Minimize protein-ligand complex with AMBER: Difference between revisions

From DISI
Jump to navigation Jump to search
 
Line 125: Line 125:
  setenv AMBERHOME /nfs/soft/amber/amber14
  setenv AMBERHOME /nfs/soft/amber/amber14
  setenv LD_LIBRARY_PATH ""
  setenv LD_LIBRARY_PATH ""
  setenv LD_LIBRARY_PATH "/usr/local/cuda-6.0/lib64/:$LD_LIBRARY_PATH"
  #setenv LD_LIBRARY_PATH "/usr/local/cuda-6.0/lib64/:$LD_LIBRARY_PATH"
setenv LD_LIBRARY_PATH "/nfs/soft/cuda-6.5/lib64/:\$LD_LIBRARY_PATH"
   
   
  cat << EOF1 > ! 01mi.in
  cat << EOF1 > ! 01mi.in
Line 172: Line 173:
  setenv AMBERHOME /nfs/soft/amber/amber14
  setenv AMBERHOME /nfs/soft/amber/amber14
  setenv LD_LIBRARY_PATH ""
  setenv LD_LIBRARY_PATH ""
  setenv LD_LIBRARY_PATH "/usr/local/cuda-6.0/lib64/:$LD_LIBRARY_PATH"
  #setenv LD_LIBRARY_PATH "/usr/local/cuda-6.0/lib64/:$LD_LIBRARY_PATH"
setenv LD_LIBRARY_PATH "/nfs/soft/cuda-6.5/lib64/:\$LD_LIBRARY_PATH"  
   
   
  cat << EOF1 > ! 01mi.lig.in
  cat << EOF1 > ! 01mi.lig.in

Latest revision as of 01:34, 15 December 2018

Written on 2018/04/02

This tutorial is for Shoichet Lab group members but we hope that it might be useful to the community.

Outlined below are five steps to minimize a complex.

This is helpful for preparing for a docking calculation (or perhaps, for minimizing a docked pose).


Step 1. Break up the complex into ligand (small molecule) and receptor.

You can do this in a number of ways. You can use your favorite molecular visualization software like UCSF Chimera. Or, you can use your favorite text editor like VIM to do it.

Here is a script that does it automatically:

Create a script called "run.001.be_balsti_reduce.csh". Put the following lines into the script using your favorite text editor like VIM.

#!/bin/csh 

set path = ( /nfs/home/tbalius/zzz.programs/msms $path )

  ln -s /nfs/home/tbalius/zzz.programs/msms/atmtypenumbers .
  #python ~/zzz.scripts/be_blasti.py --pdbcode $pdbname nocarbohydrate renumber | tee -a pdbinfo_using_biopython.log
  python ~/zzz.scripts/be_blasti.py --pdbfile complex.pdb nocarbohydrate original_numbers | tee -a pdbinfo_using_biopython.log


  if !(-s rec.pdb) then
      echo "rec.pdb is not found"
  endif

  mv rec.pdb temp.pdb
  grep -v TER temp.pdb | grep -v END  > rec.pdb

  rm temp.pdb

  if (-s lig.pdb) then
     sed -e "s/HETATM/ATOM  /g" lig.pdb > xtal-lig.pdb
  else if (-s pep.pdb) then ## if no ligand and peptide
     sed -e "s/HETATM/ATOM  /g" pep.pdb > xtal-lig.pdb
  else
     echo "Warning: No ligand or peptid."
  endif

run the script as follows:

 csh run.001.be_balsti_reduce.csh

Step 2. Prepare ligand charges and force field.

Here is a script that calls antechamber to prepare ligand charges and force field. Here are three things to be mindful of:

  • The flag "-nc" specifies the net charge of the small molecule be sure to replace this with the correct value.
  • The ligand must have hydrogen's added.
  • Also, the atom typing is dependent on the 3D coordinates.

You can call openbabel to calculate charges, for example:

# add GASTEIGER charges 
obabel -ipdb lig.pdb -omol2 -O lig.gast.mol2 -charge
cat lig.gast.mol2 | awk 'BEGIN{sum=0};{if(NF==9){print $9;sum=sum+$9}};END{printf"\nsum=%.0f\n\n",sum}'
set chrg = `cat lig.gast.mol2 | awk 'BEGIN{sum=0};{if(NF==9){sum=sum+$9}};END{printf"%.0f\n",sum}'`


Create a script called "run.002.ligprep.antechamber.csh".

#! /bin/tcsh


set workdir = `pwd`
cd $workdir


 setenv AMBERHOME /nfs/soft/amber/amber14

rm lig; mkdir lig; cd lig

#cp $workdir/xtal-lig.pdb lig.pdb
cp $workdir/33443.pdb lig.pdb
#sed -i 's/<0> /LIG/g' lig1.mol2

$AMBERHOME/bin/antechamber -i lig.pdb -fi pdb -o lig.ante.mol2 -fo mol2 

$AMBERHOME/bin/antechamber -i lig.ante.mol2 -fi mol2 -o lig.ante.charge.mol2 -fo mol2 -c bcc -at sybyl -nc 1
$AMBERHOME/bin/antechamber -i lig.ante.mol2 -fi mol2  -o lig.ante.pdb  -fo pdb
$AMBERHOME/bin/antechamber -i lig.ante.charge.mol2 -fi mol2  -o lig.ante.charge.prep -fo prepi
$AMBERHOME/bin/parmchk -i lig.ante.charge.prep -f  prepi -o lig.ante.charge.frcmod

Step 3. Prepare the coordinate file and topology-and-parameter files for amber.

Create a script called "run.003.tleap.csh".

setenv AMBERHOME /nfs/soft/amber/amber14

#grep -v 'H$' receptor.pdb > rec.pdb

