#!/bin/bash

#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

#check that all requested variables are defined
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")

if [ -e variables.txt ]; then
	echo "variables.txt exists."
	index=0
	while [ $index -lt 17 ]; do 
		if [ "${VALUES[$index]}" = "" ] ; then
			echo "${VARIABLES[$index]} is not defined in variables.txt" 
			exit 1 
		fi
	((index++))
	done
else
	echo "variables.txt does not exist" 
	exit 1
fi

WORKDIR="$(dirname "$(readlink -f "$0")")"
EM="$WORKDIR"/1_EM
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")

#counting: move past index groups for the systems, PEG, and water, as they will not be used 
#$start is the first customized index group
start=2
index=0 
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then 
		((start++))
	fi
((index++))
done

if [ ${WATERNMOL} -gt 0 ]; then
	((start++))
fi

##### ANALYSIS
# NVT production energy, mean squared displacement, and green-kubo/einstein viscosity analysis
mkdir -p "$PROD"/analysis
echo "Potential"   | gmx energy -f "$PROD"/out/energy.edr -o "$PROD"/analysis/potential.xvg
echo "Temperature" | gmx energy -f "$PROD"/out/energy.edr -o "$PROD"/analysis/temperature.xvg
echo "Pressure"    | gmx energy -f "$PROD"/out/energy.edr -o "$PROD"/analysis/pressure.xvg

###MSD
# gmx msd for the system, all individual oligomers, water, and the PolyEG as a group
#repeat calculations will occur if system is bulk/water
#run gmx msd for the system 
echo "0" | gmx msd -f "$PROD"/out/traj.xtc -s "$PROD"/md.tpr -o "$PROD"/analysis/msd_system.xvg
#run gmx msd for each oligomer
#variable concantenation for indexing: \n is used as a line break. the variable concantenates with itself until all inputs are gathered
#then the variable is used as the input to make the index file 
index=0
msdgroup=2 
msd=" "
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then 
		echo "$msdgroup" | gmx msd -f "$PROD"/out/traj.xtc -s "$PROD"/md.tpr -o "$PROD"/analysis/msd_"${POLYNAME[$index]}".xvg 
		msd="${msd}r ${POLYNAME[index]} | "
		((msdgroup++)) 
	fi 
	((index++))
done 

#run gmx msd for water 
if [ $WATERNMOL -gt 0 ]; then 
	echo "$msdgroup" | gmx msd -f "$PROD"/out/traj.xtc -s "$PROD"/md.tpr -o "$PROD"/analysis/msd_${MODEL}.xvg 
	((msdgroup++))  
fi 

#index Polyeg and run gmx msd on new group 
if [ $POLYEGNMOL -gt 0 ]; then
	if [ $MODEL != "bulk" ]; then
		echo -e  "$msd \n q" | gmx make_ndx -f "$EM"/start.gro -o "$PROD"/analysis/msd.ndx 
		echo "$msdgroup" | gmx msd -f "$PROD"/out/traj.xtc -s "$PROD"/md.tpr -n "$PROD"/analysis/msd.ndx -o "$PROD"/analysis/msd_PolyEG.xvg 
	fi
fi

### END MSD
 
###DISTANCE
#distance index 
index=0
distance=" "
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then 
	distance="${distance}a OA${POLYA[index]}1 | a OA${POLYA[index]}2 \n "
	fi
	((index++))
done 
echo -e "${distance} q" | gmx make_ndx -f "$EM"/start.gro -o "$PROD"/analysis/distance.ndx 

#run gmx distance and save to a txt 
#save the output of a command to a text file: command >> output.txt (use > to overwrite the file)
#similarly: << for an input file
index=0
group=$start
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then 
	gmx distance\
		-f "$PROD"/out/traj.xtc\
		-s "$PROD"/md.tpr\
		-n "$PROD"/analysis/distance.ndx\
		-oav "$PROD"/analysis/distance_"${POLYNAME[$index]}".xvg\
		-select $group >> "$PROD"/analysis/distance.txt
	((group++))
	fi
((index++))
done
###END DISTANCE

###RDF
mkdir $PROD/analysis/rdf
mkdir $PROD/analysis/rdf/total
#intermolecular total 
group=$start
index=0 
rdf=" "
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then 
	rdf="${rdf}a OA${POLYA[$index]}1 \n a OA${POLYA[$index]}2 \n"
	((group += 2))
		if [ ${POLYNAME[$index]} != "MEG" ]; then 
		ether="t peg_OE & r ${POLYNAME[$index]}\n name $group OE_${POLYNAME[$index]}\n"
		rdf="${rdf}${ether}"
		((group++))
		fi
	fi
((index++))
done 
echo -e "${rdf}\n q" | gmx make_ndx\
	-f "$PROD"/md.tpr\
	-o "$PROD"/analysis/rdf/rdf_total.ndx
#echo -e "rdf input:$rdf" 

