Difference between revisions of "2024 AMBER tutorial 1 with PDBID 2ITO"

From Rizzo_Lab
Jump to: navigation, search
Line 6: Line 6:
  
 
Update the directory structure to the following:
 
Update the directory structure to the following:
 +
  001_structure
 +
  002_forceFieldParameters
 +
  003_tleap
 +
  004_equilibration
 +
  005_productionRun
 +
  006a_RMSD
 +
  006b_hbonds
 +
  006c_mmgbsa
 
    
 
    
  
Line 92: Line 100:
  
 
We need to equilibrate the system before running the MD simulations. Note, for files 1-7, set the restraint mask to the entire system (1-325). For files 8-9, set the restraint mask to cover only the receptor (1-324).
 
We need to equilibrate the system before running the MD simulations. Note, for files 1-7, set the restraint mask to the entire system (1-325). For files 8-9, set the restraint mask to cover only the receptor (1-324).
 +
 +
Create:
  
 
   vi 01.min.mdin
 
   vi 01.min.mdin
 +
 +
With the following:
 +
 +
  Minmize all the hydrogens
 +
  &cntrl
 +
  imin=1,          ! Minimize the initial structure
 +
  maxcyc=5000,    ! Maximum number of cycles for minimization
 +
  ntb=1,            ! Constant volume
 +
  ntp=0,            ! No pressure scaling
 +
  ntf=1,            ! Complete force evaluation
 +
  ntwx= 1000,      ! Write to trajectory file every ntwx steps
 +
  ntpr= 1000,      ! Print to mdout every ntpr steps
 +
  ntwr= 1000,      ! Write a restart file every ntwr steps
 +
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
 +
  ntr=1,            ! Turn on restraints
 +
  restraintmask=":1-325 & !@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
 +
  /
 +
 +
Create:
 +
 
   vi 02.equil.mdin
 
   vi 02.equil.mdin
 +
 +
With the following:
 +
 +
  MD simulation
 +
  &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-325 & !@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
 +
  /
 +
 +
 +
Create:
 +
 
   vi 03.min.mdin
 
   vi 03.min.mdin
 +
 +
With the following:
 +
 +
  Minmize all the hydrogens
 +
  &cntrl
 +
  imin=1,          ! Minimize the initial structure
 +
  maxcyc=1000,    ! Maximum number of cycles for minimization
 +
  ntb=1,            ! Constant volume
 +
  ntp=0,            ! No pressure scaling
 +
  ntf=1,            ! Complete force evaluation
 +
  ntwx= 1000,      ! Write to trajectory file every ntwx steps
 +
  ntpr= 1000,      ! Print to mdout every ntpr steps
 +
  ntwr= 1000,      ! Write a restart file every ntwr steps
 +
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
 +
  ntr=1,            ! Turn on restraints
 +
  restraintmask=":1-325 & !@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
 +
  /
 +
 +
 +
Create:
 +
 
   vi 04.min.mdin
 
   vi 04.min.mdin
 +
 +
With the following:
 +
 +
  Minmize all the hydrogens
 +
  &cntrl
 +
  imin=1,          ! Minimize the initial structure
 +
  maxcyc=1000,    ! Maximum number of cycles for minimization
 +
  ntb=1,            ! Constant volume
 +
  ntp=0,            ! No pressure scaling
 +
  ntf=1,            ! Complete force evaluation
 +
  ntwx= 1000,      ! Write to trajectory file every ntwx steps
 +
  ntpr= 1000,      ! Print to mdout every ntpr steps
 +
  ntwr= 1000,      ! Write a restart file every ntwr steps
 +
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
 +
  ntr=1,            ! Turn on restraints
 +
  restraintmask=":1-325 & !@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
 +
  /
 +
 +
 +
Create:
 +
 
   vi 05.min.mdin
 
   vi 05.min.mdin
 +
 +
