diff --git a/src/gravity.c b/src/gravity.c
index 53ab6b816f964e4e2b071df1e3192d972c5567de..f0968d85f2c215a57e4a7643ab722eb276fb120b 100644
--- a/src/gravity.c
+++ b/src/gravity.c
@@ -652,21 +652,23 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
           SELF_GRAVITY_MULTIPOLE_ORDER);
 
   /* Creare files and write header */
-  const double epsilon = gravity_get_softening(0, e->gravity_properties);
   FILE *file_swift = fopen(file_name_swift, "w");
   fprintf(file_swift, "# Gravity accuracy test - SWIFT FORCES\n");
   fprintf(file_swift, "# G= %16.8e\n", e->physical_constants->const_newton_G);
   fprintf(file_swift, "# N= %d\n", SWIFT_GRAVITY_FORCE_CHECKS);
-  fprintf(file_swift, "# epsilon= %16.8e\n", epsilon);
   fprintf(file_swift, "# periodic= %d\n", s->periodic);
   fprintf(file_swift, "# theta= %16.8e\n", e->gravity_properties->theta_crit);
   fprintf(file_swift, "# Git Branch: %s\n", git_branch());
   fprintf(file_swift, "# Git Revision: %s\n", git_revision());
   fprintf(file_swift,
-          "# %16s %16s %16s %16s %16s %16s %16s %16s %16s %16s %16s %16s\n",
+          "# %16s %16s %16s %16s %16s %16s %16s %16s %16s %16s %16s %16s "
+          "%16s %16s %16s %16s %16s %16s %16s %16s %16s %16s %16s %16s %16s\n",
           "id", "pos[0]", "pos[1]", "pos[2]", "a_swift[0]", "a_swift[1]",
           "a_swift[2]", "potential", "a_PM[0]", "a_PM[1]", "a_PM[2]",
-          "potentialPM");
+          "potentialPM", "a_p2p[0]", "a_p2p[1]", "a_p2p[2]",
+          "a_m2p[0]", "a_m2p[1]", "a_m2p[2]",
+          "a_m2l[0]", "a_m2l[1]", "a_m2l[2]",
+          "n_p2p", "n_m2p", "n_m2l", "n_PM");
 
   /* Output particle SWIFT accelerations  */
   for (size_t i = 0; i < s->nr_gparts; ++i) {
@@ -688,11 +690,17 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
 
       fprintf(file_swift,
               "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e "
-              "%16.8e %16.8e %16.8e\n",
+              "%16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e "
+              "%16.8e %16.8e %16.8e %18lld %18lld %18lld %18lld\n",
               id, gpi->x[0], gpi->x[1], gpi->x[2], gpi->a_grav[0],
               gpi->a_grav[1], gpi->a_grav[2],
               gravity_get_comoving_potential(gpi), gpi->a_grav_PM[0],
-              gpi->a_grav_PM[1], gpi->a_grav_PM[2], gpi->potential_PM);
+              gpi->a_grav_PM[1], gpi->a_grav_PM[2], gpi->potential_PM,
+              gpi->a_grav_p2p[0], gpi->a_grav_p2p[1], gpi->a_grav_p2p[2],
+              gpi->a_grav_m2p[0], gpi->a_grav_m2p[1], gpi->a_grav_m2p[2],
+              gpi->a_grav_m2l[0], gpi->a_grav_m2l[1], gpi->a_grav_m2l[2],
+              gpi->num_interacted_p2p, gpi->num_interacted_m2p,
+              gpi->num_interacted_m2l, gpi->num_not_interacted);
     }
   }
 
@@ -714,7 +722,6 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
     fprintf(file_exact, "# Gravity accuracy test - EXACT FORCES\n");
     fprintf(file_exact, "# G= %16.8e\n", e->physical_constants->const_newton_G);
     fprintf(file_exact, "# N= %d\n", SWIFT_GRAVITY_FORCE_CHECKS);
