[工作范文]模擬心得MATERIAL STUDIO 中SORPTION

《[工作范文]模擬心得MATERIAL STUDIO 中SORPTION》由會員分享，可在線閱讀，更多相關(guān)《[工作范文]模擬心得MATERIAL STUDIO 中SORPTION（59頁珍藏版）》請?jiān)谘b配圖網(wǎng)上搜索。

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) %+++++++

2、++++++++++++++++++++++++++++++++++++++ %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

3、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 %+++++++

4、++++++++++++++++++++++++++++++++++++++ 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]; %+++++++++++++++++++++++++++++++++++++

5、++++++++ 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=(

6、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); for k=1:length(P) while abs(Z1(k)-Z0(k))>1e-6 Z0(k)=Z1(k); Z1(k)=1+beta(k)

7、-q.*beta(k).*(Z0(k)-beta(k))./((Z0(k)+epsilon.*beta(k)).*(Z0(k)+sigma.*beta(k))); end end I=(1./(sigma-epsilon)).*log((Z1+sigma.*beta)./(Z1+epsilon.*beta)); %+++++++++++++++++++++++++++++++++++++++++++++ %gain the density of gas %+++++++++++++++++++++++++++++++++++++++++++++ rho=(P./(Z1.*R.

8、*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 b

9、inary %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 g

10、as. 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 %+++++++++++++++

11、++++++++++++++++++++++++++++++ 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

12、.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); %

13、+++++++++++++++++++++++++++++++++++++++++++++ %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

14、).^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

15、.*beta(k)).*(Z0(k)+sigma.*beta(k))); end end I=(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

16、.*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(Z

17、1-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結(jié)構(gòu)，我們需要找到具體的實(shí)驗(yàn)文獻(xiàn)和作者，然后去CCDC下載，如圖1所示。下載中需要輸入文獻(xiàn)名和作者的姓。等一會，看看輸入的那個(gè)郵箱，就會看見CIF文件了，不過得到的是TXT文件，需要改掉擴(kuò)展名，輸入MS,在MS里面手動改成文獻(xiàn)上那樣，因?yàn)橛袝r(shí)候得到的結(jié)構(gòu)會很不規(guī)則。電荷是使用DMOL3計(jì)算得到的，但是對于某些MOF,由于含有太多

18、的過渡金屬，用DMOL3優(yōu)化得到電荷效果不好，需要使用GAUSSIAN計(jì)算電荷，在使用GAUSSVIEW生成GJF文件的時(shí)候，需要在最終結(jié)果里面加入POP=CHELPG這一行，具體請看GAUSSIAN使用手冊,分析結(jié)果里面就可以看到CHELPG電荷了。得到的結(jié)構(gòu)首先需要使用分子動力學(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)同壓力與逸度差別不大，在高壓下，就要選擇

19、逸度了，如果認(rèn)為低壓下取樣數(shù)很少，就要BUILD->SYMMETRY->SUPERCELL,建立合適的超晶胞，進(jìn)行低壓下的模擬。MOF中一般來說，UFF力場與實(shí)驗(yàn)對比比較好，選擇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得到的不同，在寫文章的時(shí)候，需要轉(zhuǎn)化一下。 沸石的電荷一般用DMOL

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

21、胞，超晶胞的具體重復(fù)倍數(shù)可以視情況而定。然后使用EDIT下面的ATOM SELECTION中RAIDIAL DISTANCE,確定中心，這個(gè)就需要幾何知識了，可以確定X和Y,變動Z,也可以確定Y和Z，變動X,變動的那個(gè)數(shù)值時(shí)從0到超晶胞邊長。對于挖孔，可以只挖一個(gè)孔，也可以挖2個(gè)以上的孔，其實(shí)可以知道這些不同的中心就是不同的線段而已。經(jīng)過嘗試每隔2，可以把孔道打通的很干凈，不過此時(shí)得在EDIT下面的子菜單FIND PATTERNS 里面尋找到一個(gè)SI與三個(gè)O連接的基團(tuán)，刪除此類基團(tuán)。 然后就是加H了，加完H,還得賦予上使用量化模塊計(jì)算得到的ESP或者CHELPG電荷，使用FORCITE模塊進(jìn)

22、行MD優(yōu)化得到希望得到的MCM-41構(gòu)型。 圖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) PIPI 0 1 O 19.14690000 19.30120000 45.38230000 Zn 18.1

