diff --git a/examples/SmoothedMetallicity/getGlass.sh b/examples/SmoothedMetallicity/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..ffd92e88deae6e91237059adac2a6c2067caee46
--- /dev/null
+++ b/examples/SmoothedMetallicity/getGlass.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5
diff --git a/examples/SmoothedMetallicity/makeIC.py b/examples/SmoothedMetallicity/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..86679d5efe897b9dfae7db94b36d74bb047661e6
--- /dev/null
+++ b/examples/SmoothedMetallicity/makeIC.py
@@ -0,0 +1,110 @@
+###############################################################################
+# 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
+
+# Parameters
+gamma = 5./3.      # Gas adiabatic index
+rho0 = 1.          # Background density
+P0 = 1.e-6         # Background pressure
+Nelem = 9          # Gear: 9, EAGLE: 9
+low_metal = -6     # Low iron fraction
+high_metal = -5    # high iron fraction
+sigma_metal = 0.1  # 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)
+
+# ---------------------------------------------------
+glass = h5py.File("glassCube_32.hdf5", "r")
+
+# Read particle positions and h from the glass
+pos = glass["/PartType0/Coordinates"][:, :]
+h = glass["/PartType0/SmoothingLength"][:]
+
+numPart = np.size(h)
+vol = 1.
+
+# Generate extra arrays
+v = np.zeros((numPart, 3))
+ids = np.linspace(1, numPart, numPart)
+m = np.zeros(numPart)
+u = np.zeros(numPart)
+r = np.zeros(numPart)
+Z = 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)))
+nber = numPart - nber
+Z[np.logical_not(select), :] = high_metal * (1 + np.random.normal(
+    loc=0., scale=sigma_metal, size=(nber, Nelem)))
+
+# --------------------------------------------------
+
+# File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [1., 1., 1.]
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+grp.attrs["Dimension"] = 3
+
+# Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+# Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+# Particle group
+grp = file.create_group("/PartType0")
+grp.create_dataset('Coordinates', data=pos, dtype='d')
+grp.create_dataset('Velocities', data=v, dtype='f')
+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')
+
+
+file.close()
diff --git a/examples/SmoothedMetallicity/plotSolution.py b/examples/SmoothedMetallicity/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..e5bca3dfb7fe1e43c836733894c9e297cdd468ca
--- /dev/null
+++ b/examples/SmoothedMetallicity/plotSolution.py
@@ -0,0 +1,208 @@
+#!/usr/bin/env python
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+#                    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/>.
+# 
+##############################################################################
+
+# Computes the analytical solution of the 3D Smoothed Metallicity example.
+
+import h5py
+import sys
+import numpy as np
+import matplotlib
+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
+
+Nelem = 9
+# 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)
+
+# ---------------------------------------------------------------
+# Don't touch anything after this.
+# ---------------------------------------------------------------
+
+# Plot parameters
+params = {
+    'axes.labelsize': 10,
+    'axes.titlesize': 10,
+    'font.size': 12,
+    'legend.fontsize': 12,
+    'xtick.labelsize': 10,
+    'ytick.labelsize': 10,
+    'text.usetex': True,
+    'figure.figsize': (9.90, 6.45),
+    'figure.subplot.left': 0.045,
+    'figure.subplot.right': 0.99,
+    'figure.subplot.bottom': 0.05,
+    'figure.subplot.top': 0.99,
+    'figure.subplot.wspace': 0.15,
+    'figure.subplot.hspace': 0.12,
+    'lines.markersize': 6,
+    'lines.linewidth': 3.,
+    'text.latex.unicode': True
+}
+
+plt.rcParams.update(params)
+plt.rc('font', **{'family': 'sans-serif', 'sans-serif': ['Times']})
+
+
+snap = int(sys.argv[1])
+
+
+# Read the simulation data
+sim = h5py.File("smoothed_metallicity_%04d.hdf5" % snap, "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"]
+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"]
+git = sim["Code"].attrs["Git Revision"]
+
+pos = sim["/PartType0/Coordinates"][:, :]
+d = pos[:, 0] - boxSize / 2
+smooth_metal = sim["/PartType0/SmoothedElementAbundance"][:, :]
+metal = sim["/PartType0/ElementAbundance"][:, :]
+h = sim["/PartType0/SmoothingLength"][:]
+h = np.mean(h)
+
+if (Nelem != metal.shape[1]):
+    print("Unexpected number of element, please check makeIC.py"
+          " and plotSolution.py")
+    exit(1)
+
+N = 1000
+d_a = np.linspace(-boxSize / 2., boxSize / 2., N)
+
+# Now, work our the solution....
+
+
+def calc_a(d, high_metal, low_metal, std_dev, h):
+    """
+    Compute analytical solution:
+                        ___ High Metallicity
+    Low Metallicity ___/
+
+    where the linear part length is given by the smoothing length.
+
+    Parameters
+    ----------
+
+    d: np.array
+        Position to compute
+    high_metal: float
+        Value on the high metallicity side
+
+    low_metal: float
+        Value on the low metallicity side
+
+    std_dev: float
+        Standard deviation of the noise added
+
+    h: float
+        Mean smoothing length
+    """
+
+    # solution
+    a = np.zeros([len(d), Nelem])
+    # function at +- 1 sigma
+    sigma = np.zeros([len(d), Nelem, 2])
+
+    for i in range(Nelem):
+        # compute low metallicity side
+        s = d < -h
+        a[s, i] = low_metal[i]
+        # compute high metallicity side
+        s = d > h
+        a[s, i] = high_metal[i]
+
+        # compute non constant parts
+        m = (high_metal[i] - low_metal[i]) / (2.0 * h)
+        c = (high_metal[i] + low_metal[i]) / 2.0
+        # compute left linear part
+        s = d < - boxSize / 2.0 + h
+        a[s, i] = - m * (d[s] + boxSize / 2.0) + c
+        # compute middle linear part
+        s = np.logical_and(d >= -h, d <= h)
+        a[s, i] = m * d[s] + c
+        # compute right linear part
+        s = d > boxSize / 2.0 - h
+        a[s, i] = - m * (d[s] - boxSize / 2.0) + c
+
+    sigma[:, :, 0] = a * (1 + std_dev)
+    sigma[:, :, 1] = a * (1 - std_dev)
+    return a, sigma
+
+
+# compute solution
+sol, sig = calc_a(d_a, high_metal, low_metal, sigma_metal, h)
+
+# Plot the interesting quantities
+plt.figure()
+
+# Metallicity --------------------------------
+plt.subplot(221)
+for e in range(Nelem):
+    plt.plot(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)
+
+# Metallicity --------------------------------
+e = 0
+plt.subplot(223)
+plt.plot(d, smooth_metal[:, e], '.', color='r', ms=0.5, alpha=0.2)
+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.xlabel("${\\rm{Distance}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Smoothed~Metallicity}}~Z_\\textrm{sm}$", labelpad=0)
+plt.xlim(-0.5, 0.5)
+plt.ylim(low_metal[e]-1, high_metal[e]+1)
+
+# Information -------------------------------------
+plt.subplot(222, frameon=False)
+
+plt.text(-0.49, 0.9, "Smoothed Metallicity in 3D at $t=%.2f$" % time,
+         fontsize=10)
+plt.plot([-0.49, 0.1], [0.82, 0.82], 'k-', lw=1)
+plt.text(-0.49, 0.7, "$\\textsc{Swift}$ %s" % git, fontsize=10)
+plt.text(-0.49, 0.6, scheme, fontsize=10)
+plt.text(-0.49, 0.5, kernel, fontsize=10)
+plt.text(-0.49, 0.4, chemistry + "'s Chemistry", fontsize=10)
+plt.text(-0.49, 0.3, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+         fontsize=10)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
+plt.savefig("SmoothedMetallicity.png", dpi=200)
diff --git a/examples/SmoothedMetallicity/run.sh b/examples/SmoothedMetallicity/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..de8c55d678bcb611934af450940d8ed8e6c15d6b
--- /dev/null
+++ b/examples/SmoothedMetallicity/run.sh
@@ -0,0 +1,19 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e glassCube_32.hdf5 ]
+then
+    echo "Fetching initial glass file for the SmoothedMetallicity example..."
+    ./getGlass.sh
+fi
+if [ ! -e smoothed_metallicity.hdf5 ]
+then
+    echo "Generating initial conditions for the SmoothedMetallicity example..."
+    python makeIC.py
+fi
+
+# Run SWIFT
+../swift -n 1 -s -t 4 smoothed_metallicity.yml 2>&1 | tee output.log
+
+# Plot the solution
+python plotSolution.py 1
diff --git a/examples/SmoothedMetallicity/smoothed_metallicity.yml b/examples/SmoothedMetallicity/smoothed_metallicity.yml
new file mode 100644
index 0000000000000000000000000000000000000000..2e37695392b12c545bbbdbe7fd94748d5b3b9ff8
--- /dev/null
+++ b/examples/SmoothedMetallicity/smoothed_metallicity.yml
@@ -0,0 +1,34 @@
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1   # Grams
+  UnitLength_in_cgs:   1   # Centimeters
+  UnitVelocity_in_cgs: 1   # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   2e-4  # The end time of the simulation (in internal units).
+  dt_min:     1e-7  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-4  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            smoothed_metallicity # Common part of the name of output files
+  time_first:          0.    # Time of the first output (in internal units)
+  delta_time:          5e-5  # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-5 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./smoothed_metallicity.hdf5          # The file to read
+
diff --git a/src/Makefile.am b/src/Makefile.am
index 63d888f6c20b77a4922461a8063d65416fa1dd2e..0c8b53de3cca4ad4836834d13c9710ccce0cce00 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -122,12 +122,15 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
                  chemistry/none/chemistry.h \
 		 chemistry/none/chemistry_io.h \
 		 chemistry/none/chemistry_struct.h \
