计算化学公社

标题: 关于cp2k的metadynamics的后处理的问题 [打印本页]

作者
Author: DwyaneWan 时间: 2022-7-8 05:20
标题: 关于cp2k的metadynamics的后处理的问题
各位老师好，学生我最近用CP2K的metadynamics功能，CV选用的是coordination number。然后r0选取的是2.5 A，指数分别是6和12。然后我跑了30 ps，用graph.psm功能算到的自由能面的图如图1.

我通过计算发现coordination number不可能为负值，然后我搜索了论坛里的相关回答：http://bbs.keinsci.com/thread-30067-1-1.html。丁越老师说这个CV通过graph计算出来的是概率密度。我想问的是如果我想通过这个CV如何才能换算成我设置的coordination number的范围的CV呢？

我在阅读别的文献，发现他通过coordination number 做CV得到的是distance的自由能面。如图二。我想请问，如果用CN设定为metadynamics的CV，如何将结果的CV转换成distance呢？是不是需要在跑模拟的过程中再设定两个distance的colvar？最后感谢各位老师的回复！！

作者
Author: 丁越 时间: 2022-7-8 09:26
本帖最后由丁越于 2022-8-21 22:49 编辑

你理解错了我说的话，你看自由能的计算公式F=-k_B*T*ln(P)，其中P是CVs的概率密度。求概率最简单的方法当然是把CVs的数值范围内划分为n个bin间隔，然后统计CVs落到每个bin间隔内的概率密度，然后就得到了每个bin的自由能，把这些bin间隔得到的自由能连起来就是一维自由能曲线了。
对于第二个问题你想通过bias配位数CVs的结果去计算以距离为CVs的自由能，这并不是简单的将已经添加有偏势的结果下直接去计算距离为CVs的自由能面，而是先需要reweighting将有偏的结果转化为无偏，然后才可以计算以距离为CVs的自由能。你可以参考这篇文献：Reconstructing the Equilibrium Boltzmann Distribution
from Well-Tempered Metadynamics
你的配位数出现了负值大概率是参数设置的有问题，你看下面的图片就因该清楚了：

作者
Author: chands 时间: 2022-7-8 09:38

丁越发表于 2022-7-8 09:26
你理解错了我说的话，你看自由能的计算公式F=-k_B*T*ln(P)，其中P是CVs的概率密度。求概率最简单的方法当然 ...

我这里借楼问一个问题，如果在NVT系综中做溶液的metadynamics，并计算自由能，根据定义，算的是Helmholtz自由能，包括你这里写得公式也是Helmholtz自由能，但是如果还用量子化学方法计算，得到的是Gibbs自由能。F = G - PV，差了一项 PV。我隐约觉得液相的delta (PV) 项是可以忽略的，这时反应的delta F 和delta G 可以不做区分，但是我想不清楚。希望得到解答。

作者
Author: 丁越 时间: 2022-7-8 09:55
我也是正在学习这些，按照我的理解，F并不是特指亥姆霍兹自由能，这就要看你在哪个系综下采样了，如果是在NPT系综下采样，那么PV的变化就体现其中，得到的F其实是吉布斯自由能。不知道我的理解是否正确？

作者
Author: DwyaneWan 时间: 2022-7-9 01:56

丁越发表于 2022-7-8 09:26
你理解错了我说的话，你看自由能的计算公式F=-k_B*T*ln(P)，其中P是CVs的概率密度。求概率最简单的方法当然 ...

感谢老师详细回答

那如果我也将这些输出的CV转化成我定义的CV，是不是也要reweighting呢？

作者
Author: 丁越 时间: 2022-7-9 09:03

DwyaneWan 发表于 2022-7-9 01:56
感谢老师详细回答

那如果我也将这些输出的CV转化成我定义的CV，是不是也要reweighting呢？

你的有偏势能是添加在了配位数CVs中，而距离并没有bias，所以要想通过bias配位数的结果计算以距离为CVs的自由能的话就需要reweighting。你可以看看那篇论文了解下。

