diff --git a/configure.ac b/configure.ac
index a74df0fc5281af08816d2bb50f6f60e2f501fe2f..a8facd73fc9132081c90a5ecb640c4e877e1bd6a 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1340,6 +1340,7 @@ with_subgrid_cooling=none
 with_subgrid_chemistry=none
 with_subgrid_tracers=none
 with_subgrid_entropy_floor=none
+with_subgrid_pressure_floor=none
 with_subgrid_stars=none
 with_subgrid_star_formation=none
 with_subgrid_feedback=none
@@ -1351,10 +1352,9 @@ case "$with_subgrid" in
    none)
    ;;
    GEAR)
-	with_subgrid_cooling=grackle
+	with_subgrid_cooling=grackle_0
 	with_subgrid_chemistry=GEAR
-	with_subgrid_tracers=none
-	with_subgrid_entropy_floor=none
+	with_subgrid_pressure_floor=GEAR
 	with_subgrid_stars=GEAR
 	with_subgrid_star_formation=GEAR
 	with_subgrid_feedback=none
@@ -1413,6 +1413,12 @@ AC_ARG_WITH([gravity],
    [with_gravity="default"]
 )
 
+if test "$with_subgrid" = "EAGLE"; then
+   if test "$with_gravity" = "default"; then
+      with_gravity="with-potential"
+   fi
+fi
+
 case "$with_gravity" in
    with-potential)
       AC_DEFINE([POTENTIAL_GRAVITY], [1], [Gravity scheme with potential calculation])
@@ -1653,7 +1659,8 @@ esac
 #  Cooling function
 AC_ARG_WITH([cooling],
    [AS_HELP_STRING([--with-cooling=<function>],
-      [cooling function @<:@none, const-du, const-lambda, EAGLE, grackle, grackle1, grackle2, grackle3 default: none@:>@]
+      [cooling function @<:@none, const-du, const-lambda, EAGLE, grackle_* default: none@:>@.
+      For Grackle, you need to provide the primordial chemistry parameter (e.g. grackle_0)]
    )],
    [with_cooling="$withval"],
    [with_cooling="none"]
@@ -1680,21 +1687,10 @@ case "$with_cooling" in
    compton)
       AC_DEFINE([COOLING_COMPTON], [1], [Compton cooling off the CMB])
    ;;
-   grackle)
-      AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
-      AC_DEFINE([COOLING_GRACKLE_MODE], [0], [Grackle chemistry network, mode 0])
-   ;;
-   grackle1)
-      AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
-      AC_DEFINE([COOLING_GRACKLE_MODE], [1], [Grackle chemistry network, mode 1])
-   ;;
-   grackle2)
+   grackle_*)
       AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
-      AC_DEFINE([COOLING_GRACKLE_MODE], [2], [Grackle chemistry network, mode 2])
-   ;;
-   grackle3)
-      AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
-      AC_DEFINE([COOLING_GRACKLE_MODE], [3], [Grackle chemistry network, mode 3])
+      primordial_chemistry=${with_cooling:8}
+      AC_DEFINE_UNQUOTED([COOLING_GRACKLE_MODE], [$primordial_chemistry], [Grackle chemistry network])
    ;;
    EAGLE)
       AC_DEFINE([COOLING_EAGLE], [1], [Cooling following the EAGLE model])
@@ -1929,6 +1925,35 @@ case "$with_entropy_floor" in
    ;;
 esac 
 
