自旋污染的基本概念以及对计算结果的影响

自旋污染的基本概念以及对计算结果的影响

关于 Hartree-Fock计算（通常采用Rootaan的自洽场（SCF）方式）介绍性的描述集中在单重态体系，在这样体系中所有电子自旋都已经配对。通过假设计算是限制在每个占据轨道有两个电子计算，计算可以相对地容易进行。这通常是指自旋限制的Hartree-Fock计算或RHF。

对于自旋多重度不为一的体系，是不可能象自旋多重度为一的体系那样用RHF方式的。 在一个非限制性的计算中，有两组完全 的轨道轨道，一组是alpha电子的，一组是beta电子的。 通常这两组轨道用同样的基函数但是不同的分子轨道参数。

非限制性计算的优点是他们计算的效率非常高。 缺点是波函数不再是总自旋的本征函数，因此计算中可能会导入一些误差。

这种误差成为自旋污染。 自选污染如何影响计算结果：

自旋污染导致在我们看来是期望的波函数中混合进了一点其他的自旋态。这种情况下由于有了较多的变分自由度，有时候会导致轻微的降低计算的总能量。更多的情况是因为较高能量态的混入导致总能量的轻微升高。然而这些变化是不正确的波函数 的结果。因为这不是一种系统误差，不同态间的能量差可能受到相反的影响。高自旋的污染可以影响稽核和布局数分析，和显著地影响自旋密度。

作为自旋污染出现的一种检查，大多ab initio程序会打印出期望的总自旋值，如果没有自旋污染，它会等于s(s+1)，这里s 等于1/2的未成对电子数倍。由有机分子计算经验得到的一个简单规律是，当和s(s+1) 的差别小于10%时，自旋污染的影响可以忽略。虽然这个规则提供了一种快速检验的方法，通常用实证据或更严格的计算来检查你的结果是明智的。

自旋污染通常在非限制的 Hartree-Fock (UHF)和非限制性的Møller-Plesset (UMP2, UMP3, UMP4) 计算中出现。在DFT计算中很少会发现显著的自旋污染，甚至在采用非限制性的Kohn-Sham 轨道的情况。

非限制性计算通常结合一步自旋湮灭步骤，这个步骤从计算的波函数中除去很大百分比的自旋污染。它有助于最小化自旋污染但是无法完全阻止自旋污染。的最后值是对出现的自旋污染的最好的检查。 在Gaussian中，选项\告诉程序用湮灭后的波函数产生布局数分析。我不知道是否有程序在执行几何优化过程中使用湮灭的波函

数。限制性开壳成计算有可能进行自旋限制性开壳层计算（ROHF）。这种方法的好处是没有自旋污染。缺点是会需要额外的CPU时间在对单占据的和双占据的轨道都进行正确处理以及它们之间的相互作用上。作为采用的数学算法的一个结果，ROHF计算给出好的总能量和波函数但是单占据的轨道不严格服从Koopman定理。

当自旋污染导致的误差不可接受时，先执行开壳层计算是给出可靠波函数的最好方法。

在Gaussian程序中，限制性开壳层计算可以执行Hartree-Fock，密度泛函理论，MP2 和一些半经验波函数计算。但是ROMP2方 式不知吃解析梯度，因此最快的方式是在其他方式优化的几何上进行单点能计算。如果几何优化一定要在这样的理论水平下进 行，必须采用一种非基于梯度的方式，如Fletcher-Powell优化。（注意在G94手册中对所有情况都还没有这种功能）

自旋投影方式

另一种近似是先运行一种非限制性计算然后在得到的波函数上投影掉自旋污染PUHF，PUMP2）。

自旋投影结果并不给出用限制性开壳层计算得到的能量。 这是因为非限制性的轨道被优化来描述被污染的自旋态而不是被优化来描述自旋投影的态。

用自旋强制的UHF方式（SUHF）可以得到类似的效果。在这种方式中，通过一个拉格朗日乘因子法，UHF波函数中的自旋污染误 差被限制。当乘的因子趋于无穷时，这种方法完全地除去了自旋污染。实际上小的正值就可以除去大多数的自旋污染。