+		 chemistry/none/chemistry_iact.h \
                  chemistry/gear/chemistry.h \
 		 chemistry/gear/chemistry_io.h \
 		 chemistry/gear/chemistry_struct.h \
+		 chemistry/gear/chemistry_iact.h \
                  chemistry/EAGLE/chemistry.h \
 		 chemistry/EAGLE/chemistry_io.h \
-		 chemistry/EAGLE/chemistry_struct.h
+		 chemistry/EAGLE/chemistry_struct.h\
+		 chemistry/EAGLE/chemistry_iact.h
 
 
 # Sources and flags for regular library
diff --git a/src/cell.c b/src/cell.c
index d04535ac3500e917cc4035d2131fcd5e65f51fd4..b415cc077b569c84c85462634458e68aa5e4833a 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -49,6 +49,7 @@
 /* Local headers. */
 #include "active.h"
 #include "atomic.h"
+#include "chemistry.h"
 #include "drift.h"
 #include "engine.h"
 #include "error.h"
@@ -2440,6 +2441,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
       /* Get ready for a density calculation */
       if (part_is_active(p, e)) {
         hydro_init_part(p, &e->s->hs);
+        chemistry_init_part(p, e->chemistry);
       }
     }
 
diff --git a/src/chemistry.h b/src/chemistry.h
index b4571e61e80004bdb64e088394f409219c18ea5c..a5cbd77efbaab99e98ca5f031512e7e5bef0a613 100644
--- a/src/chemistry.h
+++ b/src/chemistry.h
@@ -31,10 +31,13 @@
 /* Import the right chemistry definition */
 #if defined(CHEMISTRY_NONE)
 #include "./chemistry/none/chemistry.h"
