分子动力学模拟

Gromacs中文教程

淮海一粟

分子动力学（MD）模拟分为三步：首先，要准备好模拟系统；然后，对准备好的系

统进行模拟；最后，对模拟结果进行分析。虽然第二步是最耗费计算资源的，有时候需要

计算几个月，但是最耗费体力的步骤在于模拟系统准备和结果分析。本教程涉及模拟系统

准备、模拟和结果分析。

一、数据格式处理

准备好模拟系统是MD最重要的步骤之一。MD模拟原子尺度的动力学过程，可用于

理解实验现象、验证理论假说，或者为一个待验证的新假说提供基础。然而，对于上述各

种情形，都需要根据实际情况对模拟过程进行设计；这意味着模拟的时候必须十分小心。

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

本教程模拟的是蛋白质。首先需要找到蛋白质序列并选择其起始结构，见前述；然后

就要检查这个结构是否包含所有的残基和原子，这些残基和原子有时候也是模拟所必需的。

本教程假定不存在缺失，故略去。

另一个需要注意的问题是结构文件中可能包含非标准残基，被修饰过的残基或者配体，

这些基团还没有力场参数。如果有这些基团，要么被除去，要么就需要补充力场参数，这

牵涉到MD的高级技巧。本教程假定所有的蛋白质不含这类残基。

结构质量

对结构文件进行检查以了解结构文件的质量是一个很好的练习。例如，晶体结构解析

过程中，对于谷氨酰胺和天冬酰胺有可能产生不正确的构象；对于组氨酸的质子化状态和

侧链构象的解析也可能有问题。为了得到正确的结构，可以利用一些程序和服务器(如

WHATIF)。本教程假定所用的结构没有问题，我们只进行数据格式处理。

二、结构转换和拓扑化

一个分子可以由各个原子的坐标、键接情况与非键相互作用来确定。由于.pdb结构

文件只含有原子坐标，我们首先必须建立拓扑文件，该文件描述了原子类型、电荷、成键

情况等信息。拓扑文件对应着一种力场，选择何种力场对于拓扑文件的建立是一个值得仔

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

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

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

转换。

为了转换结构并建立拓扑数据，可以用pdb2gmx程序，该程序对特定结构块（如氨

基酸）建立拓扑数据。它需要使用结构块库，如果结构里面有库中没有的结构单元，转换

过程就会失败。输入以下命令来进行结构转换，选择GROMOS53a6力场和SPC水分子

模型.选项–ignh表示起始文件中的氢原子会先被除掉，然后在相应的力场中重建。

仔细阅读屏幕提示，并检查对于组氨酸质子化所做的选择，注意蛋白质总电荷数；也

要浏览输入文件（,pdb格式）和输出文件（,gromacs格式）。

注意两者的区别；最后浏览拓扑文件及其结构。

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

现在，一个适合于所选力场的正确结构文件已经建立好，同时还得到了相应的拓扑数

据。然而，由于结构转换中，牵涉到氢原子的删除和重建，这可能会带来一些相互作用力

的变化，比如由于两个原子的位置靠得太近引起的作用力。因此我们必须对结构进行一次

能量优化，使之松弛一点。这其实就是一个模拟步骤，包括两个过程：第一步，结构和拓

扑数据被合成在一起，同时还包含了一些控制参数。这个过程对产生一个文件，该文件作

为输入文件，用于第二步mdrun程序的动力学模拟。把这个过程分为两步的好处是，输

入文件可以传送到专门的（高性能）计算机网络或者超级计算机，在那里完成模拟计算。

然而，一般只在正式动力学模拟中才会这么做。

为了执行第一步，需要用到一个.mdp文件，该文件（文件名）已经被做

好并放在你们的目录下面。该文件包含一系列能量最小化的控制参数，请查看之。注意，

文件以"integrator"行开始，表示指定的数学模型。本例中，被指定采用最速下降法计算

