談談BSSE校正與Gaussian對它的處理
談談BSSE校正與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色散校正？”》（http://www.shanxitv.org/413）、《亂談DFT-D》（http://www.shanxitv.org/83）。BSSE問題會影響勢能面，因此，BSSE效應不僅影響相互作用能，也會影響幾何優化、振動分析等任務的結果，只不過BSSE對相互作用能的影響顯著大于對其它問題的影響。計算強相互作用能（如化學鍵鍵能）時不要刻意考慮BSSE問題，見《計算化學鍵鍵能時若基組合理則考慮BSSE是完全多余的》（http://www.shanxitv.org/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，見此文的介紹：《大體系弱相互作用計算的解決之道》（http://www.shanxitv.org/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彌散函數，詳見《談談彌散函數和“月份”基組》http://www.shanxitv.org/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校正并計算分子間結合能的例子》（http://www.shanxitv.org/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在符號上是等價的）。