diff --git a/src/gravity/MultiSoftening/gravity_iact.h b/src/gravity/MultiSoftening/gravity_iact.h
index 6fce3ddd512018e9ea4be21111c75904c77cb925..b0d0dcbb38ad095ff8ba531d40ce3077ae62685c 100644
--- a/src/gravity/MultiSoftening/gravity_iact.h
+++ b/src/gravity/MultiSoftening/gravity_iact.h
@@ -174,7 +174,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
   *f_y = m->M_000 * d.D_010;
   *f_z = m->M_000 * d.D_001;
 
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
 
   /* 1st order contributions */
 
@@ -279,7 +279,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
   *f_y = m->M_000 * d.D_010;
   *f_z = m->M_000 * d.D_001;
 
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
 
   /* 1st order contributions */
 
diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index 3dcffe1cc04c5e10d3b2353b7e21c532747c1475..71abcff1c9a4a325673fcca56e746dcf34597702 100644
--- a/src/gravity_derivatives.h
+++ b/src/gravity_derivatives.h
@@ -207,6 +207,10 @@ potential_derivatives_compute_M2L(const float r_x, const float r_y,
                                   const float r_s_inv,
                                   struct potential_derivatives_M2L *pot) {
 
+#ifdef SWIFT_DEBUG_CHECKS
+  if (r2 < eps * eps) error("Computing M2L derivatives below softening length");
+#endif
+
   float Dt_1;
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
   float Dt_3;
@@ -224,8 +228,8 @@ potential_derivatives_compute_M2L(const float r_x, const float r_y,
   float Dt_11;
 #endif
 
-  /* Un-softened un-truncated case (Newtonian potential) */
-  if (!periodic && r2 > eps * eps) {
+  /* Un-truncated case (Newtonian potential) */
+  if (!periodic) {
 
     Dt_1 = r_inv;
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
@@ -248,8 +252,8 @@ potential_derivatives_compute_M2L(const float r_x, const float r_y,
 #error "Missing implementation for order >5"
 #endif
 
-    /* Un-softened truncated case */
-  } else if (periodic && r2 > eps * eps) {
+    /* Truncated case */
+  } else {
 
     /* Get the derivatives of the truncated potential */
     const float r = r2 * r_inv;
@@ -291,38 +295,6 @@ potential_derivatives_compute_M2L(const float r_x, const float r_y,
 #endif
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 5
 #error "Missing implementation for order >5"
-#endif
-
-    /* Softened case */
-  } else {
-    const float r = r2 * r_inv;
-    const float u = r * eps_inv;
-    const float u_inv = r_inv * eps;
-
-    Dt_1 = eps_inv * D_soft_1(u, u_inv);
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
-    const float eps_inv2 = eps_inv * eps_inv;
-    const float eps_inv3 = eps_inv * eps_inv2;
-    Dt_3 = -eps_inv3 * D_soft_3(u, u_inv);
-#endif
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
-    const float eps_inv5 = eps_inv3 * eps_inv2;
-    Dt_5 = eps_inv5 * D_soft_5(u, u_inv);
-#endif
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
-    const float eps_inv7 = eps_inv5 * eps_inv2;
-    Dt_7 = -eps_inv7 * D_soft_7(u, u_inv);
-#endif
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
-    const float eps_inv9 = eps_inv7 * eps_inv2;
-    Dt_9 = eps_inv9 * D_soft_9(u, u_inv);
-#endif
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
-    const float eps_inv11 = eps_inv9 * eps_inv2;
-    Dt_11 = -eps_inv11 * D_soft_11(u, u_inv);
-#endif
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
-#error "Missing implementation for order >5"
 #endif
   }
 
@@ -463,6 +435,10 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
                                   const float r_s_inv,
                                   struct potential_derivatives_M2P *pot) {
 
+#ifdef SWIFT_DEBUG_CHECKS
+  if (r2 < eps * eps) error("Computing M2L derivatives below softening length");
+#endif
+
   float Dt_1;
   float Dt_3;
   float Dt_5;
@@ -471,8 +447,8 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
   float Dt_9;
 #endif
 
-  /* Un-softened un-truncated case (Newtonian potential) */
-  if (!periodic && r2 > eps * eps) {
+  /* Un-truncated case (Newtonian potential) */
+  if (!periodic) {
 
     const float r_inv2 = r_inv * r_inv;
 
@@ -484,8 +460,8 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
     Dt_9 = -7.f * Dt_7 * r_inv2; /* -105 / r^9 */
 #endif
 
-    /* Un-softened truncated case */
-  } else if (periodic && r2 > eps * eps) {
+    /* Truncated case */
+  } else if (periodic) {
 
     /* Get the derivatives of the truncated potential */
     const float r = r2 * r_inv;
@@ -512,30 +488,6 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
             45.f * r * r * d.chi_2 - 105.f * r * d.chi_1 + 105.f * d.chi_0) *
            r_inv9;
 #endif
-
-    /* Softened case */
-  } else {
-
-    const float r = r2 * r_inv;
-    const float u = r * eps_inv;
-    const float u_inv = r_inv * eps;
-    const float eps_inv2 = eps_inv * eps_inv;
-
-    Dt_1 = eps_inv * D_soft_1(u, u_inv);
-
-    const float eps_inv3 = eps_inv * eps_inv2;
-    Dt_3 = -eps_inv3 * D_soft_3(u, u_inv);
-
-    const float eps_inv5 = eps_inv3 * eps_inv2;
-    Dt_5 = eps_inv5 * D_soft_5(u, u_inv);
-
-    const float eps_inv7 = eps_inv5 * eps_inv2;
-    Dt_7 = -eps_inv7 * D_soft_7(u, u_inv);
-
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
-    const float eps_inv9 = eps_inv7 * eps_inv2;
-    Dt_9 = eps_inv9 * D_soft_9(u, u_inv);
-#endif
   }
 
   /* Compute some powers of r_x, r_y and r_z */
diff --git a/src/multipole.h b/src/multipole.h
index 956ec3e47a5c9b031c6cc396581ae21613d6f577..afe925f9b356fc68ea93a8a223332689e6e98a12 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -2039,7 +2039,7 @@ INLINE static void gravity_M2L_symmetric(
     const float rs_inv) {
 
   /* Recover some constants */
-  const float eps = m_a->max_softening;
+  const float eps = max(m_a->max_softening, m_b->max_softening);
   const float eps_inv = 1.f / eps;
 
   /* Compute distance vector */