Gromacs教程-水中的溶菌酶

GROMACS 教程 水中的溶菌酶 Justin Lemkul Department of Biochemistry, Virginia Tech YongMa2008@小木虫 译 yo ng ma 20 08 @小 木 虫 写在前面: 1 本人没有系统学习过 MD, 本教程是本人开始自学 GROMACS 的 入门教程，在课余时间的，对于很多专有词汇，不是很懂，翻译 的不尽正确，欢迎批评指正 2 由于版本问题，本教程中某些命令行或需要输入的变量可能会有不 同，请自行斟酌 3 本教程中红色为命令行，蓝色为超级链接，橙黄色为程序运行时终 端显示文字或文件内部文字 4 由于 Word 排版问题，某些命令行中的空格不是很明显，请注意；， 同时由于时间仓促难免出现排版错误，请见谅 5 原版教程链接如下，强烈建议有兴趣的童鞋学习原版 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx‐tutorials/lysozyme/index.html yo ng ma 20 08 @小 木 虫 本例将指导新用户进行设置一个蛋白质(lysozyme)加上离子在水盒子里的 模拟过程。每个步骤将包含对输入输出文件一般应用的典型设置的解释。 这个教程假定你正在用 gromacs 的 4.5.x 版本。 第一步，准备拓扑 我们必须先下载我们要用的蛋白质结构。在这个教程中，我们将用鸡蛋清溶 菌酶(PDB代码1AKI). 去 RCSB(http://www.rcsb.org/pdb/home/home.do) 网页 下载 PDB 文本的结构。 下载结构之后，可以利用 VMD，chimera, PyMOL 等可视化程序看一下蛋白质 结构。看了这个分子之后，你要去掉结晶水。用纯文本编辑器，比如 vi, emacs(Linux/Mac)或者 notepad(Windows)。不要用文字处理软件！删掉那些相 关行(PDB 文件中 residue“HOH”). 注意这个过程不是必须的(比如在水分子的 活性部位结合案例中)。我们强调我们这里不需要结晶水。 检查你的.pdb 文件看 MISSING 下面列的条目。因为这些条目显示了不存在 于该晶体结构中的原子或者残基。终端可能没有或者不显示关于动力学的问题。 任何不完整的内部序列或氨基酸残基的原子的缺失将导致 pdb2gmx 失败。这些丢 失的原子/残基必须通过其他软件来弥补。还要注意 pdb2gmx 不是万能的，不能 为任意分子产生拓扑文件，而只针对力场库中已经定义的残基(在*.rtp 文件中 的一般的蛋白质，核酸和一些有限的辅助因子，像是 NAD(H) 和 ATP). 现在结晶水已经没了，我们需要验证那些需要的原子还都在，PDB 文件应该 只含有蛋白质原子，而且准备好了迎接第一个 gromacs 工具 pdb2gmx. Pdb2gmx 的目的是产生三个文件。 1. 分子的拓扑文件 2. 位置限定文件 3. 后处理文件 拓扑文件(默认为 topol.top)包含定义分子模拟的所有信息。这些信息包含 非成键参数(原子类型和电荷)和成键信息(键长，键角和二面角)。拓扑文件形成 后，我们要仔细看一下。 用如下命令执行 pdb2gmx 命令： pdb2gmx -f 1AKI.pdb -o 1AKI_processed.gro -water spce 结构将被 pdb2gmx 处理，你要选择力场： Select the Force Field: From '/usr/local/gromacs/share/gromacs/top': 1: AMBER03 force field (Duan et al., J. Comp. Chem. 24, 1999-2012, 2003) 2: AMBER94 force field (Cornell et al., JACS 117, 5179-5197, 1995) dell 高亮 dell 高亮 yo ng ma 20 08 @小 木 虫 3: AMBER96 force field (Kollman et al., Acc. Chem. Res. 29, 461-469, 1996) 4: AMBER99 force field (Wang et al., J. Comp. Chem. 21, 1049-1074, 2000) 5: AMBER99SB force field (Hornak et al., Proteins 65, 712-725, 2006) 6: AMBER99SB-ILDN force field (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 (with CMAP) - version 2.0 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: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals) 15: [DEPRECATED] Encad all-atom force field, using full solvent charges 16: [DEPRECATED] Encad all-atom force field, using scaled-down vacuum charges 17: [DEPRECATED] Gromacs force field (see manual) 18: [DEPRECATED] Gromacs force field with hydrogens for NMR 力场包含将要写入拓扑文件的信息。这是个非常重要的选择。你必须完全了 解每个力场并且决定哪个力场对你的研究对象是最适合的。这个教程中，我们选 用全原子 OPLS 力场，在命令提示行输入 14，然后按 Enter 键。 还有些 pdb2gmx 会用的其他选项，下面列了几个： • -ignh: 忽略 PDB 文件中的氢原子，特别是对 NMR 结构有用。否则如果氢 原子存在，他们必须有正确的顺序和名称，正如 Gromacs 所期望的那样。 • -ter：为 N-和 C-末端交互式分配电荷态 • -inter：为谷氨酸，天门冬氨酸，赖氨酸，精氨酸，和 His 交互式指定电 荷态；分配二硫键到胱氨酸 现在你已经得到了三个文件，1AKI_processed.gro, topol.top, 和 posre.itp. 1AKI_processed.gro 是个 gromacs 格式结构文件，包含立场中定义 的所有原子(氢原子已经被加到蛋白质中的氨基酸上了)。 topol.top 文件是系 统拓扑文件(稍后会解释)。 posre.itp 文件包含用来限制位置的信息(后面解 释)。 第二步，检查拓扑文件 我们来看一下拓扑输出文件(topol.top)中有什么。再次使用纯文本编辑器， 检查其内容。在几句注释行之后(前面用分号；标注)，你会找到下面的语句： #include "oplsaa.ff/forcefield.itp" yo ng ma 20 08 @小 木 虫 此行调用在 OPLS- AA 的力场参数，这是这个文件的开头部分，表示所有下 面的参数都是从这个力场中来的。下面一个重要行是[ moleculetype ]，在那下 面你会找到： ; Name nrexcl Protein_A 3 “Protein A”定义了分子名字，基于这个蛋白质在 PDB 文件中被标定为 A 链。 它对成键的邻居有三个排斥。关于排斥的更多的信息可以从 gromacs 的手册上找 到。这部分内容的讨论不是本教材的范围。 下部分定义了蛋白质中的[ atoms ]，信息是如下几列： [ atoms ] ; nr type resnr residue atom cgnr charge mass typeB chargeB massB 1 opls_287 1 LYSH N 1 -0.3 14.0067 ; qtot -0.3 2 opls_290 1 LYSH H1 1 0.33 1.008 ; qtot 0.03 3 opls_290 1 LYSH H2 1 0.33 1.008 ; qtot 0.36 4 opls_290 1 LYSH H3 1 0.33 1.008 ; qtot 0.69 5 opls_293B 1 LYSH CA 1 0.25 12.011 ; qtot 0.94 6 opls_140 1 LYSH HA 1 0.06 1.008 ; qtot 1 该信息的解释如下： • nr:原子序数 • type：原子类型 • resnr: 氨基酸残基序数 • residue: 氨基酸残基名 注意这里的残基在 PDB 文件中原来是 LYS, 变成 LYSH 意味着这个残基被质子化 了。 • atom：原子名称 • cgnr：电荷组数 电荷组定义了整数电荷单位，旨在加速计算。 • charge: 无需解释 dell 矩形 yo ng ma 20 08 @小 木 虫 Qtot 描述保持对分子总电荷的持续计数 • mass: 也无需解释 • typeB,chargeB,massB：用于自由能微扰(这里不讨论) 下面的项，包括 [ bonds ], [ pairs ], [ angles ], and [ dihedrals ]. 这些项都是无需解释的。这些项中的参数和函数类型在 gromacs 手册的第 5章中 阐述。特别的，1-4 相互作用包含在“pairs“(gromacs 手册中 5.3.4 部分) 从位置限制拓扑文件开始，提示部分涉及了一些有用的/必须的拓扑文件。 "posre.itp"是从 pdb2gmx 命令产生的，它定义了一个力常数用于做平衡时能保 持原子不动(后面会很多涉及)。 ; Include Position restraint file #ifdef POSRES #include "posre.itp" #endif 这里结束了蛋白质 A分子类型的定义。拓扑文件的提示部分用来定义其他分 子并提供系统层次的描述。下面的分子类型(默认)是溶液，在这里是 SPC/E 水。 其他的对水的选择包括 SPC, TIP3P, 和 TIP4P。我们通过在 pdb2gmx 命令中使用 "-water spce"选这个水模型。对于其他不同水模型的详细描述，可以点这里 (http://www.lsbu.ac.uk/water/models.html) ，但是注意 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 正如你看到的，水分子也可以被位置限制，通过利用一个力常数(kpr) 为 1000 kJ mol-1 nm-2. 离子的信息也别包含在内： ; Include generic topology for ions #include "oplsaa.ff/ions.itp" 最后是系统层次定义： [ system ] ; Name yo ng ma 20 08 @小 木 虫 LYSOZYME [ molecules ] ; Compound #mols Protein_A 1 [ system ] 直 接 给 出 模 拟 中 将 要 写 到 输 出 文 件 中 的 系 统 名 字。 [ molecules ]直接列出系统中的所有分子(直接从[ moleculetype ]的名 字中列出)，为了他们在坐标文件(.gro)中出现。 现在我们检查完了拓扑文件中的内容，我们可以继续组建我们的系统了。 第三步，定义单位胞 & 添加溶剂 现在你对 gromacs 的拓扑文件比较熟悉了，到了继续创建系统的时候了。这 个例子中，我们将要模拟一个非常简单的水系统。模拟蛋白质或者其他分子在不 同溶液中也是可以的，只要涉及的分子的参数是合适的就行。 定义盒子和添加溶液有两步： 1. 用 editconf 定义盒子维度 2. 用 genbox 往盒子里填水 您现在面临一个如何处理的晶胞选择。对于本教程的目的，我们将使用一个 简单的立方体晶胞。当你对周期性的边界条件和箱型更了解时，我强烈推荐菱形 十二面体，因为在同一周期距离下，其体积为立方体的 71％，因此可以减少需 要加入的熔剂水分子的数量。 用 editconf 定义箱子： editconf -f 1AKI_processed.gro -o 1AKI_newbox.gro -c -d 1.0 -bt cubic 上面的命令把蛋白质放在盒子中心(-c), 并且把它放在距离盒子边缘 1.0 纳米位置(-d 1.0). 盒子类型是立方体(-bt cubic)。 盒子边缘的距离是个重要 参数。 因为我们要用到边界条件，我们必须满足最低限度的图像公约，那就是 一个蛋白质永远不能看到他的周期图像，否则计算的力是伪力。指定一个距离是 1.0 纳米的溶剂盒子意味着两个周期性蛋白质图像的距离至少是 2.0 纳米。这个 距离对于模拟中经常用到的截止模型都是足够的。 定义了盒子之后，我们要填充溶质(水)了。用 genbox 添加溶质： genbox -cp 1AKI_newbox.gro -cs spc216.gro -o 1AKI_solv.gro -p topol.top 蛋白质的配置(-cp)包含在前面一步 editconf 的输出文件中，溶剂的配置 (-cs)是 gromacs 安装的一部分。我们要用 spc216.gro，它是个通用的 3 点 平衡溶剂模型。你可以用为 SPC, SPC/E, or TIP3P 水模型。使用 spc216.gro dell 高亮 dell 高亮 yo ng ma 20 08 @小 木 虫 作为溶液配置，因为他们都是三点型的水模型。输出文件是 1AKI_solv.gro 我 们告诉 genbox 拓扑文件的名字(topol.top)，所以它能修改。注意 topol.top 中[ molecules ] 的变化： [ molecules ] ; Compound #mols Protein_A 1 SOL 10832 genbox 做的是追踪了增加了多少水分子，这个数字会写入你的拓扑文件显 示它已经做了。注意如果你用其他的(非水)溶剂，genbox 将不会在拓扑文件中 写入信息。它与水分子的兼容性更新是硬编码 第四步，添加离子 现在我们有了包含一个带电荷的蛋白质的溶液盒子。pdb2gmx 的输出文件告 诉我们这个蛋白质包含+8e 的静电荷(由于它的氨基酸组成)。如果你忽略了 pdb2gmx 输出文件的这个信息，可以看看你的 topol.top 文件中的最后一行，它 将显示“qtot 8.”. 因为生命中不存在静电荷，所以我们需要往我们的系统里 面添加离子。 gromacs程序中添加离子的命令是genion.genion做的是从拓扑文件中读取 数据并用你选用的离子去替换水分子。这个输入文件叫做运行输入文件，拓展名 是.tpr. 这个文件是用 gromacs 的 grompp(Gromacs Pre-Processor) 工具产生 的，它也将在我们随后运行第一次模拟时用到。grompp 做的是处理坐标文件和 拓扑文件(定义分子)来产生一个原子层次的输入文件(.tpr)。.tpr 文件 包含了 系统中所有原子的参数。 为了用 grompp 产生.tpr 文件我们还需要另外一个拓展名为.mdp (molecular dynamics parameter file)的输入文件，grompp 将组合 mdp 文件中 设定的参数和结构和拓扑文件来生成.tpr 文件。 .mdp 文件在能量最小化或者分子动力学模拟中常用的文件，但在这里只是简单 的用来生成一个系统中原子的描述。一个.mdp 例子(我们将用)可以在这里下载。 (http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutoria ls/lysozyme/Files/ions.mdp) 实际上，这里用的.mdp 文件可以包含任何参数的合理组合。我通常用一个 能量最小化脚本，因为它非常基础而且不用涉及任何复杂的参数组合。 下面的命令将组合你的.tpr 文件。 grompp -f ions.mdp -c 1AKI_solv.gro -p topol.top -o ions.tpr yo ng ma 20 08 @小 木 虫 现在我们有一个在原子层次描述我们系统的.tpr 文件。这个文件将用在 genion 工具中： genion -s ions.tpr -o 1AKI_solv_ions.gro -p topol.top -pname NA -nname CL -nn 8 提示行出现后，选 12“SOL”来嵌入离子。你肯定不想用离子去换掉蛋白质 的某一部分。 在 genion 命令中，我们提供了一个结构/态文件(-s)作为输入文件，产生一 个.gro 作为输出文件(-o), 处理拓扑文件(-p)来反映拿掉的水分子和增加的离 子，定义了阳离子和阴离子的名字(分别是-pname 和 -nname)， 并且告诉 genion 仅添加必须的离子来中和蛋白质中的静电荷，添加正确的负离子数目(-nn 8). 除了用-neutral 中和系统电荷你也可以添加-conc 组合起来，用 genion 添加特 定的浓度的离子。关于怎么用这些选项，可参考 genion 的手册页。 Gromacs 前期版本中离子名称是随力场而定的，到了 4.5 版本之后就完全标 准化了。原子名字一直是元素符号的大写字母，与[ moleculetype ]中一致。残 基名称可能带或者不带电荷的符号(+/-). 如果你遇到什么困难你可以参考 ions.itp 中的正确名字。 你的 [ molecules ] 指令现在看起来就像这样： [ molecules ] ; Compound #mols Protein_A 1 SOL 10824 NA 0 CL 8 这里的“NA 0”表示系统里没有钠离子，你可以安全的删除这行或者留着。 前期 gromacs 版本会报错。Gromacs-4.0 和以后版本允许 0 分子存在。 yo ng ma 20 08 @小 木 虫 第五步，能量最小化 加了溶剂的电中性的系统现在已经做好了。我们开始做动力学之前，我们必 须保证我们的系统里面没有立体几何冲突或不当。 结构要通过一个叫能量最小 化的过程弛豫一下。 能量最小化过程就跟添加离子差不多。我们再一次使用 grompp 来把结构， 拓扑和模拟参数组合到一个二进制输入文件里面(.tpr), 但是这次我们要通过 gromacs 的引擎 mdrun，而不是 genion 来进行能量最小化。 通过 grompp 用这个参数的文件，组合二进制输入文件。 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutoria ls/lysozyme/Files/minim.mdp grompp -f minim.mdp -c 1AKI_solv_ions.gro -p topol.top -o em.tpr 确保你在运行 genbox 和 genion 时已经更新了你的 topol.top 文件，否则你 将得到一堆错误信息(坐标文件中的坐标跟拓扑文件不匹配，等等). 现在我们做好了调用 mdrun 来运行能量最小化的准备： mdrun -v -deffnm em 使用- v 标志因为我们没什么耐心：它使 mdrun 冗长，使它在屏幕上打印每 一步的进展情况。使用- deffnm 标志将定义输入和输出的文件名。所以如果你 没有名字，你的 grompp 就输出“em.tpr”，你必须明确指定其与 mdrun -s 标志 的名称。就我们而言，我们将得到以下文件： • em.log: ASCII-文本的能量最小化过程日志文件 • em.edr: 二进制能量文件 • em.trr: 二进制高精度轨迹 • em.gro: 能量最小化结构 有两个重要的因素来评价能量最小化是否成功。第一个是势能(在能量最小 化过程的最后面输出，即使你没用-v 标志)。Eopt 应该是负值，根据系统大小和 水分子的多少，大约在 105-106数量级(对水中的蛋白质而言)。第二个重要的指 标是力的最大值。Fmax，在 minim.mdp 中设置的目标是“emtol=1000.0“表示 Fmax不能大于1000 kJ mol-1 nm-1。有可能当Fmax>emtol时得到一个合理的Eopt。 如果出现这样的问题，你的系统对于模拟来讲可能不够稳定。评价一下为什么会 这样，可能要改一下你的关于能量最小化的参数设置(integrator, emstep 等 等) 我们来做一些分析。em.edr 文件中包含 gromacs 在能量最小化过程中得到 的所有的能量。你可以用 gromacs 里面的 g_energy 来分析任何一个.edr 文件： yo ng ma 20 08 @小 木 虫 g_energy -f em.edr -o potential.xvg 在提示行中，输入”10 0”来选择 potential(10),zero(0)来结束输入。 Eopt 的平均值将被显示，写入一个叫做 potential.xvg 的文件。想画出这些数 据，你需要用 Xmgrace 描点工具。结果中的点看起来就像这样，展示一个漂亮的 稳定的 Eopt 收敛曲线。 现在我们的系统已经在能量最低点，我们可以开始真正的动力学模拟了。 第六步，平衡 EM 保证我们在结构和溶液取向方面有一个合理的初始结构。为了开始真正 的动力学模拟，我们必须对蛋白质周围的溶剂和离子做一些平衡。如果我们要在 这个点来尝试非限制性动力学，系统将会崩溃。原因是溶剂只在自己内部优化还 没考虑溶质。我们想要模拟的温度也要引入来确定关于溶质(蛋白质)的合适的取 向。达到正确的温度(取决于动能)之后，我们要对系统施加压力知道它达到合适 的密度。 还记得我们好久前用pdb2gmx生成的posre.itp文件么，现在要派上用场了。 Posre.itp 文件的目的就是对蛋白质中的重原子(非氢原子)实施位置限制力。除 非经过大量的能量惩罚，否则是不能动的。位置限定的用途就是在没有引起蛋白 质结构变化的前提下，允许我们来平衡蛋白质周围的溶剂。 平衡往往是分两个阶段进行。第一阶段在一个 NVT 系统（固定数目的粒子， 体积和温度）下进行。这个系统也被称为“等温等容”或“规范”。这种程序的 时间表与该系统的内容相关，但在 NVT 系，系统的温度应达到预期值的平台。如 果温度还没有稳定下来，将需要更多的时间。通常情况下，50 至 100 ps 的就足 够了，我们会为这项工作进行了 100 ps 的 NVT 系平衡。根据您的机器，这可能 需要(对双核的 MacBook 刚刚超过一小时)一段时间。从这里取得.mdp 文件 dell 高亮 yo ng ma 20 08 @小 木 虫 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorial s/lysozyme/Files/nvt.mdp 我们将要像在能量最小化过程中一样来调用 grompp 和 mdrun grompp -f nvt.mdp -c em.gro -p topol.top -o nvt.tpr mdrun -deffnm nvt 对采用的参数的全部解释可以从 gromacs 手册找到，包括一些解释。注意一 些.mdp 文件中的参数： • gen_vel = yes: 产生初始速度。利用不同的随机种子(gen_seed)，提供 不同的初始速度,因此多个(不同的)模拟可以从一个相同的初始结构开 始. • tcoupl = V-rescale: 速度重新调整恒温器是对 Berendsen 弱耦合方法的 改进，该方法不能产生正确的动能。 • pcoupl = no: 没有用压力耦合 再次用 g_energy 来分析温度进展： g_energy -f nvt.edr 在提示行输入“15 0”来选择系统温度并退出。结果显示： 从这些点中可以看出，系统温度很快就达到了目标温度(300K)，对其余的平 衡保持稳定。对于这个系统，较短的平衡时间（50 ps）已经足够。 第七步，平衡，第二部分 dell 高亮 yo ng ma 20 08 @小 木 虫 先前那步，NVT 平衡，稳定了系统的温度。采集数据之前，我们还要稳定系 统的压强（当然还有密度）。压强平衡是在 NPT 条件下进行集成，其中,粒子，压 强，和温度都不变。该条件也被称为“等温等压”，与实验条件最相近。 100ps 的 NPT 平 衡 的 .mdp 文 件 可 以 在 这 里 (http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutoria ls/lysozyme/Files/npt.mdp) 找到。该文件与 NVT 平衡时用的 mdp 文件没有太 大的不同。注意添加的压强耦合项，使用 Parrinello-rahman 恒压器。 几项其他的变动： • continuation = yes: 我们将从 NVT 平衡状态继续模拟 • gen_vel =速度从轨迹中读取 我们将调用在 NVT 平衡时用过的 grompp 和 mdrun。 注意我们现在要用到-t 来包含 NVT 平衡过程中的检查点文件。这个文件包含了继续模拟所需要的变量。 我们需要包含这个文件来保存在 NVT 过程中得到的速度。坐标文件(-c) 是 NVT 模拟的最后文件。 grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p topol.top -o npt.tpr mdrun -deffnm npt 我们来分析压强进展，再次用到 g_energy. g_energy -f npt.edr -o pressure.xvg 在提示行输入“15 0”来选择系统压强并退出。结果如下图所示： dell 高亮 dell 高亮 dell 高亮 yo ng ma 20 08 @小 木 虫 在 100ps 的平衡阶段，压强值有较大的波动，这并不意外。数据的平均值在 图中用红线标示。在该平衡过程中，压强平均值是 1.05bar. 我们再来看看密度，这次在 g_energy 的提示行输入“21 0” g_energy -f npt.edr -o density.xvg 跟压强一样，运行得到的密度平均值也用红线标示。100ps 过程中的密度平 均值为 998.3 kg m-3, 与实验值 1000 kg m-3接近。SPC/E 水模型的参数比较接近 的克隆了水的实验值。该密度值比较稳定，意味着系统在这个压强和密度下平衡 良好。 第八步，进行 MD 随着两个平衡阶段的完成，系统现在已经在需要的温度和压强下平衡良好。 我们已经准备好了进行位置限制和并为取得数据而进行 MD 了。这个过程跟我们 以前做的类似。我们还要用到检查点文件(该情况下，它包含了压强耦合信息) 来 进行 grompp。 我们要进行一个 1ns 的 MD 模拟，脚本文件可以从这里 (http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutoria ls/lysozyme/Files/md.mdp) 得到。 grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr mdrun -deffnm md_0_1 Mdrun 步骤最好在并行计算机计算，因为它要花上几小时才能完成。执行这 个任务合适的命令是： mpirun -np X mdrun_mpi -deffnm md_0_1 yo ng ma 20 08 @小 木 虫 其中 X表示模拟所需要用的处理器数目。 第九步，分析 现在我们已经模拟了我们的蛋白质，我们应该来分析一下我们的系统。什么 类型的数据才是重要的呢?这是模拟前问的一个重要问，所以你应该对你需要采 集的系统的数据类型有自己的想法。在本教程中，我们介绍了一些基本工具。 第一个是 trjconv，这是一个后期处理工具，用来处理坐标，纠正周期性或 者手动来调整轨迹(时间单位，帧频率等)。对于本练习，我们要用 trjconv 来核 算系统中的周期性。蛋白质可能会在单元格中扩散，也可能“跳”在两个盒子之 间。我们采用下面命令来处理这种情况。 trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -pbc mol -ur compact 选择 0("System")。我们要对这个“修正”了的轨迹进行我们的分析。先来 看看结构稳定性。GROMACS 有个内置工具进行 RMSD 计算，名字叫 g_rms。用这个 命令来使用这个工具： g_rms -s md_0_1.tpr -f md_0_1_noPBC.xtc -o rmsd.xvg -tu ns 对最小二乘拟合和群计算 RMSD，都选择 4("Backbone")。-tu 标志将输出结 果以 ns 显示，即使轨迹文件是用 ps 写的。这是为了输出文件更加清晰(尤其是 你做一个较长的 105ps 的模拟时)。输出文件将显示在最小化的，平衡了的系统 的结构下的 RMSD. 如果我们要计算晶体结构的 RMSD 值，我们需要用到下面命令： yo ng ma 20 08 @小 木 虫 g_rms -s em.tpr -f md_0_1_noPBC.xtc -o rmsd_xtal.xvg -tu ns 结果将如下所示： 两个数据显示 RMSD 大约在 0.1nm，表示该结构非常稳定。两者的微小不同 在于当 t=0ns 时晶体结构是不同的。这是可以预料的，因为它已经被能量最小化 了，而且如我们前面讨论的，位置限定不是 100%完美的。 蛋白质的回转半径可衡量其紧凑度。如果蛋白质稳定的折叠了，Rg 将保持 一个相对温度值。如果蛋白质展开了，它的 Rg 将随时间变化。我们来分析一下 我们模拟的溶菌酶的回转半径。 g_gyrate -s md_0_1.tpr -f md_0_1_noPBC.xtc -o gyrate.xvg yo ng ma 20 08 @小 木 虫 从这个合理的恒定的 Rg 值我们可以看到，在温度为 300k 的 1 纳秒时间内蛋 白质的紧凑型形式（折叠）非常稳定。这一结果并不出人意料，而这足以说明 GROMACS 较高的分析能力。 现在已经用 gromacs 操作了一个分子动力学模拟过程，并分析了一些结果， 本教程不是个全面的东西。用 gromacs 你可以有很多类型的模拟可以操作(自由 能计算，分平衡 MD,和正常模式分析还有一些其他的)。你应该看一些文献和 gromacs 手册来调整这里提供的.mdp 文件来达到更有效地和更精确的目标。 如果你对改善这个教程有什么建议，如果你发现了错误，或者什么东西不清 楚，请不要客气给我发邮件(jalemkul@vt.edu)。注意：这不是邀请你为了任何 gromacs 的问题而给我发邮件。我没有为了一个私人的教程或者个人帮助服务而 为自己打广告。 享受你的模拟吧！