diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index afe7b52f1e024363884b2347eaa0ac77653cd50f..49c1b77ab770731840c8d127684d6bd2f9732e25 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -784,6 +784,7 @@ INPUT		       += @top_srcdir@/src/stars/GEAR
 INPUT		       += @top_srcdir@/src/feedback/EAGLE
 INPUT		       += @top_srcdir@/src/feedback/GEAR
 INPUT		       += @top_srcdir@/src/black_holes/EAGLE
+INPUT		       += @top_srcdir@/src/sink/Default
 INPUT		       += @top_srcdir@/logger
 
 # This tag can be used to specify the character encoding of the source files
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_sink/getIC.sh b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/getIC.sh
new file mode 100755
index 0000000000000000000000000000000000000000..a3d16db27aac06abda683a7bd75e72a275f8b9d4
--- /dev/null
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/getIC.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/3e11-star-only-DM-halo-galaxy.hdf5
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_sink/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/isolated_galaxy.yml
new file mode 100644
index 0000000000000000000000000000000000000000..c83abbd55d207095655cca733fd77f596a5b6c15
--- /dev/null
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/isolated_galaxy.yml
@@ -0,0 +1,42 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.9891E43   # 10^10 solar masses 
+  UnitLength_in_cgs:   3.08567758E21   # 1 kpc 
+  UnitVelocity_in_cgs: 1E5   # km/s
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:          0.025               # Constant dimensionless multiplier for time integration.
+  theta:        0.7                 # Opening angle (Multipole acceptance criterion).
+  max_physical_DM_softening: 0.7    # Physical softening length (in internal units).
+  max_physical_baryon_softening: 0.2  # Physical softening length (in internal units)
+
+# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
+TimeIntegration:
+  time_begin:        0.    # The starting time of the simulation (in internal units).
+  time_end:          1.    # The end time of the simulation (in internal units).
+  dt_min:            1e-6  # The minimal time-step size of the simulation (in internal units).
+  dt_max:            1e-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:   output      # Common part of the name of output files
+  time_first: 0.          # (Optional) Time of the first output if non-cosmological time-integration (in internal units)
+  delta_time: 0.01        # Time difference between consecutive outputs (in internal units)
+
+Scheduler:
+  max_top_level_cells:  96
+  
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:           1e-1     # Time between statistics output
+  time_first:             0.     # (Optional) Time of the first stats output if non-cosmological time-integration (in internal units)
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  isolated_galaxy.hdf5  # The file to read
+  periodic:                    0    # Are we running with periodic ICs?
+
+ 
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_sink/makeIC.py b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..3587c9d858ab1a72a39b1b243b6d58bae7fad383
--- /dev/null
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/makeIC.py
@@ -0,0 +1,72 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (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/>.
+#
+##############################################################################
+
+import h5py
+import numpy as np
+from shutil import copyfile
+
+# Add sink particles to the isolated galaxy
+
+fileName = "3e11-star-only-DM-halo-galaxy.hdf5"
+output = "isolated_galaxy.hdf5"
+
+L = 13756  # Size of the box (internal units)
+L_sink = 1000 # Size of the sink particle area (L_sink < L)
+N = 1000  # Number of sink particles
+max_vel = 100  # Maximal velocity of the sink particles (in km / s)
+mass = 0.000142  # Mass of a sink particle (internal units)
+min_id = 360001  # minimal ids of the other particles
+
+# ---------------------------------------------------
+
+copyfile(fileName, output)
+
+# File
+file = h5py.File(output, 'a')
+
+pos = np.random.rand(N, 3) * L_sink + 0.5 * (L - L_sink)
+vel = 2 * (np .random.rand(N, 3) - 0.5) * max_vel
+m = mass * np.ones(N)
+# Generate extra arrays
+ids = np.linspace(min_id, min_id + N, N)
+
+# --------------------------------------------------
+
+# Header
+grp = file["Header"]
+numpart = grp.attrs["NumPart_Total"][:]
+numpart[3] = N
+grp.attrs["NumPart_Total"] = numpart
+grp.attrs["NumPart_ThisFile"] = numpart
+
+# Units
+grp = file["Units"]
+uv = grp.attrs["Unit length in cgs (U_L)"]
+uv /= grp.attrs["Unit time in cgs (U_t)"]
+
+vel *= 1e5 / uv
+
+# Particle group
+grp = file.create_group("/PartType3")
+grp.create_dataset('Coordinates', data=pos, dtype='d')
+grp.create_dataset('Velocities', data=vel, dtype='f')
+grp.create_dataset('Masses', data=m, dtype='f')
+grp.create_dataset('ParticleIDs', data=ids, dtype='L')
+
+file.close()
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_sink/run.sh b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..7d2cee8c29c0eadc886b3188d532567ea861bde5
--- /dev/null
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/run.sh
@@ -0,0 +1,12 @@
+#!/bin/bash
+
+if [ ! -e 3e11-star-only-DM-halo-galaxy.hdf5 ]
+then
+    echo "Fetching initial conditons for the isolated galaxy with an external potential ..."
+    ./getIC.sh
+fi
+
+echo "Generate initial conditions"
+python3 makeIC.py
+
+../../swift --sinks --external-gravity --self-gravity --threads=16 isolated_galaxy.yml 2>&1 | tee output.log
diff --git a/examples/main.c b/examples/main.c
index a24714e9090f66c25bd55b8b48c670279989139e..60de1f4015582fed82b12c2a62629741e0e6e1b3 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -105,6 +105,7 @@ int main(int argc, char *argv[]) {
   struct space s;
   struct spart *sparts = NULL;
   struct bpart *bparts = NULL;
+  struct sink *sinks = NULL;
   struct unit_system us;
   struct los_props los_properties;
 
@@ -171,6 +172,7 @@ int main(int argc, char *argv[]) {
   int with_mpole_reconstruction = 0;
   int with_structure_finding = 0;
   int with_logger = 0;
+  int with_sink = 0;
   int with_qla = 0;
   int with_eagle = 0;
   int with_gear = 0;
@@ -221,6 +223,8 @@ int main(int argc, char *argv[]) {
       OPT_BOOLEAN('S', "stars", &with_stars, "Run with stars.", NULL, 0, 0),
       OPT_BOOLEAN('B', "black-holes", &with_black_holes,
                   "Run with black holes.", NULL, 0, 0),
+      OPT_BOOLEAN('k', "sinks", &with_sink, "Run with sink particles.", NULL, 0,
+                  0),
       OPT_BOOLEAN(
           'u', "fof", &with_fof,
           "Run Friends-of-Friends algorithm to perform black hole seeding.",
@@ -439,6 +443,13 @@ int main(int argc, char *argv[]) {
   }
 #endif
 
+#ifdef WITH_MPI
+  if (with_sink) {
+    printf("Error: sink particles are not available yet with MPI.\n");
+    return 1;
+  }
+#endif
+
   /* The CPU frequency is a long long, so we need to parse that ourselves. */
   if (cpufreqarg != NULL) {
     if (sscanf(cpufreqarg, "%llu", &cpufreq) != 1) {
@@ -623,6 +634,7 @@ int main(int argc, char *argv[]) {
   if (myrank == 0) {
     message("sizeof(part)        is %4zi bytes.", sizeof(struct part));
     message("sizeof(xpart)       is %4zi bytes.", sizeof(struct xpart));
+    message("sizeof(sink)        is %4zi bytes.", sizeof(struct sink));
     message("sizeof(spart)       is %4zi bytes.", sizeof(struct spart));
     message("sizeof(bpart)       is %4zi bytes.", sizeof(struct bpart));
     message("sizeof(gpart)       is %4zi bytes.", sizeof(struct gpart));
@@ -1004,33 +1016,37 @@ int main(int argc, char *argv[]) {
     fflush(stdout);
 
     /* Get ready to read particles of all kinds */
-    size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0, Nspart = 0, Nbpart = 0;
+    size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0;
+    size_t Nspart = 0, Nbpart = 0, Nsink = 0;
     double dim[3] = {0., 0., 0.};
 
     if (myrank == 0) clocks_gettime(&tic);
 #if defined(HAVE_HDF5)
 #if defined(WITH_MPI)
 #if defined(HAVE_PARALLEL_HDF5)
-    read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
-                     &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
-                     &flag_entropy_ICs, with_hydro, with_gravity, with_stars,
-                     with_black_holes, with_cosmology, cleanup_h,
-                     cleanup_sqrt_a, cosmo.h, cosmo.a, myrank, nr_nodes,
-                     MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, dry_run);
+    read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sinks, &sparts,
+                     &bparts, &Ngas, &Ngpart, &Ngpart_background, &Nsink,
+                     &Nspart, &Nbpart, &flag_entropy_ICs, with_hydro,
+                     with_gravity, with_sink, with_stars, with_black_holes,
+                     with_cosmology, cleanup_h, cleanup_sqrt_a, cosmo.h,
+                     cosmo.a, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL,
+                     nr_threads, dry_run);
 #else
-    read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
-                   &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
-                   &flag_entropy_ICs, with_hydro, with_gravity, with_stars,
-                   with_black_holes, with_cosmology, cleanup_h, cleanup_sqrt_a,
-                   cosmo.h, cosmo.a, myrank, nr_nodes, MPI_COMM_WORLD,
-                   MPI_INFO_NULL, nr_threads, dry_run);
+    read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sinks, &sparts,
+                   &bparts, &Ngas, &Ngpart, &Ngpart_background, &Nsink, &Nspart,
+                   &Nbpart, &flag_entropy_ICs, with_hydro, with_gravity,
+                   with_sink, with_stars, with_black_holes, with_cosmology,
+                   cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank,
+                   nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
+                   dry_run);
 #endif
 #else
-    read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
-                   &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
-                   &flag_entropy_ICs, with_hydro, with_gravity, with_stars,
-                   with_black_holes, with_cosmology, cleanup_h, cleanup_sqrt_a,
-                   cosmo.h, cosmo.a, nr_threads, dry_run);
+    read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sinks, &sparts,
+                   &bparts, &Ngas, &Ngpart, &Ngpart_background, &Nsink, &Nspart,
+                   &Nbpart, &flag_entropy_ICs, with_hydro, with_gravity,
+                   with_sink, with_stars, with_black_holes, with_cosmology,
+                   cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads,
+                   dry_run);
 #endif
 #endif
 
@@ -1061,22 +1077,27 @@ int main(int argc, char *argv[]) {
       for (size_t k = 0; k < Ngpart; ++k)
         if (gparts[k].type == swift_type_gas) error("Linking problem");
     }
+    if (!with_sink && !dry_run) {
+      for (size_t k = 0; k < Ngpart; ++k)
+        if (gparts[k].type == swift_type_sink) error("Linking problem");
+    }
 
     /* Check that the other links are correctly set */
     if (!dry_run)
-      part_verify_links(parts, gparts, sparts, bparts, Ngas, Ngpart, Nspart,
-                        Nbpart, 1);
+      part_verify_links(parts, gparts, sinks, sparts, bparts, Ngas, Ngpart,
+                        Nsink, Nspart, Nbpart, 1);
 #endif
 
     /* Get the total number of particles across all nodes. */
     long long N_total[swift_type_count + 1] = {0};
-    long long Nbaryons = Ngas + Nspart + Nbpart;
+    long long Nbaryons = Ngas + Nspart + Nbpart + Nsink;
 #if defined(WITH_MPI)
     long long N_long[swift_type_count + 1] = {0};
     N_long[swift_type_gas] = Ngas;
     N_long[swift_type_dark_matter] =
         with_gravity ? Ngpart - Ngpart_background - Nbaryons : 0;
     N_long[swift_type_dark_matter_background] = Ngpart_background;
+    N_long[swift_type_sink] = Nsink;
     N_long[swift_type_stars] = Nspart;
     N_long[swift_type_black_hole] = Nbpart;
     N_long[swift_type_count] = Ngpart;
@@ -1087,6 +1108,7 @@ int main(int argc, char *argv[]) {
     N_total[swift_type_dark_matter] =
         with_gravity ? Ngpart - Ngpart_background - Nbaryons : 0;
     N_total[swift_type_dark_matter_background] = Ngpart_background;
+    N_total[swift_type_sink] = Nsink;
     N_total[swift_type_stars] = Nspart;
     N_total[swift_type_black_hole] = Nbpart;
     N_total[swift_type_count] = Ngpart;
@@ -1094,17 +1116,19 @@ int main(int argc, char *argv[]) {
 
     if (myrank == 0)
       message(
-          "Read %lld gas particles, %lld stars particles, %lld black hole "
+          "Read %lld gas particles, %lld sink particles, %lld stars particles, "
+          "%lld black hole "
           "particles, %lld DM particles and %lld DM background particles from "
           "the ICs.",
-          N_total[swift_type_gas], N_total[swift_type_stars],
-          N_total[swift_type_black_hole], N_total[swift_type_dark_matter],
+          N_total[swift_type_gas], N_total[swift_type_sink],
+          N_total[swift_type_stars], N_total[swift_type_black_hole],
+          N_total[swift_type_dark_matter],
           N_total[swift_type_dark_matter_background]);
 
     const int with_DM_particles = N_total[swift_type_dark_matter] > 0;
     const int with_baryon_particles =
         (N_total[swift_type_gas] + N_total[swift_type_stars] +
-             N_total[swift_type_black_hole] >
+             N_total[swift_type_black_hole] + N_total[swift_type_sink] >
          0) ||
         (with_DM_particles && generate_gas_in_ics);
 
@@ -1114,8 +1138,8 @@ int main(int argc, char *argv[]) {
 
     /* Initialize the space with these data. */
     if (myrank == 0) clocks_gettime(&tic);
-    space_init(&s, params, &cosmo, dim, &hydro_properties, parts, gparts,
-               sparts, bparts, Ngas, Ngpart, Nspart, Nbpart, periodic,
+    space_init(&s, params, &cosmo, dim, &hydro_properties, parts, gparts, sinks,
+               sparts, bparts, Ngas, Ngpart, Nsink, Nspart, Nbpart, periodic,
                replicate, generate_gas_in_ics, with_hydro, with_self_gravity,
                with_star_formation, with_DM_background_particles, talking,
                dry_run, nr_nodes);
@@ -1163,13 +1187,14 @@ int main(int argc, char *argv[]) {
       space_check_cosmology(&s, &cosmo, myrank);
 
     /* Also update the total counts (in case of changes due to replication) */
-    Nbaryons = s.nr_parts + s.nr_sparts + s.nr_bparts;
+    Nbaryons = s.nr_parts + s.nr_sparts + s.nr_bparts + s.nr_sinks;
 #if defined(WITH_MPI)
     N_long[swift_type_gas] = s.nr_parts;
     N_long[swift_type_dark_matter] =
         with_gravity ? s.nr_gparts - Ngpart_background - Nbaryons : 0;
     N_long[swift_type_count] = s.nr_gparts;
     N_long[swift_type_stars] = s.nr_sparts;
+    N_long[swift_type_sink] = s.nr_sinks;
     N_long[swift_type_black_hole] = s.nr_bparts;
     MPI_Allreduce(&N_long, &N_total, swift_type_count + 1, MPI_LONG_LONG_INT,
                   MPI_SUM, MPI_COMM_WORLD);
@@ -1179,6 +1204,7 @@ int main(int argc, char *argv[]) {
         with_gravity ? s.nr_gparts - Ngpart_background - Nbaryons : 0;
     N_total[swift_type_count] = s.nr_gparts;
     N_total[swift_type_stars] = s.nr_sparts;
+    N_total[swift_type_sink] = s.nr_sinks;
     N_total[swift_type_black_hole] = s.nr_bparts;
 #endif
 
@@ -1191,6 +1217,7 @@ int main(int argc, char *argv[]) {
               s.cdim[1], s.cdim[2]);
       message("%zi parts in %i cells.", s.nr_parts, s.tot_cells);
       message("%zi gparts in %i cells.", s.nr_gparts, s.tot_cells);
+      message("%zi sinks in %i cells.", s.nr_sinks, s.tot_cells);
       message("%zi sparts in %i cells.", s.nr_sparts, s.tot_cells);
       message("%zi bparts in %i cells.", s.nr_bparts, s.tot_cells);
       message("maximum depth is %d.", s.maxdepth);
@@ -1247,12 +1274,13 @@ int main(int argc, char *argv[]) {
     if (with_fof) engine_policies |= engine_policy_fof;
     if (with_logger) engine_policies |= engine_policy_logger;
     if (with_line_of_sight) engine_policies |= engine_policy_line_of_sight;
+    if (with_sink) engine_policies |= engine_policy_sink;
 
     /* Initialize the engine with the space and policies. */
     if (myrank == 0) clocks_gettime(&tic);
     engine_init(&e, &s, params, output_options, N_total[swift_type_gas],
-                N_total[swift_type_count], N_total[swift_type_stars],
-                N_total[swift_type_black_hole],
+                N_total[swift_type_count], N_total[swift_type_sink],
+                N_total[swift_type_stars], N_total[swift_type_black_hole],
                 N_total[swift_type_dark_matter_background], engine_policies,
                 talking, &reparttype, &us, &prog_const, &cosmo,
                 &hydro_properties, &entropy_floor, &gravity_properties,
@@ -1279,10 +1307,13 @@ int main(int argc, char *argv[]) {
       const long long N_DM = N_total[swift_type_dark_matter] +
                              N_total[swift_type_dark_matter_background];
       message(
-          "Running on %lld gas particles, %lld stars particles %lld black "
-          "hole particles and %lld DM particles (%lld gravity particles)",
-          N_total[swift_type_gas], N_total[swift_type_stars],
-          N_total[swift_type_black_hole], N_DM, N_total[swift_type_count]);
+          "Running on %lld gas particles, %lld sink particles, %lld stars "
+          "particles, "
+          "%lld black hole particles and %lld DM particles (%lld gravity "
+          "particles)",
+          N_total[swift_type_gas], N_total[swift_type_sink],
+          N_total[swift_type_stars], N_total[swift_type_black_hole], N_DM,
+          N_total[swift_type_count]);
       message(
           "from t=%.3e until t=%.3e with %d ranks, %d threads / rank and %d "
           "task queues / rank (dt_min=%.3e, dt_max=%.3e)...",
diff --git a/examples/main_fof.c b/examples/main_fof.c
index 25eef7b6e5bd4795566732a58435705318f5c9c6..047a0f358c5c9e78a102069f906a99081a5ed22b 100644
--- a/examples/main_fof.c
+++ b/examples/main_fof.c
@@ -95,6 +95,7 @@ int main(int argc, char *argv[]) {
   struct phys_const prog_const;
   struct space s;
   struct spart *sparts = NULL;
+  struct sink *sinks = NULL;
   struct bpart *bparts = NULL;
   struct unit_system us;
 
@@ -141,6 +142,7 @@ int main(int argc, char *argv[]) {
   int dump_threadpool = 0;
   int with_fp_exceptions = 0;
   int with_cosmology = 0;
+  int with_sinks = 0;
   int with_stars = 0;
   int with_black_holes = 0;
   int with_hydro = 0;
@@ -164,6 +166,8 @@ int main(int argc, char *argv[]) {
                   "Run with cosmological information.", NULL, 0, 0),
       OPT_BOOLEAN(0, "hydro", &with_hydro, "Read gas particles from the ICs.",
                   NULL, 0, 0),
+      OPT_BOOLEAN(0, "sinks", &with_sinks, "Read sinks from the ICs.", NULL, 0,
+                  0),
       OPT_BOOLEAN(0, "stars", &with_stars, "Read stars from the ICs.", NULL, 0,
                   0),
       OPT_BOOLEAN(0, "black-holes", &with_black_holes,
@@ -317,6 +321,7 @@ int main(int argc, char *argv[]) {
   if (myrank == 0) {
     message("sizeof(part)        is %4zi bytes.", sizeof(struct part));
     message("sizeof(xpart)       is %4zi bytes.", sizeof(struct xpart));
+    message("sizeof(sink)        is %4zi bytes.", sizeof(struct sink));
     message("sizeof(spart)       is %4zi bytes.", sizeof(struct spart));
     message("sizeof(gpart)       is %4zi bytes.", sizeof(struct gpart));
     message("sizeof(multipole)   is %4zi bytes.", sizeof(struct multipole));
@@ -447,34 +452,37 @@ int main(int argc, char *argv[]) {
 
   /* Get ready to read particles of all kinds */
   int flag_entropy_ICs = 0;
-  size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0, Nspart = 0, Nbpart = 0;
+  size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0;
+  size_t Nsink = 0, Nspart = 0, Nbpart = 0;
   double dim[3] = {0., 0., 0.};
 
   if (myrank == 0) clocks_gettime(&tic);
 #if defined(HAVE_HDF5)
 #if defined(WITH_MPI)
 #if defined(HAVE_PARALLEL_HDF5)
-  read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
-                   &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
-                   &flag_entropy_ICs, with_hydro,
-                   /*with_grav=*/1, with_stars, with_black_holes,
+  read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sinks, &sparts,
+                   &bparts, &Ngas, &Ngpart, &Ngpart_background, &Nsink, &Nspart,
+                   &Nbpart, &flag_entropy_ICs, with_hydro,
+                   /*with_grav=*/1, with_sinks, with_stars, with_black_holes,
                    with_cosmology, cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a,
                    myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
                    /*dry_run=*/0);
 #else
-  read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas,
-                 &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
-                 &flag_entropy_ICs, with_hydro,
-                 /*with_grav=*/1, with_stars, with_black_holes, with_cosmology,
-                 cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank, nr_nodes,
-                 MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, /*dry_run=*/0);
+  read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sinks, &sparts,
+                 &bparts, &Ngas, &Ngpart, &Ngpart_background, &Nsink, &Nspart,
+                 &Nbpart, &flag_entropy_ICs, with_hydro,
+                 /*with_grav=*/1, with_sinks, with_stars, with_black_holes,
+                 with_cosmology, cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a,
+                 myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
+                 /*dry_run=*/0);
 #endif
 #else
-  read_ic_single(
-      ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas, &Ngpart,
-      &Ngpart_background, &Nspart, &Nbpart, &flag_entropy_ICs, with_hydro,
-      /*with_grav=*/1, with_stars, with_black_holes, with_cosmology, cleanup_h,
-      cleanup_sqrt_a, cosmo.h, cosmo.a, nr_threads, /*dry_run=*/0);
+  read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sinks, &sparts,
+                 &bparts, &Ngas, &Ngpart, &Ngpart_background, &Nsink, &Nspart,
+                 &Nbpart, &flag_entropy_ICs, with_hydro,
+                 /*with_grav=*/1, with_sinks, with_stars, with_black_holes,
+                 with_cosmology, cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a,
+                 nr_threads, /*dry_run=*/0);
 #endif
 #endif
   if (myrank == 0) {
@@ -498,6 +506,7 @@ int main(int argc, char *argv[]) {
   N_long[swift_type_gas] = Ngas;
   N_long[swift_type_dark_matter] = Ngpart - Ngpart_background - Nbaryons;
   N_long[swift_type_dark_matter_background] = Ngpart_background;
+  N_long[swift_type_sink] = Nsink;
   N_long[swift_type_stars] = Nspart;
   N_long[swift_type_black_hole] = Nbpart;
   N_long[swift_type_count] = Ngpart;
@@ -507,6 +516,7 @@ int main(int argc, char *argv[]) {
   N_total[swift_type_gas] = Ngas;
   N_total[swift_type_dark_matter] = Ngpart - Ngpart_background - Nbaryons;
   N_total[swift_type_dark_matter_background] = Ngpart_background;
+  N_total[swift_type_sink] = Nsink;
   N_total[swift_type_stars] = Nspart;
   N_total[swift_type_black_hole] = Nbpart;
   N_total[swift_type_count] = Ngpart;
@@ -514,17 +524,19 @@ int main(int argc, char *argv[]) {
 
   if (myrank == 0)
     message(
-        "Read %lld gas particles, %lld stars particles, %lld black hole "
+        "Read %lld gas particles, %lld sink particles, %lld stars particles, "
+        "%lld black hole "
         "particles, %lld DM particles and %lld DM background particles from "
         "the ICs.",
-        N_total[swift_type_gas], N_total[swift_type_stars],
-        N_total[swift_type_black_hole], N_total[swift_type_dark_matter],
+        N_total[swift_type_gas], N_total[swift_type_sink],
+        N_total[swift_type_stars], N_total[swift_type_black_hole],
+        N_total[swift_type_dark_matter],
         N_total[swift_type_dark_matter_background]);
 
   const int with_DM_particles = N_total[swift_type_dark_matter] > 0;
   const int with_baryon_particles =
-      (N_total[swift_type_gas] + N_total[swift_type_stars] +
-       N_total[swift_type_black_hole]) > 0;
+      (N_total[swift_type_gas] + N_total[swift_type_sink] +
+       N_total[swift_type_stars] + N_total[swift_type_black_hole]) > 0;
 
   /* Do we have background DM particles? */
   const int with_DM_background_particles =
@@ -533,7 +545,8 @@ int main(int argc, char *argv[]) {
   /* Initialize the space with these data. */
   if (myrank == 0) clocks_gettime(&tic);
   space_init(&s, params, &cosmo, dim, /*hydro_props=*/NULL, parts, gparts,
-             sparts, bparts, Ngas, Ngpart, Nspart, Nbpart, periodic, replicate,
+             sinks, sparts, bparts, Ngas, Ngpart, Nsink, Nspart, Nbpart,
+             periodic, replicate,
              /*generate_gas_in_ics=*/0, /*hydro=*/N_total[0] > 0, /*gravity=*/1,
              /*with_star_formation=*/0, with_DM_background_particles, talking,
              /*dry_run=*/0, nr_nodes);
@@ -570,6 +583,7 @@ int main(int argc, char *argv[]) {
   N_long[swift_type_gas] = s.nr_parts;
   N_long[swift_type_dark_matter] = s.nr_gparts - Ngpart_background - Nbaryons;
   N_long[swift_type_count] = s.nr_gparts;
+  N_long[swift_type_sink] = s.nr_sinks;
   N_long[swift_type_stars] = s.nr_sparts;
   N_long[swift_type_black_hole] = s.nr_bparts;
   MPI_Allreduce(&N_long, &N_total, swift_type_count + 1, MPI_LONG_LONG_INT,
@@ -578,6 +592,7 @@ int main(int argc, char *argv[]) {
   N_total[swift_type_gas] = s.nr_parts;
   N_total[swift_type_dark_matter] = s.nr_gparts - Ngpart_background - Nbaryons;
   N_total[swift_type_count] = s.nr_gparts;
+  N_total[swift_type_sink] = s.nr_sinks;
   N_total[swift_type_stars] = s.nr_sparts;
   N_total[swift_type_black_hole] = s.nr_bparts;
 #endif
@@ -591,6 +606,7 @@ int main(int argc, char *argv[]) {
             s.cdim[1], s.cdim[2]);
     message("%zi parts in %i cells.", s.nr_parts, s.tot_cells);
     message("%zi gparts in %i cells.", s.nr_gparts, s.tot_cells);
+    message("%zi sinks in %i cells.", s.nr_sinks, s.tot_cells);
     message("%zi sparts in %i cells.", s.nr_sparts, s.tot_cells);
     message("maximum depth is %d.", s.maxdepth);
     fflush(stdout);
@@ -605,8 +621,8 @@ int main(int argc, char *argv[]) {
   /* Initialize the engine with the space and policies. */
   if (myrank == 0) clocks_gettime(&tic);
   engine_init(&e, &s, params, output_options, N_total[swift_type_gas],
-              N_total[swift_type_count], N_total[swift_type_stars],
-              N_total[swift_type_black_hole],
+              N_total[swift_type_count], N_total[swift_type_sink],
+              N_total[swift_type_stars], N_total[swift_type_black_hole],
               N_total[swift_type_dark_matter_background], engine_policies,
               talking, &reparttype, &us, &prog_const, &cosmo,
               /*hydro_properties=*/NULL, /*entropy_floor=*/NULL,
@@ -630,10 +646,12 @@ int main(int argc, char *argv[]) {
     const long long N_DM = N_total[swift_type_dark_matter] +
                            N_total[swift_type_dark_matter_background];
     message(
-        "Running FOF on %lld gas particles, %lld stars particles %lld black "
+        "Running FOF on %lld gas particles, %lld sink particles, %lld stars "
+        "particles %lld black "
         "hole particles and %lld DM particles (%lld gravity particles)",
-        N_total[swift_type_gas], N_total[swift_type_stars],
-        N_total[swift_type_black_hole], N_DM, N_total[swift_type_count]);
+        N_total[swift_type_gas], N_total[swift_type_sink],
+        N_total[swift_type_stars], N_total[swift_type_black_hole], N_DM,
+        N_total[swift_type_count]);
     message(
         "from t=%.3e until t=%.3e with %d ranks, %d threads / rank and %d "
         "task queues / rank (dt_min=%.3e, dt_max=%.3e)...",
diff --git a/src/cell.c b/src/cell.c
index 35c98b9219a601b1c97674ede28e73ac736258bd..5ca50b378a0126b88c9591eef09d4230a459f6e6 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1648,6 +1648,8 @@ void cell_bunlocktree(struct cell *c) {
  * @param bparts_offset Offset of the cell bparts array relative to the
  *        space's bparts array, i.e. c->black_holes.parts -
  * s->black_holes.parts.
+ * @param sinks_offset Offset of the cell sink array relative to the
+ *        space's sink array, i.e. c->sinks.parts - s->sinks.parts.
  * @param buff A buffer with at least max(c->hydro.count, c->grav.count)
  * entries, used for sorting indices.
  * @param sbuff A buffer with at least max(c->stars.count, c->grav.count)
@@ -1656,18 +1658,23 @@ void cell_bunlocktree(struct cell *c) {
  * entries, used for sorting indices for the sparts.
  * @param gbuff A buffer with at least max(c->hydro.count, c->grav.count)
  * entries, used for sorting indices for the gparts.
+ * @param sinkbuff A buffer with at least max(c->sinks.count, c->grav.count)
+ * entries, used for sorting indices for the sinks.
  */
 void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
-                ptrdiff_t bparts_offset, struct cell_buff *buff,
-                struct cell_buff *sbuff, struct cell_buff *bbuff,
-                struct cell_buff *gbuff) {
+                ptrdiff_t bparts_offset, ptrdiff_t sinks_offset,
+                struct cell_buff *buff, struct cell_buff *sbuff,
+                struct cell_buff *bbuff, struct cell_buff *gbuff,
+                struct cell_buff *sinkbuff) {
   const int count = c->hydro.count, gcount = c->grav.count,
-            scount = c->stars.count, bcount = c->black_holes.count;
+            scount = c->stars.count, bcount = c->black_holes.count,
+            sink_count = c->sinks.count;
   struct part *parts = c->hydro.parts;
   struct xpart *xparts = c->hydro.xparts;
   struct gpart *gparts = c->grav.parts;
   struct spart *sparts = c->stars.parts;
   struct bpart *bparts = c->black_holes.parts;
+  struct sink *sinks = c->sinks.parts;
   const double pivot[3] = {c->loc[0] + c->width[0] / 2,
                            c->loc[1] + c->width[1] / 2,
                            c->loc[2] + c->width[2] / 2};
@@ -1696,6 +1703,11 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
         bbuff[k].x[2] != bparts[k].x[2])
       error("Inconsistent bbuff contents.");
   }
+  for (int k = 0; k < sink_count; k++) {
+    if (sinkbuff[k].x[0] != sinks[k].x[0] ||
+        sinkbuff[k].x[1] != sinks[k].x[1] || sinkbuff[k].x[2] != sinks[k].x[2])
+      error("Inconsistent sinkbuff contents.");
+  }
 #endif /* SWIFT_DEBUG_CHECKS */
 
   /* Fill the buffer with the indices. */
@@ -1925,6 +1937,60 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
     c->progeny[k]->black_holes.parts = &c->black_holes.parts[bucket_offset[k]];
   }
 
+  /* Now do the same song and dance for the sinks. */
+  for (int k = 0; k < 8; k++) bucket_count[k] = 0;
+
+  /* Fill the buffer with the indices. */
+  for (int k = 0; k < sink_count; k++) {
+    const int bid = (sinkbuff[k].x[0] > pivot[0]) * 4 +
+                    (sinkbuff[k].x[1] > pivot[1]) * 2 +
+                    (sinkbuff[k].x[2] > pivot[2]);
+    bucket_count[bid]++;
+    sinkbuff[k].ind = bid;
+  }
+
+  /* Set the buffer offsets. */
+  bucket_offset[0] = 0;
+  for (int k = 1; k <= 8; k++) {
+    bucket_offset[k] = bucket_offset[k - 1] + bucket_count[k - 1];
+    bucket_count[k - 1] = 0;
+  }
+
+  /* Run through the buckets, and swap particles to their correct spot. */
+  for (int bucket = 0; bucket < 8; bucket++) {
+    for (int k = bucket_offset[bucket] + bucket_count[bucket];
+         k < bucket_offset[bucket + 1]; k++) {
+      int bid = sinkbuff[k].ind;
+      if (bid != bucket) {
+        struct sink sink = sinks[k];
+        struct cell_buff temp_buff = sinkbuff[k];
+        while (bid != bucket) {
+          int j = bucket_offset[bid] + bucket_count[bid]++;
+          while (sinkbuff[j].ind == bid) {
+            j++;
+            bucket_count[bid]++;
+          }
+          memswap(&sinks[j], &sink, sizeof(struct sink));
+          memswap(&sinkbuff[j], &temp_buff, sizeof(struct cell_buff));
+          if (sinks[j].gpart)
+            sinks[j].gpart->id_or_neg_offset = -(j + sinks_offset);
+          bid = temp_buff.ind;
+        }
+        sinks[k] = sink;
+        sinkbuff[k] = temp_buff;
+        if (sinks[k].gpart)
+          sinks[k].gpart->id_or_neg_offset = -(k + sinks_offset);
+      }
+      bucket_count[bid]++;
+    }
+  }
+
+  /* Store the counts and offsets. */
+  for (int k = 0; k < 8; k++) {
+    c->progeny[k]->sinks.count = bucket_count[k];
+    c->progeny[k]->sinks.parts = &c->sinks.parts[bucket_offset[k]];
+  }
+
   /* Finally, do the same song and dance for the gparts. */
   for (int k = 0; k < 8; k++) bucket_count[k] = 0;
 
@@ -1965,6 +2031,9 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
           } else if (gparts[j].type == swift_type_stars) {
             sparts[-gparts[j].id_or_neg_offset - sparts_offset].gpart =
                 &gparts[j];
+          } else if (gparts[j].type == swift_type_sink) {
+            sinks[-gparts[j].id_or_neg_offset - sinks_offset].gpart =
+                &gparts[j];
           } else if (gparts[j].type == swift_type_black_hole) {
             bparts[-gparts[j].id_or_neg_offset - bparts_offset].gpart =
                 &gparts[j];
@@ -1978,6 +2047,8 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
         } else if (gparts[k].type == swift_type_stars) {
           sparts[-gparts[k].id_or_neg_offset - sparts_offset].gpart =
               &gparts[k];
+        } else if (gparts[k].type == swift_type_sink) {
+          sinks[-gparts[k].id_or_neg_offset - sinks_offset].gpart = &gparts[k];
         } else if (gparts[k].type == swift_type_black_hole) {
           bparts[-gparts[k].id_or_neg_offset - bparts_offset].gpart =
               &gparts[k];
diff --git a/src/cell.h b/src/cell.h
index 766f3e20715bb62b694e9f6c16e7d65ba55316fa..1ea9e25f4181e11264bdb501040fb7a0292a1762 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -752,6 +752,19 @@ struct cell {
 
   } black_holes;
 
+  /*! Sink particles variables */
+  struct {
+
+    /*! Pointer to the #sink data. */
+    struct sink *parts;
+
+    /*! Nr of #sink in this cell. */
+    int count;
+
+    /*! Nr of #sink this cell can hold after addition of new one. */
+    int count_total;
+  } sinks;
+
 #ifdef WITH_MPI
   /*! MPI variables */
   struct {
@@ -843,9 +856,10 @@ struct cell {
 
 /* Function prototypes. */
 void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
-                ptrdiff_t bparts_offset, struct cell_buff *buff,
-                struct cell_buff *sbuff, struct cell_buff *bbuff,
-                struct cell_buff *gbuff);
+                ptrdiff_t bparts_offset, ptrdiff_t sinks_offset,
+                struct cell_buff *buff, struct cell_buff *sbuff,
+                struct cell_buff *bbuff, struct cell_buff *gbuff,
+                struct cell_buff *sinkbuff);
 void cell_sanitize(struct cell *c, int treated);
 int cell_locktree(struct cell *c);
 void cell_unlocktree(struct cell *c);
diff --git a/src/common_io.c b/src/common_io.c
index 457016bb700bd0d31d6310e0e9a4bd80eeed1aac..ef66a9ed4df2a904ec1b2c1360b961151b6a8d51 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -42,6 +42,7 @@
 #include "kernel_hydro.h"
 #include "part.h"
 #include "part_type.h"
+#include "sink_io.h"
 #include "star_formation_io.h"
 #include "stars_io.h"
 #include "threadpool.h"
@@ -825,6 +826,19 @@ static long long cell_count_non_inhibited_black_holes(const struct cell* c) {
   return count;
 }
 
+static long long cell_count_non_inhibited_sinks(const struct cell* c) {
+  const int total_count = c->sinks.count;
+  struct sink* sinks = c->sinks.parts;
+  long long count = 0;
+  for (int i = 0; i < total_count; ++i) {
+    if ((sinks[i].time_bin != time_bin_inhibited) &&
+        (sinks[i].time_bin != time_bin_not_created)) {
+      ++count;
+    }
+  }
+  return count;
+}
+
 /**
  * @brief Write a single 1D array to a hdf5 group.
  *
@@ -916,22 +930,24 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const double dim[3],
   /* Count of particles in each cell */
   long long *count_part = NULL, *count_gpart = NULL,
             *count_background_gpart = NULL, *count_spart = NULL,
-            *count_bpart = NULL;
+            *count_bpart = NULL, *count_sink = NULL;
   count_part = (long long*)malloc(nr_cells * sizeof(long long));
   count_gpart = (long long*)malloc(nr_cells * sizeof(long long));
   count_background_gpart = (long long*)malloc(nr_cells * sizeof(long long));
   count_spart = (long long*)malloc(nr_cells * sizeof(long long));
   count_bpart = (long long*)malloc(nr_cells * sizeof(long long));
+  count_sink = (long long*)malloc(nr_cells * sizeof(long long));
 
   /* Global offsets of particles in each cell */
   long long *offset_part = NULL, *offset_gpart = NULL,
             *offset_background_gpart = NULL, *offset_spart = NULL,
-            *offset_bpart = NULL;
+            *offset_bpart = NULL, *offset_sink = NULL;
   offset_part = (long long*)malloc(nr_cells * sizeof(long long));
   offset_gpart = (long long*)malloc(nr_cells * sizeof(long long));
   offset_background_gpart = (long long*)malloc(nr_cells * sizeof(long long));
   offset_spart = (long long*)malloc(nr_cells * sizeof(long long));
   offset_bpart = (long long*)malloc(nr_cells * sizeof(long long));
+  offset_sink = (long long*)malloc(nr_cells * sizeof(long long));
 
   /* Offsets of the 0^th element */
   offset_part[0] = 0;
@@ -939,6 +955,7 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const double dim[3],
   offset_background_gpart[0] = 0;
   offset_spart[0] = 0;
   offset_bpart[0] = 0;
+  offset_sink[0] = 0;
 
   /* Collect the cell information of *local* cells */
   long long local_offset_part = 0;
@@ -946,6 +963,7 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const double dim[3],
   long long local_offset_background_gpart = 0;
   long long local_offset_spart = 0;
   long long local_offset_bpart = 0;
+  long long local_offset_sink = 0;
   for (int i = 0; i < nr_cells; ++i) {
 
     /* Store in which file this cell will be found */
@@ -981,6 +999,7 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const double dim[3],
           cell_count_non_inhibited_background_dark_matter(&cells_top[i]);
       count_spart[i] = cell_count_non_inhibited_stars(&cells_top[i]);
       count_bpart[i] = cell_count_non_inhibited_black_holes(&cells_top[i]);
+      count_sink[i] = cell_count_non_inhibited_sinks(&cells_top[i]);
 
       /* Offsets including the global offset of all particles on this MPI rank
        * Note that in the distributed case, the global offsets are 0 such that
@@ -994,12 +1013,14 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const double dim[3],
       offset_spart[i] = local_offset_spart + global_offsets[swift_type_stars];
       offset_bpart[i] =
           local_offset_bpart + global_offsets[swift_type_black_hole];
+      offset_sink[i] = local_offset_sink + global_offsets[swift_type_sink];
 
       local_offset_part += count_part[i];
       local_offset_gpart += count_gpart[i];
       local_offset_background_gpart += count_background_gpart[i];
       local_offset_spart += count_spart[i];
       local_offset_bpart += count_bpart[i];
+      local_offset_sink += count_sink[i];
 
     } else {
 
@@ -1014,12 +1035,14 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const double dim[3],
       count_background_gpart[i] = 0;
       count_spart[i] = 0;
       count_bpart[i] = 0;
+      count_sink[i] = 0;
 
       offset_part[i] = 0;
       offset_gpart[i] = 0;
       offset_background_gpart[i] = 0;
       offset_spart[i] = 0;
       offset_bpart[i] = 0;
+      offset_sink[i] = 0;
     }
   }
 
@@ -1032,6 +1055,8 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const double dim[3],
                 MPI_COMM_WORLD);
   MPI_Allreduce(MPI_IN_PLACE, count_background_gpart, nr_cells,
                 MPI_LONG_LONG_INT, MPI_BOR, MPI_COMM_WORLD);
+  MPI_Allreduce(MPI_IN_PLACE, count_sink, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
+                MPI_COMM_WORLD);
   MPI_Allreduce(MPI_IN_PLACE, count_spart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
                 MPI_COMM_WORLD);
   MPI_Allreduce(MPI_IN_PLACE, count_bpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
@@ -1043,6 +1068,8 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const double dim[3],
                 MPI_BOR, MPI_COMM_WORLD);
   MPI_Allreduce(MPI_IN_PLACE, offset_background_gpart, nr_cells,
                 MPI_LONG_LONG_INT, MPI_BOR, MPI_COMM_WORLD);
+  MPI_Allreduce(MPI_IN_PLACE, offset_sink, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
+                MPI_COMM_WORLD);
   MPI_Allreduce(MPI_IN_PLACE, offset_spart, nr_cells, MPI_LONG_LONG_INT,
                 MPI_BOR, MPI_COMM_WORLD);
   MPI_Allreduce(MPI_IN_PLACE, offset_bpart, nr_cells, MPI_LONG_LONG_INT,
@@ -1137,6 +1164,14 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const double dim[3],
                      "PartType2", "counts");
     }
 
+    if (global_counts[swift_type_sink] > 0 && num_fields[swift_type_sink] > 0) {
+      io_write_array(h_grp_files, nr_cells, files, INT, "PartType3", "files");
+      io_write_array(h_grp_offsets, nr_cells, offset_sink, LONGLONG,
+                     "PartType3", "offsets");
+      io_write_array(h_grp_counts, nr_cells, count_sink, LONGLONG, "PartType3",
+                     "counts");
+    }
+
     if (global_counts[swift_type_stars] > 0 &&
         num_fields[swift_type_stars] > 0) {
       io_write_array(h_grp_files, nr_cells, files, INT, "PartType4", "files");
@@ -1168,11 +1203,13 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const double dim[3],
   free(count_background_gpart);
   free(count_spart);
   free(count_bpart);
+  free(count_sink);
   free(offset_part);
   free(offset_gpart);
   free(offset_background_gpart);
   free(offset_spart);
   free(offset_bpart);
+  free(offset_sink);
 }
 
 #endif /* HAVE_HDF5 */
@@ -1558,6 +1595,86 @@ void io_convert_bpart_l_mapper(void* restrict temp, int N,
     props.convert_bpart_l(e, bparts + delta + i, &temp_l[i * dim]);
 }
 
+/**
+ * @brief Mapper function to copy #sink into a buffer of floats using a
+ * conversion function.
+ */
+void io_convert_sink_f_mapper(void* restrict temp, int N,
+                              void* restrict extra_data) {
+
+  const struct io_props props = *((const struct io_props*)extra_data);
+  const struct sink* restrict sinks = props.sinks;
+  const struct engine* e = props.e;
+  const size_t dim = props.dimension;
+
+  /* How far are we with this chunk? */
+  float* restrict temp_f = (float*)temp;
+  const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
+
+  for (int i = 0; i < N; i++)
+    props.convert_sink_f(e, sinks + delta + i, &temp_f[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #sink into a buffer of ints using a
+ * conversion function.
+ */
+void io_convert_sink_i_mapper(void* restrict temp, int N,
+                              void* restrict extra_data) {
+
+  const struct io_props props = *((const struct io_props*)extra_data);
+  const struct sink* restrict sinks = props.sinks;
+  const struct engine* e = props.e;
+  const size_t dim = props.dimension;
+
+  /* How far are we with this chunk? */
+  int* restrict temp_i = (int*)temp;
+  const ptrdiff_t delta = (temp_i - props.start_temp_i) / dim;
+
+  for (int i = 0; i < N; i++)
+    props.convert_sink_i(e, sinks + delta + i, &temp_i[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #sink into a buffer of doubles using a
+ * conversion function.
+ */
+void io_convert_sink_d_mapper(void* restrict temp, int N,
+                              void* restrict extra_data) {
+
+  const struct io_props props = *((const struct io_props*)extra_data);
+  const struct sink* restrict sinks = props.sinks;
+  const struct engine* e = props.e;
+  const size_t dim = props.dimension;
+
+  /* How far are we with this chunk? */
+  double* restrict temp_d = (double*)temp;
+  const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
+
+  for (int i = 0; i < N; i++)
+    props.convert_sink_d(e, sinks + delta + i, &temp_d[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #sink into a buffer of doubles using a
+ * conversion function.
+ */
+void io_convert_sink_l_mapper(void* restrict temp, int N,
+                              void* restrict extra_data) {
+
+  const struct io_props props = *((const struct io_props*)extra_data);
+  const struct sink* restrict sinks = props.sinks;
+  const struct engine* e = props.e;
+  const size_t dim = props.dimension;
+
+  /* How far are we with this chunk? */
+  long long* restrict temp_l = (long long*)temp;
+  const ptrdiff_t delta = (temp_l - props.start_temp_l) / dim;
+
+  for (int i = 0; i < N; i++)
+    props.convert_sink_l(e, sinks + delta + i, &temp_l[i * dim]);
+}
+
 /**
  * @brief Copy the particle data into a temporary buffer ready for i/o.
  *
@@ -1735,6 +1852,54 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_spart_l_mapper, temp_l, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
+    } else if (props.convert_sink_f != NULL) {
+
+      /* Prepare some parameters */
+      float* temp_f = (float*)temp;
+      props.start_temp_f = (float*)temp;
+      props.e = e;
+
+      /* Copy the whole thing into a buffer */
+      threadpool_map((struct threadpool*)&e->threadpool,
+                     io_convert_sink_f_mapper, temp_f, N, copySize,
+                     threadpool_auto_chunk_size, (void*)&props);
+
+    } else if (props.convert_sink_i != NULL) {
+
+      /* Prepare some parameters */
+      int* temp_i = (int*)temp;
+      props.start_temp_i = (int*)temp;
+      props.e = e;
+
+      /* Copy the whole thing into a buffer */
+      threadpool_map((struct threadpool*)&e->threadpool,
+                     io_convert_sink_i_mapper, temp_i, N, copySize,
+                     threadpool_auto_chunk_size, (void*)&props);
+
+    } else if (props.convert_sink_d != NULL) {
+
+      /* Prepare some parameters */
+      double* temp_d = (double*)temp;
+      props.start_temp_d = (double*)temp;
+      props.e = e;
+
+      /* Copy the whole thing into a buffer */
+      threadpool_map((struct threadpool*)&e->threadpool,
+                     io_convert_sink_d_mapper, temp_d, N, copySize,
+                     threadpool_auto_chunk_size, (void*)&props);
+
+    } else if (props.convert_sink_l != NULL) {
+
+      /* Prepare some parameters */
+      long long* temp_l = (long long*)temp;
+      props.start_temp_l = (long long*)temp;
+      props.e = e;
+
+      /* Copy the whole thing into a buffer */
+      threadpool_map((struct threadpool*)&e->threadpool,
+                     io_convert_sink_l_mapper, temp_l, N, copySize,
+                     threadpool_auto_chunk_size, (void*)&props);
+
     } else if (props.convert_bpart_f != NULL) {
 
       /* Prepare some parameters */
@@ -1898,10 +2063,12 @@ struct duplication_data {
   struct gpart* gparts;
   struct spart* sparts;
   struct bpart* bparts;
+  struct sink* sinks;
   int Ndm;
   int Ngas;
   int Nstars;
   int Nblackholes;
+  int Nsinks;
 };
 
 void io_duplicate_hydro_gparts_mapper(void* restrict data, int Ngas,
@@ -2018,6 +2185,64 @@ void io_duplicate_stars_gparts(struct threadpool* tp,
                  sizeof(struct spart), threadpool_auto_chunk_size, &data);
 }
 
+void io_duplicate_sinks_gparts_mapper(void* restrict data, int Nsinks,
+                                      void* restrict extra_data) {
+
+  struct duplication_data* temp = (struct duplication_data*)extra_data;
+  const int Ndm = temp->Ndm;
+  struct sink* sinks = (struct sink*)data;
+  const ptrdiff_t offset = sinks - temp->sinks;
+  struct gpart* gparts = temp->gparts + offset;
+
+  for (int i = 0; i < Nsinks; ++i) {
+
+    /* Duplicate the crucial information */
+    gparts[i + Ndm].x[0] = sinks[i].x[0];
+    gparts[i + Ndm].x[1] = sinks[i].x[1];
+    gparts[i + Ndm].x[2] = sinks[i].x[2];
+
+    gparts[i + Ndm].v_full[0] = sinks[i].v[0];
+    gparts[i + Ndm].v_full[1] = sinks[i].v[1];
+    gparts[i + Ndm].v_full[2] = sinks[i].v[2];
+
+    gparts[i + Ndm].mass = sinks[i].mass;
+
+    /* Set gpart type */
+    gparts[i + Ndm].type = swift_type_sink;
+
+    /* Link the particles */
+    gparts[i + Ndm].id_or_neg_offset = -(long long)(offset + i);
+    sinks[i].gpart = &gparts[i + Ndm];
+  }
+}
+
+/**
+ * @brief Copy every #sink into the corresponding #gpart and link them.
+ *
+ * This function assumes that the DM particles, gas particles and star particles
+ * are all at the start of the gparts array and adds the sinks particles
+ * afterwards
+ *
+ * @param tp The current #threadpool.
+ * @param sinks The array of #sink freshly read in.
+ * @param gparts The array of #gpart freshly read in with all the DM, gas
+ * and star particles at the start.
+ * @param Nsinks The number of sink particles read in.
+ * @param Ndm The number of DM, gas and star particles read in.
+ */
+void io_duplicate_sinks_gparts(struct threadpool* tp, struct sink* const sinks,
+                               struct gpart* const gparts, size_t Nsinks,
+                               size_t Ndm) {
+
+  struct duplication_data data;
+  data.gparts = gparts;
+  data.sinks = sinks;
+  data.Ndm = Ndm;
+
+  threadpool_map(tp, io_duplicate_sinks_gparts_mapper, sinks, Nsinks,
+                 sizeof(struct sink), threadpool_auto_chunk_size, &data);
+}
+
 void io_duplicate_black_holes_gparts_mapper(void* restrict data,
                                             int Nblackholes,
                                             void* restrict extra_data) {
@@ -2153,6 +2378,40 @@ void io_collect_sparts_to_write(const struct spart* restrict sparts,
           count, Nsparts_written);
 }
 
+/**
+ * @brief Copy every non-inhibited #sink into the sinks_written array.
+ *
+ * @param sinks The array of #sink containing all particles.
+ * @param sinks_written The array of #sink to fill with particles we want to
+ * write.
+ * @param Nsinks The total number of #sink.
+ * @param Nsinks_written The total number of #sink to write.
+ */
+void io_collect_sinks_to_write(const struct sink* restrict sinks,
+                               struct sink* restrict sinks_written,
+                               const size_t Nsinks,
+                               const size_t Nsinks_written) {
+
+  size_t count = 0;
+
+  /* Loop over all parts */
+  for (size_t i = 0; i < Nsinks; ++i) {
+
+    /* And collect the ones that have not been removed */
+    if (sinks[i].time_bin != time_bin_inhibited &&
+        sinks[i].time_bin != time_bin_not_created) {
+
+      sinks_written[count] = sinks[i];
+      count++;
+    }
+  }
+
+  /* Check that everything is fine */
+  if (count != Nsinks_written)
+    error("Collected the wrong number of sink-particles (%zu vs. %zu expected)",
+          count, Nsinks_written);
+}
+
 /**
  * @brief Copy every non-inhibited #bpart into the bparts_written array.
  *
diff --git a/src/common_io.h b/src/common_io.h
index 1c89176e276c8c9d127b610550a8658f6957715c..920974db86dd850ca1883bc959685a8184022104 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -40,6 +40,7 @@ struct velociraptor_gpart_data;
 struct spart;
 struct bpart;
 struct xpart;
+struct sink;
 struct io_props;
 struct engine;
 struct threadpool;
@@ -134,6 +135,10 @@ void io_collect_parts_to_write(const struct part* restrict parts,
                                struct xpart* restrict xparts_written,
                                const size_t Nparts,
                                const size_t Nparts_written);
+void io_collect_sinks_to_write(const struct sink* restrict sinks,
+                               struct sink* restrict sinks_written,
+                               const size_t Nsinks,
+                               const size_t Nsinks_written);
 void io_collect_sparts_to_write(const struct spart* restrict sparts,
                                 struct spart* restrict sparts_written,
                                 const size_t Nsparts,
@@ -167,6 +172,9 @@ void io_duplicate_stars_gparts(struct threadpool* tp,
                                struct spart* const sparts,
                                struct gpart* const gparts, size_t Nstars,
                                size_t Ndm);
+void io_duplicate_sinks_gparts(struct threadpool* tp, struct sink* const sinks,
+                               struct gpart* const gparts, size_t Nsinks,
+                               size_t Ndm);
 void io_duplicate_black_holes_gparts(struct threadpool* tp,
                                      struct bpart* const bparts,
                                      struct gpart* const gparts, size_t Nstars,
diff --git a/src/distributed_io.c b/src/distributed_io.c
index 4a1b9865bb013b5291405b900ef31e92538146f8..96f9667edec513584302704d4a1356e02e738848 100644
--- a/src/distributed_io.c
+++ b/src/distributed_io.c
@@ -54,6 +54,7 @@
 #include "output_options.h"
 #include "part.h"
 #include "part_type.h"
+#include "sink_io.h"
 #include "star_formation_io.h"
 #include "stars_io.h"
 #include "tools.h"
@@ -248,6 +249,7 @@ void write_output_distributed(struct engine* e,
   const struct part* parts = e->s->parts;
   const struct xpart* xparts = e->s->xparts;
   const struct gpart* gparts = e->s->gparts;
+  const struct sink* sinks = e->s->sinks;
   const struct spart* sparts = e->s->sparts;
   const struct bpart* bparts = e->s->bparts;
   struct output_options* output_options = e->output_options;
@@ -267,6 +269,7 @@ void write_output_distributed(struct engine* e,
   /* 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 Nsinks = e->s->nr_sinks;
   const size_t Nstars = e->s->nr_sparts;
   const size_t Nblackholes = e->s->nr_bparts;
 
@@ -281,12 +284,14 @@ void write_output_distributed(struct engine* e,
       e->s->nr_gparts - e->s->nr_inhibited_gparts - e->s->nr_extra_gparts;
   const size_t Ngas_written =
       e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
+  const size_t Nsinks_written =
+      e->s->nr_sinks - e->s->nr_inhibited_sinks - e->s->nr_extra_sinks;
   const size_t Nstars_written =
       e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
   const size_t Nblackholes_written =
       e->s->nr_bparts - e->s->nr_inhibited_bparts - e->s->nr_extra_bparts;
   const size_t Nbaryons_written =
-      Ngas_written + Nstars_written + Nblackholes_written;
+      Ngas_written + Nstars_written + Nblackholes_written + Nsinks_written;
   const size_t Ndm_written =
       Ntot_written > 0 ? Ntot_written - Nbaryons_written - Ndm_background : 0;
 
@@ -331,7 +336,7 @@ void write_output_distributed(struct engine* e,
 
   /* Compute offset in the file and total number of particles */
   const long long N[swift_type_count] = {Ngas_written,   Ndm_written,
-                                         Ndm_background, 0,
+                                         Ndm_background, Nsinks_written,
                                          Nstars_written, Nblackholes_written};
 
   /* Gather the total number of particles to write */
@@ -471,6 +476,7 @@ void write_output_distributed(struct engine* e,
     struct xpart* xparts_written = NULL;
     struct gpart* gparts_written = NULL;
     struct velociraptor_gpart_data* gpart_group_data_written = NULL;
+    struct sink* sinks_written = NULL;
     struct spart* sparts_written = NULL;
     struct bpart* bparts_written = NULL;
 
@@ -632,6 +638,34 @@ void write_output_distributed(struct engine* e,
         }
       } break;
 
+      case swift_type_sink: {
+        if (Nsinks == Nsinks_written) {
+
+          /* No inhibted particles: easy case */
+          Nparticles = Nsinks;
+          sink_write_particles(sinks, list, &num_fields, with_cosmology);
+
+        } else {
+
+          /* Ok, we need to fish out the particles we want */
+          Nparticles = Nsinks_written;
+
+          /* Allocate temporary arrays */
+          if (swift_memalign("sinks_written", (void**)&sinks_written,
+                             sink_align,
+                             Nsinks_written * sizeof(struct sink)) != 0)
+            error("Error while allocating temporary memory for sinks");
+
+          /* Collect the particles we want to write */
+          io_collect_sinks_to_write(sinks, sinks_written, Nsinks,
+                                    Nsinks_written);
+
+          /* Select the fields to write */
+          sink_write_particles(sinks_written, list, &num_fields,
+                               with_cosmology);
+        }
+      } break;
+
       case swift_type_stars: {
         if (Nstars == Nstars_written) {
 
@@ -759,6 +793,7 @@ void write_output_distributed(struct engine* e,
     if (gparts_written) swift_free("gparts_written", gparts_written);
     if (gpart_group_data_written)
       swift_free("gpart_group_written", gpart_group_data_written);
+    if (sinks_written) swift_free("sinks_written", sinks_written);
     if (sparts_written) swift_free("sparts_written", sparts_written);
     if (bparts_written) swift_free("bparts_written", bparts_written);
 
diff --git a/src/engine.c b/src/engine.c
index 98c3cb433c56d3b15ee3ece1dfb8a4ea9852a842..61346629eb39e2cb70ab738f2ba67fca7674e2ad 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -128,7 +128,8 @@ const char *engine_policy_names[] = {"none",
                                      "time-step limiter",
                                      "time-step sync",
                                      "logger",
-                                     "line of sight"};
+                                     "line of sight",
+                                     "sink"};
 
 /** The rank of the engine as a global variable (for messages). */
 int engine_rank;
@@ -1513,6 +1514,7 @@ void engine_print_task_counts(const struct engine *e) {
   fflush(stdout);
   message("nr_parts = %zu.", e->s->nr_parts);
   message("nr_gparts = %zu.", e->s->nr_gparts);
+  message("nr_sink = %zu.", e->s->nr_sinks);
   message("nr_sparts = %zu.", e->s->nr_sparts);
   message("nr_bparts = %zu.", e->s->nr_bparts);
 
@@ -1770,9 +1772,10 @@ void engine_rebuild(struct engine *e, const int repartitioned,
     engine_recompute_displacement_constraint(e);
 
 #ifdef SWIFT_DEBUG_CHECKS
-  part_verify_links(e->s->parts, e->s->gparts, e->s->sparts, e->s->bparts,
-                    e->s->nr_parts, e->s->nr_gparts, e->s->nr_sparts,
-                    e->s->nr_bparts, e->verbose);
+  part_verify_links(e->s->parts, e->s->gparts, e->s->sinks, e->s->sparts,
+                    e->s->bparts, e->s->nr_parts, e->s->nr_gparts,
+                    e->s->nr_sinks, e->s->nr_sparts, e->s->nr_bparts,
+                    e->verbose);
 #endif
 
   /* Initial cleaning up session ? */
@@ -2160,6 +2163,7 @@ void engine_first_init_particles(struct engine *e) {
   space_first_init_gparts(e->s, e->verbose);
   space_first_init_sparts(e->s, e->verbose);
   space_first_init_bparts(e->s, e->verbose);
+  space_first_init_sinks(e->s, e->verbose);
 
   if (e->verbose)
     message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
@@ -2213,6 +2217,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
   space_init_gparts(s, e->verbose);
   space_init_sparts(s, e->verbose);
   space_init_bparts(s, e->verbose);
+  space_init_sinks(s, e->verbose);
 
   /* Update the cooling function */
   if ((e->policy & engine_policy_cooling) ||
@@ -2284,6 +2289,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
   space_init_gparts(e->s, e->verbose);
   space_init_sparts(e->s, e->verbose);
   space_init_bparts(e->s, e->verbose);
+  space_init_sinks(e->s, e->verbose);
 
   /* Print the number of active tasks ? */
   if (e->verbose) engine_print_task_counts(e);
@@ -2419,9 +2425,10 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
 
 #ifdef SWIFT_DEBUG_CHECKS
   space_check_timesteps(e->s);
-  part_verify_links(e->s->parts, e->s->gparts, e->s->sparts, e->s->bparts,
-                    e->s->nr_parts, e->s->nr_gparts, e->s->nr_sparts,
-                    e->s->nr_bparts, e->verbose);
+  part_verify_links(e->s->parts, e->s->gparts, e->s->sinks, e->s->sparts,
+                    e->s->bparts, e->s->nr_parts, e->s->nr_gparts,
+                    e->s->nr_sinks, e->s->nr_sparts, e->s->nr_bparts,
+                    e->verbose);
 #endif
 
   /* Ready to go */
@@ -3483,14 +3490,15 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
 
   /* Re-link everything to the gparts. */
   if (s->nr_gparts > 0)
-    part_relink_all_parts_to_gparts(s->gparts, s->nr_gparts, s->parts,
+    part_relink_all_parts_to_gparts(s->gparts, s->nr_gparts, s->parts, s->sinks,
                                     s->sparts, s->bparts, &e->threadpool);
 
 #ifdef SWIFT_DEBUG_CHECKS
 
   /* Verify that the links are correct */
-  part_verify_links(s->parts, s->gparts, s->sparts, s->bparts, s->nr_parts,
-                    s->nr_gparts, s->nr_sparts, s->nr_bparts, e->verbose);
+  part_verify_links(s->parts, s->gparts, s->sinks, s->sparts, s->bparts,
+                    s->nr_parts, s->nr_gparts, s->nr_sinks, s->nr_sparts,
+                    s->nr_bparts, e->verbose);
 #endif
 
   if (e->verbose)
@@ -3816,6 +3824,7 @@ static void engine_dumper_init(struct engine *e) {
  * @param output_options Output options for snapshots.
  * @param Ngas total number of gas particles in the simulation.
  * @param Ngparts total number of gravity particles in the simulation.
+ * @param Nsinks total number of sink particles in the simulation.
  * @param Nstars total number of star particles in the simulation.
  * @param Nblackholes total number of black holes in the simulation.
  * @param Nbackground_gparts Total number of background DM particles.
@@ -3841,9 +3850,9 @@ static void engine_dumper_init(struct engine *e) {
  */
 void engine_init(struct engine *e, struct space *s, struct swift_params *params,
                  struct output_options *output_options, long long Ngas,
-                 long long Ngparts, long long Nstars, long long Nblackholes,
-                 long long Nbackground_gparts, int policy, int verbose,
-                 struct repartition *reparttype,
+                 long long Ngparts, long long Nsinks, long long Nstars,
+                 long long Nblackholes, long long Nbackground_gparts,
+                 int policy, int verbose, struct repartition *reparttype,
                  const struct unit_system *internal_units,
                  const struct phys_const *physical_constants,
                  struct cosmology *cosmo, struct hydro_props *hydro,
@@ -3868,6 +3877,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
   e->total_nr_parts = Ngas;
   e->total_nr_gparts = Ngparts;
   e->total_nr_sparts = Nstars;
+  e->total_nr_sinks = Nsinks;
   e->total_nr_bparts = Nblackholes;
   e->total_nr_DM_background_gparts = Nbackground_gparts;
   e->proxy_ind = NULL;
diff --git a/src/engine.h b/src/engine.h
index 63492619c6e06895bb1727d7cd4e95fb82e99548..4fb53017f14f9c1b8e64ddf5cd12587b376ea3c4 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -81,8 +81,9 @@ enum engine_policy {
   engine_policy_timestep_sync = (1 << 22),
   engine_policy_logger = (1 << 23),
   engine_policy_line_of_sight = (1 << 24),
+  engine_policy_sink = (1 << 25),
 };
-#define engine_maxpolicy 25
+#define engine_maxpolicy 26
 extern const char *engine_policy_names[engine_maxpolicy + 1];
 
 /**
@@ -505,6 +506,9 @@ struct engine {
   integertime_t ti_next_los;
   int los_output_count;
 
+  /* Total number of sink particles in the system. */
+  long long total_nr_sinks;
+
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   /* Run brute force checks only on steps when all gparts active? */
   int force_checks_only_all_active;
@@ -543,9 +547,9 @@ void engine_dump_snapshot(struct engine *e);
 void engine_init_output_lists(struct engine *e, struct swift_params *params);
 void engine_init(struct engine *e, struct space *s, struct swift_params *params,
                  struct output_options *output_options, long long Ngas,
-                 long long Ngparts, long long Nstars, long long Nblackholes,
-                 long long Nbackground_gparts, int policy, int verbose,
-                 struct repartition *reparttype,
+                 long long Ngparts, long long Nsinks, long long Nstars,
+                 long long Nblackholes, long long Nbackground_gparts,
+                 int policy, int verbose, struct repartition *reparttype,
                  const struct unit_system *internal_units,
                  const struct phys_const *physical_constants,
                  struct cosmology *cosmo, struct hydro_props *hydro,
diff --git a/src/engine_drift.c b/src/engine_drift.c
index 26f8160200c5bdc51f7dab51b56e9a4d771d3d39..0e447d6235b6e3f4b1f71dfea173daa8971d56f7 100644
--- a/src/engine_drift.c
+++ b/src/engine_drift.c
@@ -357,9 +357,10 @@ void engine_drift_all(struct engine *e, const int drift_mpoles) {
   space_check_drift_point(
       e->s, e->ti_current,
       drift_mpoles && (e->policy & engine_policy_self_gravity));
-  part_verify_links(e->s->parts, e->s->gparts, e->s->sparts, e->s->bparts,
-                    e->s->nr_parts, e->s->nr_gparts, e->s->nr_sparts,
-                    e->s->nr_bparts, e->verbose);
+  part_verify_links(e->s->parts, e->s->gparts, e->s->sinks, e->s->sparts,
+                    e->s->bparts, e->s->nr_parts, e->s->nr_gparts,
+                    e->s->nr_sinks, e->s->nr_sparts, e->s->nr_bparts,
+                    e->verbose);
 #endif
 
   if (e->verbose)
diff --git a/src/engine_redistribute.c b/src/engine_redistribute.c
index ea61533f0f7dc4014042d25841a8c6adabdae4d6..7246f71ee0a0d7620e6750c3706069720acce7ac 100644
--- a/src/engine_redistribute.c
+++ b/src/engine_redistribute.c
@@ -518,6 +518,12 @@ static void engine_redistribute_relink_mapper(void *map_data, int num_elements,
 void engine_redistribute(struct engine *e) {
 
 #ifdef WITH_MPI
+#ifdef SWIFT_DEBUG_CHECKS
+  const int nr_sinks_new = 0;
+#endif
+  if (e->policy & engine_policy_sink) {
+    error("Not implemented yet");
+  }
 
   const int nr_nodes = e->nr_nodes;
   const int nodeID = e->nodeID;
@@ -868,8 +874,8 @@ void engine_redistribute(struct engine *e) {
 
   /* Sort the gparticles according to their cell index. */
   if (nr_gparts > 0)
-    space_gparts_sort(s->gparts, s->parts, s->sparts, s->bparts, g_dest,
-                      &g_counts[nodeID * nr_nodes], nr_nodes);
+    space_gparts_sort(s->gparts, s->parts, s->sinks, s->sparts, s->bparts,
+                      g_dest, &g_counts[nodeID * nr_nodes], nr_nodes);
 
 #ifdef SWIFT_DEBUG_CHECKS
   /* Verify that the gpart have been sorted correctly. */
@@ -1194,8 +1200,9 @@ void engine_redistribute(struct engine *e) {
   }
 
   /* Verify that the links are correct */
-  part_verify_links(s->parts, s->gparts, s->sparts, s->bparts, nr_parts_new,
-                    nr_gparts_new, nr_sparts_new, nr_bparts_new, e->verbose);
+  part_verify_links(s->parts, s->gparts, s->sinks, s->sparts, s->bparts,
+                    nr_parts_new, nr_gparts_new, nr_sinks_new, nr_sparts_new,
+                    nr_bparts_new, e->verbose);
 
 #endif
 
diff --git a/src/engine_split_particles.c b/src/engine_split_particles.c
index b34a37ee039aac72e41c7ff0cae27e2bd1f1ea90..04e1ff90eb14dbedf509d05bc10a5049c2601a79 100644
--- a/src/engine_split_particles.c
+++ b/src/engine_split_particles.c
@@ -356,15 +356,17 @@ void engine_split_gas_particles(struct engine *e) {
 
     /* We now need to correct all the pointers of the other particle arrays */
     part_relink_all_parts_to_gparts(gparts_new, s->nr_gparts, s->parts,
-                                    s->sparts, s->bparts, &e->threadpool);
+                                    s->sinks, s->sparts, s->bparts,
+                                    &e->threadpool);
     s->gparts = gparts_new;
   }
 
 #ifdef SWIFT_DEBUG_CHECKS
   /* Verify that whatever reallocation happened we are still having correct
    * links */
-  part_verify_links(s->parts, s->gparts, s->sparts, s->bparts, s->nr_parts,
-                    s->nr_gparts, s->nr_sparts, s->nr_bparts, e->verbose);
+  part_verify_links(s->parts, s->gparts, s->sinks, s->sparts, s->bparts,
+                    s->nr_parts, s->nr_gparts, s->nr_sinks, s->nr_sparts,
+                    s->nr_bparts, e->verbose);
 #endif
 
   /* We now have enough memory in the part array to accomodate the new
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index cba60ce5b282fd5dc8852967d9bb83b7b747e9b5..11b5139d0a143b592e503e501aced3d2fe1b8362 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -24,6 +24,7 @@
 
 /* Local includes. */
 #include "cosmology.h"
+#include "error.h"
 #include "gravity_properties.h"
 #include "kernel_gravity.h"
 #include "minmax.h"
diff --git a/src/gravity/MultiSoftening/gravity.h b/src/gravity/MultiSoftening/gravity.h
index bbbea116e4c397a03236aa6d06db82f2a537b5d2..6ed0c38b7ee90807cc412efeb765167ded77e70c 100644
--- a/src/gravity/MultiSoftening/gravity.h
+++ b/src/gravity/MultiSoftening/gravity.h
@@ -48,7 +48,6 @@ __attribute__((always_inline)) INLINE static float gravity_get_mass(
  */
 __attribute__((always_inline)) INLINE static float gravity_get_softening(
     const struct gpart* gp, const struct gravity_props* restrict grav_props) {
-
   return gp->epsilon;
 }
 
diff --git a/src/io_properties.h b/src/io_properties.h
index dc84c0ed1aadd03e4e3b4f02602732062395af3f..0e5d92f4974aec776708b1b1688474de34afca3d 100644
--- a/src/io_properties.h
+++ b/src/io_properties.h
@@ -78,6 +78,14 @@ typedef void (*conversion_func_bpart_double)(const struct engine*,
 typedef void (*conversion_func_bpart_long_long)(const struct engine*,
                                                 const struct bpart*,
                                                 long long*);
+typedef void (*conversion_func_sink_float)(const struct engine*,
+                                           const struct sink*, float*);
+typedef void (*conversion_func_sink_int)(const struct engine*,
+                                         const struct sink*, int*);
+typedef void (*conversion_func_sink_double)(const struct engine*,
+                                            const struct sink*, double*);
+typedef void (*conversion_func_sink_long_long)(const struct engine*,
+                                               const struct sink*, long long*);
 
 /**
  * @brief The properties of a given dataset for i/o
@@ -127,6 +135,7 @@ struct io_props {
   const struct gpart* gparts;
   const struct spart* sparts;
   const struct bpart* bparts;
+  const struct sink* sinks;
 
   /* Are we converting? */
   int conversion;
@@ -154,6 +163,12 @@ struct io_props {
   conversion_func_bpart_int convert_bpart_i;
   conversion_func_bpart_double convert_bpart_d;
   conversion_func_bpart_long_long convert_bpart_l;
+
+  /* Conversion function for sink */
+  conversion_func_sink_float convert_sink_f;
+  conversion_func_sink_int convert_sink_i;
+  conversion_func_sink_double convert_sink_d;
+  conversion_func_sink_long_long convert_sink_l;
 };
 
 /**
@@ -206,6 +221,9 @@ INLINE static struct io_props io_make_input_field_(
   r.convert_bpart_f = NULL;
   r.convert_bpart_d = NULL;
   r.convert_bpart_l = NULL;
+  r.convert_sink_f = NULL;
+  r.convert_sink_d = NULL;
+  r.convert_sink_l = NULL;
 
   return r;
 }
@@ -994,4 +1012,185 @@ INLINE static struct io_props io_make_output_field_convert_bpart_LONGLONG(
   return r;
 }
 
+/**
+ * @brief Constructs an #io_props (with conversion) from its parameters
+ */
+#define io_make_output_field_convert_sink(name, type, dim, units, a_exponent,  \
+                                          sink, convert, desc)                 \
+  io_make_output_field_convert_sink_##type(name, type, dim, units, a_exponent, \
+                                           sizeof(sink[0]), sink, convert,     \
+                                           desc)
+
+/**
+ * @brief Construct an #io_props from its parameters
+ *
+ * @param name Name of the field to read
+ * @param type The type of the data
+ * @param dimension Dataset dimension (1D, 3D, ...)
+ * @param units The units of the dataset
+ * @param a_exponent Exponent of the scale-factor to convert to physical units.
+ * @param sinkSize The size in byte of the particle
+ * @param sinks The particle array
+ * @param functionPtr The function used to convert a sink-particle to a float
+ * @param description Description of the field added to the meta-data.
+ *
+ * Do not call this function directly. Use the macro defined above.
+ */
+INLINE static struct io_props io_make_output_field_convert_sink_INT(
+    const char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension,
+    enum unit_conversion_factor units, float a_exponent, size_t sinkSize,
+    const struct sink* sinks, conversion_func_sink_int functionPtr,
+    const char description[DESCRIPTION_BUFFER_SIZE]) {
+
+  struct io_props r;
+  bzero(&r, sizeof(struct io_props));
+
+  strcpy(r.name, name);
+  if (strlen(description) == 0) {
+    sprintf(r.description, "No description given");
+  } else {
+    strcpy(r.description, description);
+  }
+  r.type = type;
+  r.dimension = dimension;
+  r.importance = UNUSED;
+  r.units = units;
+  r.scale_factor_exponent = a_exponent;
+  r.partSize = sinkSize;
+  r.sinks = sinks;
+  r.conversion = 1;
+  r.convert_sink_i = functionPtr;
+
+  return r;
+}
+
+/**
+ * @brief Construct an #io_props from its parameters
+ *
+ * @param name Name of the field to read
+ * @param type The type of the data
+ * @param dimension Dataset dimension (1D, 3D, ...)
+ * @param units The units of the dataset
+ * @param a_exponent Exponent of the scale-factor to convert to physical units.
+ * @param sinkSize The size in byte of the particle
+ * @param sinks The particle array
+ * @param functionPtr The function used to convert a sink-particle to a float
+ * @param description Description of the field added to the meta-data.
+ *
+ * Do not call this function directly. Use the macro defined above.
+ */
+INLINE static struct io_props io_make_output_field_convert_sink_FLOAT(
+    const char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension,
+    enum unit_conversion_factor units, float a_exponent, size_t sinkSize,
+    const struct sink* sinks, conversion_func_sink_float functionPtr,
+    const char description[DESCRIPTION_BUFFER_SIZE]) {
+
+  struct io_props r;
+  bzero(&r, sizeof(struct io_props));
+
+  strcpy(r.name, name);
+  if (strlen(description) == 0) {
+    sprintf(r.description, "No description given");
+  } else {
+    strcpy(r.description, description);
+  }
+  r.type = type;
+  r.dimension = dimension;
+  r.importance = UNUSED;
+  r.units = units;
+  r.scale_factor_exponent = a_exponent;
+  r.partSize = sinkSize;
+  r.sinks = sinks;
+  r.conversion = 1;
+  r.convert_sink_f = functionPtr;
+
+  return r;
+}
+
+/**
+ * @brief Construct an #io_props from its parameters
+ *
+ * @param name Name of the field to read
+ * @param type The type of the data
+ * @param dimension Dataset dimension (1D, 3D, ...)
+ * @param units The units of the dataset
+ * @param a_exponent Exponent of the scale-factor to convert to physical units.
+ * @param sinkSize The size in byte of the particle
+ * @param sinks The particle array
+ * @param functionPtr The function used to convert a sink-particle to a double
+ * @param description Description of the field added to the meta-data.
+ *
+ * Do not call this function directly. Use the macro defined above.
+ */
+INLINE static struct io_props io_make_output_field_convert_sink_DOUBLE(
+    const char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension,
+    enum unit_conversion_factor units, float a_exponent, size_t sinkSize,
+    const struct sink* sinks, conversion_func_sink_double functionPtr,
+    const char description[DESCRIPTION_BUFFER_SIZE]) {
+
+  struct io_props r;
+  bzero(&r, sizeof(struct io_props));
+
+  strcpy(r.name, name);
+  if (strlen(description) == 0) {
+    sprintf(r.description, "No description given");
+  } else {
+    strcpy(r.description, description);
+  }
+  r.type = type;
+  r.dimension = dimension;
+  r.importance = UNUSED;
+  r.units = units;
+  r.scale_factor_exponent = a_exponent;
+  r.partSize = sinkSize;
+  r.sinks = sinks;
+  r.conversion = 1;
+  r.convert_sink_d = functionPtr;
+
+  return r;
+}
+
+/**
+ * @brief Construct an #io_props from its parameters
+ *
+ * @param name Name of the field to read
+ * @param type The type of the data
+ * @param dimension Dataset dimension (1D, 3D, ...)
+ * @param units The units of the dataset
+ * @param a_exponent Exponent of the scale-factor to convert to physical units.
+ * @param sinkSize The size in byte of the particle
+ * @param sinks The particle array
+ * @param functionPtr The function used to convert a sink-particle to a double
+ * @param description Description of the field added to the meta-data.
+ *
+ * Do not call this function directly. Use the macro defined above.
+ */
+INLINE static struct io_props io_make_output_field_convert_sink_LONGLONG(
+    const char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension,
+    enum unit_conversion_factor units, float a_exponent, size_t sinkSize,
+    const struct sink* sinks, conversion_func_sink_long_long functionPtr,
+    const char description[DESCRIPTION_BUFFER_SIZE]) {
+
+  struct io_props r;
+  bzero(&r, sizeof(struct io_props));
+
+  strcpy(r.name, name);
+  if (strlen(description) == 0) {
+    sprintf(r.description, "No description given");
+  } else {
+    strcpy(r.description, description);
+  }
+  r.type = type;
+  r.dimension = dimension;
+  r.importance = UNUSED;
+  r.units = units;
+  r.scale_factor_exponent = a_exponent;
+  r.partSize = sinkSize;
+  r.sinks = sinks;
+  r.conversion = 1;
+  r.convert_sink_l = functionPtr;
+
+  return r;
+}
+
 #endif /* SWIFT_IO_PROPERTIES_H */
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 7b19ed7e18e7163340f305a97abef7e0db463c3b..6a1d004e493889b300930637876d402166da48e8 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -55,6 +55,7 @@
 #include "output_options.h"
 #include "part.h"
 #include "part_type.h"
+#include "sink_io.h"
 #include "star_formation_io.h"
 #include "stars_io.h"
 #include "tracers_io.h"
@@ -691,18 +692,21 @@ void write_array_parallel(struct engine* e, hid_t grp, char* fileName,
  * @param dim (output) The dimension of the volume read from the file.
  * @param parts (output) The array of #part read from the file.
  * @param gparts (output) The array of #gpart read from the file.
+ * @param sinks (output) The array of #sink read from the file.
  * @param sparts (output) The array of #spart read from the file.
  * @param bparts (output) The array of #bpart read from the file.
  * @param Ngas (output) The number of particles read from the file.
  * @param Ngparts (output) The number of particles read from the file.
  * @param Ngparts_background (output) The number of background DM particles read
  * from the file.
+ * @param Nsink (output) The number of particles read from the file.
  * @param Nstars (output) The number of particles read from the file.
  * @param Nblackholes (output) The number of particles read from the file.
  * @param flag_entropy (output) 1 if the ICs contained Entropy in the
  * InternalEnergy field
  * @param with_hydro Are we running with hydro ?
  * @param with_gravity Are we running with gravity ?
+ * @param with_sink Are we running with sink ?
  * @param with_stars Are we running with stars ?
  * @param with_black_holes Are we running with black holes ?
  * @param with_cosmology Are we running with cosmology ?
@@ -721,14 +725,15 @@ void write_array_parallel(struct engine* e, hid_t grp, char* fileName,
  */
 void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
                       double dim[3], struct part** parts, struct gpart** gparts,
-                      struct spart** sparts, struct bpart** bparts,
-                      size_t* Ngas, size_t* Ngparts, size_t* Ngparts_background,
+                      struct sink** sinks, struct spart** sparts,
+                      struct bpart** bparts, size_t* Ngas, size_t* Ngparts,
+                      size_t* Ngparts_background, size_t* Nsinks,
                       size_t* Nstars, size_t* Nblackholes, int* flag_entropy,
-                      int with_hydro, int with_gravity, int with_stars,
-                      int with_black_holes, int with_cosmology, int cleanup_h,
-                      int cleanup_sqrt_a, double h, double a, int mpi_rank,
-                      int mpi_size, MPI_Comm comm, MPI_Info info, int n_threads,
-                      int dry_run) {
+                      int with_hydro, int with_gravity, int with_sink,
+                      int with_stars, int with_black_holes, int with_cosmology,
+                      int cleanup_h, int cleanup_sqrt_a, double h, double a,
+                      int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
+                      int n_threads, int dry_run) {
 
   hid_t h_file = 0, h_grp = 0;
   /* GADGET has only cubic boxes (in cosmological mode) */
@@ -744,7 +749,7 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
 
   /* Initialise counters */
   *Ngas = 0, *Ngparts = 0, *Ngparts_background = 0, *Nstars = 0,
-  *Nblackholes = 0;
+  *Nblackholes = 0, *Nsinks = 0;
 
   /* Open file */
   /* message("Opening file '%s' as IC.", fileName); */
@@ -883,6 +888,15 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
     bzero(*parts, *Ngas * sizeof(struct part));
   }
 
+  /* Allocate memory to store black hole particles */
+  if (with_sink) {
+    *Nsinks = N[swift_type_sink];
+    if (swift_memalign("sinks", (void**)sinks, sink_align,
+                       *Nsinks * sizeof(struct sink)) != 0)
+      error("Error while allocating memory for sink particles");
+    bzero(*sinks, *Nsinks * sizeof(struct sink));
+  }
+
   /* Allocate memory to store stars particles */
   if (with_stars) {
     *Nstars = N[swift_type_stars];
@@ -895,7 +909,7 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
   /* Allocate memory to store black hole particles */
   if (with_black_holes) {
     *Nblackholes = N[swift_type_black_hole];
-    if (swift_memalign("sparts", (void**)bparts, bpart_align,
+    if (swift_memalign("bparts", (void**)bparts, bpart_align,
                        *Nblackholes * sizeof(struct bpart)) != 0)
       error("Error while allocating memory for black_holes particles");
     bzero(*bparts, *Nblackholes * sizeof(struct bpart));
@@ -909,6 +923,7 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
                N[swift_type_dark_matter] +
                N[swift_type_dark_matter_background] +
                (with_stars ? N[swift_type_stars] : 0) +
+               (with_sink ? N[swift_type_sink] : 0) +
                (with_black_holes ? N[swift_type_black_hole] : 0);
     *Ngparts_background = Ndm_background;
     if (swift_memalign("gparts", (void**)gparts, gpart_align,
@@ -966,6 +981,13 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
         }
         break;
 
+      case swift_type_sink:
+        if (with_sink) {
+          Nparticles = *Nsinks;
+          sink_read_particles(*sinks, list, &num_fields);
+        }
+        break;
+
       case swift_type_stars:
         if (with_stars) {
           Nparticles = *Nstars;
@@ -1014,15 +1036,21 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
       io_duplicate_hydro_gparts(&tp, *parts, *gparts, *Ngas,
                                 Ndm + Ndm_background);
 
+    /* Duplicate the sink particles into gparts */
+    if (with_sink)
+      io_duplicate_sinks_gparts(&tp, *sinks, *gparts, *Nsinks,
+                                Ndm + Ndm_background + *Ngas);
+
     /* Duplicate the stars particles into gparts */
     if (with_stars)
       io_duplicate_stars_gparts(&tp, *sparts, *gparts, *Nstars,
-                                Ndm + Ndm_background + *Ngas);
+                                Ndm + Ndm_background + *Ngas + *Nsinks);
 
     /* Duplicate the stars particles into gparts */
     if (with_black_holes)
-      io_duplicate_black_holes_gparts(&tp, *bparts, *gparts, *Nblackholes,
-                                      Ndm + Ndm_background + *Ngas + *Nstars);
+      io_duplicate_black_holes_gparts(
+          &tp, *bparts, *gparts, *Nblackholes,
+          Ndm + Ndm_background + *Ngas + *Nsinks + *Nstars);
 
     threadpool_clean(&tp);
   }
@@ -1060,6 +1088,8 @@ void prepare_file(struct engine* e, const char* fileName,
   const struct gpart* gparts = e->s->gparts;
   const struct spart* sparts = e->s->sparts;
   const struct bpart* bparts = e->s->bparts;
+  const struct sink* sinks = e->s->sinks;
+
   struct output_options* output_options = e->output_options;
   const int with_cosmology = e->policy & engine_policy_cosmology;
   const int with_cooling = e->policy & engine_policy_cooling;
@@ -1233,6 +1263,10 @@ void prepare_file(struct engine* e, const char* fileName,
         }
         break;
 
+      case swift_type_sink:
+        sink_write_particles(sinks, list, &num_fields, with_cosmology);
+        break;
+
       case swift_type_stars:
         stars_write_particles(sparts, list, &num_fields, with_cosmology);
         num_fields += chemistry_write_sparticles(sparts, list + num_fields);
@@ -1334,6 +1368,7 @@ void write_output_parallel(struct engine* e,
   const struct gpart* gparts = e->s->gparts;
   const struct spart* sparts = e->s->sparts;
   const struct bpart* bparts = e->s->bparts;
+  const struct sink* sinks = e->s->sinks;
   struct output_options* output_options = e->output_options;
   struct output_list* output_list = e->output_list_snapshots;
   const int with_cosmology = e->policy & engine_policy_cosmology;
@@ -1352,6 +1387,7 @@ void write_output_parallel(struct engine* e,
   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 Nsinks = e->s->nr_sinks;
   const size_t Nblackholes = e->s->nr_bparts;
   // const size_t Nbaryons = Ngas + Nstars;
   // const size_t Ndm = Ntot > 0 ? Ntot - Nbaryons : 0;
@@ -1366,18 +1402,20 @@ void write_output_parallel(struct engine* e,
       e->s->nr_gparts - e->s->nr_inhibited_gparts - e->s->nr_extra_gparts;
   const size_t Ngas_written =
       e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
+  const size_t Nsinks_written =
+      e->s->nr_sinks - e->s->nr_inhibited_sinks - e->s->nr_extra_sinks;
   const size_t Nstars_written =
       e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
   const size_t Nblackholes_written =
       e->s->nr_bparts - e->s->nr_inhibited_bparts - e->s->nr_extra_bparts;
   const size_t Nbaryons_written =
-      Ngas_written + Nstars_written + Nblackholes_written;
+      Ngas_written + Nstars_written + Nblackholes_written + Nsinks_written;
   const size_t Ndm_written =
       Ntot_written > 0 ? Ntot_written - Nbaryons_written - Ndm_background : 0;
 
   /* Compute offset in the file and total number of particles */
   size_t N[swift_type_count] = {Ngas_written,   Ndm_written,
-                                Ndm_background, 0,
+                                Ndm_background, Nsinks_written,
                                 Nstars_written, Nblackholes_written};
   long long N_total[swift_type_count] = {0};
   long long offset[swift_type_count] = {0};
@@ -1544,6 +1582,7 @@ void write_output_parallel(struct engine* e,
     struct velociraptor_gpart_data* gpart_group_data_written = NULL;
     struct spart* sparts_written = NULL;
     struct bpart* bparts_written = NULL;
+    struct sink* sinks_written = NULL;
 
     /* Write particle fields from the particle structure */
     switch (ptype) {
@@ -1703,6 +1742,33 @@ void write_output_parallel(struct engine* e,
 
       } break;
 
+      case swift_type_sink: {
+        if (Nsinks == Nsinks_written) {
+
+          /* No inhibted particles: easy case */
+          Nparticles = Nsinks;
+          sink_write_particles(sinks, list, &num_fields, with_cosmology);
+        } else {
+
+          /* Ok, we need to fish out the particles we want */
+          Nparticles = Nsinks_written;
+
+          /* Allocate temporary arrays */
+          if (swift_memalign("sinks_written", (void**)&sinks_written,
+                             sink_align,
+                             Nsinks_written * sizeof(struct sink)) != 0)
+            error("Error while allocating temporary memory for sink");
+
+          /* Collect the particles we want to write */
+          io_collect_sinks_to_write(sinks, sinks_written, Nsinks,
+                                    Nsinks_written);
+
+          /* Select the fields to write */
+          sink_write_particles(sinks_written, list, &num_fields,
+                               with_cosmology);
+        }
+      } break;
+
       case swift_type_stars: {
         if (Nstars == Nstars_written) {
 
@@ -1830,6 +1896,7 @@ void write_output_parallel(struct engine* e,
       swift_free("gpart_group_written", gpart_group_data_written);
     if (sparts_written) swift_free("sparts_written", sparts_written);
     if (bparts_written) swift_free("bparts_written", bparts_written);
+    if (sinks_written) swift_free("sinks_written", sinks_written);
 
 #ifdef IO_SPEED_MEASUREMENT
     MPI_Barrier(MPI_COMM_WORLD);
diff --git a/src/parallel_io.h b/src/parallel_io.h
index b816d1d7c6fa74cf1ec83bba2bbf768639811f17..9b5795fe2840bb97545452d106de888d1de17025 100644
--- a/src/parallel_io.h
+++ b/src/parallel_io.h
@@ -35,13 +35,14 @@ struct unit_system;
 
 void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
                       double dim[3], struct part** parts, struct gpart** gparts,
-                      struct spart** sparts, struct bpart** bparts,
-                      size_t* Ngas, size_t* Ngparts, size_t* Ngparts_background,
+                      struct sink** sinks, struct spart** sparts,
+                      struct bpart** bparts, size_t* Ngas, size_t* Ngparts,
+                      size_t* Ngparts_background, size_t* Nsinks,
                       size_t* Nsparts, size_t* Nbparts, int* flag_entropy,
-                      int with_hydro, int with_gravity, int with_stars,
-                      int with_black_holes, int with_cosmology, int cleanup_h,
-                      int cleanup_sqrt_a, double h, double a, int mpi_rank,
-                      int mpi_size, MPI_Comm comm, MPI_Info info,
+                      int with_hydro, int with_gravity, int with_sinks,
+                      int with_stars, int with_black_holes, int with_cosmology,
+                      int cleanup_h, int cleanup_sqrt_a, double h, double a,
+                      int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
                       int nr_threads, int dry_run);
 
 void write_output_parallel(struct engine* e,
diff --git a/src/part.c b/src/part.c
index be685f1d862e5fe69fd45c73998fbfc0a0b3b2ff..cfc5338a69ee0f85483b8b26591f3f9d0c7d6475 100644
--- a/src/part.c
+++ b/src/part.c
@@ -65,6 +65,22 @@ void part_relink_gparts_to_sparts(struct spart *sparts, const size_t N,
   }
 }
 
+/**
+ * @brief Re-link the #gpart%s associated with the list of #sink%s.
+ *
+ * @param sinks The list of #sink.
+ * @param N The number of sink-particles to re-link;
+ * @param offset The offset of #sink%s relative to the global sinks list.
+ */
+void part_relink_gparts_to_sinks(struct sink *sinks, const size_t N,
+                                 const ptrdiff_t offset) {
+  for (size_t k = 0; k < N; k++) {
+    if (sinks[k].gpart) {
+      sinks[k].gpart->id_or_neg_offset = -(k + offset);
+    }
+  }
+}
+
 /**
  * @brief Re-link the #gpart%s associated with the list of #bpart%s.
  *
@@ -129,12 +145,29 @@ void part_relink_bparts_to_gparts(struct gpart *gparts, const size_t N,
   }
 }
 
+/**
+ * @brief Re-link the #sink%s associated with the list of #gpart%s.
+ *
+ * @param gparts The list of #gpart.
+ * @param N The number of particles to re-link;
+ * @param sinks The global #sink array in which to find the #gpart offsets.
+ */
+void part_relink_sinks_to_gparts(struct gpart *gparts, const size_t N,
+                                 struct sink *sinks) {
+  for (size_t k = 0; k < N; k++) {
+    if (gparts[k].type == swift_type_sink) {
+      sinks[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
+    }
+  }
+}
+
 /**
  * @brief Helper structure to pass data to the liking mapper functions.
  */
 struct relink_data {
   struct part *const parts;
   struct gpart *const garts;
+  struct sink *const sinks;
   struct spart *const sparts;
   struct bpart *const bparts;
 };
@@ -156,6 +189,7 @@ void part_relink_all_parts_to_gparts_mapper(void *restrict map_data, int count,
   struct spart *const sparts = data->sparts;
   struct bpart *const bparts = data->bparts;
   struct gpart *const gparts = (struct gpart *)map_data;
+  struct sink *const sinks = data->sinks;
 
   for (int k = 0; k < count; k++) {
     if (gparts[k].type == swift_type_gas) {
@@ -164,12 +198,14 @@ void part_relink_all_parts_to_gparts_mapper(void *restrict map_data, int count,
       sparts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
     } else if (gparts[k].type == swift_type_black_hole) {
       bparts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
+    } else if (gparts[k].type == swift_type_sink) {
+      sinks[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
     }
   }
 }
 
 /**
- * @brief Re-link both the #part%s, #spart%s and #bpart%s associated
+ * @brief Re-link both the #part%s, #sink%s, #spart%s and #bpart%s associated
  * with the list of #gpart%s.
  *
  * This function uses thread parallelism and should not be called inside
@@ -179,41 +215,44 @@ void part_relink_all_parts_to_gparts_mapper(void *restrict map_data, int count,
  * @param gparts The list of #gpart.
  * @param N The number of particles to re-link;
  * @param parts The global #part array in which to find the #gpart offsets.
+ * @param sinks The global #sink array in which to find the #gpart offsets.
  * @param sparts The global #spart array in which to find the #gpart offsets.
  * @param bparts The global #bpart array in which to find the #gpart offsets.
  * @param tp The #threadpool object.
  */
 void part_relink_all_parts_to_gparts(struct gpart *gparts, const size_t N,
-                                     struct part *parts, struct spart *sparts,
-                                     struct bpart *bparts,
+                                     struct part *parts, struct sink *sinks,
+                                     struct spart *sparts, struct bpart *bparts,
                                      struct threadpool *tp) {
 
-  struct relink_data data = {parts, /*gparts=*/NULL, sparts, bparts};
+  struct relink_data data = {parts, /*gparts=*/NULL, sinks, sparts, bparts};
   threadpool_map(tp, part_relink_all_parts_to_gparts_mapper, gparts, N,
                  sizeof(struct gpart), 0, &data);
 }
 
 /**
- * @brief Verifies that the #gpart, #part and #spart are correctly linked
- * together
- * and that the particle poisitions match.
+ * @brief Verifies that the #gpart, #part, #sink, #spart and #bpart are
+ * correctly linked together and that the particle positions match.
  *
  * This is a debugging function.
  *
  * @param parts The #part array.
  * @param gparts The #gpart array.
+ * @param sinks The #sink array.
  * @param sparts The #spart array.
  * @param bparts The #bpart array.
  * @param nr_parts The number of #part in the array.
  * @param nr_gparts The number of #gpart in the array.
+ * @param nr_sinks The number of #sink in the array.
  * @param nr_sparts The number of #spart in the array.
  * @param nr_bparts The number of #bpart in the array.
  * @param verbose Do we report verbosely in case of success ?
  */
 void part_verify_links(struct part *parts, struct gpart *gparts,
-                       struct spart *sparts, struct bpart *bparts,
-                       size_t nr_parts, size_t nr_gparts, size_t nr_sparts,
-                       size_t nr_bparts, int verbose) {
+                       struct sink *sinks, struct spart *sparts,
+                       struct bpart *bparts, size_t nr_parts, size_t nr_gparts,
+                       size_t nr_sinks, size_t nr_sparts, size_t nr_bparts,
+                       int verbose) {
 
   ticks tic = getticks();
 
@@ -339,6 +378,40 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
       if (gparts[k].time_bin != bpart->time_bin)
         error("Linked particles are not at the same time !");
     }
+
+    else if (gparts[k].type == swift_type_sink) {
+
+      /* Check that it is linked */
+      if (gparts[k].id_or_neg_offset > 0)
+        error("Sink gpart not linked to anything !");
+
+      /* Find its link */
+      const struct sink *sink = &sinks[-gparts[k].id_or_neg_offset];
+
+      /* Check the reverse link */
+      if (sink->gpart != &gparts[k]) error("Linking problem !");
+
+      /* Check that the particles are at the same place */
+      if (gparts[k].x[0] != sink->x[0] || gparts[k].x[1] != sink->x[1] ||
+          gparts[k].x[2] != sink->x[2])
+        error(
+            "Linked particles are not at the same position !\n"
+            "gp->x=[%e %e %e] sink->x=[%e %e %e] diff=[%e %e %e]",
+            gparts[k].x[0], gparts[k].x[1], gparts[k].x[2], sink->x[0],
+            sink->x[1], sink->x[2], gparts[k].x[0] - sink->x[0],
+            gparts[k].x[1] - sink->x[1], gparts[k].x[2] - sink->x[2]);
+
+      /* Check that the particles have the same mass */
+      if (gparts[k].mass != sink->mass)
+        error(
+            "Linked particles do not have the same mass!\n"
+            "gp->m=%e sink->m=%e",
+            gparts[k].mass, sink->mass);
+
+      /* Check that the particles are at the same time */
+      if (gparts[k].time_bin != sink->time_bin)
+        error("Linked particles are not at the same time !");
+    }
   }
 
   /* Now check that all parts are linked */
@@ -422,6 +495,33 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
     }
   }
 
+  /* Now check that all sinks are linked */
+  for (size_t k = 0; k < nr_sinks; ++k) {
+
+    /* Ok, there is a link */
+    if (sinks[k].gpart != NULL) {
+
+      /* Check the link */
+      if (sinks[k].gpart->id_or_neg_offset != -(ptrdiff_t)k) {
+        error("Linking problem !");
+
+        /* Check that the particles are at the same place */
+        if (sinks[k].x[0] != sinks[k].gpart->x[0] ||
+            sinks[k].x[1] != sinks[k].gpart->x[1] ||
+            sinks[k].x[2] != sinks[k].gpart->x[2])
+          error("Linked particles are not at the same position !");
+
+        /* Check that the particles have the same mass */
+        if (sinks[k].mass != sinks[k].gpart->mass)
+          error("Linked particles do not have the same mass!\n");
+
+        /* Check that the particles are at the same time */
+        if (sinks[k].time_bin != sinks[k].gpart->time_bin)
+          error("Linked particles are not at the same time !");
+      }
+    }
+  }
+
   if (verbose) message("All links OK");
   if (verbose)
     message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
diff --git a/src/part.h b/src/part.h
index b01618dae8f659af93eec6c71bf2a7305e945f91..136ff23f6e1f820204541b86a7612ed5d1defdee 100644
--- a/src/part.h
+++ b/src/part.h
@@ -43,6 +43,7 @@ struct threadpool;
 #define spart_align 128
 #define gpart_align 128
 #define bpart_align 128
+#define sink_align 128
 
 /* Import the right hydro particle definition */
 #if defined(MINIMAL_SPH)
@@ -120,26 +121,38 @@ struct threadpool;
 #error "Invalid choice of black hole particle"
 #endif
 
+/* Import the right sink particle definition */
+#if defined(SINK_NONE)
+#include "./sink/Default/sink_part.h"
+#else
+#error "Invalid choice of sink particle"
+#endif
+
 void part_relink_gparts_to_parts(struct part *parts, const size_t N,
                                  const ptrdiff_t offset);
 void part_relink_gparts_to_sparts(struct spart *sparts, const size_t N,
                                   const ptrdiff_t offset);
 void part_relink_gparts_to_bparts(struct bpart *bparts, const size_t N,
                                   const ptrdiff_t offset);
+void part_relink_gparts_to_sinks(struct sink *sinks, const size_t N,
+                                 const ptrdiff_t offset);
 void part_relink_parts_to_gparts(struct gpart *gparts, const size_t N,
                                  struct part *parts);
 void part_relink_sparts_to_gparts(struct gpart *gparts, const size_t N,
                                   struct spart *sparts);
 void part_relink_bparts_to_gparts(struct gpart *gparts, const size_t N,
                                   struct bpart *bparts);
+void part_relink_sinks_to_gparts(struct gpart *gparts, const size_t N,
+                                 struct sink *sinks);
 void part_relink_all_parts_to_gparts(struct gpart *gparts, const size_t N,
-                                     struct part *parts, struct spart *sparts,
-                                     struct bpart *bparts,
+                                     struct part *parts, struct sink *sinks,
+                                     struct spart *sparts, struct bpart *bparts,
                                      struct threadpool *tp);
 void part_verify_links(struct part *parts, struct gpart *gparts,
-                       struct spart *sparts, struct bpart *bparts,
-                       size_t nr_parts, size_t nr_gparts, size_t nr_sparts,
-                       size_t nr_bparts, int verbose);
+                       struct sink *sinks, struct spart *sparts,
+                       struct bpart *bparts, size_t nr_parts, size_t nr_gparts,
+                       size_t nr_sinks, size_t nr_sparts, size_t nr_bparts,
+                       int verbose);
 
 #ifdef WITH_MPI
 /* MPI data type for the particle transfers */
diff --git a/src/part_type.c b/src/part_type.c
index 9cc8c35e7413e28d518289ccfcdd00d302ce8a28..e200a12bcf8ed72e9e1f404d84e5140190d026dc 100644
--- a/src/part_type.c
+++ b/src/part_type.c
@@ -21,4 +21,4 @@
 #include "part_type.h"
 
 const char* part_type_names[swift_type_count] = {
-    "Gas", "DM", "DMBackground", "Dummy", "Stars", "BH"};
+    "Gas", "DM", "DMBackground", "Sink", "Stars", "BH"};
diff --git a/src/part_type.h b/src/part_type.h
index 8311f5d17a7b838b9577a029387a5f01fb1a0801..ad93817304ad8ee1d6106385ed8bb77b3a564150 100644
--- a/src/part_type.h
+++ b/src/part_type.h
@@ -28,6 +28,7 @@ enum part_type {
   swift_type_gas = 0,
   swift_type_dark_matter = 1,
   swift_type_dark_matter_background = 2,
+  swift_type_sink = 3,
   swift_type_stars = 4,
   swift_type_black_hole = 5,
   swift_type_count
diff --git a/src/runner_time_integration.c b/src/runner_time_integration.c
index d966800fc989e6553de3f0fb8b7c4760a1c98bb8..0e14524aedda6276bd3e280e2c50040e9210a7fb 100644
--- a/src/runner_time_integration.c
+++ b/src/runner_time_integration.c
@@ -183,7 +183,10 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
 
       /* If the g-particle has no counterpart and needs to be kicked */
       if ((gp->type == swift_type_dark_matter ||
-           gp->type == swift_type_dark_matter_background) &&
+           gp->type == swift_type_dark_matter_background ||
+           gp->type == swift_type_sink) &&
+          // TODO loic remove this
+
           gpart_is_starting(gp, e)) {
 
         const integertime_t ti_step = get_integer_timestep(gp->time_bin);
@@ -407,7 +410,10 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
 
       /* If the g-particle has no counterpart and needs to be kicked */
       if ((gp->type == swift_type_dark_matter ||
-           gp->type == swift_type_dark_matter_background) &&
+           gp->type == swift_type_dark_matter_background ||
+           gp->type == swift_type_sink) &&
+          // TODO loic remove this
+
           gpart_is_active(gp, e)) {
 
         const integertime_t ti_step = get_integer_timestep(gp->time_bin);
@@ -672,7 +678,9 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
 
       /* If the g-particle has no counterpart */
       if (gp->type == swift_type_dark_matter ||
-          gp->type == swift_type_dark_matter_background) {
+          gp->type == swift_type_dark_matter_background ||
+          gp->type == swift_type_sink) {
+        // Loic TODO remove sink
 
         /* need to be updated ? */
         if (gpart_is_active(gp, e)) {
@@ -704,6 +712,13 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
           /* What is the next starting point for this cell ? */
           ti_gravity_beg_max = max(ti_current, ti_gravity_beg_max);
 
+          // Loic TODO remove this
+          /* Synchronize the time step */
+          if (gp->type == swift_type_sink) {
+            struct sink *sink = &e->s->sinks[-gp->id_or_neg_offset];
+            sink->time_bin = gp->time_bin;
+          }
+
         } else { /* gpart is inactive */
 
           if (!gpart_is_inhibited(gp, e)) {
diff --git a/src/serial_io.c b/src/serial_io.c
index 916d138e61193abe092d42a79fafacac63c0c4f6..16ed7e78014eedec2af41b3bcedf1a103ebc7e7e 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -55,6 +55,7 @@
 #include "output_options.h"
 #include "part.h"
 #include "part_type.h"
+#include "sink_io.h"
 #include "star_formation_io.h"
 #include "stars_io.h"
 #include "tracers_io.h"
@@ -460,12 +461,14 @@ void write_array_serial(const struct engine* e, hid_t grp, char* fileName,
  * @param dim (output) The dimension of the volume read from the file.
  * @param parts (output) The array of #part (gas particles) read from the file.
  * @param gparts (output) The array of #gpart read from the file.
+ * @param sinks (output) Array of #sink particles.
  * @param sparts (output) Array of #spart particles.
  * @param bparts (output) Array of #bpart particles.
  * @param Ngas (output) The number of #part read from the file on that node.
  * @param Ngparts (output) The number of #gpart read from the file on that node.
  * @param Ngparts_background (output) The number of background #gpart (type 2)
  * read from the file on that node.
+ * @param Nsinks (output) The number of #sink read from the file on that node.
  * @param Nstars (output) The number of #spart read from the file on that node.
  * @param Nblackholes (output) The number of #bpart read from the file on that
  * node.
@@ -473,6 +476,7 @@ void write_array_serial(const struct engine* e, hid_t grp, char* fileName,
  * InternalEnergy field
  * @param with_hydro Are we reading gas particles ?
  * @param with_gravity Are we reading/creating #gpart arrays ?
+ * @param with_sink Are we reading sink particles ?
  * @param with_stars Are we reading star particles ?
  * @param with_black_holes Are we reading black hole particles ?
  * @param with_cosmology Are we running with cosmology ?
@@ -498,13 +502,15 @@ void write_array_serial(const struct engine* e, hid_t grp, char* fileName,
  */
 void read_ic_serial(char* fileName, const struct unit_system* internal_units,
                     double dim[3], struct part** parts, struct gpart** gparts,
-                    struct spart** sparts, struct bpart** bparts, size_t* Ngas,
-                    size_t* Ngparts, size_t* Ngparts_background, size_t* Nstars,
+                    struct sink** sinks, struct spart** sparts,
+                    struct bpart** bparts, size_t* Ngas, size_t* Ngparts,
+                    size_t* Ngparts_background, size_t* Nsinks, size_t* Nstars,
                     size_t* Nblackholes, int* flag_entropy, int with_hydro,
-                    int with_gravity, int with_stars, int with_black_holes,
-                    int with_cosmology, int cleanup_h, int cleanup_sqrt_a,
-                    double h, double a, int mpi_rank, int mpi_size,
-                    MPI_Comm comm, MPI_Info info, int n_threads, int dry_run) {
+                    int with_gravity, int with_sink, int with_stars,
+                    int with_black_holes, int with_cosmology, int cleanup_h,
+                    int cleanup_sqrt_a, double h, double a, int mpi_rank,
+                    int mpi_size, MPI_Comm comm, MPI_Info info, int n_threads,
+                    int dry_run) {
 
   hid_t h_file = 0, h_grp = 0;
   /* GADGET has only cubic boxes (in cosmological mode) */
@@ -523,7 +529,7 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
 
   /* Initialise counters */
   *Ngas = 0, *Ngparts = 0, *Ngparts_background = 0, *Nstars = 0,
-  *Nblackholes = 0;
+  *Nblackholes = 0, *Nsinks = 0;
 
   /* First read some information about the content */
   if (mpi_rank == 0) {
@@ -671,6 +677,15 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
     bzero(*parts, *Ngas * sizeof(struct part));
   }
 
+  /* Allocate memory to store sinks particles */
+  if (with_sink) {
+    *Nsinks = N[swift_type_sink];
+    if (swift_memalign("sinks", (void**)sinks, sink_align,
+                       *Nsinks * sizeof(struct sink)) != 0)
+      error("Error while allocating memory for sink particles");
+    bzero(*sinks, *Nsinks * sizeof(struct sink));
+  }
+
   /* Allocate memory to store stars particles */
   if (with_stars) {
     *Nstars = N[swift_type_stars];
@@ -696,6 +711,7 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
     *Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
                N[swift_type_dark_matter] +
                N[swift_type_dark_matter_background] +
+               (with_sink ? N[swift_type_sink] : 0) +
                (with_stars ? N[swift_type_stars] : 0) +
                (with_black_holes ? N[swift_type_black_hole] : 0);
     *Ngparts_background = Ndm_background;
@@ -766,6 +782,13 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
             }
             break;
 
+          case swift_type_sink:
+            if (with_sink) {
+              Nparticles = *Nsinks;
+              sink_read_particles(*sinks, list, &num_fields);
+            }
+            break;
+
           case swift_type_stars:
             if (with_stars) {
               Nparticles = *Nstars;
@@ -823,15 +846,21 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
       io_duplicate_hydro_gparts(&tp, *parts, *gparts, *Ngas,
                                 Ndm + Ndm_background);
 
+    /* Duplicate the sink particles into gparts */
+    if (with_sink)
+      io_duplicate_sinks_gparts(&tp, *sinks, *gparts, *Nsinks,
+                                Ndm + Ndm_background + *Ngas);
+
     /* Duplicate the stars particles into gparts */
     if (with_stars)
       io_duplicate_stars_gparts(&tp, *sparts, *gparts, *Nstars,
-                                Ndm + Ndm_background + *Ngas);
+                                Ndm + Ndm_background + *Ngas + *Nsinks);
 
     /* Duplicate the black holes particles into gparts */
     if (with_black_holes)
-      io_duplicate_black_holes_gparts(&tp, *bparts, *gparts, *Nblackholes,
-                                      Ndm + Ndm_background + *Ngas + *Nstars);
+      io_duplicate_black_holes_gparts(
+          &tp, *bparts, *gparts, *Nblackholes,
+          Ndm + Ndm_background + *Ngas + *Nsinks + *Nstars);
 
     threadpool_clean(&tp);
   }
@@ -873,6 +902,7 @@ void write_output_serial(struct engine* e,
   const struct gpart* gparts = e->s->gparts;
   const struct spart* sparts = e->s->sparts;
   const struct bpart* bparts = e->s->bparts;
+  const struct sink* sinks = e->s->sinks;
   struct output_options* output_options = e->output_options;
   struct output_list* output_list = e->output_list_snapshots;
   const int with_cosmology = e->policy & engine_policy_cosmology;
@@ -892,6 +922,7 @@ void write_output_serial(struct engine* e,
   /* 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 Nsinks = e->s->nr_sinks;
   const size_t Nstars = e->s->nr_sparts;
   const size_t Nblackholes = e->s->nr_bparts;
   // const size_t Nbaryons = Ngas + Nstars;
@@ -908,12 +939,14 @@ void write_output_serial(struct engine* e,
       e->s->nr_gparts - e->s->nr_inhibited_gparts - e->s->nr_extra_gparts;
   const size_t Ngas_written =
       e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
+  const size_t Nsinks_written =
+      e->s->nr_sinks - e->s->nr_inhibited_sinks - e->s->nr_extra_sinks;
   const size_t Nstars_written =
       e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
   const size_t Nblackholes_written =
       e->s->nr_bparts - e->s->nr_inhibited_bparts - e->s->nr_extra_bparts;
   const size_t Nbaryons_written =
-      Ngas_written + Nstars_written + Nblackholes_written;
+      Ngas_written + Nstars_written + Nblackholes_written + Nsinks_written;
   const size_t Ndm_written =
       Ntot_written > 0 ? Ntot_written - Nbaryons_written - Ndm_background : 0;
 
@@ -945,7 +978,7 @@ void write_output_serial(struct engine* e,
 
   /* Compute offset in the file and total number of particles */
   size_t N[swift_type_count] = {Ngas_written,   Ndm_written,
-                                Ndm_background, 0,
+                                Ndm_background, Nsinks_written,
                                 Nstars_written, Nblackholes_written};
   long long N_total[swift_type_count] = {0};
   long long offset[swift_type_count] = {0};
@@ -1141,6 +1174,7 @@ void write_output_serial(struct engine* e,
         struct velociraptor_gpart_data* gpart_group_data_written = NULL;
         struct spart* sparts_written = NULL;
         struct bpart* bparts_written = NULL;
+        struct sink* sinks_written = NULL;
 
         /* Write particle fields from the particle structure */
         switch (ptype) {
@@ -1307,6 +1341,33 @@ void write_output_serial(struct engine* e,
 
           } break;
 
+          case swift_type_sink: {
+            if (Nsinks == Nsinks_written) {
+
+              /* No inhibted particles: easy case */
+              Nparticles = Nsinks;
+              sink_write_particles(sinks, list, &num_fields, with_cosmology);
+            } else {
+
+              /* Ok, we need to fish out the particles we want */
+              Nparticles = Nsinks_written;
+
+              /* Allocate temporary arrays */
+              if (swift_memalign("sinks_written", (void**)&sinks_written,
+                                 sink_align,
+                                 Nsinks_written * sizeof(struct sink)) != 0)
+                error("Error while allocating temporary memory for sinks");
+
+              /* Collect the particles we want to write */
+              io_collect_sinks_to_write(sinks, sinks_written, Nsinks,
+                                        Nsinks_written);
+
+              /* Select the fields to write */
+              sink_write_particles(sinks_written, list, &num_fields,
+                                   with_cosmology);
+            }
+          } break;
+
           case swift_type_stars: {
             if (Nstars == Nstars_written) {
 
@@ -1450,6 +1511,7 @@ void write_output_serial(struct engine* e,
           swift_free("gpart_group_written", gpart_group_data_written);
         if (sparts_written) swift_free("sparts_written", sparts_written);
         if (bparts_written) swift_free("bparts_written", sparts_written);
+        if (sinks_written) swift_free("sinks_written", sinks_written);
 
         /* Close particle group */
         H5Gclose(h_grp);
diff --git a/src/serial_io.h b/src/serial_io.h
index 764bcf9ca571f06a0e3c9ba389079f072bd61d0f..537340b0b1ff1d77f2b428a4eea9c4429539e12f 100644
--- a/src/serial_io.h
+++ b/src/serial_io.h
@@ -37,13 +37,15 @@ struct unit_system;
 
 void read_ic_serial(char* fileName, const struct unit_system* internal_units,
                     double dim[3], struct part** parts, struct gpart** gparts,
-                    struct spart** sparts, struct bpart** bparts, size_t* Ngas,
-                    size_t* Ngparts, size_t* Ngparts_background, size_t* Nstars,
+                    struct sink** sinks, struct spart** sparts,
+                    struct bpart** bparts, size_t* Ngas, size_t* Ngparts,
+                    size_t* Ngparts_background, size_t* Nsinks, size_t* Nstars,
                     size_t* Nblackholes, int* flag_entropy, int with_hydro,
-                    int with_gravity, int with_stars, int with_black_holes,
-                    int with_cosmology, int cleanup_h, int cleanup_sqrt_a,
-                    double h, double a, int mpi_rank, int mpi_size,
-                    MPI_Comm comm, MPI_Info info, int n_threads, int dry_run);
+                    int with_gravity, int with_sinks, int with_stars,
+                    int with_black_holes, int with_cosmology, int cleanup_h,
+                    int cleanup_sqrt_a, double h, double a, int mpi_rank,
+                    int mpi_size, MPI_Comm comm, MPI_Info info, int n_threads,
+                    int dry_run);
 
 void write_output_serial(struct engine* e,
                          const struct unit_system* internal_units,
diff --git a/src/single_io.c b/src/single_io.c
index 45e79169c99730c378c8ecf3f8503c48d78ee6d6..0657671f75a23715c7968c142b9d1f50a031285e 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -54,6 +54,7 @@
 #include "output_options.h"
 #include "part.h"
 #include "part_type.h"
+#include "sink_io.h"
 #include "star_formation_io.h"
 #include "stars_io.h"
 #include "tracers_io.h"
@@ -376,16 +377,19 @@ void write_array_single(const struct engine* e, hid_t grp, char* fileName,
  * @param dim (output) The dimension of the volume.
  * @param parts (output) Array of #part particles.
  * @param gparts (output) Array of #gpart particles.
+ * @param sinks (output) Array of #sink particles.
  * @param sparts (output) Array of #spart particles.
  * @param bparts (output) Array of #bpart particles.
  * @param Ngas (output) number of Gas particles read.
  * @param Ngparts (output) The number of #gpart read.
+ * @param Nsinks (output) The number of #sink read.
  * @param Nstars (output) The number of #spart read.
  * @param Nblackholes (output) The number of #bpart read.
  * @param flag_entropy (output) 1 if the ICs contained Entropy in the
  * InternalEnergy field
  * @param with_hydro Are we reading gas particles ?
  * @param with_gravity Are we reading/creating #gpart arrays ?
+ * @param with_sink Are we reading sink particles ?
  * @param with_stars Are we reading star particles ?
  * @param with_black_hole Are we reading black hole particles ?
  * @param with_cosmology Are we running with cosmology ?
@@ -407,12 +411,14 @@ void write_array_single(const struct engine* e, hid_t grp, char* fileName,
 void read_ic_single(const char* fileName,
                     const struct unit_system* internal_units, double dim[3],
                     struct part** parts, struct gpart** gparts,
-                    struct spart** sparts, struct bpart** bparts, size_t* Ngas,
-                    size_t* Ngparts, size_t* Ngparts_background, size_t* Nstars,
+                    struct sink** sinks, struct spart** sparts,
+                    struct bpart** bparts, size_t* Ngas, size_t* Ngparts,
+                    size_t* Ngparts_background, size_t* Nsinks, size_t* Nstars,
                     size_t* Nblackholes, int* flag_entropy, int with_hydro,
-                    int with_gravity, int with_stars, int with_black_holes,
-                    int with_cosmology, int cleanup_h, int cleanup_sqrt_a,
-                    double h, double a, int n_threads, int dry_run) {
+                    int with_gravity, int with_sink, int with_stars,
+                    int with_black_holes, int with_cosmology, int cleanup_h,
+                    int cleanup_sqrt_a, double h, double a, int n_threads,
+                    int dry_run) {
 
   hid_t h_file = 0, h_grp = 0;
   /* GADGET has only cubic boxes (in cosmological mode) */
@@ -427,7 +433,7 @@ void read_ic_single(const char* fileName,
 
   /* Initialise counters */
   *Ngas = 0, *Ngparts = 0, *Ngparts_background = 0, *Nstars = 0,
-  *Nblackholes = 0;
+  *Nblackholes = 0, *Nsinks = 0;
 
   /* Open file */
   /* message("Opening file '%s' as IC.", fileName); */
@@ -555,6 +561,15 @@ void read_ic_single(const char* fileName,
     bzero(*parts, *Ngas * sizeof(struct part));
   }
 
+  /* Allocate memory to store star particles */
+  if (with_sink) {
+    *Nsinks = N[swift_type_sink];
+    if (swift_memalign("sinks", (void**)sinks, sink_align,
+                       *Nsinks * sizeof(struct sink)) != 0)
+      error("Error while allocating memory for sink particles");
+    bzero(*sinks, *Nsinks * sizeof(struct sink));
+  }
+
   /* Allocate memory to store star particles */
   if (with_stars) {
     *Nstars = N[swift_type_stars];
@@ -569,7 +584,7 @@ void read_ic_single(const char* fileName,
     *Nblackholes = N[swift_type_black_hole];
     if (swift_memalign("bparts", (void**)bparts, bpart_align,
                        *Nblackholes * sizeof(struct bpart)) != 0)
-      error("Error while allocating memory for stars particles");
+      error("Error while allocating memory for black hole particles");
     bzero(*bparts, *Nblackholes * sizeof(struct bpart));
   }
 
@@ -580,6 +595,7 @@ void read_ic_single(const char* fileName,
     *Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
                N[swift_type_dark_matter] +
                N[swift_type_dark_matter_background] +
+               (with_sink ? N[swift_type_sink] : 0) +
                (with_stars ? N[swift_type_stars] : 0) +
                (with_black_holes ? N[swift_type_black_hole] : 0);
     *Ngparts_background = Ndm_background;
@@ -638,6 +654,13 @@ void read_ic_single(const char* fileName,
         }
         break;
 
+      case swift_type_sink:
+        if (with_sink) {
+          Nparticles = *Nsinks;
+          sink_read_particles(*sinks, list, &num_fields);
+        }
+        break;
+
       case swift_type_stars:
         if (with_stars) {
           Nparticles = *Nstars;
@@ -684,15 +707,21 @@ void read_ic_single(const char* fileName,
       io_duplicate_hydro_gparts(&tp, *parts, *gparts, *Ngas,
                                 Ndm + Ndm_background);
 
+    /* Duplicate the star particles into gparts */
+    if (with_sink)
+      io_duplicate_sinks_gparts(&tp, *sinks, *gparts, *Nsinks,
+                                Ndm + Ndm_background + *Ngas);
+
     /* Duplicate the star particles into gparts */
     if (with_stars)
       io_duplicate_stars_gparts(&tp, *sparts, *gparts, *Nstars,
-                                Ndm + Ndm_background + *Ngas);
+                                Ndm + Ndm_background + *Ngas + *Nsinks);
 
     /* Duplicate the black hole particles into gparts */
     if (with_black_holes)
-      io_duplicate_black_holes_gparts(&tp, *bparts, *gparts, *Nblackholes,
-                                      Ndm + Ndm_background + *Ngas + *Nstars);
+      io_duplicate_black_holes_gparts(
+          &tp, *bparts, *gparts, *Nblackholes,
+          Ndm + Ndm_background + *Ngas + *Nsinks + *Nstars);
 
     threadpool_clean(&tp);
   }
@@ -732,6 +761,7 @@ void write_output_single(struct engine* e,
   const struct gpart* gparts = e->s->gparts;
   const struct spart* sparts = e->s->sparts;
   const struct bpart* bparts = e->s->bparts;
+  const struct sink* sinks = e->s->sinks;
   struct output_options* output_options = e->output_options;
   struct output_list* output_list = e->output_list_snapshots;
   const int with_cosmology = e->policy & engine_policy_cosmology;
@@ -750,6 +780,7 @@ void write_output_single(struct engine* e,
   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 Nsinks = e->s->nr_sinks;
   const size_t Nblackholes = e->s->nr_bparts;
   // const size_t Nbaryons = Ngas + Nstars;
   // const size_t Ndm = Ntot > 0 ? Ntot - Nbaryons : 0;
@@ -767,17 +798,19 @@ void write_output_single(struct engine* e,
       e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
   const size_t Nstars_written =
       e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
+  const size_t Nsinks_written =
+      e->s->nr_sinks - e->s->nr_inhibited_sinks - e->s->nr_extra_sinks;
   const size_t Nblackholes_written =
       e->s->nr_bparts - e->s->nr_inhibited_bparts - e->s->nr_extra_bparts;
   const size_t Nbaryons_written =
-      Ngas_written + Nstars_written + Nblackholes_written;
+      Ngas_written + Nstars_written + Nblackholes_written + Nsinks_written;
   const size_t Ndm_written =
       Ntot_written > 0 ? Ntot_written - Nbaryons_written - Ndm_background : 0;
 
   /* Format things in a Gadget-friendly array */
   long long N_total[swift_type_count] = {
       (long long)Ngas_written,   (long long)Ndm_written,
-      (long long)Ndm_background, 0,
+      (long long)Ndm_background, (long long)Nsinks,
       (long long)Nstars_written, (long long)Nblackholes_written};
 
   /* File name */
@@ -935,6 +968,7 @@ void write_output_single(struct engine* e,
     struct velociraptor_gpart_data* gpart_group_data_written = NULL;
     struct spart* sparts_written = NULL;
     struct bpart* bparts_written = NULL;
+    struct sink* sinks_written = NULL;
 
     /* Write particle fields from the particle structure */
     switch (ptype) {
@@ -1094,6 +1128,33 @@ void write_output_single(struct engine* e,
         }
       } break;
 
+      case swift_type_sink: {
+        if (Nsinks == Nsinks_written) {
+
+          /* No inhibted particles: easy case */
+          N = Nsinks;
+          sink_write_particles(sinks, list, &num_fields, with_cosmology);
+        } else {
+
+          /* Ok, we need to fish out the particles we want */
+          N = Nsinks_written;
+
+          /* Allocate temporary arrays */
+          if (swift_memalign("sinks_written", (void**)&sinks_written,
+                             sink_align,
+                             Nsinks_written * sizeof(struct sink)) != 0)
+            error("Error while allocating temporary memory for sinks");
+
+          /* Collect the particles we want to write */
+          io_collect_sinks_to_write(sinks, sinks_written, Nsinks,
+                                    Nsinks_written);
+
+          /* Select the fields to write */
+          sink_write_particles(sinks_written, list, &num_fields,
+                               with_cosmology);
+        }
+      } break;
+
       case swift_type_stars: {
         if (Nstars == Nstars_written) {
 
@@ -1227,6 +1288,7 @@ void write_output_single(struct engine* e,
       swift_free("gpart_group_written", gpart_group_data_written);
     if (sparts_written) swift_free("sparts_written", sparts_written);
     if (bparts_written) swift_free("bparts_written", bparts_written);
+    if (sinks_written) swift_free("sinks_written", sinks_written);
 
     /* Close particle group */
     H5Gclose(h_grp);
diff --git a/src/single_io.h b/src/single_io.h
index 24e6bd4e3d0d4181a8e73ae838931256e5687f8d..5f03fb96696278b1236f9335d7bb80e48b767937 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -33,12 +33,14 @@ struct unit_system;
 void read_ic_single(const char* fileName,
                     const struct unit_system* internal_units, double dim[3],
                     struct part** parts, struct gpart** gparts,
-                    struct spart** sparts, struct bpart** bparts, size_t* Ngas,
-                    size_t* Ndm, size_t* Ndm_background, size_t* Nstars,
+                    struct sink** sinks, struct spart** sparts,
+                    struct bpart** bparts, size_t* Ngas, size_t* Ndm,
+                    size_t* Ndm_background, size_t* Nsinks, size_t* Nstars,
                     size_t* Nblackholes, int* flag_entropy, int with_hydro,
-                    int with_gravity, int with_stars, int with_black_holes,
-                    int with_cosmology, int cleanup_h, int cleanup_sqrt_a,
-                    double h, double a, int nr_threads, int dry_run);
+                    int with_gravity, int with_sinks, int with_stars,
+                    int with_black_holes, int with_cosmology, int cleanup_h,
+                    int cleanup_sqrt_a, double h, double a, int nr_threads,
+                    int dry_run);
 
 void write_output_single(struct engine* e,
                          const struct unit_system* internal_units,
diff --git a/src/sink/Default/sink.h b/src/sink/Default/sink.h
index b2c9e2a240443e08f4dfa72903a50d33c6023b45..75748d3d502ef961858ff6de4e2caa340c0ed7a3 100644
--- a/src/sink/Default/sink.h
+++ b/src/sink/Default/sink.h
@@ -22,8 +22,6 @@
 #include <float.h>
 
 /* Local includes */
-#include "dimension.h"
-#include "kernel_hydro.h"
 #include "minmax.h"
 #include "sink_part.h"
 
@@ -45,7 +43,6 @@ __attribute__((always_inline)) INLINE static float sink_compute_timestep(
  * read in to do some conversions.
  *
  * @param sp The particle to act upon
- * @param props The properties of the sink model.
  */
 __attribute__((always_inline)) INLINE static void sink_first_init_sink(
     struct sink* sp) {
diff --git a/src/space.c b/src/space.c
index 60b3c3975c4c2170156a80057dd00a04dc48153c..6aeddd5a168f1ec80494a3a25efee9162ca808aa 100644
--- a/src/space.c
+++ b/src/space.c
@@ -59,6 +59,7 @@
 #include "pressure_floor.h"
 #include "proxy.h"
 #include "restart.h"
+#include "sink.h"
 #include "sort_part.h"
 #include "space_unique_id.h"
 #include "star_formation.h"
@@ -91,6 +92,9 @@ int space_extra_bparts = space_extra_bparts_default;
 /*! Number of extra #gpart we allocate memory for per top-level cell */
 int space_extra_gparts = space_extra_gparts_default;
 
+/*! Number of extra #sink we allocate memory for per top-level cell */
+int space_extra_sinks = space_extra_sinks_default;
+
 /*! Maximum number of particles per ghost */
 int engine_max_parts_per_ghost = engine_max_parts_per_ghost_default;
 int engine_max_sparts_per_ghost = engine_max_sparts_per_ghost_default;
@@ -125,10 +129,12 @@ struct index_data {
   size_t count_inhibited_gpart;
   size_t count_inhibited_spart;
   size_t count_inhibited_bpart;
+  size_t count_inhibited_sink;
   size_t count_extra_part;
   size_t count_extra_gpart;
   size_t count_extra_spart;
   size_t count_extra_bpart;
+  size_t count_extra_sink;
 };
 
 /**
@@ -212,6 +218,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
     c->grav.count = 0;
     c->grav.count_total = 0;
     c->grav.updated = 0;
+    c->sinks.count = 0;
     c->stars.count = 0;
     c->stars.count_total = 0;
     c->stars.updated = 0;
@@ -269,6 +276,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
     c->hydro.xparts = NULL;
     c->grav.parts = NULL;
     c->grav.parts_rebuild = NULL;
+    c->sinks.parts = NULL;
     c->stars.parts = NULL;
     c->stars.parts_rebuild = NULL;
     c->black_holes.parts = NULL;
@@ -734,12 +742,14 @@ void space_allocate_extras(struct space *s, int verbose) {
   size_t nr_gparts = s->nr_gparts;
   size_t nr_sparts = s->nr_sparts;
   size_t nr_bparts = s->nr_bparts;
+  size_t nr_sinks = s->nr_sinks;
 
   /* The current number of actual particles */
   size_t nr_actual_parts = nr_parts - s->nr_extra_parts;
   size_t nr_actual_gparts = nr_gparts - s->nr_extra_gparts;
   size_t nr_actual_sparts = nr_sparts - s->nr_extra_sparts;
   size_t nr_actual_bparts = nr_bparts - s->nr_extra_bparts;
+  size_t nr_actual_sinks = nr_sinks - s->nr_extra_sinks;
 
   /* The number of particles we allocated memory for (MPI overhead) */
   size_t size_parts = s->size_parts;
@@ -765,17 +775,21 @@ void space_allocate_extras(struct space *s, int verbose) {
   const size_t expected_num_extra_gparts = nr_local_cells * space_extra_gparts;
   const size_t expected_num_extra_sparts = nr_local_cells * space_extra_sparts;
   const size_t expected_num_extra_bparts = nr_local_cells * space_extra_bparts;
+  const size_t expected_num_extra_sinks = nr_local_cells * space_extra_sinks;
 
   if (verbose) {
-    message("Currently have %zd/%zd/%zd/%zd real particles.", nr_actual_parts,
-            nr_actual_gparts, nr_actual_sparts, nr_actual_bparts);
-    message("Currently have %zd/%zd/%zd/%zd spaces for extra particles.",
-            s->nr_extra_parts, s->nr_extra_gparts, s->nr_extra_sparts,
-            s->nr_extra_bparts);
+    message("Currently have %zd/%zd/%zd/%zd/%zd real particles.",
+            nr_actual_parts, nr_actual_gparts, nr_actual_sinks,
+            nr_actual_sparts, nr_actual_bparts);
+    message("Currently have %zd/%zd/%zd/%zd/%zd spaces for extra particles.",
+            s->nr_extra_parts, s->nr_extra_gparts, s->nr_extra_sinks,
+            s->nr_extra_sparts, s->nr_extra_bparts);
     message(
-        "Requesting space for future %zd/%zd/%zd/%zd part/gpart/sparts/bparts.",
+        "Requesting space for future %zd/%zd/%zd/%zd/%zd "
+        "part/gpart/sinks/sparts/bparts.",
         expected_num_extra_parts, expected_num_extra_gparts,
-        expected_num_extra_sparts, expected_num_extra_bparts);
+        expected_num_extra_sinks, expected_num_extra_sparts,
+        expected_num_extra_bparts);
   }
 
   if (expected_num_extra_parts < s->nr_extra_parts)
@@ -786,6 +800,8 @@ void space_allocate_extras(struct space *s, int verbose) {
     error("Reduction in top-level cells number not handled.");
   if (expected_num_extra_bparts < s->nr_extra_bparts)
     error("Reduction in top-level cells number not handled.");
+  if (expected_num_extra_sinks < s->nr_extra_sinks)
+    error("Reduction in top-level cells number not handled.");
 
   /* Do we have enough space for the extra gparts (i.e. we haven't used up any)
    * ? */
@@ -977,6 +993,12 @@ void space_allocate_extras(struct space *s, int verbose) {
     s->nr_extra_parts = expected_num_extra_parts;
   }
 
+  /* Do we have enough space for the extra sinks (i.e. we haven't used up any)
+   * ? */
+  if (nr_actual_sinks + expected_num_extra_sinks > nr_sinks) {
+    error("Not implemented yet");
+  }
+
   /* Do we have enough space for the extra sparts (i.e. we haven't used up any)
    * ? */
   if (nr_actual_sparts + expected_num_extra_sparts > nr_sparts) {
@@ -1144,9 +1166,10 @@ void space_allocate_extras(struct space *s, int verbose) {
 #ifdef SWIFT_DEBUG_CHECKS
   /* Verify that the links are correct */
   if ((nr_gparts > 0 && nr_parts > 0) || (nr_gparts > 0 && nr_sparts > 0) ||
-      (nr_gparts > 0 && nr_bparts > 0))
-    part_verify_links(s->parts, s->gparts, s->sparts, s->bparts, nr_parts,
-                      nr_gparts, nr_sparts, nr_bparts, verbose);
+      (nr_gparts > 0 && nr_bparts > 0) || (nr_gparts > 0 && nr_sinks > 0))
+    part_verify_links(s->parts, s->gparts, s->sinks, s->sparts, s->bparts,
+                      nr_parts, nr_gparts, nr_sinks, nr_sparts, nr_bparts,
+                      verbose);
 #endif
 
   /* Free the list of local cells */
@@ -1231,30 +1254,35 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
   size_t nr_gparts = s->nr_gparts;
   size_t nr_sparts = s->nr_sparts;
   size_t nr_bparts = s->nr_bparts;
+  size_t nr_sinks = s->nr_sinks;
 
   /* The number of particles we allocated memory for */
   size_t size_parts = s->size_parts;
   size_t size_gparts = s->size_gparts;
   size_t size_sparts = s->size_sparts;
   size_t size_bparts = s->size_bparts;
+  size_t size_sinks = s->size_sinks;
 
   /* Counter for the number of inhibited particles found on the node */
   size_t count_inhibited_parts = 0;
   size_t count_inhibited_gparts = 0;
   size_t count_inhibited_sparts = 0;
   size_t count_inhibited_bparts = 0;
+  size_t count_inhibited_sinks = 0;
 
   /* Counter for the number of extra particles found on the node */
   size_t count_extra_parts = 0;
   size_t count_extra_gparts = 0;
   size_t count_extra_sparts = 0;
   size_t count_extra_bparts = 0;
+  size_t count_extra_sinks = 0;
 
   /* Number of particles we expect to have after strays exchange */
   const size_t h_index_size = size_parts + space_expected_max_nr_strays;
   const size_t g_index_size = size_gparts + space_expected_max_nr_strays;
   const size_t s_index_size = size_sparts + space_expected_max_nr_strays;
   const size_t b_index_size = size_bparts + space_expected_max_nr_strays;
+  const size_t sink_index_size = size_sinks + space_expected_max_nr_strays;
 
   /* Allocate arrays to store the indices of the cells where particles
      belong. We allocate extra space to allow for particles we may
@@ -1263,7 +1291,10 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
   int *g_index = (int *)swift_malloc("g_index", sizeof(int) * g_index_size);
   int *s_index = (int *)swift_malloc("s_index", sizeof(int) * s_index_size);
   int *b_index = (int *)swift_malloc("b_index", sizeof(int) * b_index_size);
-  if (h_index == NULL || g_index == NULL || s_index == NULL || b_index == NULL)
+  int *sink_index =
+      (int *)swift_malloc("sink_index", sizeof(int) * sink_index_size);
+  if (h_index == NULL || g_index == NULL || s_index == NULL ||
+      b_index == NULL || sink_index == NULL)
     error("Failed to allocate temporary particle indices.");
 
   /* Allocate counters of particles that will land in each cell */
@@ -1275,8 +1306,12 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
       (int *)swift_malloc("cell_spart_counts", sizeof(int) * s->nr_cells);
   int *cell_bpart_counts =
       (int *)swift_malloc("cell_bpart_counts", sizeof(int) * s->nr_cells);
+  int *cell_sink_counts =
+      (int *)swift_malloc("cell_sink_counts", sizeof(int) * s->nr_cells);
+
   if (cell_part_counts == NULL || cell_gpart_counts == NULL ||
-      cell_spart_counts == NULL || cell_bpart_counts == NULL)
+      cell_spart_counts == NULL || cell_bpart_counts == NULL ||
+      cell_sink_counts == NULL)
     error("Failed to allocate cell particle count buffer.");
 
   /* Initialise the counters, including buffer space for future particles */
@@ -1285,6 +1320,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
     cell_gpart_counts[i] = 0;
     cell_spart_counts[i] = 0;
     cell_bpart_counts[i] = 0;
+    cell_sink_counts[i] = 0;
   }
 
   /* Run through the particles and get their cell index. */
@@ -1304,6 +1340,10 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
     space_bparts_get_cell_index(s, b_index, cell_bpart_counts,
                                 &count_inhibited_bparts, &count_extra_bparts,
                                 verbose);
+  if (nr_sinks > 0)
+    space_sinks_get_cell_index(s, sink_index, cell_sink_counts,
+                               &count_inhibited_sinks, &count_extra_sinks,
+                               verbose);
 
 #ifdef SWIFT_DEBUG_CHECKS
   /* Some safety checks */
@@ -1315,6 +1355,8 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
     error("We just repartitioned but still found inhibited gparts.");
   if (repartitioned && count_inhibited_bparts)
     error("We just repartitioned but still found inhibited bparts.");
+  if (repartitioned && count_inhibited_sinks)
+    error("We just repartitioned but still found inhibited sinks.");
 
   if (count_extra_parts != s->nr_extra_parts)
     error(
@@ -1332,6 +1374,10 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
     error(
         "Number of extra bparts in the bpart array not matching the space "
         "counter.");
+  if (count_extra_sinks != s->nr_extra_sinks)
+    error(
+        "Number of extra sinks in the sink array not matching the space "
+        "counter.");
 #endif
 
   const ticks tic2 = getticks();
@@ -1491,6 +1537,60 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
     error("Counts of inhibited b-particles do not match!");
 #endif /* SWIFT_DEBUG_CHECKS */
 
+  /* Move non-local sinks and inhibited sinks to the end of the list. */
+  if ((with_dithering || !repartitioned) &&
+      (s->e->nr_nodes > 1 || count_inhibited_sinks > 0)) {
+
+    for (size_t k = 0; k < nr_sinks; /* void */) {
+
+      /* Inhibited particle or foreign particle */
+      if (sink_index[k] == -1 ||
+          cells_top[sink_index[k]].nodeID != local_nodeID) {
+
+        /* One fewer particle */
+        nr_sinks -= 1;
+
+        /* Swap the particle */
+        memswap(&s->sinks[k], &s->sinks[nr_sinks], sizeof(struct sink));
+
+        /* Swap the link with the gpart */
+        if (s->sinks[k].gpart != NULL) {
+          s->sinks[k].gpart->id_or_neg_offset = -k;
+        }
+        if (s->sinks[nr_sinks].gpart != NULL) {
+          s->sinks[nr_sinks].gpart->id_or_neg_offset = -nr_sinks;
+        }
+
+        /* Swap the index */
+        memswap(&sink_index[k], &sink_index[nr_sinks], sizeof(int));
+
+      } else {
+        /* Increment when not exchanging otherwise we need to retest "k".*/
+        k++;
+      }
+    }
+  }
+
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Check that all sinks are in the correct place. */
+  size_t check_count_inhibited_sinks = 0;
+  for (size_t k = 0; k < nr_sinks; k++) {
+    if (sink_index[k] == -1 ||
+        cells_top[sink_index[k]].nodeID != local_nodeID) {
+      error("Failed to move all non-local sinks to send list");
+    }
+  }
+  for (size_t k = nr_sinks; k < s->nr_sinks; k++) {
+    if (sink_index[k] != -1 &&
+        cells_top[sink_index[k]].nodeID == local_nodeID) {
+      error("Failed to remove local sinks from send list");
+    }
+    if (sink_index[k] == -1) ++check_count_inhibited_sinks;
+  }
+  if (check_count_inhibited_sinks != count_inhibited_sinks)
+    error("Counts of inhibited sink-particles do not match!");
+#endif /* SWIFT_DEBUG_CHECKS */
+
   /* Move non-local gparts and inhibited parts to the end of the list. */
   if ((with_dithering || !repartitioned) &&
       (s->e->nr_nodes > 1 || count_inhibited_gparts > 0)) {
@@ -1685,6 +1785,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
   s->nr_parts = nr_parts;
   s->nr_sparts = nr_sparts;
   s->nr_bparts = nr_bparts;
+  s->nr_sinks = nr_sinks;
 
 #endif /* WITH_MPI */
 
@@ -1779,6 +1880,36 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
   }
 #endif /* SWIFT_DEBUG_CHECKS */
 
+  /* Sort the sink according to their cells. */
+  if (nr_sinks > 0)
+    space_sinks_sort(s->sinks, sink_index, cell_sink_counts, s->nr_cells, 0);
+
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Verify that the sink have been sorted correctly. */
+  for (size_t k = 0; k < nr_sinks; k++) {
+    const struct sink *sink = &s->sinks[k];
+
+    if (sink->time_bin == time_bin_inhibited)
+      error("Inhibited particle sorted into a cell!");
+
+    /* New cell index */
+    const int new_bind =
+        cell_getid(s->cdim, sink->x[0] * s->iwidth[0],
+                   sink->x[1] * s->iwidth[1], sink->x[2] * s->iwidth[2]);
+
+    /* New cell of this sink */
+    const struct cell *c = &s->cells_top[new_bind];
+
+    if (sink_index[k] != new_bind)
+      error("sink's new cell index not matching sorted index.");
+
+    if (sink->x[0] < c->loc[0] || sink->x[0] > c->loc[0] + c->width[0] ||
+        sink->x[1] < c->loc[1] || sink->x[1] > c->loc[1] + c->width[1] ||
+        sink->x[2] < c->loc[2] || sink->x[2] > c->loc[2] + c->width[2])
+      error("sink not sorted into the right top-level cell!");
+  }
+#endif /* SWIFT_DEBUG_CHECKS */
+
   /* Extract the cell counts from the sorted indices. Deduct the extra
    * particles. */
   size_t last_index = 0;
@@ -1815,6 +1946,18 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
     }
   }
 
+  /* Extract the cell counts from the sorted indices. Deduct the extra
+   * particles. */
+  size_t last_sink_index = 0;
+  sink_index[nr_sinks] = s->nr_cells;  // sentinel.
+  for (size_t k = 0; k < nr_sinks; k++) {
+    if (sink_index[k] < sink_index[k + 1]) {
+      cells_top[sink_index[k]].sinks.count =
+          k - last_sink_index + 1 - space_extra_sinks;
+      last_sink_index = k + 1;
+    }
+  }
+
   /* We no longer need the indices as of here. */
   swift_free("h_index", h_index);
   swift_free("cell_part_counts", cell_part_counts);
@@ -1822,6 +1965,8 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
   swift_free("cell_spart_counts", cell_spart_counts);
   swift_free("b_index", b_index);
   swift_free("cell_bpart_counts", cell_bpart_counts);
+  swift_free("sink_index", sink_index);
+  swift_free("cell_sink_counts", cell_sink_counts);
 
   /* Update the slice of unique IDs. */
   space_update_unique_id(s);
@@ -1865,11 +2010,12 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
   s->nr_inhibited_gparts = 0;
   s->nr_inhibited_sparts = 0;
   s->nr_inhibited_bparts = 0;
+  s->nr_inhibited_sinks = 0;
 
   /* Sort the gparts according to their cells. */
   if (nr_gparts > 0)
-    space_gparts_sort(s->gparts, s->parts, s->sparts, s->bparts, g_index,
-                      cell_gpart_counts, s->nr_cells);
+    space_gparts_sort(s->gparts, s->parts, s->sinks, s->sparts, s->bparts,
+                      g_index, cell_gpart_counts, s->nr_cells);
 
 #ifdef SWIFT_DEBUG_CHECKS
   /* Verify that the gpart have been sorted correctly. */
@@ -1916,9 +2062,10 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
 #ifdef SWIFT_DEBUG_CHECKS
   /* Verify that the links are correct */
   if ((nr_gparts > 0 && nr_parts > 0) || (nr_gparts > 0 && nr_sparts > 0) ||
-      (nr_gparts > 0 && nr_bparts > 0))
-    part_verify_links(s->parts, s->gparts, s->sparts, s->bparts, nr_parts,
-                      nr_gparts, nr_sparts, nr_bparts, verbose);
+      (nr_gparts > 0 && nr_bparts > 0) || (nr_gparts > 0 && nr_sinks > 0))
+    part_verify_links(s->parts, s->gparts, s->sinks, s->sparts, s->bparts,
+                      nr_parts, nr_gparts, nr_sinks, nr_sparts, nr_bparts,
+                      verbose);
 #endif
 
   /* Hook the cells up to the parts. Make list of local and non-empty cells */
@@ -1928,6 +2075,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
   struct gpart *gfinger = s->gparts;
   struct spart *sfinger = s->sparts;
   struct bpart *bfinger = s->bparts;
+  struct sink *sink_finger = s->sinks;
   s->nr_cells_with_particles = 0;
   s->nr_local_cells_with_particles = 0;
   s->nr_local_cells = 0;
@@ -1945,9 +2093,9 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
 #endif
 
     const int is_local = (c->nodeID == engine_rank);
-    const int has_particles = (c->hydro.count > 0) || (c->grav.count > 0) ||
-                              (c->stars.count > 0) ||
-                              (c->black_holes.count > 0);
+    const int has_particles =
+        (c->hydro.count > 0) || (c->grav.count > 0) || (c->stars.count > 0) ||
+        (c->black_holes.count > 0) || (c->sinks.count > 0);
 
     if (is_local) {
       c->hydro.parts = finger;
@@ -1955,6 +2103,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
       c->grav.parts = gfinger;
       c->stars.parts = sfinger;
       c->black_holes.parts = bfinger;
+      c->sinks.parts = sink_finger;
 
       /* Store the state at rebuild time */
       c->stars.parts_rebuild = c->stars.parts;
@@ -1970,6 +2119,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
       gfinger = &gfinger[c->grav.count_total];
       sfinger = &sfinger[c->stars.count_total];
       bfinger = &bfinger[c->black_holes.count_total];
+      sink_finger = &sink_finger[c->sinks.count_total];
 
       /* Add this cell to the list of local cells */
       s->local_cells_top[s->nr_local_cells] = k;
@@ -2109,6 +2259,10 @@ void space_reorder_extras(struct space *s, int verbose) {
   /* Re-order the black hole particles */
   if (space_extra_bparts)
     error("Missing implementation of BH extra reordering");
+
+  /* Re-order the sink particles */
+  if (space_extra_sinks)
+    error("Missing implementation of sink extra reordering");
 }
 
 /**
@@ -2713,6 +2867,134 @@ void space_bparts_get_cell_index_mapper(void *map_data, int nr_bparts,
   atomic_add_f(&s->sum_bpart_vel_norm, sum_vel_norm);
 }
 
+/**
+ * @brief #threadpool mapper function to compute the sink-particle cell indices.
+ *
+ * @param map_data Pointer towards the sink-particles.
+ * @param nr_sinks The number of sink-particles to treat.
+ * @param extra_data Pointers to the space and index list
+ */
+void space_sinks_get_cell_index_mapper(void *map_data, int nr_sinks,
+                                       void *extra_data) {
+
+  /* Unpack the data */
+  struct sink *restrict sinks = (struct sink *)map_data;
+  struct index_data *data = (struct index_data *)extra_data;
+  struct space *s = data->s;
+  int *const ind = data->ind + (ptrdiff_t)(sinks - s->sinks);
+
+  /* Get some constants */
+  const int periodic = s->periodic;
+  const int dithering = s->e->gravity_properties->with_dithering;
+  const double delta_dithering_x =
+      s->pos_dithering[0] - s->pos_dithering_old[0];
+  const double delta_dithering_y =
+      s->pos_dithering[1] - s->pos_dithering_old[1];
+  const double delta_dithering_z =
+      s->pos_dithering[2] - s->pos_dithering_old[2];
+  const double dim_x = s->dim[0];
+  const double dim_y = s->dim[1];
+  const double dim_z = s->dim[2];
+  const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
+  const double ih_x = s->iwidth[0];
+  const double ih_y = s->iwidth[1];
+  const double ih_z = s->iwidth[2];
+
+  /* Init the local count buffer. */
+  int *cell_counts = (int *)calloc(sizeof(int), s->nr_cells);
+  if (cell_counts == NULL)
+    error("Failed to allocate temporary cell count buffer.");
+
+  /* Init the local collectors */
+  size_t count_inhibited_sink = 0;
+  size_t count_extra_sink = 0;
+
+  for (int k = 0; k < nr_sinks; k++) {
+
+    /* Get the particle */
+    struct sink *restrict sink = &sinks[k];
+
+    double old_pos_x = sink->x[0];
+    double old_pos_y = sink->x[1];
+    double old_pos_z = sink->x[2];
+
+    if (periodic && dithering && sink->time_bin != time_bin_not_created) {
+      old_pos_x += delta_dithering_x;
+      old_pos_y += delta_dithering_y;
+      old_pos_z += delta_dithering_z;
+    }
+
+#ifdef SWIFT_DEBUG_CHECKS
+    if (!periodic && sink->time_bin != time_bin_inhibited) {
+      if (old_pos_x < 0. || old_pos_x > dim_x)
+        error("Particle outside of volume along X.");
+      if (old_pos_y < 0. || old_pos_y > dim_y)
+        error("Particle outside of volume along Y.");
+      if (old_pos_z < 0. || old_pos_z > dim_z)
+        error("Particle outside of volume along Z.");
+    }
+#endif
+
+    /* Put it back into the simulation volume */
+    double pos_x = box_wrap(old_pos_x, 0.0, dim_x);
+    double pos_y = box_wrap(old_pos_y, 0.0, dim_y);
+    double pos_z = box_wrap(old_pos_z, 0.0, dim_z);
+
+    /* Treat the case where a particle was wrapped back exactly onto
+     * the edge because of rounding issues (more accuracy around 0
+     * than around dim) */
+    if (pos_x == dim_x) pos_x = 0.0;
+    if (pos_y == dim_y) pos_y = 0.0;
+    if (pos_z == dim_z) pos_z = 0.0;
+
+    /* Get its cell index */
+    const int index =
+        cell_getid(cdim, pos_x * ih_x, pos_y * ih_y, pos_z * ih_z);
+
+#ifdef SWIFT_DEBUG_CHECKS
+    if (index < 0 || index >= cdim[0] * cdim[1] * cdim[2])
+      error("Invalid index=%d cdim=[%d %d %d] p->x=[%e %e %e]", index, cdim[0],
+            cdim[1], cdim[2], pos_x, pos_y, pos_z);
+
+    if (pos_x >= dim_x || pos_y >= dim_y || pos_z >= dim_z || pos_x < 0. ||
+        pos_y < 0. || pos_z < 0.)
+      error("Particle outside of simulation box. p->x=[%e %e %e]", pos_x, pos_y,
+            pos_z);
+#endif
+
+    /* Is this particle to be removed? */
+    if (sink->time_bin == time_bin_inhibited) {
+      ind[k] = -1;
+      ++count_inhibited_sink;
+    } else if (sink->time_bin == time_bin_not_created) {
+      /* Is this a place-holder for on-the-fly creation? */
+      ind[k] = index;
+      cell_counts[index]++;
+      ++count_extra_sink;
+
+    } else {
+      /* List its top-level cell index */
+      ind[k] = index;
+      cell_counts[index]++;
+
+      /* Update the position */
+      sink->x[0] = pos_x;
+      sink->x[1] = pos_y;
+      sink->x[2] = pos_z;
+    }
+  }
+
+  /* Write the counts back to the global array. */
+  for (int k = 0; k < s->nr_cells; k++)
+    if (cell_counts[k]) atomic_add(&data->cell_counts[k], cell_counts[k]);
+  free(cell_counts);
+
+  /* Write the count of inhibited and extra bparts */
+  if (count_inhibited_sink)
+    atomic_add(&data->count_inhibited_sink, count_inhibited_sink);
+  if (count_extra_sink) atomic_add(&data->count_extra_sink, count_extra_sink);
+}
+
 /**
  * @brief Computes the cell index of all the particles.
  *
@@ -2745,10 +3027,12 @@ void space_parts_get_cell_index(struct space *s, int *ind, int *cell_counts,
   data.count_inhibited_gpart = 0;
   data.count_inhibited_spart = 0;
   data.count_inhibited_bpart = 0;
+  data.count_inhibited_sink = 0;
   data.count_extra_part = 0;
   data.count_extra_gpart = 0;
   data.count_extra_spart = 0;
   data.count_extra_bpart = 0;
+  data.count_extra_sink = 0;
 
   threadpool_map(&s->e->threadpool, space_parts_get_cell_index_mapper, s->parts,
                  s->nr_parts, sizeof(struct part), threadpool_auto_chunk_size,
@@ -2794,10 +3078,12 @@ void space_gparts_get_cell_index(struct space *s, int *gind, int *cell_counts,
   data.count_inhibited_gpart = 0;
   data.count_inhibited_spart = 0;
   data.count_inhibited_bpart = 0;
+  data.count_inhibited_sink = 0;
   data.count_extra_part = 0;
   data.count_extra_gpart = 0;
   data.count_extra_spart = 0;
   data.count_extra_bpart = 0;
+  data.count_extra_sink = 0;
 
   threadpool_map(&s->e->threadpool, space_gparts_get_cell_index_mapper,
                  s->gparts, s->nr_gparts, sizeof(struct gpart),
@@ -2842,11 +3128,13 @@ void space_sparts_get_cell_index(struct space *s, int *sind, int *cell_counts,
   data.count_inhibited_part = 0;
   data.count_inhibited_gpart = 0;
   data.count_inhibited_spart = 0;
+  data.count_inhibited_sink = 0;
   data.count_inhibited_bpart = 0;
   data.count_extra_part = 0;
   data.count_extra_gpart = 0;
   data.count_extra_spart = 0;
   data.count_extra_bpart = 0;
+  data.count_extra_sink = 0;
 
   threadpool_map(&s->e->threadpool, space_sparts_get_cell_index_mapper,
                  s->sparts, s->nr_sparts, sizeof(struct spart),
@@ -2860,6 +3148,51 @@ void space_sparts_get_cell_index(struct space *s, int *sind, int *cell_counts,
             clocks_getunit());
 }
 
+/**
+ * @brief Computes the cell index of all the sink-particles.
+ *
+ * @param s The #space.
+ * @param sink_ind The array of indices to fill.
+ * @param cell_counts The cell counters to update.
+ * @param count_inhibited_sinks (return) The number of #sink to remove.
+ * @param count_extra_sinks (return) The number of #sink for on-the-fly
+ * creation.
+ * @param verbose Are we talkative ?
+ */
+void space_sinks_get_cell_index(struct space *s, int *sink_ind,
+                                int *cell_counts, size_t *count_inhibited_sinks,
+                                size_t *count_extra_sinks, int verbose) {
+
+  const ticks tic = getticks();
+
+  /* Pack the extra information */
+  struct index_data data;
+  data.s = s;
+  data.ind = sink_ind;
+  data.cell_counts = cell_counts;
+  data.count_inhibited_part = 0;
+  data.count_inhibited_gpart = 0;
+  data.count_inhibited_spart = 0;
+  data.count_inhibited_bpart = 0;
+  data.count_inhibited_sink = 0;
+  data.count_extra_part = 0;
+  data.count_extra_gpart = 0;
+  data.count_extra_spart = 0;
+  data.count_extra_bpart = 0;
+  data.count_extra_sink = 0;
+
+  threadpool_map(&s->e->threadpool, space_sinks_get_cell_index_mapper, s->sinks,
+                 s->nr_sinks, sizeof(struct sink), threadpool_auto_chunk_size,
+                 &data);
+
+  *count_inhibited_sinks = data.count_inhibited_sink;
+  *count_extra_sinks = data.count_extra_sink;
+
+  if (verbose)
+    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+            clocks_getunit());
+}
+
 /**
  * @brief Computes the cell index of all the b-particles.
  *
@@ -2892,10 +3225,12 @@ void space_bparts_get_cell_index(struct space *s, int *bind, int *cell_counts,
   data.count_inhibited_gpart = 0;
   data.count_inhibited_spart = 0;
   data.count_inhibited_bpart = 0;
+  data.count_inhibited_sink = 0;
   data.count_extra_part = 0;
   data.count_extra_gpart = 0;
   data.count_extra_spart = 0;
   data.count_extra_bpart = 0;
+  data.count_extra_sink = 0;
 
   threadpool_map(&s->e->threadpool, space_bparts_get_cell_index_mapper,
                  s->bparts, s->nr_bparts, sizeof(struct bpart),
@@ -3095,11 +3430,72 @@ void space_bparts_sort(struct bpart *bparts, int *restrict ind,
   swift_free("bparts_offsets", offsets);
 }
 
+/**
+ * @brief Sort the sink-particles according to the given indices.
+ *
+ * @param sinks The array of #sink to sort.
+ * @param ind The indices with respect to which the #sink are sorted.
+ * @param counts Number of particles per index.
+ * @param num_bins Total number of bins (length of counts).
+ * @param sinks_offset Offset of the #sink array from the global #sink.
+ * array.
+ */
+void space_sinks_sort(struct sink *sinks, int *restrict ind,
+                      int *restrict counts, int num_bins,
+                      ptrdiff_t sinks_offset) {
+  /* Create the offsets array. */
+  size_t *offsets = NULL;
+  if (swift_memalign("sinks_offsets", (void **)&offsets, SWIFT_STRUCT_ALIGNMENT,
+                     sizeof(size_t) * (num_bins + 1)) != 0)
+    error("Failed to allocate temporary cell offsets array.");
+
+  offsets[0] = 0;
+  for (int k = 1; k <= num_bins; k++) {
+    offsets[k] = offsets[k - 1] + counts[k - 1];
+    counts[k - 1] = 0;
+  }
+
+  /* Loop over local cells. */
+  for (int cid = 0; cid < num_bins; cid++) {
+    for (size_t k = offsets[cid] + counts[cid]; k < offsets[cid + 1]; k++) {
+      counts[cid]++;
+      int target_cid = ind[k];
+      if (target_cid == cid) {
+        continue;
+      }
+      struct sink temp_sink = sinks[k];
+      while (target_cid != cid) {
+        size_t j = offsets[target_cid] + counts[target_cid]++;
+        while (ind[j] == target_cid) {
+          j = offsets[target_cid] + counts[target_cid]++;
+        }
+        memswap(&sinks[j], &temp_sink, sizeof(struct sink));
+        memswap(&ind[j], &target_cid, sizeof(int));
+        if (sinks[j].gpart)
+          sinks[j].gpart->id_or_neg_offset = -(j + sinks_offset);
+      }
+      sinks[k] = temp_sink;
+      ind[k] = target_cid;
+      if (sinks[k].gpart)
+        sinks[k].gpart->id_or_neg_offset = -(k + sinks_offset);
+    }
+  }
+
+#ifdef SWIFT_DEBUG_CHECKS
+  for (int k = 0; k < num_bins; k++)
+    if (offsets[k + 1] != offsets[k] + counts[k])
+      error("Bad offsets after shuffle.");
+#endif /* SWIFT_DEBUG_CHECKS */
+
+  swift_free("sinks_offsets", offsets);
+}
+
 /**
  * @brief Sort the g-particles according to the given indices.
  *
  * @param gparts The array of #gpart to sort.
  * @param parts Global #part array for re-linking.
+ * @param sinks Global #sink array for re-linking.
  * @param sparts Global #spart array for re-linking.
  * @param bparts Global #bpart array for re-linking.
  * @param ind The indices with respect to which the gparts are sorted.
@@ -3107,8 +3503,9 @@ void space_bparts_sort(struct bpart *bparts, int *restrict ind,
  * @param num_bins Total number of bins (length of counts).
  */
 void space_gparts_sort(struct gpart *gparts, struct part *parts,
-                       struct spart *sparts, struct bpart *bparts,
-                       int *restrict ind, int *restrict counts, int num_bins) {
+                       struct sink *sinks, struct spart *sparts,
+                       struct bpart *bparts, int *restrict ind,
+                       int *restrict counts, int num_bins) {
   /* Create the offsets array. */
   size_t *offsets = NULL;
   if (swift_memalign("gparts_offsets", (void **)&offsets,
@@ -3144,6 +3541,8 @@ void space_gparts_sort(struct gpart *gparts, struct part *parts,
           sparts[-gparts[j].id_or_neg_offset].gpart = &gparts[j];
         } else if (gparts[j].type == swift_type_black_hole) {
           bparts[-gparts[j].id_or_neg_offset].gpart = &gparts[j];
+        } else if (gparts[j].type == swift_type_sink) {
+          sinks[-gparts[j].id_or_neg_offset].gpart = &gparts[j];
         }
       }
       gparts[k] = temp_gpart;
@@ -3154,6 +3553,8 @@ void space_gparts_sort(struct gpart *gparts, struct part *parts,
         sparts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
       } else if (gparts[k].type == swift_type_black_hole) {
         bparts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
+      } else if (gparts[k].type == swift_type_sink) {
+        sinks[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
       }
     }
   }
@@ -3329,15 +3730,19 @@ void space_map_cells_pre(struct space *s, int full,
  *        c->black_holes.count or @c NULL.
  * @param gbuff A buffer for particle sorting, should be of size at least
  *        c->grav.count or @c NULL.
+ * @param sink_buff A buffer for particle sorting, should be of size at least
+ *        c->sinks.count or @c NULL.
  */
 void space_split_recursive(struct space *s, struct cell *c,
                            struct cell_buff *buff, struct cell_buff *sbuff,
-                           struct cell_buff *bbuff, struct cell_buff *gbuff) {
+                           struct cell_buff *bbuff, struct cell_buff *gbuff,
+                           struct cell_buff *sink_buff) {
 
   const int count = c->hydro.count;
   const int gcount = c->grav.count;
   const int scount = c->stars.count;
   const int bcount = c->black_holes.count;
+  const int sink_count = c->sinks.count;
   const int with_self_gravity = s->with_self_gravity;
   const int depth = c->depth;
   int maxdepth = 0;
@@ -3357,12 +3762,13 @@ void space_split_recursive(struct space *s, struct cell *c,
   struct spart *sparts = c->stars.parts;
   struct bpart *bparts = c->black_holes.parts;
   struct xpart *xparts = c->hydro.xparts;
+  struct sink *sinks = c->sinks.parts;
   struct engine *e = s->e;
   const integertime_t ti_current = e->ti_current;
 
   /* If the buff is NULL, allocate it, and remember to free it. */
-  const int allocate_buffer =
-      (buff == NULL && gbuff == NULL && sbuff == NULL && bbuff == NULL);
+  const int allocate_buffer = (buff == NULL && gbuff == NULL && sbuff == NULL &&
+                               bbuff == NULL && sink_buff == NULL);
   if (allocate_buffer) {
     if (count > 0) {
       if (swift_memalign("tempbuff", (void **)&buff, SWIFT_STRUCT_ALIGNMENT,
@@ -3428,6 +3834,23 @@ void space_split_recursive(struct space *s, struct cell *c,
         bbuff[k].x[2] = bparts[k].x[2];
       }
     }
+    if (sink_count > 0) {
+      if (swift_memalign("temp_sink_buff", (void **)&sink_buff,
+                         SWIFT_STRUCT_ALIGNMENT,
+                         sizeof(struct cell_buff) * sink_count) != 0)
+        error("Failed to allocate temporary indices.");
+      for (int k = 0; k < sink_count; k++) {
+#ifdef SWIFT_DEBUG_CHECKS
+        if (sinks[k].time_bin == time_bin_inhibited)
+          error("Inhibited particle present in space_split()");
+        if (sinks[k].time_bin == time_bin_not_created)
+          error("Extra particle present in space_split()");
+#endif
+        sink_buff[k].x[0] = sinks[k].x[0];
+        sink_buff[k].x[1] = sinks[k].x[1];
+        sink_buff[k].x[2] = sinks[k].x[2];
+      }
+    }
   }
 
   /* Check the depth. */
@@ -3456,9 +3879,11 @@ void space_split_recursive(struct space *s, struct cell *c,
       cp->hydro.count = 0;
       cp->grav.count = 0;
       cp->stars.count = 0;
+      cp->sinks.count = 0;
       cp->black_holes.count = 0;
       cp->hydro.count_total = 0;
       cp->grav.count_total = 0;
+      cp->sinks.count_total = 0;
       cp->stars.count_total = 0;
       cp->black_holes.count_total = 0;
       cp->hydro.ti_old_part = c->hydro.ti_old_part;
@@ -3504,11 +3929,13 @@ void space_split_recursive(struct space *s, struct cell *c,
 
     /* Split the cell's particle data. */
     cell_split(c, c->hydro.parts - s->parts, c->stars.parts - s->sparts,
-               c->black_holes.parts - s->bparts, buff, sbuff, bbuff, gbuff);
+               c->black_holes.parts - s->bparts, c->sinks.parts - s->sinks,
+               buff, sbuff, bbuff, gbuff, sink_buff);
 
     /* Buffers for the progenitors */
     struct cell_buff *progeny_buff = buff, *progeny_gbuff = gbuff,
-                     *progeny_sbuff = sbuff, *progeny_bbuff = bbuff;
+                     *progeny_sbuff = sbuff, *progeny_bbuff = bbuff,
+                     *progeny_sink_buff = sink_buff;
 
     for (int k = 0; k < 8; k++) {
 
@@ -3517,7 +3944,7 @@ void space_split_recursive(struct space *s, struct cell *c,
 
       /* Remove any progeny with zero particles. */
       if (cp->hydro.count == 0 && cp->grav.count == 0 && cp->stars.count == 0 &&
-          cp->black_holes.count == 0) {
+          cp->black_holes.count == 0 && cp->sinks.count == 0) {
 
         space_recycle(s, cp);
         c->progeny[k] = NULL;
@@ -3526,13 +3953,14 @@ void space_split_recursive(struct space *s, struct cell *c,
 
         /* Recurse */
         space_split_recursive(s, cp, progeny_buff, progeny_sbuff, progeny_bbuff,
-                              progeny_gbuff);
+                              progeny_gbuff, progeny_sink_buff);
 
         /* Update the pointers in the buffers */
         progeny_buff += cp->hydro.count;
         progeny_gbuff += cp->grav.count;
         progeny_sbuff += cp->stars.count;
         progeny_bbuff += cp->black_holes.count;
+        progeny_sink_buff += cp->sinks.count;
 
         /* Update the cell-wide properties */
         h_max = max(h_max, cp->hydro.h_max);
@@ -3851,6 +4279,9 @@ void space_split_recursive(struct space *s, struct cell *c,
   if (s->nr_parts > 0)
     c->owner = ((c->hydro.parts - s->parts) % s->nr_parts) * s->nr_queues /
                s->nr_parts;
+  else if (s->nr_sinks > 0)
+    c->owner = ((c->sinks.parts - s->sinks) % s->nr_sinks) * s->nr_queues /
+               s->nr_sinks;
   else if (s->nr_sparts > 0)
     c->owner = ((c->stars.parts - s->sparts) % s->nr_sparts) * s->nr_queues /
                s->nr_sparts;
@@ -3869,6 +4300,7 @@ void space_split_recursive(struct space *s, struct cell *c,
     if (gbuff != NULL) swift_free("tempgbuff", gbuff);
     if (sbuff != NULL) swift_free("tempsbuff", sbuff);
     if (bbuff != NULL) swift_free("tempbbuff", bbuff);
+    if (sink_buff != NULL) swift_free("temp_sink_buff", sink_buff);
   }
 }
 
@@ -3890,7 +4322,7 @@ void space_split_mapper(void *map_data, int num_cells, void *extra_data) {
   /* Loop over the non-empty cells */
   for (int ind = 0; ind < num_cells; ind++) {
     struct cell *c = &cells_top[local_cells_with_particles[ind]];
-    space_split_recursive(s, c, NULL, NULL, NULL, NULL);
+    space_split_recursive(s, c, NULL, NULL, NULL, NULL, NULL);
   }
 
 #ifdef SWIFT_DEBUG_CHECKS
@@ -4116,7 +4548,7 @@ void space_list_useful_top_level_cells(struct space *s) {
 
     const int has_particles =
         (c->hydro.count > 0) || (c->grav.count > 0) || (c->stars.count > 0) ||
-        (c->black_holes.count > 0) ||
+        (c->black_holes.count > 0) || (c->sinks.count > 0) ||
         (c->grav.multipole != NULL && c->grav.multipole->m_pole.M_000 > 0.f);
 
     if (has_particles) {
@@ -4245,6 +4677,41 @@ void space_synchronize_bpart_positions_mapper(void *map_data, int nr_bparts,
   }
 }
 
+void space_synchronize_sink_positions_mapper(void *map_data, int nr_sinks,
+                                             void *extra_data) {
+  /* Unpack the data */
+  const struct sink *sinks = (struct sink *)map_data;
+
+  for (int k = 0; k < nr_sinks; k++) {
+
+    /* Get the particle */
+    const struct sink *sink = &sinks[k];
+
+    /* Skip unimportant particles */
+    if (sink->time_bin == time_bin_not_created ||
+        sink->time_bin == time_bin_inhibited)
+      continue;
+
+    /* Get its gravity friend */
+    struct gpart *gp = sink->gpart;
+
+#ifdef SWIFT_DEBUG_CHECKS
+    if (gp == NULL) error("Unlinked particle!");
+#endif
+
+    /* Synchronize positions, velocities and masses */
+    gp->x[0] = sink->x[0];
+    gp->x[1] = sink->x[1];
+    gp->x[2] = sink->x[2];
+
+    gp->v_full[0] = sink->v[0];
+    gp->v_full[1] = sink->v[1];
+    gp->v_full[2] = sink->v[2];
+
+    gp->mass = sink->mass;
+  }
+}
+
 /**
  * @brief Make sure the baryon particles are at the same position and
  * have the same velocity and mass as their #gpart friends.
@@ -4272,6 +4739,11 @@ void space_synchronize_particle_positions(struct space *s) {
                    s->bparts, s->nr_bparts, sizeof(struct bpart),
                    threadpool_auto_chunk_size, /*extra_data=*/NULL);
 
+  if (s->nr_gparts > 0 && s->nr_sinks > 0)
+    threadpool_map(&s->e->threadpool, space_synchronize_sink_positions_mapper,
+                   s->sinks, s->nr_sinks, sizeof(struct sink),
+                   threadpool_auto_chunk_size, /*extra_data=*/NULL);
+
   if (s->e->verbose)
     message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
             clocks_getunit());
@@ -4644,6 +5116,71 @@ void space_first_init_bparts(struct space *s, int verbose) {
             clocks_getunit());
 }
 
+void space_first_init_sinks_mapper(void *restrict map_data, int count,
+                                   void *restrict extra_data) {
+
+  struct sink *restrict sink = (struct sink *)map_data;
+  const struct space *restrict s = (struct space *)extra_data;
+  const struct engine *e = s->e;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  const ptrdiff_t delta = sink - s->sinks;
+#endif
+
+  const struct cosmology *cosmo = e->cosmology;
+  const float a_factor_vel = cosmo->a;
+
+  /* Convert velocities to internal units */
+  for (int k = 0; k < count; k++) {
+
+    sink[k].v[0] *= a_factor_vel;
+    sink[k].v[1] *= a_factor_vel;
+    sink[k].v[2] *= a_factor_vel;
+
+#ifdef HYDRO_DIMENSION_2D
+    sink[k].x[2] = 0.f;
+    sink[k].v[2] = 0.f;
+#endif
+
+#ifdef HYDRO_DIMENSION_1D
+    sink[k].x[1] = sink[k].x[2] = 0.f;
+    sink[k].v[1] = sink[k].v[2] = 0.f;
+#endif
+  }
+
+  /* Initialise the rest */
+  for (int k = 0; k < count; k++) {
+
+    sink_first_init_sink(&sink[k]);
+
+#ifdef SWIFT_DEBUG_CHECKS
+    if (sink[k].gpart && sink[k].gpart->id_or_neg_offset != -(k + delta))
+      error("Invalid gpart -> sink link");
+
+    /* Initialise the time-integration check variables */
+    sink[k].ti_drift = 0;
+    sink[k].ti_kick = 0;
+#endif
+  }
+}
+
+/**
+ * @brief Initialises all the sink-particles by setting them into a valid state
+ *
+ * Calls stars_first_init_sink() on all the particles
+ */
+void space_first_init_sinks(struct space *s, int verbose) {
+  const ticks tic = getticks();
+  if (s->nr_sinks > 0)
+    threadpool_map(&s->e->threadpool, space_first_init_sinks_mapper, s->sinks,
+                   s->nr_sinks, sizeof(struct sink), threadpool_auto_chunk_size,
+                   s);
+
+  if (verbose)
+    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+            clocks_getunit());
+}
+
 void space_init_parts_mapper(void *restrict map_data, int count,
                              void *restrict extra_data) {
 
@@ -4767,6 +5304,33 @@ void space_init_bparts(struct space *s, int verbose) {
             clocks_getunit());
 }
 
+void space_init_sinks_mapper(void *restrict map_data, int sink_count,
+                             void *restrict extra_data) {
+
+  struct sink *restrict sinks = (struct sink *)map_data;
+  for (int k = 0; k < sink_count; k++) sink_init_sink(&sinks[k]);
+}
+
+/**
+ * @brief Calls the #sink initialisation function on all particles in the
+ * space.
+ *
+ * @param s The #space.
+ * @param verbose Are we talkative?
+ */
+void space_init_sinks(struct space *s, int verbose) {
+
+  const ticks tic = getticks();
+
+  if (s->nr_sinks > 0)
+    threadpool_map(&s->e->threadpool, space_init_sinks_mapper, s->sinks,
+                   s->nr_sinks, sizeof(struct sink), threadpool_auto_chunk_size,
+                   /*extra_data=*/NULL);
+  if (verbose)
+    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+            clocks_getunit());
+}
+
 void space_convert_quantities_mapper(void *restrict map_data, int count,
                                      void *restrict extra_data) {
   struct space *s = (struct space *)extra_data;
@@ -4814,10 +5378,12 @@ void space_convert_quantities(struct space *s, int verbose) {
  * @param hydro_properties The properties of the hydro scheme.
  * @param parts Array of Gas particles.
  * @param gparts Array of Gravity particles.
+ * @param sinks Array of sink particles.
  * @param sparts Array of stars particles.
  * @param bparts Array of black hole particles.
  * @param Npart The number of Gas particles in the space.
  * @param Ngpart The number of Gravity particles in the space.
+ * @param Nsink The number of sink particles in the space.
  * @param Nspart The number of stars particles in the space.
  * @param Nbpart The number of black hole particles in the space.
  * @param periodic flag whether the domain is periodic or not.
@@ -4839,8 +5405,8 @@ void space_convert_quantities(struct space *s, int verbose) {
 void space_init(struct space *s, struct swift_params *params,
                 const struct cosmology *cosmo, double dim[3],
                 const struct hydro_props *hydro_properties, struct part *parts,
-                struct gpart *gparts, struct spart *sparts,
-                struct bpart *bparts, size_t Npart, size_t Ngpart,
+                struct gpart *gparts, struct sink *sinks, struct spart *sparts,
+                struct bpart *bparts, size_t Npart, size_t Ngpart, size_t Nsink,
                 size_t Nspart, size_t Nbpart, int periodic, int replicate,
                 int generate_gas_in_ics, int hydro, int self_gravity,
                 int star_formation, int DM_background, int verbose, int dry_run,
@@ -4862,22 +5428,27 @@ void space_init(struct space *s, struct swift_params *params,
   s->nr_gparts = Ngpart;
   s->nr_sparts = Nspart;
   s->nr_bparts = Nbpart;
+  s->nr_sinks = Nsink;
   s->size_parts = Npart;
   s->size_gparts = Ngpart;
   s->size_sparts = Nspart;
   s->size_bparts = Nbpart;
+  s->size_sinks = Nsink;
   s->nr_inhibited_parts = 0;
   s->nr_inhibited_gparts = 0;
   s->nr_inhibited_sparts = 0;
   s->nr_inhibited_bparts = 0;
+  s->nr_inhibited_sinks = 0;
   s->nr_extra_parts = 0;
   s->nr_extra_gparts = 0;
   s->nr_extra_sparts = 0;
   s->nr_extra_bparts = 0;
+  s->nr_extra_sinks = 0;
   s->parts = parts;
   s->gparts = gparts;
   s->sparts = sparts;
   s->bparts = bparts;
+  s->sinks = sinks;
   s->min_part_mass = FLT_MAX;
   s->min_gpart_mass = FLT_MAX;
   s->min_spart_mass = FLT_MAX;
@@ -4902,8 +5473,8 @@ void space_init(struct space *s, struct swift_params *params,
 
 #ifdef SWIFT_DEBUG_CHECKS
     if (!dry_run)
-      part_verify_links(parts, gparts, sparts, bparts, Npart, Ngpart, Nspart,
-                        Nbpart, 1);
+      part_verify_links(parts, gparts, sinks, sparts, bparts, Npart, Ngpart,
+                        Nsink, Nspart, Nbpart, 1);
 #endif
   }
 
@@ -4920,14 +5491,16 @@ void space_init(struct space *s, struct swift_params *params,
     gparts = s->gparts;
     sparts = s->sparts;
     bparts = s->bparts;
+    sinks = s->sinks;
     Npart = s->nr_parts;
     Ngpart = s->nr_gparts;
     Nspart = s->nr_sparts;
     Nbpart = s->nr_bparts;
+    Nsink = s->nr_sinks;
 
 #ifdef SWIFT_DEBUG_CHECKS
-    part_verify_links(parts, gparts, sparts, bparts, Npart, Ngpart, Nspart,
-                      Nbpart, 1);
+    part_verify_links(parts, gparts, sinks, sparts, bparts, Npart, Ngpart,
+                      Nsink, Nspart, Nbpart, 1);
 #endif
   }
 
@@ -4981,6 +5554,11 @@ void space_init(struct space *s, struct swift_params *params,
       params, "Scheduler:cell_extra_gparts", space_extra_gparts_default);
   space_extra_bparts = parser_get_opt_param_int(
       params, "Scheduler:cell_extra_bparts", space_extra_bparts_default);
+  space_extra_sinks = parser_get_opt_param_int(
+      params, "Scheduler:cell_extra_sinks", space_extra_sinks_default);
+  if (space_extra_sinks != 0) {
+    error("Extra sink particles not implemented yet.");
+  }
 
   engine_max_parts_per_ghost =
       parser_get_opt_param_int(params, "Scheduler:engine_max_parts_per_ghost",
@@ -5050,6 +5628,11 @@ void space_init(struct space *s, struct swift_params *params,
       bparts[k].x[1] += shift[1];
       bparts[k].x[2] += shift[2];
     }
+    for (size_t k = 0; k < Nsink; k++) {
+      sinks[k].x[0] += shift[0];
+      sinks[k].x[1] += shift[1];
+      sinks[k].x[2] += shift[2];
+    }
   }
 
   if (!dry_run) {
@@ -5109,6 +5692,20 @@ void space_init(struct space *s, struct swift_params *params,
           if (bparts[k].x[j] < 0 || bparts[k].x[j] >= s->dim[j])
             error("Not all b-particles are within the specified domain.");
     }
+
+    /* Same for the sinks */
+    if (periodic) {
+      for (size_t k = 0; k < Nsink; k++)
+        for (int j = 0; j < 3; j++) {
+          while (sinks[k].x[j] < 0) sinks[k].x[j] += s->dim[j];
+          while (sinks[k].x[j] >= s->dim[j]) sinks[k].x[j] -= s->dim[j];
+        }
+    } else {
+      for (size_t k = 0; k < Nsink; k++)
+        for (int j = 0; j < 3; j++)
+          if (sinks[k].x[j] < 0 || sinks[k].x[j] >= s->dim[j])
+            error("Not all sink-particles are within the specified domain.");
+    }
   }
 
   /* Allocate the extra parts array for the gas particles. */
@@ -5172,18 +5769,21 @@ void space_replicate(struct space *s, int replicate, int verbose) {
   const size_t nr_gparts = s->nr_gparts;
   const size_t nr_sparts = s->nr_sparts;
   const size_t nr_bparts = s->nr_bparts;
+  const size_t nr_sinks = s->nr_sinks;
   const size_t nr_dm = nr_gparts - nr_parts - nr_sparts - nr_bparts;
 
   s->size_parts = s->nr_parts = nr_parts * factor;
   s->size_gparts = s->nr_gparts = nr_gparts * factor;
   s->size_sparts = s->nr_sparts = nr_sparts * factor;
   s->size_bparts = s->nr_bparts = nr_bparts * factor;
+  s->size_sinks = s->nr_sinks = nr_sinks * factor;
 
   /* Allocate space for new particles */
   struct part *parts = NULL;
   struct gpart *gparts = NULL;
   struct spart *sparts = NULL;
   struct bpart *bparts = NULL;
+  struct sink *sinks = NULL;
 
   if (swift_memalign("parts", (void **)&parts, part_align,
                      s->nr_parts * sizeof(struct part)) != 0)
@@ -5197,6 +5797,10 @@ void space_replicate(struct space *s, int replicate, int verbose) {
                      s->nr_sparts * sizeof(struct spart)) != 0)
     error("Failed to allocate new spart array.");
 
+  if (swift_memalign("sinks", (void **)&sinks, sink_align,
+                     s->nr_sinks * sizeof(struct sink)) != 0)
+    error("Failed to allocate new sink array.");
+
   if (swift_memalign("bparts", (void **)&bparts, bpart_align,
                      s->nr_bparts * sizeof(struct bpart)) != 0)
     error("Failed to allocate new bpart array.");
@@ -5216,6 +5820,8 @@ void space_replicate(struct space *s, int replicate, int verbose) {
                nr_bparts * sizeof(struct bpart));
         memcpy(gparts + offset * nr_gparts, s->gparts,
                nr_gparts * sizeof(struct gpart));
+        memcpy(sinks + offset * nr_sinks, s->sinks,
+               nr_sinks * sizeof(struct sink));
 
         /* Shift the positions */
         const double shift[3] = {i * s->dim[0], j * s->dim[1], k * s->dim[2]};
@@ -5240,9 +5846,14 @@ void space_replicate(struct space *s, int replicate, int verbose) {
           bparts[n].x[1] += shift[1];
           bparts[n].x[2] += shift[2];
         }
+        for (size_t n = offset * nr_sinks; n < (offset + 1) * nr_sinks; ++n) {
+          sinks[n].x[0] += shift[0];
+          sinks[n].x[1] += shift[1];
+          sinks[n].x[2] += shift[2];
+        }
 
         /* Set the correct links (recall gpart are sorted by type at start-up):
-           first DM (unassociated gpart), then gas, then stars */
+           first DM (unassociated gpart), then gas, then sinks, then stars */
         if (nr_parts > 0 && nr_gparts > 0) {
           const size_t offset_part = offset * nr_parts;
           const size_t offset_gpart = offset * nr_gparts + nr_dm;
@@ -5252,9 +5863,19 @@ void space_replicate(struct space *s, int replicate, int verbose) {
             gparts[offset_gpart + n].id_or_neg_offset = -(offset_part + n);
           }
         }
+        if (nr_sinks > 0 && nr_gparts > 0) {
+          const size_t offset_sink = offset * nr_sinks;
+          const size_t offset_gpart = offset * nr_gparts + nr_dm + nr_parts;
+
+          for (size_t n = 0; n < nr_sinks; ++n) {
+            sinks[offset_sink + n].gpart = &gparts[offset_gpart + n];
+            gparts[offset_gpart + n].id_or_neg_offset = -(offset_sink + n);
+          }
+        }
         if (nr_sparts > 0 && nr_gparts > 0) {
           const size_t offset_spart = offset * nr_sparts;
-          const size_t offset_gpart = offset * nr_gparts + nr_dm + nr_parts;
+          const size_t offset_gpart =
+              offset * nr_gparts + nr_dm + nr_parts + nr_sinks;
 
           for (size_t n = 0; n < nr_sparts; ++n) {
             sparts[offset_spart + n].gpart = &gparts[offset_gpart + n];
@@ -5264,7 +5885,7 @@ void space_replicate(struct space *s, int replicate, int verbose) {
         if (nr_bparts > 0 && nr_gparts > 0) {
           const size_t offset_bpart = offset * nr_bparts;
           const size_t offset_gpart =
-              offset * nr_gparts + nr_dm + nr_parts + nr_sparts;
+              offset * nr_gparts + nr_dm + nr_parts + nr_sinks + nr_sparts;
 
           for (size_t n = 0; n < nr_bparts; ++n) {
             bparts[offset_bpart + n].gpart = &gparts[offset_gpart + n];
@@ -5280,10 +5901,12 @@ void space_replicate(struct space *s, int replicate, int verbose) {
   swift_free("gparts", s->gparts);
   swift_free("sparts", s->sparts);
   swift_free("bparts", s->bparts);
+  swift_free("sinks", s->sinks);
   s->parts = parts;
   s->gparts = gparts;
   s->sparts = sparts;
   s->bparts = bparts;
+  s->sinks = sinks;
 
   /* Finally, update the domain size */
   s->dim[0] *= replicate;
@@ -5292,8 +5915,9 @@ void space_replicate(struct space *s, int replicate, int verbose) {
 
 #ifdef SWIFT_DEBUG_CHECKS
   /* Verify that everything is correct */
-  part_verify_links(s->parts, s->gparts, s->sparts, s->bparts, s->nr_parts,
-                    s->nr_gparts, s->nr_sparts, s->nr_bparts, verbose);
+  part_verify_links(s->parts, s->gparts, s->sinks, s->sparts, s->bparts,
+                    s->nr_parts, s->nr_gparts, s->nr_sinks, s->nr_sparts,
+                    s->nr_bparts, verbose);
 #endif
 }
 
@@ -5343,6 +5967,9 @@ void space_generate_gas(struct space *s, const struct cosmology *cosmo,
   if (s->nr_bparts != 0)
     error("Generating gas particles from DM but BHs already exists!");
 
+  if (s->nr_sinks != 0)
+    error("Generating gas particles from DM but sink already exists!");
+
   /* Pull out information about particle splitting */
   const int particle_splitting = hydro_properties->particle_splitting;
   const float splitting_mass_threshold =
@@ -5801,6 +6428,7 @@ void space_clean(struct space *s) {
   swift_free("gparts", s->gparts);
   swift_free("sparts", s->sparts);
   swift_free("bparts", s->bparts);
+  swift_free("sinks", s->sinks);
 #ifdef WITH_MPI
   swift_free("parts_foreign", s->parts_foreign);
   swift_free("sparts_foreign", s->sparts_foreign);
@@ -5847,6 +6475,8 @@ void space_struct_dump(struct space *s, FILE *stream) {
                        "space_extra_parts", "space_extra_parts");
   restart_write_blocks(&space_extra_gparts, sizeof(int), 1, stream,
                        "space_extra_gparts", "space_extra_gparts");
+  restart_write_blocks(&space_extra_sinks, sizeof(int), 1, stream,
+                       "space_extra_sinks", "space_extra_sinks");
   restart_write_blocks(&space_extra_sparts, sizeof(int), 1, stream,
                        "space_extra_sparts", "space_extra_sparts");
   restart_write_blocks(&space_extra_bparts, sizeof(int), 1, stream,
@@ -5878,6 +6508,10 @@ void space_struct_dump(struct space *s, FILE *stream) {
     restart_write_blocks(s->gparts, s->nr_gparts, sizeof(struct gpart), stream,
                          "gparts", "gparts");
 
+  if (s->nr_sinks > 0)
+    restart_write_blocks(s->sinks, s->nr_sinks, sizeof(struct sink), stream,
+                         "sinks", "sinks");
+
   if (s->nr_sparts > 0)
     restart_write_blocks(s->sparts, s->nr_sparts, sizeof(struct spart), stream,
                          "sparts", "sparts");
@@ -5920,6 +6554,8 @@ void space_struct_restore(struct space *s, FILE *stream) {
                       "space_extra_parts");
   restart_read_blocks(&space_extra_gparts, sizeof(int), 1, stream, NULL,
                       "space_extra_gparts");
+  restart_read_blocks(&space_extra_sinks, sizeof(int), 1, stream, NULL,
+                      "space_extra_sinks");
   restart_read_blocks(&space_extra_sparts, sizeof(int), 1, stream, NULL,
                       "space_extra_sparts");
   restart_read_blocks(&space_extra_bparts, sizeof(int), 1, stream, NULL,
@@ -5986,6 +6622,16 @@ void space_struct_restore(struct space *s, FILE *stream) {
                         NULL, "gparts");
   }
 
+  s->sinks = NULL;
+  if (s->nr_sinks > 0) {
+    if (swift_memalign("sinks", (void **)&s->sinks, sink_align,
+                       s->size_sinks * sizeof(struct sink)) != 0)
+      error("Failed to allocate restore sink array.");
+
+    restart_read_blocks(s->sinks, s->nr_sinks, sizeof(struct sink), stream,
+                        NULL, "sinks");
+  }
+
   s->sparts = NULL;
   if (s->nr_sparts > 0) {
     if (swift_memalign("sparts", (void **)&s->sparts, spart_align,
@@ -6013,6 +6659,10 @@ void space_struct_restore(struct space *s, FILE *stream) {
   if (s->nr_parts > 0 && s->nr_gparts > 0)
     part_relink_parts_to_gparts(s->gparts, s->nr_gparts, s->parts);
 
+  /* Re-link the sinks. */
+  if (s->nr_sinks > 0 && s->nr_gparts > 0)
+    part_relink_sinks_to_gparts(s->gparts, s->nr_gparts, s->sinks);
+
   /* Re-link the sparts. */
   if (s->nr_sparts > 0 && s->nr_gparts > 0)
     part_relink_sparts_to_gparts(s->gparts, s->nr_gparts, s->sparts);
@@ -6023,8 +6673,9 @@ void space_struct_restore(struct space *s, FILE *stream) {
 
 #ifdef SWIFT_DEBUG_CHECKS
   /* Verify that everything is correct */
-  part_verify_links(s->parts, s->gparts, s->sparts, s->bparts, s->nr_parts,
-                    s->nr_gparts, s->nr_sparts, s->nr_bparts, 1);
+  part_verify_links(s->parts, s->gparts, s->sinks, s->sparts, s->bparts,
+                    s->nr_parts, s->nr_gparts, s->nr_sinks, s->nr_sparts,
+                    s->nr_bparts, 1);
 #endif
 }
 
diff --git a/src/space.h b/src/space.h
index aafc52d97dfe4a9d854c9b4a0b3020002aec4b8c..189d2cee3c4728c3e14f5cd11d50e0a115a67d76 100644
--- a/src/space.h
+++ b/src/space.h
@@ -52,6 +52,7 @@ struct hydro_props;
 #define space_extra_gparts_default 0
 #define space_extra_sparts_default 100
 #define space_extra_bparts_default 0
+#define space_extra_sinks_default 0
 #define space_expected_max_nr_strays_default 100
 #define space_subsize_pair_hydro_default 256000000
 #define space_subsize_self_hydro_default 32000
@@ -282,6 +283,21 @@ struct space {
   /*! Structure dealing with the computation of a unique ID */
   struct unique_id unique_id;
 
+  /*! The total number of #sink in the space. */
+  size_t nr_sinks;
+
+  /*! The total number of #sink we allocated memory for. */
+  size_t size_sinks;
+
+  /*! Number of inhibted sinks in the space */
+  size_t nr_inhibited_sinks;
+
+  /*! Number of extra #sink we allocated (for on-the-fly creation) */
+  size_t nr_extra_sinks;
+
+  /*! The particle data (cells have pointers to this). */
+  struct sink *sinks;
+
 #ifdef WITH_MPI
 
   /*! Buffers for parts that we will receive from foreign cells. */
@@ -308,18 +324,21 @@ void space_free_buff_sort_indices(struct space *s);
 void space_parts_sort(struct part *parts, struct xpart *xparts, int *ind,
                       int *counts, int num_bins, ptrdiff_t parts_offset);
 void space_gparts_sort(struct gpart *gparts, struct part *parts,
-                       struct spart *sparts, struct bpart *bparts, int *ind,
-                       int *counts, int num_bins);
+                       struct sink *sinks, struct spart *sparts,
+                       struct bpart *bparts, int *ind, int *counts,
+                       int num_bins);
 void space_sparts_sort(struct spart *sparts, int *ind, int *counts,
                        int num_bins, ptrdiff_t sparts_offset);
 void space_bparts_sort(struct bpart *bparts, int *ind, int *counts,
                        int num_bins, ptrdiff_t bparts_offset);
+void space_sinks_sort(struct sink *sinks, int *ind, int *counts, int num_bins,
+                      ptrdiff_t sinks_offset);
 void space_getcells(struct space *s, int nr_cells, struct cell **cells);
 void space_init(struct space *s, struct swift_params *params,
                 const struct cosmology *cosmo, double dim[3],
                 const struct hydro_props *hydro_properties, struct part *parts,
-                struct gpart *gparts, struct spart *sparts,
-                struct bpart *bparts, size_t Npart, size_t Ngpart,
+                struct gpart *gparts, struct sink *sinks, struct spart *sparts,
+                struct bpart *bparts, size_t Npart, size_t Ngpart, size_t Nsink,
                 size_t Nspart, size_t Nbpart, int periodic, int replicate,
                 int generate_gas_in_ics, int hydro, int gravity,
                 int star_formation, int DM_background, int verbose, int dry_run,
@@ -357,15 +376,20 @@ void space_sparts_get_cell_index(struct space *s, int *sind, int *cell_counts,
 void space_bparts_get_cell_index(struct space *s, int *sind, int *cell_counts,
                                  size_t *count_inhibited_bparts,
                                  size_t *count_extra_bparts, int verbose);
+void space_sinks_get_cell_index(struct space *s, int *sind, int *cell_counts,
+                                size_t *count_inhibited_sinks,
+                                size_t *count_extra_sinks, int verbose);
 void space_synchronize_particle_positions(struct space *s);
 void space_first_init_parts(struct space *s, int verbose);
 void space_first_init_gparts(struct space *s, int verbose);
 void space_first_init_sparts(struct space *s, int verbose);
 void space_first_init_bparts(struct space *s, int verbose);
+void space_first_init_sinks(struct space *s, int verbose);
 void space_init_parts(struct space *s, int verbose);
 void space_init_gparts(struct space *s, int verbose);
 void space_init_sparts(struct space *s, int verbose);
 void space_init_bparts(struct space *s, int verbose);
+void space_init_sinks(struct space *s, int verbose);
 void space_convert_quantities(struct space *s, int verbose);
 void space_link_cleanup(struct space *s);
 void space_check_drift_point(struct space *s, integertime_t ti_drift,
diff --git a/tests/testReading.c b/tests/testReading.c
index f1013d50dfe951eb93591a8767e0f13e6d2c04ea..550ad1a1a0e266fb91fccdbcaa7d913ec73dbdd7 100644
--- a/tests/testReading.c
+++ b/tests/testReading.c
@@ -28,7 +28,8 @@
 
 int main(int argc, char *argv[]) {
 
-  size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0, Nspart = 0, Nbpart = 0;
+  size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0, Nspart = 0, Nbpart = 0,
+         Nsink = 0;
   int flag_entropy_ICs = -1;
   int i, j, k;
   double dim[3];
@@ -36,6 +37,7 @@ int main(int argc, char *argv[]) {
   struct gpart *gparts = NULL;
   struct spart *sparts = NULL;
   struct bpart *bparts = NULL;
+  struct sink *sinks = NULL;
 
   /* Default unit system */
   struct unit_system us;
@@ -50,11 +52,12 @@ int main(int argc, char *argv[]) {
 #endif
 
   /* Read data */
-  read_ic_single("input.hdf5", &us, dim, &parts, &gparts, &sparts, &bparts,
-                 &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
-                 &flag_entropy_ICs,
+  read_ic_single("input.hdf5", &us, dim, &parts, &gparts, &sinks, &sparts,
+                 &bparts, &Ngas, &Ngpart, &Ngpart_background, &Nsink, &Nspart,
+                 &Nbpart, &flag_entropy_ICs,
                  /*with_hydro=*/1,
                  /*with_gravity=*/1,
+                 /*with_sink=*/0,
                  /*with_stars=*/0,
                  /*with_black_holes=*/0,
                  /*with_cosmology=*/0,
diff --git a/tests/testSelectOutput.c b/tests/testSelectOutput.c
index 063d788e7f00ed6f01a5d793b84b47ea74bb8f6d..a4354493bff6683cb51090ca936e19872b53e42e 100644
--- a/tests/testSelectOutput.c
+++ b/tests/testSelectOutput.c
@@ -86,7 +86,8 @@ int main(int argc, char *argv[]) {
   clocks_set_cpufreq(cpufreq);
 
   // const char *base_name = "testSelectOutput";
-  size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0, Nspart = 0, Nbpart = 0;
+  size_t Ngas = 0, Ngpart = 0, Ngpart_background = 0, Nspart = 0, Nbpart = 0,
+         Nsink = 0;
   int flag_entropy_ICs = -1;
   int periodic = 1;
   double dim[3];
@@ -94,6 +95,7 @@ int main(int argc, char *argv[]) {
   struct gpart *gparts = NULL;
   struct spart *sparts = NULL;
   struct bpart *bparts = NULL;
+  struct sink *sinks = NULL;
 
   /* parse parameters */
   message("Reading parameters.");
@@ -116,11 +118,12 @@ int main(int argc, char *argv[]) {
 
   /* Read data */
   message("Reading initial conditions.");
-  read_ic_single("input.hdf5", &us, dim, &parts, &gparts, &sparts, &bparts,
-                 &Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
-                 &flag_entropy_ICs,
+  read_ic_single("input.hdf5", &us, dim, &parts, &gparts, &sinks, &sparts,
+                 &bparts, &Ngas, &Ngpart, &Ngpart_background, &Nsink, &Nspart,
+                 &Nbpart, &flag_entropy_ICs,
                  /*with_hydro=*/1,
                  /*with_gravity=*/0,
+                 /*with_sink=*/0,
                  /*with_stars=*/0,
                  /*with_black_holes=*/0,
                  /*with_cosmology=*/0,