diff --git a/src/Makefile.am b/src/Makefile.am
index ce68e2f464953f969eab8b16dcb30a148071a568..46871618a7abd138ceb61c8a0a725c4e185d1fd6 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -65,11 +65,13 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
 nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
 		 kernel_long_gravity.h vector.h cache.h runner_doiact.h runner_doiact_vec.h runner_doiact_grav.h runner_doiact_fft.h \
                  runner_doiact_nosort.h units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
-		 dimension.h equation_of_state.h part_type.h periodic.h memswap.h dump.h logger.h sign.h \
+		 dimension.h part_type.h periodic.h memswap.h dump.h logger.h sign.h \
 		 gravity.h gravity_io.h gravity_cache.h \
 		 gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
 		 gravity/Default/gravity_debug.h gravity/Default/gravity_part.h  \
 		 sourceterms.h \
+		 equation_of_state.h \
+		 equation_of_state/ideal_gas/equation_of_state.h equation_of_state/isothermal/equation_of_state.h \
 	 	 hydro.h hydro_io.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 \
diff --git a/src/equation_of_state.c b/src/equation_of_state.c
index a792bc58cf13fa5f710dc1cd80713ff90e700b18..0953ee7114314f3baf6e86a14a6ddb0e98442dd9 100644
--- a/src/equation_of_state.c
+++ b/src/equation_of_state.c
@@ -20,48 +20,6 @@
 /* This object's header. */
 #include "equation_of_state.h"
 
-/* local headers */
-#include "common_io.h"
-
 /* Equation of state for the physics model
  * (temporary ugly solution as a global variable) */
 struct eos_parameters eos;
-
-void eos_init(struct eos_parameters *e, const struct swift_params *params) {
-
-#if defined(EOS_IDEAL_GAS)
-/* nothing to do here */
-#elif defined(EOS_ISOTHERMAL_GAS)
-  e->isothermal_internal_energy =
-      parser_get_param_float(params, "EoS:isothermal_internal_energy");
-#endif
-}
-
-void eos_print(const struct eos_parameters *e) {
-
-#if defined(EOS_IDEAL_GAS)
-  message("Equation of state: Ideal gas.");
-#elif defined(EOS_ISOTHERMAL_GAS)
-  message(
-      "Equation of state: Isothermal with internal energy "
-      "per unit mass set to %f.",
-      e->isothermal_internal_energy);
-#endif
-
-  message("Adiabatic index gamma: %f.", hydro_gamma);
-}
-
-#if defined(HAVE_HDF5)
-void eos_print_snapshot(hid_t h_grpsph, const struct eos_parameters *e) {
-
-  io_write_attribute_f(h_grpsph, "Adiabatic index", hydro_gamma);
-
-#if defined(EOS_IDEAL_GAS)
-  io_write_attribute_s(h_grpsph, "Equation of state", "Ideal gas");
-#elif defined(EOS_ISOTHERMAL_GAS)
-  io_write_attribute_s(h_grpsph, "Equation of state", "Isothermal gas");
-  io_write_attribute_f(h_grpsph, "Thermal energy per unit mass",
-                       e->isothermal_internal_energy);
-#endif
-}
-#endif
diff --git a/src/equation_of_state.h b/src/equation_of_state.h
index 76b8bb922133c9c4f29453a14068fe3f9044d66f..195b52514f2acc0c40959e09c088a06f0a411869 100644
--- a/src/equation_of_state.h
+++ b/src/equation_of_state.h
@@ -29,315 +29,13 @@
 /* Config parameters. */
 #include "../config.h"
 
-/* Some standard headers. */
-#include <math.h>
-
-/* Local headers. */
-#include "adiabatic_index.h"
-#include "debug.h"
-#include "inline.h"
-
-extern struct eos_parameters eos;
-
-/* ------------------------------------------------------------------------- */
+/* Import the right functions */
 #if defined(EOS_IDEAL_GAS)
