diff --git a/examples/plot.py b/examples/plot.py
index 6f395138d81a2d7448eb9afe045d33496a03a571..f89c44525014c645f4a0f9cf164c9fb3119c4d7a 100644
--- a/examples/plot.py
+++ b/examples/plot.py
@@ -64,7 +64,6 @@ accx_new = accx_new[rank]
 accy_new = accy_new[rank]
 accz_new = accz_new[rank]
 
-
 # # Read Gadget accelerations
 # data=loadtxt("particle_dump_gadget.dat")
 # id = data[:,0]
diff --git a/examples/test_bh_3.c b/examples/test_bh_3.c
index 693a1b51fc4131702fb1d8b4ef32740c2096a4cc..820ffe03399c1cacfa88dfc3e47965b20831b1cb 100644
--- a/examples/test_bh_3.c
+++ b/examples/test_bh_3.c
@@ -815,6 +815,7 @@ void legacy_comp_com(struct cell *c, int *countCoMs) {
     c->legacy.com[2] = 0.0;
     c->legacy.mass = 0.0;
   }
+
 }
 
 /**
@@ -968,14 +969,13 @@ void interact_exact(int N, struct part *parts, int monitor) {
     /* Some things to store locally. */
     double pix[3] = {parts[i].x[0], parts[i].x[1], parts[i].x[2]};
     double mi = parts[i].mass;
-    double ai[3] = {0.0, 0.0, 0.0};
 
     /* Loop over every other particle. */
     for (j = i + 1; j < N; ++j) {
 
       /* Compute the pairwise distance. */
       for (r2 = 0.0, k = 0; k < 3; k++) {
-        dx[k] = parts[i].x[k] - pix[k];
+        dx[k] = parts[j].x[k] - pix[k];
         r2 += dx[k] * dx[k];
       }
 
@@ -985,15 +985,11 @@ void interact_exact(int N, struct part *parts, int monitor) {
 
       for (k = 0; k < 3; k++) {
         double wdx = w * dx[k];
-        parts[j].a_exact[k] += wdx * mi;
-        ai[k] -= wdx * parts[j].mass;
+        parts[j].a_exact[k] -= wdx * mi;
+        parts[i].a_exact[k] += wdx * parts[j].mass;
       }
     }
 
-    /* Store the particle acceleration. */
-    for (k = 0; k < 3; k++) {
-      parts[i].a_exact[k] = ai[k];
-    }
   }
 
   for (i = 0; i < N; ++i)
@@ -1098,8 +1094,6 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
   root->loc[1] = 0.0;
   root->loc[2] = 0.0;
   root->h = 1.0;
-  root->h = 1.0;
-  root->h = 1.0;
   root->count = N;
   root->parts = parts;
   cell_split(root, &s);
@@ -1109,7 +1103,7 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
   /* Do a N^2 interactions calculation */
 
   tic_exact = getticks();
-  // interact_exact( N , parts , ICHECK );
+  interact_exact( N , parts , ICHECK );
   toc_exact = getticks();
 
   printf("Exact calculation (1 thread) took %lli (= %e) ticks\n",
@@ -1210,7 +1204,7 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
 
   /* Dump the particles to a file */
   file = fopen("particle_dump.dat", "w");
-  fprintf(file, "# a_exact.x   a_exact.y    a_exact.z    a_legacy.x    "
+  fprintf(file, "# x y z a_exact.x   a_exact.y    a_exact.z    a_legacy.x    "
                 "a_legacy.y    a_legacy.z    a_new.x     a_new.y    a_new.z\n");
   for (k = 0; k < N; ++k)
     fprintf(file, "%d %e %e %e %e %e %e %e %e %e %e %e %e\n", parts[k].id,