diff --git a/src/engine.c b/src/engine.c
index 242832c67d4fdb810d3e9218a3fbc68f23d3862a..2f532870c8d4201358da9df5e78a06373b4e8d4e 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -167,6 +167,7 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
       scheduler_addunlock(s, c->drift, c->init);
 
       if (is_self_gravity) {
+
         /* Gravity non-neighbouring pm calculations */
         c->grav_long_range = scheduler_addtask(
             s, task_type_grav_long_range, task_subtype_none, 0, 0, c, NULL, 0);
@@ -175,8 +176,15 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
         c->grav_top_level = scheduler_addtask(
             s, task_type_grav_top_level, task_subtype_none, 0, 0, c, NULL, 0);
 
+        /* Gravity recursive down-pass */
+        c->grav_down = scheduler_addtask(s, task_type_grav_down,
+                                         task_subtype_none, 0, 0, c, NULL, 0);
+
         scheduler_addunlock(s, c->init, c->grav_long_range);
         scheduler_addunlock(s, c->init, c->grav_top_level);
+        scheduler_addunlock(s, c->grav_long_range, c->grav_down);
+        scheduler_addunlock(s, c->grav_top_level, c->grav_down);
+        scheduler_addunlock(s, c->grav_down, c->kick2);
       }
 
       /* Generate the ghost task. */
@@ -1799,8 +1807,7 @@ static inline void engine_make_self_gravity_dependencies(
 
   /* init --> gravity --> grav_down --> kick */
   scheduler_addunlock(sched, c->super->init, gravity);
-  // scheduler_addunlock(sched, gravity, c->super->grav_down);
-  scheduler_addunlock(sched, gravity, c->super->kick2);
+  scheduler_addunlock(sched, gravity, c->super->grav_down);
 }
 
 /**
@@ -2169,27 +2176,28 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
  */
 void engine_make_gravityrecursive_tasks(struct engine *e) {
 
-  struct space *s = e->s;
-  struct scheduler *sched = &e->sched;
-  const int nodeID = e->nodeID;
-  const int nr_cells = s->nr_cells;
-  struct cell *cells = s->cells_top;
+  /* struct space *s = e->s; */
+  /* struct scheduler *sched = &e->sched; */
+  /* const int nodeID = e->nodeID; */
+  /* const int nr_cells = s->nr_cells; */
+  /* struct cell *cells = s->cells_top; */
 
-  for (int k = 0; k < nr_cells; k++) {
+  /* for (int k = 0; k < nr_cells; k++) { */
 
-    /* Only do this for local cells containing gravity particles */
-    if (cells[k].nodeID == nodeID && cells[k].gcount > 0) {
+  /*   /\* Only do this for local cells containing gravity particles *\/ */
+  /*   if (cells[k].nodeID == nodeID && cells[k].gcount > 0) { */
 
-      /* Create tasks at top level. */
-      struct task *up = NULL;
-      struct task *down =
-          scheduler_addtask(sched, task_type_grav_down, task_subtype_none, 0, 0,
-                            &cells[k], NULL, 0);
+  /*     /\* Create tasks at top level. *\/ */
+  /*     struct task *up = NULL; */
+  /*     struct task *down = NULL; */
+  /*         /\* scheduler_addtask(sched, task_type_grav_down,
+   * task_subtype_none, 0, 0, *\/ */
+  /*         /\*                   &cells[k], NULL, 0); *\/ */
 
-      /* Push tasks down the cell hierarchy. */
-      engine_addtasks_grav(e, &cells[k], up, down);
-    }
-  }
+  /*     /\* Push tasks down the cell hierarchy. *\/ */
+  /*     engine_addtasks_grav(e, &cells[k], up, down); */
+  /*   } */
+  /* } */
 }
 
 /**
@@ -2842,6 +2850,8 @@ void engine_skip_force_and_kick(struct engine *e) {
     if (t->type == task_type_drift || t->type == task_type_kick1 ||
         t->type == task_type_kick2 || t->type == task_type_timestep ||
         t->subtype == task_subtype_force || t->subtype == task_subtype_grav ||
+        t->type == task_type_grav_long_range ||
+        t->type == task_type_grav_top_level || t->type == task_type_grav_down ||
         t->type == task_type_cooling || t->type == task_type_sourceterms)
       t->skip = 1;
   }
diff --git a/src/runner.c b/src/runner.c
index 02774b9c587349a0554e9534d9e2364948d525f9..afd14db0206cc7fad92b07c6bc83a6eb8c1fcc2d 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -113,7 +113,7 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
 
 /* Import the gravity loop functions. */
 #include "runner_doiact_fft.h"
