From 3c418b81bd972a7d3d57c8e532fa34c055cff4c0 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Thu, 8 Aug 2019 09:11:40 +0200 Subject: [PATCH 01/13] replaced local_particles in p3m and icc --- .../electrostatics_magnetostatics/coulomb.cpp | 19 ++++++++------- .../electrostatics_magnetostatics/coulomb.hpp | 8 +++---- .../electrostatics_magnetostatics/icc.cpp | 20 ++++++++-------- .../electrostatics_magnetostatics/icc.hpp | 2 +- .../electrostatics_magnetostatics/p3m.cpp | 23 ++++++++++--------- .../electrostatics_magnetostatics/p3m.hpp | 3 ++- src/core/energy.cpp | 6 ++--- src/core/energy.hpp | 2 +- src/core/forces.cpp | 8 +++---- src/core/forces.hpp | 2 +- src/core/pressure.cpp | 8 +++---- 11 files changed, 51 insertions(+), 50 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/coulomb.cpp b/src/core/electrostatics_magnetostatics/coulomb.cpp index 59e81c6b0c6..e4cce0c4305 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.cpp +++ b/src/core/electrostatics_magnetostatics/coulomb.cpp @@ -40,8 +40,7 @@ void pressure_n(int &n_coulomb) { } } -void calc_pressure_long_range(Observable_stat &virials, - Observable_stat &p_tensor) { +void calc_pressure_long_range(Observable_stat &virials, Observable_stat &p_tensor, const ParticleRange &particles) { switch (coulomb.method) { #ifdef P3M case COULOMB_ELC_P3M: @@ -54,9 +53,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(); + p3m_charge_assign(particles); virials.coulomb[1] = p3m_calc_kspace_forces(0, 1); - p3m_charge_assign(); + p3m_charge_assign(particles); p3m_calc_kspace_stress(p_tensor.coulomb + 9); break; } @@ -287,7 +286,7 @@ void init() { } } -void calc_long_range_force() { +void calc_long_range_force(const ParticleRange &particles) { switch (coulomb.method) { #ifdef P3M case COULOMB_ELC_P3M: @@ -296,7 +295,7 @@ void calc_long_range_force() { ELC_p3m_charge_assign_both(); ELC_P3M_self_forces(); } else - p3m_charge_assign(); + p3m_charge_assign(particles); p3m_calc_kspace_forces(1, 0); @@ -320,7 +319,7 @@ 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); @@ -351,7 +350,7 @@ 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,13 +358,13 @@ void calc_energy_long_range(Observable_stat &energy) { << "long range energy calculation not implemented for GPU P3M"; break; case COULOMB_P3M: - p3m_charge_assign(); + p3m_charge_assign(particles); energy.coulomb[1] = p3m_calc_kspace_forces(0, 1); 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); else { diff --git a/src/core/electrostatics_magnetostatics/coulomb.hpp b/src/core/electrostatics_magnetostatics/coulomb.hpp index 8fa5ccb9795..e7b8d6a6bb3 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.hpp +++ b/src/core/electrostatics_magnetostatics/coulomb.hpp @@ -7,6 +7,7 @@ #include "Observable_stat.hpp" #include +#include /** \name Type codes for the type of Coulomb interaction Enumeration of implemented methods for the electrostatic @@ -50,8 +51,7 @@ extern Coulomb_parameters coulomb; namespace Coulomb { void pressure_n(int &n_coulomb); -void calc_pressure_long_range(Observable_stat &virials, - Observable_stat &p_tensor); +void calc_pressure_long_range(Observable_stat &virials, Observable_stat &p_tensor, const ParticleRange &particles); void sanity_checks(int &state); double cutoff(const Utils::Vector3d &box_l); @@ -64,9 +64,9 @@ void on_resort_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/icc.cpp b/src/core/electrostatics_magnetostatics/icc.cpp index ace28c81ad8..dc8f3892014 100644 --- a/src/core/electrostatics_magnetostatics/icc.cpp +++ b/src/core/electrostatics_magnetostatics/icc.cpp @@ -55,7 +55,7 @@ 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 +63,7 @@ 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 +86,7 @@ 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,7 +105,7 @@ 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 + force_calc_iccp3m(particles, ghost_particles); /* Calculate electrostatic forces (SR+LR) excluding source source interaction*/ ghost_communicator(&cell_structure.collect_ghost_force_comm); @@ -188,8 +188,8 @@ 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 +197,15 @@ 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..1a2fb4dd593 100644 --- a/src/core/electrostatics_magnetostatics/icc.hpp +++ b/src/core/electrostatics_magnetostatics/icc.hpp @@ -95,7 +95,7 @@ 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/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp index e6b581dd2dd..f110c2e77a2 100644 --- a/src/core/electrostatics_magnetostatics/p3m.cpp +++ b/src/core/electrostatics_magnetostatics/p3m.cpp @@ -251,7 +251,7 @@ 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 +493,42 @@ 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++; diff --git a/src/core/electrostatics_magnetostatics/p3m.hpp b/src/core/electrostatics_magnetostatics/p3m.hpp index de1c993ba4d..069299376c3 100644 --- a/src/core/electrostatics_magnetostatics/p3m.hpp +++ b/src/core/electrostatics_magnetostatics/p3m.hpp @@ -61,6 +61,7 @@ #include #include +#include /************************************************ * data types @@ -200,7 +201,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..97aa02170af 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,10 +123,10 @@ 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 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/forces.cpp b/src/core/forces.cpp index a030809ed42..ef522bde8a3 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) { @@ -177,11 +177,11 @@ 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 */ 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/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 From 39f6d6c9443f4f7c2766ec74ac7897dd9c1c1068 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Thu, 8 Aug 2019 09:51:52 +0200 Subject: [PATCH 02/13] replaced local_particles in elc and mmm2d --- src/core/cells.cpp | 2 +- .../electrostatics_magnetostatics/coulomb.cpp | 32 ++-- .../electrostatics_magnetostatics/coulomb.hpp | 2 +- .../electrostatics_magnetostatics/elc.cpp | 163 +++++++++--------- .../electrostatics_magnetostatics/elc.hpp | 18 +- .../electrostatics_magnetostatics/mmm2d.cpp | 36 ++-- .../electrostatics_magnetostatics/mmm2d.hpp | 9 +- src/core/event.cpp | 4 +- src/core/event.hpp | 3 +- 9 files changed, 137 insertions(+), 132 deletions(-) diff --git a/src/core/cells.cpp b/src/core/cells.cpp index c120a8ae6a8..72a6b86e509 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -402,7 +402,7 @@ void cells_resort_particles(int global_flag) { resort_particles = Cells::RESORT_NONE; rebuild_verletlist = 1; - 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 e4cce0c4305..187488607c6 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.cpp +++ b/src/core/electrostatics_magnetostatics/coulomb.cpp @@ -219,7 +219,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: @@ -227,7 +227,7 @@ void on_resort_particles() { break; #endif case COULOMB_MMM2D: - MMM2D_on_resort_particles(); + MMM2D_on_resort_particles(particles); break; default: break; @@ -291,18 +291,18 @@ void calc_long_range_force(const ParticleRange &particles) { #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(particles); p3m_calc_kspace_forces(1, 0); 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 @@ -329,7 +329,7 @@ void calc_long_range_force(const ParticleRange &particles) { break; #endif case COULOMB_MMM2D: - MMM2D_add_far_force(); + MMM2D_add_far_force(particles); MMM2D_dielectric_layers_force_contribution(); break; #ifdef SCAFACOS @@ -369,24 +369,24 @@ void calc_energy_long_range(Observable_stat &energy, const ParticleRange &partic energy.coulomb[1] = p3m_calc_kspace_forces(0, 1); 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 * 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); // 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); // 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 @@ -396,7 +396,7 @@ void calc_energy_long_range(Observable_stat &energy, const ParticleRange &partic 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 e7b8d6a6bb3..68fc7372088 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.hpp +++ b/src/core/electrostatics_magnetostatics/coulomb.hpp @@ -60,7 +60,7 @@ 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(); diff --git a/src/core/electrostatics_magnetostatics/elc.cpp b/src/core/electrostatics_magnetostatics/elc.cpp index 1bc8fd3e62a..465ce008c75 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,25 @@ 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 +157,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 +165,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 +174,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 +182,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 +269,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 +279,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 +303,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 +337,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 +347,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 +364,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 +439,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 +451,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 +471,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 +517,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 +539,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 +562,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 +572,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 +595,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 +679,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 +701,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 +786,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 +819,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 +855,7 @@ 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 +877,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 +984,7 @@ 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 +992,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 +1043,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 +1083,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 +1124,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 +1319,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 +1348,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 +1356,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 +1384,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 +1514,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 +1546,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 +1554,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 +1582,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 +1590,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 +1616,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 +1624,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/mmm2d.cpp b/src/core/electrostatics_magnetostatics/mmm2d.cpp index 604e448413d..9f96ad4a058 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,7 @@ 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 +1130,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 +1170,7 @@ 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 +1179,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 +1219,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 +1262,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 +1277,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 +1699,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 +1707,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 +1835,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 +1848,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..390b00e626c 100644 --- a/src/core/electrostatics_magnetostatics/mmm2d.hpp +++ b/src/core/electrostatics_magnetostatics/mmm2d.hpp @@ -34,6 +34,7 @@ #include "config.hpp" #include +#include #ifdef ELECTROSTATICS @@ -88,13 +89,13 @@ 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 +113,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/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(); From c50f1ef359ac10f49eb021bee0917c4acec1b52d Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Thu, 8 Aug 2019 10:08:35 +0200 Subject: [PATCH 03/13] replaced local_particles in p3m --- .../electrostatics_magnetostatics/coulomb.cpp | 18 +++++------ .../electrostatics_magnetostatics/p3m.cpp | 30 +++++++++---------- .../electrostatics_magnetostatics/p3m.hpp | 2 +- 3 files changed, 25 insertions(+), 25 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/coulomb.cpp b/src/core/electrostatics_magnetostatics/coulomb.cpp index 187488607c6..1b7a6e6385b 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.cpp +++ b/src/core/electrostatics_magnetostatics/coulomb.cpp @@ -54,7 +54,7 @@ void calc_pressure_long_range(Observable_stat &virials, Observable_stat &p_tenso break; case COULOMB_P3M: { p3m_charge_assign(particles); - virials.coulomb[1] = p3m_calc_kspace_forces(0, 1); + virials.coulomb[1] = p3m_calc_kspace_forces(0, 1, particles); p3m_charge_assign(particles); p3m_calc_kspace_stress(p_tensor.coulomb + 9); break; @@ -297,7 +297,7 @@ void calc_long_range_force(const ParticleRange &particles) { } else 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(particles); @@ -322,10 +322,10 @@ void calc_long_range_force(const ParticleRange &particles) { 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: @@ -359,29 +359,29 @@ void calc_energy_long_range(Observable_stat &energy, const ParticleRange &partic break; case COULOMB_P3M: p3m_charge_assign(particles); - energy.coulomb[1] = p3m_calc_kspace_forces(0, 1); + 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(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 * 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(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(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(particles); diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp index f110c2e77a2..a5e47bef968 100644 --- a/src/core/electrostatics_magnetostatics/p3m.cpp +++ b/src/core/electrostatics_magnetostatics/p3m.cpp @@ -132,7 +132,7 @@ 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 @@ -674,7 +674,7 @@ 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 @@ -688,7 +688,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 @@ -768,7 +768,7 @@ 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 */ @@ -881,32 +881,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; @@ -914,7 +914,7 @@ 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); @@ -924,7 +924,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]); @@ -941,7 +941,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 069299376c3..f3345347ae7 100644 --- a/src/core/electrostatics_magnetostatics/p3m.hpp +++ b/src/core/electrostatics_magnetostatics/p3m.hpp @@ -183,7 +183,7 @@ 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); From c9f34aaeea7d67930412d653ca05a9951746435a Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Thu, 8 Aug 2019 10:23:06 +0200 Subject: [PATCH 04/13] replaced local_particles in dp3m --- .../electrostatics_magnetostatics/dipole.cpp | 18 +++++------ .../electrostatics_magnetostatics/dipole.hpp | 5 +-- .../p3m-dipolar.cpp | 32 +++++++++---------- .../p3m-dipolar.hpp | 5 +-- src/core/energy.cpp | 2 +- src/core/forces.cpp | 2 +- 6 files changed, 33 insertions(+), 31 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/dipole.cpp b/src/core/electrostatics_magnetostatics/dipole.cpp index f802536ce57..8fb3755a041 100644 --- a/src/core/electrostatics_magnetostatics/dipole.cpp +++ b/src/core/electrostatics_magnetostatics/dipole.cpp @@ -155,21 +155,21 @@ 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(); // 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 @@ -205,16 +205,16 @@ 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); + dp3m_dipole_assign(particles); + energy.dipolar[1] = dp3m_calc_kspace_forces(0, 1, particles); energy.dipolar[2] = add_mdlc_energy_corrections(); break; #endif diff --git a/src/core/electrostatics_magnetostatics/dipole.hpp b/src/core/electrostatics_magnetostatics/dipole.hpp index d21fb88f9c7..9009367d240 100644 --- a/src/core/electrostatics_magnetostatics/dipole.hpp +++ b/src/core/electrostatics_magnetostatics/dipole.hpp @@ -11,6 +11,7 @@ extern double dipolar_cutoff; #include #include +#include /** \name Compounds for Dipole interactions */ /*@{*/ @@ -69,9 +70,9 @@ 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/p3m-dipolar.cpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp index 40b16b92107..1d642589c62 100644 --- a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp +++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp @@ -181,7 +181,7 @@ 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 +636,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 +645,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 +761,7 @@ 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 +773,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 +825,7 @@ 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 +837,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 +866,7 @@ 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 */ @@ -1043,7 +1043,7 @@ double dp3m_calc_kspace_forces(int force_flag, int energy_flag) { 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); + 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 +1129,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 +1141,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 +1151,14 @@ 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 +1167,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 +1214,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..ca598f1f6e5 100644 --- a/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp +++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp @@ -46,6 +46,7 @@ #include #include +#include struct dp3m_data_struct { dp3m_data_struct(); @@ -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,7 @@ 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/energy.cpp b/src/core/energy.cpp index 97aa02170af..8ade508d33d 100644 --- a/src/core/energy.cpp +++ b/src/core/energy.cpp @@ -130,7 +130,7 @@ void calc_long_range_energies(const ParticleRange &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/forces.cpp b/src/core/forces.cpp index ef522bde8a3..62ea3141a86 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -187,6 +187,6 @@ void calc_long_range_forces(const ParticleRange &particles) { #ifdef DIPOLES /* calculate k-space part of the magnetostatic interaction. */ - Dipole::calc_long_range_force(); + Dipole::calc_long_range_force(particles); #endif /*ifdef DIPOLES */ } From 95b3391ee48f022759705172ef03110fa6a1f55b Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Thu, 8 Aug 2019 10:29:51 +0200 Subject: [PATCH 05/13] replaced local_particles in mdlc_correction --- .../electrostatics_magnetostatics/dipole.cpp | 8 ++-- .../mdlc_correction.cpp | 39 ++++++++++--------- .../mdlc_correction.hpp | 5 ++- 3 files changed, 27 insertions(+), 25 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/dipole.cpp b/src/core/electrostatics_magnetostatics/dipole.cpp index 8fb3755a041..d5ec6394a19 100644 --- a/src/core/electrostatics_magnetostatics/dipole.cpp +++ b/src/core/electrostatics_magnetostatics/dipole.cpp @@ -159,7 +159,7 @@ 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(particles); @@ -178,7 +178,7 @@ void calc_long_range_force(const ParticleRange &particles) { break; #ifdef DP3M case DIPOLAR_MDLC_DS: - add_mdlc_force_corrections(); + add_mdlc_force_corrections(particles); // fall through #endif case DIPOLAR_DS: @@ -215,7 +215,7 @@ void calc_energy_long_range(Observable_stat &energy, const ParticleRange &partic case DIPOLAR_MDLC_P3M: dp3m_dipole_assign(particles); energy.dipolar[1] = dp3m_calc_kspace_forces(0, 1, particles); - energy.dipolar[2] = add_mdlc_energy_corrections(); + energy.dipolar[2] = add_mdlc_energy_corrections(particles; break; #endif case DIPOLAR_ALL_WITH_ALL_AND_NO_REPLICA: @@ -224,7 +224,7 @@ void calc_energy_long_range(Observable_stat &energy, const ParticleRange &partic #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[2] = add_mdlc_energy_corrections(particles); break; #endif case DIPOLAR_DS: diff --git a/src/core/electrostatics_magnetostatics/mdlc_correction.cpp b/src/core/electrostatics_magnetostatics/mdlc_correction.cpp index 2f59fc84e30..4a7499007aa 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,11 @@ 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(); } @@ -106,8 +107,8 @@ double slab_dip_count_mu(double *mt, double *mx, double *my) { Algorithm implemented accordingly to the paper of A. Brodka, Chem. Phys. Lett. 400, 62-67, (2004). */ -double get_DLC_dipolar(int kcut, std::vector &fs, - std::vector &ts) { +double get_DLC_dipolar(int kcut, std::vector &fs, std::vector &ts, + const ParticleRange &particles) { int ip; @@ -142,7 +143,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 +183,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 +235,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 +267,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 +277,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 +303,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 +347,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 +364,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 +375,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 +410,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 +422,7 @@ 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 +436,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..8ffc30c097b 100644 --- a/src/core/electrostatics_magnetostatics/mdlc_correction.hpp +++ b/src/core/electrostatics_magnetostatics/mdlc_correction.hpp @@ -39,6 +39,7 @@ #ifndef _DLC_DIPOLAR_H #define _DLC_DIPOLAR_H +#include #include "config.hpp" #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 From 42f7828b1db276a996edf7dc51379346d801928b Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Thu, 8 Aug 2019 10:33:29 +0200 Subject: [PATCH 06/13] extracted local_cells in magnetic_non_p3m_methods --- src/core/electrostatics_magnetostatics/dipole.cpp | 12 ++++++------ .../magnetic_non_p3m_methods.cpp | 11 +++++------ .../magnetic_non_p3m_methods.hpp | 6 +++--- 3 files changed, 14 insertions(+), 15 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/dipole.cpp b/src/core/electrostatics_magnetostatics/dipole.cpp index d5ec6394a19..6d0a4c30201 100644 --- a/src/core/electrostatics_magnetostatics/dipole.cpp +++ b/src/core/electrostatics_magnetostatics/dipole.cpp @@ -174,7 +174,7 @@ void calc_long_range_force(const ParticleRange &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: @@ -182,7 +182,7 @@ void calc_long_range_force(const ParticleRange &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 @@ -215,20 +215,20 @@ void calc_energy_long_range(Observable_stat &energy, const ParticleRange &partic case DIPOLAR_MDLC_P3M: dp3m_dipole_assign(particles); energy.dipolar[1] = dp3m_calc_kspace_forces(0, 1, particles); - energy.dipolar[2] = add_mdlc_energy_corrections(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[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/magnetic_non_p3m_methods.cpp b/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.cpp index 5f617e3b8c8..76cdcc090a9 100644 --- a/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.cpp +++ b/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.cpp @@ -127,7 +127,7 @@ 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 +142,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 +184,7 @@ 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 +224,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 +344,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..dffbe8a8941 100644 --- a/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.hpp +++ b/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.hpp @@ -32,6 +32,7 @@ * sum, * */ +#include #include "config.hpp" #ifdef DIPOLES @@ -47,7 +48,7 @@ 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 +65,7 @@ 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 From 4c4703d98efc3cd24822bf68937acf9d307a6904 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Thu, 8 Aug 2019 10:36:48 +0200 Subject: [PATCH 07/13] replaced local_particles in virtual_sites_com (unused function) --- src/core/virtual_sites/virtual_sites_com.cpp | 4 ++-- src/core/virtual_sites/virtual_sites_com.hpp | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) 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..47252204d27 100644 --- a/src/core/virtual_sites/virtual_sites_com.hpp +++ b/src/core/virtual_sites/virtual_sites_com.hpp @@ -35,7 +35,7 @@ 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 From 84e05fc990cbe16a6d567bc59a7621076e0f6840 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Thu, 8 Aug 2019 10:47:33 +0200 Subject: [PATCH 08/13] replaced local_particles in lb_particle_coupling --- src/core/forces.cpp | 2 +- src/core/grid_based_algorithms/lb_particle_coupling.cpp | 4 ++-- src/core/grid_based_algorithms/lb_particle_coupling.hpp | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core/forces.cpp b/src/core/forces.cpp index 62ea3141a86..f83d177cd59 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -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 */ diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.cpp b/src/core/grid_based_algorithms/lb_particle_coupling.cpp index a8202225a8c..40090334d10 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.cpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.cpp @@ -182,7 +182,7 @@ 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 +244,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..8444785662d 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.hpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.hpp @@ -10,7 +10,7 @@ * 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); From b9928c34c71f51016a4bdc8ec8f20937f0cddf1a Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Thu, 8 Aug 2019 11:01:18 +0200 Subject: [PATCH 09/13] fixed missing header in lb_particle_coupling --- src/core/grid_based_algorithms/lb_particle_coupling.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.hpp b/src/core/grid_based_algorithms/lb_particle_coupling.hpp index 8444785662d..26c805aac0a 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.hpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.hpp @@ -4,6 +4,7 @@ #include #include +#include "ParticleRange.hpp" /** Calculate particle lattice interactions. * So far, only viscous coupling with Stokesian friction is implemented. From 51af0925d722453091a901840d79dcf047f0a23f Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Thu, 8 Aug 2019 11:44:33 +0200 Subject: [PATCH 10/13] local_cells in h5md --- src/core/io/writer/h5md_core.cpp | 4 ++-- src/core/io/writer/h5md_core.hpp | 2 +- src/script_interface/h5md/h5md.cpp | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core/io/writer/h5md_core.cpp b/src/core/io/writer/h5md_core.cpp index 30d88712555..12b5c81c57a 100644 --- a/src/core/io/writer/h5md_core.cpp +++ b/src/core/io/writer/h5md_core.cpp @@ -392,7 +392,7 @@ 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 +441,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/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") From 4ad8144dc389b4f2d467e033c0e35a6ba67af5da Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Thu, 8 Aug 2019 11:52:30 +0200 Subject: [PATCH 11/13] formatting --- src/config/myconfig-default.hpp | 2 +- .../electrostatics_magnetostatics/coulomb.cpp | 12 ++++++---- .../electrostatics_magnetostatics/coulomb.hpp | 9 ++++--- .../electrostatics_magnetostatics/dipole.cpp | 11 +++++---- .../electrostatics_magnetostatics/dipole.hpp | 5 ++-- .../electrostatics_magnetostatics/elc.cpp | 12 ++++++---- .../electrostatics_magnetostatics/icc.cpp | 19 +++++++++------ .../electrostatics_magnetostatics/icc.hpp | 3 ++- .../magnetic_non_p3m_methods.cpp | 7 ++++-- .../magnetic_non_p3m_methods.hpp | 8 ++++--- .../mdlc_correction.cpp | 11 +++++---- .../mdlc_correction.hpp | 2 +- .../electrostatics_magnetostatics/mmm2d.cpp | 6 +++-- .../electrostatics_magnetostatics/mmm2d.hpp | 10 +++++--- .../p3m-dipolar.cpp | 24 ++++++++++++------- .../p3m-dipolar.hpp | 5 ++-- .../electrostatics_magnetostatics/p3m.cpp | 18 ++++++++------ .../electrostatics_magnetostatics/p3m.hpp | 5 ++-- .../lb_particle_coupling.cpp | 3 ++- .../lb_particle_coupling.hpp | 5 ++-- src/core/io/writer/h5md_core.cpp | 3 ++- src/core/virtual_sites/virtual_sites_com.hpp | 6 ++--- 22 files changed, 118 insertions(+), 68 deletions(-) diff --git a/src/config/myconfig-default.hpp b/src/config/myconfig-default.hpp index 2626c147b92..f68c22c1694 100644 --- a/src/config/myconfig-default.hpp +++ b/src/config/myconfig-default.hpp @@ -84,4 +84,4 @@ // Further featuers #define VIRTUAL_SITES_RELATIVE #define VIRTUAL_SITES_INERTIALESS_TRACERS -#define COLLISION_DETECTION +#define COLLISION_DETECTION \ No newline at end of file diff --git a/src/core/electrostatics_magnetostatics/coulomb.cpp b/src/core/electrostatics_magnetostatics/coulomb.cpp index 1b7a6e6385b..bf71440af35 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.cpp +++ b/src/core/electrostatics_magnetostatics/coulomb.cpp @@ -40,7 +40,9 @@ void pressure_n(int &n_coulomb) { } } -void calc_pressure_long_range(Observable_stat &virials, Observable_stat &p_tensor, const ParticleRange &particles) { +void calc_pressure_long_range(Observable_stat &virials, + Observable_stat &p_tensor, + const ParticleRange &particles) { switch (coulomb.method) { #ifdef P3M case COULOMB_ELC_P3M: @@ -325,7 +327,7 @@ void calc_long_range_force(const ParticleRange &particles) { nptiso.p_vir[0] += p3m_calc_kspace_forces(1, 1, particles); else #endif - p3m_calc_kspace_forces(1, 0, particles); + p3m_calc_kspace_forces(1, 0, particles); break; #endif case COULOMB_MMM2D: @@ -350,7 +352,8 @@ void calc_long_range_force(const ParticleRange &particles) { #endif } -void calc_energy_long_range(Observable_stat &energy, const ParticleRange &particles) { +void calc_energy_long_range(Observable_stat &energy, + const ParticleRange &particles) { switch (coulomb.method) { #ifdef P3M case COULOMB_P3M_GPU: @@ -369,7 +372,8 @@ void calc_energy_long_range(Observable_stat &energy, const ParticleRange &partic energy.coulomb[1] = p3m_calc_kspace_forces(0, 1, particles); else { energy.coulomb[1] = 0.5 * p3m_calc_kspace_forces(0, 1, particles); - energy.coulomb[1] += 0.5 * ELC_P3M_dielectric_layers_energy_self(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(particles); diff --git a/src/core/electrostatics_magnetostatics/coulomb.hpp b/src/core/electrostatics_magnetostatics/coulomb.hpp index 68fc7372088..6166573fc99 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.hpp +++ b/src/core/electrostatics_magnetostatics/coulomb.hpp @@ -6,8 +6,8 @@ #ifdef ELECTROSTATICS #include "Observable_stat.hpp" -#include #include +#include /** \name Type codes for the type of Coulomb interaction Enumeration of implemented methods for the electrostatic @@ -51,7 +51,9 @@ extern Coulomb_parameters coulomb; namespace Coulomb { void pressure_n(int &n_coulomb); -void calc_pressure_long_range(Observable_stat &virials, Observable_stat &p_tensor, const ParticleRange &particles); +void calc_pressure_long_range(Observable_stat &virials, + Observable_stat &p_tensor, + const ParticleRange &particles); void sanity_checks(int &state); double cutoff(const Utils::Vector3d &box_l); @@ -66,7 +68,8 @@ void init(); void calc_long_range_force(const ParticleRange &particles); -void calc_energy_long_range(Observable_stat &energy, const ParticleRange &particles); +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 6d0a4c30201..11d59b2937d 100644 --- a/src/core/electrostatics_magnetostatics/dipole.cpp +++ b/src/core/electrostatics_magnetostatics/dipole.cpp @@ -169,7 +169,7 @@ void calc_long_range_force(const ParticleRange &particles) { fprintf(stderr, "dipolar_P3M at this moment is added to p_vir[0]\n"); } else #endif - dp3m_calc_kspace_forces(1, 0, particles); + dp3m_calc_kspace_forces(1, 0, particles); break; #endif @@ -205,7 +205,8 @@ void calc_long_range_force(const ParticleRange &particles) { } } -void calc_energy_long_range(Observable_stat &energy, const ParticleRange &particles) { +void calc_energy_long_range(Observable_stat &energy, + const ParticleRange &particles) { switch (dipole.method) { #ifdef DP3M case DIPOLAR_P3M: @@ -223,12 +224,14 @@ void calc_energy_long_range(Observable_stat &energy, const ParticleRange &partic break; #ifdef DP3M case DIPOLAR_MDLC_DS: - energy.dipolar[1] = magnetic_dipolar_direct_sum_calculations(0, 1, particles); + 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, particles); + 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 9009367d240..4883919fb9b 100644 --- a/src/core/electrostatics_magnetostatics/dipole.hpp +++ b/src/core/electrostatics_magnetostatics/dipole.hpp @@ -10,8 +10,8 @@ extern double dipolar_cutoff; #include -#include #include +#include /** \name Compounds for Dipole interactions */ /*@{*/ @@ -72,7 +72,8 @@ void init(); void calc_long_range_force(const ParticleRange &particles); -void calc_energy_long_range(Observable_stat &energy, const ParticleRange &particles); +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 465ce008c75..b5086eadf8f 100644 --- a/src/core/electrostatics_magnetostatics/elc.cpp +++ b/src/core/electrostatics_magnetostatics/elc.cpp @@ -133,8 +133,10 @@ static double Q_energy(double omega); /*@}*/ /** \name p,q <> 0 per frequency code */ /*@{*/ -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 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(const ParticleRange &particles); static double dipole_energy(const ParticleRange &particles); @@ -855,7 +857,8 @@ static double Q_energy(double omega) { /* PQ particle blocks */ /*****************************************************************/ -static void setup_PQ(int p, int q, double omega, const ParticleRange &particles) { +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])); @@ -984,7 +987,8 @@ 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 void add_PQ_force(int p, int q, double omega, + const ParticleRange &particles) { int ic; double pref_x = C_2PI * ux * p / omega; diff --git a/src/core/electrostatics_magnetostatics/icc.cpp b/src/core/electrostatics_magnetostatics/icc.cpp index dc8f3892014..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(const ParticleRange &particles, const ParticleRange &ghosts_particles); +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(const ParticleRange &particles, const ParticleRange &gho * 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(const ParticleRange &particles, const ParticleRange &ghost_particles); +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, const ParticleRange &ghost_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, const ParticleRange &ghost_ for (int j = 0; j < iccp3m_cfg.num_iteration; j++) { double hmax = 0.; - force_calc_iccp3m(particles, ghost_particles); /* 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,7 +191,8 @@ int iccp3m_iteration(const ParticleRange &particles, const ParticleRange &ghost_ return iccp3m_cfg.citeration; } -void force_calc_iccp3m(const ParticleRange &particles, const ParticleRange &ghost_particles) { +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) { @@ -200,7 +204,8 @@ void force_calc_iccp3m(const ParticleRange &particles, const ParticleRange &ghos Coulomb::calc_long_range_force(particles); } -void init_forces_iccp3m(const ParticleRange &particles, const ParticleRange &ghosts_particles) { +void init_forces_iccp3m(const ParticleRange &particles, + const ParticleRange &ghosts_particles) { for (auto &p : particles) { p.f = ParticleForce{}; } diff --git a/src/core/electrostatics_magnetostatics/icc.hpp b/src/core/electrostatics_magnetostatics/icc.hpp index 1a2fb4dd593..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, const ParticleRange &ghost_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 76cdcc090a9..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, const ParticleRange &particles) { +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"); @@ -184,7 +185,9 @@ int magnetic_dipolar_direct_sum_sanity_checks() { /************************************************************/ -double magnetic_dipolar_direct_sum_calculations(int force_flag, int energy_flag, const ParticleRange &particles) { +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; diff --git a/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.hpp b/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.hpp index dffbe8a8941..eafab7a1774 100644 --- a/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.hpp +++ b/src/core/electrostatics_magnetostatics/magnetic_non_p3m_methods.hpp @@ -32,8 +32,8 @@ * sum, * */ -#include #include "config.hpp" +#include #ifdef DIPOLES #include "particle_data.hpp" @@ -48,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, const ParticleRange &particles); +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 @@ -65,7 +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, const ParticleRange &particles); +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 4a7499007aa..7612daf4f54 100644 --- a/src/core/electrostatics_magnetostatics/mdlc_correction.cpp +++ b/src/core/electrostatics_magnetostatics/mdlc_correction.cpp @@ -54,7 +54,7 @@ static double mu_max; void calc_mu_max() { auto local_particles = local_cells.particles(); mu_max = std::accumulate( - local_particles.begin(), local_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); @@ -77,7 +77,8 @@ 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, const ParticleRange &particles) { +double slab_dip_count_mu(double *mt, double *mx, double *my, + const ParticleRange &particles) { Utils::Vector3d node_sums{}; Utils::Vector3d tot_sums{}; @@ -107,7 +108,8 @@ double slab_dip_count_mu(double *mt, double *mx, double *my, const ParticleRange Algorithm implemented accordingly to the paper of A. Brodka, Chem. Phys. Lett. 400, 62-67, (2004). */ -double get_DLC_dipolar(int kcut, std::vector &fs, std::vector &ts, +double get_DLC_dipolar(int kcut, std::vector &fs, + std::vector &ts, const ParticleRange &particles) { int ip; @@ -422,7 +424,8 @@ double add_mdlc_energy_corrections(const ParticleRange &particles) { //---- Compute the corrections ---------------------------------- // First the DLC correction - dip_DLC_energy += dipole.prefactor * get_DLC_energy_dipolar(dip_DLC_kcut, particles); + dip_DLC_energy += + dipole.prefactor * get_DLC_energy_dipolar(dip_DLC_kcut, particles); // printf("Energy DLC = %20.15le // \n",dip_DLC_energy); diff --git a/src/core/electrostatics_magnetostatics/mdlc_correction.hpp b/src/core/electrostatics_magnetostatics/mdlc_correction.hpp index 8ffc30c097b..7dc37864aaf 100644 --- a/src/core/electrostatics_magnetostatics/mdlc_correction.hpp +++ b/src/core/electrostatics_magnetostatics/mdlc_correction.hpp @@ -39,8 +39,8 @@ #ifndef _DLC_DIPOLAR_H #define _DLC_DIPOLAR_H -#include #include "config.hpp" +#include #if defined(DIPOLES) && defined(DP3M) diff --git a/src/core/electrostatics_magnetostatics/mmm2d.cpp b/src/core/electrostatics_magnetostatics/mmm2d.cpp index 9f96ad4a058..9074ea48050 100644 --- a/src/core/electrostatics_magnetostatics/mmm2d.cpp +++ b/src/core/electrostatics_magnetostatics/mmm2d.cpp @@ -1118,7 +1118,8 @@ static double PQ_energy(double omega) { /* main loops */ /*****************************************************************/ -static void add_force_contribution(int p, int q, const ParticleRange &particles) { +static void add_force_contribution(int p, int q, + const ParticleRange &particles) { double omega, fac; if (q == 0) { @@ -1170,7 +1171,8 @@ static void add_force_contribution(int p, int q, const ParticleRange &particles) } } -static double energy_contribution(int p, int q, const ParticleRange &particles) { +static double energy_contribution(int p, int q, + const ParticleRange &particles) { double eng; double omega, fac; diff --git a/src/core/electrostatics_magnetostatics/mmm2d.hpp b/src/core/electrostatics_magnetostatics/mmm2d.hpp index 390b00e626c..015f57ef8d5 100644 --- a/src/core/electrostatics_magnetostatics/mmm2d.hpp +++ b/src/core/electrostatics_magnetostatics/mmm2d.hpp @@ -33,8 +33,8 @@ #include "config.hpp" -#include #include +#include #ifdef ELECTROSTATICS @@ -92,10 +92,14 @@ int MMM2D_set_params(double maxPWerror, double far_cut, double delta_top, double MMM2D_add_far(int f, int e, const ParticleRange &particles); /** the actual long range force calculation */ -inline void MMM2D_add_far_force(const ParticleRange &particles) { MMM2D_add_far(1, 0, particles); } +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(const ParticleRange &particles) { return MMM2D_add_far(0, 1, particles); } +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, diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp index 1d642589c62..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, const ParticleRange &particles); +static double calc_surface_term(int force_flag, int energy_flag, + const ParticleRange &particles); /** \name P3M Tuning Functions */ /************************************************************/ @@ -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, const ParticleRange &particles) { +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 */ @@ -825,7 +827,8 @@ static void P3M_assign_torques(double prefac, int d_rs, const ParticleRange &par #endif /* Assign the dipolar forces obtained from k-space */ -static void dp3m_assign_forces_dip(double prefac, int d_rs, const ParticleRange &particles) { +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 */ @@ -866,7 +869,8 @@ static void dp3m_assign_forces_dip(double prefac, int d_rs, const ParticleRange /*****************************************************************************/ -double dp3m_calc_kspace_forces(int force_flag, int energy_flag, const ParticleRange &particles) { +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, const ParticleRa /* 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, particles); + 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, const ParticleRa 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, particles); + dipole_prefac * pow(2 * Utils::pi() / box_geo.length()[0], 2), d_rs, + particles); } P3M_TRACE(fprintf(stderr, @@ -1151,7 +1156,8 @@ double dp3m_calc_kspace_forces(int force_flag, int energy_flag, const ParticleRa /************************************************************/ -double calc_surface_term(int force_flag, int energy_flag, const ParticleRange &particles) { +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); diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp index ca598f1f6e5..ac2211c209c 100644 --- a/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp +++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp @@ -44,9 +44,9 @@ #include "p3m-common.hpp" #include "particle_data.hpp" +#include #include #include -#include struct dp3m_data_struct { dp3m_data_struct(); @@ -194,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, const ParticleRange &particles); +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 a5e47bef968..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, const ParticleRange &particles); +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(const ParticleRange &particles); +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); @@ -520,8 +522,7 @@ void p3m_charge_assign(const ParticleRange &particles) { } /** Assign the charges */ -template -void p3m_do_charge_assign(const ParticleRange &particles) { +template void p3m_do_charge_assign(const ParticleRange &particles) { /* charged particle counter, charge fraction counter */ int cp_cnt = 0; /* prepare local FFT mesh */ @@ -674,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, const ParticleRange &particles) { +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 @@ -768,7 +770,8 @@ static void P3M_assign_forces(double force_prefac, int d_rs, const ParticleRange } } -double p3m_calc_kspace_forces(int force_flag, int energy_flag, const ParticleRange &particles) { +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 */ @@ -914,7 +917,8 @@ double p3m_calc_kspace_forces(int force_flag, int energy_flag, const ParticleRan /************************************************************/ -double p3m_calc_dipole_term(int force_flag, int energy_flag, const ParticleRange &particles) { +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); diff --git a/src/core/electrostatics_magnetostatics/p3m.hpp b/src/core/electrostatics_magnetostatics/p3m.hpp index f3345347ae7..6067092267b 100644 --- a/src/core/electrostatics_magnetostatics/p3m.hpp +++ b/src/core/electrostatics_magnetostatics/p3m.hpp @@ -59,9 +59,9 @@ #include "fft.hpp" #include "p3m-common.hpp" +#include #include #include -#include /************************************************ * data types @@ -183,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, const ParticleRange &particles); +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); diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.cpp b/src/core/grid_based_algorithms/lb_particle_coupling.cpp index 40090334d10..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, const ParticleRange &particles) { +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 diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.hpp b/src/core/grid_based_algorithms/lb_particle_coupling.hpp index 26c805aac0a..387ff734b60 100644 --- a/src/core/grid_based_algorithms/lb_particle_coupling.hpp +++ b/src/core/grid_based_algorithms/lb_particle_coupling.hpp @@ -3,15 +3,16 @@ #include -#include #include "ParticleRange.hpp" +#include /** Calculate particle lattice interactions. * So far, only viscous coupling with Stokesian friction is implemented. * 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, const ParticleRange &particles); +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 12b5c81c57a..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, const ParticleRange &particles) { +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 diff --git a/src/core/virtual_sites/virtual_sites_com.hpp b/src/core/virtual_sites/virtual_sites_com.hpp index 47252204d27..13b9e624e29 100644 --- a/src/core/virtual_sites/virtual_sites_com.hpp +++ b/src/core/virtual_sites/virtual_sites_com.hpp @@ -37,9 +37,9 @@ void update_mol_vel_particle(Particle *); // associated real particles void distribute_mol_force(const ParticleRange &particles) -// Gets the (first) virtual particle of the same molecule as the given (real) -// particle -Particle *get_mol_com_particle(Particle *calling_p); + // Gets the (first) virtual particle of the same molecule as the given + // (real) particle + Particle *get_mol_com_particle(Particle *calling_p); #endif From 33010b48c85f3b5024ce5edb6977995e940bd619 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Thu, 8 Aug 2019 11:57:28 +0200 Subject: [PATCH 12/13] reverted accidental newline changes --- src/config/myconfig-default.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/config/myconfig-default.hpp b/src/config/myconfig-default.hpp index f68c22c1694..2626c147b92 100644 --- a/src/config/myconfig-default.hpp +++ b/src/config/myconfig-default.hpp @@ -84,4 +84,4 @@ // Further featuers #define VIRTUAL_SITES_RELATIVE #define VIRTUAL_SITES_INERTIALESS_TRACERS -#define COLLISION_DETECTION \ No newline at end of file +#define COLLISION_DETECTION From a3ac1ac34e7e8c3eaba4a8e17bed0618db8996a7 Mon Sep 17 00:00:00 2001 From: Alexander Reinauer Date: Fri, 9 Aug 2019 07:08:40 +0200 Subject: [PATCH 13/13] missing semicolon --- src/core/virtual_sites/virtual_sites_com.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core/virtual_sites/virtual_sites_com.hpp b/src/core/virtual_sites/virtual_sites_com.hpp index 13b9e624e29..05b4c8bfb3c 100644 --- a/src/core/virtual_sites/virtual_sites_com.hpp +++ b/src/core/virtual_sites/virtual_sites_com.hpp @@ -35,11 +35,11 @@ void update_mol_vel_particle(Particle *); // Distribute forces that have accumulated on virtual particles to the // associated real particles -void distribute_mol_force(const ParticleRange &particles) +void distribute_mol_force(const ParticleRange &particles); - // Gets the (first) virtual particle of the same molecule as the given - // (real) particle - Particle *get_mol_com_particle(Particle *calling_p); +// Gets the (first) virtual particle of the same molecule as the given +// (real) particle +Particle *get_mol_com_particle(Particle *calling_p); #endif