+#include "./chemistry/none/chemistry_iact.h"
 #elif defined(CHEMISTRY_GEAR)
 #include "./chemistry/gear/chemistry.h"
+#include "./chemistry/gear/chemistry_iact.h"
 #elif defined(CHEMISTRY_EAGLE)
 #include "./chemistry/EAGLE/chemistry.h"
+#include "./chemistry/EAGLE/chemistry_iact.h"
 #else
 #error "Invalid choice of chemistry function."
 #endif
diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h
index df0a4599a0f81eaf619b078b86c0ade4e1eba7b8..0ee4d2a2bfc54b1fe4b9b64f321cc23862841952 100644
--- a/src/chemistry/EAGLE/chemistry.h
+++ b/src/chemistry/EAGLE/chemistry.h
@@ -50,6 +50,32 @@ chemistry_get_element_name(enum chemistry_element elem) {
   return chemistry_element_names[elem];
 }
 
+/**
+ * @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_data containing chemistry informations.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_init_part(
+    struct part* restrict p, const struct chemistry_data* cd) {}
+
+/**
+ * @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_data containing chemistry informations.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_end_density(
+    struct part* restrict p, const struct chemistry_data* cd) {}
+
 /**
  * @brief Sets the chemistry properties of the (x-)particles to a valid start
  * state.
diff --git a/src/chemistry/EAGLE/chemistry_iact.h b/src/chemistry/EAGLE/chemistry_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..a00938e3557e68b0007e1464be2f05baefaf554a
--- /dev/null
+++ b/src/chemistry/EAGLE/chemistry_iact.h
@@ -0,0 +1,59 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_EAGLE_CHEMISTRY_IACT_H
+#define SWIFT_EAGLE_CHEMISTRY_IACT_H
+
+/**
+ * @file EAGLE/chemistry_iact.h
+ * @brief Smooth metal interaction functions following the EAGLE model.
+ */
+
+#include "chemistry_struct.h"
+
+/**
+ * @brief Do chemistry computation after the runner_iact_density (symmetric
+ * version)
+ *
+ * @param r2 Distance squared between particles
+ * @param dx Distance between particles
+ * @param hi Smoothing length of i
+ * @param hj Smoothing length of j
+ * @param pi #part i
+ * @param pj #part j
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_chemistry(
+    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+}
+
+/**
+ * @brief Do chemistry computation after the runner_iact_density (non-symmetric
+ * version)
+ *
+ * @param r2 Distance squared between particles
+ * @param dx Distance between particles
+ * @param hi Smoothing length of i
+ * @param hj Smoothing length of j
+ * @param pi #part i
+ * @param pj #part j
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
+    float r2, float *dx, float hi, float hj, struct part *pi,
+    const struct part *pj) {}
+
+#endif /* SWIFT_EAGLE_CHEMISTRY_IACT_H */
diff --git a/src/chemistry/gear/chemistry.h b/src/chemistry/gear/chemistry.h
index 0718cfd59fd2e86c4eb968241d5418dfd52045d9..5dc7a48939cdb298233abbe358473f59fc7e72f4 100644
--- a/src/chemistry/gear/chemistry.h
+++ b/src/chemistry/gear/chemistry.h
@@ -38,18 +38,17 @@
 #include "units.h"
 
 /**
- * @brief Sets the chemistry properties of the (x-)particles to a valid start
- * state.
- *
- * Nothing to do here.
- *
- * @param p Pointer to the particle data.
- * @param xp Pointer to the extended particle data.
- * @param data The global chemistry information.
+ * @brief Return a string containing the name of a given #chemistry_element.
  */
-__attribute__((always_inline)) INLINE static void chemistry_first_init_part(
-    const struct part* restrict p, struct xpart* restrict xp,
-    const struct chemistry_data* data) {}
+__attribute__((always_inline)) INLINE static const char*
+chemistry_get_element_name(enum chemistry_element elem) {
+
+  static const char* chemistry_element_names[chemistry_element_count] = {
+      "Oxygen",    "Magnesium", "Sulfur", "Iron",    "Zinc",
+      "Strontium", "Yttrium",   "Barium", "Europium"};
+
+  return chemistry_element_names[elem];
+}
 
 /**
  * @brief Initialises the chemistry properties.
@@ -72,7 +71,74 @@ static INLINE void chemistry_init_backend(
  */
 static INLINE void chemistry_print_backend(const struct chemistry_data* data) {
 
-  message("Chemistry function is 'gear'.");
+  message("Chemistry function is 'Gear'.");
+}
+
+/**
+ * @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_data containing chemistry informations.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_init_part(
+    struct part* restrict p, const struct chemistry_data* cd) {
+
+  struct chemistry_part_data* cpd = &p->chemistry_data;
+
+  for (int i = 0; i < chemistry_element_count; i++) {
+    cpd->smoothed_metal_mass_fraction[i] = 0.f;
+  }
+}
+
+/**
+ * @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_data containing chemistry informations.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_end_density(
+    struct part* restrict p, const struct chemistry_data* cd) {
+
+  /* Some smoothing length multiples. */
+  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 < 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;
+
+    /* Finish the calculation by inserting the missing h-factors */
+    cpd->smoothed_metal_mass_fraction[i] *= factor;
+  }
+}
+
+/**
+ * @brief Sets the chemistry properties of the (x-)particles to a valid start
+ * state.
+ *
+ * Nothing to do here.
+ *
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ * @param data The global chemistry information.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_first_init_part(
+    struct part* restrict p, struct xpart* restrict xp,
+    const struct chemistry_data* data) {
+  chemistry_init_part(p, data);
 }
 
 #endif /* SWIFT_CHEMISTRY_GEAR_H */