半电子近似对自由基计算，半经验程序常用半电子近似。半电子方式是RHF计算中处理一个单占据轨道的数学技术。它在拥有一个近似的 波函数和轨道能量的情况下能得到一致性的总能量。因为采用的是单行列式计算，因此没有自旋污染。

一致性的总能量使得这种方法可以计算单态-三态能隙：对单重态用RHF，对三重态用半电子近似。半电子近似计算不服从 Koopman定理。也不能得到自旋密度。Mulliken布居分析通常是相当合理的。

结合到具体的计算方法，使用 CASSCF，CASPT2，MRMP2，MRCI等multi-reference方法可以比较有效的消除spin contamination的问题。根据个人体会使用PUMP2(projected UMP2)似乎也是不错的选择（在Gaussian中的实现为在route section直接使用UMP2（对开壳层体系）关键词，然后在结果中选用projected后的能量即可），而且速度比基于

multi- reference的方法更快（CASSCF除外）。ab initio的计算主要基于error cancelling 比如反应热，反应势垒等，没有人对绝 对能量感兴趣。anyway，虽然UHF有spin contamination的问题，但是使用UHF计算的ground state的能量一般比ROHF更低，因为前者中LCAO的自由度更多。

补充信息

一些讨论和结果在：

W. J. Hehre, L. Radom, P. v.R. Schleyer, J. A. Pople \Wiley (1986)

An article that compares unrestricted, restricted and projected results is M. W. Wong, L. Radom J. Phys. Chem. 99, 8582 (1995)

Some specific examples and a discussion of the half electron method are given in T. Clark \

A more mathematical treatment can be found in the paper

J. S. Andrews, D. Jayatilake, R. G. A. Bone, N. C. Handy, R. D. Amos Chem. Phys. Lett. 183, 423 (1991)

SUHF results are examined in

P. K. Nandi, T. Kar, A. B. Sannigrahi Journal of Molecular Structure (Theochem) 362, 69 (1996)

谁能解释一下什么是自旋污染？ 就我的理解：

自旋污染表现在计算结果中操作（2指平方）不等于其本征值：s(s+1) 例如有一个单电子的体系，=0.75。但是计算结果可能要大于0.75，若比0.75大1

0％，那么计算结果就不可靠了。这也是对结果的影响。

只有开壳层才有自旋污染，闭壳层没有此问题。直观上来看，造成自旋污染的原因是 混入高自旋的态，所以总自旋变大。从计算原理上说，开壳层的a，b轨道能级不一样，会

自旋极化；若极化过多，则会产生自旋污染（呵呵，这话我也不太懂，不知道是否有错误

。哪位大虾指点一下）。

根据文献报导，UMP2是自旋污染最大的。若有未成对电子，需要用PMP2（即投影MP

2的能量）。其他的，UHF,UMP3,UMP4的自旋污染比较大。

要避免的话，HF就选择ROHF方法。不过此法有不少缺点，目前基本不用。第二个方法

就是采用自旋投影方法，gaussian基本使用此法。（我没用过，哪位用过说说怎么用的） 。第三个方法就是半电子近似的方法，这个用于半经验方法，象半经验程序mopac中有。第

四个办法，就是用casscf方法了，它没有自旋污染问题。还有......应该还有其他办法： ）

其实我没有处理过这个问题，但是我对此也感兴趣。不知我说的对不对，与大家讨论 讨论。

自旋污染看“before annihilation”一行

MP2能量看“EUMP2”一行有ROMP2算法的，可能能解决自旋污染问题。LZ可以试试看。

怎样消除自旋污染

Re：QC强制收敛-一般是没有办法了再强制，所用时间很长。 conconver=8是默认 可以调整=5 降低收敛标准

IOP里面有IOP（5/13=1）就是防止自选污染 还可以算下去

Re：用DFT方法进行优化自旋污染会小一些。 个人意见，

仅供参考

Re：我用的就是DFT 方法。但是也挺严重的。因为是高自旋，没有办法。 Re：如果是高自旋的体系是不是用MCSCF更为合理一些。 CASSCF优化， CASPT2校正能量。

Re：RO方法没有自旋污染，如RO-DFT（ROB3LYP，ROOLYP...），ROMP2，...但是可能会不收敛

