四种计算自由能方法的原理示例教程|Jerkwin

原始文档 Calculating Free Energies
2017-05-24 07:18:53

【按】本教程以一个简单实例对计算自由能的四种方法进行了示例. 教程没有过度讲解原理, 而是从实际操作出发, 教你如何手动计算. 知晓了具体的计算操作过程对理解其他自由能计算程序的计算原理和计算结果很有帮助. 故翻译此教程供参考.

在本教程中, 我们将要计算Cu(I)和Cu(II)离子之间转移一个电子导致的自由能变化.

问题: 这个自由能差值有多大呢?

由于我们已经知道要计算的自由能差值的参考值, 所有我们可以验证模拟结果的合理性和精确性. 我们将逐步进行, 使用互不相同却又有联系的几种方法进行计算:

自由能微扰(Zwanzig)
接受比例(Bennet)
热力学积分(Kirkwood)
非平衡近似(Jarzynski)

体系设置

首先, 使用你熟悉的分子编辑软件创建两个铜离子, 将其保存为pdb文件. 然后, 修改文件使两个铜原子处于不同的组中(两个离子具有不同的组编号. 如果你不确定, 请查看pdb文件的格式说明). 此外, 将两个组都命名为CU1+, 其原子名称都是CU. 使用相同组名称的原因在于, 在我们将要使用的GROMOS力场中, 这两个离子所用的Lennard-Jones参数不同, 这会导致后面的处理变得比较复杂. 在这里我们不关心原子类型的精确性, 但对于更精确的工作, 你应该也考虑改变原子类型.

将两个离子置于4x4x4 nm^3立方盒子的中心. 假定你的pdb文件为2Cu.pdb, 可使用如下命令:

gmx editconf -f 2Cu.pdb -c -box 4 4 4 -o 2Cu_in_box.gro

接下来, 我们要使用pdb2gmx来自动生成体系的拓扑文件. 在此之前, 要检查下原子名称是否如我们所要求的那样. 使用熟悉的文本编辑器打开你创建的坐标文件2Cu_in_box.gro, 其内容类似下面:

two copper ions in a box
2
    1CU1+    CU    1   0.000   0.000   0.000
    2CU1+    CU    2   0.000   0.000   2.000
4.00000   4.00000   4.00000

如果你的文件如上面所示(否则请修改), 可以使用pdb2gmx命令

gmx pdb2gmx -f 2Cu_in_box.gro

运行命令后会提示你选择力场和水模型, 我们选择GROMOS96 54a7力场和spc水. 查看得到的文件, 特别是topol.top文件, 其中包含了体系的拓扑信息, 两个铜离子所带的电荷. 对于第一个铜离子带1+电荷而第二个铜离子带2+电荷的情况, 我们称其为状态A; 当第一个铜离子带2+电荷而第二个铜离子带1+电荷时, 我们称其为状态B. 本教程的目的是计算状态A和状态B之间的自由能差值. 试着理解拓扑文件中有关电荷的设置. 我们将通过更改这些电荷来创建状态A和B对应的拓扑文件.

由于我们要计算的是两个离子在水中进行电子转移导致的自由能差值, 我们使用GROMACS的solvate命令自动向盒子中添加水分子:

gmx solvate -cp conf.gro -cs spc216.gro -p -o solvated.gro

最后一个选项保证自动更新拓扑文件. 查看修改后的拓扑文件. 你也可以使用熟悉的查看程序查看溶剂化后的体系solvated.gro(或先使用editconf其转换为pdb文件再查看). 你可以看到加入的水看起来比较有序, 我们首先要进行能量最小化消除体系中过大的相互作用力. 你可以使用下面的mdp文件. 查看GROMACS手册中mdp的有关关键词, 学习并理解每个选项的含义.

integrator           = steep
nsteps               = 500
comm-mode            = Linear
nstcomm              = 1
comm-grps            =

emtol                = 100
emstep               = 0.01

nstxout              = 0
nstvout              = 0
nstfout              = 0

nstlog               = 500
nstenergy            = 500

nstlist              = 5

ns_type              = grid
pbc                  = xyz
rlist                = 0.9

coulombtype          = PME
rcoulomb-switch      = 0
rcoulomb             = 0.9
epsilon-r            = 1
vdw-type             = Cut-off
rvdw-switch          = 0
rvdw                 = 0.9

DispCorr             = EnerPres

fourierspacing       = 0.12
pme_order            = 4
ewald_rtol           = 1e-05
ewald_geometry       = 3d
epsilon_surface      = 0
optimize_fft         = no

Tcoupl               = berendsen
tc-grps              = System
tau_t                = 0.1
ref_t                = 300