23、0790000 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.2

24、6360000 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.2

25、3230000 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 2

26、74 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 281 Zn 0 LANL2DZ **** C O H 0 6-31G(d) **** Zn 0 LANL2DZ Zn 1.35 具體的各參數(shù)解釋，請看GAUSSIAN使用手冊。 對于如何在分子篩的骨架添加CH2CH2NH2這樣的基團(tuán)，因?yàn)榉肿雍Y是無定形結(jié)構(gòu)，在MS里面肯定是不能自動全部添加好，只能使用編程語言添加?？梢园袽

27、CM-41的結(jié)構(gòu)保存為 CAR格式文件（為什么非要保存為CAR文件，因?yàn)镃AR文件里面有原子的電荷），然后利用隨機(jī)數(shù)發(fā)生器，在內(nèi)部隨機(jī)生成位置，只要次位置與原來骨架之間的距離小于C-SI鍵長的話，那么就認(rèn)為這個(gè)位置是可以接受的，并且把此位置命名為C原子，剩下來的C和N照樣按照這個(gè)添加，可以寫一個(gè)添加原子的子程序，調(diào)用三次就好。然后把得到的CAR文件導(dǎo)入到MS中，自動加氫就好。在MCM-41中添加胺基的源程序是這樣的： integer natom,natom0,nho,namino *********************************************

28、*************** * 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=999

29、2,nho=2319,namino=435) character a(natom)*5,fx(natom)*4,fft(natom)*5,atomname(natom)*2 integer occupation(natom),kjishu,ron,nOtemp,kstop,natom_add,Tatom integer Ohydroxy_SN(nho),Temp_SN double precision xc(natom),yc(natom),zc(natom) double precision rox,roy,roz,Ohydroxy(nho,4),Otemp(nho,

30、4) double precision OTa(4),OTb(4),temp(3),list3(nho*3,12) integer templist3(nho*3,3),Nra double 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*1 real distance,search_step,dis1,dis2,dis3,lbond,char

31、ge(natom) ****************************************************************** * define global variables ****************************************************************** common charge common xc,yc,zc ***********************************************************

32、******* * Read the input file ****************************************************************** open(10,file=MCM41-final.car,status=old) do 20 i=1,natom0 read(10,*)a(i),xc(i),yc(i),zc(i),fx(i),occupation(i),fft(i), & atomname(i),charge(i) 20 continue close(

33、10) ******************************************************************* * write the initial file to check whether the initial structure is read correctly. ******************************************************************* open(30,file=output.car,access=append) do 40 i=1,natom0 write

34、(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=0 Tatom=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=0 do 45 itt=1,namino ****

35、***************************************************** * 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) yte

36、mp=Ohydroxy(Nra,3) ztemp=Ohydroxy(Nra,4) lbond=1.90 call addatom(xtemp,ytemp,ztemp,lbond,Tatom,Temp_SN, & xfinal,yfinal,zfinal) natom_add=natom_add+1 Tatom=natom0+natom_add * NSiT=NSiT+1 * NSi_C=INT((NSiT+60)/99) * NSi_S=NSiT+60-99*NSi_C

37、* * 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)=Si xc(Tatom)=xfinal yc(Tatom)=yfinal zc(Tatom)=zfinal fx(Tat

38、om)=XXXX occupation(Tatom)=1 fft(Tatom)=Si3 atomname(Tatom)=Si charge(Tatom)=5.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),ato

39、mname(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+

40、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.templ

41、ist(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) * *

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

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

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

45、))=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 co

46、ntinue *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), * &

47、 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) ***************************************************

48、* 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 * kjis

49、hu=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,*)(Ohyd

50、roxy(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) * r

51、oz=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)

52、=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+ *

53、 & (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)=Ot

54、emp(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 ****************

55、**************************** * 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

56、 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=sqr

57、t((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. * &

58、 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(

59、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+na

60、tom_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) * wri

61、te(140,888)a(Tatom),xc(Tatom), * & yc(Tatom),zc(Tatom), * & fx(Tatom),occupation(Tatom), * & fft(Tatom),atomname(Tatom), * & charge(Tatom) * close(140) ************************************************************

62、****** * 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, &

63、 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_C

64、C=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

65、 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), &

66、 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, &