mkdir "$PROD"/analysis/rdf/"$BASENAME"_pngs

group=$start 
index=0
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then 
		gmx rdf\
			-f "$PROD"/out/traj.xtc\
			-s "$PROD"/md.tpr\
			-n "$PROD"/analysis/rdf/rdf_total.ndx\
			-o "$PROD"/analysis/rdf/total/rdf_${POLYNAME[$index]}_total.xvg\
			-ref $group -sel $(expr $group + 1) $(expr $group + 2)
		((group += 3))
	fi
((index++))
done

mkdir $PROD/analysis/rdf/intra
#intramolecular  
group=$start
index=0 
rdf=" "
offset=0
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ]; then 
		mol=1
		if [ ${POLYNUM[$index]} -lt 10 ]; then
			while [ $mol -le ${POLYNUM[$index]} ]; do 
			num=$(($mol + $offset))
			rdf="${rdf}a OA${POLYA[$index]}1 & ri "$num"\n name $group OA"${POLYA[$index]}"1_m"$num"\n"
			((group++))
			rdf="${rdf}a OA${POLYA[$index]}2 & ri "$num"\n name $group OA"${POLYA[$index]}"2_m"$num"\n"
			((group++))		
				if [ ${POLYNAME[$index]} != "MEG" ]; then
				ether="t peg_OE & r ${POLYNAME[$index]} &ri "$num"\n name $group OE_"${POLYNAME[$index]}"_m"$num"\n"
				((group++))
				rdf="${rdf}${ether}"
				fi			
			((mol++))
			done
		else
			while [ $mol -le 10 ]; do
			num=$(($mol + $offset))
			rdf="${rdf}a OA${POLYA[$index]}1 & ri "$num"\n name $group OA"${POLYA[$index]}"1_m"$num"\n"
			((group++))
			rdf="${rdf}a OA${POLYA[$index]}2 & ri "$num"\n name $group OA"${POLYA[$index]}"2_m"$num"\n"
			((group++))		
				if [ ${POLYNAME[$index]} != "MEG" ]; then
				ether="t peg_OE & r ${POLYNAME[$index]} &ri "$num"\n name $group OE_"${POLYNAME[$index]}"_m"$num"\n"
				((group++))
				rdf="${rdf}${ether}"
				fi			
			((mol++))
			done
		fi
	offset=(${offset}+${POLYNUM[$index]})
	fi 
((index++))
done 
echo -e "${rdf}\n q" | gmx make_ndx -f "$PROD"/md.tpr -o "$PROD"/analysis/rdf/rdf_intra.ndx
#echo -e "rdf intra input:\n$rdf"

group=$start 
index=0
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then 
	mol=1
		if [ ${POLYNUM[$index]} -lt 10 ]; then
			while [ $mol -le ${POLYNUM[$index]} ]; do
			gmx rdf\
				-f "$PROD"/out/traj.xtc\
				-s "$PROD"/md.tpr\
				-n "$PROD"/analysis/rdf/rdf_intra.ndx\
				-o $PROD/analysis/rdf/intra/rdf_${POLYNAME[$index]}_intra_${mol}.xvg\
				-ref $group -sel $(expr $group + 1) $(expr $group + 2)
			((group += 3))
			((mol++))
			done 
		else
			while [ $mol -le 10 ]; do
			gmx rdf\
				-f "$PROD"/out/traj.xtc\
				-s "$PROD"/md.tpr\
				-n "$PROD"/analysis/rdf/rdf_intra.ndx\
				-o $PROD/analysis/rdf/intra/rdf_${POLYNAME[$index]}_intra_${mol}.xvg\
				-ref $group -sel $(expr $group + 1) $(expr $group + 2)
			((group += 3))
			((mol++))
			done 
		fi
	fi
((index++))
done

#intermolecular interresidual

#makes a list of all OA1 atoms present
#ie OAD1 OAR1 etc 
index=0
OA1=" "
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then 
	OA1="${OA1}a OA${POLYA[$index]}1|"
	fi
((index++))
done 

#ref group OA1 for the residue
#sel group 1 ethers in other residues
#sel group 2 OA1 atoms in other residues
group=$start
index=0 
rdf=" "
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then 
	rdf="${rdf}a OA${POLYA[$index]}1 \n"
	((group++))		
	ether="t peg_OE &! r ${POLYNAME[$index]}\n name $group OE_!"${POLYNAME[$index]}"\n"
	((group++))
	rdf="${rdf}${ether}${OA1}&! r ${POLYNAME[$index]}\n name $group OA1_!"${POLYNAME[$index]}"\n"
	((group++))
	fi
((index++))
done
#echo -e "rdf inter_dif input:\n$rdf"
echo -e "${rdf}\n q" | gmx make_ndx -f "$PROD"/md.tpr -o "$PROD"/analysis/rdf/rdf_inter_dif.ndx

