From 2739c50c1e953085ea2291710ddef3c5a8afb17d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 11 Dec 2024 17:07:19 +0100 Subject: [PATCH] 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 | 2 +- src/core/electrostatics/elc.cpp | 8 +- src/core/electrostatics/p3m.cpp | 4 +- src/core/magnetostatics/dp3m.cpp | 9 +- src/core/p3m/common.hpp | 11 ++- src/core/p3m/data_struct.hpp | 2 +- src/core/p3m/influence_function.hpp | 73 ++++++++++----- src/core/p3m/influence_function_dipolar.hpp | 98 +++++++++++++-------- src/core/unit_tests/p3m_test.cpp | 4 +- 9 files changed, 131 insertions(+), 80 deletions(-) diff --git a/doc/doxygen/Doxyfile.in b/doc/doxygen/Doxyfile.in index 67333f21af..b6c0a0e1b2 100644 --- a/doc/doxygen/Doxyfile.in +++ b/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 diff --git a/src/core/electrostatics/elc.cpp b/src/core/electrostatics/elc.cpp index a7988e84db..4e28654ae1 100644 --- a/src/core/electrostatics/elc.cpp +++ b/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}; diff --git a/src/core/electrostatics/p3m.cpp b/src/core/electrostatics/p3m.cpp index 0826e3871b..99edb97c4c 100644 --- a/src/core/electrostatics/p3m.cpp +++ b/src/core/electrostatics/p3m.cpp @@ -127,7 +127,7 @@ void CoulombP3MImpl::count_charged_particles() { template void CoulombP3MImpl::calc_influence_function_force() { auto const [KX, KY, KZ] = p3m.fft->get_permutations(); - p3m.g_force = grid_influence_function( + p3m.g_force = grid_influence_function( p3m.params, p3m.mesh.start, p3m.mesh.stop, KX, KY, KZ, get_system().box_geo->length_inv()); } @@ -138,7 +138,7 @@ void CoulombP3MImpl::calc_influence_function_force() { template void CoulombP3MImpl::calc_influence_function_energy() { auto const [KX, KY, KZ] = p3m.fft->get_permutations(); - p3m.g_energy = grid_influence_function( + p3m.g_energy = grid_influence_function( p3m.params, p3m.mesh.start, p3m.mesh.stop, KX, KY, KZ, get_system().box_geo->length_inv()); } diff --git a/src/core/magnetostatics/dp3m.cpp b/src/core/magnetostatics/dp3m.cpp index f59613943c..1310ff24a7 100644 --- a/src/core/magnetostatics/dp3m.cpp +++ b/src/core/magnetostatics/dp3m.cpp @@ -116,8 +116,9 @@ double DipolarP3MImpl::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( + 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::calc_surface_term( template void DipolarP3MImpl::calc_influence_function_force() { - dp3m.g_force = grid_influence_function( + dp3m.g_force = grid_influence_function( dp3m.params, dp3m.mesh.start, dp3m.mesh.stop, get_system().box_geo->length_inv()); } template void DipolarP3MImpl::calc_influence_function_energy() { - dp3m.g_energy = grid_influence_function( + dp3m.g_energy = grid_influence_function( dp3m.params, dp3m.mesh.start, dp3m.mesh.stop, get_system().box_geo->length_inv()); } diff --git a/src/core/p3m/common.hpp b/src/core/p3m/common.hpp index 3cb9c1cb45..1f0a0a1739 100644 --- a/src/core/p3m/common.hpp +++ b/src/core/p3m/common.hpp @@ -53,6 +53,7 @@ auto constexpr P3M_EPSILON_METALLIC = 0.0; #include #include +/** @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 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, 3> inline calc_meshift( +std::array, 3> inline calc_p3m_mesh_shift( Utils::Vector3i const &mesh_size, bool zero_out_midpoint = false) { std::array, 3> ret{}; @@ -262,4 +262,3 @@ std::array, 3> inline calc_meshift( return ret; } -} // namespace detail diff --git a/src/core/p3m/data_struct.hpp b/src/core/p3m/data_struct.hpp index 9bd81ad7c2..35e1cffc2f 100644 --- a/src/core/p3m/data_struct.hpp +++ b/src/core/p3m/data_struct.hpp @@ -66,7 +66,7 @@ template 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) */ diff --git a/src/core/p3m/influence_function.hpp b/src/core/p3m/influence_function.hpp index c3373a88d1..c3487578cb 100644 --- a/src/core/p3m/influence_function.hpp +++ b/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 /** - * @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 -double G_opt(int cao, double alpha, Utils::Vector3d const &k, - Utils::Vector3d const &h) { +double G_opt(P3MParameters const ¶ms, 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(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 +template std::vector grid_influence_function( P3MParameters const ¶ms, 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 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(params.cao, params.alpha, k, params.a)); + g[index] = FloatType(G_opt(params, k)); } ++index; }); return g; } + +#endif // defined(P3M) diff --git a/src/core/p3m/influence_function_dipolar.hpp b/src/core/p3m/influence_function_dipolar.hpp index 74af200a20..d08a807069 100644 --- a/src/core/p3m/influence_function_dipolar.hpp +++ b/src/core/p3m/influence_function_dipolar.hpp @@ -39,31 +39,54 @@ #include #include -/** Calculate the aliasing sums for the optimal influence function. +/** + * @brief Calculate the aliasing sums for the optimal influence function. + * + * Evaluate the aliasing sums in the numerator and denominator of + * the expression for the optimal influence function, see + * @cite hockney88a eq (8-22), p. 275. + * + * The full equation is: + * @f{align*}{ + * G_{\text{opt}}(\vec{n}, S, \text{cao}) &= + * \displaystyle\frac{1}{\sum_\vec{m} U^2} + * \sum_\vec{m}U^2 + * \displaystyle\frac{\left(\vec{d}_{\text{op}}[\vec{n}](\vec{s}[\vec{n}] + + * \vec{m}\odot\vec{N})\right)^S} + * {\left(\vec{s}[\vec{n}] + \vec{m}\odot\vec{N}\right)^2} + * e^{-\pi^2\left(\alpha\vec{N}\odot\vec{a}\right)^{-2}\left(\vec{s}[\vec{n}] + * +\vec{m}\odot\vec{N}\right)^2} \\ * - * Calculates the aliasing sums in the numerator and denominator of - * the expression for the optimal influence function (see - * @cite hockney88a : 8-22, p. 275). + * U &= \operatorname{det}\left[I_3\cdot\operatorname{sinc}\left( + * \vec{s}[\vec{n}]\oslash{\vec{N}}+\vec{m}\odot\vec{N}\right) + * \right]^{\text{cao}} + * @f} * - * \tparam S order (2 for energy, 3 for forces) - * \param params DP3M parameters - * \param shift shift for a given n-vector - * \param d_op differential operator for a given n-vector - * \return The result of the fraction. + * 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, + * @f$ \oslash @f$ the Hadamard division. + * + * @tparam S Order of the differential operator (2 for energy, 3 for forces) + * @tparam m Number of Brillouin zones that contribute to the aliasing sums + * + * @param params DP3M parameters + * @param shift shifting integer vector for a specific k-vector + * @param d_op differential operator for a specific k-vector + * @return The optimal influence function for a specific k-vector. */ -template +template double G_opt_dipolar(P3MParameters const ¶ms, Utils::Vector3i const &shift, Utils::Vector3i const &d_op) { - auto constexpr limit = P3M_BRILLOUIN; auto constexpr exp_limit = 30.; - auto constexpr m_start = Utils::Vector3i::broadcast(-limit); - auto constexpr m_stop = Utils::Vector3i::broadcast(limit + 1); + auto constexpr m_start = Utils::Vector3i::broadcast(-m); + auto constexpr m_stop = Utils::Vector3i::broadcast(m + 1); auto const cao = params.cao; auto const mesh = params.mesh[0]; auto const offset = static_cast(shift) / static_cast(mesh); - auto const f2 = Utils::sqr(std::numbers::pi / params.alpha_L); + auto const exp_prefactor = Utils::sqr(std::numbers::pi / params.alpha_L); auto indices = Utils::Vector3i{}; auto nm = Utils::Vector3i{}; auto fnm = Utils::Vector3d{}; @@ -73,14 +96,14 @@ double G_opt_dipolar(P3MParameters const ¶ms, Utils::Vector3i const &shift, for_each_3d( m_start, m_stop, indices, [&]() { - auto const norm_sq = nm.norm2(); - auto const sz = std::pow(Utils::product(fnm), 2 * cao); - auto const exp_term = f2 * norm_sq; + auto const U2 = std::pow(Utils::product(fnm), 2 * cao); + auto const nm2 = nm.norm2(); + auto const exp_term = exp_prefactor * nm2; if (exp_term < exp_limit) { - auto const f3 = sz * std::exp(-exp_term) / norm_sq; + auto const f3 = U2 * std::exp(-exp_term) / nm2; numerator += f3 * Utils::int_pow(d_op * nm); } - denominator += sz; + denominator += U2; }, [&](unsigned dim, int n) { nm[dim] = shift[dim] + n * mesh; @@ -98,14 +121,15 @@ double G_opt_dipolar(P3MParameters const ¶ms, Utils::Vector3i const &shift, * over a regular grid of k vectors, and returns the values as a vector. * * @tparam S Order of the differential operator, e.g. 2 for energy, 3 for force + * @tparam m Number of Brillouin zones that contribute to the aliasing sums * * @param params DP3M 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 inv_box_l Inverse box length * @return Values of the influence function at regular grid points. */ -template +template std::vector grid_influence_function( P3MParameters const ¶ms, Utils::Vector3i const &n_start, Utils::Vector3i const &n_stop, Utils::Vector3d const &inv_box_l) { @@ -124,8 +148,8 @@ std::vector grid_influence_function( auto prefactor = Utils::int_pow<3>(static_cast(params.mesh[0])) * 2. * Utils::int_pow<2>(inv_box_l[0]); - auto const offset = detail::calc_meshift(params.mesh, false)[0]; - auto const d_op = detail::calc_meshift(params.mesh, true)[0]; + auto const offset = calc_p3m_mesh_shift(params.mesh, false)[0]; + auto const d_op = calc_p3m_mesh_shift(params.mesh, true)[0]; auto const half_mesh = params.mesh[0] / 2; auto indices = Utils::Vector3i{}; auto shift_off = Utils::Vector3i{}; @@ -137,8 +161,8 @@ std::vector grid_influence_function( [&]() { if (((indices[0] % half_mesh != 0) or (indices[1] % half_mesh != 0) or (indices[2] % half_mesh != 0))) { - g[index] = FloatType(prefactor * - G_opt_dipolar(params, shift_off, d_op_off)); + g[index] = FloatType( + prefactor * G_opt_dipolar(params, shift_off, d_op_off)); } ++index; }, @@ -150,12 +174,12 @@ std::vector grid_influence_function( return g; } +template inline double G_opt_dipolar_self_energy(P3MParameters const ¶ms, Utils::Vector3i const &shift) { - auto constexpr limit = P3M_BRILLOUIN + 1; - auto constexpr m_start = Utils::Vector3i::broadcast(-limit); - auto constexpr m_stop = Utils::Vector3i::broadcast(limit + 1); + auto constexpr m_start = Utils::Vector3i::broadcast(-m); + auto constexpr m_stop = Utils::Vector3i::broadcast(m + 1); auto const cao = params.cao; auto const mesh = params.mesh[0]; auto const offset = @@ -177,19 +201,21 @@ inline double G_opt_dipolar_self_energy(P3MParameters const ¶ms, /** * @brief Calculate self-energy of the influence function. * + * @tparam m Number of Brillouin zones that contribute to the aliasing sums + * * @param params DP3M 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 g Energies on the grid. * @return Total self-energy. */ -template +template inline double grid_influence_function_self_energy( P3MParameters const ¶ms, Utils::Vector3i const &n_start, Utils::Vector3i const &n_stop, std::vector const &g) { - auto const offset = detail::calc_meshift(params.mesh, false)[0]; - auto const d_op = detail::calc_meshift(params.mesh, true)[0]; + auto const offset = calc_p3m_mesh_shift(params.mesh, false)[0]; + auto const d_op = calc_p3m_mesh_shift(params.mesh, true)[0]; auto const half_mesh = params.mesh[0] / 2; auto indices = Utils::Vector3i{}; auto shift_off = Utils::Vector3i{}; @@ -202,8 +228,8 @@ inline double grid_influence_function_self_energy( [&]() { if (((indices[0] % half_mesh != 0) or (indices[1] % half_mesh != 0) or (indices[2] % half_mesh != 0))) { - auto const U2 = G_opt_dipolar_self_energy(params, shift_off); - energy += double(g[index]) * U2 * d_op_off.norm2(); + auto const U2 = G_opt_dipolar_self_energy(params, shift_off); + energy += static_cast(g[index]) * U2 * d_op_off.norm2(); } ++index; }, diff --git a/src/core/unit_tests/p3m_test.cpp b/src/core/unit_tests/p3m_test.cpp index 3f091d565c..eea67334c0 100644 --- a/src/core/unit_tests/p3m_test.cpp +++ b/src/core/unit_tests/p3m_test.cpp @@ -35,7 +35,7 @@ BOOST_AUTO_TEST_CASE(calc_meshift_false) { std::vector{0, 1, 2, 3, -3, -2, -1}}}; auto const mesh = Utils::Vector3i{{1, 4, 7}}; - auto const val = detail::calc_meshift(mesh, false); + auto const val = calc_p3m_mesh_shift(mesh, false); for (std::size_t i = 0u; i < 3u; ++i) { for (std::size_t j = 0u; j < ref[i].size(); ++j) { @@ -50,7 +50,7 @@ BOOST_AUTO_TEST_CASE(calc_meshift_true) { std::vector{0, 1, 2, 0, -3, -2, -1}}}; auto const mesh = Utils::Vector3i{{1, 4, 7}}; - auto const val = detail::calc_meshift(mesh, true); + auto const val = calc_p3m_mesh_shift(mesh, true); for (std::size_t i = 0u; i < 3u; ++i) { for (std::size_t j = 0u; j < ref[i].size(); ++j) {