计算化学公社

标题: 关于如何导出柔性扫描结果的问题 [打印本页]

作者
Author: 594yy 时间: 2021-11-30 11:16
标题: 关于如何导出柔性扫描结果的问题
本人现在对二聚水分子关于某一键长进行柔性扫描，扫描后想要导出各个扫描点的电子能（单点能），熵，吉布斯自由能等。按照sob老师在帖子（http://sobereva.com/474）中导出方法，只能导出Total Energy，这个Total Energy是指代什么能量？电子能（单点能），熵，吉布斯自由能等应当如何导出？望各位老师不吝赐教，谢谢！

作者
Author: wzkchem5 时间: 2021-11-30 16:15
熵、吉布斯自由能只在梯度为0的点（平衡位置、过渡态）处有定义，所以只能输出电子能（也就是total energy），不能输出每一点的熵、吉布斯自由能。
一定需要自由能的话，只能算potential of mean force，但是很多软件不支持这个

作者
Author: 594yy 时间: 2021-11-30 20:34

wzkchem5 发表于 2021-11-30 16:15
熵、吉布斯自由能只在梯度为0的点（平衡位置、过渡态）处有定义，所以只能输出电子能（也就是total energy ...

那可否将每个键长下的scan结果输出，再进行freq计算，得到熵、吉布斯自由能的值？

作者
Author: wzkchem5 时间: 2021-11-30 20:42

594yy 发表于 2021-11-30 13:34
那可否将每个键长下的scan结果输出，再进行freq计算，得到熵、吉布斯自由能的值？

可以做freq计算，但结果是没有意义的，不能用。热力学量的计算公式的前提就是当前结构的梯度为零，如果不满足这个前提，程序算出来的数没有任何物理意义

作者
Author: 594yy 时间: 2021-11-30 20:53
这个是我前两天向您请教“如何计算无能垒反应的反应速率”，您给我的文章（Baik）关于计算过程的图，图中指出计算速率常数需要计算不同键长条件下的各种热力学数据，那这些数据是不是就无法通过高斯柔性扫描获得

作者
Author: 594yy 时间: 2021-11-30 23:14

wzkchem5 发表于 2021-11-30 20:42
可以做freq计算，但结果是没有意义的，不能用。热力学量的计算公式的前提就是当前结构的梯度为零，如果不 ...

好吧。我前两天向您请教过如何计算无能垒化学反应速率，您推荐我看Baik（10.1021/acs.organomet.8b00456）的文章，我看里面需要先计算不同键长条件下的熵，电子能量等，然后才能得出无能垒过程的“吉布斯自由能垒”，最后计算出反应速率。
那么现在看来，是无法通过Gaussine来使用Baik文章的方法计算无能垒过程的反应速率，因为无法得到各个键长条件下的熵，电子能量等，对吗？

作者
Author: wzkchem5 时间: 2021-11-30 23:38

594yy 发表于 2021-11-30 16:14
好吧。我前两天向您请教过如何计算无能垒化学反应速率，您推荐我看Baik（10.1021/acs.organomet.8b00456 ...

Baik的方法是一个特例。注意在自由能垒附近，Hessian应该是有且仅有一个虚频的，而虚频的方向基本就是你扫描的那个键长。如果虚频刚好是沿着你扫描的键长的方向的话，此时可以发现freq算出的自由能刚好等于potential of mean force (PMF)，因为PMF的定义就是算Gibbs自由能的时候不考虑你的扫描坐标对自由能的贡献，算出来的数就是PMF，而虚频即使考虑进来也对自由能没有贡献，所以Baik的方法扫出的自由能曲线，在靠近平衡键长附近、没有虚频的地方是不能用、没有明确物理意义的，但是一旦扫到键开始明显解离、出现虚频的区域，freq读到的自由能可以认为是PMF，因此是有物理意义的。其中只有后者会影响计算出来的自由能垒。
所以确实我上一个帖子说得不严格，严格的说法应该是：当体系有一个虚频，且虚频和扫描的键长的方向严格一致的时候，即便不是梯度为0的结构，freq算出来的自由能也有物理意义，也就是等于PMF。如果虚频和扫描的键长的方向只是近似一致，那么freq算出来的自由能不够严格，但是这个不严格性有可能是可以接受的，至于虚频和扫描的键长的方向差多少，自由能就不可接受了，业内没有统一标准。

