Gromacs教程-水中的溶菌酶

admin管理员组

文章数量:1568433

2024年3月30日发(作者：)

GROMACS 教程

水中的溶菌酶

Justin Lemkul

Department of Biochemistry, Virginia Tech

YongMa2008@小木虫 译

写在前面:

1 本人没有系统学习过MD, 本教程是本人开始自学GROMACS的

入门教程，在课余时间翻译的，对于很多专有词汇，不是很懂，翻译

的不尽正确，欢迎批评指正

2 由于版本问题，本教程中某些命令行或需要输入的变量可能会有不

同，请自行斟酌

端显示文字或文件内部文字

/Pages/Personal/justin/gmx‐tutorials/lysozyme/

y

o

5 原版教程链接如下，强烈建议有兴趣的童鞋学习原版

n

g

m

a

同时由于时间仓促难免出现排版错误，请见谅

2

0

4 由于Word排版问题，某些命令行中的空格不是很明显，请注意；，

0

8

@

小

木

3 本教程中红色为命令行，蓝色为超级链接，橙黄色为程序运行时终

虫

本例将指导新用户进行设置一个蛋白质(lysozyme)加上离子在水盒子里的

模拟过程。每个步骤将包含对输入输出文件一般应用的典型设置的解释。

这个教程假定你正在用gromacs的4.5.x版本。

第一步，准备拓扑

我们必须先下载我们要用的蛋白质结构。在这个教程中，我们将用鸡蛋清溶

菌酶(PDB代码1AKI). 去RCSB

(/pdb/home/) 网页

下载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. 后处理文件

拓扑文件(默认为)包含定义分子模拟的所有信息。这些信息包含

非成键参数(原子类型和电荷)和成键信息(键长，键角和二面角)。拓扑文件形成

后，我们要仔细看一下。

用如下命令执行pdb2gmx命令：

pdb2gmx -f -o 1AKI_ -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)

y

o

n

g

m

a

2

0

0

8

@

小

木

虫

力场包含将要写入拓扑文件的信息。这是个非常重要的选择。你必须完全了

解每个力场并且决定哪个力场对你的研究对象是最适合的。这个教程中，我们选

用全原子OPLS力场，在命令提示行输入14，然后按Enter键。

-ignh: 忽略PDB文件中的氢原子，特别是对NMR结构有用。否则如果氢

原子存在，他们必须有正确的顺序和名称，正如Gromacs所期望的那样。

•

-ter：为N-和C-末端交互式分配电荷态

•

-inter：为谷氨酸，天门冬氨酸，赖氨酸，精氨酸，和His交互式指定电

荷态；分配二硫键到胱氨酸

•

现在你已经得到了三个文件，1AKI_, , 和

. 1AKI_是个gromacs格式结构文件，包含立场中定义

的所有原子(氢原子已经被加到蛋白质中的氨基酸上了)。 文件是系

统拓扑文件(稍后会解释)。 文件包含用来限制位置的信息(后面解

释)。

第二步，检查拓扑文件

我们来看一下拓扑输出文件()中有什么。再次使用纯文本编辑器，

检查其内容。在几句注释行之后(前面用分号；标注)，你会找到下面的语句：

#include "/"

y

o

n

g

m

a

2

0

还有些pdb2gmx会用的其他选项，下面列了几个：

0

8

@

小

木

虫

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- 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: 无需解释

y

o

n

g

m

a

2

0

0

8

@

小

木

虫

Qtot 描述保持对分子总电荷的持续计数

mass: 也无需解释

•

typeB,chargeB,massB：用于自由能微扰(这里不讨论)

•

下面的项，包括 [ bonds ], [ pairs ], [ angles ], and [ dihedrals ].

这些项都是无需解释的。这些项中的参数和函数类型在gromacs手册的第5章中

阐述。特别的，1-4相互作用包含在“pairs“(gromacs手册中5.3.4部分)

从位置限制拓扑文件开始，提示部分涉及了一些有用的/必须的拓扑文件。

""是从pdb2gmx命令产生的，它定义了一个力常数用于做平衡时能保

持原子不动(后面会很多涉及)。

; Include Position restraint file

#ifdef POSRES

#include ""

#endif

这里结束了蛋白质A分子类型的定义。拓扑文件的提示部分用来定义其他分

子并提供系统层次的描述。下面的分子类型(默认)是溶液，在这里是SPC/E水。

其他的对水的选择包括SPC, TIP3P, 和TIP4P。我们通过在pdb2gmx命令中使用

