模拟心得MATERIAL STUDIO 中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);

%+++++++++++++++++++++++++++++++++++++++++++++ êlculate 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; 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))); 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);

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);

%+++++++++++++++++++++++++++++++++++++++++++++ êlculate 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))); 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.*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,由

3

于含有太多的过渡金属，用DMOL优化得到电荷效果不好，需要使用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得到的不同，在写文章的时候，需要转化一下。

沸石的电荷一般用DMOL计算得到就可以了，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构型。

3

图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.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 281

Zn 0 LANL2DZ **** C O H 0 6-31G(d) ****

Zn 0 LANL2DZ

Zn 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)

* *

****************************************************************** * define global variables

******************************************************************

****************************************************************** * 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), common charge common xc,yc,zc

real distance,search_step,dis1,dis2,dis3,lbond,charge(natom) 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,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

& atomname(i),charge(i) 20 continue close(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(30,888)a(i),xc(i),yc(i),zc(i),fx(i),occupation(i),fft(i),

& atomname(i),charge(i) 40 continue close(30) * * * * * * * * *

do 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)

call addatom(xtemp,ytemp,ztemp,lbond,Tatom,Temp_SN,

NNT=0 NN_S=0 NN_C=0 NCT=0 NC_S=0 NC_C=0 NSiT=0 NSi_S=0 NSi_C=0 natom_add=0

Tatom=natom0+natom_add

lbond=1.90

& xfinal,yfinal,zfinal) * * * * * *

if(NSi_C.eq.0)NSi_CC='T' if(NSi_C.eq.1)NSi_CC='U' NSiT=NSiT+1

NSi_C=INT((NSiT+60)/99) NSi_S=NSiT+60-99*NSi_C natom_add=natom_add+1

Tatom=natom0+natom_add

* * * * * *

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' xc(Tatom)=xfinal yc(Tatom)=yfinal zc(Tatom)=zfinal fx(Tatom)='XXXX' occupation(Tatom)=1 fft(Tatom)='Si3' charge(Tatom)=5.000

open(140,file='amino.car',access='append') write(140,888)a(Tatom),xc(Tatom),

a(Tatom)='Si'

atomname(Tatom)='Si'

& yc(Tatom),zc(Tatom), & fx(Tatom),occupation(Tatom), & fft(Tatom),atomname(Tatom), & charge(Tatom)

* do 2060 i=1,nho-2 * * * * * * * * * *

do 2075 ix=-2,kjishu+2

if(Ohydroxy(j,1).gt.templist(ix)-0.1.and. do 2065 ix=-2,kjishu+2

if(Ohydroxy(i,1).gt.templist(ix)-0.1.and.

goto 2060 endif close(140)

* & Ohydroxy(i,1).lt.templist(ix)+0.1)then

*2065 continue

do 2070 j=i+1,nho-1

* & 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.

goto 2080 endif

* & Ohydroxy(k,1).lt.templist(ix)+0.1)then

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

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 kjishu=kjishu+1

* list3(kjishu,imm+8)=Ohydroxy(k,imm) *2110 continue

* templist3(kjishu,3)=int(Ohydroxy(k,1)) * charge(Ohydroxy_SN(k))=2.0 *

*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 goto 2061 * endif

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

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

* kstop=0

* do 70 i=1,nOtemp-1 * *

do 80 j=i+1,nOtemp

distance=sqrt((Ohydroxy(i,2)-Ohydroxy(j,2))**2+

distance=sqrt((rox-Ohydroxy(i,2))**2+(roy-Ohydroxy(i,3))**2+ if(distance.gt.1.0.and.distance.lt.2.6)then Otemp(nOtemp,1)=nOtemp Otemp(nOtemp,2)=Ohydroxy(i,2)

Otemp(nOtemp,3)=Ohydroxy(i,3) Otemp(nOtemp,4)=Ohydroxy(i,4)

* & (roz-Ohydroxy(i,4))**2) * nOtemp=nOtemp+1

Ohydroxy(kjishu,2)=xc(i) Ohydroxy(kjishu,3)=yc(i) Ohydroxy(kjishu,4)=zc(i) endif

*2140 continue

endif

*60 continue

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

************************************************* * * * * * * * * * *

goto 90 * endif *80 continue *70 continue

******************************************** * warning 1

******************************************** * if(kstop.eq.0)then * * *

write(*,*)'Warning',kstop,': no tri-grafting any more' stop endif OTb(1)=j

OTb(2)=Otemp(j,2) OTb(3)=Otemp(j,3) OTb(4)=Otemp(j,4) OTa(1)=i

OTa(2)=Otemp(i,2) OTa(3)=Otemp(i,3) OTa(4)=Otemp(i,4) kstop=1

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

open(140,file='amino.car',access='append') write(140,888)a(Tatom),xc(Tatom), 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' charge(Tatom)=3.000

* atomname(Tatom)='Si'

* & 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) * * * * * * * * * * *

