NAMD自由能计算教程—1、用eABF和meta-eABF进行多维自由能计算

JCenter · 发表于 Post on 2024-6-23 21:39:34

本帖最后由 JCenter 于 2024-6-23 21:42 编辑

请教下付博：我看完了colvars在gmx的使用手册，有几个疑问，想请您有时间指点一下

：
1.内置在gmx的colvars的距离单位是跟随gmx的nm还是默认的埃。我看手册中一会儿是nm,一会儿是埃。
width 0.1
  lowerboundary 0.01
  upperboundary 1.40
2.关于eABF、well-tempered-meta-eABF等方法的使用。耦合了well-tempered算法和Mtd的eabf可以用于一维PMF的计算吗。我看文章里，这种耦合多用于二维自由能面的计算。
3.eabf偏置方法中的fullSamples N我没有理解的。每个width区间内的平衡步数。是指每个width平衡N步后采集后面的数据，还是说仅采集平衡后的N步的数据呢
4.well-tempered-meta-eABF方法中，  高斯峰的峰宽（3）和和峰高（0.1KJ/mol）值会明显影响PMF曲线吗,因为well-tempered-meta-eABF方法，自由能计算是填谷（well-tempered-meta）和削谷（eabf）同时进行。如果单单使用Mtd,必然影响结果得数值精度。如果是两种方法呢。
5.简谐墙的上下边界值与变量的上下阈值一致，会导致上下两个值采集不到样吗，还是说，可以把力常数设的小一点，比如5 KJ/mol?
6.关于总时间的设定，我该怎么考虑呢。假设每个width20ns,我有14个width，因此总时间可以设为280ns?
感谢您的指导

fhh2626 · 发表于 Post on 2024-6-26 18:56:00

1、印象中好像是nm
2、可以，几维并没有什么本质区别
3、没看懂你的问题，你说的两句话不是一个意思吗？
4、不会影响收敛后的结果，有可能会影响效率，单位是kJ的话可能用0.4更好
5、不会，建议用reflecting boundary
6、越长越好

fhh2626 · 发表于 Post on 2024-6-26 18:57:05

JCenter 发表于 2024-6-23 21:39
请教下付博：我看完了colvars在gmx的使用手册，有几个疑问，想请您有时间指点一下：
1.内置在gmx的co ...

见上面的回复

JCenter · 发表于 Post on 2024-6-27 13:50:10

fhh2626 发表于 2024-6-26 18:56
1、印象中好像是nm
2、可以，几维并没有什么本质区别
3、没看懂你的问题，你说的两句话不是一个意思吗？
...

好嘞，感谢付博的详细的指导

。
问题3可能没表述好，我重新换下表述：假设在每个width下设了fullSamples 5000，这个5000是指在CV的每个width范围内先跑个5000步，以平衡该区间的过程，然后从5000步的平衡后开始计算平均力，其计算平均力的步数并非5000这个值，从而计算自由能；还是指在CV的每个width范围内，系统自己判定平衡过程，默认平衡完后，再采样5000次，对这5000步进行平均力计算，从而计算自由能呢。主要是没理解每个width区间内的平衡步数，这个值在模拟计算过程中的物理意义。
问题6，这个自由能有无判定的收敛标准吗。还是说，对于一个体系，我就跑个不同时间模拟下的PMF,对比曲线变化？感觉这种参考判据不太好掌控，得不断地换时间。

fhh2626 · 发表于 Post on 2024-6-28 14:47:38

JCenter 发表于 2024-6-27 13:50
好嘞，感谢付博的详细的指导。
问题3可能没表述好，我重新换下表述：假设在每个width下设了fullSam ...

在每个bin里面平衡5000步，这5000步仅收集力，5000步以后才开始施加偏置力（实际实现略有不同）

设置historyfreq以后会输出不同时刻的自由能曲线/面，可以用来判断收敛性

JCenter · 发表于 Post on 2024-7-3 02:55:52

fhh2626 发表于 2024-6-28 14:47
在每个bin里面平衡5000步，这5000步仅收集力，5000步以后才开始施加偏置力（实际实现略有不同）

设置h ...

感谢付博，言简意赅

docshen777 · 发表于 Post on 2024-7-8 18:02:35

付老师，想请教您一些问题，我想用meta-eabf方法采样我的体系，我的plumed的输入文件如下

