GROMACS Tutorial 1-Lysozyme in Water-translated with notes

一些GROMACS基础知识

略

pdb2gmx

我们必须下载我们将要研究的蛋白质结构文件。在本教程中，我们将使用鸡蛋蛋白溶菌酶(PDB:1AKI)，去RCSB网站下载其PDB。
一旦您下载了结构，您就可以使用可视化程序(如VMD、Chimera、PyMOL等)将结构可视化。当你可视化结构之后，你就会想要去掉结晶水。要删除水分子(PDB文件中的“HOH”)，可以使用纯文本编辑器，如vi、emacs (Linux/Mac)或记事本(Windows)。不要使用文字处理软件!或者，你可以使用grep非常容易地删除这些行:

grep -v HOH 1aki.pdb > 1AKI_clean.pdb

注意，这种方法并不是普遍适用的(例如形成紧密结合或其他功能活性位点水分子，grep命令可参见此处)。在本篇教程中我们不需要结晶水。

记得检查您.pdb文件中 MISSING 注释下列出的条目，这些条目表明原子或整个残基不存在于晶体结构中。终端区域（terminal regions）的缺失可能不会对模拟造成影响，而不完整的内部序列或任何缺少原子的氨基酸残基都会导致 pdb2gmx 指令的失败。这些缺失的原子或残基必须通过使用其他软件包模键出来（例如 ChimeraX 等）。还要注意，pdb2gmx 并不是万能的，它只能生成由力场定义的残基（位于.rtp文件中，通常是蛋白质、核酸和数量非常有限的辅酶，如 NAD(H) 和 ATP）而不是任意分子的拓扑文件（这也就意味着如果你需要使用自定义的非标准氨基酸或者其他小分子，都需要一些额外的操作）。

现在结晶水已经没有了，并且我们已经验证了所有必要的原子都已经存在，此时的PDB文件应该只包含蛋白质原子，然后就可以将其输入到第一个GROMACS module pdb2gmx 中。pdb2gmx 的目的是生成三个文件:

1. 分子的拓扑结构文件（.top）。
2. 位置限制文件（posre.itp）。
3. 处理后的结构文件（.gro）。

拓扑文件(默认生成名为：topol.top)包含了在模拟中定义分子所需的所有信息。这些信息包括非键参数(原子类型和电荷)以及键参数(键、角度和二面体)。具体内容将在后续介绍(也可以参考此处 note1)。

执行 pdb2gmx 命令:

gmx pdb2gmx -f 1AKI_clean.pdb -o 1AKI_processed.gro -water spce

这个结构文件（1AKI_clean.pdb）将由 pdb2gmx 处理，并且会提示您选择一个力场:

Select the Force Field:
From '/usr/local/gromacs/share/gromacs/top':
 1: AMBER03 protein, nucleic AMBER94 (Duan et al., J. Comp. Chem. 24, 1999-2012, 2003)
 2: AMBER94 force field (Cornell et al., JACS 117, 5179-5197, 1995)
 3: AMBER96 protein, nucleic AMBER94 (Kollman et al., Acc. Chem. Res. 29, 461-469, 1996)
 4: AMBER99 protein, nucleic AMBER94 (Wang et al., J. Comp. Chem. 21, 1049-1074, 2000)
 5: AMBER99SB protein, nucleic AMBER94 (Hornak et al., Proteins 65, 712-725, 2006)
 6: AMBER99SB-ILDN protein, nucleic AMBER94 (Lindorff-Larsen et al., Proteins 78, 1950-58, 2010)
 7: AMBERGS force field (Garcia & Sanbonmatsu, PNAS 99, 2782-2787, 2002)
 8: CHARMM27 all-atom force field (CHARM22 plus CMAP for proteins)
 9: GROMOS96 43a1 force field
10: GROMOS96 43a2 force field (improved alkane dihedrals)
11: GROMOS96 45a3 force field (Schuler JCC 2001 22 1205)
12: GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656)
13: GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656)
14: GROMOS96 54a7 force field (Eur. Biophys. J. (2011), 40,, 843-856, DOI: 10.1007/s00249-011-0700-9)
15: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)

