1、 模擬心得 material studio 中sorption第一個(gè)課題是模擬金屬有機(jī)框架和共價(jià)有機(jī)框架吸附co2以及分離co2/ch4,使用的軟件是material studio，使用的是sorption模塊，輸入的是逸度。單組份求逸度的matlab程序，只需要在主程序窗 口輸入function rho,f =pengrobinson(p1,t,n)（p1,t,n是具體的數(shù)值）就可以得到不同的壓力和溫度下的逸度。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);雙組份的求逸度的程序：求逸度的過(guò)程和單組份的一樣。雙組份的逸度求解程序：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);對(duì)于mof結(jié)構(gòu)，我們需要找到具體的實(shí)驗(yàn)文獻(xiàn)和作者，然后去ccdc下載，如圖1所示。下載中需要輸入文獻(xiàn)名和作者的姓。等一會(huì)，看看輸入的那個(gè)郵箱，就會(huì)看見(jiàn)cif文件了，不過(guò)得到的是txt文件，需要改掉擴(kuò)展名，輸入ms,在ms里面手動(dòng)改成文獻(xiàn)上那樣，因?yàn)橛袝r(shí)候得到的結(jié)構(gòu)會(huì)很不規(guī)則。電荷是使用dmol3計(jì)算得到的，但是對(duì)于某些mof,由于含

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

15、要選擇逸度了，如果認(rèn)為低壓下取樣數(shù)很少，就要build-symmetry-supercell,建立合適的超晶胞，進(jìn)行低壓下的模擬。mof中一般來(lái)說(shuō)，uff力場(chǎng)與實(shí)驗(yàn)對(duì)比比較好，選擇uff的比較多。 計(jì)算mof和沸石的connolly surface需要用到ms中的atom volume & surfaces這個(gè)任務(wù)，co2的connoly radius是0.165nm,可以再變化vdw scale factor的時(shí)候，得到一些不同的自由體積。文獻(xiàn)中 specific area 和 free volume的單位與ms得到的不同，在寫(xiě)文章的時(shí)候，需要轉(zhuǎn)化一下。沸石的電荷一般用dmol3計(jì)算得到就可

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

17、數(shù)可以視情況而定。然后使用edit下面的atom selection中raidial distance,確定中心，這個(gè)就需要幾何知識(shí)了，可以確定x和y,變動(dòng)z,也可以確定y和z，變動(dòng)x,變動(dòng)的那個(gè)數(shù)值時(shí)從0到超晶胞邊長(zhǎng)。對(duì)于挖孔，可以只挖一個(gè)孔，也可以挖2個(gè)以上的孔，其實(shí)可以知道這些不同的中心就是不同的線段而已。經(jīng)過(guò)嘗試每隔2，可以把孔道打通的很干凈，不過(guò)此時(shí)得在edit下面的子菜單find patterns 里面尋找到一個(gè)si與三個(gè)o連接的基團(tuán)，刪除此類(lèi)基團(tuán)。然后就是加h了，加完h,還得賦予上使用量化模塊計(jì)算得到的esp或者chelpg電荷，使用forcite模塊進(jìn)行md優(yōu)化得到希望得到的m

18、cm-41構(gòu)型。圖1 ccdc要求cif文件的界面附錄：用gaussview寫(xiě)的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具體的各參數(shù)解釋，請(qǐng)看gaussian使用手冊(cè)。 對(duì)于如何在分子篩的骨架添加ch2ch2nh2這樣的基團(tuán)，因?yàn)榉肿雍Y是無(wú)定形結(jié)構(gòu)，在ms里面肯定是不能自動(dòng)全部添加好，只能使用編程語(yǔ)言添加?？梢园裮cm-41的結(jié)構(gòu)保存為 car格式文件（為什么非要保存為car文件，因?yàn)閏ar文件里面有原子的電荷），然后利用隨機(jī)數(shù)發(fā)生器，在內(nèi)部隨機(jī)生成位置，只要次位置與原來(lái)骨架之間的距離小于c-si鍵長(zhǎng)的話，那么就認(rèn)為這個(gè)位置是可以接受的，并且把此位置命名為c原子，剩下來(lái)的c和n照樣按照這個(gè)添加，可以寫(xiě)一個(gè)添加原子的子

22、程序，調(diào)用三次就好。然后把得到的car文件導(dǎo)入到ms中，自動(dòng)加氫就好。在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*