




版權說明:本文檔由用戶提供并上傳,收益歸屬內容提供方,若內容存在侵權,請進行舉報或認領
文檔簡介
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*
溫馨提示
- 1. 本站所有資源如無特殊說明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
- 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請聯系上傳者。文件的所有權益歸上傳用戶所有。
- 3. 本站RAR壓縮包中若帶圖紙,網頁內容里面會有圖紙預覽,若沒有圖紙預覽就沒有圖紙。
- 4. 未經權益所有人同意不得將文件中的內容挪作商業或盈利用途。
- 5. 人人文庫網僅提供信息存儲空間,僅對用戶上傳內容的表現方式做保護處理,對用戶上傳分享的文檔內容本身不做任何修改或編輯,并不能對任何下載內容負責。
- 6. 下載文件中如有侵權或不適當內容,請與我們聯系,我們立即糾正。
- 7. 本站不保證下載資源的準確性、安全性和完整性, 同時也不承擔用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。
最新文檔
- 股份保本協議書
- 花店入伙協議書
- 租賃終止協議書
- 玉器鑒定協議書
- 統一安裝協議書
- 土地入股合作社協議書
- 破除陰婚協議書
- 職工貸款協議書
- 資產調出協議書
- 藥店代銷協議書
- 2025年商法知識競賽考試試卷及答案
- 2025年山東省臨沂市平邑縣中考一模語文試題(含答案)
- 2025年電子信息工程專業考試試題及答案
- 【威海】2025年山東省威海技師學院公開招聘工作人員29人筆試歷年典型考題及考點剖析附帶答案詳解
- 2025年第六屆全國國家版圖知識競賽題庫及答案
- 機械租賃投標服務方案
- 2025年烘焙師職業資格考試真題卷:烘焙師職業競賽與評價試題
- 2025年北京市朝陽區九年級初三一模英語試卷(含答案)
- GB 7718-2025食品安全國家標準預包裝食品標簽通則
- Unit1-Unit2重點短語(背誦版+默寫版)外研版英語新七年級下冊
- 《抗休克藥物治療》課件
評論
0/150
提交評論