2018 AMBER tutorial with 1c87

From Rizzo_Lab
Revision as of 15:50, 28 March 2018 by AMS536 2018 02 (talk | contribs) (IV. Simulation Analysis)
Jump to: navigation, search

In this tutorial, we will learn how to run a MD simulation of a protein-ligand complex and learn how to calculate structural fluctuations (with RMSD) and free energies of binding (MM-GBSA).

I. Prepare the files and directories

We can create the directories to make all the files organized.

 ~username/AMS536/2018/Amber_Tutorial/000.parameters
                                                                    /001.tleap_build
                                                                    /002.equilibration       
                                                                    /003.production
                                                                    /004.analysis
                                                                    /zzz.master

Once the directories are created, we are going to convert 1c87_lig_withH_charged.mol2 file and 1c87_rec_withH_charged.mol2 file to pdb file in chimera, then copy the converted files into zzz.master directory.

II. Structural Preparation

Before we submit our jobs for equilibration or production, we must parametrize the ligands by using antechamber and check for any missing force field parameters. Make sure we have the proper partial charge. Since there are two acids in the 1c87 structure, -nc will be -2.

cd 000.parameters
# Parameterize the ligand
antechamber -i ../zzz.master/2nnq.lig.withH.charged.for_antechamber.pdb -fi pdb -o 2nnq_lig.am1bcc.mol2 -fo mol2 -at gaff2 -c bcc -rn LIG -nc -2
# Check for missing FF parameters
parmchk2 -i 2nnq_lig.am1bcc.mol2 -f mol2 -o 2nnq_lig.am1bcc.frcmod

Then, we will create the topology and coordinate files for the simulation through tleap. Create a file called tleap.in with the following:

#!/usr/bin/sh

###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/rec_withH_1c87.pdb

###Load Ligand frcmod/mol2
loadamberparams ../000.parameters/1c87_lig.am1bcc.frcmod
lig=loadmol2 ../000.parameters/1c87_lig.am1bcc.mol2

###Create gas-phase complex
gascomplex= combine {rec lig}
###Write gas-phase pdb
savepdb gascomplex 1c87.gas.complex.pdb
 
###Write gas-phase toplogy and coord files for MMGBSA calc 
saveamberparm gascomplex 1c87.gas.complex.prmtop 1c87.gas.complex.rst7
saveamberparm rec 1c87.gas.receptor.prmtop 1c87.gas.receptor.rst7
saveamberparm lig 1c87.gas.ligand.prmtop 1c87.gas.ligand.rst7

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

###Solvate the systemsolvateoct 
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 1c87.wet.complex.pdb

###Check the system
charge solvcomplex
check solvcomplex
###Write Solvated topology and coord file
saveamberparm solvcomplex 1c87.wet.complex.prmtop 1c87.wet.complex.rst7

Create Amber topology and coordinates files for the MD simulation.

tleap -f tleap.in

III. Equilibration

We can now create the input files for the equilibration of the system in 002.equilbration directory. We are going to run minimization in order to remove steric clashes and allow the receptor and ligand to equilibrate into a low energy state (MD simulation).

01.min.mdin:

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-298 & !@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
/

02.equil.mdin:

 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-298 & !@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
 /

03.min.mdin:

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-298 & !@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
 /

04.min.mdin:

 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-298 & !@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
 /

05.min.mdin:

  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-298 & !@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
 /

06.equil.mdin:

 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-298 & !@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
 /

07.equil.mdin:

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-298 & !@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
 /

08.equil.mdin:

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-298 & !@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
 iwrap=1,          ! iwrap is turned on
 /

09.equil.mdin:

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-298 & !@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
 iwrap=1,          ! iwrap is turned on
 /

Create a submission script md.equilibration.sh in order to start equilibration. We are going to run it for 24 hours in order to have enough time to complete all equilibration jobs.

vi md.equilibration.sh
#!/bin/sh
#PBS -N 1c87_equilibration
#PBS -l walltime=24:00:00
#PBS -l nodes=2:ppn=28
#PBS -j oe
#PBS -q long

cd $PBS_O_WORKDIR

echo "Started Equilibration on `date` "
do_parallel="sander"

prmtop="../001.tleap_build/1c87.wet.complex.prmtop"
coords="../001.tleap_build/1c87.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` "

Submit the job

qsub md.equilibration.sh

Once submitted the job, info, rst7, trj files will be created for each input file. We can use VMD to open up the prmtop file and trajectory file. Determine file type: AMBER coordinates.

Minimization


IV. Production

Move to 003.production directory and open up 001.restrained directory. We are going to create the file for running production. Create the input file for MD simulation of production.

vi 10.prod.mdin 
 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-298@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 are going to create the bash script.

vi prod.sh
#!/bin/sh
#PBS -N 1c87_production
#PBS -l walltime=24:00:00
#PBS -l nodes=2:ppn=28
#PBS -j oe
#PBS -q long

cd $PBS_O_WORKDIR

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

prmtop="/gpfs/projects/AMS536/2018/seung/AMBER_directory_arch/001.tleap_build/1c87.wet.complex.prmtop"
coords="../../002.equilbration/07.equil"


MDINPUTS=(10.prod)

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

 $do_cuda -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 Production on `date` "

Submit the job.

qsub prod.sh

V. Simulation Analysis

After completing production simulations, we can now start analyzing and apply post processing. We can create three directories:

mkdir 001.rmsd  002.hbond  003.mmgbsa

We will now remove the water molecules and create a trajectory aligned to the CA atoms (backbone) of the crystal structure. In order to do this, we can create "cpptraj.strip.wat.in" file :

#!/usr/bin/sh
#read in trajectory
trajin /gpfs/projects/AMS536/2018/seung/AMBER_directory_arch/003.production/001.restrained/10.prod.trj 
#read in reference
reference /gpfs/projects/AMS536/2018/seung/AMBER_directory_arch/001.tleap_build/1c87.wet.complex.rst7 
#compute rmsd and align CA to the crystal structure
rmsd rms1 reference :1-298@CA
#strip Solvent
strip :WAT:Na+:Cl-
#create gas-phase trajectory
trajout 1c87.stripfit.restrained.gas.trj nobox

and we are going to run it:

 cpptraj -p ../../001.tleap_build/1c87_sunitinib.wet.complex.prmtop -i cpptraj.strip.wat.in

We can take a look at the trajectory file in VMD.

Next is to compute ligand rmsd vs. the crystal structure pose and histogram the RMSD.

#!/usr/bin/sh
#trajin the trajectory
trajin 1c87.stripfit.restrained.gas.trj
#read in the reference
reference ../../001.tleap_build/1c87.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 ":298&!(@H=)" nofit mass out 1c87.lig.restrained.rmsd.nofit.dat time .005 
#histogram the nofit rmsd
histogram rms1,*,*,.1,* norm out 1c87.lig.restrained.rmsd.nofit.histogram.dat

and again, we will run it:

   cpptraj -p ../../001.tleap_build/1c87_sunitinib.gas.complex.prmtop -i cpptraj.rmsd.lig.in