本文将以乙酰基噻唑分子为例，简单讲解GAFF力场小分子转换为Gromacs拓扑格式ITP/TOP应当包含的内容。
O S N cc cd c o c 3 ss cd nd GAFF2 首先我们使用以下命令，构建GAFF力场的Gromacs拓扑文件。
1 2 3 4 5 6 7 8 obabel -:"c1c(C(=O)C)scn1" --gen3d -addH -omol2 -O MOL.mol2 acpype -i MOL.mol2 -n 0 -a gaff2 -k "maxcyc=10"
然后在目录中会显示包括如下文件在内的力场文件：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 # 重要的Amber中间文件 MOL_AC.frcmod # AmberTools的pharmchk产生文件，主要关注成键参数的penalty score MOL_AC.prmtop # Amber格式的小分子力场文件 MOL_AC.inpcrd # Amber格式的小分子坐标文件 MOL_bcc_gaff2.mol2 # Amber原子类型和AM1-BCC电荷的MOL2格式 # 其他力场或软件格式 MOL_GMX.itp # Gromacs ITP拓扑文件 MOL_GMX.top # Gromacs TOP拓扑文件 posre_MOL.itp # Gromacs 的位置限制ITP拓扑文件 MOL_GMX.gro # Gromacs GRO坐标文件 MOL_GMX_OPLS.itp # Gromacs ITP拓扑文件，仅改用OPLS力场的原子类型和范德华参数 MOL_CNS.par # CNS 拓扑文件 MOL_CNS.top # CNS 拓扑文件 MOL_CNS.prm # CHARMM 拓扑文件 MOL_CNS.rtf # CHARMM 拓扑文件
其中MOL_GMX.itp和MOL_GMX.top就是我们需要的Gromacs格式的GAFF力场小分子拓扑。
Gromacs拓扑文件格式官方文档 对拓扑文件有着非常详尽的释义，我们将结合官方文档对MOL_GMX.itp和MOL_GMX.top进行解释。
首先是TOP文件全文：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ; MOL_GMX.itp created by acpype (v: 2022.6.6) on Fri Feb 2 12:27:46 2024 [ defaults ] ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ 1 2 yes 0.5 0.8333 #include MOL_GMX.itp [ system ] MOL [ molecules ] ; Compound nmols MOL 1
TOP文件包含体系中所有原子的子拓扑参数，由力场参数，子拓扑，和体系分子数三个主要部分组成。第一部分为力场参数，首先是[ defaults ]词条，词条下有5列参数，第一列为nbfunc，范德华势能函数类型，1代表Lennard-Jones势，2代表Buckingham势或EXP势，AMBER力场使用表Lennard-Jones势；第二列是范德华势的参数组合规则，AMBER使用的组合规则见第二章中GAFF力场参数说明；第三列gen-pairs为1,4相互作用对的参数来源，gen-pairs为yes时会开启Fudge LJ，为no时关闭Fudge LJ，Fudge QQ总是开启，不受次参数调整，第五和第六列为范德华和库仑势的具体缩放参数。如果TOP文件中含有大分子力场，[ defaults ]参数后会#include amber99sb-ildn/ffnonbonded.itp引入AMBER力场中大分子非键参数。
第二部分是各体系分子子拓扑文件，即ITP文件，使用#include方法引入，具体解析见下文。第三部分为体系分子数。首先是系统名[ system ]，命名可由用户自定义，然后是[ molecules ]的具体分子组成，每行代表一类分子的数量，如本体系中UNK分子数为1。对于复杂体系，可能会有多个蛋白或配体，排列顺序务必和坐标文件中的分子顺序保持一致。
GAFF力场小分子参数文件(MOL_GMX.itp)全文如下：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 [ atomtypes ] ;name bond_type mass charge ptype sigma epsilon Amb ss ss 0.00000 0.00000 A 3.53241e-01 1.18156e+00 ; 1.98 0.2824 cc cc 0.00000 0.00000 A 3.31521e-01 4.13379e-01 ; 1.86 0.0988 nd nd 0.00000 0.00000 A 3.38417e-01 3.93714e-01 ; 1.90 0.0941 cd cd 0.00000 0.00000 A 3.31521e-01 4.13379e-01 ; 1.86 0.0988 c c 0.00000 0.00000 A 3.31521e-01 4.13379e-01 ; 1.86 0.0988 o o 0.00000 0.00000 A 3.04812e-01 6.12119e-01 ; 1.71 0.1463 c3 c3 0.00000 0.00000 A 3.39771e-01 4.51035e-01 ; 1.91 0.1078 h4 h4 0.00000 0.00000 A 2.53639e-01 6.73624e-02 ; 1.42 0.0161 hc hc 0.00000 0.00000 A 2.60018e-01 8.70272e-02 ; 1.46 0.0208 [ moleculetype ] ;name nrexcl MOL 3 [ atoms ] ; nr type resi res atom cgnr charge mass ; qtot bond_type 1 ss 1 UNL S1 1 0.048700 32.06000 ; qtot 0.049 2 cc 1 UNL C2 2 0.267900 12.01000 ; qtot 0.317 3 nd 1 UNL N3 3 -0.605000 14.01000 ; qtot -0.288 4 cd 1 UNL C4 4 0.323200 12.01000 ; qtot 0.035 5 cc 1 UNL C5 5 -0.275400 12.01000 ; qtot -0.241 6 c 1 UNL C6 6 0.503600 12.01000 ; qtot 0.263 7 o 1 UNL O7 7 -0.534100 16.00000 ; qtot -0.271 8 c3 1 UNL C8 8 -0.181100 12.01000 ; qtot -0.452 9 h4 1 UNL H9 9 0.048100 1.00800 ; qtot -0.404 10 h4 1 UNL H10 10 0.186000 1.00800 ; qtot -0.218 11 hc 1 UNL H11 11 0.072700 1.00800 ; qtot -0.145 12 hc 1 UNL H12 12 0.072700 1.00800 ; qtot -0.073 13 hc 1 UNL H13 13 0.072700 1.00800 ; qtot 0.000 [ bonds ] ; ai aj funct r k 1 2 1 1.7560e-01 1.9414e+05 ; S1 - C2 1 5 1 1.7560e-01 1.9414e+05 ; S1 - C5 2 3 1 1.3170e-01 3.7715e+05 ; C2 - N3 2 6 1 1.4680e-01 2.4719e+05 ; C2 - C6 3 4 1 1.3690e-01 3.0886e+05 ; N3 - C4 4 5 1 1.3730e-01 3.4819e+05 ; C4 - C5 4 9 1 1.0820e-01 3.3798e+05 ; C4 - H9 5 10 1 1.0820e-01 3.3798e+05 ; C5 - H10 6 7 1 1.2180e-01 5.4610e+05 ; C6 - O7 6 8 1 1.5240e-01 2.0351e+05 ; C6 - C8 8 11 1 1.0970e-01 3.1455e+05 ; C8 - H11 8 12 1 1.0970e-01 3.1455e+05 ; C8 - H12 8 13 1 1.0970e-01 3.1455e+05 ; C8 - H13 [ pairs ] ; ai aj funct 1 7 1 ; S1 - O7 1 8 1 ; S1 - C8 1 9 1 ; S1 - H9 2 9 1 ; C2 - H9 2 10 1 ; C2 - H10 2 11 1 ; C2 - H11 2 12 1 ; C2 - H12 2 13 1 ; C2 - H13 3 7 1 ; N3 - O7 3 8 1 ; N3 - C8 3 10 1 ; N3 - H10 4 6 1 ; C4 - C6 5 6 1 ; C5 - C6 7 11 1 ; O7 - H11 7 12 1 ; O7 - H12 7 13 1 ; O7 - H13 9 10 1 ; H9 - H10 [ angles ] ; ai aj ak funct theta cth 1 2 3 1 1.1451e+02 6.9705e+02 ; S1 - C2 - N3 1 2 6 1 1.2197e+02 5.2718e+02 ; S1 - C2 - C6 1 5 4 1 1.1155e+02 5.5982e+02 ; S1 - C5 - C4 1 5 10 1 1.1997e+02 3.5564e+02 ; S1 - C5 - H10 2 1 5 1 9.0240e+01 6.5270e+02 ; C2 - S1 - C5 2 3 4 1 1.0549e+02 6.1840e+02 ; C2 - N3 - C4 2 6 7 1 1.2393e+02 7.2551e+02 ; C2 - C6 - O7 2 6 8 1 1.1729e+02 5.4392e+02 ; C2 - C6 - C8 3 2 6 1 1.2188e+02 7.1379e+02 ; N3 - C2 - C6 3 4 5 1 1.1165e+02 7.6232e+02 ; N3 - C4 - C5 3 4 9 1 1.2114e+02 5.1631e+02 ; N3 - C4 - H9 4 5 10 1 1.2848e+02 3.9999e+02 ; C4 - C5 - H10 5 4 9 1 1.2848e+02 3.9999e+02 ; C5 - C4 - H9 6 8 11 1 1.0877e+02 3.9664e+02 ; C6 - C8 - H11 6 8 12 1 1.0877e+02 3.9664e+02 ; C6 - C8 - H12 6 8 13 1 1.0877e+02 3.9664e+02 ; C6 - C8 - H13 7 6 8 1 1.2320e+02 7.0793e+02 ; O7 - C6 - C8 11 8 12 1 1.0758e+02 3.2635e+02 ; H11 - C8 - H12 11 8 13 1 1.0758e+02 3.2635e+02 ; H11 - C8 - H13 12 8 13 1 1.0758e+02 3.2635e+02 ; H12 - C8 - H13 [ dihedrals ] ; propers ; for gromacs 4.5 or higher, using funct 9 ; i j k l func phase kd pn 1 2 3 4 9 180.00 19.87400 2 ; S1- C2- N3- C4 1 2 6 7 9 180.00 12.02900 2 ; S1- C2- C6- O7 1 2 6 8 9 180.00 12.02900 2 ; S1- C2- C6- C8 1 5 4 3 9 180.00 16.73600 2 ; S1- C5- C4- N3 1 5 4 9 9 180.00 16.73600 2 ; S1- C5- C4- H9 2 1 5 4 9 180.00 4.60240 2 ; C2- S1- C5- C4 2 1 5 10 9 180.00 4.60240 2 ; C2- S1- C5- H10 2 3 4 5 9 180.00 19.87400 2 ; C2- N3- C4- C5 2 3 4 9 9 180.00 19.87400 2 ; C2- N3- C4- H9 2 6 8 11 9 180.00 0.00000 2 ; C2- C6- C8- H11 2 6 8 12 9 180.00 0.00000 2 ; C2- C6- C8- H12 2 6 8 13 9 180.00 0.00000 2 ; C2- C6- C8- H13 3 2 6 7 9 180.00 12.02900 2 ; N3- C2- C6- O7 3 2 6 8 9 180.00 12.02900 2 ; N3- C2- C6- C8 3 4 5 10 9 180.00 16.73600 2 ; N3- C4- C5- H10 4 3 2 6 9 180.00 19.87400 2 ; C4- N3- C2- C6 5 1 2 3 9 180.00 4.60240 2 ; C5- S1- C2- N3 5 1 2 6 9 180.00 4.60240 2 ; C5- S1- C2- C6 7 6 8 11 9 0.00 3.47272 1 ; O7- C6- C8- H11 7 6 8 11 9 180.00 0.16736 3 ; O7- C6- C8- H11 7 6 8 12 9 0.00 3.47272 1 ; O7- C6- C8- H12 7 6 8 12 9 180.00 0.16736 3 ; O7- C6- C8- H12 7 6 8 13 9 0.00 3.47272 1 ; O7- C6- C8- H13 7 6 8 13 9 180.00 0.16736 3 ; O7- C6- C8- H13 9 4 5 10 9 180.00 16.73600 2 ; H9- C4- C5- H10 [ dihedrals ] ; impropers ; treated as propers in GROMACS to use correct AMBER analytical function ; i j k l func phase kd pn 1 2 3 6 4 180.00 4.60240 2 ; S1- C2- N3- C6 1 5 10 4 4 180.00 4.60240 2 ; S1- C5- H10- C4 5 9 4 3 4 180.00 4.60240 2 ; C5- H9- C4- N3 8 2 6 7 4 180.00 43.93200 2 ; C8- C2- C6- O7
该ITP文件可分为两段，第一段为基础力场参数[ atomtypes ]，第二段为分子力场参数[ moleculetype ]及后续所有参数块。
[ atomtypes ]第一段是基础力场参数[ atomtypes ]，即原子定义和对应的范德华参数。每行代表一种原子类型，Gromacs标准格式中，第一列为原子类型，第二列为原子序号，第三列为原子质量（单位a.m.u），第四列为原子电荷（单位e），第五列为原子类型，第六列为$\sigma$（单位nm），第七列为$\epsilon$（单位kJ/mol）。从GAFF生成的力场文件中，第一列与GAFF.prm原子类型一致，第二列可忽略，原子质量和原子电荷为0.0，因为分子中的具体原子的质量和电荷会覆盖这两个参数，第五列和第六列则由GAFF.dat文件中的范德华参数$r_0$和$\epsilon$转化而来。这一部分的原子类型只有小分子用到的原子类型，事实上读者也可以做一份完全版的GAFF力场基础参数的ITP文件来替换这段，对于有多种不同小分子的体系减少去除重复原子类型的工作。
[ moleculetype ]第二段则是具体的分子类型，本段由[ moleculetype ]开始直到结束，严格地讲，[ moleculetype ]的内容直到读到下一个[ moleculetype ]或[ system ]为止。ITP文件中没有[ system ]内容块，所以实际上Gromacs程序是在读取TOP文件时#include MOL.itp后读到[ system ]或另一个[ moleculetype ]才停止。
[ moleculetype ]下的分子属性包括原子[ atoms ]，键长[ bonds ]，键角[ angles ]，二面角和异二面角[ dihedrals ]，1,4作用对[ pairs ]。
分子类型[ moleculetype ]只有一行，两列，第一列为分子名，用于区分不同分子，TOP文件的[ ystem ]参数块会使用这个名字；第二列为nrexcl，非键作用排除原子数，这里是3，代表连接再3个原子以内的弱相互作用全部排除，即不考虑1,2相互作用和1,3相互作用，考虑1,4相互作用。
1 2 3 [ moleculetype ] ;name nrexcl MOL 3
[ atoms ] [ bonds ] [ angles ]原子类型[ atoms ]行数与坐标中的原子数一致，原子类型也与坐标中的原子一一对应。每一行有8列，第一列为原子序号，严格从1开始，第二列为GAFF原子类型，第三列为残基编号，严格从1开始，第四列残基名，可与分子名不同，第五列cgnr为电荷组，对大分子的重复单元如残基或碱基有意义，对小分子无意义，第六列偏电荷，单位$e$，第七列为质量，单位$a.m.u$，;后都是注释，本文件注释主要计算各个原子累加电荷是否为分子静电荷，仅作提示用。Gromacs的grompp程序会优先读取此处的电荷和质量，而忽略[ atomtypes ]中的电荷和质量。
1 2 3 [ atoms ] ; nr type resi res atom cgnr charge mass ; qtot bond_type 1 ss 1 UNL S1 1 0.048700 32.06000 ; qtot 0.049
键参数[ bonds ]为分子内的所有键连关系，共5列，第一列和第二列为键连接的原子序号，第三列funct为键长力场函数类型，此处是1代表简谐势，第四列平衡键长，单位nm，第五列为力常数$k_l$，单位$kJ/(mol\cdot\text{\AA}^2)$，分号;后为ACPYPE提供的键连原子名注释，用于人工检查，程序不会读取。
1 2 3 [ bonds ] ; ai aj funct r k 1 2 1 1.7560e-01 1.9414e+05 ; S1 - C2
键角参数[ angles ]为分子内所有的键角关系，共6列，第一列，第二列和第三列为键角连接的三个原子序号，第四列funct为键角力场函数类型，此处是1代表简谐势，第五列平衡键角，单位°，第六列为力常数$k_\Theta$，单位$kJ/(mol\cdot Degree^2)$，分号;后为ACPYPE提供的键角原子名注释，用于人工检查，程序不会读取。
1 2 3 [ angles ] ; ai aj ak funct theta cth 1 2 3 1 1.1451e+02 6.9705e+02 ; S1 - C2 - N3
[ dihedrals ] [ pairs ]二面角和异二面角参数[ dihedrals ]为分子内的二面角关系，共8列，前四列为二面角连接的三个原子序号，第五列funct为二面角力场函数类型，此处是4代表周期函数，第六列为相位，单位°，第七列为势垒高度，单位$kJ$，第八列为周期，分号;后为ACPYPE提供的二面角原子名注释，用于人工检查，程序不会读取。
1 2 3 4 5 6 7 8 9 10 11 [ dihedrals ] ; propers ; for gromacs 4.5 or higher, using funct 9 ; i j k l func phase kd pn 1 2 3 4 9 180.00 19.87400 2 ; S1- C2- N3- C4 ... [ dihedrals ] ; impropers ; treated as propers in GROMACS to use correct AMBER analytical function ; i j k l func phase kd pn 1 2 3 6 4 180.00 4.60240 2 ; S1- C2- N3- C6 ...
1,4相互作用[ pairs ]，共3列，第一列和第二列为成对的两个1,4原子序号，第三列funct为函数类型，此处是1代表使用缩放的非键相互作用，分号;后为ACPYPE提供的1,4作用对原子名注释，用于人工检查，程序不会读取。
根据以往经验，1,4相互作用关系数量应该和二面角参数数量相等，但实际上该分子二面角参数有25个，1,4相互作用关系却只有20个。原因在于AMBER/GAFF力场不计算某些1,4作用，如芳环内的1,4作用，因为对于这样的体系成键项几乎不变，所以可忽略1,4相互作用。
1 2 3 [ pairs ] ; ai aj funct 1 7 1 ; S1 - O7
除此以外还有样例中没提到的平衡过程中常用的位置限制[ position_restraints ]，构象或距离约束[ constraints ]，水模型特有的限制方法[ settles ]，以及用作哑原子的多种[ virtual_sites ]等，可以查阅官方文档 进一步学习。
对Amber程序熟悉的读者可能知道，AmberTools在使用tleap导入小分子力场文件时是可以使用.mol2和.frcmod文件的。其中.frcmod文件是由AmberTools中的parmchk或parmchk2程序生成的力场检查文件。例如在tleap中使用如下命令导入GAFF2小分子力场参数：
1 2 3 4 5 6 7 8 9 ## Load Force Fields source leaprc.gaff2 loadAmberParams MOL_AC.frcmod MOL = loadmol2 MOL_bcc_gaff2.mol2 saveamberparm MOL MOL.prmtop MOL.inpcrd quit
就可以获取.prmtop和.inpcrd文件。mol2文件读者可能比较熟悉，SYBYL的.mol2只存储原子类型，原子坐标，原子电荷，以及键连关系，缺少的键长，键角，二面角参数则在leaprc.gaff2和.frcmod文件中。
我们打开文章第一部分生成的MOL_AC.frcmod文件：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 Remark line goes here MASS BOND ANGLE DIHE cc-cd-ss-cd 2 2.200 180.000 2.000 same as X -c2-ss-X , penalty score=232.0 nc-cd-ss-cd 2 2.200 180.000 2.000 same as X -c2-ss-X , penalty score=232.0 h5-cd-ss-cd 2 2.200 180.000 2.000 same as X -c2-ss-X , penalty score=232.0 c -cd-ss-cd 2 2.200 180.000 2.000 same as X -c2-ss-X , penalty score=232.0 IMPROPER cd-h4-cc-nc 1.1 180.0 2.0 Same as X -X -ca-ha, penalty score= 67.2 (use general term)) c -cc-cd-ss 1.1 180.0 2.0 Using the default value c3-cd-c -o 10.5 180.0 2.0 Using general improper torsional angle X- X- c- o, penalty score= 6.0) h5-nc-cd-ss 1.1 180.0 2.0 Same as X -X -ca-ha, penalty score= 67.5 (use general term)) NONBON
文件中只列出了4个二面角和4个异二面角参数。原因是除这8个成键项，其余的成键项和非键项在前一章描述的GAFF2.dat文件中存在，因此AmberTools的tleap程序使用source leaprc.gaff2命令引入GAFF2参数，即可补完小分子力场参数。
对于我们使用Gromacs格式的ITP文件来说，.frcmod文件还有意义么？答案是有意义的，这里单开一节描述.frcmod文件仍是有意义的，因为.frcmod会告诉软件或读者，这里的8个力场参数是缺失的，是parmchk通过规则猜测 的。比如说：
1 nc-cd-ss-cd 2 2.200 180.000 2.000 same as X -c2-ss-X , penalty score=232.0
二面角nc-cd-ss-cd在GAFF力场参数中是缺失的，因此parmchk查找到X -c2-ss-X这个二面角（S-烯烃C），并给出罚分[enalty score]=232.0。这个罚分相当高，原则上是要人工拟合调整的，不过读者可以基于常识进行判断，我们的噻唑环上的二面角和烯烃的二面角一样不可翻转，因此直接借用烯烃二面角参数不影响噻唑环的稳定性，因此这一组二面角参数是可以直接借鉴使用的。
异二面角同理，比如cd-h4-cc-nc噻唑的H的平面性问题使用了芳环的X -X -ca-ha通用异二面角参数，罚分中等，为67.5；c3-cd-c -o羰基平面性使用了通用羰基异二面角X- X- c- o，罚分最低，只有6.0。
不论是使用Amber格式还是Gromacs格式，.frcmod文件中的补全参数都需要小心判断，特别是高罚分的力场参数。比如本例中的4个二面角参数 X-cd-ss- X，罚分高达232，尽管我们根据化学经验判断对分子构象影响不大，但仍要注意。
本次教程介绍了Gromacs ITP文件格式下的小分子力场的描写方式，读者可根据第一章定制小分子力场(一)：分子力场基础 的AMBER力场形式对应学习。从下一章开始，我们将正式讨论定制化力场参数的拟合原理和拟合过程。