guass03出错信息及解决办法整理（待续） 收敛问题的调整

如果SCF计算收敛失败，可以采用以下技巧：

1. 在Guess关键字中使用Core，Huckel或Mix选项，试验不同的初始猜测。 2. 对开壳层体系，尝试收敛到同一分子的闭壳层离子，接下来用作开壳层计算的初始猜测。添加电子可以给出更合理的虚轨道，但是作为普遍的经验规则，阳离子比阴离子更容易收敛。选项Guess=Read定义初始猜测从Gaussian计算生成的checkpoint文件中读取。

3. 另一个初始猜测方法是首先用小基组进行计算，由前一个波函得到用于大基组计算的初始猜测(Guess=Read自动进行)。 4. 尝试能级移动(SCF=Vshift)。

5. 如果接近SCF但未达到，收敛标准就会放松或者忽略收敛标准。这通常用于不是在初始猜测而是在平衡结构收敛的几何优化。SCF=Sleazy 放松收敛标准，Conver选项给出更多的控制。

6. 一些程序通过减小积分精度加速SCF。对于使用弥散函数，长程作用或者低能量激发态的体系，必须使用高积分精度：SCF=NoVarAcc。

7. 尝试改变结构。首先略微减小键长，接下来略微增加键长，接下来再对结构作一点改变。

8. 考虑使用不同的基组。

9. 考虑使用不同理论级别的计算。这并不总是实用的，但除此之外，增加迭代数量

总是使得计算时间和使用更高理论级别差不多。

10. 关闭DIIS外推(SCF=NoDIIS)。同时进行更多的迭代( SCF=(MaxCycle=N) )。 11. 更多的SCF迭代( SCF(MaxCycle=N)，其中N是迭代数)。这很少有帮助，但值得一试。

12. 使用强制的收敛方法。SCF=QC通常最佳，但在极少数情况下SCF=DM 更快。不要忘记给计算额外增加一千个左右的迭代。应当测试这个方法获得的波函，保证它最小，并且正好不是稳定点(使用Stable关键字)。

13. 试着改用DIIS之外其它方法(SCF=SD或SCF=SSD)。

自旋污染

自旋污染是在非限制从头计算中的误差，通常在总能量中轻微出现。这是由于基态波函和更高能态发生混合。高的自旋污染会影响结构和布居数分析，并极大地影响自旋密度。它还会影响收敛，特别是MPn任务。可以通过比较总自旋S的预期值和计算值验证自旋污染。S应当等于s(s+1)。对于有机分子，偏差不超过正确值10%的计算通常是可接受的。 未成对电子数 s 预期S值 0 0 0 1 0.5 0.75 2 1.0 2.00 3 1.5 3.75 4 2.0 6.00 5 2.5 8.75 。

自旋污染经常在非限制Hartree-Fock (UHF)计算以及非限制M?ller-Plesset (UMP2, UMP3, UMP4)计算中遇到。而在DFT，CI 或CC计算中几乎见不到任何显著的自旋污染。在包含金属原子的体系中和在过渡结构中，出现高自旋污染是很平常的。

非限制计算经常与自旋猝灭步骤结合，它能从某些点计算的波函中删除相当比例的自旋污染。这使得自旋污染最小，但不能彻底防止。在Gaussian中，选项\通知程序使用猝灭波函产生布居数分析。

可以进行自旋限制开壳层计算(ROHF, ROMP2)。它的优势是没有自旋污染。缺点是需要额外的CPU时间，用于正确控制单占据轨道和双占据轨道以及它们的相互作用。 这个数学方法的结果是，ROHF计算给出很好的总能量和波函，但是单占据轨道能量不严格服从 Koopman理论。前缀RO可以用在HF，所有的DFT方法，AM1，MINDO3，MNDO，和MP2波函中。ROMP2方法尚不支持分析梯度。

另一个方法是运行非限制计算，在获得波函(即Gaussian输出文件中打印的 PUHF和PMP2值)之后，接着投影除去自旋污染。这给出的是对能量而不是对波函的修正。

自旋投影结果不会给出使用限制开壳层计算得到的能量。这是因为优化的非限制轨道描述污染态而不是自旋投影态。

