diff --git a/examples/ZeldovichPancake_3D/zeldovichPancake.yml b/examples/ZeldovichPancake_3D/zeldovichPancake.yml
index 39e02e3aab685bc4435d80e1d756b0b0c843c879..ca415706eefdcf679d0d5437dbf320b64334c997 100644
--- a/examples/ZeldovichPancake_3D/zeldovichPancake.yml
+++ b/examples/ZeldovichPancake_3D/zeldovichPancake.yml
@@ -46,6 +46,7 @@ Cosmology:
 Gravity:
   mesh_side_length:   16
   eta: 0.025
-  theta: 0.8
+  theta: 0.5
+  r_cut_max: 6.
   comoving_softening: 0.0001
   max_physical_softening: 0.0001
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 62f4f4a5be68924f8c1b159e0e20f50297190316..3a7b0e05f1187e2fce8f6a3753d4d0d010a2835a 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -149,6 +149,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
     const float h, const float h_inv, const struct multipole *m, float *f_x,
     float *f_y, float *f_z, float *pot) {
 
+  error("www");
+
 /* In the case where the order is < 3, then there is only a monopole term left.
  * We can default to the normal P-P interaction with the mass of the multipole
  * and its CoM as the "particle" property */
@@ -235,14 +237,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
 /* In the case where the order is < 3, then there is only a monopole term left.
  * We can default to the normal P-P interaction with the mass of the multipole
  * and its CoM as the "particle" property */
-#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
+#if 1  // SELF_GRAVITY_MULTIPOLE_ORDER < 3
 
   float f_ij, pot_ij;
   runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
                                 m->M_000, r_s_inv, &f_ij, &pot_ij);
-  *f_x = f_ij * r_x;
-  *f_y = f_ij * r_y;
-  *f_z = f_ij * r_z;
+  *f_x = -f_ij * r_x;
+  *f_y = -f_ij * r_y;
+  *f_z = -f_ij * r_z;
   *pot = pot_ij;
 
 #else
diff --git a/src/gravity_cache.h b/src/gravity_cache.h
index 5d2338127acabc076e27ccaee9c1bc8748dbd661..7eebc396b16b85205ebc8d8d423734f3813e781f 100644
--- a/src/gravity_cache.h
+++ b/src/gravity_cache.h
@@ -135,6 +135,35 @@ static INLINE void gravity_cache_init(struct gravity_cache *c,
   c->count = padded_count;
 }
 
