分子动力学模拟I

Gromacs中文教程

淮海一粟

分子动力学（MD）模拟分为三步：首先，要准备好模拟系统；然后，对准备好的系统进行模拟；最后，对模拟结果进行分析。虽然第二步是最耗费计算资源的，有时候需要计算几个月，但是最耗费体力的步骤在于模拟系统准备和结果分析。本教程涉及模拟系统准备、模拟和结果分析。

一、数据格式处理

准备好模拟系统是MD最重要的步骤之一。MD模拟原子尺度的动力学过程，可用于理解实验现象、验证理论假说，或者为一个待验证的新假说提供基础。然而，对于上述各种情形，都需要根据实际情况对模拟过程进行设计；这意味着模拟的时候必须十分小心。

丢失的残基、原子和非标准基团

本教程模拟的是蛋白质。首先需要找到蛋白质序列并选择其起始结构，见前述；然后就要检查这个结构是否包含所有的残基和原子，这些残基和原子有时候也是模拟所必需的。本教程假定不存在缺失，故略去。

另一个需要注意的问题是结构文件中可能包含非标准残基，被修饰过的残基或者配体，这些基团还没有力场参数。如果有这些基团，要么被除去，要么就需要补充力场参数，这牵涉到MD的高级技巧。本教程假定所有的蛋白质不含这类残基。

结构质量

对结构文件进行检查以了解结构文件的质量是一个很好的练习。例如，晶体结构解析过程中，对于谷氨酰胺和天冬酰胺有可能产生不正确的构象；对于组氨酸的质子化状态和侧链构象的解析也可能有问题。为了得到正确的结构，可以利用一些程序和服务器(如WHATIF)。本教程假定所用的结构没有问题，我们只进行数据格式处理。

二、结构转换和拓扑化

一个分子可以由各个原子的坐标、键接情况与非键相互作用来确定。由于.pdb结构文件只含有原子坐标，我们首先必须建立拓扑文件，该文件描述了原子类型、电荷、成键情况等信息。拓扑文件对应着一种力场，选择何种力场对于拓扑文件的建立是一个值得仔

细考虑的问题。这里我们用的是GROMOS96 53a6连接原子力场，该力场对于氨基酸侧

链的自由能预测较好，并且与NMR试验结果较吻合。

拓扑数据要与结构吻合，这点很重要；这意味着结构数据也需要在所用的力场中进行转换。

为了转换结构并建立拓扑数据，可以用pdb2gmx程序，该程序对特定结构块（如氨基酸）建立拓扑数据。它需要使用结构块库，如果结构里面有库中没有的结构单元，转换过程就会失败。输入以下命令来进行结构转换，选择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.gro -o XXXXX.pdb

仔细阅读屏幕提示，并检查对于组氨酸质子化所做的选择，注意蛋白质总电荷数；也要浏览输入文件（protein.pdb, pdb格式）和输出文件（protein.gro, gromacs格式）。注意两者的区别；最后浏览拓扑文件及其结构。

三、结构的能量最小化（真空中）

现在，一个适合于所选力场的正确结构文件已经建立好，同时还得到了相应的拓扑数据。然而，由于结构转换中，牵涉到氢原子的删除和重建，这可能会带来一些相互作用力的变化，比如由于两个原子的位置靠得太近引起的作用力。因此我们必须对结构进行一次能量优化，使之松弛一点。这其实就是一个模拟步骤，包括两个过程：第一步，结构和拓扑数据被合成在一起，同时还包含了一些控制参数。这个过程对产生一个文件，该文件作为输入文件，用于第二步mdrun程序的动力学模拟。把这个过程分为两步的好处是，输入文件可以传送到专门的（高性能）计算机网络或者超级计算机，在那里完成模拟计算。然而，一般只在正式动力学模拟中才会这么做。

为了执行第一步，需要用到一个.mdp文件，该文件（文件名minim.mdp）已经被做好并放在你们的目录下面。该文件包含一系列能量最小化的控制参数，请查看之。注意，文件以"integrator"行开始，表示指定的数学模型。本例中，被指定采用最速下降法计算5000步。现在，用以下命令合并结构和拓扑文件，以及一些控制参数：

grompp -v -f minim.mdp -c protein.gro -p protein.top -o protein-EM-vacuum.tpr

grompp程序将会标明此体系存在非零电荷，并将写入有关系统和控制参数的其他信息。该程序还会产生一个额外的输出文件(mdout.mdp)，该文件包含所有的控制参数设置. 下一步，运行mdrun：

mdrun -v -deffnm protein-EM-vacuum -c protein-EM-vacuum.gro

由于该蛋白只含有1500个原子，能量优化非常快。选项-v 将把每步的势能和最大势能力打印出来，以便于跟踪能量优化过程。-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 number 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时，仅有少数几个通用的形状。我们将采用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格式文件(.gro)中，最后一行表示溶剂盒子的形状。我们常常用三斜晶