diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index acf3ad79058ab82069ef5113843e61f2d196b426..1fa789032bb464b8887c6fc48d643bff35f6b6d0 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -148,10 +148,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
     float *restrict f_x, float *restrict f_y, float *restrict f_z,
     float *restrict 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 < 3
+#if SELF_GRAVITY_MULTIPOLE_ORDER < 2
 
   float f_ij, pot_ij;
   runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,
@@ -186,8 +186,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
   /* *f_z = m->M_001 * d.D_002 + m->M_010 * d.D_011 + m->M_100 * d.D_101 ; */
 
 #endif
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
 
   /* 2nd order contributions */
   *f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
@@ -198,8 +197,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
           m->M_101 * d.D_102 + m->M_110 * d.D_111 + m->M_200 * d.D_201;
 
 #endif
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
 
   /* 3rd order contributions */
   *f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
@@ -215,6 +213,25 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
           m->M_120 * d.D_121 + m->M_201 * d.D_202 + m->M_210 * d.D_211 +
           m->M_300 * d.D_301;
 
+#endif
+#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;
 #endif
 
   /* Take care of the the sign convention */
@@ -254,10 +271,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
     const struct multipole *m, float *restrict f_x, float *restrict f_y,
     float *restrict f_z, float *restrict 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 < 3
+#if SELF_GRAVITY_MULTIPOLE_ORDER < 2
 
   float f_ij, pot_ij;
   runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
@@ -292,8 +309,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
   /* *f_z = m->M_001 * d.D_002 + m->M_010 * d.D_011 + m->M_100 * d.D_101 ; */
 
 #endif
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
 
   /* 2nd order contributions */
   *f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
@@ -304,8 +320,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
           m->M_101 * d.D_102 + m->M_110 * d.D_111 + m->M_200 * d.D_201;
 
 #endif
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
 
   /* 3rd order contributions */
   *f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
@@ -321,7 +336,25 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
           m->M_120 * d.D_121 + m->M_201 * d.D_202 + m->M_210 * d.D_211 +
           m->M_300 * d.D_301;
 #endif
+#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;
+#endif
   /* Take care of the the sign convention */
   *f_x *= -1.f;
   *f_y *= -1.f;
diff --git a/src/gravity/MultiSoftening/gravity_iact.h b/src/gravity/MultiSoftening/gravity_iact.h
index 6432a6550eb3402faba385e57330fdbbadd88b17..f5020646e1bdc693005d01476735f9ad8485c843 100644
--- a/src/gravity/MultiSoftening/gravity_iact.h
+++ b/src/gravity/MultiSoftening/gravity_iact.h
@@ -150,10 +150,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
     float *restrict f_x, float *restrict f_y, float *restrict f_z,
     float *restrict 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 < 3
+#if SELF_GRAVITY_MULTIPOLE_ORDER < 2
 
   float f_ij, pot_ij;
   runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,
@@ -179,7 +179,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
   *f_z = m->M_000 * d.D_001;
   *pot = m->M_000 * d.D_000;
 
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
 
   /* 1st order contributions */
 
@@ -191,8 +191,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
   /* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
 
 #endif
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
 
   /* 2nd order contributions */
   *f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
@@ -205,8 +204,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
           m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
 
 #endif
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
 
   /* 3rd order contributions */
   *f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
@@ -227,7 +225,31 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
           m->M_300 * d.D_300;
 
 #endif
+#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;
+
+#endif
   /* Take care of the the sign convention */
   *f_x *= -1.f;
   *f_y *= -1.f;
@@ -263,10 +285,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
     const struct multipole *m, float *restrict f_x, float *restrict f_y,
     float *restrict f_z, float *restrict 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 < 3
+#if SELF_GRAVITY_MULTIPOLE_ORDER < 2
 
   float f_ij, pot_ij;
   runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
@@ -292,7 +314,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
   *f_z = m->M_000 * d.D_001;
   *pot = m->M_000 * d.D_000;
 
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
 
   /* 1st order contributions */
 
@@ -304,8 +326,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
   /* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
 
 #endif
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
 
   /* 2nd order contributions */
   *f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
@@ -318,8 +339,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
           m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
 
 #endif
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
 
   /* 3rd order contributions */
   *f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
@@ -340,7 +360,31 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
           m->M_300 * d.D_300;
 
 #endif
+#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;
+
+#endif
   /* Take care of the the sign convention */
   *f_x *= -1.f;
   *f_y *= -1.f;
diff --git a/src/gravity/Potential/gravity_iact.h b/src/gravity/Potential/gravity_iact.h
index c220982827ffb54e475dfe83ad53a0cce891a636..8d36ab666650a5ae6239657b247a0f2dad1f7756 100644
--- a/src/gravity/Potential/gravity_iact.h
+++ b/src/gravity/Potential/gravity_iact.h
@@ -150,10 +150,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
     float *restrict f_x, float *restrict f_y, float *restrict f_z,
     float *restrict 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 < 3
+#if SELF_GRAVITY_MULTIPOLE_ORDER < 2
 
   float f_ij, pot_ij;
   runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,
