【3.2.6】gromacs教程--自由能
自由能计算：水中的甲烷 Free Energy Calculations: Methane in Water
source /opt/fastone/softwares/hpc-kits/hpc-kits.sh
source /data/software/gromacs/gromacs-2020.3/bin/GMXRC
cd /data/user/sam/project/gromacs/test6
一、理论
本教程将假定您对什么是自由能计算，存在的不同类型以及该技术的基本理论有一个合理的了解。 在这里提供全面的教育既不可行。 相反，我将本教程的重点放在GROMACS中运行自由能计算的实践方面，并在整个教程中重点介绍了相关的理论要点。 我将在这里为那些对免费能源计算不熟悉的人提供有用的论文清单。 不要认为此列表是详尽无遗的。 它还不应取代对统计力学以及对相关主题的许多书籍和论文的正确研究。
C. D. Christ, A. E. Mark, and W. F. van Gunsteren (2010) J. Comput. Chem. 31: 1569-1582. DOI
A. Pohorille, C. Jarzynski, and C. Chipot (2010) J. Phys. Chem. B 114: 10235-10253. DOI
A. Villa and A. E. Mark (2002) J. Comput. Chem. 23: 548-553. DOI
本教程的目的是重现非常简单的系统的结果，对于该系统，存在准确的自由能估算值。 选择系统（水中的甲烷）由Shirts等人处理。 在力场和氨基酸侧链类似物的水合自由能的系统研究中。 完整的出版物可以在这里找到。 本教程将假定您已阅读并理解本文的更广泛内容。
此处进行的数据分析将使用GROMACS 4.5版（以前称为g_bar）中引入的GROMACS“ bar”模块，而不是使用热力学积分方法来评估自由能差。 它使用Bennett接受率（BAR，因此是模块的名称）方法来计算自由能差。 BAR的相应论文可以在这里找到。 还假定了此方法的知识，此处将不对其进行详细讨论。
自由能计算有许多实际应用，其中一些更常见的应用包括溶剂化/水合的自由能以及小分子与一些较大的受体生物分子（通常是蛋白质）结合的自由能。这两个过程都涉及需要从系统中添加（引入/偶联）或去除（解偶联/消除）所需的小分子，并计算所得的自由能变化。
在自由能计算期间，可以转换两种类型的非键相互作用，即库仑相互作用和范德华相互作用。绑定的交互也可以进行操作，但为简单起见，此处不再赘述。在本教程中，我们将计算一个非常简单的过程的自由能：关闭水中目标分子（在本例中为甲烷）的原子位置之间的Lennard-Jones相互作用。 Shirts等人非常精确地计算了该数量。因此，这就是我们要在此处复制的数量。
二、分析步骤
2.2 Examine the Topology 检查拓扑
下载该系统的坐标文件和拓扑。 这些文件是该系统的David Mobley教程的一部分（不再在线），它们是Michael Shirts在前一页引用的论文中使用的原始文件（为了与最新的GROMACS版本兼容而略作修改）。
该系统在一个由596个TIP3P水分子组成的盒子中包含一个甲烷分子（在坐标文件中，丙氨酸的β-碳称为“ ALAB”）。 查看拓扑，我们发现：
; Topology for methane in TIP3P
#include "oplsaa.ff/forcefield.itp"
[ moleculetype ]
; Name nrexcl
Methane 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB massB
1 opls_138 1 ALAB CB 1 0.000 12.011
2 opls_140 1 ALAB HB1 2 0.000 1.008
3 opls_140 1 ALAB HB2 3 0.000 1.008
4 opls_140 1 ALAB HB3 4 0.000 1.008
5 opls_140 1 ALAB HB4 5 0.000 1.008
[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1
1 3 1
1 4 1
1 5 1
[ angles ]
; ai aj ak funct c0 c1 c2 c3
2 1 3 1
2 1 4 1
2 1 5 1
3 1 4 1
3 1 5 1
4 1 5 1
; water topology
#include "oplsaa.ff/tip3p.itp"
[ system ]
; Name
Methane in water
[ molecules ]
; Compound #mols
Methane 1
SOL 596
您会注意到，所有charge都设置为零。此设置背后有一个实际的原因。通常，溶质与水之间的电荷相互作用在范德华方程之前被关闭。如果在Lennard-Jones项关闭时保持电荷相互作用，则正电荷和负电荷将被允许无限近地相互接近，从而导致一个非常不稳定的系统，可能会爆炸。本教程中的过程基本上假设在此之前已正确关闭了电荷，并且我们将仅关闭溶质和溶剂之间的范德华相互作用。
请注意，在以前的GROMACS版本中，B型，chargeB和massB标题的内容必须与分子的B状态相对应。从GROMACS 4.0版开始，不再需要对简单（去耦合）进行拓扑修改（万岁！），但是在将一个分子突变为另一个分子（其中键合和非键合的术语可能会发生变化）的情况下，会进行旧式修改（因此-称为“双重拓扑方法”）。 GROMACS手册的第5.8.4节提供了这种转换的示例。
2.3 工作流程
将系统从状态A转换为状态B的自由能变化ΔGAB是耦合参数λ的函数，该参数表示状态A和状态B之间发生的变化的程度，即程度 哈密顿量受到干扰，系统也发生了变化。 在不同的λ值下进行的仿真使我们能够绘制∂H/∂λ曲线，从中得出ΔGAB。 计划自由能计算的第一步是要使用多少个λ点来描述从状态A（λ= 0）到状态B（λ= 1）的转换，目的是收集足够的数据以详尽地采样相空间和 产生可靠的∂H/∂λ曲线。
对于线性依赖于λ的库仑相互作用的解耦，通常应从0到1等距的λ间距就足够了，通常0.05到0.1的λ间距值即可。请注意，这是一个广义的概括，实际上，许多分子将需要更多的技巧，例如那些通过氢键与周围环境发生强烈相互作用的分子。在这种情况下，可能需要减小λ间距，从而可能甚至以非对称方式使用更多点。
对于范德华相互作用的去耦，λ间距可能会更加棘手。由于Shirts et al.描述的原因，在其他地方，可能需要大量的λ点才能正确描述转换。可能需要在∂H/∂λ曲线的斜率较陡的区域中聚集λ点。在本教程中，我们将使用等距的λ间距0.05，但是在许多情况下，可能需要使用不同的间距，尤其是在0.6≤λ≤0.8范围内的聚类值。
对于每个λ值，必须执行完整的工作流程（能量最小化，平衡和数据收集）。通常，我发现批量运行这些作业是最有效的，将一步的输出直接传递到下一步。这里使用的工作流程将是：
- 最陡的下降最小化
- NVT平衡
- NPT平衡
- Data collection under an NPT ensemble
以下链接提供了教程不同步骤的.mdp文件（仅适用于λ= 0）：
我还提供了一个附件Perl脚本write_mdp.pl，您可以使用该脚本非常快速地生成运行这些模拟所需的所有输入文件。
您还可以根据我将描述的工作流程下载用于运行这些作业的Shell脚本（job.sh）。 如果您发出：
perl write_mdp.pl em_steep.mdp
perl write_mdp.pl npt.mdp
perl write_mdp.pl nvt.mdp
perl write_mdp.pl md.mdp
您将收到名为em_steep_0.mdp，em_steep_1.mdp，… em_steep_20.mdp的文件，并将init_lambda_state的相关值插入文件名中。 类似地，您可以用相同的方式（nvt.mdp，npt.mdp和md.mdp）产生其余的.mdp文件。
在job.sh中，您将需要更改 $FREE_ENERGY 和 $MDP 的值，否则该脚本肯定无法工作。
有关.mdp设置的注意事项
在继续之前，重要的是要了解.mdp文件中与自由能有关的设置（以em_steep_0.mdp为例）。 对于工作流中的所有步骤，它们都是相同的。 我将假定您熟悉.mdp文件中的其他设置。 如果不是这种情况，请在继续之前参考一些更基本的教程材料。
| 参数 | 说明 |
|---|---|
| free_energy = yes | 表示我们正在做自由能计算，并且将在所选分子的A和B状态（定义如下）之间进行插值。 |
| init_lambda_state = 0 | 在以前的GROMACS版本中，“ init_lambda”关键字直接指定单个λ值。 从5.0版开始，λ现在是一个向量，可以转换键和非键相互作用。 使用init_lambda_state关键字，我们指定要在仿真中使用的向量的索引（从零开始）（稍后将对此进行详细介绍）。 |
| delta_lambda = 0 | 可以在称为“缓慢增长”的技术中，每个时间步长将λ的值增加一些量（即δλ/δt）。该方法可能会伴有很大的误差，因此我们将不会随时间改变λ值。 |
| fep_lambdas =（nothing） | 您将注意到未指定此关键字。在以前的GROMACS版本中，init_lambda和foreign_lambda的使用控制λi的值以及λi的附加值，对于这些值，将针对λi的配置评估能量差。这已不再是这种情况。可以在fep_lambdas关键字中显式设置λ的值，但相反，我们允许calc_lambda_neighbors关键字（请参见下文）自动确定这些其他值。 |
| calc_lambda_neighbors = 1 | 相对于λi计算出能量差的相邻窗口的数量。例如，如果init_lambda_state设置为10，则在运行期间以calc_lambda_neighbors = 1计算与λ状态9和11有关的能量差。 |
| vdw_lambdas = … | 用于范德华相互作用转换的λ值数组。 |
| coul_lambdas = … | 用于库仑（静电）相互作用转换的λ值数组。 |
| bonding_lambdas = … | 用于键合相互作用转换的λ值数组。 |
| inhibitort_lambdas = … | 用于约束转换的λ值数组。 |
| mass_lambdas = … | 用于原子质量转换的λ值数组； 用于改变分子的化学性质。 |
| temperature_lambdas = … | 用于温度转换的λ值数组； 仅用于模拟回火。 |
| sc-alpha = 0.5 | “软核” Lennard-Jones方程中使用的α比例因子的值。 有关此术语以及后续三个术语的完整说明，请参见GROMACS手册第4.5.1节中的公式4.124-4.126。 |
| sc-coul =no | 线性变换库伦相互作用。 这是默认行为，但为清楚起见而写出来。 |
| sc-power = 1.0 | 软核方程中λ的幂。 |
| sc-sigma = 0.3 | 分配给具有C6或C12参数等于零或σ<sc-sigma（通常为H原子）的任何原子类型的σ值。 软核Lennard-Jones公式中使用了该值。 |
| couple-moltype = Methane | 其中的[分子类型]名称将从状态A插入到状态B。请注意，此处给出的名称必须与[分子类型]名称匹配，而不与残基名称匹配。在某些情况下，这些可能是相同的，但出于指导目的，我为[分子类型]和组分残基保留了不同的名称。 |
| couple-lambda0 = vdw | 插值[分子类型]与系统其余部分之间在状态A中存在的非键相互作用的类型。值“ vdw”表示只有范德华项在甲烷和水之间有效。没有溶质-溶剂库仑相互作用。 |
| couple-lambda1 = none | 插值[分子类型]与系统其余部分之间在状态B中存在的非键相互作用的类型。值“无”表示状态B中的范德华相互作用和库仑相互作用都处于关闭状态。相对于couple-lambda0，这表明仅范德华项已关闭。 |
| couple-intramol =no | 不解耦分子内相互作用。即，λ因子仅应用于溶质-溶剂非键相互作用，而不应用于溶质-溶质非键相互作用。设置“ couple-intramol = yes”对可能在较长距离内发生分子内相互作用的较大分子（例如肽或其他聚合物）很有用。否则，分子的远端部分将以不自然的强烈方式通过显式对相互作用进行相互作用（因为溶质与溶剂之间的相互作用会随着λ的作用而减弱，但分子内术语将不会），从而导致构型错误地偏向由此产生的自由能。 |
| nstdhdl = 10 | 将∂H/∂λ和ΔH写入dhdl.xvg的频率。值100可能就足够了，因为结果值将高度相关并且文件将变得非常大。您可能希望自己的工作将此值增加到100。 |
与Shirts等人使用的设置相比，此处将使用的.mdp设置也存在一些差异：
- rlist = rcoulomb = rvdw = 1.2。为了使用PME，rlist必须等于rcoulomb，这是在3.3.1版发行后的某个时间在grompp中强制执行的检查。在版本4.6中引入的Verlet截止方案允许缓冲邻居列表，也需要rvdw = rcoulomb。为了节省能源，将在运行期间调整rlist的值，从而为切换的范德华相互作用提供必要的缓冲。转换范围（从1.0-1.2 nm）与Shirts等人的方法一致。用于溶质-溶剂相互作用的处理。
- 通过使用Langevin积分器（“ sd”设置指定“随机动力学”，例如具有随机摩擦力）来处理温度耦合，而不是使用Andersen恒温器。
- tau_t = 1.0。重要的是要说明此特定值，因为当使用Langevin积分器时，该值对应于反摩擦系数，因此为ps-1。值1.0避免过度抑制水的动力。
让我们花一点时间仔细剖析λ向量。例如，对于init_lambda_state = 0，这意味着我们指定转换中的状态与* _lambdas关键字中索引为0的向量相对应，如下所示：
; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
vdw_lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
coul_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
bonded_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
restraint_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
mass_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
temperature_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
设置initial_lambda_state = 1将对应于右边的下一列（范德瓦尔斯的λ= 0.05），而initial_lambda_state = 20将指定最后一列（λ向量），其中范德华斯的λ值= 1.0。 在本教程中，我们仅转换范德华相互作用，而与电荷，键合参数，约束等相关的所有事物都保留下来。 请注意，在此示例中，当拓扑中的电荷为零且.mdp设置从不指示应该存在电荷时（couple-lambda0或couple-lambda1都不包含’q'），coul_lambdas的值选择无关紧要 ，但出于学究的目的，.mdp文件清楚地表明没有进行电荷转换。
2.4 能源最小化 Energy Minimization
我提供的用于运行这些计算的job.sh脚本将创建以下目录层次结构：
Lambda_0/
Lambda_0/EM/
Lambda_0/NVT/
Lambda_0/NPT/
Lambda_0/Production_MD/
这样，对于init_lambda_state的每个值，工作流中的所有步骤都在单个目录内执行。 我发现这是组织作业及其输出的便捷方法。
该脚本还假定.mdp文件也在$ MDP目录中按层次结构组织，该目录在脚本中设置为环境变量：
$MDP/
$MDP/EM/
$MDP/NVT/
$MDP/NPT/
$MDP/Production_MD/
如前所述，将使用最速下降法进行能量最小化。 job.sh脚本中的相关部分是：
mkdir Lambda_$LAMBDA
cd Lambda_$LAMBDA
#################################
# ENERGY MINIMIZATION 1: STEEP #
#################################
echo "Starting minimization for lambda = $LAMBDA..."
mkdir EM
cd EM
# Iterative calls to grompp and mdrun to run the simulations
gmx grompp -f $MDP/em_steep_$LAMBDA.mdp -c $FREE_ENERGY/methane_water.gro
-p $FREE_ENERGY/topol.top -o min$LAMBDA.tpr
gmx mdrun -deffnm min$LAMBDA
echo "Minimization complete."
2.5 Equilibration 平衡
系统将分两阶段进行平衡，第一阶段为恒定体积（NVT），第二阶段为恒定压力（NPT）。
#####################
# NVT EQUILIBRATION #
#####################
echo "Starting constant volume equilibration..."
cd ../
mkdir NVT
cd NVT
gmx grompp -f $MDP/nvt_$LAMBDA.mdp -c ../EM/min$LAMBDA.gro
-p $FREE_ENERGY/topol.top -o nvt$LAMBDA.tpr
gmx mdrun -deffnm nvt$LAMBDA
echo "Constant volume equilibration complete."
#####################
# NPT EQUILIBRATION #
#####################
echo "Starting constant pressure equilibration..."
cd ../
mkdir NPT
cd NPT
gmx grompp -f $MDP/npt_$LAMBDA.mdp -c ../NVT/nvt$LAMBDA.gro
-p $FREE_ENERGY/topol.top -t ../NVT/nvt$LAMBDA.cpt -o npt$LAMBDA.tpr
gmx mdrun -deffnm npt$LAMBDA
echo "Constant pressure equilibration complete."
2.6 Production MD
既然系统已经达到平衡，我们就可以开始生产MD模拟，在此期间我们将收集∂H/∂λ数据。
#################
# PRODUCTION MD #
#################
echo "Starting production MD simulation..."
cd ../
mkdir Production_MD
cd Production_MD
gmx grompp -f $MDP/md_$LAMBDA.mdp -c ../NPT/npt$LAMBDA.gro
-p $FREE_ENERGY/topol.top -t ../NPT/npt$LAMBDA.cpt -o md$LAMBDA.tpr
gmx mdrun -deffnm md$LAMBDA
echo "Production MD complete."
在GPU工作站上，整个工作流程大约需要2.5个小时。 由于要运行许多模拟（21组λ向量），因此最好在有更多处理器可用的集群上运行这些作业，以便可以同时运行作业。 一旦完成所有生产模拟，我们就可以分析结果数据。
2.7 分析
GROMACS的条形模块使ΔGAB的计算非常简单。 只需在工作目录中收集mdrun生成的所有md * .xvg文件（每个λ值一个），然后运行gmx bar：
gmx bar -f md*.xvg -o -oi
该程序除了产生三个输出文件：bar.xvg，barint.xvg和histogram.xvg外，还将在屏幕上打印许多有用的信息。 gmx栏的屏幕输出应类似于：
Detailed results in kT (see help for explanation):
lam_A lam_B DG +/- s_A +/- s_B +/- stdev +/-
0 1 0.08 0.00 0.03 0.00 0.03 0.00 0.25 0.00
1 2 0.07 0.00 0.03 0.00 0.03 0.00 0.25 0.00
2 3 0.06 0.00 0.03 0.00 0.04 0.00 0.26 0.00
3 4 0.04 0.00 0.04 0.00 0.05 0.00 0.28 0.00
4 5 0.02 0.00 0.03 0.00 0.04 0.00 0.29 0.00
5 6 0.00 0.00 0.05 0.00 0.06 0.00 0.30 0.00
6 7 -0.02 0.00 0.05 0.00 0.06 0.00 0.32 0.00
7 8 -0.06 0.01 0.06 0.00 0.07 0.00 0.35 0.01
8 9 -0.09 0.01 0.06 0.00 0.07 0.00 0.37 0.00
9 10 -0.14 0.01 0.08 0.00 0.10 0.01 0.41 0.01
10 11 -0.22 0.01 0.09 0.01 0.12 0.01 0.46 0.00
11 12 -0.31 0.01 0.12 0.01 0.15 0.01 0.51 0.01
12 13 -0.45 0.01 0.17 0.01 0.20 0.01 0.57 0.00
13 14 -0.59 0.01 0.16 0.02 0.17 0.02 0.58 0.01
14 15 -0.62 0.01 0.11 0.01 0.11 0.01 0.47 0.00
15 16 -0.54 0.01 0.06 0.01 0.06 0.01 0.35 0.00
16 17 -0.42 0.00 0.03 0.01 0.03 0.01 0.25 0.00
17 18 -0.28 0.00 0.02 0.00 0.02 0.00 0.18 0.00
18 19 -0.16 0.00 0.01 0.00 0.01 0.00 0.13 0.00
19 20 -0.05 0.00 0.00 0.00 0.00 0.00 0.10 0.00
Final results in kJ/mol:
point 0 - 1, DG 0.19 +/- 0.00
point 1 - 2, DG 0.17 +/- 0.01
point 2 - 3, DG 0.14 +/- 0.00
point 3 - 4, DG 0.09 +/- 0.01
point 4 - 5, DG 0.06 +/- 0.00
point 5 - 6, DG 0.01 +/- 0.01
point 6 - 7, DG -0.06 +/- 0.01
point 7 - 8, DG -0.14 +/- 0.02
point 8 - 9, DG -0.23 +/- 0.01
point 9 - 10, DG -0.35 +/- 0.01
point 10 - 11, DG -0.53 +/- 0.01
point 11 - 12, DG -0.77 +/- 0.03
point 12 - 13, DG -1.12 +/- 0.03
point 13 - 14, DG -1.46 +/- 0.02
point 14 - 15, DG -1.53 +/- 0.04
point 15 - 16, DG -1.35 +/- 0.02
point 16 - 17, DG -1.04 +/- 0.01
point 17 - 18, DG -0.70 +/- 0.01
point 18 - 19, DG -0.40 +/- 0.00
point 19 - 20, DG -0.13 +/- 0.00
total 0 - 20, DG -9.13 +/- 0.09
获得的ΔGABI值为-9.13±0.09 kJ mol-1或-2.18±0.02 kcal mol-1。 由于我为该演示进行的过程是甲烷的解偶联，因此反向过程（将不带电荷的甲烷引入水中，因此该过程的实际水合能量）对应于2.18±0.02 kcal mol-1的ΔGAB（假设可逆性） ），与Shirts等人获得的价值非常一致。 尽管我使用的范德华变换的λ向量数量大约是原始作者的一半，但我仍然使用2.24±0.01 kcal mol-1（根据该论文的支持信息表II）进行计算。 前面描述的其他更改。
从技术上讲，这就是我们要做的所有事情，但是gmx bar还可以打印一些有用的输出文件（所有这些文件都是可选的）。 它们的内容在这里值得探讨，因为它们提供了去耦过程和采样成功的详细信息。
2.7.1 输出文件1：bar.xvg
让我们看一下gmx bar给我们的其他输出文件，从bar.xvg开始。 该输出文件绘制了每个λ间隔（即相邻的汉密尔顿之间的相对自由能）的相对自由能差，并且应看起来像这样（重新绘制为条形图而不是默认线表示，并添加一些网格线后） 透视）：
垂直的灰线表示在此进行的去耦过程中使用的λ值。因此，每个黑条表示相邻的λ值之间的自由能差。在bar.xvg中，我们首先了解calc_lambda_neighbors在模拟过程中的工作。例如，在我们对init_lambda_state = 0的仿真中，每nstdhdl（10）个步骤评估λ= 0.05（由calc_lambda_neighbors = 1指定的最接近的相邻窗口）处的自由能。这些“外来”λ计算导致λ值之间发生能量重叠，从而使我们具有相空间重叠和足够的采样。有兴趣的读者应参考gmx条描述中引用的Wu和Kofke的论文，以讨论此概念以及对熵值sA和sB的解释。
因此，从λ= 0到λ= 1的自由能变化仅是每对相邻的λ模拟的自由能变化之和，以bar.xvg表示。
ΔG的值对应于由gmx条打印到屏幕上的输出的前半部分。屏幕输出的后半部分包含这些相同的ΔG值，已转换为kJ mol-1（在298 K时为1 kT = 2.478 kJ mol-1）。
2.7.2 输出文件2：barint.xvg
barint.xvg文件将累积ΔG绘制为λ的函数。 在barint.xvg中，λ= 1处的点对应于从λ向量0到λ向量1的ΔG之和，如上面的屏幕输出所示：
lam_A lam_B DG +/- s_A +/- s_B +/- stdev +/-
0 1 0.08 0.00 0.03 0.00 0.03 0.00 0.25 0.00
1 2 0.07 0.00 0.03 0.00 0.03 0.00 0.25 0.00
因此，在barint.xvg中在λ= 0处的累积ΔG的值为零，在λ= 0.1处的值为0.0 + 0.08 = 0.08，这就是我们在barint.xvg中找到的值。 下面是barint.xvg（蓝色）与bar.xvg的内容（黑色，λ值之间的ΔG值）旁边的图，用于说明累积ΔG。 通过取λ= 20（-3.69 kT）时的ΔG值，我们可以恢复上面显示的ΔG值（-9.14 kJ mol-1）。 此处的值（-9.14）与上方的值（-9.13）之差是由于仅使用两位小数而产生的舍入误差，而bar程序在上面的计算中使用了全精度，但是所有输出仅打印两位小数 。
2.8 进阶应用
计算自由能的一种常见应用是确定配体与受体（例如蛋白质）之间的结合ΔG，您需要将配体与受体复合并在本体溶液中进行（去偶联）， 因为在这种情况下ΔG是配体和受体络合的自由能变化（即将配体引入结合位点）和使配体解溶剂（即从溶液中除去）的自由能之和：
ΔGbind = ΔGcomplexation + ΔGdesolvation
ΔGsolvation = -ΔGdesolvation
∴ ΔGbind = ΔGcomplexation - ΔGsolvation
结合的ΔG将是这两个状态之间的差，并且可以通过使用以下（略微简化的）热力学循环来计算，其中在以上等式中，ΔG1为ΔG溶剂化，而ΔG3为ΔG复合。
可以在本文的图1中找到针对此类过程的特别详细的热力学循环。
完全相互作用的物种（例如配体）转变为“虚拟”（配体构型中的一组原子中心，不与周围环境发生任何未键合的相互作用）需要同时关闭两个范德华力 相互作用（如本教程中所示）以及感兴趣的溶质与其周围环境之间的库仑相互作用。 甲烷的完整转化系列如下所示：
在上述热力学循环中，ΔG1和ΔG3实际上代表配体的引入（从虚拟分子开始），即偶联而不是去偶联。 上面的热力学循环是极其基本的，并且没有考虑诸如方向限制，蛋白质构象变化等因素。本文将不讨论任何文献。 对于大多数小分子（以下通常称为“ LIG”），以下设置似乎很有吸引力：
couple-moltype = LIG
couple-intramol = no
couple-lambda0 = none
couple-lambda1 = vdw-q
在这种情况下，必须格外小心。 在以前版本的GROMACS中，这种方法可能非常不稳定，因为在系统保留一些电荷的同时除去范德华项，可以使带相反电荷的原子非常紧密地相互作用（在没有电子云的全部作用的情况下） 。 结果将是不稳定的配置和不可靠的能量，即使系统没有崩溃也是如此。 因此，顺次进行（解）耦合会更好。 在版本5.0中，使用新的λ向量功能可以很容易地做到这一点：
couple-moltype = LIG
couple-intramol = no
couple-lambda0 = none
couple-lambda1 = vdw-q
init_lambda_state = 0
calc_lambda_neighbors = 1
vdw_lambdas = 0.00 0.05 0.10 ... 1.00 1.00 1.00 ... 1.00
coul_lambdas = 0.00 0.00 0.00 ... 0.00 0.05 0.10 ... 1.00
在这种情况下，库仑相互作用的λ值始终为零，而变换范德华相互作用的λ值则发生变化。 然后，范德华相互作用完全打开（λ= 1.00），而库仑相互作用逐渐打开。 跟踪λ= 0和λ= 1对应的状态是此过程的关键。 在上述情况下，couple-lambda0表示交互已关闭，而couple-lambda1表示交互已打开。
这些类型的转换可用于计算水合/溶剂化自由能，因为该量通常在力场参数化中用作目标数据。
参考资料
