GROMACS使用教程要点.docx

GROMACS使用教程要点.docx

《GROMACS使用教程要点.docx》由会员分享，可在线阅读，更多相关《GROMACS使用教程要点.docx（19页珍藏版）》请在冰豆网上搜索。

GROMACS使用教程要点

GROMACS教程

GROMACS是一个使用经典分子动力学理论研究蛋白质动力学的高端的高效的工具。

GROMACS是遵守GNU许可的免费软件，可以从以下站点下载：

http:

//www.gromacs.org，并且可以在linux和Windows上使用。

在本教程中，将研究一个从漏斗形蜘蛛的毒液中分离的毒素。

我们将使用显性溶剂动力学的方法来进行研究。

首先比较真空中和溶解的模型。

我们将把毒素肽溶在水盒子里，紧接着用牛顿运动定律加以平衡。

我们还将比较偿离子在显性溶剂动力学中的影响。

更全面的用法指导请参考官网的GROMACS用户手册http:

//www.gromacs.org

注意：

在本教程中，将要生成的gromacs（*.gro）结构文件，可以用VMD（下http:

//www.ks.uiuc.edu/Research/vmd/）查看。

一Gromacs基本模拟流程

1下载pdb文件

1OMB.pdb（http:

//www.rcsb.org/pdb/）

2用pdb2gmx处理pdb文件

pdb2gmx–ignh–ffG43a1–f1OMB.pdb–ofws.pdb–pfws.top–waterspce

pdb2gmx此命令将pdb文件转换成gromacs文件并产生拓扑文件。

-ignh因为本pdb文件是由NMR产生的，含有氢原子，因此用-ignh选项忽略文件中的氢原子。

-ff指定力场（G43a1是Gromos96力场，一个通用原子力场）。

-f读入pdb文件，

-o指定一个新产生的pdb文件（也可以是其它多种类型文件）的文件名。

-p指定新产生的拓扑文件名。

拓扑文件包含了所有力场参数（基于一开始选择的力场），因此非常重要。

-water来指定水模型研究表明SPC/E水模型在水盒子模拟中表现最好。

用SPC/E水模型研究长程静电相互作用较好。

#注：

对于下面将要用到的任何命令，都可以使用“-h”查看该命令的使用方法，比如，对于命令pdb2gmx可以使用：

pdb2gmx–h

3建立盒子

editconf-btcubic–ffws.pdb–ofws.pdb–d0.9

用上面的命令建立了一个简单的立方体盒子.

-d决定了盒子的尺寸，即盒子边缘距离分子边缘0.9nm（9Å）。

理论上在绝大多数系统中，-d都不能小于0.85nm。

注：

editconf也可以用来进行gromacs文件（*.gro）和pdb文件（*.pdb）的相互转化。

例如：

editconf–ffile.gro–ofile.pdb则将file.gro转换为file.pdb

现在就可以用产生的文件进行真空模拟了。

真空模拟就是先能量最小化，然后进行动态模拟。

4在盒子中放入溶剂

genbox–cpfws.pdb–csspc216.gro–ofws_b4em.pdb–pfws.top

genbox命令在editconf产生的盒子基础上生成水盒子。

上面的命令行指定了SPC水盒子。

genbox命令可以在给定尺寸的盒子中加入正确数目的水分子。

5设置能量最小化

em.mdp文件：

Gromacs用*.mdp文件指定所有计算的参数。

它用最速下降法消除原子位置碰撞。

编辑文件，将nsteps变成400。

如果最小化不能收敛，就用nsteps=500再做一次。

（最小化在400步内一般是能收敛的，但不同的平台可能结果会不一样。

）要重做的话，必须重新运行grompp（注意：

预处理器的位置在你的机器上可能不同，用which命令来定位，即whichcpp）

em.mdp文件内容：

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

cpp–指定预处理器的位置

define–传递给预处理器的一些定义。

–DFLEXIBLE告诉grompp将spc水模型而非刚性SPC包含进拓扑文件，以便用最陡下降法进一步最小化能量。

constraints–设置模型约束

integrator–steep，告诉gompp使用最速下降法进行能量最小化。

cg则代表使用共轭梯度法。