"-water spce"选这个水模型。对于其他不同水模型的详细描述，可以点这里

(/water/) ，但是注意gromacs没有包含所

有的模型。

; Include water topology

#include "/"

#ifdef POSRES_WATER

; Position restraint for each water oxygen

[ position_restraints ]

; i funct fcx fcy fcz

1 1 1000 1000 1000

#endif

正如你看到的，水分子也可以被位置限制，通过利用一个力常数(

k

pr

) 为

1000 kJ mol

-1

nm

-2

. 离子的信息也别包含在内：

; Include generic topology for ions

#include "/"

最后是系统层次定义：

[ system ]

; Name

y

o

n

g

m

a

2

0

0

8

@

小

木

虫

LYSOZYME

[ molecules ]

; Compound #mols

Protein_A 1

[ system ] 直接给出模拟中将要写到输出文件中的系统名

字。 [ molecules ]直接列出系统中的所有分子(直接从[ moleculetype ]的名

字中列出)，为了他们在坐标文件(.gro)中出现。

现在我们检查完了拓扑文件中的内容，我们可以继续组建我们的系统了。

第三步，定义单位胞 & 添加溶剂

现在你对gromacs的拓扑文件比较熟悉了，到了继续创建系统的时候了。这

个例子中，我们将要模拟一个非常简单的水系统。模拟蛋白质或者其他分子在不

同溶液中也是可以的，只要涉及的分子的参数是合适的就行。

定义盒子和添加溶液有两步：

1. 用editconf定义盒子维度

2. 用genbox往盒子里填水

用editconf定义箱子：

editconf -f 1AKI_ -o 1AKI_ -c -d 1.0 -bt cubic

上面的命令把蛋白质放在盒子中心(-c), 并且把它放在距离盒子边缘1.0

纳米位置(-d 1.0). 盒子类型是立方体(-bt cubic)。 盒子边缘的距离是个重要

参数。 因为我们要用到边界条件，我们必须满足最低限度的图像公约，那就是

一个蛋白质永远不能看到他的周期图像，否则计算的力是伪力。指定一个距离是

1.0 纳米的溶剂盒子意味着两个周期性蛋白质图像的距离至少是2.0纳米。这个

距离对于模拟中经常用到的截止模型都是足够的。

定义了盒子之后，我们要填充溶质(水)了。用genbox添加溶质：

genbox -cp 1AKI_ -cs -o 1AKI_ -p

蛋白质的配置(-cp)包含在前面一步editconf的输出文件中，溶剂的配置

(-cs)是标准gromacs安装的一部分。我们要用，它是个通用的3点

平衡溶剂模型。你可以用为SPC, SPC/E, or TIP3P水模型。使用

y

o

n

g

您现在面临一个如何处理的晶胞选择。对于本教程的目的，我们将使用一个

简单的立方体晶胞。当你对周期性的边界条件和箱型更了解时，我强烈推荐菱形

十二面体，因为在同一周期距离下，其体积为立方体的71％，因此可以减少需

要加入的熔剂水分子的数量。

m

a

2

0

0

8

@

小

木

虫

作为溶液配置，因为他们都是三点型的水模型。输出文件是1AKI_ 我

们告诉genbox拓扑文件的名字()，所以它能修改。注意

中[ molecules ] 的变化：

[ molecules ]

; Compound #mols

Protein_A 1

SOL 10832

genbox做的是追踪了增加了多少水分子，这个数字会写入你的拓扑文件显

示它已经做了。注意如果你用其他的(非水)溶剂，genbox将不会在拓扑文件中

写入信息。它与水分子的兼容性更新是硬编码

第四步，添加离子

现在我们有了包含一个带电荷的蛋白质的溶液盒子。pdb2gmx的输出文件告

诉我们这个蛋白质包含+8e的静电荷(由于它的氨基酸组成)。如果你忽略了

pdb2gmx输出文件的这个信息，可以看看你的文件中的最后一行，它

将显示“qtot 8.”. 因为生命中不存在静电荷，所以我们需要往我们的系统里

面添加离子。

gromacs程序中添加离子的命令是做的是从拓扑文件中读取

数据并用你选用的离子去替换水分子。这个输入文件叫做运行输入文件，拓展名

是.tpr. 这个文件是用gromacs的grompp(Gromacs Pre-Processor) 工具产生

的，它也将在我们随后运行第一次模拟时用到。grompp做的是处理坐标文件和

