Gaussian98使用经验谈

    经过一段时间对GAUSSIAN的使用，发现不仅有助于加深对量子化学的了解，而且可以巩固光谱学方面的理论知识。至少在研究光谱学理论方面，GAUSSIAN是确实一个相当不错的东西。在使用中，也遇到很多问题，而且在GAUSSIAN的帮助文件中没有涉及。这些问题有的是靠自己解决的，有的是和别人探讨的，还有的是从网上或者文献中查到的。
    受本人研究方向所限，文中涉及的大多数都是处理分子激发态方面的问题。如果在分子光谱学（包括红外、拉曼、可见-紫外）方面的知识有欠缺的话，建议先看一看哈里斯和伯特卢西写的《对称性与光谱学》（胡玉才，戴寰译，高等教育出版社，1988），这是一本浅显易懂的速成书籍。
    本人水平有限，对有些问题的理解可能有误，欢迎批评指正。

（一）一些杂七杂八的小问题
§1-1.内存与速度问题
    G98W默认的内存是128兆，我以前一直用64兆，计算机是P233，所以奇慢无比。我后来发现，在计算之前把所有的杀毒程序都关闭（我的计算机上装了三个），可以节省时间20%。后来又加了一条64兆内存，时间就减少到原来的一半。如果把计算机换成P300，用128兆内存，只需要从前五分之一的时间。而且大的内存支持更大的基组。所以计算机配置可马虎不得。
§1-2.GAUSSIAN输入的一个注意事项
    输入的分子坐标（包括键长和键角）必须是实型，不能是整型，否则会报错。例如，键长2.0埃，可以写成r=2.或者r=2.0，不能写成r=2。
§1-3.GAUSSIAN的缺陷
    对于含有重原子的分子来说，光谱项的自旋、轨道角动量已经失去了意义，这个时候的光谱项应该用自旋-轨道耦合总角动量的各个分量Ω表示。例如，一氯化铊分子的3Π电子态，在实验上观察到的是两个分得很开的Ω分量：A0+态和B1态。但是GAUSSIAN仅能处理含有从H原子到Cl原子的分子（因为不是重原子，所以没有必要）。目前处理这类问题，一般是使用免费的MOLFDIR，GAMESS(USA)，NWChem，DIRAC，COLUMBUS或商业的GAMESS(UK)，MOLPRO，ADF等程序。

（二）关于基组的设置
§2-1.基组
    基组有两种，一种是全电子基组，另一种是价电子基组。价电子基组对内层的电子用包含了相对论效应的赝势（缩写PP，也称为模型势，有效核势，缩写ECP）进行近似。以铟原子为例，可以输入下面的命令来观察铟原子3-21G基组的形式：

# 3-21G GFPrint

In

0,2
In

    输出结果为：

3-21G (6D, 7F)
S 3 1.00 N=1层，s轨道
.1221454700E+05 .6124760000E-01
.1848913600E+04 .3676754000E+00
.4063683300E+03 .6901359000E+00
SP 3 1.00 N=2层，s和p轨道
.5504422550E+03 -.1127094330E+00 .1523702990E+00
.1197743540E+03 .8344349530E-01 .6096507530E+00
.3866926950E+02 .9696880360E+00 .3970249500E+00
SP 3 1.00 N=3层，s和p轨道
.4702931320E+02 -.2758954250E+00 -.1408484750E+00
.2249642350E+02 .5977348150E-01 .5290866940E+00
.6697116970E+01 .1082147540E+01 .6620681110E+00
D 3 1.00 N=3层，d轨道
.1021735600E+03 .1205559000E+00
.2839463200E+02 .4884976000E+00
.8924804500E+01 .5850190000E+00
SP 3 1.00 N=4层，s和p轨道
.6572360380E+01 .4284830560E+00 .1091304620E-01
.2502157560E+01 -.4633643610E+00 .5036758870E+00
.9420245940E+00 -.8219679320E+00 .5581808650E+00
D 3 1.00 N=4层，d轨道
.4535363700E+01 .2508574000E+00
.1537148100E+01 .5693113000E+00
.4994922600E+00 .3840635000E+00

    以上都是内层电子轨道，每个轨道用三个高斯型函数拟合，也就是3-21G中的3。

