Fix thermostats

Properly handle null values for kT and friction coefficients.
Make friction coefficients required parameters.

doc/tutorials/electrokinetics/electrokinetics.ipynb

@@ -1042,7 +1042,7 @@
     "                lattice=lattice, density=DENSITY_FLUID, kinematic_viscosity=VISCOSITY_KINEMATIC,\n",
     "                tau=TAU, ext_force_density=EXT_FORCE_DENSITY, kT=KT, seed=42)\n",
     "system.lb = lbf\n",
-    "system.thermostat.set_lb(LB_fluid=lbf, seed=42)\n",
+    "system.thermostat.set_lb(LB_fluid=lbf, seed=42, gamma=0.)\n",
     "eksolver = espressomd.electrokinetics.EKNone(lattice=lattice)\n",
     "system.ekcontainer = espressomd.electrokinetics.EKContainer(tau=TAU, solver=eksolver)"
    ]

samples/lb_profile.py

@@ -45,7 +45,6 @@
     agrid=1.0, density=1.0, kinematic_viscosity=1.0, tau=0.01,
     ext_force_density=[0, 0, 0.15], kT=0.0)
 system.lb = lb_fluid
-system.thermostat.set_lb(LB_fluid=lb_fluid, seed=23)
 ctp = espressomd.math.CylindricalTransformationParameters(
     center=[5.0, 5.0, 0.0],
     axis=[0, 0, 1],

src/core/lb/particle_coupling.cpp

@@ -291,18 +291,23 @@ void ParticleCoupling::kernel(std::vector<Particle *> const &particles) {
 #endif
     Utils::Vector3d force_on_particle = {};
     if (coupling_mode == particle_force) {
-      auto &v_fluid = *it_interpolated_velocities;
-      if (m_box_geo.type() == BoxType::LEES_EDWARDS) {
-        // Account for the case where the interpolated velocity has been read
-        // from a ghost of the particle across the LE boundary (or vice verssa)
-        // Then the particle velocity is shifted by +,- the LE shear velocity
-        auto const vel_correction = lees_edwards_vel_shift(
-            *it_positions_velocity_coupling, p.pos(), m_box_geo);
-        v_fluid += vel_correction;
+#ifndef THERMOSTAT_PER_PARTICLE
+      if (m_thermostat.gamma > 0.)
+#endif
+      {
+        auto &v_fluid = *it_interpolated_velocities;
+        if (m_box_geo.type() == BoxType::LEES_EDWARDS) {
+          // Account for the case where the interpolated velocity has been read
+          // from a ghost of the particle across the LE boundary (or vice versa)
+          // Then the particle velocity is shifted by +,- the LE shear velocity
+          auto const vel_correction = lees_edwards_vel_shift(
+              *it_positions_velocity_coupling, p.pos(), m_box_geo);
+          v_fluid += vel_correction;
+        }
+        auto const drag_force = lb_drag_force(p, m_thermostat.gamma, v_fluid);
+        auto const random_force = get_noise_term(p);
+        force_on_particle = drag_force + random_force;
       }
-      auto const drag_force = lb_drag_force(p, m_thermostat.gamma, v_fluid);
-      auto const random_force = get_noise_term(p);
-      force_on_particle = drag_force + random_force;
       ++it_interpolated_velocities;
       ++it_positions_velocity_coupling;
     }

src/core/thermostat.hpp

@@ -97,6 +97,9 @@ inline bool are_kT_equal(double old_kT, double new_kT) {
   }
   auto const large_kT = (old_kT > new_kT) ? old_kT : new_kT;
   auto const small_kT = (old_kT > new_kT) ? new_kT : old_kT;
+  if (small_kT == 0.) {
+    return false;
+  }
   return (large_kT / small_kT - 1. < relative_tolerance);
 }
 } // namespace Thermostat

src/script_interface/thermostat/thermostat.hpp

@@ -33,6 +33,7 @@
 #include <cassert>
 #include <limits>
 #include <memory>
+#include <span>
 #include <stdexcept>
 #include <string>
 
@@ -92,7 +93,7 @@ class Interface : public AutoParameters<Interface<CoreClass>, System::Leaf> {
   }
 
   virtual bool invalid_rng_state(VariantMap const &params) const {
-    return (not params.count("seed") or is_none(params.at("seed"))) and
+    return (not params.contains("seed") or is_none(params.at("seed"))) and
            is_seed_required();
   }
 
@@ -101,14 +102,15 @@ class Interface : public AutoParameters<Interface<CoreClass>, System::Leaf> {
   get_member_handle(::Thermostat::Thermostat &thermostat) = 0;
 
   void set_new_parameters(VariantMap const &params) {
-    if (params.count("__check_rng_state") and invalid_rng_state(params)) {
-      context()->parallel_try_catch([]() {
+    context()->parallel_try_catch([&]() {
+      if (params.contains("__check_rng_state") and invalid_rng_state(params)) {
         throw std::invalid_argument("Parameter 'seed' is needed on first "
                                     "activation of the thermostat");
-      });
-    }
+      }
+      check_required_parameters(params);
+    });
     for (auto const &key : get_parameter_insertion_order()) {
-      if (params.count(key)) {
+      if (params.contains(key)) {
         auto const &v = params.at(key);
         if (key == "is_active") {
           is_active = get_value<bool>(v);
@@ -119,6 +121,17 @@ class Interface : public AutoParameters<Interface<CoreClass>, System::Leaf> {
     }
   }
 
+  virtual std::span<std::string_view const> get_required_parameters() const = 0;
+
+  void check_required_parameters(VariantMap const &params) const {
+    for (auto const &required : get_required_parameters()) {
+      auto name = std::string(required);
+      if (not params.contains(name)) {
+        throw std::runtime_error("Parameter '" + name + "' is missing");
+      }
+    }
+  }
+
 protected:
   template <typename T>
   auto make_autoparameter(T CoreThermostat::*member, char const *name) {
@@ -169,7 +182,7 @@ class Interface : public AutoParameters<Interface<CoreClass>, System::Leaf> {
     auto params = parameters;
     if (not is_seed_required()) {
       for (auto key : {std::string("seed"), std::string("philox_counter")}) {
-        if (params.count(key) == 0ul) {
+        if (not params.contains(key)) {
           params[key] = get_parameter(key);
         }
       }
@@ -216,7 +229,7 @@ class Interface : public AutoParameters<Interface<CoreClass>, System::Leaf> {
   }
 
   virtual std::optional<double> extract_kT(VariantMap const &params) const {
-    if (params.count("kT")) {
+    if (params.contains("kT")) {
       auto const value = get_value<double>(params, "kT");
       sanity_checks_positive(value, "kT");
       return value;
@@ -301,6 +314,11 @@ class Langevin : public Interface<::LangevinThermostat> {
     return thermostat.langevin;
   }
 
+  std::span<std::string_view const> get_required_parameters() const override {
+    static constexpr std::array names{std::string_view("gamma")};
+    return names;
+  }
+
 public:
   Langevin() {
     add_parameters({
@@ -319,7 +337,7 @@ class Langevin : public Interface<::LangevinThermostat> {
         Interface<::LangevinThermostat>::extend_parameters(parameters);
 #ifdef ROTATION
     // If gamma_rotation is not set explicitly, use the translational one.
-    if (params.count("gamma_rotation") == 0ul and params.count("gamma")) {
+    if (not params.contains("gamma_rotation") and params.contains("gamma")) {
       params["gamma_rotation"] = params.at("gamma");
     }
 #endif // ROTATION
@@ -333,6 +351,11 @@ class Brownian : public Interface<::BrownianThermostat> {
     return thermostat.brownian;
   }
 
+  std::span<std::string_view const> get_required_parameters() const override {
+    static constexpr std::array names{std::string_view("gamma")};
+    return names;
+  }
+
 public:
   Brownian() {
     add_parameters({
@@ -351,7 +374,7 @@ class Brownian : public Interface<::BrownianThermostat> {
         Interface<::BrownianThermostat>::extend_parameters(parameters);
 #ifdef ROTATION
     // If gamma_rotation is not set explicitly, use the translational one.
-    if (params.count("gamma_rotation") == 0ul and params.count("gamma")) {
+    if (not params.contains("gamma_rotation") and params.contains("gamma")) {
       params["gamma_rotation"] = params.at("gamma");
     }
 #endif // ROTATION
@@ -366,6 +389,12 @@ class IsotropicNpt : public Interface<::IsotropicNptThermostat> {
     return thermostat.npt_iso;
   }
 
+  std::span<std::string_view const> get_required_parameters() const override {
+    static constexpr std::array names{std::string_view("gamma0"),
+                                      std::string_view("gammav")};
+    return names;
+  }
+
 public:
   IsotropicNpt() {
     add_parameters({
@@ -418,10 +447,15 @@ class LBThermostat : public Interface<::LBThermostat> {
 
 protected:
   bool invalid_rng_state(VariantMap const &params) const override {
-    return (not params.count("seed") or is_none(params.at("seed"))) and
-           params.count("__global_kT") and is_seed_required() and
+    return (not params.contains("seed") or is_none(params.at("seed"))) and
+           params.contains("__global_kT") and is_seed_required() and
            get_value<double>(params, "__global_kT") != 0.;
   }
+
+  std::span<std::string_view const> get_required_parameters() const override {
+    static constexpr std::array names{std::string_view("gamma")};
+    return names;
+  }
 };
 #endif // WALBERLA
 
@@ -432,6 +466,10 @@ class DPDThermostat : public Interface<::DPDThermostat> {
     return thermostat.dpd;
   }
 
+  std::span<std::string_view const> get_required_parameters() const override {
+    return {};
+  }
+
 public:
   ::ThermostatFlags get_thermo_flag() const final { return THERMO_DPD; }
 };
@@ -444,6 +482,10 @@ class Stokesian : public Interface<::StokesianThermostat> {
     return thermostat.stokesian;
   }
 
+  std::span<std::string_view const> get_required_parameters() const override {
+    return {};
+  }
+
 public:
   ::ThermostatFlags get_thermo_flag() const final { return THERMO_SD; }
 };
@@ -455,6 +497,10 @@ class ThermalizedBond : public Interface<::ThermalizedBondThermostat> {
     return thermostat.thermalized_bond;
   }
 
+  std::span<std::string_view const> get_required_parameters() const override {
+    return {};
+  }
+
 public:
   ::ThermostatFlags get_thermo_flag() const final { return THERMO_BOND; }
 };

testsuite/python/lb.py

@@ -54,7 +54,6 @@ class LBTest:
     system.cell_system.skin = 1.0
     if espressomd.gpu_available():
         system.cuda_init_handle.call_method("set_device_id_per_rank")
-    interpolation = False
     n_nodes = system.cell_system.get_state()["n_nodes"]
 
     def setUp(self):
@@ -75,15 +74,13 @@ def test_properties(self):
         # activated actor
         lbf = self.lb_class(kT=1.0, seed=42, **self.params, **self.lb_params)
         self.system.lb = lbf
-        self.system.thermostat.set_lb(LB_fluid=lbf, seed=1)
         self.assertTrue(lbf.is_active)
         self.check_properties(lbf)
         self.system.lb = None
 
         # deactivated actor
         lbf = self.lb_class(kT=1.0, seed=42, **self.params, **self.lb_params)
         self.system.lb = lbf
-        self.system.thermostat.set_lb(LB_fluid=lbf, seed=1)
         self.system.lb = None
         self.assertFalse(lbf.is_active)
         self.check_properties(lbf)
@@ -383,7 +380,7 @@ def test_pressure_tensor_observable(self):
         lbf = self.lb_class(kT=1., seed=1, ext_force_density=[0, 0, 0],
                             **self.params, **self.lb_params)
         system.lb = lbf
-        system.thermostat.set_lb(LB_fluid=lbf, seed=1)
+        system.thermostat.set_lb(LB_fluid=lbf, seed=1, gamma=1.)
         system.integrator.run(10)
         pressure_tensor = np.copy(
             np.mean(lbf[:, :, :].pressure_tensor, axis=(0, 1, 2)))

testsuite/python/lb_interpolation.py

@@ -128,26 +128,33 @@ def test_mach_limit_check(self):
         self.assertIsNone(self.lbf[0, 0, 0].boundary)
 
     def test_interpolated_force(self):
+        def sample(kernel):
+            for i, j, k in itertools.product(range(3), range(3), range(3)):
+                system.lb[:, :, :].velocity = lb_vel
+                p.pos = (
+                    (i + 0.5) * AGRID,
+                    (j + 0.5) * AGRID,
+                    (k + 0.5) * AGRID)
+                p.v = [0., 0., 0.]
+                system.integrator.run(1)
+                np.testing.assert_allclose(np.copy(p.f), kernel(i, j, k))
         system = self.system
-        system.thermostat.set_lb(LB_fluid=system.lb, seed=42, gamma=1.)
         system.integrator.run(1)
         lb_vel = np.zeros([12, 12, 12, 3], dtype=float)
-        for i in range(12):
-            for j in range(12):
-                for k in range(12):
-                    lb_vel[i, j, k] = 1e-3 * np.array([i + 1, j + 1, k + 1])
+        for i, j, k in itertools.product(range(12), range(12), range(12)):
+            lb_vel[i, j, k] = 1e-3 * np.array([i + 1, j + 1, k + 1])
         p = system.part.add(pos=3 * [1.5])
-        for i in range(3):
-            for j in range(3):
-                for k in range(3):
-                    system.lb[:, :, :].velocity = lb_vel
-                    p.pos = (
-                        (i + 0.5) * AGRID,
-                        (j + 0.5) * AGRID,
-                        (k + 0.5) * AGRID)
-                    p.v = [0., 0., 0.]
-                    system.integrator.run(1)
-                    np.testing.assert_allclose(np.copy(p.f), lb_vel[i, j, k])
+        # enable particle coupling
+        system.thermostat.set_lb(LB_fluid=system.lb, seed=42, gamma=1.)
+        sample(kernel=lambda i, j, k: lb_vel[i, j, k])
+        # disable particle coupling
+        system.thermostat.set_lb(LB_fluid=system.lb, seed=42, gamma=0.)
+        sample(kernel=lambda i, j, k: [0., 0., 0.])
+        if espressomd.has_features("THERMOSTAT_PER_PARTICLE"):
+            # disable particle coupling globally & enable per-particle coupling
+            system.thermostat.set_lb(LB_fluid=system.lb, seed=42, gamma=0.)
+            p.gamma = 1.5
+            sample(kernel=lambda i, j, k: 1.5 * lb_vel[i, j, k])
 
 
 @utx.skipIfMissingFeatures(["WALBERLA"])

testsuite/python/lb_pressure_tensor.py

@@ -46,7 +46,6 @@ class TestLBPressureTensor:
     @classmethod
     def tearDownClass(cls):
         cls.system.lb = None
-        cls.system.thermostat.turn_off()
 
     @classmethod
     def setUpClass(cls):

testsuite/python/lb_thermostat.py

@@ -208,6 +208,8 @@ def test_exceptions(self):
         with self.assertRaisesRegex(RuntimeError, r"set_lb\(\) got an unexpected keyword argument 'act_on_virtual'"):
             self.system.thermostat.set_lb(
                 LB_fluid=self.lbf, act_on_virtual=False)
+        with self.assertRaisesRegex(RuntimeError, "Parameter 'gamma' is missing"):
+            self.system.thermostat.set_lb(LB_fluid=self.lbf)
         self.system.part.add(pos=[0., 0., 0.], gamma=[1., 2., 3.], id=2)
         with self.assertRaisesRegex(Exception, r"ERROR: anisotropic particle \(id 2\) coupled to LB"):
             self.system.integrator.run(1)