【3.3.3】gromacs力场

“力场”，请不要被“场”这个听起来像是十分高深的物理名词给吓坏了。

分子动力学模拟中使用的力场，包含两个重要的部分：

1。 模拟粒子之间相互作用的方程（即经典力学的相互作用力方程如：库仑定律，范德华作用方程等）。 2. 方程的参数（即各个不同粒子，原子本身的参数，如带点量等等）。

可以想想，模计算机模拟好多成键或者不成键的粒子的运动，总要让它们互相推推拉拉吧，于是力场就是定义它们推推拉拉的方式（按照物理定律）。

一、介绍

1.1 力场类型

一般分类为三种：

1. 全原子力场：精确定义每一个原子的参数。
2. 联合原子力场：省略非极性氢原子，同时把其参数整合到与他们成键的相邻原子上（比如甲基，只由一个碳原子表示）。
3. 粗颗粒力场：进一步精简分子结构的力场参数，种类比较多，比如有讲蛋白侧链看作一个颗粒的力场，或者甚至将整个氨基酸残基看成一个颗粒的力场等等。

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

1.2 目前比较流行的力场有：

1. AMBER： 包含好几个版本的力场，为全原子力场；
2. CHARMM：全原子力场，是软件CHARMM的一部分；
3. GROMOS：GROMOS软件使用的力场，版本较多，为联合原子力场；OPLS：包含全原子和联合原子力场两个版本； 粗颗粒力场：种类较多，没有固定版本或者种类，一般根据研究需要开发。

二、下载

2.1 在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

三、讨论

3.1 gromacs决定力场参数的方式

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的情况下有两种情况：

1. 对于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却发现连接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还是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)。