【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的目的是生成三个文件:

  1. 分子的拓扑。
  2. 位置限制文件。
  3. 后处理的结构文件。

拓扑(默认为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

有关[分子]指令的一些关键说明:

  1. 列出的分子的顺序必须与坐标(在本例中为.gro)文件中的分子顺序完全匹配。
  2. 列出的名称必须与每个物种的[分子类型]名称匹配,而不是残基名称或其他任何名称。

如果您无法随时满足这些具体要求,则会从grompp(稍后讨论)中遇到致命错误,涉及名称不匹配,找不到分子或其他一些错误。

现在,我们已经检查了拓扑文件的内容,我们可以继续构建系统。

1.3 Defining the Unit Cell & Adding Solvent 定义晶胞并添加溶剂

现在您已经熟悉了GROMACS拓扑的内容,是时候继续构建我们的系统了。 在这个例子中,我们将模拟一个简单的水系统。 前提是可以为所有涉及的物种提供良好的参数,但是可以在不同的溶剂中模拟蛋白质和其他分子。

定义盒子并将其填充溶剂有两个步骤:

  1. 使用editconf模块定义盒子尺寸。
  2. 使用溶剂化物模块(以前称为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 。 我可以在这里为您提供帮助,但仅限于为整个社区提供服务,而不仅仅是最终用户。

参考资料

个人公众号,比较懒,很少更新,可以在上面提问题,如果回复不及时,可发邮件给我: tiehan@sina.cn

Sam avatar
About Sam
专注生物信息 专注转化医学