文件目录
最近在研究蛋白质复合物的时候遇到了一系列与力场有关的麻烦事，所以决定找个时间研究一下力场文件究竟写了些什么，免得以后再看不懂报错出来的都是啥。
本文以 gromos54a7.ff 力场（你可以在 GROMACS安装目录/share/gromacs/top 的子目录找到）为例，对力场文件进行了简要的解读。在目录中，括号内的序号标志着该文件及其拓展内容可以在GROMACS手册对应章节找到。
我们的 gromos54a7.ff 下文件有这些：
atomtypes.atp(5.2.1)
具体内容如下：
O 15.99940 ; carbonyl oxygen (C=O)
OM 15.99940 ; carboxyl oxygen (CO-)
OA 15.99940 ; hydroxyl, sugar or ester oxygen
OE 15.99940 ; ether or ester oxygen
OW 15.99940 ; water oxygen
N 14.00670 ; peptide nitrogen (N or NH)
NT 14.00670 ; terminal nitrogen (NH2)
NL 14.00670 ; terminal nitrogen (NH3)
NR 14.00670 ; aromatic nitrogen
NZ 14.00670 ; Arg NH (NH2)
NE 14.00670 ; Arg NE (NH)
.atp 是 atom types 的简写，它是力场对不同环境的原子的命名与标注。第一列为原子名称，第二列为原子质量，第三列（一般存在）为简要的解释。力场中所有文件（如后文提到的 .rtp 文件）对于原子的命名必须与 .atp 文件相统一，否则力场就无法识别该原子。
aminoacid related files
ff目录下前六个文件是与氨基酸相关的参数文件，为了照顾逻辑关系，我们没有按照上图的顺序进行解释：
aminoacids.rtp(5.6.1)
基本内容如下（省略了大部分对其他氨基酸的定义）：
[ bondedtypes ]
; bonds angles dihedrals impropers
2 2 1 2
[ ALA ]
[ atoms ]
; name type charge chargegroup
N N -0.31000 0
H H 0.31000 0
CA CH1 0.00000 1
CB CH3 0.00000 1
C C 0.450 2
O O -0.450 2
[ bonds ]
; ai aj gromos type
N H gb_2
N CA gb_21
CA CB gb_27
CA C gb_27
C O gb_5
C +N gb_10
[ exclusions ]
; ai aj
[ angles ]
; ai aj ak gromos type
-C N H ga_32
-C N CA ga_31
H N CA ga_18
N CA CB ga_13
N CA C ga_13
CB CA C ga_13
CA C O ga_30
CA C +N ga_19
O C +N ga_33
[ impropers ]
; ai aj ak al gromos type
N -C CA H gi_1
CA N C CB gi_2
C CA +N O gi_1
[ dihedrals ]
; ai aj ak al gromos type
-CA -C N CA gd_14
-C N CA C gd_44
-C N CA C gd_43
N CA C +N gd_45
N CA C +N gd_42
.rtp 文件是 residue topology parameters 的缩写。它主要包含了对各种常见氨基酸（以及部分小分子物质如CCl4）的结构参数。
文件的的第一行均为[ bondedtypes ]，后面的四个数字分别表示着键、键角、二面角、异常二面角的相互作用类型（目前笔者尚未查询到有关各数字所对应作用类型的信息）。此后各残基以块的结构书写在文件里，每个块包含[ name ]、[ atoms ]（必需）和[ bonds ]、[ exclusions ]、[ angles ]、[ impropers ]、[ dihedrals ]（非必需）的信息。在非必需的信息中，出现原子名称前有-、+意味着该原子来自此残基的后一个或前一个残基。其余各参数释义标注在了代码片内（关于gromos type的信息请见本文 ff files-ffnonbonded.itp 部分）。
aminoacids.r2b(5.6.2)
该文件的基本内容如下：
; rtp residue to rtp building block table
;GMX Force-field
CYS CYSH
HISD HISA
HISE HISB
LYS LYSH
LYSN LYS
HEM HEME
.r2b是 residue to building block 的缩写（正如文件第一行注释所写），用来转换名称以便于力场识别。第一列为残基名称（具体释义可参考手册），第二列为构建单元名称。
aminoacids.hdb(5.6.4)
文件基本内容如下：
ALA 1
1 1 H N -C CA
ARG 4
1 1 H N -C CA
1 1 HE NE CD CZ
2 3 HH1 NH1 CZ NE
2 3 HH2 NH2 CZ NE
ARGN 4
1 1 H N -C CA
1 1 HE NE CD CZ
1 2 HH1 NH1 CZ NE
2 3 HH2 NH2 CZ NE
ASN 2
1 1 H N -C CA
2 3 HD2 ND2 CG CB
.hdb是 hydrogen database 的缩写，储存着如何将H原子连接到已有原子的信息。我们将.rtp文件和.hdb文件的基本结构块放在一起看（以ARG为例）如下：
在.hdb文件中：第一行是残基名称和需要添加的H的原子的数目（可以对比.rtp文件发现，只有与主链上的重原子C相连的H没有相应信息，与O/N/S等相连的均已写在了 .rtp 文件中）；剩余每行对应需要添加H原子的具体信息：第一列至右依次代表需要添加的H原子数目、添加方式（稍后讲解）、H原子的名称、三个（或四个）用于定位H（n）的辅助原子（i/j/k）。
;from .hdb
ARG 4
1 1 H N -C CA ;例子，见代码块下
1 1 HE NE CD CZ
2 3 HH1 NH1 CZ NE
2 3 HH2 NH2 CZ NE
;from .rtp
[ ARG ]
[ atoms ]
N N -0.31000 0
H H 0.31000 0
CA CH1 0.00000 1
CB CH2 0.00000 1
CG CH2 0.00000 1
CD CH2 0.09000 2
NE NE -0.11000 2
HE H 0.24000 2
CZ C 0.34000 2
NH1 NZ -0.26000 2
HH11 H 0.24000 2
HH12 H 0.24000 2
NH2 NZ -0.26000 2
HH21 H 0.24000 2
HH22 H 0.24000 2
C C 0.450 3
O O -0.450 3
……
具体的添加方式代号代表的含义请参考手册，这里举一个例子：ARG氢原子生成的第一个是在α-C上生成的，1代表产生一个 n，置于 ijk 形成的平面内并位于 j-i-k 角平分线距 i 0.1nm 处，同时 n-i-j 和 n-i-k 均大于90°；这是一个典型的肽主链上的单原子H。
aminoacids.c.tdb(5.6.5)
基本内容如下：
[ None ]
[ COO- ]
[ replace ]
C C C 12.011 0.27
O O1 OM 15.9994 -0.635
OXT O2 OM 15.9994 -0.635
[ add ]
2 8 O C CA N
OM 15.9994 -0.635
[ bonds ]
C O1 gb_6
C O2 gb_6
[ angles ]
O1 C O2 ga_38
CA C O1 ga_22
CA C O2 ga_22
[ dihedrals ]
N CA C O2 gd_45
N CA C O2 gd_42
[ impropers ]
C CA O2 O1 gi_1
[ COOH ]
[ replace ]
C C C 12.011 0.33
O O OA 15.9994 -0.288
OXT OT O 15.9994 -0.45
[ add ]
1 2 OT C CA N
O 15.9994 -0.45
1 2 HO O C CA
H 1.008 0.408
[ bonds ]
C O gb_13
C OT gb_5
O HO gb_1
[ angles ]
O C OT ga_33
C O HO ga_12
CA C O ga_19
CA C OT ga_30
[ dihedrals ]
N CA C O gd_45
N CA C O gd_42
CA C O HO gd_12
[ impropers ]
C CA OT O gi_1
文件名中 c 则意味着该文件与肽链C端氨基酸参数有关；而 .tdb 为 terminal database（末端数据库）的意思。这个文件会被 pdb2gmx 调用，用来告诉程序如何将末端原子连接到已有原子、应该删除或更改那些原子、还应该添加那些键相互作用。
可以看出，文件的基本结构单元是长这样的：
[ COO- ]
[ replace ]
C C C 12.011 0.27
O O1 OM 15.9994 -0.635
OXT O2 OM 15.9994 -0.635
[ add ]
2 8 O C CA N
OM 15.9994 -0.635
[ bonds ]
C O1 gb_6
C O2 gb_6
[ angles ]
O1 C O2 ga_38
CA C O1 ga_22
CA C O2 ga_22
[ dihedrals ]
N CA C O2 gd_45
N CA C O2 gd_42
[ impropers ]
C CA O2 O1 gi_1
手册里面将这种单元称之为块。当 pdb2gmx 在读取肽链上的每一个氨基酸时，它会自动判断氨基酸所处位置与需要设定的类型。例如位于肽链中间部分的氨基酸会被判定到[ None ]下而不进行任何改变；对于C末端的氨基酸则会判定到[ COOH ]或[ COO- ]进行相应更改。每个块下又有细分的信息：
- [ replace ]：用于提供替换原子所需的信息：
[ replace ]
O（要被替换的原子） O1（新的原子名称） OM（新原子的类型） 15.9994（新原子的质量） -0.635（新原子的电荷）
- [ add ]：用于添加新的原子；第一行与.hdb文件基本格式相似，只不过添加类型多了8/9两种（且为C末端专用，具体请参见手册），第二行用于声明所添加原子的具体信息：
[ add ]
2（需要添加的原子数） 8（添加类型） O（i/j/k/l 原子） C CA N
OM（被添加原子的类型） 15.9994（质量） -0.635（电荷） （可选，电荷组）
- [ delete ]：用于删除已有原子；每个原子一行：
[ delete ]
H
- 其余键参数、角参数与.rtp文件中格式相同；用于添加额外的键参数。
aminoacids.n.tdb
与 aminoacids.c.tdb 相同，只不过为N端。
aminoacids.vsd(5.6.6)
基本内容如下：
; no CH3 groups - Gromos96 is a united atom forcefield
[ NH3 ]
NL C MNH3
NL CH1 MNH3
NL CH2 MNH3
[ NH2 ]
NT planar
NZ planar
NL C MNH3
NL CH1 MNH3
NL CH2 MNH3
; Data for generating dummy aromatic rings.
; Actually we dont need all these bonds and angles,
; but by specifying them here it is easier to improve
; the dummy generation code later.
[ PHE ]
CG CD1 0.139
CG CD2 0.139
CD1 CE1 0.139
CD2 CE2 0.139
CE1 CZ 0.139
CE2 CZ 0.139
CD1 HD1 0.109
CD2 HD2 0.109
CE1 HE1 0.109
CE2 HE2 0.109
CZ HZ 0.109
CG CD1 CE1 120.0
CD1 CE1 CZ 120.0
CE1 CZ CE2 120.0
CZ CE2 CD2 120.0
CE2 CD2 CG 120.0
CD2 CG CD1 120.0
CG CD1 HD1 120.0
CG CD2 HD2 120.0
HD1 CD1 CE1 120.0
CD1 CE1 HE1 120.0
HE1 CE1 CZ 120.0
CE1 CZ HZ 120.0
HZ CZ CE2 120.0
CZ CE2 HE2 120.0
HE2 CE2 CD2 120.0
HD2 CD2 CG 120.0
.vsd 是 virtual site database 的缩写，用于创建拓扑中的虚拟位点（可以参考手册4.7）。文件分为两个部分：第一部分确定了用于基团的质心，第一列向右依次为连接到 2/3 氢原子的原子类型、连接下一个重原子的原子的类型、所使用的质心类型（MNH3为经典四面体类型，planar为平面类型）；第二部分为芳香族侧链的参数：2或3列原子名称、指定键长（2列）或键角（3列）参数。
summarizing briefly
在 pdb2gmx 构建拓扑文件过程中，对于蛋白质而言，其主要应用上述文件对蛋白质进行拓扑构建。首先是读取 aminoacids.rtp 对各个氨基酸残基进行基本的处理，顺带用 aminoacids.r2b 进行名称规整，而后借助 aminoacids.hdb 添加缺失的H原子并用 aminoacids.c/n.tdb 对末端、aminoacids.vsd 对刚性结构进行处理；整个过程并非顺序执行，而有可能是并行存在的。
ff files
ffbonded.itp
.itp 是 independent topology parameter 的简写，而 ffbonded 则意味着其对力场键参数做出了定义。
参考下文的代码块（为节省空间删去了一些空的注释行）可以发现：
- 第一个主类是键拉伸参数（Table 2.5.2.1，bond-stretching parameters），代码简记为 gb，#define 语句对键进行了详细的定义：第一个参数为键的代号（Bond type code）、第二个为力常数（Force constant）、第三个为理论键长（Ideal bond length）、语句的下一行注释为符合这一键类型的例子；
- 第二个主类是键角参数（Table 2.5.3.1，bond-angle bending parameters），代码简记为 ga，在 table 下也有对 #define 参数的解释；
- 第三个主类是异常二面角参数（Table 2.5.4.1，improper (harmonic) dihedral angle parameters），代码简记为 gi，同样有对参数的注释；
- 第四个主类为二面角扭转参数（Table 2.5.5.1，(trigonometric) dihedral torsional angle parameters），代码简记为 gd，同样有对参数的注释；
- 最后还有一部分对虚拟位点参数的调用（
#include "ff_dum.itp"，ff_dum.itp 文件将会在后文解释）、对特殊结构的约束和对特殊键型的补充定义，在定义前均有注释行进行解释；其中[ constrainttypes ]，（还没研究明白，以后再写）
; Table 2.5.2.1
; GROMOS bond-stretching parameters
;
; Bond type code
; Force constant
; Ideal bond length
; Examples of usage in terms of non-bonded atom types
;
; ICB(H)[N] CB[N] B0[N]
;
#define gb_1 0.1000 1.5700e+07
; H - OA 750
;
; ……
;
; Table 2.5.3.1
; GROMOS bond-angle bending parameters
;
; Bond-angle type code
; Force constant
; Ideal bond angle
; Example of usage in terms of non-bonded atom types
;
; ICT(H)[N] CT[N] (T0[N])
;
#define ga_1 90.00 380.00
; NR(heme) - FE - C 90
;
; ……
;
; Table 2.5.4.1
; GROMOS improper (harmonic) dihedral angle parameters
;
; Improper dihedral-angle type code
; Force constant
; Ideal improper dihedral angle
; Example of usage
;
; ICQ(H)[N] CQ[N] (Q0[N])
;
#define gi_1 0.0 167.42309
; planar groups 40
;
; ……
;
; Table 2.5.5.1 (Note: changes with respect to the 43A1 table)
; GROMOS (trigonometric) dihedral torsional angle parameters
;
; Dihedral-angle type code
; Force constant
; Phase shift
; Multiplicity
; Example of usage in terms of non-bonded atom types
;
; ICP(H)[N] CP[N] PD[N] NP[N]
;
#define gd_1 180.000 2.67 1
; CHn-CHn-CHn-OA (sugar) 0.6
;
; ……
;
; get the constraint distances for dummy atom constructions
#include "ff_dum.itp"
[ constrainttypes ]
; now the constraints for the rigid NH3 groups
MNH3 C 2 DC_MNC1
MNH3 CH1 2 DC_MNC2
MNH3 CH2 2 DC_MNC2
MNH3 MNH3 2 DC_MNMN
; and the angle-constraints for OH and SH groups in proteins:
CH2 H 2 DC_CO
CH1 H 2 DC_CO
C H 2 DC_CO
P H 2 DC_PO
; bond-, angle- and dihedraltypes for specbonds:
[ bondtypes ]
S S 2 gb_36
NR FE 2 gb_34
[ angletypes ]
CH1 CH2 S 2 ga_16
CH2 S S 2 ga_6
CR1 NR FE 2 ga_34
NR FE NR 2 ga_17
[ dihedraltypes ]
S S 1 gd_21
NR FE 1 gd_38
CH2 S 1 gd_26
ffnonbonded.itp(5.3.2/5.3.4)
- 先看[ atomtypes ]部分：由于原子的静态性质中，质量来自 atomtypes.atp 文件，电荷来自于相应的 .rtp 文件，所以在非键参数的定义中质量与电荷被标为了0；第一列向右依次为为原子名称、原子序号、质量、电荷、粒子类型（A代表原子，D或V代表虚拟相互作用位点，详细见手册5.2）、范德华参数V（或 σ σ σ）、W（或 ε ε ε）。
- 第二部分为非键参数定义：其中 func 代表非键函数类型（1为Lennard-Jones，2为Buckingham），参数 c6，c12 的含义取决于拓扑文件（.top）中 [ defaults ] 的 comb-rule 定义（详细请见手册5.3.2与5.7.1[ defaults ]部分）。
- 第三部分为对相互作用的定义（详细解释参考手册5.3.4）
[ atomtypes ]
;name at.num mass charge ptype c6 c12
O 8 0.000 0.000 A 0.0022619536 1e-06
OM 8 0.000 0.000 A 0.0022619536 7.4149321e-07
OA 8 0.000 0.000 A 0.0022619536 1.505529e-06
OE 8 0.000 0.000 A 0.0022619536 1.21e-06
OW 8 0.000 0.000 A 0.0026173456 2.634129e-06
;
; ……
;
[ nonbond_params ]
; i j func c6 c12
OM O 1 2.261954E-03 8.611000E-07
OA O 1 2.261954E-03 1.386510E-06
OA OM 1 2.261954E-03 2.258907E-06
OE O 1 2.261954E-03 1.100000E-06
OE OM 1 2.261954E-03 9.472100E-07
OE OA 1 2.261954E-03 1.505529E-06
OW O 1 2.433170E-03 1.833990E-06
OW OM 1 2.433170E-03 2.987943E-06
OW OA 1 2.433170E-03 1.991421E-06
OW OE 1 2.433170E-03 1.991421E-06
;
; ……
;
[ pairtypes ]
; i j func c6 c12
O O 1 2.261954E-03 7.414932E-07
OM O 1 2.261954E-03 7.414932E-07
OM OM 1 2.261954E-03 7.414932E-07
OA O 1 2.261954E-03 9.687375E-07
OA OM 1 2.261954E-03 9.687375E-07
OA OA 1 2.261954E-03 1.265625E-06
OE O 1 2.261954E-03 9.687375E-07
OE OM 1 2.261954E-03 9.687375E-07
OE OA 1 2.261954E-03 1.265625E-06
OE OE 1 2.261954E-03 1.265625E-06
ff_dum.itp
该文件定义了一些用于处理虚拟位点的参数，如下：
; These constraints are used for vsite constructions as generated by pdb2gmx.
; Values depend on the details of the forcefield, vis. bondlengths and angles
; These parameters are designed to be used with the GROMOS96 forcefields
; G43a1, G43a2 and G43b1.
; Constraints for the rigid NH3/CH3 groups depend on the hygrogen mass,
; since an increased hydrogen mass translates into increased momentum of
; inertia which translates into a larger distance between the dummy masses.
#ifdef HEAVY_H
; now the constraints for the rigid NH3 groups
#define DC_MNC1 0.175695
#define DC_MNC2 0.188288
#define DC_MNMN 0.158884
; now the constraints for the rigid CH3 groups
#define DC_MCN 0.198911
#define DC_MCS 0.226838
#define DC_MCC 0.204247
#define DC_MCNR 0.199798
#define DC_MCMC 0.184320
#else
; now the constraints for the rigid NH3 groups
#define DC_MNC1 0.144494
#define DC_MNC2 0.158002
#define DC_MNMN 0.079442
; now the constraints for the rigid CH3 groups
#define DC_MCN 0.161051
#define DC_MCS 0.190961
#define DC_MCC 0.166809
#define DC_MCNR 0.162009
#define DC_MCMC 0.092160
#endif
; and the angle-constraints for OH and SH groups in proteins:
#define DC_CS 0.23721
#define DC_CO 0.19849
#define DC_PO 0.21603
其中有一个gromacs文件中常见的条件语句：
#ifdef XXX
;#……
#else
;#……
#endif
我们便创建了一个名为 XXX 的选择功能，只需要在进入模拟之前的 .mdp 文件中使用define = -DXXX或在相应的 .top 文件前些部分使用#define XXX，就可以使模拟调用#ifdef XXX部分的参数；如果没有事先声明，则会调用#else部分的参数。
forcefield.itp
#define _FF_GROMOS96
#define _FF_GROMOS54A7
[ defaults ]
;非键函数类型 组合规则 对生成与否 L-J1-4修正因子 静电1-4修正因子
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 1 no 1.0 1.0
#include "ffnonbonded.itp"
#include "ffbonded.itp"
这个文件实际上标明了力场的[ defaults ]和#include设置。中文注释为笔者添加。
Molecule.itp files
ions.itp
该文件被用于生成离子的拓扑（主要是用在genion）。其中nrexcl表示排除不远于x条键的非键相互作用力。
[ moleculetype ]
; molname nrexcl
CU1 1
[ atoms ]
; id at-type res-nr residu name at-name cg-nr charge mass
1 CU1+ 1 CU1 CU 1 1 63.54600
etc
此外ff 内还包括了其他用于solvate或脂质生成的 .itp 文件：
水分子：
脂质分子：
其基本格式与 ions.itp 相同，且由于这部分与 .top 文件非常相似，我们将会在另一篇文章里面细说。
扫一扫
1351
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
