WIP

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

maintainer/CI/build_cmake.sh

@@ -404,7 +404,7 @@ if [ "${with_coverage}" = true ] || [ "${with_coverage_python}" = true ]; then
         lcov --gcov-tool "${GCOV:-gcov}" -q --remove coverage.info '/usr/*' --output-file coverage.info # filter out system
         lcov --gcov-tool "${GCOV:-gcov}" -q --remove coverage.info '*/doc/*' --output-file coverage.info # filter out docs
         if [ -d _deps/ ]; then
-          lcov --gcov-tool "${GCOV:-gcov}" -q --remove coverage.info $(realpath _deps/)'/*' --output-file coverage.info # filter out docs
+          lcov --gcov-tool "${GCOV:-gcov}" -q --remove coverage.info $(realpath _deps/)'/*' --output-file coverage.info # filter out external projects
         fi
     fi
     if [ "${with_coverage_python}" = true ]; then

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/CMakeLists.txt

@@ -23,7 +23,4 @@ target_sources(
           ${CMAKE_CURRENT_SOURCE_DIR}/angle_cossquare.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/bonded_interaction_data.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/bonded_tab.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/fene.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/rigid_bond.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/thermalized_bond.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/thermalized_bond_utils.cpp)
+          ${CMAKE_CURRENT_SOURCE_DIR}/fene.cpp)

src/core/bonded_interactions/bonded_interaction_data.cpp

@@ -17,7 +17,9 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 #include "bonded_interaction_data.hpp"
+#include "rigid_bond.hpp"
 #include "system/System.hpp"
+#include "thermalized_bond.hpp"
 
 #include <boost/range/numeric.hpp>
 #include <boost/variant.hpp>
@@ -60,6 +62,22 @@ double maximal_cutoff_bonded() {
 
 void BondedInteractionsMap::on_ia_change() {
   if (System::is_system_set()) {
-    System::get_system().on_short_range_ia_change();
+    auto &system = System::get_system();
+    system.on_short_range_ia_change();
+    system.on_thermostat_param_change(); // thermalized bonds
+  }
+  n_thermalized_bonds = 0;
+#ifdef BOND_CONSTRAINT
+  n_rigid_bonds = 0;
+#endif
+  for (auto &kv : *this) {
+    if (boost::get<ThermalizedBond>(&(*kv.second)) != nullptr) {
+      ++n_thermalized_bonds;
+    }
+#ifdef BOND_CONSTRAINT
+    if (boost::get<RigidBond>(&(*kv.second)) != nullptr) {
+      ++n_rigid_bonds;
+    }
+#endif
   }
 }

src/core/bonded_interactions/bonded_interaction_data.hpp

@@ -143,10 +143,18 @@ class BondedInteractionsMap {
   auto get_zero_based_type(int bond_id) const {
     return contains(bond_id) ? at(bond_id)->which() : 0;
   }
+  auto get_n_thermalized_bonds() const { return n_thermalized_bonds; }
+#ifdef BOND_CONSTRAINT
+  auto get_n_rigid_bonds() const { return n_rigid_bonds; }
+#endif
 
 private:
   container_type m_params = {};
   key_type next_key = static_cast<key_type>(0);
+  int n_thermalized_bonds = 0;
+#ifdef BOND_CONSTRAINT
+  int n_rigid_bonds = 0;
+#endif
   void on_ia_change();
 };
 

src/core/bonded_interactions/rigid_bond.hpp

@@ -21,19 +21,13 @@
 #ifndef RIGID_BOND_HPP
 #define RIGID_BOND_HPP
 /** \file
- *  Definition of the rigid bond data type. It is utilized by the
- *  Rattle algorithm.
- *
- *  Implementation in \ref rigid_bond.cpp.
+ *  Definition of the rigid bond data type for the Rattle algorithm.
  */
 
 #include <utils/Vector.hpp>
 
 #include <cmath>
 
-/** Number of rigid bonds. */
-extern int n_rigidbonds;
-
 /** Parameters for the rigid_bond/SHAKE/RATTLE ALGORITHM */
 struct RigidBond {
   /** Square of the length of Constrained Bond */
@@ -51,7 +45,11 @@ struct RigidBond {
 
   static constexpr int num = 1;
 
-  RigidBond(double d, double p_tol, double v_tol);
+  RigidBond(double d, double p_tol, double v_tol) {
+    this->d2 = d * d;
+    this->p_tol = 2.0 * p_tol;
+    this->v_tol = v_tol;
+  }
 
 private:
   friend boost::serialization::access;

src/core/bonded_interactions/thermalized_bond.hpp

@@ -1,4 +1,3 @@
-
 /*
  * Copyright (C) 2010-2022 The ESPResSo project
  * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
@@ -24,19 +23,15 @@
 #define THERMALIZED_DIST_H
 /** \file
  *  Routines to thermalize the center of mass and distance of a particle pair.
- *
- *  Implementation in \ref thermalized_bond.cpp.
  */
 
-/** number of thermalized bonds */
-extern int n_thermalized_bonds;
-
 #include "Particle.hpp"
 
 #include <utils/Vector.hpp>
 
 #include <boost/optional.hpp>
 
+#include <cmath>
 #include <tuple>
 
 /** Parameters for Thermalized bond */
@@ -56,7 +51,26 @@ struct ThermalizedBond {
   static constexpr int num = 1;
 
   ThermalizedBond(double temp_com, double gamma_com, double temp_distance,
-                  double gamma_distance, double r_cut);
+                  double gamma_distance, double r_cut) {
+    this->temp_com = temp_com;
+    this->gamma_com = gamma_com;
+    this->temp_distance = temp_distance;
+    this->gamma_distance = gamma_distance;
+    this->r_cut = r_cut;
+
+    pref1_com = -1.;
+    pref2_com = -1.;
+    pref1_dist = -1.;
+    pref2_dist = -1.;
+  }
+
+  void recalc_prefactors(double time_step) {
+    pref1_com = gamma_com;
+    pref2_com = std::sqrt(24. * gamma_com / time_step * temp_com);
+    pref1_dist = gamma_distance;
+    pref2_dist = std::sqrt(24. * gamma_distance / time_step * temp_distance);
+  }
+
   boost::optional<std::tuple<Utils::Vector3d, Utils::Vector3d>>
   forces(Particle const &p1, Particle const &p2,
          Utils::Vector3d const &dx) const;

src/core/bonded_interactions/thermalized_bond_kernel.hpp

@@ -26,6 +26,7 @@
 
 #include "Particle.hpp"
 #include "random.hpp"
+#include "system/System.hpp"
 #include "thermostat.hpp"
 
 #include <utils/Vector.hpp>
@@ -55,8 +56,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>(