GROMACS Tutorial 3-SMD & Umbrella Sampling

Introduction

本教程将指导用户设置和运行一个牵引模拟（pulling simulations），以计算两个物质之间的结合能。在此我们假设你已经成功地完成了溶菌酶教程、其他一些教程或者精通基本的GROMACS模拟方法和拓扑组织。本教程特别关注于正确构建系统的方法和牵引代码本身的设置。

结合能（ΔG_bind）来自于从一系列伞形采样模拟中提取的平均力势（potential of mean force，PMF）。在伞形采样中，我们生成了一系列的初始构型，每个构型对应于一个 PMF 曲线上的位点，在这个位点上，我们所研究的分子（通常称为配体）在 受体（或受体的一部分）—配体 质心（center-of-mass，COM）距离增加时会受到其他基团，经由伞形偏置势的谐波抑制。这种约束使我们可以采集到配体在一个事先定义好的反应坐标（某两个基团之间的）上的某一个确定的区域内的构象（也即历遍某个构象空间）。这些（构象空间的）窗口必须允许配体位置有轻微的重叠，才能正确地重建PMF曲线。（品红色部分为译者添加）

伞形采样的步骤（以及本教程中使用的步骤）如下：

1. 沿单自由度（反应坐标）生成一系列构象；
2. 从步骤1的轨迹中提取与所需 COM 间距对应的帧；
3. 对每个配置进行伞形采样模拟，将其限制在与所选COM距离对应的窗口内；
4. 采用加权直方图分析法（WHAM）提取 PMF，计算 ΔG_bind

本教程假设读者使用的是GROMACS 2018或更高版本。作者最初的工作（该 workflow 的来源）是在版本4.0.5中完成的，但原则上可以应用于4.0.x或4.5.x系列中的任何版本。牵引代码在3.3.3版本之后被完全重写，在5.1版本中再次重写，所以使用GROMACS 2018及以上版本很重要。

在开始前，请务必阅读相关文章以便于理解此类模拟的大体思路。简单地将本教程使用的设置应用于其他系统是不合适的，一些需要具体考虑的因素将在后面讨论。本教程的目的是为用户提供伞形采样的理论理解和一个直接从文献中拿来的实际例子。

The Topology

伞形采样模拟生成分子拓扑和其他模拟一样：获取结构的坐标文件，并使用 pdb2gmx 生成拓扑。有些系统需要特殊考虑（例如，蛋白-配体复合物，膜蛋白等）。对于 蛋白—配体 系统，请参考本教程，对于膜蛋白，我推荐这个教程。伞形采样的原理很容易扩展到这些系统中。

本教程使用的系统是单个肽从 Aβ42 原纤维的生长端的解离，基于作者最近发表的模拟。在这些模拟中使用了野生型Aβ42原纤维的结构文件（每个链的 N 端乙酰化），可以在这里找到（原始文件PDB号为2BEG）。

通过pdb2gmx得到该结构：

gmx pdb2gmx -f 2BEG_model1_capped.pdb -ignh -ter -o complex.gro

选择 GROMOS96 53A6 参数集，SPC 为水模型，N 端为 None，C端为添加 COO^-。在topol_Protein_chain_B.itp 的末尾添加：

#ifdef POSRES_B
#include "posre_Protein_chain_B.itp"
#endif

在后面的牵引模拟中，我们将使用 B 链 作为固定的参考，因此只需要对这条链进行特殊的位置约束，其他链都不需要（ 在这里可以将 B 链理解为受体，A 链理解为配体，我们要做的就是固定受体来拉远配体，以达到模拟解离的效果）。

Defining the Unit Cell

为牵引模拟定义胞元并无特殊之处。然而，有一点需要注意：必须在牵引方向上留出足够的空间，以便在不与系统周期性图像相互作用的情况下进行连续拉动。也就是说，必须始终满足最小图像约定（the minimum image convention），并且拉动距离必须始终小于进行拉动的盒向量长度的一半。你可能会问，为什么？

