diff --git a/configure.ac b/configure.ac
index 3c5e43b579e8cd6065d962294f535467f8fd1441..3d5ab71e707a414eecd36894c4d9db0c9dfe33ef 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1954,7 +1954,7 @@ esac
 #  chemistry function
 AC_ARG_WITH([chemistry],
    [AS_HELP_STRING([--with-chemistry=<function>],
-      [chemistry function @<:@none, GEAR_*, QLA, EAGLE default: none@:>@
+      [chemistry function @<:@none, GEAR_*, GEAR_DIFFUSION_*, QLA, EAGLE default: none@:>@
       For GEAR, you need to provide the number of elements (e.g. GEAR_10)]
    )],
    [with_chemistry="$withval"],
@@ -1973,6 +1973,11 @@ case "$with_chemistry" in
    none)
       AC_DEFINE([CHEMISTRY_NONE], [1], [No chemistry function])
    ;;
+   GEAR_DIFFUSION_*)
+      AC_DEFINE([CHEMISTRY_GEAR_DIFFUSION], [1], [Chemistry taken from the GEAR model including the metal diffusion])
+      number_element=${with_chemistry:15}
+      AC_DEFINE_UNQUOTED([GEAR_CHEMISTRY_ELEMENT_COUNT], [$number_element], [Number of element to follow])
+   ;;
    GEAR_*)
       AC_DEFINE([CHEMISTRY_GEAR], [1], [Chemistry taken from the GEAR model])
       number_element=${with_chemistry:5}
diff --git a/examples/SubgridTests/SmoothedMetallicity/makeIC.py b/examples/SubgridTests/SmoothedMetallicity/makeIC.py
index 542b4c5911c942015d16595f42e73ca8978d20da..ce5a2b857763aff6bf690995353d98f84da53dbc 100644
--- a/examples/SubgridTests/SmoothedMetallicity/makeIC.py
+++ b/examples/SubgridTests/SmoothedMetallicity/makeIC.py
@@ -1,40 +1,45 @@
 ###############################################################################
 # 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
 
-# Generates a swift IC file for the Smoothed Metallicity test in a periodic cubic box
+# Generates a swift IC file for the Smoothed Metallicity test
+# in a periodic cubic box
 
 # Parameters
 gamma = 5./3.      # Gas adiabatic index
 rho0 = 1.          # Background density
-P0 = 1.e-6         # Background pressure
-Nelem = 9          # Gear: 9, EAGLE: 9
+P0 = 1e-6          # Background pressure
+Nelem = 10         # Gear: 10, EAGLE: 9
 low_metal = -6     # Low iron fraction
-high_metal = -5    # high iron fraction
-sigma_metal = 0.1  # relative standard deviation for the metallicities
+high_metal = -5.5  # high iron fraction
+max_shift = 1      # Shift between the different elements
+sigma_metal = 0.2  # relative standard deviation for the metallicities
 fileName = "smoothed_metallicity.hdf5"
 
 # shift all metals in order to obtain nicer plots
-low_metal = [low_metal] * Nelem + np.linspace(0, 3, Nelem)
-high_metal = [high_metal] * Nelem + np.linspace(0, 3, Nelem)
+low_metal = [low_metal] * Nelem + np.linspace(0, max_shift, Nelem)
+low_metal = 10**low_metal
+
+high_metal = [high_metal] * Nelem + np.linspace(0, max_shift, Nelem)
+high_metal = 10**high_metal
 
 # ---------------------------------------------------
 glass = h5py.File("glassCube_32.hdf5", "r")
@@ -52,21 +57,24 @@ ids = np.linspace(1, numPart, numPart)
 m = np.zeros(numPart)
 u = np.zeros(numPart)
 r = np.zeros(numPart)
-Z = np.zeros((numPart, Nelem))
+mass_frac = np.zeros((numPart, Nelem))
 
-r = np.sqrt((pos[:, 0] - 0.5)**2 + (pos[:, 1] - 0.5)**2 + (pos[:, 2] - 0.5)**2)
 m[:] = rho0 * vol / numPart
 u[:] = P0 / (rho0 * (gamma - 1))
 
 # set metallicities
 select = pos[:, 0] < 0.5
 nber = sum(select)
-Z[select, :] = low_metal * (1 + np.random.normal(loc=0., scale=sigma_metal,
-                                                 size=(nber, Nelem)))
+mass_frac[select, :] = low_metal * (
+    1 + np.random.normal(loc=0., scale=sigma_metal,
+                         size=(nber, Nelem)))
 nber = numPart - nber
