量子化學研究中切換泛函應當注意的問題

量子化學研究中切換泛函應當注意的問題

文/Sobereva @北京科音

 First release: 2018-May-2  Last update: 2021-Jun-20

量子化學研究者經常會用一種做法，即計算一種問題用A泛函，計算另一種問題則切換為B泛函。筆者在線上和線下答疑的時候，都常看到有人不恰當地切換泛函，故此文專門說說這個問題。

一定要明白，如果你要切換泛函，必須必須必須有非常正當的理由。一般有這三種正當理由：
(1) A泛函和B泛函適用的問題不同，因此在計算不同問題時用不同泛函。什么問題適合什么泛函，在此文里有匯總《簡談量子化學計算中DFT泛函的選擇》（http://www.shanxitv.org/272）。比如對某個有機大共軛體系，用PBE0優化基態結構，而用CAM-B3LYP計算電子光譜，這就是很合理的。因為PBE0對于優化有機體系基態的結構公認是適合的，在J. Chem. Theory Comput., 12, 459 (2016)的橫測中也充分證明了這點。但PBE0對于大共軛體系的電子光譜計算往往不理想（HF成分偏低所致），而CAM-B3LYP則沒這個問題，已被大量文獻用于大共軛體系電子光譜的計算且得到了不錯的結果。但對于基態幾何結構優化問題，用CAM-B3LYP則不及PBE0，文獻里也很少這么做。再加上CAM-B3LYP是范圍分離泛函，耗時也比PBE0高，因此基態不用CAM-B3LYP去優化。
(2) A泛函的耗時明顯比B泛函低，而精度也明顯比B泛函差。因此用A泛函計算那些對計算級別要求低而且計算量較大的任務，而用B泛函較準確地計算那些對計算級別要求高的問題，即好鋼用在刀刃上。最典型的情況就是用便宜的泛函優化，而用明顯更好且明顯更貴的泛函算單點能。例如計算弱相互作用，我們可以用B3LYP-D3優化結構，而用精度和耗時都顯著高一個檔次的雙雜化泛函PWPB95-D3去算單點能。
(3) 與研究內容、目的有關。比如要對比研究一系列泛函A、B、C、D對于某類分子的極化率、NMR、偶極矩之類問題的計算精度，那就沒必要用這四種泛函分別去優化結構和計算性質了，只需要用一個對當前體系合理的泛函去優化即可，之后計算性質的時候再切換成不同泛函去分別計算。

如果沒有正當理由，絕對不要切換泛函，否則只會在發文章時給自己找麻煩！內行的審稿人看到你在文章中不同問題用了不同泛函，幾乎總會思考一下你這么做的必要性，如果你切換泛函顯得沒有任何必要，就會產生疑點，被審稿人懷疑并詢問這么做的目的。如果你有正當理由，就如實解釋即可；但如果你沒有正當理由，沒法給出令審稿人信服的解釋，這就尷尬了。你如果說我就樂意這么干，會把人家惹怒，最終可能被拒；如果你想方設法東拼西湊點理由來應付，可能有的審稿人會被你糊弄過去，也可能有的審稿人依然覺得你的做法不妥但懶得繼續追究，而如果碰上特別較真的審稿人，那恐怕最后你只能把大量數據重算來使得所用泛函統一。因此，為了避免潛在的麻煩，如果不是有十分確切的原因、能確保被審稿人問及為什么換泛函時能給出有理有據的回復，就絕對不要切換泛函！之前在一次北京科音的量子化學培訓班里，有個學員是做過渡金屬配合物的，導師讓他用PBE優化結構，然后用M06算能量，他讓我提提看法。我問他為什么這兩步用不同的泛函，他卻說不出來理由。像這種情況，自己都不清楚為什么要換泛函，而且研究經驗也不多，就更不能切換泛函了。本來不切換泛函，審稿人萬一問起來為什么用那個泛函，可能回復起來都有點吃力；這回又切換泛函，吸引了審稿人的注意力，增加了被提問的幾率，到時候還得解釋那兩個泛函分別用在那兩個場合的原因，這無疑是自找麻煩。

