Particle property ranges, remove ParticleRanges from P3M (#4737)

Partial fix for #4721, closes #4696

Description of changes:
- Implements the `ParticlePropertyRange` struct containing static methods to transform a `ParticleRange` to a `boost::iterator_range` of properties of the particles, such as `ForceRange` or `ParticleRange`, which can then be used in range-based for-loops.

The refactoring of the P3M part is not finished yet, but I think this PR could also be extended to remove the ParticleRanges from the remaining Coulomb interactions.

src/core/ParticlePropertyIterator.hpp

@@ -0,0 +1,64 @@
+/*
+ * 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/>.
+ */
+
+#pragma once
+
+#include "BoxGeometry.hpp"
+#include "Particle.hpp"
+#include "ParticleRange.hpp"
+
+#include <boost/iterator/transform_iterator.hpp>
+#include <boost/range/iterator_range.hpp>
+
+#include <utils/Vector.hpp>
+
+namespace ParticlePropertyRange {
+
+namespace detail {
+template <class Kernel>
+auto create_transform_range(ParticleRange const &particles, Kernel kernel) {
+  auto transform_iterator_begin =
+      boost::make_transform_iterator(particles.begin(), kernel);
+  auto transform_iterator_end =
+      boost::make_transform_iterator(particles.end(), kernel);
+  return boost::make_iterator_range<decltype(transform_iterator_begin)>(
+      transform_iterator_begin, transform_iterator_end);
+}
+} // namespace detail
+
+auto unfolded_pos_range(ParticleRange const &particles,
+                        BoxGeometry const &box) {
+  auto return_unfolded_pos = [&box](Particle &p) {
+    return ::detail::unfolded_position(p.pos(), p.image_box(), box.length());
+  };
+  return detail::create_transform_range(particles, return_unfolded_pos);
+}
+
+auto charge_range(ParticleRange const &particles) {
+  auto return_charge = [](Particle &p) -> double & { return p.q(); };
+  return detail::create_transform_range(particles, return_charge);
+}
+auto force_range(ParticleRange const &particles) {
+  auto return_force = [](Particle &p) -> Utils::Vector3d & {
+    return p.force();
+  };
+  return detail::create_transform_range(particles, return_force);
+}
+
+} // namespace ParticlePropertyRange

src/core/electrostatics/p3m.cpp

@@ -41,6 +41,7 @@
 #include "BoxGeometry.hpp"
 #include "LocalBox.hpp"
 #include "Particle.hpp"
+#include "ParticlePropertyIterator.hpp"
 #include "ParticleRange.hpp"
 #include "actor/visitors.hpp"
 #include "cell_system/CellStructure.hpp"
@@ -64,6 +65,7 @@
 #include <boost/mpi/collectives/broadcast.hpp>
 #include <boost/mpi/collectives/reduce.hpp>
 #include <boost/mpi/communicator.hpp>
+#include <boost/range/combine.hpp>
 #include <boost/range/numeric.hpp>
 
 #include <algorithm>
@@ -354,8 +356,9 @@ void CoulombP3M::assign_charge(double q, Utils::Vector3d const &real_pos) {
 }
 
 template <int cao> struct AssignForces {
+  template <typename combined_ranges>
   void operator()(p3m_data_struct &p3m, double force_prefac,
-                  ParticleRange const &particles) const {
+                  combined_ranges const &p_q_force_range) const {
     using Utils::make_const_span;
     using Utils::Span;
     using Utils::Vector;
@@ -365,9 +368,11 @@ template <int cao> struct AssignForces {
     /* charged particle counter */
     auto p_index = std::size_t{0ul};
 
-    for (auto &p : particles) {
-      if (p.q() != 0.0) {
-        auto const pref = p.q() * force_prefac;
+    for (auto zipped : p_q_force_range) {
+      auto p_q = boost::get<0>(zipped);
+      auto &p_force = boost::get<1>(zipped);
+      if (p_q != 0.0) {
+        auto const pref = p_q * force_prefac;
         auto const w = p3m.inter_weights.load<cao>(p_index);
 
         Utils::Vector3d force{};
@@ -376,22 +381,23 @@ template <int cao> struct AssignForces {
                                        p3m.E_mesh[2][ind]};
         });
 
-        p.force() -= pref * force;
+        p_force -= pref * force;
         ++p_index;
       }
     }
   }
 };
 
+template <typename combined_ranges>
 static auto calc_dipole_moment(boost::mpi::communicator const &comm,
-                               ParticleRange const &particles,
-                               BoxGeometry const &box_geo) {
-  auto const local_dip = boost::accumulate(
-      particles, Utils::Vector3d{},
-      [&box_geo](Utils::Vector3d const &dip, auto const &p) {
-        return dip + p.q() * box_geo.unfolded_position(p.pos(), p.image_box());
-      });
-
+                               combined_ranges const &p_q_unfolded_pos_range) {
+  auto const local_dip =
+      boost::accumulate(p_q_unfolded_pos_range, Utils::Vector3d{},
+                        [](Utils::Vector3d const &dip, auto const &q_pos) {
+                          auto const p_q = boost::get<0>(q_pos);
+                          auto const &p_unfolded_pos = boost::get<1>(q_pos);
+                          return dip + p_q * p_unfolded_pos;
+                        });
   return boost::mpi::all_reduce(comm, local_dip, std::plus<>());
 }
 
@@ -467,13 +473,19 @@ double CoulombP3M::long_range_kernel(bool force_flag, bool energy_flag,
   p3m.sm.gather_grid(p3m.rs_mesh.data(), comm_cart, p3m.local_mesh.dim);
   fft_perform_forw(p3m.rs_mesh.data(), p3m.fft, comm_cart);
 
+  auto p_q_range = ParticlePropertyRange::charge_range(particles);
+  auto p_force_range = ParticlePropertyRange::force_range(particles);
+  auto p_unfolded_pos_range =
+      ParticlePropertyRange::unfolded_pos_range(particles, box_geo);
+
   // Note: after these calls, the grids are in the order yzx and not xyz
   // anymore!!!
   /* The dipole moment is only needed if we don't have metallic boundaries. */
   auto const box_dipole =
       (p3m.params.epsilon != P3M_EPSILON_METALLIC)
-          ? calc_dipole_moment(comm_cart, particles, box_geo)
-          : Utils::Vector3d::broadcast(0.);
+          ? std::make_optional(calc_dipole_moment(
+                comm_cart, boost::combine(p_q_range, p_unfolded_pos_range)))
+          : std::nullopt;
   auto const volume = box_geo.volume();
   auto const pref = 4. * Utils::pi() / volume / (2. * p3m.params.epsilon + 1.);
 
@@ -520,14 +532,17 @@ double CoulombP3M::long_range_kernel(bool force_flag, bool energy_flag,
                        p3m.local_mesh.dim);
 
     auto const force_prefac = prefactor / volume;
-    Utils::integral_parameter<int, AssignForces, 1, 7>(p3m.params.cao, p3m,
-                                                       force_prefac, particles);
+    Utils::integral_parameter<int, AssignForces, 1, 7>(
+        p3m.params.cao, p3m, force_prefac,
+        boost::combine(p_q_range, p_force_range));
 
     // add dipole forces
-    if (p3m.params.epsilon != P3M_EPSILON_METALLIC) {
-      auto const dm = prefactor * pref * box_dipole;
-      for (auto &p : particles) {
-        p.force() -= p.q() * dm;
+    if (box_dipole) {
+      auto const dm = prefactor * pref * box_dipole.value();
+      for (auto zipped : boost::combine(p_q_range, p_force_range)) {
+        auto p_q = boost::get<0>(zipped);
+        auto &p_force = boost::get<1>(zipped);
+        p_force -= p_q * dm;
       }
     }
   }
@@ -551,8 +566,8 @@ double CoulombP3M::long_range_kernel(bool force_flag, bool energy_flag,
       energy -= p3m.square_sum_q * Utils::pi() /
                 (2. * volume * Utils::sqr(p3m.params.alpha));
       /* dipole correction */
-      if (p3m.params.epsilon != P3M_EPSILON_METALLIC) {
-        energy += pref * box_dipole.norm2();
+      if (box_dipole) {
+        energy += pref * box_dipole.value().norm2();
       }
     }
     return prefactor * energy;