diff --git a/src/multipole.h b/src/multipole.h
index 8e8ecae54e5dedde1a789d2a5469b6e39d853b18..a0728e47b805c0311740cb74a3a43ef829ef26f3 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -1498,23 +1498,28 @@ INLINE static void gravity_M2M(struct multipole *m_a,
  * @param pos_a The position of the multipole.
  * @param props The #gravity_props of this calculation.
  * @param periodic Is the calculation periodic ?
+ * @param dim The size of the simulation box.
  */
 INLINE static void gravity_M2L(struct grav_tensor *l_b,
                                const struct multipole *m_a,
                                const double pos_b[3], const double pos_a[3],
-                               const struct gravity_props *props,
-                               int periodic) {
+                               const struct gravity_props *props, int periodic,
+                               const double dim[3]) {
 
   /* Recover some constants */
   const double eps2 = props->epsilon2;
 
   /* Compute distance vector */
-  const double dx =
-      periodic ? box_wrap(pos_b[0] - pos_a[0], 0., 1.) : pos_b[0] - pos_a[0];
-  const double dy =
-      periodic ? box_wrap(pos_b[1] - pos_a[1], 0., 1.) : pos_b[1] - pos_a[1];
-  const double dz =
-      periodic ? box_wrap(pos_b[2] - pos_a[2], 0., 1.) : pos_b[2] - pos_a[2];
+  double dx = pos_b[0] - pos_a[0];
+  double dy = pos_b[1] - pos_a[1];
+  double dz = pos_b[2] - pos_a[2];
+
+  /* Apply BC */
+  if (periodic) {
+    dx = nearest(dx, dim[0]);
+    dy = nearest(dy, dim[1]);
+    dz = nearest(dz, dim[2]);
+  }
 
   /* Compute distance */
   const double r2 = dx * dx + dy * dy + dz * dz;
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 0a5915cd10a170171378d708f3f0cda54c34e9b0..daf0914dacade073e600d7b5befdd9890d5773e2 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -110,10 +110,12 @@ void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
 void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
                            struct cell *restrict cj) {
 
+  /* Some constants */
   const struct engine *e = r->e;
+  const struct space *s = e->s;
+  const int periodic = s->periodic;
+  const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
   const struct gravity_props *props = e->gravity_properties;
-  const int periodic = e->s->periodic;
-  const struct multipole *multi_j = &cj->multipole->m_pole;
   // const float a_smooth = e->gravity_properties->a_smooth;
   // const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
 
@@ -122,6 +124,9 @@ void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
   /* Anything to do here? */
   if (!cell_is_active(ci, e)) return;
 
+  /* Short-cut to the multipole */
+  const struct multipole *multi_j = &cj->multipole->m_pole;
+
 #ifdef SWIFT_DEBUG_CHECKS
   if (ci == cj) error("Interacting a cell with itself using M2L");
 
@@ -136,7 +141,7 @@ void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
 
   /* Let's interact at this level */
   gravity_M2L(&ci->multipole->pot, multi_j, ci->multipole->CoM,
-              cj->multipole->CoM, props, periodic * 0);
+              cj->multipole->CoM, props, periodic, dim);
 
   TIMER_TOC(timer_dopair_grav_mm);
 }