平均力势----以环糊精-苯酚体系为例

以环糊精-苯酚体系为例,得到两者之间的平均力势(PMF),即自由能势能面。

理论知识

网上有比较多的讲解关于自由能计算的文章,也有一些例子,见参考资料[1234],过多的东西也就不在此处讲了,本文主要讲解计算步骤。计算自由能的方法有很多,如下图:

本例采用伞形采样

实例

本文以计算环糊精和苯酚之间的结合自由能为例,阐明如何计算得到PMF。

材料及软件准备

  • 环糊精和苯酚分子结构文件从ChemSpider下载
  • 通过avogadro或者GaussView等软件给环糊精分子加氢,并检查结构是否正确
  • acpype在线或者单机均可。建议单机,因为对于较大的分子(环糊精)速度快(2 h 29 min vs 36 min
  • 可视化软件vmd和模拟软件GROMACS

拓扑获取

经过加氢和结构检查,那么我们首先需要得到单个分子的top文件,方法有很多,由于此处是想用amber力场,因此使用acpype去得到拓扑文件最为方便快捷。具体使用请看博文:分子动力学模拟三聚氰胺溶液性质中的部分。

体系构建

  • 方法也有很多,首先将环糊精通过gmx editconf构建一个适当的盒子,在这里为了减少计算量,所以盒子比较小,为2 x 12 x 2 $nm^3$。此处是以y轴为拉伸方向,当然你可以通过gmx editconf中的-rotate选项旋转分子,使用常规的z轴为拉伸方向。
  • 复制苯酚gro文件分子坐标,粘贴到上一步获得的gro文件后面,盒子参数之前。然后vmd可视化,按快捷键6,选取苯酚分子将其拖动到环糊精上方,然后按住shift键鼠标拖动调整苯酚分子的方向,最终模型如下图: 然后通过gmx solvate命令添加spc216.gro水,这里我采用tip3p水模型。

拓扑文件修改

因为top文件需要和gro文件对应,因此需要进行修改,即需要组合前面获得的两个分子的单独itp文件,注意要把它们两个的[ atomtypes ]部分剪切到前面来,并合并相同原子类型。最终见下图:

此处的posre_cyclodextrin.itpposre_phenol.itp是位置限制,通过gmx genrestr -f xxx.gro -o posre_xxx.itp -fc 1000 1000 1000获得。

模拟

前面几步分别为:

  • EM能量极小化
  • NVT预平衡
  • NPT预平衡
  • 成品模拟
    这些都比较常规,我就不说了。。

构象产生

这一步是通过拉伸力沿着反应坐标产生一系列的构象。
参数文件pull.mdp前面部分与NPT部分差别不大,但是在这里限制了环糊精的位置。主要差别在于后面的拉伸代码部分,用于产生一系列的构象。拉伸部分的mdp参数(5.1以上版本)如下:

1
2
3
4
5
6
7
8
9
10
11
12
pull                    = yes
pull_ncoords = 1 ; only one reaction coordinate
pull_ngroups = 2 ; two groups defining one reaction coordinate
pull_group1_name = UNL ; Cyclodextrin
pull_group2_name = MOL ; pull group, phenol
pull_coord1_type = umbrella ; harmonic potential
pull_coord1_geometry = distance ; simple distance increase
pull_coord1_dim = N Y N ; pull along Y
pull_coord1_groups = 1 2 ; groups 1 (Cyclodextrin) and 2 (phenol) define the reaction coordinate
pull_coord1_start = yes ; define initial COM distance > 0
pull_coord1_rate = 0.0025 ; 2.5 nm per ns
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2

很好理解,沿着Y轴拉伸分子苯酚逐渐远离环糊精,以此产生一系列构象。

构象选择

沿着反应坐标Y,会产生许多的构象,我们需要通过每个构象的距离选取部分典型构象来研究。

  • 计算每一个构象的距离,比如我产生了500个构象,首先需要计算每一个构象中环糊精和苯酚分子的质心距离,可以通过下面的bash脚本产生:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    echo 0 | gmx trjconv -s pull.tpr -f pull.trr -o conf.gro -sep

    for (( i=0; i<501; i++ ))
    do
    gmx distance -s pull.tpr -f conf${i}.gro -n index.ndx -select 'com of group "UNL" plus com of group "MOL"' -oall dist${i}.xvg
    done

    touch summary_distances.dat
    for (( i=0; i<501; i++ ))
    do
    d=`tail -n 1 dist${i}.xvg | awk '{print $2}'`
    echo "${i} ${d}" >> summary_distances.dat
    rm dist${i}.xvg
    done

    exit;
  • 选取构象
    这里还是按照那个蛋白的伞形采样方式,选取距离差为0.2 nm左右的构象,脚本去这里下载,使用方式为:

    1
    python setupUmbrella.py summary_distances.dat 0.2 run-umbrella.sh &> caught-output.txt

这样就得到了所有在0.2 nm左右的构象,然后对这些构象通过伞形采样模拟即可。亦或者试试我这个umbrella.bsh脚本去产生运行文件:
使用方式为:

1
bash umbrella.bsh 0.2 summary_distances.dat

  • 伞形采样

    • 上一步的构象,可以用来模拟了,首先在NPT下进行一个短暂的平衡,注意此时的NPT的参数文件mdp中的拉伸代码部分需要将拉伸速率设置为0,即pull_coord1_rate = 0.0,也包括continuation = no, gen_vel = yes, gen_temp = 298refcoord_scaling = com
    • 正式伞形采样模拟。进行长时间的模拟,比如10 ns以上,可以通过设置多个时间梯度分别完成自由能曲线图来比较是否收敛。注意此时的mdp中的参数设置包括continuation = yes, gen_vel = norefcoord_scaling = com,即利用上一步的npt.gro中的速度。
      最好写个脚本去完成以上步骤。
  • 提取PMF

    • 一种常用的PMF提取方法是加权直方图分析法(WHAM)。分别准备tpr-files.datpullf-files.dat文本文件,内容为上一步参数的所有.tpr文件名(包括后缀)和*_pullf.xvg文件名(包括后缀)。例如:
    1
    2
    3
    4
    5
    umbrella0.tpr
    umbrella11.tpr
    umbrella38.tpr
    ...
    umbrella494.tpr
    1
    2
    3
    4
    5
    umbrella0_pullf.xvg
    umbrella11_pullf.xvg
    umbrella38_pullf.xvg
    ...
    umbrella494_pullf.xvg
    • WHAM计算PMF。通过命令:
    1
    gmx wham -it tpr-files.dat -if pullf-files.dat -o -hist

就可以得到直方图和PMF自由能曲线图。

由于模拟时间较短,并且从直方图中也可以看出在1.5 nm附近重叠并不好,需要在这里多取几个构象计算,并延长模拟时间直到收敛。我这里就不做了,懒得算。

补充

如果你对顺反异构转变的自由能感兴趣,你可以参看AMBER高级教程A17。如果你想用GROMACS进行该过程的模拟,也是可以的。首先需要定义拉伸部分,给个参考如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
; Pull code
pull = yes
pull_ncoords = 1 ; only one reaction coordinate
pull_ngroups = 6 ; six groups defining one reaction coordinate
pull_group1_name = O
pull_group2_name = C
pull_group3_name = C
pull_group4_name = N
pull_group5_name = N
pull_group6_name = H
pull_coord1_type = umbrella ; harmonic potential
pull_coord1_geometry = dihedral ; dihedral
pull_coord1_groups = 1 2 3 4 5 6
pull_coord1_start = yes ; define initial COM distance > 0
pull_coord1_rate = 0.2 ; deg/ps
pull_coord1_k = 1000 ; kJ mol-1 rad-2

测试后发现产生构象和伞形采样过程是没有问题的,但是当利用gmx wham -it tpr-files.dat -if pullf-files.dat -cycl去产生PMF曲线时出现-nan值,查阅资料发现并没有太多关于这方面的问题,仅仅在1,2这些地方提及过,但好像并没有实质性的解决,虽然GROMACS 2019最新版本的wham代码改了一部分。感觉还是单位转换有问题,比如这里的力常数1000单位是kJ mol-1 rad-2,而gmx wham使用的单位是kJ mol-1 deg-2。因此最简单的方式是改一下\src\gromacs\gmxana\gmx_wham源码中的单位,然后重新编译,个人认为需要将力常数乘以系数$(pi/180)^2$,对应于源代码中的:

测试发现这样产生的PMF曲线是正常的,与Amber教程的结果保持一致。

注意:这么更改以后的gmx wham只能用于计算角度相关的PMF,不能用于距离拉伸的情况。
结果对比如下,,有点奇怪。

总结

利用伞形采样来计算自由能曲线是一种非常常见的方法,许多论文中的PMF曲线图都是这么来计算的。一般需要注意的就是重叠模拟时间要足够好才行。

参考资料

[1] 杂谈自由能计算,PMF,伞形抽样,WHAM
[2] Umbrella Sampling
[3] Umbrella sampling DOI:10.1002/wcms.66
[4] GROMACS Tutorial 5 - Methane-methane PMF from window sampling