如何看待《关于高斯计算精度参数的常见误解》这篇文章的观点

北大-陶豫

关于高斯计算精度参数的常见误解https://mp.weixin.qq.com/s/EC3ZO8oz5A6UnQE9_Q99yA
文中的观点可取吗？实际计算中还有哪些要注意的呢？

sobereva

这篇文章作者有明显的误区，希望此文的读者都能看到以下笔者的comment避免误入歧途。（笔者comment的是那篇文章2022-01-15 09:18发布的版本，不晓得之后会不会再被修改）

以下文字太长不看版：
G16写int=ultrafine关键词完全是多余的
没特殊情况不要随便用opt=tight

-----------------

文中说“手动添加int=ultrafine关键词可以提高频率计算的精度，对G16来说并非冗余关键词”这是完全错误的。

首先要知道DFT（这里不考虑双杂化泛函）计算振动频率、极化率等等>=2阶能量导数属性，需要求解CPHF方程，这对KS-DFT计算来说也被特称为CPKS方程，这用来得到轨道展开系数对扰动的导数。求解CPKS方程在Gaussian里专门有一套格点，可以叫做“CPHF格点”。这和我们平时说的DFT积分格点（下文称作“DFT格点”）可以相同也可以不同。DFT格点影响能量、受力、激发能等等各种涉及到DFT计算的任务的结果（和耗时），而CPHF格点仅仅影响求解CPHF方程的求解结果（和耗时）。

G09里默认的是：DFT格点用fine，CPHF格点用coarse
G16里默认的是：DFT格点用ultrafine，CPHF格点用SG1

文中说的coarse < SG1 < fine < ultrafine < superfine关系并没错，属于Gaussian用户的常识性知识。也因此，从积分格点角度来说，G16算的各种与DFT有关的问题的数值层面的精度比G09要强，具体强多少，看体系、基组和泛函对格点的敏感度。比如当弥散函数特别特别多，或体系势能面特别平缓，或使用M06系列等对积分格点很敏感的泛函的时候，由于G09的格点略廉价，G16的结果可能比G09有不可忽视的改进。而文中不加限定条件地说“G09算频率的格点是远远不够的”这是大错特错，例如用B3LYP/def2-SVP算柔性不是很强的体系的opt freq任务，G09和G16的差别甚微。读者可以自行测试。

文中作者没有认识到CPHF求解对格点的敏感性本来就没有计算能量用的DFT格点那么高。这有点类似于opt freq任务用的计算级别可以比专门算单点时要低一两个档次（见《浅谈为什么优化和振动分析不需要用大基组》http://sobereva.com/387）。Gaussian的开发者不是傻子，无论是G09还是G16，默认的CPHF格点都比DFT格点低两个档正是出于这种考虑。而且作为非常流行的商业程序的开发者，必然也用很多体系做过专门的测试以保证默认设置下的数值精度对绝大多数情况来说都是充分的。

文中作者貌似嫌G16用的ultrafine/SG1这种格点搭配还不够，似乎以为用了int=ultrafine关键词的时候CPHF格点也会升到ultrafine去，并因此觉得这可以提升频率计算精度，这是明显错误的。实际上G16写不写int=ultrafine，对CPHF用的积分格点都是SG1。如果需要直接指定CPHF格点成为比如fine，则需要写CPHF=grid=fine。G16默认的SG1格点对于绝大多数情况已经足够好了，仅当对于前述的“弥散函数特别特别多，或体系势能面特别平缓，或使用M06系列等对积分格点很敏感的泛函”的时候，刻意用CPHF=grid=fine可能对改进数值精度有可查觉的益处，对于解决CPHF方程迭代求解难收敛的问题奏效几率也非常高。根据我以往的经验，没有任何情况需要用非常贵的CPHF=grid=ultrafine。

