diff --git a/examples/main.c b/examples/main.c
index 9c6fd9655b172763502e0ac713ee3c5ce3b29d97..23746ba83913b0a7384ef74c90525ea25c6e47ec 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -1201,14 +1201,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 49fa279b64ca2cedf8d0a51ed3fd0a41bb1dcac3..5343c7f6c106f9cde4867acd6098987e22198c28 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -84,6 +84,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)
@@ -160,11 +161,11 @@ DomainDecomposition:
 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          # Scale-factor of the first snaphot (cosmological run)
-  time_first:           0.01          # Time of the first structure finding output (in internal units).
-  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)
+  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 ------------------------------------------
 
diff --git a/src/engine.c b/src/engine.c
index b8e679ba4e1a8602ddef9ca9e32cd93159a10c5b..24370eb103bbd488375c3d8b9dd46a0d645e8799 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3111,13 +3111,6 @@ void engine_step(struct engine *e) {
 #endif
 }
 
-enum output_type {
-  output_none,
-  output_snapshot,
-  output_statistics,
-  output_stf
-};
-
 /**
  * @brief Check whether any kind of i/o has to be performed during this
  * step.
@@ -3132,7 +3125,17 @@ 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);
 
-  /* What kind of output do we want? And at which 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;
 
@@ -3186,7 +3189,22 @@ void engine_check_for_dumps(struct engine *e) {
     switch (type) {
       case output_snapshot:
 
-        /* Dump... */
+        /* Do we want a corresponding VELOCIraptor output? */
+        if (e->snapshot_invoke_stf) {
+
+#ifdef HAVE_VELOCIRAPTOR
+
+          /* Unleash the raptor! */
+          velociraptor_init(e);
+          velociraptor_invoke(e, /*linked_with_snap=*/1);
+#else
+          error(
+              "Asking for a VELOCIraptor output but SWIFT was compiled without "
+              "the interface!");
+#endif
+        }
+
+          /* Dump... */
 #ifdef WITH_LOGGER
         /* Write a file containing the offsets in the particle logger. */
         engine_dump_index(e);
@@ -3212,17 +3230,16 @@ void engine_check_for_dumps(struct engine *e) {
 
 #ifdef HAVE_VELOCIRAPTOR
 
-	/* Unleash the raptor! */
+        /* Unleash the raptor! */
         velociraptor_init(e);
-        velociraptor_invoke(e);
+        velociraptor_invoke(e, /*linked_with_snap=*/1);
 
         /* ... 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!");
+            "the interface!");
 #endif
         break;
 
@@ -3261,7 +3278,8 @@ void engine_check_for_dumps(struct engine *e) {
         }
       }
     }
-  }
+
+  } /* While loop over output types */
 
   /* Restore the information we stored */
   e->ti_current = ti_current;
@@ -4021,9 +4039,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;
@@ -4089,13 +4110,13 @@ 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);
     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->delta_time_stf =
-        parser_get_param_double(params, "StructureFinding:delta_time");
+        parser_get_opt_param_double(params, "StructureFinding:delta_time", -1.);
   }
 
   engine_init_output_lists(e, params);