open(150,file='amino.car',access='append') a(Tatom)='C'//'NC_S'//'NC_CC' xc(Tatom)=xfinal zc(Tatom)=zfinal fx(Tatom)='XXXX' occupation(Tatom)=1 fft(Tatom)='C_3' charge(Tatom)=6.000 a(Tatom)='C' yc(Tatom)=yfinal

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' NCT=NCT+1 NC_C=INT(NCT/99) NC_S=NCT-99*NC_C natom_add=natom_add+1

Tatom=natom0+natom_add

atomname(Tatom)='C'

write(150,888)a(Tatom),xc(Tatom),

& yc(Tatom),zc(Tatom), & fx(Tatom),occupation(Tatom), & fft(Tatom),atomname(Tatom), & charge(Tatom)

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

close(150)

& xfinal,yfinal,zfinal) * * * * * * * * * * *

a(Tatom)='C'//'NC_S'//'NC_CC' xc(Tatom)=xfinal yc(Tatom)=yfinal zc(Tatom)=zfinal fx(Tatom)='XXXX' occupation(Tatom)=1 fft(Tatom)='C_3' charge(Tatom)=7.000 a(Tatom)='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' NCT=NCT+1 NC_C=INT(NCT/99) NC_S=NCT-99*NC_C natom_add=natom_add+1

Tatom=natom0+natom_add

atomname(Tatom)='C'

open(160,file='amino.car',access='append') write(160,888)a(Tatom),xc(Tatom),

& yc(Tatom),zc(Tatom), & fx(Tatom),occupation(Tatom), & fft(Tatom),atomname(Tatom), & charge(Tatom)

****************************************************************** * add C3 atom on the chosen C2 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,

close(160)

& xfinal,yfinal,zfinal) * * * * * * * * * * *

a(Tatom)='C'//'NC_S'//'NC_CC' xc(Tatom)=xfinal yc(Tatom)=yfinal zc(Tatom)=zfinal fx(Tatom)='XXXX' occupation(Tatom)=1 fft(Tatom)='C_3' a(Tatom)='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' NCT=NCT+1 NC_C=INT(NCT/99) NC_S=NCT-99*NC_C natom_add=natom_add+1

Tatom=natom0+natom_add

atomname(Tatom)='C'

open(170,file='amino.car',access='append') write(170,888)a(Tatom),xc(Tatom), charge(Tatom)=8.000

& yc(Tatom),zc(Tatom), & fx(Tatom),occupation(Tatom), & fft(Tatom),atomname(Tatom), & charge(Tatom)

****************************************************************** * add N atom on the chosen C3 atom

******************************************************************

lbond=1.53

xtemp=xc(Tatom) ytemp=yc(Tatom) ztemp=zc(Tatom) Temp_SN=Tatom

call addatom(xtemp,ytemp,ztemp,lbond,Tatom,Temp_SN,

close(170)

& xfinal,yfinal,zfinal) * * * * * * * * * * *

a(Tatom)='N'//'NN_S'//'NN_CC' xc(Tatom)=xfinal yc(Tatom)=yfinal zc(Tatom)=zfinal a(Tatom)='N'

if(NN_C.eq.0)NN_CC='' if(NN_C.eq.1)NN_CC='A' if(NN_C.eq.2)NN_CC='B' if(NN_C.eq.3)NN_CC='C' if(NN_C.eq.4)NN_CC='D' if(NN_C.eq.5)NN_CC='E' if(NN_C.eq.6)NN_CC='F' NNT=NNT+1 NN_C=INT(NNT/99) NN_S=NNT-99*NN_C natom_add=natom_add+1

