本文主要是介绍分子动力学模拟学习-Gromacs工具链，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
1、总体流程
在gromacs的使用说明中有一个flow chart，比较简略。以下针对一般体系（非蛋白等领域）进行了一些调整，通用性更强。
在做分子动力学模拟时，其复杂性除了以上各种输入输出文件的操作，另一点就是力场参数、原子电荷数据的获取。从本质上讲，MD模拟软件就是读取系统内原子的位置、质量、电荷、成键和非键关系，在不同的环境条件下进行计算，而我们做的就是正确地构建输入文件。
2、几何建模
各种软件都可以，最终目标是输出.pdb。从此处开始贯彻一个原则，即当系统中存在多种分子时，最好对每种分子从建模到生成拓扑文件都单独处理，最后再汇总。
对于几何文件，汇总的方式可以是通过建模软件在盒子里加入所需的分子，或者通过gmx的内置工具如editconf，solvate和insert-molucules进行。最终汇总的pdb（或者gro）文件内，各原子的编号需和汇总的拓扑文件（.top）中的原子编号一致，gmx通过该编号来对应各个原子的空间位置。
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.215n2t中一定要包含所有的原子类型，否则会出现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进行格式转换，不过比较大的分子似乎转换会失败。
这篇关于分子动力学模拟学习-Gromacs工具链的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