SP 2 1.00 N=5层，s和p轨道
.1001221380E+01 -.4364171950E+00 -.2316333510E-01
.1659704190E+00 .1189893480E+01 -.9903308880E+00

    这是价电子轨道，每个轨道用二个高斯型函数拟合，也就是3-21G中的2。

SP 1 1.00 N=6层，s和p轨道
.5433974090E-01 .1000000000E+01 .1000000000E+01

    这是空轨道，每个轨道用一个高斯型函数拟合，也就是3-21G中的1。
    可见3-21G是全电子基组。再看看铟原子LanL2DZ基组的结果：

LANL2DZ (5D, 7F)
S 2 1.00
.4915000000E+00 -.4241868100E+01
.3404000000E+00 .4478482600E+01
S 1 1.00
.7740000000E-01 .1000000000E+01
P 2 1.00
.9755000000E+00 -.1226473000E+00
.1550000000E+00 .1041757100E+01
P 1 1.00
.4740000000E-01 .1000000000E+01

    这些都是价电子轨道。在ECP上加上了双zeta（DZ）基组。所以LanL2DZ是价电子基组。
    在GAUSSIAN中，CEP也是很重要的价电子基组，并且98版中对元素的适用范围扩大到了86号元素Rn。需要注意的是，超过了Ar元素后，每一种元素的CEP-4G，CEP-31G，CEP-121G是等价的（但和加了极化*的CEP还是有区别的）。SDD适用于除了Fr和Ra以外的所有元素，但是结果的准确性并不太好。
    需要注意的是有些方法不需要输入基组，像各种半经验方法，G1，G2等等。
§2-2.基组的选择
    《中国科学》B30（5），2000，419页的文章给出了一种选择基组的方案，认为基函数的选择应该根据体系中不同原子的性质及实际的化学环境：
1。描述一般体系时可根据该原子在元素周期表的位置从左到右依次增加基函数用量。
2。对核负电原子，使用更多的基函数以及适当的极化函数或弥散函数。对核正电原子，基函数用量可适当减少。
3。对正常饱和共价键原子不需要增加极化或弥散函数，对氢键、弱相互作用体系、官能团、零价或低价态金属原子等敏感体系，要增加极化或弥散函数。
    这样就可以用适中的基组和可承受的计算量的到相当可靠的计算结果。方案适用于HF，MPn，DFT等计算中。
    在GAUSSIAN中有关的关键字是ExtraBasis和Gen。二者都是在routesection指定基组的基础上，再加上附加的基函数。不同的是Gen可以自己定义基组。实际上可并不那么简单。下面是一个使用ExtraBasis的例子：

# HF/3-21G ExtraBasis

AgCl

molecule specification

Cl 0
6-311G(d,p)
****

    这里对分子中所有的原子都使用3-21G，但是对Cl原子使用更大的基组6-311G(d,p)。
如果计算过程中报错，一般是由于下面两点原因造成的：
1。并不是所有的基组都可以用（哪些基组适用于哪些原子需要去参考帮助文件），即使是分别适用于不同原子的各个基组，放在一起也不一定就行，可能有匹配或收敛的问题，也有可能是BUG。所以需要反复多试。
2。基组设置的顺序问题。还是看上面的例子，我认为route section一行中设置的基组，应当适用于该分子中所有的原子（一般来说是小一些的，虽然在计算中并不对所有出现的原子都使用这个基组），而在分子坐标下面的基组只要对所定义的原子适用就可以了（一般来说是大一些的基组）。比如，如果把上面的输入变成：

# HF/6-311G(d,p) ExtraBasis

AgCl

molecule specification

Ag 0
3-21G
****

    虽然都是对Ag原子用3-21G基组，对Cl原子用6-311G(d,p)基组，但是可能会报错。当然实际情况并不总是如此，换个分子，基组哪个在前哪个在后可能并没有关系，但是注意基组顺序确实可以防止出错。这可能是一个BUG。

（三）用GAUSSIAN处理激发态
§3-1.输入文件中的电子态多重度
    输入文件中的多重度，指的是具有这种多重度的最低能量的电子态，不一定就是基态，不过一般输入文件都用基态。激发态的多重度是在输入文件的route section行控制的，如：singlets，triplets。虽然在GAUSSIAN使用帮助中没有明说，但是从给出的例子中可以体会出来。另外，在HyperChem的帮助文件中有类似的说明，虽然是针对HyperChem的，但我想对GAUSSIAN一样适用。
