Rewrite non-bonded interactions as script interface objects (#4558)

Fixes #4555 and fixes #4569

Description of changes:
- rewrite non-bonded interactions as script interface objects
   - convert Cython file `interaction.pyx` to Python file `interaction.py`
   - 1600 lines of Cython code were removed
- bugfix: the `set_params()` method now raises an error if required arguments are not provided
   - scripts may silently break when a potential is set and later updated in a way that is affected by the bug described in #4569
- new feature: possibility to deactivate/reset non-bonded interactions
- reduce code duplication in the script interface of bonded interactions

.gitlab-ci.yml

@@ -438,6 +438,8 @@ empty:
      with_static_analysis: 'true'
      with_scafacos: 'false'
      with_stokesian_dynamics: 'false'
+     with_coverage: 'false'
+     with_coverage_python: 'true'
   script:
     - bash maintainer/CI/build_cmake.sh
   tags:

doc/sphinx/inter_non-bonded.rst

@@ -18,23 +18,41 @@ explicitly set. A bonded interaction between a set of particles has to
 be specified explicitly by the command, while the command is used to
 define the interaction parameters.
 
-.. _Isotropic non-bonded interactions:
+Non-bonded interaction are configured via the
+:class:`espressomd.interactions.NonBondedInteraction` class,
+which is a member of :class:`espressomd.system.System`::
 
-Isotropic non-bonded interactions
----------------------------------
+    system.non_bonded_inter[type1, type2]
 
-Non-bonded interaction are configured via the :class:`espressomd.interactions.NonBondedInteraction` class, which is a member of :class:`espressomd.system.System`::
+This command defines an interaction between all particles of type ``type1``
+and ``type2``. All available interaction potentials and their parameters are
+listed below. For example, the following adds a WCA potential between
+particles of type 0::
 
-    system.non_bonded_inter[type1, type2]
+    system.non_bonded_inter[0, 0].wca.set_params(epsilon=1., sigma=2.)
+
+Each type pair can have multiple potentials active at the same time.
+To deactivate a specific potential for a given type pair, do::
+
+    system.non_bonded_inter[0, 0].wca.deactivate()
+
+To deactivate all potentials for a given type pair, do::
 
-This command defines an interaction between all particles of type ``type1`` and
-``type2``. Possible interaction types and their parameters are
-listed below.
+    system.non_bonded_inter[0, 0].reset()
+
+To deactivate all potentials between all type pairs, do::
+
+    system.non_bonded_inter.reset()
 
 For many non-bonded interactions, it is possible to artificially cap the
 forces, which often allows to equilibrate the system much faster. See
 the subsection :ref:`Capping the force during warmup` for more details.
 
+.. _Isotropic non-bonded interactions:
+
+Isotropic non-bonded interactions
+---------------------------------
+
 .. _Tabulated interaction:
 
 Tabulated interaction

src/core/CMakeLists.txt

@@ -55,7 +55,8 @@ add_library(
   virtual_sites.cpp
   exclusions.cpp
   PartCfg.cpp
-  EspressoSystemStandAlone.cpp)
+  EspressoSystemStandAlone.cpp
+  TabulatedPotential.cpp)
 add_library(espresso::core ALIAS espresso_core)
 
 if(CUDA)

src/core/TabulatedPotential.cpp

@@ -0,0 +1,55 @@
+/*
+ * Copyright (C) 2010-2022 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "TabulatedPotential.hpp"
+
+#include <stdexcept>
+#include <vector>
+
+TabulatedPotential::TabulatedPotential(double minval, double maxval,
+                                       std::vector<double> const &force,
+                                       std::vector<double> const &energy)
+    : minval{minval}, maxval{maxval} {
+
+  if (minval > maxval) {
+    throw std::domain_error("TabulatedPotential parameter 'max' must be "
+                            "larger than or equal to parameter 'min'");
+  }
+  if (minval != -1.) {
+    if (minval == maxval and force.size() != 1) {
+      throw std::domain_error(
+          "TabulatedPotential parameter 'force' must contain 1 element");
+    }
+    if (force.empty()) {
+      throw std::domain_error("TabulatedPotential parameter 'force' must "
+                              "contain at least 1 element");
+    }
+    if (force.size() != energy.size()) {
+      throw std::invalid_argument("TabulatedPotential parameter 'force' must "
+                                  "have the same size as parameter 'energy'");
+    }
+    invstepsize = static_cast<double>(force.size() - 1) / (maxval - minval);
+  } else {
+    invstepsize = 0.;
+  }
+  force_tab = force;
+  energy_tab = energy;
+}

src/core/TabulatedPotential.hpp

@@ -25,7 +25,6 @@
 #include <boost/serialization/access.hpp>
 #include <boost/serialization/vector.hpp>
 
-#include <cassert>
 #include <vector>
 
 /** Evaluate forces and energies using a custom potential profile.
@@ -46,6 +45,11 @@ struct TabulatedPotential {
   /** Tabulated energies. */
   std::vector<double> energy_tab;
 
