解决SCF不收敛问题的方法

解决SCF不收敛问题的方法

文/Sobereva(3)
First release: 2010-May-17    Last update: 2017-Oct-25


1 前言

量子化学计算离不开SCF(自洽场)迭代,如半经验方法、HF、DFT等。在SCF迭代中,由Fock矩阵F对角化获得新的系数矩阵C和轨道能{ε},然后构造密度矩阵D=C'C'^(T),其中C'为不含虚轨道的C矩阵,再由D构造新的Fock矩阵,反复进行直到收敛,可以写为F_(1)->C_(1)->D_(1)->F_(2)->C_(2)->D_(2)...。收敛判据不是唯一的,比如Gaussian中用的判据是当前步与上一步的密度矩阵元变化量的最大值和方均根(RMS)以及能量的变化,当数值都小于一定范围就认为已经收敛了。默认判据和设定下多数情况在25步以内能达到收敛。但慢收敛甚至完全不能收敛的情况是经常会遇到的。常见的是迭代后期能量随迭代呈现震荡,直到达到默认的最大步数仍未收敛。也可能能量虽然震荡但总趋势是慢慢降低的,一直迭代下去能收敛,但震荡行为明显拖慢了收敛速度。也可能震荡的规律性不显著,迭代过程能量曲线看上去有随机性,但就是很难达到收敛限。不收敛、难收敛情况的出现有很多数值巧合因素,但有些情况不好收敛是众所周知的,例如:基组含弥散函数、体系处于明显的非平衡构型、前线区域轨道能级密集、HOMO-LUMO能隙较小(过渡金属化合物、成键方式古怪等静态相关较强体系中容易出现此情况)、限制性开壳层(RO)计算。

人们提出过一些方案加速收敛或试图解决不收敛,如DIIS、阻尼方法、温度展宽、能级移动、二次收敛方法等。主流的量化程序,如Gaussian、ORCA中默认就会使用一部分这样的帮助收敛。而这些主流的量化程序为了加速SCF计算耗时,会引入一些数值近似,比如一开始先默认用较低的DFT积分格点、忽略数值较小的积分、Incremental fock方式加速Fock矩阵的构建等等,这些近似有时候会阻碍收敛。

SCF不收敛怎么办(Gaussian中对应的错误提示是L502出错)是初学者最最最常问的问题,本文第2节将解决不收敛问题的可行办法进行了汇总,并给出在Gaussian中对应的关键词。不建议读者听信网上其它一些地方的关于解决SCF不收敛的说法,那90%都是误人子弟。如果读者对其中的DIIS、能级移动、阻尼、温度展宽和二次收敛方法的原理感兴趣,可以阅读本文第3节的详细介绍。


2 Gaussian中解决SCF不收敛问题的常用做法

