View source for 2016 AMBER tutorial with Beta Trypsin


For additional Rizzo Lab tutorials see [[AMBER Tutorials]].

In this tutorial, we will learn how to run a molecular dynamics simulation of a protein-ligand complex.  We will then post-process that simulation by calculating structural fluctuations (with RMSD) and free energies of binding (MM-GBSA).

==I. Introduction==

Yaping

'''AMBER'''


[http://ambermd.org/ Amber] -Assisted Model Building with Energy Refinement - is a multi-program suite for macromolecular simulations.Amber is distributed in two parts: AmberTools15 and Amber14. You can use AmberTools15 without Amber14, but not vice versa. Amber14 is the most recent version of the software and it includes new force fields such as ff14SB. In addition, in this release, more features from sander have been added to pmemd for both CPU and GPU platforms, including performance improvements, and support for extra points, multi-dimension replica exchange, a Monte Carlo barostat, ScaledMD, Jarzynski sampling, explicit solvent constant pH, GBSA, and hydrogen mass repartitioning. Support is also included for the latest Kepler, Titan and GTX7xx GPUs expanded options for Poisson-Boltzmann solvation calculations, accelerated molecular dynamics, additional features in sander pmemd code, and expanded methods for free energy calculations. Our lab is set up with Ambe r14 and the latest update of AmberTools15 which contains the programs such as antechamber and tleap to set up your simulation. 
The [http://ambermd.org/doc12/Amber14.pdf Amber 14 Manual]is available to get started with using Amber14. You can search the document for keywords such as "tleap" if you use Adobe Acrobat to view the file. Additionally, [http://ambermd.org/doc12/AmberTools15.pdf AmberTools Reference Manual]is another reference for the programs available under Amber tools. 

Here below are some of the programs available in both Amber and AmberTools: 
:  1.'''LEaP''': LEaP is an X-windows-based program that provides for basic model building and Amber coordinate parameter/topology input file creation. It includes a molecular editor which allows for building residues and manipulating molecules. 
:  2.'''ANTECHAMBER''': This program suite automates the process of developing force field descriptors for most organic molecules. It starts with structures (usually in PDB format), and generates files that can be read into LEaP for use in molecular modeling. The force field description that is generated is designed to be compatible with the usual Amber force fields for proteins and nucleic acids. 
:  3.'''SANDER''': Sander is short for Simulated annealing with NMR-derived energy restraints. This allows for NMR refinement based on NOE-derived distance restraints, torsion angle restraints, and penalty functions based on chemical shifts and NOESY volumes. Sander is also the "main" program used for molecular dynamics simulations, and is also used for replica-exchange, thermodynamic integration, and potential of mean force (PMF) calculations. Sander also includes QM/MM capability. 
:  4.'''PMEMD''': This is an extensively-modified version (originally by Bob Duke) of the sander program, optimized for periodic, PME simulations, and for GB simulations. It is faster than sander and scales better on parallel machines. 
:  5.'''PTRAJ'''and'''CPPTRAJ''': These are used to analyze MD trajectories, computing a variety of things, like RMS deviation from a reference structure, hydrogen bonding analysis, time-correlation functions, diffusional behavior, and so on. 
:  6.'''MM_PBSA''' and'''MM_PBSA.py''': These are scripts that automate post-processing of MD trajectories, to analyze energetics using continuum solvent ideas. It can be used to break energies energies into "pieces" arising from different residues, and to estimate free energy differences between conformational basins. 
:  7.'''NAB''': Originally named as "nucleic acid builder", NAB is a specialized language for writing programs that manipulate molecules and carry out molecular mechanics or distance-geometry based modeling. NAB provides and interface to Poisson-Boltzmann and RISM integral-equation solvent models. The "amberlite" package uses NAB to study protein-ligand interaction energetics. 

There is also a [http://lists.ambermd.org/mailman/listinfo/amber mailing list] available as an additional resource. What you can do with it is: you document your questions and sent to this mail address, some specialists of Amber will be assigned to reply your email and help you.

'''Beta Trypsin'''

Trypsin is a proteolytic enzyme, important for the digestion of proteins.In humans, the protein is produced in its inactive form, trypsinogen, within the pancrease. Trypsinogen enters the small intestine, via the common bile duct, where it converted to active trypsin.  Trypsin cleaves a terminal hexapeptide from trypsinogen to yield a single-chain [beta]-trypsin.  Subsequent autolysis produces other active forms having two or more peptide chains.  The two predominant forms of trypsin are [alpha]-trypsin, which has two peptide chains bound by disulfide bonds, and [beta]-trypsin which is a novel class of mechanism-based inhibitors of the serine proteases is developed using epitaxial selection. 

'''Organizing Directories'''

It makes things easier to organize your files in a clean and logical way. The following directory structure and naming scheme is a convenient way to organize your files. We could make these directories first before doing anything further
<pre>~username/AMS536-Spring2016/Amber_Tutorial/001.MOL.PREP/
                                           002.ANTE.TLEAP/ 
                                           003.PMEMD/       
                                           004.PTRAJ/
                                           005.MMGBSA/</pre>

==II. Structural Preparation==


====Antechamber, Parmchk, tLeap====

Omar, Katie

====Antechamber, Parmchk, tLeap====

Before beginning the Molecular Dynamics protocol using AMBER, you must first set up your files.
In your 001.mol.prep folder, add the following files from your docktutorial directory:

  1BJU.lig.mol2
  1BJU.rec.noH.pdb
 


IMPORTANT: delete any headers before the atoms/helix information.

Open 1BJU.rec.noH.pdb 



Before running tleap, make a new directory: mkdir 002.ANTE.TLEAP

  cp ~/../../tleap.lig.in
 
*tleap.lig.in
 source leaprc.ff14SB                                               #Load a force field
 source leaprc.gaff                  
 LIG = loadmol2 1BJU.lig.mol2                                       #Conformation file
 saveamberparm LIG 1BJU.lig.gas.leap.prm7 1BJU.lig.gas.leap.rst7    #Save the ligand gas phase AMBER topology and coordinate file
 solvatebox LIG TIP3PBOX 10.0                                       #Solvate the lig using TIP3P, solvent box radii 10 angstroms
 saveamberparm LIG 1BJU.lig.wat.leap.prm7 1BJU.lig.wat.leap.rst7    #Save the ligand water phase AMBER topology and coordinate files
 quit

  tleap –s –f tleap.lig.in > tleap.lig.out
The atom types which used in '1BJU.lig.mol2' file are not recognized. 

So the name of the atoms has to be redefined based on the gaff force field, '''antechamber''' can be used to solve the error.
  antechamber –i 1BJU.lig.mol2 –fi mol2 –o 1BJU.lig.gaff.mol2 –fo mol2
Change the content of tleap.lig.in to
(gaff file)

Re-run tleap:
  tleap –s –f tleap.lig.in > tleap.lig.out


So now run parmchk to fix the missing parameters
  parmchk –i 1BJU.lig.gaff.mol2 –f mol2 –o 1BJU.lig.ante.frcmod
Change the content of tleap.lig.in to
(load amber params)
run tleap.lig.in, and all the errors are fixed.
[[Image:error1&error2 fixed.png‎|thumb|center|600px]]
  cp ~/../../tleap.rec.in
*tleap.rec.in
 source leaprc.ff14SB                                                  #Load a force field
 REC = loadpdb ../001.chimera/1BJU.rec.noH.pdb                         #Conformation file for receptor 
 saveamberparm REC 1BJU.rec.gas.leap.prm7 4TKG.rec.gas.leap.rst7       #Save the receptor gas phase AMBER topology and coordinate file
 solvateBox REC TIP3PBOX 10.0                                          #Solvate the lig using TIP3P, solvent box radii 10 angstroms
 saveamberparm REC 1BJU.rec.wat.leap.prm7 4TKG.rec.wat.leap.rst7       #Save the receptor water phase AMBER topology and coordinate file
 quit
 
tleap –s –f tleap.rec.in > tleap.rec.out



run tleap.rec.in, and all the errors are fixed.
  
  cp ~/../../tleap.com.in
  tleap –s –f tleap.com.in > tleap.com.out

==III. Simulation using pmemd==

====PMEMD====
 01mi.in: equilibration
 &cntrl
 imin = 1, maxcyc = 1000, ntmin = 2,
 ntx = 1, ntc = 1, ntf = 1,
 ntb = 1, ntp = 0,
 ntwx = 1000, ntwe = 0, ntpr = 1000,
 cut = 8.0,
 ntr = 1,
 restraintmask = ':980-1181 & !@H=',
 restraint_wt = 5.0,



We first start by performing a series of minimization and short MD simulation steps for the receptor and ligand. The purpose of these steps is to remove steric clashes (minimization) and allow the receptor:ligand to equilibrate into a low energy state (MD simulation). After minimization, we can perform MD simulations in which we could record the trajectories to predict the folded state of protein.


1) In the first minimization step, we use a large restraint_wt on all heavy atoms (5.0 kcal/molA^2). Thus, only hydrogens are allowed to relax and only steric clashes involving hydrogens are removed.  

Create the input file 01mi.in

 01mi.in: equilibration 
 &cntrl
   imin = 1, maxcyc = 1000, ntmin = 2,
   ntx = 1, ntc = 1, ntf = 1,
   ntb = 1, ntp = 0,
   ntwx = 1000, ntwe = 0, ntpr = 1000,
   cut = 8.0,
   ntr = 1,
   restraintmask = ':1-224 & !@H=',
   restraint_wt = 5.0,

imin = 1:  flag to do minimization only

ntx = 1: read coordinates only from rst file

ntc = 1: flag for shake to perform bond length constraints. ntc =1 is the default, meaning SHAKE is not performed. In other words, minimization is performed on every bond. 

ntb = 1: flag for constant volume

ntp = 0: NO constant pressure

ntwx = 1000: saving snapshot to a trajectory file every 1000 steps.

ntwe = 0: human readable energy file is appended every 0 steps.

ntpr = 1000: newest calculated energy is stored in .info file every 1000 steps.

cut: cutoff for energy evaluation (Angstroms)

ntr = 1: flag for restraining specified atoms

restraintmask = ':1-224 & !@H=';  String that specifies restrained atoms (all atoms except hydrogens are restrained)

restraint_wt = 5.0;  restraint strength (kcal/molA^2)


2) The second step is a short MD simulation to bring the receptor:ligand system into a more equilibrated state. Create an input file 02md.in by copying and editing 01mi.in.
 cp 01mi.in 02md.in

 02md.in: equilibration
 &cntrl
   imin = 0, ntx = 1, irest = 0, nstlim = 50000,
   temp0 = 298.15, tempi = 298.15, ig = 71287,
   ntc = 2, ntf = 1, ntt = 1, dt = 0.001,
   ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5,
   ntwx = 1000, ntwe = 0, ntwr = 1000, ntpr = 1000,
   cut = 8.0, iwrap = 1,
   ntr = 1, nscm = 100,
   restraintmask = ':1-224 & !@H=', restraint_wt = 5.0,

imin = 0:  flag to do MD simulation

ntc = 2: flag for SHAKE to perform bond length constraints involving hydrogen

ntb = 2: flag for constant pressure

ntp = 1: constant pressure

ntr = 1: flag for restraining specified atoms

restraintmask = ':1-224 & !@H=';  String that specifies restrained atoms (all atoms except hydrogens are restrained)

restraint_wt = 5.0;  restraint strength (kcal/molA^2)


3) Create more minimization input files with decreasing restraint_wt (03mi.in, 04mi.in, 05mi.in). Create more MD simulation input files with different restraints (06md.in, 07md.in, 08md.in, 09md.in, 10md.in, 11md.in):

 03mi.in: equilibration
 &cntrl
   imin = 1, maxcyc = 1000, ntmin = 2,
   ntx = 1, ntc = 1, ntf = 1,
   ntb = 1, ntp = 0,
   ntwx = 1000, ntwe = 0, ntpr = 1000,
   cut = 8.0,
   ntr = 1,
   restraintmask = ':1-224 & !@H=',
   restraint_wt = 2.0,

 04mi.in: equilibration
 &cntrl
   imin = 1, maxcyc = 1000, ntmin = 2,
   ntx = 1, ntc = 1, ntf = 1,
   ntb = 1, ntp = 0,
   ntwx = 1000, ntwe = 0, ntpr = 1000,
   cut = 8.0,
   ntr = 1,
   restraintmask = ':1-224 & !@H=',
   restraint_wt = 0.1,

 05mi.in: equilibration
 &cntrl
   imin = 1, maxcyc = 1000, ntmin = 2,
   ntx = 1, ntc = 1, ntf = 1,
   ntb = 1, ntp = 0,
   ntwx = 1000, ntwe = 0, ntpr = 1000,
   cut = 8.0,
   ntr = 1,
   restraintmask = ':1-224 & !@H=',
   restraint_wt = 0.05,

 06md.in: equilibration
 &cntrl
   imin = 0, ntx = 1, irest = 0, nstlim = 50000,
   temp0 = 298.15, tempi = 298.15, ig = 71287,
   ntc = 2, ntf = 1, ntt = 1, dt = 0.001,
   ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5,
   ntwx = 1000, ntwe = 0, ntwr = 1000, ntpr = 1000,
   cut = 8.0, iwrap = 1,
   ntr = 1, nscm = 100,
   restraintmask = ':1-224 & !@H=', restraint_wt = 1.0,

 07md.in: equilibration
 &cntrl
   imin = 0, ntx = 5, irest = 1, nstlim = 50000,
   temp0 = 298.15, tempi = 298.15, ig = 71287,
   ntc = 2, ntf = 1, ntt = 1, dt = 0.001,
   ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5,
   ntwx = 1000, ntwe = 0, ntwr = 1000, ntpr = 1000,
   cut = 8.0, iwrap = 1,
   ntr = 1, nscm = 100,
   restraintmask = ':1-224 & !@H=', restraint_wt = 0.5,

 08md.in: equilibration
 &cntrl
   imin = 0, ntx = 5, irest = 1, nstlim = 50000,
   temp0 = 298.15, tempi = 298.15, ig = 71287,
   ntc = 2, ntf = 1, ntt = 1, dt = 0.001,
   ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5,
   ntwx = 1000, ntwe = 0, ntwr = 1000, ntpr = 1000,
   cut = 8.0, iwrap = 1,
   ntr = 1, nscm = 100,
   restraintmask = ':1-223 & @CA,C,N', restraint_wt = 0.1,

 09md.in: equilibration
 &cntrl
   imin = 0, ntx = 5, irest = 1, nstlim = 50000,
   temp0 = 298.15, tempi = 298.15, ig = 71287,
   ntc = 2, ntf = 1, ntt = 1, dt = 0.001,
   ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5,
   ntwx = 1000, ntwe = 0, ntwr = 1000, ntpr = 1000,
   cut = 8.0, iwrap = 1,
   ntr = 1, nscm = 100,
   restraintmask = ':1-223 & @CA,C,N', restraint_wt = 0.1,

 10md.in: production
 &cntrl
   imin = 0, ntx = 5, irest = 1, nstlim = 500000,
   temp0 = 298.15, tempi = 298.15, ig = 71287,
   ntc = 2, ntf = 1, ntt = 1, dt = 0.002,
   ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5,
   ntwx = 500, ntwe = 0, ntwr = 500, ntpr = 500,
   cut = 8.0, iwrap = 1,
   ntr = 1, nscm = 100,
   restraintmask = ':1-223 & @CA,C,N', restraint_wt = 0.1,

 11md.in: production
 &cntrl
   imin = 0, ntx = 5, irest = 1, nstlim = 500000,
   temp0 = 298.15, tempi = 298.15, ig = 71287,
   ntc = 2, ntf = 1, ntt = 1, dt = 0.002,
   ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5,
   ntwx = 500, ntwe = 0, ntwr = 500, ntpr = 500,
   cut = 8.0, iwrap = 1,
   ntr = 1, nscm = 100,
   restraintmask = ':1-223 & @CA,C,N', restraint_wt = 0.1,


