Difference between revisions of "2023 AMBER tutorial 3 with PDBID 2P16"

Latest revision as of 15:04, 28 April 2023

Introduction

Howdy! If you've gotten this far with your MD project, you've now moved out of the realm of DOCK6 and into the world of molecular dynamics (MD). This tutorial will cover how to use the AMBER 16 software package to simulate your protein-ligand complex on a sub-microsecond timeframe and use the resulting data to calculate binding energies.

Preparing the Structure

Before you start any preparation, you need to keep yourself organized to make sure that you're running the right files and jobs for your MD simulation, otherwise, it may explode.

To do this, you can either copy the example AMBER directory to your scratch directory (this may change per year):

/gpfs/projects/AMS536/2023

Or, you can create the following directories:

mkdir 000.parameters
mkdir 001.tleap_build
mkdir 002.equilibration
mkdir 003.production
mkdir 004.analysis
mkdir zzz.master

For this tutorial, we will be copying the example AMBER setup into our folder to run our complex.

Simulation Parameters

If you're looking at this tutorial after the DOCK6 tutorial and are looking to use the original ligand, then you may want to fetch a fresh PDB file for this. Otherwise, make sure you Cartesian minimize any new ligands before sending them through the gauntlet that is MD simulation and preparation.

All of your original (stock/DOCK6 generated) files should be put into the zzz.master folder, but make sure you're in the proper directory for each step!

Receptor File Generation

To generate a fresh copy of the receptor file, go to the Protein Data Bank (PDB) website and download the PDB structure for 2P16, 2p16.pdb, which will contain the ligand-receptor complex. Open the pdb file in Chimera and go to

 Select -> Structure -> Main

Alternatively, you can instead for this protein:

 Select -> Residues -> All Standard Residues

which will select the receptor protein, and then:

 Select -> Invert (all models)

to select everything except the receptor. These will all be deleted:

 Actions -> Atoms/Bonds -> delete

Now, save your receptor as a PDB. You can save it as a mol2, but you are much more likely to bump into problems down the road:

 File -> Save PDB... 

and save as 2P16_rec.pdb. Close Chimera. You may have thought you were done, but there are still issues with your protein. Specifically, AMBER will not read disulfide bonds between cysteines properly. To fix this, you will need to manually edit your PDB file (using Notepad, vi, nano, etc.) and edit all disulfide-bonded cysteine by changing the residue from CYS to CYX. In addition, there will be some lines to add to your tleap.build.in file later on to make sure that AMBER properly recognizes the disulfide bonds. Save these changes.

Ligand File Generation

To generate the ligand file, open the 2P16.pdb file again, and go to

 Select -> Structure -> Main

Alternatively, you can instead for this protein:

 Select -> Residues -> All Standard Residues

Then you will delete the protein:

 Actions -> Atoms/Bonds -> delete

You should also delete all non-essential cofactors, ions, and waters:

 Select -> Residue -> HOH -> Actions -> Atoms/Bonds -> delete

 Select -> Residue -> CA -> Actions -> Atoms/Bonds -> delete

to delete everything except for the GG2 residue, which is the ligand. After this is completed, you will be adding hydrogens and partial charges:

 Tools > Structure editing > DockPrep

Make sure to select the same parameters as in the DOCK6 prep (ff14SB and AM1BCC). Save this file:

 File -> Save Mol2... 

and save as 2P16_lig_withH.mol2. Close Chimera. Transfer them back into your Seawulf directory (./zzz.master). The next step involved getting your ligands ready for tleap.

Force Field Parameterization

Before running TLeap, we need to generate appropriate force field parameters for the ligand. You can achieve this using antechamber- a program built specifically for parametrization of small molecules, followed by parmchk to verify if it ran correctly. Make sure that the -nc flag corresponds to the net charge of the molecule.

If your ligand happens to be a peptide with standard residues, you can skip this step.

Switch to the 000.parameters directory and run the following commands:

antechamber -i ../zzz.master/2p16_lig_wH.mol2 -fi pdb -o 2p16_ligand_antechamber.mol2 -fo mol2 -at gaff2 -c bcc -rn LIG -nc 0
parmchk2 -i 2p16_ligand_antechamber.mol2 -f mol2 -o 2p16_ligand.am1bcc.frcmod

TLeap

TLeap will generate files describing the topology (prmtop) and initial parameters (rst7) for the protein, ligand, and complex in gas phase, as well as the solvated (wet) complex. Switch to the 001.tleap_build directory and open a file named tleap.build.in and type the following:

