各種后HF方法精度簡單橫測

1 前言

此文對現有的五花八門的后HF方法做一個集中橫測。限于時間精力，測試很粗糙、很不系統，任務非常簡單，但結論還是有一定普遍意義的。測試任務有四：
(1)CO在cc-pVDZ下的相關能。以CCSDTQ結果作為參考標準衡量其它方法對相關能的考慮。之所以用這么小基組是因為哪怕再大一點，用CCSDTQ都會很吃力。
(2)CNH在cc-pVTZ的相關能。以CCSDT(2)Q的結果作為參考標準。此測試檢驗用更大的基組、對另一種體系時，各種方法對相關能的考慮程度。
雖然考察不同方法對相關能的描述充分程度應當用起碼5個分子才有較可靠的統計結果，但只考慮CO和CNH，還是能說明一定問題的。由于化學上感興趣的都是相對能量，不同方法誤差抵消程度不同，所以不能只看對相關能的描述充分程度衡量方法精確程度，而必須考察實際化學問題。因此又做了一個非常簡單的實際問題的測試：
(3)HCN異構化過程的勢壘和能量變化，在cc-pVDZ級別下計算。所涉及的方法最高到CCSDT(2)Q。同時也順便考察了幾種泛函的結果。
(4)同上，但是在平時高精度計算常用的cc-pVQZ下計算，此時方法最高用到CCSD(T)。

本文數據使用NWChem 6.6和Gaussian 09 D.01計算。后HF中除了MP4(SDQ)、BD(T)、BD(TQ)、QCISD(T)、QCISD(TQ)由于NWChem不支持故用Gaussian來做以外，其它都是NWChem6.6做的（SCS-CCSD是CCSD模塊做的，其它是TCE模塊做的）。計算都用凍核近似，凡是兩個程序都支持的方法結果是精確一致的。DFT有的是NWChem有的是Gaussian算的。Gaussian計算時用4核并行，NWChem計算時是串行運行。耗時統計的是wall clock time。

2 CO在cc-pVDZ下的相關能

以CCSDTQ的相關能作為參考，得到的各種方法對相關能的考慮程度如下


如果單純以相關能論英雄的話（這樣做很片面），可見CISD爛得很，CISDT依然不咋樣，CISDTQ雖然挺理想但耗時已經高得要命，所以單參考CI算基態幾乎沒人用，這是眾所周知的。MP2作為最廉價的后HF方法也比CISD強（但僅限于當前體系，其它體系往往情況反過來）。CCSD比MP2強得多，和QCISD半斤八兩，這也是眾所周知的。CC2是CCSD的近似，不過從當前數據來看和CCSD相仿佛。MP3比MP2從此例看沒改進，還貴不少，這也是MP3根本沒人用的原因。到了CCSD(T)、QCISD(T)、BD(T)這個檔，結果就已經很精確了，很接近CCSDTQ對相關能的考慮，這三種方法的精度半斤八兩。CCSD[T]在原理上比CCSD(T)對CCSDT的近似程度更高，此例其相關能卻與CCSDTQ幾乎精確一樣，但這多半只是巧合。CCSDT(2)Q類似于通俗命名的CCSDT(Q)，精度幾乎和CCSDTQ基本沒有任何差異，但耗時低得多，因此可以作為極高精度計算比金標準CCSD(T)更高一檔的鉆石標準來用，很多對CCSD(T)尚不知足的文章都是這么做的，如JPCA,121,1693。QCISD(TQ)、BD(TQ)、LR-CCSD(TQ)-1比起QCISD(T)、BD(T)、CCSD(T)又通過微擾方式把四重激發算符考慮了，但結果并沒看出有什么改進，甚至精度還略降低了一絲。CCSDT的結果介于CCSD(T)和CCSDTQ之間，而不及CCSDT(2)Q，這和這個方法的定位正好一致。對比MP2-MP3-MP4-MP5，可以看到對相關能考慮按照這個順序是波折變化的，這是MP系列著名的特點。MP4明顯高估了相關能，MP5又明顯低估，都不理想，遠不及同檔次的耦合簇。理論方法名字后面有個橫杠再接泛函名字的是以DFT波函數為參考態做的后HF計算，可看到結果都不理想，和參考值偏離很大。

