教程翻译----SPCE水模型的PMF

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

原文链接:https://github.com/chen3262/Chem7590-waterPulling

SPCE水模型的PMF

沿所选坐标的自由能面称为平均力势(PMF)。通常,PMF模拟与平衡牵引方法(umbrella sampling)或非平衡牵引方法(Jarzynski equility)结合使用。两种方法都允许在高能量屏障存在时进行增强采样。PMF也可以由沿反应坐标的径向分布函数g(r)构成。然而若存在高能势垒,g(r)可能无法提供准确PMF。在本练习中,您将分析由两种牵引方法和g(r)产生的PMF。具体地说是水的氧-氧之间的自由能面。

起始文件

本练习使用的是Owens上的GROMACS软件包。您可以从/users/PAS1314/osu9080/Chem7590/waterPulling复制必要的文件。您将使用以下批处理脚本运行所有作业:

1
2
3
4
Batch NTP.sh
Batch GenConfig.sh
Batch Umbrella.sh
Batch Jarzynski.sh

在使用脚本之前,必须替换input文件夹的路径(第24行)。

SPC/E水的径向分布函数

相关命令:

1
gmx rdf -f md.trr -s md.tpr -n index.ndx -o rdf_OW_OW.xvg

然后根据下面PMFg(r)关系得到PMF曲线。
$$W(r)=-k_{B}Tln(g(r))$$

略,详见原文。

Umbrella sampling方法

  • 首先产生一系列构象。拉伸部分的mdp文件内容如下:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    pull			= yes
    pull_coord1_type = constraint ;constant-force
    pull_coord1_geometry = distance ;direction-periodic ;distance
    pull_coord1_dim = Y Y Y
    pull-group1-name = target
    pull-group2-name = ref
    pull-ngroups = 2
    pull-ncoords = 1
    pull-coord1-groups = 1 2
    pull_coord1_rate = 0.005 ; nm per ps
    pull_coord1_start = yes
    pull-nstxout = 100
    pull-nstfout = 100
  • 然后分别运行单个窗口模拟,这里是根据拉伸时间来确定要单独模拟的窗口构象,并进行单独的简谐力限制下的模拟。计算脚本和相应的mdp中的拉伸部分代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
#!/bin/bash
# get configure
for i in 0 40 60 80 120 140
do
echo 0| gmx trjconv -f md.trr -s md.tpr -b $i -e $i -o Conf_${i}.gro
done

# do each simulation
for i in 0 40 60 80 120 140
do
gmx grompp -f waterUmbrella.mdp -c Conf_${i}.gro -n water.ndx -p water.top -o md_${i}.tpr
gmx mdrun -deffnm md_${i} -v
done
1
2
3
4
5
6
7
8
9
10
11
12
13
14
pull			= yes
pull_coord1_type = umbrella
pull_coord1_geometry = distance
pull_coord1_dim = Y Y Y
pull-group1-name = target
pull-group2-name = ref
pull-ngroups = 2
pull-ncoords = 1
pull-coord1-groups = 1 2
pull-coord1-k = 1000 ; kJ mol^-1 nm^-2
pull_coord1_rate = 0.0 ; nm per ps
pull_coord1_start = yes
pull-nstxout = 1
pull-nstfout = 1
  • 最终得到校正以后的PMF曲线如下图: 略,详见原文。
    注意:原文提及到的C++编译好的Nearest.exe,在Windows平台上没法运行,目的是为了找到目标氧原子周围最近距离,约为0.25 nm的参考氧原子索引号。我曾想过用bash脚本直接去处理water.gro坐标得到距离,但是直接这么做需要考虑体系周期性问题,虽然我这里也提供了一个对应的bash脚本直接执行bash nearest.bsh -i water.gro -n 1108就可以。而这里我引入了一个tpr然后用gmx distance去计算距离,这样就不用了自己去考虑周期性问题了。具体脚本如下,如文中所言找到索引号1108的氧原子距离约0.25 nm的另一个氧原子索引号(274):
    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
    27
    28
    29
    30
    31
    #!/bin/bash
    # a script to find the nearest OW for target OW you given
    # 定义
    tpr=$1 # md.tpr
    gro=$2 # water.gro
    targetndx=$3 # 1108
    rm -rf distance.txt
    # 打印表头
    printf "%-8s%8s\n" "#Index" "#Distance(nm)" > distance.txt

    # .gro中氧原子编号从1开始,每隔3为氧原子,最后一个氧原子编号为1528
    for refndx in {1..1528..3}
    do
    {
    # 排除自身
    if [ $refndx != $targetndx ];then
    # 创建ndx文件并计算距离
    echo -e "a $targetndx $refndx\n q\n" |gmx make_ndx -f $gro -o OW_${refndx}.ndx > /dev/null 2>&1
    echo -e "a_${targetndx}_${refndx}" |gmx distance -s $tpr -f $gro -n OW_${refndx}.ndx -oall dist_$refndx.xvg > /dev/null 2>&1
    # 提取距离并打印
    dist=`tail -n 1 dist_${refndx}.xvg |awk '{print $2}'`
    #echo "Write $refndx into distance.txt "
    printf "%-8d%8.3f\n" $refndx $dist >> distance.txt
    rm -rf OW_${refndx}.ndx dist_${refndx}.xvg
    fi
    } &
    done

    wait
    # 距离按照大小排序
    sort -n -k2 distance.txt > distance_sort.txt

