diff --git a/src/cell.c b/src/cell.c
index 0fe4a64d8e19391c49ab91695a78c20b4e11ffbc..d12007c3c339268a9d4d23d939f5d37451e3d890 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -4010,6 +4010,10 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
       /* Maximal smoothing length */
       cell_h_max = max(cell_h_max, sp->h);
 
+      /* Get ready for a density calculation */
+      if (spart_is_active(sp, e)) {
+        stars_init_spart(sp);
+      }
     } /* Note: no need to compute dx_max as all spart have a gpart */
 
     /* Now, get the maximal particle motion from its square */
diff --git a/src/engine.c b/src/engine.c
index f32bc0dbe7e2346f9d15984575f1128a0b589d21..54d62d2389bf6b21298f0d5c95d85d7e865337a0 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3900,6 +3900,73 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
 #endif
 }
 
+#ifdef DEBUG_INTERACTIONS_STARS
+/**
+ * @brief Exchange the feedback counters between stars
+ * @param e The #engine.
+ */
+void engine_collect_stars_counter(struct engine *e) {
+
+#ifdef WITH_MPI
+  if (e->total_nr_sparts > 1e5) {
+    message("WARNING: too many sparts, skipping exchange of counters");
+    return;
+  }
+
+  /* Get number of particles for each rank */
+  size_t *n_parts = (size_t *)malloc(e->nr_nodes * sizeof(size_t));
+
+#define MPI_size MPI_UNSIGNED_LONG
+  int err = MPI_Allgather(&e->s->nr_sparts_foreign, 1, MPI_size, n_parts, 1, MPI_size, MPI_COMM_WORLD);
+  if (err != MPI_SUCCESS)
+    error("Communication failed");
+
+  /* Compute derivated quantities */
+  int total = 0;
+  int *n_parts_int = (int *) malloc(e->nr_nodes * sizeof(int));
+  int *displs = (int *) malloc(e->nr_nodes * sizeof(int));
+  for(int i = 0; i < e->nr_nodes; i++) {
+    displs[i] = total;
+    total += n_parts[i];
+    n_parts_int[i] = n_parts[i];
+  }
+
+  /* Get all particles */
+  struct spart *parts = (struct spart *) malloc(total * sizeof(struct spart));
+  err = MPI_Allgatherv(e->s->sparts_foreign, e->s->nr_sparts_foreign, spart_mpi_type,
+		       parts, n_parts_int, displs, spart_mpi_type, MPI_COMM_WORLD);
+  if (err != MPI_SUCCESS)
+    error("Communication failed");
+
+  /* Reset counters */
+  for(size_t i = 0; i < e->s->nr_sparts_foreign; i++) {
+    e->s->sparts_foreign[i].num_ngb_force = 0;
+  }
+
+  /* Update counters */
+  struct spart *local_parts = e->s->sparts;
+  for(int i = 0; i < total; i++) {
+    /* Avoids counting twice local interactions */
+    if (i >= displs[engine_rank] &&
+	i < displs[engine_rank] + n_parts_int[engine_rank])
+      continue;
+
+    const long long id_i = parts[i].id;
+
+    for(size_t j = 0; j < e->s->nr_sparts; j++) {
+      const long long id_j = local_parts[j].id;
+      if (id_j == id_i && parts[i].num_ngb_force != 0) {
+	local_parts[j].num_ngb_force += parts[i].num_ngb_force;
+	break;
+      }
+    }
+  }
+#endif  
+}
+
+#endif
+
+
 /**
  * @brief Writes a snapshot with the current state of the engine
  *
@@ -3936,6 +4003,10 @@ void engine_dump_snapshot(struct engine *e) {
   }
 #endif
 
+#ifdef DEBUG_INTERACTIONS_STARS
+  engine_collect_stars_counter(e);
+#endif
+  
 /* Dump... */
 #if defined(HAVE_HDF5)
 #if defined(WITH_MPI)
diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c
index 68754276e1917de50d2dc1bf80e5fef159041265..d9b9a234873bf00d77ffb38e5f44a73c4549fb9f 100644
--- a/src/engine_marktasks.c
+++ b/src/engine_marktasks.c
@@ -358,6 +358,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
 
 	const int should_do = t_subtype == task_subtype_stars_density ||
 	  cj->nodeID != ci->nodeID;
+
         /* Set the correct sorting flags */
         if (t_type == task_type_pair) {
 
diff --git a/src/runner.c b/src/runner.c
index d3f7635d520ef9f8738eafbf9986f4ff254fd807..4856cb7b7f09a8a83b20175b9eabdefd97dae343 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -2882,7 +2882,7 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int clear_sorts,
 
 #ifdef WITH_MPI
 
-  const struct spart *restrict sparts = c->stars.parts;
+  struct spart *restrict sparts = c->stars.parts;
   const size_t nr_sparts = c->stars.count;
   const integertime_t ti_current = r->e->ti_current;
 
@@ -2907,6 +2907,9 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int clear_sorts,
 
     /* Collect everything... */
     for (size_t k = 0; k < nr_sparts; k++) {
+#ifdef DEBUG_INTERACTIONS_STARS
+      sparts[k].num_ngb_force = 0;
+#endif
       if (sparts[k].time_bin == time_bin_inhibited) continue;
       time_bin_min = min(time_bin_min, sparts[k].time_bin);
       time_bin_max = max(time_bin_max, sparts[k].time_bin);
diff --git a/src/stars/Default/stars.h b/src/stars/Default/stars.h
index 67c665c76691321b9f06b84775278464752f7269..c6b9fe82f67e8f5623b91596055ff74808dbc5cb 100644
--- a/src/stars/Default/stars.h
+++ b/src/stars/Default/stars.h
@@ -76,8 +76,6 @@ __attribute__((always_inline)) INLINE static void stars_reset_predicted_values(
 /**
  * @brief Finishes the calculation of (non-gravity) forces acting on stars
  *
- * Multiplies the forces and accelerations by the appropiate constants
- *
  * @param sp The particle to act upon
  */
 __attribute__((always_inline)) INLINE static void stars_end_force(
@@ -110,6 +108,7 @@ __attribute__((always_inline)) INLINE static void stars_end_density(
   /* Finish the calculation by inserting the missing h-factors */
   sp->density.wcount *= h_inv_dim;
   sp->density.wcount_dh *= h_inv_dim_plus_one;
+
 }
 
 /**