Amber小分子-蛋白复合体分子动力学模拟

以前经常用GROMACS进行分子动力学模拟，后来试了一下Amber后发现，在我当前配置的GPU资源上，果然还是Amber更快一些，GROMACS太吃CPU资源了，导致多卡GPU没法高效的跑多任务，于是开始摇摇晃晃的Amber分子动力学模拟。事实证明Amber除了快以外，在AmberTools的加持下，还更简单上手。以下只是一个我自己摸索的能够跑通"小分子-蛋白复合体"的Amber分子模拟过程。详细原理后面有空再补（挖坑ing）

1. 小分子、蛋白预处理

1.1 小分子预处理

被别人准备的小分子害的跑了很多次，找了很多原因都没跑通，最终发现大部分原因是初始小分子结构不准确。小分子的预处理的唯一准则是得到一个质子化后结构正确的三维坐标即可。

晶体结构中的小分子

当你的模拟体系的小分子来源于RCSB PDB数据库的复合物中的小分子时，可以直接从PDB数据库中获取小分子，推荐格式为MOL2。
不要随便加氢后就直接拿去跑模拟，无数教训告诉我，质子化以后为了稳妥起见，还是做一下能量最小化比较好。

对接复合物种的小分子
当你的模拟体系的小分子需要先对接，对接后的小分子构象作为模拟的初始结构时，先做小分子与蛋白的分子对接，挑选最佳的对接构象作为模拟对象，将对接后的小分子保存为MOL2格式。

1.2 蛋白预处理

蛋白的预处理比较简单，只需要将从RCSB PDB数据库中获取的蛋白部分进行质子化即可，在剔除非蛋白部分的原子后，可以使用AmberTools中的pdb4amber完成蛋白质的加氢。具体的加氢程序是由杜克大学的Richardson实验室开发的Reduce程序完成。

pdb4amber -i test.pdb -o rec.pdbpdb -y --reduce

-i：输入蛋白质的pdb格式文件
-y：该参数是指去除所有氢原
--reduce ：运行Reduce来加氢
之所以用了-y和–reduce选项，是因为pdb4amber会检查每个组氨酸（HIE、HIP、HID）的类型，如果用其他软件加氢则不一定是amber所支持的原子类型，所以还是用pdb4amber一劳永逸吧。

如果你的蛋白质需要做一些修改，例如某些氨基酸位点的突变，则需要做突变后结构的能量最小化，否则初始结构中某些原子比较靠近，导致系统能量过高从而使得模拟崩溃。
我常用Schrödinger的Maestro中的Protein Prepare Wizard做蛋白的常态化处理，可以一次性解决质子化、侧链缺失、loop缺失、能量最小化等问题。

2. 模拟体系预处理

2.1 小分子预处理

利用Acpype为小分子生成拓扑参数文件。Acpype基于Antechamber，能够生成CNS/XPLOR、GROMACS、CHARMM和Amber等MD软件的拓扑参数文件。对于我们这种既使用GROMACS又使用Amber的人来说十分友好。

acpype -i test.mol2 -a gaff2 -c bcc -n 0

-i：输入文件，建议为mol2格式，其他格式将通过openbabel转换
-a：小分子力场参数
-c：电荷类型
-n：电荷数
之所以选择gaff2力场，是因为该力场涵盖了大多数的元素。bcc电荷是一种半经验的拟合静电势电荷，计算速度极快，操作简单，包含在Antechamber中，虽然比起RESP电荷而言精度较低，但无需依赖额外的程序。等有空再研究下精度更高的RESP电荷如何简单获取。

2.2 小分子-蛋白复合物体系

利用Leap/Tleap程序完成“小分子-蛋白”复合物amber参数体系的构建。

tleap -f tleap.leaprc

tleap.leaprc的参数如下所示：

source leaprc.protein.ff14SB #蛋白ff14SB力场
source leaprc.water.tip3p #水的tip3p力场
source leaprc.gaff2 #小分子的gaff2力场
set default pbradii bondi #原子半径
loadamberparams frcmod.ionsjc_tip3p #溶剂的力场参数文件
loadamberparams test.acpype/test_AC.frcmod #修改后的小分子力场参数文件
rec = loadpdb rec.pdb #加载蛋白文件
lig = loadmol2 test.acpype/test_bcc_gaff2.mol2 #加载小分子文件
complex = combine{rec,lig} #行程复合体
addions complex Na+ 0 #电荷平衡
addions complex Cl- 0 #电荷平衡
solvateoct complex TIP3PBOX 10 #形成一个TIP3P水模型的盒子，半径是10埃米
saveamberparm complex complex_ions.prmtop complex_ions.inpcrd #保存amber格式的拓扑和坐标文件
quit