Pcoupl               = berendsen
Pcoupltype           = isotropic
tau_p                = 0.5
compressibility      = 4.5e-5
ref_p                = 1.0

gen_vel              = yes
gen_temp             = 300
gen_seed             = 1993

将以上设置保存为em.mdp, 然后使用grompp创建运行mdrun所需的输入文件:

gmx grompp -f em.mdp -c solvated.gro -p topol.top

正常情况下会得到topol.tpr. 如果出现错误或警告, 检查并修正它们.

得到topol.tpr后, 我们可以使用mdrun进行能量最小化(最速下降方法)

gmx mdrun -v

能量最小化完成后检查得到的结构. 这一步中, 设定的能量收敛阈值达到与否并不重要, 因为我们只是为了消除由未正确放置的水分子导致的过大相互作用力. 你可以使用熟悉的分子查看软件来检查结构. 如果一切都看起来正常, 添加3个Cl-作为抗衡离子使体系满足电中性. 在Ewald和PME方法中, 电中性是很重要的. 我们使用genion命令来自动完成这个工作. 运行这个命令需要一个tpr文件, 因此, 我们首先需要使用grompp创建一个tpr文件, 我们使用confout.gro(即能量最小化后的结构)作为输入结构文件

gmx grompp -c confout.gro -f em.mdp -p topol.top

然后使用genion

gmx genion -s topol.tpr -p topol.top -nname CL -nn 3 -nq -1 -o neutral.gro

提示时选择溶剂SOL组.

对得到的结构neutral.gro进行能量最小化可能会更好, 但由于体系比较简单, 不进行这一步也无大碍.

接下来, 进行100 ps的平衡模拟, 我们不知道在这段时间内体系能否充分达到热力学平衡, 但我们先忽略这个问题吧.

问题: 根据哪些信息你可判定体系已经充分平衡了?

我们进行等温模拟, 使用所谓的弱耦合方法. 所用的npt.mdp文件如下. 请确保你理解其中选项的作用以及使用的原因.

integrator               = md
tinit                    = 0
dt                       = 0.002
nsteps                   = 50000
init_step                = 0
comm-mode                = Linear
nstcomm                  = 1
nstxout                  = 500
nstvout                  = 0
nstfout                  = 0
nstcheckpoint            = 1000
nstlog                   = 500
nstenergy                = 500
nstlist                  = 5
ns_type                  = grid
pbc                      = xyz
rlist                    = 0.9
domain-decomposition     = no
coulombtype              = PME
rcoulomb-switch          = 0
rcoulomb                 = 0.9
epsilon-r                = 1
vdw-type                 = Cut-off
rvdw-switch              = 0
rvdw                     = 0.9
DispCorr                 = EnerPres
fourierspacing           = 0.12
fourier_nx               = 0
fourier_ny               = 0
fourier_nz               = 0
pme_order                = 4
ewald_rtol               = 1e-05
ewald_geometry           = 3d
epsilon_surface          = 0
optimize_fft             = no
Tcoupl                   = berendsen
tc-grps                  = System
tau_t                    = 0.1
ref_t                    = 300
Pcoupl                   = berendsen
Pcoupltype               = isotropic
tau_p                    = 0.5
compressibility          = 4.5e-5
ref_p                    = 1.0
gen_vel                  = yes
gen_temp                 = 300
gen_seed                 = 1993

使用这一输入文件, 你可以创建tpr文件, 然后使用mdrun运行得到的tpr文件:

gmx grompp -c neutral.gro -f npt.mdp
gmx mdrun

自由能微扰(FEP)

平衡后, 我们进行1 ns的成品模拟, 并使用得到的轨迹进行FEP计算. 将mdp文件中的模拟步数改为1 ns对应的值, 然后进行模拟, 每ps保存一次构型. 这一模拟所用时间大约是平衡模拟的10倍.

模拟完成后, 我们使用mdrun的-rerun选项来计算轨迹中每一帧的 $V_B-V_A$. 因而-rerun所用的拓扑文件应对应于状态B(即第一个铜离子带+2电荷, 第二个铜离子带+1电荷). 将我们前面使用的状态A的拓扑文件另存为其他文件名, 如topol_B.itp, 并修改[ atom ]段的电荷使其对应于状态B.

使用状态B的拓扑文件创建新的tpr文件, 最好使用与以前不同的文件名(state_B.tpr). 更好的做法是创建一个新的目录并在其中进行新的计算. 你需要将输出选项设置为1, 即每一步都输出, 因为前面得到的轨迹仅包含每500步(即每ps)一次的快照.