With the following:
 +
 +
  Minmize all the hydrogens
 +
  &cntrl
 +
  imin=1,          ! Minimize the initial structure
 +
  maxcyc=1000,    ! Maximum number of cycles for minimization
 +
  ntb=1,            ! Constant volume
 +
  ntp=0,            ! No pressure scaling
 +
  ntf=1,            ! Complete force evaluation
 +
  ntwx= 1000,      ! Write to trajectory file every ntwx steps
 +
  ntpr= 1000,      ! Print to mdout every ntpr steps
 +
  ntwr= 1000,      ! Write a restart file every ntwr steps
 +
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
 +
  ntr=1,            ! Turn on restraints
 +
  restraintmask=":1-325 & !@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
 +
  /
 +
 +
 +
Create:
 +
 
   vi 06.equil.mdin
 
   vi 06.equil.mdin
 +
 +
With the following:
 +
 +
  MD simulation
 +
  &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-325 & !@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
 +
  /
 +
 +
 +
Create:
 +
 
   vi 07.equil.mdin
 
   vi 07.equil.mdin
 +
 +
With the following:
 +
 +
  MD simulation
 +
  &cntrl
 +
  imin=0,          ! Perform MD
 +
  nstlim=50000      ! Number of MD steps
 +
  ntx=5,            ! Positions and velocities read formatted
 +
  irest=1,          ! Restart calculation
 +
  ntc=1,            ! No SHAKE on for bonds with hydrogen
 +
  dt=0.001,        ! Timestep (ps)
 +
  ntb=2,            ! Constant Pressure
 +
  ntp=1,            ! Isotropic pressure scaling
 +
  barostat=1        ! Berendsen
 +
  taup=0.5          ! Pressure relaxtion time (ps)
 +
  ntf=1,            ! Complete force evaluation
 +
  ntt=3,            ! Langevin thermostat
 +
  gamma_ln=2.0      ! Collision Frequency for thermostat
 +
  ig=-1,            ! Random seed for thermostat
 +
  temp0=298.15      ! Simulation temperature (K)
 +
  ntwx= 1000,      ! Write to trajectory file every ntwx steps
 +
  ntpr= 1000,      ! Print to mdout every ntpr steps
 +
  ntwr= 1000,      ! Write a restart file every ntwr steps
 +
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
 +
  ntr=1,            ! Turn on restraints
 +
  restraintmask=":1-325 & !@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:
 +
 
   vi 08.equil.mdin
 
   vi 08.equil.mdin
 +
 +
With the following:
 +
 +
  MD simulations
 +
  &cntrl
 +
  imin=0,          ! Perform MD
 +
  nstlim=50000      ! Number of MD steps
 +
  ntx=5,            ! Positions and velocities read formatted
 +
  irest=1,          ! Restart calculation
 +
  ntc=1,            ! No SHAKE on for bonds with hydrogen
 +
  dt=0.001,        ! Timestep (ps)
 +
  ntb=2,            ! Constant Pressure
 +
  ntp=1,            ! Isotropic pressure scaling
 +
  barostat=1        ! Berendsen
 +
  taup=0.5          ! Pressure relaxtion time (ps)
 +
  ntf=1,            ! Complete force evaluation
 +
  ntt=3,            ! Langevin thermostat
 +
  gamma_ln=2.0      ! Collision Frequency for thermostat
 +
  ig=-1,            ! Random seed for thermostat
 +
  temp0=298.15      ! Simulation temperature (K)
 +
  ntwx= 1000,      ! Write to trajectory file every ntwx steps
 +
  ntpr= 1000,      ! Print to mdout every ntpr steps
 +
  ntwr= 1000,      ! Write a restart file every ntwr steps
 +
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
 +
  ntr=1,            ! Turn on restraints
 +
  restraintmask=":1-324@CA.C.N", ! atoms to be restrained, only the backbone
 +
  restraint_wt=0.1, ! force constant for restraint
 +
  ntxo=1,          ! Write coordinate file in ASCII format
 +
  ioutfm=0,        ! Write trajectory file in ASCII format
 +
  iwrap=1,          ! iwrap is turned on
 +
  /
 +
 +
 +