生成 complex complex_ions.prmtop 和complex_ions.inpcrd 。

2.4 Amber模拟可以直接跳过

其实当处理到这一步时，只需要再利用ParmEd转换成GROMAC的输入格式即能立即进入模拟体系。

import parmed as pmd
amber=pmd.load_file('complex_ions.prmtop','complex_ions.inpcrd')
amber.save('complex_ions.top')
amber.save('complex_ions.gro')

3. Amber MD

首先需要各个阶段的参数文件，min1.in，min2.in，heat.in，press.in，eq.in，md.in，详见4. Amber参数文件。
主要的模拟步骤如下脚本所示：

#expose a single GPU to the pmemd.cuda executable
export CUDA_VISIBLE_DEVICES="6"#1. Minimize the system with restraints just on the backbone of the molecule
pmemd.cuda -O -i min1.in -o complex_ions_min1.out -p complex_ions.prmtop -c complex_ions.inpcrd -r complex_ions_min1.rst -ref complex_ions.inpcrd#2. Relax the system with no restraints
pmemd.cuda -O -i min2.in -o complex_ions_min2.out -p complex_ions.prmtop -c complex_ions_min1.rst -r complex_ions_min2.rst#3. Heat up the system from 100 K under constant volume
pmemd.cuda -O -i heat.in -o complex_ions_heat.out -p complex_ions.prmtop -c complex_ions_min2.rst -r complex_ions_heat.rst -x complex_ions_heat.nc -ref complex_ions_min2.rst#4. Relax the system at a constant pressure
pmemd.cuda -O -i press.in -o complex_ions_press.out -p complex_ions.prmtop -c complex_ions_heat.rst -r complex_ions_press.rst -x complex_ions_press.nc -ref complex_ions_min2.rst#5. eq, 10ns
pmemd.cuda -O -i eq.in -o complex_ions_eq.out -p complex_ions.prmtop -c complex_ions_press.rst -r complex_ions_eq.rst -x complex_ions_eq.nc#6. md, 100ns
pmemd.cuda -O -i md.in -o complex_ions_md.out -p complex_ions.prmtop -c complex_ions_eq.rst -r complex_ions_md.rst -x complex_ions_md.nc

4. Amber参数文件

min1.in

&cntrl
imin=1,
maxcyc=10000,
ncyc=5000,
ntb=1,
ntr=1,
restraintmask=‘@CA,N,C’,
restraint_wt=10,
cut=8.0
/
END

min2.in

&cntrl
imin=1,
maxcyc=100000,
ncyc=5000,
ntb=1,
ntc=1,
ntf=1,
ntpr=10,
cut=8.0
/
END

heat.in

explicit solvent initial heating: 50ps
&cntrl
imin=0,
irest=0,
nstlim=25000, dt=0.002,
ntc=2, ntf=2, ntx=1,
cut=8.0, ntb=1,
ntpr=500, ntwx=500,
ntt=3, gamma_ln=2.0,
tempi=0.0, temp0=300.0, ig=-1,
ntr=1,
restraintmask=‘:1-1104’,
restraint_wt=2.0,
iwrap=1
nmropt=1
/
&wt TYPE=‘TEMP0’, ISTEP1=0, ISTEP2=25000,
VALUE1=0.0, VALUE2=300.0, /
&wt TYPE = ‘END’ /
END

press.in

explicit solvent density: 50ps
&cntrl
imin=0,
irest=1,
ntx=5,
nstlim=25000, dt=0.002,
ntc=2, ntf=2,
cut=8.0, ntb=2, ntp=1, taup=2.0,
ntpr=500, ntwx=500,
ntt=3, gamma_ln=2.0,
temp0=300.0, ig=-1,
ntr=0,
/
END

eq.in（10 ns）

&cntrl
imin=0, irest=1, ntx=5,
nstlim=5000000, dt=0.002,
ntc=2, ntf=2,
cut=10.0, ntb=2, ntp=1, taup=2.0,
ntpr=500, ntwx=500, ntwr=5000,
ntt=3, gamma_ln=2.0,
temp0=300.0,
/
END

md.in （100 ns）

explicit solvent production run: 100ns
&cntrl
imin=0,
irest=1,
ntx=5,
nstlim=50000000, dt=0.002,
ntc=2, ntf=2,
cut=8.0, ntb=2, ntp=1, taup=2.0,
ntpr=5000, ntwx=5000, ntwr=50000,
ntt=3, gamma_ln=2.0,
temp0=300.0, ig=-1,
iwrap=1
/
END