談談量子化學研究中什么時候用B3LYP泛函優化幾何結構是適當的

談談量子化學研究中什么時候用B3LYP泛函優化幾何結構是適當的

First release: 2020-Jun-11 Last update: 2020-Nov-17

1 前言

量子化學研究中什么場合用什么DFT泛函合適，在《簡談量子化學計算中DFT泛函的選擇》（http://www.shanxitv.org/272）中我已經說過了。量子化學研究中最常見的兩類任務就是單點能的計算和幾何優化。B3LYP這個1994年提出的老泛函用于計算能量方面已經嚴重落伍了，不管是什么體系的單點能計算，在主流量子化學程序里總有結果更好甚至能吊打B3LYP的泛函可以用，對有機體系堅持用B3LYP算能量甚至可能導致文章被拒，見《堅持使用B3LYP算有機體系的人的下場》（http://bbs.keinsci.com/thread-12773-1-1.html）。

然而，對于幾何優化來說，B3LYP至今依然有其實用價值，絕不是說在精度方面沒有比它更好的選擇，而是因為用B3LYP仍然有一定實際好處：
1 在雜化泛函里，至少對于量子化學領域用得最多的Gaussian來說，B3LYP的計算速度是最快之一，這得益于B3LYP泛函的形式較為簡單。相對于形式復雜的M06-2X，以及wB97XD、CAM-B3LYP等范圍分離泛函，B3LYP明顯便宜很多。
2 B3LYP對積分格點依賴性低，因此在不是很高質量的積分格點下勢能面就已經是較為平滑的，因此可以用比較便宜的積分格點（因此對于Gaussian 16用戶，我一般都建議使用int=fine關鍵詞把昂貴的ultrafine積分格點降低到fine來明顯節約時間）。相反，較現代的泛函大多是meta或meta雜化類型的泛函，對積分格點更為敏感，尤其是M06-2X等一些明尼蘇達系列泛函對積分格點更是相當敏感，若積分格點質量不夠高，容易優化不收斂、收斂了之后有小虛頻。對于M06-2X優化某些弱相互作用體系，連Gaussian中int=ultrafine檔次的格點都不夠用。如果對積分格點方面的知識不懂的話，看《密度泛函計算中的格點積分方法》（http://www.shanxitv.org/69）。
3 B3LYP歷史悠久，有海量量子化學研究的文章都用的是B3LYP優化結構，而且在優化大多數體系結構方面并沒有其它泛函比它有特別顯著優勢，因此用B3LYP比較不容易被審稿人質疑。

在本文，就具體對于幾種情況說一下B3LYP什么時候可以用、什么時候不適合用，以及為什么。

2 弱相互作用體系

簡單來說，涉及弱相互作用的體系絕對不要用B3LYP優化，否則如今極容易被審稿人批評；然而只要給B3LYP加上諸如DFT-D3等主流形式的色散校正，就可以萬事大吉，內行人一般都不會質疑。目前最常用的是就是B3LYP結合DFT-D3(BJ)色散校正，稱為B3LYP-D3(BJ)，在幾乎所有主流量子化學程序里都支持。如果不了解色散校正，看以下三篇
談談“計算時是否需要加DFT-D3色散校正？”
http://www.shanxitv.org/413
DFT-D色散校正的使用
http://www.shanxitv.org/210
亂談DFT-D
http://www.shanxitv.org/83

這里所謂的弱相互作用體系就是指弱相互作用可能明顯影響體系構型/構象的體系。例如，DNA雙螺旋結構當中，平行的堿基對之間有顯著的氫鍵作用，每個堿基與它上、下方的堿基之間有顯著的pi-pi堆積作用，都屬于弱相互作用。更具體來說，B3LYP描述弱相互作用失敗是因為其描述弱相互作用當中的色散作用完全失敗，下面這張出自Chem. Rev., 116, 5105 (2016)的圖對這點體現得很清楚：

可見對于完全靠色散作用結合的Ar二聚體，B3LYP都給不出極小點結構，而加了DFT-D3色散校正后極小點位置就和精確位置很接近了。