Create:
 +
 
   vi 09.equil.mdin
 
   vi 09.equil.mdin
 +
 +
With the following:
 +
 +
  MD simulations
 +
  &cntrl
 +
  imin=0,          ! Perform MD
 +
  nstlim=50000      ! Number of MD steps
 +
  ntx=5,            ! Positions and velocities read formatted
 +
  irest=1,          ! Restart calculation
 +
  ntc=1,            ! No SHAKE on for bonds with hydrogen
 +
  dt=0.001,        ! Timestep (ps)
 +
  ntb=2,            ! Constant Pressure
 +
  ntp=1,            ! Isotropic pressure scaling
 +
  barostat=1        ! Berendsen
 +
  taup=0.5          ! Pressure relaxtion time (ps)
 +
  ntf=1,            ! Complete force evaluation
 +
  ntt=3,            ! Langevin thermostat
 +
  gamma_ln=2.0      ! Collision Frequency for thermostat
 +
  ig=-1,            ! Random seed for thermostat
 +
  temp0=298.15      ! Simulation temperature (K)
 +
  ntwx= 1000,      ! Write to trajectory file every ntwx steps
 +
  ntpr= 1000,      ! Print to mdout every ntpr steps
 +
  ntwr= 1000,      ! Write a restart file every ntwr steps
 +
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
 +
  ntr=1,            ! Turn on restraints
 +
  restraintmask=":1-324@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 the following slurm file with the following and run:
 +
 +
  2ITO.equil.slurm
 +
 +
  #!/bin/bash
 +
  #
 +
  #SBATCH --job-name=2ITO_equilibration
 +
  #SBATCH --output=equilibration_output.txt
 +
  #SBATCH --ntasks-per-node=24
 +
  #SBATCH --nodes=1
 +
  #SBATCH --time=48:00:00
 +
  #SBATCH -p long-24core
 +
 +
  module load amber/16
 +
 +
  do_parallel="mpirun -np 80 pmemd.MPI"
 +
 +
  prmtop="../003_tleap/2ITO.wet.complex.prmtop"
 +
  coords="../003_tleap/2ITO.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
 +
  
 
=Production=
 
=Production=
Line 110: Line 412:
  
 
and add the following:
 
and add the following:
 +
 +
  MD simulations
 +
  &cntrl
 +
  imin=0,          ! Perform MD
 +
  nstlim=5000000,  ! Number of MD steps
 +
  ntx=5,            ! Positions and velocities read formatted
 +
  irest=1,          ! Restart calculation
 +
  ntc=2,            ! SHAKE on for bonds with hydrogen
 +
  dt=0.002,        ! Timestep (ps)
 +
  ntb=2,            ! Constant Pressure
 +
  ntp=1,            ! Isotropic pressure scaling
 +
  barostat=1        ! Berendsen
 +
  taup=0.5          ! Pressure relaxtion time (ps)
 +
  ntf=2,            ! No force evaluation for bonds with hydrogen
 +
  ntt=3,            ! Langevin thermostat
 +
  gamma_ln=2.0      ! Collision Frequency for thermostat
 +
  ig=-1,            ! Random seed for thermostat
 +
  temp0=298.15      ! Simulation temperature (K)
 +
  ntwx= 2500,      ! Write to trajectory file every ntwx steps
 +
  ntpr= 2500,      ! Print to mdout every ntpr steps
 +
  ntwr= 5000000,    ! Write a restart file every ntwr steps
 +
  cut=8.0,          ! Nonbonded cutoff in Angstroms
 +
  ntr=1,            ! Turn on restraints
 +
  restraintmask=":1-324@CA,C,N", ! atoms to be restrained
 +
  restraint_wt=0.1, ! force constant for restraint
 +
  ntxo=1,          ! Write coordinate file in ASCII format
 +
  ioutfm=0,        ! Write trajectory file in ASCII format
 +
  iwrap=1,          ! iwrap is turned on
 +
  $end
  
 
