diff --git a/examples/CoolingHaloWithSpin/cooling_halo.yml b/examples/CoolingHaloWithSpin/cooling_halo.yml
index 3c2266f6fc28a2989c2b2a1909adddfb8e41abf6..2e5f09b9cae199090f29ae2341c232218f8099d0 100644
--- a/examples/CoolingHaloWithSpin/cooling_halo.yml
+++ b/examples/CoolingHaloWithSpin/cooling_halo.yml
@@ -10,7 +10,7 @@ InternalUnitSystem:
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
   time_end:   10.    # The end time of the simulation (in internal units).
-  dt_min:     1e-7  # The minimal time-step size of the simulation (in internal units).
+  dt_min:     1e-10 # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-1  # The maximal time-step size of the simulation (in internal units).
 
 # Parameters governing the conserved quantities statistics
@@ -35,12 +35,12 @@ InitialConditions:
  
 # External potential parameters
 IsothermalPotential:
-  position_x:      0.     # location of centre of isothermal potential in internal units
+  position_x:      0.     # Location of centre of isothermal potential in internal units
   position_y:      0.
   position_z:      0.	
-  vrot:            200.     # rotation speed of isothermal potential in internal units
-  timestep_mult:   0.03     # controls time step
-  epsilon:         1.0      #softening for the isothermal potential
+  vrot:            200.   # Rotation speed of isothermal potential in internal units
+  timestep_mult:   0.03   # Controls time step
+  epsilon:         1.0    # Softening for the isothermal potential
 
 # Cooling parameters
 LambdaCooling:
@@ -48,4 +48,4 @@ LambdaCooling:
   minimum_temperature:         1.0e4  # Minimal temperature (Kelvin)
   mean_molecular_weight:       0.59   # Mean molecular weight
   hydrogen_mass_abundance:     0.75   # Hydrogen mass abundance (dimensionless)
-  cooling_tstep_mult:          0.1    # Dimensionless pre-factor for the time-step condition
+  cooling_tstep_mult:          1.0    # Dimensionless pre-factor for the time-step condition
diff --git a/examples/CoolingHaloWithSpin/run.sh b/examples/CoolingHaloWithSpin/run.sh
index 3a0d9c02000e760b030a96107038d3c6163f3227..131fbf3cb10d2014546683b5f43194840544fd55 100755
--- a/examples/CoolingHaloWithSpin/run.sh
+++ b/examples/CoolingHaloWithSpin/run.sh
@@ -4,7 +4,8 @@
 echo "Generating initial conditions for the isothermal potential box example..."
 python makeIC.py 10000 
 
-../swift -g -s -C -t 16 cooling_halo.yml 2>&1 | tee output.log
+# Run SWIFT with external potential, SPH and cooling
+../swift -g -s -C -t 1 cooling_halo.yml 2>&1 | tee output.log
 
 # python radial_profile.py 10
 
