Amber进行分子动力学模拟以及计算mmpbsa

1.文件处理

• 将处理后的共晶结构或者对接结果导入，以蛋白与分子相互作用分子
  • 经过对接软件处理后的结果
  • 以蛋白为pdb文件，分子为sdf文件为例

2.蛋白与分子处理

(1) 前处理

1. 使分子带电荷

antechamber -i L.sdf -fi sdf -o L.mol2 -fo mol2 -c bcc -s 2
parmchk2 -i L.mol2 -f mol2 -o L.frcmod

2. 蛋白重加氢

pdb4amber -i R1.pdb -o R_noH.pdb -y --dry
reduce R_noH.pdb > R_H.pdb
pdb4amber -i R_H.pdb -o R.pdb

(2) 生成crd与prm文件

1. tleap.in 文件

source leaprc.gaff
source leaprc.water.tip3p
loadAmberParams frcmod.ionsjc_tip3p 
L = loadmol2 L.mol2 
loadamberparams L.frcmod
saveoff L L.lib
R = loadpdb R.pdb
saveamberparm L L.prmtop L.inpcrd
saveamberparm R R.prmtop R.inpcrd
complex=combine {R L}
saveamberparm complex complex.prmtop complex.inpcrd
charge complex
addions complex Na+ 0 
addions complex Cl- 0
solvatebox complex TIP3PBOX 10.0
saveamberparm complex complex_solvated.prmtop complex_solvated.inpcrd
quit

2. 运行命令

tleap -s -f tleap.in

3.分子动力学模拟

(1)能量最小化

• 分子模拟的初始状态中，体系中原子间可能相距较近使得体系能量过高造成模拟崩溃或原子移步过大移出盒子造成原子丢失，因此一般均要进行能量最小化。此过程将会使得原子排开，减少原子重叠造成的体系能量过高；并限制原子在一个时间内的最大位移。如果输入的是经过处理的共晶构象，也可避免此步骤。
• 约束主链最小化 : min.in文件

 minimise
 &cntrl
  imin=1,maxcyc=1000,ncyc=500, 
  cut=8.0,ntb=1,
  ntc=2,ntf=2,
  ntpr=10,
  ntr=1, restraintmask=':1-110', restraint_wt=2.0
 /

• imin=1: 表示进行能量最小化
  • maxcyc=1000：最小化的最大循环数
  • ncyc=500：在500个周期后最小化方法从最陡下降法转换为共轭梯度法
• cut=8.0: 以埃为单位的非键截断距离(对于PME而言, 表示直接空间加和的截断。不使用低于8.0的值；较高的数字略微提高精度, 但大大增加计算成本)
• ntb=1，定容的周期边界
• ntc=2，不计算震动约束键的力
• ntf=2，允许震动来约束所有和氢键相连的键
• ntpr=10，默认值50，能量和温度等信息写入mdout和mdinfo文件的步长
• ntr=1，约束原子参数(0不约束，>0约束): 由restraintmask限定，straint_wt为力限制。
  • restraintmask=‘:1-110’，约束原子字符
  • restraint_wt=2，位置约束的重量，单位kcal/mol-A2
• 执行命令

pmemd.cuda -O -i min.in -o min.out -p complex_solvated.prmtop -c complex_solvated.inpcrd -r min.rst -ref complex_solvated.inpcrd

(2)体系加热

• heat.in文件

heat 50ps
 &cntrl
  imin=0,irest=0,ntx=1,
  nstlim=25000,dt=0.002,
  ntc=2,ntf=2,
  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-435',
  restraint_wt=10.0,
  nmropt=1
 /
 &wt TYPE='TEMP0', istep1=0, istep2=100000,
  value1=0.1, value2=300.0/
 &wt TYPE='END' /

• imin=0：不进行最小化；irest=0：run as a new MD；ntx=1：无初始速度信息
• nstlim=50000：跑MD的步数；dt=0.002：时间步长（ps）
• ntc=2：允许震动来约束所有和氢键相连的键；ntf=2：不计算震动约束键的力
• cut=8.0：以埃为单位的非键截断距离；ntb=1：定容的周期边界
• ntpr=500：每500步打印一次输出；ntwx=500：每500步写一次amber轨迹文件
• ntt=3：Langevin恒温器温度控制； gamma_ln=2.0：Langevin恒温器碰撞频率
• tempi=0.0：初始温度；temp0=300.0：最终温度；ig=-1：为伪随机数生成程序随机化种子
• ntr=1，约束原子参数(0不约束，>0约束): 由restraintmask限定，straint_wt为力限制。
  • restraintmask=‘:1-110’：约束原子字符
  • restraint_wt=10：位置约束的重量，单位kcal/mol-A2
