谈谈能量的基组外推

谈谈能量的基组外推

文/Sobereva @北京科音
First release: 2012-Dec-30   Last update: 2017-Jun-8

 

1 前言

量子化学计算中一个主要的误差来源是基组不够大,或者说距离完备基组极限差距较大。用越大的基组计算耗时越多。对于一些性质,如能量、优化出的结构、偶极矩等,随着基组的增大计算结果会逐渐收敛。只要知道收敛趋势的数学形式,就可以通过较小的基组的计算结果外推出较大基组下的结果,乃至完备基组(Complete basis set, CBS)下的结果。对于很高精度计算,例如小体系高精度弱相互作用计算,外推到CBS已经成了司空见惯的做法。本文将介绍单点能外推至CBS的方法,并给出实例。
 

2 HF/DFT能量的外推

基组外推的前提条件是所用的基组序列必须是以系统方式构建的。对于这样的基组序列,随着基组的增大,可以确保结果的误差逐步、平稳地降低,有明确的收敛趋势。对于Dunning相关一致性基组(cc-pVnZ系列、aug-cc-pVnZ系列等)、ANO系列基组、Jensen的极化一致性基组(pc-n系列)、Ahlrich的def2系列基组等,都符合这个要求。比如可以用cc-pVDZ、cc-pVTZ、cc-pVQZ计算的三个能量外推到CBS的能量。而诸如Pople系列基组则不符合要求,例如不适合使用诸如3-21G*、6-31G*、6-311+G(2d,p)的计算结果外推至CBS,虽说他们的误差肯定越来越小,但误差降低的趋势并不平滑,没有可循的收敛趋势以用于外推。

基组外推的公式并不唯一,有的简单有的复杂。对于不同基组序列,最佳的外推公式也可能有微小的不同。

对于Dunning相关一致性基组(或ANO系列),SCF能量(指HF、DFT能量)收敛趋势符合以下指数关系,可以通过它们来外推至CBS:
E_SCF(L)=E_SCF(∞)+A*exp(-α*L)

E_SCF(L)=E_SCF(∞)+A*exp(-α*√L)
这两种用谁都可以,结果差异甚小。其中A、α是需要确定的参数,E_SCF(∞)就是我们最终想要得到的CBS下的SCF能量。L是基组所含基函数所具有的最高角量子数。例如对于Dunning相关一致性基组,L=2对应cc-pVDZ,因为它对于主族重原子最高角动量为d,角量子数是2。类似地,L=3对应cc-pVTZ、L=4对应cc-pVQZ...由于式中包含E_SCF(∞)、A、α这三个常数,所以原则上需要按基组由低到高依次计算三次能量才能获得它们。显然最便宜的选择就是用DZ/TZ/QZ来外推,此时结果已经很准了,一般没必要用更高贵的诸如TZ/QZ/5Z来外推。

为了避免外推时的麻烦,对于不同的基组序列有人事先拟合了α参数,因此实际外推时只需要计算两次能量就行了。以下数值取自ORCA手册2.9版6.1.3.4节
对于cc-pVnZ系列,L=2/3的外推α=4.42,L=3/4的外推α=5.46
对于def2系列,L=2/3的外推α=10.39,L=3/4的外推α=7.88

假定α是已知量,我们就可以直接写出E_SCF(∞)的表达式。设
E_SCF(N)=E_SCF(∞)+A*exp(-α*√N)
E_SCF(M)=E_SCF(∞)+A*exp(-α*√M)
2式可写为A=[E_SCF(M)-E_SCF(∞)]/exp(-α*√M),代入1式,经过整理得:
E_SCF(∞) = [ E_SCF(M)*exp(-α*√N) - E_SCF(N)*exp(-α*√M) ] / [exp(-α*√N) - exp(-α*√M)]

由于HF、DFT形式上是基于单电子近似的,所以基组只要能够合理表现出单电子性质就够了,这对基组质量要求并不很高,因此能量随基组增大收敛得也比较快。SCF计算通常在3-zeta级别结果已经很好了,4-zeta级别就挺接近CBS极限了。因此对于SCF计算做CBS外推的意义不是非常大。另外顺带一提的是,单电子性质实际上并不需要很高角动量来描述,比如对于主族元素,角动量函数达到f就已经足够了,引入g及更高角动量函数对结果几乎不会有什么改进,远不如把计算量投入在增加低角动量的基函数上来得有用。
 

3 相关能的外推