diff --git a/src/debug.c b/src/debug.c
index 21f539b62eef3b0b36f52520486f1e725abc0cda..f5f2f4974a6f2d0e8da8fce71e98233a2ed3deeb 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -194,7 +194,8 @@ int checkSpacehmax(struct space *s) {
 
 /**
  * @brief Check if the h_max and dx_max values of a cell's hierarchy are
- * consistent with the particles. Report verbosely if not.
+ * consistent with the particles. Also checks if particles are correctly
+ * in a cell. Report verbosely if not.
  *
  * @param c the top cell of the hierarchy.
  * @param depth the recursion depth for use in messages. Set to 0 initially.
@@ -206,43 +207,62 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
 
   float h_max = 0.0f;
   float dx_max = 0.0f;
-  if (!c->split) {
-    const size_t nr_parts = c->count;
-    struct part *parts = c->parts;
-    for (size_t k = 0; k < nr_parts; k++) {
-      h_max = (h_max > parts[k].h) ? h_max : parts[k].h;
+  int result = 1;
+
+  const double loc_min[3] = {c->loc[0], c->loc[1], c->loc[2]};
+  const double loc_max[3] = {c->loc[0] + c->width[0], c->loc[1] + c->width[1],
+                             c->loc[2] + c->width[2]};
+
+  const size_t nr_parts = c->count;
+  struct part *parts = c->parts;
+  struct xpart *xparts = c->xparts;
+  for (size_t k = 0; k < nr_parts; k++) {
+
+    struct part *const p = &parts[k];
+    struct xpart *const xp = &xparts[k];
+
+    if (p->x[0] < loc_min[0] || p->x[0] > loc_max[0] || p->x[1] < loc_min[1] ||
+        p->x[1] > loc_max[1] || p->x[2] < loc_min[2] || p->x[2] > loc_max[2]) {
+
+      message(
+          "Inconsistent part position p->x=[%e %e %e], c->loc=[%e %e %e] "
+          "c->width=[%e %e %e]",
+          p->x[0], p->x[1], p->x[2], c->loc[0], c->loc[1], c->loc[2],
+          c->width[0], c->width[1], c->width[2]);
+
+      result = 0;
     }
-  } else {
-    for (int k = 0; k < 8; k++)
+
+    const float dx2 = xp->x_diff[0] * xp->x_diff[0] +
+                      xp->x_diff[1] * xp->x_diff[1] +
+                      xp->x_diff[2] * xp->x_diff[2];
+
+    h_max = max(h_max, p->h);
+    dx_max = max(dx_max, sqrt(dx2));
+  }
+
+  if (c->split) {
+    for (int k = 0; k < 8; k++) {
       if (c->progeny[k] != NULL) {
         struct cell *cp = c->progeny[k];
         checkCellhdxmax(cp, depth);
-        dx_max = max(dx_max, cp->dx_max);
-        h_max = max(h_max, cp->h_max);
       }
+    }
   }
 
   /* Check. */
-  int result = 1;
   if (c->h_max != h_max) {
     message("%d Inconsistent h_max: cell %f != parts %f", *depth, c->h_max,
             h_max);
-    error("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
+    message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
     result = 0;
   }
   if (c->dx_max != dx_max) {
     message("%d Inconsistent dx_max: %f != %f", *depth, c->dx_max, dx_max);
-    error("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
+    message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
     result = 0;
   }
 
-  /* Check rebuild criterion. */
-  /* if (h_max > c->dmin) {
-    message("%d Inconsistent c->dmin: %f > %f", *depth, h_max, c->dmin);
-    error("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
-    result = 0;
-  } */
-
   return result;
 }
 
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index 3fd357a2d8778f5ca8b014935d538350eccb99c6..bfb5cd1ce39a9908573c66406f41b56561a870d6 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -33,7 +33,7 @@
  * @param dt Time since the last kick
  */
 __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
-    const struct part *restrict p, float dt) {
+    const struct part *restrict p) {
 
   return p->u;
 }
@@ -45,7 +45,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
  * @param dt Time since the last kick
  */
 __attribute__((always_inline)) INLINE static float hydro_get_pressure(
-    const struct part *restrict p, float dt) {
+    const struct part *restrict p) {
 
   return gas_pressure_from_internal_energy(p->rho, p->u);
 }
@@ -57,7 +57,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_pressure(
  * @param dt Time since the last kick
  */
 __attribute__((always_inline)) INLINE static float hydro_get_entropy(
-    const struct part *restrict p, float dt) {
+    const struct part *restrict p) {
 
   return gas_entropy_from_internal_energy(p->rho, p->u);
 }
@@ -69,7 +69,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
  * @param dt Time since the last kick
  */
 __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
-    const struct part *restrict p, float dt) {
+    const struct part *restrict p) {
 
   return p->force.soundspeed;
 }
@@ -97,34 +97,30 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
 }
 
 /**
- * @brief Modifies the thermal state of a particle to the imposed internal
- * energy
+ * @brief Returns the time derivative of internal energy of a particle
  *
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * We assume a constant density.
  *
- * @param p The particle
- * @param u The new internal energy
+ * @param p The particle of interest
  */
