gromacs文件介绍and一些杂知识

（1）gromacs(GMX) 各种文件格式详细，可以查阅GROMACS 手册第5章第6小节，以下为简要介绍。

CPT文件：该文件为模拟断点文件（check point，.cpt）。该文件为模拟过程固定时间间隔产生，保存模拟系统所有信息。该文件一部分可以在能量文件（.edr）找到，一部分可以在双精度轨迹文件（.trr）中找到。如果模拟不幸因为外界条件中断（如断电，模拟人发脾气砸电脑等），可以使用该文件重新在断点处开始模拟，以节省模拟时间。同时也可以依靠该断点文件开始，并延长模拟计算（见tpbconv）。

EDR文件：系统能量文件（energy，.edr）。该文件记录模拟输入文件中定义的能量组的各种相互作用能量等。

EPS文件：封装文件格式（.eps），并不是GROMACS自身文件格式，可以当图片打开。LINUX系统下一般已经有默认打开程序，WINDOWS要安装其他打开程序（可以GOOGLE以下）。GROMACS的DSSP和罗麽占陀罗图等通过xpm2ps处理后都是这个文件格式。习惯就好。

G87文件：分子坐标文件（.g87）。该文件记录并只记录原子坐标和速度，不含原子序号。并只记录常压强模拟系统的盒子信息。

G96文件：分子坐标文件（.g96）。GROMOS96程序的分子坐标文件，模拟程序以15.9的C语言格式写入，精度较高，但是会比较大。包含有文件头，时间步，原子坐标，原子速度，以及盒子信息等。

GRO文件：分子坐标文件（.gro）。GROMACS的最主要分子坐标文件，明白这个文件，就基本明白使用GROMACS了。该文件类型的各个文本列

字

数

固

定

，

C

语

言

的

写

入

格

式

为

：

\。具体固定文本列有：残基序号，5位数；残基名称，5字母；原子名称，5字母；原子序号，5

1

为数；原子坐标三列，X，Y，Z坐标各8位数，含3个小数位；速度同坐标，速度单位为nm/ps（km/s）。

ITP文件：分子拓扑文件（.itp）。被主拓扑文件（.top）包含的分拓扑文件，一般包含某个特定分子的类型。于主拓扑文件区别有它不引用其他力场文件，同时包含[system]，[molecule]等拓扑字节。

M2P文件：xpm2ps程序配置文件，定义输出eps文件中颜色，字体种类及大小等。

MDP文件：GROMACS的模拟配置文件（.mdp）。该文件所含定义较多，各关键字的含义可以查阅GROMACS手册。（这是使用GROMACS进行分子动力学模拟最最最最（10个最）重要的文件，no mdp文件，no GROMACS模拟。好好看书，以明白各个关键字的含义。因为它太重要，所以不在此简要描述。

N2T文件：原子名称及类型对照文件（.n2t）。x2top程序可以按照原子名称得到该原子的原子类型力场参数，N2T就是x2top程序扫描的数据库，文件很小。文件中文本行有原子名称，原子类型，原子电量，原子质量，该原子与其他原子成键距离等。

NDX文件：原子索引文件（.ndx）。该文件含原子的序号，当使用make_ndx程序生成索引文件时，可以定义不同的原子组，每组名下即是该组所含各个原子的序号。

PDB文件：分子坐标文件（.pdb）。这个就不用说了（说真的，如果真没有听过这个文件类型的话，看这篇文章有点浪费时间。）

RTP文件：残基力场参数文件（.rtp）。该文件包含常见残基的力场信息，包括残基所含原子，成键种类等。使用pdb2gmx处理PDB文件时，程序按照PDB文件信息，在RTP文件中寻找对应的残基力场信息。

2

TOP文件：模拟系统的拓扑文件（.top）。该文件就是所谓十分及其著名的系统拓扑文件啦，其包含各个关键字都十分易懂；一般其还包含引用其他力场文件（#include）。TOP文件一般由pdb2gmx产生，grompp程序生成模拟TPR文件时使用。

TPR文件：模拟打包文件（.tpr）。该文件打包模拟需要各种信息，包括模拟系统，模拟控制等。

TRJ文件：全精度轨迹文件（.trj）。该文件包含模拟系统模拟各个时间下的原子坐标，速度和受力等。所含帧数频率由MDP文件控制，文件较大。 TRR文件：以上同，一般为默认格式。由于所含信息多，可以也EDR文件一起使用，重新开始模拟程序。

XPM文件：数据矩阵文件（.xpm）。该文件矩阵中每个值即是矩阵点所表示的物理量大小（也可以是布尔值）。该文件其实就是二维图，可以失踪xpm2ps转换为图片。

XTC文件：模拟轨迹单精度文件（.xtc）。单精度轨迹文件，文件较TRR和TRJ小，为常用分析文件。包含模拟系统中原子坐标，模拟时间，和模拟盒子信息。

XVG文件：二维图标文件（.xvg）。二维画图工具xmgrace的默认文件，可以使用xmgrace打开。 （2）Gromacs中几个特殊文件 aminoacids.dat

该文件保存GMX默认的蛋白质和核算的默认残基名称。如果计算过程要建立一个新的蛋白质或者核算残基，可以将新的残基名称加到该文件中，并增加文件第一个的整数即可。有时候可以将该文件拷贝到当前工作文件夹进行编辑，以不影响其他计算的命名（GMX的文件搜索总是从当前目

3

录开始的。） FF.dat

GMX默认力场列表，即pdb2gmx处理PDB文件时可以选择的立场列表。增加新的力场，可以编辑该文件，并修改文件第一行的整数，使其与力场种类熟目一致。 specbond.dat

GMX处理特殊化学键的文件，特殊化学键包括二硫键，血红素铁原子于其他原子成键等。该文件第一行指明特殊键对的数目，第二行开始即为各个特殊键对的信息，其中第一列为键对第一个残基的名称，第二列为该残基成键原子的名称，第三列为该原子可以成键的数目，第四到第六列为成键另一个残基的信息，第七列为该化学键的平衡长度，此后两列为成键后残基的新名称。 vdwradii.dat

原子范德华半径数据库。使用genbox为系统添加水分子，或者使用genion为系统添加离子时，各个原子间的距离要大于两个原子范德华半径之和，否则则为原子重叠 （3）常见水分子模型

进行分子动力学模拟，水分子十分重要，除非选择使用连续介质模型（implictit water model）。水分子模型较多，选择这些模型要结合使用的力场，并参考别人已经的数据。一下简单介绍几种常见的水分子模型，希望对了解它们有点帮助。

按照一般化学常识，水分子由三个原子构成，主要的参数应该有各个原子的质量，电量，氢氧键的长度以及H-O-H的键角。没有错，最简单的水分子模型就是这些参数都固定的刚性水分子模型。如SPC模型和TIP3P

4

模型。这两种模型中，原子质量和电量都在同一个质点上。唯一不同的是TIP3P的H-O-H键角比理论值109.47小，为104.52度。这两种水模型只有氧原子具有范德华作用系数，氢原子的范德华系数为0。

以上两种模型有对应的改进模型，SPC的改进模型为SPC/E，起主要改进其实就是使溶液系统的总能量乘以5.22 kJ/mol。这样可以使SPC溶液属性更加接近实验值。TIP3P在CHARMM力场中的改进是给氢原子一定的范德华系数，这样做的结果的计算根据复杂。。。（很无奈，因为结果好，所以也没有办法。）

由于真是情况下水分子的电量分布并不是完全在原子上的，如氧原子的一部分负电量就在H-O-H的对角线上，还有两个电子对处在H-O化学键的延长线上。为了得到更加真实的水分子模型，四个粒子以上的模型就被应用到分子动力学模拟中。其中最著名的有TIP4P模型。该模型在三个原子中间，H-O-H化学键的对角线上多了一个不含质量，只带电量的点。很多蛋白质模拟计算中，TIP4P和OPLS力场结合使用都得到很好的效果。

以上提到，水分子的氧原子在H-O化学键延长线上有两个电子对，于是有的人就在这两处添加了两个只带电量的粒子。2000年报道的TIP5P模型，计算结果也很好。还有一些牛人，结合TIP4P和TIP5P，要研制TIP6P，很好很强大。。。

不得不说，并不是模型的所含粒子越多越好。粒子越多，就算付出越大，因为要计算的相互作用更多 （4）力场

“力场”，请不要被“场”这个听起来像是十分高深的物理名词给吓坏了。 分子动力学模拟中使用的力场，包含两个重要的部分：

5

1）模拟粒子之间相互作用的方程（即经典力学的相互作用力方程如库仑定律，范德华作用方程等）。

2）方程的参数（即各个不同粒子，原子本身的参数，如带点量等等）。 可以想想，模计算机模拟好多成键或者不成键的粒子的运动，总要让它们互相推推拉拉吧，于是力场就是定义它们推推拉拉的方式（按照物理定律）。

力场类型，一般分类为三种：

i）全原子力场：精确定义每一个原子的参数。

ii）联合原子力场：省略非极性氢原子，同时把其参数整合到与他们成键的相邻原子上（比如甲基，只由一个碳原子表示）。

iii）粗颗粒力场：进一步精简分子结构的力场参数，种类比较多，比如有讲蛋白侧链看作一个颗粒的力场，或者甚至将整个氨基酸残基看成一个颗粒的力场等等。

一般来说力场的方程和参数是自成一个系统的，所以一般不能在一个系统中使用两个力场的参数。更具体的将，同一个原子在力场一中的带电量与起在第二个力场中是不一样的，化学键也一样。一般来讲，也不能特定修改力场中模一个原子的参数，因为原子之间是互相交叠依赖（比如未来保证整个氨基酸残基电量为0，各个原子电量加和必须为0）。但是，这并不是说一定不行，相反的，为了模拟一些不常见的分子，经常需要根据已有的参数（力场里面的，其他论文等）来构建新的分子参数。具体方法可以参考Mr. Google等著名老师。 目前比较流行的力场有：

AMBER： 包含好几个版本的力场，为全原子力场；

6

CHARMM：全原子力场，是软件CHARMM的一部分；

GROMOS：GROMOS软件使用的力场，版本较多，为联合原子力场；OPLS：包含全原子和联合原子力场两个版本；

粗颗粒力场：种类较多，没有固定版本或者种类，一般根据研究需要开发。 （5）Gromacs重启模拟计算

以前介绍过如果使用 GMX 3.x 重新由于种种原因停止的模拟，以下为 GMX 4.x 下重启模拟的方法。

