Refactor P3M code

Use the same variable names in P3M and DP3M influence functions,
move their Brillouin limit to a template parameter, and improve
their documentation. Avoid division by zero in ELC with constant
potential.

doc/doxygen/Doxyfile.in

@@ -1511,7 +1511,7 @@ MATHJAX_RELPATH        = https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/
 # MATHJAX_EXTENSIONS = TeX/AMSmath TeX/AMSsymbols
 # This tag requires that the tag USE_MATHJAX is set to YES.
 
-MATHJAX_EXTENSIONS     =
+MATHJAX_EXTENSIONS     = TeX/AMSmath
 
 # The MATHJAX_CODEFILE tag can be used to specify a file with javascript pieces
 # of code that will be used on startup of the MathJax code. See the MathJax site

src/core/electrostatics/elc.cpp

@@ -368,10 +368,6 @@ ElectrostaticLayerCorrection::z_energy(ParticleRange const &particles) const {
   auto const &box_geo = *get_system().box_geo;
   auto const xy_area_inv = box_geo.length_inv()[0] * box_geo.length_inv()[1];
   auto const pref = prefactor * 2. * std::numbers::pi * xy_area_inv;
-  auto const delta = elc.delta_mid_top * elc.delta_mid_bot;
-  auto const fac_delta_mid_bot = elc.delta_mid_bot / (1. - delta);
-  auto const fac_delta_mid_top = elc.delta_mid_top / (1. - delta);
-  auto const fac_delta = delta / (1. - delta);
 
   /* for non-neutral systems, this shift gives the background contribution
    * (rsp. for this shift, the DM of the background is zero) */
@@ -397,6 +393,10 @@ ElectrostaticLayerCorrection::z_energy(ParticleRange const &particles) const {
       }
     } else {
       // metallic boundaries
+      auto const delta = elc.delta_mid_top * elc.delta_mid_bot;
+      auto const fac_delta_mid_bot = elc.delta_mid_bot / (1. - delta);
+      auto const fac_delta_mid_top = elc.delta_mid_top / (1. - delta);
+      auto const fac_delta = delta / (1. - delta);
       clear_vec(gblcblk, size);
       auto const h = elc.box_h;
       ImageSum const image_sum{delta, shift, lz};

src/core/electrostatics/p3m.cpp

@@ -127,7 +127,7 @@ void CoulombP3MImpl<FloatType, Architecture>::count_charged_particles() {
 template <typename FloatType, Arch Architecture>
 void CoulombP3MImpl<FloatType, Architecture>::calc_influence_function_force() {
   auto const [KX, KY, KZ] = p3m.fft->get_permutations();
-  p3m.g_force = grid_influence_function<FloatType, 1>(
+  p3m.g_force = grid_influence_function<FloatType, 1, P3M_BRILLOUIN>(
       p3m.params, p3m.mesh.start, p3m.mesh.stop, KX, KY, KZ,
       get_system().box_geo->length_inv());
 }
@@ -138,7 +138,7 @@ void CoulombP3MImpl<FloatType, Architecture>::calc_influence_function_force() {
 template <typename FloatType, Arch Architecture>
 void CoulombP3MImpl<FloatType, Architecture>::calc_influence_function_energy() {
   auto const [KX, KY, KZ] = p3m.fft->get_permutations();
-  p3m.g_energy = grid_influence_function<FloatType, 0>(
+  p3m.g_energy = grid_influence_function<FloatType, 0, P3M_BRILLOUIN>(
       p3m.params, p3m.mesh.start, p3m.mesh.stop, KX, KY, KZ,
       get_system().box_geo->length_inv());
 }

src/core/magnetostatics/dp3m.cpp

@@ -116,8 +116,9 @@ double
 DipolarP3MImpl<FloatType, Architecture>::calc_average_self_energy_k_space()
     const {
   auto const &box_geo = *get_system().box_geo;
-  auto const node_phi = grid_influence_function_self_energy(
-      dp3m.params, dp3m.mesh.start, dp3m.mesh.stop, dp3m.g_energy);
+  auto const node_phi =
+      grid_influence_function_self_energy<FloatType, P3M_BRILLOUIN>(
+          dp3m.params, dp3m.mesh.start, dp3m.mesh.stop, dp3m.g_energy);
 
   double phi = 0.;
   boost::mpi::reduce(comm_cart, node_phi, phi, std::plus<>(), 0);
@@ -519,14 +520,14 @@ double DipolarP3MImpl<FloatType, Architecture>::calc_surface_term(
 
 template <typename FloatType, Arch Architecture>
 void DipolarP3MImpl<FloatType, Architecture>::calc_influence_function_force() {
-  dp3m.g_force = grid_influence_function<FloatType, 3>(
+  dp3m.g_force = grid_influence_function<FloatType, 3, P3M_BRILLOUIN>(
       dp3m.params, dp3m.mesh.start, dp3m.mesh.stop,
       get_system().box_geo->length_inv());
 }
 
 template <typename FloatType, Arch Architecture>
 void DipolarP3MImpl<FloatType, Architecture>::calc_influence_function_energy() {
-  dp3m.g_energy = grid_influence_function<FloatType, 2>(
+  dp3m.g_energy = grid_influence_function<FloatType, 2, P3M_BRILLOUIN>(
       dp3m.params, dp3m.mesh.start, dp3m.mesh.stop,
       get_system().box_geo->length_inv());
 }

src/core/p3m/common.hpp

@@ -53,6 +53,7 @@ auto constexpr P3M_EPSILON_METALLIC = 0.0;
 #include <span>
 #include <stdexcept>
 
+/** @brief P3M kernel architecture. */
 enum class Arch { CPU, GPU };
 
 /** @brief Structure to hold P3M parameters and some dependent variables. */
@@ -65,8 +66,8 @@ struct P3MParameters {
   /** cutoff radius for real space electrostatics (>0), rescaled to
    *  @p r_cut_iL = @p r_cut * @p box_l_i. */
   double r_cut_iL;
-  /** number of mesh points per coordinate direction (>0). */
-  Utils::Vector3i mesh = {};
+  /** number of mesh points per coordinate direction (>0), in real space. */
+  Utils::Vector3i mesh;
   /** offset of the first mesh point (lower left corner) from the
    *  coordinate origin ([0,1[). */
   Utils::Vector3d mesh_off;
@@ -237,15 +238,14 @@ template <typename FloatType> struct P3MFFTMesh {
 
 #endif // defined(P3M) or defined(DP3M)
 
-namespace detail {
-/** Calculate indices that shift @ref P3MParameters::mesh "mesh" by `mesh/2`.
+/** @brief Calculate indices that shift @ref P3MParameters::mesh by `mesh/2`.
  *  For each mesh size @f$ n @f$ in @c mesh_size, create a sequence of integer
  *  values @f$ \left( 0, \ldots, \lfloor n/2 \rfloor, -\lfloor n/2 \rfloor,
  *  \ldots, -1\right) @f$ if @c zero_out_midpoint is false, otherwise
  *  @f$ \left( 0, \ldots, \lfloor n/2 - 1 \rfloor, 0, -\lfloor n/2 \rfloor,
  *  \ldots, -1\right) @f$.
  */
-std::array<std::vector<int>, 3> inline calc_meshift(
+std::array<std::vector<int>, 3> inline calc_p3m_mesh_shift(
     Utils::Vector3i const &mesh_size, bool zero_out_midpoint = false) {
   std::array<std::vector<int>, 3> ret{};
 
@@ -262,4 +262,3 @@ std::array<std::vector<int>, 3> inline calc_meshift(
 
   return ret;
 }
-} // namespace detail

src/core/p3m/data_struct.hpp

@@ -66,7 +66,7 @@ template <typename FloatType> struct p3m_data_struct {
    *  i.e. the prefactor @f$ 2i\pi/L @f$ is missing!
    */
   void calc_differential_operator() {
-    d_op = detail::calc_meshift(params.mesh, true);
+    d_op = calc_p3m_mesh_shift(params.mesh, true);
   }
 
   /** @brief Force optimised influence function (k-space) */

src/core/p3m/influence_function.hpp

@@ -19,6 +19,10 @@
 
 #pragma once
 
+#include "config/config.hpp"
+
+#if defined(P3M)
+
 #include "p3m/common.hpp"
 #include "p3m/for_each_3d.hpp"
 #include "p3m/math.hpp"
@@ -36,34 +40,54 @@
 #include <vector>
 
 /**
- * @brief Hockney/Eastwood/Ballenegger optimal influence function.
+ * @brief Calculate the aliasing sums for the optimal influence function.
  *
  * This implements Eq. 30 of @cite cerda08d, which can be used
  * for monopole and dipole P3M by choosing the appropriate S factor.
  *
- * @tparam S Order of the differential operator, e.g. 0 for potential,
- *          1 for electric field, ...
- * @tparam m Number of aliasing terms to take into account.
+ * The full equation is:
+ * @f{align*}{
+ * G_{\text{opt}}(\vec{n}, S, \text{cao}) &=
+ *   \displaystyle\frac{4\pi}{\sum_\vec{m} U^2}
+ *   \sum_\vec{m}U^2
+ *   \displaystyle\frac{\left(\vec{k}\odot\vec{k}_m\right)^S}
+ *                     {|\vec{k}_m|^2}
+ *   e^{-1/(2\alpha)^2|\vec{k}_m|^2} \\
+ *
+ * U &= \operatorname{det}\left[I_3\cdot\operatorname{sinc}\left(
+ *          \frac{\vec{k}_m\odot\vec{a}}{2\pi}\right)\right]^{\text{cao}} \\
+ *
+ * \vec{k} &= \frac{2\pi}{\vec{N}\odot\vec{a}}\vec{s}[\vec{n}] \\
  *
- * @param cao Charge assignment order.
- * @param alpha Ewald splitting parameter.
- * @param k k-vector to evaluate the function for.
- * @param h Grid spacing.
+ * \vec{k}_m &= \frac{2\pi}{\vec{N}\odot\vec{a}}
+ *              \left(\vec{s}[\vec{n}]+\vec{m}\odot\vec{N}\right)
+ * @f}
+ *
+ * with @f$ I_3 @f$ the 3x3 identity matrix, @f$ \vec{N} @f$ the global mesh
+ * size in k-space coordinates, @f$ \vec{a} @f$ the mesh size, @f$ \vec{m} @f$
+ * the Brillouin zone coordinates, @f$ \odot @f$ the Hadamard product.
+ *
+ * @tparam S Order of the differential operator (0 for potential, 1 for E-field)
+ * @tparam m Number of Brillouin zones that contribute to the aliasing sums
+ *
+ * @param params   P3M parameters
+ * @param k        k-vector to evaluate the function for.
+ * @return The optimal influence function for a specific k-vector.
  */
 template <std::size_t S, std::size_t m>
-double G_opt(int cao, double alpha, Utils::Vector3d const &k,
-             Utils::Vector3d const &h) {
+double G_opt(P3MParameters const &params, Utils::Vector3d const &k) {
 
   auto const k2 = k.norm2();
   if (k2 == 0.) {
     return 0.;
   }
 
-  auto constexpr limit = 30.;
+  auto constexpr exp_limit = 30.;
   auto constexpr m_start = Utils::Vector3i::broadcast(-m);
   auto constexpr m_stop = Utils::Vector3i::broadcast(m + 1);
-  auto const exponent_prefactor = Utils::sqr(1. / (2. * alpha));
-  auto const wavevector = (2. * std::numbers::pi) / h;
+  auto const cao = params.cao;
+  auto const exponent_prefactor = Utils::sqr(1. / (2. * params.alpha));
+  auto const wavevector = (2. * std::numbers::pi) / params.a;
   auto const wavevector_i = 1. / wavevector;
   auto indices = Utils::Vector3i{};
   auto km = Utils::Vector3d{};
@@ -76,9 +100,9 @@ double G_opt(int cao, double alpha, Utils::Vector3d const &k,
       [&]() {
         auto const U2 = std::pow(Utils::product(fnm), 2 * cao);
         auto const km2 = km.norm2();
-        auto const exponent = exponent_prefactor * km2;
-        if (exponent < limit) {
-          auto const f3 = std::exp(-exponent) * (4. * std::numbers::pi / km2);
+        auto const exp_term = exponent_prefactor * km2;
+        if (exp_term < exp_limit) {
+          auto const f3 = std::exp(-exp_term) * (4. * std::numbers::pi / km2);
           numerator += U2 * f3 * Utils::int_pow<S>(k * km);
         }
         denominator += U2;
@@ -97,26 +121,25 @@ double G_opt(int cao, double alpha, Utils::Vector3d const &k,
  * This evaluates the optimal influence function @ref G_opt
  * over a regular grid of k vectors, and returns the values as a vector.
  *
- * @tparam S Order of the differential operator, e.g. 0 for potential,
- *          1 for electric field...
- * @tparam m Number of aliasing terms to take into account.
+ * @tparam S Order of the differential operator (0 for potential, 1 for E-field)
+ * @tparam m Number of Brillouin zones that contribute to the aliasing sums
  *
  * @param params P3M parameters.
- * @param n_start Lower left corner of the grid.
- * @param n_stop Upper right corner of the grid.
+ * @param n_start Lower left corner of the grid in k-space.
+ * @param n_stop Upper right corner of the grid in k-space.
  * @param KX k-space x-axis index.
  * @param KY k-space y-axis index.
  * @param KZ k-space z-axis index.
  * @param inv_box_l Inverse box length.
  * @return Values of G_opt at regular grid points.
  */
-template <typename FloatType, std::size_t S, std::size_t m = 0>
+template <typename FloatType, std::size_t S, std::size_t m>
 std::vector<FloatType> grid_influence_function(
     P3MParameters const &params, Utils::Vector3i const &n_start,
     Utils::Vector3i const &n_stop, int const KX, int const KY, int const KZ,
     Utils::Vector3d const &inv_box_l) {
 
-  auto const shifts = detail::calc_meshift(params.mesh);
+  auto const shifts = calc_p3m_mesh_shift(params.mesh);
   auto const size = n_stop - n_start;
 
   /* The influence function grid */
@@ -141,10 +164,12 @@ std::vector<FloatType> grid_influence_function(
           Utils::Vector3d{{shifts[0u][indices[KX]] * wavevector[0u],
                            shifts[1u][indices[KY]] * wavevector[1u],
                            shifts[2u][indices[KZ]] * wavevector[2u]}};
-      g[index] = FloatType(G_opt<S, m>(params.cao, params.alpha, k, params.a));
+      g[index] = FloatType(G_opt<S, m>(params, k));
     }
     ++index;
   });
 
   return g;
 }
+
+#endif // defined(P3M)