Gaussian中非内置的理论方法和泛函的用法
Gaussian中非内置的理论方法和泛函的用法
Usage of nonbuilt-in theoretical methods and functionals in Gaussian
文/Sobereva@北京科音
First release: 2016-Sep-4 Last update: 2022-Feb-14
本文把Gaussian中没有内置关键词，但能利用IOp或其它手段等效实现的一些理论方法的用法进行介绍。
首先要提醒一点，对于诸如opt freq这样多步任务，输入文件里设的IOp只对第一步生效，不会自动传递给后面的任务。所以通过IOp使用非内置方法的时候不要为省事把opt和freq放在一起用。
1 微扰相关的方法
1.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。
1.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低。
1.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)。
1.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代码实现。
1.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.2 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.3 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.4 明尼苏达系列老泛函
Truhlar组之前搞过一大堆乱七八糟的泛函，什么MPW1K、PBE1KCIS，也没多少人用，就他们自己玩得happy。从M05起我看才算真正步入正轨。那些烂七八糟的老明尼苏阿达系列泛函我强烈不推荐使用，早被更好的泛函完爆了，初学者也甭老效仿文献试图去用。如果非要用，去看他们网页http://comp.chem.umn.edu/info/DFT.htm，能通过Gaussian实现的泛函都给了用法。而在Gaussian中即便通过IOp也使不了的，死活非要用就必须找他们要修改版的Gaussian了。2.5 对长程校正泛函调节ω
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)。
对于范围分离泛函还可以通过调节α、β参数来让短程极限和长程极限的HF成份对应于特定的数值，笔者专门在这篇博文里进行了介绍：《在Gaussian中自定义范围分离泛函的方法》（http://sobereva.com/550）。
2.6 QTP17泛函
在《正确地认识分子的能隙(gap)、HOMO和LUMO》（http://sobereva.com/543）中笔者介绍了QTP系列泛函。其中QTP17泛函结合aug-cc-pVTZ基组计算出的价层和内层轨道能量都能较好满足Koopmans定理，因此轨道能量的负值直接就可以近似当做光电子谱测定的外层、深层、内核电离能了。而且还可以直接结合Multiwfn绘制光电子谱而不需要对轨道能量做任何校正，参见《使用Multiwfn绘制光电子谱》（http://sobereva.com/478）。
在Gaussian中QTP17的用法是写# SLYP/aug-cc-pVTZ IOp(3/76=1000006200,3/77=0000003800,3/78=0800010000)。注意由于QTP17泛函参数是在aug-cc-pVTZ基组下训练的，若改成其它基组精度会下降。
3 非内置的双杂化泛函的用法
3.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节的表达式就知道这两种写法是等价的。
3.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)精度名列前茅，而且难得的是它能直接在许多主流量化程序里不需要修改源代码就能使用。在此文的补充材料里可以找到各种纯泛函对应的双杂化泛函的参数，在“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。
DSD-PBEP86-D3(BJ)从Gaussian 16开始已经内置了，写DSDPBEP86关键词即可。这里说一下怎么在Gaussian 09中通过特殊方式等价地使用。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)的话整体误差会有所增加，但个别项目可能反倒误差会降低（不排除巧合可能）。
2020-Sep-19补充：在J. Phys. Chem. A, 123, 5129 (2019)中DSD-PBEP86-D3(BJ)的作者又提出了其修改版revDSD-PBEP86-D3(BJ)，精度有了明显提升，在所有双杂化泛函中精度已经是名列前茅了，精度和目前（2020年）最佳的双杂化泛函ωB97M(2)已经很接近了。在Gaussian 16中可以通过以下关键词使用
# DSDPBEP86/基组 IOp(3/125=0079905785,3/78=0429604296,3/76=0310006900,3/74=1004) em=gd3bj IOp(3/174=0437700,3/175=-1,3/176=0,3/177=-1,3/178=5500000)
3.3 其它双杂化泛函
18碳环是非常独特的体系，笔者做过大量研究，汇总见http://sobereva.com/carbon_ring.html。在Phys. Chem. Chem. Phys. (2022) DOI: 10.1039/d1cp04996h中作者表明范围分离双杂化泛函RSX-QIDH（也叫RSX-PBE-QIDH）可以比较好地计算碳环键长均等的过渡态结构和非均等的极小点结构的能量差。此泛函在Gaussian 16中的用法为：
PBEQIDH IOP(3/74=1001009,3/78=0666606666,3/107=0270000000,
3/108=0270000000,3/119=0310000000,3/120=0310000000,
3/125=0333403334,3/130=06900,3/131=06900)