谈谈BSSE校正与Gaussian对它的处理

谈谈BSSE校正与Gaussian对它的处理
On the BSSE correction and its processing in Gaussian

文/Sobereva @北京科音
First release 2009-Aug-27 Last update: 2018-Jul-5

计算A、B分子间的弱相互作用能时，一般不能简单地通过E_interaction = E_AB - E(A) - E(B)来计算，因为E_AB能量相对于E(A) + E(B)的降低来自两方面，一方面是真实的A、B分子间的相互作用能，这是我们要求的；另一方面来自于A、B分子的基函数在复合物体系中重叠，相当于增大了复合物的基组而使E(AB)能量降低（严格来说前提是所用的理论方法是基于变分原理的），这个部分贡献如果也掺入E_interaction，则高估了相互作用能（即实际上结合能没有算出来的那么负），所以要去掉，它称为基组重叠误差（Basis Set Superposition Error, BSSE）。所以双分子的相互作用能应该表述为E_interaction = E_AB - E(A) - E(B) + E_BSSE。对于弱相互作用，E_BSSE所占E_interaction的比例往往不小，甚至超过它，如果不进行校正，可能正负号都不对。

一般来说，只有计算弱相互作用才需要考虑BSSE问题，相关信息看《谈谈“计算时是否需要加DFT-D3色散校正？”》（//www.umsyar.com/413）、《乱谈DFT-D》（//www.umsyar.com/83）。BSSE问题会影响势能面，因此，BSSE效应不仅影响相互作用能，也会影响几何优化、振动分析等任务的结果，只不过BSSE对相互作用能的影响显著大于对其它问题的影响。计算强相互作用能（如化学键键能）时不要刻意考虑BSSE问题，见《计算化学键键能时若基组合理则考虑BSSE是完全多余的》（//www.umsyar.com/381）。

用2-Zeta基组如6-31G*、def2-SVP算弱相互作用能时BSSE问题是极为显著的，必须通过适当方法计算出E_BSSE来解决，否则结果根本没法用。3-zeta基组如6-311G**、TZVP算弱相互作用能的BSSE问题比2-Zeta基组的轻，但依然必须解决，否则结果定量不准。E_BSSE会随着基组趋于完备而逐渐减小至0，即基组越大BSSE问题越轻。对于4-Zeta基组，BSSE问题就很不显著了。而给基组加上弥散函数也能显著减小BSSE问题，此时不专门计算E_BSSE往往也可以接受。对于氢键、卤键等静电主导的强度不很弱的作用，用了带弥散的3-Zeta级别的基组后，不做BSSE校正算出来的相互作用能也不错。而范德华吸引、pi-pi堆积作用能的计算，一般强度很弱，即便用了aug-cc-pVTZ级别的基组较大的带弥散的3-Zeta基组，BSSE问题仍然不可忽略，若是仅仅用aug-cc-pVDZ的话则更得明确考虑BSSE校正。

计算E_BSSE有多种方法，目前使用最广泛的是Boys和Bernardi发展的counterpoise方法，被Gaussian等很多程序直接支持。应注意这种方法计算出来的只是实际E_BSSE的近似，并非完全严格、精确，而本来也没有完全严格的方法来计算E_BSSE。设E_i为第i个分子在自身基组下的能量，E_i'为第i个分子在全部n个分子的基函数都出现下的能量，则n个分子间的相互作用能的E_BSSE=∑[i]( E_i - E_i' )。对于变分的方法，由于基组越大能量越低，因此E_i'比E_i更负，故E_BSSE必为正值。注意计算E_i与E_i'时的分子几何结构必须与i处在复合物时的一致，Gaussian做counterpoise任务时会自动这样处理。每个几何结构下都可以计算E_BSSE，数值不同。由于我们计算弱相互作用时用的复合物结构一般是经过优化的，因此一般做counterpoise计算的时候一般也是用的优化后的复合物结构。

若要计算A、B两个分子的E_BSSE，并且顺带计算校正后的相互作用能，在Gaussian中使用counterpoise=2关键字（可简写为counter=2）。会自动依次计算5个体系，依次输出以下能量：
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的能量

计算过程中会输出类似这样的语句Counterpoise: doing DCBS calculation for fragment 1。这里就是说明接下来计算的是E_A,bAB（假设A分子为fragment 1），其中DCBS代表dimer centered basis set，说明以A、B分子为中心的基函数都出现，但是计算中并不纳入B的电子和原子核，这称为计算A的能量时添加了B的基函数；如果是doing MCBS calculation for fragment 1，就是要计算E_A，MCBS代表monomer centered basis set，计算中只出现属于A分子的基函数。