-__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
-    struct part *restrict p, float u) {
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
+    const struct part *restrict p) {
 
-  p->u = u;
+  return p->force.u_dt;
 }
 
 /**
- * @brief Modifies the thermal state of a particle to the imposed entropy
+ * @brief Returns the time derivative of internal energy of a particle
  *
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * We assume a constant density.
  *
- * @param p The particle
- * @param S The new entropy
+ * @param p The particle of interest.
+ * @param du_dt The new time derivative of the internal energy.
  */
-__attribute__((always_inline)) INLINE static void hydro_set_entropy(
-    struct part *restrict p, float S) {
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
+    struct part *restrict p, float du_dt) {
 
-  p->u = gas_internal_energy_from_entropy(p->rho, S);
+  p->force.u_dt = du_dt;
 }
 
 /**
@@ -152,26 +148,6 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
   return min(dt_cfl, dt_u_change);
 }
 
-/**
- * @brief Initialises the particles for the first time
- *
- * This function is called only once just after the ICs have been
- * read in to do some conversions.
- *
- * @param p The particle to act upon
- * @param xp The extended particle data to act upon
- */
-__attribute__((always_inline)) INLINE static void hydro_first_init_part(
-    struct part *restrict p, struct xpart *restrict xp) {
-
-  p->ti_begin = 0;
-  p->ti_end = 0;
-  xp->v_full[0] = p->v[0];
-  xp->v_full[1] = p->v[1];
-  xp->v_full[2] = p->v[2];
-  xp->u_full = p->u;
-}
-
 /**
  * @brief Prepares a particle for the density calculation.
  *
@@ -244,8 +220,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
  * @param time The current time
  */
 __attribute__((always_inline)) INLINE static void hydro_prepare_force(
-    struct part *restrict p, struct xpart *restrict xp, int ti_current,
-    double timeBase) {
+    struct part *restrict p, struct xpart *restrict xp) {
 
   /* Some smoothing length multiples. */
   const float h = p->h;
@@ -270,17 +245,18 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   p->force.balsara = normDiv_v / (normDiv_v + normRot_v + 0.0001f * fc * h_inv);
 
   /* Viscosity parameter decay time */
-  const float tau = h / (2.f * const_viscosity_length * p->force.soundspeed);
+  /* const float tau = h / (2.f * const_viscosity_length * p->force.soundspeed);
+   */
 
   /* Viscosity source term */
-  const float S = max(-normDiv_v, 0.f);
+  /* const float S = max(-normDiv_v, 0.f); */
 
   /* Compute the particle's viscosity parameter time derivative */
-  const float alpha_dot = (const_viscosity_alpha_min - p->alpha) / tau +
-                          (const_viscosity_alpha_max - p->alpha) * S;
+  /* const float alpha_dot = (const_viscosity_alpha_min - p->alpha) / tau + */
+  /*                         (const_viscosity_alpha_max - p->alpha) * S; */
 
   /* Update particle's viscosity paramter */
-  p->alpha += alpha_dot * (p->ti_end - p->ti_begin) * timeBase;
+  /* p->alpha += alpha_dot * (p->ti_end - p->ti_begin) * timeBase; */  // MATTHIEU
 }
 
 /**
@@ -305,6 +281,22 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
   p->force.v_sig = 0.0f;
 }
 
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param p The particle.
+ * @param xp The extended data of this particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
+    struct part *restrict p, const struct xpart *restrict xp) {
+
+  /* Re-set the predicted velocities */
+  p->v[0] = xp->v_full[0];
+  p->v[1] = xp->v_full[1];
+  p->v[2] = xp->v_full[2];
+}
+
 /**
  * @brief Predict additional particle fields forward in time when drifting
  *
@@ -316,8 +308,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
  * @param timeBase The minimal time-step size
  */
 __attribute__((always_inline)) INLINE static void hydro_predict_extra(
-    struct part *restrict p, struct xpart *restrict xp, float dt, int t0,
-    int t1, double timeBase) {
+    struct part *restrict p, struct xpart *restrict xp, float dt) {
   float u, w;
 
   const float h_inv = 1.f / p->h;
@@ -368,8 +359,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
  * @param half_dt The half time-step for this kick
  */
 __attribute__((always_inline)) INLINE static void hydro_kick_extra(
-    struct part *restrict p, struct xpart *restrict xp, float dt,
-    float half_dt) {}
+    struct part *restrict p, struct xpart *restrict xp, float dt) {}
 
 /**
  *  @brief Converts hydro quantity of a particle at the start of a run
@@ -379,6 +369,28 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
  * @param p The particle to act upon
  */
 __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