作者
Author: 594yy 时间: 2021-11-30 23:54

wzkchem5 发表于 2021-11-30 23:38
Baik的方法是一个特例。注意在自由能垒附近，Hessian应该是有且仅有一个虚频的，而虚频的方向基本就是你 ...

万分感谢wzk老师解答，目前我对Gaussine了解不够深，有些专业名词还是不太懂

，我能否这么理解您的回答：
也就是说，还是可以用Gaussine进行无能垒反应速率计算，过程如下：
1、寻找键解离点，采用freq保证有唯一虚频，虚频方向与键解离方向一致。
2、以上述点为起点，进行柔性扫描，获得扫描结果。
3、将扫描结果逐个进行freq计算，此时计算出来的各种热力学数据就可以用于计算“吉布斯自由能垒”了对吧？

作者
Author: wzkchem5 时间: 2021-12-1 00:40

594yy 发表于 2021-11-30 16:54
万分感谢wzk老师解答，目前我对Gaussine了解不够深，有些专业名词还是不太懂，我能否这么理解您的 ...

第一步不需要，直接做柔性扫描即可。这样扫出来的曲线，有的点没有虚频，因此算出来的自由能没有物理意义，但是这些结果本来在计算自由能垒的时候也用不着，所以画在图上也不碍事。

作者
Author: 594yy 时间: 2021-12-1 09:12

wzkchem5 发表于 2021-12-1 00:40
第一步不需要，直接做柔性扫描即可。这样扫出来的曲线，有的点没有虚频，因此算出来的自由能没有物理意义 ...

我计算出来有些点虚频不止一个可以吗？

作者
Author: wzkchem5 时间: 2021-12-1 16:20

594yy 发表于 2021-12-1 02:12
我计算出来有些点虚频不止一个可以吗？

那样不影响算出来的自由能有意义，但是影响你找到的是不是自由能垒。自由能垒的定义和过渡态类似，自由能垒处PMF必须有且仅有一个虚频，假如你的能量不止一个虚频的话，PMF大概率也不止一个虚频，因而不是自由能垒。应该检查你扫描的坐标是否合理，可能你扫描的坐标不能确切代表真正的反应坐标

作者
Author: 594yy 时间: 2021-12-1 19:36

wzkchem5 发表于 2021-12-1 16:20
那样不影响算出来的自由能有意义，但是影响你找到的是不是自由能垒。自由能垒的定义和过渡态类似，自由能 ...

好的，我尝试重新寻找路径，再次感谢wzk老师耐心的解答

作者
Author: 594yy 时间: 2021-12-2 11:05

wzkchem5 发表于 2021-12-1 16:20
那样不影响算出来的自由能有意义，但是影响你找到的是不是自由能垒。自由能垒的定义和过渡态类似，自由能 ...

wzk老师，我把自由能垒计算出来之后，可以通过sob老师发的shermo程序计算速率吗？那虚频应该设置为多少呢？

作者
Author: wzkchem5 时间: 2021-12-2 16:22

594yy 发表于 2021-12-2 04:05
wzk老师，我把自由能垒计算出来之后，可以通过sob老师发的shermo程序计算速率吗？那虚频应该设置为多少呢 ...

可以。计算的时候不要考虑隧穿，这样计算就不需要用到虚频

作者
Author: 594yy 时间: 2021-12-2 16:43

wzkchem5 发表于 2021-12-2 16:22
可以。计算的时候不要考虑隧穿，这样计算就不需要用到虚频

好的好的，谢谢老师

作者
Author: Dolantin 时间: 2023-6-2 12:13

wzkchem5 发表于 2021-12-2 16:22
可以。计算的时候不要考虑隧穿，这样计算就不需要用到虚频

