diff --git a/configure.ac b/configure.ac
index 539dea934771fcc125a0a13e280379e399274ea6..238e29468dd1a06d054b44d6d78b73afb3555407 100644
--- a/configure.ac
+++ b/configure.ac
@@ -668,9 +668,6 @@ case "$with_hydro" in
    minimal)
       AC_DEFINE([MINIMAL_SPH], [1], [Minimal SPH])
    ;; 
-   debug)
-      AC_DEFINE([DEBUG_INTERACTIONS_SPH], [1], [Debug SPH])
-   ;; 
    hopkins)
       AC_DEFINE([HOPKINS_PE_SPH], [1], [Pressure-Entropy SPH])
    ;; 
@@ -689,6 +686,29 @@ case "$with_hydro" in
    ;;
 esac
 
+# Check if debugging interactions is switched on.
+AC_ARG_ENABLE([debug-interactions],
+   [AS_HELP_STRING([--enable-debug-interactions],
+     [Activate interaction debugging, logging a maximum of @<:@N@:>@ neighbours. Defaults to 256 if no value set.]
+   )],
+   [enable_debug_interactions="$enableval"],
+   [enable_debug_interactions="no"]
+)
+if test "$enable_debug_interactions" != "no"; then
+  if test "$with_hydro" = "gadget2"; then
+      AC_DEFINE([DEBUG_INTERACTIONS_SPH],1,[Enable interaction debugging])
+    if test "$enable_debug_interactions" == "yes"; then
+      AC_DEFINE([MAX_NUM_OF_NEIGHBOURS],256,[The maximum number of particle neighbours to be logged])
+      [enable_debug_interactions="yes (Logging up to 256 neighbours)"]
+    else
+      AC_DEFINE_UNQUOTED([MAX_NUM_OF_NEIGHBOURS], [$enableval] ,[The maximum number of particle neighbours to be logged])
+      [enable_debug_interactions="yes (Logging up to $enableval neighbours)"]
+    fi
+  else 
+    [enable_debug_interactions="no (only available for gadget2 hydro scheme)"]
+  fi
+fi
+
 # SPH Kernel function
 AC_ARG_WITH([kernel],
    [AS_HELP_STRING([--with-kernel=<kernel>],
@@ -945,11 +965,12 @@ AC_MSG_RESULT([
    Multipole order     : $with_multipole_order
    No gravity below ID : $no_gravity_below_id
 
-   Individual timers    : $enable_timers
-   Task debugging       : $enable_task_debugging
-   Threadpool debugging : $enable_threadpool_debugging
-   Debugging checks     : $enable_debugging_checks
-   Naive interactions   : $enable_naive_interactions
-   Gravity checks       : $gravity_force_checks
+   Individual timers     : $enable_timers
+   Task debugging        : $enable_task_debugging
+   Threadpool debugging  : $enable_threadpool_debugging
+   Debugging checks      : $enable_debugging_checks
+   Interaction debugging : $enable_debug_interactions
+   Naive interactions    : $enable_naive_interactions
+   Gravity checks        : $gravity_force_checks
 
  ------------------------])
diff --git a/examples/check_interactions.sh b/examples/check_interactions.sh
new file mode 100755
index 0000000000000000000000000000000000000000..c61daa6f231e0a066c525acd4a522dc7ea4b2bcd
--- /dev/null
+++ b/examples/check_interactions.sh
@@ -0,0 +1,128 @@
+#!/bin/bash -l
+#
+# Runs the SedovBlast_3D and EAGLE_12 examples using naive, serial and vectorised particle interactions. Then compares the output between versions using check_ngbs.py
+
+echo
+echo "# Running SedovBlast_3D and EAGLE_12 with naive interactions and neighbour logging, 16 thread"
+echo
+
+cd ../
+
+./autogen.sh
+
+# Naive interactions run
+./configure --disable-mpi --enable-debug-interactions=1024 --disable-vec --enable-naive-interactions
+make clean; make -j 6
+
+cd examples/SedovBlast_3D/
+
+./getGlass.sh
+python makeIC.py
+
+../swift -s -t 16 -n 5 sedov.yml -P SPH:h_tolerance:10
+
+mv sedov_0000.hdf5 sedov_naive.hdf5
+
+cd ../EAGLE_12/
+
+# Link to ICs
+ln -s /gpfs/data/Swift/web-storage/ICs/EAGLE_ICs_12.hdf5 EAGLE_ICs_12.hdf5
+
+../swift -s -t 16 -n 5 eagle_12.yml -P SPH:h_tolerance:10  
+
+mv eagle_0000.hdf5 eagle_12_naive.hdf5
+
+cd ../../
+
+echo
+echo "# Running SedovBlast_3D and EAGLE_12 with serial interactions and neighbour logging, 16 thread"
+echo
+
+# Serial interactions run
+./configure --disable-mpi --enable-debug-interactions=1024 --disable-vec
+make clean; make -j 6
+
+cd examples/SedovBlast_3D/
+
+../swift -s -t 16 -n 5 sedov.yml -P SPH:h_tolerance:10
+
+mv sedov_0000.hdf5 sedov_serial.hdf5
+
+cd ../EAGLE_12/
+
+../swift -s -t 16 -n 5 eagle_12.yml -P SPH:h_tolerance:10  
+
+mv eagle_0000.hdf5 eagle_12_serial.hdf5
+
+cd ../../
+
+echo
+echo "# Running SedovBlast_3D and EAGLE_12 with vectorised interactions and neighbour logging, 16 thread"
+echo
+
+# Vectorised interactions run
+./configure --disable-mpi --enable-debug-interactions=1024
+make clean; make -j 6
+
+cd examples/SedovBlast_3D/
+
+../swift -s -t 16 -n 5 sedov.yml -P SPH:h_tolerance:10
+
+mv sedov_0000.hdf5 sedov_vec.hdf5
+
+# Compare outputs
+if python ../check_ngbs.py sedov_naive.hdf5 sedov_serial.hdf5 
+then
+  echo "SedovBlast_3D comparison between naive and serial passed"
+else
+  echo "SedovBlast_3D comparison between naive and serial failed"
+  exit 1
+fi
+
+if python ../check_ngbs.py sedov_naive.hdf5 sedov_vec.hdf5 
+then
+  echo "SedovBlast_3D comparison between naive and vectorised passed"
+else
+  echo "SedovBlast_3D comparison between naive and vectorised failed"
+  exit 1
+fi
+
+if python ../check_ngbs.py sedov_serial.hdf5 sedov_vec.hdf5 
+then
+  echo "SedovBlast_3D comparison between serial and vectorised passed"
+else
+  echo "SedovBlast_3D comparison between serial and vectorised failed"
+  exit 1
+fi
+
+cd ../EAGLE_12/
+
+../swift -s -t 16 -n 5 eagle_12.yml -P SPH:h_tolerance:10  
+
+mv eagle_0000.hdf5 eagle_12_vec.hdf5
+
+# Compare outputs
+if python ../check_ngbs.py eagle_12_naive.hdf5 eagle_12_serial.hdf5 
+then
+  echo "EAGLE_12 comparison between naive and serial passed"
+else
+  echo "EAGLE_12 comparison between naive and serial failed"
+  exit 1
+fi
+
+if python ../check_ngbs.py eagle_12_naive.hdf5 eagle_12_vec.hdf5 
+then
+  echo "EAGLE_12 comparison between naive and vectorised passed"
+else
+  echo "EAGLE_12 comparison between naive and vectorised failed"
+  exit 1
+fi
+
+if python ../check_ngbs.py eagle_12_serial.hdf5 eagle_12_vec.hdf5 
+then
+  echo "EAGLE_12 comparison between serial and vectorised passed"
+  exit 0
+else
+  echo "EAGLE_12 comparison between serial and vectorised failed"
+  exit 1
+fi
diff --git a/examples/check_ngbs.py b/examples/check_ngbs.py
new file mode 100644
index 0000000000000000000000000000000000000000..6b403c6a4555c442e3004921748c57890cdb2e0f
--- /dev/null
+++ b/examples/check_ngbs.py
@@ -0,0 +1,270 @@
+import h5py as h
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import os.path
+
+kernel_gamma = 1.825742
+kernel_gamma2 = kernel_gamma * kernel_gamma
+kernel_gamma_dim = np.power(kernel_gamma,3)
+hydro_dimension_unit_sphere = 4. * np.pi / 3.
+kernel_norm = hydro_dimension_unit_sphere * kernel_gamma_dim
+error = False
+
+inputFile1 = ""
+inputFile2 = ""
+
+# Check list of density neighbours and check that they are correct.
+def check_density_neighbours(pids, ngb_ids_naive, ngb_ids_sort, mask, pos, h, num_invalid, acc):
+
+    error_val = False
+
+    for k in range(0,num_invalid):
+
+        # Filter neighbour lists for valid particle ids
+        filter_neigh_naive = [i for i in ngb_ids_naive[mask][k] if i > -1]
+        filter_neigh_sort = [i for i in ngb_ids_sort[mask][k] if i > -1]
+
+        # Check neighbour lists for differences
+        id_list = set(filter_neigh_naive).symmetric_difference(set(filter_neigh_sort))
+        
+        pid = pids[mask][k]
+
+        # Loop over discrepancies and check if they are actually neighbours
+        for pjd in id_list:
+            pi_pos = pos[np.where(pids == pid)]
+            pj_pos = pos[np.where(pids == pjd)]
+            
+            hi = h[np.where(pids == pid)]
+            
+            dx = pi_pos[0][0] - pj_pos[0][0]
+            dy = pi_pos[0][1] - pj_pos[0][1]
+            dz = pi_pos[0][2] - pj_pos[0][2]
+            
+            r2 = dx*dx + dy*dy + dz*dz
+            
+            hig2 = hi*hi*kernel_gamma2
+            
+            diff = abs(r2 - hig2)
+            
+            print "Particle {} is missing {}, hig2: {}, r2: {}, |r2 - hig2|: {}".format(pid,pjd,hig2, r2, diff)
+            
+            if diff < acc * hig2:
+                print "Missing interaction due to precision issue will be ignored."
+            else:
+                error_val = True
+
+    return error_val
+
+# Check list of force neighbours and check that they are correct.
+def check_force_neighbours(pids, ngb_ids_naive, ngb_ids_sort, mask, pos, h, num_invalid, acc):
+
+    error_val = False
+
+    for k in range(0,num_invalid):
+
+        # Filter neighbour lists for valid particle ids
+        filter_neigh_naive = [i for i in ngb_ids_naive[mask][k] if i > -1]
+        filter_neigh_sort = [i for i in ngb_ids_sort[mask][k] if i > -1]
+
+        # Check neighbour lists for differences
+        id_list = set(filter_neigh_naive).symmetric_difference(set(filter_neigh_sort))
+        
+        pid = pids[mask][k]
+
+        # Loop over discrepancies and check if they are actually neighbours
+        for pjd in id_list:
+            pi_pos = pos[np.where(pids == pid)]
+            pj_pos = pos[np.where(pids == pjd)]
+            
+            hi = h[np.where(pids == pid)]
+            hj = h[np.where(pids == pjd)]
+            
+            dx = pi_pos[0][0] - pj_pos[0][0]
+            dy = pi_pos[0][1] - pj_pos[0][1]
+            dz = pi_pos[0][2] - pj_pos[0][2]
+            
+            r2 = dx*dx + dy*dy + dz*dz
+            
+            hig2 = hi*hi*kernel_gamma2
+            hjg2 = hj*hj*kernel_gamma2
+            
+            diff = abs(r2 - max(hig2, hjg2))
+            
+            print "Particle {} is missing {}, hig2: {}, hjg2: {}, r2: {}, |r2 - max(hig2,hjg2)|: {}".format(pid,pjd,hig2, hjg2, r2, diff)
+
+            
+            if diff < acc * max(hig2,hjg2):
+                print "Missing interaction due to precision issue will be ignored."
+            else:
+                error_val = True
+
+    return error_val
+
+# Parse command line arguments
+if len(sys.argv) < 3:
+    print "Error: pass input files as arguments"
+    sys.exit()
+else:
+    inputFile1 = sys.argv[1]
+    inputFile2 = sys.argv[2]
+    if os.path.exists(inputFile1) != 1:
+        print "\n{} does not exist!\n".format(inputFile1)
+        sys.exit()
+    if os.path.exists(inputFile2) != 1:
+        print "\n{} does not exist!\n".format(inputFile2)
+        sys.exit()
+
+# Open input files    
+file_naive = h.File(inputFile1, "r")
+file_sort = h.File(inputFile2, "r")
+
+# Read input file fields
+ids_naive = file_naive["/PartType0/ParticleIDs"][:]
+ids_sort = file_sort["/PartType0/ParticleIDs"][:]
+
+h_naive = file_naive["/PartType0/SmoothingLength"][:]
+h_sort = file_sort["/PartType0/SmoothingLength"][:]
+
+pos_naive = file_naive["/PartType0/Coordinates"][:,:]
+pos_sort = file_sort["/PartType0/Coordinates"][:,:]
+
+num_density_naive = file_naive["/PartType0/Num_ngb_density"][:]
+num_density_sort = file_sort["/PartType0/Num_ngb_density"][:]
+
+num_force_naive = file_naive["/PartType0/Num_ngb_force"][:]
+num_force_sort = file_sort["/PartType0/Num_ngb_force"][:]
+
+neighbour_ids_density_naive = file_naive["/PartType0/Ids_ngb_density"][:]
+neighbour_ids_density_sort = file_sort["/PartType0/Ids_ngb_density"][:]
+
+neighbour_ids_force_naive = file_naive["/PartType0/Ids_ngb_force"][:]
+neighbour_ids_force_sort = file_sort["/PartType0/Ids_ngb_force"][:]
+
+#wcount_naive = file_naive["/PartType0/Wcount"][:]
+#wcount_sort = file_sort["/PartType0/Wcount"][:]
+#
+#wcount_naive = wcount_naive * np.power(h_naive,3) * kernel_norm
+#wcount_sort = wcount_sort * np.power(h_sort,3) * kernel_norm
+
+# Cross check
+max_density_ngbs_naive = np.max(num_density_naive)
+max_density_ngbs_sort = np.max(num_density_sort)
+max_force_ngbs_naive = np.max(num_force_naive)
+max_force_ngbs_sort = np.max(num_force_sort)
+
+print "                   Min     Mean     Max "
+print "                   ---------------------"
+print "Ngbs density naiv: ", np.min(num_density_naive), np.mean(num_density_naive), max_density_ngbs_naive
+print "Ngbs density sort: ", np.min(num_density_sort), np.mean(num_density_sort), max_density_ngbs_sort
+print "Ngbs force naiv:   ", np.min(num_force_naive), np.mean(num_force_naive), max_force_ngbs_naive
+print "Ngbs force sort:   ", np.min(num_force_sort), np.mean(num_force_sort), max_force_ngbs_sort
+#print "Wcount naiv:   ", np.min(wcount_naive), np.mean(wcount_naive), np.max(wcount_naive)
+#print "Wcount sort:   ", np.min(wcount_sort), np.mean(wcount_sort), np.max(wcount_sort)
+
+# Sort
+index_naive = np.argsort(ids_naive)
+index_sort = np.argsort(ids_sort)
+
+num_density_naive = num_density_naive[index_naive]
+num_density_sort = num_density_sort[index_sort]
+num_force_naive = num_force_naive[index_naive]
+num_force_sort = num_force_sort[index_sort]
+ids_naive = ids_naive[index_naive]
+ids_sort = ids_sort[index_sort]
+neighbour_ids_density_naive = neighbour_ids_density_naive[index_naive]
+neighbour_ids_density_sort = neighbour_ids_density_sort[index_sort]
+neighbour_ids_force_naive = neighbour_ids_force_naive[index_naive]
+neighbour_ids_force_sort = neighbour_ids_force_sort[index_sort]
+#wcount_naive = wcount_naive[index_naive]
+#wcount_sort = wcount_sort[index_sort]
+h_naive = h_naive[index_naive]
+h_sort = h_sort[index_sort]
+pos_naive = pos_naive[index_naive]
+pos_sort = pos_sort[index_sort]
+
+neighbour_length_naive = len(neighbour_ids_density_naive[0])
+neighbour_length_sort = len(neighbour_ids_density_sort[0])
+
+# Check that input files are logging the same number of neighbours
+if neighbour_length_naive != neighbour_length_sort:
+    print "Input files have logged different numbers of neighbour lengths!"
+    print "{} has logged: {} neighbours".format(inputFile1, neighbour_length_naive)
+    print "{} has logged: {} neighbours".format(inputFile2, neighbour_length_sort)
+    exit(1)
+
+if (max_density_ngbs_naive > neighbour_length_naive or max_force_ngbs_naive > neighbour_length_naive or
+    max_density_ngbs_sort > neighbour_length_sort or max_force_ngbs_sort > neighbour_length_sort):
+    print "The number of neighbours has exceeded the number of neighbours logged."
+    print "Modify NUM_OF_NEIGHBOURS in hydro_part.h to log more neighbours."
+    print "The highest neighbour count is: ", max(max_density_ngbs_naive,max_force_ngbs_naive, max_density_ngbs_sort,max_force_ngbs_sort)
+    exit(1)
+
+# First check
+print "\n                         Min    Max"
+print "                         ----------"
+print "Differences for density:  ", min(num_density_naive - num_density_sort), max(num_density_naive - num_density_sort)
+print "Differences for force:    ", min(num_force_naive - num_force_sort), max(num_force_naive - num_force_sort)
+
+# Get the IDs that are different
+mask_density = num_density_naive != num_density_sort
+mask_force = num_force_naive != num_force_sort
+num_invalid_density = np.sum(mask_density)
+num_invalid_force = np.sum(mask_force)
+
+print "\nNum non-zero density: ", num_invalid_density
+print "Num non-zero force:   ", num_invalid_force
+
+print "\nParticle IDs with incorrect densities"
+print "----------------------------------------"
+print ids_naive[mask_density]
+
+# Check density neighbour lists
+error += check_density_neighbours(ids_naive, neighbour_ids_density_naive,
+        neighbour_ids_density_sort, mask_density, pos_naive, h_naive,
+        num_invalid_density, 2e-6)
+
+print "Num of density interactions", inputFile1
+print num_density_naive[mask_density]
+
+print "Num of density interactions", inputFile2
+print num_density_sort[mask_density]
+
+print "\nParticle IDs with incorrect forces"
+print "------------------------------------"
+print ids_naive[mask_force]
+
+# Check force neighbour lists
+error += check_force_neighbours(ids_naive, neighbour_ids_force_naive,
+        neighbour_ids_force_sort, mask_force, pos_naive, h_naive,
+        num_invalid_force, 2e-6)
+
+print "Num of force interactions", inputFile1
+print num_force_naive[mask_force]
+
+#print "Smoothing lengths", inputFile1
+#print h_naive[mask_force]
+
+print "Num of force interactions", inputFile2
+print num_force_sort[mask_force]
+
+#print "Smoothing lengths", inputFile2
+#print h_sort[mask_force]
+
+# Statistics of h difference
+h_relative = (h_naive - h_sort) / h_naive
+print "h statistics: {} {} (Min, 1st Percentile)".format(np.min(h_relative), np.percentile(h_relative,1))
+print "h statistics: {} {} (Mean, Median)".format(np.mean(h_relative), np.median(h_relative))
+print "h statistics: {} {} (Max, 99th Percentile)".format(np.max(h_relative), np.percentile(h_relative, 99))
+
+if error:
+    print "\n------------------"
+    print "Differences found."
+    print "------------------"
+    exit(1)
+else:
+    print "\n---------------------"
+    print "No differences found."
+    print "---------------------"
+    exit(0)
diff --git a/src/debug.c b/src/debug.c
index 9042a5ebe6f1ba6b45e1c113219a816d8eafaffb..8675355ca3a97eaa29d79a7bc0e064bb577168d3 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -40,9 +40,7 @@
 #include "space.h"
 
 /* Import the right hydro definition */
-#if defined(DEBUG_INTERACTIONS_SPH)
-#include "./hydro/DebugInteractions/hydro_debug.h"
-#elif defined(MINIMAL_SPH)
+#if defined(MINIMAL_SPH)
 #include "./hydro/Minimal/hydro_debug.h"
 #elif defined(GADGET2_SPH)
 #include "./hydro/Gadget2/hydro_debug.h"
diff --git a/src/hydro.h b/src/hydro.h
index cbe26761e2470bea73028d1e2ed12204ac730d0e..abb49d35b204bbaf986f502d796883e7eb778e7f 100644
--- a/src/hydro.h
+++ b/src/hydro.h
@@ -27,11 +27,7 @@
 #include "kernel_hydro.h"
 
 /* Import the right functions */
-#if defined(DEBUG_INTERACTIONS_SPH)
-#include "./hydro/DebugInteractions/hydro.h"
-#include "./hydro/DebugInteractions/hydro_iact.h"
-#define SPH_IMPLEMENTATION "Debug SELF/PAIR"
-#elif defined(MINIMAL_SPH)
+#if defined(MINIMAL_SPH)
 #include "./hydro/Minimal/hydro.h"
 #include "./hydro/Minimal/hydro_iact.h"
 #define SPH_IMPLEMENTATION "Minimal version of SPH (e.g. Price 2010)"
diff --git a/src/hydro/DebugInteractions/hydro.h b/src/hydro/DebugInteractions/hydro.h
deleted file mode 100644
index 194c9e52bcc549f6633e34e571aee2a9012e7d0f..0000000000000000000000000000000000000000
--- a/src/hydro/DebugInteractions/hydro.h
+++ /dev/null
@@ -1,314 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Coypright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#ifndef SWIFT_DEBUG_INTERACTIONS_HYDRO_H
-#define SWIFT_DEBUG_INTERACTIONS_HYDRO_H
-
-/**
- * @file DebugInteractions/hydro.h
- * @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"
-
-#include <float.h>
-
-/**
- * @brief Returns the internal energy of a particle
- *
- * @param p The particle of interest
- */
-__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
-    const struct part *restrict p) {
-
-  return 0.f;
-}
-
-/**
- * @brief Returns the pressure of a particle
- *
- * @param p The particle of interest
- */
-__attribute__((always_inline)) INLINE static float hydro_get_pressure(
-    const struct part *restrict p) {
-
-  return 0.f;
-}
-
-/**
- * @brief Returns the entropy of a particle
- *
- * @param p The particle of interest
- */
-__attribute__((always_inline)) INLINE static float hydro_get_entropy(
-    const struct part *restrict p) {
-
-  return 0.f;
-}
-
-/**
- * @brief Returns the sound speed of a particle
- *
- * @param p The particle of interest
- */
-__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
-    const struct part *restrict p) {
-
-  return 0.f;
-}
-
-/**
- * @brief Returns the density of a particle
- *
- * @param p The particle of interest
- */
-__attribute__((always_inline)) INLINE static float hydro_get_density(
-    const struct part *restrict p) {
-
-  return 0.f;
-}
-
-/**
- * @brief Returns the mass of a particle
- *
- * @param p The particle of interest
- */
-__attribute__((always_inline)) INLINE static float hydro_get_mass(
-    const struct part *restrict p) {
-  return 0.f;
-}
-
-/**
- * @brief Returns the time derivative of internal energy of a particle
- *
- * We assume a constant density.
- *
- * @param p The particle of interest
- */
-__attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
-    const struct part *restrict p) {
-  return 0.f;
-}
-
-/**
- * @brief Returns the time derivative of internal energy of a particle
- *
- * We assume a constant density.
- *
- * @param p The particle of interest.
- * @param du_dt The new time derivative of the internal energy.
- */
-__attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
-    struct part *restrict p, float du_dt) {}
-
-/**
- * @brief Computes the hydro time-step of a given particle
- *
- * @param p Pointer to the particle data
- * @param xp Pointer to the extended particle data
- *
- */
-__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
-    const struct part *restrict p, const struct xpart *restrict xp,
-    const struct hydro_props *restrict hydro_properties) {
-
-  return FLT_MAX;
-}
-
-/**
- * @brief Does some extra hydro operations once the actual physical time step
- * for the particle is known.
- *
- * @param p The particle to act upon.
- * @param dt Physical time step of the particle during the next step.
- */
-__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
-    struct part *p, float dt) {}
-
-/**
- * @brief Prepares a particle for the density calculation.
- *
- * Zeroes all the relevant arrays in preparation for the sums taking place in
- * the variaous density tasks
- *
- * @param p The particle to act upon
- * @param hs #hydro_space containing hydro specific space information.
- */
-__attribute__((always_inline)) INLINE static void hydro_init_part(
-    struct part *restrict p, const struct hydro_space *hs) {
-
-  for (int i = 0; i < 256; ++i) p->ids_ngbs_density[i] = -1;
-  p->num_ngb_density = 0;
-
-  p->density.wcount = 0.f;
-  p->density.wcount_dh = 0.f;
-}
-
-/**
- * @brief Finishes the density calculation.
- *
- * Multiplies the density and number of neighbours by the appropiate constants
- * and add the self-contribution term.
- *
- * @param p The particle to act upon
- */
-__attribute__((always_inline)) INLINE static void hydro_end_density(
-    struct part *restrict p) {
-
-  /* Some smoothing length multiples. */
-  const float h = p->h;
-  const float h_inv = 1.0f / h;                       /* 1/h */
-  const float h_inv_dim = pow_dimension(h_inv);       /* 1/h^d */
-  const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
-
-  /* Final operation on the density (add self-contribution). */
-  p->density.wcount += kernel_root;
-  p->density.wcount_dh -= hydro_dimension * kernel_root;
-
-  /* Finish the calculation by inserting the missing h-factors */
-  p->density.wcount *= h_inv_dim;
-  p->density.wcount_dh *= h_inv_dim_plus_one;
-}
-
-/**
- * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
- *
- * @param p The particle to act upon
- * @param xp The extended particle data to act upon
- */
-__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
-    struct part *restrict p, struct xpart *restrict xp) {
-
-  /* Some smoothing length multiples. */
-  const float h = p->h;
-  const float h_inv = 1.0f / h;                 /* 1/h */
-  const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
-
-  /* Re-set problematic values */
-  p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
-  p->density.wcount_dh = 0.f;
-}
-
-/**
- * @brief Prepare a particle for the force calculation.
- *
- * Computes viscosity term, conduction term and smoothing length gradient terms.
- *
- * @param p The particle to act upon
- * @param xp The extended particle data to act upon
- * @param ti_current The current time (on the timeline)
- * @param timeBase The minimal time-step size
- */
-__attribute__((always_inline)) INLINE static void hydro_prepare_force(
-    struct part *restrict p, struct xpart *restrict xp) {}
-
-/**
- * @brief Reset acceleration fields of a particle
- *
- * Resets all hydro acceleration and time derivative fields in preparation
- * for the sums taking place in the variaous force tasks
- *
- * @param p The particle to act upon
- */
-__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
-    struct part *restrict p) {
-
-  for (int i = 0; i < 256; ++i) p->ids_ngbs_force[i] = -1;
-  p->num_ngb_force = 0;
-
-  p->force.h_dt = 0.f;
-}
-
-/**
- * @brief Sets the values to be predicted in the drifts to their values at a
- * kick time
- *
- * @param p The particle.
- * @param xp The extended data of this particle.
- */
-__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
-    struct part *restrict p, const struct xpart *restrict xp) {}
-
-/**
- * @brief Predict additional particle fields forward in time when drifting
- *
- * @param p The particle
- * @param xp The extended data of the particle
- * @param dt The drift time-step.
- * @param t0 The time at the start of the drift
- * @param t1 The time at the end of the drift
- * @param timeBase The minimal time-step size
- */
-__attribute__((always_inline)) INLINE static void hydro_predict_extra(
-    struct part *restrict p, const struct xpart *restrict xp, float dt) {}
-
-/**
- * @brief Finishes the force calculation.
- *
- * Multiplies the forces and accelerationsby the appropiate constants
- *
- * @param p The particle to act upon
- */
-__attribute__((always_inline)) INLINE static void hydro_end_force(
-    struct part *restrict p) {}
-
-/**
- * @brief Kick the additional variables
- *
- * @param p The particle to act upon
- * @param xp The particle extended data to act upon
- * @param dt The time-step for this kick
- */
-__attribute__((always_inline)) INLINE static void hydro_kick_extra(
-    struct part *restrict p, struct xpart *restrict xp, float dt) {}
-
-/**
- *  @brief Converts hydro quantity of a particle at the start of a run
- *
- * Requires the density to be known
- *
- * @param p The particle to act upon
- */
-__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
-    struct part *restrict p, struct xpart *restrict xp) {}
-
-/**
- * @brief Initialises the particles for the first time
- *
- * This function is called only once just after the ICs have been
- * read in to do some conversions.
- *
- * @param p The particle to act upon
- * @param xp The extended particle data to act upon
- */
-__attribute__((always_inline)) INLINE static void hydro_first_init_part(
-    struct part *restrict p, struct xpart *restrict xp) {
-
-  p->time_bin = 0;
-  xp->v_full[0] = p->v[0];
-  xp->v_full[1] = p->v[1];
-  xp->v_full[2] = p->v[2];
-
-  hydro_reset_acceleration(p);
-  hydro_init_part(p, NULL);
-}
-
-#endif /* SWIFT_DEBUG_INTERACTIONS_HYDRO_H */
diff --git a/src/hydro/DebugInteractions/hydro_debug.h b/src/hydro/DebugInteractions/hydro_debug.h
deleted file mode 100644
index 5420f9cece2415fb5e2438b551eb605187741c13..0000000000000000000000000000000000000000
--- a/src/hydro/DebugInteractions/hydro_debug.h
+++ /dev/null
@@ -1,39 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Coypright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#ifndef SWIFT_DEBUG_INTERACTIONS_HYDRO_DEBUG_H
-#define SWIFT_DEBUG_INTERACTIONS_HYDRO_DEBUG_H
-
-/**
- * @file DebugInteractions/hydro_debug.h
- * @brief Empty SPH implementation used solely to test the SELF/PAIR routines.
- */
-
-__attribute__((always_inline)) INLINE static void hydro_debug_particle(
-    const struct part* p, const struct xpart* xp) {
-  printf(
-      "x=[%.3e,%.3e,%.3e], "
-      "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
-      "h=%.3e, wcount=%.3f, wcount_dh=%.3e, dh/dt=%.3e time_bin=%d\n",
-      p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
-      xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
-      p->h, p->density.wcount, p->density.wcount_dh, p->force.h_dt,
-      p->time_bin);
-}
-
-#endif /* SWIFT_DEBUG_INTERACTIONS_HYDRO_DEBUG_H */
diff --git a/src/hydro/DebugInteractions/hydro_iact.h b/src/hydro/DebugInteractions/hydro_iact.h
deleted file mode 100644
index 89613f4902f0f9ab476110a7532664135b7b74ed..0000000000000000000000000000000000000000
--- a/src/hydro/DebugInteractions/hydro_iact.h
+++ /dev/null
@@ -1,113 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#ifndef SWIFT_DEBUG_INTERACTIONS_HYDRO_IACT_H
-#define SWIFT_DEBUG_INTERACTIONS_HYDRO_IACT_H
-
-/**
- * @file DebugInteractions/hydro_iact.h
- * @brief Empty SPH implementation used solely to test the SELF/PAIR routines.
- */
-
-/**
- * @brief Density loop
- */
-__attribute__((always_inline)) INLINE static void runner_iact_density(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
-
-  float wi, wi_dx;
-  float wj, wj_dx;
-
-  /* Get r and r inverse. */
-  const float r = sqrtf(r2);
-
-  /* Compute the kernel function for pi */
-  const float hi_inv = 1.f / hi;
-  const float ui = r * hi_inv;
-  kernel_deval(ui, &wi, &wi_dx);
-
-  /* Compute contribution to the number of neighbours */
-  pi->density.wcount += wi;
-  pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
-
-  /* Compute the kernel function for pj */
-  const float hj_inv = 1.f / hj;
-  const float uj = r * hj_inv;
-  kernel_deval(uj, &wj, &wj_dx);
-
-  /* Compute contribution to the number of neighbours */
-  pj->density.wcount += wj;
-  pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
-
-  /* Update ngb counters */
-  pi->ids_ngbs_density[pi->num_ngb_density] = pj->id;
-  ++pi->num_ngb_density;
-  pj->ids_ngbs_density[pj->num_ngb_density] = pi->id;
-  ++pj->num_ngb_density;
-}
-
-/**
- * @brief Density loop (non-symmetric version)
- */
-__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
-
-  float wi, wi_dx;
-
-  /* Get r and r inverse. */
-  const float r = sqrtf(r2);
-
-  /* Compute the kernel function */
-  const float hi_inv = 1.0f / hi;
-  const float ui = r * hi_inv;
-  kernel_deval(ui, &wi, &wi_dx);
-
-  /* Compute contribution to the number of neighbours */
-  pi->density.wcount += wi;
-  pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
-
-  /* Update ngb counters */
-  pi->ids_ngbs_density[pi->num_ngb_density] = pj->id;
-  ++pi->num_ngb_density;
-}
-
-/**
- * @brief Force loop
- */
-__attribute__((always_inline)) INLINE static void runner_iact_force(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
-
-  /* Update ngb counters */
-  pi->ids_ngbs_force[pi->num_ngb_force] = pj->id;
-  ++pi->num_ngb_force;
-  pj->ids_ngbs_force[pj->num_ngb_force] = pi->id;
-  ++pj->num_ngb_force;
-}
-
-/**
- * @brief Force loop (non-symmetric version)
- */
-__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
-    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
-
-  /* Update ngb counters */
-  pi->ids_ngbs_force[pi->num_ngb_force] = pj->id;
-  ++pi->num_ngb_force;
-}
-
-#endif /* SWIFT_DEBUG_INTERACTIONS_HYDRO_IACT_H */
diff --git a/src/hydro/DebugInteractions/hydro_io.h b/src/hydro/DebugInteractions/hydro_io.h
deleted file mode 100644
index 61b391a0c711edd5ad8acdbb6b914c070470992a..0000000000000000000000000000000000000000
--- a/src/hydro/DebugInteractions/hydro_io.h
+++ /dev/null
@@ -1,114 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Coypright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#ifndef SWIFT_DEBUG_INTERACTIONS_HYDRO_IO_H
-#define SWIFT_DEBUG_INTERACTIONS_HYDRO_IO_H
-
-/**
- * @file DebugInteractions/hydro_io.h
- * @brief Empty SPH implementation used solely to test the SELF/PAIR routines.
- */
-
-#include "adiabatic_index.h"
-#include "hydro.h"
-#include "io_properties.h"
-#include "kernel_hydro.h"
-
-/**
- * @brief Specifies which particle fields to read from a dataset
- *
- * @param parts The particle array.
- * @param list The list of i/o properties to read.
- * @param num_fields The number of i/o fields to read.
- */
-void hydro_read_particles(struct part* parts, struct io_props* list,
-                          int* num_fields) {
-
-  *num_fields = 4;
-
-  /* List what we want to read */
-  list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
-                                UNIT_CONV_LENGTH, parts, x);
-  list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
-                                UNIT_CONV_SPEED, parts, v);
-  list[2] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
-                                UNIT_CONV_LENGTH, parts, h);
-  list[3] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
-                                UNIT_CONV_NO_UNITS, parts, id);
-}
-
-float convert_u(struct engine* e, struct part* p) {
-
-  return hydro_get_internal_energy(p);
-}
-
-float convert_P(struct engine* e, struct part* p) {
-
-  return hydro_get_pressure(p);
-}
-
-/**
- * @brief Specifies which particle fields to write to a dataset
- *
- * @param parts The particle array.
- * @param list The list of i/o properties to write.
- * @param num_fields The number of i/o fields to write.
- */
-void hydro_write_particles(struct part* parts, struct io_props* list,
-                           int* num_fields) {
-
-  *num_fields = 8;
-
-  /* List what we want to write */
-  list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
-                                 parts, x);
-  list[1] =
-      io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v);
-  list[2] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
-                                 parts, h);
-  list[3] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
-                                 UNIT_CONV_NO_UNITS, parts, id);
-  list[4] = io_make_output_field("Num_ngb_density", INT, 1, UNIT_CONV_NO_UNITS,
-                                 parts, num_ngb_density);
-  list[5] = io_make_output_field("Num_ngb_force", INT, 1, UNIT_CONV_NO_UNITS,
-                                 parts, num_ngb_force);
-  list[6] = io_make_output_field("Ids_ngb_density", LONGLONG, 256,
-                                 UNIT_CONV_NO_UNITS, parts, ids_ngbs_density);
-  list[7] = io_make_output_field("Ids_ngb_force", LONGLONG, 256,
-                                 UNIT_CONV_NO_UNITS, parts, ids_ngbs_force);
-}
-
-/**
- * @brief Writes the current model of SPH to the file
- * @param h_grpsph The HDF5 group in which to write
- */
-void writeSPHflavour(hid_t h_grpsph) {
-
-  /* Viscosity and thermal conduction */
-  io_write_attribute_s(h_grpsph, "Thermal Conductivity Model", "None");
-  io_write_attribute_s(h_grpsph, "Viscosity Model", "None");
-}
-
-/**
- * @brief Are we writing entropy in the internal energy field ?
- *
- * @return 1 if entropy is in 'internal energy', 0 otherwise.
- */
-int writeEntropyFlag() { return 0; }
-
-#endif /* SWIFT_DEBUG_INTERACTIONS_HYDRO_IO_H */
diff --git a/src/hydro/DebugInteractions/hydro_part.h b/src/hydro/DebugInteractions/hydro_part.h
deleted file mode 100644
index 85410b2f3f9f11409398aeebc7bfe824c71cf217..0000000000000000000000000000000000000000
--- a/src/hydro/DebugInteractions/hydro_part.h
+++ /dev/null
@@ -1,110 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#ifndef SWIFT_DEBUG_INTERACTIONS_HYDRO_PART_H
-#define SWIFT_DEBUG_INTERACTIONS_HYDRO_PART_H
-
-/**
- * @file DebugInteractions/hydro_part.h
- * @brief Empty SPH implementation used solely to test the SELF/PAIR routines.
- */
-
-#include "cooling_struct.h"
-
-/* Extra particle data not needed during the SPH loops over neighbours. */
-struct xpart {
-
-  /* Offset between current position and position at last tree rebuild. */
-  float x_diff[3];
-
-  /* Offset between the current position and position at the last sort. */
-  float x_diff_sort[3];
-
-  /* Velocity at the last full step. */
-  float v_full[3];
-
-  /* Additional data used to record cooling information */
-  struct cooling_xpart_data cooling_data;
-
-} SWIFT_STRUCT_ALIGN;
-
-/* Data of a single particle. */
-struct part {
-
-  /* Particle ID. */
-  long long id;
-
-  /* Pointer to corresponding gravity part. */
-  struct gpart* gpart;
-
-  /* Particle position. */
-  double x[3];
-
-  /* Particle predicted velocity. */
-  float v[3];
-
-  /* Particle acceleration. */
-  float a_hydro[3];
-
-  /* Particle cutoff radius. */
-  float h;
-
-  struct {
-
-    /* Number of neighbours. */
-    float wcount;
-
-    /* Number of neighbours spatial derivative. */
-    float wcount_dh;
-
-  } density;
-
-  struct {
-
-    float h_dt;
-
-  } force;
-
-  /* Time-step length */
-  timebin_t time_bin;
-
-#ifdef SWIFT_DEBUG_CHECKS
-
-  /* Time of the last drift */
-  integertime_t ti_drift;
-
-  /* Time of the last kick */
-  integertime_t ti_kick;
-
-#endif
-
-  /*! List of interacting particles in the density SELF and PAIR */
-  long long ids_ngbs_density[256];
-
-  /*! List of interacting particles in the force SELF and PAIR */
-  long long ids_ngbs_force[256];
-
-  /*! Number of interactions in the density SELF and PAIR */
-  int num_ngb_density;
-
-  /*! Number of interactions in the force SELF and PAIR */
-  int num_ngb_force;
-
-} SWIFT_STRUCT_ALIGN;
-
-#endif /* SWIFT_DEBUG_INTERACTIONS_HYDRO_PART_H */
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 66a475f32ec06eb40ff2bc890bc156f76e3b7b9f..fdd0795cfcb084cd470403095b63d29b8f55c158 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -175,6 +175,11 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
 __attribute__((always_inline)) INLINE static void hydro_init_part(
     struct part *restrict p, const struct hydro_space *hs) {
 
+#ifdef DEBUG_INTERACTIONS_SPH
+  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS; ++i) p->ids_ngbs_density[i] = -1;
+  p->num_ngb_density = 0;
+#endif
+
   p->rho = 0.f;
   p->density.wcount = 0.f;
   p->density.wcount_dh = 0.f;