4) Submitting the job to LIRED. Make a file called 1BJU.MD.qsub.csh:
 vim 1BJU.MD.qsub.csh
Writing information in it:
 #!/bin/tcsh
 #PBS -l walltime=04:00:00
 #PBS -l nodes=4:ppn=24 
 #PBS -q short
 #PBS -N 1BJU.MD
 #PBS -V
 #PBS -j oe
 cd /gpfs/home/guest40/AMS536-bjiang/Amber-tutorial/003.PMEMD
 /gpfs/software/intel_2016_update1/impi/5.1.2.150/bin64/mpirun -np 96 /gpfs/home/tmcgee/local/amber/bin/pmemd.MPI -O -i 01mi.in -o 01mi.out \
 -p ../002.ANTE.TLEAP/1BJU.com.wat.leap.prm7 -c ../002.ANTE.TLEAP/1BJU.com.wat.leap.rst7 \
 -ref ../002.ANTE.TLEAP/1BJU.com.wat.leap.rst7 -x 01mi.trj -inf 01mi.info -r 01mi.rst7
 /gpfs/software/intel_2016_update1/impi/5.1.2.150/bin64/mpirun -np 96 /gpfs/home/tmcgee/local/amber/bin/pmemd.MPI -O -i 02md.in -o 02md.out \
 -p ../002.ANTE.TLEAP/1BJU.com.wat.leap.prm7 -c 01mi.rst7 -ref 01mi.rst7 \
 -x 02md.trj -inf 02md.info -r 02md.rst7
 /gpfs/software/intel_2016_update1/impi/5.1.2.150/bin64/mpirun -np 96 /gpfs/home/tmcgee/local/amber/bin/pmemd.MPI -O -i 03mi.in -o 03mi.out \
 -p ../002.ANTE.TLEAP/1BJU.com.wat.leap.prm7 -c 02md.rst7 -ref 02md.rst7 \
 -x 03mi.trj -inf 03mi.info -r 03mi.rst7
 /gpfs/software/intel_2016_update1/impi/5.1.2.150/bin64/mpirun -np 96 /gpfs/home/tmcgee/local/amber/bin/pmemd.MPI -O -i 04mi.in -o 04mi.out \
 -p ../002.ANTE.TLEAP/1BJU.com.wat.leap.prm7 -c 03mi.rst7 -ref 03mi.rst7 \
 -x 04mi.trj -inf 04mi.info -r 04mi.rst7
 /gpfs/software/intel_2016_update1/impi/5.1.2.150/bin64/mpirun -np 96 /gpfs/home/tmcgee/local/amber/bin/pmemd.MPI -O -i 05mi.in -o 05mi.out \
 -p ../002.ANTE.TLEAP/1BJU.com.wat.leap.prm7 -c 04mi.rst7 -ref 04mi.rst7 \
 -x 05mi.trj -inf 05mi.info -r 05mi.rst7
 /gpfs/software/intel_2016_update1/impi/5.1.2.150/bin64/mpirun -np 96 /gpfs/home/tmcgee/local/amber/bin/pmemd.MPI -O -i 06md.in -o 06md.out \
 -p ../002.ANTE.TLEAP/1BJU.com.wat.leap.prm7 -c 05mi.rst7 -ref 05mi.rst7 \
 -x 06md.trj -inf 06md.info -r 06md.rst7
 /gpfs/software/intel_2016_update1/impi/5.1.2.150/bin64/mpirun -np 96 /gpfs/home/tmcgee/local/amber/bin/pmemd.MPI -O -i 07md.in -o 07md.out \
 -p ../002.ANTE.TLEAP/1BJU.com.wat.leap.prm7 -c 06md.rst7 -ref 06md.rst7 \
 -x 07md.trj -inf 07md.info -r 07md.rst7
 /gpfs/software/intel_2016_update1/impi/5.1.2.150/bin64/mpirun -np 96 /gpfs/home/tmcgee/local/amber/bin/pmemd.MPI -O -i 08md.in -o 08md.out \
 -p ../002.ANTE.TLEAP/1BJU.com.wat.leap.prm7 -c 07md.rst7 -ref 07md.rst7 \
 -x 07md.trj -inf 08md.info -r 08md.rst7
 /gpfs/software/intel_2016_update1/impi/5.1.2.150/bin64/mpirun -np 96 /gpfs/home/tmcgee/local/amber/bin/pmemd.MPI -O -i 09md.in -o 09md.out \
 -p ../002.ANTE.TLEAP/1BJU.com.wat.leap.prm7 -c 08md.rst7 -ref 08md.rst7 \
 -x 07md.trj -inf 09md.info -r 09md.rst7
 /gpfs/software/intel_2016_update1/impi/5.1.2.150/bin64/mpirun -np 96 /gpfs/home/tmcgee/local/amber/bin/pmemd.MPI -O -i 10md.in -o 10md.out \
 -p ../002.ANTE.TLEAP/1BJU.com.wat.leap.prm7 -c 09md.rst7 -ref 09md.rst7 \
 -x 10md.trj -inf 10md.info -r 10md.rst7
 /gpfs/software/intel_2016_update1/impi/5.1.2.150/bin64/mpirun -np 96 /gpfs/home/tmcgee/local/amber/bin/pmemd.MPI -O -i 11md.in -o 11md.out \
 -p ../002.ANTE.TLEAP/1BJU.com.wat.leap.prm7 -c 10md.rst7 -ref 10md.rst7 \
 -x 11md.trj -inf 11md.info -r 11md.rst7
