Methane in Water
Introduction
本教程将指导用户通过计算一个简单的自由能变化的过程（中性甲烷和盒内水分子之间的范德华相互作用的去除解耦过程）。 这个量曾由 Shirts 等人在系统地评估氨基酸侧链的溶剂化自由能时用不同的原子力场计算出来。
本教程需要 GROMACS 2018.x 系列的版本。
Theory
本教程将假设您对什么是自由能计算、存在的不同类型以及该技术的基本理论有一个清除的认识。在这里提供完整的讲解既不现实也不可能。相反，我将把本教程的重点放在通过实际操作实现在GROMACS中运行自由能的计算，并在整个教程中强调相关且必要的理论要点。作者在这里提供一个有用的论文列表，以供那些不太熟悉自由能的新手参考。不要认为这个列表是详尽的，它也不应取代那些关于统计力学的研究论著以及类似主题的书籍和论文。
- C. D. Christ, A. E. Mark, and W. F. van Gunsteren (2010) J. Comput. Chem. 31: 1569-1582.
- A. Pohorille, C. Jarzynski, and C. Chipot (2010) J. Phys. Chem. B 114: 10235-10253.
- A. Villa and A. E. Mark (2002) J. Comput. Chem. 23: 548-553.
本教程的目的是重现一个非常简单的系统的结果，该系统存在一个精确的自由能判断结果。所选择的系统(methane in water)是由 Shirts 等人在系统研究氨基酸侧链类似物的力场和水合自由能时所处理的系统。完整的文章可以在这里找到。本教程将假设您已经阅读并理解了此文甚至更广泛的观点。
这里进行的数据分析将使用GROMACS“bar”模块，而不是使用热力学积分方法来评估自由能差异，该模块在GROMACS 4.5版本中引入(以前称为g_bar)。它使用 Bennett Acceptance Ratio （BAR，模块名称的由来）方法来计算自由能差。BAR对应的论文在这里。本文还假设您对这种方法有一定的了解，这里将不进行详细讨论。
自由能计算有许多实际应用，其中一些更常见的包括溶剂化/水合自由能，以及小分子与一些较大的受体生物分子（通常是蛋白质）结合的自由能。这两种过程都需要从系统中添加（引入/耦合）或移除（解耦/湮灭）相应的小分子，并计算由此产生的自由能变化。
在自由能计算中，有两种类型的非键相互作用可以转换，即库仑相互作用和范德华相互作用。键相互作用也可以被操控，但是为了简单起见，这里不讨论。在本教程中，我们将计算一个非常简单的过程的自由能：消解一个的分子（在本例中是甲烷）的原子位点之间的 Lennard-Jones 相互作用。这个数值是由 Shirts 等人非常精确地计算出来的，因此这是我们在这里想要复刻这个结果。
The Topology
下载该系统的坐标文件和拓扑结构。这些文件是作为 David Mobley 关于这个系统的教程的一部分提供的（网上已经无了），并且是 Michael Shirts 在前页引用的论文中使用的原始文件（为了与最近的GROMACS版本兼容做了些修改）。
该系统包含一个单分子甲烷（坐标文件中称为“ALAB”，代表丙氨酸的β-碳）以及同一个盒中596个TIP3P水分子。查看拓扑结构，我们发现:
; 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
你会注意到所有的电荷都被设为零。这是因为，通常溶质和水之间的电荷相互作用在范德华力之前是不起作用的。如果 Lennard-Jones 被屏蔽，而电荷相互作用仍在继续，正电荷和负电荷将会无限近地相互靠近，从而导致系统非常不稳定，很可能会使计算量爆炸。本教程假设电荷在此之前已经被适当地屏蔽了，我们将只屏蔽溶质和溶剂之间的 van der Waals 相互作用。
请注意，在GROMACS的早期版本中，typeB、chargeB和massB标题的内容必须对应于分子的 B-state。从GROMACS 4.0版本开始，不再需要为简单（解）耦合修改拓扑（hooray!），但在一个分子突变为另一个分子的情况下，其中有键和非键项可能会改变，仍然需要老式的修改（所谓的“双拓扑方法”）。GROMACS手册第5.8.4节提供了一个转换的示例。
Free Energy Workflow
由 A 状态转变为 B 状态的自由能变 ΔGAB 是是一个耦合参数 λ 的函数，它代表了状态 A 和状态 B 之间发生的变化程度，即哈密顿量被扰动的程度和系统被改变的程度。通过在 λ 的不同值上进行模拟画出相应 ∂H/∂λ 曲线，然后我们就可以由此得到 ΔGAB。 规划自由能计算的第一步是确定我们要选取多少个λ来描述从状态 A(λ= 0) 到状态 B(λ= 1) 的转换过程，进而收集到足够多的数据来彻底对所有相空间采样并生产可靠 ∂H/∂λ 曲线。
-
对于解耦的库仑相互作用（decoupling Coulombic interactions），它与 λ 线性相关；因此以0.05-0.1为间距等间地选取 λ 通常足以满足我们的需要。请注意，这只是一个泛泛的概括，事实上，许多分子需要更灵巧的处理，例如那些通过氢键与周围环境强烈相互作用的分子。在这种情况下，λ 的间距可能需要减少以使用更多的点，甚至可能需要采用不对称取点的方式。
-
对于解耦范德华相互作用（decoupling van der Waals interactions），他与 λ 关系更加复杂（间距选取更棘手）。根据 Shirts 等人和其他文章所描述的原因，我们可能需要大量的 λ 点来恰当地描述这种转换。在 ∂H/∂λ 曲线的斜率陡峭的区域增加 λ 采样点的数量可能是很关键的。为了本教程的目的，我们将使用等距的 λ 间隔 0.05，但在许多情况下，可能需要使用不同的间隔，特别是在0.6≤λ≤0.8的范围内增加 λ 采样（clustering values）。
对于每一个λ值，都必须进行完整的工作流程 Workflow（能量最小化、平衡、数据收集）。作者发现将这些作业按批次处理运行，将一个步骤的输出直接传递给下一个步骤是最有效的。这里使用的 Workflow 将是：
- 最速下降法的能量最小化
- NVT平衡
- NPT平衡
- 在NPT集合下收集数据
对于教程的不同步骤(仅针对λ = 0)的.mdp文件在以下链接：
作者还附加了一个 Perl 脚本 write_mdp.pl，您可以使用它来快速地生成运行这些模拟所需的所有输入文件（.mdp）（这个脚本仅可按照间隔 0.05 来生成 .mdp 文件，当然你可以修改脚本的相应参数来满足你自己的需要）。您还可以下载这个 shell 脚本（job.sh，自动运行上述 workflow 的脚本，建议使用前先看一看，熟悉一下），以便根据作者将描述的 workflow 运行这些工作。
输入：
perl write_mdp.pl em_steep.mdp
您将得到名为 em_steep_0.mdp, em_steep_1.mdp, … em_steep_20.mdp 的一系列文件，并且相应的 init_lambda_state 值已经插入到文件名中。类似地，您可以使用相同的方法得到其他步骤所需的 .mdp文件（nvt.mdp, npt.mdp, and md.mdp）。
在job.sh中，您需要更改 $FREE_ENERGY 和 $MDP 的值，否则脚本无法工作。
Notes on .mdp settings
在继续之前，务必理解.mdp文件中与自由能量相关的设置（使用 em_steep_0.mdp 为例）。这些参数对于 workflow 中的所有步骤都是相同的。假设您熟悉 .mdp 文件中的其他设置（如果您不熟悉，请在继续之前参考一些更基本的教程材料）。
| free_energy = yes | 声明我们正在进行自由能计算并将在所选分子（稍后解释）的 A 态和 B 态之间进行插值 |
| init_lambda_state = 0 | 在以前的GROMACS版本中，“init_lambda”关键字直接指定了 λ 的单个值。自5.0版，λ 成为了一个矢量并且允许键和非键相互作用的转换。通过 init_lambda_state，我们指定了要在模拟中使用的向量的索引（从0开始）(稍后详细介绍) |
| delta_lambda = 0 | 在称为“缓慢增长”的技术中，λ 的值可以在每个时间步中有一定增加(例如 δλ/δt 值不恒为0时)。这种方法可能出现严重的错误，因此我们将没有时间依赖性的变化应用到我们的 λ 值上 |
| fep_lambdas = (nothing) | 注意，这个关键字没有指定任何参数。在以前的GROMACS版本中，init_lambda 和 foreign_lambda 的使用控制了 λi 的值和 λ 的附加值，然后评估在 λi 处构型的能量差异。而现在可以显式地在 fep_lambdas 中设置 λ 的值。但在此我们通过设置 calc_lambda_neighbors（见下）来自动确定这些附加值 |
| calc_lambda_neighbors = 1 | 根据 λi 计算能差的窗口间隔。例如，如果 init_lambda_state 设置为10，那么在 calc_lambda_neighbors = 1 的运行期间，将会计算 λ=9 和 λ=11 的能量差 |
| vdw_lambdas = ... | 用于 van der Waals interactions 变换的 λ 值数组 |
| coul_lambdas = ... | 用于 Coulombic (electrostatic) interactions 变换的 λ 值数组 |
| bonded_lambdas = ... | 用于 bonded interactions 变换的 λ 值数组 |
| restraint_lambdas = ... | 用于 restraints 变换的 λ 值数组 |
| mass_lambdas = ... | 用于原子质量变换的 λ 值数组（如果有分子化学性质改变的操作） |
| temperature_lambdas = ... | 温度变化的 λ 值数组（仅用于模拟退火） |
| sc-alpha = 0.5 | 在 soft-core Lennard-Jones 方程中使用的α标度因子的值。参见GROMACS手册第4.5.1节中的公式4.124 - 4.126以及接下来的三个参数 |
| sc-coul = no | 线性变换库仑相互作用。这个是默认值，但为了清楚起见写出来 |
| sc-power = 1.0 | soft-core 等式中 λ 的值 |
| sc-sigma = 0.3 | 当某一原子类型具有C6或C12参数且其值位于区间[0,σ)时，原C6或C12参数将会被替换为 σ（0.3 为典型的H原子）。用于 soft-core Lennard-Jones方程 |
| couple-moltype = Methane | 这个参数为那个拓扑结构将会从状态 A 插入到状态 B 的 [moleculetype] 的名称。注意这里给出的名称必须匹配 .top 中 [moleculetype] 的名称，而不是残基名称。在某些情况下，它们可能是相同的，但出于指导目的，作者对 [moleculetype] 和成分残基保留了不同的名称 |
| couple-lambda0 = vdw | 在插入的 [moleculetype] 和系统的其余部分之间，处于 A 态的非键相互作用的类型。vdw 表示甲烷和水之间只有范德华项是活跃的；没有溶质-溶剂的库仑相互作用 |
| couple-lambda1 = none | 在插入的 [moleculetype] 和系统其余部分之间，处于 B 态的非键相互作用的类型。none 表示状态 B 的范德华和库仑相互作用都是关闭的 |
| couple-intramol = no | 该参数设置意味着不要解耦分子内的相互作用。也就是说，λ 值只应用于 溶质-溶剂 非键相互作用，而不应用于 溶质-溶质 非键相互作用；设置 “couple-intramol = yes” 对于较大的分子来说是有用的，因为它们可能在较长的距离上发生分子内相互作用（例如多肽或其他聚合物）。否则，分子的远端部分将通过显式的对相互作用以一种不自然的强方式发生相互作用（因为 溶质-溶剂 的相互作用作为 λ 的函数被削弱了，但分子内的不会），从而产生不正确的自由能的构型。 |
| nstdhdl = 10 | ∂H/∂λ 和 ΔH 被写入 dhdl.xvg 的频率。值为100可能就足够使结果值高度相关，并且再高将使文件变得非常大。（我们在这仅作为一个教程，不需要100这么高的值）对于您自己的工作，您可能希望将此值增加到100。 |
这里使用的.mdp设置与Shirts等人使用的设置也有一些不同：
- rlist=rcoulomb=rvdw=1.2：为了使用PME, rlist 必须等于 rcoulomb，在3.3.1版本发布后的某个时候，grompp 强制执行了这个设置。在4.6版本中引入的 Verlet 截止方案允许缓存邻近列表（buffered neighbor lists）也需要 rvdw=rcoulomb。rlist 的值将在运行过程中进行调整，以节省能量，从而为交换范德瓦尔斯相互作用提供必要的缓存。转换范围(1.0-1.2 nm)与 Shirts 等人处理 溶质-溶剂 相互作用的方法一致。
- 温度耦合是通过使用 Langevin 积分器来处理的(“sd”设置指定了“随机动力学”，例如随机摩擦力)，而不是 Andersen 恒温器。
- tau_t = 1.0：这个特定的值很重要，因为当使用 Langevin 积分器时，这个值对应于ps-1中一个逆摩擦系数。1.0 的值避免了水的动力学过度阻尼。
让我们花一点时间更仔细地剖析 λ vector。例如，对于 init_lambda_state = 0，这意味着我们指定处于转换中的状态对应于在 *_lambdas 中索引为0的 vector ，如下所示:
; 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 将对应于右边的下一列（van der Waals 的 λ = 0.05），而 initial_lambda_state = 20 将指定最后一列，其中van der Waals 的 λ 值为 1.0。出于教学目的，我们只转换范德华相互作用，而不考虑与电荷、键参数、约束等其他有关的一切。请注意，在这个例子中，拓扑电荷为零，并且 .mdp 设置从未表明电荷应该存在(couple-lambda0 和 couple-lambda1都不包含’q’)，coul_lambdas 值的选择是无关紧要的，但处于某种学究式的传统，.mdp文件需要明确表示没有进行电荷转换。
Energy Minimization
job.sh 将会创建以下目录等级：
Lambda_0/
Lambda_0/EM/
Lambda_0/NVT/
Lambda_0/NPT/
Lambda_0/Production_MD/
这样，工作流中的所有步骤都在针对 init_lambda_state 每个值的单个目录中执行。作者发现这是组织作业及其输出的一种方便的方法。
脚本还假设.mdp文件也按层次结构组织在$MDP目录中，该目录在脚本中被设置为一个环境变量:
如前所述，能量最小化将使用最速下降法进行。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."
Equilibration
我们将要进行 1.NVT 与 2.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."
Data Collection
系统平衡完成后，我们将要进行动力学模拟并收集 ∂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.5h。因为有太多模拟需要进行，最好集群上并行这些程序。当所有模拟结束后，我们就可以分析数据了。
Analysis
GROMACS 中的 bar 模块使得计算 ΔGAB 变得非常简单。将运行所得所有 .xvg 文件整理然后运行 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
所得ΔGAB的值为 -9.13±0.09 kJ mol-1，或 -2.18±0.02 kcal mol-1。由于这个演示所进行的过程是甲烷的解耦，因此反向过程（将不带电荷的甲烷引入水中，从而得到该过程的实际水合能）对应ΔGAB 为2.18±0.02 kcal mol-1(假设可逆性)，与Shirts等人得到的 2.24±0.01 kcal mol-1 的值很好地吻合(根据该论文的补充材料表II)，更何况除了前面描述的其他变化，我们在范德华变换中仅使用了大约一半的λ。
从技术上讲，这是我们得到答案所需要做的一切，但 gmx 还输出了许多有用的文件(所有这些文件都是可选的)。它们的内容值得在这里探讨，因为它们提供了关于解耦过程和采样成功的大量细节。
Output file 1: bar.xvg
让我们看看 gmx ba r提供的其他输出文件。bar.xvg 绘制了λ每个区间的相对自由能差异(即，相邻的哈密顿量之间)，并应该看起来像这样(重新绘制为条形图而不是默认的直线表示，并添加了一些网格线的角度):
横轴上垂直的灰色线表示在这里进行的解耦过程中使用的λ的值。因此，每个黑条表示λ相邻值之间的自由能差。在 bar.xvg中，我们首先查看 calc_lambda_neighbors 在模拟期间做了什么。例如，在我们对 init_lambda_state = 0 的模拟中，λ = 0.05(最近邻窗口，由 calc_lambda_neighbors = 1指定)处的自由能在每一个 nstdhdl(10) 步中计算。这些“外来的” λ 计算导致λ值之间的能量重叠，这样我们就有相空间重叠和充分的采样。有兴趣的读者应该参考 gmx bar 描述中引用的Wu和Kofke的论文，以讨论这个概念以及熵值sA和sB的解释。
因此，bar.xvg 中绘制的是自由能从 λ = 0 到 λ = 1 的变化仅仅是每对相邻λ模拟的自由能变化的和。
ΔG 的值对应于按 gmx bar 打印到屏幕上的输出的前半部分。屏幕输出的第二部分包含单位转换为kJ mol-1 的 ΔG(1 kT = 2.478 kJ mol-1在298 K)。
Output file 2: barint.xvg
barint.xvg 文件将累积 ΔG 绘制为 λ 的函数。在barint.xvg， λ = 1 处的点对应于 ΔG 从 λ = 0到 λ = 1 的和，如上面的屏幕输出所示：
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 中发现累积 ΔG 在 barint.xvg 的 λ = 0 是零，而在 λ = 0.1 的值是 0.0 + 0.08 = 0.08（DG列）。下面是 barint.xvg (蓝色)与 bar.xvg (黑色，值ΔG之间的λ值)的叠合图，以说明累积ΔG的效果。通过在λ = 20 (-3.69 kT)处取ΔG的值，我们可以恢复上面所示的ΔG的值(-9.14 kJ mol-1)。这里的值(-9.14)和上面的值(-9.13)之间的差异是由于只使用两位小数点后的舍入误差造成的，而 bar 程序在上面的计算中使用了完全精度，但对所有输出只打印两位小数点后的数字。
Advanced Applications
我们计算自由能通常是为了确定配体和受体之间结合(如蛋白质)的ΔG，您需要完成受体-配体复合物的（解）结合以及他们在溶液中的表现，因为ΔG(在本例中)是受体-配体耦合（complexation，例如配体进入活性位点的过程）以及配体去溶剂化（desolvation，例如配体从水溶液中被剔除的过程）过程中的自由能变总和。
结合的 ΔG将是这两种状态的差，并可以通过使用以下(稍微简化)热力学循环计算，其中ΔG1代表ΔGsolvation，ΔG3代表ΔGcomplexation。
这类过程的一个特别详细的热力学循环可以在本文的图1中找到。
要将完全相互作用（fully-interacting）的物质（例如配体）转变为“假体”（dummy，配体构型中的一组原子中心，不参与与其周围环境的任何非键相互作用)，需要同时关闭我们所研究的溶质与其周围环境之间的范德华相互作用(如本教程所演示的)和库仑相互作用。完整的甲烷转化过程如下图所示：
在上述热力学循环中，Δ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 表示相互作用是开启的。
这些类型的转换对于计算水化/溶剂化自由能是有用的，因为这个自由能经常被用作力场参数化的目标数据。
Summary
你现在已经计算了一个简单变换的自由能变化，这个变换之前已经被高精度地计算过，简单溶质(未带电的甲烷)和溶剂(水)之间的范德华相互作用的解耦。希望本教程能够帮助您理解更复杂的系统。
Happy simulating!
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
