diff --git a/examples/EAGLE_100/README b/examples/EAGLE_100/README
index e3af3c0e1281f8e9ba9e0aae3fa6dd8475359a47..d3ef58755240a73436c9b126a881412f1abbac16 100644
--- a/examples/EAGLE_100/README
+++ b/examples/EAGLE_100/README
@@ -13,4 +13,4 @@ The particle load of the main EAGLE simulation can be reproduced by
 running these ICs on 4096 cores.
 
 MD5 checksum of the ICs:
-2301ea73e14207b541bbb04163c5269e  EAGLE_ICs_100.hdf5
+bb5c3531e8573739ff4ee35049750ac9  EAGLE_ICs_100.hdf5
diff --git a/examples/EAGLE_100/eagle_100.yml b/examples/EAGLE_100/eagle_100.yml
index 651c906fb33634009f30f818e7a45d806c87fe46..f1817df36eb31a3d1ba3120afee93febb8db5e64 100644
--- a/examples/EAGLE_100/eagle_100.yml
+++ b/examples/EAGLE_100/eagle_100.yml
@@ -40,4 +40,3 @@ SPH:
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  ./EAGLE_ICs_100.hdf5     # The file to read
-  cleanup_h:  1
diff --git a/examples/main.c b/examples/main.c
index 716a3922a022f2d05c447f554444f9e4a277c073..4b4c36c68cd4420a0d831b9709eff4fd25256a70 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -32,6 +32,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#include <sys/stat.h>
 #include <unistd.h>
 
 /* MPI headers. */
@@ -92,6 +93,7 @@ void print_help_message() {
   printf("  %2s %14s %s\n", "-P", "{sec:par:val}",
          "Set parameter value and overwrites values read from the parameters "
          "file. Can be used more than once.");
+  printf("  %2s %14s %s\n", "-r", "", "Continue using restart files.");
   printf("  %2s %14s %s\n", "-s", "", "Run with hydrodynamics.");
   printf("  %2s %14s %s\n", "-S", "", "Run with stars.");
   printf("  %2s %14s %s\n", "-t", "{int}",
@@ -118,6 +120,22 @@ void print_help_message() {
 int main(int argc, char *argv[]) {
 
   struct clocks_time tic, toc;
+  struct engine e;
+
+  /* Structs used by the engine. Declare now to make sure these are always in
+   * scope.  */
+  struct chemistry_data chemistry;
+  struct cooling_function_data cooling_func;
+  struct external_potential potential;
+  struct gpart *gparts = NULL;
+  struct gravity_props gravity_properties;
+  struct hydro_props hydro_properties;
+  struct part *parts = NULL;
+  struct phys_const prog_const;
+  struct sourceterms sourceterms;
+  struct space s;
+  struct spart *sparts = NULL;
+  struct unit_system us;
 
   int nr_nodes = 1, myrank = 0;
 
@@ -161,6 +179,7 @@ int main(int argc, char *argv[]) {
   int dump_tasks = 0;
   int dump_threadpool = 0;
   int nsteps = -2;
+  int restart = 0;
   int with_cosmology = 0;
   int with_external_gravity = 0;
   int with_sourceterms = 0;
@@ -177,11 +196,12 @@ int main(int argc, char *argv[]) {
   int nparams = 0;
   char *cmdparams[PARSER_MAX_NO_OF_PARAMS];
   char paramFileName[200] = "";
+  char restart_file[200] = "";
   unsigned long long cpufreq = 0;
 
   /* Parse the parameters */
   int c;
-  while ((c = getopt(argc, argv, "acCdDef:FgGhMn:P:sSt:Tv:y:Y:")) != -1)
+  while ((c = getopt(argc, argv, "acCdDef:FgGhMn:P:rsSt:Tv:y:Y:")) != -1)
     switch (c) {
       case 'a':
 #if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA)
@@ -242,6 +262,9 @@ int main(int argc, char *argv[]) {
         cmdparams[nparams] = optarg;
         nparams++;
         break;
+      case 'r':
+        restart = 1;
+        break;
       case 's':
         with_hydro = 1;
         break;
@@ -458,227 +481,318 @@ int main(int argc, char *argv[]) {
   }
 #endif
 
-  /* Initialize unit system and constants */
-  struct unit_system us;
-  struct phys_const prog_const;
-  units_init(&us, params, "InternalUnitSystem");
-  phys_const_init(&us, &prog_const);
-  if (myrank == 0 && verbose > 0) {
-    message("Internal unit system: U_M = %e g.", us.UnitMass_in_cgs);
-    message("Internal unit system: U_L = %e cm.", us.UnitLength_in_cgs);
-    message("Internal unit system: U_t = %e s.", us.UnitTime_in_cgs);
-    message("Internal unit system: U_I = %e A.", us.UnitCurrent_in_cgs);
-    message("Internal unit system: U_T = %e K.", us.UnitTemperature_in_cgs);
-    phys_const_print(&prog_const);
+  /* Common variables for restart and IC sections. */
+  int clean_h_values = 0;
+  int flag_entropy_ICs = 0;
+
+  /* Work out where we will read and write restart files. */
+  char restart_dir[PARSER_MAX_LINE_SIZE];
+  parser_get_opt_param_string(params, "Restarts:subdir", restart_dir,
+                              "restart");
+
+  /* The directory must exist. */
+  if (myrank == 0) {
+    if (access(restart_dir, W_OK | X_OK) != 0) {
+      if (restart) {
+        error("Cannot restart as no restart subdirectory: %s (%s)", restart_dir,
+              strerror(errno));
+      } else {
+        if (mkdir(restart_dir, 0777) != 0)
+          error("Failed to create restart directory: %s (%s)", restart_dir,
+                strerror(errno));
+      }
+    }
   }
 
-  /* Initialise the hydro properties */
-  struct hydro_props hydro_properties;
-  if (with_hydro) hydro_props_init(&hydro_properties, params);
-  if (with_hydro) eos_init(&eos, params);
+  /* Basename for any restart files. */
+  char restart_name[PARSER_MAX_LINE_SIZE];
+  parser_get_opt_param_string(params, "Restarts:basename", restart_name,
+                              "swift");
 
-  /* Initialise the gravity properties */
-  struct gravity_props gravity_properties;
-  if (with_self_gravity) gravity_props_init(&gravity_properties, params);
-
-  /* Read particles and space information from (GADGET) ICs */
-  char ICfileName[200] = "";
-  parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
-  const int replicate =
-      parser_get_opt_param_int(params, "InitialConditions:replicate", 1);
-  const int clean_h_values =
-      parser_get_opt_param_int(params, "InitialConditions:cleanup_h", 0);
-  if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
-  fflush(stdout);
+  /* How often to check for the stop file and dump restarts and exit the
+   * application. */
+  int restart_stop_steps =
+      parser_get_opt_param_int(params, "Restarts:stop_steps", 100);
 
-  /* Get ready to read particles of all kinds */
-  struct part *parts = NULL;
-  struct gpart *gparts = NULL;
-  struct spart *sparts = NULL;
-  size_t Ngas = 0, Ngpart = 0, Nspart = 0;
-  double dim[3] = {0., 0., 0.};
-  int periodic = 0;
-  int flag_entropy_ICs = 0;
-  if (myrank == 0) clocks_gettime(&tic);
+  /* If restarting, look for the restart files. */
+  if (restart) {
+
+    /* Attempting a restart. */
+    char **restart_files = NULL;
+    int restart_nfiles = 0;
+
+    if (myrank == 0) {
+      message("Restarting SWIFT");
+
+      /* Locate the restart files. */
+      restart_files =
+          restart_locate(restart_dir, restart_name, &restart_nfiles);
+      if (restart_nfiles == 0)
+        error("Failed to locate any restart files in %s", restart_dir);
+
+      /* We need one file per rank. */
+      if (restart_nfiles != nr_nodes)
+        error("Incorrect number of restart files, expected %d found %d",
+              nr_nodes, restart_nfiles);
+
+      if (verbose > 0)
+        for (int i = 0; i < restart_nfiles; i++)
+          message("found restart file: %s", restart_files[i]);
+    }
+
+#ifdef WITH_MPI
+    /* Distribute the restart files, need one for each rank. */
+    if (myrank == 0) {
+
+      for (int i = 1; i < nr_nodes; i++) {
+        strcpy(restart_file, restart_files[i]);
+        MPI_Send(restart_file, 200, MPI_BYTE, i, 0, MPI_COMM_WORLD);
+      }
+
+      /* Keep local file. */
+      strcpy(restart_file, restart_files[0]);
+
+      /* Finished with the list. */
+      restart_locate_free(restart_nfiles, restart_files);
+
+    } else {
+      MPI_Recv(restart_file, 200, MPI_BYTE, 0, 0, MPI_COMM_WORLD,
+               MPI_STATUS_IGNORE);
+    }
+    if (verbose > 1) message("local restart file = %s", restart_file);
+#else
+
+    /* Just one restart file. */
+    strcpy(restart_file, restart_files[0]);
+#endif
+
+    /* Now read it. */
+    restart_read(&e, restart_file);
+
+    /* And initialize the engine with the space and policies. */
+    if (myrank == 0) clocks_gettime(&tic);
+    engine_config(1, &e, params, nr_nodes, myrank, nr_threads, with_aff,
+                  talking, restart_file);
+    if (myrank == 0) {
+      clocks_gettime(&toc);
+      message("engine_config took %.3f %s.", clocks_diff(&tic, &toc),
+              clocks_getunit());
+      fflush(stdout);
+    }
+
+    /* Check if we are already done when given steps on the command-line. */
+    if (e.step >= nsteps && nsteps > 0)
+      error("Not restarting, already completed %d steps", e.step);
+
+  } else {
+
+    /* Not restarting so look for the ICs. */
+    /* Initialize unit system and constants */
+    units_init(&us, params, "InternalUnitSystem");
+    phys_const_init(&us, &prog_const);
+    if (myrank == 0 && verbose > 0) {
+      message("Internal unit system: U_M = %e g.", us.UnitMass_in_cgs);
+      message("Internal unit system: U_L = %e cm.", us.UnitLength_in_cgs);
+      message("Internal unit system: U_t = %e s.", us.UnitTime_in_cgs);
+      message("Internal unit system: U_I = %e A.", us.UnitCurrent_in_cgs);
+      message("Internal unit system: U_T = %e K.", us.UnitTemperature_in_cgs);
+      phys_const_print(&prog_const);
+    }
+
+    /* Initialise the hydro properties */
+    if (with_hydro) hydro_props_init(&hydro_properties, params);
+    if (with_hydro) eos_init(&eos, params);
+
+    /* Initialise the gravity properties */
+    if (with_self_gravity) gravity_props_init(&gravity_properties, params);
+
+    /* Read particles and space information from (GADGET) ICs */
+    char ICfileName[200] = "";
+    parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
+    const int replicate =
+        parser_get_opt_param_int(params, "InitialConditions:replicate", 1);
+    clean_h_values =
+        parser_get_opt_param_int(params, "InitialConditions:cleanup_h", 0);
+    if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
+    fflush(stdout);
+
+    /* Get ready to read particles of all kinds */
+    size_t Ngas = 0, Ngpart = 0, Nspart = 0;
+    double dim[3] = {0., 0., 0.};
+    int periodic = 0;
+    if (myrank == 0) clocks_gettime(&tic);
 #if defined(WITH_MPI)
 #if defined(HAVE_PARALLEL_HDF5)
-  read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
+    read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
+                     &Ngpart, &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
+                     (with_external_gravity || with_self_gravity), with_stars,
+                     myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL,
+                     nr_threads, dry_run);
+#else
+    read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
                    &Ngpart, &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
                    (with_external_gravity || with_self_gravity), with_stars,
                    myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
                    dry_run);
-#else
-  read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas, &Ngpart,
-                 &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
-                 (with_external_gravity || with_self_gravity), with_stars,
-                 myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
-                 dry_run);
 #endif
 #else
-  read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas, &Ngpart,
-                 &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
-                 (with_external_gravity || with_self_gravity), with_stars,
-                 nr_threads, dry_run);
+    read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
+                   &Ngpart, &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
+                   (with_external_gravity || with_self_gravity), with_stars,
+                   nr_threads, dry_run);
 #endif
-  if (myrank == 0) {
-    clocks_gettime(&toc);
-    message("Reading initial conditions took %.3f %s.", clocks_diff(&tic, &toc),
-            clocks_getunit());
-    fflush(stdout);
-  }
+    if (myrank == 0) {
+      clocks_gettime(&toc);
+      message("Reading initial conditions took %.3f %s.",
+              clocks_diff(&tic, &toc), clocks_getunit());
+      fflush(stdout);
+    }
 
 #ifdef SWIFT_DEBUG_CHECKS
-  /* Check once and for all that we don't have unwanted links */
-  if (!with_stars) {
-    for (size_t k = 0; k < Ngpart; ++k)
-      if (gparts[k].type == swift_type_star) error("Linking problem");
-  }
-  if (!with_hydro) {
-    for (size_t k = 0; k < Ngpart; ++k)
-      if (gparts[k].type == swift_type_gas) error("Linking problem");
-  }
+    /* Check once and for all that we don't have unwanted links */
+    if (!with_stars) {
+      for (size_t k = 0; k < Ngpart; ++k)
+        if (gparts[k].type == swift_type_star) error("Linking problem");
+    }
+    if (!with_hydro) {
+      for (size_t k = 0; k < Ngpart; ++k)
+        if (gparts[k].type == swift_type_gas) error("Linking problem");
+    }
 #endif
 
-  /* Get the total number of particles across all nodes. */
-  long long N_total[3] = {0, 0, 0};
+    /* Get the total number of particles across all nodes. */
+    long long N_total[3] = {0, 0, 0};
 #if defined(WITH_MPI)
-  long long N_long[3] = {Ngas, Ngpart, Nspart};
-  MPI_Allreduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM,
-                MPI_COMM_WORLD);
+    long long N_long[3] = {Ngas, Ngpart, Nspart};
+    MPI_Allreduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM,
+                  MPI_COMM_WORLD);
 #else
-  N_total[0] = Ngas;
-  N_total[1] = Ngpart;
-  N_total[2] = Nspart;
+    N_total[0] = Ngas;
+    N_total[1] = Ngpart;
+    N_total[2] = Nspart;
 #endif
-  if (myrank == 0)
-    message(
-        "Read %lld gas particles, %lld star particles and %lld gparts from the "
-        "ICs.",
-        N_total[0], N_total[2], N_total[1]);
 
-  /* Initialize the space with these data. */
-  if (myrank == 0) clocks_gettime(&tic);
-  struct space s;
-  space_init(&s, params, dim, parts, gparts, sparts, Ngas, Ngpart, Nspart,
-             periodic, replicate, with_self_gravity, talking, dry_run);
-  if (myrank == 0) {
-    clocks_gettime(&toc);
-    message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
-            clocks_getunit());
-    fflush(stdout);
-  }
-
-  /* Say a few nice things about the space we just created. */
-  if (myrank == 0) {
-    message("space dimensions are [ %.3f %.3f %.3f ].", s.dim[0], s.dim[1],
-            s.dim[2]);
-    message("space %s periodic.", s.periodic ? "is" : "isn't");
-    message("highest-level cell dimensions are [ %i %i %i ].", s.cdim[0],
-            s.cdim[1], s.cdim[2]);
-    message("%zi parts in %i cells.", s.nr_parts, s.tot_cells);
-    message("%zi gparts in %i cells.", s.nr_gparts, s.tot_cells);
-    message("%zi sparts in %i cells.", s.nr_sparts, s.tot_cells);
-    message("maximum depth is %d.", s.maxdepth);
-    fflush(stdout);
-  }
-
-  /* Verify that each particle is in it's proper cell. */
-  if (talking && !dry_run) {
-    int icount = 0;
-    space_map_cells_pre(&s, 0, &map_cellcheck, &icount);
-    message("map_cellcheck picked up %i parts.", icount);
-  }
-
-  /* Verify the maximal depth of cells. */
-  if (talking && !dry_run) {
-    int data[2] = {s.maxdepth, 0};
-    space_map_cells_pre(&s, 0, &map_maxdepth, data);
-    message("nr of cells at depth %i is %i.", data[0], data[1]);
-  }
+    if (myrank == 0)
+      message(
+          "Read %lld gas particles, %lld star particles and %lld gparts from "
+          "the "
+          "ICs.",
+          N_total[0], N_total[2], N_total[1]);
+
+    /* Initialize the space with these data. */
+    if (myrank == 0) clocks_gettime(&tic);
+    space_init(&s, params, dim, parts, gparts, sparts, Ngas, Ngpart, Nspart,
+               periodic, replicate, with_self_gravity, talking, dry_run);
+
+    if (myrank == 0) {
+      clocks_gettime(&toc);
+      message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
+              clocks_getunit());
+      fflush(stdout);
+    }
 
-/* Initialise the table of Ewald corrections for the gravity checks */
-#ifdef SWIFT_GRAVITY_FORCE_CHECKS
-  if (periodic) gravity_exact_force_ewald_init(dim[0]);
+/* Also update the total counts (in case of changes due to replication) */
+#if defined(WITH_MPI)
+    N_long[0] = s.nr_parts;
+    N_long[1] = s.nr_gparts;
+    N_long[2] = s.nr_sparts;
+    MPI_Allreduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM,
+                  MPI_COMM_WORLD);
+#else
+    N_total[0] = s.nr_parts;
+    N_total[1] = s.nr_gparts;
+    N_total[2] = s.nr_sparts;
 #endif
 