GROMACS计算距离时同时考虑周期性。如果你有一个10纳米的盒子，你设定了超过5.0纳米的牵引距离，牵引的距离将参考周期性长度而（在实际模拟中）小于5.0纳米！这一事实将显著影响模拟结果，因为你认为你正在拉的距离不是实际计算的距离。

我们将在一个12.0纳米的盒子中拉动总距离为5.0纳米的物体，以避免上述现象。原纤维的质心将放置在尺寸为6.560 x 4.362 x 12的盒子中（3.280,2.181,2.4775）的位置。使用editconf将原纤维放置在此位置：

gmx editconf -f complex.gro -o newbox.gro -center 3.280 2.181 2.4775 -box 6.560 4.362 12

您可以使用VMD查看原纤维在盒子中的位置：加载结构并打开Tk控制台，输入以下内容(注意 > 表示Tk提示符，而不是你实际输入的东西)：

> pbc box

他理应长这样：

Solvation and Ions

首先，我们将加入溶剂水：

gmx solvate -cp newbox.gro -cs spc216.gro -o solv.gro -p topol.top

接下来，我们将使用 genion 添加离子，使用这个.mdp文件。添加100 mM NaCl并中和多余电荷：

gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral -conc 0.1

Energy Minimization and Equilibration

能量最小化和平衡的步骤也同其他系统一样。在这里，我们将执行EM，然后NPT平衡。最小化的 .mdp 文件可以在这里找到，NPT平衡的可以在这里找到。

像往常一样调用 grompp 和 mdrun：

gmx grompp -f minim.mdp -c solv_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em

gmx grompp -f npt.mdp -c em.gro -p topol.top -r em.gro -o npt.tpr
gmx mdrun -deffnm npt

Generating Configurations

为了进行伞形采样，必须沿着反应坐标 ζ 生成一系列构型。其中一些配置将作为伞形采样窗口的起始配置，它们在相互独立的模拟过程中进行。下图解释了这样做的原理：上面的图像显示了我们现在将运行的牵引模拟，该模拟的目的是沿着反应坐标生成一系列构型。这些配置在模拟完成后将被提取（顶部和中间图像之间的虚线箭头）。中间的图像对应于在每个采样窗口内进行的独立模拟，自由肽的质心在其特定窗口内受到伞状偏置电位的限制。底部的图像显示了不同配置下模拟得到的直方图（不过是理想的结果），各个相邻配置的窗口重叠，以便以后可以从这些模拟中导出连续的能量函数。

对于本教程，反应坐标是z轴。为了产生这些结构，我们必须把A肽从原纤维中拉出来。我们将使用500ps的MD，每1ps保存一次快照，以获得配置的合理分布，这个设置是基于试错法建立的。在其他系统中，可能需要更频繁或更节俭地保存配置。这一步的整体思路是沿着反应坐标捕获足够多的构型，以根据肽A和肽B之间的质心距离来获得伞样窗口的规则间距。

此操作的.mdp文件可以在这里找到（在GROMACS 5.x 版本往上对 pull_code 的名称进行了较大的调整（详细可参考GROMACS手册 7.3.21 部分），上文提供的 .mdp 文件已经无法用于最新版本，这里译者进行了一定的修改，可以在这里下载最新的）。下面简要说明一下所使用的牵引选项：

; Pull code
pull                    = yes
pull_ncoords            = 1         ; only one reaction coordinate
pull_ngroups            = 2         ; two groups defining one reaction coordinate
pull_group1_name        = Chain_A
pull_group2_name        = Chain_B
pull_coord1_type        = umbrella  ; harmonic potential
pull_coord1_geometry    = distance  ; simple distance increase
pull_coord1_dim         = N N Y     ; pull along z
pull_coord1_groups      = 1 2       ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
pull_coord1_start       = yes       ; define initial COM distance > 0
pull_coord1_rate        = 0.01      ; 0.01 nm per ps = 10 nm per ns
pull_coord1_k           = 1000      ; kJ mol^-1 nm^-2

• pull = yes：告诉grompp读取 COM 牵引设置。如果省略这一句，所有的牵引设置将被忽略；

