diff --git a/src/engine.c b/src/engine.c
index 05b0739c5555a9843f9a72d13b39a7e7831386fa..d5fe48bf1f75d3afe17efef60391d9a9f099af48 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1340,28 +1340,31 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
 void engine_make_gravityinteraction_tasks(struct engine *e) {
 
   struct space *s = e->s;
-  struct scheduler *sched = &e->sched;
+  /* struct scheduler *sched = &e->sched; */
   const int nr_cells = s->nr_cells;
-  struct cell *cells = s->cells;
+  /* struct cell *cells = s->cells; */
 
   /* Loop over all cells. */
   for (int i = 0; i < nr_cells; i++) {
 
-    /* If it has gravity particles, add a self-task */
-    if (cells[i].gcount > 0) {
-      scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, -1, 0,
-                        &cells[i], NULL, 0);
-
-      /* Loop over all remainding cells */
-      for (int j = i + 1; j < nr_cells; j++) {
-
-        /* If that other cell has gravity parts, add a pair interaction */
-        if (cells[j].gcount > 0) {
-          scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, -1, 0,
-                            &cells[i], &cells[j], 0);
-        }
-      }
-    }
+    /* /\* If it has gravity particles, add a self-task *\/ */
+    /* if (cells[i].gcount > 0) { */
+    /*   scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, -1, 0,
+     */
+    /*                     &cells[i], NULL, 0); */
+
+    /*   /\* Loop over all remainding cells *\/ */
+    /*   for (int j = i + 1; j < nr_cells; j++) { */
+
+    /*     /\* If that other cell has gravity parts, add a pair interaction *\/
+     */
+    /*     if (cells[j].gcount > 0) { */
+    /*       scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, -1,
+     * 0, */
+    /*                         &cells[i], &cells[j], 0); */
+    /*     } */
+    /*   } */
+    /* } */
   }
 }
 