-  /* Initialise the external potential properties */
-  struct external_potential potential;
-  if (with_external_gravity)
-    potential_init(params, &prog_const, &us, &s, &potential);
-  if (myrank == 0) potential_print(&potential);
-
-  /* Initialise the cooling function properties */
-  struct cooling_function_data cooling_func;
-  if (with_cooling) cooling_init(params, &us, &prog_const, &cooling_func);
-  if (myrank == 0) cooling_print(&cooling_func);
-
-  /* Initialise the chemistry */
-  struct chemistry_data chemistry;
-  chemistry_init(params, &us, &prog_const, &chemistry);
-  if (myrank == 0) chemistry_print(&chemistry);
+    /* Say a few nice things about the space we just created. */
+    if (myrank == 0) {
+      message("space dimensions are [ %.3f %.3f %.3f ].", s.dim[0], s.dim[1],
+              s.dim[2]);
+      message("space %s periodic.", s.periodic ? "is" : "isn't");
+      message("highest-level cell dimensions are [ %i %i %i ].", s.cdim[0],
+              s.cdim[1], s.cdim[2]);
+      message("%zi parts in %i cells.", s.nr_parts, s.tot_cells);
+      message("%zi gparts in %i cells.", s.nr_gparts, s.tot_cells);
+      message("%zi sparts in %i cells.", s.nr_sparts, s.tot_cells);
+      message("maximum depth is %d.", s.maxdepth);
+      fflush(stdout);
+    }
 
-  /* Initialise the feedback properties */
-  struct sourceterms sourceterms;
-  if (with_sourceterms) sourceterms_init(params, &us, &sourceterms);
-  if (with_sourceterms && myrank == 0) sourceterms_print(&sourceterms);
-
-  /* Construct the engine policy */
-  int engine_policies = ENGINE_POLICY | engine_policy_steal;
-  if (with_drift_all) engine_policies |= engine_policy_drift_all;
-  if (with_mpole_reconstruction)
-    engine_policies |= engine_policy_reconstruct_mpoles;
-  if (with_hydro) engine_policies |= engine_policy_hydro;
-  if (with_self_gravity) engine_policies |= engine_policy_self_gravity;
-  if (with_external_gravity) engine_policies |= engine_policy_external_gravity;
-  if (with_cosmology) engine_policies |= engine_policy_cosmology;
-  if (with_cooling) engine_policies |= engine_policy_cooling;
-  if (with_sourceterms) engine_policies |= engine_policy_sourceterms;
-  if (with_stars) engine_policies |= engine_policy_stars;
-
-  /* Initialize the engine with the space and policies. */
-  if (myrank == 0) clocks_gettime(&tic);
-  struct engine e;
-  engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, N_total[0],
-              N_total[1], with_aff, engine_policies, talking, &reparttype, &us,
-              &prog_const, &hydro_properties, &gravity_properties, &potential,
-              &cooling_func, &chemistry, &sourceterms);
-  if (myrank == 0) {
-    clocks_gettime(&toc);
-    message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
-            clocks_getunit());
-    fflush(stdout);
-  }
+    /* Verify that each particle is in it's proper cell. */
+    if (talking && !dry_run) {
+      int icount = 0;
+      space_map_cells_pre(&s, 0, &map_cellcheck, &icount);
+      message("map_cellcheck picked up %i parts.", icount);
+    }
 
-/* Init the runner history. */
-#ifdef HIST
-  for (k = 0; k < runner_hist_N; k++) runner_hist_bins[k] = 0;
-#endif
+    /* Verify the maximal depth of cells. */
+    if (talking && !dry_run) {
+      int data[2] = {s.maxdepth, 0};
+      space_map_cells_pre(&s, 0, &map_maxdepth, data);
+      message("nr of cells at depth %i is %i.", data[0], data[1]);
+    }
 
-#if defined(WITH_MPI)
-  N_long[0] = s.nr_parts;
-  N_long[1] = s.nr_gparts;
-  N_long[2] = s.nr_sparts;
-  MPI_Reduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM, 0,
-             MPI_COMM_WORLD);
-#else
-  N_total[0] = s.nr_parts;
-  N_total[1] = s.nr_gparts;
-  N_total[2] = s.nr_sparts;
-#endif
+    /* Initialise the external potential properties */
+    if (with_external_gravity)
+      potential_init(params, &prog_const, &us, &s, &potential);
+    if (myrank == 0) potential_print(&potential);
+
+    /* Initialise the cooling function properties */
+    if (with_cooling) cooling_init(params, &us, &prog_const, &cooling_func);
+    if (myrank == 0) cooling_print(&cooling_func);
+
+    /* Initialise the chemistry */
+    chemistry_init(params, &us, &prog_const, &chemistry);
+    if (myrank == 0) chemistry_print(&chemistry);
+
+    /* Initialise the feedback properties */
+    if (with_sourceterms) sourceterms_init(params, &us, &sourceterms);
+    if (with_sourceterms && myrank == 0) sourceterms_print(&sourceterms);
+
+    /* Construct the engine policy */
+    int engine_policies = ENGINE_POLICY | engine_policy_steal;
+    if (with_drift_all) engine_policies |= engine_policy_drift_all;
+    if (with_mpole_reconstruction)
+      engine_policies |= engine_policy_reconstruct_mpoles;
+    if (with_hydro) engine_policies |= engine_policy_hydro;
+    if (with_self_gravity) engine_policies |= engine_policy_self_gravity;
+    if (with_external_gravity)
+      engine_policies |= engine_policy_external_gravity;
+    if (with_cosmology) engine_policies |= engine_policy_cosmology;
+    if (with_cooling) engine_policies |= engine_policy_cooling;
+    if (with_sourceterms) engine_policies |= engine_policy_sourceterms;
+    if (with_stars) engine_policies |= engine_policy_stars;
+
+    /* Initialize the engine with the space and policies. */
+    if (myrank == 0) clocks_gettime(&tic);
+    engine_init(&e, &s, params, N_total[0], N_total[1], engine_policies,
+                talking, &reparttype, &us, &prog_const, &hydro_properties,
+                &gravity_properties, &potential, &cooling_func, &chemistry,
+                &sourceterms);
+    engine_config(0, &e, params, nr_nodes, myrank, nr_threads, with_aff,
+                  talking, restart_file);
+    if (myrank == 0) {
+      clocks_gettime(&toc);
+      message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
+              clocks_getunit());
+      fflush(stdout);
+    }
 
-  /* Get some info to the user. */
-  if (myrank == 0) {
-    long long N_DM = N_total[1] - N_total[2] - N_total[0];
-    message(
-        "Running on %lld gas particles, %lld star particles and %lld DM "
-        "particles (%lld gravity particles)",
-        N_total[0], N_total[2], N_total[1] > 0 ? N_DM : 0, N_total[1]);
-    message(
-        "from t=%.3e until t=%.3e with %d threads and %d queues (dt_min=%.3e, "
-        "dt_max=%.3e)...",
-        e.timeBegin, e.timeEnd, e.nr_threads, e.sched.nr_queues, e.dt_min,
-        e.dt_max);
-    fflush(stdout);
+    /* Get some info to the user. */
+    if (myrank == 0) {
+      long long N_DM = N_total[1] - N_total[2] - N_total[0];
+      message(
+          "Running on %lld gas particles, %lld star particles and %lld DM "
+          "particles (%lld gravity particles)",
+          N_total[0], N_total[2], N_total[1] > 0 ? N_DM : 0, N_total[1]);
+      message(
+          "from t=%.3e until t=%.3e with %d threads and %d queues "
+          "(dt_min=%.3e, "
+          "dt_max=%.3e)...",
+          e.timeBegin, e.timeEnd, e.nr_threads, e.sched.nr_queues, e.dt_min,
+          e.dt_max);
+      fflush(stdout);
+    }
   }
 
   /* Time to say good-bye if this was not a serious run. */
@@ -694,18 +808,31 @@ int main(int argc, char *argv[]) {
     return 0;
   }
 
+/* Initialise the table of Ewald corrections for the gravity checks */
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  if (periodic) gravity_exact_force_ewald_init(e.s->dim[0]);
+#endif
+
+/* Init the runner history. */
+#ifdef HIST
+  for (k = 0; k < runner_hist_N; k++) runner_hist_bins[k] = 0;
+#endif
+
+  if (!restart) {
+
 #ifdef WITH_MPI
-  /* Split the space. */
-  engine_split(&e, &initial_partition);
-  engine_redistribute(&e);
+    /* Split the space. */
+    engine_split(&e, &initial_partition);
+    engine_redistribute(&e);
 #endif
 
-  /* Initialise the particles */
-  engine_init_particles(&e, flag_entropy_ICs, clean_h_values);
+    /* Initialise the particles */
+    engine_init_particles(&e, flag_entropy_ICs, clean_h_values);
 
-  /* Write the state of the system before starting time integration. */
-  engine_dump_snapshot(&e);
-  engine_print_stats(&e);
+    /* Write the state of the system before starting time integration. */
+    engine_dump_snapshot(&e);
+    engine_print_stats(&e);
+  }
 
   /* Legend */
   if (myrank == 0)
