计算化学公社

标题: CP2K软件metadynamics模块的自由能后处理问题 [打印本页]

作者
Author: xingyuan 时间: 2018-8-5 21:02
标题: CP2K软件metadynamics模块的自由能后处理问题
本帖最后由 xingyuan 于 2018-8-5 21:11 编辑

各位老师大家好。最近在学习CP2K软件上面的Metadynamics模块，看到一篇文献，我想实现的是类似的功能，即想得到的是如文献上面的图3，为了方便给大家截图下来了。
所以我想问MTD模拟后如何得到free-energy的变化情况，如何后处理（我没有找到关于free energy的数据在哪里）？望各位老师教我。感激不尽！！！
以下是我查到的资料，不知道对不，望各位老师指正。

我在CP2K官网上看到一个例子，是这样做的，关于设置MTD的关键词如下：
&COLVAR
   &COORDINATION
      KINDS_FROM  O
      KINDS_TO C
      R_0 [angstrom]  1.8
      NN  8
      ND  14
   &END COORDINATION
&END COLVAR
.......

&FREE_ENERGY
&METADYN
   DO_HILLS
   NT_HILLS 100
   WW 3.0e-3
   &METAVAR
      SCALE 0.2
      COLVAR 1
   &END METAVAR
结果会得到一个-colvar.metadynLog，一个-hills.metadynLog文件，分别表示CV值的变化和HILLS的坐标等。但是没有提到free-energy的值，希望老师们教导。感谢。

作者
Author: youyno 时间: 2018-8-6 09:48
graph.popt是cp2k编译时产生的用于metadynamics后处理自由能的程序，后缀根据版本不同有所区别，具体命令可以去到安装包内cp2k-4.1/src/metadyn_tools中的文件查看，graph.F文件中包含了详细的参数解释。基本的命令如下
graph.popt -ndim 需要计算自由能的CV个数 -ndw 那几个cv(如 1 2或1 3) -file 实时更新的restart文件（通常是作业名-1.restart） -cp2k 执行完毕后会生成fes.dat的文件，最后一列为自由能单位是a.u.前面几列为CV

作者
Author: xingyuan 时间: 2018-8-6 17:58
本帖最后由 xingyuan 于 2018-8-18 21:09 编辑

youyno 发表于 2018-8-6 09:48
graph.popt是cp2k编译时产生的用于metadynamics后处理自由能的程序，后缀根据版本不同有所区别，具体命令可 ...

非常感谢您的解答！！
按照您的说的我遇到了点小麻烦，MTD模拟后产生了COLVAR.metadynLog文件和HILLS.metadynLog文件，但是用graph.popt命令需要有HILLS文件，我尝试把这几个文件重命名为HILLS仍然没有成功。目前结果是把.restart文件重命名为HILLS，得到的fes.dat是空文件
下面几张图分别是：MTD模拟后该文件夹下面所含的文件（有一部分删掉了），另外两张图分别是按照您说的提交命令文件以及生成的.out文件。您能帮忙看一下应该如何做吗？
[attach]14948[/attach]

作者
Author: youyno 时间: 2018-8-6 20:05
老哥，restart文件前要加-file命令，这么简单的程序没有必要提交到服务器上算，主节点就可以算也不需要并行

作者
Author: xingyuan 时间: 2018-8-6 20:22
本帖最后由 xingyuan 于 2018-8-6 20:28 编辑

youyno 发表于 2018-8-6 20:05
老哥，restart文件前要加-file命令，这么简单的程序没有必要提交到服务器上算，主节点就可以算也不需要并行 ...

好的好的，感谢.可是加了-file命令，out文件中还是有 forrtl: severe (174): SIGSEGV, segmentation fault occurred错误，然后fes.dat是空文件。这个是什么错误啊

作者
Author: youyno 时间: 2018-8-7 16:52
你没有理解我的意思
命令应该是 graph.popt -ndim 2 -ndw 1 2 3 -file gr2hno3_mtd_3cv_p1-1_2000.restart -cp2k

作者
Author: xingyuan 时间: 2018-8-10 15:12

youyno 发表于 2018-8-7 16:52
你没有理解我的意思
命令应该是 graph.popt -ndim 2 -ndw 1 2 3 -file gr2hno3_mtd_3cv_p1-1_2000.restart ...

非常感谢！！