總之，以上數據顯示耦合簇的可靠性是不可動搖的。

下面是計算耗時的對比


注意圖中是對數坐標。對于對當前體系耗時很低的那些方法，耗時數據沒什么意義，因為統計時候的不確定度就已經達到耗時的數量級了，另外有些數據是Gaussian算的有些是NWChem算的，沒有可比性。圖中，耗時最高的是CCSDTQ，僅僅算一個CO耗時就高達6622秒，還是cc-pVDZ這么低級的基組，所以幾乎沒絲毫實用性。CCSDT還算有一定實用性，結合cc-pVDZ耗時82秒，比CCSDTQ低兩個數量級。在CCSDT基礎上微擾考慮四重激發算符得到的CCSDT(2)Q使耗時高了很多，達到394s，但仍比CCSDTQ耗時低一個數量級以上。CISDT和CCSDT耗時基本持平。

顯然，cc-pVDZ哪怕結合MP2都是沒有任何實際意義的，后HF方法必須靠高角動量基函數方能充分描述電子相關作用。由于當前例子表現出CCSDT(2)Q已經可以代替CCSDTQ作為金標準，而耗時低得多，因此不考慮CCSDTQ時我們可以用更大的基組算更大的體系，見下。

3 CNH在cc-pVTZ的相關能

一系列方法結合cc-pVTZ對CNH考慮的相關能程度如下，是以CCSDT(2)Q的相關能作為參照。


大趨勢和上一節類似。由于實際中感興趣的是相對能量，而以上相當于比較絕對能量，而且只是選取個別體系，所以對相關能的考慮相差在1%~2%左右的方法，并不能說明誰在實際問題中精度更好。但是對于相差比較大的情況，光靠相關能的描述程度還是一定程度可以反映方法精度的。

這一節的測試納入了SCS-MP2和SCS-CCSD。SCS-MP2眾所周知對于除了弱相互作用等問題，多數情況比MP2更好，不過從對相關能的考慮這點上倒沒有體現出來。SCS-CCSD更是夸張地高估了相關能。當前測試CC2對相關能考慮比CCSD差，這正常體現出CC2是對CCSD的近似。MP4依然高估相關能，但是從誤差大小來看比其近似版本MP4(SDQ)改進不少，而MP5又貴表現又很一般，沒什么實用性。同檔次CC、BD、QCI從相關能上來看依然沒體現出明顯差異。(TQ)比起便宜得多的(T)并沒體現出優勢。CCSD(2)T是一種在CCSD上以微擾方式考慮T3算符的做法，和CCSD(T)是同一檔但是做法上有一定差異，從本節和上一節的測試上看比最常用的CCSD(T)略遜色一絲。

當前體系在cc-pVTZ下，幾種方法的耗時為：
CCSDT(2)Q：9435s
CCSDT：885s
CCSD(T)：35.3s
CCSD：27.8s
CISDT：770s
CISD：31s
MP2：5.2s
QCISD(T)：8s（G09計算）
QCISD(TQ)：577s（G09計算）
BD(T)：19s（G09計算）
BD(TQ)：661s（G09計算）
MP5：541s（G09計算）

可見，結合有實際意義的cc-pVTZ，CCSDT(2)Q也就能算得動幾個原子的體系，用于極高精度計算很適合。相對來說，CCSD(T)的耗時只是它的零頭。CCSDT比CCSD(T)貴得多得多，但是從相關能數據來看改進甚微。CISDT比CCSDT便宜一點，畢竟原理和算法上簡單得多得多，但精度實在太爛，沒實用價值。(TQ)比(T)貴得不是一星半點，而是以數量級衡量。BD比同級別QCI和CC系列都貴，但通常并不帶來額外好處，再加上支持它的程序又少，所以也沒什么意義，并不流行。MP5則是又貴又臭，耗時直逼QCISD(TQ)。