其實B3LYP本身對于描述靜電主導的弱相互作用定性正確，優化也可以給出靠譜的結構。氫鍵有強有弱，本質各有不同，見《透徹認識氫鍵本質、簡單可靠地估計氫鍵強度：一篇2019年JCC上的重要研究文章介紹》（http://www.shanxitv.org/513）。如文中所示，對于比如水二聚體中的氫鍵來說，靜電作用占主要成分，所以B3LYP優化水二聚體的結構實際上是沒有太大問題的。但由于DFT-D3是“免費”的校正，無需額外耗時，而且加了之后對各種類型弱相互作用體系的優化都能變得很理想，對于優化靜電主導型的氫鍵作用體系還能有定量改進，所以不加DFT-D3完全說不通，也沒有理由如今還直接用B3LYP來優化這種體系。

下面看J. Am. Chem. Soc., 130, 16055 (2008)里面的一個例子，對比了TPSS（右圖）和加了DFT-D2形式的色散校正后（左圖）優化后的結果，這和B3LYP與B3LYP-D3(BJ)的情形類似

由圖可見，考慮色散校正后優化后的DNA結構和實際很相符，而直接用和B3LYP一樣對色散作用描述能力為零的TPSS優化的話，雖然氫鍵距離倒是沒明顯問題，但由于本質是色散作用的pi-pi堆積完全不能描述，使得堿基間的間距完全錯亂。

前一陣子思想家公社QQ群里有人問了如下問題，被我當做反面教材放到了北京科音基礎量化班的ppt里了，讀過前文的讀者想必知道他遇到的問題是因為什么了。如今還直接用B3LYP優化這種體系完全是缺乏量子化學計算最基本常識的表現。

有些初學者由于缺乏理論化學直覺，不知道什么時候弱相互作用會明顯影響體系結構，索性優化的時候就都加上DFT-D3色散校正就完了，有益無害，還不用擔心審稿人會質疑。若真碰上外行審稿人質疑你，屆時也總能解釋得通。

判斷什么時候弱相互作用可能影響構型/構象其實是很容易的事，比如《使用Molclus結合xtb做的動力學模擬對瑞德西韋(Remdesivir)做構象搜索》（http://bbs.keinsci.com/thread-16255-1-1.html）一文里對瑞德西韋做構象搜索，一看體系結構就知道有很多可旋轉的鍵，故構象非常豐富，不同構象下原子間相互接觸情況不同，顯然色散作用表現與否可能顯著影響構象的結構，更嚴重影響相對能量，因此此文里用了B3LYP-D3(BJ)來優化。一般來說，只要是做構型搜索、構象搜索，我都建議用B3LYP-D3(BJ)作為最終優化的級別，結果合理，相對來說又便宜又不容易出虛頻。

對于陰陽離子作用的情況用B3LYP優化是否需要考慮色散校正看具體情況。比如乙酸鈉，鈉離子和乙酸根的吸引相互作用幾乎完全是靜電和極化作用，B3LYP本身就能很好描述，因此色散校正可加可不加。而對于離子液體，雖然靜電作用占陰陽離子之間作用的絕大部分，但畢竟這類離子并不算小，特別是有的還帶有柔性的鏈，因此色散作用對結合能的影響不僅不可以忽略，還可能影響離子的構象，因此應當加上色散校正。

值得一提的是，有一篇JCTC上的文章對不同泛函優化弱相互作用體系做了橫測，發現B3LYP-D3(BJ)優化精度幾乎是最高的。B3LYP-D3(BJ)做優化不僅便宜，即便從精度角度來看，用B3LYP-D3(BJ)優化也是極為理想的選擇，詳見《證明B3LYP-D3(BJ)非常適合優化弱相互作用體系的一篇文章》（http://bbs.keinsci.com/thread-19495-1-1.html）。

3 主族元素構成的體系

這一節說的有機體系特指不含大共軛特征、不顯著涉及弱相互作用的情況。

