diff --git a/src/Makefile.am b/src/Makefile.am
index f1f9157b1881a673456f7dbf394aa32e09b35379..db86d524cecd7a263e2ddf32b2545ad00397b678 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -114,7 +114,7 @@ AM_SOURCES = space.c runner_main.c runner_doiact_hydro.c runner_doiact_limiter.c
 
 # Include files for distribution, not installation.
 nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
-		 gravity_iact.h kernel_long_gravity.h vector.h cache.h \
+		 gravity_iact.h kernel_long_gravity.h vector.h accumulate.h cache.h \
 	         runner_doiact_nosort.h runner_doiact_hydro.h runner_doiact_stars.h runner_doiact_black_holes.h runner_doiact_grav.h \
                  runner_doiact_functions_hydro.h runner_doiact_functions_stars.h runner_doiact_functions_black_holes.h \
 		 runner_doiact_functions_limiter.h runner_doiact_limiter.h units.h intrinsics.h minmax.h \
diff --git a/src/accumulate.h b/src/accumulate.h
new file mode 100644
index 0000000000000000000000000000000000000000..4968b331ef713fe42f89e15e1ec7371639018040
--- /dev/null
+++ b/src/accumulate.h
@@ -0,0 +1,92 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * 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_ACCUMULATE_H
+#define SWIFT_ACCUMULATE_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local includes */
+#include "atomic.h"
+
+/**
+ * @file accumulate.h
+ * @brief Defines a series of functions used to update the fields of a particle
+ * atomically. These functions should be used in any task that can run in
+ * parallel to another one.
+ */
+
+/**
+ * @brief Add x to the value stored at the location address (int version)
+ *
+ * When SWIFT_TASKS_WITHOUT_ATOMICS is *not* defined this function uses an
+ * atomic operation.
+ *
+ * @param address The address to update.
+ * @param x The value to add to *address.
+ */
+__attribute__((always_inline)) INLINE static void accumulate_add_i(
+    volatile int *const address, const int x) {
+
+#ifdef SWIFT_TASKS_WITHOUT_ATOMICS
+  *address += x;
+#else
+  atomic_add(address, x);
+#endif
+}
+
+/**
+ * @brief Add x to the value stored at the location address (float version)
+ *
+ * When SWIFT_TASKS_WITHOUT_ATOMICS is *not* defined this function uses an
+ * atomic operation.
+ *
+ * @param address The address to update.
+ * @param x The value to add to *address.
+ */
+__attribute__((always_inline)) INLINE static void accumulate_add_f(
+    volatile float *const address, const float x) {
+
+#ifdef SWIFT_TASKS_WITHOUT_ATOMICS
+  *address += x;
+#else
+  atomic_add_f(address, x);
+#endif
+}
+
+/**
+ * @brief Add x to the value stored at the location address (double version)
+ *
+ * When SWIFT_TASKS_WITHOUT_ATOMICS is *not* defined this function uses an
+ * atomic operation.
+ *
+ * @param address The address to update.
+ * @param x The value to add to *address.
+ */
+__attribute__((always_inline)) INLINE static void accumulate_add_d(
+    volatile double *const address, const double x) {
+
+#ifdef SWIFT_TASKS_WITHOUT_ATOMICS
+  *address += x;
+#else
+  atomic_add_d(address, x);
+#endif
+}
+
+#endif /* SWIFT_ACCUMULATE_H */
diff --git a/src/gravity/MultiSoftening/gravity.h b/src/gravity/MultiSoftening/gravity.h
index ce47dbe5ce13f03df8836a58972b544bc028fd6b..d7716cf712254054e3cf364a807f3cb1c12355f3 100644
--- a/src/gravity/MultiSoftening/gravity.h
+++ b/src/gravity/MultiSoftening/gravity.h
@@ -23,7 +23,7 @@
 #include <float.h>
 
 /* Local includes. */
-#include "atomic.h"
+#include "accumulate.h"
 #include "cosmology.h"
 #include "error.h"
 #include "gravity_properties.h"
