模拟心得
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