DRR ...
LABEL=eabf
ARG=Ob,Nc
FULLSAMPLES=500
GRID_MIN=-5,-5
GRID_MAX=40,20
GRID_BIN=180,100
OUTPUTFREQ=1000
HISTORYFREQ=1000
TEXTOUTPUT
... DRR

METAD ...
LABEL=mtd
ARG=eabf.Ob_fict,eabf.Nc_fict
HEIGHT=15
SIGMA=1,1
BIASFACTOR=1500
TEMP=300
PACE=1000
FILE=HILLS.txt
GRID_MIN=-5,-5
GRID_MAX=40,20
GRID_BIN=180,100
GRID_WFILE=GRIDS
GRID_WSTRIDE=1000000
... METAD

PRINT STRIDE=1000 FILE=Colvar.data ARG=*

我对此有一些疑问：
1. 我参照您的讲座和博文的内容，设定的eabf和meta的间隔分别为500步和1000步，这是否存在不协同的问题？eabf进行了2次，meta进行一次这样是可以的吗？是否会有影响？
2. 输入文件中，是把eabf偏差后的集合变量拿给meta去做添加高斯峰，而我们统计势能面时，却只使用eabf的偏置势，这是为什么呢？这也是最困惑我的地方。
3. 若是只需要eabf的偏置势，META部分是否还需要输出grids？还是在eabf的地方输出？

谢谢付老师，希望您有空能解答一下学生的疑问。

fhh2626 · 发表于 Post on 2024-7-8 20:28:43

docshen777 发表于 2024-7-8 18:02
付老师，想请教您一些问题，我想用meta-eabf方法采样我的体系，我的plumed的输入文件如下

DRR ...

1. fullsamples不是你理解的意思，eABF每步都会更新内部的平均力
2. 不是只使用eABF的偏置势，而是使用eABF模块中的estimator，严格来说这个estimator适用于所有使用扩展自由度的自由能计算方法

docshen777 · 发表于 Post on 2024-7-8 21:13:04

fhh2626 发表于 2024-7-8 20:28
1. fullsamples不是你理解的意思，eABF每步都会更新内部的平均力
2. 不是只使用eABF的偏置势，而是使用e ...

付老师，谢谢您的解答。关于您提到的estimator这些概念，请问哪里有资料可以学习到吗？我想去深入了解

fhh2626 · 发表于 Post on 2024-7-9 10:05:23

docshen777 发表于 2024-7-8 21:13
付老师，谢谢您的解答。关于您提到的estimator这些概念，请问哪里有资料可以学习到吗？我想去深入了解

https://pubs.aip.org/aip/jcp/art ... rectedFrom=fulltext
https://pubs.acs.org/doi/10.1021/acs.jctc.6b00447
https://pubs.acs.org/doi/10.1021/acs.jpcb.6b10055

docshen777 · 发表于 Post on 2024-7-9 13:01:27

fhh2626 发表于 2024-7-9 10:05
https://pubs.aip.org/aip/jcp/article-abstract/123/14/144104/905291/Bridging-the-gap-between-thermo ...

谢谢付老师。我采用plumed2.8.4版本+lammps来使用WTM-eabf，但报错如下：

ERROR
I cannot understand line: DRR LABEL=eabf ARG=Ob,Nc FULLSAMPLES=500 GRID_MIN=-5,-5 GRID_MAX=40,20 GRID_BIN=180,100 OUTPUTFREQ=1000 HISTORYFREQ=1000 TEXTOUTPUT

但是我编译时的确装了drr包，使用plumed config show命令显示如下：# list of modules
# syntax: module name on/off (default-on/default-off/always)
module adjmat on (default-off)
module analysis on (default-on)
module annfunc on (default-off)
module asmjit on (always)
module bias on (default-on)
module blas on (always)
module cltools on (default-on)
module colvar on (default-on)
module core on (always)
module crystallization on (default-off)
module dimred on (default-off)
module drr on (default-off)
module eds on (default-off)
module fisst on (default-off)
module function on (default-on)
module funnel on (default-off)
module generic on (default-on)
module gridtools on (always)
module isdb on (default-on)
module lapack on (always)
module lepton on (always)
module logmfd on (default-off)
module manyrestraints on (default-off)
module mapping on (default-on)
module maze on (default-off)
module molfile on (default-on)
module multicolvar on (default-on)
module opes on (default-off)
module pamm on (default-off)
module piv on (default-off)
module reference on (always)
module s2cm on (default-off)
module sasa on (default-off)
module secondarystructure on (default-on)
module setup on (default-on)
module tools on (always)
module vatom on (default-on)
module ves on (default-off)
module vesselbase on (always)
module xdrfile on (default-on)

