在上一篇我们讲了GROMACS的简介和安装，在这一篇中我们讲一下GROMACS的使用。GROMACS中含有精准且全面的立场（尤其是氨基酸），这使其被广泛用于生物模拟中，它不仅可以进行全原子模拟计算，还可以进行粗粒化计算。在这我们将分三个板块举例介绍GROMACS的使用：GROMACS的简单使用，GROMACS计算自由能、加电场体系以及非生物体系的使用，GROMACS用于Martini立场粗粒化计算。GROMACS在生物的研究中较为专业，故针对GROMACS的使用以多肽与生物膜为例。在第一部分，先来细致的介绍一下GROMACS的基础应用。

 

首先我们先来了解一下GROMACS的输入输出文件：

 

1. 构型文件：

       要模拟一个体系必须有相应物质的构型，比如磷脂双分子层，蛋白质等。对于蛋白质，多肽等可以从蛋白质数据库中搜索下载（网址http://www.rcsb.org/pdb/home/home.do），新手一般根据文献中模拟过得多肽代称进行搜索（比如图1左边，HIV-Tat代称是1JFW）。磷脂双分子层膜可以在lipidbook中下载（网址https://lipidbook.bioch.ox.ac.uk/lipid/或者http://wcm.ucalgary.ca/tieleman/downloads）；含胆固醇的磷脂膜可以在Computational Molecular Biophysics Group中下载（网址http://cmb.bio.uni-goettingen.de/cholmembranes.html）；此外，在MemBuilder中也可以根据自己的要求免费生成相应的磷脂膜（网址为http://bioinf.modares.ac.ir/software/mb2/builder.php），包括混合膜，如图1右边。其他的一些分子可以在ATB中下载。或者用一些软件也可以实现，比如vmd，MS等。对于碳纳米管，可以在gromacs主页上搜索Carbon Nanotube，有相应的模拟解说，个人觉得最简单的方法是在网站http://turin.nss.udel.edu/research/tubegenonline.html上生成。关于水分子和离子等一些小分子可在模拟过程中通过命令来直接添加。

图1.  蛋白质数据库

 

图2.  Membuilder界面

 

2. 拓扑文件：

       只有构型文件是不够的，因为构型文件里面只有这些分子的各个原子的顺序和坐标，并没有它们之间的连键情况等，故进行模拟还需要体系中各成分的拓扑文件。对于蛋白质的拓扑文件，可以在gromacs中通过命令pdb2gmx来生成，因为gromacs中包含很多种氨基酸的拓扑文件和转换脚本。对于其他分子的拓扑文件，在下载构型文件时，往往会包含相应的itp文件，即为分子的拓扑文件。水分子和离子等小分子的拓扑文件在gromacs的立场里面都有，可直接调用。

 

3. 参数文件：

       一般在gromacs tutorial中都可以直接下载（后缀为mdp的文件，包括ions.mdp, minim.mdp, nvt.mdp, npt.mdp, md.mdp），然后设置一些参数，包含模拟方法、时间、步长、输出时间间隔、键参数设置、相邻单元、静电、温度、压力、速度、周期边界条件等，有时还包括一些限制、电场、力等。对于温度的设置，应当设置在磷脂的相转变温度以上，如下是一部分磷脂的相转变温度。

 表1.  部分磷脂的相转变温度

 

4. 立场文件：

       GROMACS安装后的文件夹内部路径/pwd/gromacs-x1.x2.x3/share/gromacs/top中包含很多立场，有amber、gromos、opls/aa、charmm，martini立场等，但是有时候我们下载的磷脂拓扑是基于Berger立场而言的，因此运用gromacs中的立场时，应在立场文件中添加相应的非键参数。但是有网站下载的磷脂构型则不需要。对于多肽，我们选取gromacs53a7立场，避免了多肽在水中长时间模拟后解螺旋的问题。此外，立场文件中还包含了水和离子等相关参数。

 

注意事项：我们应当保证拓扑文件和构型文件的原子名称一致，拓扑文件和立场文件中原子类型一致。

 

在准备好输入文件之后就可以进行模拟了，具体步骤如下：

1. 整个体系的构造

（1）将多肽的pdb文件转换成gro文件，同时生成其拓扑文件。

gmx  pdb2gmx  –f  protein.pdb  –o  protein.gro  –water  spc (其中水的类型可以选择spc、tip3p等)

       在模拟蛋白质的时候，要考虑氢原子、蛋白质的两端以及带电残基等，可以用-ingh(忽略氢原子), -ter(选择多肽的两端氨基和羧基的带电情况), -inter(交替分配电荷状态)来实现。在这一步中将会选择立场，蛋白质的两端。完成该步之后生成protein.gro, topol.top, posre.itp，其中topol.top为多肽的拓扑文件，在其中引入了立场，多肽的原子性质，连键情况、形成的角等，并且引入了水和离子的拓扑文件；posre.itp是蛋白质的重原子，用于后面对蛋白质进行位置限制。实际上总的拓扑文件里只要包含立场、各个分子的结构性质即可，无论是直接写在里面也好，还是通过include来引入也好。

 

