gromacs_tutor_翻译_

GROMACS Introductory Tutorial

Gromacs ver 4.0

Author: John E. Kerrigan, Ph.D.

Associate Director, Bioinformatics

The Cancer Institute of NJ

195 Little Albany Street

New Brunswick, NJ 08903

Phone: (732) 235-4473

Fax: (732) 235-6267

Email: kerrigje@https://www.360docs.net/doc/207262485.html,

翻译：梁（leunglm@https://www.360docs.net/doc/207262485.html,）

1

GROMACS教程：蜘蛛毒素肽的研究

Yu, H., Rosen, M. K., Saccomano, N. A., Phillips, D., Volkmann, R. A., Schreiber, S. L.: Sequential assignment and structure determination of spider toxin omega-Aga-IVB. Biochemistry 32 pp. 13123 (1993)

GROMACS是一个使用经典分子动力学理论研究蛋白质动力学的高端的高效的工具[1]。这个软件包是遵守GNU许可的免费软件，可以从以下站点下载：https://www.360docs.net/doc/207262485.html,。GROMACS可以在linux，unix，和Windows（新开发的）上使用

摘要：在本教程中，你将研究一个从漏斗形蜘蛛的毒液中分离的毒素。过去毒液毒素用来鉴定阳离子通道。钙离子通道调节这种离子进入细胞。神经信号受到神经细胞中离子平衡的高度控制。人们认为象蜘蛛毒素这类毒液中暴露的带正电的残基会倾向于与细胞离子通道入口的带负电的残基结合。本教程中的蜘蛛毒素中带正电的残基主要朝向肽链的一侧。离子通道的堵塞导致了神经信号的中断，最终导致麻痹和死亡（通过呼吸判断）。

我们将使用显性溶剂动力学的方法研究这个小肽。首先比较真空中和溶解的模型。我们将把小肽溶在水盒子里，紧接着用牛顿运动定律加以平衡。我们还将比较和对比补偿离子在显性溶剂动力学中的影响。我们希望解决如下问题：

小肽的二级结构在动力学条件下是否稳定？

带正电的残基侧链是否主要朝向肽结构的一侧？补偿离子是保持在正电残基附近还是在

四周移动？

水在维持蛋白结构中扮演什么角色？