哪怕你切換泛函有一定道理，但是理由并沒有充分到特別有必要這樣做，而且自己對理論方面的自信不足、沒能耐回擊審稿人的一些錯誤的看法，那么為了免得審稿人找碴，最好別切換泛函。考慮一種情況，B3LYP做優化和振動分析，M06-2X算單點，假設當前研究的是有機體系。這倆泛函這么搭配，確實有一定理由，因為M06-2X的耗時比起B3LYP要高一些（尤其是在振動分析上），M06-2X對積分格點要求比起B3LYP高很多，優化以及振動分析如果用M06-2X做的話，若不用較高積分格點精度（比如Gaussian里用ultrafine格點），出現虛頻、幾何優化難收斂的概率還是挺大的，而用較高積分格點精度則會進一步提升M06-2X的耗時。而且B3LYP結合DFT-D3后對大多數體系都能優化得不錯，見《談談量子化學研究中什么時候用B3LYP泛函優化幾何結構是適當的》（http://www.shanxitv.org/557）。因此，B3LYP做opt freq，單點用M06-2X，這么切換泛函確實有實際意義。但如果希望減少被某些審稿人comment的可能性，則建議把泛函統一，即opt freq和單點計算全都用M06-2X。對有機體系，M06-2X做幾何優化也是比B3LYP略微好一點點的，特別是存在分子內弱相互作用的時候（當然，對B3LYP可以通過加DFT-D校正彌補這點）以及存在大共軛的時候，何況M06-2X的計算耗時終究還是和B3LYP在一個檔次，差距遠沒有普通泛函和雙雜化泛函之間那么懸殊。當然絕不是說筆者不建議用B3LYP優化而M06-2X算單點這種組合，但這么做的時候應當在文章當中交代清楚切換泛函的理由，否則容易令某些缺乏經驗的讀者和審稿人起疑。至于前面提到的對過渡金屬體系，PBE優化，M06算單點，這種做法則是不可取的，因為根本沒有切換的意義，M06對這類體系的結果未必比PBE更合理，而且耗時還會提高，因此要用PBE就都用PBE，要用M06就都用M06，真要想結果更好，算單點時候應當用雙雜化。

補充：有本文的讀者表示，我看那知名的做有機反應計算的xxx大牛，人家就是B3LYP優化，M06L算單點，在高檔次期刊發了很多文章。言下之意泛函這么切換才是正確、最佳的。在本文上面一段我剛說完B3LYP優化，單點切換成M06-2X的理由就不是非常顯著，而M06L算單點，對有機體系本來就明顯不適合，更沒有什么依據表明結果比B3LYP更好，干嘛換泛函？這不是沒事找事么？人家名氣大，人家的做法就是對的？況且這大牛只不過是應用性計算做得多才出名，不是專門搞理論方法的，存在認識誤區、有不恰當的做法是很正常的。而且人家聲名在外，審稿人往往還是大牛的熟人，即便審稿人對其切換泛函有看法、覺得并不妥當，只要結果不存在明顯、嚴重不合理性，審稿人看在作者的面子上一般也不會追究。而你和那位大牛不同，你如果盲目效仿這大牛的并不科學的做法，就會有極高的幾率被審稿人盯上，還可能不依不饒，到時候你怎么對自己都說不清楚的切換泛函的做法予以合理解釋？就算你把大牛這么做的一些文章給審稿人看，審稿人也未必會放過你。要牢記：在不懂基本原理的情況下盲目效仿他人的做法是極度危險的！