#!/bin/bash

###Load Protein force field
source leaprc.protein.ff14SB
###Load GAFF force field (for our ligand)
source leaprc.gaff
###Load TIP3P (water) force field
source leaprc.water.tip3p
####Load Ions frcmod for the tip3p model
loadamberparams frcmod.ionsjc_tip3p
###Needed so we can use igb=8 model
set default PBradii mbondi3


###Load Protein pdb file
rec=loadpdb ../zzz.master/2p16_rec.pdb

#bond rec.7.SG rec.12.SG 
#bond rec.27.SG rec.43.SG 
#bond rec.156.SG rec.170.SG  
#bond rec.181.SG rec.209.SG 
#bond rec.96.SG rec.109.SG
#bond rec.111.SG rec.124.SG
#bond rec.42.SG rec.58.SG
#bond rec.22.SG rec.27.SG

###Load Ligand frcmod/mol2
loadamberparams ../000.parameters/2p16_ligand.am1bcc.frcmod
lig=loadmol2 ../000.parameters/2p16_ligand.am1bcc.mol2

###Create gas-phase complex
gascomplex= combine {rec lig}

###Write gas-phase pdb
savepdb gascomplex 2p16.gas.complex.pdb

###Write gas-phase toplogy and coord files for MMGBSA calc
saveamberparm gascomplex 2p16.gas.complex.prmtop 2p16.gas.complex.rst7
saveamberparm rec 2p16.gas.receptor.prmtop 2p16.gas.receptor.rst7
saveamberparm lig 2p16.gas.ligand.prmtop 2p16.gas.ligand.rst7

###Create solvated complex (albeit redundant)
solvcomplex= combine {rec lig}

###Solvate the system
solvateoct solvcomplex TIP3PBOX 12.0

###Neutralize system (it will add either Na or Cl depending on net charge)
addions solvcomplex Cl- 0
addions solvcomplex Na+ 0

###Write solvated pdb file
savepdb solvcomplex 2p16.wet.complex.pdb

###Check the system
charge solvcomplex
check solvcomplex

###Write Solvated topology and coord file
saveamberparm solvcomplex 2p16.wet.complex.prmtop 2p16.wet.complex.rst7

You can then run this script and copy over the wet complex prmtop and rst7 files to your local machine:

tleap -f tleap.build.in
scp 'username@login.seawulf.stonybrook.edu:path_to_amber_directory/001.tleap_build/2p16.wet.complex.*' .

You can then visualize these files using the following Chimera commands to ensure the previous steps ran properly:

Tools → MD/Ensemble Analysis → MD Movie

Select file 2p16.wet.complex.prmtop for the prmtop section, and then hit add and select the 2p16.wet.complex.rst7 file to visualize the complex. Make sure to visualize the waters and ensure they encapsulate the complex.

Equilibration

Before running AMBER, we will need to energy minimize the system to ensure it is at equilibrium. This helps remove steric clashes seen in higher energy conformations. For this we will need 9 input files, each representing a step in the process. There are two key parameters in each file that are important to note- the restraintmask, which determines which atoms are restrained (i.e. which atoms an energy penalty of movement is applied to), and the restraint_wt which determines the extent of the restraint placed (i.e. the force constant of restraint).

Switch to the 002.equilibration directory, open a file entitled 01.min.mdin and type the following:

 Minmize all the hydrogens
  &cntrl
  imin=1,           ! Minimize the initial structure
  maxcyc=5000,    ! Maximum number of cycles for minimization
  ntb=1,            ! Constant volume
  ntp=0,            ! No pressure scaling
  ntf=1,            ! Complete force evaluation
  ntwx= 1000,       ! Write to trajectory file every ntwx steps
  ntpr= 1000,       ! Print to mdout every ntpr steps
  ntwr= 1000,       ! Write a restart file every ntwr steps
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
  ntr=1,            ! Turn on restraints
  restraintmask=":1-235 & !@H=", ! atoms to be restrained
  restraint_wt=5.0, ! force constant for restraint
  ntxo=1,           ! Write coordinate file in ASCII format
  ioutfm=0,         ! Write trajectory file in ASCII format
  /

Ensure that the number in the restraintmask parameter corresponds to the number of residues in the protein + 1 for the ligand, since we are restraining both for the first few iterations. You can check this number at the bottom of the 2p16.gas.complex.pdb file generated in tleap step.