group=$start 
index=0
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then 
		gmx rdf\
			-f "$PROD"/out/traj.xtc\
			-s "$PROD"/md.tpr\
			-n "$PROD"/analysis/rdf/rdf_inter_dif.ndx\
			-o $PROD/analysis/rdf/rdf_${POLYNAME[$index]}_inter_dif.xvg\
			-ref $group -sel $(expr $group + 1) $(expr $group + 2)
		((group += 3))
		xmgrace\
			-nxy $PROD/analysis/rdf/rdf_${POLYNAME[$index]}_inter_dif.xvg\
			-hdevice PNG -hardcopy\
			-printfile $PROD/analysis/rdf/"$BASENAME"_pngs/rdf_${POLYNAME[$index]}_inter_dif.png
	fi
((index++))
done

#begin python analysis
#rdf_ave.py:
#creates an xvg that is the average of the 10 intra files
#also subtracts intra from total to give the inter_same graph
#rdf_norm.py
#creates a normalized intra average from intra_ave.xvg

#export variables python needs from this script
export PROD

#the python code asks the user to input the residue acronym. 
#input only residues which are present or the code will error
index=0
while [ $index -lt 9 ]; do 
export ${POLYS[index]}
	if [ "${POLYNUM[$index]}" != "0" ]; then
		echo "${POLYNAME[$index]}" | python3 ~/BrockportMDFiles/python_scripts/rdf_ave.py 
		echo "${POLYNAME[$index]}" | python3 ~/BrockportMDFiles/python_scripts/rdf_norm.py 
	fi
((index++))
done

#the newly constructed xvg files copy the header from an existing file. change the data labels. 
index=0
num=1
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ]; then
		sed -i "/subtitle/s/OA${POLYA[$index]}1_m$num/OA${POLYA[$index]}1/g" $PROD/analysis/rdf/rdf_${POLYNAME[$index]}_intra_ave.xvg 
		sed -i "/s0/s/OA${POLYA[$index]}2_m$num/OA${POLYA[$index]}2/g" $PROD/analysis/rdf/rdf_${POLYNAME[$index]}_intra_ave.xvg 
		sed -i "/s1/s/OE_${POLYNAME[$index]}_m$num/OE_${POLYNAME[$index]}/g" $PROD/analysis/rdf/rdf_${POLYNAME[$index]}_intra_ave.xvg 
		sed -i "/subtitle/s/OA${POLYA[$index]}1_m$num/OA${POLYA[$index]}1/g" $PROD/analysis/rdf/rdf_${POLYNAME[$index]}_intra_norm.xvg 
		sed -i "/s0/s/OA${POLYA[$index]}2_m$num/OA${POLYA[$index]}2/g" $PROD/analysis/rdf/rdf_${POLYNAME[$index]}_intra_norm.xvg 
		sed -i "/s1/s/OE_${POLYNAME[$index]}_m$num/OE_${POLYNAME[$index]}/g" $PROD/analysis/rdf/rdf_${POLYNAME[$index]}_intra_norm.xvg 
		num=$(( $num + ${POLYNUM[$index]} ))
		xmgrace -nxy $PROD/analysis/rdf/rdf_${POLYNAME[$index]}_intra_ave.xvg\
			-hdevice PNG -hardcopy\
			-printfile $PROD/analysis/rdf/"$BASENAME"_pngs/rdf_${POLYNAME[$index]}_intra_ave.png
		xmgrace -nxy $PROD/analysis/rdf/rdf_${POLYNAME[$index]}_inter_same.xvg\
			-hdevice PNG -hardcopy\
			-printfile $PROD/analysis/rdf/"$BASENAME"_pngs/rdf_${POLYNAME[$index]}_inter_same.png
		xmgrace -nxy $PROD/analysis/rdf/rdf_${POLYNAME[$index]}_intra_norm.xvg\
			-hdevice PNG -hardcopy\
			-printfile $PROD/analysis/rdf/"$BASENAME"_pngs/rdf_${POLYNAME[$index]}_intra_norm.png
	fi
((index++))
done 

#water rdf
if [ $WATERNMOL -gt 0 ]; then

	group=$start
	index=0 
	rdf=" "
	water="a OW\n "
	((group++))

	while [ $index -lt 9 ]; do 
		if [ "${POLYNUM[$index]}" != "0" ] ; then 
			rdf="${rdf}t peg_OA & r ${POLYNAME[$index]}\n name $group OA_${POLYNAME[$index]}\n"
			((group++))
			ether="t peg_OE & r ${POLYNAME[$index]}\n name $group OE_${POLYNAME[$index]}\n"
			((group++))
			rdf="${rdf}${ether}"
		fi
		((index++))
	done

	rdf=${water}${rdf}	

	echo -e "${rdf}\n q" | gmx make_ndx\
		-f "$PROD"/md.tpr\
		-o "$PROD"/analysis/rdf/rdf_water.ndx
