Difference between revisions of "2019 AMBER tutorial with PDBID 2BXF"

From Rizzo_Lab
Jump to: navigation, search
(TLeap)
(H Bond)
 
(13 intermediate revisions by the same user not shown)
Line 29: Line 29:
 
Create tleap.build.in file
 
Create tleap.build.in file
  
  #!/usr/bin/sh
+
#!/usr/bin/sh
  
 
  ###Load Protein force field
 
  ###Load Protein force field
Line 41: Line 41:
 
  ###Needed so we can use igb=8 model
 
  ###Needed so we can use igb=8 model
 
  set default PBradii mbondi3
 
  set default PBradii mbondi3
 
 
 
  ###Load Protein pdb file
 
  ###Load Protein pdb file
 
  rec=loadpdb ../zzz.master/backup.pdb
 
  rec=loadpdb ../zzz.master/backup.pdb
Line 62: Line 60:
 
  bond rec.554.SG rec.563.SG
 
  bond rec.554.SG rec.563.SG
 
  bond rec.356.SG rec.365.SG
 
  bond rec.356.SG rec.365.SG
 
 
  ###Load Ligand frcmod/mol2
 
  ###Load Ligand frcmod/mol2
 
  loadamberparams ../000.parameters/2BXF_lig.am1bcc.frcmod
 
  loadamberparams ../000.parameters/2BXF_lig.am1bcc.frcmod
 
  lig=loadmol2 ../000.parameters/2BXF_lig.am1bcc.mol2
 
  lig=loadmol2 ../000.parameters/2BXF_lig.am1bcc.mol2
 
 
  ###Create gas-phase complex
 
  ###Create gas-phase complex
 
  gascomplex= combine {rec lig}
 
  gascomplex= combine {rec lig}
 
 
  ###Write gas-phase pdb
 
  ###Write gas-phase pdb
 
  savepdb gascomplex 2BXF.gas.complex.pdb
 
  savepdb gascomplex 2BXF.gas.complex.pdb
 
 
  ###Write gas-phase toplogy and coord files for MMGBSA calc
 
  ###Write gas-phase toplogy and coord files for MMGBSA calc
 
  saveamberparm gascomplex 2BXF.gas.complex.prmtop 2BXF.gas.complex.rst7
 
  saveamberparm gascomplex 2BXF.gas.complex.prmtop 2BXF.gas.complex.rst7
 
  saveamberparm rec 2BXF.gas.receptor.prmtop 2BXF.gas.receptor.rst7
 
  saveamberparm rec 2BXF.gas.receptor.prmtop 2BXF.gas.receptor.rst7
 
  saveamberparm lig 2BXF.gas.ligand.prmtop 2BXF.gas.ligand.rst7
 
  saveamberparm lig 2BXF.gas.ligand.prmtop 2BXF.gas.ligand.rst7
 
 
  ###Create solvated complex (albeit redundant)
 
  ###Create solvated complex (albeit redundant)
 
  solvcomplex= combine {rec lig}
 
  solvcomplex= combine {rec lig}
 
 
  ###Solvate the system
 
  ###Solvate the system
 
  solvateoct solvcomplex TIP3PBOX 12.0
 
  solvateoct solvcomplex TIP3PBOX 12.0
 
 
  ###Neutralize system (it will add either Na or Cl depending on net charge)
 
  ###Neutralize system (it will add either Na or Cl depending on net charge)
 
  addions solvcomplex Cl- 0
 
  addions solvcomplex Cl- 0
 
  addions solvcomplex Na+ 0
 
  addions solvcomplex Na+ 0
 
 
  ###Write solvated pdb file
 
  ###Write solvated pdb file
 
  savepdb solvcomplex 2BXF.wet.complex.pdb
 
  savepdb solvcomplex 2BXF.wet.complex.pdb
 
 
  ###Check the system
 
  ###Check the system
 
  charge solvcomplex
 
  charge solvcomplex
 
  check solvcomplex
 
  check solvcomplex
 
 
  ###Write Solvated topology and coord file
 
  ###Write Solvated topology and coord file
 
  saveamberparm solvcomplex 2BXF.wet.complex.prmtop 2BXF.wet.complex.rst7
 
  saveamberparm solvcomplex 2BXF.wet.complex.prmtop 2BXF.wet.complex.rst7
     
 
  
