詳解Gaussian中混合基組、自定義基組和贗勢基組的輸入

Gaussian的初學者經常在混合基組、自定義基組和贗勢的輸入方面犯錯，此文就具體說說怎么正確用它們。
 

1 混合基組的輸入

混合基組就是指不同的原子用不同的基組。寫gen關鍵詞代表從坐標部分后面讀入基組定義。用IOp(3/24=1)可以顯示每個原子實際用的基組，檢查是否設對了；也可以用GFPrint關鍵詞輸出相似內容，但沒那么清楚；也可以用GFInput關鍵詞來輸出實際所用基組，這些信息可以直接作為基組輸入的信息。

下面用的例子是CH4，例如想讓C是6-311G*，H是6-31G**，則輸入文件為
# m062x/gen
[空行]
test
[空行]
0 1
 C                 -0.00000000    0.00000000    0.00000000
 H                 -0.00000000    0.00000000    1.09000000
 H                 -0.00000000   -1.02766186   -0.36333333
 H                 -0.88998127    0.51383093   -0.36333333
 H                  0.88998127    0.51383093   -0.36333333
[空行]
C 0   <--元素或者原子序號寫完之后寫個0示意寫完了，但不寫也沒關系
6-311G*
****
H 0   <--等價于寫成2 3 4 5，也可以寫成2-5或者2-3 4-5或者2 3-5等，用-表示范圍在原子多時會方便很多
6-31G**
****  <--即便是基組定義的最后一部分，也要寫這個分隔符
定義基組時，寫的順序完全是無所謂的，可以先定義H，也可以先定義C。

必須確保每個原子都定義了基組。另外也不能多定義基組，比如將C 0改成C F 0來說明F也用6-311G*基組，由于體系里沒F，會因此報錯。但可以在元素符號前寫上-，這說明分子中有這個元素則基組定義生效，如果沒有這個元素也不報錯。例如可以將上面C 0改成C -F 0，也等價于在末尾添加這一項
-F 0
6-311G*
****
都不會引起報錯。

如果只想讓H2、H3、H4是6-31G**，H5是用6-311G*，不能這樣寫：
C H 0
6-311G*
****
2 3 4 0
6-31G**
****
這樣寫后果是H2、H3、H4同時被賦予6-311G*和6-31G**定義的基函數，是錯誤的。因為對同一個原子重復定義基組的話基函數是累加的關系，而不是覆蓋的關系。為了正確達到目的，應當這樣寫：
C 5 0      <--元素符號和原子序號可以同時出現，等價于1 5 0
6-311G*
****
2 3 4 0
6-31G**
****
 
 

2 自定義基組的形式

這里說的自定義基組，是指的自行輸入、修改基組的具體定義，而不是直接用Gaussian內置的現成的基組。這需要對基組的定義形式有基本的了解，如果尚一無所知的話，建議看看《基組入門資料小合集》（http://bbs.keinsci.com/forum.php?mod=viewthread&tid=2190）。

基組的定義可以從BSE基組數據庫上獲得，也就是進入https://www.basissetexchange.org，點擊元素符號，左側點擊要用的基組，然后Format選Gaussian，然后點Get Basis Set按鈕，把里面的數據拷下來即可。Gaussian自己也有基組庫，對于Linux版本是能直接看到的，也就是Gaussian目錄下的basis子目錄。比如用文本編輯器打開其中的6311.gbs，就會看到6-311G對各種元素的定義，而6311s.gbs里記錄了6-311G的極化函數對各種元素的定義，若將它們和從BSE上拷下來的基組定義相對比，會發現是一致的。

這里以碳原子的6-31+G*基組為例說明基組在Gaussian中的定義格式。
C     0
S   6   1.00     <--下面定義S型基函數殼層，其收縮度為6，即這個基函數殼層由下面寫的6個primitive殼層收縮而成。1.00是刻度因子，不用去管它，總是輸入1.00就行了（具體含義參見手冊里gen關鍵詞）
   3047.5249000              0.0018347  <---- 第1個primitive殼層的指數和收縮系數
    457.3695100              0.0140373  <---- 第2個primitive殼層...
    103.9486900              0.0688426        
     29.2101550              0.2321844        
      9.2866630              0.4679413        
      3.1639270              0.3623120        
SP   3   1.00     <--SP型基函數殼層，僅在Pople系列基組及其它少數基組中才出現，其中S和P型基函數共享相同指數但收縮系數不同，在Gaussian中用SP殼層比起單獨用一個S殼層和一個P殼層計算起來更快
      7.8682724             -0.1193324              0.0689991   <--共用的指數，S的收縮系數，P的收縮系數
      1.8812885             -0.1608542              0.3164240        
      0.5442493              1.1434564              0.7443083        