请教一下W老师，在对柔性扫描的每一帧构型做freq任务然后把【虚频和扫描键长方向一致的构型的自由能】和键长拟合成曲线，取最高点作为PMF，再代入到sob老师基于TST的Excel表格中，为什么计算的时候不考虑隧穿呢？前面的freq任务不是已经得出虚频数值了吗？至少可以用Wigner法计算隧穿系数吧，如果涉及氢转移不考虑隧穿误差应该会比较大。望您不吝赐教，感谢！

作者
Author: wzkchem5 时间: 2023-6-2 13:46

Dolantin 发表于 2023-6-2 05:13
请教一下W老师，在对柔性扫描的每一帧构型做freq任务然后把【虚频和扫描键长方向一致的构型的自由能】和 ...

原因有以下几个：
（1）Wigner等方法推导时用到了“过渡态是能量势能面上的鞍点”的条件，当你用的结构是PMF上的最高点时，这个前提条件不成立，很难讲对于隧穿系数的影响多大。
（2）假如PMF最高点和过渡态非常接近，那么Wigner等方法仍然可用，但是此时你也没必要费劲算那么多自由能、然后取PMF最高点。
（3）你回复的我的那个帖子讲的是用Baik法计算无垒反应速率，无垒反应的自由能垒一般在解离区，因此虚频小，隧穿校正小，即使忽略也没什么影响，估计小于Baik法本身的误差。如果算有垒反应，可以用Polyrate做基于IRC的VTST计算，比Baik法这类基于柔性扫描的方法严格，且此时可以用SCT、LCT等方法计算隧穿。

作者
Author: Dolantin 时间: 2023-6-2 14:49

wzkchem5 发表于 2023-6-2 13:46
原因有以下几个：
（1）Wigner等方法推导时用到了“过渡态是能量势能面上的鞍点”的条件，当你用的结构 ...

感谢w老师的细致解答。无垒反应的自由能垒一般在解离区，是因为处于解离区状态的分子比较稳定，所以虚频小吗？这一点，还不是很明白。

作者
Author: wzkchem5 时间: 2023-6-3 01:28

Dolantin 发表于 2023-6-2 07:49
感谢w老师的细致解答。无垒反应的自由能垒一般在解离区，是因为处于解离区状态的分子比较稳定，所以虚频 ...

不是。你把解离曲线画出来，由图易知随着键长的增加，虚频逐渐趋于0。所以一般而言，同等条件下，键长越长的过渡态，虚频的绝对值越小。
这个和稳定性没关系，平衡位置附近的分子，稳定性比解离区更强，但是离平衡键长近的过渡态，虚频反而比接近解离的过渡态要大得多。

作者
Author: Dolantin 时间: 2023-6-9 09:54

wzkchem5 发表于 2023-6-3 01:28
不是。你把解离曲线画出来，由图易知随着键长的增加，虚频逐渐趋于0。所以一般而言，同等条件下，键长越 ...

学习了，谢谢W老师

作者
Author: didiao 时间: 2025-6-2 19:32

wzkchem5 发表于 2021-11-30 23:38
Baik的方法是一个特例。注意在自由能垒附近，Hessian应该是有且仅有一个虚频的，而虚频的方向基本就是你 ...

老师您好，我想计算一个两个原子结合成产物的自由能，看了baik的文章，我有几个问题
1.我看他的-TS计算的前几个点是单调减小的，他去拟合了sigmoidal function，我柔性扫描间隔0.1埃，我计算了前9个点，前5个从-16.65kcal/mol降到了-17.46，第6、7、8点出现了反应方向的虚频，-TS是从-16.57— -16.61— -16.64，然后第九个点有没有虚频了，-TS是-18.129，我这个就不是单调减小的了，怎么去拟合sigmoidal function呢？
2.我看您说只有虚频方向对的，算出来的自由能才是有用的，（1）那是不是我在5-9点之间再减小0.05埃步长扫描一下，算出自由能最高的点，（2）再单独拿出两个原子分解计算能量再加和，（3）在把优化好的产物计算出能量，这三个点的能量就可以算出正向逆向的自由能能垒了？