Create Amber topology and coordinates files for the MD simulation
 
  
 
  tleap -f tleap.build.in
 
  tleap -f tleap.build.in
  
  #!/usr/bin/sh
+
= Equilibration =
 +
 
 +
Create the following input files for the equilibration of the system.
 +
  vim 01.min.mdin
 +
 
 +
Minmize all the hydrogens
 +
&cntrl
 +
  imin=1,          ! Minimize the initial structure
 +
  ntmin=2,        ! Use steepest descent Ryota Added
 +
  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-579 & !@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
 +
/
 +
 
 +
vim 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-579 & !@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
 +
/
 +
 
 +
vim 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-579 & !@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
 +
/
 +
 
 +
vim 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-579 & !@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
 +
/
 +
 
 +
vim 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-579 & !@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
 +
/
  
  ###Load Protein force field
+
  vim 06.equil.mdin
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
 
  
 +
  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-579 & !@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
 +
/
  
  ###Load Protein pdb file
+
  vim 07.equil.mdin
rec=loadpdb ../zzz.master/backup.pdb
 
bond rec.49.SG  rec.58.SG
 
bond rec.241.SG rec.249.SG
 
bond rec.71.SG  rec.87.SG
 
bond rec.86.SG  rec.97.SG
 
bond rec.120.SG rec.165.SG
 
bond rec.164.SG rec.173.SG
 
bond rec.510.SG rec.555.SG
 
bond rec.472.SG rec.482.SG
 
bond rec.457.SG rec.473.SG
 
bond rec.196.SG rec.242.SG
 
bond rec.261.SG rec.275.SG
 
bond rec.312.SG rec.357.SG
 
bond rec.274.SG rec.285.SG
 
bond rec.388.SG rec.434.SG
 
bond rec.433.SG rec.444.SG
 
bond rec.554.SG rec.563.SG
 
bond rec.356.SG rec.365.SG
 
  
  ###Load Ligand frcmod/mol2
+
  MD simulation
  loadamberparams ../000.parameters/2BXF_lig.am1bcc.frcmod
+
  &cntrl
lig=loadmol2 ../000.parameters/2BXF_lig.am1bcc.mol2
+
  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-579 & !@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
 +
/
  
  ###Create gas-phase complex
+
  vim 08.equil.mdin
gascomplex= combine {rec lig}
 
  
  ###Write gas-phase pdb
+
  MD simulations
  savepdb gascomplex 2BXF.gas.complex.pdb
+
&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-578@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
 +
/
  
  ###Write gas-phase toplogy and coord files for MMGBSA calc
+
  vim 09.equil.mdin
saveamberparm gascomplex 2BXF.gas.complex.prmtop 2BXF.gas.complex.rst7
+
   
saveamberparm rec 2BXF.gas.receptor.prmtop 2BXF.gas.receptor.rst7
 
  saveamberparm lig 2BXF.gas.ligand.prmtop 2BXF.gas.ligand.rst7
 
  
  ###Create solvated complex (albeit redundant)
+
  MD simulations
  solvcomplex= combine {rec lig}
+
&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-578@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
 +
/
  
###Solvate the system
+
Create a submission script that utilizes, above input files to equilibrate the biological system.
  solvateoct solvcomplex TIP3PBOX 12.0
+
  vim mdequilibration.sh
  
  ###Neutralize system (it will add either Na or Cl depending on net charge)
+
  #!/bin/sh                                                                                                                                                                                                                                   
  addions solvcomplex Cl- 0
+
#PBS -N 2BXF_equilibration                                                                                                                                                                                                                   
  addions solvcomplex Na+ 0
+
#PBS -l walltime=04: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/2BXF.wet.complex.prmtop"
 +
coords="../001.tleap_build/2BXF.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` "
  
###Write solvated pdb file
+
Submit the job
savepdb solvcomplex 2BXF.wet.complex.pdb
 
  
  ###Check the system
+
  qsub md.equilibration.sh