执行bash nearest.bsh md.tpr water.gro 1108以后会得到distance_sort.txt文件,内容为:

1
2
3
4
5
6
#Index  #Distance(nm)
274 0.251
907 0.268
1168 0.299
1291 0.316
...

得到参考原子编号274,与原文一致。

非平衡拉伸(Jarzynski)方法

另一种获得PMF的方法是进行非平衡拉伸,利用Jarzynski等式将自由能差与不可逆功联系起来:

选取5个目标氧原子,并确定其最近距离的参考氧原子编号,写进ndx。然后对每一个进行成品模拟,即5个拉伸模拟。得到一系列输出文件,然后采用有限差分近似方法进行评估:

可以通过如下命令结合上面的脚本得到五个模拟的结果:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# define folders for input and output
mkdir Output_Jarzynski

# move data and executable to tmp folder
mkdir workdir
cp md.tpr spce.itp water_tmp.ndx waterJarzynski.mdp water.gro water.top nearest.bsh workdir

# enter tmp folder
cd workdir

# do five simulation
for i in 1108 1198 1036 1516 1243 #loop over chosen OW atoms
do
bash nearest.bsh md.tpr water.gro $i
wait
ref=`awk 'NR==2 {print $1}' distance_sort.txt`
cat water_tmp.ndx | sed -e s/tttt/${i}/ | sed -e s/rrrr/${ref}/ > water_${i}.ndx
gmx grompp -f waterJarzynski.mdp -c water.gro -n water_${i}.ndx -p water.top -o OW_${i}.tpr > /dev/null 2>&1
gmx mdrun -v -deffnm OW_${i} -nt 2 > /dev/null 2>&1
echo "Index $i has been finished!"
done

# retrieve output
cp OW_* water_*.ndx ../Output_Jarzynski

然后根据产生的OW_*_pullx.xvgOW_*_pullf.xvg文件去得到拉力所做的功(work (kJ/mol))。这里可以采用一个bash小脚本去处理数据,使用方式例如:

1
bash xf_to_work.bsh OW_1108_pullx.xvg OW_1108_pullf.xvg 1

其中OW_1108_pullx.xvg OW_1108_pullf.xvg为模拟输出文件,分别为距离和力;1表示输出文件名为work1.dat。重复命令并得到五次拉伸情况下的距离与功的关系(work1.dat, ..., work5.dat)。绘图如下:

与原教程所给的图类似。
接下来将5次平均,得到下面的图:

最后对曲线按照如下公式进行校正,可参考前面写的一篇教程。这里我用的是S1.Y = S0.Y+2*8.314E-3*298.15*ln(S0.X)+10.5

校正以后的曲线如下图:

结论

比较不同方法得到的PMF曲线的重叠性如何?