4 HCN異構化過程，cc-pVDZ基組

這里比較一下不同方法在cc-pVDZ基組下計算的極為簡單的HCN->CNH異構化過程的勢壘誤差和反應能誤差，用的是電子能量。極小點結構、過渡態結構都是B3LYP/6-31G*優化得到的。誤差都是絕對誤差。勢壘誤差具體指的是正向勢壘誤差和逆向勢壘誤差的平均值。結果如下，順序是按照勢壘誤差從小到大排的，參考值是CCSDT(2)Q的結果。之所以是按照勢壘誤差來排序，是因為其挑戰性比計算反應能更大，誤差更難以抵消。


由圖可見，CCSD(T)表現極好，勝過所有其它方法，和CCSDT(2)Q的誤差僅有不到0.05kcal/mol。BD(T)、QCISD(T)都略遜于CCSD(T)一絲。特別值得一提的是CCSDT的誤差反倒比它的近似CCSD(T)更大。其原因在于，CCSD(T)雖然更形式上沒CCSDT精確，但是對T的微擾方式近似考慮導致的誤差，和忽略了Q的貢獻產生了一定抵消，因此反倒比嚴格考慮T時候的結果更好。實際上在一些文章對其它體系的測試中也有類似的情況。所以，不要盲目以為CCSDT一定比CCSD(T)更好，從而干出砸鍋賣鐵也非要攀CCSDT的這種事。還有很值得一說的是，也不要盲目認為(TQ)就比(T)更好。從圖中可見，(TQ)比(T)的誤差更大，而耗時則高得多得多，對當前體系耗時則和CCSDT接近。一些文章使用(TQ)，9成是認為精度比起用得臭大街的(T)更好，而且計算資源又有余裕，所以想搞個(TQ)來顯擺一下，但實際上精度往往還不如已經用爛的(T)。簡而言之，對于普通有機體系，要是嫌CCSD(T)還不夠高大上，計算資源富得流油，追求絕對的更高精度，能用的也就CCSDT(2)Q或者同檔的CCSDT(Q)了（都昂貴到炸，結合有實際意義的基組下能算得動的實際體系少之又少），而用CCSDT或(TQ)級別的方法可能還反倒費力不討好。雖然前面兩節的測試中CCSD[T]的相關能和參考值幾乎精確一致，但到了實際問題中，可見其誤差比CCSD(T)是明顯大一點的，而耗時和CCSD(T)相仿佛，所以沒有什么意義（不過有的文章如dx.doi.org/10.1021/ct3008777指出其計算弱相互作用好于CCSD(T)）。由圖可見MP5的精度和不考慮(T)的CC/QCI/BD半斤八兩，但耗時高得多，因此完全不值得用。CCSD的變體SCS-CCSD雖然眾所周知算弱相互作用精度極好，但是算當前問題，比CCSD并無優勢，異構化能的精度遠低于CCSD。MP4算的異構化能不錯，但勢壘算得比較糟糕。CI系列對相關能考慮效率低也直接體現在當前體系計算結果差上，很昂貴的CISDT誤差依然挺大。B3LYP對于當前體系意外地不錯（巧合所致），而M06-2X則令人比較失望，根據大量測試，其計算有機體系精度明顯好于B3LYP，但在當前體系卻栽了，也反映出DFT的可靠性終究還是和中高級別后HF有著顯著差距。圖中xfine是指在NWChem里使用比默認好得多的積分格點的情況，由此比較積分格點對結果的影響。雖然眾所周知明尼蘇達系列泛函對積分格點敏感，但對當前體系并未發現其結果對積分格點有顯著依賴性。MP2雖然通常對HF的結果有顯著改進，但對于當前體系表現令人極為失望，比發揮失常的M06-2X誤差還顯著更大，可見MP2算化學反應早已過時，沒有絲毫繼續使用的價值。SCS-MP2雖然算化學反應能一般認為不錯，原文里的測試發現有時能趕上QCISD的水準，但是對當前問題，表現得也很糟糕。CC2的誤差非常大，比CCSD差得遠，因此CC2也就以LR-CC2形式算激發能有點用，對于算化學反應毛用也沒有。CCD算反應能對此例恰好不錯，但勢壘遠不如CCSD，所以考慮T1算符的重要性極大。MP4(SDQ)是完整的MP4，即MP4(SDTQ)的近似，耗時降低了很多，但由圖可見誤差也有增加。最后，看基于DFT波函數做的后HF，即CCSD(T)-M062X、CCSD(T)-B3LYP和MP2-B3LYP，可見它們的誤差都很大，毫無用處，精度比標準的基于HF軌道的形式差得遠。