绝对没有哪个或哪些关键词只要加上就一定能解不收敛问题,解决不收敛问题需要在理解关键词含义的基础上反复尝试。以下方法越靠前的越值得优先考虑,大多数可以结合在一起使用以得到更好效果。

  • (1)对于M06、M06-2X等明尼苏达系列泛函可尝试增加泛函积分格点的精度,对于其它泛函这个做法见效几率则比较有限。G09默认是int=fine,相当于(75,302),可以提升到int=ultrafine,相当于(99,590),耗时也会增加,收敛的几率也会增加,结果精度也会提升。
    注:从G16开始,默认就是int=ultrafine了。
  • (2)Gaussian会自动在计算初期时降低积分计算精度以加快计算,但有可能因此阻碍收敛。对于使用了弥散函数的情况(其它情况不算)出现不收敛时建议先尝试SCF=NoVarAcc来避免因自动降低积分精度而导致的这个问题。
  • (3)用int=acc2e=12可以增加积分精度,对因为使用大量弥散函数(其它情况不算)导致的大体系难收敛问题奏效几率较大。G09默认的是int=acc2e=10。
    注:从G16开始,默认就是int=acc2e=12了。
  • (4)Gaussian默认使用Incremental Fock方式以近似方式构建Fock矩阵来显著节约迭代过程每一步的时间,但可能因此阻碍收敛。可以用SCF=noincfock避免这个做法。
    Hint:含有很大量弥散函数的基组(如aug-cc-pVTZ及更大)计算中、大体系时,往往很难收敛,建议SCF(novaracc,noincfock) int=acc2e=12三管齐下,解决问题的几率较大,而耗时增加不算很大。但这样并非总是奏效,把弥散去掉后将收敛的波函数作为初猜也是很有效的,见下文的(9)。如果你的计算里用的弥散函数不多,或者根本都没加弥散函数,靠SCF(novaracc,noincfock) int=acc2e=12解决不收敛的概率较低,应当用本文其它方法。
  • (5)用能级移动法提升虚轨道能量以增大gap,避免虚轨道和占据轨道间过度混合:SCF=vshift=x。x一般为300~500。这并不影响最终输出结果中的轨道能,只影响收敛过程。对于gap较小的含过渡金属的体系值得优先尝试此做法。
  • (6)用SCF=conver=N关键词改变收敛限,这代表令密度矩阵RMS收敛限为1E-N,密度矩阵最大变化和能量收敛限都为1E(-N+2)。G09/16对单点计算默认是SCF=tight,相当于SCF=conver=8,然而此时对密度矩阵的收敛要求有点太严了,往往不容易达到,而且达到的时候能量通常都早已经收敛到极高精度了。对于计算单点能的目的可放心降到conver=6,相当于把收敛标准放宽100倍,通常收敛时能量变化已经非常小了。但是做几何优化、振动分析的时候不建议降低默认的SCF收敛限,否则可能结果不准确,还可能阻碍几何优化的收敛。
    不过,对于优化来说,有些体系的初猜结构可能构建得和实际结构偏差较大,导致优化初期SCF收敛难,此时可降低SCF收敛限到=6(但不要设得更松)使结构粗略收敛,再在最后的结构下用默认的SCF收敛限进一步严格优化,是很有帮助的做法。但此法对于G16不适用,因为G16中后HF和涉及能量导数的任务都强制要求收敛到至少SCF=conver=8的水平,否则无法继续进行。
  • (7)改变默认的初猜方法,可使用比如guess=core、huckel(指的是扩展Huckel方法)、INDO(实指INDO、CNDO、Huckel的混合)等关键词。默认的是使用Harris泛函的DFT方法的结果作为初猜,是通过自由原子密度叠加作为分子电子密度来构建KS算符,然后将变分得到的轨道用做后续计算的初猜,这步只做一次而不像普通DFT方法继续做迭代。core的初猜是过于粗糙的,说明初始的Fock矩阵元只含单电子项而不含双电子项,即密度矩阵用空矩阵,这默认被用作一些半经验方法计算的初猜。扩展Huckel计算也不需要迭代,也默认被用作一些半经验方法的初猜。INDO由于是迭代的方法,给高等级计算提供初猜波函数之前自身也需要初猜(扩展Huckel方法),本身亦有不收敛的可能。guess=AM1一般不能直接使用。
  • (8)使用不同的计算方法得到收敛波函数用作初猜,比如某个泛函下计算不收敛,可以尝试其它泛函或HF,如果能收敛,再用guess=read读取其收敛的波函数作为初猜。
  • (9)小基组往往比大基组易收敛,可将小基组收敛的波函数作为大基组计算时的初猜。但两个基组尺度差异太大时也往往不能收敛,可以试图逐次提升基组水平反复如此操作,如STO-3G->3-21G->6-31G*->6-311G**。虽然看上去要反复做很多次计算而很花时间,但不应忘记计算时间是随基函数数目呈指数型增长,这样做可以减少在大基组下迭代的次数,所以从总时间来看往往并不吃亏。这如同几何优化时经常先用低水平方法预优化,再转到高水平方法。
    注:实际上,用小基组(包括使用半经验方法的STO)作为大基组计算时的初猜是经过波函数的投影对系数矩阵进行转换,一般依靠小基组基函数与大基组基函数的重叠积分。对于极小基,半经验的收敛结果就是较好初猜,所以容易收敛,而对于大基组时这就是粗糙初猜,不易收敛。所以需要按照如上方法在保证每次都容易收敛的情况下,一步步提高基组等级,从而一点一点提高波函数质量。半经验方法由于只含价层轨道,所以其波函数是投影到高水平方法的价层MO做初猜,内层MO用原子的SCF收敛的轨道波函数做初猜。
  • (10)使用二次收敛方法(Quadratic convergence)。这种方法收敛所需步数比默认情况更少,有较大程度解决不收敛。SCF=QC代表使用这种方法,SCF=XQC代表先用普通方式迭代到64轮(或者自己用scf=maxcyc设的值),如果不收敛再自动切换为QC。虽然此方法解决不收敛的可能性较大,但笔者很不推荐一上来就用,而是应当先尝试其它方法。因为二次收敛方法耗时高,SCF每一步的耗时是平常的好几倍,而且如果必须借助此方法才能收敛的话,有可能最终并没有收敛到最稳定的波函数,导致得到的能量高于真实值(为求稳妥可以同时加stable关键词检测一下所得波函数稳定性)。
  • (11)使用Fermi展宽:SCF=Fermi。此时默认也会择情自动使用能级移动和阻尼。
  • (12)DIIS是默认的加速SCF收敛的方法,使得收敛所需的步数普遍减少许多。但极个别情况也可能反倒导致不收敛,可尝试关闭之:SCF=noDIIS。
  • (13)略微改变分子几何结构,比如稍微缩短/伸长键长、改变键角,若能得到收敛的波函数,可作为原来的几何构型下计算的初猜。Gaussian在几何优化中一般会自动将当前步结构的收敛波函数作为计算下一步结构的波函数的初猜,与这个策略有相似之处。
  • (14)对于开壳层体系,可以先计算相应的闭壳层离子体系得到收敛的波函数,然后读取其收敛的波函数作为初猜。注意,RO方式计算开壳层体系比非限制性计算(U)难收敛得多得多,如果RO收敛不了,先改成U再说!
  • (15)电子数少的时候往往更容易收敛,比如阳离子比阴离子一般容易收敛。所以可以先计算电离态,将得到收敛波函数用作初猜。
  • (16)用SCF=Maxcycle=N加大迭代次数上限到几百、上千。这是解决SCF不收敛问题的最蠢的方法,但由于以讹传讹,却成了初学者们最爱用的做法。如果一个搞量化的碰见SCF不收敛时第一反应是使用这个办法解决,或者看到某人的输入文件里总是带着比如SCF=maxcyc=500这样的关键词,那么他一定是个菜鸟!初学者总喜欢将N设得很大,这种做法>99.9%的情况下都绝对没用。在Gaussian默认的128次迭代次数内若不能收敛,N加大到多大几乎都是白搭,还白白浪费计算时间。仅当你观察迭代过程能量变化,发现随着迭代能量有缓慢降低,超过128步后继续迭代下去有希望收敛,才有必要加大循环次数上限,但极少会有这种情况出现。其实稍微有点基本的逻辑就能明白,倘若迭代次数上限设到几百上千就能解决不收敛,那量化程序干嘛不默认就设到几百上千?
  • (17)上述做法若都不奏效,索性改用其它方法或基组来计算当前体系。对于DFT,HF成份的越高的泛函往往越容易收敛。越小的基组也往往越容易收敛。