作者
Author: xingyuan 时间: 2018-8-18 21:09
本帖最后由 xingyuan 于 2018-8-18 21:10 编辑

youyno 发表于 2018-8-7 16:52
你没有理解我的意思
命令应该是 graph.popt -ndim 2 -ndw 1 2 3 -file gr2hno3_mtd_3cv_p1-1_2000.restart ...

@youyno 老师，您好，又来麻烦您了。还是关于CP2K的问题，我有以下几个疑问：
1、CP2K Manual上关于计算自由能有三种方法：Alchemical Change 、Metadynamics、Umbrella Integration。其中另外两种有类似于MTD后处理的方法吗？或者应该怎样去处理数据？我利用Umbrella的方法计算了一下，并没有输出有关自由能的文件。输入文件附上。
2、文献上面提到进行Potential of Mean Force的计算用的是Blue Moon的方法，这个在CP2K中是否支持？还是说CP2K自带的限制计算就是用的Blue moon的原理呢？我的做法关键词是以下形式，这是否正确？非常感谢！
   &COLVAR
   &DISTANCE
      ATOMS 90 61
   &END DISTANCE

&MOTION
&CONSTRAINT
&COLLECTIVE
   COLVAR 1
   INTERMOLECULAR
   TARGET [angstrom] 8.0
&END COLLECTIVE
&LAGRANGE_MULTIPLIERS
   COMMON_ITERATION_LEVELS 1
&END
&END CONSTRAINT

麻烦您了。

作者
Author: youyno 时间: 2018-8-21 21:32

xingyuan 发表于 2018-8-18 21:09
@youyno 老师，您好，又来麻烦您了。还是关于CP2K的问题，我有以下几个疑问：
1、CP2K Manual上关于计算 ...

你好，原理方面我不是很清楚，你可以去cp2k的google group中去提问，软件开发者会给予答案，另外我用的限制方法一般是QUADRATIC

作者
Author: xingyuan 时间: 2018-8-22 09:19

youyno 发表于 2018-8-21 21:32
你好，原理方面我不是很清楚，你可以去cp2k的google group中去提问，软件开发者会给予答案，另外我用的限 ...

好的，感谢

作者
Author: zyj19831206 时间: 2018-9-3 15:38

youyno 发表于 2018-8-21 21:32
你好，原理方面我不是很清楚，你可以去cp2k的google group中去提问，软件开发者会给予答案，另外我用的限 ...

cp2k貌似学习资源不多。。

作者
Author: yjr 时间: 2018-9-3 15:57

xingyuan 发表于 2018-8-18 21:09
@youyno 老师，您好，又来麻烦您了。还是关于CP2K的问题，我有以下几个疑问：
1、CP2K Manual上关于计算 ...

你这个也叫thermodynamics intergration方法，自由能是grep Shake *.LagrangeMultLog | awk '{c++ ; s=s+$4}END{print s/c}'。每一步改变TARGET的值得到对应的自由能，就可以画自由能曲线了。

作者
Author: xingyuan 时间: 2018-9-13 14:28

yjr 发表于 2018-9-3 15:57
你这个也叫thermodynamics intergration方法，自由能是grep Shake *.LagrangeMultLog | awk '{c++ ; s=s+ ...

感谢感谢。manual上面有一个例子是这样做的，一直不知道叫什么方法，学到了。那您知道计算自由能的umbrella integration的后处理方法吗？感谢

作者
Author: xingyuan 时间: 2018-9-13 14:29

zyj19831206 发表于 2018-9-3 15:38
cp2k貌似学习资源不多。。

嗯嗯，所以一个问题就能卡住好些天。

作者
Author: liujodan 时间: 2019-10-16 09:15

youyno 发表于 2018-8-6 09:48
graph.popt是cp2k编译时产生的用于metadynamics后处理自由能的程序，后缀根据版本不同有所区别，具体命令可 ...

youyno老师您好，请问向我这个文件夹中，没有找到graph.popt,是不是没有编译上呢？还请指点，谢谢

作者
Author: youyno 时间: 2019-10-16 11:18

liujodan 发表于 2019-10-16 09:15
youyno老师您好，请问向我这个文件夹中，没有找到graph.popt,是不是没有编译上呢？还请指点，谢谢

graph.popt应该在和cp2k.popt所在的同一个文件夹内