得到新的tpr文件后, 将前面属于状态A的轨迹文件复制到工作目录下重命名为traj_A.trr, 然后运行:

gmx mdrun -rerun traj_A.trr -s state_B.tpr

使用energy命令抽取edr文件中的能量并选择势能(为什么?), 分别抽取状态A和状态B对应的势能. 利用这些能量, 根据Zwanzig方程计算自由能,

\[\D F=-RT\ln \left[ \left<\exp\left(-{V_B-V_A \over RT}\right)\right>_A \right]\]

问题: 自由能差值为多少? 误差来源于哪里? 为什么这么大?

误差来源于重叠不够. 我们对构型空间的抽样不充分, 其中的水分子取向倾向于产物状态.

为减小误差, 我们要使用更小的步骤进行FEP计算. 考虑根据更小的步骤来计算自由能, 如 Cu⁺+Cu²⁺ → Cu^1.1++Cu^1.9+, Cu^1.1++Cu^1.9+ → Cu^1.2++Cu^1.8+ 等等. 当然, 这种步骤在现实中是没有的, 因为不存在非整数电荷, 但由于自由能是状态函数, 所以我们可以选择任意的路径来计算它. 总的自由能是所有这些小步骤的自由能差值之和:

\[\D F=\Sum_{i=1}^n \d F_i\]

其中

\[\d F_i=-RT\ln \left[ \left<\exp\left[-{V(\text{Cu}^{(1+i/n)+}+\text{Cu}^{(2-i/n)+})-V(\text{Cu}^{(1+(i-1)/n)+}+\text{Cu}^{(2-(i-1)/n)+}) \over RT}\right]\right>_{i-1} \right]\]

如果你不擅长编程, 可以创建 $n$ 个子目录(或更多), 将状态A的拓扑文件复制到每个子目录中. 然后, 在子目录 $i$ 中将第一个铜离子的电荷改为 $1+i/n$, 第二个铜离子的电荷改为 $2-i/n$. 然后mdrun -rerun所有的模拟, 并计算总的 $\D F$.

问题: 误差变小了么? 为什么(回想分布拖尾现象)?

为了便于进行后面的教程, 保存系综(即轨迹) $i$ 中计算的 $V_{i+1}-V_i$ 到文件中.

接受比例

FEP方法是非对称的, 而且我们已经看到, 即便使用非常小的微扰步骤, 我们得到的自由能结果与参考值0仍有差距. 为解决这一问题, 我们将使用 Bennet 接受比例方法处理我们得到的FEP数据.

当满足下面的条件时, 每一FEP步骤的统计误差最小

\[\Sum_\text{sim. B} {1 \over 1+\exp[\b(V_A(r)-V_B(r)+C)]}=\Sum_\text{sim. A} {1 \over 1+\exp[\b(V_B(r)-V_A(r)-C)]}\]

累加对所有快照进行, 对应的势能保存在edr文件中(每nstenergy一次). 自由能的最佳值为满足上面方程的 $C$, 即 $\D F=C$, 由于 $C$ 出现于方程的两边, 我们只能对其进行数值求解. 最直接但可能并非最高效的求解方法是, 从小到大逐渐增加 $C$ 的值, 计算方程两边的值, 直至 $C$ 满足方程. 我们可以先使用大的间隔对 $C$ 进行扫描, 即相邻 $C$ 值的间距较大. 根据所得的初始扫描结果, 我们可以取使方程左右两边差值最小的 $C$ 值, 然后在第二次扫描中围绕此值进行更精细的扫描.

为了使用多步FEP计算中的数据, 我们也需要后向的FEP数据, 即, 从最后一点 $n$ 开始, 对每一步骤 $i$ 我们需要计算

\[\d F_i=-RT\ln \left[ \left<\exp\left[-{V(\text{Cu}^{(2-(i-1)/n)+}+\text{Cu}^{(1+(i-1)/n)+})-V(\text{Cu}^{(2-i/n)+}+\text{Cu}^{(1+i/n)+}) \over RT}\right]\right>_i \right]\]

其中的累加 $i$ 遍及 $n$ 到1. 现在我们只使用系综 $i$ 中 $V_{i-1}-V_i$ 的数据. 你需要创建并保存这些数据.

现在你有了每一步骤 $i$ 中正向和反向FEP的 $V_A$ 和 $V_B$ 数据, 查找每一步骤 $i$ 中的最佳 $C$ 值. 个人建议用Excel, MatLab或Mathematica进行. 或者, 如果你知道如何写脚本, 并知道一点awk的话, 会更容易.

问题: Bennet接受比例方法估计的自由能差值为多少?

热力学积分

