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