charge solvcomplex
 
check solvcomplex
 
  
###Write Solvated topology and coord file
+
= Production =
saveamberparm solvcomplex 2BXF.wet.complex.prmtop 2BXF.wet.complex.rst7
 
  
#!/usr/bin/sh
+
Move into 003.production
  
###Load Protein force field
+
Move into 001.restrained
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
 
  
 +
Create 10.prod.mdin
  
  ###Load Protein pdb file
+
  MD simulations
rec=loadpdb ../zzz.master/backup.pdb
+
  &cntrl
bond rec.49.SG  rec.58.SG
+
  imin=0,          ! Perform MD
bond rec.241.SG rec.249.SG
+
  nstlim=5000000,  ! Number of MD steps
bond rec.71.SG  rec.87.SG
+
  ntx=5,            ! Positions and velocities read formatted
bond rec.86.SG  rec.97.SG
+
  irest=1,          ! Restart calculation
bond rec.120.SG rec.165.SG
+
  ntc=2,            ! SHAKE on for bonds with hydrogen
bond rec.164.SG rec.173.SG
+
  dt=0.002,        ! Timestep (ps)
bond rec.510.SG rec.555.SG
+
  ntb=2,            ! Constant Pressure
bond rec.472.SG rec.482.SG
+
  ntp=1,            ! Isotropic pressure scaling
bond rec.457.SG rec.473.SG
+
  barostat=1        ! Berendsen
bond rec.196.SG rec.242.SG
+
  taup=0.5          ! Pressure relaxtion time (ps)
bond rec.261.SG rec.275.SG
+
  ntf=2,            ! No force evaluation for bonds with hydrogen
bond rec.312.SG rec.357.SG
+
  ntt=3,            ! Langevin thermostat
bond rec.274.SG rec.285.SG
+
  gamma_ln=2.0      ! Collision Frequency for thermostat
bond rec.388.SG rec.434.SG
+
  ig=-1,            ! Random seed for thermostat
bond rec.433.SG rec.444.SG
+
  temp0=298.15      ! Simulation temperature (K)
bond rec.554.SG rec.563.SG
+
  ntwx= 2500,      ! Write to trajectory file every ntwx steps
bond rec.356.SG rec.365.SG
+
  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-579@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
  
  ###Load Ligand frcmod/mol2
+
Create prod.sh
  loadamberparams ../000.parameters/2BXF_lig.am1bcc.frcmod
+
 
  lig=loadmol2 ../000.parameters/2BXF_lig.am1bcc.mol2
+
  #!/bin/sh                                                                                                                                                                                                                                   
 +
#PBS -N 2BXF_production                                                                                                                                                                                                                     
 +
#PBS -l walltime=24:00:00                                                                                                                                                                                                                   
 +
#PBS -l nodes=2:ppn=24                                                                                                                                                                                                                       
 +
#PBS -j oe                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 +
  #PBS -o production.output
 +
#PBS -q long                                                                                                                                                                                                                   
 +
 +
cd $PBS_O_WORKDIR
 +
 +
echo "Started Production on `date` "
 +
#do_parallel="mpirun pmemd.MPI"                                                                                                                                                                                                             
 +
do_cuda="mpirun pmemd.MPI"
 +
 +
  prmtop="../../001.tleap_build/2BXF.wet.complex.prmtop"
 +
  coords="../../002.equilbration/09.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` "
  
###Create gas-phase complex
 
gascomplex= combine {rec lig}
 
  
###Write gas-phase pdb
+
Submit the job
  savepdb gascomplex 2BXF.gas.complex.pdb
+
  qsub prod.sh
  
###Write gas-phase toplogy and coord files for MMGBSA calc
+
Move into 002.unrestrained
saveamberparm gascomplex 2BXF.gas.complex.prmtop 2BXF.gas.complex.rst7
 
saveamberparm rec 2BXF.gas.receptor.prmtop 2BXF.gas.receptor.rst7
 
saveamberparm lig 2BXF.gas.ligand.prmtop 2BXF.gas.ligand.rst7
 
  
###Create solvated complex (albeit redundant)
+
Create 10.prod.mdin
solvcomplex= combine {rec lig}
 
  
  ###Solvate the system
