diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index fb1e6ff6931b2fdee00792eed7f178872ee6d950..b76b127f02f8253d5738ebdaae2a685dc74f96e0 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -33,7 +33,9 @@ Snapshots:
 
 # Parameters governing the conserved quantities statistics
 Statistics:
-  delta_time:          1e-2 # Time between statistics output
+  delta_time:          1e-2      # Time between statistics output
+  energy_file_name:    energy    # (Optional) File name for energy output
+  timestep_file_name:  timesteps # (Optional) File name for timing information output. Note: No underscores "_" allowed in file name 
 
 # Parameters for the hydrodynamics scheme
 SPH:
diff --git a/examples/plot_scaling_results.py b/examples/plot_scaling_results.py
new file mode 100755
index 0000000000000000000000000000000000000000..66365bfad2edb38efa4b90c0c1602fd38bd750fd
--- /dev/null
+++ b/examples/plot_scaling_results.py
@@ -0,0 +1,207 @@
+#!/usr/bin/env python
+#
+# Usage:
+#  python plot_scaling_results.py input-file1-ext input-file2-ext ...
+#
+# Description:
+# Plots speed up, parallel efficiency and time to solution given a "timesteps" output file generated by SWIFT.
+# 
+# Example:
+# python plot_scaling_results.py _hreads_cosma_stdout.txt _threads_knl_stdout.txt
+# 
+# The working directory should contain files 1_threads_cosma_stdout.txt - 64_threads_cosma_stdout.txt and 1_threads_knl_stdout.txt - 64_threads_knl_stdout.txt, i.e wall clock time for each run using a given number of threads
+
+import sys
+import glob
+import re
+import numpy as np
+import matplotlib.pyplot as plt
+
+version = []
+branch = []
+revision = []
+hydro_scheme = []
+hydro_kernel = []
+hydro_neighbours = []
+hydro_eta = []
+threadList = []
+linestyle = ('ro-','bo-','go-','yo-','mo-')
+#cmdLine = './swift_fixdt -s -t 16 cosmoVolume.yml'
+#platform = 'KNL'
+
+# Work out how many data series there are
+if len(sys.argv) == 2:
+  inputFileNames = (sys.argv[1],"")
+  numOfSeries = 1
+elif len(sys.argv) == 3:
+  inputFileNames = (sys.argv[1],sys.argv[2])
+  numOfSeries = 2
+elif len(sys.argv) == 4:
+  inputFileNames = (sys.argv[1],sys.argv[2],sys.argv[3])
+  numOfSeries = 3
+elif len(sys.argv) == 5:
+  inputFileNames = (sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4])
+  numOfSeries = 4
+elif len(sys.argv) == 6:
+  inputFileNames = (sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4],sys.argv[5])
+  numOfSeries = 5
+
+# Get the names of the branch, Git revision, hydro scheme and hydro kernel
+def parse_header(inputFile):
+  with open(inputFile, 'r') as f:
+    found_end = False
+    for line in f:
+      if 'Branch:' in line:
+        s = line.split()
+        branch.append(s[2])
+      elif 'Revision:' in line:
+        s = line.split() 
+        revision.append(s[2])
+      elif 'Hydrodynamic scheme:' in line:
+        line = line[2:-1]
+        s = line.split()
+        line = s[2:]
+        hydro_scheme.append(" ".join(line))
+      elif 'Hydrodynamic kernel:' in line:
+        line = line[2:-1]
+        s = line.split()
+        line = s[2:5]
+        hydro_kernel.append(" ".join(line))
+      elif 'neighbours:' in line:
+        s = line.split() 
+        hydro_neighbours.append(s[4])
+      elif 'Eta:' in line:
+        s = line.split() 
+        hydro_eta.append(s[2])
+  return
+
+# Parse file and return total time taken, speed up and parallel efficiency
+def parse_files():
+  
+  times = []
+  totalTime = []
+  serialTime = []
+  speedUp = []
+  parallelEff = []
+
+  for i in range(0,numOfSeries): # Loop over each data series
+ 
+    # Get each file that starts with the cmd line arg
+    file_list = glob.glob(inputFileNames[i] + "*")
+    
+    threadList.append([])
+
+    # Create a list of threads using the list of files
+    for fileName in file_list:
+      s = re.split(r'[_.]+',fileName)
+      threadList[i].append(int(s[1]))
+
+    # Sort the thread list in ascending order and save the indices
+    sorted_indices = np.argsort(threadList[i])
+    threadList[i].sort()
+
+    # Sort the file list in ascending order acording to the thread number
+    file_list = [ file_list[j] for j in sorted_indices]
+
+    parse_header(file_list[0])
+    
+    version.append(branch[i] + " " + revision[i] + "\n" + hydro_scheme[i] + 
+                   "\n" + hydro_kernel[i] + r", $N_{neigh}$=" + hydro_neighbours[i] + 
+                   r", $\eta$=" + hydro_eta[i] + "\n")                  
+    times.append([])
+    totalTime.append([])
+    speedUp.append([])
+    parallelEff.append([])
+
+    # Loop over all files for a given series and load the times
+    for j in range(0,len(file_list)):
+      times[i].append([])
+      times[i][j].append(np.loadtxt(file_list[j],usecols=(5,)))
+      totalTime[i].append(np.sum(times[i][j]))
+
+    serialTime.append(totalTime[i][0])
+    
+    # Loop over all files for a given series and calculate speed up and parallel efficiency
+    for j in range(0,len(file_list)):
+      speedUp[i].append(serialTime[i] / totalTime[i][j])
+      parallelEff[i].append(speedUp[i][j] / threadList[i][j])
+
+  return (times,totalTime,speedUp,parallelEff)
+
+def print_results(times,totalTime,parallelEff,version):
+ 
+  for i in range(0,numOfSeries):
+    print " "
+    print "------------------------------------"
+    print version[i]
+    print "------------------------------------"
+    print "Wall clock time for: {} time steps".format(len(times[0][0][0]))
+    print "------------------------------------"
+    
+    for j in range(0,len(threadList[i])):
+      print str(threadList[i][j]) + " threads: {}".format(totalTime[i][j])
+    
+    print " "
+    print "------------------------------------"
+    print "Parallel Efficiency for: {} time steps".format(len(times[0][0][0]))
+    print "------------------------------------"
+    
+    for j in range(0,len(threadList[i])):
+      print str(threadList[i][j]) + " threads: {}".format(parallelEff[i][j])
+
+  return
+
+def plot_results(times,totalTime,speedUp,parallelEff):
+  
+  fig, axarr = plt.subplots(2, 2,figsize=(15,15))
+  speedUpPlot = axarr[0, 0]
+  parallelEffPlot = axarr[0, 1]
+  totalTimePlot = axarr[1, 0]
+  emptyPlot = axarr[1, 1]
+  
+  # Plot speed up
+  for i in range(0,numOfSeries):
+    speedUpPlot.plot(threadList[i],speedUp[i],linestyle[i],label=version[i])
+  
+  speedUpPlot.plot(threadList[i],threadList[i],color='k',linestyle='--')
+  speedUpPlot.set_ylabel("Speed Up")
+  speedUpPlot.set_xlabel("No. of Threads")
+
+  # Plot parallel efficiency
+  for i in range(0,numOfSeries):
+    parallelEffPlot.plot(threadList[i],parallelEff[i],linestyle[i])
+  
+  parallelEffPlot.set_xscale('log')
+  parallelEffPlot.set_ylabel("Parallel Efficiency")
+  parallelEffPlot.set_xlabel("No. of Threads")
+  parallelEffPlot.set_ylim([0,1.1])
+
+  # Plot time to solution     
+  for i in range(0,numOfSeries):
+    totalTimePlot.loglog(threadList[i],totalTime[i],linestyle[i],label=version[i])
+  
+  totalTimePlot.set_xscale('log')
+  totalTimePlot.set_xlabel("No. of Threads")
+  totalTimePlot.set_ylabel("Time to Solution (ms)")
+  
+  totalTimePlot.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.,prop={'size':14})
+  emptyPlot.axis('off')
+  
+  for i, txt in enumerate(threadList[0]):
+    speedUpPlot.annotate(txt, (threadList[0][i],speedUp[0][i]))
+    parallelEffPlot.annotate(txt, (threadList[0][i],parallelEff[0][i]))
+    totalTimePlot.annotate(txt, (threadList[0][i],totalTime[0][i]))
+
+  #fig.suptitle("Thread Speed Up, Parallel Efficiency and Time To Solution for {} Time Steps of Cosmo Volume\n Cmd Line: {}, Platform: {}".format(len(times[0][0][0]),cmdLine,platform))
+  fig.suptitle("Thread Speed Up, Parallel Efficiency and Time To Solution for {} Time Steps".format(len(times[0][0][0])))
+
+  return
+
+# Calculate results
+(times,totalTime,speedUp,parallelEff) = parse_files()
+
+plot_results(times,totalTime,speedUp,parallelEff)
+
+print_results(times,totalTime,parallelEff,version)
+
+plt.show()
diff --git a/src/engine.c b/src/engine.c
index 07332635f66f4acf4f0c4f99146b82cf0d84c5b2..2d7903f147f3948dd5590ad11811bce91c6d073a 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -65,6 +65,7 @@
 #include "single_io.h"
 #include "timers.h"
 #include "units.h"