Now open a file entitled 02.equil.mdin and type in the following:

 MD simualation
  &cntrl
  imin=0,           ! Perform MD
  nstlim=50000      ! Number of MD steps
  ntb=2,            ! Constant Pressure
  ntc=1,            ! No SHAKE on bonds between hydrogens
  dt=0.001,         ! Timestep (ps)
  ntp=1,            ! Isotropic pressure scaling
  barostat=1        ! Berendsen
  taup=0.5          ! Pressure relaxtion time (ps)
  ntf=1,            ! Complete force evaluation
  ntt=3,            ! Langevin thermostat
  gamma_ln=2.0      ! Collision Frequency for thermostat
  ig=-1,            ! Random seed for thermostat
  temp0=298.15      ! Simulation temperature (K)
  ntwx= 1000,       ! Write to trajectory file every ntwx steps
  ntpr= 1000,       ! Print to mdout every ntpr steps
  ntwr= 1000,       ! Write a restart file every ntwr steps
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
  ntr=1,            ! Turn on restraints
  restraintmask=":1-235 & !@H=", ! atoms to be restrained
  restraint_wt=5.0, ! force constant for restraint
  ntxo=1,           ! Write coordinate file in ASCII format
  ioutfm=0,         ! Write trajectory file in ASCII format
  iwrap=1,          ! iwrap is turned on
  /

Open a file entitled 03.min.mdin and type in the following:

  Minmize all the hydrogens
  &cntrl
  imin=1,           ! Minimize the initial structure
  maxcyc=1000,    ! Maximum number of cycles for minimization
  ntb=1,            ! Constant volume
  ntp=0,            ! No pressure scaling
  ntf=1,            ! Complete force evaluation
  ntwx= 1000,       ! Write to trajectory file every ntwx steps
  ntpr= 1000,       ! Print to mdout every ntpr steps
  ntwr= 1000,       ! Write a restart file every ntwr steps
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
  ntr=1,            ! Turn on restraints
  restraintmask=":1-235 & !@H=", ! atoms to be restrained
  restraint_wt=2.0, ! force constant for restraint
  ntxo=1,           ! Write coordinate file in ASCII format
  ioutfm=0,         ! Write trajectory file in ASCII format
  /

Open a file entitled 04.min.mdin and type in the following:

 Minmize all the hydrogens
 &cntrl
 imin=1,           ! Minimize the initial structure
 maxcyc=1000,    ! Maximum number of cycles for minimization
 ntb=1,            ! Constant volume
 ntp=0,            ! No pressure scaling
 ntf=1,            ! Complete force evaluation
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-235 & !@H=", ! atoms to be restrained
 restraint_wt=0.1, ! force constant for restraint
 ntxo=1,           ! Write coordinate file in ASCII format
 ioutfm=0,         ! Write trajectory file in ASCII format
 /

Open a file entitled 05.min.mdin and type in the following:

 Minmize all the hydrogens
 &cntrl
 imin=1,           ! Minimize the initial structure
 maxcyc=1000,    ! Maximum number of cycles for minimization
 ntb=1,            ! Constant volume
 ntp=0,            ! No pressure scaling
 ntf=1,            ! Complete force evaluation
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-235 & !@H=", ! atoms to be restrained
 restraint_wt=0.05, ! force constant for restraint
 ntxo=1,           ! Write coordinate file in ASCII format
 ioutfm=0,         ! Write trajectory file in ASCII format
 /

Open a file entitled 06.equil.mdin and type in the following:

 MD simualation
 &cntrl
 imin=0,           ! Perform MD
 nstlim=50000      ! Number of MD steps
 ntb=2,            ! Constant Pressure
 ntc=1,            ! No SHAKE on bonds between hydrogens
 dt=0.001,         ! Timestep (ps)
 ntp=1,            ! Isotropic pressure scaling
 barostat=1        ! Berendsen
 taup=0.5          ! Pressure relaxtion time (ps)
 ntf=1,            ! Complete force evaluation
 ntt=3,            ! Langevin thermostat
 gamma_ln=2.0      ! Collision Frequency for thermostat
 ig=-1,            ! Random seed for thermostat
 temp0=298.15      ! Simulation temperature (K)
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-235 & !@H=", ! atoms to be restrained
 restraint_wt=1.0, ! force constant for restraint
 ntxo=1,           ! Write coordinate file in ASCII format
 ioutfm=0,         ! Write trajectory file in ASCII format
 iwrap=1,          ! iwrap is turned on
 /