Tatom=natom0+natom_add

fx(Tatom)='XXXX' occupation(Tatom)=1 fft(Tatom)='N_3' charge(Tatom)=9.000

open(180,file='amino.car',access='append') write(180,888)a(Tatom),xc(Tatom),

atomname(Tatom)='N'

& yc(Tatom),zc(Tatom), & fx(Tatom),occupation(Tatom), & fft(Tatom),atomname(Tatom), & charge(Tatom)

45 continue

**************************************************** * rewrite the atom serial.

**************************************************** open(46,file='Satom.car',status='old')

do 47 i=natom0+1,Tatom read(46,1888)a(i) close(180)

47 continue close(46) 1888 format(A5) * *

*120 continue *110 continue *100 continue

endif

*************************************************** * write the output file

***************************************************

open(190,file='output-final.car',access='append') do 200 i=1,Tatom

write(190,888)a(i),xc(i),yc(i),zc(i),fx(i),occupation(i),fft(i),

& atomname(i),charge(i) 200 continue close(190)

888 format(A5,3X,F12.9,2X,F13.9,3X,F12.9,1X,A4,1X,I1,6X,A5,3X,A2,F7.3)

**************************************************** ****************************************************

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

* To make a list of hydroxy oxygen.

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

subroutine HO_list(Ohydroxy,Ohydroxy_SN,kjishu)

integer nho,natom,natom0

parameter (natom=15000,natom0=9992,nho=2319) integer kjishu,Ohydroxy_SN(nho) double precision Ohydroxy(nho,4)

end

******************************************************************* * declare the global varialbe

*******************************************************************

kjishu=0 do 50 i=1,natom0

if (charge(i).eq.1.0) then kjishu=kjishu+1 Ohydroxy_SN(kjishu)=i Ohydroxy(kjishu,2)=xc(i) Ohydroxy(kjishu,3)=yc(i) Ohydroxy(kjishu,4)=zc(i) endif

real charge(natom)

double precision xc(natom),yc(natom),zc(natom) common charge common xc,yc,zc

Ohydroxy(kjishu,1)=kjishu

50 continue

****************************************************************** * write out the hydroxy oxygen list

* *

if(dis2.lt.3.0)then goto 1022 endif endif

if(dis1.gt.minlbond.and.dis1.lt.maxlbond)then dis2=sqrt((targetx-xc(im))**2+

* do 1032 im=1,Tatom

* & (targety-yc(im))**2+ * & (targetz-zc(im))**2) * * *

1032 continue goto 1042 endif 1022 continue 1012 continue 1002 continue

1042 xfinal=targetx yfinal=targety

endsubroutine

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

endif zfinal=targetz

if(im.ne.Temp_SN.and.dis2.lt.3.0)then goto 1022 endif

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

* Random Number Generator

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV REAL ran2,AM,EPS,RNMX

PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1, *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791, *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)

INTEGER idum2,j,k,iv(NTAB),iy SAVE iv,iy,idum2

DATA idum2/123456789/, iv/NTAB*0/, iy/0/ if (idum.le.0) then idum=max(-idum,1) idum2=idum do 11 j=NTAB+8,1,-1 k=idum/IQ1

idum=IA1*(idum-k*IQ1)-k*IR1 if (idum.lt.0) idum=idum+IM1 if (j.le.NTAB) iv(j)=idum 11 continue iy=iv(1) endif k=idum/IQ1

idum=IA1*(idum-k*IQ1)-k*IR1 if (idum.lt.0) idum=idum+IM1 k=idum2/IQ2

idum2=IA2*(idum2-k*IQ2)-k*IR2 if (idum2.lt.0) idum2=idum2+IM2 j=1+iy/NDIV

FUNCTION RAN2(IDUM)

iy=iv(j)-idum2 iv(j)=idum

if(iy.lt.1)iy=iy+IMM1 ran2=min(AM*iy,RNMX) return END

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

各段程序大概信息已经在程序内部给出，调用添加原子的子程序3次，得到的CAR文件输入MS就可以自动加氢了。这个程序中使用的MCM-41模型还是单个晶胞4孔模型（由3*7*2超晶胞挖出来的），对于单孔模型的MCM-41需要自己改写，在CENTER.TXT那里修改，后面的搜索步骤也要相应的修改。其他种类的沸石例如MFI，LTA等等加胺基也是按照这样做的。对于怎么计算孔径分布，需要到处CAR文件中所有的H原子，其源程序是这样的：

