#!/bin/bash

# help and usage (called from command line with -h)
usage() {
	echo
	echo "Script file made for GROMACS 2020.4 version compatibility."
	echo "Individual simulation runs for determining self-diffusion in thermodynamic limit and shear viscosity via finite-sized effects."
	echo "See J. Chem. Theory Comput. 2018, 14, 11, 5959–5968"
	echo "This bash file must be in the desired working directory."
	echo
	echo "Start initial run:"
	echo "  - pass no arguments to command line (user input will be asked for)"
	echo "  - pass all flags to command line "
	echo
	echo "NOTE 1: Please choose of the above or the script will not run."
	echo "NOTE 2: If using '-m bulk', then use '-w 0' when specifying molecule numbers."
	echo "NOTE 3: If using '-b water', then use '-p 0' when specifying molecule numbers."
	echo "NOTE 4: You must use either -s or all integer flags."
	echo "NOTE 5: If -s is not used all integer flags must have a value. Use 0 for 0. Do not skip any integer flags or leave blank."
	echo
	echo "For all simulations make sure the following files are in the mixtures files directory"
	echo "  - all itp files for desired Poly EGs"
	echo "  - all single molecule gro files for desired Poly EGs"
	echo "  - [MODEL].itp (water itp file) (where [MODEL] is chosen with '-m spce' or '-m tip4p2005')"
	echo "  - single_[MODEL].gro"
	echo "  - top_mixtures.top"
	echo
	echo "and the following files in working directory:"
	echo "  - em.mdp"
	echo "  - eq_NPT_[TEMP]K.mdp"
	echo "  - md_NVT_[TEMP]K.mdp"
	echo "  - This shell file"
	echo 
	echo "Usage:"
	echo "  -b [string]     basename: identifying name for the simulation (use underscore instead of space if necessary e.g. PEG200_neat)"
	echo "  -m [string]     model to simulate (enter 'bulk' for bulk PolyEG, 'spce' for spc/e water, or 'tip4p2005' for tip4p/2005 water)"
	echo "  -t [value]      temperature in K (use underscore instead of decimal if necessary e.g. 298_15)"
	echo "  -f [string]     forcefield (enter 'OPLS' for unmodified OPLS-AA, or 'OPLSMOD' for modified OPLS-AA"
	echo "  -d [value]      density in g/cm^3 at desired temperature"
	echo "  -n [value]      total number of molecules in system"
	echo "  -s [string]     stock mixture. if this option is used integer flags should not be used. enter 'none' if not using a stock mixture. current options are: PEG200, none."
	echo "  -1 [value]      how many MEG (monoethylene glycol) molecules. enter 0 if none"
	echo "  -2 [value]      how many DEG (diethylene glycol) molecules. enter 0 if none"
	echo "  -3 [value]      how many REG (triethylene glycol) molecules. enter 0 if none"
	echo "  -4 [value]      how many TEG (tetraethylene glycol) molecules. enter 0 if none"
	echo "  -5 [value]      how many PEG (pentaethylene glycol) molecules. enter 0 if none"
	echo "  -6 [value]      how many XEG (hexaethylene glycol) molecules. enter 0 if none"
	echo "  -7 [value]      how many HEG (heptaethylene glycol) molecules. enter 0 if none"
	echo "  -8 [value]      how many OEG (octaethylene glycol) molecules. enter 0 if none"
	echo "  -9 [value]      how many NEG (nonaethylene glycol) molecules. enter 0 if none"
	echo "  -p [value]      total number of PolyEG (polyethylene glycol) molecules"
	echo "  -w [value]      how many molecules are water"
	echo    
	exit 1
}

# setting working directory and getting gromacs stuff.
WORKDIR="$(dirname "$(readlink -f "$0")")"
source /usr/local/gromacs/bin/GMXRC
export GMX_MAXCONSTRWARN=-1;   

##### Begin taking in parameters
# command line inputs
while getopts ":b:1:2:3:4:5:6:7:8:9:p:m:t:d:n:w:s:f:h" options; do
	case ${options} in
		b)
			BASENAME=$OPTARG ;;
		m)
			MODEL=$OPTARG ;;
		t)
			TEMP=$OPTARG ;;
		f)
			FORCEFIELD=$OPTARG ;;
		d)
			DENSITY=$OPTARG ;;
		n)
			TOTALNMOL=$OPTARG ;;
		s)
			STOCKMIX=$OPTARG ;;
		1)
			MEGNMOL=$OPTARG ;;
		2)
			DEGNMOL=$OPTARG ;;
		3)
			REGNMOL=$OPTARG ;;
		4)
			TEGNMOL=$OPTARG ;;
		5)
			PEGNMOL=$OPTARG ;;
		6)
			XEGNMOL=$OPTARG ;;
		7)
			HEGNMOL=$OPTARG ;;
		8)
			OEGNMOL=$OPTARG ;;
		9)
			NEGNMOL=$OPTARG ;;
		p)
			POLYEGNMOL=$OPTARG ;;
		w)
			WATERNMOL=$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

# user input if nothing is specified to command line
if [ $OPTIND -eq 1 ]; then
	read -p "Simulation name: " BASENAME
	read -p "Model to simulate (enter 'bulk' for bulk PolyEG, 'spce' for spc/e water, or 'tip4p2005' for tip4p/2005 water): " MODEL
	read -p "Temperature (in K) (Use underscore instead of decimal if necessary, e.g. 298_15): " TEMP
	read -p "Forcefield (enter 'OPLS' for unmodified OPLS-AA, or 'OPLSMOD' for modified OPLS-AA" FORCEFIELD
	read -p "${BASENAME}'s density at ${TEMP}K in g/cm^3 (if unknown, choose something reasonable): " DENSITY
	read -p "Number of total molecules in system: " TOTALNMOL
	read -p "Stock mixture (PEG200 or none):" STOCKMIX
	if [ $STOCKMIX = "none" ]; then 
    		read -p "How many monoethylene glycol molecules: " MEGNMOL
    		read -p "How many diethylene glycol molecules: " DEGNMOL
    		read -p "How many triethylene glycol molecules: " REGNMOL
		read -p "How many tetraethylene glycol molecules: " TEGNMOL
		read -p "How many pentaethylene glycol molecules: " PEGNMOL
		read -p "How many hexaethylene glycol molecules: " XEGNMOL
		read -p "How many heptaethylene glycol molecules: " HEGNMOL
		read -p "How many octaethylene glycol molecules: " OEGNMOL
		read -p "How many nonaethylene glycol molecules: " NEGNMOL
	fi 
	read -p "How many total polyethylene glycol molecules: " POLYEGNMOL
	read -p "How many of those molecules are water: " WATERNMOL
fi
##### END OF STEP

####calculating nmol values based on stock mixtures
if [ $STOCKMIX = "PEG200" ]; then
	MEGNMOL=0
	DEGNMOL=$(printf "%.0f" $(echo "($POLYEGNMOL * 0.03453)" | bc -l))
	REGNMOL=$(printf "%.0f" $(echo "($POLYEGNMOL * 0.22270)" | bc -l))
	TEGNMOL=$(printf "%.0f" $(echo "($POLYEGNMOL * 0.31989)" | bc -l))
	PEGNMOL=$(printf "%.0f" $(echo "($POLYEGNMOL * 0.24644)" | bc -l))
	XEGNMOL=$(printf "%.0f" $(echo "($POLYEGNMOL * 0.12901)" | bc -l))
	HEGNMOL=$(printf "%.0f" $(echo "($POLYEGNMOL * 0.04727)" | bc -l))
	OEGNMOL=0
	NEGNMOL=0

	NMOLCHECK=$(($MEGNMOL + $DEGNMOL + $REGNMOL + $TEGNMOL + $PEGNMOL + $XEGNMOL + $HEGNMOL + $OEGNMOL + $NEGNMOL)) 

	if [ $NMOLCHECK \> $POLYEGNMOL ]; then
		MORE=$(($NMOLCHECK - $POLYEGNMOL))
		HEGNMOL=$(($HEGNMOL - $MORE))
		NMOLCHECK=$(($MEGNMOL + $DEGNMOL + $REGNMOL + $TEGNMOL + $PEGNMOL + $XEGNMOL + $HEGNMOL + $OEGNMOL + $NEGNMOL)) 
	echo "removed $MORE molecules of HEG to achieve requested number of molecules"
	
	elif [ $NMOLCHECK \< $POLYEGNMOL ]; then
		LESS=$(($POLYEGNMOL - $NMOLCHECK))
		DEGNMOL=$(($DEGNMOL + $LESS))
		NMOLCHECK=$(($MEGNMOL + $DEGNMOL + $REGNMOL + $TEGNMOL + $PEGNMOL + $XEGNMOL + $HEGNMOL + $OEGNMOL + $NEGNMOL))
	echo "added $LESS molecules of DEG to achieve requested number of molecules" 
	
	else 
		echo "successfully calculated molecule numbers from mole fractions with no changes"
	fi 
	
	if [ $NMOLCHECK != $POLYEGNMOL ]; then
		echo "error with mole fractions computation" 
		exit 1 
	fi 

	echo "MEGNMOL = $MEGNMOL, DEGNMOL = $DEGNMOL, REGNMOL = $REGNMOL, TEGNMOL = $TEGNMOL, PEGNMOL = $PEGNMOL, XEGNMOL = $XEGNMOL, HEGNMOL = $HEGNMOL, OEGNMOL = $OEGNMOL, NEGNMOL = $NEGNMOL, POLYEGNMOL = $NMOLCHECK"