半经验程序经常对原子团计算使用半电子方法。半电子方法是一种数学技巧，用于处理RHF计算中的单占据轨道。这个结果依靠近似波函和轨道能量，与总能量一致。由于使用了单重行列式计算，所以没有自旋污染。总能量的一致使得它可以对单重态用RHF，对三重态用半电子计算单重-三重裂距。半电子计算不服从Koopman定理。同样，也不能得到自旋密度。而Mulliken布居数分析通常还算是合理的。

如果自选污染小，对自旋猝灭波函和自旋投影能量更适宜使用非限制方法。不要对DFT使用自旋投影。其中自旋污染的总和对使用限制开壳层方法更加重要。

参见： J. Baker, A. Scheiner, J. Andzelm, Chem. Phys. Lett. 216, 380 (1993). M. W. Wong, L. Radom, J. Phys. Chem. 99, 8582 (1995). H. B. Schlegel, \Computational Chemistry\使用赝势

从头计算不能仅仅忽略核电子，但它们可以被哈密顿量中的势能项代替。这称为核势，有效核势(ECP)或相对论有效核势(RECP)。核势必须连同相应的价电子基组一起使用。这不仅减少了计算时间，核势还可以包含相对论质量亏损和对重原子近核区非常重要的自旋耦合项的影响。这对Rb和更高的重元素是首选。 有些使用赝势的计算结果不同于全电子计算的结果，输出的总能量是价电子能量。这个能量的大小依赖于核势模拟多少个电子壳层，其中电子明确描述为价电子和价电子-1壳层，计算一般更准确。同样，维里定理的检查是没有意义的。

怎么解决S2的问题

# freq ub3lyp/gen scf=(maxcyc=120) geom=check guess=read 这个 的计算命令

但是最后算出 HF=-316.6140566 S2=1.702056 按道理这个S2=0.75 2重态的？

想请教下这个怎么解决？ 这个算的东西还可靠，吗 ？ 如何消？

高斯使用手册 没有这方面的 提示 就是什么投影方法 ？ 这个有谁用过吗 ？

能不能给点实例 介绍 学习！！下 ！！ 盼高手解决！！！

问题是这样的 我算的这个体系 作为初始步 4重态 能算出来 而且 结果很好 但是2重态 却这么也算不出来 ，好不容易计算个S2=1.70几 自旋污染严重 ， 我看别人文献 里面说 他们说这种结构不存在？！！！ 我是该放弃寻找还是 消除自旋污染 继续找？ 迷茫中

？ 望高手解决中！

我觉得自选污染显然就是四重态混杂引起的

你算出来的2重态能量是不是明显比四重态高，如果是，似乎没有什么必要进一步研究了 现在是2重态能量 比四重态高 ，问题就是四重态于2重态在这里会有交叉 这个时候怎么办？ 用QC IOP命令强制算出来的东西 都是自旋污染的（S2=1.7）的 这个时候我怎么办？

扫描这方法我也想了下 但是觉得不实际 而且太费时 我该怎么办？ 这个点算出来 该不改放弃 ？ 用rob3lyp

可能不如rdft/udft收敛快，甚至根本不收敛 方法已定 这个时候 有什么办法没找这个这个点 ？ 或者跟别人样说这个点不存在 换方法，用bp86

换方法 能解决自旋污染？

如果是自选污染 我可不可以认为这里有自旋交叉？ 放弃寻找这个点？

消除自旋污染只能用多参考态方法

Density Functional Theory

In density functional theory (DFT) the energy of a system is given as a sum of six components:

EDFT = ENN + ET + Ev + Ecoul + Eexch + Ecorr

The definitions for the nuclear-nuclear repulsion ENN, the nuclear-electron attraction Ev, and the classical electron-electron Coulomb repulsion Ecoul energies are the same as those used in Hartree-Fock theory. The kinetic energy of the electrons ET as well as the non-classical electron-electron exchange energies Eexch are, however, different from those used in Hartree-Fock theory. The last term Ecorr describes the correlated movement of electrons of different spin and is not accounted for in Hartree-Fock theory. Due to these differences, the exchange energies calculated exactly in Hartree-Fock theory cannot be used in density functional theory.

