Encapsulate thermostats

API change: thermalized bond random number generator seed moved to
the thermostat class to avoid accidentally modifying the global seed.

doc/sphinx/integration.rst

@@ -401,15 +401,14 @@ To add a thermostat, call the appropriate setter::
 The different thermostats available in |es| will be described in the following
 subsections.
 
-You may combine different thermostats at your own risk by turning them on
-one by one. The list of active thermostats can be cleared at any time with
+You may combine different thermostats by turning them on sequentially.
+Not all combinations of thermostats are sensible, though, and some
+integrators only work with a specific thermostat. The list of possible
+combinations of integrators and thermostats is hardcoded and automatically
+check against at the start of integration.
+Note that there is only one temperature for all thermostats.
+The list of active thermostats can be cleared at any time with
 :py:meth:`system.thermostat.turn_off() <espressomd.thermostat.Thermostat.turn_off>`.
-Not all combinations of thermostats are allowed, though (see
-:py:func:`espressomd.thermostat.AssertThermostatType` for details).
-Some integrators only work with a specific thermostat and throw an
-error otherwise. Note that there is only one temperature for all
-thermostats, although for some thermostats like the Langevin thermostat,
-particles can be assigned individual temperatures.
 
 Since |es| does not enforce a particular unit system, it cannot know about
 the current value of the Boltzmann constant. Therefore, when specifying

doc/sphinx/inter_bonded.rst

@@ -197,8 +197,9 @@ A thermalized bond can be instantiated via
     thermalized_bond = espressomd.interactions.ThermalizedBond(
         temp_com=<float>, gamma_com=<float>,
         temp_distance=<float>, gamma_distance=<float>,
-        r_cut=<float>, seed=<int>)
+        r_cut=<float>)
     system.bonded_inter.add(thermalized_bond)
+    system.thermostat.set_thermalized_bond(seed=<int>)
 
 This bond can be used to apply Langevin thermalization on the centre of mass
 and the distance of a particle pair.  Each thermostat can have its own

samples/dancing.py

@@ -49,6 +49,7 @@
 system.cell_system.skin = 0.4
 system.periodicity = [False, False, False]
 
+system.thermostat.set_stokesian(kT=0.)
 system.integrator.set_stokesian_dynamics(
     viscosity=1.0, radii={0: 1.0}, approximation_method=sd_method)
 

samples/drude_bmimpf6.py

@@ -258,10 +258,11 @@ def combination_rule_sigma(rule, sig1, sig2):
 
 if args.drude:
     print("-->Adding Drude related bonds")
+    system.thermostat.set_thermalized_bond(seed=123)
     thermalized_dist_bond = espressomd.interactions.ThermalizedBond(
         temp_com=temperature_com, gamma_com=gamma_com,
         temp_distance=temperature_drude, gamma_distance=gamma_drude,
-        r_cut=min(lj_sigmas.values()) * 0.5, seed=123)
+        r_cut=min(lj_sigmas.values()) * 0.5)
     harmonic_bond = espressomd.interactions.HarmonicBond(
         k=k_drude, r_0=0.0, r_cut=1.0)
     system.bonded_inter.add(thermalized_dist_bond)

samples/load_checkpoint.py

@@ -60,10 +60,6 @@
 print("\n### system.part test ###")
 print(f"system.part.all().pos = {system.part.all().pos}")
 
-# test "system.thermostat"
-print("\n### system.thermostat test ###")
-print(f"system.thermostat.get_state() = {system.thermostat.get_state()}")
-
 # test "p3m"
 print("\n### p3m test ###")
 print(f"p3m.get_params() = {p3m.get_params()}")

src/core/CMakeLists.txt