力场包含将被写入拓扑文件的信息。这是一个非常重要的选择！你应该通读每一个力场，然后决定哪一个最适合你所需要的的情况（对于各个力场的使用场景，可以参考李继存科学网博客）。在本教程中，我们将使用 all-atom OPLS force field，所以在命令提示符下输入 15，回车。

pdb2gmx 的一些常用附加选项如下:

• -ignh：忽略pdb文件中的氢原子，对NMR结构尤其有效。否则，如果氢原子存在，它们的命名必须与GROMACS中的力场所要求的的完全一致。由于存在着不同的命名惯例，所以处理氢原子有时会非常恶心！如果您需要保留初始的氢原子坐标，但是需要重命名，可以使用 sed 命令（详细用法请见这里）。
• -ter：分配N-、C-端电荷。
• -inter：分配Glu, Asp, Lys, Arg和His的电荷，选择参与生成二硫键的Cys。

您现在已经生成了三个新文件:1AKI_processed.gro，topol.top，posre.itp。1 AKI_processed.gro 是一个gromacs格式的结构文件，包含力场中定义的所有原子(即H原子被添加到蛋白质的氨基酸中)。topol.top 是系统拓扑(稍后会详细介绍)。posre.itp 文件包含用于限制重原子位置的信息(稍后详细介绍)。

注意: .gro文件并非必须，GROMACS可以处理许多不同的文件格式，.gro只是编写坐标文件命令的默认值。它是一种非常紧凑的格式，但精度有限。如果您喜欢使用 .pdb 格式，那么您所需要做的就是指定一个适当的文件名，将. pdb 作为扩展名输出。pdb2gmx 的目的是生成一个符合力场的拓扑文件；输出结构在很大程度上是辅助作用，目的是为了方便用户。关于不同的格式，请参阅GROMACS手册。

The Topology

让我们看看输出拓扑(topol.top)中有什么内容。再次使用纯文本编辑器查阅其内容。在几个注释行之后(前面有;)，你会发现如下:

#include "oplsaa.ff/forcefield.itp"

这一行表明力场参数调用自OPLS-AA力场。它位于文件的开头，表示所有后续参数都来自这个力场。另一比较重要的行是 [moleculetype]，在它下面你可以看到

; Name       nrexcl
Protein_A    3

名称“Protein_A”定义了分子名称，这是根据 .pdb 文件中蛋白编号为 chain A 自动命名的。对键合近邻的排除数为3，关于排除的更多信息可从GROMACS手册上找到（也可以查看 note1 来了解大概），对此信息的讨论超出了本教程的范围。

下一节定义了蛋白质中的原子 [atom]，以列的形式显示:

