diff --git a/src/approx_math.h b/src/approx_math.h
index e27418ce4d8bf0689219c842a99b8c1da8115d7f..f347bab44790d1e3120675bcbd6e7a457ca09821 100644
--- a/src/approx_math.h
+++ b/src/approx_math.h
@@ -21,7 +21,6 @@
 
 #include "inline.h"
 
-
 /**
  * @brief Approximate version of the complementay error function erfcf(x).
  *
@@ -31,7 +30,8 @@
  * Returns garbage for x < 0.
  * @param x The number to compute erfc for.
  */
-__attribute__((always_inline, const)) INLINE static float approx_erfcf(float x) {
+__attribute__((always_inline, const)) INLINE static float approx_erfcf(
+    float x) {
 
   /* 1 + 0.278393*x + 0.230389*x^2 + 0.000972*x^3 + 0.078108*x^4 */
   float arg = 0.078108f;
@@ -42,7 +42,7 @@ __attribute__((always_inline, const)) INLINE static float approx_erfcf(float x)
 
   /* 1 / arg^4 */
   const float arg2 = arg * arg;
-  const float arg4 = arg2 * arg2;  
+  const float arg4 = arg2 * arg2;
   return 1.f / arg4;
 }
 
diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index e6ab3729e979caaf3b6ec05dec92411636bff35f..3dcffe1cc04c5e10d3b2353b7e21c532747c1475 100644
--- a/src/gravity_derivatives.h
+++ b/src/gravity_derivatives.h
@@ -128,7 +128,7 @@ struct potential_derivatives_M2P {
 /**
  * @brief Converts the derivatives from a distance vector to its opposite.
  *
- * From a series of tensors D_xxx(\vec{r}), compute D_xxx(-\vec{r}).
+ * From a series of tensors D_xxx(r), compute D_xxx(-r).
  * This can be computed efficiently by flipping the sign of all the odd
  * derivative terms.
  *
diff --git a/src/multipole.h b/src/multipole.h
index bd00cd0840f614e4e57d4f35704af03c061acbec..2909eb5ca94266803650954e9b1cd0d97dec1f7e 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -293,8 +293,8 @@ INLINE static void gravity_field_tensors_init(struct grav_tensor *l,
  * @param la The gravity tensors to add to.
  * @param lb The gravity tensors to add.
  */
-INLINE static void gravity_field_tensors_add(struct grav_tensor *la,
-                                             const struct grav_tensor *lb) {
+INLINE static void gravity_field_tensors_add(
+    struct grav_tensor *restrict la, const struct grav_tensor *restrict lb) {
 #ifdef SWIFT_DEBUG_CHECKS
   if (lb->num_interacted == 0) error("Adding tensors that did not interact");
   la->num_interacted += lb->num_interacted;
@@ -502,8 +502,8 @@ INLINE static void gravity_multipole_print(const struct multipole *m) {
  * @param ma The multipole to add to.
  * @param mb The multipole to add.
  */
-INLINE static void gravity_multipole_add(struct multipole *ma,
-                                         const struct multipole *mb) {
+INLINE static void gravity_multipole_add(struct multipole *restrict ma,
+                                         const struct multipole *restrict mb) {
 
   /* Add 0th order term */
   ma->M_000 += mb->M_000;
@@ -1307,8 +1307,8 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
  * @param pos_a The position to which m_b will be shifted.
  * @param pos_b The current postion of the multipole to shift.
  */
-INLINE static void gravity_M2M(struct multipole *m_a,
-                               const struct multipole *m_b,
+INLINE static void gravity_M2M(struct multipole *restrict m_a,
+                               const struct multipole *restrict m_b,
                                const double pos_a[3], const double pos_b[3]) {
 
   /* Shift 0th order term */
@@ -1560,9 +1560,9 @@ INLINE static void gravity_M2M(struct multipole *m_a,
  * @param m_a The multipole creating the field.
  * @param pot The derivatives of the potential.
  */
-INLINE static void gravity_M2L_apply(struct grav_tensor *l_b,
-                                     const struct multipole *m_a,
-                                     const struct potential_derivatives_M2L *pot) {
+INLINE static void gravity_M2L_apply(
+    struct grav_tensor *restrict l_b, const struct multipole *restrict m_a,
+    const struct potential_derivatives_M2L *pot) {
 
 #ifdef SWIFT_DEBUG_CHECKS
   /* Count interactions */
@@ -1908,8 +1908,8 @@ INLINE static void gravity_M2L_apply(struct grav_tensor *l_b,
  */
 INLINE static void gravity_M2L_nonsym(
     struct grav_tensor *l_b, const struct multipole *m_a, const double pos_b[3],
-    const double pos_a[3], const struct gravity_props *props, int periodic,
-    const double dim[3], float rs_inv) {
+    const double pos_a[3], const struct gravity_props *props,
+    const int periodic, const double dim[3], const float rs_inv) {
 
   /* Recover some constants */
   const float eps = props->epsilon_cur;
@@ -1956,11 +1956,11 @@ INLINE static void gravity_M2L_nonsym(
  * @param rs_inv The inverse of the gravity mesh-smoothing scale.
  */
 INLINE static void gravity_M2L_symmetric(
-    struct grav_tensor *l_a, struct grav_tensor *l_b,
-    const struct multipole *m_a, const struct multipole *m_b,
+    struct grav_tensor *restrict l_a, struct grav_tensor *restrict l_b,
+    const struct multipole *restrict m_a, const struct multipole *restrict m_b,
     const double pos_b[3], const double pos_a[3],
-    const struct gravity_props *props, int periodic, const double dim[3],
-    float rs_inv) {
+    const struct gravity_props *props, const int periodic, const double dim[3],
+    const float rs_inv) {
 
   /* Recover some constants */
   const float eps = props->epsilon_cur;
@@ -2007,8 +2007,8 @@ INLINE static void gravity_M2L_symmetric(
  * @param pos_a The position to which m_b will be shifted.
  * @param pos_b The current postion of the multipole to shift.
  */
-INLINE static void gravity_L2L(struct grav_tensor *la,
-                               const struct grav_tensor *lb,
+INLINE static void gravity_L2L(struct grav_tensor *restrict la,
+                               const struct grav_tensor *restrict lb,
                                const double pos_a[3], const double pos_b[3]) {
 
   /* Initialise everything to zero */
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 29a00c7833cb9d0c39a492e90ea3dc3f816e1c87..0f36fe6a5ccafe263427f979c313dc8dc20d0c68 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -1198,9 +1198,9 @@ static INLINE void runner_dopair_grav_mm_symmetric(struct runner *r,
  * @param ci The #cell with field tensor to interact.
  * @param cj The #cell with the multipole.
  */
-static INLINE void runner_dopair_grav_mm_nonsym(struct runner *r,
-                                                struct cell *restrict ci,
-                                                struct cell *restrict cj) {
+static INLINE void runner_dopair_grav_mm_nonsym(
+    struct runner *r, struct cell *restrict ci,
+    const struct cell *restrict cj) {
 
   /* Some constants */
   const struct engine *e = r->e;
@@ -1247,13 +1247,16 @@ static INLINE void runner_dopair_grav_mm_nonsym(struct runner *r,
  * @param ci The first #cell.
  * @param cj The second #cell.
  */
-static INLINE void runner_dopair_grav_mm(struct runner *r, struct cell *ci,
-                                         struct cell *cj) {
+static INLINE void runner_dopair_grav_mm(struct runner *r,
+                                         struct cell *restrict ci,
+                                         struct cell *restrict cj) {
 
   const struct engine *e = r->e;
 
-  int do_i = cell_is_active_gravity_mm(ci, e) && (ci->nodeID == e->nodeID);
-  int do_j = cell_is_active_gravity_mm(cj, e) && (cj->nodeID == e->nodeID);
+  const int do_i =
+      cell_is_active_gravity_mm(ci, e) && (ci->nodeID == e->nodeID);
+  const int do_j =
+      cell_is_active_gravity_mm(cj, e) && (cj->nodeID == e->nodeID);
 
   if (do_i && do_j)
     runner_dopair_grav_mm_symmetric(r, ci, cj);