教程翻译----使用GROMACS计算乙醇分子溶剂化自由能

翻译并重复了一个有关GROMACS计算乙醇溶剂化自由能的实例。

背景

在这个教程中,我们将计算一个小分子(乙醇)溶剂化的自由能。这种计算可以单独进行,也可以作为结合自由能计算的一部分。这样的计算很重要,因为自由能是热系统中最重要的静态量:它的符号决定了一个分子是否可溶,或者是否能够与另一个分子结合。
我们将从关于如何计算自由能,以及溶剂化的自由能与结合自由能之间的关系的一些背景知识开始本教程学习。然后我们将集中讨论在GROMACS中进行这种计算的实用性。在本教程中,您应该使用GROMACS 5.1(或更高版本)。

计算结合自由能

计算自由能通常只能使用很小的步骤和从一个状态到另一个状态的完整的路径。例如,为了计算配体与蛋白质的结合自由能,我们最终需要比较配体与蛋白质结合的情况,以及配体与蛋白质分别在溶液中的情况:

这可以直接计算出来,例如将配体从蛋白质中拉出,并对平均力势(平均力,并对其积分)进行积分。然而,力有非常大的波动,这比我们将在本教程中使用的自由能量扰动(free energy perturbation,FEP)方法(如Bennett Acceptance Ratio(BAR))昂贵得多。
记住,AB两种状态之间的自由能差决定了它们的相对概率$p_A$和$p_B$。
$$\frac{p_A}{p_B}=exp \frac{F_B-F_A}{k_BT}$$
这里$k_B$是玻尔兹曼常数,将热能与温度联系起来($1.38*10^{-23} J/K$),T是温度。原则上,我们可以通过等待足够长的时间,并测量系统处于这种状态的频率来计算自由能差。然而,自由能的差异通常是几十kJ/mol

例如,298 K时乙醇溶剂化的自由能为-20.1 kJ/mol,相当于$-8.1 k_BT$,相对概率为$3*10^{-4}$。
我们需要等待很长一段时间才能自发地发生这种转变,甚至需要更长的时间才能得到关于这种转变的良好统计数据。
由于这个概率问题,自由能方法依赖于一个基本概念:强迫系统到它不想去的地方,然后用它不想去的程度来衡量。在FEP方法中,我们通耦合感兴趣的分子与系统其余部分相互作用强度来强迫系统到变量$\lambda$:
$$E_{total}=E_{ligand-ligand}+E_{rest-rest}+\lambda E_{ligand-rest}$$
我们慢慢的将$\lambda$从1变到0。这意味着我们可以有效地关闭一个分子,假装它在真空中($\lambda = 0$):我们强迫系统到它不想去的地方(无论是在溶剂化状态还是在真空状态,这取决于自由能差的符号)。然后我们会用BAR方法来计算它不想在那的程度。
耦合和解耦合的方式帮助我们计算结合自由能,因为我们现在可以创建一个两步路径:

我们首先将配体从溶剂中分离出来,然后在蛋白质存在下将配体重新溶剂化。因此结合自由能是这样的:
$$\Delta G_{binding}=\Delta G_1+\Delta G_2$$
该模拟分为两部分:一是计算解溶剂化自由能,计算分子与蛋白质耦合进入系统的自由能;二是模拟耦合配体从$\lambda = 0$(它不与系统相互作用)到$\lambda =1$(它与蛋白结合)。第一个模拟是溶剂化自由能的逆过程。这是我们将在本教程中重点关注的,计算性能也是一个原因:由于不涉及蛋白质,模拟盒的大小可以很小,模拟速度也会很快。

溶剂化自由能

为了计算溶剂化自由能,我们计算上图所示的$-\Delta G_1$,或者等价于下图的$\Delta G_{solv}$