[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  typeB    chargeB      massB
; residue   1 LYS rtp LYSH q +2.0
     1   opls_287      1   LYS       N      1       -0.3    14.0067   ; qtot -0.3
     2   opls_290      1   LYS      H1      1       0.33      1.008   ; qtot 0.03
     3   opls_290      1   LYS      H2      1       0.33      1.008   ; qtot 0.36
     4   opls_290      1   LYS      H3      1       0.33      1.008   ; qtot 0.69
     5  opls_293B      1   LYS      CA      1       0.25     12.011   ; qtot 0.94
     6   opls_140      1   LYS      HA      1       0.06      1.008   ; qtot 1

注释

• nr:原子编号
• type: 原子类型
• resnr: 氨基酸残基编号
• residue: 氨基酸残基名称
• atom: 原子名
• cgnr: 电荷组编号，有助于加速计算。gromacs中组的定义可以参见手册。
• charge:电荷量，qtot注释为分子总电荷量
• mass: 相对原子质量
• typeB, chargeB, massB: 用于自由能微扰（在此不讨论，可参考自由能计算教程）

后面的部分包括[ bond ]、[ pairs ]、[ angles ]和[ dihedrals ]。其中一些部分是不言自明的(键、角和二面角)。与这些部分相关的参数和函数类型在 GROMACS 手册的第5章中有详细说明。特殊的1-4相互作用包含在“pairs”中(GROMACS手册的5.3.4节)。

文件的其余部分包括定义其他一些有用/必要的拓扑，如“posre.itp "文件由 pdb2gmx 生成；它定义了一个力常数，用于在平衡过程中保持原子的位置(稍后详细介绍)。

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

上述语句标志着“Protein_A”分子类型定义的结束。拓扑文件的其余部分专门用于定义其他分子并提供系统级的描述。下一个分子类型(默认)是溶剂，在本实例中为 SPC/E water。其他典型的水模型的选择包括SPC, TIP3P 和 TIP4P。我们通过“-water spce”在进行pdb2gmx选择相应的水分子模型。不同的水模型的优点总结见此，请注意，并不是所有这些模型都存在于GROMACS中。

; Include water topology
#include "oplsaa.ff/spce.itp"

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
;  i funct       fcx        fcy        fcz
   1    1       1000       1000       1000
#endif

如你所见，水也可以被位置约束，使用1000 kJ mol-1 nm-2的力常数(kpr)。离子参数如下:

; Include generic topology for ions
#include "oplsaa.ff/ions.itp"

最后是系统级的定义：

[ system ]
; Name
LYSOZYME

[ molecules ]
; Compound        #mols
Protein_A           1

[ system ]给出了模拟时输出文件中的系统名，[ molecules ]列出了系统中所包含的所有分子。

对于[ molecules ]有几点值得关注的：

1. 所列出的分子顺序必须与坐标文件（例如.gro）中完全一致；
2. 分子名称必须与其[ moleculetype ]名称一致。

如果您不能满足这些要求，那么您将从 grompp 过程(稍后讨论)中得到名称不匹配、分子未被找到或其他报错。现在我们已经检查了拓扑文件的内容，我们可以继续构建我们的系统了。

Solvation

现在您已经熟悉了GROMACS拓扑的内容，现在可以继续构建我们的系统了。在这个例子中，我们将模拟一个简单的水体系。在不同溶剂中模拟蛋白质和其他分子也是可行的，只要对所涉及的所有分子种类都有良好的参数。

有两个步骤来定义盒子和填充溶剂:

1. 使用 editconf 模块定义盒子大小；
2. 使用 solvate 模块在盒子中填充水（溶剂）。

你现在面临着如何处理这个元胞（unit cell）的选择。出于本教程的目的，我们将使用一个简单的立方体作为单元格。随着你对周期性边界条件和盒子类型的熟悉，作者强烈推荐菱形十二面体，因为它的体积是相同周期距离的立方体盒子的71%，从而节省了溶解蛋白质所需添加的水分子的数量。

让我们用editconf来定义这个box:

gmx editconf -f 1AKI_processed.gro -o 1AKI_newbox.gro -c -d 1.0 -bt cubic

上述命令将蛋白质置于盒子中心(-c)，并使其距离盒子边缘至少1.0 nm(-d 1.0)。盒子类型定义为立方体(-bt cubic)。到盒子边缘的距离是一个重要的参数。由于我们将使用周期边界条件（详细请参考使用手册第三章），我们必须满足最小图像约定。也就是说，蛋白质永远不会看到它的周期图像，否则计算出来的力将是假的。指定溶质盒距离为1.0 nm意味着蛋白质的任意两个周期图像之间至少有2.0 nm。这个距离对于模拟中常用的任何截止方案都是足够的。

既然我们已经定义了一个盒子，我们可以把它装满溶剂(水)。溶剂化是用溶剂来完成的:

gmx solvate -cp 1AKI_newbox.gro -cs spc216.gro -o 1AKI_solv.gro -p topol.top

蛋白质的配置(-cp)由之前 editconf 步骤生成，溶剂的配置(-cs)是标准 GROMACS 安装的一部分。我们使用的是 spc216.gro，是一个通用的平衡的三点溶剂模型。您可以使用 spc216.gro 作为 SPC、SPC/E 或 TIP3P 的配置，因为它们都是三点水模型。输出被称为 1AKI_solv.gro。我们还需告诉 solvate 拓扑文件的名称(topol.top)，以便修改。注意 topol.top 中[ molecules ]的变化:

[ molecules ]
; Compound  #mols
Protein_A       1 
SOL         10832 

solvate 指令所做的就是记录它加入了多少水分子，然后把它写在你的拓扑结构上，以反映已经发生的变化。请注意，如果您使用任何其他(非水)溶剂，solvate 不会对您的拓扑结构做出这些更改！它与更新水分子的兼容性是强制编码的（意味着你需要手动修改 .top 文件）。

Adding Ions

我们现在有一个包含带电蛋白质的溶剂化系统。pdb2gmx 的输出结果告诉我们，该蛋白质的净电荷为 +8e (基于其氨基酸组成)。如果您在 pdb2gmx 输出中漏看了这个信息，请查看 topol.top 中的[ atoms ]的最后一行“qtot 8.”。由于系统不以净电荷存在，我们必须向我们的系统中加入离子（来平衡电荷）。

在GROMACS中添加离子的指令为 genion。genion 所做的是读取拓扑结构，用用户指定的离子替换水分子。该输入称为运行输入文件，扩展名为 .tpr；这个文件是由GROMACS的 grompp 模块(GROMACS pre -processor ，gromacs预处理)生成的，稍后我们运行第一个模拟时也将使用这个模块（有关 .tpr 文件的内容请参考 note1）。

我们先来介绍 grompp 指令。grompp 所做的是处理坐标文件和拓扑文件(用来描述分子的那个)以生成原子级输入(.tpr)。.tpr文件包含系统中所有原子的所有参数。

要用 grompp 生成一个.tpr文件，我们需要一个附加的输入文件，扩展名为 .mdp (分子动力学参数文件)；grompp 将把 .mdp 文件中指定的参数与坐标和拓扑信息组合在一起，生成一个 .tpr 文件。.mdp 文件通常用于运行能量最小化或MD模拟，但在本例中，它只是用于生成系统的原子描述。可以在这里下载本示例 .mdp 文件。

实际上，在此步骤中使用的 .mdp 文件可以包含任何合理的参数组合。我通常使用能量最小化脚本，因为它们非常基本，不涉及任何复杂的参数组合。请注意，本教程提供的文件仅用于 OPLS-AA 力场。其中的某些设置（特别是非键相互作用设置）对于其他力场来说是不同的。

.tpr生成指令:

gmx grompp -f ions.mdp -c 1AKI_solv.gro -p topol.top -o ions.tpr

现在我们在二进制文件 ions.tpr 中写入了系统的原子级别描述。然后我们可以将它用于 genion 指令：

gmx genion -s ions.tpr -o 1AKI_solv_ions.gro -p topol.top -pname NA -nname CL -neutral

选择13组“SOL”嵌入离子。你绝对不会想要用离子取代部分蛋白质。

在 genion 命令中，我们提供结构/状态文件(-s)作为输入；生成 .gro 文件作为输出(-o)；处理拓扑文件(-p)进行水分子的移除和离子的添加；定义正负离子名称(分别为 -pname 和 -nname)；告诉 genion 只添加必要的离子来中和蛋白质上的净电荷，通过添加正确数量的负离子(-neutral，在这种情况下将添加8个Cl-离子来抵消蛋白质上的+8电荷)。除了简单地通过指定 -neutral 和 -conc 选项来中和系统之外，您还可以使用 genion 添加指定浓度的离子。有关如何使用这些选项，请参阅genion部分的手册。

用 -pname 和 -nname 指定的离子名称在以前的GROMACS版本中是与力场相关的，但在4.5版本中进行了标准化。指定的离子名称总是元素符号的全部大写字母（这也是[ moleculetype ]的名称），然后写入拓扑文件。根据力场的不同，残留物或原子名可能会附加或不附加电荷的符号(+/-)。不要在 genion 命令中使用 atom 或 residual 名称，否则您将在后续步骤中遇到错误。

现在你的[ molecules ]应该长这样:

[ molecules ]
; Compound      #mols
Protein_A         1
SOL           10636
CL                8

Energy Minimization

溶剂化且电荷中性的系统现在已经组装好了。在我们开始动力学之前，我们必须确保系统没有空间冲突或不适当的几何形状。我们需要通过一种称为能量最小化(EM)的过程来放松结构。

能量最小化的过程很像离子的添加。我们将再次使用grompp将结构、拓扑和模拟参数组装到一个二进制输入文件(.tpr)中，但这一次，我们将不将.tpr传递给genion，而是通过GROMACS MD引擎 mdrun 运行能量最小化。

生成.tpr文件所需.mdp（在此）和指令：

gmx grompp -f minim.mdp -c 1AKI_solv_ions.gro -p topol.top -o em.tpr

在运行 genbox 和 genion 时前，确保你已经更新了你的 topol.top（指的是已经在文件中添加了相应的分子的数据），否则你会得到很多糟糕的错误消息(“坐标文件中的坐标数与拓扑文件不匹配”等)。

我们现在准备调用 mdrun 来执行EM:

gmx mdrun -v -deffnm em

-v是为不耐烦的人准备的：它使 mdrun 变得冗长，这样它在每一步都将进度输出到到屏幕上。-deffnm 将定义输入和输出的文件名。因此，如果没有将grompp输出命名为“em.tpr”，则必须使用 mdrun -s 来明确指定其名称（详细请参见手册 mdrun 部分）。在我们的例子中，我们将得到以下文件:

• em.log：EM进程的ASCII文本日志文件
• em.edr：二进制的能量文件
• em.trr：二进制全精度轨迹文件
• em.gro：能量最小化结构文件

有两个非常重要的因素来评估EM是否成功。第一个是势能(在EM过程结束时打印，即使没有-v)。E_pot 应该是负的，并且(对于水中的简单蛋白质)在10⁵-10⁶的数量级，这取决于系统的大小和水分子的数量。第二个重要的特征是最大力，F_max；该参数在minim.mdp中相应的设置为 “emtol = 1000.0”，表示目标 F_max 不大于1000 kJ mol^-1*nm^-1。有可能到达一个合理的 E_pot 与 F_max > emtol。如果发生这种情况，您的系统可能不够稳定来进行模拟。评估它可能发生的原因，也许改变你的最小化参数(积分算法，步数等)来优化它。

接下来让我们做一些分析。em.edr 文件包含了 GROMACS 在EM期间收集的所有能量项。你可以使用 GROMACS 的 energy 指令分析任何 .edr 文件:

gmx energy -f em.edr -o potential.xvg

在提示符处，键入“10 0”以选择Potential (10); zero (0)终止输入。你会看到E_pot的平均值，以及一个叫做“potential.xvg”。要绘制这些数据，需要使用 Xmgrace 绘图工具（或者Qtgrace）。一个Epot良好的、稳定的收敛的图像应该是这样的:

现在我们的能量最小化已经完成啦，接下来就开始真正的模拟了。

Equilibration

EM 确保了我们在几何形状和溶剂取向方面有一个合理的起始结构。要开始真正的动力学模拟，我们必须平衡蛋白质周围的溶剂和离子。如果我们在这一点上尝试无约束动力学，系统可能会崩溃。原因是在没有约束的情况下，溶剂分子在溶剂内必然是最优的，而对于溶质来说则不是。它需要达到我们希望模拟的温度并使得溶质(蛋白质)处于正确的方向。在我们达到正确的温度(基于动能的)后，我们将对系统施加压力，直到它达到适当的密度。

还记得我们用 pdb2gmx 生成的 posre.itp 文件吗?我们现在要用它对蛋白质的重原子(任何不是氢的原子)施加一种位置约束力（position restraining force），使运动只有在克服大量的能垒之后才可以发生。位置限制的效用是，它们允许我们在蛋白质周围平衡溶剂，而不需要增加蛋白质结构变化的变量。位置约束的原点(约束势为零的坐标)是由 grompp 的-r选项的坐标文件提供的。

平衡常分两阶段进行。第一阶段是在NVT系综下进行的(恒定的粒子数，体积和温度)。这个系综也称为“等温等容系综”或“正则系综”。这样一个过程的时间框架依赖于系统的内容，但在NVT中，系统的温度应该在期望的值上达到一个稳定的水平。如果温度还没有稳定下来，则需要额外的时间（来平衡温度）。通常，50-100 ps就足够了，本教程将进行100 ps的NVT平衡。根据您的计算机，这可能需要一段时间(如果16核左右并行运行，大约需要一个小时)。在这里获取 .mdp 文件。

我们将调用grompp和mdrun，就像我们在EM步骤中所做的那样:

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr

gmx mdrun -deffnm nvt

除了提供的附加指令之外，在GROMACS手册中还可以找到对所使用参数的完整解释。注意 .mdp 文件中的几个参数:

• gen_vel = yes：初速度的生成。使用不同的随机种子(gen_seed)给出不同的初始速度，从而在相同的起始结构进行多次(不同)模拟。
• tcoupl = V-rescale：速度重缩放温度耦合（velocity rescaling thermostat）是对Berendsen弱耦合方法的改进，该方法相比Berendsen温度耦合于能重现正确的动力学整体（详细请见手册3.4.8）。
• pcoupl = no：不使用压力耦合。

通过energy分析恒温过程：

gmx energy -f nvt.edr -o temperature.xvg

在提示符下输入“16 0”，选择系统温度并退出。结果图应该如下所示:

从图中可以清楚地看出，系统的温度很快达到目标值(300 K)，并且在剩余的平衡过程中保持稳定。对于这个系统，一个较短的平衡期(大约50 ps)可能已经足够了。

Equilibration, Phase 2

前一步的NVT平衡稳定了系统的温度。在收集数据之前，我们还必须稳定系统的压力(以及密度)。压力平衡是在NPT系综下进行的，其中粒子数、压力和温度都是恒定的。该系综也称为“等温等压系综”，与实验条件最为相似。

用于 100ps NPT 平衡的.mdp文件可以在这里找到。它与用于NVT平衡的参数文件没有很大的不同。注意增加了压力耦合部分，使用的是Parrinello-Rahman压力耦合方法（详细请参见手册3.4.9）。

其他的一些变化:

• continuation = yes：继续从NVT平衡阶段进行模拟；
• gen_vel = no：速度是从轨迹中读取的(见下)。

我们将调用 grompp 和 mrun，就像NVT平衡中进行的一样。注意，我们现在使用了 -t 以包括来自NVT平衡的检查点文件（这个文件包含继续模拟所需的所有状态变量）。为了保持NVT期间产生的速度，我们必须使用这个文件。坐标文件(-c)是NVT仿真的最终输出文件。

gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr

gmx mdrun -deffnm npt

分析恒压过程：

gmx energy -f npt.edr -o pressure.xvg

在提示符下输入“18 0”，选择系统压力并退出。结果图应该如下所示:

压力值在100 ps平衡阶段的过程中波动很大，但这种行为是意料之中的。这些数据的运行平均值在图中用红线表示。在平衡过程中，压力平均值为7.5±160.5 bar（注意参考压力设置为1 bar，那么这个结果可以接受吗？）。从巨大的均方根波动(160.5 bar)中可以清楚地看出，压力在MD模拟过程中是一个波动很大的量，因此从统计学上讲，我们无法区分得到的平均值(7.5±160.5 bar)和目标/参考值(1 bar)之间的差异。

分析密度变化，随后输入"24 0"：

gmx energy -f npt.edr -o density.xvg

和压强一样，密度的运行平均值也用红色表示。100 ps 过程的平均值为 1019±3 kg·m^-3，接近 1000 kg·m^-3 的实验值和 1008 kg·m^-3 的 SPC/E 模型的期望密度。SPC/E 水模型的参数与水的实验值非常接近。随着时间的推移，密度值非常稳定，这表明系统现在相对于压力和密度来说是很平衡的。

请注意:作者经常被问到为什么得到的密度值与教程显示的结果不匹配。这与压力相关项收敛的缓慢有关，因此您可能需要运行NPT平衡比此处指定的时间稍长一些。

Production MD

在两个平衡阶段完成后，系统现在就在期望的温度和压力下很好地平衡了。我们现在已经准备好解除位置限制，并运行MD进行数据收集。这个过程就像我们之前看到的一样，我们将利用检查点文件(在本例中，它现在包含保存压力耦合信息)来进行 grompp。我们将运行一个 1-ns MD 模拟，脚本可以在这里找到。

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr

grompp 将显示 PME 负载的估计值，这将指定有多少处理器应该用于PME计算，有多少处理器应该用于 PP 计算（详情请参阅 GROMACS 出版物和手册）：

Estimate for the relative computational load of the PME mesh part: 0.22

对于一个立方体盒子，最优设置将有一个 0.25 的 PME 负载(3:1=PP:PME，我们非常接近最优!)；对于十二面体箱体，最佳 PME 负载为0.33 (2:1=PP:PME)。当运行mdrun时，程序应该自动确定并分配 PP 和 PME 计算所需最佳处理器数量。因此，请确保为计算指定了适当的线程/核数（-nt X的值），以便获得最佳性能。

现在,执行mdrun:

gmx mdrun -deffnm md_0_1

Gromacs 模拟中断或正常结束后的续算方法

如果模拟中途停止（并不是崩溃，没有任何报错的那种停止（比如你的terminal突然shutdown了或者断电了）），可以使用前一部分模拟产生的日志文件续跑（其中-cpi表示使用检查点文件，-append默认为 yes，表示将模拟结果续写到原文件中）：

gmx mdrun -s md_0_1.tpr -cpi md_0_1.cpt -append

如果模拟正常结束后希望延长模拟时间，并输出到新的一系列名为md_0_2的文件中，可以：

gmx convert-tpr -s md_0_1.tpr -extend 1000 -o md_0_2.tpr

gmx mdrun -v -deffnm md_0_2 -cpi md_0_1.cpt -noappend

gmx convert-tpr为更改 .tpr 文件的指令；其中-extend后数值单位为 ps，表示将模拟时间延续多长；使用-noappend将打开新的输出文件，并将模拟序号名添加到所有输出文件名中。使用检查点文件（-cpi）时，输出将默认续写到以前写入的输出文件，除非使用-noappend或以前的输出文件都不存在（检查点文件除外）。

如果想要将输出接续到先前的输出文件中：

gmx convert-tpr -s md_0_1.tpr -extend 1000 -o md_0_1.tpr

gmx mdrun -v -deffnm md_0_1 -cpi md_0_1.cpt

只需修改并覆盖掉原先的 .tpr 文件并开始模拟（此处未使用-append是因为使用-cpi默认续写到先前的文件中）。

另外如果想让程序在后台运行可以尝试：

#Ctrl+z
bg %1
disown -h %1

Running GROMACS on GPU

在GROMACS 2018中，PME计算可以转移到图形处理单元(GPU)，这大大加快了模拟速度。使用Titan Xp GPU，这个系统可以以惊人的295秒/天的速度进行模拟。

从4.6版本开始，GROMACS支持使用GPU加速器来运行MD模拟。随着2018版的发布，非键交互和PME在GPU上计算，只在CPU核上计算键力。当搭建GROMACS时(安装说明请参见www.gromacs.org)，如果GPU硬件存在，它将被自动检测到。使用GPU加速的最低要求是有CUDA库和SDK，以及一个计算能力为 >= 2.0 的GPU。可以在这里找到一些比较常见的GPU及其规范的列表。

假设你有一个可用的GPU，mdrun命令如下:

gmx mdrun -deffnm md_0_1 -nb gpu

如果你有多个可用的 GPU，或者需要自定义如何通过混合并行化方案在 GROMACS 中分配工作，请参阅 GROMACS 手册和网页。这些技术细节超出了本教程的范围。

译者注：有关于 GPU 的调试可以参考http://bbs.keinsci.com/thread-13861-1-1.html

Analysis

既然我们已经模拟了蛋白质，现在应该对系统进行一些分析。什么类型的数据是重要的？这是在运行模拟之前要问的一个重要问题，因此您应该对希望在系统中收集到的数据类型有一些想法。在本教程中，将介绍一些基本的工具。

第一个是 trjconv，它被用作一种后处理工具，用于去除坐标、纠正周期性或手动改变轨迹(时间单位、帧频率等)。在这个练习中，我们将使用 trjconv 来修正系统中的任何周期性问题。例如蛋白质会扩散到整个元胞，并可能出现“断裂”或“跳跃”到盒子的另一边。为了调整这些行为，键入:

gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -pbc mol -center

选择1(“Protein”)作为居中分组，0(“System”)作为输出分组。我们将对这一“修正”的轨迹进行所有的分析。让我们先看看结构稳定性。GROMACS 有一个内置的用于 RMSD 计算的实用程序，称为 rms。要使用 rms，发出以下命令:

gmx rms -s md_0_1.tpr -f md_0_1_noPBC.xtc -o rmsd.xvg -tu ns

对于最小二乘拟合和RMSD计算组，选择4(“Backbone”)。-tu 将输出以ns表示的结果，即将轨迹用ps表示。这样做是为了使输出更加清晰（特别是如果你有一个很长的模拟- 1e+05 ps看起来不如100 ns好）。输出图将显示RMSD相对于最小化的平衡系统中的结构:

如果我们希望计算相对于晶体结构的 RMSD，使用以下指令:

gmx rms -s em.tpr -f md_0_1_noPBC.xtc -o rmsd_xtal.xvg -tu ns

把结果放在一起，看起来像这样:

两个时间序列的RMSD水平均降至~0.1 nm (1 Å)，表明该结构非常稳定。图之间的细微差别表明，在t = 0 ns 时的结构与晶体结构略有不同。这是意料之中的，因为它已经被能量最小化了，而且因为位置限制不是100%完美的，正如前面讨论过的。

More Analysis

蛋白质的旋转半径（the radius of gyration，Rg）是其紧密度的一种度量。如果一个蛋白质折叠很稳定，它理应保持一个相对稳定的Rg值。如果蛋白质失去了三级结构，它的Rg会随时间变化。让我们在模拟中分析溶菌酶的旋转半径:

gmx gyrate -s md_0_1.tpr -f md_0_1_noPBC.xtc -o gyrate.xvg

选择 group 1 (Protein)：

我们可以从合理且不变的Rg值中看出，在300k的1ns过程中，蛋白质在其紧凑(折叠)形式中保持得非常稳定。这个结果并不出乎意料，但说明了内置的GROMACS分析的高级能力。

Summary

现在，您已经使用GROMACS进行了分子动力学模拟，并分析了一些结果。本教程不应被视为万能的。我们可以用GROMACS进行更多类型的模拟(自由能计算、非平衡MD和正态模态分析，仅举几个例子)。为了提高效率和准确性，您还应该查阅文献和GROMACS手册，以便对这里提供的 .mdp 文件进行调整（有关 .mdp 选项的内容可以参考这里）。

略

Happy simulating!

tips

如果预估模拟时间太长而且忘记提前放入后台，可以用 bg 和 disown 的方法置入后台运行：

首先Ctrl+z暂停，

bg %1

disown -h %1