后HF计算在HF基础上考虑了电子相关作用(主要是库仑相关)以改进结果。合理描述电子相关作用就要能够较好地描述出电子间的相关穴。由于这是双电子性质,明显比描述单电子性质对基组质量有更高的要求。尤其是相关穴的cusp(歧点)特征,对于非R12/F12方式的后HF计算,需要很大基组才能充分准确描述(高角动量函数起着关键性作用),这是由于通过单电子基函数乘积的方式来描述cusp效率很低所致。由于内在本质的不同,相关能的收敛趋势和HF能量的收敛趋势也有显著不同,因此在外推时需要分别进行,最后将CBS下HF能量和CBS下相关能加和得到总的CBS能量。
(注:虽然HF波函数能够在当前级别下精确表现Fermi相关,这是双粒子性质无误,但实际上它是Slater行列式的反对称化的要求而自然而然产生的,并非是靠大基组来表现的)

Klopper-Helgaker的形式是最常用的外推相关能的方法,理论依据也最明确,即:E_corr(L)=E_corr(∞)+A*L^(-3)

实际上,自旋相同电子间(singlet pair)的相关能和自旋相反电子间(triplet pair)的相关能收敛趋势也是不同的,前者随L^(-3)而收敛,后者随L^(-5)而收敛,严格的做法是分别拟合这两种相关能。但由于后者随L增加明显收敛得更快,因此通常在外推时不考虑后者,而直接用上面给出的外推公式。

Klopper-Helgaker公式只含两个未知量E_corr(∞)和A,因此只需要用基组序列中L相邻的两个基组各算一次能量即可。假设一个是L=M,一个是L=N,则
E_corr(M)=E_corr(∞)+A*M^(-3)
E_corr(N)=E_corr(∞)+A*N^(-3)
第二个式子可写为A=[E_corr(N)-E_corr(∞)]/N^(-3),代入到第一个式子,得
[E_corr(M)-E_corr(∞)]*M^3=[E_corr(N)-E_corr(∞)]*N^3
整理一下,得到了明确的CBS相关能计算公式:
E_corr(∞) = [E_corr(N)*N^3 - E_corr(M)*M^3] / (N^3 - M^3)

由于相关能随基组收敛得较慢,所以外推是很有必要的。选用基组序列中哪两个基组外推相关能是很关键的问题,直接影响外推结果准确度。对于实在算不动的情况,就用DZ/TZ来外推,但由于DZ级别质量太低,不能指望外推出的结果真的能接近CBS极限,实际上也就是比TZ强一点罢了,甚至更糟也有可能。通常都用TZ/QZ来外推,这样的结果比较可靠,但也仍不能指望达到CBS极限。可以期望外推出的结果精度会比外推时所用的最大基组高出一级,如TZ/QZ推出来的约是5Z的精度(视具体情况而定,可能比5Z更好或更差),这样的精度已经相当高了。

Truhlar等人的相关能外推公式更为广义:
E_corr(M)=E_corr(∞)+A*M^(-β)
当β=3,就回到了Klopper-Helgaker公式。为了使外推更准确,对不同基组序列应拟合不同的β值,下面列举一些(更多的见ORCA手册里基组外推部分)
对于cc-pVnZ系列,L=2/3的外推β=2.46,L=3/4的外推β=3.05
对于def2系列,L=2/3的外推β=2.4,L=3/4的外推β=2.97
由于对于L=3/4的外推β很接近3,所以此时或更高级别的外推建议都用3,理论依据更强。而L=2/3的外推时β建议用拟合值。

另外也有别的小众一些的外推方法,比如Varandas提出的,见JCP, 126, 244105 (2007),结果和上面介绍的没什么区别,形式还更为复杂,这里就不再多提了。需要外推至CBS都是高精度计算的场合,因此必须计算方法精度也必须比较高才能让总能量误差较小。然而高级别方法结合大基组的计算往往无福消受,这时候可以运用一些技巧。例如,一种常见的近似获得CCSD(T)/CBS下相关能的方法是:
E_corr(∞)_CCSD(T) ≈ E_corr(∞)_MP2 + [ E_corr(Small)_CCSD(T) - E_corr(Small)_MP2 ]
也就是说在MP2级别下将相关能外推到CBS,然后加上小基组下CCSD(T)和MP2相关能的差值,就近似获得了CCSD(T)下外推到CBS的相关能,这比起直接在CCSD(T)下外推省时多了。当然,这个近似成立的条件是CCSD(T)和MP2的相关能差异受基组影响不大,实际上这个条件通常是合理的。
 

4 实例:N2的能量外推

这里以N2作为具体例子来外推HF能量和MP2相关能。以下是cc-pVDZ/TZ/QZ/5Z/6Z计算的能量(hartree),键长为1.100314埃,Δ代表相关能部分。
E_HF(DZ)=-108.953748406
E_HF(TZ)=-108.982940693(-0.029192287)
E_HF(QZ)=-108.990532392(-0.007591699)
E_HF(5Z)=-108.992208035(-0.001675643)
E_HF(6Z)=-108.992531280(-0.000323245)
ΔE_MP2(DZ)=-0.3070859654
ΔE_MP2(TZ)=-0.3743967513(-0.067310786)
ΔE_MP2(QZ)=-0.3994430246(-0.025046273)
ΔE_MP2(5Z)=-0.4098026746(-0.010359650)
ΔE_MP2(6Z)=-0.4145085175(-0.004705843)