GMX 4.x 的模拟程序 mdrun 较以往版本有不少不同。在模拟过程中，mdrun 按照 mdp 文件在一定时间间隔保存一个断点文件（checkpoint file, .cpt文件），该文件保存了该时刻模拟系统的所有物理量信息。如果由于不可预见原因，模拟中断，则可以使用该文件重新在该时刻开始进行模拟。

重启模拟的命令如下： -------

mdrun -s topol.tpr -cpi state.cpt -append -------

以上 state.cpt文件为最新生产的断点文件（ mdrun 会保存另外一个断点文件：state_prev.cpt，为上一个时刻保存的断点文件，双保险。）使用 “-append \的作用是将模拟输出添加到已有文件中，包括轨迹文件，记录文件，能量文件等，相同帧的信息将被后生产的信息覆盖。

当然，也可以继续像 GMX 3.x 一样使用 tpbconv生产新的 tpr 文件继续模拟，详细请参见旧文或手册。 （6）Gromacs多链模拟

进行模拟计算时，如果模拟分子由两条以上的链组成，一般都要明确告诉

7

模拟软件区分两条链。模拟软件一般没有那么聪明，除非明确定义，否则它会把两条以上的化学链（如肽链，DNA，其他聚酰胺等）看成一条链。在建立模拟文件是，上一条链尾端会于下一条链头部加一个共价化学键（如肽键）。由于该化学键一般很长，开始模拟时系统就“爆炸”了。 AMBER软件在处理这样的问题的，需要编辑原始的PDB文件，在每一条链结尾处添加“TER”。在GMX中，这种做法行不通（其实开发人员应该考虑这个问题）。解决的办法要在原始PDB文件中给每一条链添加链标识符，如“A”，“B”等等。（如果26个字母不够用，那就使用数字1到9，然后还可以使用特殊字符，如\，”￥“等等）。这样，使用 pdb2gmx 处理PDB文件的时候，就会得到各个链的拓扑文件，如topol_A.itp，topol_B.itp等等，并都被topol.top包含。

以上所述使用一个字符标识PDB文件中不同的链，是因为PDB文件只使用第22字符列作为链标识位，两个字符以上不认。（即AA，AB标识的链都被认为是A链。）

那么如果拿到一个没有链标识符的PDB坐标文件或GRO文件，该怎么办呢？那么要先使用make_ndx将不同链的残基选作不同的分子组（group），然后使用editconf将不同组输出成带链标识符的PDB文件，命令如： editconf -f File.pdb(/File.gro) -n indenx.ndx -o chian_A.pdb -lable A （以上可以等等A链，以此类推得到不同的链的PDB文件）。

最后将这些PDB文件组合成一个PDB文件，再由pdb2gmx处理即可。甚是麻烦。。。 （7）GMX.5 eneconv

GMX 分子模拟，有一个非常重要的能量输出文件，即 edr 文件。eneconv 就是对 GMX 能量输出文件进行处理的程序。

8

一 个模拟可以分对次进行，于是得到很多 edr 文件。使用 eneconv 的 “ -f ” 参数，然后把这些能量文件罗列出来，那么就可以对这些能量文件进行合并，输出一个完整的能量文件。如果另个能量文件中有重复的模拟步骤，那么后一个读入的 能量文件将覆盖前一个文件。也可以使用 “-settime” 参数对每一个输入文件的开始时间进行设置，以免互相覆盖。如下就是一个程序运行例子： ++++++++++++++++ eneconv -o fixed.edr -f *.edr ++++++++++++++++

即对当前目录下所有 edr 文件进行合并，然后输出为 fixed.edr文件。 当用 “ -f ” 参数读入单一一个能量 edr 文件时，也不是没有用，可以和其他参数配合，对能量文件进行编辑, 如 “ -dt ” 参数可以设定对原来能量文件进行规定时间间隔输出到新能量文件中；“ -offset ” 参数设定写出输入能量文件的时间帧（从哪一个模拟时间开始写入新的能量文件）等。这些参数，还有上面说到的 “-settime” 参数，都可以一起使用，加上 “-b” 和 “-e” 设定开始和结束读取模拟时间帧，就能得到新的称心如意的新能量文件。 程序输入文件： ----------- -f:

输入能量文件，即 edr 文件。 -o:

输出文件，也是 edr 文件。 其他参数：

9

--------- -b:

设定从哪一个模拟时间帧对输入文件进行读取。 -e:

设定从哪一个模拟时间帧对输入文件结束读取。 -dt:

设定输出文件的模拟时间间隔，比如 “-dt 10” 表示每 10ps 输出一次。 -offset:

设定从哪一个时间帧开始输出到新的能量文件中。 -settime:

交换式设定每一个输入文件在新输出文件中的开始时间。 -sort:

自动排序输入文件。 -scalefac:

该参数输入为一个实数，程序会将能量文件中的每一个能量项乘以这一个实数。 -error:

如果输入文件中有错误，程序自动退出。

注意：新的输出文件中，只有能量项是正确的，用于统计的 sigma 和 E^2 并没有更新，所有需要使用其他工具，如 g_analysis 进行新的统计

（8）GMX.4 editconf

editconf 是 GMX 最重要的程序之一。 它的主要功能是对系统结构进行编辑，同时它也把系统结构文件保存或者转换到不同的文件格式中，如 gro、g96、pdb 文件等。

10

在 分子动力学模拟中，通常给模拟系统添加一个模拟盒子（周期性的盒子，原子从这边出去了，就从那边进来，通俗吧）。 editconf 使用 “ -box ”、“ -d ” 和 “ -angle ” 等参数对模拟系统的盒子进行设定。在为系统设定盒子的时候，“ -box ” 和 “ -d ” 都把系统放置在盒子的中间，除非 editconf 明确使用另外一个参数 “ -noc ”。 (也就是 not center 的意思啦。)

使 用 “ -bt ” 参数，editconf 可以设定使用盒子的类型。editconf 支持以下几种盒子类型：triclinic （斜方体）、cubic（正方体）、dodecahedron (等边十二面体)， octahedron (等边八面体，两头被剪切的那种，即两个四面体放在一起，切去方向相反的尖尖的两头，同时保证所有的边长相等）。等边十二面体和等边八面体的体积分别是同“ 周期映像距离”的正方体体积的 0.71 和 0.77 倍，越接近球形，体积越小，计算代价越小。

当使用立方体、等边十二面体或者等边八面体时，“ -box ” 的参数可以只为一个实数值，该数值即为边长。“ -box ”参数也可以是三个实数，结果为以这三个实数为边长的长方体（如不指定 “ -angle ”，默认盒子矢量间夹角为90度）。

参数 “ -angle ” 只同参数 “ -box ” 和参数 “ -box ” 一起使用，以指定盒子矢量间的夹角；不能和参数 “ -d ”（用于指定系统原子到盒子边界的最小距离）一起使用。 如果使用 “ -n ” 或者 “ -ndef ” 参数为程序指定一个索引文件，则 editconf 可以选择系统中某一个组计算盒子大小或者几何中心，否则整个系统都会被考虑在内。 参数 “ -rotate” 用来旋转系统，按照给定数值在该坐标系上进行旋转。如 “ -rotate 0 30 0 ”，则表示将系统绕 Y 轴顺时间方向旋转 30 度。

参数 “ -princ ” 用来将系统（或者系统某一部分）固有坐标对齐到坐标轴上。这样做可以缩小盒子的体积，特别是当分子为一长条形状的时候。但是分子在模拟过程中会平移或者转动，所以使用的时候要特别小心。

扫 描参数 “ -scale ” 在其他参数被读入前可以扫描系统中各个原子的坐标，与其他

11

参数配合使用可以对系统中原子坐标进行修改以得到不同的系统性质。如果 “ -density ” 一起使用，可以得到不同的系统密度（这样做会改变系统盒子大小，要特别小心！）。使用该参数时，如果输入为 .gro 文件，可能结果会不精确。 当 “ -scale ” 参数为单一 “ -1 ” 时，输出的系统结果为该方向上的镜面映像，当三个方向上的输入全部为 “ -1 ”（即 “-scale -1 -1 -1”）时，系统的结果机构为原来结构的坐标原点对称映像。

在程序输出时，可以只输出系统的某一个组，或者某一个部分。可以建立划分更为细致的索引文件，这样可以进行细致的选择。

系统结构的周期性可以在程序中进行粗略去除，但是必须保证结构文件中周期性盒子的信息绝对正确，因为 editconf 去除周期性的算法十分简单，只是将原子坐标直接减去盒子边长。

当 程序输出文件为 .pdb 文件时，可以使用 “ -bf ” 参数为输出文件添加 B-factor。B-factor 要在一定的格式的文件中读取，这个格式要求文件第一行为文件中所含 B-factor 数值个数，文件第二行开始没行要有一个索引号，然后就是该索引后面的 B-factor 值。 B-factor 默认是一个残基一个数值，如果每一个原子有一直 B-factor 值，则要使用 “ -atom ” 参数。 如果使用 “ -legend ” 参数，那么程序在结构文件中会生成一列 CA 原子，这些原子带有 B-factor数字，并从最小 B-factor 值排列到最高，便于可视化软件显示。

如果使用 “ -mead ” 参数，那么程序将输出一个被静电统计软件 MEAD（玻松-波尔兹曼方程求解软件）使用的 pdb 文件（或者 pqr 文件）。使用这个参数有一个前提条件，就是输入文件必须是模拟打包文件（如 tpr 文件），因为这样的文件包含了力场参数。这个输出 pdb (pqr) 文件中，B-factor 数值列为该行原子的范德华半径，而 occupancy 列则为该原子的电量。“ -grasp ” 参数的作用很类似，但是 B-factor 列和 occupancy 列的数值互换。

12

另外一个十分的有用的参数就是 “ -label ” 参数，它可以为 pdb 文件加上一个链标记。如果一个文件里面不同残基属于不同肽链，那么这个参数可以帮这些残基指定肽链归属，这样不但可以帮助可视化，在建立模拟系统时也十分方便。

editconf 还可以修改结构文件的盒子类型，比如以下就是立方体盒子修改成一个等边八面体:

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

editconf -f -rotate 0 45 35.264 -bt octahedron -box -o +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

其中， 的数值为正方形变成乘于根号三除以二: a*sqrt(3)/2 。 程序输入文件： ------------- -f:

输入结构文件，文件格式可以是 gro、g96、pdb、tpr、tpb 或者 tpa。 -n: 索引文件。 -o:

输出文件，文件格式为 gro、g96 或者 pdb 。 -mead:

输出文件，MEAD程序的坐标文件，包含力场参数。 -bf:

B-factor数值文件，文件格式见上文。 其他参数：

13

--------- -w:

程序结束自动打开输出文件。 -ndef:

在索引文件的默认组中选择输出范围。 -bt:

指定盒子类型，支持的盒子类型有 triclinic（斜方体）、cubic（正方体）、dodecahedron（等边十二面体）或者 octahedron（等边八面体）。 -box:

指定盒子三个矢径，（a，b，c）。 -angle:

盒子矢径间夹角，默认都是 90 度（bc，ac，ab）。 -d:

溶质分子（蛋白，核酸等）到盒子边界的距离。 -c:

把溶质分子放在盒子中间，“-box” 和 “-d” 参数的默认参数。 -center:

为输出结构指定几何坐标中心，默认为 （0，0，0）。 -translate:

平移分子机构，比如 “-translate 0 0 2.5” 表示将分子结构向 Z 轴正方向平移 2.5 纳米。 -rotate:

旋转分子结构，见上文。

14

-princ:

对齐分子结构，见上文。 -scale:

扫描分子结构，“-scale 1 1 1” 表示在三个坐标方向上进行扫描。 -density:

指定输出结构文件的分子密度。 -pbc:

处理结构文件的周期性。 -grasp:

将原子的力场参数存储在输出文件中，见上文。 -rvdw:

给定默认原子的默认范德华半径。该参数可用于处理缺少力场参数的原子。 -sig56:

使用范德华能量最低点到原子中心距离的一半距离，而不使用简单使用范德华半径的一半。（特殊使用，多特殊？） -vdwread:

从默认 vdwradii.dat 文件中读取原子的范德华半径，而不使用力场提供的范德华半径。 -atom:

为各个原子添加 B-factor 数值。 -legend:

在结构文件中添加 B-factor 参考值，见上文。 -label:

15

为各个残基指定肽链归属。

注意：对于复杂的分子，最好使用 trjconv 去除结构周期性。 ：） (9)朗格文动力学 ( Langevin Dynamics ) 2009-06-19 14:10

朗格文动力学 ( Langevin Dynamics ) 是控制模拟系统能量的一种常用算法，在多种分子模拟软件中都可以看到。

分子模拟在一定的系综下进行，所以要保持系统状态不变，比如控制系统温度，压强等。由于计算机不是百分百精准，微小的误差在长时间的模拟过程中可能被不停积累和放大，于是需要不同的方法对系统进行不停调整。这些调整的方法有很多，比如 Berendersen, Nose-hoover, Langevin。

模拟系统中的原子受到起周围原子相互作用力，依据所受作用力根据牛顿第二定律运动。朗格文动力学的实现方法是给原子添加两个额外的作用力，即摩擦力和随机力。该摩擦力大小为原子速度乘于其质量再乘于一个摩擦因子（gamma)，其方向与原子速度相反。而随即力可以理解为来自溶液分子的随机相互作用等。这个两个力一起调节系统中各个原子的运动，以达到对整个系统能量的调控，即调控系统温度，压强等。

