Use global counter in thermostats

doc/sphinx/system_setup.rst

@@ -210,6 +210,15 @@ 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 a monotonically increasing global counter that is
+incremented after each integration step. It is the user's responsibility to
+decide whether the thermostats should be deterministic (by using a fixed seed)
+or not (by using a randomized seed).
+
 .. _Langevin thermostat:
 
 Langevin thermostat
@@ -251,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
@@ -299,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::
 
     system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1.5)
 
@@ -465,12 +466,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
@@ -487,8 +482,7 @@ needs to be activated via::
     system.thermostat.set_stokesian(kT=1.0, seed=43)
 
 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.hpp

@@ -77,7 +77,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);
+      integrator_counter.value(), 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/constraints/ShapeBasedConstraint.cpp

@@ -74,10 +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);
-        // Additional use of DPD here requires counter increase
-        dpd_rng_counter_increment();
+        force1 += dpd_pair_force(p, part_rep, ia_params, dist_vec, dist,
+                                 dist * dist, integrator_counter.value());
       }
 #endif
     } else if (m_penetrable && (dist <= 0)) {
@@ -87,9 +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);
-          // Additional use of DPD here requires counter increase
-          dpd_rng_counter_increment();
+                                   dist * dist, integrator_counter.value());
         }
 #endif
       }

src/core/dpd.cpp

@@ -46,10 +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<RNGSalt::SALT_DPD>(
-      dpd.rng_get(), (pid1 < pid2) ? pid2 : pid1, (pid1 < pid2) ? pid1 : pid2);
+  return Random::noise_uniform<RNGSalt::SALT_DPD>(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,
@@ -116,15 +117,15 @@ Vector3d dpd_pair_force(DPDParameters const &params, 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 {};
   }
 
   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);

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

src/core/forces.cpp

@@ -49,7 +49,8 @@
 
 ActorList forceActors;
 
-void init_forces(const ParticleRange &particles) {
+void init_forces(const ParticleRange &particles, uint64_t counter,
+                 double time_step) {
   ESPRESSO_PROFILER_CXX_MARK_FUNCTION;
   /* The force initialization depends on the used thermostat and the
      thermodynamic ensemble */
@@ -63,7 +64,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, time_step);
   }
 
   /* initialize ghost forces with zero
@@ -80,7 +81,8 @@ void init_forces_ghosts(const ParticleRange &particles) {
   }
 }
 
-void force_calc(CellStructure &cell_structure) {
+void force_calc(CellStructure &cell_structure, uint64_t const counter,
+                double time_step) {
   ESPRESSO_PROFILER_CXX_MARK_FUNCTION;
 
   espressoSystemInterface.update();
@@ -94,7 +96,7 @@ void force_calc(CellStructure &cell_structure) {
 #ifdef ELECTROSTATICS
   iccp3m_iteration(particles, cell_structure.ghost_particles());
 #endif
-  init_forces(particles);
+  init_forces(particles, counter, time_step);
 
   for (auto &forceActor : forceActors) {
     forceActor->computeForces(espressoSystemInterface);
@@ -119,8 +121,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);

src/core/forces.hpp

@@ -38,15 +38,10 @@
 
 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,
+                 double time_step);
 
 /** Set forces of all ghosts to zero */
 void init_forces_ghosts(const ParticleRange &particles);
@@ -61,10 +56,10 @@ 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, uint64_t counter,
+                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,29 @@ 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,
+                                      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, counter, time_step),
           p.p.rotation ? convert_vector_body_to_space(
-                             p, friction_thermo_langevin_rotation(langevin, p))
+                             p, friction_thermo_langevin_rotation(
+                                    langevin, p, counter, time_step))
                        : Utils::Vector3d{}};
 #else
-  return friction_thermo_langevin(langevin, p);
+  return friction_thermo_langevin(langevin, p, counter, 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,
+                                               uint64_t counter,
+                                               double time_step) {
+  return thermostat_force(part, counter, time_step) + external_force(part);
 }
 
 inline Utils::Vector3d calc_non_bonded_pair_force_parts(
@@ -234,7 +238,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 +293,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
 

src/core/grid_based_algorithms/lb.cpp

@@ -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 {};
         };