dt–能量最小化用不用。

只在动力学模拟中用（如md）。

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

nstlist–更新邻居列表的频率。

nstlist=10表示每10步更新一次。

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

coulombtype–告诉gromacs如何计算静电。

PME为particlemeshewald法（参见Gromacs用户手册）

rcoulomb–指定库仑力阈值

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

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

EMStuff

emtol–最大的力如果小于此值则能量最小化收敛（结束）（单位kJmol–1nm–1）

emstep–初始步长（nm）

6用grompp程序进行文件处理

grompp是预处理程序（即thegromacspre-processor的缩写）

grompp–fem.mdp–cfws_b4em.pdb–pfws.top–ofws_em.tpr

-f标签指定输入参数文件（*.mdp）。

-c输入结构文件（pdb文件，*.pdb）;

-p输入拓扑文件

-o输出mdrun的输入文件（*.tpr）。

7使用genion和tpr文件添加离子

对生成的tpr文件加入补偿离子以中和系统中的净电荷。

我们的模型中有+2.00静电，因此加入两个氯离子。

将fws_em.tpr文件拷贝到“ionwet”子目录，并且将fws.top和posre.itp拷贝到这个目录。

用genion命令添加氯离子：

genion–sfws_em.tpr–ofws_ion.pdb–nnameCL-–nn2–gfws_ion.log

-nname指定阴离子名称（在GromosG43a1力场中，用CL-表示氯离子。

参见ions.itp查看wrt力场中离子详细信息）

-nn是指定加入的阴离子数目。

-g输出genion的log文件。

运行这个命令时，提示提供一个连续的溶剂组，应该是组12（SOL）。

输入12，回车。

程序会告知你有两个溶剂分子被氯离子代替。

现在你必须修改fws.top文件：

添加

#include“ions.itp”（注意：

3.2及以后版本会自动添加）

经过包含声明后，力场在最后减掉两分子SOL，加入两分子Cl。

8用fws_ion.pdb来产生能量最小化的输入文件

你还需要修改pr_md.mdp和md.mdp两个文件中的温度耦合参数。

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

;Berendsentemperaturecouplingusingvelrescalingison

Tcoupl=v-rescale

tau_t=0.10.1

tc_grps=proteinnon-protein

ref_t=300300

记住：

如果要加入氯离子，需要重新运行第6步的grompp。

首先删除旧的fws_em.tpr文件，然后运行下面的grompp命令：

grompp–fem.mdp–cfws_ion.pdb–pfws.top–ofws_em.tpr

9在后台运行能量最小化

nohupmdrun–v–sfws_em.tpr–ofws_em.trr–cfws_b4pr.pdb–eem.edr–gem.log&

nohup...&使任务后台运行

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

tail–15em.log

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

用tail-50em.log：

二设置位置限制性动力学模拟

什么是位置限制性模拟？

你限制（或部分冻结）大分子中的原子位置，而允许溶剂分子运动。

这样做像是将水分子浸入大分子。

水分子松弛时间约为10ps。

因此我们要进行超过10ps的位置限制性模拟。

本实例中用20ps，大的模型（大蛋白或脂）可能需要更长的平衡时间，50ps或100ps或更长。

下面的设置在这个gromacs力场中运行良好。

其他力场请参考用户手册（例如在GROMOS96力场中，建议nstlist=10andrvdw=1.4）。

在coulombtype,PME代表“ParticleMeshEwald”静电势。

PME是计算长程静电势的最优算法（给出最可信的能量评估，尤其在用Na+,Cl-,Ca2+等作为补偿离子的体系）。

由于这个蛋白具有暴露的带电残基，使系统带有+2静电荷，所以适用PME算法，更为有益的是用补偿离子使系统处于电中性。

constraints中的all-bonds选项可以应用线形限制算法确定系统中的所有键长（当dt>0.001ps时尤为重要）。

学习一下下面的mdp文件。

pr.mdp:

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

constraints声明如前所述。

all-bonds设定LINCS算法限制所有键。

integrator告诉gromacs进行何种动态算法（另外的选项“sd”代表stochasticdynamics）

dt是每步的时间（我们选择了2fs；但此处的单位一定是ps！

