详解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://sobereva.com/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系列非常方便。