下载pdb文件1OMB.PDB(https://www.360docs.net/doc/207262485.html,/pdb/)

建议先用Deep View软件观察下载的文件，看是否有混乱的现象（缺失侧链）Deep View 可以补充缺失的侧链（但是，注意Deep View可能会在添加的侧链前加上特殊符

号，而且这些符号只能手工去掉！）这个pdb文件不存在侧链的缺失，我们就不必担心

了。

在你的unix主目录下建立一个名为“fwspider”的文件夹。在此文件夹下新建三个子

文件夹，分别命名为：“invacuo”“wet”和“ionwet”。用sftp将1OMB.pdb拷贝到

2

fwspider的每个子文件夹中。

（重要！无论何时将一个文本文件从windows系统拷贝到unix系统，一定要转成unix

文本文件？？可以用to_unix命令（如to_unix filename filename，将filename文件转成unix

文本文件。在RedHat Linux中，用dos2unix命令。）象MS Word这类Windows文本编辑器

加入的控制符可能在unix程序中产生错误。）

用pdb2gmx处理pdb文件。pdb2gmx（查看选项用pdb2gmx –h; 其实可以用-h选项查看所有Gromacs命令的帮助文档）命令将pdb文件转换成gromacs文件并产生拓扑文件。因为本pdb 文件是由NMR产生的，含有氢原子，因此用-ignh选项忽略文件中的氢原子。-ff选项用来指定力场（G43a1 是Gromos96力场，一个通用原子力场）。-f读入pdb文件，-o指定一个新产生的pdb文件（也可以是其它多种类型文件）的文件名。-p指定新产生的拓扑文件名。拓扑文件包含了所有力场参数（基于一开始选择的力场），因此非常重要。研究表明SPC/E 水模型[2]在水盒子模拟中表现最好。用spce水模型研究长程静电相互作用较好。因此我们用-water来指定水模型。

pdb2gmx –ignh –ff G43a1 –f 1OMB.pdb –o fws.pdb –p fws.top –water spce

建立盒子。

editconf -bt cubic –f fws.pdb –o fws.pdb –d 0.9

用上面的命令建立了一个简单的立方体盒子。-d决定了盒子的尺寸，即盒子边缘距离分子边缘0.9nm（9?）。理论上在绝大多数系统中，-d都不能小于0.85nm。[4]

[注释：editconf也可以用来进行gromacs文件（*.gro）和pdb文件（*.pdb）的相互转化。例如：editconf –f file.gro –o file.pdb 将file.gro 转换为 file.pdb]

现在就可以用产生的文件进行真空模拟了。真空模拟就是先能量最小化，然后进行动态模拟。（思考真空模拟为何不需要位置限制动态模拟？）

在盒子中放入溶剂

genbox –cp fws.pdb –cs spc216.gro –o fws_b4em.pdb –p fws.top

genbox命令在editconf产生的盒子基础上生成水盒子。上面的命令行指定了spc水盒子。genbox命令可以在给定尺寸的盒子中加入正确数目的水分子。

设置能量最小化

用em.mdp文件。Gromacs用*.mdp文件指定所有计算的参数。看一下文件的内容。它确定

3

用最陡下降法消除范德华位置冲突。编辑文件，将nsteps变成400。如果最小化不能收敛，就用nsteps=500再做一次。（最小化在400步内是能收敛的，但不同的平台可能结果会不一样。）要重做的话，必须重新运行grompp（注意：预处理器的位置在你的机器上可能不同，用which命令来定位[即 which cpp]！）

em.mdp内容：

em.mdp文件用的重要概念：

title –标题随便取（最长64个字，简单点好）

cpp –指定预处理器的位置

define –传递给预处理器的一些定义。–DFLEXIBLE 告诉grompp将spc水模型而非刚性SPC包含进拓扑文件，以便用最陡下降法进一步最小化能量。

constraints-设置模型约束

integrator – steep告诉gompp运行最陡下降法进行能量最小化。用cg进行共轭梯度法。

dt –能量最小化用不用。只在动力学模拟中用（如md）。

nsteps – 在能量最小化中，指定最大反应步数。

nstlist –更新邻居列表的频率。nstlist = 10表示每10步更新一次。

rlist – 短程邻居列表的阈值。

4

5

coulombtype –告诉gromacs 如何计算静电。PME 为particle mesh ewald 法（参见Gromacs 用户手册）

rcoulomb –指定库仑力阈值

vdwtype – 告诉Gromacs 如何计算范德华作用(cut-off, Shift 等）

rvdw –指定LJ 或Buckingham 势能距离阈值

EM Stuff emtol –最大的力如果小于此值则能量最小化收敛（结束）（单位是kJ mol

–1 nm –1） emstep –初始步长（nm ）

下面用grompp 程序进行文件处理。grompp 是预处理程序（即the gromacs pre-processor 的缩写）

grompp –f em.mdp –c fws_b4em.pdb –p fws.top –o fws_em.tpr

grompp 中-f 标签制定输入参数文件（*.mdp ）。-c 输入结构文件（pdb 文件，*.pdb ）; -p 输入拓扑文件，-o 输出mdrun 的输入文件（*.tpr ）。

加氯离子后的pr_md.mdp 和md.mdp 文件的温度耦合参数

; Berendsen temperature coupling using vel rescaling is on

Tcoupl = v-rescale

tau_t = 0.1 0.1

tc_grps = protein non-protein

ref_t = 300 300

记住：如果要加入氯离子，需要重新运行grompp。首先删除旧的fws_em.tpr文件，然后运行下面的grompp命令。我们加入氯离子是为了中和模型中的总静电荷。

grompp –f em.mdp –c fws_ion.pdb –p fws.top –o fws_em.tpr

在后台运行能量最小化（在命令后加&）

nohup mdrun –v –s fws_em.tpr –o fws_em.trr –c fws_b4pr.pdb –e em.edr –g em.log &

用tail命令检查最小化的进程

tail –15 em.log

当能量最小化结束，你将看到log文件中有如下总结文字，表明最陡下降收敛了。

用tail -50 em.log

设置位置限制性模拟。什么是位置限制性模拟？你限制（或部分冻结）大分子中的原子位置，而允许溶剂分子运动。这样做像是将水分子浸入大分子。水分子松弛时间约为10ps。因此我们要进行超过10ps的位置限制性模拟。本实例中用20ps（最少高出一个数量级？？）大的模型（大蛋白或脂）可能需要更长的平衡时间，50ps或100ps或更长。下面的设置在这个gromacs力场中运行良好。其他力场请参考用户手册（例如在GROMOS96力场中，建议nstlist = 10 and rvdw = 1.4）。在coulombtype, PME代表“Particle Mesh Ewald”静电势。[5，6] PME是计算长程静电势的最优算法（给出最可信的能量评估，尤其在用Na+, Cl-, Ca 2+等作为补偿离子的体系。）由于这个蛋白具有暴露的带电残基，使系统带有+2静电荷，所以适用PME算法，更为有益的是用补偿离子使系统处于电中性。constraints中的all-bonds选项可以应用线形限制算法确定系统中的所有键长（当dt > 0.001 ps时尤为重要）[7]。学习一下下面的mdp文件。

pr.mdp:

6

7

mdp 文件中的重要参数

define 声明中的–DPOSRE 告诉Gromacs 运行位置限制动力学模拟

constraints 声明如前所述。all-bonds 设定LINCS 算法限制所有键。[7]

Integrator 告诉gromacs 进行何种动态算法（另外的选项“sd ”代表stochastic dynamics ） dt 是每步的时间（我们选择了2fs ；但此处的单位一定是ps ！）

nsteps 是运行的步数（总模拟时间= nsteps * dt ）。

nstxout 告诉gromacs轨迹文件收集模拟快照（坐标）的频率（nstxout = 250 且dt = 0.002，所以每0.5ps收集一张快照）

coulombtype选择gromacs计算原子静电相互作用方法（PME 代表 particle mesh ewald；另

外还可以用cut-off）。

rcoulomb 和rvdw是计算静电和范德华作用的阈值（单位nm，1.0nm=10.0埃）

温度耦合部分非常重要，必须正确填写。

Tcoupl = v-rescale [8, 9]（用随机条件重新调解速度的温度耦合类型。）

tau_t =温度耦合的时间常数（单位ps）。必须每个tc_grps指定一个，且顺序对应。

tc_grps =与调温器耦合的组（模型中的每个原子或残基都用一定的索引组表示）

ref_t =代表耦合的参照温度（即动力学模拟的温度，单位K）。每个tc_grp对应一个ref_t.

预处理pr.mdp文件

grompp –f pr.mdp –c fws_b4pr.pdb –p fws.top –o pr.tpr –maxwarn 3

nohup mdrun –s pr.tpr –o pr.trr –c fws_b4md.pdb –e pr.edr –g pr.log &

用tail命令检查pr.log文件。

md.mdp文件和pr.mdp文件相仿。有几处不同。define声明不再需要，因我们不再做位置限

制模拟。

用于explicit solvation 的md.mdp文件内容（特殊注释：做真空模拟，去掉温度耦合中的“sol”部分。在有补偿离子的模拟中，为离子加入相应的温度耦合参数）

8

9

grompp –f md.mdp –c fws_b4md.pdb –p fws.top –o md.tpr –maxwarn 3

nohup mdrun –s md.tpr –o md.trr –c fws_md.pdb –e md.edr –g md.log &

用tail 命令查看md.log 文件（我们在2.66 GHz Core 2 Quad 机器上用了51分钟跑完100ps 模拟）

可以用trjconv 命令压缩轨迹文件以节省硬盘空间。

trjconv –f md.trr –o md.xtc

得到*.xtc 文件后就可以删除*.trr 文件了。

用ngmx 观察轨迹文件（也可以用VMD 观察轨迹文件）

ngmx –f md.trr (or md.xtc) –s md.tpr

10

当观察器启动后，将看到一个多选项的对话框。选择标“protein ”的多选框，点击OK 。

ngmx 的初始启动对话框

选择“protein”可以只看蛋白分子，而不受盒子中另外约3000水分子的影响。

用X-Rotate上下旋转盒子（鼠标左键向上，右键向下）。用Y-Rotate左右旋转盒子（左键向左，右键向右）。最下面的Scale用来放大或缩小视图（左键放大，右键缩小）。

观察模型中的其它组，点击Display > Filter …，初始对话框就会出现，允许选择观察另外的索引组（如backbone）

要观察模拟轨迹动画，点击Display > Animate。动画播放控制在窗口的底部。点击中间的箭头按钮逐帧观看。点向前的双箭头观看整个轨迹动画，点暂停按钮停止动画。点向左的双箭

头按钮重置动画。

不幸的是，file菜单下的save as pdb选项还不能用。因此，保存并观察*.pdb文件最好的方法是用visual molecular dynamics (VMD)（下载地址：https://www.360docs.net/doc/207262485.html,/Research/vmd/，它是学术免费，且在unix和windows下运行）.

结果分析

Gromacs的一个主要优点（不是指它遵循GNU公共协议而免费！）是有一系列分析轨迹文件的小程序。我们下面就讨论一些比较常用的分析工具。

组

make_ndx

程序make_ndx用来生成组（你想分析的某些特定原子或残基的ID标签）。Gromacs缺省已

经定义了一些组，普通分析可能够用了。但如果你想深入分析，则要用make_ndx程序标注

模型中的特定项。

如何使用make_ndx建立索引文件（ndx）。为了固定某些特定组，或获得一些能量信

息，可以用make_ndx指定这些组。我们来看一个固定蛋白N端和C端的例子。通常用

make_ndx建立索引组供grompp程序调用。

在本例中，我们有一个三螺旋的胶原蛋白结构文件，进行位置限制动力学模拟，我们想固

定N端和C端来进行模拟。首先，确定结构文件（clg_b4md.pdb）的N端和C端残基号。用如

下简单命令。

make_ndx –f clg_b4md.pdb –o clg_ter.ndx

11

12

你将看到如下输出信息（我们省略了开头的一些描述性信息），后面是命令提示符（>）

用“r ”命令输入代表三螺旋N 端和C 端的残基号。

注意：你也可以用连接符指定残基范围（如确定残基1到36，用> r 1-36）

新建索引组的缺省名字(r_1_36_37_72_73_108 )很繁琐，可以用name 命令修改。我们在命令中用索引号#（15）。

13

用“v ”命令查看名字是否改成功了。用“q ”保存并退出。

现在怎么固定组呢？简单，在md.mdp 文件中加入下面几行：

记住当用新的mdp 文件时，首先用grompp 将新索引文件加入tpr 文件。用grompp 的-n 标签，例如：

grompp –f md.mdp –c pr.gro –p clg.top –n clg_ter.ndx –o md.tpr

特性研究

g_confrms

14

要比较最后结构和初始PDB 文件的差异，用g_confrms (用g_confrms –h 查看详细信息)。此程序计算两个结构的最小二乘拟合。

g_confrms –f1 1OMB.pdb –f2 md.gro –o fit.pdb

你将被提示选择一个组（两次都选（组4））。程序将报告RMSD 值，并产生一个输出文件（fit.pdb ）。输出文件中包含两个位置重叠的结构。

g_covar

计算斜方差（参见手册）。也可用于从动态轨迹计算平均结构。如计算1ns 动态模拟的后200ps 的平均结构：

g_covar –f traj.xtc –s topol.tpr –b 801 –e 1000 –av traj_avg.pdb

警告-平均结构往往较粗糙，需进一步执行能量最小化。

g_energy

能量数据作图，如压力、体积、密度等。

g_energy –f md.edr –o fws_pe.xvg

首先要选择输出（*.xvg ）的数据。输出文件是一个电子数据表，可以用Xmgr 或Grace 打开。它是一个文本文件，在进行一些小的改动后可以用Microsoft Excel 打开。

如用上面的命令，你将看到如下结果（你的可能不同）：

15

如计算势能，输入“Potential”，回车

再按一次回车

我们得到一个平均势能和RMSD 的总结（单位kJ/mol ）

输入如下命令用Grace 打开*.xvg 文件：

xmgrace -nxy fws_pe.xvg

可以到以下地址下载Grace http://plasma-gate.weizmann.ac.il/Grace/ .Grace 只能在linux 和unix 上运行。如果没有Grace 或Xmgr ，可以作为空格分隔文件导入MS Excel 。

g_gyrate

用g_gyrate 测量回旋半径。这个指标用于度量结构的紧密度。此程序计算某（些）原子质量与分子重心的关系。

g_gyrate –f md.trr –s md.tpr –o fws_gyrate.xvg

g_rms 与 g_rmsdist

这两个程序用于计算结构的RMSD 值。用g_rms 计算动态模拟过程中的结构与初始结构的结构偏差。（-dt 10 选项告诉程序每10帧计算一次）

g_rms –s *.tpr –f *.xtc –dt 10

计算与NMR 结构的rmsd 值用如下命令：

g_rms –s em.tpr –f md.trr –o fws_rmsd.xvg

选择4组（Backbone ）计算最小二乘拟合。程序生成一个rmsd 随时间变化的图（rmsd.xvg ）。以空格分隔文件的形式导入Excel 。

g_rmsf

计算原子位置的根均方波动（rmsf）。与g_covar相似，此程序也可以计算平均结构。例如，计算一个2ns（2000ps）模拟的后500ps的平均结构，用如下命令：

g_rmsf –f traj.xtc –s topol.tpr –b 1501 –e 2000 –o traj_rmsf.xvg –ox traj_avg.pdb

选择一个RMSD图上平衡的范围（用g_rms计算）。例如：

上面的实例是一个1ns模拟（你的结果可能不同）。这个模拟需要延长到完全平衡。上面

的例子我们用200-500ps范围计算平均结构是因为我们看到这段比较稳定（与原始结构

比）。命令如下：

g_rmsf –s md.tpr –f md.trr –b 200 –e 500 –ox fws_avg.pdb

提示时选择组1“Protein”。

警告-平均结构往往较粗糙，需进一步执行能量最小化。

推荐的真空能量最小化的em.mdp文件。先用最陡下降法，再用共轭梯度法。警告！需要

用pdb2gmx重新生成拓扑文件，尤其是当你选择特定组（而非整体系统）计算平均结构时。

16

17

程序g_rmsf也可以用来计算温度因子。计算的温度因子可以和X光晶体结构的温度因子比较。

g_rmsf –s md.tpr –f md.xtc –o rmsf.xvg –oq fws_bfac.pdb

仍选择“Backbone”组。

do_dssp

用do_dssp命令计算模型的二级结构。前提是你必须在电脑中（/usr/local/bin）安装了

dssp程序(http://swift.cmbi.ru.nl/gv/dssp/ )。

do_dssp –s md.tpr –f md.trr –o fws_ss.xpm

选择计算组1（Protein）。用xpm2ps将xpm文件转成eps格式。然后用ImageMagick转化程序将eps文件转成png文件或其他格式文件。

xpm2ps –f fws_ss.xpm –o fws_ss.eps

convert fws_ss.eps fws_ss.png

18

19

残基数在y 轴，时间（ps ）在x 轴。看下面的NMR 结构：

从上面的dssp 图上，我们看到3个红色区域代表3个beta 片层（左图

中的黄色部分）。中间的较短区是最不稳定的。左边的图是用

pymol 程序(https://www.360docs.net/doc/207262485.html,/ )做的。

g_hbond

程序g_hbond 用于计算模拟过程中分子间或组间的氢键的数目，距离或角度。

g_hbond –f md.trr –s md.tpr –num fws_hnum.xvg

g_hbond 程序中的几何学关系

Gromacs 3.3.1的缺省值为：

r ≤0.35nm

α≤30°

用-r 和-a 标签设定其它限值。默认地，g_hbond 计算共体受体距离（r DA ）。你也可以用-da no 标签改变这种计算，设置 da 为no 则计算r HA 距离。

更多信息请参考Gromacs手册。

g_saltbr

用g_saltbr分析模拟中残基间的盐桥。程序会输出一系列*.xvg文件。给出-/-，+/-（最关注的）和+/+残基间的距离。

g_saltbr –f md.trr –s md.tpr

如何将特定时间点的轨迹保存成*.pdb文件：

要用特定帧（此例中为3000ps）代替整个轨迹，用-dump选项，如

trjconv -f traj.xtc -s file.tpr -o time_3000ps.pdb -dump 3000

遵照导师进行其它分析。

附录

如何重新进行损坏的计算

tpbconv -s prev.tpr -f prev.trr -e prev.edr -o restart.tpr

mdrun -s restart.tpr -deffnm myrestart (–deffnm将mdrun中的所有文件名设成默认名字。)

如何延长一个计算

tpbconv -f traj.trr -s topol.tpr -e ener.edr -o tpxout.tpr –time $VALUE –until $VALUE

其中$VALUE = ps（例如你要将2ns模拟延长到5ns，则$VALUE = 5000）

如何设置并行计算

grompp –np # –f md.mdp –c pr.gro –p fws.top –o md.tpr

20