diff --git a/src/Makefile.am b/src/Makefile.am
index 82b7aa73e10c1fa0b369886c48322985aadbc618..7b216af3ad083379f182fb6e8d887b4767837bc5 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/chemistry.h b/src/chemistry.h
index b7a17bdb313f3b5ddc2416cdc7a729e4f8916ffe..96745c805a77f1ea384179d18a94bff84a642c17 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/gear/chemistry_iact.h b/src/chemistry/gear/chemistry_iact.h
index 076afe86ef8229617a208fc35ad5505cfaba00f6..ea56fd8ad52ead50e92f726bdaee4dddd404ac0c 100644
--- a/src/chemistry/gear/chemistry_iact.h
+++ b/src/chemistry/gear/chemistry_iact.h
@@ -36,10 +36,10 @@
  */
 __attribute__((always_inline)) INLINE static void runner_iact_chemistry(
 float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj,
-const chemistry_data chem_data) {
+const struct chemistry_data *chem_data) {
 
-  struct chemistry_part_data *chi = &pi->Chemistry_data;
-  struct chemistry_part_data *chj = &pj->Chemistry_data;
+  struct chemistry_part_data *chi = &pi->chemistry_data;
+  struct chemistry_part_data *chj = &pj->chemistry_data;
 
   float wi, wi_dx;
   float wj, wj_dx;
@@ -61,8 +61,8 @@ const chemistry_data chem_data) {
 
   /* Compute contribution to the smooth metallicity */
   for(size_t i=0; i < chemistry_element_count; i++) {
-    chi.smoothed_metal_mass_fraction[i] += mj * chj.smoothed_metal_mass_fraction[i] * wi;
-    chj.smoothed_metal_mass_fraction[i] += mi * chi.smoothed_metal_mass_fraction[i] * wj;
+    chi->smoothed_metal_mass_fraction[i] += mj * chj->smoothed_metal_mass_fraction[i] * wi;
+    chj->smoothed_metal_mass_fraction[i] += mi * chi->smoothed_metal_mass_fraction[i] * wj;
   }
 }
 
@@ -71,15 +71,14 @@ const chemistry_data chem_data) {
  */
 __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,
-const chemistry_data chem_data) {
+const struct chemistry_data *chem_data) {
 
-  struct chemistry_part_data *chi = &pi->Chemistry_data;
-  struct chemistry_part_data *chj = &pj->Chemistry_data;
+  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 mi = pi->mass;
   const float mj = pj->mass;
 
   /* Get r */
@@ -91,7 +90,7 @@ const chemistry_data chem_data) {
 
   /* Compute contribution to the smooth metallicity */
   for(size_t i=0; i < chemistry_element_count; i++) {
-    chi.smoothed_metal_mass_fraction[i] += mj * chj.smoothed_metal_mass_fraction[i] * wi;
+    chi->smoothed_metal_mass_fraction[i] += mj * chj->smoothed_metal_mass_fraction[i] * wi;
   }
 }
 
diff --git a/src/chemistry/gear/chemistry_io.h b/src/chemistry/gear/chemistry_io.h
index 1e63faf2dd775ee8ea8090cd4b62966d905c1327..bd83659b6a83ef403e903a9f0546e36908fa2b50 100644
--- a/src/chemistry/gear/chemistry_io.h
+++ b/src/chemistry/gear/chemistry_io.h
@@ -59,8 +59,6 @@ int chemistry_write_particles(const struct part* parts, struct io_props* list) {
   list[0] = io_make_output_field("SmoothedElementAbundance", FLOAT, chemistry_element_count,
 				 UNIT_CONV_NO_UNITS,
 				 parts, chemistry_data.smoothed_metal_mass_fraction);
-  }
-
   return 1;
 }
 
diff --git a/src/runner.c b/src/runner.c
index 9561faa9a3efc1f7842f96e51f604d4967460748..23431e9964eab97957d98be0474daaa38398bac0 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -701,7 +701,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
 
           /* Finish the density calculation */
           hydro_end_density(p);
-	  chemistry_end(p);
+	  chemistry_end(p, e->chemistry);
 
           /* Compute one step of the Newton-Raphson scheme */
           const float n_sum = p->density.wcount * h_old_dim;
@@ -739,7 +739,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
 
             /* Re-initialise everything */
             hydro_init_part(p, &s->hs);
-	    chemistry_init_part(p, &e->chemistry);
+	    chemistry_init_part(p, e->chemistry);
 
             /* Off we go ! */
             continue;
