搜索
模拟心得MATERIALSTUDIO中SORPTION第一个课题是模拟金属有机框架和共价有机框架吸附CO2以及分离CO2/CH4,使用的软件是Materialstudio，使用的是Sorption模块，输入的是逸度。单组份求逸度的MATLAB程序，只需要在主程序窗口输入function[rho,f]=PengRobinson(P1,T,N)（P1,T,N是具体的数值）就可以得到不同的压力和温度下的逸度。function[rho,f]=PengRobinson(P1,T,N)%+++++++++++++++++++++++++++++++++++++++++++++%PengRobinsonisusedtocalculatethedensityandfugacityofsingle%componentgasatgivenpressurewithPeng-Robinsonequation.%PengRobinsonv1.00betaincludetheparameterofn-alkanes(1-5),CO2(6)%andCO(7).%WhereP1meansinputpressure(kPa),Tistemperature(K),Nmeanstheserialnumberofgas.rho%isdensity,fisfugacity.%e.g.Ifyouwannacalculatedensityandfugacityofmethaneat300kPa,298k,you%needinput[rho,f]=PengRobinson(300,298,1).%+++++++++++++++++++++++++++++++++++++++++++++%setphysicalparameters%+++++++++++++++++++++++++++++++++++++++++++++P=P1./100;t_M=[16.04330.07044.09758.12372.15044.0128.01];t_omiga=[0.0120.1000.1520.20.2520.2240.048];t_Tc=[190.6305.3369.8425.1469.7304.2132.9];t_Pc=[45.9948.7242.4837.9633.7073.8334.99];%+++++++++++++++++++++++++++++++++++++++++++++Tc=t_Tc(N);Pc=t_Pc(N);omiga=t_omiga(N);M=t_M(N);%+++++++++++++++++++++++++++++++++++++++++++++R=83.14;epsilon=1-2.^(0.5);sigma=1+2.^(0.5);%+++++++++++++++++++++++++++++++++++++++++++++%calculatetheZofPRequation%+++++++++++++++++++++++++++++++++++++++++++++alpha=(1+(0.37464+1.54226*omiga-0.26992*omiga.^2)*(1-(T/Tc)^(0.5))).^2;a=((0.45724*R^2*Tc^2)/Pc)*alpha;b=0.07779.*R.*Tc./Pc;beta=b.*P./(R.*T);q=a./(b.*R.*T);Z0=zeros(length(P),1);Z1=ones(length(P),1);fork=1:length(P)whileabs(Z1(k)-Z0(k))1e-6Z0(k)=Z1(k);Z1(k)=1+beta(k)-q.*beta(k).*(Z0(k)-beta(k))./((Z0(k)+epsilon.*beta(k)).*(Z0(k)+sigma.*beta(k)));endendI=(1./(sigma-epsilon)).*log((Z1+sigma.*beta)./(Z1+epsilon.*beta));%+++++++++++++++++++++++++++++++++++++++++++++%gainthedensityofgas%+++++++++++++++++++++++++++++++++++++++++++++rho=(P./(Z1.*R.*T)).*M.*1e6;rho=vpa(rho,6);phi=exp(Z1-1-log(Z1-beta)-q.*I);f=phi.*P1;f=vpa(f,5);双组份的求逸度的程序：求逸度的过程和单组份的一样。双组份的逸度求解程序：function[rho1,rho2,f1,f2]=PengRobinson_Binary(P1,T,N,y)%+++++++++++++++++++++++++++++++++++++++++++++%PengRobinsonisusedtocalculatethedensityandfugacityofbinary%componentgasatgivenpressurewithPeng-Robinsonequation.%PengRobinsonv1.00betaincludetheparameterofn-alkanes(1-5),%isoButane(6)isoPentane(7),neoPentane(8)hydrogen(9)carbondioxide(10)%WhereP1meansinputpressure(kPa),Tistemperature(K),Nmeanstheserialnumberofgas.rho%isdensity(g/m3),fisfugacity.%e.g.Ifyouwannacalculatedensityandfugacityofmixtureofmethaneandethaneat300kPa,298k,you%needinput[rho,f]=PengRobinson(300,298,[1,2],[0.5,0.5]).%+++++++++++++++++++++++++++++++++++++++++++++%setphysicalparameters%+++++++++++++++++++++++++++++++++++++++++++++P=P1./100;t_M=[16.04330.07044.09758.12372.15058.12372.15072.1502.01644.01];t_omiga=[0.0120.1000.1520.20.2520.1810.2290.197-0.2160.224];t_Tc=[190.6305.3369.8425.1469.7408.1460.39433.7533.19304.2];t_Pc=[45.9948.7242.4837.9633.7036.4833.8131.9913.1373.83];%+++++++++++++++++++++++++++++++++++++++++++++Tc=[t_Tc(N(1));t_Tc(N(2))];Pc=[t_Pc(N(1));t_Pc(N(2))];omiga=[t_omiga(N(1));t_omiga(N(2))];M=[t_M(N(1));t_M(N(2))];%+++++++++++++++++++++++++++++++++++++++++++++R=83.14;epsilon=1-2.^(0.5);sigma=1+2.^(0.5);%+++++++++++++++++++++++++++++++++++++++++++++%calculatetheZofPRequation%+++++++++++++++++++++++++++++++++++++++++++++alpha=(1+(0.37464+1.54226.*omiga-0.26992.*omiga.^2).*(1-(T./Tc).^(0.5))).^2;a=((0.45724.*R.^2.*Tc.^2)./Pc).*alpha;b=0.07779.*R.*Tc./Pc;a12=(a(1).*a(2)).^0.5;am=(y(1).^2).*a(1)+2.*y(1).*y(2).*a12+(y(2).^2)*a(2);bm=y(1).*b(1)+y(2).*b(2);beta=bm.*P./(R.*T);q=am./(bm.*R.*T);Z0=zeros(length(P),1);Z1=ones(length(P),1);fork=1:length(P)whileabs(Z1(k)-Z0(k))1e-6Z0(k)=Z1(k);Z1(k)=1+beta(k)-q.*beta(k).*(Z0(k)-beta(k))./((Z0(k)+epsilon.*beta(k)).*(Z0(k)+sigma.*beta(k)));endendI=(1./(sigma-epsilon)).*log((Z1+sigma.*beta)./(Z1+epsilon.*beta));%+++++++++++++++++++++++++++++++++++++++++++++%gainthedensityofgas%+++++++++++++++++++++++++++++++++++++++++++++%rho=(P./(Z1.*R.*T)).*M.*1e6;%rho=vpa(rho,6);rho=(P./(Z1.*R.*T)).*1e6;rho1=vpa(y(1).*rho.*M(1),6);rho2=vpa(y(2).*rho.*M(2),6);phi1=zeros(length(P),1);phi2=zeros(length(P),1);f1=zeros(length(P),1);f2=zeros(length(P),1);phi1=exp((b(1)./bm).*(Z1-1)-log(Z1-beta)-q.*I.*(2.*(y(1).*a(1)+y(2).*a12)./am-b(1)./bm));phi2=exp((b(2)./bm).*(Z1-1)-log(Z1-beta)-q.*I.*(2.*(y(2).*a(2)+y(1).*a12)./am-b(2)./bm));f1=phi1.*P1.*y(1);f2=phi2.*P1.*y(2);f1=vpa(f1,5);f2=vpa(f2,5);对于MOF结构，我们需要找到具体的实验文献和作者，然后去CCDC下载，如图1所示。下载中需要输入文献名和作者的姓。等一会，看看输入的那个邮箱，就会看见CIF文件了，不过得到的是TXT文件，需要改掉扩展名，输入MS,在MS里面手动改成文献上那样，因为有时候得到的结构会很不规则。电荷是使用DMOL3计算得到的，但是对于某些MOF,由于含有太多的过渡金属，用DMOL3优化得到电荷效果不好，需要使用GAUSSIAN计算电荷，在使用GAUSSVIEW生成GJF文件的时候，需要在最终结果里面加入POP=CHELPG这一行，具体请看GAUSSIAN使用手册,分析结果里面就可以看到CHELPG电荷了。得到的结构首先需要使用分子动力学进行优化，使用FORCITE模块，选择GEOMETRICALOPTIMIZATION这个任务，电荷加和方式最好用EWALD方法，VANDERWAALS选择ATOMBASED.得到的结构就可以进行SORPTION模拟了，选择FIXEDPRESSURE任务，在低压下，可以认同压力与逸度差别不大，在高压下，就要选择逸度了，如果认为低压下取样数很少，就要BUILD-SYMMETRY-SUPERCELL,建立合适的超晶胞，进行低压下的模拟。MOF中一般来说，UFF力场与实验对比比较好，选择UFF的比较多。计算MOF和沸石的CONNOLLYSURFACE需要用到MS中的atomvolume&surface