计算化学公社
标题: NAMD自由能计算教程—1、用eABF和meta-eABF进行多维自由能计算 [打印本页]
作者Author: fhh2626 时间: 2018-9-4 10:59
标题: NAMD自由能计算教程—1、用eABF和meta-eABF进行多维自由能计算
本帖最后由 fhh2626 于 2020-5-28 16:18 编辑
NAMD自由能计算教程—1、用eABF和meta-eABF进行多维自由能计算
fhh2626 Updated in 2020/5/28
本文首发计算化学公社分子模拟板块，禁止转载。
ChangeLog:
2020/5/28
well-tempered-meta-eABF的用法（建议使用这个）
NAMD 2.14b1及后续版本中边界力常数的设置方法
expandBoundaries选项
目前几大MD软件当中，Gromacs的中文资料比较丰富，而其他软件的则较少。尤其是NAMD，在网上很难找到中文教程。因此，作为NAMD的developer之一，我决定开一个坑，向大家介绍一下NAMD的中高级功能，可能包括但不限于TclForce，QM/MM，自由能计算，可极化力场，广义系综等等。希望我的教程能让大家觉得NAMD比别的MD软件要强大：）
这一系列教程假定读者已经具有一定的MD基础，能自己建立输入文件并且跑简单的平衡模拟（对于志在科研的你来说应该做得到吧。。。），由于我比较擅长自由能计算，所以我准备先介绍自由能计算这个主题，暂定有以下几个内容：
1、用eABF和meta-eABF进行多维自由能计算
2、谈谈metadynamics和其estimator
3、umbrella sampling、选择反应坐标和自由能计算当中的tricks
4、Colvars模块的高级应用
5、“炼金术”方法——FEP和TI
6、精确结合自由能计算和BFEE插件
下面进入主题。我把自由能计算方法大致分为两类：基于“炼金术”（alchemistry）的自由能计算方法和基于反应坐标的自由能计算方法。前者主要包括FEP和TI，由于其可以直接消失生长原子而不顾及化学原理，形似“炼金术”而得名。后者则包括ABF、Metadynamics、Umbrella Sampling和SMD四大方法。SMD是非平衡模拟方法，这里不多介绍，这个主题主要介绍其他的自由能计算方法。
如果要使用NAMD的话，则推荐使用ABF类方法，原因是NAMD的作者之一Chris Chipot和NAMD中自由能计算模块Colvars的作者Jérôme Hénin都是ABF阵营的兄弟。所以NAMD对ABF的支持较好。
在这个主题中，我会向大家介绍两种ABF的衍生方法，eABF和meta-eABF方法。其中eABF方法的基本原理是由Tony Lelièvre和Yang Wei首先提出的，由我最先实现在NAMD当中，之后Jerome提出了eABF的两个改进（czar estimator和pABF）。eABF目前是ABF的公认上位方法，即在任何情况下都应该使用eABF和不是传统ABF。meta-eABF是我最近提出的，（私货注意）是目前最快的自由能计算方法之一，（私货注意）我个人认为在任何情况下都可以使用meta-eABF代替eABF和metadynamics。
这里不详细讨论两种方法的原理，我会在第二讲和第三讲中涉及到各个方法的原理。此外，要详细了解原理或者引用参考文献的可以参考这几篇文献：
eABF：
Fu et al. J. Chem. Theory Comput., 2016, 12(8), pp 3506–3513
Lesage et al. J. Phys. Chem. B, 2017, 121(15), pp 3676–3685
meta-eABF:
Fu et al. J. Phys. Chem. Lett., 2018, 9(16), pp 4738–4745
Well-tempered-meta-eABF
Fu et al. Acc. Chem. Res. 2019, 52 (11), 3254–3264
下面以NANMA为例来介绍这两种方法的使用。NANMA又称alanine dipeptide，其构象变化通常可以用两个二面角来描述，这里我们将计算沿着这两个二面角的自由能变化：
(, 下载次数 Times of downloads: 164)
这里假定读者具有能自己生成大部分输入文件的基础，生成好的pdb，psf和namd文件已经提供在附件当中。
namd文件中有以下几句话：
colvars on
这句话表示激活NAMD当中的Colvars模块
colvarsConfig colvarsA.in
这句话表示Colvars模块的输入文件名为colvarsA.in
然后打开colvarsA.in文件
colvarsTrajFrequency 20000
colvarsRestartFrequency 1000000
上面两句分别定义了colvars文件记录数据的频率和输出的频率
colvar {
name phi
这里定义了一个叫phi的几何变量
width 5.0
该几何变量记录数据的精度
lowerboundary -180
upperboundary 180
记录的上下界
extendedlagrangian on
extendedFluctuation 5.0
extendedTimeConstant 200
定义体系中的扩展自由度，即eABF中的“e”，这里不用深究含义，对于任何体系均可设置extendedFluctuation=width和extendedTimeConstant=200
dihedral {
group1 {
atomnumbers 13
}
group2 {
atomnumbers 15
}
group3 {
atomnumbers 17
}
group4 {
atomnumbers 1
}
}
对于dihedral，即二面角colvar，需要定义4个原子（团），这里定义了13-15-17-1的二面角。定义原子（团）可以直接指定pdb当中的原子序号，也可以通过给一个reference pdb，在里边将某一列由0改为1实现，后面那种方法通常用于定义大原子团
}
colvar {
name psi
width 5.0
lowerboundary -180
upperboundary 180
extendedlagrangian on
extendedFluctuation 5.0
extendedTimeConstant 200
dihedral {
group1 {
atomnumbers 3
}
group2 {
atomnumbers 1
}
group3 {
atomnumbers 17
}
group4 {
atomnumbers 15
}
}
}
NAMD 2.14及之后的版本修改了边界力常数的定义方式，现在需要额外定义一个harmonicWalls block
#harmonicWalls {
# name mywalls
# colvars phi psi
# lowerWalls -180 -180
# upperWalls 180 180
# forceConstant 1.0
#}
分别定义phi psi的边界墙，力常数为1.0，注意单位是对应变量的width值，也就是说力常数会随着反应坐标种类/width不同而自动scale。1.0的力常数适用于大多数情况，在这里由于二面角是周期性的，所以无需定义边界力常数
abf {
定义跑的作业为ABF，因为之前设了扩展自由度，因此跑的是eABF
colvars phi psi
表示针对上面定义的phi和psi colvars进行自由能能计算
fullSamples 500
在每一个width区间内的平衡步数，500-2000适用于绝大多数情况
historyfreq 1000000
写history文件的频率
writeCZARwindowFile
使用czar estimator时产生一些中间文件
}
提交作业，跑50 ns以后（至少有一台工作站吧！），对eabf.czar.pmf作图。
(, 下载次数 Times of downloads: 124)
接下来尝试一下meta-eABF，meta-eABF是eABF和metadynamics的结合，打开colvarsA.in可以看到，配置文件绝大部分都跟eABF一样，多出了以下几句话，
subtractAppliedForce on
expandboundaries on
这是处理metadynamics高斯峰的一个选项，不写这个也不会影响计算结果，但是计算效率比较低
这是计算eABF偏置力的一个内部选项，使用meta-eABF时必须打开
metadynamics {
colvars phi psi
这里同时向定义的两个几何变量的扩展自由度上施加metadynamics偏置力，即meta-eABF
hillwidth 5
metadynamics的峰宽（单位为width），在meta-eABF中均可使用5
hillweight 0.1
metadynamics的峰高（kcal/mol），在meta-eABF中均可使用0.1
welltempered on
biastemperature 4000
最早的meta-eABF并没有引入well-tempered算法，因为当时我觉得没啥用。后来在复杂体系计算中发现加上well-tempered算法会是的计算收敛性更好，因为eABF本身是渐近收敛的，well-tempered-MtD也是，但是单纯的metadynamics并不是渐近收敛的。
Biastemperature定义了MtD渐近收敛的速率，4000对大多数情况适用
}
提交作业，运行10 ns，对metaeabf.abf1.czar.pmf作图。
(, 下载次数 Times of downloads: 156)
作者Author: sobereva 时间: 2018-9-4 11:26
“分子模拟论坛”是哪个论坛？
作者Author: fhh2626 时间: 2018-9-4 11:58
脑残了，我想说的是计算化学公社的分子模拟板块
作者Author: abdoman 时间: 2018-9-4 14:42
这个比较好！
学习学习。
作者Author: ruanyang 时间: 2018-9-4 19:25
顶顶，NAMD很不错
作者Author: pyscf 时间: 2018-9-5 00:49
请问fu博士：这个meta-eABF在namd的哪个版本之后可以用？谢谢
作者Author: fhh2626 时间: 2018-9-5 16:30
2.12好像不行，只有最新的Nightly version可以，需要自己编译，编译的方法在我的新帖里面
作者Author: 霜晨月 时间: 2018-11-13 14:06
谢谢楼主，感觉很有意思。
请教一下，如果是比较复杂的CV，比如 多个距离的平方和 或 RMSD的加权平均，可以直接用吗？如果不行的话，可以利用Plumed实现吗？
作者Author: fhh2626 时间: 2018-11-14 22:10
这些都不算复杂的CV，前者可以用colvars里面的custom function实现，后者可以用Linear and polynomial combinations of components实现。
如果是复杂的CV，可以用Scripted functions实现
作者Author: 霜晨月 时间: 2018-11-17 20:42
谢谢fu博士，看来NAMD比我想像的强大得多。
作者Author: 16aDream 时间: 2019-3-1 20:05
本帖最后由 16aDream 于 2019-3-1 20:17 编辑
Amber里面的ABMD模块是用abp计算的，这与abf或者eabf相比哪种更推荐？(好像amber没法用abf)另外我在使用ABMD的时候做了几个测试发现体系还是被局部极小值困住了，这是因为flooding time取的不合理吗，还是因为模拟时间不够长？
作者Author: fhh2626 时间: 2019-3-2 09:14
abmd是metadynamics的一个变种，陷入极小值有可能是模拟时间不够长，也有可能是反应坐标选得不好或者MtD的参数设得不好
作者Author: sixofgad 时间: 2019-3-14 15:08
目前还属于在学习阶段，感觉代码操作还是有点困难，不知道有没有图形界面的MD软件
学习使用呢
作者Author: fhh2626 时间: 2019-3-14 16:45
其实现在跑模拟已经很简单了，拿NAMD为例，跑模拟最少只需要3个文件，pdb（定义初始原子坐标），psf（定义原子连接方式、电荷和原子类型），prm（力场参数）和conf（模拟参数）文件，如果你跑的是生物体系（比如蛋白质、膜）的话，前面三种文件都可以用charmm-gui自动生成，你需要的就只是conf文件，这个文件的格式就类似下面的样子：
Structure xxx.psf
Coordinates xxx.pdb
Temperature 300 ;定义模拟温度为300K
......
你只要做完NAMD官网上任意一个tutorial，就可以用里面的conf文件作为模板，改几个参数适应你自己的体系就行了，任何图形界面的软件都只会把事情弄得更复杂
作者Author: sixofgad 时间: 2019-3-14 18:46
哦哦哦，这样的，不过NAMD的中文教程是真的少，有点基础教程什么的也好啊，楼主有木有兴趣做一个，感觉你是个大神哈哈哈
作者Author: fhh2626 时间: 2019-3-15 15:45
我倒是有这个想法。。。不过最近太忙了，没空弄
其实NAMD相比于其他软件，在教程方面做得已经很不错了，http://www.ks.uiuc.edu/Training/Tutorials/ ，这里就几乎有各方面的教程，建议从NAMD tutorial做起，很容易就能上手了
作者Author: pyscf 时间: 2019-3-22 06:33
希望楼主不要断更啊 这么有用的专题 顶上去！
作者Author: yjmaxpayne 时间: 2019-3-22 16:12
最近在玩Plumed+Gromacs，学习一下NAMD开发者的大作
作者Author: lonemen 时间: 2019-3-22 16:49
太好了，谢谢！也想学习
作者Author: xpyp 时间: 2019-3-23 00:06
请问，上述方法如eABF和meta-eABF能用GPU加速吗？
作者Author: fhh2626 时间: 2019-3-23 07:55
你用namd gpu版本运行就行啦，自由能计算本身是在cpu上完成的，但是md部分是在gpu上计算的，自由能计算所消耗的时间和md相比可以忽略不计
作者Author: xpyp 时间: 2019-3-23 09:54
好的，谢谢！顺便问下，NAMD中的FEP是否也是如此？
作者Author: fhh2626 时间: 2019-3-23 10:37
不fep不一样，fep目前发布的版本还不支持gpu，但是实际上gpu的fep代码已经开发完成了，你进namd的git，下里面的alchemy分支编译，就可以跑gpu版本的fep了
作者Author: xpyp 时间: 2019-3-23 10:56
many thanks
作者Author: springcute 时间: 2019-4-12 16:58
这个帖子太有用了，非常感谢楼主大神！
我最近一篇稿子，审稿人让补自由能计算，是一个70氨基酸的多肽，但三个螺旋之间有二硫键连接，请问：1）应该如何选择cv呢？2）二硫键的存在是否对自由能计算有影响？
作者Author: fhh2626 时间: 2019-4-12 17:56
你想计算描述什么过程的自由能？
作者Author: springcute 时间: 2019-4-12 18:05
抱歉没描述清楚。我想计算的是解折叠过程的自由能，分子的结构主体是3个螺旋和螺旋之间的loop，但螺旋间有二硫键来稳定结构。我曾经模拟过加热解折叠的过程，但不知道如何用eabf来计算其解折叠自由能呢？如果用多肽两端的原子距离作为cv的话，担心拉不断二硫键从而影响自由能的计算。
作者Author: fhh2626 时间: 2019-4-12 19:47
这种复杂情况需要试，感觉的话二维自由能计算应该可以，一维用end-to-end distance，另外一维用相对于完美结构的RMSD
作者Author: springcute 时间: 2019-4-13 21:01
谢谢您的建议，我尝试一下。
作者Author: xcandy 时间: 2019-5-11 22:00
请问老师，您这个图是拿什么软件做的，很养眼啊
作者Author: fhh2626 时间: 2019-5-12 17:20
结构图是vmd，自由能图是origin
作者Author: xcandy 时间: 2019-5-18 15:24
谢谢老师，我想算一节DNA的自由能（就是这个物质的能量），之前是用炼金术做的，做forward跟backward，一个体系得跑10来天，不知道能不能用meta-eABF。看了ABF的说明书，理解的不是很清楚，是要在Z轴上加牵引力吗
作者Author: xcandy 时间: 2019-5-22 09:55
已解决，自己摸索很多坑，看您之前说想总结教程，水平很菜，需要苦力我可以帮忙
作者Author: 朝歌夜弦 时间: 2019-5-27 16:59
请问fu博士，您好：我在跑完namd模拟之后，得到日志、坐标、轨迹文件，如何利用这些计算体系的结合自由能，体系是用namd和vmd完成的。
作者Author: fhh2626 时间: 2019-5-28 13:12
结合自由能就比较麻烦了，你具体是什么体系？一般来说有MMPBSA、FEP和PMF计算三种途径
作者Author: 朝歌夜弦 时间: 2019-5-28 19:35
您好，我做的是一个用VMD建立的一个全原子给予外部环境下某蛋白质体系分子动力学模拟，namd计算的，我可以通过什么方式计算它的结合能呢
作者Author: fhh2626 时间: 2019-5-29 16:41
要看你具体的体系，比如小分子大小
作者Author: robin 时间: 2019-7-25 16:40
不好意思想請教一下
如果要用NAMD QM/MM 計算蛋白質反應機制
目前想採用的是二維SMD將反應座標取出日後做Umbrella sampling
然而當我想對其中一個維度的距離固定時
發現到NAMD restraint.constraint的指令無法使用在QM區域的原子
請問我該怎麼做呢?
作者Author: fhh2626 时间: 2019-7-25 20:45
用colvars做约束是有效的，另外我觉得受限于模拟的时间，QM/MM+US不太可能做出好的自由能曲线
作者Author: xcandy 时间: 2019-8-26 11:26
老师，我看你的conf文件里，rigidbonds 设置的是 all；说明书里是 none。 请问这两个值对结果有影响吗
作者Author: fhh2626 时间: 2019-8-26 14:00
rigidbonds all表示用SHAKE/RATTLE和SETTLE算法限制氢原子的长度，否则氢原子振动频率过快，不能使用2fs的时间步长
作者Author: xcandy 时间: 2019-8-26 14:13
ABF说明书里有一句："be fully decoupled from any other constrained degrees of freedom",因为我自己的体系算出来PMF值很小（不正常的小）而且也加了SHAKE算法，所以想会不会是这个影响的
作者Author: fhh2626 时间: 2019-8-26 14:39
把完整的句子发出来我才能看懂是什么意思
不建议用ABF，ABF对反应坐标的要求比较tricky，如果你用了SHAKE/RATTLE的话，选择反应坐标时需要包含所有对应的氢原子，建议用eABF，不会出现这些问题
作者Author: xcandy 时间: 2019-8-26 14:55
本帖最后由 xcandy 于 2019-8-26 14:57 编辑
sorry..
"The ABF scheme assumes that the reaction coordinate, ξ, is fully unconstrained.The concept of unconstrained reaction coordinate imposes that ξ be fully decoupled from any other constrained degrees of freedom — viz. by means of the SHAKE or RATTLE algorithm."
上面所提到的反应坐标的限制只针对于ABF吗？eABF跟 meta-eABF不受影响吗？
作者Author: fhh2626 时间: 2019-8-26 17:50
ABF要计算沿反应坐标的力场力，要求反应坐标方向上不能有其他相互干扰的力（constraints, restraints, 其他维度的偏置力等等）
eABF和meta-eABF添加了扩展自由度，并且计算的时扩展自由度上的受力，所以没有这个问题
另外，伞状采样方法是一个基于histogram的方法，也不存在这个问题
作者Author: Daniel_Arndt 时间: 2019-10-5 10:05
我目前研究一个大环内酯类药物在丙酮中的构象，它有两个醚和两个羟基，形成了一个分子内氢键的“网”，我跑了几个ns，发现这个氢键的“网”始终没有被破坏掉。但结合之前的实验，我们推测它在溶剂中还应该有另一个构象，其中这个氢键的“网”被破坏了，取而代之的是它和溶剂形成的分子间氢键。我也刚刚开始接触enhanced sampling，我想问一下寻找这种构象适不适合用eABF处理？如果适合的话，有没有推荐的collective variable？因为我看过的文献中常常处理的是一些生物大分子，它们适用的collective variable不大方便照搬到这种小分子药物上。
作者Author: ene 时间: 2019-10-5 11:04
你这种跑跑加速动力学或者模拟退火估计也可以
作者Author: wuzhiyi 时间: 2019-10-5 19:10
分子动力学软件慢就是原罪 要是速度能敢上gromacs马上资料就多多了
再牛逼的enhance sampling也不能抵消因为跑的慢带来的收敛问题
作者Author: fhh2626 时间: 2019-10-5 20:14
本帖最后由 fhh2626 于 2019-10-5 20:17 编辑
NAMD+colvars比gmx+plumed快一些，主要是colvars比plumed快
增强采样和平衡模拟完全没有可比性，平衡模拟再快也就是几倍的事情，增强采样对收敛效率的提升是几个数量级的问题
另外就算不做自由能计算，gmx的速度在MD软件中也算不上快的，现在第一梯队是amber和openmm
作者Author: fhh2626 时间: 2019-10-5 20:15
我觉得你这种情况更适合用aMD，REMD这一类不需要选反应坐标的增强采样方法
作者Author: Daniel_Arndt 时间: 2019-10-6 09:24
多谢你的建议！
作者Author: s345538842 时间: 2019-10-22 21:12
付老师，您好，在使用eabf进行多维度的pmf计算时，虽然能够继续计算，但是总会出现以下提示的错误：
colvars: Integrated in 144 steps, error: 9.96932e-07
colvars: Integrated in 128 steps, error: 9.25099e-07
colvars: Integrated in 1 steps, error: 7.02279e-07
colvars: Integrated in 1 steps, error: 6.7063e-07
colvars: Synchronizing (emptying the buffer of) trajectory file "emeta_d-nc_573.colvars.traj".
colvars: Saving collective variables state to "rest.colvars.state".
colvars: Integrated in 1 steps, error: 5.93777e-07
colvars: Integrated in 1 steps, error: 5.7674e-07
用的是Lammps计算的，
版本号为：Initializing the collective variables module, version "2019-08-05".
Using LAMMPS interface, version "2019-08-01".
# collective
indexFile OH_wat.ndx
colvarsTrajFrequency 1000
colvarsRestartFrequency 10000
colvar {
name dist
width 0.2
lowerboundary 2.0
upperboundary 12.0
extendedlagrangian on
extendedFluctuation 0.2
extendedTimeConstant 200
subtractAppliedForce on
distance {
group1 { indexGroup Na }
group2 { indexGroup F }
}
}
colvar {
name ncoord
width 0.5
lowerboundary 4.5
upperboundary 13
outputValue on
extendedlagrangian on
extendedFluctuation 0.5
extendedTimeConstant 200
subtractAppliedForce on
coordNum {
group1 { indexGroup Na }
group2 { indexGroup OW }
cutoff 3.5
expNumer 6
expDenom 12
}
}
abf {
colvars dist ncoord
fullSamples 500
historyfreq 10000
writeCZARwindowFile
}
metadynamics {
colvars dist ncoord
hillwidth 5
hillweight 0.1
}
作者Author: fhh2626 时间: 2019-10-22 22:00
这个不是错误，error是误差的意思
作者Author: s345538842 时间: 2019-10-22 22:31
了解，非常感谢
作者Author: s345538842 时间: 2019-10-23 15:18
fu博士，您好，再请教个问题：利用eabf-meta 进行多维度PMF计算时，是否需要用abf_integrate进行后续的数据处理？
之前看到Giacomo回复过：If you have biased and sampled two collective variables, then by definition the potential of mean force is two-dimensional and it is the final result, requiring no post-processing :-) From this you can extract an infinite number of 1D PMFs, depending on the ones that you need. In most cases, you will calculate the Boltzmann exponential of the 2D PMF, which will give you the partition function of the system in the reduced two-dimensional space (assuming that the calculation was statistically converged).
https://lammps.sandia.gov/threads/msg81923.html
但是在colvars的colvars-refman-lammps文件中写明：In dimension 2 or greater, integrating the discretized gradient becomes non-trivial. The standalone utility abf_integrate is provided to perform that task.
且[size=13.3333px]Prof. Giacomo Fiorin在回答另外一个人时，也说是需要的：
That depends on what you are using: metadynamics will write the 2D PMF automatically; ABF will need post-processing for anything above 1 dimension; umbrella sampling will also require post-processing.https://lammps.sandia.gov/threads/msg70294.html
关于本篇文章中，您提供的eabf-meta的文档中也有abf_integrate的程序，但是您用的是metaeabf.abf1.czar.pmf作图。
我自己也尝试了 abf_integrate计算的结果与abf1.czar.pmf结果的对比，结果显示两者类似，但是在数值上存在大小的差异。
因此请问fu博士，利用eabf进行多维度PMF计算时，是否需要用abf_integrate进行后续的数据处理？
非常期待和感谢您的回复。
作者Author: fhh2626 时间: 2019-10-23 15:31
abf_integrate是一个老的积分器，现在colvars里面自带了一个新的积分器，所以直接用输出文件里面的pmf就可以，不过你想要用abf_integrate也没问题
作者Author: s345538842 时间: 2019-10-23 15:40
了解，非常感谢
作者Author: s345538842 时间: 2019-11-21 00:28
fu 博士，您好，再次打扰。
我现在想做分子在气液两相中的分配系数，因此需要使溶质分子从溶液中运动到气相里面。采用的是distanceZ的colvar方法进行了粗略的计算，具体文件如下：
计算过程中发现液相的水分子会随着溶质分子一起运动，不知为何，还请指教。
翻看论坛中的博文发现您以前提出过“1、建议约束水的质心，并且把PBC的盒子放大一些，不要让水跨过PBC边界，要不然colvars计算质心会出错”（http://bbs.keinsci.com/thread-12240-1-1.html）
但是因为水分子也可能会进入到气相中，所以我采用的是dummyAtom，使溶质分子分子的绝对坐标发生变化。
ColvarsTrajfrequency 100
ColvarsRestartfrequency 1000
colvar {
name distZ
width 2
lowerboundary 70
upperboundary 100
lowerwallconstant 20.0
upperwallconstant 20.0
extendedlagrangian on
extendedFluctuation 0.2
extendedTimeConstant 200
distanceZ {
main { atomnumbers 5998 }
ref { dummyAtom (20.0, 20.0, 0.0) }
axis (0,0,1)
}
}
abf {
colvars distZ
fullSamples 500
historyfreq 1000
writeCZARwindowFile
}
作者Author: fhh2626 时间: 2019-11-21 00:47
溶质如果是亲水的话，就有可能带着水分子走的啊，宏观上就是形成了共沸物
这种情况建议你用FEP计算溶解自由能
作者Author: s345538842 时间: 2019-11-22 00:32
fu博士，您好，非常感谢您的解答。我今天试了下把溶质分子放在水分子上面的真空层跑nvt，发现水分子是会被强烈的吸引过去的。
作者Author: minishotgun 时间: 2019-11-29 12:44
付老师您好，我想向您请教一个关于meta-eABF参数——hillwidth的问题
您在帖子里说
“metadynamics {
colvars phi psi
这里同时向定义的两个几何变量的扩展自由度上施加metadynamics偏置力，即meta-eABF
hillwidth 5
metadynamics的峰宽（单位为width），在meta-eABF中均可使用5
hillweight 0.1
metadynamics的峰高（kcal/mol），在meta-eABF中均可使用0.1
}”
我对”metadynamics的峰宽（单位为width）“这句话不是很理解，然后就去查了些资料
namd user's guide里说：
“hillWidth < Width of a Gaussian hill, measured in number of grid points >
Context: metadynamics
Acceptable Values: positive decimal
Description: The value of this keyword is the Gaussian width 2σξi, expressed in number of
grid points: the grid spacing is determined by width (see 9.3.8). The values of the parameter
σξi for each variable ξi inits physical units are printed by NAMD at initialization time. Values
between 1 and 3 are recommended: smaller values fail to adequately interpolate each Gaussian
[30], while larger values may be unable to account for steep free-energy gradients.”
namd官网上又说:
"hillWidth <Relative width of a Gaussian hill with respect to the colvar's width>
Context: metadynamics
Default value: √(2Π) / 2
Acceptable values: positive decimal
Description: The Gaussian width along each colvar, 2σ, is set as the product between this number and the colvar's parameter wi given by width (see 13.2.1); such product is printed in the standard output of VMD. The default value of this number gives a Gaussian hill function whose volume is equal to the product of Π×wi, the volume of one histogram bin (see 13.5.7), and W.
我的理解是如果cv的width = 5（单位°），hillwidth = 5，则添加的每个能量的hillwidth的值为5×5=25（单位°），请问这样理解对吗
我在做meta-eABF的时候，cv选的是distance z，但不慎把hillwidth设成0.5了，请问这个对结果影响大吗，是否可以通过延长模拟时间来弥补这个失误呀...
感谢老师指导
作者Author: sebastian6186 时间: 2019-12-16 15:10
谢谢分享
作者Author: chw985550192 时间: 2020-1-15 10:38
版主您好，非常好的帖子，非常感谢，有个问题想咨询下，我想以距离为反应坐标来计算自由能，那么 这两个值extendedFluctuation extendedTimeConstant 要怎么设置呢。
您也说了，对于任何体系均可设置extendedFluctuation=width和extendedTimeConstant=200
我的width=0.1，我设置了 extendedFluctuation 0.1 ，extendedTimeConstant 200，发现自由能不收敛，是不是这个参数设置不太好，还是有其它的参数设置需要注意，导致求出的能量不收敛，您可以回答下吗，谢谢。
作者Author: fhh2626 时间: 2020-1-15 17:05
自由能计算不收敛具体是指什么？应该不是这两个参数的问题
作者Author: lonewolf659 时间: 2020-3-31 10:56
本帖最后由 lonewolf659 于 2020-3-31 11:01 编辑
楼主，你好！请教一个问题。如果我想计算一个protein在2个conformation间变化的自由能势能面，应该如何选择free energy的计算方式，应该怎样定义CV啊？
NAMD是否有类似的功能，使用某个structure作为target reference，然后从任意初始structure开始，计算自由能变化？
谢谢！
作者Author: fhh2626 时间: 2020-3-31 14:59
可以选RMSD作为反应坐标，但是可能不太容易收敛
作者Author: lonewolf659 时间: 2020-3-31 17:39
谢谢楼主，我试试看
作者Author: youyno 时间: 2020-4-1 12:19
http://bioweb.cbm.uam.es/software/MEPSA/
这个程序更好，可以找到最小能量路径，局部最小值等等
作者Author: fhh2626 时间: 2020-4-1 21:59
这个是自由能计算完了，得到自由能面以后找LFEP的东西
作者Author: huangjin 时间: 2020-4-1 22:29
这个很好，学习一下，非常感谢！
作者Author: azero 时间: 2020-4-27 23:59
本帖最后由 azero 于 2020-4-28 09:09 编辑
请教在plumed中，怎么设置eABF的width
1. 是不是用公式 (upperboundary - lowerboundary) / width计算，然后把该值填进 GRID_BIN项？
2. 此外，还看到 ZGRID_BIN，搞不清 ZGRID_BIN和 GRID_BIN的区别3. KAPPA是不是就是力常数？看看定义似乎是，但有些例子是KAPPA=418.4 ，418.4我不常见，然后就不太确定了
作者Author: fhh2626 时间: 2020-4-28 09:14
可以不填写GRID_BIN项，width填写到GRID_SPACING就行了，要用GRID_BIN可以用(上边界-下边界)/bin宽度计算，不需要平方
一般不需要填写ZGRID_*项，会自动和GRID_*保持一致，GRID是给平均力用的，ZGRID是给estimator用的，边界设置可以不同，但一般都保持一致。
作者Author: azero 时间: 2020-4-28 23:12
谢谢付老师的解答~
今天尝试用gromacs+plumed 跑了一下eABF,然后有很多问题请教付老师你:
1.怎么判断PMF已经收敛？
查了一下，说延长模拟时间，如果PMF曲线变化不大即可；或者grad文件作的图是连续的即可(这个没看明白什么连续)。
这里是eabf_winall.abf系列文件作的图，请较一下付老师, 模拟是否收敛了。
(, 下载次数 Times of downloads: 85)
2.假如我想续跑一个模拟，怎么对drrstate 文件合并？
看到手册里有plumed drr_tool --merge，但似乎不是合并轨迹的。
作者Author: fhh2626 时间: 2020-5-2 18:02
最好是count比较均匀，grad连续的意思是相邻的两个点的梯度不要有太大的突变
作者Author: Lacrimosa 时间: 2020-5-26 13:44
Fu老师，您好，我想请教一下，我看到有些文献和教程中用rdf根据公式w=-kTln（g）算出来了PMF，请问这种方法应该归为您提到的几类中的哪一类呢？一般适用于什么样的问题呢？我想计算矿物孔道中吸附离子的自由能变化用哪些方法比较好呢？
作者Author: fhh2626 时间: 2020-5-26 18:32
那是用平衡模拟就能很好的研究的情况，直接统计平衡模拟中探测到各种态的概率，只有及其简单的情况才会用
作者Author: Lacrimosa 时间: 2020-5-26 19:54
感谢您的指导~
作者Author: minishotgun 时间: 2020-6-2 10:57
付老师您好，目前我们在用meta-eABF方法进行模拟，但自由能面一直不收敛（如图所示）。请问我们可以在脚本里加上well-tempered-meta-eABF的参数，让我们现在在跑的meta-eABF变成well-tempered-meta-eABF吗（我们担心这样做会引起计算误差，导致结果不精确）。
感谢老师指导。
作者Author: fhh2626 时间: 2020-6-2 16:28
不能直接变，但是你可以把hillheight和hillwidth改小一些（比如0.01和3）
一般来说不收敛是因为你跑的体系本身就非常复杂，需要很长时间的模拟
作者Author: minishotgun 时间: 2020-6-2 17:34
谢谢老师，那我把这两个参数改小一点。我跑的体系是比较复杂，我再多跑一段时间看看。
作者Author: s345538842 时间: 2020-9-20 09:24
本帖最后由 s345538842 于 2020-9-20 09:35 编辑
(, 下载次数 Times of downloads: 126)
DRR ...
LABEL=eabf
ARG=ncoord
FULLSAMPLES=500
GRID_MAX=1.2
GRID_MIN=0.0
# GRID_BIN=60
GRID_SPACING=0.01
FRICTION=10.0
TEXTOUTPUT
OUTPUTFREQ=5000
HISTORYFREQ=50000
TAU=0.5
TEMP=300
...
METAD ...
LABEL=mtd
ARG=eabf.ncoord_fict
SIGMA=0.1
GRID_MAX=2.0
GRID_MIN=-0.2
GRID_SPACING=0.01
# GRID_BIN=60
HEIGHT=1.3
PACE=500
TEMP=300
BIASFACTOR=8
...
uwalls: UPPER_WALLS ARG=eabf.ncoord_fict AT=1.2 KAPPA=100000 EPS=3
lwalls: LOWER_WALLS ARG=eabf.ncoord_fict AT=0 KAPPA=100000 EPS=3
付老师，您好，我用plumed+lammps计算配位数变化的自由能变化，但是在计算过程中 COLVAR存在剧烈波动，结果就是eabf.czar.pmf 存在巨大问题，参数设置见以上，请问问题出在哪里？单独使用eabf也存在这个问题，使用WTmeatd 就没有这个问题。
作者Author: fhh2626 时间: 2020-9-22 11:20
我感觉参数设置都挺对的，没看出来问题在哪，你用WT-MtD算出来的结果合理吗？
话说Lammps用自带的Colvars不香吗，怎么还外挂了一个Plumed
作者Author: eagletyr 时间: 2020-11-2 21:22
求教Fu老师，我主要用gromacs来做蛋白质和小分子的体系。1.我用拉伸动力学产生伞状采样的初始构象，拉伸时对弹簧常数以及拉伸速度选取有什么规则么？多谢。由于我采用的是恒速拉伸，RC vs time做图时，看到的图并不是恒速的（刚开始由于小分子与蛋白质作用强，反应坐标RC变化比设的拉伸速度慢很多，后面比较正常），此时是否需要加大弹簧常数以使得RC vs time做图近似一条直线（表示拉伸速度恒定）？2. 我在伞状采样时对于弹簧常数的选取又有什么要注意的么？只知道各相邻窗口要重叠的好，这个是要多试几次？多谢。
作者Author: fhh2626 时间: 2020-11-3 10:28
1、比较随意，你只是要个初始结构而已。另外恒速拉伸不代表体系的真实RC是恒速的
2、只能是感觉了
作者Author: eagletyr 时间: 2020-11-3 11:14
本帖最后由 eagletyr 于 2020-11-3 11:18 编辑
好的，多谢老师。另外想请教一下，假设我从 classical MD（100 ns）的轨迹里取了5个结构出来，从这5个结构出发分别来拉伸并伞形采样，请问的是这5个结构从轨迹取出来后，是要经过能量最小化、预平衡、拉伸的过程还是直接拉伸就可以？多谢。
作者Author: fhh2626 时间: 2020-11-3 13:43
直接拉伸
作者Author: eagletyr 时间: 2020-11-3 21:46
好的，多谢FU老师
作者Author: s345538842 时间: 2020-12-1 11:12
总是充满了无奈
，最近又试了试，发现把MD的计算步长调小一点就可以了，谢谢傅老师
作者Author: yang001 时间: 2020-12-1 17:12
你好，我按照 https://pubs.acs.org/doi/abs/10.1021%2Facs.jcim.7b00695 和 https://github.com/fhh2626/Binding-free-energy-estimator-BFEE-要求安装好了BFEE，但是总是出现这个错误，请问该怎样解决呢？error copying "C:/Users/yangk/Downloads/ci7b00695_si_002/p41-abl_my/bound.psf" to "./bound.psf": permission denied
error copying "C:/Users/yangk/Downloads/ci7b00695_si_002/p41-abl_my/bound.psf" to "./bound.psf": permission denied
while executing
"file copy -force $psfDir ./bound.psf"
(procedure "copy_files" line 19)
invoked from within
"copy_files"
(procedure "::BFEE::generateFiles" line 3)
invoked from within
"::BFEE::generateFiles"
invoked from within
".bfee.setupTab.setup.generate invoke "
invoked from within
".bfee.setupTab.setup.generate instate {pressed !disabled} { .bfee.setupTab.setup.generate state !pressed; .bfee.setupTab.setup.generate invoke } "
(command bound to event)
作者Author: fhh2626 时间: 2020-12-1 18:25
VMD的当前路径不要和输入文件放一个文件夹
作者Author: beyond 时间: 2020-12-11 06:28
刚刚用Gromacs+Colvars 跑了一个meta-eABF, 请教楼主一个问题：
如何判断模拟收敛了呢？貌似模拟过程中，直接输出了PMF，并且备份前面输出的一个。
但是再往前的PMF就丢失了，并不能像PLUMED那样，可以给出任意模拟时间的PMF。
不知道楼主是如何监测模拟进程的 ？谢谢！
作者Author: fhh2626 时间: 2020-12-11 10:04
如果你设置了historyFreq的话，colvars会输出一个.hist.czar.pmf文件，记录所有的pmf
作者Author: beyond 时间: 2020-12-11 19:14
谢谢楼主！第一眼看到hist，还以为是histogram呢
我再看看Colvars的manual！
作者Author: yang001 时间: 2020-12-16 11:27
本帖最后由 yang001 于 2020-12-16 11:32 编辑
你好，我按照“BFEE: A User-Friendly Graphical Interface Facilitating Absolute Binding Free-Energy Calculations”这篇文章附件中“p41-abl”例子采用bfee v1.1生成了模拟初始文件并顺利模拟（constant T控制采用langevin动力学，NAMDv2.13）完毕，但是计算的自由能（采用默认参数）为+1.23kcal/mol，请问是什么原因呢？
作者Author: fhh2626 时间: 2020-12-16 13:40
BFEE只是把输入文件做好，但是很多具体的模拟参数还是要自己调一下，比如默认的模拟步数第1、7、8步（很可能是第7步）很可能收敛得不是很好，要分窗口并且延长模拟时间
作者Author: yang001 时间: 2020-12-16 15:26
请问这种PMF随时间的变化曲线是怎么做出来的
作者Author: yang001 时间: 2020-12-17 09:49
请问您说的分窗口是如何实现的？比如对于007_r来说，将原来的 {lowerboundary 10.0； upperboundary 33.0}改为 {lowerboundary 10.0； upperboundary 20.0} ， {lowerboundary 16.0； upperboundary 26.0} ， {lowerboundary 22.0； upperboundary 33.0}三个分别模拟吗？那么后续结果如何处理呢？
作者Author: fhh2626 时间: 2020-12-17 12:30
输出有个.hist.pmf文件，你算每一帧和最后一帧的RMSD差别，或者算每一帧和0向量的差别，就可以得到PMF RMSD随时间的变化
分窗口的话，比如第7步距离范围是10-35，你可以分别跑范围10-13,13-16,16-19,19-35，然后把输出的PMF加减一个常数，头尾连接起来，因为靠近蛋白质的部分比较复杂，所以靠近蛋白质的部分可以把窗口分小一些
作者Author: yang001 时间: 2020-12-17 15:15
还有一个问题，我在模拟自己的小分子和蛋白质之间结合自由能时，在170ns的模拟时长内，从RMSD_unbound2.abf1.count中可以看到0-1.175A之内的count数为0，而在RMSD_bound2.abf1.count中的数目也比其他RMSD区间少20倍，这种情况对结果影响大吗？需要采取什么措施来解决这个问题吗？
作者Author: fhh2626 时间: 2020-12-17 15:45
这个正常的，因为几何熵的作用导致RMSD=0的自由能接近于无限大，所以那个点不会有采样的
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |