diff --git a/.gitignore b/.gitignore
index 272c5328c768bb2e021cf9be22496ce9f8b8a2e7..1fd274009591c7ac98d71e05d536c98d4c17c485 100644
--- a/.gitignore
+++ b/.gitignore
@@ -40,12 +40,14 @@ examples/*/*/partition_fixed_costs.h
 examples/*/*/memuse_report-step*.dat
 examples/*/*/memuse_report-step*.log
 examples/*/*/restart/*
-examples/*/*/used_parameters.yml
 examples/*/stf_output*
 examples/*/stf.*
 examples/*/fof_output*
 examples/*/log*
+examples/*/*/used_parameters.yml
 examples/*/*/unused_parameters.yml
+examples/*/*/fof_used_parameters.yml
+examples/*/*/fof_unused_parameters.yml
 examples/*/*.mpg
 examples/*/*/gravity_checks_*.dat
 examples/*/*/coolingtables.tar.gz
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/fof.c b/src/fof.c
index 2f374ed1b1c4236b5cb5dcd67e1a2eaef3f6c939..de32d06d956db334a286414221781d7de0c5569c 100644
--- a/src/fof.c
+++ b/src/fof.c
@@ -236,7 +236,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 +968,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 +998,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 +1271,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 +1362,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 +1370,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 +1417,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 +1425,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 +1517,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 +1528,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 +1565,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 +1665,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 +1954,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 +2098,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 +2541,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 +2551,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 +2647,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 +2661,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 +2677,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)
diff --git a/src/fof_io.h b/src/fof_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..9863b02ed41571bcfdf421d17cf57582a2ac6345
--- /dev/null
+++ b/src/fof_io.h
@@ -0,0 +1,117 @@
+/*******************************************************************************
+ * 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"
+
+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;
+}
+
+/**
+ * @brief Specifies which FoF-related particle fields to write to a dataset
+ *
+ * @param parts The particle array.
+ * @param parts 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) {
+
+  list[0] = io_make_output_field_convert_part("GroupIDs", LONGLONG, 1,
+                                              UNIT_CONV_NO_UNITS, parts, xparts,
+                                              convert_part_group_id);
+
+  return 1;
+}
+
+/**
+ * @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) {
+
+  list[0] = io_make_output_field("GroupIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS,
+                                 gparts, fof_data.group_id);
+
+  return 1;
+}
+
+/**
+ * @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) {
+
+  list[0] = io_make_output_field_convert_spart("GroupIDs", LONGLONG, 1,
+                                               UNIT_CONV_NO_UNITS, sparts,
+                                               convert_spart_group_id);
+
+  return 1;
+}
+
+/**
+ * @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) {
+
+  list[0] = io_make_output_field_convert_bpart("GroupIDs", LONGLONG, 1,
+                                               UNIT_CONV_NO_UNITS, bparts,
+                                               convert_bpart_group_id);
+
+  return 1;
+}
+
+#endif /* SWIFT_FOF_IO_H */
diff --git a/src/fof_struct.h b/src/fof_struct.h
new file mode 100644
index 0000000000000000000000000000000000000000..a4f7f327073db7219a6dd60d48daf4cc4fe7ad7a
--- /dev/null
+++ b/src/fof_struct.h
@@ -0,0 +1,37 @@
+/*******************************************************************************
+ * 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"
+
+/**
+ * @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;
+};
+
+#endif /* SWIFT_FOF_STRUCT_H */
diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h
index 2e443e8fdc2479f3f2feff30281dccad9a7b6519..2150b3c445eb8fbd109e209508269e81d2e773c3 100644
--- a/src/gravity/Default/gravity_io.h
+++ b/src/gravity/Default/gravity_io.h
@@ -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 db1fbe2426a3e858dd429aa495641a4e0d4ed631..0585ba6f057d4efe1bf78efce375b43ba39930a4 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/serial_io.c b/src/serial_io.c
index c05c69e1e86b737e1e7f06b995f89f6bb548901a..936b30da7cf24d95f3de9a7869a37f4931d1383a 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 7a7856d29ce35435f750f2f71d66e0269d114261..819566df7e15f70e735d207c4f73dc4001983852 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);