@@ -1390,9 +1393,12 @@ void engine_make_gravityrecursive_tasks(struct engine *e) {
       struct task *up =
           scheduler_addtask(sched, task_type_grav_up, task_subtype_none, 0, 0,
                             &cells[k], NULL, 0);
-      struct task *down =
-          scheduler_addtask(sched, task_type_grav_down, task_subtype_none, 0, 0,
-                            &cells[k], NULL, 0);
+
+      struct task *down = NULL;
+      /* struct task *down = */
+      /*     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);
@@ -1428,8 +1434,8 @@ void engine_maketasks(struct engine *e) {
   engine_make_hydroloop_tasks(e);
 
   /* Add the gravity mm tasks. */
-  if (e->policy & engine_policy_self_gravity)
-    engine_make_gravityinteraction_tasks(e);
+  /* if (e->policy & engine_policy_self_gravity) */
+  /*   engine_make_gravityinteraction_tasks(e); */
 
   /* Split the tasks. */
   scheduler_splittasks(sched);
diff --git a/src/runner.c b/src/runner.c
index 77153d0608ef086fe592a261e08dd4c701eb3b13..c096d2e480a8a5e468cb2ca79971b419a6ee801e 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1340,8 +1340,8 @@ void *runner_main(void *data) {
             runner_dosub1_density(r, ci, cj, t->flags, 1);
           else if (t->subtype == task_subtype_force)
             runner_dosub2_force(r, ci, cj, t->flags, 1);
-          else if (t->subtype == task_subtype_grav)
-            runner_dosub_grav(r, ci, cj, 1);
+          /* else if (t->subtype == task_subtype_grav) */
+          /*   runner_dosub_grav(r, ci, cj, 1); */
           else
             error("Unknown task subtype.");
           break;
@@ -1362,21 +1362,21 @@ void *runner_main(void *data) {
         case task_type_recv:
           runner_dorecv_cell(r, ci, 1);
           break;
-        case task_type_grav_pp:
-          if (t->cj == NULL)
-            runner_doself_grav(r, t->ci);
-          else
-            runner_dopair_grav(r, t->ci, t->cj);
-          break;
-        case task_type_grav_mm:
-          runner_dograv_mm(r, t->ci, t->cj);
-          break;
+        /* case task_type_grav_pp: */
+        /*   if (t->cj == NULL) */
+        /*     runner_doself_grav(r, t->ci); */
+        /*   else */
+        /*     runner_dopair_grav(r, t->ci, t->cj); */
+        /*   break; */
+        /* case task_type_grav_mm: */
+        /*   runner_dograv_mm(r, t->ci, t->cj); */
+        /*   break; */
         case task_type_grav_up:
           runner_dograv_up(r, t->ci);
           break;
-        case task_type_grav_down:
-          runner_dograv_down(r, t->ci);
-          break;
+        /* case task_type_grav_down: */
+        /*   runner_dograv_down(r, t->ci); */
+        /*   break; */
         case task_type_grav_external:
           runner_dograv_external(r, t->ci, 1);
           break;
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 02626295a49f314fef840bc044a476f5c9cf332d..19e77879501c1e34a54e0d436e9a140a9b273d1a 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -24,177 +24,177 @@
 #include "clocks.h"
 #include "part.h"
 
-/**
- * @brief Compute the sorted gravity interactions between a cell pair.
- *
- * @param r The #runner.
- * @param ci The first #cell.
- * @param cj The second #cell.
- */
-
-void runner_dopair_grav_new(struct runner *r, struct cell *ci,
-                            struct cell *cj) {
-
-  struct engine *restrict e = r->e;
-  int pid, pjd, k, sid;
-  double rshift, shift[3] = {0.0, 0.0, 0.0}, nshift[3];
-  struct entry *restrict sort_i, *restrict sort_j;
-  struct gpart *restrict pi, *restrict pj, *restrict parts_i, *restrict parts_j;
-  double pix[3];
-  float dx[3], r2, h_max, di, dj;
-  int count_i, count_j, cnj, cnj_new;
-  const int ti_current = e->ti_current;
-  struct multipole m;
-#ifdef VECTORIZE
-  int icount = 0;
-  float r2q[VEC_SIZE] __attribute__((aligned(16)));
-  float dxq[3 * VEC_SIZE] __attribute__((aligned(16)));
-  struct gpart *piq[VEC_SIZE], *pjq[VEC_SIZE];
-#endif
-  TIMER_TIC
-
-  /* Anything to do here? */
-  if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return;
-
-  /* Get the sort ID. */
-  sid = space_getsid(e->s, &ci, &cj, shift);
-
-  /* Make sure the cells are sorted. */
-  runner_dogsort(r, ci, (1 << sid), 0);
-  runner_dogsort(r, cj, (1 << sid), 0);
-
-  /* Have the cells been sorted? */
-  if (!(ci->gsorted & (1 << sid)) || !(cj->gsorted & (1 << sid)))
-    error("Trying to interact unsorted cells.");
-
-  /* Get the cutoff shift. */
-  for (rshift = 0.0, k = 0; k < 3; k++)
-    rshift += shift[k] * runner_shift[3 * sid + k];
-
-  /* Pick-out the sorted lists. */
-  sort_i = &ci->gsort[sid * (ci->count + 1)];
-  sort_j = &cj->gsort[sid * (cj->count + 1)];
-
-  /* Get some other useful values. */
-  h_max =
-      sqrtf(ci->h[0] * ci->h[0] + ci->h[1] * ci->h[1] + ci->h[2] * ci->h[2]) *
-      const_theta_max;
-  count_i = ci->gcount;
-  count_j = cj->gcount;
-  parts_i = ci->gparts;
-  parts_j = cj->gparts;
-  cnj = count_j;
-  multipole_reset(&m);
-  nshift[0] = -shift[0];
-  nshift[1] = -shift[1];
-  nshift[2] = -shift[2];
-
-  /* Loop over the parts in ci. */
-  for (pid = count_i - 1; pid >= 0; pid--) {
-
-    /* Get a hold of the ith part in ci. */
-    pi = &parts_i[sort_i[pid].i];
-    if (pi->ti_end > ti_current) continue;
-    di = sort_i[pid].d + h_max - rshift;
-
-    for (k = 0; k < 3; k++) pix[k] = pi->x[k] - shift[k];
-
-    /* Loop over the parts in cj. */
-    for (pjd = 0; pjd < cnj && sort_j[pjd].d < di; pjd++) {
-
-      /* Get a pointer to the jth particle. */
-      pj = &parts_j[sort_j[pjd].i];
-
-      /* Compute the pairwise distance. */
-      r2 = 0.0f;
-      for (k = 0; k < 3; k++) {
-        dx[k] = pix[k] - pj->x[k];
-        r2 += dx[k] * dx[k];
-      }
-
-#ifndef VECTORIZE
-
-      // if ( pi->part->id == 3473472412525 || pj->part->id == 3473472412525 )
-      //     message( "interacting particles pi=%lli and pj=%lli with r=%.3e in
-      // cells %lli/%lli." , pi->part->id , pj->part->id , sqrtf(r2) , ((long
-      // long int)ci) / sizeof(struct cell) , ((long long int)cj) /
-      // sizeof(struct cell) );
+/* /\** */
+/*  * @brief Compute the sorted gravity interactions between a cell pair. */
+/*  * */
+/*  * @param r The #runner. */
+/*  * @param ci The first #cell. */
+/*  * @param cj The second #cell. */
+/*  *\/ */
+
+/* void runner_dopair_grav_new(struct runner *r, struct cell *ci, */
+/*                             struct cell *cj) { */
+
+/*   struct engine *restrict e = r->e; */
+/*   int pid, pjd, k, sid; */
+/*   double rshift, shift[3] = {0.0, 0.0, 0.0}, nshift[3]; */
+/*   struct entry *restrict sort_i, *restrict sort_j; */
+/*   struct gpart *restrict pi, *restrict pj, *restrict parts_i, *restrict parts_j; */
+/*   double pix[3]; */
+/*   float dx[3], r2, h_max, di, dj; */
+/*   int count_i, count_j, cnj, cnj_new; */
+/*   const int ti_current = e->ti_current; */
+/*   struct multipole m; */
+/* #ifdef VECTORIZE */
+/*   int icount = 0; */
+/*   float r2q[VEC_SIZE] __attribute__((aligned(16))); */
+/*   float dxq[3 * VEC_SIZE] __attribute__((aligned(16))); */
+/*   struct gpart *piq[VEC_SIZE], *pjq[VEC_SIZE]; */
+/* #endif */
+/*   TIMER_TIC */
+
+/*   /\* Anything to do here? *\/ */
+/*   if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return; */
+
+/*   /\* Get the sort ID. *\/ */
+/*   sid = space_getsid(e->s, &ci, &cj, shift); */
+
+/*   /\* Make sure the cells are sorted. *\/ */
+/*   runner_dogsort(r, ci, (1 << sid), 0); */
+/*   runner_dogsort(r, cj, (1 << sid), 0); */
+
+/*   /\* Have the cells been sorted? *\/ */
+/*   if (!(ci->gsorted & (1 << sid)) || !(cj->gsorted & (1 << sid))) */
+/*     error("Trying to interact unsorted cells."); */
+
+/*   /\* Get the cutoff shift. *\/ */
+/*   for (rshift = 0.0, k = 0; k < 3; k++) */
+/*     rshift += shift[k] * runner_shift[3 * sid + k]; */
+
+/*   /\* Pick-out the sorted lists. *\/ */
+/*   sort_i = &ci->gsort[sid * (ci->count + 1)]; */
+/*   sort_j = &cj->gsort[sid * (cj->count + 1)]; */
+
+/*   /\* Get some other useful values. *\/ */
+/*   h_max = */
+/*       sqrtf(ci->h[0] * ci->h[0] + ci->h[1] * ci->h[1] + ci->h[2] * ci->h[2]) * */
+/*       const_theta_max; */
+/*   count_i = ci->gcount; */
+/*   count_j = cj->gcount; */
+/*   parts_i = ci->gparts; */
+/*   parts_j = cj->gparts; */
+/*   cnj = count_j; */
+/*   multipole_reset(&m); */
+/*   nshift[0] = -shift[0]; */
+/*   nshift[1] = -shift[1]; */
+/*   nshift[2] = -shift[2]; */
+
+/*   /\* Loop over the parts in ci. *\/ */
+/*   for (pid = count_i - 1; pid >= 0; pid--) { */
+
+/*     /\* Get a hold of the ith part in ci. *\/ */
+/*     pi = &parts_i[sort_i[pid].i]; */
+/*     if (pi->ti_end > ti_current) continue; */
+/*     di = sort_i[pid].d + h_max - rshift; */
+
+/*     for (k = 0; k < 3; k++) pix[k] = pi->x[k] - shift[k]; */
+
+/*     /\* Loop over the parts in cj. *\/ */
+/*     for (pjd = 0; pjd < cnj && sort_j[pjd].d < di; pjd++) { */
+
+/*       /\* Get a pointer to the jth particle. *\/ */
+/*       pj = &parts_j[sort_j[pjd].i]; */
+
+/*       /\* Compute the pairwise distance. *\/ */
+/*       r2 = 0.0f; */
+/*       for (k = 0; k < 3; k++) { */
+/*         dx[k] = pix[k] - pj->x[k]; */
+/*         r2 += dx[k] * dx[k]; */
+/*       } */
+
+/* #ifndef VECTORIZE */
+
+/*       // if ( pi->part->id == 3473472412525 || pj->part->id == 3473472412525 ) */
+/*       //     message( "interacting particles pi=%lli and pj=%lli with r=%.3e in */
+/*       // cells %lli/%lli." , pi->part->id , pj->part->id , sqrtf(r2) , ((long */
+/*       // long int)ci) / sizeof(struct cell) , ((long long int)cj) / */
+/*       // sizeof(struct cell) ); */
+
+/*       runner_iact_grav(r2, dx, pi, pj); */
+
+/* #else */
+
+/*       /\* Add this interaction to the queue. *\/ */
+/*       r2q[icount] = r2; */
+/*       dxq[3 * icount + 0] = dx[0]; */
+/*       dxq[3 * icount + 1] = dx[1]; */
+/*       dxq[3 * icount + 2] = dx[2]; */
+/*       piq[icount] = pi; */
+/*       pjq[icount] = pj; */
+/*       icount += 1; */
+
+/*       /\* Flush? *\/ */
+/*       if (icount == VEC_SIZE) { */
+/*         runner_iact_vec_grav(r2q, dxq, piq, pjq); */
+/*         icount = 0; */
+/*       } */
 
-      runner_iact_grav(r2, dx, pi, pj);
+/* #endif */
 
-#else
+/*     } /\* loop over the parts in cj. *\/ */
 
-      /* Add this interaction to the queue. */
-      r2q[icount] = r2;
-      dxq[3 * icount + 0] = dx[0];
-      dxq[3 * icount + 1] = dx[1];
-      dxq[3 * icount + 2] = dx[2];
-      piq[icount] = pi;
-      pjq[icount] = pj;
-      icount += 1;
-
-      /* Flush? */
-      if (icount == VEC_SIZE) {
-        runner_iact_vec_grav(r2q, dxq, piq, pjq);
-        icount = 0;
-      }
-
-#endif
+/*     /\* Set the new limit. *\/ */
+/*     cnj_new = pjd; */
 
-    } /* loop over the parts in cj. */
+/*     /\* Add trailing parts to the multipole. *\/ */
+/*     for (pjd = cnj_new; pjd < cnj; pjd++) { */
 
-    /* Set the new limit. */
-    cnj_new = pjd;
+/*       /\* Add the part to the multipole. *\/ */
+/*       multipole_addpart(&m, &parts_j[sort_j[pjd].i]); */
 
-    /* Add trailing parts to the multipole. */
-    for (pjd = cnj_new; pjd < cnj; pjd++) {
+/*     } /\* add trailing parts to the multipole. *\/ */
 
-      /* Add the part to the multipole. */
-      multipole_addpart(&m, &parts_j[sort_j[pjd].i]);
+/*     /\* Set the new cnj. *\/ */
+/*     cnj = cnj_new; */
 
-    } /* add trailing parts to the multipole. */
+/*     /\* Interact the ith particle with the multipole. *\/ */
+/*     multipole_iact_mp(&m, pi, nshift); */
 
-    /* Set the new cnj. */
-    cnj = cnj_new;
+/*   } /\* loop over the parts in ci. *\/ */
 
-    /* Interact the ith particle with the multipole. */
-    multipole_iact_mp(&m, pi, nshift);
+/* #ifdef VECTORIZE */
+/*   /\* Pick up any leftovers. *\/ */
+/*   if (icount > 0) */
+/*     for (k = 0; k < icount; k++) */
+/*       runner_iact_grav(r2q[k], &dxq[3 * k], piq[k], pjq[k]); */
+/* #endif */
 
-  } /* loop over the parts in ci. */
+/*   /\* Re-set the multipole. *\/ */
+/*   multipole_reset(&m); */
 
-#ifdef VECTORIZE
-  /* Pick up any leftovers. */
-  if (icount > 0)
-    for (k = 0; k < icount; k++)
-      runner_iact_grav(r2q[k], &dxq[3 * k], piq[k], pjq[k]);
-#endif
+/*   /\* Loop over the parts in cj and interact with the multipole in ci. *\/ */
+/*   for (pid = count_i - 1, pjd = 0; pjd < count_j; pjd++) { */
 
-  /* Re-set the multipole. */
-  multipole_reset(&m);
+/*     /\* Get the position of pj along the axis. *\/ */
+/*     dj = sort_j[pjd].d - h_max + rshift; */
 
-  /* Loop over the parts in cj and interact with the multipole in ci. */
-  for (pid = count_i - 1, pjd = 0; pjd < count_j; pjd++) {
+/*     /\* Add any left-over parts in cell_i to the multipole. *\/ */
+/*     while (pid >= 0 && sort_i[pid].d < dj) { */
 
-    /* Get the position of pj along the axis. */
-    dj = sort_j[pjd].d - h_max + rshift;
+/*       /\* Add this particle to the multipole. *\/ */
+/*       multipole_addpart(&m, &parts_i[sort_i[pid].i]); */
 
-    /* Add any left-over parts in cell_i to the multipole. */
-    while (pid >= 0 && sort_i[pid].d < dj) {
+/*       /\* Decrease pid. *\/ */
+/*       pid -= 1; */
+/*     } */
 
-      /* Add this particle to the multipole. */
-      multipole_addpart(&m, &parts_i[sort_i[pid].i]);
+/*     /\* Interact pj with the multipole. *\/ */
+/*     multipole_iact_mp(&m, &parts_j[sort_j[pjd].i], shift); */
 
-      /* Decrease pid. */
-      pid -= 1;
-    }
+/*   } /\* loop over the parts in cj and interact with the multipole. *\/ */
 
-    /* Interact pj with the multipole. */
-    multipole_iact_mp(&m, &parts_j[sort_j[pjd].i], shift);
-
-  } /* loop over the parts in cj and interact with the multipole. */
-
-  TIMER_TOC(TIMER_DOPAIR);
-}
+/*   TIMER_TOC(TIMER_DOPAIR); */
+/* } */
 
 /**
  * @brief Compute the recursive upward sweep, i.e. construct the
@@ -231,367 +231,367 @@ void runner_dograv_up(struct runner *r, struct cell *c) {
     multipole_init(&c->multipole, c->gparts, c->gcount);
 }
 
-/**
- * @brief Compute the recursive downward sweep, i.e. apply the multipole
- *        acceleration on all the particles.
- *
- * @param r The #runner.
- * @param c The top-level #cell.
- */
-
-void runner_dograv_down(struct runner *r, struct cell *c) {
-
-  struct multipole *m = &c->multipole;
+/* /\** */
+/*  * @brief Compute the recursive downward sweep, i.e. apply the multipole */
+/*  *        acceleration on all the particles. */
+/*  * */
+/*  * @param r The #runner. */
+/*  * @param c The top-level #cell. */
+/*  *\/ */
+
+/* void runner_dograv_down(struct runner *r, struct cell *c) { */
+
+/*   struct multipole *m = &c->multipole; */
+
+/*   /\* Split? *\/ */
+/*   if (c->split) { */
+
+/*     /\* Apply this cell's acceleration on the multipoles below. *\/ */
+/*     for (int k = 0; k < 8; k++) */
+/*       if (c->progeny[k] != NULL) { */
+/*         struct multipole *mp = &c->progeny[k]->multipole; */
+/*         mp->a[0] += m->a[0]; */
+/*         mp->a[1] += m->a[1]; */
+/*         mp->a[2] += m->a[2]; */
+/*       } */
+
+/*     /\* Recurse. *\/ */
+/*     for (int k = 0; k < 8; k++) */
+/*       if (c->progeny[k] != NULL) runner_dograv_down(r, c->progeny[k]); */
+
+/*   } */
+
+/*   /\* No, leaf node. *\/ */
+/*   else { */
+
+/*     /\* Apply the multipole acceleration to all gparts. *\/ */
+/*     for (int k = 0; k < c->gcount; k++) { */
+/*       struct gpart *p = &c->gparts[k]; */
+/*       p->a_grav[0] += m->a[0]; */
+/*       p->a_grav[1] += m->a[1]; */
+/*       p->a_grav[2] += m->a[2]; */
+/*     } */
+/*   } */
+/* } */
+
+/* /\** */
+/*  * @brief Compute the multipole-multipole interaction between two cells. */
+/*  * */
+/*  * @param r The #runner. */
+/*  * @param ci The first #cell. */
+/*  * @param cj The second #cell. */
+/*  *\/ */
+
+/* void runner_dograv_mm(struct runner *r, struct cell *restrict ci, */
+/*                       struct cell *restrict cj) { */
+
+/*   struct engine *e = r->e; */
+/*   int k; */
+/*   double shift[3] = {0.0, 0.0, 0.0}; */
+/*   float dx[3], theta; */
+
+/*   /\* Compute the shift between the cells. *\/ */
+/*   for (k = 0; k < 3; k++) { */
+/*     dx[k] = cj->loc[k] - ci->loc[k]; */
+/*     if (r->e->s->periodic) { */
+/*       if (dx[k] < -e->s->dim[k] / 2) */
+/*         shift[k] = e->s->dim[k]; */
+/*       else if (dx[k] > e->s->dim[k] / 2) */
+/*         shift[k] = -e->s->dim[k]; */
+/*       dx[k] += shift[k]; */
+/*     } */
+/*   } */
+/*   theta = */
+/*       sqrt((dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]) / */
+/*            (ci->h[0] * ci->h[0] + ci->h[1] * ci->h[1] + ci->h[2] * ci->h[2])); */
+
+/*   /\* Do an MM or an MP/PM? *\/ */
+/*   if (theta > const_theta_max * 4) { */
+
+/*     /\* Update the multipoles. *\/ */
+/*     multipole_iact_mm(&ci->multipole, &cj->multipole, shift); */
+
+/*   } else { */
+
+/*     /\* Interact the multipoles via their parts. *\/ */
+/*     for (k = 0; k < ci->gcount; k++) */
+/*       multipole_iact_mp(&cj->multipole, &ci->gparts[k], shift); */
+/*     for (k = 0; k < cj->gcount; k++) */
+/*       multipole_iact_mp(&ci->multipole, &cj->gparts[k], shift); */
+/*   } */
+/* } */
+
+/* /\** */
+/*  * @brief Compute the interactions between a cell pair. */
+/*  * */
+/*  * @param r The #runner. */
+/*  * @param ci The first #cell. */
+/*  * @param cj The second #cell. */
+/*  *\/ */
+
+/* void runner_dopair_grav(struct runner *r, struct cell *restrict ci, */
+/*                         struct cell *restrict cj) { */
+
+/*   struct engine *e = r->e; */
+/*   int pid, pjd, k, count_i = ci->gcount, count_j = cj->gcount; */
+/*   double shift[3] = {0.0, 0.0, 0.0}; */
+/*   struct gpart *restrict parts_i = ci->gparts, *restrict parts_j = cj->gparts; */
+/*   struct gpart *restrict pi, *restrict pj; */
+/*   double pix[3]; */
+/*   float dx[3], r2; */
+/*   const int ti_current = r->e->ti_current; */
+/* #ifdef VECTORIZE */
+/*   int icount = 0; */
+/*   float r2q[VEC_SIZE] __attribute__((aligned(16))); */
+/*   float dxq[3 * VEC_SIZE] __attribute__((aligned(16))); */
+/*   struct gpart *piq[VEC_SIZE], *pjq[VEC_SIZE]; */
+/* #endif */
+/*   TIMER_TIC */
+
+/*   /\* Anything to do here? *\/ */
+/*   if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return; */
+
+/*   /\* Get the relative distance between the pairs, wrapping. *\/ */
+/*   if (e->s->periodic) */
+/*     for (k = 0; k < 3; k++) { */
+/*       if (cj->loc[k] - ci->loc[k] < -e->s->dim[k] / 2) */
+/*         shift[k] = e->s->dim[k]; */
+/*       else if (cj->loc[k] - ci->loc[k] > e->s->dim[k] / 2) */
+/*         shift[k] = -e->s->dim[k]; */
+/*     } */
+
+/*   /\* Loop over the parts in ci. *\/ */
+/*   for (pid = 0; pid < count_i; pid++) { */
+
+/*     /\* Get a hold of the ith part in ci. *\/ */
+/*     pi = &parts_i[pid]; */
+/*     for (k = 0; k < 3; k++) pix[k] = pi->x[k] - shift[k]; */
+
+/*     /\* Loop over the parts in cj. *\/ */
+/*     for (pjd = 0; pjd < count_j; pjd++) { */
+
+/*       /\* Get a pointer to the jth particle. *\/ */
+/*       pj = &parts_j[pjd]; */
+
+/*       /\* Compute the pairwise distance. *\/ */
+/*       r2 = 0.0f; */
+/*       for (k = 0; k < 3; k++) { */
+/*         dx[k] = pix[k] - pj->x[k]; */
+/*         r2 += dx[k] * dx[k]; */
+/*       } */
+
+/* /\* Compute the interaction. *\/ */
+/* #ifndef VECTORIZE */
+
+/*       // if ( pi->part->id == 3473472412525 || pj->part->id == 3473472412525 ) */
+/*       //     message( "interacting particles pi=%lli and pj=%lli with r=%.3e in */
+/*       // cells %lli/%lli." , pi->part->id , pj->part->id , sqrtf(r2) , ((long */
+/*       // long int)ci) / sizeof(struct cell) , ((long long int)cj) / */
+/*       // sizeof(struct cell) ); */
+
+/*       runner_iact_grav(r2, dx, pi, pj); */
+
+/* #else */
+
+/*       /\* Add this interaction to the queue. *\/ */
+/*       r2q[icount] = r2; */
+/*       dxq[3 * icount + 0] = dx[0]; */
+/*       dxq[3 * icount + 1] = dx[1]; */
+/*       dxq[3 * icount + 2] = dx[2]; */
+/*       piq[icount] = pi; */
+/*       pjq[icount] = pj; */
+/*       icount += 1; */
+
+/*       /\* Flush? *\/ */
+/*       if (icount == VEC_SIZE) { */
+/*         runner_iact_vec_grav(r2q, dxq, piq, pjq); */
+/*         icount = 0; */
+/*       } */
+
+/* #endif */
+
+/*     } /\* loop over the parts in cj. *\/ */
+
+/*   } /\* loop over the parts in ci. *\/ */
+
+/* #ifdef VECTORIZE */
+/*   /\* Pick up any leftovers. *\/ */
+/*   if (icount > 0) */
+/*     for (k = 0; k < icount; k++) */
+/*       runner_iact_grav(r2q[k], &dxq[3 * k], piq[k], pjq[k]); */
+/* #endif */
+
+/*   TIMER_TOC(timer_dopair_grav); */
+/* } */
+
+/* /\** */
+/*  * @brief Compute the interactions within a cell. */
+/*  * */
+/*  * @param r The #runner. */
+/*  * @param c The #cell. */
+/*  *\/ */
+
+/* void runner_doself_grav(struct runner *r, struct cell *restrict c) { */
+
+/*   int pid, pjd, k, count = c->gcount; */
+/*   struct gpart *restrict parts = c->gparts; */
+/*   struct gpart *restrict pi, *restrict pj; */
+/*   double pix[3] = {0.0, 0.0, 0.0}; */
+/*   float dx[3], r2; */
+/*   const int ti_current = r->e->ti_current; */
+/* #ifdef VECTORIZE */
+/*   int icount = 0; */
+/*   float r2q[VEC_SIZE] __attribute__((aligned(16))); */
+/*   float dxq[3 * VEC_SIZE] __attribute__((aligned(16))); */
+/*   struct gpart *piq[VEC_SIZE], *pjq[VEC_SIZE]; */
+/* #endif */
+/*   TIMER_TIC */
+
+/*   /\* Anything to do here? *\/ */
+/*   if (c->ti_end_min > ti_current) return; */
+
+/*   /\* Loop over every part in c. *\/ */
+/*   for (pid = 0; pid < count; pid++) { */
+
+/*     /\* Get a hold of the ith part in ci. *\/ */
+/*     pi = &parts[pid]; */
+/*     for (k = 0; k < 3; k++) pix[k] = pi->x[k]; */
+
+/*     /\* Loop over every other part in c. *\/ */
+/*     for (pjd = pid + 1; pjd < count; pjd++) { */
+
+/*       /\* Get a pointer to the jth particle. *\/ */
+/*       pj = &parts[pjd]; */
+
+/*       /\* Compute the pairwise distance. *\/ */
+/*       r2 = 0.0f; */
+/*       for (k = 0; k < 3; k++) { */
+/*         dx[k] = pix[k] - pj->x[k]; */
+/*         r2 += dx[k] * dx[k]; */
+/*       } */
+
+/* /\* Compute the interaction. *\/ */
+/* #ifndef VECTORIZE */
+
+/*       // if ( pi->part->id == 3473472412525 || pj->part->id == 3473472412525 ) */
+/*       //     message( "interacting particles pi=%lli and pj=%lli with r=%.3e." , */
+/*       // pi->part->id , pj->part->id , sqrtf(r2) ); */
+
+/*       runner_iact_grav(r2, dx, pi, pj); */
+
+/* #else */
+
+/*       /\* Add this interaction to the queue. *\/ */
+/*       r2q[icount] = r2; */
+/*       dxq[3 * icount + 0] = dx[0]; */
+/*       dxq[3 * icount + 1] = dx[1]; */
+/*       dxq[3 * icount + 2] = dx[2]; */
+/*       piq[icount] = pi; */
+/*       pjq[icount] = pj; */
+/*       icount += 1; */
+
+/*       /\* Flush? *\/ */
+/*       if (icount == VEC_SIZE) { */
+/*         runner_iact_vec_grav(r2q, dxq, piq, pjq); */
+/*         icount = 0; */
+/*       } */
+
+/* #endif */
+
+/*     } /\* loop over the remaining parts in c. *\/ */
+
+/*   } /\* loop over the parts in c. *\/ */
+
+/* #ifdef VECTORIZE */
+/*   /\* Pick up any leftovers. *\/ */
+/*   if (icount > 0) */
+/*     for (k = 0; k < icount; k++) */
+/*       runner_iact_grav(r2q[k], &dxq[3 * k], piq[k], pjq[k]); */
+/* #endif */
+
+/*   TIMER_TOC(timer_doself_grav); */
+/* } */
+
+/* /\** */
+/*  * @brief Compute a gravity sub-task. */
+/*  * */
+/*  * @param r The #runner. */
+/*  * @param ci The first #cell. */
+/*  * @param cj The second #cell. */
+/*  * @param gettimer Flag to record timer or not. */
+/*  *\/ */
+
+/* void runner_dosub_grav(struct runner *r, struct cell *ci, struct cell *cj, */
+/*                        int gettimer) { */
+
+/*   int j, k, periodic = r->e->s->periodic; */
+/*   struct space *s = r->e->s; */
+
+/*   TIMER_TIC */
+
+/*   /\* Self-interaction? *\/ */
+/*   if (cj == NULL) { */
+
+/*     /\* If the cell is split, recurse. *\/ */
+/*     if (ci->split) { */
+
+/*       /\* Split this task into tasks on its progeny. *\/ */
+/*       for (j = 0; j < 8; j++) */
+/*         if (ci->progeny[j] != NULL) { */
+/*           runner_dosub_grav(r, ci->progeny[j], NULL, 0); */
+/*           for (k = j + 1; k < 8; k++) */
+/*             if (ci->progeny[k] != NULL) */
+/*               runner_dosub_grav(r, ci->progeny[j], ci->progeny[k], 0); */
+/*         } */
 
-  /* Split? */
-  if (c->split) {
+/*     } */
 
-    /* Apply this cell's acceleration on the multipoles below. */
-    for (int k = 0; k < 8; k++)
-      if (c->progeny[k] != NULL) {
-        struct multipole *mp = &c->progeny[k]->multipole;
-        mp->a[0] += m->a[0];
-        mp->a[1] += m->a[1];
-        mp->a[2] += m->a[2];
-      }
+/*     /\* Otherwise, just make a pp task out of it. *\/ */
+/*     else */
+/*       runner_doself_grav(r, ci); */
 
-    /* Recurse. */
-    for (int k = 0; k < 8; k++)
-      if (c->progeny[k] != NULL) runner_dograv_down(r, c->progeny[k]);
+/*   } */
 
-  }
+/*   /\* Nope, pair. *\/ */
+/*   else { */
+
+/*     /\* Get the opening angle theta. *\/ */
+/*     float dx[3], theta; */
+/*     for (k = 0; k < 3; k++) { */
+/*       dx[k] = fabs(ci->loc[k] - cj->loc[k]); */
+/*       if (periodic && dx[k] > 0.5 * s->dim[k]) dx[k] = -dx[k] + s->dim[k]; */
+/*       if (dx[k] > 0.0f) dx[k] -= ci->h[k]; */
+/*     } */
+/*     theta = (dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]) / */
+/*             (ci->h[0] * ci->h[0] + ci->h[1] * ci->h[1] + ci->h[2] * ci->h[2]); */
+
+/*     /\* Split the interaction? *\/ */
+/*     if (theta < const_theta_max * const_theta_max) { */
 
-  /* No, leaf node. */
-  else {
-
-    /* Apply the multipole acceleration to all gparts. */
-    for (int k = 0; k < c->gcount; k++) {
-      struct gpart *p = &c->gparts[k];
-      p->a_grav[0] += m->a[0];
-      p->a_grav[1] += m->a[1];
-      p->a_grav[2] += m->a[2];
-    }
-  }
-}
+/*       /\* Are both ci and cj split? *\/ */
+/*       if (ci->split && cj->split) { */
+
+/*         /\* Split this task into tasks on its progeny. *\/ */
+/*         for (j = 0; j < 8; j++) */
+/*           if (ci->progeny[j] != NULL) { */
+/*             for (k = 0; k < 8; k++) */
+/*               if (cj->progeny[k] != NULL) */
+/*                 runner_dosub_grav(r, ci->progeny[j], cj->progeny[k], 0); */
+/*           } */
+
+/*       } */
+
+/*       /\* Otherwise, make a pp task out of it. *\/ */
+/*       else */
+/*         runner_dopair_grav(r, ci, cj); */
+
+/*     } */
+
+/*     /\* Otherwise, mm interaction is fine. *\/ */
+/*     else */
+/*       runner_dograv_mm(r, ci, cj); */
+/*   } */
 
-/**
- * @brief Compute the multipole-multipole interaction between two cells.
- *
- * @param r The #runner.
- * @param ci The first #cell.
- * @param cj The second #cell.
- */
-
-void runner_dograv_mm(struct runner *r, struct cell *restrict ci,
-                      struct cell *restrict cj) {
-
-  struct engine *e = r->e;
-  int k;
-  double shift[3] = {0.0, 0.0, 0.0};
-  float dx[3], theta;
-
-  /* Compute the shift between the cells. */
-  for (k = 0; k < 3; k++) {
-    dx[k] = cj->loc[k] - ci->loc[k];
-    if (r->e->s->periodic) {
-      if (dx[k] < -e->s->dim[k] / 2)
-        shift[k] = e->s->dim[k];
-      else if (dx[k] > e->s->dim[k] / 2)
-        shift[k] = -e->s->dim[k];
-      dx[k] += shift[k];
-    }
-  }
-  theta =
-      sqrt((dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]) /
-           (ci->h[0] * ci->h[0] + ci->h[1] * ci->h[1] + ci->h[2] * ci->h[2]));
-
-  /* Do an MM or an MP/PM? */
-  if (theta > const_theta_max * 4) {
-
-    /* Update the multipoles. */
-    multipole_iact_mm(&ci->multipole, &cj->multipole, shift);
-
-  } else {
-
-    /* Interact the multipoles via their parts. */
-    for (k = 0; k < ci->gcount; k++)
-      multipole_iact_mp(&cj->multipole, &ci->gparts[k], shift);
-    for (k = 0; k < cj->gcount; k++)
-      multipole_iact_mp(&ci->multipole, &cj->gparts[k], shift);
-  }
-}
-
-/**
- * @brief Compute the interactions between a cell pair.
- *
- * @param r The #runner.
- * @param ci The first #cell.
- * @param cj The second #cell.
- */
-
-void runner_dopair_grav(struct runner *r, struct cell *restrict ci,
-                        struct cell *restrict cj) {
-
-  struct engine *e = r->e;
-  int pid, pjd, k, count_i = ci->gcount, count_j = cj->gcount;
-  double shift[3] = {0.0, 0.0, 0.0};
-  struct gpart *restrict parts_i = ci->gparts, *restrict parts_j = cj->gparts;
-  struct gpart *restrict pi, *restrict pj;
-  double pix[3];
-  float dx[3], r2;
-  const int ti_current = r->e->ti_current;
-#ifdef VECTORIZE
-  int icount = 0;
-  float r2q[VEC_SIZE] __attribute__((aligned(16)));
-  float dxq[3 * VEC_SIZE] __attribute__((aligned(16)));
-  struct gpart *piq[VEC_SIZE], *pjq[VEC_SIZE];
-#endif
-  TIMER_TIC
-
-  /* Anything to do here? */
-  if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return;
-
-  /* Get the relative distance between the pairs, wrapping. */
-  if (e->s->periodic)
-    for (k = 0; k < 3; k++) {
-      if (cj->loc[k] - ci->loc[k] < -e->s->dim[k] / 2)
-        shift[k] = e->s->dim[k];
-      else if (cj->loc[k] - ci->loc[k] > e->s->dim[k] / 2)
-        shift[k] = -e->s->dim[k];
-    }
-
-  /* Loop over the parts in ci. */
-  for (pid = 0; pid < count_i; pid++) {
-
-    /* Get a hold of the ith part in ci. */
-    pi = &parts_i[pid];
-    for (k = 0; k < 3; k++) pix[k] = pi->x[k] - shift[k];
-
-    /* Loop over the parts in cj. */
-    for (pjd = 0; pjd < count_j; pjd++) {
-
-      /* Get a pointer to the jth particle. */
-      pj = &parts_j[pjd];
-
-      /* Compute the pairwise distance. */
-      r2 = 0.0f;
-      for (k = 0; k < 3; k++) {
-        dx[k] = pix[k] - pj->x[k];
-        r2 += dx[k] * dx[k];
-      }
-
-/* Compute the interaction. */
-#ifndef VECTORIZE
-
-      // if ( pi->part->id == 3473472412525 || pj->part->id == 3473472412525 )
-      //     message( "interacting particles pi=%lli and pj=%lli with r=%.3e in
-      // cells %lli/%lli." , pi->part->id , pj->part->id , sqrtf(r2) , ((long
-      // long int)ci) / sizeof(struct cell) , ((long long int)cj) /
-      // sizeof(struct cell) );
-
-      runner_iact_grav(r2, dx, pi, pj);
-
-#else
-
-      /* Add this interaction to the queue. */
-      r2q[icount] = r2;
-      dxq[3 * icount + 0] = dx[0];
-      dxq[3 * icount + 1] = dx[1];
-      dxq[3 * icount + 2] = dx[2];
-      piq[icount] = pi;
-      pjq[icount] = pj;
-      icount += 1;
-
-      /* Flush? */
-      if (icount == VEC_SIZE) {
-        runner_iact_vec_grav(r2q, dxq, piq, pjq);
-        icount = 0;
-      }
-
-#endif
-
-    } /* loop over the parts in cj. */
-
-  } /* loop over the parts in ci. */
-
-#ifdef VECTORIZE
-  /* Pick up any leftovers. */
-  if (icount > 0)
-    for (k = 0; k < icount; k++)
-      runner_iact_grav(r2q[k], &dxq[3 * k], piq[k], pjq[k]);
-#endif
-
-  TIMER_TOC(timer_dopair_grav);
-}
-
-/**
- * @brief Compute the interactions within a cell.
- *
- * @param r The #runner.
- * @param c The #cell.
- */
-
-void runner_doself_grav(struct runner *r, struct cell *restrict c) {
-
-  int pid, pjd, k, count = c->gcount;
-  struct gpart *restrict parts = c->gparts;
-  struct gpart *restrict pi, *restrict pj;
-  double pix[3] = {0.0, 0.0, 0.0};
-  float dx[3], r2;
-  const int ti_current = r->e->ti_current;
-#ifdef VECTORIZE
-  int icount = 0;
-  float r2q[VEC_SIZE] __attribute__((aligned(16)));
-  float dxq[3 * VEC_SIZE] __attribute__((aligned(16)));
-  struct gpart *piq[VEC_SIZE], *pjq[VEC_SIZE];
-#endif
-  TIMER_TIC
-
-  /* Anything to do here? */
-  if (c->ti_end_min > ti_current) return;
-
-  /* Loop over every part in c. */
-  for (pid = 0; pid < count; pid++) {
-
-    /* Get a hold of the ith part in ci. */
-    pi = &parts[pid];
-    for (k = 0; k < 3; k++) pix[k] = pi->x[k];
-
-    /* Loop over every other part in c. */
-    for (pjd = pid + 1; pjd < count; pjd++) {
-
-      /* Get a pointer to the jth particle. */
-      pj = &parts[pjd];
-
-      /* Compute the pairwise distance. */
-      r2 = 0.0f;
-      for (k = 0; k < 3; k++) {
-        dx[k] = pix[k] - pj->x[k];
-        r2 += dx[k] * dx[k];
-      }
-
-/* Compute the interaction. */
-#ifndef VECTORIZE
-
-      // if ( pi->part->id == 3473472412525 || pj->part->id == 3473472412525 )
-      //     message( "interacting particles pi=%lli and pj=%lli with r=%.3e." ,
-      // pi->part->id , pj->part->id , sqrtf(r2) );
-
-      runner_iact_grav(r2, dx, pi, pj);
-
-#else
-
-      /* Add this interaction to the queue. */
-      r2q[icount] = r2;
-      dxq[3 * icount + 0] = dx[0];
-      dxq[3 * icount + 1] = dx[1];
-      dxq[3 * icount + 2] = dx[2];
-      piq[icount] = pi;
-      pjq[icount] = pj;
-      icount += 1;
-
-      /* Flush? */
-      if (icount == VEC_SIZE) {
-        runner_iact_vec_grav(r2q, dxq, piq, pjq);
-        icount = 0;
-      }
-
-#endif
-
-    } /* loop over the remaining parts in c. */
-
-  } /* loop over the parts in c. */
-
-#ifdef VECTORIZE
-  /* Pick up any leftovers. */
-  if (icount > 0)
-    for (k = 0; k < icount; k++)
-      runner_iact_grav(r2q[k], &dxq[3 * k], piq[k], pjq[k]);
-#endif
-
-  TIMER_TOC(timer_doself_grav);
-}
-
-/**
- * @brief Compute a gravity sub-task.
- *
- * @param r The #runner.
- * @param ci The first #cell.
- * @param cj The second #cell.
- * @param gettimer Flag to record timer or not.
- */
-
-void runner_dosub_grav(struct runner *r, struct cell *ci, struct cell *cj,
-                       int gettimer) {
-
-  int j, k, periodic = r->e->s->periodic;
-  struct space *s = r->e->s;
-
-  TIMER_TIC
-
-  /* Self-interaction? */
-  if (cj == NULL) {
-
-    /* If the cell is split, recurse. */
-    if (ci->split) {
-
-      /* Split this task into tasks on its progeny. */
-      for (j = 0; j < 8; j++)
-        if (ci->progeny[j] != NULL) {
-          runner_dosub_grav(r, ci->progeny[j], NULL, 0);
-          for (k = j + 1; k < 8; k++)
-            if (ci->progeny[k] != NULL)
-              runner_dosub_grav(r, ci->progeny[j], ci->progeny[k], 0);
-        }
-
-    }
-
-    /* Otherwise, just make a pp task out of it. */
-    else
-      runner_doself_grav(r, ci);
-
-  }
-
-  /* Nope, pair. */
-  else {
-
-    /* Get the opening angle theta. */
-    float dx[3], theta;
-    for (k = 0; k < 3; k++) {
-      dx[k] = fabs(ci->loc[k] - cj->loc[k]);
-      if (periodic && dx[k] > 0.5 * s->dim[k]) dx[k] = -dx[k] + s->dim[k];
-      if (dx[k] > 0.0f) dx[k] -= ci->h[k];
-    }
-    theta = (dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]) /
-            (ci->h[0] * ci->h[0] + ci->h[1] * ci->h[1] + ci->h[2] * ci->h[2]);
-
-    /* Split the interaction? */
-    if (theta < const_theta_max * const_theta_max) {
-
-      /* Are both ci and cj split? */
-      if (ci->split && cj->split) {
-
-        /* Split this task into tasks on its progeny. */
-        for (j = 0; j < 8; j++)
-          if (ci->progeny[j] != NULL) {
-            for (k = 0; k < 8; k++)
-              if (cj->progeny[k] != NULL)
-                runner_dosub_grav(r, ci->progeny[j], cj->progeny[k], 0);
-          }
-
-      }
-
-      /* Otherwise, make a pp task out of it. */
-      else
-        runner_dopair_grav(r, ci, cj);
-
-    }
-
-    /* Otherwise, mm interaction is fine. */
-    else
-      runner_dograv_mm(r, ci, cj);
-  }
-
-  if (gettimer) TIMER_TOC(timer_dosub_grav);
-}
+/*   if (gettimer) TIMER_TOC(timer_dosub_grav); */
+/* } */
 #endif /* SWIFT_RUNNER_DOIACT_GRAV_H */
diff --git a/src/scheduler.c b/src/scheduler.c
index d1d343240b37f5afd5f41fecacf106b0e85f726f..887e8d54fd1fda941fe0024c45733226ce8d4b43 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -518,130 +518,138 @@ void scheduler_splittasks(struct scheduler *s) {
     } /* pair interaction? */
 
     /* Gravity interaction? */
-    else if (t->type == task_type_grav_mm) {
-
-      /* Get a handle on the cells involved. */
-      struct cell *ci = t->ci;
-      struct cell *cj = t->cj;
-
-      /* Self-interaction? */
-      if (cj == NULL) {
-
-        /* Ignore this task if the cell has no gparts. */
-        if (ci->gcount == 0) t->type = task_type_none;
-
-        /* If the cell is split, recurse. */
-        else if (ci->split) {
-
-          /* Make a single sub-task? */
-          if (scheduler_dosub && ci->gcount < space_subsize / ci->gcount) {
-
-            t->type = task_type_sub;
-            t->subtype = task_subtype_grav;
-
-          }
-
-          /* Otherwise, just split the task. */
-          else {
-
-            /* Split this task into tasks on its progeny. */
-            t->type = task_type_none;
-            for (int j = 0; j < 8; j++)
-              if (ci->progeny[j] != NULL && ci->progeny[j]->gcount > 0) {
-                if (t->type == task_type_none) {
-                  t->type = task_type_grav_mm;
-                  t->ci = ci->progeny[j];
-                  t->cj = NULL;
-                } else
-                  t = scheduler_addtask(s, task_type_grav_mm, task_subtype_none,
-                                        0, 0, ci->progeny[j], NULL, 0);
-                for (int k = j + 1; k < 8; k++)
-                  if (ci->progeny[k] != NULL && ci->progeny[k]->gcount > 0) {
-                    if (t->type == task_type_none) {
-                      t->type = task_type_grav_mm;
-                      t->ci = ci->progeny[j];
-                      t->cj = ci->progeny[k];
-                    } else
-                      t = scheduler_addtask(s, task_type_grav_mm,
-                                            task_subtype_none, 0, 0,
-                                            ci->progeny[j], ci->progeny[k], 0);
-                  }
-              }
-            redo = (t->type != task_type_none);
-          }
-
-        }
-
-        /* Otherwise, just make a pp task out of it. */
-        else
-          t->type = task_type_grav_pp;
-
-      }
-
-      /* Nope, pair. */
-      else {
-
-        /* Make a sub-task? */
-        if (scheduler_dosub && ci->gcount < space_subsize / cj->gcount) {
-
-          t->type = task_type_sub;
-          t->subtype = task_subtype_grav;
-
-        }
-
-        /* Otherwise, split the task. */
-        else {
-
-          /* Get the opening angle theta. */
-          float dx[3], theta;
-          for (int k = 0; k < 3; k++) {
-            dx[k] = fabs(ci->loc[k] - cj->loc[k]);
-            if (s->space->periodic && dx[k] > 0.5 * s->space->dim[k])
-              dx[k] = -dx[k] + s->space->dim[k];
-            if (dx[k] > 0.0f) dx[k] -= ci->h[k];
-          }
-          theta =
-              (dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]) /
-              (ci->h[0] * ci->h[0] + ci->h[1] * ci->h[1] + ci->h[2] * ci->h[2]);
-
-          /* Ignore this task if the cell has no gparts. */
-          if (ci->gcount == 0 || cj->gcount == 0) t->type = task_type_none;
-
-          /* Split the interaction? */
-          else if (theta < const_theta_max * const_theta_max) {
-
-            /* Are both ci and cj split? */
-            if (ci->split && cj->split) {
-
-              /* Split this task into tasks on its progeny. */
-              t->type = task_type_none;
-              for (int j = 0; j < 8; j++)
-                if (ci->progeny[j] != NULL && ci->progeny[j]->gcount > 0) {
-                  for (int k = 0; k < 8; k++)
-                    if (cj->progeny[k] != NULL && cj->progeny[k]->gcount > 0) {
-                      if (t->type == task_type_none) {
-                        t->type = task_type_grav_mm;
-                        t->ci = ci->progeny[j];
-                        t->cj = cj->progeny[k];
-                      } else
-                        t = scheduler_addtask(
-                            s, task_type_grav_mm, task_subtype_none, 0, 0,
-                            ci->progeny[j], cj->progeny[k], 0);
-                    }
-                }
-              redo = (t->type != task_type_none);
-
-            }
-
-            /* Otherwise, make a pp task out of it. */
-            else
-              t->type = task_type_grav_pp;
-          }
-        }
-
-      } /* gravity pair interaction? */
-
-    } /* gravity interaction? */
+    /* else if (t->type == task_type_grav_mm) { */
+
+    /*   /\* Get a handle on the cells involved. *\/ */
+    /*   struct cell *ci = t->ci; */
+    /*   struct cell *cj = t->cj; */
+
+    /*   /\* Self-interaction? *\/ */
+    /*   if (cj == NULL) { */
+
+    /*     /\* Ignore this task if the cell has no gparts. *\/ */
+    /*     if (ci->gcount == 0) t->type = task_type_none; */
+
+    /*     /\* If the cell is split, recurse. *\/ */
+    /*     else if (ci->split) { */
+
+    /*       /\* Make a single sub-task? *\/ */
+    /*       if (scheduler_dosub && ci->gcount < space_subsize / ci->gcount) {
+     */
+
+    /*         t->type = task_type_sub; */
+    /*         t->subtype = task_subtype_grav; */
+
+    /*       } */
+
+    /*       /\* Otherwise, just split the task. *\/ */
+    /*       else { */
+
+    /*         /\* Split this task into tasks on its progeny. *\/ */
+    /*         t->type = task_type_none; */
+    /*         for (int j = 0; j < 8; j++) */
+    /*           if (ci->progeny[j] != NULL && ci->progeny[j]->gcount > 0) { */
+    /*             if (t->type == task_type_none) { */
+    /*               t->type = task_type_grav_mm; */
+    /*               t->ci = ci->progeny[j]; */
+    /*               t->cj = NULL; */
+    /*             } else */
+    /*               t = scheduler_addtask(s, task_type_grav_mm,
+     * task_subtype_none, */
+    /*                                     0, 0, ci->progeny[j], NULL, 0); */
+    /*             for (int k = j + 1; k < 8; k++) */
+    /*               if (ci->progeny[k] != NULL && ci->progeny[k]->gcount > 0) {
+     */
+    /*                 if (t->type == task_type_none) { */
+    /*                   t->type = task_type_grav_mm; */
+    /*                   t->ci = ci->progeny[j]; */
+    /*                   t->cj = ci->progeny[k]; */
+    /*                 } else */
+    /*                   t = scheduler_addtask(s, task_type_grav_mm, */
+    /*                                         task_subtype_none, 0, 0, */
+    /*                                         ci->progeny[j], ci->progeny[k],
+     * 0); */
+    /*               } */
+    /*           } */
+    /*         redo = (t->type != task_type_none); */
+    /*       } */
+
+    /*     } */
+
+    /*     /\* Otherwise, just make a pp task out of it. *\/ */
+    /*     else */
+    /*       t->type = task_type_grav_pp; */
+
+    /*   } */
+
+    /*   /\* Nope, pair. *\/ */
+    /*   else { */
+
+    /*     /\* Make a sub-task? *\/ */
+    /*     if (scheduler_dosub && ci->gcount < space_subsize / cj->gcount) { */
+
+    /*       t->type = task_type_sub; */
+    /*       t->subtype = task_subtype_grav; */
+
+    /*     } */
+
+    /*     /\* Otherwise, split the task. *\/ */
+    /*     else { */
+
+    /*       /\* Get the opening angle theta. *\/ */
+    /*       float dx[3], theta; */
+    /*       for (int k = 0; k < 3; k++) { */
+    /*         dx[k] = fabs(ci->loc[k] - cj->loc[k]); */
+    /*         if (s->space->periodic && dx[k] > 0.5 * s->space->dim[k]) */
+    /*           dx[k] = -dx[k] + s->space->dim[k]; */
+    /*         if (dx[k] > 0.0f) dx[k] -= ci->h[k]; */
+    /*       } */
+    /*       theta = */
+    /*           (dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]) / */
+    /*           (ci->h[0] * ci->h[0] + ci->h[1] * ci->h[1] + ci->h[2] *
+     * ci->h[2]); */
+
+    /*       /\* Ignore this task if the cell has no gparts. *\/ */
+    /*       if (ci->gcount == 0 || cj->gcount == 0) t->type = task_type_none;
+     */
+
+    /*       /\* Split the interaction? *\/ */
+    /*       else if (theta < const_theta_max * const_theta_max) { */
+
+    /*         /\* Are both ci and cj split? *\/ */
+    /*         if (ci->split && cj->split) { */
+
+    /*           /\* Split this task into tasks on its progeny. *\/ */
+    /*           t->type = task_type_none; */
+    /*           for (int j = 0; j < 8; j++) */
+    /*             if (ci->progeny[j] != NULL && ci->progeny[j]->gcount > 0) {
+     */
+    /*               for (int k = 0; k < 8; k++) */
+    /*                 if (cj->progeny[k] != NULL && cj->progeny[k]->gcount > 0)
+     * { */
+    /*                   if (t->type == task_type_none) { */
+    /*                     t->type = task_type_grav_mm; */
+    /*                     t->ci = ci->progeny[j]; */
+    /*                     t->cj = cj->progeny[k]; */
+    /*                   } else */
+    /*                     t = scheduler_addtask( */
+    /*                         s, task_type_grav_mm, task_subtype_none, 0, 0, */
+    /*                         ci->progeny[j], cj->progeny[k], 0); */
+    /*                 } */
+    /*             } */
+    /*           redo = (t->type != task_type_none); */
+
+    /*         } */
+
+    /*         /\* Otherwise, make a pp task out of it. *\/ */
+    /*         else */
+    /*           t->type = task_type_grav_pp; */
+    /*       } */
+    /*     } */
+
+    /*   } /\* gravity pair interaction? *\/ */
+
+    /* } /\* gravity interaction? *\/ */
 
   } /* loop over all tasks. */
 }
