分子动力学模拟学习-Gromacs工具链-编程知识

1、总体流程

在gromacs的使用说明中有一个flow chart，比较简略。以下针对一般体系（非蛋白等领域）进行了一些调整，通用性更强。

在做分子动力学模拟时，其复杂性除了以上各种输入输出文件的操作，另一点就是力场参数、原子电荷数据的获取。从本质上讲，MD模拟软件就是读取系统内原子的位置、质量、电荷、成键和非键关系，在不同的环境条件下进行计算，而我们做的就是正确地构建输入文件。

2、几何建模

各种软件都可以，最终目标是输出.pdb。从此处开始贯彻一个原则，即当系统中存在多种分子时，最好对每种分子从建模到生成拓扑文件都单独处理，最后再汇总。

对于几何文件，汇总的方式可以是通过建模软件在盒子里加入所需的分子，或者通过gmx的内置工具如editconf，solvate和insert-molucules进行。最终汇总的pdb（或者gro）文件内，原子的编号需和汇总的拓扑文件（.top）中的原子编号一致，gmx通过该编号来对应各个原子的空间位置。

04.28补充：上述的“原子编号一致”说法不严谨。在几何文件（pdb或gro）中的原子编号是从1到n，所有分子统一编号，例如编号1~a是分子甲，a+1~b是分子乙，b+1~n是分子丙。而在top文件内，每一段落的原子编号都是从1开始的（具体的“段落”的含义见下方第3部分）。在top内，各分子的“段落”排序可以随意，但在最后的[molecules]中，分子的顺序要与几何文件内一致。此外，在[atoms]内部的原子顺序理应也要与几何文件中，对应分子的原子顺序相同。在调用grompp时，gmx会检查top中的[atoms]的第五列和几何文件的第3列是否一致，如果不一致会产生warning。

3、生成拓扑文件

3.1、x2top用于小分子或晶体

在gmx中和拓扑文件相关的文件类型有很多，有.top，.itp，.rtp，.atp。如果采用x2top工具生成拓扑文件，我们只需要关注.top和.itp类型，另外两个是用于另一个工具pdb2gmx生成拓扑文件。pdb2gmx本身适合于蛋白质、核酸等生物大分子，有许多预置的特性，不适合小分子或者晶体，但经过修改可以用于高分子聚合物，但终究是麻烦。

以下预设的情况是对于某个分子，我们并不能使用gmx的自带力场参数或原子电荷（指partial charge，以下的“原子电荷”都是这个含义），需要自行填入文献或实验数据。

.top文件有较严格的组织方式，几个部分必须按照顺序排列。一种典型的格式为：

