【3.2.1】gromacs教程--水中的溶菌酶
使用
source /opt/fastone/softwares/hpc-kits/hpc-kits.sh
source /data/software/gromacs/gromacs-2020.3/bin/GMXRC
查看帮助：
gmx help commands
或者：
gmx help selections
一、水中的溶菌酶 Lysozyme in Water
cd /data/user/sam/project/gromacs/test
从PDB上下载1aki.pdb
1.1 准备拓扑结构
注意，这种方法并非普遍适用（例如，紧密结合或以其他方式起作用的活性位点水分子的情况）。 对于我们这里的意图，我们不需要结晶水。
grep -v HOH 1aki.pdb > 1AKI_clean.pdb
始终检查.pdb文件中是否有注释MISSING下列出的条目，因为这些条目表示晶体结构中不存在的原子或完整残基。 终端区域可能不存在，并且可能不会出现动态问题。 内部序列不完整或原子缺失的任何氨基酸残基都会导致pdb2gmx失败。 这些缺失的原子/残基必须使用其他软件包进行建模。 另请注意，pdb2gmx并不是魔术。 它不能生成任意分子的拓扑，只能生成由力场定义的残基（在* .rtp文件中-通常是蛋白质，核酸和非常有限量的辅助因子，如NAD（H）和ATP）。
现在，晶体水不见了，并且我们已经验证了所有必需的原子都存在，PDB文件应该仅包含蛋白质原子，并且可以输入到第一个GROMACS模块pdb2gmx中。 pdb2gmx的目的是生成三个文件：
- 分子的拓扑。
- 位置限制文件。
- 后处理的结构文件。
拓扑（默认为topol.top）包含在模拟中定义分子所需的所有信息。 此信息包括非键合参数（原子类型和电荷）以及键合参数（键合，角度和二面体）。 拓扑生成后，我们将对其进行更详细的研究。
gmx pdb2gmx -f 1AKI_clean.pdb -o 1AKI_processed.gro -water spce
接着选择力场
Select the Force Field:
From '/data/software/gromacs/gromacs-2020.3/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)
力字段将包含将写入拓扑的信息。 这是一个非常重要的选择！ 您应该始终仔细阅读每个力场，并确定最适合您的情况。 在本教程中，我们将使用全原子OPLS力场，因此在命令提示符下键入15，然后单击“ Enter”。
还有许多其他选项可以传递给pdb2gmx。 这里列出了一些常用的：
-ignh：忽略PDB文件中的H原子； 特别适用于NMR结构。 否则，如果存在H原子，则它们必须精确地位于GROMACS中力场所期望的名称中。 存在不同的约定，因此处理H原子有时会令人头疼！ 如果您需要保留初始H坐标，但需要重命名，那么Linux sed命令就是您的朋友。
-ter：为N和C终端交互式分配电荷状态。
-inter：交互分配Glu，Asp，Lys，Arg和His的电荷状态； 选择哪些Cys参与二硫键。
现在，您已经生成了三个新文件：1AKI_processed.gro，topol.top和posre.itp。
- 1AKI_processed.gro是GROMACS格式的结构文件，其中包含在力场内定义的所有原子（即，H原子已添加到蛋白质中的氨基酸上）。
- topol.top文件是系统拓扑（稍后将对此进行详细介绍）。
- posre.itp文件包含用于限制重原子位置的信息（稍后会详细介绍）。
最后一点：许多用户认为.gro文件是强制性的。这不是真的。 GROMACS可以处理许多不同的文件格式，.gro只是编写坐标文件的命令的默认格式。这是一种非常紧凑的格式，但是精度有限。例如，如果您喜欢使用PDB格式，则只需指定一个扩展名为.pdb的文件名作为输出即可。 pdb2gmx的目的是产生一个符合力场的拓扑。输出结构在很大程度上是此目的的副作用，目的是为了方便用户。该格式几乎可以是任何您喜欢的格式（有关不同格式，请参阅GROMACS手册）。
1.2 检查拓扑结构
让我们看一下输出拓扑中的内容（topol.top）。 同样，使用纯文本编辑器检查其内容。 在几个注释行（以;开头）之后，您将找到以下内容：
#include "oplsaa.ff/forcefield.itp"
该行调用OPLS-AA力字段内的参数。 它在文件的开头，指示所有后续参数均从该力场派生。 下一个重要的行是[Moleculartype]，您将在下面找到
; Name nrexcl
Protein_A 3
根据蛋白质在PDB文件中被标记为链A的事实，名称“ Protein_A”定义了分子名称。 bonded neighbors有3个例外。 有关排除的更多信息，请参见GROMACS手册。 关于此信息的讨论超出了本教程的范围。
下一节定义了蛋白质中的[原子]。 该信息以列的形式显示：
[ 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：氨基酸残基名称。请注意，该残基在PDB文件中为“ LYS”。 .rtp条目“ LYSH”的使用表示残基已质子化（中性pH处于主要状态）。
atom：原子名称
cgnr：Charge group number。 电荷组定义整数电荷的单位； 它们有助于加快计算速度
charge：Self-explanatory 。 “ qtot”描述符是分子上电荷的运行总计
mass：lso self-explanatory
typeB，chargeB，massB：用于自由能微扰（此处不讨论）
随后的部分包括[键]，[对]，[角度]和[二面体，dihedrals]。 这些部分中的一些是不言自明的（键，角和二面角）。 与这些部分相关的参数和功能类型在GROMACS手册的第5章中进行了详细说明。 “成对”下包含特殊的1-4交互（GROMACS手册的5.3.4节）。
该文件的其余部分涉及从位置限制开始定义其他一些有用/必要的拓扑。 “ posre.itp”文件是由pdb2gmx生成的； 它定义了一个力常数，用于在平衡过程中将原子保持在适当的位置（稍后会详细介绍）。
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
这样就结束了“ Protein_A”分子类型的定义。 拓扑文件的其余部分专用于定义其他分子并提供系统级描述。 下一个分子类型（默认情况下）是溶剂，在这种情况下为SPC / E水。 水的其他典型选择包括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]伪指令给出了在仿真过程中将被写入输出文件的系统的名称。 [分子]指令列出了系统中的所有分子。
[ system ]
; Name
LYSOZYME
[ molecules ]
; Compound #mols
Protein_A 1
有关[分子]指令的一些关键说明：
- 列出的分子的顺序必须与坐标（在本例中为.gro）文件中的分子顺序完全匹配。
- 列出的名称必须与每个物种的[分子类型]名称匹配，而不是残基名称或其他任何名称。
如果您无法随时满足这些具体要求，则会从grompp（稍后讨论）中遇到致命错误，涉及名称不匹配，找不到分子或其他一些错误。
现在，我们已经检查了拓扑文件的内容，我们可以继续构建系统。
1.3 Defining the Unit Cell & Adding Solvent 定义晶胞并添加溶剂
现在您已经熟悉了GROMACS拓扑的内容，是时候继续构建我们的系统了。 在这个例子中，我们将模拟一个简单的水系统。 前提是可以为所有涉及的物种提供良好的参数，但是可以在不同的溶剂中模拟蛋白质和其他分子。
定义盒子并将其填充溶剂有两个步骤：
- 使用editconf模块定义盒子尺寸。
- 使用溶剂化物模块（以前称为genbox）向盒子装满水。
现在，您可以选择如何处理晶胞。 就本教程而言，我们将使用一个简单的立方体框作为单位单元。 随着您对周期性边界条件和盒类型的适应性提高，我强烈建议使用菱形十二面体，因为其体积约为相同周期距离的立方盒的71％，从而节省了需要添加的水分子数量 溶解蛋白质。
让我们使用editconf定义盒子：
gmx editconf -f 1AKI_processed.gro -o 1AKI_newbox.gro -c -d 1.0 -bt cubic
上面的命令将蛋白质放在方框（-c）的中央，并将其放置在距方框边缘至少1.0 nm（-d 1.0）的位置。 框类型定义为多维数据集（-bt立方）。 到盒子边缘的距离是一个重要的参数。 由于我们将使用周期性边界条件，因此我们必须满足最小图像约定。 也就是说，蛋白质永远都不会看到其周期性图像（ periodic image），否则计算出的力将是虚假的。 指定溶质箱距离为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，我们告诉溶剂化拓扑文件的名称（topol.top），以便可以对其进行修改。 请注意topol.top的[分子]指令的更改：
[ molecules ]
; Compound #mols
Protein_A 1
SOL 10832
溶剂化物的作用是跟踪已添加的水分子数量，然后将其写入拓扑以反映所做的更改。 请注意，如果使用任何其他（非水）溶剂，则溶剂化物不会对拓扑结构做出这些更改！ 它与更新水分子的兼容性是硬编码的。
1.4 Adding Ions 添加离子
现在，我们有了包含带电蛋白质的溶剂化系统。 pdb2gmx的输出告诉我们，该蛋白质的净电荷为+ 8e（基于其氨基酸组成）。如果您在pdb2gmx输出中错过了此信息，请查看topol.top中[atoms]指令的最后一行；它应该（部分）显示为“ qtot 8”。由于生命不存在净电荷，因此我们必须向系统中添加离子。
在GROMACS中添加离子的工具称为genion。 Genion所做的工作是通过拓扑读取的，并用用户指定的离子替换水分子。该输入称为运行输入文件，其扩展名为.tpr;。该文件是由GROMACS grompp模块（GROMACS预处理器）生成的，稍后在我们运行第一个模拟时也将使用该文件。 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通过添加正确数量的负离子（中性，仅添加中和蛋白质上净电荷所必需的离子）。添加8个Cl-离子以抵消蛋白质上的+8电荷）。除了简单地通过结合指定-neutral和-conc选项来中和系统之外，还可以使用genion添加指定的离子浓度。有关如何使用这些选项的信息，请参考genion手册页。
用-pname和-nname指定的离子名称在以前的GROMACS版本中是特定于力场的，但在4.5版中已标准化。指定的离子名称始终是所有大写字母中的元素符号，即随后写入拓扑中的[分子类型]名称。残留物或原子名称可以或可以不附加电荷的符号（+/-），具体取决于力场。不要在genion命令中使用原子或残基名称，否则在后续步骤中会遇到错误。
您的[分子]指令现在应如下所示：
[ molecules ]
; Compound #mols
Protein_A 1
SOL 10636
CL 8
（这个图用什么看的？ MOE?）
1.5 能量最小化 Energy Minimization
现在组装溶剂化的电子中性系统。 在开始动力学之前，必须确保系统没有空间冲突或不适当的几何形状。 通过称为能量最小化（EM）的过程来放松结构。
EM的过程非常类似于离子的添加。 我们将再次使用grompp将结构，拓扑和模拟参数组装到二进制输入文件（.tpr）中，但是这一次，我们将通过GROMACS MD引擎，mdrun 运行能量最小化，而不是将.tpr传递给genion。 。
使用以下输入参数文件，使用grompp组装二进制输入：
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标志显式指定其名称。 在我们的例子中，我们将获得以下文件：
em.log：EM进程的ASCII文本日志文件
em.edr：二进制能量文件
em.trr：二进制全精度轨迹
em.gro：节能结构
有两个非常重要的因素需要评估以确定EM是否成功。
- 第一个是势能（在EM过程结束时打印，即使没有-v也是如此）。 Epot应为负值，（对于水中的简单蛋白质）应在10^5-10^6左右，具体取决于系统大小和水分子数量。
- 第二个重要特征是最大力Fmax，为其设置的目标值设置为minim.mdp-“ emtol = 1000.0”-表示目标Fmax不大于1000 kJ mol-1 nm-1。 Fmax> emtol可以得出一个合理的Epot。 如果发生这种情况，您的系统可能不够稳定，无法进行仿真。 评估可能发生这种情况的原因，并可能更改最小化参数（积分器，emstep等）。
让我们做一点分析。 em.edr文件包含GROMACS在EM期间收集的所有能量项。 您可以使用GROMACS能量模块分析任何.edr文件：
gmx energy -f em.edr -o potential.xvg
在提示符下，键入“ 10 0”以选择“电位”（10）； 零（0）终止输入。 将显示Epot的平均值，并且将写入一个名为“ potential.xvg”的文件。 要绘制此数据，您将需要Xmgrace绘制工具。 生成的图看起来应该像这样，说明Epot很好，稳定地融合：
1.6 Equilibration 平衡
EM确保我们在几何形状和溶剂取向方面具有合理的起始结构。要开始真正的动力学，我们必须平衡蛋白质周围的溶剂和离子。如果我们此时尝试不受限制的动力学，系统可能会崩溃。原因是溶剂大部分在自身内部进行了优化，而不一定需要使用溶质。需要将其加热到我们希望模拟的温度，并建立关于溶质（蛋白质）的正确方向。在达到正确的温度（基于动能）之后，我们将向系统施加压力，直到达到适当的密度。
还记得pdb2gmx很久以前生成的posre.itp文件吗？我们现在要使用它。 posre.itp的目的是对蛋白质的重原子（任何不是氢的原子）施加位置限制力。允许运动，但前提是要克服大量的能量消耗。位置限制的用途是它们使我们能够平衡蛋白质周围的溶剂，而不会增加蛋白质结构变化的变量。通过传递给grompp的-r选项的坐标文件提供位置约束的原点（约束势为零的坐标）。
平衡通常分两个阶段进行。第一阶段是在NVT集成下进行的（恒定的粒子数，体积和温度）。该合奏也称为“等温-等速”或“规范”。此过程的时间范围取决于系统的内容，但在NVT中，系统的温度应达到所需值的平稳状态。如果温度尚未稳定，则需要更多时间。通常，50-100 ps就足够了，我们将为此练习进行100 ps的NVT平衡。根据您的计算机，这可能需要一段时间（如果在16个内核左右并行运行，则不到一个小时）。在此处获取.mdp文件。
就像在EM步骤中一样，我们将调用grompp和mdrun：
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：速度重定标恒温器是对Berendsen弱耦合方法的改进，该方法没有重现正确的动力学集合。
pcoupl =否：不应用压力耦合。
让我们再次使用能量来分析温度变化：
gmx energy -f nvt.edr -o temperature.xvg
在提示符下键入“ 16 0”以选择系统温度并退出。 结果图应如下所示：
从图中可以明显看出，系统温度迅速达到目标值（300 K），并且在其余的平衡过程中保持稳定。 对于该系统，较短的平衡周期（大约50 ps）可能已足够。
上一步NVT平衡可稳定系统温度。在收集数据之前，我们还必须稳定系统的压力（以及密度）。在NPT集成下进行压力平衡，其中颗粒数，压力和温度均恒定。该合奏也称为“等温-等压”合奏，与实验条件最相似。
可以在此处找到用于100 ps NPT平衡的.mdp文件。它与用于NVT平衡的参数文件没有太大不同。注意使用Parrinello-Rahman恒压器增加了压力耦合部分。
其他一些变化：
continuation = yes：我们正在从NVT平衡阶段继续进行仿真
gen_vel = no：从轨迹读取速度（见下文）
就像调用NVT平衡一样，我们将调用grompp和mdrun。请注意，我们现在包括-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，这样的结果可以接受吗？ 压力是在MD模拟过程中波动很大的量，从较大的均方根波动（160.5 bar）可以清楚地看出，因此从统计学上讲，人们无法区分所得平均值之间的差异（7.5±160.5 bar） ）和目标/参考值（1巴）。
让我们看一下密度，这次使用能量并在提示符下输入“ 24 0”。
gmx energy -f npt.edr -o density.xvg
与压力一样，密度的移动平均值也以红色绘制。 在100 ps的过程中，平均值为1019±3 kg m-3，接近实验值1000 kg m-3和SPC / E模型的预期密度1008 kg m-3。 SPC / E水模型的参数紧密复制了水的实验值。 密度值随时间变化非常稳定，这表明该系统现在在压力和密度方面达到了很好的平衡。
请注意：我经常会问到为什么所获得的密度值与我的结果不匹配的问题。 与压力有关的项收敛较慢，因此您可能需要运行NPT平衡稍长于此处指定的时间。
1.8 Production MD 生产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 4出版物和手册。
Estimate for the relative computational load of the PME mesh part: 0.22
对于立方盒，最佳设置的PME负载为0.25（3：1 PP：PME-我们非常接近最佳值！）； 对于十二面体盒，最佳PME负载为0.33（2：1 PP：PME）。 执行mdrun时，程序应自动确定为PP和PME计算分配的最佳处理器数量。 因此，请确保为计算指定适当数量的线程/核（-nt X的值），以便获得最佳性能。
现在，执行mdrun：
gmx mdrun -deffnm md_0_1
在GROMACS 2018中，可以将PME计算卸载到图形处理单元（GPU）中，从而大大加快了仿真速度。 使用Titan Xp GPU，该系统可以惊人的295 ns /天的速度进行仿真！
在GPU上运行GROMACS
从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手册和网页。这些技术细节超出了本教程的范围。
1.9 数据分析
现在我们已经模拟了蛋白质，我们应该在系统上进行一些分析。 哪些类型的数据很重要？ 在运行模拟之前，这是一个非常重要的问题，因此您应该对要在自己的系统中收集的数据类型有一些想法。 对于本教程，将介绍一些基本工具。
第一个是trjconv，它用作后处理工具以去除坐标，校正周期性或手动更改轨迹（时间单位，帧频等）。 在本练习中，我们将使用trjconv考虑系统中的任何周期性。 蛋白质将扩散穿过单位细胞，并可能出现“折断”或“跳跃”到盒子的另一侧。 要解决此类行为，请发出以下命令：
gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -pbc mol -center
选择1（“蛋白质”）作为要居中的组，选择0（“系统”）进行输出。 我们将对这种“校正”的轨迹进行所有分析。 首先让我们看一下结构稳定性。 GROMACS具有一个用于RMSD计算的内置实用程序，称为rms。 要使用rms，请发出以下命令：
gmx rms -s md_0_1.tpr -f md_0_1_noPBC.xtc -o rmsd.xvg -tu ns
对于最小二乘拟合和RMSD计算组，请选择4（“主干”）。 即使轨迹是用ps编写的，-tu标志也将以ns的形式输出结果。 这样做是为了使输出清晰（特别是在您进行长时间仿真时-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 levels低至〜0.1 nm（1Å），表明该结构非常稳定。 曲线之间的细微差异表明，t = 0 ns时的结构与该晶体结构略有不同。 如前所述，这是可以预期的，因为它已实现了能量最小化，并且因为位置约束不是100％完美的。
蛋白质的回转半径是其紧密度的量度。 如果蛋白质被稳定折叠，它将可能保持相对稳定的Rg值。 如果蛋白质解折叠，其Rg将随着时间变化。 让我们在模拟中分析溶菌酶的回转半径：
gmx gyrate -s md_0_1.tpr -f md_0_1_noPBC.xtc -o gyrate.xvg
选择组1（蛋白质）进行分析。
从合理不变的Rg值可以看出，该蛋白质在300 K下在1 ns的过程中以紧密（折叠）形式保持非常稳定的状态。这个结果并不出乎意料，但说明了已建立的GROMACS分析的先进能力 。
1.10 总结
您现在已经使用GROMACS进行了分子动力学模拟，并分析了一些结果。 本教程不应视为全面的教程。 可以使用GROMACS进行更多类型的仿真（自由能计算，非平衡MD和正常模式分析，仅举几例）。 您还应该阅读文献和GROMACS手册，以调整此处提供的.mdp文件，以提高效率和准确性。
如果您有改进本教程的建议，发现错误或其他不清楚的地方，请随时给我发送电子邮件。 请注意：这不是邀请我发送有关GROMACS问题的电子邮件。 我不会宣传自己为私人家教或个人帮助服务。 这就是 https://mailman-1.sys.kth.se/mailman/listinfo/gmx-users 。 我可以在这里为您提供帮助，但仅限于为整个社区提供服务，而不仅仅是最终用户。
参考资料