+#  Pressure floor
+AC_ARG_WITH([pressure-floor], 
+    [AS_HELP_STRING([--with-pressure-floor=<floor>],
+       [pressure floor @<:@none, GEAR, default: none@:>@
+       The hydro model needs to be compatible.]
+    )],
+    [with_pressure_floor="$withval"],
+    [with_pressure_floor="none"]
+)
+if test "$with_subgrid" != "none"; then
+   if test "$with_pressure_floor" != "none"; then
+      AC_MSG_ERROR([Cannot provide with-subgrid and with-pressure-floor together])
+   else
+      with_pressure_floor="$with_subgrid_pressure_floor"
+   fi
+fi
+
+case "$with_pressure_floor" in
+   none)
+      AC_DEFINE([PRESSURE_FLOOR_NONE], [1], [No pressure floor])
+   ;;
+   GEAR)
+      AC_DEFINE([PRESSURE_FLOOR_GEAR], [1], [GEAR pressure floor])
+   ;;
+   *)
+      AC_MSG_ERROR([Unknown pressure floor model])
+   ;;
+esac 
+
 #  Star formation
 AC_ARG_WITH([star-formation], 
     [AS_HELP_STRING([--with-star-formation=<sfm>],
@@ -1974,6 +1999,13 @@ AC_DEFINE_UNQUOTED([SELF_GRAVITY_MULTIPOLE_ORDER], [$with_multipole_order], [Mul
 AC_PATH_PROG([GIT_CMD], [git])
 AC_SUBST([GIT_CMD])
 
+# Check that the gravity model is compatible with the subgrid
+if test $with_black_holes = "EAGLE"; then
+   if test $with_gravity != "with-potential"; then
+      AC_MSG_ERROR([The EAGLE BH model needs the gravity scheme to provide potentials. The code must be compile with --with-gravity=with-potential.])
+   fi
+fi
+
 # Make the documentation. Add conditional to handle disable option.
 DX_INIT_DOXYGEN(libswift,doc/Doxyfile,doc/)
 AM_CONDITIONAL([HAVE_DOXYGEN], [test "$ac_cv_path_ac_pt_DX_DOXYGEN" != ""])
@@ -2056,6 +2088,7 @@ AC_MSG_RESULT([
    External potential  : $with_potential
 
    Entropy floor        : $with_entropy_floor
+   Pressure floor       : $with_pressure_floor
    Cooling function     : $with_cooling
    Chemistry            : $with_chemistry
    Tracers              : $with_tracers
diff --git a/doc/RTD/source/SubgridModels/EAGLE/index.rst b/doc/RTD/source/SubgridModels/EAGLE/index.rst
index 56736cc4cca4bd9042a2c84d07dcf535f48c566c..cac6716c7568b7c7d2168a77586cc80123df7127 100644
--- a/doc/RTD/source/SubgridModels/EAGLE/index.rst
+++ b/doc/RTD/source/SubgridModels/EAGLE/index.rst
@@ -577,6 +577,11 @@ Black-hole creation
 Black-hole accretion
 ~~~~~~~~~~~~~~~~~~~~
 
+.. _EAGLE_black_hole_reposition:
+
+Black-hole repositioning
+~~~~~~~~~~~~~~~~~~~~~~~~
+
 .. _EAGLE_black_hole_feedback:
 
 AGN feedback
diff --git a/doc/RTD/source/SubgridModels/GEAR/index.rst b/doc/RTD/source/SubgridModels/GEAR/index.rst
index 2a211759bfc4895fd07279b72f78200d6ea47546..152056ef2eb60da189bf18b3b3c9909f88ace6de 100644
--- a/doc/RTD/source/SubgridModels/GEAR/index.rst
+++ b/doc/RTD/source/SubgridModels/GEAR/index.rst
@@ -5,6 +5,29 @@
 GEAR model
 ===========
 
+Pressure Floor
+~~~~~~~~~~~~~~
+
+In order to avoid the artificial collapse of unresolved clumps, a minimum in pressure is applied to the particles.
+This additional pressure can be seen as the pressure due to unresolved hydrodynamics turbulence and is given by:
+
+.. math::
+    P_\textrm{Jeans} = \frac{\rho}{\gamma} \left( \frac{4}{\pi} G h^2 \rho N_\textrm{Jeans}^{2/3} - \sigma^2 \right)
+
+where :math:`\rho` is the density, :math:`\gamma` the adiabatic index, :math:`G` is the gravitational constant,
+:math:`h` the smoothing length, :math:`N_\textrm{Jeans}` is the number of particle required in order to resolve a clump and
+:math:`\sigma` the velocity dispersion.
+
+
+This must be directly implemented into the hydro schemes, therefore only a subset of schemes (Gadget-2 and Pressure-Energy) are currently implemented.
+In order to implement it, you need equation 12 in `Hopkins 2013 <https://arxiv.org/abs/1206.5006>`_:
+
+.. math::
+   m_i \frac{\mathrm{d}v_i}{\mathrm{d}t} = - \sum_j x_i x_j \left[ \frac{P_i}{y_i^2} f_{ij} \nabla_i W_{ij}(h_i) + \frac{P_j}{y_j^2} f_{ji} \nabla_j W_{ji}(h_j) \right]
+
+and simply replace the :math:`P_i, P_j` by the pressure with the floor.
+Here the :math:`x, y` are simple weights that should never have the pressure floor included even if they are related to the pressure (e.g. pressure-entropy).
+
 
 Cooling: Grackle
 ~~~~~~~~~~~~~~~~
diff --git a/doc/RTD/source/conf.py b/doc/RTD/source/conf.py
index 8cda2a011e2e12d8d1184224ff77c1739e8aacf8..30780800e7829caccb5c443c4bfd97649479b590 100644
--- a/doc/RTD/source/conf.py
+++ b/doc/RTD/source/conf.py
@@ -25,7 +25,7 @@ author = 'SWIFT Team'
 # The short X.Y version
 version = '0.8'
 # The full version, including alpha/beta/rc tags
-release = '0.8.4
+release = '0.8.4'
 
 
 # -- General configuration ---------------------------------------------------
diff --git a/examples/Cooling/ConstantCosmoTempEvolution/plot_thermal_history.py b/examples/Cooling/ConstantCosmoTempEvolution/plot_thermal_history.py
index 1494102531104b252e3edaa467920db7383ac6e6..1f9ae6557d1a55bd75d3440679ae42d653da3399 100644
--- a/examples/Cooling/ConstantCosmoTempEvolution/plot_thermal_history.py
+++ b/examples/Cooling/ConstantCosmoTempEvolution/plot_thermal_history.py
@@ -68,21 +68,21 @@ for snap in snap_list:
     z = np.append(z, data.metadata.z)
         
     # Convert gas temperatures and densities to right units
-    data.gas.temperature.convert_to_cgs()
+    data.gas.temperatures.convert_to_cgs()
 
     # Get mean and standard deviation of temperature
-    T_mean.append(np.mean(data.gas.temperature) * data.gas.temperature.units)
-    T_std.append(np.std(data.gas.temperature) * data.gas.temperature.units)
+    T_mean.append(np.mean(data.gas.temperatures) * data.gas.temperatures.units)
+    T_std.append(np.std(data.gas.temperatures) * data.gas.temperatures.units)
 
     # Get mean and standard deviation of density
-    rho_mean.append(np.mean(data.gas.density) * data.gas.density.units)
-    rho_std.append(np.std(data.gas.density) * data.gas.density.units)
+    rho_mean.append(np.mean(data.gas.densities) * data.gas.densities.units)
+    rho_std.append(np.std(data.gas.densities) * data.gas.densities.units)
     
 ## Turn into numpy arrays
-T_mean = np.array(T_mean) * data.gas.temperature.units
-T_std = np.array(T_std) * data.gas.temperature.units
-rho_mean = np.array(rho_mean) * data.gas.density.units
-rho_std = np.array(rho_std) * data.gas.density.units
+T_mean = np.array(T_mean) * data.gas.temperatures.units
+T_std = np.array(T_std) * data.gas.temperatures.units
+rho_mean = np.array(rho_mean) * data.gas.densities.units
+rho_std = np.array(rho_std) * data.gas.densities.units
 
 ## Put Density into units of mean baryon density today
 
diff --git a/examples/Cooling/CoolingBox/plotEnergy.py b/examples/Cooling/CoolingBox/plotEnergy.py
index 9c7af57d3d9dfdcfa222e9f77701f230d25f9ddc..6234c319a3001fa4e364639bb2c5480a31cf99dd 100644
--- a/examples/Cooling/CoolingBox/plotEnergy.py
+++ b/examples/Cooling/CoolingBox/plotEnergy.py
@@ -87,7 +87,7 @@ temp_snap = np.zeros(N)
 time_snap_cgs = np.zeros(N)
 for i in range(N):
     snap = File(files[i], 'r')
-    u = snap["/PartType0/InternalEnergy"][:] * snap["/PartType0/Masses"][:]
+    u = snap["/PartType0/InternalEnergies"][:] * snap["/PartType0/Masses"][:]
     u = sum(u) / total_mass[0]
     temp_snap[i] = energyUnits(u)
     time_snap_cgs[i] = snap["/Header"].attrs["Time"] * unit_time
diff --git a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py
index 9cde36cb05b88e60b3bf3527f3857a3774bf5dca..cb624be3cfb1f0f803cfce7a14fb1a772cccf515 100644
--- a/examples/Cooling/CoolingRedshiftDependence/plotSolution.py
+++ b/examples/Cooling/CoolingRedshiftDependence/plotSolution.py
@@ -94,15 +94,17 @@ def get_data(handle: float, n_snaps: int):
                 t0 = data.metadata.t.to(Myr).value
 
             times.append(data.metadata.t.to(Myr).value - t0)
-            temps.append(np.mean(data.gas.temperature.to(K).value))
+            temps.append(np.mean(data.gas.temperatures.to(K).value))
             densities.append(
-                np.mean(data.gas.density.to(mh / cm ** 3).value)
+                np.mean(data.gas.densities.to(mh / cm ** 3).value)
                 / (data.metadata.scale_factor ** 3)
             )
 
             try:
                 energies.append(
-                    np.mean((data.gas.internal_energy * data.gas.masses).to(erg).value)
+                    np.mean(
+                        (data.gas.internal_energies * data.gas.masses).to(erg).value
+                    )
                     * data.metadata.scale_factor ** (2)
                 )
                 radiated_energies.append(
diff --git a/examples/Cooling/FeedbackEvent_3D/plotSolution.py b/examples/Cooling/FeedbackEvent_3D/plotSolution.py
index fe6f93996dafee2b6e81a70d4786374d86355f6f..6631fff2244e73244d0311f4e823c42dfcea7609 100644
--- a/examples/Cooling/FeedbackEvent_3D/plotSolution.py
+++ b/examples/Cooling/FeedbackEvent_3D/plotSolution.py
@@ -50,7 +50,6 @@ except:
     plt.rcParams.update(rcParams)
 
 
-
 def get_data_dump(metadata):
     """
     Gets a big data dump from the SWIFT metadata
@@ -104,10 +103,10 @@ r = r.to(kpc).value
 data = dict(
     x=r,
     v_r=v_r,
-    u=sim.gas.temperature.to(K).value,
-    S=sim.gas.entropy.to(keV / K).value,
-    P=sim.gas.pressure.to(kPa).value,
-    rho=sim.gas.density.to(mh / (cm ** 3)).value,
+    u=sim.gas.temperatures.to(K).value,
+    S=sim.gas.entropies.to(keV / K).value,
+    P=sim.gas.pressures.to(kPa).value,
+    rho=sim.gas.densities.to(mh / (cm ** 3)).value,
 )
 
 # Try to add on the viscosity and diffusion.
@@ -156,8 +155,13 @@ log = dict(
     v_r=False, v_phi=False, u=False, S=False, P=False, rho=False, visc=False, diff=False
 )
 ylim = dict(
-    v_r=[-4, 25], u=[4750, 6000], rho=[0.09, 0.16], visc=[0, 2.0], diff=[0, 0.25],
-    P=[3e-18, 10e-18], S=[-0.5e60, 4e60] 
+    v_r=[-4, 25],
+    u=[4750, 6000],
+    rho=[0.09, 0.16],
+    visc=[0, 2.0],
+    diff=[0, 0.25],
+    P=[3e-18, 10e-18],
+    S=[-0.5e60, 4e60],
 )
 
 current_axis = 0
diff --git a/examples/Cosmology/ComovingSodShock_1D/plotSolution.py b/examples/Cosmology/ComovingSodShock_1D/plotSolution.py
index 95674c04bfafd0cd549b69814df82f9a4f80a949..b09873339fc4ded024c70f66301c5cce762b97fc 100644
--- a/examples/Cosmology/ComovingSodShock_1D/plotSolution.py
+++ b/examples/Cosmology/ComovingSodShock_1D/plotSolution.py
@@ -82,12 +82,12 @@ git = str(sim["Code"].attrs["Git Revision"])
 
 x = sim["/PartType0/Coordinates"][:,0]
 v = sim["/PartType0/Velocities"][:,0] * anow
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
 try:
-    alpha = sim["/PartType0/Viscosity"][:]
+    alpha = sim["/PartType0/ViscosityParameters"][:]
     plot_alpha = True 
 except:
     plot_alpha = False
diff --git a/examples/Cosmology/ComovingSodShock_2D/plotSolution.py b/examples/Cosmology/ComovingSodShock_2D/plotSolution.py
index 8adb3cf5c550ab9724f6a8f34c1a1260a25712e1..ec28b9449bfa6b00d3088952a6b1bf871462f86c 100644
--- a/examples/Cosmology/ComovingSodShock_2D/plotSolution.py
+++ b/examples/Cosmology/ComovingSodShock_2D/plotSolution.py
@@ -83,10 +83,10 @@ git = sim["Code"].attrs["Git Revision"]
 
 x = sim["/PartType0/Coordinates"][:,0]
 v = sim["/PartType0/Velocities"][:,0] * anow
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
 
 N = 1000  # Number of points
 x_min = -1.
diff --git a/examples/Cosmology/ComovingSodShock_3D/plotSolution.py b/examples/Cosmology/ComovingSodShock_3D/plotSolution.py
index d85f4be9a49d108d133928a81ea4482fa9099792..34abae364f5421a0436352b9564537f2f31e8742 100644
--- a/examples/Cosmology/ComovingSodShock_3D/plotSolution.py
+++ b/examples/Cosmology/ComovingSodShock_3D/plotSolution.py
@@ -83,19 +83,19 @@ git = sim["Code"].attrs["Git Revision"]
 
 x = sim["/PartType0/Coordinates"][:,0]
 v = sim["/PartType0/Velocities"][:,0] * anow
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
 
 try:
-    diffusion = sim["/PartType0/Diffusion"][:]
+    diffusion = sim["/PartType0/DiffusionParameters"][:]
     plot_diffusion = True
 except:
     plot_diffusion = False
 
 try:
-    viscosity = sim["/PartType0/Viscosity"][:]
+    viscosity = sim["/PartType0/ViscosityParameters"][:]
     plot_viscosity = True
 except:
     plot_viscosity = False
diff --git a/examples/Cosmology/ConstantCosmoVolume/plotSolution.py b/examples/Cosmology/ConstantCosmoVolume/plotSolution.py
index f77889d7cb19c230accf25290b88a321e0f59616..6df0bfa4e1f7c4932f652d82a8ca95f5c54b0e9f 100644
--- a/examples/Cosmology/ConstantCosmoVolume/plotSolution.py
+++ b/examples/Cosmology/ConstantCosmoVolume/plotSolution.py
@@ -109,19 +109,19 @@ for i in range(119):
     z[i] = sim["/Cosmology"].attrs["Redshift"][0]
     a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
     
-    S = sim["/PartType0/Entropy"][:]
+    S = sim["/PartType0/Entropies"][:]
     S_mean[i] = np.mean(S)
     S_std[i] = np.std(S)
     
-    u = sim["/PartType0/InternalEnergy"][:]
+    u = sim["/PartType0/InternalEnergies"][:]
     u_mean[i] = np.mean(u)
     u_std[i] = np.std(u)
 
-    P = sim["/PartType0/Pressure"][:]
+    P = sim["/PartType0/Pressures"][:]
     P_mean[i] = np.mean(P)
     P_std[i] = np.std(P)
 
-    rho = sim["/PartType0/Density"][:]
+    rho = sim["/PartType0/Densities"][:]
     rho_mean[i] = np.mean(rho)
     rho_std[i] = np.std(rho)
 
diff --git a/examples/Cosmology/ZeldovichPancake_3D/plotSolution.py b/examples/Cosmology/ZeldovichPancake_3D/plotSolution.py
index eef247fb761e75f8dde8e8abe84075efbd7cb46a..285ed3ae218549468c309e9109ef53f10a2f66ab 100644
--- a/examples/Cosmology/ZeldovichPancake_3D/plotSolution.py
+++ b/examples/Cosmology/ZeldovichPancake_3D/plotSolution.py
@@ -78,10 +78,10 @@ gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
 
 x = sim["/PartType0/Coordinates"][:,0]
 v = sim["/PartType0/Velocities"][:,0]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
 m = sim["/PartType0/Masses"][:]
 try:
     phi = sim["/PartType0/Potential"][:]
@@ -98,8 +98,8 @@ if os.path.exists(filename_g):
     sim_g = h5py.File(filename_g, "r")
     x_g = sim_g["/PartType0/Coordinates"][:,0]
     v_g = sim_g["/PartType0/Velocities"][:,0]
-    u_g = sim_g["/PartType0/InternalEnergy"][:]
-    rho_g = sim_g["/PartType0/Density"][:]
+    u_g = sim_g["/PartType0/InternalEnergies"][:]
+    rho_g = sim_g["/PartType0/Densities"][:]
     phi_g = sim_g["/PartType0/Potential"][:]
     a_g = sim_g["/Header"].attrs["Time"]
     print("Gadget Scale-factor:", a_g, "redshift:", 1/a_g - 1.)
diff --git a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
index 3270af86afffba13cfac6be7c7c21dd84980adc5..b5122afdba28ffaed1e52cdd8bd2db8a4e566c03 100644
--- a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
@@ -170,3 +170,4 @@ EAGLEAGN:
   coupling_efficiency:              0.15       # Fraction of the radiated energy that couples to the gas in feedback events.
   AGN_delta_T_K:                    3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
   AGN_num_ngb_to_heat:              1.         # Target number of gas neighbours to heat in an AGN feedback event.
+  max_reposition_mass:              2e8        # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
diff --git a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
index 6461982f70c9635652b8253887a9998777f37c69..888009f27761feda91d3cb70df2366b9f04fd3eb 100644
--- a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
@@ -171,3 +171,4 @@ EAGLEAGN:
   coupling_efficiency:              0.15       # Fraction of the radiated energy that couples to the gas in feedback events.
   AGN_delta_T_K:                    3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
   AGN_num_ngb_to_heat:              1.         # Target number of gas neighbours to heat in an AGN feedback event.
+  max_reposition_mass:              2e8        # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
diff --git a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
index c6f7f0e01bd292b03b0eb31c6c1b57d5163dafe7..f2e3b0656541ca698c20a85a308e538871d4fdc4 100644
--- a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
@@ -171,5 +171,4 @@ EAGLEAGN:
   coupling_efficiency:              0.15       # Fraction of the radiated energy that couples to the gas in feedback events.
   AGN_delta_T_K:                    3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
   AGN_num_ngb_to_heat:              1.         # Target number of gas neighbours to heat in an AGN feedback event.
-
-  
+  max_reposition_mass:              2e8        # Maximal BH mass considered for BH repositioning in solar masses.
diff --git a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
index 9205b9f6f560398b046625f6c65fcc887492991d..685d24d20f14f567175e1d7086a65455017199e0 100644
--- a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
@@ -162,3 +162,4 @@ EAGLEAGN:
   coupling_efficiency:              0.15       # Fraction of the radiated energy that couples to the gas in feedback events.
   AGN_delta_T_K:                    3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
   AGN_num_ngb_to_heat:              1.         # Target number of gas neighbours to heat in an AGN feedback event.
+  max_reposition_mass:              2e8        # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
diff --git a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
index ccd89efbead0c511316bd664f6754ed0797c7c50..b782631d7e264e582455f70b4fb399f7db1f4731 100644
--- a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
@@ -169,3 +169,4 @@ EAGLEAGN:
   coupling_efficiency:              0.15       # Fraction of the radiated energy that couples to the gas in feedback events.
   AGN_delta_T_K:                    3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
   AGN_num_ngb_to_heat:              1.         # Target number of gas neighbours to heat in an AGN feedback event.
+  max_reposition_mass:              2e8        # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
diff --git a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
index e2a98725a94c9fd34b45a95530f480562c0aaec8..53afda9445b4f8ee7c258b5c42af2a605eba4755 100644
--- a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
@@ -161,3 +161,4 @@ EAGLEAGN:
   coupling_efficiency:              0.15       # Fraction of the radiated energy that couples to the gas in feedback events.
   AGN_delta_T_K:                    3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
   AGN_num_ngb_to_heat:              1.         # Target number of gas neighbours to heat in an AGN feedback event.
+  max_reposition_mass:              2e8        # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
diff --git a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
index 9136f85424beaf70b614cde0e70b67c134056831..694c1c111d8f5b2455d4cfcdc567b8854b12ff2c 100644
--- a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
@@ -171,3 +171,4 @@ EAGLEAGN:
   coupling_efficiency:              0.15       # Fraction of the radiated energy that couples to the gas in feedback events.
   AGN_delta_T_K:                    3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
   AGN_num_ngb_to_heat:              1.         # Target number of gas neighbours to heat in an AGN feedback event.
+  max_reposition_mass:              2e8        # Maximal BH mass considered for BH repositioning in solar masses.
\ No newline at end of file
diff --git a/examples/HydroTests/BlobTest_3D/makeMovie.py b/examples/HydroTests/BlobTest_3D/makeMovie.py
index 9ae4a538e000fa7006f093b67d325ff433e97089..55acdaab01e22ecf8bae0dd24967a348ceabf34b 100644
--- a/examples/HydroTests/BlobTest_3D/makeMovie.py
+++ b/examples/HydroTests/BlobTest_3D/makeMovie.py
@@ -24,7 +24,7 @@ resolution = 1024
 snapshot_name = "blob"
 cmap = "Spectral_r"
 text_args = dict(color="black")
-# plot = "pressure"
+# plot = "pressures"
 # name = "Pressure $P$"
 plot = "density"
 name = "Fluid Density $\\rho$"
diff --git a/examples/HydroTests/Diffusion_1D/plotSolution.py b/examples/HydroTests/Diffusion_1D/plotSolution.py
index 66c8ffc6418f06589a2918ae4d8ed460b0081972..a394e919b4e2e80728c97499c3c3544256a6ad2c 100644
--- a/examples/HydroTests/Diffusion_1D/plotSolution.py
+++ b/examples/HydroTests/Diffusion_1D/plotSolution.py
@@ -24,6 +24,7 @@ import numpy as np
 
 try:
     from scipy.integrate import solve_ivp
+
     solve_ode = True
 except:
     solve_ode = False
@@ -32,6 +33,7 @@ from swiftsimio import load
 
 matplotlib.use("Agg")
 
+
 def solve_analytic(u_0, u_1, t_0, t_1, alpha=0.1):
     """
     Solves the analytic equation:
@@ -63,7 +65,12 @@ def solve_analytic(u_0, u_1, t_0, t_1, alpha=0.1):
 
         return np.array([-1.0 * common, 1.0 * common])
 
-    ret = solve_ivp(gradient, t_span=[t_0.value, t_1.value], y0=[u_0.value, u_1.value], t_eval=np.linspace(t_0.value, t_1.value, 100))
+    ret = solve_ivp(
+        gradient,
+        t_span=[t_0.value, t_1.value],
+        y0=[u_0.value, u_1.value],
+        t_eval=np.linspace(t_0.value, t_1.value, 100),
+    )
 
     t = ret.t
     high = ret.y[1]
@@ -81,7 +88,7 @@ def get_data_dump(metadata):
         viscosity = metadata.viscosity_info
     except:
         viscosity = "No info"
-    
+
     try:
         diffusion = metadata.diffusion_info
     except:
@@ -136,6 +143,7 @@ def setup_axes(size=[8, 8], dpi=300):
 
     return fig, ax
 
+
 def mean_std_max_min(data):
     """
     Returns:
@@ -157,7 +165,7 @@ def extract_plottables_u(data_list):
     """
 
     data = [
-        np.diff(x.gas.internal_energy.value) / np.mean(x.gas.internal_energy.value)
+        np.diff(x.gas.internal_energies.value) / np.mean(x.gas.internal_energies.value)
         for x in data_list
     ]
 
@@ -175,8 +183,10 @@ def extract_plottables_x(data_list):
     dx = boxsize / n_part
 
     original_x = np.arange(n_part, dtype=float) * dx + (0.5 * dx)
-    
-    deviations = [1e6 * abs(original_x - x.gas.coordinates.value[:, 0]) / dx for x in data_list]
+
+    deviations = [
+        1e6 * abs(original_x - x.gas.coordinates.value[:, 0]) / dx for x in data_list
+    ]
 
     return mean_std_max_min(deviations)
 
@@ -187,7 +197,7 @@ def extract_plottables_rho(data_list):
     mean, stdev, max, min * 1e6 deviations from mean density 
     """
 
-    P = [x.gas.density.value for x in data_list]
+    P = [x.gas.densities.value for x in data_list]
     mean_P = [np.mean(x) for x in P]
     deviations = [1e6 * (x - y) / x for x, y in zip(mean_P, P)]
 
@@ -241,28 +251,34 @@ def make_plot(start: int, stop: int, handle: str):
 
     if solve_ode:
         times_ode, diff = solve_analytic(
-            u_0=data_list[0].gas.internal_energy.min(),
-            u_1=data_list[0].gas.internal_energy.max(),
+            u_0=data_list[0].gas.internal_energies.min(),
+            u_1=data_list[0].gas.internal_energies.max(),
             t_0=t[0],
             t_1=t[-1],
             alpha=(
-                np.sqrt(5.0/3.0 * (5.0/3.0 - 1.0)) * 
-                alpha / data_list[0].gas.smoothing_length[0].value
-            )
+                np.sqrt(5.0 / 3.0 * (5.0 / 3.0 - 1.0))
+                * alpha
+                / data_list[0].gas.smoothing_length[0].value
+            ),
         )
 
         ax[1].plot(
             times_ode,
-            (diff) / np.mean(data_list[0].gas.internal_energy),
+            (diff) / np.mean(data_list[0].gas.internal_energies),
             label="Analytic",
             linestyle="dotted",
-            c="C3"
+            c="C3",
         )
 
-        #import pdb;pdb.set_trace()
+        # import pdb;pdb.set_trace()
 
     ax[2].fill_between(
-        t, x_means - x_stdevs, x_means + x_stdevs, color="C0", alpha=0.5, edgecolor="none"
+        t,
+        x_means - x_stdevs,
+        x_means + x_stdevs,
+        color="C0",
+        alpha=0.5,
+        edgecolor="none",
     )
     ax[2].plot(t, x_means, label="Mean", c="C0")
     ax[2].plot(t, x_maxs, label="Max", linestyle="dashed", c="C1")
@@ -270,11 +286,18 @@ def make_plot(start: int, stop: int, handle: str):
 
     try:
         # Give diffusion info a go; this may not be present
-        diff_means, diff_stdevs, diff_maxs, diff_mins = extract_plottables_diff(data_list)
+        diff_means, diff_stdevs, diff_maxs, diff_mins = extract_plottables_diff(
+            data_list
+        )
 
         ax[3].set_ylabel(r"Diffusion parameter $\alpha_{diff}$")
         ax[3].fill_between(
-            t, diff_means - diff_stdevs, diff_means + diff_stdevs, color="C0", alpha=0.5, edgecolor="none"
+            t,
+            diff_means - diff_stdevs,
+            diff_means + diff_stdevs,
+            color="C0",
+            alpha=0.5,
+            edgecolor="none",
         )
         ax[3].plot(t, diff_means, label="Mean", c="C0")
         ax[3].plot(t, diff_maxs, label="Max", linestyle="dashed", c="C1")
@@ -284,9 +307,16 @@ def make_plot(start: int, stop: int, handle: str):
         # Diffusion info must not be present.
         rho_means, rho_stdevs, rho_maxs, rho_mins = extract_plottables_rho(data_list)
 
-        ax[3].set_ylabel("Deviation from mean density $(\\rho_i - \\bar{\\rho}) / \\bar{\\rho}$ [$\\times 10^{6}$]")
+        ax[3].set_ylabel(
+            "Deviation from mean density $(\\rho_i - \\bar{\\rho}) / \\bar{\\rho}$ [$\\times 10^{6}$]"
+        )
         ax[3].fill_between(
-            t, rho_means - rho_stdevs, rho_means + rho_stdevs, color="C0", alpha=0.5, edgecolor="none"
+            t,
+            rho_means - rho_stdevs,
+            rho_means + rho_stdevs,
+            color="C0",
+            alpha=0.5,
+            edgecolor="none",
         )
         ax[3].plot(t, rho_means, label="Mean", c="C0")
         ax[3].plot(t, rho_maxs, label="Max", linestyle="dashed", c="C1")
diff --git a/examples/HydroTests/EvrardCollapse_3D/plotSolution.py b/examples/HydroTests/EvrardCollapse_3D/plotSolution.py
index 8422b9c45fd573f3d0ae36324d6e39ab23cceb25..b405771a4792e5ca4fef181b3787952bf0078d67 100644
--- a/examples/HydroTests/EvrardCollapse_3D/plotSolution.py
+++ b/examples/HydroTests/EvrardCollapse_3D/plotSolution.py
@@ -75,10 +75,10 @@ x = sqrt((coords[:,0] - 0.5 * boxSize)**2 + (coords[:,1] - 0.5 * boxSize)**2 + \
          (coords[:,2] - 0.5 * boxSize)**2)
 vels = sim["/PartType0/Velocities"]
 v = sqrt(vels[:,0]**2 + vels[:,1]**2 + vels[:,2]**2)
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
 
 # Bin the data
 x_bin_edge = logspace(-3., log10(2.), 100)
diff --git a/examples/HydroTests/Gradients/plot.py b/examples/HydroTests/Gradients/plot.py
index d6750ffc581437ebbf402ec44bcb1d14cb82a698..7b4248c9fc9c9d6b9c5b15410185d6489849bed8 100644
--- a/examples/HydroTests/Gradients/plot.py
+++ b/examples/HydroTests/Gradients/plot.py
@@ -30,7 +30,7 @@ inputfile = sys.argv[1]
 outputfile = "gradients_{0}.png".format(sys.argv[2])
 
 f = h5py.File(inputfile, "r")
-rho = np.array(f["/PartType0/Density"])
+rho = np.array(f["/PartType0/Densities"])
 gradrho = np.array(f["/PartType0/GradDensity"])
 coords = np.array(f["/PartType0/Coordinates"])
 
diff --git a/examples/HydroTests/GreshoVortex_2D/plotSolution.py b/examples/HydroTests/GreshoVortex_2D/plotSolution.py
index d497a6b297bf38b39cf85a9107a769c20f815b77..2d4697b6ffaac0639da67ee90d824c75791ea573 100644
--- a/examples/HydroTests/GreshoVortex_2D/plotSolution.py
+++ b/examples/HydroTests/GreshoVortex_2D/plotSolution.py
@@ -100,10 +100,10 @@ r = sqrt(x**2 + y**2)
 v_r = (x * vel[:,0] + y * vel[:,1]) / r
 v_phi = (-y * vel[:,0] + x * vel[:,1]) / r
 v_norm = sqrt(vel[:,0]**2 + vel[:,1]**2)
-rho = sim["/PartType0/Density"][:]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Densities"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
 
 # Bin te data
 r_bin_edge = np.arange(0., 1., 0.02)
diff --git a/examples/HydroTests/GreshoVortex_3D/plotSolution.py b/examples/HydroTests/GreshoVortex_3D/plotSolution.py
index 545440c997d9ebc3ab11562d0a7d9fa143e23ed2..20beab7514759c764f5ca7c379183506b764a819 100644
--- a/examples/HydroTests/GreshoVortex_3D/plotSolution.py
+++ b/examples/HydroTests/GreshoVortex_3D/plotSolution.py
@@ -103,19 +103,19 @@ r = sqrt(x**2 + y**2)
 v_r = (x * vel[:,0] + y * vel[:,1]) / r
 v_phi = (-y * vel[:,0] + x * vel[:,1]) / r
 v_norm = sqrt(vel[:,0]**2 + vel[:,1]**2)
-rho = sim["/PartType0/Density"][:]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Densities"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
 
 try:
-    diffusion = sim["/PartType0/Diffusion"][:]
+    diffusion = sim["/PartType0/DiffusionParameters"][:]
     plot_diffusion = True
 except:
     plot_diffusion = False
 
 try:
-    viscosity = sim["/PartType0/Viscosity"][:]
+    viscosity = sim["/PartType0/ViscosityParameters"][:]
     plot_viscosity = True
 except:
     plot_viscosity = False
diff --git a/examples/HydroTests/InteractingBlastWaves_1D/plotSolution.py b/examples/HydroTests/InteractingBlastWaves_1D/plotSolution.py
index 1719162dec34626d6f4ecb8267c4d06f85b3db26..d617fb239ce21acab73b5cb057dd3cdf4b260d59 100644
--- a/examples/HydroTests/InteractingBlastWaves_1D/plotSolution.py
+++ b/examples/HydroTests/InteractingBlastWaves_1D/plotSolution.py
@@ -55,11 +55,11 @@ snap = int(sys.argv[1])
 # Open the file and read the relevant data
 file = h5py.File("interactingBlastWaves_{0:04d}.hdf5".format(snap), "r")
 x = file["/PartType0/Coordinates"][:,0]
-rho = file["/PartType0/Density"]
+rho = file["/PartType0/Densities"]
 v = file["/PartType0/Velocities"][:,0]
-u = file["/PartType0/InternalEnergy"]
-S = file["/PartType0/Entropy"]
-P = file["/PartType0/Pressure"]
+u = file["/PartType0/InternalEnergies"]
+S = file["/PartType0/Entropies"]
+P = file["/PartType0/Pressures"]
 time = file["/Header"].attrs["Time"][0]
 
 scheme = file["/HydroScheme"].attrs["Scheme"]
diff --git a/examples/HydroTests/KelvinHelmholtz_2D/plotSolution.py b/examples/HydroTests/KelvinHelmholtz_2D/plotSolution.py
index 77ab6fb244da25d13760f90653fac7eac11a0ee7..f599fcb784633b2d6765ea79767fc658196faa5f 100644
--- a/examples/HydroTests/KelvinHelmholtz_2D/plotSolution.py
+++ b/examples/HydroTests/KelvinHelmholtz_2D/plotSolution.py
@@ -77,10 +77,10 @@ x = pos[:,0] - boxSize / 2
 y = pos[:,1] - boxSize / 2
 vel = sim["/PartType0/Velocities"][:,:]
 v_norm = sqrt(vel[:,0]**2 + vel[:,1]**2)
-rho = sim["/PartType0/Density"][:]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Densities"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
 
 # Plot the interesting quantities
 figure()
diff --git a/examples/HydroTests/Noh_1D/plotSolution.py b/examples/HydroTests/Noh_1D/plotSolution.py
index 25b9b2f16b24cba5def592a5cf00dbae82195ef7..7f0b5d403ef816b0dda57823010472476a7ecc32 100644
--- a/examples/HydroTests/Noh_1D/plotSolution.py
+++ b/examples/HydroTests/Noh_1D/plotSolution.py
@@ -69,10 +69,10 @@ git = sim["Code"].attrs["Git Revision"]
 
 x = sim["/PartType0/Coordinates"][:,0]
 v = sim["/PartType0/Velocities"][:,0]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
 
 N = 1001  # Number of points
 x_min = -1
diff --git a/examples/HydroTests/Noh_2D/plotSolution.py b/examples/HydroTests/Noh_2D/plotSolution.py
index 775ddf4e8a7954c14034ad51a6b66622c41a6996..b53212c4688ec790ad8f3f83f81243f9ec52266d 100644
--- a/examples/HydroTests/Noh_2D/plotSolution.py
+++ b/examples/HydroTests/Noh_2D/plotSolution.py
@@ -71,10 +71,10 @@ x = sim["/PartType0/Coordinates"][:,0]
 y = sim["/PartType0/Coordinates"][:,1]
 vx = sim["/PartType0/Velocities"][:,0]
 vy = sim["/PartType0/Velocities"][:,1]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
 
 r = np.sqrt((x-1)**2 + (y-1)**2)
 v = -np.sqrt(vx**2 + vy**2)
diff --git a/examples/HydroTests/Noh_3D/plotSolution.py b/examples/HydroTests/Noh_3D/plotSolution.py
index 386b9f728b5e8d8e38fb7ec9aeaa336d194e35dd..20e8ca805de1cb700b8b462ae27495080f5d3268 100644
--- a/examples/HydroTests/Noh_3D/plotSolution.py
+++ b/examples/HydroTests/Noh_3D/plotSolution.py
@@ -74,10 +74,10 @@ z = sim["/PartType0/Coordinates"][:,2]
 vx = sim["/PartType0/Velocities"][:,0]
 vy = sim["/PartType0/Velocities"][:,1]
 vz = sim["/PartType0/Velocities"][:,2]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
 
 r = np.sqrt((x-1)**2 + (y-1)**2 + (z-1)**2)
 v = -np.sqrt(vx**2 + vy**2 + vz**2)
diff --git a/examples/HydroTests/Rayleigh-Taylor_2D/plotInitialProfile.py b/examples/HydroTests/Rayleigh-Taylor_2D/plotInitialProfile.py
index e89c4e8525fe4c88e517acbd453b0941f8f573c8..8928a6e597adbe4e52f905ba95fa68f424b6cabb 100644
--- a/examples/HydroTests/Rayleigh-Taylor_2D/plotInitialProfile.py
+++ b/examples/HydroTests/Rayleigh-Taylor_2D/plotInitialProfile.py
@@ -12,8 +12,8 @@ f = load(filename)
 
 # Get data from snapshot
 x, y, _ = f.gas.coordinates.value.T
-rho = f.gas.density.value
-a = f.gas.entropy.value
+rho = f.gas.densities.value
+a = f.gas.entropies.value
 
 # Get analytical solution
 y_an = np.linspace(0, makeIC.box_size[1], N)
diff --git a/examples/HydroTests/SedovBlast_1D/plotSolution.py b/examples/HydroTests/SedovBlast_1D/plotSolution.py
index c6d4a989da7493f7b500946610eea8832696bf6f..d82ad9e94610a916d900ea93863b2881757e73b3 100644
--- a/examples/HydroTests/SedovBlast_1D/plotSolution.py
+++ b/examples/HydroTests/SedovBlast_1D/plotSolution.py
@@ -78,13 +78,13 @@ x = pos[:,0] - boxSize / 2
 vel = sim["/PartType0/Velocities"][:,:]
 r = abs(x)
 v_r = x * vel[:,0] / r
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
 
 try:
-    alpha = sim["/PartType0/Viscosity"][:]
+    alpha = sim["/PartType0/ViscosityParameters"][:]
     plot_alpha = True 
 except:
     plot_alpha = False
diff --git a/examples/HydroTests/SedovBlast_2D/plotSolution.py b/examples/HydroTests/SedovBlast_2D/plotSolution.py
index 2b5de6f32b8673bbc825fbb5236f4e2ab3b4f408..6f504a09c9432368ce141ec0d28c28699f5ba7f3 100644
--- a/examples/HydroTests/SedovBlast_2D/plotSolution.py
+++ b/examples/HydroTests/SedovBlast_2D/plotSolution.py
@@ -80,10 +80,10 @@ y = pos[:,1] - boxSize / 2
 vel = sim["/PartType0/Velocities"][:,:]
 r = sqrt(x**2 + y**2)
 v_r = (x * vel[:,0] + y * vel[:,1]) / r
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
 
 # Bin te data
 r_bin_edge = np.arange(0., 0.5, 0.01)
diff --git a/examples/HydroTests/SedovBlast_3D/plotSolution.py b/examples/HydroTests/SedovBlast_3D/plotSolution.py
index b0f2e08441b3fa550e61602ba852228a04362fbc..fec4f1101406be5803a3f1601812d1cb85275409 100644
--- a/examples/HydroTests/SedovBlast_3D/plotSolution.py
+++ b/examples/HydroTests/SedovBlast_3D/plotSolution.py
@@ -81,19 +81,19 @@ z = pos[:,2] - boxSize / 2
 vel = sim["/PartType0/Velocities"][:,:]
 r = sqrt(x**2 + y**2 + z**2)
 v_r = (x * vel[:,0] + y * vel[:,1] + z * vel[:,2]) / r
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
 
 try:
-    diffusion = sim["/PartType0/Diffusion"][:]
+    diffusion = sim["/PartType0/DiffusionParameters"][:]
     plot_diffusion = True
 except:
     plot_diffusion = False
 
 try:
-    viscosity = sim["/PartType0/Viscosity"][:]
+    viscosity = sim["/PartType0/ViscosityParameters"][:]
     plot_viscosity = True
 except:
     plot_viscosity = False
diff --git a/examples/HydroTests/SineWavePotential_1D/plotSolution.py b/examples/HydroTests/SineWavePotential_1D/plotSolution.py
index 3bb889aaabd3cdac0274afb09647d0e3aebb02cc..ae99d98aaaa11d9c68473b74106054963a075895 100644
--- a/examples/HydroTests/SineWavePotential_1D/plotSolution.py
+++ b/examples/HydroTests/SineWavePotential_1D/plotSolution.py
@@ -43,9 +43,9 @@ fileName = sys.argv[1]
 
 file = h5py.File(fileName, 'r')
 coords = np.array(file["/PartType0/Coordinates"])
-rho = np.array(file["/PartType0/Density"])
-P = np.array(file["/PartType0/Pressure"])
-u = np.array(file["/PartType0/InternalEnergy"])
+rho = np.array(file["/PartType0/Densities"])
+P = np.array(file["/PartType0/Pressures"])
+u = np.array(file["/PartType0/InternalEnergies"])
 m = np.array(file["/PartType0/Masses"])
 vs = np.array(file["/PartType0/Velocities"])
 ids = np.array(file["/PartType0/ParticleIDs"])
diff --git a/examples/HydroTests/SineWavePotential_2D/plotSolution.py b/examples/HydroTests/SineWavePotential_2D/plotSolution.py
index ee02f59c404db87a790465d2786e6296525e36b0..5c87b0f4f3682d486063522715517763b1035567 100644
--- a/examples/HydroTests/SineWavePotential_2D/plotSolution.py
+++ b/examples/HydroTests/SineWavePotential_2D/plotSolution.py
@@ -38,8 +38,8 @@ fileName = sys.argv[1]
 
 file = h5py.File(fileName, 'r')
 coords = np.array(file["/PartType0/Coordinates"])
-rho = np.array(file["/PartType0/Density"])
-u = np.array(file["/PartType0/InternalEnergy"])
+rho = np.array(file["/PartType0/Densities"])
+u = np.array(file["/PartType0/InternalEnergies"])
 agrav = np.array(file["/PartType0/GravAcceleration"])
 m = np.array(file["/PartType0/Masses"])
 ids = np.array(file["/PartType0/ParticleIDs"])
diff --git a/examples/HydroTests/SineWavePotential_3D/plotSolution.py b/examples/HydroTests/SineWavePotential_3D/plotSolution.py
index 13cae037b64eff4ad4fec0003bf0f5ad3ce94896..7bfa82a5990572c478976614c50107e5254f0e00 100644
--- a/examples/HydroTests/SineWavePotential_3D/plotSolution.py
+++ b/examples/HydroTests/SineWavePotential_3D/plotSolution.py
@@ -38,8 +38,8 @@ fileName = sys.argv[1]
 
 file = h5py.File(fileName, 'r')
 coords = np.array(file["/PartType0/Coordinates"])
-rho = np.array(file["/PartType0/Density"])
-u = np.array(file["/PartType0/InternalEnergy"])
+rho = np.array(file["/PartType0/Densities"])
+u = np.array(file["/PartType0/InternalEnergies"])
 agrav = np.array(file["/PartType0/GravAcceleration"])
 m = np.array(file["/PartType0/Masses"])
 ids = np.array(file["/PartType0/ParticleIDs"])
diff --git a/examples/HydroTests/SodShockSpherical_2D/plotSolution.py b/examples/HydroTests/SodShockSpherical_2D/plotSolution.py
index 57b7f7ddc64bc25df031eb0cba7547f40d46b31a..61060631eea1a7320ef207457d35031318bceccf 100644
--- a/examples/HydroTests/SodShockSpherical_2D/plotSolution.py
+++ b/examples/HydroTests/SodShockSpherical_2D/plotSolution.py
@@ -74,10 +74,10 @@ coords = sim["/PartType0/Coordinates"]
 x = sqrt((coords[:,0] - 0.5)**2 + (coords[:,1] - 0.5)**2)
 vels = sim["/PartType0/Velocities"]
 v = sqrt(vels[:,0]**2 + vels[:,1]**2)
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
 
 # Bin the data
 rho_bin,x_bin_edge,_ = \
diff --git a/examples/HydroTests/SodShockSpherical_3D/plotSolution.py b/examples/HydroTests/SodShockSpherical_3D/plotSolution.py
index 539bfba799e3231bd26ae2eb39c271baa1fa6a4b..0a92f3aaf1831d67cb59dd71bc08cd8d973d9def 100644
--- a/examples/HydroTests/SodShockSpherical_3D/plotSolution.py
+++ b/examples/HydroTests/SodShockSpherical_3D/plotSolution.py
@@ -75,10 +75,10 @@ x = sqrt((coords[:,0] - 0.5)**2 + (coords[:,1] - 0.5)**2 + \
          (coords[:,2] - 0.5)**2)
 vels = sim["/PartType0/Velocities"]
 v = sqrt(vels[:,0]**2 + vels[:,1]**2 + vels[:,2]**2)
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
 
 # Bin the data
 rho_bin,x_bin_edge,_ = \
diff --git a/examples/HydroTests/SodShock_1D/plotSolution.py b/examples/HydroTests/SodShock_1D/plotSolution.py
index a7e6d374bac616440dace666b85c3e7ade479bcd..770d05ac0493323c2cdd0f5905e409113d1a9eae 100644
--- a/examples/HydroTests/SodShock_1D/plotSolution.py
+++ b/examples/HydroTests/SodShock_1D/plotSolution.py
@@ -79,12 +79,12 @@ git = str(sim["Code"].attrs["Git Revision"])
 
 x = sim["/PartType0/Coordinates"][:,0]
 v = sim["/PartType0/Velocities"][:,0]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
 try:
-    alpha = sim["/PartType0/Viscosity"][:]
+    alpha = sim["/PartType0/ViscosityParameters"][:]
     plot_alpha = True 
 except:
     plot_alpha = False
diff --git a/examples/HydroTests/SodShock_2D/plotSolution.py b/examples/HydroTests/SodShock_2D/plotSolution.py
index 19cbe0ffb766845c051ffb6cea81bd918d890e36..769079da8824d58535e239bd8b54b592ce981a37 100644
--- a/examples/HydroTests/SodShock_2D/plotSolution.py
+++ b/examples/HydroTests/SodShock_2D/plotSolution.py
@@ -79,10 +79,10 @@ git = sim["Code"].attrs["Git Revision"]
 
 x = sim["/PartType0/Coordinates"][:,0]
 v = sim["/PartType0/Velocities"][:,0]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
 
 N = 1000  # Number of points
 x_min = -1.
diff --git a/examples/HydroTests/SodShock_3D/plotSolution.py b/examples/HydroTests/SodShock_3D/plotSolution.py
index 69b2fe4887e986156ed01e0f4177d01ccbed6035..c2028cedcea820e56c07962fe3ca2f1fe1347d40 100644
--- a/examples/HydroTests/SodShock_3D/plotSolution.py
+++ b/examples/HydroTests/SodShock_3D/plotSolution.py
@@ -79,19 +79,19 @@ git = sim["Code"].attrs["Git Revision"]
 
 x = sim["/PartType0/Coordinates"][:,0]
 v = sim["/PartType0/Velocities"][:,0]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
 
 try:
-    diffusion = sim["/PartType0/Diffusion"][:]
+    diffusion = sim["/PartType0/DiffusionParameters"][:]
     plot_diffusion = True
 except:
     plot_diffusion = False
 
 try:
-    viscosity = sim["/PartType0/Viscosity"][:]
+    viscosity = sim["/PartType0/ViscosityParameters"][:]
     plot_viscosity = True
 except:
     plot_viscosity = False
diff --git a/examples/HydroTests/SodShock_BCC_3D/plotSolution.py b/examples/HydroTests/SodShock_BCC_3D/plotSolution.py
index 365b679991e9a3a5bbb9e9d5108066c04e497c2f..9660e12a3d226bd6b0e3c152031c93cedb345933 100644
--- a/examples/HydroTests/SodShock_BCC_3D/plotSolution.py
+++ b/examples/HydroTests/SodShock_BCC_3D/plotSolution.py
@@ -94,10 +94,10 @@ time = sim.metadata.t.value
 data = dict(
     x=sim.gas.coordinates.value[:, 0],
     v=sim.gas.velocities.value[:, 0],
-    u=sim.gas.internal_energy.value,
-    S=sim.gas.entropy.value,
-    P=sim.gas.pressure.value,
-    rho=sim.gas.density.value,
+    u=sim.gas.internal_energies.value,
+    S=sim.gas.entropies.value,
+    P=sim.gas.pressures.value,
+    rho=sim.gas.densities.value,
     y=sim.gas.coordinates.value[:, 1],
     z=sim.gas.coordinates.value[:, 2],
 )
@@ -164,12 +164,9 @@ for key, label in plot.items():
             zorder=-1,
         )
 
-        mask_noraster = np.logical_and.reduce([
-            data["y"] < 0.52,
-            data["y"] > 0.48,
-            data["z"] < 0.52,
-            data["z"] > 0.48
-        ])
+        mask_noraster = np.logical_and.reduce(
+            [data["y"] < 0.52, data["y"] > 0.48, data["z"] < 0.52, data["z"] > 0.48]
+        )
 
         axis.plot(
             data["x"][mask_noraster],
diff --git a/examples/HydroTests/SquareTest_2D/plotSolutionLegacy.py b/examples/HydroTests/SquareTest_2D/plotSolutionLegacy.py
index 956da800c9096232d2e82cf4cff4c780672e0a8f..d8701c3d44390f1d2637f798c0e9af23531c4600 100644
--- a/examples/HydroTests/SquareTest_2D/plotSolutionLegacy.py
+++ b/examples/HydroTests/SquareTest_2D/plotSolutionLegacy.py
@@ -88,10 +88,10 @@ while centre_y < 0.:
 pos = sim["/PartType0/Coordinates"][:,:]
 vel = sim["/PartType0/Velocities"][:,:]
 v_norm = sqrt(vel[:,0]**2 + vel[:,1]**2)
-rho = sim["/PartType0/Density"][:]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Densities"][:]
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
 x = pos[:,0] - centre_x
 y = pos[:,1] - centre_y
 
diff --git a/examples/HydroTests/VacuumSpherical_2D/plotSolution.py b/examples/HydroTests/VacuumSpherical_2D/plotSolution.py
index 6a65206ae20ccf79392054d047ba6be04f169f3e..de551c0ef4a5f4e027402c881069b7b4780f43d8 100644
--- a/examples/HydroTests/VacuumSpherical_2D/plotSolution.py
+++ b/examples/HydroTests/VacuumSpherical_2D/plotSolution.py
@@ -63,12 +63,12 @@ snap = int(sys.argv[1])
 file = h5py.File("vacuum_{0:04d}.hdf5".format(snap), "r")
 coords = file["/PartType0/Coordinates"]
 x = np.sqrt((coords[:,0] - 0.5)**2 + (coords[:,1] - 0.5)**2)
-rho = file["/PartType0/Density"][:]
+rho = file["/PartType0/Densities"][:]
 vels = file["/PartType0/Velocities"]
 v = np.sqrt(vels[:,0]**2 + vels[:,1]**2)
-u = file["/PartType0/InternalEnergy"][:]
-S = file["/PartType0/Entropy"][:]
-P = file["/PartType0/Pressure"][:]
+u = file["/PartType0/InternalEnergies"][:]
+S = file["/PartType0/Entropies"][:]
+P = file["/PartType0/Pressures"][:]
 time = file["/Header"].attrs["Time"][0]
 
 scheme = file["/HydroScheme"].attrs["Scheme"]
diff --git a/examples/HydroTests/VacuumSpherical_3D/plotSolution.py b/examples/HydroTests/VacuumSpherical_3D/plotSolution.py
index c73e48ee2d311692cdf4aa3b0e52f4766b339df8..4a04acde575eae46de67c70b536b8774befd96ae 100644
--- a/examples/HydroTests/VacuumSpherical_3D/plotSolution.py
+++ b/examples/HydroTests/VacuumSpherical_3D/plotSolution.py
@@ -64,12 +64,12 @@ file = h5py.File("vacuum_{0:04d}.hdf5".format(snap), "r")
 coords = file["/PartType0/Coordinates"]
 x = np.sqrt((coords[:,0] - 0.5)**2 + (coords[:,1] - 0.5)**2 + \
             (coords[:,2] - 0.5)**2)
-rho = file["/PartType0/Density"][:]
+rho = file["/PartType0/Densities"][:]
 vels = file["/PartType0/Velocities"]
 v = np.sqrt(vels[:,0]**2 + vels[:,1]**2 + vels[:,2]**2)
-u = file["/PartType0/InternalEnergy"][:]
-S = file["/PartType0/Entropy"][:]
-P = file["/PartType0/Pressure"][:]
+u = file["/PartType0/InternalEnergies"][:]
+S = file["/PartType0/Entropies"][:]
+P = file["/PartType0/Pressures"][:]
 time = file["/Header"].attrs["Time"][0]
 
 scheme = file["/HydroScheme"].attrs["Scheme"]
diff --git a/examples/HydroTests/Vacuum_1D/plotSolution.py b/examples/HydroTests/Vacuum_1D/plotSolution.py
index fceac10c25fd58b5bbcb6e31884cd62b4cfd61f5..eac7dc9e3ac43822ad167372f9f33bf2f5af0e2a 100644
--- a/examples/HydroTests/Vacuum_1D/plotSolution.py
+++ b/examples/HydroTests/Vacuum_1D/plotSolution.py
@@ -61,11 +61,11 @@ snap = int(sys.argv[1])
 # Open the file and read the relevant data
 file = h5py.File("vacuum_{0:04d}.hdf5".format(snap), "r")
 x = file["/PartType0/Coordinates"][:,0]
-rho = file["/PartType0/Density"]
+rho = file["/PartType0/Densities"]
 v = file["/PartType0/Velocities"][:,0]
-u = file["/PartType0/InternalEnergy"]
-S = file["/PartType0/Entropy"]
-P = file["/PartType0/Pressure"]
+u = file["/PartType0/InternalEnergies"]
+S = file["/PartType0/Entropies"]
+P = file["/PartType0/Pressures"]
 time = file["/Header"].attrs["Time"][0]
 
 scheme = file["/HydroScheme"].attrs["Scheme"]
diff --git a/examples/HydroTests/Vacuum_2D/plotSolution.py b/examples/HydroTests/Vacuum_2D/plotSolution.py
index 4d197234237df10b8cdbf197048a65991da023cf..ffd0eb1cdd857764f7ecc4e2d0c93fee3c5f29e8 100644
--- a/examples/HydroTests/Vacuum_2D/plotSolution.py
+++ b/examples/HydroTests/Vacuum_2D/plotSolution.py
@@ -62,11 +62,11 @@ snap = int(sys.argv[1])
 # Open the file and read the relevant data
 file = h5py.File("vacuum_{0:04d}.hdf5".format(snap), "r")
 x = file["/PartType0/Coordinates"][:,0]
-rho = file["/PartType0/Density"][:]
+rho = file["/PartType0/Densities"][:]
 v = file["/PartType0/Velocities"][:,0]
-u = file["/PartType0/InternalEnergy"][:]
-S = file["/PartType0/Entropy"][:]
-P = file["/PartType0/Pressure"][:]
+u = file["/PartType0/InternalEnergies"][:]
+S = file["/PartType0/Entropies"][:]
+P = file["/PartType0/Pressures"][:]
 time = file["/Header"].attrs["Time"][0]
 
 scheme = file["/HydroScheme"].attrs["Scheme"]
diff --git a/examples/HydroTests/Vacuum_3D/plotSolution.py b/examples/HydroTests/Vacuum_3D/plotSolution.py
index 4d197234237df10b8cdbf197048a65991da023cf..ffd0eb1cdd857764f7ecc4e2d0c93fee3c5f29e8 100644
--- a/examples/HydroTests/Vacuum_3D/plotSolution.py
+++ b/examples/HydroTests/Vacuum_3D/plotSolution.py
@@ -62,11 +62,11 @@ snap = int(sys.argv[1])
 # Open the file and read the relevant data
 file = h5py.File("vacuum_{0:04d}.hdf5".format(snap), "r")
 x = file["/PartType0/Coordinates"][:,0]
-rho = file["/PartType0/Density"][:]
+rho = file["/PartType0/Densities"][:]
 v = file["/PartType0/Velocities"][:,0]
-u = file["/PartType0/InternalEnergy"][:]
-S = file["/PartType0/Entropy"][:]
-P = file["/PartType0/Pressure"][:]
+u = file["/PartType0/InternalEnergies"][:]
+S = file["/PartType0/Entropies"][:]
+P = file["/PartType0/Pressures"][:]
 time = file["/Header"].attrs["Time"][0]
 
 scheme = file["/HydroScheme"].attrs["Scheme"]
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plotSolution.py b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plotSolution.py
index 89a87923148cb6872ab17b6d7229aef597ef3287..1ff8df3569f25590e5acb8046edeca0a1333d556 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plotSolution.py
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plotSolution.py
@@ -29,7 +29,7 @@ rcParams.update(params)
 rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
 
 snap = int(sys.argv[1])
-filename = "output_%.4d.hdf5"%snap
+filename = "output_%.4d.hdf5" % snap
 
 f = h5.File(filename, "r")
 
@@ -40,7 +40,7 @@ year_in_cgs = 3600.0 * 24 * 365.0
 Msun_in_cgs = 1.98848e33
 G_in_cgs = 6.67259e-8
 pc_in_cgs = 3.08567758e18
-Msun_p_pc2 = Msun_in_cgs / pc_in_cgs**2
+Msun_p_pc2 = Msun_in_cgs / pc_in_cgs ** 2
 
 # Gemoetry info
 boxsize = f["/Header"].attrs["BoxSize"]
@@ -52,66 +52,94 @@ unit_mass_in_cgs = f["/Units"].attrs["Unit mass in cgs (U_M)"]
 unit_time_in_cgs = f["/Units"].attrs["Unit time in cgs (U_t)"]
 
 # Calculate Gravitational constant in internal units
-G = G_in_cgs * ( unit_length_in_cgs**3 / unit_mass_in_cgs / unit_time_in_cgs**2)**(-1)
+G = G_in_cgs * (unit_length_in_cgs ** 3 / unit_mass_in_cgs / unit_time_in_cgs ** 2) ** (
+    -1
+)
 
 # Read parameters of the SF model
 KS_law_slope = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_exponent"])
 KS_law_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_normalisation"])
 KS_thresh_Z0 = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_Z0"])
 KS_thresh_slope = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_slope"])
-KS_thresh_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_norm_H_p_cm3"])
+KS_thresh_norm = float(
+    f["/Parameters"].attrs["EAGLEStarFormation:threshold_norm_H_p_cm3"]
+)
 KS_gas_fraction = float(f["/Parameters"].attrs["EAGLEStarFormation:gas_fraction"])
-KS_thresh_max_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_max_density_H_p_cm3"])
-KS_high_den_thresh = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_threshold_H_p_cm3"])
-KS_law_slope_high_den = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_exponent"])
-EOS_gamma_effective = float(f["/Parameters"].attrs["EAGLEStarFormation:EOS_gamma_effective"])                           
-EOS_density_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:EOS_density_norm_H_p_cm3"])                           
-EOS_temp_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:EOS_temperature_norm_K"])                           
+KS_thresh_max_norm = float(
+    f["/Parameters"].attrs["EAGLEStarFormation:threshold_max_density_H_p_cm3"]
+)
+KS_high_den_thresh = float(
+    f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_threshold_H_p_cm3"]
+)
+KS_law_slope_high_den = float(
+    f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_exponent"]
+)
+EOS_gamma_effective = float(
+    f["/Parameters"].attrs["EAGLEStarFormation:EOS_gamma_effective"]
+)
+EOS_density_norm = float(
+    f["/Parameters"].attrs["EAGLEStarFormation:EOS_density_norm_H_p_cm3"]
+)
+EOS_temp_norm = float(
+    f["/Parameters"].attrs["EAGLEStarFormation:EOS_temperature_norm_K"]
+)
 
 # Read reference metallicity
 EAGLE_Z = float(f["/Parameters"].attrs["EAGLEChemistry:init_abundance_metal"])
 
 # Read parameters of the entropy floor
-EAGLEfloor_Jeans_rho_norm = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_density_threshold_H_p_cm3"])
-EAGLEfloor_Jeans_temperature_norm_K = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_temperature_norm_K"])
-EAGLEfloor_Jeans_gamma_effective = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_gamma_effective"])
-EAGLEfloor_cool_rho_norm = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_density_threshold_H_p_cm3"])
-EAGLEfloor_cool_temperature_norm_K = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_temperature_norm_K"])
-EAGLEfloor_cool_gamma_effective = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_gamma_effective"])
+EAGLEfloor_Jeans_rho_norm = float(
+    f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_density_threshold_H_p_cm3"]
+)
+EAGLEfloor_Jeans_temperature_norm_K = float(
+    f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_temperature_norm_K"]
+)
+EAGLEfloor_Jeans_gamma_effective = float(
+    f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_gamma_effective"]
+)
+EAGLEfloor_cool_rho_norm = float(
+    f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_density_threshold_H_p_cm3"]
+)
+EAGLEfloor_cool_temperature_norm_K = float(
+    f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_temperature_norm_K"]
+)
+EAGLEfloor_cool_gamma_effective = float(
+    f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_gamma_effective"]
+)
 
 # Properties of the KS law
-KS_law_norm_cgs = KS_law_norm * Msun_in_cgs / ( 1e6 * pc_in_cgs**2 * year_in_cgs )
-gamma = 5./3.
+KS_law_norm_cgs = KS_law_norm * Msun_in_cgs / (1e6 * pc_in_cgs ** 2 * year_in_cgs)
+gamma = 5.0 / 3.0
 EOS_press_norm = k_in_cgs * EOS_temp_norm * EOS_density_norm
 
 # Star formation threshold
-SF_thresh = KS_thresh_norm * (EAGLE_Z / KS_thresh_Z0)**(KS_thresh_slope)
+SF_thresh = KS_thresh_norm * (EAGLE_Z / KS_thresh_Z0) ** (KS_thresh_slope)
 
 # Read gas properties
 gas_pos = f["/PartType0/Coordinates"][:, :]
 gas_mass = f["/PartType0/Masses"][:]
-gas_rho = f["/PartType0/Density"][:]
+gas_rho = f["/PartType0/Densities"][:]
 gas_T = f["/PartType0/Temperature"][:]
 gas_SFR = f["/PartType0/SFR"][:]
-gas_XH = f["/PartType0/ElementAbundance"][:, 0]
-gas_Z = f["/PartType0/Metallicity"][:]
-gas_hsml = f["/PartType0/SmoothingLength"][:]
+gas_XH = f["/PartType0/ElementMassFractions"][:, 0]
+gas_Z = f["/PartType0/Metallicities"][:]
+gas_hsml = f["/PartType0/SmoothingLengths"][:]
 gas_sSFR = gas_SFR / gas_mass
 
 # Read the Star properties
 stars_pos = f["/PartType4/Coordinates"][:, :]
 stars_BirthDensity = f["/PartType4/BirthDensity"][:]
 stars_BirthTime = f["/PartType4/BirthTime"][:]
-stars_XH = f["/PartType4/ElementAbundance"][:,0]
+stars_XH = f["/PartType4/ElementAbundance"][:, 0]
 
 # Centre the box
 gas_pos[:, 0] -= centre[0]
 gas_pos[:, 1] -= centre[1]
 gas_pos[:, 2] -= centre[2]
 
-stars_pos[:,0] -= centre[0]
-stars_pos[:,1] -= centre[1]
-stars_pos[:,2] -= centre[2]
+stars_pos[:, 0] -= centre[0]
+stars_pos[:, 1] -= centre[1]
+stars_pos[:, 2] -= centre[2]
 
 # Turn the mass into better units
 gas_mass *= unit_mass_in_cgs / Msun_in_cgs
@@ -132,9 +160,13 @@ stars_BirthDensity *= stars_XH
 
 # Equations of state
 eos_cool_rho = np.logspace(-5, 5, 1000)
-eos_cool_T = EAGLEfloor_cool_temperature_norm_K * (eos_cool_rho / EAGLEfloor_cool_rho_norm) ** ( EAGLEfloor_cool_gamma_effective - 1.0 )
+eos_cool_T = EAGLEfloor_cool_temperature_norm_K * (
+    eos_cool_rho / EAGLEfloor_cool_rho_norm
+) ** (EAGLEfloor_cool_gamma_effective - 1.0)
 eos_Jeans_rho = np.logspace(-1, 5, 1000)
-eos_Jeans_T = EAGLEfloor_Jeans_temperature_norm_K * (eos_Jeans_rho / EAGLEfloor_Jeans_rho_norm) ** (EAGLEfloor_Jeans_gamma_effective - 1.0 ) 
+eos_Jeans_T = EAGLEfloor_Jeans_temperature_norm_K * (
+    eos_Jeans_rho / EAGLEfloor_Jeans_rho_norm
+) ** (EAGLEfloor_Jeans_gamma_effective - 1.0)
 
 ########################################################################3
 
@@ -156,7 +188,15 @@ subplot(111, xscale="log", yscale="log")
 plot(eos_cool_rho, eos_cool_T, "k--", lw=0.6)
 plot(eos_Jeans_rho, eos_Jeans_T, "k--", lw=0.6)
 plot([SF_thresh, SF_thresh], [1, 1e10], "k:", lw=0.6)
-text(SF_thresh*0.9, 2e4, "$n_{\\rm H, thresh}=%.3f~{\\rm cm^{-3}}$"%SF_thresh, fontsize=8, rotation=90, ha="right", va="bottom")
+text(
+    SF_thresh * 0.9,
+    2e4,
+    "$n_{\\rm H, thresh}=%.3f~{\\rm cm^{-3}}$" % SF_thresh,
+    fontsize=8,
+    rotation=90,
+    ha="right",
+    va="bottom",
+)
 scatter(gas_nH[gas_SFR > 0.0], gas_T[gas_SFR > 0.0], s=0.2)
 xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
 ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=2)
@@ -188,37 +228,61 @@ star_mask = (
     & (stars_pos[:, 2] > -1.0)
 )
 
-stars_BirthDensity = stars_BirthDensity[star_mask] 
-#stars_BirthFlag = stars_BirthFlag[star_mask]
+stars_BirthDensity = stars_BirthDensity[star_mask]
+# stars_BirthFlag = stars_BirthFlag[star_mask]
 stars_BirthTime = stars_BirthTime[star_mask]
 
 # Histogram of the birth density
 figure()
 subplot(111, xscale="linear", yscale="linear")
-hist(np.log10(stars_BirthDensity),density=True,bins=20,range=[-2,5])
+hist(np.log10(stars_BirthDensity), density=True, bins=20, range=[-2, 5])
 xlabel("${\\rm Stellar~birth~density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
 ylabel("${\\rm Probability}$", labelpad=-7)
 savefig("BirthDensity.png", dpi=200)
 
 # Plot of the specific star formation rate in the galaxy
-rhos = 10**np.linspace(-1,np.log10(KS_high_den_thresh),100)
-rhoshigh = 10**np.linspace(np.log10(KS_high_den_thresh),5,100)
+rhos = 10 ** np.linspace(-1, np.log10(KS_high_den_thresh), 100)
+rhoshigh = 10 ** np.linspace(np.log10(KS_high_den_thresh), 5, 100)
 
-P_effective = EOS_press_norm * ( rhos / EOS_density_norm)**(EOS_gamma_effective)
-P_norm_high = EOS_press_norm * (KS_high_den_thresh  / EOS_density_norm)**(EOS_gamma_effective)
-sSFR = KS_law_norm_cgs * (Msun_p_pc2)**(-KS_law_slope) * (gamma/G_in_cgs * KS_gas_fraction *P_effective)**((KS_law_slope-1.)/2.)
-KS_law_norm_high_den_cgs = KS_law_norm_cgs * (Msun_p_pc2)**(-KS_law_slope) * (gamma/G_in_cgs * KS_gas_fraction * P_norm_high)**((KS_law_slope-1.)/2.)
-sSFR_high_den = KS_law_norm_high_den_cgs * ((rhoshigh/KS_high_den_thresh)**EOS_gamma_effective)**((KS_law_slope_high_den-1)/2.)
+P_effective = EOS_press_norm * (rhos / EOS_density_norm) ** (EOS_gamma_effective)
+P_norm_high = EOS_press_norm * (KS_high_den_thresh / EOS_density_norm) ** (
+    EOS_gamma_effective
+)
+sSFR = (
+    KS_law_norm_cgs
+    * (Msun_p_pc2) ** (-KS_law_slope)
+    * (gamma / G_in_cgs * KS_gas_fraction * P_effective) ** ((KS_law_slope - 1.0) / 2.0)
+)
+KS_law_norm_high_den_cgs = (
+    KS_law_norm_cgs
+    * (Msun_p_pc2) ** (-KS_law_slope)
+    * (gamma / G_in_cgs * KS_gas_fraction * P_norm_high) ** ((KS_law_slope - 1.0) / 2.0)
+)
+sSFR_high_den = KS_law_norm_high_den_cgs * (
+    (rhoshigh / KS_high_den_thresh) ** EOS_gamma_effective
+) ** ((KS_law_slope_high_den - 1) / 2.0)
 
 # density - sSFR plane
 figure()
 subplot(111)
-hist2d(np.log10(gas_nH), np.log10(gas_sSFR), bins=50,range=[[-1.5,5],[-.5,2.5]])
-plot(np.log10(rhos),np.log10(sSFR)+np.log10(year_in_cgs)+9.,'k--',label='sSFR low density EAGLE')
-plot(np.log10(rhoshigh),np.log10(sSFR_high_den)+np.log10(year_in_cgs)+9.,'k--',label='sSFR high density EAGLE')
+hist2d(np.log10(gas_nH), np.log10(gas_sSFR), bins=50, range=[[-1.5, 5], [-0.5, 2.5]])
+plot(
+    np.log10(rhos),
+    np.log10(sSFR) + np.log10(year_in_cgs) + 9.0,
+    "k--",
+    label="sSFR low density EAGLE",
+)
+plot(
+    np.log10(rhoshigh),
+    np.log10(sSFR_high_den) + np.log10(year_in_cgs) + 9.0,
+    "k--",
+    label="sSFR high density EAGLE",
+)
 xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=2)
 ylabel("${\\rm sSFR}~[{\\rm Gyr^{-1}}]$", labelpad=0)
-xticks([-1, 0, 1, 2, 3, 4], ["$10^{-1}$", "$10^0$", "$10^1$", "$10^2$", "$10^3$", "$10^4$"])
+xticks(
+    [-1, 0, 1, 2, 3, 4], ["$10^{-1}$", "$10^0$", "$10^1$", "$10^2$", "$10^3$", "$10^4$"]
+)
 yticks([0, 1, 2], ["$10^0$", "$10^1$", "$10^2$"])
 xlim(-1.4, 4.9)
 ylim(-0.5, 2.2)
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plot_box_evolution.py b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plot_box_evolution.py
index 67da3c390be1240323941b835e056dcd6e27feed..94f27c87ff46c48ffc6b0df8c3e02c7abb6df875 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plot_box_evolution.py
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plot_box_evolution.py
@@ -1,24 +1,25 @@
 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
- #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- # 
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- # 
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- # GNU General Public License for more details.
- # 
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program.  If not, see <http://www.gnu.org/licenses/>.
- # 
- ##############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+#                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
 
 import matplotlib
+
 matplotlib.use("Agg")
 from pylab import *
 from scipy import stats
@@ -27,40 +28,42 @@ import numpy as np
 import glob
 import os.path
 
-def find_indices(a,b):
-        result = np.zeros(len(b))
-        for i in range(len(b)):
-                result[i] = ((np.where(a == b[i]))[0])[0]
 
-        return result
+def find_indices(a, b):
+    result = np.zeros(len(b))
+    for i in range(len(b)):
+        result[i] = ((np.where(a == b[i]))[0])[0]
+
+    return result
 
 
 # Plot parameters
-params = {'axes.labelsize': 10,
-'axes.titlesize': 10,
-'font.size': 12,
-'legend.fontsize': 12,
-'xtick.labelsize': 10,
-'ytick.labelsize': 10,
-'text.usetex': True,
- 'figure.figsize' : (9.90,6.45),
-'figure.subplot.left'    : 0.1,
-'figure.subplot.right'   : 0.99,
-'figure.subplot.bottom'  : 0.1,
-'figure.subplot.top'     : 0.95,
-'figure.subplot.wspace'  : 0.2,
-'figure.subplot.hspace'  : 0.2,
-'lines.markersize' : 6,
-'lines.linewidth' : 3.,
-'text.latex.unicode': True
+params = {
+    "axes.labelsize": 10,
+    "axes.titlesize": 10,
+    "font.size": 12,
+    "legend.fontsize": 12,
+    "xtick.labelsize": 10,
+    "ytick.labelsize": 10,
+    "text.usetex": True,
+    "figure.figsize": (9.90, 6.45),
+    "figure.subplot.left": 0.1,
+    "figure.subplot.right": 0.99,
+    "figure.subplot.bottom": 0.1,
+    "figure.subplot.top": 0.95,
+    "figure.subplot.wspace": 0.2,
+    "figure.subplot.hspace": 0.2,
+    "lines.markersize": 6,
+    "lines.linewidth": 3.0,
+    "text.latex.unicode": True,
 }
 rcParams.update(params)
-rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
 
 
 # Number of snapshots and elements
-newest_snap_name = max(glob.glob('output_*.hdf5'), key=os.path.getctime)
-n_snapshots = int(newest_snap_name.replace('output_','').replace('.hdf5','')) + 1
+newest_snap_name = max(glob.glob("output_*.hdf5"), key=os.path.getctime)
+n_snapshots = int(newest_snap_name.replace("output_", "").replace(".hdf5", "")) + 1
 n_elements = 9
 
 # Read the simulation data
@@ -84,10 +87,10 @@ unit_energy_in_cgs = unit_mass_in_cgs * unit_vel_in_cgs * unit_vel_in_cgs
 unit_length_in_si = 0.01 * unit_length_in_cgs
 unit_mass_in_si = 0.001 * unit_mass_in_cgs
 unit_time_in_si = unit_time_in_cgs
-unit_density_in_cgs = unit_mass_in_cgs*unit_length_in_cgs**-3
-unit_pressure_in_cgs = unit_mass_in_cgs/unit_length_in_cgs*unit_time_in_cgs**-2
-unit_int_energy_in_cgs = unit_energy_in_cgs/unit_mass_in_cgs
-unit_entropy_in_cgs = unit_energy_in_cgs/unit_temp_in_cgs
+unit_density_in_cgs = unit_mass_in_cgs * unit_length_in_cgs ** -3
+unit_pressure_in_cgs = unit_mass_in_cgs / unit_length_in_cgs * unit_time_in_cgs ** -2
+unit_int_energy_in_cgs = unit_energy_in_cgs / unit_mass_in_cgs
+unit_entropy_in_cgs = unit_energy_in_cgs / unit_temp_in_cgs
 Gyr_in_cgs = 3.155e16
 Msun_in_cgs = 1.989e33
 
@@ -95,61 +98,111 @@ box_energy = zeros(n_snapshots)
 box_mass = zeros(n_snapshots)
 box_star_mass = zeros(n_snapshots)
 box_metal_mass = zeros(n_snapshots)
-element_mass = zeros((n_snapshots,n_elements))
+element_mass = zeros((n_snapshots, n_elements))
 t = zeros(n_snapshots)
 
 # Read data from snapshots
 for i in range(n_snapshots):
-	print("reading snapshot "+str(i))
-	# Read the simulation data
-	sim = h5py.File("output_%04d.hdf5"%i, "r")
-	t[i] = sim["/Header"].attrs["Time"][0]
-	#ids = sim["/PartType0/ParticleIDs"][:]
-	
-	masses = sim["/PartType0/Masses"][:]
-	box_mass[i] = np.sum(masses)
-	
-	star_masses = sim["/PartType4/Masses"][:]
-	box_star_mass[i] = np.sum(star_masses)
-	
-	metallicities = sim["/PartType0/Metallicity"][:]
-	box_metal_mass[i] = np.sum(metallicities * masses)
-	
-	internal_energies = sim["/PartType0/InternalEnergy"][:]
-	box_energy[i] = np.sum(masses * internal_energies)
+    print("reading snapshot " + str(i))
+    # Read the simulation data
+    sim = h5py.File("output_%04d.hdf5" % i, "r")
+    t[i] = sim["/Header"].attrs["Time"][0]
+    # ids = sim["/PartType0/ParticleIDs"][:]
+
+    masses = sim["/PartType0/Masses"][:]
+    box_mass[i] = np.sum(masses)
+
+    star_masses = sim["/PartType4/Masses"][:]
+    box_star_mass[i] = np.sum(star_masses)
+
+    metallicities = sim["/PartType0/Metallicities"][:]
+    box_metal_mass[i] = np.sum(metallicities * masses)
+
+    internal_energies = sim["/PartType0/InternalEnergies"][:]
+    box_energy[i] = np.sum(masses * internal_energies)
 
 # Plot the interesting quantities
 figure()
 
 # Box mass --------------------------------
 subplot(221)
-plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_mass[1:] - box_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', marker = "*", ms=0.5, label='swift')
+plot(
+    t[1:] * unit_time_in_cgs / Gyr_in_cgs,
+    (box_mass[1:] - box_mass[0]) * unit_mass_in_cgs / Msun_in_cgs,
+    linewidth=0.5,
+    color="k",
+    marker="*",
+    ms=0.5,
+    label="swift",
+)
 xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
 ylabel("Change in total gas particle mass (Msun)", labelpad=2)
-ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
 
 # Box metal mass --------------------------------
 subplot(222)
-plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_metal_mass[1:] - box_metal_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', ms=0.5, label='swift')
+plot(
+    t[1:] * unit_time_in_cgs / Gyr_in_cgs,
+    (box_metal_mass[1:] - box_metal_mass[0]) * unit_mass_in_cgs / Msun_in_cgs,
+    linewidth=0.5,
+    color="k",
+    ms=0.5,
+    label="swift",
+)
 xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
 ylabel("Change in total metal mass of gas particles (Msun)", labelpad=2)
-ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
 
 # Box energy --------------------------------
 subplot(223)
-plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_energy[1:] - box_energy[0])* unit_energy_in_cgs, linewidth=0.5, color='k', ms=0.5, label='swift')
+plot(
+    t[1:] * unit_time_in_cgs / Gyr_in_cgs,
+    (box_energy[1:] - box_energy[0]) * unit_energy_in_cgs,
+    linewidth=0.5,
+    color="k",
+    ms=0.5,
+    label="swift",
+)
 xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
 ylabel("Change in total energy of gas particles (erg)", labelpad=2)
-ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
 
 # Box mass --------------------------------
 subplot(224)
-plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_mass[1:] - box_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', marker = "*", ms=0.5, label='gas')
-plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_star_mass[1:] - box_star_mass[n_snapshots-1])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='r', marker = "*", ms=0.5, label='stars')
-plot(t[1:] * unit_time_in_cgs / Gyr_in_cgs, (box_star_mass[1:] - box_star_mass[n_snapshots-1] + box_mass[1:] - box_mass[0])* unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='g', marker = "*", ms=0.5, label='total')
+plot(
+    t[1:] * unit_time_in_cgs / Gyr_in_cgs,
+    (box_mass[1:] - box_mass[0]) * unit_mass_in_cgs / Msun_in_cgs,
+    linewidth=0.5,
+    color="k",
+    marker="*",
+    ms=0.5,
+    label="gas",
+)
+plot(
+    t[1:] * unit_time_in_cgs / Gyr_in_cgs,
+    (box_star_mass[1:] - box_star_mass[n_snapshots - 1])
+    * unit_mass_in_cgs
+    / Msun_in_cgs,
+    linewidth=0.5,
+    color="r",
+    marker="*",
+    ms=0.5,
+    label="stars",
+)
+plot(
+    t[1:] * unit_time_in_cgs / Gyr_in_cgs,
+    (box_star_mass[1:] - box_star_mass[n_snapshots - 1] + box_mass[1:] - box_mass[0])
+    * unit_mass_in_cgs
+    / Msun_in_cgs,
+    linewidth=0.5,
+    color="g",
+    marker="*",
+    ms=0.5,
+    label="total",
+)
 xlabel("${\\rm{Time}} (Gyr)$", labelpad=0)
 ylabel("Change in total gas particle mass (Msun)", labelpad=2)
-ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
 legend()
 
 savefig("box_evolution.png", dpi=200)
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/plotSolution.py b/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/plotSolution.py
index 73e4878e8e00a35fe19c359652be0d57153dea62..044ad2bc78958cbf669c7257122b1ff80a94ba1a 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/plotSolution.py
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/plotSolution.py
@@ -49,7 +49,7 @@ rcParams.update(params)
 rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
 
 snap = int(sys.argv[1])
-filename = "output_%.4d.hdf5"%snap
+filename = "output_%.4d.hdf5" % snap
 
 f = h5.File(filename, "r")
 
@@ -60,7 +60,7 @@ year_in_cgs = 3600.0 * 24 * 365.0
 Msun_in_cgs = 1.98848e33
 G_in_cgs = 6.67259e-8
 pc_in_cgs = 3.08567758e18
-Msun_p_pc2 = Msun_in_cgs / pc_in_cgs**2
+Msun_p_pc2 = Msun_in_cgs / pc_in_cgs ** 2
 
 # Gemoetry info
 boxsize = f["/Header"].attrs["BoxSize"]
@@ -72,57 +72,85 @@ unit_mass_in_cgs = f["/Units"].attrs["Unit mass in cgs (U_M)"]
 unit_time_in_cgs = f["/Units"].attrs["Unit time in cgs (U_t)"]
 
 # Calculate Gravitational constant in internal units
-G = G_in_cgs * ( unit_length_in_cgs**3 / unit_mass_in_cgs / unit_time_in_cgs**2)**(-1)
+G = G_in_cgs * (unit_length_in_cgs ** 3 / unit_mass_in_cgs / unit_time_in_cgs ** 2) ** (
+    -1
+)
 
 # Read parameters of the SF model
 KS_law_slope = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_exponent"])
 KS_law_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_normalisation"])
 KS_thresh_Z0 = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_Z0"])
 KS_thresh_slope = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_slope"])
-KS_thresh_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_norm_H_p_cm3"])
+KS_thresh_norm = float(
+    f["/Parameters"].attrs["EAGLEStarFormation:threshold_norm_H_p_cm3"]
+)
 KS_gas_fraction = float(f["/Parameters"].attrs["EAGLEStarFormation:gas_fraction"])
-KS_thresh_max_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:threshold_max_density_H_p_cm3"])
-KS_high_den_thresh = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_threshold_H_p_cm3"])
-KS_law_slope_high_den = float(f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_exponent"])
-EOS_gamma_effective = float(f["/Parameters"].attrs["EAGLEStarFormation:EOS_gamma_effective"])                           
-EOS_density_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:EOS_density_norm_H_p_cm3"])                           
-EOS_temp_norm = float(f["/Parameters"].attrs["EAGLEStarFormation:EOS_temperature_norm_K"])                           
+KS_thresh_max_norm = float(
+    f["/Parameters"].attrs["EAGLEStarFormation:threshold_max_density_H_p_cm3"]
+)
+KS_high_den_thresh = float(
+    f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_threshold_H_p_cm3"]
+)
+KS_law_slope_high_den = float(
+    f["/Parameters"].attrs["EAGLEStarFormation:KS_high_density_exponent"]
+)
+EOS_gamma_effective = float(
+    f["/Parameters"].attrs["EAGLEStarFormation:EOS_gamma_effective"]
+)
+EOS_density_norm = float(
+    f["/Parameters"].attrs["EAGLEStarFormation:EOS_density_norm_H_p_cm3"]
+)
+EOS_temp_norm = float(
+    f["/Parameters"].attrs["EAGLEStarFormation:EOS_temperature_norm_K"]
+)
 
 # Read reference metallicity
 EAGLE_Z = float(f["/Parameters"].attrs["EAGLEChemistry:init_abundance_metal"])
 
 # Read parameters of the entropy floor
-EAGLEfloor_Jeans_rho_norm = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_density_threshold_H_p_cm3"])
-EAGLEfloor_Jeans_temperature_norm_K = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_temperature_norm_K"])
-EAGLEfloor_Jeans_gamma_effective = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_gamma_effective"])
-EAGLEfloor_cool_rho_norm = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_density_threshold_H_p_cm3"])
-EAGLEfloor_cool_temperature_norm_K = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_temperature_norm_K"])
-EAGLEfloor_cool_gamma_effective = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_gamma_effective"])
+EAGLEfloor_Jeans_rho_norm = float(
+    f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_density_threshold_H_p_cm3"]
+)
+EAGLEfloor_Jeans_temperature_norm_K = float(
+    f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_temperature_norm_K"]
+)
+EAGLEfloor_Jeans_gamma_effective = float(
+    f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_gamma_effective"]
+)
+EAGLEfloor_cool_rho_norm = float(
+    f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_density_threshold_H_p_cm3"]
+)
+EAGLEfloor_cool_temperature_norm_K = float(
+    f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_temperature_norm_K"]
+)
+EAGLEfloor_cool_gamma_effective = float(
+    f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_gamma_effective"]
+)
 
 # Properties of the KS law
-KS_law_norm_cgs = KS_law_norm * Msun_in_cgs / ( 1e6 * pc_in_cgs**2 * year_in_cgs )
-gamma = 5./3.
+KS_law_norm_cgs = KS_law_norm * Msun_in_cgs / (1e6 * pc_in_cgs ** 2 * year_in_cgs)
+gamma = 5.0 / 3.0
 EOS_press_norm = k_in_cgs * EOS_temp_norm * EOS_density_norm
 
 # Star formation threshold
-SF_thresh = KS_thresh_norm * (EAGLE_Z / KS_thresh_Z0)**(KS_thresh_slope)
+SF_thresh = KS_thresh_norm * (EAGLE_Z / KS_thresh_Z0) ** (KS_thresh_slope)
 
 # Read gas properties
 gas_pos = f["/PartType0/Coordinates"][:, :]
 gas_mass = f["/PartType0/Masses"][:]
-gas_rho = f["/PartType0/Density"][:]
+gas_rho = f["/PartType0/Densities"][:]
 gas_T = f["/PartType0/Temperature"][:]
 gas_SFR = f["/PartType0/SFR"][:]
-gas_XH = f["/PartType0/ElementAbundance"][:, 0]
-gas_Z = f["/PartType0/Metallicity"][:]
-gas_hsml = f["/PartType0/SmoothingLength"][:]
+gas_XH = f["/PartType0/ElementMassFractions"][:, 0]
+gas_Z = f["/PartType0/Metallicities"][:]
+gas_hsml = f["/PartType0/SmoothingLengths"][:]
 gas_sSFR = gas_SFR / gas_mass
 
 # Read the Star properties
 stars_pos = f["/PartType4/Coordinates"][:, :]
 stars_BirthDensity = f["/PartType4/BirthDensity"][:]
 stars_BirthTime = f["/PartType4/BirthTime"][:]
-stars_XH = f["/PartType4/ElementAbundance"][:,0]
+stars_XH = f["/PartType4/ElementAbundance"][:, 0]
 stars_BirthTemperature = f["/PartType4/BirthTemperature"][:]
 
 # Centre the box
@@ -130,9 +158,9 @@ gas_pos[:, 0] -= centre[0]
 gas_pos[:, 1] -= centre[1]
 gas_pos[:, 2] -= centre[2]
 
-stars_pos[:,0] -= centre[0]
-stars_pos[:,1] -= centre[1]
-stars_pos[:,2] -= centre[2]
+stars_pos[:, 0] -= centre[0]
+stars_pos[:, 1] -= centre[1]
+stars_pos[:, 2] -= centre[2]
 
 # Turn the mass into better units
 gas_mass *= unit_mass_in_cgs / Msun_in_cgs
@@ -156,9 +184,13 @@ stars_BirthDensity *= stars_XH
 
 # Equations of state
 eos_cool_rho = np.logspace(-5, 5, 1000)
-eos_cool_T = EAGLEfloor_cool_temperature_norm_K * (eos_cool_rho / EAGLEfloor_cool_rho_norm) ** ( EAGLEfloor_cool_gamma_effective - 1.0 )
+eos_cool_T = EAGLEfloor_cool_temperature_norm_K * (
+    eos_cool_rho / EAGLEfloor_cool_rho_norm
+) ** (EAGLEfloor_cool_gamma_effective - 1.0)
 eos_Jeans_rho = np.logspace(-1, 5, 1000)
-eos_Jeans_T = EAGLEfloor_Jeans_temperature_norm_K * (eos_Jeans_rho / EAGLEfloor_Jeans_rho_norm) ** (EAGLEfloor_Jeans_gamma_effective - 1.0 ) 
+eos_Jeans_T = EAGLEfloor_Jeans_temperature_norm_K * (
+    eos_Jeans_rho / EAGLEfloor_Jeans_rho_norm
+) ** (EAGLEfloor_Jeans_gamma_effective - 1.0)
 
 ########################################################################3
 
@@ -180,7 +212,15 @@ subplot(111, xscale="log", yscale="log")
 plot(eos_cool_rho, eos_cool_T, "k--", lw=0.6)
 plot(eos_Jeans_rho, eos_Jeans_T, "k--", lw=0.6)
 plot([SF_thresh, SF_thresh], [1, 1e10], "k:", lw=0.6)
-text(SF_thresh*0.9, 2e4, "$n_{\\rm H, thresh}=%.3f~{\\rm cm^{-3}}$"%SF_thresh, fontsize=8, rotation=90, ha="right", va="bottom")
+text(
+    SF_thresh * 0.9,
+    2e4,
+    "$n_{\\rm H, thresh}=%.3f~{\\rm cm^{-3}}$" % SF_thresh,
+    fontsize=8,
+    rotation=90,
+    ha="right",
+    va="bottom",
+)
 scatter(gas_nH[gas_SFR > 0.0], gas_T[gas_SFR > 0.0], s=0.2)
 xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
 ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=2)
@@ -199,66 +239,92 @@ star_mask = (
     & (stars_pos[:, 2] > -1.0)
 )
 
-stars_BirthDensity = stars_BirthDensity[star_mask] 
-#stars_BirthFlag = stars_BirthFlag[star_mask]
+stars_BirthDensity = stars_BirthDensity[star_mask]
+# stars_BirthFlag = stars_BirthFlag[star_mask]
 stars_BirthTime = stars_BirthTime[star_mask]
 
 # Histogram of the birth density
 figure()
 subplot(111, xscale="linear", yscale="linear")
-hist(np.log10(stars_BirthDensity),density=True,bins=20,range=[-2,5])
+hist(np.log10(stars_BirthDensity), density=True, bins=20, range=[-2, 5])
 xlabel("${\\rm Stellar~birth~density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
 ylabel("${\\rm Probability}$", labelpad=3)
 savefig("BirthDensity.png", dpi=200)
 
-# Histogram of the birth temperature 
+# Histogram of the birth temperature
 figure()
 subplot(111, xscale="linear", yscale="linear")
-hist(np.log10(stars_BirthTemperature),density=True,bins=20,range=[3.5,5.0])
+hist(np.log10(stars_BirthTemperature), density=True, bins=20, range=[3.5, 5.0])
 xlabel("${\\rm Stellar~birth~temperature}~[{\\rm K}]$", labelpad=0)
 ylabel("${\\rm Probability}$", labelpad=3)
 savefig("BirthTemperature.png", dpi=200)
 
 # Plot of the specific star formation rate in the galaxy
-rhos = 10**np.linspace(-1,np.log10(KS_high_den_thresh),100)
-rhoshigh = 10**np.linspace(np.log10(KS_high_den_thresh),5,100)
+rhos = 10 ** np.linspace(-1, np.log10(KS_high_den_thresh), 100)
+rhoshigh = 10 ** np.linspace(np.log10(KS_high_den_thresh), 5, 100)
 
-P_effective = EOS_press_norm * ( rhos / EOS_density_norm)**(EOS_gamma_effective)
-P_norm_high = EOS_press_norm * (KS_high_den_thresh  / EOS_density_norm)**(EOS_gamma_effective)
-sSFR = KS_law_norm_cgs * (Msun_p_pc2)**(-KS_law_slope) * (gamma/G_in_cgs * KS_gas_fraction *P_effective)**((KS_law_slope-1.)/2.)
-KS_law_norm_high_den_cgs = KS_law_norm_cgs * (Msun_p_pc2)**(-KS_law_slope) * (gamma/G_in_cgs * KS_gas_fraction * P_norm_high)**((KS_law_slope-1.)/2.)
-sSFR_high_den = KS_law_norm_high_den_cgs * ((rhoshigh/KS_high_den_thresh)**EOS_gamma_effective)**((KS_law_slope_high_den-1)/2.)
+P_effective = EOS_press_norm * (rhos / EOS_density_norm) ** (EOS_gamma_effective)
+P_norm_high = EOS_press_norm * (KS_high_den_thresh / EOS_density_norm) ** (
+    EOS_gamma_effective
+)
+sSFR = (
+    KS_law_norm_cgs
+    * (Msun_p_pc2) ** (-KS_law_slope)
+    * (gamma / G_in_cgs * KS_gas_fraction * P_effective) ** ((KS_law_slope - 1.0) / 2.0)
+)
+KS_law_norm_high_den_cgs = (
+    KS_law_norm_cgs
+    * (Msun_p_pc2) ** (-KS_law_slope)
+    * (gamma / G_in_cgs * KS_gas_fraction * P_norm_high) ** ((KS_law_slope - 1.0) / 2.0)
+)
+sSFR_high_den = KS_law_norm_high_den_cgs * (
+    (rhoshigh / KS_high_den_thresh) ** EOS_gamma_effective
+) ** ((KS_law_slope_high_den - 1) / 2.0)
 
 # density - sSFR plane
 figure()
 subplot(111)
-hist2d(np.log10(gas_nH), np.log10(gas_sSFR), bins=50,range=[[-1.5,5],[-.5,2.5]])
-plot(np.log10(rhos),np.log10(sSFR)+np.log10(year_in_cgs)+9.,'k--',label='sSFR low density EAGLE')
-plot(np.log10(rhoshigh),np.log10(sSFR_high_den)+np.log10(year_in_cgs)+9.,'k--',label='sSFR high density EAGLE')
+hist2d(np.log10(gas_nH), np.log10(gas_sSFR), bins=50, range=[[-1.5, 5], [-0.5, 2.5]])
+plot(
+    np.log10(rhos),
+    np.log10(sSFR) + np.log10(year_in_cgs) + 9.0,
+    "k--",
+    label="sSFR low density EAGLE",
+)
+plot(
+    np.log10(rhoshigh),
+    np.log10(sSFR_high_den) + np.log10(year_in_cgs) + 9.0,
+    "k--",
+    label="sSFR high density EAGLE",
+)
 xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=2)
 ylabel("${\\rm sSFR}~[{\\rm Gyr^{-1}}]$", labelpad=0)
-xticks([-1, 0, 1, 2, 3, 4], ["$10^{-1}$", "$10^0$", "$10^1$", "$10^2$", "$10^3$", "$10^4$"])
+xticks(
+    [-1, 0, 1, 2, 3, 4], ["$10^{-1}$", "$10^0$", "$10^1$", "$10^2$", "$10^3$", "$10^4$"]
+)
 yticks([0, 1, 2], ["$10^0$", "$10^1$", "$10^2$"])
 xlim(-1.4, 4.9)
 ylim(-0.5, 2.2)
 savefig("density-sSFR.png", dpi=200)
 
-SFR_low = 10**(np.log10(sSFR)+np.log10(year_in_cgs)+np.log10(median_gas_mass))
-SFR_high = 10**(np.log10(sSFR_high_den)+np.log10(year_in_cgs)+np.log10(median_gas_mass))
-SFR_low_min = np.floor(np.log10(.75*np.min(SFR_low)))
-SFR_high_max = np.ceil(np.log10(1.25*np.max(SFR_high)))
+SFR_low = 10 ** (np.log10(sSFR) + np.log10(year_in_cgs) + np.log10(median_gas_mass))
+SFR_high = 10 ** (
+    np.log10(sSFR_high_den) + np.log10(year_in_cgs) + np.log10(median_gas_mass)
+)
+SFR_low_min = np.floor(np.log10(0.75 * np.min(SFR_low)))
+SFR_high_max = np.ceil(np.log10(1.25 * np.max(SFR_high)))
 
 # 3D Density vs SFR
 rcParams.update({"figure.subplot.left": 0.18})
 figure()
 subplot(111, xscale="log", yscale="log")
 scatter(gas_nH, gas_SFR, s=0.2)
-plot(rhos,SFR_low,'k--',lw=1,label='SFR low density EAGLE')
-plot(rhoshigh,SFR_high,'k--',lw=1,label='SFR high density EAGLE')
+plot(rhos, SFR_low, "k--", lw=1, label="SFR low density EAGLE")
+plot(rhoshigh, SFR_high, "k--", lw=1, label="SFR high density EAGLE")
 xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
 ylabel("${\\rm SFR}~[{\\rm M_\\odot~\\cdot~yr^{-1}}]$", labelpad=2)
 xlim(1e-2, 1e5)
-ylim(10**SFR_low_min, 10**(SFR_high_max+0.1))
+ylim(10 ** SFR_low_min, 10 ** (SFR_high_max + 0.1))
 savefig("rho_SFR.png", dpi=200)
 rcParams.update({"figure.subplot.left": 0.15})
 ########################################################################3
diff --git a/examples/SantaBarbara/SantaBarbara-256/plotTempEvolution.py b/examples/SantaBarbara/SantaBarbara-256/plotTempEvolution.py
index dab4b2c90a7b751c8d143ed38c614473c951988a..63e46ccaee7be9ea18090e13ae15bb0a1fae4bef 100644
--- a/examples/SantaBarbara/SantaBarbara-256/plotTempEvolution.py
+++ b/examples/SantaBarbara/SantaBarbara-256/plotTempEvolution.py
@@ -128,7 +128,7 @@ for i in range(n_snapshots):
     z[i] = sim["/Cosmology"].attrs["Redshift"][0]
     a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
 
-    u = sim["/PartType0/InternalEnergy"][:]
+    u = sim["/PartType0/InternalEnergies"][:]
 
     # Compute the temperature
     u *= unit_length_in_si ** 2 / unit_time_in_si ** 2
diff --git a/examples/SantaBarbara/SantaBarbara-256/rhoTPlot.py b/examples/SantaBarbara/SantaBarbara-256/rhoTPlot.py
index c290268eaa548e188bb652104ea9e726ea88a267..3bcf01d2a49bc1c53b243ffcff12359201d26d87 100644
--- a/examples/SantaBarbara/SantaBarbara-256/rhoTPlot.py
+++ b/examples/SantaBarbara/SantaBarbara-256/rhoTPlot.py
@@ -28,10 +28,10 @@ def get_data(filename):
 
     data = SWIFTDataset(filename)
 
-    data.gas.density.convert_to_units(mh / (cm ** 3))
-    data.gas.temperature.convert_to_cgs()
+    data.gas.densities.convert_to_units(mh / (cm ** 3))
+    data.gas.temperatures.convert_to_cgs()
 
-    return data.gas.density, data.gas.temperature
+    return data.gas.densities, data.gas.temperatures
 
 
 def make_hist(filename, density_bounds, temperature_bounds, bins):
@@ -155,10 +155,8 @@ def make_movie(args, density_bounds, temperature_bounds, bins):
     def format_metadata(metadata: SWIFTMetadata):
         t = metadata.t * units.units["Unit time in cgs (U_t)"]
         t.convert_to_units(Gyr)
-        
-        x = "$a$: {:2.2f}\n$z$: {:2.2f}\n$t$: {:2.2f}".format(
-            metadata.a, metadata.z, t
-        )
+
+        x = "$a$: {:2.2f}\n$z$: {:2.2f}\n$t$: {:2.2f}".format(metadata.a, metadata.z, t)
 
         return x
 
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotRhoT.py b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotRhoT.py
index 4ba8ad66daca1d9614be8917a77407dd99209dea..4f02213ec2a66700d28ad5f8e57e00c30f3019d7 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotRhoT.py
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotRhoT.py
@@ -112,8 +112,8 @@ def T(u, H_frac=H_mass_fraction, T_trans=H_transition_temp):
     return ret
 
 
-rho = sim["/PartType0/Density"][:]
-u = sim["/PartType0/InternalEnergy"][:]
+rho = sim["/PartType0/Densities"][:]
+u = sim["/PartType0/InternalEnergies"][:]
 
 # Compute the temperature
 u *= unit_length_in_si ** 2 / unit_time_in_si ** 2
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotTempEvolution.py b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotTempEvolution.py
index a3458ac1598e5657f3f597dfb10b36a7a641e68f..1e8cf9ea1082372d8e395c352f908c7ce693d99f 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotTempEvolution.py
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/plotTempEvolution.py
@@ -126,7 +126,7 @@ for i in range(n_snapshots):
     z[i] = sim["/Cosmology"].attrs["Redshift"][0]
     a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
 
-    u = sim["/PartType0/InternalEnergy"][:]
+    u = sim["/PartType0/InternalEnergies"][:]
 
     # Compute the temperature
     u *= (unit_length_in_si**2 / unit_time_in_si**2)
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml
index b108290fcd146461827d5858742bd0971ac66945..5050f9e05e278df25b050c10d9897e29d4e21804 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml
@@ -85,6 +85,16 @@ EAGLEChemistry:
   init_abundance_Silicon:   0.0
   init_abundance_Iron:      0.0
 
+EAGLEEntropyFloor:
+  Jeans_density_threshold_H_p_cm3: 0.1       # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+  Jeans_over_density_threshold:    10.       # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
+  Jeans_temperature_norm_K:        8000      # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
+  Jeans_gamma_effective:           1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
+  Cool_density_threshold_H_p_cm3:  1e-5      # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+  Cool_over_density_threshold:     10.       # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
+  Cool_temperature_norm_K:         8000      # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
+  Cool_gamma_effective:            1.        # Slope the of the EAGLE Cool limiter entropy floor
+  
 # Cooling with Grackle 3.0
 GrackleCooling:
   CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
@@ -99,3 +109,6 @@ GrackleCooling:
 
 GearChemistry:
   InitialMetallicity: 0.01295
+
+GEARPressureFloor:
+  Jeans_factor: 10.       # Number of particles required to suppose a resolved clump and avoid the pressure floor.
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/plotTempEvolution.py b/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/plotTempEvolution.py
index aa6c5df5fe5ff5c7d0944a45bb11344f70c57844..d707f70450471f2d2fc589dbc382366280e0e7f3 100644
--- a/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/plotTempEvolution.py
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_hydro/plotTempEvolution.py
@@ -123,7 +123,7 @@ for i in range(n_snapshots):
     z[i] = sim["/Cosmology"].attrs["Redshift"][0]
     a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
 
-    u = sim["/PartType0/InternalEnergy"][:]
+    u = sim["/PartType0/InternalEnergies"][:]
 
     # Compute the temperature
     u *= (unit_length_in_si**2 / unit_time_in_si**2)
diff --git a/examples/SubgridTests/BlackHoleSwallowing/check_masses.py b/examples/SubgridTests/BlackHoleSwallowing/check_masses.py
index c7f1d7b2c1f90efa285c13c75f0e5243f36e49ea..a5b55ed00df0851f989858ddffd00ea34df88a27 100644
--- a/examples/SubgridTests/BlackHoleSwallowing/check_masses.py
+++ b/examples/SubgridTests/BlackHoleSwallowing/check_masses.py
@@ -66,5 +66,5 @@ for i in range(np.size(ids_removed)):
     result = np.where(ids_gas == ids_removed)
     print result
 
-#rho_gas = f["/PartType0/Density"][:]
+#rho_gas = f["/PartType0/Densities"][:]
 #print np.mean(rho_gas), np.std(rho_gas)
diff --git a/examples/SubgridTests/SmoothedMetallicity/plotSolution.py b/examples/SubgridTests/SmoothedMetallicity/plotSolution.py
index e5bca3dfb7fe1e43c836733894c9e297cdd468ca..068fe5378e19c34ee8a68398f4e0ed096d0982e0 100644
--- a/examples/SubgridTests/SmoothedMetallicity/plotSolution.py
+++ b/examples/SubgridTests/SmoothedMetallicity/plotSolution.py
@@ -3,20 +3,20 @@
 # This file is part of SWIFT.
 # Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
 #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
-# 
+#
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published
 # by the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
-# 
+#
 # This program is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 # GNU General Public License for more details.
-# 
+#
 # You should have received a copy of the GNU Lesser General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
-# 
+#
 ##############################################################################
 
 # Computes the analytical solution of the 3D Smoothed Metallicity example.
@@ -25,12 +25,13 @@ import h5py
 import sys
 import numpy as np
 import matplotlib
+
 matplotlib.use("Agg")
 import matplotlib.pyplot as plt
 
 # Parameters
-low_metal = -6     # low metal abundance
-high_metal = -5    # High metal abundance
+low_metal = -6  # low metal abundance
+high_metal = -5  # High metal abundance
 sigma_metal = 0.1  # relative standard deviation for Z
 
 Nelem = 9
@@ -44,27 +45,27 @@ high_metal = [high_metal] * Nelem + np.linspace(0, 3, Nelem)
 
 # Plot parameters
 params = {
-    'axes.labelsize': 10,
-    'axes.titlesize': 10,
-    'font.size': 12,
-    'legend.fontsize': 12,
-    'xtick.labelsize': 10,
-    'ytick.labelsize': 10,
-    'text.usetex': True,
-    'figure.figsize': (9.90, 6.45),
-    'figure.subplot.left': 0.045,
-    'figure.subplot.right': 0.99,
-    'figure.subplot.bottom': 0.05,
-    'figure.subplot.top': 0.99,
-    'figure.subplot.wspace': 0.15,
-    'figure.subplot.hspace': 0.12,
-    'lines.markersize': 6,
-    'lines.linewidth': 3.,
-    'text.latex.unicode': True
+    "axes.labelsize": 10,
+    "axes.titlesize": 10,
+    "font.size": 12,
+    "legend.fontsize": 12,
+    "xtick.labelsize": 10,
+    "ytick.labelsize": 10,
+    "text.usetex": True,
+    "figure.figsize": (9.90, 6.45),
+    "figure.subplot.left": 0.045,
+    "figure.subplot.right": 0.99,
+    "figure.subplot.bottom": 0.05,
+    "figure.subplot.top": 0.99,
+    "figure.subplot.wspace": 0.15,
+    "figure.subplot.hspace": 0.12,
+    "lines.markersize": 6,
+    "lines.linewidth": 3.0,
+    "text.latex.unicode": True,
 }
 
 plt.rcParams.update(params)
-plt.rc('font', **{'family': 'sans-serif', 'sans-serif': ['Times']})
+plt.rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
 
 
 snap = int(sys.argv[1])
@@ -83,18 +84,17 @@ git = sim["Code"].attrs["Git Revision"]
 
 pos = sim["/PartType0/Coordinates"][:, :]
 d = pos[:, 0] - boxSize / 2
-smooth_metal = sim["/PartType0/SmoothedElementAbundance"][:, :]
-metal = sim["/PartType0/ElementAbundance"][:, :]
-h = sim["/PartType0/SmoothingLength"][:]
+smooth_metal = sim["/PartType0/SmoothedElementMassFractions"][:, :]
+metal = sim["/PartType0/ElementMassFractions"][:, :]
+h = sim["/PartType0/SmoothingLengths"][:]
 h = np.mean(h)
 
-if (Nelem != metal.shape[1]):
-    print("Unexpected number of element, please check makeIC.py"
-          " and plotSolution.py")
+if Nelem != metal.shape[1]:
+    print("Unexpected number of element, please check makeIC.py" " and plotSolution.py")
     exit(1)
 
 N = 1000
-d_a = np.linspace(-boxSize / 2., boxSize / 2., N)
+d_a = np.linspace(-boxSize / 2.0, boxSize / 2.0, N)
 
 # Now, work our the solution....
 
@@ -142,14 +142,14 @@ def calc_a(d, high_metal, low_metal, std_dev, h):
         m = (high_metal[i] - low_metal[i]) / (2.0 * h)
         c = (high_metal[i] + low_metal[i]) / 2.0
         # compute left linear part
-        s = d < - boxSize / 2.0 + h
-        a[s, i] = - m * (d[s] + boxSize / 2.0) + c
+        s = d < -boxSize / 2.0 + h
+        a[s, i] = -m * (d[s] + boxSize / 2.0) + c
         # compute middle linear part
         s = np.logical_and(d >= -h, d <= h)
         a[s, i] = m * d[s] + c
         # compute right linear part
         s = d > boxSize / 2.0 - h
-        a[s, i] = - m * (d[s] - boxSize / 2.0) + c
+        a[s, i] = -m * (d[s] - boxSize / 2.0) + c
 
     sigma[:, :, 0] = a * (1 + std_dev)
     sigma[:, :, 1] = a * (1 - std_dev)
@@ -165,7 +165,7 @@ plt.figure()
 # Metallicity --------------------------------
 plt.subplot(221)
 for e in range(Nelem):
-    plt.plot(metal[:, e], smooth_metal[:, e], '.', ms=0.5, alpha=0.2)
+    plt.plot(metal[:, e], smooth_metal[:, e], ".", ms=0.5, alpha=0.2)
 
 xmin, xmax = metal.min(), metal.max()
 ymin, ymax = smooth_metal.min(), smooth_metal.max()
@@ -178,27 +178,28 @@ plt.ylabel("${\\rm{Smoothed~Metallicity}}~Z_\\textrm{sm}$", labelpad=0)
 # Metallicity --------------------------------
 e = 0
 plt.subplot(223)
-plt.plot(d, smooth_metal[:, e], '.', color='r', ms=0.5, alpha=0.2)
-plt.plot(d_a, sol[:, e], '--', color='b', alpha=0.8, lw=1.2)
-plt.fill_between(d_a, sig[:, e, 0], sig[:, e, 1], facecolor="b",
-                 interpolate=True, alpha=0.5)
+plt.plot(d, smooth_metal[:, e], ".", color="r", ms=0.5, alpha=0.2)
+plt.plot(d_a, sol[:, e], "--", color="b", alpha=0.8, lw=1.2)
+plt.fill_between(
+    d_a, sig[:, e, 0], sig[:, e, 1], facecolor="b", interpolate=True, alpha=0.5
+)
 plt.xlabel("${\\rm{Distance}}~r$", labelpad=0)
 plt.ylabel("${\\rm{Smoothed~Metallicity}}~Z_\\textrm{sm}$", labelpad=0)
 plt.xlim(-0.5, 0.5)
-plt.ylim(low_metal[e]-1, high_metal[e]+1)
+plt.ylim(low_metal[e] - 1, high_metal[e] + 1)
 
 # Information -------------------------------------
 plt.subplot(222, frameon=False)
 
-plt.text(-0.49, 0.9, "Smoothed Metallicity in 3D at $t=%.2f$" % time,
-         fontsize=10)
-plt.plot([-0.49, 0.1], [0.82, 0.82], 'k-', lw=1)
+plt.text(-0.49, 0.9, "Smoothed Metallicity in 3D at $t=%.2f$" % time, fontsize=10)
+plt.plot([-0.49, 0.1], [0.82, 0.82], "k-", lw=1)
 plt.text(-0.49, 0.7, "$\\textsc{Swift}$ %s" % git, fontsize=10)
 plt.text(-0.49, 0.6, scheme, fontsize=10)
 plt.text(-0.49, 0.5, kernel, fontsize=10)
 plt.text(-0.49, 0.4, chemistry + "'s Chemistry", fontsize=10)
-plt.text(-0.49, 0.3, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
-         fontsize=10)
+plt.text(
+    -0.49, 0.3, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta), fontsize=10
+)
 plt.xlim(-0.5, 0.5)
 plt.ylim(0, 1)
 plt.xticks([])
diff --git a/examples/SubgridTests/StellarEvolution/check_continuous_heating.py b/examples/SubgridTests/StellarEvolution/check_continuous_heating.py
index f3c1b5d7fd682d914f2dbc05259c2dab0baf1e32..b5940eba43c89b8af4a883cc8f7022e33293b869 100644
--- a/examples/SubgridTests/StellarEvolution/check_continuous_heating.py
+++ b/examples/SubgridTests/StellarEvolution/check_continuous_heating.py
@@ -93,7 +93,7 @@ for i in range(n_snapshots):
 	sim = h5py.File("stellar_evolution_%04d.hdf5"%i, "r")
 	print('reading snapshot '+str(i))
 	masses[:,i] = sim["/PartType0/Masses"]
-	internal_energy[:,i] = sim["/PartType0/InternalEnergy"]
+	internal_energy[:,i] = sim["/PartType0/InternalEnergies"]
 	velocity_parts[:,:,i] = sim["/PartType0/Velocities"]
 	time[i] = sim["/Header"].attrs["Time"][0]
 
diff --git a/examples/SubgridTests/StellarEvolution/check_stellar_evolution.py b/examples/SubgridTests/StellarEvolution/check_stellar_evolution.py
index 02c1e9343de7b58cfddc8dee3bf0215a4b80ccf4..5680eb4d64f29ab32831d995e31ae9c27de82a71 100644
--- a/examples/SubgridTests/StellarEvolution/check_stellar_evolution.py
+++ b/examples/SubgridTests/StellarEvolution/check_stellar_evolution.py
@@ -1,4 +1,5 @@
 import matplotlib
+
 matplotlib.use("Agg")
 from pylab import *
 import h5py
@@ -7,32 +8,35 @@ import numpy as np
 import glob
 
 # Number of snapshots and elements
-newest_snap_name = max(glob.glob('stellar_evolution_*.hdf5'), key=os.path.getctime)
-n_snapshots = int(newest_snap_name.replace('stellar_evolution_','').replace('.hdf5','')) + 1
+newest_snap_name = max(glob.glob("stellar_evolution_*.hdf5"), key=os.path.getctime)
+n_snapshots = (
+    int(newest_snap_name.replace("stellar_evolution_", "").replace(".hdf5", "")) + 1
+)
 n_elements = 9
 
 # Plot parameters
-params = {'axes.labelsize': 10,
-'axes.titlesize': 10,
-'font.size': 9,
-'legend.fontsize': 9,
-'xtick.labelsize': 10,
-'ytick.labelsize': 10,
-'text.usetex': True,
- 'figure.figsize' : (3.15,3.15),
-'figure.subplot.left'    : 0.3,
-'figure.subplot.right'   : 0.99,
-'figure.subplot.bottom'  : 0.18,
-'figure.subplot.top'     : 0.92,
-'figure.subplot.wspace'  : 0.21,
-'figure.subplot.hspace'  : 0.19,
-'lines.markersize' : 6,
-'lines.linewidth' : 2.,
-'text.latex.unicode': True
+params = {
+    "axes.labelsize": 10,
+    "axes.titlesize": 10,
+    "font.size": 9,
+    "legend.fontsize": 9,
+    "xtick.labelsize": 10,
+    "ytick.labelsize": 10,
+    "text.usetex": True,
+    "figure.figsize": (3.15, 3.15),
+    "figure.subplot.left": 0.3,
+    "figure.subplot.right": 0.99,
+    "figure.subplot.bottom": 0.18,
+    "figure.subplot.top": 0.92,
+    "figure.subplot.wspace": 0.21,
+    "figure.subplot.hspace": 0.19,
+    "lines.markersize": 6,
+    "lines.linewidth": 2.0,
+    "text.latex.unicode": True,
 }
 
 rcParams.update(params)
-rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
 
 # Read the simulation data
 sim = h5py.File("stellar_evolution_0000.hdf5", "r")
@@ -43,7 +47,9 @@ neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"][0]
 eta = sim["/HydroScheme"].attrs["Kernel eta"][0]
 alpha = sim["/HydroScheme"].attrs["Alpha viscosity"][0]
 H_mass_fraction = sim["/HydroScheme"].attrs["Hydrogen mass fraction"][0]
-H_transition_temp = sim["/HydroScheme"].attrs["Hydrogen ionization transition temperature"][0]
+H_transition_temp = sim["/HydroScheme"].attrs[
+    "Hydrogen ionization transition temperature"
+][0]
 T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
 T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
 git = sim["Code"].attrs["Git Revision"]
@@ -68,136 +74,245 @@ n_parts = sim["/Header"].attrs["NumPart_Total"][0]
 n_sparts = sim["/Header"].attrs["NumPart_Total"][4]
 
 # Declare arrays for data
-masses = zeros((n_parts,n_snapshots))
-star_masses = zeros((n_sparts,n_snapshots))
-mass_from_AGB = zeros((n_parts,n_snapshots))
-metal_mass_frac_from_AGB = zeros((n_parts,n_snapshots))
-mass_from_SNII = zeros((n_parts,n_snapshots))
-metal_mass_frac_from_SNII = zeros((n_parts,n_snapshots))
-mass_from_SNIa = zeros((n_parts,n_snapshots))
-metal_mass_frac_from_SNIa = zeros((n_parts,n_snapshots))
-iron_mass_frac_from_SNIa = zeros((n_parts,n_snapshots))
-metallicity = zeros((n_parts,n_snapshots))
-abundances = zeros((n_parts,n_elements,n_snapshots))
-internal_energy = zeros((n_parts,n_snapshots))
-coord_parts = zeros((n_parts,3))
-velocity_parts = zeros((n_parts,3,n_snapshots))
+masses = zeros((n_parts, n_snapshots))
+star_masses = zeros((n_sparts, n_snapshots))
+mass_from_AGB = zeros((n_parts, n_snapshots))
+metal_mass_frac_from_AGB = zeros((n_parts, n_snapshots))
+mass_from_SNII = zeros((n_parts, n_snapshots))
+metal_mass_frac_from_SNII = zeros((n_parts, n_snapshots))
+mass_from_SNIa = zeros((n_parts, n_snapshots))
+metal_mass_frac_from_SNIa = zeros((n_parts, n_snapshots))
+iron_mass_frac_from_SNIa = zeros((n_parts, n_snapshots))
+metallicity = zeros((n_parts, n_snapshots))
+abundances = zeros((n_parts, n_elements, n_snapshots))
+internal_energy = zeros((n_parts, n_snapshots))
+coord_parts = zeros((n_parts, 3))
+velocity_parts = zeros((n_parts, 3, n_snapshots))
 coord_sparts = zeros(3)
 time = zeros(n_snapshots)
 
 # Read fields we are checking from snapshots
 for i in range(n_snapshots):
-	sim = h5py.File("stellar_evolution_%04d.hdf5"%i, "r")
-	print('reading snapshot '+str(i))
-	abundances[:,:,i] = sim["/PartType0/ElementAbundance"]
-	metallicity[:,i] = sim["/PartType0/Metallicity"]
-	masses[:,i] = sim["/PartType0/Masses"]
-	star_masses[:,i] = sim["/PartType4/Masses"]
-	mass_from_AGB[:,i] = sim["/PartType0/TotalMassFromAGB"]
-	metal_mass_frac_from_AGB[:,i] = sim["/PartType0/MetalMassFracFromAGB"]
-	mass_from_SNII[:,i] = sim["/PartType0/TotalMassFromSNII"]
-	metal_mass_frac_from_SNII[:,i] = sim["/PartType0/MetalMassFracFromSNII"]
-	mass_from_SNIa[:,i] = sim["/PartType0/TotalMassFromSNIa"]
-	metal_mass_frac_from_SNIa[:,i] = sim["/PartType0/MetalMassFracFromSNIa"]
-	iron_mass_frac_from_SNIa[:,i] = sim["/PartType0/IronMassFracFromSNIa"]
-	internal_energy[:,i] = sim["/PartType0/InternalEnergy"]
-	velocity_parts[:,:,i] = sim["/PartType0/Velocities"]
-	time[i] = sim["/Header"].attrs["Time"][0]
+    sim = h5py.File("stellar_evolution_%04d.hdf5" % i, "r")
+    print("reading snapshot " + str(i))
+    abundances[:, :, i] = sim["/PartType0/ElementMassFractions"]
+    metallicity[:, i] = sim["/PartType0/Metallicities"]
+    masses[:, i] = sim["/PartType0/Masses"]
+    star_masses[:, i] = sim["/PartType4/Masses"]
+    mass_from_AGB[:, i] = sim["/PartType0/TotalMassFromAGB"]
+    metal_mass_frac_from_AGB[:, i] = sim["/PartType0/MetalMassFracFromAGB"]
+    mass_from_SNII[:, i] = sim["/PartType0/TotalMassFromSNII"]
+    metal_mass_frac_from_SNII[:, i] = sim["/PartType0/MetalMassFracFromSNII"]
+    mass_from_SNIa[:, i] = sim["/PartType0/TotalMassFromSNIa"]
+    metal_mass_frac_from_SNIa[:, i] = sim["/PartType0/MetalMassFracFromSNIa"]
+    iron_mass_frac_from_SNIa[:, i] = sim["/PartType0/IronMassFracFromSNIa"]
+    internal_energy[:, i] = sim["/PartType0/InternalEnergies"]
+    velocity_parts[:, :, i] = sim["/PartType0/Velocities"]
+    time[i] = sim["/Header"].attrs["Time"][0]
 
 # Define ejecta factor
 ejecta_factor = 1.0e-2
-ejecta_factor_metallicity = 1.0 - 2.0/n_elements
-ejecta_factor_abundances = 1.0/n_elements
+ejecta_factor_metallicity = 1.0 - 2.0 / n_elements
+ejecta_factor_abundances = 1.0 / n_elements
 ejected_mass = star_initial_mass
-energy_per_SNe = 1.0e51/unit_energy_in_cgs
+energy_per_SNe = 1.0e51 / unit_energy_in_cgs
 
 # Check that the total amount of enrichment is as expected.
 # Define tolerance
 eps = 0.01
 
 # Total mass
-total_part_mass = np.sum(masses,axis = 0)
-if abs((total_part_mass[n_snapshots-1] - total_part_mass[0])/total_part_mass[0] - ejected_mass/total_part_mass[0])*total_part_mass[0]/ejected_mass < eps:
-	print("total mass released consistent with expectation")
+total_part_mass = np.sum(masses, axis=0)
+if (
+    abs(
+        (total_part_mass[n_snapshots - 1] - total_part_mass[0]) / total_part_mass[0]
+        - ejected_mass / total_part_mass[0]
+    )
+    * total_part_mass[0]
+    / ejected_mass
+    < eps
+):
+    print("total mass released consistent with expectation")
 else:
-	print("mass increase "+str(total_part_mass[n_snapshots-1]/total_part_mass[0])+" expected "+ str(1.0+ejected_mass/total_part_mass[0]))
+    print(
+        "mass increase "
+        + str(total_part_mass[n_snapshots - 1] / total_part_mass[0])
+        + " expected "
+        + str(1.0 + ejected_mass / total_part_mass[0])
+    )
 
 # Check that mass is conserved (i.e. total star mass decreases by same amount as total gas mass increases)
-total_spart_mass = np.sum(star_masses,axis = 0)
-if abs((total_part_mass[n_snapshots-1] + total_spart_mass[n_snapshots-1]) / (total_part_mass[0] + total_spart_mass[0]) - 1.0) < eps**3:
-	print("total mass conserved")
+total_spart_mass = np.sum(star_masses, axis=0)
+if (
+    abs(
+        (total_part_mass[n_snapshots - 1] + total_spart_mass[n_snapshots - 1])
+        / (total_part_mass[0] + total_spart_mass[0])
+        - 1.0
+    )
+    < eps ** 3
+):
+    print("total mass conserved")
 else:
-	print("initial part, spart mass " + str(total_part_mass[0]) + " " + str(total_spart_mass[0]) + " final mass " + str(total_part_mass[n_snapshots-1]) + " " + str(total_spart_mass[n_snapshots-1]))
+    print(
+        "initial part, spart mass "
+        + str(total_part_mass[0])
+        + " "
+        + str(total_spart_mass[0])
+        + " final mass "
+        + str(total_part_mass[n_snapshots - 1])
+        + " "
+        + str(total_spart_mass[n_snapshots - 1])
+    )
 
 # Total metal mass from AGB
-total_metal_mass_AGB = np.sum(np.multiply(metal_mass_frac_from_AGB,masses),axis = 0)
-expected_metal_mass_AGB = ejecta_factor*ejected_mass
-if abs(total_metal_mass_AGB[n_snapshots-1] - expected_metal_mass_AGB)/expected_metal_mass_AGB < eps:
-	print("total AGB metal mass released consistent with expectation")
+total_metal_mass_AGB = np.sum(np.multiply(metal_mass_frac_from_AGB, masses), axis=0)
+expected_metal_mass_AGB = ejecta_factor * ejected_mass
+if (
+    abs(total_metal_mass_AGB[n_snapshots - 1] - expected_metal_mass_AGB)
+    / expected_metal_mass_AGB
+    < eps
+):
+    print("total AGB metal mass released consistent with expectation")
 else:
-	print("total AGB metal mass "+str(total_metal_mass_AGB[n_snapshots-1])+" expected "+ str(expected_metal_mass_AGB))
+    print(
+        "total AGB metal mass "
+        + str(total_metal_mass_AGB[n_snapshots - 1])
+        + " expected "
+        + str(expected_metal_mass_AGB)
+    )
 
 # Total mass from AGB
-total_AGB_mass = np.sum(mass_from_AGB,axis = 0)
-expected_AGB_mass = ejecta_factor*ejected_mass
-if abs(total_AGB_mass[n_snapshots-1] - expected_AGB_mass)/expected_AGB_mass < eps:
-	print("total AGB mass released consistent with expectation")
+total_AGB_mass = np.sum(mass_from_AGB, axis=0)
+expected_AGB_mass = ejecta_factor * ejected_mass
+if abs(total_AGB_mass[n_snapshots - 1] - expected_AGB_mass) / expected_AGB_mass < eps:
+    print("total AGB mass released consistent with expectation")
 else:
-	print("total AGB mass "+str(total_AGB_mass[n_snapshots-1])+" expected "+ str(expected_AGB_mass))
+    print(
+        "total AGB mass "
+        + str(total_AGB_mass[n_snapshots - 1])
+        + " expected "
+        + str(expected_AGB_mass)
+    )
 
 # Total metal mass from SNII
-total_metal_mass_SNII = np.sum(np.multiply(metal_mass_frac_from_SNII,masses),axis = 0)
-expected_metal_mass_SNII = ejecta_factor*ejected_mass
-if abs(total_metal_mass_SNII[n_snapshots-1] - expected_metal_mass_SNII)/expected_metal_mass_SNII < eps:
-	print("total SNII metal mass released consistent with expectation")
+total_metal_mass_SNII = np.sum(np.multiply(metal_mass_frac_from_SNII, masses), axis=0)
+expected_metal_mass_SNII = ejecta_factor * ejected_mass
+if (
+    abs(total_metal_mass_SNII[n_snapshots - 1] - expected_metal_mass_SNII)
+    / expected_metal_mass_SNII
+    < eps
+):
+    print("total SNII metal mass released consistent with expectation")
 else:
-	print("total SNII metal mass "+str(total_metal_mass_SNII[n_snapshots-1])+" expected "+ str(expected_metal_mass_SNII))
+    print(
+        "total SNII metal mass "
+        + str(total_metal_mass_SNII[n_snapshots - 1])
+        + " expected "
+        + str(expected_metal_mass_SNII)
+    )
 
 # Total mass from SNII
-total_SNII_mass = np.sum(mass_from_SNII,axis = 0)
-expected_SNII_mass = ejecta_factor*ejected_mass
-if abs(total_SNII_mass[n_snapshots-1] - expected_SNII_mass)/expected_SNII_mass < eps:
-	print("total SNII mass released consistent with expectation")
+total_SNII_mass = np.sum(mass_from_SNII, axis=0)
+expected_SNII_mass = ejecta_factor * ejected_mass
+if (
+    abs(total_SNII_mass[n_snapshots - 1] - expected_SNII_mass) / expected_SNII_mass
+    < eps
+):
+    print("total SNII mass released consistent with expectation")
 else:
-	print("total SNII mass "+str(total_SNII_mass[n_snapshots-1])+" expected "+ str(expected_SNII_mass))
+    print(
+        "total SNII mass "
+        + str(total_SNII_mass[n_snapshots - 1])
+        + " expected "
+        + str(expected_SNII_mass)
+    )
 
 # Total metal mass from SNIa
-total_metal_mass_SNIa = np.sum(np.multiply(metal_mass_frac_from_SNIa,masses),axis = 0)
-expected_metal_mass_SNIa = ejecta_factor*ejected_mass
-if abs(total_metal_mass_SNIa[n_snapshots-1] - expected_metal_mass_SNIa)/expected_metal_mass_SNIa < eps:
-	print("total SNIa metal mass released consistent with expectation")
+total_metal_mass_SNIa = np.sum(np.multiply(metal_mass_frac_from_SNIa, masses), axis=0)
+expected_metal_mass_SNIa = ejecta_factor * ejected_mass
+if (
+    abs(total_metal_mass_SNIa[n_snapshots - 1] - expected_metal_mass_SNIa)
+    / expected_metal_mass_SNIa
+    < eps
+):
+    print("total SNIa metal mass released consistent with expectation")
 else:
-	print("total SNIa metal mass "+str(total_metal_mass_SNIa[n_snapshots-1])+" expected "+ str(expected_metal_mass_SNIa))
+    print(
+        "total SNIa metal mass "
+        + str(total_metal_mass_SNIa[n_snapshots - 1])
+        + " expected "
+        + str(expected_metal_mass_SNIa)
+    )
 
 # Total iron mass from SNIa
-total_iron_mass_SNIa = np.sum(np.multiply(iron_mass_frac_from_SNIa,masses),axis = 0)
-expected_iron_mass_SNIa = ejecta_factor*ejected_mass
-if abs(total_iron_mass_SNIa[n_snapshots-1] - expected_iron_mass_SNIa)/expected_iron_mass_SNIa < eps:
-	print("total SNIa iron mass released consistent with expectation")
+total_iron_mass_SNIa = np.sum(np.multiply(iron_mass_frac_from_SNIa, masses), axis=0)
+expected_iron_mass_SNIa = ejecta_factor * ejected_mass
+if (
+    abs(total_iron_mass_SNIa[n_snapshots - 1] - expected_iron_mass_SNIa)
+    / expected_iron_mass_SNIa
+    < eps
+):
+    print("total SNIa iron mass released consistent with expectation")
 else:
-	print("total SNIa iron mass "+str(total_iron_mass_SNIa[n_snapshots-1])+" expected "+ str(expected_iron_mass_SNIa))
+    print(
+        "total SNIa iron mass "
+        + str(total_iron_mass_SNIa[n_snapshots - 1])
+        + " expected "
+        + str(expected_iron_mass_SNIa)
+    )
 
 # Total mass from SNIa
-total_SNIa_mass = np.sum(mass_from_SNIa,axis = 0)
-expected_SNIa_mass = ejecta_factor*ejected_mass
-if abs(total_SNIa_mass[n_snapshots-1] - expected_SNIa_mass)/expected_SNIa_mass < eps:
-	print("total SNIa mass released consistent with expectation")
+total_SNIa_mass = np.sum(mass_from_SNIa, axis=0)
+expected_SNIa_mass = ejecta_factor * ejected_mass
+if (
+    abs(total_SNIa_mass[n_snapshots - 1] - expected_SNIa_mass) / expected_SNIa_mass
+    < eps
+):
+    print("total SNIa mass released consistent with expectation")
 else:
-	print("total SNIa mass "+str(total_SNIa_mass[n_snapshots-1])+" expected "+ str(expected_SNIa_mass))
+    print(
+        "total SNIa mass "
+        + str(total_SNIa_mass[n_snapshots - 1])
+        + " expected "
+        + str(expected_SNIa_mass)
+    )
 
 # Total metal mass
-total_metal_mass = np.sum(np.multiply(metallicity,masses),axis = 0)
-expected_metal_mass = ejecta_factor_metallicity*ejected_mass
-if abs(total_metal_mass[n_snapshots-1] - expected_metal_mass)/expected_metal_mass < eps:
-	print("total metal mass released consistent with expectation")
+total_metal_mass = np.sum(np.multiply(metallicity, masses), axis=0)
+expected_metal_mass = ejecta_factor_metallicity * ejected_mass
+if (
+    abs(total_metal_mass[n_snapshots - 1] - expected_metal_mass) / expected_metal_mass
+    < eps
+):
+    print("total metal mass released consistent with expectation")
 else:
-	print("total metal mass "+str(total_metal_mass[n_snapshots-1])+" expected "+ str(expected_metal_mass))
+    print(
+        "total metal mass "
+        + str(total_metal_mass[n_snapshots - 1])
+        + " expected "
+        + str(expected_metal_mass)
+    )
 
 # Total mass for each element
-expected_element_mass = ejecta_factor_abundances*ejected_mass
+expected_element_mass = ejecta_factor_abundances * ejected_mass
 for i in range(n_elements):
-	total_element_mass = np.sum(np.multiply(abundances[:,i,:],masses),axis = 0)
-	if abs(total_element_mass[n_snapshots-1] - expected_element_mass)/expected_element_mass < eps:
-		print("total element mass released consistent with expectation for element "+str(i))
-	else:
-		print("total element mass "+str(total_element_mass[n_snapshots-1])+" expected "+ str(expected_element_mass) + " for element "+ str(i))
+    total_element_mass = np.sum(np.multiply(abundances[:, i, :], masses), axis=0)
+    if (
+        abs(total_element_mass[n_snapshots - 1] - expected_element_mass)
+        / expected_element_mass
+        < eps
+    ):
+        print(
+            "total element mass released consistent with expectation for element "
+            + str(i)
+        )
+    else:
+        print(
+            "total element mass "
+            + str(total_element_mass[n_snapshots - 1])
+            + " expected "
+            + str(expected_element_mass)
+            + " for element "
+            + str(i)
+        )
+
diff --git a/examples/SubgridTests/StellarEvolution/check_stochastic_heating.py b/examples/SubgridTests/StellarEvolution/check_stochastic_heating.py
index da837540041a9295a33b55e16b5e996394576cd7..1cacc13653d821da3abd2a09566be347608c64f7 100644
--- a/examples/SubgridTests/StellarEvolution/check_stochastic_heating.py
+++ b/examples/SubgridTests/StellarEvolution/check_stochastic_heating.py
@@ -93,7 +93,7 @@ for i in range(n_snapshots):
 	sim = h5py.File("stellar_evolution_%04d.hdf5"%i, "r")
 	print('reading snapshot '+str(i))
 	masses[:,i] = sim["/PartType0/Masses"]
-	internal_energy[:,i] = sim["/PartType0/InternalEnergy"]
+	internal_energy[:,i] = sim["/PartType0/InternalEnergies"]
 	velocity_parts[:,:,i] = sim["/PartType0/Velocities"]
 	time[i] = sim["/Header"].attrs["Time"][0]
 
diff --git a/examples/SubgridTests/StellarEvolution/plot_box_evolution.py b/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
index a46db721e153ade2d386a5790e26a290cead4f90..bcfa85a1afac021e280a797641dd677bccd291d3 100644
--- a/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
+++ b/examples/SubgridTests/StellarEvolution/plot_box_evolution.py
@@ -112,16 +112,16 @@ for i in range(n_snapshots):
 	star_masses = sim["/PartType4/Masses"][:]
 	swift_box_star_mass[i] = np.sum(star_masses)
 
-	metallicities = sim["/PartType0/Metallicity"][:]
+	metallicities = sim["/PartType0/Metallicities"][:]
 	swift_box_gas_metal_mass[i] = np.sum(metallicities * masses)
 
-	element_abundances = sim["/PartType0/ElementAbundance"][:][:]
+	element_abundances = sim["/PartType0/ElementMassFractions"][:][:]
 	for j in range(n_elements):
 		swift_element_mass[i,j] = np.sum(element_abundances[:,j] * masses)
 
         v = sim["/PartType0/Velocities"][:,:]
         v2 = v[:,0]**2 + v[:,1]**2 + v[:,2]**2
-        u = sim["/PartType0/InternalEnergy"][:]
+        u = sim["/PartType0/InternalEnergies"][:]
         swift_internal_energy[i] = np.sum(masses * u)
         swift_kinetic_energy[i] = np.sum(0.5 * masses * v2)
         swift_total_energy[i] = swift_kinetic_energy[i] + swift_internal_energy[i]
diff --git a/examples/SubgridTests/StellarEvolution/plot_particle_evolution.py b/examples/SubgridTests/StellarEvolution/plot_particle_evolution.py
index 8b935e537b14a9d1d9cc4eec7c5cd0794c6fc489..be1588f9b707d448b2905611defd9e760c5f91de 100644
--- a/examples/SubgridTests/StellarEvolution/plot_particle_evolution.py
+++ b/examples/SubgridTests/StellarEvolution/plot_particle_evolution.py
@@ -1,29 +1,30 @@
 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
- #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- # 
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- # 
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- # GNU General Public License for more details.
- # 
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program.  If not, see <http://www.gnu.org/licenses/>.
- # 
- ##############################################################################
-
-# Assuming output snapshots contain evolution of box of gas with star at its 
-# centre, this script will plot the evolution of the radial velocities, internal 
-# energies, mass and metallicities of the nearest n particles to the star over 
-# the duration of the simulation. 
+# This file is part of SWIFT.
+# Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+#                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
+
+# Assuming output snapshots contain evolution of box of gas with star at its
+# centre, this script will plot the evolution of the radial velocities, internal
+# energies, mass and metallicities of the nearest n particles to the star over
+# the duration of the simulation.
 
 import matplotlib
+
 matplotlib.use("Agg")
 from pylab import *
 from scipy import stats
@@ -33,38 +34,42 @@ import glob
 import os.path
 
 # Function to find index in array a for each element in array b
-def find_indices(a,b):
-        result = np.zeros(len(b))
-        for i in range(len(b)):
-                result[i] = ((np.where(a == b[i]))[0])[0]
-        return result
+def find_indices(a, b):
+    result = np.zeros(len(b))
+    for i in range(len(b)):
+        result[i] = ((np.where(a == b[i]))[0])[0]
+    return result
+
 
 # Plot parameters
-params = {'axes.labelsize': 10,
-'axes.titlesize': 10,
-'font.size': 12,
-'legend.fontsize': 12,
-'xtick.labelsize': 10,
-'ytick.labelsize': 10,
-'text.usetex': True,
- 'figure.figsize' : (9.90,6.45),
-'figure.subplot.left'    : 0.1,
-'figure.subplot.right'   : 0.99,
-'figure.subplot.bottom'  : 0.1,
-'figure.subplot.top'     : 0.95,
-'figure.subplot.wspace'  : 0.2,
-'figure.subplot.hspace'  : 0.2,
-'lines.markersize' : 6,
-'lines.linewidth' : 3.,
-'text.latex.unicode': True
+params = {
+    "axes.labelsize": 10,
+    "axes.titlesize": 10,
+    "font.size": 12,
+    "legend.fontsize": 12,
+    "xtick.labelsize": 10,
+    "ytick.labelsize": 10,
+    "text.usetex": True,
+    "figure.figsize": (9.90, 6.45),
+    "figure.subplot.left": 0.1,
+    "figure.subplot.right": 0.99,
+    "figure.subplot.bottom": 0.1,
+    "figure.subplot.top": 0.95,
+    "figure.subplot.wspace": 0.2,
+    "figure.subplot.hspace": 0.2,
+    "lines.markersize": 6,
+    "lines.linewidth": 3.0,
+    "text.latex.unicode": True,
 }
 rcParams.update(params)
-rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+rc("font", **{"family": "sans-serif", "sans-serif": ["Times"]})
 
 
 # Number of snapshots and elements
-newest_snap_name = max(glob.glob('stellar_evolution_*.hdf5'), key=os.path.getctime)
-n_snapshots = int(newest_snap_name.replace('stellar_evolution_','').replace('.hdf5','')) + 1
+newest_snap_name = max(glob.glob("stellar_evolution_*.hdf5"), key=os.path.getctime)
+n_snapshots = (
+    int(newest_snap_name.replace("stellar_evolution_", "").replace(".hdf5", "")) + 1
+)
 n_particles_to_plot = 500
 
 # Read the simulation data
@@ -87,25 +92,25 @@ unit_energy_in_cgs = unit_mass_in_cgs * unit_vel_in_cgs * unit_vel_in_cgs
 unit_length_in_si = 0.01 * unit_length_in_cgs
 unit_mass_in_si = 0.001 * unit_mass_in_cgs
 unit_time_in_si = unit_time_in_cgs
-unit_density_in_cgs = unit_mass_in_cgs*unit_length_in_cgs**-3
-unit_pressure_in_cgs = unit_mass_in_cgs/unit_length_in_cgs*unit_time_in_cgs**-2
-unit_int_energy_in_cgs = unit_energy_in_cgs/unit_mass_in_cgs
-unit_entropy_in_cgs = unit_energy_in_cgs/unit_temp_in_cgs
+unit_density_in_cgs = unit_mass_in_cgs * unit_length_in_cgs ** -3
+unit_pressure_in_cgs = unit_mass_in_cgs / unit_length_in_cgs * unit_time_in_cgs ** -2
+unit_int_energy_in_cgs = unit_energy_in_cgs / unit_mass_in_cgs
+unit_entropy_in_cgs = unit_energy_in_cgs / unit_temp_in_cgs
 Myr_in_cgs = 3.154e13
 Msun_in_cgs = 1.989e33
 
 # Read data of zeroth snapshot
-pos = sim["/PartType0/Coordinates"][:,:]
-x = pos[:,0] - boxSize / 2
-y = pos[:,1] - boxSize / 2
-z = pos[:,2] - boxSize / 2
-vel = sim["/PartType0/Velocities"][:,:]
-r = sqrt(x**2 + y**2 + z**2)
-v_r = (x * vel[:,0] + y * vel[:,1] + z * vel[:,2]) / r
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-rho = sim["/PartType0/Density"][:]
+pos = sim["/PartType0/Coordinates"][:, :]
+x = pos[:, 0] - boxSize / 2
+y = pos[:, 1] - boxSize / 2
+z = pos[:, 2] - boxSize / 2
+vel = sim["/PartType0/Velocities"][:, :]
+r = sqrt(x ** 2 + y ** 2 + z ** 2)
+v_r = (x * vel[:, 0] + y * vel[:, 1] + z * vel[:, 2]) / r
+u = sim["/PartType0/InternalEnergies"][:]
+S = sim["/PartType0/Entropies"][:]
+P = sim["/PartType0/Pressures"][:]
+rho = sim["/PartType0/Densities"][:]
 mass = sim["/PartType0/Masses"][:]
 IDs = sim["/PartType0/ParticleIDs"][:]
 
@@ -123,34 +128,34 @@ t = zeros(n_snapshots)
 
 # Read data from rest of snapshots
 for i in range(n_snapshots):
-	print("reading snapshot "+str(i))
-	# Read the simulation data
-	sim = h5py.File("stellar_evolution_%04d.hdf5"%i, "r")
-	t[i] = sim["/Header"].attrs["Time"][0]
-	
-	pos = sim["/PartType0/Coordinates"][:,:]
-	x = pos[:,0] - boxSize / 2
-	y = pos[:,1] - boxSize / 2
-	z = pos[:,2] - boxSize / 2
-	vel = sim["/PartType0/Velocities"][:,:]
-	r = sqrt(x**2 + y**2 + z**2)
-	v_r = (x * vel[:,0] + y * vel[:,1] + z * vel[:,2]) / r
-	u = sim["/PartType0/InternalEnergy"][:]
-	S = sim["/PartType0/Entropy"][:]
-	P = sim["/PartType0/Pressure"][:]
-	rho = sim["/PartType0/Density"][:]
-	mass = sim["/PartType0/Masses"][:]
-	metallicity = sim["/PartType0/Metallicity"][:]
-	internal_energy = sim["/PartType0/InternalEnergy"][:]
-	IDs = sim["/PartType0/ParticleIDs"][:]
-
-	# Find which particles we want to plot and store their data
-	indices = (find_indices(IDs,part_IDs_to_plot)).astype(int)
-	masses_to_plot[:,i] = mass[indices[:]]
-	v_r_to_plot[:,i] = v_r[indices[:]]
-	metallicities_to_plot[:,i] = metallicity[indices[:]]
-	internal_energies_to_plot[:,i] = internal_energy[indices[:]]
-	
+    print("reading snapshot " + str(i))
+    # Read the simulation data
+    sim = h5py.File("stellar_evolution_%04d.hdf5" % i, "r")
+    t[i] = sim["/Header"].attrs["Time"][0]
+
+    pos = sim["/PartType0/Coordinates"][:, :]
+    x = pos[:, 0] - boxSize / 2
+    y = pos[:, 1] - boxSize / 2
+    z = pos[:, 2] - boxSize / 2
+    vel = sim["/PartType0/Velocities"][:, :]
+    r = sqrt(x ** 2 + y ** 2 + z ** 2)
+    v_r = (x * vel[:, 0] + y * vel[:, 1] + z * vel[:, 2]) / r
+    u = sim["/PartType0/InternalEnergies"][:]
+    S = sim["/PartType0/Entropies"][:]
+    P = sim["/PartType0/Pressures"][:]
+    rho = sim["/PartType0/Densities"][:]
+    mass = sim["/PartType0/Masses"][:]
+    metallicity = sim["/PartType0/Metallicities"][:]
+    internal_energy = sim["/PartType0/InternalEnergies"][:]
+    IDs = sim["/PartType0/ParticleIDs"][:]
+
+    # Find which particles we want to plot and store their data
+    indices = (find_indices(IDs, part_IDs_to_plot)).astype(int)
+    masses_to_plot[:, i] = mass[indices[:]]
+    v_r_to_plot[:, i] = v_r[indices[:]]
+    metallicities_to_plot[:, i] = metallicity[indices[:]]
+    internal_energies_to_plot[:, i] = internal_energy[indices[:]]
+
 
 # Plot the interesting quantities
 figure()
@@ -158,33 +163,61 @@ figure()
 # Radial velocity --------------------------------
 subplot(221)
 for j in range(n_particles_to_plot):
-	plot(t * unit_time_in_cgs / Myr_in_cgs, v_r_to_plot[j,:] * unit_vel_in_cgs, linewidth=0.5, color='k', ms=0.5, alpha=0.1)
+    plot(
+        t * unit_time_in_cgs / Myr_in_cgs,
+        v_r_to_plot[j, :] * unit_vel_in_cgs,
+        linewidth=0.5,
+        color="k",
+        ms=0.5,
+        alpha=0.1,
+    )
 xlabel("Time (Myr)", labelpad=0)
 ylabel("Radial velocity $(\\rm{cm} \cdot \\rm{s}^{-1})$", labelpad=0)
-ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
 
 # Internal energy --------------------------------
 subplot(222)
 for j in range(n_particles_to_plot):
-	plot(t * unit_time_in_cgs / Myr_in_cgs, internal_energies_to_plot[j,:] * unit_energy_in_cgs / unit_mass_in_cgs, linewidth=0.5, color='k', ms=0.5, alpha=0.1)
+    plot(
+        t * unit_time_in_cgs / Myr_in_cgs,
+        internal_energies_to_plot[j, :] * unit_energy_in_cgs / unit_mass_in_cgs,
+        linewidth=0.5,
+        color="k",
+        ms=0.5,
+        alpha=0.1,
+    )
 xlabel("Time (Myr)", labelpad=0)
 ylabel("Internal energy $(\\rm{erg} \cdot \\rm{g}^{-1})$", labelpad=2)
-ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
 
 # Masses --------------------------------
 subplot(223)
 for j in range(n_particles_to_plot):
-	plot(t * unit_time_in_cgs / Myr_in_cgs, masses_to_plot[j,:] * unit_mass_in_cgs / Msun_in_cgs, linewidth=0.5, color='k', ms=0.5, alpha=0.1)
+    plot(
+        t * unit_time_in_cgs / Myr_in_cgs,
+        masses_to_plot[j, :] * unit_mass_in_cgs / Msun_in_cgs,
+        linewidth=0.5,
+        color="k",
+        ms=0.5,
+        alpha=0.1,
+    )
 xlabel("Time (Myr)", labelpad=0)
 ylabel("Mass (Msun)", labelpad=2)
-ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
 
 # Metallicities --------------------------------
 subplot(224)
 for j in range(n_particles_to_plot):
-	plot(t * unit_time_in_cgs / Myr_in_cgs, metallicities_to_plot[j,:] , linewidth=0.5, color='k', ms=0.5, alpha=0.1)
+    plot(
+        t * unit_time_in_cgs / Myr_in_cgs,
+        metallicities_to_plot[j, :],
+        linewidth=0.5,
+        color="k",
+        ms=0.5,
+        alpha=0.1,
+    )
 xlabel("Time (Myr)", labelpad=0)
 ylabel("Metallicity", labelpad=2)
-ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ticklabel_format(style="sci", axis="y", scilimits=(0, 0))
 
 savefig("particle_evolution.png", dpi=200)
diff --git a/examples/main.c b/examples/main.c
index b3577b0e2adec7650bb6acd73b128552071f47d8..491b1ac1e72895ba7f9028d7c0b1bf9181f70360 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -777,6 +777,13 @@ int main(int argc, char *argv[]) {
     else
       bzero(&entropy_floor, sizeof(struct entropy_floor_properties));
 
+    /* Initialise the pressure floor */
+    if (with_hydro)
+      pressure_floor_init(&pressure_floor_props, &prog_const, &us,
+                          &hydro_properties, params);
+    else
+      bzero(&pressure_floor_props, sizeof(struct pressure_floor_properties));
+
     /* Initialise the stars properties */
     if (with_stars)
       stars_props_init(&stars_properties, &prog_const, &us, params,
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 6aeeaad3e41498b48e215f0a761ecdffb486a118..5ef0679d072816b7ee0e042d5e473ec6d02a7599 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -285,6 +285,11 @@ EAGLEEntropyFloor:
   Cool_temperature_norm_K:         8000      # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
   Cool_gamma_effective:            1.        # Slope the of the EAGLE Cool limiter entropy floor
   
+# Parameters related to pressure floors    ----------------------------------------------
+
+GEARPressureFloor:
+  Jeans_factor: 10.       # Number of particles required to suppose a resolved clump and avoid the pressure floor.
+
 # Parameters related to cooling function  ----------------------------------------------
 
 # Constant du/dt cooling function
@@ -408,4 +413,4 @@ EAGLEAGN:
   coupling_efficiency:              0.15       # Fraction of the radiated energy that couples to the gas in feedback events.
   AGN_delta_T_K:                    3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
   AGN_num_ngb_to_heat:              1.         # Target number of gas neighbours to heat in an AGN feedback event.
-
+  max_reposition_mass:              2e8        # Maximal BH mass considered for BH repositioning in solar masses.
diff --git a/src/Makefile.am b/src/Makefile.am
index 0f7cef5e3c6b861ebde639d90624a360f1be7047..3d748b135c93cf409464a199232b243c699fabe1 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -79,7 +79,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c engine_maketasks.c
     collectgroup.c hydro_space.c equation_of_state.c \
     chemistry.c cosmology.c restart.c mesh_gravity.c velociraptor_interface.c \
     outputlist.c velociraptor_dummy.c logger_io.c memuse.c fof.c \
-    hashmap.c \
+    hashmap.c pressure_floor.c \
     $(EAGLE_COOLING_SOURCES) $(EAGLE_FEEDBACK_SOURCES)
 
 # Include files for distribution, not installation.
@@ -98,28 +98,28 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
 	 	 hydro.h hydro_io.h hydro_parameters.h \
 		 hydro/Minimal/hydro.h hydro/Minimal/hydro_iact.h hydro/Minimal/hydro_io.h \
                  hydro/Minimal/hydro_debug.h hydro/Minimal/hydro_part.h \
-				 hydro/Minimal/hydro_parameters.h \
+		 hydro/Minimal/hydro_parameters.h \
 		 hydro/Default/hydro.h hydro/Default/hydro_iact.h hydro/Default/hydro_io.h \
                  hydro/Default/hydro_debug.h hydro/Default/hydro_part.h \
-				 hydro/Default/hydro_parameters.h \
+		 hydro/Default/hydro_parameters.h \
 		 hydro/Gadget2/hydro.h hydro/Gadget2/hydro_iact.h hydro/Gadget2/hydro_io.h \
                  hydro/Gadget2/hydro_debug.h hydro/Gadget2/hydro_part.h \
-				 hydro/Gadget2/hydro_parameters.h \
+		 hydro/Gadget2/hydro_parameters.h \
 		 hydro/PressureEntropy/hydro.h hydro/PressureEntropy/hydro_iact.h hydro/PressureEntropy/hydro_io.h \
                  hydro/PressureEntropy/hydro_debug.h hydro/PressureEntropy/hydro_part.h \
-				 hydro/PressureEntropy/hydro_parameters.h \
+		 hydro/PressureEntropy/hydro_parameters.h \
 		 hydro/PressureEnergy/hydro.h hydro/PressureEnergy/hydro_iact.h hydro/PressureEnergy/hydro_io.h \
                  hydro/PressureEnergy/hydro_debug.h hydro/PressureEnergy/hydro_part.h \
-				 hydro/PressureEnergy/hydro_parameters.h \
+		 hydro/PressureEnergy/hydro_parameters.h \
 		 hydro/PressureEnergyMorrisMonaghanAV/hydro.h hydro/PressureEnergyMorrisMonaghanAV/hydro_iact.h hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h \
                  hydro/PressureEnergyMorrisMonaghanAV/hydro_debug.h hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h \
-				 hydro/PressureEnergyMorrisMonaghanAV/hydro_parameters.h \
+		 hydro/PressureEnergyMorrisMonaghanAV/hydro_parameters.h \
 		 hydro/AnarchyPU/hydro.h hydro/PressureEnergy/hydro_iact.h hydro/PressureEnergy/hydro_io.h \
                  hydro/AnarchyPU/hydro_debug.h hydro/PressureEnergy/hydro_part.h \
-				 hydro/AnarchyPU/hydro_parameters.h \
+		 hydro/AnarchyPU/hydro_parameters.h \
 		 hydro/AnarchyDU/hydro.h hydro/PressureEnergy/hydro_iact.h hydro/PressureEnergy/hydro_io.h \
                  hydro/AnarchyDU/hydro_debug.h hydro/PressureEnergy/hydro_part.h \
-				 hydro/AnarchyDU/hydro_parameters.h \
+		 hydro/AnarchyDU/hydro_parameters.h \
 		 hydro/GizmoMFV/hydro.h hydro/GizmoMFV/hydro_iact.h \
                  hydro/GizmoMFV/hydro_io.h hydro/GizmoMFV/hydro_debug.h \
                  hydro/GizmoMFV/hydro_part.h \
@@ -131,7 +131,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
                  hydro/GizmoMFV/hydro_slope_limiters.h \
                  hydro/GizmoMFV/hydro_unphysical.h \
                  hydro/GizmoMFV/hydro_velocities.h \
-				 hydro/GizmoMFV/hydro_parameters.h \
+		 hydro/GizmoMFV/hydro_parameters.h \
 		 hydro/GizmoMFM/hydro.h hydro/GizmoMFM/hydro_iact.h \
                  hydro/GizmoMFM/hydro_io.h hydro/GizmoMFM/hydro_debug.h \
                  hydro/GizmoMFM/hydro_part.h \
@@ -142,7 +142,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
                  hydro/GizmoMFM/hydro_slope_limiters_face.h \
                  hydro/GizmoMFM/hydro_slope_limiters.h \
                  hydro/GizmoMFM/hydro_unphysical.h \
-				 hydro/GizmoMFM/hydro_parameters.h \
+		 hydro/GizmoMFM/hydro_parameters.h \
                  hydro/Shadowswift/hydro_debug.h \
                  hydro/Shadowswift/hydro_gradients.h hydro/Shadowswift/hydro.h \
                  hydro/Shadowswift/hydro_iact.h \
@@ -190,7 +190,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
                  cooling/const_lambda/cooling_io.h \
                  cooling/grackle/cooling.h cooling/grackle/cooling_struct.h \
                  cooling/grackle/cooling_io.h \
-		 cooling/EAGLE/cooling.h cooling/EAGLE/cooling_struct.h \
+		 cooling/EAGLE/cooling.h cooling/EAGLE/cooling_struct.h cooling/EAGLE/cooling_tables.h \
                  cooling/EAGLE/cooling_io.h cooling/EAGLE/interpolate.h cooling/EAGLE/cooling_rates.h \
                  chemistry/none/chemistry.h \
 		 chemistry/none/chemistry_io.h \
@@ -221,7 +221,10 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
                  black_holes/Default/black_holes_struct.h \
                  black_holes/EAGLE/black_holes.h black_holes/EAGLE/black_holes_io.h \
 		 black_holes/EAGLE/black_holes_part.h black_holes/EAGLE/black_holes_iact.h \
-                 black_holes/EAGLE/black_holes_properties.h 
+                 black_holes/EAGLE/black_holes_properties.h \
+                 black_holes/EAGLE/black_holes_struct.h \
+                 pressure_floor.h \
+		 pressure_floor/GEAR/pressure_floor.h pressure_floor/none/pressure_floor.h
 
 
 # Sources and flags for regular library
diff --git a/src/black_holes/Default/black_holes.h b/src/black_holes/Default/black_holes.h
index 56866764fe4b65603309856bfb2c5ae6b8d82f00..082134cc78f5a7f5e0aa9adeb231c92b9d66c233 100644
--- a/src/black_holes/Default/black_holes.h
+++ b/src/black_holes/Default/black_holes.h
@@ -187,6 +187,20 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
     const struct phys_const* constants, const struct cosmology* cosmo,
     const double dt) {}
 
+/**
+ * @brief Finish the calculation of the new BH position.
+ *
+ * Nothing to do here.
+ *
+ * @param bp The black hole particle.
+ * @param props The properties of the black hole scheme.
+ * @param constants The physical constants (in internal units).
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void black_holes_end_reposition(
+    struct bpart* restrict bp, const struct black_holes_props* props,
+    const struct phys_const* constants, const struct cosmology* cosmo) {}
+
 /**
  * @brief Reset acceleration fields of a particle
  *
diff --git a/src/black_holes/Default/black_holes_iact.h b/src/black_holes/Default/black_holes_iact.h
index ccb1e8c90cd8fef42dd0bb5b6e2d9c6eb9382226..b26b8a50f40599402f4ea9886dbbd7495e8efd59 100644
--- a/src/black_holes/Default/black_holes_iact.h
+++ b/src/black_holes/Default/black_holes_iact.h
@@ -30,16 +30,15 @@
  * @param pj Second particle (gas, not updated).
  * @param xpj The extended data of the second particle (not updated).
  * @param cosmo The cosmological model.
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
  * @param ti_current Current integer time value (for random numbers).
  */
 __attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_bh_gas_density(const float r2, const float *dx,
-                                  const float hi, const float hj,
-                                  struct bpart *restrict bi,
-                                  const struct part *restrict pj,
-                                  const struct xpart *restrict xpj,
-                                  const struct cosmology *cosmo,
-                                  const integertime_t ti_current) {
+runner_iact_nonsym_bh_gas_density(
+    const float r2, const float *dx, const float hi, const float hj,
+    struct bpart *restrict bi, const struct part *restrict pj,
+    const struct xpart *restrict xpj, const struct cosmology *cosmo,
+    const struct gravity_props *grav_props, const integertime_t ti_current) {
 
   float wi, wi_dx;
 
@@ -80,16 +79,15 @@ runner_iact_nonsym_bh_gas_density(const float r2, const float *dx,
  * @param pj Second particle (gas)
  * @param xpj The extended data of the second particle.
  * @param cosmo The cosmological model.
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
  * @param ti_current Current integer time value (for random numbers).
  */
 __attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
-                                  const float hi, const float hj,
-                                  const struct bpart *restrict bi,
-                                  struct part *restrict pj,
-                                  struct xpart *restrict xpj,
-                                  const struct cosmology *cosmo,
-                                  const integertime_t ti_current) {}
+runner_iact_nonsym_bh_gas_swallow(
+    const float r2, const float *dx, const float hi, const float hj,
+    const struct bpart *restrict bi, struct part *restrict pj,
+    struct xpart *restrict xpj, const struct cosmology *cosmo,
+    const struct gravity_props *grav_props, const integertime_t ti_current) {}
 
 /**
  * @brief Swallowing interaction between two BH particles (non-symmetric).
@@ -104,7 +102,7 @@ runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
  * @param bi First particle (black hole).
  * @param bj Second particle (black hole)
  * @param cosmo The cosmological model.
- * @param grav_props The properties of the gravity scheme (softening, G, ...)
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
  * @param ti_current Current integer time value (for random numbers).
  */
 __attribute__((always_inline)) INLINE static void
@@ -127,16 +125,15 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
  * @param pj Second particle (gas)
  * @param xpj The extended data of the second particle.
  * @param cosmo The cosmological model.
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
  * @param ti_current Current integer time value (for random numbers).
  */
 __attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_bh_gas_feedback(const float r2, const float *dx,
-                                   const float hi, const float hj,
-                                   struct bpart *restrict bi,
-                                   struct part *restrict pj,
-                                   struct xpart *restrict xpj,
-                                   const struct cosmology *cosmo,
-                                   const integertime_t ti_current) {
+runner_iact_nonsym_bh_gas_feedback(
+    const float r2, const float *dx, const float hi, const float hj,
+    struct bpart *restrict bi, struct part *restrict pj,
+    struct xpart *restrict xpj, const struct cosmology *cosmo,
+    const struct gravity_props *grav_props, const integertime_t ti_current) {
 #ifdef DEBUG_INTERACTIONS_BH
   /* Update ngb counters */
   if (si->num_ngb_force < MAX_NUM_OF_NEIGHBOURS_BH)
diff --git a/src/black_holes/Default/black_holes_io.h b/src/black_holes/Default/black_holes_io.h
index 5288ae595663871bb3187782712d52708fac7931..53322da1dd90a98fb0f1f2529774a20dc718b30c 100644
--- a/src/black_holes/Default/black_holes_io.h
+++ b/src/black_holes/Default/black_holes_io.h
@@ -102,6 +102,7 @@ INLINE static void convert_bpart_vel(const struct engine* e,
  * @param bparts The b-particle array.
  * @param list The list of i/o properties to write.
  * @param num_fields The number of i/o fields to write.
+ * @param with_cosmology Are we running a cosmological simulation?
  */
 INLINE static void black_holes_write_particles(const struct bpart* bparts,
                                                struct io_props* list,
diff --git a/src/black_holes/EAGLE/black_holes.h b/src/black_holes/EAGLE/black_holes.h
index 3bf4e2295e1e5f13993e94f76ef6b83882a5e2e7..a1b41f0954e1bec890329779d0edb7fcd31265f1 100644
--- a/src/black_holes/EAGLE/black_holes.h
+++ b/src/black_holes/EAGLE/black_holes.h
@@ -24,6 +24,7 @@
 #include "black_holes_struct.h"
 #include "cosmology.h"
 #include "dimension.h"
+#include "gravity.h"
 #include "kernel_hydro.h"
 #include "minmax.h"
 #include "physical_constants.h"
@@ -89,16 +90,43 @@ __attribute__((always_inline)) INLINE static void black_holes_init_bpart(
   bp->circular_velocity_gas[2] = 0.f;
   bp->ngb_mass = 0.f;
   bp->num_ngbs = 0;
+  bp->reposition.x[0] = -FLT_MAX;
+  bp->reposition.x[1] = -FLT_MAX;
+  bp->reposition.x[2] = -FLT_MAX;
+  bp->reposition.min_potential = FLT_MAX;
 }
 
 /**
  * @brief Predict additional particle fields forward in time when drifting
  *
+ * The fields do not get predicted but we move the BH to its new position
+ * if a new one was calculated in the repositioning loop.
+ *
  * @param bp The particle
  * @param dt_drift The drift time-step for positions.
  */
 __attribute__((always_inline)) INLINE static void black_holes_predict_extra(
-    struct bpart* restrict bp, float dt_drift) {}
+    struct bpart* restrict bp, float dt_drift) {
+
+  /* Are we doing some repositioning? */
+  if (bp->reposition.min_potential != FLT_MAX) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+    if (bp->reposition.x[0] == -FLT_MAX || bp->reposition.x[1] == -FLT_MAX ||
+        bp->reposition.x[2] == -FLT_MAX) {
+      error("Something went wrong with the new repositioning position");
+    }
+#endif
+
+    bp->x[0] = bp->reposition.x[0];
+    bp->x[1] = bp->reposition.x[1];
+    bp->x[2] = bp->reposition.x[2];
+
+    bp->gpart->x[0] = bp->reposition.x[0];
+    bp->gpart->x[1] = bp->reposition.x[1];
+    bp->gpart->x[2] = bp->reposition.x[2];
+  }
+}
 
 /**
  * @brief Sets the values to be predicted in the drifts to their values at a
@@ -430,6 +458,35 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
   }
 }
 
+/**
+ * @brief Finish the calculation of the new BH position.
+ *
+ * Here, we check that the BH should indeed be moved in the next drift.
+ *
+ * @param bp The black hole particle.
+ * @param props The properties of the black hole scheme.
+ * @param constants The physical constants (in internal units).
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void black_holes_end_reposition(
+    struct bpart* restrict bp, const struct black_holes_props* props,
+    const struct phys_const* constants, const struct cosmology* cosmo) {
+
+  const float potential = gravity_get_comoving_potential(bp->gpart);
+
+  /* Is the potential lower (i.e. the BH is at the bottom already)
+   * OR is the BH massive enough that we don't reposition? */
+  if (potential < bp->reposition.min_potential ||
+      bp->subgrid_mass > props->max_reposition_mass) {
+
+    /* No need to reposition */
+    bp->reposition.min_potential = FLT_MAX;
+    bp->reposition.x[0] = -FLT_MAX;
+    bp->reposition.x[1] = -FLT_MAX;
+    bp->reposition.x[2] = -FLT_MAX;
+  }
+}
+
 /**
  * @brief Reset acceleration fields of a particle
  *
diff --git a/src/black_holes/EAGLE/black_holes_iact.h b/src/black_holes/EAGLE/black_holes_iact.h
index e747060a76e0a5bb027a2b1b553edb45aa2542bc..7e9852dee582ab49643368c0a7b0d79ecd6b5047 100644
--- a/src/black_holes/EAGLE/black_holes_iact.h
+++ b/src/black_holes/EAGLE/black_holes_iact.h
@@ -20,6 +20,7 @@
 #define SWIFT_EAGLE_BH_IACT_H
 
 /* Local includes */
+#include "gravity.h"
 #include "hydro.h"
 #include "random.h"
 #include "space.h"
@@ -35,16 +36,15 @@
  * @param pj Second particle (gas, not updated).
  * @param xpj The extended data of the second particle (not updated).
  * @param cosmo The cosmological model.
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
  * @param ti_current Current integer time value (for random numbers).
  */
 __attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_bh_gas_density(const float r2, const float *dx,
-                                  const float hi, const float hj,
-                                  struct bpart *restrict bi,
-                                  const struct part *restrict pj,
-                                  const struct xpart *restrict xpj,
-                                  const struct cosmology *cosmo,
-                                  const integertime_t ti_current) {
+runner_iact_nonsym_bh_gas_density(
+    const float r2, const float *dx, const float hi, const float hj,
+    struct bpart *restrict bi, const struct part *restrict pj,
+    const struct xpart *restrict xpj, const struct cosmology *cosmo,
+    const struct gravity_props *grav_props, const integertime_t ti_current) {
 
   float wi, wi_dx;
 
@@ -117,16 +117,15 @@ runner_iact_nonsym_bh_gas_density(const float r2, const float *dx,
  * @param pj Second particle (gas)
  * @param xpj The extended data of the second particle.
  * @param cosmo The cosmological model.
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
  * @param ti_current Current integer time value (for random numbers).
  */
 __attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
-                                  const float hi, const float hj,
-                                  const struct bpart *restrict bi,
-                                  struct part *restrict pj,
-                                  struct xpart *restrict xpj,
-                                  const struct cosmology *cosmo,
-                                  const integertime_t ti_current) {
+runner_iact_nonsym_bh_gas_swallow(
+    const float r2, const float *dx, const float hi, const float hj,
+    struct bpart *restrict bi, struct part *restrict pj,
+    struct xpart *restrict xpj, const struct cosmology *cosmo,
+    const struct gravity_props *grav_props, const integertime_t ti_current) {
 
   float wi;
 
@@ -140,6 +139,44 @@ runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
   const float ui = r * hi_inv;
   kernel_eval(ui, &wi);
 
+  /* Start by checking the repositioning criteria */
+
+  /* Note the factor 9 is taken from EAGLE. Will be turned into a parameter */
+  const float max_dist_repos2 =
+      kernel_gravity_softening_plummer_equivalent_inv *
+      kernel_gravity_softening_plummer_equivalent_inv * 9.f *
+      grav_props->epsilon_cur2;
+
+  /* This gas neighbour is close enough that we can consider it's potential
+     for repositioning */
+  if (r2 < max_dist_repos2) {
+
+    /* Compute relative peculiar velocity between the two BHs
+     * Recall that in SWIFT v is (v_pec * a) */
+    const float delta_v[3] = {bi->v[0] - pj->v[0], bi->v[1] - pj->v[1],
+                              bi->v[2] - pj->v[2]};
+    const float v2 = delta_v[0] * delta_v[0] + delta_v[1] * delta_v[1] +
+                     delta_v[2] * delta_v[2];
+
+    const float v2_pec = v2 * cosmo->a2_inv;
+
+    /* Check the velocity criterion */
+    if (v2_pec < 0.25f * bi->sound_speed_gas * bi->sound_speed_gas) {
+
+      const float potential = gravity_get_comoving_potential(pj->gpart);
+
+      /* Is the potential lower? */
+      if (potential < bi->reposition.min_potential) {
+
+        /* Store this as our new best */
+        bi->reposition.min_potential = potential;
+        bi->reposition.x[0] = pj->x[0];
+        bi->reposition.x[1] = pj->x[1];
+        bi->reposition.x[2] = pj->x[2];
+      }
+    }
+  }
+
   /* Is the BH hungry? */
   if (bi->subgrid_mass > bi->mass) {
 
@@ -188,13 +225,13 @@ runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
  * @param bi First particle (black hole).
  * @param bj Second particle (black hole)
  * @param cosmo The cosmological model.
- * @param grav_props The properties of the gravity scheme (softening, G, ...)
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
  * @param ti_current Current integer time value (for random numbers).
  */
 __attribute__((always_inline)) INLINE static void
 runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
                                  const float hi, const float hj,
-                                 const struct bpart *restrict bi,
+                                 struct bpart *restrict bi,
                                  struct bpart *restrict bj,
                                  const struct cosmology *cosmo,
                                  const struct gravity_props *grav_props,
@@ -209,6 +246,33 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
 
   const float v2_pec = v2 * cosmo->a2_inv;
 
+  /* Note the factor 9 is taken from EAGLE. Will be turned into a parameter */
+  const float max_dist_repos2 =
+      kernel_gravity_softening_plummer_equivalent_inv *
+      kernel_gravity_softening_plummer_equivalent_inv * 9.f *
+      grav_props->epsilon_cur2;
+
+  /* This gas neighbour is close enough that we can consider it's potential
+     for repositioning */
+  if (r2 < max_dist_repos2) {
+
+    /* Check the velocity criterion */
+    if (v2_pec < 0.25f * bi->sound_speed_gas * bi->sound_speed_gas) {
+
+      const float potential = gravity_get_comoving_potential(bj->gpart);
+
+      /* Is the potential lower? */
+      if (potential < bi->reposition.min_potential) {
+
+        /* Store this as our new best */
+        bi->reposition.min_potential = potential;
+        bi->reposition.x[0] = bj->x[0];
+        bi->reposition.x[1] = bj->x[1];
+        bi->reposition.x[2] = bj->x[2];
+      }
+    }
+  }
+
   /* Find the most massive of the two BHs */
   float M = bi->subgrid_mass;
   float h = hi;
@@ -218,9 +282,10 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
   }
 
   /* Note the factor 9 is taken from EAGLE. Will be turned into a parameter */
-  const float max_distance2 = kernel_gravity_softening_plummer_equivalent_inv *
-                              kernel_gravity_softening_plummer_equivalent_inv *
-                              9.f * grav_props->epsilon_cur2;
+  const float max_dist_merge2 =
+      kernel_gravity_softening_plummer_equivalent_inv *
+      kernel_gravity_softening_plummer_equivalent_inv * 9.f *
+      grav_props->epsilon_cur2;
 
   const float G_Newton = grav_props->G_Newton;
 
@@ -234,7 +299,7 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
     /* Merge if gravitationally bound AND if within max distance
      * Note that we use the kernel support here as the size and not just the
      * smoothing length */
-    if (v2_pec < G_Newton * M / (kernel_gamma * h) && (r2 < max_distance2)) {
+    if (v2_pec < G_Newton * M / (kernel_gamma * h) && (r2 < max_dist_merge2)) {
 
       /* This particle is swallowed by the BH with the largest ID of all the
        * candidates wanting to swallow it */
@@ -269,16 +334,15 @@ runner_iact_nonsym_bh_bh_swallow(const float r2, const float *dx,
  * @param pj Second particle (gas)
  * @param xpj The extended data of the second particle.
  * @param cosmo The cosmological model.
+ * @param grav_props The properties of the gravity scheme (softening, G, ...).
  * @param ti_current Current integer time value (for random numbers).
  */
 __attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_bh_gas_feedback(const float r2, const float *dx,
-                                   const float hi, const float hj,
-                                   const struct bpart *restrict bi,
-                                   struct part *restrict pj,
-                                   struct xpart *restrict xpj,
-                                   const struct cosmology *cosmo,
-                                   const integertime_t ti_current) {
+runner_iact_nonsym_bh_gas_feedback(
+    const float r2, const float *dx, const float hi, const float hj,
+    const struct bpart *restrict bi, struct part *restrict pj,
+    struct xpart *restrict xpj, const struct cosmology *cosmo,
+    const struct gravity_props *grav_props, const integertime_t ti_current) {
 
   /* Get the heating probability */
   const float prob = bi->to_distribute.AGN_heating_probability;
diff --git a/src/black_holes/EAGLE/black_holes_io.h b/src/black_holes/EAGLE/black_holes_io.h
index ea8b6f70d2847e65128cb23051ee059a1460b52d..cd6c4599f93b3a324ed2be16e6e3fb4ee06ec02c 100644
--- a/src/black_holes/EAGLE/black_holes_io.h
+++ b/src/black_holes/EAGLE/black_holes_io.h
@@ -104,6 +104,7 @@ INLINE static void convert_bpart_vel(const struct engine* e,
  * @param bparts The b-particle array.
  * @param list The list of i/o properties to write.
  * @param num_fields The number of i/o fields to write.
+ * @param with_cosmology Are we running a cosmological simulation?
  */
 INLINE static void black_holes_write_particles(const struct bpart* bparts,
                                                struct io_props* list,
diff --git a/src/black_holes/EAGLE/black_holes_part.h b/src/black_holes/EAGLE/black_holes_part.h
index 0948ff9b009c44c1a5d2a9b37673e4d96db08f50..cbfa0fcaf560d1c771d8595ec84ad08e9ba3608e 100644
--- a/src/black_holes/EAGLE/black_holes_part.h
+++ b/src/black_holes/EAGLE/black_holes_part.h
@@ -121,6 +121,16 @@ struct bpart {
 
   } to_distribute;
 
+  struct {
+
+    /*! Value of the minimum potential across all neighbours. */
+    float min_potential;
+
+    /*! New position of the BH following the reposition procedure */
+    double x[3];
+
+  } reposition;
+
   /*! Chemistry information (e.g. metal content at birth, swallowed metal
    * content, etc.) */
   struct chemistry_bpart_data chemistry_data;
diff --git a/src/black_holes/EAGLE/black_holes_properties.h b/src/black_holes/EAGLE/black_holes_properties.h
index 50981b43e1ec23927d9958c95ec396fe4a4e4d3a..db50e5b3e08b1a66497aaf75da9d3acdea3287ee 100644
--- a/src/black_holes/EAGLE/black_holes_properties.h
+++ b/src/black_holes/EAGLE/black_holes_properties.h
@@ -74,6 +74,11 @@ struct black_holes_props {
   /*! Number of gas neighbours to heat in a feedback event */
   float num_ngbs_to_heat;
 
+  /* ---- Properties of the repositioning model --- */
+
+  /*! Maximal mass of BH to reposition */
+  float max_reposition_mass;
+
   /* ---- Common conversion factors --------------- */
 
   /*! Conversion factor from temperature to internal energy */
@@ -157,6 +162,14 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
   bp->num_ngbs_to_heat =
       parser_get_param_float(params, "EAGLEAGN:AGN_num_ngb_to_heat");
 
+  /* Reposition parameters --------------------------------- */
+
+  bp->max_reposition_mass =
+      parser_get_param_float(params, "EAGLEAGN:max_reposition_mass");
+
+  /* Convert to internal units */
+  bp->max_reposition_mass *= phys_const->const_solar_mass;
+
   /* Common conversion factors ----------------------------- */
 
   /* Calculate temperature to internal energy conversion factor (all internal
diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h
index a04d7f94e7340a68d25914828de40bb2e1e8020c..c06b29f4417673bbdac79ab76226770c9c327406 100644
--- a/src/chemistry/EAGLE/chemistry.h
+++ b/src/chemistry/EAGLE/chemistry.h
@@ -170,9 +170,10 @@ __attribute__((always_inline)) INLINE static void chemistry_first_init_part(
     p->chemistry_data.metal_mass_fraction_total =
         data->initial_metal_mass_fraction_total;
 
-    for (int elem = 0; elem < chemistry_element_count; ++elem)
+    for (int elem = 0; elem < chemistry_element_count; ++elem) {
       p->chemistry_data.metal_mass_fraction[elem] =
           data->initial_metal_mass_fraction[elem];
+    }
   }
   chemistry_init_part(p, data);
 }
@@ -241,4 +242,106 @@ static INLINE void chemistry_print_backend(
           chemistry_element_count);
 }
 
+/**
+ * @brief Updates the metal mass fractions after diffusion at the end of the
+ * force loop.
+ *
+ * @param p The particle to act upon.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_end_force(
+    struct part* restrict p, const struct cosmology* cosmo) {}
+
+/**
+ * @brief Computes the chemistry-related time-step constraint.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param cosmo The current cosmological model.
+ * @param us The internal system of units.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float chemistry_timestep(
+    const struct phys_const* restrict phys_const,
+    const struct cosmology* restrict cosmo,
+    const struct unit_system* restrict us,
+    const struct hydro_props* hydro_props,
+    const struct chemistry_global_data* cd, const struct part* restrict p) {
+  return FLT_MAX;
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * star particle to be used in feedback/enrichment related routines.
+ *
+ * EAGLE uses smooth abundances for everything.
+ *
+ * @param sp Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_total_metal_mass_fraction_for_feedback(
+    const struct spart* restrict sp) {
+
+  return sp->chemistry_data.smoothed_metal_mass_fraction_total;
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * gas particle to be used in cooling related routines.
+ *
+ * EAGLE uses smooth abundances for everything.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_total_metal_mass_fraction_for_cooling(
+    const struct part* restrict p) {
+
+  return p->chemistry_data.smoothed_metal_mass_fraction_total;
+}
+
+/**
+ * @brief Returns the abundance array (metal mass fractions) of the
+ * gas particle to be used in cooling related routines.
+ *
+ * EAGLE uses smooth abundances for everything.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float const*
+chemistry_get_metal_mass_fraction_for_cooling(const struct part* restrict p) {
+
+  return p->chemistry_data.smoothed_metal_mass_fraction;
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * gas particle to be used in star formation related routines.
+ *
+ * EAGLE uses smooth abundances for everything.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_total_metal_mass_fraction_for_star_formation(
+    const struct part* restrict p) {
+
+  return p->chemistry_data.smoothed_metal_mass_fraction_total;
+}
+
+/**
+ * @brief Returns the abundance array (metal mass fractions) of the
+ * gas particle to be used in star formation related routines.
+ *
+ * EAGLE uses smooth abundances for everything.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float const*
+chemistry_get_metal_mass_fraction_for_star_formation(
+    const struct part* restrict p) {
+
+  return p->chemistry_data.smoothed_metal_mass_fraction;
+}
+
 #endif /* SWIFT_CHEMISTRY_EAGLE_H */
diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c
index 3fda0c80040b6ad139b9100064b937e695e6030f..c5595ef514d5a145dd6e9c326c1ff1f2bc9da196 100644
--- a/src/cooling/EAGLE/cooling.c
+++ b/src/cooling/EAGLE/cooling.c
@@ -207,7 +207,7 @@ INLINE static double bisection_iter(
     const float d_He, const double Lambda_He_reion_cgs,
     const double ratefact_cgs,
     const struct cooling_function_data *restrict cooling,
-    const float abundance_ratio[chemistry_element_count + 2],
+    const float abundance_ratio[eagle_cooling_N_abundances],
     const double dt_cgs, const long long ID) {
 
   /* Bracketing */
@@ -419,14 +419,14 @@ void cooling_cool_part(const struct phys_const *phys_const,
    * Note that we need to add S and Ca that are in the tables but not tracked
    * by the particles themselves.
    * The order is [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe] */
-  float abundance_ratio[chemistry_element_count + 2];
+  float abundance_ratio[eagle_cooling_N_abundances];
   abundance_ratio_to_solar(p, cooling, abundance_ratio);
 
   /* Get the Hydrogen and Helium mass fractions */
-  const float XH =
-      p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_H];
-  const float XHe =
-      p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_He];
+  const float *const metal_fraction =
+      chemistry_get_metal_mass_fraction_for_cooling(p);
+  const float XH = metal_fraction[chemistry_element_H];
+  const float XHe = metal_fraction[chemistry_element_He];
 
   /* Get the Helium mass fraction. Note that this is He / (H + He), i.e. a
    * metal-free Helium mass fraction as per the Wiersma+08 definition */
@@ -600,10 +600,10 @@ float cooling_get_temperature(
   const double u_cgs = u * cooling->internal_energy_to_cgs;
 
   /* Get the Hydrogen and Helium mass fractions */
-  const float XH =
-      p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_H];
-  const float XHe =
-      p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_He];
+  const float *const metal_fraction =
+      chemistry_get_metal_mass_fraction_for_cooling(p);
+  const float XH = metal_fraction[chemistry_element_H];
+  const float XHe = metal_fraction[chemistry_element_He];
 
   /* Get the Helium mass fraction. Note that this is He / (H + He), i.e. a
    * metal-free Helium mass fraction as per the Wiersma+08 definition */
diff --git a/src/cooling/EAGLE/cooling_rates.h b/src/cooling/EAGLE/cooling_rates.h
index 6707a9a9d70b58eb218cfc8ebeca9b1f08fdcb12..39978c9bfa8cf2fbea25e8856c2939a8d6ee69e4 100644
--- a/src/cooling/EAGLE/cooling_rates.h
+++ b/src/cooling/EAGLE/cooling_rates.h
@@ -50,57 +50,50 @@
  * @param cooling #cooling_function_data struct.
  * @param ratio_solar (return) Array of ratios to solar abundances.
  */
-__attribute__((always_inline)) INLINE void abundance_ratio_to_solar(
+__attribute__((always_inline)) INLINE static void abundance_ratio_to_solar(
     const struct part *p, const struct cooling_function_data *cooling,
-    float ratio_solar[chemistry_element_count + 2]) {
+    float ratio_solar[eagle_cooling_N_abundances]) {
 
-  ratio_solar[0] =
-      p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_H] *
-      cooling->SolarAbundances_inv[0 /* H */];
+  /* Get the individual metal mass fractions from the particle */
+  const float *const metal_fraction =
+      chemistry_get_metal_mass_fraction_for_cooling(p);
 
-  ratio_solar[1] =
-      p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_He] *
-      cooling->SolarAbundances_inv[1 /* He */];
+  ratio_solar[0] = metal_fraction[chemistry_element_H] *
+                   cooling->SolarAbundances_inv[0 /* H */];
 
-  ratio_solar[2] =
-      p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_C] *
-      cooling->SolarAbundances_inv[2 /* C */];
+  ratio_solar[1] = metal_fraction[chemistry_element_He] *
+                   cooling->SolarAbundances_inv[1 /* He */];
 
-  ratio_solar[3] =
-      p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_N] *
-      cooling->SolarAbundances_inv[3 /* N */];
+  ratio_solar[2] = metal_fraction[chemistry_element_C] *
+                   cooling->SolarAbundances_inv[2 /* C */];
 
-  ratio_solar[4] =
-      p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_O] *
-      cooling->SolarAbundances_inv[4 /* O */];
+  ratio_solar[3] = metal_fraction[chemistry_element_N] *
+                   cooling->SolarAbundances_inv[3 /* N */];
 
-  ratio_solar[5] =
-      p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_Ne] *
-      cooling->SolarAbundances_inv[5 /* Ne */];
+  ratio_solar[4] = metal_fraction[chemistry_element_O] *
+                   cooling->SolarAbundances_inv[4 /* O */];
 
-  ratio_solar[6] =
-      p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_Mg] *
-      cooling->SolarAbundances_inv[6 /* Mg */];
+  ratio_solar[5] = metal_fraction[chemistry_element_Ne] *
+                   cooling->SolarAbundances_inv[5 /* Ne */];
 
-  ratio_solar[7] =
-      p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_Si] *
-      cooling->SolarAbundances_inv[7 /* Si */];
+  ratio_solar[6] = metal_fraction[chemistry_element_Mg] *
+                   cooling->SolarAbundances_inv[6 /* Mg */];
+
+  ratio_solar[7] = metal_fraction[chemistry_element_Si] *
+                   cooling->SolarAbundances_inv[7 /* Si */];
 
   /* For S, we use the same ratio as Si */
-  ratio_solar[8] =
-      p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_Si] *
-      cooling->SolarAbundances_inv[7 /* Si */] *
-      cooling->S_over_Si_ratio_in_solar;
+  ratio_solar[8] = metal_fraction[chemistry_element_Si] *
+                   cooling->SolarAbundances_inv[7 /* Si */] *
+                   cooling->S_over_Si_ratio_in_solar;
 
   /* For Ca, we use the same ratio as Si */
-  ratio_solar[9] =
-      p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_Si] *
-      cooling->SolarAbundances_inv[7 /* Si */] *
-      cooling->Ca_over_Si_ratio_in_solar;
-
-  ratio_solar[10] =
-      p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_Fe] *
-      cooling->SolarAbundances_inv[10 /* Fe */];
+  ratio_solar[9] = metal_fraction[chemistry_element_Si] *
+                   cooling->SolarAbundances_inv[7 /* Si */] *
+                   cooling->Ca_over_Si_ratio_in_solar;
+
+  ratio_solar[10] = metal_fraction[chemistry_element_Fe] *
+                    cooling->SolarAbundances_inv[10 /* Fe */];
 }
 
 /**
@@ -313,7 +306,7 @@ __attribute__((always_inline)) INLINE double eagle_Compton_cooling_rate(
  */
 INLINE static double eagle_metal_cooling_rate(
     const double log10_u_cgs, const double redshift, const double n_H_cgs,
-    const float solar_ratio[chemistry_element_count + 2], const int n_H_index,
+    const float solar_ratio[eagle_cooling_N_abundances], const int n_H_index,
     const float d_n_H, const int He_index, const float d_He,
     const struct cooling_function_data *cooling, double *element_lambda) {
 
@@ -537,7 +530,7 @@ INLINE static double eagle_metal_cooling_rate(
  */
 INLINE static double eagle_cooling_rate(
     const double log10_u_cgs, const double redshift, const double n_H_cgs,
-    const float abundance_ratio[chemistry_element_count + 2],
+    const float abundance_ratio[eagle_cooling_N_abundances],
     const int n_H_index, const float d_n_H, const int He_index,
     const float d_He, const struct cooling_function_data *cooling) {
 
diff --git a/src/cooling/EAGLE/cooling_tables.c b/src/cooling/EAGLE/cooling_tables.c
index cddc8d50e0c2f00c30846ebaf65b2bf250e9cc8a..4261e9ac0a6fee9f77c03afe22b7a9b66ade487d 100644
--- a/src/cooling/EAGLE/cooling_tables.c
+++ b/src/cooling/EAGLE/cooling_tables.c
@@ -213,10 +213,6 @@ void read_cooling_header(const char *fname,
   if (N_SolarAbundances != eagle_cooling_N_abundances)
     error("Invalid solar abundances array length.");
 
-  /* Check value */
-  if (N_SolarAbundances != chemistry_element_count + 2)
-    error("Number of abundances not compatible with the chemistry model.");
-
   dataset = H5Dopen(tempfile_id, "/Header/Number_of_metals", H5P_DEFAULT);
   status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
                    &N_Elements);
diff --git a/src/cooling/EAGLE/newton_cooling.c b/src/cooling/EAGLE/newton_cooling.c
index c68a77614425f85da011ab45d2ec488710e068dc..50af3294c31dfc2d44a679c26a78629d53d03420 100644
--- a/src/cooling/EAGLE/newton_cooling.c
+++ b/src/cooling/EAGLE/newton_cooling.c
@@ -214,7 +214,7 @@ __attribute__((always_inline)) INLINE double eagle_convert_u_to_temp(
  */
 INLINE static double eagle_metal_cooling_rate(
     double log10_u_cgs, double redshift, double n_H_cgs,
-    const float solar_ratio[chemistry_element_count + 2], int n_H_index,
+    const float solar_ratio[eagle_cooling_N_abundances], int n_H_index,
     float d_n_H, int He_index, float d_He,
     const struct cooling_function_data *restrict cooling, double *dlambda_du,
     double *element_lambda) {
@@ -569,7 +569,7 @@ INLINE static double eagle_metal_cooling_rate(
    */
   INLINE static double eagle_cooling_rate(
       double log_u_cgs, double redshift, double n_H_cgs,
-      const float abundance_ratio[chemistry_element_count + 2], int n_H_index,
+      const float abundance_ratio[eagle_cooling_N_abundances], int n_H_index,
       float d_n_H, int He_index, float d_He,
       const struct cooling_function_data *restrict cooling,
       double *dLambdaNet_du) {
@@ -609,7 +609,7 @@ INLINE static double eagle_metal_cooling_rate(
                     const struct cosmology *restrict cosmo,
                     const struct cooling_function_data *restrict cooling,
                     const struct phys_const *restrict phys_const,
-                    const float abundance_ratio[chemistry_element_count + 2],
+                    const float abundance_ratio[eagle_cooling_N_abundances],
                     float dt, int *bisection_flag) {
 
     double logu, logu_old;
@@ -621,8 +621,9 @@ INLINE static double eagle_metal_cooling_rate(
     const float log_table_bound_low = (cooling->Therm[0] + 0.05) / M_LOG10E;
 
     /* convert Hydrogen mass fraction in Hydrogen number density */
-    const float XH =
-        p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_H];
+    float const *metal_fraction =
+        chemistry_get_metal_mass_fraction_for_cooling(p);
+    const float XH = metal_fraction[chemistry_element_H];
     const double n_H = hydro_get_physical_density(p, cosmo) * XH /
                        phys_const->const_proton_mass;
     const double n_H_cgs = n_H * cooling->number_density_to_cgs;
diff --git a/src/cooling/grackle/cooling.h b/src/cooling/grackle/cooling.h
index 57f4b146c451557365e08d260ddd45276844af9c..98d34dc630aa0c83198ddd60b6ee9e2775778d0f 100644
--- a/src/cooling/grackle/cooling.h
+++ b/src/cooling/grackle/cooling.h
@@ -52,7 +52,7 @@ static gr_float cooling_time(
     const struct unit_system* restrict us,
     const struct cosmology* restrict cosmo,
     const struct cooling_function_data* restrict cooling,
-    const struct part* restrict p, const struct xpart* restrict xp);
+    const struct part* restrict p, struct xpart* restrict xp);
 static gr_float cooling_new_energy(
     const struct phys_const* restrict phys_const,
     const struct unit_system* restrict us,
@@ -573,7 +573,7 @@ __attribute__((always_inline)) INLINE static gr_float cooling_time(
     const struct unit_system* restrict us,
     const struct cosmology* restrict cosmo,
     const struct cooling_function_data* restrict cooling,
-    const struct part* restrict p, const struct xpart* restrict xp) {
+    const struct part* restrict p, struct xpart* restrict xp) {
 
   /* set current time */
   code_units units = cooling->units;
diff --git a/src/engine.c b/src/engine.c
index a367fb43567c31023699d848f47f416fc6d57bca..1f701d7b7e5699a7dd377a75ec914df54db9747e 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -133,6 +133,9 @@ int engine_rank;
 /** The current step of the engine as a global variable (for messages). */
 int engine_current_step;
 
+extern int engine_max_parts_per_ghost;
+extern int engine_max_sparts_per_ghost;
+
 /**
  * @brief Data collected from the cells at the end of a time-step
  */
@@ -5696,6 +5699,51 @@ void engine_config(int restart, int fof, struct engine *e,
   e->sched.mpi_message_limit =
       parser_get_opt_param_int(params, "Scheduler:mpi_message_limit", 4) * 1024;
 
+  if (restart) {
+
+    /* Overwrite the constants for the scheduler */
+    space_maxsize = parser_get_opt_param_int(params, "Scheduler:cell_max_size",
+                                             space_maxsize_default);
+    space_subsize_pair_hydro =
+        parser_get_opt_param_int(params, "Scheduler:cell_sub_size_pair_hydro",
+                                 space_subsize_pair_hydro_default);
+    space_subsize_self_hydro =
+        parser_get_opt_param_int(params, "Scheduler:cell_sub_size_self_hydro",
+                                 space_subsize_self_hydro_default);
+    space_subsize_pair_stars =
+        parser_get_opt_param_int(params, "Scheduler:cell_sub_size_pair_stars",
+                                 space_subsize_pair_stars_default);
+    space_subsize_self_stars =
+        parser_get_opt_param_int(params, "Scheduler:cell_sub_size_self_stars",
+                                 space_subsize_self_stars_default);
+    space_subsize_pair_grav =
+        parser_get_opt_param_int(params, "Scheduler:cell_sub_size_pair_grav",
+                                 space_subsize_pair_grav_default);
+    space_subsize_self_grav =
+        parser_get_opt_param_int(params, "Scheduler:cell_sub_size_self_grav",
+                                 space_subsize_self_grav_default);
+    space_splitsize = parser_get_opt_param_int(
+        params, "Scheduler:cell_split_size", space_splitsize_default);
+    space_subdepth_diff_grav =
+        parser_get_opt_param_int(params, "Scheduler:cell_subdepth_diff_grav",
+                                 space_subdepth_diff_grav_default);
+    space_extra_parts = parser_get_opt_param_int(
+        params, "Scheduler:cell_extra_parts", space_extra_parts_default);
+    space_extra_sparts = parser_get_opt_param_int(
+        params, "Scheduler:cell_extra_sparts", space_extra_sparts_default);
+    space_extra_gparts = parser_get_opt_param_int(
+        params, "Scheduler:cell_extra_gparts", space_extra_gparts_default);
+    space_extra_bparts = parser_get_opt_param_int(
+        params, "Scheduler:cell_extra_bparts", space_extra_bparts_default);
+
+    engine_max_parts_per_ghost =
+        parser_get_opt_param_int(params, "Scheduler:engine_max_parts_per_ghost",
+                                 engine_max_parts_per_ghost_default);
+    engine_max_sparts_per_ghost = parser_get_opt_param_int(
+        params, "Scheduler:engine_max_sparts_per_ghost",
+        engine_max_sparts_per_ghost_default);
+  }
+
   /* Allocate and init the threads. */
   if (swift_memalign("runners", (void **)&e->runners, SWIFT_CACHE_ALIGNMENT,
                      e->nr_threads * sizeof(struct runner)) != 0)
diff --git a/src/feedback/EAGLE/feedback.c b/src/feedback/EAGLE/feedback.c
index cf6848760fb9827718befe48e43ab53cfe77f8ef..69f92fe992a306d018294009abefedb49a7deeae 100644
--- a/src/feedback/EAGLE/feedback.c
+++ b/src/feedback/EAGLE/feedback.c
@@ -105,15 +105,15 @@ double eagle_feedback_energy_fraction(const struct spart* sp,
 
   /* Star properties */
 
-  /* Smoothed metallicity (metal mass fraction) at birth time of the star */
-  const double Z_smooth = sp->chemistry_data.smoothed_metal_mass_fraction_total;
+  /* Metallicity (metal mass fraction) at birth time of the star */
+  const double Z = chemistry_get_total_metal_mass_fraction_for_feedback(sp);
 
   /* Physical density of the gas at the star's birth time */
   const double rho_birth = sp->birth_density;
   const double n_birth = rho_birth * props->rho_to_n_cgs;
 
   /* Calculate f_E */
-  const double Z_term = pow(max(Z_smooth, 1e-6) / Z_0, n_Z);
+  const double Z_term = pow(max(Z, 1e-6) / Z_0, n_Z);
   const double n_term = pow(n_birth / n_0, -n_n);
   const double denonimator = 1. + Z_term * n_term;
 
@@ -333,7 +333,7 @@ INLINE static void evolve_SNIa(const float log10_min_mass,
    * and use updated values for the star's age and timestep in this function */
   if (log10_max_mass > props->log10_SNIa_max_mass_msun) {
 
-    const float Z = sp->chemistry_data.metal_mass_fraction_total;
+    const float Z = chemistry_get_total_metal_mass_fraction_for_feedback(sp);
     const float max_mass = exp10f(props->log10_SNIa_max_mass_msun);
     const float lifetime_Gyr = lifetime_in_Gyr(max_mass, Z, props);
 
@@ -401,6 +401,9 @@ INLINE static void evolve_SNII(float log10_min_mass, float log10_max_mass,
 
   int low_imf_mass_bin_index, high_imf_mass_bin_index, mass_bin_index;
 
+  /* Metallicity (metal mass fraction) at birth time of the star */
+  const double Z = chemistry_get_total_metal_mass_fraction_for_feedback(sp);
+
   /* If mass at beginning of step is less than tabulated lower bound for IMF,
    * limit it.*/
   if (log10_min_mass < props->log10_SNII_min_mass_msun)
@@ -423,9 +426,7 @@ INLINE static void evolve_SNII(float log10_min_mass, float log10_max_mass,
   int iz_low = 0, iz_high = 0, low_index_3d, high_index_3d, low_index_2d,
       high_index_2d;
   float dz = 0.;
-  determine_bin_yield_SNII(&iz_low, &iz_high, &dz,
-                           log10(sp->chemistry_data.metal_mass_fraction_total),
-                           props);
+  determine_bin_yield_SNII(&iz_low, &iz_high, &dz, log10(Z), props);
 
   /* compute metals produced */
   float metal_mass_released[chemistry_element_count], metal_mass_released_total;
@@ -557,6 +558,9 @@ INLINE static void evolve_AGB(const float log10_min_mass, float log10_max_mass,
 
   int low_imf_mass_bin_index, high_imf_mass_bin_index, mass_bin_index;
 
+  /* Metallicity (metal mass fraction) at birth time of the star */
+  const double Z = chemistry_get_total_metal_mass_fraction_for_feedback(sp);
+
   /* If mass at end of step is greater than tabulated lower bound for IMF, limit
    * it.*/
   if (log10_max_mass > props->log10_SNII_min_mass_msun)
@@ -574,9 +578,7 @@ INLINE static void evolve_AGB(const float log10_min_mass, float log10_max_mass,
   int iz_low = 0, iz_high = 0, low_index_3d, high_index_3d, low_index_2d,
       high_index_2d;
   float dz = 0.f;
-  determine_bin_yield_AGB(&iz_low, &iz_high, &dz,
-                          log10(sp->chemistry_data.metal_mass_fraction_total),
-                          props);
+  determine_bin_yield_AGB(&iz_low, &iz_high, &dz, log10(Z), props);
 
   /* compute metals produced */
   float metal_mass_released[chemistry_element_count], metal_mass_released_total;
@@ -718,7 +720,7 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props,
   const double star_age_Gyr = age * conversion_factor;
 
   /* Get the metallicity */
-  const float Z = sp->chemistry_data.metal_mass_fraction_total;
+  const float Z = chemistry_get_total_metal_mass_fraction_for_feedback(sp);
 
   /* Properties collected in the stellar density loop. */
   const float ngb_gas_mass = sp->feedback_data.to_collect.ngb_mass;
diff --git a/src/hydro/AnarchyDU/hydro_iact.h b/src/hydro/AnarchyDU/hydro_iact.h
index 19489f7460753b15961603e68c22c039f6de27d3..09ecde01d44cfef1c0edbfe0d3f1aa909d553a95 100644
--- a/src/hydro/AnarchyDU/hydro_iact.h
+++ b/src/hydro/AnarchyDU/hydro_iact.h
@@ -428,6 +428,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   /* Get the time derivative for h. */
   pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
   pj->force.h_dt -= mi * dvdr * r_inv / rhoi * wj_dr;
+
+  /* Update if we need to; this should be guaranteed by the gradient loop but
+   * due to some possible synchronisation problems this is here as a _quick
+   * fix_. Added: 14th August 2019. To be removed by 1st Jan 2020. (JB) */
+  pi->viscosity.v_sig = max(pi->viscosity.v_sig, v_sig);
+  pj->viscosity.v_sig = max(pj->viscosity.v_sig, v_sig);
 }
 
 /**
@@ -547,6 +553,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Get the time derivative for h. */
   pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
+
+  /* Update if we need to; this should be guaranteed by the gradient loop but
+   * due to some possible synchronisation problems this is here as a _quick
+   * fix_. Added: 14th August 2019. To be removed by 1st Jan 2020. (JB) */
+  pi->viscosity.v_sig = max(pi->viscosity.v_sig, v_sig);
 }
 
 /**
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index afc9c17d5b6ed691685bf42df2562ecb6f4409d8..93c0e58068341ab29a7ed0fb6b0d21f76980d952 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -41,6 +41,7 @@
 #include "hydro_space.h"
 #include "kernel_hydro.h"
 #include "minmax.h"
+#include "pressure_floor.h"
 
 #include "./hydro_parameters.h"
 
@@ -109,7 +110,8 @@ hydro_get_drifted_physical_internal_energy(const struct part *restrict p,
 __attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
     const struct part *restrict p) {
 
-  return gas_pressure_from_entropy(p->rho, p->entropy);
+  const float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+  return pressure_floor_get_pressure(p, p->rho, comoving_pressure);
 }
 
 /**
@@ -121,7 +123,10 @@ __attribute__((always_inline)) INLINE static float hydro_get_comoving_pressure(
 __attribute__((always_inline)) INLINE static float hydro_get_physical_pressure(
     const struct part *restrict p, const struct cosmology *cosmo) {
 
-  return gas_pressure_from_entropy(p->rho * cosmo->a3_inv, p->entropy);
+  const float phys_pressure =
+      gas_pressure_from_entropy(p->rho * cosmo->a3_inv, p->entropy);
+  const float phys_rho = hydro_get_physical_density(p, cosmo);
+  return pressure_floor_get_pressure(p, phys_rho, phys_pressure);
 }
 
 /**
@@ -388,13 +393,15 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
   const float rho_inv = 1.f / p->rho;
 
   /* Compute the pressure */
-  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+  float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+  comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure);
 
   /* Compute the sound speed */
-  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+  const float soundspeed =
+      gas_soundspeed_from_pressure(p->rho, comoving_pressure);
 
   /* Divide the pressure by the density squared to get the SPH term */
-  const float P_over_rho2 = pressure * rho_inv * rho_inv;
+  const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv;
 
   /* Update variables. */
   p->force.P_over_rho2 = P_over_rho2;
@@ -591,13 +598,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   const float abs_div_physical_v = fabsf(div_physical_v);
 
   /* Compute the pressure */
-  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+  float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+  comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure);
 
   /* Compute the sound speed */
-  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+  const float soundspeed =
+      gas_soundspeed_from_pressure(p->rho, comoving_pressure);
 
   /* Divide the pressure by the density squared to get the SPH term */
-  const float P_over_rho2 = pressure * rho_inv * rho_inv;
+  const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv;
 
   /* Compute the Balsara switch */
   /* Pre-multiply in the AV factor; hydro_props are not passed to the iact
@@ -665,14 +674,16 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
   p->entropy = xp->entropy_full;
 
   /* Re-compute the pressure */
-  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+  float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+  comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure);
 
   /* Compute the new sound speed */
-  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+  const float soundspeed =
+      gas_soundspeed_from_pressure(p->rho, comoving_pressure);
 
   /* Divide the pressure by the density squared to get the SPH term */
   const float rho_inv = 1.f / p->rho;
-  const float P_over_rho2 = pressure * rho_inv * rho_inv;
+  const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv;
 
   /* Update variables */
   p->force.soundspeed = soundspeed;
@@ -730,14 +741,16 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
     p->rho *= expf(w2);
 
   /* Re-compute the pressure */
-  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+  float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+  comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure);
 
   /* Compute the new sound speed */
-  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+  const float soundspeed =
+      gas_soundspeed_from_pressure(p->rho, comoving_pressure);
 
   /* Divide the pressure by the density squared to get the SPH term */
   const float rho_inv = 1.f / p->rho;
-  const float P_over_rho2 = pressure * rho_inv * rho_inv;
+  const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv;
 
   /* Update variables */
   p->force.soundspeed = soundspeed;
@@ -840,14 +853,16 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
   }
 
   /* Compute the pressure */
-  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+  float comoving_pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+  comoving_pressure = pressure_floor_get_pressure(p, p->rho, comoving_pressure);
 
   /* Compute the sound speed */
-  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+  const float soundspeed =
+      gas_soundspeed_from_pressure(p->rho, comoving_pressure);
 
   /* Divide the pressure by the density squared to get the SPH term */
   const float rho_inv = 1.f / p->rho;
-  const float P_over_rho2 = pressure * rho_inv * rho_inv;
+  const float P_over_rho2 = comoving_pressure * rho_inv * rho_inv;
 
   p->force.soundspeed = soundspeed;
   p->force.P_over_rho2 = P_over_rho2;
diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index e3a4f98160150717e50651760d8af68af2a70662..57e0091a7d79fb5404a30c73bca31b2c555e01b1 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -46,6 +46,7 @@
 #include "hydro_space.h"
 #include "kernel_hydro.h"
 #include "minmax.h"
+#include "pressure_floor.h"
 
 #include "./hydro_parameters.h"
 
@@ -221,7 +222,9 @@ hydro_get_comoving_soundspeed(const struct part *restrict p) {
 
   /* Compute the sound speed -- see theory section for justification */
   /* IDEAL GAS ONLY -- P-U does not work with generic EoS. */
-  const float square_rooted = sqrtf(hydro_gamma * p->pressure_bar / p->rho);
+  const float comoving_pressure =
+      pressure_floor_get_pressure(p, p->rho, p->pressure_bar);
+  const float square_rooted = sqrtf(hydro_gamma * comoving_pressure / p->rho);
 
   return square_rooted;
 }
@@ -236,7 +239,10 @@ __attribute__((always_inline)) INLINE static float
 hydro_get_physical_soundspeed(const struct part *restrict p,
                               const struct cosmology *cosmo) {
 
-  return cosmo->a_factor_sound_speed * p->force.soundspeed;
+  const float phys_rho = hydro_get_physical_density(p, cosmo);
+
+  return pressure_floor_get_pressure(
+      p, phys_rho, cosmo->a_factor_sound_speed * p->force.soundspeed);
 }
 
 /**
@@ -497,7 +503,6 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
   p->density.wcount = 0.f;
   p->density.wcount_dh = 0.f;
   p->rho = 0.f;
-  p->density.rho_dh = 0.f;
   p->pressure_bar = 0.f;
   p->density.pressure_bar_dh = 0.f;
 
@@ -531,7 +536,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
 
   /* Final operation on the density (add self-contribution). */
   p->rho += p->mass * kernel_root;
-  p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
   p->pressure_bar += p->mass * p->u * kernel_root;
   p->density.pressure_bar_dh -= hydro_dimension * p->mass * p->u * kernel_root;
   p->density.wcount += kernel_root;
@@ -539,7 +543,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
 
   /* Finish the calculation by inserting the missing h-factors */
   p->rho *= h_inv_dim;
-  p->density.rho_dh *= h_inv_dim_plus_one;
   p->pressure_bar *= (h_inv_dim * hydro_gamma_minus_one);
   p->density.pressure_bar_dh *= (h_inv_dim_plus_one * hydro_gamma_minus_one);
   p->density.wcount *= h_inv_dim;
@@ -584,7 +587,6 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
   p->pressure_bar =
       p->mass * p->u * hydro_gamma_minus_one * kernel_root * h_inv_dim;
   p->density.wcount = kernel_root * h_inv_dim;
-  p->density.rho_dh = 0.f;
   p->density.wcount_dh = 0.f;
   p->density.pressure_bar_dh = 0.f;
 
@@ -640,10 +642,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
                              hydro_one_over_gamma_minus_one) /
                             (1.f + common_factor * p->density.wcount_dh);
 
+  /* Get the pressures */
+  const float comoving_pressure_with_floor =
+      pressure_floor_get_pressure(p, p->rho, p->pressure_bar);
+
   /* Update variables. */
   p->force.f = grad_h_term;
   p->force.soundspeed = soundspeed;
   p->force.balsara = balsara;
+  p->force.pressure_bar_with_floor = comoving_pressure_with_floor;
 }
 
 /**
@@ -755,6 +762,11 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
   const float soundspeed = hydro_get_comoving_soundspeed(p);
 
   p->force.soundspeed = soundspeed;
+
+  /* update the required variables */
+  const float comoving_pressure_with_floor =
+      pressure_floor_get_pressure(p, p->rho, p->pressure_bar);
+  p->force.pressure_bar_with_floor = comoving_pressure_with_floor;
 }
 
 /**
diff --git a/src/hydro/PressureEnergy/hydro_debug.h b/src/hydro/PressureEnergy/hydro_debug.h
index 7ffc370ed4d6abd273fc3d8d5b887f5ccf8e001c..861b16299f824abbee759c80e79d99d449a348c5 100644
--- a/src/hydro/PressureEnergy/hydro_debug.h
+++ b/src/hydro/PressureEnergy/hydro_debug.h
@@ -30,13 +30,13 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       "x=[%.3e,%.3e,%.3e], "
       "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
       "u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
-      "h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, \n"
+      "h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, rho=%.3e, \n"
       "p_dh=%.3e, p_bar=%.3e \n"
       "time_bin=%d wakeup=%d\n",
       p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
       xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
       p->u, p->u_dt, p->force.v_sig, hydro_get_comoving_pressure(p), p->h,
-      p->force.h_dt, (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho,
+      p->force.h_dt, (int)p->density.wcount, p->mass, p->rho,
       p->density.pressure_bar_dh, p->pressure_bar, p->time_bin, p->wakeup);
 }
 
diff --git a/src/hydro/PressureEnergy/hydro_iact.h b/src/hydro/PressureEnergy/hydro_iact.h
index b59fdbbe418d0442fd290b3596fd08531fae88e6..62106bb35cd09a0d67d56fb7b8ab58ba7cac7491 100644
--- a/src/hydro/PressureEnergy/hydro_iact.h
+++ b/src/hydro/PressureEnergy/hydro_iact.h
@@ -69,7 +69,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   kernel_deval(ui, &wi, &wi_dx);
 
   pi->rho += mj * wi;
-  pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
 
   pi->pressure_bar += mj * wi * pj->u;
   pi->density.pressure_bar_dh -=
@@ -83,7 +82,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   kernel_deval(uj, &wj, &wj_dx);
 
   pj->rho += mi * wj;
-  pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
   pj->pressure_bar += mi * wj * pi->u;
   pj->density.pressure_bar_dh -=
       mi * pi->u * (hydro_dimension * wj + uj * wj_dx);
@@ -149,7 +147,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   kernel_deval(ui, &wi, &wi_dx);
 
   pi->rho += mj * wi;
-  pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
 
   pi->pressure_bar += mj * wi * pj->u;
 
@@ -259,11 +256,18 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   /* Convolve with the kernel */
   const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
 
+  /* Compute the ratio of pressures */
+  const float pressure_inverse_i =
+      pi->force.pressure_bar_with_floor / (pi->pressure_bar * pi->pressure_bar);
+  const float pressure_inverse_j =
+      pj->force.pressure_bar_with_floor / (pj->pressure_bar * pj->pressure_bar);
+
   /* SPH acceleration term */
-  const float sph_acc_term =
-      pj->u * pi->u * hydro_gamma_minus_one * hydro_gamma_minus_one *
-      ((f_ij / pi->pressure_bar) * wi_dr + (f_ji / pj->pressure_bar) * wj_dr) *
-      r_inv;
+  const float sph_acc_term = pj->u * pi->u * hydro_gamma_minus_one *
+                             hydro_gamma_minus_one *
+                             ((f_ij * pressure_inverse_i) * wi_dr +
+                              (f_ji * pressure_inverse_j) * wj_dr) *
+                             r_inv;
 
   /* Assemble the acceleration */
   const float acc = sph_acc_term + visc_acc_term;
@@ -278,11 +282,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   pj->a_hydro[2] += mi * acc * dx[2];
 
   /* Get the time derivative for u. */
+
   const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one *
-                              pj->u * pi->u * (f_ij / pi->pressure_bar) *
+                              pj->u * pi->u * (f_ij * pressure_inverse_i) *
                               wi_dr * dvdr * r_inv;
+
   const float sph_du_term_j = hydro_gamma_minus_one * hydro_gamma_minus_one *
-                              pi->u * pj->u * (f_ji / pj->pressure_bar) *
+                              pi->u * pj->u * (f_ji * pressure_inverse_j) *
                               wj_dr * dvdr * r_inv;
 
   /* Viscosity term */
@@ -386,11 +392,18 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   /* Convolve with the kernel */
   const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
 
+  /* Compute the ratio of pressures */
+  const float pressure_inverse_i =
+      pi->force.pressure_bar_with_floor / (pi->pressure_bar * pi->pressure_bar);
+  const float pressure_inverse_j =
+      pj->force.pressure_bar_with_floor / (pj->pressure_bar * pj->pressure_bar);
+
   /* SPH acceleration term */
-  const float sph_acc_term =
-      pj->u * pi->u * hydro_gamma_minus_one * hydro_gamma_minus_one *
-      ((f_ij / pi->pressure_bar) * wi_dr + (f_ji / pj->pressure_bar) * wj_dr) *
-      r_inv;
+  const float sph_acc_term = pj->u * pi->u * hydro_gamma_minus_one *
+                             hydro_gamma_minus_one *
+                             ((f_ij * pressure_inverse_i) * wi_dr +
+                              (f_ji * pressure_inverse_j) * wj_dr) *
+                             r_inv;
 
   /* Assemble the acceleration */
   const float acc = sph_acc_term + visc_acc_term;
@@ -402,7 +415,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Get the time derivative for u. */
   const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one *
-                              pj->u * pi->u * (f_ij / pi->pressure_bar) *
+                              pj->u * pi->u * (f_ij * pressure_inverse_i) *
                               wi_dr * dvdr * r_inv;
 
   /* Viscosity term */
diff --git a/src/hydro/PressureEnergy/hydro_part.h b/src/hydro/PressureEnergy/hydro_part.h
index ded83d70176409332b669517ed203fdf95ecfd5e..500569533a5328f3cefb681f5dfbfaf8ab4460c0 100644
--- a/src/hydro/PressureEnergy/hydro_part.h
+++ b/src/hydro/PressureEnergy/hydro_part.h
@@ -132,9 +132,6 @@ struct part {
       /*! Derivative of the neighbour number with respect to h. */
       float wcount_dh;
 
-      /*! Derivative of density with respect to h */
-      float rho_dh;
-
       /*! Derivative of the weighted pressure with respect to h */
       float pressure_bar_dh;
 
@@ -168,6 +165,9 @@ struct part {
 
       /*! Balsara switch */
       float balsara;
+
+      /*! Pressure term including the pressure floor */
+      float pressure_bar_with_floor;
     } force;
   };
 
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 0082546ce14030a37fa63d87ea88bc99153b8213..d272fb7381664d91f763119ae68043ce59a77e1f 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -34,6 +34,7 @@
 #include "hydro.h"
 #include "kernel_hydro.h"
 #include "parser.h"
+#include "pressure_floor.h"
 #include "units.h"
 
 #define hydro_props_default_max_iterations 30
@@ -186,6 +187,9 @@ void hydro_props_print(const struct hydro_props *p) {
   /* Print equation of state first */
   eos_print(&eos);
 
+  /* Then the pressure floor */
+  pressure_floor_print(&pressure_floor_props);
+
   /* Now SPH */
   message("Hydrodynamic scheme: %s in %dD.", SPH_IMPLEMENTATION,
           (int)hydro_dimension);
@@ -238,6 +242,7 @@ void hydro_props_print(const struct hydro_props *p) {
 void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
 
   eos_print_snapshot(h_grpsph, &eos);
+  pressure_floor_print_snapshot(h_grpsph);
 
   io_write_attribute_i(h_grpsph, "Dimension", (int)hydro_dimension);
   io_write_attribute_s(h_grpsph, "Scheme", SPH_IMPLEMENTATION);
diff --git a/src/pressure_floor.c b/src/pressure_floor.c
new file mode 100644
index 0000000000000000000000000000000000000000..2901241f1a2ff42d526a0122008dcf6ab2d8ea6f
--- /dev/null
+++ b/src/pressure_floor.c
@@ -0,0 +1,24 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* This object's header. */
+#include "pressure_floor.h"
+
+/* Pressure floor for the physics model. */
+struct pressure_floor_properties pressure_floor_props;
diff --git a/src/pressure_floor.h b/src/pressure_floor.h
new file mode 100644
index 0000000000000000000000000000000000000000..4389dfe0891dd6bfbc1e6730cac1c6ddcc920547
--- /dev/null
+++ b/src/pressure_floor.h
@@ -0,0 +1,54 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_PRESSURE_FLOOR_H
+#define SWIFT_PRESSURE_FLOOR_H
+
+/**
+ * @file src/pressure_floor.h
+ * @brief Branches between the different pressure floor models
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local includes */
+#include "common_io.h"
+#include "cosmology.h"
+#include "error.h"
+#include "inline.h"
+
+extern struct pressure_floor_properties pressure_floor_props;
+
+/* Check if pressure floor is implemented in hydro */
+#ifndef PRESSURE_FLOOR_NONE
+#if defined(GADGET2_SPH) || defined(HOPKINS_PU_SPH)
+/* Implemented */
+#else
+#error Pressure floor not implemented with this hydro scheme
+#endif
+
+#endif
+/* Import the right pressure floor definition */
+#if defined(PRESSURE_FLOOR_NONE)
+#include "./pressure_floor/none/pressure_floor.h"
+#elif defined(PRESSURE_FLOOR_GEAR)
+#include "./pressure_floor/GEAR/pressure_floor.h"
+#endif
+
+#endif /* SWIFT_PRESSURE_FLOOR_H */
diff --git a/src/pressure_floor/GEAR/pressure_floor.h b/src/pressure_floor/GEAR/pressure_floor.h
new file mode 100644
index 0000000000000000000000000000000000000000..dd715c155a095f9ad5f97b1d14cabfc94d9b11b0
--- /dev/null
+++ b/src/pressure_floor/GEAR/pressure_floor.h
@@ -0,0 +1,124 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_PRESSURE_FLOOR_GEAR_H
+#define SWIFT_PRESSURE_FLOOR_GEAR_H
+
+#include "adiabatic_index.h"
+#include "cosmology.h"
+#include "equation_of_state.h"
+#include "hydro_properties.h"
+#include "parser.h"
+#include "part.h"
+#include "units.h"
+
+/**
+ * @file src/pressure_floor/GEAR/pressure_floor.h
+ * @brief Pressure floor used in the GEAR model
+ */
+
+/**
+ * @brief Properties of the pressure floor in the GEAR model.
+ */
+struct pressure_floor_properties {
+
+  /*! Jeans factor */
+  float n_jeans;
+
+  /*! The constants in internal units (4 G N_jeans^(2/3) / PI) */
+  float constants;
+};
+
+/**
+ * @brief Compute the physical pressure floor of a given #part.
+ *
+ * Note that the particle is not updated!!
+ *
+ * The inputs can be either in physical or comoving coordinates (the output is
+ * in the same coordinates).
+ *
+ * @param p The #part.
+ * @param rho The physical or comoving density.
+ * @param pressure The physical or comoving pressure without any pressure floor.
+ *
+ * @return The physical or comoving pressure with the floor.
+ */
+static INLINE float pressure_floor_get_pressure(const struct part *p,
+                                                const float rho,
+                                                const float pressure) {
+
+  /* Compute pressure floor */
+  float floor = p->h * p->h * rho * pressure_floor_props.constants;
+  // TODO add sigma (will be done once the SF is merged)
+  floor *= rho * hydro_one_over_gamma;
+
+  return fmax(pressure, floor);
+}
+
+/**
+ * @brief Initialise the pressure floor by reading the parameters and converting
+ * to internal units.
+ *
+ * The input temperatures and number densities are converted to pressure and
+ * density assuming a neutral gas of primoridal abundance.
+ *
+ * @param params The YAML parameter file.
+ * @param us The system of units used internally.
+ * @param phys_const The physical constants.
+ * @param hydro_props The propoerties of the hydro scheme.
+ * @param props The pressure floor properties to fill.
+ */
+static INLINE void pressure_floor_init(struct pressure_floor_properties *props,
+                                       const struct phys_const *phys_const,
+                                       const struct unit_system *us,
+                                       const struct hydro_props *hydro_props,
+                                       struct swift_params *params) {
+
+  /* Read the Jeans factor */
+  props->n_jeans =
+      parser_get_param_float(params, "GEARPressureFloor:Jeans_factor");
+
+  /* Compute the constants */
+  props->constants =
+      4.0 * M_1_PI * phys_const->const_newton_G * pow(props->n_jeans, 2. / 3.);
+}
+
+/**
+ * @brief Print the properties of the pressure floor to stdout.
+ *
+ * @param props The pressure floor properties.
+ */
+static INLINE void pressure_floor_print(
+    const struct pressure_floor_properties *props) {
+
+  message("Pressure floor is 'GEAR' with:");
+  message("Jeans factor: %g", props->n_jeans);
+}
+
+#ifdef HAVE_HDF5
+
+/**
+ * @brief Writes the current model of pressure floor to the file
+ * @param h_grp The HDF5 group in which to write
+ */
+INLINE static void pressure_floor_print_snapshot(hid_t h_grp) {
+
+  io_write_attribute_s(h_grp, "Pressure floor", "GEAR");
+}
+#endif
+#endif /* SWIFT_PRESSURE_FLOOR_GEAR_H */
diff --git a/src/pressure_floor/none/pressure_floor.h b/src/pressure_floor/none/pressure_floor.h
new file mode 100644
index 0000000000000000000000000000000000000000..2d0b95a71f9ccce6697e79760aa43b560933e7bd
--- /dev/null
+++ b/src/pressure_floor/none/pressure_floor.h
@@ -0,0 +1,99 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_PRESSURE_FLOOR_NONE_H
+#define SWIFT_PRESSURE_FLOOR_NONE_H
+
+#include "adiabatic_index.h"
+#include "cosmology.h"
+#include "equation_of_state.h"
+#include "hydro_properties.h"
+#include "parser.h"
+#include "part.h"
+#include "units.h"
+
+/**
+ * @file src/pressure_floor/none/pressure_floor.h
+ * @brief Model without pressure floor
+ */
+
+/**
+ * @brief Properties of the pressure floor in the NONE model.
+ */
+struct pressure_floor_properties {};
+
+/**
+ * @brief Compute the physical pressure floor of a given #part.
+ *
+ * Note that the particle is not updated!!
+ *
+ * The inputs can be either in physical or comoving coordinates (the output is
+ * in the same coordinates).
+ *
+ * @param p The #part.
+ * @param rho The physical or comoving density.
+ * @param pressure The physical or comoving pressure without any pressure floor.
+ *
+ * @return The physical or comoving pressure with the floor.
+ */
+static INLINE float pressure_floor_get_pressure(const struct part *p,
+                                                const float rho,
+                                                const float pressure) {
+  return pressure;
+}
+
+/**
+ * @brief Initialise the pressure floor by reading the parameters and converting
+ * to internal units.
+ *
+ * The input temperatures and number densities are converted to pressure and
+ * density assuming a neutral gas of primoridal abundance.
+ *
+ * @param params The YAML parameter file.
+ * @param us The system of units used internally.
+ * @param phys_const The physical constants.
+ * @param hydro_props The propoerties of the hydro scheme.
+ * @param props The pressure floor properties to fill.
+ */
+static INLINE void pressure_floor_init(struct pressure_floor_properties *props,
+                                       const struct phys_const *phys_const,
+                                       const struct unit_system *us,
+                                       const struct hydro_props *hydro_props,
+                                       struct swift_params *params) {}
+
+/**
+ * @brief Print the properties of the pressure floor to stdout.
+ *
+ * @param props The pressure floor properties.
+ */
+static INLINE void pressure_floor_print(
+    const struct pressure_floor_properties *props) {}
+
+#ifdef HAVE_HDF5
+
+/**
+ * @brief Writes the current model of pressure floor to the file
+ * @param h_grp The HDF5 group in which to write
+ */
+INLINE static void pressure_floor_print_snapshot(hid_t h_grp) {
+
+  io_write_attribute_s(h_grp, "Pressure floor", "none");
+}
+#endif
+
+#endif /* SWIFT_PRESSURE_FLOOR_NONE_H */
diff --git a/src/runner.c b/src/runner.c
index cf3d15d72d647dee8ca3518dec8fa598b4f71075..318a9143c83bdc1b56e4d7876900fd998e3ebb5a 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -905,6 +905,10 @@ void runner_do_black_holes_swallow_ghost(struct runner *r, struct cell *c,
 
       if (bpart_is_active(bp, e)) {
 
+        /* Compute the final operations for repositioning of this BH */
+        black_holes_end_reposition(bp, e->black_holes_properties,
+                                   e->physical_constants, e->cosmology);
+
         /* Get particle time-step */
         double dt;
         if (with_cosmology) {
diff --git a/src/runner_doiact_black_holes.h b/src/runner_doiact_black_holes.h
index 71ce7aefcfc3689cc8476884729e4e7cd77bef66..ce159c7ac24a508bc625070ed50b3aad7dd9fa8d 100644
--- a/src/runner_doiact_black_holes.h
+++ b/src/runner_doiact_black_holes.h
@@ -159,7 +159,8 @@ void DOSELF1_BH(struct runner *r, struct cell *c, int timer) {
 #endif
 
         if (r2 < hig2) {
-          IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, cosmo, ti_current);
+          IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, cosmo, e->gravity_properties,
+                      ti_current);
         }
       } /* loop over the parts in ci. */
     }   /* loop over the bparts in ci. */
@@ -310,7 +311,8 @@ void DO_NONSYM_PAIR1_BH_NAIVE(struct runner *r, struct cell *restrict ci,
 #endif
 
         if (r2 < hig2) {
-          IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, cosmo, ti_current);
+          IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, cosmo, e->gravity_properties,
+                      ti_current);
         }
       } /* loop over the parts in cj. */
     }   /* loop over the bparts in ci. */
@@ -475,7 +477,8 @@ void DOPAIR1_SUBSET_BH_NAIVE(struct runner *r, struct cell *restrict ci,
 #endif
       /* Hit or miss? */
       if (r2 < hig2) {
-        IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, cosmo, ti_current);
+        IACT_BH_GAS(r2, dx, hi, hj, bi, pj, xpj, cosmo, e->gravity_properties,
+                    ti_current);
       }
     } /* loop over the parts in cj. */
   }   /* loop over the parts in ci. */
@@ -550,7 +553,8 @@ void DOSELF1_SUBSET_BH(struct runner *r, struct cell *restrict ci,
 
       /* Hit or miss? */
       if (r2 < hig2) {
-        IACT_BH_GAS(r2, dx, hi, pj->h, bi, pj, xpj, cosmo, ti_current);
+        IACT_BH_GAS(r2, dx, hi, pj->h, bi, pj, xpj, cosmo,
+                    e->gravity_properties, ti_current);
       }
     } /* loop over the parts in cj. */
   }   /* loop over the parts in ci. */
diff --git a/src/space.c b/src/space.c
index 250af5efafa506408c9fd91c1adffb7cd96a1a21..d8e37fd5a43a151daf8501b9011922983244d99e 100644
--- a/src/space.c
+++ b/src/space.c
@@ -4591,8 +4591,9 @@ void space_init(struct space *s, struct swift_params *params,
     Ngpart = s->nr_gparts;
 
 #ifdef SWIFT_DEBUG_CHECKS
-    part_verify_links(parts, gparts, sparts, bparts, Npart, Ngpart, Nspart,
-                      Nbpart, 1);
+    if (!dry_run)
+      part_verify_links(parts, gparts, sparts, bparts, Npart, Ngpart, Nspart,
+                        Nbpart, 1);
 #endif
   }
 
diff --git a/src/space.h b/src/space.h
index c294bfae3612c699345bd4a267c40b67cffc2bf1..1b31f0495890a68074225cc68c6acabd1abb53e8 100644
--- a/src/space.h
+++ b/src/space.h
@@ -112,9 +112,6 @@ struct space {
   /*! The minimum top-level cell width allowed. */
   double cell_min;
 
-  /*! Current maximum displacement for particles. */
-  float dx_max;
-
   /*! Space dimensions in number of top-cells. */
   int cdim[3];
 
diff --git a/src/star_formation/EAGLE/star_formation.h b/src/star_formation/EAGLE/star_formation.h
index f2c77e036842b4ca040c58a6bcad1513b03a42bd..bdd6dcb414d2a559e113b54de290fb47b82e38f8 100644
--- a/src/star_formation/EAGLE/star_formation.h
+++ b/src/star_formation/EAGLE/star_formation.h
@@ -232,8 +232,11 @@ INLINE static int star_formation_is_star_forming(
    * because we also need to check if the physical density exceeded
    * the appropriate limit */
 
-  const double Z = p->chemistry_data.smoothed_metal_mass_fraction_total;
-  const double X_H = p->chemistry_data.smoothed_metal_mass_fraction[0];
+  const double Z =
+      chemistry_get_total_metal_mass_fraction_for_star_formation(p);
+  const float* const metal_fraction =
+      chemistry_get_metal_mass_fraction_for_star_formation(p);
+  const double X_H = metal_fraction[chemistry_element_H];
   const double n_H = physical_density * X_H;
 
   /* Get the density threshold */
@@ -279,7 +282,9 @@ INLINE static void star_formation_compute_SFR(
 
   /* Hydrogen number density of this particle */
   const double physical_density = hydro_get_physical_density(p, cosmo);
-  const double X_H = p->chemistry_data.smoothed_metal_mass_fraction[0];
+  const float* const metal_fraction =
+      chemistry_get_metal_mass_fraction_for_star_formation(p);
+  const double X_H = metal_fraction[chemistry_element_H];
   const double n_H = physical_density * X_H / phys_const->const_proton_mass;
 
   /* Are we above the threshold for automatic star formation? */
diff --git a/src/star_formation/EAGLE/star_formation_logger.h b/src/star_formation/EAGLE/star_formation_logger.h
index a843f63ce518aa9be6b1d389d42cf71c8608b903..fea15242e9a83b81cf5c7801c0546a1bbb5a8eae 100644
--- a/src/star_formation/EAGLE/star_formation_logger.h
+++ b/src/star_formation/EAGLE/star_formation_logger.h
@@ -100,13 +100,13 @@ INLINE static void star_formation_logger_add_to_accumulator(
     const struct star_formation_history *sf_add) {
 
   /* Update the SFH structure */
-  sf_update->new_stellar_mass += sf_add->new_stellar_mass;
+  sf_update->new_stellar_mass = sf_add->new_stellar_mass;
 
-  sf_update->SFR_active += sf_add->SFR_active;
+  sf_update->SFR_active = sf_add->SFR_active;
 
-  sf_update->SFRdt_active += sf_add->SFRdt_active;
+  sf_update->SFRdt_active = sf_add->SFRdt_active;
 
-  sf_update->SFR_inactive += sf_add->SFR_inactive;
+  sf_update->SFR_inactive = sf_add->SFR_inactive;
 }
 
 /**
diff --git a/src/swift.h b/src/swift.h
index b9f8818d8b833231971abb1afb36ee4507648488..fe9196a8fcf6d1845c9446c480c7961504a4756f 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -55,6 +55,7 @@
 #include "map.h"
 #include "memuse.h"
 #include "mesh_gravity.h"
+#include "minmax.h"
 #include "multipole.h"
 #include "outputlist.h"
 #include "parallel_io.h"
@@ -64,6 +65,7 @@
 #include "periodic.h"
 #include "physical_constants.h"
 #include "potential.h"
+#include "pressure_floor.h"
 #include "profiler.h"
 #include "queue.h"
 #include "random.h"
diff --git a/tests/testRandom.c b/tests/testRandom.c
index a550273c63d1cac048e511c8c48908ec8376441e..036f58b29115a2e646cea35e873ccdc9a4164e4e 100644
--- a/tests/testRandom.c
+++ b/tests/testRandom.c
@@ -255,30 +255,42 @@ int main(int argc, char* argv[]) {
 
     /* Verify that the mean and variance match the expected values for a uniform
      * distribution */
-    const double tolmean = 2e-4;
-    const double tolvar = 1e-3;
-    const double tolcorr = 4e-4;
-
-    if ((fabs(mean - 0.5) / 0.5 > tolmean) ||
-        (fabs(var - 1. / 12.) / (1. / 12.) > tolvar) ||
-        (correlation > tolcorr) || (correlationID > tolcorr) ||
-        (fabs(meanID - 0.5) / 0.5 > tolmean) ||
-        (fabs(varID - 1. / 12.) / (1. / 12.) > tolvar) ||
-        (corr_star_sf > tolcorr) || (corr_star_se > tolcorr) ||
-        (corr_star_bh > tolcorr) || (corr_sf_se > tolcorr) ||
-        (corr_sf_bh > tolcorr) || (corr_se_bh > tolcorr) ||
-        (fabs(mean_sf - 0.5) / 0.5 > tolmean) ||
-        (fabs(mean_se - 0.5) / 0.5 > tolmean) ||
-        (fabs(mean_bh - 0.5) / 0.5 > tolmean) ||
-        (fabs(var_sf - 1. / 12.) / (1. / 12.) > tolvar) ||
-        (fabs(var_se - 1. / 12.) / (1. / 12.) > tolvar) ||
-        (fabs(var_bh - 1. / 12.) / (1. / 12.) > tolvar)) {
+
+    /* Set the allowed standard deviation */
+    const double std_check = 5.;
+
+    /* The mean is expected to deviate a maximum of std_check * std / sqrt(N) */
+    const double tolmean = std_check / sqrtf(12.f * count);
+
+    /* the variance is expected to deviate a maximum of std_check * variance
+     * * sqrt(2/(n-1)) */
+    const double tolvar =
+        std_check * sqrtf(2.f / (12.f * ((double)count - 1.f)));
+
+    /* The correlation coefficient is expected to deviate sqrt(1-R^2)
+     * / sqrt(n-2), in our case <R> = 0, so we get 1/sqrt(n-2) */
+    const double tolcorr = std_check / sqrtf((double)count - 2.);
+
+    if ((fabs(mean - 0.5) > tolmean) || (fabs(var - 1. / 12.) > tolvar) ||
+        (fabs(correlation) > tolcorr) || (fabs(correlationID) > tolcorr) ||
+        (fabs(meanID - 0.5) > tolmean) || (fabs(varID - 1. / 12.) > tolvar) ||
+        (fabs(corr_star_sf) > tolcorr) || (fabs(corr_star_se) > tolcorr) ||
+        (fabs(corr_star_bh) > tolcorr) || (fabs(corr_sf_se) > tolcorr) ||
+        (fabs(corr_sf_bh) > tolcorr) || (fabs(corr_se_bh) > tolcorr) ||
+        (fabs(mean_sf - 0.5) > tolmean) || (fabs(mean_se - 0.5) > tolmean) ||
+        (fabs(mean_bh - 0.5) > tolmean) || (fabs(var_sf - 1. / 12.) > tolvar) ||
+        (fabs(var_se - 1. / 12.) > tolvar) ||
+        (fabs(var_bh - 1. / 12.) > tolvar)) {
       message("Test failed!");
       message("Global result:");
       message("Result:    count=%d mean=%f var=%f, correlation=%f", count, mean,
               var, correlation);
       message("Expected:  count=%d mean=%f var=%f, correlation=%f", count, 0.5f,
               1. / 12., 0.);
+      message("Max difference:           mean=%f var=%f, correlation=%f",
+              tolmean, tolvar, tolcorr);
+      message("Difference:               mean=%f var=%f, correlation=%f",
+              fabs(mean - 0.5f), fabs(var - 1. / 12.), fabs(correlation));
       message("ID part");
       message(
           "Result:    count=%d mean=%f var=%f"
@@ -288,23 +300,47 @@ int main(int argc, char* argv[]) {
           "Expected:  count=%d mean=%f var=%f"
           " correlation=%f",
           count, .5f, 1. / 12., 0.);
+      message("Max difference:           mean=%f var=%f, correlation=%f",
+              tolmean, tolvar, tolcorr);
+      message("Difference:               mean=%f var=%f, correlation=%f",
+              fabs(meanID - 0.5f), fabs(varID - 1. / 12.), fabs(correlation));
       message("Different physical processes:");
       message(
           "Means:    stars=%f stellar feedback=%f stellar "
-          " enrichement=%f black holes=%f",
+          " enrichment=%f black holes=%f",
           mean, mean_sf, mean_se, mean_bh);
       message(
           "Expected: stars=%f stellar feedback=%f stellar "
-          " enrichement=%f black holes=%f",
+          " enrichment=%f black holes=%f",
           .5f, .5f, .5f, .5f);
+      message(
+          "Max diff: stars=%f stellar feedback=%f stellar "
+          " enrichment=%f black holes=%f",
+          tolmean, tolmean, tolmean, tolmean);
+      message(
+          "Diff:     stars=%f stellar feedback=%f stellar "
+          " enrichment=%f black holes=%f",
+          fabs(mean - .5f), fabs(mean_sf - .5f), fabs(mean_se - .5f),
+          fabs(mean_bh - .5f));
+      message(" ");
       message(
           "Var:      stars=%f stellar feedback=%f stellar "
-          " enrichement=%f black holes=%f",
+          " enrichment=%f black holes=%f",
           var, var_sf, var_se, var_bh);
       message(
           "Expected: stars=%f stellar feedback=%f stellar "
-          " enrichement=%f black holes=%f",
+          " enrichment=%f black holes=%f",
           1. / 12., 1. / 12., 1 / 12., 1. / 12.);
+      message(
+          "Max diff: stars=%f stellar feedback=%f stellar "
+          " enrichment=%f black holes=%f",
+          tolvar, tolvar, tolvar, tolvar);
+      message(
+          "Diff:     stars=%f stellar feedback=%f stellar "
+          " enrichment=%f black holes=%f",
+          fabs(var - 1. / 12.), fabs(var_sf - 1. / 12.),
+          fabs(var_se - 1. / 12.), fabs(var_bh - 1. / 12.));
+      message(" ");
       message(
           "Correlation: stars-sf=%f stars-se=%f stars-bh=%f "
           "sf-se=%f sf-bh=%f se-bh=%f",
@@ -314,6 +350,15 @@ int main(int argc, char* argv[]) {
           "Expected:    stars-sf=%f stars-se=%f stars-bh=%f "
           "sf-se=%f sf-bh=%f se-bh=%f",
           0., 0., 0., 0., 0., 0.);
+      message(
+          "Max diff:    stars-sf=%f stars-se=%f stars-bh=%f "
+          "sf-se=%f sf-bh=%f se-bh=%f",
+          tolcorr, tolcorr, tolcorr, tolcorr, tolcorr, tolcorr);
+      message(
+          "Diff:        stars-sf=%f stars-se=%f stars-bh=%f "
+          "sf-se=%f sf-bh=%f se-bh=%f",
+          fabs(corr_star_sf), fabs(corr_star_se), fabs(corr_star_bh),
+          fabs(corr_sf_se), fabs(corr_sf_bh), fabs(corr_se_bh));
       return 1;
     }
   }
diff --git a/tests/testRandomSpacing.c b/tests/testRandomSpacing.c
index 0c30ffcf5c8b493665d4cc1153e573a5836f4eec..0d2777ee702458ccaa6170483c48b83ce1a4fc7e 100644
--- a/tests/testRandomSpacing.c
+++ b/tests/testRandomSpacing.c
@@ -75,6 +75,7 @@ int main(int argc, char* argv[]) {
       const double r =
           random_unit_interval(id, ti_current, random_number_star_formation);
 
+      /* Count the number of random numbers below the boundaries */
       if (r < boundary[0]) count[0] += 1;
       if (r < boundary[1]) count[1] += 1;
       if (r < boundary[2]) count[2] += 1;
@@ -83,6 +84,7 @@ int main(int argc, char* argv[]) {
       if (r < boundary[5]) count[5] += 1;
     }
 
+    /* Print counted number of random numbers below the boundaries */
     message("Categories           | %6.0e %6.0e %6.0e %6.0e %6.0e %6.0e",
             boundary[0], boundary[1], boundary[2], boundary[3], boundary[4],
             boundary[5]);
@@ -105,21 +107,40 @@ int main(int argc, char* argv[]) {
     expected_result_int[4] = (int)expected_result[4];
     expected_result_int[5] = (int)expected_result[5];
 
+    /* Print the expected numbers */
     message("expected             | %6d %6d %6d %6d %6d %6d",
             expected_result_int[0], expected_result_int[1],
             expected_result_int[2], expected_result_int[3],
             expected_result_int[4], expected_result_int[5]);
 
     int std_expected_result[6];
-    std_expected_result[0] = (int)max(sqrt(expected_result[0]), 1);
-    std_expected_result[1] = (int)max(sqrt(expected_result[1]), 1);
-    std_expected_result[2] = (int)max(sqrt(expected_result[2]), 1);
-    std_expected_result[3] = (int)max(sqrt(expected_result[3]), 1);
-    std_expected_result[4] = (int)max(sqrt(expected_result[4]), 1);
-    std_expected_result[5] = (int)max(sqrt(expected_result[5]), 1);
 
+    /* Calculate the allowed standard error deviation the maximum of:
+     * 1. the standard error of the expected number doing sqrt(N_expected)
+     * 2. The standard error of the counted number doing sqrt(N_count)
+     * 3. 1 to prevent low number statistics to crash for 1 while expected
+     *    close to zero.
+     *
+     * 1 and 2 are for large numbers essentially the same but for small numbers
+     * it becomes imporatant (e.g. count=6 expected=.9, allowed 5+.9 so 6
+     * fails, but sqrt(6) ~ 2.5 so it should be fine) */
+    std_expected_result[0] =
+        (int)max3(sqrt(expected_result[0]), 1, sqrt(count[0]));
+    std_expected_result[1] =
+        (int)max3(sqrt(expected_result[1]), 1, sqrt(count[1]));
+    std_expected_result[2] =
+        (int)max3(sqrt(expected_result[2]), 1, sqrt(count[2]));
+    std_expected_result[3] =
+        (int)max3(sqrt(expected_result[3]), 1, sqrt(count[3]));
+    std_expected_result[4] =
+        (int)max3(sqrt(expected_result[4]), 1, sqrt(count[4]));
+    std_expected_result[5] =
+        (int)max3(sqrt(expected_result[5]), 1, sqrt(count[5]));
+
+    /* We want 5 sigma (can be changed if necessary) */
     const int numb_sigma = 5;
 
+    /* Print the differences and the 5 sigma differences */
     message("Difference           | %6d %6d %6d %6d %6d %6d",
             abs(expected_result_int[0] - count[0]),
             abs(expected_result_int[1] - count[1]),
@@ -136,6 +157,7 @@ int main(int argc, char* argv[]) {
             numb_sigma * std_expected_result[4],
             numb_sigma * std_expected_result[5]);
 
+    /* Fail if it is not within numb_sigma (5) of the expected difference. */
     if (count[0] >
             expected_result_int[0] + numb_sigma * std_expected_result[0] ||
         count[0] <
diff --git a/theory/SPH/Derivation/sph_derivation.tex b/theory/SPH/Derivation/sph_derivation.tex
index c770035df6050afe836d943dbe6538d4eafc262b..af39fb6833ab7a878d4328f8ba7645b5f07bf7ae 100644
--- a/theory/SPH/Derivation/sph_derivation.tex
+++ b/theory/SPH/Derivation/sph_derivation.tex
@@ -186,7 +186,7 @@ The spatial differential is fairly straightforward and is given by
 The differential with respect to the smoothing length is also straightforward
 after remembering that $W_{ij}(h_j) = w(|r_{ij}|/h_j)/h_j^{n_d}$. Then,
 \begin{align}
-  \frac{\partial y_i}{\partial h_i} = -\sum_{j=1}^N \frac{x_j}{h_j}
+  \frac{\partial y_i}{\partial h_i} = -\sum_{j=1}^N \frac{x_j}{h_i}
   \left[
     n_d W_{ij}(h_i) + \frac{|r_{ij}|}{h_i}
     \left.