diff --git a/src/chemistry/gear/chemistry_iact.h b/src/chemistry/gear/chemistry_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..269d33e35fc7a2ba4f93c018c0c942fe99fb944c
--- /dev/null
+++ b/src/chemistry/gear/chemistry_iact.h
@@ -0,0 +1,115 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GEAR_CHEMISTRY_IACT_H
+#define SWIFT_GEAR_CHEMISTRY_IACT_H
+
+/**
+ * @file GEAR/chemistry_iact.h
+ * @brief Smooth metal interaction functions following the GEAR version of
+ * smooth metalicity.
+ *
+ * The interactions computed here are the ones presented in Wiersma, Schaye et
+ * al. 2009
+ */
+
+#include "chemistry_struct.h"
+
+/**
+ * @brief do chemistry computation after the runner_iact_density (symmetric
+ * version)
+ *
+ * @param r2 Distance squared between particles
+ * @param dx Distance between particles
+ * @param hi Smoothing length of i
+ * @param hj Smoothing length of j
+ * @param pi #part i
+ * @param pj #part j
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_chemistry(
+    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+  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;
+
+  /* Get r */
+  const float r = sqrtf(r2);
+
+  /* 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);
+
+  /* Compute contribution to the smooth metallicity */
+  for (int i = 0; i < 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;
+  }
+}
+
+/**
+ * @brief do chemistry computation after the runner_iact_density (non symmetric
+ * version)
+ *
+ * @param r2 Distance squared between particles
+ * @param dx Distance between particles
+ * @param hi Smoothing length of i
+ * @param hj Smoothing length of j
+ * @param pi #part i
+ * @param pj #part j
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
+    float r2, float *dx, float hi, float hj, struct part *pi,
+    const struct part *pj) {
+
+  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;
+
+  /* Get r */
+  const float r = sqrtf(r2);
+
+  /* Compute the kernel function for pi */
+  const float ui = r / hi;
+  kernel_deval(ui, &wi, &wi_dx);
+
+  /* Compute contribution to the smooth metallicity */
+  for (int i = 0; i < chemistry_element_count; i++) {
+    chi->smoothed_metal_mass_fraction[i] +=
+        mj * chj->metal_mass_fraction[i] * wi;
+  }
+}
+
+#endif /* SWIFT_GEAR_CHEMISTRY_IACT_H */
diff --git a/src/chemistry/gear/chemistry_io.h b/src/chemistry/gear/chemistry_io.h
index 8a6cb2279ce444d47bf732b8f89b5f161dbb6994..81ee8d59f58eef67097b38dec6ce04f7319b04d1 100644
--- a/src/chemistry/gear/chemistry_io.h
+++ b/src/chemistry/gear/chemistry_io.h
@@ -19,6 +19,8 @@
 #ifndef SWIFT_CHEMISTRY_IO_GEAR_H
 #define SWIFT_CHEMISTRY_IO_GEAR_H
 
+#include "chemistry.h"
+#include "chemistry_struct.h"
 #include "io_properties.h"
 
 /**
@@ -32,10 +34,9 @@
 int chemistry_read_particles(struct part* parts, struct io_props* list) {
 
   /* List what we want to read */
-  list[0] =
-      io_make_input_field("HeDensity", FLOAT, 1, COMPULSORY, UNIT_CONV_DENSITY,
-                          parts, chemistry_data.he_density);
-
+  list[0] = io_make_input_field(
+      "ElementAbundance", FLOAT, chemistry_element_count, OPTIONAL,
+      UNIT_CONV_NO_UNITS, parts, chemistry_data.metal_mass_fraction);
   return 1;
 }
 
