diff --git a/examples/main.c b/examples/main.c
index 0fce7500af2e3d4fe98feb7e4d3f484658bd2ded..19c1011d5fa7bce54c3bac18c92b6e36b7d3883d 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -58,8 +58,8 @@
 
 int main(int argc, char *argv[]) {
 
-  int c, icount, j, k, N, periodic = 1;
-  long long N_total = -1;
+  int c, icount, j, k, Ngas = 0, Ngpart = 0, periodic = 1;
+  long long N_total[2] = {0};
   int nr_threads = 1, nr_queues = -1;
   int dump_tasks = 0;
   int data[2];
@@ -67,6 +67,7 @@ int main(int argc, char *argv[]) {
   double h_max = -1.0, scaling = 1.0;
   double time_end = DBL_MAX;
   struct part *parts = NULL;
+  struct gpart *gparts = NULL;
   struct space s;
   struct engine e;
   struct UnitSystem us;
@@ -361,7 +362,7 @@ int main(int argc, char *argv[]) {
                  MPI_COMM_WORLD, MPI_INFO_NULL);
 #endif
 #else
-  read_ic_single(ICfileName, dim, &parts, &N, &periodic);
+  read_ic_single(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic);
 #endif
 
   if (myrank == 0) {
@@ -372,24 +373,34 @@ int main(int argc, char *argv[]) {
   }
 
 #if defined(WITH_MPI)
-  long long N_long = N;
-  MPI_Reduce(&N_long, &N_total, 1, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
+  long long N_long[2] = {Ngas, Ngpart};
+  MPI_Reduce(&N_long, &N_total, 2, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
+  N_total[1] -= N_total[0]
 #else
-  N_total = N;
+  N_total[0] = Ngas;
+  N_total[1] = Ngpart - Ngas;
 #endif
-  if (myrank == 0) message("Read %lld particles from the ICs", N_total);
+      if (myrank == 0)
+      message("Read %lld gas particles and %lld DM particles from the ICs",
+              N_total[0], N_total[1]);
 
   /* Apply h scaling */
   if (scaling != 1.0)
-    for (k = 0; k < N; k++) parts[k].h *= scaling;
+    for (k = 0; k < Ngas; k++) parts[k].h *= scaling;
 
   /* Apply shift */
-  if (shift[0] != 0 || shift[1] != 0 || shift[2] != 0)
-    for (k = 0; k < N; k++) {
+  if (shift[0] != 0 || shift[1] != 0 || shift[2] != 0) {
+    for (k = 0; k < Ngas; k++) {
       parts[k].x[0] += shift[0];
       parts[k].x[1] += shift[1];
       parts[k].x[2] += shift[2];
     }
+    for (k = 0; k < Ngpart; k++) {
+      gparts[k].x[0] += shift[0];
+      gparts[k].x[1] += shift[1];
+      gparts[k].x[2] += shift[2];
+    }
+  }
 
   /* Set default number of queues. */
   if (nr_queues < 0) nr_queues = nr_threads;
@@ -399,7 +410,8 @@ int main(int argc, char *argv[]) {
 
   /* Initialize the space with this data. */
   if (myrank == 0) clocks_gettime(&tic);
-  space_init(&s, dim, parts, N, periodic, h_max, myrank == 0);
+  space_init(&s, dim, parts, gparts, Ngas, Ngpart, periodic, h_max,
+             myrank == 0);
   if (myrank == 0 && verbose) {
     clocks_gettime(&toc);
     message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
@@ -482,9 +494,11 @@ int main(int argc, char *argv[]) {
 
   if (myrank == 0) {
     message(
-        "Running on %lld particles until t=%.3e with %i threads and %i "
+        "Running on %lld gas particles and %lld DM particles until t=%.3e with "
+        "%i threads and %i "
         "queues (dt_min=%.3e, dt_max=%.3e)...",
-        N_total, time_end, e.nr_threads, e.sched.nr_queues, e.dt_min, e.dt_max);
+        N_total[0], N_total[1], time_end, e.nr_threads, e.sched.nr_queues,
+        e.dt_min, e.dt_max);
     fflush(stdout);
   }
 
diff --git a/src/cell.c b/src/cell.c
index 3242c6b0b8226719d3f5311534f7945bb6c20a70..df11782048dfa80c697f53feefe8fabc104eb23b 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -345,9 +345,9 @@ void cell_split(struct cell *c) {
 
   int i, j;
   const int count = c->count, gcount = c->gcount;
-  struct part* parts = c->parts;
-  struct xpart* xparts = c->xparts;
-  struct gpart* gparts = c->gparts;
+  struct part *parts = c->parts;
+  struct xpart *xparts = c->xparts;
+  struct gpart *gparts = c->gparts;
   int left[8], right[8];
   double pivot[3];
 
diff --git a/src/common_io.c b/src/common_io.c
index 0e9c5b6cd51ec64e04995cbc47be64b1f36bc131..f3d4b024d20bea91e60d0f730efcd5bf14d0e2a1 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -359,7 +359,9 @@ FILE* prepareXMFfile() {
 
   if (tempFile == NULL) error("Unable to open temporary file.");
 
-  for (int i = 0; fgets(buffer, 1024, tempFile) != NULL && i < counter - 3; i++) {
+  int i = 0;
+  while (fgets(buffer, 1024, tempFile) != NULL && i < counter - 3) {
+    i++;
     fprintf(xmfFile, "%s", buffer);
   }
   fprintf(xmfFile, "\n");
@@ -471,4 +473,85 @@ void writeXMFline(FILE* xmfFile, char* fileName, char* name, long long N,
   fprintf(xmfFile, "</Attribute>\n");
 }
 
+/**
+ * @brief Prepare the DM particles (in gparts) read in for the addition of the
+ * other particle types
+ *
+ * This function assumes that the DM particles are all at the start of the
+ * gparts array
+ *
+ * @param gparts The array of #gpart freshly read in.
+ * @param Ndm The number of DM particles read in.
+ */
+void prepare_dm_gparts(struct gpart* gparts, int Ndm) {
+
+  /* Let's give all these gparts a negative id */
+  for (int i = 0; i < Ndm; ++i) {
+    gparts[i].id = -abs(gparts[i].id);
+  }
+}
+
+/**
+ * @brief Copy every #part into the corresponding #gpart and link them.
+ *
+ * This function assumes that the DM particles are all at the start of the
+ * gparts array and adds the hydro particles afterwards
+ *
+ * @param parts The array of #part freshly read in.
+ * @param gparts The array of #gpart freshly read in with all the DM particles
+ *at the start
+ * @param Ngas The number of gas particles read in.
+ * @param Ndm The number of DM particles read in.
+ */
+void duplicate_hydro_gparts(struct part* parts, struct gpart* gparts, int Ngas,
+                            int Ndm) {
+
+  for (int i = 0; i < Ngas; ++i) {
+
+    /* Duplicate the crucial information */
+    gparts[i + Ndm].x[0] = parts[i].x[0];
+    gparts[i + Ndm].x[1] = parts[i].x[1];
+    gparts[i + Ndm].x[2] = parts[i].x[2];
+
+    gparts[i + Ndm].v_full[0] = parts[i].v[0];
+    gparts[i + Ndm].v_full[1] = parts[i].v[1];
+    gparts[i + Ndm].v_full[2] = parts[i].v[2];
+
+    gparts[i + Ndm].mass = parts[i].mass;
+
+    /* Link the particles */
+    gparts[i + Ndm].part = &parts[i];
+    parts[i].gpart = &gparts[i + Ndm];
+  }
+}
+
+/**
+ * @brief Copy every DM #gpart into the dmparts array.
+ *
+ * @param gparts The array of #gpart containing all particles.
+ * @param Ntot The number of #gpart.
+ * @param dmparts The array of #gpart containg DM particles to be filled.
+ * @param Ndm The number of DM particles.
+ */
+void collect_dm_gparts(struct gpart* gparts, int Ntot, struct gpart* dmparts,
+                       int Ndm) {
+
+  int count = 0;
+
+  /* Loop over all gparts */
+  for (int i = 0; i < Ntot; ++i) {
+
+    /* And collect the DM ones */
+    if (gparts[i].id < 0) {
+      memcpy(&dmparts[count], &gparts[i], sizeof(struct gpart));
+      count++;
+    }
+  }
+
+  /* Check that everything is fine */
+  if (count != Ndm)
+    error("Collected the wrong number of dm particles (%d vs. %d expected)",
+          count, Ndm);
+}
+
 #endif
diff --git a/src/common_io.h b/src/common_io.h
index c221ad3aaf84deb83c0067ee4ece729b98003430..f48fe85a0e651d971209c55469f05e200ee5759d 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -24,6 +24,7 @@
 #include "../config.h"
 
 /* Includes. */
+#include "part.h"
 #include "units.h"
 
 #if defined(HAVE_HDF5)
@@ -55,9 +56,26 @@ enum DATA_IMPORTANCE {
   OPTIONAL = 0
 };
 
+/**
+ * @brief The different particle types present in a GADGET IC file
+ *
+ */
+enum PARTICLE_TYPE {
+  GAS = 0,
+  DM = 1,
+  STAR = 4,
+  BH = 5
+};
+
 hid_t hdf5Type(enum DATA_TYPE type);
 size_t sizeOfType(enum DATA_TYPE type);
 
+void collect_dm_gparts(struct gpart* gparts, int Ntot, struct gpart* dmparts,
+                       int Ndm);
+void prepare_dm_gparts(struct gpart* gparts, int Ndm);
+void duplicate_hydro_gparts(struct part* parts, struct gpart* gparts, int Ngas,
+                            int Ndm);
+
 void readAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data);
 
 void writeAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data,
diff --git a/src/debug.c b/src/debug.c
index 48a0995a5d2992c2b23c7294833309508ee5b59d..4c1434118c98aab7def28d3a53493767d249d774 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -38,6 +38,8 @@
 #error "Invalid choice of SPH variant"
 #endif
 
+#include "./gravity/Default/gravity_debug.h"
+
 /**
  * @brief Looks for the particle with the given id and prints its information to
  *the standard output.
@@ -67,21 +69,20 @@ void printParticle(struct part *parts, struct xpart *xparts, long long int id,
   if (!found) printf("## Particles[???] id=%lld not found\n", id);
 }
 
-void printgParticle(struct gpart *parts, long long int id, size_t N) {
+void printgParticle(struct gpart *gparts, long long int id, size_t N) {
 
   int found = 0;
 
   /* Look for the particle. */
   for (size_t i = 0; i < N; i++)
-    if (parts[i].id == -id || (parts[i].id > 0 && parts[i].part->id == id)) {
-      printf(
-          "## gParticle[%zd]: id=%lld, x=[%.16e,%.16e,%.16e], "
-          "v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e], m=%.3e, t_begin=%d, "
-          "t_end=%d\n",
-          i, parts[i].part->id, parts[i].x[0], parts[i].x[1], parts[i].x[2],
-          parts[i].v[0], parts[i].v[1], parts[i].v[2], parts[i].a[0],
-          parts[i].a[1], parts[i].a[2], parts[i].mass, parts[i].ti_begin,
-          parts[i].ti_end);
+    if (gparts[i].id == -id) {
+      printf("## gParticle[%zd] (DM) :\n id=%lld", i, -gparts[i].id);
+      gravity_debug_particle(&gparts[i]);
+      found = 1;
+      break;
+    } else if (gparts[i].id > 0 && gparts[i].part->id == id) {
+      printf("## gParticle[%zd] (hydro) :\n id=%lld", i, gparts[i].id);
+      gravity_debug_particle(&gparts[i]);
       found = 1;
       break;
     }
diff --git a/src/engine.c b/src/engine.c
index 0739609170eceeeae08f1523330e1d7436688a0b..c06c65fb7a1d2e50016e4370c553769c8c72ac45 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -595,7 +595,8 @@ void engine_exchange_cells(struct engine *e) {
  * @return The number of arrived parts copied to parts and xparts.
  */
 
-int engine_exchange_strays(struct engine *e, int offset, size_t *ind, size_t N) {
+int engine_exchange_strays(struct engine *e, int offset, size_t *ind,
+                           size_t N) {
 
 #ifdef WITH_MPI
 
@@ -1165,8 +1166,8 @@ void engine_print_task_counts(struct engine *e) {
     else
       counts[task_type_count] += 1;
 #ifdef WITH_MPI
-  printf("[%04i] %s engine_print_task_counts: task counts are [ %s=%i", e->nodeID,
-         clocks_get_timesincestart(), taskID_names[0], counts[0]);
+  printf("[%04i] %s engine_print_task_counts: task counts are [ %s=%i",
+         e->nodeID, clocks_get_timesincestart(), taskID_names[0], counts[0]);
 #else
   printf("%s engine_print_task_counts: task counts are [ %s=%i",
          clocks_get_timesincestart(), taskID_names[0], counts[0]);
diff --git a/src/gravity/Default/gravity_debug.h b/src/gravity/Default/gravity_debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..98e0c40a5700b4da70f27fb0955592bb5d2287c3
--- /dev/null
+++ b/src/gravity/Default/gravity_debug.h
@@ -0,0 +1,28 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 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/>.
+ *
+ ******************************************************************************/
+
+__attribute__((always_inline))
+    INLINE static void gravity_debug_particle(struct gpart* p) {
+  printf(
+      "x=[%.3e,%.3e,%.3e], "
+      "v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
+      "mass=%.3e t_begin=%d, t_end=%d\n",
+      p->x[0], p->x[1], p->x[2], p->v_full[0], p->v_full[1], p->v_full[2],
+      p->a[0], p->a[1], p->a[2], p->mass, p->ti_begin, p->ti_end);
+}
diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..a493b581768bf3e9110670435e651c18982b00fb
--- /dev/null
+++ b/src/gravity/Default/gravity_io.h
@@ -0,0 +1,75 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 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/>.
+ *
+ ******************************************************************************/
+
+/**
+ * @brief Reads the different particles to the HDF5 file
+ *
+ * @param h_grp The HDF5 group in which to read the arrays.
+ * @param N The number of particles on that MPI rank.
+ * @param N_total The total number of particles (only used in MPI mode)
+ * @param offset The offset of the particles for this MPI rank (only used in MPI
+ *mode)
+ * @param parts The particle array
+ *
+ */
+__attribute__((always_inline)) INLINE static void darkmatter_read_particles(
+    hid_t h_grp, int N, long long N_total, long long offset,
+    struct gpart* gparts) {
+
+  /* Read arrays */
+  readArray(h_grp, "Coordinates", DOUBLE, N, 3, gparts, N_total, offset, x,
+            COMPULSORY);
+  readArray(h_grp, "Masses", FLOAT, N, 1, gparts, N_total, offset, mass,
+            COMPULSORY);
+  readArray(h_grp, "Velocities", FLOAT, N, 3, gparts, N_total, offset, v_full,
+            COMPULSORY);
+  readArray(h_grp, "ParticleIDs", ULONGLONG, N, 1, gparts, N_total, offset, id,
+            COMPULSORY);
+}
+
+/**
+ * @brief Writes the different particles to the HDF5 file
+ *
+ * @param h_grp The HDF5 group in which to write the arrays.
+ * @param fileName The name of the file (unsued in MPI mode).
+ * @param xmfFile The XMF file to write to (unused in MPI mode).
+ * @param Ndm The number of DM particles on that MPI rank.
+ * @param Ndm_total The total number of g-particles (only used in MPI mode)
+ * @param mpi_rank The MPI rank of this node (only used in MPI mode)
+ * @param offset The offset of the particles for this MPI rank (only used in MPI
+ *mode)
+ * @param gparts The #gpart array
+ * @param us The unit system to use
+ *
+ */
+__attribute__((always_inline)) INLINE static void darkmatter_write_particles(
+    hid_t h_grp, char* fileName, FILE* xmfFile, int Ndm, long long Ndm_total,
+    int mpi_rank, long long offset, struct gpart* gparts,
+    struct UnitSystem* us) {
+
+  /* Write arrays */
+  writeArray(h_grp, fileName, xmfFile, "Coordinates", DOUBLE, Ndm, 3, gparts,
+             Ndm_total, mpi_rank, offset, x, us, UNIT_CONV_LENGTH);
+  writeArray(h_grp, fileName, xmfFile, "Masses", FLOAT, Ndm, 1, gparts,
+             Ndm_total, mpi_rank, offset, mass, us, UNIT_CONV_MASS);
+  writeArray(h_grp, fileName, xmfFile, "Velocities", FLOAT, Ndm, 3, gparts,
+             Ndm_total, mpi_rank, offset, v_full, us, UNIT_CONV_SPEED);
+  writeArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, Ndm, 1, gparts,
+             Ndm_total, mpi_rank, offset, id, us, UNIT_CONV_NO_UNITS);
+}
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index 17262f9d413b1e801b71b456eb4e6baeb925febc..7ce7b81892582f2a90f7dd07f7f244c0d4ed8afb 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -26,7 +26,7 @@ struct gpart {
   double x[3];
 
   /* Particle velocity. */
-  float v[3];
+  float v_full[3];
 
   /* Particle acceleration. */
   float a[3];
@@ -44,7 +44,7 @@ struct gpart {
   union {
 
     /* Particle ID. */
-    size_t id;
+    long long id;
 
     /* Pointer to corresponding SPH part. */
     struct part* part;
diff --git a/src/gravity_io.h b/src/gravity_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..6276e50473de64abd1014bdb36a63a14e02ca8cf
--- /dev/null
+++ b/src/gravity_io.h
@@ -0,0 +1,26 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 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_GRAVITY_IO_H
+#define SWIFT_GRAVITY_IO_H
+
+#include "./const.h"
+
+#include "./gravity/Default/gravity_io.h"
+
+#endif /* SWIFT_GRAVITY_IO_H */
diff --git a/src/hydro/Default/hydro_part.h b/src/hydro/Default/hydro_part.h
index 4dbc0388ac20cd7f354007e2e80aa6fccccdc5ad..60453d0c7995f7af2a3166502a24aa590873a043 100644
--- a/src/hydro/Default/hydro_part.h
+++ b/src/hydro/Default/hydro_part.h
@@ -111,7 +111,7 @@ struct part {
   float mass;
 
   /* Particle ID. */
-  unsigned long long id;
+  long long id;
 
   /* Pointer to corresponding gravity part. */
   struct gpart* gpart;
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index 298cfc5ac1d801f03f7ba104a819fd3de531d55c..05754d07dd70bed071e99c86b95eb17eb2194012 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/src/hydro/Gadget2/hydro_part.h
@@ -104,7 +104,7 @@ struct part {
   float div_v;
 
   /* Particle ID. */
-  unsigned long long id;
+  long long id;
 
   /* Pointer to corresponding gravity part. */
   struct gpart* gpart;
diff --git a/src/hydro/Minimal/hydro_part.h b/src/hydro/Minimal/hydro_part.h
index 8875e737c6f3c0ab0b7a390b56e9b946c145fc1b..173397ef2c72ee99f4d10742f3645afd1e706218 100644
--- a/src/hydro/Minimal/hydro_part.h
+++ b/src/hydro/Minimal/hydro_part.h
@@ -101,7 +101,7 @@ struct part {
     } force;
   };
 
-  unsigned long long id; /*!< Particle unique ID. */
+  long long id; /*!< Particle unique ID. */
 
   struct gpart* gpart; /*!< Pointer to corresponding gravity part. */
 
diff --git a/src/part.h b/src/part.h
index 168e80b68bb211d1d773f1f80f1fa9cc757edfcb..865403e8c2c157dc5a8ff7a32bc41be676d7919b 100644
--- a/src/part.h
+++ b/src/part.h
@@ -35,6 +35,7 @@
 
 /* Some constants. */
 #define part_align 64
+#define gpart_align 32
 #define xpart_align 32
 
 /* Import the right particle definition */
diff --git a/src/single_io.c b/src/single_io.c
index f3225ecf487d223c3bea560cbae626a03da36928..270ba429a5ef54c9501fab9c12159bd4a83c3a77 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -199,9 +199,6 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
     chunk_shape[1] = 0;
   }
 
-  /* Make sure the chunks are not larger than the dataset */
-  if (chunk_shape[0] > N) chunk_shape[0] = N;
-
   /* Change shape of data space */
   h_err = H5Sset_extent_simple(h_space, rank, shape, NULL);
   if (h_err < 0) {
@@ -302,6 +299,8 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
 
 /* Import the right hydro definition */
 #include "hydro_io.h"
+/* Import the right gravity definition */
+#include "gravity_io.h"
 
 /**
  * @brief Reads an HDF5 initial condition file (GADGET-3 type)
@@ -322,13 +321,15 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
  * Calls #error() if an error occurs.
  *
  */
-void read_ic_single(char* fileName, double dim[3], struct part** parts, int* N,
+void read_ic_single(char* fileName, double dim[3], struct part** parts,
+                    struct gpart** gparts, int* Ngas, int* Ngparts,
                     int* periodic) {
   hid_t h_file = 0, h_grp = 0;
-  double boxSize[3] = {0.0, -1.0, -1.0};
   /* GADGET has only cubic boxes (in cosmological mode) */
+  double boxSize[3] = {0.0, -1.0, -1.0};
+  /* GADGET has 6 particle types. We only keep the type 0 & 1 for now...*/
   int numParticles[6] = {0};
-  /* GADGET has 6 particle types. We only keep the type 0*/
+  int Ndm;
 
   /* Open file */
   /* message("Opening file '%s' as IC.", fileName); */
@@ -357,35 +358,81 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, int* N,
   readAttribute(h_grp, "BoxSize", DOUBLE, boxSize);
   readAttribute(h_grp, "NumPart_Total", UINT, numParticles);
 
-  *N = numParticles[0];
+  *Ngas = numParticles[0];
+  Ndm = numParticles[1];
   dim[0] = boxSize[0];
   dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1];
   dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2];
 
   /* message("Found %d particles in a %speriodic box of size [%f %f %f].",  */
-  /* 	 *N, (periodic ? "": "non-"), dim[0], dim[1], dim[2]); */
+  /* 	  *N, (periodic ? "": "non-"), dim[0], dim[1], dim[2]);  */
 
   /* Close header */
   H5Gclose(h_grp);
 
-  /* Allocate memory to store particles */
-  if (posix_memalign((void*)parts, part_align, *N * sizeof(struct part)) != 0)
-    error("Error while allocating memory for particles");
-  bzero(*parts, *N * sizeof(struct part));
+  /* Total number of particles */
+  *Ngparts = *Ngas + Ndm;
+
+  /* Allocate memory to store SPH particles */
+  if (posix_memalign((void*)parts, part_align, *Ngas * sizeof(struct part)) !=
+      0)
+    error("Error while allocating memory for SPH particles");
+  bzero(*parts, *Ngas * sizeof(struct part));
+
+  /* Allocate memory to store all particles */
+  if (posix_memalign((void*)gparts, gpart_align,
+                     *Ngparts * sizeof(struct gpart)) != 0)
+    error("Error while allocating memory for gravity particles");
+  bzero(*gparts, *Ngparts * sizeof(struct gpart));
 
   /* message("Allocated %8.2f MB for particles.", *N * sizeof(struct part) /
    * (1024.*1024.)); */
 
   /* Open SPH particles group */
   /* message("Reading particle arrays..."); */
-  h_grp = H5Gopen(h_file, "/PartType0", H5P_DEFAULT);
-  if (h_grp < 0) error("Error while opening particle group.\n");
+  message("BoxSize = %lf\n", dim[0]);
+  message("NumPart = [%d, %d] Total = %d\n", *Ngas, Ndm, *Ngparts);
 
-  /* Read particle fields into the particle structure */
-  hydro_read_particles(h_grp, *N, *N, 0, *parts);
+  /* Loop over all particle types */
+  for (int ptype = 0; ptype < 6; ptype++) {
 
-  /* Close particle group */
-  H5Gclose(h_grp);
+    /* Don't do anything if no particle of this kind */
+    if (numParticles[ptype] == 0) continue;
+
+    /* Open the particle group in the file */
+    char partTypeGroupName[15];
+    sprintf(partTypeGroupName, "/PartType%d", ptype);
+    h_grp = H5Gopen(h_file, partTypeGroupName, H5P_DEFAULT);
+    if (h_grp < 0) {
+      error("Error while opening particle group %s.", partTypeGroupName);
+    }
+
+    /* message("Group %s found - reading...", partTypeGroupName); */
+
+    /* Read particle fields into the particle structure */
+    switch (ptype) {
+
+      case GAS:
+        hydro_read_particles(h_grp, *Ngas, *Ngas, 0, *parts);
+        break;
+
+      case DM:
+        darkmatter_read_particles(h_grp, Ndm, Ndm, 0, *gparts);
+        break;
+
+      default:
+        error("Particle Type %d not yet supported. Aborting", ptype);
+    }
+
+    /* Close particle group */
+    H5Gclose(h_grp);
+  }
+
+  /* Prepare the DM particles */
+  prepare_dm_gparts(*gparts, Ndm);
+
+  /* Now duplicate the hydro particle into gparts */
+  duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
 
   /* message("Done Reading particles..."); */
 
@@ -410,15 +457,21 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, int* N,
 void write_output_single(struct engine* e, struct UnitSystem* us) {
 
   hid_t h_file = 0, h_grp = 0, h_grpsph = 0;
-  int N = e->s->nr_parts;
+  const int Ngas = e->s->nr_parts;
+  const int Ntot = e->s->nr_gparts;
   int periodic = e->s->periodic;
-  int numParticles[6] = {N, 0};
-  int numParticlesHighWord[6] = {0};
   int numFiles = 1;
   struct part* parts = e->s->parts;
+  struct gpart* gparts = e->s->gparts;
+  struct gpart* dmparts = NULL;
   FILE* xmfFile = 0;
   static int outputCount = 0;
 
+  /* Number of particles of each type */
+  const int Ndm = Ntot - Ngas;
+  int numParticles[6] = {Ngas, Ndm, 0};
+  int numParticlesHighWord[6] = {0};
+
   /* File name */
   char fileName[200];
   sprintf(fileName, "output_%03i.hdf5", outputCount);
@@ -430,7 +483,7 @@ void write_output_single(struct engine* e, struct UnitSystem* us) {
   xmfFile = prepareXMFfile();
 
   /* Write the part corresponding to this specific output */
-  writeXMFheader(xmfFile, N, fileName, e->time);
+  writeXMFheader(xmfFile, Ngas, fileName, e->time);
 
   /* Open file */
   /* message("Opening file '%s'.", fileName); */
@@ -486,17 +539,56 @@ void write_output_single(struct engine* e, struct UnitSystem* us) {
   /* Print the system of Units */
   writeUnitSystem(h_file, us);
 
-  /* Create SPH particles group */
-  /* message("Writing particle arrays..."); */
-  h_grp =
-      H5Gcreate(h_file, "/PartType0", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-  if (h_grp < 0) error("Error while creating particle group.\n");
+  /* Loop over all particle types */
+  for (int ptype = 0; ptype < 6; ptype++) {
 
-  /* Write particle fields from the particle structure */
-  hydro_write_particles(h_grp, fileName, xmfFile, N, N, 0, 0, parts, us);
+    /* Don't do anything if no particle of this kind */
+    if (numParticles[ptype] == 0) continue;
 
-  /* Close particle group */
-  H5Gclose(h_grp);
+    /* Open the particle group in the file */
+    char partTypeGroupName[15];
+    sprintf(partTypeGroupName, "/PartType%d", ptype);
+    h_grp = H5Gcreate(h_file, partTypeGroupName, H5P_DEFAULT, H5P_DEFAULT,
+                      H5P_DEFAULT);
+    if (h_grp < 0) {
+      error("Error while creating particle group.\n");
+    }
+
+    /* message("Writing particle arrays..."); */
+
+    /* Write particle fields from the particle structure */
+    switch (ptype) {
+
+      case GAS:
+        hydro_write_particles(h_grp, fileName, xmfFile, Ngas, Ngas, 0, 0, parts,
+                              us);
+        break;
+
+      case DM:
+        /* Allocate temporary array */
+        if (posix_memalign((void*)dmparts, gpart_align,
+                           Ndm * sizeof(struct gpart)) != 0)
+          error("Error while allocating temporart memory for DM particles");
+        bzero(&dmparts, Ndm * sizeof(struct gpart));
+
+        /* Collect the DM particles from gpart */
+        collect_dm_gparts(gparts, Ntot, dmparts, Ndm);
+
+        /* Write DM particles */
+        darkmatter_write_particles(h_grp, fileName, xmfFile, Ndm, Ndm, 0, 0,
+                                   dmparts, us);
+
+        /* Free temporary array */
+        free(dmparts);
+        break;
+
+      default:
+        error("Particle Type %d not yet supported. Aborting", ptype);
+    }
+
+    /* Close particle group */
+    H5Gclose(h_grp);
+  }
 
   /* Write LXMF file descriptor */
   writeXMFfooter(xmfFile);
diff --git a/src/single_io.h b/src/single_io.h
index f6689901106a2cd5d85a873d4047e1d21edd3547..4f4aa10536a890a9a46ba752a17463557f27d5e7 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -26,8 +26,8 @@
 
 #if defined(HAVE_HDF5) && !defined(WITH_MPI)
 
-void read_ic_single(char* fileName, double dim[3], struct part** parts, int* N,
-                    int* periodic);
+void read_ic_single(char* fileName, double dim[3], struct part** parts,
+                    struct gpart** gparts, int* Ngas, int* Ndm, int* periodic);
 
 void write_output_single(struct engine* e, struct UnitSystem* us);
 
diff --git a/src/space.c b/src/space.c
index 5938101d94b17097277d7fc06a1cd4d3c1181727..4d1ebfca87a933a8516ab144fc06624f46c0cc6c 100644
--- a/src/space.c
+++ b/src/space.c
@@ -1188,8 +1188,9 @@ struct cell *space_getcell(struct space *s) {
  * recursively.
  */
 
-void space_init(struct space *s, double dim[3], struct part *parts, size_t N,
-                int periodic, double h_max, int verbose) {
+void space_init(struct space *s, double dim[3], struct part *parts,
+                struct gpart *gparts, size_t N, size_t Ngpart, int periodic,
+                double h_max, int verbose) {
 
   /* Store everything in the space. */
   s->dim[0] = dim[0];
diff --git a/src/space.h b/src/space.h
index 87b618b141b7871aeb5ff76eae5bdc24f063bdef..91485ff7e2ebe9da8ab927748589ae9f71320803 100644
--- a/src/space.h
+++ b/src/space.h
@@ -127,12 +127,14 @@ extern struct parallel_sort space_sort_struct;
 /* function prototypes. */
 void space_parts_sort(struct space *s, size_t *ind, size_t N, int min, int max,
                       int verbose);
-void space_gparts_sort(struct gpart *gparts, size_t *ind, size_t N, int min, int max);
+void space_gparts_sort(struct gpart *gparts, size_t *ind, size_t N, int min,
+                       int max);
 struct cell *space_getcell(struct space *s);
 int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
                  double *shift);
-void space_init(struct space *s, double dim[3], struct part *parts, size_t N,
-                int periodic, double h_max, int verbose);
+void space_init(struct space *s, double dim[3], struct part *parts,
+                struct gpart *gparts, size_t N, size_t Ngpart, int periodic,
+                double h_max, int verbose);
 void space_map_cells_pre(struct space *s, int full,
                          void (*fun)(struct cell *c, void *data), void *data);
 void space_map_parts(struct space *s,