+/**
+ * @param Zero all the output fields (acceleration and potential) of a
+ * #gravity_cache.
+ *
+ * @param c The #gravity_cache to zero.
+ * @param gcount_padded The padded size of the cache arrays.
+ */
+__attribute__((always_inline)) INLINE static void gravity_cache_zero_output(
+    struct gravity_cache *c, const int gcount_padded) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (gcount_padded % VEC_SIZE != 0)
+    error("Padded gcount size not a multiple of the vector length");
+#endif
+
+  /* Make the compiler understand we are in happy vectorization land */
+  swift_declare_aligned_ptr(float, a_x, c->a_x, SWIFT_CACHE_ALIGNMENT);
+  swift_declare_aligned_ptr(float, a_y, c->a_y, SWIFT_CACHE_ALIGNMENT);
+  swift_declare_aligned_ptr(float, a_z, c->a_z, SWIFT_CACHE_ALIGNMENT);
+  swift_declare_aligned_ptr(float, pot, c->pot, SWIFT_CACHE_ALIGNMENT);
+  swift_assume_size(gcount_padded, VEC_SIZE);
+
+  /* Zero everything */
+  bzero(a_x, gcount_padded * sizeof(float));
+  bzero(a_y, gcount_padded * sizeof(float));
+  bzero(a_z, gcount_padded * sizeof(float));
+  bzero(pot, gcount_padded * sizeof(float));
+}
+
 /**
  * @brief Fills a #gravity_cache structure with some #gpart and shift them.
  *
@@ -222,6 +251,9 @@ __attribute__((always_inline)) INLINE static void gravity_cache_populate(
     active[i] = 0;
     use_mpole[i] = 0;
   }
+
+  /* Zero the output as well */
+  gravity_cache_zero_output(c, gcount_padded);
 }
 
 /**
@@ -284,6 +316,9 @@ gravity_cache_populate_no_mpole(const timebin_t max_active_bin,
     m[i] = 0.f;
     active[i] = 0;
   }
+
+  /* Zero the output as well */
+  gravity_cache_zero_output(c, gcount_padded);
 }
 
 /**
@@ -301,12 +336,16 @@ gravity_cache_populate_no_mpole(const timebin_t max_active_bin,
  */
 __attribute__((always_inline)) INLINE static void
 gravity_cache_populate_all_mpole(const timebin_t max_active_bin,
+                                 const int periodic, const float dim[3],
                                  struct gravity_cache *c,
                                  const struct gpart *restrict gparts,
                                  const int gcount, const int gcount_padded,
-                                 const struct cell *cell,
+                                 const struct cell *cell, const float CoM[3],
+                                 const float r_max2,
                                  const struct gravity_props *grav_props) {
 
+  const float theta_crit2 = grav_props->theta_crit2;
+
   /* Make the compiler understand we are in happy vectorization land */
   swift_declare_aligned_ptr(float, x, c->x, SWIFT_CACHE_ALIGNMENT);
   swift_declare_aligned_ptr(float, y, c->y, SWIFT_CACHE_ALIGNMENT);
@@ -327,6 +366,24 @@ gravity_cache_populate_all_mpole(const timebin_t max_active_bin,
     m[i] = gparts[i].mass;
     active[i] = (int)(gparts[i].time_bin <= max_active_bin);
     use_mpole[i] = 1;
+
+#ifdef SWIFT_DEBUG_CHECKS
+    /* Distance to the CoM of the other cell. */
+    float dx = x[i] - CoM[0];
+    float dy = y[i] - CoM[1];
+    float dz = z[i] - CoM[2];
+
+    /* Apply periodic BC */
+    if (periodic) {
+      dx = nearestf(dx, dim[0]);
+      dy = nearestf(dy, dim[1]);
+      dz = nearestf(dz, dim[2]);
+    }
+    const float r2 = dx * dx + dy * dy + dz * dz;
+
+    if (!gravity_M2P_accept(r_max2, theta_crit2, r2))
+      error("Using m-pole where the test fails");
+#endif
   }
 
 #ifdef SWIFT_DEBUG_CHECKS
@@ -350,6 +407,9 @@ gravity_cache_populate_all_mpole(const timebin_t max_active_bin,
     active[i] = 0;
     use_mpole[i] = 0;
   }
+
+  /* Zero the output as well */
+  gravity_cache_zero_output(c, gcount_padded);
 }
 
 /**
diff --git a/src/multipole.h b/src/multipole.h
index d75704db2514ba6f1f175a47678a42a85a40eeaa..f85a2fe63d690c23841920cf56d06c13d153d04d 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -2427,7 +2427,7 @@ __attribute__((always_inline)) INLINE static int gravity_M2P_accept(
   // MATTHIEU: Make this mass-dependent ?
 
   /* Multipole acceptance criterion (Dehnen 2002, eq.10) */
-  return (r2 * theta_crit2 * 0.1 > r_max2);
+  return (r2 * theta_crit2 > r_max2);
 }
 
 #endif /* SWIFT_MULTIPOLE_H */
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 5109cc77a670481dcca37f080efe21c0a0c9bbf8..7e3933f185568444a4f545fe53aab481e51efa1b 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -92,7 +92,7 @@ static INLINE void runner_do_grav_down(struct runner *r, struct cell *c,
     /* Leaf case */
 
     /* We can abort early if no interactions via multipole happened */
-    if (!c->multipole->pot.interacted) return;
+    // if (!c->multipole->pot.interacted) return;
 
     if (!cell_are_gpart_drifted(c, e)) error("Un-drifted gparts");
 
@@ -242,10 +242,10 @@ static INLINE void runner_dopair_grav_pp_full(
     }
 
     /* Store everything back in cache */
-    ci_cache->a_x[pid] = a_x;
-    ci_cache->a_y[pid] = a_y;
-    ci_cache->a_z[pid] = a_z;
-    ci_cache->pot[pid] = pot;
+    ci_cache->a_x[pid] += a_x;
+    ci_cache->a_y[pid] += a_y;
+    ci_cache->a_z[pid] += a_z;
+    ci_cache->pot[pid] += pot;
   }
 }
 