(10)GMX.3 do_dssp

本程序调用第三发软件 DSSP，从轨迹文件中读出蛋白质的二级结构信息。要使用该程序，必先安装好 DSSP 软件（恩，用 google 或者 baidu 小搜一下）。安装 DSSP 完毕后(/home/user/soft/bin/dssp)，在 shell 中定义 DSSP 变量：

setenv DSSP /home/user/soft/bin/dssp

如果使用 bash，则 export DSSP='/home/user/soft/bin/dssp'。可以该变量直接加到 shell 的配置文件中。

16

程序输出文件为 .xpm 文件，该文件其实是一个文本文件，可以一般文本编辑器打开。文件中用不同字符表示蛋白质每一个残基属于什么二级结构，并随时间变化，同时定义每个字符的颜色。.xpm 文件可视化有一点麻烦，Gromacs 的另外一个程序 xpm2ps 可以将其转化为 PostScript 文件，为 .eps 文件后缀，可以直接打开或直接使用到 Latex 文件中。

do_dssp可以统计属于某二级结构类型的残基数目，“-sss” 参数用于选择二级结构类型，“-sc“ 参数则把统计结果输出到 scount.xvg 文件中（其实这个文件包含了所有不同二级结构氨基酸残基数目，可以用 xmgrace 的 “-nxy” 参数打开）。本程序还可以计算每一个残基的溶液可及化面积(Solvent Accesible Surface, SAS)，以绝对平方埃（A^2)表示，或者最大 SAS 面积的分数表示（最大 SAS 定义为该一残基在甘氨酸链中的 SAS 面积）。Gromacs 另外一个程序 g_sas 也可以求取 SAS， 效率较高并且相对简单。

程序同时可以使用输出一个 ssdump.dat 文件。该文件包含了各个残基的二级结构信息，可以被另外一个程序 g_chi 使用，以分析残基直接二面角性质与二级结构直接的关系。（其实这个文件就是一个文本文件，里面用字符代表残基的二级结构，如 H 表示螺旋，B 表示折叠等。） 程序输入输出文件： --------------- -f:

轨迹文件，xtc、trr等。 -s :

模拟系统 tpr 文件。 -n: 索引文件。

17

-ssdump:

输出 ssdump.dat 文件。 -map:

程序输出 xmp 文件的原色库文件，如无则默认输出。（让其默认挺好的啦） -o:

输出文件，xmp 文件。各个残基属于某二级结构信息并随模拟时间变化。 -sc:

统计某二级结构残基数目，输出文件 scount.xvg。 -a:

各个残基的溶液可及化面积，xmp文件。 -ta:

总可及化面积，包括疏水和亲水面积，xvg 文件。 -aa：

平均可及化面积，xvg 文件。 其他参数： -------- -b:

指定可是统计时间帧。 -e:

指定结束统计时间帧。（“-b” 和 “-e” 可以用来读取轨迹中某一时间断的信息。）-dt: 指定读取帧的时间间隔。 -tu:

设定输出时间单位。

18

-w:

统计结束，自动打开输出文件。 -xvgr:

在输出 xvg 文件中添加其他信息，坐标标题等。 -sss:

设定统计二级结构类型。 (11)GMX.2 anadock

anadock是一个分析分子对接（docking）软件 Autodock 计算结果的程序，其作用是根据距离或者 RMSD 对分子对接结果的分子结构进行归簇，然后依据分子对接能量和自由能，并进行统计和输出。（本人没有用过 AutoDock，只用过 Dock，故不作任何评论。）

对分子结构进行归簇，Gromacs 也提供 g_cluster 程序，故可以使用该程序对结构进行提前归簇在计算能量，同样可以达到目的。 程序输入输出文件： ----------------------- -f:

输入文件，PDB 文件格式。 -ox:

归簇输出文件，PDB文件格式。 -od:

能量输出文件，为 xmgrace 的 xvg 文件格式。 -of:

自由能输出文件，xvg 文件格式。 -g:

19

程序记录文件。 其他参数： ------------ -xvgr:

在 xvg 输出文件中添加文字说明，坐标名称等。 -free:

使用 AutoDock 估算的自由能对结构进行结构排序。 -rms:

使用 RMS 进行归簇。 -cutoff:

归簇尺度，偏离大于给出数值时认为不同簇，默认为 0.2 nm。 （13）Restraint & Constraint

在分子模拟中，这两个词都表示限制，但是意义不同。通俗点讲，受到 Restraint 的分子或原子原则上可以运动；受到 Constraint 的分子原子原则上不能运动。

Restraint 通常的实现方法是给目标施加一个弹性势能，当目标偏离平衡位置时，弹性势能使目标恢复到平衡位置。Constraint 则直接使目标固定住，则施加无限大的弹性势能。

在 MD 软件如Gromacs中，Restraint 定义的弹性势能直接施加到指定原子上，并在 pores.itp 文件中列出。当 MDP （Molecular Dynamics Parameters) 文件使用关键字 “define = -DPOSRES” 时，引用该 pores.itp 文件。 Constraint 的实现则在MPD文件中定义 ”freezegrps “ 关键字，把需要进行 Constraint 的分子或原子直接添加到关键字中。

BTW: GMX 中 Restraint 和 Constraint 都具有方向性，即限制可以在某一个方向上

20

实现。其他模拟软件雷同。 （14）GMX.1 参数概述

根据GMX手册预告，所有GMX的程序都有6个标准参数，其中有的是隐藏的（参数太多，吓死人）。此六个参数及其作用罗列如下： -h:

在屏幕上打印该GMX程序的帮助信息，并退出程序。就是不使用这个参数，任何GMX程序运行时，都会打印很多信息，包括GMX的软件消息，简短的协议，还有各个参数的输入等。 -hidden:

打印隐藏的命令参数。 -quiet:

打印尽可能少的信息，所有参数输出被忽略，所以很quiet。 -man:

以指定文件格式，打印程序的手册到指定文件中。支持的文件格式有html, tex, nroff, acsii, completion, py 和 xml. （真怀疑开发人员怎么不被这么多文件格式给烦死，牛人。） -debug:

输出程序的调试信息。这个对终端用户不是很有用，输出文件一般为core加一串数字。 -nice:

设置程序的优先级。一般忽略，相信聪明可爱的计算机。 以下为GMX一些其他介绍：

1、GMX编译过程中，如果使用Motif或者Lesstif图形包，那么可以使用 -X 参数得到一个编辑各个参数输入的图形界面。如果没有使用以上图形包，使用该参数会

21

报错。

2、在SGI-IRIX系统上编译的GMX，可以使用 -npri 参数设定无阻塞（non-blocking）优先级。

3、可选参数的如无指定文件，那么就不会有输出文件（不然垃圾文件太多）。但是必选参数的输出文件一定会输出，如果指定按照默认文件明输出。

4、所有的GMX程序都接受无文件后缀名的文件，并自行补全文件名。如果不指明输入文件命，程序会自动在当前目录下搜索默认文件名。如果有多个同名不同文件类型的输入，那么第一个支持类型文件会被使用。（相信很多人一般不至于懒惰到这个程度，连个文件名都不给人家程序。）

5、几乎所有GMX程序（不包括 mdrun、nmrun 和 eneconv）都会先检查输入参数是否合法（重复参数等），如出错程序自动退出。 6、枚举参数（enum）必须为指定类型之一。

7、矢量输入必须文一个或者三个数值，如果指定一个数值，那么其他两个与指定数值一致。

8、多数GMX程序都可以使用 -tu 参数来指定时间单位，如 ps、fs、ns等等。这个多数应用与分析输出，成图等。

9、所有GMX程序可以直接读取 g-zipped 压缩文件，xtc、trr 和 trj 等压缩轨迹文件可能会报错，建议只打包不压缩。