integer hno

open(10,file='H.car',status='old')

do 20 i=1,hno

read(10,*)a(i),xc(i),yc(i),zc(i),fx(i),occupation(i),fft(i), character a(hno)*5,fx(hno)*4,fft(hno)*5,natom(hno)*1 double precision xc(hno),yc(hno),zc(hno) integer occupation(hno) real charge(hno) parameter(hno=2356)

& natom(i),charge(i) 20 continue close(10)

jmm=0

open(30,file='output.car',access='append') *40 amm=RAN2(IDUM) 40 imm=int(RAN2(IDUM)*hno)

if(charge(imm).lt.1.0)then jmm=jmm+1

write(30,888)a(imm),xc(imm),yc(imm),zc(imm),fx(imm), charge(imm)=2.0 endif

& occupation(imm),fft(imm),natom(imm),charge(imm)

888 format(A5,3X,F12.9,2X,F13.9,3X,F12.9,1X,A4,1X,I1,6X,A5,3X,A1,F8.3)

end

if(jmm.lt.415)then else goto 50 endif goto 40

if(jmm.eq.50.or.jmm.eq.100.or.jmm.eq.150.or.jmm.eq.200.or.

jmm.eq.250.or.jmm.eq.300.or.jmm.eq.350.or.jmm.eq.400)then

write(30,*)'********************',jmm,'**************************' write(30,*)'********************',jmm,'**************************' endif &

50 close(30)

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

* Random Number Generator

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV REAL ran2,AM,EPS,RNMX

PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1, *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791, *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)

INTEGER idum2,j,k,iv(NTAB),iy SAVE iv,iy,idum2

DATA idum2/123456789/, iv/NTAB*0/, iy/0/ if (idum.le.0) then idum=max(-idum,1) idum2=idum do 11 j=NTAB+8,1,-1 k=idum/IQ1

idum=IA1*(idum-k*IQ1)-k*IR1 if (idum.lt.0) idum=idum+IM1 if (j.le.NTAB) iv(j)=idum 11 continue

FUNCTION RAN2(IDUM)

iy=iv(1) endif k=idum/IQ1

idum=IA1*(idum-k*IQ1)-k*IR1 if (idum.lt.0) idum=idum+IM1 k=idum2/IQ2

idum2=IA2*(idum2-k*IQ2)-k*IR2 if (idum2.lt.0) idum2=idum2+IM2 j=1+iy/NDIV iy=iv(j)-idum2 iv(j)=idum

if(iy.lt.1)iy=iy+IMM1 ran2=min(AM*iy,RNMX) return END

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

导出的氢原子的坐标需要按照孔径分布的函数关系式，输入孔径分布的自编程序，自编程序是这样的：

Towhee的安装：

首先要保证机器中有gfortran编译器，没有的话，我用的是FEDORA12,可以使用yum install gcc-gfortran下载一个，安装好了后，下载towhee的源代码，就可以开始编译了。 ln -s /home/zyj/zhuomian/towhee-6.2.7 /towheebase（这个必须得cd到最上层才行） cd /home/zyj/zhuomian/towhee-6.2.7 ./configure --enable-fix-GNU cd Source

make towhee

然后把里面do_test拷贝到任何一个例子里面，然后就运行./do_test，例子就可以开始跑了。Towhee能够做GBMC,并且input文件里面，主要是分子的构型信息。

GCMC程序MUSIC的编译以及使用过程

要使用Music，首先就得安装，不过安装过程肯定不想Windows操作系统下面那么方便，点点鼠标就OK了，在这里需要使用命令行。由于music使用了Fortran 90编写，因而不能和towhee一样，采用Linux自带的gcc就能编译，也许有人会问了，我把gcc更新或者采用最新的gcc也不行吗，答曰不行，因为据我了解，gcc强大在于对C语言的编译，而对于Fortran，只能支持F77，对于F90只能爱莫能助。因此我们在编译music之前，需要安装一个F90编译器。Music网页上推荐了几款，这里我们推荐使用intel Fortran compiler（以下简称ifort），因为它有非商业版本下载免费使用，只要你不是用于商业用途，就可以，但是你也不会得到

