diff --git a/src/engine.c b/src/engine.c
index 5d8dd1ac16de66715a278428d8c4c9ce54b03246..727f9fa8f8e415e66ca4bc019a19d48a5ccacf55 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1343,10 +1343,11 @@ static inline void engine_make_gravity_dependencies(struct scheduler *sched,
                                                     struct cell *c) {
 
   /* init --> gravity --> kick */
-  /* grav_up --> gravity ( --> kick) */
   scheduler_addunlock(sched, c->super->init, gravity);
-  scheduler_addunlock(sched, c->super->grav_up, gravity);
   scheduler_addunlock(sched, gravity, c->super->kick);
+
+  /* grav_up --> gravity ( --> kick) */
+  scheduler_addunlock(sched, c->super->grav_up, gravity);
 }
 
 /**
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 7349ad8e224ccfee8504d65cef2d0f5ac47d9998..e3a181cebee7acb52c5ebb0d1dc671e928fcc646 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -77,7 +77,7 @@ gravity_compute_timestep_self(const struct phys_const* const phys_const,
  */
 __attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
     struct gpart* gp) {
-  gp->epsilon = 0.1;  // MATTHIEU
+  gp->epsilon = 0.;  // MATTHIEU
 }
 
 /**
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 2b3dfc8555c8acb3d21f11772a47b4e3bcc88b45..048d47e737ce70634850c3172eab7dd5ca12414d 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -147,14 +147,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
   /* Apply the gravitational acceleration. */
   const float r = sqrtf(r2);
   const float ir = 1.f / r;
+  const float u = r * rlr_inv;
   const float mrinv3 = multi->mass * ir * ir * ir;
 
+  /* Get long-range correction */
+  float f_lr;
+  kernel_long_grav_eval(u, &f_lr);
+
 #if const_gravity_multipole_order < 2
 
   /* 0th and 1st order terms */
-  gp->a_grav[0] += mrinv3 * dx[0];
-  gp->a_grav[1] += mrinv3 * dx[1];
-  gp->a_grav[2] += mrinv3 * dx[2];
+  gp->a_grav[0] += mrinv3 * f_lr * dx[0];
+  gp->a_grav[1] += mrinv3 * f_lr * dx[1];
+  gp->a_grav[2] += mrinv3 * f_lr * dx[2];
 
   gp->mass_interacted += multi->mass;
 #elif const_gravity_multipole_order == 2
@@ -178,9 +183,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
   const float qRR = qRx * dx[0] + qRy * dx[1] + qRz * dx[2];
   const float C = D1 + 0.5f * D2 * q + 0.5f * D3 * qRR;
 
-  gp->a_grav[0] -= C * dx[0] + D2 * qRx;
-  gp->a_grav[1] -= C * dx[1] + D2 * qRy;
-  gp->a_grav[2] -= C * dx[2] + D2 * qRz;
+  gp->a_grav[0] -= f_lr * (C * dx[0] + D2 * qRx);
+  gp->a_grav[1] -= f_lr * (C * dx[1] + D2 * qRy);
+  gp->a_grav[2] -= f_lr * (C * dx[2] + D2 * qRz);
 
   gp->mass_interacted += multi->mass;
 #else
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index cf41c735e669bb73fa2f4d2ddba3a064bbce175a..f6bced124c2335e306e0fa4c2fad54809d1895e6 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -76,6 +76,8 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
 
   TIMER_TIC;
 
+// message("rlr_inv= %f", rlr_inv);
+
 #ifdef SANITY_CHECKS
   if (gcount == 0) error("Empty cell!");  // MATTHIEU sanity check
 
@@ -487,9 +489,9 @@ static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) {
                           cj->loc[2] - pos_i[2]};  // z
     const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
 
-    if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_pm(r, ci, cj);
-
     if (r2 > max_d2) continue;
+
+    if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_pm(r, ci, cj);
   }
 }
 
diff --git a/src/tools.c b/src/tools.c
index f655e5e14c3d04d2cb21a5f7156aeee6f2c29467..a60c785104a0bb1bb5e508954fb4ee68f5106668 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -496,7 +496,7 @@ void shuffle_particles(struct part *parts, const int count) {
 void gravity_n2(struct gpart *gparts, const int gcount,
                 const struct phys_const *constants, float rlr) {
 
-  const float rlr_inv = 0.;  // 1. / rlr;
+  const float rlr_inv = 1. / rlr;
   const float max_d = const_gravity_r_cut * rlr;
   const float max_d2 = max_d * max_d;