我们把分子和耦合到变量$\lambda$(上面的等式),然后通过GROMACS内置的BAR计算。
BAR方法依赖于在状态$\lambda_A$和$\lambda_B$成对模拟的输出,如果$\lambda_A$和$\lambda_B$足够接近,自由能差能够直接计算出来(见Bennett’s原始文章:Bennett, J. Comp. Phys, (1976) vol. 22 p. 245),通过计算从$\lambda_A$到$\lambda_B$的Monte Carlo acceptance rates,反之亦然。这里的“足够接近”一词意味着在这两种状态之间的切换应该在两个方向上都是可能的:在两个端点都应该允许一些相同的构象(例如它们应该共享相空间的某些部分)。
从$\lambda_A$到$\lambda_B$过程中最明显的点就是$\lambda_A=0$和$\lambda_B=1$。然而,这些端点通常很少有共同的状态:它们共享很少的相空间。正因为如此,自由能永远不会收敛到一个可用的值。这就是为什么我们要把问题分成:

使用多个$\lambda$是需要的。我们将有效地“缓慢”打开(或关闭)配体与溶剂之间的相互作用。这意味着我们需要运行尽可能多的模拟$\lambda$点,我们需要告诉每一个模拟邻近$\lambda$点,我们将结合许多模拟结果对结果进行后处理(我们使用7个$\lambda$点:0, 0.2, 0.4, 0.6, 0.8, 0.9和1)。例如我们运行一个在$\lambda = 0.4$的模拟,这个模拟将计算出这个$\lambda$和附近$\lambda = 0.2$和$\lambda = 0.6$之间的能量差。
我们要走一条捷径:同时关闭静电(Coulomb)相互作用和范德华(Lennard-Jones)相互作用。为了得到高质量的结果,这些阶段通常是分开的,但是为了方便起见,我们将同时进行这两个阶段。Gromacs使用“软核”相互作用来确保正常的相互作用(Lennard-Jones和Coulomb)关闭,永远不会有两个点电荷相互叠加:这是通过在中间$\lambda$点开启一种有效排斥粒子的相互作用来实现的(这样它就和自由能差抵消了)。

准备体系

我们将从一个拓扑文件开始,能够在http://www.gromacs.org/Documentation/Tutorials/Free_energy_of_solvation_tutorial下载,得到gromacs-freeenergy-tutorial.tar.gz。或者直接:

1
wget http://www.gromacs.org/@api/deki/files/261/=gromacs-free-energy-tutorial.tgz

然后解压文件夹即可:

1
tar xzvf gromacs-free-energy-tutorial.tar.gz

然后可以看到topol.top文件和一个非常基本的坐标文件,名为ethanol.gro。这个拓扑使用的是OPLS力场,定义了一个乙醇分子,包括定义了SPC/E水分子。
问题:
查看拓扑文件topol.top。对于乙醇分子的定义,你能找到有哪些原子,它们是如何连接的吗?我们通过借用苏氨酸侧链得到了乙醇的参数。
我们首先准备模拟:原始构象文件有一个与之关联的模拟盒子(您可以通过查看文件ethanol.gro看到这一点)。我们这样做:

1
gmx editconf -f ethanol.gro -o box.gro -bt dodecahedron -d 1

设置模拟盒子。在这种情况下,它将使模拟盒子做成一个斜方十二面体,溶质(乙醇分子)与盒子边距最小为1 nm。这个盒子是一个菱形的十二面体,因为它提供了一个比矩形盒子更有效的周期镜像:我们可以用更少的水来处理相同距离的乙醇分子周期成像。请参阅Gromacs手册,了解这个盒子的形状及其周期性成像的排列方式。
接下来,我们把系统溶解在水中:

1
gmx solvate -cp box.gro -cs -o solvated.gro -p topol.top

这将生成一个系统,其中包含~300个水分子,这些水分子来自-cs选项的默认文件名:平衡了的水分子。
为了使结构适合于模拟,我们将首先最小化它的能量,两次:一次用柔性键,一次用约束键。为了最小化,我们将使用以下设置(请参阅包含的文件em.mdp):

