diff --git a/examples/main.c b/examples/main.c
index 10e592109b075069286c74ba27508f5ba73818b1..b97f254acf1268362639c9d8c87431d1895c82b1 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -693,8 +693,8 @@ int main(int argc, char *argv[]) {
 
     /* And initialize the engine with the space and policies. */
     if (myrank == 0) clocks_gettime(&tic);
-    engine_config(1, &e, params, nr_nodes, myrank, nr_threads, with_aff,
-                  talking, restart_file);
+    engine_config(/*restart=*/1, /*fof=*/0, &e, params, nr_nodes, myrank,
+                  nr_threads, with_aff, talking, restart_file);
     if (myrank == 0) {
       clocks_gettime(&toc);
       message("engine_config took %.3f %s.", clocks_diff(&tic, &toc),
@@ -924,8 +924,7 @@ int main(int argc, char *argv[]) {
     if (myrank == 0)
       message(
           "Read %lld gas particles, %lld stars particles, %lld black hole "
-          "particles"
-          " and %lld gparts from the ICs.",
+          "particles and %lld gparts from the ICs.",
           N_total[0], N_total[2], N_total[3], N_total[1]);
 
     /* Verify that the fields to dump actually exist */
@@ -1054,8 +1053,8 @@ int main(int argc, char *argv[]) {
                 &stars_properties, &black_holes_properties,
                 &feedback_properties, &mesh, &potential, &cooling_func,
                 &starform, &chemistry, &fof_properties);
-    engine_config(0, &e, params, nr_nodes, myrank, nr_threads, with_aff,
-                  talking, restart_file);
+    engine_config(/*restart=*/0, /*fof=*/0, &e, params, nr_nodes, myrank,
+                  nr_threads, with_aff, talking, restart_file);
 
     if (myrank == 0) {
       clocks_gettime(&toc);
diff --git a/examples/main_fof.c b/examples/main_fof.c
index e99fa13a6210617288cd5a7f05dbeb0fcb146fab..927b0fc660eb7956560ce1a1c42957c3b5870ab8 100644
--- a/examples/main_fof.c
+++ b/examples/main_fof.c
@@ -290,7 +290,7 @@ int main(int argc, char *argv[]) {
   if (myrank == 0)
     message("WARNING: Debugging checks activated. Code will be slower !");
 #endif
-  
+
   /* Do we choke on FP-exceptions ? */
   if (with_fp_exceptions) {
 #ifdef HAVE_FE_ENABLE_EXCEPT
@@ -339,7 +339,7 @@ int main(int argc, char *argv[]) {
   /* Broadcast the parameter file */
   MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
 #endif
-  
+
   /* Check that we can write the snapshots by testing if the output
    * directory exists and is searchable and writable. */
   char basename[PARSER_MAX_LINE_SIZE];
@@ -391,30 +391,30 @@ int main(int argc, char *argv[]) {
     message("Internal unit system: U_T = %e K.", us.UnitTemperature_in_cgs);
     phys_const_print(&prog_const);
   }
-  
+
   /* Read particles and space information from ICs */
   char ICfileName[200] = "";
   parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
   const int periodic =
-    parser_get_param_int(params, "InitialConditions:periodic");
+      parser_get_param_int(params, "InitialConditions:periodic");
   const int replicate =
-    parser_get_opt_param_int(params, "InitialConditions:replicate", 1);
+      parser_get_opt_param_int(params, "InitialConditions:replicate", 1);
   clean_smoothing_length_values = parser_get_opt_param_int(
-							   params, "InitialConditions:cleanup_smoothing_lengths", 0);
+      params, "InitialConditions:cleanup_smoothing_lengths", 0);
   const int cleanup_h = parser_get_opt_param_int(
-						 params, "InitialConditions:cleanup_h_factors", 0);
+      params, "InitialConditions:cleanup_h_factors", 0);
   const int cleanup_sqrt_a = parser_get_opt_param_int(
-						      params, "InitialConditions:cleanup_velocity_factors", 0);  
+      params, "InitialConditions:cleanup_velocity_factors", 0);
   /* Initialise the cosmology */
   cosmology_init(params, &us, &prog_const, &cosmo);
 
   /* Initialise the gravity scheme */
   gravity_props_init(&gravity_properties, params, &cosmo, /*with_cosmology=*/1,
-		     periodic);
+                     periodic);
 
   /* Initialise the hydro scheme */
   hydro_props_init(&hydro_properties, &prog_const, &us, params);
-  
+
   /* Be verbose about what happens next */
   if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
   if (myrank == 0 && cleanup_h)
@@ -422,7 +422,7 @@ int main(int argc, char *argv[]) {
   if (myrank == 0 && cleanup_sqrt_a)
     message("Cleaning up a-factors from velocity (a=%f)", cosmo.a);
   fflush(stdout);
-  
+
   /* Get ready to read particles of all kinds */
   size_t Ngas = 0, Ngpart = 0, Nspart = 0, Nbpart = 0;
   double dim[3] = {0., 0., 0.};
@@ -431,48 +431,46 @@ int main(int argc, char *argv[]) {
 #if defined(WITH_MPI)
 #if defined(HAVE_PARALLEL_HDF5)
   read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
-		   &Ngas, &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs,
-		   /*with_hydro=*/1, /*with_grav=*/ 1,
-		   /*with_stars=*/1, /*with_black_holes=*/1, cleanup_h, cleanup_sqrt_a,
-		   cosmo.h, cosmo.a, myrank, nr_nodes, MPI_COMM_WORLD,
-		   MPI_INFO_NULL, nr_threads, /*dry_run=*/0);
+                   &Ngas, &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs,
+                   /*with_hydro=*/1, /*with_grav=*/1,
+                   /*with_stars=*/1, /*with_black_holes=*/1, cleanup_h,
+                   cleanup_sqrt_a, cosmo.h, cosmo.a, myrank, nr_nodes,
+                   MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, /*dry_run=*/0);
 #else
-  read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
-		 &Ngas, &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs,
-		 /*with_hydro=*/1, /*with_grav=*/ 1,
-		 /*with_stars=*/1, /*with_black_holes=*/1,
-		 cleanup_h, cleanup_sqrt_a,
-		 cosmo.h, cosmo.a, myrank, nr_nodes, MPI_COMM_WORLD,
-		 MPI_INFO_NULL, nr_threads, /*dry_run=*/0);
+  read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas,
+                 &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs,
+                 /*with_hydro=*/1, /*with_grav=*/1,
+                 /*with_stars=*/1, /*with_black_holes=*/1, cleanup_h,
+                 cleanup_sqrt_a, cosmo.h, cosmo.a, myrank, nr_nodes,
+                 MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, /*dry_run=*/0);
 #endif
 #else
-  read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
-		 &Ngas, &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs,
-		 /*with_hydro=*/1, /*with_grav=*/ 1,
-		 /*with_stars=*/1, /*with_black_holes=*/1,
-		 cleanup_h, cleanup_sqrt_a,
-		 cosmo.h, cosmo.a, nr_threads, /*dry_run=*/0);
+  read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas,
+                 &Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs,
+                 /*with_hydro=*/1, /*with_grav=*/1,
+                 /*with_stars=*/1, /*with_black_holes=*/1, cleanup_h,
+                 cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads, /*dry_run=*/0);
 #endif
 #endif
   if (myrank == 0) {
     clocks_gettime(&toc);
-    message("Reading initial conditions took %.3f %s.",
-	    clocks_diff(&tic, &toc), clocks_getunit());
+    message("Reading initial conditions took %.3f %s.", clocks_diff(&tic, &toc),
+            clocks_getunit());
     fflush(stdout);
   }
