搜索
谈谈BSSE校正与Gaussian对它的处理文/Sobereva(3)Firstrelease2009-Aug-27Lastupdate:2016-Apr-4计算A、B分子间的弱相互作用能时，一般不能简单地通过E_interaction=E_AB-E(A)-E(B)来计算，因为E_AB能量相对于E(A)+E(B)的降低来自两方面，一方面是真实的A、B分子间的相互作用能，这是我们要求的；另一方面来自于A、B分子的基函数在复合物体系中重叠，相当于增大了复合物的基组而使E(AB)能量降低（严格来说前提是所用的理论方法是基于变分原理的），这个部分贡献如果也掺入E_interaction，则高估了相互作用能（即实际上结合能没有算出来的那么负），所以要去掉，它称为BasisSetSuperpositionError(BSSE)。所以双分子的相互作用能应该表述为E_interaction=E_AB-E(A)-E(B)+E_BSSE。对于弱相互作用，E_BSSE所占E_interaction的比例往往不小，甚至超过它，如果不进行校正，可能正负号都不对。基组越小，单体间相互作用越弱则E_BSSE越大（也有人认为对于很小的基组，由于基函数不容易延伸到相邻原子上，E_BSSE反倒不大）。E_BSSE会随基组趋于完备而逐渐减小至0，给基组加上弥散函数能有效减小E_BSSE。对于氢键复合物，由于相互作用不算很弱，所以用了带弥散的中上等基组后，不做BSSE校正无妨。而pi-pi相互作用，即便用了aug-cc-pVTZ级别的基组，BSSE仍然不很小，若是用aug-cc-pVDZ的话十分建议加上BSSE校正。计算E_BSSE有多种方法，Gaussian用的是目前使用最广泛的Boys和Bernardi发展的counterpoisecorrection方法，应注意这种方法计算出来的只是实际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的能量Gaussian最后会输出Counterpoise:correctedenergy(记为E_corrected)和Counterpoise:BSSEenergy(记为E_BSSE)。E_BSSE是CP校正能，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:doingDCBScalculationforfragment1。这里就是说明接下来计算的是E_A,bAB（假设A分子为fragment1），其中DCBS代表dimercenteredbasisset，说明以A、B分子为中心的基函数都出现，但是计算中并不纳入B的电子和原子核，这称为计算A的能量时添加了B的ghost轨道；如果是doingMCBScalculationforfragment1，就是要计算E_A，MCBS代表monomercenteredbasisset，计算中只出现属于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。如果有特殊原因，需要手动进行上述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总是存在的，但常被忽略，尤其是分子内BSSE。计算柔性分子构象能量差人们往往还会注意这一点，但计算普通小分子时很多人会认为没有BSSE，实际上此时也存在着BSSE。如计算乙醇和它的同分异构体二甲基醚，都存在分子内BSSE，并且数值不一样。可以考察其中的氧原子，在乙醇中它连着一个C一个H，而二甲基醚中连着两个C，其它原子带来的BSSE显然是不一样的，这会影响到两种异构体的相对能量差，尽管微小到被忽略也无妨。哪怕对于乙醇和乙烷的甲基碳，虽然都是连着三个H和一个C，但是键长是有略微差别的，BSSE仍存在极微小的差异。再比如计算乙烷交错式与重叠式的能量差，显然两种构象BSSE也是不同的。有时这不再是个小问题，例如2006年报道的用MP2等后HF方法结合一些pople基组，优化出的苯环结构是略微弯曲的，显然不对，经过分子内BSSE校正后避免了这个问题，有兴趣可以去看JCP,128,144108和JCTC,5,2574。分子内BSSE校正的方法不如分子间校正的方法广泛、成熟，虽然按上述方法切割成分子片段（比如把苯环切成一堆C-H）是可行的，但是如何切割过于任意，操作也麻烦，不易普适化。有人提出了atomiccounterpoise(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的耗时都可能大于计算复合物的耗时。2Counterpoise有时有过校正问题，例如计算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怎么用可以达到最佳效果，各有不同的说法。比如JCTC,10,49(2013)中作者认为对于aug-cc-pVDZ/TZ，在计算弱相互作用时counterpoise校正能只用一半效果最佳，比起不用counterpoise或用完整的counterpoise更好。而对于aug-cc-pVQZ及以上级别的基组，或者进行基组外推，则应当用完整的counterpoise，比不用counterpoise或只用一部分counterpoise更好。4使用Counterpoise时能量没有解析导数，只能通过有限差分以数值方式获得导数，因此优化，尤其是频率计算都非常慢，因此强烈不建议在counterpoise下进行优化和计算频率。实际上BSSE问题对于优化出的几何结构的影响很小。因此通常的弱相互作用计算，都是在不使用counterpoise的情况下先做优化，然后再带着counterpoise计算相互作用能。整个过程中基组应该带着弥散函数。5溶剂模型下考虑BSSE问题没有严格的办法，建议先在气相下做CP计算得到BSSE校正值，然后加到溶剂下计算的相互作用能上。附：Gaussian中counterpoise输入文件的写法每个原子说明最后需要有一个整数说明这个原子属于第几个片段。例如：#MP2/cc-pVTZCounterpoise=2CounterpoisewithCartesian0,1-整体、片段1、片段2的电荷和自旋多重度都一样，只需写一次即可。写全就是010101H0.000.000.921F0.170.002.732H0.770.003.432F0.000.000.001使用Z矩阵时写法有些特殊，例如：#B3LYP/6-31G*counter=2AdaptedfromtestJob5620,1O,0.0,0.0,0.0,1-第