GROMACS 是一个使用经典分子动力学理论研究蛋白质动力学的高端的高效的工具。GROMACS是遵守GNU许可的免费软件，可以从以下站点下载：http://www.gromacs.org，并且可以在linux和 Windows上使用。

在本教程中，将研究一个从漏斗形蜘蛛的毒液中分离的毒素。我们将使用显性溶剂动力学的方法来进行研究。首先比较真空中和溶解的模型。我们将把毒素肽溶在水盒子里，紧接着用牛顿运动定律加以平衡。我们还将比较偿离子在显性溶剂动力学中的影响。

1 下载pdb文件

1OMB.pdb (http://www.rcsb.org/pdb/)

2 用pdb2gmx 处理 pdb 文件

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

 

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  -bt  cubic –f fws.pdb  –o fws.pdb  –d 0.9

 

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

-d      决定了盒子的尺寸，即盒子边缘距离分子边缘 0.9nm （9Å）。理论上在绝大多数系统中，-d 都不能小于0.85nm。

 

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

例如：editconf  –f file.gro –o file.pdb  则将file.gro  转换为 file.pdb

 

现在就可以用产生的文件进行真空模拟了。真空模拟就是先能量最小化，然后进行动态模拟。 

4 在盒子中放入溶剂  

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

 

genbox命令在editconf产生的盒子基础上生成水盒子。上面的命令行指定了SPC水盒子。

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

5 设置能量最小化

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

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

 

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为particle mesh ewald 法（参见 Gromacs用户手册） 

rcoulomb   –指定库仑力阈值

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

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

EM Stuff

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

emstep     –初始步长（nm）

6 用grompp程序进行文件处理

grompp是预处理程序（即the gromacs pre-processor 的缩写） 

 

grompp –f em.mdp –c fws_b4em.pdb –p fws.top –o fws_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 –s fws_em.tpr –o fws_ion.pdb –nname CL- –nn 2 –g fws_ion.log 

-nname   指定阴离子名称（在Gromos G43a1力场中，用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文件的温度耦合参数

; 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

 

记住：如果要加入氯离子，需要重新运行第6步的grompp。首先删除旧的 fws_em.tpr 文件，然后运行下面的 grompp命令：

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

9 在后台运行能量最小化

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

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

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

tail  –15 em.log

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

用tail -50 em.log ：

设置位置限制性动力学模拟、非限制性动力学模拟以及对结果分析详见收费内容。

下图是基于Pymol程序制作的NRM结构图：