diff --git a/src/core/ParticlePropertyIterator.hpp b/src/core/ParticlePropertyIterator.hpp new file mode 100644 index 0000000000..40a8cacf43 --- /dev/null +++ b/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 . + */ +#ifndef ESPRESSO_SRC_CORE_PARTICLE_EXTRACT_PROPERTIES_HPP +#define ESPRESSO_SRC_CORE_PARTICLE_EXTRACT_PROPERTIES_HPP + +#include "Particle.hpp" +#include "ParticleRange.hpp" +#include +#include +#include +#include + +struct ParticlePropertyRange { + + template + 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( + 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 diff --git a/src/core/electrostatics/p3m.cpp b/src/core/electrostatics/p3m.cpp index c954125b5e..805f4eaee7 100644 --- a/src/core/electrostatics/p3m.cpp +++ b/src/core/electrostatics/p3m.cpp @@ -49,6 +49,7 @@ #include "grid.hpp" #include "integrate.hpp" #include "tuning.hpp" +#include #include #include @@ -63,6 +64,7 @@ #include #include #include +#include #include #include @@ -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; } } }