计算化学公社
标题: Gaussian中几何优化收敛后Freq时出现NO或虚频的原因和解决方法 [打印本页]
作者Author: sobereva 时间: 2015-1-17 00:58
标题: Gaussian中几何优化收敛后Freq时出现NO或虚频的原因和解决方法
Gaussian中几何优化收敛后Freq时出现NO或虚频的原因和解决方法
Reasons and solutions for solving appearance of NO or imaginary frequency when convergence of geometric optimization is reached in Gaussian
文/Sobereva @北京科音 First release: 2015-Jan-16 Last update: 2024-Feb-20
1 前言
众所周知，做振动分析之前必须先在相同级别下优化。总有Gaussian使用者问，优化收敛了，也显示四个YES了，怎么做freq之后四个收敛判断标准却有的出现NO，甚至还出现虚频了？此文件就专门说一下。Gaussian中opt=TS方式的搜索过渡态的算法和几何优化很相似，因此下面的讨论对于过渡态的优化也完全适用，只不过要消的虚频显然不是感兴趣的那个（看振动动画可判断），而是额外存在的多余的虚频。
特别要强调的是，按本文的做法，绝对不存在消不掉的虚频，笔者从来没见过例外（除非犯了低级错误，诸如振动分析之前没做几何优化，或几何优化和振动分析不在严格相同级别下进行等等）。不要在网上问我怎么按照此文的做法还是没消掉虚频，消不掉一定是没认真领会此文并充分尝试。不存在本文没提到的其它能解决虚频的做法。
当遇到数值非常小的虚频时，对于计算能量相关问题，并非必须不惜一切代价得把所有虚频消掉不可。关于这点我有博文做了深入讨论，非常建议阅读：《谈谈有小虚频时热力学量的计算》（http://sobereva.com/699）。
2 关于几何优化收敛的判断和出现虚频的原因
牛顿法是高效、最常用的寻找多元函数驻点的方法，在量化研究中，这个方法被广泛用来搜索过渡态和优化几何结构。然而，牛顿法的每一步需要计算一次Hessian矩阵（更具体来说是其逆矩阵），也就是能量对几何坐标的二阶导数矩阵，如果进一步考虑质量权重后就是所谓的力常数矩阵。对于较高计算级别或大体系，计算一次Hessian矩阵是相当耗时的，另外很多方法在Gaussian里没有二阶解析导数（如CCSD、MP4等），对于这样的方法计算大体系的Hessian矩阵更是耗时甚巨甚至无法承受。因此，为了节省优化时间，一般量子化学程序默认用的都是准牛顿法，也就是每一步不精确计算Hessian，而是只计算受力，通过受力和上一步的Hessian矩阵来近似得到这一步的Hessian矩阵。虽然准牛顿法由于用的是近似的Hessian矩阵，比起基于精确Hessian的牛顿法需要更多步数才能找到过渡态/极小点，但由于每一步的耗时低得多（比如准牛顿法的10步耗时可能只相当于牛顿法的1、2步），所以整体还是比牛顿法便宜得多得多，所以是大多数量化程序默认的优化方法。当Gaussian发现当前结构的受力以及位移（下一步结构相当于当前结构的变化）都收敛了，优化就停止。受力和位移都按照最大值和方均根（RMS）来检测，因此共有四个判断标准，即四个都YES时就宣告收敛。
做频率计算，则必须先精确计算一次Hessian矩阵（或者读取精确的Hessian矩阵），才能经过一番处理得到频率，参见《基于fch中的Hessian矩阵计算振动频率的简单程序Hess2freq》（http://sobereva.com/328）里的介绍。freq任务之后，Gaussian会自动根据当前结构的受力和现成的精确Hessian矩阵也用上述四个判断标准考察一下当前结构是否真的准确收敛到了极小点。
正因为几何优化默认用的是近似的Hessian矩阵，而频率计算用的是精确Hessian矩阵，因此判断收敛时，在后两项，即最大位移和均方根位移上的判断会存在差异（受力那两项是一致的）。常见的情况是优化的最后一步四个标准都显示YES了，但是到了freq任务时，基于精确Hessian矩阵判断则发现后两个并非都YES，甚至出现了虚频，这说明优化的精度不是很高。特别是虚频如果较大的话，更是体现当前结构偏离极小点可能挺明显。
下面是个典型的例子，opt最后的输出为
Maximum Force 0.000401 0.000450 YES
RMS Force 0.000197 0.000300 YES
Maximum Displacement 0.001791 0.001800 YES
RMS Displacement 0.000951 0.001200 YES
Freq最后的输出为
Maximum Force 0.000401 0.000450 YES
RMS Force 0.000197 0.000300 YES
Maximum Displacement 0.002231 0.001800 NO
RMS Displacement 0.001325 0.001200 NO
那么，opt最后都YES了但freq最后没有都YES，结果还能用么？笔者认为，如果数值比默认情况的收敛限高不太多，不超过2倍，且没有虚频出现，一般还是可以接受的。比如上面的情况就基本可以接受。但是，如果你对频率/结构的准确度要求高，或者都出现了虚频，则应当做更精确的优化。
3 提升优化和数值精度解决虚频的做法
下面介绍让优化更精确来避免freq最后出现NO乃至虚频的做法。
3.1 几何优化时用更严的收敛限
在opt里面写上tight可以将默认的优化收敛限设严一个数量级以上。由于极小点定位得更精准了，出现虚频的可能性自然大大降低了。但代价是达到收敛限的难度大大增加，因此优化达到收敛所需的步数会大大增加、耗时大大提升，还增加了因为震荡而一直达不到较严格的收敛限的概率。
注意如果写了opt=tight freq，则freq最后也会用tight标准来判断是否收敛，而实际上freq最后只要能满足默认收敛限判断标准就够了，不需要非得满足tight那么严苛的判断标准。
值得一提的是优化过程中即便位移没收敛，但只要受力小于收敛限的100倍也自动被Gaussian算作收敛。我认为这主要是考虑到势能面非常平缓的柔性大分子或分子团簇，相对于其体系的尺度，把结构收敛到那么精确意义也不大。这种情况下，连opt最后都有NO，显然freq最后也肯定有NO。若对结构精度要求高的话可以把优化的收敛限设严。
3.2 提升DFT积分格点精度
对于DFT计算要考虑交换-相关泛函的积分格点质量问题。如果不了解什么是DFT的积分格点，可参看《密度泛函计算中的格点积分方法》（http://sobereva.com/69）。如果积分格点质量不够理想，DFT计算的受力和Hessian的精度就不够高，这往往是造成虚频出现的原因。此时即便用了tight收敛限也解决不了问题，而且改用tight后可能还极难收敛、容易震荡。如果你是Gaussian 09的用户，用opt=tight时应同时搭配int=ultrafine使用比默认的int=fine更高精度的积分格点。从Gaussian 16开始int=ultrafine已成为默认，因此就不需要手动写了。
对积分格点的要求与使用的泛函、体系都有密切联系。一些明尼苏达系列泛函（比如M06-2X等，但不是所有）对于积分格点精度要求比B3LYP、PBE0、TPSS等其它常见泛函高得多。因此比如M06-2X在int=ultrafine下依然有时会因为积分格点精度不够而导致很微小的虚频的出现，此时可尝试更高档次的积分格点int=superfine，但这非常昂贵。
对于DFT计算极个别色散作用主导的弱相互作用体系，有时候更需要用int=superfine。比如笔者曾经在ωB97XD/def2-TZVP(-f)下计算C30嵌套结构，老是有二点几波数的微小虚频死活消不掉。最终有一次我虽然没设更严收敛限，但是用了int=superfine，可算把虚频给搞掉了。但这种情况极度罕见，切勿视为普遍情况。只是说对于某些弱相互作用体系，哪怕不是明尼苏达系列泛函，当你碰到极微小虚频且实在干不掉时，还是值得得试一把int=superfine。
3.3 使用精确的Hessian矩阵
在opt里写上calcall，那么几何优化过程中每一步都会精确计算Hessian，自然最终能优化得比基于近似Hessian的准牛顿法更准确。此时几何优化和freq所用的Hessian都是一致的，因此对收敛的判断结果也是完全相同的。也就是说，只要opt最后是YES则freq最后也必然是YES，并且多数情况可以确保无虚频（值得一提的是，用calcall的时候，优化任务结束后会自动做振动分析，因此完全没有必要再单独写freq关键词要求做振动分析了）。
不过，每一步优化都精确计算Hessian太耗时。如果你之前已经在默认情况下优化过，可以取最后的结构重新优化，并在opt里写上calcfc，这样只在优化第一步的时候精确计算Hessian，而之后还是近似估计Hessian，这样收敛后也往往可以避免freq之后出现NO和虚频。但是如果初始结构离实际极小点太远则这么做无效，因为等到结构收敛时Hessian矩阵可能又偏离精确Hessian比较远了。
从Gaussian 16开始，opt里加入了一个很有用的选项recalc。比如opt=recalc=n就代表第一步精确计算一次Hessian矩阵，之后每n步重新计算一次Hessian矩阵。这种做法的效果和耗时都介于calcfc和calcall之间。对计算Hessian成本不很高的情况，n取3~5左右比较合适，对计算Hessian成本较高的情况，n取5~10比较合适，能达到效果和耗时的权衡。读者请随机应变。
4 其它导致虚频的因素和解决方法
4.1 结构对称性问题
有时候虚频是由于优化所用的初始结构对称性太高导致的。比如联苯在基态下实际上两个苯环间是有一定二面角的，但优化联苯所用的初始结构如果是纯平面的，则优化过程就会一直保持平面状态，最终也得到平面构型，显然会发现有对应于两个平面彼此间扭转的虚频。碰到这种情况，最佳的解决方法是在GaussView观看振动模式的界面中选中虚频模式，选上Manual Displacement并拉动滑条来沿着虚频模式略微调整结构，然后点save structure得到新的结构，再拿这个新结构保存新的输入文件，重新优化和做振动分析后通常就已经没虚频了。有时候对于大体系，也可能因为某些局部区域的初始对称性太高，一直维持到最后而出现虚频，这时候也应按照虚频调节结构。对于其它情况，这么按照虚频调节结构的做法解决虚频的几率比较有限，虽然也可以尝试一下。而对于有多个虚频的情况，想通过这种做法一次性就把所有虚频都解决没太大可能。我经常看到网上有很多人但凡只要看到虚频就总是妄图通过这种按照虚频调结构的办法试图消掉，这是明显不对的！！！真正应该优先尝试的是本文前面说的做法！
4.2 SMD隐式溶剂模型问题
SMD隐式模型流行度很高，在《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》（http://sobereva.com/327）我专门介绍过。我强烈不建议对opt和freq任务用SMD模型，因为此模型下很容易由于数值噪音而产生虚假的虚频，而且还可能由此造成几何优化难收敛。因此当优化和振动分析必须在溶剂模型下做时，强烈建议用Gaussian的scrf关键词默认的IEFPCM溶剂模型。从几何优化结果来讲，SMD相较IEFPCM也没有丝毫额外的好处，如327号博文所述，SMD比IEFPCM长处在于额外多考虑了非极性部分，对计算溶剂环境下的能量很重要，而非极性部分对优化的结构影响甚微。所以，opt freq在IEFPCM下做，而之后算能量时再在SMD下做，是比较得当的做法。
4.3 IOp定义泛函的问题
在Gaussian中可以自定义泛函，见《Gaussian中非内置的理论方法和泛函的用法》（http://sobereva.com/344）和《在Gaussian中自定义范围分离泛函的方法》（http://sobereva.com/550），这需要用到IOp关键词。很多人用自定义泛函的时候还是按习惯写上opt freq关键词，这是不行的。要知道，IOp是没法传递到下一个任务的，结合opt freq用的时候只对opt有效，而对freq无效，自然会导致振动分析和优化时用的泛函定义不同、对应的势能面不同，故很可能因此出虚频。所以用自定义泛函的时候opt和freq必须分开做。
4.4 程序bug问题
Gaussian从G09 D.01版开始加入DFT-D3，相关介绍见《谈谈“计算时是否需要加DFT-D3色散校正？”》（http://sobereva.com/413）、《DFT-D色散校正的使用》（http://sobereva.com/210）、《乱谈DFT-D）《http://sobereva.com/83》。DFT-D3有零阻尼和BJ阻尼两个版本，分别用em=gd3和em=gd3bj关键词使用。后者最常用，因为算相互作用能的精度整体更好一点。BJ阻尼版本在G09 D.01里的解析二阶导数存在bug，会导致你结合DFT-D3(BJ)做优化完全精确的情况下，结合DFT-D3(BJ)做的振动分析也可能出现虚频。因此解决办法是用G09 E.01及以后的版本，或者用零阻尼版本DFT-D3(0)做opt和freq。
PS1：准牛顿法、牛顿法在不同量子化学程序里具体实现不同，有很多数值方面的技巧和额外考虑，本文中没有细谈。想了解更多的话可以看看《过渡态、反应路径的计算方法及相关问题》（http://sobereva.com/44）以及《几何优化、过渡态搜索、IRC综述与原文合集》（http://bbs.keinsci.com/thread-105-1-1.html）当中的详细讨论。Gaussian做opt默认用的方法并不是简简单单的准牛顿法，而是牵扯到很多细节，在手册的opt部分里有专门的描述。
PS2：有人问，为什么有的时候优化时明明用了calcall，但是基于其结构再做freq任务的时候显示的判断标准却依然和优化最后一步不同，甚至又出现了NO。可能原因有三 (1)优化时用的级别、数值设定（如DFT积分格点）等因素和几何优化时不完全一致 (2)用GaussView读取优化的输出文件，然后又保存成freq任务的输入文件的时候，由于小数位数有限，因此造成了一点数值误差 (3)几何优化的时候用了GDIIS（默认的GEDIIS也有一定GDIIS的成份），由于这种方法预测下一步位移的时候还会参考之前步的信息，因此会和freq时基于牛顿法判断出的位移有所不同。
作者Author: Shannon 时间: 2015-1-17 05:21
有没办法让opt=calcall每一步都输出Frequency呢？ 如果让gaussian每步都算hessian却不计算出每一步的Freq，有点太浪费资源了。
如果能算出每一步的freq的话，在gaussianview 里面看每一步虚频的方向，这样万一跑歪了也能比较方便的人工矫正。
作者Author: sobereva 时间: 2015-1-17 06:36
有些人喜欢手动去消虚频什么的，其实只要Hessian是精确的，优化算法自动会向着能够消除虚频的方向走。而所谓的手动去消虚频，则难以把握好尺度，移动小了没意义，移动大了则当前虚频没消掉反倒还产生了新的虚频（况且变量间还是耦合的，手动按照某个振动模式移动考虑不了这点）。所以calcall的话就不建议这么人为移动结构了。
至于每步都输出频率，利用IOp实现非标准计算路径可能能行。但没法指定gview读取哪一步的，而且如果不对输出文件做一些手工修改来帮助gview识别，可能一个都不认。
作者Author: Shannon 时间: 2015-1-17 07:12
用M062X优化过渡态结构，学校没好电脑，舍不得用int=ultrafine，常被m62x的 虚假势能面极小点坑害，所以想用calcall来找到到底是哪一步的时候过渡态优化跑歪了。
作者Author: jiangning198511 时间: 2015-1-17 12:04
Sob大侠每次出手都是精品之作 
作者Author: sobereva 时间: 2015-1-17 14:56
int=ultrafine所付出的代价远远远远小于calcall。
calcall用了之后计算量和原来就不是一个数量级了，但int=ultrafine不会比原先多耗时太多，原先算得动那么用了这个也肯定能算得动。
作者Author: nkallwar 时间: 2015-1-26 12:37
如果是比较不同结构的能量，那么其中一个结构用了 int=ultrafine ，其余的结构是不是都得再用相同精度的积分格子重新算一边才能进行比较？
作者Author: sobereva 时间: 2015-1-26 12:49
是的
作者Author: nkallwar 时间: 2015-1-26 12:51
那是只需要计算一遍能量，还是要重新优化一些呢，前面的结构里并没有出现 opt收敛而freq不收敛的情况
作者Author: sobereva 时间: 2015-1-26 13:09
默认的int=fine和int=ultrafine对结构影响甚微，所以只需要在相同格点设定下算能量就够了
作者Author: nkallwar 时间: 2015-1-28 13:45
遇到一个奇怪的问题，opt和freq的结果让我ft了， 首先是opt的log文件：
******************************************
%chk=M06-Cr-Sym-simpleMe-S7-int.chk
%nprocshared=12
Will use up to 12 processors via shared memory.
%mem=800MW
---------------------------------------------------------------------
# opt=(calcfc,tight) genecp guess=mix geom=connectivity int=ultrafine
使用了opt=tight 和int=ultrafine ，基组是 6-31G(d) for C N O H, Lanl08(f) for Cr,
然后收敛看着很好：
Item Value Threshold Converged?
Maximum Force 0.000000 0.000015 YES
RMS Force 0.000000 0.000010 YES
Maximum Displacement 0.000040 0.000060 YES
RMS Displacement 0.000009 0.000040 YES
但得到的结构看上去上感觉很别扭，优化得到的坐标和freq作业输入文件如下：
%chk=M06-Cr-Sym-simpleMe-S7-int-freq.chk
%nprocshared=12
%mem=800MW
# freq genecp geom=connectivity int=ultrafine um06
dfstructure
0 7
Cr -1.35587600 -0.13969200 0.22924400
Cr 1.35588000 0.13972300 0.22925000
N -2.70474300 1.33577800 0.47363400
N -2.82567500 -1.37891500 -0.54018000
N 2.70469100 -1.33578800 0.47364300
N 2.82572500 1.37892200 -0.54014300
O -0.00000400 -0.00009000 1.47906000
O 0.00000700 0.00011500 -0.99049600
C -2.20390300 2.59099800 1.02088400
C -3.97177500 1.26454000 0.16519300
C -4.65512700 0.18810100 -0.38935800
H -5.70338300 0.30686600 -0.58204600
C -4.06571400 -1.03367800 -0.72061300
C -2.39533600 -2.68415300 -1.01949400
C 2.20382200 -2.59099100 1.02090800
C 3.97172200 -1.26459300 0.16518300
C 4.65510900 -0.18818000 -0.38937200
H 5.70335600 -0.30698700 -0.58208000
C 4.06574400 1.03363200 -0.72060100
C 2.39543700 2.68418700 -1.01942800
H -4.56949100 2.14311400 0.36652500
H -4.72192000 -1.76414600 -1.17580300
H -1.98432800 -3.27211800 -0.20383400
H -3.21748400 -3.24523900 -1.45342100
H -1.62214600 -2.56620400 -1.76991400
H -1.65099300 3.14209900 0.26548000
H -3.01377300 3.22057200 1.37660500
H -1.53482600 2.38698900 1.84585300
H 3.21759600 3.24522900 -1.45339200
H 1.62220400 2.56628500 -1.76981100
H 1.98450300 3.27217200 -0.20374400
H 4.72197600 1.76407500 -1.17579200
H 4.56940900 -2.14319000 0.36650500
H 3.01367900 -3.22057800 1.37663700
H 1.53474600 -2.38695400 1.84586900
H 1.65090300 -3.14209100 0.26550900
频率分析结果一大堆虚频，收敛情况一塌糊涂：
Item Value Threshold Converged?
Maximum Force 0.023783 0.000450 NO
RMS Force 0.007623 0.000300 NO
Maximum Displacement 0.068771 0.001800 NO
RMS Displacement 0.028796 0.001200 NO
这是什么情况？
作者Author: sobereva 时间: 2015-1-28 14:14
你的几何优化和频率计算明显不是同一个级别进行的，出现一堆虚频是必然的
作者Author: nkallwar 时间: 2015-1-28 14:47
丢人了，做opt的时候居然忘写UM06了， 居然还屁颠屁颠发出来被大神看笑话，不写的话默认就是HF了吧
作者Author: sobereva 时间: 2015-1-28 23:33
是的，默认是HF/STO-3G
作者Author: 虎王 时间: 2015-1-29 22:14
如果优化加了int=ultrafine，是否计算能量、轨道、跃迁等都需要加int=ultrafine？
作者Author: sobereva 时间: 2015-1-30 00:24
没这个必要
作者Author: nkallwar 时间: 2015-6-3 20:43
本帖最后由 nkallwar 于 2015-6-3 20:44 编辑
最近总出现这样的问题
两个displacement判据早早就变成了0
然后两个force判据就不动了，最后link 9999 错误
Item Value Threshold Converged?
Maximum Force 0.022910 0.000450 NO
RMS Force 0.002915 0.000300 NO
Maximum Displacement 0.000000 0.001800 YES
RMS Displacement 0.000000 0.001200 YES
****
Item Value Threshold Converged?
Maximum Force 0.022910 0.000450 NO
RMS Force 0.002915 0.000300 NO
Maximum Displacement 0.000000 0.001800 YES
RMS Displacement 0.000000 0.001200 YES
*****
后面都是这样
这种情况是什么原因造成的？应该用什么方法帮助收敛呢？
作者Author: sobereva 时间: 2015-6-3 20:53
gdiis和calcall都试试。
注意考察优化过程结构的变化
作者Author: dxyan 时间: 2016-2-25 09:06
本帖最后由 dxyan 于 2016-2-25 18:16 编辑
我的体系就是这样，优化算完后4个YES，然后算频率第三项就变成了NO，而且奇怪的是我每个体系算下来都是这样的。昨天按照楼主说的方式用了opt=tight int=ultrafine由在之前的最后结构上重新做了优化，优化后依然4个YES，但是算出的能量比之前稍微高了一些。然后用这个优化的结果重新做了光计算频率的输入文件，此时没选tight，结果频率计算下来反而后两项都变成了NO，而且还出现了一个虚频。不明白为什么啊？
作者Author: sobereva 时间: 2016-2-25 18:28
你优化时候用opt=tight int=ultrafine，在freq的时候就必须也得用int=ultrafine
作者Author: dxyan 时间: 2016-2-25 20:38
注意：如果写opt=tight freq，则freq最后也会用tight标准来判断是否收敛。实际上只要能满足默认判断标准即可，不需要非得满足tight的判断标准
------------------------------------------那前面说算频率不小用tight啊，这算频率为什么要加int=ultrafine呢？
作者Author: sobereva 时间: 2016-2-25 22:15
你几何优化用什么积分格点，频率计算就必须用什么积分格点，否则不是同一个级别
“算频率不小用tight啊”是病句，理解不能
作者Author: dxyan 时间: 2016-2-26 08:43
打错字了，不好意思。
。我意思是之前帖子说算频率不需要非得满足tight的判断标准，我理解就是算频率不需要加上tight。因为之前帖子说opt=tight，必须加上int=ultrafine。所以后来算频率这2个都没加。
作者Author: dxyan 时间: 2016-3-5 21:35
请问如何在gaussianview 里面看每一步虚频的方向，并且人工矫正啊？
作者Author: dxyan 时间: 2016-3-5 21:52
sob老师，请教一下，比如我想优化A,B两个分子，A由于结构优化4个YES，再算频率时出现了2个NO，然后用opt=tight int=ultrafine重新优化了结构，频率计算的时候只加了int=ultrafine来解决这个问题。这样算下来，能量的绝对值比不加opt=tight int=ultrafine时大。那我要比较A,B结构的能量，如果B没有出现出现opt收敛而freq不收敛的情况，是不是B也要用opt=tight int=ultrafine重新优化，并且加int=ultrafine再做频率计算。
另外我要算A和C，B和C的相互作用能，那么C也要用opt=tight int=ultrafine优化吗？AC和BC的复合物在优化的时候是否也要用opt=tight int=ultrafine？谢谢！
作者Author: sobereva 时间: 2016-3-6 02:13
没法获得每一步虚频方向。只有freq任务处理的那个结构才能观看虚频模式。
作者Author: sobereva 时间: 2016-3-6 02:17
只要比较相互的能量、求差值，积分格点必须完全一致，int=ultrafine要加就都加。
其它分子倒不一定非得opt=tight，当前情况是否用tight收敛限对结构和能量的影响只是极为轻微的。不过若你是完美主义，还是都用tight为宜。
作者Author: dxyan 时间: 2016-3-6 17:38
哦哦，那就是出现了一个用了int=ultrafine，其他都要加。谢谢SOb老师。
作者Author: quora 时间: 2016-4-20 22:14
学习了，谢谢sob老师
作者Author: ayjaybfv 时间: 2016-6-29 01:59
学习了
作者Author: ter20 时间: 2016-7-5 08:35
最近在做激发态优化，用的Opt=Vtight, Int=UltraFine，有几个结构总是有一个虚频，怎么都消不掉，试了maxstep=1也不行（Sob也不推荐这个法子）。但是TD计算没有解析二阶导数，无法使用Opt=CalcAll，不知道该如何处理比较合适？是换优化算法吗？
作者Author: sobereva 时间: 2016-7-5 09:39
注意做振动分析的时候此时也得用ultrafine
拿最后结构，opt里写readfc读取freq产生的精确Hessian再优化几步
也可以尝试gdiis。
或按照虚频模式位移一点点再用calcfc并继续优化
作者Author: ter20 时间: 2016-7-5 11:04
感谢Sob老师的回复，我挨个试试。还有些问题：
1. 我试了CalcFC，会报错，说TD没有解析梯度，我想应该是您的笔误吧
2. 我看了结果文件，里面有Using GEDIIS/GDIIS optimiser这样的语句，应该是已经用了GDIIS算法了
3. 如果关键字opt=(vtight,maxstep=1) Int=UltraFine freq这样肯定是没问题的对吧？你说的Freq时也要用Int=UltraFine是指Freq和Opt分开做是吧
作者Author: sobereva 时间: 2016-7-5 12:54
1 你就freq之后readfc就完了。应该也有选项强行让高斯用数值Hessian。
2 不是，高斯默认是自动把GEDIIS、RFO通过一定规则组合在一起用的，和直接用GDIIS不同。有时候默认的优化方式不收敛，改用GDIIS能解决。
3 没问题，但opt里还应该增加notrust，否则优化过程中程序还是会自动根据情况自动改变步长上限。通常步长上限用3、5就够了，1太费时间。是。
作者Author: ter20 时间: 2016-7-5 13:15
多谢Sob老师！又学习了不少
作者Author: whtu 时间: 2016-7-11 16:29
sob老师，我也遇到了类似的情况：
我的输入文件如下：
%chk=cation172.chk
%nprocshared=24
%mem=60GB
# opt=tight freq b3lyp/6-311G(d)
cation172
1 1
C 2.14795413 1.19505188 0.11503324
C 1.42130080 -0.06121699 -0.40010118
N 0.03601540 -0.06227507 0.09171507
C -0.65757445 1.13905339 -0.39474465
C -0.65959441 1.14242752 -1.93473963
C -1.38621227 2.40096214 -2.44436410
C -1.38823228 2.40433624 -3.98435908
C -2.11485005 3.66287090 -4.49398356
C 0.03794348 -0.06549587 1.56171027
C 0.76456142 -1.32403044 2.07133474
C -0.65760824 -1.26144080 -0.40000417
C -2.10885962 -1.26254927 0.11523189
H 1.41989737 -0.05887260 -1.47009770
H 1.92615872 -0.93565334 -0.04601145
H -1.66591143 1.13828326 -0.03675591
H -0.15269190 2.01191552 -0.03682661
H -0.65901166 -1.25909641 -1.47000068
H -0.15275032 -2.13587716 -0.04591444
H -0.97039352 -0.06626604 1.91969897
H 0.54282598 0.80736626 1.91962835
H -2.61374211 -2.13541141 -0.24268618
H -2.10745620 -1.26489365 1.18522841
H -2.61371754 -0.38811292 -0.23885785
H 3.15629112 1.19582207 -0.24295548
H 1.64309619 2.06948824 -0.23905647
H 2.14935758 1.19270748 1.18502976
H 0.76596484 -1.32637484 3.14133125
H 0.25967891 -2.19689258 1.71341665
H 1.77289841 -1.32326027 1.71334604
H 0.34874256 1.14319759 -2.29272838
H -1.16447702 0.26956542 -2.29265766
H -2.39454924 2.40019209 -2.08637533
H -0.88132964 3.27382424 -2.08644610
H -0.37989531 2.40510622 -4.34234786
H -1.89311496 1.53147417 -4.34227708
C -2.11687009 3.66624499 -6.03397854
H -3.12318701 3.66210093 -4.13599477
H -1.60996736 4.53573296 -4.13606558
C -2.84348780 4.92477969 -6.54360301
H -1.10853312 3.66701491 -6.39196733
H -2.62175282 2.79338294 -6.39189652
H -2.84489134 4.92712402 -7.61359952
H -3.85182476 4.92400977 -6.18561421
H -2.33860506 5.79764173 -6.18568503
得到的log文件中，最后一个的最大受力、均方根位移等是4个YES：
Item Value Threshold Converged?
Maximum Force 0.000004 0.000015 YES
RMS Force 0.000000 0.000010 YES
Maximum Displacement 0.000041 0.000060 YES
RMS Displacement 0.000007 0.000040 YES
而最后显示有一个一个虚频，这个是因为做freq时上面的4项其实并不都是YES的么？只是没有写入到log文件？
作者Author: sobereva 时间: 2016-7-11 17:35
帖子里已经对原理分析讨论得很详细了，再仔细看。
关键是Hessian不同所导致的，优化时候是近似的，freq时候是精确的。
作者Author: 菩城 时间: 2016-7-11 19:19
好贴，感谢！
作者Author: 小丫姐姐 时间: 2017-7-14 09:56
老师好。结构优化不收敛后采用了帖子里的方法 SCF(novaracc,noincfock) int=acc2e=12 opt freq ，出现opt收敛freq后两项不收敛还差距较大的情况。http://gaussian.com/faq3/建议读取chk文件再优化，这个输入文件该如何写呢？本身计算花费了四天，如果calcall的话耗时会更多吧？感谢。
作者Author: sobereva 时间: 2017-7-14 10:07
SCF(novaracc,noincfock) int=acc2e=12这是处理SCF不收敛的，和这里说的帮助结构收敛完全是两码事：
量子化学计算中帮助几何优化收敛的常用方法
http://sobereva.com/164
作者Author: 小丫姐姐 时间: 2017-7-14 10:27
嗯嗯，理解错误。那么我现在的freq不收敛用calcall的话可能耗时太久，我尝试一下用opt=tight freq，int=ultrafine。另外，如果按照http://gaussian.com/faq3/说的读取chk文件的话输入文件该怎么写呢？菜鸟一枚，想着多试几种方法，看哪一种比较合适。
作者Author: sobereva 时间: 2017-7-14 10:41
甭管高斯页面上的，特别是对于没有判断能力的初学者来说
光从chk里读取结构什么的根本不可能解决问题
作者Author: 小丫姐姐 时间: 2017-7-14 10:44
好的，那我试试opt=tight freq,int=ultrafine.谢谢老师。
作者Author: 曼竹w 时间: 2017-9-12 18:24
sob老师，请问优化分子结构的时候出现了一个虚频，请问是提取出最后一步加opt=tight呢，还是在刚开始的输入文件中加opt=tight 呢？
作者Author: sobereva 时间: 2017-9-12 18:51
最后一步
作者Author: 曼竹w 时间: 2017-9-12 18:53
好的，谢谢sob老师。
作者Author: whtu 时间: 2017-9-21 10:47
sob老师，请问输入文件的关键词为：# opt freq b3lyp/cc-pvdz empiricaldispersion=gd3bj，log文件最后满足4个YES，但是用gview打开log文件，却发现有一个虚频，其大小为-9.01，请问这样是否可以接受？
作者Author: sobereva 时间: 2017-9-21 12:03
虽然不一定对你要研究的问题有什么明显不良影响，但还是应当消掉虚频，省得被审稿人吐槽之类。
PS：上传较大文件前请先压缩以节约论坛空间
作者Author: xylz6188 时间: 2017-9-21 13:59
sobereva老师，您好。请教一下，我用DMol3来作反应机理的时候，由于是周期性模型，偶尔会有一些小的虚频，实在消不掉，我一般把虚频小于40以下的全部忽律了，不知道这样可不可以。
作者Author: whtu 时间: 2017-9-21 14:10
好的。谢谢sob老师
作者Author: sobereva 时间: 2017-9-21 18:40
最好都消掉，否则被问到，不好解释
作者Author: xylz6188 时间: 2017-9-21 19:58
各种方法都用了，就是消不掉，MS做过渡态确实有很多问题
作者Author: 柒月小鱼 时间: 2017-9-22 10:08
刚好解决问题 谢谢sob老师
作者Author: qcn 时间: 2017-9-22 10:23
请问sob老师，A和B两分子，A在默认情况下几何收敛，B在opt=tight int=ultrafine或calcall下才收敛，随后它们在相同条件下计算单点能，算出来的单点能有可比性吗？还是A必须和B在相同条件下重新几何优化？
作者Author: sobereva 时间: 2017-9-22 15:50
两个结构的精度一般不会差多少，可忽略不计（但也有反例）。计算单点能时只需确保积分格点相同即可
作者Author: li447fan 时间: 2020-12-16 22:12
请问如果用orca计算出现虚频是不是就不能用第二种方法了？用B3LYPD3BJ/def2SVP算的上百原子的有机分子，出现2-3个虚频，应该opt时就采取严格收敛方式，还是，像http://bbs.keinsci.com/thread-2139-1-1.html评论中说的，改变原子坐标呢（这个原子振动坐标怎么加其实没看懂，在原坐标的基础上，根据振动频率加吗）
作者Author: sobereva 时间: 2020-12-17 09:31
不管什么程序，处理思想都大同小异，ORCA没有特殊性。
帖子里说了，通常只有破坏对称性的虚频才需要考虑按照虚频模式调节结构（原坐标加上正则坐标乘上人为指定的系数）
作者Author: li447fan 时间: 2020-12-18 10:36
如果用了tightopt，tightSCF，Grid5，依然有两处较小的虚频，应该继续使用更严格的verytight吗
作者Author: sobereva 时间: 2020-12-19 09:22
Grid6
如果用了RIJCOSX，也用gridx6
作者Author: li447fan 时间: 2020-12-21 00:30
社长，用了gridx6计算后，虚频反而比之前误用grid5更严重了，发现opt结束时不是全yes，是有一个no的，只是最后normal结束的。是不是应该重新更好地优化结构呢？还是绝对值小于10的虚频就不管了呢
作者Author: sobereva 时间: 2020-12-21 00:39
分清楚gridx和grid，各司其职
作者Author: li447fan 时间: 2021-3-29 10:14
在加入recalc=3之后，结构反而越来越难收敛了，基本只有1个yes或都是no，跑了快一个星期了。请问是继续等呢，还是尝试换个方法。（一开始，B3LYP-D3/6-311G**下opt四个yes，freq两个yes）
作者Author: 北大-陶豫 时间: 2021-3-29 16:52
优化的啥，能贴一下输入输出文件不
作者Author: HJC 时间: 2021-4-9 14:00
想请问老师，计算频率的时候虽然出现了两个NO，但无虚频，这样优化后的结构是否能用？谢谢您
Item Value Threshold Converged?
Maximum Force 0.000021 0.000450 YES
RMS Force 0.000003 0.000300 YES
Maximum Displacement 0.008541 0.001800 NO
RMS Displacement 0.002329 0.001200 NO
作者Author: 喵星大佬 时间: 2021-4-9 19:47
说明势能面比较平滑，力已经充分收敛了但是位移还没有，影响不大，如果介意也可以用这次振动分析得到的heseeian再优化一下，基本就OK了
作者Author: jiuyangli 时间: 2022-10-19 01:25
请问老师，我在做激发态优化时，opt的4个都是yes，但是用opt后的结构做频率计算，4个都是no，均是在同一计算级别下做的，这是怎么回事？
作者Author: sobereva 时间: 2022-10-19 01:26
freq的时候guess=read读取opt最后的波函数再试。若还有这个问题，把opt和freq的输入、输出文件都压缩后上传
作者Author: jiuyangli 时间: 2022-10-19 09:00
老师，这是guess=read重新计算的，还是这样。
(, 下载次数 Times of downloads: 2)
作者Author: sobereva 时间: 2022-10-19 18:58
opt的gjf也上传，否则都不知道你用modredundant读的是什么
作者Author: jiuyangli 时间: 2022-10-22 09:34
老师，这是重新上传后的。
作者Author: sobereva 时间: 2022-10-23 04:17
写modredundant是完全多余的
使用IOp自定义泛函的时候不要用opt freq，因为IOp没法传递到后续子任务，导致freq部分和opt部分用的不是同一个泛函
作者Author: Sui1999 时间: 2024-3-15 09:46
老师您好，我在做一个由CHONP元素组成的体系，使用orca501优化结构，关键词是! b3lyp/G cc-pVTZ d3 opt freq，计算的结果总是有虚频，使用tihgtopt，最后的结果仍然有虚频，利用“Recalc_Hess 3”则无法收敛；使用高斯用相同的关键词最后得到的结构没有虚频，拿该结构再用orca仍然会有虚频；使用cc-pvdz的话则可以收敛且没有虚频，但是使用该结构再去优化，仍然会有虚频；且上述虚频都是4个，大小和振动模式十分相似，想问一下老师这种情况有什么好的解决办法吗。
作者Author: sobereva 时间: 2024-3-15 17:35
提高积分格点精度再试
我自己基本不用ORCA做opt freq任务
cc-pVTZ给B3LYP用不是好主意，下文说了
谈谈量子化学中基组的选择
http://sobereva.com/336（http://bbs.keinsci.com/thread-3545-1-1.html）
作者Author: Sui1999 时间: 2024-3-15 18:42
感谢老师的回复，已经试过关键词“defgrid3”，已经是orca最高的格点精度了，还是会产生虚频，且虚频和之前的还是很类似。
作者Author: sobereva 时间: 2024-3-15 20:06
懒得折腾就用Gaussian吧
作者Author: wzkchem5 时间: 2024-3-15 23:33
可以试试关掉RIJCOSX。
对于杂化泛函，ORCA默认打开RIJCOSX，其中的RI部分需要辅助基组，如果用户没有给定辅助基组，ORCA会默认对非相对论基组用def2/J，对相对论基组用SARC/J。但是def2/J恐怕未必适合cc-pVTZ。高斯没有类似RIJCOSX的方法，因此不会引入这种近似，但与此同时也会算得更慢。
当然正如sob老师所说，cc-pVTZ不太适合给DFT，尤其是非双杂化泛函用，一定要用triple zeta基组优化的话，不如用def-TZVP、def2-TZVP(-f)或def2-TZVP。这样一来也就自动解决了辅助基组的问题，因为def2/J完全适合搭配这些基组使用
作者Author: qingmo 时间: 2024-6-1 20:53
本帖最后由 qingmo 于 2024-6-1 22:00 编辑
老师您好！我对一个分子结构做常规优化后，再做freq计算。没有虚频，但其中的最大位移是NO。因此我取最后结构又做了opt＝calcall计算，查看log文件结果是没有虚频，四个也都是YES。
之后我打开opt＝calcall计算后的log文件取出最后结构再做freq计算，虽然还是没虚频，但关于位移的两项却变成NO了。。。
请问老师这是什么原因，calcall不是和freq一样都是精确计算Hessian计算吗
作者Author: sobereva 时间: 2024-6-2 09:29
优化过程计算位移的算法不仅依赖于当前步的信息，还依赖于之间一些步。而且你若要100%精确基于opt最后的结构算，应该用geom=check从chk里读取全精度的坐标。
作者Author: qingmo 时间: 2024-6-2 16:26
本帖最后由 qingmo 于 2024-6-2 16:29 编辑
感谢老师您的答复！已改成geom＝cheak提取chk坐标计算。
老师您说的第一句话“依赖于之间一些步”，我不是太懂。。请问可以理解为此时体系结构和上一步优化的结构之间的关系吗？
另外，geom计算是读取了全精度坐标，如果把全精度坐标提取出来应怎么做？我在log和fchk文件里都没找到全精度的。。之前都是GV里打开out文件计算，确实菜鸟行为了
作者Author: qingmo 时间: 2024-6-3 11:41
本帖最后由 qingmo 于 2024-6-3 12:49 编辑
按照老师您所说的用geom=check从chk里读取全精度的坐标步骤后，freq得到的确实如opt=calcall步骤一样，没有出现虚频和NO了，看来应该就是坐标问题。
之后用GV观察到fchk文件坐标小数点后有效数字达八位，而原来out文件只有六位，因此复制该坐标做freq计算，此时结果又出现了NO。
不知道是因为GV观察到的fchk文件坐标其实也并不是精确坐标？还是因为即使是，但gaussian要想读取到精确的，只能根据geom=check从chk步骤才能读取到？
（然后分别用out文件坐标，fchk文件坐标，和geom=check读取chk文件计算激发能，激发能和振子强度基本一致。是不是只要没有虚频，就算位移相关项出现了NO，关于算跃迁偶极矩，激发能这类物理量，此时的结构也是可以接受的？）
作者Author: lurensan 时间: 2024-6-3 12:04
各位老师，请问这个怎样把虚频消掉，尝试了很多方法都不好使
(, 下载次数 Times of downloads: 1)
作者Author: sobereva 时间: 2024-6-3 23:51
本帖里该说的都说了，没有可补充的
作者Author: sobereva 时间: 2024-6-3 23:52
gview载入、保存时的有效位数问题。gjf里的数据记录精度比chk里的低
作者Author: qingmo 时间: 2024-6-4 13:32
好的，谢谢老师
作者Author: blue11 时间: 2024-9-5 10:13
请问sob老师，如果A和B优化过程中，适用的坐标不一致，需要该成一致的坐标吗？
如果A的不同构型A1和A2使用的坐标也不一致，需要放在一样的坐标下进行比较吗？
还是只需要一样的积分格点就可以了？谢谢sob老师！
作者Author: sobereva 时间: 2024-9-5 10:49
看不懂
A和B都不是一个物质，不可能坐标一致
作者Author: blue11 时间: 2024-9-5 19:34
不好意思，sob老师，没说明白我的意思。我想问一下，我的A在高斯默认的冗余内坐标下优化成功了一个构型，但是在笛卡尔坐标下优化成功了另一个构型。请问一下，我可以计算一下这两个构型的的能量，然后选出一个能量最稳定的吗？
作者Author: sobereva 时间: 2024-9-6 00:05
如果你是指在两种坐标系下优化得到了不同构型，应当用能量更低的（确切来说是自由能更低的）。但注意未必这就是真正能量最低结构，有可能需要用molclus做构象/构型搜索来确定真正的能量最低结构。
作者Author: 虎王 时间: 2024-10-8 18:55
请问一下，用了opt(gdiis)后可以优化收敛，但log出现了Warning! Use of Opt=GDIIS is deprecated since it is seldom a good choice. 是否有什么影响？例如用了opt(gdiis)是否会导致优化构型不精确？
作者Author: zjxitcc 时间: 2024-10-8 18:58
本帖最后由 zjxitcc 于 2024-10-8 18:59 编辑
如果正常结束，没有影响，不会不精确。可以google了解一下优化方法的基本原理，便不会有此困惑。
作者Author: 虎王 时间: 2024-10-8 19:04
好的。非常感谢。
作者Author: sobereva 时间: 2024-10-9 04:20
这个warning是多余的
默认的GEDIIS很多时候并不如GDIIS
能决定最终精度的是几何优化收敛限，而不是优化算法
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |