diff --git a/examples/SmallCosmoVolume/stf_input_6dfof_dmonly_sub.cfg b/examples/SmallCosmoVolume/stf_input_6dfof_dmonly_sub.cfg
index e2859da17f0a71e26394d0ba40ed744d8f6ad124..872e0ad6f44d8092ce1da6ac030a949dc4dba5d5 100644
--- a/examples/SmallCosmoVolume/stf_input_6dfof_dmonly_sub.cfg
+++ b/examples/SmallCosmoVolume/stf_input_6dfof_dmonly_sub.cfg
@@ -7,19 +7,9 @@
 ################################
 #input is from a cosmological so can use parameters like box size, h, Omega_m to calculate length and density scales
 Cosmological_input=1
-#sets the size of the chunk input particle is read in number of particles (ie: memory allocated to store all useful/desired input data)
-#Typically not need to be set, default value is
-#Input_chunk_size=100000
-#sets the total buffer size in bytes used to store temporary particle information
-#of mpi read threads before they are broadcast to the appropriate waiting non-read threads
-#if not set, default value of -1 is equivalent to 1e6 particles per mpi process, quite large
-#but significantly minimises the number of send/receives
-#MPI_particle_total_buf_size=-1
-
-#gadget input related
-#NSPH_extra_blocks=0 #read extra sph blocks
-#NStar_extra_blocks=0 #read extra star blocks
-#NBH_extra_blocks=0 #read extra black hole blocks
+
+#Type of snapshot to read. Ignored when using within SWIFT.
+HDF_name_convention=6 # ILLUSTRIS 0, GADGETX 1, EAGLE 2, GIZMO 3, SIMBA 4, MUFASA 5, SWIFTEAGLE 6
 
 Particle_search_type=2 #search all particles, see allvars for other types
 Baryon_searchflag=0 #if 1 search for baryons separately using phase-space search when identifying substructures, 2 allows special treatment in field FOF linking and phase-space substructure search, 0 treat the same as dark matter particles
@@ -30,35 +20,33 @@ Halo_core_search=1
 Significance_level=1.0 #how significant a substructure is relative to Poisson noise. Values >= 1 are fine.
 
 ################################
-#unit options, should always be provided
+# unit options, should always be provided
 ################################
-#conversion of output length units to kpc
-Length_unit_to_kpc=1.0
-#conversion of output velocity units to km/s
-Velocity_to_kms=1.0
-#conversion of output mass units to solar masses
-Mass_to_solarmass=1.0
-#units conversion from input input to desired internal unit
-Length_unit=1.0 #default code unit,
-Velocity_unit=1.0 #default velocity unit,
-Mass_unit=1.0 #default mass unit,
-Gravity=4.302051e+01 #for 1e10 Msun, km/s and Mpc
-Hubble_unit=100.0 # assuming units are km/s and Mpc, then value of Hubble in km/s/Mpc
 
-################################
-#compute related options
-################################
-#if this is set, the code will store (or search for) local velocity density distribution data in this file.
-#the most time consuming part of the code is the calculation of the local density distribution, therefore if
-#one wishes to repeat searches, saving and reading this file can significantly reduce computation time
-#Output_den=1
+# This is only for i/o. Specifies what units the code was running in.
+# These should be set to whatever internal units we use.
+# They have no impact on the way the code runs.
+Length_unit_to_kpc=1000. #conversion of output length units to kpc
+Velocity_to_kms=1.0      #conversion of output velocity units to km/s
+Mass_to_solarmass=1e+10  #conversion of output mass units to solar masses
+
+# units conversion from input to desired internal unit.
+# These should be set to 1 unless a conversion is expected.
+Length_unit=1.0    #default length unit
+Velocity_unit=1.0  #default velocity unit
+Mass_unit=1.0      #default mass unit
+
+# These are ignored when running within SWIFT.
+# When using standalone code, G and H must match the value used in the run. 
+Gravity=4.300927e+01   # In internal units (here 10^10 Msun, km/s, Mpc)
+Hubble_unit=100.0      # This is H0 / h in internal units. 
 
 ################################
 #search related options
 ################################
 
 #how to search a simulation