@@ -4437,14 +4458,16 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
 
     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 {
@@ -4472,7 +4495,10 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
 
     if (e->policy & engine_policy_structure_finding) {
 
-      if (e->delta_time_stf <= 0.)
+      if (e->delta_time_stf == -1. && !e->snapshot_invoke_stf)
+        error("A value for `StructureFinding:delta_time` must be specified");
+
+      if (e->delta_time_stf <= 0. && e->delta_time_stf != -1.)
         error("Time between STF (%e) must be positive.", e->delta_time_stf);
 
       if (e->time_first_stf_output < e->time_begin)
@@ -4510,6 +4536,14 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
     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);
 
diff --git a/src/engine.h b/src/engine.h
index ad8aa7ce33fd01cc18678884067bc8224d71e347..ee7d5aca6892b7c73e152c075a8a980a1f2ea204 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -240,11 +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 */
-  int delta_step_stf;
   double a_first_stf_output;
   double time_first_stf_output;
   double delta_time_stf;
@@ -255,7 +255,8 @@ struct engine {
   /* Integer time of the next stf output */
   integertime_t ti_next_stf;
 
-  char stfBaseName[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/velociraptor_interface.c b/src/velociraptor_interface.c
index ccf1e95352c9e69c39ed454ba50fea652bd319f5..7026ed89805493c9e888e83cdef3665f202d9440 100644
--- a/src/velociraptor_interface.c
+++ b/src/velociraptor_interface.c
@@ -134,6 +134,8 @@ void velociraptor_init(struct engine *e) {
   struct unitinfo unit_info;
   struct siminfo sim_info;
 
+  const ticks tic = getticks();
+
   /* Set cosmological constants. */
   cosmo_info.atime = e->cosmology->a;
   cosmo_info.littleh = e->cosmology->h;
@@ -220,7 +222,7 @@ void velociraptor_init(struct engine *e) {
   parser_get_param_string(e->parameter_file,
                           "StructureFinding:config_file_name", configfilename);
   snprintf(outputFileName, PARSER_MAX_LINE_SIZE + 128, "%s.VELOCIraptor",
-           e->stfBaseName);
+           e->stf_base_name);
 
   message("Config file name: %s", configfilename);
   message("Period: %e", sim_info.period);
@@ -241,6 +243,10 @@ void velociraptor_init(struct engine *e) {
   if (!InitVelociraptor(configfilename, outputFileName, cosmo_info, unit_info,
                         sim_info))
     error("Exiting. 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,9 +256,9 @@ 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;
@@ -263,7 +269,8 @@ void velociraptor_invoke(struct engine *e) {
   const size_t nr_hydro_parts = s->nr_parts;
   const int nr_cells = s->nr_cells;
   int *cell_node_ids = NULL;
-  static int stf_output_count = 0;
+
+  const ticks tic = getticks();
 
   /* Allow thread to run on any core for the duration of the call to
    * VELOCIraptor so that
@@ -279,23 +286,22 @@ void velociraptor_invoke(struct engine *e) {
 
   pthread_setaffinity_np(thread, sizeof(cpu_set_t), &cpuset);
 
-  ticks tic = getticks();
-
   /* Allocate and populate array of cell node IDs. */
   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);
+  if (e->verbose)
+    message("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. */
   char outputFileName[PARSER_MAX_LINE_SIZE + 128];
 
   snprintf(outputFileName, PARSER_MAX_LINE_SIZE + 128, "%s_%04i.VELOCIraptor",
-           e->stfBaseName, stf_output_count);
+           e->stf_base_name, e->stf_output_count);
 
   /* Allocate and populate an array of swift_vel_parts to be passed to
    * VELOCIraptor. */
@@ -310,8 +316,8 @@ void velociraptor_invoke(struct engine *e) {
   const float energy_scale = 1.0;
   const float a = e->cosmology->a;
 
-  message("Energy scaling factor: %f", energy_scale);
-  message("a: %f", a);
+  /* message("Energy scaling factor: %f", energy_scale); */
+  /* message("a: %f", a); */
 
   /* Convert particle properties into VELOCIraptor units */
   for (size_t i = 0; i < nr_gparts; i++) {
@@ -343,7 +349,7 @@ void velociraptor_invoke(struct engine *e) {
   }
 
   /* Call VELOCIraptor. */
-  if (!InvokeVelociraptor(nr_gparts, nr_hydro_parts, stf_output_count,
+  if (!InvokeVelociraptor(nr_gparts, nr_hydro_parts, e->stf_output_count,
                           swift_parts, cell_node_ids, outputFileName))
     error("Exiting. Call to VELOCIraptor failed on rank: %d.", e->nodeID);
 
@@ -354,10 +360,11 @@ void velociraptor_invoke(struct engine *e) {
   free(cell_node_ids);
   free(swift_parts);
 
-  stf_output_count++;
+  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 67ca9153d3301dedc1f7529eedce8e63f8b392c0..2547fa56c1677e93b1c59a1435e9a6ab92c1f308 100644
--- a/src/velociraptor_interface.h
+++ b/src/velociraptor_interface.h
@@ -27,6 +27,6 @@ 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 */