5 HCN異構化過程，cc-pVQZ基組

在實際中除了菜鳥外誰也不會用cc-pVDZ去試圖獲得定量準確的結果，即便對于DFT也太low。這里改用高精度計算常用的cc-pVQZ基組再做測試。由于基組較大，用CCSDT(2)Q對于筆者的計算條件已經很吃力了，而CCSDT、幾種(TQ)級別的方法又往往并不比CCSD(T)準，故這里以公認優秀的，而且前面測試中也表現優異的CCSD(T)的結果作為參考值，來衡量一批更便宜的方法的精度，結果如下（注意，實際上在CCSD(T)/cc-pVQZ這個級別時，誤差的很大程度來源已經是凍核近似了，試圖達到kJ/mol的精度應當考慮核電子相關，不過當前不考慮這個問題，并不顯著影響我們橫測的結論）。


由圖可見，大基組下B3LYP對這個體系表現依舊出色，但這只是僥幸而已，說明不了什么。SCS-CCSD在大基組下也依然沒比CCSD體現出優勢，反應能精度還不增反降。按照前人的測試，B2PLYP是肯定不如DSD-PBEP86-D3(BJ)的，不過當前體系前者也僥幸贏了后者。不過就算B2PLYP對當前體系打了雞血，還是沒法和CCSD一拼。是否考慮DFT-D3(BJ)，對B3LYP、B2PLYP的精度影響甚微，基本可以忽略不計，所以不涉及弱相互作用，加D3也沒什么壞處，但也別指望能帶來多大益處。Gaussian一派人搞的APFD在此體系表現差強人意。有趣的是，wB97XD和M062X對當前體系表現出來的精度很相近，都狀態不好。MP4沒什么用處，精度<=雙雜化，而耗時高得多得多，可以被埋汰了。SCS-MP2、SCS-MP3、MP2.5都是Grimme提出來的，然而從當前體系來看，都不算什么好idea，MP2.5比MP3更渣，SCS-MP3則是渣上加渣，居然反倒比MP3還爛。也隨便取了一個GGA泛函PBE納入測試，算的反應能很不錯，但必定只是巧合。PBE算勢壘的精度不理想，但也比比它昂貴得多的MP2強，再次說明MP2如今對于化學反應研究無用武之地，精度比其好的泛函一大把，更是被雙雜化完爆。隨手還測試了一下PM6，半經驗算出來的是焓變，雖說沒法和電子能量變化精確嚴格對比，但如果姑且當成一樣，PM6的反應能計算誤差是1.48kcal/mol，碰巧結果倒也不賴。不過其計算的反應能壘大約是CCSD(T)的兩倍，非常離譜！雖然一方面來自于焓變和電子能量變的差異，但更關鍵的還是在于一般的半經驗根本沒考慮勢壘的計算，誤差這么離譜倒也能理解。唯一一個專門算勢壘的半經驗方法是PM7-TS，讀者有興趣可以用MOPAC算算試試。

