diff --git a/src/engine.c b/src/engine.c
index b68b456ee24bcdd7a4b37c30bc49987fcd852691..508f1e28582556d87cdb973300660130ee17a24c 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1617,19 +1617,14 @@ void engine_exchange_top_multipoles(struct engine *e) {
   /* Each node (space) has constructed its own top-level multipoles.
    * We now need to make sure every other node has a copy of everything.
    *
-   * WARNING: Adult stuff ahead: don't do this at home!
-   *
-   * Since all nodes have their top-level multi-poles computed
-   * and all foreign ones set to 0 (all bytes), we can gather all the m-poles
-   * by doing a bit-wise OR reduction across all the nodes directly in
-   * place inside the multi-poles_top array.
-   * This only works if the foreign m-poles on every nodes are zeroed and no
-   * multi-pole is present on more than one node (two things guaranteed by the
-   * domain decomposition).
+   * We use our home-made reduction operation that simply performs a XOR
+   * operation on the multipoles. Since only local multipoles are non-zero and
+   * each multipole is only present once, the bit-by-bit XOR will
+   * create the desired result.
    */
-  int err = MPI_Allreduce(MPI_IN_PLACE, e->s->multipoles_top,
-                          e->s->nr_cells * sizeof(struct gravity_tensors),
-                          MPI_BYTE, MPI_BOR, MPI_COMM_WORLD);
+  int err = MPI_Allreduce(MPI_IN_PLACE, e->s->multipoles_top, e->s->nr_cells,
+                          multipole_mpi_type, multipole_mpi_reduce_op,
+                          MPI_COMM_WORLD);
   if (err != MPI_SUCCESS)
     mpi_error(err, "Failed to all-reduce the top-level multipoles.");
 
@@ -4688,6 +4683,7 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
 /* Construct types for MPI communications */
 #ifdef WITH_MPI
   part_create_mpi_types();
+  multipole_create_mpi_types();
   stats_create_mpi_type();
   proxy_create_mpi_type();
   task_create_mpi_comms();
diff --git a/src/multipole.c b/src/multipole.c
index bd5c6d6546fa0546108dcd53d7fe4060293c37a7..a77e6fce297802fb4118b7ac3d4c6a9bf4ecfd22 100644
--- a/src/multipole.c
+++ b/src/multipole.c
@@ -20,3 +20,70 @@
 
 /* Config parameters. */
 #include "../config.h"
+
+/* This object's header. */
+#include "multipole.h"
+
+/* MPI headers. */
+#ifdef WITH_MPI
+#include <mpi.h>
+#endif
+
+#ifdef WITH_MPI
+
+/* MPI data type for the multipole transfer and reduction */
+MPI_Datatype multipole_mpi_type;
+MPI_Op multipole_mpi_reduce_op;
+
+/**
+ * @brief Apply a bit-by-bit XOR operattion on #gravity_tensors (i.e. does
+ * a^=b).
+ *
+ * @param a The #gravity_tensors to add to.
+ * @param b The #gravity_tensors to add.
+ */
+void gravity_binary_xor(struct gravity_tensors *a,
+                        const struct gravity_tensors *b) {
+
+  char *aa = (char *)a;
+  const char *bb = (const char *)b;
+
+  for (size_t i = 0; i < sizeof(struct gravity_tensors); ++i) {
+    aa[i] ^= bb[i];
+  }
+}
+
+/**
+ * @brief MPI reduction function for the #gravity_tensors.
+ *
+ * @param invec Array of #gravity_tensors to read.
+ * @param inoutvec Array of #gravity_tensors to read and do the reduction into.
+ * @param len The length of the array.
+ * @param datatype The MPI type this function acts upon (unused).
+ */
+void gravity_tensors_mpi_reduce(void *invec, void *inoutvec, int *len,
+                                MPI_Datatype *datatype) {
+
+  for (int i = 0; i < *len; ++i) {
+    gravity_binary_xor(&((struct gravity_tensors *)inoutvec)[i],
+                       &((const struct gravity_tensors *)invec)[i]);
+  }
+}
+
+void multipole_create_mpi_types(void) {
+
+  /* Create the datatype for multipoles */
+  /* We just consider each structure to be a byte field disregarding their */
+  /* detailed content */
+  if (MPI_Type_contiguous(
+          sizeof(struct gravity_tensors) / sizeof(unsigned char), MPI_BYTE,
+          &multipole_mpi_type) != MPI_SUCCESS ||
+      MPI_Type_commit(&multipole_mpi_type) != MPI_SUCCESS) {
+    error("Failed to create MPI type for multipole.");
+  }
+
+  /* And the reduction operator */
+  MPI_Op_create(gravity_tensors_mpi_reduce, 1, &multipole_mpi_reduce_op);
+}
+
+#endif
diff --git a/src/multipole.h b/src/multipole.h
index dc22adfa311db0ada7cede4324564d862b43a496..8139dc0548bb94b108d6e32da4b19808998f48d3 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -207,6 +207,13 @@ struct gravity_tensors {
   };
 } SWIFT_STRUCT_ALIGN;
 
+#ifdef WITH_MPI
+/* MPI datatypes for transfers */
+extern MPI_Datatype multipole_mpi_type;
+extern MPI_Op multipole_mpi_reduce_op;
+void multipole_create_mpi_types(void);
+#endif
+
 /**
  * @brief Reset the data of a #multipole.
  *
diff --git a/src/part.c b/src/part.c
index f45223d7cb8cec3ccec690f046cf2e0181bb5fa9..3a626e652cf28f0376cadc1d9a40ab85b752e6c1 100644
--- a/src/part.c
+++ b/src/part.c
@@ -26,8 +26,10 @@
 #endif
 
 /* This object's header. */
-#include "error.h"
 #include "multipole.h"
+
+/* Local headers */
+#include "error.h"
 #include "part.h"
 
 /**
@@ -259,7 +261,6 @@ MPI_Datatype part_mpi_type;
 MPI_Datatype xpart_mpi_type;
 MPI_Datatype gpart_mpi_type;
 MPI_Datatype spart_mpi_type;
-MPI_Datatype multipole_mpi_type;
 
 /**
  * @brief Registers MPI particle types.
@@ -292,11 +293,5 @@ void part_create_mpi_types(void) {
       MPI_Type_commit(&spart_mpi_type) != MPI_SUCCESS) {
     error("Failed to create MPI type for sparts.");
   }
-  if (MPI_Type_contiguous(
-          sizeof(struct gravity_tensors) / sizeof(unsigned char), MPI_BYTE,
-          &multipole_mpi_type) != MPI_SUCCESS ||
-      MPI_Type_commit(&multipole_mpi_type) != MPI_SUCCESS) {
-    error("Failed to create MPI type for multipole.");
-  }
 }
 #endif
diff --git a/src/part.h b/src/part.h
index 36d4cc5ba2051aad1952cbb1a81102d62c076c17..625b928106cc18c219a250bf3ddac9bd36ddd3d2 100644
--- a/src/part.h
+++ b/src/part.h
@@ -108,7 +108,6 @@ extern MPI_Datatype part_mpi_type;
 extern MPI_Datatype xpart_mpi_type;
 extern MPI_Datatype gpart_mpi_type;
 extern MPI_Datatype spart_mpi_type;
-extern MPI_Datatype multipole_mpi_type;
 
 void part_create_mpi_types(void);
 #endif