-
-struct eos_parameters {};
-
-/**
- * @brief Returns the internal energy given density and entropy
- *
- * Computes \f$u = \frac{S\rho^{\gamma-1} }{\gamma - 1}\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param entropy The entropy \f$S\f$.
- */
-__attribute__((always_inline)) INLINE static float
-gas_internal_energy_from_entropy(float density, float entropy) {
-
-  return entropy * pow_gamma_minus_one(density) *
-         hydro_one_over_gamma_minus_one;
-}
-
-/**
- * @brief Returns the pressure given density and entropy
- *
- * Computes \f$P = S\rho^\gamma\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param entropy The entropy \f$S\f$.
- */
-__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
-    float density, float entropy) {
-
-  return entropy * pow_gamma(density);
-}
-
-/**
- * @brief Returns the entropy given density and pressure.
- *
- * Computes \f$A = \frac{P}{\rho^\gamma}\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param pressure The pressure \f$P\f$.
- * @return The entropy \f$A\f$.
- */
-__attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
-    float density, float pressure) {
-
-  return pressure * pow_minus_gamma(density);
-}
-
-/**
- * @brief Returns the sound speed given density and entropy
- *
- * Computes \f$c = \sqrt{\gamma S \rho^{\gamma-1}}\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param entropy The entropy \f$S\f$.
- */
-__attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
-    float density, float entropy) {
-
-  return sqrtf(hydro_gamma * pow_gamma_minus_one(density) * entropy);
-}
-
-/**
- * @brief Returns the entropy given density and internal energy
- *
- * Computes \f$S = \frac{(\gamma - 1)u}{\rho^{\gamma-1}}\f$.
- *
- * @param density The density \f$\rho\f$
- * @param u The internal energy \f$u\f$
- */
-__attribute__((always_inline)) INLINE static float
-gas_entropy_from_internal_energy(float density, float u) {
-
-  return hydro_gamma_minus_one * u * pow_minus_gamma_minus_one(density);
-}
-
-/**
- * @brief Returns the pressure given density and internal energy
- *
- * Computes \f$P = (\gamma - 1)u\rho\f$.
- *
- * @param density The density \f$\rho\f$
- * @param u The internal energy \f$u\f$
- */
-__attribute__((always_inline)) INLINE static float
-gas_pressure_from_internal_energy(float density, float u) {
-
-  return hydro_gamma_minus_one * u * density;
-}
-
-/**
- * @brief Returns the internal energy given density and pressure.
- *
- * Computes \f$u = \frac{1}{\gamma - 1}\frac{P}{\rho}\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param pressure The pressure \f$P\f$.
- * @return The internal energy \f$u\f$.
- */
-__attribute__((always_inline)) INLINE static float
-gas_internal_energy_from_pressure(float density, float pressure) {
-  return hydro_one_over_gamma_minus_one * pressure / density;
-}
-
-/**
- * @brief Returns the sound speed given density and internal energy
- *
- * Computes \f$c = \sqrt{\gamma (\gamma - 1) u }\f$.
- *
- * @param density The density \f$\rho\f$
- * @param u The internal energy \f$u\f$
- */
-__attribute__((always_inline)) INLINE static float
-gas_soundspeed_from_internal_energy(float density, float u) {
-
-  return sqrtf(u * hydro_gamma * hydro_gamma_minus_one);
-}
-
-/**
- * @brief Returns the sound speed given density and pressure
- *
- * Computes \f$c = \sqrt{\frac{\gamma P}{\rho} }\f$.
- *
- * @param density The density \f$\rho\f$
- * @param P The pressure \f$P\f$
- */
-__attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
-    float density, float P) {
-
-  const float density_inv = 1.f / density;
-  return sqrtf(hydro_gamma * P * density_inv);
-}
-
-/* ------------------------------------------------------------------------- */
+#include "./equation_of_state/ideal_gas/equation_of_state.h"
 #elif defined(EOS_ISOTHERMAL_GAS)