#water-water rdf	
	gmx rdf\
		-f "$PROD"/out/traj.xtc\
		-s "$PROD"/md.tpr\
		-n "$PROD"/analysis/rdf/rdf_water.ndx\
		-o "$PROD"/analysis/rdf/rdf_W_W.xvg\
		-ref $start -sel $start
		>> "$PROD"/analysis/hbond/possible.txt
	xmgrace\
		-nxy "$PROD"/analysis/rdf/rdf_W_W.xvg\
		-hdevice PNG -hardcopy\
		-printfile $PROD/analysis/rdf/"$BASENAME"_pngs/rdf_W_W.png	

#water-OH and water_OE rdf
	group=$start 
	index=0
	((group++))
	while [ $index -lt 9 ]; do 
		if [ "${POLYNUM[$index]}" != "0" ] ; then 
		gmx rdf\
			-f "$PROD"/out/traj.xtc\
			-s "$PROD"/md.tpr\
			-n "$PROD"/analysis/rdf/rdf_water.ndx\
			-o "$PROD"/analysis/rdf/rdf_W_"${POLYNAME[$index]}".xvg\
			-ref $start -sel $group $(expr $group + 1)
			>> "$PROD"/analysis/hbond/possible.txt
		((group+= 2))		
		xmgrace\
			-nxy "$PROD"/analysis/rdf/rdf_W_"${POLYNAME[$index]}".xvg\
			-hdevice PNG -hardcopy\
			-printfile $PROD/analysis/rdf/"$BASENAME"_pngs/rdf_W_"${POLYNAME[$index]}".png
		fi
	((index++))
	done
		

fi

###END RDF

###HBOND
mkdir $PROD/analysis/hbond
mkdir $PROD/analysis/hbond/life

#water-water hbond
if [ $WATERNMOL -gt 0 ]; then
	echo -e "$(expr $start - 1)\n$(expr $start - 1)\n" | gmx hbond\
		-f "$PROD"/out/traj.xtc\
		-s "$PROD"/md.tpr\
		-num "$PROD"/analysis/hbond/hbond_W_W.xvg\
		-ac "$PROD"/analysis/hbond/life/life_W_W.xvg\
		>> "$PROD"/analysis/hbond/possible.txt
fi 

#intermolecular (within the same residue: includes intramolecular, which need to be subtracted away)
#index
group=$start
index=0 
hbond=" "
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then 
	hbond="${hbond}t peg_OA | t peg_HA & r "${POLYNAME[$index]}" \n name $group OH_"${POLYNAME[$index]}"\n"
	((group++))
		if [ "${POLYNAME[$index]}" != "MEG" ]; then 
		ether="t peg_OE & r "${POLYNAME[$index]}"\n name $group OE_"${POLYNAME[$index]}"\n"
		((group++))
		hbond="${hbond}${ether}"
		fi
	fi 
((index++))
done 
echo -e "${hbond}\n q" | gmx make_ndx -f "$PROD"/md.tpr -o "$PROD"/analysis/hbond/hbond_total.ndx
#echo -e "hbond total input:\n$hbond"  

#hbond
group=$start 
index=0
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then 
		echo -e "$group\n$group\n" | gmx hbond\
			-f "$PROD"/out/traj.xtc\
			-s "$PROD"/md.tpr\
			-n "$PROD"/analysis/hbond/hbond_total.ndx\
			-num "$PROD"/analysis/hbond/hbond_"${POLYNAME[$index]}"_total_OH.xvg\
			-ac "$PROD"/analysis/hbond/life/life_"${POLYNAME[$index]}"_total_OH.xvg\
			>> "$PROD"/analysis/hbond/possible.txt
		echo -e "$group\n$(expr $group + 1)\n" | gmx hbond\
			-f "$PROD"/out/traj.xtc\
			-s "$PROD"/md.tpr\
			-n "$PROD"/analysis/hbond/hbond_total.ndx\
			-num "$PROD"/analysis/hbond/hbond_"${POLYNAME[$index]}"_total_OE.xvg\
			-ac "$PROD"/analysis/hbond/life/life_"${POLYNAME[$index]}"_total_OE.xvg\
			>> "$PROD"/analysis/hbond/possible.txt
		let group="$group + 2"
	fi
((index++))
done

#intramolecular  
#index
group=$start
index=0 
hbond=" "
offset=0
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ]; then 
		mol=1
			if [ ${POLYNUM[$index]} -lt 10 ]; then 
				while [ $mol -le ${POLYNUM[$index]} ]; do 
				let num="$mol + $offset"
				hbond="${hbond}t peg_OA | t peg_HA & r "${POLYNAME[$index]}" &ri $num \n name $group OH_"${POLYNAME[$index]}"_m"$num"\n"
				((group++))
					if [ "${POLYNAME[$index]}" != "MEG" ]; then
					ether="t peg_OE & r "${POLYNAME[$index]}" &ri $num\n name $group OE_"${POLYNAME[$index]}"_m"$num"\n"
					((group++))
					hbond="${hbond}${ether}"
					fi			
				((mol++))
				done
			else
				while [ $mol -le 10 ]; do
				let num="$mol + $offset"
				hbond="${hbond}t peg_OA | t peg_HA & r "${POLYNAME[$index]}" &ri $num \n name $group OH_"${POLYNAME[$index]}"_m"$num"\n"
				((group++))
					if [ "${POLYNAME[$index]}" != "MEG" ]; then
					ether="t peg_OE & r "${POLYNAME[$index]}" &ri $num\n name $group OE_"${POLYNAME[$index]}"_m"$num"\n"
					((group++))
					hbond="${hbond}${ether}"
					fi			
				((mol++))
				done
			fi
	offset=(${offset}+"${POLYNUM[$index]}")
	fi 