-    fprintf(file_exact, "# epsilon=%16.8e\n", epsilon);
     fprintf(file_exact, "# periodic= %d\n", s->periodic);
     fprintf(file_exact, "# Git Branch: %s\n", git_branch());
     fprintf(file_exact, "# Git Revision: %s\n", git_revision());
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index a8ace72d2b2bb8cb92d2d386ae3ed5402c4f76e0..1b176cf593574c63817059ce943a9373556d1ca8 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -152,13 +152,27 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   gp->potential_PM = 0.f;
-  gp->a_grav_PM[0] = 0.f;
-  gp->a_grav_PM[1] = 0.f;
-  gp->a_grav_PM[2] = 0.f;
+
+  /* Track accelerations of each component. */
+  for (int i=0; i<3; i++) {
+    gp->a_grav_PM[i] = 0.f;
+    gp->a_grav_p2p[i] = 0.f;
+    gp->a_grav_m2p[i] = 0.f;
+    gp->a_grav_m2l[i] = 0.f;
+  }
+
+  /* Interaction counters. */  
+  gp->num_interacted_m2p = 0;
+  gp->num_interacted_m2l = 0;
+  gp->num_interacted_p2p = 0;
+  gp->num_not_interacted = 0;
 #endif
 
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
   gp->num_interacted = 0;
+#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
   gp->initialised = 1;
 #endif
 }
