解决SCF不收敛问题的方法
文/Sobereva
First release: 2010-May-17 Last update: 2014-Aug-18
量子化学常涉及SCF计算，如基于迭代的半经验方法、HF、DFT等，本文只以HF迭代为例来讨论。在HF迭代中，由Fock矩阵F对角化获得新的系数矩阵C和轨道能ε，然后构造密度矩阵D=C'C'^(T)，其中C'为不含虚轨道系数的C矩阵，再由D构造新的Fock矩阵，反复进行直到收敛，可以写为F_(1)->C_(1)->D_(1)->F_(2)->C_(2)->D_(2)...。收敛的判据不是唯一的，一般是指当前步与上一步的所有的密度矩阵元变化量的方均根以及能量的变化在很小范围内。多数收敛在20步以内能达到，但慢收敛甚至完全不能收敛的情况是容易遇到的。常见的是迭代后期能量随迭代进行在几个值上反复地循环呈现震荡，直到达到最大步数仍未收敛。也可能能量虽然震荡而总趋势是慢慢降低的，但这明显减慢了收敛速度。也可能规律性没这么强，看上去有随机性。当体系含弥散轨道、处于非平衡构型、前线区域能级密集、HOMO-LUMO间能隙较小（过渡金属化合物中容易出现）、开壳层时往往不易或不能收敛。人们为此提出了一些解决方案加速收敛和避免不收敛，本文对其中一些常见的方法进行介绍。在最后，给出解决不收敛问题的实际处理办法的简要汇总，并给出在Gaussian中的实现方法。如果懒得看理论介绍，只想立刻解决实际问题，就直接看文末的汇总就行了。
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)就削弱了当前步与上一步密度矩阵之间的差异，使密度矩阵随迭代变化更为平滑，帮助收敛。高斯在迭代初期会使用动态阻尼，可以用SCF(NDamp=N)设定直到第N次迭代都用阻尼。这个方法对收敛很有帮助，但也并非总能奏效，比如前线轨道能级密集的情况。
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更为随意，都能很快收敛，但初猜的密度矩阵离收敛比较近时DII
S收敛更快。将EDIIS 与DIIS组合使用，误差矩阵最大值较大时使用EDIIS，较小时使用DIIS，可以成为兼具二者优点的既高效又稳健的方法，这也是高斯默认的。
虽然DIIS对加速收敛、解决不能收敛的情形很有好处，但有时却导致收敛困难，此时可关掉DIIS再试，在高斯中用SCF=nodiis，若用#P查看迭代细节，会看到每次迭代输出的一大段DIIS信息不再显示。另外应注意的是DIIS收敛的结果不一定是能量全局最小点。
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)依靠ψoc c_(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)。
在高斯中使用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绝对值减小，使收敛更慢。
4 直接最小化(Direct minimization)
这是一类方法而不特指某种具体方法。总能量E是系数矩阵C（或密度矩阵D）的函数，可写作E(C)，变化轨道组合系数会引起总能量发生变化，所以通过SCF寻总能量最低的波函数实际上可视为常规的非线性优化问题。这十分类似于分子几何结构优化，都是在能量面上寻极小点。几何优化是不断调整结构到总能量最低的构型；而SCF过程则是在特定构型下，不断优化各轨道组合系数到总能量最低的体系波函数。所以几何优化常用的最陡下降法、共轭梯度法、牛顿法、准牛顿法等都可以应用到求解最优电子波函数问题上，MCSCF 也利用了牛顿法。而上述的DIIS、EDIIS则反过来也被用在几何优化上，分别称GDIIS和GEDIIS，G代表Geometry。使用直接最小化方法每步不再通过对角化Fock矩阵得到新的
系数，而是根据不同方法的定义来不断调整系数。直接最小化能保证一定收敛到一个能量极小点，对于DIIS收敛困难的体系很管用。
直接最小化方法用在求解波函数上虽然要做一些特殊考虑，但大体思路都是一样的，常见方法包括下列：
4.1 最陡下降法：每一步按照当前位置的能量负梯度方向走步，也就是能量下降最快的方向，直到走到这个方向能量最低点。最陡下降法往往在初期能量降低得较快，但由于每步之间方向是正交的，进入能量面呈山谷形状的区域后容易反复震荡，收敛缓慢。如果将步幅乘上一个调节因子来减小，避免走到此方向能量最低点，对避免震荡有好处，称调整步幅的最陡下降法(scaled steepest descent)。
4.2 共轭梯度法：令每次最陡下降法移动的方向通过前一步的方向进行校正，可缓解最陡下降法的震荡问题，加快收敛。
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)的二阶导数，每一步运算耗时很长。
4.4 准牛顿法：使用估算方法快速得到E(C)''而不是直接计算，比牛顿法大大节约时间，效率甚至接近DIIS。
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)，但速度不如DII
S快。
下面介绍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]区间内取得极

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。