用IOp(5/13=1)是最脑残的“解决”不收敛的做法,文末的附录专门对此进行了大批判。

值得一提的是,化学意义越强的结构、电子态通常越容易收敛。诸如成键方式普通的有机分子在平衡结构下计算基态时SCF极少会碰到不收敛的情况。遇到非常难收敛的情况应注意是否有以下问题:

  1. 原子间距离太近
  2. 初猜的结构相对于真实的结构有严重变形
  3. 搭建的结构异想天开
  4. 对于过渡金属配合物,自旋多重度设定得不符合实际基态

3 促进SCF收敛的一些方法的原理

3.1 阻尼方法(Damping)

设正常解得的第n步的密度矩阵为D_(n),阻尼方法使实际用于构建第n+1步Fock矩阵F_(n+1)的密度矩阵变为D'_(n),D'_(n)=w*D_(n-1)+(1-w)*D_(n)。w是权重系数,既可以设为常数也可以根据迭代过程动态调整。在迭代出现震荡时,D(n),D(n+1),D(n+2)...的变化呈锯齿状,用w参数平均化后的密度矩阵D'_(n)代替D_(n)就削弱了当前步与上一步密度矩阵之间的差异,使密度矩阵随迭代变化更为平滑,帮助收敛。Gaussian在迭代初期会使用动态阻尼,这个方法对收敛很有帮助,但也并非总能奏效,比如前线轨道能级密集的情况。