fi
####END OF STEP

####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")
VARIABLES=("BASENAME" "MODEL" "TEMP" "FORCEFIELD" "DENSITY" "TOTALNMOL" "STOCKMIX" "MEGNMOL" "DEGNMOL" "REGNMOL" "TEGNMOL" "PEGNMOL" "XEGNMOL" "HEGNMOL" "OEGNMOL" "NEGNMOL" "POLYEGNMOL" "WATERNMOL")
VALUES=("$BASENAME" "$MODEL" "$TEMP" "$FORCEFIELD" "$DENSITY" "$TOTALNMOL" "$STOCKMIX" "$MEGNMOL" "$DEGNMOL" "$REGNMOL" "$TEGNMOL" "$PEGNMOL" "$XEGNMOL" "$HEGNMOL" "$OEGNMOL" "$NEGNMOL" "$POLYEGNMOL" "$WATERNMOL")
#18 variables
####END OF STEP

#####check inputted variables
#are all variables defined
index=0 
while [ $index -lt 18 ]; do 
	if [ "${VALUES[$index]}" = "" ] ; then
		echo "${VARIABLES[$index]} was not defined"
		exit 1
	fi
((index++))
done

#was a valid model selected
if [ $MODEL != "bulk" ] && [ $MODEL != "spce" ] && [ $MODEL != "tip4p2005" ]; then
	echo "You have chosen an unsupported model. Please correct -m input"
	exit 1
fi

#was a valid forcefield selected
if [ $FORCEFIELD != "OPLS" ] && [ $FORCEFIELD != "OPLSMOD" ]; then 
	echo "You have chosen an unsupported forcefield. Please correct -f input"
	exit 1
fi

#density cannot be 0 
if [ $DENSITY = 0 ]; then
	echo "density must be nonzero. Please correct -d input"
	exit 1
fi

#do molecule totals add up correctly
if [ $(( $POLYEGNMOL + $WATERNMOL )) -ne $TOTALNMOL ]; then
        echo "Sum of PolyEG and water molecules is not equal to total number of molecules specified. Please correct input."
        exit 1
fi

#was a valid stock mixture selected
if [ $STOCKMIX != "none" ] && [ $STOCKMIX != "PEG200" ]; then
	echo "You have chosen an unsupported stock mixture. Please correct -s input."
	exit 1 
fi

#do molecule totals add up correctly
if [ $(( $MEGNMOL + $DEGNMOL + $REGNMOL + $TEGNMOL + $PEGNMOL + $XEGNMOL + $HEGNMOL + $OEGNMOL + $NEGNMOL )) -ne $POLYEGNMOL ]; then
	echo "Sum of polyethylene glycol molecules does not equal total specified with -p. Please correct input."
	exit 1
fi

#does number of water molecules match model
if [ $MODEL = "bulk" ] && [ $WATERNMOL -ne 0 ]; then
	echo "bulk model indicates no water molecules. Please correct either -m or -w input."
	exit 1
elif [ $MODEL != "bulk" ] && [ $WATERNMOL = 0 ]; then
	echo "conflicting information from -m and -w"
	exit 1