用B3LYP優化這類體系一般都是能得到合理的結果的，有大量文獻可以提供支撐。比如J. Chem. Theory Comput., 8, 2165 (2012)選了一批由兩、三個原子構成的五花八門的體系，在普通泛函（即雙雜化以外的泛函）里B3LYP做優化、振動分析結果幾乎是最好的，比M06-2X明顯更好。根據J. Chem. Theory Comput., 12, 459 (2016)的測試，對于普通有機體系，在普通泛函中非常常用的PBE0表現得幾乎最好，B3LYP雖然不是拔尖的，但也在很合理的范圍內，表現得不錯。在J. Phys. Chem. Lett., 11, 9957 (2020)里作者定義了geometry energy offset (GEO)參數從結構誤差對能量影響的角度考察了一批泛函優化小分子的精度，發現B3LYP的表現基本上是最好的。

正因為B3LYP做優化比較穩，在一些高精度熱力學組合方法里都用B3LYP來優化。比如G4、G3//B3LYP、CBS-QB3都是基于B3LYP來優化的。在J. Phys. Chem. A, 121, 4379 (2017)里提出的組合方法里是用B3LYP-D3(BJ)做優化和振動分析，再結合DLPNO-CCSD(T)算的高精度單點和一些經驗的后校正來得到高精度熱力學數據的。

總的來說，即便被優化的體系不顯著涉及弱相互作用，我還是建議加上DFT-D3校正，畢竟這相當于增加了B3LYP的普適性。向外行人推薦優化級別的時候，或者給他們寫輸入文件模板的時候，建議也用B3LYP-D3(BJ)而非B3LYP。有人問我為什么《計算RESP原子電荷的超級懶人腳本（一行命令就算出結果）》（http://www.shanxitv.org/476）這個文章的腳本里用的優化方法是B3LYP-D3(BJ)而不是B3LYP，其實就是出于這個考慮，誰知道用這個腳本的人會優化什么樣的體系？都統一加上DFT-D3就省得解釋了。

4 含過渡金屬體系

B3LYP在大量文章里被用于優化過渡金屬配合物。這種做法可以接受，但如果對配位鍵鍵長很關注，B3LYP絕對不是好的選擇。在任何涉及到過渡金屬配位鍵鍵長的橫測中，B3LYP沒有一次表現特別優異過，基本都是表現中庸，見比如J. Chem. Theory Comput., 2, 1282 (2006)、J. Chem. Theory Comput., 13, 5291 (2017)里的泛函橫測。適合研究過渡金屬配合物的泛函在《簡談量子化學計算中DFT泛函的選擇》里寫了，比如TPSSh在優化過渡金屬配位鍵鍵長的測試文章里經常表現得非常出眾。

我每次講完北京科音量子化學培訓班的理論方法那部分、提到不太建議拿B3LYP算過渡金屬配合物之后，老有學員在答疑時問我他都已經用B3LYP對他的這種體系做了計算了，現在該怎么辦。我的意見是，既然都已經算了，這次就先這么算完，以后就別再用了。

如果你對含有過渡金屬體系的計算關注的重點只是配體部分，比如單純是配體部分發生的反應等，那么用B3LYP優化完全沒關系。

如果你算的體系牽扯到d族過渡金屬之間直接成鍵，此時靜態相關很強，一定要用純泛函優化，而不要用B3LYP以及任何其它的雜化泛函。比如J. Chem. Theory Comput., 8, 908 (2012)里研究了[V(C5H5)]2Pn體系，其中兩個V之間直接成鍵，實驗V-V鍵長是2.538埃，純泛函BP86優化出來是2.568埃，非常理想，TPSS、PBE等純泛函表現得也都不錯，而B3LYP優化出來為2.994埃，明顯和實驗差得非常遠。

5 大共軛體系

B3LYP不適合優化大共軛體系，這點很多人都不知道。

輪烯（annulene）是典型的大共軛體系，18個pi電子離域在整體。在Angew. Chem. Int. Ed., 43, 4200 (2004)中指出[18]annulene的C-C鍵鍵長應當是不均等的，而且體系是非平面的，而B3LYP算出來的則是C-C鍵鍵長都相等，是平面構型（D6h點群），明顯和實際不符。我用B3LYP/def2-TZVP優化出的無虛頻結構如下圖左側所示。下方右圖是我用wB97XD/def2-TZVP優化出的無虛頻結果，是D3點群，可見鍵長確實是不均等的，而且由于H-H的位阻作用導致體系輕微偏離了平面結構，比B3LYP明顯與實際情況更相符。

