diff --git a/configure.ac b/configure.ac
index c715fb34baeb469e87977da779d82f48f564999e..6f47f86d304bf666fc261bd6c941ec081cc17a6f 100644
--- a/configure.ac
+++ b/configure.ac
@@ -225,16 +225,6 @@ if test "x$enable_debug" = "xyes"; then
    fi
 fi
 
-# Check if stand-alone FoF is on.
-AC_ARG_ENABLE([stand-alone-fof],
-   [AS_HELP_STRING([--enable-stand-alone-fof],
-     [Activate the compilation of the stand-alone friends-of-friends post-processing tool.],
-   )],
-   [enable_standalone_fof="$enableval"],
-   [enable_standalone_fof="no"]
-)
-AM_CONDITIONAL([HAVEFOF],[test $enable_standalone_fof = "yes"])
-
 # Check if task debugging is on.
 AC_ARG_ENABLE([task-debugging],
    [AS_HELP_STRING([--enable-task-debugging],
@@ -1341,6 +1331,7 @@ case "$with_subgrid" in
 	with_subgrid_star_formation=GEAR
 	with_subgrid_feedback=none
 	with_subgrid_black_holes=none
+	enable_fof=no
    ;;
    EAGLE)
 	with_subgrid_cooling=EAGLE
@@ -1351,12 +1342,40 @@ case "$with_subgrid" in
 	with_subgrid_star_formation=EAGLE
 	with_subgrid_feedback=EAGLE
 	with_subgrid_black_holes=EAGLE
+	enable_fof=yes
    ;;
    *)
       AC_MSG_ERROR([Unknown subgrid choice: $with_subgrid])
    ;;
 esac
 
+# Check if FoF is on.
+AC_ARG_ENABLE([fof],
+   [AS_HELP_STRING([--enable-fof],
+     [Activate the friends-of-friends (FoF) code.],
+   )],
+   [enable_fof="$enableval"],
+   [enable_fof="no"]
+)
+if test "$enable_fof" = "yes"; then
+   AC_DEFINE([WITH_FOF], 1, [Enable FoF])
+fi
+
+# Check if stand-alone FoF is on.
+AC_ARG_ENABLE([stand-alone-fof],
+   [AS_HELP_STRING([--enable-stand-alone-fof],
+     [Activate the compilation of the stand-alone friends-of-friends (FoF) post-processing tool.],
+   )],
+   [enable_standalone_fof="$enableval"],
+   [enable_standalone_fof="no"]
+)
+if test "$enable_standalone_fof" = "yes"; then
+   enable_fof="yes + stand-alone tool"
+   AC_DEFINE([WITH_FOF], 1, [Enable FoF])
+   AC_DEFINE([WITH_STAND_ALONE_FOF], 1, [Enable stand-alone FoF])
+fi
+AM_CONDITIONAL([HAVESTANDALONEFOF],[test $enable_standalone_fof = "yes"])
+
 # Gravity scheme.
 AC_ARG_WITH([gravity],
    [AS_HELP_STRING([--with-gravity=<scheme>],
@@ -1990,6 +2009,7 @@ AC_MSG_RESULT([
    Pthread barriers     : $have_pthread_barrier
    VELOCIraptor enabled : $have_velociraptor
    Particle Logger      : $with_logger
+   FoF activated:       : $enable_fof
 
    Hydro scheme       : $with_hydro
    Dimensionality     : $with_dimension
@@ -2013,7 +2033,6 @@ AC_MSG_RESULT([
    Star feedback model  : $with_feedback
    Black holes model    : $with_black_holes
 
-   Stand-alone FoF tool:       : $enable_standalone_fof
    Individual timers           : $enable_timers
    Task debugging              : $enable_task_debugging
    Threadpool debugging        : $enable_threadpool_debugging
diff --git a/examples/Makefile.am b/examples/Makefile.am
index ba92dd4cd4c2b1a4852714edc1e6c4aadd8b6c32..d8161ec28ec4b6ffc14ce0b0f9861c1b8fde2ce2 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -37,14 +37,14 @@ MPI_FLAGS = -DWITH_MPI $(PARMETIS_INCS) $(METIS_INCS)
 bin_PROGRAMS = swift
 
 # Also build the FOF tool?
-if HAVEFOF
+if HAVESTANDALONEFOF
 bin_PROGRAMS += fof
 endif
 
 # Build MPI versions as well?
 if HAVEMPI
 bin_PROGRAMS += swift_mpi
-if HAVEFOF
+if HAVESTANDALONEFOF
 bin_PROGRAMS += fof_mpi
 endif
 endif
diff --git a/examples/main.c b/examples/main.c
index 09deb8fe9cd6f205fcdadf48c25c454b3fe12bb7..b3577b0e2adec7650bb6acd73b128552071f47d8 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -368,6 +368,12 @@ int main(int argc, char *argv[]) {
     return 1;
   }
 
+  if (with_fof) {
+#ifndef WITH_FOF
+    error("Running with FOF but compiled without it!");
+#endif
+  }
+
   if (with_fof && !with_self_gravity) {
     if (myrank == 0)
       printf(
@@ -842,7 +848,9 @@ int main(int argc, char *argv[]) {
 
     /* Initialise the FOF properties */
     bzero(&fof_properties, sizeof(struct fof_props));
+#ifdef WITH_FOF
     if (with_fof) fof_init(&fof_properties, params, &prog_const, &us);
+#endif
 
     /* Be verbose about what happens next */
     if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
diff --git a/src/Makefile.am b/src/Makefile.am
index 7556089b828cb901f0ba2b1e576cf4793e3e1e8e..0f7cef5e3c6b861ebde639d90624a360f1be7047 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -49,7 +49,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     gravity_softened_derivatives.h vector_power.h collectgroup.h hydro_space.h sort_part.h \
     chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h space_getsid.h utilities.h \
     mesh_gravity.h cbrt.h exp10.h velociraptor_interface.h swift_velociraptor_part.h outputlist.h \
-    logger_io.h tracers_io.h tracers.h tracers_struct.h star_formation_io.h fof.h \
+    logger_io.h tracers_io.h tracers.h tracers_struct.h star_formation_io.h fof.h fof_struct.h fof_io.h \
     star_formation_struct.h star_formation.h star_formation_iact.h \
     star_formation_logger.h star_formation_logger_struct.h \
     velociraptor_struct.h velociraptor_io.h random.h memuse.h black_holes.h black_holes_io.h \
diff --git a/src/common_io.c b/src/common_io.c
index ad6faecd0ded12b4280de9ee8d9665c10b7af220..124a122de7cba6b1b6fde36e0ecfea663f2acb31 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -1329,6 +1329,18 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_gpart_i_mapper, temp_i, N, copySize, 0,
                      (void*)&props);
 
+    } else if (props.convert_gpart_i != NULL) {
+
+      /* Prepare some parameters */
+      int* temp_i = (int*)temp;
+      props.start_temp_i = (int*)temp;
+      props.e = e;
+
+      /* Copy the whole thing into a buffer */
+      threadpool_map((struct threadpool*)&e->threadpool,
+                     io_convert_gpart_i_mapper, temp_i, N, copySize, 0,
+                     (void*)&props);
+
     } else if (props.convert_gpart_d != NULL) {
 
       /* Prepare some parameters */
@@ -1388,6 +1400,18 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_spart_i_mapper, temp_i, N, copySize, 0,
                      (void*)&props);
 
+    } else if (props.convert_spart_i != NULL) {
+
+      /* Prepare some parameters */
+      int* temp_i = (int*)temp;
+      props.start_temp_i = (int*)temp;
+      props.e = e;
+
+      /* Copy the whole thing into a buffer */
+      threadpool_map((struct threadpool*)&e->threadpool,
+                     io_convert_spart_i_mapper, temp_i, N, copySize, 0,
+                     (void*)&props);
+
     } else if (props.convert_spart_d != NULL) {
 
       /* Prepare some parameters */
diff --git a/src/engine.c b/src/engine.c
index 2e63356b729a81598295d49b2837b3b7080b7c01..8f14cfd0088cff9d47bee38662a5aeff357dd8e0 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -5575,8 +5575,10 @@ void engine_config(int restart, int fof, struct engine *e,
   stats_create_mpi_type();
   proxy_create_mpi_type();
   task_create_mpi_comms();
+#ifdef WITH_FOF
   fof_create_mpi_types();
-#endif
+#endif /* WITH_FOF */
+#endif /* WITH_MPI */
 
   if (!fof) {
 
@@ -6256,7 +6258,9 @@ void engine_struct_dump(struct engine *e, FILE *stream) {
   feedback_struct_dump(e->feedback_props, stream);
   black_holes_struct_dump(e->black_holes_properties, stream);
   chemistry_struct_dump(e->chemistry, stream);
+#ifdef WITH_FOF
   fof_struct_dump(e->fof_properties, stream);
+#endif
   parser_struct_dump(e->parameter_file, stream);
   if (e->output_list_snapshots)
     output_list_struct_dump(e->output_list_snapshots, stream);
@@ -6374,10 +6378,12 @@ void engine_struct_restore(struct engine *e, FILE *stream) {
   chemistry_struct_restore(chemistry, stream);
   e->chemistry = chemistry;
 
+#ifdef WITH_FOF
   struct fof_props *fof_props =
       (struct fof_props *)malloc(sizeof(struct fof_props));
   fof_struct_restore(fof_props, stream);
   e->fof_properties = fof_props;
+#endif
 
   struct swift_params *parameter_file =
       (struct swift_params *)malloc(sizeof(struct swift_params));
@@ -6493,6 +6499,8 @@ void engine_activate_fof_tasks(struct engine *e) {
 void engine_fof(struct engine *e, const int dump_results,
                 const int seed_black_holes) {
 
+#ifdef WITH_FOF
+
   ticks tic = getticks();
 
   /* Compute number of DM particles */
@@ -6530,4 +6538,7 @@ void engine_fof(struct engine *e, const int dump_results,
   if (engine_rank == 0)
     message("Complete FOF search took: %.3f %s.",
             clocks_from_ticks(getticks() - tic), clocks_getunit());
+#else
+  error("SWIFT was not compiled with FOF enabled!");
+#endif
 }
diff --git a/src/fof.c b/src/fof.c
index 2f374ed1b1c4236b5cb5dcd67e1a2eaef3f6c939..c7c9323a3ed0e5c52cf09c2644c3754e6126a24f 100644
--- a/src/fof.c
+++ b/src/fof.c
@@ -20,6 +20,8 @@
 /* Config parameters. */
 #include "../config.h"
 
+#ifdef WITH_FOF
+
 /* Some standard headers. */
 #include <errno.h>
 #include <libgen.h>
@@ -236,7 +238,7 @@ void fof_allocate(const struct space *s, const long long total_nr_DM_particles,
   /* Set initial group ID of the gparts */
   const size_t group_id_default = props->group_id_default;
   for (size_t i = 0; i < nr_local_gparts; i++) {
-    gparts[i].group_id = group_id_default;
+    gparts[i].fof_data.group_id = group_id_default;
   }
 
   /* Set initial group index and group size */
@@ -968,7 +970,7 @@ void fof_search_pair_cells_foreign(
         /* Check that the links have not already been added to the list. */
         for (int l = 0; l < local_link_count; l++) {
           if ((local_group_links)[l].group_i == root_i &&
-              (local_group_links)[l].group_j == pj->group_id) {
+              (local_group_links)[l].group_j == pj->fof_data.group_id) {
             found = 1;
             break;
           }
@@ -998,8 +1000,9 @@ void fof_search_pair_cells_foreign(
           local_group_links[local_link_count].group_i_size =
               group_size[root_i - node_offset];
 
-          local_group_links[local_link_count].group_j = pj->group_id;
-          local_group_links[local_link_count].group_j_size = pj->group_size;
+          local_group_links[local_link_count].group_j = pj->fof_data.group_id;
+          local_group_links[local_link_count].group_j_size =
+              pj->fof_data.group_size;
 
           local_link_count++;
         }
@@ -1270,9 +1273,9 @@ void fof_calc_group_mass_mapper(void *map_data, int num_elements,
   for (int ind = 0; ind < num_elements; ind++) {
 
     /* Only check groups above the minimum size. */
-    if (gparts[ind].group_id != group_id_default) {
+    if (gparts[ind].fof_data.group_id != group_id_default) {
 
-      hashmap_key_t index = gparts[ind].group_id - group_id_offset;
+      hashmap_key_t index = gparts[ind].fof_data.group_id - group_id_offset;
       hashmap_value_t *data = hashmap_get(&map, index);
 
       /* Update group mass */
@@ -1361,7 +1364,7 @@ void fof_calc_group_mass(struct fof_props *props, const struct space *s,
   for (size_t i = 0; i < nr_gparts; i++) {
 
     /* Check if the particle is in a group above the threshold. */
-    if (gparts[i].group_id != group_id_default) {
+    if (gparts[i].fof_data.group_id != group_id_default) {
 
       const size_t root = fof_find_global(i, group_index, nr_gparts);
 
@@ -1369,7 +1372,7 @@ void fof_calc_group_mass(struct fof_props *props, const struct space *s,
       if (is_local(root, nr_gparts)) {
 
         const size_t index =
-            gparts[i].group_id - group_id_offset - num_groups_prev;
+            gparts[i].fof_data.group_id - group_id_offset - num_groups_prev;
 
         /* Update group mass */
         group_mass[index] += gparts[i].mass;
@@ -1416,7 +1419,7 @@ void fof_calc_group_mass(struct fof_props *props, const struct space *s,
   for (size_t i = 0; i < nr_gparts; i++) {
 
     /* Only check groups above the minimum size and mass threshold. */
-    if (gparts[i].group_id != group_id_default) {
+    if (gparts[i].fof_data.group_id != group_id_default) {
 
       size_t root = fof_find_global(i, group_index, nr_gparts);
 
@@ -1424,7 +1427,7 @@ void fof_calc_group_mass(struct fof_props *props, const struct space *s,
       if (is_local(root, nr_gparts)) {
 
         const size_t index =
-            gparts[i].group_id - group_id_offset - num_groups_prev;
+            gparts[i].fof_data.group_id - group_id_offset - num_groups_prev;
 
         /* Only seed groups above the mass threshold. */
         if (group_mass[index] > seed_halo_mass) {
@@ -1516,7 +1519,8 @@ void fof_calc_group_mass(struct fof_props *props, const struct space *s,
         (fof_mass_recv[i].global_root >= node_offset + nr_gparts)) {
       error("Received global root index out of range!");
     }
-    group_mass[gparts[fof_mass_recv[i].global_root - node_offset].group_id -
+    group_mass[gparts[fof_mass_recv[i].global_root - node_offset]
+                   .fof_data.group_id -
                group_id_offset - num_groups_prev] +=
         fof_mass_recv[i].group_mass;
   }
@@ -1526,7 +1530,7 @@ void fof_calc_group_mass(struct fof_props *props, const struct space *s,
   for (size_t i = 0; i < nrecv; i++) {
 
     const int offset =
-        gparts[fof_mass_recv[i].global_root - node_offset].group_id -
+        gparts[fof_mass_recv[i].global_root - node_offset].fof_data.group_id -
         group_id_offset - num_groups_prev;
 
     /* Only seed groups above the mass threshold. */
@@ -1563,7 +1567,7 @@ void fof_calc_group_mass(struct fof_props *props, const struct space *s,
     }
 
     const int offset =
-        gparts[fof_mass_recv[i].global_root - node_offset].group_id -
+        gparts[fof_mass_recv[i].global_root - node_offset].fof_data.group_id -
         group_id_offset - num_groups_prev;
 
     /* If the densest particle found locally is not the global max, make sure we
@@ -1663,10 +1667,10 @@ void fof_calc_group_mass(struct fof_props *props, const struct space *s,
   /* JSW TODO: Parallelise with threadpool*/
   for (size_t i = 0; i < nr_gparts; i++) {
 
-    const size_t index = gparts[i].group_id - group_id_offset;
+    const size_t index = gparts[i].fof_data.group_id - group_id_offset;
 
     /* Only check groups above the minimum mass threshold. */
-    if (gparts[i].group_id != group_id_default) {
+    if (gparts[i].fof_data.group_id != group_id_default) {
 
       if (group_mass[index] > seed_halo_mass) {
 
@@ -1952,12 +1956,12 @@ void fof_dump_group_data(const struct fof_props *props,
                                   : -1;
 #ifdef WITH_MPI
     fprintf(file, "  %8zu %12zu %12e %12e %18lld %18lld\n",
-            gparts[group_offset - node_offset].group_id,
+            gparts[group_offset - node_offset].fof_data.group_id,
             group_size[group_offset - node_offset], group_mass[i],
             max_part_density[i], max_part_density_index[i], part_id);
 #else
     fprintf(file, "  %8zu %12zu %12e %12e %18lld %18lld\n",
-            gparts[group_offset].group_id, group_size[group_offset],
+            gparts[group_offset].fof_data.group_id, group_size[group_offset],
             group_mass[i], max_part_density[i], max_part_density_index[i],
             part_id);
 #endif
@@ -2096,8 +2100,8 @@ void fof_search_foreign_cells(struct fof_props *props, const struct space *s) {
       for (int k = 0; k < local_cell->grav.count; k++) {
         const size_t root =
             fof_find_global(offset[k] - node_offset, group_index, nr_gparts);
-        gparts[k].group_id = root;
-        gparts[k].group_size = group_size[root - node_offset];
+        gparts[k].fof_data.group_id = root;
+        gparts[k].fof_data.group_size = group_size[root - node_offset];
       }
     }
   }
@@ -2539,7 +2543,8 @@ void fof_search_tree(struct fof_props *props,
         cmp_func_group_size);
 
   /* Set default group ID for all particles */
-  for (size_t i = 0; i < nr_gparts; i++) gparts[i].group_id = group_id_default;
+  for (size_t i = 0; i < nr_gparts; i++)
+    gparts[i].fof_data.group_id = group_id_default;
 
   /*
     Assign final group IDs to local root particles where the global root is on
@@ -2548,10 +2553,10 @@ void fof_search_tree(struct fof_props *props,
   */
   for (size_t i = 0; i < num_groups_local; i++) {
 #ifdef WITH_MPI
-    gparts[high_group_sizes[i].index - node_offset].group_id =
+    gparts[high_group_sizes[i].index - node_offset].fof_data.group_id =
         group_id_offset + i + num_groups_prev;
 #else
-    gparts[high_group_sizes[i].index].group_id = group_id_offset + i;
+    gparts[high_group_sizes[i].index].fof_data.group_id = group_id_offset + i;
 #endif
   }
 
@@ -2644,7 +2649,7 @@ void fof_search_tree(struct fof_props *props,
       error("Received global root index out of range!");
     }
     fof_index_recv[i].global_root =
-        gparts[fof_index_recv[i].global_root - node_offset].group_id;
+        gparts[fof_index_recv[i].global_root - node_offset].fof_data.group_id;
   }
 
   /* Send the result back */
@@ -2658,7 +2663,7 @@ void fof_search_tree(struct fof_props *props,
         (fof_index_send[i].local_root >= node_offset + nr_gparts)) {
       error("Sent local root index out of range!");
     }
-    gparts[fof_index_send[i].local_root - node_offset].group_id =
+    gparts[fof_index_send[i].local_root - node_offset].fof_data.group_id =
         fof_index_send[i].global_root;
   }
 
@@ -2674,7 +2679,7 @@ void fof_search_tree(struct fof_props *props,
   /* Assign every particle the group_id of its local root. */
   for (size_t i = 0; i < nr_gparts; i++) {
     const size_t root = fof_find_local(i, nr_gparts, group_index);
-    gparts[i].group_id = gparts[root].group_id;
+    gparts[i].fof_data.group_id = gparts[root].fof_data.group_id;
   }
 
   if (verbose)
@@ -2769,3 +2774,5 @@ void fof_struct_restore(struct fof_props *props, FILE *stream) {
   restart_read_blocks((void *)props, sizeof(struct fof_props), 1, stream, NULL,
                       "fof_props");
 }
+
+#endif /* WITH_FOF */
diff --git a/src/fof_io.h b/src/fof_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..e6e91c577ad2676cc9bb2d731b2250c9c05cafd4
--- /dev/null
+++ b/src/fof_io.h
@@ -0,0 +1,138 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_FOF_IO_H
+#define SWIFT_FOF_IO_H
+
+/* Config parameters. */
+#include "../config.h"
+
+#ifdef WITH_FOF
+
+INLINE static void convert_part_group_id(const struct engine* e,
+                                         const struct part* p,
+                                         const struct xpart* xp,
+                                         long long* ret) {
+  ret[0] = p->gpart->fof_data.group_id;
+}
+
+INLINE static void convert_spart_group_id(const struct engine* e,
+                                          const struct spart* sp,
+                                          long long* ret) {
+  ret[0] = sp->gpart->fof_data.group_id;
+}
+
+INLINE static void convert_bpart_group_id(const struct engine* e,
+                                          const struct bpart* bp,
+                                          long long* ret) {
+  ret[0] = bp->gpart->fof_data.group_id;
+}
+
+#endif /* WITH_FOF */
+
+/**
+ * @brief Specifies which FoF-related particle fields to write to a dataset
+ *
+ * @param parts The particle array.
+ * @param xparts The extended particle array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return The number of fields to write.
+ */
+INLINE static int fof_write_parts(const struct part* parts,
+                                  const struct xpart* xparts,
+                                  struct io_props* list) {
+
+#ifdef WITH_FOF
+
+  list[0] = io_make_output_field_convert_part("GroupIDs", LONGLONG, 1,
+                                              UNIT_CONV_NO_UNITS, parts, xparts,
+                                              convert_part_group_id);
+  return 1;
+#else
+  return 0;
+#endif
+}
+
+/**
+ * @brief Specifies which FoF-related g-particle fields to write to a dataset
+ *
+ * @param gparts The g-particle array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return The number of fields to write.
+ */
+INLINE static int fof_write_gparts(const struct gpart* gparts,
+                                   struct io_props* list) {
+
+#ifdef WITH_FOF
+
+  list[0] = io_make_output_field("GroupIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS,
+                                 gparts, fof_data.group_id);
+
+  return 1;
+#else
+  return 0;
+#endif
+}
+
+/**
+ * @brief Specifies which FoF-related s-particle fields to write to a dataset
+ *
+ * @param sparts The s-particle array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return The number of fields to write.
+ */
+INLINE static int fof_write_sparts(const struct spart* sparts,
+                                   struct io_props* list) {
+
+#ifdef WITH_FOF
+
+  list[0] = io_make_output_field_convert_spart("GroupIDs", LONGLONG, 1,
+                                               UNIT_CONV_NO_UNITS, sparts,
+                                               convert_spart_group_id);
+  return 1;
+#else
+  return 0;
+#endif
+}
+
+/**
+ * @brief Specifies which FoF-related b-particle fields to write to a dataset
+ *
+ * @param bparts The b-particle array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return The number of fields to write.
+ */
+INLINE static int fof_write_bparts(const struct bpart* bparts,
+                                   struct io_props* list) {
+
+#ifdef WITH_FOF
+
+  list[0] = io_make_output_field_convert_bpart("GroupIDs", LONGLONG, 1,
+                                               UNIT_CONV_NO_UNITS, bparts,
+                                               convert_bpart_group_id);
+  return 1;
+#else
+  return 0;
+#endif
+}
+
+#endif /* SWIFT_FOF_IO_H */
diff --git a/src/fof_struct.h b/src/fof_struct.h
new file mode 100644
index 0000000000000000000000000000000000000000..38022b56db7d273291e2fa783c48adb7887f516f
--- /dev/null
+++ b/src/fof_struct.h
@@ -0,0 +1,48 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_FOF_STRUCT_H
+#define SWIFT_FOF_STRUCT_H
+
+/* Config parameters. */
+#include "../config.h"
+
+#ifdef WITH_FOF
+
+/**
+ * @brief Particle-carried fields for the FoF scheme.
+ */
+struct fof_gpart_data {
+
+  /*! Particle group ID */
+  size_t group_id;
+
+  /*! Size of the FOF group of this particle */
+  size_t group_size;
+};
+
+#else
+
+/**
+ * @brief Particle-carried fields for the FoF scheme.
+ */
+struct fof_gpart_data {};
+
+#endif
+
+#endif /* SWIFT_FOF_STRUCT_H */
diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h
index 2e443e8fdc2479f3f2feff30281dccad9a7b6519..1ba3899e7ecc346227c10679bb8b704937c625b2 100644
--- a/src/gravity/Default/gravity_io.h
+++ b/src/gravity/Default/gravity_io.h
@@ -104,7 +104,7 @@ INLINE static void darkmatter_write_particles(const struct gpart* gparts,
                                               int* num_fields) {
 
   /* Say how much we want to write */
-  *num_fields = 5;
+  *num_fields = 4;
 
   /* List what we want to write */
   list[0] = io_make_output_field_convert_gpart(
@@ -115,8 +115,6 @@ INLINE static void darkmatter_write_particles(const struct gpart* gparts,
       io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, gparts, mass);
   list[3] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
                                  UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset);
-  list[4] = io_make_output_field("GroupIDs", INT, 1, UNIT_CONV_NO_UNITS, gparts,
-                                 group_id);
 }
 
 #endif /* SWIFT_DEFAULT_GRAVITY_IO_H */
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index 2c16a34a11d3a7674f5a6c122da780aeaff4a9dc..29f6bb6a64202c3fc3407407932e8f2e31f994af 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -19,7 +19,11 @@
 #ifndef SWIFT_DEFAULT_GRAVITY_PART_H
 #define SWIFT_DEFAULT_GRAVITY_PART_H
 
-/* Gravity particle. */
+#include "fof_struct.h"
+
+/**
+ * @brief Gravity particle.
+ */
 struct gpart {
 
   /*! Particle ID. If negative, it is the negative offset of the #part with
@@ -38,15 +42,15 @@ struct gpart {
   /*! Particle mass. */
   float mass;
 
+  /*! Particle FoF properties (group ID, group size, ...) */
+  struct fof_gpart_data fof_data;
+
   /*! Time-step length */
   timebin_t time_bin;
 
   /*! Type of the #gpart (DM, gas, star, ...) */
   enum part_type type;
 
-  /* Particle group ID and size in the FOF. */
-  size_t group_id, group_size;
-
 #ifdef SWIFT_DEBUG_CHECKS
 
   /* Numer of gparts this gpart interacted with */
diff --git a/src/gravity/Potential/gravity_part.h b/src/gravity/Potential/gravity_part.h
index 7b47a6198333edf7e6acbd0daf37c6b9c1eb6ce1..97126731f0dda7484600b05c0fd04faf72165d63 100644
--- a/src/gravity/Potential/gravity_part.h
+++ b/src/gravity/Potential/gravity_part.h
@@ -19,7 +19,11 @@
 #ifndef SWIFT_POTENTIAL_GRAVITY_PART_H
 #define SWIFT_POTENTIAL_GRAVITY_PART_H
 
-/* Gravity particle. */
+#include "fof_struct.h"
+
+/**
+ * @brief Gravity particle.
+ */
 struct gpart {
 
   /*! Particle ID. If negative, it is the negative offset of the #part with
@@ -41,15 +45,15 @@ struct gpart {
   /*! Gravitational potential */
   float potential;
 
+  /*! Particle FoF properties (group ID, group size, ...) */
+  struct fof_gpart_data fof_data;
+
   /*! Time-step length */
   timebin_t time_bin;
 
   /*! Type of the #gpart (DM, gas, star, ...) */
   enum part_type type;
 
-  /* Particle group ID and size in the FOF. */
-  size_t group_id, group_size;
-
 #ifdef SWIFT_DEBUG_CHECKS
 
   /* Numer of gparts this gpart interacted with */
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 5d68dc5d9f0f4ac0d12f7829f2c2a12d9964fe4a..54154ce970a9fccf0f2f6324d86faaebed392555 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -130,16 +130,6 @@ INLINE static void convert_part_potential(const struct engine* e,
     ret[0] = 0.f;
 }
 
-INLINE static void convert_part_group_id(const struct engine* e,
-                                         const struct part* p,
-                                         const struct xpart* xp, int* ret) {
-
-  if (p->gpart != NULL)
-    ret[0] = p->gpart->group_id;
-  else
-    ret[0] = 0;
-}
-
 /**
  * @brief Specifies which particle fields to write to a dataset
  *
@@ -153,7 +143,7 @@ INLINE static void hydro_write_particles(const struct part* parts,
                                          struct io_props* list,
                                          int* num_fields) {
 
-  *num_fields = 11;
+  *num_fields = 10;
 
   /* List what we want to write */
   list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
@@ -180,9 +170,6 @@ INLINE static void hydro_write_particles(const struct part* parts,
   list[9] = io_make_output_field_convert_part("Potential", FLOAT, 1,
                                               UNIT_CONV_POTENTIAL, parts,
                                               xparts, convert_part_potential);
-  list[10] =
-      io_make_output_field_convert_part("GroupIDs", INT, 1, UNIT_CONV_NO_UNITS,
-                                        parts, xparts, convert_part_group_id);
 
 #ifdef DEBUG_INTERACTIONS_SPH
 
diff --git a/src/parallel_io.c b/src/parallel_io.c
index b90a3848d9acc88a8f41956a27f40c72277b7524..3b43e623784fd60735ac759e89fcf5536094076f 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -45,6 +45,7 @@
 #include "engine.h"
 #include "entropy_floor.h"
 #include "error.h"
+#include "fof_io.h"
 #include "gravity_io.h"
 #include "gravity_properties.h"
 #include "hydro_io.h"
@@ -1206,6 +1207,7 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
                                               with_cosmology);
         num_fields +=
             star_formation_write_particles(parts, xparts, list + num_fields);
+        num_fields += fof_write_parts(parts, xparts, list + num_fields);
         if (with_stf) {
           num_fields +=
               velociraptor_write_parts(parts, xparts, list + num_fields);
@@ -1214,6 +1216,7 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
 
       case swift_type_dark_matter:
         darkmatter_write_particles(gparts, list, &num_fields);
+        num_fields += fof_write_gparts(gparts, list + num_fields);
         if (with_stf) {
           num_fields += velociraptor_write_gparts(e->s->gpart_group_data,
                                                   list + num_fields);
@@ -1225,6 +1228,7 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
         num_fields += chemistry_write_sparticles(sparts, list + num_fields);
         num_fields +=
             tracers_write_sparticles(sparts, list + num_fields, with_cosmology);
+        num_fields += fof_write_sparts(sparts, list + num_fields);
         if (with_stf) {
           num_fields += velociraptor_write_sparts(sparts, list + num_fields);
         }
@@ -1233,6 +1237,7 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
       case swift_type_black_hole:
         black_holes_write_particles(bparts, list, &num_fields);
         num_fields += chemistry_write_bparticles(bparts, list + num_fields);
+        num_fields += fof_write_bparts(bparts, list + num_fields);
         if (with_stf) {
           num_fields += velociraptor_write_bparts(bparts, list + num_fields);
         }
@@ -1535,6 +1540,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
             num_fields += cooling_write_particles(
                 parts, xparts, list + num_fields, e->cooling_func);
           }
+          num_fields += fof_write_parts(parts, xparts, list + num_fields);
           if (with_stf) {
             num_fields +=
                 velociraptor_write_parts(parts, xparts, list + num_fields);
@@ -1573,6 +1579,8 @@ void write_output_parallel(struct engine* e, const char* baseName,
                 cooling_write_particles(parts_written, xparts_written,
                                         list + num_fields, e->cooling_func);
           }
+          num_fields +=
+              fof_write_parts(parts_written, xparts_written, list + num_fields);
           if (with_stf) {
             num_fields += velociraptor_write_parts(
                 parts_written, xparts_written, list + num_fields);
@@ -1590,6 +1598,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
           /* This is a DM-only run without inhibited particles */
           Nparticles = Ntot;
           darkmatter_write_particles(gparts, list, &num_fields);
+          num_fields += fof_write_gparts(gparts, list + num_fields);
           if (with_stf) {
             num_fields += velociraptor_write_gparts(e->s->gpart_group_data,
                                                     list + num_fields);
@@ -1622,11 +1631,10 @@ void write_output_parallel(struct engine* e, const char* baseName,
 
           /* Select the fields to write */
           darkmatter_write_particles(gparts_written, list, &num_fields);
+          num_fields += fof_write_gparts(gparts_written, list + num_fields);
           if (with_stf) {
-#ifdef HAVE_VELOCIRAPTOR
             num_fields += velociraptor_write_gparts(gpart_group_data_written,
                                                     list + num_fields);
-#endif
           }
         }
       } break;
@@ -1640,6 +1648,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
           num_fields += chemistry_write_sparticles(sparts, list + num_fields);
           num_fields += tracers_write_sparticles(sparts, list + num_fields,
                                                  with_cosmology);
+          num_fields += fof_write_sparts(sparts, list + num_fields);
           if (with_stf) {
             num_fields += velociraptor_write_sparts(sparts, list + num_fields);
           }
@@ -1663,6 +1672,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
           num_fields += chemistry_write_sparticles(sparts, list + num_fields);
           num_fields += tracers_write_sparticles(sparts, list + num_fields,
                                                  with_cosmology);
+          num_fields += fof_write_sparts(sparts_written, list + num_fields);
           if (with_stf) {
             num_fields +=
                 velociraptor_write_sparts(sparts_written, list + num_fields);
@@ -1677,7 +1687,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
           Nparticles = Nblackholes;
           black_holes_write_particles(bparts, list, &num_fields);
           num_fields += chemistry_write_bparticles(bparts, list + num_fields);
-
+          num_fields += fof_write_bparts(bparts, list + num_fields);
           if (with_stf) {
             num_fields += velociraptor_write_bparts(bparts, list + num_fields);
           }
@@ -1699,7 +1709,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
           /* Select the fields to write */
           black_holes_write_particles(bparts_written, list, &num_fields);
           num_fields += chemistry_write_bparticles(bparts, list + num_fields);
-
+          num_fields += fof_write_bparts(bparts_written, list + num_fields);
           if (with_stf) {
             num_fields +=
                 velociraptor_write_bparts(bparts_written, list + num_fields);
diff --git a/src/runner.c b/src/runner.c
index b7422b524b6a68a37a60eea580fb35edb30665c7..7f645665eb786ab578fd8e0cfccaa72420c08b1b 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -4973,6 +4973,8 @@ void runner_do_logger(struct runner *r, struct cell *c, int timer) {
  */
 void runner_do_fof_self(struct runner *r, struct cell *c, int timer) {
 
+#ifdef WITH_FOF
+
   TIMER_TIC;
 
   const struct engine *e = r->e;
@@ -4985,6 +4987,10 @@ void runner_do_fof_self(struct runner *r, struct cell *c, int timer) {
   rec_fof_search_self(e->fof_properties, dim, search_r2, periodic, gparts, c);
 
   if (timer) TIMER_TOC(timer_fof_self);
+
+#else
+  error("SWIFT was not compiled with FOF enabled!");
+#endif
 }
 
 /**
@@ -4998,6 +5004,8 @@ void runner_do_fof_self(struct runner *r, struct cell *c, int timer) {
 void runner_do_fof_pair(struct runner *r, struct cell *ci, struct cell *cj,
                         int timer) {
 
+#ifdef WITH_FOF
+
   TIMER_TIC;
 
   const struct engine *e = r->e;
@@ -5011,4 +5019,7 @@ void runner_do_fof_pair(struct runner *r, struct cell *ci, struct cell *cj,
                       cj);
 
   if (timer) TIMER_TOC(timer_fof_pair);
+#else
+  error("SWIFT was not compiled with FOF enabled!");
+#endif
 }
diff --git a/src/serial_io.c b/src/serial_io.c
index 220f7b3eb403a991c8cdc8f3a4ebf6add8f4d261..7a6b73935d7d500594ec920c23594aa538029253 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -45,6 +45,7 @@
 #include "engine.h"
 #include "entropy_floor.h"
 #include "error.h"
+#include "fof_io.h"
 #include "gravity_io.h"
 #include "gravity_properties.h"
 #include "hydro_io.h"
@@ -1172,6 +1173,7 @@ void write_output_serial(struct engine* e, const char* baseName,
                 num_fields += cooling_write_particles(
                     parts, xparts, list + num_fields, e->cooling_func);
               }
+              num_fields += fof_write_parts(parts, xparts, list + num_fields);
               if (with_stf) {
                 num_fields +=
                     velociraptor_write_parts(parts, xparts, list + num_fields);
@@ -1210,6 +1212,8 @@ void write_output_serial(struct engine* e, const char* baseName,
                     cooling_write_particles(parts_written, xparts_written,
                                             list + num_fields, e->cooling_func);
               }
+              num_fields += fof_write_parts(parts_written, xparts_written,
+                                            list + num_fields);
               if (with_stf) {
                 num_fields += velociraptor_write_parts(
                     parts_written, xparts_written, list + num_fields);
@@ -1228,6 +1232,7 @@ void write_output_serial(struct engine* e, const char* baseName,
               /* This is a DM-only run without inhibited particles */
               Nparticles = Ntot;
               darkmatter_write_particles(gparts, list, &num_fields);
+              num_fields += fof_write_gparts(gparts_written, list + num_fields);
               if (with_stf) {
                 num_fields += velociraptor_write_gparts(e->s->gpart_group_data,
                                                         list + num_fields);
@@ -1261,6 +1266,7 @@ void write_output_serial(struct engine* e, const char* baseName,
 
               /* Select the fields to write */
               darkmatter_write_particles(gparts_written, list, &num_fields);
+              num_fields += fof_write_gparts(gparts_written, list + num_fields);
               if (with_stf) {
                 num_fields += velociraptor_write_gparts(
                     gpart_group_data_written, list + num_fields);
@@ -1278,6 +1284,7 @@ void write_output_serial(struct engine* e, const char* baseName,
                   chemistry_write_sparticles(sparts, list + num_fields);
               num_fields += tracers_write_sparticles(sparts, list + num_fields,
                                                      with_cosmology);
+              num_fields += fof_write_sparts(sparts, list + num_fields);
               if (with_stf) {
                 num_fields +=
                     velociraptor_write_sparts(sparts, list + num_fields);
@@ -1300,9 +1307,10 @@ void write_output_serial(struct engine* e, const char* baseName,
               /* Select the fields to write */
               stars_write_particles(sparts_written, list, &num_fields);
               num_fields +=
-                  chemistry_write_sparticles(sparts, list + num_fields);
-              num_fields += tracers_write_sparticles(sparts, list + num_fields,
-                                                     with_cosmology);
+                  chemistry_write_sparticles(sparts_written, list + num_fields);
+              num_fields += tracers_write_sparticles(
+                  sparts_written, list + num_fields, with_cosmology);
+              num_fields += fof_write_sparts(sparts_written, list + num_fields);
               if (with_stf) {
                 num_fields += velociraptor_write_sparts(sparts_written,
                                                         list + num_fields);
@@ -1318,7 +1326,7 @@ void write_output_serial(struct engine* e, const char* baseName,
               black_holes_write_particles(bparts, list, &num_fields);
               num_fields +=
                   chemistry_write_bparticles(bparts, list + num_fields);
-
+              num_fields += fof_write_bparts(bparts, list + num_fields);
               if (with_stf) {
                 num_fields +=
                     velociraptor_write_bparts(bparts, list + num_fields);
@@ -1342,7 +1350,7 @@ void write_output_serial(struct engine* e, const char* baseName,
               black_holes_write_particles(bparts_written, list, &num_fields);
               num_fields +=
                   chemistry_write_bparticles(bparts, list + num_fields);
-
+	      num_fields += fof_write_bparts(bparts_written, list + num_fields);
               if (with_stf) {
                 num_fields += velociraptor_write_bparts(bparts_written,
                                                         list + num_fields);
diff --git a/src/single_io.c b/src/single_io.c
index 0a5d20236d85b7a707054e856f86323172d2ff8c..20d7f98e1bbbe6761b15378508f50104dd14f92d 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -44,6 +44,7 @@
 #include "engine.h"
 #include "entropy_floor.h"
 #include "error.h"
+#include "fof_io.h"
 #include "gravity_io.h"
 #include "gravity_properties.h"
 #include "hydro_io.h"
@@ -971,6 +972,7 @@ void write_output_single(struct engine* e, const char* baseName,
             num_fields += cooling_write_particles(
                 parts, xparts, list + num_fields, e->cooling_func);
           }
+          num_fields += fof_write_parts(parts, xparts, list + num_fields);
           if (with_stf) {
             num_fields +=
                 velociraptor_write_parts(parts, xparts, list + num_fields);
@@ -1009,6 +1011,8 @@ void write_output_single(struct engine* e, const char* baseName,
                 cooling_write_particles(parts_written, xparts_written,
                                         list + num_fields, e->cooling_func);
           }
+          num_fields +=
+              fof_write_parts(parts_written, xparts_written, list + num_fields);
           if (with_stf) {
             num_fields += velociraptor_write_parts(
                 parts_written, xparts_written, list + num_fields);
@@ -1026,6 +1030,7 @@ void write_output_single(struct engine* e, const char* baseName,
           /* This is a DM-only run without inhibited particles */
           N = Ntot;
           darkmatter_write_particles(gparts, list, &num_fields);
+          num_fields += fof_write_gparts(gparts, list + num_fields);
           if (with_stf) {
             num_fields += velociraptor_write_gparts(e->s->gpart_group_data,
                                                     list + num_fields);
@@ -1058,6 +1063,7 @@ void write_output_single(struct engine* e, const char* baseName,
 
           /* Select the fields to write */
           darkmatter_write_particles(gparts_written, list, &num_fields);
+          num_fields += fof_write_gparts(gparts_written, list + num_fields);
           if (with_stf) {
             num_fields += velociraptor_write_gparts(gpart_group_data_written,
                                                     list + num_fields);
@@ -1074,6 +1080,7 @@ void write_output_single(struct engine* e, const char* baseName,
           num_fields += chemistry_write_sparticles(sparts, list + num_fields);
           num_fields += tracers_write_sparticles(sparts, list + num_fields,
                                                  with_cosmology);
+          num_fields += fof_write_sparts(sparts, list + num_fields);
           if (with_stf) {
             num_fields += velociraptor_write_sparts(sparts, list + num_fields);
           }
@@ -1098,6 +1105,7 @@ void write_output_single(struct engine* e, const char* baseName,
               chemistry_write_sparticles(sparts_written, list + num_fields);
           num_fields += tracers_write_sparticles(
               sparts_written, list + num_fields, with_cosmology);
+          num_fields += fof_write_sparts(sparts_written, list + num_fields);
           if (with_stf) {
             num_fields +=
                 velociraptor_write_sparts(sparts_written, list + num_fields);
@@ -1112,7 +1120,7 @@ void write_output_single(struct engine* e, const char* baseName,
           N = Nblackholes;
           black_holes_write_particles(bparts, list, &num_fields);
           num_fields += chemistry_write_bparticles(bparts, list + num_fields);
-
+          num_fields += fof_write_bparts(bparts, list + num_fields);
           if (with_stf) {
             num_fields += velociraptor_write_bparts(bparts, list + num_fields);
           }
@@ -1135,6 +1143,7 @@ void write_output_single(struct engine* e, const char* baseName,
           black_holes_write_particles(bparts_written, list, &num_fields);
           num_fields +=
               chemistry_write_bparticles(bparts_written, list + num_fields);
+          num_fields += fof_write_bparts(bparts_written, list + num_fields);
           if (with_stf) {
             num_fields +=
                 velociraptor_write_bparts(bparts_written, list + num_fields);