-    struct part *restrict p) {}
+    struct part *restrict p, struct xpart *restrict xp) {}
+
+/**
+ * @brief Initialises the particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_first_init_part(
+    struct part *restrict p, struct xpart *restrict xp) {
+
+  p->time_bin = 0;
+  xp->v_full[0] = p->v[0];
+  xp->v_full[1] = p->v[1];
+  xp->v_full[2] = p->v[2];
+  xp->u_full = p->u;
+
+  hydro_reset_acceleration(p);
+  hydro_init_part(p);
+}
 
 #endif /* SWIFT_DEFAULT_HYDRO_H */
diff --git a/src/hydro/Default/hydro_debug.h b/src/hydro/Default/hydro_debug.h
index d02d3ef82c1b3d751731f49850c06df4b146b164..3be9c9e1760591423edbd218d19b46ddf9aad01e 100644
--- a/src/hydro/Default/hydro_debug.h
+++ b/src/hydro/Default/hydro_debug.h
@@ -25,11 +25,10 @@ __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],\n "
       "h=%.3e, "
-      "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, t_begin=%d, t_end=%d\n",
+      "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, time_bin=%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->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->ti_begin,
-      p->ti_end);
+      p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->time_bin);
 }
 
 #endif /* SWIFT_DEFAULT_HYDRO_DEBUG_H */
diff --git a/src/hydro/Default/hydro_part.h b/src/hydro/Default/hydro_part.h
index c7464bcf338b1c5b81ffa91d92264c2bd35e9313..332eecb27fb65a6b4da48cbb595450a432c44615 100644
--- a/src/hydro/Default/hydro_part.h
+++ b/src/hydro/Default/hydro_part.h
@@ -55,12 +55,6 @@ struct part {
   /* Particle cutoff radius. */
   float h;
 
-  /* Particle time of beginning of time-step. */
-  int ti_begin;
-
-  /* Particle time of end of time-step. */
-  int ti_end;
-
   /* Particle internal energy. */
   float u;
 
@@ -125,6 +119,9 @@ struct part {
   /* Pointer to corresponding gravity part. */
   struct gpart* gpart;
 
+  /* Particle time-bin */
+  timebin_t time_bin;
+
 } SWIFT_STRUCT_ALIGN;
 
 #endif /* SWIFT_DEFAULT_HYDRO_PART_H */
diff --git a/src/runner.c b/src/runner.c
index 3b17fb8a953e873b05a1a59c72bf54483427ca28..803ef5e3b06f453de2161cfb2a9f82c3daf85926 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -734,8 +734,16 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
       }
     }
 
+#ifdef SWIFT_DEBUG_CHECKS
+    if (count) {
+      message("Smoothing length failed to converge on %i particles.", count);
+
+      error("Aborting....");
+    }
+#else
     if (count)
       message("Smoothing length failed to converge on %i particles.", count);
+#endif
 
     /* Be clean */
     free(pid);