Open a file entitled 07.equil.mdin and type in the following:

 MD simulation
 &cntrl
 imin=0,           ! Perform MD
 nstlim=50000      ! Number of MD steps
 ntx=5,            ! Positions and velocities read formatted
 irest=1,          ! Restart calculation
 ntc=1,            ! No SHAKE on for bonds with hydrogen
 dt=0.001,         ! Timestep (ps)
 ntb=2,            ! Constant Pressure
 ntp=1,            ! Isotropic pressure scaling
 barostat=1        ! Berendsen
 taup=0.5          ! Pressure relaxtion time (ps)
 ntf=1,            ! Complete force evaluation
 ntt=3,            ! Langevin thermostat
 gamma_ln=2.0      ! Collision Frequency for thermostat
 ig=-1,            ! Random seed for thermostat
 temp0=298.15      ! Simulation temperature (K)
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-235 & !@H=", ! atoms to be restrained
 restraint_wt=0.5, ! force constant for restraint
 ntxo=1,           ! Write coordinate file in ASCII format
 ioutfm=0,         ! Write trajectory file in ASCII format
 iwrap=1,          ! iwrap is turned on
 /

Open a file entitled 08.equil.mdin and type in the following:

 MD simulations
 &cntrl
 imin=0,           ! Perform MD
 nstlim=50000      ! Number of MD steps
 ntx=5,            ! Positions and velocities read formatted
 irest=1,          ! Restart calculation
 ntc=1,            ! No SHAKE on for bonds with hydrogen
 dt=0.001,         ! Timestep (ps)
 ntb=2,            ! Constant Pressure
 ntp=1,            ! Isotropic pressure scaling
 barostat=1        ! Berendsen
 taup=0.5          ! Pressure relaxtion time (ps)
 ntf=1,            ! Complete force evaluation
 ntt=3,            ! Langevin thermostat
 gamma_ln=2.0      ! Collision Frequency for thermostat
 ig=-1,            ! Random seed for thermostat
 temp0=298.15      ! Simulation temperature (K)
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-234@CA.C.N", ! atoms to be restrained, only the backbone
 restraint_wt=0.1, ! force constant for restraint
 ntxo=1,           ! Write coordinate file in ASCII format
 ioutfm=0,         ! Write trajectory file in ASCII format
 iwrap=1,          ! iwrap is turned on
 /

Ensure the restraintmask number is 1 lower for this step and the next one, since we're only restraining the protein and not the ligand for the last two rounds of the minimization.

Open a file entitled 09.equil.mdin and type in the following:

 MD simulations
 &cntrl
 imin=0,           ! Perform MD
 nstlim=50000      ! Number of MD steps
 ntx=5,            ! Positions and velocities read formatted
 irest=1,          ! Restart calculation
 ntc=1,            ! No SHAKE on for bonds with hydrogen
 dt=0.001,         ! Timestep (ps)
 ntb=2,            ! Constant Pressure
 ntp=1,            ! Isotropic pressure scaling
 barostat=1        ! Berendsen
 taup=0.5          ! Pressure relaxtion time (ps)
 ntf=1,            ! Complete force evaluation
 ntt=3,            ! Langevin thermostat
 gamma_ln=2.0      ! Collision Frequency for thermostat
 ig=-1,            ! Random seed for thermostat
 temp0=298.15      ! Simulation temperature (K)
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-234@CA.C.N", ! atoms to be restrained
 restraint_wt=0.1, ! force constant for restraint
 ntxo=1,           ! Write coordinate file in ASCII format
 ioutfm=0,         ! Write trajectory file in ASCII format
 iwrap=1,          ! iwrap is turned on
 /

We will need to run the minimization on the cluster using slurm, and so will need to create the following script entitled equil.slurm:

#! /bin/sh
#SBATCH --job-name=sys_equilibration
#SBATCH --output=equil_output.txt
#SBATCH --ntasks-per-node=24
#SBATCH --nodes=6
#SBATCH --time=48:00:00
#SBATCH -p long-24core
module load amber/16

cd $SLURM_SUBMIT_DIR

echo "Started Equilibration on `date` "
do_parallel="mpirun pmemd.MPI"

prmtop="../001.tleap_build/2p16.wet.complex.prmtop"
coords="../001.tleap_build/2p16.wet.complex"


MDINPUTS=(01.min 02.equil 03.min 04.min 05.min 06.equil 07.equil 08.equil 09.equil)

for input in ${MDINPUTS[@]}; do

  $do_parallel -O -i ${input}.mdin -o ${input}.mdout -p $prmtop -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf ${input}.info -r ${input}.rst7
  coords=$input