-Z[np.logical_not(select), :] = high_metal * (1 + np.random.normal(
+mass_frac[~select, :] = high_metal * (1 + np.random.normal(
     loc=0., scale=sigma_metal, size=(nber, Nelem)))
 
+v[select, 2] = 1
+v[~select, 2] = -1
+
 # --------------------------------------------------
 
 # File
@@ -100,7 +108,8 @@ grp.create_dataset('Masses', data=m, dtype='f')
 grp.create_dataset('SmoothingLength', data=h, dtype='f')
 grp.create_dataset('InternalEnergy', data=u, dtype='f')
 grp.create_dataset('ParticleIDs', data=ids, dtype='L')
-grp.create_dataset('ElementAbundance', data=Z, dtype='f')
+grp.create_dataset('MetalMassFraction', data=mass_frac, dtype='d')
+grp.create_dataset('ElementAbundance', data=mass_frac, dtype='d')
 
 
 file.close()
diff --git a/examples/SubgridTests/SmoothedMetallicity/plotSolution.py b/examples/SubgridTests/SmoothedMetallicity/plotSolution.py
index 068fe5378e19c34ee8a68398f4e0ed096d0982e0..ab5e18d7d9d9bb15b154693e61c6fb93bb8bcc79 100644
--- a/examples/SubgridTests/SmoothedMetallicity/plotSolution.py
+++ b/examples/SubgridTests/SmoothedMetallicity/plotSolution.py
@@ -30,14 +30,17 @@ matplotlib.use("Agg")
 import matplotlib.pyplot as plt
 
 # Parameters
-low_metal = -6  # low metal abundance
-high_metal = -5  # High metal abundance
-sigma_metal = 0.1  # relative standard deviation for Z
+low_metal = -6     # low metal abundance
+high_metal = -5.5  # High metal abundance
+sigma_metal = 0.2  # relative standard deviation for Z
+max_shift = 1      # Shift between the different elements
 
-Nelem = 9
+Nelem = 10
 # shift all metals in order to obtain nicer plots
-low_metal = [low_metal] * Nelem + np.linspace(0, 3, Nelem)
-high_metal = [high_metal] * Nelem + np.linspace(0, 3, Nelem)
+low_metal = [low_metal] * Nelem + np.linspace(0, max_shift, Nelem)
+low_metal = 10**low_metal
+high_metal = [high_metal] * Nelem + np.linspace(0, max_shift, Nelem)
+high_metal = 10**high_metal
 
 # ---------------------------------------------------------------
 # Don't touch anything after this.
@@ -80,12 +83,20 @@ kernel = sim["/HydroScheme"].attrs["Kernel function"]
 neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
 eta = sim["/HydroScheme"].attrs["Kernel eta"]
 chemistry = sim["/SubgridScheme"].attrs["Chemistry Model"]
+chemistry = str(chemistry)
 git = sim["Code"].attrs["Git Revision"]
 
 pos = sim["/PartType0/Coordinates"][:, :]
 d = pos[:, 0] - boxSize / 2
-smooth_metal = sim["/PartType0/SmoothedElementMassFractions"][:, :]
-metal = sim["/PartType0/ElementMassFractions"][:, :]
+gear_name = "/PartType0/MetalMassFractions"
+eagle_name = "/PartType0/ElementAbundances"
+if gear_name in sim:
+    smooth_metal = sim["/PartType0/SmoothedMetalMassFractions"][:, :]
+    metal = sim[gear_name][:, :]
+else:
+    smooth_metal = sim["/PartType0/SmoothedElementAbundances"][:, :]
+    metal = sim[eagle_name][:, :]
+
 h = sim["/PartType0/SmoothingLengths"][:]
 h = np.mean(h)
 
@@ -165,15 +176,15 @@ plt.figure()
 # Metallicity --------------------------------
 plt.subplot(221)
 for e in range(Nelem):
-    plt.plot(metal[:, e], smooth_metal[:, e], ".", ms=0.5, alpha=0.2)
+    plt.loglog(metal[:, e], smooth_metal[:, e], ".", ms=0.5, alpha=0.2)
 
 xmin, xmax = metal.min(), metal.max()
 ymin, ymax = smooth_metal.min(), smooth_metal.max()
 x = max(xmin, ymin)
 y = min(xmax, ymax)
-plt.plot([x, y], [x, y], "--k", lw=1.0)
-plt.xlabel("${\\rm{Metallicity}}~Z_\\textrm{part}$", labelpad=0)
-plt.ylabel("${\\rm{Smoothed~Metallicity}}~Z_\\textrm{sm}$", labelpad=0)
+plt.loglog([x, y], [x, y], "--k", lw=1.0)
+plt.xlabel("${\\rm{Metal~mass~fraction}}$", labelpad=0)
+plt.ylabel("${\\rm{Smoothed~metal~mass~fraction}}$", labelpad=0)
 
 # Metallicity --------------------------------
 e = 0
@@ -183,10 +194,11 @@ plt.plot(d_a, sol[:, e], "--", color="b", alpha=0.8, lw=1.2)
 plt.fill_between(
     d_a, sig[:, e, 0], sig[:, e, 1], facecolor="b", interpolate=True, alpha=0.5
 )
+plt.yscale("log")
 plt.xlabel("${\\rm{Distance}}~r$", labelpad=0)
-plt.ylabel("${\\rm{Smoothed~Metallicity}}~Z_\\textrm{sm}$", labelpad=0)
+plt.ylabel("${\\rm{Smoothed~metal~mass~fraction}}$", labelpad=0)
 plt.xlim(-0.5, 0.5)
-plt.ylim(low_metal[e] - 1, high_metal[e] + 1)
+# plt.ylim(low_metal[e] - 1, high_metal[e] + 1)
 
 # Information -------------------------------------
 plt.subplot(222, frameon=False)
diff --git a/examples/SubgridTests/SmoothedMetallicity/smoothed_metallicity.yml b/examples/SubgridTests/SmoothedMetallicity/smoothed_metallicity.yml
index f6841c6bd0744b4bbeacbe136a126b4ed5631f6f..aeb725235f59f9b0d2c2126d3af26ab3ec7e50da 100644
--- a/examples/SubgridTests/SmoothedMetallicity/smoothed_metallicity.yml
+++ b/examples/SubgridTests/SmoothedMetallicity/smoothed_metallicity.yml
@@ -33,3 +33,7 @@ InitialConditions:
   file_name:  ./smoothed_metallicity.hdf5          # The file to read
   periodic:   1
 
+GEARChemistry:
+  scale_initial_metallicity: 0 # Scale the initial metallicity with solar abundances
+  initial_metallicity: -1 # Read from the snapshot
+  diffusion_coefficient: 1e-3 # Coefficient for the diffusion (see Shen et al. 2010)
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index d89758b4af0a0456f324b0fc682609cb47f8a839..dd54e48be256ae095fd994e9c706b438f384f4b1 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -455,6 +455,7 @@ EAGLEChemistry:
 GEARChemistry:
   initial_metallicity: 1         # Initial metallicity of the gas (mass fraction)
   scale_initial_metallicity: 1   # Should we scale the initial metallicity with the solar one?
+  diffusion_coefficient: 1e-3    # Coefficient for the diffusion (see Shen et al. 2010; differs by gamma^2 due to h^2).
 
 
 # Parameters related to star formation models  -----------------------------------------------
diff --git a/src/Makefile.am b/src/Makefile.am
index f1558a1504b352752f57c1b1946e3ab380ace12d..e52c63e60c93a8cb7ac297bcec7487e9ca4249d8 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -139,7 +139,7 @@ nobase_noinst_HEADERS += gravity_iact.h kernel_long_gravity.h vector.h accumulat
 nobase_noinst_HEADERS += runner_doiact_nosort.h runner_doiact_hydro.h runner_doiact_stars.h runner_doiact_black_holes.h runner_doiact_grav.h 
 nobase_noinst_HEADERS += runner_doiact_functions_hydro.h runner_doiact_functions_stars.h runner_doiact_functions_black_holes.h 
 nobase_noinst_HEADERS += runner_doiact_functions_limiter.h runner_doiact_limiter.h units.h intrinsics.h minmax.h 
-nobase_noinst_HEADERS += runner_doiact_rt.h runner_doiact_functions_rt.h runner_doiact_sinks.h
+nobase_noinst_HEADERS += runner_doiact_rt.h runner_doiact_functions_rt.h runner_doiact_sinks.h runner_doiact_functions_sinks.h
 nobase_noinst_HEADERS += kick.h timestep.h drift.h adiabatic_index.h io_properties.h dimension.h part_type.h periodic.h memswap.h 
 nobase_noinst_HEADERS += timestep_limiter.h timestep_limiter_iact.h timestep_sync.h timestep_sync_part.h timestep_limiter_struct.h 
 nobase_noinst_HEADERS += dump.h logger.h sign.h logger_io.h hashmap.h gravity.h gravity_io.h gravity_logger.h  gravity_cache.h output_options.h
@@ -324,8 +324,11 @@ nobase_noinst_HEADERS += pressure_floor/GEAR/pressure_floor_iact.h pressure_floo
 nobase_noinst_HEADERS += pressure_floor/GEAR/pressure_floor_struct.h pressure_floor/none/pressure_floor_struct.h 
 nobase_noinst_HEADERS += sink/Default/sink.h sink/Default/sink_io.h sink/Default/sink_part.h sink/Default/sink_properties.h 
 nobase_noinst_HEADERS += sink/Default/sink_iact.h
-nobase_noinst_HEADERS += sink.h sink_io.h sink_properties.h runner_doiact_functions_sinks.h
-
+nobase_noinst_HEADERS += sink.h sink_io.h sink_properties.h
+nobase_noinst_HEADERS += chemistry/GEAR_DIFFUSION/chemistry.h
+nobase_noinst_HEADERS += chemistry/GEAR_DIFFUSION/chemistry_io.h
+nobase_noinst_HEADERS += chemistry/GEAR_DIFFUSION/chemistry_struct.h
+nobase_noinst_HEADERS += chemistry/GEAR_DIFFUSION/chemistry_iact.h
 
 # Sources and special flags for the gravity library
 libgrav_la_SOURCES = runner_doiact_grav.c
diff --git a/src/chemistry.h b/src/chemistry.h
index c6113cbbc64e3d1dcbd4c0cea0fdef25b1aaa75f..45db33d1f746294af017c123d7cbf02ac818b995 100644
--- a/src/chemistry.h
+++ b/src/chemistry.h
@@ -35,6 +35,9 @@
 #elif defined(CHEMISTRY_GEAR)
 #include "./chemistry/GEAR/chemistry.h"
 #include "./chemistry/GEAR/chemistry_iact.h"
+#elif defined(CHEMISTRY_GEAR_DIFFUSION)
+#include "./chemistry/GEAR_DIFFUSION/chemistry.h"
+#include "./chemistry/GEAR_DIFFUSION/chemistry_iact.h"
 #elif defined(CHEMISTRY_QLA)
 #include "./chemistry/QLA/chemistry.h"
 #include "./chemistry/QLA/chemistry_iact.h"
diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h
index cab8916bce17fb5beda5f01ef679bf2fdbdaf0ad..f4e109eda5ac14ad8b5e496d76a3f89cb294c89e 100644
--- a/src/chemistry/EAGLE/chemistry.h
+++ b/src/chemistry/EAGLE/chemistry.h
@@ -294,9 +294,13 @@ static INLINE void chemistry_print_backend(
  *
  * @param p The particle to act upon.
  * @param cosmo The current cosmological model.
+ * @param with_cosmology Are we running with the cosmology?
+ * @param time Current time of the simulation.
+ * @param dt Time step (in physical units).
  */
 __attribute__((always_inline)) INLINE static void chemistry_end_force(
-    struct part* restrict p, const struct cosmology* cosmo) {}
+    struct part* restrict p, const struct cosmology* cosmo,
+    const int with_cosmology, const double time, const double dt) {}
 
 /**
  * @brief Computes the chemistry-related time-step constraint.
diff --git a/src/chemistry/EAGLE/chemistry_iact.h b/src/chemistry/EAGLE/chemistry_iact.h
index 573291637d66e39a973e662edec084d3a4687050..a6f0be62c7182162bd8127aa8ec5e56e25922c76 100644
--- a/src/chemistry/EAGLE/chemistry_iact.h
+++ b/src/chemistry/EAGLE/chemistry_iact.h
@@ -130,4 +130,54 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
       mj * chj->iron_mass_fraction_from_SNIa * wi;
 }
 
+/**
+ * @brief do metal diffusion computation in the <FORCE LOOP>
+ * (symmetric version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ * @param time_base The time base used in order to convert integer to float
+ * time.
+ * @param ti_current The current time (in integer)
+ * @param cosmo The #cosmology.
+ * @param with_cosmology Are we running with cosmology?
+ *
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_diffusion(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H, float time_base,
+    integertime_t t_current, const struct cosmology *cosmo,
+    const int with_cosmology) {}
+
+/**
+ * @brief do metal diffusion computation in the <FORCE LOOP>
+ * (nonsymmetric version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ * @param time_base The time base used in order to convert integer to float
+ * time.
+ * @param ti_current The current time (in integer)
+ * @param cosmo The #cosmology.
+ * @param with_cosmology Are we running with cosmology?
+ *
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_diffusion(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H, float time_base,
+    integertime_t t_current, const struct cosmology *cosmo,
+    const int with_cosmology) {}
+
 #endif /* SWIFT_EAGLE_CHEMISTRY_IACT_H */
diff --git a/src/chemistry/GEAR/chemistry.h b/src/chemistry/GEAR/chemistry.h
index 754d1636ddf9ecf30931e3ecf2db7ad3025822b0..9e47c0edb926e19a6de708b02a12a77e668bce7c 100644
--- a/src/chemistry/GEAR/chemistry.h
+++ b/src/chemistry/GEAR/chemistry.h
@@ -145,6 +145,12 @@ static INLINE void chemistry_init_backend(struct swift_params* parameter_file,
   const float initial_metallicity = parser_get_param_float(
       parameter_file, "GEARChemistry:initial_metallicity");
 
+  if (initial_metallicity < 0) {
+    message("Setting the initial metallicity from the snapshot.");
+  } else {
+    message("Setting the initial metallicity from the parameter file.");
+  }
+
   /* Set the initial metallicities */
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
     data->initial_metallicities[i] = initial_metallicity;
@@ -177,10 +183,6 @@ __attribute__((always_inline)) INLINE static void chemistry_init_part(
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
     /* Reset the smoothed metallicity */
     cpd->smoothed_metal_mass_fraction[i] = 0.f;
-
-    /* Convert the total mass into mass fraction */
-    /* Now the metal mass is not available anymore */
-    cpd->metal_mass_fraction[i] = cpd->metal_mass[i] / p->mass;
   }
 }
 
@@ -204,21 +206,15 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density(
   const float h = p->h;
   const float h_inv = 1.0f / h;                       /* 1/h */
   const float factor = pow_dimension(h_inv) / p->rho; /* 1 / h^d * rho */
-  const float m = p->mass;
 
   struct chemistry_part_data* cpd = &p->chemistry_data;
 
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
     /* Final operation on the density (add self-contribution). */
-    cpd->smoothed_metal_mass_fraction[i] +=
-        m * cpd->metal_mass_fraction[i] * kernel_root;
+    cpd->smoothed_metal_mass_fraction[i] += cpd->metal_mass[i] * kernel_root;
 
     /* Finish the calculation by inserting the missing h-factors */
     cpd->smoothed_metal_mass_fraction[i] *= factor;
-
-    /* Convert the mass fraction into a total mass */
-    /* Now the metal mass fraction is not available anymore */
-    cpd->metal_mass[i] = m * cpd->metal_mass_fraction[i];
   }
 }
 
@@ -229,7 +225,8 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density(
  * @param cosmo The current cosmological model.
  */
 __attribute__((always_inline)) INLINE static void chemistry_end_force(
-    struct part* restrict p, const struct cosmology* cosmo) {}
+    struct part* restrict p, const struct cosmology* cosmo,
+    const int with_cosmology, const double time, const double dt) {}
 
 /**
  * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
@@ -248,9 +245,7 @@ chemistry_part_has_no_neighbours(struct part* restrict p,
   /* Set the smoothed fractions with the non smoothed fractions */
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
     p->chemistry_data.smoothed_metal_mass_fraction[i] =
-        p->chemistry_data.metal_mass_fraction[i];
-    p->chemistry_data.metal_mass[i] =
-        p->chemistry_data.metal_mass_fraction[i] * p->mass;
+        p->chemistry_data.metal_mass[i] / hydro_get_mass(p);
   }
 }
 
@@ -294,7 +289,14 @@ __attribute__((always_inline)) INLINE static void chemistry_first_init_part(
     struct xpart* restrict xp) {
 
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
-    p->chemistry_data.metal_mass[i] = data->initial_metallicities[i] * p->mass;
+    if (data->initial_metallicities[i] < 0) {
+      /* Use the value from the IC. We are reading the metal mass fraction. */
+      p->chemistry_data.metal_mass[i] *= hydro_get_mass(p);
+    } else {
+      /* Use the value from the parameter file */
+      p->chemistry_data.metal_mass[i] =
+          data->initial_metallicities[i] * hydro_get_mass(p);
+    }
   }
 
   chemistry_init_part(p, data);