Create the slurm file:
 
Create the slurm file:
Line 116: Line 447:
  
 
and add the following:
 
and add the following:
 +
 +
  #!/bin/bash
 +
  #
 +
  #SBATCH --job-name=2ITO_production
 +
  #SBATCH --output=2ITO_production.txt
 +
  #SBATCH --ntasks-per-node=28
 +
  #SBATCH --nodes=4
 +
  #SBATCH --time=48:00:00
 +
  #SBATCH -p long-28core
 +
 +
  module load amber/16
 +
 +
  do_parallel="mpirun pmemd.MPI"
 +
 +
  prmtop="../003_tleap/2ITO.wet.complex.prmtop"
 +
  coords="../004_equil/09.equil"
 +
 +
  MDINPUTS=(10.prod)
 +
 +
  for input in ${MDINPUTS[@]}; do
 +
  $do_parallel -O -i ${input}.mdin -o ${input}.mdout -p $prmtop -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf ${input}.info -r ${input}.rst7
 +
  coords=$input
 +
  done
  
 
Run the slurm script
 
Run the slurm script
Line 127: Line 481:
 
First we want to strip all the water:
 
First we want to strip all the water:
  
−
vi cpp_strip_water.in
+
  vi cpp_strip_water.in
  
 
with the following:
 
with the following:
  
 +
  #!/usr/bin/sh
 +
  parm ../003_tleap/2ITO.wet.complex.prmtop
 +
  #read in trajectory
 +
  trajin ../005_production/10.prod.trj
 +
  #read in reference
 +
  reference ../003_tleap/2ITO.wet.complex.rst7
 +
  #compute RMSD + align CA to crystal structure
 +
  rmsd rms1 reference :1-325@CA
 +
  #strip solvent
 +
  strip :WAT:Na+:Cl-
 +
  #create gas-phase trajectory
 +
  trajout 2ITO_stripfit_restrained_gas.trj nobox
 +
  cpp_strip_water.in (END)
 +
 +
Run it with the following command:
 +
 +
  cpptraj -i cpp_strip_water.in
  
 
To look at the RMSD of the ligand, we only want to look at the residue that includes the ligand:
 
To look at the RMSD of the ligand, we only want to look at the residue that includes the ligand:
  
−
vi cpp_rmsd_lig.in
+
  vi cpp_rmsd_lig.in
  
 
with the following:
 
with the following:
  
−
 
+
  #!/usr/bin/sh
−
Run it with the following command:
+
  #trajin the trajectory
−
 
+
  trajin 2ITO_stripfit_restrained_gas.trj
−
   cpptraj -i cpp_strip_water.in
+
  #read in reference
−
 
+
  reference ../003_tleap/2ITO.gas.complex.rst7
 +
  #compute RMSD (do not fit internal geometris first, included rigid body motions, convert frames to ns (framenum*.005)
 +
  rmsd rms1 ":325&!(@H=)" nofit mass out 2ITO_lig_restrained_rmsd_nofit.dat time .005
 +
  #histogram the nofit rmsd
 +
   histogram rms1,*,*,.1,* norm out 2ITO_lig_restrained_rmsd_nofit_histogram.dat
  
 
Run it with the following command:
 
Run it with the following command:
Line 153: Line 528:
  
 
with the following:
 
with the following:
 +
  #!/usr/bin/sh
 +
  #trajin the trajectory
 +
  trajin 2ITO_stripfit_restrained_gas.trj
 +
  #read in reference
 +
  reference ../003_tleap/2ITO.gas.complex.rst7
 +
  #compute RMSD (do not fit internal geometries first, included rigid body motions, convert frames to ns (framenum*.005)
 +
  rmsd rms1 ":1-324&!(@H=)" nofit mass out 2ITO_receptor_restrained_rmsd_nofit.dat time .005
 +
  #histogram the nofit rmsd
 +
  histogram rms1,*,*,.1,* norm out 2ITO_receptor_restrained_rmsd_nofit_histogram.dat
  
  
