#!/bin/bash

usage () {
	echo 
	echo "Script file made for GROMACS 2020.4 version compatibility."
	echo "This bash file must be in the desired working directory."
	echo "Usage:"
	echo "   -c [string]     convergence time: time in ps density converges (overestimate)"
	echo
	exit 1
	}

while getopts ":c:h" options; do
	case ${options} in
		c)
			CONVERGE=$OPTARG ;;
		h)
			usage ;;
		\?)
			echo "Invalid input. Unknown flag. Refer to the help with -h"
			exit 1 ;;
		:)
			echo "Invalid input: -$OPTARG requires an argument. Refer to the help with -h"
			exit 1 ;;
	esac
done
	
if [ $OPTIND -eq 1 ]; then
	echo "please specify convergence time with -c"
	exit 1 
fi

#fetch variables from variables.txt
#if variables.txt was not created by another program it must be created and placed in the workdir
#copy and edit the file from the shell scripts folder

while read -r flag value ;do
	case $flag in
		b)
			BASENAME=$value ;;
		m)
			MODEL=$value ;;
		t)
			TEMP=$value ;;
		f)
			FORCEFIELD=$value ;;
		n)
			TOTALNMOL=$value ;;
		s)
			STOCKMIX=$value ;;
		1)
			MEGNMOL=$value ;;
		2)
			DEGNMOL=$value ;;
		3)
			REGNMOL=$value ;;
		4)
			TEGNMOL=$value ;;
		5)
			PEGNMOL=$value ;;
		6)
			XEGNMOL=$value ;;
		7)
			HEGNMOL=$value ;;
		8)
			OEGNMOL=$value ;;
		9)
			NEGNMOL=$value ;;
		p)
			POLYEGNMOL=$value ;;
		w)
			WATERNMOL=$value ;;
	esac
done < variables.txt

VALUES=("$BASENAME" "$MODEL" "$TEMP" "$FORCEFIELD" "$TOTALNMOL" "$STOCKMIX" "$MEGNMOL" "$DEGNMOL" "$REGNMOL" "$TEGNMOL" "$PEGNMOL" "$XEGNMOL" "$HEGNMOL" "$OEGNMOL" "$NEGNMOL" "$POLYEGNMOL" "$WATERNMOL")
VARIABLES=("BASENAME" "MODEL" "TEMP" "FORCEFIELD" "TOTALNMOL" "STOCKMIX" "MEGNMOL" "DEGNMOL" "REGNMOL" "TEGNMOL" "PEGNMOL" "XEGNMOL" "HEGNMOL" "OEGNMOL" "NEGNMOL" "POLYEGNMOL" "WATERNMOL")
#check that all requested variables are defined
if [ -e variables.txt ]; then
	echo "variables.txt exists."
	index=0
	while [ $index -lt 17 ]; do 
		echo " ${VARIABLES[$index]} = ${VALUES[$index]}"
		if [ "${VALUES[$index]}" = "" ] ; then
			echo "${VARIABLES[$index]} was not defined"
			exit 1
		fi
	((index++))
	done
else
	echo "variables.txt does not exist" 
	exit 1
fi

WORKDIR="$(dirname "$(readlink -f "$0")")"
EM="$WORKDIR"/1_EM
EQUIL="$WORKDIR"/2_equil_NPT_${TEMP}K
PROD="$WORKDIR"/3_prod_NVT_${TEMP}K

#lists for indexed loops 
POLYA=("M" "D" "R" "T" "P" "X" "H" "O" "N")
POLYNUM=("$MEGNMOL" "$DEGNMOL" "$REGNMOL" "$TEGNMOL" "$PEGNMOL" "$XEGNMOL" "$HEGNMOL" "$OEGNMOL" "$NEGNMOL") 
POLYS=("MEGNMOL" "DEGNMOL" "REGNMOL" "TEGNMOL" "PEGNMOL" "XEGNMOL" "HEGNMOL" "OEGNMOL" "NEGNMOL")
POLYNAME=("MEG" "DEG" "REG" "TEG" "PEG" "XEG" "HEG" "OEG" "NEG")

#get average density
gmx analyze -f "$EQUIL"/analysis/density.xvg -b $CONVERGE >> "$EQUIL"/analysis/density.txt

#fluctuation properties
(echo "Volume" && echo "Temperature" && echo "Enthalpy") | gmx energy\
	-f "$EQUIL"/out/energy.edr\
	-o "$EQUIL"/analysis/fluct.xvg\
	-nmol $TOTALNMOL\
	-fluct_props -driftcorr \
	-b $CONVERGE >> "$EQUIL"/analysis/fluctuation.txt

#dielectric constant
echo "0" | gmx dipoles\
	-f "$EQUIL"/out/traj.xtc\
	-s "$EQUIL"/npt.tpr\
	-b $CONVERGE\
        -temp $TEMP\
	-o "$EQUIL"/analysis/Mtot.xvg\
	-d "$EQUIL"/analysis/dipsist.xvg\
	-eps "$EQUIL"/analysis/epsilon.xvg\
	-a "$EQUIL"/analysis/aver.xvg >> "$EQUIL"/analysis/epsilon.txt

#dipole moment
mkdir $EQUIL/analysis/dipole 

echo "q" | gmx make_ndx \
	-f $EQUIL/npt.tpr\
	-o $EQUIL/analysis/dipole/dipole.ndx

dip_group=2 
index=0
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then 
		echo "$dip_group" | gmx dipoles\
			-f $EQUIL/out/traj.xtc\
			-s $EQUIL/npt.tpr\
			-n $EQUIL/analysis/dipole/dipole.ndx\
			-b $CONVERGE\
			-temp $TEMP\
			-eps $EQUIL/analysis/dipole/epsilon_"${POLYNAME[$index]}".xvg\
			-o $EQUIL/analysis/dipole/Mtot_"${POLYNAME[$index]}".xvg\
			-d $EQUIL/analysis/dipole/dipdist_"${POLYNAME[$index]}".xvg\
			-a $EQUIL/analysis/dipole/aver_"${POLYNAME[$index]}".xvg
		echo "${POLYNAME[$index]}" >> $EQUIL/analysis/dipole.txt
		gmx analyze -f $EQUIL/analysis/dipole/epsilon_"${POLYNAME[$index]}".xvg >> $EQUIL/analysis/dipole.txt
		((dip_group++))
	fi
((index++))
done

if [ $WATERNMOL != 0 ]; then
	echo "$dip_group" | gmx dipoles\
		-f $EQUIL/out/traj.xtc\
		-s $EQUIL/npt.tpr\
		-n $EQUIL/analysis/dipole/dipole.ndx\
		-b $CONVERGE\
		-temp $TEMP\
		-eps $EQUIL/analysis/dipole/epsilon_water.xvg\
		-d $EQUIL/analysis/dipole/dipdist_water.xvg\
		-o $EQUIL/analysis/dipole/Mtot_water.xvg\
		-a $EQUIL/analysis/dipole/aver_water.xvg
	echo "water" >> $EQUIL/analysis/dipole.txt
	gmx analyze -f $EQUIL/analysis/dipole/epsilon_water.xvg >> $EQUIL/analysis/dipole.txt
	((dip_group++))
fi