-
-struct eos_parameters {
-
-  /*! Thermal energy per unit mass */
-  float isothermal_internal_energy;
-};
-
-/**
- * @brief Returns the internal energy given density and entropy
- *
- * Since we are using an isothermal EoS, the entropy and density values are
- * ignored.
- * Computes \f$u = u_{cst}\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param entropy The entropy \f$S\f$.
- */
-__attribute__((always_inline)) INLINE static float
-gas_internal_energy_from_entropy(float density, float entropy) {
-
-  return eos.isothermal_internal_energy;
-}
-
-/**
- * @brief Returns the pressure given density and entropy
- *
- * Since we are using an isothermal EoS, the entropy value is ignored.
- * Computes \f$P = (\gamma - 1)u_{cst}\rho\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param entropy The entropy \f$S\f$.
- */
-__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
-    float density, float entropy) {
-
-  return hydro_gamma_minus_one * eos.isothermal_internal_energy * density;
-}
-
-/**
- * @brief Returns the entropy given density and pressure.
- *
- * Since we are using an isothermal EoS, the pressure value is ignored.
- * Computes \f$A = (\gamma - 1)u_{cst}\rho^{-(\gamma-1)}\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param pressure The pressure \f$P\f$ (ignored).
- * @return The entropy \f$A\f$.
- */
-__attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
-    float density, float pressure) {
-
-  return hydro_gamma_minus_one * eos.isothermal_internal_energy *
-         pow_minus_gamma_minus_one(density);
-}
-
-/**
- * @brief Returns the sound speed given density and entropy
- *
- * Since we are using an isothermal EoS, the entropy and density values are
- * ignored.
- * Computes \f$c = \sqrt{u_{cst} \gamma (\gamma-1)}\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param entropy The entropy \f$S\f$.
- */
-__attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
-    float density, float entropy) {
-
-  return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
-               hydro_gamma_minus_one);
-}
-
-/**
- * @brief Returns the entropy given density and internal energy
- *
- * Since we are using an isothermal EoS, the energy value is ignored.
- * Computes \f$S = \frac{(\gamma - 1)u_{cst}}{\rho^{\gamma-1}}\f$.
- *
- * @param density The density \f$\rho\f$
- * @param u The internal energy \f$u\f$
- */
-__attribute__((always_inline)) INLINE static float
-gas_entropy_from_internal_energy(float density, float u) {
-
-  return hydro_gamma_minus_one * eos.isothermal_internal_energy *
-         pow_minus_gamma_minus_one(density);
-}
-
-/**
- * @brief Returns the pressure given density and internal energy
- *
- * Since we are using an isothermal EoS, the energy value is ignored.
- * Computes \f$P = (\gamma - 1)u_{cst}\rho\f$.
- *
- * @param density The density \f$\rho\f$
- * @param u The internal energy \f$u\f$
- */
-__attribute__((always_inline)) INLINE static float
-gas_pressure_from_internal_energy(float density, float u) {
-
-  return hydro_gamma_minus_one * eos.isothermal_internal_energy * density;
-}
-
-/**
- * @brief Returns the internal energy given density and pressure.
- *
- * Just returns the constant internal energy.
- *
- * @param density The density \f$\rho\f$ (ignored).
- * @param pressure The pressure \f$P\f$ (ignored).
- * @return The internal energy \f$u\f$ (which is constant).
- */
-__attribute__((always_inline)) INLINE static float
-gas_internal_energy_from_pressure(float density, float pressure) {
-  return eos.isothermal_internal_energy;
-}
-
-/**
- * @brief Returns the sound speed given density and internal energy
- *
- * Since we are using an isothermal EoS, the energy and density values are
- * ignored.
- * Computes \f$c = \sqrt{u_{cst} \gamma (\gamma-1)}\f$.
- *
- * @param density The density \f$\rho\f$
- * @param u The internal energy \f$u\f$
- */
-__attribute__((always_inline)) INLINE static float
-gas_soundspeed_from_internal_energy(float density, float u) {
-
-  return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
-               hydro_gamma_minus_one);
-}
-
-/**
- * @brief Returns the sound speed given density and pressure
- *
- * Since we are using an isothermal EoS, the pressure and density values are
- * ignored.
- * Computes \f$c = \sqrt{u_{cst} \gamma (\gamma-1)}\f$.
- *
- * @param density The density \f$\rho\f$
- * @param P The pressure \f$P\f$
- */
-__attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
-    float density, float P) {
-
-  return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
-               hydro_gamma_minus_one);
-}
-
-/* ------------------------------------------------------------------------- */
+#include "./equation_of_state/isothermal/equation_of_state.h"
 #else
-
-#error "An Equation of state needs to be chosen in const.h !"
-
-#endif
-
-void eos_init(struct eos_parameters *e, const struct swift_params *params);
-void eos_print(const struct eos_parameters *e);
-
-#if defined(HAVE_HDF5)
-void eos_print_snapshot(hid_t h_grpsph, const struct eos_parameters *e);
+#error "Invalid choice of equation of state"
 #endif
 
 #endif /* SWIFT_EQUATION_OF_STATE_H */