（2）给生成的多肽构型定义一个盒子（根据整个体系的盒子尺寸来构造），并且将多肽放在合适的位置上。

gmx  editconf  –f  protein.gro  –o  pro_newbox.gro  –box  *  *  *  -center  *  * *

       盒子的大小要取决与膜的大小（见下面步骤）。为了将多肽放在膜上方来进行模拟，所以在此定义的位置应该与后面生物膜的位置有一定距离。该步生成pro_newbox.gro。

 

（3）给生物膜定义一个同样的盒子，并放置在相同的位置上。

gmx  editconf  –f  bilayer.pdb  –o  bilayer_newbox.gro  –box  *  *  *  -center  * *  *

       注意双层膜的位置不要与多肽重叠。并且盒子的长和宽要与膜的大小相匹配，否则后面加水后磷脂的疏水尾部也将暴露在水中。如果不知道下载的膜的长宽，可以先根据经验定义一个盒子，然后用vmd打开后，输入指令pbc box，就可以显示盒子的框架了，看看是否适合，不合适就进行调整。关于vmd的使用，在本文的最后会有详细介绍。如果下载的磷脂膜模型较小，可以用cat命令将两个放在一起，具体使用方法类比下面cat命令介绍。或者用如下命令：

genconf  -f  bilayer_newbox  -nbox  2  1  1  -o   bilayer_big.gro

其中-nbox后面的数字分别是在x、y、z方向上复制的倍数。得到的新的磷脂膜最好经过一段时间的平衡，使其连接处更加连贯。

图3.   vmd中显示盒子

 

（4）将多肽与膜放在一起，并且添加水分子。

cat  pro_newbox.gro  bilayer_newbox.gro  >  whole.gro

vi  whole.gro  (打开whole.gro编辑，将蛋白质最下方的盒子大小，以及双层膜上方的粒子总数等均删除，使蛋白质与磷脂膜原子相衔接)

vi  topol.top  (编辑拓扑文件，将磷脂分子加入)

       在加水之前，应当将vdwradii.dat中的C的值从0.15改为0.375，这样是为了防止加水后生物膜的空隙中也有水，如果此时还有极少量的水分子在膜空隙中，可以用写字板等直接编辑生成的gro文件，删除其中的水分子，具体步骤祥见下面的vmd部分。删除水分子之后，别忘了top文件里面相应的水分子数目要减少。

gmx  solvate  –cp  whole.gro  –cs  spc216.gro  –o  sol.gro  –p  topol.top 

  

图4.  模拟中有四个多肽，288个DPPC磷脂分子，并且引入了磷脂限制的判定

 

（5）添加离子，用于平衡电荷或模拟更真实环境

gmx  grompp  -f  ions.mdp  -c  sol.gro  -p  topol.top  -o  ions.tpr

gmx  genion  -s  ions.tpr  -o  system.gro  -p  topol.top  -pname  NA  -nname  CL  –np  10  –nn  18

       其中-nn是添加负电荷CL-的数目，-np是添加正离子NA+的数目，应该根据具体情况而定，但是必须保证整个体系最终是电中性的，否则模拟结果不可信。

 

2. 开始模拟

整个体系的构型已经搭建完毕，接下来就开始步入模拟了。

（1） 能量最小化

gmx  grompp  –f  minim.mdp  –c  system.gro  –p  topol.top  –o  min.tpr

        接下来就可以任务提交了，前面的步骤也可以用windows版gromacs实现。在提交任务时，需要用到脚本，内容如下：

 

       输出文件包括构型文件（min.gro）、日志文件（min.log）、轨迹文件（min.trr、min.xtc）。可以通过修改参数文件内的输出时间间隔来控制输出，最好不要过于细致（占内存较大），也不要过于粗略（不利于分析）。

 

（2） 建立index文件

gmx  make_ndx  -f  min.gro  -o  index.ndx

       在建立过程中，会选择分组，比如，我们可以将水和离子分为一组，根据提示可以将水和离子前面的序号用“|”连接即可，如14|16。

 

（3） nvt模拟

       在nvt模拟之前，我们会发现模拟体系中磷脂头部与水分子是有一定距离的，为了更好的使磷脂膜溶在水中，我们进行同时限制蛋白质和磷脂分子的nvt模拟。在mdp文件中的define一项中写入-DPOSRES  –DPOSRES_LIPID，与之呼应的是在top文件中应当添加相应的内容：

#ifdef  POSRES_LIPID

#include  "lipid_posre.itp"

#endif

       其中磷脂的限制文件lipid_posre.itp中只需限制磷原子的位置即可。此外，nvt平衡中通常使用的热浴为Berendsen或者V-rescale，使最终温度能基本达到稳定。

