diff --git a/src/engine.c b/src/engine.c
index 6c11d9f970cda8acb48c165b33cdf9bbf185a908..fce807675ad71f176e597b585fe42b0e425ef720 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -5767,6 +5767,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
   e->delta_time_snapshot =
       parser_get_param_double(params, "Snapshots:delta_time");
   e->snaplist_snapshots = NULL;
+  e->snaplist_stats = NULL;
   e->ti_next_snapshot = 0;
   parser_get_param_string(params, "Snapshots:basename", e->snapshot_base_name);
   e->snapshot_compression =
@@ -6496,12 +6497,77 @@ void engine_compute_next_snapshot_time(struct engine *e) {
   }
 }
 
+/**
+ * @brief Read the next time (on the time line) for a dump
+ *
+ * @param e The #engine.
+ */
+void engine_read_next_snapshot_time(struct engine *e) {
+  int is_cosmo = e->policy & engine_policy_cosmology;
+  const struct snaplist *t = e->snaplist_snapshots;
+  
+  /* Find upper-bound on last output */
+  double time_end;
+  if (is_cosmo)
+    time_end = e->cosmology->a_end;
+  else
+    time_end = e->time_end;
+
+  /* Find next snasphot above current time */
+  double time;
+  time = t->times[0];
+  size_t ind = 0;
+  while (time < time_end) {
+    
+    /* Output time on the integer timeline */
+    if (is_cosmo)
+      e->ti_next_snapshot = log(time / e->cosmology->a_begin) / e->time_base;
+    else
+      e->ti_next_snapshot = (time - e->time_begin) / e->time_base;
+
+    /* Found it? */
+    if (e->ti_next_snapshot > e->ti_current) break;
+
+    ind += 1;
+    if (ind == t->size)
+      break;
+    
+    time = t->times[ind];
+  }
+
+  /* Deal with last snapshot */
+  if (e->ti_next_snapshot >= max_nr_timesteps ||
+      ind == t->size || time >= time_end) {
+    e->ti_next_snapshot = -1;
+    if (e->verbose) message("No further output time.");
+  } else {
+
+    /* Be nice, talk... */
+    if (is_cosmo) {
+      const double next_snapshot_time =
+          exp(e->ti_next_snapshot * e->time_base) * e->cosmology->a_begin;
+      if (e->verbose)
+        message("Next snapshot time set to a=%e.", next_snapshot_time);
+    } else {
+      const double next_snapshot_time =
+          e->ti_next_snapshot * e->time_base + e->time_begin;
+      if (e->verbose)
+        message("Next snapshot time set to t=%e.", next_snapshot_time);
+    }
+  }
+  
+}
+
 /**
  * @brief Computes the next time (on the time line) for a statistics dump
  *
  * @param e The #engine.
  */
 void engine_compute_next_statistics_time(struct engine *e) {
+  if (e->snaplist_stats) {
+    engine_read_next_statistics_time(e);
+    return;
+  }
 
   /* Find upper-bound on last output */
   double time_end;
@@ -6611,6 +6677,133 @@ void engine_compute_next_stf_time(struct engine *e) {
   }
 }
 
