分享自用的amber脚本

lrf1980

最近在学习amber，分享一下自己用的amber脚本，主要用MMPBSA做free binding energy的计算。两个脚本，一个是准备文件用的脚本，另外一个是做模拟的脚本
文件准备脚本，可以自己copy下面的代码，然后命名成＊.sh， 记得加上可执行权限，chmod +x *.sh, 当然也可以自己一行行的在终端输入
这里用到的计算配体电荷是用的amber自带的方法，bcc 电荷。 ${charge} 这部分是表示配体的电荷。假设配体带一个负电荷，那么这里的输入就是：
antechamber -i ligand.mol2 -fi mol2 -o lig.mol2 -c bcc -nc -1 -rn LIG -df 2 -pf y -s 0

#! /bin/bash
# A script for MD simulation by Amber
#
# use ligand.mol2 for ligand molecule and receptor as target protein
# step 1. prepare the protein-ligand complex either from docking results or PDB database
#
# step 2. use MOE／Yasara or any other related software to do the protanate3D to get ligand.mol2 and receptor.pdb file
#
# step 3. save ligand in mol2 and receptor in pdb format as mentioned above
#
# step 4. calculate the bcc charge
#
antechamber -i ligand.mol2 -fi mol2 -o lig.mol2 -fo mol2 -c bcc -nc ${charge} -rn LIG -df 2 -pf y -s 0
#
# step 5. get topology parameter for ligand (indicate the ligand name by using LIG)
#
parmchk -i lig.mol2 -f mol2 -o lig.frcmod
#
# step 6. tleap for other files and parameters generation
cat >leap.in <<EOF
source leaprc.ff99SB
source leaprc.gaff
lig = loadmol2 lig.mol2
loadamberparams lig.frcmod
check lig
rec = loadpdb receptor.pdb
com=combine {rec lig}
check com
saveamberparm lig lig.prmtop lig.inpcrd
saveamberparm rec rec.prmtop rec.inpcrd
saveamberparm com com.prmtop com.inpcrd
savepdb com com.pdb
solvateoct com TIP3PBOX 10.0
addions com Na+ 0
addions com Cl- 0
saveamberparm com com_wat.prmtop com_wat.inpcrd
savepdb com com_wat.pdb
quit
END
END
EOF
tleap -f leap.in

模拟用的脚本
在使用该脚本之前，首先需要知道加了水盒子的配体-蛋白质文件里的一些基本情况。
脚本中标注下划线的两个地方在每次使用前，都需要根据实际情况修改
1. RES 326 13200， 这部分是复合物中起始水和最后一个水的编号。这一步是在做蛋白质的minimization，固定水
2. RES 1 308， 这部分是复合物中的所有氨基酸及配体的编号。这一步是在做溶剂及离子的minimization，固定蛋白质
根据自己体系改好这两个部分的参数，就可以直接运行脚本文件了。

#! /bin/bash
####################### Minimize Process ##########################
# step 1. initial protein minimization
cat > min0.in <<EOF
Complex: initial protein minimization
&cntrl
imin = 1,
ntpr = 25,
ntb =1, cut = 12,
ntmin = 2, maxcyc = 6000, ncyc = 600,
/
Water molecules fixed
500.0
RES 326 13200
END
END
EOF

# step 2. initial solvent + ions minimization
cat > min1.in <<EOF
Complex: initial slovent + ions minimization
&cntrl
imin =1,
ntc = 2, ntf = 2,
ntpr = 250, ntwx = 1000,
ntb = 1, cut = 12,
ntr = 1,
maxcyc = 10000, ncyc = 1000,
/
Protein fixed
500.0
RES 1 308
END
END
EOF

# step 3. whole system minimization
cat > min2.in <<EOF
Complex: whole system minimization
&cntrl
imin = 1,
ntc = 2, ntf = 2,
ntpr = 250, ntwx = 1000,
ntb =1, cut = 12,
maxcyc = 10000, ncyc = 1000
/
END
END
EOF

