Difference between revisions of "2022 AMBER tutorial 1 with PDBID 6ME2"

From Rizzo_Lab
Jump to: navigation, search
(→‎Receptor File Generation)
(→‎MD Analysis)
 
(7 intermediate revisions by the same user not shown)
Line 33: Line 33:
 
   Actions -> Atoms/Bonds -> delete
 
   Actions -> Atoms/Bonds -> delete
  
−
which will isolate the receptor. Go to  
+
which will isolate the receptor. Go to  
  
 
   File -> Save Mol2...  
 
   File -> Save Mol2...  
Line 41: Line 41:
 
=== '''Ligand File Generation''' ===
 
=== '''Ligand File Generation''' ===
  
−
To generate the ligand file, open the '''6me2.pdb''' file again, and go to '''Select -> Structure -> protein''' then '''Actions -> Atoms/Bonds -> delete''' to delete the receptor. Then, go to '''Select -> Residue -> HOH''', '''Select -> Residue -> OLA''', and '''Select -> Residue -> PEG''' to delete everything except for the '''JEV''' residue, which is the ligand of interest.
+
To generate the ligand file, open the '''6me2.pdb''' file again, and go to  
  
−
Note: there are problems with the file containing the ligand as-is. In order to proceed properly, it is necessary to make some changed to this file. While the file is still open, go to '''Tools > Structure Editing > Rotamers'' and change the non-standard YCM residue to the most probable rotamers of CYS.
+
  Select -> Structure -> protein
  
−
After this is completed, go to '''Tools > Structure editing > Add H''' to add hydrogens to the ligand. Then, go to '''Tools > Structure editing > Add Charge > (have Amber ff14SB and AM1-BCC selected) -> Ok''' and then go to '''File -> Save Mol2...''' and save as '''6me2_lig_dockprep.mol2'''. Close Chimera.
+
then
 +
 
 +
  Actions -> Atoms/Bonds -> delete
 +
 
 +
to delete the receptor. Then, go to
 +
 
 +
  Select -> Residue -> HOH -> Actions -> Atoms/Bonds -> delete
 +
 
 +
  Select -> Residue -> OLA -> Actions -> Atoms/Bonds -> delete
 +
 
 +
  Select -> Residue -> PEG -> Actions -> Atoms/Bonds -> delete
 +
 
 +
to delete everything except for the '''JEV''' residue, which is the ligand of interest here.
 +
 
 +
Note: there are problems with the file containing the ligand as-is. In order to proceed properly, it is necessary to make some changes to this file. While the file is still open, go to  
 +
 
 +
  Tools > Structure Editing > Rotamers
 +
 
 +
and change the non-standard '''YCM''' residue to the most probable rotamers of '''CYS'''.
 +
 
 +
After this is completed, go to
 +
 
 +
  Tools > Structure editing > Add H
 +
 
 +
to add hydrogens to the ligand. Then, go to  
 +
 
 +
  Tools > Structure editing > Add Charge > (have Amber ff14SB and AM1-BCC selected) -> Ok
 +
 
 +
and then go to  
 +
 
 +
  File -> Save Mol2...  
 +
 
 +
and save as '''6me2_lig_dockprep.mol2'''. Close Chimera.
  
 
== '''Generating Simulation Parameters''' ==
 
== '''Generating Simulation Parameters''' ==
Line 63: Line 95:
  
 
== '''Using TLeap''' ==
 
== '''Using TLeap''' ==
 +
Now, switch to the following directory:
 +
 +
  cd 003_leap
 +
 +
Where files will be generated to simulate the ligand and receptor, such that further calculations involving both structures--such as binding affinities, free energies, and trajectories. Two types of files will be generated that can be used to set up the system: parameter ('''parm7''') and restart ('''rst7''') files. Start the process by creating the following input file:
 +
 +
  vi leap.in
 +
 +
And put the following into the input 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 ../001_structure/6me2_fresh.pdb
 +
  ###load ligand frcmod/mol2
 +
  loadamberparams ../002_parameters/6me2_ligand.am1bcc.frcmod
 +
  lig=loadmol2 ../002_parameters/1HW9_ligand_antechamber.mol2
 +
  ###create gase-phase complex
 +
  gascomplex= combine {rec lig}
 +
  ###write gas-phase pdb
 +
  savepdb gascomplex 6me2.gas.complex.pdb
 +
  ###write gase-phase toplogy and coord files for MMGBSA calc
 +
  saveamberparm gascomplex 6me2.complex.parm7 6me2.gas.complex.rst7
 +
  saveamberparm rec 6me2.gas.receptor.parm7 6me2.gas.receptor.rst7
 +
  saveamberparm lig 6me2.gas.ligand.parm7 6me2.gas.ligand.rst7
 +
  ###create solvated complex (albeit redundant)
 +
  solvcomplex= combine {rec lig}
 +
  ###solvate the system
 +
  solvateoct solvcomplex TIP3PBOX 12.0
 +
  ###Neutralize system
 +
  addions solvcomplex Cl- 0
 +
  addions solvcomplex Na+ 0
 +
  #write solvated pdb file
 +
  savepdb solvcomplex 6me2.wet.complex.pdb
 +
  ###check the system
 +
  charge solvcomplex
 +
  check solvcomplex
 +
  ###write solvated toplogy and coordinate file
 +
  saveamberparm solvcomplex 6me2.wet.complex.parm7 6me2.wet.complex.rst7
 +
  quit
 +
 +
After this runs, copy the '''parm7''' (parameter) and '''rst7''' (coordinate) files to the local machine, where they can be opened in Chimera. To open in Chimera, go to
 +
 +
  Tools -> MD/Ensemble Analysis -> MD Movie --> Choose appropriate parameter and coordinate files to load
  
 
== '''Equilibrating the system''' ==
 
== '''Equilibrating the system''' ==
 +
Now, move to the following directory:
 +
 +
  cd 004_equil
 +
 +
In this process, energy minimization and equilibration will be performed on the system. There are 9 input files that need to be generated for this file, as follows:
 +
 +
  vi 01.min.mdin
 +
 +
With the following input:
 +
 +
  Minimize 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="!@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
 +
  /
 +
 +
Next,
 +
 +
  vi 02.equil.mdin
 +
 +
with input
 +
 +
  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=":!@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 
 +
  /
 +
 +
Next,
 +
 +
  vi 03.min.mdin
 +
 +
with input
 +
 +
  Minimize 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="!@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
 +
  /
 +
 +
Next,
 +
 +
  vi 04.min.mdin
 +
 +
with input
 +
 +
  Minimize 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="!@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
 +
  /
 +
 +
Next,
 +
 +
  vi 05.min.mdin
 +
 +
with input
 +
 +
  Minimize 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="!@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
 +
  /
 +
 +
Next,
 +
 +
  vi 06.equil.mdin
 +
 +
with input
 +
 +
  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="!@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
 +
  /
 +
 +
Next,
 +
 +
  vi 07.equil.mdin
 +
 +
with input
 +
 +
  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="!@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
 +
  /
 +
 +
Next,
 +
 +
  vi 08.equil.mdin *See below for specific comments on this input
 +
 +
with input
 +
 +
  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-433@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
 +
  /
 +
 +
And finally,
 +
 +
  vi 09.equil.mdin *See below for comments on this input
 +
 +
with input
 +
 +
  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-433@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
 +
  /
 +
 +
Note that the restraint mask in files 8 and 9 have to be specified to the relevant system being simulated. These are residue numbers which can be found in the '''wet.complex.pdb''' file generated in the '''003_leap''' directory.
 +
 +
To submit all 9 input files as a job on Seawulf, create the following script:
 +
 +
  vi mdequilibration.sh
 +
 +
with input
 +
 +
  #!/bin/sh
 +
  #SBATCH --job-name=6me2_equilibration
 +
  #SBATCH --ntasks-per-node=40
 +
  #SBATCH --nodes=2
 +
  #SBATCH --time=8:00:00
 +
  #SBATCH -p long-40core 
 +
  echo "started Equilibration on 'date' "
 +
  parm7="../003_leap/6me2.wet.complex.parm7"
 +
  coords="../003_leap/6me2.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
 +
  pmemd -O -i ${input}.mdin -o ${input}.mdout -p $parm7 -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf ${input}.info -r
 +
  ${input}.rst7
 +
  coords=$input
 +
  done
 +
  echo " Finished equilibration on 'date' "
 +
 +
And submit this to the queue:
 +
 +
  sbatch mdequilibration.sh
  
 
== '''MD Production''' ==
 
== '''MD Production''' ==
 +
Now, enter the following directory:
 +
 +
  cd 005_production
 +
 +
and create the following production input file:
 +
 +
  vi 10.prod.mdin
 +
 +
with inputs
 +
 +
  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-432@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
 +
  /
 +
 +
The restraint mask should be modified on this file as well to only include the receptor, not the ligand (i.e., it should be one fewer than the restraint mask value in the 8th input file from the previous directory).
 +
 +
Create a file that will be submitted to Seawulf:
 +
 +
  vi mdproduction.sh
 +
 +
with the following:
 +
 +
  !/bin/sh
 +
  #SBATCH --job-name=6me2_prod
 +
  #SBATCH --ntasks-per-node=24
 +
  #SBATCH --nodes=2
 +
  #SBATCH --time=7-00:00:00
 +
  #SBATCH -p extended-24core
 +
  echo "started Equilibration on 'date' "
 +
  parm7="../003_leap/6me2.wet.complex.parm7"
 +
  coords="../004_equil/09.equil"
 +
  MDINPUTS=( 10.prod)
 +
  for input in ${MDINPUTS[@]}; do
 +
    pmemd -O -i ${input}.mdin -o ${input}.mdout -p ${parm7} -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf ${input}.info -r
 +
  ${input}.rst7
 +
    coords=$input
 +
  done
 +
  echo "Finished Equilibration on `date` "
 +
 +
And now submit this job to the queue:
 +
 +
  sbatch mdproduction.sh
 +
 +
And note that this job will likely take a very long time, perhaps days.
  
 
== '''MD Analysis''' ==
 
== '''MD Analysis''' ==
 +
 +
Since the trajectories have been generated, they can now be inspected using RMSD calculations, hydrogen bond analysis, and MMGBSA calculations.
  
 
=== '''RMSD Calculation''' ===
 
=== '''RMSD Calculation''' ===
 +
 +
To perform the RMSD calculations, first create an input file which will delete the water and charges from the production trajectory.
 +
 +
  vi cpptraj_strip_wat.in
 +
 +
with the following inputs
 +
 +
  #!/usr/bin/sh
 +
  parm ../003_leap/6me2.wet.complex.parm7
 +
  #read in trajectory
 +
  trajin ../004_production/10.prod.trj
 +
  #read in reference
 +
  reference ../003_leap/6me2.wet.complex.rst7
 +
  #compute RMSD + align CA to crystal structure
 +
  rmsd rms1 reference :1-263@CA
 +
  #strip solvent
 +
  strip :WAT:Na+:Cl-
 +
  #create gas-phase trajectory
 +
  trajout 6me2_stripfit_restrained_gas.trj nobox
 +
 +
The reference mask should contain both the ligand and receptor in its value.
 +
 +
Now run this file using the following command:
 +
 +
  cpptraj -i cpptraj_strip_wat.in
 +
 +
Next, create another input file which will be used to calculate the RMSD of the ligand by typing
 +
 +
  vi cpptraj_rmsd_lig.in
 +
 +
with inputs
 +
 +
  #!/usr/bin/sh
 +
  #trajin the trajectory
 +
  trajin 6me2_stripfit_restrained_gas.trj
 +
  #read in reference
 +
  reference ../003_leap/6me2.gas.complex.rst7
 +
  #compute RMSD (do not fit internal geometris first, included rigid body motions, convert frames to ns (framenum*.005)
 +
  rmsd rms1 ":254&!(@H=)" nofit mass out 6me2_lig_restrained_rmsd_nofit.dat time .005
 +
  #histogram the nofit rmsd
 +
  histogram rms1,*,*,.1,* norm out 6me2_lig_restrained_rmsd_nofit_histogram.dat
 +
 +
And check that the reference mask only factors in the ligand in this file. Run this input file using the following command:
 +
 +
  cpptraj -p ../003_leap/6me2.complex.parm7 -i cpptraj_rmsd_lig.in
 +
 +
Next, the receptor RMSD can be calculated by creating the following input file:
 +
 +
  vi cpptraj_rmsd_rec.in
 +
 +
with inputs
 +
 +
  #!/usr/bin/sh
 +
  #trajin the trajectory
 +
  trajin 6me2_stripfit_restrained_gas.trj
 +
  #read in reference
 +
  reference ../003_leap/6me2.gas.complex.rst7
 +
  #compute RMSD (do not fit internal geometries first, included rigid body motions, convert frames to ns (framenum*.005)
 +
  rmsd rms1 ":1-253&!(@H=)" nofit mass out 6me2_rec_restrained_rmsd_nofit.dat time .005
 +
  #histogram the nofit rmsd
 +
  histogram rms1,*,*,.1,* norm out 6me2_rec_restrained_rmsd_nofit_histogram.dat
 +
 +
And check that the reference mask only takes into account the receptor. Now run this file using the following command:
 +
 +
  cpptraj -p ../003_leap/6me2.gas.complex.parm7 -i cpptraj_rmsd_rec.in
  
 
=== '''Calculating Number of Hydrogen Bonds''' ===
 
=== '''Calculating Number of Hydrogen Bonds''' ===
 +
 +
Next, the number of hydrogen bonds will be analyzed to see if any exist between the receptor and ligand. Create the following input file:
 +
 +
  vi cpptraj_hbond.in
 +
 +
with inputs
 +
 +
  !/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-263 out 6me2_calcipotriol_hbond.out avgout 6me2_calcipotriol_hbond_avg.dat solventdonor :WAT solventacceptor :WAT@-O
 +
  nointramol brid\
 +
  geout 6me2_calcipotriol_bridge-water.dat dist 3.0 angle 140
 +
 +
Be sure that the '''hbobnd hb1''' parameter corresponds to the residue numbers of the relevant ligand and receptor. Now, ru this input file with the following command:
 +
 +
  cpptraj -p ../003_leap/6me2_wetcomplex.parm7 -i cpptraj_hbond.in
  
 
=== '''MMGBSA Calculation''' ===
 
=== '''MMGBSA Calculation''' ===
 +
 +
Next, MMGBSA will be used to calculate free energies of binding in this system. First, create the following input file:
 +
 +
  vi mmgbsa.in
 +
 +
with inputs
 +
 +
  mmgbsa 6me2 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,
 +
  /
 +
 +
To run this on Seawulf, create the following script:
 +
 +
  vi mmgbsa.sh
 +
 +
and input the following
 +
 +
  #bin/bash
 +
  #SBATCH --time=2-00:00:00
 +
  #SBATCH --nodes=2
 +
  #SBATCH --ntasks=28
 +
  #SBATCH --job-name=6me2_MMGBSA
 +
  #SBATCH --output=6me2_MMGBSA.log
 +
  #SBATCH -p extended-28core
 +
  #Define topology files
 +
  solv_parm="../003_leap/6me2.wet.complex.parm7"
 +
  complex_parm="../003_leap/6me2.complex.parm7"
 +
  receptor_parm="../003_leap/6me2.gas.receptor.parm7"
 +
  lig_parm="../003_leap/6me2.gas.ligand.parm7"
 +
  trajectory="10_prod.trj"
 +
  MMPBSA.py -O -i mmgbsa.in \
 +
      -o 6me2_mmgbsa_results.dat \
 +
      -eo 6me2_mmgbsa_perframe.dat \
 +
      -sp ${solv_parm} \
 +
      -cp ${complex_parm} \
 +
      -rp ${receptor_parm} \
 +
      -lp ${lig_parm} \
 +
        -y ${trajectory}
 +
 +
And finally, submit this as a job to the queue by typing
 +
 +
  sbatch mmgbsa.sh

Latest revision as of 19:33, 3 May 2022

Introduction

AMBER is a program designed for computing biomolecular simulations. In this tutorial, AMBER will be used to simulate the dynamics between the ligand and receptor of the 6ME2 PDB file.

Directory Setup

For this tutorial, it is recommended to define the following directories to stay organized: 

 mkdir 001_structure
 mkdir 002_parameters
 mkdir 003_leap
 mkdir 004_equil
 mkdir 005_production

6ME2 Structure Files

If the 6ME2 DOCK tutorials were followed before this AMBER tutorial, the user has some experience generating and prepping receptor and ligand filed from the original 6ME2 PDB file. However, it is recommended that the user NOT use the files generated in the 6ME2 DOCK tutorials and to instead start from scratch here, as this allows the chance to catch mistakes that may have been made but not recognized while following the initial DOCK tutorials.

To begin this tutorial, search the Protein Data Bank (PDB) website using the code 6ME2, or use the Fetch... button on the home screen of UCSF Chimera for this four-letter code. All the following structures should be saved to the following directory: 

 cd 001_structure

Receptor File Generation

To generate a fresh copy of the receptor file, open the 6me2.pdb file containing the ligand-receptor complex. Go to 

 Select -> Structure -> protein

which will select the receptor, and then go to 

 Select -> Invert (all models)

to select everything other than the receptor in the file. Then, go to 

 Actions -> Atoms/Bonds -> delete

which will isolate the receptor. Go to 

 File -> Save Mol2...

and save as 6me2_rec.mol2. Close Chimera.

Ligand File Generation

To generate the ligand file, open the 6me2.pdb file again, and go to 

 Select -> Structure -> protein

then 

 Actions -> Atoms/Bonds -> delete

to delete the receptor. Then, go to 

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

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

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

to delete everything except for the JEV residue, which is the ligand of interest here.

Note: there are problems with the file containing the ligand as-is. In order to proceed properly, it is necessary to make some changes to this file. While the file is still open, go to 

 Tools > Structure Editing > Rotamers

and change the non-standard YCM residue to the most probable rotamers of CYS.

After this is completed, go to 

 Tools > Structure editing > Add H

to add hydrogens to the ligand. Then, go to 

 Tools > Structure editing > Add Charge > (have Amber ff14SB and AM1-BCC selected) -> Ok

and then go to 

 File -> Save Mol2...

and save as 6me2_lig_dockprep.mol2. Close Chimera.

Generating Simulation Parameters

In order to generate parameters for this tutorial, switch to the following directory: 

 cd 002_parameters

To generate the parameters, run the following command: 

 antechamber -i ../001_structure/6me2_lig_wH.mol2 -fi mol2 -o 6me2_ligand_antechamber.mol2 -fo mol2 -at gaff2 -c bcc -rn LIG -nc 0

Note the -nc flag is set to 0 in the above line; this should correspond to the protonation state of the ligand. If a user is following this tutorial to simulate dynamics on a structure other than 6ME2, be sure to double check and change this value accordingly.

After the 6me2_ligand_antechamber.mol2 output file is generated, run the following command to run parmch2: 

 parmchk2 -i 6me2_ligand_antechamber.mol2 -f mol2 -o 6me2_ligand.am1bcc.frcmod

Using TLeap

Now, switch to the following directory: 

 cd 003_leap

Where files will be generated to simulate the ligand and receptor, such that further calculations involving both structures--such as binding affinities, free energies, and trajectories. Two types of files will be generated that can be used to set up the system: parameter (parm7) and restart (rst7) files. Start the process by creating the following input file: 

 vi leap.in

And put the following into the input 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 ../001_structure/6me2_fresh.pdb
 ###load ligand frcmod/mol2
 loadamberparams ../002_parameters/6me2_ligand.am1bcc.frcmod
 lig=loadmol2 ../002_parameters/1HW9_ligand_antechamber.mol2
 ###create gase-phase complex
 gascomplex= combine {rec lig}
 ###write gas-phase pdb
 savepdb gascomplex 6me2.gas.complex.pdb
 ###write gase-phase toplogy and coord files for MMGBSA calc
 saveamberparm gascomplex 6me2.complex.parm7 6me2.gas.complex.rst7
 saveamberparm rec 6me2.gas.receptor.parm7 6me2.gas.receptor.rst7
 saveamberparm lig 6me2.gas.ligand.parm7 6me2.gas.ligand.rst7
 ###create solvated complex (albeit redundant)
 solvcomplex= combine {rec lig}
 ###solvate the system
 solvateoct solvcomplex TIP3PBOX 12.0
 ###Neutralize system
 addions solvcomplex Cl- 0
 addions solvcomplex Na+ 0
 #write solvated pdb file
 savepdb solvcomplex 6me2.wet.complex.pdb
 ###check the system
 charge solvcomplex
 check solvcomplex
 ###write solvated toplogy and coordinate file
 saveamberparm solvcomplex 6me2.wet.complex.parm7 6me2.wet.complex.rst7
 quit

After this runs, copy the parm7 (parameter) and rst7 (coordinate) files to the local machine, where they can be opened in Chimera. To open in Chimera, go to 

 Tools -> MD/Ensemble Analysis -> MD Movie --> Choose appropriate parameter and coordinate files to load

Equilibrating the system

Now, move to the following directory: 

 cd 004_equil

In this process, energy minimization and equilibration will be performed on the system. There are 9 input files that need to be generated for this file, as follows: 

 vi 01.min.mdin

With the following input: 

 Minimize 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="!@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
 /

Next, 

 vi 02.equil.mdin

with input 

 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=":!@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  
 /

Next, 

 vi 03.min.mdin

with input 

 Minimize 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="!@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
 /

Next, 

 vi 04.min.mdin

with input 

 Minimize 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="!@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
 /

Next, 

 vi 05.min.mdin

with input 

 Minimize 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="!@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 
 /

Next, 

 vi 06.equil.mdin

with input 

 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="!@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
 /

Next, 

 vi 07.equil.mdin

with input 

 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="!@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
 /

Next, 

 vi 08.equil.mdin *See below for specific comments on this input

with input 

 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-433@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
 /

And finally, 

 vi 09.equil.mdin *See below for comments on this input

with input 

 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-433@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
 /

Note that the restraint mask in files 8 and 9 have to be specified to the relevant system being simulated. These are residue numbers which can be found in the wet.complex.pdb file generated in the 003_leap directory.

To submit all 9 input files as a job on Seawulf, create the following script: 

 vi mdequilibration.sh

with input 

 #!/bin/sh
 #SBATCH --job-name=6me2_equilibration
 #SBATCH --ntasks-per-node=40
 #SBATCH --nodes=2
 #SBATCH --time=8:00:00
 #SBATCH -p long-40core  
 echo "started Equilibration on 'date' "
 parm7="../003_leap/6me2.wet.complex.parm7"
 coords="../003_leap/6me2.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
 pmemd -O -i ${input}.mdin -o ${input}.mdout -p $parm7 -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf ${input}.info -r 
 ${input}.rst7 
 coords=$input
 done
 echo " Finished equilibration on 'date' "

And submit this to the queue: 

 sbatch mdequilibration.sh

MD Production

Now, enter the following directory: 

 cd 005_production

and create the following production input file: 

 vi 10.prod.mdin

with inputs 

 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-432@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 
 /

The restraint mask should be modified on this file as well to only include the receptor, not the ligand (i.e., it should be one fewer than the restraint mask value in the 8th input file from the previous directory).

Create a file that will be submitted to Seawulf: 

 vi mdproduction.sh

with the following: 

 !/bin/sh
 #SBATCH --job-name=6me2_prod
 #SBATCH --ntasks-per-node=24
 #SBATCH --nodes=2
 #SBATCH --time=7-00:00:00
 #SBATCH -p extended-24core
 echo "started Equilibration on 'date' "
 parm7="../003_leap/6me2.wet.complex.parm7"
 coords="../004_equil/09.equil"
 MDINPUTS=( 10.prod)
 for input in ${MDINPUTS[@]}; do
   pmemd -O -i ${input}.mdin -o ${input}.mdout -p ${parm7} -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf ${input}.info -r 
 ${input}.rst7
   coords=$input
 done
 echo "Finished Equilibration on `date` "

And now submit this job to the queue: 

 sbatch mdproduction.sh

And note that this job will likely take a very long time, perhaps days.

MD Analysis

Since the trajectories have been generated, they can now be inspected using RMSD calculations, hydrogen bond analysis, and MMGBSA calculations.

RMSD Calculation

To perform the RMSD calculations, first create an input file which will delete the water and charges from the production trajectory. 

 vi cpptraj_strip_wat.in

with the following inputs 

 #!/usr/bin/sh
 parm ../003_leap/6me2.wet.complex.parm7
 #read in trajectory
 trajin ../004_production/10.prod.trj
 #read in reference
 reference ../003_leap/6me2.wet.complex.rst7
 #compute RMSD + align CA to crystal structure
 rmsd rms1 reference :1-263@CA
 #strip solvent
 strip :WAT:Na+:Cl-
 #create gas-phase trajectory
 trajout 6me2_stripfit_restrained_gas.trj nobox

The reference mask should contain both the ligand and receptor in its value.

Now run this file using the following command: 

 cpptraj -i cpptraj_strip_wat.in

Next, create another input file which will be used to calculate the RMSD of the ligand by typing 

 vi cpptraj_rmsd_lig.in

with inputs 

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

And check that the reference mask only factors in the ligand in this file. Run this input file using the following command: 

 cpptraj -p ../003_leap/6me2.complex.parm7 -i cpptraj_rmsd_lig.in

Next, the receptor RMSD can be calculated by creating the following input file: 

 vi cpptraj_rmsd_rec.in

with inputs 

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

And check that the reference mask only takes into account the receptor. Now run this file using the following command: 

 cpptraj -p ../003_leap/6me2.gas.complex.parm7 -i cpptraj_rmsd_rec.in

Calculating Number of Hydrogen Bonds

Next, the number of hydrogen bonds will be analyzed to see if any exist between the receptor and ligand. Create the following input file: 

 vi cpptraj_hbond.in

with inputs 

 !/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-263 out 6me2_calcipotriol_hbond.out avgout 6me2_calcipotriol_hbond_avg.dat solventdonor :WAT solventacceptor :WAT@-O
 nointramol brid\
 geout 6me2_calcipotriol_bridge-water.dat dist 3.0 angle 140

Be sure that the hbobnd hb1 parameter corresponds to the residue numbers of the relevant ligand and receptor. Now, ru this input file with the following command: 

 cpptraj -p ../003_leap/6me2_wetcomplex.parm7 -i cpptraj_hbond.in

MMGBSA Calculation

Next, MMGBSA will be used to calculate free energies of binding in this system. First, create the following input file: 

 vi mmgbsa.in

with inputs 

 mmgbsa 6me2 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,
  /

To run this on Seawulf, create the following script: 

 vi mmgbsa.sh

and input the following 

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

And finally, submit this as a job to the queue by typing 

 sbatch mmgbsa.sh
Retrieved from "https://ringo.ams.stonybrook.edu/index.php?title=2022_AMBER_tutorial_1_with_PDBID_6ME2&oldid=15026"