View source for 2019 AMBER tutorial with PDBID 2BXF

2BXF with an explicit solvent model

= Prepare the files =

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

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

Copy into zzz.master

= Parameters =

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

Move into 000.programs

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

Check for missing force field parameters
 
 parmchk2 -i 2BXF_lig.am1bcc.mol2 -f mol2 -o 2BXF_lig.am1bcc.frcmod

= TLeap =

Move into 001.tleap_build

Create tleap.build.in file

 #!/usr/bin/sh

 ###Load Protein force field
 source leaprc.protein.ff14SB
 ###Load GAFF force field (for our ligand)
 source leaprc.gaff
 ###Load TIP3P (water) force field
 source leaprc.water.tip3p
 ####Load Ions frcmod for the tip3p model
 loadamberparams frcmod.ionsjc_tip3p
 ###Needed so we can use igb=8 model
 set default PBradii mbondi3
 ###Load Protein pdb file
 rec=loadpdb ../zzz.master/backup.pdb
 bond rec.49.SG  rec.58.SG
 bond rec.241.SG rec.249.SG
 bond rec.71.SG  rec.87.SG
 bond rec.86.SG  rec.97.SG
 bond rec.120.SG rec.165.SG
 bond rec.164.SG rec.173.SG
 bond rec.510.SG rec.555.SG
 bond rec.472.SG rec.482.SG
 bond rec.457.SG rec.473.SG
 bond rec.196.SG rec.242.SG
 bond rec.261.SG rec.275.SG
 bond rec.312.SG rec.357.SG
 bond rec.274.SG rec.285.SG
 bond rec.388.SG rec.434.SG
 bond rec.433.SG rec.444.SG
 bond rec.554.SG rec.563.SG
 bond rec.356.SG rec.365.SG
 ###Load Ligand frcmod/mol2
 loadamberparams ../000.parameters/2BXF_lig.am1bcc.frcmod
 lig=loadmol2 ../000.parameters/2BXF_lig.am1bcc.mol2
 ###Create gas-phase complex
 gascomplex= combine {rec lig}
 ###Write gas-phase pdb
 savepdb gascomplex 2BXF.gas.complex.pdb
 ###Write gas-phase toplogy and coord files for MMGBSA calc
 saveamberparm gascomplex 2BXF.gas.complex.prmtop 2BXF.gas.complex.rst7
 saveamberparm rec 2BXF.gas.receptor.prmtop 2BXF.gas.receptor.rst7
 saveamberparm lig 2BXF.gas.ligand.prmtop 2BXF.gas.ligand.rst7
 ###Create solvated complex (albeit redundant)
 solvcomplex= combine {rec lig}
 ###Solvate the system
 solvateoct solvcomplex TIP3PBOX 12.0
 ###Neutralize system (it will add either Na or Cl depending on net charge)
 addions solvcomplex Cl- 0
 addions solvcomplex Na+ 0
 ###Write solvated pdb file
 savepdb solvcomplex 2BXF.wet.complex.pdb
 ###Check the system
 charge solvcomplex
 check solvcomplex
 ###Write Solvated topology and coord file
 saveamberparm solvcomplex 2BXF.wet.complex.prmtop 2BXF.wet.complex.rst7


 tleap -f tleap.build.in

= Equilibration =