- # searches for separate 6dfof cores in field haloes, and then more than just flags halo as merging, assigns particles to each merging "halo". 2 is full separation, 1 is flagging, 0 is off
+# searches for separate 6dfof cores in field haloes, and then more than just flags halo as merging, assigns particles to each merging "halo". 2 is full separation, 1 is flagging, 0 is off
 #also useful for zoom simulations or simulations of individual objects, setting this flag means no field structure search is run
 Singlehalo_search=0 #if file is single halo in which one wishes to search for substructure
 #additional option for field haloes
@@ -69,7 +57,6 @@ Minimum_size=20 #min 20 particles
 Minimum_halo_size=-1 #if field halos have different minimum sizes, otherwise set to -1.
 
 #for field fof halo search
-
 Halo_linking_length_factor=2.0 #factor by which Physical_linking_length is changed when searching for field halos. Typical values are ~2 when using iterative substructure search.
 Halo_velocity_linking_length_factor=5.0 #for 6d fof halo search increase ellv from substructure search
 
@@ -88,7 +75,6 @@ Velocity_opening_angle=0.10 #angle between velocities. 18 degrees here, typical
 Physical_linking_length=0.10 #physical linking length. IF reading periodic volumes in gadget/hdf/ramses, in units of the effective inter-particle spacing. Otherwise in user defined code units. Here set to 0.10 as iterative flag one, values of 0.1-0.3 are typical.
 Velocity_linking_length=0.20 #where scaled by structure dispersion
 
-
 #for iterative substructure search, rarely ever need to change this
 Iterative_threshold_factor=1.0 #change in threshold value when using iterative search. Here no increase in threshold if iterative or not
 Iterative_linking_length_factor=2.0 #increase in final linking final iterative substructure search will be sqrt(2.25)*this factor
@@ -106,9 +92,6 @@ Halo_core_adaptive_sigma_fac=2.0 #used when running fully adaptive core search w
 Halo_core_num_loops=3 #allows the core search to iterate, shrinking the velocity linking length to used till the number of cores identified decreases or this limit is reached. Allows apative search with larger linking length to be robust. Typically values are 3-5
 Halo_core_loop_ellv_fac=0.75 #Factor by which velocity linking length is decreased when running loops for core search.  Typically values are 0.75
 
-#for zoom simulations, alter the effective resolution, allowing quick scaling of linking lenghts passed
-#Effective_resolution=1024.0 #here effective resolution of 1024^3
-
 ################################
 #Unbinding options (VELOCIraptor is able to accurately identify tidal debris so particles need not be bound to a structure)
 ################################
@@ -125,17 +108,6 @@ Softening_length=0.
 #don't keep background potential when unbinding
 Keep_background_potential=0
 
-################################
-#Cosmological parameters
-#this is typically overwritten by information in the gadget/hdf header if those input file types are read
-################################
-h_val=1.0
-Omega_m=0.3
-Omega_Lambda=0.7
-Critical_density=1.0
-Virial_density=200 #so-called virial overdensity value
-Omega_b=0. #no baryons
-
 ################################
 #Calculation of properties related options
 ################################
@@ -151,7 +123,6 @@ Comoving_units=0
 Write_group_array_file=0 #write a group array file
 Separate_output_files=0 #separate output into field and substructure files similar to subfind
 Binary_output=2 #binary output 1, ascii 0, and HDF 2
-HDF_name_convention=6
 
 #halo ids are adjusted by this value * 1000000000000 (or 1000000 if code compiled with the LONGINTS option turned off)
 #to ensure that halo ids are temporally unique. So if you had 100 snapshots, for snap 100 set this to 100 and 100*1000000000000 will
diff --git a/examples/main.c b/examples/main.c
index e8577a353c8a4de48e26064c7dbd324feffe6119..773accc461fa26c819f0b603e7de689fa50e6cc6 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -943,8 +943,13 @@ int main(int argc, char *argv[]) {
     engine_dump_snapshot(&e);
     engine_print_stats(&e);
 
+#ifdef HAVE_VELOCIRAPTOR
     /* Call VELOCIraptor for the first time after the first snapshot dump. */
-    // if (e.policy & engine_policy_structure_finding) velociraptor_invoke(&e);
+    // if (e.policy & engine_policy_structure_finding) {
+    // velociraptor_init(&e);
+    // velociraptor_invoke(&e);
+    //}
+#endif
   }
 
   /* Legend */
