diff --git a/src/const.h b/src/const.h
index 737a6870fcac3777c954cfbf68e2e5ca1b16ba0b..f41034101be3738f7e956c4becf047a50f007476 100644
--- a/src/const.h
+++ b/src/const.h
@@ -45,7 +45,7 @@
 /* Try to keep cells regular by adding a correction velocity. */
 //#define GIZMO_STEER_MOTION
 /* Use the total energy instead of the thermal energy as conserved variable. */
-#define GIZMO_TOTAL_ENERGY
+//#define GIZMO_TOTAL_ENERGY
 
 /* Options to control handling of unphysical values (GIZMO_SPH only). */
 /* In GIZMO, mass and energy (and hence density and pressure) can in principle
diff --git a/src/hydro/GizmoMFM/hydro.h b/src/hydro/GizmoMFM/hydro.h
index d57fd61ed9883ea6797a4da792edd5566369785d..dad33486edfb3c06063627db01f6f33993e65e58 100644
--- a/src/hydro/GizmoMFM/hydro.h
+++ b/src/hydro/GizmoMFM/hydro.h
@@ -156,6 +156,10 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
      is to have a way of remembering that we need more neighbours for this
      particle */
   p->geometry.wcorr = 1.0f;
+
+  /* set the entropy to a non-zero initial value (to make sure we
+     don't zero out the density before we reach hydro_convert_quantities) */
+  p->conserved.entropy = 1.0f;
 }
 
 /**
@@ -334,6 +338,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
   /* energy contains the total thermal energy, we want the specific energy.
      this is why we divide by the volume, and not by the density */
   p->P = hydro_gamma_minus_one * energy * volume_inv;
+  p->A = p->conserved.entropy / m;
+  p->P = p->A * pow_gamma(p->rho);
 #endif
 
   /* sanity checks */
@@ -518,6 +524,11 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
     const struct hydro_props* hydro_props) {
 
   p->conserved.energy /= cosmo->a_factor_internal_energy;
+
+  /* initialize the entropy */
+  p->conserved.entropy =
+    p->conserved.energy * pow_minus_gamma_minus_one(p->rho) *
+      hydro_gamma_minus_one;
 }
 
 /**
@@ -570,6 +581,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
     const float u = (p->conserved.energy + p->flux.energy * dt_drift) * m_inv;
 #endif
     p->P = hydro_gamma_minus_one * u * p->rho;
+    p->P = p->A * pow_gamma(p->rho);
 #endif
   }
 
@@ -635,12 +647,14 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
   p->conserved.momentum[0] += p->flux.momentum[0] * dt_therm;
   p->conserved.momentum[1] += p->flux.momentum[1] * dt_therm;
   p->conserved.momentum[2] += p->flux.momentum[2] * dt_therm;
+
 #if defined(EOS_ISOTHERMAL_GAS)
   /* We use the EoS equation in a sneaky way here just to get the constant u */
   p->conserved.energy =
       p->conserved.mass * gas_internal_energy_from_entropy(0.0f, 0.0f);
 #else
   p->conserved.energy += p->flux.energy * dt_therm;
+  p->conserved.energy = hydro_one_over_gamma_minus_one * p->conserved.entropy * pow_gamma_minus_one(p->rho);
 #endif
 
 #ifndef HYDRO_GAMMA_5_3
diff --git a/src/hydro/GizmoMFM/hydro_gradients.h b/src/hydro/GizmoMFM/hydro_gradients.h
index 6f751d970287ca7ba137c1138ca6f3bf00e6c4cb..805cebf68d2ef039b9e4dd156fe953f046bc4801 100644
--- a/src/hydro/GizmoMFM/hydro_gradients.h
+++ b/src/hydro/GizmoMFM/hydro_gradients.h
@@ -95,7 +95,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
     struct part* restrict pi, struct part* restrict pj, float hi, float hj,
     const float* dx, float r, const float* xij_i, float* Wi, float* Wj) {
 
-  /* perform gradient reconstruction in space and time */
+  /* perform gradient reconstruction in space */
   /* Compute interface position (relative to pj, since we don't need the actual
    * position) eqn. (8) */
   const float xij_j[3] = {xij_i[0] + dx[0], xij_i[1] + dx[1], xij_i[2] + dx[2]};
diff --git a/src/hydro/GizmoMFM/hydro_part.h b/src/hydro/GizmoMFM/hydro_part.h
index 8097b1b2560f24f78636bbb855700054524fe0bb..0adfed995252e0c0bf8223bf25fd032ded8d377a 100644
--- a/src/hydro/GizmoMFM/hydro_part.h
+++ b/src/hydro/GizmoMFM/hydro_part.h
@@ -77,6 +77,9 @@ struct part {
   /* Pressure. */
   float P;
 
+  /* Entropic function. */
+  float A;
+
   union {
     /* Quantities used during the volume (=density) loop. */
     struct {
@@ -169,6 +172,9 @@ struct part {
     /* Fluid thermal energy (not per unit mass!). */
     float energy;
 
+    /* Fluid entropy. */
+    float entropy;
+
   } conserved;
 
   /* Geometrical quantities used for hydro. */