1、Gromacs 中文教程淮海一粟分子动力学（MD ）模拟分为三步：首先，要准备好模拟系统；然后，对准备好的系统进行模拟；最后，对模拟结果进行分析。虽然第二步是最耗费计算资源的，有时候需要计算几个月，但是最耗费体力的步骤在于模拟系统准备和结果分析。本教程涉及模拟系统准备、模拟和结果分析。一、数据格式处理准备好模拟系统是 MD 最重要的步骤之一。MD 模拟原子尺度的动力学过程，可用于理解实验现象、验证理论假说，或者为一个待验证的新假说提供基础。然而，对于上述各种情形，都需要根据实际情况对模拟过程进行设计；这意味着模拟的时候必须十分小心。丢失的残基、原子和非标准基团本教程模拟的是蛋白质。首先需要找到
2、蛋白质序列并选择其起始结构，见前述；然后就要检查这个结构是否包含所有的残基和原子，这些残基和原子有时候也是模拟所必需的。本教程假定不存在缺失，故略去。另一个需要注意的问题是结构文件中可能包含非标准残基，被修饰过的残基或者配体，这些基团还没有力场参数。如果有这些基团，要么被除去，要么就需要补充力场参数，这牵涉到 MD 的高级技巧。本教程假定所有的蛋白质不含这类残基。结构质量对结构文件进行检查以了解结构文件的质量是一个很好的练习。例如，晶体结构解析过程中，对于谷氨酰胺和天冬酰胺有可能产生不正确的构象；对于组氨酸的质子化状态和侧链构象的解析也可能有问题。为了得到正确的结构，可以利用一些程序和服务器(
3、如 WHATIF)。本教程假定所用的结构没有问题，我们只进行数据格式处理。二、结构转换和拓扑化一个分子可以由各个原子的坐标、键接情况与非键相互作用来确定。由于.pdb 结构文件只含有原子坐标，我们首先必须建立拓扑文件，该文件描述了原子类型、电荷、成键情况等信息。拓扑文件对应着一种力场，选择何种力场对于拓扑文件的建立是一个值得仔细考虑的问题。这里我们用的是 GROMOS96 53a6 连接原子力场，该力场对于氨基酸侧链的自由能预测较好，并且与 NMR 试验结果较吻合。拓扑数据要与结构吻合，这点很重要；这意味着结构数据也需要在所用的力场中进行转换。为了转换结构并建立拓扑数据，可以用 pdb2gmx
4、 程序，该程序对特定结构块（如氨基酸）建立拓扑数据。它需要使用结构块库，如果结构里面有库中没有的结构单元，转换过程就会失败。输入以下命令来进行结构转换，选择 GROMOS 53a6 力场和 SPC 水分子模型. 选项 ignh 表示起始文件中的氢原子会先被除掉，然后在相应的力场中重建。pdb2gmx -f protein.pdb -o protein.gro -p protein.top -ignh You can always generate a pdb file from a gro file using the following command:editconf -f XXXXX.g
5、ro -o XXXXX.pdb仔细阅读屏幕提示，并检查对于组氨酸质子化所做的选择，注意蛋白质总电荷数；也要浏览输入文件（protein.pdb, pdb 格式）和输出文件（protein.gro, gromacs 格式）。注意两者的区别；最后浏览拓扑文件及其结构。三、结构的能量最小化（真空中）现在，一个适合于所选力场的正确结构文件已经建立好，同时还得到了相应的拓扑数据。然而，由于结构转换中，牵涉到氢原子的删除和重建，这可能会带来一些相互作用力的变化，比如由于两个原子的位置靠得太近引起的作用力。因此我们必须对结构进行一次能量优化，使之松弛一点。这其实就是一个模拟步骤，包括两个过程：第一步，结构和
6、拓扑数据被合成在一起，同时还包含了一些控制参数。这个过程对产生一个文件，该文件作为输入文件，用于第二步 mdrun 程序的动力学模拟。 把这个过程分为两步的好处是，输入文件可以传送到专门的（高性能）计算机网络或者超级计算机，在那里完成模拟计算。然而，一般只在正式动力学模拟中才会这么做。 为了执行第一步，需要用到一个.mdp 文件，该文件（文件名 minim.mdp）已经被做好并放在你们的目录下面。该文件包含一系列能量最小化的控制参数，请查看之。注意，文件以“integrator“行开始，表示指定的数学模型。本例中，被指定采用最速下降法计算5000 步。现在，用以下命令合并结构和拓扑文件，以及一
7、些控制参数：grompp -v -f minim.mdp -c protein.gro -p protein.top -o protein-EM-vacuum.tprgrompp 程序将会标明此体系存在非零电荷，并将写入有关系统和控制参数的其他信息。该程序还会产生一个额外的输出文件(mdout.mdp)，该文件包含所有的控制参数设置. 下一步，运行 mdrun：mdrun -v -deffnm protein-EM-vacuum -c protein-EM-vacuum.gro由于该蛋白只含有 1500 个原子，能量优化非常快。选项 -v 将把每步的势能和最大势能力打印出来，以便于跟踪能量优化
8、过程。-deffnm 选项设定了输出文件的默认文件名，而不需要对每个输出文件都进行命名，以减少选项设置。现在，结构松弛下来了，我们该设定周期性边界并添加溶剂。Which method was used for the energy minimization? How many steps were specified in the parameter file and how many steps did the minimization take? What could cause the energy minimization to stop before the specified nu
9、mber of steps was reached? What is the final potential energy of the system? Compare the structures before and after energy minimization by loading them into PyMol四、周期性边界条件设定在添加溶剂之前，需要选定模拟体系的大致外形。大多数常见的模拟都是在周期性边界条件(PBC)下进行的，这意味着要定义一个外形单元框，该框可用于空间填充堆积。这样，就可以对一个无限、周期性系统进行模拟，以避免由于边框造成的边界效应。建立PBC 时，仅有少数
10、几个通用的形状。我们将采用 rhombic dodecahedron，因为其对应于最优堆积空间，因此是自由旋转原子的最佳选择。为了不产生周期性影像的直接接触，我们设定了蛋白质距离边框的最小距离为 1.0，这样的话，两个相邻的堆积单元的距离就会大于 2.0 nm。周期性边界条件用 editconf 命令设定：editconf -f protein-EM-vacuum.gro -o protein-PBC.gro -bt dodecahedron -d 1.2看看 editconf 的输出，注意体积的变化。同时，也看看文件 protein-PBC.gro 的最后一行。在 gromacs 格式文件
11、(.gro)中，最后一行表示溶剂盒子的形状。我们常常用三斜晶系矩阵表示法，其中前三个数字表示对角线元素 (xx, yy, zz)，最后六个表示非对角线元素(xy, xz, yx, yz, zx, zy)。What is the new unit cell volume? 五、溶剂条件设定溶剂盒子设定好了之后，就可添加溶剂了。溶剂模型有好几种，每种模型多少都与一种力场相对应。GROMOS96 力场常用于 SPC 水分子模型。溶剂添加不需要写入拓扑数据，但是也可以在拓扑数据中包含溶剂模型。向盒子中添加溶剂用以下命令：genbox -cp protein-PBC.gro -cs spc216.gro
12、 -p protein.top -o protein-water.gro看看 protein.top 的结尾溶剂添加的情况，添加了多少溶剂？How many water molecules are added to the system and what volume does that correspond to? How many atoms are in the system now? 用以前学过的命令，将 .gro 文件转换为 .pdb 文件。用 Pymol 程序读入溶剂化的结构文件 (protein-water.pdb)。Pymol 中，可用“show cell”命令显示盒子形状。s
13、how cell 这将画出一个框架线来显示盒子的三维空间边界。这个三斜晶图形对应菱形十二面体，但可能不太明显。另一个值得注意的事情是，并非所有的溶剂都在盒子内，但以砖形排布。非常类似地，蛋白质也有部分伸出水外面。 Why is it not a problem to have the protein sticking out of the water? 六 添加相反的离子以平衡体系电荷目前我们已经有了溶剂化的蛋白质了，但是体系的净电荷数不为零。为了使体系电中性，我们需要添加一定数量的相反离子。添加一定离子的步骤，是一个很好的练习，程序genion 能完成这些任务，但是它需要一个输入文件，例如，
14、一个包含了结构和拓扑数据的文件。就像前面一样，我们可以用 grompp 来产生这样的文件：grompp -v -f minim.mdp -c protein-water.gro -p protein.top -o protein-water.tpr文件 protein-water.tpr 既可作为 genion 程序的输入文件。我们设定了 NA+/CL-(-pname NA+ -nname CL-)的浓度为 0.15 M (-conc 0.15)，并指定特定离子的加入以是体系电中性 (-neutral)。Genion 程序将询问用离子取代何种分子，选择 SOL。genion -s protei
15、n-water.tpr -o protein-solvated.gro -conc 0.15 -neutral -pname NA+ -nname CL-注意添加离子的数量，并注意一个额外的氯离子被加入，以中和体系电荷。通过将水分子置换为离子，体系的拓扑数据文件（protein.top）将不正确了（你可以修改拓扑文件，减少溶剂分子数目；同时加入一行命令，说明添加 Na 离子的数量并写入另一行命令说明Cl 离子的数量）。Edit the topology file and decrease the number of solvent molecules. Also add a line spec
16、ifying the number of NA ions and a line specifying the amount of CL. How many sodium and chloride ions are added to the system?七、溶剂化体系的能量最小化到此为止，模拟系统已经准备好了。在进行正式模拟前，还需要使体系再次松弛。因为溶剂的添加以及离子的置换，可能产生不利的相互作用力，例如，原子重叠、同种电荷的接近，等，因此需要对体系再次能量优化，其步骤与前面相同：先运行 grompp，再运行 mdrun：grompp -v -f minim.mdp -c protein-
17、solvated.gro -p protein.top -o protein-EM-solvated.tpr mdrun -v -deffnm protein-EM-solvated -c protein-EM-solvated.gro How many steps were specified in the parameter file and how many steps did the minimization take? What is the final potential energy of the system?八、溶剂和氢原子位置的松弛：位置限制性 MD为了缓解体系的随机分布的
18、张力，我们是溶剂与蛋白质相适应，比如，我们使溶剂自由移动，而使得蛋白质上面的非氢原子或多或少地固定在相对位置。这么做的目的是为了 确保溶剂的构象“适合”蛋白质。这一步是正式 MD 的第一步。这一步的控制参数在pr.mdp 文件中，执行时需要被读入。看看这个文件， 注意 integrate 行和 define 行。后者用于允许拓扑文件的流动控制（to allow flow control in the topology file）： define = -DPOSRES 将导致定义关键字 “POSRES“。查看拓扑文件的最后一行，看看怎样进行位置限制。采用 grompp 和 mdrun 程序运行模
19、拟：Remark: edit the file pr.mdp and change the value of gen_seed grompp -v -f pr.mdp -c protein-EM-solvated.gro -p protein.top -o protein-PR.tpr mdrun -v -deffnm protein-PR What is the length of the simulation in picoseconds? How is the inclusion of the position restraints handled in the topology fil
20、e? 你会很有趣地看到总能量、势能和动能的变化。在前一步，粒子没有速度，所以没有动能，因此也没有温度。现在，模拟一开始，原子就被赋予速度因此获得动能。模拟过程的能量信息被保存为不可读的二进制文件格式（.edr 文件）。该文件的信息可用 gromacs工具 g_energy 程序提取。注：在下一步模拟运行时，你能用以下命令并在执行过程中回答有关问题，以节省时间。如果这一步不奏效，让该程序继续在这里运行，而打开另外一个终端，运行下一步的grompp 和 mdrun 程序，但要注意两个终端的都要在正确的文件夹下运行。g_energy -f protein-PR.edr -o thermodynami
21、cs-PR.xvg你将要在一列能量参数中进行选择，选择好的参数结果将会被输出。输入与温度、势能、动能和总能量相对应的数字，然后输入 0（零），回车。输出文件 (.xvg) 可以用 xmgr 或 xmgrace 程序查看，这些程序需要预先安装好。也可以用 Python 语言脚本 xvg2ascii.py 查看，它能在终端窗口上显示出图线。xmgrace -nxy thermodynamics-PR.xvg注：绿色曲线是总能量，黑色是势能，蓝色是温度，红色是动能。你也可以点击曲线并填写“图标”窗口以改变图标。别忘了点击 accept 以接受改变。What happens with the temp
22、erature? What happens to the potential/kinetic/total energy and how can this be explained?九、非限制性 MD，打开温度耦合系统现在应该相对松弛了一些，我们就引入温度并开始对系统进行热浴耦合吧。换句话说，将在打开温度耦合的情况下运行一小段时间，让系统松弛到一个新的状态。下载控制参数文件 nvt.mdp，然后检查一下其中的温度耦合参数；最后用 grompp 和 mdrun 运行模拟：grompp -v -f nvt.mdp -c protein-PR.gro -p protein.top -o protein
23、-NVT.tpr mdrun -v -deffnm protein-NVTThere are two groups which are separately coupled to a heat bath. Which groups are that and what do you think is in each of these groups? 运行结束后，查看一下能量和温度。再次提醒，你可以在查看能量的同时，运行下一步 。数据提取方法同前：g_energy -f protein-NVT.edr -o thermodynamics-NVT.xvg观察能量图。What happens to t
24、he temperature? What happens to the potential/kinetic/total energy and how can this be explained? 十、非限制性 MD 模拟：打开压力耦合当引入温度并松弛后，就该引入压力了。导入控制参数文件 npt.mdp，短时间运行有压力下的模拟。先查看一下此控制参数文件，找到控制引入压力的句柄。此处也应该注意，因为我们需要重新赋予离子速度，通过 gen_vel 行实现。在这行下面，设定了速度分布（Maxwell 分布）的温度，并通过 gen_seed 产生初始值。 gen_seed 现在设定为 9999，但是会
25、被改为一个新的值。MD 模拟是最重要的一步，所以每次都以相同的坐标体系、速度和相同的参数开始模拟时，就会导致模拟的结果一模一样，这不是我们所希望的。当你设定好了控制参数后，运行以下命令进行模拟。Remark: edit the file npt.mdp and change the value of gen_seed grompp -v -f npt.mdp -c protein-NVT.gro -p protein.top -o protein-NPT.tpr mdrun -v -deffnm protein-NPT 运行完成后，再次查看能量和温度，同时注意观察压力变化。数据提取方法同前：g
26、_energy -f protein-NPT.edr -o thermodynamics-NPT.xvg看看能量曲线、温度和压力的变化。What happens to the temperature? What happens to the potential/kinetic/total energy and how can this be explained? 最后一个问题与从有限数量的粒子系统中提取热动力特性直接有关。粗略地说，热动力学指的是大量粒子（例如，几十亿个而不是几千个）的行为。取大量粒子的平均特性能够使数据波动减少，相反，仅对少量粒子进行平均，波动就较大。 由于我们引入了压力，系
27、统的密度会发生变化，所以从能量文件中提取密度数据，方法如下： g_energy -f protein-NPT.edr -o density-NPT.xvg How does the density behave over time? Why does the density of the system change if pressure coupling is turned on? 十一、正式模拟到了最后一步了。我们已经有了好赖平衡了的含溶剂系统，系统中包含我们感兴趣的蛋白质，所以该进行正式模拟了。记住，正式模拟并不表示整个模拟都可以用来分析感兴趣的特性。虽然一些初始参数值被擦除了，系统看起
28、来还不太可能就已经达到平衡状态了。在分析阶段，我们应当检查哪一部分模拟结果可以认为是在平衡状态下得到的，只有这部分数据才适合用于分析。先设定运行参数。这里只需要运行模拟步骤，类似于系统准备的第二步。然而，这一步也是需要考虑模拟目的的，因此控制参数里面应当包含能反映分析内容的有关参数。考虑以下问题： 在什么时间尺度上能够反映所研究的问题？需要运行多久？ 需要得到多少帧数据，时间解析度是多少？ 需不需要存储速度数据？ 是否需要输出所有的原子数据，还是只要有蛋白质坐标数据就够了？ 每隔多久记录一次能量文件和日志文件？ 每隔多久记录一次坐标和速度文件？ .我们将运行 10 纳秒的模拟。下载或拷贝参数文
29、件 md.mdp，先检查一下文件内容。下载或拷贝 md.mdp 文件并查看其中的内容。How many steps are needed to get a total simulation time of 10 ns?你要算一下，为了完成 10 纳秒的模拟，需要多少步数，并编辑好参数文件；然后用Grompp 命令将最终的结构文件和系统准备步骤中得到的拓扑文件合并为一个输入文件。 grompp -v -f md.mdp -c protein-NPT.gro -p protein.top -o topol.tpr虽然输入文件是二进制格式，我们还是可以看一看其内容。某些情况下，会发生一些意外，情况，其中有些意外情况可能与内部控制和力场参数有关，在这些情况下，这么做就尤为有用。利用 gmxdump 程序，可以读取这个文件。读出的文件可能会有很多页，建议这个时候用 more 或 less命令分屏显示文件。输入以下命令读取输入文件（注意，在实际操作时，需要把 topol.tpr 文件名换为运行 grompp 程序所得到的文件名）：gmxdump -s topol.tpr | less