+/**
+ * @brief Read the next time (on the time line) for a statistics dump
+ *
+ * @param e The #engine.
+ */
+void engine_read_next_statistics_time(struct engine *e) {
+  int is_cosmo = e->policy & engine_policy_cosmology;
+  const struct snaplist *t = e->snaplist_stats;
+
+  /* Find upper-bound on last output */
+  double time_end;
+  if (is_cosmo)
+    time_end = e->cosmology->a_end;
+  else
+    time_end = e->time_end;
+
+  /* Find next snasphot above current time */
+  double time = t->times[0];
+  size_t ind = 0;
+  while (time < time_end) {
+
+    /* Output time on the integer timeline */
+    if (is_cosmo)
+      e->ti_next_stats = log(time / e->cosmology->a_begin) / e->time_base;
+    else
+      e->ti_next_stats = (time - e->time_begin) / e->time_base;
+
+    /* Found it? */
+    if (e->ti_next_stats > e->ti_current) break;
+
+    ind += 1;
+    if (ind == t->size)
+      break;
+    
+    time = t->times[ind];
+  }
+
+  /* Deal with last statistics */
+  if (e->ti_next_stats >= max_nr_timesteps ||
+      ind == t->size || time >= time_end) {
+    e->ti_next_stats = -1;
+    if (is_cosmo) {
+      const double next_statistics_time =
+          exp(e->ti_next_stats * e->time_base) * e->cosmology->a_begin;
+      if (e->verbose)
+        message("Next output time for stats set to a=%e.",
+                next_statistics_time);
+    } else {
+      const double next_statistics_time =
+          e->ti_next_stats * e->time_base + e->time_begin;
+      if (e->verbose)
+        message("Next output time for stats set to t=%e.",
+                next_statistics_time);
+    }
+  }
+}
+
+/**
+  * @brief Read all the snaplist files
+ *
+ * @param e The #engine.
+ * @param params The #swift_params.
+ */
+void engine_read_snaplist_files(struct engine *e, const struct swift_params *params) {
+  char filename[PARSER_MAX_LINE_SIZE];
+  struct snaplist *list;
+
+  /* get cosmo */
+  struct cosmology *cosmo = NULL;
+  if (e->policy & engine_policy_cosmology)
+    cosmo = e->cosmology;
+  
+  /* Deal with snapshots */
+  e->snaplist_snapshots = (struct snaplist*) malloc(sizeof(struct snaplist));
+  list = e->snaplist_snapshots;
+  
+  strcpy(filename, "");
+  parser_get_opt_param_string(params, "Snapshots:snaplist",
+			      filename, "");
+
+  /* Read snaplist for snapshots */
+  if (strcmp(filename, "")) {
+    message("Reading snaplist file.");
+    snaplist_read_file(list, filename, cosmo, engine_max_snaplist_snapshots);
+
+    if (list->size < 2)
+      error("You need to provide more snapshots in '%s'", filename);
+
+    /* Set data for later checks */
+    if (cosmo) {
+      e->delta_time_snapshot = list->times[1] / list->times[0];
+      e->a_first_snapshot = list->times[0];
+    }
+    else {
+      e->delta_time_snapshot = list->times[1] - list->times[0];
+      e->time_first_snapshot = list->times[0];
+    }
+  }
+
+    /* Deal with stats */
+  e->snaplist_stats = (struct snaplist*) malloc(sizeof(struct snaplist));
+  list = e->snaplist_stats;
+  
+  strcpy(filename, "");
+  parser_get_opt_param_string(params, "Statistics:snaplist",
+			      filename, "");
+
+  /* Read snaplist for stats */
+  if (strcmp(filename, "")) {
+    message("Reading snaplist file.");
+    snaplist_read_file(list, filename, cosmo, engine_max_snaplist_stats);
+
+    if (list->size < 2)
+      error("You need to provide more snapshots in '%s'", filename);
+
+    /* Set data for later checks */
+    if (cosmo) {
+      e->delta_time_statistics = list->times[1] / list->times[0];
+      e->a_first_statistics = list->times[0];
+    }
+    else {
+      e->delta_time_statistics = list->times[1] - list->times[0];
+      e->time_first_statistics = list->times[0];
+    }
+  }
+}
+
 /**
  * @brief Computes the maximal time-step allowed by the max RMS displacement
  * condition.
@@ -6748,6 +6941,10 @@ void engine_clean(struct engine *e) {
     snaplist_clean(e->snaplist_snapshots);
     free(e->snaplist_snapshots);
   }
+  if (e->snaplist_stats) {
+    snaplist_clean(e->snaplist_stats);
+    free(e->snaplist_stats);
+  }
   free(e->links);
   free(e->cell_loc);
   scheduler_clean(&e->sched);
@@ -6792,6 +6989,8 @@ void engine_struct_dump(struct engine *e, FILE *stream) {
   parser_struct_dump(e->parameter_file, stream);
   if (e->snaplist_snapshots)
     snaplist_struct_dump(e->snaplist_snapshots, stream);
+  if (e->snaplist_stats)
+    snaplist_struct_dump(e->snaplist_stats, stream);
 }
 
 /**
@@ -6894,106 +7093,14 @@ void engine_struct_restore(struct engine *e, FILE *stream) {
     e->snaplist_snapshots = snaplist_snapshots;
   }
 
+  if (e->snaplist_stats) {
+    struct snaplist *snaplist_stats =
+      (struct snaplist *) malloc(sizeof(struct snaplist));
+    snaplist_struct_restore(snaplist_stats, stream);
+    e->snaplist_stats = snaplist_stats;
+  }
+
   /* Want to force a rebuild before using this engine. Wait to repartition.*/
   e->forcerebuild = 1;
   e->forcerepart = 0;
 }