@@ -311,6 +316,11 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
 __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
     struct part *restrict p) {
 
+#ifdef DEBUG_INTERACTIONS_SPH
+  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS; ++i) p->ids_ngbs_force[i] = -1;
+  p->num_ngb_force = 0;
+#endif
+
   /* Reset the acceleration. */
   p->a_hydro[0] = 0.0f;
   p->a_hydro[1] = 0.0f;
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 666e0be835f114959f2f2a074900700b383e10fe..d30fd7b80318465cb654d769accd0e185336a878 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -103,6 +103,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   pj->density.rot_v[0] += facj * curlvr[0];
   pj->density.rot_v[1] += facj * curlvr[1];
   pj->density.rot_v[2] += facj * curlvr[2];
+
+#ifdef DEBUG_INTERACTIONS_SPH
+  /* Update ngb counters */
+  if (pi->num_ngb_density < MAX_NUM_OF_NEIGHBOURS)
+    pi->ids_ngbs_density[pi->num_ngb_density] = pj->id;
+  ++pi->num_ngb_density;
+
+  if (pj->num_ngb_density < MAX_NUM_OF_NEIGHBOURS)
+    pj->ids_ngbs_density[pj->num_ngb_density] = pi->id;
+  ++pj->num_ngb_density;
+#endif
 }
 
 /**
@@ -151,6 +162,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   pi->density.rot_v[0] += fac * curlvr[0];
   pi->density.rot_v[1] += fac * curlvr[1];
   pi->density.rot_v[2] += fac * curlvr[2];
+
+#ifdef DEBUG_INTERACTIONS_SPH
+  /* Update ngb counters */
+  if (pi->num_ngb_density < MAX_NUM_OF_NEIGHBOURS)
+    pi->ids_ngbs_density[pi->num_ngb_density] = pj->id;
+  ++pi->num_ngb_density;
+#endif
 }
 
 #ifdef WITH_VECTORIZATION