§3-2.GAUSSIAN中研究激发态的方法
    GAUSSIAN可以用CIS，CASSCF和ZINDO方法计算激发态。GAUSSIAN中的CIS实际是对分子激发态的零级组态相关处理。有些从头计算程序（如GAMESS）包含计算激发态的一级组态相关（FOCI）和二级组态相关（SOCI），我不清楚是不是和CIS相对应。从结果来看，FOCI和SOCI比CIS略有改善。ZINDO方法适用范围窄，只能用于周期表前48种原子构成的分子，这里不谈。98版还新增加了含时（Time Depend）方法处理激发态（包括含时Hatree-Fock(TDHF)和含时密度泛函(TDDFT)）。
    我们一般计算的激发态称为Low-lying Excited States，也就是由原子价电子构成的分子电子激发态。更高的电子态是Rydberg态。从头算法和其他半经验方法对空轨道，特别是较高空轨道的计算结果不好，因而得到的Rydberg态的结果缺乏明确的物理意义。
§3-3.含时方法（TD）输入方式
    帮助文件对TDDFT这么重要的方法说得真是太简单了，让人去参考CIS的输入方法，连个例子都没有。CIS的输入方法是：

# CIS(Triplets,NStates=10)/3-21G Test

    那就用B3LYP泛函试试看吧：

# TDB3LYP(Triplets,NStates=10)/3-21G Test

    但是不行。其实正确的输入方法是：

# B3LYP/3-21G TD(Triplets,NStates=10)Test

    这个问题浪费了我很长时间。为了表示它的重要，单独分为一节。
§3-4.CIS，TDHF，TDDFT方法激发能比较
    一般来说，CIS和TDHF方法结果的准确度是比较差的，TDDFT的结果是最好的。而各种TDDFT相比，B3LYP，B3P86，MPW1PW91是最好的（依具体的分子和使用的基组而不同）。可以参考文献：
J. Chem. Phys., 111(7), 1999, 2889
Chem. Phys. Lett., 297(1-2), 1998, 60
J. Chem. Phys., 109(23), 1998, 10180
    这些文章都是计算C，H，O组成的有机物在基态平衡位置的垂直激发能。我对含有重元素In的分子激发态进行计算，发现如果选好基组，很多TDDFT的结果和实验值的差距甚至小于500个波数，而CIS的结果则可能要差5000多个波数。所以现在使用TDDFT方法的比较多。另外还有人用TDDFT方法计算DNA这样的大分子。
§3-5.GAUSSIAN处理激发态的BUG
1。用94版进行CIS计算，50-50选项不支持LanL2DZ和3-21G*基组。98版则不会出现这个问题。这个问题可能还跟计算的分子有关。
2。同时计算单重激发态和三重激发态（也就是使用50-50这个选项），有时候会出现致命错误：某些电子态在某些位置的能量有不规则起伏，而单独计算单重态或三重态（使用Singlets或者Triplets选项）则不会出现这个问题。我在计算中曾经遇见过两次。下面这个例子是从网上看到的，我想这位德国的大哥也遇见了同样的问题。
    他对同一个任务执行了三种计算：

1) #Bp86/6-31G* td=(nstates=5,50-50) scf=direct density
2) #Bp86/6-31G* td=(triplets,nstates=1) scf=direct density
--link1--
#Bp86/6-31G* td=(nstates=10) scf=direct guess=read geom=check
--link1--
#Bp86/6-31G* td=(triplets,nstates=10) scf=direct guess=read
geom=check
3) #Bp86/6-31G* td=(nstates=10) scf=direct
#Bp86/6-31G* td=(triplets,nstates=10) scf=direct density

    第一种计算产生10个态（5S和5T），后两种产生20个。关键字Density对结果没有影响。结果是：

             1)     2)      3)
Triplet-B1U  3.1077 3.1077  3.1077
Triplet-B3G  3.9568 3.9568  3.9568
Triplet-AG   4.1612 4.1612  4.1612
Triplet-B2U  4.1769 4.1769  4.1769
Triplet-B3G  4.2101 4.2101  4.2101
Singlet-B3G  4.3609 4.3609  4.3609
Singlet-B1U  4.3788 >4.4864 >4.4864
Singlet-B2U  4.4864 >4.6409 >4.6409
Triplet-AG   4.5752 >4.7144 >4.7144
Singlet-B3G 4.6409 > -- >5.0010

    通过比较发现，对前六个态，三种处理的结果相同，但是后面的态，使用了50-50的方法就和另外两种方法的结果不一样了。
    我个人认为，尽量不要使用50-50这个选项。推荐单重态、三重态分开计算，虽谈多花一倍时间，但是结果保险。