拓扑文件(定义分子)来产生一个原子层次的输入文件(.tpr)。.tpr文件 包含了

系统中所有原子的参数。

为了用grompp产生.tpr文件我们还需要另外一个拓展名为.mdp

(molecular dynamics

parameter file)的输入文件，grompp 将组合mdp文件中

设定的参数和结构和拓扑文件来生成.tpr文件。

.mdp文件在能量最小化或者分子动力学模拟中常用的文件，但在这里只是简单

的用来生成一个系统中原子的描述。一个.mdp例子(我们将用)可以在这里

下载。

(/Pages/Personal/justin/gmx-tutoria

ls/lysozyme/Files/)

实际上，这里用的.mdp文件可以包含任何参数的合理组合。我通常用一个

能量最小化脚本，因为它非常基础而且不用涉及任何复杂的参数组合。

下面的命令将组合你的.tpr文件。

grompp -f -c 1AKI_ -p -o

y

o

n

g

m

a

2

0

0

8

@

小

木

虫

现在我们有一个在原子层次描述我们系统的.tpr文件。这个文件将用在

genion工具中：

genion -s -o 1AKI_solv_ -p -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 ]中一致。残

基名称可能带或者不带电荷的符号(+/-). 如果你遇到什么困难你可以参考

中的正确名字。

你的 [ molecules ] 指令现在看起来就像这样：

[ molecules ]

; Compound #mols

Protein_A 1

SOL 10824

NA 0

CL 8

这里的“NA 0”表示系统里没有钠离子，你可以安全的删除这行或者留着。

前期gromacs版本会报错。Gromacs-4.0和以后版本允许0分子存在。

y

o

n

g

m

a

2

0

0

8

@

小

木

虫

第五步，能量最小化

加了溶剂的电中性的系统现在已经做好了。我们开始做动力学之前，我们必

须保证我们的系统里面没有立体几何冲突或不当。 结构要通过一个叫能量最小

化的过程弛豫一下。

能量最小化过程就跟添加离子差不多。我们再一次使用grompp来把结构，

拓扑和模拟参数组合到一个二进制输入文件里面(.tpr), 但是这次我们要通过

gromacs的引擎mdrun，而不是genion来进行能量最小化。

通过grompp用这个

参数的文件，组合二进制输入文件。

/Pages/Personal/justin/gmx-tutoria

ls/lysozyme/Files/

grompp -f -c 1AKI_solv_ -p -o

确保你在运行genbox和genion时已经更新了你的文件，否则你

将得到一堆错误信息(坐标文件中的坐标跟拓扑文件不匹配，等等).

mdrun -v -deffnm em

: ASCII-文本的能量最小化过程日志文件

•

: 二进制能量文件

•

: 二进制高精度轨迹

•

: 能量最小化结构

•

有两个重要的因素来评价能量最小化是否成功。第一个是势能(在能量最小

化过程的最后面输出，即使你没用-v标志)。Eopt应该是负值，根据系统大小和

56

水分子的多少，大约在10

-10数量级(对水中的蛋白质而言)。第二个重要的指

标是力的最大值。Fmax，在中设置的目标是“emtol=1000.0“表示

有可能当Fmax>emtol时得到一个合理的Eopt。

Fmax不能大于1000 kJ mol

-1

nm

-1

。

如果出现这样的问题，你的系统对于模拟来讲可能不够稳定。评价一下为什么会

这样，可能要改一下你的关于能量最小化的参数设置(integrator, emstep 等

等)

我们来做一些分析。文件中包含gromacs在能量最小化过程中得到

的所有的能量。你可以用gromacs里面的g_energy来分析任何一个.edr文件：

y

o

n

g

使用- v标志因为我们没什么耐心：它使mdrun冗长，使它在屏幕上打印每

一步的进展情况。使用- deffnm标志将定义输入和输出的文件名。所以如果你

没有名字，你的grompp就输出“”，你必须明确指定其与mdrun -s标志

的名称。就我们而言，我们将得到以下文件：

m

a

2

0

0

8

@

小

木

现在我们做好了调用mdrun来运行能量最小化的准备：

虫

g_energy -f -o

在提示行中，输入”10 0”来选择potential(10),zero(0)来结束输入。

Eopt的平均值将被显示，写入一个叫做的文件。想画出这些数

据，你需要用Xmgrace描点工具。结果中的点看起来就像这样，展示一个漂亮的

稳定的Eopt收敛曲线。

