大家好,又见面了,我是你们的朋友全栈君。
与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