else
	echo ""
fi

#export variables to be used by the analysis script
cp ~/BrockportMDFiles/shell_scripts/variables.txt $WORKDIR 
index=0
while [ $index -lt 18 ]; do
	sed -i "/${VARIABLES[$index]}/s/${VARIABLES[$index]}/${VALUES[$index]}/g" "$WORKDIR"/variables.txt
	((index++))
done
####END OF STEP

####create top file
# to understand sed replacing, see point 33 at https://catonmat.net/sed-one-liners-explained-part-one
# note sed -i means editing existing files (no backup)

if [ $FORCEFIELD = "OPLS" ]; then
	MOLEDIR=/home/mmh/BrockportMDFiles/mixture_start_files
	cp "$MOLEDIR"/top_mixtures.top "$WORKDIR"/top_${BASENAME}.top
else
	MOLEDIR=/home/mmh/BrockportMDFiles/modified_mixture_start_files
	cp "$MOLEDIR"/top_mod_mixtures.top "$WORKDIR"/top_${BASENAME}.top
fi

#set molecule numbers for oligomers and remove unused oligomers
index=0
while [ $index -lt 9 ]; do
	if [ "${POLYNUM[$index]}" != "0" ]; then
		sed -i "/"${POLYS[$index]}"/s/"${POLYS[$index]}"/"${POLYNUM[$index]}"/g" "$WORKDIR"/top_${BASENAME}.top
	else
		sed -i "/"${POLYS[$index]}"/d" "$WORKDIR"/top_${BASENAME}.top
	fi
	((index++))
done

#set water molecules for model used and remove unused models
if [ $WATERNMOL -eq 0 ]; then
	sed -i "/WATERNMOL/d" "$WORKDIR"/top_${BASENAME}.top
elif [ $MODEL = "tip4p2005" ]; then
	sed -i "/TIP4P/s/WATERNMOL/$WATERNMOL/g" "$WORKDIR"/top_${BASENAME}.top
	sed -i "/SPCE/d" "$WORKDIR"/top_${BASENAME}.top
else
	sed -i "/SPCE/s/WATERNMOL/$WATERNMOL/g" "$WORKDIR"/top_${BASENAME}.top
	sed -i "/TIP4P/d" "$WORKDIR"/top_${BASENAME}.top
fi

####END OF STEP

####use insert-molecules to creat start.gro
# setting molecular weights and scaling of density
# molecules inserted into a box such that [BOXDENSITY] = [DENSITY] * [SCALE]
AVAGADRO=$(echo "6.02214076 * 10^23" | bc)

#molecular weight and density scale for each PolyEG and water
#MEG
	MWM=62.0676
	DSCALEM=0.65
#DEG
	MWD=106.12
	DSCALED=0.60
#REG
	MWR=150.1724
	DSCALER=0.55
#TEG
	MWT=194.2248
	DSCALET=0.50
#PEG
	MWP=238.2772
	DSCALEP=0.45
#XEG
	MWX=282.3296
	DSCALEX=0.40
#HEG
	MWH=326.3820
	DSCALEH=0.35
#OEG
	MWO=370.4344
	DSCALEO=0.30
#NEG
	MWN=414.4868
	DSCALEN=0.25
#water
	MWW=18.02
	DSCALEW=0.9

# bc does not handle fractional powers; use exp(ln(ans)/3) to get ans^(1/3) with math library (-l)
#formula for single oligomers:
#BOXSIDE=$(echo "scale=10; e( l(($TOTALNMOL * $MW * 10^21)/($AVAGADRO * $DENSITY * $DENSITYSCALE))/3 )" | bc -l )  # in nm
BOXSIDE=$(echo "scale=10; e( l( ((($MEGNMOL * $MWM)/$DSCALEM) + (($DEGNMOL * $MWD)/$DSCALED) + (($REGNMOL * $MWR)/$DSCALER) + (($TEGNMOL * $MWT)/$DSCALET) + (($PEGNMOL * $MWP)/$DSCALEP) + (($XEGNMOL * $MWX)/$DSCALEX) + (($HEGNMOL * $MWH)/$DSCALEH) + (($OEGNMOL * $MWO)/$DSCALEO) + (($NEGNMOL * $MWN)/$DSCALEN) + (($WATERNMOL * $MWW)/$DSCALEW)) * (10^21/($AVAGADRO * $DENSITY)) )/3 )" | bc -l)
echo "The dimensions of the box for gmx insert-molecules is $BOXSIDE nm x $BOXSIDE nm x $BOXSIDE nm."