10、大多数GMX程序可以处理比输入文件中原子数少的轨迹文件（轨迹文件中原子数目少于输入文件），但是必须保重轨迹文件中的原子处于输入文件的开头（即轨迹文件中没有的原子在输入文件尾端）。（这个貌似很容易产生歧义，语文没学好，见谅！） 造成这个的原因是GMX系统中，输入文件中所有原子按照序号进行，只有一个数字，没有名称。 （15）进行模拟的基本步骤行

22

以下是进行分子动力学模拟的最基本步骤，来自Gromacs维基。具体的步骤可能会因为模拟目标的不同而存在差异，这里只做基本的概述。 十分清楚要什么样的模拟系统性质和现象，明白模拟的目标。

选择适当的工具以实现模拟目标。这需要熟知其他研究人员对相似模拟系统的模拟方法等资料，要好好读别人的论文，主要包括：

-模拟使用的软件，当然还要考虑使用的力场，有时因为力场关系，需要更改使用的软件。

-力场描述了系统中原子的性质和相互作用，选择一个适合自己研究目标的力场至关重要，还是需要多看论文。根据sensenbobo血泪史和别人的经验，模拟核酸的话，可以采用amber力场，水分子可以使用TIP3P；模拟蛋白质的话，可以采用OPLS-AA，水分子使用TIP4P（重申一下，我本人没有跑过测试）。

获取模拟分子的坐标文件（PDB），或者自己生成一个。自己生成很好啊，强烈鼓励，但是要十分明白有一定的危险性，要明白每一步模拟系统的现象是为什么。 根据模拟软件，生成模拟系统初步坐标文件。

获取或者生成系统的拓扑文件，比如Gromacs的pdb2gmx命令，或者使用网络服务器PRODRG（google一下这几个字母），amber的leap命令和antechamber，NAMD的psfgen或者VMD。请参考相关手册。

给模拟系统添加一个盒子（盒子的翻译方法我十分不喜欢，箱子显得太大，边界显得太难理解，但是盒子听起来想到蛋糕，我们是严肃的科学家。娃哈哈。）在盒子中添加水分子，添加离子如NA和CL之类（添油加醋，生活就是这样进行的）。添加任何东西到你的系统中，就要相应的修改一下拓扑文件，以保持与系统对应，不然会出错哦。

进行能量最优话。这样做是消去系统中原子之间距离太近造成的高能量。必要的时候，分子可以现在真空中能量最优话，然后再做溶剂环境中能量最优话。

23

选择恰当的模拟参数模拟平衡模拟系统。需要慎重考虑力场的作用，必要的时候保持系统体积温度不变（NVT系综）进行位置限制性平衡系统，然后在保持系统温度压强不变（NPT系综）修正系统的密度。然后再选择正确的系综（NPT，NVT，NET，BT？？）进行数据收集。（真烦！）

跑足够长时间的模拟，使系统达到理想的平衡。使用模拟软件的工具或者其他软件观察系统的性质，如温度，能量瞪瞪。

使用适当的参数进行数据采集模拟，注意模拟衔接过程中各个系统性质保持不变，如果原子速度之类。这里还要考虑力场的作用，同时要考虑怎么测量描述模拟目标。 模拟足够长的时间使模拟目标明显出现（可能也出现跟你期望相反的结果，出现率应该百分之五十，加上个人感情因素，可能会达到百分之八十。保重！）。 分析模拟结果，解释得到结果。要好好利用视图工具，同时要十分注意漂亮的三维图说明不了什么问题。图表，比较，统计，一个都不能少。 就这样吧。

（16）Gromacs使用－新手入门篇】

其实gromacs的思路比较简单，大体就是以下的流程,在此先对em进行一下总结：

1. convert the pdb-file to a gromacs structure file(.gro / .pdb) and a gromacs topology file(.top) 这一步需要gromacs的预处理程序pdb2gmx 标准的输入命令应该是：

＃ pdb2gmx -ignh -f speptide.pdb -p speptide.top -o speptide.gro

这个命令之后有一个force field的选择 系统给出了11种力场，具体以后实际计算中如何选择力场还是应该多参阅文献和gromacs说明书的。 这步之后可以在屏幕看到：

--------- PLEASE NOTE ------------

You have succesfully generated a topology from: speptide.pdb. The G43a1 force field and the spc water model are used.

Note that the default mechanism for selecting a force fields has changed, starting from GROMACS version 3.2.0

--------- ETON ESAELP ------------至此，完成了从pdb到gro和top的转变。另外，我也看过有的作这一步命令时候并不进行top文件的转换，因为在genbox（对盒子加水溶剂化时还要对top文件进行修改）

2. solvate the peptide in water 这一步需要用到editconf和genbox两个程序，命令的标准书写应该是： # editconf -bt cubic -f speptide.gro -o speptide_box.gro -c -d 0.5 可以在屏幕上看到 Read 191 atoms

Volume: 8.17315 nm^3, corresponds to roughly 3600 electrons

24

No velocities found

system size : 1.774 3.372

1.367 (nm)

diameter : 3.517 (nm)

center : 2.650 1.453

2.417 (nm)

box vectors : 1.774 3.372

1.366 (nm)

box angles : 90.00 90.00

90.00 (degrees)

box volume : 8.17 (nm^3) shift : -0.391

0.806 -0.158 (nm) new center : 2.259 2.259

2.259 (nm)

new box vectors : 4.517 4.517

4.517 (nm) new box angles :

25

90.00 90.00

90.00 (degrees) new box volume : 92.18 (nm^3)

因为程序的默认，所以这一步命令可以写作： # editconf –f speptide –o –d 0.5

这里－f之后的sp文件是默认gro文件或者pdb文件，－o后面生成的output file是放入盒子的体系的gro文件 按照这样的格式，输出的文件被默认命名为－out.gro 体系放入盒子之后在盒子中注入水

# genbox -cp speptide_box.gro -cs -p speptide.top -o speptide_water.gro 也可简写为 ：

# genbox –cp out –cs –p sprptide –o b4em Screen--

Output configuration contains 9035 atoms in 2967 residues Volume :

92.1813 (nm^3) Density :

994.449 (g/l)

Number of SOL molecules: 2948

之后看加水之后的top文件 用命令 # tail speptide.top 可以看到 [ system ] ; Name

Protein in water

[ molecules ] ; Compound #mols Protein 1 SOL 2948

其中SOL后面的数字2948就是genbox加入盒子的水分子的数目。

Editconf程序的另一个用途是讲gro文件转化回pdb这时可以讲speptide_water.gro转化回pdb观察 ＃ editconf –f speptide_water.gro –o speptide.pdb 拖回本机 用spbdv或者vmd观察 。

The next step is to generate index file . in the tutor of

gromacs , we are told that there are a set of index groups to select, unfortunately, I didi not find them ,so I have to use make_ndx to generate one . # make_ndx –f b4em

3. perform an enery minimization of the peptide in solvent

26

Now the simulation system is almost ready. Before we can start the dynamics, we must perform an energy minimization, to alleviate any bad contacts (atoms overlapping such that a significant repulsion would result, causing numerical problems in the simulation) that might be present in the system.

# grompp –v –f em.mdp –c speptide_water.gro –p speptide.top –o speptide_em.tpr Or use this command for short:

# grompp –v –f em –c b4em –p speptide –o em

After this command ,the bad contacs have been removed .so we can do enery minimization now. # mdrun –v –s speptide_em.tpr –o speptide_em.trr –c after_em.gro –g emlog.log Or

# mdrun –v –s em –o em –c after_em –g emlog From screen ,we can see these ------

Steepest Descents:

Tolerance (Fmax) =

2.00000e+03 Number of steps = 100 Step=

0, Dmax= 1.0e-02 nm, Epot= -8.02562e+04 Fmax= 2.70468e+04, atom= 6792 Step=

1, Dmax= 1.0e-02 nm, Epot= -8.69641e+04 Fmax= 1.24531e+04, atom= 6657 Step=

2, Dmax= 1.2e-02 nm, Epot= -9.48067e+04 Fmax= 5.54072e+03, atom= 9009 Step=

3, Dmax= 1.4e-02 nm, Epot= -1.02743e+05 Fmax= 2.86957e+03, atom= 8766 Step=

4, Dmax= 1.7e-02 nm, Epot= -1.08941e+05 Fmax= 1.53075e+04, atom= 150 Step=

5, Dmax= 2.1e-02 nm, Epot= -1.10083e+05 Fmax= 1.31653e+04, atom= 150 Step=

6, Dmax= 2.5e-02 nm, Epot= -1.10564e+05 Fmax= 2.16454e+04, atom= 150 Step=

7, Dmax= 3.0e-02 nm, Epot= -1.11528e+05 Fmax= 1.82398e+04, atom= 150 Step=

9, Dmax= 1.8e-02 nm, Epot= -1.12767e+05 Fmax= 3.87548e+03, atom= 120 Step=

11, Dmax= 1.1e-02 nm, Epot= -1.13598e+05 Fmax= 1.23013e+04, atom= 120 Step=

12, Dmax= 1.3e-02 nm, Epot= -1.14464e+05 Fmax= 5.37656e+03, atom= 120 Step=

13, Dmax= 1.5e-02 nm, Epot= -1.14787e+05 Fmax= 1.74494e+04, atom= 120 Step=

14, Dmax= 1.9e-02 nm, Epot= -1.15874e+05 Fmax= 8.35820e+03, atom= 120 Step=

16, Dmax= 1.1e-02 nm, Epot= -1.16420e+05 Fmax= 6.65421e+03, atom= 120

27

Step=

17, Dmax= 1.3e-02 nm, Epot= -1.16807e+05 Fmax= 1.08106e+04, atom= 120 Step=

18, Dmax= 1.6e-02 nm, Epot= -1.17248e+05 Fmax= 1.12834e+04, atom= 120 Step=

19, Dmax= 1.9e-02 nm, Epot= -1.17425e+05 Fmax= 1.45145e+04, atom= 120 Step=

20, Dmax= 2.3e-02 nm, Epot= -1.17590e+05 Fmax= 1.67089e+04, atom= 120 Step=

22, Dmax= 1.4e-02 nm, Epot= -1.18845e+05 Fmax= 2.65542e+03, atom= 121 Step=

23, Dmax= 1.7e-02 nm, Epot= -1.19198e+05 Fmax= 2.39072e+04, atom= 120 Step=

24, Dmax= 2.0e-02 nm, Epot= -1.20588e+05 Fmax= 6.17310e+03, atom= 120 Step=

26, Dmax= 1.2e-02 nm, Epot= -1.20901e+05 Fmax= 1.04224e+04, atom= 120 Step=

27, Dmax= 1.4e-02 nm, Epot= -1.21229e+05 Fmax= 8.89909e+03, atom= 120 Step=