-  
+
 #ifdef SWIFT_DEBUG_CHECKS
   /* Check that the other links are correctly set */
-  part_verify_links(parts, gparts, sparts, bparts, Ngas, Ngpart, Nspart,
-		    Nbpart, 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[4] = {0, 0, 0};
 #if defined(WITH_MPI)
   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);
+                MPI_COMM_WORLD);
 #else
   N_total[0] = Ngas;
   N_total[1] = Ngpart;
@@ -482,146 +480,146 @@ int main(int argc, char *argv[]) {
 
   if (myrank == 0)
     message(
-	    "Read %lld gas particles, %lld stars particles, %lld black hole "
-	    "particles and %lld gparts from the ICs.",
-	    N_total[0], N_total[2], N_total[3], N_total[1]);
-  
+        "Read %lld gas particles, %lld stars particles, %lld black hole "
+        "particles and %lld gparts from the ICs.",
+        N_total[0], N_total[2], N_total[3], N_total[1]);
+
   /* Initialize the space with these data. */
   if (myrank == 0) clocks_gettime(&tic);
   space_init(&s, params, &cosmo, dim, parts, gparts, sparts, bparts, Ngas,
-	     Ngpart, Nspart, Nbpart, periodic, replicate, /*generate_gas_in_ics=*/0,
-	     N_total[0] > 0, 1, /*with_star_formation=*/0, talking,
-	     /*dry_run=*/0);
-  
+             Ngpart, Nspart, Nbpart, periodic, replicate,
+             /*generate_gas_in_ics=*/0, N_total[0] > 0, 1,
+             /*with_star_formation=*/0, talking,
+             /*dry_run=*/0);
+
   if (myrank == 0) {
     clocks_gettime(&toc);
     message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
-	    clocks_getunit());
+            clocks_getunit());
     fflush(stdout);
   }
-  
+
   /* Initialise the long-range gravity mesh */
   if (periodic) {
 #ifdef HAVE_FFTW
     pm_mesh_init(&mesh, &gravity_properties, s.dim, nr_threads);
 #else
     /* Need the FFTW library if periodic and self gravity. */
-    error(
-          "No FFTW library found. Cannot compute periodic long-range forces.");
+    error("No FFTW library found. Cannot compute periodic long-range forces.");
 #endif
   } else {
     pm_mesh_init_no_mesh(&mesh, s.dim);
   }