§3-6.画激发态的势能曲线
    最笨的方法莫过于算完一个点后，更改输入文件的坐标，计算下一个点了。这样费时费力。在优化基态几何构型的时候，可以使用Scan关键字进行扫描，实际上Scan对激发态一样适用。这样开上一晚上计算机，到早上所有激发态的势能曲线数据就全有了。例如计算从R=0.6埃到1.2埃氢分子的10个单重激发态：

# RCIS(NStates=10)/6-311G Scan

H2 Excited States

0,1
H1
H2 H1 r

r 0.6 60 0.01

    需要注意三点：
1。结果给的是垂直激发能，必须换算成相对于基态平衡位置的能量。
2。扫面范围不要离基态平衡位置太远，一般从0.75Re到2Re就可以了，太远了会报错。
3。某些理论和某些基组的搭配可能不支持Scan，这时候就只好手工控制了。
    最后是数据处理，这一步是最枯燥的了。数据处理需要用到Excel、Oringin、Matlab。这些软件用不着钻研很深，知道几个最基本的命令就足够用了。
§3-7.高对称性的分子使用CIS方法
    默认的CIS计算求出最低的三个激发态。以它为例，GAUSSIAN在初始的猜测中选择电子态数目两倍的向量（也就是6个）进行迭代，直到有三个收敛。对于大多数分子体系来说，这种默认的向量数目足够了。但是对于高对称性的分子，则需要特殊的方案，否则不会得到比较高的收敛态。方法是增加NStates的数量，以增加初始向量的数量。推荐NStates值为最大阿贝尔点群操作的数量，也就是输出文件在Standardorientation之前的Symmetry部分（在输入文件中使用#控制打印输出就可以得到，#T不行）。这是一个计算苯分子的例子：

