Difference between revisions of "2018 AMBER tutorial with 1c87"
From Rizzo_Lab
m |
|||
| (2 intermediate revisions by the same user not shown) | |||
| Line 423: | Line 423: | ||
We will now remove the water molecules and create a trajectory aligned to the CA atoms (backbone) of the crystal structure. In order to do this, we can create "cpptraj.strip.wat.in" file : | We will now remove the water molecules and create a trajectory aligned to the CA atoms (backbone) of the crystal structure. In order to do this, we can create "cpptraj.strip.wat.in" file : | ||
| + | With the RMSD values, we are able to observe how far the ligand and receptor are moved. | ||
#!/usr/bin/sh | #!/usr/bin/sh | ||
| Line 457: | Line 458: | ||
==H-Bond== | ==H-Bond== | ||
| + | Move to 002.hbond directory. We can find out hydrogen bonding between the ligand and the enzyme residues. | ||
| + | We are going to make a bash script: | ||
| + | #!/bin/sh | ||
| + | #PBS -N 1c87_production | ||
| + | #PBS -l walltime=07:00:00 | ||
| + | #PBS -l nodes=1:ppn=28 | ||
| + | #PBS -j oe | ||
| + | #PBS -q long | ||
| + | #PBS -V | ||
| + | |||
| + | cd $PBS_O_WORKDIR | ||
| + | |||
| + | echo "Started Analysis on `date` " | ||
| + | #do_parallel="mpirun pmemd.MPI" | ||
| + | |||
| + | cpptraj -p ../../001.tleap_build/1c87.wet.complex.prmtop -i cpptraj.hbond.in | ||
| + | |||
| + | echo "Finished Analysis on `date` " | ||
| + | |||
| + | Once this script is create, make sure to ''cpptraj'' in the command line to see if cpptraj is working fine. | ||
| + | Then, make an input file named 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-297 out 1c87.hbond.out avgout 1c87_sunitinib.hbond.avg.dat solventdonor :WAT solventacceptor :WAT@O nointramol bridgeout 1c87_sunitinib.bridge-water.dat dist 3.0 angle 140 | ||
| + | |||
| + | Submit the job. | ||
| + | qsub hbond.sh | ||
| + | |||
| + | Once the job is completed, ''1c87.hbond.out, 1c87_sunitinib.bridge-water.dat, 1c87_sunitinib.hbond.avg.dat'' are created. | ||
| + | We can take a look at the h-bonds with ligand. | ||
| + | [[File:Screen Shot 2018-04-04 at 4.21.40 PM.png|thumb|center|1000px| hydrogen bond formation]] | ||
| + | |||
| + | |||
==MM-GBSA== | ==MM-GBSA== | ||
| + | Move to 003.mmgbsa directory. | ||
| + | MM-GBSA values tell us the relative binding free energy of ligand-receptor complex and we are able to see how strongly they bind to each other. | ||
| + | |||
| + | We are also going to create the bash script for MM-GBSA calculation, '''mmgbsa.sh'': | ||
| + | |||
| + | #!/bin/bash | ||
| + | #PBS -l walltime=35:00:00 | ||
| + | #PBS -l nodes=1:ppn=24 | ||
| + | #PBS -N 1c87_mmgbsa | ||
| + | #PBS -V | ||
| + | #PBS -j oe | ||
| + | #PBS -q long-24core | ||
| + | |||
| + | cd $PBS_O_WORKDIR | ||
| + | |||
| + | #Define topology files | ||
| + | solv_prmtop="../../001.tleap_build/1c87.wet.complex.prmtop" | ||
| + | complex_prmtop="../../001.tleap_build/1c87.gas.complex.prmtop" | ||
| + | receptor_prmtop="../../001.tleap_build/1c87.gas.receptor.prmtop" | ||
| + | ligand_prmtop="../../001.tleap_build/1c87.gas.ligand.prmtop" | ||
| + | trajectory="../../003.production/001.restrained/10.prod.trj" | ||
| + | |||
| + | |||
| + | #create mmgbsa input file | ||
| + | cat >mmgbsa.in<<EOF | ||
| + | mmgbsa 1c87 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 1c87.mmgbsa.results.dat \ | ||
| + | -eo 1c87.mmgbsa.per-frame.dat \ | ||
| + | -sp ${solv_prmtop} \ | ||
| + | -cp ${complex_prmtop} \ | ||
| + | -rp ${receptor_prmtop} \ | ||
| + | -lp ${ligand_prmtop} \ | ||
| + | -y ${trajectory} | ||
| + | |||
| + | Then, make an input file, ''mmgbsa.in''. | ||
| + | |||
| + | mmgbsa 1c87 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, | ||
| + | / | ||
| + | |||
| + | Submit the job. | ||
| + | qsub mmgbsa.sh | ||
| + | |||
| + | The following output files will be created when the job is completed. | ||
| + | ''_MMPBSA_ligand.nc.0'' | ||
| + | ''_MMPBSA_ligand_nm.nc.0'' | ||
| + | ''_MMPBSA_com_nm_traj_cpptraj.out _MMPBSA_ligand.pdb'' | ||
| + | ''_MMPBSA_complex_gb.mdout.0 _MMPBSA_lig_nm_traj_cpptraj.out'' | ||
| + | ''_MMPBSA_complex.nc.0 _MMPBSA_normal_traj_cpptraj.out'' | ||
| + | ''_MMPBSA_complex_nm.nc.0 _MMPBSA_receptor.nc.0'' | ||
| + | ''_MMPBSA_complex.pdb _MMPBSA_receptor_nm.nc.0'' | ||
| + | ''_MMPBSA_dummycomplex.inpcrd _MMPBSA_receptor.pdb'' | ||
| + | ''_MMPBSA_dummyligand.inpcrd _MMPBSA_rec_nm_traj_cpptraj.out'' | ||
| + | ''_MMPBSA_dummyreceptor.inpcrd reference.frc'' | ||
| + | ''_MMPBSA_gb.mdin'' | ||
Latest revision as of 15:26, 4 April 2018
In this tutorial, we will learn how to run a MD simulation of a protein-ligand complex and learn how to calculate structural fluctuations (with RMSD) and free energies of binding (MM-GBSA).
Contents
I. Prepare the files and directories
We can create the directories to make all the files organized.
~username/AMS536/2018/Amber_Tutorial/000.parameters
/001.tleap_build
/002.equilibration
/003.production
/004.analysis
/zzz.master
Once the directories are created, we are going to convert 1c87_lig_withH_charged.mol2 file and 1c87_rec_withH_charged.mol2 file to pdb file in chimera, then copy the converted files into zzz.master directory.
II. Structural Preparation
Before we submit our jobs for equilibration or production, we must parametrize the ligands by using antechamber and check for any missing force field parameters. Make sure we have the proper partial charge. Since there are two acids in the 1c87 structure, -nc will be -2.
cd 000.parameters # Parameterize the ligand antechamber -i ../zzz.master/2nnq.lig.withH.charged.for_antechamber.pdb -fi pdb -o 2nnq_lig.am1bcc.mol2 -fo mol2 -at gaff2 -c bcc -rn LIG -nc -2 # Check for missing FF parameters parmchk2 -i 2nnq_lig.am1bcc.mol2 -f mol2 -o 2nnq_lig.am1bcc.frcmod
Then, we will create the topology and coordinate files for the simulation through tleap. Create a file called tleap.in with the following:
#!/usr/bin/sh
###Load Protein force field
source leaprc.protein.ff14SB
###Load GAFF force field (for our ligand)
source leaprc.gaff
###Load TIP3P (water) force field
source leaprc.water.tip3p
####Load Ions frcmod for the tip3p model
loadamberparams frcmod.ionsjc_tip3p
###Needed so we can use igb=8 model
set default PBradii mbondi3
###Load Protein pdb file
rec=loadpdb ../zzz.master/rec_withH_1c87.pdb
###Load Ligand frcmod/mol2
loadamberparams ../000.parameters/1c87_lig.am1bcc.frcmod
lig=loadmol2 ../000.parameters/1c87_lig.am1bcc.mol2
###Create gas-phase complex
gascomplex= combine {rec lig}
###Write gas-phase pdb
savepdb gascomplex 1c87.gas.complex.pdb
###Write gas-phase toplogy and coord files for MMGBSA calc
saveamberparm gascomplex 1c87.gas.complex.prmtop 1c87.gas.complex.rst7
saveamberparm rec 1c87.gas.receptor.prmtop 1c87.gas.receptor.rst7
saveamberparm lig 1c87.gas.ligand.prmtop 1c87.gas.ligand.rst7
###Create solvated complex (albeit redundant)
solvcomplex= combine {rec lig}
###Solvate the systemsolvateoct
solvcomplex TIP3PBOX 12.0
###Neutralize system (it will add either Na or Cl depending on net charge)
addions solvcomplex Cl- 0
addions solvcomplex Na+ 0
###Write solvated pdb file
savepdb solvcomplex 1c87.wet.complex.pdb
###Check the system
charge solvcomplex
check solvcomplex
###Write Solvated topology and coord file
saveamberparm solvcomplex 1c87.wet.complex.prmtop 1c87.wet.complex.rst7
Create Amber topology and coordinates files for the MD simulation.
tleap -f tleap.in
III. Equilibration
We can now create the input files for the equilibration of the system in 002.equilbration directory. We are going to run minimization in order to remove steric clashes and allow the receptor and ligand to equilibrate into a low energy state (MD simulation).
01.min.mdin:
Minmize all the hydrogens &cntrl imin=1, ! Minimize the initial structure maxcyc=5000, ! Maximum number of cycles for minimization ntb=1, ! Constant volume ntp=0, ! No pressure scaling ntf=1, ! Complete force evaluation ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-298 & !@H=", ! atoms to be restrained restraint_wt=5.0, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format /
02.equil.mdin:
MD simualation &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntb=2, ! Constant Pressure ntc=1, ! No SHAKE on bonds between hydrogens dt=0.001, ! Timestep (ps) ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-298 & !@H=", ! atoms to be restrained restraint_wt=5.0, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
03.min.mdin:
Minmize all the hydrogens &cntrl imin=1, ! Minimize the initial structure maxcyc=1000, ! Maximum number of cycles for minimization ntb=1, ! Constant volume ntp=0, ! No pressure scaling ntf=1, ! Complete force evaluation ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-298 & !@H=", ! atoms to be restrained restraint_wt=2.0, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format /
04.min.mdin:
Minmize all the hydrogens &cntrl imin=1, ! Minimize the initial structure maxcyc=1000, ! Maximum number of cycles for minimization ntb=1, ! Constant volume ntp=0, ! No pressure scaling ntf=1, ! Complete force evaluation ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-298 & !@H=", ! atoms to be restrained restraint_wt=0.1, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format /
05.min.mdin:
Minmize all the hydrogens &cntrl imin=1, ! Minimize the initial structure maxcyc=1000, ! Maximum number of cycles for minimization ntb=1, ! Constant volume ntp=0, ! No pressure scaling ntf=1, ! Complete force evaluation ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-298 & !@H=", ! atoms to be restrained restraint_wt=0.05, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format /
06.equil.mdin:
MD simualation &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntb=2, ! Constant Pressure ntc=1, ! No SHAKE on bonds between hydrogens dt=0.001, ! Timestep (ps) ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-298 & !@H=", ! atoms to be restrained restraint_wt=1.0, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
07.equil.mdin:
MD simulation &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntx=5, ! Positions and velocities read formatted irest=1, ! Restart calculation ntc=1, ! No SHAKE on for bonds with hydrogen dt=0.001, ! Timestep (ps) ntb=2, ! Constant Pressure ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-298 & !@H=", ! atoms to be restrained restraint_wt=0.5, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
08.equil.mdin:
MD simulations &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntx=5, ! Positions and velocities read formatted irest=1, ! Restart calculation ntc=1, ! No SHAKE on for bonds with hydrogen dt=0.001, ! Timestep (ps) ntb=2, ! Constant Pressure ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-298 & !@H=", ! atoms to be restrained restraint_wt=0.1, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
09.equil.mdin:
MD simulations &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntx=5, ! Positions and velocities read formatted irest=1, ! Restart calculation ntc=1, ! No SHAKE on for bonds with hydrogen dt=0.001, ! Timestep (ps) ntb=2, ! Constant Pressure ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-298 & !@H=", ! atoms to be restrained restraint_wt=0.1, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
Create a submission script md.equilibration.sh in order to start equilibration. We are going to run it for 24 hours in order to have enough time to complete all equilibration jobs.
vi md.equilibration.sh
#!/bin/sh
#PBS -N 1c87_equilibration
#PBS -l walltime=24:00:00
#PBS -l nodes=2:ppn=28
#PBS -j oe
#PBS -q long
cd $PBS_O_WORKDIR
echo "Started Equilibration on `date` "
do_parallel="sander"
prmtop="../001.tleap_build/1c87.wet.complex.prmtop"
coords="../001.tleap_build/1c87.wet.complex"
MDINPUTS=(01.min 02.equil 03.min 04.min 05.min 06.equil 07.equil 08.equil 09.equil)
for input in ${MDINPUTS[@]}; do
$do_parallel -O -i ${input}.mdin -o ${input}.mdout -p $prmtop -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf ${input}.info -r ${input}.rst7
coords=$input
done
echo "Finished Equilibration on `date` "
Submit the job
qsub md.equilibration.sh
Once submitted the job, info, rst7, trj files will be created for each input file. We can use VMD to open up the prmtop file and trajectory file. Determine file type: AMBER coordinates.
IV. Production
Move to 003.production directory and open up 001.restrained directory. We are going to create the file for running production. Create the input file for MD simulation of production.
vi 10.prod.mdin
MD simulations &cntrl imin=0, ! Perform MD nstlim=5000000, ! Number of MD steps ntx=5, ! Positions and velocities read formatted irest=1, ! Restart calculation ntc=2, ! SHAKE on for bonds with hydrogen dt=0.002, ! Timestep (ps) ntb=2, ! Constant Pressure ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=2, ! No force evaluation for bonds with hydrogen ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 2500, ! Write to trajectory file every ntwx steps ntpr= 2500, ! Print to mdout every ntpr steps ntwr= 5000000, ! Write a restart file every ntwr steps cut=8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-298@CA,C,N", ! atoms to be restrained restraint_wt=0.1, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
We are going to create the bash script.
vi prod.sh
#!/bin/sh
#PBS -N 1c87_production
#PBS -l walltime=24:00:00
#PBS -l nodes=2:ppn=28
#PBS -j oe
#PBS -q long
cd $PBS_O_WORKDIR
echo "Started Production on `date` "
#do_parallel="mpirun pmemd.MPI"
do_cuda="mpirun pmemd.MPI"
prmtop="/gpfs/projects/AMS536/2018/seung/AMBER_directory_arch/001.tleap_build/1c87.wet.complex.prmtop"
coords="../../002.equilbration/07.equil"
MDINPUTS=(10.prod)
for input in ${MDINPUTS[@]}; do
$do_cuda -O -i ${input}.mdin -o ${input}.mdout -p $prmtop -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf ${input}.info -r ${input}.rst7
coords=$input
done
echo "Finished Production on `date` "
Submit the job.
qsub prod.sh
V. Simulation Analysis
RMSD
After completing production simulations, we can now start analyzing and apply post processing. We can create three directories:
mkdir 001.rmsd 002.hbond 003.mmgbsa
We will now remove the water molecules and create a trajectory aligned to the CA atoms (backbone) of the crystal structure. In order to do this, we can create "cpptraj.strip.wat.in" file : With the RMSD values, we are able to observe how far the ligand and receptor are moved.
#!/usr/bin/sh #read in trajectory trajin /gpfs/projects/AMS536/2018/seung/AMBER_directory_arch/003.production/001.restrained/10.prod.trj #read in reference reference /gpfs/projects/AMS536/2018/seung/AMBER_directory_arch/001.tleap_build/1c87.wet.complex.rst7 #compute rmsd and align CA to the crystal structure rmsd rms1 reference :1-298@CA #strip Solvent strip :WAT:Na+:Cl- #create gas-phase trajectory trajout 1c87.stripfit.restrained.gas.trj nobox
and we are going to run it:
cpptraj -p ../../001.tleap_build/1c87_sunitinib.wet.complex.prmtop -i cpptraj.strip.wat.in
We can take a look at the trajectory file in VMD.
Next is to compute ligand rmsd vs. the crystal structure pose and histogram the RMSD.
#!/usr/bin/sh #trajin the trajectory trajin 1c87.stripfit.restrained.gas.trj #read in the reference reference ../../001.tleap_build/1c87.gas.complex.rst7 #compute the RMSD (do not fit the internal geometries first, included rigid body motions #and convert the frames to ns (framenum*.005) rmsd rms1 ":298&!(@H=)" nofit mass out 1c87.lig.restrained.rmsd.nofit.dat time .005 #histogram the nofit rmsd histogram rms1,*,*,.1,* norm out 1c87.lig.restrained.rmsd.nofit.histogram.dat
and again, we will run it:
cpptraj -p ../../001.tleap_build/1c87_sunitinib.gas.complex.prmtop -i cpptraj.rmsd.lig.in
H-Bond
Move to 002.hbond directory. We can find out hydrogen bonding between the ligand and the enzyme residues. We are going to make a bash script:
#!/bin/sh #PBS -N 1c87_production #PBS -l walltime=07:00:00 #PBS -l nodes=1:ppn=28 #PBS -j oe #PBS -q long #PBS -V cd $PBS_O_WORKDIR echo "Started Analysis on `date` " #do_parallel="mpirun pmemd.MPI" cpptraj -p ../../001.tleap_build/1c87.wet.complex.prmtop -i cpptraj.hbond.in echo "Finished Analysis on `date` "
Once this script is create, make sure to cpptraj in the command line to see if cpptraj is working fine. Then, make an input file named 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-297 out 1c87.hbond.out avgout 1c87_sunitinib.hbond.avg.dat solventdonor :WAT solventacceptor :WAT@O nointramol bridgeout 1c87_sunitinib.bridge-water.dat dist 3.0 angle 140
Submit the job.
qsub hbond.sh
Once the job is completed, 1c87.hbond.out, 1c87_sunitinib.bridge-water.dat, 1c87_sunitinib.hbond.avg.dat are created. We can take a look at the h-bonds with ligand.
MM-GBSA
Move to 003.mmgbsa directory. MM-GBSA values tell us the relative binding free energy of ligand-receptor complex and we are able to see how strongly they bind to each other.
We are also going to create the bash script for MM-GBSA calculation, 'mmgbsa.sh:
#!/bin/bash
#PBS -l walltime=35:00:00
#PBS -l nodes=1:ppn=24
#PBS -N 1c87_mmgbsa
#PBS -V
#PBS -j oe
#PBS -q long-24core
cd $PBS_O_WORKDIR
#Define topology files
solv_prmtop="../../001.tleap_build/1c87.wet.complex.prmtop"
complex_prmtop="../../001.tleap_build/1c87.gas.complex.prmtop"
receptor_prmtop="../../001.tleap_build/1c87.gas.receptor.prmtop"
ligand_prmtop="../../001.tleap_build/1c87.gas.ligand.prmtop"
trajectory="../../003.production/001.restrained/10.prod.trj"
#create mmgbsa input file
cat >mmgbsa.in<<EOF
mmgbsa 1c87 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 1c87.mmgbsa.results.dat \
-eo 1c87.mmgbsa.per-frame.dat \
-sp ${solv_prmtop} \
-cp ${complex_prmtop} \
-rp ${receptor_prmtop} \
-lp ${ligand_prmtop} \
-y ${trajectory}
Then, make an input file, mmgbsa.in.
mmgbsa 1c87 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, /
Submit the job.
qsub mmgbsa.sh
The following output files will be created when the job is completed. _MMPBSA_ligand.nc.0 _MMPBSA_ligand_nm.nc.0 _MMPBSA_com_nm_traj_cpptraj.out _MMPBSA_ligand.pdb _MMPBSA_complex_gb.mdout.0 _MMPBSA_lig_nm_traj_cpptraj.out _MMPBSA_complex.nc.0 _MMPBSA_normal_traj_cpptraj.out _MMPBSA_complex_nm.nc.0 _MMPBSA_receptor.nc.0 _MMPBSA_complex.pdb _MMPBSA_receptor_nm.nc.0 _MMPBSA_dummycomplex.inpcrd _MMPBSA_receptor.pdb _MMPBSA_dummyligand.inpcrd _MMPBSA_rec_nm_traj_cpptraj.out _MMPBSA_dummyreceptor.inpcrd reference.frc _MMPBSA_gb.mdin