-//#include "runner_doiact_grav.h"
+#include "runner_doiact_grav.h"
 
 /**
  * @brief Perform source terms
@@ -1350,11 +1350,11 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
 #ifdef SWIFT_DEBUG_CHECKS
       if (e->policy & engine_policy_self_gravity) {
         gp->mass_interacted += gp->mass;
-        if (gp->mass_interacted != e->s->total_mass)
+        if (fabs(gp->mass_interacted - e->s->total_mass) > gp->mass)
           error(
               "g-particle did not interact gravitationally with all other "
-              "particles gp->mass_interacted=%e, total_mass=%e",
-              gp->mass_interacted, e->s->total_mass);
+              "particles gp->mass_interacted=%e, total_mass=%e, gp->mass=%e",
+              gp->mass_interacted, e->s->total_mass, gp->mass);
       }
 #endif
     }
@@ -1686,8 +1686,7 @@ void *runner_main(void *data) {
           else if (t->subtype == task_subtype_force)
             runner_doself2_force(r, ci);
           else if (t->subtype == task_subtype_grav)
-            // runner_doself_grav(r, ci, 1);
-            ;
+            runner_doself_grav(r, ci, 1);
           else if (t->subtype == task_subtype_external_grav)
             runner_do_grav_external(r, ci, 1);
           else
@@ -1704,8 +1703,7 @@ void *runner_main(void *data) {
           else if (t->subtype == task_subtype_force)
             runner_dopair2_force(r, ci, cj);
           else if (t->subtype == task_subtype_grav)
-            // runner_dopair_grav(r, ci, cj, 1);#
-            ;
+            runner_dopair_grav(r, ci, cj, 1);
           else
             error("Unknown/invalid task subtype (%d).", t->subtype);
           break;
@@ -1720,8 +1718,7 @@ void *runner_main(void *data) {
           else if (t->subtype == task_subtype_force)
             runner_dosub_self2_force(r, ci, 1);
           else if (t->subtype == task_subtype_grav)
-            // runner_dosub_grav(r, ci, cj, 1);#
-            ;
+            runner_dosub_grav(r, ci, cj, 1);
           else if (t->subtype == task_subtype_external_grav)
             runner_do_grav_external(r, ci, 1);
           else
@@ -1738,8 +1735,7 @@ void *runner_main(void *data) {
           else if (t->subtype == task_subtype_force)
             runner_dosub_pair2_force(r, ci, cj, t->flags, 1);
           else if (t->subtype == task_subtype_grav)
-            // runner_dosub_grav(r, ci, cj, 1);
-            ;
+            runner_dosub_grav(r, ci, cj, 1);
           else
             error("Unknown/invalid task subtype (%d).", t->subtype);
           break;
@@ -1806,7 +1802,7 @@ void *runner_main(void *data) {
           // runner_do_grav_top_level(r);
           break;
         case task_type_grav_long_range:
-          // runner_do_grav_fft(r);
+          runner_do_grav_long_range(r, t->ci, 1);
           break;
         case task_type_cooling:
           if (e->policy & engine_policy_cooling) runner_do_end_force(r, ci, 1);
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 9d2606ceb06fd6d32592010376e867a6ae582bf0..9498caeb0b4b15cc849160375a50ef984393d847 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -34,23 +34,23 @@
  */
 void runner_do_grav_up(struct runner *r, struct cell *c) {
 
-  if (c->split) { /* Regular node */
-
-    /* Recurse. */
-    for (int k = 0; k < 8; k++)
-      if (c->progeny[k] != NULL) runner_do_grav_up(r, c->progeny[k]);
-
-    /* Collect the multipoles from the progeny. */
-    multipole_reset(&c->multipole);
-    for (int k = 0; k < 8; k++) {
-      if (c->progeny[k] != NULL)
-        multipole_add(&c->multipole, &c->progeny[k]->multipole);
-    }
-
-  } else { /* Leaf node. */
-    /* Just construct the multipole from the gparts. */
-    multipole_init(&c->multipole, c->gparts, c->gcount);
-  }
+  /* if (c->split) { /\* Regular node *\/ */
+
+  /*   /\* Recurse. *\/ */
+  /*   for (int k = 0; k < 8; k++) */
+  /*     if (c->progeny[k] != NULL) runner_do_grav_up(r, c->progeny[k]); */
+
+  /*   /\* Collect the multipoles from the progeny. *\/ */
+  /*   multipole_reset(&c->multipole); */
+  /*   for (int k = 0; k < 8; k++) { */
+  /*     if (c->progeny[k] != NULL) */
+  /*       multipole_add(&c->multipole, &c->progeny[k]->multipole); */
+  /*   } */
+
+  /* } else { /\* Leaf node. *\/ */
+  /*   /\* Just construct the multipole from the gparts. *\/ */
+  /*   multipole_init(&c->multipole, c->gparts, c->gcount); */
+  /* } */
 }
 
 /**
@@ -68,16 +68,15 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
   const struct engine *e = r->e;
   const int gcount = ci->gcount;
   struct gpart *restrict gparts = ci->gparts;
-  const struct multipole multi = cj->multipole;
+  const struct multipole *multi = cj->multipole;
   const float rlr_inv = 1. / (const_gravity_a_smooth * ci->super->width[0]);
 
   TIMER_TIC;
 
 #ifdef SWIFT_DEBUG_CHECKS
-  if (gcount == 0) error("Empty cell!");  // MATTHIEU sanity check
+  if (gcount == 0) error("Empty cell!");
 
-  if (multi.mass == 0.0)  // MATTHIEU sanity check
-    error("Multipole does not seem to have been set.");
+  if (multi->mass == 0.0) error("Multipole does not seem to have been set.");
 #endif
 
   /* Anything to do here? */