1
2
3
4
5
6
------em.mdp------
integrator = steep
nsteps = 500
coulombtype = pme
vdw-type = pme
------------

通过将输入文件预处理为一个运行文件来运行平衡:

1
gmx grompp -f em.mdp -c solvated.gro -o em.tpr

它生成运行文件em.tpr。使用以下命令运行此文件:

1
gmx mdrun -v -deffnm em

由于这是一个非常小的系统,您可能需要处理MPI级别的数量(用于域分解)或OpenMP线程的数量。例如,您可以添加-nt 4来将线程(处理器核心)的数量限制为4个,因为具有大约1000个原子的模拟系统太小,无法支持高度并行运行。您设法在多少个内核上运行它?如果你发现这是一个麻烦,并且有多余的资源,你也可以先创建一个更大的水盒子。

预平衡

温度耦合:我们试着计算吉布斯自由能的差,为了达到这个目的,当乙醇分子解耦时,系统必须保持温度和压强不变。全局平衡(即在我们施加几个不同的值之前所做的平衡)是用equil.mdp来完成的:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
------equil.mdp------
integrator = md
nsteps = 20000
dt = 0.002
nstenergy = 100
rlist = 1.0
nstlist = 10
vdw-type = pme
rvdw = 1.0
coulombtype = pme
rcoulomb = 1.0
fourierspacing = 0.12
constraints = all-bonds
tcoupl = v-rescale
tc-grps = system
tau-t = 0.2
ref-t = 300
pcoupl = berendsen
ref-p = 1
compressibility = 4.5e-5
tau-p = 5
gen-vel = yes
gen-temp = 300
------------

我们将使用v-rescale恒温器和Berendsen恒压器。我们运行:

1
2
gmx grompp -f equil.mdp -c em.gro -o equil.tpr
gmx mdrun -deffnm equil -v

与以前一样,您可能需要限制使用的核心数量。在这之后,我们应该在一分钟内准备好乙醇在水中的平衡结构。输出构象的名称是equil.gro
问题:
你应该检查一下系统是否平衡了。你怎么能这么做?

创建$\lambda$点

平衡后,我们准备把系统分割成不同$\lambda$点。为了简化这一过程,我们准备了一个名为mklambdas.sh的脚本,可以在前面提到的压缩文件中找到。对于最后一次运行,我们将使用run设置run.mdp

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
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
------run.mdp------
; we'll use the sd integrator with 100000 time steps (200ps)
integrator = sd
nsteps = 100000
dt = 0.002
nstenergy = 1000
nstlog = 5000
; cut-offs at 1.0nm
rlist = 1.0
dispcorr = EnerPres
vdw-type = pme
rvdw = 1.0
; Coulomb interactions
coulombtype = pme
rcoulomb = 1.0
fourierspacing = 0.12
; Constraints
constraints = all-bonds
; set temperature to 300K
tcoupl = v-rescale
tc-grps = system
tau-t = 0.2
ref-t = 300
; set pressure to 1 bar with a thermostat that gives a correct
; thermodynamic ensemble
pcoupl = parrinello-rahman
ref-p = 1
compressibility = 4.5e-5
tau-p = 5
; and set the free energy parameters
free-energy = yes
couple-moltype = ethanol
; these 'soft-core' parameters make sure we never get overlapping
; charges as lambda goes to 0
sc-power = 1
sc-sigma = 0.3
sc-alpha = 1.0
; we still want the molecule to interact with itself at lambda=0
couple-intramol = no
couple-lambda1 = vdwq
couple-lambda0 = none
init-lambda-state = $LAMBDA$
; These are the lambda states at which we simulate
; for separate LJ and Coulomb decoupling, use
fep-lambdas = 0.0 0.2 0.4 0.6 0.8 0.9 1.0
------------