文中作者对比G16默认情况和加上int=ultrafine的情况在M062X/def2TZVP下计算环戊烷构象异构化的过渡态的opt freq任务，发现两个情况算的频率有细微偏差，得到结论说G16写int=ultrafine也会影响结果，在我来要么是弄错了数据，要么只是极特定情况下的巧合，例如此文说的情况《数值误差对计算化学结果重现性的影响》（http://sobereva.com/88）。我在G16 B.01和G16 C.01下皆未发现G16写int=ultrafine会影响振动频率计算结果的现象。我还专门亲自用G16 C.01在加和不加int=ultrafine的情况下在M062X/def2TZVP级别下对环戊烷构象异构化的过渡态做了opt freq任务，输出文件如下，可见频率根本没有丝毫差别。

opt freq.rar (161.92 KB, 下载次数 Times of downloads: 15)

2022-1-15 16:45 上传 Uploaded

点击下载
Click to download

有趣的是，作者在末尾说写上int=ultrafine“几乎不增加任何计算量”，这不仅不是说明加这个关键词有益，而恰恰说明加这个关键词根本不影响G16默认用的格点，要不然怎么可能不增加计算量？如果真是加了这个关键词既不增加计算量还能改进数值精度，都是内行的Gaussian开发者干嘛不把这当做默认设置？

文中说“考虑到积分在freq任务的耗时占比很低”这也是错误的，实际对比一下不同DFT格点、不同CPHF格点下的耗时便知。G09用默认的int=fine和int=ultrafine，我以前实测用M06-2X/def-TZVP对多巴胺小分子做freq任务耗时能相差一倍。

此文文末作者说“笔者的建议是无论是G09还是G16的用户，计算时都加上opt=tight以及int=ultrafine这两个关键词”。这个建议切勿效仿。前面我已经明确说了G16用户写int=ultrafine是完全多余的，而且实际上对于势能面不是特别平缓、用的只是B3LYP、PBE0、CAM-B3LYP等对积分格点要求并不高的泛函，用ultrafine档次的DFT格点都完全没必要，因此在很多情况下我用G16的时候还刻意写上int=fine来节约时间（但不建议初学者效仿）。

至于作者建议加上opt=tight，这更是绝对不可取，这个关键词仅当该加的时候才应当加。比起默认收敛限，用tight需要花费很多步才能达到收敛，大大增加耗时（绝不是一星半点），而绝大多数情况比收敛到默认收敛限没有任何可查觉的益处；实际上在默认收敛限下，对于绝大多数体系，最后几步肉眼都根本看不出来差异（在gview里快速切换最后几帧，一般原子球在屏幕上也就有几个像素的位置改变而已）。而且，在默认收敛限下，尚且许多人日常计算中还频繁碰上几何优化不收敛，用了tight显然更是严重雪上加霜，几何优化难收敛的问题出现概率会暴增，为了能达到tight收敛限可能用户需要花极大精力来回折腾，花费巨量时间。真正需要用tight的情况是优化完了有虚频，这在《Gaussian中几何优化收敛后Freq时出现NO或虚频的原因和解决方法》（http://sobereva.com/278）我专门说了，以及势能面极其平缓的时候（如《静电效应主导了氢气、氮气二聚体的构型》http://sobereva.com/209一文研究中做的H2二聚体研究）。Gaussian的优化收敛规则是受力小于收敛限的100倍即便位移还没收敛也正常结束任务，但势能面极其平缓的时候可能达到这个标准时此时位移收敛还不够充分，这种情况可以靠tight要求让结构也能收敛得比较准确。文中作者还说“...初始坐标没有进行对称化，在默认的收敛标准下分子因为过于柔性导致肉眼可见的偏离对称结构，分子轨道系数中原本对称的部分的贡献有明显的不同，进而导致激发态的计算出现问题。后来将初始结构做了对称化后就没问题了”，这例子完全不能体现用tight收敛限带来的益处，体系有对称性时必须做对称化再计算，这属于最基本的常识，用tight收敛限只是令不恰当的计算暴露的问题没那么严重而已，而并不解决问题。有对称性的体系计算时不考虑对称性即便用tight收敛限也救不回来，依然可能导致定性错误。例如《我对一篇存在大量错误的J.Mol.Model.期刊上的18碳环研究文章的comment》（http://sobereva.com/584）里我提到的计算，作者计算碳环体系没考虑对称性问题，导致用优化后的结构算ECD出现了根本不该有的吸收峰。若初始结构不满足实际对称性，哪怕优化用了tight收敛限照样还是会有虚假的吸收峰出现（只不过转子强度会小得多而已）。