對于Gaussian，雖然支持密度擬合方法加速純泛函的計算，但是由于做得非常差，一般也沒人去利用，所以在Gaussian里，純泛函和雜化泛函的耗時基本沒有區別。但是，在ORCA、Turbomole等程序里，由于DFT的密度擬合做得很出色，純泛函比雜化泛函在耗時上對大體系能低一個數量級以上（見《大體系弱相互作用計算的解決之道》http://www.shanxitv.org/214），因此如果你是這些程序的用戶，推薦用純泛函來充分節約高耗時步驟的計算時間。“xx 步驟使用純泛函是為了節約時間”也是一個切換泛函的非常正當的理由，足以說服懂行的審稿人。對于ORCA用戶，比如使用BLYP-D3(BJ)做優化和振動分析，而用B3LYP-D3(BJ)或M06-2X算單點，這是完全可以接受的而且很恰當的做法。但如果你是Gaussian用戶，這么切換泛函是不合適的，沒有正當理由（而且多數情況純泛函精度比雜化泛函差）。純泛函有一些顯著的軟肋，比如計算(超)極化率會明顯高估、計算激發能會明顯低估等等，算這些任務可千萬別為了圖便宜用純泛函。

有些人，算不同問題的時候切換泛函不是出于計算量的考慮，也不是出于計算精度的考慮，而是為了湊實驗數據，覺得計算結果與實驗越接近就越好、越容易發表。其實這種做法是不好的，畢竟理論計算并不是實驗的陪襯，沒必要刻意去一味地迎合實驗。盲目去讓結果貼合實驗，那做理論計算還有什么意義？每種泛函對于某類體系的某種問題的計算終究是有誤差的，只要選用的泛函在原理上合適、計算過程各方面要素都已經考慮周全，存在一定誤差應當坦然接受，只要文章中的理論計算能把要考察的問題分析、討論清楚就行了。眾所周知，泛函的HF成分越高，算出來的激發能就越高，于是有很多人利用這個趨勢，偏要挑一個激發能和實驗值最接近的泛函來算激發能，然后在文章里得意地說，本文計算結果和實驗相符非常好。這樣的做法其實真沒什么意思，而且如果你為了達到這個目的，用了一個對當前問題很少使用的泛函，那么你湊數據的意圖很可能會被有經驗的審稿人識破。比如算個有機體系的UV-Vis光譜，之前優化的時候你用的B3LYP（HF成份20%），但可能你發現B3LYP略微高估了激發能，于是嘗試使用HF成份更低（10%）的TPSSh泛函來算，并發現其結果誤差相對于實驗來說比B3LYP更小一點點，于是就在文章中用TPSSh算光譜了。這么做就給人感覺在動機上非常可疑，畢竟這個泛函極少用于有機體系的激發能計算，若是審稿人追究用這個泛函的原因，那么你很難給出正當理由，如果坦白交代就是單純為了湊實驗數據，這必然不會給審稿人好印象。因此，對當前這個例子，索性優化和激發能計算都用B3LYP，這就不會令審稿人覺得有什么疑點，就算誤差大一點也比起惹質疑好。還有一種挺常見的情況是研究者比較了幾種泛函，發現計算吸收光譜時候A泛函最接近實驗，但計算發射光譜的時候B泛函最接近實驗，如果因此就在文章中算吸收用A，算發射用B，那也是非常不合適的，這會讓審稿人覺得動機不純，畢竟不管是吸收還是發射都是激發能計算問題，用倆泛函是什么鬼？這種情況，就應該用一個對當前類型體系用得比較多的，而且實際發現算吸收和發射的誤差都在可容忍范圍內的泛函，別湊實驗湊得太露骨。

前些天有人在思想家公社QQ群里問“請教用m062x優化，b3lyp-d3做單點，這個配置合理嗎，能量的方法比優化的小了，主要覺得m062x算能量有點低估”（我不知道他算的什么體系，姑且視為是有機體系吧）。像這種做法，顯然非常不合理。有經驗的人都知道B3LYP-D3本身就比M06-2X便宜，通常精度還更低，這么做明顯違背一般邏輯。若碰到懂行的審稿人，要么會說你缺乏常識而批你一通，要么會懷疑你在刻意湊數據而要求你給個說法，說不定還會覺得你在計算時由于忽略了什么重要因素，導致本應該精度較好的方法M06-2X算出來的結果差，逼得你把M06-2X的數據隱藏起來，于是要求做更多解釋...假設對于當前問題，由于M06-2X自身存在的不足，可能算能量時表現得就是比B3LYP-D3要差，就是應當用B3LYP-D3才最合適，那你何不在優化時候也用B3LYP-D3，使得整個計算過程用的泛函統一？這樣從邏輯角度上顯得更自然，也免得因切換泛函被別人質疑。