括号里是相对于上一个数值的变化,将它作图,如下所示



可见,HF能量收敛得明显要比相关能快得多,QZ级别的HF能量精度已经很高了,而相关能收敛得却仍不十分充分。虽然相关能占总能量比重很小,但是相关能计算的不准确却是总能量误差的主要来源。

利用第二节的HF能量外推公式,以及L=2/3时α=4.42、L=3/4时α=5.46的参数,可以得到以下结果,这可以方便地用Excel来做。括号里是相对于外推时用的最大基组下结果的差值。
E_HF(DZ/TZ->CBS)=-108.9924345(-0.014689896)
E_HF(TZ/QZ->CBS)=-108.9928198(-0.002287408)
通过比较可见,用DZ/TZ外推的结果颇接近直接用6Z计算的结果。而TZ/QZ外推的结果明显比起6Z更为精确,可以认为就是CBS极限的结果了。

利用第三节的Klopper-Helgaker相关能外推公式,可以得到以下结果
ΔE_MP2(DZ/TZ->CBS)=-0.402738135(-0.0283413837)
ΔE_MP2(TZ/QZ->CBS)=-0.417720035(-0.0182770104)
通过比较可见,用DZ/TZ外推的结果介于QZ和5Z之间(更接近于QZ)。TZ/QZ外推的结果比6Z稍好点。和外推HF能量相比,外推相关能对总能量的影响明显更大,这凸显了相关能外推的重要性。另外也看出,相关能外推的效果比起HF能量外推要弱,没法在较低级别基组下就一下子外推出很高级别基组下的结果。又由于相关能随基组收敛慢,所以如果对结果精度要求很高,外推相关能尽量要在TZ/QZ甚至更高级别下做,DZ/TZ外推的结果离CBS极限的差距还较大。

如果L=2/3外推时用Truhlar的公式,β用推荐的2.46,则
ΔE_MP2(DZ/TZ->CBS)=-0.413728888(-0.0393321367)
明显比Klopper-Helgaker公式外推的结果更接近于高级别基组外推的结果。所以L=2/3外推时推荐用Truhlar公式结合最佳的β参数。

想得到由TZ/QZ外推出的MP2/CBS能量,只要将-108.9928198和-0.417720035加和即可。

要注意的是,化学上感兴趣的并非总能量,而是相对能量,外推对结果的影响还是会在求差值过程中所大幅抵消的。求差值时,要么所有体系都做外推(而且要用相同基组外推),要么都不外推;如果有的外推有的不外推,结果误差会巨大。
 

5 考虑BSSE时的外推

分子间相互作用能一般这样计算,也就是所谓超分子方式计算:
E_AB=E_A-E_B
E_A和E_B是两个单体的能量。单体的几何结构用自由状态的还是复合物状态的依情况而定,这和本文内容无关。

想要得到CBS下的相互作用能,就是对E_AB、E_A和E_B都按照上述过程外推至CBS能量。但实际上问题没有这么简单。因为计算E_AB时需要考虑基组重叠误差(BSSE),详见http://sobereva.com/46。对于以色散作用主导的高精度计算考虑BSSE尤为重要,此时哪怕是aug-cc-pVQZ级别的计算,BSSE的影响仍不可忽视。通常就是用Counterpoise来解决BSSE,这使得复合物能量的外推过程变得略麻烦些,所以专门在这一节说说。

经过Counterpoise方式对BSSE校正后的AB能量E_AB'这样计算(下文中的符号带一撇皆代表做了BSSE校正):
E_AB'=E_AB+E_BSSE
E_BSSE = (E_A - E_A,bAB) + (E_B - E_B,bAB)
其中
E_AB:A、B基组下AB复合物的能量
E_A,bAB:A、B基组下A的能量
E_B,bAB:A、B基组下B的能量
E_A:A基组下A的能量
E_B:B基组下B的能量

对于后HF计算,BSSE校正能分为对HF能量的校正E_BSSE_HF和对相关能的校正E_BSSE_corr:
E_BSSE_HF = (E_A_HF - E_A,bAB_HF) + (E_B_HF - E_B,bAB_HF)
E_BSSE_corr = (E_A_corr - E_A,bAB_corr) + (E_B_corr - E_B,bAB_corr)

BSSE校正后的HF能量和相关能即为
E_AB'_HF=E_AB_HF+E_BSSE_HF
E_AB'_corr=E_AB_corr+E_BSSE_corr