• nmropt=1：核磁共振限制和读取重量改
• TYPE=‘TEMP0’, istep1=0, istep2=100000, value1=0.1, value2=300.0：允许恒温器在整个模拟过程中改变目标温度。一共100000步，温度将从0K增加到300K。
  因为这个模拟非常简短，所以NTPR和NTWX的设置非常低。使用这些设置进行长时间的MD模拟将产生非常大的输出和轨迹文件，比普通的MD要慢。对于真正的跑MD，你需要增加NTPR和NTWX，不然时间太长
• 执行命令

pmemd.cuda -O -i heat.in -p complex_solvated.prmtop -c min.rst -o heat.out -r heat.rst -ref min.rst -x heat.crd

(3)均匀密度

• density.in文件

density 50ps
 &cntrl
  imin=0,irest=1,ntx=5,
  nstlim=50000,dt=0.002,
  ntc=2,ntf=2,
  cut=4.0, ntb=2, ntp=1, taup=1.0,
  ntpr=500, ntwx=500,
  ntt=3, gamma_ln=2.0,
  temp0=300.0, ig=-1,
  ntr=1, restraintmask=':1-110',
  restraint_wt=2.0,
 /

• imin=0：不进行最小化；irest=1：重启之前的MD（意味着inpcrd文件提供初始原子速度）；ntx=5：从未格式化的inpcrd坐标文件中读取坐标和速度
• nstlim=50000：跑MD的步数；dt=0.002：时间步长（ps）
• ntc=2：允许震动来约束所有和氢键相连的键；ntf=2：不计算震动约束键的力
• cut=4.0：以埃为单位的非键截断距离，ntb=2：定容的周期边界,
• ntp=1, taup=1.0：压力松弛时间 1.0ps
• ntpr=500：每500步打印一次输出，ntwx=500：每500步写一次amber轨迹文件
• ntt=3：Langevin恒温器温度控制,，gamma_ln=2.0：Langevin恒温器碰撞频率
• temp0=300.0：最终温度；ig=-1：为伪随机数生成程序随机化种子
• ntr=1，约束原子参数(0不约束，>0约束): 由restraintmask限定，straint_wt为力限制。
  • restraintmask=‘:1-110’：约束原子字符
  • restraint_wt=2：位置约束的重量，单位kcal/mol-A2
• 执行命令

pmemd.cuda -O -i density.in -p complex_solvated.prmtop -c heat.rst -o density.out -r density.rst -ref heat.rst -x density.crd

(4)全局平衡

• eq.in文件

eq 10ns
&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

• 执行命令

pmemd.cuda -O -i eq.in -p complex_solvated.prmtop -c density.rst -o eq.out -r eq.rst -ref density.rst -x eq.crd

(5)分子动力学模拟

• md.in文件

md 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,
  ntt=3, gamma_ln=2.0,
  temp0=300.0, ig=-1,
  iwrap=1
 /

• iwrap=1：处理镜像
• 执行命令

pmemd.cuda -O -i md.in -p complex_solvated.prmtop -c eq.rst -o md.out -r md.rst -ref eq.rst -x md.crd

4.mmpbsa计算

(1) 新建mmpbsa.in文件，并写入

&general
   startframe=2000, endframe=5000, interval=50,
   keep_files=0, debug_printlevel=2
/
&gb
   igb=2, saltcon=0.1
/
&decomp
   idecomp=1, print_res='1-110;111', csv_format=0,
/

• startframe=2000, endframe=5000, interval=50：从2000帧开始到5000帧结束，每隔50帧取一帧出来计算
• &decomp：计算分解自由能
• print_res=‘1-110; 111’：蛋白-配体分别的序号，这个要从prmtop, inpcrd生成的pdb文件中看；单个则为单个氨基酸编号

(2) mmpbsa模块计算

$AMBERHOME/bin/MMPBSA.py -O -i mmpbsa.in -o FINAL_RESULTS_MMPBSA.dat -do FINAL_DECOMP_MMPBSA.dat -sp complex_solvated.prmtop -cp complex.prmtop -rp R.prmtop -lp L.prmtop -y md.crd