请问这是为什么呢？

fhh2626 · 发表于 Post on 2024-7-10 11:36:03

docshen777 发表于 2024-7-9 13:01
谢谢付老师。我采用plumed2.8.4版本+lammps来使用WTM-eabf，但报错如下：

是不是不能多行设置？plumed有些版本有bug，这种换行的写法识别不了

docshen777 · 发表于 Post on 2024-7-11 10:32:12

fhh2626 发表于 2024-7-10 11:36
是不是不能多行设置？plumed有些版本有bug，这种换行的写法识别不了

ERROR
I cannot understand line: DRR LABEL=eabf ARG=Ob,Nc FULLSAMPLES=500 GRID_MIN=-5,-5 GRID_MAX=40,20 GRID_BIN=180,100 OUTPUTFREQ=1000 HISTORYFREQ=1000 TEXTOUTPUT
Maybe a missing space or a typo?

写成一行依然不行哎

yyf123 · 发表于 Post on 2024-8-21 13:18:53

以下是邮件询问付博士的内容，供大家讨论

问题：

colvar {

name eulerTheta

eulerTheta {

      atoms {

         indexGroup  ligand

         centerReference on

         rotateReference on

         centerToOrigin    on

      enableFitGradients on

         fittingGroup {

            indexGroup  protein

         }

         refpositionsfile  ./complex_largeBox.xyz

      }

      refpositionsfile  ./complex_largeBox.xyz

}

}

1.  以上是bfee2产生的内容，我研究了很久，看了手册，实在是不能理解centerReference（我理解的就是对齐complex_largeBox.xyz 这个文件里面的ligand的意思，但是既然是对齐这个ligand，如图（例子文件里面的膜蛋白），为什么不是重合在一起的呢），另外rotateReference， centerToOrigin，fittingGroup是什么意思我也是完全看不懂，希望能指点一二

2. 那个large_box的文件，只是增大了盒子然后多了水分子，但是磷脂分子没有增多？磷脂会跑散吧？

3. 延续第一个问题，我发现只有距离变量没有centerReference之类的参数，这是为何？

付博士回：
fittingGroup是计算变量前对体系进行一个旋转平移。如果设置为蛋白的话，就是在计算这个角度前，对体系进行一个旋转平移，使帧与帧之间蛋白不存在相对运动。因为欧拉角与小分子的取向有关，你也可以理解为设了这个选项以后就是小分子相对于蛋白的取向。

前面那几个选项定义了计算四元数/RMSD的方法，默认就是on。例如计算RMSD时，先要对小分子进行一个最优旋转，再对小分子进行一个最优平移，消除小分子的平动和转动。欧拉角就是根据这个最优旋转求出来的。

yyf123 · 发表于 Post on 2024-8-21 13:31:44

本帖最后由 yyf123 于 2024-8-21 13:33 编辑

邮件向付博士的提问内容：

我使用2024gromacs自带的colvars计算PMF，钠离子从箭头处（水里面），通过离子通道(如图)，加了圆柱势（如图），但是！不知道为什么即使钠离子在水里面的PMF十分的曲折（如图）（我看您发的wtm-eabf的那个帖子下面很多人做的pmf都不平滑, 不知道为啥），19日之后讨论结果我会上传论坛供大家学习

1. 是因为圆柱势能的影响吗？
2. 能否把圆柱势能用reflectingUpperboundary替代？

付博士回：
可能是你算的时间太短了？你跑了多久？在模拟刚开始有噪声是正常的，一般要一两百纳秒可能才会比较光滑。

另外hillWidth这个参数设为偶数的话物理意义有点不明确，建议改为3或者5，hillWeight也可以适当降低到0.5左右。
还有你这个圆柱是否太宽了？可能要确认一下。

解决办法：

FullSamples改为100000，biasTemperature为2791, hillWidth为5.0, hillWeight为0.7 kj/mol，newHillFrequency为500
仅供大家参考

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[NAMD] NAMD自由能计算教程—1、用eABF和meta-eABF进行多维自由能计算