#test 
index=0
while [ $index -lt 9 ] ; do 
if [ "${POLYNUM[$index]}" != "0" ]; then
echo " "${POLYNAME[$index]}"="${POLYNUM[$index]}" "
fi
((index++))
done

# insert molecules with final output saved at start.gro
mkdir $WORKDIR/gro_intermediates
index=0
intermediate=1
while  [ $index -lt 9 ]; do
	if [ "${POLYNUM[$index]}" != "0" ]; then
		if [ -e $WORKDIR/gro_intermediates/intermediate_$(expr $intermediate - 1).gro ]; then
			gmx insert-molecules \
			-ci $MOLEDIR/single_"${POLYNAME[$index]}".gro \
			-f $WORKDIR/gro_intermediates/intermediate_$(expr $intermediate - 1).gro \
			-o $WORKDIR/gro_intermediates/intermediate_${intermediate}.gro \
			-nmol "${POLYNUM[$index]}" \
			-try 10000
			((intermediate++))
		else
			gmx insert-molecules \
			-ci $MOLEDIR/single_"${POLYNAME[$index]}".gro \
			-o $WORKDIR/gro_intermediates/intermediate_${intermediate}.gro \
			-nmol "${POLYNUM[$index]}" \
			-box $BOXSIDE $BOXSIDE $BOXSIDE \
			-try 10000
			((intermediate++))
		fi 
	fi
	((index++))
done

if [ $WATERNMOL -gt 0 ]; then
	if [ -e $WORKDIR/gro_intermediates/intermediate_$(expr $intermediate - 1).gro ]; then
		gmx insert-molecules \
		-ci $MOLEDIR/single_${MODEL}.gro \
		-f $WORKDIR/gro_intermediates/intermediate_$(expr $intermediate - 1).gro \
		-o $WORKDIR/gro_intermediates/intermediate_${intermediate}.gro \
		-nmol $WATERNMOL \
		-try 10000
		((intermediate++))
	else
		gmx insert-molecules \
		-ci $MOLEDIR/single_${MODEL}.gro \
		-o $WORKDIR/gro_intermediates/intermediate_${intermediate}.gro \
		-nmol $WATERNMOL \
		-box $BOXSIDE $BOXSIDE $BOXSIDE \
		-try 10000
		((intermediate++))
	fi
fi
mv $WORKDIR/gro_intermediates/intermediate_$(expr $intermediate - 1).gro $WORKDIR/start.gro

####END OF STEP

if [ -e $WORKDIR/start.gro ] && [ -e $WORKDIR/top_${BASENAME}.top ] && [ -e "$WORKDIR"/em.mdp ];then
	echo "starting energy minimization"
else
	echo "missing one or more files needed to start energy minimization"
	exit 1
fi

####energy minimization 
mkdir -p "$WORKDIR"/1_EM/out
EM="$WORKDIR"/1_EM
mv "$WORKDIR"/em.mdp "$EM"
mv "$WORKDIR"/start.gro "$EM"

gmx grompp \
	-f "$EM"/em.mdp \
	-c "$EM"/start.gro \
	-p "$WORKDIR"/top_${BASENAME}.top \
	-o "$EM"/em.tpr \
	-po "$EM"/out/mdout.mdp        

gmx mdrun \
	-s "$EM"/em.tpr \
	-c "$EM"/out/em.gro \
	-o "$EM"/out/traj.trr \
	-e "$EM"/out/energy.edr \
	-g "$EM"/out/log.log \
	-pme cpu -pmefft cpu -update cpu -bonded cpu -nb gpu 

# energy minimization potential energy analysis
	mkdir -p "$EM"/analysis
	echo "Potential" | gmx energy -f "$EM"/out/energy.edr -o "$EM"/analysis/potential.xvg
##### END STEP

##### NPT equilibration run
mkdir -p "$WORKDIR"/2_equil_NPT_${TEMP}K/out
EQUIL="$WORKDIR"/2_equil_NPT_${TEMP}K
cp "$EM"/out/em.gro "$EQUIL"
mv "$WORKDIR"/eq_NPT_${TEMP}K.mdp "$EQUIL"