@@ -188,9 +202,12 @@ __attribute__((always_inline)) INLINE static void gravity_end_force(
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   gp->potential_PM *= const_G;
-  gp->a_grav_PM[0] *= const_G;
-  gp->a_grav_PM[1] *= const_G;
-  gp->a_grav_PM[2] *= const_G;
+  for (int i = 0; i < 3; i++) {
+    gp->a_grav_PM[i] *= const_G;
+    gp->a_grav_p2p[i] *= const_G;
+    gp->a_grav_m2p[i] *= const_G;
+    gp->a_grav_m2l[i] *= const_G;
+  }
 #endif
 
 #ifdef SWIFT_DEBUG_CHECKS
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index bf0d612b077579829bffd69b5e821d86c75345c1..79459d0330f50f757738a4a4b0ce6990bcae6972 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -58,9 +58,6 @@ struct gpart {
 
 #ifdef SWIFT_DEBUG_CHECKS
 
-  /* Numer of gparts this gpart interacted with */
-  long long num_interacted;
-
   /* Time of the last drift */
   integertime_t ti_drift;
 
@@ -72,6 +69,11 @@ struct gpart {
 
 #endif
 
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
+  /* Total number of gparts this gpart interacted with */
+  long long num_interacted;
+#endif
+
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
 
   /*! Acceleration taken from the mesh only */
@@ -80,11 +82,24 @@ struct gpart {
   /*! Potential taken from the mesh only */
   float potential_PM;
 
-  /* Brute-force particle acceleration. */
+  /* Acceleration taken from each component of the tree */
+  float a_grav_p2p[3];
+  float a_grav_m2p[3];
+  float a_grav_m2l[3];
+
+  /* Brute-force particle accelerations */
   double a_grav_exact[3];
+  double a_grav_exact_long[3];
+  double a_grav_exact_short[3];
 
   /* Brute-force particle potential. */
   double potential_exact;
+
+  /* Type specific interaction counters */
+  long long num_interacted_m2p;
+  long long num_interacted_m2l;
+  long long num_interacted_p2p;
+  long long num_not_interacted;
 #endif
 
 } SWIFT_STRUCT_ALIGN;
diff --git a/src/gravity/MultiSoftening/gravity.h b/src/gravity/MultiSoftening/gravity.h
index 893a835011c244a99b443edbd0f095aa5961cb2e..62f245bcc30d37a7fb64c8784a388006fc97b7e8 100644
--- a/src/gravity/MultiSoftening/gravity.h
+++ b/src/gravity/MultiSoftening/gravity.h
@@ -162,13 +162,27 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   gp->potential_PM = 0.f;
-  gp->a_grav_PM[0] = 0.f;
-  gp->a_grav_PM[1] = 0.f;
-  gp->a_grav_PM[2] = 0.f;
+
+  /* Track accelerations of each component. */
+  for (int i=0; i<3; i++) {
+    gp->a_grav_PM[i] = 0.f;
+    gp->a_grav_p2p[i] = 0.f;
+    gp->a_grav_m2p[i] = 0.f;
+    gp->a_grav_m2l[i] = 0.f;
+  }
+
+  /* Interaction counters. */
+  gp->num_interacted_m2p = 0;
+  gp->num_interacted_m2l = 0;
+  gp->num_interacted_p2p = 0;
+  gp->num_not_interacted = 0;
 #endif
 
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
   gp->num_interacted = 0;
+#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
   gp->initialised = 1;
 #endif
 }
@@ -202,9 +216,12 @@ __attribute__((always_inline)) INLINE static void gravity_end_force(
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   gp->potential_PM *= const_G;
-  gp->a_grav_PM[0] *= const_G;
-  gp->a_grav_PM[1] *= const_G;
-  gp->a_grav_PM[2] *= const_G;
+  for (int i = 0; i < 3; i++) {
+    gp->a_grav_PM[i] *= const_G;
+    gp->a_grav_p2p[i] *= const_G;
+    gp->a_grav_m2p[i] *= const_G;
+    gp->a_grav_m2l[i] *= const_G;
+  }
 #endif
 
 #ifdef SWIFT_DEBUG_CHECKS
diff --git a/src/gravity/MultiSoftening/gravity_part.h b/src/gravity/MultiSoftening/gravity_part.h
index e865d99484735820ebb81e7a10fc3b5bb3edfdb2..fcc212bbee9d30b1c1fe5c32f9a5311cc0f7ec92 100644
--- a/src/gravity/MultiSoftening/gravity_part.h
+++ b/src/gravity/MultiSoftening/gravity_part.h
@@ -59,9 +59,6 @@ struct gpart {
 
 #ifdef SWIFT_DEBUG_CHECKS
 
-  /* Numer of gparts this gpart interacted with */
-  long long num_interacted;
-
   /* Time of the last drift */
   integertime_t ti_drift;
 
@@ -73,6 +70,11 @@ struct gpart {
 
 #endif
 
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
+  /* Total number of gparts this gpart interacted with */
+  long long num_interacted;
+#endif
+
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
 
   /*! Acceleration taken from the mesh only */
@@ -81,11 +83,24 @@ struct gpart {
   /*! Potential taken from the mesh only */
   float potential_PM;
 
-  /* Brute-force particle acceleration. */
+  /* Acceleration taken from each component of the tree */
+  float a_grav_p2p[3];
+  float a_grav_m2p[3];
+  float a_grav_m2l[3];
+
+  /* Brute-force particle accelerations */
   double a_grav_exact[3];
+  double a_grav_exact_long[3];
+  double a_grav_exact_short[3];
 
   /* Brute-force particle potential. */
   double potential_exact;
+
+  /* Type specific interaction counters */
+  long long num_interacted_m2l;
+  long long num_interacted_m2p;
+  long long num_interacted_p2p;
+  long long num_not_interacted;
 #endif
 
 } SWIFT_STRUCT_ALIGN;
diff --git a/src/multipole.h b/src/multipole.h
index 4cbac0e7634982b7c25574126858ac049094b5a8..83277ef82f541d7ac69ba8db493f5051feaa1928 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -92,11 +92,18 @@ struct grav_tensor {
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 5
 #error "Missing implementation for order >5"
 #endif
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
 
   /* Total number of gpart this field tensor interacted with */
   long long num_interacted;
 
+  /* Total number of gpart this field tensor did not interact
+   * with. i.e., the distance to the cell was > r_cut_max. */
+  long long num_not_interacted;
+#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
+
   /* Last time this tensor was zeroed */
   integertime_t ti_init;
 
@@ -168,7 +175,7 @@ struct multipole {
 #error "Missing implementation for order >5"
 #endif
 
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
 
   /* Total number of gpart in this multipole */
   long long num_gpart;
@@ -307,8 +314,12 @@ INLINE static void gravity_field_tensors_init(struct grav_tensor *l,
 INLINE static void gravity_field_tensors_add(
     struct grav_tensor *restrict la, const struct grav_tensor *restrict lb) {
 #ifdef SWIFT_DEBUG_CHECKS
-  if (lb->num_interacted == 0) error("Adding tensors that did not interact");
+  if (lb->num_interacted + lb->num_not_interacted == 0)
+      error("Adding tensors that did not interact");
+#endif
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
   la->num_interacted += lb->num_interacted;
+  la->num_not_interacted += lb->num_not_interacted;
 #endif
 
   la->interacted = 1;
@@ -597,7 +608,7 @@ INLINE static void gravity_multipole_add(struct multipole *restrict ma,
 #error "Missing implementation for order >5"
 #endif
 
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
   ma->num_gpart += mb->num_gpart;
 #endif
 }
@@ -1319,7 +1330,7 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
 #error "Missing implementation for order >5"
 #endif
 
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
   multi->m_pole.num_gpart = gcount;
 #endif
 }
@@ -1576,7 +1587,7 @@ INLINE static void gravity_M2M(struct multipole *restrict m_a,
 #error "Missing implementation for order >5"
 #endif
 
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
   m_a->num_gpart = m_b->num_gpart;
 #endif
 }
@@ -1594,7 +1605,7 @@ INLINE static void gravity_M2L_apply(
     struct grav_tensor *restrict l_b, const struct multipole *restrict m_a,
     const struct potential_derivatives_M2L *pot) {
 
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
   /* Count interactions */
   l_b->num_interacted += m_a->num_gpart;
 #endif
@@ -2090,8 +2101,12 @@ INLINE static void gravity_L2L(struct grav_tensor *restrict la,
   gravity_field_tensors_init(la, 0);
 
 #ifdef SWIFT_DEBUG_CHECKS
-  if (lb->num_interacted == 0) error("Shifting tensors that did not interact");
+  if (lb->num_interacted + lb->num_not_interacted == 0)
+      error("Shifting tensors that did not interact");
+#endif
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
   la->num_interacted = lb->num_interacted;
+  la->num_not_interacted = lb->num_not_interacted;
 #endif
 
   /* Distance to shift by */
@@ -2443,8 +2458,15 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
                                const double loc[3], struct gpart *gp) {
 
 #ifdef SWIFT_DEBUG_CHECKS
-  if (lb->num_interacted == 0) error("Interacting with empty field tensor");
-  gp->num_interacted += lb->num_interacted;
+  if (lb->num_interacted + lb->num_not_interacted == 0)
+      error("Interacting with empty field tensor");
+#endif
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
+  gp->num_interacted += lb->num_interacted + lb->num_not_interacted;
+#endif
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  gp->num_interacted_m2l += lb->num_interacted;
+  gp->num_not_interacted += lb->num_not_interacted;
 #endif
 
   /* Local accumulator */
@@ -2565,6 +2587,12 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
   gp->a_grav[1] += a_grav[1];
   gp->a_grav[2] += a_grav[2];
   gravity_add_comoving_potential(gp, pot);
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  gp->a_grav_m2l[0] += a_grav[0];
+  gp->a_grav_m2l[1] += a_grav[1];
+  gp->a_grav_m2l[2] += a_grav[2];
+#endif
 }
 
 /**
diff --git a/src/runner_doiact_grav.c b/src/runner_doiact_grav.c
index a8a07162b69613354bca35b4bc72fd27fc19e0d0..4d7ba1c34cc7bdc74db1ca08e94cd1adbf49b872 100644
--- a/src/runner_doiact_grav.c
+++ b/src/runner_doiact_grav.c
@@ -266,6 +266,12 @@ static INLINE void runner_dopair_grav_pp_full(
       if (pjd < gcount_j && !gpart_is_inhibited(&gparts_j[pjd], e))
         gparts_i[pid].num_interacted++;
 #endif
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+      /* Update the p2p interaction counter if it's not a padded gpart */
+      if (pjd < gcount_j && !gpart_is_inhibited(&gparts_j[pjd], e))
+        gparts_i[pid].num_interacted_p2p++;
+#endif
     }
 
     /* Store everything back in cache */
@@ -273,6 +279,12 @@ static INLINE void runner_dopair_grav_pp_full(
     ci_cache->a_y[pid] += a_y;
     ci_cache->a_z[pid] += a_z;
     ci_cache->pot[pid] += pot;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+    gparts_i[pid].a_grav_p2p[0] += a_x;
+    gparts_i[pid].a_grav_p2p[1] += a_y;
+    gparts_i[pid].a_grav_p2p[2] += a_z;
+#endif
   }
 }
 
@@ -410,6 +422,12 @@ static INLINE void runner_dopair_grav_pp_truncated(
       if (pjd < gcount_j && !gpart_is_inhibited(&gparts_j[pjd], e))
         gparts_i[pid].num_interacted++;
 #endif
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+      /* Update the p2p interaction counter if it's not a padded gpart */
+      if (pjd < gcount_j && !gpart_is_inhibited(&gparts_j[pjd], e))
+        gparts_i[pid].num_interacted_p2p++;
+#endif
     }
 
     /* Store everything back in cache */
@@ -417,6 +435,12 @@ static INLINE void runner_dopair_grav_pp_truncated(
     ci_cache->a_y[pid] += a_y;
     ci_cache->a_z[pid] += a_z;
     ci_cache->pot[pid] += pot;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+    gparts_i[pid].a_grav_p2p[0] += a_x;
+    gparts_i[pid].a_grav_p2p[1] += a_y;
+    gparts_i[pid].a_grav_p2p[2] += a_z;
+#endif
   }
 }
 
@@ -543,6 +567,16 @@ static INLINE void runner_dopair_grav_pm_full(
     if (pid < gcount_i)
       gparts_i[pid].num_interacted += cj->grav.multipole->m_pole.num_gpart;
 #endif
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+    /* Update the M2P interaction counter and forces. */
+    if (pid < gcount_i) {
+      gparts_i[pid].num_interacted_m2p += cj->grav.multipole->m_pole.num_gpart;
+      gparts_i[pid].a_grav_m2p[0] += f_x;
+      gparts_i[pid].a_grav_m2p[1] += f_y;
+      gparts_i[pid].a_grav_m2p[2] += f_z;
+    }
+#endif
   }
 }
 
@@ -674,6 +708,16 @@ static INLINE void runner_dopair_grav_pm_truncated(
     if (pid < gcount_i)
       gparts_i[pid].num_interacted += cj->grav.multipole->m_pole.num_gpart;
 #endif
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+    /* Update the M2P interaction counter and forces. */
+    if (pid < gcount_i) {
+      gparts_i[pid].num_interacted_m2p += cj->grav.multipole->m_pole.num_gpart;
+      gparts_i[pid].a_grav_m2p[0] += f_x;
+      gparts_i[pid].a_grav_m2p[1] += f_y;
+      gparts_i[pid].a_grav_m2p[2] += f_z;
+    }
+#endif
   }
 }
 
@@ -999,6 +1043,12 @@ static INLINE void runner_doself_grav_pp_full(
       if (pjd < gcount && !gpart_is_inhibited(&gparts[pjd], e))
         gparts[pid].num_interacted++;
 #endif
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+      /* Update the P2P interaction counter if it's not a padded gpart */
+      if (pjd < gcount && !gpart_is_inhibited(&gparts[pjd], e))
+        gparts[pid].num_interacted_p2p++;
+#endif
     }
 
     /* Store everything back in cache */
@@ -1006,6 +1056,12 @@ static INLINE void runner_doself_grav_pp_full(
     ci_cache->a_y[pid] += a_y;
     ci_cache->a_z[pid] += a_z;
     ci_cache->pot[pid] += pot;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+    gparts[pid].a_grav_p2p[0] += a_x;
+    gparts[pid].a_grav_p2p[1] += a_y;
+    gparts[pid].a_grav_p2p[2] += a_z;
+#endif
   }
 }
 
@@ -1126,6 +1182,12 @@ static INLINE void runner_doself_grav_pp_truncated(
       if (pjd < gcount && !gpart_is_inhibited(&gparts[pjd], e))
         gparts[pid].num_interacted++;
 #endif
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+    /* Update the P2P interaction counter if it's not a padded gpart */
+    if (pjd < gcount && !gpart_is_inhibited(&gparts[pjd], e))
+      gparts[pid].num_interacted_p2p++;
+#endif
     }
 
     /* Store everything back in cache */
@@ -1133,6 +1195,12 @@ static INLINE void runner_doself_grav_pp_truncated(
     ci_cache->a_y[pid] += a_y;
     ci_cache->a_z[pid] += a_z;
     ci_cache->pot[pid] += pot;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+    gparts[pid].a_grav_p2p[0] += a_x;
+    gparts[pid].a_grav_p2p[1] += a_y;
+    gparts[pid].a_grav_p2p[2] += a_z;
+#endif
   }
 }
 
@@ -1568,12 +1636,12 @@ void runner_dopair_recursive_grav(struct runner *r, struct cell *ci,
   /* Are we beyond the distance where the truncated forces are 0? */
   if (periodic && r_lr_check > max_distance) {
 
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
     /* Need to account for the interactions we missed */
     if (cell_is_active_gravity(ci, e))
-      multi_i->pot.num_interacted += multi_j->m_pole.num_gpart;
+      multi_i->pot.num_not_interacted += multi_j->m_pole.num_gpart;
     if (cell_is_active_gravity(cj, e))
-      multi_j->pot.num_interacted += multi_i->m_pole.num_gpart;
+      multi_j->pot.num_not_interacted += multi_i->m_pole.num_gpart;
 #endif
     return;
   }
@@ -1777,9 +1845,9 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
       /* Are we beyond the distance where the truncated forces are 0 ?*/
       if (min_radius2 > max_distance2) {
 
-#ifdef SWIFT_DEBUG_CHECKS
+#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_GRAVITY_FORCE_CHECKS)
         /* Need to account for the interactions we missed */
-        multi_i->pot.num_interacted += multi_j->m_pole.num_gpart;
+        multi_i->pot.num_not_interacted += multi_j->m_pole.num_gpart;
 #endif
 
         /* Record that this multipole received a contribution */