Gaussian中非内置的理论方法和泛函的用法

Gaussian中非内置的理论方法和泛函的用法

文/Sobereva  2016-Sep-4


本文把Gaussian中没有内置关键词,但能利用IOp或其它手段等效实现的一些理论方法的用法进行介绍,本文对应于Gaussian 09 D.01/E.01。估计其中一些方法到了下个Gaussian大版本肯定会有直接对应的关键词。

1 微扰相关的方法

2.1 SCS-MP2

这是Grimme在J. Chem. Phys., 118, 9095 (2003)提出的对MP2的改进。MP2相关能可以写为电子对的加和,SCS-MP2把自旋相同电子对的贡献乘上1/3以削弱之,而把自旋相反的贡献乘上6/5以增强之。SCS-MP2算反应能比MP2有了很大提高,达到QCISD级别甚至个别时达到QCISD(T)级别。几何结构也比MP2略有改进。但势垒计算精度比MP2反倒略低一点,弱相互作用计算精度恶化不少。其它方面和MP2差不多。整体性能明显不如目前最好的双杂化泛函,所以如今用处也不太大。

SCS-MP2的使用方法是在MP2计算时额外加上IOp(3/125=0333312000)。3/125=MMMMMNNNNN代表将自旋平行项与自旋反平行项的贡献以MMMMM/10000比NNNNN/10000的比例组合。

MP2计算时从输出中看到
     alpha-alpha T2 =       0.1336200196D-01 E2=     -0.4731575741D-01
     alpha-beta  T2 =       0.8298753242D-01 E2=     -0.3004708246D+00
     beta-beta   T2 =       0.1336200196D-01 E2=     -0.4731575741D-01
SCS-MP2计算时为
     alpha-alpha T2 =       0.1336200196D-01 E2=     -0.1577034194D-01
     alpha-beta  T2 =       0.8298753242D-01 E2=     -0.3605649895D+00
     beta-beta   T2 =       0.1336200196D-01 E2=     -0.1577034194D-01
可见alpha-alpha、beta-beta的能量贡献都被乘上了1/3,而alpha-beta乘上了6/5。

2.2 SOS-MP2

于J. Chem. Phys., 121, 9793(2004)中提出,忽略了自旋平行贡献,反平行部分乘以1.3,结果比MP2烂,特别是弱相互作用。不过号称通过Laplace变换方法可以令标度降到O(N^4)从而使花费比MP2的O(N^5)低。不过再便宜也没普通DFT便宜,MP2精度尚且不及选得最合适的DFT泛函,所以SOS-MP2这个破玩意儿并没什么卵用。大家知道有这么个垃圾就够了。

用法是MP2 IOp(3/125=0000013000)。由于Gaussian并未对之进行优化,所以耗时不会比MP2低。

2.3 SCSN-MP2

在J. Chem. Theory Comput., 3, 80 (2007)中提出,忽略自旋反平行部分,自旋平行部分的系数为1.76,拟合自核酸碱基对相互作用能。弱相互作用计算精度挺不错,误差比MP2低两三倍,有使用价值。另外还有Mol. Phys., 105, 1073 (2007)中提出的SCS(MI)-MP2,MI代表molecular interaction,它与SCSN-MP2的提出目的和精度都很类似,并对不同基组拟合了不同的系数,这里就不提了。从SCSN-MP2和SCS(MI)-MP2的系数都能看出,对于计算弱相互作用的目的,应该削减MP2自旋反平行而增大自旋平行的贡献,这和一般目的的SCS-MP2恰好相反。

用法是MP2 IOp(3/125=1760000000)。

2.4 MP2.5

是在ChemPhysChem, 10, 282-289 (2009)中提出的,将MP2总能量加上乘以0.5的MP3三阶微扰能。计算量和MP3一样。系数0.5来自于分析计算精度、基组依赖性和理论意义。各方面性能比MP2强,特别是计算弱相互作用准确度很高,号称在中等基组下(不加弥散亦可)就能接近CCSD(T)/CBS。

在Gaussian中用MP3计算完后,会看到比如
 E2 =    -0.3951023394D+00 EUMP2 =    -0.11430702788798D+03
 Keep R2 and R3 ints in memory in symmetry-blocked form, NReq=9641770.
 DD1Dir will call FoFMem   1 times, MxPair=        42
 NAB=    21 NAA=     0 NBB=     0.
 E3=       -0.48480932D-02        EUMP3=      -0.11431187598D+03
MP2.5能量就是-114.30702788798-0.0048480932*0.5=-114.3094519

