diff --git a/doc/RTD/source/ParameterFiles/index.rst b/doc/RTD/source/ParameterFiles/index.rst
index d62cba2886697232c24b60cd649f6c53357c8a3b..0e0a6249cd942f5ffd3d856789b279e468d91023 100644
--- a/doc/RTD/source/ParameterFiles/index.rst
+++ b/doc/RTD/source/ParameterFiles/index.rst
@@ -414,12 +414,13 @@ parameter is the base name that will be used for all the outputs in the run:
 
 * The base name of the HDF5 snapshots: ``basename``.
 
-This name will then be appended by an under-score and 6 digits followed by
-``.hdf5`` (e.g. ``base_name_001234.hdf5``). The 6 digits are used to label the
-different outputs, starting at ``000000``. In the default setup the digits
-simply increase by one for each snapshot. However, if the optional parameter
-``int_time_label_on`` is switched on, then the 6-digits will the physical time
-of the simulation rounded to the nearest integer [#f3]_. 
+This name will then be appended by an under-score and 4 digits followed by
+``.hdf5`` (e.g. ``base_name_1234.hdf5``). The 4 digits are used to label the
+different outputs, starting at ``0000``. In the default setup the digits simply
+increase by one for each snapshot. However, if the optional parameter
+``int_time_label_on`` is switched on, then we use 6 digits and these will the
+physical time of the simulation rounded to the nearest integer
+(e.g. ``base_name_001234.hdf5``) [#f3]_.
 
 The time of the first snapshot is controlled by the two following options:
 
@@ -434,11 +435,28 @@ in the internal units of time. Users also have to provide the difference in time
 * Time difference between consecutive outputs: ``delta_time``.
 
 In non-cosmological runs this is also expressed in internal units. For
-cosmological runs, this value is *multiplied* to obtain the scale-factor of the
-next snapshot. This implies that the outputs are equally space in
-:math:`\log(a)` (See :ref:`Output_list_label` to have snapshots not regularly spaced in time).
-
-Users can optionally specify the level of compression used by the HDF5 libary
+cosmological runs, this value is *multiplied* to obtain the
+scale-factor of the next snapshot. This implies that the outputs are
+equally space in :math:`\log(a)` (See :ref:`Output_list_label` to have
+snapshots not regularly spaced in time).
+
+When running the code with structure finding activated, it is often
+useful to have a structure catalog written at the same simulation time
+as the snapshots. To activate this, the following parameter can be
+switched on:
+
+* Run VELOCIraptor every time a snapshot is dumped: ``invoke_stf``
+  (default: ``0``).
+
+This produces catalogs using the options specified for the stand-alone
+VELOCIraptor outputs (see the section :ref:`Parameters_structure_finding`) but
+with a base name and output number that matches the snapshot name
+(e.g. ``stf_base_name_1234.hdf5``) irrespective of the name specified in the
+section dedicated to VELOCIraptor. Note that the invocation of VELOCIraptor at
+every dump is done additionally to the stand-alone dumps that can be specified
+in the corresponding section of the YAML parameter file.
+
+Users can optionally specify the level of compression used by the HDF5 library
 using the parameter:
 
 * GZIP compression level of the HDF5 arrays: ``compression`` (default: ``0``).
@@ -464,14 +482,16 @@ When un-specified, these all take the same value as assumed by the internal
 system of units. These are rarely used but can offer a practical alternative to
 converting data in the post-processing of the simulations. 
 
-For a standard cosmological run, the full section would be:
+For a standard cosmological run with structure finding activated, the
+full section would be:
 
 .. code:: YAML
 
    Snapshots:
      basename:            output
      scale_factor_first:  0.02    # z = 49
-     delta_time:          1.02    
+     delta_time:          1.02
+     invoke_stf:          1
 	    
 Showing all the parameters for a basic hydro test-case, one would have:
 
@@ -481,6 +501,7 @@ Showing all the parameters for a basic hydro test-case, one would have:
      basename:            sedov
      time_first:          0.01
      delta_time:          0.005
+     invoke_stf:          0
      int_time_label_on:   0
      compression:         3
      UnitLength_in_cgs:   1.  # Use cm in outputs
@@ -598,6 +619,12 @@ Scheduler
 Domain Decomposition
 --------------------
 
+.. _Parameters_structure_finding:
+
+Structure finding (VELOCIraptor)
+--------------------------------
+
+
 .. [#f1] The thorough reader (or overly keen SWIFT tester) would find  that the speed of light is :math:`c=1.8026\times10^{12}\,\rm{fur}\,\rm{ftn}^{-1}`, Newton's constant becomes :math:`G_N=4.896735\times10^{-4}~\rm{fur}^3\,\rm{fir}^{-1}\,\rm{ftn}^{-2}` and Planck's constant turns into :math:`h=4.851453\times 10^{-34}~\rm{fur}^2\,\rm{fir}\,\rm{ftn}^{-1}`.
 
 
diff --git a/doc/RTD/source/VELOCIraptorInterface/stfwithswift.rst b/doc/RTD/source/VELOCIraptorInterface/stfwithswift.rst
index a663c37f93a6cede8c4528583c44183059414432..ed261b76abbcefaf5643a69069bb4b8ea1a0894c 100644
--- a/doc/RTD/source/VELOCIraptorInterface/stfwithswift.rst
+++ b/doc/RTD/source/VELOCIraptorInterface/stfwithswift.rst
@@ -50,8 +50,10 @@ HDF5 library, not a parallel build.
 Compiling SWIFT
 ---------------
 The next part is compiling SWIFT with VELOCIraptor and assumes you already
-downloaded SWIFT from the GitLab_, this can be done by running::
+downloaded SWIFT from the GitLab_, this can be done by running
 
+.. code:: bash
+  
   ./autogen.sh 
   ./configure --with-velociraptor=/path/to/VELOCIraptor-STF/src 
   make 
@@ -60,16 +62,16 @@ In which ``./autogen.sh`` only needs to be run once after the code is cloned
 from the GitLab_, and ``/path/to/`` is the path to the ``VELOCIraptor-STF``
 directory on your machine. In general ``./configure`` can be run with other
 options as desired. After this we can run SWIFT with VELOCIraptor, but for this
-we first need to add several lines to the yaml file of our simulation::
+we first need to add several lines to the yaml file of our simulation
 
     
-  #structure finding options
-  StructureFinding:
-  config_file_name:     stf_input_6dfof_dmonly_sub.cfg
-  basename:             ./stf
-  output_time_format:   1
-  scale_factor_first:   0.02
-  delta_time:           1.02
+.. code:: YAML
+
+   StructureFinding:      
+     config_file_name:     stf_input_6dfof_dmonly_sub.cfg
+     basename:             ./stf
+     scale_factor_first:   0.02
+     delta_time:           1.02
 
 In which we specify the ``.cfg`` file that is used by VELOCIraptor and the 
 other parameters which SWIFT needs to use. In the case of 
diff --git a/examples/DwarfGalaxy/dwarf_galaxy.yml b/examples/DwarfGalaxy/dwarf_galaxy.yml
index 0d815a99c42249bcbbdaf21dbaa34a55f61731aa..4c5e2a82b017725929138de011b1f3ed1fe9f1ef 100644
--- a/examples/DwarfGalaxy/dwarf_galaxy.yml
+++ b/examples/DwarfGalaxy/dwarf_galaxy.yml
@@ -10,10 +10,8 @@ InternalUnitSystem:
 StructureFinding:
   config_file_name:     stf_input.cfg # Name of the STF config file.
   basename:             ./stf         # Common part of the name of output files.
-  output_time_format:   0             # Specifies the frequency format of structure finding. 0 for simulation steps (delta_step) and 1 for simulation time intervals (delta_time).
   scale_factor_first:   0.92          # Scale-factor of the first snaphot (cosmological run)
   time_first:           0.01        # Time of the first structure finding output (in internal units).
-  delta_step:           1000          # Time difference between consecutive structure finding outputs (in internal units) in simulation steps.
   delta_time:           1.10          # Time difference between consecutive structure finding outputs (in internal units) in simulation time intervals.
 
 # Cosmological parameters
diff --git a/examples/EAGLE_25/eagle_25.yml b/examples/EAGLE_25/eagle_25.yml
index 0aec970db486a164696b23fdc1e281fbe4853486..cab0dbcd5efc0528ddc65a6dde1e5c2d7cb6b9a9 100644
--- a/examples/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_25/eagle_25.yml
@@ -10,10 +10,8 @@ InternalUnitSystem:
 StructureFinding:
   config_file_name:     stf_input.cfg    # Name of the STF config file.
   basename:             ./stf         # Common part of the name of output files.
-  output_time_format:   0             # Specifies the frequency format of structure finding. 0 for simulation steps (delta_step) and 1 for simulation time intervals (delta_time).
   scale_factor_first:   0.92          # Scale-factor of the first snaphot (cosmological run)
   time_first:           0.01        # Time of the first structure finding output (in internal units).
-  delta_step:           1000          # Time difference between consecutive structure finding outputs (in internal units) in simulation steps.
   delta_time:           1.10          # Time difference between consecutive structure finding outputs (in internal units) in simulation time intervals.
 
 # Cosmological parameters
diff --git a/examples/EAGLE_6/eagle_6.yml b/examples/EAGLE_6/eagle_6.yml
index 7c64c1cdedb6c8e9714471f4bad9611f548d05fa..e80fac8167a832c17cd10e1d2ae7cd854f314d17 100644
--- a/examples/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_6/eagle_6.yml
@@ -10,10 +10,8 @@ InternalUnitSystem:
 StructureFinding:
   config_file_name:     stf_input.cfg # Name of the STF config file.
   basename:             ./stf         # Common part of the name of output files.
-  output_time_format:   0             # Specifies the frequency format of structure finding. 0 for simulation steps (delta_step) and 1 for simulation time intervals (delta_time).
   scale_factor_first:   0.92          # Scale-factor of the first snaphot (cosmological run)
   time_first:           0.01        # Time of the first structure finding output (in internal units).
-  delta_step:           1000          # Time difference between consecutive structure finding outputs (in internal units) in simulation steps.
   delta_time:           1.10          # Time difference between consecutive structure finding outputs (in internal units) in simulation time intervals.
 
 # Cosmological parameters
diff --git a/examples/SmallCosmoVolume_DM/small_cosmo_volume_dm.yml b/examples/SmallCosmoVolume_DM/small_cosmo_volume_dm.yml
index 910137edc442c994a9f31a8c62e16818ca4ae97d..ebe3a78ee0d03eb53752b1dfa8fa749931a754a9 100644
--- a/examples/SmallCosmoVolume_DM/small_cosmo_volume_dm.yml
+++ b/examples/SmallCosmoVolume_DM/small_cosmo_volume_dm.yml
@@ -10,7 +10,6 @@ InternalUnitSystem:
 StructureFinding:
   config_file_name:     stf_input_6dfof_dmonly_sub.cfg
   basename:             ./stf
-  output_time_format:   1
   scale_factor_first:   0.02
   delta_time:           1.02
 
diff --git a/examples/SmallCosmoVolume_VELOCIraptor/small_cosmo_volume.yml b/examples/SmallCosmoVolume_VELOCIraptor/small_cosmo_volume.yml
index c8157a7a0e0065b1f58667fb8437b9e3883eda75..d6b9a78fe3c2a891492affbdea9787d62916d3ed 100644
--- a/examples/SmallCosmoVolume_VELOCIraptor/small_cosmo_volume.yml
+++ b/examples/SmallCosmoVolume_VELOCIraptor/small_cosmo_volume.yml
@@ -37,8 +37,9 @@ SPH:
 # Parameters governing the snapshots
 Snapshots:
   basename:            snap
-  delta_time:          1.02
+  delta_time:          1.05
   scale_factor_first:  0.02
+  invoke_stf:          1
   
 # Parameters governing the conserved quantities statistics
 Statistics:
@@ -52,16 +53,16 @@ Scheduler:
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  small_cosmo_volume.hdf5
+  periodic:                    1
   cleanup_h_factors:           1    
   cleanup_velocity_factors:    1  
-  generate_gas_in_ics: 1            # Generate gas particles from the DM-only ICs
-  cleanup_smoothing_lengths: 1      # Since we generate gas, make use of the (expensive) cleaning-up procedure.
+  generate_gas_in_ics:         1     # Generate gas particles from the DM-only ICs
+  cleanup_smoothing_lengths:   1     # Since we generate gas, make use of the (expensive) cleaning-up procedure.
 
 # Structure finding options (requires velociraptor)
 StructureFinding:
   config_file_name:     stfconfig_input.cfg
   basename:             ./stf
-  output_time_format:   1
   scale_factor_first:   0.02
   delta_time:           1.02
  
diff --git a/examples/main.c b/examples/main.c
index 9c6fd9655b172763502e0ac713ee3c5ce3b29d97..6d3ec54878034cfff8e989d5a30dbb4ea299b5df 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -919,6 +919,10 @@ int main(int argc, char *argv[]) {
       fflush(stdout);
     }
 
+#ifdef HAVE_VELOCIRAPTOR
+    if (with_structure_finding) velociraptor_init(&e);
+#endif
+
     /* Get some info to the user. */
     if (myrank == 0) {
       long long N_DM = N_total[1] - N_total[2] - N_total[0];
@@ -1201,14 +1205,6 @@ int main(int argc, char *argv[]) {
 #endif
     // write a final snapshot with logger, in order to facilitate a restart
     engine_dump_snapshot(&e);
-
-#ifdef HAVE_VELOCIRAPTOR
-    /* Call VELOCIraptor at the end of the run to find groups. */
-    if (e.policy & engine_policy_structure_finding) {
-      velociraptor_init(&e);
-      velociraptor_invoke(&e);
-    }
-#endif
   }
 
 #ifdef WITH_MPI
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 2fcc1545958346ed914e1feb52bdb9a829efaa07..75f7d77ac8659d8c0d8a75f622691a4ca21ae772 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -85,6 +85,7 @@ Snapshots:
   scale_factor_first: 0.1 # (Optional) Scale-factor of the first snapshot if cosmological time-integration.
   time_first: 0.          # (Optional) Time of the first output if non-cosmological time-integration (in internal units)
   delta_time: 0.01        # Time difference between consecutive outputs (in internal units)
+  invoke_stf: 0           # (Optional) Call VELOCIraptor every time a snapshot is written irrespective of the VELOCIraptor output strategy.
   compression: 0          # (Optional) Set the level of compression of the HDF5 datasets [0-9]. 0 does no compression.
   int_time_label_on:   0  # (Optional) Enable to label the snapshots using the time rounded to an integer (in internal units)
   UnitMass_in_cgs:     1  # (Optional) Unit system for the outputs (Grams)
@@ -157,6 +158,16 @@ DomainDecomposition:
   itr:              100     # When adaptive defines the ratio of inter node communication time to data redistribution time, in the range 0.00001 to 10000000.0.
                             # Lower values give less data movement during redistributions, at the cost of global balance which may require more communication.
 
+# Structure finding options (requires velociraptor)
+StructureFinding:
+  config_file_name:     stf_input.cfg # Name of the STF config file.
+  basename:             ./stf         # Common part of the name of output files.
+  scale_factor_first:   0.92          # (Optional) Scale-factor of the first snaphot (cosmological run)
+  time_first:           0.01          # (Optional) Time of the first structure finding output (in internal units).
+  delta_time:           1.10          # (Optional) Time difference between consecutive structure finding outputs (in internal units) in simulation time intervals.
+  output_list_on:       0   	      # (Optional) Enable the output list
+  output_list:          stflist.txt   # (Optional) File containing the output times (see documentation in "Parameter File" section)
+
 # Parameters related to the equation of state ------------------------------------------
 
 EoS:
@@ -287,15 +298,3 @@ EAGLEChemistry:
   init_abundance_Magnesium: 0.000        # Inital fraction of particle mass in Magnesium
   init_abundance_Silicon:   0.000        # Inital fraction of particle mass in Silicon
   init_abundance_Iron:      0.000        # Inital fraction of particle mass in Iron
-
-# Structure finding options (requires velociraptor)
-StructureFinding:
-  config_file_name:     stf_input.cfg # Name of the STF config file.
-  basename:             ./stf         # Common part of the name of output files.
-  output_time_format:   0             # Specifies the frequency format of structure finding. 0 for simulation steps (delta_step) and 1 for simulation time intervals (delta_time).
-  scale_factor_first:   0.92          # Scale-factor of the first snaphot (cosmological run)
-  time_first:           0.01          # Time of the first structure finding output (in internal units).
-  delta_step:           1000          # Time difference between consecutive structure finding outputs (in internal units) in simulation steps.
-  delta_time:           1.10          # Time difference between consecutive structure finding outputs (in internal units) in simulation time intervals.
-  output_list_on:      0   	      # (Optional) Enable the output list
-  output_list:         stflist.txt    # (Optional) File containing the output times (see documentation in "Parameter File" section)
diff --git a/src/Makefile.am b/src/Makefile.am
index db8cc00ed7dbb7d0b3b6812f06ce3daab8b9920e..6afe1c8509b77bc6eb5a2a22111673e1a17aaf15 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
+    logger_io.h tracers_io.h tracers.h tracers_struct.h velociraptor_struct.h velociraptor_io.h
 
 # source files for EAGLE cooling
 EAGLE_COOLING_SOURCES =
diff --git a/src/common_io.c b/src/common_io.c
index e6c0d016f4651e08a745ab0d399c6205c1f68415..733cf1dacac5f0c73ea401a584e2aa40eadd4a23 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -808,6 +808,28 @@ void io_convert_part_d_mapper(void* restrict temp, int N,
                          &temp_d[i * dim]);
 }
 
+/**
+ * @brief Mapper function to copy #part into a buffer of doubles using a
+ * conversion function.
+ */
+void io_convert_part_l_mapper(void* restrict temp, int N,
+                              void* restrict extra_data) {
+
+  const struct io_props props = *((const struct io_props*)extra_data);
+  const struct part* restrict parts = props.parts;
+  const struct xpart* restrict xparts = props.xparts;
+  const struct engine* e = props.e;
+  const size_t dim = props.dimension;
+
+  /* How far are we with this chunk? */
+  long long* restrict temp_l = (long long*)temp;
+  const ptrdiff_t delta = (temp_l - props.start_temp_l) / dim;
+
+  for (int i = 0; i < N; i++)
+    props.convert_part_l(e, parts + delta + i, xparts + delta + i,
+                         &temp_l[i * dim]);
+}
+
 /**
  * @brief Mapper function to copy #gpart into a buffer of floats using a
  * conversion function.
@@ -848,6 +870,26 @@ void io_convert_gpart_d_mapper(void* restrict temp, int N,
     props.convert_gpart_d(e, gparts + delta + i, &temp_d[i * dim]);
 }
 
+/**
+ * @brief Mapper function to copy #gpart into a buffer of doubles using a
+ * conversion function.
+ */
+void io_convert_gpart_l_mapper(void* restrict temp, int N,
+                               void* restrict extra_data) {
+
+  const struct io_props props = *((const struct io_props*)extra_data);
+  const struct gpart* restrict gparts = props.gparts;
+  const struct engine* e = props.e;
+  const size_t dim = props.dimension;
+
+  /* How far are we with this chunk? */
+  long long* restrict temp_l = (long long*)temp;
+  const ptrdiff_t delta = (temp_l - props.start_temp_l) / dim;
+
+  for (int i = 0; i < N; i++)
+    props.convert_gpart_l(e, gparts + delta + i, &temp_l[i * dim]);
+}
+
 /**
  * @brief Mapper function to copy #spart into a buffer of floats using a
  * conversion function.
@@ -888,6 +930,26 @@ void io_convert_spart_d_mapper(void* restrict temp, int N,
     props.convert_spart_d(e, sparts + delta + i, &temp_d[i * dim]);
 }
 
+/**
+ * @brief Mapper function to copy #spart into a buffer of doubles using a
+ * conversion function.
+ */
+void io_convert_spart_l_mapper(void* restrict temp, int N,
+                               void* restrict extra_data) {
+
+  const struct io_props props = *((const struct io_props*)extra_data);
+  const struct spart* restrict sparts = props.sparts;
+  const struct engine* e = props.e;
+  const size_t dim = props.dimension;
+
+  /* How far are we with this chunk? */
+  long long* restrict temp_l = (long long*)temp;
+  const ptrdiff_t delta = (temp_l - props.start_temp_l) / dim;
+
+  for (int i = 0; i < N; i++)
+    props.convert_spart_l(e, sparts + delta + i, &temp_l[i * dim]);
+}
+
 /**
  * @brief Copy the particle data into a temporary buffer ready for i/o.
  *
@@ -945,6 +1007,18 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_part_d_mapper, temp_d, N, copySize, 0,
                      (void*)&props);
 
+    } else if (props.convert_part_l != NULL) {
+
+      /* Prepare some parameters */
+      long long* temp_l = (long long*)temp;
+      props.start_temp_l = (long long*)temp;
+      props.e = e;
+
+      /* Copy the whole thing into a buffer */
+      threadpool_map((struct threadpool*)&e->threadpool,
+                     io_convert_part_l_mapper, temp_l, N, copySize, 0,
+                     (void*)&props);
+
     } else if (props.convert_gpart_f != NULL) {
 
       /* Prepare some parameters */
@@ -969,6 +1043,18 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_gpart_d_mapper, temp_d, N, copySize, 0,
                      (void*)&props);
 
+    } else if (props.convert_gpart_l != NULL) {
+
+      /* Prepare some parameters */
+      long long* temp_l = (long long*)temp;
+      props.start_temp_l = (long long*)temp;
+      props.e = e;
+
+      /* Copy the whole thing into a buffer */
+      threadpool_map((struct threadpool*)&e->threadpool,
+                     io_convert_gpart_l_mapper, temp_l, N, copySize, 0,
+                     (void*)&props);
+
     } else if (props.convert_spart_f != NULL) {
 
       /* Prepare some parameters */
@@ -993,6 +1079,18 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_spart_d_mapper, temp_d, N, copySize, 0,
                      (void*)&props);
 
+    } else if (props.convert_spart_l != NULL) {
+
+      /* Prepare some parameters */
+      long long* temp_l = (long long*)temp;
+      props.start_temp_l = (long long*)temp;
+      props.e = e;
+
+      /* Copy the whole thing into a buffer */
+      threadpool_map((struct threadpool*)&e->threadpool,
+                     io_convert_spart_l_mapper, temp_l, N, copySize, 0,
+                     (void*)&props);
+
     } else {
       error("Missing conversion function");
     }