gmx  grompp  –f  nvt.mdp  –c  min.gro  –p  topol.top  –n  index.ndx  –o nvt.tpr

       类似能量最小化，只须将脚本中的min改为nvt即可。输出文件比能量最小化多了一个断点文件（nvt.cpt）

 

（4） npt模拟

       nvt模拟之后温度基本达到稳定，但是还需将模拟盒子根据设定的压力进行平衡模拟以松弛磷脂膜，因此进行npt模拟只需限制蛋白质即可。但还应注意的是，在npt模拟中一般使用Nose-Hoover热浴，与nvt中的两个热浴不同，该热浴允许比较大的涨落，因此不适用于nvt模拟。对于压力运用Parrinello-Rahman方法，由于模拟盒子中有磷脂膜，所以我们在此使用semiisotropic，与膜平行的两个方向为一个压力，与膜垂直的方向为一个压力。

gmx  grompp  –f  npt.mdp  –c  nvt.gro  –t  nvt.cpt  –p  topol.top  –n  index.ndx  –o  npt.tpr

将脚本中nvt改为npt，即可提交任务。

 

（5） md模拟

       平衡了温度和压力之后，开始较长时间的md模拟，模拟过程与npt类似。

如果计算结束以后，觉得仍未平衡，或者想要更长时间的模拟，可以进行续算，一种方法是类似上面，只是将md.gro作为初始构型，进行md模拟；另一种方法是延长计算时间，用-extend 10000延长10000ns。

gmx  grompp  –f  md.mdp  –c  npt.gro  –t  npt.cpt  –p  topol.top  –n  index.ndx  –o  md.tpr

修改脚本，进行任务提交。

 

3. 分析

       关于分析，gromacs具有很多命令，可以直接进行分析，可根据模拟需求来进行选择。在分析之前，应当先建立分组索引文件，再对其中的某一组或者几组进行分析，包括能量(gmx energy)、密度(gmx density)、径向分布函数(gmx rdf)、均方位移(gmx msd)、氢键(gmx hbond)、回旋半径(gmx gyrate)、蛋白质的二级结构(gmx do_dssp)、静电势(gmx potential)等。同时，gromacs还提供了查看轨迹文件工具 gmx view，此方法无需OpenGL支持。比如对于磷脂膜密度的分析，在做分析之前要建立相应的index文件：

gmx  make_ndx  -f  md.tpr  -o  density_groups.ndx

...

 > 12 & a C1 | a C2 | a C3 | a N4 | ... | a O11

 > name 22 Headgroups

 > 12 & a C12 | a C13 | a O14 | ... | a C50

 > name 23 Tails

 > q

选择的原子取决于磷脂的构型。

       更多的分析可参见gromacs tutorial，gromacs manual，文献等。值得注意地是，在分析的时候经常会用到一段时间的平均值，在选取时间时（-b 起始时间-e终止时间）应当避免前面一小段不平衡的时间。如果是计算了两个文件，也可以将两个文件合并（gmx trjcat）之后再做分析。

 

注意事项：

       如果模拟过程中报错，应当查找错误信息再进行下一步，如果不明白错误信息，可以将错误信息复制到gromacs主页搜索栏进行搜索；如果出现警告，则应当查看警告内容，比如体系整体显电正性，则必须改正，有的警告则在接受范围内，可以忽略，在grompp是加-maxwarm 1；旧版gromacs命令前面没有gmx，另外加水的命令不是solvate，是genbox。

 

       讲到gromacs的使用，不得不提到一个用于查看gro或者pdb文件的软件vmd。vmd的功能也有很多，在此仅介绍与gromacs模拟相关的、常用的功能。

1. 显示

       用vmd打开构型文件，我们可以分组显示（graphics>representations>selection>resname）,也可以改变显示的颜色（graphics>colors或graphics>representations>draw style）和类型（graphics>representations>draw style），也可以改变背景色（graphics>colors>display）。不仅如此，vmd还可以显示一部分原子或者分子，比如存在于膜中的水分子，在分组显示之后，我们可以在水的分组中输入and z>30 and z<80等

 

2. 查找

       上面提到，当磷脂膜中还存在几个水分子时，我们可以手动删除，这时候可以借助vmd快速查找需要删除的原子或者分子。通过上述查看的方法，我们可以快速找到膜中水分子，将其显示为VDW模式（graphics>representations>draw style>drawing method），然后在mouse中选择center，此时点击某分子，在命令框中即可出现该原子的信息，记下其index，注意此index比构型文件中的原子序号小1，因为它是从0开始的。最后查找到所有的膜中水分子的氢或者氧原子的index，用写字板打开构型文件，查找对应的原子序号，将整个水分子删除，最后别忘了改一下在文件开端的原子总数。

 

3. 查看轨迹文件

       用vmd打开构型文件，并导入相应的trr文件，即可播放动画

 

4. 输出

       将最终想表达的图片文件输出 file > render