diff --git a/src/cell.c b/src/cell.c
index 645dfbd3eaa94c4266243f19309c205824ce79f8..0f8694364bc806a5ccbc63f02a57103b5fee10f9 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -3173,7 +3173,7 @@ void cell_set_super_mapper(void *map_data, int num_elements, void *extra_data) {
 
     /* All top-level cells get an MPI tag. */
 #ifdef WITH_MPI
-    if (c->mpi.tag < 0 && c->mpi.sendto) cell_tag(c);
+    cell_ensure_tagged(c);
 #endif
 
     /* Super-pointer for hydro */
diff --git a/src/cell.h b/src/cell.h
index bc3a21dbbdbb1e03b6d51a4363e0153a750340f8..cfb5778f06d2d436f3844e33adf2bb31761aeaaa 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -850,20 +850,23 @@ __attribute__((always_inline)) INLINE static int cell_need_rebuild_for_pair(
 }
 
 /**
- * @brief Add a unique tag to a cell mostly for MPI communications.
+ * @brief Add a unique tag to a cell, mostly for MPI communications.
+ *
+ * This function locks the cell so that tags can be added concurrently.
  *
  * @param c The #cell to tag.
  */
-__attribute__((always_inline)) INLINE static void cell_tag(struct cell *c) {
+__attribute__((always_inline)) INLINE static void cell_ensure_tagged(
+    struct cell *c) {
 #ifdef WITH_MPI
 
-#ifdef SWIFT_DEBUG_CHECKS
-  if (c->mpi.tag > 0) error("setting tag for already tagged cell");
-#endif
-
+  lock_lock(&c->hydro.lock);
   if (c->mpi.tag < 0 &&
       (c->mpi.tag = atomic_inc(&cell_next_tag)) > cell_max_tag)
     error("Ran out of cell tags.");
+  if (lock_unlock(&c->hydro.lock) != 0) {
+    error("Failed to unlock cell.");
+  }
 #else
   error("SWIFT was not compiled with MPI enabled.");
 #endif  // WITH_MPI
diff --git a/src/engine.c b/src/engine.c
index 123c692edfbfbe6f6d8cef027998b43c0fa79c9d..de1e269ca4ada16750e33fa5f897d716f65a7953 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2116,7 +2116,7 @@ void engine_prepare(struct engine *e) {
   /* Unskip active tasks and check for rebuild */
   if (!e->forcerebuild && !e->forcerepart && !e->restarting) engine_unskip(e);
 
-  const ticks tic2 = getticks();
+  const ticks tic3 = getticks();
 
 #ifdef WITH_MPI
   MPI_Allreduce(MPI_IN_PLACE, &e->forcerebuild, 1, MPI_INT, MPI_MAX,
@@ -2125,7 +2125,7 @@ void engine_prepare(struct engine *e) {
 
   if (e->verbose)
     message("Communicating rebuild flag took %.3f %s.",
-            clocks_from_ticks(getticks() - tic2), clocks_getunit());
+            clocks_from_ticks(getticks() - tic3), clocks_getunit());
 
   /* Do we need repartitioning ? */
   if (e->forcerepart) {
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index 54a73513935cbd73db7362711bd992e669cbe7ca..dda69319608f06bb73de58fdfa4c41155ca6fd88 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -80,8 +80,8 @@ void engine_addtasks_send_gravity(struct engine *e, struct cell *ci,
     /* Create the tasks and their dependencies? */
     if (t_grav == NULL) {
 
-      /* Create a tag for this cell. */
-      if (ci->mpi.tag < 0) cell_tag(ci);
+      /* Make sure this cell is tagged. */
+      cell_ensure_tagged(ci);
 
       t_grav = scheduler_addtask(s, task_type_send, task_subtype_gpart,
                                  ci->mpi.tag, 0, ci, cj);
@@ -139,8 +139,8 @@ void engine_addtasks_send_hydro(struct engine *e, struct cell *ci,
     /* Create the tasks and their dependencies? */
     if (t_xv == NULL) {
 
-      /* Create a tag for this cell. */
-      if (ci->mpi.tag < 0) cell_tag(ci);
+      /* Make sure this cell is tagged. */
+      cell_ensure_tagged(ci);
 
       t_xv = scheduler_addtask(s, task_type_send, task_subtype_xv, ci->mpi.tag,
                                0, ci, cj);
@@ -238,8 +238,8 @@ void engine_addtasks_send_timestep(struct engine *e, struct cell *ci,
     /* Create the tasks and their dependencies? */
     if (t_ti == NULL) {
 
-      /* Create a tag for this cell. */
-      if (ci->mpi.tag < 0) cell_tag(ci);
+      /* Make sure this cell is tagged. */
+      cell_ensure_tagged(ci);
 
       t_ti = scheduler_addtask(s, task_type_send, task_subtype_tend,
                                ci->mpi.tag, 0, ci, cj);
@@ -1816,6 +1816,63 @@ void engine_make_hydroloop_tasks_mapper(void *map_data, int num_elements,
   }
 }
 
+struct cell_type_pair {
+  struct cell *ci, *cj;
+  int type;
+};
+
+void engine_addtasks_send_mapper(void *map_data, int num_elements,
+                                 void *extra_data) {
+  struct engine *e = (struct engine *)extra_data;
+  struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data;
+
+  for (int k = 0; k < num_elements; k++) {
+    struct cell *ci = cell_type_pairs[k].ci;
+    struct cell *cj = cell_type_pairs[k].cj;
+    const int type = cell_type_pairs[k].type;
+
+    /* Add the send task for the particle timesteps. */
+    engine_addtasks_send_timestep(e, ci, cj, NULL);
+
+    /* Add the send tasks for the cells in the proxy that have a hydro
+     * connection. */
+    if ((e->policy & engine_policy_hydro) && (type & proxy_cell_type_hydro))
+      engine_addtasks_send_hydro(e, ci, cj, /*t_xv=*/NULL,
+                                 /*t_rho=*/NULL, /*t_gradient=*/NULL);
+
+    /* Add the send tasks for the cells in the proxy that have a gravity
+     * connection. */
+    if ((e->policy & engine_policy_self_gravity) &&
+        (type & proxy_cell_type_gravity))
+      engine_addtasks_send_gravity(e, ci, cj, NULL);
+  }
+}
+
+void engine_addtasks_recv_mapper(void *map_data, int num_elements,
+                                 void *extra_data) {
+  struct engine *e = (struct engine *)extra_data;
+  struct cell_type_pair *cell_type_pairs = (struct cell_type_pair *)map_data;
+
+  for (int k = 0; k < num_elements; k++) {
+    struct cell *ci = cell_type_pairs[k].ci;
+    const int type = cell_type_pairs[k].type;
+
+    /* Add the recv task for the particle timesteps. */
+    engine_addtasks_recv_timestep(e, ci, NULL);
+
+    /* Add the recv tasks for the cells in the proxy that have a hydro
+     * connection. */
+    if ((e->policy & engine_policy_hydro) && (type & proxy_cell_type_hydro))
+      engine_addtasks_recv_hydro(e, ci, NULL, NULL, NULL);
+
+    /* Add the recv tasks for the cells in the proxy that have a gravity
+     * connection. */
+    if ((e->policy & engine_policy_self_gravity) &&
+        (type & proxy_cell_type_gravity))
+      engine_addtasks_recv_gravity(e, ci, NULL);
+  }
+}
+
 /**
  * @brief Fill the #space's task list.
  *
@@ -1978,7 +2035,7 @@ void engine_maketasks(struct engine *e) {
                    sched->nr_tasks, sizeof(struct task), 0, e);
 
   if (e->verbose)
-    message("Linking stars tasks took %.3f %s (including reweight).",
+    message("Linking stars tasks took %.3f %s.",
             clocks_from_ticks(getticks() - tic2), clocks_getunit());
 
 #ifdef WITH_MPI
@@ -1992,32 +2049,34 @@ void engine_maketasks(struct engine *e) {
 
     /* Loop over the proxies and add the send tasks, which also generates the
      * cell tags for super-cells. */
+    int max_num_send_cells = 0;
+    for (int pid = 0; pid < e->nr_proxies; pid++)
+      max_num_send_cells += e->proxies[pid].nr_cells_out;
+    struct cell_type_pair *send_cell_type_pairs = NULL;
+    if ((send_cell_type_pairs = (struct cell_type_pair *)malloc(
+             sizeof(struct cell_type_pair) * max_num_send_cells)) == NULL)
+      error("Failed to allocate temporary cell pointer list.");
+    int num_send_cells = 0;
+
     for (int pid = 0; pid < e->nr_proxies; pid++) {
 
       /* Get a handle on the proxy. */
       struct proxy *p = &e->proxies[pid];
 
-      for (int k = 0; k < p->nr_cells_out; k++)
-        engine_addtasks_send_timestep(e, p->cells_out[k], p->cells_in[0], NULL);
-
-      /* Loop through the proxy's outgoing cells and add the
-         send tasks for the cells in the proxy that have a hydro connection. */
-      if (e->policy & engine_policy_hydro)
-        for (int k = 0; k < p->nr_cells_out; k++)
-          if (p->cells_out_type[k] & proxy_cell_type_hydro)
-            engine_addtasks_send_hydro(e, p->cells_out[k], p->cells_in[0], NULL,
-                                       NULL, NULL);
-
-      /* Loop through the proxy's outgoing cells and add the
-         send tasks for the cells in the proxy that have a gravity connection.
-         */
-      if (e->policy & engine_policy_self_gravity)
-        for (int k = 0; k < p->nr_cells_out; k++)
-          if (p->cells_out_type[k] & proxy_cell_type_gravity)
-            engine_addtasks_send_gravity(e, p->cells_out[k], p->cells_in[0],
-                                         NULL);
+      for (int k = 0; k < p->nr_cells_out; k++) {
+        send_cell_type_pairs[num_send_cells].ci = p->cells_out[k];
+        send_cell_type_pairs[num_send_cells].cj = p->cells_in[0];
+        send_cell_type_pairs[num_send_cells++].type = p->cells_out_type[k];
+      }
     }
 
+    threadpool_map(&e->threadpool, engine_addtasks_send_mapper,
+                   send_cell_type_pairs, num_send_cells,
+                   sizeof(struct cell_type_pair),
+                   /*chunk=*/0, e);
+
+    free(send_cell_type_pairs);
+
     if (e->verbose)
       message("Creating send tasks took %.3f %s.",
               clocks_from_ticks(getticks() - tic2), clocks_getunit());
@@ -2035,29 +2094,28 @@ void engine_maketasks(struct engine *e) {
 
     /* Loop over the proxies and add the recv tasks, which relies on having the
      * cell tags. */
+    int max_num_recv_cells = 0;
+    for (int pid = 0; pid < e->nr_proxies; pid++)
+      max_num_recv_cells += e->proxies[pid].nr_cells_in;
+    struct cell_type_pair *recv_cell_type_pairs = NULL;
+    if ((recv_cell_type_pairs = (struct cell_type_pair *)malloc(
+             sizeof(struct cell_type_pair) * max_num_recv_cells)) == NULL)
+      error("Failed to allocate temporary cell pointer list.");
+    int num_recv_cells = 0;
     for (int pid = 0; pid < e->nr_proxies; pid++) {
 
       /* Get a handle on the proxy. */
       struct proxy *p = &e->proxies[pid];
-
-      for (int k = 0; k < p->nr_cells_in; k++)
-        engine_addtasks_recv_timestep(e, p->cells_in[k], NULL);
-
-      /* Loop through the proxy's incoming cells and add the
-         recv tasks for the cells in the proxy that have a hydro connection. */
-      if (e->policy & engine_policy_hydro)
-        for (int k = 0; k < p->nr_cells_in; k++)
-          if (p->cells_in_type[k] & proxy_cell_type_hydro)
-            engine_addtasks_recv_hydro(e, p->cells_in[k], NULL, NULL, NULL);
-
-      /* Loop through the proxy's incoming cells and add the
-         recv tasks for the cells in the proxy that have a gravity connection.
-         */
-      if (e->policy & engine_policy_self_gravity)
-        for (int k = 0; k < p->nr_cells_in; k++)
-          if (p->cells_in_type[k] & proxy_cell_type_gravity)
-            engine_addtasks_recv_gravity(e, p->cells_in[k], NULL);
+      for (int k = 0; k < p->nr_cells_in; k++) {
+        recv_cell_type_pairs[num_recv_cells].ci = p->cells_in[k];
+        recv_cell_type_pairs[num_recv_cells++].type = p->cells_in_type[k];
+      }
     }
+    threadpool_map(&e->threadpool, engine_addtasks_recv_mapper,
+                   recv_cell_type_pairs, num_recv_cells,
+                   sizeof(struct cell_type_pair),
+                   /*chunk=*/0, e);
+    free(recv_cell_type_pairs);
 
     if (e->verbose)
       message("Creating recv tasks took %.3f %s.",
diff --git a/src/lock.h b/src/lock.h
index b2dd2eac9d0ca5d7807907e31cf3fa31894f9aed..39601b0c52e414dad1a507b406c54640a254df30 100644
--- a/src/lock.h
+++ b/src/lock.h
@@ -34,6 +34,7 @@
 #define lock_trylock(l) (pthread_spin_lock(l) != 0)
 #define lock_unlock(l) (pthread_spin_unlock(l) != 0)
 #define lock_unlock_blind(l) pthread_spin_unlock(l)
+#define lock_static_initializer ((pthread_spinlock_t)0)
 
 #elif defined(PTHREAD_LOCK)
 #include <pthread.h>
@@ -44,6 +45,7 @@
 #define lock_trylock(l) (pthread_mutex_trylock(l) != 0)
 #define lock_unlock(l) (pthread_mutex_unlock(l) != 0)
 #define lock_unlock_blind(l) pthread_mutex_unlock(l)
+#define lock_static_initializer PTHREAD_MUTEX_INITIALIZER
 
 #else
 #define swift_lock_type volatile int
@@ -52,12 +54,12 @@
 INLINE static int lock_lock(volatile int *l) {
   while (atomic_cas(l, 0, 1) != 0)
     ;
-  // while( *l );
   return 0;
 }
 #define lock_trylock(l) ((*(l)) ? 1 : atomic_cas(l, 0, 1))
 #define lock_unlock(l) (atomic_cas(l, 1, 0) != 1)
 #define lock_unlock_blind(l) atomic_cas(l, 1, 0)
+#define lock_static_initializer 0
 #endif
 
 #endif /* SWIFT_LOCK_H */
diff --git a/src/proxy.c b/src/proxy.c
index 325ed78644b07a497374e40bfc8518edcb018593..4a67b4b3584c43b2df63f17303eba9ec5e742cb0 100644
--- a/src/proxy.c
+++ b/src/proxy.c
@@ -269,6 +269,93 @@ void proxy_cells_exchange_second(struct proxy *p) {
 #endif
 }
 
+#ifdef WITH_MPI
+
+void proxy_cells_count_mapper(void *map_data, int num_elements,
+                              void *extra_data) {
+  struct cell *cells = (struct cell *)map_data;
+
+  for (int k = 0; k < num_elements; k++) {
+    if (cells[k].mpi.sendto) cells[k].mpi.pcell_size = cell_getsize(&cells[k]);
+  }
+}
+
+struct pack_mapper_data {
+  struct space *s;
+  int *offset;
+  struct pcell *pcells;
+  int with_gravity;
+};
+
+void proxy_cells_pack_mapper(void *map_data, int num_elements,
+                             void *extra_data) {
+  struct cell *cells = (struct cell *)map_data;
+  struct pack_mapper_data *data = (struct pack_mapper_data *)extra_data;
+
+  for (int k = 0; k < num_elements; k++) {
+    if (cells[k].mpi.sendto) {
+      ptrdiff_t ind = &cells[k] - data->s->cells_top;
+      cells[k].mpi.pcell = &data->pcells[data->offset[ind]];
+      cell_pack(&cells[k], cells[k].mpi.pcell, data->with_gravity);
+    }
+  }
+}
+
+void proxy_cells_exchange_first_mapper(void *map_data, int num_elements,
+                                       void *extra_data) {
+  struct proxy *proxies = (struct proxy *)map_data;
+
+  for (int k = 0; k < num_elements; k++) {
+    proxy_cells_exchange_first(&proxies[k]);
+  }
+}
+
+struct wait_and_unpack_mapper_data {
+  struct space *s;
+  int num_proxies;
+  MPI_Request *reqs_in;
+  struct proxy *proxies;
+  int with_gravity;
+  swift_lock_type lock;
+};
+
+void proxy_cells_wait_and_unpack_mapper(void *unused_map_data, int num_elements,
+                                        void *extra_data) {
+
+  // MATTHIEU: This is currently unused. Scalar (non-threadpool) version is
+  // faster but we still need to explore why this happens.
+
+  struct wait_and_unpack_mapper_data *data =
+      (struct wait_and_unpack_mapper_data *)extra_data;
+
+  for (int k = 0; k < num_elements; k++) {
+    int pid = MPI_UNDEFINED;
+    MPI_Status status;
+    int res;
+
+    /* We need a lock to prevent concurrent calls to MPI_Waitany on
+       the same array of requests since this is not supported in the MPI
+       standard (v3.1). This is not really a problem since the threads
+       would block inside MPI_Waitany anyway. */
+    lock_lock(&data->lock);
+    if ((res = MPI_Waitany(data->num_proxies, data->reqs_in, &pid, &status)) !=
+            MPI_SUCCESS ||
+        pid == MPI_UNDEFINED)
+      mpi_error(res, "MPI_Waitany failed.");
+    if (lock_unlock(&data->lock) != 0) {
+      error("Failed to release lock.");
+    }
+
+    // message( "cell data from proxy %i has arrived." , pid );
+    for (int count = 0, j = 0; j < data->proxies[pid].nr_cells_in; j++)
+      count += cell_unpack(&data->proxies[pid].pcells_in[count],
+                           data->proxies[pid].cells_in[j], data->s,
+                           data->with_gravity);
+  }
+}
+
+#endif  // WITH_MPI
+
 /**
  * @brief Exchange the cell structures with all proxies.
  *
@@ -294,13 +381,14 @@ void proxy_cells_exchange(struct proxy *proxies, int num_proxies,
 
   /* Run through the cells and get the size of the ones that will be sent off.
    */
+  threadpool_map(&s->e->threadpool, proxy_cells_count_mapper, s->cells_top,
+                 s->nr_cells, sizeof(struct cell), /*chunk=*/0,
+                 /*extra_data=*/NULL);
   int count_out = 0;
   int offset[s->nr_cells];
   for (int k = 0; k < s->nr_cells; k++) {
     offset[k] = count_out;
-    if (s->cells_top[k].mpi.sendto)
-      count_out +=
-          (s->cells_top[k].mpi.pcell_size = cell_getsize(&s->cells_top[k]));
+    if (s->cells_top[k].mpi.sendto) count_out += s->cells_top[k].mpi.pcell_size;
   }
 
   if (s->e->verbose)
@@ -316,19 +404,19 @@ void proxy_cells_exchange(struct proxy *proxies, int num_proxies,
   tic2 = getticks();
 
   /* Pack the cells. */
-  for (int k = 0; k < s->nr_cells; k++)
-    if (s->cells_top[k].mpi.sendto) {
-      cell_pack(&s->cells_top[k], &pcells[offset[k]], with_gravity);
-      s->cells_top[k].mpi.pcell = &pcells[offset[k]];
-    }
+  struct pack_mapper_data data = {s, offset, pcells, with_gravity};
+  threadpool_map(&s->e->threadpool, proxy_cells_pack_mapper, s->cells_top,
+                 s->nr_cells, sizeof(struct cell), /*chunk=*/0, &data);
 
   if (s->e->verbose)
     message("Packing cells took %.3f %s.", clocks_from_ticks(getticks() - tic2),
             clocks_getunit());
 
   /* Launch the first part of the exchange. */
+  threadpool_map(&s->e->threadpool, proxy_cells_exchange_first_mapper, proxies,
+                 num_proxies, sizeof(struct proxy), /*chunk=*/0,
+                 /*extra_data=*/NULL);
   for (int k = 0; k < num_proxies; k++) {
-    proxy_cells_exchange_first(&proxies[k]);
     reqs_in[k] = proxies[k].req_cells_count_in;
     reqs_out[k] = proxies[k].req_cells_count_out;
   }
diff --git a/src/runner.c b/src/runner.c
index c7e186f5cf6263199c26fb32fbe745fb48b56c08..95ffa5ca5e4cbd68a032703c6c481952496865ab 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -337,7 +337,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
     free(sid);
   }
 
-  if (timer) TIMER_TOC(timer_do_stars_ghost);
+  if (timer) TIMER_TOC(timer_dostars_ghost);
 }
 
 /**
diff --git a/src/runner_doiact_stars.h b/src/runner_doiact_stars.h
index e696e4fd10853536008d9d9fafc90e6475fd291a..7d2959da717b0b5dafd688245d32e0e611e169e5 100644
--- a/src/runner_doiact_stars.h
+++ b/src/runner_doiact_stars.h
@@ -242,7 +242,7 @@ void runner_dopair_subset_stars_density(struct runner *r,
     } /* loop over the parts in cj. */
   }   /* loop over the parts in ci. */
 
-  TIMER_TOC(timer_dopair_subset_naive);
+  TIMER_TOC(timer_dopair_subset_stars_density);
 }
 
 /**
@@ -944,7 +944,7 @@ void runner_dosub_subset_stars_density(struct runner *r, struct cell *ci,
 
   } /* otherwise, pair interaction. */
 
-  if (gettimer) TIMER_TOC(timer_dosub_subset);
+  if (gettimer) TIMER_TOC(timer_dosub_subset_pair_stars_density);
 }
 
 /**
@@ -1371,7 +1371,7 @@ void runner_dosub_pair_stars_density(struct runner *r, struct cell *ci,
     runner_dopair_branch_stars_density(r, ci, cj);
   }
 
-  if (gettimer) TIMER_TOC(TIMER_DOSUB_PAIR);
+  if (gettimer) TIMER_TOC(timer_dosub_pair_stars_density);
 }
 
 /**
diff --git a/src/timers.h b/src/timers.h
index c43f0154d2aaaa9d1c5aed3ad51b912e4fc5d751..52dc076ceb17f31c4e6c6515250faa46c4a9761e 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -93,8 +93,10 @@ enum {
   timer_dostars_ghost,
   timer_doself_subset_stars_density,
   timer_dopair_subset_stars_density,
-  timer_dosubpair_stars_density,
+  timer_dosub_pair_stars_density,
   timer_dosub_self_stars_density,
+  timer_dosub_subset_pair_stars_density,
+  timer_dosub_subset_self_stars_density,
   timer_logger,
   timer_count,
 };
diff --git a/src/tools.c b/src/tools.c
index ff33afbca3fdc97ede61ff09998d8dc27bf0154b..c0400aa7b42322fce276a5e788af7bcb9e7f3625 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -714,7 +714,7 @@ int compare_values(double a, double b, double threshold, double *absDiff,
  *
  * @return 1 if difference found, 0 otherwise
  */
-int compare_particles(struct part a, struct part b, double threshold) {
+int compare_particles(struct part *a, struct part *b, double threshold) {
 
 #ifdef GADGET2_SPH
 
@@ -722,117 +722,117 @@ int compare_particles(struct part a, struct part b, double threshold) {
   double absDiff = 0.0, absSum = 0.0, relDiff = 0.0;
 
   for (int k = 0; k < 3; k++) {
-    if (compare_values(a.x[k], b.x[k], threshold, &absDiff, &absSum,
+    if (compare_values(a->x[k], b->x[k], threshold, &absDiff, &absSum,
                        &relDiff)) {
       message(
           "Relative difference (%e) larger than tolerance (%e) for x[%d] of "
           "particle %lld.",
-          relDiff, threshold, k, a.id);
-      message("a = %e, b = %e", a.x[k], b.x[k]);
+          relDiff, threshold, k, a->id);
+      message("a = %e, b = %e", a->x[k], b->x[k]);
       result = 1;
     }
   }
   for (int k = 0; k < 3; k++) {
-    if (compare_values(a.v[k], b.v[k], threshold, &absDiff, &absSum,
+    if (compare_values(a->v[k], b->v[k], threshold, &absDiff, &absSum,
                        &relDiff)) {
       message(
           "Relative difference (%e) larger than tolerance (%e) for v[%d] of "
           "particle %lld.",
-          relDiff, threshold, k, a.id);
-      message("a = %e, b = %e", a.v[k], b.v[k]);
+          relDiff, threshold, k, a->id);
+      message("a = %e, b = %e", a->v[k], b->v[k]);
       result = 1;
     }
   }
   for (int k = 0; k < 3; k++) {
-    if (compare_values(a.a_hydro[k], b.a_hydro[k], threshold, &absDiff, &absSum,
-                       &relDiff)) {
+    if (compare_values(a->a_hydro[k], b->a_hydro[k], threshold, &absDiff,
+                       &absSum, &relDiff)) {
       message(
           "Relative difference (%e) larger than tolerance (%e) for a_hydro[%d] "
           "of particle %lld.",
-          relDiff, threshold, k, a.id);
-      message("a = %e, b = %e", a.a_hydro[k], b.a_hydro[k]);
+          relDiff, threshold, k, a->id);
+      message("a = %e, b = %e", a->a_hydro[k], b->a_hydro[k]);
       result = 1;
     }
   }
-  if (compare_values(a.rho, b.rho, threshold, &absDiff, &absSum, &relDiff)) {
+  if (compare_values(a->rho, b->rho, threshold, &absDiff, &absSum, &relDiff)) {
     message(
         "Relative difference (%e) larger than tolerance (%e) for rho of "
         "particle %lld.",
-        relDiff, threshold, a.id);
-    message("a = %e, b = %e", a.rho, b.rho);
+        relDiff, threshold, a->id);
+    message("a = %e, b = %e", a->rho, b->rho);
     result = 1;
   }
-  if (compare_values(a.density.rho_dh, b.density.rho_dh, threshold, &absDiff,
+  if (compare_values(a->density.rho_dh, b->density.rho_dh, threshold, &absDiff,
                      &absSum, &relDiff)) {
     message(
         "Relative difference (%e) larger than tolerance (%e) for rho_dh of "
         "particle %lld.",
-        relDiff, threshold, a.id);
-    message("a = %e, b = %e", a.density.rho_dh, b.density.rho_dh);
+        relDiff, threshold, a->id);
+    message("a = %e, b = %e", a->density.rho_dh, b->density.rho_dh);
     result = 1;
   }
-  if (compare_values(a.density.wcount, b.density.wcount, threshold, &absDiff,
+  if (compare_values(a->density.wcount, b->density.wcount, threshold, &absDiff,
                      &absSum, &relDiff)) {
     message(
         "Relative difference (%e) larger than tolerance (%e) for wcount of "
         "particle %lld.",
-        relDiff, threshold, a.id);
-    message("a = %e, b = %e", a.density.wcount, b.density.wcount);
+        relDiff, threshold, a->id);
+    message("a = %e, b = %e", a->density.wcount, b->density.wcount);
     result = 1;
   }
-  if (compare_values(a.density.wcount_dh, b.density.wcount_dh, threshold,
+  if (compare_values(a->density.wcount_dh, b->density.wcount_dh, threshold,
                      &absDiff, &absSum, &relDiff)) {
     message(
         "Relative difference (%e) larger than tolerance (%e) for wcount_dh of "
         "particle %lld.",
-        relDiff, threshold, a.id);
-    message("a = %e, b = %e", a.density.wcount_dh, b.density.wcount_dh);
+        relDiff, threshold, a->id);
+    message("a = %e, b = %e", a->density.wcount_dh, b->density.wcount_dh);
     result = 1;
   }
-  if (compare_values(a.force.h_dt, b.force.h_dt, threshold, &absDiff, &absSum,
+  if (compare_values(a->force.h_dt, b->force.h_dt, threshold, &absDiff, &absSum,
                      &relDiff)) {
     message(
         "Relative difference (%e) larger than tolerance (%e) for h_dt of "
         "particle %lld.",
-        relDiff, threshold, a.id);
-    message("a = %e, b = %e", a.force.h_dt, b.force.h_dt);
+        relDiff, threshold, a->id);
+    message("a = %e, b = %e", a->force.h_dt, b->force.h_dt);
     result = 1;
   }
-  if (compare_values(a.force.v_sig, b.force.v_sig, threshold, &absDiff, &absSum,
-                     &relDiff)) {
+  if (compare_values(a->force.v_sig, b->force.v_sig, threshold, &absDiff,
+                     &absSum, &relDiff)) {
     message(
         "Relative difference (%e) larger than tolerance (%e) for v_sig of "
         "particle %lld.",
-        relDiff, threshold, a.id);
-    message("a = %e, b = %e", a.force.v_sig, b.force.v_sig);
+        relDiff, threshold, a->id);
+    message("a = %e, b = %e", a->force.v_sig, b->force.v_sig);
     result = 1;
   }
-  if (compare_values(a.entropy_dt, b.entropy_dt, threshold, &absDiff, &absSum,
+  if (compare_values(a->entropy_dt, b->entropy_dt, threshold, &absDiff, &absSum,
                      &relDiff)) {
     message(
         "Relative difference (%e) larger than tolerance (%e) for entropy_dt of "
         "particle %lld.",
-        relDiff, threshold, a.id);
-    message("a = %e, b = %e", a.entropy_dt, b.entropy_dt);
+        relDiff, threshold, a->id);
+    message("a = %e, b = %e", a->entropy_dt, b->entropy_dt);
     result = 1;
   }
-  if (compare_values(a.density.div_v, b.density.div_v, threshold, &absDiff,
+  if (compare_values(a->density.div_v, b->density.div_v, threshold, &absDiff,
                      &absSum, &relDiff)) {
     message(
         "Relative difference (%e) larger than tolerance (%e) for div_v of "
         "particle %lld.",
-        relDiff, threshold, a.id);
-    message("a = %e, b = %e", a.density.div_v, b.density.div_v);
+        relDiff, threshold, a->id);
+    message("a = %e, b = %e", a->density.div_v, b->density.div_v);
     result = 1;
   }
   for (int k = 0; k < 3; k++) {
-    if (compare_values(a.density.rot_v[k], b.density.rot_v[k], threshold,
+    if (compare_values(a->density.rot_v[k], b->density.rot_v[k], threshold,
                        &absDiff, &absSum, &relDiff)) {
       message(
           "Relative difference (%e) larger than tolerance (%e) for rot_v[%d] "
           "of particle %lld.",
-          relDiff, threshold, k, a.id);
-      message("a = %e, b = %e", a.density.rot_v[k], b.density.rot_v[k]);
+          relDiff, threshold, k, a->id);
+      message("a = %e, b = %e", a->density.rot_v[k], b->density.rot_v[k]);
       result = 1;
     }
   }
diff --git a/src/tools.h b/src/tools.h
index 5ec2ca42810e1cb733e1aa674bd1c94ac5c569bd..22eba1519ebf80673cb2d8540791e6b4d092bab0 100644
--- a/src/tools.h
+++ b/src/tools.h
@@ -54,7 +54,7 @@ void gravity_n2(struct gpart *gparts, const int gcount,
                 const struct gravity_props *gravity_properties, float rlr);
 int compare_values(double a, double b, double threshold, double *absDiff,
                    double *absSum, double *relDiff);
-int compare_particles(struct part a, struct part b, double threshold);
+int compare_particles(struct part *a, struct part *b, double threshold);
 
 long get_maxrss(void);
 
diff --git a/tests/testInteractions.c b/tests/testInteractions.c
index 0241a46634bdd12fb8a89bc3912c3b160198c389..dae55a337642e1616e94119263ff8f1c2a617c89 100644
--- a/tests/testInteractions.c
+++ b/tests/testInteractions.c
@@ -197,10 +197,10 @@ int check_results(struct part serial_test_part, struct part *serial_parts,
                   struct part vec_test_part, struct part *vec_parts,
                   int count) {
   int result = 0;
-  result += compare_particles(serial_test_part, vec_test_part, ACC_THRESHOLD);
+  result += compare_particles(&serial_test_part, &vec_test_part, ACC_THRESHOLD);
 
   for (int i = 0; i < count; i++)
-    result += compare_particles(serial_parts[i], vec_parts[i], ACC_THRESHOLD);
+    result += compare_particles(&serial_parts[i], &vec_parts[i], ACC_THRESHOLD);
 
   return result;
 }
diff --git a/tests/testPeriodicBC.c b/tests/testPeriodicBC.c
index e4a5a85845f92bfbecae85d9dc9453b38b9b7646..a16710d9117037a3bd6b51dcb03b02b4fa821509 100644
--- a/tests/testPeriodicBC.c
+++ b/tests/testPeriodicBC.c
@@ -273,7 +273,7 @@ int check_results(struct part *serial_parts, struct part *vec_parts, int count,
   int result = 0;
 
   for (int i = 0; i < count; i++)
-    result += compare_particles(serial_parts[i], vec_parts[i], threshold);
+    result += compare_particles(&serial_parts[i], &vec_parts[i], threshold);
 
   return result;
 }