From 302ffc48027fb133b96599e4cddda401115b0cd8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 5 Sep 2022 23:01:18 +0200 Subject: [PATCH] Rewrite BMHTF interaction interface --- .../nonbonded_interactions/bmhtf-nacl.cpp | 47 ++++------ .../nonbonded_interactions/bmhtf-nacl.hpp | 6 +- .../nonbonded_interaction_data.cpp | 2 +- .../nonbonded_interaction_data.hpp | 4 + src/python/espressomd/interactions.pxd | 17 ---- src/python/espressomd/interactions.pyx | 88 ++++++------------- .../interactions/NonBondedInteraction.hpp | 38 ++++++++ .../interactions/initialize.cpp | 3 + .../interactions_non-bonded_interface.py | 8 ++ 9 files changed, 101 insertions(+), 112 deletions(-) diff --git a/src/core/nonbonded_interactions/bmhtf-nacl.cpp b/src/core/nonbonded_interactions/bmhtf-nacl.cpp index 74c2d0df9d0..e89fbf193df 100644 --- a/src/core/nonbonded_interactions/bmhtf-nacl.cpp +++ b/src/core/nonbonded_interactions/bmhtf-nacl.cpp @@ -25,35 +25,26 @@ #include "bmhtf-nacl.hpp" #ifdef BMHTF_NACL -#include "interactions.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" -#include - -int BMHTF_set_params(int part_type_a, int part_type_b, double A, double B, - double C, double D, double sig, double cut) { - double shift, dist2, pw6; - IA_parameters *data = get_ia_param_safe(part_type_a, part_type_b); - - if (!data) - return ES_ERROR; - - dist2 = cut * cut; - pw6 = dist2 * dist2 * dist2; - shift = -(A * exp(B * (sig - cut)) - C / pw6 - D / pw6 / dist2); - - data->bmhtf.A = A; - data->bmhtf.B = B; - data->bmhtf.C = C; - data->bmhtf.D = D; - data->bmhtf.sig = sig; - data->bmhtf.cut = cut; - data->bmhtf.computed_shift = shift; - - /* broadcast interaction parameters */ - mpi_bcast_ia_params(part_type_a, part_type_b); - - return ES_OK; +#include + +#include + +BMHTF_Parameters::BMHTF_Parameters(double a, double b, double c, double d, + double sig, double cutoff) + : A{a}, B{b}, C{c}, D{d}, sig{sig}, cut{cutoff} { + if (a < 0.) { + throw std::domain_error("BMHTF parameter 'a' has to be >= 0"); + } + if (c < 0.) { + throw std::domain_error("BMHTF parameter 'c' has to be >= 0"); + } + if (d < 0.) { + throw std::domain_error("BMHTF parameter 'd' has to be >= 0"); + } + computed_shift = C / Utils::int_pow<6>(cut) + D / Utils::int_pow<8>(cut) - + A * std::exp(B * (sig - cut)); } -#endif +#endif // BMHTF_NACL diff --git a/src/core/nonbonded_interactions/bmhtf-nacl.hpp b/src/core/nonbonded_interactions/bmhtf-nacl.hpp index c2a5b87d110..91ff2179d20 100644 --- a/src/core/nonbonded_interactions/bmhtf-nacl.hpp +++ b/src/core/nonbonded_interactions/bmhtf-nacl.hpp @@ -33,14 +33,10 @@ #include "nonbonded_interaction_data.hpp" -#include #include #include -int BMHTF_set_params(int part_type_a, int part_type_b, double A, double B, - double C, double D, double sig, double cut); - /** Calculate BMHTF force factor */ inline double BMHTF_pair_force_factor(IA_parameters const &ia_params, double dist) { @@ -67,5 +63,5 @@ inline double BMHTF_pair_energy(IA_parameters const &ia_params, double dist) { return 0.0; } -#endif +#endif // BMHTF_NACL #endif diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp index c2801e512ba..2abad44a70f 100644 --- a/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp +++ b/src/core/nonbonded_interactions/nonbonded_interaction_data.cpp @@ -166,7 +166,7 @@ static double recalc_maximal_cutoff(const IA_parameters &data) { #endif #ifdef BMHTF_NACL - max_cut_current = std::max(max_cut_current, data.bmhtf.cut); + max_cut_current = std::max(max_cut_current, data.bmhtf.max_cutoff()); #endif #ifdef MORSE diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp index 922638e3b43..a2fcc9343ed 100644 --- a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp +++ b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp @@ -134,6 +134,10 @@ struct BMHTF_Parameters { double sig = 0.0; double cut = INACTIVE_CUTOFF; double computed_shift = 0.0; + BMHTF_Parameters() = default; + BMHTF_Parameters(double A, double B, double C, double D, double sig, + double cut); + double max_cutoff() const { return cut; } }; /** Morse potential */ diff --git a/src/python/espressomd/interactions.pxd b/src/python/espressomd/interactions.pxd index 839de598f65..f30569e5f29 100644 --- a/src/python/espressomd/interactions.pxd +++ b/src/python/espressomd/interactions.pxd @@ -47,15 +47,6 @@ cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp": int n double k0 - cdef struct BMHTF_Parameters: - double A - double B - double C - double D - double sig - double cut - double computed_shift - cdef struct Morse_Parameters: double eps double alpha @@ -114,8 +105,6 @@ cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp": SmoothStep_Parameters smooth_step - BMHTF_Parameters bmhtf - Morse_Parameters morse Buckingham_Parameters buckingham @@ -149,12 +138,6 @@ IF SMOOTH_STEP: double d, int n, double eps, double k0, double sig, double cut) -IF BMHTF_NACL: - cdef extern from "nonbonded_interactions/bmhtf-nacl.hpp": - int BMHTF_set_params(int part_type_a, int part_type_b, - double A, double B, double C, - double D, double sig, double cut) - IF MORSE: cdef extern from "nonbonded_interactions/morse.hpp": int morse_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 9fba81af9a6..20ad80c9662 100644 --- a/src/python/espressomd/interactions.pyx +++ b/src/python/espressomd/interactions.pyx @@ -938,51 +938,41 @@ IF SMOOTH_STEP == 1: IF BMHTF_NACL == 1: - cdef class BMHTFInteraction(NonBondedInteraction): + @script_interface_register + class BMHTFInteraction(NewNonBondedInteraction): + """BMHTF 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["a"] < 0: - raise ValueError("BMHTF a has to be >=0") - if self._params["c"] < 0: - raise ValueError("BMHTF c has to be >=0") - if self._params["d"] < 0: - raise ValueError("BMHTF d has to be >=0") - if self._params["cutoff"] < 0: - raise ValueError("BMHTF cutoff has to be >=0") + Parameters + ---------- + a : :obj:`float` + Magnitude of exponential part of the interaction. + b : :obj:`float` + Exponential factor of the interaction. + c : :obj:`float` + Magnitude of the term decaying with the sixth power of r. + d : :obj:`float` + Magnitude of the term decaying with the eighth power of r. + sig : :obj:`float` + Shift in the exponent. + cutoff : :obj:`float` + Cutoff distance of the interaction. + """ - 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 { - "a": ia_params.bmhtf.A, - "b": ia_params.bmhtf.B, - "c": ia_params.bmhtf.C, - "d": ia_params.bmhtf.D, - "sig": ia_params.bmhtf.sig, - "cutoff": ia_params.bmhtf.cut, - } + _so_name = "Interactions::InteractionBMHTF" def is_active(self): """Check if interaction is active. """ - return (self._params["a"] > 0) and ( - self._params["c"] > 0) and (self._params["d"] > 0) - - def _set_params_in_es_core(self): - if BMHTF_set_params(self._part_types[0], - self._part_types[1], - self._params["a"], - self._params["b"], - self._params["c"], - self._params["d"], - self._params["sig"], - self._params["cutoff"]): - raise Exception("Could not set BMHTF parameters") + return self.a > 0. and self.c > 0. and self.d > 0. def default_params(self): """Python dictionary of default parameters. @@ -996,28 +986,6 @@ IF BMHTF_NACL == 1: """ return "BMHTF" - def set_params(self, **kwargs): - """ - Set parameters for the BMHTF interaction. - - Parameters - ---------- - a : :obj:`float` - Magnitude of exponential part of the interaction. - b : :obj:`float` - Exponential factor of the interaction. - c : :obj:`float` - Magnitude of the term decaying with the sixth power of r. - d : :obj:`float` - Magnitude of the term decaying with the eighth power of r. - sig : :obj:`float` - Shift in the exponent. - cutoff : :obj:`float` - Cutoff distance of the interaction. - - """ - super().set_params(**kwargs) - def valid_keys(self): """All parameters that can be set. @@ -1438,8 +1406,6 @@ class NonBondedInteractionHandle(ScriptInterfaceHelper): self.soft_sphere = SoftSphereInteraction(_type1, _type2) IF SMOOTH_STEP: self.smooth_step = SmoothStepInteraction(_type1, _type2) - IF BMHTF_NACL: - self.bmhtf = BMHTFInteraction(_type1, _type2) IF MORSE: self.morse = MorseInteraction(_type1, _type2) IF BUCKINGHAM: diff --git a/src/script_interface/interactions/NonBondedInteraction.hpp b/src/script_interface/interactions/NonBondedInteraction.hpp index 7fe6a846c30..ccf4569a8bc 100644 --- a/src/script_interface/interactions/NonBondedInteraction.hpp +++ b/src/script_interface/interactions/NonBondedInteraction.hpp @@ -357,6 +357,34 @@ class InteractionGaussian }; #endif // GAUSSIAN +#ifdef BMHTF_NACL +class InteractionBMHTF + : public InteractionPotentialInterface<::BMHTF_Parameters> { +protected: + CoreInteraction IA_parameters::*get_ptr_offset() const override { + return &::IA_parameters::bmhtf; + } + +public: + InteractionBMHTF() { + add_parameters({ + make_autoparameter(&CoreInteraction::A, "a"), + make_autoparameter(&CoreInteraction::B, "b"), + make_autoparameter(&CoreInteraction::C, "c"), + make_autoparameter(&CoreInteraction::D, "d"), + make_autoparameter(&CoreInteraction::sig, "sig"), + make_autoparameter(&CoreInteraction::cut, "cutoff"), + }); + } + + void make_new_instance(VariantMap const ¶ms) override { + m_ia_si = make_shared_from_args( + params, "a", "b", "c", "d", "sig", "cutoff"); + } +}; +#endif // BMHTF_NACL + class NonBondedInteractionHandle : public AutoParameters { std::array m_types = {-1, -1}; @@ -382,6 +410,9 @@ class NonBondedInteractionHandle #ifdef GAUSSIAN std::shared_ptr m_gaussian; #endif +#ifdef BMHTF_NACL + std::shared_ptr m_bmhtf; +#endif template auto make_autoparameter(std::shared_ptr &member, const char *key) const { @@ -420,6 +451,9 @@ class NonBondedInteractionHandle #endif #ifdef GAUSSIAN make_autoparameter(m_gaussian, "gaussian"), +#endif +#ifdef BMHTF_NACL + make_autoparameter(m_bmhtf, "bmhtf"), #endif }); } @@ -486,6 +520,10 @@ class NonBondedInteractionHandle #ifdef GAUSSIAN set_member( m_gaussian, "gaussian", "Interactions::InteractionGaussian", params); +#endif +#ifdef BMHTF_NACL + set_member(m_bmhtf, "bmhtf", + "Interactions::InteractionBMHTF", params); #endif } diff --git a/src/script_interface/interactions/initialize.cpp b/src/script_interface/interactions/initialize.cpp index b4359575250..ab4c766f575 100644 --- a/src/script_interface/interactions/initialize.cpp +++ b/src/script_interface/interactions/initialize.cpp @@ -77,6 +77,9 @@ void initialize(Utils::Factory *om) { #ifdef GAUSSIAN om->register_new("Interactions::InteractionGaussian"); #endif +#ifdef BMHTF_NACL + om->register_new("Interactions::InteractionBMHTF"); +#endif } } // namespace Interactions } // namespace ScriptInterface diff --git a/testsuite/python/interactions_non-bonded_interface.py b/testsuite/python/interactions_non-bonded_interface.py index 66f69ff7352..184a9e6ad4f 100644 --- a/testsuite/python/interactions_non-bonded_interface.py +++ b/testsuite/python/interactions_non-bonded_interface.py @@ -289,6 +289,14 @@ def test_gaussian_exceptions(self): ("eps", "sig") ) + @utx.skipIfMissingFeatures("BMHTF_NACL") + def test_bmhtf_exceptions(self): + self.check_potential_exceptions( + espressomd.interactions.BMHTFInteraction, + {"a": 3., "b": 2., "c": 1., "d": 4., "sig": 0.13, "cutoff": 1.2}, + ("a", "c", "d") + ) + if __name__ == "__main__": ut.main()