diff --git a/src/core/ParticlePropertyIterator.hpp b/src/core/ParticlePropertyIterator.hpp index 2034be8ac3..e733eaaf01 100644 --- a/src/core/ParticlePropertyIterator.hpp +++ b/src/core/ParticlePropertyIterator.hpp @@ -42,19 +42,25 @@ auto create_transform_range(ParticleRange const &particles, Kernel kernel) { } } // namespace detail -auto unfolded_pos_range(ParticleRange const &particles, - BoxGeometry const &box) { +inline 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) { +inline auto pos_range(ParticleRange const &particles) { + auto return_pos = [](Particle &p) -> Utils::Vector3d & { return p.pos(); }; + return detail::create_transform_range(particles, return_pos); +} + +inline 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) { + +inline auto force_range(ParticleRange const &particles) { auto return_force = [](Particle &p) -> Utils::Vector3d & { return p.force(); }; diff --git a/src/core/electrostatics/elc.cpp b/src/core/electrostatics/elc.cpp index 77eb2e8e73..69be07a7eb 100644 --- a/src/core/electrostatics/elc.cpp +++ b/src/core/electrostatics/elc.cpp @@ -32,6 +32,7 @@ #include "BoxGeometry.hpp" #include "Particle.hpp" +#include "ParticlePropertyIterator.hpp" #include "ParticleRange.hpp" #include "cell_system/CellStructure.hpp" #include "communication.hpp" @@ -43,6 +44,7 @@ #include #include +#include #include #include @@ -1118,9 +1120,9 @@ static void p3m_assign_image_charge(elc_data const &elc, CoulombP3M &p3m, } } -template +template void charge_assign(elc_data const &elc, CoulombP3M &solver, - ParticleRange const &particles) { + combined_ranges const &p_q_pos_range) { if (protocol == ChargeProtocol::BOTH or protocol == ChargeProtocol::IMAGE) { solver.p3m.inter_weights.reset(solver.p3m.params.cao); } @@ -1128,49 +1130,52 @@ void charge_assign(elc_data const &elc, CoulombP3M &solver, for (int i = 0; i < solver.p3m.local_mesh.size; i++) solver.p3m.rs_mesh[i] = 0.; - for (auto const &p : particles) { - if (p.q() != 0.) { + for (auto zipped : p_q_pos_range) { + auto const p_q = boost::get<0>(zipped); + auto const &p_pos = boost::get<1>(zipped); + if (p_q != 0.) { if (protocol == ChargeProtocol::BOTH or protocol == ChargeProtocol::REAL) { - solver.assign_charge(p.q(), p.pos(), solver.p3m.inter_weights); + solver.assign_charge(p_q, p_pos, solver.p3m.inter_weights); } if (protocol == ChargeProtocol::BOTH or protocol == ChargeProtocol::IMAGE) { - p3m_assign_image_charge(elc, solver, p.q(), p.pos()); + p3m_assign_image_charge(elc, solver, p_q, p_pos); } } } } -template +template void modify_p3m_sums(elc_data const &elc, CoulombP3M &solver, - ParticleRange const &particles) { + combined_range const &p_q_pos_range) { Utils::Vector3d node_sums{}; - for (auto const &p : particles) { - auto const q = p.q(); - if (q != 0.) { - auto const z = p.pos()[2]; + for (auto zipped : p_q_pos_range) { + auto const p_q = boost::get<0>(zipped); + auto const &p_pos = boost::get<1>(zipped); + if (p_q != 0.) { + auto const p_z = p_pos[2]; if (protocol == ChargeProtocol::BOTH or protocol == ChargeProtocol::REAL) { node_sums[0] += 1.; - node_sums[1] += Utils::sqr(q); - node_sums[2] += q; + node_sums[1] += Utils::sqr(p_q); + node_sums[2] += p_q; } if (protocol == ChargeProtocol::BOTH or protocol == ChargeProtocol::IMAGE) { - if (z < elc.space_layer) { + if (p_z < elc.space_layer) { node_sums[0] += 1.; - node_sums[1] += Utils::sqr(elc.delta_mid_bot * q); - node_sums[2] += elc.delta_mid_bot * q; + node_sums[1] += Utils::sqr(elc.delta_mid_bot * p_q); + node_sums[2] += elc.delta_mid_bot * p_q; } - if (z > (elc.box_h - elc.space_layer)) { + if (p_z > (elc.box_h - elc.space_layer)) { node_sums[0] += 1.; - node_sums[1] += Utils::sqr(elc.delta_mid_top * q); - node_sums[2] += elc.delta_mid_top * q; + node_sums[1] += Utils::sqr(elc.delta_mid_top * p_q); + node_sums[2] += elc.delta_mid_top * p_q; } } } @@ -1190,6 +1195,10 @@ double ElectrostaticLayerCorrection::long_range_energy( auto &solver = *solver_ptr; auto const &box_geo = *get_system().box_geo; + auto p_q_range = ParticlePropertyRange::charge_range(particles); + auto p_pos_range = ParticlePropertyRange::pos_range(particles); + auto p_q_pos_range = boost::combine(p_q_range, p_pos_range); + // assign the original charges (they may not have been assigned yet) solver.charge_assign(particles); @@ -1203,17 +1212,17 @@ double ElectrostaticLayerCorrection::long_range_energy( 0.5 * elc.dielectric_layers_self_energy(solver, box_geo, particles); // assign both original and image charges - charge_assign(elc, solver, particles); - modify_p3m_sums(elc, solver, particles); + charge_assign(elc, solver, p_q_pos_range); + modify_p3m_sums(elc, solver, p_q_pos_range); energy += 0.5 * solver.long_range_energy(particles); // assign only the image charges now - charge_assign(elc, solver, particles); - modify_p3m_sums(elc, solver, particles); + charge_assign(elc, solver, p_q_pos_range); + modify_p3m_sums(elc, solver, p_q_pos_range); energy -= 0.5 * solver.long_range_energy(particles); // restore modified sums - modify_p3m_sums(elc, solver, particles); + modify_p3m_sums(elc, solver, p_q_pos_range); return energy; }, @@ -1226,17 +1235,20 @@ void ElectrostaticLayerCorrection::add_long_range_forces( std::visit( [this, &particles](auto const &solver_ptr) { auto &solver = *solver_ptr; + auto p_q_range = ParticlePropertyRange::charge_range(particles); + auto p_pos_range = ParticlePropertyRange::pos_range(particles); + auto p_q_pos_range = boost::combine(p_q_range, p_pos_range); if (elc.dielectric_contrast_on) { auto const &box_geo = *get_system().box_geo; - modify_p3m_sums(elc, solver, particles); - charge_assign(elc, solver, particles); + modify_p3m_sums(elc, solver, p_q_pos_range); + charge_assign(elc, solver, p_q_pos_range); elc.dielectric_layers_self_forces(solver, box_geo, particles); } else { solver.charge_assign(particles); } solver.add_long_range_forces(particles); if (elc.dielectric_contrast_on) { - modify_p3m_sums(elc, solver, particles); + modify_p3m_sums(elc, solver, p_q_pos_range); } }, base_solver); diff --git a/src/core/electrostatics/p3m.cpp b/src/core/electrostatics/p3m.cpp index 1a78b36d10..7bb20f4a27 100644 --- a/src/core/electrostatics/p3m.cpp +++ b/src/core/electrostatics/p3m.cpp @@ -323,10 +323,13 @@ template struct AssignCharge { [q, &p3m](int ind, double w) { p3m.rs_mesh[ind] += w * q; }); } - void operator()(p3m_data_struct &p3m, ParticleRange const &particles) { - for (auto &p : particles) { - if (p.q() != 0.0) { - this->operator()(p3m, p.q(), p.pos(), p3m.inter_weights); + template + void operator()(p3m_data_struct &p3m, combined_ranges const &p_q_pos_range) { + for (auto zipped : p_q_pos_range) { + auto const p_q = boost::get<0>(zipped); + auto const &p_pos = boost::get<1>(zipped); + if (p_q != 0.0) { + this->operator()(p3m, p_q, p_pos, p3m.inter_weights); } } } @@ -340,8 +343,11 @@ void CoulombP3M::charge_assign(ParticleRange const &particles) { for (int i = 0; i < p3m.local_mesh.size; i++) p3m.rs_mesh[i] = 0.0; - Utils::integral_parameter(p3m.params.cao, p3m, - particles); + auto p_q_range = ParticlePropertyRange::charge_range(particles); + auto p_pos_range = ParticlePropertyRange::pos_range(particles); + + Utils::integral_parameter( + p3m.params.cao, p3m, boost::combine(p_q_range, p_pos_range)); } void CoulombP3M::assign_charge(double q, Utils::Vector3d const &real_pos, diff --git a/src/core/electrostatics/p3m.hpp b/src/core/electrostatics/p3m.hpp index 776eb57cb3..a53f408d92 100644 --- a/src/core/electrostatics/p3m.hpp +++ b/src/core/electrostatics/p3m.hpp @@ -99,10 +99,6 @@ struct CoulombP3M : public Coulomb::Actor { bool is_tuned() const { return m_is_tuned; } - /** Compute the k-space part of forces and energies. */ - double kernel(bool force_flag, bool energy_flag, - ParticleRange const &particles); - /** @brief Recalculate all derived parameters. */ void init(); void on_activation() {