diff --git a/src/equation_of_state/ideal_gas/equation_of_state.h b/src/equation_of_state/ideal_gas/equation_of_state.h
new file mode 100644
index 0000000000000000000000000000000000000000..42314e0fd87cb5ee2b81bc8cf29029ab4951c869
--- /dev/null
+++ b/src/equation_of_state/ideal_gas/equation_of_state.h
@@ -0,0 +1,203 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016   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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_IDEAL_GAS_EQUATION_OF_STATE_H
+#define SWIFT_IDEAL_GAS_EQUATION_OF_STATE_H
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "adiabatic_index.h"
+#include "common_io.h"
+#include "debug.h"
+#include "inline.h"
+
+extern struct eos_parameters eos;
+/* ------------------------------------------------------------------------- */
+
+struct eos_parameters {};
+
+/**
+ * @brief Returns the internal energy given density and entropy
+ *
+ * Computes \f$u = \frac{S\rho^{\gamma-1} }{\gamma - 1}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float
+gas_internal_energy_from_entropy(float density, float entropy) {
+
+  return entropy * pow_gamma_minus_one(density) *
+         hydro_one_over_gamma_minus_one;
+}
+
+/**
+ * @brief Returns the pressure given density and entropy
+ *
+ * Computes \f$P = S\rho^\gamma\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
+    float density, float entropy) {
+
+  return entropy * pow_gamma(density);
+}
+
+/**
+ * @brief Returns the entropy given density and pressure.
+ *
+ * Computes \f$A = \frac{P}{\rho^\gamma}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param pressure The pressure \f$P\f$.
+ * @return The entropy \f$A\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
+    float density, float pressure) {
+
+  return pressure * pow_minus_gamma(density);
+}
+
+/**
+ * @brief Returns the sound speed given density and entropy
+ *
+ * Computes \f$c = \sqrt{\gamma S \rho^{\gamma-1}}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
+    float density, float entropy) {
+
+  return sqrtf(hydro_gamma * pow_gamma_minus_one(density) * entropy);
+}
+
+/**
+ * @brief Returns the entropy given density and internal energy
+ *
+ * Computes \f$S = \frac{(\gamma - 1)u}{\rho^{\gamma-1}}\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_entropy_from_internal_energy(float density, float u) {
+
+  return hydro_gamma_minus_one * u * pow_minus_gamma_minus_one(density);
+}
+
+/**
+ * @brief Returns the pressure given density and internal energy
+ *
+ * Computes \f$P = (\gamma - 1)u\rho\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_pressure_from_internal_energy(float density, float u) {
+
+  return hydro_gamma_minus_one * u * density;
+}
+
+/**
+ * @brief Returns the internal energy given density and pressure.
+ *
+ * Computes \f$u = \frac{1}{\gamma - 1}\frac{P}{\rho}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param pressure The pressure \f$P\f$.
+ * @return The internal energy \f$u\f$.
+ */
+__attribute__((always_inline)) INLINE static float
+gas_internal_energy_from_pressure(float density, float pressure) {
+  return hydro_one_over_gamma_minus_one * pressure / density;
+}
+
+/**
+ * @brief Returns the sound speed given density and internal energy
+ *
+ * Computes \f$c = \sqrt{\gamma (\gamma - 1) u }\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_soundspeed_from_internal_energy(float density, float u) {
+
+  return sqrtf(u * hydro_gamma * hydro_gamma_minus_one);
+}
+
+/**
+ * @brief Returns the sound speed given density and pressure
+ *
+ * Computes \f$c = \sqrt{\frac{\gamma P}{\rho} }\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param P The pressure \f$P\f$
+ */
+__attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
+    float density, float P) {
+
+  const float density_inv = 1.f / density;
+  return sqrtf(hydro_gamma * P * density_inv);
+}
+
+/**
+ * @brief Initialize the eos parameters
+ *
+ * @param e The #eos_paramters
+ * @param params The parsed parameters
+ */
+__attribute__((always_inline)) INLINE static void eos_init(
+    struct eos_parameters *e, const struct swift_params *params) {}
+
+/**
+ * @brief Print the equation of state
+ *
+ * @param e The #eos_parameters
+ */
+__attribute__((always_inline)) INLINE static void eos_print(
+    const struct eos_parameters *e) {
+
+  message("Equation of state: Ideal gas.");
+
+  message("Adiabatic index gamma: %f.", hydro_gamma);
+}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Write equation of state information to the snapshot
+ *
+ * @param h_grpsph The HDF5 group in which to write
+ * @param e The #eos_parameters
+ */
+__attribute__((always_inline)) INLINE static void eos_print_snapshot(
+    hid_t h_grpsph, const struct eos_parameters *e) {
+
+  io_write_attribute_f(h_grpsph, "Adiabatic index", hydro_gamma);
+
+  io_write_attribute_s(h_grpsph, "Equation of state", "Ideal gas");
+}
+#endif
+
+#endif /* SWIFT_IDEAL_GAS_EQUATION_OF_STATE_H */
diff --git a/src/equation_of_state/isothermal/equation_of_state.h b/src/equation_of_state/isothermal/equation_of_state.h
new file mode 100644
index 0000000000000000000000000000000000000000..71890b4df656cb5f44d3cb0fbb3bd6005bd6ab6a
--- /dev/null
+++ b/src/equation_of_state/isothermal/equation_of_state.h
@@ -0,0 +1,233 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016   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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_ISOTHERMAL_EQUATION_OF_STATE_H
+#define SWIFT_ISOTHERMAL_EQUATION_OF_STATE_H
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "adiabatic_index.h"
+#include "common_io.h"
+#include "debug.h"
+#include "inline.h"
+
+extern struct eos_parameters eos;
+/* ------------------------------------------------------------------------- */
+
+struct eos_parameters {
+
+  /*! Thermal energy per unit mass */
+  float isothermal_internal_energy;
+};
+
+/**
+ * @brief Returns the internal energy given density and entropy
+ *
+ * Since we are using an isothermal EoS, the entropy and density values are
+ * ignored.
+ * Computes \f$u = u_{cst}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float
+gas_internal_energy_from_entropy(float density, float entropy) {
+
+  return eos.isothermal_internal_energy;
+}
+
+/**
+ * @brief Returns the pressure given density and entropy
+ *
+ * Since we are using an isothermal EoS, the entropy value is ignored.
+ * Computes \f$P = (\gamma - 1)u_{cst}\rho\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
+    float density, float entropy) {
+
+  return hydro_gamma_minus_one * eos.isothermal_internal_energy * density;
+}
+
+/**
+ * @brief Returns the entropy given density and pressure.
+ *
+ * Since we are using an isothermal EoS, the pressure value is ignored.
+ * Computes \f$A = (\gamma - 1)u_{cst}\rho^{-(\gamma-1)}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param pressure The pressure \f$P\f$ (ignored).
+ * @return The entropy \f$A\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
+    float density, float pressure) {
+
+  return hydro_gamma_minus_one * eos.isothermal_internal_energy *
+         pow_minus_gamma_minus_one(density);
+}
+
+/**
+ * @brief Returns the sound speed given density and entropy
+ *
+ * Since we are using an isothermal EoS, the entropy and density values are
+ * ignored.
+ * Computes \f$c = \sqrt{u_{cst} \gamma (\gamma-1)}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
+    float density, float entropy) {
+
+  return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
+               hydro_gamma_minus_one);
+}
+
+/**
+ * @brief Returns the entropy given density and internal energy
+ *
+ * Since we are using an isothermal EoS, the energy value is ignored.
+ * Computes \f$S = \frac{(\gamma - 1)u_{cst}}{\rho^{\gamma-1}}\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_entropy_from_internal_energy(float density, float u) {
+
+  return hydro_gamma_minus_one * eos.isothermal_internal_energy *
+         pow_minus_gamma_minus_one(density);
+}
+
+/**
+ * @brief Returns the pressure given density and internal energy
+ *
+ * Since we are using an isothermal EoS, the energy value is ignored.
+ * Computes \f$P = (\gamma - 1)u_{cst}\rho\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_pressure_from_internal_energy(float density, float u) {
+
+  return hydro_gamma_minus_one * eos.isothermal_internal_energy * density;
+}
+
+/**
+ * @brief Returns the internal energy given density and pressure.
+ *
+ * Just returns the constant internal energy.
+ *
+ * @param density The density \f$\rho\f$ (ignored).
+ * @param pressure The pressure \f$P\f$ (ignored).
+ * @return The internal energy \f$u\f$ (which is constant).
+ */
+__attribute__((always_inline)) INLINE static float
+gas_internal_energy_from_pressure(float density, float pressure) {
+  return eos.isothermal_internal_energy;
+}
+
+/**
+ * @brief Returns the sound speed given density and internal energy
+ *
+ * Since we are using an isothermal EoS, the energy and density values are
+ * ignored.
+ * Computes \f$c = \sqrt{u_{cst} \gamma (\gamma-1)}\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_soundspeed_from_internal_energy(float density, float u) {
+
+  return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
+               hydro_gamma_minus_one);
+}
+
+/**
+ * @brief Returns the sound speed given density and pressure
+ *
+ * Since we are using an isothermal EoS, the pressure and density values are
+ * ignored.
+ * Computes \f$c = \sqrt{u_{cst} \gamma (\gamma-1)}\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param P The pressure \f$P\f$
+ */
+__attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
+    float density, float P) {
+
+  return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
+               hydro_gamma_minus_one);
+}
+
+/* ------------------------------------------------------------------------- */
+
+/**
+ * @brief Initialize the eos parameters
+ *
+ * @param e The #eos_paramters
+ * @param params The parsed parameters
+ */
+__attribute__((always_inline)) INLINE static void eos_init(
+    struct eos_parameters *e, const struct swift_params *params) {
+
+  e->isothermal_internal_energy =
+      parser_get_param_float(params, "EoS:isothermal_internal_energy");
+}
+
+/**
+ * @brief Print the equation of state
+ *
+ * @param e The #eos_parameters
+ */
+__attribute__((always_inline)) INLINE static void eos_print(
+    const struct eos_parameters *e) {
+
+  message(
+      "Equation of state: Isothermal with internal energy "
+      "per unit mass set to %f.",
+      e->isothermal_internal_energy);
+
+  message("Adiabatic index gamma: %f.", hydro_gamma);
+}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Write equation of state information to the snapshot
+ *
+ * @param h_grpsph The HDF5 group in which to write
+ * @param e The #eos_parameters
+ */
+__attribute__((always_inline)) INLINE static void eos_print_snapshot(
+    hid_t h_grpsph, const struct eos_parameters *e) {
+
+  io_write_attribute_f(h_grpsph, "Adiabatic index", hydro_gamma);
+
+  io_write_attribute_s(h_grpsph, "Equation of state", "Isothermal gas");
+  io_write_attribute_f(h_grpsph, "Thermal energy per unit mass",
+                       e->isothermal_internal_energy);
+}
+#endif
+
+#endif /* SWIFT_ISOTHERMAL_EQUATION_OF_STATE_H */
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index b7fa7dd4aa00249d009037961888c89abc0d0313..b2af8909bed1780586a5130370222c9b8157d724 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -59,9 +59,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   const float mi = pi->mass;
   const float mj = pj->mass;
 
