关于BSSE校正

本文转载自：http://blog.sina.com.cn/zhangsummer183

 

计算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的比例往往不小，甚至超过它，如果不进行校正，可能正负号都不对。

基组越小，单体间相互作用越弱则E_BSSE越大。E_BSSE会随基组趋于完备而逐渐减小至0，给基组加上弥散函数能有效减小E_BSSE。对于氢键复合物，由于相互作用不算很弱，所以用了带弥散的中上等基组后，不做BSSE校正无妨。而pi-pi相互作用，即便用了aug-cc-pVTZ级别的基组，BSSE仍然不很小，若是用aug-cc-pVDZ的话十分建议加上BSSE校正。

计算E_BSSE有多种方法，Gaussian用的是目前使用最广泛的Boys和Bernardi发展的counterpoise correction方法，应注意这种方法计算出来的只是实际E_BSSE的近似，并非完全严格、精确。设E_i为第i个分子在自身基组下的能量，E_i'为第i个分子在全部n个分子上的基函数都出现下的能量，则计算n个分子相互作用能中的E_BSSE = ∑[i]( E_i - E_i' )，E_BSSE必为正值。注意计算E_i与E_i'时的分子几何结构必须与i处在复合物时的一致，Gaussian会自动这样处理。

要计算A、B两个分子的相互作用能，在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的能量

最后会输出corrected energy (E_corrected)和BSSE energy (E_BSSE)。E_BSSE是校正能，E_BSSE = (E_A - E_A,bAB) + (E_B - E_B,bAB) ；E_corrected就是消除了因单体基组重叠造成的能量降低后的AB复合物能量，E_corrected = E_AB + E_BSSE。

BSSE校正后的真实的相互作用能这样计算：E_interaction = E_corrected - (E_A` + E_B`)，其中E_A`和E_B`分别是A和B在孤立状态下经过优化后的能量。注意，相互作用能是复合物（考虑了BSSE时）减去单体在孤立状态时的能量差，单体在孤立状态的结构与在复合物中的结构并不相同，尤其单体结构呈柔性、分子间相互作用比较强时差异会较为明显。在实际中为了方便、省时往往做这样的近似：E_A`=E_A，E_B`=E_B，即是说计算单体能量时使用处在复合物状态下的结构，这个能量在计算E_BSSE时就已经顺便计算了，所以省得再单独算，此时相互作用能公式也有另一种等价写法E_interaction = E_AB - E_A,bAB - E_B,bAB，可见E_A和E_B并没被用到，其实可以不用在counterpoise任务中计算。

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

若计算n个分子间的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'。E_corrected = E_complex + E_BSSE。相互作用能公式即E_interaction = E_corrected - ( E_1` + E_2` + ... + E_n` )，其中`的含义与前面所述一致，也可以将E_i`近似为E_i。计算过程中也用DCBS和MCBS来说明接下来将要计算的是哪项，但此时DCBS中的D的含义就不是具体指Dimer了，而是多分子复合物。

在分子内相互作用能计算时一般也要考虑BSSE，比如一条长链分子，计算一字形和字母C形的能量差就不能忽略这个问题。但由于两个片段属于同一个分子，需要特殊处理，而不能直接用上述方法。比如可以将分子人为地切成两段，悬键用比如H来封闭，然后适当调整两个片段，即让切断的部位离得远一些（否则这部分也会对E_BSSE产生贡献），可能造成BSSE的部位相对位置保持不变，然后将两个片段当成两个分子来同上获得E_BSSE。

如果有特殊原因，需要手动进行上述BSSE计算的每步操作而不直接用Gaussian提供的Counterpoise，可通过设定Ghost原子来实现，事实上这也是老版本不支持Counterpoise任务的Gaussian获得E_BSSE所必需的。只要把某个原子名后面加上-Bq就说明它是Ghost原子，如Na-Bq。Ghost原子也就是指这个原子照常有基函数，但是没有原子核和电子。这样靠每次手工修改输入文件分别得到E_AB、E_A,bAB、E_B,bAB、E_A、E_B也就能算得E_BSSE，显然没直接用Counterpoise关键词省事。