Various approaches exist to calculate the exchange and correlation energy terms in DFT methods. These approaches differ in using either only the electron density (local methods) or the electron density as well as its gradients (gradient corrected methods or generalized gradient approximation, GGA). Aside from these \DFT methods, another group of hybrid functionals exists, in which mixtures of DFT and Hartree-Fock exchange energies are used.

1) Local methods

The only local exchange functional available in Gaussian is the Slater functional. Combination with the local VWN correlation functional by Vosko, Wilk and Nusair gives the Local Density Approximation (LDA) method. For open shell systems, using unrestricted wavefunctions, this is also referred to as the Local Spin Density Approximation (LSDA). This method is used in Gaussian with either the LSDA or SVWN keyword. A sample input for the frequency calculation for the methanol-1-yl radical using the SVWN functional and the 6-31G(d) basis set is as follows:

#SVWN/6-31G(d) freq

SVWN/6-31G(d) freq methanol-1-yl radical 0 2 C1

O2 1 r2

H3 2 r3 1 a3

H4 1 r4 2 a4 3 d4 H5 1 r5 2 a5 4 d5

r2=1.35454381 r3=0.97675734 r4=1.09713107 r5=1.09208721 a3=108.95266528

a4=119.30335589 a5=113.0329959 d4=-26.58949552 d5=-149.77421707

Please observe that the acrynom VWN is used to indicate different functionals in different software packages. While in Gaussian this keyword refers to functional III described in the VWN-paper (S. J. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 1980, 58, 1200), most other packages refer to functional V in the same paper. This latter functional can also be used in Gaussian with the VWN5 keyword.

2) Gradient corrected methods

The gradient-corrected exchange functionals available in Gaussian include (given with their abbreviations in parenthesis):

? ? ? ? ?

Becke88 (B)

Perdew-Wang (PW91)

Modified Perdew-Wang by Barone and Adamo (MPW) Gill96 (G96)

Perdew-Burke-Ernzerhof (PBE)

and the gradient-corrected correlation functionals are

? ? ? ? ?

LYP by Lee, Yang, and Parr (LYP) Perdew-Wang (PW91) Perdew 86 (P86) Becke96 (B96)

Perdew-Burke-Ernzerhof (PBE)

In all cases, the names of these functionals refer to their respective authors and the year of publication. All combinations of exchange and correlation functionals are possible, the keywords being composed of the acronyms for the two functionals. The frequently used BLYP method, for example, combines Becke's 1988 exchange functional with the correlation functional by Lee, Yang, and Parr. An input file for the BLYP frequency calculation of the methanol-1-yl radical is:

#P BLYP/6-31G(d) freq

BLYP/6-31G(d) opt min methanol radical

0 2

C1

O2 1 r2

H3 2 r3 1 a3

H4 1 r4 2 a4 3 d4

H5 1 r5 2 a5 4 d5

r2=1.38575119 r3=0.98034945 r4=1.09658704 r5=1.09123991 a3=108.13230593 a4=118.45265645 a5=112.14316572 d4=-29.40173022 d5=-145.83450795

Another frequently used GGA functional is BP86 composed of the Becke 1988 exchange functional and the Perdew 86 correlation functional. The PW91 functional combines exchange and correlation functionals developed by the same authors in 1991. The keywords used in Gaussian for a particular GGA functional are combinations of the acronyms for exchange and correlation functionals. Use of the PW91 exchange and correlation functionals thus requires the keyword PW91PW91.

3) Hybrid functionals

The basic idea behind the hybrid functionals is to mix exchange energies calculated in an exact (Hartree-Fock-like) manner with those obtained from DFT methods in order to improve performance. Frequently used methods are:

?

Becke-Half-and-Half-LYP (BHandHLYP) uses a 1:1 mixture of DFT and exact exchange energies: EXC = 0.5*EX(HF) + 0.5*EX(B88) + EC(LYP)

Becke-3-LYP (B3LYP) uses a different mixing scheme involving three mixing parameters: EXC = 0.2*EX(HF) + 0.8*EX(LSDA) + 0.72*DEX(B88) + 0.81*EC(LYP) + 0.19*EC(VWN)

In this latter case, the B88 gradient correction to the local LSDA exchange energies carries its own scaling factor of 0.72. The three scaling factors have been derived through fitting the parameters to a set of thermochemical data (G1 set).

