GROMACS使用教程要点Word文档下载推荐.docx

GROMACS使用教程要点Word文档下载推荐.docx

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

用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）（下载地址：

//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≤