-  /* Get r and r inverse. */
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and 1/r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Compute the kernel function for pi */
   const float hi_inv = 1.f / hi;
@@ -148,9 +148,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   /* Get the masses. */
   const float mj = pj->mass;
 
-  /* Get r and r inverse. */
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and 1/r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Compute the kernel function */
   const float hi_inv = 1.0f / hi;
@@ -440,8 +440,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float fac_mu = pow_three_gamma_minus_five_over_two(a);
   const float a2_Hubble = a * a * H;
 
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and 1/r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Get some values in local variables. */
   const float mi = pi->mass;
@@ -559,8 +560,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float fac_mu = pow_three_gamma_minus_five_over_two(a);
   const float a2_Hubble = a * a * H;
 
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and 1/r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Get some values in local variables. */
   // const float mi = pi->mass;
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index bfd8f1dff5febbd14b8a1e3a1cb8c01d49bc85ba..42fd93d6062cbfcea5cf5297eeda0bb6525f3cad 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -53,7 +53,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
 
   float wi, wj, wi_dx, wj_dx;
 
-  const float r = sqrtf(r2);
+  /* Get r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Get the masses. */
   const float mi = pi->mass;
@@ -101,8 +103,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   /* Get the masses. */
   const float mj = pj->mass;
 
-  /* Get r and r inverse. */
-  const float r = sqrtf(r2);
+  /* Get r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   const float h_inv = 1.f / hi;
   const float ui = r * h_inv;
@@ -134,8 +137,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float fac_mu = pow_three_gamma_minus_five_over_two(a);
   const float a2_Hubble = a * a * H;
 
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and r inverse. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Recover some data */
   const float mi = pi->mass;