+
  MD simulations
  solvateoct solvcomplex TIP3PBOX 12.0
+
  &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
 +
  ntxo=1,          ! Write coordinate file in ASCII format
 +
  ioutfm=0,        ! Write trajectory file in ASCII format
 +
  iwrap=1,          ! iwrap is turned on
 +
/
  
###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
+
Create prod.sh
  savepdb solvcomplex 2BXF.wet.complex.pdb
+
 
 +
  #!/bin/sh                                                                                                                                                                                                                                   
 +
#PBS -N 2BXF_production                                                                                                                                                                                                                     
 +
#PBS -l walltime=24:00:00                                                                                                                                                                                                                   
 +
#PBS -l nodes=2:ppn=24                                                                                                                                                                                                                       
 +
#PBS -j oe                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
 +
  #PBS -o production.output
 +
#PBS -q long                                                                                                                                                                                                                   
 +
 +
cd $PBS_O_WORKDIR
 +
 +
echo "Started Production on `date` "
 +
#do_parallel="mpirun pmemd.MPI"                                                                                                                                                                                                             
 +
do_cuda="mpirun pmemd.MPI"
 +
 +
  prmtop="../../001.tleap_build/2BXF.wet.complex.prmtop"
 +
coords="../../002.equilbration/09.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` "
  
###Check the system
 
charge solvcomplex
 
check solvcomplex
 
  
###Write Solvated topology and coord file
+
Submit the job
  saveamberparm solvcomplex 2BXF.wet.complex.prmtop 2BXF.wet.complex.rst7
+
  qsub prod.sh
  
 
= Analysis =
 
= Analysis =
Line 255: Line 510:
 
  reference ../../001.tleap_build/2BXF.wet.complex.rst7
 
  reference ../../001.tleap_build/2BXF.wet.complex.rst7
 
  #compute rmsd and align CA to the crystal structure                         
 
  #compute rmsd and align CA to the crystal structure                         
  rmsd rms1 reference :1-131@CA
+
  rmsd rms1 reference :1-578@CA
 
  #strip Solvent                                                               
 
  #strip Solvent                                                               
 
  strip :WAT:Na+:Cl-
 
  strip :WAT:Na+:Cl-
Line 272: Line 527:
 
  #compute the RMSD (do not fit the internal geometries first, included rigid body motions                                                               
 
  #compute the RMSD (do not fit the internal geometries first, included rigid body motions                                                               
 
  #and convert the frames to ns (framenum*.005)                               
 
  #and convert the frames to ns (framenum*.005)                               
  rmsd rms1 ":132&!(@H=)" nofit mass out 2BXF.lig.restrained.rmsd.nofit.dat time .005
+
  rmsd rms1 ":579&!(@H=)" nofit mass out 2BXF.lig.restrained.rmsd.nofit.dat time .005
 
  #histogram the nofit rmsd                                                   
 
  #histogram the nofit rmsd                                                   
 
  histogram rms1,*,*,.1,* norm out 2BXF.lig.restrained.rmsd.nofit.histogram.dat
 
  histogram rms1,*,*,.1,* norm out 2BXF.lig.restrained.rmsd.nofit.histogram.dat
Line 288: Line 543:
 
  #compute the RMSD (do not fit the internal geometries first, included rigid body motions         
 
  #compute the RMSD (do not fit the internal geometries first, included rigid body motions         
 
  #and convert the frames to ns (framenum*.005)                                                     
 
  #and convert the frames to ns (framenum*.005)                                                     
  rmsd rms1 ":1-131&!(@H=)" nofit mass out 2BXF.rec.restrained.rmsd.nofit.dat time .005
+
  rmsd rms1 ":1-578&!(@H=)" nofit mass out 2BXF.rec.restrained.rmsd.nofit.dat time .005
 
  #histogram the nofit rmsd                                                                         
 
  #histogram the nofit rmsd                                                                         
 
  histogram rms1,*,*,.1,* norm out 2BXF.rec.restrained.rmsd.nofit.histogram.dat
 
  histogram rms1,*,*,.1,* norm out 2BXF.rec.restrained.rmsd.nofit.histogram.dat