Changing "cd /gpfs/home/guest40/AMS536-bjiang/Amber-tutorial/003.PMEMD" to your own path.

Submit your job by command:
 qsub 1BJU.MD.qsub.csh

==IV. Simulation Analysis==

===Ptraj===

Ptraj is an abbreviation for Process Trajectory. The points or frames in the trajectory are referred to as snapshots in an MD simulation. First create a new directory called 004.PTRAJ\

 mkdir 004.PTRAJ
 cd 004.PTRAJ

Once inside the the directory, we are going to combine the two 1 nanosecond trajectories together. We will also strip away all the waters, and create a .strip output file. 

ptraj.strip.wat.in:
 trajin ../003.PMEMD/10md.trj 1 1000 1
 trajin ../003.PMEMD/11md.trj 1 1000 1
 trajout 1BJU.trj.strip nobox
 strip :WAT

The "1 1000 1" in the input files dictate which frames will be used for Ptraj, where the first two numbers "1 1000" are the beginning and ending snapshots. The last number "1" is the frequency in which your snapshots will be saved. 

ptraj.rec.rmsd.in:
 trajin 1BJU.trj.strip
 trajout 1BJU.com.trj.stripfit 
 reference ../002.ANTE.TLEAP/1BJU.com.gas.leap.rst7
 rms reference out 1BJU.rmsd.CA.dat :1-223@CA