; 行首加分号为注释语句
; gmx读取文件的逻辑是各列以空格为分隔符，连续的多个空格（不论有多少个）都会被视为是一个分隔符。所以下面标了（不重要）的项，也需要填入至少一个字符。; 该字段分别为非键作用函数类型（LJ或buckingham），混合规则，是否考虑分子内1-4非键作用（pair），1-4非键作用的范德华作用因子，1-4非键作用的库伦作用因子
; 关于非键作用函数类型和混合规则的说明见https://manual.gromacs.org/current/reference-manual/topologies/parameter-files.html
; 关于分子内1-4非键作用见http://sobereva.com/4
[ defaults ]
; nbfunc	comb-rule	gen-pairs	fudgeLJ	fudgeQQ
1		3		yes		0.5		0.5; 此处要填入系统内所有原子类型
; 各列分别是原子类型，该类型原子数或键类型（这列似乎不重要，填什么都可以，例如全0），原子质量（不重要，因为质量一般会在[atoms]中写），电荷（不重要，因为电荷一般会在[atoms]中写，并且同种原子电荷未必相等），
; 粒子类型（见https://manual.gromacs.org/current/reference-manual/topologies/particle-type.html#tab-ptype） ，范德华力的sigma、epsilon参数（根据nbfunc，也可能是C12、C6或者buckingham的ABC）
[ atomtypes ]
; atom_type  bond_type    mass    charge   ptype          sigma      epsilon; 这里用来写不同原子之间的非键作用参数，对于没写的原子对则直接使用混合规则
[ nonbond_params]
; i    j 	func          sigma      epsilon; 从这里开始直到[dihedrals]，每种分子都需要写一组
; 各列含义为分子名称（可自定义），非键排除，排除相邻nrexcl个键的非键相互作用（引自https://blog.csdn.net/CocoCream/article/details/123769268） 全称可能是number of exclusion
[ moleculetype ]
; Name            nrexcl; 分子内各原子信息，各列分别是原子编号，原子类型，残基编号（不重要），残基名（不重要），原子符号（C、H、O等，最好与几何文件内的符号一致），原子电荷组（用于批量化操作库仑力截断半径等），
; 原子电荷，原子质量，用于自由能计算的参数（见https://manual.gromacs.org/current/reference-manual/topologies/topology-file-formats.html#topologies-for-free-energy-calculations）
[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  typeB    chargeB      massB; 成键参数伸缩项，各列分别是原子i的类型，原子j的类型，函数类型（见https://manual.gromacs.org/current/reference-manual/topologies/topology-file-formats.html#tab-topfile2）
; 平衡时的键长，参数1，参数2，参数3（不一定会都用上，例如funct=1时只需要参数c1，其他两个没有意义）
[ bonds ]
;  ai    aj funct            c0            c1            c2            c3; 分子内的原子对，如果前面的gen-pairs选择了yes就会从这里读取原子对，如果这里写了参数则还会读取这里的参数
[ pairs ]
;  ai    aj funct            c0            c1            c2            c3; 成键参数弯曲项，类似伸缩项
[ angles ]
;  ai    aj    ak funct            c0            c1            c2            c3; 成键参数扭转项，类似伸缩项
[ dihedrals ]
;  ai    aj    ak    al funct            c0            c1            c2            c3            c4            c5; 整个系统的名称
[ system ]
; Name
zro2_bigger; 指定系统内各分子的数量
[ molecules ]
; Compound        #mols
zro2_bigger         1

除了全部列在此处以外，还可以通过#include语句来引用.itp文件（.itp的用处就在于此），这样可以实现对同一个分子在不同系统中的引用。不过引用时需注意将[atomtypes]字段复制到最终的.top中，因为atomtypes字段只允许出现一次。

更多关于.top的说明：

http://bbs.keinsci.com/forum.php?mod=viewthread&tid=28201&highlight=%CD%D8%C6%CB

http://bbs.keinsci.com/forum.php?mod=viewthread&tid=27657&highlight=%CD%D8%C6%CB

http://bbs.keinsci.com/thread-14723-1-13.html

http://bbs.keinsci.com/forum.php?mod=viewthread&tid=29033&extra=&highlight=top&page=1

http://bbs.keinsci.com/forum.php?mod=viewthread&tid=37240&highlight=top

http://bbs.keinsci.com/thread-19761-1-1.html

http://bbs.keinsci.com/forum.php?mod=viewthread&tid=32292&highlight=top

http://bbs.keinsci.com/forum.php?mod=viewthread&tid=27464&highlight=top

http://bbs.keinsci.com/forum.php?mod=viewthread&tid=39032&highlight=atomtypes

http://bbs.keinsci.com/forum.php?mod=viewthread&tid=19125&highlight=nonbond

下面是如何用x2top来生成上述的文件。

上面的.top文件中，可以从几何文件获取的只有原子编号，怎样将原子的位置信息转化为成键信息（两个原子是否成键），以及如何指定原子的质量和电荷呢？这就要配合.n2t文件一起使用。n2t即name to type，意义是将pdb中的各个原子映射到gmx中定义的原子类型。在原子类型中，带有该原子与谁成键，以及该类型原子的质量和电荷信息（其中电荷信息没用，因为同一类型的原子电荷未必相同）。典型的n2t文件：

; 原子符号    原子类型    电荷    质量    成键数量    成键原子1 键长1    成键原子2 键长2 ……
ZR   ZRO2_ZRO8   0     91.224    8    O  0.2195   O 0.2195   O 0.2195   O 0.2195   O 0.2195   O  0.2195   O 0.215    O 0.24
ZR   ZRO2_ZRO6   0     91.224    6    O  0.2195   O 0.2195   O 0.2195   O 0.2195   O 0.2195           O  0.2195
O    ZRO2_OZR4   0     15.9994   4    ZR 0.2195  ZR 0.2195  ZR 0.2195  ZR 0.2195
O    ZRO2_OZRH   0     15.9994   2    ZR 0.215    H 0.11
H    ZRO2_HO     0     1.008     1     O 0.11
O    ZRO2_OZR3   0     15.9994   3    ZR 0.215   ZR 0.215   ZR 0.215

n2t中一定要包含所有的原子类型，否则会出现not found错误。如果已经定义了所有类型，仍然有原子not found，则需要检查键长是否合适，是否与几何对应。不一定完全相同，但一定要接近，可以通过建模软件进行测量（gmx的长度单位是nm）。

关于获取原子电荷的方法

原子电荷应与其力场适配，通常研究力场的文献内会写明拟合力场参数时，所使用的原子电荷。如果找不到对应电荷，可能需要通过第一性原理计算软件进行计算。可以计算原子电荷的路线包括MS计算QEq或者通过dmol3或castep计算mulliken或hirshfeld电荷、gaussian和antechamber拟合RESP电荷、gaussian和multiwfn拟合RESP电荷、cp2k计算RESP电荷等。

关于获取力场参数的方法

力场参数一般通过#include对应的力场的.itp文件来进行，如果知道参数，也可以根据gmx自带的力场.itp文件手写。

如果对成键作用没有要求（例如研究固液界面，固体发生的变化不重要），也可以在x2top中直接写入成键参数，见GROMACS教程：创建周期性体系的拓扑文件：以石墨烯为例|Jerkwin

其他Interatomic Potentials - LAMMPS Tube

单位转换https://www.colby.edu/chemistry/PChem/Hartree.html

3.2 高分子的拓扑

高分子拓扑的生成可以魔改pdb2gmx，或用antechamber+actype。更简单逃课的方法是利用现成工具MD模拟中力场文件生成工具 - 知乎、自己开发的力场文件生成工具 - 知乎

经过尝试，针对聚合物，比较好的路线是在Marvin JS中画好分子式，复制至PolyParGen生成.gro和.itp，支持oplsaa和amber力场。自行上传需要cml格式文件，可以通过OPENBABEL - Chemical file format converter进行格式转换，不过比较大的分子似乎转换会失败。