diff --git a/.gitignore b/.gitignore
index 51976c1ddf004c12570565dac227c022bda92efe..27ff2cf7ff880de2007e27978fcc1cfa22d2bd75 100644
--- a/.gitignore
+++ b/.gitignore
@@ -42,6 +42,11 @@ examples/*/*/*.hdf5
 examples/*/snapshots*
 examples/*/restart/*
 examples/*/*/used_parameters.yml
+examples/*/err_file*
+examples/*/out_file*
+examples/*/stf_output*
+examples/*/stf_ouput*
+examples/*/log*
 examples/*/*/unused_parameters.yml
 examples/*/*.mpg
 examples/*/gravity_checks_*.dat
diff --git a/configure.ac b/configure.ac
index 32ef711f2e306cf16edd6ac5f58a08bd2c172dc1..56f88116c845def33300b87630524f06fec666bb 100644
--- a/configure.ac
+++ b/configure.ac
@@ -771,6 +771,38 @@ if test "$with_hdf5" = "yes"; then
 fi
 AM_CONDITIONAL([HAVEPARALLELHDF5],[test "$have_parallel_hdf5" = "yes"])
 
+# Check for VELOCIraptor.
+have_velociraptor="no"
+AC_ARG_WITH([velociraptor],
+    [AS_HELP_STRING([--with-velociraptor=PATH],
+       [Directory where velociraptor library exists @<:@yes/no@:>@]
+    )],
+    [with_velociraptor="$withval"],
+    [with_velociraptor="no"]
+)
+if test "x$with_velociraptor" != "xno"; then
+   AC_PROG_FC
+   AC_FC_LIBRARY_LDFLAGS
+   if test "x$with_velociraptor" != "xyes" -a "x$with_velociraptor" != "x"; then
+      VELOCIRAPTOR_LIBS="-L$with_velociraptor -lstf -lstdc++ -lhdf5_cpp"
+      CFLAGS="$CFLAGS -fopenmp"
+   else
+      VELOCIRAPTOR_LIBS=""
+   fi
+
+   have_velociraptor="yes"
+
+   AC_CHECK_LIB(
+      [stf],
+      [InitVelociraptor],
+      [AC_DEFINE([HAVE_VELOCIRAPTOR],1,[The VELOCIraptor library appears to be present.])],
+      [AC_MSG_ERROR(Cannot find VELOCIraptor library at $with_velociraptor)],
+      [$VELOCIRAPTOR_LIBS $HDF5_LDFLAGS $HDF5_LIBS $GSL_LIBS]
+   )
+fi
+AC_SUBST([VELOCIRAPTOR_LIBS])
+AM_CONDITIONAL([HAVEVELOCIRAPTOR],[test -n "$VELOCIRAPTOR_LIBS"])
+
 # Check for floating-point execeptions
 AC_CHECK_FUNC(feenableexcept, AC_DEFINE([HAVE_FE_ENABLE_EXCEPT],[1],
     [Defined if the floating-point exception can be enabled using non-standard GNU functions.]))
