摘抄:Gaussian L9999问题(不收敛问题的解决方法)

L9999报错运行情况

一般来讲,L9999出现时大多运行到当前运行状态的最后的part还未收敛造成的。

输出文件一般都有点大(行数会比较超规格)

比如上图,因为循环了太多次,算到20万行了(输出文件也有14M,正常情况在10M以下)。

特征行:Error termination via Lnk1e in l9999.exe

最后计算时间会偏长:

Elapsed time: 0 days 5 hours 0 minutes 34.2 seconds.

情况1 虚频不合 (*1)

1
2
3
Optimization stopped.
-- Wrong number of Negative eigenvalues: Desired= 1 Actual= 4
-- Flag reset to prevent archiving.

成因:Gaussian会在过渡态优化时自动检查虚频的数量,若不为1即终止。这一检查没有必要

解决:将NoEigenTest关键词加入到opt的选项中,如 opt=(TS, CalcFC, NoEigenTest)

情况2 步数内未收敛 (*1)

1
2
3
Optimization stopped.
-- Number of steps exceeded, NStep= 100
-- Flag reset to prevent archiving.

成因:几何优化未能在指定步数内收敛

解决:打开GaussView,File-Open,给下方的Read Intermediate Geometries(Gaussian Optimation Only)选项打勾:

然后点击Result-Optimization,可以看到优化时的能量变化曲线。

观察优化曲线应纵向放大,方法为用鼠标反复选择纵轴方向仅包含最末几个点、但横轴方向包含所有点的“宽而矮”的矩形范围,直到能看清最末几个点的变化情况为止。

根据图形判断是否震荡。如下图,需纵向放大至右边才能看清确实震荡了。(注意图中偶尔有突然的、无“周期性的”某1~2个点的spike是正常现象,不是震荡,下方小幅度但有“周期性”的才是震荡。)

(左侧看不出是否震荡。需纵向放大至右侧才能看出)

Gaussian中判断几何优化收敛有四个标准,在默认收敛设定下,这四个标准是:
最大受力<0.00045;方均根受力<0.00030;最大位移<0.00180;方均根位移<0.00120
当这四个标准都满足了就宣告收敛。方均根受力/方均根位移体现的是体系中所有原子的平均受力/位移情况。另外,优化过程中只要受力小于预定的收敛限100倍,哪怕位移还没低于收敛限,也算作已收敛。

如果没有震荡

【仅当初始结构离最终结构比较远,特别是体系比较小(此时默认的步数上限较少)或者用了较小的步长上限时,且从gview显示的优化过程受力和位移曲线上能看到结构变化和受力有整体降低趋势,那么重算或续算时才有必要将步数上限设大,通过opt=maxcyc=N来指定上限为N步。】

opt=maxcycle=当前步数的2~3倍。可以读取当前最末/最优的一个构象做为新的初猜。注意检查构象、轨迹是否合理,如果结构已经不合理了,则应适当调整初猜后再做(尤其是找过渡态的时候)。

如果有振荡

参考自此页

如果在原先设置下的优化任务最终没能收敛,可以取优化过程中受力或能量最低的一个结构(可以用gview切换到相应的帧,保存成新的输入文件。也可以用geom(check,step=n)从%chk指定的chk文件中读取第n步的坐标),然后加上下文提到的关键词接着进行优化。

尝试不同的优化方法

当使用某一种方法不收敛时,可以考虑用别的方法。这三种方法可以用opt=RFO、opt=GDIIS、opt=GEDIIS来选择。个人强烈推荐对于弱相互作用体系、柔性大体系优化遇到不收敛时优先尝试opt=GDIIS。

使用更好的Hessian矩阵

opt=calcfc可在优化的第一步使用精确计算的Hessian矩阵,但是仍可能后续优化过程中Hessian矩阵逐渐变得越来越不精确而依然收敛失败。

opt=calcall不仅在第一步,在后续的每一步中也都精确计算Hessian矩阵,这使得很多优化失败的情况都能得到解决,优化所需步数通常也会减少很多,而且能够保证最终优化结果准确(因为最终判断是否收敛时是基于精确的Hessian矩阵所得结果),但代价是每一步计算量会很大。

使用不同的坐标系和坐标定义方式

在Gaussian中,支持以下三种坐标系下的优化。输入文件里的坐标采用哪种坐标系记录和实际优化中使用哪种坐标系完全无关。如果不注明优化过程所用的坐标系,都会自动生成冗余内坐标,并在冗余内坐标下优化。

(1) 笛卡尔坐标(Opt=Cartesian)

主要适合于含有较多分子的簇体系的优化。另外,当冗余内坐标下收敛过程因为坐标问题报错时(比如构成二面角的四个原子中的三个在优化过程中排成了直线,此时无法定义二面角),改用此坐标系优化可以解决。注意,使用GDIIS时即便写了cartesian,Gaussian依然会在默认的冗余内坐标下优化。