if [ -e $EQUIL/em.gro ] && [ -e $WORKDIR/top_${BASENAME}.top ] && [ -e "$EQUIL"/eq_NPT_${TEMP}K.mdp ];then
	echo "starting equilibration"
else
	echo "missing one or more files needed to start equilibration"
	exit 1
fi

gmx grompp \
	-f "$EQUIL"/eq_NPT_${TEMP}K.mdp  \
	-c "$EQUIL"/*.gro \
	-p "$WORKDIR"/top_${BASENAME}.top \
	-o "$EQUIL"/npt.tpr \
	-po "$EQUIL"/out/mdout.mdp \
	-maxwarn 1

gmx mdrun \
	-s "$EQUIL"/npt.tpr \
	-c "$EQUIL"/out/npt.gro \
	-x "$EQUIL"/out/traj.xtc \
	-e "$EQUIL"/out/energy.edr \
	-g "$EQUIL"/out/log.log \
	-cpo "$EQUIL"/out/npt.cpt \
	-cpt 1 \
	-pme cpu -pmefft cpu -update cpu -bonded cpu -nb gpu -v

# NPT equilibration energy analysis
mkdir -p "$EQUIL"/analysis
echo "Potential"   | gmx energy -f "$EQUIL"/out/energy.edr -o "$EQUIL"/analysis/potential.xvg
echo "Temperature" | gmx energy -f "$EQUIL"/out/energy.edr -o "$EQUIL"/analysis/temperature.xvg
echo "Pressure"    | gmx energy -f "$EQUIL"/out/energy.edr -o "$EQUIL"/analysis/pressure.xvg
echo "Volume"      | gmx energy -f "$EQUIL"/out/energy.edr -o "$EQUIL"/analysis/volume.xvg
echo "Density"     | gmx energy -f "$EQUIL"/out/energy.edr -o "$EQUIL"/analysis/density.xvg

# Choose frame from NPT of average density for NVT
# use automatic feature of boxcut.sh for cutting frame
bash /home/mmh/BrockportMDFiles/shell_scripts/boxcut.sh -a -p "$EQUIL"/
# only use manual cutting of frame if automatic fails
if [ ! -f "$EQUIL"/equilibrium.gro ]; then
	bash /home/mmh/BrockportMDFiles/shell_scripts/boxcut.sh -f 0.5 -p "$EQUIL"/
fi
##### END STEP

##### NVT production run

mkdir -p "$WORKDIR"/3_prod_NVT_${TEMP}K/out
PROD="$WORKDIR"/3_prod_NVT_${TEMP}K
cp "$EQUIL"/equilibrium.gro "$PROD"
mv "$WORKDIR"/md_NVT_${TEMP}K.mdp "$PROD"

if [ -e $PROD/equilibrium.gro ] && [ -e $WORKDIR/top_${BASENAME}.top ] && [ -e "$PROD"/md_NVT_${TEMP}K.mdp ];then
	echo "starting production"
else
	echo "missing one or more files needed to start production"
	exit 1
fi

gmx grompp \
	-f "$PROD"/md_NVT_${TEMP}K.mdp \
	-c "$PROD"/equilibrium.gro \
	-p "$WORKDIR"/top_${BASENAME}.top \
	-o "$PROD"/md.tpr \
	-po "$PROD"/out/mdout.mdp 

gmx mdrun \
	-s "$PROD"/md.tpr \
	-c "$PROD"/out/md.gro \
	-x "$PROD"/out/traj.xtc \
	-e "$PROD"/out/energy.edr \
	-g "$PROD"/out/log.log \
	-cpo "$PROD"/out/md.cpt \
	-cpt 1 \
	-pme cpu -pmefft cpu -update cpu -bonded cpu -nb gpu -v

##### END NVT

##### ANALYSIS

if [ -e "$PROD"/out/md.gro ];then
	echo "starting post-production analysis"
else
	echo "production did not finish, analysis cannot be done"
	exit 1
fi

cp ~/BrockportMDFiles/shell_scripts/analysis.sh $WORKDIR
bash analysis.sh 
##### END ANALYSIS
