模拟心得MATERIALSTUDIO中SORPTION

模拟心得 MATERIAL STUDIO 中SORPTION第一个课题是模拟金属有机框架和共价有机框架吸附CO2以及分离CO2/CH4,使用的软件是Material studio，使用的是Sorption模块，输入的是逸度。单组份求逸度的MATLAB程序，只需要在主程序窗 口输入function rho,f =PengRobinson(P1,T,N)（P1,T,N是具体的数值）就可以得到不同的压力和温度下的逸度。function rho,f =PengRobinson(P1,T,N)%+%PengRobinson is used to calculate the density and fugacity of single%component gas at given pressure with Peng-Robinson equation.%PengRobinson v1.00 beta include the parameter of n-alkanes(1-5), CO2(6)%and CO(7).%Where P1 means input pressure(kPa), T is temperature(K), N means the serial number of gas. rho%is density, f is fugacity.%e.g. If you wanna calculate density and fugacity of methane at 300kPa, 298k,you%need input rho,f =PengRobinson(300,298,1). %+%set physical parameters%+P=P1./100;t_M=16.043 30.070 44.097 58.123 72.150 44.01 28.01;t_omiga=0.012 0.100 0.152 0.2 0.252 0.224 0.048;t_Tc=190.6 305.3 369.8 425.1 469.7 304.2 132.9 ;t_Pc=45.99 48.72 42.48 37.96 33.70 73.83 34.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);%+%calculate the Z of PR equation%+alpha=(1+(0.37464+1.54226*omiga-0.26992*omiga.2)*(1-(T/Tc)(0.5).2;a=(0.45724*R2*Tc2)/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);for k=1:length(P) while abs(Z1(k)-Z0(k)1e-6 Z0(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);%+%gain the density of gas%+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)%+%PengRobinson is used to calculate the density and fugacity of binary%component gas at given pressure with Peng-Robinson equation.%PengRobinson v1.00 beta include the parameter of n-alkanes(1-5),%isoButane(6) isoPentane(7), neoPentane(8) hydrogen(9) carbon dioxide(10)%Where P1 means input pressure(kPa), T is temperature(K), N means the serial number of gas. rho%is density(g/m3), f is fugacity.%e.g. If you wanna calculate density and fugacity of mixture of methane and ethane at 300kPa,298k, you%need input rho,f =PengRobinson(300,298,1,2,0.5,0.5). %+%set physical parameters%+P=P1./100;t_M=16.043 30.070 44.097 58.123 72.150 58.123 72.150 72.150 2.016 44.01;t_omiga=0.012 0.100 0.152 0.2 0.252 0.181 0.229 0.197 -0.216 0.224;t_Tc=190.6 305.3 369.8 425.1 469.7 408.1 460.39 433.75 33.19 304.2;t_Pc=45.99 48.72 42.48 37.96 33.70 36.48 33.81 31.99 13.13 73.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);%+%calculate the Z of PR equation%+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);for k=1:length(P) while abs(Z1(k)-Z0(k)1e-6 Z0(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);%+%gain the density of gas%+%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模块，选择GEOMETRICAL OPTIMIZATION这个任务，电荷加和方式最好用EWALD方法，VANDERWAALS选择ATOM BASED.得到的结构就可以进行SORPTION模拟了，选择FIXED PRESSURE任务，在低压下，可以认同压力与逸度差别不大，在高压下，就要选择逸度了，如果认为低压下取样数很少，就要BUILD-SYMMETRY-SUPERCELL,建立合适的超晶胞，进行低压下的模拟。MOF中一般来说，UFF力场与实验对比比较好，选择UFF的比较多。 计算MOF和沸石的CONNOLLY SURFACE需要用到MS中的atom volume & surfaces这个任务，CO2的CONNOLY RADIUS是0.165NM,可以再变化VDW SCALE FACTOR的时候，得到一些不同的自由体积。文献中 specific area 和 free volume的单位与MS得到的不同，在写文章的时候，需要转化一下。沸石的电荷一般用DMOL3计算得到就可以了，ESP电荷比其他两种电荷更加合适，因为计算方法比较适应真实体系。对于怎么换配体，需要点击晶体的框架，然后扩大晶体的边长，这样，就在删除配体后，有空间画新配体了，然后用FORCITE优化，优化的过程中勾选上，优化晶胞的选项。金属掺杂，是先在MOF或者分子筛中切取簇模型，然后赋予这个簇模型不同的电荷，这样再把这个得到的电荷赋予到整体的骨架中，由于此时整体的骨架电荷不为0，就需要一定数目的金属原子去平衡它，这些金属原子可以作为吸附质预先吸附进这个骨架。对于怎么构造MCM-41的骨架，可以使用MS自带的STRUCUTURE中的GLASS下面的无定形SIO2,也是通过构造超晶胞，超晶胞的具体重复倍数可以视情况而定。然后使用EDIT下面的ATOM SELECTION中RAIDIAL DISTANCE,确定中心，这个就需要几何知识了，可以确定X和Y,变动Z,也可以确定Y和Z，变动X,变动的那个数值时从0到超晶胞边长。对于挖孔，可以只挖一个孔，也可以挖2个以上的孔，其实可以知道这些不同的中心就是不同的线段而已。经过尝试每隔2，可以把孔道打通的很干净，不过此时得在EDIT下面的子菜单FIND PATTERNS 里面寻找到一个SI与三个O连接的基团，删除此类基团。然后就是加H了，加完H,还得赋予上使用量化模块计算得到的ESP或者CHELPG电荷，使用FORCITE模块进行MD优化得到希望得到的MCM-41构型。图1 CCDC要求CIF文件的界面附录：用GAUSSVIEW写的MOF簇的GAUSSIAN输入文件：% chk=1.chk%mem=100MW# b3lyp geom=connectivity gen pseudo=read sp scf=tight pop=(mk,chelpg)maxdisk=25gb iop(6/33=2)PIPI0 1 O 19.14690000 19.30120000 45.38230000 Zn 18.10790000 18.22300000 44.34000000 Zn 20.18800000 18.26520000 46.47520000 Zn 18.08020000 20.39100000 46.39280000 Zn 20.21140000 20.32720000 44.31250000 Zn 6.68970000 19.24060000 44.93580000 Zn 19.26360000 6.83520000 45.02830000 Zn 19.18610000 19.14570000 32.92150000 Zn 31.60320000 19.10710000 45.02760000 Zn 19.18580000 31.75920000 44.84940000 C 9.94100000 19.27140000 45.14170000 C 19.23230000 10.08690000 45.26670000 C 19.16500000 28.50940000 45.03930000 C 19.20710000 19.14190000 36.17600000 C 28.35370000 19.17540000 45.25230000 N 7.87110000 17.86280000 44.83600000. 271 272 274 1.5 273 1.0 273 274 276 1.5 275 276 278 1.0 277 278 1.0 280 1.0 278 279 1.0 279 281 1.0 280 281Zn 0LANL2DZ*C O H 06-31G(d)*Zn 0LANL2DZZn 1.35具体的各参数解释，请看GAUSSIAN使用手册。 对于如何在分子筛的骨架添加CH2CH2NH2这样的基团，因为分子筛是无定形结构，在MS里面肯定是不能自动全部添加好，只能使用编程语言添加。可以把MCM-41的结构保存为 CAR格式文件（为什么非要保存为CAR文件，因为CAR文件里面有原子的电荷），然后利用随机数发生器，在内部随机生成位置，只要次位置与原来骨架之间的距离小于C-SI键长的话，那么就认为这个位置是可以接受的，并且把此位置命名为C原子，剩下来的C和N照样按照这个添加，可以写一个添加原子的子程序，调用三次就好。然后把得到的CAR文件导入到MS中，自动加氢就好。在MCM-41中添加胺基的源程序是这样的： integer natom,natom0,nho,namino* Number of atoms in the original structure is 9992.* But the parameter natom should include the number of atoms added subsequently.* So here the value of natom is set to 15000* parameter (natom=15000,natom0=9992,nho=2319,namino=435)character a(natom)*5,fx(natom)*4,fft(natom)*5,atomname(natom)*2integer occupation(natom),kjishu,ron,nOtemp,kstop,natom_add,Tatominteger Ohydroxy_SN(nho),Temp_SNdouble precision xc(natom),yc(natom),zc(natom)double precision rox,roy,roz,Ohydroxy(nho,4),Otemp(nho,4)double precision OTa(4),OTb(4),temp(3),list3(nho*3,12)integer templist3(nho*3,3),Nradouble precision xtemp,ytemp,ztemp,xfinal,yfinal,zfinal*integer NSiT,NSi_S,NSi_C,NCT,NC_S,NC_C,NNT,NN_S,NN_C*character NSi_CC*1,NC_CC*1,NN_CC*1real distance,search_step,dis1,dis2,dis3,lbond,charge(natom)* define global variables*common charge common xc,yc,zc * Read the input file*open(10,file=MCM41-final.car,status=old)do 20 i=1,natom0read(10,*)a(i),xc(i),yc(i),zc(i),fx(i),occupation(i),fft(i), & atomname(i),charge(i)20 continue close(10)* write the initial check whether the initial structure is read correctly.* open(30,file=output.car,access=append)do 40 i=1,natom0write(30,888)a(i),xc(i),yc(i),zc(i),fx(i),occupation(i),fft(i), & atomname(i),charge(i)40 continue close(30)natom_add=0Tatom=natom0+natom_add*NSiT=0*NSi_S=0*NSi_C=0*NCT=0*NC_S=0*NC_C=0*NNT=0*NN_S=0*NN_C=0do 45 itt=1,namino* add Si to the chosen Oxygen atom* call HO_list(Ohydroxy,Ohydroxy_SN,kjishu) Nra=int(RAN2(IDUM)*kjishu)Temp_SN=Ohydroxy_SN(Nra)xtemp=Ohydroxy(Nra,2)ytemp=Ohydroxy(Nra,3)ztemp=Ohydroxy(Nra,4) lbond=1.90call addatom(xtemp,ytemp,ztemp,lbond,Tatom,Temp_SN, & xfinal,yfinal,zfinal) natom_add=natom_add+1Tatom=natom0+natom_add*NSiT=NSiT+1*NSi_C=INT(NSiT+60)/99)*NSi_S=NSiT+60-99*NSi_C*if(NSi_C.eq.0)NSi_CC=T*if(NSi_C.eq.1)NSi_CC=U*if(NSi_C.eq.2)NSi_CC=V*if(NSi_C.eq.3)NSi_CC=W*if(NSi_C.eq.4)NSi_CC=X*if(NSi_C.eq.5)NSi_CC=Y*if(NSi_C.eq.6)NSi_CC=Z*a(Tatom)=Si/NSi_S/NSi_CC a(Tatom)=Sixc(Tatom)=xfinalyc(Tatom)=yfinalzc(Tatom)=zfinalfx(Tatom)=XXXXoccupation(Tatom)=1fft(Tatom)=Si3 atomname(Tatom)=Sicharge(Tatom)=5.000open(140,file=amino.car,access=append)write(140,888)a(Tatom),xc(Tatom), & yc(Tatom),zc(Tatom), & fx(Tatom),occupation(Tatom), & fft(Tatom),atomname(Tatom), & charge(Tatom)close(140)* do 2060 i=1,nho-2* do 2065 ix=-2,kjishu+2* if(Ohydroxy(i,1).gt.templist(ix)-0.1.and.* & Ohydroxy(i,1).lt.templist(ix)+0.1)then*goto 2060*endif*2065 continue * * do 2070 j=i+1,nho-1* do 2075 ix=-2,kjishu+2* if(Ohydroxy(j,1).gt.templist(ix)-0.1.and.* & Ohydroxy(j,1).lt.templist(ix)+0.1)then* goto 2070* endif*2075 continue * do 2080 k=j+1,nho* do 2085 ix=-2,kjishu+2* if(Ohydroxy(k,1).gt.templist(ix)-0.1.and.* & Ohydroxy(k,1).lt.templist(ix)+0.1)then* goto 2080* endif*2085 continue * dis1=sqrt(Ohydroxy(i,2)-Ohydroxy(j,2)*2+* & (Ohydroxy(i,3)-Ohydroxy(j,3)*2+* & (Ohydroxy(i,4)-Ohydroxy(j,4)*2)* dis2=sqrt(Ohydroxy(i,2)-Ohydroxy(k,2)*2+* & (Ohydroxy(i,3)-Ohydroxy(k,3)*2+* & (Ohydroxy(i,4)-Ohydroxy(k,4)*2)* dis3=sqrt(Ohydroxy(j,2)-Ohydroxy(k,2)*2+* & (Ohydroxy(j,3)-Ohydroxy(k,3)*2+* & (Ohydroxy(j,4)-Ohydroxy(k,4)*2)* if (dis1.gt.1.8.and.dis1.lt.2.6.and.* & dis2.gt.1.8.and.dis1.lt.2.6.and.* & dis3.gt.1.8.and.dis1.lt.2.6)then* kjishu=kjishu+1* do 2090 imm=1,4* list3(kjishu,imm)=Ohydroxy(i,imm)*2090 continue* templist3(kjishu,1)=int(Ohydroxy(i,1)* charge(Ohydroxy_SN(i)=2.0* do 2100 imm=1,4* list3(kjishu,imm+4)=Ohydroxy(j,imm)*2100 continue* templist3(kjishu,2)=int(Ohydroxy(j,1)* charge(Ohydroxy_SN(j)=2.0* do 2110 imm=1,4* list3(kjishu,imm+8)=Ohydroxy(k,imm)*2110 continue* templist3(kjishu,3)=int(Ohydroxy(k,1)* charge(Ohydroxy_SN(k)=2.0* goto 2061* endif*2080 continue*2070 continue*2060 continue * open(2120,file=list3.car,access=append)*do 2130 i=1,kjishu* write(2120,777)int(list3(i,1),list3(i,2),list3(i,3),list3(i,4),* & templist3(i,1)* write(2120,777)int(list3(i,5),list3(i,6),list3(i,7),list3(i,8),* & templist3(i,2)* write(2120,777)int(list3(i,9),list3(i,10),list3(i,11),list3(i,12)* & ,templist3(i,3)*write(2120,*)*2130 continue*777 Format(I4,3X,F12.9,2X,F13.9,3X,F12.9,3X,I4)* close(2120)* renew hydroxy oxygen list* do 2135 i=1,nho*Ohydroxy_SN(i)=0* do 2136 j=1,4* Ohydroxy(i,j)=0*2136 continue*2135 continue* kjishu=0* do 2140 i=1,natom0* if (charge(i).eq.1.0) then* kjishu=kjishu+1* Ohydroxy_SN(kjishu)=i* Ohydroxy(kjishu,1)=kjishu* Ohydroxy(kjishu,2)=xc(i)* Ohydroxy(kjishu,3)=yc(i)* Ohydroxy(kjishu,4)=zc(i)* endif*2140 continue* open(2150,file=hydrogen oxygen-1.car,access=append)* do 2160 i=1,nho*write(2150,*)(Ohydroxy(i,j),j=1,4),Ohydroxy_SN(i)*2160 continue* close(2150)* stop* choose a initial oxygen atom randomly* ron=int(ran2(idum)*nho)* rox=Ohydroxy(nho,2)* roy=Ohydroxy(nho,3)* roz=Ohydroxy(nho,4)* nOtemp=0* do 60 i=1,nho* distance=sqrt(rox-Ohydroxy(i,2)*2+(roy-Ohydroxy(i,3)*2+* & (roz-Ohydroxy(i,4)*2)* if(distance.gt.1.0.and.distance.lt.2.6)then* nOtemp=nOtemp+1* Otemp(nOtemp,1)=nOtemp* Otemp(nOtemp,2)=Ohydroxy(i,2)* Otemp(nOtemp,3)=Ohydroxy(i,3)* Otemp(nOtemp,4)=Ohydroxy(i,4)* endif*60 continue* kstop=0* do 70 i=1,nOtemp-1* do 80 j=i+1,nOtemp* distance=sqrt(Ohydroxy(i,2)-Ohydroxy(j,2)*2+* & (Ohydroxy(i,3)-Ohydroxy(j,3)*2+* & (Ohydroxy(i,4)-Ohydroxy(j,4)*2)* if(distance.gt.1.0.and.distance.lt.2.6)then* Mark of interrupt service routine 2* kstop=1 * OTa(1)=i* OTa(2)=Otemp(i,2)* OTa(3)=Otemp(i,3)* OTa(4)=Otemp(i,4)* OTb(1)=j* OTb(2)=Otemp(j,2)* OTb(3)=Otemp(j,3)* OTb(4)=Otemp(j,4)* goto 90* endif*80 continue*70 continue* warning 1* if(kstop.eq.0)then*write(*,*)Warning,kstop,: no tri-grafting any more*stop*endif * warning 1* kstop=0*90 search_step=0.01* do 100 i=-300,300* do 110 j=-300,300* do 120 k=-300,300* temp(1)=rox+search_step*i* temp(2)=roy+search_step*j* temp(3)=roz+search_step*k* dis1=sqrt(temp(1)-rox)*2+* & (temp(2)-roy)*2+* & (temp(3)-roz)*2)* dis2=sqrt(temp(1)-OTa(2)*2+* & (temp(2)-OTa(3)*2+* & (temp(3)-OTa(4)*2) * dis3=sqrt(temp(1)-OTb(2)*2+* & (temp(2)-OTb(3)*2+* & (temp(3)-OTb(4)*2) * if (dis1.gt.1.8.and.dis1.lt.2.0.and.* & dis2.gt.1.8.and.dis1.lt.2.0.and.* & dis3.gt.1.8.and.dis1.lt.2.0)then* do 130 im=1,natom+natom_add* distance=sqrt(temp(1)-xc(im)*2+* & (temp(2)-yc(im)*2+* & (temp(3)-zc(im)*2)* if(distance.le.3.0)then*goto 120*endif*130 continue* Mark of interrupt service routine 2* kstop=2* natom_add=natom_add+1* Tatom=natom0+natom_add* a(Tatom)=Si* xc(Tatom)=temp(1)* yc(Tatom)=temp(2)* zc(Tatom)=temp(3)* fx(Tatom)=XXXX* occupation(Tatom)=1* fft(Tatom)=Si3* atomname(Tatom)=Si* charge(Tatom)=3.000* open(140,file=amino.car,access=append)* write(140,888)a(Tatom),xc(Tatom),* & yc(Tatom),zc(Tatom),* & fx(Tatom),occupation(Tatom),* & fft(Tatom),atomname(Tatom),* & charge(Tatom)* close(140) * add C1 atom on the chosen Si atom* lbond=1.93 xtemp=xc(Tatom) ytemp=yc(Tatom) ztemp=zc(Tatom) Temp_SN=Tatom call addatom(xtemp,ytemp,ztemp,lbond,Tatom,Temp_SN, & xfinal,yfinal,zfinal) natom_add=natom_add+1 Tatom=natom0+natom_add* NCT=NCT+1* NC_C=INT(NCT/99)* NC_S=NCT-99*NC_C* if(NC_C.eq.0)NC_CC=* if(NC_C.eq.1)NC_CC=A* if(NC_C.eq.2)NC_CC=B* if(NC_C.eq.3)NC_CC=C* if(NC_C.eq.4)NC_CC=D* if(NC_C.eq.5)NC_CC=E* if(NC_C.eq.6)NC_CC=F* a(Tatom)=C/NC_S/NC_CC a(Tatom)=C xc(Tatom)=xfinal yc(Tatom)=yfinal zc(Tatom)=zfinal fx(Tatom)=XXXX occupation(Tatom)=1 fft(Tatom)=C_3 atomname(Tatom)=C charge(Tatom)=6.000 open(150,file=amino.car,access=append) write(150,888)a(Tatom),xc(Tatom), & yc(Tatom),zc(Tatom), & fx(Tatom),occupation(Tatom), & fft(Tatom),atomname(Tatom), & charge(Tatom) close(150)* add C2 atom on the chosen C1 atom* lbond=1.50 xtemp=xc(Tatom) ytemp=yc(Tatom) ztemp=zc(Tatom) Temp_SN=Tatom call addatom(xtemp,ytemp,ztemp,lbond,Tatom,Temp_SN, &