作者
Author: wzkchem5 时间: 2025-6-2 21:02

didiao 发表于 2025-6-2 19:32
老师您好，我想计算一个两个原子结合成产物的自由能，看了baik的文章，我有几个问题
1.我看他的-TS计算 ...

也可以不拟合，简单地取自由能最大的点。我的文章就这么做过，可以在https://pubs.acs.org/doi/10.1021/jacs.4c18445的SI里搜Baik

作者
Author: didiao 时间: 2025-6-2 22:21

wzkchem5 发表于 2025-6-2 21:02
也可以不拟合，简单地取自由能最大的点。我的文章就这么做过，可以在https://pubs.acs.org/doi/10.1021/j ...

非常感谢您的解答

作者
Author: didiao 时间: 2025-6-4 17:02

wzkchem5 发表于 2025-6-2 21:02
也可以不拟合，简单地取自由能最大的点。我的文章就这么做过，可以在https://pubs.acs.org/doi/10.1021/j ...

老师好，我看您的计算中，随着键的拉长，可以找到最大的自由能点，这个点减去两个产物的自由能和是10.5，然后您又加了个反应的自由能，得到自由能垒最终估计值，这个反应自由能▲G是怎么算出来的呢？
对于我的两个物质合成一个物质的反应来说，是不是随着键拉长的自由能最高的点减去两个物质的自由能和，就是反应势垒呢？

作者
Author: wzkchem5 时间: 2025-6-4 21:00

didiao 发表于 2025-6-4 17:02
老师好，我看您的计算中，随着键的拉长，可以找到最大的自由能点，这个点减去两个产物的自由能和是10.5， ...

反应自由能的计算方法已经在SI其他章节说清楚了。就是高精度单点+低精度的Gibbs自由能校正量
是的

作者
Author: didiao 时间: 2025-6-5 15:45

wzkchem5 发表于 2025-6-4 21:00
反应自由能的计算方法已经在SI其他章节说清楚了。就是高精度单点+低精度的Gibbs自由能校正量
是的

谢谢老师

作者
Author: didiao 时间: 2025-6-30 16:36

wzkchem5 发表于 2025-6-4 21:00
反应自由能的计算方法已经在SI其他章节说清楚了。就是高精度单点+低精度的Gibbs自由能校正量
是的

老师，对于我计算KOH=K+OH，我有几个问题
1.它进性柔性扫描，得到的是K·+OH·这样的双自由基还是K+和OH-呢？
2.如果得到的是K+和OH-，那我就只能计算K+ +OH-=KOH的能垒吗？这时候减去的两个反应物能量就是K+和OH-吗？
3.如果是2的情况，有没有办法通过baik的方法得到两个自由基反应的能垒呢？

作者
Author: wzkchem5 时间: 2025-6-30 17:08

didiao 发表于 2025-6-30 16:36
老师，对于我计算KOH=K+OH，我有几个问题
1.它进性柔性扫描，得到的是K·+OH·这样的双自由基还是K+和OH ...

唯一确定的方法是实际扫出来看看（注意做broken symmetry，并留意是否收敛到了稳定波函数，broken symmetry不能保证波函数稳定）。KOH在水溶液里是异裂的，但气相相比水溶液强烈不利于异裂，不好说哪个因素更重要。
不管产物是离子还是自由基，都可以用Baik的方法算，但是考虑到你的体系很小，完全可以用VRC-VTST乃至AIMD等更精确的方法算，如果用Baik的方法算反而会被质疑，说更精确的方法你又不是用不起，为什么不用

作者
Author: didiao 时间: 2025-6-30 17:14

wzkchem5 发表于 2025-6-30 17:08
唯一确定的方法是实际扫出来看看（注意做broken symmetry，并留意是否收敛到了稳定波函数，broken symmet ...

老师，VRC-VTST应该用啥软件呢？

作者
Author: wzkchem5 时间: 2025-6-30 19:55

didiao 发表于 2025-6-30 17:14
老师，VRC-VTST应该用啥软件呢？

用polyrate

