diff --git a/src/core/nonbonded_interactions/hertzian.cpp b/src/core/nonbonded_interactions/hertzian.cpp index 025f513b9e5..698ab9499bf 100644 --- a/src/core/nonbonded_interactions/hertzian.cpp +++ b/src/core/nonbonded_interactions/hertzian.cpp @@ -25,25 +25,13 @@ #include "hertzian.hpp" #ifdef HERTZIAN -#include "interactions.hpp" #include "nonbonded_interaction_data.hpp" -#include - -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 diff --git a/src/core/nonbonded_interactions/hertzian.hpp b/src/core/nonbonded_interactions/hertzian.hpp index 767c401dc6d..89dc3bebc23 100644 --- a/src/core/nonbonded_interactions/hertzian.hpp +++ b/src/core/nonbonded_interactions/hertzian.hpp @@ -31,15 +31,10 @@ #include "nonbonded_interaction_data.hpp" -#include - #include #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) { diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp index 7ba4fdfc48a..1706bb7b0c2 100644 --- a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp +++ b/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 */ diff --git a/src/python/espressomd/interactions.pxd b/src/python/espressomd/interactions.pxd index bd801330582..8c4fef2bbe7 100644 --- a/src/python/espressomd/interactions.pxd +++ b/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, diff --git a/src/python/espressomd/interactions.pyx b/src/python/espressomd/interactions.pyx index 7faff462664..f35fcb08ead 100644 --- a/src/python/espressomd/interactions.pyx +++ 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: diff --git a/src/script_interface/interactions/NonBondedInteraction.hpp b/src/script_interface/interactions/NonBondedInteraction.hpp index bfcef3f4460..01aec3ba3d1 100644 --- a/src/script_interface/interactions/NonBondedInteraction.hpp +++ b/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 ¶ms) override { + m_ia_si = make_shared_from_args( + params, "epsilon", "sigma"); + } +}; +#endif // HERTZIAN + class NonBondedInteractionHandle : public AutoParameters { std::array m_types = {-1, -1}; @@ -299,6 +321,9 @@ class NonBondedInteractionHandle #ifdef GAUSSIAN std::shared_ptr m_gaussian; #endif +#ifdef HERTZIAN + std::shared_ptr m_hertzian; +#endif template auto make_autoparameter(std::shared_ptr &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( m_gaussian, "gaussian", "Interactions::InteractionGaussian", params); +#endif +#ifdef HERTZIAN + set_member( + m_hertzian, "hertzian", "Interactions::InteractionHertzian", params); #endif } diff --git a/src/script_interface/interactions/initialize.cpp b/src/script_interface/interactions/initialize.cpp index cb4ac03bab6..ebf8a228536 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 GAUSSIAN om->register_new("Interactions::InteractionGaussian"); #endif +#ifdef HERTZIAN + om->register_new("Interactions::InteractionHertzian"); +#endif #ifdef LENNARD_JONES om->register_new("Interactions::InteractionLJ"); #endif diff --git a/testsuite/python/interactions_non-bonded.py b/testsuite/python/interactions_non-bonded.py index ae48aff290e..e3bf29d4b4a 100644 --- a/testsuite/python/interactions_non-bonded.py +++ b/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)