Some scripts for Amber

0-AmberScripts.md

This gist contains some scripts for Amber.

• anteligand.sh: Use antechamber to prepare a ligand for Amber
• complex_tleap.sh: Use tleap to prepare a ligand-protein complex input parameter file (Amber14)
• checkPDBID_MDprep.sh: After running many preparation jobs for PDB ligand-protein complexes, you should check whether the preparation is done.
• rename_water.py: rename water to HOH O,H1,H2 for PDB/PQR

Maybe need some other tools installed, such as pdb2pqr and g09 and so on.

anteligand.sh

#! /bin/bash
# Author: Hom. 2015.11.12
# Usage: $0 inputfile chargetype
# For: Convert normal molecule file to frcmod and charged input file

if [ -z $1 ]; then
echo "No ligand file is assigned!"
exit 1
fi


if [ -z ${AMBERHOME} ]; then
#echo "The AMBERHOME var or ligand mol2 file is not set! Check it!"
#echo "Use it as './antchm.sh ligand.mol2'"
#exit

# If no define AMBERHOME, source the ambertools
source ~/AmberTools/amber.sh
fi


chargemode=$2
if [ -z $2 ];then
    chargemode="bcc"
fi

basename=${1%.*}
exdname=${1##*.}

if [ $exdname = "pqr" ];then
    cp $1 ${basename}.mpdb
    exdname="mpdb"
fi

ligbcc=${basename}_${chargemode}
if [ $chargemode = "no" ];then
    antechamber -i ${1} -fi $exdname -o ${ligbcc}.mol2 -fo mol2 -pf y
else
    antechamber -i ${1} -fi $exdname -o ${ligbcc}.mol2 -fo mol2 -c $chargemode -pf y
fi

sleep 1

antechamber -i ${ligbcc}.mol2 -fi mol2 -o ${ligbcc}.prepi -fo prepi -pf y
parmchk -i ${ligbcc}.prepi -f prepi -o ${ligbcc}.frcmod

if [ $exdname = "pqr" ];then
    rm ${basename}.mpdb
fi

rm sqm.*

checkPDBID_MDprep.sh

#! /bin/bash
# After running many preparation jobs for PDB ligand-protein complexes, you should check whether the preparation is done.
# This script help to check the input files for MD, also known as output from preparation.
#
# Need Argument: a file saving pdb id
# Cases are run in PDB id directories and save lig.top/crd, pro.top/crd, com.top/crd and complex.top/crd in each directory.
# Output to kout.log

echo "Checking wrong results...." | tee kout.log
 
for dir in `cat $1`
do
cd $dir
 
if [ ! -f lig.top -o ! -f lig.crd ];then
    echo "$dir ligand not finished"  >> ../kout.log
fi
 
if [ ! -f pro.top -o ! -f pro.crd ];then
    echo "$dir protein not finished"  >> ../kout.log
fi
 
if [ ! -f com.top -o ! -f com.crd ];then
    echo "$dir complex not finished"  >> ../kout.log
fi
 
if [ ! -f complex.top -o ! -f complex.crd ];then
    echo "$dir complex out not finished"  >> ../kout.log
fi
 
cd ..
done
 
echo "">> kout.log
 
echo "Wrong ligand: " `grep "ligand" kout.log |wc -l` > kout_tmp.tmp
echo "Wrong protein: " `grep "protein" kout.log |wc -l` >> kout_tmp.tmp
echo "Wrong complex: " `grep "complex out" kout.log |wc -l` >> kout_tmp.tmp
 
echo "#! /usr/bin/env python
import os,sys
fi=open(sys.argv[1])
wrongid=[]
for line in fi:
    tmp=line.split()
    if (len(tmp)>0 and len(tmp[0])==4):
        wrongid.append(tmp[0])
setid=set(wrongid)
print 'Total wrong cases: ',len(setid)
for i in setid:
    print i" > tmp_result.py
python tmp_result.py kout.log >> kout_tmp.tmp
cat kout_tmp.tmp >> kout.log
rm kout_tmp.tmp tmp_result.py

complex_tleap.sh

#!/bin/bash
# This script need 3 arguments:
#   1. ligand input file
#   2. protein input file
#   3. box radius

echo "#############################################################################"
echo "The script help to build a leap.in file and run the tleap for the complex."
echo "First parameter is the ligand file name ABC"
echo "Second is the protein file name"
echo "Third is the radius of TIP3P box."
# echo -n "Enter anything to continue....."
# read
echo "#############################################################################"

if [ -z ${AMBERHOME} ]; then
# echo "ERROR:The AMBERHOME var is not set! Check it!"
# echo "ERROR:Be sure the tleap could be run!"
# exit

# If no define AMBERHOME, source the ambertools
source ~/AmberTools/amber.sh
fi

ligfile=$1
if [ -z $1 ];then
ligfile="ligand.mol2"
echo "Set ligand file name to default: ligand.mol2"
fi

profile=$2
if [ -z $2 ];then
profile="protein.pdb"
echo "Set protein file name to default: protein.pdb"
fi

boxradius=$3
if [ -z $3 ]; then
boxradius="10"
echo "Set box radius to default: 10"
fi

basename1=${ligfile%.*}
exdname1=`echo "${ligfile##*.}" | tr A-Z a-z`

if [ ! -f ${basename1}.frcmod ];then
echo "No ligand parameters file (ligandfile_prefix.frcmod)!"
exit 1
fi

# other options:
#   loadamberprep ligand.prep
#    list

echo "source leaprc.ff14SB
source leaprc.gaff
loadamberparams ${basename1}.frcmod">leap.in

if [ exdname1 = "mol2" ];then
echo "lig=loadmol2 $ligfile" >>leap.in
# no recommend with pdb file for ligand, no charge
elif [ exdname1 = "pdb" ];then
echo "lig=loadpdb $ligfile" >>leap.in
fi

echo "check lig
saveamberparm lig lig.top lig.crd
pro=loadpdb $profile
saveamberparm pro pro.top pro.crd
com=combine {lig pro}
saveamberparm com com.top com.crd
check com
addions com Na+ 0
addions com Cl- 0
solvatebox com TIP3PBOX $boxradius
check com
saveamberparm com complex.top complex.crd
quit" >>leap.in

tleap -s -f leap.in
ambpdb -p complex.top < complex.crd > complex_out.pdb

echo "###############################################################################"
echo "Note: Finish!"
echo "Note: Check the residues num for the restrain in minimization and heat step in MD!"

rename_water.py

#! /usr/bin/env python
import sys
fin=sys.argv[1]; 
fout=sys.argv[2]
fi=open(fin); 
fo=open(fout,'w')
op=0
for line in fi:
    if (line[:6]=='ATOM  ' or line[:6]=='HETATM'):
        if ( line[17:20]=='WAT' or line[17:20]=='HOH'):
            if (line[13]=='O'):
                op=0
                fo.write(line[:13]+'O   HOH'+line[20:])
            elif (line[13]=='H'):
                op+=1
                if op is 1:
                    fo.write(line[:13]+'H1  HOH'+line[20:])
                else:
                    fo.write(line[:13]+'H2  HOH'+line[20:])
                    op=0
            elif (line[12]=='H'):
                op+=1
                if op is 1:
                    fo.write(line[:12]+' H1  HOH'+line[20:])
                else:
                    fo.write(line[:12]+' H2  HOH'+line[20:])
                    op=0
            else: fo.write(line)
        else: fo.write(line)
    else:fo.write(line)
fi.close(); fo.close()