后来在Phys. Chem. Chem. Phys., 13, 21121 (2011)中还提出了MP2.X,X是给三阶微扰能乘上的系数。MP2.5用在小基组上结果不如在大基组好,为了使得较小基组下就能得到不错的弱相互作用计算结果,MP2.X对从小到大的基组都通过S66测试集拟合了MP3校正能的权重系数,这使得不同基组下(乃至低至6-31G*)得到的弱相互作用能精度都相仿佛,和MP2.5/aug-cc-pVTZ下差不多。虽看似很好,但对更多的体系的可靠性还有待广泛验证。

由于MP2.5、MP2.X能量只能在MP3单点任务结束后根据输出信息手动计算,没法像SCS-MP2那样直接用IOp,所以这两种方法只能算能量,没法用来优化、计算频率等,除非自己去改Gaussian代码实现。

2.5 SCS-MP3

Grimme在J. Comput. Chem., 24, 1529 (2003)中提出的,是SCS-MP2能量加上乘以0.25的三阶微扰校正能。热化学性能比SCS-MP2好很多,往往接近QCISD(T)。但弱相互作用能精度不算理想。

用法是MP3 IOp(3/125=0333312000),然后把输出文件中EUMP2能量加上乘以了0.25的E3。


顺带一提,SCS的思路也被用在CI和耦合簇方法上。SCS-CCSD的正反自旋参数来自于拟合几十个反应能。算反应能,特别是算弱相互作用能比CCSD好很多。Phys. Chem. Chem. Phys., 12, 9611 (2010)中提出的SCS-MI-CCSD则令参数向S22弱相互作用测试集拟合,算弱相互作用极为接近金标准CCSD(T),比MP2.5、SCS-CCSD又明显更好。但可惜Gaussian09中没法实现这些SCS的耦合簇方法,因为程序也不把耦合簇的自旋平行和反平行的贡献独立输出。


2 非内置的普通DFT泛函的用法

2.1 规则

在Gaussian中可以通过IOp在一定程度上自定义泛函。简单来说,若定义IOp(3/76=ab)、IOp(3/77=cd)、IOp(3/78=ef),这里abcdef都已经除以了10000,则所用泛函的实际形式是(下文说的GGA也包含meta-GGA)
E_XC = a*[ d*E_X_LDA + c*ΔE_X_GGA ] + b*E_X_HF + f*E_C_LDA + e*ΔE_C_GGA
每一项的含义如下
E_X_LDA:LDA交换项
ΔE_X_GGA:GGA交换项对LDA交换项的梯度修正部分
E_X_HF:HF交换项。前头的系数b就是常说的HF杂化成份
E_C_LDA:LDA相关项
ΔE_C_GGA:GGA相关项对LDA相关项的梯度修正部分

无论对于交换还是相关部分,都是E_GGA=ΔE_GGA+E_LDA。比如E_X_B88=E_X_LDA+ΔE_X_B88,E_C_LYP=E_C_LDA+ΔE_C_LYP。

对于DFT,a总为1(其实这个参数纯属多余,即便有用也可以直接融合到d、c里嘛)。对于纯泛函,满足d=c=1、f=e=1、b=0;对于单参数杂化泛函,比如PBE0、BHandHLYP,满足d=c、d+b=1、f=e。对于B3LYP这样的三参数杂化泛函才可能d!=c、f!=e。而对于HF,显然b=1,其它皆为0。

用#P在计算时会有诸如这样的输出:
IExCor= 402 DFT=T Ex=B+HF Corr=LYP ExCW=0 ScaHFX= 0.200000
ScaDFX= 0.800000 0.720000 1.000000 0.810000
其中ScaHFX就是b,ScaDFX后面的值是d c f e。Ex就是所用的交换泛函,当前是B88和HF杂化。Corr是所用的相关泛函,当前是LYP。

2.1 B3LYP及变体

老不死的B3LYP定义为
E_XC_B3LYP = (1-a0)*E_X_LDA + a0*E_X_HF + aX*ΔE_X_B88 + E_C_VWN + aC*ΔE_C_LYP
其中a0=0.2、aX=0.72、aC=0.81。

显然,换算成上一节的系数表达方式,就是a=1.0,b=0.2,c=0.72,d=0.8,e=0.81,f=1.0。相应地,计算的时候关键词就写成BLYP IOp(3/76=1000002000,3/77=0720008000,3/78=0810010000)。