这里有一个关于自由能设置的扩充部分。有些值像$LAMBDA$:这些值将被脚本mklambdas.sh替换。自由能设置如下:把乙醇分子耦合到变量$\lambda$。$\lambda = 0$意味着分子解耦合,$\lambda = 1$表示分子是完全耦合的(vdwq表示Lennard-Jones + Coulomb)。sc-power, sc-sigmasc-alpha设置控制“软核”相互作用,防止系统在解耦时产生重叠粒子。
仍然需要设置实际init-lambda的$\lambda$值和foreign-lambda值:该字段决定的其他$\lambda$值模拟应该计算能量差异。就我们的目的而言,我们只会计算所有其他$\lambda$值的能量差异和对所有模拟保持不变(这是对性能的影响可以忽略不计)。
脚本mklambdas.sh将为每个lambda点创建一个目录,并用构象、拓扑和模拟设置填充它,同时替换mdp中的\$LAMBDAS\$和\$ALL_LAMBDAS\$。我们可以这样运行:

1
bash mklambdas.sh run.mdp topol.top equil.gro

它应该生成许多目录:

1
2
3
4
5
6
7
lambda_00
lambda_01
lambda_02
lambda_03
lambda_04
lambda_05
lambda_06

每一个包含了

1
2
3
conf.gro
grompp.mdp
topol.top

验证替换操作是否正确:

1
grep init-lambda lambda_*/grompp.mdp

应该会显示:

1
2
lambda_00/grompp.mdp:init-lambda-state = 0
lambda_01/grompp.mdp:init-lambda-state = 1

等等,每一个$\lambda$点都需要如下运行:

1
2
3
4
cd lambda_00
gmx grompp
cd ../lambda_01
gmx grompp

检查输出是否成功,只是为了确定。现在我们可以开始运行了。
总运行时间将会在计算机为每个$\lambda$点只要几分钟。这意味着我们可以按顺序运行这些作业,但之后我们将不得不等待一段时间,这将浪费并行化的大好机会。
 相反,我们将尝试并行运行它们,并假设我们有某种批处理系统,可以将作业提交给该系统。因为这个系统只有1000个粒子,所以缩放是很困难的。通常,现代计算集群的每个节点可能有32-64个或更多内核。因此,我们将欺骗队列系统运行多个作业。
