diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index f6eb888cf344671d7b5b6b19632b9a216c08b1ee..e808f534549867df320f08f94963f87424677a73 100644
--- a/src/gravity_derivatives.h
+++ b/src/gravity_derivatives.h
@@ -85,6 +85,20 @@ struct potential_derivatives {
 #endif
 };
 
+/**
+ * @brief Compute all the relevent derivatives of the softened and truncated
+ * gravitational potential.
+ *
+ * @param r_x x-component of distance vector
+ * @param r_y y-component of distance vector
+ * @param r_z z-component of distance vector
+ * @param r2 Square norm of distance vector
+ * @param r_inv Inverse norm of distance vector
+ * @param eps Softening length.
+ * @param eps2 Square of softening length.
+ * @param eps_inv Inverse of softening length.
+ * @param pot (return) The structure containing all the derivatives.
+ */
 __attribute__((always_inline)) INLINE static void compute_potential_derivatives(
     float r_x, float r_y, float r_z, float r2, float r_inv, float eps,
     float eps2, float eps_inv, struct potential_derivatives *pot) {
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index 8110b218eb064cdfc488afee9e0cdfb0c4f056f4..4782875ff98354f8b2a1d578524e4fa2a1807eb7 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -59,7 +59,7 @@ void gravity_props_init(struct gravity_props *p,
   /* Softening lengths */
   p->epsilon = 3. * parser_get_param_double(params, "Gravity:epsilon");
   p->epsilon2 = p->epsilon * p->epsilon;
-  p->epsilon_inv = 1. / p->epsilon;
+  p->epsilon_inv = 1.f / p->epsilon;
 }
 
 void gravity_props_print(const struct gravity_props *p) {
diff --git a/src/gravity_properties.h b/src/gravity_properties.h
index c663f628b48c1af8d701792c256d7d30db248137..886f29ff0834efba1486133c80d83aa86c317b31 100644
--- a/src/gravity_properties.h
+++ b/src/gravity_properties.h
@@ -58,13 +58,13 @@ struct gravity_props {
   double theta_crit_inv;
 
   /*! Softening length */
-  double epsilon;
+  float epsilon;
 
   /*! Square of softening length */
-  double epsilon2;
+  float epsilon2;
 
   /*! Inverse of softening length */
-  double epsilon_inv;
+  float epsilon_inv;
 };
 
 void gravity_props_print(const struct gravity_props *p);
diff --git a/src/multipole.h b/src/multipole.h
index 7eeccff73bed41503aab569b9436e77fcafbefb9..fc7866a28977e6ffa0ed3b822f5a0c014f8be6b4 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -1516,14 +1516,14 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
                                const double dim[3]) {
 
   /* Recover some constants */
-  const double eps = props->epsilon;
-  const double eps2 = props->epsilon2;
-  const double eps_inv = props->epsilon_inv;
+  const float eps = props->epsilon;
+  const float eps2 = props->epsilon2;
+  const float eps_inv = props->epsilon_inv;
 
   /* Compute distance vector */
-  double dx = pos_b[0] - pos_a[0];
-  double dy = pos_b[1] - pos_a[1];
-  double dz = pos_b[2] - pos_a[2];
+  float dx = (float)(pos_b[0] - pos_a[0]);
+  float dy = (float)(pos_b[1] - pos_a[1]);
+  float dz = (float)(pos_b[2] - pos_a[2]);
 
   /* Apply BC */
   if (periodic) {
@@ -1533,8 +1533,8 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
   }
 
   /* Compute distance */
-  const double r2 = dx * dx + dy * dy + dz * dz;
-  const double r_inv = 1. / sqrt(r2);
+  const float r2 = dx * dx + dy * dy + dz * dz;
+  const float r_inv = 1. / sqrtf(r2);
 
   /* Compute all derivatives */
   struct potential_derivatives pot;