diff --git a/README b/README
index cd2a397a18e872e7914b24fd58cc588ec1d6c8c0..2fc16c28f5b14904e566cf734f2c0e34235f93ae 100644
--- a/README
+++ b/README
@@ -22,12 +22,12 @@ Valid options are:
   -d          Dry run. Read the parameter file, allocate memory but does not read 
               the particles from ICs and exit before the start of time integration.
               Allows user to check validy of parameter and IC files as well as memory limits.
+  -D          Always drift all particles even the ones far from active particles.
   -e          Enable floating-point exceptions (debugging mode)
   -f    {int} Overwrite the CPU frequency (Hz) to be used for time measurements
   -g          Run with an external gravitational potential
   -G          Run with self-gravity
-  -n    {int} Execute a fixed number of time steps. When unset use the time_end
-              parameter to stop. 
+  -n    {int} Execute a fixed number of time steps. When unset use the time_end parameter to stop. 
   -s          Run with SPH
   -t    {int} The number of threads to use on each MPI rank. Defaults to 1 if not specified.
   -v     [12] Increase the level of verbosity 1: MPI-rank 0 writes 
diff --git a/examples/main.c b/examples/main.c
index e1eace107a42e0b7df948fd1183e32ed6b9678a4..561f20ae3fefa723c0fb2c0677f12b9eae724306 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -67,6 +67,8 @@ void print_help_message() {
   printf("  %2s %8s %s\n", "", "",
          "Allows user to check validy of parameter and IC files as well as "
          "memory limits.");
+  printf("  %2s %8s %s\n", "-D", "",
+         "Always drift all particles even the ones far from active particles.");
   printf("  %2s %8s %s\n", "-e", "",
          "Enable floating-point exceptions (debugging mode)");
   printf("  %2s %8s %s\n", "-f", "{int}",
@@ -75,7 +77,8 @@ void print_help_message() {
          "Run with an external gravitational potential");
   printf("  %2s %8s %s\n", "-G", "", "Run with self-gravity");
   printf("  %2s %8s %s\n", "-n", "{int}",
-         "Execute a fixed number of time steps");
+         "Execute a fixed number of time steps. When unset use the time_end "
+         "parameter to stop.");
   printf("  %2s %8s %s\n", "-s", "", "Run with SPH");
   printf("  %2s %8s %s\n", "-t", "{int}",
          "The number of threads to use on each MPI rank. Defaults to 1 if not "
@@ -146,6 +149,7 @@ int main(int argc, char *argv[]) {
   int with_self_gravity = 0;
   int with_hydro = 0;
   int with_fp_exceptions = 0;
+  int with_drift_all = 0;
   int verbose = 0;
   int nr_threads = 1;
   char paramFileName[200] = "";
@@ -153,7 +157,7 @@ int main(int argc, char *argv[]) {
 
   /* Parse the parameters */
   int c;
-  while ((c = getopt(argc, argv, "acdef:gGhn:st:v:y:")) != -1) switch (c) {
+  while ((c = getopt(argc, argv, "acdDef:gGhn:st:v:y:")) != -1) switch (c) {
       case 'a':
         with_aff = 1;
         break;
@@ -163,6 +167,9 @@ int main(int argc, char *argv[]) {
       case 'd':
         dry_run = 1;
         break;
+      case 'D':
+        with_drift_all = 1;
+        break;
       case 'e':
         with_fp_exceptions = 1;
         break;
@@ -438,6 +445,7 @@ int main(int argc, char *argv[]) {
 
   /* Construct the engine policy */
   int engine_policies = ENGINE_POLICY | engine_policy_steal;
+  if (with_drift_all) engine_policies |= engine_policy_drift_all;
   if (with_hydro) engine_policies |= engine_policy_hydro;
   if (with_self_gravity) engine_policies |= engine_policy_self_gravity;
   if (with_external_gravity) engine_policies |= engine_policy_external_gravity;
diff --git a/src/engine.c b/src/engine.c
index 348b8b02a7d5032d8bd8a821b56c6edf443f3df3..75cccd7eeede7ae061aad8aec2e8bd4d3fb8d26e 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -68,7 +68,7 @@
 #include "units.h"
 #include "version.h"
 
-const char *engine_policy_names[13] = {"none",
+const char *engine_policy_names[14] = {"none",
                                        "rand",
                                        "steal",
                                        "keep",
@@ -80,7 +80,8 @@ const char *engine_policy_names[13] = {"none",
                                        "hydro",
                                        "self_gravity",
                                        "external_gravity",
-                                       "cosmology_integration"};
+                                       "cosmology_integration",
+                                       "drift_all"};
 
 /** The rank of the engine as a global variable (for messages). */
 int engine_rank;
@@ -2093,7 +2094,9 @@ void engine_prepare(struct engine *e) {
     e->drift_all = 1;
     threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells,
                    e->s->nr_cells, sizeof(struct cell), 1, e);
-    e->drift_all = 0;
+
+    /* Restore the default drifting policy */
+    e->drift_all = (e->policy & engine_policy_drift_all);
 
     /* And now rebuild */
     engine_rebuild(e);
@@ -2479,7 +2482,9 @@ void engine_step(struct engine *e) {
     e->drift_all = 1;
     threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells,
                    e->s->nr_cells, sizeof(struct cell), 1, e);
-    e->drift_all = 0;
+
+    /* Restore the default drifting policy */
+    e->drift_all = (e->policy & engine_policy_drift_all);
 
     /* Dump... */
     engine_dump_snapshot(e);
diff --git a/src/engine.h b/src/engine.h
index 18055279e9c08f298a3c171ea01f9e891456bca4..3d688925e73a9a07bbb02859be96bc84cfd2ae54 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -61,7 +61,8 @@ enum engine_policy {
   engine_policy_hydro = (1 << 8),
   engine_policy_self_gravity = (1 << 9),
   engine_policy_external_gravity = (1 << 10),
-  engine_policy_cosmology = (1 << 11)
+  engine_policy_cosmology = (1 << 11),
+  engine_policy_drift_all = (1 << 12)
 };
 
 extern const char *engine_policy_names[];
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index bec6b004a483e094e8ef299956a28ddfc7007a20..0108e0663c0d84ff6b5698456f6be34d5ee08c14 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -432,7 +432,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float balsara_j = pj->force.balsara;
 
   /* Are the particles moving towards each others ? */
-  const float omega_ij = (dvdr < 0.f) ? dvdr: 0.f;
+  const float omega_ij = (dvdr < 0.f) ? dvdr : 0.f;
   const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
 
   /* Signal velocity */
@@ -707,7 +707,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float balsara_j = pj->force.balsara;
 
   /* Are the particles moving towards each others ? */
-  const float omega_ij = (dvdr < 0.f) ? dvdr: 0.f;
+  const float omega_ij = (dvdr < 0.f) ? dvdr : 0.f;
   const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
 
   /* Signal velocity */
diff --git a/src/runner.c b/src/runner.c
index 13c8fa6ec14c26eb1da5530aee9cdb5f4ab7adf0..8434b3445bf3cebf6d7f5d1110c54a9ccdefb5a1 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -487,8 +487,8 @@ void runner_do_ghost(struct runner *r, struct cell *c) {
           h_corr = (target_wcount - p->density.wcount) / p->density.wcount_dh;
 
           /* Truncate to the range [ -p->h/2 , p->h ]. */
-	  h_corr = (h_corr < p->h) ? h_corr : p->h;
-	  h_corr = (h_corr > -0.5f * p->h) ? h_corr : -0.5f * p->h;
+          h_corr = (h_corr < p->h) ? h_corr : p->h;
+          h_corr = (h_corr > -0.5f * p->h) ? h_corr : -0.5f * p->h;
         }
 
         /* Did we get the right number density? */
@@ -626,7 +626,7 @@ static void runner_do_drift(struct cell *c, struct engine *e) {
       const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
                         gp->x_diff[1] * gp->x_diff[1] +
                         gp->x_diff[2] * gp->x_diff[2];
-      dx2_max = (dx2_max  > dx2) ? dx2_max : dx2;
+      dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
     }
 
     /* Loop over all the particles in the cell (more work for these !) */
@@ -644,7 +644,7 @@ static void runner_do_drift(struct cell *c, struct engine *e) {
       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];
-      dx2_max = (dx2_max  > dx2) ? dx2_max : dx2;
+      dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
 
       /* Maximal smoothing length */
       h_max = (h_max > p->h) ? h_max : p->h;
diff --git a/src/serial_io.c b/src/serial_io.c
index 83f1f357468a80a820cfbd95a6afcc95198203bf..d0618ea93e07168274dff02f31ae95c89503859b 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -225,11 +225,11 @@ void prepareArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
   }
 
   /* Impose data compression */
-  if(e->snapshotCompression > 0) {
+  if (e->snapshotCompression > 0) {
     h_err = H5Pset_deflate(h_prop, e->snapshotCompression);
     if (h_err < 0) {
       error("Error while setting compression options for field '%s'.",
-	    props.name);
+            props.name);
     }
   }
 
diff --git a/src/timestep.h b/src/timestep.h
index 582f10a658e0be595a22e1066341de3803afb7fc..d92f88d06451892ce47db4b9468db9714bb52baa 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -74,7 +74,8 @@ __attribute__((always_inline)) INLINE static int get_gpart_timestep(
   /*     gravity_compute_timestep_self(e->physical_constants, gp); */
   const float new_dt_self = FLT_MAX;  // MATTHIEU
 
-  float new_dt = (new_dt_external < new_dt_self) ? new_dt_external : new_dt_self;
+  float new_dt =
+      (new_dt_external < new_dt_self) ? new_dt_external : new_dt_self;
 
   /* Limit timestep within the allowed range */
   new_dt = (new_dt < e->dt_max) ? new_dt : e->dt_max;
@@ -111,11 +112,12 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
     /*     gravity_compute_timestep_self(e->physical_constants, p->gpart); */
     const float new_dt_self = FLT_MAX;  // MATTHIEU
 
-    new_dt_grav = (new_dt_external < new_dt_self) ? new_dt_external : new_dt_self;
+    new_dt_grav =
+        (new_dt_external < new_dt_self) ? new_dt_external : new_dt_self;
   }
 
   /* Final time-step is minimum of hydro and gravity */
-  float new_dt = (new_dt_hydro < new_dt_grav) ? new_dt_hydro: new_dt_grav;
+  float new_dt = (new_dt_hydro < new_dt_grav) ? new_dt_hydro : new_dt_grav;
 
   /* Limit change in h */
   const float dt_h_change =