28, Dmax= 1.7e-02 nm, Epot= -1.21346e+05 Fmax= 1.51125e+04, atom= 120 Step=

29, Dmax= 2.1e-02 nm, Epot= -1.21686e+05 Fmax= 1.29390e+04, atom= 120 Step=

31, Dmax= 1.2e-02 nm, Epot= -1.22216e+05 Fmax= 3.21419e+03, atom= 120 Step=

33, Dmax= 7.5e-03 nm, Epot= -1.22523e+05 Fmax= 6.52631e+03, atom= 120 Step=

34, Dmax= 8.9e-03 nm, Epot= -1.22796e+05 Fmax= 6.31916e+03, atom= 120 Step=

35, Dmax= 1.1e-02 nm, Epot= -1.22998e+05 Fmax= 8.04922e+03, atom= 120 Step=

36, Dmax= 1.3e-02 nm, Epot= -1.23167e+05 Fmax= 9.93619e+03, atom= 120 Step=

37, Dmax= 1.5e-02 nm, Epot= -1.23317e+05 Fmax= 1.10029e+04, atom= 120 Step=

39, Dmax= 9.3e-03 nm, Epot= -1.23857e+05 Fmax= 1.88770e+03, atom= 810 writing lowest energy coordinates.

Steepest Descents converged to Fmax < 2000 in 40 steps Potential Energy = -1.2385705e+05 Maximum force =

1.8876959e+03 on atom 81 Norm of force =

1.3663420e+04 NOTICE!!!!!!!!!

If the potential enery after minimization is lower than -1.1e+05kJ/mol, it is acceptable and the structure can be used

28

for MD caculations

（17）gromac的一些笔记

用了一个星期，基本操作差不多了。套路自顶向下大概是这样：

1．核心步骤当然是mdrun了，当然后续的数据分析才是对理论水平的要求，先不讨论。这里需要一个tpr输入文件——由grompp生成的计算体系和计算控制详细信息的二进制文件。计算生成trr轨迹文件、xtc压缩的轨道文件、edr能量等状态参量记录文件、log计算记录等。后续处理靠这些输出文件和 gmx提供的命令。 2．grompp预处理文件，将三种文件整合到一个tpr文件作为mdrun的输入。分别是：

a．mdp文件：计算控制文件，里面要给md模拟时的计算参数加以描述，控制时间、精度、温度什么的。

b．gro文件：分子坐标文件，记录要模拟的体系中全部（要计算的）原子的坐标、名称等。 c．top文件：力场参数文件。记录gro文件中的原子之间相互作用立场参数和原子分子类别信息的文件。3．mdp文件一般可以看手册根据需要自己填写，或者看模版。

4．gro、top文件，用editconf加周期边界元胞，genbox加溶剂。那么最初的gro、top文件如何得到呢？ a．如果是从网站上直接当下来的pdb文件，用pdb2gmx转化。（有时需要先去掉里面的水，有时需要加-ignh不考虑多余的氢，有时需要加一些离子（genion）来平衡蛋白质上的电荷）

b．要是还要加入自己的分子，比如在MS什么的软件中先画了一个，存成pdb文件，这时不能直接转，因为残基gmx不认识，立场库里面没有描述，需要自己写top文件。不过这很累，还要自己计算，又不会，怎么办呢，网上推荐用ProDrg软件（只找到网络版的网址google一下就有），把pdb文件 copy进去，run一下就出来gro，和itp文件了（还有其他的一坨，用不到）。把自己gro文件里的内容copy到转出来的gro里，总原子数和序号改一下，再把itp文件include到老top文件中，[ molecules ] 标签下面添上加入的分子名和数，就可以grompp了。

具体的命令用法和文件格式就看手册吧，写得是相当的详细。网上相关的介绍例子也很多很清楚，尤其是那个氩气的。

在gromacs中使用Amber力场

由于种种原因，Gromacs中没有直接包含Amber力场。

参考 http://folding.stanford.edu/ffamber/ 的介绍，下面给出在gromacs中使用Amber力场的方法： Install FF:

1.下载力场 ：

http://folding.stanford.edu/ffamber/ffamber_v3.3.1.tar.gz (490 kB)

http://folding.stanford.edu/ffamber/ffamber_v3.3.1-doc.tar.gz (6.54 MB)带文献例子 2.解压缩

3.把aminoacides.dat、vdwradii.dat拷到top文件夹下覆盖，实际上一股脑考过去就完了。 4.子文件夹下的东西都拷到外面，也就是top里面。 5.修改FF.dat,修改第一行的数字，加上新添的力场数。 6.为每个新的力场增加新行，如：

ffamber94 AMBER94 Cornell protein/nucleic forcefield ffamber99 AMBER99 Wang protein/nucleic acid forcefield ffamber99p AMBER99p protein/nucleic forcefield ffamber03 AMBER03 Duan protein/nucleic forcefield Edit PDB:

