【转】高斯优化过渡态的经典总结
2012-10-16 10:51阅读:
【转】高斯优化过渡态的经典总结
一般地,优化所得驻点的性质(极小点还是过渡态)要靠频率来确定;而对过渡态,要确定反应路径(即到底是哪个反应的过渡态)必需要做IRC了,不然靠不住的(往往用QST找到的过渡态并不一定就是连接输入反应物和产物的过渡态)。
在我们用QST2或QST3来优化过渡态时,需输入反应物和产物,实际上反应物和产物的输入顺序是没有关系的。就是说,先输反应物后输产物和先输产物后输反应物得到的是同样的过渡态。这也好理解,QST2里对过渡态的初始猜测实际上是程序自动将输入的反应物和产物的各变量取个平均,所以输入顺序是没有关系的。对QST3和TS,需人为指定过渡态的初始猜测。
上面说的反应物和产物的输入顺序没有关系,有个前提条件,就是反应物和产物的自旋多重度一致。对QST3,因为是人为指定过渡态的初始猜测,所以没有影响;而对QST2,过渡态的自旋多重度默认和后面输入的一致。如果反应物和产物的自旋多重度一致,那随便先输哪个都没关系;而如果反应物和产物的自旋多重度不一致,这时该怎么办???
当反应物基态多重度为1,产物基态多重度为3时,到底将过渡态的多重度定为1还是3?还是两个都试,哪个能量低取哪个?假如取1,在下来做IRC
验证时,Reverse还好办,Forward却是按多重度为1做的,这样怎么能和多重度为3的产物连接起来?
要将其关联起来,需有一个为激发态,这样才能在自旋多重度上保持一致。可这种处理对吗?还有更好的思想吗?另外,我觉得IRC难用得很,很难控制,也不知自己钻进了哪条死胡同!!
一直对此很迷茫,希望各位大侠援助啊!
对于这点确实很迷茫,我觉得好像得用cas解决,有看过类似的文献,上面有用cas得出过反应物基态经激发态再回到基态得出产物的反映路径。我没有试过,对于cas我很是头疼,理论和实践都一无所知。
作 IRC 能说明哪些问题?
irc做完,就说明我们完成了这个反应路径地计算,
在SUMMARY OF REACTION PATH
FOLLOWING:,我们可以看到分子过渡态的键长键角与能量随着反应坐标的变化而变化,如果你将反应坐标与能量作图,就可以得到一条过渡态曲线,由于irc计算是结构沿着反应路径的方向,在每个点进行优化的,所以如果你找的过渡态是正确且优化是成功的话,TS确实是连接两个minimum的。
一般IRC的两个终点并不恰好就是反应物和生成物,除非不断加大Maxpoints。然而此时IRC又往往会出错,或者并不能达到指定的Maxpoints,例如,我想正反两个方向各算20各点,IRC计算正常完成,但结果中并没有41个构型。如果这时作能量-反应坐标曲线,如何确定反应物和生成物的反应坐标?
如果你的过渡态确实做得漂亮,那么就一定可以成功地由过渡态找到与之对应的反应物、生成物。可这只是理论上的东东,现实中,由于反应的势能面实在太复杂,你千辛万苦所得到的过渡态是很难天衣无缝地连接反应物、生成物的。不断加大Maxpoints
简直是Kill
time。
一般地,我们可以将IRC最后得到的两个结构的所有数据分别与你做QST时所用的两个结构进行对比,只要结构数据的差别、能量差别比较小,就可以基本认定(因为反应的势能面太复杂)你IRC是成功的了。
说得好!取点越多越好,但我发现往往作不到这一点,所以我也是将IRC最后得到的两个结构的数据分别生产物态和产物态的两个结构进行对比。
而且我发现,作IRC时频率计算给出的讯息很重要,呵呵!
过渡态寻找小结
刚刚做了一段时间的过渡态,期间碰到了许多的困难。寻找过渡态不是一件容易的事(对于我和大多数刚涉及量化的人来说),因此我希望通过写这个经验小结能对大家有些帮助。
1.首先遇到的问题是,用哪种方法来寻找过渡态?
GAUSSIAN提供的方法是QSTN和TSN方法。两种方法各有优点和缺点。QSTN方法特别QST3方法要求输入反应物,过渡态的猜测结构,产物这三者的结构。特别麻烦。但很管用,一般不会出现不收敛的情况。对于TSN(对应关键词为OPT=TS)方法,只要求输入过渡态的初始结构,但这个初始结构非常的关键,如果结构不好,则很容易出现不收敛的情况。所以我建议,如果是刚开始做过度态的话,用QSTN方法是好的选择,等有了“感觉”之后,再用TSN方法。
2. 怎么解决经常出现的错误?
在找过度态的时候,经常碰到的一些问题就是不收敛(1).,有一个错误的本征值(错误信息为:there
is a wrong sign eigenvalue
in hessian
matrix.....)(2),和LINK9999错误导致退出。(3)
对于不收敛的情况,可以分为两类,比如提示信息里的CONVERGENCE
FAILER
提醒收敛到了10(-5),而此时你设定的SCF循环次数也仅仅是64步,那么完全有希望通过加大SCF循环次数来达到收敛的目的。倘若只收敛到10(-3)或10(-2),此时加大循环次数可能就没用了。结果还是CONVERGE
FAILER。
此时可采用SCF=QC,来达到强制收敛的目的
。因为SCF=QC(LINK508)的计算量比默认的L502要大,所以不到万不得以就不用它了。
出现第二个错误可以直接用
关键词OPT=NOEIGEN
来实现。
LINK9999出错是因为已经走完了默认的步数,但还未完成。系统会自动跳出。出现这种情况大多数就是因为优化步数和SCF步数超过了默认值。可用OPT(MAXCYCLE=100)和SCF(MAXCYCLE=300)来改错。
3.怎么样控制过渡态的优化,使得过渡态不至于收敛到其他的分子结构中去?
我用GAUSS VIEW
可以解决这个问题
,当刚开始运行GAUSSIAN时,你用GVIEW去打开输出文件时,你可以看到你的过渡态的初始输入结构,当一个循环过后(从上一个LINK502到下一个LINK502),你再打开输出文件,你就可以清晰地看到优化一步后分子的构型,这样就可以随时监控过度态分子的结构,倘若已经有收敛到其他分子构型的趋势时,你就可以把它给KILL了,而不至于需要等全部工作结束后,打开输出文件才知道已经不是想要的过渡态了。
如果收敛到其它的构型上去,可以考虑缩小OPT的步长.iop(1/8=2或3)即可。
4.
还需要加其它的关键词吗?
建议在OPT中加入CALCFC。这样可以加大找到过渡态的几率。本人深有体会!
先写这么多了,难免有错误和不恰当之处,还希望大家来指点和补充!
我找过渡态的经验是:找过渡态的关键点,是如何根据反应体系的立体结构和反应过程化学键的断裂和形成的轨道的空间形状,提出一个合理初始的过渡态的结构,再进行优化。提出的初始过渡态结构越接近真实结构,就越容易找到,否则花一年时间也可能找不到一个过渡态。而找过渡态的QSTN方法、TSN方法的差别不是太大。
这几个帖子都很好,我也找了十几个简单的过渡态,不是直接找到的,而是先用目标过渡态的四元环结构做,然后对四元环结构逐步“修饰”直到那个四元环就是我需要的四元环为止,这个方法来得一般比较简单,但也有无法“修饰”到目标过渡态的情况,其他方法都是一样的。
10(-5)表示十的-5次方。
每个循环都会列出SCF计算的结果,和收敛值啊。可以找的到的。
iop(1/8)=2
相当于在找过度态的时候,以2A(A为基本单位长度)为单位来寻找过渡态。
就象有一百个抽屉里只放了一个苹果,当然是一个一个抽屉打开,找到苹果的把握最大了。此时意味着iop(1/8)=1。
有虚频说明结构处于不稳定状态。过渡态的初始结构猜测应该根据已有的知识和经验啊,哪种构型或构象最可能是过渡态,注意有时候猜测可能与计算结果完全相反。
反应坐标
势能面是研究化学反应历程的基础。对一个势能面来说,我们比较关心其关键点。势能面上的关键点是指势能的极值点,包括反应物
(R)、产物
(P)、中间体 (I)
与过渡态(TS)
或鞍点(saddle
point)等。连接反应物、过渡态和产物的反应途径,是一条能量最低的路径(minimum
energy
path),称为反应坐标(reaction
coordinate)。各关键点是通过反应坐标联结起来的,过渡态是其上面的一级鞍点。
*不同类型的反应坐标
我们常说的反应坐标、最小能量路径等用语,严格地说,是有一定区别的。
(1)最速下降路径:也称最小陡降路径。从过渡态出发,体系会沿着最陡的斜坡向反应物深谷和产物深谷移动,移动路径与等能量面正交,方向与势能梯度相反。这条路径的任何地方,其垂直方向有一极小点,就是垂直切割最速下降路径的交点。
(2)内禀反应坐标:内禀反应坐标就是连接反应物、过渡态和产物的最速下降路径。
(3)最小能量路径(MEP):通常选定一个反应坐标,使能量相对于其它坐标为最小。这样的反应路径,常常会引起弯曲甚至不连续,特别地与最速下降路径偏离甚远,而后者是真正的最小能量路径。举例来说,Relaxed
Scan
的结果就是一条MEP。
(4)经典轨迹:上述路径在势能面上对应着一组原子无限缓慢地移动(零动量),而轨迹是通过在包括核最初动量的合适初始条件下求解反应物的运动方程得到的。动力学效应,例如离心力,将使分子沿着不同于最速下降的路径移动。随着分子象台球一样从势能山脊上弹起又滑下山坡,轨迹可能十分繁长。
irc输入和结果
%Chk=irc
#p hf/sto-3g irc(forward,calcfc)
guess=read geom=allcheck
******************************************
%nproc=2
%Chk=irc2
#p hf/sto-3g irc(reverse,calcfc)
guess=read geom=allcheck scf(maxcyc=200)
(results:)
SUMMARY OF REACTION PATH
FOLLOWING:
与scan看结果一样是看OPti上的那个orient。
如果是单独作irc(forward)和irc(reverse),取出最后opt结果上面的那个orient来看irc到底是个什么样的过程。第一个是过渡态的orient了。
关于虚频的一个简单的理解
首先,什么是频率。
中学的时候我们学过简谐振动,对应的回复力是f=-kx,对应的能量曲线,是一个开口向上的二次函数E=kx^2/2.
这样的振动,对应的x=0的点是能量极小值点(简单情况下也就是最小值点)。这时的振动频率我们也会求:ω=2π
sqrt(k/m)。显然它是一个正的频率,也就是通常意义下的振动频率。
那么,一维情况下,如果能量曲线是一个开口向下的二次曲线呢?首先,从能量上看,这是个不稳定的点,中学的物理书上称为“不稳平衡”。用现在的观点看,就是这一点导数是零(受力为0),且是能量极大值。如果套用上面的公式,“回复力”f=-k'x(实际上已经不是回复,而是让x越来越远了),这里k'是个负数,ω=2π
sqrt(k'/m)显然就是一个虚数了,即所谓的虚频。Gaussian里面给出一个负的频率,就是对应这个虚频的。
实际情况下,分子的能量是一个高维的势能面,构型优化的时候,有时得到了极小值点,这样这个点的任意方向上,都可以近似为开口向上的二次函数,这样这里对应的振动频率就都是正的。对于极大值点,在每个方向都是开口向下的二次函数,那么频率就会都是负的——当然一般优化很少会遇到这样的情况。对于频率有正有负的情况,说明找到的点在某些方向上是极大值,有些方向上是极小值。如果要得到稳定的能量最低构型,显然需要通过微调分子的构型,消去所有的虚频。如何微调?要看虚频的振动方向。想象着虚频对应的就是开口向下的二次函数,显然,把分子坐标按照振动的方向移动一点点,分子应该就可以顺着势能面找到新的稳定点,但是也不能太小。而所谓的过渡态,则是连接反应物和产物之间的最低能量路径上的能量极大值。好比山谷中的A,B两点,它们之间的一个小土丘,就是过渡态,从A到B的反应,需要越过的是这个小土丘,而不是两边的高山。这样,过渡态就是在一个方向上是极大值,而在其它方向上都是极小值的点。因此,过渡态只有一个虚频。
频率分析只能在势能面的稳定点进行,这样,频率分析就必须在已经优化好的结构上进行。频率分析的另外一个用处是判断稳定点的本质。稳定点表述的是在势能面上力为零的点,它即可能是极小值,也可能是鞍点。极小值在势能面的各个方向都是极小的。而鞍点则是在某些方向上是极小的,但在某一个方向上是极大的,因为鞍点是连接两个极小值的点。有关鞍点的信息:
1.负的频率2.频率相应简正振动的模式
当一个结构产生负的振动频率时,可以表明在该振动方向可能存在着能量更低的结构。判断所得鞍点是不是需要的鞍点的方法,就是察看它的简正振动模式,分析是不是可以导向所需要的产物或反应物。进一步的,更好的办法是通过IRC计算来判断反应物,产物与得到的鞍点是否有关系。
向大家推荐一本书,翻译成中文的(原文我没有见过),可能是《化学反应中的电子》。
在寻找过渡态计算中opt=ts,
需要计算频率,如果有一个虚频就说明优化所得的构型是个
过渡态。有时在做部分限制优化opt=minimal时,也能得到一个虚频(一般发生在frozen
部分的化学键上)。我想问一下,虚频的数值大小本身具有什么物理意义。后一种情况是否也可以认定为过渡态。
我的理解: 虚频就是在这个振动方向上力常数是负的
过渡态之所以有一个虚频,是因为它是一个鞍点,鞍点上有一个负的梯度的方向分子沿这个方向振动时将转化为反应物或产物。就象从山上掉下来,受到的力可以认为是负的。但并不是有一个虚频就是过渡态,还要在这个虚频的振动方向上分别指向反应物和产物
Gaussian 程序频率计算中的几个问题
在Gaussian计算中,为了确定优化得到的几何结构是势能面上的局域极小点还是鞍点,或者要得到相关的热力学性质,经常需要对优化后的几何结构进行振动分析。这里我们将讨论几个频率计算中常见的一些问题。希望能对初学Gaussian的人有所帮助。
首先,原则上说,振动频率分析只对稳定结构有意义。这里所说的稳定结构包括是势能面上的局域极小点和鞍点。如下图1所示是一维自由度上的势能面,A和B处在势能面的局域极小点,而处在势能面的鞍点上。他们在都处在平衡位置(原子核受力为零),不同的是,A和B来说离开平衡位置会受到指向平衡位置处的力,而C离开平衡位置会受到远离平衡位置的力。
因此A和B处在稳定平衡点,C处在不稳定平衡点。实际上,一个分子可以有很多的自由度,如果在所有自由度上分子都处在稳定平衡,就是稳定的分子。频率分析得结果是所有频率都是正的,表明这是一个局域的极小点。如果分子只在一个自由度上处于不稳定平衡位置,其他自由度上都处在稳定平衡位置,说明该结构是一阶鞍点。分子在稳定自由度方向上的振动才是真实的振动,在不稳定自由度方向上的实际上是不会有振动的。不过我们可以对不稳定方向上的运动也按振动来做数学处理,会的到负的振动频率,我们称它为虚频。虚频的出现表明该结构为鞍点。
图1
势能面上的局域极小点和鞍点
第二,Gaussian计算中,频率的计算一定要在和分子结构优化相同的方法,基组下进行,否则计算的结果是没有意义的。我们知道,任何理论水平下的计算,都是在一定的近似下进行的,不同的理论水平的近似程度是不同的。在一种理论水平A下优化的稳定结构Geom_A会和另一种理论水平B下优化的稳定结构Geom_B有差别,也就是说Geom_A不会是理论水平B下的稳定结构。根据前面我们所讨论的,在理论水平B下对一个不稳定的结构进行频率分析是没有意义的。
图2示意说明了不同理论水平下稳定点结构的不同。
图2 不同理论水平下优化的稳定结构是不同的
第三,频率计算中可以考虑同位素效应(Freq=ReadIsotopes)。在波恩-奥本海默近似下,对于同一种元素采用不同的同位素对几何优化和电子结构计算没有影响,频率计算所需的力常数矩阵(Hessian矩阵)也不会变化,变化的只是约化质量。容易理解,重的同位素会导致低的振动频率。实际上,原子序数大的元素的同位素效应非常不明显,一般只需考虑H原子的同位素效应。
第四,各种方法计算的频率和实验结果之间存在系统误差,需要乘以一个约化因子来进行校正(Scale=f)。一般来说,理论计算的频率值会比实验结果大。下面是一些理论水平下的约化因子。注意频率和零点能的约化因子是可以不同的。更多水平下的约化因子需要查文献获得。
方法:
约化因子
约化因子
(频率)
(ZPE)
HF/3-21G
0.9085
0.9409
HF/6-31G(d)
0.8929
0.9135
MP2(Full)/6-31G(d)
0.9427
0.9646
MP2(FC)/6-31G(d)
0.9434
0.9676
BLYP/6-31G(d)
0.9940
1.0119
B3LYP/6-31G(d)
0.9613
0.9804
SVWN/6-31G(d)
0.9833
1.0079
第五,Gaussian的频率计算有时会遇到下面的警告:
Warning -- explicit consideration of
3 degrees of freedom
as vibrations may cause
significant error
这时一般有两种可能:一种可能是,优化的几何结构不够精确,还没有达到稳定点。对于这种情况,需要考虑用OPT=tight或OPT=VeryTight,结合Int=fine或者Int=VeryFine进行更加精确的优化。第二种可能是,在计算的振动模式中,有部分的低频模式对应于内转动模式,在热力学分析中应该按自由转子或受阻转子模型处理,如果按谐振子模型处理会造成较大的误差。采用受阻内转子模型,可以用Freq=Hindered-Rotor计算或者自己手动计算这些振动模式对热力学性质的贡献。
第六,对于计算出的力常数较小的振动模式,其势能面形状可能不满足谐振子模型的要求,力常数函数需要考虑非谐性的成分,这时可以用Freq=Anharmonic来进行非谐性处理。
第七,在频率输出前,会有这样的输出:
Low frequencies --- -0.0008
0.0003 0.0013 40.6275
59.3808 66.4408
Low frequencies --- 1799.1892
3809.4604 3943.3536
上面行的六个模式实际对应于平动(前三个)和转动(后三个),理论上,这六个数都应该是零。如果是对过渡态进行的频率计算,第一行会先输出虚频,然后才给出平动转动的数值。一般来说,平动的三个数都是接近于零的。如果转动的三个数不接近于零(一般解析频率10个波数以内,数值频率50个波数以内),说明需要进行OPT=Tight或OPT=VeryTight计算。第二行的几个实频率应该和随后输出中的相应频率进行对比,如果基本一样,说明频率计算结果很好;如果差别较大,说明这些频率受到的平动和转动的污染较大。
第八,频率计算后程序会进行一步优化计算,正常会得到四个判据都是YES的结果,这说明优化的结构很好。但是有时也会遇到四个判据不全是YES的情况,这时需要仔细分析,一般有三种情况:1、Force和RMS
Force都是收敛的,Displacement和RMS
Displacement虽然是NO,但是都比较接近收敛判据;
2、如果Force和RMS
Force都是收敛的,Displacement和RMS
Displacement不收敛且数值远大于收敛判据;3、Force和RMS
Force不收敛。对于第一种情况,我们可以不管它,仍然可以认为计算结果是可靠的。对于第二和第三种情况,一般说明优化过程中估算的Hessian是不准确的,优化的结构可能还没有达到稳定点,需要重新优化。如果反复优化仍然无法解决这个问题,建议在优化过程中使用OPT=CalcAll。
第九,优化一个局域极小点,收敛后从力常数本征值看应该是局域极小点(没有负的本征值),但是频率计算中会出现虚频和负的本征值。这种情况下虚频一般是由转动模式造成的,说明分子中两个基团之间的相互位置不是很合适,需要绕转动的键相对转动到一个合适的位置,重新优化。绕某个键转动两个基团,有时可以很方便地用修改二面角的方法实现:OPT=Modredundant结合分子描述后输入:*
m n *
[+=]value。其中m,n为连接两个基团的键的顶端原子。有时,用OPT=CalcAll也可以解决这个问题。
第十,Gaussian程序中,频率计算是和温度、压力无关的。因为Gaussian所考虑的是原子核在一定核构型和电子运动状态构成的势场内振动,一般来说,温度和压力的变化不会影响分子的结构和电子运动状态,所以振动的势能函数是与温度、压力无关的。Freq=Isotopes关键词要求输入的温度和压力影响的只是平动,转动和振动的平衡分布和它们对热力学性质的贡献。
分子激发态计算的理论方法一
1. 激发态计算方法的分类
按照化学模型,目前计算激发态的方法可以分为:
____DFT (TDDFT,
MC-DFT,DFT+CI...)
Chemical____|
Model |____Wave
function____|~~~~Ab initio (MR-CI,
CIS...)
|____Semiempirical (Zindo/1,
INDO/S,CNDO/S...)
这与基态计算方法的分类是一样的,从中看不出激发态计算的特点。另一种分类是按照获得激发能(或激发态能量)的方式,可以分为两类:i)
电子组态方法。选择不同的电子组态的线性组合进行迭代,得到激发
态的绝对能量(有些程序经过简单的后期处理也能给出相对能量)。
分为多参考(MR)方法和单参考(SR)方法两种。常见的多参考方法有:MCSCF(有CASSCF,RASSCF等多种形式),MR-CISD,CASPT2等。标准的
MR的计算过程一般是三步:HF—〉MCSCF(常用CASSCF)—〉MRCISD或CASPT2...R是MR的特例,不需要中间的MCSCF步骤,计算中只需要对MR程序模块指定一个参考组态即可,计算量大大减少,但是精度也会下降。对于一些无法用SR很好描述的体系(如镧系化合物),计算是失败的,这时就必须用MR方法。ii)
线性响应理论(LRT)。基态用普通的理论方法求得,然后通过求解特殊
的线性响应方程,获得激发态相对于基态的垂直激发能(相对能量)。
有的方法求解方程需要做迭代,有的则不需要。
这一方法要求作为参考态的基态必须被很好地描述。常见的有TDDFT,
CIS,CCSD-LRT等,共同点是基态都使用SR模型。基态用MR模型的有MCLR
(MCSCF-LRT),MR-AQCC-LRT等方法,用于计算那些无法用SR模型很好描述基态的体系。
2. 一些方法的特点及软件
i)
Zindo/1,INDO/S,CNDO/S
都属于半经验方法,可用于研究原子数在100左右的大分子体系。软件:Zindo/1
Gaussian,
Zindo。Zindo结合ZINDO-MN模块可以用Zindo/1方法计算溶剂对激发态的影响。INDO/S,CNDO/S
商业MOPAC的MOS-F模块。免费的MOPAC
6.0/7.0无此功能
ii) CIS
基态能量用Hartree-Fock方法得到,激发态精度也是HF级别的,得到半定量的结果。适于研究原子数在50左右的中等分子体系。此外还有精度
更高一些的CIS(D),计算量当然更大。Gaussian
03的CIS还能考虑溶剂
环境对激发态影响。程序化的这些方法一般只能计算单重态和三重态。扩展CIS
(XCIS)方法
还能计算双重态和四重态。
软件:CIS Gaussian,
Q-Chem,CaChe,商业的MOPAC,TURBOMOLE,最新版MWChemCIS(D)
Q-Chem,TURBOMOLEXCIS
Q-Chem
iii)TD-DFT和TD-HF
目前比较流行的激发态计算方法,适用于十几个原子构成的中小分子,
小分子计算也能得到令人满意的结果。TD-DFT目前的缺陷是无法研究多
电子激发和Rydberg态。TD-HF早期也叫RPA,是TD-DFT的特殊形式,精度
较差。
Gaussian
03中TDHF/TDDFT的激发态计算还能考虑溶剂环境影响。
程序化的TD-DFT一般只提供对单重态和三重态计算的正式支持。最新版
NMChem的TD-DFT除了计算闭壳层分子的单重态和三重态,还能计算开壳
层分子的双重态和四重态。
软件:
TD-HF Gaussian, ADF, NWChem,
Q-Chem,DeMon,GAMESS-US,
ACES
II,TURBOMOLE
TD-DFT Gaussian, ADF,
NWChem,
Q-Chem。DeMon,TURBOMOLE,
CADPAC,BDF。
一些计算晶体的程序也可以做分子的TD-DFT计算,但是由于使用平面波基组,结果不太好。
iv) CC-LRT
CC-LRT包括CCS-LRT
(等价于CIS),CC2-LRT
(对CCSD-LRT的近似),CCSD-
LRT,CC3-LRT
(对CCSDT-LRT的近似)等不同激发级别的形式。最常用的是
CCSD-LRT,有很多种形式,如EOM-CCSD,VOOD-CCSD等。适用于中小分子
和小分子的计算,结果与TD-DFT相当,可能还略好一些。
软件:
CC-LRT
DALTON。CCSD-LRT的其他形式可在ACES
II,Q-Chem,PSI,
MOLPRO中找到。TURBOMOLE可做CC2-LRT计算。
v) SAC-CI。
精度略高于TD-DFT的方法,适于研究十几个原子构成的有机分子,还可
以包含多电子激发的影响。
软件:SAC-CI Gaussian 03
vi) DFT+CI混合方法。
结合了DFT和CI的优点,适用于中小分子激发态计算。
软件: DFT+CI TURBOMOLE
vii) MCSCF
(有CASSCF,RASSCF等)。
做(电子组态方法)多参考计算的基础。适用于中小分子激发态计算。
常用电子组态形式的CASSCF。线性响应理论形式(MCLR)用的较少。
软件:
MCSCF GAMESS-US, GAMESS-UK,
Gaussian, DALTON, MOLPRO,
MOLCAS,
COLUMBUS,NWChem
MCLR
DALTON,GAMESS-UK
viii)MR-CI及其变体
MR-CISD是计算(几个原子构成的)小分子激发态的标准方法。由于这
一方法存在着大小一致性问题,因而又发展了很多变体。比较简单的是
对MR-CISD做事后修正(+Q),主要有Davidson修正和Pople修正两种。复杂一些的有MR-ACPF、MR-AQCC、MR-CEPA等,其中表现最好的是MR-AQCC,以及后来发展的虚轨道修正MR-AQCC
(MR-AQCC-V)。有的程序还允许对MR-CI做激发级别的修正,得到MR-CIS,MR-CISD(T)
等。
软件: MR-CISD GAMESS-US,
GAMESS-UK, DALTON, MOLPRO, MOLCAS,
DIRAC,
MOLFDIR,COLUMBUS
变体: MOLPRO,
MOLCAS,COLUMBUS。COLUMBUS还能做MR-AQCC-LRT
计算,此方法综合了MR-AQCC与线性响应理论的特点。
MOLPRO对MR-CI及其变体使用了内部收缩,速度更快。
ix) MR-MPn
多参考微扰理论。目前主要有二阶微扰和三阶微扰。在程序作了优化的前提下,它们比相同级别的MR-CI计算速度快(例如在MOLPRO中,CASPT2比MR-CISD快),但是精度略低。目前应用的难题是入侵态问题。有些程序通过与MR-CI混合可以部分地解决这一问题。
软件: CASPT2 MOLPRO,
MOLCAS,COLUMBUS,GAMESS-UK。GAMESS-US的MC-QDPT方法和Gaussian的CASSCF+MP2也属于MR-MP2。
CASPT3
MOLPRO,GAMESS-UK
x) MR-CC和FCI
MR-CI的理论极限是完全组态相互作用(FCI),常用的MR-CISD是把FCI在
CI空间的展开截断到二级的形式。
鉴于耦合簇CC能获得比CI更好的结果,而且不会有大小一致性的问题,如果把FCI展开到CC空间并在一定的级别截断,就得到MR-CC。
由于FCI计算量相当大,所以目前也只能计算较轻原子构成的双原子分子。从理论上讲,MR-CCSD是小分子激发态计算最实用、最精确的方法。但目
前还有入侵态的问题没有解决。公开发布的程序中目前还没有此方法。
xi) MC-DFT,CAS-DFT:
多参考的DFT,尚在研究中。MOLCAS声称将在6.0版包含此方法。
xii) MR-Gn
多参考的Gaussian-n理论,可想而知是相当昂贵的方法。目前只有GAMESS-
US小组的一篇文章研究过这一理论。将来能否用于激发态计算并获得高精度结果还未知。
3. 方法的选择
构成分子的原子数:
2.....10.......20........50..........100..........
|_____________________|
|______________|
半经验
CIS
|_____________________|
MCSCF(注1)
|______________|
TD-DFT(注2),
CCSD-LRT
|_____|
MR-CI及其变体,MR-MPn(注2)
|
FCI
可见,越往下精度越高,适用分子越小。
注1:CASSCF的计算量和活性空间大小2A和活性电子数E有关。
C(2A,E)=(2A)!/[(2A-E)!*E!]。
C(2A,E)越大则计算量越大。RASSCF由于对组态进行了筛选,能适当减少计算量。对于小分子,MCSCF一般只作为高级MR计算(MR-CI,MR-MPn)的中间步骤。
注2:对于含有重元素的小分子激发态计算,还要考虑自旋轨道耦合(SO)造成的电子态分裂。例如AtBr的3Pi会分裂为0-,0+,1,2四个态。能够在激发态计
算中包含SO的程序有MR-CISD及其变体
+ SO:
COLUMBUS,MOLPRO,GAMESS-US
MR-MPn + SO:
MOLPRO,GAMESS-US
MCSCF + SO:
COLUMBUS,MOLPRO,GAMESS-US,
MOLCAS,TD-DFT + SO:
ADF
2002,尚在测试中,目前只能对原子获得较好的激发能
完全相对论MR-CISD:
DIRAC, MOLFDIR
完全相对论TD-DFT:
BDF,尚在测试中
如果是用Gaussian或者GAMESS-US,我的理解是:
1.先用AM1进行几何优化。
2.用优化的结构进行Hartree-Fock计算。
3.用第二步的分子轨道进行CIS计算(Gaussian使用Guess=Read读入第二步的checkpoint文件;GAMESS使用guess=moread,并把第二步的MO复制到$VEC中)。Gaussian的CIS自动给出激发能和谐振强度,GAMESS只给出激发能,计算谐振强度需要另外的计算步骤。
4.Gaussian中的Zindo和CIS使用方法相同,速度更快,但适用的原子有限,如果计算量不是太大的话,还可以用TD,结果更好。另外Gaussian的CIS只能算单态、三态,其它的像双态、四重态、五重态可能也会给出结果,但是分析很麻烦。用GAMESS进行CIS计算可以定义任意多重度,比较直观。