Implement dipole fields calculation with DDS

Introduce dipolar direct sum CPU kernels to compute the total dipole
field on magnetic particles.

Co-authored-by: Jean-Noël Grad <jgrad@icp.uni-stuttgart.de>

doc/sphinx/installation.rst

@@ -430,6 +430,9 @@ General features
 
 -  ``DIPOLAR_DIRECT_SUM`` This activates the GPU implementation of the dipolar direct sum.
 
+-  ``DIPOLE_FIELD_TRACKING`` This enables the CPU implementation of the dipolar direct sum
+   to calculate the total dipole field at particle positions.
+
 -  ``ROTATION`` Switch on rotational degrees of freedom for the particles, as well as
    the corresponding quaternion integrator.
 

doc/sphinx/magnetostatics.rst

@@ -155,8 +155,13 @@ system's list of active actors. The only required parameter is the prefactor
 The CPU implementation has an optional argument ``n_replicas`` which
 adds periodic copies to the system along periodic directions. In that
 case, the minimum image convention is no longer used.
+Additionally, enabling the ``DIPOLE_FIELDS_TRACKING`` feature enables the CPU
+implementation to calculate the total dipole field at the position of each
+magnetic particle in the primary simulation box. These values are stored in
+the particle handle's ``dip_fld`` property and can be accessed directly or
+via an observable.
 
-Both implementations support MPI-parallelization.
+Both the CPU and GPU implementations support MPI-parallelization.
 
 
 .. _Barnes-Hut octree sum on GPU:

maintainer/configs/maxset.hpp

@@ -35,6 +35,7 @@
 #ifdef SCAFACOS
 #define SCAFACOS_DIPOLES
 #endif
+#define DIPOLE_FIELD_TRACKING
 
 #define ENGINE
 #define LB_ELECTROHYDRODYNAMICS

src/config/features.def

@@ -62,6 +62,7 @@ DIPOLES                         implies ROTATION
 DP3M                            equals DIPOLES and FFTW
 DIPOLAR_DIRECT_SUM              requires CUDA
 DIPOLAR_DIRECT_SUM              equals DIPOLES and ROTATION and CUDA
+DIPOLE_FIELD_TRACKING           implies DIPOLES
 DIPOLAR_BARNES_HUT              requires CUDA
 DIPOLAR_BARNES_HUT              equals DIPOLES and ROTATION and CUDA
 

src/core/Particle.hpp