3.2 迭代子空间中直接求逆(DIIS,Direct Inversion in the Iterative Subspace)

这是由Pulay发展的基于外推的方法。常规的SCF迭代收敛缓慢,DIIS能明显减少收敛所需迭代次数,加快SCF计算,是SCF计算中使用得十分普遍的方法,使多数分子都能在20步以内收敛。这个方法利用之前步的信息来估算出最好(最接近收敛,亦即“误差”最小)的Fock矩阵。下一步的Fock矩阵F_(n+1)由之前步的Fock矩阵线性组合而成,F_(n+1)=∑[i]c(i)*F_(i),c(i)是组合系数,i的加和从1到n,也可以设成比如从n-10到n。引入的每一步的F(i)都会引入误差,以误差矩阵Err_(i)表示,总误差矩阵Err_tot=∑[i]c(i)Err_(i),误差函数ErrF(c)为Err_tot的模。若找到一套系数能让误差函数ErrF(c)最小,则组合出的F_(n+1)就是最佳、最接近收敛的Fock矩阵。只需通过令ErrF对每个c(i)求导得0,并且用拉格朗日乘子法将归一化条件∑[i]c(i)=1作为限制,就能求解出各个组合系数。这等价于求解矩阵方程Ac=b,因此c=A^(-1)b,其中对矩阵A直接求逆的操作就是名字中Direct inversion的由来,iterative subspace就是指由之前迭代步的信息构成的子空间。若只利用当前步与此前的10步的信息构建F_(n+1),则子空间的维度就是10,之前陈旧的信息就被丢掉了,这也是常见做法,因为那些信息对当前迭代的贡献已经很小了,留着还会占内存。

有了上面的思路后还需要具体定义误差矩阵。第n步迭代求解HF方程就是求解F_(n)C_(n)=SC_(n)E_(n),E是对角矩阵,矩阵元为轨道能。显然用上一步的系数矩阵不可能令此式成立,否则当C_(n-1)=C_(n)时就说明已收敛了,即F_(n)C_(n-1)-SC_(n-1)E_(n)≠0,此式左边的结果正是误差矩阵,若将式子右乘C_(n-1)^(T)S,则可整理为Err_(n-1)=F_(n)D_(n-1)S-SD_(n-1)F_(n),这就是Pulay定义的误差矩阵。

也有使用以能量描述误差的相似方法,称为EDIIS,E代表Energy,它是DIIS与后文中RCA的思想结合的方法。EDIIS中下一步密度矩阵D_(n+1)=∑[i]c(i)D_(i),由于总能量是密度矩阵的函数,求出一套系数令总能量最小,就确定了D_(n+1)。注意EDIIS求解c(i)的时候限制它们必须都大于或等于0,而DIIS没这个要求,所以EDIIS被称为内插方法,而DIIS被称为外推方法,但实际上DIIS迭代中既有内插也有外推的情况。EDIIS初猜的密度矩阵比DIIS更为随意,都能很快收敛,但初猜的密度矩阵离收敛比较近时DIIS收敛更快。将EDIIS与DIIS组合使用,误差矩阵最大值较大时使用EDIIS,较小时使用DIIS,可以成为兼具二者优点的既高效又稳健的方法,这也是Gaussian默认的。

3.3 能级提升法(Level Shifting)

也称能级移动法。先介绍下SCF迭代过程的本质。设第F_(n)对角化得到的MO轨道波函数系为ψ_(n)(其中包括ψ_(n,1)、ψ_(n,2)...ψ_(n,N),N为基函数数目,能量由低到高排序),其中占据的轨道部分叫ψocc_(n),非占据的部分叫ψvir_(n),下一步n+1的Fock矩阵元F_(n+1)依靠ψocc_(n)的信息构建。以ψ_(n)为基时,显然F_(n)是对角矩阵,矩阵元就是各个轨道的能量。在收敛前以ψ_(n)为基时F_(n+1)不是对角矩阵,如果仍是对角矩阵,等于不用再做对角化,ψ_(n+1)就等于ψ_(n),就说明已经收敛了。对F_(n+1)进行对角化所得到的变换矩阵X,就是描述在这一次迭代中新得到的ψ_(n+1)是如何由上一步的ψ_(n)转换而来。由于ψ_(n)和ψ_(n+1)都是正交基,所以X为酉矩阵。X可看成由四个子方阵组成,分别是对角块X[occ(n)_occ(n+1)]和X[vir(n)_vir(n+1)],以及非对角块X[occ(n)_vir(n+1)]和X[vir(n)_occ(n+1)](二者互为转置关系)。由于占据轨道波函数之间的酉变换不改变体系能量,而虚轨道对能量有没有贡献,所以X[occ(n)_occ(n+1)]和X[vir(n)_vir(n+1)]的数值对迭代能量变化并没有贡献,只有对应占据轨道与虚轨道混合的非对角块X[vir(n)_occ(n+1)]才会使总能量变化。

