Revision as of 15:53, 23 March 2018

2nnq with an explicit solvent model

Prepare the files

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

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

Copy into zzz.master

Parameters

Move into 000.programs

Paramaterize the ligand

antechamber -i ../zzz.master/2nnq.lig.withH.charged.pdb -fi pdb -o 2nnq_lig.am1bcc.mol2 -fo mol2 -at gaff2 -c bcc -rn LIG -nc 1

Check for missing force field parameters

parmchk2 -i 2nnq_lig.am1bcc.mol2 -f mol2 -o 2nnq_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/2nnq.rec.withH.charged.pdb

###Load Ligand frcmod/mol2                                                      
loadamberparams ../000.parameters/2nnq_lig.am1bcc.frcmod
lig=loadmol2 ../000.parameters/2nnq_lig.am1bcc.mol2

###Create gas-phase complex                                                     
gascomplex= combine {rec lig}

###Write gas-phase pdb                                                          
savepdb gascomplex 2nnq.gas.complex.pdb

###Write gas-phase toplogy and coord files for MMGBSA calc                      
saveamberparm gascomplex 2nnq.gas.complex.prmtop 2nnq.gas.complex.rst7
saveamberparm rec 2nnq.gas.receptor.prmtop 2nnq.gas.receptor.rst7
saveamberparm lig 2nnq.gas.ligand.prmtop 2nnq.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

###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            

Create Amber topology and coordinates files for the MD simulation

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
 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-131 & !@H=", ! atoms to be restrained
 restraint_wt=5.0, ! force constant for restraint
 ntxo=1,           ! Write coordinate file in ASCII format
 ioutfm=0,         ! Write trajectory file in ASCII format

/

Create a submission script to equilibrate the biological system

#!/bin/sh                                                                                                                                                                                                                                     
#PBS -N 2nnq_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/2nnq.wet.complex.prmtop"
coords="../001.tleap_build/2nnq.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