diff --git a/src/core/cells.cpp b/src/core/cells.cpp index d4a3d1e8590..3b1a1c269a5 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -418,7 +418,7 @@ void cells_resort_particles(int global_flag) { realloc_particlelist(&displaced_parts, 0); - on_resort_particles(); + on_resort_particles(local_cells.particles()); CELL_TRACE( fprintf(stderr, "%d: leaving cells_resort_particles\n", this_node)); diff --git a/src/core/electrostatics_magnetostatics/coulomb.cpp b/src/core/electrostatics_magnetostatics/coulomb.cpp index 59e81c6b0c6..bf71440af35 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.cpp +++ b/src/core/electrostatics_magnetostatics/coulomb.cpp @@ -41,7 +41,8 @@ void pressure_n(int &n_coulomb) { } void calc_pressure_long_range(Observable_stat &virials, - Observable_stat &p_tensor) { + Observable_stat &p_tensor, + const ParticleRange &particles) { switch (coulomb.method) { #ifdef P3M case COULOMB_ELC_P3M: @@ -54,9 +55,9 @@ void calc_pressure_long_range(Observable_stat &virials, "WARNING: pressure calculated, but GPU P3M pressure not implemented\n"); break; case COULOMB_P3M: { - p3m_charge_assign(); - virials.coulomb[1] = p3m_calc_kspace_forces(0, 1); - p3m_charge_assign(); + p3m_charge_assign(particles); + virials.coulomb[1] = p3m_calc_kspace_forces(0, 1, particles); + p3m_charge_assign(particles); p3m_calc_kspace_stress(p_tensor.coulomb + 9); break; } @@ -220,7 +221,7 @@ void on_coulomb_change() { } } -void on_resort_particles() { +void on_resort_particles(const ParticleRange &particles) { switch (coulomb.method) { #ifdef P3M case COULOMB_ELC_P3M: @@ -228,7 +229,7 @@ void on_resort_particles() { break; #endif case COULOMB_MMM2D: - MMM2D_on_resort_particles(); + MMM2D_on_resort_particles(particles); break; default: break; @@ -287,23 +288,23 @@ void init() { } } -void calc_long_range_force() { +void calc_long_range_force(const ParticleRange &particles) { switch (coulomb.method) { #ifdef P3M case COULOMB_ELC_P3M: if (elc_params.dielectric_contrast_on) { - ELC_P3M_modify_p3m_sums_both(); - ELC_p3m_charge_assign_both(); - ELC_P3M_self_forces(); + ELC_P3M_modify_p3m_sums_both(particles); + ELC_p3m_charge_assign_both(particles); + ELC_P3M_self_forces(particles); } else - p3m_charge_assign(); + p3m_charge_assign(particles); - p3m_calc_kspace_forces(1, 0); + p3m_calc_kspace_forces(1, 0, particles); if (elc_params.dielectric_contrast_on) - ELC_P3M_restore_p3m_sums(); + ELC_P3M_restore_p3m_sums(particles); - ELC_add_force(); + ELC_add_force(particles); break; #endif @@ -320,17 +321,17 @@ void calc_long_range_force() { #endif #ifdef P3M case COULOMB_P3M: - p3m_charge_assign(); + p3m_charge_assign(particles); #ifdef NPT if (integ_switch == INTEG_METHOD_NPT_ISO) - nptiso.p_vir[0] += p3m_calc_kspace_forces(1, 1); + nptiso.p_vir[0] += p3m_calc_kspace_forces(1, 1, particles); else #endif - p3m_calc_kspace_forces(1, 0); + p3m_calc_kspace_forces(1, 0, particles); break; #endif case COULOMB_MMM2D: - MMM2D_add_far_force(); + MMM2D_add_far_force(particles); MMM2D_dielectric_layers_force_contribution(); break; #ifdef SCAFACOS @@ -351,7 +352,8 @@ void calc_long_range_force() { #endif } -void calc_energy_long_range(Observable_stat &energy) { +void calc_energy_long_range(Observable_stat &energy, + const ParticleRange &particles) { switch (coulomb.method) { #ifdef P3M case COULOMB_P3M_GPU: @@ -359,35 +361,36 @@ void calc_energy_long_range(Observable_stat &energy) { << "long range energy calculation not implemented for GPU P3M"; break; case COULOMB_P3M: - p3m_charge_assign(); - energy.coulomb[1] = p3m_calc_kspace_forces(0, 1); + p3m_charge_assign(particles); + energy.coulomb[1] = p3m_calc_kspace_forces(0, 1, particles); break; case COULOMB_ELC_P3M: // assign the original charges first // they may not have been assigned yet - p3m_charge_assign(); + p3m_charge_assign(particles); if (!elc_params.dielectric_contrast_on) - energy.coulomb[1] = p3m_calc_kspace_forces(0, 1); + energy.coulomb[1] = p3m_calc_kspace_forces(0, 1, particles); else { - energy.coulomb[1] = 0.5 * p3m_calc_kspace_forces(0, 1); - energy.coulomb[1] += 0.5 * ELC_P3M_dielectric_layers_energy_self(); + energy.coulomb[1] = 0.5 * p3m_calc_kspace_forces(0, 1, particles); + energy.coulomb[1] += + 0.5 * ELC_P3M_dielectric_layers_energy_self(particles); // assign both original and image charges now - ELC_p3m_charge_assign_both(); - ELC_P3M_modify_p3m_sums_both(); + ELC_p3m_charge_assign_both(particles); + ELC_P3M_modify_p3m_sums_both(particles); - energy.coulomb[1] += 0.5 * p3m_calc_kspace_forces(0, 1); + energy.coulomb[1] += 0.5 * p3m_calc_kspace_forces(0, 1, particles); // assign only the image charges now - ELC_p3m_charge_assign_image(); - ELC_P3M_modify_p3m_sums_image(); + ELC_p3m_charge_assign_image(particles); + ELC_P3M_modify_p3m_sums_image(particles); - energy.coulomb[1] -= 0.5 * p3m_calc_kspace_forces(0, 1); + energy.coulomb[1] -= 0.5 * p3m_calc_kspace_forces(0, 1, particles); // restore modified sums - ELC_P3M_restore_p3m_sums(); + ELC_P3M_restore_p3m_sums(particles); } - energy.coulomb[2] = ELC_energy(); + energy.coulomb[2] = ELC_energy(particles); break; #endif #ifdef SCAFACOS @@ -397,7 +400,7 @@ void calc_energy_long_range(Observable_stat &energy) { break; #endif case COULOMB_MMM2D: - *energy.coulomb += MMM2D_far_energy(); + *energy.coulomb += MMM2D_far_energy(particles); *energy.coulomb += MMM2D_dielectric_layers_energy_contribution(); break; default: diff --git a/src/core/electrostatics_magnetostatics/coulomb.hpp b/src/core/electrostatics_magnetostatics/coulomb.hpp index 8fa5ccb9795..6166573fc99 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.hpp +++ b/src/core/electrostatics_magnetostatics/coulomb.hpp @@ -6,6 +6,7 @@ #ifdef ELECTROSTATICS #include "Observable_stat.hpp" +#include #include /** \name Type codes for the type of Coulomb interaction @@ -51,7 +52,8 @@ extern Coulomb_parameters coulomb; namespace Coulomb { void pressure_n(int &n_coulomb); void calc_pressure_long_range(Observable_stat &virials, - Observable_stat &p_tensor); + Observable_stat &p_tensor, + const ParticleRange &particles); void sanity_checks(int &state); double cutoff(const Utils::Vector3d &box_l); @@ -60,13 +62,14 @@ void deactivate(); void integrate_sanity_check(); void on_observable_calc(); void on_coulomb_change(); -void on_resort_particles(); +void on_resort_particles(const ParticleRange &particles); void on_boxl_change(); void init(); -void calc_long_range_force(); +void calc_long_range_force(const ParticleRange &particles); -void calc_energy_long_range(Observable_stat &energy); +void calc_energy_long_range(Observable_stat &energy, + const ParticleRange &particles); int energy_n(); int iccp3m_sanity_check(); diff --git a/src/core/electrostatics_magnetostatics/dipole.cpp b/src/core/electrostatics_magnetostatics/dipole.cpp index f802536ce57..11d59b2937d 100644 --- a/src/core/electrostatics_magnetostatics/dipole.cpp +++ b/src/core/electrostatics_magnetostatics/dipole.cpp @@ -155,34 +155,34 @@ void init() { } } -void calc_long_range_force() { +void calc_long_range_force(const ParticleRange &particles) { switch (dipole.method) { #ifdef DP3M case DIPOLAR_MDLC_P3M: - add_mdlc_force_corrections(); + add_mdlc_force_corrections(particles); // fall through case DIPOLAR_P3M: - dp3m_dipole_assign(); + dp3m_dipole_assign(particles); #ifdef NPT if (integ_switch == INTEG_METHOD_NPT_ISO) { - nptiso.p_vir[0] += dp3m_calc_kspace_forces(1, 1); + nptiso.p_vir[0] += dp3m_calc_kspace_forces(1, 1, particles); fprintf(stderr, "dipolar_P3M at this moment is added to p_vir[0]\n"); } else #endif - dp3m_calc_kspace_forces(1, 0); + dp3m_calc_kspace_forces(1, 0, particles); break; #endif case DIPOLAR_ALL_WITH_ALL_AND_NO_REPLICA: - dawaanr_calculations(1, 0); + dawaanr_calculations(1, 0, particles); break; #ifdef DP3M case DIPOLAR_MDLC_DS: - add_mdlc_force_corrections(); + add_mdlc_force_corrections(particles); // fall through #endif case DIPOLAR_DS: - magnetic_dipolar_direct_sum_calculations(1, 0); + magnetic_dipolar_direct_sum_calculations(1, 0, particles); break; case DIPOLAR_DS_GPU: // Do nothing. It's an actor @@ -205,30 +205,33 @@ void calc_long_range_force() { } } -void calc_energy_long_range(Observable_stat &energy) { +void calc_energy_long_range(Observable_stat &energy, + const ParticleRange &particles) { switch (dipole.method) { #ifdef DP3M case DIPOLAR_P3M: - dp3m_dipole_assign(); - energy.dipolar[1] = dp3m_calc_kspace_forces(0, 1); + dp3m_dipole_assign(particles); + energy.dipolar[1] = dp3m_calc_kspace_forces(0, 1, particles); break; case DIPOLAR_MDLC_P3M: - dp3m_dipole_assign(); - energy.dipolar[1] = dp3m_calc_kspace_forces(0, 1); - energy.dipolar[2] = add_mdlc_energy_corrections(); + dp3m_dipole_assign(particles); + energy.dipolar[1] = dp3m_calc_kspace_forces(0, 1, particles); + energy.dipolar[2] = add_mdlc_energy_corrections(particles); break; #endif case DIPOLAR_ALL_WITH_ALL_AND_NO_REPLICA: - energy.dipolar[1] = dawaanr_calculations(0, 1); + energy.dipolar[1] = dawaanr_calculations(0, 1, particles); break; #ifdef DP3M case DIPOLAR_MDLC_DS: - energy.dipolar[1] = magnetic_dipolar_direct_sum_calculations(0, 1); - energy.dipolar[2] = add_mdlc_energy_corrections(); + energy.dipolar[1] = + magnetic_dipolar_direct_sum_calculations(0, 1, particles); + energy.dipolar[2] = add_mdlc_energy_corrections(particles); break; #endif case DIPOLAR_DS: - energy.dipolar[1] = magnetic_dipolar_direct_sum_calculations(0, 1); + energy.dipolar[1] = + magnetic_dipolar_direct_sum_calculations(0, 1, particles); break; case DIPOLAR_DS_GPU: break; diff --git a/src/core/electrostatics_magnetostatics/dipole.hpp b/src/core/electrostatics_magnetostatics/dipole.hpp index d21fb88f9c7..4883919fb9b 100644 --- a/src/core/electrostatics_magnetostatics/dipole.hpp +++ b/src/core/electrostatics_magnetostatics/dipole.hpp @@ -10,6 +10,7 @@ extern double dipolar_cutoff; #include +#include #include /** \name Compounds for Dipole interactions */ @@ -69,9 +70,10 @@ void on_coulomb_change(); void on_boxl_change(); void init(); -void calc_long_range_force(); +void calc_long_range_force(const ParticleRange &particles); -void calc_energy_long_range(Observable_stat &energy); +void calc_energy_long_range(Observable_stat &energy, + const ParticleRange &particles); void energy_n(int &n_dipolar); int set_mesh(); diff --git a/src/core/electrostatics_magnetostatics/elc.cpp b/src/core/electrostatics_magnetostatics/elc.cpp index 1bc8fd3e62a..b5086eadf8f 100644 --- a/src/core/electrostatics_magnetostatics/elc.cpp +++ b/src/core/electrostatics_magnetostatics/elc.cpp @@ -112,8 +112,8 @@ static int n_scycache; /** \name sin/cos storage */ /*@{*/ -static void prepare_scx_cache(); -static void prepare_scy_cache(); +static void prepare_scx_cache(const ParticleRange &particles); +static void prepare_scy_cache(const ParticleRange &particles); /*@}*/ /** \name common code */ /*@{*/ @@ -121,25 +121,27 @@ static void distribute(int size); /*@}*/ /** \name p=0 per frequency code */ /*@{*/ -static void setup_P(int p, double omega); -static void add_P_force(); +static void setup_P(int p, double omega, const ParticleRange &particles); +static void add_P_force(const ParticleRange &particles); static double P_energy(double omega); /*@}*/ /** \name q=0 per frequency code */ /*@{*/ -static void setup_Q(int q, double omega); -static void add_Q_force(); +static void setup_Q(int q, double omega, const ParticleRange &particles); +static void add_Q_force(const ParticleRange &particles); static double Q_energy(double omega); /*@}*/ /** \name p,q <> 0 per frequency code */ /*@{*/ -static void setup_PQ(int p, int q, double omega); -static void add_PQ_force(int p, int q, double omega); +static void setup_PQ(int p, int q, double omega, + const ParticleRange &particles); +static void add_PQ_force(int p, int q, double omega, + const ParticleRange &particles); static double PQ_energy(double omega); -static void add_dipole_force(); -static double dipole_energy(); -static double z_energy(); -static void add_z_force(); +static void add_dipole_force(const ParticleRange &particles); +static double dipole_energy(const ParticleRange &particles); +static double z_energy(const ParticleRange &particles); +static void add_z_force(const ParticleRange &particles); /*@}*/ /* COMMON */ @@ -157,7 +159,7 @@ void ELC_setup_constants() { /* SC Cache */ /************/ -static void prepare_scx_cache() { +static void prepare_scx_cache(const ParticleRange &particles) { int freq; double arg; @@ -165,7 +167,7 @@ static void prepare_scx_cache() { double pref = C_2PI * ux * freq; int o = (freq - 1) * n_localpart; int ic = 0; - for (auto const &part : local_cells.particles()) { + for (auto const &part : particles) { arg = pref * part.r.p[0]; scxcache[o + ic].s = sin(arg); scxcache[o + ic].c = cos(arg); @@ -174,7 +176,7 @@ static void prepare_scx_cache() { } } -static void prepare_scy_cache() { +static void prepare_scy_cache(const ParticleRange &particles) { int freq; double arg; @@ -182,7 +184,7 @@ static void prepare_scy_cache() { double pref = C_2PI * uy * freq; int o = (freq - 1) * n_localpart; int ic = 0; - for (auto const &part : local_cells.particles()) { + for (auto const &part : particles) { arg = pref * part.r.p[1]; scycache[o + ic].s = sin(arg); scycache[o + ic].c = cos(arg); @@ -269,9 +271,9 @@ static void checkpoint(char *text, int p, int q, int e_size) { #endif #ifdef LOG_FORCES -static void clear_log_forces(char *where) { +static void clear_log_forces(char *where, const ParticleRange &particles) { fprintf(stderr, "%s\n", where); - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { fprintf(stderr, "%d %g %g %g\n", p.p.identity, p.f.f[0], p.f.f[1], p.f.f[2]); for (int j = 0; j < 3; j++) @@ -279,17 +281,19 @@ static void clear_log_forces(char *where) { } } #else -#define clear_log_forces(w) +#define clear_log_forces(w, particles) #endif /*****************************************************************/ /* dipole terms */ /*****************************************************************/ -static void add_dipole_force() { +static void add_dipole_force(const ParticleRange &particles) { double pref = coulomb.prefactor * 4 * M_PI * ux * uy * uz; int size = 3; + auto local_particles = particles; + /* for nonneutral systems, this shift gives the background contribution (rsp. for this shift, the DM of the background is zero) */ double shift = 0.5 * box_geo.length()[2]; @@ -301,7 +305,7 @@ static void add_dipole_force() { gblcblk[1] = 0; // sum q_i z_i gblcblk[2] = 0; // sum q_i - for (auto const &p : local_cells.particles()) { + for (auto const &p : local_particles) { gblcblk[0] += p.p.q * (p.r.p[2] - shift); gblcblk[1] += p.p.q * p.r.p[2]; gblcblk[2] += p.p.q; @@ -335,7 +339,7 @@ static void add_dipole_force() { field_tot -= coulomb.field_applied + coulomb.field_induced; } - for (auto &p : local_cells.particles()) { + for (auto &p : local_particles) { p.f.f[2] -= field_tot * p.p.q; if (!elc_params.neutralize) { @@ -345,7 +349,7 @@ static void add_dipole_force() { } } -static double dipole_energy() { +static double dipole_energy(const ParticleRange &particles) { double pref = coulomb.prefactor * 2 * M_PI * ux * uy * uz; int size = 7; /* for nonneutral systems, this shift gives the background contribution @@ -362,7 +366,7 @@ static double dipole_energy() { gblcblk[5] = 0; // sum q_i (z_i - L/2)^2 boundary layers gblcblk[6] = 0; // sum q_i z_i primary box - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { gblcblk[0] += p.p.q; gblcblk[2] += p.p.q * (p.r.p[2] - shift); gblcblk[4] += p.p.q * (Utils::sqr(p.r.p[2] - shift)); @@ -437,7 +441,7 @@ inline double image_sum_t(double q, double z) { } /*****************************************************************/ -static double z_energy() { +static double z_energy(const ParticleRange &particles) { double pref = coulomb.prefactor * 2 * M_PI * ux * uy; int size = 4; @@ -449,7 +453,7 @@ static double z_energy() { if (elc_params.dielectric_contrast_on) { if (elc_params.const_pot) { clear_vec(gblcblk, size); - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { gblcblk[0] += p.p.q; gblcblk[1] += p.p.q * (p.r.p[2] - shift); if (p.r.p[2] < elc_params.space_layer) { @@ -469,7 +473,7 @@ static double z_energy() { double fac_delta = delta / (1 - delta); clear_vec(gblcblk, size); - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { gblcblk[0] += p.p.q; gblcblk[1] += p.p.q * (p.r.p[2] - shift); if (elc_params.dielectric_contrast_on) { @@ -515,15 +519,16 @@ static double z_energy() { } /*****************************************************************/ -static void add_z_force() { +static void add_z_force(const ParticleRange &particles) { double pref = coulomb.prefactor * 2 * M_PI * ux * uy; if (elc_params.dielectric_contrast_on) { + auto local_particles = particles; int size = 1; if (elc_params.const_pot) { clear_vec(gblcblk, size); /* just counter the 2 pi |z| contribution stemming from P3M */ - for (auto &p : local_cells.particles()) { + for (auto &p : local_particles) { if (p.r.p[2] < elc_params.space_layer) gblcblk[0] -= elc_params.delta_mid_bot * p.p.q; if (p.r.p[2] > (elc_params.h - elc_params.space_layer)) @@ -536,7 +541,7 @@ static void add_z_force() { double fac_delta = delta / (1 - delta); clear_vec(gblcblk, size); - for (auto &p : local_cells.particles()) { + for (auto &p : local_particles) { if (p.r.p[2] < elc_params.space_layer) { gblcblk[0] += fac_delta * (elc_params.delta_mid_bot + 1) * p.p.q; } else { @@ -559,7 +564,7 @@ static void add_z_force() { distribute(size); - for (auto &p : local_cells.particles()) { + for (auto &p : local_particles) { p.f.f[2] += gblcblk[0] * p.p.q; } } @@ -569,7 +574,7 @@ static void add_z_force() { /* PoQ exp sum */ /*****************************************************************/ -static void setup_P(int p, double omega) { +static void setup_P(int p, double omega, const ParticleRange &particles) { int ic, o = (p - 1) * n_localpart; double pref = -coulomb.prefactor * 4 * M_PI * ux * uy / (expm1(omega * box_geo.length()[2])); @@ -592,7 +597,7 @@ static void setup_P(int p, double omega) { clear_vec(gblcblk, size); ic = 0; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { double e = exp(omega * p.r.p[2]); partblk[size * ic + POQESM] = p.p.q * scxcache[o + ic].s / e; @@ -676,7 +681,7 @@ static void setup_P(int p, double omega) { } } -static void setup_Q(int q, double omega) { +static void setup_Q(int q, double omega, const ParticleRange &particles) { int ic, o = (q - 1) * n_localpart; double pref = -coulomb.prefactor * 4 * M_PI * ux * uy / (expm1(omega * box_geo.length()[2])); @@ -698,7 +703,7 @@ static void setup_Q(int q, double omega) { clear_vec(lclimge, size); clear_vec(gblcblk, size); ic = 0; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { double e = exp(omega * p.r.p[2]); partblk[size * ic + POQESM] = p.p.q * scycache[o + ic].s / e; @@ -783,12 +788,12 @@ static void setup_Q(int q, double omega) { } } -static void add_P_force() { +static void add_P_force(const ParticleRange &particles) { int ic; int size = 4; ic = 0; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { p.f.f[0] += partblk[size * ic + POQESM] * gblcblk[POQECP] - partblk[size * ic + POQECM] * gblcblk[POQESP] + partblk[size * ic + POQESP] * gblcblk[POQECM] - @@ -816,12 +821,12 @@ static double P_energy(double omega) { return eng; } -static void add_Q_force() { +static void add_Q_force(const ParticleRange &particles) { int ic; int size = 4; ic = 0; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { p.f.f[1] += partblk[size * ic + POQESM] * gblcblk[POQECP] - partblk[size * ic + POQECM] * gblcblk[POQESP] + partblk[size * ic + POQESP] * gblcblk[POQECM] - @@ -852,7 +857,8 @@ static double Q_energy(double omega) { /* PQ particle blocks */ /*****************************************************************/ -static void setup_PQ(int p, int q, double omega) { +static void setup_PQ(int p, int q, double omega, + const ParticleRange &particles) { int ic, ox = (p - 1) * n_localpart, oy = (q - 1) * n_localpart; double pref = -coulomb.prefactor * 8 * M_PI * ux * uy / (expm1(omega * box_geo.length()[2])); @@ -874,7 +880,7 @@ static void setup_PQ(int p, int q, double omega) { clear_vec(gblcblk, size); ic = 0; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { double e = exp(omega * p.r.p[2]); partblk[size * ic + PQESSM] = @@ -981,7 +987,8 @@ static void setup_PQ(int p, int q, double omega) { } } -static void add_PQ_force(int p, int q, double omega) { +static void add_PQ_force(int p, int q, double omega, + const ParticleRange &particles) { int ic; double pref_x = C_2PI * ux * p / omega; @@ -989,7 +996,7 @@ static void add_PQ_force(int p, int q, double omega) { int size = 8; ic = 0; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { p.f.f[0] += pref_x * (partblk[size * ic + PQESCM] * gblcblk[PQECCP] + partblk[size * ic + PQESSM] * gblcblk[PQECSP] - partblk[size * ic + PQECCM] * gblcblk[PQESCP] - @@ -1040,37 +1047,37 @@ static double PQ_energy(double omega) { /* main loops */ /*****************************************************************/ -void ELC_add_force() { +void ELC_add_force(const ParticleRange &particles) { int p, q; double omega; - prepare_scx_cache(); - prepare_scy_cache(); + prepare_scx_cache(particles); + prepare_scy_cache(particles); - clear_log_forces("start"); + clear_log_forces("start", particles); - add_dipole_force(); + add_dipole_force(particles); - clear_log_forces("dipole"); + clear_log_forces("dipole", particles); - add_z_force(); + add_z_force(particles); - clear_log_forces("z_force"); + clear_log_forces("z_force", particles); /* the second condition is just for the case of numerical accident */ for (p = 1; ux * (p - 1) < elc_params.far_cut && p <= n_scxcache; p++) { omega = C_2PI * ux * p; - setup_P(p, omega); + setup_P(p, omega, particles); distribute(4); - add_P_force(); + add_P_force(particles); checkpoint("************distri p", p, 0, 2); } for (q = 1; uy * (q - 1) < elc_params.far_cut && q <= n_scycache; q++) { omega = C_2PI * uy * q; - setup_Q(q, omega); + setup_Q(q, omega, particles); distribute(4); - add_Q_force(); + add_Q_force(particles); checkpoint("************distri q", 0, q, 2); } @@ -1080,37 +1087,37 @@ void ELC_add_force() { q <= n_scycache; q++) { omega = C_2PI * sqrt(Utils::sqr(ux * p) + Utils::sqr(uy * q)); - setup_PQ(p, q, omega); + setup_PQ(p, q, omega, particles); distribute(8); - add_PQ_force(p, q, omega); + add_PQ_force(p, q, omega, particles); checkpoint("************distri pq", p, q, 4); } } - clear_log_forces("end"); + clear_log_forces("end", particles); } -double ELC_energy() { +double ELC_energy(const ParticleRange &particles) { double eng; int p, q; double omega; - eng = dipole_energy(); - eng += z_energy(); - prepare_scx_cache(); - prepare_scy_cache(); + eng = dipole_energy(particles); + eng += z_energy(particles); + prepare_scx_cache(particles); + prepare_scy_cache(particles); /* the second condition is just for the case of numerical accident */ for (p = 1; ux * (p - 1) < elc_params.far_cut && p <= n_scxcache; p++) { omega = C_2PI * ux * p; - setup_P(p, omega); + setup_P(p, omega, particles); distribute(4); eng += P_energy(omega); checkpoint("E************distri p", p, 0, 2); } for (q = 1; uy * (q - 1) < elc_params.far_cut && q <= n_scycache; q++) { omega = C_2PI * uy * q; - setup_Q(q, omega); + setup_Q(q, omega, particles); distribute(4); eng += Q_energy(omega); checkpoint("E************distri q", 0, q, 2); @@ -1121,7 +1128,7 @@ double ELC_energy() { q <= n_scycache; q++) { omega = C_2PI * sqrt(Utils::sqr(ux * p) + Utils::sqr(uy * q)); - setup_PQ(p, q, omega); + setup_PQ(p, q, omega, particles); distribute(8); eng += PQ_energy(omega); checkpoint("E************distri pq", p, q, 4); @@ -1316,11 +1323,11 @@ int ELC_set_params(double maxPWerror, double gap_size, double far_cut, //////////////////////////////////////////////////////////////////////////////////// -void ELC_P3M_self_forces() { +void ELC_P3M_self_forces(const ParticleRange &particles) { Utils::Vector3d pos; double q; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { if (p.r.p[2] < elc_params.space_layer) { q = elc_params.delta_mid_bot * p.p.q * p.p.q; @@ -1345,7 +1352,7 @@ void ELC_P3M_self_forces() { //////////////////////////////////////////////////////////////////////////////////// -void ELC_p3m_charge_assign_both() { +void ELC_p3m_charge_assign_both(const ParticleRange &particles) { Utils::Vector3d pos; /* charged particle counter, charge fraction counter */ int cp_cnt = 0; @@ -1353,7 +1360,7 @@ void ELC_p3m_charge_assign_both() { for (int i = 0; i < p3m.local_mesh.size; i++) p3m.rs_mesh[i] = 0.0; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { if (p.p.q != 0.0) { p3m_assign_charge(p.p.q, p.r.p, cp_cnt); @@ -1381,13 +1388,13 @@ void ELC_p3m_charge_assign_both() { #endif } -void ELC_p3m_charge_assign_image() { +void ELC_p3m_charge_assign_image(const ParticleRange &particles) { Utils::Vector3d pos; /* prepare local FFT mesh */ for (int i = 0; i < p3m.local_mesh.size; i++) p3m.rs_mesh[i] = 0.0; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { if (p.p.q != 0.0) { if (p.r.p[2] < elc_params.space_layer) { @@ -1511,13 +1518,13 @@ double ELC_P3M_dielectric_layers_energy_contribution(const Particle *p1, ////////////////////////////////////////////////////////////////////////////////// -double ELC_P3M_dielectric_layers_energy_self() { +double ELC_P3M_dielectric_layers_energy_self(const ParticleRange &particles) { Utils::Vector3d pos; double q; double eng = 0.0; // Loop cell neighbors - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { // Loop neighbor cell particles if (p.r.p[2] < elc_params.space_layer) { @@ -1543,7 +1550,7 @@ double ELC_P3M_dielectric_layers_energy_self() { ///////////////////////////////////////////////////////////////////////////////// -void ELC_P3M_modify_p3m_sums_both() { +void ELC_P3M_modify_p3m_sums_both(const ParticleRange &particles) { double node_sums[3], tot_sums[3]; for (int i = 0; i < 3; i++) { @@ -1551,7 +1558,7 @@ void ELC_P3M_modify_p3m_sums_both() { tot_sums[i] = 0.0; } - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { if (p.p.q != 0.0) { node_sums[0] += 1.0; @@ -1579,7 +1586,7 @@ void ELC_P3M_modify_p3m_sums_both() { p3m.square_sum_q = Utils::sqr(tot_sums[2]); } -void ELC_P3M_modify_p3m_sums_image() { +void ELC_P3M_modify_p3m_sums_image(const ParticleRange &particles) { double node_sums[3], tot_sums[3]; for (int i = 0; i < 3; i++) { @@ -1587,7 +1594,7 @@ void ELC_P3M_modify_p3m_sums_image() { tot_sums[i] = 0.0; } - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { if (p.p.q != 0.0) { if (p.r.p[2] < elc_params.space_layer) { @@ -1613,7 +1620,7 @@ void ELC_P3M_modify_p3m_sums_image() { } // this function is required in force.cpp for energy evaluation -void ELC_P3M_restore_p3m_sums() { +void ELC_P3M_restore_p3m_sums(const ParticleRange &particles) { double node_sums[3], tot_sums[3]; for (int i = 0; i < 3; i++) { @@ -1621,7 +1628,7 @@ void ELC_P3M_restore_p3m_sums() { tot_sums[i] = 0.0; } - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { if (p.p.q != 0.0) { node_sums[0] += 1.0; diff --git a/src/core/electrostatics_magnetostatics/elc.hpp b/src/core/electrostatics_magnetostatics/elc.hpp index 0708fed77bf..3ae9a5b5d62 100644 --- a/src/core/electrostatics_magnetostatics/elc.hpp +++ b/src/core/electrostatics_magnetostatics/elc.hpp @@ -107,10 +107,10 @@ int ELC_set_params(double maxPWerror, double min_dist, double far_cut, bool const_pot, double pot_diff); /// the force calculation -void ELC_add_force(); +void ELC_add_force(const ParticleRange &particles); /// the energy calculation -double ELC_energy(); +double ELC_energy(const ParticleRange &particles); /// check the ELC parameters /// @retval ES_OK @@ -132,22 +132,22 @@ void ELC_P3M_dielectric_layers_force_contribution(const Particle *p1, Utils::Vector3d &force1, Utils::Vector3d &force2); /// self energies of top and bottom layers with their virtual images -double ELC_P3M_dielectric_layers_energy_self(); +double ELC_P3M_dielectric_layers_energy_self(const ParticleRange &particles); /// forces of particles in border layers with themselves -void ELC_P3M_self_forces(); +void ELC_P3M_self_forces(const ParticleRange &particles); /// assign the additional, virtual charges, used only in energy.cpp -void ELC_p3m_charge_assign_both(); +void ELC_p3m_charge_assign_both(const ParticleRange &particles); /// assign the additional, virtual charges, used only in energy.cpp -void ELC_p3m_charge_assign_image(); +void ELC_p3m_charge_assign_image(const ParticleRange &particles); /// take into account the virtual charges in the charge sums, used in energy.cpp -void ELC_P3M_modify_p3m_sums_both(); +void ELC_P3M_modify_p3m_sums_both(const ParticleRange &particles); /// take into account the virtual charges in the charge sums, used in energy.cpp -void ELC_P3M_modify_p3m_sums_image(); +void ELC_P3M_modify_p3m_sums_image(const ParticleRange &particles); /// assign the additional, virtual charges, used only in energy.cpp -void ELC_P3M_restore_p3m_sums(); +void ELC_P3M_restore_p3m_sums(const ParticleRange &particles); #endif diff --git a/src/core/electrostatics_magnetostatics/icc.cpp b/src/core/electrostatics_magnetostatics/icc.cpp index ace28c81ad8..60ed48958d4 100644 --- a/src/core/electrostatics_magnetostatics/icc.cpp +++ b/src/core/electrostatics_magnetostatics/icc.cpp @@ -55,7 +55,8 @@ iccp3m_struct iccp3m_cfg; /* functions that are used in icc* to compute the electric field acting on the * induced charges, excluding forces other than the electrostatic ones */ -void init_forces_iccp3m(); +void init_forces_iccp3m(const ParticleRange &particles, + const ParticleRange &ghosts_particles); /** Calculation of the electrostatic forces between source charges (= real * charges) and wall charges. For each electrostatic method the proper functions @@ -63,7 +64,8 @@ void init_forces_iccp3m(); * directly, short range parts need helper functions according to the particle * data organisation. A modified version of \ref force_calc in \ref forces.hpp. */ -void force_calc_iccp3m(); +void force_calc_iccp3m(const ParticleRange &particles, + const ParticleRange &ghost_particles); /** Variant of add_non_bonded_pair_force where only Coulomb * contributions are calculated */ @@ -86,7 +88,8 @@ void iccp3m_alloc_lists() { iccp3m_cfg.sigma.resize(n_ic); } -int iccp3m_iteration(const ParticleRange &particles) { +int iccp3m_iteration(const ParticleRange &particles, + const ParticleRange &ghost_particles) { if (iccp3m_cfg.n_ic == 0) return 0; @@ -105,8 +108,8 @@ int iccp3m_iteration(const ParticleRange &particles) { for (int j = 0; j < iccp3m_cfg.num_iteration; j++) { double hmax = 0.; - force_calc_iccp3m(); /* Calculate electrostatic forces (SR+LR) excluding - source source interaction*/ + force_calc_iccp3m(particles, ghost_particles); /* Calculate electrostatic + forces (SR+LR) excluding source source interaction*/ ghost_communicator(&cell_structure.collect_ghost_force_comm); double diff = 0; @@ -188,8 +191,9 @@ int iccp3m_iteration(const ParticleRange &particles) { return iccp3m_cfg.citeration; } -void force_calc_iccp3m() { - init_forces_iccp3m(); +void force_calc_iccp3m(const ParticleRange &particles, + const ParticleRange &ghost_particles) { + init_forces_iccp3m(particles, ghost_particles); short_range_loop(Utils::NoOp{}, [](Particle &p1, Particle &p2, Distance &d) { /* calc non bonded interactions */ @@ -197,15 +201,16 @@ void force_calc_iccp3m() { d.dist2); }); - Coulomb::calc_long_range_force(); + Coulomb::calc_long_range_force(particles); } -void init_forces_iccp3m() { - for (auto &p : local_cells.particles()) { +void init_forces_iccp3m(const ParticleRange &particles, + const ParticleRange &ghosts_particles) { + for (auto &p : particles) { p.f = ParticleForce{}; } - for (auto &p : ghost_cells.particles()) { + for (auto &p : ghosts_particles) { p.f = ParticleForce{}; } } diff --git a/src/core/electrostatics_magnetostatics/icc.hpp b/src/core/electrostatics_magnetostatics/icc.hpp index 75da1a9543f..2aba8ce3896 100644 --- a/src/core/electrostatics_magnetostatics/icc.hpp +++ b/src/core/electrostatics_magnetostatics/icc.hpp @@ -95,7 +95,8 @@ extern iccp3m_struct iccp3m_cfg; /* global variable with ICCP3M configuration */ /** The main iterative scheme, where the surface element charges are calculated * self-consistently. */ -int iccp3m_iteration(const ParticleRange &particles); +int iccp3m_iteration(const ParticleRange &particles, + const ParticleRange &ghost_particles); /** The allocation of ICCP3M lists for python interface */ diff --git a/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.cpp b/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.cpp index 5f617e3b8c8..cafcc649e43 100644 --- a/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.cpp +++ b/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.cpp @@ -127,7 +127,8 @@ double calc_dipole_dipole_ia(Particle *p1, const Utils::Vector3d &dip1, ============================================================================= */ -double dawaanr_calculations(int force_flag, int energy_flag) { +double dawaanr_calculations(int force_flag, int energy_flag, + const ParticleRange &particles) { if (n_nodes != 1) { fprintf(stderr, "error: DAWAANR is just for one cpu...\n"); @@ -142,7 +143,7 @@ double dawaanr_calculations(int force_flag, int energy_flag) { // Variable to sum up the energy double u = 0; - auto parts = local_cells.particles(); + auto parts = particles; // Iterate over all cells for (auto it = parts.begin(), end = parts.end(); it != end; ++it) { @@ -184,8 +185,9 @@ int magnetic_dipolar_direct_sum_sanity_checks() { /************************************************************/ -double magnetic_dipolar_direct_sum_calculations(int force_flag, - int energy_flag) { +double +magnetic_dipolar_direct_sum_calculations(int force_flag, int energy_flag, + const ParticleRange &particles) { std::vector x, y, z; std::vector mx, my, mz; std::vector fx, fy, fz; @@ -225,7 +227,7 @@ double magnetic_dipolar_direct_sum_calculations(int force_flag, } int dip_particles = 0; - for (auto const &p : local_cells.particles()) { + for (auto const &p : particles) { if (p.p.dipm != 0.0) { const Utils::Vector3d dip = p.calc_dip(); @@ -345,7 +347,7 @@ double magnetic_dipolar_direct_sum_calculations(int force_flag, int dip_particles2 = 0; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { if (p.p.dipm != 0.0) { p.f.f[0] += dipole.prefactor * fx[dip_particles2]; diff --git a/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.hpp b/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.hpp index 51d3013013b..eafab7a1774 100644 --- a/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.hpp +++ b/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.hpp @@ -33,6 +33,7 @@ * */ #include "config.hpp" +#include #ifdef DIPOLES #include "particle_data.hpp" @@ -47,7 +48,8 @@ double calc_dipole_dipole_ia(Particle *p1, Particle *p2, int force_flag); /** Core of the DAWAANR method: here you compute all the magnetic forces, * torques and the magnetic energy for the whole system*/ -double dawaanr_calculations(int force_flag, int energy_flag); +double dawaanr_calculations(int force_flag, int energy_flag, + const ParticleRange &particles); /** switch on DAWAANR magnetostatics. @return ES_ERROR, if not on a single CPU @@ -64,8 +66,8 @@ int magnetic_dipolar_direct_sum_sanity_checks(); /* Core of the method: here you compute all the magnetic forces,torques and the * energy for the whole system using direct sum*/ -double magnetic_dipolar_direct_sum_calculations(int force_flag, - int energy_flag); +double magnetic_dipolar_direct_sum_calculations(int force_flag, int energy_flag, + const ParticleRange &particles); /** switch on direct sum magnetostatics. @param n_cut cut off for the explicit summation diff --git a/src/core/electrostatics_magnetostatics/mdlc_correction.cpp b/src/core/electrostatics_magnetostatics/mdlc_correction.cpp index 2f59fc84e30..7612daf4f54 100644 --- a/src/core/electrostatics_magnetostatics/mdlc_correction.cpp +++ b/src/core/electrostatics_magnetostatics/mdlc_correction.cpp @@ -52,8 +52,9 @@ static int n_local_particles = 0; static double mu_max; void calc_mu_max() { + auto local_particles = local_cells.particles(); mu_max = std::accumulate( - local_cells.particles().begin(), local_cells.particles().end(), 0.0, + local_particles.begin(), local_particles.end(), 0.0, [](double mu, Particle const &p) { return std::max(mu, p.p.dipm); }); MPI_Allreduce(MPI_IN_PLACE, &mu_max, 1, MPI_DOUBLE, MPI_MAX, comm_cart); @@ -76,11 +77,12 @@ inline double g2_DLC_dip(double g, double x) { } /* Compute Mx, My, Mz and Mtotal */ -double slab_dip_count_mu(double *mt, double *mx, double *my) { +double slab_dip_count_mu(double *mt, double *mx, double *my, + const ParticleRange &particles) { Utils::Vector3d node_sums{}; Utils::Vector3d tot_sums{}; - for (auto const &p : local_cells.particles()) { + for (auto const &p : particles) { if (p.p.dipm != 0.0) { node_sums += p.calc_dip(); } @@ -107,7 +109,8 @@ double slab_dip_count_mu(double *mt, double *mx, double *my) { Lett. 400, 62-67, (2004). */ double get_DLC_dipolar(int kcut, std::vector &fs, - std::vector &ts) { + std::vector &ts, + const ParticleRange &particles) { int ip; @@ -142,7 +145,7 @@ double get_DLC_dipolar(int kcut, std::vector &fs, // {Re(S+), Im(S+), Re(S-), Im(S-)} ip = 0; - for (auto const &p : local_cells.particles()) { + for (auto const &p : particles) { if (p.p.dipm > 0) { Utils::Vector3d const dip = p.calc_dip(); @@ -182,7 +185,7 @@ double get_DLC_dipolar(int kcut, std::vector &fs, // ... Now we can compute the contributions to E,Fj,Ej for the current // g-value ip = 0; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { if (p.p.dipm > 0) { // We compute the contributions to the forces ............ @@ -234,7 +237,7 @@ double get_DLC_dipolar(int kcut, std::vector &fs, //-tz[0]*M_PI/(box_l[0]*box_l[1]) ); ip = 0; - for (auto const &p : local_cells.particles()) { + for (auto const &p : particles) { if (p.p.dipm > 0) { ts[ip] = vector_product(p.calc_dip(), ts[ip]); } @@ -266,7 +269,7 @@ double get_DLC_dipolar(int kcut, std::vector &fs, Algorithm implemented accordingly to the paper of A. Brodka, Chem. Phys. Lett. 400, 62-67, (2004). */ -double get_DLC_energy_dipolar(int kcut) { +double get_DLC_energy_dipolar(int kcut, const ParticleRange &particles) { int ix, iy, ip; double gx, gy, gr; @@ -276,7 +279,7 @@ double get_DLC_energy_dipolar(int kcut) { double s1; double energy, piarea, facux, facuy; - n_local_particles = local_cells.particles().size(); + n_local_particles = particles.size(); facux = 2.0 * M_PI / box_geo.length()[0]; facuy = 2.0 * M_PI / box_geo.length()[1]; @@ -302,7 +305,7 @@ double get_DLC_energy_dipolar(int kcut) { ip = 0; - for (auto const &p : local_cells.particles()) { + for (auto const &p : particles) { if (p.p.dipm > 0) { const Utils::Vector3d dip = p.calc_dip(); @@ -346,11 +349,11 @@ double get_DLC_energy_dipolar(int kcut) { /** Compute and add the terms needed to correct the 3D dipolar * methods when we have an slab geometry */ -void add_mdlc_force_corrections() { +void add_mdlc_force_corrections(const ParticleRange &particles) { int dip_DLC_kcut = dlc_params.far_cut; double mz = 0.0, mx = 0.0, my = 0.0, volume, mtot = 0.0; - n_local_particles = local_cells.particles().size(); + n_local_particles = particles.size(); volume = box_geo.length()[0] * box_geo.length()[1] * box_geo.length()[2]; @@ -363,7 +366,7 @@ void add_mdlc_force_corrections() { //---- Compute the corrections ---------------------------------- // First the DLC correction - get_DLC_dipolar(dip_DLC_kcut, dip_DLC_f, dip_DLC_t); + get_DLC_dipolar(dip_DLC_kcut, dip_DLC_f, dip_DLC_t, particles); // Now we compute the the correction like Yeh and Klapp to take into account // the fact that you are using a @@ -374,13 +377,13 @@ void add_mdlc_force_corrections() { // Shape Dependent Correction. // See Brodka, Chem. Phys. Lett. 400, 62, (2004). - mz = slab_dip_count_mu(&mtot, &mx, &my); + mz = slab_dip_count_mu(&mtot, &mx, &my, particles); // --- Transfer the computed corrections to the Forces, Energy and torques // of the particles int ip = 0; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { if ((p.p.dipm) != 0.0) { // SDC correction term is zero for the forces p.f.f += dipole.prefactor * dip_DLC_f[ip]; @@ -409,7 +412,7 @@ void add_mdlc_force_corrections() { /** Compute and add the terms needed to correct the energy of * 3D dipolar methods when we have an slab geometry */ -double add_mdlc_energy_corrections() { +double add_mdlc_energy_corrections(const ParticleRange &particles) { double dip_DLC_energy = 0.0; double mz = 0.0, mx = 0.0, my = 0.0, volume, mtot = 0.0; int dip_DLC_kcut; @@ -421,7 +424,8 @@ double add_mdlc_energy_corrections() { //---- Compute the corrections ---------------------------------- // First the DLC correction - dip_DLC_energy += dipole.prefactor * get_DLC_energy_dipolar(dip_DLC_kcut); + dip_DLC_energy += + dipole.prefactor * get_DLC_energy_dipolar(dip_DLC_kcut, particles); // printf("Energy DLC = %20.15le // \n",dip_DLC_energy); @@ -435,7 +439,7 @@ double add_mdlc_energy_corrections() { // Shape Dependent Correction. // See Brodka, Chem. Phys. Lett. 400, 62, (2004). - mz = slab_dip_count_mu(&mtot, &mx, &my); + mz = slab_dip_count_mu(&mtot, &mx, &my, particles); if (this_node == 0) { #ifdef DP3M diff --git a/src/core/electrostatics_magnetostatics/mdlc_correction.hpp b/src/core/electrostatics_magnetostatics/mdlc_correction.hpp index a248139129f..7dc37864aaf 100644 --- a/src/core/electrostatics_magnetostatics/mdlc_correction.hpp +++ b/src/core/electrostatics_magnetostatics/mdlc_correction.hpp @@ -40,6 +40,7 @@ #define _DLC_DIPOLAR_H #include "config.hpp" +#include #if defined(DIPOLES) && defined(DP3M) @@ -76,8 +77,8 @@ extern DLC_struct dlc_params; int mdlc_set_params(double maxPWerror, double gap_size, double far_cut); int mdlc_sanity_checks(); -void add_mdlc_force_corrections(); -double add_mdlc_energy_corrections(); +void add_mdlc_force_corrections(const ParticleRange &particles); +double add_mdlc_energy_corrections(const ParticleRange &particles); void calc_mu_max(); #endif diff --git a/src/core/electrostatics_magnetostatics/mmm2d.cpp b/src/core/electrostatics_magnetostatics/mmm2d.cpp index 604e448413d..9074ea48050 100644 --- a/src/core/electrostatics_magnetostatics/mmm2d.cpp +++ b/src/core/electrostatics_magnetostatics/mmm2d.cpp @@ -228,7 +228,7 @@ static void prepareBernoulliNumbers(int nmax); static int MMM2D_tune_near(double error); /** energy of all local particles with their copies */ -void MMM2D_self_energy(); +void MMM2D_self_energy(const ParticleRange &particles); /*@}*/ @@ -582,7 +582,7 @@ static void setup_z_force() { } } -static void add_z_force() { +static void add_z_force(const ParticleRange &particles) { double add; double *othcblk; int size = 2; @@ -592,7 +592,7 @@ static void add_z_force() { if (mmm2d_params.const_pot_on) { double gbl_dm_z = 0; double lcl_dm_z = 0; - for (auto const &p : local_cells.particles()) { + for (auto const &p : particles) { lcl_dm_z += p.p.q * (p.r.p[2] + p.l.i[2] * box_geo.length()[2]); } @@ -650,7 +650,7 @@ static void setup_z_energy() { } } -static double z_energy() { +static double z_energy(const ParticleRange &particles) { int np, c, i; Particle *part; double *othcblk; @@ -671,7 +671,7 @@ static double z_energy() { double gbl_dm_z = 0; double lcl_dm_z = 0; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { lcl_dm_z += p.p.q * (p.r.p[2] + p.l.i[2] * box_geo.length()[2]); } @@ -1118,7 +1118,8 @@ static double PQ_energy(double omega) { /* main loops */ /*****************************************************************/ -static void add_force_contribution(int p, int q) { +static void add_force_contribution(int p, int q, + const ParticleRange &particles) { double omega, fac; if (q == 0) { @@ -1130,7 +1131,7 @@ static void add_force_contribution(int p, int q) { distribute(1, 1.); - add_z_force(); + add_z_force(particles); checkpoint("************2piz", 0, 0, 1); } else { @@ -1170,7 +1171,8 @@ static void add_force_contribution(int p, int q) { } } -static double energy_contribution(int p, int q) { +static double energy_contribution(int p, int q, + const ParticleRange &particles) { double eng; double omega, fac; @@ -1179,7 +1181,7 @@ static double energy_contribution(int p, int q) { setup_z_energy(); clear_image_contributions(2); distribute(2, 1.); - eng = z_energy(); + eng = z_energy(particles); checkpoint("E************2piz", 0, 0, 2); } else { omega = C_2PI * ux * p; @@ -1219,7 +1221,7 @@ static double energy_contribution(int p, int q) { return eng; } -double MMM2D_add_far(int f, int e) { +double MMM2D_add_far(int f, int e, const ParticleRange &particles) { int p, q; double R, dR, q2; @@ -1262,9 +1264,9 @@ double MMM2D_add_far(int f, int e) { if (ux2 * Utils::sqr(p) + uy2 * Utils::sqr(q) < Utils::sqr(R)) break; if (f) - add_force_contribution(p, q); + add_force_contribution(p, q, particles); if (e) - eng += energy_contribution(p, q); + eng += energy_contribution(p, q, particles); } undone[p] = q; } @@ -1277,9 +1279,9 @@ double MMM2D_add_far(int f, int e) { for (; q >= 0; q--) { // printf("xxxxx %d %d\n", p, q); if (f) - add_force_contribution(p, q); + add_force_contribution(p, q, particles); if (e) - eng += energy_contribution(p, q); + eng += energy_contribution(p, q, particles); } } @@ -1699,7 +1701,7 @@ double mmm2d_coulomb_pair_energy(double charge_factor, return 0.0; } -void MMM2D_self_energy() { +void MMM2D_self_energy(const ParticleRange &particles) { Utils::Vector3d dv{}; double seng = coulomb.prefactor * calc_mmm2d_copy_pair_energy(dv); @@ -1707,7 +1709,7 @@ void MMM2D_self_energy() { in the far formula which counts everything twice and in the end divides by two*/ - auto parts = local_cells.particles(); + auto parts = particles; self_energy = std::accumulate(parts.begin(), parts.end(), 0.0, [seng](double sum, Particle const &p) { return sum + seng * Utils::sqr(p.p.q); @@ -1835,7 +1837,7 @@ void MMM2D_init() { } } -void MMM2D_on_resort_particles() { +void MMM2D_on_resort_particles(const ParticleRange &particles) { /* if we need MMM2D far formula, allocate caches */ if (cell_structure.type == CELL_STRUCTURE_LAYERED) { n_localpart = cells_get_n_particles(); @@ -1848,7 +1850,7 @@ void MMM2D_on_resort_particles() { lclcblk.resize(cells.size() * 8); gblcblk.resize(n_layers * 8); } - MMM2D_self_energy(); + MMM2D_self_energy(particles); } void MMM2D_dielectric_layers_force_contribution() { diff --git a/src/core/electrostatics_magnetostatics/mmm2d.hpp b/src/core/electrostatics_magnetostatics/mmm2d.hpp index 04c89c8d740..015f57ef8d5 100644 --- a/src/core/electrostatics_magnetostatics/mmm2d.hpp +++ b/src/core/electrostatics_magnetostatics/mmm2d.hpp @@ -33,6 +33,7 @@ #include "config.hpp" +#include #include #ifdef ELECTROSTATICS @@ -88,13 +89,17 @@ int MMM2D_set_params(double maxPWerror, double far_cut, double delta_top, double delta_bot, bool const_pot_on, double pot_diff); /** the general long range force/energy calculation */ -double MMM2D_add_far(int f, int e); +double MMM2D_add_far(int f, int e, const ParticleRange &particles); /** the actual long range force calculation */ -inline void MMM2D_add_far_force() { MMM2D_add_far(1, 0); } +inline void MMM2D_add_far_force(const ParticleRange &particles) { + MMM2D_add_far(1, 0, particles); +} /** the actual long range energy calculation */ -inline double MMM2D_far_energy() { return MMM2D_add_far(0, 1); } +inline double MMM2D_far_energy(const ParticleRange &particles) { + return MMM2D_add_far(0, 1, particles); +} /** pairwise calculated parts of MMM2D force (near neighbors) */ void add_mmm2d_coulomb_pair_force(double pref, Utils::Vector3d const &d, @@ -112,7 +117,7 @@ void MMM2D_init(); /** if the number of particles has changed (even per node), the particle buffers for the coefficients have to be resized. */ -void MMM2D_on_resort_particles(); +void MMM2D_on_resort_particles(const ParticleRange &particles); /** energy contribution from dielectric layers */ double MMM2D_dielectric_layers_energy_contribution(); diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp index 40b16b92107..a83603dff37 100644 --- a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp +++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp @@ -181,7 +181,8 @@ static double dp3m_k_space_error(double box_size, double prefac, int mesh, /*@}*/ /** Compute the dipolar surface terms */ -static double calc_surface_term(int force_flag, int energy_flag); +static double calc_surface_term(int force_flag, int energy_flag, + const ParticleRange &particles); /** \name P3M Tuning Functions */ /************************************************************/ @@ -636,7 +637,7 @@ void dp3m_interpolate_dipole_assignment_function() { } } -void dp3m_dipole_assign() { +void dp3m_dipole_assign(const ParticleRange &particles) { /* magnetic particle counter, dipole fraction counter */ int cp_cnt = 0; @@ -645,7 +646,7 @@ void dp3m_dipole_assign() { for (int j = 0; j < dp3m.local_mesh.size; j++) i[j] = 0.0; - for (auto const &p : local_cells.particles()) { + for (auto const &p : particles) { if (p.p.dipm != 0.0) { dp3m_assign_dipole(p.r.p.data(), p.p.dipm, p.calc_dip().data(), cp_cnt); cp_cnt++; @@ -761,7 +762,8 @@ void dp3m_shrink_wrap_dipole_grid(int n_dipoles) { #ifdef ROTATION /* Assign the torques obtained from k-space */ -static void P3M_assign_torques(double prefac, int d_rs) { +static void P3M_assign_torques(double prefac, int d_rs, + const ParticleRange &particles) { /* particle counter, charge fraction counter */ int cp_cnt = 0, cf_cnt = 0; /* index, index jumps for dp3m.rs_mesh array */ @@ -773,7 +775,7 @@ static void P3M_assign_torques(double prefac, int d_rs) { cp_cnt = 0; cf_cnt = 0; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { if ((p.p.dipm) != 0.0) { const Utils::Vector3d dip = p.calc_dip(); q_ind = dp3m.ca_fmp[cp_cnt]; @@ -825,7 +827,8 @@ static void P3M_assign_torques(double prefac, int d_rs) { #endif /* Assign the dipolar forces obtained from k-space */ -static void dp3m_assign_forces_dip(double prefac, int d_rs) { +static void dp3m_assign_forces_dip(double prefac, int d_rs, + const ParticleRange &particles) { /* particle counter, charge fraction counter */ int cp_cnt = 0, cf_cnt = 0; /* index, index jumps for dp3m.rs_mesh array */ @@ -837,7 +840,7 @@ static void dp3m_assign_forces_dip(double prefac, int d_rs) { cp_cnt = 0; cf_cnt = 0; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { if ((p.p.dipm) != 0.0) { const Utils::Vector3d dip = p.calc_dip(); q_ind = dp3m.ca_fmp[cp_cnt]; @@ -866,7 +869,8 @@ static void dp3m_assign_forces_dip(double prefac, int d_rs) { /*****************************************************************************/ -double dp3m_calc_kspace_forces(int force_flag, int energy_flag) { +double dp3m_calc_kspace_forces(int force_flag, int energy_flag, + const ParticleRange &particles) { int i, d, d_rs, ind, j[3]; /**************************************************************/ /* k-space energy */ @@ -1042,8 +1046,9 @@ double dp3m_calc_kspace_forces(int force_flag, int energy_flag) { /* redistribute force component mesh */ dp3m_spread_force_grid(dp3m.rs_mesh); /* Assign force component from mesh to particle */ - P3M_assign_torques( - dipole_prefac * (2 * Utils::pi() / box_geo.length()[0]), d_rs); + P3M_assign_torques(dipole_prefac * + (2 * Utils::pi() / box_geo.length()[0]), + d_rs, particles); } P3M_TRACE(fprintf(stderr, "%d: done torque calculation.\n", this_node)); #endif /*if def ROTATION */ @@ -1129,8 +1134,8 @@ double dp3m_calc_kspace_forces(int force_flag, int energy_flag) { dp3m_spread_force_grid(dp3m.rs_mesh_dip[2]); /* Assign force component from mesh to particle */ dp3m_assign_forces_dip( - dipole_prefac * pow(2 * Utils::pi() / box_geo.length()[0], 2), - d_rs); + dipole_prefac * pow(2 * Utils::pi() / box_geo.length()[0], 2), d_rs, + particles); } P3M_TRACE(fprintf(stderr, @@ -1141,7 +1146,7 @@ double dp3m_calc_kspace_forces(int force_flag, int energy_flag) { } /* of if(force_flag) */ if (dp3m.params.epsilon != P3M_EPSILON_METALLIC) { - surface_term = calc_surface_term(force_flag, energy_flag); + surface_term = calc_surface_term(force_flag, energy_flag, particles); if (this_node == 0) k_space_energy_dip += surface_term; } @@ -1151,14 +1156,15 @@ double dp3m_calc_kspace_forces(int force_flag, int energy_flag) { /************************************************************/ -double calc_surface_term(int force_flag, int energy_flag) { +double calc_surface_term(int force_flag, int energy_flag, + const ParticleRange &particles) { const double pref = dipole.prefactor * 4 * M_PI * (1. / box_geo.length()[0]) * (1. / box_geo.length()[1]) * (1. / box_geo.length()[2]) / (2 * dp3m.params.epsilon + 1); double suma, a[3]; double en; - auto const n_local_part = local_cells.particles().size(); + auto const n_local_part = particles.size(); // We put all the dipolar momenta in a the arrays mx,my,mz according to the // id-number of the particles @@ -1167,7 +1173,7 @@ double calc_surface_term(int force_flag, int energy_flag) { std::vector mz(n_local_part); int ip = 0; - for (auto const &p : local_cells.particles()) { + for (auto const &p : particles) { const Utils::Vector3d dip = p.calc_dip(); mx[ip] = dip[0]; my[ip] = dip[1]; @@ -1214,7 +1220,7 @@ double calc_surface_term(int force_flag, int energy_flag) { } ip = 0; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { p.f.torque[0] -= pref * sumix[ip]; p.f.torque[1] -= pref * sumiy[ip]; p.f.torque[2] -= pref * sumiz[ip]; diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp index 976a3be5ab8..ac2211c209c 100644 --- a/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp +++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp @@ -44,6 +44,7 @@ #include "p3m-common.hpp" #include "particle_data.hpp" +#include #include #include @@ -145,7 +146,7 @@ bool dp3m_sanity_checks(const Utils::Vector3i &grid); * If Dstore_ca_frac is true, then the charge fractions are buffered in * Dcur_ca_fmp and Dcur_ca_frac. */ -void dp3m_dipole_assign(); +void dp3m_dipole_assign(const ParticleRange &particles); /** Reset @ref dp3m core parameters */ void dp3m_deactivate(); @@ -193,7 +194,8 @@ int dp3m_adaptive_tune(char **log); /** Compute the k-space part of forces and energies for the magnetic * dipole-dipole interaction */ -double dp3m_calc_kspace_forces(int force_flag, int energy_flag); +double dp3m_calc_kspace_forces(int force_flag, int energy_flag, + const ParticleRange &particles); /** Calculate number of magnetic particles, the sum of the squared * charges and the squared sum of the charges. diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp index e6b581dd2dd..3c649442adb 100644 --- a/src/core/electrostatics_magnetostatics/p3m.cpp +++ b/src/core/electrostatics_magnetostatics/p3m.cpp @@ -132,7 +132,8 @@ static void p3m_init_a_ai_cao_cut(); static void p3m_calc_lm_ld_pos(); /** Calculates the dipole term */ -static double p3m_calc_dipole_term(int force_flag, int energy_flag); +static double p3m_calc_dipole_term(int force_flag, int energy_flag, + const ParticleRange &particles); /** Gather FFT grid. * After the charge assignment Each node needs to gather the @@ -251,7 +252,8 @@ static void p3m_tune_aliasing_sums(int nx, int ny, int nz, const int mesh[3], * wrapper. * \tparam cao charge assignment order. */ -template static void p3m_do_charge_assign(); +template +static void p3m_do_charge_assign(const ParticleRange &particles); template void p3m_do_assign_charge(double q, Utils::Vector3d &real_pos, int cp_cnt); @@ -493,41 +495,41 @@ void p3m_interpolate_charge_assignment_function() { } /* Template wrapper for p3m_do_charge_assign() */ -void p3m_charge_assign() { +void p3m_charge_assign(const ParticleRange &particles) { switch (p3m.params.cao) { case 1: - p3m_do_charge_assign<1>(); + p3m_do_charge_assign<1>(particles); break; case 2: - p3m_do_charge_assign<2>(); + p3m_do_charge_assign<2>(particles); break; case 3: - p3m_do_charge_assign<3>(); + p3m_do_charge_assign<3>(particles); break; case 4: - p3m_do_charge_assign<4>(); + p3m_do_charge_assign<4>(particles); break; case 5: - p3m_do_charge_assign<5>(); + p3m_do_charge_assign<5>(particles); break; case 6: - p3m_do_charge_assign<6>(); + p3m_do_charge_assign<6>(particles); break; case 7: - p3m_do_charge_assign<7>(); + p3m_do_charge_assign<7>(particles); break; } } /** Assign the charges */ -template void p3m_do_charge_assign() { +template void p3m_do_charge_assign(const ParticleRange &particles) { /* charged particle counter, charge fraction counter */ int cp_cnt = 0; /* prepare local FFT mesh */ for (int i = 0; i < p3m.local_mesh.size; i++) p3m.rs_mesh[i] = 0.0; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { if (p.p.q != 0.0) { p3m_do_assign_charge(p.p.q, p.r.p, cp_cnt); cp_cnt++; @@ -673,7 +675,8 @@ void p3m_shrink_wrap_charge_grid(int n_charges) { /* Assign the forces obtained from k-space */ template -static void P3M_assign_forces(double force_prefac, int d_rs) { +static void P3M_assign_forces(double force_prefac, int d_rs, + const ParticleRange &particles) { /* charged particle counter, charge fraction counter */ int cp_cnt = 0; #ifdef P3M_STORE_CA_FRAC @@ -687,7 +690,7 @@ static void P3M_assign_forces(double force_prefac, int d_rs) { /* index, index jumps for rs_mesh array */ int q_ind = 0; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { auto const q = p.p.q; if (q != 0.0) { #ifdef P3M_STORE_CA_FRAC @@ -767,7 +770,8 @@ static void P3M_assign_forces(double force_prefac, int d_rs) { } } -double p3m_calc_kspace_forces(int force_flag, int energy_flag) { +double p3m_calc_kspace_forces(int force_flag, int energy_flag, + const ParticleRange &particles) { int i, d, d_rs, ind, j[3]; /**************************************************************/ /* Prefactor for force */ @@ -880,32 +884,32 @@ double p3m_calc_kspace_forces(int force_flag, int energy_flag) { /* Assign force component from mesh to particle */ switch (p3m.params.cao) { case 1: - P3M_assign_forces<1>(force_prefac, d_rs); + P3M_assign_forces<1>(force_prefac, d_rs, particles); break; case 2: - P3M_assign_forces<2>(force_prefac, d_rs); + P3M_assign_forces<2>(force_prefac, d_rs, particles); break; case 3: - P3M_assign_forces<3>(force_prefac, d_rs); + P3M_assign_forces<3>(force_prefac, d_rs, particles); break; case 4: - P3M_assign_forces<4>(force_prefac, d_rs); + P3M_assign_forces<4>(force_prefac, d_rs, particles); break; case 5: - P3M_assign_forces<5>(force_prefac, d_rs); + P3M_assign_forces<5>(force_prefac, d_rs, particles); break; case 6: - P3M_assign_forces<6>(force_prefac, d_rs); + P3M_assign_forces<6>(force_prefac, d_rs, particles); break; case 7: - P3M_assign_forces<7>(force_prefac, d_rs); + P3M_assign_forces<7>(force_prefac, d_rs, particles); break; } } } /* if(force_flag) */ if (p3m.params.epsilon != P3M_EPSILON_METALLIC) { - k_space_energy += p3m_calc_dipole_term(force_flag, energy_flag); + k_space_energy += p3m_calc_dipole_term(force_flag, energy_flag, particles); } return k_space_energy; @@ -913,7 +917,8 @@ double p3m_calc_kspace_forces(int force_flag, int energy_flag) { /************************************************************/ -double p3m_calc_dipole_term(int force_flag, int energy_flag) { +double p3m_calc_dipole_term(int force_flag, int energy_flag, + const ParticleRange &particles) { double pref = coulomb.prefactor * 4 * M_PI * (1. / box_geo.length()[0]) * (1. / box_geo.length()[1]) * (1. / box_geo.length()[2]) / (2 * p3m.params.epsilon + 1); @@ -923,7 +928,7 @@ double p3m_calc_dipole_term(int force_flag, int energy_flag) { for (double &j : lcl_dm) j = 0; - for (auto const &p : local_cells.particles()) { + for (auto const &p : particles) { for (int j = 0; j < 3; j++) /* dipole moment with unfolded coordinates */ lcl_dm[j] += p.p.q * (p.r.p[j] + p.l.i[j] * box_geo.length()[j]); @@ -940,7 +945,7 @@ double p3m_calc_dipole_term(int force_flag, int energy_flag) { if (force_flag) { for (double &j : gbl_dm) j *= pref; - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { for (int j = 0; j < 3; j++) p.f.f[j] -= gbl_dm[j] * p.p.q; } diff --git a/src/core/electrostatics_magnetostatics/p3m.hpp b/src/core/electrostatics_magnetostatics/p3m.hpp index de1c993ba4d..6067092267b 100644 --- a/src/core/electrostatics_magnetostatics/p3m.hpp +++ b/src/core/electrostatics_magnetostatics/p3m.hpp @@ -59,6 +59,7 @@ #include "fft.hpp" #include "p3m-common.hpp" +#include #include #include @@ -182,7 +183,8 @@ void p3m_scaleby_box_l(); /** Compute the k-space part of forces and energies for the charge-charge * interaction */ -double p3m_calc_kspace_forces(int force_flag, int energy_flag); +double p3m_calc_kspace_forces(int force_flag, int energy_flag, + const ParticleRange &particles); /** Compute the k-space part of the stress tensor **/ void p3m_calc_kspace_stress(double *stress); @@ -200,7 +202,7 @@ void p3m_count_charged_particles(); * in @ref p3m_data_struct::ca_fmp "ca_fmp" and @ref p3m_data_struct::ca_frac * "ca_frac". */ -void p3m_charge_assign(); +void p3m_charge_assign(const ParticleRange &particles); /** Assign a single charge into the current charge grid. * diff --git a/src/core/energy.cpp b/src/core/energy.cpp index a5f9daf0ff7..8ade508d33d 100644 --- a/src/core/energy.cpp +++ b/src/core/energy.cpp @@ -107,7 +107,7 @@ void energy_calc(double *result) { add_single_particle_energy(&p); } } - calc_long_range_energies(); + calc_long_range_energies(local_cells.particles()); auto local_parts = local_cells.particles(); Constraints::constraints.add_energy(local_parts, sim_time, energy); @@ -123,14 +123,14 @@ void energy_calc(double *result) { /************************************************************/ -void calc_long_range_energies() { +void calc_long_range_energies(const ParticleRange &particles) { #ifdef ELECTROSTATICS /* calculate k-space part of electrostatic interaction. */ - Coulomb::calc_energy_long_range(energy); + Coulomb::calc_energy_long_range(energy, particles); #endif /* ifdef ELECTROSTATICS */ #ifdef DIPOLES - Dipole::calc_energy_long_range(energy); + Dipole::calc_energy_long_range(energy, particles); #endif /* ifdef DIPOLES */ } diff --git a/src/core/energy.hpp b/src/core/energy.hpp index 9d2627cd3ba..ce33536fdd3 100644 --- a/src/core/energy.hpp +++ b/src/core/energy.hpp @@ -56,7 +56,7 @@ void master_energy_calc(); void energy_calc(double *result); /** Calculate long range energies (P3M, MMM2d...). */ -void calc_long_range_energies(); +void calc_long_range_energies(const ParticleRange &particles); /** Calculate the total energy */ double calculate_current_potential_energy_of_system(); diff --git a/src/core/event.cpp b/src/core/event.cpp index 1041dcb26f7..7e7a876bb85 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -298,10 +298,10 @@ void on_lbboundary_change() { #endif } -void on_resort_particles() { +void on_resort_particles(const ParticleRange &particles) { EVENT_TRACE(fprintf(stderr, "%d: on_resort_particles\n", this_node)); #ifdef ELECTROSTATICS - Coulomb::on_resort_particles(); + Coulomb::on_resort_particles(particles); #endif /* ifdef ELECTROSTATICS */ /* DIPOLAR interactions so far don't need this */ diff --git a/src/core/event.hpp b/src/core/event.hpp index 329b79dcd46..cf53dc6f74b 100644 --- a/src/core/event.hpp +++ b/src/core/event.hpp @@ -37,6 +37,7 @@ * * Implementation in event.cpp. */ +#include "ParticleRange.hpp" /** \name Hook procedures * These procedures are called if several significant changes to @@ -69,7 +70,7 @@ void on_particle_change(); void on_particle_charge_change(); /** called every time the particles are resorted from node to node. */ -void on_resort_particles(); +void on_resort_particles(const ParticleRange &particles); /** called every time the Coulomb parameters are changed. */ void on_coulomb_change(); diff --git a/src/core/forces.cpp b/src/core/forces.cpp index a030809ed42..f83d177cd59 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -92,7 +92,7 @@ void force_calc(CellStructure &cell_structure) { auto particles = cell_structure.local_cells().particles(); #ifdef ELECTROSTATICS - iccp3m_iteration(particles); + iccp3m_iteration(particles, cell_structure.ghost_cells().particles()); #endif init_forces(particles); @@ -103,7 +103,7 @@ void force_calc(CellStructure &cell_structure) { #endif } - calc_long_range_forces(); + calc_long_range_forces(particles); // Only calculate pair forces if the maximum cutoff is >0 if (max_cut > 0) { @@ -143,7 +143,7 @@ void force_calc(CellStructure &cell_structure) { // Must be done here. Forces need to be ghost-communicated immersed_boundaries.volume_conservation(); - lb_lbcoupling_calc_particle_lattice_ia(thermo_virtual); + lb_lbcoupling_calc_particle_lattice_ia(thermo_virtual, particles); #ifdef METADYNAMICS /* Metadynamics main function */ @@ -177,16 +177,16 @@ void force_calc(CellStructure &cell_structure) { recalc_forces = 0; } -void calc_long_range_forces() { +void calc_long_range_forces(const ParticleRange &particles) { ESPRESSO_PROFILER_CXX_MARK_FUNCTION; #ifdef ELECTROSTATICS /* calculate k-space part of electrostatic interaction. */ - Coulomb::calc_long_range_force(); + Coulomb::calc_long_range_force(particles); #endif /*ifdef ELECTROSTATICS */ #ifdef DIPOLES /* calculate k-space part of the magnetostatic interaction. */ - Dipole::calc_long_range_force(); + Dipole::calc_long_range_force(particles); #endif /*ifdef DIPOLES */ } diff --git a/src/core/forces.hpp b/src/core/forces.hpp index 19b41a2ba18..670b42e9ee0 100644 --- a/src/core/forces.hpp +++ b/src/core/forces.hpp @@ -64,7 +64,7 @@ void init_forces_ghosts(const ParticleRange &particles); void force_calc(CellStructure &cell_structure); /** Calculate long range forces (P3M, MMM2d...). */ -void calc_long_range_forces(); +void calc_long_range_forces(const ParticleRange &particles); /*@}*/ #endif diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.cpp b/src/core/grid_based_algorithms/lb_particle_coupling.cpp index a8202225a8c..fca486b67a0 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.cpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.cpp @@ -182,7 +182,8 @@ void add_swimmer_force(Particle &p) { #endif } // namespace -void lb_lbcoupling_calc_particle_lattice_ia(bool couple_virtual) { +void lb_lbcoupling_calc_particle_lattice_ia(bool couple_virtual, + const ParticleRange &particles) { ESPRESSO_PROFILER_CXX_MARK_FUNCTION; if (lattice_switch == ActiveLB::GPU) { #ifdef CUDA @@ -244,7 +245,7 @@ void lb_lbcoupling_calc_particle_lattice_ia(bool couple_virtual) { }; /* local cells */ - for (auto &p : local_cells.particles()) { + for (auto &p : particles) { if (!p.p.is_virtual or couple_virtual) { auto const force = lb_viscous_coupling( &p, noise_amplitude * f_random(p.identity())); diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.hpp b/src/core/grid_based_algorithms/lb_particle_coupling.hpp index 2ff1cbbe3f0..387ff734b60 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.hpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.hpp @@ -3,6 +3,7 @@ #include +#include "ParticleRange.hpp" #include /** Calculate particle lattice interactions. @@ -10,7 +11,8 @@ * Include all particle-lattice forces in this function. * The function is called from \ref force_calc. */ -void lb_lbcoupling_calc_particle_lattice_ia(bool couple_virtual); +void lb_lbcoupling_calc_particle_lattice_ia(bool couple_virtual, + const ParticleRange &particles); void lb_lbcoupling_propagate(); uint64_t lb_lbcoupling_get_rng_state(); void lb_lbcoupling_set_rng_state(uint64_t counter); diff --git a/src/core/io/writer/h5md_core.cpp b/src/core/io/writer/h5md_core.cpp index 30d88712555..814d91f0994 100644 --- a/src/core/io/writer/h5md_core.cpp +++ b/src/core/io/writer/h5md_core.cpp @@ -392,7 +392,8 @@ void File::fill_arrays_for_h5md_write_with_particle_property( } } -void File::Write(int write_dat, PartCfg &partCfg) { +void File::Write(int write_dat, PartCfg &partCfg, + const ParticleRange &particles) { #ifdef H5MD_DEBUG std::cout << "Called " << __func__ << " on node " << this_node << std::endl; #endif @@ -441,7 +442,7 @@ void File::Write(int write_dat, PartCfg &partCfg) { /* loop over all local cells. */ int particle_index = 0; - for (auto ¤t_particle : local_cells.particles()) { + for (auto ¤t_particle : particles) { fill_arrays_for_h5md_write_with_particle_property( particle_index, id, typ, mass, pos, image, vel, f, charge, current_particle, write_dat, bond); diff --git a/src/core/io/writer/h5md_core.hpp b/src/core/io/writer/h5md_core.hpp index c5bf838152b..6eba42a2be4 100644 --- a/src/core/io/writer/h5md_core.hpp +++ b/src/core/io/writer/h5md_core.hpp @@ -73,7 +73,7 @@ class File { * write methods. * Boolean values for position, velocity, force and mass. */ - void Write(int write_dat, PartCfg &partCfg); + void Write(int write_dat, PartCfg &partCfg, const ParticleRange &particles); std::string &filename() { return m_filename; }; std::string &scriptname() { return m_scriptname; }; diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index 0c75cf0fdc5..10d1ad6f4cf 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -72,7 +72,7 @@ nptiso_struct nptiso = {0.0, /************************************************************/ /** Calculate long range virials (P3M, MMM2d...). */ -void calc_long_range_virials(); +void calc_long_range_virials(const ParticleRange &particles); /** Initializes a virials Observable stat. */ void init_virials(Observable_stat *stat); @@ -138,7 +138,7 @@ void pressure_calc(double *result, double *result_t, double *result_nb, /* rescale kinetic energy (=ideal contribution) */ virials.data.e[0] /= (3.0 * volume * time_step * time_step); - calc_long_range_virials(); + calc_long_range_virials(local_cells.particles()); #ifdef VIRTUAL_SITES virtual_sites()->pressure_and_stress_tensor_contribution( @@ -175,10 +175,10 @@ void pressure_calc(double *result, double *result_t, double *result_nb, /************************************************************/ -void calc_long_range_virials() { +void calc_long_range_virials(const ParticleRange &particles) { #ifdef ELECTROSTATICS /* calculate k-space part of electrostatic interaction. */ - Coulomb::calc_pressure_long_range(virials, p_tensor); + Coulomb::calc_pressure_long_range(virials, p_tensor, particles); #endif /*ifdef ELECTROSTATICS */ #ifdef DIPOLES diff --git a/src/core/virtual_sites/virtual_sites_com.cpp b/src/core/virtual_sites/virtual_sites_com.cpp index ade86a14951..9e36f20e583 100644 --- a/src/core/virtual_sites/virtual_sites_com.cpp +++ b/src/core/virtual_sites/virtual_sites_com.cpp @@ -46,8 +46,8 @@ void update_mol_pos_particle(Particle *p) { } } -void distribute_mol_force() { - for (auto &p : local_cells.particles()) { +void distribute_mol_force(const ParticleRange &particles) { + for (auto &p : particles) { if (p.p.is_virtual) { if (p.f.f.norm2() != 0) { put_mol_force_on_parts(&p); diff --git a/src/core/virtual_sites/virtual_sites_com.hpp b/src/core/virtual_sites/virtual_sites_com.hpp index 0f1f41add8d..05b4c8bfb3c 100644 --- a/src/core/virtual_sites/virtual_sites_com.hpp +++ b/src/core/virtual_sites/virtual_sites_com.hpp @@ -35,10 +35,10 @@ void update_mol_vel_particle(Particle *); // Distribute forces that have accumulated on virtual particles to the // associated real particles -void distribute_mol_force(); +void distribute_mol_force(const ParticleRange &particles); -// Gets the (first) virtual particle of the same molecule as the given (real) -// particle +// Gets the (first) virtual particle of the same molecule as the given +// (real) particle Particle *get_mol_com_particle(Particle *calling_p); #endif diff --git a/src/script_interface/h5md/h5md.cpp b/src/script_interface/h5md/h5md.cpp index 51ac5616640..1a16e864314 100644 --- a/src/script_interface/h5md/h5md.cpp +++ b/src/script_interface/h5md/h5md.cpp @@ -34,7 +34,7 @@ Variant H5mdScript::call_method(const std::string &name, if (name == "init_file") m_h5md->InitFile(); else if (name == "write") - m_h5md->Write(m_h5md->what(), partCfg()); + m_h5md->Write(m_h5md->what(), partCfg(), local_cells.particles()); else if (name == "flush") m_h5md->Flush(); else if (name == "close")