關于加DFT-D3的問題，在《談談“計算時是否需要加DFT-D3色散校正？”》（http://www.shanxitv.org/413）中有專門討論，本文不再詳述。諸如優化用M06-2X而單點用M06-2X-D3這種做法是純屬自找麻煩的，雖然通常來說對M06-2X加不加D3影響很小，但你似乎都知道加D3會帶來好處于是在算單點時候加了，有什么正當理由在優化時候不加D3？何況是完全免費的校正，目前主流程序又不是不支持優化時候帶D3。而諸如優化弱相互作用復合物用B3LYP，算相互作用能的時候才改用B3LYP-D3(BJ)，那文章被秒拒一點都不冤枉。眾所周知B3LYP完全不能描述色散作用，而弱相互作用問題又和色散作用密切相關，B3LYP優化這種體系往往連結構的定性正確都沒法保證。

如果你用的后校正方法會明顯增加計算量，或者用于優化過程有顯著障礙，而且在優化時候其必要性又沒那么高，那么只在計算能量的時候加校正才是恰當的。比如解決BSSE問題的Counterpoise (CP)校正，用了它會導致失去解析梯度而令優化任務計算耗時極大地提升，而且從實際上看只要基組不是很小且理論方法合理，不用CP也不至于令優化結果誤差明顯，因此CP只適合在最后算能量的階段才考慮。而至于Grimme提出的gCP方法解決BSSE問題，它類似于DFT-D3一樣是免費的后校正，因此在用的時候也應該像DFT-D3一樣，優化和單點要加就都加，要不加就都不加。

眾所周知，幾何優化、振動分析、產生IRC這三個任務，用的計算級別必須嚴格相同，顯然更是絕對不能切換泛函，不管你有什么理由都不行。比如在ORCA 4.0里用M06-2X優化了分子，但是在ORCA 4.0里面M06-2X沒有解析Hessian而只能用數值Hessian，因此算頻率巨慢，故想改用B3LYP做振動分析，這顯然是絕對不行的。

凡是需要橫向對比的計算，在計算時也不要切換泛函。不同泛函給出的絕對能量沒有絲毫可比性，沒法互相求差，比如過渡態用M06-2X，極小點用B3LYP，它們的能量求差算出來的勢壘毫無意義。如果比較的問題是基于相對能量算的，比如算一批體系的反應能壘來比較誰的反應活性強，有的體系用wB97XD，有的體系用M06-2X，這也是不行的，因為不同泛函都存在著不同程度的系統性誤差，只有相同泛函的計算結果之間對比，從原理上講系統誤差抵消得才較好，從而能反映出研究的不同體系自身特點的不同。但為了遵從這個原則，有時候可能會遇到一些麻煩的情況。比如研究一個含過渡金屬的復合反應，有的階段的基元反應只涉及主族原子，而有的階段則涉及過渡金屬。眾所周知，M06-2X適合主族問題，而反應直接牽扯到過渡金屬時，相對于M06-2X而言，其同門的M06則更為適宜。那么計算這種體系的整個反應過程，能不能有的用M06-2X，有的用M06？這樣做是不合適的，給審稿人解釋起來也困難。像這種情況，索性用一個傾向性沒那么強的泛函，比如后來提出的明尼蘇達系列泛函MN15，它號稱對主族和過渡金屬考慮得比較均衡；或者用傳統的且普適性較強的B3LYP、PBE0之類。用它們來算，有可能整個反應中涉及主族的基元反應精度不如用M06-2X，涉及過渡金屬的基元反應精度不如用M06，這沒關系，最后都用精度高普適性也好的雙雜化算個單點就完了（ORCA里雙雜化計算由于能利用密度擬合技術，并不很昂貴）。