@@ -246,8 +250,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float fac_mu = pow_three_gamma_minus_five_over_two(a);
   const float a2_Hubble = a * a * H;
 
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and r inverse. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Recover some data */
   // const float mi = pi->mass;
diff --git a/src/hydro/PressureEntropy/hydro_iact.h b/src/hydro/PressureEntropy/hydro_iact.h
index c9e3c4ec7728df8a8e7d59b8c21f22bd613463d8..b8f8c1983a3b1fb67781f7228194deb770273988 100644
--- a/src/hydro/PressureEntropy/hydro_iact.h
+++ b/src/hydro/PressureEntropy/hydro_iact.h
@@ -54,9 +54,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   const float mi = pi->mass;
   const float mj = pj->mass;
 
-  /* Get r and r inverse. */
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and 1/r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Compute the kernel function for pi */
   const float hi_inv = 1.f / hi;
@@ -142,9 +142,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   /* Get the masses. */
   const float mj = pj->mass;
 
-  /* Get r and r inverse. */
-  const float r = sqrtf(r2);
-  const float ri = 1.0f / r;
+  /* Get r and 1/r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Compute the kernel function */
   const float h_inv = 1.0f / hi;
@@ -164,7 +164,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   pi->density.pressure_dh -=
       mj * pj->entropy_one_over_gamma * (hydro_dimension * wi + ui * wi_dx);
 
-  const float fac = mj * wi_dx * ri;
+  const float fac = mj * wi_dx * r_inv;
 
   /* Compute dv dot r */
   dv[0] = pi->v[0] - pj->v[0];
@@ -205,8 +205,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float fac_mu = pow_three_gamma_minus_five_over_two(a);
   const float a2_Hubble = a * a * H;
 
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and 1/r . */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Get some values in local variables. */
   const float mi = pi->mass;
@@ -318,8 +319,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float fac_mu = pow_three_gamma_minus_five_over_two(a);
   const float a2_Hubble = a * a * H;
 
-  const float r = sqrtf(r2);
-  const float r_inv = 1.0f / r;
+  /* Get r and 1/r. */
+  const float r_inv = 1.0f / sqrtf(r2);
+  const float r = r2 * r_inv;
 
   /* Get some values in local variables. */
   // const float mi = pi->mass;