当前步解出来的虚轨道与占据轨道能量差值越小,说明它们越相似,HF方程的近似解也越容易是由它们为基组合而成,所以下一步所得的占据轨道就会混进很多当前步的虚轨道成分,使总能量变化剧烈,当然总能量下降得快是好事,但也因此可能出现总能量不降反升的情况。若每步虚轨道与占据轨道能量差都不大,则每步它们都有大规模的混合,可能能量降低与升高会混杂出现而产生振荡的情况,必然使收敛缓慢。Saunders等人提出的能级提升法就是人为地将每步得到的虚轨道能量提升(对应的虚轨道波函数也因此改变),这样每步占据轨道与虚轨道的混合就会减小,如果提升得足够多,可以证明一定能在迭代中令总能量不断降低,而不会出现升高的情况,如此可以使能量最终收敛于极小点。虽然虚轨道能量提升越多,收敛越稳定,总能量变化越缓慢,越能保证总能量每步迭代都能降低,但如果原始的迭代过程本来就能让每步的总能量都降低,提升虚轨道能量则会使降低的量也减小,导致收敛速度减慢,还有可能出现低于占据轨道能量的空轨道,形成所谓的“空穴”,这样的非基态结果没什么意义,所以能量也不能肆意提升。能量提升值的选取有一定任意性,对不同虚轨道可以统一提升同一个值也可以提升不同值。

将第n步虚轨道的能量值提升x也就是使得E_(n)对应虚轨道能量的矩阵元都增加了x成为E'_(n),因此需要修改F_(n)成为F'_(n),这样通过求解F'_(n)C'_(n)=SC'_(n)E'_(n)就能得到所期望的E'_(n)以及相应的轨道波函数C'_(n)。用于限制性闭壳层SCF,可以推出F'_(n)=F_(n)+x*(S-0.5*SD_(n)S)。

在Gaussian中使用SCF=Vshift=n关键词可以设定将虚轨道能量提升n*0.001hartree,一般n设为几百,若仍不收敛可尝试提高更多,这个方法对解决不收敛问题很有效,可以尝试多提升一些。如果体系本身就容易收敛则不要用这个方法,会使收敛更慢。注意能级移动只是帮助收敛的方法,在最终输出的能量中会从所得虚轨道能量中减去这个值,修正回实际的虚轨道能量,所以对结果没影响。

PS:这个问题可以通过二阶微扰理论做更细致的分析。将f_(n+1)与f_(n)的差作为微扰项f~,对第i个占据轨道的一级校正波函数δψ(i)=∑[j∈vir]H(i,j)/(e(i)-e(j))*ψ_(n,j),其中H(i,j)=<ψ_(n,i)|f~|ψ_(n,j)>,j为虚轨道,e是相应轨道的能量。所以ψ_(n+1,i)≈ψ_(n,i)+δψ(i),因此轨道i能量变化为δe(i)=∑[j∈vir]H(i,j)^2/(e(i)-e(j)),体系总能量的变化δE=∑[i∈occ]δe(i)。可见轨道能级差e(i)-e(j)越小,从虚轨道j混进占据轨道i的量越多,δE的值也会越负。

当F~较大时,由于高阶项不可忽略,加上后可能实际的δE是正值,对收敛不利,如果将H(i,j)乘上一个数值较小的阻尼因子,令高阶项可忽略,就能让δE仍是负值,这是上述阻尼方法另外形式的表述。