如果您有一个适当的MPI库(即,而不是内置的thread-MPI),您可以使用-multidir选项来gmx mdrun并指定所有不同的目录。多玩一会儿,看看会发生什么!您可以使用多个MPI级别(每个模拟一个)和多个核心来混合/匹配吗?
 另一种选择(特别是如果没有任何GROMACS的MPI版本)是在一个批处理脚本中启动多个作业。要在批处理系统设置中使用的确切设置(通常#PBS#SBATCH取决于正在使用的批处理系统)字段严重依赖于正在运行的系统,因此这里的示例可能无法工作(为您的系统修改它)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#!/bin/sh
#SBATCH -J fe
#SBATCH -e run1.stderr
#SBATCH -o run1.stdout
#SBATCH --mem-per-cpu=1000
#SBATCH -t 00:10:00
# One node of 12 cores per job:
#SBATCH -N 1
#SBATCH -n 12
module load gromacs # add suitable modules such as cuda, fftw, etc.
( cd lambda_0; gmx mdrun -nt 4 >& run.log ) &
( cd lambda_01; gmx mdrun -nt 4 >& run.log ) &
( cd lambda_02; gmx mdrun -nt 4 >& run.log ) &
wait # wait for all background tasks to finish

它将运行三个4线程版本的mdrun(设置为-nt 4),每个版本都在自己的目录中,最长运行10分钟。让尽可能多的这些脚本的所有不同的$\lambda$-values(例如3),并将它们提交到队列中。

后处理:提取自由能

模拟完成后,我们可以从输出数据中提取出全自由能差。
检查目录lambda_00lambda_06中名为dhdl.xvg的文件。这些包含了计算自由能差的能量差。使用Gromacs BAR工具gmx bar将它们组合成自由能:

1
gmx bar -b 100 -f lambda_*/dhdl.xvg -o -oh -oi

这是纯粹的魔法。作为mdrun中的自由能代码的一部分,我们已经计算了抵消焓变到邻近$\lambda$点,所以这在dhdl.xvg中是可用的(连同点本身的所有信息,它是什么模拟,等等)。不需要重新运行任何模拟或存储整个轨迹。其次,gmx bar命令将完成BAR自由能所需的所有复杂处理,并只给出结果。其中-b 100表示前100 ps应该被忽略:它们作为另一种平衡,这一次是在模拟的条件下。你应该得到一个近似的自由能差-20.6 +/- 2.4 kJ/mol(如果在不同的硬件上运行,情况可能有所不同:这个结果来自标准的x86_64集群)。这应该与实验值-20.9 kJ/mol进行比较。

问题:
随着标准误差估计值的减小,较长时间的运行会稍微改变自由能值(试试)。您得到的确切值将取决于许多设置,比如是否使用LJ-PME。为什么有时实验结果与模拟结果之间会有显著的差异(即大于估计误差)?如何改进这一点呢?
问题:
查看各个$\lambda$点的误差:单个点之间的差异很多。这对整体计算的效率意味着什么?如何改进呢?

下一步该去哪儿

 在计算了溶剂化自由能后,我们解出了早期方程结合自由能的第一部分。第二部分是将分子耦合到(或脱离)与蛋白质结合的状态。这就引入了一个额外的复杂性:我们最终会遇到弱耦合配体在我们的系统中游荡的情况:

 这很糟糕,因为这是一种不可逆的情况:突然间很少有状态能从弱耦合映射到更强耦合的分子,这将大大降低自由能计算的准确性。
 这种情况可以通过迫使配体停留在相对于蛋白质的特定位置来补救。这可以通过Gromacs的pull code来实现,该代码允许指定任意的力或约束作用于任意一组原子的质心上,这些力或约束作用于任意一组原子上。使用伞形的拉式,我们可以指定我们想要这个指定位置的二次电位,即使配体完全解耦,也会迫使配体停留在其原始位置。
 一种确定力的中心位置的方法是选择蛋白质中靠近配体的一组原子,并在完全配体耦合的情况下进行模拟,在这种情况下拉伸代码是启用的,但力为零。然后,拉伸代码将频繁地输出配体的坐标,从中可以计算出平均位置和预期偏差。然后,这可以作为生产运行期间拉代码的力中心位置和拉伸代码的力常数的参考点。
 一旦自由能被计算出来,就必须注意纠正我们已经捕获了分子的事实。这很容易进行分析。

可选问题:
给定配体质心位置的测量标准偏差,我们如何为拉伸代码选择力常数?
可选问题:
我们如何正确使用拉伸代码:对分子施加二次势能对自由能的贡献是多少?

教程改进之处

  • 设置的lambda点太少,误差太大。
  • 单机批量运行脚本能够更方便,比如:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    #!/bin/bash
    for i in {0..6}
    do
    folder=`printf "lambda_%02d" $i`
    cd $folder
    gmx grompp
    gmx mdrun -v
    cd ..
    done
  • 同时关闭静电(Coulomb)和范德华(VDW)相互作用的方法欠佳。

  • 模拟时间太短。

重复教程

结果与上面教程保持一致。增加$\lambda$点会使得误差变小。比如求得的溶剂化自由能为-18.05 +/- 1.68 kJ/mol

其他参考资料

[1] Guidelines for the analysis of free energy calculations
[2] GROMACS Tutorial 3 - Several methanes in water
[3] Free Energy Calculations: Methane in Water
[4] Numerical modelling of free energy for methanol and water mixtures for biodiesel production