Added files for new gravity implementation that ignores the calculation of the potential.

configure.ac

@@ -917,7 +917,26 @@ case "$with_subgrid" in
    ;;
 esac
 
+# Gravity scheme.
+AC_ARG_WITH([gravity],
+   [AS_HELP_STRING([--with-gravity=<scheme>],
+      [Gravity scheme to use @<:@default, with-potential, default: default@:>@]
+   )],
+   [with_gravity="$withval"],
+   [with_gravity="default"]
+)
 
+case "$with_gravity" in
+   with-potential)
+      AC_DEFINE([POTENTIAL_GRAVITY], [1], [Gravity scheme with potential calculation])
+   ;;
+   default)
+      AC_DEFINE([DEFAULT_GRAVITY], [1], [Default gravity scheme])
+   ;;
+   *)
+      AC_MSG_ERROR([Unknown gravity scheme: $with_gravity])
+   ;;
+esac
 
 # Hydro scheme.
 AC_ARG_WITH([hydro],
@@ -1351,12 +1370,14 @@ AC_MSG_RESULT([
    Equation of state  : $with_eos
    Adiabatic index    : $with_gamma
    Riemann solver     : $with_riemann
-   Cooling function   : $with_cooling
-   Chemistry          : $with_chemistry
 
-   External potential  : $with_potential
+   Gravity scheme      : $with_gravity
    Multipole order     : $with_multipole_order
    No gravity below ID : $no_gravity_below_id
+   External potential  : $with_potential
+
+   Cooling function   : $with_cooling
+   Chemistry          : $with_chemistry
 
    Individual timers     : $enable_timers
    Task debugging        : $enable_task_debugging

src/debug.c

@@ -62,7 +62,14 @@
 #error "Invalid choice of SPH variant"
 #endif
 
+/* Import the right gravity definition */
+#if defined(DEFAULT_GRAVITY)
 #include "./gravity/Default/gravity_debug.h"
+#elif defined(POTENTIAL_GRAVITY)
+#include "./gravity/Potential/gravity_debug.h"
+#else
+#error "Invalid choice of gravity variant"
+#endif
 
 /**
  * @brief Looks for the particle with the given id and prints its information to

src/gravity.h

@@ -27,10 +27,19 @@
 #include "inline.h"
 #include "part.h"
 
-/* So far only one model here */
-/* Straight-forward import */
+/* Import the right functions */
+#if defined(DEFAULT_GRAVITY)
 #include "./gravity/Default/gravity.h"
 #include "./gravity/Default/gravity_iact.h"
+#define GRAVITY_IMPLEMENTATION "Default (no potential)"
+#elif defined(POTENTIAL_GRAVITY)
+#include "./gravity/Potential/gravity.h"
+#include "./gravity/Potential/gravity_iact.h"
+#define GRAVITY_IMPLEMENTATION "With potential calculation"
+#else
+#error "Invalid choice of gravity variant"
+#endif
+
 
 struct engine;
 struct space;

src/gravity/Potential/gravity.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *               2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_DEFAULT_GRAVITY_H
+#define SWIFT_DEFAULT_GRAVITY_H
+
+#include <float.h>
+
+#include "cosmology.h"
+#include "gravity_properties.h"
+#include "kernel_gravity.h"
+#include "minmax.h"
+
+/**
+ * @brief Returns the mass of a particle
+ *
+ * @param gp The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float gravity_get_mass(
+    const struct gpart* restrict gp) {
+
+  return gp->mass;
+}
+
+/**
+ * @brief Returns the softening of a particle
+ *
+ * @param gp The particle of interest
+ * @param grav_props The global gravity properties.
+ */
+__attribute__((always_inline)) INLINE static float gravity_get_softening(
+    const struct gpart* gp, const struct gravity_props* restrict grav_props) {
+
+  return grav_props->epsilon_cur;
+}
+
+/**
+ * @brief Returns the comoving potential of a particle
+ *
+ * @param gp The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float
+gravity_get_comoving_potential(const struct gpart* restrict gp) {
+
+  return gp->potential;
+}
+
+/**
+ * @brief Returns the physical potential of a particle
+ *
+ * @param gp The particle of interest.
+ * @param cosmo The cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+gravity_get_physical_potential(const struct gpart* restrict gp,
+                               const struct cosmology* cosmo) {
+
+  return gp->potential * cosmo->a_inv;
+}
+
+/**
+ * @brief Computes the gravity time-step of a given particle due to self-gravity
+ *
+ * We use Gadget-2's type 0 time-step criterion.
+ *
+ * @param gp Pointer to the g-particle data.
+ * @param a_hydro The accelerations coming from the hydro scheme (can be 0).
+ * @param grav_props Constants used in the gravity scheme.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static float
+gravity_compute_timestep_self(const struct gpart* const gp,
+                              const float a_hydro[3],
+                              const struct gravity_props* restrict grav_props,
+                              const struct cosmology* cosmo) {
+
+  /* Get physical acceleration (gravity contribution) */
+  float a_phys_x = gp->a_grav[0] * cosmo->a_factor_grav_accel;
+  float a_phys_y = gp->a_grav[1] * cosmo->a_factor_grav_accel;
+  float a_phys_z = gp->a_grav[2] * cosmo->a_factor_grav_accel;
+
+  /* Get physical acceleration (hydro contribution) */
+  a_phys_x += a_hydro[0] * cosmo->a_factor_hydro_accel;
+  a_phys_y += a_hydro[1] * cosmo->a_factor_hydro_accel;
+  a_phys_z += a_hydro[2] * cosmo->a_factor_hydro_accel;
+
+  const float ac2 =
+      a_phys_x * a_phys_x + a_phys_y * a_phys_y + a_phys_z * a_phys_z;
+
+  const float ac_inv = (ac2 > 0.f) ? 1.f / sqrtf(ac2) : FLT_MAX;
+
+  const float epsilon = gravity_get_softening(gp, grav_props);
+
+  const float dt = sqrtf(2. * kernel_gravity_softening_plummer_equivalent_inv *
+                         cosmo->a * grav_props->eta * epsilon * ac_inv);
+
+  return dt;
+}
+
+/**
+ * @brief Prepares a g-particle for the gravity calculation
+ *
+ * Zeroes all the relevant arrays in preparation for the sums taking place in
+ * the variaous tasks
+ *
+ * @param gp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void gravity_init_gpart(
+    struct gpart* gp) {
+
+  /* Zero the acceleration */
+  gp->a_grav[0] = 0.f;
+  gp->a_grav[1] = 0.f;
+  gp->a_grav[2] = 0.f;
+  gp->potential = 0.f;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  gp->potential_PM = 0.f;
+  gp->a_grav_PM[0] = 0.f;
+  gp->a_grav_PM[1] = 0.f;
+  gp->a_grav_PM[2] = 0.f;
+#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
+  gp->num_interacted = 0;
+#endif
+}
+
+/**
+ * @brief Finishes the gravity calculation.
+ *
+ * Multiplies the forces and accelerations by the appropiate constants
+ *
+ * @param gp The particle to act upon
+ * @param const_G Newton's constant in internal units
+ */
+__attribute__((always_inline)) INLINE static void gravity_end_force(
+    struct gpart* gp, float const_G) {
+
+  /* Let's get physical... */
+  gp->a_grav[0] *= const_G;
+  gp->a_grav[1] *= const_G;
+  gp->a_grav[2] *= const_G;
+  gp->potential *= const_G;
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+  gp->potential_PM *= const_G;
+  gp->a_grav_PM[0] *= const_G;
+  gp->a_grav_PM[1] *= const_G;
+  gp->a_grav_PM[2] *= const_G;
+#endif
+}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * @param gp The particle to act upon
+ * @param dt The time-step for this kick
+ */
+__attribute__((always_inline)) INLINE static void gravity_kick_extra(
+    struct gpart* gp, float dt) {}
+
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param gp The particle.
+ */
+__attribute__((always_inline)) INLINE static void
+gravity_reset_predicted_values(struct gpart* gp) {}
+
+/**
+ * @brief Initialises the g-particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param gp The particle to act upon
+ * @param grav_props The global properties of the gravity calculation.
+ */
+__attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
+    struct gpart* gp, const struct gravity_props* grav_props) {
+
+  gp->time_bin = 0;
+
+  gravity_init_gpart(gp);
+}
+
+#endif /* SWIFT_DEFAULT_GRAVITY_H */