能级移动法则相应于将e(i)-e(j)项替换为e(i)-e(j)-x。在一般情况下由于占据轨道能量低于虚轨道,因此e(i)-e(j)总是负值,此时能量总是降低的,这没问题。但迭代过程中aufbau定律(即电子按轨道能量从低往上排)并不总是遵守的,这种情况经常出现在过渡金属络合物的体系。原因是波函数离收敛比较远时,分子的平均势场与真实势场偏差较大(“真实”是指相对当前步更真实的下一次迭代的平均势场),此时SCF给出的轨道能量因此与真实偏差较大,甚至顺序是倒转的,这样程序根据此时SCF轨道能量由低到高判断出的轨道占据顺序就与真实顺序不符。当顺序相反时,e(i)-e(j)就成为正值,可能δE也会成为正值,导致迭代后总能量不降反升,这光靠乘上阻尼因子已无法解决,此时若用较大的能级提升量x使e(i)-e(j)-x项为负,就能使能量仍然降低。若e(i)-e(j)本来为负值,从中减去了x就会使其绝对值变大,造成δE绝对值减小,使收敛更慢。

3.4 直接最小化(Direct minimization)

这是一类方法而不特指某种具体方法。总能量E是系数矩阵C(或密度矩阵D)的函数,可写作E(C),变化轨道组合系数会引起总能量发生变化,所以通过SCF寻找总能量最低的波函数实际上可视为常规的非线性优化问题。这十分类似于分子几何结构优化,都是在能量面上寻找极小点。几何优化是不断调整结构找到总能量最低的构型;而SCF过程则是在特定构型下,不断优化各轨道组合系数找到总能量最低的体系波函数。所以几何优化常用的最陡下降法、共轭梯度法、牛顿法、准牛顿法等都可以应用到求解最优电子波函数问题上,MCSCF也利用了牛顿法。而上述的DIIS、EDIIS则反过来也被用在几何优化上,分别称GDIIS和GEDIIS,G代表Geometry。使用直接最小化方法每步不再通过对角化Fock矩阵得到新的系数,而是根据不同方法的定义来不断调整系数。直接最小化能保证一定收敛到一个能量极小点,对于DIIS收敛困难的体系很管用。

直接最小化方法用在求解波函数上虽然要做一些特殊考虑,但大体思路都是一样的,常见方法包括下列:
3.4.1 最陡下降法:每一步按照当前位置的能量负梯度方向走步,也就是能量下降最快的方向,直到走到这个方向能量最低点。最陡下降法往往在初期能量降低得较快,但由于每步之间方向是正交的,进入能量面呈山谷形状的区域后容易反复震荡,收敛缓慢。如果将步幅乘上一个调节因子来减小,避免走到此方向能量最低点,对避免震荡有好处,称调整步幅的最陡下降法(scaled steepest descent)。
3.4.2 共轭梯度法:令每次最陡下降法移动的方向通过前一步的方向进行校正,可缓解最陡下降法的震荡问题,加快收敛。
3.4.3 牛顿法:令E(C)对C进行Taylor展开,舍去三阶及以上项作为近似,并根据dE/dC=0的条件,推得C_(n+1)=C_(n)-E(C)'|C_(n) * (E(C)''|C_(n))^(-1),这里|C_(n)代表将C_(n)代入左边的变量。显然对于纯二次型区域(是指可以用多变量的二次函数准确描述的区域,三阶及以上导数为0)这个表达式是准确的,可以一步达到极小点。在接近收敛时能量面可近似视为二次区域,所以牛顿法收敛很快,但是在优化初期效率很低,一方面是二次区域的近似不再那么合理,相对于DIIS迭代次数可能需要更多;另一方面是牛顿法需要计算E(C)的二阶导数,每一步运算耗时很长。
3.4.4 准牛顿法:使用估算方法快速得到E(C)''而不是直接计算,比牛顿法大大节约时间,效率甚至接近DIIS。
3.4.5 增强型牛顿法(Augmented Newton):单纯的牛顿法并不稳定,需要特殊的处理,其中通过修改原始E(C)''来调整步进长度和方向的方法称为增强型牛顿法。

往往不将C或D作为E的变量,而是将前后两次MO的变换矩阵X作为变量,好处是MO是正交归一集,相互独立利于优化。前面已说到,在对能量有影响的只是X的非对角块,所以更具体来说E是X[vir(n)_occ(n+1)]的函数。对E(X)的求导可写为以MO为基的单电子和双电子积分,由于需要从AO的双电子积分变换过来,所以很耗时。