diff --git a/src/space.c b/src/space.c
index 267eca1989a43caa53827695df24aa8c8f81db12..d87107a80add1b04674b74a8af55004427ddeac3 100644
--- a/src/space.c
+++ b/src/space.c
@@ -452,8 +452,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
       const int t = ind[k];
       ind[k] = ind[nr_parts];
       ind[nr_parts] = t;
-    }
-    else {
+    } else {
       /* Increment when not exchanging otherwise we need to retest "k".*/
       k++;
     }
@@ -488,8 +487,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
       const int t = gind[k];
       gind[k] = gind[nr_gparts];
       gind[nr_gparts] = t;
-    }
-    else {
+    } else {
       /* Increment when not exchanging otherwise we need to retest "k".*/
       k++;
     }
@@ -508,7 +506,6 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
     }
   }*/
 
-
   /* Exchange the strays, note that this potentially re-allocates
      the parts arrays. */
   size_t nr_parts_exchanged = s->nr_parts - nr_parts;
diff --git a/src/task.c b/src/task.c
index 31c2bfcb3462b9f91547ebd631f645e190c07143..77503d2207ae1a80e1c16f9eecc9fd59ed93f871 100644
--- a/src/task.c
+++ b/src/task.c
@@ -47,10 +47,10 @@
 
 /* Task type names. */
 const char *taskID_names[task_type_count] = {
-    "none",          "sort",      "self",       "pair",       "sub",
-    "init",          "ghost",     "drift",      "kick",       "send",
-    "recv",          "grav_pp",   "grav_mm",    "grav_up",    "grav_down",
-    "grav_external", "part_sort", "gpart_sort", "split_cell", "rewait"};
+    "none",       "sort",    "self",          "pair",      "sub",
+    "init",       "ghost",   "drift",         "kick",      "send",
+    "recv",       "grav_up", "grav_external", "part_sort", "gpart_sort",
+    "split_cell", "rewait"};
 
 const char *subtaskID_names[task_type_count] = {"none",  "density",
                                                 "force", "grav"};