ptraj.lig.rmsd.in:
 trajin 1BJU.com.trj.stripfit
 reference ../002.ANTE.TLEAP/1BJU.com.gas.leap.rst7
 rms reference out 1BJU.lig.rmsd.dat :224@C*,N*,O*,S* nofit

ptraj.strip.rec.in:
 trajin 1BJU.com.trj.stripfit 1 2000 1
 trajout 1BJU.lig.trj.stripfit
 strip :1-223

ptraj.strip.lig.in:
 trajin 1BJU.com.trj.stripfit 1 2000 1
 trajout 1BJU.rec.trj.stripfit
 strip :224

ptraj.dist.in:
 trajin 1BJU.com.trj.stripfit
 distance 171OD2_224N1 :171@OD2 :224@N1 out 171OD2_224N1.dat
 distance 171OD1_224N2 :171@OD1 :224@N2 out 171OD1_224N2.dat

To run the file
  cpptraj ../*.in > *.log

===MM-GBSA Energy Calculation===
Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) is the standard method to calculate binding affinities between ligands and receptors. The more negative the value the  stronger the binding affinity. Topology files are needed for the ligand, receptor, and ligand-receptor complex. 

Begin by creating a new directory

   mkdir 005.MMGBSA/

as well as an input file

   vim.gbrescore.in

the following parameters will provide a good estimate of the energy of binding based on solvation enthalpy and binding enthalpy

   Single point GB energy calc
   &cntrl
 ntf    = 1,        ntb    = 0,        ntc    = 2,
 idecomp= 0,
 igb    = 5,        saltcon= 0.00,
 gbsa   = 2,        surften= 1.0,
 offset = 0.09,     extdiel= 78.5,
 cut    = 99999.0,  nsnb   = 99999,
 imin   = 5,        maxcyc = 1,        ncyc   = 0,
 /

We can now submit a job to the cluster we are working on. This script, inside a file named run.sander.rescore.csh , is designed for LIRED, however it can be modified for whatever supercomputer you are working on

   #! /bin/tcsh
 #PBS -l nodes=1:ppn=1
 #PBS -l walltime=48:00:00
 #PBS -o zzz.qsub.out
 #PBS -e zzz.qsub.err
 #PBS -V
 #PBS -N mmgbsa
 
 set workdir = /nfs/guest46/amber_tutorial/005.mmgbsa
 
 cd $workdir
 
 sander -O -i gb.rescore.in \
 -o gb.rescore.out.com \
 -p ../002.tleap/4TKG.com.gas.leap.prm7 \
 -c ../002.tleap/4TKG.com.gas.leap.rst7 \
 -y ../004.ptraj/4TKG.com.trj.stripfit \
 -r restrt.com \
 -ref ../002.tleap/4TKG.com.gas.leap.rst7 \
 -x mdcrd.com \
 -inf mdinfo.com

 sander -O -i gb.rescore.in \
 -o gb.rescore.out.lig \
 -p ../002.tleap/4TKG.lig.gas.leap.prm7 \
 -c ../002.tleap/4TKG.lig.gas.leap.rst7 \
 -y ../004.ptraj/4TKG.lig.trj.stripfit \
 -r restrt.lig \
 -ref ../002.tleap/4TKG.lig.gas.leap.rst7 \
 -x mdcrd.lig \
 -inf mdinfo.lig 
 
 sander -O -i gb.rescore.in \
 -o gb.rescore.out.test.rec \
 -p ../002.tleap/4TKG.rec.gas.leap.prm7 \
 -c ../002.tleap/4TKG.rec.gas.leap.rst7 \
 -y ../004.ptraj/4TKG.rec.trj.stripfit \
 -r restrt.rec \
 -ref ../002.tleap/4TKG.rec.gas.leap.rst7 \
 -x mdcrd.rec \
 -inf mdinfo.rec 
 
 exit

Then submit the job as a .csh 

   qsub run.sander.rescore.csh

Your ligand, receptor, and receptor-ligand complex will all have .out files. The format should look 

   FINAL RESULTS
  
  
   
   NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER
      1       5.9132E+03     2.0005E+01     1.2640E+02     C         159
  
 BOND    =      661.8980  ANGLE   =     1751.7992  DIHED      =     2581.7692
 VDWAALS =    -1696.6585  EEL     =   -13958.9335  EGB        =    -3125.9524
 1-4 VDW =      747.0185  1-4 EEL =     7750.8118  RESTRAINT  =        0.0000
 ESURF   =    11201.4791
   minimization completed, ENE= 0.59132314E+04 RMS= 0.200047E+02

We can now use a system of equations to elucidate various energies related to our system. 
VDWAALS = ΔGvdw 
E(ELS) = ΔGcoulombic
E(GB) = ΔGpolar
SASA = ESURF


ΔGnonpolar can be estimated by:

ΔGnonpolar = (Solvent Accessible Surface Area)*0.00542 + 0.92

This equation is based on empirical observations relating a molecules surface area and solvation energy

ΔGnonpolar and ΔGmmgbsa are related by anther equation:

ΔGmmgbsa = ΔGvdw + ΔGcoul + ΔGpolar + ΔGnonpolar

all the information needed is in our output files.

ΔΔGbind can now be solved via the equation:

ΔΔGbind = ΔG(mmgbsa complex) – (ΔG(mmgbsa ligand) + ΔG(mmgbsa receptor)	(3)

ΔΔGbind is a function of ligand position, so it is useful to plot it. For publications, it is standard protocol to include the mean and standard deviation for your ΔΔGbind.

Lets make a script named get.mmgbsa.csh to extract this sophisticated energetic estimation from the three initial output files of the calculations obtained solve for ΔΔGbind:

    bash get.mmgbsa.csh

The output file will be .dat and we can easily use the program xmgrace to plot the information visually
    
   xmgrace MMGBSA_vs_time.dat

==V. Frequently Encountered Problems==