PLUMED MetaDynamics Simulation with GROMACS

该教程原文链接为 https://www.plumed.org/doc-v2.4/user-doc/html/belfast-6.html

本教程的目的是向用户介绍使用 PLUMED 运行元动力学模拟。我们将在真空中建立丙氨酸二肽的简单模拟，分析输出，并从模拟中估计自由能。我们还将学习如何运行回火元动力学模拟（well-tempered metadynamics simulation），并检测与集体变量选择不当有关的问题。“😇”内的内容为译者添加。

Summary of theory

在元动力学中，在几个选定的自由度 s(q) 的空间中构建了一个外部历史相关的偏置势，通常称为集体变量（CV）。该势为轨迹在CV空间中的高斯沉积的总和：$$V(\vec{s},t)=\sum _{k\tau <t}W(k\tau )\exp \left(-\sum _{i=1}^d\frac{(s_i-s_i(q(k\tau )){)}^2}{2{\sigma }_i^2}\right).$$其中，τ 是高斯沉积步长，σ_i 是第i个CV的高斯宽度，W(kτ) 是高斯高度。元动力学偏置势的作用是将系统从局部极小值推到相空间的新区域。此外，在长时间限制中，偏置势收敛为减去作为CV的函数的自由能：$$V(\vec{s},t\to \infty )=-F(\vec{s})+C.$$在标准元动力学中，在整个模拟过程中都会添加恒定高度的高斯。结果是系统最终被推动并探索高自由能区域，并且根据偏置势计算的自由能的估计值围绕真实值振荡。在回火元动力学中，高斯高度随着模拟时间的推移而降低，具体如下：$$W(k\tau )=W_0\exp \left(-\frac{V(\vec{s}(q(k\tau )),k\tau )}{k_B\Delta T}\right),$$其中 W₀ 是初始高斯高度，ΔT 是具有温度维度的输入参数，k_B 是玻尔兹曼常数。通过对高斯高度的重新缩放，偏置电在长时间限制内平滑收敛，但它不能完全补偿潜在的自由能：$$V(\vec{s},t\to \infty )=-\frac{\Delta T}{T+\Delta T}F(\vec{s})+C.$$其中 T 是系统的温度。因此，在长时间限制中，CV 在高于系统温度 T 的温度 T+ΔT 下对系综进行采样。可以选择参数 ΔT 来调节自由能探索的程度：ΔT=0 对应标准分子动力学，ΔT → ∞ 对应标准元动力学。在回火元动力学文献和 PLUMED 中，您经常会遇到术语“偏因子”，即CV的温度（T+ΔT）和系统温度（T）之间的比率：$$\gamma =\frac{T+\Delta T}{T}.$$因此，应该仔细选择偏置因子，以便在模拟的时间尺度上有效地跨越相关的自由能垒。

其他信息可以在关于元动力学的几篇综述论文中找到

完成本教程后，你将知道如何：

• 使用PLUMED运行元动力学模拟 😇 Ex1.
• 分析模拟的输出 😇 Ex1.
• 重新启动元动力学模拟 😇 Ex2. 续跑元动力学
• 从元动力学模拟计算自由能 😇 Ex3.
• 使用PLUMED运行一个温和的元动力学模拟 😇 Ex4.
• 检测集体变量选择的问题 😇 Ex5.

此项目的tarball（可点击进行下载）包含以下目录：

• TOPO：它包含在真空中模拟丙氨酸二肽的gromacs拓扑和配置文件
• Exercise_1：使用2个CV、二面体phi和psi运行元动力学模拟，并分析输出
• Exercise_2：重新启动元动力学模拟
• Exercise_3：从元动力学模拟计算自由能并监测收敛性
• Exercise_4：使用2个CV、二面体phi和psi进行回火的元动力学模拟
• Exercise_5：用1个CV，二面体psi进行回火的元动力学模拟

Instructions