diff --git a/src/engine.c b/src/engine.c
index 89af0b0966dd3683782eb219a1123f0a315228d4..1808de3b5c422480f7a63b1ad771160332d8abf5 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -4983,13 +4983,16 @@ void engine_step(struct engine *e) {
   /* Perform structure finding? */
   if (run_stf) {
 
-    // MATTHIEU: Add a drift_all here. And check the order with the order i/o
-    // options.
+  // MATTHIEU: Add a drift_all here. And check the order with the order i/o
+  // options.
 
+#ifdef HAVE_VELOCIRAPTOR
+    velociraptor_init(e);
     velociraptor_invoke(e);
 
     /* ... and find the next output time */
     if (e->stf_output_freq_format == TIME) engine_compute_next_stf_time(e);
+#endif
   }
 
   /* Restore the information we stored */
@@ -5799,6 +5802,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
   e->chemistry = chemistry;
   e->sourceterms = sourceterms;
   e->parameter_file = params;
+  e->cell_loc = NULL;
 #ifdef WITH_MPI
   e->cputime_last_step = 0;
   e->last_repartition = 0;
@@ -5892,7 +5896,6 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
         parser_get_param_double(params, "StructureFinding:time_first");
     e->a_first_stf = parser_get_opt_param_double(
         params, "StructureFinding:scale_factor_first", 0.1);
-    // velociraptor_init(e);
     e->stf_output_freq_format =
         parser_get_param_int(params, "StructureFinding:output_time_format");
     if (e->stf_output_freq_format == STEPS) {
diff --git a/src/velociraptor_interface.c b/src/velociraptor_interface.c
index 3d9dab8ee68dd85eb54e4df3e05622854aef1aa0..d7331ce49f102f52adafff1364dce173fc247586 100644
--- a/src/velociraptor_interface.c
+++ b/src/velociraptor_interface.c
@@ -30,6 +30,7 @@
 /* Local includes. */
 #include "common_io.h"
 #include "engine.h"
+#include "hydro.h"
 #include "swift_velociraptor_part.h"
 
 #ifdef HAVE_VELOCIRAPTOR
@@ -77,13 +78,16 @@ void velociraptor_init(struct engine *e) {
   message("Omega_cdm: %e", cosmo_info.Omega_cdm);
   message("w_de: %e", cosmo_info.w_de);
 
+  if (e->cosmology->w != -1.)
+    error("w_de is not 1. It is: %lf", e->cosmology->w);
+
   /* Set unit conversions. */
   unit_info.lengthtokpc = 1.0;
   unit_info.velocitytokms = 1.0;
   unit_info.masstosolarmass = 1.0;
   unit_info.energyperunitmass = 1.0;
   unit_info.gravity = e->physical_constants->const_newton_G;
-  unit_info.hubbleunit = e->cosmology->H; /* TODO: double check this. */
+  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);
@@ -122,15 +126,18 @@ void velociraptor_init(struct engine *e) {
   sim_info.icellwidth[1] = s->iwidth[1] / unit_info.lengthtokpc;
   sim_info.icellwidth[2] = s->iwidth[2] / unit_info.lengthtokpc;
 
-  /* 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];
+  /* 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;
@@ -230,23 +237,21 @@ void velociraptor_invoke(struct engine *e) {
   bzero(swift_parts, nr_gparts * sizeof(struct swift_vel_part));
 
   const float energy_scale = 1.0;
-  const float a2 = e->cosmology->a * e->cosmology->a;
+  const float a = e->cosmology->a;
 
   message("Energy scaling factor: %f", energy_scale);
-  message("a^2: %f", a2);
+  message("a: %f", a);
 
   /* Convert particle properties into VELOCIraptor units */
   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] / a2;  // MATTHIEU: Check this a^2 !!
-    swift_parts[i].v[1] = gparts[i].v_full[1] / a2;
-    swift_parts[i].v[2] = gparts[i].v_full[2] / a2;
+    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].mass = gravity_get_mass(&gparts[i]);
-    swift_parts[i].potential = gravity_get_comoving_potential(
-        &gparts[i]);  // MATTHIEU: Need factors here?
+    swift_parts[i].potential = gravity_get_comoving_potential(&gparts[i]);
     swift_parts[i].type = gparts[i].type;
 
     /* Set gas particle IDs from their hydro counterparts and set internal