@@ -368,10 +368,10 @@ static INLINE void runner_dopair_grav_pp_truncated(
     }
 
     /* Store everything back in cache */
-    ci_cache->a_x[pid] = a_x;
-    ci_cache->a_y[pid] = a_y;
-    ci_cache->a_z[pid] = a_z;
-    ci_cache->pot[pid] = pot;
+    ci_cache->a_x[pid] += a_x;
+    ci_cache->a_y[pid] += a_y;
+    ci_cache->a_z[pid] += a_z;
+    ci_cache->pot[pid] += pot;
   }
 }
 
@@ -474,10 +474,10 @@ static INLINE void runner_dopair_grav_pm_full(
                              &f_z, &pot_ij);
 
     /* Store it back */
-    a_x[pid] = f_x;
-    a_y[pid] = f_y;
-    a_z[pid] = f_z;
-    pot[pid] = pot_ij;
+    a_x[pid] += f_x;
+    a_y[pid] += f_y;
+    a_z[pid] += f_z;
+    pot[pid] += pot_ij;
 
 #ifdef SWIFT_DEBUG_CHECKS
     /* Update the interaction counter */
@@ -572,7 +572,7 @@ static INLINE void runner_dopair_grav_pm_truncated(
 
     const float r2 = dx * dx + dy * dy + dz * dz;
 
-#ifdef SWIFT_DEBUG_CHECKSa
+#ifdef SWIFT_DEBUG_CHECKS
     const float r_max_j = cj->multipole->r_max;
     const float r_max2 = r_max_j * r_max_j;
     const float theta_crit2 = e->gravity_properties->theta_crit2;
@@ -591,10 +591,10 @@ static INLINE void runner_dopair_grav_pm_truncated(
                                   multi_j, &f_x, &f_y, &f_z, &pot_ij);
 
     /* Store it back */
-    a_x[pid] = f_x;
-    a_y[pid] = f_y;
-    a_z[pid] = f_z;
-    pot[pid] = pot_ij;
+    a_x[pid] += f_x;
+    a_y[pid] += f_y;
+    a_z[pid] += f_z;
+    pot[pid] += pot_ij;
 
 #ifdef SWIFT_DEBUG_CHECKS
     /* Update the interaction counter */
@@ -692,11 +692,11 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
 
   /* Fill the caches */
   gravity_cache_populate(e->max_active_bin, periodic, dim, ci_cache, ci->gparts,
-                         gcount_i, gcount_padded_i, shift_i, CoM_j, rmax2_j, ci,
-                         e->gravity_properties);
+                         gcount_i, gcount_padded_i, shift_i, CoM_j,
+                         rmax2_j, ci, e->gravity_properties);
   gravity_cache_populate(e->max_active_bin, periodic, dim, cj_cache, cj->gparts,
-                         gcount_j, gcount_padded_j, shift_j, CoM_i, rmax2_i, cj,
-                         e->gravity_properties);
+                         gcount_j, gcount_padded_j, shift_j, CoM_i,
+                         rmax2_i, cj, e->gravity_properties);
 
   /* Can we use the Newtonian version or do we need the truncated one ? */
   if (!periodic) {
@@ -751,9 +751,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
                                         ci->gparts, cj->gparts);
 
         /* Then the M2P */
-        runner_dopair_grav_pm_truncated(ci_cache, gcount_padded_i, CoM_j,
-                                        multi_j, dim, r_s_inv, e, ci->gparts,
-                                        gcount_i, cj);
+        if (0)
+          runner_dopair_grav_pm_truncated(ci_cache, gcount_padded_i, CoM_j,
+                                          multi_j, dim, r_s_inv, e, ci->gparts,
+                                          gcount_i, cj);
       }
       if (cj_active && symmetric) {
 
@@ -763,9 +764,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
                                         cj->gparts, ci->gparts);
 
         /* Then the M2P */