###################### Heat Process ##############################
# step 4. heat from 0K to 100K in 50 ps
cat > heat.in<<EOF
50 ps NVT heat from 0K to 300K
&cntrl
  imin = 0, irest = 0,
  ntx = 1, ntb = 1,
  iwrap = 1,
  cut = 12,
  ntc = 2, ntf = 2,
  tempi = 0.0, temp0 = 300.0,
  ntt = 3,
  gamma_ln = 1.0,
  nstlim = 25000, dt = 0.002,
  ntpr = 500, ntwx = 500, ntwr = 500,
  nmropt = 1,
/
&wt TYPE='TEMP0', istep1 = 0, istep2 = 25000,
value1 = 0.1, value2 = 300.0, /
&wt TYPE='END' /
END
END
EOF
#
#
########################## Density ############################
# step 5. start to equillibrate the system
cat > dens.in<<EOF
100 ps NPT
&cntrl
  imin = 0, irest = 1, ntx = 5,
  ntb = 2, pres0 = 1.0, ntp = 1,
  taup = 2.0,
  iwrap = 1,
  cut = 12.0,
  ntc = 2, ntf = 2,
  tempi = 300.0, temp0 = 300.0,
  ntt = 3, gamma_ln = 1.0,
  nstlim = 50000, dt = 0.002,
  ntpr = 500, ntwx = 500, ntwr = 500,
/
END
END
EOF
#
######################## MD run ##############################
# step 6. MD run for the system
cat > equil.in<<EOF
1ns NVT MD
&cntrl
  imin = 0,
  irest = 1, ntx = 5,
  ntb = 1, ntp = 0,
  cut = 12.0,
  ntr = 0,
  iwrap = 1,
  ntc = 2, ntf = 2,
  tempi = 300.0, temp0 = 300.0,
  ntt= 3, gamma_ln = 1.0,
  nstlim = 500000, dt = 0.002,
  ntpr = 1000, ntwx = 1000, ntwr = 10000,
/
END
END
EOF
#
pmemd.cuda -O -i min0.in -p com_wat.prmtop -c com_wat.inpcrd -ref com_wat.inpcrd -r min0.rst -o min0.out
pmemd.cuda -O -i min1.in -p com_wat.prmtop -c min0.rst -ref min0.rst -r min1.rst -o min1.out
pmemd.cuda -O -i min2.in -p com_wat.prmtop -c min1.rst -ref min1.rst -r min2.rst -o min2.out
#
ambpdb -p com_wat.prmtop <min0.rst>min0.pdb
ambpdb -p com_wat.prmtop <min1.rst>min1.pdb
ambpdb -p com_wat.prmtop <min2.rst>min2.pdb
#
pmemd.cuda -O -i heat.in -p com_wat.prmtop -c min2.rst -ref min2.rst -r heat.rst -o heat.out -x heat.mdcrd
#
pmemd.cuda -O -i dens.in -p com_wat.prmtop -c heat.rst -ref heat.rst -r dens.rst -o dens.out -x dens.mdcrd
#
pmemd.cuda -O -i equil.in -p com_wat.prmtop -c dens.rst -ref dens.rst -r equil1.rst -o equil1.out -x equil1.mdcrd
pmemd.cuda -O -i equil.in -p com_wat.prmtop -c equil1.rst -ref equil1.rst -r equil2.rst -o equil2.out -x equil2.mdcrd
pmemd.cuda -O -i equil.in -p com_wat.prmtop -c equil2.rst -ref equil2.rst -r equil3.rst -o equil3.out -x equil3.mdcrd
pmemd.cuda -O -i equil.in -p com_wat.prmtop -c equil3.rst -ref equil3.rst -r equil4.rst -o equil4.out -x equil4.mdcrd
pmemd.cuda -O -i equil.in -p com_wat.prmtop -c equil4.rst -ref equil4.rst -r equil5.rst -o equil5.out -x equil5.mdcrd