实际上无论计算何种分子，BSSE总是存在的，但常被忽略，尤其是分子内BSSE。计算柔性分子构象能量差人们往往还会注意这一点，但计算普通小分子时很多人会认为没有BSSE，实际上此时也存在着BSSE。如计算乙醇和它的同分异构体二甲基醚，都存在分子内BSSE，并且数值不一样。可以考察其中的氧原子，在乙醇中它连着一个C一个H，而二甲基醚中连着两个C，其它原子带来的BSSE显然是不一样的，这会影响到两种异构体的相对能量差，尽管微小到被忽略也无妨。哪怕对于乙醇和乙烷的甲基碳，虽然都是连着三个H和一个C，但是键长是有略微差别的，BSSE仍存在极微小的差异。再比如计算乙烷交错式与重叠式的能量差，显然两种构象BSSE也是不同的。有时这不再是个小问题，例如2006年报道的用MP2等后HF方法结合一些pople基组，优化出的苯环结构是略微弯曲的，显然不对，经过分子内BSSE校正后避免了这个问题。

分子内BSSE校正的方法不如分子间校正的方法广泛、成熟，虽然按上述方法切割成分子片段（比如把苯环切成一堆C-H）是可行的，但是如何切割过于任意，操作也麻烦，不易普适化。有人提出了atomic counterpoise(ACP)方法，将体系的总BSSE由每个原子的BSSE加和得到，即E_BSSE=∑[A]E_A,bA-E_A,bS，其中E_A,bA是A原子孤立存在时的能量，E_A,bS是A原子在周围Ghost原子基组下的能量。显然不可能计算每个原子的BSSE时都要带着全部其它原子的轨道，这样太慢，实际计算E_A,bS时可以只让A原子附近7埃以内的其它Ghost原子的轨道出现，再靠外的原子的轨道由于伸展不到A所以不用考虑，这样效率是比较高的。

Counterpoise方法有几个问题值得注意：
1 由于Ghost原子和真实原子是不同的，所以会破坏体系对称性。故计算高对称性二聚体时E_A,bAB和E_B,bAB的耗时都可能大于计算复合物的耗时。
2 Counterpoise有时有过校正问题，例如计算E_A,bAB时，A感受到的是完全“空闲”的B的基函数，能被A完全利用；而在复合物中A感受到的B的基函数已经有一定占据了，不能被A完全利用。由于两种状况B的基函数的状态不同，用E_A - E_A,bAB作为复合物中A的校正显得校正过头了，故所得相互作用能会被低估，有时结果还不如不用Counterpoise，对于最需要E_BSSE的小基组，过校正问题反倒容易出现。故Counterpoise该不该用其实没有明确答案，这和体系有一定关系。过校正容易在氢键体系出现，所以建议用个好点的基组就行了，而不用Counterpoise。对于范德华复合物，尤其是pi-pi作用，建议用aug-cc-pVDZ及以上基组+Counterpoise。
3 在Counterpoise下进行优化时由于无法获能量的得解析导数，只能通过数值方式进行，因此优化很慢。

附：高斯中counterpoise输入文件的写法

每个原子说明最后需要有一个整数说明这个原子属于第几个片段，例如：
# MP2/6-31G Counterpoise=2 Opt

Counterpoise with Cartesian

0,1 <-整体、片段1、片段2的电荷和自旋多重度都一样，只需写一次即可。
1 0.00 0.00 0.92 1
9 0.17 0.00 2.73 2
1 0.77 0.00 3.43 2
9 0.00 0.00 0.00 1

使用Z矩阵时写法有些特殊，下面是test562的输入文件：
#p uhf/3-21g counter=2 test scf=tight opt freq

Gaussian Test Job 562:
Counterpoise with different charges and multiplicites

0,1,0,3,1,2  <-总电荷，总自旋；片段1的电荷和自旋多重度；片段2的电荷和自旋多重度。
O,0.0,0.0,0.0,1  <-第一个原子必须用笛卡尔坐标。
O,1,ROO,2  <-片段2的原子。
X,1,1.,2,X3O
H,1,RO1H,3,HOX3,2,90.,0,1
H,1,RO1H,3,HOX3,2,-90.,0,1  <-片段1的原子。-90.与1之间必须写0来代表这是二面角
X,2,1.,1,52.5,3,180.,0  <-定位用的虚原子，不写所属片段号。由于是二面角，仍得写0。
H,2,RO2H1,6,H7OX,1,180.,0,2
H,2,RO2H2,6,H8OX,1,0.,0,2
 
ROO=2.98308
RO1H=0.94839
X3O=120.2827
HOX3=52.90868
RO2H1=0.94686
RO2H2=0.95173
H7OX=52.98178
H8OX=51.9632