src/gravity/Potential/gravity_debug.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_DEFAULT_GRAVITY_DEBUG_H
+#define SWIFT_DEFAULT_GRAVITY_DEBUG_H
+
+__attribute__((always_inline)) INLINE static void gravity_debug_particle(
+    const struct gpart* p) {
+  printf(
+      "mass=%.3e time_bin=%d\n"
+      "x=[%.5e,%.5e,%.5e], v_full=[%.5e,%.5e,%.5e], a=[%.5e,%.5e,%.5e]\n",
+      p->mass, p->time_bin, p->x[0], p->x[1], p->x[2], p->v_full[0],
+      p->v_full[1], p->v_full[2], p->a_grav[0], p->a_grav[1], p->a_grav[2]);
+#ifdef SWIFT_DEBUG_CHECKS
+  printf("num_interacted=%lld ti_drift=%lld ti_kick=%lld\n", p->num_interacted,
+         p->ti_drift, p->ti_kick);
+#endif
+}
+
+#endif /* SWIFT_DEFAULT_GRAVITY_DEBUG_H */

src/gravity/Potential/gravity_iact.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_DEFAULT_GRAVITY_IACT_H
+#define SWIFT_DEFAULT_GRAVITY_IACT_H
+
+/* Includes. */
+#include "kernel_gravity.h"
+#include "kernel_long_gravity.h"
+#include "multipole.h"
+
+/* Standard headers */
+#include <float.h>
+
+/**
+ * @brief Computes the intensity of the force at a point generated by a
+ * point-mass.
+ *
+ * The returned quantity needs to be multiplied by the distance vector to obtain
+ * the force vector.
+ *
+ * @param r2 Square of the distance to the point-mass.
+ * @param h2 Square of the softening length.
+ * @param h_inv Inverse of the softening length.
+ * @param h_inv3 Cube of the inverse of the softening length.
+ * @param mass Mass of the point-mass.
+ * @param f_ij (return) The force intensity.
+ * @param pot_ij (return) The potential.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
+    const float r2, const float h2, const float h_inv, const float h_inv3,
+    const float mass, float *f_ij, float *pot_ij) {
+
+  /* Get the inverse distance */
+  const float r_inv = 1.f / sqrtf(r2 + FLT_MIN);
+
+  /* Should we soften ? */
+  if (r2 >= h2) {
+
+    /* Get Newtonian gravity */
+    *f_ij = mass * r_inv * r_inv * r_inv;
+    *pot_ij = -mass * r_inv;
+
+  } else {
+
+    const float r = r2 * r_inv;
+    const float ui = r * h_inv;
+
+    float W_f_ij, W_pot_ij;
+    kernel_grav_force_eval(ui, &W_f_ij);
+    kernel_grav_pot_eval(ui, &W_pot_ij);
+
+    /* Get softened gravity */
+    *f_ij = mass * h_inv3 * W_f_ij;
+    *pot_ij = mass * h_inv * W_pot_ij;
+  }
+}
+
+/**
+ * @brief Computes the intensity of the force at a point generated by a
+ * point-mass truncated for long-distance periodicity.
+ *
+ * The returned quantity needs to be multiplied by the distance vector to obtain
+ * the force vector.
+ *
+ * @param r2 Square of the distance to the point-mass.
+ * @param h2 Square of the softening length.
+ * @param h_inv Inverse of the softening length.
+ * @param h_inv3 Cube of the inverse of the softening length.
+ * @param mass Mass of the point-mass.
+ * @param r_s_inv Inverse of the mesh smoothing scale.
+ * @param f_ij (return) The force intensity.
+ * @param pot_ij (return) The potential.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
+    const float r2, const float h2, const float h_inv, const float h_inv3,
+    const float mass, const float r_s_inv, float *f_ij, float *pot_ij) {
+
+  /* Get the inverse distance */
+  const float r_inv = 1.f / sqrtf(r2 + FLT_MIN);
+  const float r = r2 * r_inv;
+
+  /* Should we soften ? */
+  if (r2 >= h2) {
+
+    /* Get Newtonian gravity */
+    *f_ij = mass * r_inv * r_inv * r_inv;
+    *pot_ij = -mass * r_inv;
+
+  } else {
+
+    const float ui = r * h_inv;
+    float W_f_ij, W_pot_ij;
+
+    kernel_grav_force_eval(ui, &W_f_ij);
+    kernel_grav_pot_eval(ui, &W_pot_ij);
+
+    /* Get softened gravity */
+    *f_ij = mass * h_inv3 * W_f_ij;
+    *pot_ij = mass * h_inv * W_pot_ij;
+  }
+
+  /* Get long-range correction */
+  const float u_lr = r * r_s_inv;
+  float corr_f_lr, corr_pot_lr;
+  kernel_long_grav_force_eval(u_lr, &corr_f_lr);
+  kernel_long_grav_pot_eval(u_lr, &corr_pot_lr);
+  *f_ij *= corr_f_lr;
+  *pot_ij *= corr_pot_lr;
+}
+
+/**
+ * @brief Computes the forces at a point generated by a multipole.
+ *
+ * This assumes M_100 == M_010 == M_001 == 0.
+ * This uses the quadrupole and trace of the octupole terms only and defaults to
+ * the monopole if the code is compiled with low-order gravity only.
+ *
+ * @param r_x x-component of the distance vector to the multipole.
+ * @param r_y y-component of the distance vector to the multipole.
+ * @param r_z z-component of the distance vector to the multipole.
+ * @param r2 Square of the distance vector to the multipole.
+ * @param h The softening length.
+ * @param h_inv Inverse of the softening length.
+ * @param m The multipole.
+ * @param f_x (return) The x-component of the acceleration.
+ * @param f_y (return) The y-component of the acceleration.
+ * @param f_z (return) The z-component of the acceleration.
+ * @param pot (return) The potential.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
+    const float r_x, const float r_y, const float r_z, const float r2,
+    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.
+ * 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
+
+  float f_ij, pot_ij;
+  runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,
+                           &f_ij, &pot_ij);
+  *f_x = f_ij * r_x;
+  *f_y = f_ij * r_y;
+  *f_z = f_ij * r_z;
+  *pot = pot_ij;
+
+#else
+
+  /* Get the inverse distance */
+  const float r_inv = 1.f / sqrtf(r2);
+
+  /* Compute the derivatives of the potential */
+  struct potential_derivatives_M2P d;
+  compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 0, 0.f,
+                                    &d);
+
+  /* 0th order contributions */
+  *f_x = m->M_000 * d.D_100;
+  *f_y = m->M_000 * d.D_010;
+  *f_z = m->M_000 * d.D_001;
+  *pot = m->M_000 * d.D_000;
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+
+/* 1st order contributions */
+
+/* 1st order contributions are all 0 since the dipole is 0 */
+
+/* *f_x = m->M_001 * d.D_101 + m->M_010 * d.D_110 + m->M_100 * d.D_200 ; */
+/* *f_y = m->M_001 * d.D_011 + m->M_010 * d.D_020 + m->M_100 * d.D_110 ; */
+/* *f_z = m->M_001 * d.D_002 + m->M_010 * d.D_011 + m->M_100 * d.D_101 ; */
+/* *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
+
+  /* 2nd order contributions */
+  *f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
+          m->M_101 * d.D_201 + m->M_110 * d.D_210 + m->M_200 * d.D_300;
+  *f_y += m->M_002 * d.D_012 + m->M_011 * d.D_021 + m->M_020 * d.D_030 +
+          m->M_101 * d.D_111 + m->M_110 * d.D_120 + m->M_200 * d.D_210;
+  *f_z += m->M_002 * d.D_003 + m->M_011 * d.D_012 + m->M_020 * d.D_021 +
+          m->M_101 * d.D_102 + m->M_110 * d.D_111 + m->M_200 * d.D_201;
+  *pot += m->M_002 * d.D_002 + m->M_011 * d.D_011 + m->M_020 * d.D_020 +
+          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
+
+  /* 3rd order contributions */
+  *f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
+          m->M_030 * d.D_130 + m->M_102 * d.D_202 + m->M_111 * d.D_211 +
+          m->M_120 * d.D_220 + m->M_201 * d.D_301 + m->M_210 * d.D_310 +
+          m->M_300 * d.D_400;
+  *f_y += m->M_003 * d.D_013 + m->M_012 * d.D_022 + m->M_021 * d.D_031 +
+          m->M_030 * d.D_040 + m->M_102 * d.D_112 + m->M_111 * d.D_121 +
+          m->M_120 * d.D_130 + m->M_201 * d.D_211 + m->M_210 * d.D_220 +
+          m->M_300 * d.D_310;
+  *f_z += m->M_003 * d.D_004 + m->M_012 * d.D_013 + m->M_021 * d.D_022 +
+          m->M_030 * d.D_031 + m->M_102 * d.D_103 + m->M_111 * d.D_112 +
+          m->M_120 * d.D_121 + m->M_201 * d.D_202 + m->M_210 * d.D_211 +
+          m->M_300 * d.D_301;
+  *pot += m->M_003 * d.D_003 + m->M_012 * d.D_012 + m->M_021 * d.D_021 +
+          m->M_030 * d.D_030 + m->M_102 * d.D_102 + m->M_111 * d.D_111 +
+          m->M_120 * d.D_120 + m->M_201 * d.D_201 + m->M_210 * d.D_210 +
+          m->M_300 * d.D_300;
+
+#endif
+
+  /* Take care of the the sign convention */
+  *f_x *= -1.f;
+  *f_y *= -1.f;
+  *f_z *= -1.f;
+  *pot *= -1.f;
+#endif
+}
+
+/**
+ * @brief Computes the forces at a point generated by a multipole, truncated for
+ * long-range periodicity.
+ *
+ * This assumes M_100 == M_010 == M_001 == 0.
+ * This uses the quadrupole term and trace of the octupole terms only and
+ * defaults to the monopole if the code is compiled with low-order gravity only.
+ *
+ * @param r_x x-component of the distance vector to the multipole.
+ * @param r_y y-component of the distance vector to the multipole.
+ * @param r_z z-component of the distance vector to the multipole.
+ * @param r2 Square of the distance vector to the multipole.
+ * @param h The softening length.
+ * @param h_inv Inverse of the softening length.
+ * @param r_s_inv The inverse of the gravity mesh-smoothing scale.
+ * @param m The multipole.
+ * @param f_x (return) The x-component of the acceleration.
+ * @param f_y (return) The y-component of the acceleration.
+ * @param f_z (return) The z-component of the acceleration.
+ * @param pot (return) The potential.
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
+    const float r_x, const float r_y, const float r_z, const float r2,
+    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.
+ * 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
+
+  float f_ij, pot_ij;
+  runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
+                                m->M_000, r_s_inv, &f_ij, &pot_ij);
+  *f_x = f_ij * r_x;
+  *f_y = f_ij * r_y;
+  *f_z = f_ij * r_z;
+  *pot = -pot_ij;
+
+#else
+
+  /* Get the inverse distance */
+  const float r_inv = 1.f / sqrtf(r2);
+
+  /* Compute the derivatives of the potential */
+  struct potential_derivatives_M2P d;
+  compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 1,
+                                    r_s_inv, &d);
+
+  /* 0th order contributions */
+  *f_x = m->M_000 * d.D_100;
+  *f_y = m->M_000 * d.D_010;
+  *f_z = m->M_000