1.DNA残基名改成DA、DT、DC、DG，分子为A、B。 2.所有的(*)字符换成(')。

3.端基的残基名改为DX5/DX3,或NXXX/CXXX。 4.DT中：O2->O 5.DT中：C5M->C7 6.DT中：H5M->H7

29

7.DC中：O2->O

8.比较麻烦:两个等位的H，形如nHX,->HXn 9.蛋白质的话，LYS->LYP

残基中的原子修改方法是将标准的DNA文件和gromacs里Amber的rtp文件逐一比较得到的，希望会对这方面的研究者有所帮助。

（18）gromacs读取碳纳米管

有人问我leap怎么读入碳纳米管，这里说一下，顺便也说一下gromacs读入的方法。只说操作方法，不涉及选用力场等问题。

碳纳米管结构比较大，而结构很有规律，一般仅有C原子。若不考虑边界问题，可以将全部电荷当成0，这种近似并不会对结果真实性有多少影响。对于碳纳米管这样的特征，适合当成一般的小分子来处理，但不宜用antechamber而且也没必要用。

先用Nanotube modeler生成一段碳纳米管，含240个C原子，导出结构为tube.pdb。但是会发现，原子名都是相同的，而leap要求结构中相同残基中每个原子都有独立的名字，所以需要先改名。比较方便的方法是用ultraedit，打开后开启列模式，选定从第1个原子到最后一个原子的15和16列(即原子名C后面的两列)，选Column-Insert Number，OK，即把原子名改为了C1至C240，保存为tubename.pdb。 这里假定C原子使用gaff力场的c2原子类型对应的力场参数，依次运行 xleap -f leaprc.gaff

a=loadpdb /sob/tubename.pdb

bondbydistance a 默认是2埃内成键，故碳纳米管的碳原子都会正常成键，可用edit a 检查。

然后把碳纳米管内全部240个原子都设成c2原子类型，可以用批处理脚本。由于leap只能一条一条运行批处理文件中的内容而不支持shell脚本的循环，需要将循环转化成普通命令脚本再在leap运行。写脚本： for ((i=1;i<=240;i=i+1)) do

echo \done

在shell下运行，得到leapdo，然后在leap里运行source leapdo即可

check a 进行参数检查，应该OK。如果想改力场参数或缺少某些参数，可以edit gaff往里添加和修改。 如果要加溶剂和普通情况无异，这里略过。 最后saveamberparm a tube.top tube.inpcrd =============gromacs读入碳纳米管

还是用上面的tube.pdb为例子。小分子在gromacs中常用prodrg来处理，但是理由和antechamber一样并不适合，而且有更方便的办法。gromacs提供了一个构建拓扑文件的工具x2top，专适用于构建结构规律性很强的体系。x2top中有一些bug和“规矩”，而且在不同gromacs版本中bug和“规矩”还不一样，这里使用gromacs4.0.4的x2top。

首先需要写n2t文件，比如使用ffgmx2力场，就在力场文件所在文件夹里写一个ffgmx2.n2t，内容如下： C CX 0.0 12.011 3 C 0.14 C 0.14 C 0.14 C CX 0.0 12.011 2 C 0.14 C 0.14 C CX 0.0 12.011 1 C 0.14 第一行说明如果体系中任何一个C原子，与周围3个C原子的距离都在0.14nm左右（判断标准大概为正负10%），就把它当作CX原子类型（原子类型名字自定，可以不属于力场包含的原子类型），电荷为0.0，相对原子质量为12.011。第二行和第三行与第一行意义类似，即代表与两个C原子和与一个C原子成键的C原子被当成什么类型、多少电荷、多少原子质量，此例中将碳纳米管所有原子都当成CX。 运行x2top -f tube.pdb -o tube.top -ff gmx2

这样就得到tube.top，里面包含这个碳纳米管全部键结项。

-ff代表用什么力场，此处即ffgmx2，如果-ff select则是出现列表自行选择。这里如果使用ffgmx力场即便正确写了ffgmx.n2t也可能无法正确执行，即不能正确根据距离判断成键，此时应换用别的力场来执行x2top。 如果运行后程序卡住不动，可尝试加-nopbc解决，这是一个bug。

默认情况下，x2top会自动在.top里相应键结项加上力场参数，这种方式自动加入的力场参数并不是根据力场中

30

相应原子类型间的键结参数加入的。平衡距离/键角/二面角就是tube.pdb中的相应项的距离/键角/二面角，键的力常数/键角力常数/二面角力常数默认是400000/400/5，可以分别用-kb、-kt、-kp设定。如果不想让其自动加入，可以加上-noparam。

如果n2t中把C原子转换为ffgmx2力场中已有的类型，比如CB，那么此时这个.top搭配tube.pdb已经可以直接用于模拟了。

此例将C设为力场中没有的原子类型CX是为了假设需要加入专门的力场参数的情况，此时需要定义新的原子类型。首先加入VDW和1-4VDW参数，分别在ffgmx2nb.itp里的[atomtypes]和[nonbond_params]里面添加。 如果前面x2top时加了-noparam，还需要设定键结参数。分两种情况：

1 直接在.top里面每一个键结项后面加入相应力场参数，类似于x2top自动在.top里添加力场参数的方式，用ultraedit等软件的列模式批量加入效率最高。

2 不在.top里面加，而是在ff???bon.itp里面的[bondtypes] [angletypes] [dihedraltypes]里面加入相应的项，或者把这些内容单独写进一个.itp文件，之后include进.top里。这种方法比较省事、易于修改。 关于这两种参数定义方式的讨论，详见《gromacs决定键结参数的方式》 http://hi.http://www.wodefanwen.com//sobereva/blog/item/dcc7ef253cc72e6735a80f53.html 之后加盒子、填充溶剂等与一般体系无异。

（18）gromacs决定力场参数的方式 文/sobereva last updated: 2009-Apr-5 Q:

43a1力场中，有的键长、键角、二面角都是空白的，只有原子的信息，这种情况是不是constraint啊？ 比如这样：

51 91 2 52 53 2 gb_4 A:

[bonds]里面指定了原子间连接关系，不管后面力场参数是不是空的，都说明原子是相连的。此时，如果在mdp里面设了constraint=all-bonds，所有[bonds]里面的键就会被约束住(不是按照距离远近判断，而是看[bonds]里面怎么连接的)。比如51 92 2后面虽然没写力场参数，此时也会被约束住。如果写了力场参数，则力场参数视为无效，也按约束来处理，用gmxdump -s察看.tpr，看不到它们处于成键项(比如810 type=366 (G96BONDS) 796 797)，而会看到处于约束项(比如672 type=429 (CONSTR) 796 797)。也就是说constraint的优先级最高。

不同的力场决定成键参数的方式不同。对于设constraint=none的情况下有两种情况：

对于GROMOS96力场来说，用什么参数只取决于grompp时读入的所要研究体系的.top(或者被include的.itp)，里面写了哪些力场参数就发挥哪些效果，没有写的力场参数就当作没有(被设为0)。比如[bonds]里写51 91 2，虽然指定了连接关系，但没写力场参数，又没有设约束，等跑起来之后就可能会逐渐远离。 对于OPLSAA力场来说，哪些键用哪个力场参数是在grompp时根据原子类型来决定的(类似amber中leap程序的方式)，而不像GROMOS96那样在.rtp里面事先定义。OPLSAA力场在所要研究体系的.top(或者被include的.itp)里面只定义有键结项(本文中键结项泛指键长、键角、二面角)并不直接定义参数，例如： [ bonds ]

1 2 1 ...... [ pairs ]

1 9 1 ......

[ angles ]

2 1 3 1 ......

[ dihedrals ]

2 1 4 5 3 ......

这和GROMOS96力场处理方式截然不同。虽然[bonds]里1 2 1没定义参数，但是用gmxdump -s看tpr却发现连

31

接1和2原子的键存在，模拟过程中也保持一定距离。这说明对于OPLSAA力场，grompp从.top(或者被include的.itp)获知的仅仅是有哪些键结项，这些键结项涉及的原子会通过开头[atoms]段里面的定义将原子序号转化为原子的type(如opls_136)，再由ffoplsaanb.itp的定义转化为bond_type(如CT)，再从ffopsabon.itp里面找到对应参数，放进.tpr里面，在模拟过程中生效。

对于[bonds]里的内容处理方法简单来说就是，设了constraint=all-bonds，不管什么力场，[bonds]里写了哪些项就约束哪些。如果constraint=none，用GROMOS96力场时有项有参数的就起作用，有项没参数就等于没写。OPLSAA里面有项没参数也起作用。

(关于自定义constraints的情况可参考《解析gromacs的restraint、constraint和freeze》) 善用gmxdump -s察看.tpr是最可靠的方法。

在这里我将问题扩展一下，说一下gromacs决定键结参数的机制。实际上OPLSAA力场与GROMOS96力场都是遵循同样的方法来处理，但是如上所述，决定键结参数有很大区别，并非因为程序对它们处理机制的差异，而是力场参数文件的差异。下面都不考虑约束问题。

例如用某个力场，研究的体系的.top里面的一个普通的成键定义： [ bonds ]

1 2 1 ga_1 //注意第二个1代表函数类型，不是原子号 首先gmx会用C预处理程序(一般就是cpp)将ga_1还原为实际的参数，转换关系定义在ff????bon.itp里，在grompp时加上-pp，从得到的processed.top中就会看到ga_1被还原了。当然如果ga_1位置直接写的就是参数，就不必转换了。这样1和2之间的键参数就定下来了。 另一种情况，虽然定义了成键项，但没写参数： [ bonds ]

1 2 1

那么程序会首先将根据原子序号，在[ atoms ]里面查找到对应的原子类型，比如转换为CH2 N，然后按照原子类型到ff????bon.itp里面的[ bondtypes ]里面找原子类型相对应的参数，找到了就定下来了，如果没找到，在grompp的时候就会出Warning提示这个键的参数被设为0，相当于不存在。

GROMOS96与OPLSAA定义力场参数的差异就在于它们的ff????bon.itp的内容。

比如GROMOS96的ffG43a1bon.itp里面[ bondtypes ]几乎没有什么内容，只是二硫键和血红素用的几个键的定义。所以必须在[ bonds ]里面直接写出参数(这步由pdb2gmx完成，主要将.rtp里面的相应残基的键/角/二面角项和参数搬过来)，否则在ffG43a1bon.itp的[ bondtypes ]里也不会找到参数，那个成键项等于没有被定义。 而比如OPLSAA的ffoplsaabon.itp，[ bondtypes ]里面包含了全部OPLSAA力场的成键项，所以就不必在体系.top的[ bonds ]里面直接写参数，grompp时在这里一般都能找到对应的。在ffoplsaa.rtp里面定义的除了原子电荷外，主要就是残基的连接关系，也就是[ bonds ]，没有参数，在[ angles ]和[ dihedrals ]里面几乎没有内容或者根本没有。这种情况下pdb2gmx的主要作用除了把.rtp的对应内容直接搬到.top里面以外，另一个主要工作就是根据连接关系自动补全角和二面角项，所以.top里面看到的键/角/二面角项是全的，如上所述没有写参数。

从大体来说，实际上GROMOS96力场是有些特殊的，它直接在rtp中事先定义好了全部键结项和参数，top也包含全部键结项和参数，grompp调用的参数直接取自top。而OPLSAA则算是比较普通的情况，rtp只定义连接关系，在top中更进一步把全部的键结项表达出来，该用什么参数等grompp时再去从力场参数库中搜。amber的leap也是如此，载入结构载入的只是连接关系，各个键结项用的参数在saveamberparm的时候对应去搜。 决定angle、dihedral的参数与决定bond的方式相同。另外如果比如[ bonds ]里面已经写了参数，同时在[ bondtypes ]里面也能找到原子类型对应的参数，则优先使用在[ bonds ]里面直接写的。

这里再说一下对二硫键等特殊键的处理方式。pdb2gmx读取结构时，也载入specbond.dat，specbond.dat内如如下 6

CYS SG 1 CYS SG 1 0.2 CYS2 CYS2 ......

这就说明如果发现在同一个链上有两个CYS上的SG间的距离在0.2nm附近(注意不是仅仅指的在0.2nm以内)，就认为这两个SG成键，并且残基名由CYS变成CYS2。这时这个半胱氨酸的参数就用的.rtp里面CYS2残基的参数，加氢时也使用.hdb里面CYS2残基的加氢方式而非默认的CYSH，也就是硫上面不加氢。在.top里面可以看到出现了这两个硫原子之间的成键项以及相关的角、二面角项，但是没有参数(无论是GROMOS还是

32

OPLSAA)，这就说明要从比如[ bondtypes ]里面读对应的参数，这就是为什么比如ffG43a1bon.itp的[ bondtypes ]里面有S S 2 gb_33这样的二硫键参数定义。硫原子的原子名是SG但是原子类型是S，所以二硫键和这个参数项对应。对于血红素的一些特殊成键也是这么处理的。

也可以自行添加、修改specbond.dat的内容来处理类似的特殊情况。 这里也说说非键作用和1-4作用，下文指的非键作用不包括1-4作用。 非键作用：

对于VDW参数的确定，GROMOS96力场在ffG43a1nb.itp里面的[ nonbond_params ]中提供了各种各样原子类型的组合，计算i、j两个原子间的VDW作用时会根据原子类型对应获得C(ij)参数(包括C6和C12项)。若计算两个原子间VDW作用时，发现[ nonbond_params ]里面没有对应的，就通过组合方法确定，从[ atomtypes ]中得到原子i的C(i)与原子j的C(j)，由于ffG43a1.itp里面comb-rule设的是1，会根据(C(i)*C(j))^1/2公式来获得C(ij)。由于OPLSAA并没有[ nonbond_params ]段，故所有C(ij)项都是如上算出来的(OPLSAA的comb-rule是3，详见手册)。

对于静电作用，只依赖于原子电荷，没什么特殊的。

要注意，在拓扑文件中的[ moleculetype ]段中，nrexcl决定了相隔多少个键以内不计算非键相互作用，例如nrexcl=5时，与某个原子相隔5个及以内的原子间都不计算非键作用，苯分子若设nrexcl=5，则导致完全不计算原子间的非键作用。 1-4作用：

1-4作用项虽然看起来也属于非键作用，常被称为1-4非键作用，但gromacs对它的处理与非键作用完全不同。1-4作用包括库仑1-4和1-4VDW相互作用，哪些原子间有1-4相互作用都定义在体系.top里面[ pairs ]，里面每一项就是一个1-4作用项，如果没写就说明这两个原子间不拥有1-4相互作用（即便从成键关系上属于1-4关系的原子也不算）。如果两个原子包含在1-4作用项中，则不计算它们的非键作用，相当于从非键列表中排除了。nrexcl对1-4作用完全无效，根据nrexcl的设定，即便某两个原子间不计算非键作用，若存在于1-4作用项中，也会照常计算1-4作用。即nrexcl只管非键作用。 1-4作用项中的库仑1-4作用，就是照常计算静电作用然后乘以FudgeQQ(在比如ffoplsaa.itp的[ defaults ]中定义)。 1-4作用项中的1-4VDW作用，所用C(ij)参数与非键VDW中的不同，需要额外的参数，但是GROMOS96和OPLSAA的体系.top在[ pairs ]里默认都没写参数，这是因为：

GROMOS96力场，明确定义了各种原子类型组合的1-4VDW参数，全部储存在比如ffG43a1nb.itp里面的[ pairtypes ]当中，体系.top里面[ pairs ]当中定义的每一个1-4VDW项的参数都从这里对照原子类型自动提取。在ffG43a1.itp当中会看到gen-pairs设为了no，意思是如果读不到对应的1-4VDW参数，就出现warning并将参数用0值代替。

OPLSAA力场，各种原子类型之间的1-4VDW作用都是普通VDW作用乘以FudgeLJ因子0.5计算得到，也就是进行削弱(GROMOS96定义的1-4VDW参数也是被削弱的)，在ffoplsaa.itp里面可以看到gen-pairs被设为yes，意思是读不到的1-4VDW参数都通过上述方法计算生成。实际上ffoplsaanb.itp里面根本就没有[ pairtypes ]，所以如果不做人为修改，1-4VDW参数都是通过计算生成。

如果直接在[ pairs ]当中的项的后面给出参数，这些项就直接用给的参数当作1-4VDW的参数。如果直接给的参数为0，或者没直接给参数、在ff????nb.itp里的[ pairtypes ]也没有对应项、而且把gen-pairs设为了no，则参数值会默认为0，这两种情况下1-4VDW作用效果皆为0。但此时1-4静电力仍会照常计算，除非把FudgeQQ也设为0，则彻底不计算1-4作用了。

不计算1-4作用实际上有两种情况：一种是不写[ pairs ]的情况，1-4作用完全忽视掉，程序根本不算。另一种是写了[ pairs ]，程序也算了，但是因为上面一段所述的参数和设定原因，导致算出来的1-4作用为0，相当于没算，这种情况下仍然会花费计算时间，在mdrun的结尾会看到1-4作用计算的开销。

PS: 基于GROMOS87的ffgmx力场，决定VDW参数方式同GROMOS96，决定键结参数方式同OPLSAA(但OPLSAA多了一步原子名的转换过程，即type到bond_type)。 （19）gromacs基础知识

Gromacs的pdb2gmx命令使用gromacs做分子动力学模拟时，第一个要用到的命令一般都是pdb2gmx。这个命令吧pdb分子文件转化成gromacs独特的gro分子结构文件类型，同时产生分子拓扑文件。

Gromacs是典型的GPL软件，每一个命令都有很多命令参数。这对熟悉windows环境的人来说有一点烦，但是如果熟悉了Linux环境，也就慢慢喜欢啦。（建议多使用命令，就像VMD， Pymol，rasmol和Chimera等等分

33

子可视化软件，如果接合命令使用，功能都非常强大。另外一个比较bt的软件叫做WHATIF的，完全建立在bt的命令菜单上，心理承受能力不强者多半吐血而终。

使用“pdb2gmx -h”可以得到pdb2gmx的所有参数及简单说明（gromacs的任何命令都可以使用-h参数得到类似帮助）。pdb2gmx的参数很多，但是常用的只有以下几个： -----------------------------------------------------------------------

-f 指定你的坐标文件，可以是pdb、gro、tpr等等包含有分子坐标的文件;

-o 输出文件，也就是处理过的分子坐标文件，同样可以是pdb、gro、g96等文件类型; -p 输出拓扑文件。pdb2gmx读入力场文件，根据坐标文件建立分子系统的拓扑;

-water 指定使用的水模型，使用pdb2gmx的时候最好加这个参数，不然后面会吃苦头。它会提前在拓扑文件中添加水分子模型文件;

-ff 指定力场文件（下文讨论），也可以不用这个参数，再自行选择;

-ignh 舍弃分子文件中的H原子，因为H原子命名规则多，有的力场不认; -his 独个指定HIS残基的质子化位置。

------------------------------------------------------------------------

其他的参数还不少，可以好好看一个pdb2gmx的帮助文件，一般的pdb2gmx的执行格式如下（假设你的分子坐标文件为sen.pdb）：

pdb2gmx -f sen.pdb -o sen.gro -p sen.top -water tip4p -ignh -his

该命令读入分子文件，使用tipp水模型，等等，然后pdb2gmx会让你选择力场文件。然后它就很聪明的帮你建立初始模拟系统啦。

gromacs自带的力场有很多：

------------------------------------------------------------------------ 0: GROMOS96 43a1 force field

1: GROMOS96 43b1 vacuum force field

2: GROMOS96 43a2 force field (improved alkane dihedrals) 3: GROMOS96 45a3 force field (Schuler JCC 2001 22 1205) 4: GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656) 5: GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656) 6: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals) 7: [DEPRECATED] Gromacs force field (see manual)