）

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

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

coulombtype选择gromacs计算原子静电相互作用方法（PME代表particlemeshewald；另外还可以用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.

当你改变温度时，别忘了改变gen_temp变量以生成速度。

pcoupl–Parrinello-Rahman恒压器。

pcoupltype–isotropic指“box”可以平均地向各个方向（x,y,z）膨胀或压缩，来维持一定的压力。

注意：

进行膜模拟时用semiisotropic。

tau_p–压力耦合的时间常数（单位ps）。

compressibility–溶剂在每bar的可压缩性（上面的设置是水在300K和1大气压下的可压缩性）。

ref_p–压力耦合的参照压力（单位bar，1大气压~0.983bar）。

grompp–fpr.mdp–cfws_b4pr.pdb–pfws.top–opr.tpr–maxwarn3

nohupmdrun–spr.tpr–opr.trr–cfws_b4md.pdb–epr.edr–gpr.log&

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

三设置非限制性动力学模拟

md.mdp文件和pr.mdp文件相仿。

有几处不同，define声明不再需要，因我们不再做位置限制模拟。

用于explicitsolvation的md.mdp文件内容（特殊注释：

做真空模拟，去掉温度耦合中的“sol”部分。

在有补偿离子的模拟中，为离子加入相应的温度耦合参数）

grompp–fmd.mdp–cfws_b4md.pdb–pfws.top–omd.tpr–maxwarn3

nohupmdrun–smd.tpr–omd.trr–cfws_md.pdb–emd.edr–gmd.log&

用tail命令查看md.log文件。

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

trjconv–fmd.trr–omd.xtc

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

四并行计算、重启、延长模拟计算

1如何重启一个计算

tpbconv-sprev.tpr-fprev.trr-eprev.edr-orestart.tpr

mdrun-srestart.tpr-deffnmmyrestart（–deffnm将mdrun中的所有文件名设成默认名字）

2如何延长一个计算

tpbconv-ftraj.trr-stopol.tpr-eener.edr-otpxout.tpr–time$VALUE–until$VALUE

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

3如何设置并行计算

grompp–npN–fmd.mdp–cpr.gro–pfws.top–omd.tpr

-np标签设定并行计算的节点数，N代表节点数，自行设定。

然后用mdrun_mpi进行并行任务：

mpirun–np#/products/gromacs/bin/mdrun_mpi–deffnmmd

五模拟结果分析

1如何将特定帧的轨迹保存成*.pdb文件

用-dump选项，如

trjconv-ftraj.xtc-sfile.tpr-otime_3000ps.pdb-dump3000

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

ngmx–fmd.trr（ormd.xtc）–smd.tpr

当观察器启动后，将看到一个多选项的对话框。

选择标“protein”的多选框，点击OK。

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

用X-Rotate上下旋转盒子（鼠标左键向上，右键向下）。

用Y-Rotate左右旋转盒子（左键向左，右键向右）。

最下面的Scale用来放大或缩小视图（左键放大，右键缩小）。

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

要观察模拟轨迹动画，点击Display>Animate。

动画播放控制在窗口的底部。

点击中间的箭头按钮逐帧观看。

点向前的双箭头观看整个轨迹动画，点暂停按钮停止动画。

点向左的双箭头按钮重置动画。

保存并观察*.pdb文件最好的方法是用visualmoleculardynamics（VMD）（下载地址：

http:

//www.ks.uiuc.edu/Research/vmd/，它是学术免费，且在linux和windows下运行）结果分析Gromacs的一个主要优点是有一系列分析轨迹文件的小程序。

3比较常用的分析工具

3.1组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–fclg_b4md.pdb–oclg_ter.ndx

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

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

注意：

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

新建索引组的缺省名字（r_1_36_37_72_73_108）很繁琐，可以用name命令修改。

我们在

命令中用索引号#（15）。

用“v”命令查看名字是否改成功了。

用“q”保存并退出。

现在怎么固定组呢？

简单，在md.mdp文件中加入下面几行：

记住当用新的mdp文件时，首先用grompp将新索引文件加入tpr文件。

用grompp的-n标签，例如：