@@ -405,7 +407,7 @@ __attribute__((always_inline)) INLINE static void chemistry_split_part(
  */
 __attribute__((always_inline)) INLINE static float const*
 chemistry_get_metal_mass_fraction_for_feedback(const struct part* restrict p) {
-
+  error("Not implemented");
   return NULL;
 }
 
@@ -420,6 +422,7 @@ chemistry_get_metal_mass_fraction_for_feedback(const struct part* restrict p) {
 __attribute__((always_inline)) INLINE static float
 chemistry_get_total_metal_mass_fraction_for_feedback(
     const struct part* restrict p) {
+  error("Not implemented");
   return 0.f;
 }
 
@@ -537,7 +540,7 @@ chemistry_get_star_total_metal_mass_for_stats(const struct spart* restrict sp) {
  */
 __attribute__((always_inline)) INLINE static float
 chemistry_get_bh_total_metal_mass_for_stats(const struct bpart* restrict bp) {
-
+  error("Not implemented");
   return 0.f;
 }
 
diff --git a/src/chemistry/GEAR/chemistry_iact.h b/src/chemistry/GEAR/chemistry_iact.h
index ff49af4abb141ee69dea48bbaf99f6b38db4ce45..e2ea733d44c0855ebe2bdc38f85a1fab2fa6ce45 100644
--- a/src/chemistry/GEAR/chemistry_iact.h
+++ b/src/chemistry/GEAR/chemistry_iact.h
@@ -49,30 +49,24 @@ __attribute__((always_inline)) INLINE static void runner_iact_chemistry(
   struct chemistry_part_data *chi = &pi->chemistry_data;
   struct chemistry_part_data *chj = &pj->chemistry_data;
 
-  float wi, wi_dx;
-  float wj, wj_dx;
-
-  /* Get the masses. */
-  const float mi = pi->mass;
-  const float mj = pj->mass;
+  float wi;
+  float wj;
 
   /* Get r */
   const float r = sqrtf(r2);
 
   /* Compute the kernel function for pi */
   const float ui = r / hi;
-  kernel_deval(ui, &wi, &wi_dx);
+  kernel_eval(ui, &wi);
 
   /* Compute the kernel function for pj */
   const float uj = r / hj;
-  kernel_deval(uj, &wj, &wj_dx);
+  kernel_eval(uj, &wj);
 
   /* Compute contribution to the smooth metallicity */
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
-    chi->smoothed_metal_mass_fraction[i] +=
-        mj * chj->metal_mass_fraction[i] * wi;
-    chj->smoothed_metal_mass_fraction[i] +=
-        mi * chi->metal_mass_fraction[i] * wj;
+    chi->smoothed_metal_mass_fraction[i] += chj->metal_mass[i] * wi;
+    chj->smoothed_metal_mass_fraction[i] += chi->metal_mass[i] * wj;
   }
 }
 
@@ -96,23 +90,69 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
   struct chemistry_part_data *chi = &pi->chemistry_data;
   const struct chemistry_part_data *chj = &pj->chemistry_data;
 
-  float wi, wi_dx;
-
-  /* Get the masses. */
-  const float mj = pj->mass;
+  float wi;
 
   /* Get r */
   const float r = sqrtf(r2);
 
   /* Compute the kernel function for pi */
   const float ui = r / hi;
-  kernel_deval(ui, &wi, &wi_dx);
+  kernel_eval(ui, &wi);
 
   /* Compute contribution to the smooth metallicity */
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
-    chi->smoothed_metal_mass_fraction[i] +=
-        mj * chj->metal_mass_fraction[i] * wi;
+    chi->smoothed_metal_mass_fraction[i] += chj->metal_mass[i] * wi;
   }
 }
 
+/**
+ * @brief do metal diffusion computation in the <FORCE LOOP>
+ * (symmetric version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ * @param time_base The time base used in order to convert integer to float
+ * time.
+ * @param ti_current The current time (in integer)
+ * @param cosmo The #cosmology.
+ * @param with_cosmology Are we running with cosmology?
+ *
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_diffusion(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H, float time_base,
+    integertime_t t_current, const struct cosmology *cosmo,
+    const int with_cosmology) {}
+
+/**
+ * @brief do metal diffusion computation in the <FORCE LOOP>
+ * (nonsymmetric version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ * @param time_base The time base used in order to convert integer to float
+ * time.
+ * @param ti_current The current time (in integer)
+ * @param cosmo The #cosmology.
+ * @param with_cosmology Are we running with cosmology?
+ *
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_diffusion(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H, float time_base,
+    integertime_t t_current, const struct cosmology *cosmo,
+    const int with_cosmology) {}
+
 #endif /* SWIFT_GEAR_CHEMISTRY_IACT_H */
diff --git a/src/chemistry/GEAR/chemistry_io.h b/src/chemistry/GEAR/chemistry_io.h
index 8112bfb13b3203d75905167f3e75c8b06f05fb83..b485329a36733226113343b942b50f28578dc0ad 100644
--- a/src/chemistry/GEAR/chemistry_io.h
+++ b/src/chemistry/GEAR/chemistry_io.h
@@ -42,8 +42,8 @@ INLINE static int chemistry_read_particles(struct part* parts,
 
   /* List what we want to read */
   list[0] = io_make_input_field(
-      "ElementAbundance", DOUBLE, GEAR_CHEMISTRY_ELEMENT_COUNT, OPTIONAL,
-      UNIT_CONV_NO_UNITS, parts, chemistry_data.metal_mass_fraction);
+      "MetalMassFraction", DOUBLE, GEAR_CHEMISTRY_ELEMENT_COUNT, OPTIONAL,
+      UNIT_CONV_NO_UNITS, parts, chemistry_data.metal_mass);
 
   return 1;
 }
@@ -53,7 +53,7 @@ INLINE static void convert_gas_metals(const struct engine* e,
                                       const struct xpart* xp, double* ret) {
 
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
-    ret[0] = p->chemistry_data.metal_mass[i] / hydro_get_mass(p);
+    ret[i] = p->chemistry_data.metal_mass[i] / hydro_get_mass(p);
   }
 }
 
@@ -74,13 +74,13 @@ INLINE static int chemistry_write_particles(const struct part* parts,
 
   /* List what we want to write */
   list[0] = io_make_output_field(
-      "SmoothedElementAbundances", DOUBLE, GEAR_CHEMISTRY_ELEMENT_COUNT,
+      "SmoothedMetalMassFractions", DOUBLE, GEAR_CHEMISTRY_ELEMENT_COUNT,
       UNIT_CONV_NO_UNITS, 0.f, parts,
       chemistry_data.smoothed_metal_mass_fraction,
-      "Element abundances smoothed over the neighbors");
+      "Mass fraction of each element smoothed over the neighbors");
 
   list[1] = io_make_output_field_convert_part(
-      "ElementAbundances", DOUBLE, GEAR_CHEMISTRY_ELEMENT_COUNT,
+      "MetalMassFractions", DOUBLE, GEAR_CHEMISTRY_ELEMENT_COUNT,
       UNIT_CONV_NO_UNITS, 0.f, parts, xparts, convert_gas_metals,
       "Mass fraction of each element");
 
