Tutorial on running Molecular Dynamics for GIST grid generation with scripts 2

Jump to navigation Jump to search

here is a script (003md.tleap_reduce.csh) for preparing systems of standard protein residues (without co-factors):

# This script uses first reduce and then tleap to prepare a receptor for amber.
# The outputs are: parameter topology file (prm7) and a coordinate file (rst7) 
# Written by Trent balius
# TEB/ MF comments Feb2017

setenv AMBERHOME /nfs/soft/amber/amber14
setenv DOCKBASE "/nfs/home/tbalius/zzz.github/DOCK"
# CUDA for GPU
#setenv LD_LIBRARY_PATH ""
#setenv LD_LIBRARY_PATH "/usr/local/cuda-6.0/lib64/:$LD_LIBRARY_PATH"

set mountdir = `pwd`
set filedir = $mountdir/MDrundir/prep/002md_align
set workdir = $mountdir/MDrundir/prep/003md_tleap

# check if workdir exists
if ( -s $workdir ) then
   echo "$workdir exits"
# if it exists do this:     rm -rf ${workdir}
mkdir -p ${workdir}
cd ${workdir}

#rm -f leap.log  tleap.rec.out tleap.rec.in

## uncomment the next line to include Xtal waters into simulation
#cp $filedir/nearby_waters_aligned.pdb water.pdb
cp $filedir/apo_ref.pdb rec.pdb

# 1st: run leap to renumber residues.

# produces tleap input file -- renumbers rec
cat << EOF >! tleap.rec.1.in 
set default PBradii mbondi2
# load the protein force field
source leaprc.ff14SB
REC = loadpdb rec.pdb
saveamberparm REC rec.1.leap.prm7 rec.1.leap.rst7

# runs tleap and converts parameter and coordinate restart file back into pdb file
$AMBERHOME/bin/tleap -s -f tleap.rec.1.in > ! tleap.rec.1.out
$AMBERHOME/bin/ambpdb -p rec.1.leap.prm7 < rec.1.leap.rst7 >! rec.1.leap.ori.pdb 

# Remove hydrogens before running leap.
grep -v ' H$' rec.1.leap.ori.pdb >! rec.1.leap.pdb

# nomenclature clean-up
$DOCKBASE/proteins/Reduce/reduce -HIS -FLIPs rec.1.leap.pdb >! rec.nowat.reduce.pdb
sed -i 's/HETATM/ATOM  /g' rec.nowat.reduce.pdb 
# grep -v means "everything but"; last grep statement removes inconsistently named HEM hydrogens (for leap to be added back in)
grep "^ATOM  " rec.nowat.reduce.pdb | sed -e 's/   new//g' | sed 's/   flip//g' | sed 's/   std//g' | grep -v "OXT" | grep -v " 0......HEM" >! rec.nowat.reduce_clean.pdb
# grep -v means "everything but"; last grep statement removes inconsistently named HEM hydrogens (for leap to be added back in)

curl docking.org/~tbalius/code/waterpaper2017/scripts/replace_his_with_hie_hid_hip.py > replace_his_with_hie_hid_hip.py
curl docking.org/~tbalius/code/waterpaper2017/scripts/replace_cys_to_cyx.py > replace_cys_to_cyx.py
curl docking.org/~tbalius/code/waterpaper2017/scripts/add.ters.py > add.ters.py

# python scripts do these three things: 1) checks his to give protonation specific names; 2) checks for disulphide bonds; 3) checks for missing residues and adds TER flag
python replace_his_with_hie_hid_hip.py rec.nowat.reduce_clean.pdb rec.nowat.1his.pdb
python replace_cys_to_cyx.py rec.nowat.1his.pdb rec.nowat.2cys.pdb
python add.ters.py rec.nowat.2cys.pdb rec.nowat.3ter.pdb

# add TER to file rec.nowat.H_mod2.pdb before HEM
#grep -v "HEM" rec.nowat.H_mod3.pdb >! rec.nowat.H_mod4.pdb
#echo "TER" >> rec.nowat.H_mod4.pdb
#grep "HEM" rec.nowat.H_mod2.pdb >> rec.nowat.H_mod4.pdb
mv rec.nowat.3ter.pdb rec.nowat.final.pdb
#mv rec.nowat.H_mod3.pdb rec.nowat.H_final.pdb

# produces tleap input file in 3 steps, first script based instruction file (loading FF and rec file), then pipes disulphide info into middle and final instructions (write just rec, solvated in TIP3P water box of 10A about rec boundaries, write rec in water-box parameters and coordinates) for complete tleap input file
cat << EOF >! tleap.rec.in 
set default PBradii mbondi2
# load the protein force field
source leaprc.ff14SB
REC = loadpdb rec.nowat.final.pdb

if (-e rec.nowat.2cys.pdb.for.leap ) then
  cat rec.nowat.2cys.pdb.for.leap >> tleap.rec.in

cat << EOF >> tleap.rec.in
saveamberparm REC rec.leap.prm7 rec.leap.rst7
solvateBox REC TIP3PBOX 10.0
saveamberparm REC rec.watbox.leap.prm7 rec.watbox.leap.rst7

$AMBERHOME/bin/tleap -s -f tleap.rec.in > ! tleap.rec.out

# for ease of visualization in pymol 
$AMBERHOME/bin/ambpdb -p rec.leap.prm7 < rec.leap.rst7 > rec.leap.pdb

echo "Look at rec.leap.pdb in pymol. May have to delete last column in pdb file for pymol. "
# may have to remove element column so pymol does not get confused

# inspect tleap.rec.out and the leap.log file
# visually inspect (VMD) rec.10wat.leap.prm7, rec.10wat.leap.rst7 and rec.watbox.leap.prm7, rec.watbox.leap.rst7

This script can easily be modified to work with metal ions. Also it may be useful to add in capping groups if breaks in the protein (do to missing residues) are near the binding site (site of interest).