diff --git a/INSTALL.swift b/INSTALL.swift
index 8e1635b0715426512503fd9dcde32f59a7ad1b62..90096bc1b88d34ab86c04f623f1d3f04ca5a2997 100644
--- a/INSTALL.swift
+++ b/INSTALL.swift
@@ -143,6 +143,8 @@ before you can build it.
  - DOXYGEN: the doxygen library is required to create the SWIFT API
             documentation.
 
+ - python:  Examples and solution script use python and rely on the
+   	    numpy library version 1.8.2 or higher.
 
 
                              SWIFT Coding style
diff --git a/examples/BigCosmoVolume/makeIC.py b/examples/BigCosmoVolume/makeIC.py
index c141337c06fb28aa4049e2823fcc7cd3e9d5513c..8f83b564a3f747d164fd03b7dddf3cd9c609d9eb 100644
--- a/examples/BigCosmoVolume/makeIC.py
+++ b/examples/BigCosmoVolume/makeIC.py
@@ -107,7 +107,7 @@ if n_copy > 1:
     for i in range(n_copy):
         for j in range(n_copy):
             for k in range(n_copy):
-                coords = np.append(coords, coords_tile + np.array([ i * boxSize, j * boxSize, k * boxSize ]), axis=0)
+                coords = np.append(coords, coords_tile + np.array([ i * boxSize[0], j * boxSize[1], k * boxSize[2] ]), axis=0)
                 v = np.append(v, v_tile, axis=0)
                 m = np.append(m, m_tile)
                 h = np.append(h, h_tile)
diff --git a/examples/CoolingHalo/makeIC.py b/examples/CoolingHalo/makeIC.py
index 0b542e200da709e2cc7f668ab8b62b94e0bf95ee..3ec1be6f7b5e568ebe8e0fefe508ef8287edb29c 100644
--- a/examples/CoolingHalo/makeIC.py
+++ b/examples/CoolingHalo/makeIC.py
@@ -227,7 +227,7 @@ ds[()] = u
 u = np.zeros(1)
 
 # Particle IDs
-ids = 1 + np.linspace(0, N, N, endpoint=False, dtype='L')
+ids = 1 + np.linspace(0, N, N, endpoint=False)
 ds = grp.create_dataset('ParticleIDs', (N, ), 'L')
 ds[()] = ids
 
diff --git a/examples/CoolingHaloWithSpin/makeIC.py b/examples/CoolingHaloWithSpin/makeIC.py
index a6d57868ad7542498b27007a5c3ef9234b9feb84..2cf3127c743f61756b3ff6c4a7738c83d185f9cd 100644
--- a/examples/CoolingHaloWithSpin/makeIC.py
+++ b/examples/CoolingHaloWithSpin/makeIC.py
@@ -233,7 +233,7 @@ ds[()] = u
 u = np.zeros(1)
 
 # Particle IDs
-ids = 1 + np.linspace(0, N, N, endpoint=False, dtype='L')
+ids = 1 + np.linspace(0, N, N, endpoint=False)
 ds = grp.create_dataset('ParticleIDs', (N, ), 'L')
 ds[()] = ids
 
diff --git a/examples/DiscPatch/GravityOnly/makeIC.py b/examples/DiscPatch/GravityOnly/makeIC.py
index 42cd26e235deb17a899a65050ef5caa9c832c59c..6e4e148392eb7ca2fbf8c29c3f737d029916c59b 100644
--- a/examples/DiscPatch/GravityOnly/makeIC.py
+++ b/examples/DiscPatch/GravityOnly/makeIC.py
@@ -150,7 +150,7 @@ ds[()] = m
 m = numpy.zeros(1)
 
 
-ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False, dtype='L')
+ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False)
 ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L')
 ds[()] = ids
 