作者
Author: liujodan 时间: 2019-10-16 14:57

youyno 发表于 2018-8-6 09:48
graph.popt是cp2k编译时产生的用于metadynamics后处理自由能的程序，后缀根据版本不同有所区别，具体命令可 ...

youyno老师您好,graph.pop 文件问题已经解决，请问自由能单位a.u.是 Ha吗？

作者
Author: youyno 时间: 2019-10-16 20:19

liujodan 发表于 2019-10-16 14:57
youyno老师您好,graph.pop 文件问题已经解决，请问自由能单位a.u.是 Ha吗？

是的！

作者
Author: liujodan 时间: 2019-10-17 08:35

youyno 发表于 2019-10-16 20:19
是的！

感谢回复，另外想请教一下，使用restart输出文件中的 -1.rastart 文件和-2000.restart （2000步输出）文件有什么区别吗? 从结果上看好像比较相似。

作者
Author: highlight 时间: 2019-10-17 08:59

liujodan 发表于 2019-10-17 08:35
感谢回复，另外想请教一下，使用restart输出文件中的 -1.rastart 文件和-2000.restart （2000步输出）文 ...

当然不一样，你到第 2000 步和到最后一步，加过的 hills 数量肯定不一样吧，高斯展开出来的结果能完全一样才奇了怪了。

作者
Author: liujodan 时间: 2019-10-19 20:49

highlight 发表于 2019-10-17 08:59
当然不一样，你到第 2000 步和到最后一步，加过的 hills 数量肯定不一样吧，高斯展开出来的结果能完全一 ...

感谢

，初学METAD 请多指教

作者
Author: liujodan 时间: 2019-10-19 20:54

youyno 发表于 2019-10-16 20:19
是的！

谢谢 youyno 老师，另外，还有一个问题想请教，我在做METAD时，设置了6 个 COLVAR，跑完后，想通过graph.popt进行FES分析，graph。popt最多能输入 CV个数为2，当CV大于2时就出错，还请请教是什么原因

作者
Author: youyno 时间: 2019-10-28 10:01

liujodan 发表于 2019-10-19 20:54
谢谢 youyno 老师，另外，还有一个问题想请教，我在做METAD时，设置了6 个 COLVAR，跑完后，想通过graph ...

你打开COLVAR.log和HILLs.Log文件看看是不是有6个CV

作者
Author: 413 时间: 2019-10-30 14:42

youyno 发表于 2019-10-16 11:18
graph.popt应该在和cp2k.popt所在的同一个文件夹内

您好，请教一个问题哈。
我现在在跑一个分子的metadynamics，variable是两个键长，我设置了数值为正的scale之后发现两个键长均往较大的方向变化。请问scale的正负是否代表了键长变化的方向？或者怎么处理才能使键长只往减小的方向变化？

作者
Author: youyno 时间: 2019-10-30 15:22

413 发表于 2019-10-30 14:42
您好，请教一个问题哈。
我现在在跑一个分子的metadynamics，variable是两个键长，我设置了数值为正的sc ...

调节R0控制键长变化趋势，如果两个原子已经成键需要断掉就把R0设置成比当前键长大的距离，如C-C键1.5埃，需要断掉就设置成2.1埃，相反，两个没有成键的原子需要成键就可以将R0设置成成键键长。

作者
Author: 413 时间: 2019-10-30 20:42

youyno 发表于 2019-10-30 15:22
调节R0控制键长变化趋势，如果两个原子已经成键需要断掉就把R0设置成比当前键长大的距离，如C-C键1.5埃， ...

调节R0？就是初始结构的意思吗？
还有别的方法吗？

作者
Author: highlight 时间: 2019-10-30 20:59
本帖最后由 highlight 于 2019-10-31 08:03 编辑

413 发表于 2019-10-30 20:42
调节R0？就是初始结构的意思吗？
还有别的方法吗？

他说的应该是使用配位数作 CV 的情况：

https://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS/COLVAR/COORDINATION.html

你的情况或许是使用了键长作 CV，所以没有 R0。
至于通过 SCALE 的正负影响键长的想法，我认为从理论上看就不可能。
如果你确实需要限制 CV 的变化范围，WALL 关键词可能会有点用。

https://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/FREE_ENERGY/METADYN/METAVAR/WALL.html