@@ -50,19 +51,29 @@ int chemistry_read_particles(struct part* parts, struct io_props* list) {
 int chemistry_write_particles(const struct part* parts, struct io_props* list) {
 
   /* List what we want to write */
-  list[0] = io_make_output_field("HeDensity", FLOAT, 1, UNIT_CONV_DENSITY,
-                                 parts, chemistry_data.he_density);
+  list[0] = io_make_output_field(
+      "SmoothedElementAbundance", FLOAT, chemistry_element_count,
+      UNIT_CONV_NO_UNITS, parts, chemistry_data.smoothed_metal_mass_fraction);
 
-  return 1;
+  list[1] = io_make_output_field("ElementAbundance", FLOAT,
+                                 chemistry_element_count, UNIT_CONV_NO_UNITS,
+                                 parts, chemistry_data.metal_mass_fraction);
+
+  return 2;
 }
 
 /**
  * @brief Writes the current model of SPH to the file
- * @param h_grpsph The HDF5 group in which to write
+ * @param h_grp The HDF5 group in which to write
  */
-void chemistry_write_flavour(hid_t h_grpsph) {
+void chemistry_write_flavour(hid_t h_grp) {
 
-  io_write_attribute_s(h_grpsph, "Chemistry Model", "GEAR");
+  io_write_attribute_s(h_grp, "Chemistry Model", "GEAR");
+  for (size_t i = 0; i < chemistry_element_count; i++) {
+    char buffer[20];
+    sprintf(buffer, "Element %lu", i);
+    io_write_attribute_s(h_grp, buffer, chemistry_get_element_name(i));
+  }
 }
 
 #endif /* SWIFT_CHEMISTRY_IO_GEAR_H */
diff --git a/src/chemistry/gear/chemistry_struct.h b/src/chemistry/gear/chemistry_struct.h
index 401a68cdbf09cc59028c66d10833e19a95b5a297..b0ddee2b37c2d1126523ff7a4283f5c18ae3a7b4 100644
--- a/src/chemistry/gear/chemistry_struct.h
+++ b/src/chemistry/gear/chemistry_struct.h
@@ -22,7 +22,18 @@
 /**
  * @brief The individual elements traced in the model.
  */
-enum chemistry_element { chemistry_element_count = 0 };
+enum chemistry_element {
+  chemistry_element_O = 0,
+  chemistry_element_Mg,
+  chemistry_element_S,
+  chemistry_element_Fe,
+  chemistry_element_Zn,
+  chemistry_element_Sr,
+  chemistry_element_Y,
+  chemistry_element_Ba,
+  chemistry_element_Eu,
+  chemistry_element_count
+};
 
 /**
  * @brief Global chemical abundance information.
@@ -33,7 +44,12 @@ struct chemistry_data {};
  * @brief Properties of the chemistry function.
  */
 struct chemistry_part_data {
-  float he_density;
+
+  /*! Fraction of the particle mass in a given element */
+  float metal_mass_fraction[chemistry_element_count];
+
+  /*! Smoothed fraction of the particle mass in a given element */
+  float smoothed_metal_mass_fraction[chemistry_element_count];
 };
 
 #endif /* SWIFT_CHEMISTRY_STRUCT_GEAR_H */
diff --git a/src/chemistry/none/chemistry.h b/src/chemistry/none/chemistry.h
index ec183f7701dd6da0bcea12eff23cad2ee30ab97e..f6f826d580502f3c65bf56480ef0c04ddff1ee37 100644
--- a/src/chemistry/none/chemistry.h
+++ b/src/chemistry/none/chemistry.h
@@ -38,18 +38,15 @@
 #include "units.h"
 
 /**
- * @brief Sets the chemistry properties of the (x-)particles to a valid start
- * state.
- *
- * Nothing to do here.
- *
- * @param p Pointer to the particle data.
- * @param xp Pointer to the extended particle data.
- * @param data The global chemistry information used for this run.
+ * @brief Return a string containing the name of a given #chemistry_element.
  */
-__attribute__((always_inline)) INLINE static void chemistry_first_init_part(
-    const struct part* restrict p, struct xpart* restrict xp,
-    const struct chemistry_data* data) {}
+__attribute__((always_inline)) INLINE static const char*
+chemistry_get_element_name(enum chemistry_element elem) {
+
+  static const char* chemistry_element_names[chemistry_element_count] = {};
+
+  return chemistry_element_names[elem];
+}
 
 /**
  * @brief Initialises the chemistry properties.
@@ -75,4 +72,39 @@ static INLINE void chemistry_print_backend(const struct chemistry_data* data) {
   message("Chemistry function is 'No chemistry'.");
 }
 
+/**
+ * @brief Finishes the density calculation.
+ *
+ * @param p The particle to act upon
+ * @param cd The global chemistry information.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_end_density(
+    struct part* restrict p, const struct chemistry_data* cd) {}
+
+/**
+ * @brief Sets the chemistry properties of the (x-)particles to a valid start
+ * state.
+ *
+ * Nothing to do here.
+ *
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ * @param data The global chemistry information used for this run.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_first_init_part(
+    const struct part* restrict p, struct xpart* restrict xp,
+    const struct chemistry_data* data) {}
+
+/**
+ * @brief Sets the chemistry properties of the (x-)particles to a valid start
+ * state.
+ *
+ * Nothing to do here.
+ *
+ * @param p Pointer to the particle data.
+ * @param data The global chemistry information.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_init_part(
+    struct part* restrict p, const struct chemistry_data* data) {}
+
 #endif /* SWIFT_CHEMISTRY_NONE_H */
diff --git a/src/chemistry/none/chemistry_iact.h b/src/chemistry/none/chemistry_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..d963e4310e529f1a0f3202bc1f3ee68f605e811e
--- /dev/null
+++ b/src/chemistry/none/chemistry_iact.h
@@ -0,0 +1,62 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_NONE_CHEMISTRY_IACT_H
+#define SWIFT_NONE_CHEMISTRY_IACT_H
+
+/**
+ * @file none/chemistry_iact.h
+ * @brief Density computation
+ */
+
+#include "cache.h"
+#include "chemistry_struct.h"
+#include "minmax.h"
+
+/**
+ * @brief do chemistry computation after the runner_iact_density (symmetric
+ * version)
+ *
+ * @param r2 Distance squared between particles
+ * @param dx Distance between particles
+ * @param hi Smoothing length of i
+ * @param hj Smoothing length of j
+ * @param pi #part i
+ * @param pj #part j
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_chemistry(
+    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+}
+
+/**
+ * @brief do chemistry computation after the runner_iact_density (non symmetric
+ * version)
+ *
+ * @param r2 Distance squared between particles
+ * @param dx Distance between particles
+ * @param hi Smoothing length of i
+ * @param hj Smoothing length of j
+ * @param pi #part i
+ * @param pj #part j
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
+    float r2, float *dx, float hi, float hj, struct part *pi,
+    const struct part *pj) {}
+
+#endif /* SWIFT_NONE_CHEMISTRY_IACT_H */
diff --git a/src/runner.c b/src/runner.c
index 74367a7fa321a38e7e52001453b8c5af151ed2db..4dbadd278bac5d8ad9a3c340409ba62d15c55f19 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -42,6 +42,7 @@
 #include "approx_math.h"
 #include "atomic.h"
 #include "cell.h"
+#include "chemistry.h"
 #include "const.h"
 #include "cooling.h"
 #include "debug.h"
@@ -716,6 +717,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
 
           /* Finish the density calculation */
           hydro_end_density(p, cosmo);
+          chemistry_end_density(p, e->chemistry);
 
           /* Compute one step of the Newton-Raphson scheme */
           const float n_sum = p->density.wcount * h_old_dim;
@@ -753,6 +755,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
 
             /* Re-initialise everything */
             hydro_init_part(p, hs);
+            chemistry_init_part(p, e->chemistry);
 
             /* Off we go ! */
             continue;
diff --git a/src/runner_doiact.h b/src/runner_doiact.h
index 4303fbc235c55ac5a9ab977011c1dcc198d0d59d..7ed6db614a5ca5a390b58498c539b8e87fffa558 100644
--- a/src/runner_doiact.h
+++ b/src/runner_doiact.h
@@ -202,6 +202,9 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci,
       if (r2 < hig2 && pi_active) {
 
         IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+        runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#endif
       }
       if (r2 < hjg2 && pj_active) {
 
@@ -210,6 +213,9 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci,
         dx[2] = -dx[2];
 
         IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+        runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#endif
       }
 
     } /* loop over the parts in cj. */
@@ -297,9 +303,15 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
         if (pi_active && pj_active) {
 
           IACT(r2, dx, hi, hj, pi, pj, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+          runner_iact_chemistry(r2, dx, hi, hj, pi, pj);
+#endif
         } else if (pi_active) {
 
           IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+          runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#endif
         } else if (pj_active) {
 
           dx[0] = -dx[0];
@@ -307,6 +319,9 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
           dx[2] = -dx[2];
 
           IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+          runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#endif
         }
       }
     } /* loop over the parts in cj. */