@@ -32,7 +32,6 @@ add_library(
   errorhandling.cpp
   forces.cpp
   ghosts.cpp
-  global_ghost_flags.cpp
   immersed_boundaries.cpp
   integrate.cpp
   npt.cpp

src/core/PropagationMode.hpp

@@ -20,9 +20,7 @@
 #pragma once
 
 namespace PropagationMode {
-/**
- * @brief Flags to create bitmasks for propagation modes.
- */
+/** @brief Flags to create bitmasks for propagation modes. */
 enum PropagationMode : int {
   NONE = 0,
   SYSTEM_DEFAULT = 1 << 0,
@@ -42,11 +40,23 @@ enum PropagationMode : int {
 };
 } // namespace PropagationMode
 
-/** \name Integrator switches */
+/** @brief Integrator identifier. */
 enum IntegratorSwitch : int {
   INTEG_METHOD_NPT_ISO = 0,
   INTEG_METHOD_NVT = 1,
   INTEG_METHOD_STEEPEST_DESCENT = 2,
   INTEG_METHOD_BD = 3,
-  INTEG_METHOD_SD = 7,
+  INTEG_METHOD_SD = 4,
+};
+
+/** @brief Thermostat flags. */
+enum ThermostatFlags : int {
+  THERMO_OFF = 0,
+  THERMO_LANGEVIN = 1 << 0,
+  THERMO_BROWNIAN = 1 << 1,
+  THERMO_NPT_ISO = 1 << 2,
+  THERMO_LB = 1 << 3,
+  THERMO_SD = 1 << 4,
+  THERMO_DPD = 1 << 5,
+  THERMO_BOND = 1 << 6,
 };

src/core/bonded_interactions/thermalized_bond_kernel.hpp

@@ -25,6 +25,7 @@
 
 #include "Particle.hpp"
 #include "random.hpp"
+#include "system/System.hpp"
 #include "thermostat.hpp"
 
 #include <utils/Vector.hpp>
@@ -54,8 +55,9 @@ ThermalizedBond::forces(Particle const &p1, Particle const &p2,
   auto const sqrt_mass_red = sqrt(p1.mass() * p2.mass() / mass_tot);
   auto const com_vel = mass_tot_inv * (p1.mass() * p1.v() + p2.mass() * p2.v());
   auto const dist_vel = p2.v() - p1.v();
+  auto const &thermalized_bond =
+      *::System::get_system().thermostat->thermalized_bond;
 
-  extern ThermalizedBondThermostat thermalized_bond;
   Utils::Vector3d force1{};
   Utils::Vector3d force2{};
   auto const noise = Random::noise_uniform<RNGSalt::THERMALIZED_BOND>(

src/core/cell_system/CellStructure.cpp

@@ -171,33 +171,33 @@ unsigned map_data_parts(unsigned data_parts) {
 
   /* clang-format off */
   return GHOSTTRANS_NONE
-         | ((DATA_PART_PROPERTIES & data_parts) ? GHOSTTRANS_PROPRTS : 0u)
-         | ((DATA_PART_POSITION & data_parts) ? GHOSTTRANS_POSITION : 0u)
-         | ((DATA_PART_MOMENTUM & data_parts) ? GHOSTTRANS_MOMENTUM : 0u)
-         | ((DATA_PART_FORCE & data_parts) ? GHOSTTRANS_FORCE : 0u)
+         | ((data_parts & DATA_PART_PROPERTIES) ? GHOSTTRANS_PROPRTS : 0u)
+         | ((data_parts & DATA_PART_POSITION) ? GHOSTTRANS_POSITION : 0u)
+         | ((data_parts & DATA_PART_MOMENTUM) ? GHOSTTRANS_MOMENTUM : 0u)
+         | ((data_parts & DATA_PART_FORCE) ? GHOSTTRANS_FORCE : 0u)
 #ifdef BOND_CONSTRAINT
-         | ((DATA_PART_RATTLE & data_parts) ? GHOSTTRANS_RATTLE : 0u)
+         | ((data_parts & DATA_PART_RATTLE) ? GHOSTTRANS_RATTLE : 0u)
 #endif
-         | ((DATA_PART_BONDS & data_parts) ? GHOSTTRANS_BONDS : 0u);
+         | ((data_parts & DATA_PART_BONDS) ? GHOSTTRANS_BONDS : 0u);
   /* clang-format on */
 }
 
 void CellStructure::ghosts_count() {
   ghost_communicator(decomposition().exchange_ghosts_comm(),
-                     GHOSTTRANS_PARTNUM);
+                     *get_system().box_geo, GHOSTTRANS_PARTNUM);
 }
 void CellStructure::ghosts_update(unsigned data_parts) {
   ghost_communicator(decomposition().exchange_ghosts_comm(),
-                     map_data_parts(data_parts));
+                     *get_system().box_geo, map_data_parts(data_parts));
 }
 void CellStructure::ghosts_reduce_forces() {
   ghost_communicator(decomposition().collect_ghost_force_comm(),
-                     GHOSTTRANS_FORCE);
+                     *get_system().box_geo, GHOSTTRANS_FORCE);
 }
 #ifdef BOND_CONSTRAINT
 void CellStructure::ghosts_reduce_rattle_correction() {
   ghost_communicator(decomposition().collect_ghost_force_comm(),
-                     GHOSTTRANS_RATTLE);
+                     *get_system().box_geo, GHOSTTRANS_RATTLE);
 }
 #endif
 
