量子化学计算中帮助几何优化收敛的常用方法几何优化，也就是寻势能面极小点结构的过程。量子化学计算中几何优化不收敛是个老生常谈的问题，在各种论坛里、里都已经反复讨论过很多遍了，但是还是时常看到有人问，而且现有的讨论也都不怎么全面，所以觉得有必要撰文谈一下。所谓几何优化不收敛，也就是始终，或者很难达到收敛要求。通常会伴随着震荡行为，即受力、几何结构变化随优化步数呈现周期性趋势。解决这种问题必须在结合经验和理论知识的前提下，通过考察实际收敛的趋势，尝试各种可能奏效处理办法。本文列举一些常用的解决不收敛，也包括加速收敛的办法。其中很多方法可以相互结合使用以达到更好的效果。这里假定用户是用Gaussian，很多方法在其它程序中也可以类似地使用。先说一下收敛标准。Gaussian中判断几何优化收敛有四个标准，在默认收敛设定下，这四个标准是：最大受力<0.00045；方均根受力<0.00030；最大位移<0.00180；方均根位移<0.00120当这四个标准都满足了，达成四个YES，就宣告收敛。另外，优化过程中只要受力小于预定 的收敛限100倍，哪怕位移还没低于收敛限，则也算作已收敛。这主要考虑到势能面非常非常缓的大的柔性分子，相对于这样尺度的分子，几何结构收敛到那么精确意义不大，放宽位移收敛限避免了收敛太慢。
有时候优化出错，不是因为几何收敛问题，而是因为每一步优化中连能量计算都没能完成。优化也可能朝着明显错误的方向进行而导致难以收敛，这极有可能是理论方法、基组、电子态及其它诸多选项的设定不合理。这些方面和优化不收敛问题本身没关系，所以不会在本文提到。
1 尝试不同的优化方法
优化几何结构的方法有很多，以前我在《过渡态、反应路径的计算方法及相关问题》（hi.baidu/sobereva/item/c5a8de0b88de1c10acdc70d7）当中详细介绍过的很多搜索过渡态的方法其实和搜索势能面极小点（即几何优化）的方法本质是一致的。对于
从头算，在Gaussian03中默认的是Berny方法，这方法实际上相当于在一般的RFO方法上添油加醋，而RFO法可理解为是Newton-Raphson方法（一种很常用的非线性优化方法）的改良。G03也支持GDIIS方法（改进版，与RFO方法相结合），它对半经验方法优化是默认的，对于大体系、势能面较平坦的情况收敛得一般比RFO要快（势能面平坦表现出的特征是每一步受力小但位移大）。在G09中新支持了GEDIIS方法，并且作为了默认的优化方法，号称是目前效率最高的优化方法，自称比RFO或GDIIS在各种条件下都能更快收敛。当使用某一种方法不收敛时，可以考虑用别的方法。这三种方法可以用opt=RFO、opt=GDIIS、opt=GEDIIS来选择。
2 使用更好的Hessian矩阵
RFO、GDIIS、GEDIIS这三种优化方法在走步时都需要Hessian矩阵（力常数矩阵），然而计算Hessian矩阵是很耗时的。因此，在默认情况下，Gaussian会通过价力场的方法近似估计出初始的Hessian矩阵，在每一步优化中只精确计算梯度，利用梯度对原先Hessian
矩阵进行修正不断得到新的Hessian矩阵。所以，在默认情况下，优化从头到尾使用的Hessian矩阵都只是近似的（这也是为什么在优化出来的结构上做Freq往往会显示还没收敛，因为Freq用的是精确的Hessian矩阵）。当Hessian矩阵离精确值偏差较大，就会造成收敛缓慢，或者始终不收敛。为解决这个问题，可以用opt=calcfc，这会在优化的第一步使用精确计算的Hessian矩阵，但是仍可能后续优化过程中Hessian矩阵逐渐变得越来越不精确而依然收敛失败。opt=calcall则不仅在第一步，在后续的每一步中也都精确计算Hessian矩阵，这使得很多优化失败的情况都能得到解决，优化所需步数通常也会减少很多，而且能够保证最终优化结果准确（因为最终判断是否收敛时是基于精确的Hessian矩阵所得结果），但代价是每一步计算量会很大。如果原先设置下的优化最终由于震荡没能收敛，可以取优化过程中受力、能量最低的一个结构（可以提取坐标重新写输入文件，也可以用geom=(check,step=n)读取第n步的坐标），在calcfc乃至calcall的条件下接着进行优化。有时候，在当前理论方法下，哪怕只是calcfc都因为计算量太大（或程序限制）而难以做到， 此时也可以考虑在更低级别下计算出Hessian矩阵作为高水平优化时的初始Hessian矩阵。具体来说，比如可以在低水平方法下做freq得到精确的Hessian矩阵，然后在高水平下优化时用opt=readfc来读取它。
maxstep3 增加迭代次数上限
这个方法用处不大。一般情况下，在默认步数下如果还不能收敛，那么多数情况是发生了震荡，增大步数上限一般也无济于事。如果初始结构离最终结构比较远，而且从Gview显示的优化过程受力和位移曲线上能看到结构变化和受力有整体降低趋势，那么可以将步长上限设大，通过opt=maxcyc=N来指定为N步。但实际上Gaussian程序有内定上限（与几何结构自由度有关），即便N设得很大，但实际的步长上限并没那么大。如果到了程序内定上限还没优化完，可以用geom=check读取最后一步的结构在新的任务中接着优化。
4 使用不同的坐标系和坐标定义方式
坐标系的选择，以及坐标的定义方式，对于优化收敛速度乃至于能否收敛都有很大关系。各种坐标系各有所长，某种坐标系对于某个体系可能比较快就能收敛，但是换到另一个体系可能收敛会比较慢。当一种坐标系无法收敛时，可尝试换另外的坐标系。一般来说，对于一个体系，哪种坐标设定下能够使各个几何变量之间耦合程度越弱，就越容易收敛。这容易理解，假设每个几何变量之间都完全没有耦合（Hessian矩阵非对角元为0），那么比如3N-6维下的几何优化就等于3N-6个一维的优化，显然一维优化是很容易完成的，所以整体优化也容易完成。如果各个坐标之间耦合比较强烈，每个坐标的微小改变又会导致其它坐标上受力的明显改变，问题就比较复杂而难以拆分，收敛会变得相对困难，对优化算法的能力要求就越高。
在Gaussian中，支持以下三种坐标系下的优化。输入文件里的坐标采用哪种坐标系记录和实际优化中使用哪种坐标系完全无关。如果不注明优化过程所用的坐标系，都会自动生成冗余内坐标，并在冗余内坐标下优化。
(1) 笛卡尔坐标（Opt=Cartesian）：通常的优化不用笛卡尔坐标，因为对多数体系这种坐标下收敛都比在冗余内坐标下要慢。主要适合于含有较多分子的簇体系的优化。另外，当冗余内坐标下收敛过程因为坐标问题报错时，改用此坐标系优化可以解决。
(2) 内坐标（Opt=Z-Matrix）：对于单分子体系，在Gview自动建立的内坐标系下优化一般比在笛卡尔坐标系下优化收敛会快。这不难理解。比如对应化学键长度的坐标之间，尤其是相隔较远这样的坐标之间，没有明显的耦合，一个键稍微伸缩一下不至于使得另一个键长坐标上的受力受到明显影响。而对于笛卡尔坐标，一个原子朝笛卡尔坐标某个轴移动一下，这很可能使周围很多原子受力发生不小的改变，因此耦合较大。但注意，如果内坐标建得很烂，胡乱定义，比如坐标跨着好几个原子，则会造成耦合严重，反倒比笛卡尔坐标下收敛更慢。
内坐标不太适合有较多分子的簇体系，因为每个单体之间内坐标怎么去链接比较合适往往说不清楚，很有可能会弄得很糟。内坐标尤为不适合环状体系的优化（是指不利用虚原子
进行特殊设定的情况下），这是因为内坐标对于环状体系是以一个原子为起点绕一圈定义的，首尾虽然在几何结构上是挨着的，但是在内坐标下却没有连着。这就导致了首尾原子间本该有的键长、键角、二面角项，要被其它各个内坐标所一起等效地“凑合”着描述，这势必造成了那些内坐标之间强烈耦合。
内坐标对于高对称性体系的优化特别有用。通过巧妙地自行设定要内坐标，并结合一些虚原子作为辅助，可以大大减少要优化的坐标数目，不仅易于收敛也显著减少了计算梯度的耗时。对于一些昂贵的，而且只能通过有限差分方法计算梯度的理论方法来说，有时只有靠这种做法才能优化得动。
另外，Gaussian也支持内坐标和笛卡尔坐标的混合，可以结合二者所长。不过需要用到这种情况的场合不多，就不多提了。
(3) 冗余内坐标：N个原子的体系的结构至少由3N-6个坐标才能描述，也即内坐标。如果在内坐标基础上，额外添加一些几何变量，使总坐标数多于3N-6个，就叫冗余内坐标。这些
多出来的变量虽然对于描述体系结构是多余（冗余）的，但是它们的存在对于帮助几何优化收敛往往是很有意义的。在默认情况下，Gaussian程序会根据输入的几何坐标自动构建出冗余内坐标，构建的方法参见J.Comp.Chem,17,49。这样自动构建的冗余内坐标中只要是离得近的原子之间都有键长项（并可能伴随产生角度、二面角项）。对于单分子体系冗余内坐标下优化通常总比在内坐标下和笛卡尔坐标下优化收敛得要容易。特别是对于环状体系，在这样的冗余内坐标下优化明显比在内坐标下好得多，这是由于环中首尾原子离得很近，因此在冗余内坐标下它们之间会有对应的键长、角度、二面角项，故而环上的每个原子在优化中都能够同等地对待。
对于氢键，在自动设定的冗余内坐标中一般也会被认为是成键的。氢键长度在优化过程中有专们对应的键长项来描述，这对收敛是有益的。然而，弱一些的氢键（或者卤键之类偏弱的相互作用），以及范德华作用的原子间由于相距较远而超过了阈值，往往不会自动添加对应的键长项。为了帮助收敛或试图解决不收敛，这时应当自行添加冗余内坐标让它们之间存在键长项。在Gaussian中，可以用opt=modredundant，然后在分子坐标末尾空一行写上要添加键长项的原子序号即可，例如
3 7
12 14
8 22
这样，在优化过程中输出的冗余内坐标的信息中，就会看到新增了R(3,7)、R(12,14)、R(8,22)键长项了。
如果是很多分子的团簇的优化，例如水簇的优化，若发生震荡难以收敛，可以对所有H---O氢键都加上对应的键长项。据说给O---O也加上键长项的话还能令收敛变得更容易一些。
尽管冗余内坐标比起内坐标增加了一些变量，因此计算梯度、Hessian时原则上会略慢于用内坐标的情况，但是相对于减少了的优化步数，这点代价是不值得一提的。
5 冻结片段
一些复杂大体系有很多基团、配体、单元，复合物体系可能有很多个单体。像这样要优化的自由度非常多，相互作用复杂时往往不易收敛。在优化时，可以先冻结住一部分坐标，比如令柔性较小的配体、片段、单体中的一部分或全部内坐标不动，而对其它自由度进行优化。由于要优化的坐标数少，收敛会变得更容易，而且每一步计算梯度或Hessian矩阵的耗时也会降低很多（尤其是当使用数值梯度时）。等收敛以后，再把所有自由度都放开继续进行优化。如果对结构精度要求不是很高，那些较为刚性的自由度索性不优化也没关系。比如研究范德华作用形成的小分子二聚体时，单体通常都保持预先优化或实验测定的结构，这并不会导致明显误差。
6 调整步长上限
opt=maxstep=N（等价于IOp(1/8=N)）可以将步长上限（也称置信半径）调整到0.01N Bohr或弧度，默认是30。当按照优化方法判断出的下一步几何变量的变化超过这个数值，则实际的改变量会调节到恰好等于这个上限值（具体来说，是在置信半径对应的球面上做
最小化寻最佳的位置）。如果优化发生震荡，且每步位移不是很大而受力较大的话，往往表明优化过程陷在比较深、比较窄的势阱里，此时可以降低步长上限，比如降低到N=3乃至N=1，以避免每次都走过头而始终达不到势能面极小点。减小步长的做法不要在优化过程中过早地使用，否则由于优化过程走得太慢，所需的优化步数会增加。另外，如果用了calcall，减小步长上限不一定有好处，因为在精确Hessian下进行优化，每一步走的是比较准确的，人为篡改了步长有可能还阻碍收敛。如果初猜结构有可能偏离实际势能面极小点结构比较远，那么可以将步长上限设得比默认更大，有助于更快接近极小点结构。
注意，maxstep设的只是最初的步长上限。随着优化的进行，步长上限其实是动态调整的。如果想一直用自己的设定，需要在opt中加上NoTrustUpdate关键词。
7 调整对称性
如果优化一开始的初猜结构具有较高对称性，那么优化过程中这个对称性会一直保持。如
果实际上体系本身势能面极小点的结构并没有那么高的对称性，则很容易造成优化不能收敛。此时应当尝试人为地随意扰动一下初始结构，比如让某个二面角稍微扭转一点来破坏对称性，并且在优化中加上nosymm，则优化过程中结构不会再因为对称性被强行限制住而不能收敛。nosymm是彻底取消对称性，实际上也可以考虑用symm(PG=X)来使用指定的比当前更低阶的X点，以在消掉虚假的、造成不优化收敛的对称元素的同时依然保留本应有的对称元素。
反之，如果一个体系看起来有对称性，而且优化过程中也确实明显看到体系逐渐变得对称性越来越高，但是始终没收敛，则建议将结构取出，在Gview里打开Edit-Point Group工具，将tolerance设低一些，使Gview能判断出此分子应有的最高阶点，然后点Symmetrize，即调整结构使得体系的结构严格满足这样的点。接下来再用Gaussian优化时，就会在较高阶点下优化，这大大节省了计算能量及其导数的时间，结构还会更精确（无对称性下优化的结构虽然看上去有了对称性，但由于数值上的微小偏差其坐标并没有完全精确满足对称性），也由于要优化的自由度显著变少了，使得优化过程容易完成。有时候在Gview里已经按照指定的点对初始结构做了对称化，但是Gaussian在算的时候却
认不出相应的点，尤其是对于高阶点（比如Gaussian很难认出C60的Ih点），这是由于程序自身毛病以及输入文件中记录精度的问题导致的。可尝试在输入文件中使用内坐标而非笛卡尔坐标记录结构，也可以用symm=loose来放宽判断点的标准。
8 增加DFT泛函积分精度
DFT计算中交换相关泛函由于形式复杂，而且种类奇多，所以它涉及到的项不是利用解析积分方法精确计算而是利用数值格点积分方法来近似计算的，详见《密度泛函计算中的格点积分方法》（hi.baidu/sobereva/item/0d3734337e60bbc02f8ec2dc）。积分格点数目越多，积分就越准确，DFT计算的能量、受力、Hessian矩阵就算得越准确。如果积分精度不够高，不仅可能造成优化收敛困难，最终优化的结果也可能不很精确。Gaussian默认的泛函积分精度一般来说也够了，但如果碰上不收敛，可以试试提高积分格点数。比如使用int=fine，乃至更高级别的int=ultrafine。对于明尼苏达系列泛函，如M06-2X，对于积分格点数要求较高，因此建议至少使用int=fine。
9 改变收敛标准
改变收敛标准是没办法的办法。如果当前研究对结构精度要求不需要那么高，比如只要有个大概结构就可以，或者是作为其它什么任务（比如动力学）的初始结构，若发现结构反复微小地震荡，而且离收敛限也不远了，索性就直接停掉优化，当做已经收敛就行了。在Gaussian中也可以用opt=loose放宽收敛限。反之，如果对几何优化精度要求高，比如一些弱相互作用的研究，或是需要计算准确的振动频率，那么应该将收敛限设严，用opt=tight乃至verytight。
10 尝试不同的初始结构
比如，优化过程中某个二面角反复在两个角度来回晃荡就是不收敛，那么可以试试将初始
结构中这个二面角设在这个两个角度的平均值，以此作为初始结构（最好同时结合calcfc）。
11 尝试不同理论方法和基组
某个计算级别下如果优化很难收敛，可以尝试换个相近的基组，或换个相近的理论方法（比如B3LYP换成B3PW91之类），如果优化能收敛，可将此结构作为原先级别下优化的初始结构（也可以将最后的波函数作为初猜，Hessian矩阵也一起读入）。理论方法和基组不应和当前计算级别相差太多，否则势能面极小点位置可能相差得不少，对解决收敛问题将没多大用处。或者，你的文章中所用的理论方法干脆就换成一个精度相仿佛，对当前体系又容易收敛的方法。也可以考虑用其它适宜的方法优化得到收敛的结构，而用你原先打算用的方法在优化好的结构上讨论能量、极化率、电荷分布等性质。
如果主要目的是想节省整个优化过程的计算耗时，那么建议先在级别低一些的理论方法和
基组下（或者用比较好的半经验方法，或适宜的分子力场）做预优化。不过这对解决原先级别下不收敛的问题倒未必有多大帮助，因为低级别下与原先级别下势能面的极小点结构往往有不可忽视的偏离。另外，较低的计算级别选择得必须有意义，比如想在高精度下优化一个范德华复合物，事先用诸如B3LYP这种根本没法合理描述散作用的泛函预优化不仅白费时间，预优化出的结构还很可能比你最初搭的结构偏离极小点结构更多。
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议