将D作为变量时,根据密度矩阵的性质可将能量写为E(D)=Tr(DF)=Tr(Dh)+Tr(DG(D)),其中h是和G是Fock矩阵中单、双电子部分。D并不能自由变化,由于MO需要正交归一,D需要受制于等幂条件DSD=D,否则结果无意义。然而在优化中满足这个要求并不容易,一种方法是通过3*D^2-2*D^3运算将不等幂的D纯化。等幂条件限制对应于轨道占据数为0或者1,而RCA(relaxed constraints algorithms)过程允许在优化过程中放松这个要求,占据数可以在0至1之间,相当于将等幂条件放松为DSD<=D,而到收敛时自动恢复等幂条件,RCA当中最简单的实现为ODA(Optimal Damping Algorithm),但速度不如DIIS快。

下面介绍ODA的大致过程。在ODA迭代过程中的任意一段可以这么表示:...D_(n)~-->F_(n+1)~-->D_(n+1)-->D_(n+1)~-->F_(n+2)~-->D_(n+2)...。其中D_(n+1)~=D_(n)~+λ(D_(n+1)-D_(n)~),λ是通过令能量E(D_(n+1)~)在λ∈[0,1]区间内取得极小值来得到的。由D_(n)~构建F_(n+1)~,再由其本征向量矩阵C_(n+1)构建D_(n+1),如上方法计算D_(n+1)~,就形成了一次迭代。每一步可以视作由当前位置D_(n)~向D_(n+1)-D_(n)~方向走到能量最低点,这个方向正是能量下降最快的方向。

在Gaussian中可以用SCF选项中的SD和SSD来使用最陡下降法和调整步幅的最陡下降法,但容易出错,不建议使用。SCF=QC是线搜索法和加入了阻尼的准牛顿法的结合,称二次收敛法(quadratically convergent),虽然往往能解决收敛问题,但速度非常慢,另外可能在收敛末期出现振荡而无法收敛。建议尝试SCF=XQC,在迭代初期仍用效率高的DIIS方法,难以收敛时才转为速度慢但收敛性好的SCF=QC,使收敛性和速度都较好,但往往在DIIS的阶段就已经报错停止。

3.5 温度展宽

也叫分数占据法(FON, Fractional occupation number),在Gaussian中通过SCF=Fermi使用,默认是关闭的。这个方法的特点是使轨道占据数n出现非整数,其数值根据费米-狄拉克分布计算,即n(i)=1/(1+exp((e(i)-e(F))/kT)),e(F)是费米能级,T是虚构的温度,密度矩阵D(i,j)=2*∑[k]n(k)C(i,k)C(j,k),由于不再分清轨道是否占据,k要对所有轨道加和。此方法用的费米-狄拉克分布函数并没有实际的物理意义,只是由于此方法意图将一部分占据轨道的电子挪到虚轨道上而借用了这个函数的形式。这样的结果就是每一步占据轨道和虚轨道之间混合加剧,前面也已提到,这种混合是使能量降低的驱动力,所以此方法可以加快收敛速度。实际应用中,此方法对于容易收敛的体系收敛速度没什么影响,对一些收敛慢的体系可以大大加快收敛速度,有时也能解决个别体系不能收敛的问题。温度展宽可以结合DIIS、阻尼一起使用以得到更好效果。

e(F)和T都是可变参数,影响FON的效果。e(F)一般设为HOMO与LUMO能量的平均值,实际上设在HOMO之下LUMO之上都可以。T越高,在费米能级附近的分数占据的轨道越多,轨道混合程度越大;T越低,则越趋于还原为整数占据。对于HOMO-LUMO能隙大的分子,T应设大来以起到明显的加速效果;而能隙小的分子由于轨道混合已比较重,不需要太高的T。由于单slater行列式要求轨道占据数必须是整数,T在迭代完成前必须回归为0以满足整数占据条件,因此T是动态变化的。在迭代开始,视分子的能隙或DIIS误差矩阵或将温度设在1000~2000K范围内,之后逐步降低T。可以线性降低,比如每步100K,也可以将T作为DIIS误差矩阵的函数,这样随迭代进行DIIS误差矩阵逐渐减小最终趋近0,T也随之逐渐降至0,对迭代的影响也逐渐减弱,这样在迭代完成时占据数就自动恢复为整数了。


