Rewrite Hertzian interaction interface, add docs

src/core/nonbonded_interactions/hertzian.cpp

@@ -25,25 +25,13 @@
 #include "hertzian.hpp"
 
 #ifdef HERTZIAN
-#include "interactions.hpp"
 #include "nonbonded_interaction_data.hpp"
 
-#include <utils/constants.hpp>
-
-int hertzian_set_params(int part_type_a, int part_type_b, double eps,
-                        double sig) {
-  IA_parameters *data = get_ia_param_safe(part_type_a, part_type_b);
-
-  if (!data)
-    return ES_ERROR;
-
-  data->hertzian.eps = eps;
-  data->hertzian.sig = sig;
-
-  /* broadcast interaction parameters */
-  mpi_bcast_ia_params(part_type_a, part_type_b);
-
-  return ES_OK;
+Hertzian_Parameters::Hertzian_Parameters(double eps, double sig)
+    : eps{eps}, sig{sig} {
+  if (eps < 0.) {
+    throw std::domain_error("Hertzian parameter 'epsilon' has to be >= 0");
+  }
 }
 
-#endif
+#endif // HERTZIAN

src/core/nonbonded_interactions/hertzian.hpp

@@ -31,15 +31,10 @@
 
 #include "nonbonded_interaction_data.hpp"
 
-#include <utils/Vector.hpp>
-
 #include <cmath>
 
 #ifdef HERTZIAN
 
-int hertzian_set_params(int part_type_a, int part_type_b, double eps,
-                        double sig);
-
 /** Calculate Hertzian force factor */
 inline double hertzian_pair_force_factor(IA_parameters const &ia_params,
                                          double dist) {

src/core/nonbonded_interactions/nonbonded_interaction_data.hpp

@@ -94,6 +94,8 @@ struct SmoothStep_Parameters {
 struct Hertzian_Parameters {
   double eps = 0.0;
   double sig = INACTIVE_CUTOFF;
+  Hertzian_Parameters() = default;
+  Hertzian_Parameters(double eps, double sig);
 };
 
 /** Gaussian potential */

src/python/espressomd/interactions.pxd

@@ -60,10 +60,6 @@ cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp":
         int n
         double k0
 
-    cdef struct Hertzian_Parameters:
-        double eps
-        double sig
-
     cdef struct BMHTF_Parameters:
         double A
         double B
@@ -139,8 +135,6 @@ cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp":
 
         Buckingham_Parameters buckingham
 
-        Hertzian_Parameters hertzian
-
         DPDParameters dpd_radial
         DPDParameters dpd_trans
 
@@ -207,11 +201,6 @@ IF SOFT_SPHERE:
         int soft_sphere_set_params(int part_type_a, int part_type_b,
                                    double a, double n, double cut, double offset)
 
-IF HERTZIAN:
-    cdef extern from "nonbonded_interactions/hertzian.hpp":
-        int hertzian_set_params(int part_type_a, int part_type_b,
-                                double eps, double sig)
-
 IF DPD:
     cdef extern from "dpd.hpp":
         int dpd_set_params(int part_type_a, int part_type_b,

src/python/espressomd/interactions.pyx

@@ -598,6 +598,9 @@ IF LJCOS == 1:
             return self.epsilon > 0.
 
         def default_params(self):
+            """Python dictionary of default parameters.
+
+            """
             return {"offset": 0.}
 
         def type_name(self):
@@ -647,9 +650,15 @@ IF LJCOS2 == 1:
         _so_name = "Interactions::InteractionLJcos2"
 
         def is_active(self):
+            """Check if interactions is active.
+
+            """
             return self.epsilon > 0.
 
         def default_params(self):
+            """Python dictionary of default parameters.
+
+            """
             return {"offset": 0.}
 
         def type_name(self):
@@ -1355,42 +1364,35 @@ IF SOFT_SPHERE == 1:
             """
             return {"a", "n", "cutoff"}
 
-
 IF HERTZIAN == 1:
 
-    cdef class HertzianInteraction(NonBondedInteraction):
+    @script_interface_register
+    class HertzianInteraction(NewNonBondedInteraction):
+        """Hertzian interaction.
 
-        def validate_params(self):
-            """Check that parameters are valid.
+        Methods
+        -------
+        set_params()
+            Set or update parameters for the interaction.
+            Parameters marked as required become optional once the
+            interaction has been activated for the first time;
+            subsequent calls to this method update the existing values.
 
-            """
-            if self._params["eps"] < 0:
-                raise ValueError("Hertzian eps a has to be >=0")
-            if self._params["sig"] < 0:
-                raise ValueError("Hertzian sig has to be >=0")
+            Parameters
+            ----------
+            epsilon : :obj:`float`
+                Magnitude of the interaction.
+            sigma : :obj:`float`
+                Interaction length scale.
+        """
 
-        def _get_params_from_es_core(self):
-            cdef IA_parameters * ia_params
-            ia_params = get_ia_param_safe(
-                self._part_types[0],
-                self._part_types[1])
-            return {
-                "eps": ia_params.hertzian.eps,
-                "sig": ia_params.hertzian.sig
-            }
+        _so_name = "Interactions::InteractionHertzian"
 
         def is_active(self):
             """Check if interaction is active.
 
             """
-            return (self._params["eps"] > 0)
-
-        def _set_params_in_es_core(self):
-            if hertzian_set_params(self._part_types[0],
-                                   self._part_types[1],
-                                   self._params["eps"],
-                                   self._params["sig"]):
-                raise Exception("Could not set Hertzian parameters")
+            return self.epsilon > 0.
 
         def default_params(self):
             """Python dictionary of default parameters.
@@ -1404,32 +1406,17 @@ IF HERTZIAN == 1:
             """
             return "Hertzian"
 
-        def set_params(self, **kwargs):
-            """
-            Set parameters for the Hertzian interaction.
-
-            Parameters
-            ----------
-            eps : :obj:`float`
-                Magnitude of the interaction.
-            sig : :obj:`float`
-                Parameter sigma. Determines the length over which the potential
-                decays.
-
-            """
-            super().set_params(**kwargs)
-
         def valid_keys(self):
             """All parameters that can be set.
 
             """
-            return {"eps", "sig"}
+            return {"epsilon", "sigma"}
 
         def required_keys(self):
             """Parameters that have to be set.
 
             """
-            return {"eps", "sig"}
+            return {"epsilon", "sigma"}
 
 IF GAUSSIAN == 1:
 
@@ -1529,8 +1516,6 @@ class NonBondedInteractionHandle(ScriptInterfaceHelper):
             self.morse = MorseInteraction(_type1, _type2)
         IF BUCKINGHAM:
             self.buckingham = BuckinghamInteraction(_type1, _type2)
-        IF HERTZIAN:
-            self.hertzian = HertzianInteraction(_type1, _type2)
         IF TABULATED:
             self.tabulated = TabulatedNonBonded(_type1, _type2)
         IF GAY_BERNE:

src/script_interface/interactions/NonBondedInteraction.hpp

@@ -280,6 +280,28 @@ class InteractionGaussian
 };
 #endif // GAUSSIAN
 
+#ifdef HERTZIAN
+class InteractionHertzian
+    : public InteractionPotentialInterface<::Hertzian_Parameters> {
+protected:
+  CoreInteraction IA_parameters::*get_ptr_offset() const override {
+    return &::IA_parameters::hertzian;
+  }
+
+public:
+  InteractionHertzian() {
+    add_parameters({
+        make_autoparameter(&CoreInteraction::eps, "epsilon"),
+        make_autoparameter(&CoreInteraction::sig, "sigma"),
+    });
+  }
+  void make_new_instance(VariantMap const &params) override {
+    m_ia_si = make_shared_from_args<CoreInteraction, double, double>(
+        params, "epsilon", "sigma");
+  }
+};
+#endif // HERTZIAN
+
 class NonBondedInteractionHandle
     : public AutoParameters<NonBondedInteractionHandle> {
   std::array<int, 2> m_types = {-1, -1};
@@ -299,6 +321,9 @@ class NonBondedInteractionHandle
 #ifdef GAUSSIAN
   std::shared_ptr<InteractionGaussian> m_gaussian;
 #endif
+#ifdef HERTZIAN
+  std::shared_ptr<InteractionHertzian> m_hertzian;
+#endif
 
   template <class T>
   auto make_autoparameter(std::shared_ptr<T> &member, const char *key) const {
@@ -331,6 +356,9 @@ class NonBondedInteractionHandle
 #endif
 #ifdef GAUSSIAN
         make_autoparameter(m_gaussian, "gaussian"),
+#endif
+#ifdef HERTZIAN
+        make_autoparameter(m_hertzian, "hertzian"),
 #endif
     });
   }
@@ -389,6 +417,10 @@ class NonBondedInteractionHandle
 #ifdef GAUSSIAN
     set_member<InteractionGaussian>(
         m_gaussian, "gaussian", "Interactions::InteractionGaussian", params);
+#endif
+#ifdef HERTZIAN
+    set_member<InteractionHertzian>(
+        m_hertzian, "hertzian", "Interactions::InteractionHertzian", params);
 #endif
   }
 

src/script_interface/interactions/initialize.cpp

@@ -59,6 +59,9 @@ void initialize(Utils::Factory<ObjectHandle> *om) {
 #ifdef GAUSSIAN
   om->register_new<InteractionGaussian>("Interactions::InteractionGaussian");
 #endif
+#ifdef HERTZIAN
+  om->register_new<InteractionHertzian>("Interactions::InteractionHertzian");
+#endif
 #ifdef LENNARD_JONES
   om->register_new<InteractionLJ>("Interactions::InteractionLJ");
 #endif

testsuite/python/interactions_non-bonded.py

@@ -196,17 +196,17 @@ def soft_sphere_force(r, a, n, cutoff, offset=0):
 # Hertzian
 
 
-def hertzian_potential(r, eps, sig):
+def hertzian_potential(r, epsilon, sigma):
     V = 0.
-    if r < sig:
-        V = eps * np.power(1 - r / sig, 5. / 2.)
+    if r < sigma:
+        V = epsilon * np.power(1 - r / sigma, 5. / 2.)
     return V
 
 
-def hertzian_force(r, eps, sig):
+def hertzian_force(r, epsilon, sigma):
     f = 0.
-    if r < sig:
-        f = 5. / 2. * eps / sig * np.power(1 - r / sig, 3. / 2.)
+    if r < sigma:
+        f = 5. / 2. * epsilon / sigma * np.power(1 - r / sigma, 3. / 2.)
     return f
 
 # Gaussian
@@ -467,8 +467,8 @@ def test_soft_sphere(self):
     def test_hertzian(self):
 
         self.run_test("hertzian",
-                      {"eps": 6.92,
-                       "sig": 2.432},
+                      {"epsilon": 6.92,
+                       "sigma": 2.432},
                       force_kernel=hertzian_force,
                       energy_kernel=hertzian_potential,
                       n_steps=244)