@@ -100,7 +100,7 @@ INLINE static int chemistry_write_sparticles(const struct spart* sparts,
 
   /* List what we want to write */
   list[0] = io_make_output_field(
-      "ElementAbundances", DOUBLE, GEAR_CHEMISTRY_ELEMENT_COUNT,
+      "MetalMassFractions", DOUBLE, GEAR_CHEMISTRY_ELEMENT_COUNT,
       UNIT_CONV_NO_UNITS, 0.f, sparts, chemistry_data.metal_mass_fraction,
       "Mass fraction of each element");
 
diff --git a/src/chemistry/GEAR/chemistry_struct.h b/src/chemistry/GEAR/chemistry_struct.h
index 1332deea87e74a874c3b0a76339d88f8edcef13a..9ad831b27d84409c131fd5aa8df1a4204dbb73ef 100644
--- a/src/chemistry/GEAR/chemistry_struct.h
+++ b/src/chemistry/GEAR/chemistry_struct.h
@@ -24,7 +24,7 @@
  */
 struct chemistry_global_data {
 
-  /* Initial metallicity Z */
+  /* Initial mass fraction */
   double initial_metallicities[GEAR_CHEMISTRY_ELEMENT_COUNT];
 };
 
@@ -33,15 +33,8 @@ struct chemistry_global_data {
  */
 struct chemistry_part_data {
 
-  union {
-    /*! Fraction of the particle mass in a given element.
-      This field is available only during the density hydro loop. */
-    double metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT];
-
-    /*! Total mass of element in a particle.
-      This field is available only outside the density hydro loop. */
-    double metal_mass[GEAR_CHEMISTRY_ELEMENT_COUNT];
-  };
+  /*! Total mass of element in a particle. */
+  double metal_mass[GEAR_CHEMISTRY_ELEMENT_COUNT];
 
   /*! Smoothed fraction of the particle mass in a given element */
   double smoothed_metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT];
diff --git a/src/chemistry/GEAR_DIFFUSION/chemistry.h b/src/chemistry/GEAR_DIFFUSION/chemistry.h
new file mode 100644
index 0000000000000000000000000000000000000000..80cbcbc4e40909e07a6727dd7b8cc6fe72926dca
--- /dev/null
+++ b/src/chemistry/GEAR_DIFFUSION/chemistry.h
@@ -0,0 +1,595 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_CHEMISTRY_GEAR_DIFFUSION_H
+#define SWIFT_CHEMISTRY_GEAR_DIFFUSION_H
+
+/**
+ * @file src/chemistry/GEAR_DIFFUSION/chemistry.h
+ * Follows Shen et al. 2010
+ */
+
+/* Some standard headers. */
+#include <float.h>
+#include <math.h>
+#include <string.h>
+
+/* Local includes. */
+#include "chemistry_struct.h"
+#include "error.h"
+#include "hydro.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "units.h"
+
+/**
+ * @brief Copies the chemistry properties of the gas particle over to the
+ * star particle.
+ *
+ * @param p the gas particles.
+ * @param xp the additional properties of the gas particles.
+ * @param sp the new created star particle with its properties.
+ */
+INLINE static void chemistry_copy_star_formation_properties(
+    struct part* p, const struct xpart* xp, struct spart* sp) {
+
+  float mass = hydro_get_mass(p);
+
+  /* Store the chemistry struct in the star particle */
+  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+    sp->chemistry_data.metal_mass_fraction[i] =
+        p->chemistry_data.smoothed_metal_mass_fraction[i];
+
+    /* Remove the metals taken by the star. */
+    p->chemistry_data.metal_mass[i] *= mass / (mass + sp->mass);
+  }
+}
+
+/**
+ * @brief Prints the properties of the chemistry model to stdout.
+ *
+ * @brief The #chemistry_global_data containing information about the current
+ * model.
+ */
+static INLINE void chemistry_print_backend(
+    const struct chemistry_global_data* data) {
+
+  message("Chemistry function is 'Gear with diffusion'.");
+}
+
+/**
+ * @brief Read the solar abundances and scale with them the initial
+ * metallicities.
+ *
+ * @param parameter_file The parsed parameter file.
+ * @param data The properties to initialise.
+ */
+static INLINE void chemistry_scale_initial_metallicities(
+    struct swift_params* parameter_file, struct chemistry_global_data* data) {
+#ifdef HAVE_HDF5
+
+  /* Get the yields table */
+  char filename[DESCRIPTION_BUFFER_SIZE];
+  parser_get_param_string(parameter_file, "GEARFeedback:yields_table",
+                          filename);
+
+  /* Open file. */
+  hid_t file_id = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
+  if (file_id < 0) error("unable to open file %s.\n", filename);
+
+  /* Open group. */
+  hid_t group_id = H5Gopen(file_id, "Data", H5P_DEFAULT);
+  if (group_id < 0) error("unable to open group Data.\n");
+
+  /* Read the data */
+  float* sol_ab = (float*)malloc(sizeof(float) * GEAR_CHEMISTRY_ELEMENT_COUNT);
+  io_read_array_attribute(group_id, "SolarMassAbundances", FLOAT, sol_ab,
+                          GEAR_CHEMISTRY_ELEMENT_COUNT);
+
+  /* Close group */
+  hid_t status = H5Gclose(group_id);
+  if (status < 0) error("error closing group.");
+
+  /* Close file */
+  status = H5Fclose(file_id);
+  if (status < 0) error("error closing file.");
+
+  /* Scale the initial metallicities */
+  char txt[DESCRIPTION_BUFFER_SIZE] = "Scaling initial metallicities by:";
+  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+    data->initial_metallicities[i] *= sol_ab[i];
+    char tmp[10];
+    sprintf(tmp, " %.2g", sol_ab[i]);
+    strcat(txt, tmp);
+  }
+
+  if (engine_rank == 0) {
+    message("%s", txt);
+  }
+#else
+  error("Cannot scale the solar abundances without HDF5");
+#endif
+}
+
+/**
+ * @brief Initialises the chemistry properties.
+ *
+ * Nothing to do here.
+ *
+ * @param parameter_file The parsed parameter file.
+ * @param us The current internal system of units.
+ * @param phys_const The physical constants in internal units.
+ * @param data The properties to initialise.
+ */
+static INLINE void chemistry_init_backend(struct swift_params* parameter_file,
+                                          const struct unit_system* us,
+                                          const struct phys_const* phys_const,
+                                          struct chemistry_global_data* data) {
+
+  /* read parameters */
+  const float initial_metallicity = parser_get_param_float(
+      parameter_file, "GEARChemistry:initial_metallicity");
+
+  if (initial_metallicity < 0) {
+    message("Setting the initial metallicity from the snapshot.");
+  } else {
+    message("Setting the initial metallicity from the parameter file.");
+  }
+
+  /* Set the initial metallicities */
+  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+    data->initial_metallicities[i] = initial_metallicity;
+  }
+
+  /* Read the diffusion coefficient */
+  data->C = parser_get_param_float(parameter_file,
+                                   "GEARChemistry:diffusion_coefficient");
+
+  /* Check if need to scale the initial metallicity */
+  const int scale_metallicity = parser_get_opt_param_int(
+      parameter_file, "GEARChemistry:scale_initial_metallicity", 0);
+
+  /* Scale the metallicities if required */
+  if (scale_metallicity) {
+    chemistry_scale_initial_metallicities(parameter_file, data);
+  }
+}
+
+/**
+ * @brief Prepares a particle for the smooth metal calculation.
+ *
+ * Zeroes all the relevant arrays in preparation for the sums taking place in
+ * the various smooth metallicity tasks
+ *
+ * @param p The particle to act upon
+ * @param cd #chemistry_global_data containing chemistry informations.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_init_part(
+    struct part* restrict p, const struct chemistry_global_data* cd) {
+
+  struct chemistry_part_data* cpd = &p->chemistry_data;
+
+  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+    /* Reset the smoothed metallicity */
+    cpd->smoothed_metal_mass_fraction[i] = 0.f;
+
+    /* Reset the diffusion equation */
+    cpd->metal_mass_dt[i] = 0;
+  }
+
+  /* Reset the shear tensor */
+  for (int i = 0; i < 3; i++) {
+    cpd->S[i][0] = 0;
+    cpd->S[i][1] = 0;
+    cpd->S[i][2] = 0;
+  }
+
+  /* Reset the diffusion. */
+  cpd->diff_coef = 0;
+}
+
+/**
+ * @brief Finishes the smooth metal calculation.
+ *
+ * Multiplies the smoothed metallicity and number of neighbours by the
+ * appropiate constants and add the self-contribution term.
+ *
+ * This function requires the #hydro_end_density to have been called.
+ *
+ * @param p The particle to act upon.
+ * @param cd #chemistry_global_data containing chemistry informations.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_end_density(
+    struct part* restrict p, const struct chemistry_global_data* cd,
+    const struct cosmology* cosmo) {
+
+  /* Some smoothing length multiples. */
+  const float h = p->h;
+  const float h_inv = 1.0f / h;                       /* 1/h */
+  const float h_inv_dim = pow_dimension(h_inv);       /* 1/h^d */
+  const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
+  const float rho = hydro_get_comoving_density(p);
+  const float factor = h_inv_dim / rho; /* 1 / h^d * rho */
+  const float rho_inv = 1.0f / rho;     /* 1 / rho */
+
+  struct chemistry_part_data* cpd = &p->chemistry_data;
+
+  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+    /* Final operation on the density (add self-contribution). */
+    cpd->smoothed_metal_mass_fraction[i] += cpd->metal_mass[i] * kernel_root;
+
+    /* Finish the calculation by inserting the missing h-factors */
+    cpd->smoothed_metal_mass_fraction[i] *= factor;
+  }
+
+  /* convert the shear factor into physical */
+  const float factor_shear = h_inv_dim_plus_one * rho_inv * cosmo->a2_inv;
+  for (int k = 0; k < 3; k++) {
+    cpd->S[k][0] *= factor_shear;
+    cpd->S[k][1] *= factor_shear;
+    cpd->S[k][2] *= factor_shear;
+  }
+
+  /* Compute the trace over 3 and add the hubble flow. */
+  float trace_3 = 0;
+  for (int i = 0; i < 3; i++) {
+    cpd->S[i][i] += cosmo->H;
+    trace_3 += cpd->S[i][i];
+  }
+  trace_3 /= 3.;
+
+  float S_t[3][3];
+  for (int i = 0; i < 3; i++) {
+    /* Make the tensor symmetric. */
+    float avg = 0.5 * (cpd->S[i][0] + cpd->S[0][i]);
+    S_t[i][0] = avg;
+    S_t[0][i] = avg;
+
+    avg = 0.5 * (cpd->S[i][1] + cpd->S[1][i]);
+    S_t[i][1] = avg;
+    S_t[1][i] = avg;
+
+    avg = 0.5 * (cpd->S[i][2] + cpd->S[2][i]);
+    S_t[i][2] = avg;
+    S_t[2][i] = avg;
+
+    /* Remove the trace. */
+    S_t[i][i] -= trace_3;
+  }
+
+  /* Compute the norm. */
+  float norm = 0;
+  for (int i = 0; i < 3; i++) {
+    norm += S_t[i][0] * S_t[i][0];
+    norm += S_t[i][1] * S_t[i][1];
+    norm += S_t[i][2] * S_t[i][2];
+  }
+  norm = sqrt(norm);
+
+  /* Compute the diffusion coefficient in physical coordinates.
+   * The norm is already in physical coordinates.
+   * We do not include kernel_gamma on purpose. */
+  const float h_phys = cosmo->a * p->h;
+  cpd->diff_coef = cd->C * norm * h_phys * h_phys;
+}
+
+/**
+ * @brief Updates to the chemistry data after the hydro force loop.
+ *
+ * @param p The particle to act upon.
+ * @param cosmo The current cosmological model.
+ * @param with_cosmology Are we running with the cosmology?
+ * @param time Current time of the simulation.
+ * @param dt Time step (in physical units).
+ */
+__attribute__((always_inline)) INLINE static void chemistry_end_force(
+    struct part* restrict p, const struct cosmology* cosmo,
+    const int with_cosmology, const double time, const double dt) {
+  if (dt == 0) {
+    return;
+  }
+
+  struct chemistry_part_data* ch = &p->chemistry_data;
+  const float h_inv = cosmo->a / p->h;
+  const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+  /* Missing factors in iact. */
+  const float factor = h_inv_dim * h_inv;
+
+  const double sum = 0;
+  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+    ch->metal_mass[i] += ch->metal_mass_dt[i] * dt * factor;
+    /* Make sure that the metallicity is 0 <= x <= 1 */
+    if (ch->metal_mass[i] < 0 || ch->metal_mass[i] > hydro_get_mass(p)) {
+      error("Negative mass or mass fraction larger than 1.");
+    }
+    /* Make sure that we do not have more metals than the sum. */
+    if (i != GEAR_CHEMISTRY_ELEMENT_COUNT) {
+      sum += ch->metal_mass[i];
+    } else if (sum > ch->metal_mass[i]) {
+      error("Found more individual elements than the sum of all of them.");
+    }
+  }
+}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param cd #chemistry_global_data containing chemistry informations.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void
+chemistry_part_has_no_neighbours(struct part* restrict p,
+                                 struct xpart* restrict xp,
+                                 const struct chemistry_global_data* cd,
+                                 const struct cosmology* cosmo) {
+
+  /* Set the smoothed fractions with the non smoothed fractions */
+  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+    p->chemistry_data.smoothed_metal_mass_fraction[i] =
+        p->chemistry_data.metal_mass[i] / hydro_get_mass(p);
+  }
+}
+
+/**
+ * @brief Computes the chemistry-related time-step constraint.
+ *
+ * No constraints in the GEAR model (no diffusion) --> FLT_MAX
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param cosmo The current cosmological model.
+ * @param us The internal system of units.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param cd The global properties of the chemistry scheme.
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float chemistry_timestep(
+    const struct phys_const* restrict phys_const,
+    const struct cosmology* restrict cosmo,
+    const struct unit_system* restrict us,
+    const struct hydro_props* hydro_props,
+    const struct chemistry_global_data* cd, const struct part* restrict p) {
+  return FLT_MAX;
+}
+
+/**
+ * @brief Sets the chemistry properties of the (x-)particles to a valid start
+ * state.
+ *
+ * @param phys_const The #phys_const.
+ * @param us The #unit_system.
+ * @param cosmo The #cosmology.
+ * @param data The global chemistry information.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_first_init_part(
+    const struct phys_const* restrict phys_const,
+    const struct unit_system* restrict us,
+    const struct cosmology* restrict cosmo,
+    const struct chemistry_global_data* data, struct part* restrict p,
+    struct xpart* restrict xp) {
+
+  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+    if (data->initial_metallicities[i] < 0) {
+      /* Use the value from the IC. We are reading the metal mass fraction. */
+      p->chemistry_data.metal_mass[i] *= hydro_get_mass(p);
+    } else {
+      /* Use the value from the parameter file */
+      p->chemistry_data.metal_mass[i] =
+          data->initial_metallicities[i] * hydro_get_mass(p);
+    }
+  }
+
+  chemistry_init_part(p, data);
+}
+
+/**
+ * @brief Sets the chemistry properties of the sparticles to a valid start
+ * state.
+ *
+ * @param data The global chemistry information.
+ * @param sp Pointer to the sparticle data.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_first_init_spart(
+    const struct chemistry_global_data* data, struct spart* restrict sp) {
+
+  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+    sp->chemistry_data.metal_mass_fraction[i] =
+        data->initial_metallicities[i] * sp->mass;
+  }
+}
+
+/**
+ * @brief Initialise the chemistry properties of a black hole with
+ * the chemistry properties of the gas it is born from.
+ *
+ * Nothing to do here.
+ *
+ * @param bp_data The black hole data to initialise.
+ * @param p_data The gas data to use.
+ * @param gas_mass The mass of the gas particle.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_bpart_from_part(
+    struct chemistry_bpart_data* bp_data,
+    const struct chemistry_part_data* p_data, const double gas_mass) {
+  error("Loic: to be implemented");
+}
+
+/**
+ * @brief Add the chemistry data of a gas particle to a black hole.
+ *
+ * Nothing to do here.
+ *
+ * @param bp_data The black hole data to add to.
+ * @param p_data The gas data to use.
+ * @param gas_mass The mass of the gas particle.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_add_part_to_bpart(
+    struct chemistry_bpart_data* bp_data,
+    const struct chemistry_part_data* p_data, const double gas_mass) {
+  error("Loic: to be implemented");
+}
+
+/**
+ * @brief Add the chemistry data of a black hole to another one.
+ *
+ * Nothing to do here.
+ *
+ * @param bp_data The black hole data to add to.
+ * @param swallowed_data The black hole data to use.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_add_bpart_to_bpart(
+    struct chemistry_bpart_data* bp_data,
+    const struct chemistry_bpart_data* swallowed_data) {
+  error("Loic: to be implemented");
+}
+
+/**
+ * @brief Split the metal content of a particle into n pieces
+ *
+ * @param p The #part.
+ * @param n The number of pieces to split into.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_split_part(
+    struct part* p, const double n) {
+  error("Loic: to be implemented");
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * star particle to be used in feedback/enrichment related routines.
+ *
+ * @param sp Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static double
+chemistry_get_total_metal_mass_fraction_for_feedback(
+    const struct spart* restrict sp) {
+
+  return sp->chemistry_data
+      .metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT - 1];
+}
+
+/**
+ * @brief Returns the abundances (metal mass fraction) of the
+ * star particle to be used in feedback/enrichment related routines.
+ *
+ * @param sp Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static double const*
+chemistry_get_metal_mass_fraction_for_feedback(
+    const struct spart* restrict sp) {
+
+  return sp->chemistry_data.metal_mass_fraction;
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * gas particle to be used in cooling related routines.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static double
+chemistry_get_total_metal_mass_fraction_for_cooling(
+    const struct part* restrict p) {
+
+  return p->chemistry_data
+      .smoothed_metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT - 1];
+}
+
+/**
+ * @brief Returns the abundance array (metal mass fractions) of the
+ * gas particle to be used in cooling related routines.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static double const*
+chemistry_get_metal_mass_fraction_for_cooling(const struct part* restrict p) {
+
+  return p->chemistry_data.smoothed_metal_mass_fraction;
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * gas particle to be used in star formation related routines.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static double
+chemistry_get_total_metal_mass_fraction_for_star_formation(
+    const struct part* restrict p) {
+
+  return p->chemistry_data
+      .smoothed_metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT - 1];
+}
+
+/**
+ * @brief Returns the abundance array (metal mass fractions) of the
+ * gas particle to be used in star formation related routines.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static double const*
+chemistry_get_metal_mass_fraction_for_star_formation(
+    const struct part* restrict p) {
+
+  return p->chemistry_data.smoothed_metal_mass_fraction;
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * black hole particle to be used in the stats related routines.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_bh_total_metal_mass_for_stats(const struct bpart* restrict bp) {
+  error("Not implemented");
+  return 0.f;
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * gas particle to be used in the stats related routines.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_total_metal_mass_for_stats(const struct part* restrict p) {
+
+  return p->chemistry_data.metal_mass[GEAR_CHEMISTRY_ELEMENT_COUNT - 1];
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * star particle to be used in the stats related routines.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_star_total_metal_mass_for_stats(const struct spart* restrict sp) {
+
+  return sp->chemistry_data
+             .metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT - 1] *
+         sp->mass;
+}
+
+#endif /* SWIFT_CHEMISTRY_GEAR_DIFFUSION_H */
diff --git a/src/chemistry/GEAR_DIFFUSION/chemistry_iact.h b/src/chemistry/GEAR_DIFFUSION/chemistry_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..896307c011e2e501e03a076e8662bd9e7bdd101c
--- /dev/null
+++ b/src/chemistry/GEAR_DIFFUSION/chemistry_iact.h
@@ -0,0 +1,297 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GEAR_DIFFUSION_CHEMISTRY_IACT_H
+#define SWIFT_GEAR_DIFFUSION_CHEMISTRY_IACT_H
+
+/**
+ * @file GEAR/chemistry_iact.h
+ * @brief Smooth metal interaction + diffusion functions following the GEAR
+ * version of smooth metalicity.
+ */
+
+/**
+ * @brief do chemistry computation after the runner_iact_density (symmetric
+ * version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_chemistry(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H) {
+
+  struct chemistry_part_data *chi = &pi->chemistry_data;
+  struct chemistry_part_data *chj = &pj->chemistry_data;
+
+  float wi, wi_dx;
+  float wj, wj_dx;
+
+  /* Get the masses. */
+  const float mj = hydro_get_mass(pj);
+  const float mi = hydro_get_mass(pi);
+
+  /* Get r */
+  const float r = sqrtf(r2);
+  const float r_inv = 1.f / r;
+
+  /* Compute the kernel function for pi */
+  const float ui = r / hi;
+  kernel_deval(ui, &wi, &wi_dx);
+
+  /* Compute the kernel function for pj */
+  const float uj = r / hj;
+  kernel_deval(uj, &wj, &wj_dx);
+
+  const float wi_dr = wi_dx * r_inv;
+  const float mj_wi_dr = mj * wi_dr;
+
+  const float wj_dr = wj_dx * r_inv;
+  const float mi_wj_dr = mi * wj_dr;
+
+  /* Compute contribution to the smooth metallicity */
+  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+    chi->smoothed_metal_mass_fraction[i] += chj->metal_mass[i] * wi;
+    chj->smoothed_metal_mass_fraction[i] += chi->metal_mass[i] * wj;
+  }
+
+  /* Compute the shear tensor */
+  for (int i = 0; i < 3; i++) {
+    chi->S[i][0] += (pj->v[0] - pi->v[0]) * dx[i] * mj_wi_dr;
+    chi->S[i][1] += (pj->v[1] - pi->v[1]) * dx[i] * mj_wi_dr;
+    chi->S[i][2] += (pj->v[2] - pi->v[2]) * dx[i] * mj_wi_dr;
+
+    chj->S[i][0] -= (pj->v[0] - pi->v[0]) * dx[i] * mi_wj_dr;
+    chj->S[i][1] -= (pj->v[1] - pi->v[1]) * dx[i] * mi_wj_dr;
+    chj->S[i][2] -= (pj->v[2] - pi->v[2]) * dx[i] * mi_wj_dr;
+  }
+}
+
+/**
+ * @brief do chemistry computation after the runner_iact_density (non symmetric
+ * version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle (not updated).
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    const struct part *restrict pj, float a, float H) {
+
+  struct chemistry_part_data *chi = &pi->chemistry_data;
+  const struct chemistry_part_data *chj = &pj->chemistry_data;
+
+  float wi, wi_dx;
+
+  /* Get the masses. */
+  const float mj = hydro_get_mass(pj);
+
+  /* Get r */
+  const float r = sqrtf(r2);
+  const float r_inv = 1.f / r;
+
+  /* Compute the kernel function for pi */
+  const float ui = r / hi;
+  kernel_deval(ui, &wi, &wi_dx);
+
+  const float wi_dr = wi_dx * r_inv;
+  const float mj_wi_dr = mj * wi_dr;
+
+  /* Compute contribution to the smooth metallicity */
+  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+    chi->smoothed_metal_mass_fraction[i] += chj->metal_mass[i] * wi;
+  }
+
+  /* Compute the shear tensor */
+  for (int i = 0; i < 3; i++) {
+    chi->S[i][0] += (pj->v[0] - pi->v[0]) * dx[i] * mj_wi_dr;
+    chi->S[i][1] += (pj->v[1] - pi->v[1]) * dx[i] * mj_wi_dr;
+    chi->S[i][2] += (pj->v[2] - pi->v[2]) * dx[i] * mj_wi_dr;
+  }
+}
+
+/**
+ * @brief do metal diffusion computation in the <FORCE LOOP>
+ * (symmetric version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ * @param time_base The time base used in order to convert integer to float
+ * time.
+ * @param ti_current The current time (in integer)
+ * @param cosmo The #cosmology.
+ * @param with_cosmology Are we running with cosmology?
+ *
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_diffusion(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H, float time_base,
+    integertime_t t_current, const struct cosmology *cosmo,
+    const int with_cosmology) {
+
+  struct chemistry_part_data *chi = &pi->chemistry_data;
+  struct chemistry_part_data *chj = &pj->chemistry_data;
+
+  /* No need to diffuse if both particles are not diffusing. */
+  if (chj->diff_coef > 0 && chi->diff_coef > 0) {
+
+    /* Get mass */
+    const float mj = hydro_get_mass(pj);
+    const float mi = hydro_get_mass(pi);
+    const float rhoj = hydro_get_physical_density(pj, cosmo);
+    const float rhoi = hydro_get_physical_density(pi, cosmo);
+
+    float wi, wj, dwi_dx, dwj_dx;
+
+    /* Get r */
+    const float r = sqrtf(r2);
+
+    /* part j */
+    /* Get the kernel for hj */
+    const float hj_inv = 1.0f / hj;
+
+    /* Compute the kernel function for pj */
+    const float xj = r * hj_inv;
+    kernel_deval(xj, &wj, &dwj_dx);
+
+    /* part i */
+    /* Get the kernel for hi */
+    const float hi_inv = 1.0f / hi;
+
+    /* Compute the kernel function for pi */
+    const float xi = r * hi_inv;
+    kernel_deval(xi, &wi, &dwi_dx);
+
+    /* Get 1/r */
+    const float r_inv = 1.f / sqrtf(r2);
+
+    const float wi_dr = dwi_dx * r_inv;
+    const float wj_dr = dwj_dx * r_inv;
+
+    const float mj_dw_r = mj * wi_dr;
+    const float mi_dw_r = mi * wj_dr;
+
+    /* Compute the diffusion coefficient <D> / <rho> in physical units. */
+    const float coef = 2. * (chi->diff_coef + chj->diff_coef) / (rhoi + rhoj);
+
+    const float coef_i = coef * mj_dw_r;
+    const float coef_j = coef * mi_dw_r;
+
+    /* Compute the time derivative */
+    const float inv_mi = 1.f / hydro_get_mass(pi);
+    const float inv_mj = 1.f / hydro_get_mass(pj);
+    for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+      const double mi_frac = chi->metal_mass[i] * inv_mi;
+      const double mj_frac = chj->metal_mass[i] * inv_mj;
+      const double dm = mi_frac - mj_frac;
+      chi->metal_mass_dt[i] += coef_i * dm;
+      chj->metal_mass_dt[i] -= coef_j * dm;
+    }
+  }
+}
+
+/**
+ * @brief do metal diffusion computation in the <FORCE LOOP>
+ * (nonsymmetric version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ * @param time_base The time base used in order to convert integer to float
+ * time.
+ * @param ti_current The current time (in integer)
+ * @param cosmo The #cosmology.
+ * @param with_cosmology Are we running with cosmology?
+ *
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_diffusion(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H, float time_base,
+    integertime_t t_current, const struct cosmology *cosmo,
+    const int with_cosmology) {
+
+  struct chemistry_part_data *chi = &pi->chemistry_data;
+  struct chemistry_part_data *chj = &pj->chemistry_data;
+
+  if (chj->diff_coef > 0 && chi->diff_coef > 0) {
+
+    /* Get mass */
+    const float mj = hydro_get_mass(pj);
+    const float rhoj = hydro_get_physical_density(pj, cosmo);
+    const float rhoi = hydro_get_physical_density(pi, cosmo);
+
+    float wi, dwi_dx;
+
+    /* Get r */
+    const float r = sqrtf(r2);
+
+    /* part i */
+    /* Get the kernel for hi */
+    const float hi_inv = 1.0f / hi;
+
+    /* Compute the kernel function for pi */
+    const float xi = r * hi_inv;
+    kernel_deval(xi, &wi, &dwi_dx);
+
+    /* Get 1/r */
+    const float r_inv = 1.f / sqrtf(r2);
+    const float wi_dr = dwi_dx * r_inv;
+
+    const float mj_dw_r = mj * wi_dr;
+
+    /* Compute the diffusion coefficient <D> / <rho> in physical units. */
+    const float coef = 2. * (chi->diff_coef + chj->diff_coef) / (rhoi + rhoj);
+
+    const float coef_i = coef * mj_dw_r;
+
+    /* Compute the time derivative */
+    const float inv_mi = 1.f / hydro_get_mass(pi);
+    const float inv_mj = 1.f / hydro_get_mass(pj);
+    for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+      const double mi_frac = chi->metal_mass[i] * inv_mi;
+      const double mj_frac = chj->metal_mass[i] * inv_mj;
+      const double dm = mi_frac - mj_frac;
+      chi->metal_mass_dt[i] += coef_i * dm;
+    }
+  }
+}
+
+#endif /* SWIFT_GEAR_DIFFUSION_CHEMISTRY_IACT_H */
diff --git a/src/chemistry/GEAR_DIFFUSION/chemistry_io.h b/src/chemistry/GEAR_DIFFUSION/chemistry_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..29761c0efc759533b7cf53f105f7d6575ee9d790
--- /dev/null
+++ b/src/chemistry/GEAR_DIFFUSION/chemistry_io.h
@@ -0,0 +1,160 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_CHEMISTRY_IO_GEAR_DIFFUSION_H
+#define SWIFT_CHEMISTRY_IO_GEAR_DIFFUSION_H
+
+#include "chemistry_struct.h"
+#include "engine.h"
+#include "error.h"
+#include "feedback.h"
+#include "io_properties.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "units.h"
+
+/**
+ * @brief Specifies which particle fields to read from a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to read.
+ *
+ * @return Returns the number of fields to read.
+ */
+INLINE static int chemistry_read_particles(struct part* parts,
+                                           struct io_props* list) {
+
+  /* List what we want to read */
+  list[0] = io_make_input_field(
+      "MetalMassFraction", DOUBLE, GEAR_CHEMISTRY_ELEMENT_COUNT, OPTIONAL,
+      UNIT_CONV_NO_UNITS, parts, chemistry_data.metal_mass);
+
+  return 1;
+}
+
+INLINE static void convert_gas_metals(const struct engine* e,
+                                      const struct part* p,
+                                      const struct xpart* xp, double* ret) {
+
+  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+    ret[i] = p->chemistry_data.metal_mass[i] / hydro_get_mass(p);
+  }
+}
+
+/**
+ * @brief Specifies which particle fields to write to a dataset
+ *
+ * @param parts The particle array.
+ * @param xparts The extra particle array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return Returns the number of fields to write.
+ */
+INLINE static int chemistry_write_particles(const struct part* parts,
+                                            const struct xpart* xparts,
+                                            struct io_props* list,
+                                            const int with_cosmology) {
+
+  /* List what we want to write */
+  list[0] = io_make_output_field(
+      "SmoothedMetalMassFractions", DOUBLE, GEAR_CHEMISTRY_ELEMENT_COUNT,
+      UNIT_CONV_NO_UNITS, 0.f, parts,
+      chemistry_data.smoothed_metal_mass_fraction,
+      "Mass fraction of each element smoothed over the neighbors");
+
+  list[1] = io_make_output_field_convert_part(
+      "MetalMassFractions", DOUBLE, GEAR_CHEMISTRY_ELEMENT_COUNT,
+      UNIT_CONV_NO_UNITS, 0.f, parts, xparts, convert_gas_metals,
+      "Mass fraction of each element");
+
+  return 2;
+}
+
+/**
+ * @brief Specifies which sparticle fields to write to a dataset
+ *
+ * @param sparts The sparticle array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return Returns the number of fields to write.
+ */
+INLINE static int chemistry_write_sparticles(const struct spart* sparts,
+                                             struct io_props* list) {
+
+  /* List what we want to write */
+  list[0] = io_make_output_field(
+      "MetalMassFractions", DOUBLE, GEAR_CHEMISTRY_ELEMENT_COUNT,
+      UNIT_CONV_NO_UNITS, 0.f, sparts, chemistry_data.metal_mass_fraction,
+      "Mass fraction of each element");
+
+  return 1;
+}
+
+/**
+ * @brief Specifies which black hole particle fields to write to a dataset
+ *
+ * @param bparts The black hole particle array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return Returns the number of fields to write.
+ */
+INLINE static int chemistry_write_bparticles(const struct bpart* bparts,
+                                             struct io_props* list) {
+
+  /* No fields to write here */
+  return 0;
+}
+
+#ifdef HAVE_HDF5
+
+/**
+ * @brief Writes the current model of SPH to the file
+ * @param h_grp The HDF5 group in which to write
+ * @param h_grp_columns The HDF5 group containing named columns
+ * @param e The #engine.
+ */
+INLINE static void chemistry_write_flavour(hid_t h_grp, hid_t h_grp_columns,
+                                           const struct engine* e) {
+
+  io_write_attribute_s(h_grp, "Chemistry Model", "GEAR with diffusion");
+  io_write_attribute_d(h_grp, "Chemistry element count",
+                       GEAR_CHEMISTRY_ELEMENT_COUNT);
+#ifdef FEEDBACK_GEAR
+  const char* element_names = e->feedback_props->stellar_model.elements_name;
+
+  /* Add to the named columns */
+  hsize_t dims[1] = {GEAR_CHEMISTRY_ELEMENT_COUNT};
+  hid_t type = H5Tcopy(H5T_C_S1);
+  H5Tset_size(type, GEAR_LABELS_SIZE);
+  hid_t space = H5Screate_simple(1, dims, NULL);
+  hid_t dset = H5Dcreate(h_grp_columns, "ElementAbundances", type, space,
+                         H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+  H5Dwrite(dset, type, H5S_ALL, H5S_ALL, H5P_DEFAULT, element_names);
+  H5Dclose(dset);
+  dset = H5Dcreate(h_grp_columns, "SmoothedElementAbundances", type, space,
+                   H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+  H5Dwrite(dset, type, H5S_ALL, H5S_ALL, H5P_DEFAULT, element_names);
+  H5Dclose(dset);
+
+  H5Tclose(type);
+#endif
+}
+#endif
+
+#endif /* SWIFT_CHEMISTRY_IO_GEAR_DIFFUSION_H */
diff --git a/src/chemistry/GEAR_DIFFUSION/chemistry_struct.h b/src/chemistry/GEAR_DIFFUSION/chemistry_struct.h
new file mode 100644
index 0000000000000000000000000000000000000000..eb8b1e255e4d4130ef843251dff34bd199e42295
--- /dev/null
+++ b/src/chemistry/GEAR_DIFFUSION/chemistry_struct.h
@@ -0,0 +1,70 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_CHEMISTRY_STRUCT_GEAR_DIFFUSION_H
+#define SWIFT_CHEMISTRY_STRUCT_GEAR_DIFFUSION_H
+
+/**
+ * @brief Global chemical abundance information.
+ */
+struct chemistry_global_data {
+
+  /* Initial mass fraction */
+  double initial_metallicities[GEAR_CHEMISTRY_ELEMENT_COUNT];
+
+  /*! Diffusion coefficent (no unit) */
+  float C;
+};
+
+/**
+ * @brief Properties of the chemistry function for #part.
+ */
+struct chemistry_part_data {
+
+  /*! Total mass of element in a particle.
+    This field is available only outside the density hydro loop. */
+  double metal_mass[GEAR_CHEMISTRY_ELEMENT_COUNT];
+
+  /*! Smoothed fraction of the particle mass in a given element */
+  double smoothed_metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT];
+
+  /*! Diffusion coefficient */
+  float diff_coef;
+
+  /*! Variation of the metal mass */
+  double metal_mass_dt[GEAR_CHEMISTRY_ELEMENT_COUNT];
+
+  /*! shear tensor in internal and physical units. */
+  float S[3][3];
+};
+
+/**
+ * @brief Properties of the chemistry function for #spart.
+ */
+struct chemistry_spart_data {
+
+  /*! Fraction of the particle mass in a given element */
+  double metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT];
+};
+
+/**
+ * @brief Chemical abundances traced by the #bpart in the GEAR model.
+ */
+struct chemistry_bpart_data {};
+
+#endif /* SWIFT_CHEMISTRY_STRUCT_GEAR_DIFFUSION_H */
diff --git a/src/chemistry/QLA/chemistry.h b/src/chemistry/QLA/chemistry.h
index c9823214570a82649c4db77049e9d611c26e94d7..7d6ffba8facf8d1e3b1e6d2abafbd55a2f3c2ba8 100644
--- a/src/chemistry/QLA/chemistry.h
+++ b/src/chemistry/QLA/chemistry.h
@@ -94,11 +94,17 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density(
 /**
  * @brief Updates to the chemistry data after the hydro force loop.
  *
+ * Nothing to do here.
+ *
  * @param p The particle to act upon.
  * @param cosmo The current cosmological model.
+ * @param with_cosmology Are we running with the cosmology?
+ * @param time Current time of the simulation.
+ * @param dt Time step (in physical units).
  */
 __attribute__((always_inline)) INLINE static void chemistry_end_force(
-    struct part* restrict p, const struct cosmology* cosmo) {}
+    struct part* restrict p, const struct cosmology* cosmo,
+    const int with_cosmology, const double time, const double dt) {}
 
 /**
  * @brief Computes the chemistry-related time-step constraint.
diff --git a/src/chemistry/QLA/chemistry_iact.h b/src/chemistry/QLA/chemistry_iact.h
index 89f0d4a3c8a2dd4fbde40b4525196d439ba66f85..1c3dcfed72fd31312f5b7bc337e599991c017afe 100644
--- a/src/chemistry/QLA/chemistry_iact.h
+++ b/src/chemistry/QLA/chemistry_iact.h
@@ -58,4 +58,54 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
     float r2, const float *dx, float hi, float hj, struct part *restrict pi,
     const struct part *restrict pj, float a, float H) {}
 
+/**
+ * @brief do metal diffusion computation in the <FORCE LOOP>
+ * (symmetric version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ * @param time_base The time base used in order to convert integer to float
+ * time.
+ * @param ti_current The current time (in integer)
+ * @param cosmo The #cosmology.
+ * @param with_cosmology Are we running with cosmology?
+ *
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_diffusion(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H, float time_base,
+    integertime_t t_current, const struct cosmology *cosmo,
+    const int with_cosmology) {}
+
+/**
+ * @brief do metal diffusion computation in the <FORCE LOOP>
+ * (nonsymmetric version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ * @param time_base The time base used in order to convert integer to float
+ * time.
+ * @param ti_current The current time (in integer)
+ * @param cosmo The #cosmology.
+ * @param with_cosmology Are we running with cosmology?
+ *
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_diffusion(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H, float time_base,
+    integertime_t t_current, const struct cosmology *cosmo,
+    const int with_cosmology) {}
+
 #endif /* SWIFT_QLA_CHEMISTRY_IACT_H */
