Merge #3070 #3080 #3096

3070: Thermalized bond philox r=jngrad a=KonradBreitsprecher

Fixes the TB part of #2683 

- Common philox counter for all thermalized bonds. Decorrelation via the bond partner IDs (not taking into account the bond ID). This was much easier to implement and is valid, as two (different) thermalized bonds on the same particle pair makes no sense.

3080: testsuite: add test for is_valid_type. r=jngrad a=KaiSzuttor

Currently it is possible to give ```nan``` as values to the particle properties because in python a ```nan``` is of type ```float```. This PR adds a check for ```None```, ```inf``` and ```nan``` to the ```is_valid_type``` and tests for the correct behavior.

3096: Core: Correct maximum permissible skin in tune_skin() r=jngrad a=RudolfWeeber

Fixes #2924 
Fix according to the discussion in #2924 




Co-authored-by: KonradBreitsprecher <[email protected]>
Co-authored-by: Konrad Breitsprecher <[email protected]>
Co-authored-by: Konrad Breitsprecher <[email protected]>
Co-authored-by: Jean-Noël Grad <[email protected]>
Co-authored-by: Kai Szuttor <[email protected]>
Co-authored-by: Kai Szuttor <[email protected]>
Co-authored-by: Rudolf Weeber <[email protected]>

doc/sphinx/inter_bonded.rst

@@ -221,25 +221,18 @@ is named ``vtol``.
 Thermalized distance bond
 ~~~~~~~~~~~~~~~~~~~~~~~~~
 
-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
-temperature and friction coefficient.
-
-The bond is configured with::
+A thermalized bond can be instantiated via
+:class:`espressomd.interactions.ThermalizedBond`::
 
     from espressomd.interactions import ThermalizedBond
     thermalized_bond = ThermalizedBond(temp_com=<float>, gamma_com=<float>,
                                        temp_distance=<float>, gamma_distance=<float>,
-                                       r_cut=<float>)
+                                       r_cut=<float>, seed=<int>)
     system.bonded_inter.add(thermalized_bond)
 
-The parameters are:
-
-    * ``temp_com``: Temperature of the Langevin thermostat for the COM of the particle pair.
-    * ``gamma_com``: Friction coefficient of the Langevin thermostat for the COM of the particle pair.
-    * ``temp_distance``: Temperature of the Langevin thermostat for the distance vector of the particle pair.
-    * ``gamma_distance``: Friction coefficient of the Langevin thermostat for the distance vector of the particle pair.
-    * ``r_cut``:  Specifies maximum distance beyond which the bond is considered broken.
+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
+temperature and friction coefficient.
 
 The bond is closely related to simulating :ref:`Particle polarizability with
 thermalized cold Drude oscillators`.

samples/drude_bmimpf6.py

@@ -279,7 +279,7 @@ def combination_rule_sigma(rule, sig1, sig2):
     thermalized_dist_bond = 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)
+        r_cut=min(lj_sigmas.values()) * 0.5, seed=123)
     harmonic_bond = HarmonicBond(k=k_drude, r_0=0.0, r_cut=1.0)
     system.bonded_inter.add(thermalized_dist_bond)
     system.bonded_inter.add(harmonic_bond)

src/core/bonded_interactions/bonded_interaction_data.hpp

@@ -3,6 +3,7 @@
 
 #include "TabulatedPotential.hpp"
 #include "particle_data.hpp"
