分子模拟软件amber_Gromacs联用GAFF力场处理水溶剂下小分子动力学

与GROMACS偏重生物大分子模拟的力场不同，AMBER支持很多方便处理有机小分子的力场（详见http://sobereva.com/115），如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) geom=connectivity
empiricaldispersion=gd3 pop=mk iop(6/33=2,6/42=6,6/50=1)
Symmetric GC5AY
5 1
*中略*
gc5ay.gesp

末行保留两行空行

这里使用b3lyp泛函代替过时的HF方法。以往考虑HF方法计算水溶剂下静电势是因为HF没有考虑电子相关作用，这会导致高估键的极性。而水溶剂效应本身也会极化溶质的电荷分布使键的极性增加，因此用HF在气相下拟合RESP电荷，可以等效地反映出水溶剂的这个效应。但这种方法并不优雅（难道高估的极性就恰好和水溶剂效应一样吗？），这里采用b3lyp加隐式水溶剂模型计算。

最后只需要得到的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.frcmod

parmchk2检查输入分子构型中GAFF的缺失参数, 并生成相应的补充参数文件gc5ay.frcmod.

使用tleap生成AMBER参数文件及坐标文件

命令行输入tleap，然后输入以下内容（针对AmberTools2018版本）
source oldff/leaprc.ff14SB
source leaprc.gaff
loadamberparams gc5ay.frcmod
gc5ay=loadmol2 gc5ay.mol2
source leaprc.water.tip3p
solvatebox gc5ay TIP3PBOX 10
addions gc5ay Cl- 0
check gc5ay
saveamberparm gc5ay gc5ay.prmtop gc5ay.inpcrd
quit

此命令包含了载入Amber14SB力场（针对三点水）/GAFF力场，创造10*10*10溶剂化盒子，补充Cl-反离子中性化，检查体系，生成拓扑文件

这样就拿到了分子的AMBER参数文件gc5ay.prmtop, 结构文件gc5ay.inpcrd.

如果要处理的是复合物，计算完ESP电荷后，应该拆开gesp文件，分别计算RESP电荷、用antechamber和parmchk2处理，在tleap中分别加载再用combine指令合并，可以用pdb文件保存检视合并状况
source oldff/leaprc.ff14SB
source leaprc.gaff
loadamberparams cortisol-cp.frcmod
loadamberparams gc5ay-cp.frcmod
host=loadmol2 gc5ay-cp.mol2
guest=loadmol2 cortisol-cp.mol2
complex = combine {host guest}
source leaprc.water.tip3p
solvatebox complex TIP3PBOX 10
addions complex Cl- 0
check complex
saveamberparm complex complex.prmtop complex.inpcrd
savepdb complex complex.pdb
quit

然后运行第三方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.tpr

gmx mdrun -v -deffnm em

脚本文件em.mdp

define = -DFLEXIBLE
integrator = cg
nsteps = 500
emtol = 100.0
emstep = 0.01
;
nstxout = 50
nstlog = 50
nstenergy = 50
;
pbc = xyz
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb = 1.0
vdwtype = Cut-off
rvdw = 1.0
DispCorr = EnerPres
;
constraints = none

执行10ns模拟

gmx grompp -f md.mdp -c em.gro -p gc5ay_GMX.top -o md.tpr

gmx mdrun -v -deffnm md

供参考的脚本文件md.mdp
define =
integrator = md
dt = 0.002 ; ps
nsteps = 5000000 ; 10ns
comm-grps = system
energygrps =
;
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 5000
nstenergy = 1000
nstxout-compressed = 1000
compressed-x-grps = system
;
pbc = xyz
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb = 1.0
vdwtype = cut-off
rvdw = 1.0
DispCorr = EnerPres
;
Tcoupl = V-rescale
tau_t = 0.2
tc_grps = system
ref_t = 298.15
;
Pcoupl = parrinello-rahman
pcoupltype = isotropic
tau_p = 2.0
ref_p = 1.0
compressibility = 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 mol

gmx trjconv -f cent.xtc -s md.tpr -fit rot+trans -o fit.xtc

gmx trjconv -f md.gro -s md.tpr -o cent.gro -center -pbc mol

gmx trjconv -f cent.gro -s md.tpr -fit rot+trans -o fit.gro

当然分析作业远远不止观察分子运动那么简单，还可以调用各种脚本考察希望的量，这里暂不讨论。