@@ -474,6 +492,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   /* Change in entropy */
   pi->entropy_dt += mj * visc_term * dvdr;
   pj->entropy_dt += mi * visc_term * dvdr;
+
+#ifdef DEBUG_INTERACTIONS_SPH
+  /* Update ngb counters */
+  if (pi->num_ngb_force < MAX_NUM_OF_NEIGHBOURS)
+    pi->ids_ngbs_force[pi->num_ngb_force] = pj->id;
+  ++pi->num_ngb_force;
+
+  if (pj->num_ngb_force < MAX_NUM_OF_NEIGHBOURS)
+    pj->ids_ngbs_force[pj->num_ngb_force] = pi->id;
+  ++pj->num_ngb_force;
+#endif
 }
 
 /**
@@ -563,6 +592,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Change in entropy */
   pi->entropy_dt += mj * visc_term * dvdr;
+
+#ifdef DEBUG_INTERACTIONS_SPH
+  /* Update ngb counters */
+  if (pi->num_ngb_force < MAX_NUM_OF_NEIGHBOURS)
+    pi->ids_ngbs_force[pi->num_ngb_force] = pj->id;
+  ++pi->num_ngb_force;
+#endif
 }
 
 #ifdef WITH_VECTORIZATION
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index d89fe5c277b0273e53304d77a241e14d2c71fc5d..34a9f33daffb6d087e41e22370a4475371459319 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -91,6 +91,10 @@ void hydro_write_particles(const struct part* parts, struct io_props* list,
 
   *num_fields = 10;
 
+#ifdef DEBUG_INTERACTIONS_SPH
+  *num_fields += 4;
+#endif
+
   /* List what we want to write */
   list[0] = io_make_output_field_convert_part(
       "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
@@ -113,6 +117,18 @@ void hydro_write_particles(const struct part* parts, struct io_props* list,
                                               parts, convert_u);
   list[9] = io_make_output_field_convert_part(
       "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, convert_P);
+#ifdef DEBUG_INTERACTIONS_SPH
+  list[10] = io_make_output_field("Num_ngb_density", INT, 1, UNIT_CONV_NO_UNITS,
+                                  parts, num_ngb_density);
+  list[11] = io_make_output_field("Num_ngb_force", INT, 1, UNIT_CONV_NO_UNITS,
+                                  parts, num_ngb_force);
+  list[12] =
+      io_make_output_field("Ids_ngb_density", LONGLONG, MAX_NUM_OF_NEIGHBOURS,
+                           UNIT_CONV_NO_UNITS, parts, ids_ngbs_density);
+  list[13] =
+      io_make_output_field("Ids_ngb_force", LONGLONG, MAX_NUM_OF_NEIGHBOURS,
+                           UNIT_CONV_NO_UNITS, parts, ids_ngbs_force);
+#endif
 }
 
 /**
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index 571aaf39ed2509d95e9f68c34ab8e6c09ded3fbb..c7deaf73608356c8a956d9130701baac1e740feb 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/src/hydro/Gadget2/hydro_part.h
@@ -143,6 +143,20 @@ struct part {
 
 #endif
 
+#ifdef DEBUG_INTERACTIONS_SPH
+  /*! List of interacting particles in the density SELF and PAIR */
+  long long ids_ngbs_density[MAX_NUM_OF_NEIGHBOURS];
+
+  /*! List of interacting particles in the force SELF and PAIR */
+  long long ids_ngbs_force[MAX_NUM_OF_NEIGHBOURS];
+
+  /*! Number of interactions in the density SELF and PAIR */
+  int num_ngb_density;
+
+  /*! Number of interactions in the force SELF and PAIR */
+  int num_ngb_force;
+#endif
+
 } SWIFT_STRUCT_ALIGN;
 
 #endif /* SWIFT_GADGET2_HYDRO_PART_H */
diff --git a/src/hydro_io.h b/src/hydro_io.h
index 29044de30960d1e80590d63609063c9552c79308..639c2f3ae640d7b74e6a2507bd4e3d5ad5625171 100644
--- a/src/hydro_io.h
+++ b/src/hydro_io.h
@@ -23,9 +23,7 @@
 #include "../config.h"
 
 /* Import the right functions */
-#if defined(DEBUG_INTERACTIONS_SPH)
-#include "./hydro/DebugInteractions/hydro_io.h"
-#elif defined(MINIMAL_SPH)
+#if defined(MINIMAL_SPH)
 #include "./hydro/Minimal/hydro_io.h"
 #elif defined(GADGET2_SPH)
 #include "./hydro/Gadget2/hydro_io.h"
diff --git a/src/part.h b/src/part.h
index ebe0ae41c6a292e5ec140dcaf32c4620c7f8586f..1b40aee0db3deb4790e07e3da9807060900d0c55 100644
--- a/src/part.h
+++ b/src/part.h
@@ -42,10 +42,7 @@
 #define gpart_align 128
 
 /* Import the right hydro particle definition */
-#if defined(DEBUG_INTERACTIONS_SPH)
-#include "./hydro/DebugInteractions/hydro_part.h"
-#define hydro_need_extra_init_loop 0
-#elif defined(MINIMAL_SPH)
+#if defined(MINIMAL_SPH)
 #include "./hydro/Minimal/hydro_part.h"
 #define hydro_need_extra_init_loop 0
 #elif defined(GADGET2_SPH)
diff --git a/src/runner.c b/src/runner.c
index c3d81a2cc829c753e235227e92880ead1c5e2a64..a8ea47e8ce449e716d5f454b4af4bd618dc4fcc9 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1935,24 +1935,15 @@ void *runner_main(void *data) {
       /* Different types of tasks... */
       switch (t->type) {
         case task_type_self:
-          if (t->subtype == task_subtype_density) {
-#if defined(WITH_VECTORIZATION) && defined(GADGET2_SPH)
-            runner_doself1_density_vec(r, ci);
-#else
-            runner_doself1_density(r, ci);
-#endif
-          }
+          if (t->subtype == task_subtype_density)
+            runner_doself1_branch_density(r, ci);
 #ifdef EXTRA_HYDRO_LOOP
           else if (t->subtype == task_subtype_gradient)
             runner_doself1_gradient(r, ci);
 #endif
-          else if (t->subtype == task_subtype_force) {
-#if defined(WITH_VECTORIZATION) && defined(GADGET2_SPH)
-            runner_doself2_force_vec(r, ci);
-#else
-            runner_doself2_force(r, ci);
-#endif
-          } else if (t->subtype == task_subtype_grav)
+          else if (t->subtype == task_subtype_force)
+            runner_doself2_branch_force(r, ci);
+          else if (t->subtype == task_subtype_grav)
             runner_doself_grav(r, ci, 1);
           else if (t->subtype == task_subtype_external_grav)
             runner_do_grav_external(r, ci, 1);
diff --git a/src/runner_doiact.h b/src/runner_doiact.h
index ef2bf1ab0a8b9e0621f1239dbb6ec6a943c5c4f5..5ec1457de9428bed73f26dab41d834489763b636 100644
--- a/src/runner_doiact.h
+++ b/src/runner_doiact.h
@@ -59,9 +59,15 @@
 #define _DOSELF2_NAIVE(f) PASTE(runner_doself2_naive, f)
 #define DOSELF2_NAIVE _DOSELF2_NAIVE(FUNCTION)
 
+#define _DOSELF1_BRANCH(f) PASTE(runner_doself1_branch, f)
+#define DOSELF1_BRANCH _DOSELF1_BRANCH(FUNCTION)
+
 #define _DOSELF1(f) PASTE(runner_doself1, f)
 #define DOSELF1 _DOSELF1(FUNCTION)
 
+#define _DOSELF2_BRANCH(f) PASTE(runner_doself2_branch, f)
+#define DOSELF2_BRANCH _DOSELF2_BRANCH(FUNCTION)
+
 #define _DOSELF2(f) PASTE(runner_doself2, f)
 #define DOSELF2 _DOSELF2(FUNCTION)
 
@@ -765,11 +771,6 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
 
   const struct engine *restrict e = r->e;
 
-#ifdef SWIFT_USE_NAIVE_INTERACTIONS
-  DOPAIR1_NAIVE(r, ci, cj);
-  return;
-#endif
-
   TIMER_TIC;
 
   /* Get the cutoff shift. */
@@ -1033,8 +1034,9 @@ void DOPAIR1_BRANCH(struct runner *r, struct cell *ci, struct cell *cj) {
   }
 #endif /* SWIFT_DEBUG_CHECKS */
 
-#if defined(WITH_VECTORIZATION) && defined(GADGET2_SPH) && \
-    (DOPAIR1_BRANCH == runner_dopair1_density_branch)
+#if defined(SWIFT_USE_NAIVE_INTERACTIONS)
+  DOPAIR1_NAIVE(r, ci, cj);
+#elif defined(WITH_VECTORIZATION) && defined(GADGET2_SPH)
   if (!sort_is_corner(sid))
     runner_dopair1_density_vec(r, ci, cj, sid, shift);
   else
@@ -1058,11 +1060,6 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
 
   struct engine *restrict e = r->e;
 
-#ifdef SWIFT_USE_NAIVE_INTERACTIONS
-  DOPAIR2_NAIVE(r, ci, cj);
-  return;
-#endif
-
   TIMER_TIC;
 
   /* Get the cutoff shift. */
@@ -1522,8 +1519,9 @@ void DOPAIR2_BRANCH(struct runner *r, struct cell *ci, struct cell *cj) {
   }
 #endif /* SWIFT_DEBUG_CHECKS */
 
-#if defined(WITH_VECTORIZATION) && defined(GADGET2_SPH) && \
-    (DOPAIR2_BRANCH == runner_dopair2_force_branch)
+#ifdef SWIFT_USE_NAIVE_INTERACTIONS
+  DOPAIR2_NAIVE(r, ci, cj);
+#elif defined(WITH_VECTORIZATION) && defined(GADGET2_SPH)
   if (!sort_is_corner(sid))
     runner_dopair2_force_vec(r, ci, cj, sid, shift);
   else
@@ -1543,17 +1541,8 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
 
   const struct engine *e = r->e;
 
-#ifdef SWIFT_USE_NAIVE_INTERACTIONS
-  DOSELF1_NAIVE(r, c);
-  return;
-#endif
-
   TIMER_TIC;
 
-  if (!cell_is_active_hydro(c, e)) return;
-
-  if (!cell_are_part_drifted(c, e)) error("Interacting undrifted cell.");
-
   struct part *restrict parts = c->parts;
   const int count = c->count;
 
@@ -1670,6 +1659,33 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
   TIMER_TOC(TIMER_DOSELF);
 }
 
+/**
+ * @brief Determine which version of DOSELF1 needs to be called depending on the
+ * optimisation level.
+ *
+ * @param r #runner
+ * @param c #cell c
+ *
+ */
+void DOSELF1_BRANCH(struct runner *r, struct cell *c) {
+
+  const struct engine *restrict e = r->e;
+
+  /* Anything to do here? */
+  if (!cell_is_active_hydro(c, e)) return;
+
+  /* Check that cells are drifted. */
+  if (!cell_are_part_drifted(c, e)) error("Interacting undrifted cell.");
+
+#if defined(SWIFT_USE_NAIVE_INTERACTIONS)
+  DOSELF1_NAIVE(r, c);
+#elif defined(WITH_VECTORIZATION) && defined(GADGET2_SPH)
+  runner_doself1_density_vec(r, c);
+#else
+  DOSELF1(r, c);
+#endif
+}
+
 /**
  * @brief Compute the cell self-interaction (symmetric).
  *
@@ -1680,16 +1696,8 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
 
   const struct engine *e = r->e;
 
-#ifdef SWIFT_USE_NAIVE_INTERACTIONS
-  DOSELF2_NAIVE(r, c);
-  return;
-#endif
-
   TIMER_TIC;
 
-  if (!cell_is_active_hydro(c, e)) return;
-  if (!cell_are_part_drifted(c, e)) error("Cell is not drifted");
-
   struct part *restrict parts = c->parts;
   const int count = c->count;
 
@@ -1798,6 +1806,33 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
   TIMER_TOC(TIMER_DOSELF);
 }
 
+/**
+ * @brief Determine which version of DOSELF2 needs to be called depending on the
+ * optimisation level.
+ *
+ * @param r #runner
+ * @param c #cell c
+ *
+ */
+void DOSELF2_BRANCH(struct runner *r, struct cell *c) {
+
+  const struct engine *restrict e = r->e;
+
+  /* Anything to do here? */
+  if (!cell_is_active_hydro(c, e)) return;
+
+  /* Check that cells are drifted. */
+  if (!cell_are_part_drifted(c, e)) error("Interacting undrifted cell.");
+
+#if defined(SWIFT_USE_NAIVE_INTERACTIONS)
+  DOSELF2_NAIVE(r, c);
+#elif defined(WITH_VECTORIZATION) && defined(GADGET2_SPH)
+  runner_doself2_force_vec(r, c);
+#else
+  DOSELF2(r, c);
+#endif
+}
+
 /**
  * @brief Compute grouped sub-cell interactions for pairs
  *
@@ -2082,12 +2117,7 @@ void DOSUB_SELF1(struct runner *r, struct cell *ci, int gettimer) {
     /* Drift the cell to the current timestep if needed. */
     if (!cell_are_part_drifted(ci, r->e)) error("Interacting undrifted cell.");
 
-#if (DOSELF1 == runner_doself1_density) && defined(WITH_VECTORIZATION) && \
-    defined(GADGET2_SPH)
-    runner_doself1_density_vec(r, ci);
-#else
-    DOSELF1(r, ci);
-#endif
+    DOSELF1_BRANCH(r, ci);
   }
 
   if (gettimer) TIMER_TOC(TIMER_DOSUB_SELF);
@@ -2374,12 +2404,7 @@ void DOSUB_SELF2(struct runner *r, struct cell *ci, int gettimer) {
 
   /* Otherwise, compute self-interaction. */
   else {
-#if (DOSELF2 == runner_doself2_force) && defined(WITH_VECTORIZATION) && \
-    defined(GADGET2_SPH)
-    runner_doself2_force_vec(r, ci);
-#else
-    DOSELF2(r, ci);
-#endif
+    DOSELF2_BRANCH(r, ci);
   }
   if (gettimer) TIMER_TOC(TIMER_DOSUB_SELF);
 }
diff --git a/src/runner_doiact_vec.c b/src/runner_doiact_vec.c
index 480fb6a3a4282703790adf87b2bc1872c394f2e2..05ad1ba7c1d4f171e993b870f11acbbe8707a3ef 100644
--- a/src/runner_doiact_vec.c
+++ b/src/runner_doiact_vec.c
@@ -664,6 +664,24 @@ __attribute__((always_inline)) INLINE void runner_doself1_density_vec(
       const int doi_mask2 = vec_is_mask_true(v_doi_mask2) &
                             vec_is_mask_true(v_doi_mask2_self_check);
 
+#ifdef DEBUG_INTERACTIONS_SPH
+      for (int bit_index = 0; bit_index < VEC_SIZE; bit_index++) {
+        if (doi_mask & (1 << bit_index)) {
+          if (pi->num_ngb_density < MAX_NUM_OF_NEIGHBOURS)
+            pi->ids_ngbs_density[pi->num_ngb_density] =
+                parts[pjd + bit_index].id;
+          ++pi->num_ngb_density;
+        }
+
+        if (doi_mask2 & (1 << bit_index)) {
+          if (pi->num_ngb_density < MAX_NUM_OF_NEIGHBOURS)
+            pi->ids_ngbs_density[pi->num_ngb_density] =
+                parts[pjd + VEC_SIZE + bit_index].id;
+          ++pi->num_ngb_density;
+        }
+      }
+#endif
+
       /* If there are any interactions left pack interaction values into c2
        * cache. */
       if (doi_mask) {
@@ -872,6 +890,24 @@ __attribute__((always_inline)) INLINE void runner_doself_subset_density_vec(
       const int doi_mask2 = vec_is_mask_true(v_doi_mask2) &
                             vec_is_mask_true(v_doi_mask2_self_check);
 
+#ifdef DEBUG_INTERACTIONS_SPH
+      for (int bit_index = 0; bit_index < VEC_SIZE; bit_index++) {
+        if (doi_mask & (1 << bit_index)) {
+          if (pi->num_ngb_density < MAX_NUM_OF_NEIGHBOURS)
+            pi->ids_ngbs_density[pi->num_ngb_density] =
+                parts[pjd + bit_index].id;
+          ++pi->num_ngb_density;
+        }
+
+        if (doi_mask2 & (1 << bit_index)) {
+          if (pi->num_ngb_density < MAX_NUM_OF_NEIGHBOURS)
+            pi->ids_ngbs_density[pi->num_ngb_density] =
+                parts[pjd + VEC_SIZE + bit_index].id;
+          ++pi->num_ngb_density;
+        }
+      }
+#endif
+
       /* If there are any interactions left pack interaction values into c2
        * cache. */
       if (doi_mask) {
@@ -1076,6 +1112,16 @@ __attribute__((always_inline)) INLINE void runner_doself2_force_vec(
       /* Combine all 3 masks. */
       vec_combine_masks(v_doi_mask, v_doi_mask_self_check);
 
+#ifdef DEBUG_INTERACTIONS_SPH
+      for (int bit_index = 0; bit_index < VEC_SIZE; bit_index++) {
+        if (vec_is_mask_true(v_doi_mask) & (1 << bit_index)) {
+          if (pi->num_ngb_force < MAX_NUM_OF_NEIGHBOURS)
+            pi->ids_ngbs_force[pi->num_ngb_force] = parts[pjd + bit_index].id;
+          ++pi->num_ngb_force;
+        }
+      }
+#endif
+
       /* If there are any interactions perform them. */
       if (vec_is_mask_true(v_doi_mask)) {
         vector v_hj_inv = vec_reciprocal(hj);
@@ -1316,6 +1362,17 @@ void runner_dopair1_density_vec(struct runner *r, struct cell *ci,
         /* Form r2 < hig2 mask. */
         vec_create_mask(v_doi_mask, vec_cmp_lt(v_r2.v, v_hig2.v));
 
+#ifdef DEBUG_INTERACTIONS_SPH
+        for (int bit_index = 0; bit_index < VEC_SIZE; bit_index++) {
+          if (vec_is_mask_true(v_doi_mask) & (1 << bit_index)) {
+            if (pi->num_ngb_density < MAX_NUM_OF_NEIGHBOURS)
+              pi->ids_ngbs_density[pi->num_ngb_density] =
+                  parts_j[sort_j[pjd + bit_index].i].id;
+            ++pi->num_ngb_density;
+          }
+        }
+#endif
+
         /* If there are any interactions perform them. */
         if (vec_is_mask_true(v_doi_mask))
           runner_iact_nonsym_1_vec_density(
@@ -1431,6 +1488,17 @@ void runner_dopair1_density_vec(struct runner *r, struct cell *ci,
         /* Form r2 < hig2 mask. */
         vec_create_mask(v_doj_mask, vec_cmp_lt(v_r2.v, v_hjg2.v));
 
+#ifdef DEBUG_INTERACTIONS_SPH
+        for (int bit_index = 0; bit_index < VEC_SIZE; bit_index++) {
+          if (vec_is_mask_true(v_doj_mask) & (1 << bit_index)) {
+            if (pj->num_ngb_density < MAX_NUM_OF_NEIGHBOURS)
+              pj->ids_ngbs_density[pj->num_ngb_density] =
+                  parts_i[sort_i[ci_cache_idx + first_pi + bit_index].i].id;
+            ++pj->num_ngb_density;
+          }
+        }
+#endif
+
         /* If there are any interactions perform them. */
         if (vec_is_mask_true(v_doj_mask))
           runner_iact_nonsym_1_vec_density(
@@ -1681,6 +1749,17 @@ void runner_dopair2_force_vec(struct runner *r, struct cell *ci,
         v_h2.v = vec_fmax(v_hig2.v, v_hjg2.v);
         vec_create_mask(v_doi_mask, vec_cmp_lt(v_r2.v, v_h2.v));
 
+#ifdef DEBUG_INTERACTIONS_SPH
+        for (int bit_index = 0; bit_index < VEC_SIZE; bit_index++) {
+          if (vec_is_mask_true(v_doi_mask) & (1 << bit_index)) {
+            if (pi->num_ngb_force < MAX_NUM_OF_NEIGHBOURS)
+              pi->ids_ngbs_force[pi->num_ngb_force] =
+                  parts_j[sort_j[pjd + bit_index].i].id;
+            ++pi->num_ngb_force;
+          }
+        }
+#endif
+
         /* If there are any interactions perform them. */
         if (vec_is_mask_true(v_doi_mask)) {
           vector v_hj_inv = vec_reciprocal(v_hj);
@@ -1806,6 +1885,17 @@ void runner_dopair2_force_vec(struct runner *r, struct cell *ci,
         v_h2.v = vec_fmax(v_hjg2.v, v_hig2.v);
         vec_create_mask(v_doj_mask, vec_cmp_lt(v_r2.v, v_h2.v));
 
+#ifdef DEBUG_INTERACTIONS_SPH
+        for (int bit_index = 0; bit_index < VEC_SIZE; bit_index++) {
+          if (vec_is_mask_true(v_doj_mask) & (1 << bit_index)) {
+            if (pj->num_ngb_force < MAX_NUM_OF_NEIGHBOURS)
+              pj->ids_ngbs_force[pj->num_ngb_force] =
+                  parts_i[sort_i[ci_cache_idx + first_pi + bit_index].i].id;
+            ++pj->num_ngb_force;
+          }
+        }
+#endif
+
         /* If there are any interactions perform them. */
         if (vec_is_mask_true(v_doj_mask)) {
           vector v_hi_inv = vec_reciprocal(v_hi);