5000步。现在，用以下命令合并结构和拓扑文件，以及一些控制参数：

grompp程序将会标明此体系存在非零电荷，并将写入有关系统和控制参数的其他信

息。该程序还会产生一个额外的输出文件()，该文件包含所有的控制参数设置.

下一步，运行mdrun：

由于该蛋白只含有1500个原子，能量优化非常快。选项-v将把每步的势能和最大势

能力打印出来，以便于跟踪能量优化过程。-deffnm选项设定了输出文件的默认文件名，

而不需要对每个输出文件都进行命名，以减少选项设置。现在，结构松弛下来了，我们该

设定周期性边界并添加溶剂。

Whichmethodwasudfortheenergyminimization?

Howmanystepswerespecifiedintheparameterfileandhowmanystepsdidthe

minimizationtake?

Whatcouldcautheenergyminimizationtostopbeforethespecifiednumberof

stepswasreached?

Whatisthefinalpotentialenergyofthesystem?

Comparethestructuresbeforeandafterenergyminimizationbyloadingthem

intoPyMol

四、周期性边界条件设定

在添加溶剂之前，需要选定模拟体系的大致外形。大多数常见的模拟都是在周期性边

界条件(PBC)下进行的，这意味着要定义一个外形单元框，该框可用于空间填充堆积。这

样，就可以对一个无限、周期性系统进行模拟，以避免由于边框造成的边界效应。建立

PBC时，仅有少数几个通用的形状。我们将采用rhombicdodecahedron，因为其对应于

最优堆积空间，因此是自由旋转原子的最佳选择。为了不产生周期性影像的直接接触，我

们设定了蛋白质距离边框的最小距离为1.0，这样的话，两个相邻的堆积单元的距离就会

大于2.0nm。周期性边界条件用editconf命令设定：

-btdodecahedron-d1.2

看看editconf的输出，注意体积的变化。同时，也看看文件的最后

一行。在gromacs格式文件(.gro)中，最后一行表示溶剂盒子的形状。我们常常用三斜晶

系矩阵表示法，其中前三个数字表示对角线元素(xx,yy,zz)，最后六个表示非对角线元素

(xy,xz,yx,yz,zx,zy)。

Whatisthenewunitcellvolume?

五、溶剂条件设定

溶剂盒子设定好了之后，就可添加溶剂了。溶剂模型有好几种，每种模型多少都与一

种力场相对应。GROMOS96力场常用于SPC水分子模型。溶剂添加不需要写入拓扑数

据，但是也可以在拓扑数据中包含溶剂模型。向盒子中添加溶剂用以下命令：

-o

看看的结尾溶剂添加的情况，添加了多少溶剂？

Howmanywatermoleculesareaddedtothesystemandwhatvolumedoesthat

correspondto?

Howmanyatomsareinthesystemnow?

用以前学过的命令，将.gro文件转换为.pdb文件。用Pymol程序读入溶剂化的结构

文件()。Pymol中，可用“showcell”命令显示盒子形状。

showcell

这将画出一个框架线来显示盒子的三维空间边界。这个三斜晶图形对应菱形十二面体，

但可能不太明显。另一个值得注意的事情是，并非所有的溶剂都在盒子内，但以砖形排布。

非常类似地，蛋白质也有部分伸出水外面。

Whyisitnotaproblemtohavetheproteinstickingoutofthewater?

六添加相反的离子以平衡体系电荷

目前我们已经有了溶剂化的蛋白质了，但是体系的净电荷数不为零。为了使体系电中

性，我们需要添加一定数量的相反离子。添加一定离子的步骤，是一个很好的练习，程序

genion能完成这些任务，但是它需要一个输入文件，例如，一个包含了结构和拓扑数据

的文件。就像前面一样，我们可以用grompp来产生这样的文件：

-o