+#include "version.h"
 
 const char *engine_policy_names[13] = {"none",
                                        "rand",
@@ -2392,6 +2393,10 @@ void engine_step(struct engine *e) {
     printf("  %6d %14e %14e %10zd %10zd %21.3f\n", e->step, e->time,
            e->timeStep, e->updates, e->g_updates, e->wallclock_time);
     fflush(stdout);
+
+    fprintf(e->file_timesteps, "  %6d %14e %14e %10zd %10zd %21.3f\n", e->step,
+            e->time, e->timeStep, e->updates, e->g_updates, e->wallclock_time);
+    fflush(e->file_timesteps);
   }
 
   /* Save some statistics */
@@ -2801,6 +2806,7 @@ void engine_init(struct engine *e, struct space *s,
   e->dt_min = parser_get_param_double(params, "TimeIntegration:dt_min");
   e->dt_max = parser_get_param_double(params, "TimeIntegration:dt_max");
   e->file_stats = NULL;
+  e->file_timesteps = NULL;
   e->deltaTimeStatistics =
       parser_get_param_double(params, "Statistics:delta_time");
   e->timeLastStatistics = e->timeBegin - e->deltaTimeStatistics;
@@ -2952,12 +2958,41 @@ void engine_init(struct engine *e, struct space *s,
 
   /* Open some files */
   if (e->nodeID == 0) {
-    e->file_stats = fopen("energy.txt", "w");
+    char energyfileName[200] = "";
+    parser_get_opt_param_string(params, "Statistics:energy_file_name",
+                                energyfileName,
+                                engine_default_energy_file_name);
+    sprintf(energyfileName + strlen(energyfileName), ".txt");
+    e->file_stats = fopen(energyfileName, "w");
     fprintf(e->file_stats,
             "# %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s\n",
             "Time", "Mass", "E_tot", "E_kin", "E_int", "E_pot", "p_x", "p_y",
             "p_z", "ang_x", "ang_y", "ang_z");
     fflush(e->file_stats);
+
+    char timestepsfileName[200] = "";
+    parser_get_opt_param_string(params, "Statistics:timestep_file_name",
+                                timestepsfileName,
+                                engine_default_timesteps_file_name);
+
+    sprintf(timestepsfileName + strlen(timestepsfileName), "_%d.txt",
+            nr_nodes * nr_threads);
+    e->file_timesteps = fopen(timestepsfileName, "w");
+    fprintf(e->file_timesteps,
+            "# Branch: %s\n# Revision: %s\n# Compiler: %s, Version: %s \n# "
+            "Number of threads: %d\n# Number of MPI ranks: %d\n# Hydrodynamic "
+            "scheme: %s\n# Hydrodynamic kernel: %s\n# No. of neighbours: %.2f "
+            "+/- %.2f\n# Eta: %f\n",
+            git_branch(), git_revision(), compiler_name(), compiler_version(),
+            e->nr_threads, e->nr_nodes, SPH_IMPLEMENTATION, kernel_name,
+            e->hydro_properties->target_neighbours,
+            e->hydro_properties->delta_neighbours,
+            e->hydro_properties->eta_neighbours);
+
+    fprintf(e->file_timesteps, "# %6s %14s %14s %10s %10s %16s [%s]\n", "Step",
+            "Time", "Time-step", "Updates", "g-Updates", "Wall-clock time",
+            clocks_getunit());
+    fflush(e->file_timesteps);
   }
 
   /* Print policy */
diff --git a/src/engine.h b/src/engine.h
index 3dfe22845a66ba02fcc270e13647efd4390c7256..8aa95b689e7cd8774cf8ddd57d46a4e838b9095a 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -72,6 +72,8 @@ extern const char *engine_policy_names[];
 #define engine_tasksreweight 10
 #define engine_parts_size_grow 1.05
 #define engine_redistribute_alloc_margin 1.2
+#define engine_default_energy_file_name "energy"
+#define engine_default_timesteps_file_name "timesteps"
 
 /* The rank of the engine as a global variable (for messages). */
 extern int engine_rank;
@@ -149,6 +151,9 @@ struct engine {
   double timeLastStatistics;
   double deltaTimeStatistics;
 
+  /* Timesteps information */
+  FILE *file_timesteps;
+
   /* The current step number. */
   int step;
 
diff --git a/tests/testReading.c b/tests/testReading.c
index 2fa88855a70a12265f180cd97528dda855322d1d..600c69c47b615e9b3f3a4aaa7a85056e74976d4f 100644
--- a/tests/testReading.c
+++ b/tests/testReading.c
@@ -27,6 +27,7 @@ int main() {
 
   size_t Ngas = 0, Ngpart = 0;
   int periodic = -1;
+  int flag_entropy_ICs = -1;
   int i, j, k, n;
   double dim[3];
   struct part *parts = NULL;
@@ -39,7 +40,7 @@ int main() {
 
   /* Read data */
   read_ic_single("input.hdf5", dim, &parts, &gparts, &Ngas, &Ngpart, &periodic,
-                 0);
+                 &flag_entropy_ICs, 0);
 
   /* Check global properties read are correct */
   assert(dim[0] == boxSize);