IRMOF1.Methane.res # Restart File to write to

IRMOF1.Methane.con # Configuration File

首先，这个文件里出现的，#后面的内容为标注内容，可有可无，但是推荐使用，因为便于日后检查，其次，里面有的空行是必须存在的，不知道理由，反正空了就不会错，不空有时会影响程序运行，反正按照这个贴出来的格式就不会错，最后，--------XXX-----------也是必须的，是用来标注不同的数据，以便于程序读取。此外，每个块的顺序最好不要改变。

总体信息： 1、任务的描述，这行随便写什么，自己能看懂就行，但是必须写，不然影响整体数据读入。（1，表示对应这个部分的第1行，下同）

2、循环的次数，这个数值是对应于文献里提到的“总共使用XX步MC尝试来计算”，数值越大，计算得到的结果越稳定，当你尝试程序是否可行时，使用较小的数字比较好；当用来计算写论文时，至少100万步以上比较合适。

3、每XX步写入一次日志文件，这里的日志文件使用来记录计算过程中的一些信息，个人感觉像飞机上的黑匣子，如果计算报错，或者结果匪夷所思，就可以查看日志文件。这个数值不宜过小，因为信息太多，写入硬盘很耗时，影响整体计算时间。一般每个计算记录1～10次即可，即总步数/记录次数=这个数值。

4、每XX步写入一次crash文件，这个文件是用于记录写入文件时正在进行计算的各个参数，如果你的计算机碰到了断电，死机，或者对计算结果不满意（由于总计算步数不足导致），均可以使用以这个文件为初始，继续算下去。这样可以节省很多时间。如果你不需要，那么和总步数一致就可以。

5、每XX步写入一次config文件，这个文件是用来计算迭代过程中，每若干步时整个分子的构型，是一个很重要的文件，后期计算吸附热，换算等温线都会使用到他，建议这个数值取为总步数/100000比较合适，当然，构型记录的越多，也费时。这个数值对应于文献中提到的“每XX步输出一次构型。”

6、这个数值是config文件与crash文件生成文件名的后缀，如果设置为1，那么运算时，生成XXX.1 YYY.1。在计算多个点的时候，这个数字比较有用。 7、随机函数的种子，一般随便娶个数字，固定不变就行。

8、这个数值在这个版本的music里已经没用了，但是需要写，随便放着就OK了。 9-10、生成的config文件与crash文件的名称的基础，实际生成的文件名会在后面加上6提到的数值.

------ Atomic Types -------------------------------------------------- 5 # number of atomic types

Methane # atom type

Methane.atm # basic atom info file

Carbon # atom type

Carbon.atm # basic atom info file

Oxygen # atom type Oxygen.atm # basic atom info file

Hydrogen # atom type

Hydrogen.atm # basic atom info file

Zinc Zinc.atm

原子类型：这个部分比较简单

1、这里的数字是计算所需要使用到的所有的原子文件的个数，包括吸附质和吸附剂，这里用的是5，甲烷，碳，氢，氧，锌。

2、原子的名称，这个名称得和前面提到的原子文件里的Atom_Name保持一致，不然会报错。

3、该原子文件的文件名，程序会在前面提到路径文件里设定路径里寻找，如果找到了调用，找不到就报错。

------ Molecule Types -------------------------------------------------

2 # number of sorbate types Methane # sorbate Methane.mol # sorbate coordinates file

IRMOF1 # sorbate

IRMOF1.mol # sorbate coordinates file 和原子类型一样，不再赘述。

------ Simulation Cell Information ------------------------------------

IRMOF1 # Fundamental cell file 2, 2, 2 # No. of unit cells in x, y, z direction 1, 1, 1 # (1 = Periodic) in x, y, z 模拟晶胞信息：

1、这个得和前面介绍的吸附剂的分子文件里Molecule Name保持一致。

2、计算时采用的晶胞个数，一般是越多，计算结果的稳定性越好，但是越费时，目前主流的水平是采用2×2×2，也就是这里的设置。 3、边界周期性设置，这个概念不详细说，只有1和0两种选择，1是采用周期性边界设置，0是不采用，如果采用，会减小计算量，一般推荐在3个方向均使用。