diff --git a/src/chemistry/none/chemistry.h b/src/chemistry/none/chemistry.h
index 3cd78a2cc63af5caf62f439ad8df5cf05cc33683..fa6e0c892d1d36d472306bb228ba803cbdbe112c 100644
--- a/src/chemistry/none/chemistry.h
+++ b/src/chemistry/none/chemistry.h
@@ -101,11 +101,17 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density(
 /**
  * @brief Updates to the chemistry data after the hydro force loop.
  *
+ * Nothing to do here.
+ *
  * @param p The particle to act upon.
  * @param cosmo The current cosmological model.
+ * @param with_cosmology Are we running with the cosmology?
+ * @param time Current time of the simulation.
+ * @param dt Time step (in physical units).
  */
 __attribute__((always_inline)) INLINE static void chemistry_end_force(
-    struct part* restrict p, const struct cosmology* cosmo) {}
+    struct part* restrict p, const struct cosmology* cosmo,
+    const int with_cosmology, const double time, const double dt) {}
 
 /**
  * @brief Computes the chemistry-related time-step constraint.
diff --git a/src/chemistry/none/chemistry_iact.h b/src/chemistry/none/chemistry_iact.h
index 52499c5e5cc3d5749904251cc967193d875579b9..13233f6a0f0d4657f49a51b2da5bc456bd829759 100644
--- a/src/chemistry/none/chemistry_iact.h
+++ b/src/chemistry/none/chemistry_iact.h
@@ -59,4 +59,54 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
     float r2, const float *dx, float hi, float hj, struct part *restrict pi,
     const struct part *restrict pj, float a, float H) {}
 
+/**
+ * @brief do metal diffusion computation in the <FORCE LOOP>
+ * (symmetric version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ * @param time_base The time base used in order to convert integer to float
+ * time.
+ * @param ti_current The current time (in integer)
+ * @param cosmo The #cosmology.
+ * @param with_cosmology Are we running with cosmology?
+ *
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_diffusion(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H, float time_base,
+    integertime_t t_current, const struct cosmology *cosmo,
+    const int with_cosmology) {}
+
+/**
+ * @brief do metal diffusion computation in the <FORCE LOOP>
+ * (nonsymmetric version)
+ *
+ * @param r2 Comoving square distance between the two particles.
+ * @param dx Comoving vector separating both particles (pi - pj).
+ * @param hi Comoving smoothing-length of particle i.
+ * @param hj Comoving smoothing-length of particle j.
+ * @param pi First particle.
+ * @param pj Second particle.
+ * @param a Current scale factor.
+ * @param H Current Hubble parameter.
+ * @param time_base The time base used in order to convert integer to float
+ * time.
+ * @param ti_current The current time (in integer)
+ * @param cosmo The #cosmology.
+ * @param with_cosmology Are we running with cosmology?
+ *
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_diffusion(
+    float r2, const float *dx, float hi, float hj, struct part *restrict pi,
+    struct part *restrict pj, float a, float H, float time_base,
+    integertime_t t_current, const struct cosmology *cosmo,
+    const int with_cosmology) {}
+
 #endif /* SWIFT_NONE_CHEMISTRY_IACT_H */