以这种方式得到不同基组下的E_AB'_HF和E_AB'_corr,就可以按照前述方法直接外推出CBS下的复合物的HF能量和相关能。这样比起在不考虑BSSE条件下,即使用E_AB_HF和E_AB_corr来外推更为合理。如果在aug-cc-pVTZ/QZ级别外推,由于BSSE的影响不是很大,所以外推时不考虑BSSE问题也不大,但如果要求精度很高,仍建议将BSSE考虑进去。

虽然诸如Gaussian程序直接提供了counterpoise关键词,从而可以一步就得出E_BSSE,这对于HF/DFT很方便,但是对于后HF外推,必须得分别得到E_BSSE_HF和E_BSSE_corr,这就只能自行手算。我建议此时使用Ghost原子的方法手动依次计算A、B、A,bAB、B,bAB体系并提取各自的HF能量和相关能,然后计算E_BSSE_HF和E_BSSE_corr,而不要靠counterpoise关键词,否则从输出文件中找数据时容易弄混。如果整体有对称性,这样比起直接用counterpoise关键词在效率上也有好处,因为counterpoise关键词会对涉及的所有计算都关闭对称性,而手动算的话,计算E_AB、E_A、E_B时还能开着对称性以节约时间。
 

6 实例:考虑BSSE时外推氢气-氮气二聚体能量

下面就以氢气-氮气二聚体平行构型为例来说明怎么得到二聚体的考虑了BSSE的外推至CBS下的能量,外推通过aug-cc-pVDZ/TZ来实现,理论方法用CCSD(T),程序用G09 B.01。这是笔者在J. Mol. Model.19,5387的研究中涉及到的计算。

先在CCSD(T)/aug-cc-pVTZ下计算,以下是route section和坐标部分
计算E_AB:
#p ccsd(T)/aug-cc-pvtz
 N                  0.00000000    0.55705000   -0.42673900
 N                  0.00000000   -0.55705000   -0.42673900
 H                  0.00000000   -0.36881700    2.98717200
 H                  0.00000000    0.36881700    2.98717200
计算E_A:
#p ccsd(T)/aug-cc-pvtz
 N                  0.00000000    0.55705000   -0.42673900
 N                  0.00000000   -0.55705000   -0.42673900
计算E_A,bAB:
#p ccsd(T)/aug-cc-pvtz nosymm
 N                  0.00000000    0.55705000   -0.42673900
 N                  0.00000000   -0.55705000   -0.42673900
 H-Bq               0.00000000   -0.36881700    2.98717200
 H-Bq               0.00000000    0.36881700    2.98717200
计算E_B:
#p ccsd(T)/aug-cc-pvtz
 H                  0.00000000   -0.36881700    2.98717200
 H                  0.00000000    0.36881700    2.98717200
计算E_B,bAB:
#p ccsd(T)/aug-cc-pvtz nosymm
 N-Bq               0.00000000    0.55705000   -0.42673900
 N-Bq               0.00000000   -0.55705000   -0.42673900
 H                  0.00000000   -0.36881700    2.98717200
 H                  0.00000000    0.36881700    2.98717200

HF能量就是输出文件中SCF Done:后面的,CCSD(T)能量搜索CCSD(T)=就能找到,它减去HF能量就是CCSD(T)的相关能。计算结果如下:
E_AB_HF=-110.11354966
E_A_HF=-108.98075536
E_A,bAB_HF=-108.98078589
E_B_HF=-1.13304859
E_B,bAB_HF=-1.13305499
E_AB_corr=-0.43998843
E_A_corr=-0.39983017
E_A,bAB_corr=-0.39989990
E_B_corr=-0.03956792
E_B,bAB_corr=-0.03957422
故aug-cc-pVTZ级别下:
E_BSSE_HF=0.00003692
E_BSSE_corr=0.00007603
E_AB'_HF=-110.11351274
E_AB'_corr=-0.43991240

把上述计算在aug-cc-pVQZ下再算一遍,可得:
E_BSSE_HF=0.00001463
E_BSSE_corr=0.00003224
E_AB'_HF=-110.12069678
E_AB'_corr=-0.46017610
可见在aug-cc-pVQZ下BSSE的影响明显比aug-cc-pVTZ要小。

由于没有aug-cc-pVnZ的α参数,三点外推又麻烦,所以这里就不外推HF能量了,直接用aug-cc-pVQZ下的HF能量作为CBS下的。如果要外推也可以凑合用相应L=3/4的cc-pVnZ的α参数,和理想的aug-cc-pVnZ的α参数差异不会很大。

利用aug-cc-pVTZ/QZ下的E_AB'_corr,按Klopper-Helgaker公式外推得到CBS下的CCSD(T)相关能-0.47496313。将它与aug-cc-pVQZ下的E_AB'_HF=-110.12069678加和,最终得到二聚体的考虑了BSSE的CBS下的能量-110.59565991。令它减去外推到CBS的单体的CCSD(T)能量即是单体间相互作用能。