Difference between revisions of "2019 AMBER tutorial with PDBID 2BXF"
From Rizzo_Lab
Stonybrook (talk | contribs) (→TLeap) |
Stonybrook (talk | contribs) (→H Bond) |
||
| (12 intermediate revisions by the same user not shown) | |||
| Line 29: | Line 29: | ||
Create tleap.build.in file | Create tleap.build.in file | ||
| − | + | #!/usr/bin/sh | |
###Load Protein force field | ###Load Protein force field | ||
| Line 41: | Line 41: | ||
###Needed so we can use igb=8 model | ###Needed so we can use igb=8 model | ||
set default PBradii mbondi3 | set default PBradii mbondi3 | ||
| − | |||
###Load Protein pdb file | ###Load Protein pdb file | ||
rec=loadpdb ../zzz.master/backup.pdb | rec=loadpdb ../zzz.master/backup.pdb | ||
| Line 61: | Line 60: | ||
bond rec.554.SG rec.563.SG | bond rec.554.SG rec.563.SG | ||
bond rec.356.SG rec.365.SG | bond rec.356.SG rec.365.SG | ||
| − | |||
###Load Ligand frcmod/mol2 | ###Load Ligand frcmod/mol2 | ||
loadamberparams ../000.parameters/2BXF_lig.am1bcc.frcmod | loadamberparams ../000.parameters/2BXF_lig.am1bcc.frcmod | ||
lig=loadmol2 ../000.parameters/2BXF_lig.am1bcc.mol2 | lig=loadmol2 ../000.parameters/2BXF_lig.am1bcc.mol2 | ||
| − | |||
###Create gas-phase complex | ###Create gas-phase complex | ||
gascomplex= combine {rec lig} | gascomplex= combine {rec lig} | ||
| − | |||
###Write gas-phase pdb | ###Write gas-phase pdb | ||
savepdb gascomplex 2BXF.gas.complex.pdb | savepdb gascomplex 2BXF.gas.complex.pdb | ||
| − | |||
###Write gas-phase toplogy and coord files for MMGBSA calc | ###Write gas-phase toplogy and coord files for MMGBSA calc | ||
saveamberparm gascomplex 2BXF.gas.complex.prmtop 2BXF.gas.complex.rst7 | saveamberparm gascomplex 2BXF.gas.complex.prmtop 2BXF.gas.complex.rst7 | ||
saveamberparm rec 2BXF.gas.receptor.prmtop 2BXF.gas.receptor.rst7 | saveamberparm rec 2BXF.gas.receptor.prmtop 2BXF.gas.receptor.rst7 | ||
saveamberparm lig 2BXF.gas.ligand.prmtop 2BXF.gas.ligand.rst7 | saveamberparm lig 2BXF.gas.ligand.prmtop 2BXF.gas.ligand.rst7 | ||
| − | |||
###Create solvated complex (albeit redundant) | ###Create solvated complex (albeit redundant) | ||
solvcomplex= combine {rec lig} | solvcomplex= combine {rec lig} | ||
| − | |||
###Solvate the system | ###Solvate the system | ||
solvateoct solvcomplex TIP3PBOX 12.0 | solvateoct solvcomplex TIP3PBOX 12.0 | ||
| − | |||
###Neutralize system (it will add either Na or Cl depending on net charge) | ###Neutralize system (it will add either Na or Cl depending on net charge) | ||
addions solvcomplex Cl- 0 | addions solvcomplex Cl- 0 | ||
addions solvcomplex Na+ 0 | addions solvcomplex Na+ 0 | ||
| − | |||
###Write solvated pdb file | ###Write solvated pdb file | ||
savepdb solvcomplex 2BXF.wet.complex.pdb | savepdb solvcomplex 2BXF.wet.complex.pdb | ||
| − | |||
###Check the system | ###Check the system | ||
charge solvcomplex | charge solvcomplex | ||
check solvcomplex | check solvcomplex | ||
| − | |||
###Write Solvated topology and coord file | ###Write Solvated topology and coord file | ||
saveamberparm solvcomplex 2BXF.wet.complex.prmtop 2BXF.wet.complex.rst7 | saveamberparm solvcomplex 2BXF.wet.complex.prmtop 2BXF.wet.complex.rst7 | ||
| − | |||
| − | |||
tleap -f tleap.build.in | tleap -f tleap.build.in | ||
| − | + | = Equilibration = | |
| + | |||
| + | Create the following input files for the equilibration of the system. | ||
| + | vim 01.min.mdin | ||
| + | |||
| + | Minmize all the hydrogens | ||
| + | &cntrl | ||
| + | imin=1, ! Minimize the initial structure | ||
| + | ntmin=2, ! Use steepest descent Ryota Added | ||
| + | maxcyc=5000, ! Maximum number of cycles for minimization | ||
| + | ntb=1, ! Constant volume | ||
| + | ntp=0, ! No pressure scaling | ||
| + | ntf=1, ! Complete force evaluation | ||
| + | ntwx= 1000, ! Write to trajectory file every ntwx steps | ||
| + | ntpr= 1000, ! Print to mdout every ntpr steps | ||
| + | ntwr= 1000, ! Write a restart file every ntwr steps | ||
| + | cut= 8.0, ! Nonbonded cutoff in Angstroms | ||
| + | ntr=1, ! Turn on restraints | ||
| + | restraintmask=":1-579 & !@H=", ! atoms to be restrained | ||
| + | restraint_wt=5.0, ! force constant for restraint | ||
| + | ntxo=1, ! Write coordinate file in ASCII format | ||
| + | ioutfm=0, ! Write trajectory file in ASCII format | ||
| + | / | ||
| + | |||
| + | vim 02.equil.mdin | ||
| + | |||
| + | MD simualation | ||
| + | &cntrl | ||
| + | imin=0, ! Perform MD | ||
| + | nstlim=50000 ! Number of MD steps | ||
| + | ntb=2, ! Constant Pressure | ||
| + | ntc=1, ! No SHAKE on bonds between hydrogens | ||
| + | dt=0.001, ! Timestep (ps) | ||
| + | ntp=1, ! Isotropic pressure scaling | ||
| + | barostat=1 ! Berendsen | ||
| + | taup=0.5 ! Pressure relaxtion time (ps) | ||
| + | ntf=1, ! Complete force evaluation | ||
| + | ntt=3, ! Langevin thermostat | ||
| + | gamma_ln=2.0 ! Collision Frequency for thermostat | ||
| + | ig=-1, ! Random seed for thermostat | ||
| + | temp0=298.15 ! Simulation temperature (K) | ||
| + | ntwx= 1000, ! Write to trajectory file every ntwx steps | ||
| + | ntpr= 1000, ! Print to mdout every ntpr steps | ||
| + | ntwr= 1000, ! Write a restart file every ntwr steps | ||
| + | cut= 8.0, ! Nonbonded cutoff in Angstroms | ||
| + | ntr=1, ! Turn on restraints | ||
| + | restraintmask=":1-579 & !@H=", ! atoms to be restrained | ||
| + | restraint_wt=5.0, ! force constant for restraint | ||
| + | ntxo=1, ! Write coordinate file in ASCII format | ||
| + | ioutfm=0, ! Write trajectory file in ASCII format | ||
| + | iwrap=1, ! iwrap is turned on | ||
| + | / | ||
| + | |||
| + | vim 03.min.mdin | ||
| + | |||
| + | Minmize all the hydrogens | ||
| + | &cntrl | ||
| + | imin=1, ! Minimize the initial structure | ||
| + | maxcyc=1000, ! Maximum number of cycles for minimization | ||
| + | ntb=1, ! Constant volume | ||
| + | ntp=0, ! No pressure scaling | ||
| + | ntf=1, ! Complete force evaluation | ||
| + | ntwx= 1000, ! Write to trajectory file every ntwx steps | ||
| + | ntpr= 1000, ! Print to mdout every ntpr steps | ||
| + | ntwr= 1000, ! Write a restart file every ntwr steps | ||
| + | cut= 8.0, ! Nonbonded cutoff in Angstroms | ||
| + | ntr=1, ! Turn on restraints | ||
| + | restraintmask=":1-579 & !@H=", ! atoms to be restrained | ||
| + | restraint_wt=2.0, ! force constant for restraint | ||
| + | ntxo=1, ! Write coordinate file in ASCII format | ||
| + | ioutfm=0, ! Write trajectory file in ASCII format | ||
| + | / | ||
| + | |||
| + | vim 04.min.mdin | ||
| + | |||
| + | Minmize all the hydrogens | ||
| + | &cntrl | ||
| + | imin=1, ! Minimize the initial structure | ||
| + | maxcyc=1000, ! Maximum number of cycles for minimization | ||
| + | ntb=1, ! Constant volume | ||
| + | ntp=0, ! No pressure scaling | ||
| + | ntf=1, ! Complete force evaluation | ||
| + | ntwx= 1000, ! Write to trajectory file every ntwx steps | ||
| + | ntpr= 1000, ! Print to mdout every ntpr steps | ||
| + | ntwr= 1000, ! Write a restart file every ntwr steps | ||
| + | cut= 8.0, ! Nonbonded cutoff in Angstroms | ||
| + | ntr=1, ! Turn on restraints | ||
| + | restraintmask=":1-579 & !@H=", ! atoms to be restrained | ||
| + | restraint_wt=0.1, ! force constant for restraint | ||
| + | ntxo=1, ! Write coordinate file in ASCII format | ||
| + | ioutfm=0, ! Write trajectory file in ASCII format | ||
| + | / | ||
| + | |||
| + | vim 05.min.mdin | ||
| + | |||
| + | Minmize all the hydrogens | ||
| + | &cntrl | ||
| + | imin=1, ! Minimize the initial structure | ||
| + | maxcyc=1000, ! Maximum number of cycles for minimization | ||
| + | ntb=1, ! Constant volume | ||
| + | ntp=0, ! No pressure scaling | ||
| + | ntf=1, ! Complete force evaluation | ||
| + | ntwx= 1000, ! Write to trajectory file every ntwx steps | ||
| + | ntpr= 1000, ! Print to mdout every ntpr steps | ||
| + | ntwr= 1000, ! Write a restart file every ntwr steps | ||
| + | cut= 8.0, ! Nonbonded cutoff in Angstroms | ||
| + | ntr=1, ! Turn on restraints | ||
| + | restraintmask=":1-579 & !@H=", ! atoms to be restrained | ||
| + | restraint_wt=0.05, ! force constant for restraint | ||
| + | ntxo=1, ! Write coordinate file in ASCII format | ||
| + | ioutfm=0, ! Write trajectory file in ASCII format | ||
| + | / | ||
| − | + | vim 06.equil.mdin | |
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| + | MD simualation | ||
| + | &cntrl | ||
| + | imin=0, ! Perform MD | ||
| + | nstlim=50000 ! Number of MD steps | ||
| + | ntb=2, ! Constant Pressure | ||
| + | ntc=1, ! No SHAKE on bonds between hydrogens | ||
| + | dt=0.001, ! Timestep (ps) | ||
| + | ntp=1, ! Isotropic pressure scaling | ||
| + | barostat=1 ! Berendsen | ||
| + | taup=0.5 ! Pressure relaxtion time (ps) | ||
| + | ntf=1, ! Complete force evaluation | ||
| + | ntt=3, ! Langevin thermostat | ||
| + | gamma_ln=2.0 ! Collision Frequency for thermostat | ||
| + | ig=-1, ! Random seed for thermostat | ||
| + | temp0=298.15 ! Simulation temperature (K) | ||
| + | ntwx= 1000, ! Write to trajectory file every ntwx steps | ||
| + | ntpr= 1000, ! Print to mdout every ntpr steps | ||
| + | ntwr= 1000, ! Write a restart file every ntwr steps | ||
| + | cut= 8.0, ! Nonbonded cutoff in Angstroms | ||
| + | ntr=1, ! Turn on restraints | ||
| + | restraintmask=":1-579 & !@H=", ! atoms to be restrained | ||
| + | restraint_wt=1.0, ! force constant for restraint | ||
| + | ntxo=1, ! Write coordinate file in ASCII format | ||
| + | ioutfm=0, ! Write trajectory file in ASCII format | ||
| + | iwrap=1, ! iwrap is turned on | ||
| + | / | ||
| − | + | vim 07.equil.mdin | |
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | + | MD simulation | |
| − | + | &cntrl | |
| − | + | imin=0, ! Perform MD | |
| + | nstlim=50000 ! Number of MD steps | ||
| + | ntx=5, ! Positions and velocities read formatted | ||
| + | irest=1, ! Restart calculation | ||
| + | ntc=1, ! No SHAKE on for bonds with hydrogen | ||
| + | dt=0.001, ! Timestep (ps) | ||
| + | ntb=2, ! Constant Pressure | ||
| + | ntp=1, ! Isotropic pressure scaling | ||
| + | barostat=1 ! Berendsen | ||
| + | taup=0.5 ! Pressure relaxtion time (ps) | ||
| + | ntf=1, ! Complete force evaluation | ||
| + | ntt=3, ! Langevin thermostat | ||
| + | gamma_ln=2.0 ! Collision Frequency for thermostat | ||
| + | ig=-1, ! Random seed for thermostat | ||
| + | temp0=298.15 ! Simulation temperature (K) | ||
| + | ntwx= 1000, ! Write to trajectory file every ntwx steps | ||
| + | ntpr= 1000, ! Print to mdout every ntpr steps | ||
| + | ntwr= 1000, ! Write a restart file every ntwr steps | ||
| + | cut= 8.0, ! Nonbonded cutoff in Angstroms | ||
| + | ntr=1, ! Turn on restraints | ||
| + | restraintmask=":1-579 & !@H=", ! atoms to be restrained | ||
| + | restraint_wt=0.5, ! force constant for restraint | ||
| + | ntxo=1, ! Write coordinate file in ASCII format | ||
| + | ioutfm=0, ! Write trajectory file in ASCII format | ||
| + | iwrap=1, ! iwrap is turned on | ||
| + | / | ||
| − | + | vim 08.equil.mdin | |
| − | |||
| − | + | MD simulations | |
| − | + | &cntrl | |
| + | imin=0, ! Perform MD | ||
| + | nstlim=50000 ! Number of MD steps | ||
| + | ntx=5, ! Positions and velocities read formatted | ||
| + | irest=1, ! Restart calculation | ||
| + | ntc=1, ! No SHAKE on for bonds with hydrogen | ||
| + | dt=0.001, ! Timestep (ps) | ||
| + | ntb=2, ! Constant Pressure | ||
| + | ntp=1, ! Isotropic pressure scaling | ||
| + | barostat=1 ! Berendsen | ||
| + | taup=0.5 ! Pressure relaxtion time (ps) | ||
| + | ntf=1, ! Complete force evaluation | ||
| + | ntt=3, ! Langevin thermostat | ||
| + | gamma_ln=2.0 ! Collision Frequency for thermostat | ||
| + | ig=-1, ! Random seed for thermostat | ||
| + | temp0=298.15 ! Simulation temperature (K) | ||
| + | ntwx= 1000, ! Write to trajectory file every ntwx steps | ||
| + | ntpr= 1000, ! Print to mdout every ntpr steps | ||
| + | ntwr= 1000, ! Write a restart file every ntwr steps | ||
| + | cut= 8.0, ! Nonbonded cutoff in Angstroms | ||
| + | ntr=1, ! Turn on restraints | ||
| + | restraintmask=":1-578@CA,C,N", ! atoms to be restrained | ||
| + | restraint_wt=0.1, ! force constant for restraint | ||
| + | ntxo=1, ! Write coordinate file in ASCII format | ||
| + | ioutfm=0, ! Write trajectory file in ASCII format | ||
| + | iwrap=1, ! iwrap is turned on | ||
| + | / | ||
| − | + | vim 09.equil.mdin | |
| − | + | ||
| − | |||
| − | |||
| − | + | MD simulations | |
| − | + | &cntrl | |
| + | imin=0, ! Perform MD | ||
| + | nstlim=50000 ! Number of MD steps | ||
| + | ntx=5, ! Positions and velocities read formatted | ||
| + | irest=1, ! Restart calculation | ||
| + | ntc=1, ! No SHAKE on for bonds with hydrogen | ||
| + | dt=0.001, ! Timestep (ps) | ||
| + | ntb=2, ! Constant Pressure | ||
| + | ntp=1, ! Isotropic pressure scaling | ||
| + | barostat=1 ! Berendsen | ||
| + | taup=0.5 ! Pressure relaxtion time (ps) | ||
| + | ntf=1, ! Complete force evaluation | ||
| + | ntt=3, ! Langevin thermostat | ||
| + | gamma_ln=2.0 ! Collision Frequency for thermostat | ||
| + | ig=-1, ! Random seed for thermostat | ||
| + | temp0=298.15 ! Simulation temperature (K) | ||
| + | ntwx= 1000, ! Write to trajectory file every ntwx steps | ||
| + | ntpr= 1000, ! Print to mdout every ntpr steps | ||
| + | ntwr= 1000, ! Write a restart file every ntwr steps | ||
| + | cut= 8.0, ! Nonbonded cutoff in Angstroms | ||
| + | ntr=1, ! Turn on restraints | ||
| + | restraintmask=":1-578@CA,C,N", ! atoms to be restrained | ||
| + | restraint_wt=0.1, ! force constant for restraint | ||
| + | ntxo=1, ! Write coordinate file in ASCII format | ||
| + | ioutfm=0, ! Write trajectory file in ASCII format | ||
| + | iwrap=1, ! iwrap is turned on | ||
| + | / | ||
| − | + | Create a submission script that utilizes, above input files to equilibrate the biological system. | |
| − | + | vim mdequilibration.sh | |
| − | ### | + | #!/bin/sh |
| − | + | #PBS -N 2BXF_equilibration | |
| − | + | #PBS -l walltime=04:00:00 | |
| + | #PBS -l nodes=2:ppn=28 | ||
| + | #PBS -j oe | ||
| + | #PBS -q long | ||
| + | |||
| + | cd $PBS_O_WORKDIR | ||
| + | |||
| + | echo "Started Equilibration on `date` " | ||
| + | do_parallel="sander" | ||
| + | |||
| + | prmtop="../001.tleap_build/2BXF.wet.complex.prmtop" | ||
| + | coords="../001.tleap_build/2BXF.wet.complex" | ||
| + | |||
| + | |||
| + | MDINPUTS=(01.min 02.equil 03.min 04.min 05.min 06.equil 07.equil 08.equil 09.equil) | ||
| + | |||
| + | for input in ${MDINPUTS[@]}; do | ||
| + | |||
| + | $do_parallel -O -i ${input}.mdin -o ${input}.mdout -p $prmtop -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf | ||
| + | ${input}.info -r ${input}.rst7 | ||
| + | coords=$input | ||
| + | done | ||
| + | |||
| + | echo "Finished Equilibration on `date` " | ||
| − | + | Submit the job | |
| − | |||
| − | + | qsub md.equilibration.sh | |
| − | |||
| − | |||
| − | + | = Production = | |
| − | |||
| − | + | Move into 003.production | |
| − | + | Move into 001.restrained | |
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| − | |||
| + | Create 10.prod.mdin | ||
| − | + | MD simulations | |
| − | + | &cntrl | |
| − | + | imin=0, ! Perform MD | |
| − | + | nstlim=5000000, ! Number of MD steps | |
| − | + | ntx=5, ! Positions and velocities read formatted | |
| − | + | irest=1, ! Restart calculation | |
| − | + | ntc=2, ! SHAKE on for bonds with hydrogen | |
| − | + | dt=0.002, ! Timestep (ps) | |
| − | + | ntb=2, ! Constant Pressure | |
| − | + | ntp=1, ! Isotropic pressure scaling | |
| − | + | barostat=1 ! Berendsen | |
| − | + | taup=0.5 ! Pressure relaxtion time (ps) | |
| − | + | ntf=2, ! No force evaluation for bonds with hydrogen | |
| − | + | ntt=3, ! Langevin thermostat | |
| − | + | gamma_ln=2.0 ! Collision Frequency for thermostat | |
| − | + | ig=-1, ! Random seed for thermostat | |
| − | + | temp0=298.15 ! Simulation temperature (K) | |
| − | + | ntwx= 2500, ! Write to trajectory file every ntwx steps | |
| − | + | ntpr= 2500, ! Print to mdout every ntpr steps | |
| + | ntwr= 5000000, ! Write a restart file every ntwr steps | ||
| + | cut=8.0, ! Nonbonded cutoff in Angstroms | ||
| + | ntr=1, ! Turn on restraints | ||
| + | restraintmask=":1-579@CA,C,N", ! atoms to be restrained | ||
| + | restraint_wt=0.1, ! force constant for restraint | ||
| + | ntxo=1, ! Write coordinate file in ASCII format | ||
| + | ioutfm=0, ! Write trajectory file in ASCII format | ||
| + | iwrap=1, ! iwrap is turned on | ||
| − | ### | + | Create prod.sh |
| − | + | ||
| − | + | #!/bin/sh | |
| + | #PBS -N 2BXF_production | ||
| + | #PBS -l walltime=24:00:00 | ||
| + | #PBS -l nodes=2:ppn=24 | ||
| + | #PBS -j oe | ||
| + | #PBS -o production.output | ||
| + | #PBS -q long | ||
| + | |||
| + | cd $PBS_O_WORKDIR | ||
| + | |||
| + | echo "Started Production on `date` " | ||
| + | #do_parallel="mpirun pmemd.MPI" | ||
| + | do_cuda="mpirun pmemd.MPI" | ||
| + | |||
| + | prmtop="../../001.tleap_build/2BXF.wet.complex.prmtop" | ||
| + | coords="../../002.equilbration/09.equil" | ||
| + | |||
| + | |||
| + | MDINPUTS=(10.prod) | ||
| + | |||
| + | for input in ${MDINPUTS[@]}; do | ||
| + | |||
| + | $do_cuda -O -i ${input}.mdin -o ${input}.mdout -p $prmtop -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf | ||
| + | ${input}.info -r ${input}.rst7 | ||
| + | coords=$input | ||
| + | done | ||
| + | |||
| + | echo "Finished Production on `date` " | ||
| − | |||
| − | |||
| − | + | Submit the job | |
| − | + | qsub prod.sh | |
| − | + | Move into 002.unrestrained | |
| − | |||
| − | |||
| − | |||
| − | + | Create 10.prod.mdin | |
| − | |||
| − | + | MD simulations | |
| − | + | &cntrl | |
| + | imin=0, ! Perform MD | ||
| + | nstlim=5000000 ! Number of MD steps | ||
| + | ntx=5, ! Positions and velocities read formatted | ||
| + | irest=1, ! Restart calculation | ||
| + | ntc=2, ! SHAKE on for bonds with hydrogen | ||
| + | dt=0.002, ! Timestep (ps) | ||
| + | ntb=2, ! Constant Pressure | ||
| + | ntp=1, ! Isotropic pressure scaling | ||
| + | barostat=1 ! Berendsen | ||
| + | taup=0.5 ! Pressure relaxtion time (ps) | ||
| + | ntf=2, ! No force evaluation for bonds with hydrogen | ||
| + | ntt=3, ! Langevin thermostat | ||
| + | gamma_ln=2.0 ! Collision Frequency for thermostat | ||
| + | ig=-1, ! Random seed for thermostat | ||
| + | temp0=298.15 ! Simulation temperature (K) | ||
| + | ntwx= 2500, ! Write to trajectory file every ntwx steps | ||
| + | ntpr= 2500, ! Print to mdout every ntpr steps | ||
| + | ntwr= 5000000, ! Write a restart file every ntwr steps | ||
| + | cut= 8.0, ! Nonbonded cutoff in Angstroms | ||
| + | ntxo=1, ! Write coordinate file in ASCII format | ||
| + | ioutfm=0, ! Write trajectory file in ASCII format | ||
| + | iwrap=1, ! iwrap is turned on | ||
| + | / | ||
| − | |||
| − | |||
| − | |||
| − | ### | + | Create prod.sh |
| − | + | ||
| + | #!/bin/sh | ||
| + | #PBS -N 2BXF_production | ||
| + | #PBS -l walltime=24:00:00 | ||
| + | #PBS -l nodes=2:ppn=24 | ||
| + | #PBS -j oe | ||
| + | #PBS -o production.output | ||
| + | #PBS -q long | ||
| + | |||
| + | cd $PBS_O_WORKDIR | ||
| + | |||
| + | echo "Started Production on `date` " | ||
| + | #do_parallel="mpirun pmemd.MPI" | ||
| + | do_cuda="mpirun pmemd.MPI" | ||
| + | |||
| + | prmtop="../../001.tleap_build/2BXF.wet.complex.prmtop" | ||
| + | coords="../../002.equilbration/09.equil" | ||
| + | |||
| + | |||
| + | MDINPUTS=(10.prod) | ||
| + | |||
| + | for input in ${MDINPUTS[@]}; do | ||
| + | |||
| + | $do_cuda -O -i ${input}.mdin -o ${input}.mdout -p $prmtop -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf | ||
| + | ${input}.info -r ${input}.rst7 | ||
| + | coords=$input | ||
| + | done | ||
| + | |||
| + | echo "Finished Production on `date` " | ||
| − | |||
| − | |||
| − | |||
| − | + | Submit the job | |
| − | + | qsub prod.sh | |
= Analysis = | = Analysis = | ||
| Line 254: | Line 510: | ||
reference ../../001.tleap_build/2BXF.wet.complex.rst7 | reference ../../001.tleap_build/2BXF.wet.complex.rst7 | ||
#compute rmsd and align CA to the crystal structure | #compute rmsd and align CA to the crystal structure | ||
| − | rmsd rms1 reference :1- | + | rmsd rms1 reference :1-578@CA |
#strip Solvent | #strip Solvent | ||
strip :WAT:Na+:Cl- | strip :WAT:Na+:Cl- | ||
| Line 271: | Line 527: | ||
#compute the RMSD (do not fit the internal geometries first, included rigid body motions | #compute the RMSD (do not fit the internal geometries first, included rigid body motions | ||
#and convert the frames to ns (framenum*.005) | #and convert the frames to ns (framenum*.005) | ||
| − | rmsd rms1 ": | + | rmsd rms1 ":579&!(@H=)" nofit mass out 2BXF.lig.restrained.rmsd.nofit.dat time .005 |
#histogram the nofit rmsd | #histogram the nofit rmsd | ||
histogram rms1,*,*,.1,* norm out 2BXF.lig.restrained.rmsd.nofit.histogram.dat | histogram rms1,*,*,.1,* norm out 2BXF.lig.restrained.rmsd.nofit.histogram.dat | ||
| Line 287: | Line 543: | ||
#compute the RMSD (do not fit the internal geometries first, included rigid body motions | #compute the RMSD (do not fit the internal geometries first, included rigid body motions | ||
#and convert the frames to ns (framenum*.005) | #and convert the frames to ns (framenum*.005) | ||
| − | rmsd rms1 ":1- | + | rmsd rms1 ":1-578&!(@H=)" nofit mass out 2BXF.rec.restrained.rmsd.nofit.dat time .005 |
#histogram the nofit rmsd | #histogram the nofit rmsd | ||
histogram rms1,*,*,.1,* norm out 2BXF.rec.restrained.rmsd.nofit.histogram.dat | histogram rms1,*,*,.1,* norm out 2BXF.rec.restrained.rmsd.nofit.histogram.dat | ||
| Line 305: | Line 561: | ||
#autoimage | #autoimage | ||
#compute intra and water mediated hydrogen bonds | #compute intra and water mediated hydrogen bonds | ||
| − | hbond hb1 :1- | + | hbond hb1 :1-579 out 2BXF_sunitinib.hbond.out avgout 2BXF_sunitinib.hbond.avg.dat solventdonor :WAT solventacceptor :WAT@O |
nointramol brid\ | nointramol brid\ | ||
geout 2BXF_sunitinib.bridge-water.dat dist 3.0 angle 140 | geout 2BXF_sunitinib.bridge-water.dat dist 3.0 angle 140 | ||
Latest revision as of 15:14, 19 April 2019
2BXF with an explicit solvent model
Contents
Prepare the files
Convert 2BXF.lig.withH.charged.mol2 to pdb in chimera
Convert 2BXF.rec.withH.charged.mol2 to pdb in chimera
Copy into zzz.master
Parameters
The system we are working on has two main components (Protein receptor & ligand). The usual forcefield "ff14SB" contains all the parameters needed for calculations of the protein. However, the ligand is a non-protein component. Therefore, ff14SB forcefield does not contain the parameters needed for the calculations regarding the ligand. Therefore, we need to generate parameters needed for the ligand. The following steps will be taken using antechamber in order to generate the ligand parameters.
Move into 000.programs
Paramaterize the ligand
antechamber -i ../zzz.master/2BXF.lig.withH.charged.pdb -fi pdb -o 2BXF_lig.am1bcc.mol2 -fo mol2 -at gaff2 -c bcc -rn LIG -nc -1
Check for missing force field parameters
parmchk2 -i 2BXF_lig.am1bcc.mol2 -f mol2 -o 2BXF_lig.am1bcc.frcmod
TLeap
Move into 001.tleap_build
Create tleap.build.in file
#!/usr/bin/sh
###Load Protein force field
source leaprc.protein.ff14SB
###Load GAFF force field (for our ligand)
source leaprc.gaff
###Load TIP3P (water) force field
source leaprc.water.tip3p
####Load Ions frcmod for the tip3p model
loadamberparams frcmod.ionsjc_tip3p
###Needed so we can use igb=8 model
set default PBradii mbondi3
###Load Protein pdb file
rec=loadpdb ../zzz.master/backup.pdb
bond rec.49.SG rec.58.SG
bond rec.241.SG rec.249.SG
bond rec.71.SG rec.87.SG
bond rec.86.SG rec.97.SG
bond rec.120.SG rec.165.SG
bond rec.164.SG rec.173.SG
bond rec.510.SG rec.555.SG
bond rec.472.SG rec.482.SG
bond rec.457.SG rec.473.SG
bond rec.196.SG rec.242.SG
bond rec.261.SG rec.275.SG
bond rec.312.SG rec.357.SG
bond rec.274.SG rec.285.SG
bond rec.388.SG rec.434.SG
bond rec.433.SG rec.444.SG
bond rec.554.SG rec.563.SG
bond rec.356.SG rec.365.SG
###Load Ligand frcmod/mol2
loadamberparams ../000.parameters/2BXF_lig.am1bcc.frcmod
lig=loadmol2 ../000.parameters/2BXF_lig.am1bcc.mol2
###Create gas-phase complex
gascomplex= combine {rec lig}
###Write gas-phase pdb
savepdb gascomplex 2BXF.gas.complex.pdb
###Write gas-phase toplogy and coord files for MMGBSA calc
saveamberparm gascomplex 2BXF.gas.complex.prmtop 2BXF.gas.complex.rst7
saveamberparm rec 2BXF.gas.receptor.prmtop 2BXF.gas.receptor.rst7
saveamberparm lig 2BXF.gas.ligand.prmtop 2BXF.gas.ligand.rst7
###Create solvated complex (albeit redundant)
solvcomplex= combine {rec lig}
###Solvate the system
solvateoct solvcomplex TIP3PBOX 12.0
###Neutralize system (it will add either Na or Cl depending on net charge)
addions solvcomplex Cl- 0
addions solvcomplex Na+ 0
###Write solvated pdb file
savepdb solvcomplex 2BXF.wet.complex.pdb
###Check the system
charge solvcomplex
check solvcomplex
###Write Solvated topology and coord file
saveamberparm solvcomplex 2BXF.wet.complex.prmtop 2BXF.wet.complex.rst7
tleap -f tleap.build.in
Equilibration
Create the following input files for the equilibration of the system.
vim 01.min.mdin
Minmize all the hydrogens &cntrl imin=1, ! Minimize the initial structure ntmin=2, ! Use steepest descent Ryota Added maxcyc=5000, ! Maximum number of cycles for minimization ntb=1, ! Constant volume ntp=0, ! No pressure scaling ntf=1, ! Complete force evaluation ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-579 & !@H=", ! atoms to be restrained restraint_wt=5.0, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format /
vim 02.equil.mdin
MD simualation &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntb=2, ! Constant Pressure ntc=1, ! No SHAKE on bonds between hydrogens dt=0.001, ! Timestep (ps) ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-579 & !@H=", ! atoms to be restrained restraint_wt=5.0, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
vim 03.min.mdin
Minmize all the hydrogens &cntrl imin=1, ! Minimize the initial structure maxcyc=1000, ! Maximum number of cycles for minimization ntb=1, ! Constant volume ntp=0, ! No pressure scaling ntf=1, ! Complete force evaluation ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-579 & !@H=", ! atoms to be restrained restraint_wt=2.0, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format /
vim 04.min.mdin
Minmize all the hydrogens &cntrl imin=1, ! Minimize the initial structure maxcyc=1000, ! Maximum number of cycles for minimization ntb=1, ! Constant volume ntp=0, ! No pressure scaling ntf=1, ! Complete force evaluation ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-579 & !@H=", ! atoms to be restrained restraint_wt=0.1, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format /
vim 05.min.mdin Minmize all the hydrogens &cntrl imin=1, ! Minimize the initial structure maxcyc=1000, ! Maximum number of cycles for minimization ntb=1, ! Constant volume ntp=0, ! No pressure scaling ntf=1, ! Complete force evaluation ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-579 & !@H=", ! atoms to be restrained restraint_wt=0.05, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format /
vim 06.equil.mdin
MD simualation &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntb=2, ! Constant Pressure ntc=1, ! No SHAKE on bonds between hydrogens dt=0.001, ! Timestep (ps) ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-579 & !@H=", ! atoms to be restrained restraint_wt=1.0, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
vim 07.equil.mdin
MD simulation &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntx=5, ! Positions and velocities read formatted irest=1, ! Restart calculation ntc=1, ! No SHAKE on for bonds with hydrogen dt=0.001, ! Timestep (ps) ntb=2, ! Constant Pressure ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-579 & !@H=", ! atoms to be restrained restraint_wt=0.5, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
vim 08.equil.mdin
MD simulations &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntx=5, ! Positions and velocities read formatted irest=1, ! Restart calculation ntc=1, ! No SHAKE on for bonds with hydrogen dt=0.001, ! Timestep (ps) ntb=2, ! Constant Pressure ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-578@CA,C,N", ! atoms to be restrained restraint_wt=0.1, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
vim 09.equil.mdin
MD simulations &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntx=5, ! Positions and velocities read formatted irest=1, ! Restart calculation ntc=1, ! No SHAKE on for bonds with hydrogen dt=0.001, ! Timestep (ps) ntb=2, ! Constant Pressure ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-578@CA,C,N", ! atoms to be restrained restraint_wt=0.1, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
Create a submission script that utilizes, above input files to equilibrate the biological system.
vim mdequilibration.sh
#!/bin/sh
#PBS -N 2BXF_equilibration
#PBS -l walltime=04:00:00
#PBS -l nodes=2:ppn=28
#PBS -j oe
#PBS -q long
cd $PBS_O_WORKDIR
echo "Started Equilibration on `date` "
do_parallel="sander"
prmtop="../001.tleap_build/2BXF.wet.complex.prmtop"
coords="../001.tleap_build/2BXF.wet.complex"
MDINPUTS=(01.min 02.equil 03.min 04.min 05.min 06.equil 07.equil 08.equil 09.equil)
for input in ${MDINPUTS[@]}; do
$do_parallel -O -i ${input}.mdin -o ${input}.mdout -p $prmtop -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf
${input}.info -r ${input}.rst7
coords=$input
done
echo "Finished Equilibration on `date` "
Submit the job
qsub md.equilibration.sh
Production
Move into 003.production
Move into 001.restrained
Create 10.prod.mdin
MD simulations &cntrl imin=0, ! Perform MD nstlim=5000000, ! Number of MD steps ntx=5, ! Positions and velocities read formatted irest=1, ! Restart calculation ntc=2, ! SHAKE on for bonds with hydrogen dt=0.002, ! Timestep (ps) ntb=2, ! Constant Pressure ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=2, ! No force evaluation for bonds with hydrogen ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 2500, ! Write to trajectory file every ntwx steps ntpr= 2500, ! Print to mdout every ntpr steps ntwr= 5000000, ! Write a restart file every ntwr steps cut=8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-579@CA,C,N", ! atoms to be restrained restraint_wt=0.1, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on
Create prod.sh
#!/bin/sh
#PBS -N 2BXF_production
#PBS -l walltime=24:00:00
#PBS -l nodes=2:ppn=24
#PBS -j oe
#PBS -o production.output
#PBS -q long
cd $PBS_O_WORKDIR
echo "Started Production on `date` "
#do_parallel="mpirun pmemd.MPI"
do_cuda="mpirun pmemd.MPI"
prmtop="../../001.tleap_build/2BXF.wet.complex.prmtop"
coords="../../002.equilbration/09.equil"
MDINPUTS=(10.prod)
for input in ${MDINPUTS[@]}; do
$do_cuda -O -i ${input}.mdin -o ${input}.mdout -p $prmtop -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf
${input}.info -r ${input}.rst7
coords=$input
done
echo "Finished Production on `date` "
Submit the job
qsub prod.sh
Move into 002.unrestrained
Create 10.prod.mdin
MD simulations &cntrl imin=0, ! Perform MD nstlim=5000000 ! Number of MD steps ntx=5, ! Positions and velocities read formatted irest=1, ! Restart calculation ntc=2, ! SHAKE on for bonds with hydrogen dt=0.002, ! Timestep (ps) ntb=2, ! Constant Pressure ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=2, ! No force evaluation for bonds with hydrogen ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 2500, ! Write to trajectory file every ntwx steps ntpr= 2500, ! Print to mdout every ntpr steps ntwr= 5000000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
Create prod.sh
#!/bin/sh
#PBS -N 2BXF_production
#PBS -l walltime=24:00:00
#PBS -l nodes=2:ppn=24
#PBS -j oe
#PBS -o production.output
#PBS -q long
cd $PBS_O_WORKDIR
echo "Started Production on `date` "
#do_parallel="mpirun pmemd.MPI"
do_cuda="mpirun pmemd.MPI"
prmtop="../../001.tleap_build/2BXF.wet.complex.prmtop"
coords="../../002.equilbration/09.equil"
MDINPUTS=(10.prod)
for input in ${MDINPUTS[@]}; do
$do_cuda -O -i ${input}.mdin -o ${input}.mdout -p $prmtop -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf
${input}.info -r ${input}.rst7
coords=$input
done
echo "Finished Production on `date` "
Submit the job
qsub prod.sh
Analysis
RMSD
How far our ligand moved and how far our receptor moved
Create cpptraj.strip.wat.in
#!/usr/bin/sh parm ../../001.tleap_build/2BXF.wet.complex.prmtop #read in trajectory trajin ../../003.production/001.restrained/10.prod.trj #read in reference reference ../../001.tleap_build/2BXF.wet.complex.rst7 #compute rmsd and align CA to the crystal structure rmsd rms1 reference :1-578@CA #strip Solvent strip :WAT:Na+:Cl- #create gas-phase trajectory trajout 2BXF.stripfit.restrained.gas.trj nobox
Run cpptraj
cpptraj -i cpptraj.strip.wat.in
Create cpptraj.rmsd.lig.in
#!/usr/bin/sh #trajin the trajectory trajin 2BXF.stripfit.restrained.gas.trj #read in the reference reference ../../001.tleap_build/2BXF.gas.complex.rst7 #compute the RMSD (do not fit the internal geometries first, included rigid body motions #and convert the frames to ns (framenum*.005) rmsd rms1 ":579&!(@H=)" nofit mass out 2BXF.lig.restrained.rmsd.nofit.dat time .005 #histogram the nofit rmsd histogram rms1,*,*,.1,* norm out 2BXF.lig.restrained.rmsd.nofit.histogram.dat
Run cpptraj
cpptraj -p ../../001.tleap_build/2BXF.gas.complex.prmtop -i cpptraj.rmsd.lig.in
Create cpptraj.rmsd.rec.in
#!/usr/bin/sh #trajin the trajectory trajin 2BXF.stripfit.restrained.gas.trj #read in the reference reference ../../001.tleap_build/2BXF.gas.complex.rst7 #compute the RMSD (do not fit the internal geometries first, included rigid body motions #and convert the frames to ns (framenum*.005) rmsd rms1 ":1-578&!(@H=)" nofit mass out 2BXF.rec.restrained.rmsd.nofit.dat time .005 #histogram the nofit rmsd histogram rms1,*,*,.1,* norm out 2BXF.rec.restrained.rmsd.nofit.histogram.dat
Run cpptraj
cpptraj -p ../../001.tleap_build/2BXF.gas.complex.prmtop -i cpptraj.rmsd.rec.in
H Bond
Hydrogen bonding
Create cpptraj.hbond.in
#!/usr/bin/sh #read in trajectory trajin ../../003.production/001.restrained/10.prod.trj #wrap everything into one periodic cell #autoimage #compute intra and water mediated hydrogen bonds hbond hb1 :1-579 out 2BXF_sunitinib.hbond.out avgout 2BXF_sunitinib.hbond.avg.dat solventdonor :WAT solventacceptor :WAT@O nointramol brid\ geout 2BXF_sunitinib.bridge-water.dat dist 3.0 angle 140
Run cpptraj
cpptraj -p ../../001.tleap_build/2BXF.wet.complex.prmtop -i cpptraj.hbond.in
MMGBSA
This will analyze how strongly our small molecule binds to our receptor.
Create mmgbsa.in
mmgbsa 2BXF analysis
&general
interval=1, netcdf=1,
keep_files=0,
/
&gb
igb=8,
saltcon=0.0, surften=0.0072,
surfoff=0.0, molsurf=0,
/
&nmode
drms=0.001, maxcyc=10000,
nminterval=250, nmendframe=2000,
nmode_igb=1,
/
Create mmgbsa.sh
#!/bin/bash
#PBS -l walltime=35:00:00
#PBS -l nodes=1:ppn=24
#PBS -N 4qmz_mmgbsa
#PBS -V
#PBS -j oe
#PBS -q long-24core
cd $PBS_O_WORKDIR
#Define topology files
solv_prmtop="../../001.tleap_build/2BXF.wet.complex.prmtop"
complex_prmtop="../../001.tleap_build/2BXF.gas.complex.prmtop"
receptor_prmtop="../../001.tleap_build/2BXF.gas.receptor.prmtop"
ligand_prmtop="../../001.tleap_build/2BXF.gas.ligand.prmtop"
trajectory="../../003.production/001.restrained/10.prod.trj"
#create mmgbsa input file
cat >mmgbsa.in<<EOF
mmgbsa HIVgp41 analysis
&general
interval=1, netcdf=1,
keep_files=0,
/
&gb
igb=8,
saltcon=0.0, surften=0.0072,
surfoff=0.0, molsurf=0,
/
&nmode
drms=0.001, maxcyc=10000,
nminterval=250, nmendframe=2000,
nmode_igb=1,
/
EOF
MMPBSA.py -O -i mmgbsa.in \
-o 2BXF.mmgbsa.results.dat \
-eo 2BXF.mmgbsa.per-frame.dat \
-sp ${solv_prmtop} \
-cp ${complex_prmtop} \
-rp ${receptor_prmtop} \
-lp ${ligand_prmtop} \
-y ${trajectory}
Submit your script
qsub mmgbsa.sh
The last line of our 2BXF.mmgbsa.results.dat file shows that our delta G binding is -24.4431 +/- 5.4851.
