diff --git a/src/hydro/DebugInteractions/hydro.h b/src/hydro/DebugInteractions/hydro.h
index dfffab642ae99e3bbb1928984e8daaf44d76c8eb..194c9e52bcc549f6633e34e571aee2a9012e7d0f 100644
--- a/src/hydro/DebugInteractions/hydro.h
+++ b/src/hydro/DebugInteractions/hydro.h
@@ -24,6 +24,7 @@
  * @brief Empty SPH implementation used solely to test the SELF/PAIR routines.
  */
 
+#include "equation_of_state.h"
 #include "hydro_properties.h"
 #include "hydro_space.h"
 #include "part.h"
diff --git a/src/hydro/Shadowswift/hydro_iact.h b/src/hydro/Shadowswift/hydro_iact.h
index 28c3ce5521c8f15a971d47dbb3e48cdb2a57e77d..15219bf879ea7d34a8a45be217ef81107d3392e6 100644
--- a/src/hydro/Shadowswift/hydro_iact.h
+++ b/src/hydro/Shadowswift/hydro_iact.h
@@ -345,4 +345,3 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0);
 }
-
diff --git a/src/runner_doiact.h b/src/runner_doiact.h
index 5812b9641dbb371eff74c1b9a1c4c68b8f1d3f79..f73f44d7db005d6e7689c8d9e134d4e448b889cb 100644
--- a/src/runner_doiact.h
+++ b/src/runner_doiact.h
@@ -940,9 +940,17 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
   const int count_j = cj->count;
   struct part *restrict parts_i = ci->parts;
   struct part *restrict parts_j = cj->parts;
+
+  /* Maximal displacement since last rebuild */
+  const double dx_max = (ci->dx_max_sort + cj->dx_max_sort);
+
+  /* Position on the axis of the particles closest to the interface */
   const double di_max = sort_i[count_i - 1].d;
   const double dj_min = sort_j[0].d;
-  const double dx_max = (ci->dx_max_sort + cj->dx_max_sort);
+
+  /* Upper bound on maximal distance in the other cell */
+  const double max_i = max(hi_max, hj_max) * kernel_gamma + dx_max - rshift;
+  const double max_j = max(hi_max, hj_max) * kernel_gamma + dx_max;
 
   if (cell_is_active(ci, e)) {
 
@@ -950,6 +958,9 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
      * We start from the centre and move outwards. */
     for (int pid = count_i - 1; pid >= 0; pid--) {
 
+      /* Did we go beyond the upper bound ? */
+      if (sort_i[pid].d + max_i < dj_min) break;
+
       /* Get a hold of the ith part in ci. */
       struct part *pi = &parts_i[sort_i[pid].i];
       const float hi = pi->h;
@@ -963,11 +974,8 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
       if (di < dj_min) continue;
 
       /* Where do we have to stop when looping over cell j? */
-      int last_pj = count_j - 1;
-      while (sort_j[last_pj].d > di && last_pj > 0) last_pj--;
-
-      last_pj += 1;
-      if (last_pj >= count_j) last_pj = count_j - 1;
+      int last_pj = 0;
+      while (last_pj < count_j - 1 && sort_j[last_pj].d < di) last_pj++;
 
       /* Get some additional information about pi */
       const float hig2 = hi * hi * kernel_gamma2;
@@ -1004,6 +1012,9 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
      * We start from the centre and move outwards. */
     for (int pjd = 0; pjd < count_j; pjd++) {
 
+      /* Did we go beyond the upper bound ? */
+      if (sort_j[pjd].d - max_j > di_max - rshift) break;
+
       /* Get a hold of the jth part in cj. */
       struct part *pj = &parts_j[sort_j[pjd].i];
       const float hj = pj->h;
@@ -1016,12 +1027,8 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
       if (dj > di_max - rshift) continue;
 
       /* Where do we have to stop when looping over cell i? */
-      int first_pi = 0;
-      while (sort_i[first_pi].d - rshift < dj && first_pi < count_i - 1)
-        first_pi++;
-
-      first_pi -= 1;
-      if (first_pi < 0) first_pi = 0;
+      int first_pi = count_i;
+      while (first_pi > 0 && sort_i[first_pi].d - rshift > dj) first_pi--;
 
       /* Get some additional information about pj */
       const float hjg2 = hj * hj * kernel_gamma2;