【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. 下载力场 :
  1. 解压缩

  2. 把aminoacides.dat、vdwradii.dat拷到top文件夹下覆盖,实际上一股脑考过去就完了。

  3. 子文件夹下的东西都拷到外面,也就是top里面。

  4. 修改FF.dat,修改第一行的数字,加上新添的力场数。

  5. 为每个新的力场增加新行,如: 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)。

药企,独角兽,苏州。团队长期招人,感兴趣的都可以发邮件聊聊:tiehan@sina.cn
个人公众号,比较懒,很少更新,可以在上面提问题,如果回复不及时,可发邮件给我: tiehan@sina.cn