#!/bin/sh
#set -xv

#  Automate all the steps needed to take a system through af-nmr.  

#------------------------------------------------------------------------------
#  set up usage statement:
#------------------------------------------------------------------------------
usage(){
cat<<EOD

Usage: `basename $0` [flags] basename

   where the file "basename.pdb" must be present

Optional flags:
   -lastres        final protein or na residue; things after that should
                   be ligands or water; default is 0, which implies that
                   there are no additional residues
   -mixedb         use a mixed basis with the primary residue at TZP,
                   others at DZP (default)
   -nomixedb       all quantum residues are treated at the TZP level
   -solinprot      use the solinprot program to get surface charges (default)
   -pbsa           use the PBSA program to get surface charges
   -rism, -3drism  use 3D-RISM to get surface charges
   -demon, -deMon  Set up quantum input files to use with deMon (default)
   -orca           Set up quantum input files to use with Orca
   -gau            Set up quantum input files to use with Gaussian
   -qchem          Set up quantum input files to use with Q-Chem
   -workdir        Do everything inside of a working directory named <basename>
   -nobuild        Skip the tleap & txmin steps; uses <basename>.pqr
   -nomin          Skip the minimization steps, but still build from 
                       <basename>.pdb
   -frcmod <file>  An additional frcmod file to load in to tleap
   -offlib <file>  An additional library file to load in to tleap

===========================================================
EOD

}

check_amberhome() {
   if [ -z "$AMBERHOME" ]; then
      echo ""
      echo "Your AMBERHOME environment variable is not set!"
      exit 1
   fi
}

check_shiftshome() {
   if [ -z "$SHIFTSHOME" ]; then
      echo ""
      echo "Your SHIFTSHOME environment variable is not set!"
      exit 1
   fi
}

implicit_solvent="solinprot"
solinprot="T"
lastres=0
program="D"
mixedb="P"
workdir="no"
min="yes"
build="yes"