@@ -265,8 +265,9 @@ void CellStructure::set_hybrid_decomposition(double cutoff_regular,
   auto &local_geo = *system.local_geo;
   auto const &box_geo = *system.box_geo;
   set_particle_decomposition(std::make_unique<HybridDecomposition>(
-      ::comm_cart, cutoff_regular, m_verlet_skin, box_geo, local_geo,
-      n_square_types));
+      ::comm_cart, cutoff_regular, m_verlet_skin,
+      [&system]() { return system.get_global_ghost_flags(); }, box_geo,
+      local_geo, n_square_types));
   m_type = CellStructureType::HYBRID;
   local_geo.set_cell_structure_type(m_type);
   system.on_cell_structure_change();

src/core/cell_system/HybridDecomposition.cpp

@@ -26,7 +26,6 @@
 
 #include "BoxGeometry.hpp"
 #include "LocalBox.hpp"
-#include "global_ghost_flags.hpp"
 
 #include <utils/Vector.hpp>
 #include <utils/mpi/sendrecv.hpp>
@@ -36,20 +35,23 @@
 
 #include <algorithm>
 #include <cstddef>
+#include <functional>
 #include <iterator>
 #include <set>
 #include <utility>
 
 HybridDecomposition::HybridDecomposition(boost::mpi::communicator comm,
                                          double cutoff_regular, double skin,
+                                         std::function<bool()> get_ghost_flags,
                                          BoxGeometry const &box_geo,
                                          LocalBox const &local_box,
                                          std::set<int> n_square_types)
     : m_comm(std::move(comm)), m_box(box_geo), m_cutoff_regular(cutoff_regular),
       m_regular_decomposition(RegularDecomposition(
           m_comm, cutoff_regular + skin, m_box, local_box)),
       m_n_square(AtomDecomposition(m_comm, m_box)),