-
-/**
-  * @brief Read all the snaplist files
- *
- * @param e The #engine.
- * @param params The #swift_params.
- */
-void engine_read_snaplist_files(struct engine *e, const struct swift_params *params) {
-  char filename[PARSER_MAX_LINE_SIZE];
-  struct snaplist *list;
-
-  /* get cosmo */
-  struct cosmology *cosmo = NULL;
-  if (e->policy & engine_policy_cosmology)
-    cosmo = e->cosmology;
-  
-  /* Deal with snapshots */
-  e->snaplist_snapshots = (struct snaplist*) malloc(sizeof(struct snaplist));
-  list = e->snaplist_snapshots;
-  
-  strcpy(filename, "");
-  parser_get_opt_param_string(params, "Snapshots:snaplist",
-			      filename, "");
-
-  /* Read snaplist for snapshots */
-  if (strcmp(filename, "")) {
-    message("Reading snaplist file.");
-    snaplist_read_file(list, filename, cosmo, engine_max_snaplist_snapshots);
-
-    if (list->size < 2)
-      error("You need to provide more snapshots in '%s'", filename);
-
-    /* Set data for later checks */
-    if (cosmo) {
-      e->delta_time_snapshot = list->times[1] / list->times[0];
-      e->a_first_snapshot = list->times[0];
-    }
-    else {
-      e->delta_time_snapshot = list->times[1] - list->times[0];
-      e->time_first_snapshot = list->times[0];
-    }
-  }
-}
-
-
-void engine_read_next_snapshot_time(struct engine *e) {
-  int is_cosmo = e->policy & engine_policy_cosmology;
-  const struct snaplist *t = e->snaplist_snapshots;
-  
-  /* Find upper-bound on last output */
-  double time_end;
-  if (is_cosmo)
-    time_end = e->cosmology->a_end;
-  else
-    time_end = e->time_end;
-
-  /* Find next snasphot above current time */
-  double time;
-  time = t->times[0];
-  size_t ind = 0;
-  while (time < time_end) {
-    
-    /* Output time on the integer timeline */
-    if (is_cosmo)
-      e->ti_next_snapshot = log(time / e->cosmology->a_begin) / e->time_base;
-    else
-      e->ti_next_snapshot = (time - e->time_begin) / e->time_base;
-
-    /* Found it? */
-    if (e->ti_next_snapshot > e->ti_current) break;
-
-    ind += 1;
-    if (ind == t->size)
-      break;
-    
-    time = t->times[ind];
-  }
-
-  /* Deal with last snapshot */
-  if (e->ti_next_snapshot >= max_nr_timesteps || ind == t->size) {
-    e->ti_next_snapshot = -1;
-    if (e->verbose) message("No further output time.");
-  } else {
-
-    /* Be nice, talk... */
-    if (is_cosmo) {
-      const double next_snapshot_time =
-          exp(e->ti_next_snapshot * e->time_base) * e->cosmology->a_begin;
-      if (e->verbose)
-        message("Next snapshot time set to a=%e.", next_snapshot_time);
-    } else {
-      const double next_snapshot_time =
-          e->ti_next_snapshot * e->time_base + e->time_begin;
-      if (e->verbose)
-        message("Next snapshot time set to t=%e.", next_snapshot_time);
-    }
-  }
-  
-}
diff --git a/src/engine.h b/src/engine.h
index 55c3493ca8ac02994fbef032c71a5776bdcc8c54..3bcdcbdb0ca4e9b2bf13e89499bbd9482f274518 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -217,6 +217,7 @@ struct engine {
   double time_first_snapshot;
   double delta_time_snapshot;
   struct snaplist *snaplist_snapshots;
+  struct snaplist *snaplist_stats;
 
   /* Integer time of the next snapshot */
   integertime_t ti_next_snapshot;
@@ -370,6 +371,7 @@ void engine_compute_next_snapshot_time(struct engine *e);
 void engine_compute_next_stf_time(struct engine *e);
 void engine_read_next_snapshot_time(struct engine *e);
 void engine_compute_next_statistics_time(struct engine *e);
+void engine_read_next_statistics_time(struct engine *e);
 void engine_recompute_displacement_constraint(struct engine *e);
 void engine_unskip(struct engine *e);
 void engine_drift_all(struct engine *e);