-  
-  /* Also update the total counts (in case of changes due to replication) */
+
+    /* Also update the total counts (in case of changes due to replication) */
 #if defined(WITH_MPI)
   N_long[0] = s.nr_parts;
   N_long[1] = s.nr_gparts;
   N_long[2] = s.nr_sparts;
   MPI_Allreduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM,
-		MPI_COMM_WORLD);
+                MPI_COMM_WORLD);
 #else
   N_total[0] = s.nr_parts;
   N_total[1] = s.nr_gparts;
   N_total[2] = s.nr_sparts;
 #endif
-  
+
   /* Say a few nice things about the space we just created. */
   if (myrank == 0) {
     message("space dimensions are [ %.3f %.3f %.3f ].", s.dim[0], s.dim[1],
-	    s.dim[2]);
+            s.dim[2]);
     message("space %s periodic.", s.periodic ? "is" : "isn't");
     message("highest-level cell dimensions are [ %i %i %i ].", s.cdim[0],
-	    s.cdim[1], s.cdim[2]);
+            s.cdim[1], s.cdim[2]);
     message("%zi parts in %i cells.", s.nr_parts, s.tot_cells);
     message("%zi gparts in %i cells.", s.nr_gparts, s.tot_cells);
     message("%zi sparts in %i cells.", s.nr_sparts, s.tot_cells);
     message("maximum depth is %d.", s.maxdepth);
     fflush(stdout);
   }
-  
+
   /* Verify that each particle is in it's proper cell. */
   if (talking) {
     int icount = 0;
     space_map_cells_pre(&s, 0, &map_cellcheck, &icount);
     message("map_cellcheck picked up %i parts.", icount);
   }
-  
+
   /* Verify the maximal depth of cells. */
   if (talking) {
     int data[2] = {s.maxdepth, 0};
     space_map_cells_pre(&s, 0, &map_maxdepth, data);
     message("nr of cells at depth %i is %i.", data[0], data[1]);
   }
-  
+
   /* Initialise the FOF properties */
   bzero(&fof_properties, sizeof(struct fof_props));
   if (with_fof) fof_init(&fof_properties, params, &prog_const, &us);
-  
+
   /* Construct the engine policy */
   int engine_policies = ENGINE_POLICY | engine_policy_steal;
   engine_policies |= engine_policy_self_gravity;
   engine_policies |= engine_policy_cosmology;
   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], N_total[3],
-	      engine_policies, talking, &reparttype, &us, &prog_const, &cosmo,
-	      &hydro_properties, /*entropy_floor=*/NULL, &gravity_properties,
-	      /*stars_properties=*/NULL, /*black_holes_properties=*/NULL,
-	      /*feedback_properties=*/NULL, &mesh, /*potential=*/NULL, /*cooling_func=*/NULL,
-	      /*starform=*/NULL, /*chemistry=*/NULL, &fof_properties);
-  engine_config(0, &e, params, nr_nodes, myrank, nr_threads, with_aff,
-		talking, NULL);
-  
+              engine_policies, talking, &reparttype, &us, &prog_const, &cosmo,
+              &hydro_properties, /*entropy_floor=*/NULL, &gravity_properties,
+              /*stars_properties=*/NULL, /*black_holes_properties=*/NULL,
+              /*feedback_properties=*/NULL, &mesh, /*potential=*/NULL,
+              /*cooling_func=*/NULL,
+              /*starform=*/NULL, /*chemistry=*/NULL, &fof_properties);
+  engine_config(/*restart=*/0, /*fof=*/1, &e, params, nr_nodes, myrank,
+                nr_threads, with_aff, talking, NULL);
+
   if (myrank == 0) {
     clocks_gettime(&toc);
     message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
-	    clocks_getunit());
+            clocks_getunit());
     fflush(stdout);
   }
-  
+
   /* Get some info to the user. */
   if (myrank == 0) {
     long long N_DM = N_total[1] - N_total[2] - N_total[3] - N_total[0];
     message(
-	    "Running on %lld gas particles, %lld stars particles %lld black "
-	    "hole particles and %lld DM particles (%lld gravity particles)",
-	    N_total[0], N_total[2], N_total[3], N_total[1] > 0 ? N_DM : 0,
-	    N_total[1]);
+        "Running on %lld gas particles, %lld stars particles %lld black "
+        "hole particles and %lld DM particles (%lld gravity particles)",
+        N_total[0], N_total[2], N_total[3], N_total[1] > 0 ? N_DM : 0,
+        N_total[1]);
     message(
-	    "from t=%.3e until t=%.3e with %d ranks, %d threads / rank and %d "
-	    "task queues / rank (dt_min=%.3e, dt_max=%.3e)...",
-	    e.time_begin, e.time_end, nr_nodes, e.nr_threads, e.sched.nr_queues,
-	    e.dt_min, e.dt_max);
+        "from t=%.3e until t=%.3e with %d ranks, %d threads / rank and %d "
+        "task queues / rank (dt_min=%.3e, dt_max=%.3e)...",
+        e.time_begin, e.time_end, nr_nodes, e.nr_threads, e.sched.nr_queues,
+        e.dt_min, e.dt_max);
     fflush(stdout);
   }