@@ -381,9 +396,15 @@ void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) {
       if (doi && doj) {
 
         IACT(r2, dx, hi, hj, pi, pj, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+        runner_iact_chemistry(r2, dx, hi, hj, pi, pj);
+#endif
       } else if (doi) {
 
         IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+        runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#endif
       } else if (doj) {
 
         dx[0] = -dx[0];
@@ -391,6 +412,9 @@ void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) {
         dx[2] = -dx[2];
 
         IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+        runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#endif
       }
     } /* loop over the parts in cj. */
   }   /* loop over the parts in ci. */
@@ -464,9 +488,15 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) {
       if (doi && doj) {
 
         IACT(r2, dx, hi, hj, pi, pj, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+        runner_iact_chemistry(r2, dx, hi, hj, pi, pj);
+#endif
       } else if (doi) {
 
         IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+        runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#endif
       } else if (doj) {
 
         dx[0] = -dx[0];
@@ -474,6 +504,9 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) {
         dx[2] = -dx[2];
 
         IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+        runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#endif
       }
     } /* loop over the parts in cj. */
   }   /* loop over the parts in ci. */
@@ -553,6 +586,9 @@ void DOPAIR_SUBSET_NAIVE(struct runner *r, struct cell *restrict ci,
       if (r2 < hig2) {
 
         IACT_NONSYM(r2, dx, hi, pj->h, pi, pj, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+        runner_iact_nonsym_chemistry(r2, dx, hi, pj->h, pi, pj);
+#endif
       }
     } /* loop over the parts in cj. */
   }   /* loop over the parts in ci. */
@@ -637,6 +673,9 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
         if (r2 < hig2) {
 
           IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+          runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#endif
         }
       } /* loop over the parts in cj. */
     }   /* loop over the parts in ci. */
@@ -684,6 +723,9 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
         if (r2 < hig2) {
 
           IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+          runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#endif
         }
       } /* loop over the parts in cj. */
     }   /* loop over the parts in ci. */
@@ -813,6 +855,9 @@ void DOSELF_SUBSET(struct runner *r, struct cell *restrict ci,
       if (r2 > 0.f && r2 < hig2) {
 
         IACT_NONSYM(r2, dx, hi, pj->h, pi, pj, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+        runner_iact_nonsym_chemistry(r2, dx, hi, pj->h, pi, pj);
+#endif
       }
     } /* loop over the parts in cj. */
   }   /* loop over the parts in ci. */
@@ -966,6 +1011,9 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
         if (r2 < hig2) {
 
           IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+          runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#endif
         }
       } /* loop over the parts in cj. */
     }   /* loop over the parts in ci. */
@@ -1046,6 +1094,9 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
         if (r2 < hjg2) {
 
           IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+          runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#endif
         }
       } /* loop over the parts in ci. */
     }   /* loop over the parts in cj. */
@@ -1314,6 +1365,9 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
            (note that we will do the other condition in the reverse loop) */
         if (r2 < hig2) {
           IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+          runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#endif
         }
       } /* loop over the active parts in cj. */
     }
@@ -1374,10 +1428,17 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
         if (r2 < hig2) {
 
           /* Does pj need to be updated too? */
-          if (part_is_active(pj, e))
+          if (part_is_active(pj, e)) {
             IACT(r2, dx, hi, hj, pi, pj, a, H);
-          else
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+            runner_iact_chemistry(r2, dx, hi, hj, pi, pj);
+#endif
+          } else {
             IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+            runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#endif
+          }
         }
       } /* loop over the parts in cj. */
     }   /* Is pi active? */