PBE0 uses a 1:3 mixture of DFT and exact exchange energies: EXC = 0.25*EX(HF) + 0.75*EX(PBE) + EC(mPW91)

B98 (Becke98) is based on a 10-parameter equation scaling components of exact exchange, GGA exchange and GGA correlation. All parameters have been optimized simultaneously to fit thermochemica data collected in the extended G2 data set: EXC = 0.2198*EX(HF) + EX + EC

?

?

?

The Becke-half-and-half-LYP method is used with the BHandHLYP keyword, the very popular Becke-3-LYP method is used with the B3LYP keyword, and the PBE0 method with the PBE1PBE keyword. A sample input file for the B3LYP frequency calculation on the methanol-1-yl radical is:

#P B3LYP/6-31G(d) freq

Becke-3-LYP/6-31G(d) opt min methanol radical 0 2 C1

O2 1 r2

H3 2 r3 1 a3

H4 1 r4 2 a4 3 d4 H5 1 r5 2 a5 4 d5

r2=1.37007444 r3=0.96903862 r4=1.08886325 r5=1.08381987 a3=108.87657926 a4=118.49962766 a5=112.61039767 d4=-29.19863769 d5=-146.73450089

It is important to note that implementations of the B3LYP and BHLYP functionals in different programs (Gaussian, Turbomole, Jaguar, Q-Chem) vary somewhat. Other hybrid functionals available in Gaussian are: B3P86 uses the P86 correlation functional instead of LYP, but retains the three parameters derived for B3LYP; B3PW91 uses the PW91 correlation functional instead of LYP, but retains the three parameters derived for B3LYP; B1B96, a one parameter functional developed by Becke using the Becke 96 correlation functional; MPW1PW91 developed by Barone and Adamo using a modified version of the PW91 exchange functional in combination with the original PW91 correlation functional and a mixing ratio of exact and DFT exchange of 0.25 : 0.75.

4) User-defined models

Calculation of reaction and activation energies in either radical reactions or transition metal mediated reactions often shows that the mixing ratio of HF- and DFT-exchange is not optimal for a given purpose. The definition of user-defined models tailored to a certain purpose is therefore sometimes necessary. This can be achieved by first choosing gradient corrected exchange and correlation functionals from the list of available choices and then specifying the mixing ratios of the various exchange and correlation energy terms. For this latter step Gaussian uses the general formula:

EXC = P2EXHF + P1(P4EXSlater + P3dEXnon-local) + P6EClocal + P5dECnon-local

The values for the parameters P1 - P6 are provided through the IOP switches IOP(5/45=xxxxyyyy), IOP(5/46=xxxxyyyy), and IOP(5/47=xxxxyyyy) with P1 = xxxx/1000 from IOP(5/45) etc. The use of this facility will be illustrated with the MPW1K functional by Truhlar and coworkers (B. J. Lynch, P. L. Fast, M.

Harris, D. G. Truhlar, J. Phys. Chem. A 2000, 104, 4811) developed to optimize reaction and activation energies of free radical reactions. This method is very similar to the MPW1PW91 hybrid functional, but uses a mixing ratio of exact and DFT exchange energy terms of 0.428 : 0.572. An input file for the methanol-1-yl radical frequency calculation at the MPW1K/6-31G(d) level is:

#P MPWPW91/6-31G(d) freq

IOP(5/45=10000428) IOP(5/46=05720572) IOP(5/47=10001000)

MPW1K/6-31G(d) opt min methanol radical 0 2 C1

O2 1 r2

H3 2 r3 1 a3

H4 1 r4 2 a4 3 d4

H5 1 r5 2 a5 4 d5

r2=1.35355619 r3=0.95690342 r4=1.08170148 r5=1.07685321 a3=109.38001654 a4=118.56662303 a5=112.99211737 d4=-29.64928765 d5=-147.42512309

The gradient corrected MPWPW91 functional is chosen here as the basis of the hybrid method and the amount of exact exchange is specified through parameter P2. The composition of the functional is also described in the output file in a compact format:

IExCor= 908 DFT=T Ex=mPW+HF Corr=PW91 ScaHFX= 0.4280 ScaDFX= 0.5720 0.5720 1.0000 1.0000