done

echo "Finished Equilibration on `date` "

Now, we can submit the job using:

sbatch equil.slurm

This should generate a .trj file for each step in the process, which you can scp over and visualize in Chimera just like in the tleap analysis.

Production

After minimizing the system and making sure its at equilibrium, we can finally run a production MD simulation which our analysis will be based on.

Switch to the 003.production directory and create a file called 10.prod.mdin and type in the following:

  MD simulations
  &cntrl
   imin=0,           ! Perform MD
   nstlim=5000000,   ! Number of MD steps
   ntx=5,            ! Positions and velocities read formatted
   irest=1,          ! Restart calculation
   ntc=2,            ! SHAKE on for bonds with hydrogen
   dt=0.002,         ! Timestep (ps)
   ntb=2,            ! Constant Pressure
   ntp=1,            ! Isotropic pressure scaling
   barostat=1        ! Berendsen
   taup=0.5          ! Pressure relaxtion time (ps)
   ntf=2,            ! No force evaluation for bonds with hydrogen
   ntt=3,            ! Langevin thermostat
   gamma_ln=2.0      ! Collision Frequency for thermostat
   ig=-1,            ! Random seed for thermostat
   temp0=298.15      ! Simulation temperature (K)
   ntwx= 2500,       ! Write to trajectory file every ntwx steps
   ntpr= 2500,       ! Print to mdout every ntpr steps
   ntwr= 5000000,    ! Write a restart file every ntwr steps
   cut=8.0,          ! Nonbonded cutoff in Angstroms
   ntr=1,            ! Turn on restraints
   restraintmask=":1-234@CA,C,N", ! atoms to be restrained
   restraint_wt=0.1, ! force constant for restraint
   ntxo=1,           ! Write coordinate file in ASCII format
   ioutfm=0,         ! Write trajectory file in ASCII format
   iwrap=1,          ! iwrap is turned on
  $end

Note that, just like in the last two runs of the equilibration, only the protein is restrained and not the ligand.

We will also need a slurm script to run this job on the cluster, so open a file called prod_slurm.sh and type the following:

 #!/bin/bash
 #
 #SBATCH --job-name=2p16_production
 #SBATCH --output=2p16_production.txt
 #SBATCH --ntasks-per-node=28
 #SBATCH --nodes=4
 #SBATCH --time=48:00:00
 #SBATCH -p long-28core 

 module load amber/16

 do_parallel="mpirun pmemd.MPI"

 prmtop="../003.tleap/2p16.wet.complex.prmtop"
 coords="../004.equil/09.equil"

 MDINPUTS=(10.prod)

 for input in ${MDINPUTS[@]}; do
   $do_parallel -O -i ${input}.mdin -o ${input}.mdout -p $prmtop -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf ${input}.info -r ${input}.rst7
   coords=$input
 done

Submit this job using sbatch and visualize the resulting trj file in Chimera.

Analysis

We can now use the finalized production run to analyze stability, hydrogen bond formation, and energetics of the system. Switch to the 004.analysis directory and create the following sub directories:

001.rmsd
002.hbond
003.mmgbsa

RMSD

Switch to the 001.rmsd directory, open a file entitled cpptraj.strip.wat.in and type in the following:

#!/usr/bin/sh
parm ../../001.tleap_build/2p16.wet.complex.prmtop
#read in trajectory
trajin ../../003.production/10.prod.trj
#read in reference
reference ../../001.tleap_build/2p16.wet.complex.rst7
#compute rmsd and align CA to the crystal structure
rmsd rms1 reference :1-235@CA
#strip Solvent
strip :WAT:Na+:Cl-
#create gas-phase trajectory
trajout 2p16.stripfit.restrained.gas.trj nobox

This file will delete the solvent and charges from the production run using the following command:

cpptraj -i cpptraj.strip.wat.in

We can now open an input file entitled cpptraj.rmsd.lig.in that calculates ligand RMSD from the .trj file generated in the previous step:

#!/usr/bin/sh
#trajin the trajectory
trajin 2p16.stripfit.restrained.gas.trj
#read in the reference
reference ../../001.tleap_build/2p16.gas.complex.rst7
#compute the RMSD (do not fit the internal geometries first, included rigid body motions
#and convert the frames to ns (framenum*.005)
rmsd rms1 ":235&!(@H=)" nofit mass out 2p16.lig.restrained.rmsd.nofit.dat time .005
#histogram the nofit rmsd
histogram rms1,*,*,.1,* norm out 2p16.lig.restrained.rmsd.nofit.histogram.dat