diff --git a/tests/tolerance_125_normal.dat b/tests/tolerance_125_normal.dat
index 0f11d03507b23c76b5703e118eede1359fe2afba..ff3d42e22b5b91378114294c4c81405492eafad2 100644
--- a/tests/tolerance_125_normal.dat
+++ b/tests/tolerance_125_normal.dat
@@ -1,4 +1,4 @@
 #   ID    pos_x    pos_y    pos_z      v_x      v_y      v_z        h      rho    div_v        S        u        P        c      a_x      a_y      a_z     h_dt    v_sig    dS/dt    du/dt
     0	  1e-4	   1e-4	    1e-4       1e-4	1e-4	 1e-4	    1e-4   1e-4	  1e-4	       1e-4	1e-4	 1e-4	  1e-4	 1e-4	  1e-4	   1e-4	   1e-4	   1e-4	    1e-4     1e-4
     0	  1e-4	   1e-4	    1e-4       1e-4	1e-4	 1e-4	    1e-4   1e-4	  1e-4	       1e-4	1e-4	 1e-4	  1e-4	 1e-4	  1e-4	   1e-4	   1e-4	   1e-4	    1e-4     1e-4
-    0	  1e-6	   1e-6	    1e-6       1e-6	1e-6	 1e-6	    1e-6   1e-6	  1e-6	       1e-6	1e-6	 1e-6	  1e-6	 1e-5	  1e-5	   1e-5	   1e-5	   1e-5	    1e-5     1e-5
+    0	  1e-6	   1e-6	    1e-6       1e-6	1e-6	 1e-6	    1e-6   1e-6	  1e-6	       1e-6	1e-6	 1e-6	  1e-6	 2e-5	  2e-5	   2e-5	   1e-5	   1e-5	    1e-5     1e-5
diff --git a/tests/tolerance_27_normal.dat b/tests/tolerance_27_normal.dat
index 0fe55e84a42e7541068744e1e554afff1731ed3f..713e6c20c39ae9be0a4069a17c1ce45728fb5a34 100644
--- a/tests/tolerance_27_normal.dat
+++ b/tests/tolerance_27_normal.dat
@@ -1,4 +1,4 @@
 #   ID      pos_x      pos_y      pos_z        v_x        v_y        v_z           rho        rho_dh        wcount     wcount_dh         div_v       curl_vx       curl_vy       curl_vz
-    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   2e-6	      4e-5	    4e-4       1e-2		 1e-5	     6e-6	   6e-6		 6e-6
-    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1e-6	      1.2e-4	    1e-4       2e-4		 2e-4	     1e-4	   1e-4	 	 1e-4
+    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   2e-6	      6.2e-5	    4e-4       1.2e-2		 1e-5	     6e-6	   6e-6		 8e-6
+    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1e-6	      2e-4	    1e-4       2e-4		 2e-4	     1e-4	   1e-4	 	 1e-4
     0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1e-6	      1e-6	    1e-6       1e-6		 1e-6	     1e-6	   1e-6		 1e-6
diff --git a/tests/tolerance_27_perturbed.dat b/tests/tolerance_27_perturbed.dat
index e1fd7518e55811386967938dfff0d62f04325cb0..4b7534f4ba43e7117452af819937f2d2b0451352 100644
--- a/tests/tolerance_27_perturbed.dat
+++ b/tests/tolerance_27_perturbed.dat
@@ -1,4 +1,4 @@
 #   ID      pos_x      pos_y      pos_z        v_x        v_y        v_z           rho        rho_dh        wcount     wcount_dh         div_v       curl_vx       curl_vy       curl_vz
     0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   2e-6       1e-4	    2e-4       1e-2		 1e-5	     3e-6	   3e-6		 7e-6
-    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1e-6	      1.3e-3	    1e-5       2e-3		 6e-5	     2e-3	   2e-3	 	 2e-3
+    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1e-6	      1.3e-3	    1e-5       2e-3		 6e-5	     3e-3	   2e-3	 	 2e-3
     0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1e-6	      2e-3	    1e-6       1e0		 1e-6	     2e-6	   2e-6		 2e-6
diff --git a/tests/tolerance_27_perturbed_h.dat b/tests/tolerance_27_perturbed_h.dat
index 5d6710c3946c0515f94c58eb3e0ecbfe14135b71..a1703f5b296d16f045fd28462fc9c61fc9358585 100644
--- a/tests/tolerance_27_perturbed_h.dat
+++ b/tests/tolerance_27_perturbed_h.dat
@@ -1,4 +1,4 @@
 #   ID      pos_x      pos_y      pos_z        v_x        v_y        v_z           rho        rho_dh        wcount     wcount_dh         div_v       curl_vx       curl_vy       curl_vz
