| 查看: 5220 | 回复: 22 | ||||||
| 【奖励】 本帖被评价17次，作者zzgyb增加金币 16 个 | ||||||
| 当前主题已经存档。 | ||||||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||||||
[资源]
guassian使用经验
|
||||||
|
经过一段时间对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原子的分子（因为不是重原子，所以没有必要），只能算是80年代 初期的水平。目前处理这类问题，除了自己编程，一般是使用MOLFDIR 程序。 （二）关于基组的设置 §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以外的所有元素，但是结果的 准确性并不太好。另外，在下面的两个网址，有几乎整个周期表元素的 赝势下载，不过我想不出那些短寿命的超铀元素有什么用： http://www.mapr.ucl.ac.be/ABINIT http://www.theochem.uni-stuttgart.de 需要注意的是有些方法不需要输入基组，像各种半经验方法，G1， G2等等。 §2-2.基组的选择 《中国科学》B30（5），2000，419页的文章给出了一种选择基组 的方案，认为基函数的选择应该根据体系中不同原子的性质及实际的化 学环境： 1。描述一般体系时可根据该原子在元素周期表的位置从左到右依 次增加基函数用量。 2。对核负电原子，使用更多的基函数以及适当的极化函数或弥散 函数。对核正电原子，基函数用量可适当减少。 3。对正常饱和共价键原子不需要增加极化或弥散函数，对氢键、 弱相互作用体系、官能团、零价或低价态金属原子等敏感体系，要增加 极化或弥散函数。 这样就可以用适中的基组和可承受的计算量的到相当可靠的计算结 果。方案适用于HF，MPn，DFT等计算中。 在GAUSSIAN中有关的关键字是ExtraBasis和Gen。二者都是在route section指定基组的基础上，再加上附加的基函数。不同的是Gen可以自 己定义基组。实际上可并不那么简单。下面是一个使用ExtraBasis的例 子： # HF/3-21G ExtraBasis AgCl molecule specification Cl 0 6-311G(d,p) **** 这里对分子中所有的原子都使用3-21G，但是对Cl原子使用更大的 基组6-311G(d,p)。 如果计算过程中报错，一般是由于下面两点原因造成的： 1。并不是所有的基组都可以用（哪些基组适用于哪些原子需要去参考 帮助文件），即使是分别适用于不同原子的各个基组，放在一起也不 一定就行，可能有个匹配的问题。所以需要多试。 2。基组设置的顺序问题。还是看上面的例子，我认为route section 一行中设置的基组应该适用于分子中所有的原子（或者应该是比较小 的基组）——虽然在计算中并不对所有的原子都使用这个基组，而在 分子坐标下面定义的基组只要对该原子适用就可以了（一般来说是大 一些的基组）。比如，如果把上面的输入变成： # HF/6-311G(d,p) ExtraBasis AgCl molecule specification Ag 0 3-21G **** 虽然都是对Ag原子用3-21G基组，对Cl原子用6-311G(d,p)基组，但是可 能会报错。当然实际情况并不总是如此，Ag和Cl的基组哪个在前哪个在 后可能并没有太大关系，但是有时候注意基组顺序确实可以解决问题。 （三）用GAUSSIAN处理激发态 §3-1.输入文件中的电子态多重度 输入文件中的多重度，指的是具有这种多重度的最低能量的电子态， 不一定就是基态，不过一般输入文件都用基态。激发态的多重度是在输 入文件的route section行控制的，如：singlets，triplets。虽然在 GAUSSIAN使用帮助中没有明说，但是从给出的例子中可以体会出来。另 外，在HyperChem的帮助文件中有类似的说明，虽然是针对HyperChem的， 但我想对GAUSSIAN一样适用。 §3-2.GAUSSIAN中研究激发态的方法 GAUSSIAN可以用CIS或者ZINDO方法计算激发态，CIS能够研究的分 子种类与使用的基组有关，GAUSSIAN中的CIS实际是对分子激发态的零 级处理。有些从头计算程序包含计算激发态的一级组态相关（FOCI）和 二级组态相关（SOCI），我不清楚是不是和CIS相对应。但是，从结果 来看，FOCI和SOCI并没有比CIS改善多少。ZINDO方法适用范围窄，只能 用于周期表前48种元素构成的分子，这里不谈。98版还新增加了含时 （TD）方法处理激发态（包括含时Hatree-Fock和含时密度泛函）。 我们一般计算的激发态称为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和TDDFT方法结果的准确度是比较差的，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值为最大阿贝尔点群操作的数量，也就是输出文件在Standard orientation之前的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 [ Last edited by zzgyb on 2008-1-10 at 12:03 ] |
回复此楼» 收录本帖的淘帖专辑推荐
 Gaussian | 化学 | 基组的介绍总结 | gauss |
» 猜你喜欢
26年储能、电池方向博士申请
已经有1人回复
-
求调剂，本科吉大，一志愿华工339
已经有2人回复
-
物理化学论文润色/翻译怎么收费?
已经有176人回复
-
意大利博洛尼亚大学招聘研究助理2名（2D材料/锂离子与钠离子电池方向）
已经有1人回复
-
关于碳量子点抗菌的研究
已经有15人回复
-
南大化学调剂
已经有11人回复
-
0703调剂
已经有2人回复
-
1
已经有20人回复
-
0703 调剂
已经有5人回复
-
西南某211招27届博士
已经有23人回复
» 本主题相关商家推荐: (我也要在这里推广)
5楼2008-01-28 09:52:18
★ ★
小红豆(金币+2,VIP+0):非常感谢！辛苦了！
小红豆(金币+2,VIP+0):非常感谢！辛苦了！
|
（四）求分子的光谱常数 §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Π态。再往下看，第二激发态和第一激发态的电子组态、 能量完全相同，说明这是一个简并态，而实际上除了Σ态以外的态都是 二度简并的。所以第一激发态和第二激发态实际都是InCl的第一激发态 1Π态。 第一个组态可以耦合出的单重态光谱项有：1Δ，1Σ+，1Σ-。因 为这不是一个简并态，所以只能是1Σ+和1Σ-其中之一。这就需要进一 步知道+、-对称性。方法是看跃迁系数： 3 -> 6 -.28038 4 -> 7 -.28038 由于正负号相同，所以是1Σ+态。那么第四激发态必然是1Σ-了。 另外还可以看第二个组态耦合出来的光谱项，只有1Σ+，因此第三 个激发态必然是1Σ+。两个结果一致。 [ Last edited by zzgyb on 2008-1-8 at 14:17 ] |
2楼2008-01-08 14:16:37
3楼2008-01-08 23:21:33
6楼2008-02-03 10:25:29
简单回复
小红豆14楼
2008-04-05 19:39
回复