B3LYP在提出的时候没有明确说清楚LYP用的VWN到底是VWN3还是VWN5,二者结果有一点差异。Gaussian里默认用的是前者,比后者更好一点点,见Chem. Phys. Lett., 268, 345 (1997)的对比。但有的程序默认用的则是VWN5,如GAMESS-US。如果想用VWN5的B3LYP的话,就可以写BV5LYP IOp(3/76=1000002000,3/77=0720008000,3/78=0810010000),这里V5LYP相关泛函是指把VWN5作为LYP里的LDA部分。

在J. Chem. Phys., 117, 4729 (2002)中提出了B3LYP*,它把B3LYP的HF成份从20%降低到15%,这大大解决了B3LYP高估过渡金属配合物高自旋态稳定性的问题(即低估了高、低自旋态间的能量差)。用B3LYP*就写BLYP IOp(3/76=1000001500,3/77=0720008500,3/78=0810010000),即把X_HF的系数b降到0.15,相应地把X_LDA的系数d提升到0.85,使二者加和仍保持为1。

HF成份对结果有很多系统性的影响,比如HF成份越大TDDFT算的激发能越高。因此可以容易地调节b并相应地修改d,使得算出的光谱和实际尽量接近(虽然洒家鄙夷这种刻意往实验结果上凑的勾当)。

2.2 PBE0及变体

PBE0在Gaussian里写作PBE1PBE,是把PBE泛函中引入25% HF成份,即a=1.0,b=0.25,c=d=0.75,e=f=1.0。相当于PBEPBE IOp(3/76=1000002500,3/77=0750007500)。由于IOp(3/78=1000010000)是默认的,所以可以不用写。

J. Chem. Phys., 138, 021104 (2013)定义了一个无聊的泛函PBE0-1/3,就是把HF成份增加到1/3,故对应于PBEPBE IOp(3/76=1000003333) IOp(3/77=0666706667)。

DFT-D3原文里顺带定义了PBE38,是把HF成份改为3/8=37.5%,测试表明算含频极化率是最好的,对应于PBEPBE IOp(3/76=1000003750) IOp(3/77=0625006250)。

2.3 明尼苏达系列老泛函

Truhlar组之前搞过一大堆乱七八糟的泛函,什么MPW1K、PBE1KCIS,也没多少人用,就他们自己玩得happy。从M05起我看才算真正步入正轨。那些烂七八糟的老明尼苏阿达系列泛函我强烈不推荐使用,早被更好的泛函完爆了,初学者也甭老效仿文献试图去用。如果非要用,去看他们网页http://comp.chem.umn.edu/info/DFT.htm,能通过Gaussian实现的泛函都给了用法。而像PW6B95这样的在Gaussian中即便通过IOp也使不了的,死活非要用就必须找他们要修改版的Gaussian了。

2.4 对长程校正泛函调节ω

Gaussian中对纯泛函可以用LC-前缀做长程校正,比如LC-BLYP。此方法将1/r12算符划分为短程和长程部分,短程不用HF交换项,而长程用100% HF交换项。此做法有一个参数ω很重要,此值越大长程部分涵盖范围就越广,会明显影响结果。LC-加到各种纯泛函上默认是ω=0.47,ωB97X默认是ω=0.3,LC-wPBE默认是0.4。如果要调节这个参数,就把3/107和3/108都设为MMMMM00000,代表ω值为MMMMM/10000。比如要把LC-wPBE的ω设0.3就写LC-wPBE IOp(3/107=0300000000,3/108=0300000000)。



3 非内置的双杂化泛函的用法

2.1 老式双杂化泛函

比较老的双杂化泛函,如B2PLYP、mPW2PLYP、B2GP-PLYP没有考虑色散校正也没有对掺入的二阶微扰能E2以SCS方式考虑。这类老双杂化泛函的交换相关能通式为
E_XC = (1-c1)*E_X_GGA + c1*E_X_HF + (1-c2)*E_C_GGA + c2*E2

最老、最知名,但也是最差的双杂化泛函B2PLYP是J. Chem. Phys., 124, 034108 (2006)中提出的,也是Gaussian09中内置的,定义为
E_XC=0.47*E_X_B88+0.53*E_X_HF+0.73*E_C_LYP+0.27*E2
计算时会看到相应的泛函定义信息,只要看懂了2.1节就自然能理解:
 IExCor=  419 DFT=T Ex+Corr=B2PLYP ExCW=0 ScaHFX=  0.530000
 ScaDFX=  0.470000  0.470000  0.730000  0.730000 ScalE2=  0.270000  0.270000
这里比2.1节看到的多了ScalE2,后面两个值是给E2的自旋平行和反平行部分分别乘的系数,这里都是0.27,所以说没有以SCS方式考虑E2。