+#include <utils/Counter.hpp>
 
 /** @file
  *  Data structures for bonded interactions.

src/core/bonded_interactions/thermalized_bond.cpp

@@ -27,9 +27,50 @@
 #include "bonded_interaction_data.hpp"
 #include "communication.hpp"
 #include "global.hpp"
+#include "random.hpp"
 
 #include <utils/constants.hpp>
 
+/** Common TB RNG Counter */
+std::unique_ptr<Utils::Counter<uint64_t>> thermalized_bond_rng_counter;
+
+/** Communication */
+void mpi_bcast_thermalized_bond_rng_counter_slave(const uint64_t counter) {
+  thermalized_bond_rng_counter =
+      std::make_unique<Utils::Counter<uint64_t>>(counter);
+}
+
+REGISTER_CALLBACK(mpi_bcast_thermalized_bond_rng_counter_slave)
+
+void mpi_bcast_thermalized_bond_rng_counter(const uint64_t counter) {
+  mpi_call(mpi_bcast_thermalized_bond_rng_counter_slave, counter);
+}
+
+void thermalized_bond_rng_counter_increment() {
+  thermalized_bond_rng_counter->increment();
+}
+
+/** Interface */
+bool thermalized_bond_is_seed_required() {
+  /* Seed is required if rng is not initialized */
+  return thermalized_bond_rng_counter == nullptr;
+}
+
+uint64_t thermalized_bond_get_rng_state() {
+  return thermalized_bond_rng_counter->value();
+}
+
+void thermalized_bond_set_rng_state(const uint64_t counter) {
+  mpi_bcast_thermalized_bond_rng_counter(counter);
+  thermalized_bond_rng_counter =
+      std::make_unique<Utils::Counter<uint64_t>>(counter);
+}
+
+/** Called each integration step */
+void thermalized_bonds_rng_counter_increment() {
+  thermalized_bond_rng_counter->increment();
+}
+
 int n_thermalized_bonds = 0;
 
 int thermalized_bond_set_params(int bond_type, double temp_com,

src/core/bonded_interactions/thermalized_bond.hpp

@@ -36,6 +36,12 @@ extern int n_thermalized_bonds;
 #include "integrate.hpp"
 #include "random.hpp"
 
+#include <Random123/philox.h>
+#include <utils/u32_to_u64.hpp>
+#include <utils/uniform.hpp>
+
+extern std::unique_ptr<Utils::Counter<uint64_t>> thermalized_bond_rng_counter;
+
 /** Set the parameters of a thermalized bond
  *
  *  @retval ES_OK on success
@@ -45,12 +51,57 @@ int thermalized_bond_set_params(int bond_type, double temp_com,
                                 double gamma_com, double temp_distance,
                                 double gamma_distance, double r_cut);
 
+/* Used in integration step */
+void thermalized_bond_rng_counter_increment();
+
+/* Interface */
+bool thermalized_bond_is_seed_required();
+uint64_t thermalized_bond_get_rng_state();
+void thermalized_bond_set_rng_state(uint64_t counter);
+
+/* Setup */
 void thermalized_bond_heat_up();
 void thermalized_bond_cool_down();
 void thermalized_bond_update_params(double pref_scale);
 void thermalized_bond_init();
 
-/** Separately thermalize the COM and distance of a particle pair.
+/** Return a random 3d vector with the philox thermostat.
+    Random numbers depend on
+    1. rng_counter (initialized by seed) which is increased on
+   integration
+    2. Salt (decorrelates different counter)
+    3. Particle ID, Partner ID
+*/
+
+inline Utils::Vector3d v_noise(int particle_id, int partner_id) {
+
+  using rng_type = r123::Philox4x64;
+  using ctr_type = rng_type::ctr_type;
+  using key_type = rng_type::key_type;
+
+  ctr_type c{{thermalized_bond_rng_counter->value(),
+              static_cast<uint64_t>(RNGSalt::THERMALIZED_BOND)}};
+
+  /** Bond is stored on particle_id, so concatenation with the
+      partner_id is unique. This will give the same RN for
+          multiple thermalized bonds on the same particle pair, which
+          should not be allowed.
+   */
+  uint64_t merged_ids;
+  auto const id1 = static_cast<uint32_t>(particle_id);
+  auto const id2 = static_cast<uint32_t>(partner_id);
+  merged_ids = Utils::u32_to_u64(id1, id2);
+  key_type k{merged_ids};
+
+  auto const noise = rng_type{}(c, k);
+
+  using Utils::uniform;
+  return Utils::Vector3d{uniform(noise[0]), uniform(noise[1]),
+                         uniform(noise[2])} -
+         Utils::Vector3d::broadcast(0.5);
+}
+
+/** Separately thermalizes the com and distance of a particle pair.
  *  @param[in]  p1        First particle.
  *  @param[in]  p2        Second particle.
  *  @param[in]  iaparams  Bonded parameters for the pair interaction.
@@ -79,23 +130,25 @@ calc_thermalized_bond_forces(Particle const *const p1, Particle const *const p2,
       mass_tot_inv * (p1->p.mass * p1->m.v + p2->p.mass * p2->m.v);
   auto const dist_vel = p2->m.v - p1->m.v;
 
+  Utils::Vector3d noise = v_noise(p1->p.identity, p2->p.identity);
+
   for (int i = 0; i < 3; i++) {
     double force_lv_com, force_lv_dist;
 
     // Langevin thermostat for center of mass
     if (iaparams->p.thermalized_bond.pref2_com > 0.0) {
-      force_lv_com = -iaparams->p.thermalized_bond.pref1_com * com_vel[i] +
-                     sqrt_mass_tot * iaparams->p.thermalized_bond.pref2_com *
-                         (d_random() - 0.5);
+      force_lv_com =
+          -iaparams->p.thermalized_bond.pref1_com * com_vel[i] +
+          sqrt_mass_tot * iaparams->p.thermalized_bond.pref2_com * noise[i];
     } else {
       force_lv_com = -iaparams->p.thermalized_bond.pref1_com * com_vel[i];
     }
 
     // Langevin thermostat for distance p1->p2
     if (iaparams->p.thermalized_bond.pref2_dist > 0.0) {
-      force_lv_dist = -iaparams->p.thermalized_bond.pref1_dist * dist_vel[i] +
-                      sqrt_mass_red * iaparams->p.thermalized_bond.pref2_dist *
-                          (d_random() - 0.5);
+      force_lv_dist =
+          -iaparams->p.thermalized_bond.pref1_dist * dist_vel[i] +
+          sqrt_mass_red * iaparams->p.thermalized_bond.pref2_dist * noise[i];
     } else {
       force_lv_dist = -iaparams->p.thermalized_bond.pref1_dist * dist_vel[i];
     }

src/core/integrate.cpp

@@ -29,6 +29,7 @@
 #include "integrate.hpp"
 #include "accumulators.hpp"
 #include "bonded_interactions/bonded_interaction_data.hpp"
+#include "bonded_interactions/thermalized_bond.hpp"
 #include "cells.hpp"
 #include "collision.hpp"
 #include "communication.hpp"
@@ -401,6 +402,8 @@ void philox_counter_increment() {
     dpd_rng_counter_increment();
 #endif
   }
+  if (n_thermalized_bonds)
+    thermalized_bond_rng_counter_increment();
 }
 
 void propagate_vel_finalize_p_inst(const ParticleRange &particles) {

src/core/random.hpp

@@ -41,7 +41,7 @@
  * noise on the particle coupling and the fluid
  * thermalization.
  */
-enum class RNGSalt { FLUID, PARTICLES, LANGEVIN, SALT_DPD };
+enum class RNGSalt { FLUID, PARTICLES, LANGEVIN, SALT_DPD, THERMALIZED_BOND };
 
 namespace Random {
 extern std::mt19937 generator;

src/core/tuning.cpp

@@ -22,6 +22,7 @@
     Implementation of tuning.hpp .
 */
 #include "communication.hpp"
+#include "domain_decomposition.hpp"
 #include "errorhandling.hpp"
 #include "global.hpp"
 #include "grid.hpp"
@@ -103,9 +104,9 @@ void tune_skin(double min_skin, double max_skin, double tol, int int_steps,
   double a = min_skin;
   double b = max_skin;
   double time_a, time_b;
-
-  auto const min_local_box_l = *boost::min_element(local_geo.length());
-  double const max_permissible_skin = 0.5 * min_local_box_l - max_cut;
+  double min_cell_size =
+      std::min(std::min(dd.cell_size[0], dd.cell_size[1]), dd.cell_size[2]);
+  double const max_permissible_skin = min_cell_size - max_cut;
 
   if (adjust_max_skin and max_skin > max_permissible_skin)
     b = max_permissible_skin;

src/python/espressomd/interactions.pxd

@@ -20,6 +20,7 @@
 # Handling of interactions
 
 from libcpp.string cimport string
+from libcpp cimport bool as cbool
 from libc cimport stdint
 
 include "myconfig.pxi"
@@ -561,6 +562,10 @@ cdef extern from "object-in-fluid/out_direction.hpp":
     int oif_out_direction_set_params(int bond_type)
 cdef extern from "bonded_interactions/thermalized_bond.hpp":
     int thermalized_bond_set_params(int bond_type, double temp_com, double gamma_com, double temp_distance, double gamma_distance, double r_cut)
+    void thermalized_bond_set_rng_state(stdint.uint64_t counter)
+    cbool thermalized_bond_is_seed_required()
+    stdint.uint64_t thermalized_bond_get_rng_state()
+
 cdef extern from "bonded_interactions/bonded_coulomb.hpp":
     int bonded_coulomb_set_params(int bond_type, double prefactor)
 cdef extern from "bonded_interactions/quartic.hpp":

src/python/espressomd/interactions.pyx

@@ -2166,13 +2166,19 @@ class ThermalizedBond(BondedInteraction):
         Friction coefficient of the Langevin thermostat for the center of mass
         of the particle pair.
     temp_distance: :obj:`float`
-        Tmperature of the Langevin thermostat for the distance vector
+        Temperature of the Langevin thermostat for the distance vector
         of the particle pair.
     gamma_distance: :obj:`float`
         Friction coefficient of the Langevin thermostat for the
         distance vector of the particle pair.
     r_cut: :obj:`float`, optional
         Maximum distance beyond which the bond is considered broken.
+    seed : :obj:`int`
+        Initial counter value (or seed) of the philox RNG.
+        Required on the first thermalized bond in the system. Must be positive.
+        If prompted, it does not return the initially set counter value
+        (the seed) but the current state of the RNG.
+
     """
 
     def __init__(self, *args, **kwargs):
@@ -2185,13 +2191,13 @@ class ThermalizedBond(BondedInteraction):
         return "THERMALIZED_DIST"
 
     def valid_keys(self):
-        return {"temp_com", "gamma_com", "temp_distance", "gamma_distance", "r_cut"}
+        return {"temp_com", "gamma_com", "temp_distance", "gamma_distance", "r_cut", "seed"}
 
     def required_keys(self):
         return {"temp_com", "gamma_com", "temp_distance", "gamma_distance"}
 
     def set_default_params(self):
-        self._params = {"r_cut": 0.}
+        self._params = {"r_cut": 0., "seed": None}
 
     def _get_params_from_es_core(self):
         return \
@@ -2204,9 +2210,21 @@ class ThermalizedBond(BondedInteraction):
              "gamma_distance":
                  bonded_ia_params[
                      self._bond_id].p.thermalized_bond.gamma_distance,
-             "r_cut": bonded_ia_params[self._bond_id].p.thermalized_bond.r_cut}
+             "r_cut": bonded_ia_params[self._bond_id].p.thermalized_bond.r_cut,
+             "seed": thermalized_bond_get_rng_state()
+             }
 
     def _set_params_in_es_core(self):
+        if self.params["seed"] is None and thermalized_bond_is_seed_required():
+            raise ValueError(
+                "A seed has to be given as keyword argument on first activation of the thermalized bond")
+        if self.params["seed"] is not None:
+            utils.check_type_or_throw_except(
+                self.params["seed"], 1, int, "seed must be a positive integer")
+            if self.params["seed"] < 0:
+                raise ValueError("seed must be a positive integer")
+            thermalized_bond_set_rng_state(self.params["seed"])
+
         thermalized_bond_set_params(
             self._bond_id, self._params["temp_com"], self._params["gamma_com"],
             self._params["temp_distance"], self._params["gamma_distance"], self._params["r_cut"])

src/python/espressomd/particle_data.pyx

@@ -151,6 +151,8 @@ cdef class ParticleHandle:
 
         def __set__(self, _pos):
             cdef double mypos[3]
+            if np.isnan(_pos).any() or np.isinf(_pos).any():
+                raise ValueError("invalid particle position")
             check_type_or_throw_except(
                 _pos, 3, float, "Position must be 3 floats")
             for i in range(3):