From 8576bf891fe4259f3081501b1268b8777b4cb7bb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 5 Sep 2022 22:31:04 +0200 Subject: [PATCH] Rewrite generic Lennard-Jones interaction interface --- src/core/nonbonded_interactions/ljgen.cpp | 55 +++---- src/core/nonbonded_interactions/ljgen.hpp | 23 +-- .../nonbonded_interaction_data.cpp | 3 +- .../nonbonded_interaction_data.hpp | 20 ++- src/python/espressomd/interactions.pxd | 29 ---- src/python/espressomd/interactions.pyx | 140 +++++------------- .../interactions/NonBondedInteraction.hpp | 65 ++++++++ .../interactions/initialize.cpp | 3 + testsuite/python/interactions_non-bonded.py | 29 ++-- .../interactions_non-bonded_interface.py | 19 +++ 10 files changed, 191 insertions(+), 195 deletions(-) diff --git a/src/core/nonbonded_interactions/ljgen.cpp b/src/core/nonbonded_interactions/ljgen.cpp index 862d0dfef81..51feb1bc448 100644 --- a/src/core/nonbonded_interactions/ljgen.cpp +++ b/src/core/nonbonded_interactions/ljgen.cpp @@ -25,43 +25,36 @@ #include "ljgen.hpp" #ifdef LENNARD_JONES_GENERIC -#include "interactions.hpp" #include "nonbonded_interaction_data.hpp" -#include - -int ljgen_set_params(int part_type_a, int part_type_b, double eps, double sig, - double cut, double shift, double offset, double a1, - double a2, double b1, double b2 +#include +LJGen_Parameters::LJGen_Parameters(double epsilon, double sigma, double cutoff, + double shift, double offset, #ifdef LJGEN_SOFTCORE - , - double lambda, double softrad + double lam, double delta, #endif -) { - IA_parameters *data = get_ia_param_safe(part_type_a, part_type_b); - - if (!data) - return ES_ERROR; - - data->ljgen.eps = eps; - data->ljgen.sig = sig; - data->ljgen.cut = cut; - data->ljgen.shift = shift; - data->ljgen.offset = offset; - data->ljgen.a1 = a1; - data->ljgen.a2 = a2; - data->ljgen.b1 = b1; - data->ljgen.b2 = b2; + double e1, double e2, double b1, double b2) + : eps{epsilon}, sig{sigma}, cut{cutoff}, shift{shift}, offset{offset}, #ifdef LJGEN_SOFTCORE - data->ljgen.lambda1 = lambda; - data->ljgen.softrad = softrad; + lambda{lam}, softrad{delta}, #endif - - /* broadcast interaction parameters */ - mpi_bcast_ia_params(part_type_a, part_type_b); - - return ES_OK; + a1{e1}, a2{e2}, b1{b1}, b2{b2} { + if (epsilon < 0.) { + throw std::domain_error("Generic LJ parameter 'epsilon' has to be >= 0"); + } + if (sigma < 0.) { + throw std::domain_error("Generic LJ parameter 'sigma' has to be >= 0"); + } +#ifdef LJGEN_SOFTCORE + if (delta < 0.) { + throw std::domain_error("Generic LJ parameter 'delta' has to be >= 0"); + } + if (lam < 0. or lam > 1.) { + throw std::domain_error( + "Generic LJ parameter 'lam' has to be in the range [0, 1]"); + } +#endif // LJGEN_SOFTCORE } -#endif +#endif // LENNARD_JONES_GENERIC diff --git a/src/core/nonbonded_interactions/ljgen.hpp b/src/core/nonbonded_interactions/ljgen.hpp index b0a734ff514..c753ba1db20 100644 --- a/src/core/nonbonded_interactions/ljgen.hpp +++ b/src/core/nonbonded_interactions/ljgen.hpp @@ -47,24 +47,15 @@ #include -int ljgen_set_params(int part_type_a, int part_type_b, double eps, double sig, - double cut, double shift, double offset, double a1, - double a2, double b1, double b2 -#ifdef LJGEN_SOFTCORE - , - double lambda, double softrad -#endif -); - /** Calculate Lennard-Jones force factor */ inline double ljgen_pair_force_factor(IA_parameters const &ia_params, double dist) { - if (dist < (ia_params.ljgen.cut + ia_params.ljgen.offset)) { + if (dist < ia_params.ljgen.max_cutoff()) { auto r_off = dist - ia_params.ljgen.offset; #ifdef LJGEN_SOFTCORE r_off *= r_off; - r_off += Utils::sqr(ia_params.ljgen.sig) * (1.0 - ia_params.ljgen.lambda1) * + r_off += Utils::sqr(ia_params.ljgen.sig) * (1.0 - ia_params.ljgen.lambda) * ia_params.ljgen.softrad; r_off = sqrt(r_off); #else @@ -73,7 +64,7 @@ inline double ljgen_pair_force_factor(IA_parameters const &ia_params, auto const frac = ia_params.ljgen.sig / r_off; auto const fac = ia_params.ljgen.eps #ifdef LJGEN_SOFTCORE - * ia_params.ljgen.lambda1 * + * ia_params.ljgen.lambda * (dist - ia_params.ljgen.offset) / r_off #endif * (ia_params.ljgen.b1 * ia_params.ljgen.a1 * @@ -88,11 +79,11 @@ inline double ljgen_pair_force_factor(IA_parameters const &ia_params, /** Calculate Lennard-Jones energy */ inline double ljgen_pair_energy(IA_parameters const &ia_params, double dist) { - if (dist < (ia_params.ljgen.cut + ia_params.ljgen.offset)) { + if (dist < ia_params.ljgen.max_cutoff()) { auto r_off = dist - ia_params.ljgen.offset; #ifdef LJGEN_SOFTCORE r_off *= r_off; - r_off += pow(ia_params.ljgen.sig, 2) * (1.0 - ia_params.ljgen.lambda1) * + r_off += pow(ia_params.ljgen.sig, 2) * (1.0 - ia_params.ljgen.lambda) * ia_params.ljgen.softrad; r_off = sqrt(r_off); #else @@ -101,7 +92,7 @@ inline double ljgen_pair_energy(IA_parameters const &ia_params, double dist) { auto const frac = ia_params.ljgen.sig / r_off; auto const fac = ia_params.ljgen.eps #ifdef LJGEN_SOFTCORE - * ia_params.ljgen.lambda1 + * ia_params.ljgen.lambda #endif * (ia_params.ljgen.b1 * pow(frac, ia_params.ljgen.a1) - ia_params.ljgen.b2 * pow(frac, ia_params.ljgen.a2) + @@ -112,4 +103,4 @@ inline double ljgen_pair_energy(IA_parameters const &ia_params, double dist) { } #endif // LENNARD_JONES_GENERIC -#endif // CORE_NB_IA_LJGEN_HPP +#endif diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp index 68ca1341910..c2801e512ba 100644 --- a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp +++ b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp @@ -150,8 +150,7 @@ static double recalc_maximal_cutoff(const IA_parameters &data) { #endif #ifdef LENNARD_JONES_GENERIC - max_cut_current = - std::max(max_cut_current, (data.ljgen.cut + data.ljgen.offset)); + max_cut_current = std::max(max_cut_current, data.ljgen.max_cutoff()); #endif #ifdef SMOOTH_STEP diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp index c48b4823a82..922638e3b43 100644 --- a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp +++ b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp @@ -31,6 +31,7 @@ #include #include +#include #include #include #include @@ -72,12 +73,27 @@ struct LJGen_Parameters { double cut = INACTIVE_CUTOFF; double shift = 0.0; double offset = 0.0; + double lambda = 1.0; + double softrad = 0.0; double a1 = 0.0; double a2 = 0.0; double b1 = 0.0; double b2 = 0.0; - double lambda1 = 1.0; - double softrad = 0.0; + LJGen_Parameters() = default; + LJGen_Parameters(double epsilon, double sigma, double cutoff, double shift, + double offset, +#ifdef LJGEN_SOFTCORE + double lam, double delta, +#endif + double e1, double e2, double b1, double b2); + double get_auto_shift() const { + auto auto_shift = 0.; + if (cut != 0.) { + auto_shift = b2 * std::pow(sig / cut, a2) - b1 * std::pow(sig / cut, a1); + } + return auto_shift; + } + double max_cutoff() const { return cut + offset; } }; /** smooth step potential */ diff --git a/src/python/espressomd/interactions.pxd b/src/python/espressomd/interactions.pxd index 53edb5645c4..839de598f65 100644 --- a/src/python/espressomd/interactions.pxd +++ b/src/python/espressomd/interactions.pxd @@ -39,19 +39,6 @@ cdef extern from "TabulatedPotential.hpp": vector[double] force_tab cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp": - cdef struct LJGen_Parameters: - double eps - double sig - double cut - double shift - double offset - double a1 - double a2 - double b1 - double b2 - double lambda1 - double softrad - cdef struct SmoothStep_Parameters: double eps double sig @@ -119,8 +106,6 @@ cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp": cdef struct IA_parameters: - LJGen_Parameters ljgen - SoftSphere_Parameters soft_sphere TabulatedPotential tab @@ -158,20 +143,6 @@ IF THOLE: cdef extern from "nonbonded_interactions/thole.hpp": int thole_set_params(int part_type_a, int part_type_b, double scaling_coeff, double q1q2) -IF LENNARD_JONES_GENERIC: - cdef extern from "nonbonded_interactions/ljgen.hpp": - IF LJGEN_SOFTCORE: - cdef int ljgen_set_params(int part_type_a, int part_type_b, - double eps, double sig, double cut, - double shift, double offset, - double a1, double a2, double b1, double b2, - double genlj_lambda, double softrad) - ELSE: - cdef int ljgen_set_params(int part_type_a, int part_type_b, - double eps, double sig, double cut, - double shift, double offset, - double a1, double a2, double b1, double b2) - IF SMOOTH_STEP: cdef extern from "nonbonded_interactions/smooth_step.hpp": int smooth_step_set_params(int part_type_a, int part_type_b, diff --git a/src/python/espressomd/interactions.pyx b/src/python/espressomd/interactions.pyx index a0c9673495a..9fba81af9a6 100644 --- a/src/python/espressomd/interactions.pyx +++ b/src/python/espressomd/interactions.pyx @@ -409,109 +409,18 @@ IF WCA == 1: IF LENNARD_JONES_GENERIC == 1: - cdef class GenericLennardJonesInteraction(NonBondedInteraction): - - def validate_params(self): - """Check that parameters are valid. - - Raises - ------ - ValueError - If not true. - """ - if self._params["epsilon"] < 0: - raise ValueError("Generic Lennard-Jones epsilon has to be >=0") - if self._params["sigma"] < 0: - raise ValueError("Generic Lennard-Jones sigma has to be >=0") - if self._params["cutoff"] < 0: - raise ValueError("Generic Lennard-Jones cutoff has to be >=0") - IF LJGEN_SOFTCORE: - if self._params["delta"] < 0: - raise ValueError( - "Generic Lennard-Jones delta has to be >=0") - if self._params["lam"] < 0 or self._params["lam"] > 1: - raise ValueError( - "Generic Lennard-Jones lam has to be in the range [0,1]") - - 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.ljgen.eps, - "sigma": ia_params.ljgen.sig, - "cutoff": ia_params.ljgen.cut, - "shift": ia_params.ljgen.shift, - "offset": ia_params.ljgen.offset, - "e1": ia_params.ljgen.a1, - "e2": ia_params.ljgen.a2, - "b1": ia_params.ljgen.b1, - "b2": ia_params.ljgen.b2, - "lam": ia_params.ljgen.lambda1, - "delta": ia_params.ljgen.softrad - } - - def is_active(self): - """Check if interaction is active. - - """ - return (self._params["epsilon"] > 0) - - def _set_params_in_es_core(self): - # Handle the case of shift="auto" - if self._params["shift"] == "auto": - self._params["shift"] = -( - self._params["b1"] * (self._params["sigma"] / self._params["cutoff"])**self._params["e1"] - - self._params["b2"] * (self._params["sigma"] / self._params["cutoff"])**self._params["e2"]) - IF LJGEN_SOFTCORE: - if ljgen_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["e1"], - self._params["e2"], - self._params["b1"], - self._params["b2"], - self._params["lam"], - self._params["delta"]): - raise Exception( - "Could not set Generic Lennard-Jones parameters") - ELSE: - if ljgen_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["e1"], - self._params["e2"], - self._params["b1"], - self._params["b2"], - ): - raise Exception( - "Could not set Generic Lennard-Jones parameters") - - def default_params(self): - """Python dictionary of default parameters. - - """ - IF LJGEN_SOFTCORE: - return {"delta": 0., "lam": 1.} - ELSE: - return {} - - def type_name(self): - """Name of interaction type. - - """ - return "GenericLennardJones" + @script_interface_register + class GenericLennardJonesInteraction(NewNonBondedInteraction): + """ + Generalized Lennard-Jones potential. - def set_params(self, **kwargs): - """ - Set parameters for the generic 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 ---------- @@ -542,8 +451,30 @@ IF LENNARD_JONES_GENERIC == 1: ``LJGEN_SOFTCORE`` parameter lambda. Tune the strength of the interaction. + """ + + _so_name = "Interactions::InteractionLJGen" + + def is_active(self): + """Check if interaction is active. + """ - super().set_params(**kwargs) + return self.epsilon > 0. + + def default_params(self): + """Python dictionary of default parameters. + + """ + IF LJGEN_SOFTCORE: + return {"delta": 0., "lam": 1.} + ELSE: + return {} + + def type_name(self): + """Name of interaction type. + + """ + return "GenericLennardJones" def valid_keys(self): """All parameters that can be set. @@ -1505,9 +1436,6 @@ class NonBondedInteractionHandle(ScriptInterfaceHelper): # Here, add one line for each nonbonded ia IF SOFT_SPHERE: self.soft_sphere = SoftSphereInteraction(_type1, _type2) - IF LENNARD_JONES_GENERIC: - self.generic_lennard_jones = GenericLennardJonesInteraction( - _type1, _type2) IF SMOOTH_STEP: self.smooth_step = SmoothStepInteraction(_type1, _type2) IF BMHTF_NACL: diff --git a/src/script_interface/interactions/NonBondedInteraction.hpp b/src/script_interface/interactions/NonBondedInteraction.hpp index 4929a67f822..7fe6a846c30 100644 --- a/src/script_interface/interactions/NonBondedInteraction.hpp +++ b/src/script_interface/interactions/NonBondedInteraction.hpp @@ -196,6 +196,61 @@ class InteractionLJ : public InteractionPotentialInterface<::LJ_Parameters> { }; #endif // LENNARD_JONES +#ifdef LENNARD_JONES_GENERIC +class InteractionLJGen + : public InteractionPotentialInterface<::LJGen_Parameters> { +protected: + CoreInteraction IA_parameters::*get_ptr_offset() const override { + return &::IA_parameters::ljgen; + } + +public: + InteractionLJGen() { + 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"), +#ifdef LJGEN_SOFTCORE + make_autoparameter(&CoreInteraction::lambda, "lam"), + make_autoparameter(&CoreInteraction::softrad, "delta"), +#endif + make_autoparameter(&CoreInteraction::a1, "e1"), + make_autoparameter(&CoreInteraction::a2, "e2"), + make_autoparameter(&CoreInteraction::b1, "b1"), + make_autoparameter(&CoreInteraction::b2, "b2"), + }); + } + + void make_new_instance(VariantMap const ¶ms) override { + auto new_params = params; + auto const *shift_string = boost::get(¶ms.at("shift")); + if (shift_string != nullptr) { + if (*shift_string != "auto") { + throw std::invalid_argument( + "Generic LJ parameter 'shift' has to be 'auto' or a float"); + } + new_params["shift"] = 0.; + } + m_ia_si = make_shared_from_args( + new_params, "epsilon", "sigma", "cutoff", "shift", "offset", +#ifdef LJGEN_SOFTCORE + "lam", "delta", +#endif + "e1", "e2", "b1", "b2"); + if (shift_string != nullptr) { + m_ia_si->shift = m_ia_si->get_auto_shift(); + } + } +}; +#endif // LENNARD_JONES_GENERIC + #ifdef LJCOS class InteractionLJcos : public InteractionPotentialInterface<::LJcos_Parameters> { @@ -312,6 +367,9 @@ class NonBondedInteractionHandle #ifdef LENNARD_JONES std::shared_ptr m_lj; #endif +#ifdef LENNARD_JONES_GENERIC + std::shared_ptr m_ljgen; +#endif #ifdef LJCOS std::shared_ptr m_ljcos; #endif @@ -348,6 +406,9 @@ class NonBondedInteractionHandle #ifdef LENNARD_JONES make_autoparameter(m_lj, "lennard_jones"), #endif +#ifdef LENNARD_JONES_GENERIC + make_autoparameter(m_ljgen, "generic_lennard_jones"), +#endif #ifdef LJCOS make_autoparameter(m_ljcos, "lennard_jones_cos"), #endif @@ -406,6 +467,10 @@ class NonBondedInteractionHandle set_member(m_lj, "lennard_jones", "Interactions::InteractionLJ", params); #endif +#ifdef LENNARD_JONES_GENERIC + set_member(m_ljgen, "generic_lennard_jones", + "Interactions::InteractionLJGen", params); +#endif #ifdef LJCOS set_member(m_ljcos, "lennard_jones_cos", "Interactions::InteractionLJcos", params); diff --git a/src/script_interface/interactions/initialize.cpp b/src/script_interface/interactions/initialize.cpp index bf1208c89c2..b4359575250 100644 --- a/src/script_interface/interactions/initialize.cpp +++ b/src/script_interface/interactions/initialize.cpp @@ -59,6 +59,9 @@ void initialize(Utils::Factory *om) { #ifdef LENNARD_JONES om->register_new("Interactions::InteractionLJ"); #endif +#ifdef LENNARD_JONES_GENERIC + om->register_new("Interactions::InteractionLJGen"); +#endif #ifdef LJCOS om->register_new("Interactions::InteractionLJcos"); #endif diff --git a/testsuite/python/interactions_non-bonded.py b/testsuite/python/interactions_non-bonded.py index b5eaec6ec55..ba55bfced47 100644 --- a/testsuite/python/interactions_non-bonded.py +++ b/testsuite/python/interactions_non-bonded.py @@ -291,21 +291,32 @@ def assertItemsFractionAlmostEqual(self, a, b): @utx.skipIfMissingFeatures("LENNARD_JONES_GENERIC") def test_lj_generic(self): + params = {"epsilon": 2.12, + "sigma": 1.37, + "cutoff": 2.122, + "offset": 0.185, + "b1": 4.22, + "b2": 3.63, + "e1": 10.32, + "e2": 5.81, + "shift": -0.13} self.run_test("generic_lennard_jones", - {"epsilon": 2.12, - "sigma": 1.37, - "cutoff": 2.122, - "offset": 0.185, - "b1": 4.22, - "b2": 3.63, - "e1": 10.32, - "e2": 5.81, - "shift": -0.13}, + params, force_kernel=tests_common.lj_generic_force, energy_kernel=tests_common.lj_generic_potential, n_steps=231, force_kernel_needs_espressomd=True) + params["shift"] = "auto" + obj = espressomd.interactions.GenericLennardJonesInteraction(**params) + ref_shift = obj.b2 * (obj.sigma / obj.cutoff)**obj.e2 - \ + obj.b1 * (obj.sigma / obj.cutoff)**obj.e1 + self.assertAlmostEqual(obj.shift, ref_shift, delta=1e-10) + + params["cutoff"] = 0. + obj = espressomd.interactions.GenericLennardJonesInteraction(**params) + self.assertEqual(obj.shift, 0.) + # Test WCA Potential @utx.skipIfMissingFeatures("WCA") def test_wca(self): diff --git a/testsuite/python/interactions_non-bonded_interface.py b/testsuite/python/interactions_non-bonded_interface.py index b128fb2b9f6..3efde27b40d 100644 --- a/testsuite/python/interactions_non-bonded_interface.py +++ b/testsuite/python/interactions_non-bonded_interface.py @@ -247,6 +247,25 @@ def test_lj_exceptions(self): {"epsilon": 1., "sigma": 1., "cutoff": 1.5, "shift": 0.2}, ("epsilon", "sigma")) + @utx.skipIfMissingFeatures("LENNARD_JONES_GENERIC") + def test_ljgen_exceptions(self): + check_keys = ("epsilon", "sigma") + params = {"epsilon": 1., "sigma": 1., "cutoff": 1.5, "shift": 0.2, + "offset": 0.1, "b1": 1., "b2": 1., "e1": 12., "e2": 6.} + if espressomd.has_features(["LJGEN_SOFTCORE"]): + params["delta"] = 0.1 + params["lam"] = 0.3 + check_keys = ("delta", "lam", *check_keys) + self.check_potential_exceptions( + espressomd.interactions.GenericLennardJonesInteraction, + params, + check_keys) + with self.assertRaisesRegex(ValueError, f"Generic LJ parameter 'shift' has to be 'auto' or a float"): + invalid_params = params.copy() + invalid_params["shift"] = "unknown" + espressomd.interactions.GenericLennardJonesInteraction( + **invalid_params) + @utx.skipIfMissingFeatures("LJCOS") def test_lj_cos_exceptions(self): self.check_potential_exceptions(