@@ -124,12 +124,12 @@ void task_unlock(struct task *t) {
       cell_unlocktree(t->ci);
       if (t->cj != NULL) cell_unlocktree(t->cj);
       break;
-    case task_type_grav_pp:
-    case task_type_grav_mm:
-    case task_type_grav_down:
-      cell_gunlocktree(t->ci);
-      if (t->cj != NULL) cell_gunlocktree(t->cj);
-      break;
+    /*   case task_type_grav_pp: */
+    /*   case task_type_grav_mm: */
+    /*   case task_type_grav_down: */
+    /*     cell_gunlocktree(t->ci); */
+    /*     if (t->cj != NULL) cell_gunlocktree(t->cj); */
+    /*     break; */
     default:
       break;
   }
@@ -184,15 +184,15 @@ int task_lock(struct task *t) {
   }
 
   /* Gravity tasks? */
-  else if (type == task_type_grav_mm || type == task_type_grav_pp ||
-           type == task_type_grav_down) {
-    if (ci->ghold || (cj != NULL && cj->ghold)) return 0;
-    if (cell_glocktree(ci) != 0) return 0;
-    if (cj != NULL && cell_glocktree(cj) != 0) {
-      cell_gunlocktree(ci);
-      return 0;
-    }
-  }
+  /* else if (type == task_type_grav_mm || type == task_type_grav_pp || */
+  /*          type == task_type_grav_down) { */
+  /*   if (ci->ghold || (cj != NULL && cj->ghold)) return 0; */
+  /*   if (cell_glocktree(ci) != 0) return 0; */
+  /*   if (cj != NULL && cell_glocktree(cj) != 0) { */
+  /*     cell_gunlocktree(ci); */
+  /*     return 0; */
+  /*   } */
+  /* } */
 
   /* If we made it this far, we've got a lock. */
   return 1;
diff --git a/src/task.h b/src/task.h
index 25cc886f4b38456a0431fb6c7d0b7b1864053dd9..4956860c96020a211db5137c96ad4d260582e1c7 100644
--- a/src/task.h
+++ b/src/task.h
@@ -44,10 +44,10 @@ enum task_types {
   task_type_kick,
   task_type_send,
   task_type_recv,
-  task_type_grav_pp,
-  task_type_grav_mm,
+  /* task_type_grav_pp, */
+  /* task_type_grav_mm, */
   task_type_grav_up,
-  task_type_grav_down,
+  /* task_type_grav_down, */
   task_type_grav_external,
   task_type_part_sort,
   task_type_gpart_sort,