diff --git a/src/parallel_io.c b/src/parallel_io.c
index cb02d1806c6ede5edfeaeb13d61cbecb09ac57e6..600eee21cece8692f839b94a53f8de819b450ab4 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -1171,22 +1171,28 @@ void write_output_parallel(struct engine* e, const char* baseName,
                            int mpi_rank, int mpi_size, MPI_Comm comm,
                            MPI_Info info) {
 
-  const size_t Ngas = e->s->nr_parts;
-  const size_t Nstars = e->s->nr_sparts;
-  const size_t Ntot = e->s->nr_gparts;
   const struct part* parts = e->s->parts;
   const struct xpart* xparts = e->s->xparts;
   const struct gpart* gparts = e->s->gparts;
-  struct gpart* dmparts = NULL;
   const struct spart* sparts = e->s->sparts;
-  const struct cooling_function_data* cooling = e->cooling_func;
   struct swift_params* params = e->parameter_file;
 
-  /* Number of unassociated gparts */
-  const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
+  /* Number of particles currently in the arrays */
+  const size_t Ntot = e->s->nr_gparts;
+  const size_t Ngas = e->s->nr_parts;
+  const size_t Nstars = e->s->nr_sparts;
+  //const size_t Nbaryons = Ngas + Nstars;
+  //const size_t Ndm = Ntot > 0 ? Ntot - Nbaryons : 0;
+
+  /* Number of particles that we will write */
+  const size_t Ntot_written = e->s->nr_gparts - e->s->nr_inhibited_sparts;  
+  const size_t Ngas_written = e->s->nr_parts - e->s->nr_inhibited_parts;
+  const size_t Nstars_written = e->s->nr_sparts - e->s->nr_inhibited_gparts;
+  const size_t Nbaryons_written = Ngas_written + Nstars_written;
+  const size_t Ndm_written = Ntot_written > 0 ? Ntot_written - Nbaryons_written : 0;
 
   /* Compute offset in the file and total number of particles */
-  size_t N[swift_type_count] = {Ngas, Ndm, 0, 0, Nstars, 0};
+  size_t N[swift_type_count] = {Ngas_written, Ndm_written, 0, 0, Nstars_written, 0};
   long long N_total[swift_type_count] = {0};
   long long offset[swift_type_count] = {0};
   MPI_Exscan(&N, &offset, swift_type_count, MPI_LONG_LONG_INT, MPI_SUM, comm);
@@ -1337,37 +1343,96 @@ void write_output_parallel(struct engine* e, const char* baseName,
     struct io_props list[100];
     size_t Nparticles = 0;
 
+    struct part* parts_written = NULL;
+    struct xpart* xparts_written = NULL;
+    struct gpart* gparts_written = NULL;
+    struct spart* sparts_written = NULL;
+    
     /* Write particle fields from the particle structure */
     switch (ptype) {
 
       case swift_type_gas:
-        Nparticles = Ngas;
-        hydro_write_particles(parts, xparts, list, &num_fields);
-        num_fields += chemistry_write_particles(parts, list + num_fields);
-        num_fields +=
-            cooling_write_particles(xparts, list + num_fields, cooling);
+	{
+	  if (Ngas == Ngas_written) {
+
+	    /* No inhibted particles: easy case */
+	    Nparticles = Ngas;
+	    hydro_write_particles(parts, xparts, list, &num_fields);
+	    num_fields += chemistry_write_particles(parts, list + num_fields);
+	    num_fields += cooling_write_particles(xparts, list + num_fields, e->cooling_func);
+	  }
+	  else {
+
+	    /* Ok, we need to fish out the particles we want */
+	    Nparticles = Ngas_written;
+
+	    /* Allocate temporary arrays */
+	    if (posix_memalign((void**)&parts_written, part_align, Ngas_written * sizeof(struct part)) != 0)
+	      error("Error while allocating temporart memory for parts");
+	    if (posix_memalign((void**)&xparts_written, xpart_align, Ngas_written * sizeof(struct xpart)) != 0)
+	      error("Error while allocating temporart memory for xparts");
+
+	    /* Collect the particles we want to write */
+	    io_collect_parts_to_write(parts, xparts, parts_written, xparts_written, Ngas, Ngas_written);
+
+	    /* Select the fields to write */
+	    hydro_write_particles(parts_written, xparts_written, list, &num_fields);
+	    num_fields += chemistry_write_particles(parts_written, list + num_fields);
+	    num_fields += cooling_write_particles(xparts_written, list + num_fields, e->cooling_func);
+	  }	  
+	}
         break;
 
       case swift_type_dark_matter:
-        /* 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 */
-        io_collect_dm_gparts(gparts, Ntot, dmparts, Ndm);
-
-        /* Write DM particles */
-        Nparticles = Ndm;
-        darkmatter_write_particles(dmparts, list, &num_fields);
+	{
+	  if(Ntot == Ndm_written) {
+
+	    /* This is a DM-only run without inhibited particles */
+	    Nparticles = Ntot;
+	    darkmatter_write_particles(gparts, list, &num_fields);
+	  }
+	  else {
+	    
+	    /* Ok, we need to fish out the particles we want */
+	    Nparticles = Ndm_written;
+	    
+	    /* Allocate temporary array */
+	    if (posix_memalign((void**)&gparts_written, gpart_align, Ndm_written * sizeof(struct gpart)) != 0)
+	      error("Error while allocating temporart memory for gparts");
+	    
+	    /* Collect the non-inhibited DM particles from gpart */
+	    io_collect_gparts_to_write(gparts, gparts_written, Ntot, Ndm_written);
+	    
+	    /* Write DM particles */
+	    darkmatter_write_particles(gparts_written, list, &num_fields);
+	  }
+	}
         break;
 
       case swift_type_stars:
-        Nparticles = Nstars;
-        stars_write_particles(sparts, list, &num_fields);
+	{
+	  if(Nstars == Nstars_written) {
+
+	    /* No inhibted particles: easy case */
+	    Nparticles = Nstars;
+	    stars_write_particles(sparts, list, &num_fields);
+	  }
+	  else {
+
+	    /* Ok, we need to fish out the particles we want */
+	    Nparticles = Nstars_written;
+
+	    /* Allocate temporary arrays */
+	    if (posix_memalign((void**)&sparts_written, spart_align, Nstars_written * sizeof(struct spart)) != 0)
+	      error("Error while allocating temporart memory for sparts");
+
+	    /* Collect the particles we want to write */
+	    io_collect_sparts_to_write(sparts, sparts_written, Nstars, Nstars_written);
+
+	    /* Select the fields to write */
+	    stars_write_particles(sparts, list, &num_fields);
+	  }
+	}
         break;
 
       default:
@@ -1390,10 +1455,10 @@ void write_output_parallel(struct engine* e, const char* baseName,
     }
 
     /* Free temporary array */
-    if (dmparts) {
-      free(dmparts);
-      dmparts = 0;
-    }
+    if(parts_written) free(parts_written);
+    if(xparts_written) free(xparts_written);
+    if(gparts_written) free(gparts_written);
+    if(sparts_written) free(sparts_written);
 
 #ifdef IO_SPEED_MEASUREMENT
     MPI_Barrier(MPI_COMM_WORLD);
diff --git a/src/serial_io.c b/src/serial_io.c
index 32329a296f98687f4a44ad806b52eb259e65207c..5ae2062482a2ed0be62de006031e4ec135d329b3 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -739,22 +739,28 @@ void write_output_serial(struct engine* e, const char* baseName,
                          int mpi_size, MPI_Comm comm, MPI_Info info) {
 
   hid_t h_file = 0, h_grp = 0;
-  const size_t Ngas = e->s->nr_parts;
-  const size_t Nstars = e->s->nr_sparts;
-  const size_t Ntot = e->s->nr_gparts;
   int periodic = e->s->periodic;
   int numFiles = 1;
   const struct part* parts = e->s->parts;
   const struct xpart* xparts = e->s->xparts;
   const struct gpart* gparts = e->s->gparts;
-  struct gpart* dmparts = NULL;
   const struct spart* sparts = e->s->sparts;
-  const struct cooling_function_data* cooling = e->cooling_func;
   struct swift_params* params = e->parameter_file;
   FILE* xmfFile = 0;
 
-  /* Number of unassociated gparts */
-  const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
+  /* Number of particles currently in the arrays */
+  const size_t Ntot = e->s->nr_gparts;
+  const size_t Ngas = e->s->nr_parts;
+  const size_t Nstars = e->s->nr_sparts;
+  //const size_t Nbaryons = Ngas + Nstars;
+  //const size_t Ndm = Ntot > 0 ? Ntot - Nbaryons : 0;
+
+  /* Number of particles that we will write */
+  const size_t Ntot_written = e->s->nr_gparts - e->s->nr_inhibited_sparts;  
+  const size_t Ngas_written = e->s->nr_parts - e->s->nr_inhibited_parts;
+  const size_t Nstars_written = e->s->nr_sparts - e->s->nr_inhibited_gparts;
+  const size_t Nbaryons_written = Ngas_written + Nstars_written;
+  const size_t Ndm_written = Ntot_written > 0 ? Ntot_written - Nbaryons_written : 0;
 
   /* File name */
   char fileName[FILENAME_BUFFER_SIZE];
@@ -766,7 +772,7 @@ void write_output_serial(struct engine* e, const char* baseName,
              e->snapshot_output_count);
 
   /* Compute offset in the file and total number of particles */
-  size_t N[swift_type_count] = {Ngas, Ndm, 0, 0, Nstars, 0};
+  size_t N[swift_type_count] = {Ngas_written, Ndm_written, 0, 0, Nstars_written, 0};
   long long N_total[swift_type_count] = {0};
   long long offset[swift_type_count] = {0};
   MPI_Exscan(&N, &offset, swift_type_count, MPI_LONG_LONG_INT, MPI_SUM, comm);
@@ -1012,35 +1018,96 @@ void write_output_serial(struct engine* e, const char* baseName,
         struct io_props list[100];
         size_t Nparticles = 0;
 
+	struct part* parts_written = NULL;
+	struct xpart* xparts_written = NULL;
+	struct gpart* gparts_written = NULL;
+	struct spart* sparts_written = NULL;
+
         /* Write particle fields from the particle structure */
         switch (ptype) {
 
           case swift_type_gas:
-            Nparticles = Ngas;
-            hydro_write_particles(parts, xparts, list, &num_fields);
-            num_fields += chemistry_write_particles(parts, list + num_fields);
-            num_fields +=
-                cooling_write_particles(xparts, list + num_fields, cooling);
-            break;
+	    {
+	      if (Ngas == Ngas_written) {
+		
+		/* No inhibted particles: easy case */
+		Nparticles = Ngas;
+		hydro_write_particles(parts, xparts, list, &num_fields);
+		num_fields += chemistry_write_particles(parts, list + num_fields);
+		num_fields += cooling_write_particles(xparts, list + num_fields, e->cooling_func);
+	      }
+	      else {
+		
+		/* Ok, we need to fish out the particles we want */
+		Nparticles = Ngas_written;
+		
+		/* Allocate temporary arrays */
+		if (posix_memalign((void**)&parts_written, part_align, Ngas_written * sizeof(struct part)) != 0)
+		  error("Error while allocating temporart memory for parts");
+		if (posix_memalign((void**)&xparts_written, xpart_align, Ngas_written * sizeof(struct xpart)) != 0)
+		  error("Error while allocating temporart memory for xparts");
+		
+		/* Collect the particles we want to write */
+		io_collect_parts_to_write(parts, xparts, parts_written, xparts_written, Ngas, Ngas_written);
+		
+		/* Select the fields to write */
+		hydro_write_particles(parts_written, xparts_written, list, &num_fields);
+		num_fields += chemistry_write_particles(parts_written, list + num_fields);
+		num_fields += cooling_write_particles(xparts_written, list + num_fields, e->cooling_func);
+	      }	  
+	    }
+	    break;
 
           case swift_type_dark_matter:
-            /* 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 */
-            io_collect_dm_gparts(gparts, Ntot, dmparts, Ndm);
-
-            /* Write DM particles */
-            Nparticles = Ndm;
-            darkmatter_write_particles(dmparts, list, &num_fields);
+	    {
+	      if(Ntot == Ndm_written) {
+		
+		/* This is a DM-only run without inhibited particles */
+		Nparticles = Ntot;
+		darkmatter_write_particles(gparts, list, &num_fields);
+	      }
+	      else {
+		
+		/* Ok, we need to fish out the particles we want */
+		Nparticles = Ndm_written;
+		
+		/* Allocate temporary array */
+		if (posix_memalign((void**)&gparts_written, gpart_align, Ndm_written * sizeof(struct gpart)) != 0)
+		  error("Error while allocating temporart memory for gparts");
+		
+		/* Collect the non-inhibited DM particles from gpart */
+		io_collect_gparts_to_write(gparts, gparts_written, Ntot, Ndm_written);
+		
+		/* Write DM particles */
+		darkmatter_write_particles(gparts_written, list, &num_fields);
+	      }
+	    }
             break;
-
+	    
           case swift_type_stars:
-            Nparticles = Nstars;
-            stars_write_particles(sparts, list, &num_fields);
+	    {
+	      if(Nstars == Nstars_written) {
+		
+		/* No inhibted particles: easy case */
+		Nparticles = Nstars;
+		stars_write_particles(sparts, list, &num_fields);
+	      }
+	      else {
+		
+		/* Ok, we need to fish out the particles we want */
+		Nparticles = Nstars_written;
+		
+		/* Allocate temporary arrays */
+		if (posix_memalign((void**)&sparts_written, spart_align, Nstars_written * sizeof(struct spart)) != 0)
+		  error("Error while allocating temporart memory for sparts");
+		
+		/* Collect the particles we want to write */
+		io_collect_sparts_to_write(sparts, sparts_written, Nstars, Nstars_written);
+		
+		/* Select the fields to write */
+		stars_write_particles(sparts, list, &num_fields);
+	      }
+	    }
             break;
 
           default:
@@ -1063,10 +1130,10 @@ void write_output_serial(struct engine* e, const char* baseName,
         }
 
         /* Free temporary array */
-        if (dmparts) {
-          free(dmparts);
-          dmparts = 0;
-        }
+	if(parts_written) free(parts_written);
+	if(xparts_written) free(xparts_written);
+	if(gparts_written) free(gparts_written);
+	if(sparts_written) free(sparts_written);
 
         /* Close particle group */
         H5Gclose(h_grp);
diff --git a/src/single_io.c b/src/single_io.c
index 4a8e9ed25e7ede2432baffeeb38d38cdb52471cc..2586b2909a7c6df8d21c2e9ddbcb3c5dd29e6833 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -609,7 +609,6 @@ void write_output_single(struct engine* e, const char* baseName,
   const struct xpart* xparts = e->s->xparts;
   const struct gpart* gparts = e->s->gparts;
   const struct spart* sparts = e->s->sparts;
-  const struct cooling_function_data* cooling = e->cooling_func;
   struct swift_params* params = e->parameter_file;
   
   /* Number of particles currently in the arrays */
@@ -852,7 +851,7 @@ void write_output_single(struct engine* e, const char* baseName,
 	    N = Ngas;
 	    hydro_write_particles(parts, xparts, list, &num_fields);
 	    num_fields += chemistry_write_particles(parts, list + num_fields);
-	    num_fields += cooling_write_particles(xparts, list + num_fields, cooling);
+	    num_fields += cooling_write_particles(xparts, list + num_fields, e->cooling_func);
 	  }
 	  else {
 
@@ -871,7 +870,7 @@ void write_output_single(struct engine* e, const char* baseName,
 	    /* Select the fields to write */
 	    hydro_write_particles(parts_written, xparts_written, list, &num_fields);
 	    num_fields += chemistry_write_particles(parts_written, list + num_fields);
-	    num_fields += cooling_write_particles(xparts_written, list + num_fields, cooling);
+	    num_fields += cooling_write_particles(xparts_written, list + num_fields, e->cooling_func);
 	  }	  
 	}
         break;