1、 模擬心得 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 fu

2、gacity 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

3、 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.9

4、9 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)

5、;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

6、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

7、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 yo

8、u 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.2

9、29 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);%+

10、%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

11、(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)

12、.*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)./b

13、m).*(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,由于含

14、有太多的過渡金屬，用dmol3優化得到電荷效果不好，需要使用gaussian計算電荷，在使用gaussview生成gjf文件的時候，需要在最終結果里面加入pop=chelpg這一行，具體請看gaussian使用手冊,分析結果里面就可以看到chelpg電荷了。得到的結構首先需要使用分子動力學進行優化，使用forcite模塊，選擇geometrical optimization這個任務，電荷加和方式最好用ewald方法，vanderwaals選擇atom based.得到的結構就可以進行sorption模擬了，選擇fixed pressure任務，在低壓下，可以認同壓力與逸度差別不大，在高壓下，就

15、要選擇逸度了，如果認為低壓下取樣數很少，就要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計算得到就可

16、以了，esp電荷比其他兩種電荷更加合適，因為計算方法比較適應真實體系。對于怎么換配體，需要點擊晶體的框架，然后擴大晶體的邊長，這樣，就在刪除配體后，有空間畫新配體了，然后用forcite優化，優化的過程中勾選上，優化晶胞的選項。金屬摻雜，是先在mof或者分子篩中切取簇模型，然后賦予這個簇模型不同的電荷，這樣再把這個得到的電荷賦予到整體的骨架中，由于此時整體的骨架電荷不為0，就需要一定數目的金屬原子去平衡它，這些金屬原子可以作為吸附質預先吸附進這個骨架。對于怎么構造mcm-41的骨架，可以使用ms自帶的strucuture中的glass下面的無定形sio2,也是通過構造超晶胞，超晶胞的具體重復倍

17、數可以視情況而定。然后使用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優化得到希望得到的m

18、cm-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

19、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.1417

20、0000 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 2

21、81zn 0lanl2dz*c o h 06-31g(d)*zn 0lanl2dzzn 1.35具體的各參數解釋，請看gaussian使用手冊。 對于如何在分子篩的骨架添加ch2ch2nh2這樣的基團，因為分子篩是無定形結構，在ms里面肯定是不能自動全部添加好，只能使用編程語言添加?？梢园裮cm-41的結構保存為 car格式文件（為什么非要保存為car文件，因為car文件里面有原子的電荷），然后利用隨機數發生器，在內部隨機生成位置，只要次位置與原來骨架之間的距離小于c-si鍵長的話，那么就認為這個位置是可以接受的，并且把此位置命名為c原子，剩下來的c和n照樣按照這個添加，可以寫一個添加原子的子

22、程序，調用三次就好。然后把得到的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

23、,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

24、 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*com

25、mon 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 file to check whether the initial structure is read correctly.* open(3

26、0,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

27、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=nsi

28、t+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(ta

29、tom)=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,kjis

30、hu+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

31、,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+

32、* & (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

33、=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(

34、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),* & t

35、emplist3(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*

36、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)* end

37、if*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)* no

38、temp=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 continu

39、e* 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)=otem

40、p(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

41、,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*