diff --git a/.gitignore b/.gitignore
index f054f6c825c3e34904e7b6161b49a00917755055..4c704f0665ed0b73d3abb0ca576bfd4709e080e2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -42,6 +42,7 @@ tests/testSingle
 tests/testTimeIntegration
 tests/testSPHStep
 tests/testKernel
+tests/testParser
 
 theory/latex/swift.pdf
 theory/kernel/kernels.pdf
diff --git a/examples/main.c b/examples/main.c
index f770265c1b59dff93513a59ff7246258539325e3..9523af49ed30c54d256d287ea2846a854650bc05 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -1,4 +1,3 @@
-
 /*******************************************************************************
  * This file is part of SWIFT.
  * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
@@ -56,7 +55,6 @@
  * @brief Main routine that loads a few particles and generates some output.
  *
  */
-
 int main(int argc, char *argv[]) {
 
   int c, icount, periodic = 1;
@@ -80,7 +78,8 @@ int main(int argc, char *argv[]) {
   int nr_nodes = 1, myrank = 0;
   FILE *file_thread;
   int with_outputs = 1;
-  int with_gravity = 0;
+  int with_external_gravity = 0;
+  int with_self_gravity = 0;
   int engine_policies = 0;
   int verbose = 0, talking = 0;
   unsigned long long cpufreq = 0;
@@ -159,7 +158,7 @@ int main(int argc, char *argv[]) {
   bzero(&s, sizeof(struct space));
 
   /* Parse the options */
-  while ((c = getopt(argc, argv, "a:c:d:e:f:gh:m:oP:q:R:s:t:v:w:y:z:")) != -1)
+  while ((c = getopt(argc, argv, "a:c:d:e:f:gGh:m:oP:q:R:s:t:v:w:y:z:")) != -1)
     switch (c) {
       case 'a':
         if (sscanf(optarg, "%lf", &scaling) != 1)
@@ -189,7 +188,10 @@ int main(int argc, char *argv[]) {
         if (!strcpy(ICfileName, optarg)) error("Error parsing IC file name.");
         break;
       case 'g':
-        with_gravity = 1;
+        with_external_gravity = 1;
+        break;
+      case 'G':
+        with_self_gravity = 1;
         break;
       case 'h':
         if (sscanf(optarg, "%llu", &cpufreq) != 1)
@@ -349,10 +351,6 @@ int main(int argc, char *argv[]) {
     message("CPU frequency used for tick conversion: %llu Hz", cpufreq);
   }
 
-  /* Check we have sensible time step bounds */
-  if (dt_min > dt_max)
-    error("Minimal time step size must be large than maximal time step size ");
-
   /* Check whether an IC file has been provided */
   if (strcmp(ICfileName, "") == 0)
     error("An IC file name must be provided via the option -f");
@@ -394,7 +392,7 @@ int main(int argc, char *argv[]) {
 #endif
 
   /* MATTHIEU: Temporary fix to preserve master */
-  if (!with_gravity) {
+  if (!with_external_gravity && !with_self_gravity) {
     free(gparts);
     gparts = NULL;
     for (size_t k = 0; k < Ngas; ++k) parts[k].gpart = NULL;
@@ -482,7 +480,8 @@ int main(int argc, char *argv[]) {
 
   /* Construct the engine policy */
   engine_policies = ENGINE_POLICY | engine_policy_steal | engine_policy_hydro;
-  if (with_gravity) engine_policies |= engine_policy_external_gravity;
+  if (with_external_gravity) engine_policies |= engine_policy_external_gravity;
+  if (with_self_gravity) engine_policies |= engine_policy_self_gravity;
 
   /* Initialize the engine with this space. */
   if (myrank == 0) clocks_gettime(&tic);
@@ -545,8 +544,8 @@ int main(int argc, char *argv[]) {
   /* Legend */
   if (myrank == 0)
     printf(
-        "# Step  Time  time-step  Number of updates    CPU Wall-clock time "
-        "[%s]\n",
+        "# Step  Time  time-step  Number of updates  Number of updates "
+        "CPU Wall-clock time [%s]\n",
         clocks_getunit());
 
   /* Let loose a runner on the space. */
diff --git a/src/Makefile.am b/src/Makefile.am
index 833f07e9e6cec94a7f1f3c01dfb0f975982ff931..a96f35b3cf0d8a23aec4f8c0f8d16bec8638cbcd 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -35,13 +35,13 @@ endif
 # List required headers
 include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     engine.h swift.h serial_io.h timers.h debug.h scheduler.h proxy.h parallel_io.h \
-    common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h
+    common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h parser.h
 
 # Common source files
 AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
     serial_io.c timers.c debug.c scheduler.c proxy.c parallel_io.c \
     units.c common_io.c single_io.c multipole.c version.c map.c \
-    kernel_hydro.c kernel_gravity.c tools.c part.c partition.c clocks.c
+    kernel_hydro.c kernel_gravity.c tools.c part.c partition.c clocks.c parser.c
 
 # Include files for distribution, not installation.
 nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
diff --git a/src/cell.c b/src/cell.c
index b0e299999e9073cc14d92b769b3e157f50d64c8f..61acfaaea7a0af01a78ab773541564e9a2723f4e 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -45,6 +45,7 @@
 /* Local headers. */
 #include "atomic.h"
 #include "error.h"
+#include "gravity.h"
 #include "hydro.h"
 #include "space.h"
 #include "timers.h"
@@ -605,6 +606,27 @@ void cell_init_parts(struct cell *c, void *data) {
   c->ti_end_max = 0;
 }
 
+/**
+ * @brief Initialises all g-particles to a valid state even if the ICs were
+ *stupid
+ *
+ * @param c Cell to act upon
+ * @param data Unused parameter
+ */
+void cell_init_gparts(struct cell *c, void *data) {
+
+  struct gpart *gp = c->gparts;
+  const int gcount = c->gcount;
+
+  for (int i = 0; i < gcount; ++i) {
+    gp[i].ti_begin = 0;
+    gp[i].ti_end = 0;
+    gravity_first_init_gpart(&gp[i]);
+  }
+  c->ti_end_min = 0;
+  c->ti_end_max = 0;
+}
+
 /**
  * @brief Converts hydro quantities to a valid state after the initial density
  *calculation
diff --git a/src/cell.h b/src/cell.h
index 05b870b144dd425ee9e72ee3b529122b6e721523..a471eac44bfd3533c4220ab8c5ff2ddec724e87f 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -141,7 +141,7 @@ struct cell {
   double mass, e_pot, e_int, e_kin;
 
   /* Number of particles updated in this cell. */
-  int updated;
+  int updated, g_updated;
 
   /* Linking pointer for "memory management". */
   struct cell *next;
@@ -178,6 +178,7 @@ int cell_getsize(struct cell *c);
 int cell_link_parts(struct cell *c, struct part *parts);
 int cell_link_gparts(struct cell *c, struct gpart *gparts);
 void cell_init_parts(struct cell *c, void *data);
+void cell_init_gparts(struct cell *c, void *data);
 void cell_convert_hydro(struct cell *c, void *data);
 void cell_clean_links(struct cell *c, void *data);
 
diff --git a/src/engine.c b/src/engine.c
index 28d75bba5378b83d4d0c8f6bf5c7705b40f6edb8..9e854554806bac25d0c20268991a05b59b16161e 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -442,10 +442,10 @@ void engine_redistribute(struct engine *e) {
       if (gparts_new[k].part->gpart != &gparts_new[k])
         error("Linking problem !");
 
-      /* if (gparts_new[k].x[0] != gparts_new[k].part->x[0] || */
-      /*     gparts_new[k].x[1] != gparts_new[k].part->x[1] || */
-      /*     gparts_new[k].x[2] != gparts_new[k].part->x[2]) */
-      /*   error("Linked particles are not at the same position !"); */
+      if (gparts_new[k].x[0] != gparts_new[k].part->x[0] ||
+          gparts_new[k].x[1] != gparts_new[k].part->x[1] ||
+          gparts_new[k].x[2] != gparts_new[k].part->x[2])
+        error("Linked particles are not at the same position !");
     }
   }
   for (size_t k = 0; k < nr_parts; ++k) {
@@ -1558,6 +1558,7 @@ int engine_marktasks(struct engine *e) {
       else if (t->type == task_type_kick) {
         t->skip = (t->ci->ti_end_min > ti_end);
         t->ci->updated = 0;
+        t->ci->g_updated = 0;
       }
 
       /* Drift? */
@@ -1614,6 +1615,7 @@ void engine_print_task_counts(struct engine *e) {
   printf(" skipped=%i ]\n", counts[task_type_count]);
   fflush(stdout);
   message("nr_parts = %zi.", e->s->nr_parts);
+  message("nr_gparts = %zi.", e->s->nr_gparts);
 }
 
 /**
@@ -1743,7 +1745,7 @@ void engine_collect_kick(struct cell *c) {
   if (c->kick != NULL) return;
 
   /* Counters for the different quantities. */
-  int updated = 0;
+  int updated = 0, g_updated = 0;
   double e_kin = 0.0, e_int = 0.0, e_pot = 0.0;
   float mom[3] = {0.0f, 0.0f, 0.0f}, ang[3] = {0.0f, 0.0f, 0.0f};
   int ti_end_min = max_nr_timesteps, ti_end_max = 0;
@@ -1766,6 +1768,7 @@ void engine_collect_kick(struct cell *c) {
         ti_end_min = min(ti_end_min, cp->ti_end_min);
         ti_end_max = max(ti_end_max, cp->ti_end_max);
         updated += cp->updated;
+        g_updated += cp->g_updated;
         e_kin += cp->e_kin;
         e_int += cp->e_int;
         e_pot += cp->e_pot;
@@ -1783,6 +1786,7 @@ void engine_collect_kick(struct cell *c) {
   c->ti_end_min = ti_end_min;
   c->ti_end_max = ti_end_max;
   c->updated = updated;
+  c->g_updated = g_updated;
   c->e_kin = e_kin;
   c->e_int = e_int;
   c->e_pot = e_pot;
@@ -1846,7 +1850,15 @@ void engine_init_particles(struct engine *e) {
 
   /* Make sure all particles are ready to go */
   /* i.e. clean-up any stupid state in the ICs */
-  space_map_cells_pre(s, 1, cell_init_parts, NULL);
+  if ((e->policy & engine_policy_hydro) == engine_policy_hydro) {
+    space_map_cells_pre(s, 1, cell_init_parts, NULL);
+  }
+  if (((e->policy & engine_policy_self_gravity) ==
+       engine_policy_self_gravity) ||
+      ((e->policy & engine_policy_external_gravity) ==
+       engine_policy_external_gravity)) {
+    space_map_cells_pre(s, 1, cell_init_gparts, NULL);
+  }
 
   engine_prepare(e);
 
@@ -1920,7 +1932,7 @@ void engine_init_particles(struct engine *e) {
  */
 void engine_step(struct engine *e) {
 
-  int updates = 0;
+  int updates = 0, g_updates = 0;
   int ti_end_min = max_nr_timesteps, ti_end_max = 0;
   double e_pot = 0.0, e_int = 0.0, e_kin = 0.0;
   float mom[3] = {0.0, 0.0, 0.0};
@@ -1947,6 +1959,7 @@ void engine_step(struct engine *e) {
       e_int += c->e_int;
       e_pot += c->e_pot;
       updates += c->updated;
+      g_updates += c->g_updated;
       mom[0] += c->mom[0];
       mom[1] += c->mom[1];
       mom[2] += c->mom[2];
@@ -1972,18 +1985,20 @@ void engine_step(struct engine *e) {
     ti_end_max = in_i[0];
   }
   {
-    double in_d[4], out_d[4];
+    double in_d[5], out_d[5];
     out_d[0] = updates;
-    out_d[1] = e_kin;
-    out_d[2] = e_int;
-    out_d[3] = e_pot;
-    if (MPI_Allreduce(out_d, in_d, 4, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) !=
+    out_d[1] = g_updates;
+    out_d[2] = e_kin;
+    out_d[3] = e_int;
+    out_d[4] = e_pot;
+    if (MPI_Allreduce(out_d, in_d, 5, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) !=
         MPI_SUCCESS)
       error("Failed to aggregate energies.");
     updates = in_d[0];
-    e_kin = in_d[1];
-    e_int = in_d[2];
-    e_pot = in_d[3];
+    g_updates = in_d[1];
+    e_kin = in_d[2];
+    e_int = in_d[3];
+    e_pot = in_d[4];
   }
 #endif
 
@@ -2008,8 +2023,8 @@ void engine_step(struct engine *e) {
   if (e->nodeID == 0) {
 
     /* Print some information to the screen */
-    printf("%d %e %e %d %.3f\n", e->step, e->time, e->timeStep, updates,
-           e->wallclock_time);
+    printf("%d %e %e %d %d %.3f\n", e->step, e->time, e->timeStep, updates,
+           g_updates, e->wallclock_time);
     fflush(stdout);
 
     /* Write some energy statistics */
@@ -2251,6 +2266,28 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
   for (size_t k = 0; k < s->nr_gparts; k++)
     if (s->gparts[k].id > 0) s->gparts[k].part->gpart = &s->gparts[k];
 
+  /* Verify that the links are correct */
+  /* MATTHIEU: To be commented out once we are happy */
+  for (size_t k = 0; k < s->nr_gparts; ++k) {
+
+    if (s->gparts[k].id > 0) {
+
+      if (s->gparts[k].part->gpart != &s->gparts[k]) error("Linking problem !");
+
+      if (s->gparts[k].x[0] != s->gparts[k].part->x[0] ||
+          s->gparts[k].x[1] != s->gparts[k].part->x[1] ||
+          s->gparts[k].x[2] != s->gparts[k].part->x[2])
+        error("Linked particles are not at the same position !");
+    }
+  }
+  for (size_t k = 0; k < s->nr_parts; ++k) {
+
+    if (s->parts[k].gpart != NULL) {
+
+      if (s->parts[k].gpart->part != &s->parts[k]) error("Linking problem !");
+    }
+  }
+
 #else
   error("SWIFT was not compiled with MPI support.");
 #endif
@@ -2426,6 +2463,7 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
 
   /* Print information about the hydro scheme */
   if (e->nodeID == 0) message("Hydrodynamic scheme: %s", SPH_IMPLEMENTATION);
+  if (e->nodeID == 0) message("Hydrodynamic kernel: %s", kernel_name);
 
   /* Check we have sensible time bounds */
   if (timeBegin >= timeEnd)
@@ -2434,10 +2472,12 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
         "(t_beg = %e)",
         timeEnd, timeBegin);
 
-  /* Check we have sensible time step bounds */
+  /* Check we have sensible time-step values */
   if (e->dt_min > e->dt_max)
     error(
-        "Minimal time step size must be smaller than maximal time step size ");
+        "Minimal time-step size (%e) must be smaller than maximal time-step "
+        "size (%e)",
+        e->dt_min, e->dt_max);
 
   /* Deal with timestep */
   e->timeBase = (timeEnd - timeBegin) / max_nr_timesteps;
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 82bc52ad3e05794c8c05896075edc463a69197ff..92a9f64c1f84a9e949f4c0e9485f892b5c808cdc 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -22,14 +22,61 @@
 /**
  * @brief Computes the gravity time-step of a given particle
  *
- * @param p Pointer to the particle data
- * @param xp Pointer to the extended particle data
+ * @param gp Pointer to the g-particle data
  *
  */
 
-__attribute__((always_inline)) INLINE static float gravity_compute_timestep(
-    struct part* p, struct xpart* xp) {
+__attribute__((always_inline))
+    INLINE static float gravity_compute_timestep(struct gpart* gp) {
 
   /* Currently no limit is imposed */
   return FLT_MAX;
 }
+
+/**
+ * @brief Initialises the g-particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param gp The particle to act upon
+ */
+__attribute__((always_inline))
+    INLINE static void gravity_first_init_gpart(struct gpart* gp) {}
+
+/**
+ * @brief Prepares a g-particle for the gravity calculation
+ *
+ * Zeroes all the relevant arrays in preparation for the sums taking place in
+ * the variaous tasks
+ *
+ * @param gp The particle to act upon
+ */
+__attribute__((always_inline))
+    INLINE static void gravity_init_part(struct gpart* gp) {
+
+  /* Zero the acceleration */
+  gp->a_grav[0] = 0.f;
+  gp->a_grav[1] = 0.f;
+  gp->a_grav[2] = 0.f;
+}
+
+/**
+ * @brief Finishes the gravity calculation.
+ *
+ * Multiplies the forces and accelerations by the appropiate constants
+ *
+ * @param gp The particle to act upon
+ */
+__attribute__((always_inline))
+    INLINE static void gravity_end_force(struct gpart* gp) {}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * @param gp The particle to act upon
+ * @param dt The time-step for this kick
+ * @param half_dt The half time-step for this kick
+ */
+__attribute__((always_inline)) INLINE static void gravity_kick_extra(
+    struct gpart* gp, float dt, float half_dt) {}
diff --git a/src/gravity/Default/gravity_debug.h b/src/gravity/Default/gravity_debug.h
index 98e0c40a5700b4da70f27fb0955592bb5d2287c3..654745bfeb70dddba772af9e23797713376377a7 100644
--- a/src/gravity/Default/gravity_debug.h
+++ b/src/gravity/Default/gravity_debug.h
@@ -24,5 +24,5 @@ __attribute__((always_inline))
       "v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
       "mass=%.3e t_begin=%d, t_end=%d\n",
       p->x[0], p->x[1], p->x[2], p->v_full[0], p->v_full[1], p->v_full[2],
-      p->a[0], p->a[1], p->a[2], p->mass, p->ti_begin, p->ti_end);
+      p->a_grav[0], p->a_grav[1], p->a_grav[2], p->mass, p->ti_begin, p->ti_end);
 }
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 9158bf6ec288e78b5d3b73640a29520f3a00678b..62023345f174eb8cb9bae4d4438bdd50c9969494 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -49,8 +49,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav(
   /* Aggregate the accelerations. */
   for (k = 0; k < 3; k++) {
     w = acc * dx[k];
-    pi->a[k] -= w * mj;
-    pj->a[k] += w * mi;
+    pi->a_grav[k] -= w * mj;
+    pj->a_grav[k] += w * mi;
   }
 }
 
@@ -100,8 +100,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_grav(
     ai.v = w.v * mj.v;
     aj.v = w.v * mi.v;
     for (j = 0; j < VEC_SIZE; j++) {
-      pi[j]->a[k] -= ai.f[j];
-      pj[j]->a[k] += aj.f[j];
+      pi[j]->a_grav[k] -= ai.f[j];
+      pj[j]->a_grav[k] += aj.f[j];
     }
   }
 
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index aaaab2055d8ce05d77a37d3f40e096bd91d87926..0dfdb82e4ec11c9153f77439027d7e4451ded7f4 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -29,7 +29,7 @@ struct gpart {
   float v_full[3];
 
   /* Particle acceleration. */
-  float a[3];
+  float a_grav[3];
 
   /* Particle mass. */
   float mass;
diff --git a/src/multipole.h b/src/multipole.h
index d5939878d86ff7f7540870942f00e7fd7da7f238..85ba44d3ce95d958b721d435ccd26b72e30a79c1 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -127,7 +127,7 @@ __attribute__((always_inline)) INLINE static void multipole_iact_mp(
 
 /* Compute the forces on both multipoles. */
 #if multipole_order == 1
-  for (k = 0; k < 3; k++) p->a[k] += dx[k] * acc;
+  for (k = 0; k < 3; k++) p->a_grav[k] += dx[k] * acc;
 #else
 #error( "Multipoles of order %i not yet implemented." , multipole_order )
 #endif
diff --git a/src/parser.c b/src/parser.c
new file mode 100644
index 0000000000000000000000000000000000000000..06dc819842d54d952704e4e0c40ebec5b561f691
--- /dev/null
+++ b/src/parser.c
@@ -0,0 +1,265 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 James Willis (james.s.willis@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+/* Needs to be included so that strtok returns char * instead of a int *. */
+#include <string.h>
+#include <stdlib.h>
+
+/* This object's header. */
+#include "parser.h"
+
+/* Local headers. */
+#include "error.h"
+
+/* Private functions. */
+static int count_char(char *str, char val);
+static void parse_line(FILE *fp, struct swift_params *params);
+
+/**
+ * @brief Reads an input file and stores each parameter in a structure.
+ *
+ * @param file_name Name of file to be read
+ * @param params Structure to be populated from file
+ */
+
+void parser_read_file(const char *file_name, struct swift_params *params) {
+
+  FILE *fp;
+
+  params->count = 0;
+
+  /* Open file for reading */
+  fp = fopen(file_name, "r");
+
+  if (fp == NULL) {
+    error("Error opening parameter file: %s", file_name);
+  }
+
+  /* Read until the end of the file is reached.*/
+  while (!feof(fp)) {
+    parse_line(fp, params);
+  }
+
+  fclose(fp);
+}
+
+/**
+ * @brief Counts the number of times a specific character appears in a string.
+ *
+ * @param str String to be checked
+ * @param val Character to be counted
+ */
+
+static int count_char(char *str, char val) {
+
+  int count = 0;
+
+  /* Check if the line contains the character */
+  while (*str) {
+    if (*str++ == val) ++count;
+  }
+
+  return count;
+}
+
+/**
+ * @brief Parses a line from a file and stores any parameters in a structure.
+ *
+ * @param fp File pointer to file to be read
+ * @param params Structure to be populated from file
+ *
+ */
+
+static void parse_line(FILE *fp, struct swift_params *params) {
+
+  char line[PARSER_MAX_LINE_SIZE];
+  char trim_line[PARSER_MAX_LINE_SIZE];
+
+  /* Read a line of the file */
+  if (fgets(line, PARSER_MAX_LINE_SIZE, fp) != NULL) {
+
+    char *token;
+    /* Remove comments */
+    token = strtok(line, PARSER_COMMENT_CHAR);
+    strcpy(trim_line, token);
+
+    /* Check if the line contains a value */
+    if (strchr(trim_line, PARSER_VALUE_CHAR)) {
+      /* Check for more than one parameter on the same line. */
+      if (count_char(trim_line, PARSER_VALUE_CHAR) > 1) {
+        error("Found more than one parameter in '%s', only one allowed.", line);
+      } else {
+        /* Take first token as the parameter name. */
+        token = strtok(trim_line, PARSER_VALUE_STRING);
+        strcpy(params->data[params->count].name, token);
+
+        /* Take second token as the parameter value. */
+        token = strtok(NULL, " #\n");
+        strcpy(params->data[params->count++].value, token);
+      }
+    }
+  }
+}
+
+/**
+ * @brief Retrieve integer parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param retParam Value of the parameter found
+ *
+ */
+
+void parser_get_param_int(struct swift_params *params, char *name,
+                          int *retParam) {
+
+  char str[128];
+
+  for (int i = 0; i < params->count; i++) {
+
+    /*strcmp returns 0 if both strings are the same.*/
+    if (!strcmp(name, params->data[i].name)) {
+
+      /* Check that exactly one number is parsed. */
+      if (sscanf(params->data[i].value, "%d%s", retParam, str) != 1) {
+        error(
+            "Tried parsing int '%s' but found '%s' with illegal integer "
+            "characters '%s'.",
+            params->data[i].name, params->data[i].value, str);
+      }
+
+      return;
+    }
+  }
+
+  message("Cannot find '%s' in the structure.", name);
+}
+
+/**
+ * @brief Retrieve float parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param retParam Value of the parameter found
+ *
+ */
+
+void parser_get_param_float(struct swift_params *params, char *name,
+                            float *retParam) {
+
+  char str[128];
+
+  for (int i = 0; i < params->count; i++) {
+
+    /*strcmp returns 0 if both strings are the same.*/
+    if (!strcmp(name, params->data[i].name)) {
+
+      /* Check that exactly one number is parsed. */
+      if (sscanf(params->data[i].value, "%f%s", retParam, str) != 1) {
+        error(
+            "Tried parsing float '%s' but found '%s' with illegal float "
+            "characters '%s'.",
+            params->data[i].name, params->data[i].value, str);
+      }
+
+      return;
+    }
+  }
+
+  message("Cannot find '%s' in the structure.", name);
+}
+
+/**
+ * @brief Retrieve double parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param retParam Value of the parameter found
+ *
+ */
+
+void parser_get_param_double(struct swift_params *params, char *name,
+                             double *retParam) {
+
+  char str[128];
+
+  for (int i = 0; i < params->count; i++) {
+
+    /*strcmp returns 0 if both strings are the same.*/
+    if (!strcmp(name, params->data[i].name)) {
+
+      /* Check that exactly one number is parsed. */
+      if (sscanf(params->data[i].value, "%lf", retParam) != 1) {
+        error(
+            "Tried parsing double '%s' but found '%s' with illegal double "
+            "characters '%s'.",
+            params->data[i].name, params->data[i].value, str);
+      }
+
+      return;
+    }
+  }
+
+  message("Cannot find '%s' in the structure.", name);
+}
+
+/**
+ * @brief Retrieve string parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param retParam Value of the parameter found
+ *
+ */
+
+void parser_get_param_string(struct swift_params *params, char *name,
+                             char *retParam) {
+
+  for (int i = 0; i < params->count; i++) {
+
+    /*strcmp returns 0 if both strings are the same.*/
+    if (!strcmp(name, params->data[i].name)) {
+      strcpy(retParam, params->data[i].value);
+      return;
+    }
+  }
+}
+
+/**
+ * @brief Prints the contents of the parameter structure.
+ *
+ * @param params Structure that holds the parameters
+ *
+ */
+
+void parser_print_params(struct swift_params *params) {
+
+  printf("\n--------------------------\n");
+  printf("|  SWIFT Parameter File  |\n");
+  printf("--------------------------\n");
+
+  for (int i = 0; i < params->count; i++) {
+    printf("Parameter name: %s\n", params->data[i].name);
+    printf("Parameter value: %s\n", params->data[i].value);
+  }
+}
diff --git a/src/parser.h b/src/parser.h
new file mode 100644
index 0000000000000000000000000000000000000000..2fb4148944cd423da016341744cb6d58e222182e
--- /dev/null
+++ b/src/parser.h
@@ -0,0 +1,54 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 James Willis (james.s.willis@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_PARSER_H
+#define SWIFT_PARSER_H
+
+#include <stdio.h>
+
+#define PARSER_MAX_LINE_SIZE 128
+#define PARSER_MAX_NO_OF_PARAMS 512
+
+#define PARSER_COMMENT_CHAR "#"
+#define PARSER_VALUE_CHAR ':'
+#define PARSER_VALUE_STRING ":"
+#define PARSER_END_OF_FILE "..."
+
+struct parameter {
+  char name[PARSER_MAX_LINE_SIZE];
+  char value[PARSER_MAX_LINE_SIZE];
+};
+
+struct swift_params {
+  struct parameter data[PARSER_MAX_NO_OF_PARAMS];
+  int count;
+};
+
+/* Public API. */
+void parser_read_file(const char *file_name, struct swift_params *params);
+void parser_print_params(struct swift_params *params);
+void parser_get_param_int(struct swift_params *params, char *name,
+                          int *retParam);
+void parser_get_param_float(struct swift_params *params, char *name,
+                            float *retParam);
+void parser_get_param_double(struct swift_params *params, char *name,
+                             double *retParam);
+void parser_get_param_string(struct swift_params *params, char *name,
+                             char *retParam);
+
+#endif /* SWIFT_PARSER_H */
diff --git a/src/runner.c b/src/runner.c
index 5e1b6497d757c8c25230823e7f8d60668f947618..d382e4802d3f1a5afd22761e97fea9f5128e02b1 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -469,8 +469,10 @@ void runner_dogsort(struct runner *r, struct cell *c, int flags, int clock) {
 
 void runner_doinit(struct runner *r, struct cell *c, int timer) {
 
-  struct part *p, *parts = c->parts;
+  struct part *const parts = c->parts;
+  struct gpart *const gparts = c->gparts;
   const int count = c->count;
+  const int gcount = c->gcount;
   const int ti_current = r->e->ti_current;
 
   TIMER_TIC;
@@ -486,7 +488,7 @@ void runner_doinit(struct runner *r, struct cell *c, int timer) {
     for (int i = 0; i < count; i++) {
 
       /* Get a direct pointer on the part. */
-      p = &parts[i];
+      struct part *const p = &parts[i];
 
       if (p->ti_end <= ti_current) {
 
@@ -494,6 +496,19 @@ void runner_doinit(struct runner *r, struct cell *c, int timer) {
         hydro_init_part(p);
       }
     }
+
+    /* Loop over the gparts in this cell. */
+    for (int i = 0; i < gcount; i++) {
+
+      /* Get a direct pointer on the part. */
+      struct gpart *const gp = &gparts[i];
+
+      if (gp->ti_end <= ti_current) {
+
+        /* Get ready for a density calculation */
+        gravity_init_part(gp);
+      }
+    }
   }
 
   if (timer) TIMER_TOC(timer_init);
@@ -649,7 +664,7 @@ void runner_doghost(struct runner *r, struct cell *c) {
 }
 
 /**
- * @brief Drift particles forward in time
+ * @brief Drift particles and g-particles forward in time
  *
  * @param r The runner thread.
  * @param c The cell.
@@ -658,26 +673,39 @@ void runner_doghost(struct runner *r, struct cell *c) {
 void runner_dodrift(struct runner *r, struct cell *c, int timer) {
 
   const int nr_parts = c->count;
+  const int nr_gparts = c->gcount;
   const double timeBase = r->e->timeBase;
   const double dt = (r->e->ti_current - r->e->ti_old) * timeBase;
-  const float ti_old = r->e->ti_old;
-  const float ti_current = r->e->ti_current;
-  struct part *restrict p, *restrict parts = c->parts;
-  struct xpart *restrict xp, *restrict xparts = c->xparts;
+  const int ti_old = r->e->ti_old;
+  const int ti_current = r->e->ti_current;
+  struct part *const parts = c->parts;
+  struct xpart *const xparts = c->xparts;
+  struct gpart *const gparts = c->gparts;
   float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
-  float w;
 
   TIMER_TIC
 
   /* No children? */
   if (!c->split) {
 
-    /* Loop over all the particles in the cell */
+    /* Loop over all the g-particles in the cell */
+    for (int k = 0; k < nr_gparts; ++k) {
+
+      /* Get a handle on the gpart. */
+      struct gpart *const gp = &gparts[k];
+
+      /* Drift... */
+      gp->x[0] += gp->v_full[0] * dt;
+      gp->x[1] += gp->v_full[1] * dt;
+      gp->x[2] += gp->v_full[2] * dt;
+    }
+
+    /* Loop over all the particles in the cell (more work for these !) */
     for (int k = 0; k < nr_parts; k++) {
 
       /* Get a handle on the part. */
-      p = &parts[k];
-      xp = &xparts[k];
+      struct part *const p = &parts[k];
+      struct xpart *const xp = &xparts[k];
 
       /* Useful quantity */
       const float h_inv = 1.0f / p->h;
@@ -693,18 +721,18 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
       p->v[2] += p->a_hydro[2] * dt;
 
       /* Predict smoothing length */
-      w = p->h_dt * h_inv * dt;
-      if (fabsf(w) < 0.2f)
-        p->h *= approx_expf(w); /* 4th order expansion of exp(w) */
+      const float w1 = p->h_dt * h_inv * dt;
+      if (fabsf(w1) < 0.2f)
+        p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
       else
-        p->h *= expf(w);
+        p->h *= expf(w1);
 
       /* Predict density */
-      w = -3.0f * p->h_dt * h_inv * dt;
-      if (fabsf(w) < 0.2f)
-        p->rho *= approx_expf(w); /* 4th order expansion of exp(w) */
+      const float w2 = -3.0f * p->h_dt * h_inv * dt;
+      if (fabsf(w2) < 0.2f)
+        p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
       else
-        p->rho *= expf(w);
+        p->rho *= expf(w2);
 
       /* Predict the values of the extra fields */
       hydro_predict_extra(p, xp, ti_old, ti_current, timeBase);
@@ -760,37 +788,97 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
   const double timeBase = r->e->timeBase;
   const double timeBase_inv = 1.0 / r->e->timeBase;
   const int count = c->count;
+  const int gcount = c->gcount;
+  struct part *const parts = c->parts;
+  struct xpart *const xparts = c->xparts;
+  struct gpart *const gparts = c->gparts;
   const int is_fixdt =
       (r->e->policy & engine_policy_fixdt) == engine_policy_fixdt;
 
-  int new_dti;
-  int dti_timeline;
-
-  int updated = 0;
+  int updated = 0, g_updated = 0;
   int ti_end_min = max_nr_timesteps, ti_end_max = 0;
   double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0;
   float mom[3] = {0.0f, 0.0f, 0.0f};
   float ang[3] = {0.0f, 0.0f, 0.0f};
-  float x[3], v_full[3];
-  struct part *restrict p, *restrict parts = c->parts;
-  struct xpart *restrict xp, *restrict xparts = c->xparts;
 
   TIMER_TIC
 
   /* No children? */
   if (!c->split) {
 
+    /* Loop over the g-particles and kick the active ones. */
+    for (int k = 0; k < gcount; k++) {
+
+      /* Get a handle on the part. */
+      struct gpart *const gp = &gparts[k];
+
+      /* If the g-particle has no counterpart and needs to be kicked */
+      if (gp->id < 0 && (is_fixdt || gp->ti_end <= ti_current)) {
+
+        /* First, finish the force calculation */
+        gravity_end_force(gp);
+
+        /* Now we are ready to compute the next time-step size */
+        int new_dti;
+
+        if (is_fixdt) {
+
+          /* Now we have a time step, proceed with the kick */
+          new_dti = global_dt_max * timeBase_inv;
+
+        } else {
+
+          /* Compute the next timestep (gravity condition) */
+          float new_dt = gravity_compute_timestep(gp);
+
+          /* Limit timestep within the allowed range */
+          new_dt = fminf(new_dt, global_dt_max);
+          new_dt = fmaxf(new_dt, global_dt_min);
+
+          /* Convert to integer time */
+          new_dti = new_dt * timeBase_inv;
+
+          /* Recover the current timestep */
+          const int current_dti = gp->ti_end - gp->ti_begin;
+
+          /* Limit timestep increase */
+          if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
+
+          /* Put this timestep on the time line */
+          int dti_timeline = max_nr_timesteps;
+          while (new_dti < dti_timeline) dti_timeline /= 2;
+
+          /* Now we have a time step, proceed with the kick */
+          new_dti = dti_timeline;
+        }
+
+        /* Compute the time step for this kick */
+        const int ti_start = (gp->ti_begin + gp->ti_end) / 2;
+        const int ti_end = gp->ti_end + new_dti / 2;
+        const double dt = (ti_end - ti_start) * timeBase;
+        const double half_dt = (ti_end - gp->ti_end) * timeBase;
+
+        /* Kick particles in momentum space */
+        gp->v_full[0] += gp->a_grav[0] * dt;
+        gp->v_full[1] += gp->a_grav[1] * dt;
+        gp->v_full[2] += gp->a_grav[2] * dt;
+
+        /* Extra kick work */
+        gravity_kick_extra(gp, dt, half_dt);
+
+        /* Number of updated g-particles */
+        g_updated++;
+      }
+    }
+
+    /* Now do the hydro ones... */
+
     /* Loop over the particles and kick the active ones. */
     for (int k = 0; k < count; k++) {
 
       /* Get a handle on the part. */
-      p = &parts[k];
-      xp = &xparts[k];
-
-      const float m = p->mass;
-      x[0] = p->x[0];
-      x[1] = p->x[1];
-      x[2] = p->x[2];
+      struct part *const p = &parts[k];
+      struct xpart *const xp = &xparts[k];
 
       /* If particle needs to be kicked */
       if (is_fixdt || p->ti_end <= ti_current) {
@@ -800,8 +888,10 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
 
         /* And do the same of the extra variable */
         hydro_end_force(p);
+        if (p->gpart != NULL) gravity_end_force(p->gpart);
 
         /* Now we are ready to compute the next time-step size */
+        int new_dti;
 
         if (is_fixdt) {
 
@@ -810,9 +900,13 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
 
         } else {
 
-          /* Compute the next timestep */
+          /* Compute the next timestep (hydro condition) */
           const float new_dt_hydro = hydro_compute_timestep(p, xp);
-          const float new_dt_grav = gravity_compute_timestep(p, xp);
+
+          /* Compute the next timestep (gravity condition) */
+          float new_dt_grav = FLT_MAX;
+          if (p->gpart != NULL)
+            new_dt_grav = gravity_compute_timestep(p->gpart);
 
           float new_dt = fminf(new_dt_hydro, new_dt_grav);
 
@@ -837,7 +931,7 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
           if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
 
           /* Put this timestep on the time line */
-          dti_timeline = max_nr_timesteps;
+          int dti_timeline = max_nr_timesteps;
           while (new_dti < dti_timeline) dti_timeline /= 2;
 
           /* Now we have a time step, proceed with the kick */
@@ -847,34 +941,51 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
         /* Compute the time step for this kick */
         const int ti_start = (p->ti_begin + p->ti_end) / 2;
         const int ti_end = p->ti_end + new_dti / 2;
-        const float dt = (ti_end - ti_start) * timeBase;
-        const float half_dt = (ti_end - p->ti_end) * timeBase;
+        const double dt = (ti_end - ti_start) * timeBase;
+        const double half_dt = (ti_end - p->ti_end) * timeBase;
 
         /* Move particle forward in time */
         p->ti_begin = p->ti_end;
         p->ti_end = p->ti_begin + new_dti;
 
+        /* Get the acceleration */
+        float a_tot[3] = {p->a_hydro[0], p->a_hydro[1], p->a_hydro[2]};
+        if (p->gpart != NULL) {
+          a_tot[0] += p->gpart->a_grav[0];
+          a_tot[1] += p->gpart->a_grav[1];
+          a_tot[1] += p->gpart->a_grav[2];
+        }
+
         /* Kick particles in momentum space */
-        xp->v_full[0] += p->a_hydro[0] * dt;
-        xp->v_full[1] += p->a_hydro[1] * dt;
-        xp->v_full[2] += p->a_hydro[2] * dt;
+        xp->v_full[0] += a_tot[0] * dt;
+        xp->v_full[1] += a_tot[1] * dt;
+        xp->v_full[2] += a_tot[2] * dt;
+
+        if (p->gpart != NULL) {
+          p->gpart->v_full[0] = xp->v_full[0];
+          p->gpart->v_full[1] = xp->v_full[1];
+          p->gpart->v_full[2] = xp->v_full[2];
+        }
 
-        p->v[0] = xp->v_full[0] - half_dt * p->a_hydro[0];
-        p->v[1] = xp->v_full[1] - half_dt * p->a_hydro[1];
-        p->v[2] = xp->v_full[2] - half_dt * p->a_hydro[2];
+        /* Go back by half-step for the hydro velocity */
+        p->v[0] = xp->v_full[0] - half_dt * a_tot[0];
+        p->v[1] = xp->v_full[1] - half_dt * a_tot[1];
+        p->v[2] = xp->v_full[2] - half_dt * a_tot[2];
 
         /* Extra kick work */
         hydro_kick_extra(p, xp, dt, half_dt);
+        if (p->gpart != NULL) gravity_kick_extra(p->gpart, dt, half_dt);
 
         /* Number of updated particles */
         updated++;
+        if (p->gpart != NULL) g_updated++;
       }
 
       /* Now collect quantities for statistics */
 
-      v_full[0] = xp->v_full[0];
-      v_full[1] = xp->v_full[1];
-      v_full[2] = xp->v_full[2];
+      const double x[3] = {p->x[0], p->x[1], p->x[2]};
+      const float v_full[3] = {xp->v_full[0], xp->v_full[1], xp->v_full[2]};
+      const float m = p->mass;
 
       /* Collect mass */
       mass += m;
@@ -908,13 +1019,14 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
     /* Loop over the progeny. */
     for (int k = 0; k < 8; k++)
       if (c->progeny[k] != NULL) {
-        struct cell *cp = c->progeny[k];
+        struct cell *const cp = c->progeny[k];
 
         /* Recurse */
         runner_dokick(r, cp, 0);
 
         /* And aggregate */
         updated += cp->updated;
+        g_updated += cp->g_updated;
         e_kin += cp->e_kin;
         e_int += cp->e_int;
         e_pot += cp->e_pot;
@@ -932,6 +1044,7 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
 
   /* Store the values. */
   c->updated = updated;
+  c->g_updated = g_updated;
   c->e_kin = e_kin;
   c->e_int = e_int;
   c->e_pot = e_pot;
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index f374339da75e31b39a5295fcd8bbc23c34d8d67d..02626295a49f314fef840bc044a476f5c9cf332d 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -267,9 +267,9 @@ void runner_dograv_down(struct runner *r, struct cell *c) {
     /* Apply the multipole acceleration to all gparts. */
     for (int k = 0; k < c->gcount; k++) {
       struct gpart *p = &c->gparts[k];
-      p->a[0] += m->a[0];
-      p->a[1] += m->a[1];
-      p->a[2] += m->a[2];
+      p->a_grav[0] += m->a[0];
+      p->a_grav[1] += m->a[1];
+      p->a_grav[2] += m->a[2];
     }
   }
 }
@@ -594,5 +594,4 @@ void runner_dosub_grav(struct runner *r, struct cell *ci, struct cell *cj,
 
   if (gettimer) TIMER_TOC(timer_dosub_grav);
 }
-
 #endif /* SWIFT_RUNNER_DOIACT_GRAV_H */
diff --git a/src/space.c b/src/space.c
index ddb6d67781c2a9cd2bde917ae0dfe350280af649..0d8c171b7ea0fe9226df72153a7d3dcacca7f405 100644
--- a/src/space.c
+++ b/src/space.c
@@ -499,6 +499,28 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
   /* We no longer need the indices as of here. */
   free(gind);
 
+  /* Verify that the links are correct */
+  /* MATTHIEU: To be commented out once we are happy */
+  for (size_t k = 0; k < nr_gparts; ++k) {
+
+    if (s->gparts[k].id > 0) {
+
+      if (s->gparts[k].part->gpart != &s->gparts[k]) error("Linking problem !");
+
+      if (s->gparts[k].x[0] != s->gparts[k].part->x[0] ||
+          s->gparts[k].x[1] != s->gparts[k].part->x[1] ||
+          s->gparts[k].x[2] != s->gparts[k].part->x[2])
+        error("Linked particles are not at the same position !");
+    }
+  }
+  for (size_t k = 0; k < nr_parts; ++k) {
+
+    if (s->parts[k].gpart != NULL) {
+
+      if (s->parts[k].gpart->part != &s->parts[k]) error("Linking problem !");
+    }
+  }
+
   /* Hook the cells up to the parts. */
   // tic = getticks();
   struct part *finger = s->parts;
diff --git a/src/swift.h b/src/swift.h
index 9ab090dccd195ff4927d3e614e446b36d273f824..e568a28c888295affc9ec45b6d059d34f5b4bf04 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -27,7 +27,6 @@
 #include "cell.h"
 #include "clocks.h"
 #include "const.h"
-#include "const.h"
 #include "cycle.h"
 #include "debug.h"
 #include "engine.h"
@@ -38,7 +37,9 @@
 #include "map.h"
 #include "multipole.h"
 #include "parallel_io.h"
+#include "parser.h"
 #include "part.h"
+#include "partition.h"
 #include "queue.h"
 #include "runner.h"
 #include "scheduler.h"
@@ -47,9 +48,8 @@
 #include "space.h"
 #include "task.h"
 #include "timers.h"
-#include "units.h"
 #include "tools.h"
-#include "partition.h"
+#include "units.h"
 #include "version.h"
 
 #endif /* SWIFT_SWIFT_H */
diff --git a/src/tools.c b/src/tools.c
index a384974fdc94452079838ae0eebf30e1815f644b..d25b7401a1e0515c650333b41193d54b5e155d39 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -300,9 +300,9 @@ void pairs_single_grav(double *dim, long long int pid,
       break;
   if (k == N) error("Part not found.");
   pi = parts[k];
-  pi.a[0] = 0.0f;
-  pi.a[1] = 0.0f;
-  pi.a[2] = 0.0f;
+  pi.a_grav[0] = 0.0f;
+  pi.a_grav[1] = 0.0f;
+  pi.a_grav[2] = 0.0f;
 
   /* Loop over all particle pairs. */
   for (k = 0; k < N; k++) {
@@ -320,15 +320,15 @@ void pairs_single_grav(double *dim, long long int pid,
     }
     r2 = fdx[0] * fdx[0] + fdx[1] * fdx[1] + fdx[2] * fdx[2];
     runner_iact_grav(r2, fdx, &pi, &pj);
-    a[0] += pi.a[0];
-    a[1] += pi.a[1];
-    a[2] += pi.a[2];
-    aabs[0] += fabsf(pi.a[0]);
-    aabs[1] += fabsf(pi.a[1]);
-    aabs[2] += fabsf(pi.a[2]);
-    pi.a[0] = 0.0f;
-    pi.a[1] = 0.0f;
-    pi.a[2] = 0.0f;
+    a[0] += pi.a_grav[0];
+    a[1] += pi.a_grav[1];
+    a[2] += pi.a_grav[2];
+    aabs[0] += fabsf(pi.a_grav[0]);
+    aabs[1] += fabsf(pi.a_grav[1]);
+    aabs[2] += fabsf(pi.a_grav[2]);
+    pi.a_grav[0] = 0.0f;
+    pi.a_grav[1] = 0.0f;
+    pi.a_grav[2] = 0.0f;
   }
 
   /* Dump the result. */
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 242341ddf53e8e5b4d0e9cc855229cc4aa469e2b..586a86b984add85f6b3ed1413b85573829a28ecb 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -22,11 +22,11 @@ AM_LDFLAGS = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
 
 # List of programs and scripts to run in the test suite
 TESTS = testGreetings testReading.sh testSingle testPair.sh testPairPerturbed.sh \
-	test27cells.sh test27cellsPerturbed.sh
+	test27cells.sh test27cellsPerturbed.sh testParser.sh
 
 # List of test programs to compile
 check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
-		 testSPHStep testPair test27cells testKernel
+		 testSPHStep testPair test27cells testParser testKernel
 
 # Sources for the individual programs
 testGreetings_SOURCES = testGreetings.c
@@ -43,8 +43,11 @@ testPair_SOURCES = testPair.c
 
 test27cells_SOURCES = test27cells.c
 
+testParser_SOURCES = testParser.c
+
 testKernel_SOURCES = testKernel.c
 
 # Files necessary for distribution
 EXTRA_DIST = testReading.sh makeInput.py testPair.sh testPairPerturbed.sh \
-	     test27cells.sh test27cellsPerturbed.sh tolerance.dat
+	     test27cells.sh test27cellsPerturbed.sh tolerance.dat testParser.sh \ 
+	     testParserInput.yaml
diff --git a/tests/testParser.c b/tests/testParser.c
new file mode 100644
index 0000000000000000000000000000000000000000..a4b8789fca056fef659bca78eae9d0effb2ceb66
--- /dev/null
+++ b/tests/testParser.c
@@ -0,0 +1,67 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2016 James Willis (james.s.willis@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#include "parser.h"
+#include <assert.h>
+#include <string.h>
+#include <math.h>
+
+int main(int argc, char *argv[]) {
+
+  const char *input_file = argv[1];
+
+  /* Create a structure to read file into. */
+  struct swift_params param_file;
+
+  /* Create variables that will be set from the parameter file. */
+  int no_of_threads = 0;
+  int no_of_time_steps = 0;
+  float max_h = 0.0f;
+  double start_time = 0.0;
+  char ic_file[PARSER_MAX_LINE_SIZE];
+
+  /* Read the parameter file. */
+  parser_read_file(input_file, &param_file);
+
+  /* Print the contents of the structure. */
+  parser_print_params(&param_file);
+
+  /* Retrieve parameters and store them in variables defined above.
+   * Have to specify the name of the parameter as it appears in the
+   * input file: testParserInput.yaml.*/
+  parser_get_param_int(&param_file, "no_of_threads", &no_of_threads);
+  parser_get_param_int(&param_file, "no_of_time_steps", &no_of_time_steps);
+  parser_get_param_float(&param_file, "max_h", &max_h);
+  parser_get_param_double(&param_file, "start_time", &start_time);
+  parser_get_param_string(&param_file, "ic_file", ic_file);
+
+  /* Print the variables to check their values are correct. */
+  printf(
+      "no_of_threads: %d, no_of_time_steps: %d, max_h: %f, start_time: %lf, "
+      "ic_file: %s\n",
+      no_of_threads, no_of_time_steps, max_h, start_time, ic_file);
+
+  assert(no_of_threads == 16);
+  assert(no_of_time_steps == 10);
+  assert(fabs(max_h - 1.1255) < 0.00001);
+  assert(fabs(start_time - 1.23456789) < 0.00001);
+  assert(strcmp(ic_file, "ic_file.ini") == 0); /*strcmp returns 0 if correct.*/
+
+  return 0;
+}
diff --git a/tests/testParser.sh b/tests/testParser.sh
new file mode 100755
index 0000000000000000000000000000000000000000..3dad7f386f792ff2beb6e94eb093bad4085023a4
--- /dev/null
+++ b/tests/testParser.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+./testParser testParserInput.yaml
diff --git a/tests/testParserInput.yaml b/tests/testParserInput.yaml
new file mode 100644
index 0000000000000000000000000000000000000000..d695e6a8ddd327e31224f36a6e34767ea8d36408
--- /dev/null
+++ b/tests/testParserInput.yaml
@@ -0,0 +1,9 @@
+---
+no_of_threads:      16 # The number of threads that will be used. 
+no_of_time_steps:   10
+max_h:              1.1255
+start_time:         1.23456789
+#Input file
+ic_file:            ic_file.ini
+
+...