再看一個情況，優化一個體系基態用B3LYP，而做電子激發計算發現它的某個激發態是電荷轉移激發態，由于B3LYP算電荷轉移激發態不好（而且有文章證明其優化電荷轉移激發態結構也同樣不好），于是在優化這個激發態的時候改用了長程校正泛函wB97XD。這種做法正當不？畢竟這是同一個體系，又都是優化任務，卻用不同泛函來做，是不是這樣不妥？實際上這樣做是沒有明顯問題的，因為激發態和基態的研究首先就不屬于同一類問題，而且你切換泛函的理由也比較明確。而如果你就是希望統一，完全避免潛在的麻煩，那你基態和激發態都用wB97XD也完全可以，本身wB97XD算基態（至少對于有機體系來說）也是很適用的。注意如果牽扯到基態和激發態電子結構特征直接比較，特別是要求差的時候，那么切換泛函是不行的。比如，要通過Multiwfn繪制垂直躍遷過程的密度差圖，用wB97XD算激發態密度，而用B3LYP算基態密度，這樣的結果沒有任何意義，因為不同泛函下電子結構還是存在不小差異的，這么求的密度差就沒法如實體現出電子躍遷本身特征了，而會引入顯著的人為因素，甚至使得密度差圖是完全誤導性的。

上面雖然說的都是切換泛函，但上述討論也可以推廣到切換計算級別上，即DFT與波函數方法間的切換、基組的切換。比如優化時候用6-31G**，計算單點時候用def2-SVP，這就不合適了，本來都是同一檔次的基組，切換干啥？意圖何在？再比如，某有機體系，優化時候用M06-2X，算單點時候用MP2，這幾乎是作死行為。在有機體系能量計算上，MP2的精度明顯不及M06-2X，耗時還比M06-2X高很多，這和“算單點用的精度應該至少不低于幾何優化”這個量化計算常識完全相違背！如果你在優化和單點計算時都用MP2，可能有的內行審稿人只是在心里吐槽一下：都什么年代了還用MP2，這作者有沒有看過近10年的研究文章啊。但未必會寫comment；如果優化和單點都用M06-2X，有經驗的審稿人可能心想：其實單點用精度更高的雙雜化更理想，也未必會寫comment；而如果你優化用M06-2X算單點卻用MP2，這位審稿人看到作者連基本常識都沒有，絕對忍不住要在審稿意見里寫comment批一頓了。

以上討論還可以進一步推廣到切換溶劑模型、DFT積分格點等方面，同樣也是如果沒有特別原因，建議在整個計算中保持這些設定的統一。對于溶劑環境中的問題的研究，我們多數情況用隱式溶劑模型，而隱式溶劑模型一般也就增加耗時的15~20%左右（用線性響應溶劑模型時），即曰只要你在真空中算得動，你在溶劑模型下也肯定算得動。雖然大多數情況，溶劑效應對分子結構的影響是很小的，在氣相下優化就夠了，但是畢竟優化時考慮溶劑模型對耗時增加才那么丁點，何不在優化時也用溶劑模型，令整個計算保持統一，省得被審稿人找碴？而且，有些體系（特別是局部帶顯著電荷的體系），溶劑效應對結構影響非常明顯，甚至是定性的，你始終都帶溶劑模型，也消除了因為某些情況沒考慮溶劑效應導致結果不合理的風險。也有另一種情況，比如我要研究某個分子在8種溶劑下的垂直電離能，我根據經驗可以判定溶劑效應對這個分子結構的影響是可以忽略不計的，于是我只在氣相下優化這個分子，省得在各個溶劑下都優化一遍了，大大節約了計算量。這樣做是完全可以的，理由是非常正當的。

顯然，本文不可能把各種各樣理論研究過程中切換泛函的情況全都包羅萬象討論一遍，對其它的情況大家應當根據自己積累的理論知識判斷切換泛函的合理性、必要性，如果實在難以判斷，應求助真正懂行的人。

