关于cp2k的metadynamics的后处理的问题

DwyaneWan

各位老师好，学生我最近用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？ 最后感谢各位老师的回复！！

丁越

你理解错了我说的话，你看自由能的计算公式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
你的配位数出现了负值大概率是参数设置的有问题，你看下面的图片就因该清楚了：

chands

丁越 发表于 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 可以不做区分，但是我想不清楚。希望得到解答。

丁越

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

DwyaneWan

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

感谢老师详细回答

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

丁越

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

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

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

Daniel_Arndt

在我看来，如果想要把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，里面用到了一些条件概率的知识。

DwyaneWan

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

非常感谢您的回复！

chands

丁越 发表于 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

Daniel_Arndt

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

DwyaneWan

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

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

JCenter

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

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

丁越

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

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

JCenter

丁越 发表于 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的配位数负值了。

丁越

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

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