談談彌散函數和“月份”基組
談談彌散函數和“月份”基組
文/Sobereva @北京科音
First release: 2012-Jan-13 Last update: 2021-Sep-28
如果你對基組了解甚少，強烈建議閱讀本文前先閱讀《談談量子化學中基組的選擇》（http://sobereva.com/336）。
1 何時需要彌散函數？
量子化學計算中的彌散函數是指指數很小的基函數，有很廣的空間分布范圍。加彌散函數的必要性，筆者根據大量理論計算文章和實踐經驗總結如下：
? 不加彌散結果一定沒法用，必須加彌散：
計算偶極矩、多極矩
計算極化率、超極化率
計算里德堡激發態
計算陰離子體系能量、電子親和能
? 加彌散很重要，強烈建議加彌散：
計算弱相互作用能
計算Raman、ROA強度
? 加彌散對改進結果有益，計算條件允許可以加彌散：
計算反應勢壘
優化陰離子體系或弱相互作用體系的結構、對其進行振動分析
如果體系只有部分原子帶十分顯著的負電荷或者其電子十分容易被極化，或者弱相互作用只涉及到體系中的一部分，限于計算能力沒法給所有原子加彌散的話，那么至少要給這些區域的原子加上彌散函數。
當彌散函數必不可少時，為了能加上彌散函數寧可降低zeta數，比如計算偶極矩、極化率時寧可從3-zeta降到2-zeta也要把彌散函數加上。但對于彌散函數不起到關鍵性作用但是加它有益處的時候，那么至少先把基組用到不低于3-zeta的檔次再考慮加彌散，這樣才最劃算，因為此時增加基組價層分裂數帶來的益處明顯比加彌散函數大得多。
對于中性體系且其中不涉及顯著的弱相互作用、沒有局部帶顯著負電荷的時候，計算其上述未提到的特征，比如原子化能、價層電子激發、鍵能、電離能、幾何優化、NMR等等，不要加彌散，否則對結果精度沒什么改進，還徒增巨額計算時間、導致很多問題。記住亂用彌散函數只會自取其辱！
PS：筆者時不時看到居然有人算陽離子的普通問題（如能量、幾何優化）時還加彌散函數，簡直難以理喻，陽離子又不像陰離子那樣有被束縛得較弱的電子，加什么彌散函數！？審稿人一看就知道這作者缺乏計算化學常識，好感度驟降！（后來我發現某些小白給陽離子加彌散的一個可能原因：有一次有人在我的群里貼了張講基組的中文幻燈片里的圖，背景是深藍的，里面說對離子體系應該加彌散，居然也不區分是陽離子還是陰離子，這ppt真是超級誤人子弟！似乎這ppt在中文的低水平的計算化學討論場所流傳還挺廣，導致坑害了一堆初學者。所以千萬別隨便信網上來路不明的中文資料，這在《談談學量子化學如何正確地入門》http://www.shanxitv.org/355里專門說了）
由于給氫亂用彌散函數的情況特別普遍，很多人還屢教不改，所以這里專門強調一下。給氫加彌散函數絕大多數情況沒有任何意義。因為氫本身只有一個電子，電負性又小，在分子環境中一般都會被吸走不少電子，給它加彌散函數顯然對結果不會有可察覺的改進，而有機體系里氫的數目通常很多，因此給氫加彌散還令計算量增加不少。此外，凡是給氫加彌散函數的時候肯定也會給重原子加，重原子的彌散函數會很大程度延展到氫那里去，也一定程度上表現了氫的彌散函數起到的作用，因此更沒必要給氫加彌散函數。算普通問題給氫加彌散函數是多余的這點早就被大量文章的測試所證實，也被內行研究者所廣泛認同。所以，需要對體系加彌散函數時，強烈不建議用諸如6-311++G**，而建議用6-311+G**。對于Dunning系列基組要加彌散的情況，也建議用后文說的月份基組而不是aug-cc-pVnZ。即便研究氫鍵鍵能也不需要給氫加彌散函數，測試表明這對結果的影響微乎其微。但對于計算(超)極化率、氫帶明顯負電（如LiH）等情況，給氫加彌散函數確實有重要意義，此時需要給氫加彌散。
2 彌散函數造成的問題
加了彌散函數會使基函數的完備性增加，從原理上講似乎總應當對結果有益無害，但在實際計算中卻會導致如下問題：- 計算量猛增。這是眾所周知的，尤其是cc-pVnZ系列，往往加了aug-后就算不動了，第四節將討論怎么解決這個困難。
- 加入了彌散函數后SCF經常比未加時難收斂得多。
- 彌散函數的化學意義很差，與原子軌道缺乏對應關系，會對希爾伯特空間下的波函數分析方法造成嚴重不良影響。比如Mulliken電荷會變得極爛，Mayer鍵級會頗不可靠。原因不難理解，比如A原子一大堆彌散函數延伸到B原子空間內，因此B附近的電子分布有一部分會被這些彌散函數所描述，那么Mulliken布居分析就會把很多本應屬于B的電子劃歸到A原子上，導致A的電荷過負而B的過正。
- 虛軌道（即非占據軌道）的化學意義變得更含糊。尤其是彌散函數下的Hartree-Fock計算，虛軌道空間分布范圍往往會變得特別廣，導致前線軌道理論分析完全不再適用。
- 導致出現基函數線性依賴問題引起數值問題。不過，一般量化程序中都會自動檢驗基函數重疊矩陣的本征值來適當砍掉些基函數以解決這個問題。
- 可能因數值計算精度問題導致體系結構和波函數對稱性下降。比如你優化個有對稱性的體系，一開始的結構程序是能判斷出實際點群的，但優化之后，程序有時只能判斷出更低階點群了。如果你當前用了彌散函數，那么8成是彌散函數導致的，去掉彌散函數往往就能維持住對稱性。
- 如果原始基組完備性不高，卻增加了過多彌散函數，則本應該由價層基函數表現的效應會轉而被彌散函數所表現，在一些問題的研究中可能引起不合理的結果，這其實屬于分子內BSSE范疇。一個例子是JACS,128,9342發現后HF結合某些彌散版本的Pople基組（如6-31++G**）算出來的苯的穩定結構竟然是彎的，或者說本該穩定的平面結構卻有虛頻。這是因為這些Pople基組中沒有后HF計算較依賴的更高角動量基函數（尤其是f），延伸過去的較為彌散的s和p基函數為了能充分等效展現出碳的更高角動量基函數的效應而引起了結構的彎曲。
3 常見的含有彌散函數的基組
先說一下彌散函數的一般特征。每種角動量的彌散函數的指數小于基組中其它同等角動量函數的最小指數的數倍。各種基組中的彌散函數數目、所涉及的角動量不同。彌散函數一般是非收縮的。由于彌散函數對許多情況很重要，所以主流的基組大多數都有帶彌散函數的版本。帶彌散函數的版本有些是原基組的作者搞出來的，也有的是其它研究者提出來的。下面列舉一些常見的：
1 Pople系列基組：只給重原子加上一層sp（即指數相同的一層s和一層p）彌散函數就在G前頭添上一個加號，如6-31+G*；若同時還給氫、氦原子加上一層s彌散函數就在G前頭再添上一個加號，如6-311++G(2df,2p)。Pople系列基組有很多，但是彌散函數的指數都是共用的，并未單獨優化。（值得一提的是，6-311G*的耗時是要低于6-31+G*的，對于彌散函數不起到關鍵作用時，前者的結果還比后者好。無數人用6-31+G*的時候殊不知，對于他們所研究的問題，用6-311G*劃算得多得多！）
2 Dunning相關一致性基組（cc-pVnZ系列）：加上彌散函數的版本是aug-cc-pVnZ系列（aug=augmented），是給相應的cc-pVnZ基組的每種角動量函數上都添加一層同等角動量的彌散函數。例如cc-pVTZ對于C是4s,3p,2d,1f，因此aug-版本會增加一層s、一層p、一層d和一層f彌散函數。而cc-pVDZ對于氫是2s,1p，因此aug-版本會給它增加一層s和一層p彌散函數。
與aug-cc-pVnZ完全相同的彌散函數加到考慮描述核相關的cc-pCVnZ和cc-pwCVnZ基組上成為aug-cc-pCVnZ和aug-cc-pwCVnZ；加到適合DKH計算的cc-pVnZ-DK上成為aug-cc-pVnZ-DK；加到cc-pV(n+d)Z系列（這套基組是在cc-pVnZ基礎上加入一層比較緊的d函數以改善外推時的收斂性）成為aug-cc-pV(n+d)Z，其中僅d彌散函數的指數經過了重新優化。相關一致性贗勢基組cc-pVnZ-PP也有彌散函數的版本aug-cc-pVnZ-PP，彌散函數加入的類型和數量同aug-cc-pVnZ系列，但指數都重新優化了。另外，還可以對cc-pVnZ系列每個角動量加上多層彌散函數，d-aug-cc-pVnZ、t-aug-cc-pVnZ就是給每個角動量上分別增加兩層和三層彌散函數，極為昂貴，一般用于精確計算里德堡激發態、超極化率等目的。
3 Ahlrichs的def2-系列基組：目前這類基組并沒有官方的帶彌散函數的版本，這不能不說是個遺憾，但是有其它許多方式加彌散，見《給ahlrichs的def2系列基組加彌散的方法》（http://sobereva.com/340）。
4 Jensen的極化一致性基組pc-n：在JCP,117,9234中Jensen提出了給他的pc-n系列基組增加彌散函數的方法，彌散函數加到多高角動量可以根據需要自行選擇。結果表明加了s和p彌散能讓DFT計算的電子親合勢精度大有提高，而進一步改善響應屬性的計算結果還同時需要加更高角動量的彌散函數。
5 Lanl贗勢基組：在JPCA,105,8111中，作者給主族Lanl2DZ添加了一層d極化和一層p彌散函數，命名為LANL2DZdp，在計算電子親和能、振動頻率、鍵長方面都比Lanl2DZ大有提高。LANL2TZ+和LANL08+分別是給LANL2TZ和LANL08對于第一周期過渡金屬增加了一層彌散，這是考慮到d殼層充滿的第一周期過渡金屬有時容易被極化。
除了LANL2TZ+和LANL08+等基組的彌散函數的指數用的是even-tempered方式推導出的以外，上面提到的彌散函數的指數大多都是來自于不同方式最小化陰離子計算的能量（注：原基組的收縮系數和指數保持不變，彌散函數只是單純地累加到原基組）。然而，計算能量好的彌散函數基組用于計算其它屬性，尤其是強烈依賴彌散函數的（超）極化率未必好，或者說性價比不高。有不少人為了能讓（超）極化率等響應屬性在較低的計算成本下有滿意的計算結果，令彌散函數，乃至整個基組都直接來自于優化響應屬性的計算，這里也舉幾個例子：
Sadlej POL：也叫Sadlej pVTZ基組，從1988年陸續提出，參數是優化計算極化率得到的。大小接近cc-pVTZ，算極化率精度則接近昂貴得多的aug-cc-pVTZ。
Sadlej ZPOL：2004年陸續提出。對POL基組進行簡化。適合計算很大體系偶極矩和極化率，和6-311+G*耗時相仿佛，但算極化率精度比6-311++G(2df,2p)都好。
Sadlej LPol-ds：2009年提出。LPol系列基組分為LPol-ds/dl/fs/fl，依次增大。LPol-ds是其中最小的，但也比POL大得多得多，算第一超極化率精度極好，和d-aug-cc-pVTZ相仿佛，耗時則是同檔次精度中最低的。只對C、H、O、N、F有定義。
def2-SVPD、TZVPD、TZVPPD、QZVPD、QZVPPD基組：在JCP,133,134105中提出，是分別對def2-系列的SVP、TZVP、TZVPP、QZVP、QZVPP基組加上彌散函數。彌散函數的指數通過優化原子的HF極化率得到。
LFK贗勢基組：是在SBKJC贗勢基組基礎上對39個主族元素增加彌散和極化函數，目的是在贗勢基組的計算量下達到全電子Sadlej基組下的極化率計算精度。指數和收縮系數是通過逼近實驗原子極化率或高精度相對論計算值得到的。
如果某種基組目前沒有人提出帶彌散函數的版本，而計算卻又需要彌散函數，那么可以考慮將其它尺寸差不多的基組的彌散函數挪過來用。也可以通過even-tempered方式基于現有基組來構建，另外還有人提出可以將基組中指數最小的s和p的指數除以3作為s和p彌散函數的指數，見JCTC,7,3027。雖說這樣算出來的指數對特定問題肯定沒專門優化出來的指數那么好，但對于在優化指數過程中沒考慮到的更廣泛的體系和問題的類型，直接算出來的指數也有可能表現得更好。更多討論見前面提到的《給ahlrichs的def2系列基組加彌散的方法》。
4 降低彌散函數的使用量以節約時間
Pople系列基組和cc-pVnZ系列基組是目前最流行的基組。從前面的介紹可知，后者的帶彌散函數的版本額外加入的函數量比Pople系列基組的帶彌散的版本多得多，因為包含了更高角動量殼層（同時注意，角動量越大的殼層內包含的基函數數目越多，且積分越耗時），此外后者一律給氫加上彌散函數，而不是像Pople系列基組那樣通過+和++來區分。根據JCTC,7,10的統計，一般應用中aug-cc-pVnZ基函數數目比n相同的cc-pVnZ基組多了一半有余，在Gaussian下的MP2計算中，前者計算耗時是后者的2.4倍（n=D），5倍（n=T）和6.3倍（n=Q）。aug-往往帶來過大的計算負擔而使研究者難以承擔，還使收斂更困難。一些研究者發現高角動量彌散函數的價值并不大，因此考慮能否適當刪掉一些高角動量的彌散函數來緩解這些問題。不同的人提出了不同的刪減aug-cc-pVnZ彌散函數的方法。但是不同的刪減方法往往只是在個別應用中順帶進行討論，缺乏系統性。Truhlar等人在JCTC,7,10和JCTC,7,3027中給出了一套系統的刪減方法，由不同的刪減程度所得到的“月份”系列基組彌補了cc-pVnZ和aug-cc-pVnZ之間的空白，下面將具體介紹。首先說一下minimal augmentation(maug)概念，它是指對一個沒有彌散函數的基組加入最低限度的彌散函數，具體來說，就是只給重原子加一層s和p彌散函數，而不給氫加任何彌散函數。maug-cc-pVnZ基組，就是將aug-cc-pVnZ的重原子的s和p彌散函數挪到cc-pVnZ上，或者說，是把aug-cc-pVnZ的氫的所有彌散函數和重原子的角動量高于p的彌散函數都砍掉的結果。
月份基組像maug-基組一樣砍掉所有氫的彌散函數，但是重原子的高角動量彌散函數不是都砍掉，而是根據月份順序依次砍掉最高的，直到等價于maug-基組為止：
jul-cc-pVnZ：將aug-cc-pVnZ的氫的所有彌散砍掉。重元素的彌散函數不變。
jun-cc-pVnZ：將jul-cc-pVnZ的重元素的角動量最高的彌散函數砍掉。
may-cc-pVnZ：將jun-cc-pVnZ的重元素的角動量最高的彌散函數砍掉。
apr-cc-pVnZ：將may-cc-pVnZ的重元素的角動量最高的彌散函數砍掉。
...
比如對于aug-cc-pVQZ下的碳原子，最高角動量是g，所以jul-、jun-、may-和apr-分別對應于：不砍掉任何彌散函數；砍掉g彌散函數；砍掉f,g彌散函數；砍掉d,f,g彌散函數。由于apr-cc-pVQZ已經等于maug-cc-pVQZ，所以n=Q時的月份基組就沒有mar-cc-pVQZ了。如果是aug-cc-pV5Z，由于最高角動量增加了1，所以相應的月份基組的下限也會降一個月成為mar。
之所以月份基組是從jul開始，是因為完整的基組開頭的aug-可看作august。月份基組只是從中簡單地刪除函數，其它的函數的指數和收縮系數并沒有重新做優化。
眾所周知，HF/KS-DFT這樣的基于單電子有效勢的方法對于基組的角動量要求并不高，需要有彌散函數的情況下，通常使用maug-級別就夠了。aug-級別里的高角動量彌散函數給HF/KS-DFT用實在太浪費。maug-cc-pVDZ是推薦用于大體系的基組，而精度要求高一些的話，用maug-cc-pVTZ（等價于may-cc-pVTZ）就已經很好了。
后HF計算對于高角動量函數比HF/KS-DFT敏感得多，maug-的彌散函數級別顯得偏小了。對于必須需要彌散函數的情況，中低精度的計算建議用jul-cc-pVDZ，中等精度的計算建議用may-cc-pVTZ，中高精度計算建議用jul-cc-pVTZ，高精度的計算建議用jun-cc-pVQZ。之所以zeta數越多，所需精度越高，卻可以將月份適當降低，是因為隨著價層分裂數目的增加，cc-pVnZ中最外層的基函數越彌散，越能等效起到彌散函數的作用，相同zeta數下不同級別的月份基組的精度差異越小。值得一說的是，非常需要彌散函數的任務中，給cc-pVnZ哪怕僅加入maug-級別的彌散對結果的改進都要比提高zeta數要顯著得多，即彌散函數該有的時候必須有，但不需要角動量太高。
可以說，只要適當使用月份基組，就沒必要再用aug-cc-pVnZ基組了。可能有人覺得，將氫原子的全部彌散函數去掉是不是有點過分了，實際上這沒有問題。因為氫原子電負性小，電子并不富集、彌散，它周邊的重原子的彌散函數也能起到補充作用。就連高精度弱相互作用測試集S66在外推完備基組下的結果時都刪去了氫的彌散函數，更何況一般的計算。但對于個別情況，比如LiH，由于Li的電負性很小，氫會帶有很大負電荷，這時候氫也得加彌散。
正如這一節強調的，氫的彌散函數絕大多數情況沒什么用，所以很多人用比如6-311++G**的時候都應該改為6-311+G**，這樣能節省大量時間又不會使精度有什么損失。
5 月份基組在Gaussian中的使用方法
從Gaussian09 D.01版開始，直接支持月份基組。直接寫比如may-cc-pVTZ、jul-cc-pVQZ這樣就可以了。另外，Gaussian里還有個spAug-cc-pV*Z的寫法，也就是相對于cc-pV*Z給重原子加上s和p彌散（相當于maug-cc-pV*Z），但是還給氫加上s彌散。如果用的是Gaussian老版本，就得通過自定義基組方式來使用。這里以對甲烷使用may-cc-pVTZ（等價于maug-cc-pVTZ）為例進行說明。
進入BSE基組庫（https://www.basissetexchange.org），周期表中點擊碳，在左側列表里選aug-cc-pVTZ，然后Format選Gaussian94，然后點Get Basis Set按鈕，彈出的窗口有如下信息
****
C 0
S 8 1.00
8236.0000000 0.0005310
1235.0000000 0.0041080
280.8000000 0.0210870
79.2700000 0.0818530
25.5900000 0.2348170
8.9970000 0.4344010
3.3190000 0.3461290
0.3643000 -0.0089830
S 8 1.00
8236.0000000 -0.0001130
1235.0000000 -0.0008780
280.8000000 -0.0045400
79.2700000 -0.0181330
25.5900000 -0.0557600
8.9970000 -0.1268950
3.3190000 -0.1703520
0.3643000 0.5986840
S 1 1.00
0.9059000 1.0000000
S 1 1.00
0.1285000 1.0000000
S 1 1.00
0.0440200 1.0000000
P 3 1.00
18.7100000 0.0140310
4.1330000 0.0868660
1.2000000 0.2902160
P 1 1.00
0.3827000 1.0000000
P 1 1.00
0.1209000 1.0000000
P 1 1.00
0.0356900 1.0000000
D 1 1.00
1.0970000 1.0000000
D 1 1.00
0.3180000 1.0000000
D 1 1.00
0.1000000 1.0000000
F 1 1.00
0.7610000 1.0000000
F 1 1.00
0.2680000 1.0000000
****
這些是aug-cc-pVTZ中對碳原子的彌散函數的定義，沒字母的那些行中第一列是指數，第二列是收縮系數。may-是砍掉兩個最高角動量彌散函數，所以d和f彌散函數就不要了，只保留s和p的彌散函數。同等角動量中指數最小的基函數就是彌散函數，因此
D 1 1.00
0.1000000 1.0000000
和
F 1 1.00
0.2680000 1.0000000
應該被舍棄。
may-cc-pVTZ下計算CH4分子的Gaussian輸入文件應為
#P wb97xd/gen
Divokej Bill
0 1
C -2.32704417 -0.72320845 -0.00946732
H -1.97038974 -1.73201845 -0.00946732
H -1.97037133 -0.21881026 0.86418419
H -1.97037133 -0.21881026 -0.88311882
H -3.39704417 -0.72319527 -0.00946732
C 0
S 8 1.00
8236.0000000 0.0005310
1235.0000000 0.0041080
280.8000000 0.0210870
79.2700000 0.0818530
25.5900000 0.2348170
8.9970000 0.4344010
3.3190000 0.3461290
0.3643000 -0.0089830
S 8 1.00
8236.0000000 -0.0001130
1235.0000000 -0.0008780
280.8000000 -0.0045400
79.2700000 -0.0181330
25.5900000 -0.0557600
8.9970000 -0.1268950
3.3190000 -0.1703520
0.3643000 0.5986840
S 1 1.00
0.9059000 1.0000000
S 1 1.00
0.1285000 1.0000000
S 1 1.00
0.0440200 1.0000000
P 3 1.00
18.7100000 0.0140310
4.1330000 0.0868660
1.2000000 0.2902160
P 1 1.00
0.3827000 1.0000000
P 1 1.00
0.1209000 1.0000000
P 1 1.00
0.0356900 1.0000000
D 1 1.00
1.0970000 1.0000000
D 1 1.00
0.3180000 1.0000000
F 1 1.00
0.7610000 1.0000000
****
H 0
cc-pVTZ
****
這里對H使用cc-pVTZ，是由于月份基組對于氫來說相當于砍掉了所有aug-cc-pVnZ基組的彌散函數，和cc-pVnZ是等價的。
另外，BSE上目前也有月份基組的定義可以直接用。比如在BSE界面里點擊碳，直接找may-cc-pV(T+d)Z來用即可，這和我們上面自行修改得到的相同。