------ Forcefield Information ------------------------------------------- BASIC SPC

atom_atom_file_UFF # atom-atom interaction file sorb_sorb_file_UFF # sorbate-sorbate interaction file

intramolecular_file # intramolecular interaction file/specification 力场信息：

1、不做解释，一般都是用BASIC这个关键词。

2、表明计算时储存的等级，至今也不明白具体含义，反正SPC就行。

3-5、分别写出前面提到的原子-原子文件、分子-分子文件，分子内部作用文件的文件名。

------ Ideal Parameters -----------------------------------------------

Ideal # Equation of State 1 # no. of sorbates Methane # Sorbate Name

理想气体参数：

1、采用的状态方程的名称，目前似乎只有Ideal可用

2、吸附质的个数，这里是1，如果计算多组分，有几个组分就填几，记住，不包括吸附剂。

3、吸附质的名称，多组分情况会日后介绍。

------ GCMC Information ----------------------------------------------- 1 # No. of iterations

298. # temperature

Ideal Parameters # Tag for the equation of state (NULL = Ideal Gas) 15 # No. of simulation points 5000 # Block size for statistics 1 # no. of sorbates -------------------------

Methane # Sorbate Name

1,1000 # pressure

Null # sitemap filename (Null = no sitemap) 3 # no of gcmc movetypes 1.0, 1.0, 1.0 # move type weights RINSERT # type of move.1 RDELETE # type of move.2 RTRANSLATE # type of move.4

0.2, 0 # Delta Translate, adjust delta option (0=NO, 1=YES) GCMC信息：

1、迭代的次数，一般就是1，这个和总体信息的迭代次数不同，设置1就行了。 2、吸附时的温度，单位是K。

3、状态方程的标签，默认就是用Ideal Parameters。

4、计算的点数，如果是1，就是计算单点吸附情况，如果是多个，就是计算吸附等温线。

5、统计的块，块分的多，有利于对迭代过程中的数值进行数值统计，有利于统计平均更加稳定，但是前提是足够多的迭代步数。

6、吸附质的个数

7、这段横线是必须的，保留就行。

8、吸附剂的名称

9、这里的压力，其实是逸度，关于逸度解释可以看我的这个帖子：http://emuch.net/bbs/viewthread.php?tid=1112118。

这里不做多数，当前面的计算点数是1时，这里就应该是一个数，如果前面的点是N个的话，这里如果给出两个数，那么程序会自动以这两个数为下限和上限，进行对数插值，生成N个对应的数据点。一般如果点太多的话，建议使用fugacity.dat文件，这个以后再讲。

10、一般就写NULL，个人也没弄清具体含义。 11、采用GCMC尝试的种类，有几种就写几。前提是必须和下面的对应上，不然会报错。

12、采用的随机插入尝试，MUSIC提供了3中插入尝试，随机插入RINSERT，偏倚插入BINSERT，构型偏倚插入CBINSERT，第一种适用于小分子吸附；第二种适合比较大的分子，如苯这些原子较多的分子，但是需要使用map，这个之后讲；第三种常用于长链烷烃，

包括直链烷烃和支链烷烃。（后两种插入方式以后讲，这次用RINSERT足够，书写格式就是上面所示的）

13、对应的删除分子的尝试，和插入方式要对应，上面用什么方式，下面就是什么方式，有RDELETE，BDELETE，CBDELETE三种。

14、平移一个分子，这里有RTRANSLATE和FTRANSLATE可用。后者是只能移动一定范围，前者是随机位置的平移。针对于具体体系使用不同的尝试方式。

15、移动的幅度，是否矫正。一般移动的幅度不易过大，不宜收敛稳定，移动幅度过小，费时。0.2比较合适，随机平移时，一般不采用矫正移动。

因为甲烷被假设为一个伪原子，所以不存在旋转。这里一共是3中MC尝试。

------ Configuration Initialization ------------------------------------- Methane # Sorbate_Type GCMC NULL

IRMOF1 # Sorbate_Type FIXED NULL 构型初始化：

这里是写入构型初始化信息，需要罗列每个物质。

1、吸附质的名称

2、如果是全新开始算，吸附质对应的就填写GCMC，表示用GCMC插入获得构型，如果接着之前的结果继续算（由于断电等因素），这里就填写

