教程翻译----来自窗口采样的甲烷-甲烷之间的PMF

翻译、修改并重复了一篇关于计算PMF的教程。。

原文链接:GROMACS Tutorial 5 - Methane-methane PMF from window sampling

注意:原网站需要梯子才可访问。

译文如下

通过窗口采样计算甲烷-甲烷分子之间的PMF

设置

创建盒子

首先得到甲烷分子结构和拓扑文件,这里我也不再详述,前博文已经多次提及。在这里,我将采用amber99sb-ildn去模拟。方形盒子大小至少为边长3.1 nm, 可以通过gmx insert-molecules将两个甲烷分子插入到盒子中,注意这两个甲烷分子不要相离太远,可以用vmd将分子靠近。

创建索引文件

可以通过命令gmx make_ndx -f conf.gro得到,分别依次输如:a 2name X CA可以将原子编号为2的碳原子(其中一个甲烷碳)选中并命名为CA,这里的X指的是输入a 2返回后的组号。自己视情况而定,也不一定2就是你的碳原子编号。同理选中第二个甲烷碳并命名为CB

参数文件

mdp中添加一个关于质心拉伸代码部分。此代码是我们将甲烷保持一定距离的方法。可能有几种不同的设置方法(比如用拉力),但是对于这个系统,我们将手动指定两个甲烷之间的距离。点击下载获取所需的mdp文件(适用于gromacs2019
参数文件的拉伸部分做以下简单说明:

参数 解释
pull yes 使用拉伸代码
pull-ngroups 2 定义两个拉伸组
pull-group1-name CA 第一个组名
pull-group2-name CB 第二个组名
pull-ncoords 1 使用一个反应坐标
pull-coord1-geometry distance 沿着连接两个组的向量
pull-coord1-type umbrella 对这个坐标使用伞势
pull-coord1-groups 1 2 有两个组定义的拉伸坐标(定义如下)。实际上你可以有更多的拉伸坐标,使得穿过不同的分子,但在这里不适用
pull-coord1-k 5000.0 伞势中所用的力常数(kJ/(mol nm)
pull-coord1-init WINDOW 这是我们想要的两个组之间的距离。我在这里放置了这个关键字,我将在每个窗口的bash脚本中替换它
pull-coord1-rate 0.0 我们不想让组在坐标上移动,所以这里是0
pull-coord1-start no 我们手动指定每个窗口的距离,所以我们不想在计算中添加质心距离

设置了100 ps NVT平衡、1 ns NPT平衡和5 ns成品模拟运行的参数文件。我们正在计划甲烷达到正确的距离时,伞势是应用在平衡过程中。对于其他一些系统,在为每个窗口生成初始构象时,您可能必须更有条理。

模拟

对于模拟,我们将使用bash脚本替换mdp文件中的WINDOW关键字,这与我们在自由能模拟中所做的非常相似。脚本如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
#!/bin/bash
set -e
for ((i=0;i<27;i++)); do

x=$(echo "0.05*$(($i+1))" | bc);

sed 's/WINDOW/'$x'/g' mdp/min.mdp > grompp.mdp
gmx grompp -o min.$i -pp min.$i -po min.$i -n index.ndx
gmx mdrun -deffnm min.$i -pf pullf-min.$i -px pullx-min.$i

sed 's/WINDOW/'$x'/g' mdp/min2.mdp > grompp.mdp
gmx grompp -o min2.$i -c min.$i -t min.$i -pp min2.$i -po min2.$i -maxwarn 1 -n index.ndx
gmx mdrun -deffnm min2.$i -pf pullf-min2.$i -px pullx-min2.$i

sed 's/WINDOW/'$x'/g' mdp/eql.mdp > grompp.mdp
gmx grompp -o eql.$i -c min2.$i -t min2.$i -pp eql.$i -po eql.$i -n index.ndx
gmx mdrun -deffnm eql.$i -pf pullf-eql.$i -px pullx-eql.$i

sed 's/WINDOW/'$x'/g' mdp/eql2.mdp > grompp.mdp
gmx grompp -o eql2.$i -c eql.$i -t eql.$i -pp eql2.$i -po eql2.$i -n index.ndx
gmx mdrun -deffnm eql2.$i -pf pullf-eql2.$i -px pullx-eql2.$i

sed 's/WINDOW/'$x'/g' mdp/prd.mdp > grompp.mdp
gmx grompp -o prd.$i -c eql2.$i -t eql2.$i -pp prd.$i -po prd.$i -n index.ndx
gmx mdrun -deffnm prd.$i -pf pullf-prd.$i -px pullx-prd.$i
done

我们模拟了26个窗口,距离从0.051.3 nm左右。注意,我为每一步的拉力和拉力距离添加了-pf-px标志。这是因为使用-deffnm GROMACS时,将尝试将这两个文件都写入同一个文件(其实在gromacs2019中不设置也会产生)。我还用-n指定了索引文件,因为gmx grompp需要获得我们用pull参数指定的组。注意,使用了一个小技巧bc命令,以便在bash中对浮点数进行数学运算。

分析

我们要用gmx wham来得到PMF。该程序输入一个包含.tpr文件列表的文件和另一个包含.xvg文件列表的文件作为参数。
通过如下方式创建这两个文件:

1
2
ls prd.*.tpr > tpr.dat
ls pullf-prd.*.xvg > pullf.dat

然后运行gmx wham:

1
gmx wham -it tpr.dat -if pullf.dat

运行gmx wham之后,您将在名为profile.xvg的文件中获得平均力势。如果你马上画出来,它应该是这样的:

我们期望相互作用在更长的距离上趋于零。因为我们使用了一个三维偏置势,但是,我们需要一个校正。想象一个甲烷作为参考点。另一种甲烷可以在距离为r的点周围采样,覆盖半径为r的球体的表面。这给我们的采样增加了额外的构型空间,减小了熵。这个额外的熵对PMF的贡献需要去掉。回想一下等温等压系综中的吉布斯自由能是-kTln(W)其中W是配分函数。甲烷在球体表面跳舞的情况下,W与球体的表面积成正比。因此,需要添加2kTln(r)的校正。此外,我们还需要将图向上移动,使其尾部为零。我发现加大约77个对我的特定系统有效,但是您的系统可能有所不同。要在gnuplot中绘图,请在gnuplot终端中执行以下操作:

1
plot 'profile.xvg' u 1:($2+2*8.314e-3*298.15*log($1)+77) w l

你的PMF现在应该是这样的:

与教程3中的PMF相比,我们可以看到它们几乎是相同的:

一个不同之处在于,直接法永远不会像窗口采样那样接近。两个甲烷不会自然地想要靠近彼此,这就是为什么我们要增加伞势来保持它们在那里。
另一个输出是histo.xvg。这有助于确定窗口之间是否有足够的重叠。下面是本次模拟的每个直方图:

很明显,我们的窗口有足够的重叠。如果没有,我们可能不得不选择一个更小的窗口大小,或者选择一些缺少的特定点来模拟。


重复

自己按照上面的教程重复了一下,我这里为了节约时间,每个窗口的成品只进行了1 ns,原文为5 nm。基本分析结果如下,并教大家如何在xmgrace/qtgrace软件中绘图和转换。

  • histo.xvg中查看窗口重合度。读入数据到qtgrace中,具体操作可以参考Qtgrace绘图程序学习小记一文。对于多列数据记得选择NXY或者命令行指定-nxy。适当调整坐标轴可以得到下图: 可以看出窗口重叠程度良好,但是由于模拟时间比较短,相比原文曲线显得比较粗糙。
  • profile.xvg文件中的PMF。同理读入PMF数据,调整一下绘制下图:
  • 如原文所说,需要校正一下曲线。
    在上图的基础上,点击Data -> Transformations -> Evaluate expression,选中数据源Set,在下方输入:S1.Y = S0.Y+2*8.314E-3*298.15*ln(S0.X)+68.5并点击应用即可出现校正后的曲线: 表达式中的S0.Y表示源曲线Y轴数据,S1.Y为新建曲线。经过调整坐标轴,达到最终的PMF曲线如下: 结果与原文相符,重复基本成功。

    总结

    在本教程中,我们使用GROMACS的质心拉伸代码对水中的两种甲烷进行窗口采样。在此基础上,利用gmx wham算法提取了平均力势。