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

[工作范文]模拟心得MATERIAL STUDIO 中SORPTION 模拟心得 MATERIAL STUDIO 中SORPTION 第一个课题是模拟金属有机框架和共价有机框架吸附CO以及分离CO/CH,使用的软件是224 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); %+++++++++++++++++++++++++++++++++++++++++++++ %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; 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); %+++++++++++++++++++++++++++++++++++++++++++++ %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(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里面手动改成文献上那样，因 3为有时候得到的结构会很不规则。电荷是使用DMOL计算得到的，但是对于某些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 这个任务，CO的CONNOLY RADIUS是0.165NM,可以再变化VDW SCALE FACTOR的时2 候，得到一些不同的自由体积。文献中 specific area 和 free volume的单位与MS得到的不同，在写文章的时候，需要转化一下。 3沸石的电荷一般用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构型。 图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) 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 real distance,search_step,dis1,dis2,dis3,lbond,charge(natom) ****************************************************************** * define global variables ****************************************************************** common charge common xc,yc,zc ****************************************************************** * 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(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) 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 ********************************************************* * 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) 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 * * 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(Tatom)='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),atomname(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+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.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+ * & (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=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(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), * & 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 * 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) * endif *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) * 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)=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+ * & (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)=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 ******************************************** * 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,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 * 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' * atomname(Tatom)='Si' * charge(Tatom)=3.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),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) 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_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' * 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' 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), & 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, & 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_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' * 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' atomname(Tatom)='C' charge(Tatom)=7.000 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) close(160) ****************************************************************** * 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, & 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_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' * 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' atomname(Tatom)='C' charge(Tatom)=8.000 open(170,file='amino.car',access='append') write(170,888)a(Tatom),xc(Tatom), & yc(Tatom),zc(Tatom), & fx(Tatom),occupation(Tatom), & fft(Tatom),atomname(Tatom), & charge(Tatom) close(170) ****************************************************************** * 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, & xfinal,yfinal,zfinal) natom_add=natom_add+1 Tatom=natom0+natom_add * NNT=NNT+1 * NN_C=INT(NNT/99) * NN_S=NNT-99*NN_C * 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' * a(Tatom)='N'//'NN_S'//'NN_CC' a(Tatom)='N' xc(Tatom)=xfinal yc(Tatom)=yfinal zc(Tatom)=zfinal fx(Tatom)='XXXX' occupation(Tatom)=1 fft(Tatom)='N_3' atomname(Tatom)='N' charge(Tatom)=9.000 open(180,file='amino.car',access='append') write(180,888)a(Tatom),xc(Tatom), & yc(Tatom),zc(Tatom), & fx(Tatom),occupation(Tatom), & fft(Tatom),atomname(Tatom), & charge(Tatom) close(180) 45 continue **************************************************** * rewrite the atom serial. **************************************************** open(46,file='Satom.car',status='old') do 47 i=natom0+1,Tatom read(46,1888)a(i) 47 continue close(46) 1888 format(A5) * endif * *120 continue *110 continue *100 continue *************************************************** * 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) **************************************************** **************************************************** end *&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& * 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) ******************************************************************* * declare the global varialbe ******************************************************************* real charge(natom) double precision xc(natom),yc(natom),zc(natom) common charge common xc,yc,zc kjishu=0 do 50 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) endif 50 continue ****************************************************************** * write out the hydroxy oxygen list ****************************************************************** open(51,file='HO_list.car',access='append') do 52 i=1,kjishu write(51,999)int(Ohydroxy(i,1)),Ohydroxy(i,2), & Ohydroxy(i,3),Ohydroxy(i,4),Ohydroxy_SN(i) 52 continue write(51,*)'***********************************************' write(51,*)'' close(51) 999 Format(I4,3X,F12.9,2X,F13.9,3X,F12.9,3X,I4) endsubroutine *&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& *&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& *&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& * add new atom *&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine addatom(xtemp,ytemp,ztemp,lbond,Tatom,Temp_SN, & xfinal,yfinal,zfinal) integer natom parameter(natom=15000) real search_step,lbond,minlbond,maxlbond double precision xtemp,ytemp,ztemp,targetx,targety,targetz double precision dis1,dis2,xc(natom),yc(natom),zc(natom) double precision xfinal,yfinal,zfinal double precision center(14,2),xcenter,ycenter,lc,lct integer ia,ib,ic,ja,jb,jc integer Tatom,Temp_SN ***************************************************************** * difine global variable ***************************************************************** common xc,yc,zc ***************************************************************** ***************************************************************** open(3000,file='center.txt',status='old') do 3010 i=1,14 read(3000,*)(center(i,j),j=1,2) 3010 continue close(3000) * open(3011,file='center-1.txt',access='append') * do 3012 i=1,14 * write(3011,*)(center(i,j),j=1,2) *3012 continue * close(3011) lc=1000 xcenter=0 ycenter=0 do 3020 i=1,14 lct=sqrt((xtemp-center(i,1))**2+(ytemp-center(i,2))**2) if (lct.lt.lc)then lc=lct xcenter=center(i,1) ycenter=center(i,2) endif 3020 continue if(ytemp.le.ycenter)then ja=20 jb=0 jc=-1 jd=j else ja=-20 jb=0 jc=1 jd=j endif if(xtemp.le.xcenter-5)then ia=20 ib=0 ic=-1 minlbond=lbond-0.1 maxlbond=lbond+0.1 search_step=0.1 do 1000 i=ia,ib,ic do 1010 j=ja,jb,jc do 1020 k=-20,20 targetx=xtemp+search_step*i targety=ytemp+search_step*j targetz=ztemp+search_step*k dis1=sqrt((targetx-xtemp)**2+ & (targety-ytemp)**2+ & (targetz-ztemp)**2) if(dis1.gt.minlbond.and.dis1.lt.maxlbond)then do 1030 im=1,Tatom if(im.ne.Temp_SN)then dis2=sqrt((targetx-xc(im))**2+ & (targety-yc(im))**2+ & (targetz-zc(im))**2) if(dis2.lt.3.0)then goto 1020 endif endif * if(im.ne.Temp_SN.and.dis2.lt.3.0)then * goto 1020 * endif 1030 continue goto 1040 endif 1020 continue 1010 continue 1000 continue 1040 xfinal=targetx yfinal=targety zfinal=targetz elseif(xtemp.gt.xcenter+5)then ia=-20 ib=0 ic=1 minlbond=lbond-0.1 maxlbond=lbond+0.1 search_step=0.1 do 1001 i=ia,ib,ic do 1011 j=ja,jb,jc do 1021 k=-20,20 targetx=xtemp+search_step*i targety=ytemp+search_step*j targetz=ztemp+search_step*k dis1=sqrt((targetx-xtemp)**2+ & (targety-ytemp)**2+ & (targetz-ztemp)**2) if(dis1.gt.minlbond.and.dis1.lt.maxlbond)then do 1031 im=1,Tatom if(im.ne.Temp_SN)then dis2=sqrt((targetx-xc(im))**2+ & (targety-yc(im))**2+ & (targetz-zc(im))**2) if(dis2.lt.3.0)then goto 1021 endif endif * if(dis1.gt.minlbond.and.dis1.lt.maxlbond)then * do 1031 im=1,Tatom * dis2=sqrt((targetx-xc(im))**2+ * & (targety-yc(im))**2+ * & (targetz-zc(im))**2) * if(im.ne.Temp_SN.and.dis2.lt.3.0)then * goto 1021 * endif 1031 continue goto 1041 endif 1021 continue 1011 continue 1001 continue 1041 xfinal=targetx yfinal=targety zfinal=targetz else ia=20 ib=-20 ic=-1 minlbond=lbond-0.1 maxlbond=lbond+0.1 search_step=0.1 do 1002 j=ja,jb,jc do 1012 i=ia,ib,ic do 1022 k=-20,20 targetx=xtemp+search_step*i targety=ytemp+search_step*j targetz=ztemp+search_step*k dis1=sqrt((targetx-xtemp)**2+ & (targety-ytemp)**2+ & (targetz-ztemp)**2) if(dis1.gt.minlbond.and.dis1.lt.maxlbond)then do 1032 im=1,Tatom if(im.ne.Temp_SN)then dis2=sqrt((targetx-xc(im))**2+ & (targety-yc(im))**2+ & (targetz-zc(im))**2) if(dis2.lt.3.0)then goto 1022 endif endif * if(dis1.gt.minlbond.and.dis1.lt.maxlbond)then * do 1032 im=1,Tatom * dis2=sqrt((targetx-xc(im))**2+ * & (targety-yc(im))**2+ * & (targetz-zc(im))**2) * if(im.ne.Temp_SN.and.dis2.lt.3.0)then * goto 1022 * endif 1032 continue goto 1042 endif 1022 continue 1012 continue 1002 continue 1042 xfinal=targetx yfinal=targety zfinal=targetz endif endsubroutine *&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& &&&&&&&&&&& *&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& &&&&&&&&&&& *&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& * Random Number Generator *&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& FUNCTION RAN2(IDUM) 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 iy=iv(j)-idum2 iv(j)=idum if(iy.lt.1)iy=iy+IMM1 ran2=min(AM*iy,RNMX) return END *&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& *&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 次，得到的CAR文件输入MS各段程序大概信息已经在程序内部给出，调用添加原子的子程序3 就可以自动加氢了。这个程序中使用的MCM-41模型还是单个晶胞4孔模型(由3*7*2超晶胞挖出来的)，对于单孔模型的MCM-41需要自己改写，在CENTER.TXT那里修改，后面的搜索步骤也要相应的修改。其他种类的沸石例如MFI，LTA等等加胺基也是按照这样做的。对于怎么计算孔径分布，需要到处CAR文件中所有的H原子，其源程序是这样的: integer hno parameter(hno=2356) 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) 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), & 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), & occupation(imm),fft(imm),natom(imm),charge(imm) charge(imm)=2.0 endif 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 if(jmm.lt.415)then goto 40 else goto 50 endif 50 close(30) 888 format(A5,3X,F12.9,2X,F13.9,3X,F12.9,1X,A4,1X,I1,6X,A5,3X,A1,F8.3) end *&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& * Random Number Generator *&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& FUNCTION RAN2(IDUM) 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 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)，因为它有非商业版本下载免费使用，只要你不是用于商业用途，就可以，但是你也不会得到 技术支持。网址如下 只要按照要求填写一个可用的邮箱，就可以免费下载。邮箱是用来接收许可证文件和序列号的，两者只要使用一个就行。 这里我下载的是11.1.059版本，虽然比较大，但是支持32位和64位。注意的是，如果你使用的是AMD的处理器，下载的时候得注意选择正确合适的版本。 下载好后，我们开始安装ifort。首先解压缩，因为下载下来的是一个tgz的包，可以使用命令 $ tar -zvxf l_cprof_p_11.1.059_ia32.tgz ($表示一般用户权限，#表示需要root权限) 进行解压缩。然后进入该文件夹，使用命令 $ ./install.sh 开始安装操作，安装过程比较简单，中间会提示你输入序列号，收到的邮件里有。安装好了后，为了确保我们是否安装成功，我们检验下。首先在终端里输入 $ source /opt/intel/Compiler/11.1/059/bin/ia32/ifortvars_ia32.sh (我使用的是32位的系统，具体情况自己修改即可)然后输入 $ which ifort 如果正确，系统会显示ifort的安装路径，如果有兴趣的朋友可以对一个小程序进行试编译。 安装好后，我们就开始对music进行编译 下载下来的music也是tgz的压缩包，我们采用同样的方式进行解压缩，解压好后，我们将其放在用户的家目录下/home/XXX/ Music目录下有以下几个目录/ctrlfiles /drivers /src /tests /utils还有一个GPL的授权书。其中/ctrlfiles主要存放control文件的 模板 个人简介word模板免费下载关于员工迟到处罚通告模板康奈尔office模板下载康奈尔 笔记本 模板 下载软件方案模板免费下载 (control文件干嘛后面讲。);/drivers主要存放运算模块的F90文件，包括了gcmc和MD等;/src下包含了程序中的源代码，格式基本为F90，也有少量的F;/tests存放的是大量的实例，用于检验编译的可执行文件是否正确，以及对使用方法的验证;/utils包含了一些实用的小程序，辅助使用music。 这里我们开始编译，首先确保打开了ifort编译器，每次打开终端，如果要编译，最好确认下。如果没有 $ source /opt/intel/Compiler/11.1/059/bin/ia32/ifortvars_ia32.sh 即可，然后开始。 $ cd music-4.0/src $ cp ../drivers/music_gcmc.F90 music.F90 这里需要修改下makemake文件和Makefile文件，这两个文件都在/src下，红色为修改部分。具体原因不做解释。 ===================makemake 文件================ $compiler = getoldvalue('Makefile','FC'); 44行 if ($compiler eq '0') {$compiler = 'ifort'; } =============================================== 改好后，输入 $ ./makemake 再修改MakeFile文件 ===================Makefile 文件================ # Executable name 4行 CMD = gcmc.exe ifeq ($(FORTRAN_COMPILER),) 27行 FC = ifort else FC= $(FORTRAN_COMPILER) endif ifeq ($(FC),ifort) 60行 F90FLAGS = -g -debug all -check noarg_temp_created # -O3 # -DSPME # untested options ============================================ 修改好后 $ make 这样等待编译好后，就会在/music-4.0/下生成一个gcmc.exe文件，这个就是用用来计算的执行文件。 编译好后，还要记得删除编译过程中删除的.mod和.o文件，因为这些文件的存在，会影响其他模块的编译。因为每个F90文件都包含了不同模块的代码，每次编译都运行了其中一部分代码，如果不删除，下次编译另外一个模块时，这些存在的.mod和.o文件会使ifort不再对相应的F90文件进行编译，因为它认为编译过了，这样就会导致其他模块编译错误。 运行命令 $ rm *.mod *.o 即可删除对应的文件。 这样一个模块就编译好了，不过对于music-4.0来说，还有部分模块没有加工制作好，比如NVTMC，会出现编译错误，我写信给Dr. Chempath S.和Dr. Düren T.(他们都参与了music的编写制作，并使用music发表过高质量的文章) 他们都说可以使用GCMC来替代，只要取消插入和删除两种尝试即可。 好了，编译方面就写这么多，下面会一个个介绍计算所需要的每个文件。 文件 MuSiC程序作计算时，是需要一个个不同类型的文件为编译好的可执行文件提供数据的，这里以GCMC计算为例，计算过程中最少需要原子文件、分子文件、原子-原子作用文件、分子-分子作用文件、分子内部作用文件，控制文件以及环境变量文件，如果还有更多的需求，需要添加更多的文件。不过不要害怕，如果熟悉了这些文件的作用，你会发现一点都不难(用过towhee，这个程序里使用的文件让人无语，呵呵)下面我们一个个来解析。 原子文件:主要是提供你计算所需要调用的原子的信息，包括原子名，原子量等，我们以C原子为例: ##### Basic Atom Information Atom_Name: Carbon Atom_Symbol:C Atom_SS_Charge:0.0 Atom_SZ_Charge:0.0 Atom_Mass: 12.0 Atom_Valency: 4 第一行是原子的名字，这个名字很重要，是后面MuSiC程序计算的时候所需要调用的名字，所以最好不要随便定义，以该元素的英文单词为佳，以便以后使用时便于识别。这里 注意的是，如果对于不同物质里的同一个元素，建议使用不同的元素名，这样以便于区分使 用，好处以后讲。比如CO2中的C，建议命名为Carbon_CO2，其他以此类推。 第二行是该元素的符号，建议使用元素周期表中的标准写法，因为是为了以后模拟的 结果可以被一些分子结构显示软件正确识别。 第三行和第四行分别是关于该元素的电荷，具体是什么含义不详，不过在MuSiC程序 计算中不会被使用到，但是写还是得写在上面。 第五行是原子的原子量。第六行是该元素的最大化合价。 最后，该文件保存为XXX.atm的格式，虽说Linux的文件不讲究什么后缀，但是加 上.atm便于识别，可能也便于程序调用。 分子文件:分子文件是计算的基础，是将我们事先写好的原子组成一个可以计算的分 子，这里以正丁烷举例: ##### Basic Molecule Information Molecule Name: Butane Coord_Info: Listed Cartesian 4 # number of atoms in molecule 1 0.000000 0.000000 0.000000 Methyl 0.0 0 0 2 1.268427 -0.855565 0.000000 Methylene 0.0 0 0 3 2.536855 0.000000 0.000000 Methylene 0.0 0 0 4 3.805282 -0.855565 0.000000 Methyl 0.0 0 0 Connect_Info: Generate 2 # number of connection types listed Methyl Methylene 1.530000 Methylene Methylene 1.530000 Connect_Info: Generate ##### Intramolecular forces 2 # number of bond types to generate Methyl Methylene 1000.00 1.530000 Methylene Methylene 1000.00 1.530000 Bond_Constraints: Listed Evans Fast 3 #number of bond constraints (0.0d0 bondlength ==> to be calculated from above) 1 2 1.530000 2 3 1.530000 3 4 1.530000 Bond_Bending: Listed Harmonica Fast 2 # number of bending centers (angle = 0.0d0 ==> to be calculated from above) 1 2 3 124.400000 112.000000 #atoms, ktheta (kcal/mol), kthetaeq (degrees) 2 3 4 124.400000 112.000000 #atoms, ktheta (kcal/mol), kthetaeq (degrees) Bond\_Torsion: Listed Cosexpansion Fast 1 # number of torsion angles 1 2 3 4 2.217500 2.905100 -3.135600 -0.731200 6.271200 -7.527100 # end of butane.mol 第一行还是分子名，建议使用英文命名，便于识别记忆。 第二行是分子的坐标信息，这里有多种类型可以选择，这里是采用了列出了笛卡尔坐标 的方式，也是最常用的一种方式。后面还可以添加不同的参数以表示优化分子的计算，比如rigid多用于刚性的分子，比如CH和H等，或者branched多用于支链烷烃，等等。 42 第三行列出该分子所包含原子的个数。这里是4。 第四行到第七行分别依次列出了正丁烷的四个原子的xyz坐标、原子类型、电荷、集和类型，其中集和类型一般为0，我也不知道具体含义，不过似乎不影响计算，也没有看到过不为0的情况，所以暂时搁置。这里有人可能会问，正丁烷怎么是四个原子，没错，是四个，不过是四个“伪原子”，这个是基于TraPPE联合原子力场，将甲烷、甲基、亚甲基等模拟为一个原子，这样整个分子结构被简化，计算效率大大提高，但是计算精度能够保证和全原子模型相差无几。 后面的部分是关于分子内部的作用力，如果一个分子被认为是刚性的，那么一个文件只要写到刚才就OK了，如果像正丁烷这样的分子，考虑为刚性是不合理的，吸附过程中肯定为弯曲，伸缩等等，所以需要添加这部分参数。首先是键的长度，其次是键伸缩的参数，这里键伸缩采取了虎克定律来表达，1000和1.53都是参数，具体这些参数在TraPPE联合原子力场里都是可以查到的，所以不做详细介绍，后面的键角变化，键的扭转等项均是类似，具体可以参照文献:Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1998, 102, 2569. 最近比较忙，所以一直没有更新，现在抓紧时间更新吧。前面讲了原子文件和分子文件，这个是计算吸附的最基础的文件，所有的计算都是建立在这个基础上的。后面所用到的文件，就是为计算提供一些数据的文件，这些文件如果拿出来单独解释不太合适，所以还是放入例子中解释比较好，这样容易记住。 这里选用了一个最简单的例子，金属有机骨架IRMOF-1在298K下吸附甲烷的吸附等温线的计算实例，我们会一步步分析每个文件的作用，并且运行计算后，将计算结果和文献结果进行对比。(IRMOF-1如果你不知道是什么不用管，一种吸附剂而已) 计算中用到的输入文件主要是两种，第一种就是关于力场的参数文件，第二种就是计算使用的控制文件。第一种文件包括原子-原子作用文件;分子-分子作用文件;分子内部作用文件: 原子-原子文件:先把实例中的该文件列出来: Methane Methane LJ SIG@3.730 EPS@148.0 HICUT@12.8 Methane Carbon LJ SIG@3.581 EPS@88.462 HICUT@12.8 Methane Hydrogen LJ SIG@3.151 EPS@57.265 HICUT@12.8 Methane Oxygen LJ SIG@3.424 EPS@66.872 HICUT@12.8 Methane Zinc LJ SIG@3.096 EPS@96.099 HICUT@12.8 可以看出来，这个文件主要是用来描述原子与原子之间的交互作用力，一般分为Lenard-Jones和Buckingham两种形式来描述，这里由于选择了使用UFF力场描述IRMOF-1，TraPPE-UA力场描述甲烷分子，所以均使用LJ这个参数，后面的数值分别是各个力场提供的参数经过LB混合规则计算过以后获得的，最后的HICUT表示了截断半径。相信了解LJ的人，都不会对这个陌生。 这里可以看出来，原子-原子文件的书写最常用的形式是: Atom Atom LJ SIG@XXX EPS@XXX HICUT@XXX 如果两个原子之间不存在交互作用，则是: Atom Atom OFF 由于甲烷被描述为一个原子的形式，所以是整体电中性的，因此在这里没有考虑库伦电荷作用，因此静电力的书写方式暂时不做介绍。 分子-分子作用文件:还是先把实例中的该文件列出来: Methane Methane NCOUL BASIC LJ FAST Methane Methane COUL OFF Methane IRMOF1 NCOUL BASIC LJ FAST Methane IRMOF1 COUL OFF IRMOF1 IRMOF1 NCOUL OFF IRMOF1 IRMOF1 COUL OFF 分子-分子文件主要是用来描述分子之间的相互作用力，比如这个实例中，存在VDW 力作用的只有甲烷-甲烷，甲烷-IRMOF-1两种，由于吸附剂自身之间不会存在作用，因此不 考虑。在这里，NCOUL表示非库仑力，即VDW力。COUL表示库仑作用力，在这个实例 中也是不考虑的，因此都是OFF。 分子内作用文件: Intra: Methane Intra: IRMOF1 这个文件是用来描述吸附剂与吸附质分子内部的作用，比如，分子的弯曲，扭 转等等，在这里，由于IRMOF-1和甲烷都被近似为刚性分子，因此，后面的参数都被简化 了。再比如正丁烷，后面会加上: Intra: Butane BENDING TORSION 这里就标出了正丁烷的分子内部的作用，包括键角的变化和整个分子的扭转。 今天收到一个虫友的短信息，问我说按照我写的东西，将music运行在fedora 12下， 也编译出来gcmc.exe这个执行文件，但是不知道如何执行，连wine都用上了，这里就让我 汗颜了，更新不够快，后面的还没讲呢，今天赶紧把这个例子完成，这样想实践的朋友就可 以先算起来了。 环境变量文件:计算时，gcmc会调用到前面提到的原子文件和分子文件，因此这里必 须设定环境变量，我们写环境变量如下: #!/bin/sh export ATOMSDIR=/home/tina/tutorial/gcmc/atoms export MOLSDIR=/home/tina/tutorial/gcmc/molecules export PMAPDIR=/home/tina/tutorial/gcmc/pmaps 这里可以看出，这是一个BASH文件，需要将分子文件、原子文件的路径给明，这里 也许你会发现PMAPDIR，这个是用来制定Pmap文件的位置，之前没讲到，但是暂时放到 这里，以后会详细讲明的。这里，路径一定要使用绝对路径，避免发生错误。 控制文件:这个文件是整个计算过程中的重中之重，因为它的信息量最大，包括了计 算的各方面设置。先按次序把这个文件中的内容分块贴出来，挨个儿分析。 ------ General Information ------------------------------------------ Methane molecule in IRMOF1 1000000 # No. of iterations 100000 # No. of steps between writes to output/log file 100000 # No. of steps between writes to crash file 5000 # No. of steps between writes to config. file 1 # Start numbering simulations from . 10283 # Iseed 3 # specifies contents of config file, 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、这里的压力，其实是逸度，关于逸度解释可以看我的这个帖子: 。 这里不做多数，当前面的计算点数是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. 至于我说的这个实例的文件包，大家可以到这个网盘上下载 文件名称: 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 distribution一般选择no，因为输出多个截图，骨架同样也输出多次，后期不好处理。 5、Movie文件所包含的晶胞个数，如果你算的是单胞，那么肯定是1,1,1，如果是多胞，那么这个根据你自己研究的情况而定。 输出的这个movie文件可以直接在MS、rasmol、gopenmol等可视化软件中显示。 B、中止条件 如果你做GCMC，一般这个不会用到，但是如果计算NVTMC，这个就很有必要了，NVTMC只吸附在骨架里的吸附质分子不变的情况下进行MC尝试移动，但是我们如何让骨架内预先有我们想要的吸附质分子个数呢。如果用GCMC获得，一般不会是整数，那么就要使用到这个部分，同样，也是只算一个压力，这个压力最好比较大，可以保证你的吸附剂可以吸附到你想要的数量，比如你想让IRMOF1吸附30个甲烷分子，如果你只输入100kPa显然这个计算是永远不可能结束的，那么就需要你将这个压力提高，比如5000kPa。加入这部分内容如下: ----------------------- Stop Criteria ------------------------------ Nmolecs, Methane, 64, EQ # stop simulation if Nmolecs of Methane is EQUAL to 64 ----------------------------------------------------------------------# this empty line is needed 这个很好理解，当甲烷的个数等于64的时候即停止gcmc运算。这里好像还可以填GE和LE等，就是大于等于/小于等于，但是似乎没什么用。 这个算好后，就可以获得一个名称为finalconfig.xyz的文件，里面的构型就包含了64个甲烷分子，后续你就可以做NVTMC了。 这里需要说明的是，music4.0的NVTMC模块还没有开发完毕，所以你是不可能将NVTMC编译成功的，但是我们要计算NVTMC怎么办呢，有替代方法，就是用GCMC，两者的不同就是系综不同，我们写一个类似于gcmc.ctr，但是添加的MC trial尝试中不要含有insert和delete，比如只有translate和rotate，那就在MC trail上和NVTMC一致了，其次由于运行的还是gcmc.exe，所以在ctr文件里，压力还是需要填上一个数字，但是这个数字在运行时不参与计算，所以随便填个就行。 呵呵，又翻页了，你说的这种是氢气的三位模型，是一种能够考虑H2弱的四极距的模型，分子文件如下: Molecule_Name: Hydrogen_3s Coord_Info: Listed Cartersian Rigid 3 # Number of atom in molecule 1 0.0000 0.0000 XXXXX Hydrogen 0.468 0 0 2 0.0000 0.0000 -XXXXX Hydrogen 0.468 0 0 3 0.0000 0.0000 0.0000 PsedoAtom -0.936 0 0 Molecule_DOF: 6 这里的PsedoAtom就是自己建立的虚原子，这种虚原子跟一般原子的原子文件一样，只是原子量设定为0。因为有电荷，所以要建立emap来加速，但是没有pmap，因为它没有LJ作用。 XXXXX为你建立的H2分子中两个H原子之间键的一半长度，这个文献里应该有，自己写上去就行了。 离子液体的MD模拟(GROMACS) 我使用的是GROMACS，由于PRODRG网站不能自动生成离子液体的GRO和TOP以 及ITP文件，需要自己建立ITP文件，给一个例子: ; FORCEFIELD of Ionic liquids V1.0 ; SO2, CO2, etc. ; LJ parameters from OPLSAA FF ; Intramolecular force field constants from OPLSAA ; 2004-07-16, by Wang Yong, Pan Haihua ; This force field uses a format that requires Gromacs 3.1.4 or later. ; ; References for the OPLS-AA force field: ; ; W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, ; J. Am. Chem. Soc. 118, 11225-11236 (1996). ; W. L. Jorgensen and N. A. McDonald, Theochem 424, 145-155 (1998). ; W. L. Jorgensen and N. A. McDonald, J. Phys. Chem. B 102, 8049-8059 (1998). ; R. C. Rizzo and W. L. Jorgensen, J. Am. Chem. Soc. 121, 4827-4836 (1999). ; M. L. Price, D. Ostrovsky, and W. L. Jorgensen, J. Comp. Chem. (2001). ; E. K. Watkins and W. L. Jorgensen, J. Phys. Chem. A 105, 4118-4125 (2001). ; G. A. Kaminski, R.A. Friesner, J.Tirado-Rives and W.L. Jorgensen, J. Phys. Chem. B 105, 6474 (2001). ; [ defaults ] ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ 1 3 yes 0.5 0.5 [ atomtypes ] ; name bond_type mass charge ptype sigma epsilon opls_302 CA 12.01100 0.000 A 2.25000e-01 2.09200e-01 ; C=NH2 opls_300 N2 14.00670 -0.000 A 3.25000e-01 7.11280e-01 ; C=NH2 opls_301 H3 1.00800 0.000 A 0.00000e+00 0.00000e+00 ; C=NH2 opls_305 CT 12.01100 0.000 A 3.50000e-01 2.76144e-01 ; CH3-NR2 opls_911 HC 1.00800 0.000 A 2.50000e-01 6.27600e-02 ; CH3-NR2 opls_271 C_3 12.01100 -0.100 A 3.75000e-01 4.39320e-01 ; C in COO- carboxylate opls_272 O2 15.99940 -0.800 A 2.96000e-01 8.78640e-01 ; O: O in COO- carboxylate,peptide terminus opls_275 CT 12.01100 -0.160 A 3.50000e-01 2.76144e-01 ; C: CH, carboxylate ion opls_282 HC 1.00800 0.060 A 2.42000e-01 6.27600e-02 ; AA H on C-alpha in ketone & aldehyde ? opls_154 OH 15.99940 -1.300 A 3.12000e-01 7.11280e-01 ; all-atom O: mono alcohols opls_155 HO 1.00800 0.300 A 0.00000e+00 0.00000e+00 ; all-atom H(O): mono alcohols, OP(=O)2 JACS_000 C 12.01100 0.700 A 2.80000e-01 0.22260e+00 ; JACS,2004,126,5300. JACS_001 O 15.99940 -0.350 A 3.05000e-01 0.65940e+00 ; JACS,2004,126,5300. ;;;;;; TETRAMETHYLGUANIDINE TMG [ moleculetype ] ; molname nrexcl TMG 3 [ atoms ] ; nr type resnr residue atom cgnr charge mass typeB chargeB massB 1 opls_300 1 TMG NT1 1 -0.826 14.0067 ; qtot -0.826 N2 2 opls_301 1 TMG HT2 1 0.420 1.0080 ; qtot -0.406 H3 3 opls_301 1 TMG HT3 1 0.420 1.0080 ; qtot 0.014 H3 4 opls_302 1 TMG CT4 2 0.511 12.0110 ; qtot 0.525 CA 5 opls_300 1 TMG NT5 2 -0.046 14.0067 ; qtot 0.479 N2 6 opls_300 1 TMG NT6 2 -0.046 14.0067 ; qtot 0.433 N2 7 opls_305 1 TMG CT7 3 -0.218 12.0110 ; qtot 0.215 CT 8 opls_911 1 TMG HT8 3 0.108 1.0080 ; qtot 0.323 HC 9 opls_911 1 TMG HT9 3 0.142 1.0080 ; qtot 0.465 HC 10 opls_911 1 TMG HT10 3 0.119 1.0080 ; qtot 0.584 HC 11 opls_305 1 TMG CT11 4 -0.256 12.0110 ; qtot 0.328 CT 12 opls_911 1 TMG HT12 4 0.107 1.0080 ; qtot 0.435 HC 13 opls_911 1 TMG HT13 4 0.143 1.0080 ; qtot 0.578 HC 14 opls_911 1 TMG HT14 4 0.138 1.0080 ; qtot 0.716 HC 15 opls_305 1 TMG CT15 5 -0.218 12.0110 ; qtot 0.498 CT 16 opls_911 1 TMG HT16 5 0.108 1.0080 ; qtot 0.606 HC 17 opls_911 1 TMG HT17 5 0.142 1.0080 ; qtot 0.748 HC 18 opls_911 1 TMG HT18 5 0.119 1.0080 ; qtot 0.867 HC 19 opls_305 1 TMG CT19 6 -0.256 12.0110 ; qtot 0.611 CT 20 opls_911 1 TMG HT20 6 0.107 1.0080 ; qtot 0.718 HC 21 opls_911 1 TMG HT21 6 0.144 1.0080 ; qtot 0.862 HC 22 opls_911 1 TMG HT22 6 0.138 1.0080 ; qtot 1.000 HC [bonds] ; i j func b0 kb 1 2 1 0.1007 363171.2 ; N2-H3 1 3 1 0.1007 363171.2 ; N2-H3 1 4 1 0.1347 402500.8 ; N2-CA 4 5 1 0.1344 402500.8 ; N2-CA 4 6 1 0.1344 402500.8 ; N2-CA 5 7 1 0.1469 282001.6 ; N2-CT 5 11 1 0.1471 282001.6 ; N2-CT 6 15 1 0.1469 282001.6 ; N2-CT 6 19 1 0.1471 282001.6 ; N2-CT 7 8 1 0.1089 284512.0 ; CT-HC 7 9 1 0.1090 284512.0 ; CT-HC 7 10 1 0.1095 284512.0 ; CT-HC 11 12 1 0.1087 284512.0 ; CT-HC 11 13 1 0.1092 284512.0 ; CT-HC 11 14 1 0.1092 284512.0 ; CT-HC 15 16 1 0.1089 284512.0 ; CT-HC 15 17 1 0.1090 284512.0 ; CT-HC 15 18 1 0.1095 284512.0 ; CT-HC 19 20 1 0.1087 284512.0 ; CT-HC 19 21 1 0.1092 284512.0 ; CT-HC 19 22 1 0.1092 284512.0 ; CT-HC [ pairs ] ; ai aj funct c0 c1 c2 c3 1 7 1 1 11 1 1 15 1 1 19 1 2 5 1 2 6 1 3 5 1 3 6 1 4 8 1 4 9 1 4 10 1 4 12 1 4 13 1 4 14 1 4 16 1 4 17 1 4 18 1 4 20 1 4 21 1 4 22 1 5 15 1 5 19 1 6 7 1 6 11 1 7 12 1 7 13 1 7 14 1 8 11 1 9 11 1 10 11 1 15 20 1 15 21 1 15 22 1 16 19 1 17 19 1 18 19 1 [angles] ; ai aj ak funct c0 c1 c2 c3 2 1 3 1 117.000 292.880 ; H3-N2-H3 from ADE,CYT,GUA,GLN,ASN,ARG 4 1 2 1 121.500 292.880 ; CA-N2-H3 from ADE,CYT,GUA,ARG 4 1 3 1 121.500 292.880 ; CA-N2-H3 from ADE,CYT,GUA,ARG 5 4 6 1 121.700 585.760 ; N2-CA-N2 from ARG(OL) 5 4 1 1 119.100 585.760 ; N2-CA-N2 from ARG(OL) 6 4 1 1 119.100 585.760 ; N2-CA-N2 from ARG(OL) 7 5 4 1 121.400 418.400 ; CA-N2-CT from ARG(OL) 11 5 7 1 115.300 418.400 ; CA-N2-CT from ARG(OL) ?? 11 5 4 1 122.300 418.400 ; CA-N2-CT from ARG(OL) 19 6 4 1 122.300 418.400 ; CA-N2-CT from ARG(OL) ?? 15 6 4 1 121.400 418.400 ; CA-N2-CT from ARG(OL) ?? 19 6 15 1 115.300 418.400 ; CA-N2-CT from ARG(OL) ?? 5 7 8 1 110.600 292.880 ; HC-CT-N2 from ARG 5 7 9 1 108.400 292.880 ; HC-CT-N2 from ARG 5 7 10 1 111.700 292.880 ; HC-CT-N2 from ARG 8 7 9 1 108.100 276.144 ; HC-CT-HC 8 7 10 1 109.900 276.144 ; HC-CT-HC 9 7 10 1 108.100 276.144 ; HC-CT-HC 5 11 12 1 111.600 292.880 ; HC-CT-N2 from ARG 5 11 13 1 108.500 292.880 ; HC-CT-N2 from ARG 5 11 14 1 111.300 292.880 ; HC-CT-N2 from ARG 12 11 13 1 108.200 276.144 ; HC-CT-HC 12 11 14 1 109.600 276.144 ; HC-CT-HC 13 11 14 1 108.600 276.144 ; HC-CT-HC 6 15 16 1 110.600 292.880 ; HC-CT-N2 6 15 17 1 108.400 292.880 ; HC-CT-N2 6 15 18 1 111.700 292.880 ; HC-CT-N2 16 15 17 1 108.400 276.144 ; HC-CT-HC 16 15 18 1 109.900 276.144 ; HC-CT-HC 17 15 18 1 108.100 276.144 ; HC-CT-HC 6 19 20 1 110.600 292.880 ; HC-CT-N2 from ARG 6 19 21 1 108.500 292.880 ; HC-CT-N2 from ARG 6 19 22 1 111.300 292.880 ; HC-CT-N2 from ARG 20 19 21 1 108.200 276.144 ; HC-CT-HC 20 19 22 1 109.600 276.144 ; HC-CT-HC 21 19 22 1 108.600 276.144 ; HC-CT-HC [ dihedrals ] ; ai aj ak al funct c0 c1 c2 c3 c4 c5 2 1 4 6 3 16.31760 0.00000 -16.31760 0.00000 0.00000 0.00000 ; guanidinium ion, H3-N2-CA-N2 3 1 4 6 3 16.31760 0.00000 -16.31760 0.00000 0.00000 0.00000 ; guanidinium ion, H3-N2-CA-N2 2 1 4 5 3 16.31760 0.00000 -16.31760 0.00000 0.00000 0.00000 ; guanidinium ion, H3-N2-CA-N2 3 1 4 5 3 16.31760 0.00000 -16.31760 0.00000 0.00000 0.00000 ; guanidinium ion, H3-N2-CA-N2 7 5 4 6 3 33.20422 0.00000 -33.20422 0.00000 0.00000 0.00000 ; methylguanidinium ion,CT-N2-CA-N2 7 5 4 1 3 33.20422 0.00000 -33.20422 0.00000 0.00000 0.00000 ; methylguanidinium ion,CT-N2-CA-N2 11 5 4 6 3 33.20422 0.00000 -33.20422 0.00000 0.00000 0.00000 ; methylguanidinium ion,CT-N2-CA-N2 11 5 4 1 3 33.20422 0.00000 -33.20422 0.00000 0.00000 0.00000 ; methylguanidinium ion,CT-N2-CA-N2 4 5 7 8 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC 4 5 7 9 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC 4 5 7 10 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC 11 5 7 8 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC ?? 11 5 7 9 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC ?? 11 5 7 10 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC ?? 19 6 4 5 3 33.20422 0.00000 -33.20422 0.00000 0.00000 0.00000 ; methylguanidinium ion, CT-N2-CA-N2 19 6 4 1 3 33.20422 0.00000 -33.20422 0.00000 0.00000 0.00000 ; methylguanidinium ion, CT-N2-CA-N2 15 6 4 5 3 33.20422 0.00000 -33.20422 0.00000 0.00000 0.00000 ; methylguanidinium ion, CT-N2-CA-N2 15 6 4 1 3 33.20422 0.00000 -33.20422 0.00000 0.00000 0.00000 ; methylguanidinium ion, CT-N2-CA-N2 4 5 11 12 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC 4 5 11 13 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC 4 5 11 14 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC 7 5 11 12 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC ?? 7 5 11 13 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC ?? 7 5 11 14 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC ?? 19 6 15 16 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC ?? 19 6 15 17 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC ?? 19 6 15 18 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC ?? 4 6 15 16 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC 4 6 15 17 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC 4 6 15 18 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC 4 6 19 20 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC 4 6 19 21 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC 4 6 19 22 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC 15 6 19 20 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC ?? 15 6 19 21 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC ?? 15 6 19 22 3 0.37028 1.11086 0.00000 -1.48114 0.00000 0.00000 ; methylguanidinium ion, CA-N2-CT-HC ?? [ dihedrals ] ; ai aj ak al funct q0 cq ; 1 2 3 4 2 0.000000e+00 1.673600e+02 ; ; 4 5 6 1 2 0.000000e+00 1.673600e+02 ; ; 5 4 7 11 2 0.000000e+00 1.673600e+02 ; ; 6 4 15 19 2 0.000000e+00 1.673600e+02 ; ;;;;;; LACTIC ACID LAC [ moleculetype ] ; molname nrexcl LAC 3 [ atoms ] ; nr type resnr residue atom cgnr charge mass typeB chargeB massB 1 opls_271 1 LAC CL1 1 0.756 12.01100 ; qtot 0.756 C_3 2 opls_272 1 LAC OL2 1 -0.796 15.99940 ; qtot -0.04 O2 3 opls_272 1 LAC OL3 1 -0.806 15.99940 ; qtot -0.846 O2 4 opls_275 1 LAC CL4 1 0.557 12.01100 ; qtot -0.289 CT 5 opls_282 1 LAC HL5 1 -0.129 1.00800 ; qtot -0.418 HC 6 opls_305 1 LAC CL6 2 -0.127 12.01100 ; qtot -0.545 CT 7 opls_911 1 LAC HL7 2 -0.020 1.00800 ; qtot -0.565 HC 8 opls_911 1 LAC HL8 2 0.003 1.00800 ; qtot -0.562 HC 9 opls_911 1 LAC HL9 2 -0.005 1.00800 ; qtot -0.567 HC 10 opls_154 1 LAC OL10 3 -0.822 15.99940 ; qtot -1.389 OH 11 opls_155 1 LAC HL11 3 0.389 1.00800 ; qtot -1.000 HO [bonds] ; i j func b0 kb 1 2 1 0.12630 548940.8 ; C_3-O2,GLU,ASP 1 3 1 0.12470 548940.8 ; C_3-O2,GLU,ASP 1 4 1 0.15690 265265.6 ; C_3-CT 4 5 1 0.11000 284512.0 ; CT-HC, CHARMM 22 parameter file 4 6 1 0.15270 224262.4 ; CT-CT, CHARMM 22 parameter file 4 10 1 0.14280 267776.0 ; CT-OH 6 7 1 0.10930 284512.0 ; CT-HC, CHARMM 22 parameter file 6 8 1 0.10960 284512.0 ; CT-HC, CHARMM 22 parameter file 6 9 1 0.10950 284512.0 ; CT-HC, CHARMM 22 parameter file 10 11 1 0.09830 462750.4 ; HO-OH, SUG(OL) wlj mod 0.96-> 0.945 [ pairs ] ; ai aj funct c0 c1 c2 c3 1 7 1 1 8 1 1 9 1 1 11 1 2 5 1 2 6 1 2 10 1 3 5 1 3 6 1 3 10 1 5 7 1 5 8 1 5 9 1 5 11 1 6 11 1 7 10 1 8 10 1 9 10 1 [angles] ; ai aj ak funct c0 c1 c2 c3 4 1 2 1 114.000 585.760 ; CT-C_3-O2 GLU(OL) SCH JPC 79,2379 4 1 3 1 116.600 585.760 ; CT-C_3-O2 GLU(OL) SCH JPC 79,2379 2 1 3 1 129.400 669.440 ; O2-C_3-O2 GLU(OL) SCH JPC 79,2379 6 4 10 1 110.300 418.400 ; CT-CT-OH 1 4 6 1 112.000 527.184 ; C_3-CT-CT, AA 6 4 5 1 108.800 313.800 ; CT-CT-HC, CHARMM 22 parameter file 1 4 10 1 109.800 418.400 ; C_3-CT-OH 5 4 10 1 108.600 292.880 ; HC-CT-OH 1 4 5 1 107.300 292.880 ; C_3-CT-HC 4 6 7 1 110.200 313.800 ; CT-CT-HC, CHARMM 22 parameter file 4 6 8 1 110.800 313.800 ; CT-CT-HC, CHARMM 22 parameter file 4 6 9 1 110.200 313.800 ; CT-CT-HC, CHARMM 22 parameter file 7 6 8 1 109.100 276.144 ; HC-CT-HC 7 6 9 1 108.500 276.144 ; HC-CT-HC 8 6 9 1 107.900 276.144 ; HC-CT-HC 4 10 11 1 101.400 460.240 ; CT-OH-HO [ dihedrals ] ; ai aj ak al funct c0 c1 c2 c3 c4 c5 5 4 1 2 3 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 ; HC-CT-C_3-O2, carboxylate ion 5 4 1 3 3 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 ; HC-CT-C_3-O2, carboxylate ion 6 4 1 2 3 3.43088 0.00000 -3.43088 0.00000 0.00000 0.00000 ; CT-CT-C_3-O2, carboxylate ion 6 4 1 3 3 3.43088 0.00000 -3.43088 0.00000 0.00000 0.00000 ; CT-CT-C_3-O2, carboxylate ion 10 4 1 2 3 2.28446 0.00000 -2.28446 0.00000 0.00000 0.00000 ; COOH terminus 10 4 1 3 3 2.28446 0.00000 -2.28446 0.00000 0.00000 0.00000 ; COOH terminus 1 4 6 7 3 -0.47070 -1.41210 0.00000 1.88280 0.00000 0.00000 ; C_3-CT-CT-HC, carboxylate ion 5 4 6 7 3 0.62760 1.88280 0.00000 -2.51040 0.00000 0.00000 ; HC-CT-CT-HC, hydrocarbon *new* 11/99 1 4 6 9 3 -0.47070 -1.41210 0.00000 1.88280 0.00000 0.00000 ; C_3-CT-CT-HC, carboxylate ion 5 4 6 9 3 0.62760 1.88280 0.00000 -2.51040 0.00000 0.00000 ; HC-CT-CT-HC, hydrocarbon *new* 11/99 1 4 10 11 3 -0.44350 3.83255 0.72801 -4.11705 0.00000 0.00000 ; CT-CT-OH-HO, alcohols AA 5 4 10 11 3 0.94140 2.82420 0.00000 -3.76560 0.00000 0.00000 ; HC-CT-OH-HO, alcohols AA 6 4 10 11 3 -0.44350 3.83255 0.72801 -4.11705 0.00000 0.00000 ; CT-CT-OH-HO, alcohols AA 7 6 4 10 3 0.97905 2.93716 0.00000 -3.91622 0.00000 0.00000 ; HC-CT-CT-OH, alcohols, ethers AA 8 6 4 10 3 0.97905 2.93716 0.00000 -3.91622 0.00000 0.00000 ; HC-CT-CT-OH, alcohols, ethers AA 9 6 4 10 3 0.97905 2.93716 0.00000 -3.91622 0.00000 0.00000 ; HC-CT-CT-OH, alcohols, ethers AA [ dihedrals ] ; ai aj ak al funct q0 cq ; 1 2 3 4 2 0.000000e+00 1.673600e+02 ; ;;;;;; LACTIC ACID CO2 [ moleculetype ] ; molname nrexcl CO2 2 [ atoms ] ; nr type resnr residue atom cgnr charge mass typeB chargeB massB 1 JACS_000 1 CO2 CC1 1 0.700 12.01100 ; 2 JACS_001 1 CO2 OC2 1 -0.350 15.99940 ; 3 JACS_001 1 CO2 OC3 1 -0.350 15.99940 ; [bonds] ; i j func b0 kb 1 2 1 0.11600 432600.0 ; C_3-O2,GLU,ASP 1 3 1 0.11600 432600.0 ; C_3-O2,GLU,ASP [ pairs ] ; ai aj funct c0 c1 c2 c3 [angles] ; ai aj ak funct c0 c1 c2 c3 2 1 3 1 180.000 235.2 ; CT-C_3-O2 GLU(OL) SCH JPC 79,2379 [ dihedrals ] ; ai aj ak al funct c0 c1 c2 c3 c4 c5 [ dihedrals ] ; ai aj ak al funct q0 cq 年发在JPCB上面的一篇文章的力场，可以仿照这个建立适此力场是王勇(浙江大学)在2007 合自己的离子液体的力场。在建立后，最好在TOP文件里面INCLUDE的时候，应该把ITP文件里 面的DEFAULTS那几行删除，这两个只能留一个。对于BMIMBF系列的离子液体，就得参考汪4 文川2004年的一篇文献，自己仿照上面的格式，自己建立力场，建立力场的过程需要参考 MANUAL。 也需要自己建立GRO文件，在GROMACS里面首先使用EDITCONF扩大盒子的半径，然后使 用GENCONF建立超盒子，使用GENBOX加水，进行分子动力学模拟。 对于建立好的盒子，如果想插入CH和CO什么的，就得编程自己插入，源程序是这样的: 42 integer nsurf,natom0,nco2,tatom,nch4 real xc,yc,zc,dl,blx,bly,blz * blx,bly,blz are the length of simulation box in x,y,z direction, respectively. parameter(nsurf=3456,tatom=1,natom0=3456,nch4=0,nco2=100, & blx=4.04680,bly=4.04590,blz=4.02660) character resid(natom0)*8,atomname(natom0)*3,co2name(3)*3 integer atomid(natom0),atomsn,molsn real coor(natom0,6),co2(nco2*3,3),ch4(nch4,3) open(10,file='sorbent.gro',status='old') do 20 i=1,natom0 read(10,*)resid(i),atomname(i),atomid(i),coor(i,1),coor(i,2), & coor(i,3) 20 continue close(10) do 25 i=1,nco2*3 do 26 j=1,3 co2(i,j)=0.0 26 continue 25 continue *************** add CO2 ************************* do 30 i=1,nco2 50 xc=RAN2(IDUM)*blx yc=RAN2(IDUM)*bly zc=RAN2(IDUM)*blz do 40 j=1,natom0 dl=sqrt((xc-coor(j,1))**2+(yc-coor(j,2))**2+(zc-coor(j,3))**2) if(dl.lt.0.25) goto 50 40 continue if(i.ge.2)then do 60 k=1,i-1 kp=(k-1)*3+2 dl=sqrt((xc-co2(kp,1))**2+(yc-co2(kp,2))**2+(zc-co2(kp,3))**2) if(dl.lt.0.25) goto 50 60 continue endif *************** add CO2 ************************* co2((i-1)*3+2,1)=xc co2((i-1)*3+2,2)=yc co2((i-1)*3+2,3)=zc co2((i-1)*3+1,1)=xc-0.115 co2((i-1)*3+1,2)=yc co2((i-1)*3+1,3)=zc co2((i-1)*3+3,1)=xc+0.115 co2((i-1)*3+3,2)=yc co2((i-1)*3+3,3)=zc 30 continue *************** add CH4 ************************* do 62 i=1,nch4 63 xc=RAN2(IDUM)*blx yc=RAN2(IDUM)*bly zc=RAN2(IDUM)*blz do 64 j=1,natom0 dl=sqrt((xc-coor(j,1))**2+(yc-coor(j,2))**2+(zc-coor(j,3))**2) if(dl.lt.0.25) goto 63 64 continue do 65 kp=1,nco2*3 dl=sqrt((xc-co2(kp,1))**2+(yc-co2(kp,2))**2+(zc-co2(kp,3))**2) if(dl.lt.0.25) goto 63 65 continue ch4(i,1)=xc ch4(i,2)=yc ch4(i,3)=zc 62 continue ****************** add CH4 ************************ ******************OUTPUT THE NEW GRO FILE ************************ open(70,file='system.gro',access='append') write(70,'(A6)')'system' ktotal=natom0+3*nco2+nch4 write(70,'(I8)')ktotal do 120 i=1,nch4 molsn=i atomsn=i write(70,666)molsn,'MET','CME',atomsn, & ch4(atomsn,1), ch4(atomsn,2),ch4(atomsn,3) 120 continue co2name(1)='OAA' co2name(2)='CAC' co2name(3)='OAB' do 100 i=1,nco2 do 110 j=1,3 molsn=i+nch4 atomsn=(i-1)*3+j+nch4 kjishu=(i-1)*3+j write(70,666)molsn,'DRG',co2name(j),atomsn, & co2(kjishu,1),co2(kjishu,2),co2(kjishu,3) 110 continue 100 continue do 80 i=1,nsurf do 90 j=1,tatom molsn=nco2+i+nch4 atomsn=nco2*3+(i-1)*tatom+j+nch4 kjishu=(i-1)*tatom+j write(70,666)molsn,'MOR',atomname(kjishu),atomsn,coor(kjishu,1), & coor(kjishu,2),coor(kjishu,3) 90 continue 80 continue ******************OUTPUT THE NEW GRO FILE ************************ 666 format(I5,A3,4X,A3,I5,3F8.3) write(70,'(3F10.5)')blx,bly,blz write(70,*)'' close(70) end FUNCTION RAN2(IDUM) 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 iy=iv(j)-idum2 iv(j)=idum if(iy.lt.1)iy=iy+IMM1 ran2=min(AM*iy,RNMX) return END 其中输入输出文件是这样 Sorbent.gro: 1SLC SIZ 1 0.859 0.105 0.875 2SLC SIZ 2 0.607 0.049 1.059 3SLC SIZ 3 0.565 0.131 0.013 4SLC SIZ 4 0.250 0.133 0.025 5SLC SIZ 5 0.150 0.052 1.091 6SLC SIZ 6 0.359 0.104 0.864 7SLC SIZ 7 0.860 1.676 0.901 8SLC SIZ 8 0.617 1.754 1.081 9SLC SIZ 9 0.565 1.670 0.039 10SLC SIZ 10 0.252 1.670 0.043 11SLC SIZ 11 0.154 1.761 1.101 12SLC SIZ 12 0.373 1.673 0.891 13SLC SIZ 13 0.153 1.918 0.204 14SLC SIZ 14 0.404 1.974 0.388 **************************** system 3756 1DRG OAA 1 3.477 1.506 1.399 1DRG CAC 2 3.592 1.506 1.399 1DRG OAB 3 3.707 1.506 1.399 2DRG OAA 4 3.782 1.647 0.485 2DRG CAC 5 3.897 1.647 0.485 2DRG OAB 6 4.012 1.647 0.485 3DRG OAA 7 3.082 1.639 2.761 3DRG CAC 8 3.197 1.639 2.761 3DRG OAB 9 3.312 1.639 2.761 4DRG OAA 10 1.775 3.039 2.157 4DRG CAC 11 1.890 3.039 2.157 4DRG OAB 12 2.005 3.039 2.157 5DRG OAA 13 0.973 2.895 1.369 5DRG CAC 14 1.088 2.895 1.369 5DRG OAB 15 1.203 2.895 1.369 6DRG OAA 16 0.826 3.432 2.987 6DRG CAC 17 0.941 3.432 2.987 6DRG OAB 18 1.056 3.432 2.987 ******************************** 此段程序，是为了遍历到所有的CH分子和CO中C和O原子而已，MOLSN是所有的分子个数，42 而ATOMSN是所有的原子个数，KJISHU是所有的CO中原子个数，其中(I-1)*3+J是为了遍历2 到所有的原子而故意这么设置的。 离子液体的MD模拟(DLPOLY) DLPOLY安装，很多帖子已经说明了，这里主要是讨论怎么生成DLPOLY需要的CONFIG和FIELD文件，其他的控制文件按照规定比较好弄。CONFIG文件，可以使用其他的软件生成类似于PDB格式这样的文件，MS就可以生成离子液体分子动力学盒子，但是MS不能直接生成含有F元素的离子液体的盒子，比如说生成BMIMBF4的盒子，可以先生成BMIMPH4的盒子， 然后选择盒子中所有的P，换成B,然后选择ATOM SELECTION中的CONNETED,就选择了所有的H,换成F就成。然后使用程序来自动生成CONFIG文件，MS可以保存为PDB格式的文件，可以通过EXCEL分列，去除其中杂乱的因素，只留下数据，然后操作。MS的pdb格式结构文件转化为DLPOLY的CONFIG的FORTRAN程序 program dlpolyinput implicit none integer i, numcoordinate, pbc, numatom real, PARAMETER :: K=0 character *20 atomname, systemname real x(200),y(200),z(200), lbox open(20,file='CONFIG',status='unknown') !打开一个输出文件 open(1,file='input.txt',status='unknown')!打开文件，并读取 read(1,*) systemname write(20,"(A80)") systemname read(1,*) numcoordinate, pbc, numatom write(20,"(I10,I10,I10)") numcoordinate, pbc, numatom read(1,*) lbox write(20,"(F20.3,F20.3,F20.3)") lbox, K, K write(20,"(F20.3,F20.3,F20.3)") K, lbox, K write(20,"(F20.3,F20.3,F20.3)") K, K, lbox do i=1, 200 read(1,*)atomname,x(i),y(i),z(i) write(20,"(A8,I10)") atomname, i write(20,"(F20.3,F20.3,F20.3)")x(i), y(i) ,z(i) end do end 这是input格式.. input-turn-config 0 1 200 28.850 C1 -7.554 -6.98 -10.036 C2 -8.917 -7.191 -10.626 C2 -8.954 -8.558 -11.308 C2 -10.315 -8.796 -11.93 S14 -7.47 -5.332 -9.259 C1 1.742 -9.347 11.89 C2 2.961 -8.475 11.609 C2 3.55 -7.894 12.896 C2 4.758 -7.009 12.613 S14 1.068 -10.075 10.362 C1 0.999 -2.023 -3.01 C2 -0.235 -2.861 -3.327 C2 -0.654 -2.719 -4.793 C2 -1.898 -3.543 -5.105 S14 1.527 -2.224 -1.278 C1 12.82 8.949 -13.352 C2 12.856 10.457 -13.568 C2 11.448 11.052 -13.642 C2 11.484 12.558 -13.876 S14 14.495 8.247 -13.21 但是很多离子液体结构使用MS不能生成立方盒子结构，这个时候，可以使用GROMACS里面的EDITCONF命令建立合适大小的盒子。得到的最终数据文件，可以通过EXCEL分列得到需要的XYZ数据，然后也是通过这个程序运行生成DLPOLY需要的CONFIG文件。GROMACS需要的GRO文件，需要手动编写，因为离子液体的原子数不算太多。而对于FIELD文件，确定好了一个离子液体的FIELD文件以后，只要在分子数上添加一个随意的数字，就可以得到含有这个数字的盒子的FIELD文件。GUI的编译在DLPOLY编译好了以后，已经编译成功了，此时不必使用javac *.java了，直接在java子文件夹下，输入启动GUI命令就可以得到GUI. MS生成的文件转化为VASP的格式