Create the following input files for the equilibration of the system.
 vim 01.min.mdin

 Minmize all the hydrogens
 &cntrl
  imin=1,           ! Minimize the initial structure
  ntmin=2,         ! Use steepest descent Ryota Added
  maxcyc=5000,    ! Maximum number of cycles for minimization
  ntb=1,            ! Constant volume
  ntp=0,            ! No pressure scaling
  ntf=1,            ! Complete force evaluation
  ntwx= 1000,       ! Write to trajectory file every ntwx steps
  ntpr= 1000,       ! Print to mdout every ntpr steps
  ntwr= 1000,       ! Write a restart file every ntwr steps
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
  ntr=1,            ! Turn on restraints
  restraintmask=":1-579 & !@H=", ! atoms to be restrained
  restraint_wt=5.0, ! force constant for restraint
  ntxo=1,           ! Write coordinate file in ASCII format
  ioutfm=0,         ! Write trajectory file in ASCII format
 /

 vim 02.equil.mdin

  MD simualation
 &cntrl
  imin=0,           ! Perform MD
  nstlim=50000      ! Number of MD steps
  ntb=2,            ! Constant Pressure
  ntc=1,            ! No SHAKE on bonds between hydrogens
  dt=0.001,         ! Timestep (ps)
  ntp=1,            ! Isotropic pressure scaling
  barostat=1        ! Berendsen
  taup=0.5          ! Pressure relaxtion time (ps)
  ntf=1,            ! Complete force evaluation
  ntt=3,            ! Langevin thermostat
  gamma_ln=2.0      ! Collision Frequency for thermostat
  ig=-1,            ! Random seed for thermostat
  temp0=298.15      ! Simulation temperature (K)
  ntwx= 1000,       ! Write to trajectory file every ntwx steps
  ntpr= 1000,       ! Print to mdout every ntpr steps
  ntwr= 1000,       ! Write a restart file every ntwr steps
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
  ntr=1,            ! Turn on restraints
  restraintmask=":1-579 & !@H=", ! atoms to be restrained
  restraint_wt=5.0, ! force constant for restraint
  ntxo=1,           ! Write coordinate file in ASCII format
  ioutfm=0,         ! Write trajectory file in ASCII format
  iwrap=1,          ! iwrap is turned on
 /

 vim 03.min.mdin

 Minmize all the hydrogens
 &cntrl
  imin=1,           ! Minimize the initial structure
  maxcyc=1000,    ! Maximum number of cycles for minimization
  ntb=1,            ! Constant volume
  ntp=0,            ! No pressure scaling
  ntf=1,            ! Complete force evaluation
  ntwx= 1000,       ! Write to trajectory file every ntwx steps
  ntpr= 1000,       ! Print to mdout every ntpr steps
  ntwr= 1000,       ! Write a restart file every ntwr steps
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
  ntr=1,            ! Turn on restraints
  restraintmask=":1-579 & !@H=", ! atoms to be restrained
  restraint_wt=2.0, ! force constant for restraint
  ntxo=1,           ! Write coordinate file in ASCII format
  ioutfm=0,         ! Write trajectory file in ASCII format
 /

 vim 04.min.mdin

  Minmize all the hydrogens
 &cntrl
  imin=1,           ! Minimize the initial structure
  maxcyc=1000,    ! Maximum number of cycles for minimization
  ntb=1,            ! Constant volume
  ntp=0,            ! No pressure scaling
  ntf=1,            ! Complete force evaluation
  ntwx= 1000,       ! Write to trajectory file every ntwx steps
  ntpr= 1000,       ! Print to mdout every ntpr steps
  ntwr= 1000,       ! Write a restart file every ntwr steps
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
  ntr=1,            ! Turn on restraints
  restraintmask=":1-579 & !@H=", ! atoms to be restrained
  restraint_wt=0.1, ! force constant for restraint
  ntxo=1,           ! Write coordinate file in ASCII format
  ioutfm=0,         ! Write trajectory file in ASCII format
 /

 vim 05.min.mdin
 
   Minmize all the hydrogens
 &cntrl
  imin=1,           ! Minimize the initial structure
  maxcyc=1000,    ! Maximum number of cycles for minimization
  ntb=1,            ! Constant volume
  ntp=0,            ! No pressure scaling
  ntf=1,            ! Complete force evaluation
  ntwx= 1000,       ! Write to trajectory file every ntwx steps
  ntpr= 1000,       ! Print to mdout every ntpr steps
  ntwr= 1000,       ! Write a restart file every ntwr steps
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
  ntr=1,            ! Turn on restraints
  restraintmask=":1-579 & !@H=", ! atoms to be restrained
  restraint_wt=0.05, ! force constant for restraint
  ntxo=1,           ! Write coordinate file in ASCII format
  ioutfm=0,         ! Write trajectory file in ASCII format
 /

 vim 06.equil.mdin

  MD simualation
 &cntrl
  imin=0,           ! Perform MD
  nstlim=50000      ! Number of MD steps
  ntb=2,            ! Constant Pressure
  ntc=1,            ! No SHAKE on bonds between hydrogens
  dt=0.001,         ! Timestep (ps)
  ntp=1,            ! Isotropic pressure scaling
  barostat=1        ! Berendsen
  taup=0.5          ! Pressure relaxtion time (ps)
  ntf=1,            ! Complete force evaluation
  ntt=3,            ! Langevin thermostat
  gamma_ln=2.0      ! Collision Frequency for thermostat
  ig=-1,            ! Random seed for thermostat
  temp0=298.15      ! Simulation temperature (K)
  ntwx= 1000,       ! Write to trajectory file every ntwx steps
  ntpr= 1000,       ! Print to mdout every ntpr steps
  ntwr= 1000,       ! Write a restart file every ntwr steps
  cut=  8.0,        ! Nonbonded cutoff in Angstroms
  ntr=1,            ! Turn on restraints
  restraintmask=":1-579 & !@H=", ! atoms to be restrained
  restraint_wt=1.0, ! force constant for restraint
  ntxo=1,           ! Write coordinate file in ASCII format
  ioutfm=0,         ! Write trajectory file in ASCII format
  iwrap=1,          ! iwrap is turned on
 /

 vim 07.equil.mdin

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

 vim 08.equil.mdin

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

 vim 09.equil.mdin
 

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

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

 #!/bin/sh                                                                                                                                                                                                                                     
 #PBS -N 2BXF_equilibration                                                                                                                                                                                                                    
 #PBS -l walltime=04:00:00                                                                                                                                                                                                                     
 #PBS -l nodes=2:ppn=28                                                                                                                                                                                                                        
 #PBS -j oe                                                                                                                                                                                                                                    
 #PBS -q long                                                                                                                                                                                                                                  
 
 cd $PBS_O_WORKDIR
 
 echo "Started Equilibration on `date` "
 do_parallel="sander" 
 
 prmtop="../001.tleap_build/2BXF.wet.complex.prmtop"
 coords="../001.tleap_build/2BXF.wet.complex" 
 
 
 MDINPUTS=(01.min 02.equil 03.min 04.min 05.min 06.equil 07.equil 08.equil 09.equil)
 
 for input in ${MDINPUTS[@]}; do
 
  $do_parallel -O -i ${input}.mdin -o ${input}.mdout -p $prmtop -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf 
 ${input}.info -r ${input}.rst7
  coords=$input
 done 
 
 echo "Finished Equilibration on `date` "