diff --git a/src/chemistry_io.h b/src/chemistry_io.h
index c94461c84bdf243990fa8f28f0049178e180d9ec..d7386f8c68a42b797a8d201e8a68329f00b4f6fc 100644
--- a/src/chemistry_io.h
+++ b/src/chemistry_io.h
@@ -27,6 +27,8 @@
 #include "./chemistry/none/chemistry_io.h"
 #elif defined(CHEMISTRY_GEAR)
 #include "./chemistry/GEAR/chemistry_io.h"
+#elif defined(CHEMISTRY_GEAR_DIFFUSION)
+#include "./chemistry/GEAR_DIFFUSION/chemistry_io.h"
 #elif defined(CHEMISTRY_QLA)
 #include "./chemistry/QLA/chemistry_io.h"
 #elif defined(CHEMISTRY_EAGLE)
diff --git a/src/chemistry_struct.h b/src/chemistry_struct.h
index 05074d8759be660876544aaa20cffec6cddb8bce..0b1cb59e592d9dbce8e77fb01cf45e1be846a55b 100644
--- a/src/chemistry_struct.h
+++ b/src/chemistry_struct.h
@@ -32,6 +32,8 @@
 #include "./chemistry/none/chemistry_struct.h"
 #elif defined(CHEMISTRY_GEAR)
 #include "./chemistry/GEAR/chemistry_struct.h"