Stoichiometry C6H6
Framework group D6H[3C2'(HC,CH)]
Deg. of freedom 2
Full point group D2H NOp 8
Largest concise Abelian subgroup D2 NOp 4
Standard orientation
------------------------------
    因此在route section中使用CIS(NStates=8)。

（四）求分子的光谱常数
§4-1.求离解能
    首先必须注意，原子化能和离解能不是同一个概念。帮助文件和Exploring Chemistry with Electronic Structure Methods一书中的解释并不确切。计算原子化能和离解能是不同的两个过程。以三原子分子GaCl2为例，计算离解能的过程是：
GaCl2(X2A1)---->Cl(2P)+GaCl(X1Σ+)
过程中的能量变化称为离解能。很明显，离解能还可以细分为第一步离解能、第二步离解能......，上面的是计算第一步离解能，计算第二步离解能的过程是：
GaCl(X1Σ+)---->Cl(2P)+Ga(2P)
而计算原子化能的过程是：
GaCl2(X2A1)---->Cl(2P)+2Ga(2P)
可见还是有区别的。不过对于双原子分子来说，这两个概念是相同的。
    第二个需要注意的是，离解能一般用符号D0或者De表示（前者不包括零点能，后者包括），但是在分子光谱中，还有一组D0和De。为避免混淆，这里需要作特别说明。
    第一组D0和De，也就是离解能，数值比较大，数量级一般在0.1到1eV。这个参数多出现在真空-紫外、振动光谱中。有一个公式：
          ωe^2
ωeχe = -------             （^是指数符号）
           4De
其中的De就是离解能。
    另一组D0和De称为离心畸变常数，数值很小，作为转动光谱中的微扰项。如：
                     h^3
De = -----------------------------------
      128(π^6)(μ^3)(c^3)(re^6)(ωe^2)
其中的De就是离心畸变常数。
    GAUSSIAN中有个求PH2分子的原子化能的例子，可以作为求离解能的参考。也就是文件e7-01。结果是
原子化能 = E(P)(Hartree)+2E(H)(Hartree)-E(PH2)(Hartree)-ZPE(Hartree)
         = (-341.25930)+(-0.50027)*2-(-342.50942)-(0.01322)(Hartree)  其中的零点能ZPE经过了矫正
         = 148.3 Kcal/mol
和实验值144.7 Kcal/mol的差距相当小了。
§4-2.求电离能
    电离能缩写为I.P.。在GAUSSIAN中有两种求I.P.的方法。
法一：见文件e7-03，这个例子求PH2分子的电离能。结果为：
     E          ZPE
PH2+ -342.14416 0.01347 （ZPE经过矫正）
PH2  -342.50942 0.01322 （这是上面求原子化能的结果）
I.P. = (E(PH2+)+ZPE(PH2+))-(E(PH2)+ZPE(PH2))
     = (-342.14416+0.01347)-(-342.50942+0.01322)(Hartree)
     = 0.36551(Hartree)
     = 9.95(eV)
法二：OVGF方法。帮助文件里有详细的使用说明。
§4-3.求激发态的光谱学参数Te，ωe，Re
    有两种方法。
1。CIS可以用结构优化的方法。Starl的文章里说得很详细了，这里不再细讲。仅仅需要补充的是，如何强制保留设置后的波函：
    在使用Guess=Alter设置波函之后，程序总会回到设置前的波函。这个时候需要使用SCF关键字中的Symm选项，强制程序保留设置以后的波函。使用方式是在route section部分添加SCF=Symm。
2。势能曲线扫描法。
    TDDFT方法不能计算激发态的波函，所以不能用来研究激发态的平衡构型和分子属性。上面的结构优化方法用在TDDFT中就不适合了。但是可以用势能曲线扫描的方法。但是，多原子分子能量是势能曲面，控制多个输入坐标变量，扫描整个曲面是困难的，所以这个方法用于双原子分子。如何扫描激发态的势能曲线见前面的§3-6节。
    有了势能曲线以后，就可以找到势能最低点Re（为了精确，这一区域的点要取密一些），Re位置的能量相对于基态平衡位置能量的差值就是电子态高度Te。势阱的深度就是离解能De。
    这些都简单，关键是求振动参数ωe。求ωe也有两种方法。一种是用Morse势函数拟合的方法。Morse势函数形式为：
V(R) = Te+De{1-EXP[-α(R-Re)]}^2 （对基态Te是0）
Te、De、Re是已知的，拟合出了α后，用公式
α = ωe ((2(π^2)c*μ)/(h*De))^0.5 （高斯单位制）
就可以得到ωe。
    这种方法太麻烦，我们不这样做。实际只要把扫描出来的Morse势函数对R求二阶导数，在Re处的导数值f2和ωe的关系有：
       4(π^2)(ωe^2)c*μ
f2 = --------------------
          h*(10^18)
    其中c、h、μ是国际单位制，ωe的单位是波数cm^-1。另外在扫描的势函数中，De的单位是波数cm^-1，R和Re的单位是埃。否则上面的f2需要乘上个系数。
    求f2用Origin软件就可以做。注意要取Re位置的f2。有了f2，就可以算出ωe，省去了计算De和拟合α的步骤，避免引入更大的误差。
§4-4.如何获得光谱项的对称性
    如果用GAUSSIAN计算比较轻的分子，一般都能给出正确的激发态光谱项对称性，但是对含有重原子的分子，可能得不到对称性。这个时候就需要根据轨道跃迁的的信息获得对称性。
    下面是用CIS方法LanL2DZ基组计算InCl分子的例子，输出文件给出InCl的电子轨道信息：

Occupied (SG) (SG) (PI) (PI) (SG)
Virtual (PI) (PI) (SG) (PI) (PI) (SG) (PI) (PI) (SG) (SG) (SG)

    这里结果只给出前4个单重激发态：

CIS wavefunction symmetry could not be determined.
Excited State 1: Singlet-?Sym 4.9226 eV 251.87 nm
f=0.3243
5 -> 6 .43587
5 -> 7 .54909
This state for optimization and/or second-order correction.
Copying the Cisingles density for this state as the 1-particle
RhoCI density.

CIS wavefunction symmetry could not be determined.
Excited State 2: Singlet-?Sym 4.9226 eV 251.87 nm
f=0.3243
5 -> 6 .54909
5 -> 7 -.43587

Excited State 3: Singlet-SG 6.7609 eV 183.38 nm
f=0.2830
3 -> 6 -.28038
4 -> 7 -.28038
5 -> 8 .57818

CIS wavefunction symmetry could not be determined.
Excited State 4: Singlet-?Sym 6.8241 eV 181.68 nm
f=0.0000
3 -> 7 .49690
4 -> 6 -.49690

    由于In原子是重原子，所以大多数电子态并没有给出对称性（第三个激发态除外）。我们发现第一激发态是从第五个轨道跃迁到第六、七个轨道，而第六、七两个轨道是简并的（PI、DELT轨道都是二重简并的）。也就是说第一激发态的电子组态是：

InCl: (内层电子)(1σ)2 (2σ)2 (1π)4 (3σ)1 (2π)1

    这个组态耦合出来的单重态光谱项有一种（当然还有一个三重的）：1Π，所以就是1Π态。再往下看，第二激发态和第一激发态的电子组态、能量完全相同，说明这是一个简并态，而实际上除了Σ态以外的态都是二度简并的。所以第一激发态和第二激发态实际都是InCl的第一激发态1Π态。
    第三个激发态给出对称性了。如果没给，方法和上面的一样，先找到它的电子组态，它是由两个组态叠加而成的：

InCl: (内层电子)(1σ)2 (2σ)2 (1π)3 (3σ)2 (2π)1
InCl: (内层电子)(1σ)2 (2σ)2 (1π)4 (3σ)1 (4σ)1

    第一个组态可以耦合出的单重态光谱项有：1Δ，1Σ+，1Σ-。因为这不是一个简并态，所以只能是1Σ+和1Σ-其中之一。这就需要进一步知道+、-对称性。方法是看跃迁系数：

3 -> 6 -.28038
4 -> 7 -.28038

    由于正负号相同，所以是1Σ+态。那么第四激发态必然是1Σ-了。另外还可以看第二个组态耦合出来的光谱项，只有1Σ+，因此第三个激发态必然是1Σ+。两个结果一致。

    以上是自旋限制的计算。非限制计算（每个轨道再分成自旋Alpha和自旋Beta）分析类似。

***************************************************************************
Gaussian 98的补充

1. 谐振强度f
    谐振强度f出现在计算激发态的结果中，如：

Excited State 3: Singlet-?Sym 5.0236 eV 246.80 nm f=0.3428
5 -> 7 0.68918

    在Gaussian 98早期版本中，TD方法计算激发态的谐振强度f是错误的（其它计算激发态的方法没有这个问题）。从A7版开始做了修正。这点要注意。
    谐振强度能干什么呢？可以预测电子激发态跃迁的寿命。用公式
t=1.499/(f*E^2)
    这里的E是电子态之间的能量间隔，单位是波数，f是电子态谐振强度，t是电子态寿命，单位是秒。默认的下能级是基态。如果计算各激发态之间的f，在Route Section中加入关键字AllTransitionDensities，就得到各电子态之间的偶极矩m(x)，m(y)，m(z)。求解f用公式：
f=(2/3)*E*[m(x)^2+m(y)^2+m(z)^2]
    因为Gaussian给出的m(x)，m(y)，m(z)是au单位，所以激发态之间的能量间隔E要换算成au单位，才能代入上面的公式。
2. 用TDDFT做几何结构优化
    一般认为，因为TDDFT不能计算激发态的波函，所以不能用来研究激发态的平衡构型和分子属性。所以用Gaussian 98提供的TDDFT方法做研究激发态的研究工作，只能得到垂直激发能，对于计算平衡键长Re、键角和振动频率ωe无能为力。
    不过，最近的报道说，TDDFT做激发态几何优化在理论上证明可行。已经有消息说，Gaussian公司准备在新版本的Gaussian中加入TDDFT几何优化了。可以参考文献：
Van Caillie C, Amos R. D., Geometric derivatives of density functional theory excitation energies using gradient-correct functionals, Chem. Phys. Lett., 317(1-2), 2000, p.159-164
和
Hsu, Hirata, Head-Gorden, J. Phys. Chem., A105, 2001, 451-458 关于线性多烯烃的TDDFT激发能和谐振强度。
3. 自定义赝势（也称为模型势，有效核势）
    Gaussian给出的Pseudo示例如下：

# Becke3LYP/Gen Pseudo=Read Opt Test

HF/6-31G(*)Opt of Cr(CO)6

0 1
Cr 0.0 0.0 0.0
molecule specification continues ...

C O 0
6-31G(d)
****

Cr 0
LANL2DZ
****

Cr 0
LANL2DZ

    这里，Cr原子的基组LanL2DZ输入了两次。这两个LanL2DZ是不同的。第一个描述的是价电子，或者核势之外的电子；第二个描述内层电子，也就是核势部分。如果自定义ECP基组，也这么输入。当然还有其他输入方法。
    这是在CIS计算中自定义In原子的基组，In的ECP使用程序内部定义的SDD赝势，Cl原子使用Route Section中的LanL2DZ。

# CIS(NStates=3,50-50)/LanL2DZ ExtraBasis pseudo=read

InCl Excited States

0,1
In
Cl In r

r=2.401

IN 0
S 3 1.00
1.53256200 0.261774000
0.944269000 -0.625406000
0.176733000 0.651304000
S 1 1.00
0.643430000E-01 1.00000000
P 3 1.00
0.641474000 0.211103000
0.342739000 -0.243029000
0.109546000 -0.560191000
P 1 1.00
0.389580000E-01 1.00000000
****

In 0
SDD

    如果自定义ECP的话，可以把上面的SDD换成：

IN 0
IN-ECP 4 46
G POTENTIAL
1
2 1.00000000 0.00000000
S-G POTENTIAL
2
2 1.43509100 29.16521900
2 0.69580500 -4.19080600
P-G POTENTIAL
2
2 1.44083200 36.99054200
2 0.70139200 -3.36582000
D-G POTENTIAL
1
2 0.96123600 20.00053100
F-G POTENTIAL
1
2 0.88436900 -6.01909200

    价电子基组和核势顺序不要颠倒。下载更多的ECP基组，可以到：

http://www.emsl.pnl.gov/forms/basisform.html

    区分核势和价电子基组的方法是：核势第一行有一个被核势代替的内壳层电子个数，而且就是固定的那么几个偶数。像上面的例子，核势代替In的46个内壳层电子。
    如何得到Gaussian中ECP基组的赝势呢？只要在Route Section行添加IOp(3/18=1)就行了。IOp的功能是极其强大的，可以参考Gaussian公司web页。这只是其中之一。例如要输出Cl原子完整的ECP基组CEP-121G，就必须用到GFPrint和IOp(3/18=1)关键字，像下面这样：

# HF/CEP-121G GFPrint IOp(3/18=1)

    其中GFPrint输出CEP-121G的价电子基组部分。当然，如果观察cc-pV[DT]Z，6-31G等非ECP基组，就不用加IOp(3/18=1)了，仅用GFPrint就行。这是输出结果中Cl原子的CEP-121G基组部分，第一部分是GFPrint结果，第二部分是IOp(3/18=1)的结果：

****************************************************************************
ATOMIC ORBITAL * GAUSSIAN FUNCTIONS
****************************************************************************
FUNCTION SHELL SCALE *
NUMBER TYPE FACTOR * EXPONENT S-COEF P-COEF D-COEF F-COEF
****************************************************************************

1- 4 SP 1.00
0.222500D+010.100000D+010.100000D+010.000000D+000.000000D+00
5- 8 SP 1.00
0.117300D+010.115280D+000.299520D+000.000000D+000.000000D+00
0.385100D+000.847170D+000.583570D+000.000000D+000.000000D+00
9- 12 SP 1.00
0.130100D+000.100000D+010.100000D+010.000000D+000.000000D+00
****************************************************************************
============================================================================
PSEUDOPOTENTIAL PARAMETERS
============================================================================
CENTER ATOMIC VALENCE ANGULAR POWER
NUMBER NUMBER ELECTRONS MOMENTUM OF R EXPONENT COEFFICIENT
============================================================================
1 17 7
D and up
1 4.8748300 -3.40738000
S - D
0 17.0036700 6.50966000
2 4.1038000 42.27785000
P - D
0 8.9002900 3.42860000
2 3.5264800 22.15256000
===========================================================================

    以上内容大部分来自www.ccl.net讨论区。里面内容除少部分是讨论使用GAMESS及其它程序外，大多是关于Gaussian的，而且Gaussian公司的几个人也经常参加讨论。建议大家有机会去看看。