18輪烯體現的B3LYP優化大共軛體系不很合理的問題不是個例，在Chem. Mater., 29, 477 (2017)中作者Bredas專門提到：B3LYP functional suffers from a significant electron self-interaction error...The consequence is a strong tendency for B3LYP to overdelocalize the wave functions; this is the opposite of Hartree?Fock methods, which overlocalize the wave functions. Thus, in extended π-conjugated systems, B3LYP favors fully coplanar conformations。即B3LYP對大共軛體系，由于此泛函的自相互作用誤差(SIE)問題，傾向于把體系優化成平面的。它傾向于把鍵長描述得均衡化也同樣是由于SIE問題所導致的。而CAM-B3LYP、wB97XD這樣長程HF成份很高的泛函，以及M06-2X這樣全局高HF成份泛函，都沒有這種問題。關于不同泛函的HF成份可參看《不同DFT泛函的HF成份一覽》（http://www.shanxitv.org/282）。遇到有機大共軛體系，在優化的時候筆者比較建議用wB97XD，而用M06-2X也完全可以，不過它對積分格點比wB97XD更敏感，用ultrafine積分格點精度的必要性明顯更大。

最近，在筆者對18碳環的研究中，B3LYP也展現了同樣的問題。18碳環有in-plane和out-plane兩套pi共軛作用，每套pi共軛作用都是18個電子在整個體系上離域，相關討論見Carbon, 165, 468 (2020)和《談談18碳環的幾何結構和電子結構》（http://www.shanxitv.org/515）。對此體系不同泛函優化得到的C-C鍵鍵長如下所示

可見HF成份不是太高的泛函，比如B3LYP、PBE0、TPSSh，優化出來都是所有C-C鍵鍵長相等。相反，wB97XD以及其它一些HF成份較高的泛函優化出來的都是長-短鍵交替的情況，這不僅與高精度的CCSD優化的結果相符，而且也與AFM實驗上觀測到的一致。這體現出不僅是B3LYP，包括與它HF成份相近的泛函都因為SIE問題傾向于把大共軛體系共軛路徑上的鍵長描述得均衡化。

6 激發態優化

一個泛函對基態優化的好壞和它對激發態優化的好壞完全沒有必然聯系，激發態勢能面的描述難度也遠高于對基態的描述。B3LYP優化激發態表現一般，至少不如它優化基態穩妥。比如J. Chem. Theory Comput., 14, 3715 (2018)里面測試了48個泛函對41個中、小分子激發態的優化精度，雖然B3LYP優化激發態的平均誤差較小，但有個別體系定性不正確。而文中發現PBE0和wB97XD優化激發態至少對被測試的體系來說都是合理的，盡管平均絕對誤差比B3LYP更大。

眾所周知，B3LYP這種HF成份不高的泛函是無法合理描述電荷轉移激發態（CT態）的激發能的，見《亂談激發態的計算方法》（http://www.shanxitv.org/265）。如果不知道什么叫CT態，看《圖解電子激發的分類》（http://www.shanxitv.org/284）和《使用Multiwfn做空穴-電子分析全面考察電子激發特征》（http://www.shanxitv.org/434）。B3LYP對CT態的勢能面的描述也同樣是不合理的，所以不能用于這種激發態的優化。在Chem. Soc. Rev., 42, 845 (2013)的5.2節給出了一個例子，下圖中GH（global hybrid）代表的是B3LYP，RSH（range-separated hybrid）代表的是CAM-B3LYP，GS和EES分別代表基態和激發態。