@@ -1348,21 +1380,22 @@ AC_MSG_RESULT([
 
    $PACKAGE_NAME v.$PACKAGE_VERSION
 
-   Compiler           : $CC
-    - vendor          : $ax_cv_c_compiler_vendor
-    - version         : $ax_cv_c_compiler_version
-    - flags           : $CFLAGS
-   MPI enabled        : $enable_mpi
-   HDF5 enabled       : $with_hdf5
-    - parallel        : $have_parallel_hdf5
-   Metis enabled      : $have_metis
-   FFTW3 enabled      : $have_fftw
-   GSL enabled        : $have_gsl
-   libNUMA enabled    : $have_numa
-   GRACKLE enabled    : $have_grackle
-   Special allocators : $have_special_allocator
-   CPU profiler       : $have_profiler
-   Pthread barriers   : $have_pthread_barrier
+   Compiler             : $CC
+    - vendor            : $ax_cv_c_compiler_vendor
+    - version           : $ax_cv_c_compiler_version
+    - flags             : $CFLAGS
+   MPI enabled          : $enable_mpi
+   HDF5 enabled         : $with_hdf5
+    - parallel          : $have_parallel_hdf5
+   Metis enabled        : $have_metis
+   FFTW3 enabled        : $have_fftw
+   GSL enabled          : $have_gsl
+   libNUMA enabled      : $have_numa
+   GRACKLE enabled      : $have_grackle
+   Special allocators   : $have_special_allocator
+   CPU profiler         : $have_profiler
+   Pthread barriers     : $have_pthread_barrier
+   VELOCIraptor enabled : $have_velociraptor
 
    Hydro scheme       : $with_hydro
    Dimensionality     : $with_dimension
diff --git a/examples/EAGLE_25/eagle_25.yml b/examples/EAGLE_25/eagle_25.yml
index 904f8c58168e0bae0ca36c5c1a71f1d006d4e3f2..d6f9ad2474cb4fc207145c73a1c1c694f2f11386 100644
--- a/examples/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_25/eagle_25.yml
@@ -6,6 +6,16 @@ InternalUnitSystem:
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
 
+# Structure finding options
+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
 Cosmology:
   h:              0.6777        # Reduced Hubble constant
diff --git a/examples/EAGLE_6/eagle_6.yml b/examples/EAGLE_6/eagle_6.yml
index 03bb82b48394525c12b7a8df0e5c5df245fa2669..eb374df964e8b021ef2b7d90caf8a1824cf3a833 100644
--- a/examples/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_6/eagle_6.yml
@@ -6,6 +6,16 @@ InternalUnitSystem:
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
 
+# Structure finding options
+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
 Cosmology:
   h:              0.6777        # Reduced Hubble constant
diff --git a/examples/EAGLE_6/testVELOCIraptor.sh b/examples/EAGLE_6/testVELOCIraptor.sh
new file mode 100755
index 0000000000000000000000000000000000000000..14ec30487006f0b7e86356837c9a801950c15c83
--- /dev/null
+++ b/examples/EAGLE_6/testVELOCIraptor.sh
@@ -0,0 +1,109 @@
+#!/bin/bash
+
+VELOCIRAPTOR_PATH=/cosma7/data/dp004/dc-will2/VELOCIraptor-STF/stf
+SCRIPT_PATH=$VELOCIRAPTOR_PATH/examples
+TREEFROG_PATH=$VELOCIRAPTOR_PATH/bin
+
+# List each test that should be run
+declare -a DM_TEST_LIST=(6dfof_dmonly_sub)
+declare -a GAS_TEST_LIST=(6dfof_gas_sub)
+INFO_FILE_TEXT="$TREEFROG_PATH/treefrog \n-s 2 -d 2 -F 0 -B 2 -m -1 -N 1 -T -1\n2 # No. of outputs to compare"
+NUM_MPI_PROC=1
+
+# Check for command line arguments
+while getopts 'g:s:' opt ; do
+  case $opt in
+    g) RUN_DM=$OPTARG ;;
+    s) RUN_GAS=$OPTARG ;;
+  esac
+done
+
+# Run catalog comparison for DM only runs
+if [ "$RUN_DM" = "1" ]; then
+  for TEST in "${DM_TEST_LIST[@]}"
+  do
+
+    # Create output directory
+    OUTPUT=halo_$NUM_MPI_PROC"_mpi_1_"$TEST
+    VEL_OUTPUT=vel_output_$NUM_MPI_PROC"_mpi_1_"$TEST
+
+    if [ ! -d $OUTPUT ]; then mkdir $OUTPUT; fi
+    if [ ! -d $VEL_OUTPUT ]; then mkdir $VEL_OUTPUT; fi
+
+    # Remove old comparison files
+    rm catcomp*
+    rm $OUTPUT/stf* 
+    rm $VEL_OUTPUT/vel_$TEST*
+
+    # Run test using SWIFT + VELOCIraptor
+    echo "Running: mpirun -np $NUM_MPI_PROC ../swift_mpi -G -t 8 eagle_6.yml -x -n 5 -P StructureFinding:basename:./$OUTPUT/stf -P StructureFinding:config_file_name:./stf_input_$TEST.cfg -P Snapshots:basename:./eagle_dmonly"
+    mpirun -np $NUM_MPI_PROC ../swift_mpi -G -t 8 eagle_6.yml -x -n 5 -P StructureFinding:basename:./$OUTPUT/stf -P StructureFinding:config_file_name:./stf_input_$TEST.cfg -P Snapshots:basename:./eagle_dmonly
+
+    # Run test using VELOCIraptor
+    echo "Running: mpirun -np $NUM_MPI_PROC $VELOCIRAPTOR_PATH/bin/stf-gas -I 2 -i eagle_dmonly_0000 -C $VELOCIRAPTOR_PATH/stf_input_$TEST.cfg -o ./$VEL_OUTPUT/vel_$TEST"
+    mpirun -np $NUM_MPI_PROC $VELOCIRAPTOR_PATH/bin/stf-gas -I 2 -i eagle_dmonly_0000 -C ./stf_input_$TEST.cfg -o ./$VEL_OUTPUT/vel_$TEST
+
+    # Create info file for python comparison script
+    echo -e $INFO_FILE_TEXT > infoFile_$TEST.txt
+    echo "vel_$TEST $VEL_OUTPUT/vel_$TEST" >> infoFile_$TEST.txt
+    echo "stf_$TEST $OUTPUT/stf_0000.VELOCIraptor" >> infoFile_$TEST.txt
+
+    # Run comparison script on VELOCIraptor output
+    if python $SCRIPT_PATH/catalogcomparisontolerancecheck.py infoFile_$TEST.txt toIfile.txt
+    then
+      echo "Catalog comparison passed: "$TEST
+    else
+      echo "Catalog comparison failed: "$TEST
+      exit 1
+    fi
+
+    echo "------------"
+
+  done
+fi
+
+# Run catalog comparison for DM + GAS runs
+if [ "$RUN_GAS" = "1" ]; then
+  for TEST in "${GAS_TEST_LIST[@]}"
+  do
+
+    # Create output directory
+    OUTPUT=halo_$NUM_MPI_PROC"_mpi_1_"$TEST
+    VEL_OUTPUT=vel_output_$NUM_MPI_PROC"_mpi_1_"$TEST
+
+    if [ ! -d $OUTPUT ]; then mkdir $OUTPUT; fi
+    if [ ! -d $VEL_OUTPUT ]; then mkdir $VEL_OUTPUT; fi
+
+    # Remove old comparison files
+    rm catcomp*
+    rm $OUTPUT/stf* 
+    rm $VEL_OUTPUT/vel_$TEST*
+
+    # Run test using SWIFT + VELOCIraptor
+    echo "Running: mpirun -np $NUM_MPI_PROC ../swift_mpi -s -G -t 8 eagle_6.yml -x -n 5 -P StructureFinding:basename:./$OUTPUT/stf -P StructureFinding:config_file_name:./stf_input_$TEST.cfg -P Snapshots:basename:./eagle_gas"
+    mpirun -np $NUM_MPI_PROC ../swift_mpi -s -G -t 8 eagle_6.yml -x -n 5 -P StructureFinding:basename:./$OUTPUT/stf -P StructureFinding:config_file_name:./stf_input_$TEST.cfg -P Snapshots:basename:./eagle_gas
+
+    # Run test using VELOCIraptor
+    echo "Running: mpirun -np $NUM_MPI_PROC $VELOCIRAPTOR_PATH/bin/stf-gas -I 2 -i eagle_gas_0000 -C ./stf_input_$TEST.cfg -o ./$VEL_OUTPUT/vel_$TEST"
+    mpirun -np $NUM_MPI_PROC $VELOCIRAPTOR_PATH/bin/stf-gas -I 2 -i eagle_gas_0000 -C ./stf_input_$TEST.cfg -o ./$VEL_OUTPUT/vel_$TEST
+
+    # Create info file for python comparison script
+    echo -e $INFO_FILE_TEXT > infoFile_$TEST.txt
+    echo "vel_$TEST $VEL_OUTPUT/vel_$TEST" >> infoFile_$TEST.txt
+    echo "stf_$TEST $OUTPUT/stf_0000.VELOCIraptor" >> infoFile_$TEST.txt
+
+    # Run comparison script on VELOCIraptor output
+    if python $SCRIPT_PATH/catalogcomparisontolerancecheck.py infoFile_$TEST.txt toIfile.txt
+    then
+      echo "Catalog comparison passed: "$TEST
+    else
+      echo "Catalog comparison failed: "$TEST
+      exit 1
+    fi
+
+    echo "------------"
+
+  done
+fi
+
+exit $?
diff --git a/examples/Makefile.am b/examples/Makefile.am
index 9b9f76294f6727ac174a34dc45d993c43708fd92..e4a12fc1a90c0f52c767efd378f4bea646c0d425 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -24,7 +24,7 @@ AM_CFLAGS = -I$(top_srcdir)/src $(HDF5_CPPFLAGS) $(GSL_INCS) $(FFTW_INCS)
 AM_LDFLAGS = $(HDF5_LDFLAGS)
 
 # Extra libraries.
