自适应偏置力的分解

xiaochenxi123

请问下面的每一步怎么做啊，实在是搞不懂，谢谢
看文章像是用adaptive biasing force (ABF)算的PMF。我查了一下，ABF方法的作者Chipot以前在NAMD mailing list 里回复了该怎样对结果进行处理来分解PMF，具体在这里
http://www.ks.uiuc.edu/Research/ ... 2009-2010/1566.html

他说
(1) loop over your trajectory and calculate the pair interactions
       that you are interested in, e.g. solute 1 and solute 2.
(2) collect the forces at play, van der Waals and electrostatics.
(3) project these forces onto your reaction coordinate, which for
       each configuration, should be determined independently.
(4) average the projected forces.
翻译成中文大概就是
1. 读取轨迹文件，对于每一帧，计算粒子1和2之间相互作用力
2. 把1和2之间的相互作用力分解成vdw和静电两部分
3. 把计算的力投影到反应坐标上。这句话的实际意思是，你会把
反应坐标分成多个bin（比如在0到10之间，每隔0.1设一个bin），假如
你读取的这帧构型反应坐标为2.0， 就把vdw和静电力的数值放到2.0 这个bin里
4. 如此反复直到处理完这个轨迹，然后对每个bin里力的数值加和求平均，就得到
不同反应坐标值时的vdw和静电力
5. 再对这些vdw和静电力沿反应坐标积分，即得到两部分力对PMF的贡献。

我不清楚NAMD或者VMD里是否有现成的命令做这个，否则的话你得自己写个小程序来做这个。
第四步中求平均值好像用matlab，其他的也不懂

fhh2626

大概是这样的：
先跑ABF(推荐用我们和Chipot新开发的eABF代替），得到一个在反应坐标空间均匀分布的轨迹文件
然后对dcd做pair interaction，NAMD里面有现成的命令，网上可以google到教程
注意，pair interaction可以分多次跑，比如第一次算主体和客体的interaction，第二次算主体和溶剂，第三次算客体和溶剂，这样就把体系受力分成了三个组分
pair interacton会输出log文件，里面会输出体系所受的静电力和范德华力，将这些力对应到每个小bin中，对每个小bin中的力求平均，然后积分
上述过程可以把自由能分解为6个部分：
主体-客体静电，主体-客体VDW，主体-溶剂静电/VDW，客体-溶剂静电/VDW

做自由能分解还是很复杂的，需要对MD和PMF计算都有一定的了解才行，不推荐新手上来就做这个（上面每一步都可以写几页的教程.。。）

xiaochenxi123

fhh2626 发表于 2017-3-28 11:32
大概是这样的：
先跑ABF(推荐用我们和Chipot新开发的eABF代替），得到一个在反应坐标空间均匀分布的轨迹文 ...

主要是对每个bin中的力求平均不会

fhh2626

xiaochenxi123 发表于 2017-3-28 16:12
主要是对每个bin中的力求平均不会

比如说你的轨迹里面有10帧，那pair interaction的结果是
1: elec: xxx, vdw: xxx
2: elec: xxx, vdw: xxx
...
10: elec: xxx, vdw: xxx

然后你的反应坐标范围为0-3，分成了3个bin（宽度就是你跑自由能计算的width）: 0-1,1-2,2-3

你需要用你的VMD写一个脚本，判断每一帧归属于哪个bin，比如你的反应坐标是距离，你通过脚本判断出第1帧A-B分子的距离是1.6，那就归属到第二个bin中去

这样你的结果就是
bin1（0-1）： (frame 1, elec:xxx, vdw:xxx) (frame 4, elec:xxx, vdw:xxx) (frame 10, elec:xxx, vdw:xxx)
bin2（1-2）: .......
bin3:.............

下一步就是求对力在反应坐标的投影，比如第一帧，vdw力向量是(1,2,3)，但是A和B的距离向量是(0,0,1)，你就要求(1,2,3)在(0,0,1)方向上投影的长度

下一步就是求每个bin里面力的平均，积分

说得有点混乱。。不知道你看懂了没

xiaochenxi123

fhh2626 发表于 2017-3-28 17:37
比如说你的轨迹里面有10帧，那pair interaction的结果是
1: elec: xxx, vdw: xxx
2: elec: xxx, vdw: x ...

你能详细些吗，实在是笨鸟，不是很懂，能不能有个其他联系方式，谢谢

xiaochenxi123

xiaochenxi123 发表于 2017-3-29 10:30
你能详细些吗，实在是笨鸟，不是很懂，能不能有个其他联系方式，谢谢

主要是方便我请教问题谢谢

fhh2626

xiaochenxi123 发表于 2017-3-29 11:09
主要是方便我请教问题谢谢

就在论坛里面说吧，这样别人也能看见。你ABF和pair interaction都没问题，就是下一步不会对吧