+#elif defined(CHEMISTRY_GEAR_DIFFUSION)
+#include "./chemistry/GEAR_DIFFUSION/chemistry_struct.h"
 #elif defined(CHEMISTRY_QLA)
 #include "./chemistry/QLA/chemistry_struct.h"
 #elif defined(CHEMISTRY_EAGLE)
diff --git a/src/runner_doiact_functions_hydro.h b/src/runner_doiact_functions_hydro.h
index b0dbbc1b97ba91975fc2669e547ee9ddf6fd10f4..4bdbcec05cf8fd9e5d5b06b13737d02886c94164 100644
--- a/src/runner_doiact_functions_hydro.h
+++ b/src/runner_doiact_functions_hydro.h
@@ -40,6 +40,11 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci,
 
   const struct engine *e = r->e;
   const struct cosmology *cosmo = e->cosmology;
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
+  const double time_base = e->time_base;
+  const integertime_t t_current = e->ti_current;
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
+#endif
 
   TIMER_TIC;
 
@@ -119,6 +124,8 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci,
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
         runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H);
+        runner_iact_nonsym_diffusion(r2, dx, hi, hj, pi, pj, a, H, time_base,
+                                     t_current, cosmo, with_cosmology);
 #endif
       }
       if (r2 < hjg2 && pj_active) {
@@ -135,6 +142,8 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci,
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
         runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H);
+        runner_iact_nonsym_diffusion(r2, dx, hj, hi, pj, pi, a, H, time_base,
+                                     t_current, cosmo, with_cosmology);
 #endif
       }
     } /* loop over the parts in cj. */
@@ -157,6 +166,11 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
 
   const struct engine *e = r->e;
   const struct cosmology *cosmo = e->cosmology;
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
+  const double time_base = e->time_base;
+  const integertime_t t_current = e->ti_current;
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
+#endif
 
   TIMER_TIC;
 
@@ -238,6 +252,8 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
           runner_iact_timebin(r2, dx, hi, hj, pi, pj, a, H);
+          runner_iact_diffusion(r2, dx, hi, hj, pi, pj, a, H, time_base,
+                                t_current, cosmo, with_cosmology);
 #endif
         } else if (pi_active) {
 
@@ -249,6 +265,8 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
           runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H);
+          runner_iact_nonsym_diffusion(r2, dx, hi, hj, pi, pj, a, H, time_base,
+                                       t_current, cosmo, with_cosmology);
 #endif
         } else if (pj_active) {
 
@@ -264,6 +282,8 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
           runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H);
+          runner_iact_nonsym_diffusion(r2, dx, hj, hi, pj, pi, a, H, time_base,
+                                       t_current, cosmo, with_cosmology);
 #endif
         }
       }
@@ -285,6 +305,11 @@ void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) {
 
   const struct engine *e = r->e;
   const struct cosmology *cosmo = e->cosmology;
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
+  const double time_base = e->time_base;
+  const integertime_t t_current = e->ti_current;
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
+#endif
 
   TIMER_TIC;
 
@@ -356,6 +381,8 @@ void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) {
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
         runner_iact_timebin(r2, dx, hi, hj, pi, pj, a, H);
+        runner_iact_diffusion(r2, dx, hi, hj, pi, pj, a, H, time_base,
+                              t_current, cosmo, with_cosmology);
 #endif
       } else if (doi) {
 
@@ -367,6 +394,8 @@ void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) {
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
         runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H);
+        runner_iact_nonsym_diffusion(r2, dx, hi, hj, pi, pj, a, H, time_base,
+                                     t_current, cosmo, with_cosmology);
 #endif
       } else if (doj) {
 
@@ -382,6 +411,8 @@ void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) {
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
         runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H);
+        runner_iact_nonsym_diffusion(r2, dx, hj, hi, pj, pi, a, H, time_base,
+                                     t_current, cosmo, with_cosmology);
 #endif
       }
     } /* loop over the parts in cj. */
@@ -402,6 +433,11 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) {
 
   const struct engine *e = r->e;
   const struct cosmology *cosmo = e->cosmology;
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
+  const double time_base = e->time_base;
+  const integertime_t t_current = e->ti_current;
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
+#endif
 
   TIMER_TIC;
 
@@ -473,6 +509,8 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) {
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
         runner_iact_timebin(r2, dx, hi, hj, pi, pj, a, H);
+        runner_iact_diffusion(r2, dx, hi, hj, pi, pj, a, H, time_base,
+                              t_current, cosmo, with_cosmology);
 #endif
       } else if (doi) {
 
@@ -484,6 +522,8 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) {
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
         runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H);
+        runner_iact_nonsym_diffusion(r2, dx, hi, hj, pi, pj, a, H, time_base,
+                                     t_current, cosmo, with_cosmology);
 #endif
       } else if (doj) {
 
@@ -499,6 +539,8 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) {
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
         runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H);
+        runner_iact_nonsym_diffusion(r2, dx, hj, hi, pj, pi, a, H, time_base,
+                                     t_current, cosmo, with_cosmology);
 #endif
       }
     } /* loop over the parts in cj. */
@@ -528,6 +570,11 @@ void DOPAIR_SUBSET_NAIVE(struct runner *r, struct cell *restrict ci,
 
   const struct engine *e = r->e;
   const struct cosmology *cosmo = e->cosmology;
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
+  const double time_base = e->time_base;
+  const integertime_t t_current = e->ti_current;
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
+#endif
 
   TIMER_TIC;
 
@@ -589,6 +636,8 @@ void DOPAIR_SUBSET_NAIVE(struct runner *r, struct cell *restrict ci,
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
         runner_iact_nonsym_timebin(r2, dx, hi, pj->h, pi, pj, a, H);
+        runner_iact_nonsym_diffusion(r2, dx, hi, pj->h, pi, pj, a, H, time_base,
+                                     t_current, cosmo, with_cosmology);
 #endif
       }
     } /* loop over the parts in cj. */
