计算化学公社

标题: 关于分子有假虚频时的热力学校正问题 [打印本页]

作者
Author: wzkchem5 时间: 2020-12-9 13:13
标题: 关于分子有假虚频时的热力学校正问题
最近想到一个问题，我觉得是很多做计算的人忽视的，想和大家讨论一下，看看大家有没有什么好的意见。
众所周知，一个分子如果计算时出现特别小的虚频（比如-50 cm-1以内），可能不一定代表这个分子真的有虚频，没准只是计算的数值误差（格点误差，SCF或几何优化的收敛阈值不够严，以及其他的各种prescreening的阈值）导致的。尤其是比如说通过把分子在虚频的方向上做扰动的方法，怎么消虚频都消不掉，那么这时一般认为（或者说怀疑）这个虚频是假虚频。当然严格来说，这时应该加大格点，把某些重要的阈值设严，直到没有虚频为止。但是这个方法不一定总是可行的，因为：
（1）量化程序的某些阈值可能是没办法改变的；
（2）加大格点、阈值设严以后，计算时间增加，有时增加的开销可能不能接受；
（3）有时候一个研究已经基本做完了，只有最后一个分子的假虚频总是消不干净，如果加大格点就意味着要把之前的计算全部用大格点重做一遍，会觉得得不偿失。
在这样的背景下，有些人（比如说ORCA团队）会建议“带瘤生存”，也就是如果有充分的理由认为一个频率是假虚频，那么可以把它看作实频，不去管它。但是这样也有问题，因为一般程序在计算热力学校正的时候，会忽略虚频的贡献。假如这个虚频是某个过渡态的真虚频，用户要算活化能/活化焓/活化Gibbs自由能，那么这样做是正确的，也是统计热力学理论要求的。但是如果这个虚频是假虚频，这样做就有问题了。
比如我有一个特别小的实频，几乎等于0 cm-1。假设我把这个实频算准了，并且用Grimme的QRRHO方法计算它在298 K下对Gibbs自由能的贡献，那么可以算出来贡献是-2.68 kcal/mol。这个数比0小很多，也很好理解，因为RRHO模型下0 cm-1的振动频率对Gibbs自由能的贡献是负无穷，QRRHO在频率很小的时候把RRHO的Gibbs自由能校正damp掉了一些，所以等于一个负的有限值。现在假设我把这个频率算成了-1 cm-1，结果它立刻就对Gibbs自由能校正没有贡献了，这就引入了多达2.68 kcal/mol的误差，都跟DFT的泛函误差差不多了，如果分子有不止一个假虚频，误差可以多达5~10 kcal/mol，后果是灾难性的。等于一个无穷小的频率误差就会导致一个并非无穷小，而且还算是比较大的Gibbs自由能误差。
那么我现在的一个想法就是，假如已经基本确定了某个虚频是假虚频，但是又没办法把它消掉，那么能不能直接把算出来的Gibbs自由能减去2.68 kcal/mol，也就是假设那个假虚频实际上是0 cm-1？这样做的考虑是，越接近0的频率越容易被错误地算成虚频，所以如果已知一个频率被算成了虚频，那么某种意义上讲，它最可能的值就是0 cm-1。我最近算的一个体系上测试结果还不错，算一个分子量200左右的有机分子在石墨上的吸附，用团簇模型描述石墨，计算发现有机分子没有虚频，石墨基底有1个小虚频，吸附了有机分子的石墨有4个小虚频，而且都消不掉，或者越消越多。如果直接算吸附Gibbs自由能，发现在-2 kcal/mol左右，这明显是不可能的，如果分子量200的分子的吸附都这么弱，那活性炭就净化不了空气了。做了校正（-3*2.68 kcal/mol）以后，吸附Gibbs自由能降到-10 kcal/mol，明显合理了。
不知道大家觉得这种校正方法是否合理，以及有没有文献研究这方面的？

作者
Author: sobereva 时间: 2020-12-9 14:03
ORCA团队的带瘤生存的主张，在我来看是程序这方面做得不济，所以给用户这么一个歪招。换做是Gaussian的开发者我估计没有赞同的。

以前我也经常遇到特别微小的虚频，都是势能面超级平缓的情况，每次我都是无一例外硬着头皮消掉，哪怕是wB97XD这样对积分格点不是特敏感的泛函也往往要上到int=superfine才能解决。

把微小虚频视为0频极限的做法，确实有道理，但如果碰上较真的审稿人，可能最终说服不了。另外，有一个微小虚频，也暗示着其它非常低的频率可能也是不精确的，可能对自由能热校正量贡献的误差也挺大。