SP   1   1.00
      0.1687144              1.0000000              1.0000000        
SP   1   1.00   <--此殼層指數很小，所以顯然是6-31+G*中的+號對應的彌散函數。彌散函數、極化函數的收縮度總是為1，即未收縮，這樣的殼層的收縮系數必定是1.0。
      0.0438000              1.0000000              1.0000000        
D   1   1.00   <--D型極化函數，對應6-31+G*中的*。
      0.8000000              1.0000000        
****

用自定義基組時，也是寫上gen關鍵詞，然后在分子坐標末尾空一行再寫上如上形式的基組定義即可。這里我們把BSE上拷下來STO-3G的定義用到甲烷上
# m062x/gen
[空行]
test
[空行]
0 1
 C                 -0.00000000    0.00000000    0.00000000
 H                 -0.00000000    0.00000000    1.09000000
 H                 -0.00000000   -1.02766186   -0.36333333
 H                 -0.88998127    0.51383093   -0.36333333
 H                  0.88998127    0.51383093   -0.36333333
[空行]
H     0
S   3   1.00
      3.42525091             0.15432897       
      0.62391373             0.53532814       
      0.16885540             0.44463454       
****
C     0
S   3   1.00
     71.6168370              0.15432897       
     13.0450960              0.53532814       
      3.5305122              0.44463454       
SP   3   1.00
      2.9412494             -0.09996723             0.15591627       
      0.6834831              0.39951283             0.60768372       
      0.2222899              0.70011547             0.39195739       
****

 

3 自定義基組與混合基組聯用

實際上，Gaussian會先把每個元素符號展開成原子序號，把每個基組名字展成基組的具體信息。第一節的例子的基組定義部分也可以等價地寫成下面這樣，這里直接寫明C的6-311G*基組的具體定義，而H則是直接寫基組名。
C     0
S   6   1.00
   4563.2400000              0.00196665       
    682.0240000              0.0152306        
    154.9730000              0.0761269        
     44.4553000              0.2608010        
     13.0290000              0.6164620        
      1.8277300              0.2210060        
SP   3   1.00
     20.9642000              0.1146600              0.0402487        
      4.8033100              0.9199990              0.2375940        
      1.4593300             -0.00303068             0.8158540        
SP   1   1.00
      0.4834560              1.0000000              1.0000000        
SP   1   1.00
      0.1455850              1.0000000              1.0000000        
D   1   1.00
      0.6260000              1.0000000
****
H 0
6-31G**
****
 

4 給已有基組額外添加函數

需要額外添加函數的情況主要是用來添加彌散和極化函數。

對于第一節的CH4的例子，如果給原本使用6-311G*的C加上彌散函數，當然最簡單的辦法就是改寫成6-311+G*就行了，但也可以通過自定義基組的方式給它額外加上彌散函數。我們可以從BSE上找到C的6-311+G*基組中對應彌散函數的那部分，也可以從basis目錄下的plus.gbs文件（Pople基組的重原子的彌散函數的定義文件）中把C的那部分信息拷下來，然后添到C的6-311G*基組定義的后面，即
C 0
6-311G*
SP   1 1.00
 0.4380000000D-01  0.1000000000D+01  0.1000000000D+01
****
H 0
6-31G**
****

也可以這么寫，和上面是等價的，因為定義了兩遍C的基組，相當于累加在了一起
C 0
6-311G*
****
H 0
6-31G**
****
C 0
SP   1 1.00
 0.4380000000D-01  0.1000000000D+01  0.1000000000D+01
****

類似地，我們可以給指定的原子增加極化函數。

很多文獻里都給出了對于某元素的某基組上額外添加的極化/彌散函數的指數。由于我們已經知道了其殼層類型是什么，而且已知其收縮度必為1，收縮系數和刻度因子都為1.0，因此就能直接按照自定義基組標準格式將這些極化/彌散添加到當前基組中。

PS：在定義基組時，往往會看到像上面這樣Dxxx的后綴，這是Gaussian這樣的Fortran程序用的雙精度浮點數的表示法，D相當于科學計數法里的E，比如0.4380000000D-01就代表0.438*10^-1=0.0438。對于小數，強烈建議總以雙精度表示法來書寫，因為Gaussian都是用雙精度浮點。如果在基組定義里直接寫0.0438，那么讀入Gaussian時會被當成單精度浮點數來讀入到雙精度變量里，由于單精度浮點的有效位數比較有限，實際讀入后可能就成了比如0.04380003215這樣，在末尾額外冒出來一些“噪音”影響計算精度。