Line 168: Line 552:
  
 
with the following:
 
with the following:
 +
 +
  !/usr/bin/sh
 +
  #read in trajectory
 +
  trajin ../004_production/10_prod.trj
 +
  #wrap everything in one periodic cell - could cause problems, may comment out #autoimage if problems later
 +
  autoimage
 +
  #compute intra + water-mediated H-bonds
 +
  hbond hb1 :1-325 out 2ITO_hbond.out avgout 2ITO_hbond_avg.dat solventdonor :WAT solventacceptor :WAT@-O
 +
  nointramol bridgeout 2ITO_bridge-water.dat dist 3.0 angle 140
 +
  
 
Run it with the following command:
 
Run it with the following command:
Line 181: Line 575:
 
with the following:
 
with the following:
  
−
and run.
+
  mmgbsa 2ITO 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 the slurm file:
 +
 
 +
  vi 2ITO_mmgbsa.slurm
 +
 
 +
with the following:
 +
 
 +
  #bin/bash
 +
  #SBATCH --time=2-00:00:00
 +
  #SBATCH --nodes=2
 +
  #SBATCH --ntasks=28
 +
  #SBATCH --job-name=2ITO_MMGBSA
 +
  #SBATCH --output=2ITO_MMGBSA.log
 +
  #SBATCH -p extended-28core
 +
  #Define topology files
 +
  solv_parm="../003_tleap/2ITO.wet.complex.prmtop"
 +
  complex_parm="../003_tleap/2ITO.complex.parm7"
 +
  receptor_parm="../003_tleap/2ITO.gas.receptor.parm7"
 +
  lig_parm="../003_tleap/2ITO.gas.ligand.parm7"
 +
  trajectory="../005_production/10_prod.trj"
 +
  MMPBSA.py -O -i mmgbsa.in \
 +
      -o 2ITO_mmgbsa_results.dat \
 +
      -eo 2ITO_mmgbsa_perframe.dat \
 +
      -sp ${solv_parm} \
 +
      -cp ${complex_parm} \
 +
      -rp ${receptor_parm} \
 +
      -lp ${lig_parm} \
 +
      -y ${trajectory}
 +
 
 +
 
 +
and run the slurm file.
 +
 
 +
sbatch 2ITO_mmgbsa.slurm

Revision as of 02:21, 6 May 2024

Introduction

AMBER is a molecular dynamics program that simulate the behavior of a protein/ligand complex of interest. AMBER considers all atoms and models their interactions over time, continuously updating positions and forces in small time increments. For this tutorial, we will focus on 2ITO.

Getting Started

Update the directory structure to the following: 

 001_structure
 002_forceFieldParameters
 003_tleap
 004_equilibration
 005_productionRun
 006a_RMSD
 006b_hbonds
 006c_mmgbsa

Structure

Similarly to the DOCK tutorial, we need to obtain the structure of the protein and ligand. For the protein, we want to fill in all loops even if it is not near the binding pocket, since AMBER considers all atoms in its calculations. Follow the DOCK tutorial for 2ITO to do this.

For the protein, do not add H or add charge. Save as 2ITO_protein_only_for_AMBER.pdb and 2ITO_protein_only_for_AMBER.mol2. For the ligand, add H and add charge. Save as 2ITO_ligand_only_for_AMBER.pdb and 2ITO_ligand_only_for_AMBER.mol2.

Force Field Parameters

Force field parameters are required to calculation the interactions between all atoms. This needs to be done for the ligand. Create the following slurm script as follows: 

 vi 2ITO_ffp.slurm

and have the following within the file:

