diff --git a/src/gravity/MultiSoftening/gravity_iact.h b/src/gravity/MultiSoftening/gravity_iact.h
index 0d1a949ea808187b216c2f45a60010fce5951da0..c3f416fdee5af046197ad660712cd415090b4ae5 100644
--- a/src/gravity/MultiSoftening/gravity_iact.h
+++ b/src/gravity/MultiSoftening/gravity_iact.h
@@ -149,7 +149,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
     const float h, const float h_inv, const struct multipole *m, float *f_x,
     float *f_y, float *f_z, float *pot) {
 
-/* In the case where the order is < 3, then there is only a monopole term left.
+/* In the case where the order is < 2, then there is only a monopole term left.
  * We can default to the normal P-P interaction with the mass of the multipole
  * and its CoM as the "particle" property */
 #if SELF_GRAVITY_MULTIPOLE_ORDER < 2
@@ -227,26 +227,26 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 3
 
   /* 4th order contributions */
-  *f_x = m->M_004 * d.D_104 + m->M_013 * d.D_113 + m->M_022 * d.D_122 +
-         m->M_031 * d.D_131 + m->M_040 * d.D_140 + m->M_103 * d.D_203 +
-         m->M_112 * d.D_212 + m->M_121 * d.D_221 + m->M_130 * d.D_230 +
-         m->M_202 * d.D_302 + m->M_211 * d.D_311 + m->M_220 * d.D_320 +
-         m->M_301 * d.D_401 + m->M_310 * d.D_410 + m->M_400 * d.D_500;
-  *f_y = m->M_004 * d.D_014 + m->M_013 * d.D_023 + m->M_022 * d.D_032 +
-         m->M_031 * d.D_041 + m->M_040 * d.D_050 + m->M_103 * d.D_113 +
-         m->M_112 * d.D_122 + m->M_121 * d.D_131 + m->M_130 * d.D_140 +
-         m->M_202 * d.D_212 + m->M_211 * d.D_221 + m->M_220 * d.D_230 +
-         m->M_301 * d.D_311 + m->M_310 * d.D_320 + m->M_400 * d.D_410;
-  *f_z = m->M_004 * d.D_005 + m->M_013 * d.D_014 + m->M_022 * d.D_023 +
-         m->M_031 * d.D_032 + m->M_040 * d.D_041 + m->M_103 * d.D_104 +
-         m->M_112 * d.D_113 + m->M_121 * d.D_122 + m->M_130 * d.D_131 +
-         m->M_202 * d.D_203 + m->M_211 * d.D_212 + m->M_220 * d.D_221 +
-         m->M_301 * d.D_302 + m->M_310 * d.D_311 + m->M_400 * d.D_401;
-  *pot = m->M_004 * d.D_004 + m->M_013 * d.D_013 + m->M_022 * d.D_022 +
-         m->M_031 * d.D_031 + m->M_040 * d.D_040 + m->M_103 * d.D_103 +
-         m->M_112 * d.D_112 + m->M_121 * d.D_121 + m->M_130 * d.D_130 +
-         m->M_202 * d.D_202 + m->M_211 * d.D_211 + m->M_220 * d.D_220 +
-         m->M_301 * d.D_301 + m->M_310 * d.D_310 + m->M_400 * d.D_400;
+  *f_x += m->M_004 * d.D_104 + m->M_013 * d.D_113 + m->M_022 * d.D_122 +
+          m->M_031 * d.D_131 + m->M_040 * d.D_140 + m->M_103 * d.D_203 +
+          m->M_112 * d.D_212 + m->M_121 * d.D_221 + m->M_130 * d.D_230 +
+          m->M_202 * d.D_302 + m->M_211 * d.D_311 + m->M_220 * d.D_320 +
+          m->M_301 * d.D_401 + m->M_310 * d.D_410 + m->M_400 * d.D_500;
+  *f_y += m->M_004 * d.D_014 + m->M_013 * d.D_023 + m->M_022 * d.D_032 +
+          m->M_031 * d.D_041 + m->M_040 * d.D_050 + m->M_103 * d.D_113 +
+          m->M_112 * d.D_122 + m->M_121 * d.D_131 + m->M_130 * d.D_140 +
+          m->M_202 * d.D_212 + m->M_211 * d.D_221 + m->M_220 * d.D_230 +
+          m->M_301 * d.D_311 + m->M_310 * d.D_320 + m->M_400 * d.D_410;
+  *f_z += m->M_004 * d.D_005 + m->M_013 * d.D_014 + m->M_022 * d.D_023 +
+          m->M_031 * d.D_032 + m->M_040 * d.D_041 + m->M_103 * d.D_104 +
+          m->M_112 * d.D_113 + m->M_121 * d.D_122 + m->M_130 * d.D_131 +
+          m->M_202 * d.D_203 + m->M_211 * d.D_212 + m->M_220 * d.D_221 +
+          m->M_301 * d.D_302 + m->M_310 * d.D_311 + m->M_400 * d.D_401;
+  *pot += m->M_004 * d.D_004 + m->M_013 * d.D_013 + m->M_022 * d.D_022 +
+          m->M_031 * d.D_031 + m->M_040 * d.D_040 + m->M_103 * d.D_103 +
+          m->M_112 * d.D_112 + m->M_121 * d.D_121 + m->M_130 * d.D_130 +
+          m->M_202 * d.D_202 + m->M_211 * d.D_211 + m->M_220 * d.D_220 +
+          m->M_301 * d.D_301 + m->M_310 * d.D_310 + m->M_400 * d.D_400;
 
 #endif
   /* Take care of the the sign convention */
@@ -283,7 +283,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
     const float h, const float h_inv, const float r_s_inv,
     const struct multipole *m, float *f_x, float *f_y, float *f_z, float *pot) {
 
-/* In the case where the order is < 3, then there is only a monopole term left.
+/* In the case where the order is < 2, then there is only a monopole term left.
  * We can default to the normal P-P interaction with the mass of the multipole
  * and its CoM as the "particle" property */
 #if SELF_GRAVITY_MULTIPOLE_ORDER < 2
@@ -361,26 +361,26 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 3
 
   /* 4th order contributions */
-  *f_x = m->M_004 * d.D_104 + m->M_013 * d.D_113 + m->M_022 * d.D_122 +
-         m->M_031 * d.D_131 + m->M_040 * d.D_140 + m->M_103 * d.D_203 +
-         m->M_112 * d.D_212 + m->M_121 * d.D_221 + m->M_130 * d.D_230 +
-         m->M_202 * d.D_302 + m->M_211 * d.D_311 + m->M_220 * d.D_320 +
-         m->M_301 * d.D_401 + m->M_310 * d.D_410 + m->M_400 * d.D_500;
-  *f_y = m->M_004 * d.D_014 + m->M_013 * d.D_023 + m->M_022 * d.D_032 +
-         m->M_031 * d.D_041 + m->M_040 * d.D_050 + m->M_103 * d.D_113 +
-         m->M_112 * d.D_122 + m->M_121 * d.D_131 + m->M_130 * d.D_140 +
-         m->M_202 * d.D_212 + m->M_211 * d.D_221 + m->M_220 * d.D_230 +
-         m->M_301 * d.D_311 + m->M_310 * d.D_320 + m->M_400 * d.D_410;
-  *f_z = m->M_004 * d.D_005 + m->M_013 * d.D_014 + m->M_022 * d.D_023 +
-         m->M_031 * d.D_032 + m->M_040 * d.D_041 + m->M_103 * d.D_104 +
-         m->M_112 * d.D_113 + m->M_121 * d.D_122 + m->M_130 * d.D_131 +
-         m->M_202 * d.D_203 + m->M_211 * d.D_212 + m->M_220 * d.D_221 +
-         m->M_301 * d.D_302 + m->M_310 * d.D_311 + m->M_400 * d.D_401;
-  *pot = m->M_004 * d.D_004 + m->M_013 * d.D_013 + m->M_022 * d.D_022 +
-         m->M_031 * d.D_031 + m->M_040 * d.D_040 + m->M_103 * d.D_103 +
-         m->M_112 * d.D_112 + m->M_121 * d.D_121 + m->M_130 * d.D_130 +
-         m->M_202 * d.D_202 + m->M_211 * d.D_211 + m->M_220 * d.D_220 +
-         m->M_301 * d.D_301 + m->M_310 * d.D_310 + m->M_400 * d.D_400;
+  *f_x += m->M_004 * d.D_104 + m->M_013 * d.D_113 + m->M_022 * d.D_122 +
+          m->M_031 * d.D_131 + m->M_040 * d.D_140 + m->M_103 * d.D_203 +
+          m->M_112 * d.D_212 + m->M_121 * d.D_221 + m->M_130 * d.D_230 +
+          m->M_202 * d.D_302 + m->M_211 * d.D_311 + m->M_220 * d.D_320 +
+          m->M_301 * d.D_401 + m->M_310 * d.D_410 + m->M_400 * d.D_500;
+  *f_y += m->M_004 * d.D_014 + m->M_013 * d.D_023 + m->M_022 * d.D_032 +
+          m->M_031 * d.D_041 + m->M_040 * d.D_050 + m->M_103 * d.D_113 +
+          m->M_112 * d.D_122 + m->M_121 * d.D_131 + m->M_130 * d.D_140 +
+          m->M_202 * d.D_212 + m->M_211 * d.D_221 + m->M_220 * d.D_230 +
+          m->M_301 * d.D_311 + m->M_310 * d.D_320 + m->M_400 * d.D_410;
+  *f_z += m->M_004 * d.D_005 + m->M_013 * d.D_014 + m->M_022 * d.D_023 +
+          m->M_031 * d.D_032 + m->M_040 * d.D_041 + m->M_103 * d.D_104 +
+          m->M_112 * d.D_113 + m->M_121 * d.D_122 + m->M_130 * d.D_131 +
+          m->M_202 * d.D_203 + m->M_211 * d.D_212 + m->M_220 * d.D_221 +
+          m->M_301 * d.D_302 + m->M_310 * d.D_311 + m->M_400 * d.D_401;
+  *pot += m->M_004 * d.D_004 + m->M_013 * d.D_013 + m->M_022 * d.D_022 +
+          m->M_031 * d.D_031 + m->M_040 * d.D_040 + m->M_103 * d.D_103 +
+          m->M_112 * d.D_112 + m->M_121 * d.D_121 + m->M_130 * d.D_130 +
+          m->M_202 * d.D_202 + m->M_211 * d.D_211 + m->M_220 * d.D_220 +
+          m->M_301 * d.D_301 + m->M_310 * d.D_310 + m->M_400 * d.D_400;
 
 #endif
   /* Take care of the the sign convention */
diff --git a/src/gravity/Potential/gravity_iact.h b/src/gravity/Potential/gravity_iact.h
index 37db5730569ec891b6e8dfce6e62424ed89438d8..5a3a8a11f91de837694c231bcc04871b6b1c5cd4 100644
--- a/src/gravity/Potential/gravity_iact.h
+++ b/src/gravity/Potential/gravity_iact.h
@@ -227,26 +227,26 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 3
 
   /* 4th order contributions */
-  *f_x = m->M_004 * d.D_104 + m->M_013 * d.D_113 + m->M_022 * d.D_122 +
-         m->M_031 * d.D_131 + m->M_040 * d.D_140 + m->M_103 * d.D_203 +
-         m->M_112 * d.D_212 + m->M_121 * d.D_221 + m->M_130 * d.D_230 +
-         m->M_202 * d.D_302 + m->M_211 * d.D_311 + m->M_220 * d.D_320 +
-         m->M_301 * d.D_401 + m->M_310 * d.D_410 + m->M_400 * d.D_500;
-  *f_y = m->M_004 * d.D_014 + m->M_013 * d.D_023 + m->M_022 * d.D_032 +
-         m->M_031 * d.D_041 + m->M_040 * d.D_050 + m->M_103 * d.D_113 +
-         m->M_112 * d.D_122 + m->M_121 * d.D_131 + m->M_130 * d.D_140 +
-         m->M_202 * d.D_212 + m->M_211 * d.D_221 + m->M_220 * d.D_230 +
-         m->M_301 * d.D_311 + m->M_310 * d.D_320 + m->M_400 * d.D_410;
-  *f_z = m->M_004 * d.D_005 + m->M_013 * d.D_014 + m->M_022 * d.D_023 +
-         m->M_031 * d.D_032 + m->M_040 * d.D_041 + m->M_103 * d.D_104 +
-         m->M_112 * d.D_113 + m->M_121 * d.D_122 + m->M_130 * d.D_131 +
-         m->M_202 * d.D_203 + m->M_211 * d.D_212 + m->M_220 * d.D_221 +
-         m->M_301 * d.D_302 + m->M_310 * d.D_311 + m->M_400 * d.D_401;
-  *pot = m->M_004 * d.D_004 + m->M_013 * d.D_013 + m->M_022 * d.D_022 +
-         m->M_031 * d.D_031 + m->M_040 * d.D_040 + m->M_103 * d.D_103 +
-         m->M_112 * d.D_112 + m->M_121 * d.D_121 + m->M_130 * d.D_130 +
-         m->M_202 * d.D_202 + m->M_211 * d.D_211 + m->M_220 * d.D_220 +
-         m->M_301 * d.D_301 + m->M_310 * d.D_310 + m->M_400 * d.D_400;
+  *f_x += m->M_004 * d.D_104 + m->M_013 * d.D_113 + m->M_022 * d.D_122 +
+          m->M_031 * d.D_131 + m->M_040 * d.D_140 + m->M_103 * d.D_203 +
+          m->M_112 * d.D_212 + m->M_121 * d.D_221 + m->M_130 * d.D_230 +
+          m->M_202 * d.D_302 + m->M_211 * d.D_311 + m->M_220 * d.D_320 +
+          m->M_301 * d.D_401 + m->M_310 * d.D_410 + m->M_400 * d.D_500;
+  *f_y += m->M_004 * d.D_014 + m->M_013 * d.D_023 + m->M_022 * d.D_032 +
+          m->M_031 * d.D_041 + m->M_040 * d.D_050 + m->M_103 * d.D_113 +
+          m->M_112 * d.D_122 + m->M_121 * d.D_131 + m->M_130 * d.D_140 +
+          m->M_202 * d.D_212 + m->M_211 * d.D_221 + m->M_220 * d.D_230 +
+          m->M_301 * d.D_311 + m->M_310 * d.D_320 + m->M_400 * d.D_410;
+  *f_z += m->M_004 * d.D_005 + m->M_013 * d.D_014 + m->M_022 * d.D_023 +
+          m->M_031 * d.D_032 + m->M_040 * d.D_041 + m->M_103 * d.D_104 +
+          m->M_112 * d.D_113 + m->M_121 * d.D_122 + m->M_130 * d.D_131 +
+          m->M_202 * d.D_203 + m->M_211 * d.D_212 + m->M_220 * d.D_221 +
+          m->M_301 * d.D_302 + m->M_310 * d.D_311 + m->M_400 * d.D_401;
+  *pot += m->M_004 * d.D_004 + m->M_013 * d.D_013 + m->M_022 * d.D_022 +
+          m->M_031 * d.D_031 + m->M_040 * d.D_040 + m->M_103 * d.D_103 +
+          m->M_112 * d.D_112 + m->M_121 * d.D_121 + m->M_130 * d.D_130 +
+          m->M_202 * d.D_202 + m->M_211 * d.D_211 + m->M_220 * d.D_220 +
+          m->M_301 * d.D_301 + m->M_310 * d.D_310 + m->M_400 * d.D_400;
 
 #endif
   /* Take care of the the sign convention */
@@ -361,26 +361,26 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 3
 
   /* 4th order contributions */
-  *f_x = m->M_004 * d.D_104 + m->M_013 * d.D_113 + m->M_022 * d.D_122 +
-         m->M_031 * d.D_131 + m->M_040 * d.D_140 + m->M_103 * d.D_203 +
-         m->M_112 * d.D_212 + m->M_121 * d.D_221 + m->M_130 * d.D_230 +
-         m->M_202 * d.D_302 + m->M_211 * d.D_311 + m->M_220 * d.D_320 +
-         m->M_301 * d.D_401 + m->M_310 * d.D_410 + m->M_400 * d.D_500;
-  *f_y = m->M_004 * d.D_014 + m->M_013 * d.D_023 + m->M_022 * d.D_032 +
-         m->M_031 * d.D_041 + m->M_040 * d.D_050 + m->M_103 * d.D_113 +
-         m->M_112 * d.D_122 + m->M_121 * d.D_131 + m->M_130 * d.D_140 +
-         m->M_202 * d.D_212 + m->M_211 * d.D_221 + m->M_220 * d.D_230 +
-         m->M_301 * d.D_311 + m->M_310 * d.D_320 + m->M_400 * d.D_410;
-  *f_z = m->M_004 * d.D_005 + m->M_013 * d.D_014 + m->M_022 * d.D_023 +
-         m->M_031 * d.D_032 + m->M_040 * d.D_041 + m->M_103 * d.D_104 +
-         m->M_112 * d.D_113 + m->M_121 * d.D_122 + m->M_130 * d.D_131 +
-         m->M_202 * d.D_203 + m->M_211 * d.D_212 + m->M_220 * d.D_221 +
-         m->M_301 * d.D_302 + m->M_310 * d.D_311 + m->M_400 * d.D_401;
-  *pot = m->M_004 * d.D_004 + m->M_013 * d.D_013 + m->M_022 * d.D_022 +
-         m->M_031 * d.D_031 + m->M_040 * d.D_040 + m->M_103 * d.D_103 +
-         m->M_112 * d.D_112 + m->M_121 * d.D_121 + m->M_130 * d.D_130 +
-         m->M_202 * d.D_202 + m->M_211 * d.D_211 + m->M_220 * d.D_220 +
-         m->M_301 * d.D_301 + m->M_310 * d.D_310 + m->M_400 * d.D_400;
+  *f_x += m->M_004 * d.D_104 + m->M_013 * d.D_113 + m->M_022 * d.D_122 +
+          m->M_031 * d.D_131 + m->M_040 * d.D_140 + m->M_103 * d.D_203 +
+          m->M_112 * d.D_212 + m->M_121 * d.D_221 + m->M_130 * d.D_230 +
+          m->M_202 * d.D_302 + m->M_211 * d.D_311 + m->M_220 * d.D_320 +
+          m->M_301 * d.D_401 + m->M_310 * d.D_410 + m->M_400 * d.D_500;
+  *f_y += m->M_004 * d.D_014 + m->M_013 * d.D_023 + m->M_022 * d.D_032 +
+          m->M_031 * d.D_041 + m->M_040 * d.D_050 + m->M_103 * d.D_113 +
+          m->M_112 * d.D_122 + m->M_121 * d.D_131 + m->M_130 * d.D_140 +
+          m->M_202 * d.D_212 + m->M_211 * d.D_221 + m->M_220 * d.D_230 +
+          m->M_301 * d.D_311 + m->M_310 * d.D_320 + m->M_400 * d.D_410;
+  *f_z += m->M_004 * d.D_005 + m->M_013 * d.D_014 + m->M_022 * d.D_023 +
+          m->M_031 * d.D_032 + m->M_040 * d.D_041 + m->M_103 * d.D_104 +
+          m->M_112 * d.D_113 + m->M_121 * d.D_122 + m->M_130 * d.D_131 +
+          m->M_202 * d.D_203 + m->M_211 * d.D_212 + m->M_220 * d.D_221 +
+          m->M_301 * d.D_302 + m->M_310 * d.D_311 + m->M_400 * d.D_401;
+  *pot += m->M_004 * d.D_004 + m->M_013 * d.D_013 + m->M_022 * d.D_022 +
+          m->M_031 * d.D_031 + m->M_040 * d.D_040 + m->M_103 * d.D_103 +
+          m->M_112 * d.D_112 + m->M_121 * d.D_121 + m->M_130 * d.D_130 +
+          m->M_202 * d.D_202 + m->M_211 * d.D_211 + m->M_220 * d.D_220 +
+          m->M_301 * d.D_301 + m->M_310 * d.D_310 + m->M_400 * d.D_400;
 
 #endif
   /* Take care of the the sign convention */