diff --git a/src/runner_doiact.h b/src/runner_doiact.h
index 3c923bbdee20624bee9b2d9fb4df3a891b431deb..cebe480e855fc3df4b14fd4b3027a01a31a96510 100644
--- a/src/runner_doiact.h
+++ b/src/runner_doiact.h
@@ -24,6 +24,10 @@
    and runner_dosub_FUNCTION calling the pairwise interaction function
    runner_iact_FUNCTION. */
 
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+#include "chemistry.h"
+#endif
+
 #define PASTE(x, y) x##_##y
 
 #define _DOPAIR1_BRANCH(f) PASTE(runner_dopair1_branch, f)
@@ -197,8 +201,8 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci,
       if (r2 < hig2 && pi_active) {
 
         IACT_NONSYM(r2, dx, hi, hj, pi, pj);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, e->chemistry);
 #endif
       }
       if (r2 < hjg2 && pj_active) {
@@ -208,8 +212,8 @@ void DOPAIR1_NAIVE(struct runner *r, struct cell *restrict ci,
         dx[2] = -dx[2];
 
         IACT_NONSYM(r2, dx, hj, hi, pj, pi);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, e->chemistry);
 #endif
       }
 
@@ -293,14 +297,14 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
         if (pi_active && pj_active) {
 
           IACT(r2, dx, hi, hj, pi, pj);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	  runner_iact_chemistry(r2, dx, hi, hj, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	  runner_iact_chemistry(r2, dx, hi, hj, pi, pj, e->chemistry);
 #endif
         } else if (pi_active) {
 
           IACT_NONSYM(r2, dx, hi, hj, pi, pj);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	  runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	  runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, e->chemistry);
 #endif
         } else if (pj_active) {
 
@@ -309,8 +313,8 @@ void DOPAIR2_NAIVE(struct runner *r, struct cell *restrict ci,
           dx[2] = -dx[2];
 
           IACT_NONSYM(r2, dx, hj, hi, pj, pi);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	  runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	  runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, e->chemistry);
 #endif
         }
       }
@@ -381,14 +385,14 @@ void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) {
       if (doi && doj) {
 
         IACT(r2, dx, hi, hj, pi, pj);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	runner_iact_chemistry(r2, dx, hi, hj, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	runner_iact_chemistry(r2, dx, hi, hj, pi, pj, e->chemistry);
 #endif
       } else if (doi) {
 
         IACT_NONSYM(r2, dx, hi, hj, pi, pj);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, e->chemistry);
 #endif
       } else if (doj) {
 
@@ -397,8 +401,8 @@ void DOSELF1_NAIVE(struct runner *r, struct cell *restrict c) {
         dx[2] = -dx[2];
 
         IACT_NONSYM(r2, dx, hj, hi, pj, pi);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, e->chemistry);
 #endif
       }
     } /* loop over the parts in cj. */
@@ -468,14 +472,14 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) {
       if (doi && doj) {
 
         IACT(r2, dx, hi, hj, pi, pj);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	runner_iact_chemistry(r2, dx, hi, hj, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	runner_iact_chemistry(r2, dx, hi, hj, pi, pj, e->chemistry);
 #endif
       } else if (doi) {
 
         IACT_NONSYM(r2, dx, hi, hj, pi, pj);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, e->chemistry);
 #endif
       } else if (doj) {
 
@@ -484,8 +488,8 @@ void DOSELF2_NAIVE(struct runner *r, struct cell *restrict c) {
         dx[2] = -dx[2];
 
         IACT_NONSYM(r2, dx, hj, hi, pj, pi);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, e->chemistry);
 #endif
       }
     } /* loop over the parts in cj. */
@@ -560,8 +564,8 @@ void DOPAIR_SUBSET_NAIVE(struct runner *r, struct cell *restrict ci,
       if (r2 < hig2) {
 
         IACT_NONSYM(r2, dx, hi, pj->h, pi, pj);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	runner_iact_nonsym_chemistry(r2, dx, hi, pj->h, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	runner_iact_nonsym_chemistry(r2, dx, hi, pj->h, pi, pj, r->e->chemistry);
 #endif
       }
     } /* loop over the parts in cj. */
@@ -641,8 +645,8 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
         if (r2 < hig2) {
 
           IACT_NONSYM(r2, dx, hi, hj, pi, pj);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	  runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	  runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, r->e->chemistry);
 #endif
         }
       } /* loop over the parts in cj. */