作者
Author: biogon 时间: 2020-12-9 14:52
假设微小虚频为0频率极限确实是个不错方法，但是怎么说服审稿人是个问题，毕竟没有证明过这个问题

作者
Author: wzkchem5 时间: 2020-12-9 17:39
本帖最后由 wzkchem5 于 2020-12-9 21:04 编辑

sobereva 发表于 2020-12-9 14:03
ORCA团队的带瘤生存的主张，在我来看是程序这方面做得不济，所以给用户这么一个歪招。换做是Gaussian的开发 ...

确实，有小虚频可能意味着实频也不准，但是我怀疑这个误差应该一般比QRRHO本身的误差要小。比如QRRHO里有个参数omega0 = 100 cm-1，是比较arbitrary地取的，Grimme说这个参数在50~150 cm-1范围内结果都差不多，100 cm-1的取值没有什么理论上的理由，也不是拟合出来的。那么现在假设有一个频率是50 cm-1，如果按omega0 = 100 cm-1算自由能校正量，结果跟频率为90 cm-1、按omega0 = 50 cm-1算的自由能校正量差不多。所以如果omega0取100 cm-1跟取50 cm-1都可以接受的话，那也就意味着90 cm-1附近的振动频率算成50 cm-1都可以接受，只要误差是系统的。
而且我其实觉得QRRHO的数学形式也很有问题，比如0 cm-1的频率的自由能校正是-2.68 kcal/mol，但是1 cm-1的校正量是-1.61 kcal/mol。也就是有极小的概率，即使你把每个频率都算准到1 cm-1以内，你也可能有1 kcal/mol的误差（尤其频率越小，算准就越困难，因为频率涉及把mass-weighted Hessian的本征值开根号，所以频率越小就对Hessian的误差越敏感）。直觉上很难接受0 cm-1和1 cm-1的热力学校正量有这么大差别。Grimme的数学形式是假设频率很小的时候，把这个振动看作一个频率相同的转子的转动，但是即使这样计算结果也是以对数速度发散的（和RRHO一样），所以又人为把这个转子的转动惯量μ改成μBav/(μ+Bav)，其中Bav=10^-44 kg m^2，这样振动频率趋于0的时候，转动惯量趋于Bav而不是趋于无穷大，这样0 cm-1的校正量就不是负无穷了。但是有意思的是，只要当频率小于0.03 cm-1的时候，μ才开始比Bav大，所以在频率大于0.1 cm-1的时候，QRRHO的自由能校正值表现得跟一个对数发散的函数基本没有区别。但是一个人要做多少计算才能碰上一次有频率（实频）小于0.1 cm-1的情况呢？我自己反正是没碰到过的。所以我觉得QRRHO还是没有真正解决低频自由能校正发散的问题，它本质上只是把发散速度减慢了，然后在频率非常接近0、一个人做一辈子计算不一定碰得上一次的区域，把发散性regularize掉，但是这个regularization对于大部分不涉及超级小的频率的计算都和没regularize一样。

要是我来开发QRRHO方法的话，第一我会让低频时的自由能校正随频率变化慢一些，第二我会让自由能校正是频率平方的解析函数，方便解析延拓到虚频。

作者
Author: QuantumicGuy 时间: 2022-6-7 10:20
最近用GFN1/(2)-xTB来优化一个机械互锁轮烷分子（442个原子，C、H、O、N、P、Pt），级别从normal到extreme都试了，总是有4个-2~ -5 cm-1的小虚频，这是一篇实验组的文章，目的是去吻合分子在电镜下的不同构象，这样的小虚频对于此研究目的是否可以忽略呢？

作者
Author: dantevinsky 时间: 2022-6-7 11:20

wzkchem5 发表于 2020-12-9 17:39
确实，有小虚频可能意味着实频也不准，但是我怀疑这个误差应该一般比QRRHO本身的误差要小。比如QRRHO里有 ...

说实话，无论是omega还是Grimme的那个SRRHO，都是一种简易的近似计算方法。本来计算精度就在那里，再在他的模型上添东西也难以摆脱他们简易近似的本质。低频矫正从上世纪以来就有一堆的方法可以做，如果需要更高精度的热力学量，就不应该用这两个方法。
话说，我在看到上世纪90年代的文章中Truhlar就提出过SRRHO的渐进策略[1]，只不过没有进一步阐述罢了。。

[1] 10.1002/jcc.540120217

作者
Author: wzkchem5 时间: 2022-6-7 17:53