另外,还有一种通过Car-Parrinello动力学来获得基态收敛波函数的方法,当其它各种解决波函数不收敛的方法都不灵的时候,这个方法就是最后选择。但这种方法不能在Gaussian中实现。也就是一开始将电子温度设得较高,在动力学过程中慢慢降温,降温到0时就得到了收敛的基态波函数,也就是电子退火过程。




附:对Gaussian中IOp(5/13=1)的大批判

通过IOp(5/13=1)来“解决”SCF不收敛是计算化学界最最最最脑残的行为,没有之一!我不吝啬任何恶毒的词语来批判用5/13=1的做法。洒家一看到一些人输入文件里的5/13=1就特别来气,简直就像看到大街上有一坨恶翔!

5/13=1意思是,SCF到达指定圈数不收敛也不报错,继续做后面的事。为什么用这个IOp脑残,稍微有一点点量子化学理论常识的人都自然明白。这哪是解决不收敛啊,分明是对不收敛视而不见!完全就是饮鸩止渴!如果到最大圈数时,恰好已经离收敛限不远了(并且收敛限本身设得并不太松),那么结果可以用,但这必须自己查阅输出文件来确认这一点;然而,初学者才不知道查阅输出文件来检查收敛到什么程度了呢!如果到达最大圈数时能量还在明显震荡,由于用了5/13=1又没报错,初学者就会傻乎乎地直接用最后给出的能量,然而这能量根本就不靠谱啊!!!用这样的结果害人害己,轻则结果不准确,重则结果定性错误!可以说,对于不知道自行检查SCF收敛程度的初学者用5/13=1,完全就是对数据不负责任,甚至个别情况可以说是数据造假!

很多对理论知识一无所知的初学者发现平时不收敛会报错的任务一加5/13=1就不报错了,于是就傻乎乎地、无条件地所有任务都加上5/13=1,甚至还让别人加5/13=1,别人甚至连5/13=1都不知道是干什么的,被以讹传讹也傻乎乎地用,这种肤浅至极的做学问的态度真是不忍直视!!!不懂某个关键词是干什么的话就不要加这个关键词,不要盲目地效仿别人(尤其是网上大量计算化学的讨论都是错的),不要被以讹传讹!

5/13=1是计算化学界的大毒瘤!也不知道始作俑者是谁,如今简直就像瘟疫一样迅速传播!笔者在笔者建立的计算化学公社论坛、思想家公社QQ群都全力抵制使用5/13=1的行为。另外,只要看到有人用5/13=1,你就可以认为这人是比使用scf(maxcyc=几百)还菜的菜鸟,一定要予以强烈的鄙视!要让5/13=1无处现形,让大家都以使用5/13=1为耻!另外强烈呼吁Gaussian官方,在发现使用5/13=1的时候,一定要在输出文件里输出巨大的警告提示,甚至干脆取消掉这个IOp,以避免初学者们受害。

另外,在http://bbs.keinsci.com/forum.php?mod=viewthread&tid=3344帖子中Gaussian官方客服也已经明确表态绝不应当在任何实际计算中用5/13=1。

顺带一提,有人根据自己的一些计算经验,认为5/13=1对于解决几何优化过程的SCF不收敛是有益的。这种观点的原理是:几何优化初期由于结构不太好,SCF不容易收敛,于是加了5/13=1让程序不报错;等到优化末期,结构收敛到差不多极小点结构附近了,SCF通常也就容易收敛了,此时虽然带着5/13=1,但是照样也能收敛到了标准收敛限,大不了人工检查一下是否SCF收敛了。实际上,即便是出于这种考虑,用5/13=1也是不对的!这里给个例子:用5_13=1有害的例子.rar
从例子可见,用了5/13=1时,由于波函数未收敛从而进一步导致受力计算烂,而又不试图通过恰当解决SCF不收敛的做法去解决,导致结构也反复震荡难以收敛。而按照本文做法,加了辅助SCF收敛的关键词后,每一步SCF都收敛了,结构也顺利收敛了。所以那些人鼓吹适合在几何优化过程中5/13=1的说法也是不恰当的。总之一句话:无论是菜鸟还是有经验的用户,无论是出于何种目的,都不要用5/13=1!

已有 2 条评论

  1. zhlu929

    这个图真的好有意思

  2. wtli

    看到配图忍不住笑出声

评论已关闭