本教程使用AMBER99SB-ILDN全原子力场模拟丙氨酸二肽体系。

Exercise 1. Setup and run a metadynamics simulation

Standard metadynamics simulation

在本练习中，我们将在真空中对丙氨酸二肽进行元动力学模拟，使用两个主链二面角 phi 和 psi 作为 CV。为了运行此模拟，我们需要准备 PLUMED 输入文件（PLUMED.dat），如下所示：

# 定义 Phi 和 Psi 二面角设置
phi: TORSION ATOMS=5,7,9,15
psi: TORSION ATOMS=7,9,15,17
#
#激活使用 phi 和 psi 的元动力学
#每 500 个时间步长进行一次高斯沉积，
#高度等于 1.2 kJ/mol，
#两个CV的宽度（标准差）均为 0.35 rad
#
metad: METAD ARG=phi,psi PACE=500 HEIGHT=1.2 SIGMA=0.35,0.35 FILE=HILLS 

# 监测这两个变量和元动力学的偏置势
PRINT STRIDE=10 ARG=phi,psi,metad.bias FILE=COLVAR

METAD的语法很简单。该指令后面跟着的选项ARG，后面写明元动力学势能将作用的CV的标签。选项PACE为进行高斯沉积的时间步长，而关键字HEIGHT以kJ/mol指定高斯沉积的高度。对于每个CV，必须使用关键字SIGMA指定高斯的宽度。高斯数将被写入由关键字FILE指示的文件中。

准备好 PLUMED 输入文件后，您需要使用类似下述命令运行 Gromacs：

gmx_mpi mdrun -v -deffnm meta_1 -plumed

在元动力学模拟过程中，PLUMED 将创建两个文件，分别命名为COLVAR和HILLS。COLVAR文件包含PRINT命令指定的所有信息，在本例中为模拟每10步的CV值，以及元动力学偏置势的值。HILLS文件包含沿模拟放置的高斯图的列表。如果我们看一下这个文件的标题，我们可以找到有关其内容的相关信息：

#! FIELDS time phi psi sigma_phi sigma_psi height biasf
#! SET multivariate false
#! SET min_phi -pi
#! SET max_phi pi
#! SET min_psi -pi
#! SET max_psi pi

从FIELDS开始的行告诉我们在HILLS文件的各个列中显示了什么：模拟时间、phi和psi的值、phi和psi的高斯宽度、高斯高度和偏移因子。偏移因子的大小仅与回火元动力学模拟相关（请参见 Exercise 4. Setup and run a well-tempered metadynamics simulation, part I），它在标准元动力学模拟中等于1。我们稍后将使用HILLS文件来计算元动力学模拟的自由能，并评估其收敛性。目前，我们可以绘制模拟过程中CV的变化：

通过上图我们可以看到，系统初始为丙氨酸二肽的两种亚稳态之一。一段时间后（t=0.3ns），系统被元动力学偏置势推动，访问了另一个局部最小值。随着模拟的继续，偏置势填充了底层的自由能景观，系统能够在整个相空间中扩散。

Standard metadynamics simulation with grid

如果我们使用上述 PLUMED 输入文件，元动力学模拟的费用会随着模拟的长度而增加，因为在每一步都必须评估越来越多的高斯值。为了避免这个问题，您可以将偏差存储在网格上。为了使用网格，我们必须在METAD指令的行中添加一些附加信息，如下所示：

# 定义 Phi 和 Psi 二面角设置
phi: TORSION ATOMS=5,7,9,15
psi: TORSION ATOMS=7,9,15,17
#
# 激活使用 phi 和 psi 的元动力学
# 每 500 个时间步长进行一次高斯沉积，
# 高度等于 1.2 kJ/mol，
# 两个CV的宽度（标准差）均为 0.35 rad
# 偏置势将会储存到网格中
# 对于两个CV而言，网格窗口大小均为 0.1 rad 
# 对于两个CV而言，网格边界均为 -pi and pi
#
METAD ...
LABEL=metad
ARG=phi,psi 
PACE=500
HEIGHT=1.2
SIGMA=0.35,0.35
FILE=HILLS
GRID_MIN=-pi,-pi
GRID_MAX=pi,pi
GRID_SPACING=0.1,0.1 
... METAD