@@ -716,8 +843,16 @@ int main(int argc, char *argv[]) {
   /* File for the timers */
   if (with_verbose_timers) timers_open_file(myrank);
 
+  /* Create a name for restart file of this rank. */
+  if (restart_genname(restart_dir, restart_name, e.nodeID, restart_file, 200) !=
+      0)
+    error("Failed to generate restart filename");
+
   /* Main simulation loop */
-  for (int j = 0; !engine_is_done(&e) && e.step - 1 != nsteps; j++) {
+  /* ==================== */
+  int force_stop = 0;
+  for (int j = 0; !engine_is_done(&e) && e.step - 1 != nsteps && !force_stop;
+       j++) {
 
     /* Reset timers */
     timers_reset_all();
@@ -728,6 +863,21 @@ int main(int argc, char *argv[]) {
     /* Print the timers. */
     if (with_verbose_timers) timers_print(e.step);
 
+    /* Every so often allow the user to stop the application and dump the
+     * restart
+     * files. */
+    if (j % restart_stop_steps == 0) {
+      force_stop = restart_stop_now(restart_dir, 0);
+      if (myrank == 0 && force_stop)
+        message("Forcing application exit, dumping restart files...");
+    }
+
+    /* Also if using nsteps to exit, will not have saved any restarts on exit,
+     * make
+     * sure we do that (useful in testing only). */
+    if (force_stop || (e.restart_onexit && e.step - 1 == nsteps))
+      engine_dump_restarts(&e, 0, 1);
+
 #ifdef SWIFT_DEBUG_TASKS
     /* Dump the task data using the given frequency. */
     if (dump_tasks && (dump_tasks == 1 || j % dump_tasks == 1)) {
@@ -852,6 +1002,10 @@ int main(int argc, char *argv[]) {
     error("call to MPI_Finalize failed with error %i.", res);
 #endif
 
+  /* Remove the stop file if used. Do this anyway, we could have missed the
+   * stop file if normal exit happened first. */
+  if (myrank == 0) force_stop = restart_stop_now(restart_dir, 1);
+
   /* Clean everything */
   if (with_verbose_timers) timers_close_file();
   engine_clean(&e);
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 19c5c95be1f1a9fcaa4e790fb24fd0f60d1ef4c2..36da7a1e9971e03124927d0dbb7fe0ca52772eb1 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -75,6 +75,15 @@ InitialConditions:
   shift_z:    0.
   replicate:  2                     # (Optional) Replicate all particles along each axis a given number of times. Default 1.
 
+# Parameters controlling restarts
+Restarts:
+  enable:      1        # (Optional) whether to enable dumping restarts at fixed intervals.
+  onexit:      0        # (Optional) whether to dump restarts on exit (*needs enable*)
+  subdir:      restart  # (Optional) name of subdirectory for restart files.
+  basename:    swift    # (Optional) prefix used in naming restart files.
+  delta_hours: 6.0      # (Optional) decimal hours between dumps of restart files.
+  stop_steps:  100      # (Optional) how many steps to process before checking if the <subdir>/stop file exists. When present the application will attempt to exit early, dumping restart files first.
+
 # Parameters governing domain decomposition
 DomainDecomposition:
   initial_type:     simple_metis # (Optional) The initial decomposition strategy: "grid",
diff --git a/src/Makefile.am b/src/Makefile.am
index 690c5a0589e7cbd71ec40c2dd2bcdd53d848bd4b..82b7aa73e10c1fa0b369886c48322985aadbc618 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -47,7 +47,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.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
+    chemistry.h chemistry_io.h chemistry_struct.h restart.h
 
 # Common source files
 AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
@@ -59,8 +59,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
-
+    chemistry.c restart.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 \
diff --git a/src/clocks.c b/src/clocks.c
index d5f0e571fe0a4694c4088bb05014fafa99d60488..4b8b269143bb3dba606c77a00eaf12b3b614eb8f 100644
--- a/src/clocks.c
+++ b/src/clocks.c
@@ -45,7 +45,7 @@ static unsigned long long clocks_cpufreq = 0;
 static ticks clocks_start = 0;
 
 /* The units of any returned times. */
-static char *clocks_units[] = {"ms", "ticks"};
+static char *clocks_units[] = {"ms", "~ms"};
 static int clocks_units_index = 0;
 static double clocks_units_scale = 1000.0;
 
@@ -176,11 +176,12 @@ static void clocks_estimate_cpufreq() {
   }
 #endif
 
-  /* If all fails just report ticks in any times. */
+  /* If all fails just report ticks for a notional 2.6GHz machine. */
   if (clocks_cpufreq == 0) {
-    clocks_cpufreq = 1;
+    unsigned long long maxfreq = 2600000;
+    clocks_cpufreq = maxfreq * 1000;
     clocks_units_index = 1;
-    clocks_units_scale = 1.0;
+    clocks_units_scale = 1000.0;
   }
 }
 
@@ -216,6 +217,22 @@ double clocks_from_ticks(ticks tics) {
   return ((double)tics / (double)clocks_get_cpufreq() * clocks_units_scale);
 }
 
+/**
+ * @brief Convert a number of milli seconds into ticks, if possible.
+ *
+ * Only an approximation as based on how well we have estimated the
+ * rtc frequency. Should be good for machines that support constant_rtc
+ * and clock_gettime(), and reasonable for most Linux machines, otherwise
+ * a guess will just be returned. See clocks_getunit() for the actual units.
+ *
+ * @param ms a number of "milliseconds" to convert to ticks
+ *
+ * @result the number of ticks, if possible.
+ */
+ticks clocks_to_ticks(double ms) {
+  return (ticks)(ms * (double)clocks_get_cpufreq() / clocks_units_scale);
+}
+
 /**
  * @brief return the time units.
  *
diff --git a/src/clocks.h b/src/clocks.h
index 4c69e130d9bfc5e61fb03fc9820ffc76d2440dc4..bdb3a6651e52f5b165e644015b91f96aa5812d57 100644
--- a/src/clocks.h
+++ b/src/clocks.h
@@ -39,6 +39,7 @@ const char *clocks_getunit();
 void clocks_set_cpufreq(unsigned long long freq);
 unsigned long long clocks_get_cpufreq();
 double clocks_from_ticks(ticks tics);
+ticks clocks_to_ticks(double interval);
 double clocks_diff_ticks(ticks tic, ticks toc);
 const char *clocks_get_timesincestart();
 
diff --git a/src/cooling.c b/src/cooling.c
index 9452857117367cfc1d4d3eab262c73eba555d4f5..89a04b30f8ce3f106d90a6db1364a3caadba6deb 100644
--- a/src/cooling.c
+++ b/src/cooling.c
@@ -22,6 +22,7 @@
 
 /* This object's header. */
 #include "cooling.h"
+#include "restart.h"
 
 /**
  * @brief Initialises the cooling properties.
@@ -52,3 +53,28 @@ void cooling_print(const struct cooling_function_data* cooling) {
 
   cooling_print_backend(cooling);
 }
+
+/**
+ * @brief Write a hydro_props struct to the given FILE as a stream of bytes.
+ *
+ * @param cooling the struct
+ * @param stream the file stream
+ */
+void cooling_struct_dump(const struct cooling_function_data* cooling,
+                         FILE* stream) {
+  restart_write_blocks((void*)cooling, sizeof(struct cooling_function_data), 1,
+                       stream, "cooling", "cooling function");
+}
+
+/**
+ * @brief Restore a hydro_props struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param cooling the struct
+ * @param stream the file stream
+ */
+void cooling_struct_restore(const struct cooling_function_data* cooling,
+                            FILE* stream) {
+  restart_read_blocks((void*)cooling, sizeof(struct cooling_function_data), 1,
+                      stream, NULL, "cooling function");
+}
diff --git a/src/cooling.h b/src/cooling.h
index ff15d95d39f857b04ed08db43b488d5f4d8ed31f..9d1001d360a1816837381e9aa52b17ba47f50fce 100644
--- a/src/cooling.h
+++ b/src/cooling.h
@@ -50,4 +50,10 @@ void cooling_init(const struct swift_params* parameter_file,
 
 void cooling_print(const struct cooling_function_data* cooling);
 
+/* Dump/restore. */
+void cooling_struct_dump(const struct cooling_function_data* cooling,
+                         FILE* stream);
+void cooling_struct_restore(const struct cooling_function_data* cooling,
+                            FILE* stream);
+
 #endif /* SWIFT_COOLING_H */
diff --git a/src/engine.h b/src/engine.h
index 02103022c6d7e38ae5e8d48054b45d60edf1a45b..25d82a5705aa78757df4bf6b2c880b71f111375b 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -83,7 +83,8 @@ enum engine_step_properties {
   engine_step_prop_redistribute = (1 << 1),
   engine_step_prop_repartition = (1 << 2),
   engine_step_prop_statistics = (1 << 3),
-  engine_step_prop_snapshot = (1 << 4)
+  engine_step_prop_snapshot = (1 << 4),
+  engine_step_prop_restarts = (1 << 5)
 };
 
 /* Some constants */
@@ -194,6 +195,7 @@ struct engine {
   char snapshotBaseName[PARSER_MAX_LINE_SIZE];
   int snapshotCompression;
   struct unit_system *snapshotUnits;
+  int snapshotOutputCount;
 
   /* Statistics information */
   FILE *file_stats;
@@ -290,6 +292,24 @@ struct engine {
   /* Temporary struct to hold a group of deferable properties (in MPI mode
    * these are reduced together, but may not be required just yet). */
   struct collectgroup1 collect_group1;
+
+  /* Whether to dump restart files. */
+  int restart_dump;
+
+  /* Whether to dump restart files after the last step. */
+  int restart_onexit;
+
+  /* Name of the restart file. */
+  const char *restart_file;
+
+  /* Ticks between restart dumps. */
+  ticks restart_dt;
+
+  /* Time after which next dump will occur. */
+  ticks restart_next;
+
+  /* Maximum number of tasks needed for restarting. */
+  int restart_max_tasks;
 };
 
 /* Function prototypes. */
@@ -301,18 +321,19 @@ void engine_drift_top_multipoles(struct engine *e);
 void engine_reconstruct_multipoles(struct engine *e);
 void engine_print_stats(struct engine *e);
 void engine_dump_snapshot(struct engine *e);
-void engine_init(struct engine *e, struct space *s,
-                 const struct swift_params *params, int nr_nodes, int nodeID,
-                 int nr_threads, long long Ngas, long long Ndm, int with_aff,
-                 int policy, int verbose, struct repartition *reparttype,
-                 const struct unit_system *internal_units,
-                 const struct phys_const *physical_constants,
-                 const struct hydro_props *hydro,
-                 const struct gravity_props *gravity,
-                 const struct external_potential *potential,
-                 const struct cooling_function_data *cooling_func,
-                 const struct chemistry_data *chemistry,
-                 struct sourceterms *sourceterms);
+void engine_init(
+    struct engine *e, struct space *s, const struct swift_params *params,
+    long long Ngas, long long Ndm, int policy, int verbose,
+    struct repartition *reparttype, const struct unit_system *internal_units,
+    const struct phys_const *physical_constants,
+    const struct hydro_props *hydro, const struct gravity_props *gravity,
+    const struct external_potential *potential,
+    const struct cooling_function_data *cooling_func,
+    const struct chemistry_data *chemistry, struct sourceterms *sourceterms);
+void engine_config(int restart, struct engine *e,
+                   const struct swift_params *params, int nr_nodes, int nodeID,
+                   int nr_threads, int with_aff, int verbose,
+                   const char *restart_file);
 void engine_launch(struct engine *e);
 void engine_prepare(struct engine *e);
 void engine_init_particles(struct engine *e, int flag_entropy_ICs,
@@ -337,4 +358,9 @@ void engine_unpin();
 void engine_clean(struct engine *e);
 int engine_estimate_nr_tasks(struct engine *e);
 
+/* Struct dump/restore support. */
+void engine_struct_dump(struct engine *e, FILE *stream);
+void engine_struct_restore(struct engine *e, FILE *stream);
+void engine_dump_restarts(struct engine *e, int drifted_all, int final_step);
+
 #endif /* SWIFT_ENGINE_H */
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index 27a5de0a4102cae4ca787c10c60cf3bbc3a983ee..1df421a8b21424dc99e52bc9dbc59560e54d14a5 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -96,3 +96,26 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
   io_write_attribute_f(h_grpgrav, "Mesh r_cut_min", p->r_cut_min);
 }
 #endif
+
+/**
+ * @brief Write a gravity_props struct to the given FILE as a stream of bytes.
+ *
+ * @param p the struct
+ * @param stream the file stream
+ */
+void gravity_props_struct_dump(const struct gravity_props *p, FILE *stream) {
+  restart_write_blocks((void *)p, sizeof(struct gravity_props), 1, stream,
+                       "gravity", "gravity props");
+}
+
+/**
+ * @brief Restore a gravity_props struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param p the struct
+ * @param stream the file stream
+ */
+void gravity_props_struct_restore(const struct gravity_props *p, FILE *stream) {
+  restart_read_blocks((void *)p, sizeof(struct gravity_props), 1, stream, NULL,
+                      "gravity props");
+}
diff --git a/src/gravity_properties.h b/src/gravity_properties.h
index f7b9950052b302a003e5d128191c9dbe68fe875f..826ffd0de05d376b930523c8a5d937e457c6d795 100644
--- a/src/gravity_properties.h
+++ b/src/gravity_properties.h
@@ -28,6 +28,7 @@
 
 /* Local includes. */
 #include "parser.h"
+#include "restart.h"
 
 /**
  * @brief Contains all the constants and parameters of the self-gravity scheme
@@ -79,4 +80,8 @@ void gravity_props_print_snapshot(hid_t h_grpsph,
                                   const struct gravity_props *p);
 #endif
 
+/* Dump/restore. */
+void gravity_props_struct_dump(const struct gravity_props *p, FILE *stream);
+void gravity_props_struct_restore(const struct gravity_props *p, FILE *stream);
+
 #endif /* SWIFT_GRAVITY_PROPERTIES */
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 995610acb2d46cb62a4fa5fc7ef76b0d9474f504..31e18913b7414110db68c6202c512de5fb90a3c1 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -127,3 +127,26 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
                        p->max_smoothing_iterations);
 }
 #endif
+
+/**
+ * @brief Write a hydro_props struct to the given FILE as a stream of bytes.
+ *
+ * @param p the struct
+ * @param stream the file stream
+ */
+void hydro_props_struct_dump(const struct hydro_props *p, FILE *stream) {
+  restart_write_blocks((void *)p, sizeof(struct hydro_props), 1, stream,
+                       "hydroprops", "hydro props");
+}
+
+/**
+ * @brief Restore a hydro_props struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param p the struct
+ * @param stream the file stream
+ */
+void hydro_props_struct_restore(const struct hydro_props *p, FILE *stream) {
+  restart_read_blocks((void *)p, sizeof(struct hydro_props), 1, stream, NULL,
+                      "hydro props");
+}
diff --git a/src/hydro_properties.h b/src/hydro_properties.h
index a887ccb6df13b649cd1ef1009059c6f08908669c..6d325c35d6ae6a9bcddbfff9ccc4601e026fc4e6 100644
--- a/src/hydro_properties.h
+++ b/src/hydro_properties.h
@@ -33,6 +33,7 @@
 
 /* Local includes. */
 #include "parser.h"
+#include "restart.h"
 
 /**
  * @brief Contains all the constants and parameters of the hydro scheme
@@ -71,4 +72,8 @@ void hydro_props_init(struct hydro_props *p, const struct swift_params *params);
 void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p);
 #endif
 
+/* Dump/restore. */
+void hydro_props_struct_dump(const struct hydro_props *p, FILE *stream);
+void hydro_props_struct_restore(const struct hydro_props *p, FILE *stream);
+
 #endif /* SWIFT_HYDRO_PROPERTIES */
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 90aeddd6af7e3a65a3e50d5cbc8c12e9c0c93fa0..0a992e82565364a4ca800e83bd66a40dce5bea1b 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -504,7 +504,6 @@ void writeArray(struct engine* e, hid_t grp, char* fileName,
 
   /* Close everything */
   H5Dclose(h_data);
-// H5Pclose(h_plist_id);
 
 #ifdef IO_SPEED_MEASUREMENT
   MPI_Barrier(MPI_COMM_WORLD);
@@ -806,7 +805,6 @@ void prepare_file(struct engine* e, const char* baseName, int outputCount,
   struct part* parts = e->s->parts;
   struct gpart* gparts = e->s->gparts;
   struct spart* sparts = e->s->sparts;
-
   FILE* xmfFile = 0;
   int periodic = e->s->periodic;
   int numFiles = 1;
@@ -820,7 +818,7 @@ void prepare_file(struct engine* e, const char* baseName, int outputCount,
   /* HDF5 File name */
   char fileName[FILENAME_BUFFER_SIZE];
   snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
-           outputCount);
+           e->snapshotOutputCount);
 
   /* Open HDF5 file with the chosen parameters */
   hid_t h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
@@ -976,7 +974,6 @@ void prepare_file(struct engine* e, const char* baseName, int outputCount,
     /* Close this particle group in the XMF file as well */
     xmf_write_groupfooter(xmfFile, (enum part_type)ptype);
   }
-  /////
 
   /* Write LXMF file descriptor */
   xmf_write_outputfooter(xmfFile, outputCount, e->time);
@@ -1269,7 +1266,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
             clocks_getunit());
 #endif
 
-  ++outputCount;
+  e->snapshotOutputCount++;
 }
 
 #endif /* HAVE_HDF5 */
diff --git a/src/parser.c b/src/parser.c
index 0b608b29263342240af68fd99d2fdd3241e2a1e6..cbd63e913010066cbd39418b81609f4669404213 100644
--- a/src/parser.c
+++ b/src/parser.c
@@ -33,6 +33,7 @@
 /* Local headers. */
 #include "common_io.h"
 #include "error.h"
+#include "restart.h"
 
 #define PARSER_COMMENT_STRING "#"
 #define PARSER_COMMENT_CHAR '#'
@@ -784,3 +785,26 @@ void parser_write_params_to_hdf5(const struct swift_params *params, hid_t grp) {
     io_write_attribute_s(grp, params->data[i].name, params->data[i].value);
 }
 #endif
+
+/**
+ * @brief Write a swift_params struct to the given FILE as a stream of bytes.
+ *
+ * @param params the struct
+ * @param stream the file stream
+ */
+void parser_struct_dump(const struct swift_params *params, FILE *stream) {
+  restart_write_blocks((void *)params, sizeof(struct swift_params), 1, stream,
+                       "parameters", "parameters");
+}
+
+/**
+ * @brief Restore a swift_params struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param params the struct
+ * @param stream the file stream
+ */
+void parser_struct_restore(const struct swift_params *params, FILE *stream) {
+  restart_read_blocks((void *)params, sizeof(struct swift_params), 1, stream,
+                      NULL, "parameters");
+}
diff --git a/src/parser.h b/src/parser.h
index bab6d8b25f5334546ac2aaf39a3f25ef7fb6ff57..4e61b16ab53b3688e60f85a42e786c44b095120a 100644
--- a/src/parser.h
+++ b/src/parser.h
@@ -83,4 +83,8 @@ void parser_get_opt_param_string(const struct swift_params *params,
 void parser_write_params_to_hdf5(const struct swift_params *params, hid_t grp);
 #endif
 
+/* Dump/restore. */
+void parser_struct_dump(const struct swift_params *params, FILE *stream);
+void parser_struct_restore(const struct swift_params *params, FILE *stream);
+
 #endif /* SWIFT_PARSER_H */
diff --git a/src/partition.c b/src/partition.c
index 297ac71ac75ea4a5c5446033418b1bb4352b116b..b46172d859e1e19221ee86ecf2c28650715281ef 100644
--- a/src/partition.c
+++ b/src/partition.c
@@ -50,6 +50,7 @@
 #include "debug.h"
 #include "error.h"
 #include "partition.h"
+#include "restart.h"
 #include "space.h"
 #include "tools.h"
 
@@ -1161,6 +1162,10 @@ void partition_init(struct partition *partition,
         "Invalid DomainDecomposition:minfrac, must be greater than 0 and less "
         "than equal to 1");
 
+  /* Clear the celllist for use. */
+  repartition->ncelllist = 0;
+  repartition->celllist = NULL;
+
 #else
   error("SWIFT was not compiled with MPI support");
 #endif
@@ -1189,7 +1194,8 @@ static int check_complete(struct space *s, int verbose, int nregions) {
     if (s->cells_top[i].nodeID <= nregions)
       present[s->cells_top[i].nodeID]++;
     else
-      message("Bad nodeID: %d", s->cells_top[i].nodeID);
+      message("Bad nodeID: s->cells_top[%d].nodeID = %d", i,
+              s->cells_top[i].nodeID);
   }
   for (int i = 0; i < nregions; i++) {
     if (!present[i]) {
@@ -1250,3 +1256,87 @@ int partition_space_to_space(double *oldh, double *oldcdim, int *oldnodeIDs,
   /* Check we have all nodeIDs present in the resample. */
   return check_complete(s, 1, nr_nodes + 1);
 }
+
+/**
+ * @brief save the nodeIDs of the current top-level cells by adding them to a
+ *             repartition struct. Used when restarting application.
+ *
+ * @param s the space with the top-level cells.
+ * @param reparttype struct to update with the a list of nodeIDs.
+ *
+ */
+void partition_store_celllist(struct space *s, struct repartition *reparttype) {
+  if (reparttype->celllist != NULL) free(reparttype->celllist);
+  reparttype->celllist = malloc(sizeof(int) * s->nr_cells);
+  reparttype->ncelllist = s->nr_cells;
+  if (reparttype->celllist == NULL) error("Failed to allocate celllist");
+
+  for (int i = 0; i < s->nr_cells; i++) {
+    reparttype->celllist[i] = s->cells_top[i].nodeID;
+  }
+}
+
+/**
+ * @brief restore the saved list of nodeIDs by applying them to the
+ *        top-level cells of a space. Used when restarting application.
+ *
+ * @param s the space with the top-level cells.
+ * @param reparttype struct with the list of nodeIDs saved,
+ *
+ */
+void partition_restore_celllist(struct space *s,
+                                struct repartition *reparttype) {
+  if (reparttype->ncelllist > 0) {
+    if (reparttype->ncelllist == s->nr_cells) {
+      for (int i = 0; i < s->nr_cells; i++) {
+        s->cells_top[i].nodeID = reparttype->celllist[i];
+      }
+      if (!check_complete(s, 1, s->e->nr_nodes)) {
+        error("Not all ranks are present in the restored partition");
+      }
+    } else {
+      error(
+          "Cannot apply the saved partition celllist as the number of"
+          "top-level cells (%d) is different to the saved number (%d)",
+          s->nr_cells, reparttype->ncelllist);
+    }
+  }
+}
+
+/**
+ * @brief Write a repartition struct to the given FILE as a stream of bytes.
+ *
+ * @param reparttype the struct
+ * @param stream the file stream
+ */
+void partition_struct_dump(struct repartition *reparttype, FILE *stream) {
+  restart_write_blocks(reparttype, sizeof(struct repartition), 1, stream,
+                       "repartition", "repartition params");
+
+  /* Also save the celllist, if we have one. */
+  if (reparttype->ncelllist > 0)
+    restart_write_blocks(reparttype->celllist,
+                         sizeof(int) * reparttype->ncelllist, 1, stream,
+                         "celllist", "repartition celllist");
+}
+
+/**
+ * @brief Restore a repartition struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param reparttype the struct
+ * @param stream the file stream
+ */
+void partition_struct_restore(struct repartition *reparttype, FILE *stream) {
+  restart_read_blocks(reparttype, sizeof(struct repartition), 1, stream, NULL,
+                      "repartition params");
+
+  /* Also restore the celllist, if we have one. */
+  if (reparttype->ncelllist > 0) {
+    reparttype->celllist = malloc(sizeof(int) * reparttype->ncelllist);
+    if (reparttype->celllist == NULL) error("Failed to allocate celllist");
+    restart_read_blocks(reparttype->celllist,
+                        sizeof(int) * reparttype->ncelllist, 1, stream, NULL,
+                        "repartition celllist");
+  }
+}
diff --git a/src/partition.h b/src/partition.h
index c3eade190c9514efb4c44011e3990745e20846fd..3ad479c6b1b343106ac736e2d9c77aa9bc93cf60 100644
--- a/src/partition.h
+++ b/src/partition.h
@@ -57,6 +57,10 @@ struct repartition {
   enum repartition_type type;
   float trigger;
   float minfrac;
+
+  /* The partition as a cell list, if used. */
+  int ncelllist;
+  int *celllist;
 };
 
 /* Simple descriptions of types for reports. */
@@ -74,4 +78,11 @@ void partition_init(struct partition *partition,
                     struct repartition *repartition,
                     const struct swift_params *params, int nr_nodes);
 
+/* Dump/restore. */
+void partition_store_celllist(struct space *s, struct repartition *reparttype);
+void partition_restore_celllist(struct space *s,
+                                struct repartition *reparttype);
+void partition_struct_dump(struct repartition *reparttype, FILE *stream);
+void partition_struct_restore(struct repartition *reparttype, FILE *stream);
+
 #endif /* SWIFT_PARTITION_H */
diff --git a/src/physical_constants.c b/src/physical_constants.c
index c851578c96e146261e9512f6899c7b82a8d91097..b70b80574e7f73fe3568e18f0809350548564c4b 100644
--- a/src/physical_constants.c
+++ b/src/physical_constants.c
@@ -27,6 +27,7 @@
 /* Local headers. */
 #include "error.h"
 #include "physical_constants_cgs.h"
+#include "restart.h"
 
 /**
  * @brief Converts physical constants to the internal unit system
@@ -34,8 +35,8 @@
  * @param us The current internal system of units.
  * @param internal_const The physical constants to initialize.
  */
-void phys_const_init(struct unit_system* us,
-                     struct phys_const* internal_const) {
+void phys_const_init(struct unit_system *us,
+                     struct phys_const *internal_const) {
 
   /* Units are declared as {U_M, U_L, U_t, U_I, U_T} */
 
@@ -105,7 +106,7 @@ void phys_const_init(struct unit_system* us,
       units_general_cgs_conversion_factor(us, dimension_length);
 }
 
-void phys_const_print(struct phys_const* internal_const) {
+void phys_const_print(struct phys_const *internal_const) {
 
   message("%25s = %e", "Gravitational constant",
           internal_const->const_newton_G);
@@ -121,3 +122,28 @@ void phys_const_print(struct phys_const* internal_const) {
   message("%25s = %e", "Parsec", internal_const->const_parsec);
   message("%25s = %e", "Solar mass", internal_const->const_solar_mass);
 }
+
+/**
+ * @brief Write a phys_const struct to the given FILE as a stream of bytes.
+ *
+ * @param internal_const the struct
+ * @param stream the file stream
+ */
+void phys_const_struct_dump(const struct phys_const *internal_const,
+                            FILE *stream) {
+  restart_write_blocks((void *)internal_const, sizeof(struct phys_const), 1,
+                       stream, "physconst", "phys_const params");
+}
+
+/**
+ * @brief Restore a phys_const struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param internal_const the struct
+ * @param stream the file stream
+ */
+void phys_const_struct_restore(const struct phys_const *internal_const,
+                               FILE *stream) {
+  restart_read_blocks((void *)internal_const, sizeof(struct phys_const), 1,
+                      stream, NULL, "phys_const params");
+}
diff --git a/src/physical_constants.h b/src/physical_constants.h
index 3731c0ef565c6159592ad2de96a222efc6cf43f2..39c6c52cba311228d6152fd59e8a6001d11a3ac1 100644
--- a/src/physical_constants.h
+++ b/src/physical_constants.h
@@ -82,4 +82,10 @@ void phys_const_init(struct unit_system* us, struct phys_const* internal_const);
 
 void phys_const_print(struct phys_const* internal_const);
 
+/* Dump/restore. */
+void phys_const_struct_dump(const struct phys_const* internal_const,
+                            FILE* stream);
+void phys_const_struct_restore(const struct phys_const* internal_const,
+                               FILE* stream);
+
 #endif /* SWIFT_PHYSICAL_CONSTANTS_H */
diff --git a/src/potential.c b/src/potential.c
index 94c2a6cc9412d1fc76b70ddebb2edba7878e6209..1fda6fc8752ff626a5262d7824ea68fd3bc16d46 100644
--- a/src/potential.c
+++ b/src/potential.c
@@ -23,6 +23,7 @@
 
 /* This object's header. */
 #include "potential.h"
+#include "restart.h"
 
 /**
  * @brief Initialises the external potential properties in the internal system
@@ -51,3 +52,29 @@ void potential_print(const struct external_potential* potential) {
 
   potential_print_backend(potential);
 }
+
+/**
+ * @brief Write an external_potential struct to the given FILE as a stream of
+ * bytes.
+ *
+ * @param potential the struct
+ * @param stream the file stream
+ */
+void potential_struct_dump(const struct external_potential* potential,
+                           FILE* stream) {
+  restart_write_blocks((void*)potential, sizeof(struct external_potential), 1,
+                       stream, "externalpotential", "external potential");
+}
+
+/**
+ * @brief Restore a external_potential struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param potential the struct
+ * @param stream the file stream
+ */
+void potential_struct_restore(const struct external_potential* potential,
+                              FILE* stream) {
+  restart_read_blocks((void*)potential, sizeof(struct external_potential), 1,
+                      stream, NULL, "external potential");
+}
diff --git a/src/potential.h b/src/potential.h
index e6ab9a5bd6bd91801eae0b3f1e3d8f65778f5065..bcb3fd284021fba339ba494c90b81c91bd2ce72f 100644
--- a/src/potential.h
+++ b/src/potential.h
@@ -50,4 +50,10 @@ void potential_init(const struct swift_params* parameter_file,
 
 void potential_print(const struct external_potential* potential);
 
+/* Dump/restore. */
+void potential_struct_dump(const struct external_potential* potential,
+                           FILE* stream);
+void potential_struct_restore(const struct external_potential* potential,
+                              FILE* stream);
+
 #endif /* SWIFT_POTENTIAL_H */
diff --git a/src/restart.c b/src/restart.c
new file mode 100644
index 0000000000000000000000000000000000000000..a507555fe62277f2bec117299f9bc12d87985153
--- /dev/null
+++ b/src/restart.c
@@ -0,0 +1,285 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 Peter W. Draper (p.w.draper@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/>.
+ *
+ ******************************************************************************/
+
+/**
+ *  @file restart.c
+ *  @brief support for SWIFT restarts
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Standard headers. */
+#include <errno.h>
+#include <glob.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/stat.h>
+#include <unistd.h>
+
+#include "engine.h"
+#include "error.h"
+#include "restart.h"
+#include "version.h"
+
+/* The signature for restart files. */
+#define SWIFT_RESTART_SIGNATURE "SWIFT-restart-file"
+#define SWIFT_RESTART_END_SIGNATURE "SWIFT-restart-file:end"
+
+#define FNAMELEN 200
+#define LABLEN 20
+
+/* Structure for a dumped header. */
+struct header {
+  size_t len;             /* Total length of data in bytes. */
+  char label[LABLEN + 1]; /* A label for data */
+};
+
+/**
+ * @brief generate a name for a restart file.
+ *
+ * @param dir the directory of restart files.
+ * @param basename the basename of the restart files.
+ * @param nodeID a unique integer, usually the nodeID of the engine.
+ * @param name pointer to a string to hold the result.
+ * @param size length of name.
+ *
+ * @result 0 if the string was large enough.
+ */
+int restart_genname(const char *dir, const char *basename, int nodeID,
+                    char *name, int size) {
+  int n = snprintf(name, size, "%s/%s_%06d.rst", dir, basename, nodeID);
+  return (n >= size);
+}
+
+/**
+ * @brief locate all the restart files in the given directory with the given
+ *        basename.
+ * @param dir the directory of restart files.
+ * @param basename the basename of the restart files.
+ * @param nfiles the number of restart files located.
+ *
+ * @result pointer to an array of strings with all the filenames found,
+ *         these should be collated using the current locale, i.e. sorted
+ *         alphabetically (so make sure the filenames are zero padded to get
+ *         numeric ordering). Release by calling restart_locate_free().
+ */
+char **restart_locate(const char *dir, const char *basename, int *nfiles) {
+  *nfiles = 0;
+
+  /* Construct the glob pattern for locating files. */
+  char pattern[FNAMELEN];
+  if (snprintf(pattern, FNAMELEN, "%s/%s_[0-9]*.rst", dir, basename) <
+      FNAMELEN) {
+
+    glob_t globbuf;
+    char **files = NULL;
+    if (glob(pattern, 0, NULL, &globbuf) == 0) {
+      *nfiles = globbuf.gl_pathc;
+      files = malloc(sizeof(char *) * *nfiles);
+      for (int i = 0; i < *nfiles; i++) {
+        files[i] = strdup(globbuf.gl_pathv[i]);
+      }
+    }
+
+    globfree(&globbuf);
+    return files;
+  }
+  error("Failed to construct pattern to locate restart files");
+
+  return NULL;
+}
+
+/**
+ * @brief Release the memory allocated to hold the restart file names.
+ *
+ * @param nfiles the number of restart files located.
+ * @param files the list of filenames found in call to restart_locate().
+ */
+void restart_locate_free(int nfiles, char **files) {
+  for (int i = 0; i < nfiles; i++) {
+    free(files[i]);
+  }
+  free(files);
+}
+
+/**
+ * @brief Write a restart file for the state of the given engine struct.
+ *
+ * @param e the engine with our state information.
+ * @param filename name of the file to write the restart data to.
+ */
+void restart_write(struct engine *e, const char *filename) {
+
+  FILE *stream = fopen(filename, "w");
+  if (stream == NULL)
+    error("Failed to open restart file: %s (%s)", filename, strerror(errno));
+
+  /* Dump our signature and version. */
+  restart_write_blocks(SWIFT_RESTART_SIGNATURE, strlen(SWIFT_RESTART_SIGNATURE),
+                       1, stream, "signature", "SWIFT signature");
+  restart_write_blocks((void *)package_version(), strlen(package_version()), 1,
+                       stream, "version", "SWIFT version");
+
+  engine_struct_dump(e, stream);
+
+  /* Just an END statement to spot truncated files. */
+  restart_write_blocks(SWIFT_RESTART_END_SIGNATURE,
+                       strlen(SWIFT_RESTART_END_SIGNATURE),
+                       1, stream, "endsignature", "SWIFT end signature");
+
+  fclose(stream);
+}
+
+/**
+ * @brief Read a restart file to construct a saved engine struct state.
+ *
+ * @param e the engine to recover from the saved state.
+ * @param filename name of the file containing the staved state.
+ */
+void restart_read(struct engine *e, const char *filename) {
+
+  FILE *stream = fopen(filename, "r");
+  if (stream == NULL)
+    error("Failed to open restart file: %s (%s)", filename, strerror(errno));
+
+  /* Get our version and signature back. These should match. */
+  char signature[strlen(SWIFT_RESTART_SIGNATURE) + 1];
+  int len = strlen(SWIFT_RESTART_SIGNATURE);
+  restart_read_blocks(signature, len, 1, stream, NULL, "SWIFT signature");
+  signature[len] = '\0';
+  if (strncmp(signature, SWIFT_RESTART_SIGNATURE, len) != 0)
+    error(
+        "Do not recognise this as a SWIFT restart file, found '%s' "
+        "expected '%s'",
+        signature, SWIFT_RESTART_SIGNATURE);
+
+  char version[FNAMELEN];
+  len = strlen(package_version());
+  restart_read_blocks(version, len, 1, stream, NULL, "SWIFT version");
+  version[len] = '\0';
+
+  /* It might work! */
+  if (strncmp(version, package_version(), len) != 0)
+    message(
+        "WARNING: restoring from a different version of SWIFT.\n You have:"
+        " '%s' and the restarts files are from: '%s'. This may fail"
+        " badly.",
+        package_version(), version);
+
+  engine_struct_restore(e, stream);
+  fclose(stream);
+}
+
+/**
+ * @brief Read blocks of memory from a file stream into a memory location.
+ *        Exits the application if the read fails and does nothing if the
+ *        size is zero.
+ *
+ * @param ptr pointer to the memory
+ * @param size size of a block
+ * @param nblocks number of blocks to read
+ * @param stream the file stream
+ * @param label the label recovered for the block, needs to be at least 20
+ *              characters, set to NULL if not required
+ * @param errstr a context string to qualify any errors.
+ */
+void restart_read_blocks(void *ptr, size_t size, size_t nblocks, FILE *stream,
+                         char *label, const char *errstr) {
+  if (size > 0) {
+    struct header head;
+    size_t nread = fread(&head, sizeof(struct header), 1, stream);
+    if (nread != 1)
+      error("Failed to read the %s header from restart file (%s)", errstr,
+            strerror(errno));
+
+    /* Check that the stored length is the same as the expected one. */
+    if (head.len != nblocks * size)
+      error("Mismatched data length in restart file for %s (%zu != %zu)",
+            errstr, head.len, nblocks * size);
+
+    /* Return label, if required. */
+    if (label != NULL) strncpy(label, head.label, LABLEN);
+
+    nread = fread(ptr, size, nblocks, stream);
+    if (nread != nblocks)
+      error("Failed to restore %s from restart file (%s)", errstr,
+            ferror(stream) ? strerror(errno) : "unexpected end of file");
+  }
+}
+
+/**
+ * @brief Write blocks of memory to a file stream from a memory location.
+ *        Exits the application if the write fails and does nothing
+ *        if the size is zero.
+ *
+ * @param ptr pointer to the memory
+ * @param size the blocks
+ * @param nblocks number of blocks to write
+ * @param stream the file stream
+ * @param label a label for the content, can only be 20 characters.
+ * @param errstr a context string to qualify any errors.
+ */
+void restart_write_blocks(void *ptr, size_t size, size_t nblocks, FILE *stream,
+                          const char *label, const char *errstr) {
+  if (size > 0) {
+
+    /* Add a preamble header. */
+    struct header head;
+    head.len = nblocks * size;
+    strncpy(head.label, label, LABLEN);
+    head.label[LABLEN] = '\0';
+
+    /* Now dump it and the data. */
+    size_t nwrite = fwrite(&head, sizeof(struct header), 1, stream);
+    if (nwrite != 1)
+      error("Failed to save %s header to restart file (%s)", errstr,
+            strerror(errno));
+
+    nwrite = fwrite(ptr, size, nblocks, stream);
+    if (nwrite != nblocks)
+      error("Failed to save %s to restart file (%s)", errstr, strerror(errno));
+  }
+}
+
+/**
+ * @brief check if the stop file exists in the given directory and optionally
+ *        remove it if found.
+ *
+ * @param dir the directory of restart files.
+ * @param cleanup remove the file if found. Should only do this from one rank
+ *                once all ranks have tested this file.
+ *
+ * @result 1 if the file was found.
+ */
+int restart_stop_now(const char *dir, int cleanup) {
+  static struct stat buf;
+  char filename[FNAMELEN];
+  strcpy(filename, dir);
+  strcat(filename, "/stop");
+  if (stat(filename, &buf) == 0) {
+    if (cleanup && unlink(filename) != 0) {
+      /* May not be fatal, so press on. */
+      message("Failed to delete restart stop file (%s)", strerror(errno));
+    }
+    return 1;
+  }
+  return 0;
+}
diff --git a/src/restart.h b/src/restart.h
new file mode 100644
index 0000000000000000000000000000000000000000..80e6c6755c09be4d3bda84860774b7da9e1a07cd
--- /dev/null
+++ b/src/restart.h
@@ -0,0 +1,42 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 Peter W. Draper (p.w.draper@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_RESTART_H
+#define SWIFT_RESTART_H
+
+#include <stdio.h>
+
+struct engine;
+
+void restart_write(struct engine *e, const char *filename);
+void restart_read(struct engine *e, const char *filename);
+
+char **restart_locate(const char *dir, const char *basename, int *nfiles);
+void restart_locate_free(int nfiles, char **files);
+int restart_genname(const char *dir, const char *basename, int nodeID,
+                    char *name, int size);
+
+void restart_read_blocks(void *ptr, size_t size, size_t nblocks, FILE *stream,
+                         char *label, const char *errstr);
+
+void restart_write_blocks(void *ptr, size_t size, size_t nblocks, FILE *stream,
+                          const char *label, const char *errstr);
+
+int restart_stop_now(const char *dir, int cleanup);
+
+#endif /* SWIFT_RESTART_H */
diff --git a/src/serial_io.c b/src/serial_io.c
index 38dd3b3e4d278535fa3f042fb3302bdefb5b5793..047c3d3f88ba56eca8720fa06bba9bf5c3db07b7 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -712,7 +712,6 @@ void write_output_serial(struct engine* e, const char* baseName,
   struct gpart* gparts = e->s->gparts;
   struct gpart* dmparts = NULL;
   struct spart* sparts = e->s->sparts;
-  static int outputCount = 0;
   FILE* xmfFile = 0;
 
   /* Number of unassociated gparts */
@@ -721,7 +720,7 @@ void write_output_serial(struct engine* e, const char* baseName,
   /* File name */
   char fileName[FILENAME_BUFFER_SIZE];
   snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
-           outputCount);
+           e->snapshotOutputCount);
 
   /* Compute offset in the file and total number of particles */
   size_t N[swift_type_count] = {Ngas, Ndm, 0, 0, Nstars, 0};
@@ -741,7 +740,7 @@ void write_output_serial(struct engine* e, const char* baseName,
   if (mpi_rank == 0) {
 
     /* First time, we need to create the XMF file */
-    if (outputCount == 0) xmf_create_file(baseName);
+    if (e->snapshotOutputCount == 0) xmf_create_file(baseName);
 
     /* Prepare the XMF file for the new entry */
     xmfFile = xmf_prepare_file(baseName);
@@ -1002,10 +1001,11 @@ void write_output_serial(struct engine* e, const char* baseName,
   }
 
   /* Write footer of LXMF file descriptor */
-  if (mpi_rank == 0) xmf_write_outputfooter(xmfFile, outputCount, e->time);
+  if (mpi_rank == 0)
+    xmf_write_outputfooter(xmfFile, e->snapshotOutputCount, e->time);
 
   /* message("Done writing particles..."); */
-  ++outputCount;
+  e->snapshotOutputCount++;
 }
 
 #endif /* HAVE_HDF5 && HAVE_MPI */
diff --git a/src/single_io.c b/src/single_io.c
index e8e281b68b183dd6c9aa9aa1a4fdfe33f56b08bf..79dc6d43d645862edb3e0f8a9a6480f368265a50 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -577,7 +577,6 @@ void write_output_single(struct engine* e, const char* baseName,
   struct gpart* gparts = e->s->gparts;
   struct gpart* dmparts = NULL;
   struct spart* sparts = e->s->sparts;
-  static int outputCount = 0;
 
   /* Number of unassociated gparts */
   const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
@@ -587,10 +586,10 @@ void write_output_single(struct engine* e, const char* baseName,
   /* File name */
   char fileName[FILENAME_BUFFER_SIZE];
   snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
-           outputCount);
+           e->snapshotOutputCount);
 
   /* First time, we need to create the XMF file */
-  if (outputCount == 0) xmf_create_file(baseName);
+  if (e->snapshotOutputCount == 0) xmf_create_file(baseName);
 
   /* Prepare the XMF file for the new entry */
   FILE* xmfFile = 0;
@@ -807,14 +806,14 @@ void write_output_single(struct engine* e, const char* baseName,
   }
 
   /* Write LXMF file descriptor */
-  xmf_write_outputfooter(xmfFile, outputCount, e->time);
+  xmf_write_outputfooter(xmfFile, e->snapshotOutputCount, e->time);
 
   /* message("Done writing particles..."); */
 
   /* Close file */
   H5Fclose(h_file);
 
-  ++outputCount;
+  e->snapshotOutputCount++;
 }
 
 #endif /* HAVE_HDF5 */
diff --git a/src/sourceterms.c b/src/sourceterms.c
index f12071cf912eae3aa8d0e25f0f3b4c5e139de667..994658740a50a764edd3988ef7d6b78e00546f8f 100644
--- a/src/sourceterms.c
+++ b/src/sourceterms.c
@@ -36,8 +36,8 @@
  * @param us The current internal system of units
  * @param source the structure that has all the source term properties
  */
-void sourceterms_init(const struct swift_params* parameter_file,
-                      struct unit_system* us, struct sourceterms* source) {
+void sourceterms_init(const struct swift_params *parameter_file,
+                      struct unit_system *us, struct sourceterms *source) {
 #ifdef SOURCETERMS_SN_FEEDBACK
   supernova_init(parameter_file, us, source);
 #endif /* SOURCETERMS_SN_FEEDBACK */
@@ -47,7 +47,7 @@ void sourceterms_init(const struct swift_params* parameter_file,
  * @brief Prints the properties of the source terms to stdout
  * @param source the structure that has all the source term properties
  */
-void sourceterms_print(struct sourceterms* source) {
+void sourceterms_print(struct sourceterms *source) {
 #ifdef SOURCETERMS_NONE
   error(" no sourceterms defined yet you ran with -F");
 #ifdef SOURCETERMS_SN_FEEDBACK
@@ -58,3 +58,28 @@ void sourceterms_print(struct sourceterms* source) {
   supernova_print(source);
 #endif /* SOURCETERMS_SN_FEEDBACK */
 };
+
+/**
+ * @brief Write a sourceterms struct to the given FILE as a stream of bytes.
+ *
+ * @param sourceterms the struct
+ * @param stream the file stream
+ */
+void sourceterms_struct_dump(const struct sourceterms *sourceterms,
+                             FILE *stream) {
+  restart_write_blocks((void *)sourceterms, sizeof(struct sourceterms), 1,
+                       stream, "sourceterms", "sourceterms");
+}
+
+/**
+ * @brief Restore a sourceterms struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param sourceterms the struct
+ * @param stream the file stream
+ */
+void sourceterms_struct_restore(const struct sourceterms *sourceterms,
+                                FILE *stream) {
+  restart_read_blocks((void *)sourceterms, sizeof(struct sourceterms), 1,
+                      stream, NULL, "sourceterms");
+}
diff --git a/src/sourceterms.h b/src/sourceterms.h
index 1445bcb777ff634d1e3a2312cb0a49ac155e1020..a5d0c3c727d70d50fb3388d5e5e3cc3d6362276f 100644
--- a/src/sourceterms.h
+++ b/src/sourceterms.h
@@ -45,6 +45,10 @@ void sourceterms_init(const struct swift_params* parameter_file,
                       struct unit_system* us, struct sourceterms* source);
 void sourceterms_print(struct sourceterms* source);
 
+/* Dump/restore. */
+void sourceterms_struct_dump(const struct sourceterms* source, FILE* stream);
+void sourceterms_struct_restore(const struct sourceterms* source, FILE* stream);
+
 /**
  * @brief Routines related to source terms
  * @param cell_min: corner of cell to test
diff --git a/src/space.c b/src/space.c
index 9e3bd998d82a8dcef851f67c9a39c024a464ef94..c0381b619f1aa8b1847fda84e3b15255cdfc9f85 100644
--- a/src/space.c
+++ b/src/space.c
@@ -53,6 +53,7 @@
 #include "memswap.h"
 #include "minmax.h"
 #include "multipole.h"
+#include "restart.h"
 #include "runner.h"
 #include "sort_part.h"
 #include "stars.h"
@@ -384,6 +385,9 @@ void space_regrid(struct space *s, int verbose) {
     }
   }
 
+  /* Are we about to allocate new top level cells without a regrid?
+   * Can happen when restarting the application. */
+  int no_regrid = (s->cells_top == NULL && oldnodeIDs == NULL);
 #endif
 
   /* Do we need to re-build the upper-level cells? */
@@ -513,6 +517,21 @@ void space_regrid(struct space *s, int verbose) {
 
       /* Finished with these. */
       free(oldnodeIDs);
+
+    } else if (no_regrid && s->e != NULL) {
+      /* If we have created the top-levels cells and not done an initial
+       * partition (can happen when restarting), then the top-level cells
+       * are not assigned to a node, we must do that and then associate the
+       * particles with the cells. Note requires that
+       * partition_store_celllist() was called once before, or just before
+       * dumping the restart files.*/
+      partition_restore_celllist(s, s->e->reparttype);
+
+      /* Now re-distribute the particles, should just add to cells? */
+      engine_redistribute(s->e);
+
+      /* Make the proxies. */
+      engine_makeproxies(s->e);
     }
 #endif /* WITH_MPI */
 
@@ -2621,7 +2640,8 @@ void space_synchronize_particle_positions(struct space *s) {
  * Calls chemistry_first_init_part() on all the particles
  */
 void space_first_init_parts(struct space *s,
-                            const struct chemistry_data *chemistry) {
+                            const struct chemistry_data *chemistry,
+                            const struct cooling_function_data *cool_func) {
 
   const size_t nr_parts = s->nr_parts;
   struct part *restrict p = s->parts;
@@ -2644,6 +2664,9 @@ void space_first_init_parts(struct space *s,
     /* Also initialise the chemistry */
     chemistry_first_init_part(&p[i], &xp[i], chemistry);
 
+    /* And the cooling */
+    cooling_first_init_part(&p[i], &xp[i], cool_func);
+
 #ifdef SWIFT_DEBUG_CHECKS
     p[i].ti_drift = 0;
     p[i].ti_kick = 0;
@@ -2651,24 +2674,6 @@ void space_first_init_parts(struct space *s,
   }
 }
 
-/**
- * @brief Initialises all the extra particle data
- *
- * Calls cooling_init_xpart() on all the particles
- */
-void space_first_init_xparts(struct space *s,
-                             const struct cooling_function_data *cool_func) {
-
-  const size_t nr_parts = s->nr_parts;
-  struct part *restrict p = s->parts;
-  struct xpart *restrict xp = s->xparts;
-
-  for (size_t i = 0; i < nr_parts; ++i) {
-
-    cooling_first_init_part(&p[i], &xp[i], cool_func);
-  }
-}
-
 /**
  * @brief Initialises all the g-particles by setting them into a valid state
  *
@@ -2695,8 +2700,8 @@ void space_first_init_gparts(struct space *s,
     gravity_first_init_gpart(&gp[i], grav_props);
 
 #ifdef SWIFT_DEBUG_CHECKS
-    gp->ti_drift = 0;
-    gp->ti_kick = 0;
+    gp[i].ti_drift = 0;
+    gp[i].ti_kick = 0;
 #endif
   }
 }
@@ -2726,8 +2731,8 @@ void space_first_init_sparts(struct space *s) {
     star_first_init_spart(&sp[i]);
 
 #ifdef SWIFT_DEBUG_CHECKS
-    sp->ti_drift = 0;
-    sp->ti_kick = 0;
+    sp[i].ti_drift = 0;
+    sp[i].ti_kick = 0;
 #endif
   }
 }
@@ -3233,3 +3238,112 @@ void space_clean(struct space *s) {
   free(s->gparts);
   free(s->sparts);
 }
+
+/**
+ * @brief Write the space struct and its contents to the given FILE as a
+ * stream of bytes.
+ *
+ * @param s the space
+ * @param stream the file stream
+ */
+void space_struct_dump(struct space *s, FILE *stream) {
+
+  restart_write_blocks(s, sizeof(struct space), 1, stream, "space",
+                       "space struct");
+
+  /* More things to write. */
+  if (s->nr_parts > 0) {
+    restart_write_blocks(s->parts, s->nr_parts, sizeof(struct part), stream,
+                         "parts", "parts");
+    restart_write_blocks(s->xparts, s->nr_parts, sizeof(struct xpart), stream,
+                         "xparts", "xparts");
+  }
+  if (s->nr_gparts > 0)
+    restart_write_blocks(s->gparts, s->nr_gparts, sizeof(struct gpart), stream,
+                         "gparts", "gparts");
+
+  if (s->nr_sparts > 0)
+    restart_write_blocks(s->sparts, s->nr_sparts, sizeof(struct spart), stream,
+                         "sparts", "sparts");
+}
+
+/**
+ * @brief Re-create a space struct and its contents from the given FILE
+ *        stream.
+ *
+ * @param s the space
+ * @param stream the file stream
+ */
+void space_struct_restore(struct space *s, FILE *stream) {
+
+  restart_read_blocks(s, sizeof(struct space), 1, stream, NULL, "space struct");
+
+  /* Things that should be reconstructed in a rebuild. */
+  s->cells_top = NULL;
+  s->cells_sub = NULL;
+  s->multipoles_top = NULL;
+  s->multipoles_sub = NULL;
+  s->local_cells_top = NULL;
+  s->grav_top_level = NULL;
+#ifdef WITH_MPI
+  s->parts_foreign = NULL;
+  s->size_parts_foreign = 0;
+  s->gparts_foreign = NULL;
+  s->size_gparts_foreign = 0;
+  s->sparts_foreign = NULL;
+  s->size_sparts_foreign = 0;
+#endif
+
+  /* More things to read. */
+  s->parts = NULL;
+  s->xparts = NULL;
+  if (s->nr_parts > 0) {
+
+    /* Need the memory for these. */
+    if (posix_memalign((void *)&s->parts, part_align,
+                       s->size_parts * sizeof(struct part)) != 0)
+      error("Failed to allocate restore part array.");
+    if (posix_memalign((void *)&s->xparts, xpart_align,
+                       s->size_parts * sizeof(struct xpart)) != 0)
+      error("Failed to allocate restore xpart array.");
+
+    restart_read_blocks(s->parts, s->nr_parts, sizeof(struct part), stream,
+                        NULL, "parts");
+    restart_read_blocks(s->xparts, s->nr_parts, sizeof(struct xpart), stream,
+                        NULL, "xparts");
+  }
+  s->gparts = NULL;
+  if (s->nr_gparts > 0) {
+    if (posix_memalign((void *)&s->gparts, gpart_align,
+                       s->size_gparts * sizeof(struct gpart)) != 0)
+      error("Failed to allocate restore gpart array.");
+
+    restart_read_blocks(s->gparts, s->nr_gparts, sizeof(struct gpart), stream,
+                        NULL, "gparts");
+  }
+
+  s->sparts = NULL;
+  if (s->nr_sparts > 0) {
+    if (posix_memalign((void *)&s->sparts, spart_align,
+                       s->size_sparts * sizeof(struct spart)) != 0)
+      error("Failed to allocate restore spart array.");
+
+    restart_read_blocks(s->sparts, s->nr_sparts, sizeof(struct spart), stream,
+                        NULL, "sparts");
+  }
+
+  /* Need to reconnect the gravity parts to their hydro and star particles. */
+  /* Re-link the parts. */
+  if (s->nr_parts > 0 && s->nr_gparts > 0)
+    part_relink_parts_to_gparts(s->gparts, s->nr_gparts, s->parts);
+
+  /* Re-link the sparts. */
+  if (s->nr_sparts > 0 && s->nr_gparts > 0)
+    part_relink_sparts_to_gparts(s->gparts, s->nr_gparts, s->sparts);
+
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Verify that everything is correct */
+  part_verify_links(s->parts, s->gparts, s->sparts, s->nr_parts, s->nr_gparts,
+                    s->nr_sparts, 1);
+#endif
+}
diff --git a/src/space.h b/src/space.h
index 06faf6bb7f25cce4b4a7799d6990d1618480349c..11cbaabdc9fc3ae17024c042ae868d464954d501 100644
--- a/src/space.h
+++ b/src/space.h
@@ -224,9 +224,8 @@ void space_do_parts_sort();
 void space_do_gparts_sort();
 void space_do_sparts_sort();
 void space_first_init_parts(struct space *s,
-                            const struct chemistry_data *chemistry);
-void space_first_init_xparts(struct space *s,
-                             const struct cooling_function_data *cool_func);
+                            const struct chemistry_data *chemistry,
+                            const struct cooling_function_data *cool_func);
 void space_first_init_gparts(struct space *s,
                              const struct gravity_props *grav_props);
 void space_first_init_sparts(struct space *s);
@@ -244,4 +243,7 @@ void space_reset_task_counters(struct space *s);
 void space_clean(struct space *s);
 void space_free_cells(struct space *s);
 
+void space_struct_dump(struct space *s, FILE *stream);
+void space_struct_restore(struct space *s, FILE *stream);
+
 #endif /* SWIFT_SPACE_H */
diff --git a/src/swift.h b/src/swift.h
index eda566b19697ceeb1d66b266b44f4769be227585..bd3c7fd5741b8c93f2931e1ec0f148d659249752 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -54,6 +54,7 @@
 #include "potential.h"
 #include "profiler.h"
 #include "queue.h"
+#include "restart.h"
 #include "runner.h"
 #include "scheduler.h"
 #include "serial_io.h"
diff --git a/src/units.c b/src/units.c
index c9038924fa540f7df86f44db885cc2ec28672a4e..5c50bb063b9411f47357407678458e6094f69e2b 100644
--- a/src/units.c
+++ b/src/units.c
@@ -39,6 +39,7 @@
 /* Includes. */
 #include "adiabatic_index.h"
 #include "error.h"
+#include "restart.h"
 
 /**
  * @brief Initialises the unit_system structure with CGS system
@@ -602,3 +603,25 @@ void units_print(const struct unit_system* us) {
   message("\tUnit Current:     %g", us->UnitCurrent_in_cgs);
   message("\tUnit Temperature: %g", us->UnitTemperature_in_cgs);
 }
+
+/**
+ * @brief Write a units struct to the given FILE as a stream of bytes.
+ *
+ * @param us the units
+ * @param stream the file stream
+ */
+void units_struct_dump(const struct unit_system* us, FILE* stream) {
+  restart_write_blocks((void*)us, sizeof(struct unit_system), 1, stream,
+                       "units", "units");
+}
+
+/**
+ * @brief Restore a units struct from the given FILE as a stream of bytes.
+ *
+ * @param us the units
+ * @param stream the file stream
+ */
+void units_struct_restore(const struct unit_system* us, FILE* stream) {
+  restart_read_blocks((void*)us, sizeof(struct unit_system), 1, stream, NULL,
+                      "units");
+}
diff --git a/src/units.h b/src/units.h
index 657b29c070f3816293c18edfc46ad6b960ec9b33..5ac70a909a77146ba5f7a441d7747acfc80c3dfa 100644
--- a/src/units.h
+++ b/src/units.h
@@ -139,4 +139,8 @@ double units_conversion_factor(const struct unit_system* from,
 
 void units_print(const struct unit_system* us);
 
+/* Dump/restore. */
+void units_struct_dump(const struct unit_system* us, FILE* stream);
+void units_struct_restore(const struct unit_system* us, FILE* stream);
+
 #endif /* SWIFT_UNITS_H */
diff --git a/theory/SPH/Derivation/sph_derivation.tex b/theory/SPH/Derivation/sph_derivation.tex
new file mode 100644
index 0000000000000000000000000000000000000000..c770035df6050afe836d943dbe6538d4eafc262b
--- /dev/null
+++ b/theory/SPH/Derivation/sph_derivation.tex
@@ -0,0 +1,219 @@
+Following \citet{hopkins2013} we use the Lagrangian formalism to determine the
+equation of motion for the SPH particles. The derivation found in
+\citet{hopkins2013} is somewhat sparse, and so we reproduce it here with more
+steps for clarity.
+
+The following derivation is underpinned by the idea of there being two
+independent ways of defining the volume associated with a particle in SPH. The
+first is the volume associated with the thermodynamical system ($\Delta V$),
+from the first law of thermodynamics, and the second being the volume around the
+particle in which we conserve an effective neighbour number ($\Delta \tilde{V}$).
+These two need not necessarily be linked in any way.
+
+We begin with the SPH lagragian,
+\begin{align}
+  L(q, \dot{q}) = \frac{1}{2}\sum^N_{i=1} m_i \dot{r}^2_i
+    - \sum^N_{i=1} m_i u_i,
+  \label{eqn:sph:derivation:sphlagrangian}
+\end{align}
+and the first law of thermodynamics,
+\begin{align}
+  \left. \frac{\partial u_i}{\partial q_i} \right|_A =
+  -\frac{P_i}{m_i} \frac{\partial \Delta V_i}{\partial q_i},
+  \label{eqn:sph:derivation:firstlaw}
+\end{align}
+where $\mathbf{q} = (\mathbf{r}_1, ..., \mathbf{r}_N, h_i, ..., h_N)$ are the
+generalised coordinates of the particles, and the derivative of the internal
+energy $u_i$ is taken at fixed entropy $A$. This gives the first concept of
+`volume' $\Delta V_i$ that the particles occupy.
+
+As mentioned earlier, the particles also have a volume associated with the
+spread of their neighbours and hence their smoothing length. We can write this as
+a constraint equation,
+\begin{align}
+    \phi_i(\mathbf{q}) = \kappa h_i^{n_d}
+      \frac{1}{\Delta \tilde{V}} - N_{ngb} = 0,
+  \label{eqn:sph:derivation:constraint}
+\end{align}
+where $N_{ngb}$ is the effecitve neighbour number, $\kappa$ is the volume of
+the unit sphere ($\kappa_{3D} = 4\pi/3$), and $n_d$ is the number of spatial
+dimensions considered in the problem. It is important to note that $N_{ngb}$
+need not be an integer quantity.
+
+\subsection{Lagrange Multipliers}
+
+With these components in hand, we can use the technique of Lagrange multipliers
+to enforce the constraint equation, with
+\begin{align}
+  \frac{\mathrm{d}}{\mathrm{d}t} \frac{\partial L}{\partial \dot{q}_i} -
+  \frac{\partial L}{\partial q_i} = 
+  \sum^{N}_{j=1} \lambda_j \frac{\partial \phi_i}{\partial q_j},
+  \label{eqn:sph:derivation:lmsum}
+\end{align}
+where the $\lambda_j$ are the lagrange multipliers.  We use the second half of
+these equations (i.e. $q_i = h_i$) to constrain $\lambda_j$. The differentials
+with respect to the smoothing lengths:
+\begin{align}
+  \frac{\partial L}{\partial \dot{h}_i} = 0, 
+  \quad
+  \frac{\partial L}{\partial h_i} =
+  -\sum^N_{j=1}m_j\frac{\partial u_j}{\partial h_i} =
+  -m_i \frac{\partial u_i}{\partial h_i}.
+\end{align}
+As all terms where $i \neq j = 0$. For the constraint equation we find
+\begin{align}
+  \frac{\partial \phi_j}{\partial h_i} = 
+  \kappa n_d h_j^{n_d -1} \frac{\partial h_j}{\partial h_i}
+    \frac{1}{\Delta \tilde{V}_j}
+  + \kappa h_j^{n_d} \frac{\partial \Delta \tilde{V}_j}{\partial h_i}
+    \frac{1}{\Delta \tilde{V}_j^2},
+\end{align}
+which clearly reduces to
+\begin{align}
+  \frac{\partial \phi_j}{\partial h_i} = \kappa
+  \left(n_d h_j^{n_d - 1}  \frac{1}{\Delta \tilde{V_j}} +
+    h_j^{n_d} \frac{\partial \Delta \tilde{V_j}}{\partial h_i}
+    \frac{1}{\Delta \tilde{V_j}^2}
+  \right)\delta_{ij},
+\end{align}
+The first law of thermodynamics gives us an expression for
+${\partial u_i}/{\partial h_i}$,
+\begin{align}
+  \frac{\partial u_i}{\partial h_i} =
+  -\frac{P_i}{m_i} \frac{\partial \Delta V_i}{\partial h_i}.
+\end{align}
+Putting all this together we find that
+\begin{align}
+  P_i \frac{\partial \Delta V_i}{\partial h_i} = 
+  \sum^N_{j=1} \kappa \lambda_j
+    & \left( n_d h_j^{n_d -1}  \frac{1}{\Delta \tilde{V_j}}\right. \nonumber\\
+    & \left. + ~  h_j^{n_d} \frac{\partial \Delta \tilde{V_j}}{\partial h_i}
+      \frac{1}{\Delta \tilde{V_j}^2}
+  \right)\delta_{ij}.
+\end{align}
+This expression can now the rearranged for $\lambda_i$,
+\begin{align}
+  \lambda_i =
+  \left(\frac{P_i}{\kappa} \frac{\Delta \tilde{V}^2_i}{h^{n_d}_i}\right)
+    \frac{h_i}{n_d \Delta \tilde{V}_i}
+    \frac{\partial \Delta V_i}{\partial h_i}
+  \left(1 - \frac{h_i}{n_d \Delta \tilde{V}_i}
+    \frac{\partial \Delta \tilde{V}_i}{\partial h_i}
+  \right)^{-1}.
+\end{align}
+In \citet{hopkins2013} the lagrange multiplier is split into two parts such that
+\begin{align}
+    \lambda_i =
+    \left(\frac{P_i}{\kappa}
+      \frac{\Delta \tilde{V}^2_i}{h^{n_d}_i}
+    \right)\psi_i,
+   \label{eqn:sph:derivation:lambda_i}
+\end{align}
+with
+\begin{align}
+  \psi_i =
+    \frac{h_i}{n_d \Delta \tilde{V}_i}
+    \frac{\partial \Delta V_i}{\partial h_i}
+    \left(1 - \frac{h_i}{n_d \Delta \tilde{V}_i}
+      \frac{\partial \Delta \tilde{V}_i}{\partial h_i}
+    \right)^{-1}.
+  \label{eqn:sph:derivation:psi_I}
+\end{align}
+
+\subsection{Generalised Equation of Motion}
+
+Now that the lagrange multipliers have been constrained, the second half of the
+equations in \ref{eqn:sph:derivation:lmsum} ($q_i = \mathbf{r}_i$) can be used
+to find the equation of motion. The differentials are given as
+\begin{align}
+    \frac{\partial L}{\partial \dot{\mathbf{r}}_i} = m_i \dot{\mathbf{r}}_i,
+    \quad
+    \frac{\partial L}{\partial \mathbf{r}_i} =
+      \sum_{j=1}^N P_j \frac{\partial \Delta V_j}{\partial \mathbf{r}_i}.
+\end{align}
+The differential of \ref{eqn:sph:derivation:constraint} is
+\begin{align}
+    \frac{\partial \phi_j}{\partial \mathbf{r}_i} =
+      - \kappa h_j^{n_d} \frac{\Delta \tilde{V}}{\partial \mathbf{r}_i}
+        \frac{1}{\Delta \tilde{V}^2},
+\end{align}
+as $h_j$ is an \emph{independent coordinate} of a particle from $\mathbf{r}_j$
+implying that the partial differential $\partial h_j/\partial \mathbf{r}_i = 0$.
+We can substitute these into the equation of motion to find
+\begin{align}
+  m_i \frac{\mathrm{d} \mathbf{v}_i}{\mathrm{d}t} =  
+    \sum^N_{j=1}  P_j \nabla_i \Delta V_j + \lambda_j
+      \left(
+        \kappa \frac{h_j^{n_d}}{\Delta \tilde{V}_j^2}
+      \right)
+      \left(
+        - \frac{\partial \Delta \tilde{V}_j}{\partial \mathbf{r}_i}
+      \right).
+\end{align}
+This then gives the result from the substitution of the lagrange multipliers, 
+\begin{align}
+  m_i \frac{\mathrm{d} \mathbf{v}_i}{\mathrm{d}t} =
+  \sum^N_{j=1} P_j \nabla_i \Delta \tilde{V}_j \psi_j +
+  P_j \nabla_i \Delta V_j.
+  \label{eqn:sph:derivation:eom}
+\end{align}
+
+Now we need to \emph{specify} what we mean by volumes in terms of SPH
+quantities. An example, familiar choice would be $\Delta V_i = m_i/\rho_i$.
+Notice that in this definition of volume there is one \emph{particle-carried}
+property and one `field' or \emph{smoothed} property. This turns out to be the
+case for any smoothed property
+\begin{align}
+  y_i = \sum^N_{j=1} x_j W_{ij}(h_i)
+  \label{eqn:sph:derivation:smoothed}
+\end{align}
+for a given particle-carried scalar $x_i$. This is also true for the $\Delta
+\tilde{V}_i = \tilde{x}_i/\tilde{y}_i$. Using these specifications, we can write
+down the volume differentials,
+\begin{align}
+  \frac{\partial \Delta V_i}{\partial h_i} =
+    -\frac{x_i}{y_i^2}\frac{\partial y_i}{\partial h_i},
+    \quad
+  \nabla_i \Delta V_j = -\frac{x_i}{y_i^2} \nabla_i y_j.
+  \label{eqn:sph:derivation:volumediffs}
+\end{align}
+The spatial differential is fairly straightforward and is given by
+\begin{align}
+  \nabla_i y_j = \nabla_i W_{ij}(h_j)
+    + \delta_{ij}\sum_{k=1}^N \nabla_i W_{ik}(h_i).
+  \label{eqn:sph:derivation:nablay}
+\end{align}
+The differential with respect to the smoothing length is also straightforward
+after remembering that $W_{ij}(h_j) = w(|r_{ij}|/h_j)/h_j^{n_d}$. Then,
+\begin{align}
+  \frac{\partial y_i}{\partial h_i} = -\sum_{j=1}^N \frac{x_j}{h_j}
+  \left[
+    n_d W_{ij}(h_i) + \frac{|r_{ij}|}{h_i}
+    \left. 
+      \frac{\partial W(u)}{\partial u}
+    \right|_{u=\frac{|r_{ij}|}{h_i}}
+  \right].
+  \label{eqn:sph:derivation:dydh}
+\end{align}
+Finally, putting all of the above together, we can reach a
+formulation-independent equation of motion for SPH,
+\begin{align}
+  \frac{\mathrm{d}\vec{v}_i}{\mathrm{d}t} = -\sum_{j=1}^N x_i x_j
+  & \left[ \frac{f_{ij}P_i}{y_i^2} \nabla_i W_{ij}(h_i) \right. \nonumber \\
+  & \left. + ~ \frac{f_{ji} P_j}{y_j^2}\nabla_i W_{ji}(h_j)\right],
+  \label{eqn:sph:derivation:spheom}
+\end{align}
+with
+\begin{align}
+  f_{ij} \equiv 1 - 
+    \frac{\tilde{x}_j}{x_j}
+    \left(
+      \frac{h_i}{n_d \tilde{y}_i} \frac{\partial y_i}{\partial h_i}
+    \right)
+    \left(
+      1+\frac{h_i}{n_d \tilde{y}_i} \frac{\partial \tilde{y}_i}{\partial h_i}
+    \right)^{-1}
+  \label{eqn:sph:derivation:fij}
+\end{align}
+being the so-called `h-terms' that are essentially correction factors due to the
+fact that $h$ is allowed to vary.
diff --git a/theory/SPH/Flavours/sph_flavours.tex b/theory/SPH/Flavours/sph_flavours.tex
index 5fe1277373552d60607671299437a371e068169c..3c80fefb4989505b76cfaf8b38676ac4276b8da8 100644
--- a/theory/SPH/Flavours/sph_flavours.tex
+++ b/theory/SPH/Flavours/sph_flavours.tex
@@ -36,8 +36,9 @@ and the derivative of its density with respect to $h$:
   \rho_{\partial h_i} \equiv \dd{\rho}{h}(\vec{x}_i) = \sum_j m_j \dd{W}{h}(\vec{x}_{ij}
   , h_i).
 \end{equation}
-The gradient terms (``h-terms'') can then be computed from the density
-and its derivative:
+This corresponds to $x_i = \tilde{x}_i = m_i$, and $y_i =\tilde{y}_i = \rho_i$
+in the \citet{hopkins2013} formalism.  The gradient terms (``h-terms'') can
+then be computed from the density and its derivative:
 
 \begin{equation}
   f_i \equiv \left(1 + \frac{h_i}{3\rho_i}\rho_{\partial h_i}
@@ -300,15 +301,17 @@ smoothing length using:
 \bar P_{\partial h_i} \equiv \dd{\bar{P}}{h}(\vec{x}_i) = \sum_j m_j
 \tilde{A_j} \dd{W}{h}(\vec{x}_{ij}), \label{eq:sph:pe:P_dh}
 \end{equation}
+This corresponds to $x_i = m_i \tilde{A}_i$, $\tilde{x}_i = m_i$, $y_i =
+\bar{P}_i$, and $\tilde{y}_i = \rho_i$ in the \citet{hopkins2013} formalism.
 The gradient terms (``h-terms'') are then obtained by combining $\bar
-P_{\partial h_i}$ and $\rho_{\partial h_i}$
-(eq. \ref{eq:sph:minimal:rho_dh}):
+P_{\partial h_i}$ and $\rho_{\partial h_i}$ (eq. \ref{eq:sph:minimal:rho_dh}):
 
-\begin{equation}
-  f_i \equiv \left(\frac{h_i}{3\rho_i}\bar P_{\partial
+\begin{align}
+    f_{ij} = & ~ 1 - \tilde{A}_j^{-1} f_i \nonumber \\
+    f_i \equiv &  \left(\frac{h_i}{3\rho_i}\bar P_{\partial
     h_i}\right)\left(1 + \frac{h_i}{3\rho_i}\rho_{\partial
     h_i}\right)^{-1}. 
-\end{equation}
+\end{align}
 
 \subsubsection{Hydrodynamical accelerations (\nth{2} neighbour loop)}
 
diff --git a/theory/SPH/Kernels/kernels.py b/theory/SPH/Kernels/kernels.py
index 069bfd1ea25c8b99b894eadb46f93b9656ba9c7e..1d3a3524ba5d62ddf5a22a641985c119d8c7703e 100644
--- a/theory/SPH/Kernels/kernels.py
+++ b/theory/SPH/Kernels/kernels.py
@@ -20,12 +20,12 @@ import matplotlib
 matplotlib.use("Agg")
 from pylab import *
 from scipy import integrate
-import distinct_colours as colours
 from scipy.optimize import curve_fit
 from scipy.optimize import fsolve
 from matplotlib.font_manager import FontProperties
 import numpy
 
+
 params = {'axes.labelsize': 9,
 'axes.titlesize': 10,
 'font.size': 12,
@@ -100,12 +100,12 @@ N_H_WendlandC4 = 4./3. * PI * H_WendlandC4**3 / (dx)**3
 N_H_WendlandC6 = 4./3. * PI * H_WendlandC6**3 / (dx)**3
 
 
-print "Smoothing length: h =", h, "Cubic spline kernel support size:   H =", H_cubic, "Number of neighbours N_H =", N_H_cubic
-print "Smoothing length: h =", h, "Quartic spline kernel support size: H =", H_quartic, "Number of neighbours N_H =", N_H_quartic
-print "Smoothing length: h =", h, "Quintic spline kernel support size: H =", H_quintic, "Number of neighbours N_H =", N_H_quintic
-print "Smoothing length: h =", h, "Wendland C2 kernel support size:    H =", H_WendlandC2, "Number of neighbours N_H =", N_H_WendlandC2
-print "Smoothing length: h =", h, "Wendland C4 kernel support size:    H =", H_WendlandC4, "Number of neighbours N_H =", N_H_WendlandC4
-print "Smoothing length: h =", h, "Wendland C6 kernel support size:    H =", H_WendlandC6, "Number of neighbours N_H =", N_H_WendlandC6
+print("Smoothing length: h =", h, "Cubic spline kernel support size:   H =", H_cubic, "Number of neighbours N_H =", N_H_cubic)
+print("Smoothing length: h =", h, "Quartic spline kernel support size: H =", H_quartic, "Number of neighbours N_H =", N_H_quartic)
+print("Smoothing length: h =", h, "Quintic spline kernel support size: H =", H_quintic, "Number of neighbours N_H =", N_H_quintic)
+print("Smoothing length: h =", h, "Wendland C2 kernel support size:    H =", H_WendlandC2, "Number of neighbours N_H =", N_H_WendlandC2)
+print("Smoothing length: h =", h, "Wendland C4 kernel support size:    H =", H_WendlandC4, "Number of neighbours N_H =", N_H_WendlandC4)
+print("Smoothing length: h =", h, "Wendland C6 kernel support size:    H =", H_WendlandC6, "Number of neighbours N_H =", N_H_WendlandC6)
 
 # Get kernel constants (Dehen & Aly 2012, table 1) for 3D kernel
 C_cubic   = 16. / PI
diff --git a/theory/SPH/swift_sph.tex b/theory/SPH/swift_sph.tex
index 693df8cbc7cc9d03e98bdc80725572016ccf7bd0..4be8b965b2888bb04b5c150b41ccc38ef2ec3f95 100644
--- a/theory/SPH/swift_sph.tex
+++ b/theory/SPH/swift_sph.tex
@@ -18,7 +18,7 @@
 
 %opening
 \title{SPH implementation in \swift}
-\author{Matthieu Schaller}
+\author{Matthieu Schaller, Josh Borrow}
 
 \begin{document}
 
@@ -29,6 +29,9 @@
 \section{Equation of state}
 \input{EoS/eos}
 
+\section{Derivation of the Equation of Motion}
+\input{Derivation/sph_derivation.tex}
+
 \section{SPH flavours}
 \input{Flavours/sph_flavours}