Thermostats refactoring (#4845)

Description of changes:
- encapsulate thermostats and remove their Cython interface
   - use one C++ `ScriptInterface` class per core class (fixes #3980)
   - introduce the `ThermostatFlags` enum, allow deactivating a single thermostat, add an `is_active` flag accessible from the Python interface (fixes #4266)
- remove 20 MPI callbacks
   - partial fix for #4614
- remove thermostat and thermalized/rigid bond global variables
   - partial fix for #2628
- remove a significant amount of historical Cython code that relied on deprecated Cython features
   - partial fix for #4542
   - ESPResSo is now compatible with Cython 3.0
   - the `nogil` deprecation warning subsists
- improve testing of the thermostats interface
- API changes:
   - thermalized bond parameter `seed` was moved to `system.thermostat.set_thermalized_bond()`
      - this change better reflects the fact there is only one global seed for all thermalized bonds
      - until now this global seed was overwritten by any newly created thermalized bond, whether it was added to the system or not
   - LB thermostat `seed` parameter now gets written to the RNG seed (silent change)
      - until now, the RNG seed was hardcoded as `0` and the `seed` value was used to set the RNG counter, thus two independent simulations with a different seed would actually get the same particle coupling random noise sequence with a time lag equal to the difference between the two seeds
      -  this is a bugfix for ensemble runs
      - fixes #4848

.github/actions/build_and_check/action.yml

@@ -6,7 +6,7 @@ runs:
     - run: |
        brew install boost boost-mpi fftw
        brew install hdf5-mpi
-       pip3 install -c requirements.txt numpy cython h5py scipy
+       pip3 install -c requirements.txt numpy "cython<3.0" h5py scipy
       shell: bash
       if: runner.os == 'macOS'
     - run: |

CMakeLists.txt

@@ -247,7 +247,7 @@ endif()
 # Python interpreter and Cython interface library
 if(ESPRESSO_BUILD_WITH_PYTHON)
   find_package(Python 3.9 REQUIRED COMPONENTS Interpreter Development NumPy)
-  find_package(Cython 0.29.21...<3.0 REQUIRED)
+  find_package(Cython 0.29.21...<3.0.8 REQUIRED)
   find_program(IPYTHON_EXECUTABLE NAMES jupyter ipython3 ipython)
 endif()
 

doc/doxygen/Doxyfile.in

@@ -319,7 +319,7 @@ TOC_INCLUDE_HEADINGS   = 5
 # globally by setting AUTOLINK_SUPPORT to NO.
 # The default value is: YES.
 
-AUTOLINK_SUPPORT       = YES
+AUTOLINK_SUPPORT       = NO
 
 # If you use STL classes (i.e. std::string, std::vector, etc.) but do not want
 # to include (a tag file for) the STL sources as input, then you should set this

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

doc/tutorials/polymers/polymers.ipynb

@@ -307,7 +307,7 @@
     "                                        kinematic_viscosity=5, tau=system.time_step,\n",
     "                                        single_precision=True)\n",
     "    system.lb = lbf\n",
-    "    system.thermostat.set_lb(LB_fluid=lbf, gamma=gamma, seed=42)"
+    "    system.thermostat.set_lb(LB_fluid=lbf, gamma=gamma, seed=0)"
    ]
   },
   {
@@ -326,6 +326,7 @@
    "outputs": [],
    "source": [
     "import logging\n",
+    "import tqdm\n",
     "import sys\n",
     "\n",
     "import numpy as np\n",
@@ -419,7 +420,7 @@
     "    rhs = np.zeros(LOOPS)\n",
     "    rfs = np.zeros(LOOPS)\n",
     "    rgs = np.zeros(LOOPS)\n",
-    "    for i in range(LOOPS):\n",
+    "    for i in tqdm.trange(LOOPS):\n",
     "        system.integrator.run(STEPS)\n",
     "        rhs[i] = system.analysis.calc_rh(\n",
     "            chain_start=0,\n",

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

requirements.txt

@@ -1,5 +1,5 @@
 # build system
-cython>=0.29.21,<3.0
+cython>=0.29.21,<=3.0.7
 setuptools>=59.6.0
 # required scientific packages
 numpy>=1.23

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/TabulatedPotential.hpp

@@ -38,7 +38,7 @@ struct TabulatedPotential {
   double minval = -1.0;
   /** Position on the x-axis of the last tabulated value. */
   double maxval = -1.0;
-  /** %Distance on the x-axis between tabulated values. */
+  /** Distance on the x-axis between tabulated values. */
   double invstepsize = 0.0;
   /** Tabulated forces. */
   std::vector<double> force_tab;

src/core/bonded_interactions/CMakeLists.txt

@@ -17,13 +17,4 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 
-target_sources(
-  espresso_core
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/angle_cosine.cpp
-          ${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)
+target_sources(espresso_core PRIVATE bonded_interaction_data.cpp)

src/core/bonded_interactions/angle_common.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 ANGLE_COMMON_H
-#define ANGLE_COMMON_H
+
+#pragma once
+
 /** \file
  *  Common code for functions calculating angle forces.
  */
@@ -95,5 +96,3 @@ angle_generic_force(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2,
   auto f_mid = -(f_left + f_right);
   return std::make_tuple(f_mid, f_left, f_right);
 }
-
-#endif /* ANGLE_COMMON_H */

src/core/bonded_interactions/angle_cosine.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 ANGLE_COSINE_H
-#define ANGLE_COSINE_H
+
+#pragma once
+
 /** \file
  *  Routines to calculate the angle energy or/and and force
  *  for a particle triple using the potential described in
@@ -50,7 +51,12 @@ struct AngleCosineBond {
 
   static constexpr int num = 2;
 
-  AngleCosineBond(double bend, double phi0);
+  AngleCosineBond(double bend, double phi0) {
+    this->bend = bend;
+    this->phi0 = phi0;
+    this->cos_phi0 = cos(phi0);
+    this->sin_phi0 = sin(phi0);
+  }
 
   std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
   forces(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2) const;
@@ -98,5 +104,3 @@ inline double AngleCosineBond::energy(Utils::Vector3d const &vec1,
   // trig identity: cos(phi - phi0) = cos(phi)cos(phi0) + sin(phi)sin(phi0)
   return bend * (1 - (cos_phi * cos_phi0 + sin_phi * sin_phi0));
 }
-
-#endif /* ANGLE_COSINE_H */

src/core/bonded_interactions/angle_cossquare.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 ANGLE_COSSQUARE_H
-#define ANGLE_COSSQUARE_H
+
+#pragma once
+
 /** \file
  *  Routines to calculate the angle energy or/and and force
  *  for a particle triple using the potential described in
@@ -31,6 +32,7 @@
 #include <utils/Vector.hpp>
 #include <utils/math/sqr.hpp>
 
+#include <cmath>
 #include <tuple>
 
 /** Parameters for three-body angular potential (cossquare). */
@@ -46,7 +48,11 @@ struct AngleCossquareBond {
 
   static constexpr int num = 2;
 
-  AngleCossquareBond(double bend, double phi0);
+  AngleCossquareBond(double bend, double phi0) {
+    this->bend = bend;
+    this->phi0 = phi0;
+    this->cos_phi0 = cos(phi0);
+  }
 
   std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
   forces(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2) const;
@@ -87,5 +93,3 @@ inline double AngleCossquareBond::energy(Utils::Vector3d const &vec1,
   auto const cos_phi = calc_cosine(vec1, vec2, true);
   return 0.5 * bend * Utils::sqr(cos_phi - cos_phi0);
 }
-
-#endif /* ANGLE_COSSQUARE_H */

src/core/bonded_interactions/angle_harmonic.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 ANGLE_HARMONIC_H
-#define ANGLE_HARMONIC_H
+
+#pragma once
+
 /** \file
  *  Routines to calculate the angle energy or/and and force
  *  for a particle triple using the potential described in
@@ -91,5 +92,3 @@ inline double AngleHarmonicBond::energy(Utils::Vector3d const &vec1,
   auto const phi = acos(cos_phi);
   return 0.5 * bend * Utils::sqr(phi - phi0);
 }
-
-#endif /* ANGLE_HARMONIC_H */