From 486612d83ff0da57d7331f2456c07e4832066db8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 4 Sep 2020 18:50:05 +0200 Subject: [PATCH] core: Forward integrator_counter state by parameter Pass the current value of the integrator_counter to the force kernels as a parameter instead of relying on external linkage. There are two exceptions: thermalized bonds and constraints, which were not designed to be extensible in this way, and still rely on external linkage. --- src/core/constraints/ShapeBasedConstraint.cpp | 6 +- src/core/dpd.cpp | 12 ++-- src/core/dpd.hpp | 2 +- src/core/forces.cpp | 13 ++-- src/core/forces.hpp | 9 +-- src/core/forces_inline.hpp | 22 +++--- src/core/integrate.cpp | 14 ++-- src/core/integrators/brownian_inline.hpp | 34 ++++++---- src/core/integrators/langevin_inline.hpp | 15 +++-- .../integrators/stokesian_dynamics_inline.hpp | 11 +-- src/core/integrators/velocity_verlet_npt.cpp | 36 ++++++---- src/core/integrators/velocity_verlet_npt.hpp | 8 ++- src/core/stokesian_dynamics/sd_interface.cpp | 6 +- src/core/stokesian_dynamics/sd_interface.hpp | 2 +- src/core/thermostat.hpp | 15 +++-- src/core/unit_tests/thermostats_test.cpp | 67 +++++++++---------- 16 files changed, 145 insertions(+), 127 deletions(-) diff --git a/src/core/constraints/ShapeBasedConstraint.cpp b/src/core/constraints/ShapeBasedConstraint.cpp index 098f932ff9b..ca32c4bb063 100644 --- a/src/core/constraints/ShapeBasedConstraint.cpp +++ b/src/core/constraints/ShapeBasedConstraint.cpp @@ -74,8 +74,8 @@ ParticleForce ShapeBasedConstraint::force(Particle const &p, dist, &torque1, &torque2); #ifdef DPD if (thermo_switch & THERMO_DPD) { - force1 += - dpd_pair_force(p, part_rep, ia_params, dist_vec, dist, dist * dist); + force1 += dpd_pair_force(p, part_rep, ia_params, dist_vec, dist, + dist * dist, integrator_counter.value()); } #endif } else if (m_penetrable && (dist <= 0)) { @@ -85,7 +85,7 @@ ParticleForce ShapeBasedConstraint::force(Particle const &p, #ifdef DPD if (thermo_switch & THERMO_DPD) { force1 += dpd_pair_force(p, part_rep, ia_params, dist_vec, dist, - dist * dist); + dist * dist, integrator_counter.value()); } #endif } diff --git a/src/core/dpd.cpp b/src/core/dpd.cpp index 49debe8f78e..c11e9ede490 100644 --- a/src/core/dpd.cpp +++ b/src/core/dpd.cpp @@ -46,11 +46,11 @@ using Utils::Vector3d; * 3. Two particle IDs (order-independent, decorrelates particles, gets rid of * seed-per-node) */ -Vector3d dpd_noise(uint32_t pid1, uint32_t pid2) { +Vector3d dpd_noise(uint32_t pid1, uint32_t pid2, uint64_t counter) { extern DPDThermostat dpd; - return Random::noise_uniform( - integrator_counter.value(), dpd.rng_seed(), (pid1 < pid2) ? pid2 : pid1, - (pid1 < pid2) ? pid1 : pid2); + return Random::noise_uniform(counter, dpd.rng_seed(), + (pid1 < pid2) ? pid2 : pid1, + (pid1 < pid2) ? pid1 : pid2); } int dpd_set_params(int part_type_a, int part_type_b, double gamma, double k, @@ -117,7 +117,7 @@ Vector3d dpd_pair_force(DPDParameters const ¶ms, Vector3d const &v, Utils::Vector3d dpd_pair_force(Particle const &p1, Particle const &p2, IA_parameters const &ia_params, Utils::Vector3d const &d, double dist, - double dist2) { + double dist2, uint64_t counter) { if (ia_params.dpd_radial.cutoff <= 0.0 && ia_params.dpd_trans.cutoff <= 0.0) { return {}; } @@ -125,7 +125,7 @@ Utils::Vector3d dpd_pair_force(Particle const &p1, Particle const &p2, auto const v21 = p1.m.v - p2.m.v; auto const noise_vec = (ia_params.dpd_radial.pref > 0.0 || ia_params.dpd_trans.pref > 0.0) - ? dpd_noise(p1.p.identity, p2.p.identity) + ? dpd_noise(p1.p.identity, p2.p.identity, counter) : Vector3d{}; auto const f_r = dpd_pair_force(ia_params.dpd_radial, v21, dist, noise_vec); diff --git a/src/core/dpd.hpp b/src/core/dpd.hpp index 39f119c25d4..9b3048e0e6d 100644 --- a/src/core/dpd.hpp +++ b/src/core/dpd.hpp @@ -51,7 +51,7 @@ void dpd_update_params(double pref2_scale); Utils::Vector3d dpd_pair_force(Particle const &p1, Particle const &p2, IA_parameters const &ia_params, Utils::Vector3d const &d, double dist, - double dist2); + double dist2, uint64_t counter); Utils::Vector9d dpd_stress(); #endif diff --git a/src/core/forces.cpp b/src/core/forces.cpp index 4b5f4b34eeb..43fc5f9dd04 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -49,7 +49,7 @@ ActorList forceActors; -void init_forces(const ParticleRange &particles) { +void init_forces(const ParticleRange &particles, uint64_t counter) { ESPRESSO_PROFILER_CXX_MARK_FUNCTION; /* The force initialization depends on the used thermostat and the thermodynamic ensemble */ @@ -63,7 +63,7 @@ void init_forces(const ParticleRange &particles) { set torque to zero for all and rescale quaternions */ for (auto &p : particles) { - p.f = init_local_particle_force(p); + p.f = init_local_particle_force(p, counter); } /* initialize ghost forces with zero @@ -80,7 +80,7 @@ void init_forces_ghosts(const ParticleRange &particles) { } } -void force_calc(CellStructure &cell_structure) { +void force_calc(CellStructure &cell_structure, uint64_t const counter) { ESPRESSO_PROFILER_CXX_MARK_FUNCTION; espressoSystemInterface.update(); @@ -94,7 +94,7 @@ void force_calc(CellStructure &cell_structure) { #ifdef ELECTROSTATICS iccp3m_iteration(particles, cell_structure.ghost_particles()); #endif - init_forces(particles); + init_forces(particles, counter); for (auto &forceActor : forceActors) { forceActor->computeForces(espressoSystemInterface); @@ -119,8 +119,9 @@ void force_calc(CellStructure &cell_structure) { short_range_loop( add_bonded_force, - [](Particle &p1, Particle &p2, Distance const &d) { - add_non_bonded_pair_force(p1, p2, d.vec21, sqrt(d.dist2), d.dist2); + [&counter](Particle &p1, Particle &p2, Distance const &d) { + add_non_bonded_pair_force(p1, p2, d.vec21, sqrt(d.dist2), d.dist2, + counter); #ifdef COLLISION_DETECTION if (collision_params.mode != COLLISION_MODE_OFF) detect_collision(p1, p2, d.dist2); diff --git a/src/core/forces.hpp b/src/core/forces.hpp index f5ae3c8a380..3422a128105 100644 --- a/src/core/forces.hpp +++ b/src/core/forces.hpp @@ -38,15 +38,11 @@ extern ActorList forceActors; -/** \name Exported Functions */ -/************************************************************/ -/*@{*/ - /******************* forces.cpp *******************/ /** initialize real particle forces with thermostat forces and ghost particle forces with zero. */ -void init_forces(const ParticleRange &particles); +void init_forces(const ParticleRange &particles, uint64_t counter); /** Set forces of all ghosts to zero */ void init_forces_ghosts(const ParticleRange &particles); @@ -61,10 +57,9 @@ void init_forces_ghosts(const ParticleRange &particles); *
  • Calculate long range interaction forces * */ -void force_calc(CellStructure &cell_structure); +void force_calc(CellStructure &cell_structure, uint64_t counter); /** Calculate long range forces (P3M, ...). */ void calc_long_range_forces(const ParticleRange &particles); -/*@}*/ #endif diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp index 04571b60941..6f56411569d 100644 --- a/src/core/forces_inline.hpp +++ b/src/core/forces_inline.hpp @@ -100,25 +100,27 @@ inline ParticleForce external_force(Particle const &p) { return f; } -inline ParticleForce thermostat_force(Particle const &p) { +inline ParticleForce thermostat_force(Particle const &p, uint64_t counter) { extern LangevinThermostat langevin; if (!(thermo_switch & THERMO_LANGEVIN)) { return {}; } #ifdef ROTATION - return {friction_thermo_langevin(langevin, p), - p.p.rotation ? convert_vector_body_to_space( - p, friction_thermo_langevin_rotation(langevin, p)) - : Utils::Vector3d{}}; + return {friction_thermo_langevin(langevin, p, counter), + p.p.rotation + ? convert_vector_body_to_space( + p, friction_thermo_langevin_rotation(langevin, p, counter)) + : Utils::Vector3d{}}; #else - return friction_thermo_langevin(langevin, p); + return friction_thermo_langevin(langevin, p, counter); #endif } /** Initialize the forces for a real particle */ -inline ParticleForce init_local_particle_force(Particle const &part) { - return thermostat_force(part) + external_force(part); +inline ParticleForce init_local_particle_force(Particle const &part, + uint64_t counter) { + return thermostat_force(part, counter) + external_force(part); } inline Utils::Vector3d calc_non_bonded_pair_force_parts( @@ -234,7 +236,7 @@ inline Utils::Vector3d calc_non_bonded_pair_force(Particle const &p1, */ inline void add_non_bonded_pair_force(Particle &p1, Particle &p2, Utils::Vector3d const &d, double dist, - double dist2) { + double dist2, uint64_t counter) { IA_parameters const &ia_params = *get_ia_param(p1.p.type, p2.p.type); Utils::Vector3d force{}; Utils::Vector3d *torque1 = nullptr; @@ -289,7 +291,7 @@ inline void add_non_bonded_pair_force(Particle &p1, Particle &p2, /** The inter dpd force should not be part of the virial */ #ifdef DPD if (thermo_switch & THERMO_DPD) { - force += dpd_pair_force(p1, p2, ia_params, d, dist, dist2); + force += dpd_pair_force(p1, p2, ia_params, d, dist, dist2, counter); } #endif diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 397560c53ec..eebfebd3dfb 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -107,6 +107,7 @@ void integrator_sanity_checks() { * @return whether or not to stop the integration loop early. */ bool integrator_step_1(ParticleRange &particles) { + auto const counter = integrator_counter.value(); switch (integ_switch) { case INTEG_METHOD_STEEPEST_DESCENT: if (steepest_descent_step(particles)) @@ -117,7 +118,7 @@ bool integrator_step_1(ParticleRange &particles) { break; #ifdef NPT case INTEG_METHOD_NPT_ISO: - velocity_verlet_npt_step_1(particles); + velocity_verlet_npt_step_1(particles, counter); break; #endif case INTEG_METHOD_BD: @@ -126,7 +127,7 @@ bool integrator_step_1(ParticleRange &particles) { break; #ifdef STOKESIAN_DYNAMICS case INTEG_METHOD_SD: - stokesian_dynamics_step_1(particles); + stokesian_dynamics_step_1(particles, counter); break; #endif // STOKESIAN_DYNAMICS default: @@ -138,6 +139,7 @@ bool integrator_step_1(ParticleRange &particles) { /** Calls the hook of the propagation kernels after force calculation */ void integrator_step_2(ParticleRange &particles) { extern BrownianThermostat brownian; + auto const counter = integrator_counter.value(); switch (integ_switch) { case INTEG_METHOD_STEEPEST_DESCENT: // Nothing @@ -147,12 +149,12 @@ void integrator_step_2(ParticleRange &particles) { break; #ifdef NPT case INTEG_METHOD_NPT_ISO: - velocity_verlet_npt_step_2(particles); + velocity_verlet_npt_step_2(particles, counter); break; #endif case INTEG_METHOD_BD: // the Ermak-McCammon's Brownian Dynamics requires a single step - brownian_dynamics_propagator(brownian, particles); + brownian_dynamics_propagator(brownian, particles, counter); break; #ifdef STOKESIAN_DYNAMICS case INTEG_METHOD_SD: @@ -190,7 +192,7 @@ int integrate(int n_steps, int reuse_forces) { // Communication step: distribute ghost positions cells_update_ghosts(global_ghost_flags()); - force_calc(cell_structure); + force_calc(cell_structure, integrator_counter.value()); if (integ_switch != INTEG_METHOD_STEEPEST_DESCENT) { #ifdef ROTATION @@ -251,7 +253,7 @@ int integrate(int n_steps, int reuse_forces) { particles = cell_structure.local_particles(); - force_calc(cell_structure); + force_calc(cell_structure, integrator_counter.value()); #ifdef VIRTUAL_SITES virtual_sites()->after_force_calc(); diff --git a/src/core/integrators/brownian_inline.hpp b/src/core/integrators/brownian_inline.hpp index 84e9d00e059..2c8d33bd7bc 100644 --- a/src/core/integrators/brownian_inline.hpp +++ b/src/core/integrators/brownian_inline.hpp @@ -163,10 +163,12 @@ inline Utils::Vector3d bd_drag_vel(Thermostat::GammaType const &brownian_gamma, * From eq. (14.37) in @cite Schlick2010. * @param[in] brownian Parameters * @param[in] p %Particle + * @param[in] counter RNG counter * @param[in] dt Time interval */ inline Utils::Vector3d bd_random_walk(BrownianThermostat const &brownian, - Particle const &p, double dt) { + Particle const &p, uint64_t counter, + double dt) { // skip the translation thermalizing for virtual sites unless enabled extern bool thermo_virtual; if (p.p.is_virtual && !thermo_virtual) @@ -213,7 +215,7 @@ inline Utils::Vector3d bd_random_walk(BrownianThermostat const &brownian, // magnitude defined in the second eq. (14.38), Schlick2010. Utils::Vector3d delta_pos_body{}; auto const noise = Random::noise_gaussian( - integrator_counter.value(), brownian.rng_seed(), p.p.identity); + counter, brownian.rng_seed(), p.p.identity); for (int j = 0; j < 3; j++) { #ifdef EXTERNAL_FORCES if (!(p.p.ext_flag & COORD_FIXED(j))) @@ -265,9 +267,10 @@ inline Utils::Vector3d bd_random_walk(BrownianThermostat const &brownian, * From eq. (10.2.16) in @cite Pottier2010. * @param[in] brownian Parameters * @param[in] p %Particle + * @param[in] counter RNG counter */ inline Utils::Vector3d bd_random_walk_vel(BrownianThermostat const &brownian, - Particle const &p) { + Particle const &p, uint64_t counter) { // skip the translation thermalizing for virtual sites unless enabled extern bool thermo_virtual; if (p.p.is_virtual && !thermo_virtual) @@ -290,7 +293,7 @@ inline Utils::Vector3d bd_random_walk_vel(BrownianThermostat const &brownian, #endif // BROWNIAN_PER_PARTICLE auto const noise = Random::noise_gaussian( - integrator_counter.value(), brownian.rng_seed(), p.identity()); + counter, brownian.rng_seed(), p.identity()); Utils::Vector3d velocity = {}; for (int j = 0; j < 3; j++) { #ifdef EXTERNAL_FORCES @@ -395,10 +398,12 @@ bd_drag_vel_rot(Thermostat::GammaType const &brownian_gamma_rotation, * An analogy of eq. (14.37) in @cite Schlick2010. * @param[in] brownian Parameters * @param[in] p %Particle + * @param[in] counter RNG counter * @param[in] dt Time interval */ inline Utils::Vector4d bd_random_walk_rot(BrownianThermostat const &brownian, - Particle const &p, double dt) { + Particle const &p, uint64_t counter, + double dt) { // first, set defaults Thermostat::GammaType sigma_pos = brownian.sigma_pos_rotation; @@ -439,7 +444,7 @@ inline Utils::Vector4d bd_random_walk_rot(BrownianThermostat const &brownian, Utils::Vector3d dphi = {}; auto const noise = Random::noise_gaussian( - integrator_counter.value(), brownian.rng_seed(), p.p.identity); + counter, brownian.rng_seed(), p.p.identity); for (int j = 0; j < 3; j++) { #ifdef EXTERNAL_FORCES if (!(p.p.ext_flag & COORD_FIXED(j))) @@ -470,9 +475,11 @@ inline Utils::Vector4d bd_random_walk_rot(BrownianThermostat const &brownian, * An analogy of eq. (10.2.16) in @cite Pottier2010. * @param[in] brownian Parameters * @param[in] p %Particle + * @param[in] counter RNG counter */ inline Utils::Vector3d -bd_random_walk_vel_rot(BrownianThermostat const &brownian, Particle const &p) { +bd_random_walk_vel_rot(BrownianThermostat const &brownian, Particle const &p, + uint64_t counter) { double sigma_vel; // Override defaults if per-particle values for T and gamma are given @@ -490,7 +497,7 @@ bd_random_walk_vel_rot(BrownianThermostat const &brownian, Particle const &p) { Utils::Vector3d domega{}; auto const noise = Random::noise_gaussian( - integrator_counter.value(), brownian.rng_seed(), p.p.identity); + counter, brownian.rng_seed(), p.p.identity); for (int j = 0; j < 3; j++) { #ifdef EXTERNAL_FORCES if (!(p.p.ext_flag & COORD_FIXED(j))) @@ -504,15 +511,16 @@ bd_random_walk_vel_rot(BrownianThermostat const &brownian, Particle const &p) { #endif // ROTATION inline void brownian_dynamics_propagator(BrownianThermostat const &brownian, - const ParticleRange &particles) { + const ParticleRange &particles, + uint64_t counter) { for (auto &p : particles) { // Don't propagate translational degrees of freedom of vs extern bool thermo_virtual; if (!(p.p.is_virtual) or thermo_virtual) { p.r.p += bd_drag(brownian.gamma, p, time_step); p.m.v = bd_drag_vel(brownian.gamma, p); - p.r.p += bd_random_walk(brownian, p, time_step); - p.m.v += bd_random_walk_vel(brownian, p); + p.r.p += bd_random_walk(brownian, p, counter, time_step); + p.m.v += bd_random_walk_vel(brownian, p, counter); /* Verlet criterion check */ if ((p.r.p - p.l.p_old).norm2() > Utils::sqr(0.5 * skin)) cell_structure.set_resort_particles(Cells::RESORT_LOCAL); @@ -523,8 +531,8 @@ inline void brownian_dynamics_propagator(BrownianThermostat const &brownian, convert_torque_to_body_frame_apply_fix(p); p.r.quat = bd_drag_rot(brownian.gamma_rotation, p, time_step); p.m.omega = bd_drag_vel_rot(brownian.gamma_rotation, p); - p.r.quat = bd_random_walk_rot(brownian, p, time_step); - p.m.omega += bd_random_walk_vel_rot(brownian, p); + p.r.quat = bd_random_walk_rot(brownian, p, counter, time_step); + p.m.omega += bd_random_walk_vel_rot(brownian, p, counter); #endif // ROTATION } sim_time += time_step; diff --git a/src/core/integrators/langevin_inline.hpp b/src/core/integrators/langevin_inline.hpp index 65cd96e5d12..65e657e85ff 100644 --- a/src/core/integrators/langevin_inline.hpp +++ b/src/core/integrators/langevin_inline.hpp @@ -36,10 +36,11 @@ * LANGEVIN_PER_PARTICLE). Applies the noise and friction term. * @param[in] langevin Parameters * @param[in] p %Particle + * @param[in] counter RNG counter */ inline Utils::Vector3d -friction_thermo_langevin(LangevinThermostat const &langevin, - Particle const &p) { +friction_thermo_langevin(LangevinThermostat const &langevin, Particle const &p, + uint64_t counter) { // Early exit for virtual particles without thermostat if (p.p.is_virtual && !thermo_virtual) { return {}; @@ -85,9 +86,8 @@ friction_thermo_langevin(LangevinThermostat const &langevin, #endif // PARTICLE_ANISOTROPY return friction_op * velocity + - noise_op * - Random::noise_uniform( - integrator_counter.value(), langevin.rng_seed(), p.p.identity); + noise_op * Random::noise_uniform( + counter, langevin.rng_seed(), p.p.identity); } #ifdef ROTATION @@ -97,10 +97,11 @@ friction_thermo_langevin(LangevinThermostat const &langevin, * LANGEVIN_PER_PARTICLE). Applies the noise and friction term. * @param[in] langevin Parameters * @param[in] p %Particle + * @param[in] counter RNG counter */ inline Utils::Vector3d friction_thermo_langevin_rotation(LangevinThermostat const &langevin, - Particle const &p) { + Particle const &p, uint64_t counter) { auto pref_friction = -langevin.gamma_rotation; auto pref_noise = langevin.pref_noise_rotation; @@ -118,7 +119,7 @@ friction_thermo_langevin_rotation(LangevinThermostat const &langevin, #endif // LANGEVIN_PER_PARTICLE auto const noise = Random::noise_uniform( - integrator_counter.value(), langevin.rng_seed(), p.p.identity); + counter, langevin.rng_seed(), p.p.identity); return hadamard_product(pref_friction, p.m.omega) + hadamard_product(pref_noise, noise); } diff --git a/src/core/integrators/stokesian_dynamics_inline.hpp b/src/core/integrators/stokesian_dynamics_inline.hpp index 32ff0d24827..556a0100c18 100644 --- a/src/core/integrators/stokesian_dynamics_inline.hpp +++ b/src/core/integrators/stokesian_dynamics_inline.hpp @@ -31,12 +31,12 @@ #include -inline void -stokesian_dynamics_propagate_vel_pos(const ParticleRange &particles) { +inline void stokesian_dynamics_propagate_vel_pos(const ParticleRange &particles, + uint64_t counter) { auto const skin2 = Utils::sqr(0.5 * skin); // Compute new (translational and rotational) velocities - propagate_vel_pos_sd(particles, comm_cart, time_step); + propagate_vel_pos_sd(particles, comm_cart, time_step, counter); for (auto &p : particles) { @@ -56,8 +56,9 @@ stokesian_dynamics_propagate_vel_pos(const ParticleRange &particles) { } } -inline void stokesian_dynamics_step_1(const ParticleRange &particles) { - stokesian_dynamics_propagate_vel_pos(particles); +inline void stokesian_dynamics_step_1(const ParticleRange &particles, + uint64_t counter) { + stokesian_dynamics_propagate_vel_pos(particles, counter); sim_time += time_step; } diff --git a/src/core/integrators/velocity_verlet_npt.cpp b/src/core/integrators/velocity_verlet_npt.cpp index 797d208c215..7a5b04a4ac8 100644 --- a/src/core/integrators/velocity_verlet_npt.cpp +++ b/src/core/integrators/velocity_verlet_npt.cpp @@ -33,7 +33,10 @@ #include -void velocity_verlet_npt_propagate_vel_final(const ParticleRange &particles) { +#include + +void velocity_verlet_npt_propagate_vel_final(const ParticleRange &particles, + uint64_t counter) { extern IsotropicNptThermostat npt_iso; nptiso.p_vel[0] = nptiso.p_vel[1] = nptiso.p_vel[2] = 0.0; @@ -41,7 +44,8 @@ void velocity_verlet_npt_propagate_vel_final(const ParticleRange &particles) { // Virtual sites are not propagated during integration if (p.p.is_virtual) continue; - auto const noise = friction_therm0_nptiso<2>(npt_iso, p.m.v, p.p.identity); + auto const noise = + friction_therm0_nptiso<2>(npt_iso, p.m.v, p.p.identity, counter); for (int j = 0; j < 3; j++) { if (!(p.p.ext_flag & COORD_FIXED(j))) { if (nptiso.geometry & nptiso.nptgeom_dir[j]) { @@ -58,7 +62,7 @@ void velocity_verlet_npt_propagate_vel_final(const ParticleRange &particles) { } /** Scale and communicate instantaneous NpT pressure */ -void velocity_verlet_npt_finalize_p_inst() { +void velocity_verlet_npt_finalize_p_inst(uint64_t counter) { extern IsotropicNptThermostat npt_iso; /* finalize derivation of p_inst */ nptiso.p_inst = 0.0; @@ -75,15 +79,16 @@ void velocity_verlet_npt_finalize_p_inst() { nptiso.p_inst = p_sum / (nptiso.dimension * nptiso.volume); nptiso.p_diff = nptiso.p_diff + (nptiso.p_inst - nptiso.p_ext) * 0.5 * time_step + - friction_thermV_nptiso(npt_iso, nptiso.p_diff); + friction_thermV_nptiso(npt_iso, nptiso.p_diff, counter); } } -void velocity_verlet_npt_propagate_pos(const ParticleRange &particles) { +void velocity_verlet_npt_propagate_pos(const ParticleRange &particles, + uint64_t counter) { double scal[3] = {0., 0., 0.}, L_new = 0.0; /* finalize derivation of p_inst */ - velocity_verlet_npt_finalize_p_inst(); + velocity_verlet_npt_finalize_p_inst(counter); /* adjust \ref nptiso_struct::nptiso.volume; prepare pos- and * vel-rescaling @@ -152,7 +157,8 @@ void velocity_verlet_npt_propagate_pos(const ParticleRange &particles) { cells_re_init(cell_structure.decomposition_type()); } -void velocity_verlet_npt_propagate_vel(const ParticleRange &particles) { +void velocity_verlet_npt_propagate_vel(const ParticleRange &particles, + uint64_t counter) { extern IsotropicNptThermostat npt_iso; #ifdef NPT nptiso.p_vel[0] = nptiso.p_vel[1] = nptiso.p_vel[2] = 0.0; @@ -170,7 +176,7 @@ void velocity_verlet_npt_propagate_vel(const ParticleRange &particles) { if (!(p.p.ext_flag & COORD_FIXED(j))) { #ifdef NPT auto const noise = - friction_therm0_nptiso<1>(npt_iso, p.m.v, p.p.identity); + friction_therm0_nptiso<1>(npt_iso, p.m.v, p.p.identity, counter); if (integ_switch == INTEG_METHOD_NPT_ISO && (nptiso.geometry & nptiso.nptgeom_dir[j])) { p.m.v[j] += (p.f.f[j] * time_step / 2.0 + noise[j]) / p.p.mass; @@ -184,17 +190,19 @@ void velocity_verlet_npt_propagate_vel(const ParticleRange &particles) { } } -void velocity_verlet_npt_step_1(const ParticleRange &particles) { - velocity_verlet_npt_propagate_vel(particles); - velocity_verlet_npt_propagate_pos(particles); +void velocity_verlet_npt_step_1(const ParticleRange &particles, + uint64_t counter) { + velocity_verlet_npt_propagate_vel(particles, counter); + velocity_verlet_npt_propagate_pos(particles, counter); sim_time += time_step; } -void velocity_verlet_npt_step_2(const ParticleRange &particles) { - velocity_verlet_npt_propagate_vel_final(particles); +void velocity_verlet_npt_step_2(const ParticleRange &particles, + uint64_t counter) { + velocity_verlet_npt_propagate_vel_final(particles, counter); #ifdef ROTATION convert_torques_propagate_omega(particles); #endif - velocity_verlet_npt_finalize_p_inst(); + velocity_verlet_npt_finalize_p_inst(counter); } #endif diff --git a/src/core/integrators/velocity_verlet_npt.hpp b/src/core/integrators/velocity_verlet_npt.hpp index b1a9a8b3b18..5d07b398f1a 100644 --- a/src/core/integrators/velocity_verlet_npt.hpp +++ b/src/core/integrators/velocity_verlet_npt.hpp @@ -25,6 +25,8 @@ #include "Particle.hpp" #include "ParticleRange.hpp" +#include + /** Special propagator for NpT isotropic. * Propagate the velocities and positions. Integration steps before force * calculation of the Velocity Verlet integrator: @@ -34,14 +36,16 @@ * Propagate pressure, box_length (2 times) and positions, rescale * positions and velocities and check Verlet list criterion (only NpT). */ -void velocity_verlet_npt_step_1(const ParticleRange &particles); +void velocity_verlet_npt_step_1(const ParticleRange &particles, + uint64_t counter); /** Final integration step of the Velocity Verlet+NpT integrator. * Finalize instantaneous pressure calculation: * \f[ v(t+\Delta t) = v(t+0.5 \Delta t) * + 0.5 \Delta t \cdot F(t+\Delta t)/m \f] */ -void velocity_verlet_npt_step_2(const ParticleRange &particles); +void velocity_verlet_npt_step_2(const ParticleRange &particles, + uint64_t counter); #endif #endif diff --git a/src/core/stokesian_dynamics/sd_interface.cpp b/src/core/stokesian_dynamics/sd_interface.cpp index 7a01dc0cf88..0ab5bd1d0a2 100644 --- a/src/core/stokesian_dynamics/sd_interface.cpp +++ b/src/core/stokesian_dynamics/sd_interface.cpp @@ -180,7 +180,7 @@ int get_sd_flags() { return sd_flags; } void propagate_vel_pos_sd(const ParticleRange &particles, const boost::mpi::communicator &comm, - const double time_step) { + const double time_step, uint64_t counter) { static std::vector parts_buffer{}; parts_buffer.clear(); @@ -231,7 +231,7 @@ void propagate_vel_pos_sd(const ParticleRange &particles, case CPU: v_sd = sd_cpu(x_host, f_host, a_host, n_part, sd_viscosity, std::sqrt(sd_kT / time_step), - static_cast(integrator_counter.value()), + static_cast(counter), static_cast(stokesian.rng_seed()), sd_flags); break; #endif @@ -240,7 +240,7 @@ void propagate_vel_pos_sd(const ParticleRange &particles, case GPU: v_sd = sd_gpu(x_host, f_host, a_host, n_part, sd_viscosity, std::sqrt(sd_kT / time_step), - static_cast(integrator_counter.value()), + static_cast(counter), static_cast(stokesian.rng_seed()), sd_flags); break; #endif diff --git a/src/core/stokesian_dynamics/sd_interface.hpp b/src/core/stokesian_dynamics/sd_interface.hpp index 241ed5265b9..ec5e8b4e873 100644 --- a/src/core/stokesian_dynamics/sd_interface.hpp +++ b/src/core/stokesian_dynamics/sd_interface.hpp @@ -57,7 +57,7 @@ int get_sd_flags(); */ void propagate_vel_pos_sd(const ParticleRange &particles, const boost::mpi::communicator &comm, - double time_step); + double time_step, uint64_t counter); #endif // STOKESIAN_DYNAMICS diff --git a/src/core/thermostat.hpp b/src/core/thermostat.hpp index be434688438..83fbdf458f2 100644 --- a/src/core/thermostat.hpp +++ b/src/core/thermostat.hpp @@ -358,13 +358,15 @@ void thermo_init(); * @param npt_iso Parameters * @param vel particle velocity * @param p_identity particle identity + * @param counter RNG counter * @return noise added to the velocity, already rescaled by * dt/2 (contained in prefactors) */ template inline Utils::Vector3d friction_therm0_nptiso(IsotropicNptThermostat const &npt_iso, - Utils::Vector3d const &vel, int p_identity) { + Utils::Vector3d const &vel, int p_identity, + uint64_t counter) { static_assert(step == 1 or step == 2, "NPT only has 2 integration steps"); constexpr auto const salt = (step == 1) ? RNGSalt::NPTISO0_HALF_STEP1 : RNGSalt::NPTISO0_HALF_STEP2; @@ -372,8 +374,8 @@ friction_therm0_nptiso(IsotropicNptThermostat const &npt_iso, if (npt_iso.pref_noise_0 > 0.0) { return npt_iso.pref_rescale_0 * vel + npt_iso.pref_noise_0 * - Random::noise_uniform(integrator_counter.value(), - npt_iso.rng_seed(), p_identity); + Random::noise_uniform(counter, npt_iso.rng_seed(), + p_identity); } return npt_iso.pref_rescale_0 * vel; } @@ -384,13 +386,12 @@ friction_therm0_nptiso(IsotropicNptThermostat const &npt_iso, * nptiso_struct::p_diff */ inline double friction_thermV_nptiso(IsotropicNptThermostat const &npt_iso, - double p_diff) { + double p_diff, uint64_t counter) { if (thermo_switch & THERMO_NPT_ISO) { if (npt_iso.pref_noise_V > 0.0) { return npt_iso.pref_rescale_V * p_diff + - npt_iso.pref_noise_V * - Random::noise_uniform( - integrator_counter.value(), npt_iso.rng_seed(), 0); + npt_iso.pref_noise_V * Random::noise_uniform( + counter, npt_iso.rng_seed(), 0); } return npt_iso.pref_rescale_V * p_diff; } diff --git a/src/core/unit_tests/thermostats_test.cpp b/src/core/unit_tests/thermostats_test.cpp index b14b82c4878..08a708adfbc 100644 --- a/src/core/unit_tests/thermostats_test.cpp +++ b/src/core/unit_tests/thermostats_test.cpp @@ -23,12 +23,12 @@ #define BOOST_TEST_DYN_LINK #include #include +#include #include #include #include "Particle.hpp" -#include "integrate.hpp" #include "integrators/brownian_inline.hpp" #include "integrators/langevin_inline.hpp" #include "random_test.hpp" @@ -38,9 +38,6 @@ extern double time_step; extern double temperature; -void reset_integrator_counter() { - integrator_counter = Utils::Counter{0ul}; -} // multiply by 100 because BOOST_CHECK_CLOSE takes a percentage tolerance, // and by 6 to account for error accumulation in thermostat functions @@ -72,7 +69,6 @@ template T thermostat_factory() { BOOST_AUTO_TEST_CASE(test_brownian_dynamics) { time_step = 0.1; temperature = 3.0; - reset_integrator_counter(); auto const brownian = thermostat_factory(); auto const dispersion = hadamard_division(particle_factory().f.f, brownian.gamma); @@ -103,7 +99,7 @@ BOOST_AUTO_TEST_CASE(test_brownian_dynamics) { auto const sigma = sqrt(brownian.gamma / (2.0 * temperature)); auto const noise = Random::noise_gaussian(0, 0, 0); auto const ref = hadamard_division(noise, sigma) * sqrt(time_step); - auto const out = bd_random_walk(brownian, p, time_step); + auto const out = bd_random_walk(brownian, p, 0u, time_step); BOOST_CHECK_CLOSE(out[0], ref[0], tol); BOOST_CHECK_CLOSE(out[1], ref[1], tol); BOOST_CHECK_CLOSE(out[2], ref[2], tol); @@ -115,7 +111,7 @@ BOOST_AUTO_TEST_CASE(test_brownian_dynamics) { auto const sigma = sqrt(temperature); auto const noise = Random::noise_gaussian(0, 0, 0); auto const ref = sigma * noise / sqrt(p.p.mass); - auto const out = bd_random_walk_vel(brownian, p); + auto const out = bd_random_walk_vel(brownian, p, 0u); BOOST_CHECK_CLOSE(out[0], ref[0], tol); BOOST_CHECK_CLOSE(out[1], ref[1], tol); BOOST_CHECK_CLOSE(out[2], ref[2], tol); @@ -156,7 +152,7 @@ BOOST_AUTO_TEST_CASE(test_brownian_dynamics) { auto const noise = Random::noise_gaussian(0, 0, 0); auto const phi = hadamard_division(noise, sigma)[0] * sqrt(time_step); - auto const out = bd_random_walk_rot(brownian, p, time_step); + auto const out = bd_random_walk_rot(brownian, p, 0u, time_step); BOOST_CHECK_CLOSE(out[0], std::cos(phi / 2), tol); BOOST_CHECK_CLOSE(out[1], std::sin(phi / 2), tol); BOOST_CHECK_CLOSE(out[2], 0.0, tol); @@ -171,7 +167,7 @@ BOOST_AUTO_TEST_CASE(test_brownian_dynamics) { auto const noise = Random::noise_gaussian(0, 0, 0); auto const ref = sigma * noise; - auto const out = bd_random_walk_vel_rot(brownian, p); + auto const out = bd_random_walk_vel_rot(brownian, p, 0u); BOOST_CHECK_CLOSE(out[0], ref[0], tol); BOOST_CHECK_CLOSE(out[1], ref[1], tol); BOOST_CHECK_CLOSE(out[2], ref[2], tol); @@ -182,7 +178,6 @@ BOOST_AUTO_TEST_CASE(test_brownian_dynamics) { BOOST_AUTO_TEST_CASE(test_langevin_dynamics) { time_step = 0.1; temperature = 3.0; - reset_integrator_counter(); auto const langevin = thermostat_factory(); auto const prefactor_squared = 24.0 * temperature / time_step; @@ -194,7 +189,7 @@ BOOST_AUTO_TEST_CASE(test_langevin_dynamics) { auto const pref = sqrt(prefactor_squared * langevin.gamma); auto const ref = hadamard_product(-langevin.gamma, p.m.v) + hadamard_product(pref, noise); - auto const out = friction_thermo_langevin(langevin, p); + auto const out = friction_thermo_langevin(langevin, p, 0u); BOOST_CHECK_CLOSE(out[0], ref[0], tol); BOOST_CHECK_CLOSE(out[1], ref[1], tol); BOOST_CHECK_CLOSE(out[2], ref[2], tol); @@ -209,7 +204,7 @@ BOOST_AUTO_TEST_CASE(test_langevin_dynamics) { auto const pref = sqrt(prefactor_squared * langevin.gamma_rotation); auto const ref = hadamard_product(-langevin.gamma_rotation, p.m.omega) + hadamard_product(pref, noise); - auto const out = friction_thermo_langevin_rotation(langevin, p); + auto const out = friction_thermo_langevin_rotation(langevin, p, 0u); BOOST_CHECK_CLOSE(out[0], ref[0], tol); BOOST_CHECK_CLOSE(out[1], ref[1], tol); BOOST_CHECK_CLOSE(out[2], ref[2], tol); @@ -221,7 +216,7 @@ BOOST_AUTO_TEST_CASE(test_noise_statistics) { time_step = 1.0; temperature = 2.0; constexpr size_t const sample_size = 10'000; - reset_integrator_counter(); + uint64_t counter = 0u; auto thermostat = thermostat_factory(); auto p1 = particle_factory(); auto p2 = particle_factory(); @@ -229,11 +224,11 @@ BOOST_AUTO_TEST_CASE(test_noise_statistics) { p2.p.identity = 1; auto const correlation = std::get<3>(noise_statistics( - [&p1, &p2, &thermostat]() -> std::array { - integrator_counter.increment(); - return {friction_thermo_langevin(thermostat, p1), - -friction_thermo_langevin(thermostat, p1), - friction_thermo_langevin(thermostat, p2)}; + [&p1, &p2, &thermostat, &counter]() -> std::array { + counter++; + return {friction_thermo_langevin(thermostat, p1, counter), + -friction_thermo_langevin(thermostat, p1, counter), + friction_thermo_langevin(thermostat, p2, counter)}; }, sample_size)); for (size_t i = 0; i < correlation.size(); ++i) { @@ -255,7 +250,7 @@ BOOST_AUTO_TEST_CASE(test_brownian_randomness) { time_step = 1.0; temperature = 2.0; constexpr size_t const sample_size = 10'000; - reset_integrator_counter(); + uint64_t counter = 0u; auto thermostat = thermostat_factory(); auto p = particle_factory(); #ifdef ROTATION @@ -266,14 +261,14 @@ BOOST_AUTO_TEST_CASE(test_brownian_randomness) { #endif auto const correlation = std::get<3>(noise_statistics( - [&p, &thermostat]() -> std::array { - integrator_counter.increment(); + [&p, &thermostat, &counter]() -> std::array { + counter++; return { - bd_random_walk(thermostat, p, time_step), - bd_random_walk_vel(thermostat, p), + bd_random_walk(thermostat, p, counter, time_step), + bd_random_walk_vel(thermostat, p, counter), #ifdef ROTATION - bd_random_walk_rot(thermostat, p, time_step), - bd_random_walk_vel_rot(thermostat, p), + bd_random_walk_rot(thermostat, p, counter, time_step), + bd_random_walk_vel_rot(thermostat, p, counter), #endif }; }, @@ -289,7 +284,7 @@ BOOST_AUTO_TEST_CASE(test_langevin_randomness) { time_step = 1.0; temperature = 2.0; constexpr size_t const sample_size = 10'000; - reset_integrator_counter(); + uint64_t counter = 0u; auto thermostat = thermostat_factory(); auto p = particle_factory(); #ifdef ROTATION @@ -299,12 +294,12 @@ BOOST_AUTO_TEST_CASE(test_langevin_randomness) { #endif auto const correlation = std::get<3>(noise_statistics( - [&p, &thermostat]() -> std::array { - integrator_counter.increment(); + [&p, &thermostat, &counter]() -> std::array { + counter++; return { - friction_thermo_langevin(thermostat, p), + friction_thermo_langevin(thermostat, p, counter), #ifdef ROTATION - friction_thermo_langevin_rotation(thermostat, p), + friction_thermo_langevin_rotation(thermostat, p, counter), #endif }; }, @@ -323,7 +318,7 @@ BOOST_AUTO_TEST_CASE(test_npt_iso_randomness) { time_step = 1.0; temperature = 2.0; constexpr size_t const sample_size = 10'000; - reset_integrator_counter(); + uint64_t counter = 0u; IsotropicNptThermostat thermostat{}; thermostat.rng_initialize(0); thermostat.gamma0 = 2.0; @@ -332,12 +327,12 @@ BOOST_AUTO_TEST_CASE(test_npt_iso_randomness) { auto p = particle_factory(); auto const correlation = std::get<3>(noise_statistics( - [&p, &thermostat]() -> std::array { - integrator_counter.increment(); + [&p, &thermostat, &counter]() -> std::array { + counter++; return { - friction_therm0_nptiso<1>(thermostat, p.m.v, 0), - friction_therm0_nptiso<2>(thermostat, p.m.v, 0), - friction_thermV_nptiso(thermostat, 1.5), + friction_therm0_nptiso<1>(thermostat, p.m.v, 0, counter), + friction_therm0_nptiso<2>(thermostat, p.m.v, 0, counter), + friction_thermV_nptiso(thermostat, 1.5, counter), }; }, sample_size));