獲得CCSD(T)/大基組有一種常見的節約耗時的做法，就是通過以下公式估算：
CCSD(T)/大基組 = CCSD(T)/小基組 + (MP2/大基組-MP2/小基組)
這里也順帶測試了一下此做法的精度，還是以CCSD(T)/cc-pVQZ的結果作為參照：
小基組取cc-pVDZ：勢壘誤差0.13 kcal/mol，反應能誤差：0.11 kcal/mol
小基組取cc-pVTZ：勢壘誤差0.06 kcal/mol，反應能誤差：0.02 kcal/mol
可見這個近似公式精度很理想，實用性極強。即便小基組取cc-pVDZ，這種近似帶來的誤差也非常小，而降低的耗時以數量級論。只要能接受比雙雜化更高一些的計算花費，那么此做法無疑是最佳選擇（可惜在實際程序中沒法用于幾何優化。雖然原理上能實現，把梯度也這么組合便可）。之所以這個近似方法帶來的誤差這么小，原因在于：(1)MP2和CCSD(T)相關能本身差異就只有不到10% (2)所利用的是MP2和CCSD(T)下E(QZ)-E(DZ)的差值，差異進一步縮小 (3)化學上研究的是體系間相對能量，求差時誤差進一步抵消。

6 總結

再次強調一下，本文的測試非常粗糙，切勿根據本文的數據以偏概全說事。衡量方法的準確度，最少也得用包含5個體系的測試集進行統計分析，而如今benchmark研究都是用動輒幾十，乃至上百，甚至上千的測試體系，否則很可能某些本來精度不佳的方法由于誤差抵消的巧合恰好顯得精度不錯，比如上面測試中B3LYP精度顯得比M06-2X還好就有違一般情況。不過，此文的測試在一定程度上還是能以小見大的，結合前人的大量測試產生出的普遍共識，本文的測試反映出以下結論（大多其實都是老生常談的，但也有少數是許多人沒有注意的）：

(1)如果研究的是微小的體系，就幾個原子，不滿足于CCSD(T)，還想要更高精度（或者就是單純想顯擺），那么CCSDT(2)Q方法（或者同檔次的CCSDT(Q)，MRCC程序支持）是最佳選擇，精度和Full CI幾乎沒差別，但耗時奇高。

(2)使用耦合簇時，(T)是非常重要的，能顯著提升耦合簇的精度，對于高精度計算是必須的。凡是能用得動CCSD時，一定要努力沖一下(T)。如果帶著(T)死活就是算不動，那根本別用CCSD，還不如用前述的近似得到CCSD(T)/大基組的方法，或者嘗試ORCA的DLPNO-CCSD(T)（通過引入數值上的近似來降低CCSD(T)耗時，對于大體系效果明顯）。

(3)CCSD(T)的精度很理想。同檔次的CCSD(2)T、CCSD[T]都遜色于它，這也是為什么CCSD(T)遠比其它同檔次對CCSDT的微擾近似常見得多得多的原因之一。

(4)QCI、BD比起耦合簇沒有優勢，BD耗時還更高。doi: 10.1002/wcms.1131更是強調QCI早該被遺忘了，不應該再被使用。除非你的目的是重現其它人用QCI得到的結果，否則打死也別用QCI。

(5)不要以為(TQ)比(T)的精度一定更好，也不要以為CCSDT比CCSD(T)的精度一定更好，至少此文測試給出了反例。最常用的Gaussian支持的最高精度的單參考方法從形式上來說是QCISD(TQ)和BD(TQ)，但試圖用這倆來獲得比CCSD(T)更高的精度，往往會自取其辱。

(6)不要以為SCS的做法總能改進結果，本文給出了鮮明的反例，別被原文和一些測試文章的光鮮數據給忽悠了，實際沒有想象的那么好，用了SCS變體精度極可能不增反降。MP2.5對于一般化學反應問題也得不到什么好效果，性價比很低。用SCS-MP2遠遠不如用雙雜化泛函。

