Refactor thermostat infrastructure (#3888)

Description of changes:
- Remove RNG correlation stemming from seed offsets (fixes #3585)
    - seeds are now used as keys
    - a monotonically increasing counter is used in each thermostat
    - the only way to reset these counters is to create a new `System`
- Remove RNG correlation stemming from resetting `sim_time` or `time_step` during simulations with SD (fixes #3840)
    - the SD thermostat now uses the same RNG interface as other thermostats
- Accelerate RNG unit tests (fixes #3573)
    - they now take 2 seconds to run in coverage and sanitizer builds in CI
- Separate thermostats from integrators
    - better separation of concerns

doc/sphinx/system_setup.rst

@@ -212,6 +212,13 @@ the temperature of a thermostat, you actually do not define the
 temperature, but the value of the thermal energy :math:`k_B T` in the
 current unit system (see the discussion on units, Section :ref:`On units`).
 
+All thermostats have a ``seed`` argument that controls the state of the random
+number generator (Philox Counter-based RNG). This seed is required on first
+activation of a thermostat, unless stated otherwise. It can be omitted in
+subsequent calls of the method that activates the same thermostat. The random
+sequence also depends on the thermostats counters that are
+incremented after each integration step.
+
 .. _Langevin thermostat:
 
 Langevin thermostat
@@ -253,12 +260,6 @@ The distribution of :math:`\overline{\eta}` is uniform and satisfies
 
 .. math:: <\overline{\eta}> = 0 , <\overline{\eta}\overline{\eta}> = 1/dt
 
-The keyword ``seed`` controls the state of the random number generator (Philox
-Counter-based RNG) and is required on first activation of the thermostat. It
-can be omitted in subsequent calls of ``set_langevin()``. It is the user's
-responsibility to decide whether the thermostat should be deterministic (by
-using a fixed seed) or not (by using a randomized seed).
-
 If the feature ``ROTATION`` is compiled in, the rotational degrees of freedom are
 also coupled to the thermostat. If only the first two arguments are
 specified then the friction coefficient for the rotation is set to the
@@ -301,10 +302,8 @@ The LB fluid can be used to thermalize particles, while also including their hyd
 The LB thermostat expects an instance of either :class:`espressomd.lb.LBFluid` or :class:`espressomd.lb.LBFluidGPU`.
 Temperature is set via the ``kT`` argument of the LB fluid.
 
-Furthermore a ``seed`` has to be given for the
-thermalization of the particle coupling. The magnitude of the frictional coupling can be adjusted by
-the parameter ``gamma``.
-To enable the LB thermostat, use::
+The magnitude of the frictional coupling can be adjusted by the
+parameter ``gamma``. To enable the LB thermostat, use::
 
     import espressomd
     import espressomd.lb
@@ -471,12 +470,6 @@ Note: the rotational Brownian dynamics implementation is compatible with particl
 the isotropic moment of inertia tensor only. Otherwise, the viscous terminal angular velocity
 is not defined, i.e. it has no constant direction over the time.
 
-The keyword ``seed`` controls the state of the random number generator (Philox
-Counter-based RNG) and is required on first activation of the thermostat. It
-can be omitted in subsequent calls of ``set_brownian()``. It is the user's
-responsibility to decide whether the thermostat should be deterministic (by
-using a fixed seed) or not (by using a randomized seed).
-
 .. _Stokesian thermostat:
 
 Stokesian thermostat
@@ -501,8 +494,7 @@ needs to be activated via::
     system.integrator.run(100)
 
 where ``kT`` denotes the desired temperature of the system, and ``seed`` the
-seed for the random number generator of the Stokesian Dynamics thermostat.
-It is independent from the other random number generators in |es|.
+seed for the random number generator.
 
 
 .. _CUDA:

src/core/bonded_interactions/thermalized_bond.cpp

@@ -27,6 +27,7 @@
 #include "bonded_interaction_data.hpp"
 #include "communication.hpp"
 #include "global.hpp"
+#include "integrate.hpp"
 
 int n_thermalized_bonds = 0;
 

src/core/bonded_interactions/thermalized_bond.hpp

@@ -33,7 +33,6 @@ extern int n_thermalized_bonds;
 
 #include "Particle.hpp"
 #include "bonded_interaction_data.hpp"
-#include "integrate.hpp"
 #include "random.hpp"
 #include "thermostat.hpp"
 
@@ -77,7 +76,8 @@ thermalized_bond_forces(Particle const &p1, Particle const &p2,
   Utils::Vector3d force1{};
   Utils::Vector3d force2{};
   auto const noise = Random::noise_uniform<RNGSalt::THERMALIZED_BOND>(
-      thermalized_bond.rng_get(), p1.p.identity, p2.p.identity);
+      thermalized_bond.rng_counter(), thermalized_bond.rng_seed(),
+      p1.p.identity, p2.p.identity);
 
   for (int i = 0; i < 3; i++) {
     double force_lv_com, force_lv_dist;

src/core/communication.cpp

@@ -97,7 +97,6 @@ int n_nodes = -1;
   CB(mpi_bcast_coulomb_params_slave)                                           \
   CB(mpi_remove_particle_slave)                                                \
   CB(mpi_rescale_particles_slave)                                              \
-  CB(mpi_bcast_nptiso_geom_slave)                                              \
   CB(mpi_bcast_steepest_descent_worker)                                        \
   CB(mpi_bcast_cuda_global_part_vars_slave)                                    \
   CB(mpi_resort_particles_slave)                                               \
@@ -496,17 +495,14 @@ void mpi_set_use_verlet_lists(bool use_verlet_lists) {
 
 /*************** BCAST NPTISO GEOM *****************/
 
-void mpi_bcast_nptiso_geom() {
-  mpi_call(mpi_bcast_nptiso_geom_slave, -1, 0);
-  mpi_bcast_nptiso_geom_slave(-1, 0);
-}
+#ifdef NPT
+REGISTER_CALLBACK(mpi_bcast_nptiso_geom_worker)
 
-void mpi_bcast_nptiso_geom_slave(int, int) {
-  MPI_Bcast(&nptiso.geometry, 1, MPI_INT, 0, comm_cart);
-  MPI_Bcast(&nptiso.dimension, 1, MPI_INT, 0, comm_cart);
-  MPI_Bcast(&nptiso.cubic_box, 1, MPI_LOGICAL, 0, comm_cart);
-  MPI_Bcast(&nptiso.non_const_dim, 1, MPI_INT, 0, comm_cart);
+void mpi_bcast_nptiso_geom() {
+  mpi_call(mpi_bcast_nptiso_geom_worker, -1, 0);
+  mpi_bcast_nptiso_geom_worker(-1, 0);
 }
+#endif
 
 /*************** BCAST STEEPEST DESCENT *****************/
 

src/core/constraints/ShapeBasedConstraint.cpp

@@ -24,6 +24,7 @@
 #include "errorhandling.hpp"
 #include "forces_inline.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
+#include "thermostat.hpp"
 
 namespace Constraints {
 Utils::Vector3d ShapeBasedConstraint::total_force() const {
@@ -77,7 +78,7 @@ ParticleForce ShapeBasedConstraint::force(Particle const &p,
         force1 +=
             dpd_pair_force(p, part_rep, ia_params, dist_vec, dist, dist * dist);
         // Additional use of DPD here requires counter increase
-        dpd_rng_counter_increment();
+        dpd.rng_increment();
       }
 #endif
     } else if (m_penetrable && (dist <= 0)) {
@@ -89,7 +90,7 @@ ParticleForce ShapeBasedConstraint::force(Particle const &p,
           force1 += dpd_pair_force(p, part_rep, ia_params, dist_vec, dist,
                                    dist * dist);
           // Additional use of DPD here requires counter increase
-          dpd_rng_counter_increment();
+          dpd.rng_increment();
         }
 #endif
       }

src/core/dpd.cpp

@@ -47,9 +47,9 @@ using Utils::Vector3d;
  *     seed-per-node)
  */
 Vector3d dpd_noise(uint32_t pid1, uint32_t pid2) {
-  extern DPDThermostat dpd;
   return Random::noise_uniform<RNGSalt::SALT_DPD>(
-      dpd.rng_get(), (pid1 < pid2) ? pid2 : pid1, (pid1 < pid2) ? pid1 : pid2);
+      dpd.rng_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,

src/core/forces.cpp

@@ -49,7 +49,7 @@
 
 ActorList forceActors;
 
-void init_forces(const ParticleRange &particles) {
+void init_forces(const ParticleRange &particles, double time_step) {
   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, time_step);
   }
 
   /* 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, double time_step) {
   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, time_step);
 
   for (auto &forceActor : forceActors) {
     forceActor->computeForces(espressoSystemInterface);

src/core/forces.hpp

@@ -38,15 +38,9 @@
 
 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, double time_step);
 
 /** Set forces of all ghosts to zero */
 void init_forces_ghosts(const ParticleRange &particles);
@@ -61,10 +55,9 @@ void init_forces_ghosts(const ParticleRange &particles);
  *  <li> Calculate long range interaction forces
  *  </ol>
  */
-void force_calc(CellStructure &cell_structure);
+void force_calc(CellStructure &cell_structure, double time_step);
 
 /** Calculate long range forces (P3M, ...). */
 void calc_long_range_forces(const ParticleRange &particles);
-/*@}*/
 
 #endif

src/core/forces_inline.hpp

@@ -39,7 +39,6 @@
 #include "forces.hpp"
 #include "immersed_boundary/ibm_tribend.hpp"
 #include "immersed_boundary/ibm_triel.hpp"
-#include "integrators/langevin_inline.hpp"
 #include "nonbonded_interactions/bmhtf-nacl.hpp"
 #include "nonbonded_interactions/buckingham.hpp"
 #include "nonbonded_interactions/gaussian.hpp"
@@ -61,6 +60,7 @@
 #include "object-in-fluid/oif_local_forces.hpp"
 #include "rotation.hpp"
 #include "thermostat.hpp"
+#include "thermostats/langevin_inline.hpp"
 
 #ifdef DIPOLES
 #include "electrostatics_magnetostatics/dipole_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, double time_step) {
   extern LangevinThermostat langevin;
   if (!(thermo_switch & THERMO_LANGEVIN)) {
     return {};
   }
 
 #ifdef ROTATION
-  return {friction_thermo_langevin(langevin, p),
+  return {friction_thermo_langevin(langevin, p, time_step),
           p.p.rotation ? convert_vector_body_to_space(
-                             p, friction_thermo_langevin_rotation(langevin, p))
+                             p, friction_thermo_langevin_rotation(langevin, p,
+                                                                  time_step))
                        : Utils::Vector3d{}};
 #else
-  return friction_thermo_langevin(langevin, p);
+  return friction_thermo_langevin(langevin, p, time_step);
 #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,
+                                               double time_step) {
+  return thermostat_force(part, time_step) + external_force(part);
 }
 
 inline Utils::Vector3d calc_non_bonded_pair_force_parts(

src/core/global.cpp

@@ -29,6 +29,7 @@
 #include "event.hpp"
 #include "grid.hpp"
 #include "grid_based_algorithms/lb_interface.hpp"
+#include "integrate.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
 #include "npt.hpp"
 #include "object-in-fluid/oif_global_forces.hpp"
@@ -45,9 +46,6 @@
 #include <unordered_map>
 
 extern double force_cap;
-extern LangevinThermostat langevin;
-extern BrownianThermostat brownian;
-extern IsotropicNptThermostat npt_iso;
 
 namespace {
 

src/core/grid_based_algorithms/lb.cpp

@@ -42,10 +42,10 @@
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
 #include "random.hpp"
 
-#include "utils/u32_to_u64.hpp"
 #include <utils/Counter.hpp>
 #include <utils/index.hpp>
 #include <utils/math/matrix_vector_product.hpp>
+#include <utils/u32_to_u64.hpp>
 #include <utils/uniform.hpp>
 using Utils::get_linear_index;
 #include <utils/constants.hpp>
@@ -841,7 +841,7 @@ lb_thermalize_modes(Lattice::index_t index, const std::array<T, 19> &modes,
     using rng_type = r123::Philox4x64;
     using ctr_type = rng_type::ctr_type;
 
-    const r123::Philox4x64::ctr_type c{
+    const ctr_type c{
         {rng_counter->value(), static_cast<uint64_t>(RNGSalt::FLUID)}};
     const T rootdensity =
         std::sqrt(std::fabs(modes[0] + lb_parameters.density));

src/core/grid_based_algorithms/lb_particle_coupling.cpp

@@ -279,7 +279,7 @@ void lb_lbcoupling_calc_particle_lattice_ia(
         auto f_random = [noise_amplitude](int id) -> Utils::Vector3d {
           if (noise_amplitude > 0.0) {
             return Random::noise_uniform<RNGSalt::PARTICLES>(
-                lb_particle_coupling.rng_counter_coupling->value(), id);
+                lb_particle_coupling.rng_counter_coupling->value(), 0, id);
           }
           return {};
         };