EM保证我们在结构和溶液取向方面有一个合理的初始结构。为了开始真正

的动力学模拟，我们必须对蛋白质周围的溶剂和离子做一些平衡。如果我们要在

这个点来尝试非限制性动力学，系统将会崩溃。原因是溶剂只在自己内部优化还

没考虑溶质。我们想要模拟的温度也要引入来确定关于溶质(蛋白质)的合适的取

向。达到正确的温度(取决于动能)之后，我们要对系统施加压力知道它达到合适

的密度。

还记得我们好久前用pdb2gmx生成的文件么，现在要派上用场了。

文件的目的就是对蛋白质中的重原子(非氢原子)实施位置限制力。除

非经过大量的能量惩罚，否则是不能动的。位置限定的用途就是在没有引起蛋白

质结构变化的前提下，允许我们来平衡蛋白质周围的溶剂。

平衡往往是分两个阶段进行。第一阶段在一个NVT系统（固定数目的粒子，

体积和温度）下进行。这个系统也被称为“等温等容”或“规范”。这种程序的

时间表与该系统的内容相关，但在NVT系，系统的温度应达到预期值的平台。如

果温度还没有稳定下来，将需要更多的时间。通常情况下，50至100 ps的就足

够了，我们会为这项工作进行了100 ps的NVT系平衡。根据您的机器，这可能

需要(对双核的MacBook刚刚超过一小时)一段时间。从这里取得.mdp文件

y

o

n

g

m

a

第六步，平衡

2

0

0

8

现在我们的系统已经在能量最低点，我们可以开始真正的动力学模拟了。

@

小

木

虫

/Pages/Personal/justin/gmx-tutorial

s/lysozyme/Files/

我们将要像在能量最小化过程中一样来调用grompp和mdrun

grompp -f -c -p -o

mdrun -deffnm nvt

对采用的参数的全部解释可以从gromacs手册找到，包括一些解释。注意一

些.mdp文件中的参数：

gen_vel = yes: 产生初始速度。利用不同的随机种子(gen_seed)，提供

不同的初始速度,因此多个(不同的)模拟可以从一个相同的初始结构开

始.

•

tcoupl = V-rescale: 速度重新调整恒温器是对Berendsen弱耦合方法的

改进，该方法不能产生正确的动能。

•

pcoupl = no: 没有用压力耦合

•

g_energy -f

从这些点中可以看出，系统温度很快就达到了目标温度(300K)，对其余的平

衡保持稳定。对于这个系统，较短的平衡时间（50 ps）已经足够。

第七步，平衡，第二部分

y

o

n

g

m

a

2

0

在提示行输入“15 0”来选择系统温度并退出。结果显示：

0

8

@

小

木

再次用g_energy来分析温度进展：

虫

先前那步，NVT平衡，稳定了系统的温度。采集数据之前，我们还要稳定系

统的压强（当然还有密度）。压强平衡是在NPT条件下进行集成，其中,粒子，压

强，和温度都不变。该条件也被称为“等温等压”，与实验条件最相近。

100ps的NPT平衡的.mdp文件可以在这里

(/Pages/Personal/justin/gmx-tutoria

ls/lysozyme/Files/) 找到。该文件与NVT平衡时用的mdp文件没有太

大的不同。注意添加的压强耦合项，使用Parrinello-rahman恒压器。

几项其他的变动：

continuation = yes: 我们将从NVT平衡状态继续模拟

•

gen_vel =速度从轨迹中读取

•

我们将调用在NVT平衡时用过的grompp和mdrun。 注意我们现在要用到-t

来包含NVT平衡过程中的检查点文件。这个文件包含了继续模拟所需要的变量。

我们需要包含这个文件来保存在NVT过程中得到的速度。坐标文件(-c) 是NVT

模拟的最后文件。

g_energy -f -o

y

o

在提示行输入“15 0”来选择系统压强并退出。结果如下图所示：

n

g

m

a

2

0

我们来分析压强进展，再次用到g_energy.

0

8

@

小

木

grompp -f -c -t -p -o

mdrun -deffnm npt

虫

在100ps的平衡阶段，压强值有较大的波动，这并不意外。数据的平均值在

图中用红线标示。在该平衡过程中，压强平均值是1.05bar.

我们再来看看密度，这次在g_energy的提示行输入“21 0”

g_energy -f -o

第八步，进行MD

随着两个平衡阶段的完成，系统现在已经在需要的温度和压强下平衡良好。

