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进行格式转换,不过比较大的分子似乎转换会失败。