diff --git a/examples/main.c b/examples/main.c
index 3c1921a296c67cf577b3043202bdd23b63e1f8c1..d7a6dd3d47a874d92fa5ed162ee1d67d5a56f2e0 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -599,28 +599,30 @@ int main(int argc, char *argv[]) {
       phys_const_print(&prog_const);
     }
 
-  /* Initialise the cosmology */
-  struct cosmology cosmo;
-  if (with_cosmology) cosmology_init(params, &us, &prog_const, &cosmo);
-  if (with_cosmology) cosmology_print(&cosmo);
-
-  struct engine ee;
-  ee.ti_current = 0;
-  ee.timeBase = (log(cosmo.a_end) - log(cosmo.a_begin)) / max_nr_timesteps;
-  cosmology_update(&cosmo, &ee);
+    /* Initialise the cosmology */
+    struct cosmology cosmo;
+    if (with_cosmology) cosmology_init(params, &us, &prog_const, &cosmo);
+    if (with_cosmology) cosmology_print(&cosmo);
+
+    // MATTHIEU START
+    struct engine ee;
+    ee.ti_current = 0;
+    ee.timeBase = (log(cosmo.a_end) - log(cosmo.a_begin)) / max_nr_timesteps;
+    cosmology_update(&cosmo, &ee);
 
-  for (int i = 0; i <= 16; ++i) {
+    for (int i = 0; i <= 16; ++i) {
 
-    ee.ti_current = (max_nr_timesteps / 16) * i;
+      ee.ti_current = (max_nr_timesteps / 16) * i;
 
-    cosmology_update(&cosmo, &ee);
+      cosmology_update(&cosmo, &ee);
 
-    message("z=%e H(z)=%e w=%f t=%e [yrs] t_l=%e [yrs]", cosmo.z, cosmo.H,
-            cosmo.w, cosmo.time / prog_const.const_year,
-            cosmo.lookback_time / prog_const.const_year);
-  }
+      message("z=%e H(z)=%e w=%f t=%e [yrs] t_l=%e [yrs]", cosmo.z, cosmo.H,
+              cosmo.w, cosmo.time / prog_const.const_year,
+              cosmo.lookback_time / prog_const.const_year);
+    }
 
-  return 0;
+    return 0;
+    // MATTHIEU END
 
     /* Initialise the hydro properties */
     if (with_hydro) hydro_props_init(&hydro_properties, params);
@@ -702,13 +704,25 @@ int main(int argc, char *argv[]) {
           "ICs.",
           N_total[0], N_total[2], N_total[1]);
 
-  /* Also update the total counts (in case of changes due to replication) */
+    /* 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);
+    }
+
+/* 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);
+    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;
@@ -743,10 +757,10 @@ int main(int argc, char *argv[]) {
       message("nr of cells at depth %i is %i.", data[0], data[1]);
     }
 
-  /* Initialise the table of Ewald corrections for the gravity checks */
-#ifdef SWIFT_GRAVITY_FORCE_CHECKS
-  if (periodic) gravity_exact_force_ewald_init(dim[0]);
-#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);
diff --git a/src/cosmology.c b/src/cosmology.c
index 1149737919f8cb00507feb80a31ad04db501ce6c..44c1f493c386c4ef25b1c2690d9d8d77cbdf9524 100644
--- a/src/cosmology.c
+++ b/src/cosmology.c
@@ -82,7 +82,7 @@ static INLINE double cosmology_dark_energy_EoS(double a, double w_0,
 }
 
 /**
- * @brief Computes the integral of the dark-energy equation of state 
+ * @brief Computes the integral of the dark-energy equation of state
  * up to a scale-factor a.
  *
  * We follow the convention of Linder & Jenkins, MNRAS, 346, 573, 2003
@@ -100,7 +100,7 @@ static INLINE double w_tilde(double a, double w0, double wa) {
  * @brief Compute \f$ E(z) \f$.
  */
 static INLINE double E(double Or, double Om, double Ok, double Ol, double w0,
-		       double wa, double a) {
+                       double wa, double a) {
   const double a_inv = 1. / a;
   return sqrt(Or * a_inv * a_inv * a_inv * a_inv + Om * a_inv * a_inv * a_inv +
               Ok * a_inv * a_inv + Ol * exp(3. * w_tilde(a, w0, wa)));