在这一部分, 我们将使用GROMACS来进行热力学积分计算. 自由能的估计值为

\[\D F_{AB}=\int_1^0 \left< {\partial V \over \partial \l} \right>_\l \rmd \l \approx \Sum_i^{n_\text{steps} } \left<{\partial V \over \partial \l}\right>_\l \D \l\]

最好在新的目录中进行这一部分的计算. 如前面一样, 我们将使用10个步骤进行计算. 因此, 你需要创建10个子目录, 并将mdp文件复制到每个子目录中. 在每个mdp文件中添加下面的关键词来启用热力学积分

free-energy     = yes
init-lambda     = 0
delta-lambda    = 0
sc-alpha        = 0
sc-sigma        = 0.3

我们可以用初始init-lambda关键词控制labmda的值, 将其设置为0, 0.1, 0.2,…1, 对应于每一子目录.

使用这些设置, GROMACS会对状态A和B之间的相互作用函数进行线性内插. 因此我们同时需要设定状态B, 只要简单地在拓扑文件中添加状态B的信息(电荷)即可

[ atoms ]
;   nr    type  resnr residue  atom  cgnr charge  mass    typeB  chargeB  massB
    1     CU1+      1   CU1+   CU    1    1       63.546   CU1+   2       63.546
    2     CU1+      2   CU2+   CU    2    2       63.546   CU1+   1       63.546

在每一子目录中进行固定 $\l$ 值的模拟, mdrun的输出文件之一为dhdl.xvg, 其中包含了随模拟时间变化的 ${\partial V \over \partial \l}$. 你可以使用analyze命令来计算平均值和误差估计随 $\l$ 的变化. 当得到 $\left<{\partial V \over \partial \l}\right>_\l$ 与 $\l$ (0到1)的关系后, 你也可以使用analyze命令计算从0到1的积分.

问题: 热力学积分给出的自由能为多少? 是否好于FEP和BAR?

Jarzynski方法

创建一个新的子目录, 并复制用于热力学积分计算的mdp文件和拓扑文件. 在本部分中, 我们将进行很短时间的慢增长模拟. 慢增长是热力学积分的连续形式, 其中的 $\l$ 作为时间的函数连续变化, 因此, $\l(t_{init}) = 0$, $\l(t_{final}) = 1$. 在mdp文件中, 增加下面的选项(在这个例子中我们进行20 ps的模拟)

dt               = 0.002
nsteps           = 10000

free-energy      = yes
init-lambda      = 0
delta-lambda     = 0.0001

我们创建许多子目录, 每一个都进行短时间(几ps)的慢增长模拟. 如果你会写脚本的话, 会更容易.

执行完所有模拟后, 可以使用anlyze程序计算所有模拟的非可逆功

gmx analyze -f dhdl.xvg -integrate

注意, x值是以ps为单位的时间, 从0到20. 因此你需要将所得的功除以20使得 $\l$ 处于0到1之间, 间隔正确. 得到这些功后, 你可以使用你熟悉的数据分析程序创建一个直方图, 并试着将其拟合为高斯函数, 使用Jarzynski 方程(指数平均)来计算可逆功

\[\exp\left[-\b \D F \right]={1\over n_\text{sim.} } \Sum_{n_\text{sim.} } \exp\left[-\b \int_{t_i}^{t_f} {\partial V \over \partial \l} {\rmd \l \over \rmd t}\rmd t \right]\]

问题: 要得到准确的自由能估计值, 需要多少个慢增长步骤?

结论

我们已经使用不同方法计算了自由能. 显然, 不应该使用非对称的自由能微扰方法, 但当它与对称的Bennet接受比例方法一起使用时, 能够给出精确的自由能估计值. 热力学积分也是对称的, 当 $\l$ 间隔足够小, 相应 $\l$ 系综间有重叠时, 这种方法同样能够给出精确的自由能估计值. 前面三种方法假定体系哈密顿量的变化足够缓慢, 体系始终处于平衡态, 因此, 这些方法都被称为平衡方法. 然而, 由于采样时间的限制, 实际情况很少这样. 与平衡方法相反, Jarzynski方法用于远离平衡的体系. 尽管考虑到我们的限制, 这听起来不错, 但由于平均问题, 需要非常多的模拟才能得到精确的自由能估计值. 最后的结果就是你并不能节省多少模拟时间. 这样, 具体使用哪种方法依赖于你要研究的问题, 要由你来决定使用哪种方法.

在此教程中, 我们仅仅改变了电荷, 但所用方法是通用的, 我们也可以改变成键相互作用. 然而, 对于构型改变, 基于伞形采样的其他方法通常更合适.