8: [DEPRECATED] Gromacs force field with hydrogens for NMR 9: Encad all-atom force field, using scaled-down vacuum charges 10: Encad all-atom force field, using full solvent charges ------------------------------------------------------------------------

本人建议使用OPLS力场和tip4p水模型。tip4p水有四个粒子，分别是两个氢原子，一个氧原子和一个没有质量的电粒子。这个电粒子在其他三个原子中间靠近氧原子。tip4p水模型多了一个粒子，模拟代价高一点，但是结果要好一点。但是最近有报道说使用spce水模型和GROMOS力场的计算结果最好，嘿嘿，好一个百家争鸣的学术氛围啊，我真的号感动哦。Gromacs也可以使用其他力场，如AMBER力场等，使用方法请参考google。 pdb2gmx的输出基本可以做真空中模拟了，现在对MD模拟的要求高，一般都要有点水。为分子系统添加水环境和离子环境需要其他命令，要慢慢来。

Gromacs的editconf命令在使用pdb2gmx创建模拟分子系统之后，可以使用editconf为你的分子画一个盒子。也可以认为使用editconf把分子放进一个盒子中，这样，你就可以往盒子里面添加水分子，离子，或者其他溶剂等等了。

使用“ eidtconf -h ”可以看到editconf的参数，其中比较常用的有以下几个： ---------------------------------------------------------- -f 指定你的坐标文件。

-n 分子系统的索引文件。索引文件是gromacs一个十分突出的功能，刚接触有点复杂，但是功能相当强大。 -o 输出文件，即放进盒子里面的分子系统。

-bt 盒子类型，有正方型，长方形，八面型等等，看个人需要跟癖好啦。

34

-box 自定义盒子大小，需要三个长度，即X、Y、Z三个方向的长度。

-d 分子离盒子表面的最短距离。这个跟-bt一起使用，基本就足够了;如果蛋白在模拟过程尺寸变化很大，那就用-box吧。

-center 确定哪一部分分子放在盒子中心，如果和索引文件一起使用，可以非常详细的定义分子的位置。 -translate 平移分子，跟X、Y、Z三个方向，随便移动。可以参考我前面一篇关于加大盒子的e文文章，练习e文，不好意思。

-rotate 转动分子，还是三个方向。我觉得editconf开发者太有才，什么都想到了，我想要什么，他就有什么，以后遇到他，我会让他给我签名。

-princ 这个参数可以用来对齐分子，比如使分子沿X轴对齐。举一个例子吧，比如你想将分子中两个残基沿Y轴对齐，那么就在索引文件中将这俩个残基标记以下，然后使用-princ，根据提示走就能对齐分子啦。 ----------------------------------------------------------

由于gromacs的很多命令都可以接受不同的文件类型，editconf也有其他功能，如和trjconv一样进行gro和pdb的转化：

editconf -f sen.pdb -o sen.gro

或者

editconf -f sen.gro -o sen.pdb

使用tpbconv重启gromacs模拟在使用gromacs的mdrun进行模拟计算过程中，很多因素可以是模拟计算终止。比如突然断电，断网或者磁盘空间满，或者windows死机（^_^）等等。重启gromacs模拟计算是一件十分方便的事情，因为gromacs众多的程序里面就有一个专门（或者吧）用来修改tpr文件的，就是tpbconv。

gromacs把模拟需要的所以文件都打包成一个tpr二进制文件，里面包含了分子坐标，各个原子在给定温度下速度和能量的分布。当模拟突然终止时，只要将终止时候系统的状态，即各个原子的位置、速度、坐标等装入tpr文件即可。tpbconv的参数也不少，可以使用\-h \查看，但是制作一个重启tpr文件的参数和格式一般如下：

tpbconv -s topol.tpr -f traj.trr -e ener.edr -o newtopol.tpr

其中topol.tpr为原来的tpr文件，traj.trr为双精度坐标文件（不要用xtc文件，因为精度不够），ener.edr为系统能量输出文件，newtopol.tpr是重启模拟文件。以上的命令得到的是在计算突然终止前一个系统构象的信息。也可以在命令中加上一个\参数来指定从那一个时间重新开始，如一下指定从一纳秒处重新开始模拟：

tpbconv -s topol.tpr -f traj.trr -e ener.edr -time 1000 -o newtopol.tpr

同时，如果模拟正常结束，而模拟时间让人觉得不够长时，可以使用tpbconv写一个延长模拟的tpr文件，一般格式如下：

tpbconv -s topol.tpr -f traj.trr -e ener.edr -extend 1000 -o newtopol.tpr

其中\表示延长1000ps的模拟时间。呵呵，非常好用。

这样断了又开始，就会产生很多轨迹文件，分析的时候非常不方便，gromacs有其他常用的命令把坐标文件，能量文件连接成一个文件，其中比较常用的如trjcat和eneconv，格式分别如下：

trjcat -f traj1.trr traj2.trr.... -o traj_all.trr eneconv -f ener1.edr ener2.edr... -o ener_all.edr

即使用\读入所有轨迹或者能量文件，使用\输出完整的轨迹和能量文件。

最后说说一个tpbconv的弱点。tpbconv不能更改你原来tpr文件中并行计算的节点数，比如你原来的tpr文件是8个节点的，那么使用tpbconv得到的重启tpr文件也是8个节点的。如果想更改使用节点数，那只能用grompp

35

重新做一个了。但是使用grompp做重启模拟文件时，就算你指定了原来的轨迹文件和能量文件，它还是会根据麦克斯韦分布重新给各个原子指定速度，真气人。

嗯，如果你觉得这是一个大问题，那就伸长脖子等gromcas新版

使用genbox命令为Gromacs模拟分子添加水环境使用editconf把分子放进一个盒子之后，接下来要做的就是往盒子里面添加水分子。Gromacs中添加水环境的命令是genbox，这是一个较其他命令简单一点的命令，因为参数不多。

通常用到的genbox参数有一下几个：

-cp ：带盒子参数的分子坐标文件，也就是editconf的输出文件；

-cs ：添加的水分子模型，如spc216、spce、tip3p、tip4p等，关于各个模型的区别，请参考scholar google； -o ：输出坐标文件，就是添加水分子之后的分子坐标文件，默认是.gro文件，但是也可以输出其他文件格式，如pdb；

-p ：系统拓扑文件，genbox会往里面写入添加水分子的个数，这个不要忘记，不然在进行下一步计算时，会出现坐标文件和拓扑文件原子数不一致的错误；

-ci ：索引文件，genbox可以为分子特定部位添加水环境，这样可以减少原子数，只有研究的分子部位添加水环境，节省计算时间。

-seed ：随机种子数，添加水分子时，各个水分子的位置是随机的，可以改变这个随机数子使水分子重新分布。添加水分子后可以用VMD等软件看看结果，一般完成这个之后事情开始变得复杂。比如某一个水分子出现在蛋白结构中，而这个位置本来就是不希望水分子一开始就存在，那么可以找出这个水分的残基标号，进行删除，同时删除拓扑文件中水分子的数目等。