grompp–fmd.mdp–cpr.gro–pclg.top–nclg_ter.ndx–omd.tpr

3.2特性研究g_confrms 

要比较最后结构和初始PDB文件的差异，用g_confrms（用g_confrms–h查看详细信息）。

此程序计算两个结构的最小二乘拟合。

g_confrms–f11OMB.pdb–f2md.gro–ofit.pdb

你将被提示选择一个组（两次都选（组4））。

程序将报告RMSD值，并产生一个输出文件（fit.pdb）。

输出文件中包含两个位置重叠的结构。

3.3g_covar计算斜方差

也可用于从动态轨迹计算平均结构。

如计算1ns动态模拟的后200ps的平均结构：

g_covar–ftraj.xtc–stopol.tpr–b801–e1000–avtraj_avg.pdb

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

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

g_energy–fmd.edr–ofws_pe.xvg

首先要选择输出（*.xvg）的数据。

输出文件是一个电子数据表，可以用Xmgr或Grace打

开。

它是一个文本文件，在进行一些小的改动后可以用MicrosoftExcel打开。

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

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

再按一次回车

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

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

xmgrace-nxyfws_pe.xvg

可以到以下地址下载Gracehttp:

//plasma-gate.weizmann.ac.il/Grace/.Grace只能在linux和unix上运行。

如果没有Grace或Xmgr，可以作为空格分隔文件导入MSExcel。

3.5g_gyrate测量回旋半径

这个指标用于度量结构的紧密度。

此程序计算某（些）原子质量与分子重心的关系。

g_gyrate–fmd.trr–smd.tpr–ofws_gyrate.xvg

3.6g_rms与g_rmsdist计算结构的RMSD值

用g_rms计算动态模拟过程中的结构与初始结构的结构偏差。

（-dt10选项告诉程序每10帧计算一次）

g_rms–s*.tpr–f*.xtc–dt10

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

g_rms–sem.tpr–fmd.trr–ofws_rmsd.xvg

选择4组（Backbone）计算最小二乘拟合。

程序生成一个rmsd随时间变化的图（rmsd.xvg）。

以空格分隔文件的形式导入Excel。

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

例如：

上面的实例是一个1ns模拟（你的结果可能不同）。

这个模拟需要延长到完全平衡。

3.7g_rmsf计算原子位置的根均方波动（rmsf）

与g_covar相似，此程序也可以计算平均结构。

例如，计算一个2ns（2000ps）模拟的后500ps的平均结构，用如下命令：

g_rmsf–ftraj.xtc–stopol.tpr–b1501–e2000–otraj_rmsf.xvg–oxtraj_avg.pdb

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

命令如下：

g_rmsf–smd.tpr–fmd.trr–b200–e500–oxfws_avg.pdb

提示时选择组1“Protein”。

推荐的真空能量最小化的em.mdp文件。

先用最陡下降法，再用共轭梯度法。

警告！

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

程序g_rmsf也可以用来计算温度因子。

计算的温度因子可以和X光晶体结构的温度因子比较。

g_rmsf–smd.tpr–fmd.xtc–ormsf.xvg–oqfws_bfac.pdb

仍选择“Backbone”组。

3.8do_dssp计算模型的二级结构

前提是你必须在电脑中（/usr/local/bin）安装了dssp程序（http:

//swift.cmbi.ru.nl/gv/dssp/）。

do_dssp–smd.tpr–fmd.trr–ofws_ss.xpm

选择计算组1（Protein）。

用xpm2ps将xpm文件转成eps格式。

然后用ImageMagick转化程序将eps文件转成png文件或其他格式文件。

xpm2ps–ffws_ss.xpm–ofws_ss.eps

convertfws_ss.epsfws_ss.png

残基数在y轴，时间（ps）在x轴。

看下面的NMR结构：

从上面的dssp图上，我们看到3个红色区域代表3个beta片层（下图中的黄色部分）。

中间的较短区是最不稳定的。

下边的图是用pymol程序（）做的。

3.9g_hbond计算模拟过程中分子间的氢键的数目、距离或角度

g_hbond–fmd.trr–smd.tpr–numfws_hnum.xvg

Gromacs4.0的缺省值为：

r≤