文中说“因此G09算频率的格点是远远不够的，反而是算能量的格点看起来还可以接受...如果读者遇到了明确对于电子能来说fine格点明显不够的体系，欢迎告知笔者”。对积分格点测试的文章有的是，比如Mol. Phys., 115, 2315 (2017)里给出的下图展现了对不同泛函在不同积分格点下算弱相互作用和原子化能的积分格点误差。int=fine对应的是剪裁后的(75,302)，可见对于不少泛函称不上理想，尤其是明尼苏达系列泛函（注：Mol Phys这篇文章用的是Q-Chem，积分格点分布和Gaussian的必然有一定出入，仅供参考）



文中“较小的频率对自由能的校正更大”这种说法不准确。比如0K下，自由能热校正量就是等于ZPE而已，显然频率越小对ZPE的贡献越小。只能说较小的频率对升温导致焓的改变量、对熵的影响更大。

文中“最小频率”用词不正确，应当是“虚频”。

总之，此文的作者给出建议时应该更谨慎些，多仔细调查一下相关计算原理，以免误导他人。


PS：造谣动动嘴，辟谣跑断腿...

liyuanhe211

G16的 opt freq 任务中，加/不加 int=ultrafine 对 freq 步骤的计算路径唯一影响是多了 3/75=-5，而 IOp(3/75) 本来默认就是-5，所以加/不加 int=ultrafine 不可能对 freq 的设定有任何影响。

算比较复杂的有机分子反应无脑默认加opt=tight纯属找事儿。

wzkchem5

sobereva 发表于 2022-1-15 09:45
这篇文章作者有明显的误区，希望此文的读者都能看到以下笔者的comment避免误入歧途。（笔者comment的是那篇 ...

量子化学公众号的作者回复了
https://mp.weixin.qq.com/s?__biz ... &wx_header=0#rd

量化小菜鸡

wzkchem5 发表于 2022-1-17 22:17
量子化学公众号的作者回复了
https://mp.weixin.qq.com/s?__biz=MzU5NjMxNjkzMw==&mid=2247486728&idx=1 ...

我在想要Gaussian的默认收敛限要是都不行，那ORCA的结果不是按这个逻辑更差。。。

wzkchem5

量化小菜鸡 发表于 2022-1-17 15:48
我在想要Gaussian的默认收敛限要是都不行，那ORCA的结果不是按这个逻辑更差。。。

我看了一下，ORCA的位移收敛限（MAX 0.004, RMS 0.002）确实比高斯（MAX 0.0018, RMS 0.0012）更粗糙，但是梯度收敛限（MAX 0.0003, RMS 0.0001）比高斯（MAX 0.00045, RMS 0.0003）更严，此外ORCA还额外要求能量变化小于0.000005。所以不好说高斯和ORCA的几何优化收敛限哪个更差。
我记得ORCA 5.0推出来之前我们这边有人做过一个测试，几百个分子用默认设置做opt freq，只有一个分子有假虚频

量化小菜鸡

wzkchem5 发表于 2022-1-17 23:14
我看了一下，ORCA的位移收敛限（MAX 0.004, RMS 0.002）确实比高斯（MAX 0.0018, RMS 0.0012）更粗糙，但 ...

好的老师，我印象之前ORCA的TightOpt才对应高斯默认opt的MaxD要求。SCF来说ORCA就更低了，VeryTightSCF对应的才是高斯的Conver=7。是不是我的信息没有更新。

