diff --git a/configure.ac b/configure.ac
index fdd263bff26a1a614d7cefe0d5fc4088c3fa8893..a35ea7a1e7fc7fb7df92b2f31b4fe4fbd2012494 100644
--- a/configure.ac
+++ b/configure.ac
@@ -273,6 +273,18 @@ if test "$enable_debugging_checks" = "yes"; then
    AC_DEFINE([SWIFT_DEBUG_CHECKS],1,[Enable expensive debugging])
 fi
 
+# Check if using our custom icbrtf is enalbled.
+AC_ARG_ENABLE([custom-icbrtf],
+   [AS_HELP_STRING([--enable-custom-icbrtf],
+     [Use SWIFT's custom icbrtf function instead of the system cbrtf @<:@yes/no@:>@]
+   )],
+   [enable_custom_icbrtf="$enableval"],
+   [enable_custom_icbrtf="no"]
+)
+if test "$enable_custom_icbrtf" = "yes"; then
+   AC_DEFINE([WITH_ICBRTF],1,[Enable custom icbrtf])
+fi
+
 # Check whether we want to default to naive cell interactions
 AC_ARG_ENABLE([naive-interactions],
    [AS_HELP_STRING([--enable-naive-interactions],
@@ -1354,5 +1366,6 @@ AC_MSG_RESULT([
    Interaction debugging : $enable_debug_interactions
    Naive interactions    : $enable_naive_interactions
    Gravity checks        : $gravity_force_checks
+   Custom icbrtf         : $enable_custom_icbrtf
 
  ------------------------])
diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index bef8ce13fadeb0f61cc08f17eb1330942d98453d..5245c738e32b8f94bf58d7436d76b007b6808c1a 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -1957,7 +1957,7 @@ MACRO_EXPANSION        = YES
 # The default value is: NO.
 # This tag requires that the tag ENABLE_PREPROCESSING is set to YES.
 
-EXPAND_ONLY_PREDEF     = YES
+EXPAND_ONLY_PREDEF     = NO
 
 # If the SEARCH_INCLUDES tag is set to YES, the include files in the
 # INCLUDE_PATH will be searched if a #include is found.
diff --git a/examples/CoolingHalo/cooling_halo.yml b/examples/CoolingHalo/cooling_halo.yml
index e4e8750f46bf7f46c692c4ccd3afe8453e0f606a..68c3478b717261698ac175835fc246e134e3a6a7 100644
--- a/examples/CoolingHalo/cooling_halo.yml
+++ b/examples/CoolingHalo/cooling_halo.yml
@@ -31,15 +31,9 @@ SPH:
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  CoolingHalo.hdf5       # The file to read
-  shift_x:    0.                  # A shift to apply to all particles read from the ICs (in internal units).
-  shift_y:    0.
-  shift_z:    0.
  
 # External potential parameters
 IsothermalPotential:
-  position_x:      0.     # location of centre of isothermal potential in internal units
-  position_y:      0.
-  position_z:      0.	
   vrot:            200.     # rotation speed of isothermal potential in internal units
   timestep_mult:   0.03     # controls time step
   epsilon:         0.1    #softening for the isothermal potential
diff --git a/examples/CoolingHaloWithSpin/cooling_halo.yml b/examples/CoolingHaloWithSpin/cooling_halo.yml
index 1d548b2fb7e531b712548eb5834b2c4de575f941..f6e9fe3b124631fc2d5336db8a7ffb18f7b34a95 100644
--- a/examples/CoolingHaloWithSpin/cooling_halo.yml
+++ b/examples/CoolingHaloWithSpin/cooling_halo.yml
@@ -34,9 +34,6 @@ InitialConditions:
  
 # External potential parameters
 IsothermalPotential:
-  position_x:      0.     # Location of centre of isothermal potential in internal units
-  position_y:      0.
-  position_z:      0.	
   vrot:            200.   # Rotation speed of isothermal potential in internal units
   timestep_mult:   0.03   # Controls time step
   epsilon:         1.0    # Softening for the isothermal potential
diff --git a/examples/ExternalPointMass/externalPointMass.yml b/examples/ExternalPointMass/externalPointMass.yml
index 6e2315cc4dd5c2c354375d39e03f13a1553ed08b..79dcd6458489958622ed7d6167d7f5f214ebd8b9 100644
--- a/examples/ExternalPointMass/externalPointMass.yml
+++ b/examples/ExternalPointMass/externalPointMass.yml
@@ -31,15 +31,11 @@ SPH:
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  PointMass.hdf5        # The file to read
-  shift_x:    50.                   # A shift to apply to all particles read from the ICs (in internal units).
-  shift_y:    50.
-  shift_z:    50.
+  shift:      [50.,50.,50.]         # A shift to apply to all particles read from the ICs (in internal units).
 
 # External potential parameters
 PointMassPotential:
-  position_x:      50.     # location of external point mass in internal units
-  position_y:      50.
-  position_z:      50.	
+  position:        [50.,50.,50.]    # location of external point mass in internal units
   mass:            1e10     # mass of external point mass in internal units
   timestep_mult:   0.03     # controls time step
 
diff --git a/examples/HydrostaticHalo/hydrostatic.yml b/examples/HydrostaticHalo/hydrostatic.yml
index 38a5df2866d39523e6e4b9f6c1c97a7537bff5a4..c37a5314976c3116eb51040b1feeaa9b23ac1326 100644
--- a/examples/HydrostaticHalo/hydrostatic.yml
+++ b/examples/HydrostaticHalo/hydrostatic.yml
@@ -34,9 +34,6 @@ InitialConditions:
  
 # External potential parameters
 IsothermalPotential:
-  position_x:      0.     # location of centre of isothermal potential in internal units
-  position_y:      0.
-  position_z:      0.	
   vrot:            200.     # rotation speed of isothermal potential in internal units
   epsilon:         1.0
   timestep_mult:   0.03     # controls time step
diff --git a/examples/IsothermalPotential/isothermal.yml b/examples/IsothermalPotential/isothermal.yml
index bd59ddc5117f71c478ce5eeda0a08817eefe8624..5dd0b831c839b3307a1c118d9bd64bda29da487c 100644
--- a/examples/IsothermalPotential/isothermal.yml
+++ b/examples/IsothermalPotential/isothermal.yml
@@ -26,15 +26,10 @@ Snapshots:
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  Isothermal.hdf5       # The file to read
-  shift_x:    200.                  # Shift all particles to be in the potential
-  shift_y:    200.
-  shift_z:    200.
+  shift:      [200.,200.,200.]      # Shift all particles to be in the potential
  
 # External potential parameters
 IsothermalPotential:
-  position_x:      0.       # location of centre of isothermal potential in internal units
-  position_y:      0.
-  position_z:      0.
-  vrot:            200.     # rotation speed of isothermal potential in internal units
-  timestep_mult:   0.01     # controls time step
-  epsilon:         0.       # No softening at the centre of the halo
+  vrot:            200.       # rotation speed of isothermal potential in internal units
+  timestep_mult:   0.01       # controls time step
+  epsilon:         0.         # No softening at the centre of the halo
diff --git a/examples/KeplerianRing/keplerian_ring.yml b/examples/KeplerianRing/keplerian_ring.yml
index 421d1e89255195cd05a66975d838c59a4ad77c72..cc5db2a06adbe9678207454c6504a6fa315675cf 100644
--- a/examples/KeplerianRing/keplerian_ring.yml
+++ b/examples/KeplerianRing/keplerian_ring.yml
@@ -32,15 +32,10 @@ SPH:
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  initial_conditions.hdf5        # The file to read
-  shift_x:    0.                   # A shift to apply to all particles read from the ICs (in internal units).
-  shift_y:    0.
-  shift_z:    0.
 
 # External potential parameters
 PointMassPotential:
-  position_x:      5.     # location of external point mass in internal units
-  position_y:      5.     # here just take this to be the centre of the ring
-  position_z:      5.	
-  mass:            1.498351e7     # mass of external point mass in internal units
-  timestep_mult:   0.03     # controls time step
+  position:        [5.,5.,5.]  # location of external point mass in internal units
+  mass:            1.498351e7  # mass of external point mass in internal units
+  timestep_mult:   0.03        # controls time step
   softening:       0.05
diff --git a/examples/KeplerianRing/write_gadget.py b/examples/KeplerianRing/write_gadget.py
index d2519cdaf4d78d9c7327419283072b8001e70fcb..0f745782ef81c7500e8b965938eb5eea33ede288 100644
--- a/examples/KeplerianRing/write_gadget.py
+++ b/examples/KeplerianRing/write_gadget.py
@@ -300,7 +300,7 @@ def write_block(f, part_type, pos, vel, ids, mass, int_energy, smoothing, other=
     else:
         data = default_data
 
-    particles = f.create_group(f"PartType{part_type}")
+    particles = f.create_group("PartType" + str(part_type))
 
     for name, value in data.items():
         particles.create_dataset(name, data=value)
diff --git a/examples/MultiTypes/multiTypes.yml b/examples/MultiTypes/multiTypes.yml
index 5872634cd620e0b7e2b7edb3a38e363c34229b9a..04647f0f00e69f5baf2560aca0feeb14a26cc50a 100644
--- a/examples/MultiTypes/multiTypes.yml
+++ b/examples/MultiTypes/multiTypes.yml
@@ -35,8 +35,6 @@ InitialConditions:
 
 # External potential parameters
 PointMassPotential:
-  position_x:      50.     # location of external point mass in internal units
-  position_y:      50.
-  position_z:      50.	
+  position:        [50.,50.,50.] # location of external point mass in internal units
   mass:            1e10     # mass of external point mass in internal units
   timestep_mult:   1e-2
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index baae1e84e917267c65f84e61038efedce256322e..b6d9228dd42e231da79df47334bfa7b04b4d563b 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -97,9 +97,7 @@ InitialConditions:
   cleanup_velocity_factors:    0    # (Optional) Clean up the scale-factors used in the definition of the velocity variable in the ICs (e.g. in Gadget files).
   cleanup_smoothing_lengths:   0    # (Optional) Clean the values of the smoothing lengths that are read in to remove stupid values. Set to 1 to activate.
   smoothing_length_scaling:    1.   # (Optional) A scaling factor to apply to all smoothing lengths in the ICs.
-  shift_x:    0.                    # (Optional) A shift to apply to all particles read from the ICs (in internal units).
-  shift_y:    0.
-  shift_z:    0.
+  shift:      [0.0,0.0,0.0]         # (Optional) A shift to apply to all particles read from the ICs (in internal units).
   replicate:  2                     # (Optional) Replicate all particles along each axis a given integer number of times. Default 1.
 
 # Parameters controlling restarts
@@ -116,9 +114,7 @@ Restarts:
 DomainDecomposition:
   initial_type:     simple_metis # (Optional) The initial decomposition strategy: "grid",
                                  #            "simple_metis", "weighted_metis", or "vectorized".
-  initial_grid_x:   10      # (Optional) Grid size if the "grid" strategy is chosen.
-  initial_grid_y:   10      # ""
-  initial_grid_z:   10      # ""
+  initial_grid: [10,10,10] # (Optional) Grid sizes if the "grid" strategy is chosen.
 
   repartition_type: costs/costs # (Optional) The re-decomposition strategy, one of:
                             # "none/none", "costs/costs", "counts/none", "none/costs", "counts/costs",
@@ -149,17 +145,14 @@ EoS:
 
 # Point mass external potentials
 PointMassPotential:
-  position_x:      50.      # location of external point mass (internal units)
-  position_y:      50.
-  position_z:      50.
-  mass:            1e10     # mass of external point mass (internal units)
-  timestep_mult:   0.03     # Dimensionless pre-factor for the time-step condition
+  position:        [50.,50.0,50.]      # location of external point mass (internal units)
+  mass:            1e10                # mass of external point mass (internal units)
+  timestep_mult:   0.03                # Dimensionless pre-factor for the time-step condition
+  softening:       0.05                # For point-mass-softened option
 
 # Isothermal potential parameters
 IsothermalPotential:
-  position_x:      100.     # Location of centre of isothermal potential with respect to centre of the box (internal units)
-  position_y:      100.
-  position_z:      100.
+  position:        [100.,100.,100.]    # Location of centre of isothermal potential with respect to centre of the box (internal units)
   vrot:            200.     # Rotation speed of isothermal potential (internal units)
   timestep_mult:   0.03     # Dimensionless pre-factor for the time-step condition
   epsilon:         0.1      # Softening size (internal units)
diff --git a/src/Makefile.am b/src/Makefile.am
index 33bcd78dfe4a1f846a297f47141c2e9914f97f50..9b043c829a20e99306eac7fa7b352859248ab7ae 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -48,7 +48,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.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 cosmology.h restart.h space_getsid.h utilities.h \
-    mesh_gravity.h
+    mesh_gravity.h cbrt.h
 
 # Common source files
 AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
diff --git a/src/adiabatic_index.h b/src/adiabatic_index.h
index cb3cabe357567eb9bc0a6f99dc8a7e1701eda094..b3bb4a4c1c861a57e58814049906bcf1af8366ec 100644
--- a/src/adiabatic_index.h
+++ b/src/adiabatic_index.h
@@ -33,6 +33,7 @@
 #include <math.h>
 
 /* Local headers. */
+#include "cbrt.h"
 #include "error.h"
 #include "inline.h"
 
@@ -112,8 +113,13 @@ __attribute__((always_inline, const)) INLINE static float pow_gamma(float x) {
 
 #if defined(HYDRO_GAMMA_5_3)
 
+#ifdef WITH_ICBRTF
+  const float icbrt = icbrtf(x); /* x^(-1/3) */
+  return icbrt * x * x;          /* x^(5/3) */
+#else
   const float cbrt = cbrtf(x); /* x^(1/3) */
   return cbrt * cbrt * x;      /* x^(5/3) */
+#endif	// WITH_ICBRTF
 
 #elif defined(HYDRO_GAMMA_7_5)
 
@@ -121,7 +127,12 @@ __attribute__((always_inline, const)) INLINE static float pow_gamma(float x) {
 
 #elif defined(HYDRO_GAMMA_4_3)
 
+#ifdef WITH_ICBRTF
+  const float icbrt = icbrtf(x); /* x^(-1/3) */
+  return icbrt * icbrt * x * x;  /* x^(4/3) */
+#else
   return cbrtf(x) * x; /* x^(4/3) */
+#endif	// WITH_ICBRTF
 
 #elif defined(HYDRO_GAMMA_2_1)
 
@@ -146,8 +157,13 @@ __attribute__((always_inline, const)) INLINE static float pow_gamma_minus_one(
 
 #if defined(HYDRO_GAMMA_5_3)
 
+#ifdef WITH_ICBRTF
+  const float icbrt = icbrtf(x); /* x^(-1/3) */
+  return x * icbrt;              /* x^(2/3) */
+#else
   const float cbrt = cbrtf(x); /* x^(1/3) */
   return cbrt * cbrt;          /* x^(2/3) */
+#endif	// WITH_ICBRTF
 
 #elif defined(HYDRO_GAMMA_7_5)
 
@@ -155,7 +171,12 @@ __attribute__((always_inline, const)) INLINE static float pow_gamma_minus_one(
 
 #elif defined(HYDRO_GAMMA_4_3)
 
+#ifdef WITH_ICBRTF
+  const float icbrt = icbrtf(x); /* x^(-1/3) */
+  return x * icbrt * icbrt;      /* x^(1/3) */
+#else
   return cbrtf(x); /* x^(1/3) */
+#endif	// WITH_ICBRTF
 
 #elif defined(HYDRO_GAMMA_2_1)
 
@@ -180,8 +201,13 @@ pow_minus_gamma_minus_one(float x) {
 
 #if defined(HYDRO_GAMMA_5_3)
 
+#ifdef WITH_ICBRTF
+  const float icbrt = icbrtf(x); /* x^(-1/3) */
+  return icbrt * icbrt;          /* x^(-2/3) */
+#else
   const float cbrt_inv = 1.f / cbrtf(x); /* x^(-1/3) */
   return cbrt_inv * cbrt_inv;            /* x^(-2/3) */
+#endif	// WITH_ICBRTF
 
 #elif defined(HYDRO_GAMMA_7_5)
 
@@ -189,7 +215,11 @@ pow_minus_gamma_minus_one(float x) {
 
 #elif defined(HYDRO_GAMMA_4_3)
 
+#ifdef WITH_ICBRTF
+  return icbrtf(x); /* x^(-1/3) */
+#else
   return 1.f / cbrtf(x); /* x^(-1/3) */
+#endif	// WITH_ICBRTF
 
 #elif defined(HYDRO_GAMMA_2_1)
 
@@ -217,9 +247,15 @@ __attribute__((always_inline, const)) INLINE static float pow_minus_gamma(
 
 #if defined(HYDRO_GAMMA_5_3)
 
+#ifdef WITH_ICBRTF
+  const float icbrt = icbrtf(x);      /* x^(-1/3) */
+  const float icbrt2 = icbrt * icbrt; /* x^(-2/3) */
+  return icbrt * icbrt2 * icbrt2;     /* x^(-5/3) */
+#else
   const float cbrt_inv = 1.f / cbrtf(x);       /* x^(-1/3) */
   const float cbrt_inv2 = cbrt_inv * cbrt_inv; /* x^(-2/3) */
   return cbrt_inv * cbrt_inv2 * cbrt_inv2;     /* x^(-5/3) */
+#endif	// WITH_ICBRTF
 
 #elif defined(HYDRO_GAMMA_7_5)
 
@@ -227,7 +263,11 @@ __attribute__((always_inline, const)) INLINE static float pow_minus_gamma(
 
 #elif defined(HYDRO_GAMMA_4_3)
 
+#ifdef WITH_ICBRTF
+  const float cbrt_inv = icbrtf(x);            /* x^(-1/3) */
+#else
   const float cbrt_inv = 1.f / cbrtf(x);       /* x^(-1/3) */
+#endif 	// WITH_ICBRTF
   const float cbrt_inv2 = cbrt_inv * cbrt_inv; /* x^(-2/3) */
   return cbrt_inv2 * cbrt_inv2;                /* x^(-4/3) */
 
diff --git a/src/cbrt.h b/src/cbrt.h
new file mode 100644
index 0000000000000000000000000000000000000000..54560e41d5e6ef143b839f014ce8f0c4a7624ff6
--- /dev/null
+++ b/src/cbrt.h
@@ -0,0 +1,94 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ *                    Matthieu Schaller (matthieu.schaller@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_CBRT_H
+#define SWIFT_CBRT_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "inline.h"
+
+/**
+ * @brief Compute the inverse cube root of a single-precision floating-point
+ * number.
+ *
+ * This function does not care about non-finite inputs.
+ *
+ * @warning This function is faster than both gcc and Intel's `cbrtf()`
+ * functions on x86 systems. However, Other compilers or other architectures
+ * may have faster implementations of the standard function `cbrtf()` that
+ * will potentionally outperform this function.
+ *
+ * @param x_in The input value.
+ *
+ * @return The inverse cubic root of @c x_in (i.e. \f$x_{in}^{-1/3} \f$) .
+ */
+__attribute__((always_inline)) INLINE static float icbrtf(float x_in) {
+
+  union {
+    float as_float;
+    unsigned int as_uint;
+    int as_int;
+  } cast;
+
+  /* Extract the exponent. */
+  cast.as_float = x_in;
+  const int exponent = ((cast.as_int & 0x7f800000) >> 23) - 127;
+
+  /* Clear the exponent and sign to get the mantissa. */
+  cast.as_uint = (cast.as_uint & ~0xff800000) | 0x3f800000;
+  const float x_norm = cast.as_float;
+
+  /* Multiply by sqrt(1/2) and subtract one, should then be in the
+     range [sqrt(1/2) - 1, sqrt(2) - 1). */
+  const float x = x_norm * (float)M_SQRT1_2 - 1.0f;
+
+  /* Compute the polynomial interpolant. */
+  float res =
+      9.99976591940035e-01f +
+      x * (-3.32901212909283e-01f +
+           x * (2.24361110929912e-01f +
+                x * (-1.88913279594895e-01f + x * 1.28384036492344e-01f)));
+
+  /* Compute the new exponent and the correction factor. */
+  int exponent_new = exponent;
+  if (exponent_new < 0) exponent_new -= 2;
+  exponent_new = -exponent_new / 3;
+  const int exponent_rem = exponent + 3 * exponent_new;
+  cast.as_uint = (exponent_new + 127) << 23;
+  static const float scale[3] = {8.90898718140339e-01f, 7.07106781186548e-01f,
+                                 5.61231024154687e-01f};
+  const float exponent_scale = cast.as_float * scale[exponent_rem];
+
+  /* Scale the result and set the correct sign. */
+  res = copysignf(res * exponent_scale, x_in);
+
+  /* One step of Newton iteration to refine the result. */
+  res *= (1.0f / 3.0f) * (4.0f - x_in * res * res * res);
+
+  /* We're done. */
+  return res;
+}
+
+#endif /* SWIFT_CBRT_H */
diff --git a/src/common_io.c b/src/common_io.c
index 494a702125cf873946d06855b5683216cb2aceaf..68311107575a89ce8a2990a8e0f7a8eeb5d2d644 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -250,13 +250,18 @@ void io_write_attribute_s(hid_t grp, const char* name, const char* str) {
 
 /**
  * @brief Reads the Unit System from an IC file.
+ *
+ * If the 'Units' group does not exist in the ICs, we will use the internal
+ * system of units.
+ *
  * @param h_file The (opened) HDF5 file from which to read.
- * @param us The unit_system to fill.
+ * @param ic_units The unit_system to fill.
+ * @param internal_units The internal system of units to copy if needed.
  * @param mpi_rank The MPI rank we are on.
- *
- * If the 'Units' group does not exist in the ICs, cgs units will be assumed
  */
-void io_read_unit_system(hid_t h_file, struct unit_system* us, int mpi_rank) {
+void io_read_unit_system(hid_t h_file, struct unit_system* ic_units,
+                         const struct unit_system* internal_units,
+                         int mpi_rank) {
 
   /* First check if it exists as this is *not* required. */
   const htri_t exists = H5Lexists(h_file, "/Units", H5P_DEFAULT);
@@ -264,16 +269,12 @@ void io_read_unit_system(hid_t h_file, struct unit_system* us, int mpi_rank) {
   if (exists == 0) {
 
     if (mpi_rank == 0)
-      message("'Units' group not found in ICs. Assuming CGS unit system.");
+      message("'Units' group not found in ICs. Assuming internal unit system.");
 
-    /* Default to CGS */
-    us->UnitMass_in_cgs = 1.;
-    us->UnitLength_in_cgs = 1.;
-    us->UnitTime_in_cgs = 1.;
-    us->UnitCurrent_in_cgs = 1.;
-    us->UnitTemperature_in_cgs = 1.;
+    units_copy(ic_units, internal_units);
 
     return;
+
   } else if (exists < 0) {
     error("Serious problem with 'Units' group in ICs. H5Lexists gives %d",
           exists);
@@ -284,15 +285,15 @@ void io_read_unit_system(hid_t h_file, struct unit_system* us, int mpi_rank) {
 
   /* Ok, Read the damn thing */
   io_read_attribute(h_grp, "Unit length in cgs (U_L)", DOUBLE,
-                    &us->UnitLength_in_cgs);
+                    &ic_units->UnitLength_in_cgs);
   io_read_attribute(h_grp, "Unit mass in cgs (U_M)", DOUBLE,
-                    &us->UnitMass_in_cgs);
+                    &ic_units->UnitMass_in_cgs);
   io_read_attribute(h_grp, "Unit time in cgs (U_t)", DOUBLE,
-                    &us->UnitTime_in_cgs);
+                    &ic_units->UnitTime_in_cgs);
   io_read_attribute(h_grp, "Unit current in cgs (U_I)", DOUBLE,
-                    &us->UnitCurrent_in_cgs);
+                    &ic_units->UnitCurrent_in_cgs);
   io_read_attribute(h_grp, "Unit temperature in cgs (U_T)", DOUBLE,
-                    &us->UnitTemperature_in_cgs);
+                    &ic_units->UnitTemperature_in_cgs);
 
   /* Clean up */
   H5Gclose(h_grp);
diff --git a/src/common_io.h b/src/common_io.h
index f26a635a66f40424984238e586fcdf5bc752fc99..61feffe2ba613c339df02e8c71cbbd1d4aec7a87 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -75,7 +75,9 @@ void io_write_attribute_s(hid_t grp, const char* name, const char* str);
 void io_write_code_description(hid_t h_file);
 void io_write_engine_policy(hid_t h_file, const struct engine* e);
 
-void io_read_unit_system(hid_t h_file, struct unit_system* us, int mpi_rank);
+void io_read_unit_system(hid_t h_file, struct unit_system* ic_units,
+                         const struct unit_system* internal_units,
+                         int mpi_rank);
 void io_write_unit_system(hid_t h_grp, const struct unit_system* us,
                           const char* groupName);
 
diff --git a/src/equation_of_state/planetary/equation_of_state.h b/src/equation_of_state/planetary/equation_of_state.h
index 2a0fe66b49c4ca3aa2135dc773f4c9aae58d3b3f..61e23dc0b4eb82e9ae5c0869f7a10dfff97fc45e 100644
--- a/src/equation_of_state/planetary/equation_of_state.h
+++ b/src/equation_of_state/planetary/equation_of_state.h
@@ -54,7 +54,7 @@ enum eos_planetary_type_id {
   eos_planetary_type_HM80 = 2,
   eos_planetary_type_ANEOS = 3,
   eos_planetary_type_SESAME = 4,
-} __attribute__((packed));
+};
 
 /**
  * @brief Minor type for the planetary equation of state.
@@ -104,7 +104,7 @@ enum eos_planetary_material_id {
   /*! SESAME iron */
   eos_planetary_id_SESAME_iron =
       eos_planetary_type_SESAME * eos_planetary_type_factor,
-} __attribute__((packed));
+};
 
 /* Individual EOS function headers. */
 #include "aneos.h"
@@ -132,7 +132,8 @@ __attribute__((always_inline)) INLINE static float
 gas_internal_energy_from_entropy(float density, float entropy,
                                  enum eos_planetary_material_id mat_id) {
 
-  const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+  const enum eos_planetary_type_id type =
+      (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
 
   /* Select the material base type */
   switch (type) {
@@ -241,7 +242,8 @@ gas_internal_energy_from_entropy(float density, float entropy,
 __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
     float density, float entropy, enum eos_planetary_material_id mat_id) {
 
-  const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+  const enum eos_planetary_type_id type =
+      (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
 
   /* Select the material base type */
   switch (type) {
@@ -344,7 +346,8 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
 __attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
     float density, float P, enum eos_planetary_material_id mat_id) {
 
-  const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+  const enum eos_planetary_type_id type =
+      (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
 
   /* Select the material base type */
   switch (type) {
@@ -445,7 +448,8 @@ __attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
 __attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
     float density, float entropy, enum eos_planetary_material_id mat_id) {
 
-  const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+  const enum eos_planetary_type_id type =
+      (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
 
   /* Select the material base type */
   switch (type) {
@@ -549,8 +553,8 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
 __attribute__((always_inline)) INLINE static float
 gas_entropy_from_internal_energy(float density, float u,
                                  enum eos_planetary_material_id mat_id) {
-
-  const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+  const enum eos_planetary_type_id type =
+      (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
 
   /* Select the material base type */
   switch (type) {
@@ -654,7 +658,8 @@ __attribute__((always_inline)) INLINE static float
 gas_pressure_from_internal_energy(float density, float u,
                                   enum eos_planetary_material_id mat_id) {
 
-  const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+  const enum eos_planetary_type_id type =
+      (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
 
   /* Select the material base type */
   switch (type) {
@@ -762,7 +767,8 @@ __attribute__((always_inline)) INLINE static float
 gas_internal_energy_from_pressure(float density, float P,
                                   enum eos_planetary_material_id mat_id) {
 
-  const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+  const enum eos_planetary_type_id type =
+      (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
 
   /* Select the material base type */
   switch (type) {
@@ -867,7 +873,8 @@ __attribute__((always_inline)) INLINE static float
 gas_soundspeed_from_internal_energy(float density, float u,
                                     enum eos_planetary_material_id mat_id) {
 
-  const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+  const enum eos_planetary_type_id type =
+      (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
 
   /* Select the material base type */
   switch (type) {
@@ -975,7 +982,8 @@ gas_soundspeed_from_internal_energy(float density, float u,
 __attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
     float density, float P, enum eos_planetary_material_id mat_id) {
 
-  const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
+  const enum eos_planetary_type_id type =
+      (enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
 
   /* Select the material base type */
   switch (type) {
diff --git a/src/hydro/MinimalMultiMat/hydro_io.h b/src/hydro/MinimalMultiMat/hydro_io.h
index 61301157e12b625f3ca1fb6cf5d5238d82403c9f..7f41f5e227b6c8a8904b5546a2568b4700109abd 100644
--- a/src/hydro/MinimalMultiMat/hydro_io.h
+++ b/src/hydro/MinimalMultiMat/hydro_io.h
@@ -69,8 +69,8 @@ INLINE static void hydro_read_particles(struct part* parts,
                                 UNIT_CONV_ACCELERATION, parts, a_hydro);
   list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
                                 UNIT_CONV_DENSITY, parts, rho);
-  list[8] =
-      io_make_input_field("MaterialID", INT, 1, OPTIONAL, 1, parts, mat_id);
+  list[8] = io_make_input_field("MaterialID", INT, 1, COMPULSORY,
+                                UNIT_CONV_NO_UNITS, parts, mat_id);
 }
 
 INLINE static void convert_S(const struct engine* e, const struct part* p,
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 749d9020504be97f27f4f7da064148ce7d0ca48f..7e892154231cb09ffca1f5f93f56240b043995ba 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -730,7 +730,7 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
   struct unit_system* ic_units =
       (struct unit_system*)malloc(sizeof(struct unit_system));
   if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
-  io_read_unit_system(h_file, ic_units, mpi_rank);
+  io_read_unit_system(h_file, ic_units, internal_units, mpi_rank);
 
   /* Tell the user if a conversion will be needed */
   if (mpi_rank == 0) {
diff --git a/src/parser.c b/src/parser.c
index 78d8aef2c3194acd0a9128867e6e5867a0cbc7b0..04015fe693fbf616e531681a3b46902812bc09e0 100644
--- a/src/parser.c
+++ b/src/parser.c
@@ -1,7 +1,7 @@
 /*******************************************************************************
  * This file is part of SWIFT.
  * Copyright (c) 2016 James Willis (james.s.willis@durham.ac.uk)
- *               2017 Peter W. Draper (p.w.draper@durham.ac.uk)
+ *               2017-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
@@ -34,6 +34,7 @@
 #include "common_io.h"
 #include "error.h"
 #include "restart.h"
+#include "tools.h"
 
 #define PARSER_COMMENT_STRING "#"
 #define PARSER_COMMENT_CHAR '#'
@@ -42,8 +43,9 @@
 #define PARSER_START_OF_FILE "---"
 #define PARSER_END_OF_FILE "..."
 
+#define CHUNK 10
+
 /* Private functions. */
-static int count_char(const char *str, char val);
 static int is_empty(const char *str);
 static int count_indentation(const char *str);
 static void parse_line(char *line, struct swift_params *params);
@@ -56,6 +58,86 @@ static void find_duplicate_section(const struct swift_params *params,
                                    const char *section_name);
 static int lineNumber = 0;
 
+/**
+ * @brief parse a YAML list of strings returning a set of pointers to
+ *        the strings.
+ *
+ * It is assumed that the [] have been removed (also no lists in lists)
+ * words are separated by commas and the strings may or may not be quoted.
+ * So lines like:
+ *
+ *    'xyz', 'ABC', "ab'c", "de:f", "g,hi", "zzz", Hello World, again
+ *
+ * Are supported as expected.
+ *
+ * @param line the line to parse.
+ * @param result array of pointers to the strings.
+ * @return the number of strings
+ */
+static int parse_quoted_strings(const char *line, char ***result) {
+
+  char word[PARSER_MAX_LINE_SIZE];
+  int nchar = 0;
+  int nwords = 0;
+  char quote = '\0';
+
+  /* Preallocate a number of pointers. */
+  char **strings;
+  int count = CHUNK;
+  strings = (char **)malloc(count * sizeof(char *));
+
+  word[0] = '\0';
+  for (unsigned int i = 0; i < strlen(line); i++) {
+    char c = line[i];
+    if (c == '"' || c == '\'') {
+      if (c == quote) {
+        quote = '\0';
+      } else if (!quote) {
+        quote = c;
+      } else {
+        word[nchar++] = c;
+      }
+    } else if (c == ',') {
+      if (!quote) {
+
+        /* Save word. */
+        word[nchar++] = '\0';
+        if (count <= nwords) {
+          count += CHUNK;
+          strings = (char **)realloc(strings, count * sizeof(char *));
+        }
+        strings[nwords] = (char *)malloc((strlen(word) + 1) * sizeof(char));
+        strcpy(strings[nwords], trim_both(word));
+        nwords++;
+
+        /* Ready for next. */
+        nchar = 0;
+        word[0] = '\0';
+
+      } else {
+        word[nchar++] = c;
+      }
+    } else {
+      word[nchar++] = c;
+    }
+  }
+
+  /* Keep unfinished words. */
+  if (nchar > 0) {
+    word[nchar] = '\0';
+    if (count <= nwords) {
+      count += 1;
+      strings = (char **)realloc(strings, count * sizeof(char));
+    }
+    strings[nwords] = (char *)malloc((strlen(word) + 1) * sizeof(char));
+    strcpy(strings[nwords], trim_both(word));
+    nwords++;
+  }
+
+  *result = strings;
+  return nwords;
+}
+
 /**
  * @brief Initialize the parser structure.
  *
@@ -115,12 +197,19 @@ void parser_set_param(struct swift_params *params, const char *namevalue) {
   /* Get the various parts. */
   char name[PARSER_MAX_LINE_SIZE];
   char value[PARSER_MAX_LINE_SIZE];
+  char section[PARSER_MAX_LINE_SIZE];
   name[0] = '\0';
   value[0] = '\0';
 
   /* Name is part until second colon. */
   const char *p1 = strchr(namevalue, ':');
   if (p1 != NULL) {
+
+    /* Section is first part until a colon. */
+    memcpy(section, namevalue, p1 - namevalue);
+    section[p1 - namevalue] = ':';
+    section[p1 - namevalue + 1] = '\0';
+
     const char *p2 = strchr(p1 + 1, ':');
     if (p2 != NULL) {
       memcpy(name, namevalue, p2 - namevalue);
@@ -145,11 +234,27 @@ void parser_set_param(struct swift_params *params, const char *namevalue) {
     if (strcmp(name, params->data[i].name) == 0) {
       message("Value of '%s' changed from '%s' to '%s'", params->data[i].name,
               params->data[i].value, value);
-      strcpy(params->data[i].value, value);
+      strcpy(params->data[i].value, trim_both(value));
       updated = 1;
     }
   }
   if (!updated) {
+    /* Is this a new section? */
+    int newsection = 1;
+    for (int i = 0; i < params->sectionCount; i++) {
+      if (strcmp(section, params->section[i].name) == 0) {
+        newsection = 0;
+        break;
+      }
+    }
+    if (newsection) {
+      strcpy(params->section[params->sectionCount].name, section);
+      params->sectionCount++;
+      if (params->sectionCount == PARSER_MAX_NO_OF_SECTIONS)
+        error("Too many sections, current maximum is %d.",
+              params->sectionCount);
+    }
+
     strcpy(params->data[params->paramCount].name, name);
     strcpy(params->data[params->paramCount].value, value);
     params->data[params->paramCount].used = 0;
@@ -159,26 +264,6 @@ void parser_set_param(struct swift_params *params, const char *namevalue) {
   }
 }
 
-/**
- * @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
- *
- * @return Number of occurrences of val inside str
- */
-
-static int count_char(const char *str, char val) {
-  int count = 0;
-
-  /* Check if the line contains the character */
-  while (*str) {
-    if (*str++ == val) ++count;
-  }
-
-  return count;
-}
-
 /**
  * @brief Counts the number of white spaces that prefix a string.
  *
@@ -260,7 +345,7 @@ static void find_duplicate_section(const struct swift_params *params,
 
 static void parse_line(char *line, struct swift_params *params) {
   /* Parse line if it doesn't begin with a comment. */
-  if (*line != PARSER_COMMENT_CHAR) {
+  if (line[0] != PARSER_COMMENT_CHAR) {
     char trim_line[PARSER_MAX_LINE_SIZE];
     char tmp_str[PARSER_MAX_LINE_SIZE];
     char *token;
@@ -306,14 +391,8 @@ static void parse_value(char *line, struct swift_params *params) {
 
   char *token;
 
-  /* Check for more than one value on the same line. */
-  if (count_char(line, PARSER_VALUE_CHAR) > 1) {
-    error("Invalid line:%d '%s', only one value allowed per line.", lineNumber,
-          line);
-  }
-
   /* Check that standalone parameters have correct indentation. */
-  if (!inSection && *line == ' ') {
+  if (!inSection && line[0] == ' ') {
     error(
         "Invalid line:%d '%s', standalone parameter defined with incorrect "
         "indentation.",
@@ -321,18 +400,20 @@ static void parse_value(char *line, struct swift_params *params) {
   }
 
   /* Check that it is a parameter inside a section.*/
-  if (*line == ' ' || *line == '\t') {
+  if (line[0] == ' ' || line[0] == '\t') {
     parse_section_param(line, &isFirstParam, section, params);
-  } else { /*Else it is the start of a new section or standalone parameter. */
-    /* Take first token as the parameter name. */
-    token = strtok(line, " :\t");
-    strcpy(tmpStr, token);
+  } else {
+    /* It is the start of a new section or standalone parameter.
+     * Take first token as the parameter name. */
+    token = strtok(line, ":\t");
+    strcpy(tmpStr, trim_trailing(token));
 
     /* Take second token as the parameter value. */
-    token = strtok(NULL, " #\n");
+    token = trim_both(strtok(NULL, "#\n"));
 
-    /* If second token is NULL then the line must be a section heading. */
-    if (token == NULL) {
+    /* If second token is NULL or empty then the line must be a section
+     * heading. */
+    if (token == NULL || strlen(token) == 0) {
       strcpy(tmpSectionName, tmpStr);
       strcat(tmpSectionName, PARSER_VALUE_STRING);
 
@@ -413,11 +494,11 @@ static void parse_section_param(char *line, int *isFirstParam,
   }
 
   /* Take first token as the parameter name and trim leading white space. */
-  token = strtok(line, " :\t");
+  token = trim_both(strtok(line, ":\t"));
   strcpy(tmpStr, token);
 
   /* Take second token as the parameter value. */
-  token = strtok(NULL, " #\n");
+  token = trim_both(strtok(NULL, "#\n"));
 
   /* Prefix the parameter name with its section name and
    * copy it into the parameter structure. */
@@ -437,6 +518,54 @@ static void parse_section_param(char *line, int *isFirstParam,
   }
 }
 
+// Retrieve parameter value from structure. TYPE is the data type, float, int
+// etc. FMT the format required for that data type, i.e. %f, %d etc. and DESC
+// a one word description of the type, "float", "int" etc.
+#define PARSER_GET_VALUE(TYPE, FMT, DESC)                                     \
+  static int get_param_##TYPE(struct swift_params *params, const char *name,  \
+                              int required, TYPE *result) {                   \
+    char str[PARSER_MAX_LINE_SIZE];                                           \
+    for (int i = 0; i < params->paramCount; i++) {                            \
+      if (strcmp(name, params->data[i].name) == 0) {                          \
+        /* Check that exactly one number is parsed, capture junk. */          \
+        if (sscanf(params->data[i].value, " " FMT "%s ", result, str) != 1) { \
+          error("Tried parsing " DESC                                         \
+                " '%s' but found '%s' with "                                  \
+                "illegal trailing characters '%s'.",                          \
+                params->data[i].name, params->data[i].value, str);            \
+        }                                                                     \
+        /* This parameter has been used */                                    \
+        params->data[i].used = 1;                                             \
+        return 1;                                                             \
+      }                                                                       \
+    }                                                                         \
+    if (required)                                                             \
+      error("Cannot find '%s' in the structure, in file '%s'.", name,         \
+            params->fileName);                                                \
+    return 0;                                                                 \
+  }
+
+// Set a parameter to a value and save for dumping.
+#define PARSER_SAVE_VALUE(PREFIX, TYPE, FMT)                      \
+  static void save_param_##PREFIX(struct swift_params *params,    \
+                                  const char *name, TYPE value) { \
+    char str[PARSER_MAX_LINE_SIZE];                               \
+    sprintf(str, "%s: " FMT, name, value);                        \
+    parser_set_param(params, str);                                \
+    params->data[params->paramCount - 1].used = 1;                \
+  }
+
+/* Instantiations. */
+PARSER_GET_VALUE(char, "%c", "char");
+PARSER_GET_VALUE(int, "%d", "int");
+PARSER_GET_VALUE(float, "%f", "float");
+PARSER_GET_VALUE(double, "%lf", "double");
+PARSER_SAVE_VALUE(char, char, "%c");
+PARSER_SAVE_VALUE(int, int, "%d");
+PARSER_SAVE_VALUE(float, float, "%g");
+PARSER_SAVE_VALUE(double, double, "%g");
+PARSER_SAVE_VALUE(string, const char *, "%s");
+
 /**
  * @brief Retrieve integer parameter from structure.
  *
@@ -445,30 +574,9 @@ static void parse_section_param(char *line, int *isFirstParam,
  * @return Value of the parameter found
  */
 int parser_get_param_int(struct swift_params *params, const char *name) {
-
-  char str[PARSER_MAX_LINE_SIZE];
-  int retParam = 0;
-
-  for (int i = 0; i < params->paramCount; i++) {
-    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);
-      }
-
-      /* this parameter has been used */
-      params->data[i].used = 1;
-
-      return retParam;
-    }
-  }
-
-  error("Cannot find '%s' in the structure, in file '%s'.", name,
-        params->fileName);
-  return 0;
+  int result = 0;
+  get_param_int(params, name, 1, &result);
+  return result;
 }
 
 /**
@@ -479,30 +587,9 @@ int parser_get_param_int(struct swift_params *params, const char *name) {
  * @return Value of the parameter found
  */
 char parser_get_param_char(struct swift_params *params, const char *name) {
-
-  char str[PARSER_MAX_LINE_SIZE];
-  char retParam = 0;
-
-  for (int i = 0; i < params->paramCount; i++) {
-    if (!strcmp(name, params->data[i].name)) {
-      /* Check that exactly one number is parsed. */
-      if (sscanf(params->data[i].value, "%c%s", &retParam, str) != 1) {
-        error(
-            "Tried parsing char '%s' but found '%s' with illegal char "
-            "characters '%s'.",
-            params->data[i].name, params->data[i].value, str);
-      }
-
-      /* this parameter has been used */
-      params->data[i].used = 1;
-
-      return retParam;
-    }
-  }
-
-  error("Cannot find '%s' in the structure, in file '%s'.", name,
-        params->fileName);
-  return 0;
+  char result = 0;
+  get_param_char(params, name, 1, &result);
+  return result;
 }
 
 /**
@@ -513,30 +600,9 @@ char parser_get_param_char(struct swift_params *params, const char *name) {
  * @return Value of the parameter found
  */
 float parser_get_param_float(struct swift_params *params, const char *name) {
-
-  char str[PARSER_MAX_LINE_SIZE];
-  float retParam = 0.f;
-
-  for (int i = 0; i < params->paramCount; i++) {
-    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);
-      }
-
-      /* this parameter has been used */
-      params->data[i].used = 1;
-
-      return retParam;
-    }
-  }
-
-  error("Cannot find '%s' in the structure, in file '%s'.", name,
-        params->fileName);
-  return 0.f;
+  float result = 0;
+  get_param_float(params, name, 1, &result);
+  return result;
 }
 
 /**
@@ -547,30 +613,9 @@ float parser_get_param_float(struct swift_params *params, const char *name) {
  * @return Value of the parameter found
  */
 double parser_get_param_double(struct swift_params *params, const char *name) {
-
-  char str[PARSER_MAX_LINE_SIZE];
-  double retParam = 0.;
-
-  for (int i = 0; i < params->paramCount; i++) {
-    if (!strcmp(name, params->data[i].name)) {
-      /* Check that exactly one number is parsed. */
-      if (sscanf(params->data[i].value, "%lf%s", &retParam, str) != 1) {
-        error(
-            "Tried parsing double '%s' but found '%s' with illegal double "
-            "characters '%s'.",
-            params->data[i].name, params->data[i].value, str);
-      }
-
-      /* this parameter has been used */
-      params->data[i].used = 1;
-
-      return retParam;
-    }
-  }
-
-  error("Cannot find '%s' in the structure, in file '%s'.", name,
-        params->fileName);
-  return 0.;
+  double result = 0;
+  get_param_double(params, name, 1, &result);
+  return result;
 }
 
 /**
@@ -605,36 +650,9 @@ void parser_get_param_string(struct swift_params *params, const char *name,
  */
 int parser_get_opt_param_int(struct swift_params *params, const char *name,
                              int def) {
-
-  char str[PARSER_MAX_LINE_SIZE];
-  int retParam = 0;
-
-  for (int i = 0; i < params->paramCount; i++) {
-    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);
-      }
-
-      /* this parameter has been used */
-      params->data[i].used = 1;
-
-      return retParam;
-    }
-  }
-
-  /* Generate string for new parameter */
-  sprintf(str, "%s: %i", name, def);
-
-  /* Add it to params */
-  parser_set_param(params, str);
-
-  /* Set parameter as used */
-  params->data[params->paramCount - 1].used = 1;
-
+  int result = 0;
+  if (get_param_int(params, name, 0, &result)) return result;
+  save_param_int(params, name, def);
   return def;
 }
 
@@ -648,36 +666,9 @@ int parser_get_opt_param_int(struct swift_params *params, const char *name,
  */
 char parser_get_opt_param_char(struct swift_params *params, const char *name,
                                char def) {
-
-  char str[PARSER_MAX_LINE_SIZE];
-  char retParam = 0;
-
-  for (int i = 0; i < params->paramCount; i++) {
-    if (!strcmp(name, params->data[i].name)) {
-      /* Check that exactly one number is parsed. */
-      if (sscanf(params->data[i].value, "%c%s", &retParam, str) != 1) {
-        error(
-            "Tried parsing char '%s' but found '%s' with illegal char "
-            "characters '%s'.",
-            params->data[i].name, params->data[i].value, str);
-      }
-
-      /* this parameter has been used */
-      params->data[i].used = 1;
-
-      return retParam;
-    }
-  }
-
-  /* Generate string for new parameter */
-  sprintf(str, "%s: %c", name, def);
-
-  /* Add it to params */
-  parser_set_param(params, str);
-
-  /* Set parameter as used */
-  params->data[params->paramCount - 1].used = 1;
-
+  char result = 0;
+  if (get_param_char(params, name, 0, &result)) return result;
+  save_param_char(params, name, def);
   return def;
 }
 
@@ -691,115 +682,392 @@ char parser_get_opt_param_char(struct swift_params *params, const char *name,
  */
 float parser_get_opt_param_float(struct swift_params *params, const char *name,
                                  float def) {
+  float result = 0;
+  if (get_param_float(params, name, 0, &result)) return result;
+  save_param_float(params, name, def);
+  return def;
+}
 
-  char str[PARSER_MAX_LINE_SIZE];
-  float retParam = 0.f;
+/**
+ * @brief Retrieve optional double parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param def Default value of the parameter of not found.
+ * @return Value of the parameter found
+ */
+double parser_get_opt_param_double(struct swift_params *params,
+                                   const char *name, double def) {
+  double result = 0;
+  if (get_param_double(params, name, 0, &result)) return result;
+  save_param_double(params, name, def);
+  return def;
+}
+
+/**
+ * @brief Retrieve string parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param def Default value of the parameter of not found.
+ * @param retParam (return) Value of the parameter found
+ */
+void parser_get_opt_param_string(struct swift_params *params, const char *name,
+                                 char *retParam, const char *def) {
 
   for (int i = 0; i < params->paramCount; i++) {
     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);
-      }
+      strcpy(retParam, params->data[i].value);
 
       /* this parameter has been used */
       params->data[i].used = 1;
-
-      return retParam;
+      return;
     }
   }
+  save_param_string(params, name, def);
+  strcpy(retParam, def);
+}
 
-  /* Generate string for new parameter */
-  sprintf(str, "%s: %f", name, def);
+/* Macro defining functions that get primitive types as simple one-line YAML
+ * arrays, that is SEC: [v1,v2,v3...] format, with the extension that the []
+ * are optional. TYPE is the data type, float etc. FMT a format to parse a
+ * single value, so "%f" for a float and DESC the type description
+ * i.e. "float".
+ */
+#define PARSER_GET_ARRAY(TYPE, FMT, DESC)                             \
+  static int get_param_##TYPE##_array(struct swift_params *params,    \
+                                      const char *name, int required, \
+                                      int nval, TYPE *values) {       \
+    char str[PARSER_MAX_LINE_SIZE];                                   \
+    char cpy[PARSER_MAX_LINE_SIZE];                                   \
+                                                                      \
+    for (int i = 0; i < params->paramCount; i++) {                    \
+      if (!strcmp(name, params->data[i].name)) {                      \
+        char *cp = cpy;                                               \
+        strcpy(cp, params->data[i].value);                            \
+        cp = trim_both(cp);                                           \
+                                                                      \
+        /* Strip off [], if present. */                               \
+        if (cp[0] == '[') cp++;                                       \
+        int l = strlen(cp);                                           \
+        if (cp[l - 1] == ']') cp[l - 1] = '\0';                       \
+        cp = trim_both(cp);                                           \
+                                                                      \
+        /* Format that captures spaces and trailing junk. */          \
+        char fmt[20];                                                 \
+        sprintf(fmt, " %s%%s ", FMT);                                 \
+                                                                      \
+        /* Parse out values which should now be "v, v, v" with        \
+         * internal     whitespace variations. */                     \
+        char *p = strtok(cp, ",");                                    \
+        for (int k = 0; k < nval; k++) {                              \
+          if (p != NULL) {                                            \
+            if (sscanf(p, fmt, &values[k], str) != 1) {               \
+              error("Tried parsing " DESC                             \
+                    " '%s' but found '%s' with "                      \
+                    "illegal " DESC " characters '%s'.",              \
+                    name, p, str);                                    \
+            }                                                         \
+          } else {                                                    \
+            error(                                                    \
+                "Array '%s' with value '%s' has too few values, "     \
+                "expected %d",                                        \
+                name, params->data[i].value, nval);                   \
+          }                                                           \
+          if (k < nval - 1) p = strtok(NULL, ",");                    \
+        }                                                             \
+        params->data[i].used = 1;                                     \
+        return 1;                                                     \
+      }                                                               \
+    }                                                                 \
+    if (required)                                                     \
+      error("Cannot find '%s' in the structure, in file '%s'.", name, \
+            params->fileName);                                        \
+    return 0;                                                         \
+  }
 
-  /* Add it to params */
-  parser_set_param(params, str);
+// Set values of a default parameter so they will be saved correctly.
+#define PARSER_SAVE_ARRAY(TYPE, FMT)                                           \
+  static int save_param_##TYPE##_array(                                        \
+      struct swift_params *params, const char *name, int nval, TYPE *values) { \
+    /* Save values against the parameter. */                                   \
+    char str[PARSER_MAX_LINE_SIZE];                                            \
+    int k = sprintf(str, "%s: [", name);                                       \
+    for (int i = 0; i < nval - 1; i++)                                         \
+      k += sprintf(&str[k], FMT ", ", values[i]);                              \
+    sprintf(&str[k], FMT "]", values[nval - 1]);                               \
+    parser_set_param(params, str);                                             \
+    params->data[params->paramCount - 1].used = 1;                             \
+    return 0;                                                                  \
+  }
 
-  /* Set parameter as used */
-  params->data[params->paramCount - 1].used = 1;
+/* Instantiations. */
+PARSER_GET_ARRAY(char, "%c", "char");
+PARSER_GET_ARRAY(int, "%d", "int");
+PARSER_GET_ARRAY(float, "%f", "float");
+PARSER_GET_ARRAY(double, "%lf", "double");
+PARSER_SAVE_ARRAY(char, "%c");
+PARSER_SAVE_ARRAY(int, "%d");
+PARSER_SAVE_ARRAY(float, "%g");
+PARSER_SAVE_ARRAY(double, "%g");
 
-  return def;
+/**
+ * @brief Retrieve char array parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param nval number of values expected.
+ * @param values Values of the parameter found, of size at least nvals.
+ */
+void parser_get_param_char_array(struct swift_params *params, const char *name,
+                                 int nval, char *values) {
+  get_param_char_array(params, name, 1, nval, values);
 }
 
 /**
- * @brief Retrieve optional double parameter from structure.
+ * @brief Retrieve optional char array parameter from structure.
  *
  * @param params Structure that holds the parameters
  * @param name Name of the parameter to be found
- * @param def Default value of the parameter of not found.
- * @return Value of the parameter found
+ * @param nval number of values expected.
+ * @param values Values of the parameter found, of size at least nvals. If the
+ *               parameter is not found these values will be returned
+ *               unmodified, so should be set to the default values.
+ * @return whether the parameter has been found.
  */
-double parser_get_opt_param_double(struct swift_params *params,
-                                   const char *name, double def) {
-
-  char str[PARSER_MAX_LINE_SIZE];
-  double retParam = 0.;
-
-  for (int i = 0; i < params->paramCount; i++) {
-    if (!strcmp(name, params->data[i].name)) {
-      /* Check that exactly one number is parsed. */
-      if (sscanf(params->data[i].value, "%lf%s", &retParam, str) != 1) {
-        error(
-            "Tried parsing double '%s' but found '%s' with illegal double "
-            "characters '%s'.",
-            params->data[i].name, params->data[i].value, str);
-      }
+int parser_get_opt_param_char_array(struct swift_params *params,
+                                    const char *name, int nval, char *values) {
+  if (get_param_char_array(params, name, 0, nval, values) != 1) {
+    save_param_char_array(params, name, nval, values);
+    return 0;
+  }
+  return 1;
+}
 
-      /* this parameter has been used */
-      params->data[i].used = 1;
+/**
+ * @brief Retrieve int array parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param nval number of values expected.
+ * @param values Values of the parameter found, of size at least nvals.
+ */
+void parser_get_param_int_array(struct swift_params *params, const char *name,
+                                int nval, int *values) {
+  get_param_int_array(params, name, 1, nval, values);
+}
 
-      return retParam;
-    }
+/**
+ * @brief Retrieve optional int array parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param nval number of values expected.
+ * @param values Values of the parameter found, of size at least nvals. If the
+ *               parameter is not found these values will be returned
+ *               unmodified, so should be set to the default values.
+ * @return whether the parameter has been found.
+ */
+int parser_get_opt_param_int_array(struct swift_params *params,
+                                   const char *name, int nval, int *values) {
+  if (get_param_int_array(params, name, 0, nval, values) != 1) {
+    save_param_int_array(params, name, nval, values);
+    return 0;
   }
+  return 1;
+}
 
-  /* Generate string for new parameter */
-  sprintf(str, "%s: %lf", name, def);
+/**
+ * @brief Retrieve float array parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param nval number of values expected.
+ * @param values Values of the parameter found, of size at least nvals.
+ */
+void parser_get_param_float_array(struct swift_params *params, const char *name,
+                                  int nval, float *values) {
+  get_param_float_array(params, name, 1, nval, values);
+}
 
-  /* Add it to params */
-  parser_set_param(params, str);
+/**
+ * @brief Retrieve optional float array parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param nval number of values expected.
+ * @param values Values of the parameter found, of size at least nvals. If the
+ *               parameter is not found these values will be returned
+ *               unmodified, so should be set to the default values.
+ * @return whether the parameter has been found.
+ */
+int parser_get_opt_param_float_array(struct swift_params *params,
+                                     const char *name, int nval,
+                                     float *values) {
+  if (get_param_float_array(params, name, 0, nval, values) != 1) {
+    save_param_float_array(params, name, nval, values);
+    return 0;
+  }
+  return 1;
+}
 
-  /* Set parameter as used */
-  params->data[params->paramCount - 1].used = 1;
+/**
+ * @brief Retrieve double array parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param nval number of values expected.
+ * @param values Values of the parameter found, of size at least nvals.
+ */
+void parser_get_param_double_array(struct swift_params *params,
+                                   const char *name, int nval, double *values) {
+  get_param_double_array(params, name, 1, nval, values);
+}
 
-  return def;
+/**
+ * @brief Retrieve optional double array parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param nval number of values expected.
+ * @param values Values of the parameter found, of size at least nvals. If the
+ *               parameter is not found these values will be returned
+ *               unmodified, so should be set to the default values.
+ * @return whether the parameter has been found.
+ */
+int parser_get_opt_param_double_array(struct swift_params *params,
+                                      const char *name, int nval,
+                                      double *values) {
+  if (get_param_double_array(params, name, 0, nval, values) != 1) {
+    save_param_double_array(params, name, nval, values);
+    return 0;
+  }
+  return 1;
 }
 
 /**
- * @brief Retrieve string parameter from structure.
+ * @brief Retrieve string array parameter from structure.
  *
  * @param params Structure that holds the parameters
  * @param name Name of the parameter to be found
- * @param def Default value of the parameter of not found.
- * @param retParam (return) Value of the parameter found
+ * @param required whether the parameter is required or not.
+ * @param nval number of values located.
+ * @param values pointer to an array of [nval] pointers to the strings.
+ *        Note this must be freed by a call to
+ *        parser_free_param_string_array when no longer required.
+ * @result whether the parameter was found or not. Note if required
+ *        an error will be thrown.
  */
-void parser_get_opt_param_string(struct swift_params *params, const char *name,
-                                 char *retParam, const char *def) {
+static int get_string_array(struct swift_params *params, const char *name,
+                            int required, int *nval, char ***values) {
+
+  char cpy[PARSER_MAX_LINE_SIZE];
+  *nval = 0;
 
   for (int i = 0; i < params->paramCount; i++) {
     if (!strcmp(name, params->data[i].name)) {
-      strcpy(retParam, params->data[i].value);
+      char *cp = cpy;
+      strcpy(cp, params->data[i].value);
+      cp = trim_both(cp);
 
-      /* this parameter has been used */
-      params->data[i].used = 1;
+      /* Strip off [], if present. */
+      if (cp[0] == '[') cp++;
+      int l = strlen(cp);
+      if (cp[l - 1] == ']') cp[l - 1] = '\0';
+      cp = trim_both(cp);
 
-      return;
+      *nval = parse_quoted_strings(cp, values);
+
+      params->data[i].used = 1;
+      return 1;
     }
   }
+  if (required)
+    error("Cannot find '%s' in the structure, in file '%s'.", name,
+          params->fileName);
+  return 0;
+}
 
-  /* Generate string for new parameter */
-  char str[PARSER_MAX_LINE_SIZE];
-  sprintf(str, "%s: %s", name, def);
+/**
+ * @brief Retrieve string array parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param nval number of values located.
+ * @param values pointer to an array of [nval] pointers to the strings.
+ *        Note this must be freed by a call to
+ *        parser_free_param_string_array when no longer required.
+ */
+void parser_get_param_string_array(struct swift_params *params,
+                                   const char *name, int *nval,
+                                   char ***values) {
+  get_string_array(params, name, 1, nval, values);
+}
 
-  /* Add it to params */
-  parser_set_param(params, str);
+/**
+ * @brief Retrieve optional string array parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param nval number of values located.
+ * @param values pointer to an array of [nval] pointers to the strings.
+ *        Note this must be freed by a call to
+ *        parser_free_param_string_array when no longer required.
+ * @param ndef the number of default values.
+ * @param def the default values as an array of pointers to strings.
+ *        Note copied to values if used.
+ * @result whether the parameter was found or not.
+ */
+int parser_get_opt_param_string_array(struct swift_params *params,
+                                      const char *name, int *nval,
+                                      char ***values, int ndef,
+                                      const char *def[]) {
+  if (get_string_array(params, name, 0, nval, values) == 1) return 1;
+
+  /* Not found, so save the default values against the parameter. Look for
+   * single quotes in value and use if not found, otherwise use double
+   * quotes. We don't support having both in a string. */
+  char cpy[PARSER_MAX_LINE_SIZE];
+  int k = sprintf(cpy, "%s: [", name);
+  int i = 0;
+  for (i = 0; i < ndef - 1; i++) {
+    if (strchr(def[i], '\'') == 0)
+      k += sprintf(&cpy[k], "'%s', ", def[i]);
+    else
+      k += sprintf(&cpy[k], "\"%s\", ", def[i]);
+  }
+  if (strchr(def[i], '\'') == 0)
+    sprintf(&cpy[k], "'%s']", def[i]);
+  else
+    sprintf(&cpy[k], "\"%s\"]", def[i]);
+  parser_set_param(params, cpy);
+
+  /* Now copy to output space. */
+  char **strings;
+  strings = (char **)malloc(ndef * sizeof(char *));
+  for (int j = 0; j < ndef; j++) {
+    strings[j] = (char *)malloc((strlen(def[j]) + 1) * sizeof(char));
+    strcpy(strings[j], def[j]);
+  }
+  *values = strings;
+  *nval = ndef;
 
-  /* Set parameter as used */
   params->data[params->paramCount - 1].used = 1;
+  return 0;
+}
 
-  strcpy(retParam, def);
+/**
+ * @brief Free string array allocated by parser_get_param_string_array.
+ *
+ * @param nval number of strings returned.
+ * @param values pointer to the returned values.
+ */
+void parser_free_param_string_array(int nval, char **values) {
+  for (int i = 0; i < nval; i++) {
+    free(values[i]);
+  }
+  free(values);
+  return;
 }
 
 /**
@@ -830,51 +1098,69 @@ void parser_print_params(const struct swift_params *params) {
 void parser_write_params_to_file(const struct swift_params *params,
                                  const char *file_name, int write_used) {
   FILE *file = fopen(file_name, "w");
-  char section[PARSER_MAX_LINE_SIZE] = {0};
   char param_name[PARSER_MAX_LINE_SIZE] = {0};
+  char section[PARSER_MAX_LINE_SIZE] = {0};
   char *token;
 
   /* Start of file identifier in YAML. */
   fprintf(file, "%s\n", PARSER_START_OF_FILE);
 
-  for (int i = 0; i < params->paramCount; i++) {
-    if (write_used && !params->data[i].used) {
-#ifdef SWIFT_DEBUG_CHECKS
-      message(
-          "Parameter `%s` was not used. "
-          "Only the parameter used are written.",
-          params->data[i].name);
-#endif
-      continue;
-    } else if (!write_used && params->data[i].used)
-      continue;
-    /* Check that the parameter name contains a section name. */
-    if (strchr(params->data[i].name, PARSER_VALUE_CHAR)) {
-      /* Copy the parameter name into a temporary string and find the section
-       * name. */
-      strcpy(param_name, params->data[i].name);
-      token = strtok(param_name, PARSER_VALUE_STRING);
-
-      /* If a new section name is found print it to the file. */
-      if (strcmp(token, section)) {
-        strcpy(section, token);
-        fprintf(file, "\n%s%c\n", section, PARSER_VALUE_CHAR);
+  /* Flags to track which parameters are written. */
+  int *written = (int *)calloc(params->paramCount, sizeof(int));
+  int nwritten = 0;
+
+  /* Loop over all sections. These are not contiguous when storing optional
+   * values. */
+  for (int k = 0; k < params->sectionCount; k++) {
+    int first = 1;
+
+    /* Locate parameters in this section. */
+    for (int i = 0; i < params->paramCount; i++) {
+      if (!written[i] && ((write_used && params->data[i].used) ||
+                          (!write_used && !params->data[i].used))) {
+
+        /* Find section part of name, if have one. */
+        if (strchr(params->data[i].name, PARSER_VALUE_CHAR)) {
+          strcpy(param_name, params->data[i].name);
+          token = strtok(param_name, PARSER_VALUE_STRING);
+          strcpy(section, token);
+          strcat(section, PARSER_VALUE_STRING);
+
+          /* If in our section name print it to the file. */
+          if (strcmp(section, params->section[k].name) == 0) {
+            if (first) {
+              fprintf(file, "\n%s\n", section);
+              first = 0;
+            }
+
+            /* Remove white space from parameter name and write it to the
+             * file. */
+            token = trim_both(strtok(NULL, "#\n"));
+            fprintf(file, "  %s%c %s\n", token, PARSER_VALUE_CHAR,
+                    params->data[i].value);
+            written[i] = 1;
+            nwritten++;
+          }
+        }
       }
+    }
+  }
 
-      /* Remove white space from parameter name and write it to the file. */
-      token = strtok(NULL, " #\n");
-
-      fprintf(file, "  %s%c %s\n", token, PARSER_VALUE_CHAR,
-              params->data[i].value);
-    } else {
-      fprintf(file, "\n%s%c %s\n", params->data[i].name, PARSER_VALUE_CHAR,
-              params->data[i].value);
+  /* Write out any parameters outside of sections. */
+  if (nwritten < params->paramCount) {
+    for (int i = 0; i < params->paramCount; i++) {
+      if (!written[i] && ((write_used && params->data[i].used) ||
+                          (!write_used && !params->data[i].used))) {
+        fprintf(file, "\n%s%c %s\n", params->data[i].name, PARSER_VALUE_CHAR,
+                params->data[i].value);
+      }
     }
   }
 
   /* End of file identifier in YAML. */
   fprintf(file, "%s\n", PARSER_END_OF_FILE);
 
+  free(written);
   fclose(file);
 }
 
diff --git a/src/parser.h b/src/parser.h
index 2b06fce03cbfe2d6a73bc28960bc83bc822e4459..b4a8dd3d32310601f86db7af46d1bb6244c96f0a 100644
--- a/src/parser.h
+++ b/src/parser.h
@@ -80,6 +80,32 @@ double parser_get_opt_param_double(struct swift_params *params,
                                    const char *name, double def);
 void parser_get_opt_param_string(struct swift_params *params, const char *name,
                                  char *retParam, const char *def);
+void parser_get_param_char_array(struct swift_params *params, const char *name,
+                                 int nval, char *values);
+void parser_get_param_int_array(struct swift_params *params, const char *name,
+                                int nval, int *values);
+void parser_get_param_float_array(struct swift_params *params, const char *name,
+                                  int nval, float *values);
+void parser_get_param_double_array(struct swift_params *params,
+                                   const char *name, int nval, double *values);
+
+int parser_get_opt_param_char_array(struct swift_params *params,
+                                    const char *name, int nval, char *values);
+int parser_get_opt_param_int_array(struct swift_params *params,
+                                   const char *name, int nval, int *values);
+int parser_get_opt_param_float_array(struct swift_params *params,
+                                     const char *name, int nval, float *values);
+int parser_get_opt_param_double_array(struct swift_params *params,
+                                      const char *name, int nval,
+                                      double *values);
+
+void parser_get_param_string_array(struct swift_params *params,
+                                   const char *name, int *nval, char ***values);
+int parser_get_opt_param_string_array(struct swift_params *params,
+                                      const char *name, int *nval,
+                                      char ***values, int ndef,
+                                      const char *def[]);
+void parser_free_param_string_array(int nval, char **values);
 
 #if defined(HAVE_HDF5)
 void parser_write_params_to_hdf5(const struct swift_params *params, hid_t grp,
diff --git a/src/partition.c b/src/partition.c
index 85a51dddf2797e7d203da95abc42639c29f11aa6..25b073c3198a5e0a37e1d39709e8e44bc9fc9473 100644
--- a/src/partition.c
+++ b/src/partition.c
@@ -1095,12 +1095,8 @@ void partition_init(struct partition *partition,
 
   /* In case of grid, read more parameters */
   if (part_type[0] == 'g') {
-    partition->grid[0] = parser_get_opt_param_int(
-        params, "DomainDecomposition:initial_grid_x", partition->grid[0]);
-    partition->grid[1] = parser_get_opt_param_int(
-        params, "DomainDecomposition:initial_grid_y", partition->grid[1]);
-    partition->grid[2] = parser_get_opt_param_int(
-        params, "DomainDecomposition:initial_grid_z", partition->grid[2]);
+    parser_get_opt_param_int_array(params, "DomainDecomposition:initial_grid",
+                                   3, partition->grid);
   }
 
   /* Now let's check what the user wants as a repartition strategy */
diff --git a/src/potential/isothermal/potential.h b/src/potential/isothermal/potential.h
index 4267c9becc7251ffe276b69f078e374504c22aab..b5f8d7c39738bfe1895c73e6e59ae1279c0f74fa 100644
--- a/src/potential/isothermal/potential.h
+++ b/src/potential/isothermal/potential.h
@@ -42,7 +42,7 @@
 struct external_potential {
 
   /*! Position of the centre of potential */
-  double x, y, z;
+  double x[3];
 
   /*! Rotation velocity */
   double vrot;
@@ -73,9 +73,9 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
     const struct phys_const* restrict phys_const,
     const struct gpart* restrict g) {
 
-  const float dx = g->x[0] - potential->x;
-  const float dy = g->x[1] - potential->y;
-  const float dz = g->x[2] - potential->z;
+  const float dx = g->x[0] - potential->x[0];
+  const float dy = g->x[1] - potential->x[1];
+  const float dz = g->x[2] - potential->x[2];
 
   const float r2_plus_epsilon2_inv =
       1.f / (dx * dx + dy * dy + dz * dz + potential->epsilon2);
@@ -115,9 +115,9 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
     double time, const struct external_potential* potential,
     const struct phys_const* const phys_const, struct gpart* g) {
 
-  const float dx = g->x[0] - potential->x;
-  const float dy = g->x[1] - potential->y;
-  const float dz = g->x[2] - potential->z;
+  const float dx = g->x[0] - potential->x[0];
+  const float dy = g->x[1] - potential->x[1];
+  const float dz = g->x[2] - potential->x[2];
   const float r2_plus_epsilon2_inv =
       1.f / (dx * dx + dy * dy + dz * dz + potential->epsilon2);
 
@@ -144,9 +144,9 @@ external_gravity_get_potential_energy(
     double time, const struct external_potential* potential,
     const struct phys_const* const phys_const, const struct gpart* g) {
 
-  const float dx = g->x[0] - potential->x;
-  const float dy = g->x[1] - potential->y;
-  const float dz = g->x[2] - potential->z;
+  const float dx = g->x[0] - potential->x[0];
+  const float dy = g->x[1] - potential->x[1];
+  const float dz = g->x[2] - potential->x[2];
 
   return -0.5f * potential->vrot * potential->vrot *
          logf(dx * dx + dy * dy + dz * dz + potential->epsilon2);
@@ -166,15 +166,12 @@ static INLINE void potential_init_backend(
     const struct unit_system* us, const struct space* s,
     struct external_potential* potential) {
 
-  potential->x =
-      s->dim[0] / 2. +
-      parser_get_param_double(parameter_file, "IsothermalPotential:position_x");
-  potential->y =
-      s->dim[1] / 2. +
-      parser_get_param_double(parameter_file, "IsothermalPotential:position_y");
-  potential->z =
-      s->dim[2] / 2. +
-      parser_get_param_double(parameter_file, "IsothermalPotential:position_z");
+  parser_get_param_double_array(parameter_file, "IsothermalPotential:position",
+                                3, potential->x);
+  potential->x[0] += s->dim[0] / 2.;
+  potential->x[1] += s->dim[1] / 2.;
+  potential->x[2] += s->dim[2] / 2.;
+
   potential->vrot =
       parser_get_param_double(parameter_file, "IsothermalPotential:vrot");
   potential->timestep_mult = parser_get_param_float(
@@ -198,7 +195,7 @@ static INLINE void potential_print_backend(
       "External potential is 'Isothermal' with properties are (x,y,z) = (%e, "
       "%e, %e), vrot = %e "
       "timestep multiplier = %e, epsilon = %e",
-      potential->x, potential->y, potential->z, potential->vrot,
+      potential->x[0], potential->x[1], potential->x[2], potential->vrot,
       potential->timestep_mult, sqrtf(potential->epsilon2));
 }
 
diff --git a/src/potential/point_mass/potential.h b/src/potential/point_mass/potential.h
index db875842d51f6c0a28dd1308fd7dd1728e746ce4..f9d56a1ff165f2331c91ea828b5ffe0e0db76c2f 100644
--- a/src/potential/point_mass/potential.h
+++ b/src/potential/point_mass/potential.h
@@ -41,7 +41,7 @@
 struct external_potential {
 
   /*! Position of the point mass */
-  double x, y, z;
+  double x[3];
 
   /*! Mass */
   double mass;
@@ -66,15 +66,15 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
     const struct gpart* restrict g) {
 
   const float G_newton = phys_const->const_newton_G;
-  const float dx = g->x[0] - potential->x;
-  const float dy = g->x[1] - potential->y;
-  const float dz = g->x[2] - potential->z;
+  const float dx = g->x[0] - potential->x[0];
+  const float dy = g->x[1] - potential->x[1];
+  const float dz = g->x[2] - potential->x[2];
   const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
   const float rinv2 = rinv * rinv;
   const float rinv3 = rinv2 * rinv;
-  const float drdv = (g->x[0] - potential->x) * (g->v_full[0]) +
-                     (g->x[1] - potential->y) * (g->v_full[1]) +
-                     (g->x[2] - potential->z) * (g->v_full[2]);
+  const float drdv = (g->x[0] - potential->x[0]) * (g->v_full[0]) +
+                     (g->x[1] - potential->x[1]) * (g->v_full[1]) +
+                     (g->x[2] - potential->x[2]) * (g->v_full[2]);
   const float dota_x = G_newton * potential->mass * rinv3 *
                        (-g->v_full[0] + 3.f * rinv2 * drdv * dx);
   const float dota_y = G_newton * potential->mass * rinv3 *
@@ -109,9 +109,9 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
     double time, const struct external_potential* restrict potential,
     const struct phys_const* restrict phys_const, struct gpart* restrict g) {
 
-  const float dx = g->x[0] - potential->x;
-  const float dy = g->x[1] - potential->y;
-  const float dz = g->x[2] - potential->z;
+  const float dx = g->x[0] - potential->x[0];
+  const float dy = g->x[1] - potential->x[1];
+  const float dz = g->x[2] - potential->x[2];
   const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
   const float rinv3 = rinv * rinv * rinv;
 
@@ -134,9 +134,9 @@ external_gravity_get_potential_energy(
     double time, const struct external_potential* potential,
     const struct phys_const* const phys_const, const struct gpart* g) {
 
-  const float dx = g->x[0] - potential->x;
-  const float dy = g->x[1] - potential->y;
-  const float dz = g->x[2] - potential->z;
+  const float dx = g->x[0] - potential->x[0];
+  const float dy = g->x[1] - potential->x[1];
+  const float dz = g->x[2] - potential->x[2];
   const float rinv = 1. / sqrtf(dx * dx + dy * dy + dz * dz);
   return -phys_const->const_newton_G * potential->mass * rinv;
 }
@@ -156,12 +156,8 @@ static INLINE void potential_init_backend(
     const struct unit_system* us, const struct space* s,
     struct external_potential* potential) {
 
-  potential->x =
-      parser_get_param_double(parameter_file, "PointMassPotential:position_x");
-  potential->y =
-      parser_get_param_double(parameter_file, "PointMassPotential:position_y");
-  potential->z =
-      parser_get_param_double(parameter_file, "PointMassPotential:position_z");
+  parser_get_param_double_array(parameter_file, "PointMassPotential:position",
+                                3, potential->x);
   potential->mass =
       parser_get_param_double(parameter_file, "PointMassPotential:mass");
   potential->timestep_mult = parser_get_param_float(
@@ -179,7 +175,7 @@ static INLINE void potential_print_backend(
   message(
       "External potential is 'Point mass' with properties (x,y,z) = (%e, %e, "
       "%e), M = %e timestep multiplier = %e.",
-      potential->x, potential->y, potential->z, potential->mass,
+      potential->x[0], potential->x[1], potential->x[2], potential->mass,
       potential->timestep_mult);
 }
 
diff --git a/src/potential/point_mass_ring/potential.h b/src/potential/point_mass_ring/potential.h
index 551efe32521a5c5ee8068ba409dbb81547103e8f..b384c3df6771c9f3342f5aed6514e75c77bf2c17 100644
--- a/src/potential/point_mass_ring/potential.h
+++ b/src/potential/point_mass_ring/potential.h
@@ -41,7 +41,7 @@
 struct external_potential {
 
   /*! Position of the point mass */
-  double x, y, z;
+  double x[3];
 
   /*! Mass */
   double mass;
@@ -68,16 +68,16 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
     const struct gpart* restrict g) {
 
   const float G_newton = phys_const->const_newton_G;
-  const float dx = g->x[0] - potential->x;
-  const float dy = g->x[1] - potential->y;
-  const float dz = g->x[2] - potential->z;
+  const float dx = g->x[0] - potential->x[0];
+  const float dy = g->x[1] - potential->x[1];
+  const float dz = g->x[2] - potential->x[2];
   const float r = sqrtf(dx * dx + dy * dy + dz * dz);
   const float rinv = 1.f / r;
   const float rinv2 = rinv * rinv;
   const float rinv3 = rinv2 * rinv;
-  const float drdv = (g->x[0] - potential->x) * (g->v_full[0]) +
-                     (g->x[1] - potential->y) * (g->v_full[1]) +
-                     (g->x[2] - potential->z) * (g->v_full[2]);
+  const float drdv = (g->x[0] - potential->x[0]) * (g->v_full[0]) +
+                     (g->x[1] - potential->x[1]) * (g->v_full[1]) +
+                     (g->x[2] - potential->x[2]) * (g->v_full[2]);
   float factor;
 
   if (r < 0.175) {
@@ -129,9 +129,9 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
     double time, const struct external_potential* restrict potential,
     const struct phys_const* restrict phys_const, struct gpart* restrict g) {
 
-  const float dx = g->x[0] - potential->x;
-  const float dy = g->x[1] - potential->y;
-  const float dz = g->x[2] - potential->z;
+  const float dx = g->x[0] - potential->x[0];
+  const float dy = g->x[1] - potential->x[1];
+  const float dz = g->x[2] - potential->x[2];
   const float r = sqrtf(dx * dx + dy * dy + dz * dz);
   const float rinv = 1.f / r;
   const float rinv2 = rinv * rinv;
@@ -174,9 +174,9 @@ external_gravity_get_potential_energy(
     double time, const struct external_potential* potential,
     const struct phys_const* const phys_const, const struct gpart* g) {
 
-  const float dx = g->x[0] - potential->x;
-  const float dy = g->x[1] - potential->y;
-  const float dz = g->x[2] - potential->z;
+  const float dx = g->x[0] - potential->x[0];
+  const float dy = g->x[1] - potential->x[1];
+  const float dz = g->x[2] - potential->x[2];
   const float rinv = 1. / sqrtf(dx * dx + dy * dy + dz * dz);
   return -phys_const->const_newton_G * potential->mass * rinv;
 }
@@ -196,12 +196,8 @@ static INLINE void potential_init_backend(
     const struct unit_system* us, const struct space* s,
     struct external_potential* potential) {
 
-  potential->x =
-      parser_get_param_double(parameter_file, "PointMassPotential:position_x");
-  potential->y =
-      parser_get_param_double(parameter_file, "PointMassPotential:position_y");
-  potential->z =
-      parser_get_param_double(parameter_file, "PointMassPotential:position_z");
+  parser_get_param_double_array(parameter_file, "PointMassPotential:position",
+                                3, potential->x);
   potential->mass =
       parser_get_param_double(parameter_file, "PointMassPotential:mass");
   potential->timestep_mult = parser_get_param_float(
@@ -219,7 +215,7 @@ static INLINE void potential_print_backend(
   message(
       "External potential is 'Point mass' with properties (x,y,z) = (%e, %e, "
       "%e), M = %e timestep multiplier = %e.",
-      potential->x, potential->y, potential->z, potential->mass,
+      potential->x[0], potential->x[1], potential->x[2], potential->mass,
       potential->timestep_mult);
 }
 
diff --git a/src/potential/point_mass_softened/potential.h b/src/potential/point_mass_softened/potential.h
index 80959ec923cbedcbea5fba3293c8ae4f94f65679..0e35e7bb9870c7954b47316a3cc30bb68cde5fc4 100644
--- a/src/potential/point_mass_softened/potential.h
+++ b/src/potential/point_mass_softened/potential.h
@@ -50,7 +50,7 @@
 struct external_potential {
 
   /*! Position of the point mass */
-  double x, y, z;
+  double x[3];
 
   /*! Mass */
   double mass;
@@ -77,9 +77,9 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
     const struct phys_const* restrict phys_const,
     const struct gpart* restrict g) {
 
-  const float dx = g->x[0] - potential->x;
-  const float dy = g->x[1] - potential->y;
-  const float dz = g->x[2] - potential->z;
+  const float dx = g->x[0] - potential->x[0];
+  const float dy = g->x[1] - potential->x[1];
+  const float dz = g->x[2] - potential->x[2];
 
   const float softening2 = potential->softening * potential->softening;
   const float r2 = dx * dx + dy * dy + dz * dz;
@@ -133,9 +133,9 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
     double time, const struct external_potential* restrict potential,
     const struct phys_const* restrict phys_const, struct gpart* restrict g) {
 
-  const float dx = g->x[0] - potential->x;
-  const float dy = g->x[1] - potential->y;
-  const float dz = g->x[2] - potential->z;
+  const float dx = g->x[0] - potential->x[0];
+  const float dy = g->x[1] - potential->x[1];
+  const float dz = g->x[2] - potential->x[2];
   const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz +
                                  potential->softening * potential->softening);
   const float rinv3 = rinv * rinv * rinv;
@@ -159,9 +159,9 @@ external_gravity_get_potential_energy(
     double time, const struct external_potential* potential,
     const struct phys_const* const phys_const, const struct gpart* g) {
 
-  const float dx = g->x[0] - potential->x;
-  const float dy = g->x[1] - potential->y;
-  const float dz = g->x[2] - potential->z;
+  const float dx = g->x[0] - potential->x[0];
+  const float dy = g->x[1] - potential->x[1];
+  const float dz = g->x[2] - potential->x[2];
   const float rinv = 1. / sqrtf(dx * dx + dy * dy + dz * dz +
                                 potential->softening * potential->softening);
 
@@ -183,12 +183,8 @@ static INLINE void potential_init_backend(
     const struct unit_system* us, const struct space* s,
     struct external_potential* potential) {
 
-  potential->x =
-      parser_get_param_double(parameter_file, "PointMassPotential:position_x");
-  potential->y =
-      parser_get_param_double(parameter_file, "PointMassPotential:position_y");
-  potential->z =
-      parser_get_param_double(parameter_file, "PointMassPotential:position_z");
+  parser_get_param_double_array(parameter_file, "PointMassPotential:position",
+                                3, potential->x);
   potential->mass =
       parser_get_param_double(parameter_file, "PointMassPotential:mass");
   potential->timestep_mult = parser_get_param_float(
@@ -208,7 +204,7 @@ static INLINE void potential_print_backend(
   message(
       "External potential is 'Point mass' with properties (x,y,z) = (%e, %e, "
       "%e), M = %e timestep multiplier = %e, softening = %e.",
-      potential->x, potential->y, potential->z, potential->mass,
+      potential->x[0], potential->x[1], potential->x[2], potential->mass,
       potential->timestep_mult, potential->softening);
 }
 
diff --git a/src/serial_io.c b/src/serial_io.c
index 83bcd88c935603ccbebaee9d9065a19877d53922..dafa75ab0baacb1b5ddeee34020c9773893bced7 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -518,7 +518,7 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
 
     /* Read the unit system used in the ICs */
     if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
-    io_read_unit_system(h_file, ic_units, mpi_rank);
+    io_read_unit_system(h_file, ic_units, internal_units, mpi_rank);
 
     if (units_are_equal(ic_units, internal_units)) {
 
diff --git a/src/single_io.c b/src/single_io.c
index 9d6ace82369e18f4dc16b3e31844db3f3cca3fe6..be33bba4d7c95454892b02836f073c085e07b2ea 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -424,7 +424,7 @@ void read_ic_single(const char* fileName,
   struct unit_system* ic_units =
       (struct unit_system*)malloc(sizeof(struct unit_system));
   if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
-  io_read_unit_system(h_file, ic_units, 0);
+  io_read_unit_system(h_file, ic_units, internal_units, 0);
 
   /* Tell the user if a conversion will be needed */
   if (units_are_equal(ic_units, internal_units)) {
diff --git a/src/space.c b/src/space.c
index 5b9e592bd60b1cf49ad276eafa351d46765b8022..e80cc37190452013bdf4a29da7bd1ebae9fa2c2a 100644
--- a/src/space.c
+++ b/src/space.c
@@ -2771,12 +2771,8 @@ void space_init(struct space *s, struct swift_params *params,
 
   /* Apply shift */
   double shift[3] = {0.0, 0.0, 0.0};
-  shift[0] =
-      parser_get_opt_param_double(params, "InitialConditions:shift_x", 0.0);
-  shift[1] =
-      parser_get_opt_param_double(params, "InitialConditions:shift_y", 0.0);
-  shift[2] =
-      parser_get_opt_param_double(params, "InitialConditions:shift_z", 0.0);
+  parser_get_opt_param_double_array(params, "InitialConditions:shift", 3,
+                                    shift);
   if ((shift[0] != 0. || shift[1] != 0. || shift[2] != 0.) && !dry_run) {
     message("Shifting particles by [%e %e %e]", shift[0], shift[1], shift[2]);
     for (size_t k = 0; k < Npart; k++) {
diff --git a/src/tools.c b/src/tools.c
index 8b8b7fdf37d91547f328a6b49d69b3c5a491aed1..9c0df6012737872eef8d97521b3a7532ceb42720 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -23,6 +23,7 @@
 #include "../config.h"
 
 /* Some standard headers. */
+#include <ctype.h>
 #include <math.h>
 #include <stddef.h>
 #include <stdio.h>
@@ -728,3 +729,46 @@ long get_maxrss() {
   getrusage(RUSAGE_SELF, &usage);
   return usage.ru_maxrss;
 }
+
+/**
+ * @brief trim leading white space from a string.
+ *
+ * Returns pointer to first character.
+ *
+ * @param s the string.
+ * @result the result.
+ */
+char *trim_leading(char *s) {
+  if (s == NULL || strlen(s) < 2) return s;
+  while (isspace(*s)) s++;
+  return s;
+}
+
+/**
+ * @brief trim trailing white space from a string.
+ *
+ * Modifies the string by adding a NULL to the end.
+ *
+ * @param s the string.
+ * @result the result.
+ */
+char *trim_trailing(char *s) {
+  if (s == NULL || strlen(s) < 2) return s;
+  char *end = s + strlen(s) - 1;
+  while (isspace(*end)) end--;
+  *(end + 1) = '\0';
+  return s;
+}
+
+/**
+ * @brief trim leading and trailing white space from a string.
+ *
+ * Can modify the string by adding a NULL to the end.
+ *
+ * @param s the string.
+ * @result the result.
+ */
+char *trim_both(char *s) {
+  if (s == NULL || strlen(s) < 2) return s;
+  return trim_trailing(trim_leading(s));
+}
diff --git a/src/tools.h b/src/tools.h
index a54510000d3c2843e8d60047752b19a46bd502d9..25d024679174eabbe89908c0254651e4bbc69e15 100644
--- a/src/tools.h
+++ b/src/tools.h
@@ -54,4 +54,8 @@ int compare_particles(struct part a, struct part b, double threshold);
 
 long get_maxrss(void);
 
+char *trim_leading(char *s);
+char *trim_trailing(char *s);
+char *trim_both(char *s);
+
 #endif /* SWIFT_TOOL_H */
diff --git a/src/units.c b/src/units.c
index 48f0a3aee6e348b5df24ac41b308aebf6f70224a..6d070c9c8ebdf56dff53b9dfeecc73d68fb6b4b6 100644
--- a/src/units.c
+++ b/src/units.c
@@ -128,6 +128,21 @@ void units_init_default(struct unit_system* us, struct swift_params* params,
       parser_get_opt_param_double(params, buffer, def->UnitTemperature_in_cgs);
 }
 
+/**
+ * @brief Copy the content of a #unit_system to another one.
+ *
+ * @param dest The destination of the copy.
+ * @param src The source of the copy.
+ */
+void units_copy(struct unit_system* dest, const struct unit_system* src) {
+
+  dest->UnitMass_in_cgs = src->UnitMass_in_cgs;
+  dest->UnitLength_in_cgs = src->UnitLength_in_cgs;
+  dest->UnitTime_in_cgs = src->UnitTime_in_cgs;
+  dest->UnitCurrent_in_cgs = src->UnitCurrent_in_cgs;
+  dest->UnitTemperature_in_cgs = src->UnitTemperature_in_cgs;
+}
+
 /**
  * @brief Returns the base unit conversion factor for a given unit system
  * @param us The unit_system used
diff --git a/src/units.h b/src/units.h
index 829a1ce542500308cbc64a2463545fbd23921eef..da2c209815b07d1d5597a598ee4a61f3132e39db 100644
--- a/src/units.h
+++ b/src/units.h
@@ -103,6 +103,7 @@ void units_init_from_params(struct unit_system*, struct swift_params*,
 void units_init_default(struct unit_system* us, struct swift_params* params,
                         const char* category, const struct unit_system* def);
 
+void units_copy(struct unit_system* dest, const struct unit_system* src);
 int units_are_equal(const struct unit_system* a, const struct unit_system* b);
 
 /* Base units */
diff --git a/tests/Makefile.am b/tests/Makefile.am
index e8f2a4f0ae39254c28e981fef5e36866e56395d8..ad7dcf92d6e4133e80bdd077447841626e59d722 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -27,7 +27,8 @@ TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetr
         testMatrixInversion testThreadpool testDump testLogger testInteractions.sh \
         testVoronoi1D testVoronoi2D testVoronoi3D testGravityDerivatives \
 	testPeriodicBC.sh testPeriodicBCPerturbed.sh testPotentialSelf \
-	testPotentialPair testEOS testUtilities testSelectOutput.sh
+	testPotentialPair testEOS testUtilities testSelectOutput.sh \
+    testCbrt
 
 # List of test programs to compile
 check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
@@ -38,7 +39,7 @@ check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
                  testRiemannHLLC testMatrixInversion testDump testLogger \
 		 testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC \
 		 testGravityDerivatives testPotentialSelf testPotentialPair testEOS testUtilities \
-		 testSelectOutput
+		 testSelectOutput testCbrt
 
 # Rebuild tests when SWIFT is updated.
 $(check_PROGRAMS): ../src/.libs/libswiftsim.a
diff --git a/tests/testCbrt.c b/tests/testCbrt.c
new file mode 100644
index 0000000000000000000000000000000000000000..b608f9a00d619570c298f4123038f930584a245c
--- /dev/null
+++ b/tests/testCbrt.c
@@ -0,0 +1,129 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2016 Pedro Gonnet (pedro.gonnet@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"
+
+/* Standard includes. */
+#include <fenv.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+/* Local includes. */
+#include "cbrt.h"
+#include "clocks.h"
+#include "error.h"
+
+int main(int argc, char *argv[]) {
+
+  /* Initialize CPU frequency, this also starts time. */
+  unsigned long long cpufreq = 0;
+  clocks_set_cpufreq(cpufreq);
+
+  /* Choke on FP-exceptions */
+  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
+
+  /* Some constants for this test. */
+  const int num_vals = 200000000;
+  const float range_min = -1e6f;
+  const float range_max = 1e6f;
+  const float err_rel_tol = 1e-7;
+  message("executing %i runs of each command.", num_vals);
+
+  /* Create and fill an array of floats. */
+  float *data = (float *)malloc(sizeof(float) * num_vals);
+  for (int k = 0; k < num_vals; k++) {
+    data[k] = (float)rand() / RAND_MAX;
+    data[k] = (1.0f - data[k]) * range_min + data[k] * range_max;
+    if (data[k] == 0.f) k--; /* Skip 0 to avoid spurious mistakes */
+  }
+
+  /* First run just checks for correctness. */
+  for (int k = 0; k < num_vals; k++) {
+    const double exact = cbrt(data[k]);  // computed in double just to be sure.
+    const float ours = 1.0f / icbrtf(data[k]);
+    const float err_abs = fabsf(exact - ours);
+    const float err_rel = 0.5f * fabsf(exact - ours) / (exact + ours);
+
+    if (err_rel > err_rel_tol && data[k] != 0.f)
+      error(
+          "failed for x = %.8e, exact = %.8e, ours = %.8e, err = %.3e, rel = "
+          "%.3e",
+          data[k], exact, ours, err_abs, err_rel);
+  }
+
+  /* Second run to check the speed of the inverse cube root. */
+  double acc_exact = 0.0f;
+  ticks tic_exact = getticks();
+  for (int k = 0; k < num_vals; k++) {
+    acc_exact += 1.0f / cbrtf(data[k]);
+  }
+  message("1.0f / cbrtf took %9.3f %s (acc = %18.11e).",
+          clocks_from_ticks(getticks() - tic_exact), clocks_getunit(),
+          acc_exact);
+
+  double acc_ours = 0.0;
+  ticks tic_ours = getticks();
+  for (int k = 0; k < num_vals; k++) {
+    acc_ours += icbrtf(data[k]);
+  }
+  message("icbrtf       took %9.3f %s (acc = %18.11e).",
+          clocks_from_ticks(getticks() - tic_ours), clocks_getunit(), acc_ours);
+
+  /* Third run to check the speed of the cube root. */
+  acc_exact = 0.0f;
+  tic_exact = getticks();
+  for (int k = 0; k < num_vals; k++) {
+    acc_exact += cbrtf(data[k]);
+  }
+  message("cbrtf        took %9.3f %s (acc = %18.11e).",
+          clocks_from_ticks(getticks() - tic_exact), clocks_getunit(),
+          acc_exact);
+
+  acc_ours = 0.0f;
+  tic_ours = getticks();
+  for (int k = 0; k < num_vals; k++) {
+    const float temp = icbrtf(data[k]);
+    acc_ours += data[k] * temp * temp;
+  }
+  message("x * icbrtf^2 took %9.3f %s (acc = %18.11e).",
+          clocks_from_ticks(getticks() - tic_ours), clocks_getunit(), acc_ours);
+
+  /* Fourth run to check the speed of (.)^(2/3). */
+  acc_exact = 0.0f;
+  tic_exact = getticks();
+  for (int k = 0; k < num_vals; k++) {
+    const float temp = cbrtf(data[k]);
+    acc_exact += temp * temp;
+  }
+  message("cbrtf^2      took %9.3f %s (acc = %18.11e).",
+          clocks_from_ticks(getticks() - tic_exact), clocks_getunit(),
+          acc_exact);
+
+  acc_ours = 0.0f;
+  tic_ours = getticks();
+  for (int k = 0; k < num_vals; k++) {
+    acc_ours += data[k] * icbrtf(data[k]);
+  }
+  message("x * icbrtf   took %9.3f %s (acc = %18.11e).",
+          clocks_from_ticks(getticks() - tic_ours), clocks_getunit(), acc_ours);
+
+  return 0;
+}
diff --git a/tests/testParser.c b/tests/testParser.c
index c50e9b1a5473c1406e97894f211cc94aa7689d12..580ab873a2ffcf6825526b925610f220e8fa74ce 100644
--- a/tests/testParser.c
+++ b/tests/testParser.c
@@ -1,6 +1,7 @@
 /*******************************************************************************
  * This file is part of SWIFT.
  * Copyright (C) 2016 James Willis (james.s.willis@durham.ac.uk).
+ *               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
@@ -33,6 +34,7 @@ int main(int argc, char *argv[]) {
   parser_read_file(input_file, &param_file);
 
   /* Print the contents of the structure to stdout. */
+  printf("\n--- Values read from file:\n");
   parser_print_params(&param_file);
 
   /* Retrieve parameters and store them in variables defined above.
@@ -50,15 +52,93 @@ int main(int argc, char *argv[]) {
   const int optional =
       parser_get_opt_param_int(&param_file, "Simulation:optional", 1);
 
+  /* Optional things not mentioned in parameter file. Should be in output
+   * files.*/
+  parser_get_opt_param_int(&param_file, "Virtual:param1", 1);
+  parser_get_opt_param_int(&param_file, "Virtual:param2", 2);
+  parser_get_opt_param_int(&param_file, "Virtual:param3", 3);
+  parser_get_opt_param_int(&param_file, "Virtual:param4", 4);
+
   char ic_file[PARSER_MAX_LINE_SIZE];
   parser_get_param_string(&param_file, "IO:ic_file", ic_file);
 
+  char csides[3];
+  parser_get_param_char_array(&param_file, "Box:csides", 3, csides);
+
+  int isides[3];
+  parser_get_param_int_array(&param_file, "Box:isides", 3, isides);
+
+  float fsides[3];
+  parser_get_param_float_array(&param_file, "Box:fsides", 3, fsides);
+
+  double dsides[3];
+  parser_get_param_double_array(&param_file, "Box:dsides", 3, dsides);
+
+  int optsides[5] = {1, 2, 3, 4, 5};
+  int haveopt =
+      parser_get_opt_param_int_array(&param_file, "Box:moresides", 5, optsides);
+
+  char **var_result;
+  int nvar_result;
+  parser_get_param_string_array(&param_file, "Words:list", &nvar_result,
+                                &var_result);
+
+  printf("\nList from Words:list parameter\n");
+  for (int i = 0; i < nvar_result; i++) printf("   %d: %s\n", i, var_result[i]);
+
+  /* Get same list without []. */
+  char **var_result2;
+  int nvar_result2;
+  parser_get_param_string_array(&param_file, "Words:nonstdlist", &nvar_result2,
+                                &var_result2);
+
+  assert(nvar_result == nvar_result2);
+  for (int i = 0; i < nvar_result; i++)
+    assert(strcmp(var_result[i], var_result2[i]) == 0);
+
+  parser_free_param_string_array(nvar_result, var_result);
+  parser_free_param_string_array(nvar_result2, var_result2);
+
+  const char *optwords[4] = {"Word1", "Word2", "Word3", "Word4"};
+  int noptwords = 4;
+  int haveoptwords = parser_get_opt_param_string_array(
+      &param_file, "Simulation:optwords", &nvar_result, &var_result, noptwords,
+      optwords);
+  printf("\nList from Simulation:optwords parameter (%d of %d)\n", nvar_result,
+         noptwords);
+  assert(noptwords == nvar_result);
+  for (int i = 0; i < nvar_result; i++) {
+    assert(strcmp(optwords[i], var_result[i]) == 0);
+    printf("   %d: %s\n", i, var_result[i]);
+  }
+  parser_free_param_string_array(nvar_result, var_result);
+
+  /* Long list of values. */
+  parser_get_param_string_array(&param_file, "Words:long", &nvar_result,
+                                &var_result);
+  printf("No of letters in alphabet = %d\n", nvar_result);
+  assert(nvar_result == 26);
+  char alphabet[27];
+  for (int i = 0; i < nvar_result; i++) {
+    alphabet[i] = var_result[i][0];
+  }
+  alphabet[26] = '\0';
+  printf("The alphabet: %s\n", alphabet);
+  assert(strcmp("abcdefghijklmnopqrstuvwxyz", alphabet) == 0);
+  parser_free_param_string_array(nvar_result, var_result);
+
+  /* Print the contents of the structure to stdout now used. */
+  printf("\n--- Values after being used:\n");
+  parser_print_params(&param_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, kernel: %d optional: %d\n",
+      "\nValues read from file:\n"
+      "no_of_threads: %d, no_of_time_steps: %d, max_h: %f, start_time: %lf,"
+      " ic_file: %s, kernel: %d optional: %d\n",
       no_of_threads, no_of_time_steps, max_h, start_time, ic_file, kernel,
       optional);
+  printf("Box: [%d, %d, %d]\n", isides[0], isides[1], isides[2]);
 
   /* Print the contents of the structure to a file in YAML format. */
   parser_write_params_to_file(&param_file, "used_parser_output.yml", 1);
@@ -72,5 +152,30 @@ int main(int argc, char *argv[]) {
   assert(kernel == 4);
   assert(optional == 1);
 
+  assert(csides[0] == 'a');
+  assert(csides[1] == 'b');
+  assert(csides[2] == 'c');
+
+  assert(isides[0] == 2);
+  assert(isides[1] == 3);
+  assert(isides[2] == 4);
+
+  assert(fsides[0] == 2.0);
+  assert(fsides[1] == 3.0);
+  assert(fsides[2] == 4.0);
+
+  assert(dsides[0] == 2.0);
+  assert(dsides[1] == 3.0);
+  assert(dsides[2] == 4.0);
+
+  assert(haveopt == 0);
+  assert(optsides[0] == 1);
+  assert(optsides[1] == 2);
+  assert(optsides[2] == 3);
+  assert(optsides[3] == 4);
+  assert(optsides[4] == 5);
+
+  assert(haveoptwords == 0);
+
   return 0;
 }
diff --git a/tests/testParserInput.yaml b/tests/testParserInput.yaml
index c7fefb3242ab4e140756789aad9979d024f83906..ce1b51399e59d0420836b49d67b95952deda21dc 100644
--- a/tests/testParserInput.yaml
+++ b/tests/testParserInput.yaml
@@ -9,11 +9,27 @@ Scheduler:
 kernel: 4
 
 Simulation:    
-  no_of_time_steps:   10
+  no_of_time_steps:   10    
   max_h:              1.1255
   start_time:         1.23456789
 
 IO:
   #Input file
   ic_file:            ic_file.ini
+
+Words:
+  list: ['xyz', 'ABC', "ab'c", "de:f", "g,hi", "zzz", Hello World, once-again]
+  nonstdlist: 'xyz', 'ABC', "ab'c", "de:f", "g,hi", "zzz", Hello World, once-again
+  long: [a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z]
+
+#  Spelling mistake, should go into unused list.
+Box:
+  csides: [a, b, c]
+  isydes: [2, 3, 4]
+  isides: [2, 3, 4]
+  fsides: [2.0, 3.0, 4.0]
+  dsides: 2.0, 3.0, 4.0
+
+Unused:
+  value: unused
 ...