作者
Author: Daniel_Arndt 时间: 2022-7-10 08:12
在我看来，如果想要把coordination number作为collective variable的free energy profile转换成distance作为collective variable的free energy profile的话，首先要确认你已经跑出的轨线同时充分采样了coordination number和distance，然后这个转换才可以说得上是有意义的。在此前提下，如果你跑轨线时用的coordination number和你绘图时用的distance都适合描述你关心的过程的话，你才能在你的论文里做出有指导意义的讨论。

具体该如何操作，你可以看这篇文章 https://doi.org/10.1016/j.jcat.2020.04.015 的Supporting Information里面的7.2 Transformation of 1D free energy profiles，里面用到了一些条件概率的知识。

作者
Author: DwyaneWan 时间: 2022-7-11 05:52

Daniel_Arndt 发表于 2022-7-10 08:12
在我看来，如果想要把coordination number作为collective variable的free energy profile转换成distance作 ...

非常感谢您的回复！

作者
Author: chands 时间: 2022-7-11 08:12
本帖最后由 chands 于 2022-7-11 14:42 编辑

丁越发表于 2022-7-8 09:55
我也是正在学习这些，按照我的理解，F并不是特指亥姆霍兹自由能，这就要看你在哪个系综下采样了，如果是在N ...

谢谢。前两天有点忙没回复。查了统计力学的书，你是对的，公式F= - k_B*T*ln(Q) (我这里用配分函数Q代替概率密度），对于NVT系综，F是Helmholtz自由能A，Q用NVT系综的配分函数，对于NPT系综，F是Gibbs自由能G，Q用NPT系综的配分函数。不管NVT还是NPT，概率密度与Q相关。
在NVT系综里，delta A = - R*T*ln (Kv)， Kv是等温等容条件下的平衡常数。在NPT系综里， delta G = - R*T*ln (K）, K是等温等压条件下的平衡常数。

今天逛乎见到如下一段话（注意这里F是Helmholtz自由能）。我想NPT系综中，溶液反应前后V变化很小，delta (PV)是可以忽略，这时Helmholtz自由能变与Gibbs自由能变可以不作区分。但是气相反应的delta (PV)可能变化很大，Helmholtz自由能变与Gibbs自由能变就不是一回事了。我再查一下资料。https://zhuanlan.zhihu.com/p/28013840

作者
Author: Daniel_Arndt 时间: 2022-8-17 13:26
楼主贴的第二张图，不会是从这里 http://bbs.keinsci.com/thread-23778-1-1.html 讨论的文献中拿来的吧？那篇文章看的时候，得小心点。

作者
Author: DwyaneWan 时间: 2022-8-20 11:31

Daniel_Arndt 发表于 2022-8-17 13:26
楼主贴的第二张图，不会是从这里 http://bbs.keinsci.com/thread-23778-1-1.html 讨论的文献中拿来的吧？那 ...

您好是的。请问应该要注意点什么呢？还望不吝赐教，感谢同学

作者
Author: JCenter 时间: 2025-6-23 17:00
本帖最后由 JCenter 于 2025-6-24 01:17 编辑

丁越发表于 2022-7-8 09:26
你理解错了我说的话，你看自由能的计算公式F=-k_B*T*ln(P)，其中P是CVs的概率密度。求概率最简单的方法当然 ...

丁越老师，您好。我看了好多帖子，还是没有理解，配位数CV值出现负值的情况（，理解不了这个负号的含义。这个负号，您能给解疑答惑一下不（我的配位数设置：NN=8,ND=14,R0=1.6埃，使用了wall,低于0.015值，发挥四次势墙作用）

作者
Author: 丁越 时间: 2025-6-24 08:15
本帖最后由丁越于 2025-6-24 09:41 编辑

JCenter 发表于 2025-6-23 17:00
丁越老师，您好。我看了好多帖子，还是没有理解，配位数CV值出现负值的情况（，理解不了这个负号的含义。 ...

我不太清除你的输入文件怎么设置的，按照配位数公式不会出现负值的情况。配位数作为CV时候r0参数的取值比较重要, 一般对于化学键来说通常会选择过渡态的键长作为r0，或者选择一倍多的键长作为r0足矣。对于其他情况需要跑一段时间MD，然后根据CV的变化来设置r0参数。

作者
Author: JCenter 时间: 2025-6-24 17:00

丁越发表于 2025-6-24 08:15
我不太清除你的输入文件怎么设置的，按照配位数公式不会出现负值的情况。配位数作为CV时候r0参数的取值比 ...

我设置配位数就是按照这种思路来的，colvar.metaldynlog文件中的配位数C值也确实都是大于0的值，就是在计算自由能，得到fes.dat文件时（graph.psmp -ndim 2 -ndw 1 2 -file xxx-1.restart -cp2k）后，配位数就出来从-0.8至-0.1的配位数负值了。

作者
Author: 丁越 时间: 2025-6-24 20:14

JCenter 发表于 2025-6-24 17:00
我设置配位数就是按照这种思路来的，colvar.metaldynlog文件中的配位数C值也确实都是大于0的值，就是在计 ...

我不用cp2k自带的MTD模块，更推荐使用plumed, 输入文件写起来非常方便，比起cp2k的繁琐输入文件不知道都简化了多少倍，而且plumed的后处理模块也做得很好，同时也支持更加丰富的MTD方法变体。

作者
Author: JCenter 时间: 2025-6-24 21:58

丁越发表于 2025-6-24 20:14
我不用cp2k自带的MTD模块，更推荐使用plumed, 输入文件写起来非常方便，比起cp2k的繁琐输入文件不知道都 ...

好嘞，感谢大佬。后面我再学习一下plumed插件。最后，若使用well-tempered metadynamics 设置偏置因子后是不是要做不同的偏置因子下的样本采集效率对比呢

作者
Author: 丁越 时间: 2025-6-25 22:00

JCenter 发表于 2025-6-24 21:58
好嘞，感谢大佬。后面我再学习一下plumed插件。最后，若使用well-tempered metadynamics 设置偏置因子后 ...

不用，花那时间不如多跑跑MTD了。

偏置因子gamma=(T + deltaT) / T, 高斯峰被scale的程度正比于exp(-V/K_B * deltaT), 也就是说当gamma越大时即deltaT越大，这个指数项也就越接近1，初期高斯峰很高就能迅速的促进相空间探索。反之当gamma设置过小时，高斯峰被scale的越小，更早的进入小振幅填充，虽然自由能面填充的比较精细了，但是就需要更长的模拟时间。同时也存在gamma过小时CV跨越势垒的能力就变得很差。所以gamma的选择与CV的势垒deltaG密切相关，经验原则是选取deltaG/K_B*T。

另外也要注意hill的初始高度的设置，一般原则也是设置几分之一的K_B*T, 最好不要超过体系的势垒。

作者
Author: JCenter 时间: 2025-6-26 00:37
本帖最后由 JCenter 于 2025-6-26 04:06 编辑

丁越发表于 2025-6-25 22:00
不用，花那时间不如多跑跑MTD了。

偏置因子gamma=(T + deltaT) / T, 高斯峰被scale的程度正比于exp(-V ...

我也是这么理解的。但后面这句话是不是打错了呀

：“gamma的选择与CV的势垒deltaG密切相关，经验原则是选取deltaG/K_B*T”。应该是kB*(T+daltaT)与deltaG相当，根据这个关系，就能得出deltaT值，从而确定偏置因子gamma值吧。即玻尔兹曼常数乘偏置后温度(T+deltaT)，而不是乘T。
此外，如果使用well-tempered方法，初始的峰高hill不需要按照kBT的几分之一吧。除非是使用标准mtd,每次加峰的高度都是一定的，这种情况下峰高越小更精确也更耗时。(如果是使用well-tempered方法，我不确定这里的kBT中的T是指T,还是指T+deltaT，请您指点一下)。随着步数的增加，添加的高斯峰都在scale，deltaT越高，添加的高斯峰峰高越来越小。这么看，初始峰高WW只要不超过能垒deltaG，一个低的初始峰高都可以吧。还请大佬谈下您的理解。

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