@@ -105,13 +104,17 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
     if (!gpart_is_active(gp, e)) continue;
 
     /* Compute the pairwise distance. */
-    const float dx[3] = {multi.CoM[0] - gp->x[0],   // x
-                         multi.CoM[1] - gp->x[1],   // y
-                         multi.CoM[2] - gp->x[2]};  // z
+    const float dx[3] = {multi->CoM[0] - gp->x[0],   // x
+                         multi->CoM[1] - gp->x[1],   // y
+                         multi->CoM[2] - gp->x[2]};  // z
     const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
 
     /* Interact !*/
-    runner_iact_grav_pm(rlr_inv, r2, dx, gp, &multi);
+    runner_iact_grav_pm(rlr_inv, r2, dx, gp, multi);
+
+#ifdef SWIFT_DEBUG_CHECKS
+    gp->mass_interacted += multi->mass;
+#endif
   }
 
   TIMER_TOC(timer_dopair_grav_pm);
@@ -195,18 +198,30 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
 
         runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
 
+#ifdef SWIFT_DEBUG_CHECKS
+        gpi->mass_interacted += gpj->mass;
+        gpj->mass_interacted += gpi->mass;
+#endif
+
       } else {
 
         if (gpart_is_active(gpi, e)) {
 
           runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
 
+#ifdef SWIFT_DEBUG_CHECKS
+          gpi->mass_interacted += gpj->mass;
+#endif
         } else if (gpart_is_active(gpj, e)) {
 
           dx[0] = -dx[0];
           dx[1] = -dx[1];
           dx[2] = -dx[2];
           runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
+
+#ifdef SWIFT_DEBUG_CHECKS
+          gpj->mass_interacted += gpi->mass;
+#endif
         }
       }
     }
@@ -277,18 +292,31 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
 
         runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
 
+#ifdef SWIFT_DEBUG_CHECKS
+        gpi->mass_interacted += gpj->mass;
+        gpj->mass_interacted += gpi->mass;
+#endif
+
       } else {
 
         if (gpart_is_active(gpi, e)) {
 
           runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
 
+#ifdef SWIFT_DEBUG_CHECKS
+          gpi->mass_interacted += gpj->mass;
+#endif
+
         } else if (gpart_is_active(gpj, e)) {
 
           dx[0] = -dx[0];
           dx[1] = -dx[1];
           dx[2] = -dx[2];
           runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
+
+#ifdef SWIFT_DEBUG_CHECKS
+          gpj->mass_interacted += gpi->mass;
+#endif
         }
       }
     }
@@ -466,7 +494,8 @@ static void runner_dosub_grav(struct runner *r, struct cell *ci,
   }
 }
 
-static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) {
+static void runner_do_grav_long_range(struct runner *r, struct cell *ci,
+                                      int timer) {
 
 #if ICHECK > 0
   for (int pid = 0; pid < ci->gcount; pid++) {
@@ -485,10 +514,10 @@ static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) {
   const struct engine *e = r->e;
   struct cell *cells = e->s->cells_top;
   const int nr_cells = e->s->nr_cells;
-  const double max_d =
-      const_gravity_a_smooth * const_gravity_r_cut * ci->width[0];
-  const double max_d2 = max_d * max_d;
-  const double pos_i[3] = {ci->loc[0], ci->loc[1], ci->loc[2]};
+  /* const double max_d = */
+  /*     const_gravity_a_smooth * const_gravity_r_cut * ci->width[0]; */
+  /* const double max_d2 = max_d * max_d; */
+  // const double pos_i[3] = {ci->loc[0], ci->loc[1], ci->loc[2]};
 
   /* Anything to do here? */
   if (!cell_is_active(ci, e)) return;
@@ -501,12 +530,12 @@ static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) {
 
     if (ci == cj) continue;
 
-    const double dx[3] = {cj->loc[0] - pos_i[0],   // x
-                          cj->loc[1] - pos_i[1],   // y
-                          cj->loc[2] - pos_i[2]};  // z
-    const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+    /* const double dx[3] = {cj->loc[0] - pos_i[0],   // x */
+    /*                       cj->loc[1] - pos_i[1],   // y */
+    /*                       cj->loc[2] - pos_i[2]};  // z */
+    /* const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; */
 
-    if (r2 > max_d2) continue;
+    // if (r2 > max_d2) continue;
 
     if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_pm(r, ci, cj);
   }