Submit the job

 qsub md.equilibration.sh

= Production =

Move into 003.production

Move into 001.restrained

Create 10.prod.mdin

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

Create prod.sh
  
  #!/bin/sh                                                                                                                                                                                                                                     
 #PBS -N 2BXF_production                                                                                                                                                                                                                       
 #PBS -l walltime=24:00:00                                                                                                                                                                                                                     
 #PBS -l nodes=2:ppn=24                                                                                                                                                                                                                        
 #PBS -j oe                                                                                                                                                                                                                                                                                                                                                                                                                                                                      
 #PBS -o production.output
 #PBS -q long                                                                                                                                                                                                                     
 
 cd $PBS_O_WORKDIR
 
 echo "Started Production on `date` "
 #do_parallel="mpirun pmemd.MPI"                                                                                                                                                                                                               
 do_cuda="mpirun pmemd.MPI"
 
 prmtop="../../001.tleap_build/2BXF.wet.complex.prmtop"
 coords="../../002.equilbration/09.equil"
 
 
 MDINPUTS=(10.prod)
 
 for input in ${MDINPUTS[@]}; do
 
  $do_cuda -O -i ${input}.mdin -o ${input}.mdout -p $prmtop -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf 
 ${input}.info -r ${input}.rst7
  coords=$input
 done
 
 echo "Finished Production on `date` "


Submit the job
 qsub prod.sh

Move into 002.unrestrained

Create 10.prod.mdin

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


Create prod.sh
  
  #!/bin/sh                                                                                                                                                                                                                                     
 #PBS -N 2BXF_production                                                                                                                                                                                                                       
 #PBS -l walltime=24:00:00                                                                                                                                                                                                                     
 #PBS -l nodes=2:ppn=24                                                                                                                                                                                                                        
 #PBS -j oe                                                                                                                                                                                                                                                                                                                                                                                                                                                                      
 #PBS -o production.output
 #PBS -q long                                                                                                                                                                                                                     
 
 cd $PBS_O_WORKDIR
 
 echo "Started Production on `date` "
 #do_parallel="mpirun pmemd.MPI"                                                                                                                                                                                                               
 do_cuda="mpirun pmemd.MPI"
 
 prmtop="../../001.tleap_build/2BXF.wet.complex.prmtop"
 coords="../../002.equilbration/09.equil"
 
 
 MDINPUTS=(10.prod)
 
 for input in ${MDINPUTS[@]}; do
 
  $do_cuda -O -i ${input}.mdin -o ${input}.mdout -p $prmtop -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf 
 ${input}.info -r ${input}.rst7
  coords=$input
 done
 
 echo "Finished Production on `date` "


Submit the job
 qsub prod.sh

= Analysis =

== RMSD ==
How far our ligand moved and how far our receptor moved

Create cpptraj.strip.wat.in

 #!/usr/bin/sh                                                               
 parm ../../001.tleap_build/2BXF.wet.complex.prmtop
 #read in trajectory                                                         
 trajin ../../003.production/001.restrained/10.prod.trj
 #read in reference                                                          
 reference ../../001.tleap_build/2BXF.wet.complex.rst7
 #compute rmsd and align CA to the crystal structure                         
 rmsd rms1 reference :1-131@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 ":132&!(@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-131&!(@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-288 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.