-  
+
 #ifdef WITH_MPI
-    /* Split the space. */
+  /* Split the space. */
   engine_split(&e, &initial_partition);
   engine_redistribute(&e);
 #endif
-  
+
 #ifdef SWIFT_DEBUG_TASKS
   e.tic_step = getticks();
 #endif
-  
+
   /* Initialise the particles */
-  engine_init_particles(&e, flag_entropy_ICs, clean_smoothing_length_values,
-			0);
-  
+  engine_init_particles(&e, flag_entropy_ICs, clean_smoothing_length_values, 0);
+
 #ifdef SWIFT_DEBUG_TASKS
   e.toc_step = getticks();
 #endif
-  
+
   /* Perform the FOF search */
   engine_fof(&e, /*dump_results=*/1, /*seed_black_holes=*/0);
 
   /* Write output. */
   engine_dump_snapshot(&e);
-  
+
 #ifdef WITH_MPI
   MPI_Barrier(MPI_COMM_WORLD);
 #endif
@@ -631,13 +629,13 @@ int main(int argc, char *argv[]) {
 #ifdef SWIFT_DEBUG_TASKS
     task_dump_all(&e, 0);
 #endif
-    
+
     /* Generate the task statistics. */
     char dumpfile[40];
     snprintf(dumpfile, 40, "thread_stats-step%d.dat", 0);
     task_dump_stats(dumpfile, &e, /* header = */ 0, /* allranks = */ 1);
   }