另外，Gaussian里還有個extrabasis關鍵詞，它與gen的區別在于，gen是完全重新定義所有原子的基組，而extrabasis關鍵詞只是用來給當前基組額外添加一些基函數，只需要定義額外添加的那些基函數即可。假設當前用的是6-31G*，我們想給碳加上彌散函數，可以把關鍵詞部分寫為# M062X/6-31G* extrabasis，然后分子坐標末尾空一行寫上
C 0
SP   1 1.00
 0.4380000000D-01  0.1000000000D+01  0.1000000000D+01
****
可見，extrabasis和gen用于添加額外函數時的明顯差異在于前者不需要把原本的基組在分子坐標后面重新寫一遍，但缺點是原先的基組不能是混合基組的，因為它不能和gen一起混用，沒法同時定義原先的混合基組。
 
 

5 技巧：利用文件引用設定基組

前面提到Linux版的basis文件夾里的.gbs文件可以查到各種Gaussian內置的基組對各種元素的定義，通過文件名可大致猜到對應什么基組。例如631.gbs對應6-31G*，內容為：
-H
S    3 1.00
 0.1873113696D+02  0.3349460434D-01
 0.2825394365D+01  0.2347269535D+00
 0.6401216923D+00  0.8137573262D+00
S    1 1.00
 0.1612777588D+00  0.1000000000D+01
****
-He
S    3 1.00
 0.3842163400D+02  0.2376600000D-01
...（略）

實際上這些文件的信息不光可以被查閱、拷貝，還可以通過Gaussian的文件引用功能把它們直接弄進輸入文件里。例如：
{分子坐標信息}
{空行}
@/sob/g03/basis/sto3g.gbs  <--如果末尾寫上/N，可以避免在輸出中顯示一遍這文件的內容。

這里@代表引用外部文件，在Gaussian的Link0模塊中這個條目會被展開成實際文件內容，如同編程中的include一樣。當然直接把.gbs文件內容復制到輸入文件里也可以。.gbs里面定義的元素眾多，肯定不會同時出現在當前分子中，由于此文件里每個元素符號前都有負號，所以多定義的那些元素不會引起報錯，而是直接被忽略掉。

利用Gaussian的文件引用方法會給實際研究帶來很大的便利，比如要研究一大批分子，都要用到某自定義基組，那么就用不著把基組定義在每個文件里都完整寫一遍，而是都直接引用一個記錄了基組定義的文件即可。在《給def2以ma-方式加彌散函數的Gaussian格式的基組定義文件（含所有def2支持的元素）》（http://www.shanxitv.org/509）一文中也演示了引用外部文件來方便地用ma-def2系列基組。
 
 

6 使用贗勢基組

贗勢基組需要和對應的贗勢一起使用，如果不熟悉贗勢或想深入了解內在細節，可參見《贗勢的函數形式以及在量子化學程序中定義的方式》（http://sobereva.com/188）。使用贗勢基組需要用genecp關鍵詞（等價于gen pseudo=read），代表先從分子坐標后面讀取贗勢基組定義，再讀取贗勢定義。例如Cu(CO)+：
# B3LYP/genecp
[空行]
Cu(CO)+
[空行]
1 1
 Cu             
 C                  1            B1
 O                  2            B2    1            180.
[空行]
   B1             1.94000000
   B2             1.11540000
[空行]
Cu 0
Lanl2DZ    <-- Cu用Lanl2DZ贗勢基組
****
C O 0
6-31G*   <-- C O用的基組為6-31G*
****
[空行]
Cu 0
Lanl2   <-- Cu用的贗勢為Lanl2，這種贗勢是Lanl2DZ贗勢基組對應的贗勢。如果不知道贗勢名是什么，此處直接寫贗勢基組名就可以。比如這里也可以寫成Lanl2DZ，和寫Lanl2是等價的
[空行]
[空行]
     

