diff --git a/Hydration_freenrg/README.md b/Hydration_freenrg/README.md
index 23c04e0..2cae265 100644
--- a/Hydration_freenrg/README.md
+++ b/Hydration_freenrg/README.md
@@ -2,8 +2,12 @@
 A set of scripts used to reproduce the hydration free energy results of [Leoffler _et al_](10.1021/acs.jctc.8b00544) using BioSimSpace. Currently designed to be run specifically on the Yoko cluster.
 
 ## Usage
-* First run "run_generate.sh", this will generate smiles strings for all molecules in the "ligands" folder, placing a file named "SMILES.csv" in the "python" folder. It will then generate perturbable systems in both vacuum and solvent, reading the syste>
-* A folder named "Systems" will be created, containing all of the pertubable systems - this is where all simulation data will be stored. Currently 3 replicas of each system are generated.
-* One system generation is complete, run "MinEqallsystems.sh", this will loop over every perturbable system in "Systems", minimising the vacuum systems, and minimising and equilibrating the solvated systems.
-* Once Minimisation and equilibration are complete, run "runallsystems.sh". This will run all vacuum and solvated systems for all replicas (expect this to take a while).
-* Once all simulations are complete, run "analyse_all.sh" to analyse all simulations. This will generate a file called "Free_en_outs.out" in each perturbable system folder, as well as adding to the "results.csv" file in the "energy_results" folder.
+* First run ``run_generate.sh``, this will generate smiles strings for all molecules in the ``ligands`` folder, placing a file named ``SMILES.csv`` in the ``python`` folder. It will then generate perturbable systems in both vacuum and solvent.
+* A folder named ``Systems`` will be created, containing all of the pertubable systems - this is where all simulation data will be stored. Currently 3 replicas of each system are generated.
+* One system generation is complete, run ``min_eq_all_systems.sh``, this will loop over every perturbable system in ``Systems``, minimising the vacuum systems, and minimising and equilibrating the solvated systems. All minimisation and equilibration is perfromed using vanilla OpenMM.
+* Once Minimisation and equilibration are complete, simulations can be run using either ``SOMD1`` or ``SOMD2``. In order to run simulations using ``SOMD1`` run ``run_all_somd1.sh``. This will run all vacuum and solvated systems for all replicas using the SOMD1 protocol defined in ``mineq_solv.py`` and ``mineq_vac.py`` (expect this to take a while). In order to run simulations with ``SOMD2`` run ``run_all_somd2.sh``, this will run all vacuum and solvated legs using the protol defined in ``run_one_somd2.sh``.
+* Once all simulations are complete, run ``analyse_all_somd1.sh`` if SOMD1 was used to run simulations, or ``analyse_all_somd2`` if SOMD2 was used. This will generate a file called ``Free_en_outs.out`` in each perturbable system folder, as well as adding to the ``results.csv`` file in the ``energy_results`` folder.
+
+## Notes
+* Ensure that, prior to running scripts, the proper conda environment is listed in all bash scripts (``conda activate <your env name>``).
+* Protocols for SOMD1 and SOMD2 are defined differently: to alter SOMD1 protocols the relevant ``mineq_*.py`` script should be altered, to alter SOMD2 protocols, the ``run_one_somd2.sh`` script sohuld be altered.
diff --git a/Hydration_freenrg/bash_scripts/analyse_all.sh b/Hydration_freenrg/bash_scripts/analyse_all_somd1.sh
similarity index 92%
rename from Hydration_freenrg/bash_scripts/analyse_all.sh
rename to Hydration_freenrg/bash_scripts/analyse_all_somd1.sh
index 1d2ba4e..82f370d 100755
--- a/Hydration_freenrg/bash_scripts/analyse_all.sh
+++ b/Hydration_freenrg/bash_scripts/analyse_all_somd1.sh
@@ -22,7 +22,7 @@ for filename in ./Systems/*; do
         echo "System directory = $PWD"
 	systemfolder="$PWD"
 	for rep in ./rep*; do
-		cd $rep/solv
+		cd $rep/solvated_somd1
 		echo "$PWD"
 		cd lambda_0.0000
 		lj-tailcorrection -t somd.prm7 -c somd.rst7 -m somd.pert -C somd.cfg -l 0.00 -r traj000000001.dcd -s 1 1> ../freenrg-LJCOR-lam-0.000.dat 2> /dev/null
@@ -38,6 +38,6 @@ for filename in ./Systems/*; do
 	done
 	cd $maindir$filename
 	echo "HERE $PWD"
-	python $maindir/python/analysis.py > hydration_energy.out
+	python $maindir/python/analysis_somd1.py > hydration_energy.out
 	cd $maindir
 done
diff --git a/Hydration_freenrg/bash_scripts/analyse_all_somd2.sh b/Hydration_freenrg/bash_scripts/analyse_all_somd2.sh
new file mode 100755
index 0000000..872a381
--- /dev/null
+++ b/Hydration_freenrg/bash_scripts/analyse_all_somd2.sh
@@ -0,0 +1,27 @@
+#!/bin/bash
+#SBATCH -n 1
+#SBATCH --cpus-per-task 4
+    # allocate 5 GB of memory per job
+#SBATCH --mem 5000
+    # Set job name
+#SBATCH --job-name=analysis
+#SBATCH --output=./scriptouts/analysis.o
+#SBATCH --error=./scriptouts/analysis.err
+
+export CUDA_VISIBLE_DEVICES=""
+eval "$(conda shell.bash hook)"
+conda activate sireDEV
+python3 --version
+cd ../
+maindir="$PWD"
+echo "master dir is $maindir"
+for filename in "$maindir"/Systems/*; do
+        cd $filename
+        echo "System directory = $PWD"
+	systemfolder="$PWD"
+	cd $maindir$filename
+	echo "HERE $PWD"
+	python $maindir/python/analysis.py > hydration_energy.out
+	cd $maindir
+done
+
diff --git a/Hydration_freenrg/bash_scripts/min_eq_all_systems.sh b/Hydration_freenrg/bash_scripts/min_eq_all_systems.sh
index bf32520..0ac82b6 100755
--- a/Hydration_freenrg/bash_scripts/min_eq_all_systems.sh
+++ b/Hydration_freenrg/bash_scripts/min_eq_all_systems.sh
@@ -12,7 +12,7 @@
 #SBATCH --error=./scriptouts/MinEq.err
 
 eval "$(conda shell.bash hook)"
-conda activate openbiosim
+conda activate sireDEV
 python --version
 cd ../
 sdir="$PWD"
diff --git a/Hydration_freenrg/bash_scripts/min_eq_one_system.sh b/Hydration_freenrg/bash_scripts/min_eq_one_system.sh
index 2d565f0..36855b4 100755
--- a/Hydration_freenrg/bash_scripts/min_eq_one_system.sh
+++ b/Hydration_freenrg/bash_scripts/min_eq_one_system.sh
@@ -14,7 +14,7 @@
 #Designed to be run from within the /Systems/{currentsystem} folder
 #Directory work should be handled by MinEqallsystems.sh
 eval "$(conda shell.bash hook)"
-conda activate openbiosim
+conda activate sireDEV
 python3 --version
 reps=( rep0 rep1 rep2 )
 rep=${reps[SLURM_ARRAY_TASK_ID]}
diff --git a/Hydration_freenrg/bash_scripts/run_all_systems.sh b/Hydration_freenrg/bash_scripts/run_all_somd1.sh
similarity index 76%
rename from Hydration_freenrg/bash_scripts/run_all_systems.sh
rename to Hydration_freenrg/bash_scripts/run_all_somd1.sh
index 5d9b9ff..45f36f8 100755
--- a/Hydration_freenrg/bash_scripts/run_all_systems.sh
+++ b/Hydration_freenrg/bash_scripts/run_all_somd1.sh
@@ -12,7 +12,7 @@
 #SBATCH --error=./scriptouts/runall.err
 
 eval "$(conda shell.bash hook)"
-conda activate openbiosim
+conda activate sireDEV
 python3 --version
 cd ../
 maindir="$PWD"
@@ -22,8 +22,8 @@ for filename in ./Systems/*; do
 	for rep in ./rep*; do
 		cd $rep
 		echo "rep directory = $PWD"
-		sbatch --wait --array=0-16 $maindir/bash_scripts/run_one_system_solv.sh
-		sbatch --wait --array=0-16 $maindir/bash_scripts/run_one_system_vac.sh
+		sbatch --wait --array=0-16 $maindir/bash_scripts/run_one_solv_somd1.sh
+		sbatch --wait --array=0-16 $maindir/bash_scripts/run_one_vac_somd1.sh
 		cd ..
 	done
 	cd $maindir
diff --git a/Hydration_freenrg/bash_scripts/run_all_somd2.sh b/Hydration_freenrg/bash_scripts/run_all_somd2.sh
new file mode 100755
index 0000000..3fd5002
--- /dev/null
+++ b/Hydration_freenrg/bash_scripts/run_all_somd2.sh
@@ -0,0 +1,28 @@
+#!/bin/bash
+#SBATCH -n 1
+    # allocate 1 processor as a manager
+#SBATCH --cpus-per-task 1
+    # allocate 0.1 GB of memory per job
+#SBATCH --mem 100
+    # Set job name
+#SBATCH --job-name=runall
+#SBATCH --output=./scriptouts/runall_1fs.o
+#SBATCH --error=./scriptouts/runall_1fs.err
+
+eval "$(conda shell.bash hook)"
+conda activate sireDEV
+python3 --version
+cd ../
+maindir="$PWD"
+for filename in ./Systems/*; do
+        cd $filename
+        echo "System directory = $PWD"
+        for rep in ./rep*; do
+                cd $rep
+                echo "rep directory = $PWD"
+                sbatch --wait $maindir/bash_scripts/run_one_somd2.sh
+                cd ..
+        done
+        cd $maindir
+done
+
diff --git a/Hydration_freenrg/bash_scripts/run_generate.sh b/Hydration_freenrg/bash_scripts/run_generate.sh
index fbe9dc8..f1d1ceb 100755
--- a/Hydration_freenrg/bash_scripts/run_generate.sh
+++ b/Hydration_freenrg/bash_scripts/run_generate.sh
@@ -12,7 +12,7 @@
 #SBATCH --error=./scriptouts/systemgenerate.err
 
 eval "$(conda shell.bash hook)"
-conda activate openbiosim
+conda activate sireDEV
 python3 --version
 ./convtoSMILES.sh
 cd ../
diff --git a/Hydration_freenrg/bash_scripts/run_one_system_solv.sh b/Hydration_freenrg/bash_scripts/run_one_solv_somd1.sh
similarity index 95%
rename from Hydration_freenrg/bash_scripts/run_one_system_solv.sh
rename to Hydration_freenrg/bash_scripts/run_one_solv_somd1.sh
index efbb569..250fc3e 100755
--- a/Hydration_freenrg/bash_scripts/run_one_system_solv.sh
+++ b/Hydration_freenrg/bash_scripts/run_one_solv_somd1.sh
@@ -20,7 +20,7 @@ lam=${lamvals[SLURM_ARRAY_TASK_ID]}
 
 echo "lambda is: " $lam
 
-cd ./solv/lambda_$lam
+cd ./solvated_somd1/lambda_$lam
 echo "$PWD"
 srun somd-freenrg -C ./somd.cfg -c ./somd.rst7 -t ./somd.prm7 -m ./somd.pert -p CUDA -l $lam
 cd ..
diff --git a/Hydration_freenrg/bash_scripts/run_one_somd2.sh b/Hydration_freenrg/bash_scripts/run_one_somd2.sh
new file mode 100755
index 0000000..8fffab4
--- /dev/null
+++ b/Hydration_freenrg/bash_scripts/run_one_somd2.sh
@@ -0,0 +1,24 @@
+#!/bin/bash
+#SBATCH -n 1
+    # Request 1 gpu per job
+#SBATCH --gres=gpu:4
+    # allocate 5 processors/CPUs per GPU
+#SBATCH --cpus-per-gpu 4
+    # allocate 5 GB of memory per job
+#SBATCH --mem 5000
+    # Set job name
+#SBATCH --job-name=solv_run
+#SBATCH --output=run.%j.o
+#SBATCH --error=run.%j.er
+
+eval "$(conda shell.bash hook)"
+conda activate sireDEV
+python3 --version
+file_name_solv="mineq_solv.bss"
+file_name_vac="mineq_vac.bss"
+# Solvated leg
+somd2 ${file_name_solv} --runtime 500ps --timestep 1fs --h-mass-factor 1.0 --num-lambda 17 --temperature 298K --pressure 1atm --checkpoint-frequency 200ns --energy-frequency 0.25ps --output-directory solvated_somd2 --no-restart --cutoff-type rf --cutoff 10.0A --overwrite --perturbable-constraint none --constraint none --somd1-compatibility
+wait
+#Vacuum leg
+somd2 ${file_name_vac} --runtime 500ps --timestep 1fs --h-mass-factor 1.0 --num-lambda 17 --temperature 298K --checkpoint-frequency 1ns --output-directory vacuum_somd2 --energy-frequency 0.25ps --no-restart --cutoff-type rf --cutoff 10.0A --overwrite --perturbable-constraint none --constraint none --somd1-compatibility
+wait
diff --git a/Hydration_freenrg/bash_scripts/run_one_system_vac.sh b/Hydration_freenrg/bash_scripts/run_one_vac_somd1.sh
similarity index 95%
rename from Hydration_freenrg/bash_scripts/run_one_system_vac.sh
rename to Hydration_freenrg/bash_scripts/run_one_vac_somd1.sh
index cedbd84..de98551 100755
--- a/Hydration_freenrg/bash_scripts/run_one_system_vac.sh
+++ b/Hydration_freenrg/bash_scripts/run_one_vac_somd1.sh
@@ -20,7 +20,7 @@ lam=${lamvals[SLURM_ARRAY_TASK_ID]}
 
 echo "lambda is: " $lam
 
-cd ./vac/lambda_$lam
+cd ./vacuum_somd1/lambda_$lam
 echo "$PWD"
 srun somd-freenrg -C ./somd.cfg -c ./somd.rst7 -t ./somd.prm7 -m ./somd.pert -p CUDA -l $lam
 cd ..
diff --git a/Hydration_freenrg/bash_scripts/test.sh b/Hydration_freenrg/bash_scripts/test.sh
deleted file mode 100755
index 859d0aa..0000000
--- a/Hydration_freenrg/bash_scripts/test.sh
+++ /dev/null
@@ -1,23 +0,0 @@
-#!/bin/bash
-#SBATCH -n 1
-    # Request 1 gpu per job
-#SBATCH --gres=gpu:1
-    # allocate 5 processors/CPUs per GPU
-#SBATCH --cpus-per-gpu 5
-    # allocate 5 GB of memory per job
-#SBATCH --mem 5000
-    # Set job name
-#SBATCH --job-name=MinEQ
-#SBATCH --output=MinEq.o
-#SBATCH --error=MinEq.err
-
-source /home/matthew/mambaforge/bin/activate openbiosim
-python --version
-cd ../
-sdir="$PWD"
-for filename in ./Systems/*; do
-	cd $filename
-	echo "before running $PWD"
-	cd ../../
-	echo "after runnning $PWD"
-done
diff --git a/Hydration_freenrg/python/analysis.py b/Hydration_freenrg/python/analysis.py
deleted file mode 100644
index c9b2dbe..0000000
--- a/Hydration_freenrg/python/analysis.py
+++ /dev/null
@@ -1,88 +0,0 @@
-"""Python script designed to find free energies, with error found across multiple runs
-Designed to be run in the folder that contains the rep folders (e.g. Systems/Ethane_Methane)"""
-import BioSimSpace as BSS
-from BioSimSpace._Exceptions import *
-import pandas as pd
-import os
-import numpy as np
-import re
-import math
-from datetime import date
-
-pert_file = os.getcwd()
-list_of_dirs = pert_file.split("/")
-pert = list_of_dirs[-1].split("_")
-
-list_of_reps = []
-for root, dir, files in os.walk("./"):
-	for di in dir:
-		if di.startswith("rep"):
-			list_of_reps.append(di)
-
-results_file = os.path.join("../../","energy_results/results.csv")
-out = open("Free_en_outs.out","w")
-free_ens = []
-errs = []
-s_r=0 #counts the number of successful replicas
-for rep in list_of_reps:
-	if not os.path.isdir(os.path.join(rep,"solv")):
-        	out.write(rep + " failed due to lack of solv directory - likely that Min/Eq failed \n")
-        	continue
-	if not os.path.isdir(os.path.join(rep,"vac")):
-		out.write(rep + " failed due to lack of vac directory - likely that Min/Eq failed \n")
-		continue
-	try:
-		pmf_solv, overlap_solv = BSS.FreeEnergy.Relative.analyse(os.path.join("./",rep,"solv"))
-	except AnalysisError:
-		out.write("\n Free energy analysis failed on " + rep + " solv, likely that one or more lambda windows crashed") 
-		continue
-	print(os.path.join("./",rep,"solv"))
-	files = [os.path.join("./",rep,"solv","freenrg-LJCOR-lam-0.000.dat"),os.path.join("./",rep,"solv","freenrg-LJCOR-lam-1.000.dat")]
-	print(os.path.join("./",rep,"solv","freenrg-LJCOR-lam-0.000.dat"))
-	corrs = []
-	#gets the energies in order [0,1], needs to output 1-0
-	for file in files:
-		print(file)
-		with open(file,"r") as f:
-			lines = f.read().splitlines()
-			last_line = lines[-1]
-			ll = last_line
-			ls = re.findall(r"[-+]?(?:\d*\.*\d+)",ll)
-			corrs.append(float(ls[0]))
-	lj_correction=corrs[1]-corrs[0]
-	try:
-		pmf_vac, overlap_vac = BSS.FreeEnergy.Relative.analyse(os.path.join("./",rep,"vac"))
-	except AnalysisError:
-		out.write("\n Free energy analysis failed on " + rep + " vac, likely that one or more lambda windows crashed")
-		continue
-	free_en = BSS.FreeEnergy.Relative.difference(pmf_solv,pmf_vac)
-	free_ens.append(free_en[0].value()+lj_correction)
-	errs.append(free_en[1].value())
-	out.write("\n" + rep)
-	out.write("\n Base Free Energy = " + str(free_en))
-	out.write("\n Lennard-jones Correction = " + str(lj_correction))
-	out.write("\n Overall Hydration Free Energy for this run = " + str(free_en[0].value()+lj_correction))
-	s_r += 1
-avg = np.nanmean([f for f in free_ens])
-err = math.sqrt(np.sum([f**2 for f in errs]))
-if not np.isnan(avg):
-	out.write("\n Average Free Energy = " + str(avg))
-if err:
-	out.write("\n Error = " + str(err))
-out.close()
-energies_dict = {} #Data to append to csv containing energies
-#titles = ["mol0","mol1","free_energy","error","successful_reps","date"]
-energies_dict["mol0"] = [pert[0]]
-energies_dict["mol1"] = [pert[1]]
-if not np.isnan(avg):
-	energies_dict["free_energy"] = [avg]
-else:
-	energies_dict["free_energy"] = [np.nan]
-if not np.isnan(err):
-        energies_dict["error"] = [err]
-else:
-        energies_dict["error"] = [np.nan]
-energies_dict["successful_reps"] = [s_r]
-energies_dict["date"] = [date.today()]
-df=pd.DataFrame.from_dict(energies_dict)
-df.to_csv(results_file,mode='a',header=False)
diff --git a/Hydration_freenrg/python/analysis_somd1.py b/Hydration_freenrg/python/analysis_somd1.py
new file mode 100644
index 0000000..d209e90
--- /dev/null
+++ b/Hydration_freenrg/python/analysis_somd1.py
@@ -0,0 +1,110 @@
+"""Python script designed to find free energies, with error found across multiple runs
+Designed to be run in the folder that contains the rep folders (e.g. Systems/Ethane_Methane)"""
+import BioSimSpace as BSS
+from BioSimSpace._Exceptions import *
+import pandas as pd
+import os
+import numpy as np
+import re
+import math
+from datetime import date
+
+pert_file = os.getcwd()
+list_of_dirs = pert_file.split("/")
+pert = list_of_dirs[-1].split("_")
+
+list_of_reps = []
+for root, dir, files in os.walk("./"):
+    for di in dir:
+        if di.startswith("rep"):
+            list_of_reps.append(di)
+
+results_file = os.path.join("../../", "energy_results/results.csv")
+out = open("Free_en_outs.out", "w")
+free_ens = []
+errs = []
+s_r = 0  # counts the number of successful replicas
+for rep in list_of_reps:
+    if not os.path.isdir(os.path.join(rep, "solvated_somd1")):
+        out.write(
+            rep + " failed due to lack of solv directory - likely that Min/Eq failed \n"
+        )
+        continue
+    if not os.path.isdir(os.path.join(rep, "vacuum_somd1")):
+        out.write(
+            rep + " failed due to lack of vac directory - likely that Min/Eq failed \n"
+        )
+        continue
+    try:
+        pmf_solv, overlap_solv = BSS.FreeEnergy.Relative.analyse(
+            os.path.join("./", rep, "solvated_somd1")
+        )
+    except AnalysisError:
+        out.write(
+            "\n Free energy analysis failed on "
+            + rep
+            + " solv, likely that one or more lambda windows crashed"
+        )
+        continue
+    print(os.path.join("./", rep, "solvated_somd1"))
+    files = [
+        os.path.join("./", rep, "solvated_somd1", "freenrg-LJCOR-lam-0.000.dat"),
+        os.path.join("./", rep, "solvated_somd1", "freenrg-LJCOR-lam-1.000.dat"),
+    ]
+    print(os.path.join("./", rep, "solvated_somd1", "freenrg-LJCOR-lam-0.000.dat"))
+    corrs = []
+    # gets the energies in order [0,1], needs to output 1-0
+    for file in files:
+        print(file)
+        with open(file, "r") as f:
+            lines = f.read().splitlines()
+            last_line = lines[-1]
+            ll = last_line
+            ls = re.findall(r"[-+]?(?:\d*\.*\d+)", ll)
+            corrs.append(float(ls[0]))
+    lj_correction = corrs[1] - corrs[0]
+    try:
+        pmf_vac, overlap_vac = BSS.FreeEnergy.Relative.analyse(
+            os.path.join("./", rep, "vacuum_somd1")
+        )
+    except AnalysisError:
+        out.write(
+            "\n Free energy analysis failed on "
+            + rep
+            + " vac, likely that one or more lambda windows crashed"
+        )
+        continue
+    free_en = BSS.FreeEnergy.Relative.difference(pmf_solv, pmf_vac)
+    free_ens.append(free_en[0].value() + lj_correction)
+    errs.append(free_en[1].value())
+    out.write("\n" + rep)
+    out.write("\n Base Free Energy = " + str(free_en))
+    out.write("\n Lennard-jones Correction = " + str(lj_correction))
+    out.write(
+        "\n Overall Hydration Free Energy for this run = "
+        + str(free_en[0].value() + lj_correction)
+    )
+    s_r += 1
+avg = np.nanmean([f for f in free_ens])
+err = math.sqrt(np.sum([f**2 for f in errs]))
+if not np.isnan(avg):
+    out.write("\n Average Free Energy = " + str(avg))
+if err:
+    out.write("\n Error = " + str(err))
+out.close()
+energies_dict = {}  # Data to append to csv containing energies
+# titles = ["mol0","mol1","free_energy","error","successful_reps","date"]
+energies_dict["mol0"] = [pert[0]]
+energies_dict["mol1"] = [pert[1]]
+if not np.isnan(avg):
+    energies_dict["free_energy"] = [avg]
+else:
+    energies_dict["free_energy"] = [np.nan]
+if not np.isnan(err):
+    energies_dict["error"] = [err]
+else:
+    energies_dict["error"] = [np.nan]
+energies_dict["successful_reps"] = [s_r]
+energies_dict["date"] = [date.today()]
+df = pd.DataFrame.from_dict(energies_dict)
+df.to_csv(results_file, mode="a", header=False)
diff --git a/Hydration_freenrg/python/analysis_somd2.py b/Hydration_freenrg/python/analysis_somd2.py
new file mode 100644
index 0000000..e9212f2
--- /dev/null
+++ b/Hydration_freenrg/python/analysis_somd2.py
@@ -0,0 +1,146 @@
+"""Python script designed to find free energies, with error found across multiple runs
+Designed to be run in the folder that contains the rep folders (e.g. Systems/Ethane_Methane)"""
+import BioSimSpace as BSS
+from BioSimSpace._Exceptions import *
+import pandas as pd
+import os
+import numpy as np
+import re
+import math
+from datetime import date
+import signal, time
+
+class Timeout():
+  """Timeout class using ALARM signal"""
+  class Timeout(Exception): pass
+
+  def __init__(self, sec):
+    self.sec = sec
+
+  def __enter__(self):
+    signal.signal(signal.SIGALRM, self.raise_timeout)
+    signal.alarm(self.sec)
+
+  def __exit__(self, *args):
+    signal.alarm(0) # disable alarm
+
+  def raise_timeout(self, *args):
+    raise Timeout.Timeout()
+
+pert_file = os.getcwd()
+print(f"LOCAtiON {pert_file}")
+list_of_dirs = pert_file.split("/")
+pert = list_of_dirs[-1].split("_")
+
+list_of_reps = []
+for root, dir, files in os.walk("./"):
+    for di in dir:
+        if di.startswith("rep"):
+            list_of_reps.append(di)
+
+print(list_of_reps)
+results_file = os.path.join("../../", "energy_results/results_1fs_nopert.csv")
+out = open("Free_en_outs.out", "w")
+free_ens = []
+errs = []
+ens_vac = []
+err_vac = []
+ens_solv = []
+err_solv = []
+s_r = 0  # counts the number of successful replicas
+for rep in list_of_reps:
+    file = os.path.join(rep, "solvated_prod_nopertrest")
+    if not os.path.isdir(os.path.join(rep, "solvated_prod_nopertrest")):
+        out.write(
+            f"{rep} failed due to lack of solv directory ({file})- likely that Min/Eq failed \n"
+        )
+        continue
+    if not os.path.isdir(os.path.join(rep, "vacuum_prod_nopertrest")):
+        out.write(
+            rep + " failed due to lack of vac directory - likely that Min/Eq failed \n"
+        )
+        continue
+    try:
+        with Timeout(60):
+            pmf_solv, overlap_solv = BSS.FreeEnergy.Relative.analyse(
+                os.path.join("./", rep, "solvated_prod_nopertrest"), estimator="MBAR"
+            )
+        print(BSS.FreeEnergy.Relative.difference(pmf_solv)[0].value())
+        ens_solv.append(BSS.FreeEnergy.Relative.difference(pmf_solv)[0].value())
+        err_solv.append(BSS.FreeEnergy.Relative.difference(pmf_solv)[1].value())
+    except:
+        out.write(
+            "\n Free energy analysis failed on "
+            + rep
+            + " solv, likely that one or more lambda windows crashed"
+        )
+        continue
+    try:
+        with Timeout(60):
+            pmf_vac, overlap_vac = BSS.FreeEnergy.Relative.analyse(
+                os.path.join("./", rep, "vacuum_prod_nopertrest"), estimator="MBAR"
+            )
+        print(BSS.FreeEnergy.Relative.difference(pmf_vac)[0].value())
+        ens_vac.append(BSS.FreeEnergy.Relative.difference(pmf_vac)[0].value())
+        err_vac.append(BSS.FreeEnergy.Relative.difference(pmf_vac)[1].value())
+    except Exception as e:
+        out.write(
+            "\n Free energy analysis failed on "
+            + rep
+            + " vac, likely that one or more lambda windows crashed"
+        )
+        print(e)
+        continue
+    free_en = BSS.FreeEnergy.Relative.difference(pmf_solv, pmf_vac)
+    free_ens.append(free_en[0].value())
+    errs.append(free_en[1].value())
+    out.write("\n" + rep)
+    out.write(
+        "\n Overall Hydration Free Energy for this run = " + str(free_en[0].value())
+    )
+    s_r += 1
+avg = np.nanmean([f for f in free_ens])
+err = math.sqrt(np.sum([f**2 for f in errs]))
+avg_solv = np.nanmean([f for f in ens_solv])
+err_solv = math.sqrt(np.sum([f**2 for f in err_solv]))
+avg_vac = np.nanmean([f for f in ens_vac])
+err_vac = math.sqrt(np.sum([f**2 for f in err_vac]))
+if not np.isnan(avg):
+    out.write("\n Average Free Energy = " + str(avg))
+if err:
+    out.write("\n Error = " + str(err))
+out.close()
+energies_dict = {}  # Data to append to csv containing energies
+# titles = ["mol0","mol1","free_energy","error","successful_reps","date"]
+energies_dict["mol0"] = [pert[0]]
+energies_dict["mol1"] = [pert[1]]
+if not np.isnan(avg):
+    energies_dict["free_energy"] = [avg]
+else:
+    energies_dict["free_energy"] = [np.nan]
+if not np.isnan(err):
+    energies_dict["error"] = [err]
+else:
+    energies_dict["error"] = [np.nan]
+if not np.isnan(avg_solv):
+    energies_dict["free_energy_solv"] = [avg_solv]
+else:
+    energies_dict["free_energy_solv"] = [np.nan]
+if not np.isnan(err_solv):
+    energies_dict["error_solv"] = [err_solv]
+else:
+    energies_dict["error_solv"] = [np.nan]
+if not np.isnan(avg_vac):
+    energies_dict["free_energy_vac"] = [avg_vac]
+else:
+    energies_dict["free_energy_vac"] = [np.nan]
+if not np.isnan(err_solv):
+    energies_dict["error_vac"] = [err_vac]
+else:
+    energies_dict["error_vac"] = [np.nan]
+
+
+energies_dict["successful_reps"] = [s_r]
+energies_dict["date"] = [date.today()]
+df = pd.DataFrame.from_dict(energies_dict)
+df.to_csv(results_file, mode="a", header=False)
diff --git a/Hydration_freenrg/python/min_eq_solv.py b/Hydration_freenrg/python/min_eq_solv.py
index b67d468..fd78506 100644
--- a/Hydration_freenrg/python/min_eq_solv.py
+++ b/Hydration_freenrg/python/min_eq_solv.py
@@ -7,6 +7,7 @@
 designed to be run in the respective replica folder
 """
 
+
 class MinimisationError(Exception):
     """Exception raised when failure occured during a minimisation step
 
@@ -14,11 +15,13 @@ class MinimisationError(Exception):
         min -- minimisation step that caused failure
         message -- explanation of error
     """
+
     def __init__(self, min, message="Minimisation Failed"):
-        self.min=min
-        self.message=message
+        self.min = min
+        self.message = message
         super().__init__(self.message)
 
+
 class EquilibrationError(Exception):
     """Exception raised when failure occured during an equilibration step
 
@@ -26,11 +29,13 @@ class EquilibrationError(Exception):
         eq -- equilibration step that caused failure
         message -- explanation of error
     """
+
     def __init__(self, eq, message="Equilibration Failed"):
-        self.eq=eq
-        self.message=message
+        self.eq = eq
+        self.message = message
         super().__init__(self.message)
 
+
 class GetSystemError(Exception):
     """Exception raised when a system of value None is found
 
@@ -38,27 +43,29 @@ class GetSystemError(Exception):
         sys -- name of none system
         message -- explanation of error
     """
+
     def __init__(self, sys, message="System with value None found"):
-        self.sys=sys
-        self.message=message
+        self.sys = sys
+        self.message = message
         super().__init__(self.message)
 
-#BSS.verbose(True)
+
+# BSS.verbose(True)
 solvated = BSS.IO.readPerturbableSystem(
     top0="pertsave_solv0.prm7",
     coords0="pertsave_solv0.rst7",
     top1="pertsave_solv1.prm7",
-    coords1= "pertsave_solv1.rst7",
+    coords1="pertsave_solv1.rst7",
 )
-#Need to manually change the device index in each CUDA equilibration to 0
-#stops wierd interactions with SLURM
-#(can be done with CUDA_VISIBLE_DEVICES instead but this is easier for large scripts)
+# Need to manually change the device index in each CUDA equilibration to 0
+# stops wierd interactions with SLURM
+# (can be done with CUDA_VISIBLE_DEVICES instead but this is easier for large scripts)
 substring = "CudaDeviceIndex"
 
 minimise_1 = BSS.Protocol.Minimisation(
     steps=1000, restraint="heavy", force_constant=5.0
 )
-process_m1 = BSS.Process.OpenMM(solvated, minimise_1, work_dir="Minimise_1")
+process_m1 = BSS.Process.Amber(solvated, minimise_1, work_dir="Minimise_1")
 process_m1.start()
 process_m1.wait()
 if process_m1.isError():
@@ -76,9 +83,9 @@ def __init__(self, sys, message="System with value None found"):
     force_constant=5.0,
     thermostat_time_constant=0.5 * BSS.Units.Time.picosecond,
 )
-process_e1 = BSS.Process.OpenMM(min1, equil_1, work_dir="Equilibrate_1",platform='CUDA')
+process_e1 = BSS.Process.Amber(min1, equil_1, work_dir="Equilibrate_1", platform="CUDA")
 cfg1 = process_e1.getConfig()
-#This will break if there are more than 1 CudaDeviceIndex variables in the config, but I can't see this happening
+# This will break if there are more than 1 CudaDeviceIndex variables in the config, but I can't see this happening
 i = [cfg1.index(s) for s in cfg1 if substring in s]
 cfg1[i[0]] = "properties = {'CudaDeviceIndex': '0'}"
 process_e1.setConfig(cfg1)
@@ -94,7 +101,7 @@ def __init__(self, sys, message="System with value None found"):
 minimise_2 = BSS.Protocol.Minimisation(
     steps=1000, restraint="heavy", force_constant=2.0
 )
-process_m2 = BSS.Process.OpenMM(eq1, minimise_2, work_dir="Minimise_2")
+process_m2 = BSS.Process.Amber(eq1, minimise_2, work_dir="Minimise_2")
 process_m2.start()
 process_m2.wait()
 if process_m2.isError():
@@ -107,7 +114,7 @@ def __init__(self, sys, message="System with value None found"):
 minimise_3 = BSS.Protocol.Minimisation(
     steps=1000, restraint="heavy", force_constant=0.1
 )
-process_m3 = BSS.Process.OpenMM(min2, minimise_3, work_dir="Minimise_3")
+process_m3 = BSS.Process.Amber(min2, minimise_3, work_dir="Minimise_3")
 process_m3.start()
 process_m3.wait()
 if process_m3.isError():
@@ -118,7 +125,7 @@ def __init__(self, sys, message="System with value None found"):
 print("Third Minimisation Complete")
 
 minimise_4 = BSS.Protocol.Minimisation(steps=1000)
-process_m4 = BSS.Process.OpenMM(min3, minimise_4, work_dir="Minimise_4")
+process_m4 = BSS.Process.Amber(min3, minimise_4, work_dir="Minimise_4")
 process_m4.start()
 process_m4.wait()
 if process_m4.isError():
@@ -137,7 +144,7 @@ def __init__(self, sys, message="System with value None found"):
     force_constant=1.0,
     thermostat_time_constant=1 * BSS.Units.Time.picosecond,
 )
-process_e2 = BSS.Process.OpenMM(min4, equil_2, work_dir="Equilibrate_2",platform='CUDA')
+process_e2 = BSS.Process.Amber(min4, equil_2, work_dir="Equilibrate_2", platform="CUDA")
 cfg2 = process_e2.getConfig()
 i = [cfg2.index(s) for s in cfg2 if substring in s]
 cfg2[i[0]] = "properties = {'CudaDeviceIndex': '0'}"
@@ -160,7 +167,7 @@ def __init__(self, sys, message="System with value None found"):
     force_constant=0.5,
     thermostat_time_constant=1 * BSS.Units.Time.picosecond,
 )
-process_e3 = BSS.Process.OpenMM(eq2, equil_3, work_dir="Equilibrate_3",platform='CUDA')
+process_e3 = BSS.Process.Amber(eq2, equil_3, work_dir="Equilibrate_3", platform="CUDA")
 cfg3 = process_e3.getConfig()
 i = [cfg3.index(s) for s in cfg3 if substring in s]
 cfg3[i[0]] = "properties = {'CudaDeviceIndex': '0'}"
@@ -179,11 +186,11 @@ def __init__(self, sys, message="System with value None found"):
     runtime=10 * BSS.Units.Time.picosecond,
     temperature=300 * BSS.Units.Temperature.kelvin,
     pressure=1.0 * BSS.Units.Pressure.bar,
-#    restraint="backbone", uncomment for systems contraining protein
-#    force_constant=0.5,  "                                         "
+    #    restraint="backbone", uncomment for systems contraining protein
+    #    force_constant=0.5,  "                                         "
     thermostat_time_constant=1 * BSS.Units.Time.picosecond,
 )
-process_e4 = BSS.Process.OpenMM(eq3, equil_4, work_dir="Equilibrate_4",platform='CUDA')
+process_e4 = BSS.Process.Amber(eq3, equil_4, work_dir="Equilibrate_4", platform="CUDA")
 cfg4 = process_e4.getConfig()
 i = [cfg4.index(s) for s in cfg4 if substring in s]
 cfg4[i[0]] = "properties = {'CudaDeviceIndex': '0'}"
@@ -204,7 +211,7 @@ def __init__(self, sys, message="System with value None found"):
     pressure=1.0 * BSS.Units.Pressure.bar,
     thermostat_time_constant=1.0 * BSS.Units.Time.picosecond,
 )
-process_e5 = BSS.Process.OpenMM(eq4, equil_5, work_dir="Equilibrate_5",platform='CUDA')
+process_e5 = BSS.Process.Amber(eq4, equil_5, work_dir="Equilibrate_5", platform="CUDA")
 cfg5 = process_e5.getConfig()
 i = [cfg5.index(s) for s in cfg5 if substring in s]
 cfg5[i[0]] = "properties = {'CudaDeviceIndex': '0'}"
@@ -218,9 +225,12 @@ def __init__(self, sys, message="System with value None found"):
     raise GetSystemError(eq5)
 print("Fifth Equilibration Complete")
 
+# Save binary for SOMD2
+eq5.save("mineq_solv")
 protocol_run = BSS.Protocol.FreeEnergyProduction(
     num_lam=17,
-    runtime=BSS.Types.Time(2.0, "ns"),
+    timestep=BSS.Types.Time(1.0, "fs"),
+    runtime=BSS.Types.Time(500.0, "ps"),
     report_interval=20000,
     restart_interval=50,
     temperature=BSS.Types.Temperature(298, "K"),
@@ -230,9 +240,8 @@ def __init__(self, sys, message="System with value None found"):
     eq5,
     protocol_run,
     engine="somd",
-    work_dir="solv",
-    extra_options={"minimise": "True","gpu":"0"},
+    work_dir="solvated_somd1",
+    extra_options={"minimise": "True", "gpu": "0", "constraint": "none"},
     setup_only=True,
 )
 print("Success!")
-
diff --git a/Hydration_freenrg/python/min_eq_vac.py b/Hydration_freenrg/python/min_eq_vac.py
index c109af7..28628b6 100644
--- a/Hydration_freenrg/python/min_eq_vac.py
+++ b/Hydration_freenrg/python/min_eq_vac.py
@@ -1,6 +1,8 @@
 """Minimisation of the vacuum systems using the default BioSimSpace minimisation protocol"""
+
 import BioSimSpace as BSS
 
+
 class MinimisationError(Exception):
     """Exception raised when failure occured during a minimisation step
 
@@ -8,11 +10,13 @@ class MinimisationError(Exception):
         min -- minimisation step that caused failure
         message -- explanation of error
     """
+
     def __init__(self, min, message="Minimisation Failed"):
-        self.min=min
-        self.message=message
+        self.min = min
+        self.message = message
         super().__init__(self.message)
 
+
 class GetSystemError(Exception):
     """Exception raised when a system of value None is found
 
@@ -20,20 +24,22 @@ class GetSystemError(Exception):
         sys -- name of none system
         message -- explanation of error
     """
+
     def __init__(self, sys, message="System with value None found"):
-        self.sys=sys
-        self.message=message
+        self.sys = sys
+        self.message = message
         super().__init__(self.message)
 
+
 vac = BSS.IO.readPerturbableSystem(
     top0="pertsave_vac0.prm7",
     coords0="pertsave_vac0.rst7",
     top1="pertsave_vac1.prm7",
-    coords1= "pertsave_vac1.rst7",
+    coords1="pertsave_vac1.rst7",
 )
 
 minimise = BSS.Protocol.Minimisation()
-process_m1 = BSS.Process.OpenMM(vac, minimise, work_dir="Minimise_Vacuum")
+process_m1 = BSS.Process.Amber(vac, minimise, work_dir="Minimise_Vacuum")
 process_m1.start()
 process_m1.wait()
 if process_m1.isError():
@@ -42,10 +48,12 @@ def __init__(self, sys, message="System with value None found"):
 if min1 is None:
     raise GetSystemError(min1)
 print("Minimisation Complete")
-
+# Saving binary for somd2
+min1.save("mineq_vac")
 protocol_run = BSS.Protocol.FreeEnergyProduction(
     num_lam=17,
-    runtime=BSS.Types.Time(2.0, "ns"),
+    runtime=BSS.Types.Time(500.0, "ps"),
+    timestep=BSS.Types.Time(1.0, "fs"),
     report_interval=10000,
     restart_interval=100,
     temperature=BSS.Types.Temperature(298, "K"),
@@ -55,9 +63,8 @@ def __init__(self, sys, message="System with value None found"):
     min1,
     protocol_run,
     engine="somd",
-    work_dir="vac",
-    extra_options={"minimise": "True","gpu":"0"},
+    work_dir="vacuum_somd1",
+    extra_options={"minimise": "True", "gpu": "0", "constraint": "none"},
     setup_only=True,
 )
 print("Success!")
-