diff --git a/examples/HydrostaticHalo/hydrostatic.yml b/examples/HydrostaticHalo/hydrostatic.yml
index ad7b4abe67231db7385d6e7ff2343dd882301434..177f950ea248a8a25701e68c45e6dfe94ecb469a 100644
--- a/examples/HydrostaticHalo/hydrostatic.yml
+++ b/examples/HydrostaticHalo/hydrostatic.yml
@@ -9,7 +9,7 @@ InternalUnitSystem:
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
-  time_end:   0.1    # The end time of the simulation (in internal units).
+  time_end:   1.0    # The end time of the simulation (in internal units).
   dt_min:     1e-6  # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
 
diff --git a/examples/HydrostaticHalo/makeIC.py b/examples/HydrostaticHalo/makeIC.py
index 49bd73f4fa122e2b35e73d323bbba2a95bd3f90f..09e7734b0eb29e20220c8444486028b573e2b9cd 100644
--- a/examples/HydrostaticHalo/makeIC.py
+++ b/examples/HydrostaticHalo/makeIC.py
@@ -73,7 +73,7 @@ const_G                = ((CONST_G_CGS*const_unit_mass_in_cgs*const_unit_time_in
 print 'G=', const_G
 
 # Parameters
-periodic= 0            # 1 For periodic box
+periodic= 1            # 1 For periodic box
 boxSize = 4.          
 G       = const_G 
 N       = int(sys.argv[1])  # Number of particles
@@ -128,6 +128,11 @@ coords += np.full((N,3),boxSize/2.)
 print "x range = (%f,%f)" %(np.min(coords[:,0]),np.max(coords[:,0]))
 print "y range = (%f,%f)" %(np.min(coords[:,1]),np.max(coords[:,1]))
 print "z range = (%f,%f)" %(np.min(coords[:,2]),np.max(coords[:,2]))
+
+print np.mean(coords[:,0])
+print np.mean(coords[:,1])
+print np.mean(coords[:,2])
+
 ds = grp.create_dataset('Coordinates', (N, 3), 'd')
 ds[()] = coords
 coords = np.zeros(1)
diff --git a/examples/HydrostaticHalo/run.sh b/examples/HydrostaticHalo/run.sh
index 7b462f447b5065e2a7da5c660dabbfe0ac37e634..b37c3b4532316953eeb3a370d72ede78063240f3 100755
--- a/examples/HydrostaticHalo/run.sh
+++ b/examples/HydrostaticHalo/run.sh
@@ -9,3 +9,5 @@ python makeIC.py 1000
 python radial_profile.py 10
 
 python internal_energy_profile.py 10
+
+python test_energy_conservation.py
diff --git a/examples/main.c b/examples/main.c
index 895bb9003e04a4723b13ea1eaa360c5aa44e7228..d7080ee1f1a97b0b056c14d4f431596fa8441f27 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -435,7 +435,7 @@ int main(int argc, char *argv[]) {
   /* Initialise the external potential properties */
   struct external_potential potential;
   if (with_external_gravity)
-    potential_init(params, &prog_const, &us, &potential);
+    potential_init(params, &prog_const, &us, &s, &potential);
   if (with_external_gravity && myrank == 0) potential_print(&potential);
 
   /* Initialise the cooling function properties */
diff --git a/src/const.h b/src/const.h
index d905f636dc0f0fcd97560b516afdfa29771e87df..1855d04de763f3d458ea3362d3eb4e6d153c5d54 100644
--- a/src/const.h
+++ b/src/const.h
@@ -92,9 +92,9 @@
 
 /* External gravity properties */
 
-#define EXTERNAL_POTENTIAL_NONE
+//#define EXTERNAL_POTENTIAL_NONE
 //#define EXTERNAL_POTENTIAL_POINTMASS
-//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
 //#define EXTERNAL_POTENTIAL_DISC_PATCH
 
 
diff --git a/src/potential.c b/src/potential.c
index 5433a05e3e7886ad88021d3916cae26adfe8b954..534540ea5bc4e7d05bbd91761548411962f2c327 100644
--- a/src/potential.c
+++ b/src/potential.c
@@ -36,9 +36,10 @@
 void potential_init(const struct swift_params* parameter_file,
                     const struct phys_const* phys_const,
                     const struct UnitSystem* us,
+		    const struct space* s,
                     struct external_potential* potential) {
 
-  potential_init_backend(parameter_file, phys_const, us, potential);
+  potential_init_backend(parameter_file, phys_const, us, s, potential);
 }
 
 /**
diff --git a/src/potential.h b/src/potential.h
index 77bd41794a3a8cd244405493898d63b3f80ff3a6..de23882770c4473aabbc41da1c4f18e5d13a0c03 100644
--- a/src/potential.h
+++ b/src/potential.h
@@ -47,6 +47,7 @@
 void potential_init(const struct swift_params* parameter_file,
                     const struct phys_const* phys_const,
                     const struct UnitSystem* us,
+		    const struct space* s,
                     struct external_potential* potential);
 
 void potential_print(const struct external_potential* potential);
diff --git a/src/potential/disc_patch/potential.h b/src/potential/disc_patch/potential.h
index 21d168818e164ad3b3e18076ba824285e40956aa..116b9dab36fef1c4acf816fb929b400620851b53 100644
--- a/src/potential/disc_patch/potential.h
+++ b/src/potential/disc_patch/potential.h
@@ -34,6 +34,7 @@
 #include "part.h"
 #include "physical_constants.h"
 #include "units.h"
+#include "space.h"
 
 /**
  * @brief External Potential Properties - Disc patch case
@@ -161,6 +162,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
 static INLINE void potential_init_backend(
     const struct swift_params* parameter_file,
     const struct phys_const* phys_const, const struct UnitSystem* us,
+    const struct space* s,
     struct external_potential* potential) {
 
   potential->surface_density = parser_get_param_double(
diff --git a/src/potential/isothermal/potential.h b/src/potential/isothermal/potential.h
index a993c09a978ca3692ec3359f7633df14760f263d..b1f0ef95b27d5afcebe6759e7e30826292eed4ed 100644
--- a/src/potential/isothermal/potential.h
+++ b/src/potential/isothermal/potential.h
@@ -32,6 +32,7 @@
 #include "part.h"
 #include "physical_constants.h"
 #include "units.h"
+#include "space.h"
 
 /**
  * @brief External Potential Properties - Isothermal sphere case
@@ -126,13 +127,14 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
 static INLINE void potential_init_backend(
     const struct swift_params* parameter_file,
     const struct phys_const* phys_const, const struct UnitSystem* us,
+    const struct space* s,
     struct external_potential* potential) {
 
-  potential->x =
+  potential->x = s->dim[0]/2. + 
       parser_get_param_double(parameter_file, "IsothermalPotential:position_x");
-  potential->y =
+  potential->y = s->dim[1]/2. +
       parser_get_param_double(parameter_file, "IsothermalPotential:position_y");
-  potential->z =
+  potential->z = s->dim[2]/2. +
       parser_get_param_double(parameter_file, "IsothermalPotential:position_z");
   potential->vrot =
       parser_get_param_double(parameter_file, "IsothermalPotential:vrot");
diff --git a/src/potential/none/potential.h b/src/potential/none/potential.h
index 8b1c3e841521f3fb42fbdf5c8922cead2ea7cbcb..3f075f8fdc9b10e107bb4663b26d10e37119c9bb 100644
--- a/src/potential/none/potential.h
+++ b/src/potential/none/potential.h
@@ -31,6 +31,7 @@
 #include "part.h"
 #include "physical_constants.h"
 #include "units.h"
+#include "space.h"
 
 /**
  * @brief External Potential Properties
@@ -81,6 +82,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
 static INLINE void potential_init_backend(
     const struct swift_params* parameter_file,
     const struct phys_const* phys_const, const struct UnitSystem* us,
+    const struct space* s,
     struct external_potential* potential) {}
 
 /**
diff --git a/src/potential/point_mass/potential.h b/src/potential/point_mass/potential.h
index f718af2e2c4ff91540e1834cb2072d321ce38705..28f9e81b43b2aca8954c5343397af941f1620ff0 100644
--- a/src/potential/point_mass/potential.h
+++ b/src/potential/point_mass/potential.h
@@ -32,6 +32,7 @@
 #include "part.h"
 #include "physical_constants.h"
 #include "units.h"
+#include "space.h"
 
 /**
  * @brief External Potential Properties - Point mass case
@@ -127,6 +128,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
 static INLINE void potential_init_backend(
     const struct swift_params* parameter_file,
     const struct phys_const* phys_const, const struct UnitSystem* us,
+    const struct space* s,
     struct external_potential* potential) {
 
   potential->x =