作者
Author: 413 时间: 2019-10-30 22:21

highlight 发表于 2019-10-30 20:59
他说的应该是使用配位数作 CV 的情况：

你的情况或许是使用了键长作 CV，所以没有 R0。

请问WALL里面，POSITION指的是什么坐标？x？y？z？

作者
Author: highlight 时间: 2019-10-31 08:09

413 发表于 2019-10-30 22:21
请问WALL里面，POSITION指的是什么坐标？x？y？z？

CV，先看看这里吧
https://www.cp2k.org/exercises:2015_cecam_tutorial:mtd1

作者
Author: 413 时间: 2019-10-31 09:12

highlight 发表于 2019-10-31 08:09
CV，先看看这里吧
https://www.cp2k.org/exercises:2015_cecam_tutorial:mtd1

谢谢您提供的链接。想请教下比如2000步的fes结果和3000步的fes结果，是不是3000步的fes一定会比2000步的合理？

作者
Author: highlight 时间: 2019-10-31 10:52

413 发表于 2019-10-31 09:12
谢谢您提供的链接。想请教下比如2000步的fes结果和3000步的fes结果，是不是3000步的fes一定会比2000步的 ...

我只知道会不同，但是由于我没有跑出过成功的 MTD 轨迹，我也没法判断。但是从理论上看，符合最后结果收敛条件的话，应该是差不多的。

作者
Author: 413 时间: 2020-3-2 16:34

youyno 发表于 2018-8-6 09:48
graph.popt是cp2k编译时产生的用于metadynamics后处理自由能的程序，后缀根据版本不同有所区别，具体命令可 ...

请教下，metadynamics的作业怎么重启呢？我找了下手册，好像给的介绍都是利用restart文件得到自由能势能面的信息

作者
Author: youyno 时间: 2020-3-2 23:21

413 发表于 2020-3-2 16:34
请教下，metadynamics的作业怎么重启呢？我找了下手册，好像给的介绍都是利用restart文件得到自由能势能 ...

直接将.restart 文件复制成inp文件提交即可

作者
Author: 413 时间: 2020-3-3 14:40
本帖最后由 413 于 2020-3-3 14:43 编辑

youyno 发表于 2020-3-2 23:21
直接将.restart 文件复制成inp文件提交即可

收到

另外再请教下，从一个势阱出来，进入到另一个势阱时（此时两个势阱之间的势能面还没被填满）。
这个时候机器计算一个step花费的时间增长了很多很多，是正常现象吗？

作者
Author: youyno 时间: 2020-3-4 16:41

413 发表于 2020-3-3 14:40
收到
另外再请教下，从一个势阱出来，进入到另一个势阱时（此时两个势阱之间的势能面还没被填满）。
...

不正常，结构可能跑崩了，最好去看看轨迹检查

作者
Author: 413 时间: 2020-3-4 17:23

youyno 发表于 2020-3-4 16:41
不正常，结构可能跑崩了，最好去看看轨迹检查

谢谢，我检查了一下，结构确实崩了。请问是不是就是说整个计算都没有意义了？我是不是可以从崩之前的restart接着算？

作者
Author: youyno 时间: 2020-3-4 21:45

413 发表于 2020-3-4 17:23
谢谢，我检查了一下，结构确实崩了。请问是不是就是说整个计算都没有意义了？我是不是可以从崩之前的rest ...

崩之后的肯定没用了，重启时记得删掉COLVAR和HILLs里面崩掉的部份

作者
Author: 413 时间: 2020-3-6 14:15

youyno 发表于 2020-3-4 21:45
崩之后的肯定没用了，重启时记得删掉COLVAR和HILLs里面崩掉的部份

好的，谢谢

作者
Author: Flights 时间: 2020-3-11 10:40
借楼，求助，请问cp2k分子动力学计算之后的结构怎么进行电子结构分析呢？把结构提取出来拿vasp计算差分电荷、bader这样？cp2k可以计算这些嘛？谢谢！

作者
Author: sobereva 时间: 2021-2-10 03:40

Flights 发表于 2020-3-11 10:40
借楼，求助，请问cp2k分子动力学计算之后的结构怎么进行电子结构分析呢？把结构提取出来拿vasp计算差分电荷 ...

