#!/bin/bash
shopt -s expand_aliases
gmx 1&>/dev/null || module load gromacs/2018.3-CUDA-8.0 2>/dev/null
#|| alias gmx="/nfsopt/gromacs/gromacs-5.1-sse41/jessie/bin/gmx"

usage() { echo "Usage: $0 -p directoryPath/ -e endfrac[0.5] -a automatic -s startfrac[0.05] -f forcefrac[-1] -d difference[0.00005] [-m dt modulo in ps (for -dt of gmx traj)] -r redo box.xvg; fractions mean how much of the trajectory is skipped, 4*startfrac is standard!, if automatic fails and -f is not used no gro will be cut" 1>&2; exit 1; }

path="$(pwd)"
FORCEFRAC=-1
difference=0.00005
redo=0
STARTFRAC=0.05
ENDFRAC=0.5
AUTOMATIC=""
IGNORE=0
MODULO=""
while getopts "hf:d:p:ras:e:im:" opt; do
  case $opt in
    f ) 
	    FORCEFRAC=$OPTARG
        ;;
    d ) 
        difference=$OPTARG
        ;;
    p ) 
        path="$OPTARG"
        ;;
    a ) 
        AUTOMATIC=" -a"
        ;;
    s ) 
        STARTFRAC=$OPTARG
        ;;
    e ) 
        ENDFRAC=$OPTARG
        ;;
    m ) 
        MODULO=" -dt $OPTARG"
        ;;
    r ) 
        redo=1
        ;;
    i ) 
        IGNORE=1
        ;;
    h ) 
	usage
	;;
   \? ) 
	usage
	;;
  esac
done

if [ "${path: -1}" != "/" ] ; then
	path="${path}/"
fi

echo "going to cut optimal box for $path"
#find box.xvg and if not make it
boxfile="${path}box.xvg"

if find "${path}"*/ -maxdepth 1 -name "*.gro" | grep -q .; then
    OUTGRO="$(find "${path}"*/ -maxdepth 1 -name "*.gro")"
    OUTGRO="${OUTGRO[0]}"
    echo "$OUTGRO"
fi

if [ ! -f "$OUTGRO" ] && [ "$IGNORE" == "0" ] ; then
    echo "simulation seems to not have finished, no */out.gro found; exiting"
    echo "use -i if sure"
    exit 1
fi

if [ -f "$OUTGRO" ] && [ -f "$boxfile" ] && [ "$OUTGRO" -nt "$boxfile" ] ; then
    echo "box.xvg existed but out.gro was newer; redoing"
    redo=1
fi

if find "$path" -name *.xtc | grep -q .; then
    trajectory="$(find "$path" -name *.xtc)"
    trajectory="${trajectory[0]}"
else
    echo no trajectory found
    exit 1
fi
if find "$path" -name run*.tpr | grep -q .; then
    runfile="$(find "$path" -name run*.tpr)"
    runfile="${runfile[0]}"
elif find "$path" -name *.tpr | grep -q .; then
    runfile="$(find "$path" -name *.tpr)"
    runfile="${runfile[0]}"
else
    echo "no run*.tpr or *.tpr file found"
    exit 1
fi

if [ -f "$boxfile" ] && [ $redo == 0 ] ; then
	echo "box.xvg already existed"
else
	echo "box.xvg not found or redoing box.xvg"
	echo "using: yes 0 | gmx traj -s $runfile -f $trajectory -ob $boxfile -fp $MODULO"
	yes 0 | gmx traj -s "$runfile" -f "$trajectory" -ob "$boxfile" -fp $MODULO 2>/dev/null
fi
echo "----------------------"
echo "Determining optimal box for the trajectory"
FRAMETIME=($(python3 /home/mmh/BrockportMDFiles/python_scripts/boxsize.py -x "$boxfile" -f $FORCEFRAC -d $difference $AUTOMATIC -s $STARTFRAC -e $ENDFRAC))
echo $FRAMETIME

if [ "$FRAMETIME" != "" ]; then
    equilibriumfile="${path}"equilibrium.gro
    echo "Will cut $equilibriumfile from $trajectory at $frameTime using"
    echo "yes 0 | gmx trjconv -s $runfile -f $trajectory -o $equilibriumfile -b $FRAMETIME -e $FRAMETIME"
    yes 0 | gmx trjconv -s "$runfile" -f "$trajectory" -o "$equilibriumfile" -b $FRAMETIME -e $FRAMETIME 2>/dev/null
    echo "done"
    echo "$runfile" "$trajectory"
else
    echo "Automatic determination of cutfraction failed. change parameters, prolong simulation or force cut with -f!"
    exit 1
fi

exit 0