Stuttgart贗勢的用法相對特殊一些，這里專門提一下。Stuttgart贗勢規范的寫法是nXY這種形式，含義是
n = 被贗勢表示的核電子數
X = S/M：擬合贗勢時的參考體系。S是只考慮模型體系僅有的單個價電子的能量，M是考慮實際原子全部價電子的能量
Y = HF/WB/DF：擬合贗勢用的數據的計算級別。HF=Hartree-Fock非相對論（即此贗勢不考慮相對論效應），WB=Wood-Boring準相對論，DF=Dirac-Fock全相對論
Gaussian自帶了這些組合：SDF2,MWB2,SDF10,MWB10,MDF10,MWB28,MWB46~60,MWB78。每種元素并不是所這些都能用，比如P只能用MWB10和SDF10，La只能用MWB28/46/47和MHF46/47（MWB28對La屬于小核贗勢，其它的屬于大核贗勢），詳見手冊pseudo關鍵詞部分的列表。
這些Stuttgart贗勢的寫法往往不容易記憶，在Gaussian中人們習慣用SDD關鍵詞。寫SDD說明對序號<=Ar的元素都用D95V或6-31G*全電子基組，而對更重的元素都使用Stuttgart贗勢和相應贗勢基組，具體用的是nXY寫法中的哪種，參見pseudo關鍵詞的表格。另外，還有個SDDAll關鍵詞，與SDD的差別在于它對所有序號>2的元素都使用Stuttgart贗勢和相應贗勢基組，而非對于輕元素用D95V代替。
Stuttgart系列贗勢在不斷發展完善，而在Gaussian中所內置（包括BSE上的）的Stuttgart贗勢既不全也不新。最新版本可以從Stuttgart系列贗勢的開發者的網站上下載到：http://www.tc.uni-koeln.de/PP/clickpse.en.html。

一個使用Stuttgart贗勢的例子：
Cu 0
SDD
****
C O 0
6-31G*
****
[空行]
Cu 0
SDD
通過手冊pseudo關鍵詞部分的表格可以得知，這兩處SDD也可以都寫為MDF10，是等價的。


對于不同元素可以使用不同贗勢，在贗勢定義部分寫多種即可，例如對鑭采用MWB46，對镥采用MWB60，應寫為
C H O N 0
6-31G**
****
La
MWB46
****
Lu
MWB60
****
[空行]
La  0
MWB46
Lu 0
MWB60
[空行]
[空行]


使用前面介紹的方法，同樣可以對贗勢基組添加額外的函數。很常見的一個情況就是對Lanl2DZ贗勢基組增加極化函數。對于Cu，加極化就是指加f函數，其指數散見在一些文獻中，也可以從BSE上的Cu的Lanl2TZ(f)基組中借來，數值是3.525。加上f極化后分子坐標后面內容就成了
Cu 0
Lanl2DZ
F 1 1.0
3.525D0 1.0D0     <---D0后綴用來將數值轉化成雙精度浮點表示
****
C O 0
6-31G*
****
[空行]
Cu 0
Lanl2
[空行]
[空行]


有時候我們需要用到Gaussian沒有內置的贗勢基組。比如Lanl08、Lanl2TZ、cc-pVnZ-PP系列等，這時候我們需要從BSE上拷定義，貼到贗勢基組和贗勢定義的地方。比如對上例的Cu用cc-pVDZ-PP，就去BSE找相應條目，頁面里顯示的信息中前一半是贗勢基組定義，后一半是贗勢定義，用到此例后分子坐標后面部分就成了
Cu     0
S   7   1.00
    560.0880000              0.0006370        
     56.6486000             -0.0097350        
     35.4258000              0.0657930        
     11.0546000             -0.4150350        
      2.3068200              0.7466110        
      0.9514290              0.4621730        
      0.1451840              0.0159830        
...（略）        
D   1   1.00
      0.2836420              1.0000000        
F   1   1.00
      2.7482000              1.0000000
****
C O 0
6-31G*
****
[空行]
CU     0
CU-ECP     4     10
g-ul potential
  1
2      1.0000000              0.0000000        
...（略）      
f-ul potential
  2
2      6.1901020             -0.2272640        
2      8.1187800             -0.4687730   
[空行]
[空行]

def2系列基組也比較特殊，這里特意說一下。此系列基組對前四周期（H~Kr）是全電子基組，對第五、六周期（不包含錒系）是贗勢基組，搭配的是小核Stuttgart贗勢。在Gaussian中使用def2系列且涉及到第五周期及以后的元素時，不需要用genecp來特意給第五周期及以后的元素定義贗勢。比如說，要優化OsO4，就直接寫# B3LYP/def2SVP opt就行了，此時Gaussian自動就知道對O使用def2SVP全電子基組，對Os使用def2SVP贗勢基組并加上相應的贗勢。所以用def2系列非常方便。