把轨迹中你要做电子结构分析的帧的结构提取出来。诸如你要做密度差，就对整体和相应片段分别做个单点计算并同时要求产生molden文件，然后按照Multiwfn手册2.9.2节所述将盒子信息加入molden文件中。之后在Multiwfn里按照下文的做法计算密度差格点数据并导出成cube文件。之后可以用VMD、VESTA等观看
使用Multiwfn作电子密度差图
http://sobereva.com/113

bader电荷没什么优点，严重缺点倒是一大堆。建议用Multiwfn载入CP2K的molden文件并计算CM5电荷，比bader电荷化学意义好多了。

作者
Author: 高斯不过如此 时间: 2021-7-30 21:16
借楼问一下，为什么CP2K用graph.popt算出来的FES是负值？这是什么意义呢？

作者
Author: a52241313 时间: 2021-9-1 15:57

xingyuan 发表于 2018-8-18 21:09
@youyno 老师，您好，又来麻烦您了。还是关于CP2K的问题，我有以下几个疑问：
1、CP2K Manual上关于计算 ...

"CP2K Manual上关于计算自由能有三种方法：Alchemical Change 、Metadynamics、Umbrella Integration。"请问这一部分的知识在网页的哪里去找，多谢！

作者
Author: elpa 时间: 2021-11-4 12:57

youyno 发表于 2018-8-6 09:48
graph.popt是cp2k编译时产生的用于metadynamics后处理自由能的程序，后缀根据版本不同有所区别，具体命令可 ...

您好，请问运行命令graph.popt -ndim 2 -ndw 1 2 -file cp2k-1.restart -find-path -point-a CVs -point-b CVs -cp2k 报错是为什么？
报错提示：forrtl: severe (59): list-directed I/O syntax error, unit -5, file Internal List-Directed Read 查了一下好像是输入文件的问题，但是restart文件是cp2k自己打印出来的我也没有改动过。

作者
Author: dhdhdh 时间: 2022-2-12 10:04

liujodan 发表于 2019-10-16 09:15
youyno老师您好，请问向我这个文件夹中，没有找到graph.popt,是不是没有编译上呢？还请指点，谢谢

请问问题解决了吗？我用的cp2k是7.0版本，输入graph.popt命令显示command not found。

作者
Author: zgdong 时间: 2022-3-23 11:24
请问您是怎么编译的，我用的也是popt版本，但没有graph相关文件

作者
Author: Vivat 时间: 2022-3-30 19:26

dhdhdh 发表于 2022-2-12 10:04
请问问题解决了吗？我用的cp2k是7.0版本，输入graph.popt命令显示command not found。

运行graph.popt的绝对路径

作者
Author: chen0201 时间: 2023-5-25 10:18

xingyuan 发表于 2018-9-13 14:28
感谢感谢。manual上面有一个例子是这样做的，一直不知道叫什么方法，学到了。那您知道计算自由能的umbrel ...

您好，我也是在学习cp2k计算pmf，请问您的参数设置是怎么样的呢？方便给我看看吗？

作者
Author: chen0201 时间: 2023-6-1 16:22

youyno 发表于 2018-8-7 16:52
你没有理解我的意思
命令应该是 graph.popt -ndim 2 -ndw 1 2 3 -file gr2hno3_mtd_3cv_p1-1_2000.restart ...

请问老师，cp2k没有graph.popt,采用的是graph.psmp -ndim 2 -ndw 1 2 3 -file C-CH3-cl-1.restart -cp2k命令，但是fes.dat是空的，请问怎么回事呢？

作者
Author: gehan 时间: 2023-6-2 12:42

youyno 发表于 2018-8-6 09:48
graph.popt是cp2k编译时产生的用于metadynamics后处理自由能的程序，后缀根据版本不同有所区别，具体命令可 ...

老师，您好。请问关于通过graph命令得到的fes.dat文件中前几列是CV值吗，为什么我得到的fes.dat文件中前几列的值和cp2k产生的COLVAR-metadynLog文件中的第二、三列的值不一样.

作者
Author: Coldwaterfish 时间: 2024-2-22 16:40

chen0201 发表于 2023-6-1 16:22
请问老师，cp2k没有graph.popt,采用的是graph.psmp -ndim 2 -ndw 1 2 3 -file C-CH3-cl-1.restart -cp2k ...

您好请问这个问题您解决了吗？我现在也遇到了同样的问题，盼回复

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