# 监测这两个变量和元动力学的偏置势
PRINT STRIDE=10 ARG=phi,psi,metad.bias FILE=COLVAR

偏置势将存储在网格上，其边界由关键字GRID_MIN和GRID_MAX指定。请注意，您应该为每个 CV提供储存窗口的数量（GRID_BIN）或网格间距（GRID_SPACING）。如果您同时给出了两个选项，则 PLUMED 将对每个维度使用最保守的选择（使窗口数量最大的情况）。如果您没有提供任何关于窗口大小的信息（既没有GRID_BIN也没有GRID_SPACING），并且如果高斯宽度是固定的，PLUMED 将使用高斯宽度的1/5作为网格间距。对于大多数应用程序来说，此默认选择应该是合理的（😇指高斯宽度的1/5😇）。

Exercise 2. Restart a metadynamics simulation

如果我们尝试在已经存在COLVAR和HILLS文件的目录中使用上面的脚本再次运行元动力学模拟，PLUMED将创建旧文件的备份副本，并运行新的模拟。相反，如果我们想重新启动以前的模拟，我们必须将关键字RESTART添加到 PLUMED 输入文件中，如下所示：

# restart previous simulation
RESTART

# set up two variables for Phi and Psi dihedral angles 
phi: TORSION ATOMS=5,7,9,15
psi: TORSION ATOMS=7,9,15,17
#
# Activate metadynamics in phi and psi
# depositing a Gaussian every 500 time steps,
# with height equal to 1.2 kJoule/mol,
# and width 0.35 rad for both CVs. 
#
metad: METAD ARG=phi,psi PACE=500 HEIGHT=1.2 SIGMA=0.35,0.35 FILE=HILLS 

# monitor the two variables and the metadynamics bias potential
PRINT STRIDE=10 ARG=phi,psi,metad.bias FILE=COLVAR

😇
同时 Gromacs 也进行续跑：

gmx convert-tpr -s meta_1.tpr -extend 10000 -o meta_1.tpr

gmx mdrun -v -deffnm meta_1 -cpi meta_1.cpt -plumed

😇

通过这种方式，PLUMED 将从HILLS文件中读取旧的高斯值，并将新信息附加到COLVAR和HILLS两个文件中。

Exercise 3. Calculate free-energies and monitor convergence

Two-dimensional free energy

可以直接根据元动力学偏置势来估计自由能作为元动力学CV的函数。为了做到这一点，应使用sum_hills程序对模拟过程中存储并存储在HILLS文件中的高斯进行求和。为了计算作为phi和psi函数的二维自由能，使用以下命令行就足够了：

plumed sum_hills --hills HILLS

上面的命令生成一个名为fes.dat的文件，内含规则网格上计算的作为phi和psi函数的自由能表面。通过使用sum_hills的以下选项，可以修改自由能量文件的默认名称，以及网格的边界和窗口大小：

--outfile # 设置输出文件路径/名称
--min # 网格的最小边界
--max # 网格的最大边界
--bin # 网格的窗口数量
--spacing #网格间距，与窗口数量等效