由上圖可見，B3LYP和能夠合理描述CT態的CAM-B3LYP對于優化基態結構看不出什么差異，但是對激發態的描述就有定性不同了，B3LYP給出的兩個苯環近乎垂直的結構是嚴重錯誤的。在這篇文章中，作者還特意強調了用長程HF成份很高的范圍分離泛函優化CT態的重要性：for CT or Rydberg EES, it is mandatory to use range-separated hybrids (we recommend CAM-B3LYP or wB97X-D) to reach physically meaningful estimates. These latter functionals also circumvent the main problems when optimising the EES geometries having a strong CT character。如果你不知道什么叫范圍分離泛函的話，看《在Gaussian中自定義范圍分離泛函的方法》（http://www.shanxitv.org/550）。

順帶一提，優化基態的泛函和優化激發態的泛函并不需要必須一致，比如你用B3LYP優化基態而用wB97XD優化激發態，在原理上完全沒問題（只不過有些審稿人比較外行，這樣做有可能會被吐槽。怕被吐槽的話基態用wB97XD優化也完全可以，雖然耗時會更高，但比起優化激發態的耗時來說可以忽略不計）。

7 優化過渡態

B3LYP優化有機反應的過渡態一般來說是能給出合理結果的，即便到現在依然使用非常廣泛。對于過渡態的幾何優化精度的benchmark文章很少。雖然在Org. Biomol. Chem., 9, 689 (2011)中比較了一堆泛函優化出的過渡態的結構，但這篇文章質量不高，測試方式很不科學，也沒有給出多少有信服力、有價值的信息。作者在末尾給出的結論是B3LYP, still have an important role to play in the search for transition state geometries for organic reactions, and can generate results that are almost as reliable as much more expensive computational methods，但這個結論并沒通過文章中的數據來充分論證，不過如果有審稿人質疑你用B3LYP優化過渡態的合理性的話倒是可以把這篇文章搬出來用。

對有機體系的過渡態搜索，使用M06-2X通常更準確和可靠，但如果你想圖便宜用B3LYP，而且實際發現得到的結構完全滿足你的期望，那也不用再嘗試其它的了。對于非有機類型的反應，如果沒有其它信息幫助你判斷什么泛函可能比B3LYP更適合優化其過渡態，也大可先用B3LYP試試。

值得一提的是，對于一個反應，并非B3LYP能找出來只有一個虛頻的過渡態則這個過渡態就真的很合理。比如F+H2->HF+H這個自由基反應的例子，用不同的泛函結合def-TZVP優化出來的過渡態結構如下

雙雜化泛函B2PLYP-D3(BJ)優化出的鍵角是116.3度，F-H距離是1.600埃、H-H距離是0.765埃，這和J. Chem. Phys., 128, 034305 (2008)中在非常高級別MRCI/jul-cc-pV5Z下得到的113.7度、1.554埃、0.771埃相符得很好。上圖中M06-2X的鍵角偏大，但也定性正確，而B3LYP給出的則可以算是定性不合理了。如果繼續跑IRC，會發現B3LYP給的TS是錯的，它對應的IRC兩頭的結構都是HF+H，而M06-2X的TS對應的IRC則是正確的。

8 團簇

對于優化主族團簇，如果是碳團簇，用B3LYP優化毫無問題。對于優化其它類型的主族團簇，特別是硼團簇，PBE0或TPSSh有很大概率表現更好，見比如J. Phys. Chem. A, 123, 10454 (2019)。

B3LYP優化金團簇的結構非常差，見J. Phys. Chem. A, 121, 2410 (2017)的測試，在常見泛函里表現最好的是TPSSh，而PBE0雖然表現不突出但也明顯比B3LYP強。在《使用量子化學程序基于簇模型計算金屬表面吸附問題》（http://www.shanxitv.org/540）里涉及到了銀團簇來模擬銀表面，沒有用B3LYP-D3(BJ)優化也是有這方面的考慮。

如果你是優化d族過渡金屬的團簇，更是絕對不要用B3LYP，而應當優先考慮BP86、TPSS、SCAN等純泛函，這類似于前面提到過的V-V鍵的情況。

B3LYP優化金屬的晶格常數很垃圾，見J. Chem. Phys., 127, 024103 (2007)中的測試和分析討論，而改用PBE0或HSE明顯好得多。