@@ -1464,6 +1525,9 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
            (note that we must avoid the r2 < hig2 cases we already processed) */
         if (r2 < hjg2 && r2 >= hig2) {
           IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+          runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#endif
         }
       } /* loop over the active parts in ci. */
     }
@@ -1527,10 +1591,17 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
         if (r2 < hjg2 && r2 >= hig2) {
 
           /* Does pi need to be updated too? */
-          if (part_is_active(pi, e))
+          if (part_is_active(pi, e)) {
             IACT(r2, dx, hj, hi, pj, pi, a, H);
-          else
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+            runner_iact_chemistry(r2, dx, hj, hi, pj, pi);
+#endif
+          } else {
             IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+            runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#endif
+          }
         }
       } /* loop over the parts in ci. */
     }   /* Is pj active? */
@@ -1701,6 +1772,9 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
         if (r2 < hj * hj * kernel_gamma2) {
 
           IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+          runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#endif
         }
       } /* loop over all other particles. */
     }
@@ -1740,15 +1814,24 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
         if (r2 < hig2 || doj) {
 
           /* Which parts need to be updated? */
-          if (r2 < hig2 && doj)
+          if (r2 < hig2 && doj) {
             IACT(r2, dx, hi, hj, pi, pj, a, H);
-          else if (!doj)
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+            runner_iact_chemistry(r2, dx, hi, hj, pi, pj);
+#endif
+          } else if (!doj) {
             IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
-          else {
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+            runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#endif
+          } else {
             dx[0] = -dx[0];
             dx[1] = -dx[1];
             dx[2] = -dx[2];
             IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+            runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#endif
           }
         }
       } /* loop over all other particles. */
@@ -1866,6 +1949,9 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
         if (r2 < hig2 || r2 < hj * hj * kernel_gamma2) {
 
           IACT_NONSYM(r2, dx, hj, hi, pj, pi, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+          runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#endif
         }
       } /* loop over all other particles. */
     }