使用grompp提取上一次模拟最后速度和能量在上面提到的使用gromacs程序包中tpbconv命令制作新的tpr文件中，最后提到新制作的.tpr文件只能使用跟原来.tpr文件一样多的CPU数目。还抱怨说这是tpbconv一个不足的地方。使用grompp可以制作一个新的.tpr文件，从上一步模拟的轨迹文件中提取速度，并从上一步能量文件中提取能量，也可以无缝的链接重启模拟计算。

要做到从上一步的最后的一个系统状态开始新的模拟计算。首先要在.mdp文件中把“ gen_vel ”参数定义为\\，这样做是为了告诉grompp不要重新为系统中的原子指定随机速度。指定新模拟开始的时间，即修改\参数。然后可以使用一下命令制作一个从上一步模拟文件中提取速度和能量的.tpr文件： ----------------------------------

grompp -f [.mdp文件] -c [上一步模拟最后的系统坐标文件] -p [拓扑文件] -t [上一步的trr轨迹文件] -e [上一步能量文件] -time [坐标文件对应的模拟时间] -o [输出tpr文件] -np [CPU数目]

----------------------------------

提取上一步模拟系统的速度时使用trr文件，是因为xtc为单精度，没有trr文件精确。\参数告诉grompp在上一步模拟文件中提取该时间的能量和速度，所以该时间要和系统的坐标文件相一致。

看起来好像要比tpbconv命令复杂一点，但是可以改变CPU数目，还算十分灵活。Gromacs是灵活的人的MD工具

(20)转自 sobereva 大师姐的blog，

(21)http://hi.http://www.wodefanwen.com//sobereva/blog/item/124afdd4e5979308a18bb7aa.html

(22)这么好的东西sob 竟然不弄到这里来，看来她不管我们了，。 2009-06-20 19:10

本文主要讨论隐式溶剂模型和GB模型（广义Born模型）的各个方面。

隐式水模型优势主要有下：

减少体系的粒子数从而降低计算开销。 没有溶剂摩擦效应，可以加快构象变化。

不显著、具体地表现溶剂运动而直接表现平均效果，便于计算自由能，不需要动力学很长时间采样。 不需要显式水的预平衡过程。

适用于常PH模拟、REMD等方法。

36

适合计算能量面，结果比较光滑，真实的水会带来噪音，溶剂的运动、碰撞使主体分子产生许多能量局部极小点。等等......

总体来说，隐式水模型相对于显式水模型，是一个离散化到连续化的过程。溶剂模型的近似关系如下（越左边越真实，越右边近似程度越大）：

显式溶剂模型-->隐式溶剂模型-->明确拆分为静电、非极性贡献-->PB-->GB

分子在溶剂环境下的总能量Etot=溶剂化能deltaG(solv)+溶质在真空条件下的能量Evac。（对于极化力场，溶剂对溶质的极化作用造成的能变也算进在deltaG(solv)里。）

deltaG(solv)=deltaG(elec)+deltaG(nonpolar)=deltaG(elec)+deltaG(vdw)+deltaG(cavity) 注：有时明确包含氢键项。

其中deltaG(elec)就是溶剂环境下的静电作用能与真空下的静电作用能之差，主要包括溶剂与溶质间的静电作用，以及溶剂对溶质原子间静电作用的屏蔽效果。

其中deltaG(nonpolar)就是除此以外，即体系各原子皆无电荷时，溶剂环境与真空环境下能量之差。它又分为deltaG(vdw)和 deltaG(cavity)。deltaG(vdw)是溶质与溶剂的范德华作用，由于范德华作用随距离衰减快，所以只考虑第一溶剂化层的溶剂，此项造成能量降低。deltaG(cavity)表现的是溶质反抗溶剂压力形成洞穴的做功以及溶质周围溶剂（第一溶剂化层为主）的重构造成的熵罚，此项造成能量升高。

因为deltaG(nonpolar)的两项都主要涉及第一溶剂化层，它正比于SASA，故最简单、常用的表达是a*SASA+b，而比例系数a由拟合实验值得到，b一般为0。另有其它方案计算，比如包括溶质体积、溶质直径等信息在内的方程。也有比如deltaG(vdw)通过SAS的积分得到，而 deltaG(cavity)仍用正比于SASA的方案。在GB模型下的MD模拟中，对于原始状态的生物大分子，由于SASA变化很小，原子因非极性作用导致的受力往往可以忽略，若需计算的话则是求它对坐标的导数。由于a*SASA+b的a为正值（一般为7.2cal/mol/?^2），可见 deltaG(nonpolar)的效果是减小溶质的SASA，此项体现了疏水效应。

deltaG(elec)项一般通过GB或PB方法计算，GB速度快但一般准确度低于PB。目前常用的GB模型计算方程由Still于1990年提出（原文http://www.brsbox.com/filebox/down/fc/90f984299e8583a61e0f33942bd13260），依据是库仑定律和born方程，推导出另一种形式。当粒子相离较远时，方程等价于库仑定律+born方程能量，而相距为0时等价于born方程，而当距离较近时等价于Onsager反应场能量。方程中还可以再引入单价粒子的屏蔽校正项。

GB方程中重要参数是有效born半径。某个原子电荷与溶剂作用对deltaG(elec)的贡献实际上是把分子中其它原子电荷去掉后，溶剂化状态与真空状态能量差。如果用的有效born半径以GB方程算的结果与之一样，就是正确的born半径。或者粗略来说：某原子用有效born半径算得的GB式中的“自身项”的能量应当与真实的分子中此原子与溶剂的静电相互作用能一致。这样，有了正确的born半径，就可以用GB式得到合理的 deltaG(elec)。有效born半径对计算deltaG(elec)有重要影响，模拟过程中分子构象不断改变，依赖于它的有效born半径也不断改变，需要每一步重新计算，若构象变化不明显可每隔一定步数重新计算，这是GB计算量中的主要部分之一。计算有效born半径在不同程序中有不同方法。一般使用CFA（库仑场近似），有效born半径就可以对溶质分子表面内的空间进行积分获得。而分子内空间的范围不便于描述，简单的方法是直接将原子VDW球内区域叠加作为分子内空间，称GB-HCT，但这样就忽略了原子VDW空间之间的缝隙，可以乘以4/3来近似校正。也有人根据原子被埋程度的概念提出了基于经验的简单、快速的函数，其中引入可调参数，在amber里称为GB-OBC方法，结果明显优于GB-HCT，应注意可调参数是通过预先优化得到的，面向不同体系模拟应使用不同参数。amber中还支持更新的GBn方法，结果优于GB-OBC，适合蛋白而不适于核酸体系。

GB方程的优点是形式简单，计算快速，结果较准确，而且有简单的解析导数，可直接计算GB环境下原子在MD过程中的受力。GB可以结合分子力场及量子力学模拟各种分子体系的溶剂化效果。可以用在类似MM/PBSA的结合自由能计算的静电项中，即MM/GBSA。一些模拟方法则只能使用GB，比如 amber中可以实现的结合MC的常PH模拟，原理是根据Metropolis判据在MD过程中动态决定是否某点被质子化或去质子化。对于REMD，因为随着体系粒子数的增加，需要的副本数目飙升（否则能量重叠较差，交换概率低而起不到效果），所以总是结合GB方法来显著降低粒子数量。对于不大的体系，GB比显式溶剂模型有更快的速度。

但GB方法也有缺点。与显式水模型仍有一定差距，将水分子进行了“连续化”的近似，故不能表现与溶剂的强相互作用、特殊相互作用（如水桥）。对于膜体系，溶剂、溶质、膜都有着不同的介电常数，这样复杂的情况介电常数环境下GB也难以表达。对于电荷较多的核酸体系，也不能较好表达多价抗衡粒子与它的相互作用。一些研究也表明用GB模拟多肽会产生与显式溶剂模型不同的构象，与实验结果有一定偏差。计算deltaG(elec)项

37

时，相对于原理严格的PB 也有差距（即便使用最佳的有效born半径），但是好处是GB计算速度快得多。不过PB的结果亦不可作为绝对的金标准，因为一些误差在PB引入的近似中就已经出现了。在计算速度上，显式溶剂体系计算量是O(N^2)，但使用Ewald等方法后，可降至Nlog(N)以下，而GB模型如果在cutoff上不做近似（往往取无限大），对于大体系反倒更慢。

除GB、PB以外也有其它一些方案，最简单的隐式模型是介电常数依赖于作用距离的方法，其中库仑作用的介电常数等于粒子间距离（即静电作用能1/r 变为了1/(r^2)），以体现溶剂的屏蔽效果，此方法精度显然低于GB。近来还有一些新方法，如AGBNP（Analytical Generalized Born plus Nonpolar），在GB基础上引入另外形式的deltaG(nonpolar)项。ALPB（analytical linearized Poisson-Boltzmann），其方程类似GB（在无限介电常数下回归为GB方程），故有着基本一致的计算速度，但引入了有效分子静电尺寸参数，模拟水比GB有着更好的精度，结果更接近于PB。这些新模型理论上有着更好的结果，但仍需广泛地应用于模拟来检验。还有隐式与显式溶剂模型相结合的方法，也就是在主体分子外面包一层溶剂分子，往往以弹簧势限制住溶剂避免跑走，而更外面则用隐式溶剂模型，结合了两种方法的一些优点。 目前Amber对隐式溶剂模型支持较好，支持GB、ALPB，也有内建的PB模块pbsa（老版本挂的是第三方的delphi）。而Gromacs在最新版本4.0.5中尚不支持GB/PB，需要外接第三方软件如APBS，但在接下来的版本中即将加入GB。

对于GB模型，适用领域与不适用领域的界限尚为模糊，一般来说，对于必须用GB的地方，或者已证明有明显优势的地方可以使用，但如果不确定，尽量还是用显式溶剂模型，毕竟GB是基于多层近似后的溶剂模型。

g_rms出来的xmgrace图有波谷！

38

刚完成一个蛋白在水中的simulation，总共有四个xtc文件，分别是pr.xtc, ann.xtc (5ns), run.xtc (10ns), run.xtc (40ns)。然后我利用如下command组合了一个xtc:

trjcat -f run.xtc run1.xtc -settime -o water50.xtc, 运行后选择的settime是0， C; 然后用g_rms -s run1.tpr -f water50.xtc -n index.ndx -o water50.xvg; 最后用xmgrace做出如上图。可以发现在连接处10ns的地方有个波谷，非常奇怪，问过教授，他说是因为做trjcat联合的时候应该从pr energy minimized的文件开始。一时不知道怎么回事？

小弟是初学者，还请各位帮帮忙，究竟是怎么回事？正确的做法又是什么？谢谢！

39

本文来源：https://www.bwwdw.com/article/gvyw.html