#!/bin/bash

# echo $1

#============================================================================
#  Use tleap to create the Amber prmtop file:
#============================================================================
cat <<EOF > leap.in
set default nocenter on
set default reorder_residues off
set default PBradii bondi
source leaprc.gaff2
loadAmberParams $1.frcmod
x = loadMol2 $1.mol2
saveamberparm x $1.parm7 $1.x
quit
EOF

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

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

../txmin $1.parm7 $1.x $1.min.x > $1.min.o

if [ "$?" -gt 0 ]; then 
   echo "Error in txmin: check $1.min.o"
   exit 1
fi
echo "====== min.o =========================================================" \
     >> $1.geom.log
cat $1.min.o >> $1.geom.log
/bin/rm -f $1.min.o

#  TODO: only make pdb files if there is a large RMS shift
#  ambpdb -p $1.parm7 -c $1.min.x > $1.min.pdb


#============================================================================
#  Check how far the struture has moved from the initial conformation:
#============================================================================

cpptraj $1.parm7 <<EOF > cpptraj.out
reference $1.x
trajin $1.min.x
rms all_heavy reference !@H= out rms.dat
EOF

if [ "$?" -gt 0 ]; then 
   echo "Error in cpptraj for $1: check cpptraj.out"
   exit 1
fi
echo "====== rms.dat  ======================================================" \
     >> $1.geom.log

if [ -s rms.dat ]; then
   cat rms.dat >> $1.geom.log
else
   echo "rms.dat missing for $1"
fi


/bin/rm -f cpptraj.out rms.dat $1.parm7 $1.min.x