Make sure to change LIG -nc 0 based off of your ligand. It is the charge on the ligand. For 2ITO, this charge is +1.

Once the slurm file has run, the following will show up in your directory: 

 2ITO_ligand_antechamber.mol2  
 ANTECHAMBER_AC.AC0     
 ANTECHAMBER_AM1BCC_PRE.AC  
 ANTECHAMBER_BOND_TYPE.AC0  
 sqm.in   
 sqm.pdb
 ANTECHAMBER_AC.AC             
 ANTECHAMBER_AM1BCC.AC  
 ANTECHAMBER_BOND_TYPE.AC   
 ATOMTYPE.INF               
 sqm.out

2ITO_ligand_antechamber.mol2 is the file with the parameters we generated.

In the command line, load amber 16: 

 module load amber/16

Run the following command to modify the parameters slightly: 

parmchk2 -i 2ITO_ligand_antechamber.mol2 -f mol2 -o 2ITO_ligand.am1bcc.frcmod


You will see 2ITO_ligand.am1bcc.frcmod file in your directory.

TLeap

TLeap is a program with AMBER that will prepare the system for MD simulations. It will generate the gas-phase and solvated systems.

Create the file: 

 vi 2ITO_tleap.in

and add the following to the file:

Create the following slurm file: 

 vi 2ITO_tleap.slurm

and add the following to the file:


Run the slurm and you should get the following new files in your directory: 

 2ITO.complex.parm7    
 2ITO.gas.complex.rst7  
 2ITO.gas.ligand.rst7     
 2ITO.gas.receptor.rst7      
 leap.log
 2ITO.gas.complex.pdb  
 2ITO.gas.ligand.parm7  
 2ITO.gas.receptor.parm7       
 2ITO.wet.complex.pdb  
 2ITO.wet.complex.prmtop
 2ITO.wet.complex.rst7    
 tleap_output.txt

You can look at the files you created by moving the files to your local computer. Move 2ITO.wet.complex.rst7 and 2ITO.wet.complex.prmtop. Open chimera and follow the following commands: 

 Tools → MD/Ensemble Analysis → MD Movie

Click browse and put in the prmtop file. Add the rst7 trajectory file. You should be able to view your system.

Equilibration

Look at the 2ITO.gas.complex.pdb file that was just generated. Scroll all the way down to see the residues in the complex. In our example, the protein is residues 1-324 and the ligand is 325.

We need to equilibrate the system before running the MD simulations. Note, for files 1-7, set the restraint mask to the entire system (1-325). For files 8-9, set the restraint mask to cover only the receptor (1-324).

Create: 

 vi 01.min.mdin

With the following: 

 Minmize all the hydrogens
 &cntrl
 imin=1,           ! Minimize the initial structure
 maxcyc=5000,    ! Maximum number of cycles for minimization
 ntb=1,            ! Constant volume
 ntp=0,            ! No pressure scaling
 ntf=1,            ! Complete force evaluation
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-325 & !@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
 /

Create: 

 vi 02.equil.mdin

With the following: 

 MD simulation
 &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-325 & !@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
 /


Create: 

 vi 03.min.mdin

With the following: 

 Minmize all the hydrogens
 &cntrl
 imin=1,           ! Minimize the initial structure
 maxcyc=1000,    ! Maximum number of cycles for minimization
 ntb=1,            ! Constant volume
 ntp=0,            ! No pressure scaling
 ntf=1,            ! Complete force evaluation
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-325 & !@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
 /


Create: 

 vi 04.min.mdin

With the following: 

 Minmize all the hydrogens
 &cntrl
 imin=1,           ! Minimize the initial structure
 maxcyc=1000,    ! Maximum number of cycles for minimization
 ntb=1,            ! Constant volume
 ntp=0,            ! No pressure scaling
 ntf=1,            ! Complete force evaluation
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-325 & !@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
 /


Create: 

 vi 05.min.mdin