wzkchem5

量化小菜鸡 发表于 2022-1-17 16:20
好的老师，我印象之前ORCA的TightOpt才对应高斯默认opt的MaxD要求。SCF来说ORCA就更低了，VeryTightSCF对 ...

这个不好比，因为ORCA默认的检查SCF收敛的方法是看总能量和单电子能量，当开启DIIS的时候是看DIIS commutator，当开启SOSCF的时候是看轨道梯度，而高斯是看密度矩阵。
你有兴趣的话可以用高斯和ORCA跑同一个体系，看最后一步的能量变化、密度矩阵的RMS和MAX变化，看看到底差多少。我们因为众所周知的原因没法用高斯，没办法自己来做这个事情

量化小菜鸡

wzkchem5 发表于 2022-1-18 00:10
这个不好比，因为ORCA默认的检查SCF收敛的方法是看总能量和单电子能量，当开启DIIS的时候是看DIIS commut ...

好的，谢谢老师。

sobereva

wzkchem5 发表于 2022-1-18 00:10
这个不好比，因为ORCA默认的检查SCF收敛的方法是看总能量和单电子能量，当开启DIIS的时候是看DIIS commut ...

我经常发现ORCA基于DIIS commutator判断收敛的时候，哪怕用strongSCF甚至tightSCF，达到收敛时能量收敛精度都还不是很高，容易心里没底。我觉得要是始终也能判断能量收敛精度就好了，毕竟这是算单点时候特别关心的，比如strongSCF时至少<1E-6 Hartree。

sobereva

量化小菜鸡 发表于 2022-1-17 22:48
我在想要Gaussian的默认收敛限要是都不行，那ORCA的结果不是按这个逻辑更差。。。

Gaussian的几何优化默认收敛限对于一般目的已经足够了，在我来看是得当的。如果有人鼓吹Gaussian默认收敛限太宽，甚至觉得tight才应该是默认的，则是非常偏激、很不成熟的，对各种类型体系都做过很多研究的经验丰富的研究者一般是不会认同的。


题外话：作为程序开发者不能无脑地把收敛限设得过严，比如非要试图保证所有情况下虚频严格清零什么的，而必须权衡考虑绝大多数应用场景下结果的精度、耗时、达到收敛的难度。如果把比如int=superfine、opt=tight当成默认关键词，虽然会给极个别用户的极个别研究带来好处（比如对柔性非常强的体系出现虚频的概率可能降低为原先的<1/10），却会给绝大多数用户的体验带来显著的坏处（平均耗时和遇到几何优化不收敛的情况增加N倍）。所以对于收敛限这个问题，在开发程序，以及写经验文章给出建议的时候，切不可光从自己的角度看待问题，更得从绝大多数程序使用者的角度全面、理性看待。

mfdsrax2

作者更新了《Reply to comment to “关于高斯计算精度参数的常见误解”》
https://mp.weixin.qq.com/s/kPdcqZ4qFbn1YTeLRUiS4w

liyuanhe211

mfdsrax2 发表于 2022-1-18 09:38
作者更新了《Reply to comment to “关于高斯计算精度参数的常见误解”》
https://mp.weixin.qq.com/s/kPd ...

相比于“Reply to comment to”，这叫“Correction of”更合理

冰释之川

liyuanhe211 发表于 2022-1-18 12:58
相比于“Reply to comment to”，这叫“Correction of”更合理

233333

wzkchem5

sobereva 发表于 2022-1-18 01:53
我经常发现ORCA基于DIIS commutator判断收敛的时候，哪怕用strongSCF甚至tightSCF，达到收敛时能量收敛精 ...

从ORCA 5.0.3开始，UHF/UKS的DIIS收敛精度会变高一些，因为以前是用alpha和beta的commutator的和来判断收敛，所以在两个commutator异号的时候收敛会比较粗糙，现在改成了两者都必须收敛