Minimize protein-ligand complex with AMBER

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. 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=%6.2f\n\n",sum}'
set chrg = `cat lig.gast.mol2 | awk 'BEGIN{sum=0};{if(NF==9){sum=sum+$9}};END{printf"%6.2f\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"

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"

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