辨析計算化學中的任務類型和理論方法
辨析計算化學中的任務類型和理論方法
文/Sobereva@北京科音 2023-Aug-5
0 前言
在網上答疑時，偶爾看到有初學者搞不清楚任務類型和理論方法，比如今天有人在思想家公社3號群問“md，aimd和qmmm的區別是啥啊？什么需求下會用到這三種呢？”，這三個詞明顯都不是一個層面的東西。此文就把計算化學中的任務類型和理論方法的關系明確一下，簡明扼要地介紹一下相關基本概念，做一個科普，以令計算化學零基礎者一次性搞清楚它們的關系。更多的計算化學的總覽性的知識在《北京科音初級量子化學培訓班》（http://www.keinsci.com/workshop/KEQC_content.html）里有介紹。
1 任務類型
有N個原子的非線型的體系有3N-6個內坐標描述原子間相對位置關系，在Born-Oppenheimer近似下體系的能量是依賴于內坐標的，也即能量是個3N-6維的函數，這被稱為勢能面（potential energy surface, PES）。計算化學領域有很多常見的任務（task）都是基于勢能面做的，即有了勢能面信息（能夠求得勢能面上任意位置的能量及其導數）就能做這些任務。下面羅列一下常見的這種類型的任務：
? 優化極小點：平時說的幾何優化（geometry optimization）一般也是指這個。此任務從一個給定的初猜結構開始，根據特定算法去尋找與之最近的勢能面上極小點的精確位置。在分子動力學程序如GROMACS里這種任務也被叫做能量極小化（energy minimization, em），只不過實際目的不一樣，em的目的主要是在動力學模擬之前釋放體系中可能存在的較大斥力（自行構建的初始模型里往往有原子間距離太近、斥力太大），免得動力學模擬一開始由于過大的原子加速度造成過大的速度而導致動力學模擬異常或崩潰。
? 優化過渡態：從一個給定的初猜結構開始，根據特定算法去尋找與之最近的勢能面上的一階鞍點的精確位置。過渡態可以視為反應過程中最有代表性的一個結構，可以由此判斷或驗證反應機理，利用它和相鄰極小點的能量求差可以得到反應勢壘，對于討論反應難易有關鍵意義，見《談談如何通過勢壘判斷反應是否容易發生》（http://www.shanxitv.org/506）。優化過渡態的算法介紹見《過渡態、反應路徑的計算方法及相關問題》（http://www.shanxitv.org/44）。
? 產生反應路徑：用于把基元反應對應的勢能面上兩個相鄰極小點之間最容易相互轉變的路徑產生，也相當于得到了一個軌跡并可以觀看對應的結構變化過程，過渡態是此路徑上能量最高點。在質量權重坐標下產生的反應路徑稱為IRC（intrinsic reaction coordinate），不考慮質量權重時產生的一般稱為MEP（minimum energy path）。具體算法很多，見《過渡態、反應路徑的計算方法及相關問題》（http://www.shanxitv.org/44）。反應路徑可認為是實際化學反應過程最有代表性的路徑，自然對于理解反應機理有重要意義，還可以用Multiwfn對反應路徑上各個點做波函數分析考察反應過程中電子結構的變化以獲得更深入的信息，例如《通過鍵級曲線和ELF/LOL/RDG等值面動畫研究化學反應過程》（http://www.shanxitv.org/200）。
? 振動分析：通常是在勢能面上的駐點（所有原子受力都為0的點）做的，用于得到相應結構下的振動頻率，可以用來計算熱力學量，公式看《使用Shermo結合量子化學程序方便地計算分子的各種熱力學數據》（http://www.shanxitv.org/552）介紹的Shermo程序手冊的附錄部分，還可以用來得到振動能級、振動波函數、檢驗幾何優化的準確性（根據虛頻判斷）、繪制振動光譜（參看《使用Multiwfn繪制紅外、拉曼、UV-Vis、ECD、VCD和ROA光譜圖》http://www.shanxitv.org/224）。
? 分子動力學（molecular dynamics, MD）：上面的任務都是“靜態”的任務，即不含時間這個維度。而分子動力學則引入了時間這個維度，可以模擬體系狀態隨時間的變化，能研究的問題與那些靜態的任務有極大的互補性，在北京科音分子動力學與GROMACS培訓班（http://www.keinsci.com/workshop/KGMX_content.html）里面有全面的介紹。分子動力學模擬過程相當于體系在勢能面上不斷運動，也相當于是對勢能面的一種采樣方式。雖然其名字叫分子動力學，但研究對象并不限于分子，比如無機固體、金屬團簇的動力學模擬也可以叫分子動力學。最常見的分子動力學的實現形式是BOMD（結合量子化學方法時另有CPMD、ADMP等），也相當于根據原子的位置、速度和受力按照經典牛頓運動方程演化原子的坐標和速度。還有其它一些特殊的動力學形式，如朗之萬動力學、耗散粒子動力學等，如今用得很有限。
? 蒙特卡羅：和分子動力學并列的另一種常見的對勢能面采樣的方法，適用場合相對更少，主要是在多孔材料吸附小分子、計算氣液共存曲線等問題上用得較多。不需要像MD那樣計算受力和速度信息，也沒有明確的時間的概念、無法像MD那樣嚴格地考察含時演化。
? 構型/構象搜索：勢能面上往往有很多極小點，能量最低的那個是全局最小點，可以視為體系最穩定的構型或構象對應的結構，其它的極小點對應亞穩的構型/構象。前面說的幾何優化只能收斂到與初猜結構最近的極小點，顯然這未必是最小點。如果想確保得到最小點（全局優化，global optimization），或者能量最低的一批構型/構象，就需要做構型/構象搜索，方法很多，比如免費的molclus（http://www.keinsci.com/research/molclus.html）就是很常用的實現工具。
2 理論方法
上一節介紹的那些任務在進行過程中，至少需要計算體系在不同結構下的能量。有的還需要計算受力（勢能面對核坐標一階導數矢量，或者說梯度，的負矢量），如幾何優化、分子動力學。有的還需要計算能量對核坐標的更高階導數，比如振動分析至少需要算二階導數（Hessian矩陣由之構成）。怎么獲得各種任務中涉及的能量及其對核坐標的導數，就是理論方法層面的問題了。下面簡要列舉一些常見的理論方法：
? 分子力場（molecular forcefield）：也稱為經典力場（classical forcefield）或分子力學（molecular mechanics, MM），或者經驗勢函數等等。雖然通常帶著“分子”倆字，但實際中不限于用于分子體系。這類方法一般使用簡單的模型描述體系（通常一個原子作為一個粒子，電子與原子核不分開考慮），并用形式簡單的數學函數（勢函數）描述原子間相互作用。例如計算原子間靜電相互作用時，分子力場普遍是給每個原子核位置分配一個點電荷（數值等于原子電荷），然后基于庫侖公式計算。由于分子力場的形式簡單，因此計算耗時極低。分子力場的一個關鍵的局限性是絕大多數分子力場都無法描述化學反應，因為成鍵關系在計算從始至終是固定不變的，是一開始就定義好的。描述化學反應問題主流的是后面提及的量子化學類型的方法，或者反應力場（遠比普通分子力場要貴、支持的程序要少）。分子力場另一個關鍵局限性是依賴于大量參數，一方面影響模擬精度，一方面決定能描述的體系，尋找和構建合適的參數往往不是易事。
? 機器學習勢（machine learning potential）：基本思想是通過機器學習的思想構造依賴于原子坐標的分子描述符（如距離矩陣、能量矩陣等）與能量之間的抽象關系，這種關系比上述傳統的分子力場用的勢函數形式明顯更為復雜。
以下提及方法都是量子化學（quantum chemistry, QC）范疇的，也可以叫基于量子力學（quantum mechanism, QM）的方法。
? 密度泛函理論（DFT）：是目前在量子化學、第一性原理計算領域用得最普遍的理論方法，因為其性價比非常高，交換-相關泛函（影響DFT計算精度的決定性因素）選擇得當的話可以得到很不錯的結果，足夠滿足絕大多數研究需要，而且精度顯著強于分子力場（除非力場是對某類體系專門訓練得特別精妙的力場），普適性強得多得多，但耗時高于分子力場N個數量級，體系越大N越大。
? Hartree-Fock（HF）：在DFT興起之前用得很普遍，如今已經完全過時。耗時不比DFT顯著更低（DFT計算時利用RI技術時反倒還可以明顯更快），而精度差得多。
? 后HF：是一大類方法的統稱，如CCSD(T)、MP2、CISD等。HF方法精度爛在于同時沒考慮到動態相關和靜態相關。后HF在HF基礎上額外把動態相關在一定程度上考慮了進來，但耗時也高了很多。動態相關、靜態相關的相關概念很有學問，在北京科音基礎量子化學培訓班（http://www.keinsci.com/workshop/KBQC_content.html）里有詳細講授，這里就不多提了。
? MCSCF：主要彌補HF缺乏對靜態相關的考慮。但由于幾乎沒有或很少考慮動態相關，所以結果也說不上好。MCSCF最常見的實現是CASSCF。
? 多參考方法：在MCSCF基礎上進一步把動態相關考慮進來，精度整體很好，普適性很強，但使用相對復雜且昂貴。典型的實現如CASPT2、NEVPT2、MRCI。
? 半經驗：一大類方法的統稱，都是對HF的簡化以巨幅降低耗時，耗時只是HF的微小零頭，整體精度也有不小降低，而且引入了依賴于元素的參數而只能用于有限的元素。典型的如AM1、PM3、PM6。一些較復雜的機器學習勢的耗時已經追平了半經驗。
? GFN-xTB：相當于納入了一定DFT思想的半經驗級別的方法，耗時和半經驗相仿佛，而整體精度和可靠性>=主流的半經驗方法，但跟DFT比還有很大差距。GFN-xTB從2017年開始發展，和歷史更早的DFTB有極大的共性。
前面介紹的HF、后HF、MCSCF、多參考方法里面都完全沒有引入經驗參數（極個別的除外，諸如CASPT2可以涉及shift參數），計算中的參數只依賴于物理常數，故它們被稱為“從頭算（Ab initio）”類型的方法。分子力場、半經驗、GFN-xTB由于引入了大量經驗參數，因此明顯不屬于從頭算。DFT算不算從頭算模棱兩可。如今尚未找到的理論上嚴格的交換-相關泛函是不含任何參數的，但如今只有近似的交換-相關泛函被提出，其中雖包含經驗參數，但整體相對較少，因此爭論DFT是否屬于從頭算沒意思，只是個說法問題。有的泛函的經驗性較小，有的經驗性則較高。如今也有的文章將HF、后HF、MCSCF、多參考方法等物理思想完全基于波函數的方法統稱為wavefunction theory (WFT)，在討論時使用WFT一詞時DFT就不會和那些方法混淆了。
還有一個詞叫“第一性原理（first-principle）”，其字義和從頭算一樣，但“第一性原理計算”如今的主流指代是“量子化學方法研究周期性體系”，不僅包括嚴格的從頭算方法，也包括DFT，因為這類計算最常用的就是DFT。DFTB雖然依賴很多參數，但DFTB計算周期性體系也往往被不計較地算作第一性原理計算。
由于有了“第一性原理計算”這個專門描述把量子化學方法用于周期性體系的計算的詞，因此平時說的“量子化學計算”、“從頭算方法計算”一般是指計算分子、團簇這樣孤立（非周期性）體系的情況。那些通過量子化學方法主要研究孤立體系的程序被稱為量子化學程序，如Gaussian、ORCA、GAMESS-US等，而通過量子化學方法主要研究周期性體系的程序被稱為第一性原理程序，如CP2K、Quantum ESPRESSO等。注意“量子計算”和量子化學計算又不是一碼事，前者強調的是利用在量子計算機上實現的量子算法，可以但絕不限于做量子化學問題的研究。
接下來再說很知名的QM/MM。這個方法是指使用量子化學方法以量子力學（QM）的思想描述體系的一部分，用分子力場以分子力學（MM）的思想描述體系的另一部分，并且恰當考慮兩部分之間的耦合。對于很大的體系，這顯然比起全都用量子化學方法描述要便宜得多得多（MM部分計算耗時相當于QM部分的九豬一毛），但代價是MM區域描述的精度比QM區要差，而且絕大多數力場不能描述這部分發生的成鍵/斷鍵過程。因此，必須恰當劃分QM和MM區，讓最需要精確描述的部分或者涉及化學反應的原子納入到QM區中，而起到較次要作用的“環境”原子納入到MM區中。
還有一個詞是ONIOM，也是劃分不同區域，不同方法著重描述不同部分，它可以把任意兩種理論方法相組合。如果把量子化學方法和分子力場相組合，特稱為ONIOM(QM:MM)，其用處和QM/MM類似但實現不同，描述時絕對不能搞混。一些相關討論見《要善用簇模型，不要盲目用ONIOM算蛋白質-小分子相互作用問題》（http://www.shanxitv.org/597）。
再說一下對電子激發態的描述。前面介紹的半經驗、HF、后HF、DFT等方法通過delta-SCF方式可以算激發態，但不夠普適。量子化學領域有很多專門的方法可以計算激發態，包括計算激發態的能量及其導數，如最常用的TDDFT、過時的CIS、較高精度的EOM-CCSD等，見《亂談激發態的計算方法》（http://www.shanxitv.org/265）。幾乎所有分子力場都是描述基態電子態的，但如果專門對激發態勢能面擬合力場參數，也可以用力場描述激發態，但已有的這樣的力場都是針對極其具體的體系的激發態擬合的，不具備普適性。機器學習勢同理。QM/MM也可以描述激發態，但僅限于電子激發完全發生在QM區的情況。
3 任務類型與理論方法的結合
第2節介紹的各種理論方法原則上都可以用于第1節介紹的各種任務。例如，分子力場、DFT、QM/MM等理論方法全都可以用于幾何優化、振動分析、產生反應路徑等任務。這些組合在程序實現上一般沒有什么需要特殊考慮的，一個程序支持什么理論方法，就能基于它們做什么依賴于勢能面的任務。
搞明白了前面介紹的名詞，從頭算動力學（ab initio molecular dynamics, AIMD）和第一性原理動力學（first-principle molecular dynamics, FPMD）顧名思義自然就知道是描述什么計算了。這兩個詞其實沒有嚴格的區分，實際中經常混淆使用，但從嚴格的習俗來說，AIMD和FPMD分別最適用于描述使用量子化學方法對孤立體系和周期性體系做分子動力學模擬。所以北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）里面講CP2K做基于GFN-xTB和DFT的分子動力學的地方標題用的是第一性原理動力學，但若說成AIMD也完全可以，不會造成什么誤會。
由于基于分子力場做分子動力學的文章數目和流行程度遠高于非常昂貴的AIMD/FPMD，因此平時說的分子動力學模擬默認是指的基于經典力場的分子動力學，這也往往被叫做經典分子動力學。但要注意“經典”這個詞在不同語境下有不同含義，有的場合下是特指把原子核當做經典粒子來考慮，此時基于BOMD形式做基于量子化學方法的分子動力學也是經典分子動力學。而與之相對的是“量子動力學”，強調把原子核以量子力學方式描述，原子核有對應的波函數。顯然，此處“量子動力學”和“基于量子化學做分子動力學”又不是同一個概念。
QM/MM MD指什么也不言而喻，顯然是用的QM/MM理論方法做分子動力學任務，這結合了QM/MM對體系描述的特點和分子動力學研究問題的能力。顯然這比所有原子都用QM描述（相當于AIMD/FPMD）要便宜得多得多，當然代價就是MM區域在動力學過程中描述得相對較糙，而且預期不會發生化學反應，同時還得想辦法去找力場參數。因此能否用QM/MM MD代替AIMD/FPMD當然得根據被研究的體系和問題判斷，盲目瞎用純粹是自討苦吃（吐槽一下，很多初學者沒有基本理論常識，一看文獻里有QM/MM就覺得好像高大上似的，就整天想著QM/MM而不知道其局限性，太naive了）。
文獻里的許多用詞并不標準，必須結合語境和常識理解說的是什么。比如可能有的文章沒明確用QM/MM MD這個詞而只說了QM/MM，若實際上在此基礎上跑了動力學，當然就是QM/MM MD。再比如有的文章表面上說是做的AIMD，但從計算細節描述來看作者把一部分當做MM來描述以節約時間，顯然實際上也是QM/MM MD。理解名詞必須有基本的應變能力。
如果用的理論方法描述的是激發態，那么前述各種任務就是對激發態做的。比如用TDDFT做激發態分子動力學，那就是動力學的每一步用的原子受力都是用TDDFT算的激發態的受力。順帶一提，有些文獻里聲稱做的是激發態動力學研究，但其實內容里根本沒做分子動力學，只是對激發態計算了勢壘，文中這么說只是由于勢壘和研究反應速率的反應動力學領域里的反應速率常數密切相關。另外值得注意的是研究激發態動力學往往需要考慮內轉換、系間竄越等勢能面切換的過程，考慮這些需要做非絕熱動力學，需要額外的算法，這方面的坑很深。所以并不是一個程序支持了計算激發態的理論方法和分子動力學模擬后就能完美、嚴格地做激發態動力學研究相關問題，情況遠沒那么簡單。