新浪博客

Gaussian98使用经验谈（一）

2007-01-12 17:22阅读：

http://blog.sina.cn/dpool/blog/u/1259836332

经过一段时间对GAUSSIAN的使用，发现不仅有助于加深对量子化学的了解，而且可以巩固光谱学方面的理论知识。至少在研究光谱学理论方面，GAUSSIAN是确实一个相当不错的东西。在使用中，也遇到很多问题，而且在GAUSSIAN的帮助文件中没有涉及。这些问题有的是靠自己解决的，有的是和别人探讨的，还有的是从网上或者文献中查到的。
受本人研究方向所限，文中涉及的大多数都是处理分子激发态方面的问题。如果在分子光谱学（包括红外、拉曼、可见-紫外）方面的知识有欠缺的话，建议先看一看哈里斯和伯特卢西写的《对称性与光谱学》（胡玉才，戴寰译，高等教育出版社，1988），这是一本浅显易懂的速成书籍。
本人水平有限，对有些问题的理解可能有误，欢迎批评指正。
******************
目录
******************
（一）一些杂七杂八的小问题
§1-1.内存与速度问题
§1-2.GAUSSIAN输入的一个注意事项
§1-3.GAUSSIAN的缺陷
（二）关于基组的设置
§2-1.基组
§2-2.基组的选择
（三）用GAUSSIAN处理激发态
§3-1.输入文件中的电子态多重度
§3-2.GAUSSIAN中研究激发态的方法
§3-3.含时方法（TD）输入方式
§3-4.CIS，TDHF，TDDFT方法激发能比较
§3-5.GAUSSIAN处理激发态的BUG
§3-6.画激发态的势能曲线
§3-7.高对称性的分子使用CIS方法
（四）求分子的光谱常数
§4-1.求离解能
§4-2.求电离能
§4-3.求激发态的光谱学参数Te，ωe，Re
§4-4.如何获得光谱项的对称性
[附]Gaussia
n 98的补充
1. 谐振强度f
2. 用TDDFT做几何结构优化
3. 自定义赝势（模型势，有效核势）
（一）一些杂七杂八的小问题
§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)。

我的更多文章

下载客户端阅读体验更佳

APP专享