量子化学计算中帮助几何优化收敛的常用方法
笔者注：读者不要说诸如“此文我看了，还是不收敛”、“我用了文中的方法还是不收敛”这类话。绝对没有哪个方法能100%解决几何优化不收敛。本文绝对没有一条馊主意，所有Gaussian支持的真正有可能能解决几何优化不收敛的做法在本文都已经全面列举了，没有任何遗漏。倘若本文的做法都做了尝试还没解决，那也别指望有任何其它方法能解决。才随便试了文中一个做法，或者以不正当的方式（没有真正了解原理的基础上）瞎试几个就说没解决问题，这根本什么也说明不了。要想解决问题，必须十分仔细阅读本文，所有有可能对当前问题有用的方法都依次尝试，有必要的时候几个一起结合使用。本文列举的做法已经把最坏的情况都考虑了，所以出现所谓的“看了本文还是没解决问题”的情况一定是读者还没真正仔细看、仔细领会、充分尝试。
量子化学计算中帮助几何优化收敛的常用方法
Common methods to facilitate geometry optimization convergence in quantum chemistry calculations
文/Sobereva @北京科音
First release: 2012-Oct-13 Last update: 2022-May-5
0 前言
几何优化就是寻找势能面极小点结构的过程。而所谓几何优化不收敛，就是优化过程始终，或者很难达到收敛要求，这种情况通常会伴随着震荡行为，即原子受力、几何结构随优化步数增加呈现一定周期性变化趋势。显然盲目增加优化步数的上限试图解决是超级愚蠢的做法。解决这种问题必须在结合经验和理论知识的前提下，通过考察实际收敛的趋势，尝试各种可能奏效处理办法，本文就专门说说这个问题。更详细、更具体的东西在笔者讲授的北京科音（http://www.keinsci.com）的量子化学培训班里会系统地讲解。
经常有人问Gaussian程序做几何优化为什么有时会出现下面这种报错
Error termination request processed by link 9999.
Error termination via Lnk1e in d:\study\G09W\l9999.exe at Sat Jun 02 20:45:20 2018.
实际上往上看，就会看到更明确的错误提示
Optimization stopped.
-- Number of steps exceeded, NStep= [当前步数上限]
-- Flag reset to prevent archiving.
即达到了步数上限还没收敛，于是优化就报错结束了。这多数情况就是因为出现了上述震荡问题而造成了难收敛所导致的。另外，还经常有初学者问，怎么优化任务算了非常非常久还没算完，这很可能也是因为发生了震荡，导致优化始终达不到收敛限所致（比如都跑了好几百步了，还在震荡中）。
如果想确认是否发生了震荡，可以用GaussView载入输出文件（打开文件的界面下方必须勾选上Read intermediate geometries），然后进入Results - Optimization，看能量和受力变化曲线。如果曲线到后期没有整体降低、逐渐收敛到平坦的大趋势，而是来回上下波动，并且观看优化轨迹也发现结构在反复波动，这就说明震荡了。几何优化正在执行的过程中也可以随时对收敛情况这样进行监视。
说怎么解决几何优化难收敛之前先说一下收敛标准。Gaussian中判断几何优化收敛有四个标准，在默认收敛设定下，这四个标准是：
最大受力<0.00045；方均根受力<0.00030；最大位移<0.00180；方均根位移<0.00120
当这四个标准都满足了就宣告收敛。方均根受力/方均根位移体现的是体系中所有原子的平均受力/位移情况。另外，优化过程中只要受力小于预定的收敛限100倍，哪怕位移还没低于收敛限，也算作已收敛。这主要考虑到势能面非常非常缓的大的柔性分子，相对于这样尺度的分子，几何结构收敛到那么精确意义不大，放宽位移收敛限避免了收敛太慢。
也注意有时候优化任务中途出错，不是因为几何收敛方面的问题，而是因为优化当中某一步连能量计算都没能正常完成，很常见的情况是由于SCF不收敛导致的，这种情况的解决见《解决SCF不收敛问题的方法》（http://sobereva.com/61）。几何优化过程也可能朝着明显错误的方向进行而导致难以收敛，这可能是理论方法、基组、电子态、计算模型的构建，以及其它某些设定不合理所致。这些方面和优化不收敛问题本身没关系，所以不会在本文提到。
如果在原先设置下的优化任务最终没能收敛，可以取优化过程中受力或能量最低的一个结构（可以用gview切换到相应的帧，保存成新的输入文件。也可以用geom(check,step=n)从%chk指定的chk文件中读取第n步的坐标），然后加上下文提到的关键词接着进行优化。
下面，就详细谈谈碰到几何优化不收敛（通常由于上述震荡问题导致）应当怎么解决。本文列举各种有可能解决不收敛，也包括加速收敛的办法。其中很多方法可以相互结合使用以达到更好的效果。本文虽然说的是Gaussian的情况，但很多方法在其它程序如ORCA中也可以类似地使用。对于优化过渡态过程出现难收敛的情况，在本文末尾有专门的提及。
1 尝试不同的优化方法
优化几何结构的方法有很多，以前我在《过渡态、反应路径的计算方法及相关问题》（http://sobereva.com/44）当中详细介绍过的很多搜索过渡态的方法其实和搜索势能面极小点（即几何优化）的方法本质是一致的。对于从头算，在Gaussian03中默认的是Berny方法，这方法实际上相当于在一般的RFO方法上添油加醋，而RFO法可理解为是Newton-Raphson方法（一种很常用的非线性优化方法）的改良。G03也支持GDIIS方法（改进版，与RFO方法相结合），它对半经验方法优化是默认的，对于大体系、势能面较平坦的情况收敛得一般比RFO要快（势能面平坦表现出的特征是每一步受力小但位移大）。在G09中新支持了GEDIIS方法，并且作为了默认的优化方法，号称是目前效率最高的优化方法，自称比RFO或GDIIS在各种条件下都能更快收敛，但实际情况中往往并非如此。当使用某一种方法不收敛时，可以考虑用别的方法。这三种方法可以用opt=RFO、opt=GDIIS、opt=GEDIIS来选择。个人强烈推荐对于弱相互作用体系、柔性大体系优化遇到不收敛时优先尝试opt=GDIIS。
2 使用更好的Hessian矩阵
RFO、GDIIS、GEDIIS这三种优化方法在走步时都需要Hessian矩阵（质权坐标下也叫力常数矩阵），然而计算Hessian矩阵是很耗时的。因此，在默认情况下，Gaussian会通过价力场的方法近似估计出初始的Hessian矩阵，在每一步优化中只精确计算梯度，利用梯度对原先Hessian矩阵进行修正不断得到新的Hessian矩阵。所以，在默认情况下，优化从头到尾使用的Hessian矩阵都只是近似的（这也是为什么在优化出来的结构上做Freq往往会显示还没收敛，因为Freq用的是精确的Hessian矩阵）。当Hessian矩阵离精确值偏差较大，就会造成收敛缓慢，或者始终不收敛。为解决这个问题，可以用opt=calcfc，这会在优化的第一步使用精确计算的Hessian矩阵，但是仍可能后续优化过程中Hessian矩阵逐渐变得越来越不精确而依然收敛失败。opt=calcall则不仅在第一步，在后续的每一步中也都精确计算Hessian矩阵，这使得很多优化失败的情况都能得到解决，优化所需步数通常也会减少很多，而且能够保证最终优化结果准确（因为最终判断是否收敛时是基于精确的Hessian矩阵所得结果），但代价是每一步计算量会很大。
在G16中，增加了opt=recalc=n关键词来设置每n步重新精确计算一次Hessian矩阵。显然n越小，收敛的可能性越大，但耗时越高。n=1时对应calcall。实际中n建议设为3~5，以达到耗时和帮助收敛效果的权衡。
有时候，在当前计算级别下，哪怕只是calcfc都因为计算量太大（或程序限制）而难以做到，此时也可以考虑在更低级别下计算出Hessian矩阵作为高水平优化时的初始Hessian矩阵。具体来说，比如可以在低水平方法下做freq得到精确的Hessian矩阵，然后在高水平下优化时用opt=readfc来读取它。
3 增加步数上限
这个方法是菜鸟最爱用的方法，但实际上几乎不会起到用处！如果一个搞量化的碰见几何优化不收敛时第一反应是使用这个办法解决，那么他一定是个菜鸟！之所以这个做法这么流行，一方面是长期的以讹传讹，另一方面是初学者们想当然地以为迭代次数越多自然越可能收敛。然而，绝大多数情况下，在默认步数下如果还不能收敛，那么多数情况是早已发生了震荡，增大步数上限也无济于事，只会白费功夫，糟蹋计算资源。
仅当初始结构离最终结构比较远，特别是体系比较小（此时默认的步数上限较少）或者用了较小的步长上限时，且从gview显示的优化过程受力和位移曲线上能看到结构变化和受力有整体降低趋势，那么重算或续算时才有必要将步数上限设大，通过opt=maxcyc=N来指定上限为N步。注意实际上Gaussian程序在16版之前有内定步数上限（与坐标变量数正相关，也就是计算输出中显示的"maximum allowed number of steps"），即便N设得巨大，若超过了内定步数上限，N将自动设为内定步数上限值。
4 使用不同的坐标系和坐标定义方式
坐标系的选择，以及坐标的定义方式，对于优化收敛速度乃至于能否收敛都有很大关系。各种坐标系各有所长，某种坐标系对于某个体系可能比较快就能收敛，但是换到另一个体系可能收敛会比较慢。当一种坐标系无法收敛时，可尝试换另外的坐标系。一般来说，对于一个体系，哪种坐标设定下能够使各个几何变量之间耦合程度越弱，就越容易收敛。这容易理解，假设每个几何变量之间都完全没有耦合（Hessian矩阵非对角元为0），那么比如3N-6维下的几何优化就等于3N-6个一维的优化，显然一维优化是很容易完成的，所以整体优化也容易完成。如果各个坐标之间耦合比较强烈，每个坐标的微小改变又会导致其它坐标上受力的明显改变，问题就比较复杂而难以拆分，收敛会变得相对困难，对优化算法的能力要求就越高。
在Gaussian中，支持以下三种坐标系下的优化。输入文件里的坐标采用哪种坐标系记录和实际优化中使用哪种坐标系完全无关。如果不注明优化过程所用的坐标系，都会自动生成冗余内坐标，并在冗余内坐标下优化。
(1) 笛卡尔坐标（Opt=Cartesian）：通常的优化不用笛卡尔坐标，因为对多数体系这种坐标下收敛都比在冗余内坐标下要慢。主要适合于含有较多分子的簇体系的优化。另外，当冗余内坐标下收敛过程因为坐标问题报错时（比如构成二面角的四个原子中的三个在优化过程中排成了直线，此时无法定义二面角），改用此坐标系优化可以解决。注意，使用GDIIS时即便写了cartesian，Gaussian依然会在默认的冗余内坐标下优化。
(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 调整步长上限
opt=maxstep=N（等价于IOp(1/8=N)）可以将步长上限（也称置信半径）调整到0.01N Bohr或弧度，默认是30。当按照优化方法判断出的下一步几何变量的变化超过这个数值，则实际的改变量会调节到恰好等于这个上限值（具体来说，是在置信半径对应的球面上做最小化寻找最佳的位置）。如果优化发生震荡，且每步位移不是很大而受力较大的话，往往表明优化过程陷在比较深、比较窄的势阱里，此时可以降低步长上限，比如降低到N=5乃至N=3，以避免每次都走过头而始终达不到势能面极小点。减小步长的做法不要在优化过程中过早地使用，否则由于优化过程走得太慢，所需的优化步数会增加。另外，如果用了calcall，减小步长上限不一定有好处，因为在精确Hessian下进行优化，每一步走的是比较准确的，人为篡改了步长有可能还阻碍收敛。
注意，maxstep设的只是最初的步长上限。随着优化的进行，步长上限其实是动态调整的。如果想一直用自己的设定，需要在opt中加上NoTrust关键词。对于结构反复微小震荡的情况，在使用maxstep减小步长上限的同时总建议同时写上此关键词。
6 调整对称性
如果优化一开始的初猜结构具有较高对称性，那么优化过程中这个对称性多数情况下会一直保持。如果实际上体系本身势能面极小点的结构并没有那么高的对称性，则很容易造成优化不能收敛，而且即便收敛了也会有虚频。此时应当尝试人为地随意扰动一下初始结构，比如让某个二面角稍微扭转一点来破坏对称性。
反之，如果一个体系看起来有对称性，而且优化过程中也确实明显看到体系逐渐变得对称性越来越高，但是始终没收敛，则建议将结构取出，在Gview里打开Edit-Point Group工具，将tolerance设低一些，使Gview能判断出此分子应有的最高阶点群，然后点Symmetrize，即调整结构使得体系的结构严格满足这样的点群。接下来再用Gaussian优化时，就会在较高阶点群下优化，这大大节省了计算能量及其导数的时间，结构还会更精确（无对称性下优化的结构虽然看上去有了对称性，但由于数值上的微小偏差其坐标并没有完全精确满足对称性），也由于要优化的自由度显著变少了，使得优化过程容易完成。有时候在Gview里已经按照指定的点群对初始结构做了对称化，但是Gaussian在算的时候却认不出相应的点群，尤其是对于高阶点群（比如Gaussian很难认出C60的Ih点群），这是由于程序自身毛病以及输入文件中记录精度的问题导致的。可尝试在输入文件中使用内坐标而非笛卡尔坐标记录结构，也可以用symm=loose来放宽判断点群的标准。
7 增加DFT泛函积分精度
DFT计算中交换相关泛函由于形式复杂，而且种类奇多，所以它涉及到的项不是利用解析积分方法精确计算而是利用数值格点积分方法来近似计算的，详见《密度泛函计算中的格点积分方法》（http://sobereva.com/69）。积分格点数目越多，积分就越准确，DFT计算的能量、受力、Hessian矩阵就算得越准确。如果积分精度不够高，不仅可能造成优化收敛困难，最终优化的结果也可能不很精确。Gaussian默认的泛函积分精度一般来说也够了，但如果碰上不收敛（特别是对于明尼苏达系列泛函，如M06-2X），可以试试提高积分格点数，比如使用int=ultrafine。注：从G16开始，默认就是int=ultrafine了。
8 改变收敛标准
改变收敛标准是没办法的办法。如果当前研究对结构精度要求不需要那么高，比如只要有个大概结构就可以，或者是作为其它什么任务（比如动力学）的初始结构，若发现结构反复微小地震荡，而且离收敛限也不远了，索性就直接停掉优化，当做已经收敛就行了。在Gaussian中也可以用opt=loose放宽收敛限。反之，如果对几何优化精度要求高，比如一些存在很弱相互作用的体系的研究（如H2二聚体），或是需要计算准确的振动频率，那么应该将收敛限设严，用opt=tight乃至verytight。
9 尝试不同的初始结构
比如，优化过程中某个二面角反复在两个角度来回晃荡就是不收敛，那么可以试试将初始结构中这个二面角设在这个两个角度的平均值，以此作为初始结构（最好同时结合calcfc）。此做法解决问题的几率不大，但是偶尔碰巧能解决，毕竟本身几何优化这种迭代的数值过程能否收敛就是有很多巧合性。如果优化过程的结构变化已经向着不靠谱或者非预期的方向上走了，往往是因为初始结构明显不合理所致，需要重新调整结构重新优化再试。
10 尝试不同理论方法和基组
某个计算级别下如果优化很难收敛，可以尝试换个相近的基组，或换个相近的理论方法（比如B3LYP换成B3PW91之类），如果优化能收敛，可将此结构作为原先级别下优化的初始结构（也可以将最后的波函数作为初猜，Hessian矩阵也一起读入）。理论方法和基组不应和当前计算级别相差太多，否则势能面极小点位置可能相差得不少，对解决收敛问题将没多大用处。这种做法属于撞大运，起效几率不大，但实在没辙了可以试试。或者，你的文章中所用的理论方法干脆就换成一个精度相仿佛，对当前体系又发现相对更容易收敛的方法。
如果主要目的是想节省整个优化过程的计算耗时，那么建议先在级别低一些的理论方法和基组下（或者用比较好的半经验方法，或适宜的分子力场）在loose收敛限下做预优化。不过这对解决原先级别下不收敛的问题倒未必有多大帮助，因为低级别下与原先级别下势能面的极小点结构往往有不可忽视的偏离。另外，较低的计算级别选择得必须有意义，比如想在高精度下优化一个范德华复合物，事先用诸如B3LYP这种根本没法合理描述色散作用的泛函预优化不仅白费时间，预优化出的结构还很可能比你最初搭的结构偏离极小点结构更多。
11 关于溶剂模型
溶剂环境下的几何优化有没有必要带溶剂模型在《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》（http://sobereva.com/327）里我已经专门说过了，不再累述。这里特别强调的一点是，对于opt freq任务，不建议用常用的SMD溶剂模型，因为我和很多人都发现用SMD的时候几何优化可能变得更困难，而且收敛后做振动分析时还容易出现虚频，这大抵是SMD溶剂模型使得势能面出现一定数值噪音所致。用Gaussian的scrf关键词默认的IEFPCM溶剂模型则没有显著这种问题。但使用IEFPCM也不收敛的话，而且若当前体系的优化没有绝对必要带溶剂模型，也可以尝试在真空下优化再试。
几何优化和优化后计算单点两个过程对溶剂效应的考虑、用的溶剂模型都没有必要非得统一。在优化之后算单点能时可以放心用SMD。
具体怎么解决优化不收敛得看实际情况，特别是要用gview检查优化过程的结构变化趋势和能量、受力曲线，看是否最后开始震荡了（如果能量/受力还有降低的大趋势则应当继续跑下去）。虽然上面讨论了很多，但从实用性出发，对于大多数情况，对于初学者笔者推荐按照以下顺序去尝试解决震荡问题：
1 对于DFT，特别是明尼苏达系列泛函（如M062X、M06、M06L），首先加上int=ultrafine再试（从G16开始此为默认，不再需要尝试此项）。如果用了SMD，改用IEFPCM
2 如果体系小，尝试opt=calcall。如果体系大算得慢，此关键词会令计算量进一步激增，实在没辙再考虑这个，或者尝试opt=recalc=3或5
3 尝试opt=gdiis
4 尝试opt(gdiis,maxstep=x,notrust)，x为3~5
5 如果上述做法都不灵，干脆改成loose收敛限，但如果之后做振动分析容易有虚频
寻找过渡态，或者说优化过渡态，本身也是几何优化过程，只不过是向势能面一阶鞍点优化，而不是向极小点优化。在Gaussian里一般用opt=TS方式基于一个初猜结构来寻找过渡态，通常用这样的关键词：opt(calcfc,noeigen,TS)。寻找过渡态也经常会遇到几何优化不收敛问题。由于在主流量化程序里优化过渡态和优化极小点的算法基本没有区别，因此上面所有解决优化极小点不收敛的方法，对于优化过渡态也完全适用。但优化过渡态对初始结构的敏感性远远高于优化极小点，因此优化过渡态过程几何优化不收敛有很大可能是因为初猜的过渡态结构离实际过渡态结构不够接近，导致最终被优化到没有化学意义的结构所造成的。到底是由于初猜结构不好，还是由于最终发生了震荡导致优化过渡态没有收敛，通过gview查看优化过程的结构变化立刻就能知道。对于不是很简单的化学反应，研究者一次就能给出能最终优化到过渡态的初猜结构并不容易，故失败时不要轻易气馁，要耐心反复尝试。Gaussian还支持QST2方法优化过渡态，它是根据用户提供的反应物和产物结构通过插值产生过渡态初猜结构。但这个初猜结构往往很不合理，明显不如有化学直觉的用户自己手动摆出来的，因此一般情况不建议用QST2。而当你用QST2时发现结构向着不靠谱的方向优化，那应当立马放弃QST2而改用opt=TS。另外，优化过渡态对于Hessian矩阵的质量比起优化极小点敏感得多得多，因此解决优化过渡态不收敛的问题，在计算条件允许的情况下，建议优先考虑calcall或recalc=3（对于小体系，由于计算Hessian矩阵耗时不算特别高，对于不太好找过渡态的情况，为了减少反复尝试的次数，建议总是用calcall或recalc）。