Line 306: Line 561:
 
  #autoimage                                                                                                                                 
 
  #autoimage                                                                                                                                 
 
  #compute intra and water mediated hydrogen bonds                                                                                           
 
  #compute intra and water mediated hydrogen bonds                                                                                           
  hbond hb1 :1-288 out 2BXF_sunitinib.hbond.out avgout 2BXF_sunitinib.hbond.avg.dat solventdonor :WAT solventacceptor :WAT@O  
+
  hbond hb1 :1-579 out 2BXF_sunitinib.hbond.out avgout 2BXF_sunitinib.hbond.avg.dat solventdonor :WAT solventacceptor :WAT@O  
 
  nointramol brid\
 
  nointramol brid\
 
  geout 2BXF_sunitinib.bridge-water.dat dist 3.0 angle 140
 
  geout 2BXF_sunitinib.bridge-water.dat dist 3.0 angle 140

Latest revision as of 15:14, 19 April 2019

2BXF with an explicit solvent model

Prepare the files

Convert 2BXF.lig.withH.charged.mol2 to pdb in chimera

Convert 2BXF.rec.withH.charged.mol2 to pdb in chimera

Copy into zzz.master

Parameters

The system we are working on has two main components (Protein receptor & ligand). The usual forcefield "ff14SB" contains all the parameters needed for calculations of the protein. However, the ligand is a non-protein component. Therefore, ff14SB forcefield does not contain the parameters needed for the calculations regarding the ligand. Therefore, we need to generate parameters needed for the ligand. The following steps will be taken using antechamber in order to generate the ligand parameters.

Move into 000.programs

Paramaterize the ligand

antechamber -i ../zzz.master/2BXF.lig.withH.charged.pdb -fi pdb -o 2BXF_lig.am1bcc.mol2 -fo mol2 -at gaff2 -c bcc -rn LIG -nc -1

Check for missing force field parameters

parmchk2 -i 2BXF_lig.am1bcc.mol2 -f mol2 -o 2BXF_lig.am1bcc.frcmod

TLeap

Move into 001.tleap_build

Create tleap.build.in file

#!/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/backup.pdb
bond rec.49.SG  rec.58.SG
bond rec.241.SG rec.249.SG
bond rec.71.SG  rec.87.SG
bond rec.86.SG  rec.97.SG
bond rec.120.SG rec.165.SG
bond rec.164.SG rec.173.SG
bond rec.510.SG rec.555.SG
bond rec.472.SG rec.482.SG
bond rec.457.SG rec.473.SG
bond rec.196.SG rec.242.SG
bond rec.261.SG rec.275.SG
bond rec.312.SG rec.357.SG
bond rec.274.SG rec.285.SG
bond rec.388.SG rec.434.SG
bond rec.433.SG rec.444.SG
bond rec.554.SG rec.563.SG
bond rec.356.SG rec.365.SG
###Load Ligand frcmod/mol2
loadamberparams ../000.parameters/2BXF_lig.am1bcc.frcmod
lig=loadmol2 ../000.parameters/2BXF_lig.am1bcc.mol2
###Create gas-phase complex
gascomplex= combine {rec lig}
###Write gas-phase pdb
savepdb gascomplex 2BXF.gas.complex.pdb
###Write gas-phase toplogy and coord files for MMGBSA calc
saveamberparm gascomplex 2BXF.gas.complex.prmtop 2BXF.gas.complex.rst7
saveamberparm rec 2BXF.gas.receptor.prmtop 2BXF.gas.receptor.rst7
saveamberparm lig 2BXF.gas.ligand.prmtop 2BXF.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 2BXF.wet.complex.pdb
###Check the system
charge solvcomplex
check solvcomplex
###Write Solvated topology and coord file
saveamberparm solvcomplex 2BXF.wet.complex.prmtop 2BXF.wet.complex.rst7


tleap -f tleap.build.in

Equilibration

Create the following input files for the equilibration of the system.

vim 01.min.mdin
Minmize all the hydrogens
&cntrl
 imin=1,           ! Minimize the initial structure
 ntmin=2,         ! Use steepest descent Ryota Added
 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-579 & !@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