With the following: 

 Minmize all the hydrogens
 &cntrl
 imin=1,           ! Minimize the initial structure
 maxcyc=1000,    ! Maximum number of cycles for minimization
 ntb=1,            ! Constant volume
 ntp=0,            ! No pressure scaling
 ntf=1,            ! Complete force evaluation
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-325 & !@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
 /


Create: 

 vi 06.equil.mdin

With the following: 

 MD simulation
 &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-325 & !@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
 /


Create: 

 vi 07.equil.mdin

With the following: 

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

 vi 08.equil.mdin

With the following: 

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


Create: 

 vi 09.equil.mdin

With the following: 

 MD simulations
 &cntrl
 imin=0,           ! Perform MD
 nstlim=50000      ! Number of MD steps
 ntx=5,            ! Positions and velocities read formatted
 irest=1,          ! Restart calculation
 ntc=1,            ! No SHAKE on for bonds with hydrogen
 dt=0.001,         ! Timestep (ps)
 ntb=2,            ! Constant Pressure
 ntp=1,            ! Isotropic pressure scaling
 barostat=1        ! Berendsen
 taup=0.5          ! Pressure relaxtion time (ps)
 ntf=1,            ! Complete force evaluation
 ntt=3,            ! Langevin thermostat
 gamma_ln=2.0      ! Collision Frequency for thermostat
 ig=-1,            ! Random seed for thermostat
 temp0=298.15      ! Simulation temperature (K)
 ntwx= 1000,       ! Write to trajectory file every ntwx steps
 ntpr= 1000,       ! Print to mdout every ntpr steps
 ntwr= 1000,       ! Write a restart file every ntwr steps
 cut=  8.0,        ! Nonbonded cutoff in Angstroms
 ntr=1,            ! Turn on restraints
 restraintmask=":1-324@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 the following slurm file with the following and run: 

 2ITO.equil.slurm

 #!/bin/bash
 #
 #SBATCH --job-name=2ITO_equilibration
 #SBATCH --output=equilibration_output.txt
 #SBATCH --ntasks-per-node=24
 #SBATCH --nodes=1
 #SBATCH --time=48:00:00
 #SBATCH -p long-24core

 module load amber/16

 do_parallel="mpirun -np 80 pmemd.MPI"

 prmtop="../003_tleap/2ITO.wet.complex.prmtop"
 coords="../003_tleap/2ITO.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


Production

The system is all equilibrated and now we can run production run. 

 vi 10.prod.mdin

and add the following: 

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

Create the slurm file: 

 production.slurm

and add the following: 

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

 module load amber/16

 do_parallel="mpirun pmemd.MPI"

 prmtop="../003_tleap/2ITO.wet.complex.prmtop"
 coords="../004_equil/09.equil"

 MDINPUTS=(10.prod)

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

Run the slurm script

Analysis

Now we can analyze the simulation by looking at the RMSD, H-Bonds, and MM/GBSA results.

RMSD

First we want to strip all the water: 

 vi cpp_strip_water.in

with the following: 

 #!/usr/bin/sh
 parm ../003_tleap/2ITO.wet.complex.prmtop
 #read in trajectory
 trajin ../005_production/10.prod.trj
 #read in reference
 reference ../003_tleap/2ITO.wet.complex.rst7
 #compute RMSD + align CA to crystal structure
 rmsd rms1 reference :1-325@CA
 #strip solvent
 strip :WAT:Na+:Cl-
 #create gas-phase trajectory
 trajout 2ITO_stripfit_restrained_gas.trj nobox
 cpp_strip_water.in (END)

Run it with the following command: 

 cpptraj -i cpp_strip_water.in

To look at the RMSD of the ligand, we only want to look at the residue that includes the ligand: 

 vi cpp_rmsd_lig.in