-  
+
 #ifdef SWIFT_DEBUG_THREADPOOL
   /* Dump the task data using the given frequency. */
   if (dump_threadpool) {
diff --git a/src/engine.c b/src/engine.c
index bb5f1e032d4db258257d159584cb673b6837b467..8d0087f0f32c078786f6b8de1d0e69e13fb47df7 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -5091,9 +5091,10 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
  * @param verbose Is this #engine talkative ?
  * @param restart_file The name of our restart file.
  */
-void engine_config(int restart, struct engine *e, struct swift_params *params,
-                   int nr_nodes, int nodeID, int nr_threads, int with_aff,
-                   int verbose, const char *restart_file) {
+void engine_config(int restart, int fof, struct engine *e,
+                   struct swift_params *params, int nr_nodes, int nodeID,
+                   int nr_threads, int with_aff, int verbose,
+                   const char *restart_file) {
 
   /* Store the values and initialise global fields. */
   e->nodeID = nodeID;
@@ -5119,9 +5120,15 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
   e->run_fof = 0;
   engine_rank = nodeID;
 
+  if (restart && fof) {
+    error(
+        "Can't configure the engine to be a stand-alone FOF and restarting "
+        "from a check-point at the same time!");
+  }
+
   /* Welcome message */
   if (e->nodeID == 0) message("Running simulation '%s'.", e->run_name);
-  
+
   /* Get the number of queues */
   int nr_queues =
       parser_get_opt_param_int(params, "Scheduler:nr_queues", nr_threads);
@@ -5307,11 +5314,10 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
           hostname(), git_branch(), git_revision(), compiler_name(),
           compiler_version(), e->nr_threads, e->nr_nodes, SPH_IMPLEMENTATION,
           kernel_name,
-	  e->hydro_properties ? e->hydro_properties->target_neighbours : 0.f,
-          e->hydro_properties ? e->hydro_properties->delta_neighbours: 0.f,
+          e->hydro_properties ? e->hydro_properties->target_neighbours : 0.f,
+          e->hydro_properties ? e->hydro_properties->delta_neighbours : 0.f,
           e->hydro_properties ? e->hydro_properties->eta_neighbours : 0.f,
-	  configuration_options(),
-          compilation_cflags());
+          configuration_options(), compilation_cflags());
 
       fprintf(
           e->file_timesteps,
@@ -5346,92 +5352,93 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
   /* Print policy */
   engine_print_policy(e);
 
-  /* Print information about the hydro scheme */
-  if (e->policy & engine_policy_hydro) {
-    if (e->nodeID == 0) hydro_props_print(e->hydro_properties);
-    if (e->nodeID == 0) entropy_floor_print(e->entropy_floor);
-  }
+  if (!fof) {
 
-  /* Print information about the gravity scheme */
-  if (e->policy & engine_policy_self_gravity)
-    if (e->nodeID == 0) gravity_props_print(e->gravity_properties);
+    /* Print information about the hydro scheme */
+    if (e->policy & engine_policy_hydro) {
+      if (e->nodeID == 0) hydro_props_print(e->hydro_properties);
+      if (e->nodeID == 0) entropy_floor_print(e->entropy_floor);
+    }
 
-  if (e->policy & engine_policy_stars)
-    if (e->nodeID == 0) stars_props_print(e->stars_properties);
+    /* Print information about the gravity scheme */
+    if (e->policy & engine_policy_self_gravity)
+      if (e->nodeID == 0) gravity_props_print(e->gravity_properties);
 
-  /* Check we have sensible time bounds */
-  if (e->time_begin >= e->time_end)
-    error(
-        "Final simulation time (t_end = %e) must be larger than the start time "
-        "(t_beg = %e)",
-        e->time_end, e->time_begin);
+    if (e->policy & engine_policy_stars)
+      if (e->nodeID == 0) stars_props_print(e->stars_properties);
 
-  /* Check we don't have inappropriate time labels */
-  if ((e->snapshot_int_time_label_on == 1) && (e->time_end <= 1.f))
-    error("Snapshot integer time labels enabled but end time <= 1");
+    /* Check we have sensible time bounds */
+    if (e->time_begin >= e->time_end)
+      error(
+          "Final simulation time (t_end = %e) must be larger than the start "
+          "time "
+          "(t_beg = %e)",
+          e->time_end, e->time_begin);
 
-  /* Check we have sensible time-step values */
-  if (e->dt_min > e->dt_max)
-    error(
-        "Minimal time-step size (%e) must be smaller than maximal time-step "
-        "size (%e)",
-        e->dt_min, e->dt_max);
+    /* Check we don't have inappropriate time labels */
+    if ((e->snapshot_int_time_label_on == 1) && (e->time_end <= 1.f))
+      error("Snapshot integer time labels enabled but end time <= 1");
 
-  /* Info about time-steps */
-  if (e->nodeID == 0) {
-    message("Absolute minimal timestep size: %e", e->time_base);
+    /* Check we have sensible time-step values */
+    if (e->dt_min > e->dt_max)
+      error(
+          "Minimal time-step size (%e) must be smaller than maximal time-step "
+          "size (%e)",
+          e->dt_min, e->dt_max);
 
-    float dt_min = e->time_end - e->time_begin;
-    while (dt_min > e->dt_min) dt_min /= 2.f;
+    /* Info about time-steps */
+    if (e->nodeID == 0) {
+      message("Absolute minimal timestep size: %e", e->time_base);
 
-    message("Minimal timestep size (on time-line): %e", dt_min);
+      float dt_min = e->time_end - e->time_begin;
+      while (dt_min > e->dt_min) dt_min /= 2.f;
 
-    float dt_max = e->time_end - e->time_begin;
-    while (dt_max > e->dt_max) dt_max /= 2.f;
+      message("Minimal timestep size (on time-line): %e", dt_min);
 
-    message("Maximal timestep size (on time-line): %e", dt_max);
-  }
+      float dt_max = e->time_end - e->time_begin;
+      while (dt_max > e->dt_max) dt_max /= 2.f;
 
-  if (e->dt_min < e->time_base && e->nodeID == 0)
-    error(
-        "Minimal time-step size smaller than the absolute possible minimum "
-        "dt=%e",
-        e->time_base);
+      message("Maximal timestep size (on time-line): %e", dt_max);
+    }
 
-  if (!(e->policy & engine_policy_cosmology))
-    if (e->dt_max > (e->time_end - e->time_begin) && e->nodeID == 0)
-      error("Maximal time-step size larger than the simulation run time t=%e",
-            e->time_end - e->time_begin);
+    if (e->dt_min < e->time_base && e->nodeID == 0)
+      error(
+          "Minimal time-step size smaller than the absolute possible minimum "
+          "dt=%e",
+          e->time_base);
 
-  /* Deal with outputs */
-  if (e->policy & engine_policy_cosmology) {
+    if (!(e->policy & engine_policy_cosmology))
+      if (e->dt_max > (e->time_end - e->time_begin) && e->nodeID == 0)
+        error("Maximal time-step size larger than the simulation run time t=%e",
+              e->time_end - e->time_begin);
 
-    if (e->delta_time_snapshot <= 1.)
-      error("Time between snapshots (%e) must be > 1.", e->delta_time_snapshot);
+    /* Deal with outputs */
+    if (e->policy & engine_policy_cosmology) {
 
-    if (e->delta_time_statistics <= 1.)
-      error("Time between statistics (%e) must be > 1.",
-            e->delta_time_statistics);
+      if (e->delta_time_snapshot <= 1.)
+        error("Time between snapshots (%e) must be > 1.",
+              e->delta_time_snapshot);
 
-    if (e->a_first_snapshot < e->cosmology->a_begin)
-      error(
-          "Scale-factor of first snapshot (%e) must be after the simulation "
-          "start a=%e.",
-          e->a_first_snapshot, e->cosmology->a_begin);
+      if (e->delta_time_statistics <= 1.)
+        error("Time between statistics (%e) must be > 1.",
+              e->delta_time_statistics);
 
-    if (e->a_first_statistics < e->cosmology->a_begin)
-      error(
-          "Scale-factor of first stats output (%e) must be after the "
-          "simulation start a=%e.",
-          e->a_first_statistics, e->cosmology->a_begin);
+      if (e->a_first_snapshot < e->cosmology->a_begin)
+        error(
+            "Scale-factor of first snapshot (%e) must be after the simulation "
+            "start a=%e.",
+            e->a_first_snapshot, e->cosmology->a_begin);
 
-    if (e->policy & engine_policy_structure_finding) {
+      if (e->a_first_statistics < e->cosmology->a_begin)
+        error(
+            "Scale-factor of first stats output (%e) must be after the "
+            "simulation start a=%e.",
+            e->a_first_statistics, e->cosmology->a_begin);
 
-      if (e->delta_time_stf == -1. && !e->snapshot_invoke_stf)
-        error("A value for `StructureFinding:delta_time` must be specified");
+      if (e->policy & engine_policy_structure_finding) {
 
-      if (e->delta_time_stf <= 1. && e->delta_time_stf != -1.)
-        error("Time between STF (%e) must be > 1.", e->delta_time_stf);
+        if (e->delta_time_stf == -1. && !e->snapshot_invoke_stf)
+          error("A value for `StructureFinding:delta_time` must be specified");
 
       if (e->a_first_stf_output < e->cosmology->a_begin)
         error(
@@ -5453,74 +5460,67 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
     }
   } else {
 
-    if (e->delta_time_snapshot <= 0.)
-      error("Time between snapshots (%e) must be positive.",
-            e->delta_time_snapshot);
+      if (e->delta_time_snapshot <= 0.)
+        error("Time between snapshots (%e) must be positive.",
+              e->delta_time_snapshot);
 
-    if (e->delta_time_statistics <= 0.)
-      error("Time between statistics (%e) must be positive.",
-            e->delta_time_statistics);
+      if (e->delta_time_statistics <= 0.)
+        error("Time between statistics (%e) must be positive.",
+              e->delta_time_statistics);
 
-    /* Find the time of the first output */
-    if (e->time_first_snapshot < e->time_begin)
-      error(
-          "Time of first snapshot (%e) must be after the simulation start "
-          "t=%e.",
-          e->time_first_snapshot, e->time_begin);
+      /* Find the time of the first output */
+      if (e->time_first_snapshot < e->time_begin)
+        error(
+            "Time of first snapshot (%e) must be after the simulation start "
+            "t=%e.",
+            e->time_first_snapshot, e->time_begin);
 
-    if (e->time_first_statistics < e->time_begin)
-      error(
-          "Time of first stats output (%e) must be after the simulation start "
-          "t=%e.",
-          e->time_first_statistics, e->time_begin);
+      if (e->time_first_statistics < e->time_begin)
+        error(
+            "Time of first stats output (%e) must be after the simulation "
+            "start "
+            "t=%e.",
+            e->time_first_statistics, e->time_begin);
 
-    if (e->policy & engine_policy_structure_finding) {
+      if (e->policy & engine_policy_structure_finding) {
 
-      if (e->delta_time_stf == -1. && !e->snapshot_invoke_stf)
-        error("A value for `StructureFinding:delta_time` must be specified");
+        if (e->delta_time_stf == -1. && !e->snapshot_invoke_stf)
+          error("A value for `StructureFinding:delta_time` must be specified");
 
-      if (e->delta_time_stf <= 0. && e->delta_time_stf != -1.)
-        error("Time between STF (%e) must be positive.", e->delta_time_stf);
+        if (e->delta_time_stf <= 0. && e->delta_time_stf != -1.)
+          error("Time between STF (%e) must be positive.", e->delta_time_stf);
 
-      if (e->time_first_stf_output < e->time_begin)
-        error("Time of first STF (%e) must be after the simulation start t=%e.",
+        if (e->time_first_stf_output < e->time_begin)
+          error(
+              "Time of first STF (%e) must be after the simulation start t=%e.",
               e->time_first_stf_output, e->time_begin);
-    }
-
-    if (e->policy & engine_policy_structure_finding) {
 
-      if (e->delta_time_fof <= 0.)
-        error("Time between FOF (%e) must be positive.", e->delta_time_fof);
 
-      if (e->time_first_fof_call < e->time_begin)
-        error("Time of first FOF (%e) must be after the simulation start t=%e.",
-              e->time_first_fof_call, e->time_begin);
-    }
   }
 
-  /* Get the total mass */
-  e->total_mass = 0.;
-  for (size_t i = 0; i < e->s->nr_gparts; ++i)
-    e->total_mass += e->s->gparts[i].mass;
+    /* Get the total mass */
+    e->total_mass = 0.;
+    for (size_t i = 0; i < e->s->nr_gparts; ++i)
+      e->total_mass += e->s->gparts[i].mass;
 
 /* Reduce the total mass */
 #ifdef WITH_MPI
-  MPI_Allreduce(MPI_IN_PLACE, &e->total_mass, 1, MPI_DOUBLE, MPI_SUM,
-                MPI_COMM_WORLD);
+    MPI_Allreduce(MPI_IN_PLACE, &e->total_mass, 1, MPI_DOUBLE, MPI_SUM,
+                  MPI_COMM_WORLD);
 #endif
 
 #if defined(WITH_LOGGER)
-  if (e->nodeID == 0)
-    message(
-        "WARNING: There is currently no way of predicting the output "
-        "size, please use it carefully");
+    if (e->nodeID == 0)
+      message(
+          "WARNING: There is currently no way of predicting the output "
+          "size, please use it carefully");
 #endif
 
-  /* Find the time of the first snapshot output */
-  engine_compute_next_snapshot_time(e);
+    /* Find the time of the first snapshot output */
+    engine_compute_next_snapshot_time(e);
 
-  /* Find the time of the first statistics output */
-  engine_compute_next_statistics_time(e);
+    /* Find the time of the first statistics output */
+    engine_compute_next_statistics_time(e);
 
   /* Find the time of the first stf output */
   if (e->policy & engine_policy_structure_finding) {
@@ -5540,36 +5540,38 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
         "activated at runtime (Use --velociraptor).");
   }
 
-  /* Whether restarts are enabled. Yes by default. Can be changed on restart. */
-  e->restart_dump = parser_get_opt_param_int(params, "Restarts:enable", 1);
+    /* Whether restarts are enabled. Yes by default. Can be changed on restart.
+     */
+    e->restart_dump = parser_get_opt_param_int(params, "Restarts:enable", 1);
 
-  /* Whether to save backup copies of the previous restart files. */
-  e->restart_save = parser_get_opt_param_int(params, "Restarts:save", 1);
+    /* Whether to save backup copies of the previous restart files. */
+    e->restart_save = parser_get_opt_param_int(params, "Restarts:save", 1);
 
-  /* Whether restarts should be dumped on exit. Not by default. Can be changed
-   * on restart. */
-  e->restart_onexit = parser_get_opt_param_int(params, "Restarts:onexit", 0);
+    /* Whether restarts should be dumped on exit. Not by default. Can be changed
+     * on restart. */
+    e->restart_onexit = parser_get_opt_param_int(params, "Restarts:onexit", 0);
 
-  /* Hours between restart dumps. Can be changed on restart. */
-  float dhours =
-      parser_get_opt_param_float(params, "Restarts:delta_hours", 6.0);
-  if (e->nodeID == 0) {
-    if (e->restart_dump)
-      message("Restarts will be dumped every %f hours", dhours);
-    else
-      message("WARNING: restarts will not be dumped");
+    /* Hours between restart dumps. Can be changed on restart. */
+    float dhours =
+        parser_get_opt_param_float(params, "Restarts:delta_hours", 6.0);
+    if (e->nodeID == 0) {
+      if (e->restart_dump)
+        message("Restarts will be dumped every %f hours", dhours);
+      else
+        message("WARNING: restarts will not be dumped");
 
-    if (e->verbose && e->restart_onexit)
-      message("Restarts will be dumped after the final step");
-  }
+      if (e->verbose && e->restart_onexit)
+        message("Restarts will be dumped after the final step");
+    }
 
-  /* Internally we use ticks, so convert into a delta ticks. Assumes we can
-   * convert from ticks into milliseconds. */
-  e->restart_dt = clocks_to_ticks(dhours * 60.0 * 60.0 * 1000.0);
+    /* Internally we use ticks, so convert into a delta ticks. Assumes we can
+     * convert from ticks into milliseconds. */
+    e->restart_dt = clocks_to_ticks(dhours * 60.0 * 60.0 * 1000.0);
 
-  /* The first dump will happen no sooner than restart_dt ticks in the
-   * future. */
-  e->restart_next = getticks() + e->restart_dt;
+    /* The first dump will happen no sooner than restart_dt ticks in the
+     * future. */
+    e->restart_next = getticks() + e->restart_dt;
+  }
 
 /* Construct types for MPI communications */
 #ifdef WITH_MPI
@@ -6487,7 +6489,8 @@ void engine_activate_fof_tasks(struct engine *e) {
  * @param dump_results Are we writing group catalogues to output files?
  * @param seed_black_holes Are we seeding black holes?
  */
-void engine_fof(struct engine *e, const int dump_results, const int seed_black_holes) {
+void engine_fof(struct engine *e, const int dump_results,
+                const int seed_black_holes) {
 
   ticks tic = getticks();
 
diff --git a/src/engine.h b/src/engine.h
index fa677312ae351f95a02a9b47644a7d6c86ee4602..998cd02854ffd8d50b5f334b8299ffb78d4c31d2 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -506,9 +506,10 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
                  const struct star_formation *starform,
                  const struct chemistry_global_data *chemistry,
                  struct fof_props *fof_properties);
-void engine_config(int restart, struct engine *e, struct swift_params *params,
-                   int nr_nodes, int nodeID, int nr_threads, int with_aff,
-                   int verbose, const char *restart_file);
+void engine_config(int restart, int fof, struct engine *e,
+                   struct swift_params *params, int nr_nodes, int nodeID,
+                   int nr_threads, int with_aff, int verbose,
+                   const char *restart_file);
 void engine_dump_index(struct engine *e);
 void engine_launch(struct engine *e);
 void engine_prepare(struct engine *e);
@@ -535,7 +536,8 @@ void engine_unpin(void);
 void engine_clean(struct engine *e);
 int engine_estimate_nr_tasks(const struct engine *e);
 void engine_print_task_counts(const struct engine *e);
-void engine_fof(struct engine *e, const int dump_results, const int seed_black_holes);
+void engine_fof(struct engine *e, const int dump_results,
+                const int seed_black_holes);
 
 /* Function prototypes, engine_maketasks.c. */
 void engine_maketasks(struct engine *e);
diff --git a/src/single_io.c b/src/single_io.c
index 805536c19e8615f7ae3e145ea0ce6b1ae956e1b7..afce1e1e6ebb6a8101ccfa109efd8df9e10fbd33 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -459,7 +459,7 @@ void read_ic_single(const char* fileName,
 
   for (int ptype = 0; ptype < swift_type_count; ++ptype)
     N[ptype] = (numParticles[ptype]) + (numParticles_highWord[ptype] << 32);
-  
+
   /* Get the box size if not cubic */
   dim[0] = boxSize[0];
   dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1];
diff --git a/src/space.c b/src/space.c
index 5855d1479c82f15dc4b203338652583065e00b15..0bf349cd9d0bab26917ccaac9be58d43162b47d8 100644
--- a/src/space.c
+++ b/src/space.c
@@ -4017,8 +4017,8 @@ void space_first_init_parts_mapper(void *restrict map_data, int count,
   const struct gravity_props *grav_props = s->e->gravity_properties;
   const int with_gravity = e->policy & engine_policy_self_gravity;
 
-  //const struct chemistry_global_data *chemistry = e->chemistry;
-  //const struct star_formation *star_formation = e->star_formation;
+  // const struct chemistry_global_data *chemistry = e->chemistry;
+  // const struct star_formation *star_formation = e->star_formation;
   const struct cooling_function_data *cool_func = e->cooling_func;
 
   /* Check that the smoothing lengths are non-zero */
@@ -4067,14 +4067,15 @@ void space_first_init_parts_mapper(void *restrict map_data, int count,
 #endif
 
     /* Also initialise the chemistry */
-    //chemistry_first_init_part(phys_const, us, cosmo, chemistry, &p[k], &xp[k]);
+    // chemistry_first_init_part(phys_const, us, cosmo, chemistry, &p[k],
+    // &xp[k]);
 
     /* Also initialise the star formation */
-    //star_formation_first_init_part(phys_const, us, cosmo, star_formation, &p[k],
-    //&xp[k]);
+    // star_formation_first_init_part(phys_const, us, cosmo, star_formation,
+    // &p[k], &xp[k]);
 
     /* And the cooling */
-    //cooling_first_init_part(phys_const, us, cosmo, cool_func, &p[k], &xp[k]);
+    // cooling_first_init_part(phys_const, us, cosmo, cool_func, &p[k], &xp[k]);
 
     /* And the tracers */
     tracers_first_init_xpart(&p[k], &xp[k], us, phys_const, cosmo, hydro_props,
@@ -4175,7 +4176,7 @@ void space_first_init_sparts_mapper(void *restrict map_data, int count,
   const struct space *restrict s = (struct space *)extra_data;
   const struct engine *e = s->e;
 
-  //const struct chemistry_global_data *chemistry = e->chemistry;
+  // const struct chemistry_global_data *chemistry = e->chemistry;
 
 #ifdef SWIFT_DEBUG_CHECKS
   const ptrdiff_t delta = sp - s->sparts;
@@ -4186,7 +4187,7 @@ void space_first_init_sparts_mapper(void *restrict map_data, int count,
   const int with_feedback = (e->policy & engine_policy_feedback);
 
   const struct cosmology *cosmo = e->cosmology;
-  //const struct stars_props *stars_properties = e->stars_properties;
+  // const struct stars_props *stars_properties = e->stars_properties;
   const float a_factor_vel = cosmo->a;
 
   /* Convert velocities to internal units */
@@ -4222,10 +4223,10 @@ void space_first_init_sparts_mapper(void *restrict map_data, int count,
   /* Initialise the rest */
   for (int k = 0; k < count; k++) {
 
-    //stars_first_init_spart(&sp[k], stars_properties);
+  // stars_first_init_spart(&sp[k], stars_properties);
 
-    /* Also initialise the chemistry */
-    //chemistry_first_init_spart(chemistry, &sp[k]);
+  /* Also initialise the chemistry */
+  // chemistry_first_init_spart(chemistry, &sp[k]);
 
 #ifdef SWIFT_DEBUG_CHECKS
     if (sp[k].gpart && sp[k].gpart->id_or_neg_offset != -(k + delta))