/
vim 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-579 & !@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
/
vim 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-579 & !@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
/
vim 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-579 & !@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
/
vim 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-579 & !@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
/
vim 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-579 & !@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
/
vim 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-579 & !@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
/
vim 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-578@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
/
vim 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-578@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
/

Create a submission script that utilizes, above input files to equilibrate the biological system.

vim mdequilibration.sh
#!/bin/sh                                                                                                                                                                                                                                     
#PBS -N 2BXF_equilibration                                                                                                                                                                                                                    
#PBS -l walltime=04: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/2BXF.wet.complex.prmtop"
coords="../001.tleap_build/2BXF.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

Production

Move into 003.production

Move into 001.restrained

Create 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-579@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

Create prod.sh

 #!/bin/sh                                                                                                                                                                                                                                     
#PBS -N 2BXF_production                                                                                                                                                                                                                       
#PBS -l walltime=24:00:00                                                                                                                                                                                                                     
#PBS -l nodes=2:ppn=24                                                                                                                                                                                                                        
#PBS -j oe                                                                                                                                                                                                                                                                                                                                                                                                                                                                      
#PBS -o production.output
#PBS -q long                                                                                                                                                                                                                     

cd $PBS_O_WORKDIR

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

prmtop="../../001.tleap_build/2BXF.wet.complex.prmtop"
coords="../../002.equilbration/09.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

Move into 002.unrestrained

Create 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
 ntxo=1,           ! Write coordinate file in ASCII format
 ioutfm=0,         ! Write trajectory file in ASCII format
 iwrap=1,          ! iwrap is turned on
/


Create prod.sh

 #!/bin/sh                                                                                                                                                                                                                                     
#PBS -N 2BXF_production                                                                                                                                                                                                                       
#PBS -l walltime=24:00:00                                                                                                                                                                                                                     
#PBS -l nodes=2:ppn=24                                                                                                                                                                                                                        
#PBS -j oe                                                                                                                                                                                                                                                                                                                                                                                                                                                                      
#PBS -o production.output
#PBS -q long                                                                                                                                                                                                                     

cd $PBS_O_WORKDIR

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

prmtop="../../001.tleap_build/2BXF.wet.complex.prmtop"
coords="../../002.equilbration/09.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

Analysis

RMSD

How far our ligand moved and how far our receptor moved

Create cpptraj.strip.wat.in

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

Run cpptraj

cpptraj -i cpptraj.strip.wat.in

Create cpptraj.rmsd.lig.in

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


Run cpptraj

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

Create cpptraj.rmsd.rec.in

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

Run cpptraj

cpptraj -p ../../001.tleap_build/2BXF.gas.complex.prmtop -i cpptraj.rmsd.rec.in

H Bond

Hydrogen bonding

Create cpptraj.hbond.in

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

Run cpptraj

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

MMGBSA

This will analyze how strongly our small molecule binds to our receptor.

Create mmgbsa.in

mmgbsa 2BXF 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,                                                                                                                              
/                                                                                                                                           


Create mmgbsa.sh

#!/bin/bash                                                                                                                                
#PBS -l walltime=35:00:00                                                                                                                  
#PBS -l nodes=1:ppn=24                                                                                                                     
#PBS -N 4qmz_mmgbsa                                                                                                                        
#PBS -V                                                                                                                                    
#PBS -j oe                                                                                                                                 
#PBS -q long-24core                                                                                                                         

cd $PBS_O_WORKDIR 

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

#create mmgbsa input file                                                                                                                  
cat >mmgbsa.in<<EOF                                                                                                                        
mmgbsa HIVgp41 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  2BXF.mmgbsa.results.dat \
           -eo 2BXF.mmgbsa.per-frame.dat \
           -sp ${solv_prmtop} \
           -cp ${complex_prmtop} \
           -rp ${receptor_prmtop} \
           -lp ${ligand_prmtop} \
            -y ${trajectory}

Submit your script

qsub mmgbsa.sh

The last line of our 2BXF.mmgbsa.results.dat file shows that our delta G binding is -24.4431 +/- 5.4851.