We can also open an input file entitled cpptraj.rmsd.protein.in that calculates protein RMSD:

#!/usr/bin/sh
#trajin the trajectory
trajin 2p16.stripfit.restrained.gas.trj
#read in the reference
reference ../../001.tleap_build/2p16.gas.complex.rst7
#compute the RMSD (do not fit the internal geometries first, included rigid body motions
#and convert the frames to ns (framenum*.005)
rmsd rms1 ":1-234&!(@H=)" nofit mass out 2p16.lig.restrained.rmsd.nofit.dat time .005
#histogram the nofit rmsd
histogram rms1,*,*,.1,* norm out 2p16.protein.restrained.rmsd.nofit.histogram.dat

Note the differences in residue numbers chosen

These scripts can be run using the following commands:

cpptraj -p ../001.tleap_build/2p16.gas.complex.prmtop -i cpptraj.rmsd.lig.in
cpptraj -p ../001.tleap_build/2p16.gas.complex.prmtop -i cpptraj.rmsd.protein.in

This should generate the following data files which can be plotted using python:

2p16.lig.restrained.rmsd.nofit.dat
2p16.lig.restrained.rmsd.nofit.histogram.dat
2p16.protein.restrained.rmsd.nofit.dat
2p16.protein.restrained.rmsd.nofit.histogram.dat

Typically, an RMSD < 2 angstroms is considered a successful run.

H-Bonds

Switch to the 002.hbond directory, open a file entitled cpptraj.hbond.in and type in the following:

#!/usr/bin/sh
#read in trajectory
trajin ../../003.production/10.prod.trj
#wrap everything into one periodic cell
#autoimage
#compute intra and water mediated hydrogen bonds
hbond hb1 :1-235 out 2p16.hbond.out avgout 2p16.hbond.avg.dat solventdonor :WAT solventacceptor :WAT@O nointramol bridgeout 2p16.bridge-water.dat dist 3.0 angle 140

This will need to be run using slurm on the cluster, so create a file called hbond_slurm.sh and type in the following:

#! /bin/sh
#SBATCH --job-name=hbond_anal
#SBATCH --ntasks-per-node=24
#SBATCH --nodes=1
#SBATCH --time=15:00:00
#SBATCH -p long-24core

cd $SLURM_SUBMIT_DIR

echo "Started Analysis on `date` "
#do_parallel="mpirun pmemd.MPI"

 cpptraj -p ../../001.tleap_build/2p16.wet.complex.prmtop -i cpptraj.hbond.in

echo "Finished Analysis on `date` "

Submit it using sbatch and examine the resulting dat files to see which h-bonds played a major role across frames.

mmGBSA

Switch to the 003.mmgbsa directory, open a file entitled mmgbsa_slurm.sh and type in the following:

#! /bin/sh
#SBATCH --job-name=production_unrest
#SBATCH --ntasks-per-node=24
#SBATCH --nodes=1
#SBATCH --time=24:00:00
#SBATCH -p long-24core

cd $SLURM_SUBMIT_DIR

#Define topology files 
solv_prmtop="../../001.tleap_build/2p16.wet.complex.prmtop"
complex_prmtop="../../001.tleap_build/2p16.gas.complex.prmtop"
receptor_prmtop="../../001.tleap_build/2p16.gas.receptor.prmtop"
ligand_prmtop="../../001.tleap_build/2p16.gas.ligand.prmtop"
trajectory="../../003.production/001.restrained/10.prod.trj"

#create mmgbsa input file
cat >mmgbsa.in<<EOF
mmgbsa 2p16 analysis
&general
  interval=1, netcdf=1,
  keep_files=0,

/
&gb
  igb=8,
  saltcon=0.0, surften=0.0072,
  surfoff=0.0, molsurf=0,
/
&nmode
  drms=0.001, maxcyc=10000,
  nminterval=250, nmendframe=2000,
  nmode_igb=1, 
/
EOF

MMPBSA.py -O -i mmgbsa.in \
           -o  2p16.mmgbsa.results.dat \
           -eo 2p16.mmgbsa.per-frame.dat \
           -sp ${solv_prmtop} \
           -cp ${complex_prmtop} \
           -rp ${receptor_prmtop} \
           -lp ${ligand_prmtop} \
            -y ${trajectory}

Run it using sbatch and examine or plot the resulting dat files to see the overall as well as per-frame calculated binding energy.