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

大家好，又见面了，我是你们的朋友全栈君。

2a3df004b6d10e59e0c1a172448aab6c.png

与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

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

发布者：全栈程序员栈长，转载请注明出处：https://javaforall.cn/141389.html原文链接：https://javaforall.cn