@@ -179,7 +179,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
   *f_z = m->M_000 * d.D_001;
   *pot = m->M_000 * d.D_000;
 
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
 
   /* 1st order contributions */
 
@@ -191,8 +191,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
   /* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
 
 #endif
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
 
   /* 2nd order contributions */
   *f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
@@ -205,8 +204,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
           m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
 
 #endif
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
 
   /* 3rd order contributions */
   *f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
@@ -227,7 +225,31 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
           m->M_300 * d.D_300;
 
 #endif
+#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;
+
+#endif
   /* Take care of the the sign convention */
   *f_x *= -1.f;
   *f_y *= -1.f;
@@ -263,10 +285,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
     const struct multipole *m, float *restrict f_x, float *restrict f_y,
     float *restrict f_z, float *restrict 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 < 3
+#if SELF_GRAVITY_MULTIPOLE_ORDER < 2
 
   float f_ij, pot_ij;
   runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
@@ -292,7 +314,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
   *f_z = m->M_000 * d.D_001;
   *pot = m->M_000 * d.D_000;
 
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
 
   /* 1st order contributions */
 
@@ -304,8 +326,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
   /* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
 
 #endif
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
 
   /* 2nd order contributions */
   *f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
@@ -318,8 +339,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
           m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
 
 #endif
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
 
   /* 3rd order contributions */
   *f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
@@ -340,7 +360,31 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
           m->M_300 * d.D_300;
 
 #endif
+#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;
+
+#endif
   /* Take care of the the sign convention */
   *f_x *= -1.f;
   *f_y *= -1.f;
diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index 8bd8340890dc7ced7a0573e2dc694c2e1650fc06..674779b54b8594fe80ed3595972b3f9950fb0f9d 100644
--- a/src/gravity_derivatives.h
+++ b/src/gravity_derivatives.h
@@ -102,9 +102,13 @@ struct potential_derivatives_M2P {
   /* 1st order terms */
   float D_100, D_010, D_001;
 
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+
   /* 2nd order terms */
   float D_200, D_020, D_002;
   float D_110, D_101, D_011;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
 
   /* 3rd order terms */
   float D_300, D_030, D_003;
@@ -112,8 +116,8 @@ struct potential_derivatives_M2P {
   float D_120, D_021;
   float D_102, D_012;
   float D_111;
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
 
   /* 4th order terms */
   float D_400, D_040, D_004;
@@ -123,6 +127,17 @@ struct potential_derivatives_M2P {
   float D_220, D_202, D_022;
   float D_211, D_121, D_112;
 #endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+
+  /* 5th order terms */
+  float D_005, D_014, D_023;
+  float D_032, D_041, D_050;
+  float D_104, D_113, D_122;
+  float D_131, D_140, D_203;
+  float D_212, D_221, D_230;
+  float D_302, D_311, D_320;
+  float D_401, D_410, D_500;
+#endif
 };
 
 /**
@@ -415,6 +430,9 @@ potential_derivatives_compute_M2L(const float r_x, const float r_y,
  * @brief Compute all the relevent derivatives of the softened and truncated
  * gravitational potential for the M2P kernel.
  *
+ * For M2P, we compute the derivatives to one order higher than
+ * SELF_GRAVITY_MULTIPOLE_ORDER, as these are needed for the accelerations.
+ *
  * @param r_x x-component of distance vector
  * @param r_y y-component of distance vector
  * @param r_z z-component of distance vector
@@ -437,71 +455,106 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
     error("Computing M2P derivatives below softening length");
 #endif
 
-  float Dt_1 = 0.f;
-  float Dt_3 = 0.f;
-  float Dt_5 = 0.f;
-  float Dt_7 = 0.f;
+  float Dt_1;
+  float Dt_3;
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+  float Dt_5;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+  float Dt_7;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+  float Dt_9;
+#endif
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 3
-  float Dt_9 = 0.f;
+  float Dt_11;
 #endif
 
   /* Un-truncated case (Newtonian potential) */
   if (!periodic) {
 
     const float r_inv2 = r_inv * r_inv;
-
     Dt_1 = r_inv;
     Dt_3 = -1.f * Dt_1 * r_inv2; /* -1 / r^3 */
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
     Dt_5 = -3.f * Dt_3 * r_inv2; /* 3 / r^5 */
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
     Dt_7 = -5.f * Dt_5 * r_inv2; /* -15 / r^7 */
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+    Dt_9 = -7.f * Dt_7 * r_inv2; /* 105 / r^9 */
+#endif
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 3
-    Dt_9 = -7.f * Dt_7 * r_inv2; /* -105 / r^9 */
+    Dt_11 = -9.f * Dt_9 * r_inv2; /* -945 / r^11 */
 #endif
 
     /* Truncated case */
-  } else if (periodic) {
+  } else {
 
     /* Get the derivatives of the truncated potential */
     const float r = r2 * r_inv;
-    struct chi_derivatives d;
-    kernel_long_grav_derivatives(r, r_s_inv, &d);
+    struct chi_derivatives derivs;
+    kernel_long_grav_derivatives(r, r_s_inv, &derivs);
 
     const float r_inv2 = r_inv * r_inv;
-    Dt_1 = d.chi_0 * r_inv;
-
     const float r_inv3 = r_inv2 * r_inv;
-    Dt_3 = (r * d.chi_1 - d.chi_0) * r_inv3;
-
+    Dt_1 = derivs.chi_0 * r_inv;
+    Dt_3 = (r * derivs.chi_1 - derivs.chi_0) * r_inv3;
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
     const float r_inv5 = r_inv2 * r_inv3;
-    Dt_5 = (r * r * d.chi_2 - 3.f * r * d.chi_1 + 3.f * d.chi_0) * r_inv5;
-
+    Dt_5 =
+        (r * r * derivs.chi_2 - 3.f * r * derivs.chi_1 + 3.f * derivs.chi_0) *
+        r_inv5;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
     const float r_inv7 = r_inv2 * r_inv5;
-    Dt_7 = (r * r * r * d.chi_3 - 6.f * r * r * d.chi_2 + 15.f * r * d.chi_1 -
-            15.f * d.chi_0) *
+    Dt_7 = (r * r * r * derivs.chi_3 - 6.f * r * r * derivs.chi_2 +
+            15.f * r * derivs.chi_1 - 15.f * derivs.chi_0) *
            r_inv7;
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
     const float r_inv9 = r_inv2 * r_inv7;
-    Dt_9 = (r * r * r * r * d.chi_4 - 10.f * r * r * r * d.chi_3 +
-            45.f * r * r * d.chi_2 - 105.f * r * d.chi_1 + 105.f * d.chi_0) *
+    Dt_9 = (r * r * r * r * derivs.chi_4 - 10.f * r * r * r * derivs.chi_3 +
+            45.f * r * r * derivs.chi_2 - 105.f * r * derivs.chi_1 +
+            105.f * derivs.chi_0) *
            r_inv9;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+    const float r_inv11 = r_inv2 * r_inv9;
+    Dt_11 = (r * r * r * r * r * derivs.chi_5 -
+             15.f * r * r * r * r * derivs.chi_4 +
+             105.f * r * r * r * derivs.chi_3 - 420.f * r * r * derivs.chi_2 +
+             945.f * r * derivs.chi_1 - 945.f * derivs.chi_0) *
+            r_inv11;
 #endif
   }
 
-  /* Compute some powers of r_x, r_y and r_z */
+/* Alright, let's get the full terms */
+
+/* Compute some powers of r_x, r_y and r_z */
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
   const float r_x2 = r_x * r_x;
   const float r_y2 = r_y * r_y;
   const float r_z2 = r_z * r_z;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
   const float r_x3 = r_x2 * r_x;
   const float r_y3 = r_y2 * r_y;
   const float r_z3 = r_z2 * r_z;
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
   const float r_x4 = r_x3 * r_x;
   const float r_y4 = r_y3 * r_y;
   const float r_z4 = r_z3 * r_z;
 #endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+  const float r_x5 = r_x4 * r_x;
+  const float r_y5 = r_y4 * r_y;
+  const float r_z5 = r_z4 * r_z;
+#endif
 
-  /* 0th order derivative */
+  /* Get the 0th order term */
   pot->D_000 = Dt_1;
 
   /* 1st order derivatives */
@@ -509,6 +562,7 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
   pot->D_010 = r_y * Dt_3;
   pot->D_001 = r_z * Dt_3;
 
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
   /* 2nd order derivatives */
   pot->D_200 = r_x2 * Dt_5 + Dt_3;
   pot->D_020 = r_y2 * Dt_5 + Dt_3;
@@ -516,7 +570,8 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
   pot->D_110 = r_x * r_y * Dt_5;
   pot->D_101 = r_x * r_z * Dt_5;
   pot->D_011 = r_y * r_z * Dt_5;
-
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
   /* 3rd order derivatives */
   pot->D_300 = r_x3 * Dt_7 + 3.f * r_x * Dt_5;
   pot->D_030 = r_y3 * Dt_7 + 3.f * r_y * Dt_5;
@@ -528,8 +583,8 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
   pot->D_102 = r_z2 * r_x * Dt_7 + r_x * Dt_5;
   pot->D_012 = r_z2 * r_y * Dt_7 + r_y * Dt_5;
   pot->D_111 = r_x * r_y * r_z * Dt_7;
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
   /* 4th order derivatives */
   pot->D_400 = r_x4 * Dt_9 + 6.f * r_x2 * Dt_7 + 3.f * Dt_5;
   pot->D_040 = r_y4 * Dt_9 + 6.f * r_y2 * Dt_7 + 3.f * Dt_5;
@@ -547,6 +602,39 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
   pot->D_121 = r_y2 * r_x * r_z * Dt_9 + r_x * r_z * Dt_7;
   pot->D_112 = r_z2 * r_x * r_y * Dt_9 + r_x * r_y * Dt_7;
 #endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+  /* 5th order derivatives */
+  pot->D_500 = r_x5 * Dt_11 + 10.f * r_x3 * Dt_9 + 15.f * r_x * Dt_7;
+  pot->D_050 = r_y5 * Dt_11 + 10.f * r_y3 * Dt_9 + 15.f * r_y * Dt_7;
+  pot->D_005 = r_z5 * Dt_11 + 10.f * r_z3 * Dt_9 + 15.f * r_z * Dt_7;
+  pot->D_410 = r_x4 * r_y * Dt_11 + 6.f * r_x2 * r_y * Dt_9 + 3.f * r_y * Dt_7;
+  pot->D_401 = r_x4 * r_z * Dt_11 + 6.f * r_x2 * r_z * Dt_9 + 3.f * r_z * Dt_7;
+  pot->D_140 = r_y4 * r_x * Dt_11 + 6.f * r_y2 * r_x * Dt_9 + 3.f * r_x * Dt_7;
+  pot->D_041 = r_y4 * r_z * Dt_11 + 6.f * r_y2 * r_z * Dt_9 + 3.f * r_z * Dt_7;
+  pot->D_104 = r_z4 * r_x * Dt_11 + 6.f * r_z2 * r_x * Dt_9 + 3.f * r_x * Dt_7;
+  pot->D_014 = r_z4 * r_y * Dt_11 + 6.f * r_z2 * r_y * Dt_9 + 3.f * r_y * Dt_7;
+  pot->D_320 = r_x3 * r_y2 * Dt_11 + r_x3 * Dt_9 + 3.f * r_x * r_y2 * Dt_9 +
+               3.f * r_x * Dt_7;
+  pot->D_302 = r_x3 * r_z2 * Dt_11 + r_x3 * Dt_9 + 3.f * r_x * r_z2 * Dt_9 +
+               3.f * r_x * Dt_7;
+  pot->D_230 = r_y3 * r_x2 * Dt_11 + r_y3 * Dt_9 + 3.f * r_y * r_x2 * Dt_9 +
+               3.f * r_y * Dt_7;
+  pot->D_032 = r_y3 * r_z2 * Dt_11 + r_y3 * Dt_9 + 3.f * r_y * r_z2 * Dt_9 +
+               3.f * r_y * Dt_7;
+  pot->D_203 = r_z3 * r_x2 * Dt_11 + r_z3 * Dt_9 + 3.f * r_z * r_x2 * Dt_9 +
+               3.f * r_z * Dt_7;
+  pot->D_023 = r_z3 * r_y2 * Dt_11 + r_z3 * Dt_9 + 3.f * r_z * r_y2 * Dt_9 +
+               3.f * r_z * Dt_7;
+  pot->D_311 = r_x3 * r_y * r_z * Dt_11 + 3.f * r_x * r_y * r_z * Dt_9;
+  pot->D_131 = r_y3 * r_x * r_z * Dt_11 + 3.f * r_x * r_y * r_z * Dt_9;
+  pot->D_113 = r_z3 * r_x * r_y * Dt_11 + 3.f * r_x * r_y * r_z * Dt_9;
+  pot->D_122 = r_x * r_y2 * r_z2 * Dt_11 + r_x * r_y2 * Dt_9 +
+               r_x * r_z2 * Dt_9 + r_x * Dt_7;
+  pot->D_212 = r_y * r_x2 * r_z2 * Dt_11 + r_y * r_x2 * Dt_9 +
+               r_y * r_z2 * Dt_9 + r_y * Dt_7;
+  pot->D_221 = r_z * r_x2 * r_y2 * Dt_11 + r_z * r_x2 * Dt_9 +
+               r_z * r_y2 * Dt_9 + r_z * Dt_7;
+#endif
 }
 
 #endif /* SWIFT_GRAVITY_DERIVATIVE_H */
diff --git a/tests/testGravityDerivatives.c b/tests/testGravityDerivatives.c
index 032460d9f519c455e485e0dddfc5f8742138d60a..9dba2dc929715b8da8763cc3dcd87a53e72d6d87 100644
--- a/tests/testGravityDerivatives.c
+++ b/tests/testGravityDerivatives.c
@@ -935,13 +935,14 @@ int main(int argc, char* argv[]) {
   message("Seed = %d", seed);
   srand(seed);
 
+  /* Start by testing M2L */
   for (int i = 0; i < 100; ++i) {
 
     const double dx = 100. * ((double)rand() / (RAND_MAX));
     const double dy = 100. * ((double)rand() / (RAND_MAX));
     const double dz = 100. * ((double)rand() / (RAND_MAX));
 
-    message("Testing gravity for r=(%e %e %e)", dx, dy, dz);
+    message("Testing M2L gravity for r=(%e %e %e)", dx, dy, dz);
 
     const double r_s = 100. * ((double)rand() / (RAND_MAX));
     const double r_s_inv = 1. / r_s;
@@ -958,6 +959,7 @@ int main(int argc, char* argv[]) {
 
     /* Compute all derivatives */
     struct potential_derivatives_M2L pot;
+    bzero(&pot, sizeof(struct potential_derivatives_M2L));
     potential_derivatives_compute_M2L(dx, dy, dz, r2, r_inv, eps, periodic,
                                       r_s_inv, &pot);
 
@@ -967,89 +969,211 @@ int main(int argc, char* argv[]) {
     /* Now check everything... */
 
     /* 0th order terms */
-    test(pot.D_000, D_000(dx, dy, dz, r_inv), tol, min, "D_000");
+    test(pot.D_000, D_000(dx, dy, dz, r_inv), tol, min, "M2L D_000");
 
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
 
     /* 1st order terms */
-    test(pot.D_100, D_100(dx, dy, dz, r_inv), tol, min, "D_100");
-    test(pot.D_010, D_010(dx, dy, dz, r_inv), tol, min, "D_010");
-    test(pot.D_001, D_001(dx, dy, dz, r_inv), tol, min, "D_001");
+    test(pot.D_100, D_100(dx, dy, dz, r_inv), tol, min, "M2L D_100");
+    test(pot.D_010, D_010(dx, dy, dz, r_inv), tol, min, "M2L D_010");
+    test(pot.D_001, D_001(dx, dy, dz, r_inv), tol, min, "M2L D_001");
 #endif
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 1
 
     /* 2nd order terms */
-    test(pot.D_200, D_200(dx, dy, dz, r_inv), tol, min, "D_200");
-    test(pot.D_020, D_020(dx, dy, dz, r_inv), tol, min, "D_020");
-    test(pot.D_002, D_002(dx, dy, dz, r_inv), tol, min, "D_002");
-    test(pot.D_110, D_110(dx, dy, dz, r_inv), tol, min, "D_110");
-    test(pot.D_101, D_101(dx, dy, dz, r_inv), tol, min, "D_101");
-    test(pot.D_011, D_011(dx, dy, dz, r_inv), tol, min, "D_011");
+    test(pot.D_200, D_200(dx, dy, dz, r_inv), tol, min, "M2L D_200");
+    test(pot.D_020, D_020(dx, dy, dz, r_inv), tol, min, "M2L D_020");
+    test(pot.D_002, D_002(dx, dy, dz, r_inv), tol, min, "M2L D_002");
+    test(pot.D_110, D_110(dx, dy, dz, r_inv), tol, min, "M2L D_110");
+    test(pot.D_101, D_101(dx, dy, dz, r_inv), tol, min, "M2L D_101");
+    test(pot.D_011, D_011(dx, dy, dz, r_inv), tol, min, "M2L D_011");
 #endif
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 2
 
     tol *= 2.5;
 
     /* 3rd order terms */
-    test(pot.D_300, D_300(dx, dy, dz, r_inv), tol, min, "D_300");
-    test(pot.D_030, D_030(dx, dy, dz, r_inv), tol, min, "D_030");
-    test(pot.D_003, D_003(dx, dy, dz, r_inv), tol, min, "D_003");
-    test(pot.D_210, D_210(dx, dy, dz, r_inv), tol, min, "D_210");
-    test(pot.D_201, D_201(dx, dy, dz, r_inv), tol, min, "D_201");
-    test(pot.D_120, D_120(dx, dy, dz, r_inv), tol, min, "D_120");
-    test(pot.D_021, D_021(dx, dy, dz, r_inv), tol, min, "D_021");
-    test(pot.D_102, D_102(dx, dy, dz, r_inv), tol, min, "D_102");
-    test(pot.D_012, D_012(dx, dy, dz, r_inv), tol, min, "D_012");
-    test(pot.D_111, D_111(dx, dy, dz, r_inv), tol, min, "D_111");
+    test(pot.D_300, D_300(dx, dy, dz, r_inv), tol, min, "M2L D_300");
+    test(pot.D_030, D_030(dx, dy, dz, r_inv), tol, min, "M2L D_030");
+    test(pot.D_003, D_003(dx, dy, dz, r_inv), tol, min, "M2L D_003");
+    test(pot.D_210, D_210(dx, dy, dz, r_inv), tol, min, "M2L D_210");
+    test(pot.D_201, D_201(dx, dy, dz, r_inv), tol, min, "M2L D_201");
+    test(pot.D_120, D_120(dx, dy, dz, r_inv), tol, min, "M2L D_120");
+    test(pot.D_021, D_021(dx, dy, dz, r_inv), tol, min, "M2L D_021");
+    test(pot.D_102, D_102(dx, dy, dz, r_inv), tol, min, "M2L D_102");
+    test(pot.D_012, D_012(dx, dy, dz, r_inv), tol, min, "M2L D_012");
+    test(pot.D_111, D_111(dx, dy, dz, r_inv), tol, min, "M2L D_111");
 #endif
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 3
 
     /* 4th order terms */
-    test(pot.D_400, D_400(dx, dy, dz, r_inv), tol, min, "D_400");
-    test(pot.D_040, D_040(dx, dy, dz, r_inv), tol, min, "D_040");
-    test(pot.D_004, D_004(dx, dy, dz, r_inv), tol, min, "D_004");
-    test(pot.D_310, D_310(dx, dy, dz, r_inv), tol, min, "D_310");
-    test(pot.D_301, D_301(dx, dy, dz, r_inv), tol, min, "D_301");
-    test(pot.D_130, D_130(dx, dy, dz, r_inv), tol, min, "D_130");
-    test(pot.D_031, D_031(dx, dy, dz, r_inv), tol, min, "D_031");
-    test(pot.D_103, D_103(dx, dy, dz, r_inv), tol, min, "D_103");
-    test(pot.D_013, D_013(dx, dy, dz, r_inv), tol, min, "D_013");
-    test(pot.D_220, D_220(dx, dy, dz, r_inv), tol, min, "D_220");
-    test(pot.D_202, D_202(dx, dy, dz, r_inv), tol, min, "D_202");
-    test(pot.D_022, D_022(dx, dy, dz, r_inv), tol, min, "D_022");
-    test(pot.D_211, D_211(dx, dy, dz, r_inv), tol, min, "D_211");
-    test(pot.D_121, D_121(dx, dy, dz, r_inv), tol, min, "D_121");
-    test(pot.D_112, D_112(dx, dy, dz, r_inv), tol, min, "D_112");
+    test(pot.D_400, D_400(dx, dy, dz, r_inv), tol, min, "M2L D_400");
+    test(pot.D_040, D_040(dx, dy, dz, r_inv), tol, min, "M2L D_040");
+    test(pot.D_004, D_004(dx, dy, dz, r_inv), tol, min, "M2L D_004");
+    test(pot.D_310, D_310(dx, dy, dz, r_inv), tol, min, "M2L D_310");
+    test(pot.D_301, D_301(dx, dy, dz, r_inv), tol, min, "M2L D_301");
+    test(pot.D_130, D_130(dx, dy, dz, r_inv), tol, min, "M2L D_130");
+    test(pot.D_031, D_031(dx, dy, dz, r_inv), tol, min, "M2L D_031");
+    test(pot.D_103, D_103(dx, dy, dz, r_inv), tol, min, "M2L D_103");
+    test(pot.D_013, D_013(dx, dy, dz, r_inv), tol, min, "M2L D_013");
+    test(pot.D_220, D_220(dx, dy, dz, r_inv), tol, min, "M2L D_220");
+    test(pot.D_202, D_202(dx, dy, dz, r_inv), tol, min, "M2L D_202");
+    test(pot.D_022, D_022(dx, dy, dz, r_inv), tol, min, "M2L D_022");
+    test(pot.D_211, D_211(dx, dy, dz, r_inv), tol, min, "M2L D_211");
+    test(pot.D_121, D_121(dx, dy, dz, r_inv), tol, min, "M2L D_121");
+    test(pot.D_112, D_112(dx, dy, dz, r_inv), tol, min, "M2L D_112");
 #endif
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 4
 
     tol *= 2.5;
 
     /* 5th order terms */
-    test(pot.D_500, D_500(dx, dy, dz, r_inv), tol, min, "D_500");
-    test(pot.D_050, D_050(dx, dy, dz, r_inv), tol, min, "D_050");
-    test(pot.D_005, D_005(dx, dy, dz, r_inv), tol, min, "D_005");
-    test(pot.D_410, D_410(dx, dy, dz, r_inv), tol, min, "D_410");
-    test(pot.D_401, D_401(dx, dy, dz, r_inv), tol, min, "D_401");
-    test(pot.D_140, D_140(dx, dy, dz, r_inv), tol, min, "D_140");
-    test(pot.D_041, D_041(dx, dy, dz, r_inv), tol, min, "D_041");
-    test(pot.D_104, D_104(dx, dy, dz, r_inv), tol, min, "D_104");
-    test(pot.D_014, D_014(dx, dy, dz, r_inv), tol, min, "D_014");
-    test(pot.D_320, D_320(dx, dy, dz, r_inv), tol, min, "D_320");
-    test(pot.D_302, D_302(dx, dy, dz, r_inv), tol, min, "D_302");
-    test(pot.D_230, D_230(dx, dy, dz, r_inv), tol, min, "D_230");
-    test(pot.D_032, D_032(dx, dy, dz, r_inv), tol, min, "D_032");
-    test(pot.D_203, D_203(dx, dy, dz, r_inv), tol, min, "D_203");
-    test(pot.D_023, D_023(dx, dy, dz, r_inv), tol, min, "D_023");
-    test(pot.D_311, D_311(dx, dy, dz, r_inv), tol, min, "D_311");
-    test(pot.D_131, D_131(dx, dy, dz, r_inv), tol, min, "D_131");
-    test(pot.D_113, D_113(dx, dy, dz, r_inv), tol, min, "D_113");
-    test(pot.D_122, D_122(dx, dy, dz, r_inv), tol, min, "D_122");
-    test(pot.D_212, D_212(dx, dy, dz, r_inv), tol, min, "D_212");
-    test(pot.D_221, D_221(dx, dy, dz, r_inv), tol, min, "D_221");
+    test(pot.D_500, D_500(dx, dy, dz, r_inv), tol, min, "M2L D_500");
+    test(pot.D_050, D_050(dx, dy, dz, r_inv), tol, min, "M2L D_050");
+    test(pot.D_005, D_005(dx, dy, dz, r_inv), tol, min, "M2L D_005");
+    test(pot.D_410, D_410(dx, dy, dz, r_inv), tol, min, "M2L D_410");
+    test(pot.D_401, D_401(dx, dy, dz, r_inv), tol, min, "M2L D_401");
+    test(pot.D_140, D_140(dx, dy, dz, r_inv), tol, min, "M2L D_140");
+    test(pot.D_041, D_041(dx, dy, dz, r_inv), tol, min, "M2L D_041");
+    test(pot.D_104, D_104(dx, dy, dz, r_inv), tol, min, "M2L D_104");
+    test(pot.D_014, D_014(dx, dy, dz, r_inv), tol, min, "M2L D_014");
+    test(pot.D_320, D_320(dx, dy, dz, r_inv), tol, min, "M2L D_320");
+    test(pot.D_302, D_302(dx, dy, dz, r_inv), tol, min, "M2L D_302");
+    test(pot.D_230, D_230(dx, dy, dz, r_inv), tol, min, "M2L D_230");
+    test(pot.D_032, D_032(dx, dy, dz, r_inv), tol, min, "M2L D_032");
+    test(pot.D_203, D_203(dx, dy, dz, r_inv), tol, min, "M2L D_203");
+    test(pot.D_023, D_023(dx, dy, dz, r_inv), tol, min, "M2L D_023");
+    test(pot.D_311, D_311(dx, dy, dz, r_inv), tol, min, "M2L D_311");
+    test(pot.D_131, D_131(dx, dy, dz, r_inv), tol, min, "M2L D_131");
+    test(pot.D_113, D_113(dx, dy, dz, r_inv), tol, min, "M2L D_113");
+    test(pot.D_122, D_122(dx, dy, dz, r_inv), tol, min, "M2L D_122");
+    test(pot.D_212, D_212(dx, dy, dz, r_inv), tol, min, "M2L D_212");
+    test(pot.D_221, D_221(dx, dy, dz, r_inv), tol, min, "M2L D_221");
 
 #endif
     message("All good!");
   }
+
+  /* And now the M2P terms */
+  for (int i = 0; i < 100; ++i) {
+
+    const double dx = 100. * ((double)rand() / (RAND_MAX));
+    const double dy = 100. * ((double)rand() / (RAND_MAX));
+    const double dz = 100. * ((double)rand() / (RAND_MAX));
+
+    message("Testing M2P gravity for r=(%e %e %e)", dx, dy, dz);
+
+    const double r_s = 100. * ((double)rand() / (RAND_MAX));
+    const double r_s_inv = 1. / r_s;
+
+    const int periodic = 0;
+
+    message("Mesh scale r_s=%e periodic=%d", r_s, periodic);
+
+    /* Compute distance */
+    const double r2 = dx * dx + dy * dy + dz * dz;
+    const double r_inv = 1. / sqrt(r2);
+    const double r = r2 * r_inv;
+    const double eps = r / 10.;
+
+    /* Compute all derivatives */
+    struct potential_derivatives_M2P pot;
+    bzero(&pot, sizeof(struct potential_derivatives_M2P));
+    potential_derivatives_compute_M2P(dx, dy, dz, r2, r_inv, eps, periodic,
+                                      r_s_inv, &pot);
+
+    /* Minimal value we care about */
+    const double min = 1e-9;
+
+    /* Now check everything...
+     *
+     * Note that the M2P derivatives are computed to order
+     * SELF_GRAVITY_MULTIPOLE_ORDER + 1 by the function above. */
+
+    /* 0th order terms */
+    test(pot.D_000, D_000(dx, dy, dz, r_inv), tol, min, "M2P D_000");
+
+    /* 1st order terms */
+    test(pot.D_100, D_100(dx, dy, dz, r_inv), tol, min, "M2P D_100");
+    test(pot.D_010, D_010(dx, dy, dz, r_inv), tol, min, "M2P D_010");
+    test(pot.D_001, D_001(dx, dy, dz, r_inv), tol, min, "M2P D_001");
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+
+    /* 2nd order terms */
+    test(pot.D_200, D_200(dx, dy, dz, r_inv), tol, min, "M2P D_200");
+    test(pot.D_020, D_020(dx, dy, dz, r_inv), tol, min, "M2P D_020");
+    test(pot.D_002, D_002(dx, dy, dz, r_inv), tol, min, "M2P D_002");
+    test(pot.D_110, D_110(dx, dy, dz, r_inv), tol, min, "M2P D_110");
+    test(pot.D_101, D_101(dx, dy, dz, r_inv), tol, min, "M2P D_101");
+    test(pot.D_011, D_011(dx, dy, dz, r_inv), tol, min, "M2P D_011");
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+
+    tol *= 2.5;
+
+    /* 3rd order terms */
+    test(pot.D_300, D_300(dx, dy, dz, r_inv), tol, min, "M2P D_300");
+    test(pot.D_030, D_030(dx, dy, dz, r_inv), tol, min, "M2P D_030");
+    test(pot.D_003, D_003(dx, dy, dz, r_inv), tol, min, "M2P D_003");
+    test(pot.D_210, D_210(dx, dy, dz, r_inv), tol, min, "M2P D_210");
+    test(pot.D_201, D_201(dx, dy, dz, r_inv), tol, min, "M2P D_201");
+    test(pot.D_120, D_120(dx, dy, dz, r_inv), tol, min, "M2P D_120");
+    test(pot.D_021, D_021(dx, dy, dz, r_inv), tol, min, "M2P D_021");
+    test(pot.D_102, D_102(dx, dy, dz, r_inv), tol, min, "M2P D_102");
+    test(pot.D_012, D_012(dx, dy, dz, r_inv), tol, min, "M2P D_012");
+    test(pot.D_111, D_111(dx, dy, dz, r_inv), tol, min, "M2P D_111");
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+
+    /* 4th order terms */
+    test(pot.D_400, D_400(dx, dy, dz, r_inv), tol, min, "M2P D_400");
+    test(pot.D_040, D_040(dx, dy, dz, r_inv), tol, min, "M2P D_040");
+    test(pot.D_004, D_004(dx, dy, dz, r_inv), tol, min, "M2P D_004");
+    test(pot.D_310, D_310(dx, dy, dz, r_inv), tol, min, "M2P D_310");
+    test(pot.D_301, D_301(dx, dy, dz, r_inv), tol, min, "M2P D_301");
+    test(pot.D_130, D_130(dx, dy, dz, r_inv), tol, min, "M2P D_130");
+    test(pot.D_031, D_031(dx, dy, dz, r_inv), tol, min, "M2P D_031");
+    test(pot.D_103, D_103(dx, dy, dz, r_inv), tol, min, "M2P D_103");
+    test(pot.D_013, D_013(dx, dy, dz, r_inv), tol, min, "M2P D_013");
+    test(pot.D_220, D_220(dx, dy, dz, r_inv), tol, min, "M2P D_220");
+    test(pot.D_202, D_202(dx, dy, dz, r_inv), tol, min, "M2P D_202");
+    test(pot.D_022, D_022(dx, dy, dz, r_inv), tol, min, "M2P D_022");
+    test(pot.D_211, D_211(dx, dy, dz, r_inv), tol, min, "M2P D_211");
+    test(pot.D_121, D_121(dx, dy, dz, r_inv), tol, min, "M2P D_121");
+    test(pot.D_112, D_112(dx, dy, dz, r_inv), tol, min, "M2P D_112");
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+
+    tol *= 2.5;
+
+    /* 5th order terms */
+    test(pot.D_500, D_500(dx, dy, dz, r_inv), tol, min, "M2P D_500");
+    test(pot.D_050, D_050(dx, dy, dz, r_inv), tol, min, "M2P D_050");
+    test(pot.D_005, D_005(dx, dy, dz, r_inv), tol, min, "M2P D_005");
+    test(pot.D_410, D_410(dx, dy, dz, r_inv), tol, min, "M2P D_410");
+    test(pot.D_401, D_401(dx, dy, dz, r_inv), tol, min, "M2P D_401");
+    test(pot.D_140, D_140(dx, dy, dz, r_inv), tol, min, "M2P D_140");
+    test(pot.D_041, D_041(dx, dy, dz, r_inv), tol, min, "M2P D_041");
+    test(pot.D_104, D_104(dx, dy, dz, r_inv), tol, min, "M2P D_104");
+    test(pot.D_014, D_014(dx, dy, dz, r_inv), tol, min, "M2P D_014");
+    test(pot.D_320, D_320(dx, dy, dz, r_inv), tol, min, "M2P D_320");
+    test(pot.D_302, D_302(dx, dy, dz, r_inv), tol, min, "M2P D_302");
+    test(pot.D_230, D_230(dx, dy, dz, r_inv), tol, min, "M2P D_230");
+    test(pot.D_032, D_032(dx, dy, dz, r_inv), tol, min, "M2P D_032");
+    test(pot.D_203, D_203(dx, dy, dz, r_inv), tol, min, "M2P D_203");
+    test(pot.D_023, D_023(dx, dy, dz, r_inv), tol, min, "M2P D_023");
+    test(pot.D_311, D_311(dx, dy, dz, r_inv), tol, min, "M2P D_311");
+    test(pot.D_131, D_131(dx, dy, dz, r_inv), tol, min, "M2P D_131");
+    test(pot.D_113, D_113(dx, dy, dz, r_inv), tol, min, "M2P D_113");
+    test(pot.D_122, D_122(dx, dy, dz, r_inv), tol, min, "M2P D_122");
+    test(pot.D_212, D_212(dx, dy, dz, r_inv), tol, min, "M2P D_212");
+    test(pot.D_221, D_221(dx, dy, dz, r_inv), tol, min, "M2P D_221");
+
+#endif
+    message("All good!");
+  }
+
+  /* All happy */
   return 0;
 }