• pull_ncord = 1：定义系统中存在的反应坐标的数量。为了多维自由能计算的目的，这个设置允许同时设置多个偏差，但本教程中只存在一个反应坐标。

• pull_ngroups = 2：有两个组定义了反应坐标的端点。反应坐标是这两个基团 COM 之间的距离。

• pull_group1_name = Chain_A和pull_group2_name = Chain_B：将被应用于牵引代码的组名，从索引文件中读取。由于在这个系统中指定了两个组（上一行中的 pull_ngroups = 2），因此必须提供两个名称。

• pull_coord1_type = umbrella：使用谐波势进行牵引。重要提示：此声明不是伞形采样。谐波势允许力根据肽A与肽b相互作用的性质而变化，也就是说，力会逐渐增加，直到某些关键的相互作用被打破。详情请参阅我们的论文。为了生成伞形采样的初始配置，您可以使用不同的设置组合(pull_coord1_type和pull_coord1_geometry)，但对实际伞形采样（下一步）时，你必须使用pull_coord1_type = umbrella。非常重要的是，不要使用非常快的牵引速率或非常强的力常数，这可能使系统各要素严重扭曲。请参阅论文(特别是支持信息)，了解我们如何选择验证所使用的牵引率。还要注意，“coord1”作为设置命的关键字之一，意味着这是对于第一个反应坐标的设置。对于每个反应坐标（pull_ncord），带有该关键字的所有设置都是必须的。如果有第二个反应坐标，设置将命名为“coord2”，第三个设置将命名为“coord3”，依此类推。

• pull_coord1_geometry = distance：参见.mdp文件中的注释;您也可以使用位置、方向或方向周期（position, direction, or direction-periodic），但必须更改其他牵引参数。

• pull_coord1_dim = N N Y：我们只在z轴上应用了偏差限值。因此，x和y被设置为“否”(N)， z被设置为“是”(Y)。

• pull_coord1_groups：反应坐标连接了组1和组2。

• pull_coord1_start = yes：初始 COM 距离是第一帧的参考距离。下一阶段将进行更多的设置，那时将进行实际的伞采样模拟

• pull_coord1_rate = 0.01：连接在牵引组上的虚拟弹簧被拉长的速率。这个速度（在模拟过程中）是恒定的。

• pull_coord1_k = 1000：用于牵引的弹簧上的力常数。

我们要使用先前添加到 topol_B 中的 #ifdef POSRES_B 语句。通过抑制原纤维的肽B，我们能够更容易地将肽A拉离。由于A链和B链之间存在广泛的非共价相互作用，如果我们不限制B链，我们就只能简单地将整个复合物沿着模拟盒子拖着走，这并不能实现我们的目的。

我们需要为这个牵引模拟定义一些自定义索引组。使用 make_ndx：

gmx make_ndx -f npt.gro
...
 > r 1-27
 > name 19 Chain_A
 > r 28-54
 > name 20 Chain_B
 > q

运行MD：

gmx grompp -f md_pull.mdp -c npt.gro -p topol.top -r npt.gro -n index.ndx -t npt.cpt -o pull.tpr

gmx mdrun -deffnm pull -pf pullf.xvg -px pullx.xvg

整个模拟结果理应如下（原教程是动态的gif）：

A 链在一段时间内没有位移，因为力在假想的弹簧上逐渐增加，直到它足以克服原纤维结构内的恢复力。大约在 SMD 模拟进行到一半的时候，链 A 开始从链 B 分离，这个过程也可以在弹簧上的力随着时间推移的中看到，这部分数据存储在pullf.xvg中：

为了准备单独的伞形采样窗口，我们需要从 SMD 轨迹中提取有用的构象。作者发现最简单的方法如下:

1. 定义窗口间距(一般为0.1 - 0.2 nm)；
2. 从刚才产生的牵引轨迹中提取出所有的帧；
3. 测量两组之间每一帧的COM距离来定义反应坐标；
4. 使用选定的构象作为伞形采样的输入。

要从你的轨迹中提取帧（traj.xtc），使用 trjconv (当提示时选择 0 来保存整个系统)：