cat << EOF >! tleap.in
set default PBradii mbondi2
# load the protein force field
source leaprc.ff12SB
# load in GAFF
source leaprc.gaff

# load ligand and covalent parameters.  
loadamberparams lig/lig.ante.charge.frcmod
loadamberprep lig/lig.ante.charge.prep

# load pdb file 
REC = loadpdb rec.pdb
LIG = loadpdb lig/lig.ante.pdb 
COM  = combine {REC LIG}

saveamberparm REC rec.leap.prm7 rec.leap.rst7
saveamberparm LIG lig.leap.prm7 lig.leap.rst7
saveamberparm COM com.leap.prm7 com.leap.rst7

quit
EOF

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

Step 4. Run minimization using the amber program PMEMD.cuda on a GPU machine

Create a script called "run.004.pmemd_cuda_min.csh"

#setenv AMBERHOME /nfs/soft/amber/amber14

setenv AMBERHOME /nfs/soft/amber/amber14
setenv LD_LIBRARY_PATH ""
#setenv LD_LIBRARY_PATH "/usr/local/cuda-6.0/lib64/:$LD_LIBRARY_PATH"
setenv LD_LIBRARY_PATH "/nfs/soft/cuda-6.5/lib64/:\$LD_LIBRARY_PATH"

cat << EOF1 > ! 01mi.in
01mi.in: minimization with GB
&cntrl
 imin = 1, maxcyc = 10000, ncyc = 500,  ntmin = 1,
 igb=1,
 ntx = 1, ntc = 1, ntf = 1,
 ntb = 0, ntp = 0,
 ntwx = 1000, ntwe = 0, ntpr = 1000,
 cut = 999.9,
 ntr = 1,
 restraintmask = '!@H=', 
 restraint_wt = 0.1,
/
EOF1


#$AMBERHOME/bin/pmemd.cuda -O -i 01mi.in -o 01mi.out -p com.leap.prm7 -c com.leap.rst7 -ref com.leap.rst7 -x 01mi.mdcrd -inf 01mi.info -r 01mi.rst7
#$AMBERHOME/bin/sander -O -i 01mi.in -o 01mi.out -p com.leap.prm7 -c com.leap.rst7 -ref com.leap.rst7 -x 01mi.mdcrd -inf 01mi.info -r 01mi.rst7

set pwd = `pwd`
#cd $pwd
 
cat << EOF > ! qsub.sander.csh
#\$ -S /bin/csh
#\$ -cwd
#\$ -q gpu.q
#\$ -o stdout
#\$ -e stderr

  cd $pwd
  
  $AMBERHOME/bin/pmemd.cuda -O -i 01mi.in -o 01mi.out -p com.leap.prm7 -c com.leap.rst7 -ref com.leap.rst7 -x 01mi.mdcrd -inf 01mi.info -r 01mi.rst7

EOF

  qsub qsub.sander.csh 

You might want to also minimize the ligand on its own:

Create a script called "run.004.pmemd_cuda_min.lig.csh".

#setenv AMBERHOME /nfs/soft/amber/amber14

setenv AMBERHOME /nfs/soft/amber/amber14
setenv LD_LIBRARY_PATH ""
#setenv LD_LIBRARY_PATH "/usr/local/cuda-6.0/lib64/:$LD_LIBRARY_PATH"
setenv LD_LIBRARY_PATH "/nfs/soft/cuda-6.5/lib64/:\$LD_LIBRARY_PATH" 

cat << EOF1 > ! 01mi.lig.in
01mi.in: minimization with GB
&cntrl
 imin = 1, maxcyc = 10000, ncyc = 500,  ntmin = 1,
 igb=1,
 ntx = 1, ntc = 1, ntf = 1,
 ntb = 0, ntp = 0,
 ntwx = 1000, ntwe = 0, ntpr = 1000,
 cut = 999.9,
 ntr = 0,
/
EOF1
#restraintmask = '!@H=', 
#restraint_wt = 0.1,


#$AMBERHOME/bin/pmemd.cuda -O -i 01mi.in -o 01mi.out -p com.leap.prm7 -c com.leap.rst7 -ref com.leap.rst7 -x 01mi.mdcrd -inf 01mi.info -r 01mi.rst7
#$AMBERHOME/bin/sander -O -i 01mi.lig.in -o 01mi.lig.out -p lig.leap.prm7 -c lig.leap.rst7 -ref lig.leap.rst7 -x 01mi.lig.mdcrd -inf 01mi.lig.info -r 01mi.lig.rst7 

set pwd = `pwd`
#cd $pwd

cat << EOF > ! qsub.sander.lig.csh
#\$ -S /bin/csh
#\$ -cwd
#\$ -q gpu.q
#\$ -o stdout
#\$ -e stderr

  cd $pwd
  $AMBERHOME/bin/pmemd.cuda -O -i 01mi.lig.in -o 01mi.lig.out -p lig.leap.prm7 -c lig.leap.rst7 -ref lig.leap.rst7 -x 01mi.lig.mdcrd -inf 01mi.lig.info -r 01mi.lig.rst7
   
EOF

  qsub qsub.sander.lig.csh

Step 5. Convert to pdb for easy visualization.

cat run.005.mkpdb.csh | awk '{print " " $0}'

#setenv AMBERHOME /nfs/soft/amber/amber14

setenv AMBERHOME /nfs/soft/amber/amber14
setenv LD_LIBRARY_PATH ""
setenv LD_LIBRARY_PATH "/usr/local/cuda-6.0/lib64/:$LD_LIBRARY_PATH"

  
   $AMBERHOME/bin/ambpdb -p com.leap.prm7 < 01mi.rst7 > 01mi.pdb
#   $AMBERHOME/bin/ambpdb -p lig.leap.prm7 < 01mi.lig.rst7 > 01mi.lig.pdb