#------------------------------------------------------------------------------
#  Checking Arguments:
#------------------------------------------------------------------------------
if [ $# -lt 1 ]; then usage; exit 1;  fi

cat > leap.in << EOF
source oldff/leaprc.ff14SB
loadAmberParams frcmod.ionsjc_tip3p
EOF

while [ $# -gt 0 ]; do
    case "$1" in
        -lastres)     shift; lastres="$1";;
        -mixedb)      mixedb="P";;
        -nomixedb)    mixedb="F";;
        -demon)       program="D";;
        -deMon)       program="D";;
        -orca)        program="O";;
        -gau)         program="G";;
        -qchem)       program="Q";;
        -solinprot)   solinprot="T"; implicit_solvent="solinprot";;
        -pbsa)        solinprot="F"; implicit_solvent="pb";;
        -rism)        solinprot="F"; implicit_solvent="3drism";;
        -3drism)      solinprot="F"; implicit_solvent="3drism";;
        -workdir)     workdir="T";;
        -nobuild)     build="no";;
        -nomin)       min="no";;
        -frcmod)      shift; echo "loadAmberParams $1" >> leap.in;;
        -offlib)      shift; echo "loadOFF $1" >> leap.in;;
        -help)        usage; exit 0;;
        --help)       usage; exit 0;;
        -h)           usage; exit 0;;

        -*) echo "Error: unknown flag: $1"
            usage
            exit 1;;
        *) if [ $# -gt 1 ]; then
             echo "Too many arguments"
             usage
             exit 1
           fi
           break;;
    esac
    shift
done

if [ "$workdir" = "T" ]; then
   /bin/mkdir -p $1
   cd $1
   /bin/rm -f $1.pdb
   /bin/ln -s ../$1.pdb
   cp ../leap.in .
fi

#  Set the options to afnmr.x:
AFNMR_OPTIONS="$program $mixedb $solinprot $lastres"

#  Program to analyze h-bonds (simulaid or cpptraj or blank)
GETHB=

#  Should not need to edit anything below here; (except to perhaps change 
#  detailed parameters in minimization, pbsa or 3drism)

#############################################################################

check_amberhome
check_shiftshome

if [ $build = "yes" ]; then

echo "Running AFNMR on $1.pdb"

#============================================================================
#  Use tleap to create the Amber prmtop file:
#============================================================================
echo "Running tleap"
cat <<EOF >> leap.in
source oldff/leaprc.ff14SB
source leaprc.gaff
addAtomTypes {
        { "M1"  "Cu" "sp3" }
        { "Y1"  "O" "sp3" }
        { "Y2"  "N" "sp3" }
        { "Y3"  "S" "sp3" }
        { "Y4"  "N" "sp3" }
        { "Y5"  "S" "sp3" }
}
GY1 = loadmol2 ~/p53_project/AZURIN/CU1/NO_OPT/AZU1/GY1.mol2
HE1 = loadmol2 ~/p53_project/AZURIN/CU1/NO_OPT/AZU1/HE1.mol2
CM1 = loadmol2 ~/p53_project/AZURIN/CU1/NO_OPT/AZU1/CM1.mol2
HE2 = loadmol2 ~/p53_project/AZURIN/CU1/NO_OPT/AZU1/HE2.mol2
MT1 = loadmol2 ~/p53_project/AZURIN/CU1/NO_OPT/AZU1/MT1.mol2
CU1 = loadmol2 ~/p53_project/AZURIN/CU1/NO_OPT/AZU1/CU1.mol2
loadamberparams frcmod.ions1lm_126_tip3p
loadamberparams ~/p53_project/AZURIN/CU1/NO_OPT/AZU1/azu_mcpbpy.frcmod
mol = loadpdb $1.pdb
bond mol.3.SG mol.26.SG
bond mol.45.O mol.129.CU
bond mol.46.ND1 mol.129.CU
bond mol.112.SG mol.129.CU
bond mol.117.ND1 mol.129.CU
bond mol.121.SD mol.129.CU
bond mol.44.C mol.45.N
bond mol.45.C mol.46.N
bond mol.46.C mol.47.N
bond mol.111.C mol.112.N
bond mol.112.C mol.113.N
bond mol.116.C mol.117.N
bond mol.117.C mol.118.N
bond mol.120.C mol.121.N
bond mol.121.C mol.122.N
saveamberparm mol $1.parm7 $1.x

EOF
#ix = loadpdb $1.pdb
#saveamberparm x $1.parm7 $1.x
#quit
#EOF

tleap -f leap.in > tleap.out
if [ ! -s $1.parm7 ]; then 
   echo "Error in tleap: check tleap.out or leap.log "
   exit 1
fi
cat tleap.out > $1.afnmr.log
/bin/rm leap.in tleap.out

if [ $min = "yes" ]; then

ambpdb -p $1.parm7 < $1.x > $1.0.pdb  2> /dev/null

#============================================================================
#  Do a constrained minimization to regularize the local geometries:
#============================================================================
echo "Minimization"

cat <<EOF > mdin
  constrained minimization
 &cntrl
   imin=1, ntb=0, igb=1, cut=99.0, rgbmax=10.0, 
   ntmin=3, maxcyc=50, drms=0.001, ntpr=1, ntxo=1,
   restraint_wt=5.0, restraintmask=':*', ntr=1,
 /
 &lmod
   xmin_method='TNCG', xmin_verbosity=0,
 /
EOF

sander -O -i mdin -p $1.parm7 -c $1.x -ref $1.x -r $1.min.x -o $1.min.o
if [ "$?" -gt 0 ]; then 
   echo "Error in sander: check $1.min.o"
   exit 1
fi
cat $1.min.o >> $1.afnmr.log
/bin/rm -f $1.min.o mdinfo mdin

#============================================================================
#  Create pqr for further use, with original residue numbers
#============================================================================
#parmed.py -n -p $1.parm7 <<EOF >> $1.afnmr.log
#setOverwrite
#addPDB $1.pdb
#outparm $1.parm7
#EOF

ambpdb -p $1.parm7 -pqr -ext -c $1.min.x > $1.min.pqr 2> /dev/null
/bin/rm -f $1.pqr
ln -s $1.min.pqr $1.pqr

#============================================================================
#  Check how far the struture has moved from the initial conformation:
#============================================================================
backbone="@C,CA,N,P,O3',O5',C3',C4',C5'"  # should work for both proteins
                                          # and nucleic acids
cpptraj $1.parm7 <<EOF > cpptraj.out
reference $1.x
trajin $1.min.x
rms all_heavy reference !@H= out rms.dat
rms backbone reference $backbone out rms.dat
EOF

if [ "$?" -gt 0 ]; then 
   echo "Error in cpptraj: check cpptraj.out"
   exit 1
fi
cat cpptraj.out rms.dat >> $1.afnmr.log
/bin/rm cpptraj.out rms.dat

#  end of "if min=yes" block:
else
    echo "Skipping minimization"
    ambpdb -p $1.parm7 -pqr -bres -ext < $1.x > $1.pqr 2> /dev/null
fi

#  end of "if build=yes" block:
else
    /bin/rm -f leap.in
    echo "nobuild: Running AFNMR on $1.pqr"
fi

#============================================================================
#  write a script to get MS surface points:
#============================================================================
cat <<EOF2 > runms.solsurf
#! /usr/bin/env perl
#
#----------------------------------------------------------------
#   run the "ms" program from input pqr file,(s) to get solvent-
#   accesible surface area.  Results to STDOUT
#----------------------------------------------------------------
#

#  solvent accessible radius for solvent probe radius of 1.0:
open (MSRAD, ">ms.rad") || die "can't open ms.rad\n";
print MSRAD <<"EOF";
    1      2.85                      phosphorus;
    2      2.50                      oxygen
    4      2.55                      nitrogen
    7      2.70                      carbon, allatom
    8      2.80                      sulfur
   15      2.20                      hydrogen
   16      2.63                      potassium
EOF

close(MSRAD); 

foreach \$file (@ARGV) {

#   change a pdb or pqr file to atm format, using our standard
#     atom types
#
open (PQR, "\$file") || die "can't find input file \$file\n";
open (ATM, ">input.atm") || die "can't open input.atm\n";

\$type{"H"} = 15; \$type{"C"} = 7; \$type{"N"} = 4; \$type{"P"} = 1;
\$type{"O"} = 2; \$type{"S"} = 8; \$type{"F"} = 4; \$type{"K"} = 16;

while (<PQR>) { 
	next unless /^ATOM|^HETATM/;
	( \$atno, \$atname, \$alt, \$resname, \$resno, \$x, \$y, \$z ) =
	    unpack("x7 a4 x a4 a a3 x3 a3 x4 a8 a8 a8",\$_);

	\$atname =~ s/^\s+//;  # remove leading spaces
	\$atname =~ s/\s+\$//;  # remove trailing spaces

	if ( defined(\$type{substr(\$atname,0,1)}) ){
		\$r = \$type{substr(\$atname,0,1)};
	} else {
		print STDERR "Error in atom type\n";
        exit 1;
	}
	printf ATM "%10.3f%10.3f%10.3f%5d    2    1\n",\$x,\$y,\$z,\$r;
}

close(ATM);

system ("$SHIFTSHOME/bin/ms -atm input.atm -rad ms.rad -d 0.5 -rp 0.0 ");
unlink("input.atm");

open (CONTACT, "contact") || die "contact file not found\n";
open (SRF, ">surf.pos") || die "cannot open surf.pos\n";

#
#  convert a "contact" solvent-accessible surface file to
#  one where each line corresponds a given surface point
#
while (<CONTACT>) {
	@_ = split(' ', \$_); 
    printf SRF "%10.3f%10.3f%10.3f\n", \$_[4],\$_[5],\$_[6];
}
close(CONTACT);
unlink("contact");

}
unlink("ms.rad");
EOF2

chmod +x runms.solsurf

if [ "$implicit_solvent" = "pb" ]; then

#============================================================================
#  run PBSA to get surface charge distribution
#============================================================================
echo "MS"
./runms.solsurf $1.pqr >> $1.afnmr.log
/bin/rm runms.solsurf

echo "PBSA"
cat > min.in <<EOF
 test of pbsa
 &cntrl 
   ntx=1, imin=1, igb=10, inp=0,
 /   
 &pb
   npbverb=0, istrng=0, epsout=80.0, epsin=1.0, space=0.5,
   accept=0.001, sprob=1.6, radiopt=0, fillratio=4, eneopt=2,
   frcopt=2, isurfchg=1,
 /   
EOF

pbsa -O -p $1.parm7 -i min.in -o pbsa.out -c $1.min.x 
if [ "$?" -gt 0 ]; then 
   echo "Error in pbsa: check pbsa.out"
   exit 1
fi
cat pbsa.out >> $1.afnmr.log
/bin/mv srfchg.pos input.xyzv
/bin/rm pbsa.out min.in

echo "gridprune"
$SHIFTSHOME/bin/gridprune < $1.pqr > gridprune.out
if [ "$?" -gt 0 ]; then 
   echo "Error in gridprune: check gridprune.out"
   exit 1
fi
cat gridprune.out >> $1.afnmr.log
pair < check.dat >> $1.afnmr.log
/bin/rm input.xyzv surf.pos gridprune.out check.dat

elif [ "$implicit_solvent" = "3drism" ]; then

#============================================================================
# get rism charge distribution:
#============================================================================
echo "MS"
./runms.solsurf $1.pqr >> $1.afnmr.log
/bin/rm runms.solsurf

echo "RISM"
rism3d.snglpnt --pdb $1.pqr --prmtop $1.parm7 \
   --xvv ${SHIFTSHOME}/dat/spc-nacl-3.xvv \
   --chgdist chgdist --buffer 12. --grdspc 0.5,0.5,0.5 --verbose 1 \
   --volfmt xyzv > rism3d.out
if [ "$?" -gt 0 ]; then 
   echo "Error in 3drism: check 3drism.out"
   exit 1
fi
cat 3drism.out >> $1.afnmr.log
/bin/mv chgdist.1.xyzv input.xyzv
/bin/rm 3drism.out

echo "gridprune"
$SHIFTSHOME/bin/gridprune < $1.pqr > gridprune.out
if [ "$?" -gt 0 ]; then 
   echo "Error in gridprune: check gridprune.out"
   exit 1
fi
cat gridprune.out >> $1.afnmr.log
pair < check.dat >> $1.afnmr.log
/bin/rm input.xyzv surf.pos gridprune.out check.dat


elif [ "$implicit_solvent" = "solinprot" ]; then

cat <<EOF3 > runsolinprot
#!/bin/sh

echo "MS"
./runms.solsurf \$1.pqr
if [ "\$?" -gt 0 ]; then 
   echo "Error in ms: check afnmr.log"
   exit 1
fi

echo "solinprot"
cat <<EOF > \$1.ogm
ON_GEOM_CENT 101 2.0
ON_GEOM_CENT 101 1.0
ON_GEOM_CENT 101 0.25
EOF

grep -v MOD \$1.pqr | awk '/ATOM/{print \$6,\$7,\$8}' > \$1.fpt

$SHIFTSHOME/bin/solinprot -epsin1 1.0 -epsin2 4.0 -ProteinField -ReactionField \$1 \$1.prot 
if [ "\$?" -gt 0 ]; then 
   echo "Error in solinprot: check afnmr.log"
   exit 1
fi
/bin/rm \$1.ogm

paste \$1.fpt \$1.rf \$1.pf | awk '{print \$1,\$2,\$3,\$4+\$5}' > input.xyzv
/bin/rm \$1.fpt \$1.rf \$1.pf

echo "gridprune"
$SHIFTSHOME/bin/gridprune -solinprot < \$1.pqr 
if [ "\$?" -gt 0 ]; then 
   echo "Error in gridprune: check afnmr.log"
   exit 1
fi
pair < check.dat
/bin/rm input.xyzv surf.pos check.dat
exit 0

EOF3
chmod +x runsolinprot

fi
#  (End of if statement about implicit_solvent)

#============================================================================
#  Get hydrogen bonding information (optional):
#============================================================================

if [ "$GETHB" = "simulaid" ]; then

#============================================================================
#  Run the simulaid program to get Hbond info. needed for afnmr
#============================================================================
echo "simulaid"
cat <<EOF > simulaid.in
a
$1.pqr
p
n
WAT
n
b
h
n
y
1.90
100.0
y
q
EOF

simulaid < simulaid.in
awk '$2=="H" && substr($7,1,1)=="O"' $1.pqr.hbn > hbond.hbn
#/bin/rm $1.pqr.hbn

elif [ "$GETHB" = "cpptraj" ]; then

#============================================================================
#  use cpptraj to get hbonds from amide H to oxygens
#============================================================================

cat <<EOF > hbond.in
parm $1.parm7
trajin $1.min.x
hbond dist 3.3 avgout hbond.hbn printatomnum donormask @N acceptormask @O=
EOF

cpptraj < hbond.in > cpptraj.out
cat cpptraj.out >> $1.afnmr.log
/bin/rm hbond.in cpptraj.out

# next line could be improved to mimic simulaid output format(?)
/usr/bin/env perl -i -p -e "s/@/ /g; s/_/ /g;" hbond.hbn

fi

#============================================================================
#  Run the fragment program to create the quantum chemistry input files;
#============================================================================
echo "afnmr.x"
$SHIFTSHOME/bin/afnmr.x  $AFNMR_OPTIONS  $1 > afnmr.out
if [ "$?" -gt 0 ]; then 
   echo "Error in afnmr: check afnmr.out"
   exit 1
fi
cat afnmr.out >> $1.afnmr.log
/bin/rm afnmr.out

#============================================================================
#  Some fine file clean-ups  (optional)
#============================================================================
/bin/rm -f $1.0.pdb $1.min.x $1???.prot.pqr \
       srfchg.pos runms.solsurf runsolinprot leap.log

if [ "$build" = "yes" ]; then
	/bin/rm -f $1.x $1.parm7
fi

if [ "$workdir" = "T" ]; then
   cd ../
   /bin/rm -f leap.in
fi

echo "done"