@@ -78,11 +78,7 @@ __attribute__((always_inline)) INLINE static float gravity_get_softening(
 __attribute__((always_inline)) INLINE static void
 gravity_add_comoving_potential(struct gpart* restrict gp, float pot) {
 
-#ifdef SWIFT_TASKS_WITHOUT_ATOMICS
-  gp->potential += pot;
-#else
-  atomic_add_f(&gp->potential, pot);
-#endif
+  accumulate_add_f(&gp->potential, pot);
 }
 
 /**
diff --git a/src/gravity/Potential/gravity.h b/src/gravity/Potential/gravity.h
index e467e05e0d9d9adbeecbf2192d4d5c50b18a9b5b..f9a9502a528c161fbc82b3028f303b7f9cad49f8 100644
--- a/src/gravity/Potential/gravity.h
+++ b/src/gravity/Potential/gravity.h
@@ -23,7 +23,7 @@
 #include <float.h>
 
 /* Local includes. */
-#include "atomic.h"
+#include "accumulate.h"
 #include "cosmology.h"
 #include "gravity_properties.h"
 #include "kernel_gravity.h"
@@ -64,11 +64,7 @@ __attribute__((always_inline)) INLINE static float gravity_get_softening(
 __attribute__((always_inline)) INLINE static void
 gravity_add_comoving_potential(struct gpart* restrict gp, float pot) {
 
-#ifdef SWIFT_TASKS_WITHOUT_ATOMICS
-  gp->potential += pot;
-#else
-  atomic_add_f(&gp->potential, pot);
-#endif
+  accumulate_add_f(&gp->potential, pot);
 }
 
 /**
diff --git a/src/gravity_cache.h b/src/gravity_cache.h
index 23103f05098df06705e2a9a9ce1ce83d99f57bdc..ba3432d7e732b783547ddd2a597dbdf5723bf606 100644
--- a/src/gravity_cache.h
+++ b/src/gravity_cache.h
@@ -23,6 +23,7 @@
 #include "../config.h"
 
 /* Local headers */
+#include "accumulate.h"
 #include "align.h"
 #include "error.h"
 #include "gravity.h"
@@ -491,15 +492,9 @@ __attribute__((always_inline)) INLINE static void gravity_cache_write_back(
   /* Write stuff back to the particles */
   for (int i = 0; i < gcount; ++i) {
     if (active[i]) {
-#ifdef SWIFT_TASKS_WITHOUT_ATOMICS
-      gparts[i].a_grav[0] += a_x[i];
-      gparts[i].a_grav[1] += a_y[i];
-      gparts[i].a_grav[2] += a_z[i];
-#else
-      atomic_add_f(&gparts[i].a_grav[0], a_x[i]);
-      atomic_add_f(&gparts[i].a_grav[1], a_y[i]);
-      atomic_add_f(&gparts[i].a_grav[2], a_z[i]);
-#endif
+      accumulate_add_f(&gparts[i].a_grav[0], a_x[i]);
+      accumulate_add_f(&gparts[i].a_grav[1], a_y[i]);
+      accumulate_add_f(&gparts[i].a_grav[2], a_z[i]);
       gravity_add_comoving_potential(&gparts[i], pot[i]);
     }
   }
diff --git a/src/mesh_gravity.c b/src/mesh_gravity.c
index f36200718e952532da8c1c13453057deefb7eb51..cbfb9ef7a2970b222bf67a915bd5044bedaaaae3 100644
--- a/src/mesh_gravity.c
+++ b/src/mesh_gravity.c
@@ -28,6 +28,7 @@
 #include "mesh_gravity.h"
 
 /* Local includes. */
+#include "accumulate.h"
 #include "active.h"
 #include "debug.h"
 #include "engine.h"
@@ -334,15 +335,9 @@ void mesh_to_gparts_CIC(struct gpart* gp, const double* pot, const int N,
   /* ---- */
 
   /* Store things back */
-#ifdef SWIFT_TASKS_WITHOUT_ATOMICS
-  gp->a_grav[0] += fac * a[0];
-  gp->a_grav[1] += fac * a[1];
-  gp->a_grav[2] += fac * a[2];
-#else
-  atomic_add_f(&gp->a_grav[0], fac * a[0]);
-  atomic_add_f(&gp->a_grav[1], fac * a[1]);
-  atomic_add_f(&gp->a_grav[2], fac * a[2]);
-#endif
+  accumulate_add_f(&gp->a_grav[0], fac * a[0]);
+  accumulate_add_f(&gp->a_grav[1], fac * a[1]);
+  accumulate_add_f(&gp->a_grav[2], fac * a[2]);
   gravity_add_comoving_potential(gp, p);
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   gp->potential_PM = p;
diff --git a/src/multipole.h b/src/multipole.h
index e7dcfbeb1f7a5b4fc6c993d619ba257dca016bdd..fede08c09f424b181ed628a929020652e7a0695a 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -2589,21 +2589,14 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
 #endif
 
   /* Update the particle */
-#ifdef SWIFT_TASKS_WITHOUT_ATOMICS
-  gp->a_grav[0] += a_grav[0];
-  gp->a_grav[1] += a_grav[1];
-  gp->a_grav[2] += a_grav[2];
-#else
-  atomic_add_f(&gp->a_grav[0], a_grav[0]);
-  atomic_add_f(&gp->a_grav[1], a_grav[1]);
-  atomic_add_f(&gp->a_grav[2], a_grav[2]);
-#endif
-  gravity_add_comoving_potential(gp, pot);
+  accumulate_add_f(&gp->a_grav[0], a_grav[0]);
+  accumulate_add_f(&gp->a_grav[1], a_grav[1]);
+  accumulate_add_f(&gp->a_grav[2], a_grav[2]);
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
-  atomic_add_f(&gp->a_grav_m2l[0], a_grav[0]);
-  atomic_add_f(&gp->a_grav_m2l[1], a_grav[1]);
-  atomic_add_f(&gp->a_grav_m2l[2], a_grav[2]);
+  accumulate_add_f(&gp->a_grav_m2l[0], a_grav[0]);
+  accumulate_add_f(&gp->a_grav_m2l[1], a_grav[1]);
+  accumulate_add_f(&gp->a_grav_m2l[2], a_grav[2]);
 #endif
 }