((index++))
done 
echo -e "${hbond}\n q" | gmx make_ndx -f "$PROD"/md.tpr -o "$PROD"/analysis/hbond/hbond_intra.ndx
#echo -e "hbond input intra:\n$hbond" 

#hbond
group=$start 
index=0
mkdir $PROD/analysis/hbond/intra
mkdir $PROD/analysis/hbond/intra/life
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then 
		mol=1 
			if [ ${POLYNUM[$index]} -lt 10 ]; then
				while [ $mol -le ${POLYNUM[$index]} ]; do 
				echo -e "$group\n$group\n" | gmx hbond\
					-f "$PROD"/out/traj.xtc\
					-s "$PROD"/md.tpr\
					-n "$PROD"/analysis/hbond/hbond_intra.ndx\
					-num $PROD/analysis/hbond/intra/hbond_"${POLYNAME[$index]}"_intra_OH_$mol.xvg\
					-ac $PROD/analysis/hbond/intra/life/life_"${POLYNAME[$index]}"_intra_OH_$mol.xvg\
					>> "$PROD"/analysis/hbond/intra/possible_"${POLYNAME[$index]}"_intra_OH.txt
				echo -e "$group\n$(expr $group + 1)\n" | gmx hbond\
					-f "$PROD"/out/traj.xtc\
					-s "$PROD"/md.tpr\
					-n "$PROD"/analysis/hbond/hbond_intra.ndx\
					-num $PROD/analysis/hbond/intra/hbond_"${POLYNAME[$index]}"_intra_OE_$mol.xvg\
					-ac $PROD/analysis/hbond/intra/life/life_"${POLYNAME[$index]}"_intra_OE_$mol.xvg\
					>> $PROD/analysis/hbond/intra/life/life_"${POLYNAME[$index]}"_intra_OE.txt
				let group="$group + 2"
				((mol++))
				done 
			else
#saves the output to a seperate txt file to be averaged
				while [ $mol -le 10 ]; do
				echo -e "$group\n$group\n" | gmx hbond\
					-f "$PROD"/out/traj.xtc\
					-s "$PROD"/md.tpr\
					-n "$PROD"/analysis/hbond/hbond_intra.ndx\
					-num $PROD/analysis/hbond/intra/hbond_"${POLYNAME[$index]}"_intra_OH_$mol.xvg\
					-ac $PROD/analysis/hbond/intra/life/life_"${POLYNAME[$index]}"_intra_OH_$mol.xvg\
					>> "$PROD"/analysis/hbond/intra/possible_"${POLYNAME[$index]}"_intra_OH.txt
				echo -e "$group\n$(expr $group + 1)\n" | gmx hbond\
					-f "$PROD"/out/traj.xtc\
					-s "$PROD"/md.tpr\
					-n "$PROD"/analysis/hbond/hbond_intra.ndx\
					-num $PROD/analysis/hbond/intra/hbond_"${POLYNAME[$index]}"_intra_OE_$mol.xvg\
					-ac $PROD/analysis/hbond/intra/life/life_"${POLYNAME[$index]}"_intra_OE_$mol.xvg\
					>> $PROD/analysis/hbond/intra/possible_"${POLYNAME[$index]}"_intra_OE.txt
				let group="$group + 2"
				((mol++))
				done 
			fi
	fi
((index++))
done

#averaging for OH
index=0
count=0
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then
#copies the line containing the average number of H bonds to a seperate file to be read
# line reads "Average number of hbonds per timeframe # out of # possible"		 
		awk '/Average/' "$PROD/analysis/hbond/intra/possible_"${POLYNAME[$index]}"_intra_OH.txt"\
			> "$PROD/analysis/hbond/intra/bondline_"${POLYNAME[$index]}"_intra_OH.txt"
		sum=0
		count=0
#reads each line saving each word to the corresponding variable
		while read -r average number of hbonds per timeframe bonds out of total possible; do 
#0 is saved as undefined and then cannot be used in math. this saves 0 as 0 
			bonds=$(echo "$bonds + 0" | bc -l)
#finds the sum of the bonds, for averaging
			sum=$(echo "$sum + $bonds" | bc -l)