文件既可作为genion程序的输入文件。我们设定了NA+/CL-(-

pnameNA+-nnameCL-)的浓度为0.15M(-conc0.15)，并指定特定离子的加入以是体系

电中性(-neutral)。Genion程序将询问用离子取代何种分子，选择'SOL'。

-conc0.15-

neutral-pnameNA+-nnameCL-

注意添加离子的数量，并注意一个额外的氯离子被加入，以中和体系电荷。通过将水

分子置换为离子，体系的拓扑数据文件（）将不正确了（你可以修改拓扑文件，

减少溶剂分子数目；同时加入一行命令，说明添加Na离子的数量并写入另一行命令说明

Cl离子的数量）。

Howmanysodiumandchlorideionsareaddedtothesystem?

七、溶剂化体系的能量最小化

到此为止，模拟系统已经准备好了。在进行正式模拟前，还需要使体系再次松弛。因

为溶剂的添加以及离子的置换，可能产生不利的相互作用力，例如，原子重叠、同种电荷

的接近，等，因此需要对体系再次能量优化，其步骤与前面相同：先运行grompp，再运

行mdrun：

-o

Howmanystepswerespecifiedintheparameterfileandhowmanystepsdidthe

minimizationtake?

Whatisthefinalpotentialenergyofthesystem?

八、溶剂和氢原子位置的松弛：位置限制性MD

为了缓解体系的随机分布的张力，我们是溶剂与蛋白质相适应，比如，我们使溶剂自

由移动，而使得蛋白质上面的非氢原子或多或少地固定在相对位置。这么做的目的是为了

确保溶剂的构象“适合”蛋白质。这一步是正式MD的第一步。这一步的控制参数在

文件中，执行时需要被读入。看看这个文件，注意integrate行和define行。后

者用于允许拓扑文件的流动控制（toallowflowcontrolinthetopologyfile）：define=-

DPOSRES将导致定义关键字"POSRES"。查看拓扑文件的最后一行，看看怎样进行位置

限制。采用grompp和mdrun程序运行模拟：

Remark:changethevalueofgen_ed

-o

mdrun-v-deffnmprotein-PR

Whatisthelengthofthesimulationinpicoconds?

Howistheinclusionofthepositionrestraintshandledinthetopologyfile?

你会很有趣地看到总能量、势能和动能的变化。在前一步，粒子没有速度，所以没有

动能，因此也没有温度。现在，模拟一开始，原子就被赋予速度因此获得动能。模拟过程

的能量信息被保存为不可读的二进制文件格式（.edr文件）。该文件的信息可用gromacs

工具g_energy程序提取。

注：在下一步模拟运行时，你能用以下命令并在执行过程中回答有关问题，以节省时

间。如果这一步不奏效，让该程序继续在这里运行，而打开另外一个终端，运行下一步的

grompp和mdrun程序，但要注意两个终端的都要在正确的文件夹下运行。

g_

你将要在一列能量参数中进行选择，选择好的参数结果将会被输出。输入与温度、势

能、动能和总能量相对应的数字，然后输入0（零），回车。

输出文件(.xvg)可以用xmgr或xmgrace程序查看，这些程序需要预先安装好。也可

以用Python语言脚本查看，它能在终端窗口上显示出图线。

注：绿色曲线是总能量，黑色是势能，蓝色是温度，红色是动能。你也可以点击曲线

并填写“图标”窗口以改变图标。别忘了点击'accept'以接受改变。

Whathappenswiththetemperature?

Whathappenstothepotential/kinetic/totalenergyandhowcanthisbeexplained?

九、非限制性MD，打开温度耦合

系统现在应该相对松弛了一些，我们就引入温度并开始对系统进行热浴耦合吧。换句

话说，将在打开温度耦合的情况下运行一小段时间，让系统松弛到一个新的状态。下载控

制参数文件，然后检查一下其中的温度耦合参数；最后用grompp和mdrun运

行模拟：

roups

arethatandwhatdoyouthinkisineachofthegroups?

