diff --git a/examples/main.c b/examples/main.c
index b00ace12e96409ddd5a52aafe308a8d3a4aa9239..8ff44b2dd8d05109d516d27e4f9ff8074ce8ec8e 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -1034,7 +1034,7 @@ int main(int argc, char *argv[]) {
 
     /* Initialize the engine with the space and policies. */
     if (myrank == 0) clocks_gettime(&tic);
-    engine_init(&e, &s, params, N_total[0], N_total[1], N_total[2],
+    engine_init(&e, &s, params, N_total[0], N_total[1], N_total[2], N_total[3],
                 engine_policies, talking, &reparttype, &us, &prog_const, &cosmo,
                 &hydro_properties, &entropy_floor, &gravity_properties,
                 &stars_properties, &black_holes_properties,
diff --git a/examples/main_fof.c b/examples/main_fof.c
index 253e20fad47503c7bc7e8b9f6a0b6d7a021f0b9f..b743d7311045af73a40ba367dc0842e9d7e80301 100644
--- a/examples/main_fof.c
+++ b/examples/main_fof.c
@@ -96,11 +96,14 @@ int main(int argc, char *argv[]) {
   struct gravity_props gravity_properties;
   struct hydro_props hydro_properties;
   struct stars_props stars_properties;
+  struct feedback_props feedback_properties;
   struct entropy_floor_properties entropy_floor;
+  struct black_holes_props black_holes_properties;
   struct part *parts = NULL;
   struct phys_const prog_const;
   struct space s;
   struct spart *sparts = NULL;
+  struct bpart *bparts = NULL;
   struct unit_system us;
 
   int nr_nodes = 1, myrank = 0;
@@ -156,6 +159,7 @@ int main(int argc, char *argv[]) {
   int with_fof = 1;
   int with_star_formation = 0;
   int with_feedback = 0;
+  int with_black_holes = 0;
   int with_fp_exceptions = 0;
   int with_drift_all = 0;
   int with_mpole_reconstruction = 0;
@@ -708,10 +712,24 @@ int main(int argc, char *argv[]) {
     /* Initialise the stars properties */
     if (with_stars)
       stars_props_init(&stars_properties, &prog_const, &us, params,
-                       &hydro_properties);
+                       &hydro_properties, &cosmo);
     else
       bzero(&stars_properties, sizeof(struct stars_props));
 
+    /* Initialise the feedback properties */
+    if (with_feedback) {
+      feedback_props_init(&feedback_properties, &prog_const, &us, params,
+                          &hydro_properties, &cosmo);
+    } else
+      bzero(&feedback_properties, sizeof(struct feedback_props));
+
+    /* Initialise the black holes properties */
+    if (with_black_holes) {
+      black_holes_props_init(&black_holes_properties, &prog_const, &us, params,
+                             &hydro_properties, &cosmo);
+    } else
+      bzero(&black_holes_properties, sizeof(struct black_holes_props));
+
     /* Initialise the gravity properties */
     if (with_self_gravity)
       gravity_props_init(&gravity_properties, params, &cosmo, with_cosmology,
@@ -728,32 +746,32 @@ int main(int argc, char *argv[]) {
     fflush(stdout);
 
     /* Get ready to read particles of all kinds */
-    size_t Ngas = 0, Ngpart = 0, Nspart = 0;
+    size_t Ngas = 0, Ngpart = 0, Nspart = 0, Nbpart = 0;
     double dim[3] = {0., 0., 0.};
     if (myrank == 0) clocks_gettime(&tic);
 #if defined(HAVE_HDF5)
 #if defined(WITH_MPI)
 #if defined(HAVE_PARALLEL_HDF5)
-    read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
-                     &Ngpart, &Nspart, &flag_entropy_ICs, with_hydro,
-                     (with_external_gravity || with_self_gravity), with_stars,
-                     cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank,
-                     nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
-                     dry_run);
+    read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
+                     &Ngas, &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs,
+                     with_hydro, (with_external_gravity || with_self_gravity),
+                     with_stars, with_black_holes, cleanup_h, cleanup_sqrt_a,
+                     cosmo.h, cosmo.a, myrank, nr_nodes, MPI_COMM_WORLD,
+                     MPI_INFO_NULL, nr_threads, dry_run);
 #else
-    read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
-                   &Ngpart, &Nspart, &flag_entropy_ICs, with_hydro,
-                   (with_external_gravity || with_self_gravity), with_stars,
-                   cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank,
-                   nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
-                   dry_run);
+    read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
+                   &Ngas, &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs,
+                   with_hydro, (with_external_gravity || with_self_gravity),
+                   with_stars, with_black_holes, cleanup_h, cleanup_sqrt_a,
+                   cosmo.h, cosmo.a, myrank, nr_nodes, MPI_COMM_WORLD,
+                   MPI_INFO_NULL, nr_threads, dry_run);
 #endif
 #else
-    read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
-                   &Ngpart, &Nspart, &flag_entropy_ICs, with_hydro,
-                   (with_external_gravity || with_self_gravity), with_stars,
-                   cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads,
-                   dry_run);
+    read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
+                   &Ngas, &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs,
+                   with_hydro, (with_external_gravity || with_self_gravity),
+                   with_stars, with_black_holes, cleanup_h, cleanup_sqrt_a,
+                   cosmo.h, cosmo.a, nr_threads, dry_run);
 #endif
 #endif
     if (myrank == 0) {
@@ -776,19 +794,21 @@ int main(int argc, char *argv[]) {
 
     /* Check that the other links are correctly set */
     if (!dry_run)
-      part_verify_links(parts, gparts, sparts, Ngas, Ngpart, Nspart, 1);
+      part_verify_links(parts, gparts, sparts, bparts, Ngas, Ngpart, Nspart,
+                        Nbpart, 1);
 #endif
 
     /* Get the total number of particles across all nodes. */
-    long long N_total[3] = {0, 0, 0};
+    long long N_total[4] = {0, 0, 0};
 #if defined(WITH_MPI)
-    long long N_long[3] = {Ngas, Ngpart, Nspart};
-    MPI_Allreduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM,
+    long long N_long[4] = {Ngas, Ngpart, Nspart, Nbpart};
+    MPI_Allreduce(&N_long, &N_total, 4, MPI_LONG_LONG_INT, MPI_SUM,
                   MPI_COMM_WORLD);
 #else
     N_total[0] = Ngas;
     N_total[1] = Ngpart;
     N_total[2] = Nspart;
+    N_total[3] = Nbpart;
 #endif
 
     if (myrank == 0)
@@ -802,9 +822,10 @@ int main(int argc, char *argv[]) {
 
     /* Initialize the space with these data. */
     if (myrank == 0) clocks_gettime(&tic);
-    space_init(&s, params, &cosmo, dim, parts, gparts, sparts, Ngas, Ngpart,
-               Nspart, periodic, replicate, generate_gas_in_ics, with_hydro,
-               with_self_gravity, with_star_formation, talking, dry_run);
+    space_init(&s, params, &cosmo, dim, parts, gparts, sparts, bparts, Ngas,
+               Ngpart, Nspart, Nbpart, periodic, replicate, generate_gas_in_ics,
+               with_hydro, with_self_gravity, with_star_formation, talking,
+               dry_run);
 
     if (myrank == 0) {
       clocks_gettime(&toc);
@@ -878,7 +899,8 @@ int main(int argc, char *argv[]) {
     if (myrank == 0) potential_print(&potential);
 
     /* Initialise the cooling function properties */
-    if (with_cooling) cooling_init(params, &us, &prog_const, &cooling_func);
+    if (with_cooling)
+      cooling_init(params, &us, &prog_const, &hydro_properties, &cooling_func);
     if (myrank == 0) cooling_print(&cooling_func);
 
     /* Initialise the star formation law and its properties */
@@ -905,17 +927,19 @@ int main(int argc, char *argv[]) {
     if (with_stars) engine_policies |= engine_policy_stars;
     if (with_fof) engine_policies |= engine_policy_fof;
     if (with_feedback) engine_policies |= engine_policy_feedback;
+    if (with_black_holes) engine_policies |= engine_policy_black_holes;
     if (with_structure_finding)
       engine_policies |= engine_policy_structure_finding;
     if (with_fof) engine_policies |= engine_policy_fof;
 
     /* Initialize the engine with the space and policies. */
     if (myrank == 0) clocks_gettime(&tic);
-    engine_init(&e, &s, params, N_total[0], N_total[1], N_total[2],
+    engine_init(&e, &s, params, N_total[0], N_total[1], N_total[2], N_total[3],
                 engine_policies, talking, &reparttype, &us, &prog_const, &cosmo,
                 &hydro_properties, &entropy_floor, &gravity_properties,
-                &stars_properties, &mesh, &potential, &cooling_func, &starform,
-                &chemistry);
+                &stars_properties, &black_holes_properties,
+                &feedback_properties, &mesh, &potential, &cooling_func,
+                &starform, &chemistry);
     engine_config(0, &e, params, nr_nodes, myrank, nr_threads, with_aff,
                   talking, restart_file);
 
diff --git a/src/Makefile.am b/src/Makefile.am
index c479d8e40ab430c02954961ff36dc4c8259f6fba..3cb445ffe7db5afb1f34dac72cfcbd8b2d1bffec 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -78,15 +78,16 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c engine_maketasks.c
     part_type.c xmf.c gravity_properties.c gravity.c \
     collectgroup.c hydro_space.c equation_of_state.c \
     chemistry.c cosmology.c restart.c mesh_gravity.c velociraptor_interface.c \
-    outputlist.c velociraptor_dummy.c logger_io.c memuse.c \
+    outputlist.c velociraptor_dummy.c logger_io.c memuse.c fof.c \
+    c_hashmap/hashmap.c \
     $(EAGLE_COOLING_SOURCES) $(EAGLE_FEEDBACK_SOURCES)
 
 # Include files for distribution, not installation.
 nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
 		 gravity_iact.h kernel_long_gravity.h vector.h cache.h runner_doiact.h runner_doiact_vec.h runner_doiact_grav.h  \
-         runner_doiact_nosort.h runner_doiact_stars.h runner_doiact_black_holes.h units.h intrinsics.h minmax.h \
-         kick.h timestep.h drift.h adiabatic_index.h io_properties.h dimension.h part_type.h periodic.h memswap.h \
-         dump.h logger.h sign.h logger_io.h timestep_limiter.h \
+                 runner_doiact_nosort.h runner_doiact_stars.h runner_doiact_black_holes.h units.h intrinsics.h minmax.h \
+                 kick.h timestep.h drift.h adiabatic_index.h io_properties.h dimension.h part_type.h periodic.h memswap.h \
+                 dump.h logger.h sign.h logger_io.h timestep_limiter.h c_hashmap/hashmap.h \
 		 gravity.h gravity_io.h gravity_cache.h \
 		 gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
 		 gravity/Default/gravity_debug.h gravity/Default/gravity_part.h  \
diff --git a/src/engine.c b/src/engine.c
index 35a7a1f3bf669c09ffa2a1af58bc190ebe1b6b5d..de94817fe564b22222eca6dd6dc74a5942910dc6 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -4887,7 +4887,8 @@ void engine_unpin(void) {
  */
 void engine_init(struct engine *e, struct space *s, struct swift_params *params,
                  long long Ngas, long long Ngparts, long long Nstars,
-                 int policy, int verbose, struct repartition *reparttype,
+                 long long Nblackholes, int policy, int verbose,
+                 struct repartition *reparttype,
                  const struct unit_system *internal_units,
                  const struct phys_const *physical_constants,
                  struct cosmology *cosmo, struct hydro_props *hydro,
@@ -4910,6 +4911,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
   e->total_nr_parts = Ngas;
   e->total_nr_gparts = Ngparts;
   e->total_nr_sparts = Nstars;
+  e->total_nr_bparts = Nblackholes;
   e->proxy_ind = NULL;
   e->nr_proxies = 0;
   e->reparttype = reparttype;
diff --git a/src/engine.h b/src/engine.h
index a26c5f502c9b551322ccf26e9bcb469679955bf2..32217428e1b878dff14c6325ea3eef9ec48f1612 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -475,7 +475,8 @@ void engine_dump_snapshot(struct engine *e);
 void engine_init_output_lists(struct engine *e, struct swift_params *params);
 void engine_init(struct engine *e, struct space *s, struct swift_params *params,
                  long long Ngas, long long Ngparts, long long Nstars,
-                 int policy, int verbose, struct repartition *reparttype,
+                 long long Nblackholes, int policy, int verbose,
+                 struct repartition *reparttype,
                  const struct unit_system *internal_units,
                  const struct phys_const *physical_constants,
                  struct cosmology *cosmo, struct hydro_props *hydro,
diff --git a/src/timers.c b/src/timers.c
index 2d3a5b11f32390be2f7c00431b7c369b000c5de6..8c193f9f2e15746bd551f8a56edafe7c80753334 100644
--- a/src/timers.c
+++ b/src/timers.c
@@ -108,14 +108,6 @@ const char* timers_names[timer_count] = {
     "do_stars_sort",
     "fof_self",
     "fof_pair",
-    "dopair_limiter",
-    "doself_limiter",
-    "dosub_pair_limiter",
-    "dosub_self_limiter",
-    "drift_spart",
-    "do_limiter",
-    "end_hydro_force",
-    "end_grav_force",
 };
 
 /* File to store the timers */
diff --git a/src/timers.h b/src/timers.h
index 0301b38c6850fdc7f12d3571f8a64e9d29aa69d3..abcd1019a708cfae28094875fa094c9886804d26 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -109,14 +109,6 @@ enum {
   timer_do_stars_sort,
   timer_fof_self,
   timer_fof_pair,
-  timer_dopair_limiter,
-  timer_doself_limiter,
-  timer_dosub_pair_limiter,
-  timer_dosub_self_limiter,
-  timer_drift_spart,
-  timer_do_limiter,
-  timer_end_hydro_force,
-  timer_end_grav_force,
   timer_count,
 };