diff --git a/src/engine.c b/src/engine.c
index 93d35ef62d3f31ac20de771f3ab635ac5e0eee73..716ede5b1bd9101417f8305ebd654e712c65a305 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -795,99 +795,16 @@ void engine_make_hydroloop_tasks(struct engine *e) {
 }
 
 /**
- * @brief Constructs the top-level pair tasks for the gravity M-M interactions
+ * @brief Counts the tasks associated with one cell and constructs the links
  *
  * @param e The #engine.
  */
-void engine_make_gravityinteraction_tasks(struct engine *e) {
+void engine_count_and_link_tasks(struct engine *e) {
 
-  struct space *s = e->s;
   struct scheduler *sched = &e->sched;
-  const int nr_cells = s->nr_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++) {
+  const int nr_tasks = sched->nr_tasks;
 
-        /* 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);
-        }
-      }
-    }
-  }
-}
-
-/**
- * @brief Fill the #space's task list.
- *
- * @param e The #engine we are working with.
- */
-
-void engine_maketasks(struct engine *e) {
-
-  struct space *s = e->s;
-  struct scheduler *sched = &e->sched;
-  struct cell *cells = s->cells;
-  const int nr_cells = s->nr_cells;
-  const int nodeID = e->nodeID;
-  const ticks tic = getticks();
-
-  /* Re-set the scheduler. */
-  scheduler_reset(sched, s->tot_cells * engine_maxtaskspercell);
-
-  /* Add the space sorting tasks. */
-  for (int i = 0; i < e->nr_threads; i++)
-    scheduler_addtask(sched, task_type_psort, task_subtype_none, i, 0, NULL,
-                      NULL, 0);
-
-  /* Construct the firt hydro loop over neighbours */
-  engine_make_hydroloop_tasks(e);
-
-  /* Add the gravity mm tasks. */
-  engine_make_gravityinteraction_tasks(e);
-
-  /* Split the tasks. */
-  scheduler_splittasks(sched);
-
-  /* Allocate the list of cell-task links. The maximum number of links
-     is the number of cells (s->tot_cells) times the number of neighbours (27)
-     times the number of interaction types (2, density and force). */
-  if (e->links != NULL) free(e->links);
-  e->size_links = s->tot_cells * 27 * 2;
-  if ((e->links = malloc(sizeof(struct link) * e->size_links)) == NULL)
-    error("Failed to allocate cell-task links.");
-  e->nr_links = 0;
-
-  /* Add the gravity up/down tasks at the top-level cells and push them down. */
-  for (int k = 0; k < nr_cells; k++)
-    if (cells[k].nodeID == nodeID && cells[k].gcount > 0) {
-
-      /* Create tasks at top level. */
-      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);
-
-      /* Push tasks down the cell hierarchy. */
-      engine_addtasks_grav(e, &cells[k], up, down);
-    }
-
-  /* Count the number of tasks associated with each cell and
-     store the density tasks in each cell, and make each sort
-     depend on the sorts of its progeny. */
-  for (int k = 0; k < sched->nr_tasks; k++) {
+  for (int k = 0; k < nr_tasks; k++) {
 
     /* Get the current task. */
     struct task *t = &sched->tasks[k];
@@ -942,16 +859,21 @@ void engine_maketasks(struct engine *e) {
       }
     }
   }
+}
 
-  /* Append a ghost task to each cell, and add kick tasks to the
-     super cells. */
-  for (int k = 0; k < nr_cells; k++) engine_mkghosts(e, &cells[k], NULL);
+/**
+ * @brief Duplicates the first hydro loop and creates the corresponding
+ *dependencies using the ghost tasks.
+ *
+ * @parma e The #engine.
+ */
+void engine_make_extra_hydroloop_tasks(struct engine *e) {
 
-  /* Run through the tasks and make force tasks for each density task.
-     Each force task depends on the cell ghosts and unlocks the kick task
-     of its super-cell. */
-  int sched_nr_tasks = sched->nr_tasks;
-  for (int k = 0; k < sched_nr_tasks; k++) {
+  struct scheduler *sched = &e->sched;
+  const int nodeID = e->nodeID;
+  const int nr_tasks = sched->nr_tasks;
+
+  for (int k = 0; k < nr_tasks; k++) {
 
     /* Get a pointer to the task. */
     struct task *t = &sched->tasks[k];
@@ -1021,28 +943,153 @@ void engine_maketasks(struct engine *e) {
     /* else if (t->type == task_type_kick && t->ci->grav_down != NULL) */
     /*   scheduler_addunlock(sched, t->ci->grav_down, t); */
   }
+}
 
-/* Add the communication tasks if MPI is being used. */
-#ifdef WITH_MPI
+/**
+ * @brief Constructs the top-level pair tasks for the gravity M-M interactions
+ *
+ * @param e The #engine.
+ */
+void engine_make_gravityinteraction_tasks(struct engine *e) {
 
-  /* Loop over the proxies. */
-  for (int pid = 0; pid < e->nr_proxies; pid++) {
+  struct space *s = e->s;
+  struct scheduler *sched = &e->sched;
+  const int nr_cells = s->nr_cells;
+  struct cell *cells = s->cells;
+
+  /* Loop over all cells. */
+  for (int i = 0; i < nr_cells; i++) {
 
-    /* Get a handle on the proxy. */
-    struct proxy *p = &e->proxies[pid];
+    /* 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 through the proxy's incoming cells and add the
-       recv tasks. */
-    for (int k = 0; k < p->nr_cells_in; k++)
-      engine_addtasks_recv(e, p->cells_in[k], NULL, NULL);
+      /* Loop over all remainding cells */
+      for (int j = i + 1; j < nr_cells; j++) {
 
-    /* Loop through the proxy's outgoing cells and add the
-       send tasks. */
-    for (int k = 0; k < p->nr_cells_out; k++)
-      engine_addtasks_send(e, p->cells_out[k], p->cells_in[0]);
+        /* 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);
+        }
+      }
+    }
   }
+}
 
-#endif
+/**
+ * @brief Constructs the gravity tasks building the multipoles and propagating
+ *them to the children
+ *
+ * @param e The #engine.
+ */
+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;
+
+  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) {
+
+      /* Create tasks at top level. */
+      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);
+
+      /* Push tasks down the cell hierarchy. */
+      engine_addtasks_grav(e, &cells[k], up, down);
+    }
+  }
+}
+
+/**
+ * @brief Fill the #space's task list.
+ *
+ * @param e The #engine we are working with.
+ */
+
+void engine_maketasks(struct engine *e) {
+
+  struct space *s = e->s;
+  struct scheduler *sched = &e->sched;
+  struct cell *cells = s->cells;
+  const int nr_cells = s->nr_cells;
+  const ticks tic = getticks();
+
+  /* Re-set the scheduler. */
+  scheduler_reset(sched, s->tot_cells * engine_maxtaskspercell);
+
+  /* Add the space sorting tasks. */
+  for (int i = 0; i < e->nr_threads; i++)
+    scheduler_addtask(sched, task_type_psort, task_subtype_none, i, 0, NULL,
+                      NULL, 0);
+
+  /* Construct the firt hydro loop over neighbours */
+  engine_make_hydroloop_tasks(e);
+
+  /* Add the gravity mm tasks. */
+  if ((e->policy & engine_policy_self_gravity) == engine_policy_self_gravity)
+    engine_make_gravityinteraction_tasks(e);
+
+  /* Split the tasks. */
+  scheduler_splittasks(sched);
+
+  /* Allocate the list of cell-task links. The maximum number of links
+     is the number of cells (s->tot_cells) times the number of neighbours (27)
+     times the number of interaction types (2, density and force). */
+  if (e->links != NULL) free(e->links);
+  e->size_links = s->tot_cells * 27 * 2;
+  if ((e->links = malloc(sizeof(struct link) * e->size_links)) == NULL)
+    error("Failed to allocate cell-task links.");
+  e->nr_links = 0;
+
+  /* Add the gravity up/down tasks at the top-level cells and push them down. */
+  if ((e->policy & engine_policy_self_gravity) == engine_policy_self_gravity)
+    engine_make_gravityrecursive_tasks(e);
+
+  /* Count the number of tasks associated with each cell and
+     store the density tasks in each cell, and make each sort
+     depend on the sorts of its progeny. */
+  engine_count_and_link_tasks(e);
+
+  /* Append a ghost task to each cell, and add kick tasks to the
+     super cells. */
+  for (int k = 0; k < nr_cells; k++) engine_mkghosts(e, &cells[k], NULL);
+
+  /* Run through the tasks and make force tasks for each density task.
+     Each force task depends on the cell ghosts and unlocks the kick task
+     of its super-cell. */
+  engine_make_extra_hydroloop_tasks(e);
+
+  /* Add the communication tasks if MPI is being used. */
+  if ((e->policy & engine_policy_mpi) == engine_policy_mpi) {
+
+    /* Loop over the proxies. */
+    for (int pid = 0; pid < e->nr_proxies; pid++) {
+
+      /* Get a handle on the proxy. */
+      struct proxy *p = &e->proxies[pid];
+
+      /* Loop through the proxy's incoming cells and add the
+         recv tasks. */
+      for (int k = 0; k < p->nr_cells_in; k++)
+        engine_addtasks_recv(e, p->cells_in[k], NULL, NULL);
+
+      /* Loop through the proxy's outgoing cells and add the
+         send tasks. */
+      for (int k = 0; k < p->nr_cells_out; k++)
+        engine_addtasks_send(e, p->cells_out[k], p->cells_in[0]);
+    }
+  }
 
   /* Set the unlocks per task. */
   scheduler_set_unlocks(sched);
diff --git a/src/engine.h b/src/engine.h
index cb94da34dfc0c59c5c2b2d6da2644e41337906eb..80d836422d7d43f15c6cfd2b1b0924445e72fbb7 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -181,7 +181,10 @@ void engine_print(struct engine *e);
 void engine_init_particles(struct engine *e);
 void engine_step(struct engine *e);
 void engine_make_hydroloop_tasks(struct engine *e);
+void engine_make_extra_hydroloop_tasks(struct engine *e);
+void engine_count_and_link_tasks(struct engine *e);
 void engine_make_gravityinteraction_tasks(struct engine *e);
+void engine_make_gravityrecursive_tasks(struct engine *e);
 void engine_maketasks(struct engine *e);
 void engine_split(struct engine *e, struct partition *initial_partition);
 int engine_exchange_strays(struct engine *e, int offset, size_t *ind, size_t N);