+  TabulatedPotential() = default;
+  TabulatedPotential(double minval, double maxval,
+                     std::vector<double> const &force,
+                     std::vector<double> const &energy);
+
   /** Evaluate the force at position @p x.
    *  @param x  Bond length/angle
    *  @return Interpolated force.
@@ -67,17 +71,6 @@ struct TabulatedPotential {
   }
 
   double cutoff() const { return maxval; }
-
-private:
-  friend boost::serialization::access;
-  template <typename Archive>
-  void serialize(Archive &ar, long int /* version */) {
-    ar &minval;
-    ar &maxval;
-    ar &invstepsize;
-    ar &force_tab;
-    ar &energy_tab;
-  }
 };
 
 #endif

src/core/bonded_interactions/bonded_interaction_data.hpp

@@ -133,6 +133,9 @@ class BondedInteractionsMap {
   bool empty() const { return m_params.empty(); }
   auto size() const { return m_params.size(); }
   auto get_next_key() const { return next_key; }
+  auto get_zero_based_type(int bond_id) const {
+    return contains(bond_id) ? at(bond_id)->which() : 0;
+  }
 
 private:
   container_type m_params = {};

src/core/bonded_interactions/bonded_interactions.dox

@@ -243,37 +243,27 @@
  *
  *  @subsection bondedIA_new_interface Adding the interaction in the Python interface
  *
- *  Please note that the following is Cython code (www.cython.org), rather than
- *  pure Python.
- *
- *  * In file <tt>src/python/espressomd/interactions.pxd</tt>:
- *    - Add the bonded interaction to \c enum_bonded_interaction.
+ *  * In file <tt>src/python/espressomd/interactions.py</tt>:
+ *    - Add the bonded interaction to \c BONDED_IA.
  *      The order of the enum values must match the order of types in
  *      @ref Bonded_IA_Parameters exactly:
  *      @code{.py}
- *      cdef enum enum_bonded_interaction:
- *          BONDED_IA_NONE = 0,
- *          BONDED_IA_FENE,
- *          [...]
+ *      class BONDED_IA(enum.IntEnum):
+ *          NONE = 0
+ *          FENE = enum.auto()
+ *          HARMONIC = enum.auto()
+ *          # ...
  *      @endcode
- *  * In file <tt>src/python/espressomd/interactions.pyx</tt>:
- *    - Implement the Cython class for the bonded interaction, using the one
+ *    - Implement the Python class for the bonded interaction, using the one
  *      for the FENE bond as template. Please use pep8 naming convention.
- *    - Set the class' member
+ *    - Set the class' members
  *      @code{.py}
  *      _so_name = "Interactions::YourNewBond"
+ *      _type_number = BONDED_IA.YOURNEWBOND
  *      @endcode
  *      where you put the name of your bond instead of \c YourNewBond.
  *      This connects the %ScriptInterface class with the Python class.
- *    - Create a new entry in the dictionary \c bonded_interaction_classes to
- *      register the new class you have written:
- *      @code{.py}
- *      bonded_interaction_classes = {
- *          int(BONDED_IA_FENE): FeneBond,
- *          int(BONDED_IA_HARMONIC): HarmonicBond,
- *          [...]
- *      }
- *      @endcode
+ *      The type number is the new enum value from \c BONDED_IA.
  *  * In file <tt>testsuite/python/interactions_bonded_interface.py</tt>:
  *    - Add a test case to verify that parameters set and gotten from the
  *      interaction are consistent.

src/core/bonded_interactions/bonded_tab.cpp

@@ -30,20 +30,7 @@
 TabulatedBond::TabulatedBond(double min, double max,
                              std::vector<double> const &energy,
                              std::vector<double> const &force) {
-  assert(max >= min);
-  assert((max == min) || force.size() > 1);
-  assert(force.size() == energy.size());
-
-  auto tab_pot = this->pot = std::make_shared<TabulatedPotential>();
-
-  /* set table limits */
-  tab_pot->minval = min;
-  tab_pot->maxval = max;
-
-  tab_pot->invstepsize = static_cast<double>(force.size() - 1) / (max - min);
-
-  tab_pot->force_tab = force;
-  tab_pot->energy_tab = energy;
+  pot = std::make_shared<TabulatedPotential>(min, max, force, energy);
 }
 
 TabulatedDistanceBond::TabulatedDistanceBond(double min, double max,

src/core/dpd.cpp

@@ -64,28 +64,15 @@ Vector3d dpd_noise(int pid1, int pid2) {
       (pid1 < pid2) ? pid1 : pid2);
 }
 
-int dpd_set_params(int part_type_a, int part_type_b, double gamma, double k,
-                   double r_c, int wf, double tgamma, double tr_c, int twf) {
-  auto &ia_params = *get_ia_param_safe(part_type_a, part_type_b);
-
-  ia_params.dpd_radial = DPDParameters{gamma, k, r_c, wf, -1.};
-  ia_params.dpd_trans = DPDParameters{tgamma, k, tr_c, twf, -1.};
-
-  /* broadcast interaction parameters */
-  mpi_bcast_ia_params(part_type_a, part_type_b);
-
-  return ES_OK;
-}
-
 void dpd_init(double kT, double time_step) {
   for (int type_a = 0; type_a < max_seen_particle_type; type_a++) {
     for (int type_b = 0; type_b < max_seen_particle_type; type_b++) {
       IA_parameters &ia_params = get_ia_param(type_a, type_b);
 
-      ia_params.dpd_radial.pref =
-          sqrt(24.0 * kT * ia_params.dpd_radial.gamma / time_step);
-      ia_params.dpd_trans.pref =
-          sqrt(24.0 * kT * ia_params.dpd_trans.gamma / time_step);
+      ia_params.dpd.radial.pref =
+          sqrt(24.0 * kT * ia_params.dpd.radial.gamma / time_step);
+      ia_params.dpd.trans.pref =
+          sqrt(24.0 * kT * ia_params.dpd.trans.gamma / time_step);
     }
   }
 }
@@ -116,19 +103,19 @@ Utils::Vector3d dpd_pair_force(Particle const &p1, Particle const &p2,
                                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) {
+  if (ia_params.dpd.radial.cutoff <= 0.0 && ia_params.dpd.trans.cutoff <= 0.0) {
     return {};
   }
 
   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)
+      (ia_params.dpd.radial.pref > 0.0 || ia_params.dpd.trans.pref > 0.0)
           ? dpd_noise(p1.id(), p2.id())
           : Vector3d{};
 
-  auto const f_r = dpd_pair_force(ia_params.dpd_radial, v21, dist, noise_vec);
-  auto const f_t = dpd_pair_force(ia_params.dpd_trans, v21, dist, noise_vec);
+  auto const f_r = dpd_pair_force(ia_params.dpd.radial, v21, dist, noise_vec);
+  auto const f_t = dpd_pair_force(ia_params.dpd.trans, v21, dist, noise_vec);
 
   /* Projection operator to radial direction */
   auto const P = tensor_product(d / dist2, d);
@@ -150,8 +137,8 @@ static auto dpd_viscous_stress_local() {
         auto const &ia_params = get_ia_param(p1.type(), p2.type());
         auto const dist = std::sqrt(d.dist2);
 
-        auto const f_r = dpd_pair_force(ia_params.dpd_radial, v21, dist, {});
-        auto const f_t = dpd_pair_force(ia_params.dpd_trans, v21, dist, {});
+        auto const f_r = dpd_pair_force(ia_params.dpd.radial, v21, dist, {});
+        auto const f_t = dpd_pair_force(ia_params.dpd.trans, v21, dist, {});
 
         /* Projection operator to radial direction */
         auto const P = tensor_product(d.vec21 / d.dist2, d.vec21);
@@ -188,4 +175,5 @@ Utils::Vector9d dpd_stress() {
 
   return Utils::flatten(stress) / volume;
 }
-#endif
+
+#endif // DPD

src/core/dpd.hpp

@@ -36,8 +36,6 @@
 
 struct IA_parameters;
 
-int dpd_set_params(int part_type_a, int part_type_b, double gamma, double k,
-                   double r_c, int wf, double tgamma, double tr_c, int twf);
 void dpd_init(double kT, double time_step);
 
 Utils::Vector3d dpd_pair_force(Particle const &p1, Particle const &p2,

src/core/event.cpp

@@ -251,6 +251,11 @@ void on_dipoles_change() {
   on_short_range_ia_change();
 }
 
+void on_non_bonded_ia_change() {
+  maximal_cutoff_nonbonded();
+  on_short_range_ia_change();
+}
+
 void on_short_range_ia_change() {
   cells_re_init(cell_structure.decomposition_type());
   recalc_forces = true;

src/core/event.hpp

@@ -79,6 +79,9 @@ void on_dipoles_change();
 /** called every time short ranged interaction parameters are changed. */
 void on_short_range_ia_change();
 
+/** called every time a non-bonded interaction parameters are changed. */
+void on_non_bonded_ia_change();
+
 /** called every time a constraint is changed. */
 void on_constraint_change();
 

src/core/interactions.cpp

@@ -18,16 +18,15 @@
  * You should have received a copy of the GNU General Public License
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
-#include "communication.hpp"
 
 #include "bonded_interactions/bonded_interaction_data.hpp"
 #include "collision.hpp"
+#include "communication.hpp"
 #include "electrostatics/coulomb.hpp"
 #include "errorhandling.hpp"
 #include "event.hpp"
 #include "magnetostatics/dipoles.hpp"
-
-#include "serialization/IA_parameters.hpp"
+#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
 
 #include <boost/mpi/collectives/broadcast.hpp>
 
@@ -77,14 +76,3 @@ bool long_range_interactions_sanity_checks() {
   }
   return false;
 }
-
-void mpi_bcast_ia_params_local(int i, int j) {
-  boost::mpi::broadcast(comm_cart, get_ia_param(i, j), 0);
-  on_short_range_ia_change();
-}
-
-REGISTER_CALLBACK(mpi_bcast_ia_params_local)
-
-void mpi_bcast_ia_params(int i, int j) {
-  mpi_call_all(mpi_bcast_ia_params_local, i, j);
-}