B2GP-PLYP是在B2PLYP之后提出的,原文见J. Phys. Chem. A, 112, 12868 (2008),精度比B2PLYP好不少,但比起后来以SCS方式考虑E2的双杂化泛函还是差多了,其系数c1=0.65,c2=0.36,纯泛函部分还是BLYP,可明确写为
E_XC=0.35*E_X_B88+0.65*E_X_HF+0.64*E_C_LYP+0.36*E2

Gaussian09支持的很可惜只有最老最烂的B2PLYP和后来提出的单精度与之半斤八两的mPW2PLYP。B2GP-PLYP在Gaussian09中没有内置,好在若想用的话可以写
B2PLYP/cc-pVTZ IOp(3/125=0360003600,3/76=1000006500,3/77=0350003500,3/78=0640006400)
这里还写B2PLYP是为了让Gaussian做双杂化计算,但是参数已经被我们改了。3/125=MMMMMNNNNN代表给E2的自旋平行和反平行部分分别乘上MMMMM/10000和NNNNN/10000,其它参数前面已经讲过了。

注:也可以把3/76=1000006500,3/77=0350003500改为3/76=0350006500,3/77=1000010000,根据2.1节的表达式就知道这两种写法是等价的。

2.2 较新的双杂化泛函

J. Phys. Chem. C, 114, 20801 (2010)中提出了DSD-BLYP,它考虑了色散校正,并且以SCS方式考虑了E2的贡献,精度比B2GP-PLYP有一定提升。从此之后各种通式类似的双杂化泛函冒出来不少,比如DSD-PBEP86、DSD-PBEhB95等,此二者都比DSD-BLYP精度又有进一步提升。这些泛函通式写为这样:
E_XC = (1-cX)*E_X_GGA + cX*E_X_HF + cC*E_C_GGA + cO*E2_OS + cS*E2_SS + E_D
这里E2_OS和E2_SS分别是E2中自旋相反和相同部分的贡献,E_D是DFT-D3色散校正能,一般都用BJ阻尼。

J. Comput. Chem., 34, 2327 (2013)是一篇好文章,在其中把各种常见纯泛函组合代入到上面的公式来拟合了系数,包括色散校正参数。经过比较,发现DSD-PBEP86-D3(BJ)精度名列前茅,而且最难得的是能直接在许多主流量化程序里用,包括Gaussian09。在此文的补充材料里可以找到各种纯泛函对应的双杂化泛函的参数,在“DSD-DFT-D3BJ”表格中可见对于DSD-PBEP86-D3(BJ),cX=0.69,cC=0.44,cO=0.52,cS=0.22,DFT-D3(BJ)的参数为s6=0.48,a2=5.6。

下面说怎么在G09中用DSD-PBEP86-D3(BJ)。DFT-D3(BJ)的参数在Gaussian09输入文件里没法直接设,自定义的话必须得设定环境变量。对于Linux版,bash环境,在计算之前先把以下内容输入到控制台中
export GAUSS_DFTD3_S6=480000
export GAUSS_DFTD3_SR6=0
export GAUSS_DFTD3_S8=0
export GAUSS_DFTD3_ABJ1=0
export GAUSS_DFTD3_ABJ2=5600000

之后输入文件里的关键词就写
B2PLYP/cc-pVTZ IOp(3/125=0220005200,3/76=1000006900,3/77=0310003100,3/78=0440004400,3/74=1004) empiricaldispersion=GD3BJ
这里3/74=1004代表把X_GGA设为PBE,C_GGA设为P86。

注1:也可以像此泛函提出者的网页http://www.compchem.me/dsd-pbep86-functional中那样写为IOp(3/125=0220005200,3/78=0440004400,3/76=0310006900,3/74=1004),是等价的。

注2:虽然JCC文中测试表明DSD-PBEhB95-D3(BJ)综合性能最好,但在G09中通过上述方法使用会报错,需要对代码进行特殊修改。

类似地,还可以根据补充材料里的信息使用各种其它的双杂化泛函,但大多不如DSD-PBEP86-D3(BJ)。另外,如果不想用色散校正的版本(比如只有G09 D.01以前版本或者嫌设定环境变量麻烦),即只用DSD-PBEP86,参数是为cX=0.72,cC=0.44,cO=0.51,cS=0.36,故应当写
B2PLYP/cc-pVTZ IOp(3/125=0360005100,3/76=1000007200,3/77=0280002800,3/78=0440004400,3/74=1004)
如文中的数据所示,不用DFT-D3(BJ)的话整体误差会有所增加,但个别项目可能反倒误差会降低(不排除巧合可能)。

添加新评论