@@ -692,8 +696,8 @@ void DOPAIR_SUBSET(struct runner *r, struct cell *restrict ci,
         if (r2 < hig2) {
 
           IACT_NONSYM(r2, dx, hi, hj, pi, pj);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, r->e->chemistry);
 #endif
         }
       } /* loop over the parts in cj. */
@@ -821,8 +825,8 @@ 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);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	runner_iact_nonsym_chemistry(r2, dx, hi, pj->h, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	runner_iact_nonsym_chemistry(r2, dx, hi, pj->h, pi, pj, r->e->chemistry);
 #endif
       }
     } /* loop over the parts in cj. */
@@ -972,8 +976,8 @@ 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);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	  runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, e->chemistry);
 #endif
         }
       } /* loop over the parts in cj. */
@@ -1055,8 +1059,8 @@ 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);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	  runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	  runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, e->chemistry);
 #endif
         }
       } /* loop over the parts in ci. */
@@ -1321,8 +1325,8 @@ 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);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	  runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	  runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, e->chemistry);
 #endif
         }
       } /* loop over the active parts in cj. */
@@ -1386,14 +1390,14 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
           /* Does pj need to be updated too? */
           if (part_is_active(pj, e)) {
 	    IACT(r2, dx, hi, hj, pi, pj);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	    runner_iact_chemistry(r2, dx, hi, hj, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	    runner_iact_chemistry(r2, dx, hi, hj, pi, pj, e->chemistry);
 #endif
 	  }
           else {
             IACT_NONSYM(r2, dx, hi, hj, pi, pj);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	    runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	    runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, e->chemistry);
 #endif
 	  }
         }
@@ -1482,8 +1486,8 @@ 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);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	  runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	  runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, e->chemistry);
 #endif
         }
       } /* loop over the active parts in ci. */
@@ -1550,14 +1554,14 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj, const int sid,
           /* Does pi need to be updated too? */
           if (part_is_active(pi, e)) {
             IACT(r2, dx, hj, hi, pj, pi);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	    runner_iact_chemistry(r2, dx, hj, hi, pj, pi);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	    runner_iact_chemistry(r2, dx, hj, hi, pj, pi, e->chemistry);
 #endif
 	  }
           else {
             IACT_NONSYM(r2, dx, hj, hi, pj, pi);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	    runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	    runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, e->chemistry);
 #endif
 	  }
         }
@@ -1725,8 +1729,8 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
         if (r2 < hj * hj * kernel_gamma2) {
 
           IACT_NONSYM(r2, dx, hj, hi, pj, pi);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	  runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	  runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, e->chemistry);
 #endif
         }
       } /* loop over all other particles. */
@@ -1769,14 +1773,14 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
           /* Which parts need to be updated? */
           if (r2 < hig2 && doj) {
             IACT(r2, dx, hi, hj, pi, pj);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	    runner_iact_chemistry(r2, dx, hi, hj, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	    runner_iact_chemistry(r2, dx, hi, hj, pi, pj, e->chemistry);
 #endif
 	  }
           else if (!doj) {
             IACT_NONSYM(r2, dx, hi, hj, pi, pj);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	    runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	    runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, e->chemistry);
 #endif
 	  }
           else {
@@ -1784,8 +1788,8 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
             dx[1] = -dx[1];
             dx[2] = -dx[2];
             IACT_NONSYM(r2, dx, hj, hi, pj, pi);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	    runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	    runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, e->chemistry);
 #endif
           }
         }
@@ -1899,8 +1903,8 @@ 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);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	  runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	  runner_iact_nonsym_chemistry(r2, dx, hj, hi, pj, pi, e->chemistry);
 #endif
         }
       } /* loop over all other particles. */
@@ -1941,14 +1945,14 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
           /* Does pj need to be updated too? */
           if (part_is_active(pj, e)) {
             IACT(r2, dx, hi, hj, pi, pj);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	    runner_iact_chemistry(r2, dx, hi, hj, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	    runner_iact_chemistry(r2, dx, hi, hj, pi, pj, e->chemistry);
 #endif
 	  }
           else {
             IACT_NONSYM(r2, dx, hi, hj, pi, pj);
-#ifdef (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
-	    runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj);
+#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
+	    runner_iact_nonsym_chemistry(r2, dx, hi, hj, pi, pj, e->chemistry);
 #endif
 	  }
         }