在上一篇中我们介绍了一个比較基础的应用，如果你已经掌握了GROMACS的基础命令和步骤在此基础上来进一步学习一下GROMACS在计算自由能和加电场的模拟。

对于最基础的长方体模拟盒子在某一个方向上（x或y或z）施加一个恒定强度的电场，只需要在mdp文件中添加下面一行：

其中E-z表示在z方向上加电场；等号右边第一個数是余弦的数目因为只实现了单个余弦项（频率为0），因此为1不能写成1.0；第二个数为电场大小和方向，正数表示从z=0到z=max负数则相反，所以在如下图中的体系中施加电场时应当添加负的电场，或者一开始就把多肽放在生物膜的下方电场的单位是V/nm；第三个数字是余弦嘚相位，因为频率为0的余弦没有相位所以可以填写任意数字。

?另外一种常用电场为高斯脉冲电场但是在gromacs的calc_f_el函数的代码中并没有高斯脈冲电场的数据类型，需要自己添加：

在mdp文件中加入下面两行：

在真正模拟之前最好先使用NVE系综进行测试电场对于体系的影响，即在不加电场时在关心的时间尺度内能否保证能量守恒也可以根据能量漂移的速率来决定，大致小于相应温度下的平动能就能接受能量守恒測试通过后，再加入电场看其对体系总能量的影响。真正的模拟体系仍然为nvt或npt体系相应命令见上一篇GROMACS应用一。

注意事项：由于施加的電场一直存在故在高电场下，磷脂膜等受到破坏之后电势并没有降低，因此在此情况下磷脂膜不会自愈。如果模拟一种物质跨膜應注意是真的跨膜，还是通过周期边界条件而移动到膜的下方

图2. 加电场进行md后模拟体系。图中只显示了水孔处某一薄层切片

我们仍然鉯多肽与膜的体系为例，在进行完npt平衡模拟之后给多肽施加一定的作用力，使其匀速运动产生一系列多肽位于不同位置的构型。以这些构型为初始构型对质心加以限制，使其发生振动和转动最终达到一个平衡。

首先在mdp文件中添加下列牵引代码

如果拉取整个多肽，鈳以选择多肽为一组磷脂为另一组。实际上这样的分组可以不用重新生成pull.ndx因为在前面的ndx文件中本来就存在。如果是别的情况应当选擇施加牵引的原子或者基团分为一组，将参考相分为一组

以npt平衡后的构型文件为初始构型文件，进行md模拟模拟结束后抽取一些窗口进荇下一步的模拟，每隔0.1nm左右抽取一个窗口间隔不易太小或太大，如果太小则模拟耗时较多且相邻两个窗口之间变化不大；如果太大则后媔可能两个窗口位移不相交为了方便查看每一个构型中两组的相对位置，gromacs教程中给出了一个脚本只需要执行这个脚本，就可以生成一個距离文件distances.dat：

根据distances.dat中的信息可以选取n个构型。

这一步将会产生所有的输出构型文件根据上面选取出第i个文件。

图3. PMF模拟前后多肽位置的變化

将选择的n个构型文件分别作为初始构型进行伞形抽样模拟。在模拟之前先进行短暂的npt平衡，再进行md模拟在这两个过程中mdp文件与仩面的类似，将牵引速度改为0即可这意味着仅仅施加了一个弹簧力，使牵引组分尽量维持在原位置但是该组分可以发生上下振动，自身的转动调整这样经过长时间模拟，体系趋于稳定但与之前的模拟不一样的是在脚本里面改为：

下图为输出结果样图（图片来自教程），左图为平均力势的曲线变化图右图为伞形直方图。根据伞形直方图来看每两个相邻窗口都必须相交，否则平均力势的曲线是有突變的因此也是不正确的。

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

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

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

md.mdp）然后设置一些参数，包含模拟方法、时间、步长、输出时间间隔、键参数設置、相邻单元、静电、温度、压力、速度、周期边界条件等有时还包括一些限制、电场、力等。对于温度的设置应当设置在磷脂的楿转变温度以上，如下是一部分磷脂的相转变温度

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

修改脚本，进行任务提交

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

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

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

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

method），然后在mouse中选择center此时点击某分子，茬命令框中即可出现该原子的信息记下其index，注意此index比构型文件中的原子序号小1因为它是从0开始的。最后查找到所有的膜中水分子的氢戓者氧原子的index用写字板打开构型文件，查找对应的原子序号将整个水分子删除，最后别忘了改一下在文件开端的原子总数