Gaussian中非內置的理論方法和泛函的用法
Gaussian中非內置的理論方法和泛函的用法
文/Sobereva@北京科音
First release: 2016-Sep-4 Last update: 2022-Feb-14
本文把Gaussian中沒有內置關鍵詞，但能利用IOp或其它手段等效實現的一些理論方法的用法進行介紹。
首先要提醒一點，對于諸如opt freq這樣多步任務，輸入文件里設的IOp只對第一步生效，不會自動傳遞給后面的任務。所以通過IOp使用非內置方法的時候不要為省事把opt和freq放在一起用。
1 微擾相關的方法
1.1 SCS-MP2
這是Grimme在J. Chem. Phys., 118, 9095 (2003)提出的對MP2的改進。MP2相關能可以寫為電子對的加和，SCS-MP2把自旋相同電子對的貢獻乘上1/3以削弱之，而把自旋相反的貢獻乘上6/5以增強之。SCS-MP2算反應能比MP2有了很大提高，達到QCISD級別甚至個別時達到QCISD(T)級別。幾何結構也比MP2略有改進。但勢壘計算精度比MP2反倒略低一點，弱相互作用計算精度惡化不少。其它方面和MP2差不多。整體性能明顯不如目前最好的雙雜化泛函，所以如今用處也不太大。SCS-MP2的使用方法是在MP2計算時額外加上IOp(3/125=0333312000)。3/125=MMMMMNNNNN代表將自旋平行項與自旋反平行項的貢獻以MMMMM/10000比NNNNN/10000的比例組合。
MP2計算時從輸出中看到
alpha-alpha T2 = 0.1336200196D-01 E2= -0.4731575741D-01
alpha-beta T2 = 0.8298753242D-01 E2= -0.3004708246D+00
beta-beta T2 = 0.1336200196D-01 E2= -0.4731575741D-01
SCS-MP2計算時為
alpha-alpha T2 = 0.1336200196D-01 E2= -0.1577034194D-01
alpha-beta T2 = 0.8298753242D-01 E2= -0.3605649895D+00
beta-beta T2 = 0.1336200196D-01 E2= -0.1577034194D-01
可見alpha-alpha、beta-beta的能量貢獻都被乘上了1/3，而alpha-beta乘上了6/5。
1.2 SOS-MP2
于J. Chem. Phys., 121, 9793(2004)中提出，忽略了自旋平行貢獻，反平行部分乘以1.3，結果比MP2爛，特別是弱相互作用。不過號稱通過Laplace變換方法可以令標度降到O(N^4)從而使花費比MP2的O(N^5)低。不過再便宜也沒普通DFT便宜，MP2精度尚且不及選得最合適的DFT泛函，所以SOS-MP2這個破玩意兒并沒什么卵用。大家知道有這么個垃圾就夠了。用法是MP2 IOp(3/125=0000013000)。由于Gaussian并未對之進行優化，所以耗時不會比MP2低。
1.3 SCSN-MP2
在J. Chem. Theory Comput., 3, 80 (2007)中提出，忽略自旋反平行部分，自旋平行部分的系數為1.76，擬合自核酸堿基對相互作用能。弱相互作用計算精度挺不錯，誤差比MP2低兩三倍，有使用價值。另外還有Mol. Phys., 105, 1073 (2007)中提出的SCS(MI)-MP2，MI代表molecular interaction，它與SCSN-MP2的提出目的和精度都很類似，并對不同基組擬合了不同的系數，這里就不提了。從SCSN-MP2和SCS(MI)-MP2的系數都能看出，對于計算弱相互作用的目的，應該削減MP2自旋反平行而增大自旋平行的貢獻，這和一般目的的SCS-MP2恰好相反。用法是MP2 IOp(3/125=1760000000)。
1.4 MP2.5
是在ChemPhysChem, 10, 282-289 (2009)中提出的，將MP2總能量加上乘以0.5的MP3三階微擾能。計算量和MP3一樣。系數0.5來自于分析計算精度、基組依賴性和理論意義。各方面性能比MP2強，特別是計算弱相互作用準確度很高，號稱在中等基組下（不加彌散亦可）就能接近CCSD(T)/CBS。在Gaussian中用MP3計算完后，會看到比如
E2 = -0.3951023394D+00 EUMP2 = -0.11430702788798D+03
Keep R2 and R3 ints in memory in symmetry-blocked form, NReq=9641770.
DD1Dir will call FoFMem 1 times, MxPair= 42
NAB= 21 NAA= 0 NBB= 0.
E3= -0.48480932D-02 EUMP3= -0.11431187598D+03
MP2.5能量就是-114.30702788798-0.0048480932*0.5=-114.3094519
后來在Phys. Chem. Chem. Phys., 13, 21121 (2011)中還提出了MP2.X，X是給三階微擾能乘上的系數。MP2.5用在小基組上結果不如在大基組好，為了使得較小基組下就能得到不錯的弱相互作用計算結果，MP2.X對從小到大的基組都通過S66測試集擬合了MP3校正能的權重系數，這使得不同基組下（乃至低至6-31G*）得到的弱相互作用能精度都相仿佛，和MP2.5/aug-cc-pVTZ下差不多。雖看似很好，但對更多的體系的可靠性還有待廣泛驗證。
由于MP2.5、MP2.X能量只能在MP3單點任務結束后根據輸出信息手動計算，沒法像SCS-MP2那樣直接用IOp，所以這兩種方法只能算能量，沒法用來優化、計算頻率等，除非自己去改Gaussian代碼實現。
1.5 SCS-MP3
Grimme在J. Comput. Chem., 24, 1529 (2003)中提出的，是SCS-MP2能量加上乘以0.25的三階微擾校正能。熱化學性能比SCS-MP2好很多，往往接近QCISD(T)。但弱相互作用能精度不算理想。用法是MP3 IOp(3/125=0333312000)，然后把輸出文件中EUMP2能量加上乘以了0.25的E3。
順帶一提，SCS的思路也被用在CI和耦合簇方法上。SCS-CCSD的正反自旋參數來自于擬合幾十個反應能。算反應能，特別是算弱相互作用能比CCSD好很多。Phys. Chem. Chem. Phys., 12, 9611 (2010)中提出的SCS-MI-CCSD則令參數向S22弱相互作用測試集擬合，算弱相互作用極為接近金標準CCSD(T)，比MP2.5、SCS-CCSD又明顯更好。但可惜Gaussian09中沒法實現這些SCS的耦合簇方法，因為程序也不把耦合簇的自旋平行和反平行的貢獻獨立輸出。
2 非內置的普通DFT泛函的用法
2.1 規則
在Gaussian中可以通過IOp在一定程度上自定義泛函。簡單來說，若定義IOp(3/76=ab)、IOp(3/77=cd)、IOp(3/78=ef)，這里abcdef都已經除以了10000，則所用泛函的實際形式是（下文說的GGA也包含meta-GGA）E_XC = a*[ d*E_X_LDA + c*ΔE_X_GGA ] + b*E_X_HF + f*E_C_LDA + e*ΔE_C_GGA
每一項的含義如下
E_X_LDA：LDA交換項
ΔE_X_GGA：GGA交換項對LDA交換項的梯度修正部分
E_X_HF：HF交換項。前頭的系數b就是常說的HF雜化成份
E_C_LDA：LDA相關項
ΔE_C_GGA：GGA相關項對LDA相關項的梯度修正部分
無論對于交換還是相關部分，都是E_GGA=ΔE_GGA+E_LDA。比如E_X_B88=E_X_LDA+ΔE_X_B88，E_C_LYP=E_C_LDA+ΔE_C_LYP。
對于DFT，a總為1（其實這個參數純屬多余，即便有用也可以直接融合到d、c里嘛）。對于純泛函，滿足d=c=1、f=e=1、b=0；對于單參數雜化泛函，比如PBE0、BHandHLYP，滿足d=c、d+b=1、f=e。對于B3LYP這樣的三參數雜化泛函才可能d!=c、f!=e。而對于HF，顯然b=1，其它皆為0。
用#P在計算時會有諸如這樣的輸出：
IExCor= 402 DFT=T Ex=B+HF Corr=LYP ExCW=0 ScaHFX= 0.200000
ScaDFX= 0.800000 0.720000 1.000000 0.810000
其中ScaHFX就是b，ScaDFX后面的值是d c f e。Ex就是所用的交換泛函，當前是B88和HF雜化。Corr是所用的相關泛函，當前是LYP。
2.2 B3LYP及變體
老不死（非貶義）的B3LYP定義為E_XC_B3LYP = (1-a0)*E_X_LDA + a0*E_X_HF + aX*ΔE_X_B88 + E_C_VWN + aC*ΔE_C_LYP
其中a0=0.2、aX=0.72、aC=0.81。
顯然，換算成上一節的系數表達方式，就是a=1.0，b=0.2，c=0.72，d=0.8，e=0.81，f=1.0。相應地，計算的時候關鍵詞就寫成BLYP IOp(3/76=1000002000,3/77=0720008000,3/78=0810010000)。
B3LYP在提出的時候沒有明確說清楚LYP用的VWN到底是VWN3還是VWN5，二者結果有一點差異。Gaussian里默認用的是前者，比后者更好一點點，見Chem. Phys. Lett., 268, 345 (1997)的對比。但有的程序默認用的則是VWN5，如GAMESS-US。如果想用VWN5的B3LYP的話，就可以寫BV5LYP IOp(3/76=1000002000,3/77=0720008000,3/78=0810010000)，這里V5LYP相關泛函是指把VWN5作為LYP里的LDA部分。
在J. Chem. Phys., 117, 4729 (2002)中提出了B3LYP*，它把B3LYP的HF成份從20%降低到15%，這大大解決了B3LYP高估過渡金屬配合物高自旋態穩定性的問題（即低估了高、低自旋態間的能量差）。用B3LYP*就寫BLYP IOp(3/76=1000001500,3/77=0720008500,3/78=0810010000)，即把X_HF的系數b降到0.15，相應地把X_LDA的系數d提升到0.85，使二者加和仍保持為1。
HF成份對結果有很多系統性的影響，比如HF成份越大TDDFT算的激發能越高。因此可以容易地調節b并相應地修改d，使得算出的光譜和實際盡量接近（雖然灑家鄙夷這種刻意往實驗結果上湊的勾當）。
2.3 PBE0及變體
PBE0在Gaussian里寫作PBE1PBE，是把PBE泛函中引入25% HF成份，即a=1.0，b=0.25，c=d=0.75，e=f=1.0。相當于PBEPBE IOp(3/76=1000002500,3/77=0750007500)。由于IOp(3/78=1000010000)是默認的，所以可以不用寫。J. Chem. Phys., 138, 021104 (2013)定義了一個無聊的泛函PBE0-1/3，就是把HF成份增加到1/3，故對應于PBEPBE IOp(3/76=1000003333) IOp(3/77=0666706667)。
DFT-D3原文里順帶定義了PBE38，是把HF成份改為3/8=37.5%，測試表明算含頻極化率是最好的，對應于PBEPBE IOp(3/76=1000003750) IOp(3/77=0625006250)。
2.4 明尼蘇達系列老泛函
Truhlar組之前搞過一大堆亂七八糟的泛函，什么MPW1K、PBE1KCIS，也沒多少人用，就他們自己玩得happy。從M05起我看才算真正步入正軌。那些爛七八糟的老明尼蘇阿達系列泛函我強烈不推薦使用，早被更好的泛函完爆了，初學者也甭老效仿文獻試圖去用。如果非要用，去看他們網頁http://comp.chem.umn.edu/info/DFT.htm，能通過Gaussian實現的泛函都給了用法。而在Gaussian中即便通過IOp也使不了的，死活非要用就必須找他們要修改版的Gaussian了。2.5 對長程校正泛函調節ω
Gaussian中對純泛函可以用LC-前綴做長程校正，比如LC-BLYP。此方法將1/r12算符劃分為短程和長程部分，短程不用HF交換項，而長程用100% HF交換項。此做法有一個參數ω很重要，此值越大長程部分涵蓋范圍就越廣，會明顯影響結果。LC-加到各種純泛函上默認是ω=0.47,ωB97X默認是ω=0.3，LC-wPBE默認是0.4。如果要調節這個參數，就把3/107和3/108都設為MMMMM00000，代表ω值為MMMMM/10000。比如要把LC-wPBE的ω設0.3就寫LC-wPBE IOp(3/107=0300000000,3/108=0300000000)。
對于范圍分離泛函還可以通過調節α、β參數來讓短程極限和長程極限的HF成份對應于特定的數值，筆者專門在這篇博文里進行了介紹：《在Gaussian中自定義范圍分離泛函的方法》（http://www.shanxitv.org/550）。
2.6 QTP17泛函
在《正確地認識分子的能隙(gap)、HOMO和LUMO》（http://www.shanxitv.org/543）中筆者介紹了QTP系列泛函。其中QTP17泛函結合aug-cc-pVTZ基組計算出的價層和內層軌道能量都能較好滿足Koopmans定理，因此軌道能量的負值直接就可以近似當做光電子譜測定的外層、深層、內核電離能了。而且還可以直接結合Multiwfn繪制光電子譜而不需要對軌道能量做任何校正，參見《使用Multiwfn繪制光電子譜》（http://www.shanxitv.org/478）。
在Gaussian中QTP17的用法是寫# SLYP/aug-cc-pVTZ IOp(3/76=1000006200,3/77=0000003800,3/78=0800010000)。注意由于QTP17泛函參數是在aug-cc-pVTZ基組下訓練的，若改成其它基組精度會下降。
3 非內置的雙雜化泛函的用法
3.1 老式雙雜化泛函
比較老的雙雜化泛函，如B2PLYP、mPW2PLYP、B2GP-PLYP沒有考慮色散校正也沒有對摻入的二階微擾能E2以SCS方式考慮。這類老雙雜化泛函的交換相關能通式為E_XC = (1-c1)*E_X_GGA + c1*E_X_HF + (1-c2)*E_C_GGA + c2*E2
最老、最知名，但也是最差的雙雜化泛函B2PLYP是J. Chem. Phys., 124, 034108 (2006)中提出的，也是Gaussian09中內置的，定義為
E_XC=0.47*E_X_B88+0.53*E_X_HF+0.73*E_C_LYP+0.27*E2
計算時會看到相應的泛函定義信息，只要看懂了2.1節就自然能理解：
IExCor= 419 DFT=T Ex+Corr=B2PLYP ExCW=0 ScaHFX= 0.530000
ScaDFX= 0.470000 0.470000 0.730000 0.730000 ScalE2= 0.270000 0.270000
這里比2.1節看到的多了ScalE2，后面兩個值是給E2的自旋平行和反平行部分分別乘的系數，這里都是0.27，所以說沒有以SCS方式考慮E2。
B2GP-PLYP是在B2PLYP之后提出的，原文見J. Phys. Chem. A, 112, 12868 (2008)，精度比B2PLYP好不少，但比起后來以SCS方式考慮E2的雙雜化泛函還是差多了，其系數c1=0.65，c2=0.36，純泛函部分還是BLYP，可明確寫為
E_XC=0.35*E_X_B88+0.65*E_X_HF+0.64*E_C_LYP+0.36*E2
Gaussian09支持的很可惜只有最老最爛的B2PLYP和后來提出的單精度與之半斤八兩的mPW2PLYP。B2GP-PLYP在Gaussian09中沒有內置，好在若想用的話可以寫
B2PLYP/cc-pVTZ IOp(3/125=0360003600,3/76=1000006500,3/77=0350003500,3/78=0640006400)
這里還寫B2PLYP是為了讓Gaussian做雙雜化計算，但是參數已經被我們改了。3/125=MMMMMNNNNN代表給E2的自旋平行和反平行部分分別乘上MMMMM/10000和NNNNN/10000，其它參數前面已經講過了。
注：也可以把3/76=1000006500,3/77=0350003500改為3/76=0350006500,3/77=1000010000，根據2.1節的表達式就知道這兩種寫法是等價的。
3.2 較新的雙雜化泛函
J. Phys. Chem. C, 114, 20801 (2010)中提出了DSD-BLYP，它考慮了色散校正，并且以SCS方式考慮了E2的貢獻，精度比B2GP-PLYP有一定提升。從此之后各種通式類似的雙雜化泛函冒出來不少，比如DSD-PBEP86、DSD-PBEhB95等，此二者都比DSD-BLYP精度又有進一步提升。這些泛函通式寫為這樣：
E_XC = (1-cX)*E_X_GGA + cX*E_X_HF + cC*E_C_GGA + cO*E2_OS + cS*E2_SS + E_D
這里E2_OS和E2_SS分別是E2中自旋相反和相同部分的貢獻，E_D是DFT-D3色散校正能，一般都用BJ阻尼。
J. Comput. Chem., 34, 2327 (2013)是一篇不錯的文章，在其中把各種常見純泛函組合代入到上面的公式來擬合了系數，包括色散校正參數。經過比較，發現DSD-PBEP86-D3(BJ)精度名列前茅，而且難得的是它能直接在許多主流量化程序里不需要修改源代碼就能使用。在此文的補充材料里可以找到各種純泛函對應的雙雜化泛函的參數，在“DSD-DFT-D3BJ”表格中可見對于DSD-PBEP86-D3(BJ)，cX=0.69，cC=0.44，cO=0.52，cS=0.22，DFT-D3(BJ)的參數為s6=0.48，a2=5.6。
DSD-PBEP86-D3(BJ)從Gaussian 16開始已經內置了，寫DSDPBEP86關鍵詞即可。這里說一下怎么在Gaussian 09中通過特殊方式等價地使用。DFT-D3(BJ)的參數在Gaussian09輸入文件里沒法直接設，自定義的話必須得設定環境變量。對于Linux版，bash環境，在計算之前先把以下內容輸入到終端中
export GAUSS_DFTD3_S6=480000
export GAUSS_DFTD3_SR6=0
export GAUSS_DFTD3_S8=0
export GAUSS_DFTD3_ABJ1=0
export GAUSS_DFTD3_ABJ2=5600000
之后輸入文件里的關鍵詞就寫
B2PLYP/cc-pVTZ IOp(3/125=0220005200,3/76=1000006900,3/77=0310003100,3/78=0440004400,3/74=1004) empiricaldispersion=GD3BJ
這里3/74=1004代表把X_GGA設為PBE，C_GGA設為P86。
注1：也可以像此泛函提出者的網頁http://www.compchem.me/dsd-pbep86-functional中那樣寫為IOp(3/125=0220005200,3/78=0440004400,3/76=0310006900,3/74=1004)，是等價的。
注2：雖然JCC文中測試表明DSD-PBEhB95-D3(BJ)綜合性能最好，但在G09中通過上述方法使用會報錯，需要對代碼進行特殊修改。
類似地，還可以根據補充材料里的信息使用各種其它的雙雜化泛函，但大多不如DSD-PBEP86-D3(BJ)。另外，如果不想用色散校正的版本（比如只有G09 D.01以前版本或者嫌設定環境變量麻煩），即只用DSD-PBEP86，參數是為cX=0.72，cC=0.44，cO=0.51，cS=0.36，故應當寫
B2PLYP/cc-pVTZ IOp(3/125=0360005100,3/76=1000007200,3/77=0280002800,3/78=0440004400,3/74=1004)
如文中的數據所示，不用DFT-D3(BJ)的話整體誤差會有所增加，但個別項目可能反倒誤差會降低（不排除巧合可能）。
2020-Sep-19補充：在J. Phys. Chem. A, 123, 5129 (2019)中DSD-PBEP86-D3(BJ)的作者又提出了其修改版revDSD-PBEP86-D3(BJ)，精度有了明顯提升，在所有雙雜化泛函中精度已經是名列前茅了，精度和目前（2020年）最佳的雙雜化泛函ωB97M(2)已經很接近了。在Gaussian 16中可以通過以下關鍵詞使用
# DSDPBEP86/基組 IOp(3/125=0079905785,3/78=0429604296,3/76=0310006900,3/74=1004) em=gd3bj IOp(3/174=0437700,3/175=-1,3/176=0,3/177=-1,3/178=5500000)
3.3 其它雙雜化泛函
18碳環是非常獨特的體系，筆者做過大量研究，匯總見http://www.shanxitv.org/carbon_ring.html。在Phys. Chem. Chem. Phys. (2022) DOI: 10.1039/d1cp04996h中作者表明范圍分離雙雜化泛函RSX-QIDH（也叫RSX-PBE-QIDH）可以比較好地計算碳環鍵長均等的過渡態結構和非均等的極小點結構的能量差。此泛函在Gaussian 16中的用法為：
PBEQIDH IOP(3/74=1001009,3/78=0666606666,3/107=0270000000,
3/108=0270000000,3/119=0310000000,3/120=0310000000,
3/125=0333403334,3/130=06900,3/131=06900)