add ParticlePropertyRange struct, start p3m refactor

src/core/ParticlePropertyIterator.hpp

@@ -0,0 +1,58 @@
+/*
+ * Copyright (C) 2010-2022 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 ESPRESSO_SRC_CORE_PARTICLE_EXTRACT_PROPERTIES_HPP
+#define ESPRESSO_SRC_CORE_PARTICLE_EXTRACT_PROPERTIES_HPP
+
+#include "Particle.hpp"
+#include "ParticleRange.hpp"
+#include <boost/iterator/transform_iterator.hpp>
+#include <boost/range/iterator_range.hpp>
+#include <tuple>
+#include <utils/Vector.hpp>
+
+struct ParticlePropertyRange {
+
+  template <class Kernel>
+  static auto create_transform_range(const ParticleRange &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);
+  }
+
+  static auto pos_range(const ParticleRange &particles) {
+    auto return_pos = [](Particle &p) { return p.pos(); };
+    return create_transform_range(particles, return_pos);
+  }
+
+  static auto charge_range(const ParticleRange &particles) {
+    auto return_charge = [](Particle &p) { return p.q(); };
+    return create_transform_range(particles, return_charge);
+  }
+
+  static auto force_range(const ParticleRange &particles) {
+    auto return_force = [](Particle &p) { return p.force(); };
+    return create_transform_range(particles, return_force);
+  }
+};
+
+#endif

src/core/electrostatics/p3m.cpp

@@ -49,6 +49,7 @@
 #include "grid.hpp"
 #include "integrate.hpp"
 #include "tuning.hpp"
+#include <ParticlePropertyIterator.hpp>
 
 #include <utils/Span.hpp>
 #include <utils/Vector.hpp>
@@ -63,6 +64,7 @@
 #include <boost/mpi/collectives/reduce.hpp>
 #include <boost/mpi/communicator.hpp>
 #include <boost/optional.hpp>
+#include <boost/range/combine.hpp>
 #include <boost/range/numeric.hpp>
 
 #include <algorithm>
@@ -455,6 +457,9 @@ 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 particle_charge_range = ParticlePropertyRange::charge_range(particles);
+  auto particle_force_range = ParticlePropertyRange::force_range(particles);
+
   // 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. */
@@ -514,8 +519,12 @@ double CoulombP3M::long_range_kernel(bool force_flag, bool energy_flag,
     // add dipole forces
     if (p3m.params.epsilon != P3M_EPSILON_METALLIC) {
       auto const dm = prefactor * pref * box_dipole.value();
-      for (auto &p : particles) {
-        p.force() -= p.q() * dm;
+      Utils::Vector3d p_force;
+      double p_q;
+      for (auto zipped :
+           boost::combine(particle_force_range, particle_charge_range)) {
+        boost::tie(p_force, p_q) = zipped;
+        p_force -= p_q * dm;
       }
     }
   }