作者
Author: Uus/pMeC6H4-/キ 时间: 2025-9-25 14:52

wzkchem5 发表于 2021-11-30 23:38
Baik的方法是一个特例。注意在自由能垒附近，Hessian应该是有且仅有一个虚频的，而虚频的方向基本就是你 ...

王老师，关于这个Baik法（Pitfalls in Computational Modeling of Chemical Reactions and How To Avoid Them, Organometallics 2018, 37, 19, 3228–3239, DOI: 10.1021/acs.organomet.8b00456）处理找不到过渡态/无垒（无电子能垒）反应的话题，我读完此文理解的过程是：
1. 对优化好的反应物与产物做振动分析，获得熵对自由能的贡献–TS；
2. 对反应物与产物做单点计算，获得电子能E(SCF)；
3. 沿反应坐标做柔性扫描，获得与反应物或产物相似的几个点的结构与电子能E(SCF)；
4. 对这几个点做振动分析，获得熵对自由能的贡献–TS；
5. 绘制E(SCF)对反应坐标的曲线并用Morse势函数拟合；
6. 绘制–TS对反应坐标的曲线并用sigmoid函数拟合；
7. 将两个拟合结果合并，得到存在极大值的自由能曲线，用于算势垒。

单看文中的描述还是不太明白的几处细节，不知可否解惑：
1. 文中给的例子的反应坐标只需一个分子内坐标的距离r来表示，如果反应坐标涉及多个内坐标的变化，需要像增强采样动力学那样找一个集合变量来表示并扫描该广义内坐标吗？
2. 都说柔性扫描涉及多个变量时会有扫描顺序与结构弛豫的滞后性问题，此处是否需要关心呢？
3. 可否用完全由原子受力决定结构变化的IRC(downhill)任务代替内坐标的柔性扫描，并以质权坐标作为拟合所需的反应坐标？
4. 文中只从反应物一侧出发获得相似的几个点的结构，能否同时从反应物与产物两侧出发获得更多结构并改善拟合效果？
5. 拟合–TS对反应坐标的曲线时用的sigmoid函数是mathworld提到的logistic function 1/(1+exp(–x))吗，还是wikipedia提到的如误差函数、tanh x、arctan x、x/sqrt(1+x^2)等也能用，理论依据是什么？
6. 合并曲线时为什么直接用了电子能，而不加上（已经在振动分析中一起获得了的）零点能以及焓对自由能的贡献部分？
7. 该方法能找到能垒似乎是以“两条拟合曲线趋势相反且能量的数量级相当”为前提的，如何有信心保证这两点呢？

作者
Author: wzkchem5 时间: 2025-9-25 15:01
本帖最后由 wzkchem5 于 2025-9-25 15:04 编辑

Uus/pMeC6H4-/キ发表于 2025-9-25 14:52
王老师，关于这个Baik法（Pitfalls in Computational Modeling of Chemical Reactions and How To Avoid ...

1. 一般都能找到某两个原子/某两个片段质心的距离作为这个变量。如果反应物可以发生一个无垒反应，通过能量单调下降的途径得到产物，而反应物又没有虚频的话，容易证明反应物必然不止一个分子。而双分子反应起码在反应初期，必然是可以用两个分子的质心距离（或者某一对原子之间的距离）作为反应坐标的
2. 肯定有一些影响，个人觉得这个因素的影响小于Baik法本身的误差，但目前未见文献支持
3. downhill的问题在于：（1）初始结构必须是两个分子离得比较远，但不是无限远的结构；（2）downhill路径依赖于初始结构的选取。如何选择IRC的初始结构而不引入人为因素，是一个问题
4. 有可能，但未必。当两个分子离得足够远的时候，谐振近似会失效，使得拟合结果无意义。键长范围不是越宽越好
5. 我觉得没有理论依据。我都是不拟合，直接取自由能最大值点，也没见人质疑（https://doi.org/10.1038/s44160-025-00821-8；https://doi.org/10.1021/jacs.4c18445）
6. 7. 见5

作者
Author: Uus/pMeC6H4-/キ 时间: 2025-9-25 15:58