lrf1980

Ares 发表于 2015-10-9 17:08
有没有cpptraj H键分析脚本?

有，等会儿换台电脑贴一下,下面脚本包括基本的轨迹分析及氢键分析。请根据具体情况修改相应的轨迹及残基编号
#use ptraj com_wat.prmtop ptraj.in to analyze the traj and generate the RMSD and RMSF
# Read in the trajectory file
#trajin ../heat.mdcrd 1 2500 10
#trajin ../dens.mdcrd 1 2500 10
trajin ../equil1.mdcrd 1 500 1
trajin ../equil2.mdcrd 1 500 1
trajin ../equil3.mdcrd 1 500 1  

center origin :1-126
image center familiar


##rmsd and rmsf analysis for all system

reference min1.rst [min1]
rms ref [min1] :1-125@CA,C,O,N,H out rmsd-backbone.dat
rms ref [min1] :1-125&!@CA,C,O,N,H out rmsd-sidechain.dat
rms ref [min1] :126<@5.0 out rmsd-pocket.dat
rms ref [min1] :126 out rmsd-lig.dat

atomicfluct out rmsf-backbone.dat :1-125@CA,C,O,N,H byres
atomicfluct out rmsf-sidechain.dat :1-125&!@CA,C,O,N,H byres
atomicfluct out rmsf-pocket.dat :126<@5.0 byres


##hbond analysis

#analysis ligand hbond status
hbond acceptormask :126 out hbond_acceptor.dat avgout avg_hbond_acceptor.dat
hbond donormask :126 out hbond_donor.dat avgout avg_hbond_donor.dat
hbond out nhb.dat avgout avg_solute.dat solventacceptor :WAT@O solventdonor :WAT solvout avg_solvent.dat bridgeout bridge.dat

#strip :WAT
#strip :Cl-

## save multiple pdb files from mdcrd
#trajout pdb_equil3/equil3 pdb multi

lrf1980

laoman 发表于 2015-7-13 08:30
简单的理解就是二进制文件是偏向于机器的，对机器来说，二进制文件比文本更易读写（可以加快MD速度），体 ...

ioutfm
         The format of coordinate and velocity trajectory files (mdcrd, mdvel and inptraj).
          As of Amber 9, the binary format used in previous versions is no longer supported;
          binary output is now in NetCDF trajectory format. While not the default option,
          binary trajectory files have many advantages: they are smaller, higher precision,
          much faster to read and write, and able to accept a wider range of coordinate (or
          velocity) values than formatted trajectory files.
          = 0 (default) Formatted ASCII trajectory
          = 1 Binary NetCDF trajectory

不知道这个是不是能满足你说的要求

laoman

lrf1980 发表于 2015-7-12 20:41
详细分享一下好在哪里？
：）

简单的理解就是二进制文件是偏向于机器的，对机器来说，二进制文件比文本更易读写（可以加快MD速度），体积也没有文本的那么大。Gromacs的的轨迹就是二进制存取的。但我想知道的是Amber的轨迹要怎么写in文件或命令才能存为二进制？manual里面有说到吗？

greatzdl

轨迹存为二进制格式会更好

lrf1980

greatzdl 发表于 2015-7-12 17:36
轨迹存为二进制格式会更好

详细分享一下好在哪里？
：）

greatzdl

lrf1980 发表于 2015-7-13 09:44
ioutfm
         The format of coordinate and velocity trajectory files (mdcrd, mdvel and inptraj).
...

就是这个

lrf1980

最近有朋友提到，在做minimization的时候应该先做溶剂，然后蛋白，我觉得也有道理。大家可以改变顺序，在自己的体系里试试看，是否结果有差别

Ares

有没有cpptraj H键分析脚本?

twq0228

amber动力学跑的脚本和分析的最新的脚本都发我一下！我邮箱twq0228@163.com 谢谢！