gmx trjconv -s pull.tpr -f pull.xtc -o conf.gro -sep

一系列坐标文件将产生（conf0. gro，conf1. gro…501个！实际上我们可以根据实际需求，添加 -skip 标签来跳过一部分，例如在此教程中使用gmx trjconv -s pull.tpr -f pull.xtc -o conf.gro -sep -skip 2就只会产生251个坐标文件），对应于保存在连续牵引模拟的每一帧。为了对生成的所有这些帧迭代调用 gmx distance，作者编写了一个 Bash 脚本来处理这个任务。它将打印一个名为“summary_distance.dat”的文件，其中包含了这些信息。可以在这里找到脚本。

改脚本内容为：

#!/bin/bash

# compute distances 
for (( i=0; i<1001; i++ ))
do
    gmx distance -s ../smd_o_000.tpr -f t1_conf${i}.gro -n ../index.ndx -select 'com of group 29 plus com of group 30' -oall dist${i}.xvg 
done

# compile summary
touch summary_distances.dat
for (( i=0; i<1001; i++ ))
do
    d=`tail -n 1 dist${i}.xvg | awk '{print $2}'`
    echo "${i} ${d}" >> summary_distances.dat
    rm dist${i}.xvg
done

exit;

在应用于其他体系之前，我们需要进行一些调整：其中两个 for 语句的 1001 改成之前gmx trjconv输出出的 conf?.gro 中?的最大值；gmx distance那一行命令中，-s 后为你做 SMD 的 tpr 文件名，-f 是之前 gmx trjconv 输出出的conf?.gro 文件名（${i}代表序号），-n 是 SMD 用到的索引文件名，-select 中 group ? 的 ? 换成你做 SMD 的两个组在索引文件中的序号。

查看 summary_distance.dat 的内容，可以看到链 A 和链 B 之间的 COM 距离随时间的变化。根据所需的间距，记录伞形采样时使用的配置。也就是说，如果你想要0.2 nm的间距，你可能会在summary_distance.dat中找到以下行:

6     0.500
...
160    0.704

然后使用 conf6.gro 和 conf160.gro 作为两个相邻的伞取样窗口的起始配置。在继续之前，请记下您希望使用的所有配置。对于本教程来说，识别具有0.2 nm间距的构型就足够了，尽管在最原始的工作中，为了更彻底地定义能量最小时的自由能，使用了不同的间距。

Umbrella Sampling

在本例中，我们将以大约 0.2 nm 的间距沿 z 轴从 0.5 - 5.0 nm 采样。下面的示例命令可能是也可能不是正确的（帧数可能不同），但作为如何在单独的坐标文件上运行grompp以生成所有23个输入的示例已经很足够了（还要注意，在大约4.5 nm范围内获得0.2 nm间距所需的窗口数量为23；在我们最初的工作中，使用了31个非对称窗口）。

（原文中本段出现在最后，为了防止有人操之过急，放到了前面。请先阅读该段落之后的内容，了解基本思路后再来下载这个脚本。）Mike Harms贡献了一个Python脚本，它可以自动化这个过程，提取坐标文件并设置grompp和mdrun命令来简化这个过程。您可以在这里下载他的脚本和相关文档。

python setupUmbrella.py summary_distances.dat 0.2 run-umbrella.sh &> caught-output.txt

在确定了采样窗口的初始配置之后，我们现在可以准备伞形采样模拟了。我们将需要生成大量的输入文件，以便进行每一个必要的模拟。例如，由于我们沿着反应坐标确定了23种构型，这意味着我们将需要23种不同的输入文件来进行23次独立的模拟。首先使用这个.mdp文件在每个窗口中运行一个简短的NPT平衡（ 牵引参数的设置与下文将用到的 .mdp相似，具体解释也见下）。

第一个grompp(输入conf*.gro，根据不同的SMD模拟结果，文件名可能不同)：

gmx grompp -f npt_umbrella.mdp -c conf6.gro -p topol.top -r conf6.gro -n index.ndx -o npt0.tpr
...
gmx grompp -f npt_umbrella.mdp -c conf449.gro -p topol.top -r conf449.gro -n index.ndx -o npt22.tpr