-EXTRA_LIBS = $(HDF5_LIBS) $(FFTW_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(TBBMALLOC_LIBS) $(GRACKLE_LIBS) $(GSL_LIBS)
+EXTRA_LIBS = $(HDF5_LIBS) $(FFTW_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(TBBMALLOC_LIBS) $(GRACKLE_LIBS) $(VELOCIRAPTOR_LIBS) $(GSL_LIBS)
 
 # MPI libraries.
 MPI_LIBS = $(METIS_LIBS) $(MPI_THREAD_LIBS)
diff --git a/examples/SmallCosmoVolume/small_cosmo_volume.yml b/examples/SmallCosmoVolume/small_cosmo_volume.yml
index b0057d536677f502c5c1e84a5af923042018a33b..32ec15db6be35fed4eb0c0168f52f0ba919158ea 100644
--- a/examples/SmallCosmoVolume/small_cosmo_volume.yml
+++ b/examples/SmallCosmoVolume/small_cosmo_volume.yml
@@ -6,6 +6,16 @@ InternalUnitSystem:
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
 
+# Structure finding options
+StructureFinding:
+  config_file_name:     stf_input_6dfof_dmonly_sub.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.02                            # Time difference between consecutive structure finding outputs (in internal units) in simulation time intervals.
+
 # WMAP9 cosmology
 Cosmology:
   Omega_m:        0.276
diff --git a/examples/SmallCosmoVolume/stf_input_6dfof_dmonly_sub.cfg b/examples/SmallCosmoVolume/stf_input_6dfof_dmonly_sub.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..e2859da17f0a71e26394d0ba40ed744d8f6ad124
--- /dev/null
+++ b/examples/SmallCosmoVolume/stf_input_6dfof_dmonly_sub.cfg
@@ -0,0 +1,164 @@
+#configuration file.
+#It is suggested that you alter this file as necessary as not all options will be desired and some conflict.
+#This file is simply meant to show options available.
+
+################################
+#input related
+################################
+#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
+
+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
+Search_for_substructure=1 #if 0, end search once field objects are found
+FoF_Field_search_type=3 #5 3DFOF search for field halos, 4 for 6DFOF clean up of field halos, 3 for 6DFOF with velocity scale distinct for each halo
+Unbind_flag=1 #run unbinding
+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
+################################
+#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
+
+################################
+#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
+#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
+Keep_FOF=0 #if field 6DFOF search is done, allows to keep structures found in 3DFOF (can be interpreted as the inter halo stellar mass when only stellar search is used).\n
+
+#minimum size for structures
+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
+
+#for mean field estimates and local velocity density distribution funciton estimator related quantiites, rarely need to change this
+Cell_fraction = 0.01 #fraction of field fof halo used to determine mean velocity distribution function. Typical values are ~0.005-0.02
+Grid_type=1 #normal entropy based grid, shouldn't have to change
+Nsearch_velocity=32 #number of velocity neighbours used to calculate local velocity distribution function. Typial values are ~32
+Nsearch_physical=256 #numerof physical neighbours from which the nearest velocity neighbour set is based. Typical values are 128-512
+
+#for substructure search, rarely ever need to change this
+FoF_search_type=1 #default phase-space FOF search. Don't really need to change
+Iterative_searchflag=1 #iterative substructure search, for substructure find initial candidate substructures with smaller linking lengths then expand search region
+Outlier_threshold=2.5 #outlier threshold for a particle to be considered residing in substructure, that is how dynamically distinct a particle is. Typical values are >2
+Velocity_ratio=2.0 #ratio of speeds used in phase-space FOF
+Velocity_opening_angle=0.10 #angle between velocities. 18 degrees here, typical values are ~10-30
+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
+Iterative_Vratio_factor=1.0 #change in Vratio when using iterative search. no change in vratio
+Iterative_ThetaOp_factor=1.0 #change in velocity opening angle. no change in velocity opening angle
+
+#for checking for halo merger remnants, which are defined as large, well separated phase-space density maxima
+
+#if searching for cores, linking lengths. likely does not need to change much
+Use_adaptive_core_search=2 #calculate dispersions in configuration & vel space to determine linking lengths
+Halo_core_ellx_fac=1.0 #how linking lengths are changed when searching for local 6DFOF cores,
+Halo_core_ellv_fac=1.0 #how velocity lengths based on dispersions are changed when searching for local 6DFOF cores
+Halo_core_ncellfac=0.05 #fraction of total halo particle number setting min size of a local 6DFOF core
+Halo_core_adaptive_sigma_fac=2.0 #used when running fully adaptive core search with phase-space tensors, specifies the width of the physical linking length in configuration space dispersion (think of this as how many sigma to include). Typically values are 2
+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)
+################################
+
+#unbinding related items
+
+Min_bound_mass_frac=0.2 #minimum bound mass fraction, not yet implemented
+#alpha factor used to determine whether particle is "bound" alaph*T+W<0. For standard subhalo catalogues use >0.9 but if interested in tidal debris 0.2-0.5
+Allowed_kinetic_potential_ratio=0.2
+#run unbinding of field structures, aka halos
+Bound_halos=0
+#simple Plummer softening length when calculating gravitational energy. If cosmological simulation with period, is fraction of interparticle spacing
+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
+################################
+#when calculating properties, for field objects calculate inclusive masses
+Inclusive_halo_masses=1 #calculate inclusive masses
+#ensures that output is comoving distances per little h
+Comoving_units=0
+
+################################
+#output related
+################################
+
+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
+#be added to the halo id as set for this snapshot, so halo 1 becomes halo 100*1000000000000+1 and halo 1 of snap 0 would just have ID=1
+Snapshot_value=1
+
+################################
+#other options
+################################
+Verbose=0 #how talkative do you want the code to be, 0 not much, 1 a lot, 2 chatterbox
diff --git a/examples/main.c b/examples/main.c
index ba36bd87032eaa04143c1f2e147cc11a6e96875c..87891ed14da81bd0e5e7cb76a5c8f41709a38666 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -105,6 +105,7 @@ void print_help_message(void) {
   printf("  %2s %14s %s\n", "-v", "[12]", "Increase the level of verbosity:");
   printf("  %2s %14s %s\n", "", "", "1: MPI-rank 0 writes,");
   printf("  %2s %14s %s\n", "", "", "2: All MPI-ranks write.");
+  printf("  %2s %14s %s\n", "-x", "", "Run with structure finding.");
   printf("  %2s %14s %s\n", "-y", "{int}",
          "Time-step frequency at which task graphs are dumped.");
   printf("  %2s %14s %s\n", "-Y", "{int}",
@@ -194,6 +195,7 @@ int main(int argc, char *argv[]) {
   int with_fp_exceptions = 0;
   int with_drift_all = 0;
   int with_mpole_reconstruction = 0;
+  int with_structure_finding = 0;
   int verbose = 0;
   int nr_threads = 1;
   int with_verbose_timers = 0;
@@ -206,7 +208,7 @@ int main(int argc, char *argv[]) {
 
   /* Parse the parameters */
   int c;
-  while ((c = getopt(argc, argv, "acCdDef:FgGhMn:o:P:rsSt:Tv:y:Y:")) != -1)
+  while ((c = getopt(argc, argv, "acCdDef:FgGhMn:o:P:rsSt:Tv:xy:Y:")) != -1)
     switch (c) {
       case 'a':
 #if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA)
@@ -303,6 +305,13 @@ int main(int argc, char *argv[]) {
           return 1;
         }
         break;
+      case 'x':
+#ifdef HAVE_VELOCIRAPTOR
+        with_structure_finding = 1;
+#else
+	error("Error: (-x) needs to have the code compiled with VELOCIraptor linked in.");
+#endif
+        break;
       case 'y':
         if (sscanf(optarg, "%d", &dump_tasks) != 1) {
           if (myrank == 0) printf("Error parsing dump_tasks (-y). \n");
@@ -483,6 +492,17 @@ int main(int argc, char *argv[]) {
   if (access(dirp, W_OK | X_OK) != 0) {
     error("Cannot write snapshots in directory %s (%s)", dirp, strerror(errno));
   }
+ 
+  /* Check that we can write the structure finding catalogues by testing if the output
+   * directory exists and is searchable and writable. */
+  if(with_structure_finding) { 
+    char stfbasename[PARSER_MAX_LINE_SIZE];
+    parser_get_param_string(params, "StructureFinding:basename", stfbasename);
+    const char *stfdirp = dirname(stfbasename);
+    if (access(stfdirp, W_OK | X_OK) != 0) {
+      error("Cannot write stf catalogues in directory %s (%s)", stfdirp, strerror(errno));
+    }
+  }
 
   /* Prepare the domain decomposition scheme */
   struct repartition reparttype;
@@ -701,6 +721,11 @@ int main(int argc, char *argv[]) {
       error("Periodic self-gravity over MPI temporarily disabled.");
 #endif
 
+#if defined(WITH_MPI) && defined(HAVE_VELOCIRAPTOR)
+    if(with_structure_finding)
+      error("VEOCIraptor not yet enabled over MPI.");
+#endif
+
 #ifdef SWIFT_DEBUG_CHECKS
     /* Check once and for all that we don't have unwanted links */
     if (!with_stars && !dry_run) {
@@ -841,6 +866,7 @@ int main(int argc, char *argv[]) {
     if (with_cooling) engine_policies |= engine_policy_cooling;
     if (with_sourceterms) engine_policies |= engine_policy_sourceterms;
     if (with_stars) engine_policies |= engine_policy_stars;
+    if (with_structure_finding) engine_policies |= engine_policy_structure_finding;
 
     /* Initialize the engine with the space and policies. */
     if (myrank == 0) clocks_gettime(&tic);
@@ -912,6 +938,10 @@ int main(int argc, char *argv[]) {
     /* Write the state of the system before starting time integration. */
     engine_dump_snapshot(&e);
     engine_print_stats(&e);
+  
+    /* Call VELOCIraptor for the first time after the first snapshot dump. */
+    //if (e.policy & engine_policy_structure_finding) velociraptor_invoke(&e);
+
   }
 
   /* Legend */
@@ -1105,6 +1135,14 @@ int main(int argc, char *argv[]) {
   engine_print_stats(&e);
   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
   if ((res = MPI_Finalize()) != MPI_SUCCESS)
     error("call to MPI_Finalize failed with error %i.", res);
diff --git a/src/Makefile.am b/src/Makefile.am
index fe17d77a874f3b2780b1ae00ba7523c3e63f51e6..6b6fb33da8a19097c6f36c6639465d435b81b542 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -48,7 +48,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.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
+    mesh_gravity.h cbrt.h velociraptor_interface.h
 
 # Common source files
 AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
@@ -60,7 +60,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
     statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \
     part_type.c xmf.c gravity_properties.c gravity.c \
     collectgroup.c hydro_space.c equation_of_state.c \
-    chemistry.c cosmology.c restart.c mesh_gravity.c
+    chemistry.c cosmology.c restart.c mesh_gravity.c velociraptor_interface.c
 
 # Include files for distribution, not installation.
 nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
@@ -159,14 +159,14 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
 libswiftsim_la_SOURCES = $(AM_SOURCES)
 libswiftsim_la_CFLAGS = $(AM_CFLAGS)
 libswiftsim_la_LDFLAGS = $(AM_LDFLAGS) $(EXTRA_LIBS)
-libswiftsim_la_LIBADD = $(GRACKLE_LIBS)
+libswiftsim_la_LIBADD = $(GRACKLE_LIBS) $(VELOCIRAPTOR_LIBS)
 
 # Sources and flags for MPI library
 libswiftsim_mpi_la_SOURCES = $(AM_SOURCES)
 libswiftsim_mpi_la_CFLAGS = $(AM_CFLAGS) $(MPI_FLAGS)
 libswiftsim_mpi_la_LDFLAGS = $(AM_LDFLAGS) $(MPI_LIBS) $(EXTRA_LIBS)
 libswiftsim_mpi_la_SHORTNAME = mpi
-libswiftsim_mpi_la_LIBADD = $(GRACKLE_LIBS)
+libswiftsim_mpi_la_LIBADD = $(GRACKLE_LIBS) $(VELOCIRAPTOR_LIBS)
 
 
 # Versioning. If any sources change then update the version_string.h file with
diff --git a/src/common_io.h b/src/common_io.h
index 61feffe2ba613c339df02e8c71cbbd1d4aec7a87..fc391b4726224022a6943c3a5cfc337ede22e4de 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -56,6 +56,15 @@ enum IO_DATA_TYPE {
   CHAR
 };
 
+/**
+ * @brief The different formats for when to run structure finding.
+ *
+ */
+enum IO_STF_OUTPUT_FORMAT {
+  STEPS = 0,
+  TIME
+};
+
 #if defined(HAVE_HDF5)
 
 hid_t io_hdf5_type(enum IO_DATA_TYPE type);
diff --git a/src/engine.c b/src/engine.c
index a8079fda8eea99a02df6baeb92017c886158fc26..66acc6b3ee41635c7da08c0614ac99712c07b1af 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -85,6 +85,7 @@
 #include "tools.h"
 #include "units.h"
 #include "version.h"
+#include "velociraptor_interface.h"
 
 /* Particle cache size. */
 #define CACHE_SIZE 512
@@ -105,7 +106,8 @@ const char *engine_policy_names[] = {"none",
                                      "reconstruct multi-poles",
                                      "cooling",
                                      "sourceterms",
-                                     "stars"};
+                                     "stars",
+                                     "structure finding"};
 
 /** The rank of the engine as a global variable (for messages). */
 int engine_rank;
@@ -4859,6 +4861,15 @@ void engine_step(struct engine *e) {
   if (e->ti_end_min >= e->ti_next_snapshot && e->ti_next_snapshot > 0)
     dump_snapshot = 1;
 
+  /* Do we want to perform structure finding? */
+  int run_stf = 0;
+  if ((e->policy & engine_policy_structure_finding)) {
+    if(e->stf_output_freq_format == STEPS && e->step % e->deltaStepSTF == 0)
+      run_stf = 1;
+    else if(e->stf_output_freq_format == TIME && e->ti_end_min >= e->ti_nextSTF && e->ti_nextSTF > 0)
+      run_stf = 1;
+  }
+
   /* Store information before attempting extra dump-related drifts */
   integertime_t ti_current = e->ti_current;
   timebin_t max_active_bin = e->max_active_bin;
@@ -4974,6 +4985,17 @@ void engine_step(struct engine *e) {
     engine_compute_next_statistics_time(e);
   }
 
+  /* Perform structure finding? */
+  if (run_stf) {
+
+    // MATTHIEU: Add a drift_all here. And check the order with the order i/o options.
+
+    velociraptor_invoke(e);
+    
+    /* ... and find the next output time */
+    if(e->stf_output_freq_format == TIME) engine_compute_next_stf_time(e);
+  }
+
   /* Restore the information we stored */
   e->ti_current = ti_current;
   e->max_active_bin = max_active_bin;
@@ -5630,7 +5652,7 @@ void engine_dump_snapshot(struct engine *e) {
 /**
  * @brief Returns the initial affinity the main thread is using.
  */
-static cpu_set_t *engine_entry_affinity(void) {
+cpu_set_t *engine_entry_affinity(void) {
 
   static int use_entry_affinity = 0;
   static cpu_set_t entry_affinity;
@@ -5862,8 +5884,25 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
   e->restart_file = restart_file;
   e->restart_next = 0;
   e->restart_dt = 0;
+  e->timeFirstSTFOutput = 0;
   engine_rank = nodeID;
 
+  /* Initialise VELOCIraptor. */
+  if (e->policy & engine_policy_structure_finding) {
+    parser_get_param_string(params, "StructureFinding:basename", e->stfBaseName);
+    e->timeFirstSTFOutput = 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) {
+      e->deltaStepSTF = parser_get_param_int(params, "StructureFinding:delta_step");
+    }
+    else if(e->stf_output_freq_format == TIME) {
+      e->deltaTimeSTF = 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);
+  }
+
   /* Get the number of queues */
   int nr_queues =
       parser_get_opt_param_int(params, "Scheduler:nr_queues", nr_threads);
@@ -6129,7 +6168,7 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
     if (e->delta_time_statistics <= 1.)
       error("Time between statistics (%e) must be > 1.",
             e->delta_time_statistics);
-
+    
     if (e->a_first_snapshot < e->cosmology->a_begin)
       error(
           "Scale-factor of first snapshot (%e) must be after the simulation "
@@ -6141,6 +6180,18 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
           "Scale-factor of first stats output (%e) must be after the "
           "simulation start a=%e.",
           e->a_first_statistics, e->cosmology->a_begin);
+    
+    if ((e->policy & engine_policy_structure_finding) && (e->stf_output_freq_format == TIME)) {
+      
+      if (e->deltaTimeSTF <= 1.)
+        error("Time between STF (%e) must be > 1.", e->deltaTimeSTF);
+
+      if (e->a_first_stf < e->cosmology->a_begin)
+        error(
+            "Scale-factor of first stf output (%e) must be after the simulation "
+            "start a=%e.",
+            e->a_first_stf, e->cosmology->a_begin);
+    }
   } else {
 
     if (e->delta_time_snapshot <= 0.)
@@ -6149,8 +6200,9 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
 
     if (e->delta_time_statistics <= 0.)
       error("Time between statistics (%e) must be positive.",
-            e->delta_time_statistics);
+          e->delta_time_statistics);
 
+    /* Find the time of the first output */
     if (e->time_first_snapshot < e->time_begin)
       error(
           "Time of first snapshot (%e) must be after the simulation start "
@@ -6162,6 +6214,26 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
           "Time of first stats output (%e) must be after the simulation start "
           "t=%e.",
           e->time_first_statistics, e->time_begin);
+
+    if ((e->policy & engine_policy_structure_finding) && (e->stf_output_freq_format == TIME)) {
+
+      if (e->deltaTimeSTF <= 0.)
+        error("Time between STF (%e) must be positive.",
+            e->deltaTimeSTF);
+
+      if (e->timeFirstSTFOutput < e->time_begin)
+        error(
+            "Time of first STF (%e) must be after the simulation start t=%e.",
+            e->timeFirstSTFOutput, e->time_begin);
+    }
+  }
+
+  if (e->policy & engine_policy_structure_finding) {
+    /* Find the time of the first stf output */
+    if(e->stf_output_freq_format == TIME) { 
+      engine_compute_next_stf_time(e);
+      message("Next STF step will be: %lld", e->ti_nextSTF);
+    }
   }
 
   /* Get the total mass */
@@ -6473,6 +6545,61 @@ void engine_compute_next_statistics_time(struct engine *e) {
   }
 }
 
+/**
+ * @brief Computes the next time (on the time line) for structure finding
+ *
+ * @param e The #engine.
+ */
+void engine_compute_next_stf_time(struct engine *e) {
+
+  /* Find upper-bound on last output */
+  double time_end;
+  if (e->policy & engine_policy_cosmology)
+    time_end = e->cosmology->a_end * e->deltaTimeSTF;
+  else
+    time_end = e->time_end + e->deltaTimeSTF;
+
+  /* Find next snasphot above current time */
+  double time = e->timeFirstSTFOutput;
+  
+  while (time < time_end) {
+
+    /* Output time on the integer timeline */
+    if (e->policy & engine_policy_cosmology)
+      e->ti_nextSTF = log(time / e->cosmology->a_begin) / e->time_base;
+    else
+      e->ti_nextSTF = (time - e->time_begin) / e->time_base;
+
+    /* Found it? */
+    if (e->ti_nextSTF > e->ti_current) break;
+
+    if (e->policy & engine_policy_cosmology)
+      time *= e->deltaTimeSTF;
+    else
+      time += e->deltaTimeSTF;
+  }
+
+  /* Deal with last snapshot */
+  if (e->ti_nextSTF >= max_nr_timesteps) {
+    e->ti_nextSTF = -1;
+    if (e->verbose) message("No further output time.");
+  } else {
+
+    /* Be nice, talk... */
+    if (e->policy & engine_policy_cosmology) {
+      const float next_snapshot_time =
+          exp(e->ti_nextSTF * e->time_base) * e->cosmology->a_begin;
+      if (e->verbose)
+        message("Next output time set to a=%e.", next_snapshot_time);
+    } else {
+      const float next_snapshot_time =
+          e->ti_nextSTF * e->time_base + e->time_begin;
+      if (e->verbose)
+        message("Next output time set to t=%e.", next_snapshot_time);
+    }
+  }
+}
+
 /**
  * @brief Computes the maximal time-step allowed by the max RMS displacement
  * condition.
@@ -6607,6 +6734,7 @@ void engine_clean(struct engine *e) {
   free(e->runners);
   free(e->snapshot_units);
   free(e->links);
+  free(e->cell_loc);
   scheduler_clean(&e->sched);
   space_clean(e->s);
   threadpool_clean(&e->threadpool);
diff --git a/src/engine.h b/src/engine.h
index ae430d9f796d6cf9bd4a1973f15b3a4e986b9997..9e44ab3a732a490104db48925904bb5629c6334d 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -70,9 +70,10 @@ enum engine_policy {
   engine_policy_reconstruct_mpoles = (1 << 12),
   engine_policy_cooling = (1 << 13),
   engine_policy_sourceterms = (1 << 14),
-  engine_policy_stars = (1 << 15)
+  engine_policy_stars = (1 << 15),
+  engine_policy_structure_finding = (1 << 16)
 };
-#define engine_maxpolicy 15
+#define engine_maxpolicy 16
 extern const char *engine_policy_names[];
 
 /**
@@ -205,6 +206,9 @@ 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;
@@ -220,6 +224,18 @@ struct engine {
   struct unit_system *snapshot_units;
   int snapshot_output_count;
 
+  /* Structure finding information */
+  int stf_output_freq_format;
+  double a_first_stf;
+  double timeFirstSTFOutput;
+  double deltaTimeSTF;
+  int deltaStepSTF;
+  
+  /* Integer time of the next stf output */
+  integertime_t ti_nextSTF;
+  
+  char stfBaseName[PARSER_MAX_LINE_SIZE];
+
   /* Statistics information */
   double a_first_statistics;
   double time_first_statistics;
@@ -348,6 +364,7 @@ struct engine {
 void engine_addlink(struct engine *e, struct link **l, struct task *t);
 void engine_barrier(struct engine *e);
 void engine_compute_next_snapshot_time(struct engine *e);
+void engine_compute_next_stf_time(struct engine *e);
 void engine_compute_next_statistics_time(struct engine *e);
 void engine_recompute_displacement_constraint(struct engine *e);
 void engine_unskip(struct engine *e);
@@ -394,6 +411,10 @@ void engine_unpin(void);
 void engine_clean(struct engine *e);
 int engine_estimate_nr_tasks(struct engine *e);
 
+#ifdef HAVE_SETAFFINITY
+cpu_set_t *engine_entry_affinity(void);
+#endif
+
 /* Struct dump/restore support. */
 void engine_struct_dump(struct engine *e, FILE *stream);
 void engine_struct_restore(struct engine *e, FILE *stream);
diff --git a/src/swift.h b/src/swift.h
index 58634443e2ed5de629123ab22a9f7560e8652107..3222376888bada072114f58ac5294f4567e8659a 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -71,6 +71,7 @@
 #include "timers.h"
 #include "tools.h"
 #include "units.h"
+#include "velociraptor_interface.h"
 #include "version.h"
 
 #endif /* SWIFT_SWIFT_H */
diff --git a/src/swift_vel_part.h b/src/swift_vel_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..f1225201755a743b821810564ef204139399f8f1
--- /dev/null
+++ b/src/swift_vel_part.h
@@ -0,0 +1,48 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 James Willis (james.s.willis@durham.ac.uk)
+ *
+ * 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_PART_H
+#define SWIFT_VELOCIRAPTOR_PART_H
+
+/* SWIFT/VELOCIraptor particle. */
+struct swift_vel_part {
+
+  /*! Particle ID. */
+  long long id;
+
+  /*! Particle position. */
+  double x[3];
+
+  /*! Particle velocity. */
+  float v[3];
+
+  /*! Particle mass. */
+  float mass;
+
+  /*! Gravitational potential */
+  float potential;
+
+  /*! Internal energy of gas particle */
+  float u;
+
+  /*! Type of the #gpart (DM, gas, star, ...) */
+  enum part_type type;
+
+};
+
+#endif /* SWIFT_VELOCIRAPTOR_PART_H */
diff --git a/src/velociraptor_interface.c b/src/velociraptor_interface.c
new file mode 100644
index 0000000000000000000000000000000000000000..f16d7df896f5463c7fdf7d87895d9f394391a9d4
--- /dev/null
+++ b/src/velociraptor_interface.c
@@ -0,0 +1,246 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 James Willis (james.s.willis@durham.ac.uk)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#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"
+
+/**
+ * @brief Initialise VELOCIraptor with input and output file names along with cosmological 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;
+
+    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);
+
+    /* 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. */
+
+    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;
+    
+    /* 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) sim_info.icosmologicalsim = 1;
+    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;
+   
+    /* 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]);
+
+    /* Initialise VELOCIraptor. */
+    if(!InitVelociraptor(configfilename, outputFileName, cosmo_info, unit_info, sim_info)) error("Exiting. VELOCIraptor initialisation failed.");
+#else
+  error("SWIFT not configure to run with VELOCIraptor.");
+#endif /* HAVE_VELOCIRAPTOR */
+
+}
+
+/**
+ * @brief Run VELOCIraptor with current particle data.
+ *
+ * @param e The #engine.
+ *
+ */
+void velociraptor_invoke(struct engine *e) {
+
+#ifdef HAVE_VELOCIRAPTOR
+    struct space *s = e->s;
+    struct gpart *gparts = s->gparts;
+    struct part *parts = s->parts;
+    const size_t nr_gparts = s->nr_gparts;
+    const size_t nr_hydro_parts = s->nr_parts;
+    const int nr_cells = s->nr_cells;
+    int *cell_node_ids = NULL;
+
+    /* 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. */
+    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_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();
+
+    /* 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);
+    
+    /* 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];
+    if(e->stf_output_freq_format == STEPS) {
+        snprintf(outputFileName, PARSER_MAX_LINE_SIZE + 128, "%s_%04i.VELOCIraptor", e->stfBaseName,
+             e->step);
+    }
+    else if(e->stf_output_freq_format == TIME) {
+        snprintf(outputFileName, PARSER_MAX_LINE_SIZE + 128, "%s_%04e.VELOCIraptor", e->stfBaseName,
+             e->time);
+    }
+
+    /* 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 a2 = e->cosmology->a * e->cosmology->a;
+    
+    message("Energy scaling factor: %f", energy_scale);
+    message("a^2: %f", a2);
+   
+    /* 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].mass = gravity_get_mass(&gparts[i]);
+      swift_parts[i].potential = gravity_get_comoving_potential(&gparts[i]); //MATTHIEU: Need factors here?
+      swift_parts[i].type = gparts[i].type;
+      
+      /* 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], 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?) !");
+      }
+    }
+
+    /* Call VELOCIraptor. */
+    if(!InvokeVelociraptor(nr_gparts, nr_hydro_parts, swift_parts, cell_node_ids, outputFileName)) error("Exiting. Call to VELOCIraptor failed on rank: %d.", e->nodeID);
+
+    /* 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);
+    
+    message("VELOCIraptor took %.3f %s on rank %d.",
+            clocks_from_ticks(getticks() - tic), clocks_getunit(), engine_rank);
+#else
+    error("SWIFT not configure to run with VELOCIraptor.");
+#endif /* HAVE_VELOCIRAPTOR */
+}
diff --git a/src/velociraptor_interface.h b/src/velociraptor_interface.h
new file mode 100644
index 0000000000000000000000000000000000000000..f568dc7c5e058ff93205e409aa6a030c97be77b1
--- /dev/null
+++ b/src/velociraptor_interface.h
@@ -0,0 +1,112 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 James Willis (james.s.willis@durham.ac.uk)
+ *
+ * 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_INTERFACE_H
+#define SWIFT_VELOCIRAPTOR_INTERFACE_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Includes. */
+#include "engine.h"
+#include "hydro.h"
+#include "swift_vel_part.h"
+
+/* Structure for passing cosmological information to VELOCIraptor. */
+struct cosmoinfo {
+
+  /*! Current expansion factor of the Universe. (cosmology.a) */
+  double atime;
+
+  /*! Reduced Hubble constant (H0 / (100km/s/Mpc) (cosmology.h) */
+  double littleh;
+
+  /*! Matter density parameter (cosmology.Omega_m) */
+  double Omega_m;
+
+  /*! Baryon density parameter (cosmology.Omega_b) */
+  double Omega_b;
+  
+  /*! Radiation constant density parameter (cosmology.Omega_lambda) */
+  double Omega_Lambda;
+
+  /*! Dark matter density parameter (cosmology.Omega_m - cosmology.Omega_b) */
+  double Omega_cdm;
+
+  /*! Dark-energy equation of state at the current time (cosmology.w)*/
+  double w_de;
+};
+
+/* Structure for passing unit information to VELOCIraptor. */
+struct unitinfo {
+
+  /* Length conversion factor to kpc. */
+  double lengthtokpc;
+  
+  /* Velocity conversion factor to km/s. */
+  double velocitytokms;
+  
+  /* Mass conversion factor to solar masses. */
+  double masstosolarmass;
+  
+  /* Potential conversion factor. */
+  double energyperunitmass;
+
+  /*! Newton's gravitationl constant (phys_const.const_newton_G)*/
+  double gravity;
+
+  /*! Hubble constant at the current redshift (cosmology.H) */
+  double hubbleunit;
+};
+
+/* Structure to hold the location of a top-level cell. */
+struct cell_loc {
+    
+    /* Coordinates x,y,z */
+    double loc[3];
+
+};
+
+/* Structure for passing simulation information to VELOCIraptor. */
+struct siminfo {
+    double period, zoomhigresolutionmass, interparticlespacing, spacedimension[3];
+    
+    /* Number of top-cells. */
+    int numcells;
+
+    /*! Locations of top-level cells. */
+    struct cell_loc *cell_loc;
+    
+    /*! Top-level cell width. */
+    double cellwidth[3];
+    
+    /*! Inverse of the top-level cell width. */
+    double icellwidth[3];
+
+    int icosmologicalsim;
+};
+
+/* 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, struct swift_vel_part *swift_parts, const int *cell_node_ids, char* output_name);
+
+/* VELOCIraptor wrapper functions. */
+void velociraptor_init(struct engine *e);
+void velociraptor_invoke(struct engine *e);
+
+#endif /* SWIFT_VELOCIRAPTOR_INTERFACE_H */