如果你用的是G09早期版本或者更早的Gaussian，任务最后会输出"Counterpoise: BSSE energy"，这即是E_BSSE，即(E_A - E_A,bAB) + (E_B - E_B,bAB)。还会看到"Counterpoise: corrected energy"，记为E_corrected，这是消除了因单体基函数重叠造成的虚假能量降低后的AB复合物能量，E_corrected = E_AB + E_BSSE。

如果你用的是Gaussian09靠后期的版本，或者G16，counterpoise任务末尾会有类似下例的输出，看着更方便
Counterpoise corrected energy = -200.665575667261 //校正后的整体能量
BSSE energy = 0.001681932370 //BSSE校正能
sum of monomers = -200.659550058648 //两个单体在自己基组下的能量和
complexation energy = -4.84 kcal/mole (raw) //校正前的相互作用能
complexation energy = -3.78 kcal/mole (corrected) //校正后的相互作用能

若计算n个分子间的E_BSSE，则关键词为counterpoise=n，能量按如下顺序输出：E_complex, E_1', E_2'...E_n', E_1, E_2...E_n。E_BSSE = (E_1 - E_1') + (E_2 - E_2') + ... + (E_n - E_n')。计算过程中也用DCBS和MCBS来说明接下来将要计算的是哪项，但此时DCBS中的D的含义就不是具体指Dimer了，而是多分子复合物。

如果有特殊原因，需要手动进行上述counterpoise的每步计算，则可以通过设定Ghost原子来实现（比如很老的Gaussian版本不支持counterpoise关键词，要获得E_BSSE不得不手动计算）。只要把某个原子名后面加上-Bq就说明它是Ghost原子（如Na-Bq），即这个原子照常有它原本的基函数，但是没有原子核和电子。因此，比如要计算E_A,bAB，只需要在复合物的输入文件中把B片段的所有原子后面都加上-Bq使之成为Ghost原子即可。注意Bq原子若破坏了原有对称性，最好加上nosymm，否则可能中途报错停止。分别计算得到E_AB、E_A,bAB、E_B,bAB、E_A和E_B之后，就能按照前面的式子立刻算得E_BSSE。在能用counterpoise关键词的情况下这样手动做counterpoise显然没直接用counterpoise关键词省事，不过也有好处，就是手动做counterpoise的每个步骤比直接用counterpoise关键词可能更省时，尤其是高度对称的体系甚至能省几倍，因为用counterpoise关键词时所有任务都关掉了对称性，而手动做时，对单体、复合物仍然可以利用对称性来节省时间。

另一种如今也挺常用的解决BSSE的方法是Grimme提出的gCP，见此文的介绍：《大体系弱相互作用计算的解决之道》（//www.umsyar.com/214）。counterpoise任务耗时是复合物单点耗时的三倍多（对于DFT而言），而gCP校正则是“免费”的。但可惜到目前为止Gaussian还不支持。

涉及分子内弱相互作用时一般也要考虑BSSE。比如一条长链分子，计算一字形和字母C形的能量差就不能忽略这个问题。但由于两个片段属于同一个分子，不能直接用counterpoise方法处理。非要用counterpoise的话，可以将分子人为地切成两段，悬键用比如H来封闭，然后适当调整两个片段，即让切断的部位离得远一些（否则这部分也会对E_BSSE产生影响），原本可能造成BSSE问题的部位相对位置保持不变，然后将两个片段当成两个分子来通过counterpoise方法获得E_BSSE。而更优雅的做法是直接在计算时令基组带充分的弥散函数，或者不带弥散函数但是用gCP，这就省得切割片段了。