(2) 内坐标(Opt=Z-Matrix)

内坐标对于高对称性体系的优化特别有用。

(3) 冗余内坐标(默认)

N个原子的体系的结构至少由3N-6个坐标才能描述,也即内坐标。如果在内坐标基础上,额外添加一些几何变量,使总坐标数多于3N-6个,就叫冗余内坐标。

调整对称性

体系对称性不高:人为地随意扰动一下初始结构,比如让某个二面角稍微扭转一点来破坏对称性。

体系对称性较高:将结构取出,在Gview里打开Edit-Point Group工具,将tolerance设低一些,使Gview能判断出此分子应有的最高阶点群,然后点Symmetrize,即调整结构使得体系的结构严格满足这样的点群

增加DFT泛函积分精度

Gaussian默认的泛函积分精度一般来说也够了,但如果碰上不收敛(特别是对于明尼苏达系列泛函,如M06-2X),可以试试提高积分格点数,比如使用int=ultrafine。注:从G16开始,默认就是int=ultrafine了。

改变收敛标准

用opt=loose放宽收敛限;或用opt=tight/verytight提高几何优化精度要求。

尝试不同的初始结构

比如,优化过程中某个二面角反复在两个角度来回晃荡就是不收敛,那么可以试试将初始结构中这个二面角设在这个两个角度的平均值,以此作为初始结构(最好同时结合calcfc)。此做法解决问题的几率不大,但是偶尔碰巧能解决,毕竟本身几何优化这种迭代的数值过程能否收敛就是有很多巧合性。如果优化过程的结构变化已经向着不靠谱或者非预期的方向上走了,往往是因为初始结构明显不合理所致,需要重新调整结构重新优化再试。

尝试不同理论方法和基组

某个计算级别下如果优化很难收敛,可以尝试换个相近的基组,或换个相近的理论方法(比如B3LYP换成B3PW91之类),如果优化能收敛,可将此结构作为原先级别下优化的初始结构(也可以将最后的波函数作为初猜,Hessian矩阵也一起读入)。理论方法和基组不应和当前计算级别相差太多,否则势能面极小点位置可能相差得不少,对解决收敛问题将没多大用处。

总结

具体怎么解决优化不收敛得看实际情况,特别是要用gview检查优化过程的结构变化趋势和能量、受力曲线,看是否最后开始震荡了(如果能量/受力还有降低的大趋势则应当继续跑下去)。

虽然上面讨论了很多,但从实用性出发,综合上述内容和原作者的建议我总结的方法使用顺序:

  1. 检查初始结构是否存在问题【opt非常需要】;
  2. 扰动结构调整对称性【opt非常需要】;
  3. 如果体系小,尝试opt=calcall。体系大尝试opt=recalc=3或5;
  4. 尝试opt=gdiis;
  5. 尝试opt(gdiis,maxstep=x,notrust),x为3~5;
  6. 如果上述做法都不灵,干脆改成loose收敛限,但如果之后做振动分析容易有虚频。

另外论坛帖子作者的建议:

我个人倾向于找到能量较低、4个收敛标准离收敛限较近,且结构合理的某点为新的初猜,并首先尝试在 opt 选项中加上(MaxStep=5, NoTrustUpdate, GDIIS) 选项;如果有计算频率的资源,可再加上opt=calcfc选项。特别仅对于对称性较高的结构,如T、C3 等群的复杂分子,若反复调整关键词后仍然震荡,可以尝试构建或破坏具有对称性的初猜结构,分别结合opt=Cartesian 优化,如 opt=(MaxStep=5,NoTrustUpdate,Cartesian), 有时有效果。

过渡态时的处理方法

寻找过渡态也是几何优化过程,只不过是向势能面一阶鞍点优化,而不是向极小点优化。

在Gaussian里一般用这样的关键词:opt(calcfc,noeigen,TS)

寻找过渡态也经常会遇到几何优化不收敛问题。由于在主流量化程序里优化过渡态和优化极小点的算法基本没有区别,因此上面所有解决优化极小点不收敛的方法,对于优化过渡态也完全适用。

但优化过渡态对初始结构的敏感性远远高于优化极小点,因此优化过渡态过程几何优化不收敛有很大可能是因为初猜的过渡态结构离实际过渡态结构不够接近,导致最终被优化到没有化学意义的结构所造成的。

到底是由于初猜结构不好,还是由于最终发生了震荡导致优化过渡态没有收敛,通过gview查看优化过程的结构变化立刻就能知道。

另外,由于优化过渡态对于Hessian矩阵的质量比起优化极小点敏感得多得多,因此解决优化过渡态不收敛的问题,在计算条件允许的情况下,建议优先考虑calcall或recalc=3(对于小体系,由于计算Hessian矩阵耗时不算特别高,对于不太好找过渡态的情况,为了减少反复尝试的次数,建议总是用calcall或recalc)。

本文摘抄内容来自: