diff --git a/src/engine.c b/src/engine.c
index facffa389e98fdeeebda222f04b2953c00b97373..4ca5ea512752bb83cd88042ad7aeda986000c938 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2980,7 +2980,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
                     MPI_COMM_WORLD) != MPI_SUCCESS)
     error("Failed to all-reduce total mass in the system.");
 #endif
-  message("Total mass in the system: %e", e->s->total_mass);
+  if (e->nodeID == 0) message("Total mass in the system: %e", e->s->total_mass);
 #endif
 
   /* Now time to get ready for the first time-step */
@@ -3079,7 +3079,7 @@ void engine_step(struct engine *e) {
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   /* Check the accuracy of the gravity calculation */
-  gravity_exact_force_check(e->s, e, 1e-1);
+  gravity_exact_force_check(e->s, e, 1e1);
 #endif
 
 /* Collect the values of rebuild from all nodes. */
diff --git a/src/gravity.c b/src/gravity.c
index 86f9fa82e3eb693cc3c051420fb9c7bff277eb9f..5b89dd398fbb380777fa9ab97fdc6c5d22d8242b 100644
--- a/src/gravity.c
+++ b/src/gravity.c
@@ -149,9 +149,12 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
       /* Check that we are not above tolerance */
       if (fabsf(err_rel[0]) > rel_tol || fabsf(err_rel[1]) > rel_tol ||
           fabsf(err_rel[2]) > rel_tol)
-        error("Error too large ! gp->a_grav=[%e %e %e] gp->a_exact=[%e %e %e]",
-              gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2],
-              gpi->a_grav_exact[0], gpi->a_grav_exact[1], gpi->a_grav_exact[2]);
+        error(
+            "Error too large ! gp->a_grav=[%e %e %e] gp->a_exact=[%e %e %e], "
+            "gp->mass_interacted=%e",
+            gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2],
+            gpi->a_grav_exact[0], gpi->a_grav_exact[1], gpi->a_grav_exact[2],
+            gpi->mass_interacted);
 
       /* Construct some statistics */
       for (int k = 0; k < 3; ++k) {
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 93a65e2f5a70e09ad14280cf9c334753359fb8b5..9ff69357163ca46a3d4b21d048cac465771e047e 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -113,7 +113,7 @@ __attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
     struct gpart* gp) {
 
   gp->time_bin = 0;
-  gp->epsilon = 0.;  // MATTHIEU
+  gp->epsilon = 0.1;  // MATTHIEU
 
   gravity_init_gpart(gp);
 }
diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h
index 6952681999f833bce7755a72aaee742a7fa0ed22..7b1c5984647c3be232770dc32fc1b112ad8bee94 100644
--- a/src/kernel_long_gravity.h
+++ b/src/kernel_long_gravity.h
@@ -37,14 +37,15 @@
 __attribute__((always_inline)) INLINE static void kernel_long_grav_eval(
     float u, float *const W) {
 
-  const float arg1 = u * 0.5f;
-  const float arg2 = u * one_over_sqrt_pi;
-  const float arg3 = -arg1 * arg1;
+  /* const float arg1 = u * 0.5f; */
+  /* const float arg2 = u * one_over_sqrt_pi; */
+  /* const float arg3 = -arg1 * arg1; */
 
-  const float term1 = erfcf(arg1);
-  const float term2 = arg2 * expf(arg3);
+  /* const float term1 = erfcf(arg1); */
+  /* const float term2 = arg2 * expf(arg3); */
 
-  *W = term1 + term2;
+  /* *W = term1 + term2; */
+  *W = 1.f;
 }
 
 #endif  // SWIFT_KERNEL_LONG_GRAVITY_H
diff --git a/src/multipole.h b/src/multipole.h
index 95b7450cd9c7aa60dba5c43ebd3c77e6316845a5..d5ede8866a54016b1ef2f5119dbfd3b868b123ee 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -104,7 +104,7 @@ INLINE static void gravity_reset(struct gravity_tensors *m) {
   bzero(m, sizeof(struct gravity_tensors));
 }
 