(7)MP2如今完全不中用，廉頗老矣。原先的MP2的市場已經完全是DFT的天下。MP4、MP5的性價比都很差，跟耦合簇沒法比，不要用。MP3也毫無用武之地。

(8)CC2比CCSD精度差得遠，當用不動CCSD時，不要試圖用CC2來節約計算量。

(9)以DFT波函數做參考態結果很糟糕。雖然有些文章這么做卻得到了不錯結果，但可能是其用的程序實現形式和NWChem的不同所致。

(10)CI系列爛得很，完全無價值，比微擾還沒用。

(11)用半經驗方法算勢壘風險極大，很可能得到完全離譜的結果，其算勢壘的誤差普遍遠高于算反應能的誤差。

(12)借助MP2近似計算CCSD(T)/大基組下的能量是極為有用的方法，一定要善加利用。

(13)僅從對相關能的考慮程度來衡量方法的精度很片面，并不能可靠展現方法對實際化學問題的計算精度。但是，對實際問題總是表現良好的方法，對相關能的考慮都很充分。

另外，還要強調的是本文考察的體系并沒有很強的多參考特征，涉及多參考體系時往往會有不同的景象（比如CCSDT可能會展現比CCSD(T)明顯的優勢）。上面的總結對于一般有機體系的研究是普遍適用的。

附錄：一些輸入文件示例

考慮到一些讀者可能不會寫本文提及的一些小眾向計算方法的輸入文件，但可能又想用它們，這里就簡單說明一下。

DSD-PBEP86-D3、MP2.5、SCS-MP2、SCS-MP3的計算方法看此文：《Gaussian中非內置的理論方法和泛函的用法》（http://www.shanxitv.org/344）。BD(TQ)、QCISD(TQ)、MP4(SDQ)、MP5就直接在Gaussian輸入文件里寫這關鍵詞即可。

NWChem做CCSDT(2)_Q/cc-pVTZ計算的輸入文件為
start  //計算時完全重新開始，不基于之前的波函數(若有的話）
memory total 7000 global 6000 mb   //內存設置，詳見手冊
GEOMETRY
symmetry c2v   //點群
 C                 -0.00000000    0.00000000   -0.74384793
 N                  0.00000000   -0.00000000    0.43281671
 H                  0.00000000   -0.00000000    1.43337065
end
basis spherical   //設置基組。用球諧型基函數。默認是笛卡爾型
* library cc-pVTZ   //所有原子用的基組
end
TCE   //設置TCE模塊
freeze atomic   //凍核。默認不凍核
CCSDT(2)_Q
END
TASK TCE ENERGY    //用TCE模塊算單點
這里可以把CCSDT(2)_Q改為LR-CCSD(TQ)-1、CCSD(2)_T、CR-CCSD(T)、CCSDT、CCSDTQ、CISD、CISDTQ、CC2、CCSD(T)、MP2、MP3、MP4等來做相應計算。CCSD[T]結果會在CCSD(T)任務中一并輸出。雖然此體系是C∞v對稱性，但NWChem的TCE模塊只支持阿貝爾群，所以只能用阿貝群里最高的C2v。

如果想基于DFT波函數做后HF計算，比如基于B3LYP做CCSD(T)，用：
start
memory total 7000 global 6000 mb
GEOMETRY
symmetry c2v
 C  0. 0. 0.
 O  0. 0. 1.128323
end
basis
* library cc-pVDZ
end
DFT
XC B3LYP
END
TCE
DFT
freeze atomic
CCSD(T)
END
TASK TCE ENERGY

若要計算SCS-CCSD，需要用NWChem里的CCSD模塊來做，不能用TCE模塊。寫法：
...其它設定
CCSD
freeze atomic
END
TASK CCSD ENERGY
輸出信息里就會看到SCS-CCSD能量。

如果不會編譯NWChem，看《NWChem 6.5、6.6編譯方法》（http://www.shanxitv.org/270）。