😇对于输出文件fes.dat，译者发现qtgrace等画图软件似乎没法很好的绘制三维形貌图。所以这里给一个MATLAB绘制的方法：1. 将fes.dat中的前三列导入MATLAB中；2. 格栅化[X,Y,Z]=griddata(CV1,CV2,FES,linspace(min(CV1),max(CV1))',linspace(min(CV2),max(CV2)))；3. 绘制三维图mesh(X,Y,Z) 或等高线图contour(X,Y,Z,20)😇

One-dimensional free energy

还可以从二维元动力学模拟中计算一维自由能。例如，如果对自由能单独作为 phi 二面角的函数感兴趣，则应使用以下命令行：

plumed sum_hills --hills HILLS --idw phi --kt 2.5

Monitor convergence

为了评估元动力学模拟的收敛性，可以计算自由能的估计值作为模拟时间的函数。在收敛时，除了恒定的偏移之外，重构的轮廓应该是相似的。选项–stride应用于估计每N个高斯沉积的自由能，选项–mintozero可用于通过将全局最小值设置为零来对齐轮廓。使用以下命令行，可得下图左：

plumed sum_hills --hills HILLS --idw phi --kt 2.5 --stride 500 --mintozero

为了更定量地评估模拟的收敛性，我们可以计算一维自由能中沿 phi 的两个局部极小值之间的自由能差，作为模拟时间的函数。我们可以使用bash脚本analyze_FES.sh来积分两个盆地中的多个自由能剖面，这两个盆地由 phi 空间中的以下区间定义：basin A, -3<phi<-1, basin B, 0.5<phi<1.5；运行下行命令，可得上图右：

./analyze_FES.sh NFES -3.0 -1.0 0.5 1.5 KBT 

其中，NFES是由sum_hills的步长选项生成的轮廓数（模拟不同时间的自由能量估计），KBT是以能量单位表示的温度（在这种情况下，KBT=2.5）。

以上所有分析及 Ex1. 对CV空间中扩散行为的观察表明，模拟是收敛的。

Exercise 4. Setup and run a well-tempered metadynamics simulation, Part I

在本练习中，我们将使用两个主链二面角 phi 和 psi 作为CV，在真空中对丙氨酸二肽进行良好的元动力学模拟。要使用回火元动力学，我们需要在METAD的行中添加两个关键字，该行指定模拟的偏置因子和温度。对于第一个例子，我们将尝试一个等于6的偏置因子。以下是 PLUMED 输入文件的内容：

# set up two variables for Phi and Psi dihedral angles 
phi: TORSION ATOMS=5,7,9,15
psi: TORSION ATOMS=7,9,15,17
#
# Activate metadynamics in phi and psi
# depositing a Gaussian every 500 time steps,
# with height equal to 1.2 kJoule/mol,
# and width 0.35 rad for both CVs.
# 激活回火动力学
# 设置偏置因子为6
#
metad: METAD ARG=phi,psi PACE=500 HEIGHT=1.2 SIGMA=0.35,0.35 FILE=HILLS BIASFACTOR=6.0 TEMP=300.0

# monitor the two variables and the metadynamics bias potential
PRINT STRIDE=10 ARG=phi,psi,metad.bias FILE=COLVAR

在使用上述指令运行模拟之后，我们可以查看HILLS文件（如下）。与标准元动力学不同，HILLS文件的最后两列报告了更多有用的信息。最后一列包含所使用的偏置因子的值，而倒数第二列包含高斯高度，高斯在模拟过程中会按照回火方法重新缩放。

#! FIELDS time phi psi sigma_phi sigma_psi height biasf
#! SET multivariate false
#! SET min_phi -pi
#! SET max_phi pi
#! SET min_psi -pi
#! SET max_psi pi
      1.0000     -1.3100    0.0525          0.35            0.35      1.4400      6
      2.0000     -1.4054    1.9742          0.35            0.35      1.4400      6
      3.0000     -1.9997    2.5177          0.35            0.35      1.4302      6
      4.0000     -2.2256    2.1929          0.35            0.35      1.3622      6

如果我们仔细查看高度列，我们会注意到在开始时报告的值高于输入文件中指定的初始高度，该初始高度应为1.2 kJ/mol。事实上，这一列报道了由预因子重新缩放的高斯高度，该预因子在回火元动力学中将偏置势与自由能联系起来。这样，当我们使用sum_hills时，沉积的高斯数之和将直接提供自由能，而无需进一步重新缩放。

我们可以绘制CV和高斯高度随时间的演变（偏置因子为 6）：

系统被初始化为局部最小值之一，在该局部最小值开始累积偏置。随着模拟的进行和添加的偏置的增加，高斯高度逐渐减小。一段时间后（t=0.8ns），系统能够逃脱局部最小值，并探索相空间的新区域。一旦发生这种情况，高斯高度将恢复到初始值，并开始再次减小。随着时间的推移，当系统在整个 CV 空间中扩散时，高斯高度变得越来越小。

我们现在可以尝试不同的偏置因子，并查看模拟效果。如果我们选择一个等于1.5的偏置因子，我们可以注意到高斯高度随着模拟时间的推移会更快地减小，正如回火方法所预期的那样。如下图，该偏差因子不足以允许系统在该模拟的时间尺度上逃离初始局部最小值：

按照前面例子中为标准元动力学描述的程序，我们可以估计自由能作为时间的函数，并使用analyze_FES.sh脚本监测模拟的收敛性。我们将对偏置因子设置为 6.0 的模拟执行此操作。在这种情况下，我们会注意到（如下图），在标准元动力学中观察到的振荡在这里是阻尼的，并且偏置势更平滑地收敛到潜在的自由能景观，前提是偏置因子足够高，可以跨越所研究系统的自由能垒。

Exercise 5. Setup and run a well-tempered metadynamics simulation, Part II

在本练习中，我们将研究在选择元动力学 CV 时忽略相关自由度的影响。我们将使用以下 PLUMED 输入文件，以 psi 二面体单独作为 CV 来运行一个回火元动力学模拟：

# set up two variables for Phi and Psi dihedral angles 
phi: TORSION ATOMS=5,7,9,15
psi: TORSION ATOMS=7,9,15,17
#
# Activate metadynamics in psi
# depositing a Gaussian every 500 time steps,
# with height equal to 1.2 kJoule/mol,
# and width 0.35 rad.
# Well-tempered metadynamics is activated,
# and the biasfactor is set to 10.0
#
metad: METAD ARG=psi PACE=500 HEIGHT=1.2 SIGMA=0.35 FILE=HILLS BIASFACTOR=10.0 TEMP=300.0

# monitor the two variables and the metadynamics bias potential
PRINT STRIDE=10 ARG=phi,psi,metad.bias FILE=COLVAR

让我们看看HILLS文件，特别是 CV psi 和高斯高度随时间的变化。下图可以看出当高斯高度已经很小时，CV-psi 的扩散行为很好。这种情况一直持续到 t=3ns，此时 CV 似乎在 CV 空间的一个小区域中停留了一段时间。这种行为是 CV集合 中不包含慢变量（slow variable）的典型情况。当在这个隐藏的自由度中发生某些事情时，偏置CV通常无法访问以前访问过的相空间的任何区域。为了理解这种行为，我们需要可视化存储在COLVAR文件中的 phi 和 psi 随时间的变化。

从上面的图中可以清楚地看出，在 t=3ns 左右发生的事情是被忽视的慢变量 phi 从一个自由能池跳到另一个自由能量池。phi 的动力学不受任何势的影响，因此我们需要平衡这一自由度，即在模拟收敛之前，观察两个盆地的多次转换。或者，我们可以将 phi 添加到 CV 集。该示例演示了如何确认回火的元动力学模拟的收敛性——有必要但不足以观察到：1）高度非常小的高斯，2）CV空间中的扩散行为（如本示例的前3ns）。

我们应该做的是从不同的初始构象开始多次重复模拟。如果在所有模拟中，当高斯高度很小时，我们在偏置 CV 中观察到扩散行为，并且我们获得了非常相似的自由能表面，那么我们可以非常确信我们的模拟收敛到了正确的值。如果不是这样，手动检查运行结果（轨迹）可以帮助我们确定缺少的慢速自由度，以便添加到有偏差的CV集中。