@@ -1254,15 +1352,21 @@ void io_collect_sparts_to_write(const struct spart* restrict sparts,
  * @brief Copy every non-inhibited DM #gpart into the gparts_written array.
  *
  * @param gparts The array of #gpart containing all particles.
+ * @param vr_data The array of gpart-related VELOCIraptor output.
  * @param gparts_written The array of #gpart to fill with particles we want to
  * write.
+ * @param vr_data_written The array of gpart-related VELOCIraptor with particles
+ * we want to write.
  * @param Ngparts The total number of #part.
  * @param Ngparts_written The total number of #part to write.
+ * @param with_stf Are we running with STF? i.e. do we want to collect vr data?
  */
-void io_collect_gparts_to_write(const struct gpart* restrict gparts,
-                                struct gpart* restrict gparts_written,
-                                const size_t Ngparts,
-                                const size_t Ngparts_written) {
+void io_collect_gparts_to_write(
+    const struct gpart* restrict gparts,
+    const struct velociraptor_gpart_data* restrict vr_data,
+    struct gpart* restrict gparts_written,
+    struct velociraptor_gpart_data* restrict vr_data_written,
+    const size_t Ngparts, const size_t Ngparts_written, const int with_stf) {
 
   size_t count = 0;
 
@@ -1274,6 +1378,8 @@ void io_collect_gparts_to_write(const struct gpart* restrict gparts,
         (gparts[i].time_bin != time_bin_not_created) &&
         (gparts[i].type == swift_type_dark_matter)) {
 
+      if (with_stf) vr_data_written[count] = vr_data[i];
+
       gparts_written[count] = gparts[i];
       count++;
     }
@@ -1281,7 +1387,7 @@ void io_collect_gparts_to_write(const struct gpart* restrict gparts,
 
   /* Check that everything is fine */
   if (count != Ngparts_written)
-    error("Collected the wrong number of s-particles (%zu vs. %zu expected)",
+    error("Collected the wrong number of g-particles (%zu vs. %zu expected)",
           count, Ngparts_written);
 }
 
diff --git a/src/common_io.h b/src/common_io.h
index bdff3e37d11c59a71c08dcf3e19758019adc3a73..eb1ee0a804f324d897842fb2a0ca33fc07e769d6 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -36,6 +36,7 @@
 struct cell;
 struct part;
 struct gpart;
+struct velociraptor_gpart_data;
 struct spart;
 struct xpart;
 struct io_props;
@@ -113,9 +114,11 @@ void io_collect_sparts_to_write(const struct spart* restrict sparts,
                                 const size_t Nsparts,
                                 const size_t Nsparts_written);
 void io_collect_gparts_to_write(const struct gpart* restrict gparts,
+                                const struct velociraptor_gpart_data* vr_data,
                                 struct gpart* restrict gparts_written,
+                                struct velociraptor_gpart_data* vr_data_written,
                                 const size_t Ngparts,
-                                const size_t Ngparts_written);
+                                const size_t Ngparts_written, int with_stf);
 void io_prepare_dm_gparts(struct threadpool* tp, struct gpart* const gparts,
                           size_t Ndm);
 void io_duplicate_hydro_gparts(struct threadpool* tp, struct part* const parts,
diff --git a/src/engine.c b/src/engine.c
index 9354da36ade1211209e8482e2c96bc93cab5213a..8899c37a604880164b06329308476cd7747ee9b2 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3169,90 +3169,87 @@ void engine_step(struct engine *e) {
  */
 void engine_check_for_dumps(struct engine *e) {
 
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
   const int with_stf = (e->policy & engine_policy_structure_finding);
-  const int stf_time_output = (e->stf_output_freq_format == io_stf_time);
+
+  /* What kind of output are we getting? */
+  enum output_type {
+    output_none,
+    output_snapshot,
+    output_statistics,
+    output_stf
+  };
+
+  /* What kind of output do we want? And at which time ?
+   * Find the earliest output (amongst all kinds) that takes place
+   * before the next time-step */
+  enum output_type type = output_none;
+  integertime_t ti_output = max_nr_timesteps;
 
   /* Save some statistics ? */
-  int save_stats = 0;
-  if (e->ti_end_min > e->ti_next_stats && e->ti_next_stats > 0) save_stats = 1;
+  if (e->ti_end_min > e->ti_next_stats && e->ti_next_stats > 0) {
+    if (e->ti_next_stats < ti_output) {
+      ti_output = e->ti_next_stats;
+      type = output_statistics;
+    }
+  }
 
   /* Do we want a snapshot? */
-  int dump_snapshot = 0;
-  if (e->ti_end_min > e->ti_next_snapshot && e->ti_next_snapshot > 0)
-    dump_snapshot = 1;
+  if (e->ti_end_min > e->ti_next_snapshot && e->ti_next_snapshot > 0) {
+    if (e->ti_next_snapshot < ti_output) {
+      ti_output = e->ti_next_snapshot;
+      type = output_snapshot;
+    }
+  }
 
   /* Do we want to perform structure finding? */
-  int run_stf = 0;
-  if (with_stf && stf_time_output) {
-    if (e->ti_end_min > e->ti_next_stf && e->ti_next_stf > 0) run_stf = 1;
-  }
-  if (with_stf && !stf_time_output) {
-    if (e->step % e->delta_step_stf == 0) run_stf = 1;
+  if (with_stf) {
+    if (e->ti_end_min > e->ti_next_stf && e->ti_next_stf > 0) {
+      if (e->ti_next_stf < ti_output) {
+        ti_output = e->ti_next_stf;
+        type = output_stf;
+      }
+    }
   }
 
   /* Store information before attempting extra dump-related drifts */
-  integertime_t ti_current = e->ti_current;
-  timebin_t max_active_bin = e->max_active_bin;
-  double time = e->time;
+  const integertime_t ti_current = e->ti_current;
+  const timebin_t max_active_bin = e->max_active_bin;
+  const double time = e->time;
+
+  while (type != output_none) {
+
+    /* Let's fake that we are at the dump time */
+    e->ti_current = ti_output;
+    e->max_active_bin = 0;
+    if (with_cosmology) {
+      cosmology_update(e->cosmology, e->physical_constants, e->ti_current);
+      e->time = e->cosmology->time;
+    } else {
+      e->time = ti_output * e->time_base + e->time_begin;
+    }
 
-  while (save_stats || dump_snapshot || run_stf) {
+    /* Drift everyone */
+    engine_drift_all(e, /*drift_mpole=*/0);
 
     /* Write some form of output */
-    if (dump_snapshot && save_stats) {
+    switch (type) {
+      case output_snapshot:
 
-      /* If both, need to figure out which one occurs first */
-      if (e->ti_next_stats == e->ti_next_snapshot) {
+        /* Do we want a corresponding VELOCIraptor output? */
+        if (with_stf && e->snapshot_invoke_stf) {
 
-        /* Let's fake that we are at the common dump time */
-        e->ti_current = e->ti_next_snapshot;
-        e->max_active_bin = 0;
-        if ((e->policy & engine_policy_cosmology)) {
-          cosmology_update(e->cosmology, e->physical_constants, e->ti_current);
-          e->time = e->cosmology->time;
-        } else {
-          e->time = e->ti_next_stats * e->time_base + e->time_begin;
-        }
-
-        /* Drift everyone */
-        engine_drift_all(e, /*drift_mpole=*/0);
-
-        /* Dump everything */
-        engine_print_stats(e);
-#ifdef WITH_LOGGER
-        /* Write a file containing the offsets in the particle logger. */
-        engine_dump_index(e);
+#ifdef HAVE_VELOCIRAPTOR
+          velociraptor_invoke(e, /*linked_with_snap=*/1);
+          e->step_props |= engine_step_prop_stf;
 #else
-        engine_dump_snapshot(e);
+          error(
+              "Asking for a VELOCIraptor output but SWIFT was compiled without "
+              "the interface!");
 #endif
-
-      } else if (e->ti_next_stats < e->ti_next_snapshot) {
-
-        /* Let's fake that we are at the stats dump time */
-        e->ti_current = e->ti_next_stats;
-        e->max_active_bin = 0;
-        if ((e->policy & engine_policy_cosmology)) {
-          cosmology_update(e->cosmology, e->physical_constants, e->ti_current);
-          e->time = e->cosmology->time;
-        } else {
-          e->time = e->ti_next_stats * e->time_base + e->time_begin;
         }
 
-        /* Drift everyone */
-        engine_drift_all(e, /*drift_mpole=*/0);
-
-        /* Dump stats */
-        engine_print_stats(e);
-
-        /* Let's fake that we are at the snapshot dump time */
-        e->ti_current = e->ti_next_snapshot;
-        e->max_active_bin = 0;
-        if (!(e->policy & engine_policy_cosmology))
-          e->time = e->ti_next_snapshot * e->time_base + e->time_begin;
-
-        /* Drift everyone */
-        engine_drift_all(e, /*drift_mpole=*/0);
-
-        /* Dump snapshot */
+          /* Dump... */
 #ifdef WITH_LOGGER
         /* Write a file containing the offsets in the particle logger. */
         engine_dump_index(e);
@@ -3260,144 +3257,81 @@ void engine_check_for_dumps(struct engine *e) {
         engine_dump_snapshot(e);
 #endif
 
-      } else if (e->ti_next_stats > e->ti_next_snapshot) {
-
-        /* Let's fake that we are at the snapshot dump time */
-        e->ti_current = e->ti_next_snapshot;
-        e->max_active_bin = 0;
-        if ((e->policy & engine_policy_cosmology)) {
-          cosmology_update(e->cosmology, e->physical_constants, e->ti_current);
-          e->time = e->cosmology->time;
-        } else {
-          e->time = e->ti_next_stats * e->time_base + e->time_begin;
-        }
-
-        /* Drift everyone */
-        engine_drift_all(e, /*drift_mpole=*/0);
-
-        /* Dump snapshot */
-#ifdef WITH_LOGGER
-        /* Write a file containing the offsets in the particle logger. */
-        engine_dump_index(e);
-#else
-        engine_dump_snapshot(e);
+        /* Free the memory allocated for VELOCIraptor i/o. */
+        if (with_stf && e->snapshot_invoke_stf) {
+#ifdef HAVE_VELOCIRAPTOR
+          free(e->s->gpart_group_data);
+          e->s->gpart_group_data = NULL;
 #endif
+        }
 
-        /* Let's fake that we are at the stats dump time */
-        e->ti_current = e->ti_next_stats;
-        e->max_active_bin = 0;
-        if (!(e->policy & engine_policy_cosmology))
-          e->time = e->ti_next_stats * e->time_base + e->time_begin;
+        /* ... and find the next output time */
+        engine_compute_next_snapshot_time(e);
+        break;
 
-        /* Drift everyone */
-        engine_drift_all(e, /*drift_mpole=*/0);
+      case output_statistics:
 
-        /* Dump stats */
+        /* Dump */
         engine_print_stats(e);
-      }
-
-      /* Let's compute the time of the next outputs */
-      engine_compute_next_snapshot_time(e);
-      engine_compute_next_statistics_time(e);
 
-    } else if (dump_snapshot) {
+        /* and move on */
+        engine_compute_next_statistics_time(e);
 
-      /* Let's fake that we are at the snapshot dump time */
-      e->ti_current = e->ti_next_snapshot;
-      e->max_active_bin = 0;
-      if ((e->policy & engine_policy_cosmology)) {
-        cosmology_update(e->cosmology, e->physical_constants, e->ti_current);
-        e->time = e->cosmology->time;
-      } else {
-        e->time = e->ti_next_stats * e->time_base + e->time_begin;
-      }
+        break;
 
-      /* Drift everyone */
-      engine_drift_all(e, /*drift_mpole=*/0);
-
-      /* Dump... */
-#ifdef WITH_LOGGER
-      /* Write a file containing the offsets in the particle logger. */
-      engine_dump_index(e);
-#else
-      engine_dump_snapshot(e);
-#endif
-
-      /* ... and find the next output time */
-      engine_compute_next_snapshot_time(e);
-
-    } else if (save_stats) {
-
-      /* Let's fake that we are at the stats dump time */
-      e->ti_current = e->ti_next_stats;
-      e->max_active_bin = 0;
-      if ((e->policy & engine_policy_cosmology)) {
-        cosmology_update(e->cosmology, e->physical_constants, e->ti_current);
-        e->time = e->cosmology->time;
-      } else {
-        e->time = e->ti_next_stats * e->time_base + e->time_begin;
-      }
-
-      /* Drift everyone */
-      engine_drift_all(e, /*drift_mpole=*/0);
-
-      /* Dump */
-      engine_print_stats(e);
-
-      /* and move on */
-      engine_compute_next_statistics_time(e);
-    }
-
-    /* Perform structure finding? */
-    if (run_stf) {
+      case output_stf:
 
 #ifdef HAVE_VELOCIRAPTOR
+        /* Unleash the raptor! */
+        velociraptor_invoke(e, /*linked_with_snap=*/0);
+        e->step_props |= engine_step_prop_stf;
 
-      // MATTHIEU: Check the order with the other i/o options.
-      if (!dump_snapshot && !save_stats) {
-
-        /* Let's fake that we are at the stats dump time */
-        e->ti_current = e->ti_next_stf;
-        e->max_active_bin = 0;
-        if ((e->policy & engine_policy_cosmology)) {
-          cosmology_update(e->cosmology, e->physical_constants, e->ti_current);
-          e->time = e->cosmology->time;
-        } else {
-          e->time = e->ti_next_stats * e->time_base + e->time_begin;
-        }
-
-        /* Drift everyone */
-        engine_drift_all(e, /*drift_mpole=*/0);
-      }
-
-      velociraptor_init(e);
-      velociraptor_invoke(e);
-
-      /* ... and find the next output time */
-      if (e->stf_output_freq_format == io_stf_time)
+        /* ... and find the next output time */
         engine_compute_next_stf_time(e);
+#else
+        error(
+            "Asking for a VELOCIraptor output but SWIFT was compiled without "
+            "the interface!");
 #endif
+        break;
+
+      default:
+        error("Invalid dump type");
     }
 
     /* We need to see whether whether we are in the pathological case
      * where there can be another dump before the next step. */
 
+    type = output_none;
+    ti_output = max_nr_timesteps;
+
     /* Save some statistics ? */
-    save_stats = 0;
-    if (e->ti_end_min > e->ti_next_stats && e->ti_next_stats > 0)
-      save_stats = 1;
+    if (e->ti_end_min > e->ti_next_stats && e->ti_next_stats > 0) {
+      if (e->ti_next_stats < ti_output) {
+        ti_output = e->ti_next_stats;
+        type = output_statistics;
+      }
+    }
 
     /* Do we want a snapshot? */
-    dump_snapshot = 0;
-    if (e->ti_end_min > e->ti_next_snapshot && e->ti_next_snapshot > 0)
-      dump_snapshot = 1;
+    if (e->ti_end_min > e->ti_next_snapshot && e->ti_next_snapshot > 0) {
+      if (e->ti_next_snapshot < ti_output) {
+        ti_output = e->ti_next_snapshot;
+        type = output_snapshot;
+      }
+    }
 
     /* Do we want to perform structure finding? */
-    run_stf = 0;
-    if (with_stf && stf_time_output) {
-      if (e->ti_end_min > e->ti_next_stf && e->ti_next_stf > 0) run_stf = 1;
+    if (with_stf) {
+      if (e->ti_end_min > e->ti_next_stf && e->ti_next_stf > 0) {
+        if (e->ti_next_stf < ti_output) {
+          ti_output = e->ti_next_stf;
+          type = output_stf;
+        }
+      }
     }
-  }
+
+  } /* While loop over output types */
 
   /* Restore the information we stored */
   e->ti_current = ti_current;
@@ -4157,9 +4091,12 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
       parser_get_opt_param_int(params, "Snapshots:compression", 0);
   e->snapshot_int_time_label_on =
       parser_get_opt_param_int(params, "Snapshots:int_time_label_on", 0);
+  e->snapshot_invoke_stf =
+      parser_get_opt_param_int(params, "Snapshots:invoke_stf", 0);
   e->snapshot_units = (struct unit_system *)malloc(sizeof(struct unit_system));
   units_init_default(e->snapshot_units, params, "Snapshots", internal_units);
   e->snapshot_output_count = 0;
+  e->stf_output_count = 0;
   e->dt_min = parser_get_param_double(params, "TimeIntegration:dt_min");
   e->dt_max = parser_get_param_double(params, "TimeIntegration:dt_max");
   e->dt_max_RMS_displacement = FLT_MAX;
@@ -4184,7 +4121,6 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
   e->cooling_func = cooling_func;
   e->chemistry = chemistry;
   e->parameter_file = params;
-  e->cell_loc = NULL;
 #ifdef WITH_MPI
   e->cputime_last_step = 0;
   e->last_repartition = 0;
@@ -4225,28 +4161,16 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
   /* Initialise VELOCIraptor output. */
   if (e->policy & engine_policy_structure_finding) {
     parser_get_param_string(params, "StructureFinding:basename",
-                            e->stfBaseName);
+                            e->stf_base_name);
+    parser_get_param_string(params, "StructureFinding:config_file_name",
+                            e->stf_config_file_name);
+
     e->time_first_stf_output =
         parser_get_opt_param_double(params, "StructureFinding:time_first", 0.);
     e->a_first_stf_output = parser_get_opt_param_double(
         params, "StructureFinding:scale_factor_first", 0.1);
-    e->stf_output_freq_format = (enum io_stf_output_format)parser_get_param_int(
-        params, "StructureFinding:output_time_format");
-
-    if (e->stf_output_freq_format == io_stf_steps) {
-      e->delta_step_stf =
-          parser_get_param_int(params, "StructureFinding:delta_step");
-    } else if (e->stf_output_freq_format == io_stf_time) {
-      e->delta_time_stf =
-          parser_get_param_double(params, "StructureFinding:delta_time");
-    } else {
-      error(
-          "Invalid flag (%d) set for output time format of structure finding.",
-          e->stf_output_freq_format);
-    }
-
-    /* overwrite input if outputlist */
-    if (e->output_list_stf) e->stf_output_freq_format = io_stf_time;
+    e->delta_time_stf =
+        parser_get_opt_param_double(params, "StructureFinding:delta_time", -1.);
   }
 
   engine_init_output_lists(e, params);
@@ -4492,10 +4416,11 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
 
       fprintf(e->file_timesteps,
               "# Step Properties: Rebuild=%d, Redistribute=%d, Repartition=%d, "
-              "Statistics=%d, Snapshot=%d, Restarts=%d\n",
+              "Statistics=%d, Snapshot=%d, Restarts=%d STF=%d, logger=%d\n",
               engine_step_prop_rebuild, engine_step_prop_redistribute,
               engine_step_prop_repartition, engine_step_prop_statistics,
-              engine_step_prop_snapshot, engine_step_prop_restarts);
+              engine_step_prop_snapshot, engine_step_prop_restarts,
+              engine_step_prop_stf, engine_step_prop_logger_index);
 
       fprintf(e->file_timesteps,
               "# %6s %14s %12s %12s %14s %9s %12s %12s %12s %16s [%s] %6s\n",
@@ -4586,17 +4511,18 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
           "simulation start a=%e.",
           e->a_first_statistics, e->cosmology->a_begin);
 
-    if ((e->policy & engine_policy_structure_finding) &&
-        (e->stf_output_freq_format == io_stf_time)) {
+    if (e->policy & engine_policy_structure_finding) {
 
-      if (e->delta_time_stf <= 1.)
+      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->delta_time_stf != -1.)
         error("Time between STF (%e) must be > 1.", e->delta_time_stf);
 
       if (e->a_first_stf_output < e->cosmology->a_begin)
         error(
             "Scale-factor of first stf output (%e) must be after the "
-            "simulation "
-            "start a=%e.",
+            "simulation start a=%e.",
             e->a_first_stf_output, e->cosmology->a_begin);
     }
   } else {
@@ -4622,10 +4548,12 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
           "t=%e.",
           e->time_first_statistics, e->time_begin);
 
-    if ((e->policy & engine_policy_structure_finding) &&
-        (e->stf_output_freq_format == io_stf_time)) {
+    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 <= 0.)
+      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)
@@ -4634,12 +4562,6 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
     }
   }
 
-  if (e->policy & engine_policy_structure_finding) {
-    /* Find the time of the first stf output */
-    if (e->stf_output_freq_format == io_stf_time)
-      engine_compute_next_stf_time(e);
-  }
-
   /* Get the total mass */
   e->total_mass = 0.;
   for (size_t i = 0; i < e->s->nr_gparts; ++i)
@@ -4664,6 +4586,19 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
   /* 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) {
+    engine_compute_next_stf_time(e);
+  }
+
+  /* Check that we are invoking VELOCIraptor only if we have it */
+  if (e->snapshot_invoke_stf &&
+      !(e->policy & engine_policy_structure_finding)) {
+    error(
+        "Invoking VELOCIraptor after snapshots but structure finding wasn't "
+        "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);
 
@@ -5253,7 +5188,6 @@ void engine_clean(struct engine *e) {
   output_list_clean(&e->output_list_stf);
 
   free(e->links);
-  free(e->cell_loc);
 #if defined(WITH_LOGGER)
   logger_clean(e->logger);
   free(e->logger);
diff --git a/src/engine.h b/src/engine.h
index 7e6cc7ffe49c90fbeffdaa0249507d04f2a35c78..4a51d6a9fe864c016f1de23ae16d0e0567a66187 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -91,7 +91,8 @@ enum engine_step_properties {
   engine_step_prop_statistics = (1 << 3),
   engine_step_prop_snapshot = (1 << 4),
   engine_step_prop_restarts = (1 << 5),
-  engine_step_prop_logger_index = (1 << 6)
+  engine_step_prop_stf = (1 << 6),
+  engine_step_prop_logger_index = (1 << 7)
 };
 
 /* Some constants */
@@ -225,9 +226,6 @@ struct engine {
   /* The internal system of units */
   const struct unit_system *internal_units;
 
-  /* Top-level cell locations for VELOCIraptor. */
-  struct cell_loc *cell_loc;
-
   /* Snapshot information */
   double a_first_snapshot;
   double time_first_snapshot;
@@ -242,12 +240,11 @@ struct engine {
   char snapshot_base_name[PARSER_MAX_LINE_SIZE];
   int snapshot_compression;
   int snapshot_int_time_label_on;
+  int snapshot_invoke_stf;
   struct unit_system *snapshot_units;
   int snapshot_output_count;
 
   /* Structure finding information */
-  enum io_stf_output_format stf_output_freq_format;
-  int delta_step_stf;
   double a_first_stf_output;
   double time_first_stf_output;
   double delta_time_stf;
@@ -258,7 +255,9 @@ struct engine {
   /* Integer time of the next stf output */
   integertime_t ti_next_stf;
 
-  char stfBaseName[PARSER_MAX_LINE_SIZE];
+  char stf_config_file_name[PARSER_MAX_LINE_SIZE];
+  char stf_base_name[PARSER_MAX_LINE_SIZE];
+  int stf_output_count;
 
   /* Statistics information */
   double a_first_statistics;
diff --git a/src/io_properties.h b/src/io_properties.h
index 9e948fc3991b0178d06fdd5d83fa900a98f84d2a..c45edb2641e374e2cfaec6c3251aff7d18f361d6 100644
--- a/src/io_properties.h
+++ b/src/io_properties.h
@@ -43,14 +43,23 @@ typedef void (*conversion_func_part_float)(const struct engine*,
 typedef void (*conversion_func_part_double)(const struct engine*,
                                             const struct part*,
                                             const struct xpart*, double*);
+typedef void (*conversion_func_part_long_long)(const struct engine*,
+                                               const struct part*,
+                                               const struct xpart*, long long*);
 typedef void (*conversion_func_gpart_float)(const struct engine*,
                                             const struct gpart*, float*);
 typedef void (*conversion_func_gpart_double)(const struct engine*,
                                              const struct gpart*, double*);
+typedef void (*conversion_func_gpart_long_long)(const struct engine*,
+                                                const struct gpart*,
+                                                long long*);
 typedef void (*conversion_func_spart_float)(const struct engine*,
                                             const struct spart*, float*);
 typedef void (*conversion_func_spart_double)(const struct engine*,
                                              const struct spart*, double*);
+typedef void (*conversion_func_spart_long_long)(const struct engine*,
+                                                const struct spart*,
+                                                long long*);
 
 /**
  * @brief The properties of a given dataset for i/o
@@ -79,6 +88,7 @@ struct io_props {
   char* start_temp_c;
   float* start_temp_f;
   double* start_temp_d;
+  long long* start_temp_l;
 
   /* Pointer to the engine */
   const struct engine* e;
@@ -98,14 +108,17 @@ struct io_props {
   /* Conversion function for part */
   conversion_func_part_float convert_part_f;
   conversion_func_part_double convert_part_d;
+  conversion_func_part_long_long convert_part_l;
 
   /* Conversion function for gpart */
   conversion_func_gpart_float convert_gpart_f;
   conversion_func_gpart_double convert_gpart_d;
+  conversion_func_gpart_long_long convert_gpart_l;
 
   /* Conversion function for spart */
   conversion_func_spart_float convert_spart_f;
   conversion_func_spart_double convert_spart_d;
+  conversion_func_spart_long_long convert_spart_l;
 };
 
 /**
@@ -147,10 +160,13 @@ INLINE static struct io_props io_make_input_field_(
   r.conversion = 0;
   r.convert_part_f = NULL;
   r.convert_part_d = NULL;
+  r.convert_part_l = NULL;
   r.convert_gpart_f = NULL;
   r.convert_gpart_d = NULL;
+  r.convert_gpart_l = NULL;
   r.convert_spart_f = NULL;
   r.convert_spart_d = NULL;
+  r.convert_spart_l = NULL;
 
   return r;
 }
@@ -191,10 +207,13 @@ INLINE static struct io_props io_make_output_field_(
   r.conversion = 0;
   r.convert_part_f = NULL;
   r.convert_part_d = NULL;
+  r.convert_part_l = NULL;
   r.convert_gpart_f = NULL;
   r.convert_gpart_d = NULL;
+  r.convert_gpart_l = NULL;
   r.convert_spart_f = NULL;
   r.convert_spart_d = NULL;
+  r.convert_spart_l = NULL;
 
   return r;
 }
@@ -242,10 +261,13 @@ INLINE static struct io_props io_make_output_field_convert_part_FLOAT(
   r.conversion = 1;
   r.convert_part_f = functionPtr;
   r.convert_part_d = NULL;
+  r.convert_part_l = NULL;
   r.convert_gpart_f = NULL;
   r.convert_gpart_d = NULL;
+  r.convert_gpart_l = NULL;
   r.convert_spart_f = NULL;
   r.convert_spart_d = NULL;
+  r.convert_spart_l = NULL;
 
   return r;
 }
@@ -285,10 +307,59 @@ INLINE static struct io_props io_make_output_field_convert_part_DOUBLE(
   r.conversion = 1;
   r.convert_part_f = NULL;
   r.convert_part_d = functionPtr;
+  r.convert_part_l = NULL;
   r.convert_gpart_f = NULL;
   r.convert_gpart_d = NULL;
+  r.convert_gpart_l = NULL;
   r.convert_spart_f = NULL;
   r.convert_spart_d = NULL;
+  r.convert_spart_l = NULL;
+
+  return r;
+}
+
+/**
+ * @brief Construct an #io_props from its parameters
+ *
+ * @param name Name of the field to read
+ * @param type The type of the data
+ * @param dimension Dataset dimension (1D, 3D, ...)
+ * @param units The units of the dataset
+ * @param partSize The size in byte of the particle
+ * @param parts The particle array
+ * @param xparts The xparticle array
+ * @param functionPtr The function used to convert a particle to a double
+ *
+ * Do not call this function directly. Use the macro defined above.
+ */
+INLINE static struct io_props io_make_output_field_convert_part_LONGLONG(
+    const char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension,
+    enum unit_conversion_factor units, size_t partSize,
+    const struct part* parts, const struct xpart* xparts,
+    conversion_func_part_long_long functionPtr) {
+
+  struct io_props r;
+  strcpy(r.name, name);
+  r.type = type;
+  r.dimension = dimension;
+  r.importance = UNUSED;
+  r.units = units;
+  r.field = NULL;
+  r.partSize = partSize;
+  r.parts = parts;
+  r.xparts = xparts;
+  r.gparts = NULL;
+  r.sparts = NULL;
+  r.conversion = 1;
+  r.convert_part_f = NULL;
+  r.convert_part_d = NULL;
+  r.convert_part_l = functionPtr;
+  r.convert_gpart_f = NULL;
+  r.convert_gpart_d = NULL;
+  r.convert_gpart_l = NULL;
+  r.convert_spart_f = NULL;
+  r.convert_spart_d = NULL;
+  r.convert_spart_l = NULL;
 
   return r;
 }
@@ -334,10 +405,13 @@ INLINE static struct io_props io_make_output_field_convert_gpart_FLOAT(
   r.conversion = 1;
   r.convert_part_f = NULL;
   r.convert_part_d = NULL;
+  r.convert_part_l = NULL;
   r.convert_gpart_f = functionPtr;
   r.convert_gpart_d = NULL;
+  r.convert_gpart_l = NULL;
   r.convert_spart_f = NULL;
   r.convert_spart_d = NULL;
+  r.convert_spart_l = NULL;
 
   return r;
 }
@@ -375,10 +449,57 @@ INLINE static struct io_props io_make_output_field_convert_gpart_DOUBLE(
   r.conversion = 1;
   r.convert_part_f = NULL;
   r.convert_part_d = NULL;
+  r.convert_part_l = NULL;
   r.convert_gpart_f = NULL;
   r.convert_gpart_d = functionPtr;
+  r.convert_gpart_l = NULL;
   r.convert_spart_f = NULL;
   r.convert_spart_d = NULL;
+  r.convert_spart_l = NULL;
+
+  return r;
+}
+
+/**
+ * @brief Construct an #io_props from its parameters
+ *
+ * @param name Name of the field to read
+ * @param type The type of the data
+ * @param dimension Dataset dimension (1D, 3D, ...)
+ * @param units The units of the dataset
+ * @param gpartSize The size in byte of the particle
+ * @param gparts The particle array
+ * @param functionPtr The function used to convert a g-particle to a double
+ *
+ * Do not call this function directly. Use the macro defined above.
+ */
+INLINE static struct io_props io_make_output_field_convert_gpart_LONGLONG(
+    const char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension,
+    enum unit_conversion_factor units, size_t gpartSize,
+    const struct gpart* gparts, conversion_func_gpart_long_long functionPtr) {
+
+  struct io_props r;
+  strcpy(r.name, name);
+  r.type = type;
+  r.dimension = dimension;
+  r.importance = UNUSED;
+  r.units = units;
+  r.field = NULL;
+  r.partSize = gpartSize;
+  r.parts = NULL;
+  r.xparts = NULL;
+  r.gparts = gparts;
+  r.sparts = NULL;
+  r.conversion = 1;
+  r.convert_part_f = NULL;
+  r.convert_part_d = NULL;
+  r.convert_part_l = NULL;
+  r.convert_gpart_f = NULL;
+  r.convert_gpart_d = NULL;
+  r.convert_gpart_l = functionPtr;
+  r.convert_spart_f = NULL;
+  r.convert_spart_d = NULL;
+  r.convert_spart_l = NULL;
 
   return r;
 }
@@ -424,10 +545,13 @@ INLINE static struct io_props io_make_output_field_convert_spart_FLOAT(
   r.conversion = 1;
   r.convert_part_f = NULL;
   r.convert_part_d = NULL;
+  r.convert_part_l = NULL;
   r.convert_gpart_f = NULL;
   r.convert_gpart_d = NULL;
+  r.convert_gpart_l = NULL;
   r.convert_spart_f = functionPtr;
   r.convert_spart_d = NULL;
+  r.convert_spart_l = NULL;
 
   return r;
 }
@@ -465,10 +589,57 @@ INLINE static struct io_props io_make_output_field_convert_spart_DOUBLE(
   r.conversion = 1;
   r.convert_part_f = NULL;
   r.convert_part_d = NULL;
+  r.convert_part_l = NULL;
   r.convert_gpart_f = NULL;
   r.convert_gpart_d = NULL;
+  r.convert_gpart_l = NULL;
   r.convert_spart_f = NULL;
   r.convert_spart_d = functionPtr;
+  r.convert_spart_l = NULL;
+
+  return r;
+}
+
+/**
+ * @brief Construct an #io_props from its parameters
+ *
+ * @param name Name of the field to read
+ * @param type The type of the data
+ * @param dimension Dataset dimension (1D, 3D, ...)
+ * @param units The units of the dataset
+ * @param spartSize The size in byte of the particle
+ * @param sparts The particle array
+ * @param functionPtr The function used to convert a s-particle to a double
+ *
+ * Do not call this function directly. Use the macro defined above.
+ */
+INLINE static struct io_props io_make_output_field_convert_spart_LONGLONG(
+    const char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension,
+    enum unit_conversion_factor units, size_t spartSize,
+    const struct spart* sparts, conversion_func_spart_long_long functionPtr) {
+
+  struct io_props r;
+  strcpy(r.name, name);
+  r.type = type;
+  r.dimension = dimension;
+  r.importance = UNUSED;
+  r.units = units;
+  r.field = NULL;
+  r.partSize = spartSize;
+  r.parts = NULL;
+  r.xparts = NULL;
+  r.gparts = NULL;
+  r.sparts = sparts;
+  r.conversion = 1;
+  r.convert_part_f = NULL;
+  r.convert_part_d = NULL;
+  r.convert_part_l = NULL;
+  r.convert_gpart_f = NULL;
+  r.convert_gpart_d = NULL;
+  r.convert_gpart_l = NULL;
+  r.convert_spart_f = NULL;
+  r.convert_spart_d = NULL;
+  r.convert_spart_l = functionPtr;
 
   return r;
 }
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 2aa7a9c49642df47cc56149a745fe29f36659426..ec761d3bf9a2540fb21ed8485e9ec8a83f6e204e 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -54,6 +54,7 @@
 #include "stars_io.h"
 #include "tracers_io.h"
 #include "units.h"
+#include "velociraptor_io.h"
 #include "xmf.h"
 
 /* The current limit of ROMIO (the underlying MPI-IO layer) is 2GB */
@@ -957,6 +958,12 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
   const int with_cosmology = e->policy & engine_policy_cosmology;
   const int with_cooling = e->policy & engine_policy_cooling;
   const int with_temperature = e->policy & engine_policy_temperature;
+#ifdef HAVE_VELOCIRAPTOR
+  const int with_stf = (e->policy & engine_policy_structure_finding) &&
+                       (e->s->gpart_group_data != NULL);
+#else
+  const int with_stf = 0;
+#endif
 
   FILE* xmfFile = 0;
   int numFiles = 1;
@@ -1145,15 +1152,25 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
         }
         num_fields += tracers_write_particles(parts, xparts, list + num_fields,
                                               with_cosmology);
-
+        if (with_stf) {
+          num_fields +=
+              velociraptor_write_parts(parts, xparts, list + num_fields);
+        }
         break;
 
       case swift_type_dark_matter:
         darkmatter_write_particles(gparts, list, &num_fields);
+        if (with_stf) {
+          num_fields += velociraptor_write_gparts(e->s->gpart_group_data,
+                                                  list + num_fields);
+        }
         break;
 
       case swift_type_stars:
         stars_write_particles(sparts, list, &num_fields);
+        if (with_stf) {
+          num_fields += velociraptor_write_sparts(sparts, list + num_fields);
+        }
         break;
 
       default:
@@ -1223,6 +1240,12 @@ void write_output_parallel(struct engine* e, const char* baseName,
   const int with_cosmology = e->policy & engine_policy_cosmology;
   const int with_cooling = e->policy & engine_policy_cooling;
   const int with_temperature = e->policy & engine_policy_temperature;
+#ifdef HAVE_VELOCIRAPTOR
+  const int with_stf = (e->policy & engine_policy_structure_finding) &&
+                       (e->s->gpart_group_data != NULL);
+#else
+  const int with_stf = 0;
+#endif
 
   /* Number of particles currently in the arrays */
   const size_t Ntot = e->s->nr_gparts;
@@ -1424,6 +1447,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
     struct part* parts_written = NULL;
     struct xpart* xparts_written = NULL;
     struct gpart* gparts_written = NULL;
+    struct velociraptor_gpart_data* gpart_group_data_written = NULL;
     struct spart* sparts_written = NULL;
 
     /* Write particle fields from the particle structure */
@@ -1440,6 +1464,10 @@ void write_output_parallel(struct engine* e, const char* baseName,
             num_fields += cooling_write_particles(
                 parts, xparts, list + num_fields, e->cooling_func);
           }
+          if (with_stf) {
+            num_fields +=
+                velociraptor_write_parts(parts, xparts, list + num_fields);
+          }
           num_fields += tracers_write_particles(
               parts, xparts, list + num_fields, with_cosmology);
 
@@ -1470,6 +1498,10 @@ void write_output_parallel(struct engine* e, const char* baseName,
                 cooling_write_particles(parts_written, xparts_written,
                                         list + num_fields, e->cooling_func);
           }
+          if (with_stf) {
+            num_fields += velociraptor_write_parts(
+                parts_written, xparts_written, list + num_fields);
+          }
           num_fields += tracers_write_particles(
               parts_written, xparts_written, list + num_fields, with_cosmology);
         }
@@ -1481,6 +1513,10 @@ 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);
+          if (with_stf) {
+            num_fields += velociraptor_write_gparts(e->s->gpart_group_data,
+                                                    list + num_fields);
+          }
         } else {
 
           /* Ok, we need to fish out the particles we want */
@@ -1491,11 +1527,28 @@ void write_output_parallel(struct engine* e, const char* baseName,
                              Ndm_written * sizeof(struct gpart)) != 0)
             error("Error while allocating temporart memory for gparts");
 
+          if (with_stf) {
+            if (posix_memalign(
+                    (void**)&gpart_group_data_written, gpart_align,
+                    Ndm_written * sizeof(struct velociraptor_gpart_data)) != 0)
+              error(
+                  "Error while allocating temporart memory for gparts STF "
+                  "data");
+          }
+
           /* Collect the non-inhibited DM particles from gpart */
-          io_collect_gparts_to_write(gparts, gparts_written, Ntot, Ndm_written);
+          io_collect_gparts_to_write(gparts, e->s->gpart_group_data,
+                                     gparts_written, gpart_group_data_written,
+                                     Ntot, Ndm_written, with_stf);
 
-          /* Write DM particles */
+          /* Select the fields to write */
           darkmatter_write_particles(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;
 
@@ -1505,6 +1558,9 @@ void write_output_parallel(struct engine* e, const char* baseName,
           /* No inhibted particles: easy case */
           Nparticles = Nstars;
           stars_write_particles(sparts, list, &num_fields);
+          if (with_stf) {
+            num_fields += velociraptor_write_sparts(sparts, list + num_fields);
+          }
         } else {
 
           /* Ok, we need to fish out the particles we want */
@@ -1521,6 +1577,10 @@ void write_output_parallel(struct engine* e, const char* baseName,
 
           /* Select the fields to write */
           stars_write_particles(sparts_written, list, &num_fields);
+          if (with_stf) {
+            num_fields +=
+                velociraptor_write_sparts(sparts_written, list + num_fields);
+          }
         }
       } break;
 
@@ -1547,6 +1607,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
     if (parts_written) free(parts_written);
     if (xparts_written) free(xparts_written);
     if (gparts_written) free(gparts_written);
+    if (gpart_group_data_written) free(gpart_group_data_written);
     if (sparts_written) free(sparts_written);
 
 #ifdef IO_SPEED_MEASUREMENT
diff --git a/src/physical_constants.c b/src/physical_constants.c
index 7752f4d3130b7174863d520b3d4d3c6a3e8eb433..3e3c72812c552aba1204086353dc7d239a5c36f9 100644
--- a/src/physical_constants.c
+++ b/src/physical_constants.c
@@ -32,7 +32,8 @@
 /**
  * @brief Converts physical constants to the internal unit system
  *
- * Some constants can be overwritten by the YAML file values.
+ * Some constants can be overwritten by the YAML file values. If the
+ * param argument is NULL, no overwriting is done.
  *
  * @param us The current internal system of units.
  * @param params The parsed parameter file.
@@ -48,8 +49,10 @@ void phys_const_init(const struct unit_system *us, struct swift_params *params,
       const_newton_G_cgs / units_general_cgs_conversion_factor(us, dimension_G);
 
   /* Overwrite G if present in the file */
-  internal_const->const_newton_G = parser_get_opt_param_double(
-      params, "PhysicalConstants:G", internal_const->const_newton_G);
+  if (params != NULL) {
+    internal_const->const_newton_G = parser_get_opt_param_double(
+        params, "PhysicalConstants:G", internal_const->const_newton_G);
+  }
 
   const float dimension_c[5] = {0, 1, -1, 0, 0}; /* [cm s^-1] */
   internal_const->const_speed_light_c =
diff --git a/src/serial_io.c b/src/serial_io.c
index a90e4aceb080b4b53997ce1e8bb66c63fed5ac45..3d9ae4ee0b5a7f4641650f03aac03d5968ee92a2 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -54,6 +54,7 @@
 #include "stars_io.h"
 #include "tracers_io.h"
 #include "units.h"
+#include "velociraptor_io.h"
 #include "xmf.h"
 
 /**
@@ -785,6 +786,12 @@ void write_output_serial(struct engine* e, const char* baseName,
   const int with_cosmology = e->policy & engine_policy_cosmology;
   const int with_cooling = e->policy & engine_policy_cooling;
   const int with_temperature = e->policy & engine_policy_temperature;
+#ifdef HAVE_VELOCIRAPTOR
+  const int with_stf = (e->policy & engine_policy_structure_finding) &&
+                       (e->s->gpart_group_data != NULL);
+#else
+  const int with_stf = 0;
+#endif
 
   FILE* xmfFile = 0;
 
@@ -1090,6 +1097,7 @@ void write_output_serial(struct engine* e, const char* baseName,
         struct part* parts_written = NULL;
         struct xpart* xparts_written = NULL;
         struct gpart* gparts_written = NULL;
+        struct velociraptor_gpart_data* gpart_group_data_written = NULL;
         struct spart* sparts_written = NULL;
 
         /* Write particle fields from the particle structure */
@@ -1106,6 +1114,10 @@ void write_output_serial(struct engine* e, const char* baseName,
                 num_fields += cooling_write_particles(
                     parts, xparts, list + num_fields, e->cooling_func);
               }
+              if (with_stf) {
+                num_fields +=
+                    velociraptor_write_parts(parts, xparts, list + num_fields);
+              }
               num_fields += tracers_write_particles(
                   parts, xparts, list + num_fields, with_cosmology);
 
@@ -1136,7 +1148,10 @@ void write_output_serial(struct engine* e, const char* baseName,
                     cooling_write_particles(parts_written, xparts_written,
                                             list + num_fields, e->cooling_func);
               }
-
+              if (with_stf) {
+                num_fields += velociraptor_write_parts(
+                    parts_written, xparts_written, list + num_fields);
+              }
               num_fields +=
                   tracers_write_particles(parts_written, xparts_written,
                                           list + num_fields, with_cosmology);
@@ -1149,6 +1164,10 @@ 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);
+              if (with_stf) {
+                num_fields += velociraptor_write_gparts(e->s->gpart_group_data,
+                                                        list + num_fields);
+              }
             } else {
 
               /* Ok, we need to fish out the particles we want */
@@ -1159,12 +1178,27 @@ void write_output_serial(struct engine* e, const char* baseName,
                                  Ndm_written * sizeof(struct gpart)) != 0)
                 error("Error while allocating temporart memory for gparts");
 
+              if (with_stf) {
+                if (posix_memalign(
+                        (void**)&gpart_group_data_written, gpart_align,
+                        Ndm_written * sizeof(struct velociraptor_gpart_data)) !=
+                    0)
+                  error(
+                      "Error while allocating temporart memory for gparts STF "
+                      "data");
+              }
+
               /* Collect the non-inhibited DM particles from gpart */
-              io_collect_gparts_to_write(gparts, gparts_written, Ntot,
-                                         Ndm_written);
+              io_collect_gparts_to_write(
+                  gparts, e->s->gpart_group_data, gparts_written,
+                  gpart_group_data_written, Ntot, Ndm_written, with_stf);
 
-              /* Write DM particles */
+              /* Select the fields to write */
               darkmatter_write_particles(gparts_written, list, &num_fields);
+              if (with_stf) {
+                num_fields += velociraptor_write_gparts(
+                    gpart_group_data_written, list + num_fields);
+              }
             }
           } break;
 
@@ -1174,6 +1208,10 @@ void write_output_serial(struct engine* e, const char* baseName,
               /* No inhibted particles: easy case */
               Nparticles = Nstars;
               stars_write_particles(sparts, list, &num_fields);
+              if (with_stf) {
+                num_fields +=
+                    velociraptor_write_sparts(sparts, list + num_fields);
+              }
             } else {
 
               /* Ok, we need to fish out the particles we want */
@@ -1190,6 +1228,10 @@ void write_output_serial(struct engine* e, const char* baseName,
 
               /* Select the fields to write */
               stars_write_particles(sparts_written, list, &num_fields);
+              if (with_stf) {
+                num_fields += velociraptor_write_sparts(sparts_written,
+                                                        list + num_fields);
+              }
             }
           } break;
 
@@ -1216,6 +1258,7 @@ void write_output_serial(struct engine* e, const char* baseName,
         if (parts_written) free(parts_written);
         if (xparts_written) free(xparts_written);
         if (gparts_written) free(gparts_written);
+        if (gpart_group_data_written) free(gpart_group_data_written);
         if (sparts_written) free(sparts_written);
 
         /* Close particle group */
diff --git a/src/single_io.c b/src/single_io.c
index 9115aa0f988113f62851b1d0ea82b5fa2ba314f6..0522174dee3882674239983bc182125e893419c5 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -53,6 +53,7 @@
 #include "stars_io.h"
 #include "tracers_io.h"
 #include "units.h"
+#include "velociraptor_io.h"
 #include "xmf.h"
 
 /**
@@ -648,6 +649,12 @@ void write_output_single(struct engine* e, const char* baseName,
   const int with_cosmology = e->policy & engine_policy_cosmology;
   const int with_cooling = e->policy & engine_policy_cooling;
   const int with_temperature = e->policy & engine_policy_temperature;
+#ifdef HAVE_VELOCIRAPTOR
+  const int with_stf = (e->policy & engine_policy_structure_finding) &&
+                       (e->s->gpart_group_data != NULL);
+#else
+  const int with_stf = 0;
+#endif
 
   /* Number of particles currently in the arrays */
   const size_t Ntot = e->s->nr_gparts;
@@ -893,6 +900,7 @@ void write_output_single(struct engine* e, const char* baseName,
     struct part* parts_written = NULL;
     struct xpart* xparts_written = NULL;
     struct gpart* gparts_written = NULL;
+    struct velociraptor_gpart_data* gpart_group_data_written = NULL;
     struct spart* sparts_written = NULL;
 
     /* Write particle fields from the particle structure */
@@ -909,6 +917,10 @@ void write_output_single(struct engine* e, const char* baseName,
             num_fields += cooling_write_particles(
                 parts, xparts, list + num_fields, e->cooling_func);
           }
+          if (with_stf) {
+            num_fields +=
+                velociraptor_write_parts(parts, xparts, list + num_fields);
+          }
           num_fields += tracers_write_particles(
               parts, xparts, list + num_fields, with_cosmology);
 
@@ -939,6 +951,10 @@ void write_output_single(struct engine* e, const char* baseName,
                 cooling_write_particles(parts_written, xparts_written,
                                         list + num_fields, e->cooling_func);
           }
+          if (with_stf) {
+            num_fields += velociraptor_write_parts(
+                parts_written, xparts_written, list + num_fields);
+          }
           num_fields += tracers_write_particles(
               parts_written, xparts_written, list + num_fields, with_cosmology);
         }
@@ -950,6 +966,10 @@ 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);
+          if (with_stf) {
+            num_fields += velociraptor_write_gparts(e->s->gpart_group_data,
+                                                    list + num_fields);
+          }
         } else {
 
           /* Ok, we need to fish out the particles we want */
@@ -960,11 +980,26 @@ void write_output_single(struct engine* e, const char* baseName,
                              Ndm_written * sizeof(struct gpart)) != 0)
             error("Error while allocating temporart memory for gparts");
 
+          if (with_stf) {
+            if (posix_memalign(
+                    (void**)&gpart_group_data_written, gpart_align,
+                    Ndm_written * sizeof(struct velociraptor_gpart_data)) != 0)
+              error(
+                  "Error while allocating temporart memory for gparts STF "
+                  "data");
+          }
+
           /* Collect the non-inhibited DM particles from gpart */
-          io_collect_gparts_to_write(gparts, gparts_written, Ntot, Ndm_written);
+          io_collect_gparts_to_write(gparts, e->s->gpart_group_data,
+                                     gparts_written, gpart_group_data_written,
+                                     Ntot, Ndm_written, with_stf);
 
-          /* Write DM particles */
+          /* Select the fields to write */
           darkmatter_write_particles(gparts_written, list, &num_fields);
+          if (with_stf) {
+            num_fields += velociraptor_write_gparts(gpart_group_data_written,
+                                                    list + num_fields);
+          }
         }
       } break;
 
@@ -974,6 +1009,9 @@ void write_output_single(struct engine* e, const char* baseName,
           /* No inhibted particles: easy case */
           N = Nstars;
           stars_write_particles(sparts, list, &num_fields);
+          if (with_stf) {
+            num_fields += velociraptor_write_sparts(sparts, list + num_fields);
+          }
         } else {
 
           /* Ok, we need to fish out the particles we want */
@@ -990,6 +1028,10 @@ void write_output_single(struct engine* e, const char* baseName,
 
           /* Select the fields to write */
           stars_write_particles(sparts_written, list, &num_fields);
+          if (with_stf) {
+            num_fields +=
+                velociraptor_write_sparts(sparts_written, list + num_fields);
+          }
         }
       } break;
 
@@ -1015,6 +1057,7 @@ void write_output_single(struct engine* e, const char* baseName,
     if (parts_written) free(parts_written);
     if (xparts_written) free(xparts_written);
     if (gparts_written) free(gparts_written);
+    if (gpart_group_data_written) free(gpart_group_data_written);
     if (sparts_written) free(sparts_written);
 
     /* Close particle group */
diff --git a/src/space.h b/src/space.h
index f3e058e684561bed411af2b7fb63bf2810b7a176..98ab2523668c9789bb644f0ebe300cf73ef6f182 100644
--- a/src/space.h
+++ b/src/space.h
@@ -35,6 +35,7 @@
 #include "lock.h"
 #include "parser.h"
 #include "part.h"
+#include "velociraptor_struct.h"
 
 /* Avoid cyclic inclusions */
 struct cell;
@@ -234,6 +235,9 @@ struct space {
   /*! The associated engine. */
   struct engine *e;
 
+  /*! The group information returned by VELOCIraptor for each #gpart. */
+  struct velociraptor_gpart_data *gpart_group_data;
+
 #ifdef WITH_MPI
 
   /*! Buffers for parts that we will receive from foreign cells. */
diff --git a/src/swift_velociraptor_part.h b/src/swift_velociraptor_part.h
index adae884c2f930c44edf4d48f47f168475bc65885..700842ac5a13e5bee4af15cc0d8726fc668ce421 100644
--- a/src/swift_velociraptor_part.h
+++ b/src/swift_velociraptor_part.h
@@ -21,7 +21,13 @@
 
 #include "part_type.h"
 
-/* SWIFT/VELOCIraptor particle. */
+/**
+ * @brief SWIFT/VELOCIraptor particle.
+ *
+ * This should match the structure Swift::swift_vel_part
+ * defined in the file NBodylib/src/NBody/SwiftParticle.h
+ * of the VELOCIraptor code.
+ */
 struct swift_vel_part {
 
   /*! Particle ID. */
@@ -42,8 +48,18 @@ struct swift_vel_part {
   /*! Internal energy of gas particle */
   float u;
 
+  /*! Temperature of a gas particle */
+  float T;
+
   /*! Type of the #gpart (DM, gas, star, ...) */
   enum part_type type;
+
+  /*! MPI rank on which this #gpart lives on the SWIFT side. */
+  int task;
+
+  /*! Index of this #gpart in the global array of this rank on the SWIFT
+    side. */
+  int index;
 };
 
 #endif /* SWIFT_VELOCIRAPTOR_PART_H */
diff --git a/src/velociraptor_dummy.c b/src/velociraptor_dummy.c
index 15ad1feb19f71571713ee30f8c302e7e83953f55..36cb65bfbe6931464f33d7e4b641f8882fdf65d0 100644
--- a/src/velociraptor_dummy.c
+++ b/src/velociraptor_dummy.c
@@ -20,9 +20,6 @@
 /* Config parameters. */
 #include "../config.h"
 
-/* Some standard headers. */
-#include <stddef.h>
-
 /* Local includes. */
 #include "error.h"
 #include "swift_velociraptor_part.h"
@@ -36,19 +33,41 @@ struct unitinfo {};
 struct cell_loc {};
 struct siminfo {};
 
+/*
 int InitVelociraptor(char *config_name, char *output_name,
                      struct cosmoinfo cosmo_info, struct unitinfo unit_info,
-                     struct siminfo sim_info) {
+                     struct siminfo sim_info, const int numthreads) {
 
   error("This is only a dummy. Call the real one!");
   return 0;
 }
+
 int InvokeVelociraptor(const size_t num_gravity_parts,
                        const size_t num_hydro_parts, const int snapnum,
                        struct swift_vel_part *swift_parts,
-                       const int *cell_node_ids, char *output_name) {
+                       const int *cell_node_ids, char *output_name,
+                       const int numthreads) {
+
+  error("This is only a dummy. Call the real one!");
+  return 0;
+}
+*/
+int InitVelociraptor(char *config_name, struct unitinfo unit_info,
+                     struct siminfo sim_info, const int numthreads) {
 
   error("This is only a dummy. Call the real one!");
   return 0;
 }
+
+struct groupinfo *InvokeVelociraptor(
+    const int snapnum, char *output_name, struct cosmoinfo cosmo_info,
+    struct siminfo sim_info, const size_t num_gravity_parts,
+    const size_t num_hydro_parts, const size_t num_star_parts,
+    struct swift_vel_part *swift_parts, const int *cell_node_ids,
+    const int numthreads, const int return_group_flags,
+    int *const num_in_groups) {
+  error("This is only a dummy. Call the real one!");
+  return 0;
+}
+
 #endif /* HAVE_DUMMY_VELOCIRAPTOR */
diff --git a/src/velociraptor_interface.c b/src/velociraptor_interface.c
index f60a18d86604af2c2e4fd3706a603cbf6ce27199..1049c4730e996112c9b4dc88effad3732af9025d 100644
--- a/src/velociraptor_interface.c
+++ b/src/velociraptor_interface.c
@@ -21,21 +21,23 @@
 #include "../config.h"
 
 /* Some standard headers. */
-#include <errno.h>
 #include <unistd.h>
 
 /* This object's header. */
 #include "velociraptor_interface.h"
 
 /* Local includes. */
-#include "common_io.h"
+#include "cooling.h"
 #include "engine.h"
 #include "hydro.h"
 #include "swift_velociraptor_part.h"
+#include "velociraptor_struct.h"
 
 #ifdef HAVE_VELOCIRAPTOR
 
-/* Structure for passing cosmological information to VELOCIraptor. */
+/**
+ * @brief Structure for passing cosmological information to VELOCIraptor.
+ */
 struct cosmoinfo {
 
   /*! Current expansion factor of the Universe. (cosmology.a) */
@@ -47,6 +49,15 @@ struct cosmoinfo {
   /*! Matter density parameter (cosmology.Omega_m) */
   double Omega_m;
 
+  /*! Radiation density parameter (cosmology.Omega_r) */
+  double Omega_r;
+
+  /*! Neutrino density parameter (0 in SWIFT) */
+  double Omega_nu;
+
+  /*! Neutrino density parameter (cosmology.Omega_k) */
+  double Omega_k;
+
   /*! Baryon density parameter (cosmology.Omega_b) */
   double Omega_b;
 
@@ -60,19 +71,21 @@ struct cosmoinfo {
   double w_de;
 };
 
-/* Structure for passing unit information to VELOCIraptor. */
+/**
+ * @brief Structure for passing unit information to VELOCIraptor.
+ */
 struct unitinfo {
 
-  /* Length conversion factor to kpc. */
+  /*! Length conversion factor to kpc. */
   double lengthtokpc;
 
-  /* Velocity conversion factor to km/s. */
+  /*! Velocity conversion factor to km/s. */
   double velocitytokms;
 
-  /* Mass conversion factor to solar masses. */
+  /*! Mass conversion factor to solar masses. */
   double masstosolarmass;
 
-  /* Potential conversion factor. */
+  /*! Potential conversion factor to (km/s)^2. */
   double energyperunitmass;
 
   /*! Newton's gravitationl constant (phys_const.const_newton_G)*/
@@ -82,18 +95,34 @@ struct unitinfo {
   double hubbleunit;
 };
 
-/* Structure to hold the location of a top-level cell. */
+/**
+ * @brief Structure to hold the location of a top-level cell.
+ */
 struct cell_loc {
 
-  /* Coordinates x,y,z */
+  /*! Coordinates x,y,z */
   double loc[3];
 };
 
-/* Structure for passing simulation information to VELOCIraptor. */
+/**
+ * @brief Structure for passing simulation information to VELOCIraptor for a
+ * given call.
+ */
 struct siminfo {
-  double period, zoomhigresolutionmass, interparticlespacing, spacedimension[3];
 
-  /* Number of top-cells. */
+  /*! Size of periodic replications */
+  double period;
+
+  /*! Mass of the high-resolution DM particles in a zoom-in run. */
+  double zoomhigresolutionmass;
+
+  /*! Mean inter-particle separation of the DM particles */
+  double interparticlespacing;
+
+  /*! Spacial extent of the simulation volume */
+  double spacedimension[3];
+
+  /*! Number of top-level cells. */
   int numcells;
 
   /*! Locations of top-level cells. */
@@ -105,142 +134,135 @@ struct siminfo {
   /*! Inverse of the top-level cell width. */
   double icellwidth[3];
 
+  /*! Holds the node ID of each top-level cell. */
+  int *cellnodeids;
+
+  /*! Is this a cosmological simulation? */
   int icosmologicalsim;
+
+  /*! Is this a zoom-in simulation? */
+  int izoomsim;
+
+  /*! Do we have DM particles? */
+  int idarkmatter;
+
+  /*! Do we have gas particles? */
+  int igas;
+
+  /*! Do we have star particles? */
+  int istar;
+
+  /*! Do we have BH particles? */
+  int ibh;
+
+  /*! Do we have other particles? */
+  int iother;
 };
 
-/* VELOCIraptor interface. */
-int InitVelociraptor(char *config_name, char *output_name,
-                     struct cosmoinfo cosmo_info, struct unitinfo unit_info,
-                     struct siminfo sim_info);
-int InvokeVelociraptor(const size_t num_gravity_parts,
-                       const size_t num_hydro_parts, const int snapnum,
-                       struct swift_vel_part *swift_parts,
-                       const int *cell_node_ids, char *output_name);
+/**
+ * @brief Structure for group information back to swift
+ */
+struct groupinfo {
+
+  /*! Index of a #gpart in the global array on this MPI rank */
+  int index;
+
+  /*! Group number of the #gpart. */
+  long long groupID;
+};
+
+int InitVelociraptor(char *config_name, struct unitinfo unit_info,
+                     struct siminfo sim_info, const int numthreads);
+
+struct groupinfo *InvokeVelociraptor(
+    const int snapnum, char *output_name, struct cosmoinfo cosmo_info,
+    struct siminfo sim_info, const size_t num_gravity_parts,
+    const size_t num_hydro_parts, const size_t num_star_parts,
+    struct swift_vel_part *swift_parts, const int *cell_node_ids,
+    const int numthreads, const int return_group_flags,
+    int *const num_in_groups);
 
 #endif /* HAVE_VELOCIRAPTOR */
 
 /**
- * @brief Initialise VELOCIraptor with input and output file names along with
- * cosmological info needed to run.
+ * @brief Initialise VELOCIraptor with configuration, units,
+ * simulation info needed to run.
  *
  * @param e The #engine.
- *
  */
 void velociraptor_init(struct engine *e) {
 
 #ifdef HAVE_VELOCIRAPTOR
-  struct space *s = e->s;
-  struct cosmoinfo cosmo_info;
-  struct unitinfo unit_info;
-  struct siminfo sim_info;
-
-  /* Set cosmological constants. */
-  cosmo_info.atime = e->cosmology->a;
-  cosmo_info.littleh = e->cosmology->h;
-  cosmo_info.Omega_m = e->cosmology->Omega_m;
-  cosmo_info.Omega_b = e->cosmology->Omega_b;
-  cosmo_info.Omega_Lambda = e->cosmology->Omega_lambda;
-  cosmo_info.Omega_cdm = e->cosmology->Omega_m - e->cosmology->Omega_b;
-  cosmo_info.w_de = e->cosmology->w;
+  const ticks tic = getticks();
 
-  message("Scale factor: %e", cosmo_info.atime);
-  message("Little h: %e", cosmo_info.littleh);
-  message("Omega_m: %e", cosmo_info.Omega_m);
-  message("Omega_b: %e", cosmo_info.Omega_b);
-  message("Omega_Lambda: %e", cosmo_info.Omega_Lambda);
-  message("Omega_cdm: %e", cosmo_info.Omega_cdm);
-  message("w_de: %e", cosmo_info.w_de);
+  /* Internal SWIFT units */
+  const struct unit_system *swift_us = e->internal_units;
 
-  if (e->cosmology->w != -1.)
-    error("w_de is not 1. It is: %lf", e->cosmology->w);
+  /* CGS units and physical constants in CGS */
+  struct unit_system cgs_us;
+  units_init_cgs(&cgs_us);
+  struct phys_const cgs_pc;
+  phys_const_init(&cgs_us, /*params=*/NULL, &cgs_pc);
 
   /* Set unit conversions. */
-  unit_info.lengthtokpc = 1.0;
-  unit_info.velocitytokms = 1.0;
-  unit_info.masstosolarmass = 1.0;
-  unit_info.energyperunitmass = 1.0;
+  struct unitinfo unit_info;
+  unit_info.lengthtokpc =
+      units_cgs_conversion_factor(swift_us, UNIT_CONV_LENGTH) /
+      (1000. * cgs_pc.const_parsec);
+  unit_info.velocitytokms =
+      units_cgs_conversion_factor(swift_us, UNIT_CONV_VELOCITY) / 1.0e5;
+  unit_info.masstosolarmass =
+      units_cgs_conversion_factor(swift_us, UNIT_CONV_MASS) /
+      cgs_pc.const_solar_mass;
+  unit_info.energyperunitmass =
+      units_cgs_conversion_factor(swift_us, UNIT_CONV_ENERGY_PER_UNIT_MASS) /
+      (1.0e10);
   unit_info.gravity = e->physical_constants->const_newton_G;
   unit_info.hubbleunit = e->cosmology->H0 / e->cosmology->h;
 
-  message("Length conversion factor: %e", unit_info.lengthtokpc);
-  message("Velocity conversion factor: %e", unit_info.velocitytokms);
-  message("Mass conversion factor: %e", unit_info.masstosolarmass);
-  message("Potential conversion factor: %e", unit_info.energyperunitmass);
-  message("G: %e", unit_info.gravity);
-  message("H: %e", unit_info.hubbleunit);
-
-  /* TODO: Find the total number of DM particles when running with star
-   * particles and BHs. */
-  const int total_nr_dmparts = e->total_nr_gparts - e->total_nr_parts;
+  /* Gather some information about the simulation */
+  struct siminfo sim_info;
 
-  /* Set simulation information. */
-  if (e->s->periodic) {
-    sim_info.period =
-        unit_info.lengthtokpc *
-        s->dim[0]; /* Physical size of box in VELOCIraptor units (kpc). */
-  } else
-    sim_info.period = 0.0;
-  sim_info.zoomhigresolutionmass = -1.0; /* Placeholder. */
-  sim_info.interparticlespacing = sim_info.period / cbrt(total_nr_dmparts);
-  if (e->policy & engine_policy_cosmology)
+  /* Are we running with cosmology? */
+  if (e->policy & engine_policy_cosmology) {
     sim_info.icosmologicalsim = 1;
-  else
+  } else {
     sim_info.icosmologicalsim = 0;
-  sim_info.spacedimension[0] = unit_info.lengthtokpc * s->dim[0];
-  sim_info.spacedimension[1] = unit_info.lengthtokpc * s->dim[1];
-  sim_info.spacedimension[2] = unit_info.lengthtokpc * s->dim[2];
-  sim_info.numcells = s->nr_cells;
-
-  sim_info.cellwidth[0] = unit_info.lengthtokpc * s->cells_top[0].width[0];
-  sim_info.cellwidth[1] = unit_info.lengthtokpc * s->cells_top[0].width[1];
-  sim_info.cellwidth[2] = unit_info.lengthtokpc * s->cells_top[0].width[2];
-
-  sim_info.icellwidth[0] = s->iwidth[0] / unit_info.lengthtokpc;
-  sim_info.icellwidth[1] = s->iwidth[1] / unit_info.lengthtokpc;
-  sim_info.icellwidth[2] = s->iwidth[2] / unit_info.lengthtokpc;
-
-  /* Only allocate cell location array on first call to velociraptor_init(). */
-  if (e->cell_loc == NULL) {
-    /* Allocate and populate top-level cell locations. */
-    if (posix_memalign((void **)&(e->cell_loc), 32,
-                       s->nr_cells * sizeof(struct cell_loc)) != 0)
-      error("Failed to allocate top-level cell locations for VELOCIraptor.");
-
-    for (int i = 0; i < s->nr_cells; i++) {
-      e->cell_loc[i].loc[0] = unit_info.lengthtokpc * s->cells_top[i].loc[0];
-      e->cell_loc[i].loc[1] = unit_info.lengthtokpc * s->cells_top[i].loc[1];
-      e->cell_loc[i].loc[2] = unit_info.lengthtokpc * s->cells_top[i].loc[2];
-    }
   }
-
-  sim_info.cell_loc = e->cell_loc;
-
-  char configfilename[PARSER_MAX_LINE_SIZE],
-      outputFileName[PARSER_MAX_LINE_SIZE + 128];
-  parser_get_param_string(e->parameter_file,
-                          "StructureFinding:config_file_name", configfilename);
-  snprintf(outputFileName, PARSER_MAX_LINE_SIZE + 128, "%s.VELOCIraptor",
-           e->stfBaseName);
-
-  message("Config file name: %s", configfilename);
-  message("Period: %e", sim_info.period);
-  message("Zoom high res mass: %e", sim_info.zoomhigresolutionmass);
-  message("Inter-particle spacing: %e", sim_info.interparticlespacing);
-  message("Cosmological: %d", sim_info.icosmologicalsim);
-  message("Space dimensions: (%e,%e,%e)", sim_info.spacedimension[0],
-          sim_info.spacedimension[1], sim_info.spacedimension[2]);
-  message("No. of top-level cells: %d", sim_info.numcells);
-  message("Top-level cell locations range: (%e,%e,%e) -> (%e,%e,%e)",
-          sim_info.cell_loc[0].loc[0], sim_info.cell_loc[0].loc[1],
-          sim_info.cell_loc[0].loc[2],
-          sim_info.cell_loc[sim_info.numcells - 1].loc[0],
-          sim_info.cell_loc[sim_info.numcells - 1].loc[1],
-          sim_info.cell_loc[sim_info.numcells - 1].loc[2]);
+  sim_info.izoomsim = 0;
+
+  /* Tell VELOCIraptor what we have in the simulation */
+  sim_info.idarkmatter = (e->total_nr_gparts - e->total_nr_parts > 0);
+  sim_info.igas = (e->policy & engine_policy_hydro);
+  sim_info.istar = (e->policy & engine_policy_stars);
+  sim_info.ibh = 0;  // sim_info.ibh = (e->policy&engine_policy_bh);
+  sim_info.iother = 0;
+
+  /* Be nice, talk! */
+  if (e->verbose) {
+    message("VELOCIraptor conf: Length conversion factor: %e",
+            unit_info.lengthtokpc);
+    message("VELOCIraptor conf: Velocity conversion factor: %e",
+            unit_info.velocitytokms);
+    message("VELOCIraptor conf: Mass conversion factor: %e",
+            unit_info.masstosolarmass);
+    message("VELOCIraptor conf: Internal energy conversion factor: %e",
+            unit_info.energyperunitmass);
+    message("VELOCIraptor conf: G: %e", unit_info.gravity);
+    message("VELOCIraptor conf: H0/h: %e", unit_info.hubbleunit);
+    message("VELOCIraptor conf: Config file name: %s", e->stf_config_file_name);
+    message("VELOCIraptor conf: Cosmological Simulation: %d",
+            sim_info.icosmologicalsim);
+  }
 
   /* Initialise VELOCIraptor. */
-  if (!InitVelociraptor(configfilename, outputFileName, cosmo_info, unit_info,
-                        sim_info))
-    error("Exiting. VELOCIraptor initialisation failed.");
+  if (InitVelociraptor(e->stf_config_file_name, unit_info, sim_info,
+                       e->nr_threads) != 1)
+    error("VELOCIraptor initialisation failed.");
+
+  if (e->verbose)
+    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+            clocks_getunit());
 #else
   error("SWIFT not configure to run with VELOCIraptor.");
 #endif /* HAVE_VELOCIRAPTOR */
@@ -250,121 +272,287 @@ void velociraptor_init(struct engine *e) {
  * @brief Run VELOCIraptor with current particle data.
  *
  * @param e The #engine.
- *
+ * @param linked_with_snap Are we running at the same time as a snapshot dump?
  */
-void velociraptor_invoke(struct engine *e) {
+void velociraptor_invoke(struct engine *e, const int linked_with_snap) {
 
 #ifdef HAVE_VELOCIRAPTOR
-  struct space *s = e->s;
-  struct gpart *gparts = s->gparts;
-  struct part *parts = s->parts;
-  struct xpart *xparts = s->xparts;
+
+  const struct cosmology *cosmo = e->cosmology;
+  const struct hydro_props *hydro_props = e->hydro_properties;
+  const struct unit_system *us = e->internal_units;
+  const struct phys_const *phys_const = e->physical_constants;
+  const struct cooling_function_data *cool_func = e->cooling_func;
+
+  /* Handle on the particles */
+  const struct space *s = e->s;
+  const struct part *parts = s->parts;
+  const struct xpart *xparts = s->xparts;
+  const struct gpart *gparts = s->gparts;
+  const struct spart *sparts = s->sparts;
   const size_t nr_gparts = s->nr_gparts;
-  const size_t nr_hydro_parts = s->nr_parts;
+  const size_t nr_parts = s->nr_parts;
+  const size_t nr_sparts = s->nr_sparts;
   const int nr_cells = s->nr_cells;
-  int *cell_node_ids = NULL;
-  static int stf_output_count = 0;
-  int active_stf_output_count;
+
+  const ticks tic = getticks();
 
   /* Allow thread to run on any core for the duration of the call to
-   * VELOCIraptor so that
-   * when OpenMP threads are spawned they can run on any core on the processor.
-   */
+   * VELOCIraptor so that  when OpenMP threads are spawned
+   * they can run on any core on the processor. */
   const int nr_cores = sysconf(_SC_NPROCESSORS_ONLN);
-  cpu_set_t cpuset;
   pthread_t thread = pthread_self();
 
   /* Set affinity mask to include all cores on the CPU for VELOCIraptor. */
+  cpu_set_t cpuset;
   CPU_ZERO(&cpuset);
   for (int j = 0; j < nr_cores; j++) CPU_SET(j, &cpuset);
-
   pthread_setaffinity_np(thread, sizeof(cpu_set_t), &cpuset);
 
-  ticks tic = getticks();
+  /* Set cosmology information for this point in time */
+  struct cosmoinfo cosmo_info;
+  cosmo_info.atime = e->cosmology->a;
+  cosmo_info.littleh = e->cosmology->h;
+  cosmo_info.Omega_m = e->cosmology->Omega_m;
+  cosmo_info.Omega_b = e->cosmology->Omega_b;
+  cosmo_info.Omega_r = e->cosmology->Omega_r;
+  cosmo_info.Omega_k = e->cosmology->Omega_k;
+  cosmo_info.Omega_nu = 0.;
+  cosmo_info.Omega_Lambda = e->cosmology->Omega_lambda;
+  cosmo_info.Omega_cdm = e->cosmology->Omega_m - e->cosmology->Omega_b;
+  cosmo_info.w_de = e->cosmology->w;
+
+  /* Report the cosmo info we use */
+  if (e->verbose) {
+    message("VELOCIraptor conf: Scale factor: %e", cosmo_info.atime);
+    message("VELOCIraptor conf: Little h: %e", cosmo_info.littleh);
+    message("VELOCIraptor conf: Omega_m: %e", cosmo_info.Omega_m);
+    message("VELOCIraptor conf: Omega_b: %e", cosmo_info.Omega_b);
+    message("VELOCIraptor conf: Omega_Lambda: %e", cosmo_info.Omega_Lambda);
+    message("VELOCIraptor conf: Omega_cdm: %e", cosmo_info.Omega_cdm);
+    message("VELOCIraptor conf: w_de: %e", cosmo_info.w_de);
+  }
+
+  /* Update the simulation information */
+  struct siminfo sim_info;
+
+  /* Period of the box (Note we assume a cubic box!) */
+  if (e->s->periodic) {
+    sim_info.period = s->dim[0];
+  } else {
+    sim_info.period = 0.0;
+  }
+
+  /* Tell VELOCIraptor this is not a zoom-in simulation */
+  sim_info.zoomhigresolutionmass = -1.0;
+
+  /* Are we running with cosmology? */
+  if (e->policy & engine_policy_cosmology) {
+    sim_info.icosmologicalsim = 1;
+    sim_info.izoomsim = 0;
+    const size_t total_nr_baryons = e->total_nr_parts + e->total_nr_sparts;
+    const size_t total_nr_dmparts = e->total_nr_gparts - total_nr_baryons;
+    sim_info.interparticlespacing = sim_info.period / cbrt(total_nr_dmparts);
+  } else {
+    sim_info.icosmologicalsim = 0;
+    sim_info.izoomsim = 0;
+    sim_info.interparticlespacing = -1;
+  }
+
+  /* Set the spatial extent of the simulation volume */
+  sim_info.spacedimension[0] = s->dim[0];
+  sim_info.spacedimension[1] = s->dim[1];
+  sim_info.spacedimension[2] = s->dim[2];
+
+  /* Store number of top-level cells */
+  sim_info.numcells = s->nr_cells;
+
+  /* Size and inverse size of the top-level cells in VELOCIraptor units */
+  sim_info.cellwidth[0] = s->cells_top[0].width[0];
+  sim_info.cellwidth[1] = s->cells_top[0].width[1];
+  sim_info.cellwidth[2] = s->cells_top[0].width[2];
+  sim_info.icellwidth[0] = s->iwidth[0];
+  sim_info.icellwidth[1] = s->iwidth[1];
+  sim_info.icellwidth[2] = s->iwidth[2];
+
+  /* Copy the poisiton of the top-level cells */
+  if (posix_memalign((void **)&sim_info.cell_loc, 32,
+                     s->nr_cells * sizeof(struct cell_loc)) != 0)
+    error("Failed to allocate top-level cell locations for VELOCIraptor.");
+  for (int i = 0; i < s->nr_cells; i++) {
+    sim_info.cell_loc[i].loc[0] = s->cells_top[i].loc[0];
+    sim_info.cell_loc[i].loc[1] = s->cells_top[i].loc[1];
+    sim_info.cell_loc[i].loc[2] = s->cells_top[i].loc[2];
+  }
+
+  if (e->verbose) {
+    message("VELOCIraptor conf: Space dimensions: (%e,%e,%e)",
+            sim_info.spacedimension[0], sim_info.spacedimension[1],
+            sim_info.spacedimension[2]);
+    message("VELOCIraptor conf: No. of top-level cells: %d", sim_info.numcells);
+    message(
+        "VELOCIraptor conf: Top-level cell locations range: (%e,%e,%e) -> "
+        "(%e,%e,%e)",
+        sim_info.cell_loc[0].loc[0], sim_info.cell_loc[0].loc[1],
+        sim_info.cell_loc[0].loc[2],
+        sim_info.cell_loc[sim_info.numcells - 1].loc[0],
+        sim_info.cell_loc[sim_info.numcells - 1].loc[1],
+        sim_info.cell_loc[sim_info.numcells - 1].loc[2]);
+  }
 
   /* Allocate and populate array of cell node IDs. */
+  int *cell_node_ids = NULL;
   if (posix_memalign((void **)&cell_node_ids, 32, nr_cells * sizeof(int)) != 0)
     error("Failed to allocate list of cells node IDs for VELOCIraptor.");
-
   for (int i = 0; i < nr_cells; i++) cell_node_ids[i] = s->cells_top[i].nodeID;
 
-  message("MPI rank %d sending %zu gparts to VELOCIraptor.", engine_rank,
-          nr_gparts);
+  /* Mention the number of particles being sent */
+  if (e->verbose)
+    message(
+        "VELOCIraptor conf: MPI rank %d sending %zu gparts to VELOCIraptor.",
+        engine_rank, nr_gparts);
 
-  /* Append base name with either the step number or time depending on what
-   * format is specified in the parameter file. */
+  /* Append base name with the current output number */
   char outputFileName[PARSER_MAX_LINE_SIZE + 128];
-  if (e->stf_output_freq_format == io_stf_steps)
-    active_stf_output_count = e->step;
-  else if (e->stf_output_freq_format == io_stf_time)
-    active_stf_output_count = stf_output_count;
-  else
-    active_stf_output_count = 0;
 
-  snprintf(outputFileName, PARSER_MAX_LINE_SIZE + 128, "%s_%04i.VELOCIraptor",
-           e->stfBaseName, active_stf_output_count);
+  /* What should the filename be? */
+  if (linked_with_snap) {
+    snprintf(outputFileName, PARSER_MAX_LINE_SIZE + 128,
+             "stf_%s_%04i.VELOCIraptor", e->snapshot_base_name,
+             e->snapshot_output_count);
+  } else {
+    snprintf(outputFileName, PARSER_MAX_LINE_SIZE + 128, "%s_%04i.VELOCIraptor",
+             e->stf_base_name, e->stf_output_count);
+  }
+
+  /* What is the snapshot number? */
+  int snapnum;
+  if (linked_with_snap) {
+    snapnum = e->snapshot_output_count;
+  } else {
+    snapnum = e->stf_output_count;
+  }
 
   /* Allocate and populate an array of swift_vel_parts to be passed to
    * VELOCIraptor. */
   struct swift_vel_part *swift_parts = NULL;
-
   if (posix_memalign((void **)&swift_parts, part_align,
                      nr_gparts * sizeof(struct swift_vel_part)) != 0)
     error("Failed to allocate array of particles for VELOCIraptor.");
 
-  bzero(swift_parts, nr_gparts * sizeof(struct swift_vel_part));
-
-  const float energy_scale = 1.0;
-  const float a = e->cosmology->a;
+  const float a_inv = e->cosmology->a_inv;
 
-  message("Energy scaling factor: %f", energy_scale);
-  message("a: %f", a);
-
-  /* Convert particle properties into VELOCIraptor units */
+  /* Convert particle properties into VELOCIraptor units.
+   * VELOCIraptor wants:
+   * - Co-moving positions,
+   * - Peculiar velocities,
+   * - Co-moving potential,
+   * - Physical internal energy (for the gas),
+   * - Temperatures (for the gas).
+   */
   for (size_t i = 0; i < nr_gparts; i++) {
+
     swift_parts[i].x[0] = gparts[i].x[0];
     swift_parts[i].x[1] = gparts[i].x[1];
     swift_parts[i].x[2] = gparts[i].x[2];
-    swift_parts[i].v[0] = gparts[i].v_full[0] / a;
-    swift_parts[i].v[1] = gparts[i].v_full[1] / a;
-    swift_parts[i].v[2] = gparts[i].v_full[2] / a;
+
+    swift_parts[i].v[0] = gparts[i].v_full[0] * a_inv;
+    swift_parts[i].v[1] = gparts[i].v_full[1] * a_inv;
+    swift_parts[i].v[2] = gparts[i].v_full[2] * a_inv;
+
     swift_parts[i].mass = gravity_get_mass(&gparts[i]);
     swift_parts[i].potential = gravity_get_comoving_potential(&gparts[i]);
+
     swift_parts[i].type = gparts[i].type;
 
+    swift_parts[i].index = i;
+#ifdef WITH_MPI
+    swift_parts[i].task = e->nodeID;
+#else
+    swift_parts[i].task = 0;
+#endif
+
     /* Set gas particle IDs from their hydro counterparts and set internal
      * energies. */
-    if (gparts[i].type == swift_type_gas) {
-      swift_parts[i].id = parts[-gparts[i].id_or_neg_offset].id;
-      swift_parts[i].u =
-          hydro_get_physical_internal_energy(
-              &parts[-gparts[i].id_or_neg_offset],
-              &xparts[-gparts[i].id_or_neg_offset], e->cosmology) *
-          energy_scale;
-    } else if (gparts[i].type == swift_type_dark_matter) {
-      swift_parts[i].id = gparts[i].id_or_neg_offset;
-      swift_parts[i].u = 0.f;
-    } else {
-      error("Particle type not handled by velociraptor (yet?) !");
+    switch (gparts[i].type) {
+
+      case swift_type_gas: {
+        const struct part *p = &parts[-gparts[i].id_or_neg_offset];
+        const struct xpart *xp = &xparts[-gparts[i].id_or_neg_offset];
+
+        swift_parts[i].id = parts[-gparts[i].id_or_neg_offset].id;
+        swift_parts[i].u = hydro_get_drifted_physical_internal_energy(p, cosmo);
+        swift_parts[i].T = cooling_get_temperature(phys_const, hydro_props, us,
+                                                   cosmo, cool_func, p, xp);
+      } break;
+
+      case swift_type_stars:
+
+        swift_parts[i].id = sparts[-gparts[i].id_or_neg_offset].id;
+        swift_parts[i].u = 0.f;
+        swift_parts[i].T = 0.f;
+        break;
+
+      case swift_type_dark_matter:
+
+        swift_parts[i].id = gparts[i].id_or_neg_offset;
+        swift_parts[i].u = 0.f;
+        swift_parts[i].T = 0.f;
+        break;
+
+      default:
+        error("Particle type not handled by VELOCIraptor.");
     }
   }
 
+  /* Values returned by VELOCIRaptor */
+  int num_gparts_in_groups = -1;
+  struct groupinfo *group_info = NULL;
+
   /* Call VELOCIraptor. */
-  if (!InvokeVelociraptor(nr_gparts, nr_hydro_parts, active_stf_output_count,
-                          swift_parts, cell_node_ids, outputFileName))
+  group_info = (struct groupinfo *)InvokeVelociraptor(
+      snapnum, outputFileName, cosmo_info, sim_info, nr_gparts, nr_parts,
+      nr_sparts, swift_parts, cell_node_ids, e->nr_threads, linked_with_snap,
+      &num_gparts_in_groups);
+
+  /* Check that the ouput is valid */
+  if (linked_with_snap && group_info == NULL && num_gparts_in_groups < 0) {
     error("Exiting. Call to VELOCIraptor failed on rank: %d.", e->nodeID);
+  }
+  if (!linked_with_snap && group_info != NULL) {
+    error("VELOCIraptor returned an array whilst it should not have.");
+  }
+
+  /* Assign the group IDs back to the gparts */
+  if (linked_with_snap) {
+
+    if (posix_memalign((void **)&s->gpart_group_data, part_align,
+                       nr_gparts * sizeof(struct velociraptor_gpart_data)) != 0)
+      error("Failed to allocate array of gpart data for VELOCIraptor i/o.");
+
+    struct velociraptor_gpart_data *data = s->gpart_group_data;
+
+    /* Zero the array (gparts not in groups have an ID of 0) */
+    bzero(data, nr_gparts * sizeof(struct velociraptor_gpart_data));
+
+    /* Copy the data at the right place */
+    for (int i = 0; i < num_gparts_in_groups; i++) {
+      data[group_info[i].index].groupID = group_info[i].groupID;
+    }
+
+    /* Free the array returned by VELOCIraptor */
+    free(group_info);
+  }
 
   /* Reset the pthread affinity mask after VELOCIraptor returns. */
   pthread_setaffinity_np(thread, sizeof(cpu_set_t), engine_entry_affinity());
 
-  /* Free cell node ids after VELOCIraptor has copied them. */
-  free(cell_node_ids);
-  free(swift_parts);
-
-  stf_output_count++;
+  /* Increase output counter (if not linked with snapshots) */
+  if (!linked_with_snap) e->stf_output_count++;
 
-  message("VELOCIraptor took %.3f %s on rank %d.",
-          clocks_from_ticks(getticks() - tic), clocks_getunit(), engine_rank);
+  if (e->verbose)
+    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+            clocks_getunit());
 #else
   error("SWIFT not configure to run with VELOCIraptor.");
 #endif /* HAVE_VELOCIRAPTOR */
diff --git a/src/velociraptor_interface.h b/src/velociraptor_interface.h
index 1f29be11c9dd8e267c87201b0a438979fec3775b..2547fa56c1677e93b1c59a1435e9a6ab92c1f308 100644
--- a/src/velociraptor_interface.h
+++ b/src/velociraptor_interface.h
@@ -22,19 +22,11 @@
 /* Config parameters. */
 #include "../config.h"
 
-/**
- * @brief The different formats for when to run structure finding.
- */
-enum io_stf_output_format {
-  io_stf_steps = 0, /*!< Output every N steps */
-  io_stf_time       /*!< Output at fixed time intervals */
-};
-
 /* Forward declaration */
 struct engine;
 
 /* VELOCIraptor wrapper functions. */
 void velociraptor_init(struct engine *e);
-void velociraptor_invoke(struct engine *e);
+void velociraptor_invoke(struct engine *e, const int linked_with_snap);
 
 #endif /* SWIFT_VELOCIRAPTOR_INTERFACE_H */
diff --git a/src/velociraptor_io.h b/src/velociraptor_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..f18398219bfbc5cd6bb58a37b103f29527fa5589
--- /dev/null
+++ b/src/velociraptor_io.h
@@ -0,0 +1,78 @@
+/*******************************************************************************
+ * 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_VELOCIRAPTOR_IO_H
+#define SWIFT_VELOCIRAPTOR_IO_H
+
+/* Config parameters. */
+#include "../config.h"
+
+INLINE static void velociraptor_convert_part_groupID(const struct engine* e,
+                                                     const struct part* p,
+                                                     const struct xpart* xp,
+                                                     long long* ret) {
+  if (p->gpart == NULL)
+    ret[0] = 0.f;
+  else {
+    const ptrdiff_t offset = p->gpart - e->s->gparts;
+    *ret = (e->s->gpart_group_data + offset)->groupID;
+  }
+}
+
+INLINE static void velociraptor_convert_spart_groupID(const struct engine* e,
+                                                      const struct spart* sp,
+                                                      long long* ret) {
+  if (sp->gpart == NULL)
+    ret[0] = 0.f;
+  else {
+    const ptrdiff_t offset = sp->gpart - e->s->gparts;
+    *ret = (e->s->gpart_group_data + offset)->groupID;
+  }
+}
+
+__attribute__((always_inline)) INLINE static int velociraptor_write_parts(
+    const struct part* parts, const struct xpart* xparts,
+    struct io_props* list) {
+
+  list[0] = io_make_output_field_convert_part(
+      "GroupID", LONGLONG, 1, UNIT_CONV_NO_UNITS, parts, xparts,
+      velociraptor_convert_part_groupID);
+
+  return 1;
+}
+
+__attribute__((always_inline)) INLINE static int velociraptor_write_gparts(
+    const struct velociraptor_gpart_data* group_data, struct io_props* list) {
+
+  list[0] = io_make_output_field("GroupID", LONGLONG, 1, UNIT_CONV_NO_UNITS,
+                                 group_data, groupID);
+
+  return 1;
+}
+
+__attribute__((always_inline)) INLINE static int velociraptor_write_sparts(
+    const struct spart* sparts, struct io_props* list) {
+
+  list[0] = io_make_output_field_convert_spart(
+      "GroupID", LONGLONG, 1, UNIT_CONV_NO_UNITS, sparts,
+      velociraptor_convert_spart_groupID);
+
+  return 1;
+}
+
+#endif /* SWIFT_VELOCIRAPTOR_IO_H */
diff --git a/src/velociraptor_struct.h b/src/velociraptor_struct.h
new file mode 100644
index 0000000000000000000000000000000000000000..b998263a6ba2fe0aaa6552f274cb8f4ee85d3b1c
--- /dev/null
+++ b/src/velociraptor_struct.h
@@ -0,0 +1,34 @@
+/*******************************************************************************
+ * 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_VELOCIRAPTOR_STRUCT_H
+#define SWIFT_VELOCIRAPTOR_STRUCT_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/**
+ * @brief Data returned by VELOCIraptor for each #gpart.
+ */
+struct velociraptor_gpart_data {
+
+  /*! Group ID of that #gpart. */
+  long long groupID;
+};
+
+#endif /* SWIFT_VELOCIRAPTOR_STRUCT_H */
diff --git a/tools/analyse_runtime.py b/tools/analyse_runtime.py
index 68dc6776fd89bb90db21cee6826259984f96c9b1..a2c3dd0f201fc47518d6bb0a6a918627db2f3e96 100755
--- a/tools/analyse_runtime.py
+++ b/tools/analyse_runtime.py
@@ -53,98 +53,52 @@ threshold = 0.008
 num_files = len(sys.argv) - 1
 
 labels = [
-    "Gpart assignment",
-    "Mesh comunication",
-    "Forward Fourier transform",
-    "Green function",
-    "Backwards Fourier transform",
-    "engine_recompute_displacement_constraint:",
-    "engine_exchange_top_multipoles:",
-    "updating particle counts",
-    "engine_estimate_nr_tasks:",
-    "Making gravity tasks",
-    "Making hydro tasks",
-    "Splitting tasks",
-    "Counting and linking tasks",
-    "Setting super-pointers",
-    "Making extra hydroloop tasks",
-    "Making extra starsloop tasks",
-    "Linking gravity tasks",
-    "Creating send tasks",
-    "Exchanging cell tags",
-    "Creating recv tasks",
-    "Counting number of foreign particles",
-    "Recursively linking foreign arrays",
-    "Setting unlocks",
-    "Ranking the tasks",
-    "scheduler_reweight:",
-    "space_list_useful_top_level_cells:",
-    "space_rebuild:",
-    "engine_drift_all:",
-    "engine_unskip:",
-    "engine_collect_end_of_step:",
-    "engine_launch:",
-    "writing particle properties",
-    "engine_repartition:",
-    "engine_exchange_cells:",
-    "Dumping restart files",
-    "engine_print_stats:",
-    "engine_marktasks:",
-    "Reading initial conditions",
-    "engine_print_task_counts:",
-    "engine_drift_top_multipoles:",
-    "Communicating rebuild flag",
-    "engine_split:",
-    "space_init",
-    "engine_init",
-    "engine_repartition_trigger:"
-]
-is_rebuild = [
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    1,
-    0,
-    0,
-    0,
-    0,
-    0,
-    0,
-    1,
-    0,
-    0,
-    1,
-    0,
-    0,
-    0,
-    0,
-    0,
-    0,
-    0,
-    0
+    ["Gpart assignment", 1],
+    ["Mesh comunication", 1],
+    ["Forward Fourier transform", 1],
+    ["Green function", 1],
+    ["Backwards Fourier transform", 1],
+    ["engine_recompute_displacement_constraint:", 1],
+    ["engine_exchange_top_multipoles:", 1],
+    ["updating particle counts", 1],
+    ["engine_estimate_nr_tasks:", 1],
+    ["Making gravity tasks", 1],
+    ["Making hydro tasks", 1],
+    ["Splitting tasks", 1],
+    ["Counting and linking tasks", 1],
+    ["Setting super-pointers", 1],
+    ["Making extra hydroloop tasks", 1],
+    ["Making extra starsloop tasks", 1],
+    ["Linking gravity tasks", 1],
+    ["Creating send tasks", 1],
+    ["Exchanging cell tags", 1],
+    ["Creating recv tasks", 1],
+    ["Counting number of foreign particles", 1],
+    ["Recursively linking foreign arrays", 1],
+    ["Setting unlocks", 1],
+    ["Ranking the tasks", 1],
+    ["scheduler_reweight:", 1],
+    ["space_list_useful_top_level_cells:", 1],
+    ["space_rebuild:", 1],
+    ["engine_drift_all:", 0],
+    ["engine_unskip:", 0],
+    ["engine_collect_end_of_step:", 0],
+    ["engine_launch:", 0],
+    ["writing particle properties", 0],
+    ["engine_repartition:", 0],
+    ["engine_exchange_cells:", 1],
+    ["Dumping restart files", 0],
+    ["engine_print_stats:", 0],
+    ["engine_marktasks:", 1],
+    ["Reading initial conditions", 0],
+    ["engine_print_task_counts:", 0],
+    ["engine_drift_top_multipoles:", 0],
+    ["Communicating rebuild flag", 0],
+    ["engine_split:", 0],
+    ["space_init", 0],
+    ["engine_init", 0],
+    ["engine_repartition_trigger:", 0],
+    ["velociraptor_invoke:", 0]
 ]
 times = np.zeros(len(labels))
 counts = np.zeros(len(labels))
@@ -184,20 +138,20 @@ for i in range(num_files):
         for i in range(len(labels)):
 
             # Extract the different blocks
-            if re.search("%s took" % labels[i], line):
+            if re.search("%s took" % labels[i][0], line):
                 counts[i] += 1.0
                 times[i] += float(
                     re.findall(r"[+-]?((\d+\.?\d*)|(\.\d+))", line)[-1][0]
                 )
 
-        # Find the last line with meaningful output (avoid crash report, batch system stuf....)
+        # Find the last line with meaningful output (avoid crash report, batch system stuff....)
         if re.findall(r"\[[0-9]{4}\][ ]\[*", line) or re.findall(
             r"^\[[0-9]*[.][0-9]+\][ ]", line
         ):
             lastline = line
 
     # Total run time
-    total_time += float(re.findall(r"[+-]?([0-9]*[.])?[0-9]+", lastline)[1])
+    total_time += float(re.findall(r"[+-]?(\[[0-9]\])?(\[[0-9]*[.][0-9]*\])+", lastline)[0][1][1:-1])
 
 # Conver to seconds
 times /= 1000.0
@@ -213,35 +167,33 @@ time_ratios = times / total_time
 
 # Better looking labels
 for i in range(len(labels)):
-    labels[i] = labels[i].replace("_", " ")
-    labels[i] = labels[i].replace(":", "")
-    labels[i] = labels[i].title()
+    labels[i][0] = labels[i][0].replace("_", " ")
+    labels[i][0] = labels[i][0].replace(":", "")
+    labels[i][0] = labels[i][0].title()
 
 times = np.array(times)
 time_ratios = np.array(time_ratios)
-is_rebuild = np.array(is_rebuild)
 
 # Sort in order of importance
 order = np.argsort(-times)
 times = times[order]
 counts = counts[order]
 time_ratios = time_ratios[order]
-is_rebuild = is_rebuild[order]
-labels = np.take(labels, order)
+labels = [labels[i] for i in order]
 
 # Keep only the important components
 important_times = [0.0]
 important_ratios = [0.0]
-important_labels = ["Others (all below %.1f\%%)" % (threshold * 100)]
 important_is_rebuild = [0]
+important_labels = ["Others (all below %.1f\%%)" % (threshold * 100)]
 need_print = True
 print("Time spent in the different code sections:")
 for i in range(len(labels)):
     if time_ratios[i] > threshold:
         important_times.append(times[i])
         important_ratios.append(time_ratios[i])
-        important_labels.append(labels[i])
-        important_is_rebuild.append(is_rebuild[i])
+        important_is_rebuild.append(labels[i][1])
+        important_labels.append(labels[i][0])
     else:
         if need_print:
             print("Elements in 'Other' category (<%.1f%%):" % (threshold * 100))
@@ -249,7 +201,7 @@ for i in range(len(labels)):
         important_times[0] += times[i]
         important_ratios[0] += time_ratios[i]
 
-    print(" - '%-40s' (%5d calls, time: %.4fs): %.4f%%" % (labels[i], counts[i], times[i], time_ratios[i] * 100))
+    print(" - '%-40s' (%5d calls, time: %.4fs): %.4f%%" % (labels[i][0], counts[i], times[i], time_ratios[i] * 100))
 
 # Anything unaccounted for?
 print(
@@ -260,8 +212,8 @@ print(
 important_ratios = np.array(important_ratios)
 important_is_rebuild = np.array(important_is_rebuild)
 
-figure()
 
+figure()
 
 def func(pct):
     return "$%4.2f\\%%$" % pct