-INLINE static void gravity_field_tensor_init(struct gravity_tensors *m) {
+INLINE static void gravity_field_tensors_init(struct gravity_tensors *m) {
 
   bzero(&m->a_x, sizeof(struct acc_tensor));
   bzero(&m->a_y, sizeof(struct acc_tensor));
@@ -116,6 +116,19 @@ INLINE static void gravity_field_tensor_init(struct gravity_tensors *m) {
 #endif
 }
 
+/**
+ * @brief Adds field tensrs to other ones (i.e. does la += lb).
+ *
+ * @param la The gravity tensors to add to.
+ * @param lb The gravity tensors to add.
+ */
+INLINE static void gravity_field_tensors_add(struct gravity_tensors *la,
+                                             const struct gravity_tensors *lb) {
+  la->a_x.F_000 += lb->a_x.F_000;
+  la->a_y.F_000 += lb->a_y.F_000;
+  la->a_z.F_000 += lb->a_z.F_000;
+}
+
 /**
  * @brief Prints the content of a #multipole to stdout.
  *
@@ -322,9 +335,9 @@ INLINE static void gravity_M2L(struct gravity_tensors *l_a,
   const double r_inv = 1. / sqrt(r2);
 
   /* 1st order multipole term */
-  l_a->a_x.F_000 = D_100(dx, dy, dz, r_inv) * m_b->mass;
-  l_a->a_y.F_000 = D_010(dx, dy, dz, r_inv) * m_b->mass;
-  l_a->a_z.F_000 = D_001(dx, dy, dz, r_inv) * m_b->mass;
+  l_a->a_x.F_000 += D_100(dx, dy, dz, r_inv) * m_b->mass;
+  l_a->a_y.F_000 += D_010(dx, dy, dz, r_inv) * m_b->mass;
+  l_a->a_z.F_000 += D_001(dx, dy, dz, r_inv) * m_b->mass;
 
 #ifdef SWIFT_DEBUG_CHECKS
   l_a->mass_interacted += m_b->mass;
@@ -345,123 +358,28 @@ INLINE static void gravity_M2L(struct gravity_tensors *l_a,
 INLINE static void gravity_L2L(struct gravity_tensors *l_a,
                                const struct gravity_tensors *l_b,
                                const double pos_a[3], const double pos_b[3],
-                               int periodic) {}
+                               int periodic) {
+  error("Not implemented yet");
+}
 
 /**
- * @brief Applies the  #acc_tensor to a set of #gpart.
+ * @brief Applies the  #acc_tensor to a  #gpart.
  *
  * Corresponds to equation (28a).
  */
 INLINE static void gravity_L2P(const struct gravity_tensors *l,
-                               struct gpart *gparts, int gcount) {
-
-  for (int i = 0; i < gcount; ++i) {
+                               struct gpart *gp) {
 
 #ifdef SWIFT_DEBUG_CHECKS
-    struct gpart *gp = &gparts[i];
-
-// if(gpart_is_active(gp, e)){
-
-    gp->mass_interacted += l->mass_interacted;
+  gp->mass_interacted += l->mass_interacted;
 #endif
-    //}
-  }
-}
-
-#if 0
-
-/* Multipole function prototypes. */
-void multipole_add(struct gravity_tensors *m_sum, const struct gravity_tensors *m_term);
-void multipole_init(struct gravity_tensors *m, const struct gpart *gparts,
-                    int gcount);
-void multipole_reset(struct gravity_tensors *m);
 
-/* static void multipole_iact_mm(struct multipole *ma, struct multipole *mb, */
-/*                               double *shift); */
-/* void multipole_addpart(struct multipole *m, struct gpart *p); */
-/* void multipole_addparts(struct multipole *m, struct gpart *p, int N); */
+  // message("a=[%e %e %e]", l->a_x.F_000, l->a_y.F_000, l->a_z.F_000);
 
-/**
- * @brief Compute the pairwise interaction between two multipoles.
- *
- * @param ma The first #multipole.
- * @param mb The second #multipole.
- * @param shift The periodicity correction.
- */
-__attribute__((always_inline)) INLINE static void multipole_iact_mm(
-    struct gravity_tensors *ma, struct gravity_tensors *mb, double *shift) {
-  /*   float dx[3], ir, r, r2 = 0.0f, acc; */
-  /*   int k; */
-
-  /*   /\* Compute the multipole distance. *\/ */
-  /*   for (k = 0; k < 3; k++) { */
-  /*     dx[k] = ma->x[k] - mb->x[k] - shift[k]; */
-  /*     r2 += dx[k] * dx[k]; */
-  /*   } */
-
-  /*   /\* Compute the normalized distance vector. *\/ */
-  /*   ir = 1.0f / sqrtf(r2); */
-  /*   r = r2 * ir; */
-
-  /*   /\* Evaluate the gravity kernel. *\/ */
-  /*   kernel_grav_eval(r, &acc); */
-
-  /*   /\* Scale the acceleration. *\/ */
-  /*   acc *= const_G * ir * ir * ir; */
-
-  /* /\* Compute the forces on both multipoles. *\/ */
-  /* #if const_gravity_multipole_order == 1 */
-  /*   float mma = ma->coeffs[0], mmb = mb->coeffs[0]; */
-  /*   for (k = 0; k < 3; k++) { */
-  /*     ma->a[k] -= dx[k] * acc * mmb; */
-  /*     mb->a[k] += dx[k] * acc * mma; */
-  /*   } */
-  /* #else */
-  /* #error( "Multipoles of order %i not yet implemented." ,
-   * const_gravity_multipole_order )
-   */
-  /* #endif */
+  /* 0th order interaction */
+  gp->a_grav[0] += l->a_x.F_000;
+  gp->a_grav[1] += l->a_y.F_000;
+  gp->a_grav[2] += l->a_z.F_000;
 }
 
-/**
- * @brief Compute the interaction of a multipole on a particle.
- *
- * @param m The #multipole.
- * @param p The #gpart.
- * @param shift The periodicity correction.
- */
-__attribute__((always_inline)) INLINE static void multipole_iact_mp(
-    struct gravity_tensors *m, struct gpart *p, double *shift) {
-
-  /*   float dx[3], ir, r, r2 = 0.0f, acc; */
-  /*   int k; */
-
-  /*   /\* Compute the multipole distance. *\/ */
-  /*   for (k = 0; k < 3; k++) { */
-  /*     dx[k] = m->x[k] - p->x[k] - shift[k]; */
-  /*     r2 += dx[k] * dx[k]; */
-  /*   } */
-
-  /*   /\* Compute the normalized distance vector. *\/ */
-  /*   ir = 1.0f / sqrtf(r2); */
-  /*   r = r2 * ir; */
-
-  /*   /\* Evaluate the gravity kernel. *\/ */
-  /*   kernel_grav_eval(r, &acc); */
-
-  /*   /\* Scale the acceleration. *\/ */
-  /*   acc *= const_G * ir * ir * ir * m->coeffs[0]; */
-
-  /* /\* Compute the forces on both multipoles. *\/ */
-  /* #if const_gravity_multipole_order == 1 */
-  /*   for (k = 0; k < 3; k++) p->a_grav[k] += dx[k] * acc; */
-  /* #else */
-  /* #error( "Multipoles of order %i not yet implemented." ,
-   * const_gravity_multipole_order )
-   */
-  /* #endif */
-}
-
-#endif
-
 #endif /* SWIFT_MULTIPOLE_H */
diff --git a/src/runner.c b/src/runner.c
index 521914bfdecb18f6c5e4e0a55ec2f0e759424797..68ad4e7966359c229c7baeb469bbb240edb3b63b 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -497,7 +497,7 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
 
   /* Reset the gravity acceleration tensors */
   if (e->policy & engine_policy_self_gravity)
-    gravity_field_tensor_init(c->multipole);
+    gravity_field_tensors_init(c->multipole);
 
   /* Recurse? */
   if (c->split) {
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 9988b7d553a962ee541f20cab64cc534c991957a..71a11f38d4721dffbe584f8bae2dfc49d065249e 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -64,12 +64,22 @@ void runner_do_grav_down(struct runner *r, struct cell *c) {
       if (cp != NULL) {
         gravity_L2L(&temp, c->multipole, cp->multipole->CoM, c->multipole->CoM,
                     1);
+
+        gravity_field_tensors_add(cp->multipole, &temp);
       }
     }
 
   } else {
 
-    gravity_L2P(c->multipole, c->gparts, c->gcount);
+    const struct engine *e = r->e;
+    struct gpart *gparts = c->gparts;
+    const int gcount = c->gcount;
+
+    /* Apply accelerations to the particles */
+    for (int i = 0; i < gcount; ++i) {
+      struct gpart *gp = &gparts[i];
+      if (gpart_is_active(gp, e)) gravity_L2P(c->multipole, gp);
+    }
   }
 }
 
@@ -237,49 +247,55 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
     /* Get a hold of the ith part in ci. */
     struct gpart *restrict gpi = &gparts_i[pid];
 
+    if (!gpart_is_active(gpi, e)) continue;
+
     /* Loop over every particle in the other cell. */
     for (int pjd = 0; pjd < gcount_j; pjd++) {
 
       /* Get a hold of the jth part in cj. */
-      struct gpart *restrict gpj = &gparts_j[pjd];
+      const struct gpart *restrict gpj = &gparts_j[pjd];
 
       /* Compute the pairwise distance. */
-      float dx[3] = {gpi->x[0] - gpj->x[0],   // x
-                     gpi->x[1] - gpj->x[1],   // y
-                     gpi->x[2] - gpj->x[2]};  // z
+      const float dx[3] = {gpi->x[0] - gpj->x[0],   // x
+                           gpi->x[1] - gpj->x[1],   // y
+                           gpi->x[2] - gpj->x[2]};  // z
       const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
 
       /* Interact ! */
-      if (gpart_is_active(gpi, e) && gpart_is_active(gpj, e)) {
-
-        runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
+      runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
 
 #ifdef SWIFT_DEBUG_CHECKS
-        gpi->mass_interacted += gpj->mass;
-        gpj->mass_interacted += gpi->mass;
+      gpi->mass_interacted += gpj->mass;
 #endif
+    }
+  }
 
-      } else {
+  /* Loop over all particles in cj... */
+  for (int pjd = 0; pjd < gcount_j; pjd++) {
 
-        if (gpart_is_active(gpi, e)) {
+    /* Get a hold of the ith part in ci. */
+    struct gpart *restrict gpj = &gparts_j[pjd];
 
-          runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
+    if (!gpart_is_active(gpj, e)) continue;
 
-#ifdef SWIFT_DEBUG_CHECKS
-          gpi->mass_interacted += gpj->mass;
-#endif
-        } else if (gpart_is_active(gpj, e)) {
+    /* Loop over every particle in the other cell. */
+    for (int pid = 0; pid < gcount_i; pid++) {
 
-          dx[0] = -dx[0];
-          dx[1] = -dx[1];
-          dx[2] = -dx[2];
-          runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
+      /* Get a hold of the ith part in ci. */
+      const struct gpart *restrict gpi = &gparts_i[pid];
+
+      /* Compute the pairwise distance. */
+      const float dx[3] = {gpj->x[0] - gpi->x[0],   // x
+                           gpj->x[1] - gpi->x[1],   // y
+                           gpj->x[2] - gpi->x[2]};  // z
+      const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+      /* Interact ! */
+      runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
 
 #ifdef SWIFT_DEBUG_CHECKS
-          gpj->mass_interacted += gpi->mass;
+      gpj->mass_interacted += gpi->mass;
 #endif
-        }
-      }
     }
   }
 
@@ -468,9 +484,9 @@ static void runner_dopair_grav(struct runner *r, struct cell *ci,
 
             } else {
 
-              /* Ok, here we can go for particle-multipole interactions */
-              runner_dopair_grav_pm(r, ci->progeny[j], cj->progeny[k]);
-              runner_dopair_grav_pm(r, cj->progeny[k], ci->progeny[j]);
+              /* Ok, here we can go for multipole-multipole interactions */
+              runner_dopair_grav_mm(r, ci->progeny[j], cj->progeny[k]);
+              runner_dopair_grav_mm(r, cj->progeny[k], ci->progeny[j]);
             }
           }
         }