View source for 2018 AMBER tutorial with 1c87

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. 
 vi md.equilibration.sh

 #!/bin/sh
 #PBS -N 1c87_equilibration
 #PBS -l walltime=12: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