谈谈能量的基组外推 | 宜武汇-ag真人国际厅网站

谈谈能量的基组外推

文/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),详见

。对于以色散作用主导的高精度计算考虑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)能量即是单体间相互作用能。

原文链接:http://sobereva.com/172

网络摘文,本文作者:15h,如若转载,请注明出处:https://www.15cov.cn/2023/08/27/谈谈能量的基组外推/

发表评论

邮箱地址不会被公开。 必填项已用*标注

网站地图