diff --git a/examples/DiscPatch/HydroStatic/makeIC.py b/examples/DiscPatch/HydroStatic/makeIC.py
index e2846d82a8cfa8bf08de83632b19ae2e7818f3c1..6ba1ccd06fed84ca728aadaa5922dbba536b6881 100644
--- a/examples/DiscPatch/HydroStatic/makeIC.py
+++ b/examples/DiscPatch/HydroStatic/makeIC.py
@@ -205,7 +205,7 @@ if (entropy_flag == 1):
 else:
     ds[()] = u    
 
-ids = 1 + numpy.linspace(0, numGas, numGas, endpoint=False, dtype='L')
+ids = 1 + numpy.linspace(0, numGas, numGas, endpoint=False)
 ds = grp0.create_dataset('ParticleIDs', (numGas, ), 'L')
 ds[()] = ids
 
diff --git a/examples/HydrostaticHalo/density_profile.py b/examples/HydrostaticHalo/density_profile.py
index 52bebb9ffefa77dae66af155fb31fed539dcde13..d0afd399f951cf3b727e869ca8571a3a802c2e8d 100644
--- a/examples/HydrostaticHalo/density_profile.py
+++ b/examples/HydrostaticHalo/density_profile.py
@@ -28,7 +28,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
 unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
+v_c = float(params.attrs["IsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
 #lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
 #X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
diff --git a/examples/HydrostaticHalo/internal_energy_profile.py b/examples/HydrostaticHalo/internal_energy_profile.py
index a1b2bda314a66eb965974d34519f66c544ee8aed..ea52cf8fc5fd098a46f05eaa58494529a868000c 100644
--- a/examples/HydrostaticHalo/internal_energy_profile.py
+++ b/examples/HydrostaticHalo/internal_energy_profile.py
@@ -46,7 +46,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
 unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
+v_c = float(params.attrs["IsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
 #lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
 #X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
diff --git a/examples/HydrostaticHalo/makeIC.py b/examples/HydrostaticHalo/makeIC.py
index f33387e18dd0ab523684227f5c745b5c8b807b7f..d5081ac84473edc87857c6872278b4d0ca6389b1 100644
--- a/examples/HydrostaticHalo/makeIC.py
+++ b/examples/HydrostaticHalo/makeIC.py
@@ -227,7 +227,7 @@ ds[()] = u
 u = np.zeros(1)
 
 # Particle IDs
-ids = 1 + np.linspace(0, N, N, endpoint=False, dtype='L')
+ids = 1 + np.linspace(0, N, N, endpoint=False)
 ds = grp.create_dataset('ParticleIDs', (N, ), 'L')
 ds[()] = ids
 
diff --git a/examples/HydrostaticHalo/velocity_profile.py b/examples/HydrostaticHalo/velocity_profile.py
index f6a7350b9731d660b2092266d4d6ad3730bab48c..9133195d942233514148aa419003ee0ab7923494 100644
--- a/examples/HydrostaticHalo/velocity_profile.py
+++ b/examples/HydrostaticHalo/velocity_profile.py
@@ -46,7 +46,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
 unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
+v_c = float(params.attrs["IsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
 header = f["Header"]
 N = header.attrs["NumPart_Total"][0]
diff --git a/examples/IsothermalPotential/makeIC.py b/examples/IsothermalPotential/makeIC.py
index 7d1c5361f9a255365517226e49c55a8a50c4d6ce..27ddf15fe63d69d4172b927cfca8b8662d11ea4e 100644
--- a/examples/IsothermalPotential/makeIC.py
+++ b/examples/IsothermalPotential/makeIC.py
@@ -138,7 +138,7 @@ ds = grp1.create_dataset('Masses', (numPart,), 'f')
 ds[()] = m
 m = numpy.zeros(1)
 
-ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False, dtype='L')
+ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False)
 ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L')
 ds[()] = ids
 
diff --git a/examples/Noh_3D/makeIC.py b/examples/Noh_3D/makeIC.py
index ec8d46639eecdefd55b917a955f13efc8ffe4126..0c25a5c8b3e967185cf16bae4b1f21c215266def 100644
--- a/examples/Noh_3D/makeIC.py
+++ b/examples/Noh_3D/makeIC.py
@@ -35,8 +35,8 @@ glass = h5py.File("glassCube_64.hdf5", "r")
 
 vol = 8.
 
-pos = glass["/PartType0/Coordinates"][:,:] * cbrt(vol)
-h = glass["/PartType0/SmoothingLength"][:] * cbrt(vol)
+pos = glass["/PartType0/Coordinates"][:,:] * vol**(1./3.)
+h = glass["/PartType0/SmoothingLength"][:] * vol**(1./3.)
 numPart = size(h)
 
 # Generate extra arrays
@@ -65,7 +65,7 @@ file = h5py.File(fileName, 'w')
 
 # Header
 grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = [cbrt(vol), cbrt(vol), cbrt(vol)]
+grp.attrs["BoxSize"] = [vol**(1./3.), vol**(1./3.), vol**(1./3.)]
 grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
diff --git a/examples/SineWavePotential_3D/plotSolution.py b/examples/SineWavePotential_3D/plotSolution.py
index 7ae5dcd2a52235812c3c11cfcf4d01f4fde647af..13cae037b64eff4ad4fec0003bf0f5ad3ce94896 100644
--- a/examples/SineWavePotential_3D/plotSolution.py
+++ b/examples/SineWavePotential_3D/plotSolution.py
@@ -55,7 +55,7 @@ rho_x = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*x))
 
 P = cs2*rho
 
-n1D = int(np.cbrt(len(P)))
+n1D = np.ceil(len(P)**(1./3.))
 gradP = np.zeros(P.shape)
 for i in range(len(P)):
   iself = int(ids[i]/n1D/n1D)
diff --git a/examples/main.c b/examples/main.c
index 9ab40413b048cf08e035aeec6775bdc8ab1ab9ca..b687ab18263a79bef9cc79360e135531d7ee7f9a 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -609,7 +609,7 @@ int main(int argc, char *argv[]) {
            clocks_getunit());
 
   /* Main simulation loop */
-  for (int j = 0; !engine_is_done(&e) && e.step != nsteps; j++) {
+  for (int j = 0; !engine_is_done(&e) && e.step - 1 != nsteps; j++) {
 
     /* Reset timers */
     timers_reset(timers_mask_all);
diff --git a/src/active.h b/src/active.h
index 38020fde8ed53a638231097643476b7ef72d6d49..02e504f762735994e6c57f7e155071fede016713 100644
--- a/src/active.h
+++ b/src/active.h
@@ -105,10 +105,12 @@ __attribute__((always_inline)) INLINE static int cell_is_all_active(
 __attribute__((always_inline)) INLINE static int part_is_active(
     const struct part *p, const struct engine *e) {
 
-  const integertime_t ti_current = e->ti_current;
-  const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+  const timebin_t max_active_bin = e->max_active_bin;
+  const timebin_t part_bin = p->time_bin;
 
 #ifdef SWIFT_DEBUG_CHECKS
+  const integertime_t ti_current = e->ti_current;
+  const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
   if (ti_end < ti_current)
     error(
         "particle in an impossible time-zone! p->ti_end=%lld "
@@ -116,7 +118,7 @@ __attribute__((always_inline)) INLINE static int part_is_active(
         ti_end, ti_current);
 #endif
 
-  return (ti_end == ti_current);
+  return (part_bin <= max_active_bin);
 }
 
 /**
@@ -129,10 +131,13 @@ __attribute__((always_inline)) INLINE static int part_is_active(
 __attribute__((always_inline)) INLINE static int gpart_is_active(
     const struct gpart *gp, const struct engine *e) {
 
+  const timebin_t max_active_bin = e->max_active_bin;
+  const timebin_t gpart_bin = gp->time_bin;
+
+#ifdef SWIFT_DEBUG_CHECKS
   const integertime_t ti_current = e->ti_current;
   const integertime_t ti_end = get_integer_time_end(ti_current, gp->time_bin);
 
-#ifdef SWIFT_DEBUG_CHECKS
   if (ti_end < ti_current)
     error(
         "g-particle in an impossible time-zone! gp->ti_end=%lld "
@@ -140,7 +145,7 @@ __attribute__((always_inline)) INLINE static int gpart_is_active(
         ti_end, ti_current);
 #endif
 
-  return (ti_end == ti_current);
+  return (gpart_bin <= max_active_bin);
 }
 
 /**
@@ -153,10 +158,13 @@ __attribute__((always_inline)) INLINE static int gpart_is_active(
 __attribute__((always_inline)) INLINE static int spart_is_active(
     const struct spart *sp, const struct engine *e) {
 
+  const timebin_t max_active_bin = e->max_active_bin;
+  const timebin_t spart_bin = sp->time_bin;
+
+#ifdef SWIFT_DEBUG_CHECKS
   const integertime_t ti_current = e->ti_current;
   const integertime_t ti_end = get_integer_time_end(ti_current, sp->time_bin);
 
-#ifdef SWIFT_DEBUG_CHECKS
   if (ti_end < ti_current)
     error(
         "s-particle in an impossible time-zone! sp->ti_end=%lld "
@@ -164,7 +172,7 @@ __attribute__((always_inline)) INLINE static int spart_is_active(
         ti_end, ti_current);
 #endif
 
-  return (ti_end == ti_current);
+  return (spart_bin <= max_active_bin);
 }
 
 /* Are cells / particles active for kick1 tasks ? */
@@ -201,11 +209,14 @@ __attribute__((always_inline)) INLINE static int cell_is_starting(
 __attribute__((always_inline)) INLINE static int part_is_starting(
     const struct part *p, const struct engine *e) {
 
+  const timebin_t max_active_bin = e->max_active_bin;
+  const timebin_t part_bin = p->time_bin;
+
+#ifdef SWIFT_DEBUG_CHECKS
   const integertime_t ti_current = e->ti_current;
   const integertime_t ti_beg =
       get_integer_time_begin(ti_current + 1, p->time_bin);
 
-#ifdef SWIFT_DEBUG_CHECKS
   if (ti_beg > ti_current)
     error(
         "particle in an impossible time-zone! p->ti_beg=%lld "
@@ -213,7 +224,7 @@ __attribute__((always_inline)) INLINE static int part_is_starting(
         ti_beg, ti_current);
 #endif
 
-  return (ti_beg == ti_current);
+  return (part_bin <= max_active_bin);
 }
 
 /**
@@ -226,11 +237,14 @@ __attribute__((always_inline)) INLINE static int part_is_starting(
 __attribute__((always_inline)) INLINE static int gpart_is_starting(
     const struct gpart *gp, const struct engine *e) {
 
+  const timebin_t max_active_bin = e->max_active_bin;
+  const timebin_t gpart_bin = gp->time_bin;
+
+#ifdef SWIFT_DEBUG_CHECKS
   const integertime_t ti_current = e->ti_current;
   const integertime_t ti_beg =
       get_integer_time_begin(ti_current + 1, gp->time_bin);
 
-#ifdef SWIFT_DEBUG_CHECKS
   if (ti_beg > ti_current)
     error(
         "g-particle in an impossible time-zone! gp->ti_beg=%lld "
@@ -238,7 +252,7 @@ __attribute__((always_inline)) INLINE static int gpart_is_starting(
         ti_beg, ti_current);
 #endif
 
-  return (ti_beg == ti_current);
+  return (gpart_bin <= max_active_bin);
 }
 
 /**
@@ -251,11 +265,14 @@ __attribute__((always_inline)) INLINE static int gpart_is_starting(
 __attribute__((always_inline)) INLINE static int spart_is_starting(
     const struct spart *sp, const struct engine *e) {
 
+  const timebin_t max_active_bin = e->max_active_bin;
+  const timebin_t spart_bin = sp->time_bin;
+
+#ifdef SWIFT_DEBUG_CHECKS
   const integertime_t ti_current = e->ti_current;
   const integertime_t ti_beg =
       get_integer_time_begin(ti_current + 1, sp->time_bin);
 
-#ifdef SWIFT_DEBUG_CHECKS
   if (ti_beg > ti_current)
     error(
         "s-particle in an impossible time-zone! sp->ti_beg=%lld "
@@ -263,6 +280,6 @@ __attribute__((always_inline)) INLINE static int spart_is_starting(
         ti_beg, ti_current);
 #endif
 
-  return (ti_beg == ti_current);
+  return (spart_bin <= max_active_bin);
 }
 #endif /* SWIFT_ACTIVE_H */
diff --git a/src/engine.c b/src/engine.c
index facffa389e98fdeeebda222f04b2953c00b97373..52d4be9b3470e398efdd9669191c92790d472444 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2804,6 +2804,18 @@ void engine_print_stats(struct engine *e) {
 
   const ticks tic = getticks();
 
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Check that all cells have been drifted to the current time.
+   * That can include cells that have not
+   * previously been active on this rank. */
+  space_check_drift_point(e->s, e->ti_current);
+
+  /* Be verbose about this */
+  message("Saving statistics at t=%e.", e->time);
+#else
+  if (e->verbose) message("Saving statistics at t=%e.", e->time);
+#endif
+
   e->save_stats = 0;
 
   struct statistics stats;
@@ -3010,7 +3022,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
 #endif
 
   /* Ready to go */
-  e->step = -1;
+  e->step = 0;
   e->forcerebuild = 1;
   e->wallclock_time = (float)clocks_diff(&time1, &time2);
 
@@ -3031,14 +3043,6 @@ void engine_step(struct engine *e) {
 
   e->tic_step = getticks();
 
-  /* Move forward in time */
-  e->ti_old = e->ti_current;
-  e->ti_current = e->ti_end_min;
-  e->step += 1;
-  e->time = e->ti_current * e->timeBase + e->timeBegin;
-  e->timeOld = e->ti_old * e->timeBase + e->timeBegin;
-  e->timeStep = (e->ti_current - e->ti_old) * e->timeBase;
-
   if (e->nodeID == 0) {
 
     /* Print some information to the screen */
@@ -3053,6 +3057,15 @@ void engine_step(struct engine *e) {
     fflush(e->file_timesteps);
   }
 
+  /* Move forward in time */
+  e->ti_old = e->ti_current;
+  e->ti_current = e->ti_end_min;
+  e->max_active_bin = get_max_active_bin(e->ti_end_min);
+  e->step += 1;
+  e->time = e->ti_current * e->timeBase + e->timeBegin;
+  e->timeOld = e->ti_old * e->timeBase + e->timeBegin;
+  e->timeStep = (e->ti_current - e->ti_old) * e->timeBase;
+
   /* Prepare the tasks to be launched, rebuild or repartition if needed. */
   engine_prepare(e);
 
@@ -3566,6 +3579,7 @@ void engine_init(struct engine *e, struct space *s,
   e->time = e->timeBegin;
   e->ti_old = 0;
   e->ti_current = 0;
+  e->max_active_bin = num_time_bins;
   e->timeStep = 0.;
   e->timeBase = 0.;
   e->timeBase_inv = 0.;
diff --git a/src/engine.h b/src/engine.h
index 35b0d4a58225e4c825eb1002ad405cb7840b8bd6..80903ee78156c7e4efb2650e59e4fa832fecbfa3 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -122,6 +122,9 @@ struct engine {
   double time;
   integertime_t ti_current;
 
+  /* The highest active bin at this time */
+  timebin_t max_active_bin;
+
   /* Time step */
   double timeStep;
 
diff --git a/src/runner.c b/src/runner.c
index 6448f21d41e307a33a0eaae7f2f486f1fdcd3fd9..792f03f9b454d2afe0116f9bda194bb00ce79543 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1418,6 +1418,8 @@ void runner_do_recv_part(struct runner *r, struct cell *c, int timer) {
 
   integertime_t ti_end_min = max_nr_timesteps;
   integertime_t ti_end_max = 0;
+  timebin_t time_bin_min = num_time_bins;
+  timebin_t time_bin_max = 0;
   float h_max = 0.f;
 
   /* If this cell is a leaf, collect the particle data. */
@@ -1425,10 +1427,9 @@ void runner_do_recv_part(struct runner *r, struct cell *c, int timer) {
 
     /* Collect everything... */
     for (size_t k = 0; k < nr_parts; k++) {
-      const integertime_t ti_end =
-          get_integer_time_end(ti_current, parts[k].time_bin);
-      ti_end_min = min(ti_end_min, ti_end);
-      ti_end_max = max(ti_end_max, ti_end);
+      if (parts[k].time_bin == time_bin_inhibited) continue;
+      time_bin_min = min(time_bin_min, parts[k].time_bin);
+      time_bin_max = max(time_bin_max, parts[k].time_bin);
       h_max = max(h_max, parts[k].h);
 
 #ifdef SWIFT_DEBUG_CHECKS
@@ -1436,6 +1437,10 @@ void runner_do_recv_part(struct runner *r, struct cell *c, int timer) {
         error("Received un-drifted particle !");
 #endif
     }
+
+    /* Convert into a time */
+    ti_end_min = get_integer_time_end(ti_current, time_bin_min);
+    ti_end_max = get_integer_time_end(ti_current, time_bin_max);
   }
 
   /* Otherwise, recurse and collect. */
@@ -1490,22 +1495,27 @@ void runner_do_recv_gpart(struct runner *r, struct cell *c, int timer) {
 
   integertime_t ti_end_min = max_nr_timesteps;
   integertime_t ti_end_max = 0;
+  timebin_t time_bin_min = num_time_bins;
+  timebin_t time_bin_max = 0;
 
   /* If this cell is a leaf, collect the particle data. */
   if (!c->split) {
 
     /* Collect everything... */
     for (size_t k = 0; k < nr_gparts; k++) {
-      const integertime_t ti_end =
-          get_integer_time_end(ti_current, gparts[k].time_bin);
-      ti_end_min = min(ti_end_min, ti_end);
-      ti_end_max = max(ti_end_max, ti_end);
+      if (gparts[k].time_bin == time_bin_inhibited) continue;
+      time_bin_min = min(time_bin_min, gparts[k].time_bin);
+      time_bin_max = max(time_bin_max, gparts[k].time_bin);
 
 #ifdef SWIFT_DEBUG_CHECKS
       if (gparts[k].ti_drift != ti_current)
         error("Received un-drifted g-particle !");
 #endif
     }
+
+    /* Convert into a time */
+    ti_end_min = get_integer_time_end(ti_current, time_bin_min);
+    ti_end_max = get_integer_time_end(ti_current, time_bin_max);
   }
 
   /* Otherwise, recurse and collect. */
@@ -1558,17 +1568,27 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int timer) {
 
   integertime_t ti_end_min = max_nr_timesteps;
   integertime_t ti_end_max = 0;
+  timebin_t time_bin_min = num_time_bins;
+  timebin_t time_bin_max = 0;
 
   /* If this cell is a leaf, collect the particle data. */
   if (!c->split) {
 
     /* Collect everything... */
     for (size_t k = 0; k < nr_sparts; k++) {
-      const integertime_t ti_end =
-          get_integer_time_end(ti_current, sparts[k].time_bin);
-      ti_end_min = min(ti_end_min, ti_end);
-      ti_end_max = max(ti_end_max, ti_end);
+      if (sparts[k].time_bin == time_bin_inhibited) continue;
+      time_bin_min = min(time_bin_min, sparts[k].time_bin);
+      time_bin_max = max(time_bin_max, sparts[k].time_bin);
+
+#ifdef SWIFT_DEBUG_CHECKS
+      if (sparts[k].ti_drift != ti_current)
+        error("Received un-drifted s-particle !");
+#endif
     }
+
+    /* Convert into a time */
+    ti_end_min = get_integer_time_end(ti_current, time_bin_min);
+    ti_end_max = get_integer_time_end(ti_current, time_bin_max);
   }
 
   /* Otherwise, recurse and collect. */
diff --git a/src/space.c b/src/space.c
index 625fe944c488c871d02b14c72d15051b3e53b3c4..00e5b87ad6d3587b02d14d514c9788f8639c5dd5 100644
--- a/src/space.c
+++ b/src/space.c
@@ -2132,47 +2132,52 @@ void space_split_recursive(struct space *s, struct cell *c,
     c->split = 0;
     maxdepth = c->depth;
 
-    /* Get dt_min/dt_max. */
+    timebin_t time_bin_min = num_time_bins, time_bin_max = 0;
+
+    /* parts: Get dt_min/dt_max and h_max. */
+    for (int k = 0; k < count; k++) {
+#ifdef SWIFT_DEBUG_CHECKS
+      if (parts[k].time_bin == time_bin_inhibited)
+        error("Inhibited particle present in space_split()");
+#endif
+      time_bin_min = min(time_bin_min, parts[k].time_bin);
+      time_bin_max = max(time_bin_max, parts[k].time_bin);
+      h_max = max(h_max, parts[k].h);
+    }
+    /* parts: Reset x_diff */
     for (int k = 0; k < count; k++) {
-      struct part *p = &parts[k];
-      struct xpart *xp = &xparts[k];
-      const float h = p->h;
-      const integertime_t ti_end =
-          get_integer_time_end(e->ti_current, p->time_bin);
-      const integertime_t ti_beg =
-          get_integer_time_begin(e->ti_current + 1, p->time_bin);
-      xp->x_diff[0] = 0.f;
-      xp->x_diff[1] = 0.f;
-      xp->x_diff[2] = 0.f;
-      if (h > h_max) h_max = h;
-      if (ti_end < ti_end_min) ti_end_min = ti_end;
-      if (ti_end > ti_end_max) ti_end_max = ti_end;
-      if (ti_beg > ti_beg_max) ti_beg_max = ti_beg;
+      xparts[k].x_diff[0] = 0.f;
+      xparts[k].x_diff[1] = 0.f;
+      xparts[k].x_diff[2] = 0.f;
     }
+    /* gparts: Get dt_min/dt_max, reset x_diff. */
     for (int k = 0; k < gcount; k++) {
-      struct gpart *gp = &gparts[k];
-      const integertime_t ti_end =
-          get_integer_time_end(e->ti_current, gp->time_bin);
-      const integertime_t ti_beg =
-          get_integer_time_begin(e->ti_current + 1, gp->time_bin);
-      gp->x_diff[0] = 0.f;
-      gp->x_diff[1] = 0.f;
-      gp->x_diff[2] = 0.f;
-      if (ti_end < ti_end_min) ti_end_min = ti_end;
-      if (ti_end > ti_end_max) ti_end_max = ti_end;
-      if (ti_beg > ti_beg_max) ti_beg_max = ti_beg;
+#ifdef SWIFT_DEBUG_CHECKS
+      if (gparts[k].time_bin == time_bin_inhibited)
+        error("Inhibited s-particle present in space_split()");
+#endif
+      time_bin_min = min(time_bin_min, gparts[k].time_bin);
+      time_bin_max = max(time_bin_max, gparts[k].time_bin);
+
+      gparts[k].x_diff[0] = 0.f;
+      gparts[k].x_diff[1] = 0.f;
+      gparts[k].x_diff[2] = 0.f;
     }
+    /* sparts: Get dt_min/dt_max */
     for (int k = 0; k < scount; k++) {
-      struct spart *sp = &sparts[k];
-      const integertime_t ti_end =
-          get_integer_time_end(e->ti_current, sp->time_bin);
-      const integertime_t ti_beg =
-          get_integer_time_begin(e->ti_current + 1, sp->time_bin);
-      if (ti_end < ti_end_min) ti_end_min = ti_end;
-      if (ti_end > ti_end_max) ti_end_max = ti_end;
-      if (ti_beg > ti_beg_max) ti_beg_max = ti_beg;
+#ifdef SWIFT_DEBUG_CHECKS
+      if (sparts[k].time_bin == time_bin_inhibited)
+        error("Inhibited g-particle present in space_split()");
+#endif
+      time_bin_min = min(time_bin_min, sparts[k].time_bin);
+      time_bin_max = max(time_bin_max, sparts[k].time_bin);
     }
 
+    /* Convert into integer times */
+    ti_end_min = get_integer_time_end(e->ti_current, time_bin_min);
+    ti_end_max = get_integer_time_end(e->ti_current, time_bin_max);
+    ti_beg_max = get_integer_time_begin(e->ti_current + 1, time_bin_max);
+
     /* Construct the multipole and the centre of mass*/
     if (s->gravity) gravity_P2M(c->multipole, c->gparts, c->gcount);
   }
diff --git a/src/task.c b/src/task.c
index ff6bb9bee21de6a77d8a81a8f4f79722e483ad30..dfc3dc538c75297cbe7e79b7d64c1b7e13a015dc 100644
--- a/src/task.c
+++ b/src/task.c
@@ -283,6 +283,10 @@ void task_unlock(struct task *t) {
   /* Act based on task type. */
   switch (type) {
 
+    case task_type_init:
+    case task_type_kick1:
+    case task_type_kick2:
+    case task_type_timestep:
     case task_type_drift:
       cell_unlocktree(ci);
       cell_gunlocktree(ci);
@@ -355,6 +359,10 @@ int task_lock(struct task *t) {
 #endif
       break;
 
+    case task_type_init:
+    case task_type_kick1:
+    case task_type_kick2:
+    case task_type_timestep:
     case task_type_drift:
       if (ci->hold || ci->ghold) return 0;
       if (cell_locktree(ci) != 0) return 0;
diff --git a/src/timeline.h b/src/timeline.h
index 25ab231e24b85f2a4b6b3d44f024588420e432ea..0f38ff3e9421c122b61caf74290c170b62b2dd92 100644
--- a/src/timeline.h
+++ b/src/timeline.h
@@ -37,7 +37,8 @@ typedef char timebin_t;
 /*! The maximal number of timesteps in a simulation */
 #define max_nr_timesteps (1LL << (num_time_bins + 1))
 
-#define time_bin_inactive (num_time_bins + 2)
+/*! Fictious time-bin to hold inhibited particles */
+#define time_bin_inhibited (num_time_bins + 2)
 
 /**
  * @brief Returns the integer time interval corresponding to a time bin
diff --git a/tests/test125cells.c b/tests/test125cells.c
index aa2b577d606e97b0c2466f420870fa2855fd722e..21fc3f5407f90d9b75485d08bf1077e0a20e88b2 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -537,6 +537,7 @@ int main(int argc, char *argv[]) {
   engine.s = &space;
   engine.time = 0.1f;
   engine.ti_current = 8;
+  engine.max_active_bin = num_time_bins;
 
   struct runner runner;
   runner.e = &engine;
diff --git a/tests/test27cells.c b/tests/test27cells.c
index bece9b5e7ac5af547c8f14392ab544e371666869..5ef78c9509c545279063234992085b79808f628e 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -397,6 +397,7 @@ int main(int argc, char *argv[]) {
   engine.s = &space;
   engine.time = 0.1f;
   engine.ti_current = 8;
+  engine.max_active_bin = num_time_bins;
 
   struct runner runner;
   runner.e = &engine;
diff --git a/tests/testPair.c b/tests/testPair.c
index 6a0a53f7d4a34f4b055080f9c26a394a199b89b4..6f6e36f2d67bc4b60f9317cd5072e4affd292f9a 100644
--- a/tests/testPair.c
+++ b/tests/testPair.c
@@ -252,6 +252,7 @@ int main(int argc, char *argv[]) {
   engine.s = &space;
   engine.time = 0.1f;
   engine.ti_current = 8;
+  engine.max_active_bin = num_time_bins;
   runner.e = &engine;
 
   volume = particles * particles * particles;