with the following: 

 #!/usr/bin/sh
 #trajin the trajectory
 trajin 2ITO_stripfit_restrained_gas.trj
 #read in reference
 reference ../003_tleap/2ITO.gas.complex.rst7
 #compute RMSD (do not fit internal geometris first, included rigid body motions, convert frames to ns (framenum*.005)
 rmsd rms1 ":325&!(@H=)" nofit mass out 2ITO_lig_restrained_rmsd_nofit.dat time .005
 #histogram the nofit rmsd
 histogram rms1,*,*,.1,* norm out 2ITO_lig_restrained_rmsd_nofit_histogram.dat

Run it with the following command: 

 cpptraj -p ../003_tleap/2ITO.complex.parm7 -i cpp_rmsd_lig.in


Now we want to look at the RMSD of the receptor. Include the residues that include only the receptor:

vi cpp_rmsd_receptor.in

with the following: 

 #!/usr/bin/sh
 #trajin the trajectory
 trajin 2ITO_stripfit_restrained_gas.trj
 #read in reference
 reference ../003_tleap/2ITO.gas.complex.rst7
 #compute RMSD (do not fit internal geometries first, included rigid body motions, convert frames to ns (framenum*.005)
 rmsd rms1 ":1-324&!(@H=)" nofit mass out 2ITO_receptor_restrained_rmsd_nofit.dat time .005
 #histogram the nofit rmsd
 histogram rms1,*,*,.1,* norm out 2ITO_receptor_restrained_rmsd_nofit_histogram.dat


Run it with the following command: 

 cpptraj -p ../003_tleap/2ITO.complex.parm7 -i cpp_rmsd_receptor.in

If you want to look at the histogram, you can graph it using python.


H-Bonds

To look at the hydrogen bonds that occur in the system during the simulation, make sure to include the receptor and ligand in the mask. Create the program: 

 vi cpp_hbonds.in

with the following: 

 !/usr/bin/sh
 #read in trajectory
 trajin ../004_production/10_prod.trj
 #wrap everything in one periodic cell - could cause problems, may comment out #autoimage if problems later
 autoimage
 #compute intra + water-mediated H-bonds
 hbond hb1 :1-325 out 2ITO_hbond.out avgout 2ITO_hbond_avg.dat solventdonor :WAT solventacceptor :WAT@-O
 nointramol bridgeout 2ITO_bridge-water.dat dist 3.0 angle 140


Run it with the following command: 

 cpptraj -p ../003_tleap/2ITO.wet.complex.prmtop -i cpp_hbonds.in


MM/GBSA

Now we will calcualte the free energy of the system. Create the program: 

 vi mmgbsa_2ITO.in

with the following: 

 mmgbsa 2ITO 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 the slurm file: 

 vi 2ITO_mmgbsa.slurm

with the following: 

 #bin/bash
 #SBATCH --time=2-00:00:00
 #SBATCH --nodes=2
 #SBATCH --ntasks=28
 #SBATCH --job-name=2ITO_MMGBSA
 #SBATCH --output=2ITO_MMGBSA.log
 #SBATCH -p extended-28core
 #Define topology files
 solv_parm="../003_tleap/2ITO.wet.complex.prmtop"
 complex_parm="../003_tleap/2ITO.complex.parm7"
 receptor_parm="../003_tleap/2ITO.gas.receptor.parm7"
 lig_parm="../003_tleap/2ITO.gas.ligand.parm7"
 trajectory="../005_production/10_prod.trj"
 MMPBSA.py -O -i mmgbsa.in \
      -o 2ITO_mmgbsa_results.dat \ 
      -eo 2ITO_mmgbsa_perframe.dat \
      -sp ${solv_parm} \
      -cp ${complex_parm} \
      -rp ${receptor_parm} \
      -lp ${lig_parm} \
      -y ${trajectory}


and run the slurm file. 

sbatch 2ITO_mmgbsa.slurm
Retrieved from "https://ringo.ams.stonybrook.edu/index.php?title=2024_AMBER_tutorial_1_with_PDBID_2ITO&oldid=17383"