wzkchem5 发表于 2025-9-25 15:01
1. 一般都能找到某两个原子/某两个片段质心的距离作为这个变量。如果反应物可以发生一个无垒反应，通过能 ...

感谢解答！小声说一句，要是当初要是单独撰文详解该方法就好了，而不是现在这样藏在一篇杂谈/教程的一节里。

3. 确实，反应物是原子受力接近0的极小点的话跑不起来IRC，沿反应坐标电子能单调上升的情况也与downhill方向相反[而IRC(uphill)闻所未闻]；要产生距离较远的loose分子团簇初始结构还得老实做柔性扫描。另一边，NEB类对反应路径插点的算法也用不了。
4. 我能理解分子距离远时谐振近似无法描述好势能面局域形状这点，当时写这问题主要是在设想反应物与产物都是势能面驻点、因而柔性扫描两端都有谐振近似仍生效的点的情况。
5-7. “直接取自由能最大值点”是在做了振动分析的反应物、过渡态（若有）、产物及所有相近的点中选的么，那确实省事。不过我还想确认一下，无论是不是驻点、有几个振动虚频，一般都能认为振动分析所得简正坐标已经相互正交，因而不用额外做“把振动虚频矢量在其他实频矢量上的投影去掉”这种操作，对吗？

作者
Author: wzkchem5 时间: 2025-9-25 20:35

Uus/pMeC6H4-/キ发表于 2025-9-25 15:58
感谢解答！小声说一句，要是当初要是单独撰文详解该方法就好了，而不是现在这样藏在一篇杂谈/教程的一节 ...

对，柔性扫描曲线上所有点分别算自由能，取最大的。
振动分析的简正坐标by definition是严格正交的，这和虚频数目、是不是驻点无关，而是由矩阵对角化的基本性质决定的

作者
Author: Uus/pMeC6H4-/キ 时间: 2025-10-15 10:52
本帖最后由 Uus/pMeC6H4-/キ于 2025-10-20 17:27 编辑

wzkchem5 发表于 2025-9-25 20:35
对，柔性扫描曲线上所有点分别算自由能，取最大的。
振动分析的简正坐标by definition是严格正交的，这 ...

看Gaussian的freq关键词的projected选项才回忆起来，当时我想说的应该不是简正坐标之间的投影，而是这个由梯度求得IRC路径切线后垂直于路径方向的投影……

Projected
For a point on a mass-weighted reaction path (IRC), compute the projected frequencies for vibrations perpendicular to the path. For the projection, the gradient is used to compute the tangent to the path. Note that this computation is very sensitive to the accuracy of the structure and the path [Baboul97]. Accordingly, the geometry should be specified to at least 5 significant digits. This computation is not meaningful at a minimum.

就像普通的振动分析所得简正坐标可以分解为各种内坐标的贡献一样，IRC路径的反应坐标应该也能分解为各种简正坐标或各种内坐标的贡献。我试图问的是，对本帖讨论的柔性扫描求自由能曲线而言，虽然不是IRC但计算的各点同样有明显的非零梯度，那是否也需要做这种投影后的振动分析获得热力学校正量？

[Baboul97] J. Chem. Phys. 107, 9413–9417 (1997) DOI: https://doi.org/10.1063/1.475238

编辑：这个老帖子也相关，但是不知道算不算回答了我的问题

作者
Author: wzkchem5 时间: 2025-10-15 12:52

Uus/pMeC6H4-/キ发表于 2025-10-15 10:52
看Gaussian的freq关键词的projected选项才回忆起来，当时我想说的应该不是简正坐标之间的投影，而是这个 ...

如果梯度与虚频方向相同，那么投影不投影都是一样的，且都是可用的，因为虚频本来就对热力学校正量没有贡献。
柔性扫描曲线在能垒附近，往往近似满足这个条件，但一般不严格满足。这时做投影或许更严格一些，但是很难说这样带来的改善是否足够大。我个人怀疑改善不大，因为在解离区有很多额外的问题，例如非谐性，这些因素导致的误差可能更大

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