@@ -618,6 +667,11 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
 
   const struct engine *e = r->e;
   const struct cosmology *cosmo = e->cosmology;
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
+  const double time_base = e->time_base;
+  const integertime_t t_current = e->ti_current;
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
+#endif
 
   TIMER_TIC;
 
@@ -686,6 +740,8 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
           runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H);
+          runner_iact_nonsym_diffusion(r2, dx, hi, hj, pi, pj, a, H, time_base,
+                                       t_current, cosmo, with_cosmology);
 #endif
         }
       } /* loop over the parts in cj. */
@@ -746,6 +802,8 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
           runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H);
+          runner_iact_nonsym_diffusion(r2, dx, hi, hj, pi, pj, a, H, time_base,
+                                       t_current, cosmo, with_cosmology);
 #endif
         }
       } /* loop over the parts in cj. */
@@ -831,6 +889,11 @@ void DOSELF_SUBSET(struct runner *r, struct cell *restrict ci,
 
   const struct engine *e = r->e;
   const struct cosmology *cosmo = e->cosmology;
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
+  const double time_base = e->time_base;
+  const integertime_t t_current = e->ti_current;
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
+#endif
 
   TIMER_TIC;
 
@@ -892,6 +955,8 @@ void DOSELF_SUBSET(struct runner *r, struct cell *restrict ci,
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
         runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H);
+        runner_iact_nonsym_diffusion(r2, dx, hi, hj, pi, pj, a, H, time_base,
+                                     t_current, cosmo, with_cosmology);
 #endif
       }
     } /* loop over the parts in cj. */
@@ -935,6 +1000,11 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
 
   const struct engine *restrict e = r->e;
   const struct cosmology *restrict cosmo = e->cosmology;
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
+  const double time_base = e->time_base;
+  const integertime_t t_current = e->ti_current;
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
+#endif
 
   TIMER_TIC;
 
@@ -1060,6 +1130,8 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
           runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H);
+          runner_iact_nonsym_diffusion(r2, dx, hi, hj, pi, pj, a, H, time_base,
+                                       t_current, cosmo, with_cosmology);
 #endif
         }
       } /* loop over the parts in cj. */
@@ -1152,6 +1224,8 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
           runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H);
+          runner_iact_nonsym_diffusion(r2, dx, hj, hi, pj, pi, a, H, time_base,
+                                       t_current, cosmo, with_cosmology);
 #endif
         }
       } /* loop over the parts in ci. */
@@ -1268,6 +1342,11 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
 
   const struct engine *restrict e = r->e;
   const struct cosmology *restrict cosmo = e->cosmology;
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
+  const double time_base = e->time_base;
+  const integertime_t t_current = e->ti_current;
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
+#endif
 
   TIMER_TIC;
 
@@ -1458,6 +1537,8 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
           runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H);
+          runner_iact_nonsym_diffusion(r2, dx, hj, hi, pj, pi, a, H, time_base,
+                                       t_current, cosmo, with_cosmology);
 #endif
         }
       } /* loop over the active parts in cj. */
@@ -1533,6 +1614,8 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
             runner_iact_timebin(r2, dx, hi, hj, pi, pj, a, H);
+            runner_iact_diffusion(r2, dx, hi, hj, pi, pj, a, H, time_base,
+                                  t_current, cosmo, with_cosmology);
 #endif
           } else {
             IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
@@ -1543,6 +1626,9 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
             runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H);
+            runner_iact_nonsym_diffusion(r2, dx, hi, hj, pi, pj, a, H,
+                                         time_base, t_current, cosmo,
+                                         with_cosmology);
 #endif
           }
         }
@@ -1650,6 +1736,8 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
           runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H);
+          runner_iact_nonsym_diffusion(r2, dx, hi, hj, pi, pj, a, H, time_base,
+                                       t_current, cosmo, with_cosmology);
 #endif
         }
       } /* loop over the active parts in ci. */
@@ -1727,6 +1815,8 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
             runner_iact_timebin(r2, dx, hj, hi, pj, pi, a, H);
+            runner_iact_diffusion(r2, dx, hj, hi, pj, pi, a, H, time_base,
+                                  t_current, cosmo, with_cosmology);
 #endif
           } else {
             IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
@@ -1737,6 +1827,9 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
             runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H);
+            runner_iact_nonsym_diffusion(r2, dx, hj, hi, pj, pi, a, H,
+                                         time_base, t_current, cosmo,
+                                         with_cosmology);
 #endif
           }
         }
@@ -1856,6 +1949,11 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
 
   const struct engine *e = r->e;
   const struct cosmology *cosmo = e->cosmology;
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
+  const double time_base = e->time_base;
+  const integertime_t t_current = e->ti_current;
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
+#endif
 
   TIMER_TIC;
 
@@ -1930,6 +2028,8 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
           runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H);
+          runner_iact_nonsym_diffusion(r2, dx, hj, hi, pj, pi, a, H, time_base,
+                                       t_current, cosmo, with_cosmology);
 #endif
         }
       } /* loop over all other particles. */
@@ -1986,6 +2086,8 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
             runner_iact_timebin(r2, dx, hi, hj, pi, pj, a, H);
+            runner_iact_diffusion(r2, dx, hi, hj, pi, pj, a, H, time_base,
+                                  t_current, cosmo, with_cosmology);
 #endif
           } else if (doi) {
 
@@ -1997,6 +2099,9 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
             runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H);
+            runner_iact_nonsym_diffusion(r2, dx, hi, hj, pi, pj, a, H,
+                                         time_base, t_current, cosmo,
+                                         with_cosmology);
 #endif
           } else if (doj) {
 
@@ -2011,6 +2116,9 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
             runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H);
+            runner_iact_nonsym_diffusion(r2, dx, hj, hi, pj, pi, a, H,
+                                         time_base, t_current, cosmo,
+                                         with_cosmology);
 #endif
           }
         }
@@ -2068,6 +2176,11 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
 
   const struct engine *e = r->e;
   const struct cosmology *cosmo = e->cosmology;
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
+  const double time_base = e->time_base;
+  const integertime_t t_current = e->ti_current;
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
+#endif
 
   TIMER_TIC;
 
@@ -2142,6 +2255,8 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
           runner_iact_nonsym_timebin(r2, dx, hj, hi, pj, pi, a, H);
+          runner_iact_nonsym_diffusion(r2, dx, hj, hi, pj, pi, a, H, time_base,
+                                       t_current, cosmo, with_cosmology);
 #endif
         }
       } /* loop over all other particles. */
@@ -2193,6 +2308,8 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
             runner_iact_timebin(r2, dx, hi, hj, pi, pj, a, H);
+            runner_iact_diffusion(r2, dx, hi, hj, pi, pj, a, H, time_base,
+                                  t_current, cosmo, with_cosmology);
 #endif
           } else {
             IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
@@ -2203,6 +2320,9 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
 #endif
 #if (FUNCTION_TASK_LOOP == TASK_LOOP_FORCE)
             runner_iact_nonsym_timebin(r2, dx, hi, hj, pi, pj, a, H);
+            runner_iact_nonsym_diffusion(r2, dx, hi, hj, pi, pj, a, H,
+                                         time_base, t_current, cosmo,
+                                         with_cosmology);
 #endif
           }
         }
diff --git a/src/runner_others.c b/src/runner_others.c
index 7efb403213f15c7ad8a539cc30047d4b053e8912..f684b1714f4eaad4d5c8b9c87c162d7f9ef89b21 100644
--- a/src/runner_others.c
+++ b/src/runner_others.c
@@ -496,6 +496,7 @@ void runner_do_sink_formation(struct runner *r, struct cell *c) {
 void runner_do_end_hydro_force(struct runner *r, struct cell *c, int timer) {
 
   const struct engine *e = r->e;
+  const int with_cosmology = e->policy & engine_policy_cosmology;
 
   TIMER_TIC;
 
@@ -518,12 +519,24 @@ void runner_do_end_hydro_force(struct runner *r, struct cell *c, int timer) {
       /* Get a handle on the part. */
       struct part *restrict p = &parts[k];
 
+      double dt = 0;
       if (part_is_active(p, e)) {
 
+        if (with_cosmology) {
+          /* Compute the time step. */
+          const integertime_t ti_step = get_integer_timestep(p->time_bin);
+          const integertime_t ti_begin =
+              get_integer_time_begin(e->ti_current - 1, p->time_bin);
+
+          dt = cosmology_get_delta_time(cosmo, ti_begin, ti_begin + ti_step);
+        } else {
+          dt = get_timestep(p->time_bin, e->time_base);
+        }
+
         /* Finish the force loop */
         hydro_end_force(p, cosmo);
         timestep_limiter_end_force(p);
-        chemistry_end_force(p, cosmo);
+        chemistry_end_force(p, cosmo, with_cosmology, e->time, dt);
 
 #ifdef SWIFT_BOUNDARY_PARTICLES
 
diff --git a/tests/testSymmetry.c b/tests/testSymmetry.c
index a5a1be2ce21d42af6e86513fb7293a6376a351de..3a3359d9cc33c872758b8ddf9285859b0ac7e8da 100644
--- a/tests/testSymmetry.c
+++ b/tests/testSymmetry.c
@@ -46,6 +46,8 @@ void test(void) {
   /* Start with some values for the cosmological parameters */
   const float a = (float)random_uniform(0.8, 1.);
   const float H = 1.f;
+  const integertime_t ti_current = 1;
+  const double time_base = 1e-5;
 
   /* Create two random particles (don't do this at home !) */
   struct part pi, pj;
@@ -219,15 +221,23 @@ void test(void) {
 
   /* Call the symmetric version */
   runner_iact_force(r2, dx, pi.h, pj.h, &pi, &pj, a, H);
+  runner_iact_diffusion(r2, dx, pi.h, pj.h, &pi, &pj, a, H, time_base,
+                        ti_current, NULL, /*with_cosmology=*/0);
   runner_iact_timebin(r2, dx, pi.h, pj.h, &pi, &pj, a, H);
 
   /* Call the non-symmetric version */
   runner_iact_nonsym_force(r2, dx, pi2.h, pj2.h, &pi2, &pj2, a, H);
+  runner_iact_nonsym_diffusion(r2, dx, pi2.h, pj2.h, &pi2, &pj2, a, H,
+                               time_base, ti_current, NULL,
+                               /*with_cosmology=*/0);
   runner_iact_nonsym_timebin(r2, dx, pi2.h, pj2.h, &pi2, &pj2, a, H);
   dx[0] = -dx[0];
   dx[1] = -dx[1];
   dx[2] = -dx[2];
   runner_iact_nonsym_force(r2, dx, pj2.h, pi2.h, &pj2, &pi2, a, H);
+  runner_iact_nonsym_diffusion(r2, dx, pj2.h, pi2.h, &pj2, &pi2, a, H,
+                               time_base, ti_current, NULL,
+                               /*with_cosmology=*/0);
   runner_iact_nonsym_timebin(r2, dx, pj2.h, pi2.h, &pj2, &pi2, a, H);
 
 /* Check that the particles are the same */