#counts the number of entries, for averaging
			((count++))
			all=$total
		done < "$PROD/analysis/hbond/intra/bondline_"${POLYNAME[$index]}"_intra_OH.txt"
#takes the average and saves to a file 
		ave=$(echo " $sum / $count " | bc -l)
		ave=$(printf "%.4f" $ave)
#evaluationg standard deviaton
		diffsquare=0		
		while read -r average number of hbonds per timeframe bonds out of total possible; do 
#0 is saved as undefined and then cannot be used in math. this saves 0 as 0 
			bonds=$(echo "$bonds + 0" | bc -l)
			diffsquare=$(echo "$diffsquare + ($bonds - $ave)^2" | bc -l)
			((count++))
		done < "$PROD/analysis/hbond/intra/bondline_"${POLYNAME[$index]}"_intra_OH.txt"
		stdevbonds=$(echo "sqrt($diffsquare/($count-1))" | bc -l)
		stdevbonds=$(printf "%.4f" $stdevbonds)
		echo ""${POLYNAME[$index]}" intra OH: $ave +/- $stdevbonds out of $all" >> "$PROD/analysis/hbond/intra_averages.txt"
#now for hbond average lifetim
		awk '/Forward/ {print $3}' "$PROD/analysis/hbond/intra/possible_"${POLYNAME[$index]}"_intra_OH.txt"\
			> "$PROD/analysis/hbond/intra/lifeline_"${POLYNAME[$index]}"_intra_OH.txt"
		addlife=0
		count=0
		while read -r ; do 
			life=${REPLY}
			addlife=$(echo "$addlife + $life" | bc -l)
			((count++))
		done < "$PROD/analysis/hbond/intra/lifeline_"${POLYNAME[$index]}"_intra_OH.txt"
		avelife=$(echo " $addlife / $count " | bc -l)
		avelife=$(printf "%.1f" $avelife)
#evaluationg standard deviaton
		diffsquare=0		
		while read -r ; do
			life=${REPLY}
			diffsquare=$(echo "$diffsquare + ($life - $avelife)^2" | bc -l)
			((count++))
		done < "$PROD/analysis/hbond/intra/lifeline_"${POLYNAME[$index]}"_intra_OH.txt"
		stdevlife=$(echo "sqrt($diffsquare/($count-1))" | bc -l)
		stdevlife=$(printf "%.1f" $stdevlife)
		echo ""${POLYNAME[$index]}" intra OH: $avelife +/- $stdevlife ps average lifetime " >> "$PROD/analysis/hbond/intra_lifeaverages.txt"
	fi
((index++))
done

##averaging for OE
#3same as above

index=0
count=0
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then
#copies the line containing the average number of H bonds to a seperate file to be read
# line reads "Average number of hbonds per timeframe # out of # possible"		 
		awk '/Average/' "$PROD/analysis/hbond/intra/possible_"${POLYNAME[$index]}"_intra_OE.txt"\
			> "$PROD/analysis/hbond/intra/bondline_"${POLYNAME[$index]}"_intra_OE.txt"
		sum=0
		count=0
#reads each line saving each word to the corresponding variable
		while read -r average number of hbonds per timeframe bonds out of total possible; do 
#0 is saved as undefined and then cannot be used in math. this saves 0 as 0 
			bonds=$(echo "$bonds + 0" | bc -l)
#finds the sum of the bonds, for averaging
			sum=$(echo "$sum + $bonds" | bc -l)
#counts the number of entries, for averaging
			((count++))
			all=$total
		done < "$PROD/analysis/hbond/intra/bondline_"${POLYNAME[$index]}"_intra_OE.txt"
#takes the average and saves to a file 
		ave=$(echo " $sum / $count " | bc -l)
		ave=$(printf "%.4f" $ave)
#evaluationg standard deviaton
		diffsquare=0		
		while read -r average number of hbonds per timeframe bonds out of total possible; do 
#0 is saved as undefined and then cannot be used in math. this saves 0 as 0 
			bonds=$(echo "$bonds + 0" | bc -l)
			diffsquare=$(echo "$diffsquare + ($bonds - $ave)^2" | bc -l)
			((count++))
		done < "$PROD/analysis/hbond/intra/bondline_"${POLYNAME[$index]}"_intra_OE.txt"
		stdevbonds=$(echo "sqrt($diffsquare/($count-1))" | bc -l)
		stdevbonds=$(printf "%.4f" $stdevbonds)
		echo ""${POLYNAME[$index]}" intra OE: $ave +/- $stdevbonds out of $all" >> "$PROD/analysis/hbond/intra_averages.txt"
#now for hbond average lifetim
		awk '/Forward/ {print $3}' "$PROD/analysis/hbond/intra/possible_"${POLYNAME[$index]}"_intra_OE.txt"\
			> "$PROD/analysis/hbond/intra/lifeline_"${POLYNAME[$index]}"_intra_OE.txt"
		addlife=0
		count=0
		while read -r ; do 
			life=${REPLY}
			addlife=$(echo "$addlife + $life" | bc -l)
			((count++))
		done < "$PROD/analysis/hbond/intra/lifeline_"${POLYNAME[$index]}"_intra_OE.txt"
		avelife=$(echo " $addlife / $count " | bc -l)
		avelife=$(printf "%.1f" $avelife)