The value given after ScaHFX corresponds to P2 while the four values listed after ScaDFX are P3 - P6. Please observe that theses lines are only printed when using the extended output option of Gaussian. Before using this facility with a particular version of Gaussian, please make sure that your IOP-directives are read correctly. In some versions of Gaussian, the options P1 - P6 are set with IOP(3/76) - IOP(3/78) instead of IOP(5/45) - IOP(5/47). In this latter case, the values of P1 - P6 are given in five-digit format. The keyword line for the PBE0 frequency calculation thus becomes:

#P PBEPBE/6-31G(d) freq

IOP(3/76=1000002500) IOP(3/77=0750007500) IOP(3/78=1000010000)

Absolute and relative CPU-times for DFT- and Hartree-Fock-calculations depend dramatically on the size of the system under investigation, the computer hardware and the operating system used. The following table contains CPU times (in seconds) for all frequency calculations on the methanol-1-yl radical described above with Gaussian 03, Rev. B.03. The platforms used are the CIP-room Pentium 4/2 GHz computers under LINUX. In all cases the frequency calculation has been performed on the structure optimized at the given theoretical level. The 6-31G(d) basis set has been used in all cases.

CPU (P4/2GHz) vib(1) vib(2) vib(3) [sec] [cm-1] [cm-1] [cm-1] 31 69 71 73 11 - 3690 1455 1036 3599 1463 1042 3761 1508 1072 3950 1552 1102 4123 1627 1156 3650 1459 1056 method USVWN/6-31G(d) UBLYP/6-31G(d) UB3LYP/6-31G(d) UMPW1K/6-31G(d) UHF/6-31G(d) obsvd.

Despite all platform-dependend variations, it is clear that hybrid DFT calculations are similarly laborious as GGA calculations, while local DFT methods are somewhat and HF-calculations much cheaper. In matching calculated (harmonic) gas phase vibrational frequencies with experimentally observed anharmonic values measured in a matrix (R. D. Johnson III, J. W. Hudjens, J. Phys. Chem., 1996, 100, 19874), please remember that uniform scaling is often performed to improve correspondence (see our earlier discussion of this point). The values reported here are unscaled. A general trend visible in this example is the increase of vibrational frequency wavenumbers with an increase in the amount of Hartree-Fock exchange contributions (varying between 0 and 100% from the second entry to the bottom in Table 1).

5) Further technical considerations

Evaluation of the exchange-correlation energies of all density functional methods implemented in Gaussian involves a grid-based numerical integration step. The computational effort required for this step strongly depends on the selected grid size. The larger the number of integration points per atom, the higher is the computational cost and the better the numerical accuracy of the calculation. Several predefined grids are available in Gaussian, which can be used through the keyword

int=(gridsize) or int=(Grid=gridsize)

Alternatively, the grid size can also be specified with IOP settings, using either IOP(5/44=gridsize) or IOP(3/75=gridsize), depending on the version of Gaussian. The following are acceptable values for gridsize:

corresponding (radial shells * CPU IOP setting angular points) [sec] 2 frequencies (cm-1) gridsize Etot CoarseGrid (35*110), pruned 33.8 -115.0520127 -110.9325 -0.0015 -0.0009 -0.0008 24.6825 59.6546 SG1Grid finegrid ultrafine mmm,nnn 1 4 5 mmmnnn (50*194), pruned 45.1 -115.0520371 (75*302), pruned 77.1 -115.0520326 (99*590), pruned 192.7 -115.0520324 (mmm*nnn), not pruned - - -30.5802 0.0015 0.0015 0.0016 10.5877 29.2740 -20.7294 -8.3913 0.0009 0.0009 0.0014 11.1426 -0.0013 -0.0013 0.0011 2.6708 4.0642 14.0308 -

The CPU times and total energies listed here (in sec.) together with the grid specifications are those for the UB3LYP/6-31G(d) frequency calculation on the methanol radical mentioned before on the CIP room Pentium 4 LINUX PCs. It is clear from these results that neither the total electronic energy nor the low-frequency vibrations converge quickly with increasing grid size. Thus, for each application it must be decided, whether the substantial computational effort for a very high grid size is justified. It is also important to note that the energy variation with grid size is significant enough that all calculations for stationary points on the same potential energy surface must be performed with the same grid in order to be comparable!