QuantumicGuy 发表于 2022-6-7 03:20
最近用GFN1/(2)-xTB来优化一个机械互锁轮烷分子（442个原子，C、H、O、N、P、Pt），级别从normal到extreme ...

问题应该不太大，不过可以考虑减小数值Hessian的步长，因为xTB没有解析Hessian，是用解析梯度做数值差分得到的数值Hessian，可能会引入一些数值噪音

作者
Author: wzkchem5 时间: 2022-6-7 18:22

dantevinsky 发表于 2022-6-7 04:20
说实话，无论是omega还是Grimme的那个SRRHO，都是一种简易的近似计算方法。本来计算精度就在那里，再在他 ...

对，但是即便是近似方法和近似方法比，有的方法在一些corner case的情况下也会有问题。比如Grimme这个方法，如果用来做VTST之类的计算，而且IRC上某些地方在垂直于反应坐标的方向上有虚频的话，就会发现振动配分函数在这个虚频转换为实频的那一点不可导，或者说导数趋于无穷大。这个是因为频率的平方是坐标的连续可导函数，为了让振动配分函数也是坐标的连续可导函数，振动配分函数在频率等于0处必须是频率的平方的连续可导函数，但是Grimme的公式在频率为0处只是频率的连续可导函数（而且还只是右可导），而不是频率的平方的连续可导函数。造成的后果就是在FEP上出现一个artificial minimum，而且这个minimum的底部是个cusp，不可导。
相比之下Truhlar的这个公式就要好得多，在频率为0处是频率平方的连续可导函数，因此要physical得多。唯一问题是需要计算有效转动惯量才能用，不过也可以仿照Grimme的做法，用一个固定的数作为有效转动惯量。

作者
Author: dnlx 时间: 2023-1-29 09:40
Grimme他们组去年的文章（/10.1002/anie.202205735）提到了以前的一种做法，就是直接乘-i变成实频。他们没有评价这种方法的优劣。你认为这种做法是否可以接受？

“Note that it is quite typical for large systems to have a few low-energy imaginary frequencies <50–100 cm-1, which are often technically related to a too small DFT integration grid. This is per se not a problem if properly dealt with, that is, the frequencies should be inverted (multiplied by -i) and treated as normal real frequencies, like in the xtb program using the mRRHO approach.”

作者
Author: wzkchem5 时间: 2023-1-29 17:30

dnlx 发表于 2023-1-29 02:40
Grimme他们组去年的文章（/10.1002/anie.202205735）提到了以前的一种做法，就是直接乘-i变成实频。他们没 ...

我觉得不如处理成0合理，因为如果这个虚频确实是因为数值误差把实频错算成虚频的话，那么虚频的绝对值越大，原来的那个实频就越可能很小。而按照Grimme的做法，就是虚频的绝对值越大，对应的实频越大，这个我感觉无论如何都是不合理的。

作者
Author: mfdsrax2 时间: 2023-1-29 18:15
我不懂这个，只是咨询一下

一个几乎等于0 cm-1的虚频，可以算出来贡献是-2.68 kcal/mol。如果把这个频率算成了-1 cm-1，结果它立刻就对Gibbs自由能校正没有贡献了

这是为什么？

作者
Author: wzkchem5 时间: 2023-1-29 18:25

mfdsrax2 发表于 2023-1-29 11:15
我不懂这个，只是咨询一下

一个几乎等于0 cm-1的虚频，可以算出来贡献是-2.68 kcal/mol。如果把这个频率 ...

如果这个分子有且仅有一个虚频，那么一般程序都会认为这个结构是过渡态，从而改用过渡态的配分函数公式，因此虚频对配分函数没有贡献。所以本质上是当有虚频时，配分函数的公式发生了突变。
如果虚频数目不止一个（例如算过渡态结果多了一个假虚频，或者优化极小值点结果有两个假虚频），据我所知没法严格推导出配分函数公式，因为这种情况下配分函数怎么定义都是一个问题。但是一般程序仍然会套用计算过渡态的逻辑，把所有虚频对配分函数的贡献忽略，其他频率的贡献照常算。

作者
Author: dnlx 时间: 2023-1-29 20:47

wzkchem5 发表于 2023-1-29 17:30
我觉得不如处理成0合理，因为如果这个虚频确实是因为数值误差把实频错算成虚频的话，那么虚频的绝对值越 ...

是的，如果出现几个小虚频，消虚频的过程中可以发现，它们是“平移”而不是“翻转”到实频。

欢迎光临计算化学公社 (http://bbs.keinsci.com/)