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/）查看。

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：