In addition to the grid selection for integral evaluation using the int=(gridsize) keyword, a particular grid can also be selected for the calculation of second derivatives using the CPHF=gridsize keyword. The available options for gridsize are identical in both cases. If no explicit grid is specified through the CPHF keyword, the following relations exist between the grids used for integral and second derivative calculations:

gridsize used for grid used for integral evaluation second derivatives CoarseGrid SG1Grid finegrid ultrafine (mmm*nnn) CoarseGrid CoarseGrid CoarseGrid SG1Grid (mmm*nnn)

These combinations are usually appropriate for most systems. The default grid size for integral evaluation in the most recent versions of Gaussian is finegrid. Better grids are required for kinetic isotope effect calculations, the optimization of stationary points on very flat potential energy surfaces, and IRC-following on flat potential energy surfaces. Increasing the CPHF grid size over what is specified per default through the integral grid selection is rarely helpful. Smaller grids can be quite useful to speed up geometry optimizations of large systems, especially when still far away from the next stationary point.

6) Treating open shell systems

For practical reasons the DFT methods also use a DFT \product of one electron functions (now termed Kohn-Sham orbitals). As in Hartree-Fock theory, electrons of unlike spin can either use the same spatial orbitals (restricted DFT methods) or different spatial orbitals (unrestricted DFT methods). Specification of a restricted open shell wavefunction for the Becke3LYP hybrid functional requires the ROB3LYP keyword, while the unrestricted version is specified using UB3LYP. If no explicit choice is made in the input file, open shell systems are treated with the unrestricted approach by default. It should be noted that unrestricted DFT methods (GGA or hybrid functionals alike) are much less affected by the phenomenon of spin contamination than Hartree-Fock theory and also give much better predictions of ESR spectra of organic radicals. The UB3LYP/6-31G(d) calculation of the allyl radical is achieved with the following input file:

#UBecke3LYP/6-31G(d) scf=(tight)

UB3LYP/6-31G(d) ally radical (C2v) 0 2 H

C,1,r2

C,2,r3,1,a3

C,2,r3,1,a3,3,180. H,3,r5,2,a5,1,0. H,4,r5,2,a5,1,0.

H,3,r7,2,a7,1,180. H,4,r7,2,a7,1,180.

r2=1.09054601 r3=1.38611149 r5=1.08492921 r7=1.08690606 a3=117.45043254 a5=121.64506508 a7=121.12827199

Comparison of the expectation value of the operator calculated at the UHF/6-31G(d) level (0.9729) with that obtained at UB3LYP/6-31G(d) (0.7818) shows that the spin contamination problem is much reduced with DFT methods. From the construction principle of the hybrid DFT methods it is clear that a larger HF-exchange component will also give rise to larger amounts of spin contamination. Using the allyl radical at the UB3LYP/6-31G(d) structure as an example, the value amounts to 0.7650 at the BLYP/6-31G(d) level and to 0.8228 at the BHandHLYP/6-31G(d) level. The precise values of the

expectation value of the operator will also depend to some extend on the basis set used.

Another important point to consider in DFT calculations concerns the self interaction error resulting from an artificial interaction of the electron density with itself. While the definition of the exchange and Coulomb energies in Hartree-Fock theory leads to a perfect cancellation of these self interaction terms, the cancellation is incomplete in DFT methods. As nicely explained by Koch and Holthausen in their monograph on DFT methods, this can be easily demonstrated using the hydrogen atom. The exact energy for this one-electron system is -0.50000 Hartree. As there is no second electron in the system, no correlation energies have to be calculated and Hartree-Fock theory is able to provide the exact answer in the complete basis set limit. Using a basis set considered \for molecular systems such as the cc-pVQZ basis set, a value of -0.499946 Hartree is indeed obtained. Due to the parameterized exchange and correlation functionals, however, much more negative values can be obtained in particular with hybrid DFT functionals. Using the same basis set, the UB3LYP functional predicts a value of -0.502346 Hartree!!

本文来源：https://www.bwwdw.com/article/pwu2.html