#evaluationg standard deviaton
		diffsquare=0		
		while read -r ; do
			life=${REPLY}
			diffsquare=$(echo "$diffsquare + ($life - $avelife)^2" | bc -l)
			((count++))
		done < "$PROD/analysis/hbond/intra/lifeline_"${POLYNAME[$index]}"_intra_OE.txt"
		stdevlife=$(echo "sqrt($diffsquare/($count-1))" | bc -l)
		stdevlife=$(printf "%.1f" $stdevlife)
		echo ""${POLYNAME[$index]}" intra OE: $avelife +/- $stdevlife ps average lifetime " >> "$PROD/analysis/hbond/intra_lifeaverages.txt"
	fi
((index++))
done


#intermolecular interresidual
#index
group=$start
index=0 
hbond=" "
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then 
	hbond="${hbond}t peg_OA | t peg_HA & r "${POLYNAME[$index]}" \n name $group OH_"${POLYNAME[$index]}"\n"
	((group++))
	hbond="${hbond}t peg_OA | t peg_HA &! r "${POLYNAME[$index]}"\n name $group OH_!"${POLYNAME[$index]}"\n"
	((group++))
		if [ "${POLYNAME[$index]}" != "MEG" ]; then 
		ether="t peg_OE &! r "${POLYNAME[$index]}"\n name $group OE_!"${POLYNAME[$index]}"\n"
		((group++))
		hbond="${hbond}${ether}"
		fi 
	fi
((index++))
done
#echo -e "hbond inter input:\n$hbond"
echo -e "${hbond}\n q" | gmx make_ndx -f "$PROD"/md.tpr -o "$PROD"/analysis/hbond/hbond_inter.ndx
#hbond
group=$start 
index=0
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then 
		echo -e "$group\n$(expr $group + 1)\n" | gmx hbond\
			-f "$PROD"/out/traj.xtc\
			-s "$PROD"/md.tpr\
			-n "$PROD"/analysis/hbond/hbond_inter.ndx\
			-num "$PROD"/analysis/hbond/hbond_"${POLYNAME[$index]}"_inter_OH.xvg\
			-ac "$PROD"/analysis/hbond/life/life_"${POLYNAME[$index]}"_inter_OH.xvg\
			>> "$PROD"/analysis/hbond/possible.txt
		echo -e "$group\n$(expr $group + 2)\n" | gmx hbond\
			-f "$PROD"/out/traj.xtc\
			-s "$PROD"/md.tpr\
			-n "$PROD"/analysis/hbond/hbond_inter.ndx\
			-num "$PROD"/analysis/hbond/hbond_"${POLYNAME[$index]}"_inter_OE.xvg\
			-ac "$PROD"/analysis/hbond/life/life_"${POLYNAME[$index]}"_inter_OE.xvg\
			>> "$PROD"/analysis/hbond/possible.txt
		if [ $WATERNMOL -gt 0 ]; then
			echo -e "$group\n$(expr $start - 1)\n" | gmx hbond\
				-f "$PROD"/out/traj.xtc\
				-s "$PROD"/md.tpr\
				-n "$PROD"/analysis/hbond/hbond_inter.ndx\
				-num "$PROD"/analysis/hbond/hbond_"${POLYNAME[$index]}"_OH_W.xvg\
				-ac "$PROD"/analysis/hbond/life/life_"${POLYNAME[$index]}"_OH_W.xvg\
				>> "$PROD"/analysis/hbond/possible.txt
		fi 
		let group="$group + 3"
	fi
((index++))
done

#ethers to water hbonds
if [ $WATERNMOL != 0 ];then
	group=$start
	index=0 
	ether=" " 
	while [ $index -lt 9 ]; do 
		if [ "${POLYNUM[$index]}" != "0" ] ; then 
			if [ "${POLYNAME[$index]}" != "MEG" ]; then 
			ether="${ether}t peg_OE & r "${POLYNAME[$index]}"\n name $group OE_"${POLYNAME[$index]}"\n"
			((group++))
			fi
		fi 
	((index++))
	done 
	echo -e "${ether}\n q" | gmx make_ndx -f "$PROD"/md.tpr -o "$PROD"/analysis/hbond/hbond_water.ndx

	group=$start 
	index=0
	while [ $index -lt 9 ]; do 
		if [ "${POLYNUM[$index]}" != "0" ] ; then 
			echo -e "$group\n$(expr $start - 1)\n" | gmx hbond\
				-f "$PROD"/out/traj.xtc\
				-s "$PROD"/md.tpr\
				-n "$PROD"/analysis/hbond/hbond_water.ndx\
				-num "$PROD"/analysis/hbond/hbond_"${POLYNAME[$index]}"_OE_W.xvg\
				-ac "$PROD"/analysis/hbond/life/life_"${POLYNAME[$index]}"_OE_W.xvg\
				>> "$PROD"/analysis/hbond/possible.txt
			((group++))
		fi
	((index++))
	done