-      m_n_square_types(std::move(n_square_types)) {
+      m_n_square_types(std::move(n_square_types)),
+      m_get_global_ghost_flags(std::move(get_ghost_flags)) {
 
   /* Vector containing cells of both child decompositions */
   m_local_cells = m_regular_decomposition.get_local_cells();
@@ -155,11 +157,11 @@ void HybridDecomposition::resort(bool global,
   m_n_square.resort(global, diff);
 
   /* basically do CellStructure::ghost_count() */
-  ghost_communicator(exchange_ghosts_comm(), GHOSTTRANS_PARTNUM);
+  ghost_communicator(exchange_ghosts_comm(), m_box, GHOSTTRANS_PARTNUM);
 
   /* basically do CellStructure::ghost_update(unsigned data_parts) */
-  ghost_communicator(exchange_ghosts_comm(),
-                     map_data_parts(global_ghost_flags()));
+  ghost_communicator(exchange_ghosts_comm(), m_box,
+                     map_data_parts(m_get_global_ghost_flags()));
 }
 
 std::size_t HybridDecomposition::count_particles(

src/core/cell_system/HybridDecomposition.hpp

@@ -39,6 +39,7 @@
 #include <boost/optional.hpp>
 
 #include <cstddef>
+#include <functional>
 #include <set>
 #include <utility>
 #include <vector>
@@ -69,14 +70,17 @@ class HybridDecomposition : public ParticleDecomposition {
   /** Set containing the types that should be handled using n_square */
   std::set<int> const m_n_square_types;
 
+  std::function<bool()> m_get_global_ghost_flags;
+
   bool is_n_square_type(int type_id) const {
     return (m_n_square_types.find(type_id) != m_n_square_types.end());
   }
 
 public:
   HybridDecomposition(boost::mpi::communicator comm, double cutoff_regular,
-                      double skin, BoxGeometry const &box_geo,
-                      LocalBox const &local_box, std::set<int> n_square_types);
+                      double skin, std::function<bool()> get_ghost_flags,
+                      BoxGeometry const &box_geo, LocalBox const &local_box,
+                      std::set<int> n_square_types);
 
   auto get_cell_grid() const { return m_regular_decomposition.cell_grid; }
 

src/core/constraints/ShapeBasedConstraint.cpp

@@ -101,11 +101,12 @@ ParticleForce ShapeBasedConstraint::force(Particle const &p,
            calc_non_central_force(p, part_rep, ia_params, dist_vec, dist);
 
 #ifdef DPD
-      if (thermo_switch & THERMO_DPD) {
-        dpd_force =
-            dpd_pair_force(p, part_rep, ia_params, dist_vec, dist, dist * dist);
+      if (m_system.thermostat->thermo_switch & THERMO_DPD) {
+        dpd_force = dpd_pair_force(p, part_rep, *m_system.thermostat->dpd,
+                                   *m_system.box_geo, ia_params, dist_vec, dist,
+                                   dist * dist);
         // Additional use of DPD here requires counter increase
-        dpd.rng_increment();
+        m_system.thermostat->dpd->rng_increment();
       }
 #endif
     } else if (m_penetrable && (dist <= 0)) {
@@ -117,11 +118,12 @@ ParticleForce ShapeBasedConstraint::force(Particle const &p,
             calc_non_central_force(p, part_rep, ia_params, dist_vec, -dist);
 
 #ifdef DPD
-        if (thermo_switch & THERMO_DPD) {
-          dpd_force = dpd_pair_force(p, part_rep, ia_params, dist_vec, dist,
-                                     dist * dist);
+        if (m_system.thermostat->thermo_switch & THERMO_DPD) {
+          dpd_force = dpd_pair_force(p, part_rep, *m_system.thermostat->dpd,
+                                     *m_system.box_geo, ia_params, dist_vec,
+                                     dist, dist * dist);
           // Additional use of DPD here requires counter increase
-          dpd.rng_increment();
+          m_system.thermostat->dpd->rng_increment();
         }
 #endif
       }

src/core/dpd.cpp

@@ -55,7 +55,7 @@
  *  3. Two particle IDs (order-independent, decorrelates particles, gets rid of
  *     seed-per-node)
  */
-Utils::Vector3d dpd_noise(int pid1, int pid2) {
+Utils::Vector3d dpd_noise(DPDThermostat const &dpd, int pid1, int pid2) {
   return Random::noise_uniform<RNGSalt::SALT_DPD>(
       dpd.rng_counter(), dpd.rng_seed(), (pid1 < pid2) ? pid2 : pid1,
       (pid1 < pid2) ? pid1 : pid2);
@@ -99,21 +99,19 @@ Utils::Vector3d dpd_pair_force(DPDParameters const &params,
   return {};
 }
 
-Utils::Vector3d dpd_pair_force(Particle const &p1, Particle const &p2,
-                               IA_parameters const &ia_params,
-                               Utils::Vector3d const &d, double dist,
-                               double dist2) {
+Utils::Vector3d
+dpd_pair_force(Particle const &p1, Particle const &p2, DPDThermostat const &dpd,
+               BoxGeometry const &box_geo, IA_parameters const &ia_params,
+               Utils::Vector3d const &d, double dist, double dist2) {
   if (ia_params.dpd.radial.cutoff <= 0.0 && ia_params.dpd.trans.cutoff <= 0.0) {
     return {};
   }
 
-  auto const &box_geo = *System::get_system().box_geo;
-
   auto const v21 =
       box_geo.velocity_difference(p1.pos(), p2.pos(), p1.v(), p2.v());
   auto const noise_vec =
       (ia_params.dpd.radial.pref > 0.0 || ia_params.dpd.trans.pref > 0.0)
-          ? dpd_noise(p1.id(), p2.id())
+          ? dpd_noise(dpd, p1.id(), p2.id())
           : Utils::Vector3d{};
 
   auto const f_r = dpd_pair_force(ia_params.dpd.radial, v21, dist, noise_vec);

src/core/dpd.hpp

@@ -18,8 +18,9 @@
  * You should have received a copy of the GNU General Public License
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
-#ifndef ESPRESSO_SRC_CORE_DPD_HPP
-#define ESPRESSO_SRC_CORE_DPD_HPP
+
+#pragma once
+
 /** \file
  *  Routines to use DPD as thermostat or pair force @cite soddemann03a
  *
@@ -30,7 +31,9 @@
 
 #ifdef DPD
 
+#include "BoxGeometry.hpp"
 #include "Particle.hpp"
+#include "thermostat.hpp"
 
 #include <utils/Vector.hpp>
 
@@ -43,11 +46,10 @@ struct IA_parameters;
 
 void dpd_init(double kT, double time_step);
 
-Utils::Vector3d dpd_pair_force(Particle const &p1, Particle const &p2,
-                               IA_parameters const &ia_params,
-                               Utils::Vector3d const &d, double dist,
-                               double dist2);
+Utils::Vector3d
+dpd_pair_force(Particle const &p1, Particle const &p2, DPDThermostat const &dpd,
+               BoxGeometry const &box_geo, IA_parameters const &ia_params,
+               Utils::Vector3d const &d, double dist, double dist2);
 Utils::Vector9d dpd_stress(boost::mpi::communicator const &comm);
 
 #endif // DPD
-#endif

src/core/ek/EKNone.hpp

@@ -32,6 +32,7 @@ struct EKNone {
   void propagate() { throw NoEKActive{}; }
   double get_tau() const { throw NoEKActive{}; }
   void veto_time_step(double) const { throw NoEKActive{}; }
+  void veto_kT(double) const { throw NoEKActive{}; }
   void sanity_checks(System::System const &) const { throw NoEKActive{}; }
   void on_cell_structure_change() const { throw NoEKActive{}; }
   void on_boxl_change() const { throw NoEKActive{}; }