运行结束后，查看一下能量和温度。再次提醒，你可以在查看能量的同时，运行下一

步。数据提取方法同前：

观察能量图。

Whathappenstothetemperature?

Whathappenstothepotential/kinetic/totalenergyandhowcanthisbeexplained?

十、非限制性MD模拟：打开压力耦合

当引入温度并松弛后，就该引入压力了。导入控制参数文件，短时间运行有

压力下的模拟。先查看一下此控制参数文件，找到控制引入压力的句柄。此处也应该注意，

因为我们需要重新赋予离子速度，通过gen_vel行实现。在这行下面，设定了速度分布

（Maxwell分布）的温度，并通过gen_ed产生初始值。gen_ed现在设定为9999，

但是会被改为一个新的值。

MD模拟是最重要的一步，所以每次都以相同的坐标体系、速度和相同的参数开始模

拟时，就会导致模拟的结果一模一样，这不是我们所希望的。

当你设定好了控制参数后，运行以下命令进行模拟。

Remark:changethevalueofgen_ed

-o

mdrun-v-deffnmprotein-NPT

运行完成后，再次查看能量和温度，同时注意观察压力变化。数据提取方法同前：

g_

看看能量曲线、温度和压力的变化。

Whathappenstothetemperature?

Whathappenstothepotential/kinetic/totalenergyandhowcanthisbeexplained?

最后一个问题与从有限数量的粒子系统中提取热动力特性直接有关。粗略地说，热动

力学指的是大量粒子（例如，几十亿个而不是几千个）的行为。取大量粒子的平均特性能

够使数据波动减少，相反，仅对少量粒子进行平均，波动就较大。

由于我们引入了压力，系统的密度会发生变化，所以从能量文件中提取密度数据，方

法如下：

Howdoesthedensitybehaveovertime?

Whydoesthedensityofthesystemchangeifpressurecouplingisturnedon?

十一、正式模拟

到了最后一步了。我们已经有了好赖平衡了的含溶剂系统，系统中包含我们感兴趣的

蛋白质，所以该进行正式模拟了。记住，正式模拟并不表示整个模拟都可以用来分析感兴

趣的特性。虽然一些初始参数值被擦除了，系统看起来还不太可能就已经达到平衡状态了。

在分析阶段，我们应当检查哪一部分模拟结果可以认为是在平衡状态下得到的，只有这部

分数据才适合用于分析。先设定运行参数。这里只需要运行模拟步骤，类似于系统准备的

第二步。然而，这一步也是需要考虑模拟目的的，因此控制参数里面应当包含能反映分析

内容的有关参数。考虑以下问题：

•在什么时间尺度上能够反映所研究的问题？需要运行多久？

•需要得到多少帧数据，时间解析度是多少？

•需不需要存储速度数据？

•是否需要输出所有的原子数据，还是只要有蛋白质坐标数据就够了？

•每隔多久记录一次能量文件和日志文件？

•每隔多久记录一次坐标和速度文件？

•...

我们将运行10纳秒的模拟。下载或拷贝参数文件，先检查一下文件内容。

下载或拷贝文件并查看其中的内容。

Howmanystepsareneededtogetatotalsimulationtimeof10ns?

你要算一下，为了完成10纳秒的模拟，需要多少步数，并编辑好参数文件；然后用

Grompp命令将最终的结构文件和系统准备步骤中得到的拓扑文件合并为一个输入文件。

虽然输入文件是二进制格式，我们还是可以看一看其内容。某些情况下，会发生一些

意外，情况，其中有些意外情况可能与内部控制和力场参数有关，在这些情况下，这么做

就尤为有用。利用gmxdump程序，可以读取这个文件。读出的文件可能会有很多页，建

议这个时候用'more'或'less'命令分屏显示文件。输入以下命令读取输入文件（注意，在

实际操作时，需要把文件名换为运行grompp程序所得到的文件名）：