fi

# making a nice summary file of all the intermolecular hydrogenbonding interaction info 
# all stored in file "inter_hbondlife_analysis.txt"
awk '/Selected/ {print $3}' "$PROD"/analysis/hbond/possible.txt > selectedboth.txt
awk '/Average number/ {print $7}' "$PROD"/analysis/hbond/possible.txt > hbondnumber.txt
awk '/Forward/ {print $2, $3, $4, $5}' "$PROD"/analysis/hbond/possible.txt > hbondlife.txt
sed -n 'n;p' selectedboth.txt > selectedeven.txt #Print odd numbered lines
sed -n 'p;n' selectedboth.txt > selectedodd.txt # Print odd numbered lines
paste selectedodd.txt selectedeven.txt hbondnumber.txt hbondlife.txt > "$PROD"/analysis/hbond/inter_hbondlife_analysis.txt
rm selectedboth.txt
rm hbondnumber.txt
rm hbondlife.txt
rm selectedeven.txt
rm selectedodd.txt

###END HBOND

###DIHEDRALS
mkdir -p "$PROD"/analysis/dihedrals
gmx mk_angndx -s "$PROD"/md.tpr -n "$PROD"/analysis/dihedrals/RB_dih.ndx -type ryckaert-bellemans

dihgroup=6
index=0

while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then
		echo "$dihgroup" | gmx angle\
			-n "$PROD"/analysis/dihedrals/RB_dih.ndx\
			-f "$PROD"/out/traj.xtc\
			-type ryckaert-bellemans\
			-od "$PROD"/analysis/dihedrals/"${POLYNAME[index]}"_OaCCOe_RB.xvg 
		xmgrace -nxy "$PROD"/analysis/dihedrals/"${POLYNAME[index]}"_OaCCOe_RB.xvg\
			-hdevice PNG -hardcopy\
			-printfile "$PROD"/analysis/dihedrals/"${POLYNAME[index]}"_OaCCOe_RB.png 
		echo -e "${POLYNAME[index]} (0,120)\n" >> "$PROD"/analysis/dihedrals.txt 
		gmx analyze\
			-f "$PROD"/analysis/dihedrals/"${POLYNAME[index]}"_OaCCOe_RB.xvg\
			-b 0 -e 120 -integrate >> "$PROD"/analysis/dihedrals.txt 
		echo -e "${POLYNAME[index]} (-120,0)\n" >> "$PROD"/analysis/dihedrals.txt 
		gmx analyze\
			-f "$PROD"/analysis/dihedrals/"${POLYNAME[index]}"_OaCCOe_RB.xvg\
			-b -120 -e 0 -integrate >> "$PROD"/analysis/dihedrals.txt 
		if [ ${POLYNAME[$index]} = "DEG" ]; then
			((dihgroup += 7))
		else
			((dihgroup +=8))
		fi
	fi
((index++))
done
###END DIHEDRALS

###GYRATION
mkdir $PROD/analysis/gyration

echo "q" | gmx make_ndx -f $PROD/md.tpr -o $PROD/analysis/gyration/gyrate.ndx

gyr_group=2 
index=0
while [ $index -lt 9 ]; do 
	if [ "${POLYNUM[$index]}" != "0" ] ; then 
		echo "$gyr_group" | gmx gyrate\
			-f $PROD/out/traj.xtc\
			-s $PROD/md.tpr\
			-o $PROD/analysis/gyration/gyrate_"${POLYNAME[$index]}".xvg\
			-nmol "${POLYNUM[$index]}"\
			-n $PROD/analysis/gyration/gyrate.ndx 
		echo "${POLYNAME[$index]}" >> $PROD/analysis/gyrate.txt
		gmx analyze\
			-f $PROD/analysis/gyration/gyrate_"${POLYNAME[$index]}".xvg\
			>> $PROD/analysis/gyrate.txt
		((gyr_group++))
	fi
((index++))
done

#water gyration
if [ ${WATERNMOL} -gt 0 ]; then
	echo "$(expr $start - 1)" | gmx gyrate\
		-f $PROD/out/traj.xtc\
		-s $PROD/md.tpr\
		-o $PROD/analysis/gyration/gyrate_water.xvg\
		-nmol $WATERNMOL\
		-n $PROD/analysis/gyration/gyrate.ndx 
	echo "water" >> $PROD/analysis/gyrate.txt
	gmx analyze\
		-f $PROD/analysis/gyration/gyrate_water.xvg\
		>> $PROD/analysis/gyrate.txt
	((gyr_group++))
fi

###END GYRATION