@@ -1319,7 +1327,13 @@ void *runner_main(void *data) {
 #ifdef SWIFT_DEBUG_CHECKS
       t->ti_run = e->ti_current;
 #ifndef WITH_MPI
-      if (cj == NULL) { /* self */
+      if (ci == NULL && cj == NULL) {
+
+        if (t->type != task_type_grav_gather_m && t->type != task_type_grav_fft)
+          error("Task not associated with cells!");
+
+      } else if (cj == NULL) { /* self */
+
         if (!cell_is_active(ci, e) && t->type != task_type_sort &&
             t->type != task_type_send && t->type != task_type_recv)
           error(
diff --git a/src/scheduler.c b/src/scheduler.c
index 2c6b84b7d12dee3b440ba575ca321a685aa823a1..f9e0533b2567e3a1dfc1d46eb262cd729c6768da 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -1056,8 +1056,6 @@ void scheduler_enqueue_mapper(void *map_data, int num_elements,
  */
 void scheduler_start(struct scheduler *s) {
 
-  message("launching %i active tasks.", s->active_count);
-
   /* Re-wait the tasks. */
   if (s->active_count > 1000) {
     threadpool_map(s->threadpool, scheduler_rewait_mapper, s->tid_active,
@@ -1082,7 +1080,12 @@ void scheduler_start(struct scheduler *s) {
       /* Don't check MPI stuff */
       if (t->type == task_type_send || t->type == task_type_recv) continue;
 
-      if (cj == NULL) { /* self */
+     if (ci == NULL && cj == NULL) {	
+
+        if (t->type != task_type_grav_gather_m && t->type != task_type_grav_fft)
+          error("Task not associated with cells!");
+
+      } else if (cj == NULL) { /* self */
 
         if (ci->ti_end_min == ti_current && t->skip &&
             t->type != task_type_sort && t->type)
diff --git a/src/timestep.h b/src/timestep.h
index 63a21bf3e2c82344f9a0b1ac5f4ac769af10c7d1..b3ef96493772e7e93853ad28a0536e4f8448dceb 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -70,14 +70,15 @@ make_integer_timestep(float new_dt, timebin_t old_bin, integertime_t ti_current,
 __attribute__((always_inline)) INLINE static integertime_t get_gpart_timestep(
     const struct gpart *restrict gp, const struct engine *restrict e) {
 
-  const float new_dt_external = external_gravity_timestep(
-      e->time, e->external_potential, e->physical_constants, gp);
+  float new_dt = FLT_MAX;
 
-  /* const float new_dt_self = */
-  /*     gravity_compute_timestep_self(e->physical_constants, gp); */
-  const float new_dt_self = FLT_MAX;  // MATTHIEU
+  if (e->policy & engine_policy_external_gravity)
+    new_dt =
+        min(new_dt, external_gravity_timestep(e->time, e->external_potential,
+                                              e->physical_constants, gp));
 
-  float new_dt = min(new_dt_external, new_dt_self);
+  if (e->policy & engine_policy_self_gravity)
+    new_dt = min(new_dt, gravity_compute_timestep_self(gp));
 
   /* Limit timestep within the allowed range */
   new_dt = min(new_dt, e->dt_max);
@@ -114,14 +115,13 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
   float new_dt_grav = FLT_MAX;
   if (p->gpart != NULL) {
 
-    const float new_dt_external = external_gravity_timestep(
-        e->time, e->external_potential, e->physical_constants, p->gpart);
+    if (e->policy & engine_policy_external_gravity)
+      new_dt_grav = min(new_dt_grav, external_gravity_timestep(
+                                         e->time, e->external_potential,
+                                         e->physical_constants, p->gpart));
 
-    /* const float new_dt_self = */
-    /*     gravity_compute_timestep_self(e->physical_constants, p->gpart); */
-    const float new_dt_self = FLT_MAX;  // MATTHIEU
-
-    new_dt_grav = min(new_dt_external, new_dt_self);
+    if (e->policy & engine_policy_self_gravity)
+      new_dt_grav = min(new_dt_grav, gravity_compute_timestep_self(p->gpart));
   }
 
   /* Final time-step is minimum of hydro and gravity */