随后 mdrun：

gmx mdrun -deffnm npt0
...
gmx mdrun -deffnm npt22

为了开始伞式采样，您只需调用 grompp 来处理每个现在已经平衡的配置的 .mdp 文件。许多牵引参数与 SMD 过程中的参数相同，但有一个例外，即 pull_coord1_rate，它现在已被设置为0。我们不想让构型沿着反应坐标移动；相反，我们想把它限制在构型空间的一个定义窗口内。设置pull_coord1_start = yes 意味着初始COM距离是参考距离，我们不必为每个配置单独定义引用(pull_coord1_init)。md 所需 .mdp 可在这里下载。

gmx grompp -f md_umbrella.mdp -c npt0.gro -t npt0.cpt -p topol.top -r npt0.gro -n index.ndx -o umbrella0.tpr
...
gmx grompp -f md_umbrella.mdp -c npt22.gro -t npt22.cpt -p topol.top -r npt22.gro -n index.ndx -o umbrella22.tpr

现在，应该将每个输入文件交给 mdrun，以进行实际的数据收集模拟。一旦所有的模拟都完成了，就可以进行数据分析了。

gmx mdrun -deffnm umbrella0
...
gmx mdrun -deffnm umbrella22

Data Analysis

对伞取样模拟进行的最常见分析是提取平均力势（potential of mean force，PMF），这将得到 结合/解离 定过程的 ΔG。ΔG 的值就是PMF曲线的最大值和最小值的差值，只要 PMF 的值在大 COM 距离下收敛到一个稳定的值。提取PMF的常用方法是加权直方图分析法（WHAM），它作为WHAM 模块包含在 GROMACS 中。wham 的输入由两个文件组成，一个列出每个窗口的 .tpr 文件的名称，另一个列出每个窗口的 pullf. xvg 或 pullx. xvg的名称。例如，一个简单的 tpr-files.dat 可能包含：

umbrella0.tpr
umbrella1.tpr
...
umbrella22.tpr

类似地，对于任意一个 pullf. xvg 或 pullx. xvg文件中的 pullf-files.dat 或 pullx-files.dat 文件列表，文件必须具有唯一的名称（例如，umbrella0_pulf.xvg、umbrella1_pullf.xvg 等)，否则 wham 将报错。按与 tpr-files.dat 相同的方式列出 pullf-files.dat 中的文件:

umbrella0_pullf.xvg
umbrella1_pullf.xvg
...
umbrella22_pullf.xvg

然后运行 gmx wham：

gmx wham -it tpr-files.dat -if pullf-files.dat -o -hist -unit kCal

然后，wham模块将打开每个 umbrella*. tpr 和 umbrella*_pullf. xvg 文件，并对它们运行 WHAM 分析。-unit kCal选项表示输出的单位是 kcal mol^-1，但也可以得到 kJ mol^-1 或 kBT。得到的 PMF和相应的伞状直方图看起来像这样:

这里得到的结果对应于ΔG = 46 kcal mol^-1，与我们公布的 50.5 kcal mol^-1 相当。由于窗口间距不相同和启动配置不同，预计会有一些差异。特别要注意的是，ξ = 0.8 nm附近的PMF轮廓有一个缺陷，这反映了直方图中明显缺乏抽样。ξ 值所在区域为高能量状态，由于没有以该区域为中心的窗口，导致采样不足。为了解决这个问题，进行了一个额外的模拟，以ξ = 0.8 nm为中心的窗口。由于所有的伞取样窗口模拟都是独立的，其他的不需要重新运行，只需要包括新 umbrella.tpr 和 umbrella_pullf.xvg 文件并重新运行 WHAM 分析。

Summary

现在，通过沿反应坐标生成一系列构型，运行偏态模拟，并提取PMF，您已经成功地进行了伞式采样模拟。这里提供的 .mdp 文件仅作为示例，不应该被认为广泛适用于所有系统。回顾文献和GROMACS手册，以调整这些文件的效率和准确性的目的。

Happy simulating!