Gromacs联用GAFF力场处理水溶剂下小分子动力学
与GROMACS偏重生物大分子模拟的力场不同，AMBER支持很多方便处理有机小分子的力场，如GAFF力场，简单而又有不错的精度，适合处理有机小分子；这里将介绍用Gaussian计算RESP电荷，交由Amber生成GAFF力场下的拓扑文件，最后用GROMACS模拟的过程。
软件版本：AmberTools18，Gromacs 2019-beta-1 GPU版
1.使用G16计算RESP电荷
首先应先在
b3lyp/6-31g(d) scrf=(smd,solvent=water) empiricaldispersion=gd3下进行优化，得到大致合理的结构，再进行静电势计算
静电势计算的GJF文件如下
%nprocshared=24%mem=40GB%chk=gc5ay.chk#b3lyp/6-31g(d)scrf=(smd,solvent=water) eom=connectivityempiricaldispersion=gd3pop=mk iop(6/33=2,6/42=6,6/50=1)Symmetric GC5AY5 1*中略*gc5ay.gesp#末行保留两行空行
最后只需要得到的gc5ay.gesp文件即可。
2.生成拓扑文件
使用上一步得到的gc5ay.gesp，运行
antechamber -i gc5ay.gesp -fi gesp -o gc5ay.mol2 -fo mol2 -pf y -c resp我们只需要gc5ay.mol2输出文件进行下一步计算, 它包含了构型以及RESP电荷。
使用parmchk2检查GAFF参数并生成缺失参数文件
使用上一步得到的gc5ay.mol2文件, 运行parmchk2命令
parmchk2 -i gc5ay.mol2 -f mol2 -o gc5ay.frcmodparmchk2检查输入分子构型中GAFF的缺失参数, 并生成相应的补充参数文件gc5ay.frcmod.
使用tleap生成AMBER参数文件及坐标文件
命令行输入tleap，然后输入以下内容（针对AmberTools2018版本）
source oldff/leaprc.ff14SBsource leaprc.gaffloadamberparams gc5ay.frcmodgc5ay=loadmol2 gc5ay.mol2source leaprc.water.tip3psolvatebox gc5ay TIP3PBOX 10addions gc5ay Cl- 0check gc5aysaveamberparm gc5ay gc5ay.prmtop gc5ay.inpcrdquit
此命令包含了载入Amber14SB力场（针对三点水）/GAFF力场，创造10*10*10溶剂化盒子，补充Cl-反离子中性化，检查体系，生成拓扑文件。
这样就拿到了分子的AMBER参数文件gc5ay.prmtop, 结构文件gc5ay.inpcrd.
如果要处理的是复合物，计算完ESP电荷后，应该拆开gesp文件，分别计算RESP电荷、用antechamber和parmchk2处理，在tleap中分别加载再用combine指令合并，可以用pdb文件保存检视合并状况
source oldff/leaprc.ff14SBsource leaprc.gaffloadamberparams cortisol-cp.frcmodloadamberparams gc5ay-cp.frcmodhost=loadmol2 gc5ay-cp.mol2guest=loadmol2 cortisol-cp.mol2complex = combine {host guest}source leaprc.water.tip3psolvatebox complex TIP3PBOX 10addions complex Cl- 0check complexsaveamberparm complex complex.prmtop complex.inpcrdsavepdb complex complex.pdbquit
然后运行第三方Python脚本acpype将amber格式拓扑文件转换为gromacs支援的格式
acpype -p gc5ay.prmtop -x gc5ay.inpcrd -d这样就得到了GROMACS支持的gc5ay_GMX.gro, gc5ay_GMX.top, em.mdp, md.mdp等文件. 一般我们只需要前面两个文件。
如果想将.top文件进行处理生成.itp文件，以便在蛋白质的拓扑文件中包含, 可以除去表头, 改动原子类型, 再除去后面的附加信息。
注意的是这里要手动修改top文件，将第七行对下来的[ atomtypes ]下的Cl-修改为大写的CL-，以及最底下描述离子信息的[ atom ]下的IM改为CL-，这才和底下的离子信息对得上，否则待会gromacs运行会报错”atom type XX not found“
3. Gromacs运行作业
执行能量最小化
gmx grompp -f em.mdp -c gc5ay_GMX.gro -p gc5ay_GMX.top -o em.tprgmx mdrun -v -deffnm em
脚本文件em.mdp
define = -DFLEXIBLEintegrator = cgnsteps = 500emtol = 100.0emstep = 0.01;nstxout = 50nstlog = 50nstenergy = 50;pbc = xyzcutoff-scheme = Verletcoulombtype = PMErcoulomb = 1.0vdwtype = Cut-offrvdw = 1.0DispCorr = EnerPres;constraints = none
执行10ns模拟
gmx grompp -f md.mdp -c em.gro -p gc5ay_GMX.top -o md.tprgmx mdrun -v -deffnm md
供参考的脚本文件md.mdp
define =integrator = mddt = 0.002 ; psnsteps = 5000000 ; 10nscomm-grps = systemenergygrps = ;nstxout = 0nstvout = 0nstfout = 0nstlog = 5000nstenergy = 1000nstxout-compressed = 1000compressed-x-grps = system;pbc = xyzcutoff-scheme = Verletcoulombtype = PMErcoulomb = 1.0vdwtype = cut-offrvdw = 1.0DispCorr = EnerPres;Tcoupl = V-rescaletau_t = 0.2tc_grps = systemref_t = 298.15 ;Pcoupl = parrinello-rahmanpcoupltype = isotropictau_p = 2.0ref_p = 1.0compressibility = 4.5e-5;freezegrps =freezedim =constraints = hbonds
在E5 2678 v3+RTX 2060 平台上花费了14分钟完成模拟
4. 分析作业
将分子置于中心并消除旋转便于观察，选择考察对象时可以选择others排除掉水溶剂，便于观察也节省空间
gmx trjconv -f md.xtc -s md.tpr -o cent.xtc -center -pbc molgmx trjconv -f cent.xtc -s md.tpr -fit rot+trans -o fit.xtcgmx trjconv -f md.gro -s md.tpr -o cent.gro -center -pbc molgmx trjconv -f cent.gro -s md.tpr -fit rot+trans -o fit.gro
当然分析作业远远不止观察分子运动那么简单，还可以调用各种脚本考察希望的量，这里暂不讨论。