RESTARTFILE XXX

XXX前面提到的crash文件的文件名。

3、吸附剂的名称

4、如果把吸附剂认为是刚性的结果，就选择FIXED，一般选择FIXED。 -------- Main Datafile Information --------

Energy, position, pair_energy # contents of datafile

写在config文件里的内容，对于MC计算，包括能量，原子的位置坐标等

好了，计算这个实例所使用到的文件基本介绍完了，如果你也想尝试，我可以把这写文件压缩上传。

下面我们开始计算，打开一个终端，进入到控制文件和gcmc.exe共同存在的文件夹下,下载得到的分子分子作用文件还得改一下，上下都改成一样的。路径也得改成自己机器的路径。

>> $ source set_path

>> $ ./gcmc.exe gcmc.Methane.IRMOF1 >out.log

这样程序就会开始运行。贴上我计算的结果，和文献中的数据比较情况。

298K下IRMOF-1吸附甲烷的吸附等温线比较，Snurr就是我们说的Snurr课题组的文献，Yang and Zhong是北京化工大学的阳庆元和仲崇立课题组。

文献参考如下：

(1) Yang, Q. Y.; Zhong, C. L. J. Phys. Chem. B 2006, 110, 17776. (2) Duren, T.; Snurr, R. Q. J. Phys. Chem. B 2004, 108, 15703.

至于我说的这个实例的文件包，大家可以到这个网盘上下载 http://www.163pan.com/files/30700020y.html

文件名称:

gcmc.zip Quote:

Originally posted by chengxuan at 2010-03-09 14:39:03: 天啊，写的太详细了，太好了。忍不住顶一下。

请教你一个小问题哦！你提到在断电时可以用crash文件接着算，我这就断过两次电，结果都是从第一步开始算的，那样很浪费时间，今天看到你这个贴，真的很开心，又学到 ...

恩，改成这样，就可以从上次计算的结束时开始计算下去，比较省事。 Methane # Sorbate_Type RESTARTFILE XXX

IRMOF1 # Sorbate_Type FIXED NULL

Extra section in control file

不少虫子都问我如何用music求snapshot，那么就借着这个机会，将music控制文件中的extra部分介绍一下，当然，这部分有很多，我主要介绍和MC有关的部分。

A、movie information

这部分是用来获得MC过程中截图的，如果你只取少量的几个截图，并且一次只展示一个瞬间的构型，那么就是文献中所说的snapshot，如果你将几百几千个截图放在一起展示，那么就是文献中通常展示的density distribution图，其实原理是一样的，如果你做的吸附剂吸附量比较大，某一瞬间的截图中包含的吸附质分子比较多，便于研究优先吸附位等规律的时候，就用snapshot，如果吸附量比较小，假设一个吸附剂总共就有两个吸附质分子（在吸附大分子时尤为明显），那一个瞬间的截图肯定是无法研究规律的，那么就需要使用density distribution图，通过统计大量的截图中吸附质的位置，来研究这些性质。

如果你要计算这个，那么就不能计算一个等温线，建议在压力那行只写上一个数字，做fixed pressure计算。其他的不用改变，在gcmc.ctr文件最后放加上这么一块内容： ------ Movie Information ---------------------------------------------- movie.xyz # Movie filename

10000, 2000 # Starting step, ending step 100 # Steps between frames Yes # Include zeolite in movie

1, 1, 1 # Number of unit cells to dump in x, y, z directions 首先这个标题还是一样不可或缺，不然系统无法识别

1、Movie文件名，这个随便定义，只要自己识别就好，后缀其实无所谓，但是为了在Windows在使用可视化软件方便，推荐用xyz后缀。

2、开始做截图的起始和结束步数，一般不推荐从第一步开始，因为还没有平衡，推荐从一半后开始，取多少步随便你，哪怕是10001到10002也行。 3、每多少步写入movie.xyz文件一次，这里就是做density distribution图的人要用到的了，你所获得的截图数是(ending-starting)/steps between frames，这和前面写入con文件的道理是一样的。 4、Movie文件里是否包含骨架原子的坐标，写yes/no，一般做snapshot选择yes，做density

本文来源：https://www.bwwdw.com/article/v0tg.html