Rewrite LJ interaction interface

src/core/nonbonded_interactions/lj.cpp

@@ -25,31 +25,33 @@
 #include "lj.hpp"
 
 #ifdef LENNARD_JONES
-#include "interactions.hpp"
 #include "nonbonded_interaction_data.hpp"
 
-#include <utils/constants.hpp>
+#include <utils/math/int_pow.hpp>
 
-int lennard_jones_set_params(int part_type_a, int part_type_b, double eps,
-                             double sig, double cut, double shift,
-                             double offset, double min) {
-  IA_parameters *data = get_ia_param_safe(part_type_a, part_type_b);
+#include <algorithm>
+#include <stdexcept>
 
-  if (!data)
-    return ES_ERROR;
-
-  data->lj.eps = eps;
-  data->lj.sig = sig;
-  data->lj.cut = cut;
-  data->lj.shift = shift;
-  data->lj.offset = offset;
-  if (min > 0) {
-    data->lj.min = min;
+LJ_Parameters::LJ_Parameters(double eps, double sig, double cut, double offset,
+                             double min)
+    : LJ_Parameters(eps, sig, cut, offset, min, 0.) {
+  if (cut != 0.) {
+    auto const sig_cut = sig / cut;
+    shift = Utils::int_pow<6>(sig_cut) - Utils::int_pow<12>(sig_cut);
   }
-  /* broadcast interaction parameters */
-  mpi_bcast_ia_params(part_type_a, part_type_b);
+}
 
-  return ES_OK;
+LJ_Parameters::LJ_Parameters(double eps, double sig, double cut, double offset,
+                             double min, double shift)
+    : eps{eps}, sig{sig}, cut{cut}, shift{shift}, offset{offset}, min{std::max(
+                                                                      min,
+                                                                      0.)} {
+  if (eps < 0.) {
+    throw std::domain_error("Lennard-Jones parameter 'epsilon' has to be >= 0");
+  }
+  if (sig < 0.) {
+    throw std::domain_error("Lennard-Jones parameter 'sigma' has to be >= 0");
+  }
 }
 
 #endif /* ifdef LENNARD_JONES */

src/core/nonbonded_interactions/lj.hpp

@@ -37,10 +37,6 @@
 #include <utils/math/int_pow.hpp>
 #include <utils/math/sqr.hpp>
 
-int lennard_jones_set_params(int part_type_a, int part_type_b, double eps,
-                             double sig, double cut, double shift,
-                             double offset, double min);
-
 /** Calculate Lennard-Jones force factor */
 inline double lj_pair_force_factor(IA_parameters const &ia_params,
                                    double dist) {

src/core/nonbonded_interactions/nonbonded_interaction_data.hpp

@@ -48,6 +48,10 @@ struct LJ_Parameters {
   double shift = 0.0;
   double offset = 0.0;
   double min = 0.0;
+  LJ_Parameters() = default;
+  LJ_Parameters(double eps, double sig, double cut, double offset, double min);
+  LJ_Parameters(double eps, double sig, double cut, double offset, double min,
+                double shift);
 };
 
 /** WCA potential */

src/core/unit_tests/EspressoSystemStandAlone_test.cpp

@@ -37,6 +37,7 @@ namespace utf = boost::unit_test;
 #include "electrostatics/p3m.hpp"
 #include "electrostatics/registration.hpp"
 #include "energy.hpp"
+#include "event.hpp"
 #include "galilei/Galilei.hpp"
 #include "integrate.hpp"
 #include "nonbonded_interactions/lj.hpp"
@@ -126,6 +127,18 @@ static void mpi_set_tuned_p3m(double prefactor) {
 }
 #endif // P3M
 
+#ifdef LENNARD_JONES
+void mpi_set_lj_local(int key, double eps, double sig, double cut,
+                      double offset, double min, double shift) {
+  LJ_Parameters lj{eps, sig, cut, offset, min, shift};
+  ::nonbonded_ia_params[key]->lj = lj;
+  ::old_nonbonded_ia_params[key].lj = lj;
+  on_non_bonded_ia_change();
+}
+
+REGISTER_CALLBACK(mpi_set_lj_local)
+#endif // LENNARD_JONES
+
 BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory,
                         *utf::precondition(if_head_node())) {
   constexpr auto tol = 100. * std::numeric_limits<double>::epsilon();
@@ -219,9 +232,12 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory,
     auto const offset = 0.1;
     auto const min = 0.0;
     auto const r_off = dist - offset;
-    auto const cut = r_off + 1e-3; // LJ for only 2 pairs: AA BB
-    lennard_jones_set_params(type_a, type_b, eps, sig, cut, shift, offset, min);
-    lennard_jones_set_params(type_b, type_b, eps, sig, cut, shift, offset, min);
+    auto const cut = r_off + 1e-3; // LJ for only 2 pairs: AB BB
+    make_particle_type_exist(std::max(type_a, type_b));
+    auto const key_ab = get_ia_param_key(type_a, type_b);
+    auto const key_bb = get_ia_param_key(type_b, type_b);
+    mpi_call_all(mpi_set_lj_local, key_ab, eps, sig, cut, offset, min, shift);
+    mpi_call_all(mpi_set_lj_local, key_bb, eps, sig, cut, offset, min, shift);
 
     // matrix indices and reference energy value
     auto const max_type = type_b + 1;

src/core/unit_tests/Verlet_list_test.cpp

@@ -37,6 +37,7 @@ namespace bdata = boost::unit_test::data;
 #include "Particle.hpp"
 #include "ParticleFactory.hpp"
 #include "communication.hpp"
+#include "event.hpp"
 #include "integrate.hpp"
 #include "integrators/steepest_descent.hpp"
 #include "nonbonded_interactions/lj.hpp"
@@ -143,8 +144,19 @@ struct : public IntegratorHelper {
   char const *name() const override { return "VelocityVerletNpT"; }
 } velocity_verlet_npt;
 #endif // NPT
+
 } // namespace Testing
 
+void mpi_set_lj_local(int key, double eps, double sig, double cut,
+                      double offset, double min, double shift) {
+  LJ_Parameters lj{eps, sig, cut, offset, min, shift};
+  ::nonbonded_ia_params[key]->lj = lj;
+  ::old_nonbonded_ia_params[key].lj = lj;
+  on_non_bonded_ia_change();
+}
+
+REGISTER_CALLBACK(mpi_set_lj_local)
+
 inline double get_dist_from_last_verlet_update(Particle const &p) {
   return (p.pos() - p.pos_at_last_verlet_update()).norm();
 }
@@ -188,7 +200,9 @@ BOOST_DATA_TEST_CASE_F(ParticleFactory, verlet_list_update,
   auto const min = 0.0;
   auto const r_off = dist - offset;
   auto const cut = r_off + 1e-3;
-  lennard_jones_set_params(0, 1, eps, sig, cut, shift, offset, min);
+  make_particle_type_exist(1);
+  auto const key = get_ia_param_key(0, 1);
+  mpi_call_all(mpi_set_lj_local, key, eps, sig, cut, offset, min, shift);
 
   // set up velocity-Verlet integrator
   auto const time_step = 0.01;

src/python/espressomd/interactions.pxd

@@ -39,14 +39,6 @@ cdef extern from "TabulatedPotential.hpp":
         vector[double] force_tab
 
 cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp":
-    cdef struct LJ_Parameters:
-        double eps
-        double sig
-        double cut
-        double shift
-        double offset
-        double min
-
     cdef struct LJGen_Parameters:
         double eps
         double sig
@@ -147,7 +139,6 @@ cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp":
         double pref
 
     cdef struct IA_parameters:
-        LJ_Parameters lj
 
         LJcos_Parameters ljcos
 
@@ -185,13 +176,6 @@ cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp":
     cdef void ia_params_set_state(string)
     cdef void reset_ia_params()
 
-IF LENNARD_JONES:
-    cdef extern from "nonbonded_interactions/lj.hpp":
-        cdef int lennard_jones_set_params(int part_type_a, int part_type_b,
-                                          double eps, double sig, double cut,
-                                          double shift, double offset,
-                                          double min)
-
 IF LJCOS:
     cdef extern from "nonbonded_interactions/ljcos.hpp":
         cdef int ljcos_set_params(int part_type_a, int part_type_b,

src/python/espressomd/interactions.pyx

@@ -286,48 +286,21 @@ class NewNonBondedInteraction(ScriptInterfaceHelper):
 
 IF LENNARD_JONES == 1:
 
-    cdef class LennardJonesInteraction(NonBondedInteraction):
-
-        def validate_params(self):
-            """Check that parameters are valid.
-
-            Raises
-            ------
-            ValueError
-                If not true.
-            """
-            if self._params["epsilon"] < 0:
-                raise ValueError("Lennard-Jones epsilon has to be >=0")
-            if self._params["sigma"] < 0:
-                raise ValueError("Lennard-Jones sigma has to be >=0")
-            if self._params["cutoff"] < 0:
-                raise ValueError("Lennard-Jones cutoff has to be >=0")
-
-        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 {
-                "epsilon": ia_params.lj.eps,
-                "sigma": ia_params.lj.sig,
-                "cutoff": ia_params.lj.cut,
-                "shift": ia_params.lj.shift,
-                "offset": ia_params.lj.offset,
-                "min": ia_params.lj.min}
-
-        def is_active(self):
-            """Check if interaction is active.
-
-            """
-            return (self._params["epsilon"] > 0)
+    @script_interface_register
+    class LennardJonesInteraction(NewNonBondedInteraction):
+        """
+        Standard 6-12 Lennard-Jones potential.
 
-        def set_params(self, **kwargs):
-            """Set parameters for the Lennard-Jones interaction.
+        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.
 
             Parameters
             ----------
-
             epsilon : :obj:`float`
                 Magnitude of the interaction.
             sigma : :obj:`float`
@@ -343,25 +316,15 @@ IF LENNARD_JONES == 1:
             min : :obj:`float`, optional
                 Restricts the interaction to a minimal distance.
 
-            """
-            super().set_params(**kwargs)
+        """
 
-        def _set_params_in_es_core(self):
-            # Handle the case of shift="auto"
-            if self._params["shift"] == "auto":
-                self._params["shift"] = -((
-                    self._params["sigma"] / self._params["cutoff"])**12 - (
-                    self._params["sigma"] / self._params["cutoff"])**6)
-
-            if lennard_jones_set_params(
-                    self._part_types[0], self._part_types[1],
-                    self._params["epsilon"],
-                    self._params["sigma"],
-                    self._params["cutoff"],
-                    self._params["shift"],
-                    self._params["offset"],
-                    self._params["min"]):
-                raise Exception("Could not set Lennard-Jones parameters")
+        _so_name = "Interactions::InteractionLJ"
+
+        def is_active(self):
+            """Check if interaction is active.
+
+            """
+            return self.epsilon > 0.
 
         def default_params(self):
             """Python dictionary of default parameters.
@@ -1610,8 +1573,6 @@ class NonBondedInteractionHandle(ScriptInterfaceHelper):
         super().__init__(_types=[_type1, _type2], **kwargs)
 
         # Here, add one line for each nonbonded ia
-        IF LENNARD_JONES:
-            self.lennard_jones = LennardJonesInteraction(_type1, _type2)
         IF SOFT_SPHERE:
             self.soft_sphere = SoftSphereInteraction(_type1, _type2)
         IF LENNARD_JONES_GENERIC:

src/script_interface/interactions/NonBondedInteraction.hpp

@@ -154,13 +154,53 @@ class InteractionWCA : public InteractionPotentialInterface<::WCA_Parameters> {
 };
 #endif // WCA
 
+#ifdef LENNARD_JONES
+class InteractionLJ : public InteractionPotentialInterface<::LJ_Parameters> {
+protected:
+  CoreInteraction IA_parameters::*get_ptr_offset() const override {
+    return &::IA_parameters::lj;
+  }
+
+public:
+  InteractionLJ() {
+    add_parameters({
+        make_autoparameter(&CoreInteraction::eps, "epsilon"),
+        make_autoparameter(&CoreInteraction::sig, "sigma"),
+        make_autoparameter(&CoreInteraction::cut, "cutoff"),
+        make_autoparameter(&CoreInteraction::shift, "shift"),
+        make_autoparameter(&CoreInteraction::offset, "offset"),
+        make_autoparameter(&CoreInteraction::min, "min"),
+    });
+  }
+
+  void make_new_instance(VariantMap const &params) override {
+    if (auto const *shift = boost::get<std::string>(&params.at("shift"))) {
+      if (*shift != "auto") {
+        throw std::invalid_argument(
+            "Parameter 'shift' has to be 'auto' or a float");
+      }
+      m_ia_si = make_shared_from_args<CoreInteraction, double, double, double,
+                                      double, double>(
+          params, "epsilon", "sigma", "cutoff", "offset", "min");
+    } else {
+      m_ia_si = make_shared_from_args<CoreInteraction, double, double, double,
+                                      double, double, double>(
+          params, "epsilon", "sigma", "cutoff", "offset", "min", "shift");
+    }
+  }
+};
+#endif // LENNARD_JONES
+
 class NonBondedInteractionHandle
     : public AutoParameters<NonBondedInteractionHandle> {
   std::array<int, 2> m_types = {-1, -1};
   std::shared_ptr<::IA_parameters> m_interaction;
 #ifdef WCA
   std::shared_ptr<InteractionWCA> m_wca;
 #endif
+#ifdef LENNARD_JONES
+  std::shared_ptr<InteractionLJ> m_lj;
+#endif
 
   template <class T>
   auto make_autoparameter(std::shared_ptr<T> &member, const char *key) const {
@@ -181,6 +221,9 @@ class NonBondedInteractionHandle
     add_parameters({
 #ifdef WCA
         make_autoparameter(m_wca, "wca"),
+#endif
+#ifdef LENNARD_JONES
+        make_autoparameter(m_lj, "lennard_jones"),
 #endif
     });
   }
@@ -223,6 +266,10 @@ class NonBondedInteractionHandle
 #ifdef WCA
     set_member<InteractionWCA>(m_wca, "wca", "Interactions::InteractionWCA",
                                params);
+#endif
+#ifdef LENNARD_JONES
+    set_member<InteractionLJ>(m_lj, "lennard_jones",
+                              "Interactions::InteractionLJ", params);
 #endif
   }
 

src/script_interface/interactions/initialize.cpp

@@ -56,6 +56,9 @@ void initialize(Utils::Factory<ObjectHandle> *om) {
       "Interactions::NonBondedInteractions");
   om->register_new<NonBondedInteractionHandle>(
       "Interactions::NonBondedInteractionHandle");
+#ifdef LENNARD_JONES
+  om->register_new<InteractionLJ>("Interactions::InteractionLJ");
+#endif
 #ifdef WCA
   om->register_new<InteractionWCA>("Interactions::InteractionWCA");
 #endif