我们已经准备好了进行位置限制和并为取得数据而进行MD了。这个过程跟我们

以前做的类似。我们还要用到检查点文件(该情况下，它包含了压强耦合信息) 来

进行grompp。 我们要进行一个1ns的MD模拟，脚本文件可以从这里

(/Pages/Personal/justin/gmx-tutoria

ls/lysozyme/Files/) 得到。

grompp -f -c -t -p -o md_0_

mdrun -deffnm md_0_1

Mdrun步骤最好在并行计算机计算，因为它要花上几小时才能完成。执行这

个任务合适的命令是：

mpirun -np X mdrun_mpi -deffnm md_0_1

y

o

n

g

跟压强一样，运行得到的密度平均值也用红线标示。100ps过程中的密度平

均值为998.3 kg m

-3

, 与实验值1000 kg m

-3

接近。SPC/E水模型的参数比较接近

的克隆了水的实验值。该密度值比较稳定，意味着系统在这个压强和密度下平衡

良好。

m

a

2

0

0

8

@

小

木

虫

其中X表示模拟所需要用的处理器数目。

第九步，分析

现在我们已经模拟了我们的蛋白质，我们应该来分析一下我们的系统。什么

类型的数据才是重要的呢?这是模拟前问的一个重要问，所以你应该对你需要采

集的系统的数据类型有自己的想法。在本教程中，我们介绍了一些基本工具。

第一个是trjconv，这是一个后期处理工具，用来处理坐标，纠正周期性或

者手动来调整轨迹(时间单位，帧频率等)。对于本练习，我们要用trjconv来核

算系统中的周期性。蛋白质可能会在单元格中扩散，也可能“跳”在两个盒子之

间。我们采用下面命令来处理这种情况。

trjconv -s md_0_ -f md_0_ -o md_0_1_ -pbc mol -ur

compact

选择0("System")。我们要对这个“修正”了的轨迹进行我们的分析。先来

看看结构稳定性。GROMACS有个内置工具进行RMSD计算，名字叫g_rms。用这个

命令来使用这个工具：

g_rms -s md_0_ -f md_0_1_ -o -tu ns

对最小二乘拟合和群计算RMSD，都选择4("Backbone")。-tu 标志将输出结

果以ns显示，即使轨迹文件是用ps写的。这是为了输出文件更加清晰(尤其是

你做一个较长的10

5

ps的模拟时)。输出文件将显示在最小化的，平衡了的系统

的结构下的RMSD.

y

o

n

g

m

a

2

0

0

8

@

小

木

虫

如果我们要计算晶体结构的RMSD值，我们需要用到下面命令：

g_rms -s -f md_0_1_ -o rmsd_ -tu ns

结果将如下所示：

g_gyrate -s md_0_ -f md_0_1_ -o

y

o

n

g

蛋白质的回转半径可衡量其紧凑度。如果蛋白质稳定的折叠了，Rg将保持

一个相对温度值。如果蛋白质展开了，它的Rg将随时间变化。我们来分析一下

我们模拟的溶菌酶的回转半径。

m

a

2

0

0

8

两个数据显示RMSD大约在0.1nm，表示该结构非常稳定。两者的微小不同

在于当t=0ns时晶体结构是不同的。这是可以预料的，因为它已经被能量最小化

了，而且如我们前面讨论的，位置限定不是100%完美的。

@

小

木

虫

从这个合理的恒定的Rg值我们可以看到，在温度为300k的1纳秒时间内蛋

白质的紧凑型形式（折叠）非常稳定。这一结果并不出人意料，而这足以说明

GROMACS较高的分析能力。

总结

现在已经用gromacs操作了一个分子动力学模拟过程，并分析了一些结果，

本教程不是个全面的东西。用gromacs你可以有很多类型的模拟可以操作(自由

能计算，分平衡MD,和正常模式分析还有一些其他的)。你应该看一些文献和

gromacs手册来调整这里提供的.mdp文件来达到更有效地和更精确的目标。

如果你对改善这个教程有什么建议，如果你发现了错误，或者什么东西不清

楚，请不要客气给我发邮件

(jalemkul@)。注意：这不是邀请你为了任何

gromacs的问题而给我发邮件。我没有为了一个私人的教程或者个人帮助服务而

为自己打广告。

享受你的模拟吧！

y

o

n

g

m

a

2

0

0

8

@

小

木

虫

本文标签： 文件系统蛋白质拓扑模拟