-    0	      1e-6       1e-6	      1e-6 	       1e-6 	    1e-6	     1e-6	         3e-6       1e-4	        5e-4       1.4e-2		         1.1e-5	       3e-6	         3e-6		       8e-6
-    0	      1e-6       1e-6	      1e-6 	       1e-6 	    1e-6	     1e-6	         1.5e-6	      1.4e-2	      1e-5     2e-3		           2.5e-4	     3e-3	         3e-3	 	       3e-3
-    0	      1e-6       1e-6	      1e-6 	       1e-6 	    1e-6	     1e-6	         1e-6	      1e-6	        1e-6       1e0		           1e-6	       4e-6	         4e-6		       4e-6
+    0	      1e-6       1e-6	      1e-6       1e-6 	    1e-6     1e-6         3e-6        1e-4	     5e-4       1.4e-2	         1.1e-5	       3e-6	         3e-6		       8e-6
+    0	      1e-6       1e-6	      1e-6       1e-6 	    1e-6     1e-6         1.5e-6      1.4e-2	     1e-5       2e-3	         2.5e-4	       3e-3	         3e-3	 	       3e-3
+    0	      1e-6       1e-6	      1e-6       1e-6 	    1e-6     1e-6         1e-6	      1e-6	     1e-6       1e0	         1e-6	       4e-6	         4e-6		       4e-6
diff --git a/tests/tolerance_27_perturbed_h2.dat b/tests/tolerance_27_perturbed_h2.dat
index 882554741734aeb270427b62580d6077907056ad..0b3921f0294f03c199e2045364144756dc5b066e 100644
--- a/tests/tolerance_27_perturbed_h2.dat
+++ b/tests/tolerance_27_perturbed_h2.dat
@@ -1,4 +1,4 @@
-#   ID      pos_x      pos_y      pos_z        v_x        v_y        v_z           rho        rho_dh        wcount     wcount_dh         div_v       curl_vx       curl_vy       curl_vz
-    0	      1e-6       1e-6	      1e-6 	       1e-6 	    1e-6	     1e-6	         3e-6       1e-4	        5e-4       1.5e-2		         1.4e-5	       3e-6	         3e-6		       1e-5
-    0	      1e-6       1e-6	      1e-6 	       1e-6 	    1e-6	     1e-6	         1.5e-6	    2.5e-2	      1e-5       5.86e-3		       1.17e-3	     3e-3	         5.65e-3	 	       3e-3
-    0	      1e-6       1e-6	      1e-6 	       1e-6 	    1e-6	     1e-6	         1e-6	      1e-6	        1e-6       1e0		           1e-6	         4e-6	         4e-6		       4e-6
+#   ID        pos_x      pos_y      pos_z        v_x        v_y         v_z        rho        rho_dh    wcount     wcount_dh    div_v       curl_vx       curl_vy       curl_vz
+    0	      1e-6       1e-6	    1e-6 	 1e-6 	    1e-6	1e-6	   3e-6       1e-4      5e-4       1.5e-2	1.4e-5	     3e-6	   3e-6		 1e-5
+    0	      1e-6       1e-6	    1e-6 	 1e-6 	    1e-6	1e-6	 1.5e-6	    2.5e-2      1e-5       5.86e-3	1.17e-3	     3e-3	   8e-3	         3e-3
+    0	      1e-6       1e-6	    1e-6	 1e-6 	    1e-6	1e-6	   1e-6	      1e-6      1e-6       1e0		1e-6	     4e-6	   4e-6		 4e-6
diff --git a/tests/tolerance_pair_force_active.dat b/tests/tolerance_pair_force_active.dat
index c5b5620338b848f1b8f4658049e53a6424e93a60..f4394edf35e7287968673dfb949c726462b5d8ad 100644
--- a/tests/tolerance_pair_force_active.dat
+++ b/tests/tolerance_pair_force_active.dat
@@ -1,4 +1,4 @@
 #   ID  wcount   h_dt
     0	   1      2.4e-3
     0	   1      2.4e-3
-    0	   1      1e-5
+    0	   1      3e-5
diff --git a/tests/tolerance_periodic_BC_normal.dat b/tests/tolerance_periodic_BC_normal.dat
index 823e4af488b343f57e3c90e89ee2d4f13d3ca94b..fdbcd53a1f0233a87d657e5bc3751e685631a4c7 100644
--- a/tests/tolerance_periodic_BC_normal.dat
+++ b/tests/tolerance_periodic_BC_normal.dat
@@ -1,4 +1,4 @@
 #   ID      pos_x      pos_y      pos_z        v_x        v_y        v_z           rho        rho_dh        wcount     wcount_dh         div_v       curl_vx       curl_vy       curl_vz
-    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   4e-6	      4e-5	    1e-3       1e-2		 2e-4	     2e-4	   2e-4		 2e-4
+    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   4e-6	      7e-5	    1e-3       1.3e-2		 2e-4	     2e-4	   2e-4		 2e-4
     0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   2e-6	      2e-4	    1e-4       2e-4		 6e-4	     2e-3	   2e-3	 	 2e-3
     0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1e-6	      1e-4	    1e-6       1e-4		 5e-4	     2e-4	   2e-4	 	 2e-4