如何向審稿人證明自己用的泛函（或其它類型的理論方法）是合適的？如何選擇合適的泛函？這個問題和本文有一定聯系，在下面筆者進行了歸納（前述的《簡談量子化學計算中DFT泛函的選擇》里筆者對泛函的選用推薦也是根據以下方式來的）。
(1)看專門的對不同泛函計算某種問題的精度的對比測試（benchmark）文章。這類文章在主流的理論/計算化學刊物，如JCTC、JCC、PCCP、JCP等期刊上可謂多如牛毛，哪怕是那些很冷門的問題（比如轉動常數、原子核處的電場梯度）都有人做過橫測。從中你會知道對這類問題哪種泛函表現比較突出從而值得使用、哪些泛函精度爛或者用起來很危險。如果審稿人懷疑你用當前泛函的合理性，直接把相關的橫測文章塞給他即可。但要注意，哪怕是類似體系的類似問題，由于不同測試文章選取的測試集不同，往往出現兩篇文章結論大相徑庭的情況，這就需要你多看測試文章，看一篇往往不夠，起碼得看至少三篇，從而逐漸形成自己的認識，不要盲目被某一篇測試文章的結論牽著走。
(2)看大部分類似體系和問題的計算文章都用的什么泛函。如果某個泛函在某類體系的某種問題上的計算被廣為使用，成為了研究者們優先使用甚至默認使用的泛函，而且結果也確實不錯，那么這個泛函大可以放心去使用。雖然未必這個泛函對這類情況真的是最最理想的泛函，但至少一般不會出現定性錯誤的情況，比起用其它亂七八糟的泛函放心得多。如果一個泛函對這種情況的計算表現得經常很差，它也注定不會流行起來。這些相關計算的文章看多了，當審稿人問你為什么用當前的泛函算當前這類問題，塞給審稿人一大把較高水準的用這個泛函做的相關研究文章往往就夠了。
(3)進行實際對比測試。比如當前研究了一大批體系，如果其中某些（或者之前考察過的類似的體系）有實驗值（或者自己用高精度級別計算了其可靠的理論值），你同時用了好幾種泛函做了計算，發現其中某個泛函精度高且穩健性高，那么你用這個泛函去計算類似的體系一般也是靠譜的，也不怕被審稿人懷疑，因為論據很充分。
(4)根據研究經驗和理論分析。如果上述三類論據你都找不到，那就只能告訴審稿人，根據我的大量研究經驗，用xxx泛函是合理的（最好能給出已發表的文章）。或者，如果你有一定理論水平，你可以表示根據xxx泛函的特點（比如HF成份、擬合參數用的數據集的特點、對某些問題的計算表現出來的普遍特點等），此泛函用于此類計算在原理上是恰當的。但這種說辭主觀性太強，論據不很充分，審稿人未必會信。最終能否說服審稿人，這在很大程度上看你的語言表達能力和理論功底了。

強烈不建議大家效仿極個別文章，或者道聽途說，用一些冷門、不知名，而且缺乏充分論據以證明其適用性的泛函。比如前一陣，群里有個人問BVP86和BP86有什么區別，為什么用BP86而不用BVP86。這個問題很莫名奇妙，BVP86又沒什么特別的好處，又非常冷門，干嘛想要用它？用google學術搜索去搜BVP86，結果才不到300條，而搜索BP86，結果多達23000多條。雖然流行程度和泛函的好壞沒有必然關系，但你去用一個審稿人可能聞所未聞的極小眾泛函，負責一些的審稿人不可能不讓你解釋原因，而你又無從解釋，也沒法利用已發表的大量文章從側面體現合理性，顯然是在給自己挖坑，被要求重算或者最終被拒也是自找的。

PS：本文涉及的話題較雜，一些內容和討論與以下文章有關，有興趣可以看看
談談隱式溶劑模型下溶解自由能和體系自由能的計算
http://www.shanxitv.org/327
淺談為什么優化和振動分析不需要用大基組
http://www.shanxitv.org/387
各種后HF方法精度簡單橫測
http://www.shanxitv.org/378
大體系弱相互作用計算的解決之道
http://www.shanxitv.org/214
亂談DFT-D
http://www.shanxitv.org/83
DFT-D色散校正的使用
http://www.shanxitv.org/210
談談BSSE校正與Gaussian對它的處理
http://www.shanxitv.org/46
亂談激發態的計算方法
http://www.shanxitv.org/265
不同DFT泛函的HF成份一覽
http://www.shanxitv.org/282