@@ -151,6 +151,11 @@ struct ParticleProperties {
   double dipm = 0.;
 #endif
 
+#ifdef DIPOLE_FIELD_TRACKING
+  /** total dipole field */
+  Utils::Vector3d dip_fld = {0., 0., 0.};
+#endif
+
 #ifdef VIRTUAL_SITES_RELATIVE
   /** The following properties define, with respect to which real particle a
    *  virtual site is placed and at what distance. The relative orientation of
@@ -229,7 +234,9 @@ struct ParticleProperties {
 #ifdef DIPOLES
     ar &dipm;
 #endif
-
+#ifdef DIPOLE_FIELD_TRACKING
+    ar &dip_fld;
+#endif
 #ifdef VIRTUAL_SITES
     ar &is_virtual;
 #ifdef VIRTUAL_SITES_RELATIVE
@@ -497,6 +504,10 @@ struct Particle { // NOLINT(bugprone-exception-escape)
   auto &dipm() { return p.dipm; }
   auto calc_dip() const { return calc_director() * dipm(); }
 #endif
+#ifdef DIPOLE_FIELD_TRACKING
+  auto const &dip_fld() const { return p.dip_fld; }
+  auto &dip_fld() { return p.dip_fld; }
+#endif
 #ifdef ROTATIONAL_INERTIA
   auto const &rinertia() const { return p.rinertia; }
   auto &rinertia() { return p.rinertia; }

src/core/energy.cpp

@@ -40,6 +40,7 @@
 #include "magnetostatics/dipoles.hpp"
 
 #include <utils/Span.hpp>
+#include <utils/mpi/iall_gatherv.hpp>
 
 #include <memory>
 
@@ -154,3 +155,14 @@ double particle_short_range_energy_contribution(int pid) {
   }
   return ret;
 }
+
+#ifdef DIPOLE_FIELD_TRACKING
+void calc_long_range_fields() {
+  auto particles = cell_structure.local_particles();
+  Dipoles::calc_long_range_field(particles);
+}
+
+REGISTER_CALLBACK(calc_long_range_fields)
+
+void mpi_calc_long_range_fields() { mpi_call_all(calc_long_range_fields); }
+#endif

src/core/energy.hpp

@@ -51,3 +51,8 @@ double mpi_observable_compute_energy();
 double particle_short_range_energy_contribution(int pid);
 
 #endif
+#ifdef DIPOLE_FIELD_TRACKING
+/** Calculate dipole fields. */
+void calc_long_range_fields();
+void mpi_calc_long_range_fields();
+#endif

src/core/magnetostatics/dipolar_direct_sum.cpp

@@ -104,6 +104,25 @@ auto pair_potential(Utils::Vector3d const &d, Utils::Vector3d const &m1,
   return pe1 / r3 - 3.0 * pe2 * pe3 / r5;
 }
 
+/**
+ * @brief Dipole field contribution from a particle with dipole moment @c m1
+ * at a distance @c d.
+ *
+ * @param d Distance vector.
+ * @param m1 Dipole moment of one particle.
+ *
+ * @return Utils::Vector3d containing dipole field components.
+ */
+auto dipole_field(Utils::Vector3d const &d, Utils::Vector3d const &m1) {
+  auto const r2 = d * d;
+  auto const r = sqrt(r2);
+  auto const r3 = r2 * r;
+  auto const r5 = r3 * r2;
+  auto const pe2 = m1 * d;
+
+  return 3.0 * pe2 * d / r5 - m1 / r3;
+}
+
 /**
  * @brief Call kernel for every 3d index in a sphere around the origin.
  *
@@ -384,6 +403,45 @@ DipolarDirectSum::long_range_energy(ParticleRange const &particles) const {
   return prefactor * u;
 }
 
+/**
+ * @brief Calculate total dipole field of each particle.
+ *
+ * This employs a parallel N-square loop over all particles.
+ * Logically this is equivalent to the potential calculation
+ * in @ref DipolarDirectSum::long_range_energy, which calculates
+ * a naive N-square sum. The difference is summation range,
+ * and the kernel calculates the dipole field rather than the energy.
+ */
+#ifdef DIPOLE_FIELD_TRACKING
+void DipolarDirectSum::dipole_field_at_part(
+    ParticleRange const &particles) const {
+  auto const &box_l = ::box_geo.length();
+  /* collect particle data */
+  auto [local_particles, all_posmom, reqs, offset] =
+      gather_particle_data(particles, n_replicas);
+
+  auto const ncut = get_n_cut(n_replicas);
+  auto const with_replicas = (ncut.norm2() > 0);
+
+  boost::mpi::wait_all(reqs.begin(), reqs.end());
+
+  auto const local_posmom_begin = all_posmom.begin() + offset;
+  auto const local_posmom_end =
+      local_posmom_begin + static_cast<long>(local_particles.size());
+
+  Utils::Vector3d u_init = {0., 0., 0.};
+  auto p = local_particles.begin();
+  for (auto pi = local_posmom_begin; pi != local_posmom_end; ++pi, ++p) {
+    auto const u = image_sum(
+        all_posmom.begin(), all_posmom.end(), pi, with_replicas, ncut, box_l,
+        u_init, [](Utils::Vector3d const &rn, Utils::Vector3d const &mj) {
+          return dipole_field(rn, mj);
+        });
+    (*p)->dip_fld() = prefactor * u;
+  }
+}
+#endif
+
 DipolarDirectSum::DipolarDirectSum(double prefactor, int n_replicas)
     : prefactor{prefactor}, n_replicas{n_replicas} {
   if (prefactor <= 0.) {

src/core/magnetostatics/dipolar_direct_sum.hpp

@@ -49,6 +49,9 @@ struct DipolarDirectSum {
 
   double long_range_energy(ParticleRange const &particles) const;
   void add_long_range_forces(ParticleRange const &particles) const;
+#ifdef DIPOLE_FIELD_TRACKING
+  void dipole_field_at_part(ParticleRange const &particles) const;
+#endif
 };
 
 #endif // DIPOLES

src/core/magnetostatics/dipoles.cpp

@@ -186,6 +186,25 @@ struct LongRangeEnergy : public boost::static_visitor<double> {
 #endif
 };
 
+#ifdef DIPOLE_FIELD_TRACKING
+struct LongRangeField : public boost::static_visitor<void> {
+  ParticleRange const &m_particles;
+  explicit LongRangeField(ParticleRange const &particles)
+      : m_particles(particles) {}
+
+  void operator()(std::shared_ptr<DipolarDirectSum> const &actor) const {
+    actor->dipole_field_at_part(m_particles);
+  }
+
+  template <typename T,
+            std::enable_if_t<!traits::has_dipole_fields<T>::value> * = nullptr>
+  void operator()(std::shared_ptr<T> const &) const {
+    runtimeErrorMsg() << "Dipoles field calculation not implemented by "
+                      << "dipolar method " << Utils::demangle<T>();
+  }
+};
+#endif
+
 void calc_long_range_force(ParticleRange const &particles) {
   if (magnetostatics_actor) {
     boost::apply_visitor(LongRangeForce(particles), *magnetostatics_actor);
@@ -200,6 +219,14 @@ double calc_energy_long_range(ParticleRange const &particles) {
   return 0.;
 }
 
+#ifdef DIPOLE_FIELD_TRACKING
+void calc_long_range_field(ParticleRange const &particles) {
+  if (magnetostatics_actor) {
+    boost::apply_visitor(LongRangeField(particles), *magnetostatics_actor);
+  }
+}
+#endif
+
 namespace detail {
 bool flag_all_reduce(bool flag) {
   return boost::mpi::all_reduce(comm_cart, flag, std::logical_or<>());

src/core/magnetostatics/dipoles.hpp

@@ -88,6 +88,12 @@ namespace traits {
 template <typename T>
 using is_solver = std::is_convertible<std::shared_ptr<T>, MagnetostaticsActor>;
 
+/** @brief The dipolar method supports dipole fields calculation. */
+template <class T> struct has_dipole_fields : std::false_type {};
+#ifdef DIPOLE_FIELD_TRACKING
+template <> struct has_dipole_fields<DipolarDirectSum> : std::true_type {};
+#endif // DIPOLE_FIELD_TRACKING
+
 } // namespace traits
 
 void calc_pressure_long_range();
@@ -104,6 +110,9 @@ void on_cell_structure_change();
 
 void calc_long_range_force(ParticleRange const &particles);
 double calc_energy_long_range(ParticleRange const &particles);
+#ifdef DIPOLE_FIELD_TRACKING
+void calc_long_range_field(ParticleRange const &particles);
+#endif
 
 namespace detail {
 bool flag_all_reduce(bool flag);

src/core/observables/ParticleDipoleFields.hpp

@@ -0,0 +1,48 @@
+/*
+ * Copyright (C) 2023 The ESPResSo project
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo 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 General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#ifndef OBSERVABLES_PARTICLEDIPOLEFIELDS_HPP
+#define OBSERVABLES_PARTICLEDIPOLEFIELDS_HPP
+
+#include "PidObservable.hpp"
+#include "energy.hpp"
+
+#include <vector>
+
+namespace Observables {
+
+/** Extract particle dipole fields.
+ *  For \f$n\f$ particles, return \f$3 n\f$ dipole fields ordered as
+ *  \f$(h_d1_x, h_d1_y, h_d1_z, \dots, h_dn_x, h_dn_y, h_dn_z)\f$.
+ */
+class ParticleDipoleFields
+    : public ParticleObservable<ParticleObservables::DipoleFields> {
+public:
+  using ParticleObservable<
+      ParticleObservables::DipoleFields>::ParticleObservable;
+  std::vector<double>
+  evaluate(ParticleReferenceRange particles,
+           const ParticleObservables::traits<Particle> &traits) const override {
+    mpi_calc_long_range_fields();
+    return ParticleObservable<ParticleObservables::DipoleFields>::evaluate(
+        particles, traits);
+  }
+};
+
+} // namespace Observables
+#endif

src/core/observables/ParticleTraits.hpp

@@ -47,6 +47,13 @@ template <> struct traits<Particle> {
     return p.calc_dip();
 #else
     return Utils::Vector3d{};
+#endif
+  }
+  auto dipole_field(Particle const &p) const {
+#ifdef DIPOLE_FIELD_TRACKING
+    return p.dip_fld();
+#else
+    return Utils::Vector3d{};
 #endif
   }
   auto velocity_body(Particle const &p) const {

src/particle_observables/include/particle_observables/observable.hpp

@@ -54,6 +54,7 @@ using Forces = Map<Force>;
 using Positions = Map<Position>;
 using Velocities = Map<Velocity>;
 using Directors = Map<Director>;
+using DipoleFields = Map<DipoleField>;
 using BodyVelocities = Map<BodyVelocity>;
 using AngularVelocities = Map<AngularVelocity>;
 using BodyAngularVelocities = Map<BodyAngularVelocity>;

src/particle_observables/include/particle_observables/properties.hpp

@@ -118,6 +118,14 @@ struct DipoleMoment {
     return particle_traits.dipole_moment(p);
   }
 };
+
+struct DipoleField {
+  template <class Particle, class Traits = default_traits<Particle>>
+  decltype(auto) operator()(Particle const &p,
+                            Traits particle_traits = {}) const {
+    return particle_traits.dipole_field(p);
+  }
+};
 } // namespace ParticleObservables
 
 #endif // OBSERVABLES_PROPERTIES_HPP

src/python/espressomd/analyze.py

@@ -572,6 +572,13 @@ def particle_energy(self, particle):
         """
         return self.call_method("particle_energy", pid=particle.id)
 
+    def dipole_fields(self):
+        """
+        Calculate the total dipole field on each particle.
+        """
+        assert_features("DIPOLE_FIELD_TRACKING")
+        self.call_method("calc_long_range_fields")
+
     def dpd_stress(self):
         assert_features("DPD")
         return np.reshape(self.call_method("dpd_stress"), (3, 3))

src/python/espressomd/observables.py

@@ -509,6 +509,28 @@ class ParticleDirectors(Observable):
     _so_name = "Observables::ParticleDirectors"
 
 
+@script_interface_register
+class ParticleDipoleFields(Observable):
+
+    """Calculates the particle dipole fields for particles with given ids.
+
+    Output format: :math:`(h_d1_x,\\ h_d1_y,\\ h_d1_z),\\ (h_d2_x,\\ h_d2_y,\\ h_d2_z),\\ \\dots,\\ (h_dn_x,\\ h_dn_y,\\ h_dn_z)`.
+
+    The particles are ordered according to the list of ids passed to the observable.
+
+    Parameters
+    ----------
+    ids : array_like of :obj:`int`
+        The ids of (existing) particles to take into account.
+
+    Returns
+    -------
+    (N, 3) :obj:`ndarray` of :obj:`float`
+
+    """
+    _so_name = "Observables::ParticleDipoleFields"
+
+
 @script_interface_register
 class ParticleDistances(Observable):
 

src/python/espressomd/particle_data.py

@@ -245,6 +245,12 @@ class ParticleHandle(ScriptInterfaceHelper):
         .. note::
            This needs the feature ``DIPOLES``.
 
+     dip_fld: (3,) array_like of :obj:`float`
+        Total dipole field value at the position of the particle.
+
+        .. note::
+           This needs the feature ``DIPOLE_FIELD_TRACKING``.
+
     ext_force: (3,) array_like of :obj:`float`
         An additional external force applied to the particle.