@@ -1903,10 +1989,17 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
         if (r2 < hig2 || r2 < hj * hj * kernel_gamma2) {
 
           /* Does pj need to be updated too? */
-          if (part_is_active(pj, e))
+          if (part_is_active(pj, e)) {
             IACT(r2, dx, hi, hj, pi, pj, a, H);
-          else
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+            runner_iact_chemistry(r2, dx, hi, hj, pi, pj);
+#endif
+          } else {
             IACT_NONSYM(r2, dx, hi, hj, pi, pj, a, H);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+            runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#endif
+          }
         }
       } /* loop over all other particles. */
     }
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 1033dc96cd93a203491499188778244903108252..c00933e587269d3743aaaf26da8515897fc6e70b 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -463,7 +463,10 @@ void runner_doself2_force_vec(struct runner *r, struct cell *ci);
 /* And go... */
 int main(int argc, char *argv[]) {
 
+#ifdef HAVE_SETAFFINITY
   engine_pin();
+#endif
+
   size_t runs = 0, particles = 0;
   double h = 1.23485, size = 1., rho = 2.5;
   double perturbation = 0.;
diff --git a/tests/test27cells.c b/tests/test27cells.c
index 97581527e552b9b6c7a8aa55a9ac8e3317af46fb..76fd18565907bd3eb1c738d2e2cb46680e926cfd 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -335,7 +335,10 @@ void runner_doself_subset_branch_density(struct runner *r,
 /* And go... */
 int main(int argc, char *argv[]) {
 
+#ifdef HAVE_SETAFFINITY
   engine_pin();
+#endif
+
   size_t runs = 0, particles = 0;
   double h = 1.23485, size = 1., rho = 1.;
   double perturbation = 0., h_pert = 0.;
diff --git a/tests/testPeriodicBC.c b/tests/testPeriodicBC.c
index 0012d77ade89e3d35046f639dc669cae59c138ca..541f1a99732dac5b74123c986f16caa3ba6f7aac 100644
--- a/tests/testPeriodicBC.c
+++ b/tests/testPeriodicBC.c
@@ -380,7 +380,10 @@ void test_boundary_conditions(struct cell **cells, struct runner runner,
 /* And go... */
 int main(int argc, char *argv[]) {
 
+#ifdef HAVE_SETAFFINITY
   engine_pin();
+#endif
+
   size_t runs = 0, particles = 0;
   double h = 1.23485, size = 1., rho = 1.;
   double perturbation = 0.;
diff --git a/tests/testSymmetry.c b/tests/testSymmetry.c
index 1ba5080cc8c265f59c4ab1a599d19bb9b4803fb0..868fec6019ada071d55d82a505f2619010257cf2 100644
--- a/tests/testSymmetry.c
+++ b/tests/testSymmetry.c
@@ -26,12 +26,13 @@
 
 #include "swift.h"
 
-int main(int argc, char *argv[]) {
+void print_bytes(void *p, size_t len) {
+  printf("(");
+  for (size_t i = 0; i < len; ++i) printf("%02x", ((unsigned char *)p)[i]);
+  printf(")\n");
+}
 
-/* Choke on FPEs */
-#ifdef HAVE_FE_ENABLE_EXCEPT
-  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
-#endif
+void test() {
 
 #if defined(SHADOWFAX_SPH)
   /* Initialize the Voronoi simulation box */
@@ -56,8 +57,8 @@ int main(int argc, char *argv[]) {
   for (size_t i = 0; i < 3; ++i) pj.x[0] = random_uniform(-1., 1.);
   pi.h = 2.f;
   pj.h = 2.f;
-  pi.id = 1;
-  pj.id = 2;
+  pi.id = 1ll;
+  pj.id = 2ll;
 
 #if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
   /* Give the primitive variables sensible values, since the Riemann solver does
@@ -124,11 +125,11 @@ int main(int argc, char *argv[]) {
   memcpy(&pi2, &pi, sizeof(struct part));
   memcpy(&pj2, &pj, sizeof(struct part));
 
-  int i_ok = memcmp(&pi, &pi2, sizeof(struct part));
-  int j_ok = memcmp(&pj, &pj2, sizeof(struct part));
+  int i_not_ok = memcmp(&pi, &pi2, sizeof(struct part));
+  int j_not_ok = memcmp(&pj, &pj2, sizeof(struct part));
 
-  if (i_ok != 0) error("Particles 'pi' do not match after copy");
-  if (j_ok != 0) error("Particles 'pj' do not match after copy");
+  if (i_not_ok) error("Particles 'pi' do not match after copy");
+  if (j_not_ok) error("Particles 'pj' do not match after copy");
 
   /* Compute distance vector */
   float dx[3];
@@ -141,20 +142,23 @@ int main(int argc, char *argv[]) {
 
   /* Call the symmetric version */
   runner_iact_density(r2, dx, pi.h, pj.h, &pi, &pj, a, H);
+  runner_iact_chemistry(r2, dx, pi.h, pj.h, &pi, &pj);
 
   /* Call the non-symmetric version */
   runner_iact_nonsym_density(r2, dx, pi2.h, pj2.h, &pi2, &pj2, a, H);
+  runner_iact_nonsym_chemistry(r2, dx, pi2.h, pj2.h, &pi2, &pj2);
   dx[0] = -dx[0];
   dx[1] = -dx[1];
   dx[2] = -dx[2];
   runner_iact_nonsym_density(r2, dx, pj2.h, pi2.h, &pj2, &pi2, a, H);
+  runner_iact_nonsym_chemistry(r2, dx, pj2.h, pi2.h, &pj2, &pi2);
 
   /* Check that the particles are the same */
-  i_ok = memcmp(&pi, &pi2, sizeof(struct part));
-  j_ok = memcmp(&pj, &pj2, sizeof(struct part));
+  i_not_ok = memcmp(&pi, &pi2, sizeof(struct part));
+  j_not_ok = memcmp(&pj, &pj2, sizeof(struct part));
 
-  if (i_ok) error("Particles 'pi' do not match after density");
-  if (j_ok) error("Particles 'pj' do not match after density");
+  if (i_not_ok) error("Particles 'pi' do not match after density");
+  if (j_not_ok) error("Particles 'pj' do not match after density");
 
   /* --- Test the force loop --- */
 
@@ -170,8 +174,8 @@ int main(int argc, char *argv[]) {
 
 /* Check that the particles are the same */
 #if defined(GIZMO_SPH)
-  i_ok = 0;
-  j_ok = 0;
+  i_not_ok = 0;
+  j_not_ok = 0;
   for (size_t i = 0; i < sizeof(struct part) / sizeof(float); ++i) {
     float a = *(((float *)&pi) + i);
     float b = *(((float *)&pi2) + i);
@@ -198,24 +202,53 @@ int main(int argc, char *argv[]) {
       message("%.8e, %.8e, %lu", c, d, i);
     }
 
-    i_ok |= a_is_b;
-    j_ok |= c_is_d;
+    i_not_ok |= a_is_b;
+    j_not_ok |= c_is_d;
   }
 #else
-  i_ok = memcmp(&pi, &pi2, sizeof(struct part));
-  j_ok = memcmp(&pj, &pj2, sizeof(struct part));
+  i_not_ok =
+      strncmp((const char *)&pi, (const char *)&pi2, sizeof(struct part));
+  j_not_ok =
+      strncmp((const char *)&pj, (const char *)&pj2, sizeof(struct part));
 #endif
 
-  if (i_ok) {
+  if (i_not_ok) {
     printParticle_single(&pi, &xpi);
     printParticle_single(&pi2, &xpi);
-    error("Particles 'pi' do not match after force");
+    print_bytes(&pj, sizeof(struct part));
+    print_bytes(&pj2, sizeof(struct part));
+    error("Particles 'pi' do not match after force (byte = %d)", i_not_ok);
   }
-  if (j_ok) {
+  if (j_not_ok) {
     printParticle_single(&pj, &xpj);
     printParticle_single(&pj2, &xpj);
-    error("Particles 'pj' do not match after force");
+    print_bytes(&pj, sizeof(struct part));
+    print_bytes(&pj2, sizeof(struct part));
+    error("Particles 'pj' do not match after force (byte = %d)", j_not_ok);
+  }
+}
+
+int main(int argc, char *argv[]) {
+
+  /* Initialize CPU frequency, this also starts time. */
+  unsigned long long cpufreq = 0;
+  clocks_set_cpufreq(cpufreq);
+
+/* Choke on FPEs */
+#ifdef HAVE_FE_ENABLE_EXCEPT
+  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
+#endif
+
+  /* Get some randomness going */
+  const int seed = time(NULL);
+  message("Seed = %d", seed);
+  srand(seed);
+
+  for (int i = 0; i < 100; ++i) {
+    message("Random test %d/100", i);
+    test();
   }
+  message("All good");
 
   return 0;
 }