counterpoise方法有几个问题值得注意：
1 很多文章专门探讨counterpoise怎么用可以达到最佳效果，各有不同的说法。比如JCTC,10,49(2013)中作者认为对于aug-cc-pVDZ/TZ，在计算弱相互作用时counterpoise校正能只用一半效果最佳，比起不用counterpoise或用完整的counterpoise更好。而对于aug-cc-pVQZ及以上级别的基组，或者进行基组外推，则应当用完整的counterpoise，比不用counterpoise或只用一部分counterpoise更好。
2 使用counterpoise时能量没有解析导数，只能通过有限差分以数值方式获得导数，因此优化，尤其是频率计算都慢到吐血，因此强烈不建议在counterpoise下进行优化和计算频率！！！实际上BSSE问题对于优化出的几何结构的影响很小。因此通常的弱相互作用计算，都应当在不使用counterpoise的情况下先做优化，然后再带着counterpoise计算相互作用能。优化时候只要基组用3-zeta就够了，不加弥散函数一般误差也不大，但计算相互作用能时即便用了counterpoise，也强烈建议带着最低程度的弥散函数（即至少给重原子带一层s,p弥散函数，详见《谈谈弥散函数和“月份”基组》//www.umsyar.com/119）。如果在优化和振动分析过程中就想明确考虑BSSE问题（优化时只用得起不带弥散函数的2-zeta基组时必须如此），那么应该用gCP方法，ORCA程序里的gCP有解析一阶和二阶导数。
3 溶剂模型下考虑BSSE问题没有严格的办法，一般做法是在气相下做counterpoise计算得到E_BSSE，然后加到溶剂模型下以常规方式E(AB)-E(A)-E(B)计算的相互作用能上。
4 counterpoise有时有过校正问题，例如计算E_A,bAB时，A感受到的是完全“空闲”的B的基函数，能被A充分利用；而在复合物中A感受到的B的基函数已经有一定占据了，不能被A充分利用。由于两种状况B的基函数的状态不同，用E_A - E_A,bAB作为复合物中A的校正显得校正过头了，故所得相互作用能会被低估，个别时候结果还不如不用counterpoise。

附：Gaussian中counterpoise输入文件的写法

每个原子最后需要有一个整数说明这个原子属于第几个片段。例如：
# B2PLYPD3/jul-cc-pVTZ counterpoise=2

Counterpoise with Cartesian

0 1 //此例整体、片段1、片段2的电荷和自旋多重度都一样，只需写一次即可
H 0.00 0.00 0.92 1
F 0.17 0.00 2.73 2
H 0.77 0.00 3.43 2
F 0.00 0.00 0.00 1

若整体与片段的电荷和自旋多重度不一样，应该写：整体电荷，整体自旋多重度，片段1电荷，片段1自旋多重度，片段2电荷，片段2自旋多重度。

从G09开始，坐标部分也可以如下这么写来设置片段，效果和上面一样。
H(fragment=1) 0.00 0.00 0.92
F(fragment=2) 0.17 0.00 2.73
H(fragment=2) 0.77 0.00 3.43
F(fragment=1) 0.00 0.00 0.00

如果你是ORCA程序的用户，做counterpoise任务最方便的做法是通过Multiwfn产生输入文件，这里给了明确的例子：《在ORCA中做counterpoise校正并计算分子间结合能的例子》（//www.umsyar.com/542）。

补充：由于有网友看过此文后对BSSE的认识仍然有误，因此我将对他的回复也贴到了这里（修改了符号以与本文一致），希望读者看了之后能对BSSE问题根源了解得更明白一些。为简化讨论，这里假设不考虑单体在复合物状态下结构的改变。

对于计算相互作用能，E_interaction = E_AB - E_A - E_B这个式子没错，这是由结合能的定义直接得到的。而必须指出的是，在实际计算中，E_interaction = E_AB,bAB - E_A,bA - E_B,bB这个式子中（b代表基组）， E_AB,bAB作为复合物能量是对的，E_A,bA和E_B,bB作为单体能量也是对的，但是，这么计算相互作用能，在基组完备性差的时候，是绝对错误的！因为计算复合物和计算单体时不在相同的基组级别下。例如，A单体中靠近B单体的原子的基组完备性在复合物状态下高（因为B的基函数侵入此原子空间），而在A单体单独计算时低，这就导致能量在求差值时基组误差无法被抵消，这就是BSSE的根源。只有在基组完备的情况下，或者计算复合物和计算单体时都在相同基组级别下才能直接通过计算出的复合物能量减去单体能量来计算结合能。即正确的写法是E_interaction = E_AB,bQ - E_A,bQ - E_B,bQ，Q可以代表平面波基函数，也可以代表一套原子中心的基函数，比如可以指AB的基函数，这正是为什么本文中给出了这样的式子：E_interaction = E_AB - E_A,bAB - E_B,bAB，即计算复合物和单体的能量时都在AB基组下（文中E_AB和E_AB,bAB在符号上是等价的）。