-        runner_dopair_grav_pm_truncated(cj_cache, gcount_padded_j, CoM_i,
-                                        multi_i, dim, r_s_inv, e, cj->gparts,
-                                        gcount_j, ci);
+        if (0)
+          runner_dopair_grav_pm_truncated(cj_cache, gcount_padded_j, CoM_i,
+                                          multi_i, dim, r_s_inv, e, cj->gparts,
+                                          gcount_j, ci);
       }
 
     } else {
@@ -781,8 +783,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
                                    ci->gparts, cj->gparts);
 
         /* Then the M2P */
-        runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j,
-                                   periodic, dim, e, ci->gparts, gcount_i, cj);
+        if (0)
+          runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j,
+                                     periodic, dim, e, ci->gparts, gcount_i,
+                                     cj);
       }
       if (cj_active && symmetric) {
 
@@ -792,8 +796,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
                                    cj->gparts, ci->gparts);
 
         /* Then the M2P */
-        runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i,
-                                   periodic, dim, e, cj->gparts, gcount_j, ci);
+        if (0)
+          runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i,
+                                     periodic, dim, e, cj->gparts, gcount_j,
+                                     ci);
       }
     }
   }
@@ -899,10 +905,10 @@ static INLINE void runner_doself_grav_pp_full(
     }
 
     /* Store everything back in cache */
-    ci_cache->a_x[pid] = a_x;
-    ci_cache->a_y[pid] = a_y;
-    ci_cache->a_z[pid] = a_z;
-    ci_cache->pot[pid] = pot;
+    ci_cache->a_x[pid] += a_x;
+    ci_cache->a_y[pid] += a_y;
+    ci_cache->a_z[pid] += a_z;
+    ci_cache->pot[pid] += pot;
   }
 }
 
@@ -1009,10 +1015,10 @@ static INLINE void runner_doself_grav_pp_truncated(
     }
 
     /* Store everything back in cache */
-    ci_cache->a_x[pid] = a_x;
-    ci_cache->a_y[pid] = a_y;
-    ci_cache->a_z[pid] = a_z;
-    ci_cache->pot[pid] = pot;
+    ci_cache->a_x[pid] += a_x;
+    ci_cache->a_y[pid] += a_y;
+    ci_cache->a_z[pid] += a_z;
+    ci_cache->pot[pid] += pot;
   }
 }
 
@@ -1209,17 +1215,18 @@ static INLINE void runner_dopair_recursive_grav_pm(struct runner *r,
       error("Not enough space in the cache! gcount_i=%d", gcount_i);
 #endif
 
-    /* Fill the cache */
-    gravity_cache_populate_all_mpole(e->max_active_bin, ci_cache, ci->gparts,
-                                     gcount_i, gcount_padded_i, ci,
-                                     e->gravity_properties);
-
     /* Recover the multipole info and the CoM locations */
     const struct multipole *multi_j = &cj->multipole->m_pole;
+    const float r_max = cj->multipole->r_max;
     const float CoM_j[3] = {(float)(cj->multipole->CoM[0]),
                             (float)(cj->multipole->CoM[1]),
                             (float)(cj->multipole->CoM[2])};
 
+    /* Fill the cache */
+    gravity_cache_populate_all_mpole(
+        e->max_active_bin, periodic, dim, ci_cache, ci->gparts, gcount_i,
+        gcount_padded_i, ci, CoM_j, r_max * r_max, e->gravity_properties);
+
     /* Can we use the Newtonian version or do we need the truncated one ? */
     if (!periodic) {
 
@@ -1540,9 +1547,13 @@ static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
         continue;
       }
 
-      /* Call the PM interaction fucntion on the active sub-cells of ci */
       ++cc;
-      runner_dopair_recursive_grav_pm(r, ci, cj);
+
+      /* Call the PM interaction fucntion on the active sub-cells of ci */
+
+      // runner_dopair_recursive_grav_pm(r, ci, cj);
+      runner_dopair_grav_mm(r, ci, cj);
+      // runner_dopair_grav_pp(r, ci, cj, 0);
 
     } /* We are in charge of this pair */
   }   /* Loop over top-level cells */