计算化学公社

 找回密码 Forget password
 注册 Register
Views: 7293|回复 Reply: 12

[综合交流] 关于分子有假虚频时的热力学校正问题

[复制链接 Copy URL]

6693

帖子

0

威望

4012

eV
积分
10705

Level 6 (一方通行)

发表于 Post on 2020-12-9 13:13:12 | 显示全部楼层 Show all |阅读模式 Reading model
最近想到一个问题,我觉得是很多做计算的人忽视的,想和大家讨论一下,看看大家有没有什么好的意见。
众所周知,一个分子如果计算时出现特别小的虚频(比如-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,明显合理了。
不知道大家觉得这种校正方法是否合理,以及有没有文献研究这方面的?
BDF(https://bdf-manual.readthedocs.io/zh_CN/latest/Introduction.html)、ORCA(https://orcaforum.kofo.mpg.de/index.php)开发团队成员

4万

帖子

99

威望

4万

eV
积分
89946

管理员

公社社长+计算化学玩家

发表于 Post on 2020-12-9 14:03:13 | 显示全部楼层 Show all
ORCA团队的带瘤生存的主张,在我来看是程序这方面做得不济,所以给用户这么一个歪招。换做是Gaussian的开发者我估计没有赞同的。

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

把微小虚频视为0频极限的做法,确实有道理,但如果碰上较真的审稿人,可能最终说服不了。另外,有一个微小虚频,也暗示着其它非常低的频率可能也是不精确的,可能对自由能热校正量贡献的误差也挺大。
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班,内容介绍以及往届资料购买请点击链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的最佳途径。培训各种常见问题见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

1158

帖子

1

威望

2801

eV
积分
3979

Level 5 (御坂)

发表于 Post on 2020-12-9 14:52:59 | 显示全部楼层 Show all
假设微小虚频为0频率极限确实是个不错方法,但是怎么说服审稿人是个问题,毕竟没有证明过这个问题

6693

帖子

0

威望

4012

eV
积分
10705

Level 6 (一方通行)

 楼主 Author| 发表于 Post on 2020-12-9 17:39:50 | 显示全部楼层 Show all
本帖最后由 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方法的话,第一我会让低频时的自由能校正随频率变化慢一些,第二我会让自由能校正是频率平方的解析函数,方便解析延拓到虚频。

BDF(https://bdf-manual.readthedocs.io/zh_CN/latest/Introduction.html)、ORCA(https://orcaforum.kofo.mpg.de/index.php)开发团队成员

68

帖子

0

威望

841

eV
积分
909

Level 4 (黑子)

#未来可寄

发表于 Post on 2022-6-7 10:20:36 | 显示全部楼层 Show all
最近用GFN1/(2)-xTB来优化一个机械互锁轮烷分子(442个原子,C、H、O、N、P、Pt),级别从normal到extreme都试了,总是有4个-2~ -5 cm-1的小虚频,这是一篇实验组的文章,目的是去吻合分子在电镜下的不同构象,这样的小虚频对于此研究目的是否可以忽略呢?

44

帖子

0

威望

648

eV
积分
692

Level 4 (黑子)

发表于 Post on 2022-6-7 11:20:01 | 显示全部楼层 Show all
wzkchem5 发表于 2020-12-9 17:39
确实,有小虚频可能意味着实频也不准,但是我怀疑这个误差应该一般比QRRHO本身的误差要小。比如QRRHO里有 ...

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

[1] 10.1002/jcc.540120217

6693

帖子

0

威望

4012

eV
积分
10705

Level 6 (一方通行)

 楼主 Author| 发表于 Post on 2022-6-7 17:53:01 | 显示全部楼层 Show all
QuantumicGuy 发表于 2022-6-7 03:20
最近用GFN1/(2)-xTB来优化一个机械互锁轮烷分子(442个原子,C、H、O、N、P、Pt),级别从normal到extreme ...

问题应该不太大,不过可以考虑减小数值Hessian的步长,因为xTB没有解析Hessian,是用解析梯度做数值差分得到的数值Hessian,可能会引入一些数值噪音
BDF(https://bdf-manual.readthedocs.io/zh_CN/latest/Introduction.html)、ORCA(https://orcaforum.kofo.mpg.de/index.php)开发团队成员

6693

帖子

0

威望

4012

eV
积分
10705

Level 6 (一方通行)

 楼主 Author| 发表于 Post on 2022-6-7 18:22:18 | 显示全部楼层 Show all
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的做法,用一个固定的数作为有效转动惯量。
BDF(https://bdf-manual.readthedocs.io/zh_CN/latest/Introduction.html)、ORCA(https://orcaforum.kofo.mpg.de/index.php)开发团队成员

111

帖子

0

威望

597

eV
积分
708

Level 4 (黑子)

发表于 Post on 2023-1-29 09:40:47 | 显示全部楼层 Show all
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.”

6693

帖子

0

威望

4012

eV
积分
10705

Level 6 (一方通行)

 楼主 Author| 发表于 Post on 2023-1-29 17:30:09 | 显示全部楼层 Show all
dnlx 发表于 2023-1-29 02:40
Grimme他们组去年的文章(/10.1002/anie.202205735)提到了以前的一种做法,就是直接乘-i变成实频。他们没 ...

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

评分 Rate

参与人数
Participants 1
eV +5 收起 理由
Reason
dnlx + 5 谢谢

查看全部评分 View all ratings

BDF(https://bdf-manual.readthedocs.io/zh_CN/latest/Introduction.html)、ORCA(https://orcaforum.kofo.mpg.de/index.php)开发团队成员

256

帖子

0

威望

964

eV
积分
1220

Level 4 (黑子)

发表于 Post on 2023-1-29 18:15:34 | 显示全部楼层 Show all
我不懂这个,只是咨询一下

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

这是为什么?

6693

帖子

0

威望

4012

eV
积分
10705

Level 6 (一方通行)

 楼主 Author| 发表于 Post on 2023-1-29 18:25:19 | 显示全部楼层 Show all
mfdsrax2 发表于 2023-1-29 11:15
我不懂这个,只是咨询一下

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

如果这个分子有且仅有一个虚频,那么一般程序都会认为这个结构是过渡态,从而改用过渡态的配分函数公式,因此虚频对配分函数没有贡献。所以本质上是当有虚频时,配分函数的公式发生了突变。
如果虚频数目不止一个(例如算过渡态结果多了一个假虚频,或者优化极小值点结果有两个假虚频),据我所知没法严格推导出配分函数公式,因为这种情况下配分函数怎么定义都是一个问题。但是一般程序仍然会套用计算过渡态的逻辑,把所有虚频对配分函数的贡献忽略,其他频率的贡献照常算。
BDF(https://bdf-manual.readthedocs.io/zh_CN/latest/Introduction.html)、ORCA(https://orcaforum.kofo.mpg.de/index.php)开发团队成员

111

帖子

0

威望

597

eV
积分
708

Level 4 (黑子)

发表于 Post on 2023-1-29 20:47:45 | 显示全部楼层 Show all
wzkchem5 发表于 2023-1-29 17:30
我觉得不如处理成0合理,因为如果这个虚频确实是因为数值误差把实频错算成虚频的话,那么虚频的绝对值越 ...

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

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-6 05:32 , Processed in 0.209846 second(s), 24 queries .

快速回复 返回顶部 返回列表 Return to list