Biomolecular Simulation with Amber
(Advanced Level)
If you have any questions/suggestions feel free to contact at bhakatsoumendranath@gmail.com
Please went through the Amber manual available at http://www.ambermd.org for details about keywords and manipulations needed during simulation. Also you can subscribe to Amber mailing list.
Enzyme and ligand preparation
Open the receptor in Chimera.
If there is receptor+ligand then in first step:
1. Delete the ligand and remove all hydrogens from receptor —-save as “rec.pdb”
Open the “rec.pdb” in H++ server in order to generate correct protonation states for amino acids at a particular pH.
2. Then open the complex and delete the receptor in order to get the ligand. Then click on structure editing>add hydrogen>add charges (AM1-BCC)>save as “lig.mol2”
System preparation
Make sure you load all modules in order to run Amber.
Example: location of modules
~/amber12/bin/antechamber
~/amber12/bin/parmchk
~/amber12/bin/tleap
~/amber12/bin/cpptraj
Steps:
1. antechamber -i lig.mol2 -fi mol2 -o LIG.mol2 -fo mol2 -j 4 -at gaff -c AM1-BCC -nc 0
2. parmchk -i LIG.mol2 -f mol2 -o LIG.frcmod
3. tleap -s –f tleap.all
#Input for tleap.all
source leaprc.ff99SB
loadamberparams frcmod.ionsjc_tip3p
source leaprc.gaff
LIG = loadmol2 LIG.mol2
loadamberparams LIG.frcmod
check LIG
receptor = loadPDB rec.pdb
complex = combine {receptor LIG}
set default PBRadii mbondi2
saveAmberParm LIG LIG.top LIG.crd
receptor = loadPDB rec.pdb
saveAmberParm receptor rec.top rec.crd
saveAmberParm complex com.top com.crd
savepdb complex complex_gas.pdb
charge complex
addIons2 complex Na+ 0
solvateBox complex TIP3PBOX 8.0
saveAmberParm complex com_solvated.top com_solvated.crd
savepdb complex com_solvated.pdb
quit
New tleap Input with improved sidechain potential with TIP4P water model
#Tleap with FF14SB, improved sidechains, TIP4PEW
source leaprc.ff14SB
loadAmberParams frcmod.tip4pew
loadAmberParams frcmod.ionsjc_tip4pew
complex = loadPDB rec_final.pdb
# This part is for Cys-Cys sulphur bridges.
savepdb x pdbout
bond x.275.SG x.325.SG
bond x.365.SG x.161.SG
bond x.388.SG x.223.SG
check x
saveAmberParm receptor rec.top rec.crd
charge complex
addIons2 complex NA 0
solvateOct complex TIP4PBOX 10.0
saveAmberParm complex com_solvated.top com_solvated.crd
savepdb complex com_solvated.pdb
quit
Input of Partial minimization
# Res 1 198 describes residue number in case protein+ligand complex. E.g. in case of inhibitor bound HIV-protease the total residue number for receptor is 198 and the ligand is 1. Thus total residue count 198+1=199
Initial minimization of MMP3 (MMMM): solvent molecules and added ions
&cntrl
imin = 1,
maxcyc = 2500,
ncyc = 750,
ntb = 1,
ntr = 1,
cut = 12.0,
/
Hold the Protein fixed
10.0
RES 1 199
END
END
Input for Full minimization
Full minimization of MMP3 (MMMM): protein, ligand, solvent molecules and added ions
&cntrl
imin = 1,
maxcyc = 200,
ncyc = 50,
ntb = 1,
ntr = 0,
cut = 12.0,
Drms = 0.0001,
/
END
Input for Heating stage
#tempi(temperature input)=0 (0K); tempo(temperature output)=300 (300K)
#Res 1 199 is similar to the explanation provided above
Heating Step of MMP3 (MMMM): stage-5
&cntrl
imin= 0,
irest=1,
NTX=1,
ntb= 1,
NTPR=500,
NTWX=500,
NTWR=500,
ntr=1,
Tempi=0.0,
Temp0=300.0,
NTT=3,
gamma_ln=1.0,
NTC=2,
NTF=2,
cut= 12.0,
nstlim=2500,
dt=0.002,
/
Keep Protein and inhibitor fixed with weak restraints
10.0
RES 1 557
END
END
Input for Equilibration stage
Equilibration Step of MMP3 (MMMM): stage-1
&cntrl
imin= 0,
irest=1,
NTX=7,
ntb=2,
ntp=1,
PRES0=1.0,
TAUP=2.0,
NTPR=500,
NTWX=500,
ntr=0,
Tempi=300.0,
Temp0=300.0,
NTT=3,
gamma_ln=1.0,
NTC=2,
NTF=2,
cut=12.0,
nstlim=250000,
dt=0.002
/
Input for MD production stage
#change nstlim as per simulation time (nstlim is number of steps). 5ns=2500000 (nstlim)
Equilibration Step of MMP3 (MMMM): stage-1
&cntrl
imin= 0,
irest=0,
ig=-1,
NTX=7, ntxo=2, ioutfm=1,
ntb=2,
ntp=1,
PRES0=1.0,
TAUP=2.0,
NTPR=500,
NTWX=500,
ntr=0,
Tempi=300.0,
Temp0=300.0,
NTT=3,
gamma_ln=1.0,
NTC=2,
NTF=2,
cut=12.0,
nstlim=2500000,
dt=0.002
/
To transfer .rst file to .pdb file
~/amber14/bin/ambpdb -p com_solvated.top <Full_Mini.rst> Full_Mini.pdb
Example of MD submission script
#!/bin/bash
#PBS -N <Name>
#PBS -q gpu
#PBS -l nodes=1:ppn=1:gpus=1:exclusive_process
#PBS -l walltime=72:00:00
#PBS -V
#PBS -k oe
echo $CUDA_VISIBLE_DEVICES
nvidia-smi
# see the path of the folder
cd /cbio/jclab/projects/collab/bhakat/apo_plasmepsin
#see the path of Amber
export AMBERHOME=/cbio/jclab/projects/collab/amber14
#Run the minimization with just CPU AMBER.
#$AMBERHOME/bin/pmemd -O -i Partial_Mini.in -o Partial_Mini.out -p rec_solvated.top -c rec_solvated.crd -ref rec_solvated.crd -r Partial_Mini.rst
#$AMBERHOME/bin/pmemd -O -i Full_Mini.in -o Full_Mini.out -p rec_solvated.top -c Partial_Mini.rst -r Full_Mini.rst
#Run MD using GPU AMBER.
#$AMBERHOME/bin/pmemd.cuda -O -i Heating.in -o Heating.out -p rec_solvated.top -c Full_Mini.rst -r Heating.rst -ref Full_Mini.rst -x Heating.mdcrd
#$AMBERHOME/bin/pmemd.cuda -O -i equil.in -o equil.out -p rec_solvated.top -c Heating.rst -r equil.rst -x equil.mdcrd
$AMBERHOME/bin/pmemd.cuda -O -i md.in -o md50.out -p rec_solvated.top -c equil.rst -r md50.rst -x md50.nc
Submit with qsub command
Check the queue with qstat/qme command
To delete job qdel <job number>
Analysis Process
strip.ptrj
#specifying the input .mdcrd file
trajin md50.mdcrd
#stripping water and Cl- and creating a com_solvated.top with a prefix stripped
strip :WAT,Cl- out prefix stripped
# to generate a stripped image of com_solvated.top
autoimage
#dry .mdcrd files
trajout md50_dry.mdcrd
Execute like ~/amber12/bin/cpptraj com_solvated.top strip.ptrj
Calculation of RMSD, RMSF, average PDB
first strip water from Full_Mini.rst file by doing the following (this step will generate dry version of Full_Mini.pdb which will be compatible with stripped topology and co-ordinate)
trajin Full_Mini.rst
strip :WAT,Cl-
trajout Full_Mini_dry.rst
then generate the average PDB by following script (this will generate av.pdb)
trajin md1_dry.mdcrd
reference Full_Mini_dry.rst
rms reference out rms_min.dat @CA
average av.pdb * pdb nobox
then calculate both RMSF and B-factor (this step will calculate RMSF and bfactor aligning on av.pdb)
trajin md1_dry.mdcrd
reference av.pdb
rms reference out rms_min.dat @CA
atomicfluct out bfactor.dat :1-198@CA byres bfactor
atomicfluct out RMSF.dat :1-198@CA byres
Note: rms-min.dat is the RMSD after aligning not the RMSF.
Calculation of inter-atomic distances, angles, dihedral angles, radius of gyration, solvent accessible surface area, secondary structure evolution, generating snapshot PDBs for a particular point of time.
Vi analysis.ptrj with the following input
trajin md50_dry.mdcrd
#calculation of distance
distance DIST @1 @784 out ptraj_distance.out
distance DIST @1603 @784 out ptraj_distance1.out
# if you want radius of gyration
radgyr out RoG.dat mass nomax
radgyr out @1918 @1800 RoG1.dat mass nomax
#generating snapshot PDBs at a particular point of time
outtraj snapshot1.pdb onlyframes 12345
outtraj snapshot2.pdb onlyframes 12500
#calculation of solvent accessible surface area
surf :1-198 out surf1.dat
#secondary structure evolution
secstruct :1-560 out dssp.gnu sumout dssp.agr
#dihedral angle
dihedral phi @1205 @535 @3334 @4451 out phi1.dat
dihedral psi @1205 @535 @3334 @4451 out psi1.dat
#calculating angle between atoms
angle @1205 @535 @334 out angle.dat
#calculating hbond
hbond out hbond.dat
Execute like ~/amber12/bin/cpptraj stripped.com_solvated.top analysis.ptrj
mmpbsa.in (#MMPBSA input) (please check Amber tutorial to manipulate the input)
Input file for running PB and GB in serial
&general
startframe=1, endframe=50000, interval=50, keep_file=0, entropy=1,
/
&gb
igb=2,saltcon=0.150,
/
&pb
istrng=0.100,
/
&decomp
idecomp=1, print_res="323;113;270;99;109;110;95;112;94;91;98;260;256;255;259;252;273;274;275",
dec_verbose=1,csv_format=0,
/
vi run_mmpbsa.sub
#!/bin/bash
#PBS -N <Name>
#PBS -q gpu
#PBS -l nodes=1:ppn=1:gpus=1:exclusive_process
#PBS -l walltime=72:00:00
#PBS -V
#PBS -k oe
echo $CUDA_VISIBLE_DEVICES
nvidia-smi
#see the path of folder
cd /cbio/jclab/projects/collab/bhakat/hiv_wild
#see the path of Amber
export AMBERHOME=/cbio/jclab/projects/collab/amber14
#Run MD using GPU AMBER.
$AMBERHOME/bin/MMPBSA.py.MPI -O -i mmpbsa.in -o mmpbsa_final.dat -sp stripped.com_solvated.top -cp com.top -rp rec.top -lp LIG.top -y md50_dry.mdcrd > mmpbsa2.log
Submit it like qsub run_mmpbsa.sub
Combining trajectories
First remove water and ions from trajectory (all in interactive mode)
vi combine.ptrj
trajin md50_dry.mdcrd
trajin md51_dry.mdcrd
trajin md52_dry.mdcrd
trajout md_total.mdcrd
Execute like ~/amber14/bin/cpptraj stripped.com_solvated.top combine.ptrj
