Skip to content

Commit

Permalink
Rewrite LJcos2 interaction interface
Browse files Browse the repository at this point in the history
  • Loading branch information
jhossbach authored and jngrad committed Aug 26, 2022
1 parent dbb08ef commit a9f5a24
Show file tree
Hide file tree
Showing 9 changed files with 95 additions and 86 deletions.
39 changes: 15 additions & 24 deletions src/core/nonbonded_interactions/ljcos2.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -25,33 +25,24 @@
#include "ljcos2.hpp"

#ifdef LJCOS2
#include "interactions.hpp"
#include "nonbonded_interaction_data.hpp"

#include <utils/constants.hpp>

#include <cmath>

int ljcos2_set_params(int part_type_a, int part_type_b, double eps, double sig,
double offset, double w) {
IA_parameters *data = get_ia_param_safe(part_type_a, part_type_b);

if (!data)
return ES_ERROR;

data->ljcos2.eps = eps;
data->ljcos2.sig = sig;
data->ljcos2.offset = offset;
data->ljcos2.w = w;

/* calculate dependent parameters */
data->ljcos2.rchange = pow(2, 1 / 6.) * sig;
data->ljcos2.cut = w + data->ljcos2.rchange;

/* broadcast interaction parameters */
mpi_bcast_ia_params(part_type_a, part_type_b);

return ES_OK;
#include <stdexcept>

LJcos2_Parameters::LJcos2_Parameters(double eps, double sig, double offset,
double w)
: eps{eps}, sig{sig}, offset{offset}, w{w} {
if (eps < 0.) {
throw std::domain_error("LJcos2 parameter 'epsilon' has to be >= 0");
}
if (sig < 0.) {
throw std::domain_error("LJcos2 parameter 'sigma' has to be >= 0");
}
if (sig != 0.) {
rchange = std::pow(2., 1. / 6.) * sig;
cut = w + rchange;
}
}

#endif /* ifdef LJCOS2 */
4 changes: 0 additions & 4 deletions src/core/nonbonded_interactions/ljcos2.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -36,16 +36,12 @@

#include "nonbonded_interaction_data.hpp"

#include <utils/Vector.hpp>
#include <utils/constants.hpp>
#include <utils/math/int_pow.hpp>
#include <utils/math/sqr.hpp>

#include <cmath>

int ljcos2_set_params(int part_type_a, int part_type_b, double eps, double sig,
double offset, double w);

/** Calculate Lennard-Jones cosine squared force factor */
inline double ljcos2_pair_force_factor(IA_parameters const &ia_params,
double dist) {
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -192,8 +192,7 @@ static double recalc_maximal_cutoff(const IA_parameters &data) {
#endif

#ifdef LJCOS2
max_cut_current =
std::max(max_cut_current, (data.ljcos2.cut + data.ljcos2.offset));
max_cut_current = std::max(max_cut_current, data.ljcos2.max_cutoff());
#endif

#ifdef GAY_BERNE
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -172,6 +172,9 @@ struct LJcos2_Parameters {
double offset = 0.0;
double w = 0.0;
double rchange = 0.0;
LJcos2_Parameters() = default;
LJcos2_Parameters(double eps, double sig, double offset, double w);
double max_cutoff() const { return cut + offset; }
};

/** Gay-Berne potential */
Expand Down
14 changes: 0 additions & 14 deletions src/python/espressomd/interactions.pxd
Original file line number Diff line number Diff line change
Expand Up @@ -106,12 +106,6 @@ cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp":
double Fmax
double r

cdef struct LJcos2_Parameters:
double eps
double sig
double offset
double w

cdef struct GayBerne_Parameters:
double eps
double sig
Expand All @@ -134,8 +128,6 @@ cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp":

cdef struct IA_parameters:

LJcos2_Parameters ljcos2

LJGen_Parameters ljgen

SoftSphere_Parameters soft_sphere
Expand Down Expand Up @@ -168,12 +160,6 @@ cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp":
cdef void ia_params_set_state(string)
cdef void reset_ia_params()

IF LJCOS2:
cdef extern from "nonbonded_interactions/ljcos2.hpp":
cdef int ljcos2_set_params(int part_type_a, int part_type_b,
double eps, double sig, double offset,
double w)

IF GAY_BERNE:
cdef extern from "nonbonded_interactions/gay_berne.hpp":
int gay_berne_set_params(int part_type_a, int part_type_b,
Expand Down
63 changes: 21 additions & 42 deletions src/python/espressomd/interactions.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -618,38 +618,22 @@ IF LJCOS == 1:
"""
return {"epsilon", "sigma", "cutoff"}

IF LJCOS2:
IF LJCOS2 == 1:

cdef class LennardJonesCos2Interaction(NonBondedInteraction):

def validate_params(self):
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["width"] < 0:
raise ValueError("Parameter width 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.ljcos2.eps,
"sigma": ia_params.ljcos2.sig,
"offset": ia_params.ljcos2.offset,
"width": ia_params.ljcos2.w}

def is_active(self):
return (self._params["epsilon"] > 0)
@script_interface_register
class LennardJonesCos2Interaction(NewNonBondedInteraction):
"""Second variant of the Lennard-Jones cosine interaction.
def set_params(self, **kwargs):
"""Set parameters for the Lennard-Jones cosine squared 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`
Expand All @@ -658,18 +642,12 @@ IF LJCOS2:
Offset distance of the interaction.
width : :obj:`float`
Width of interaction.
"""
super().set_params(**kwargs)
"""

def _set_params_in_es_core(self):
if ljcos2_set_params(self._part_types[0],
self._part_types[1],
self._params["epsilon"],
self._params["sigma"],
self._params["offset"],
self._params["width"]):
raise Exception(
"Could not set Lennard-Jones Cosine2 parameters")
_so_name = "Interactions::InteractionLJcos2"

def is_active(self):
return self.epsilon > 0.

def default_params(self):
return {"offset": 0.}
Expand All @@ -684,14 +662,18 @@ IF LJCOS2:
"""All parameters that can be set.
"""
return {"epsilon", "sigma", "offset", "width"}
return {"epsilon", "sigma", "width", "offset"}

def required_keys(self):
"""Parameters that have to be set.
"""
return {"epsilon", "sigma", "width"}

@property
def cutoff(self):
return self.call_method("get_cutoff")

IF HAT == 1:

cdef class HatInteraction(NonBondedInteraction):
Expand Down Expand Up @@ -1564,9 +1546,6 @@ class NonBondedInteractionHandle(ScriptInterfaceHelper):
IF LENNARD_JONES_GENERIC:
self.generic_lennard_jones = GenericLennardJonesInteraction(
_type1, _type2)
IF LJCOS2:
self.lennard_jones_cos2 = LennardJonesCos2Interaction(
_type1, _type2)
IF SMOOTH_STEP:
self.smooth_step = SmoothStepInteraction(_type1, _type2)
IF BMHTF_NACL:
Expand Down
45 changes: 45 additions & 0 deletions src/script_interface/interactions/NonBondedInteraction.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -222,6 +222,41 @@ class InteractionLJcos
};
#endif // LJCOS

#ifdef LJCOS2
class InteractionLJcos2
: public InteractionPotentialInterface<::LJcos2_Parameters> {
protected:
CoreInteraction IA_parameters::*get_ptr_offset() const override {
return &::IA_parameters::ljcos2;
}

public:
InteractionLJcos2() {
add_parameters({
make_autoparameter(&CoreInteraction::eps, "epsilon"),
make_autoparameter(&CoreInteraction::sig, "sigma"),
make_autoparameter(&CoreInteraction::offset, "offset"),
make_autoparameter(&CoreInteraction::w, "width"),
});
}

void make_new_instance(VariantMap const &params) override {
m_ia_si =
make_shared_from_args<CoreInteraction, double, double, double, double>(
params, "epsilon", "sigma", "offset", "width");
}

Variant do_call_method(std::string const &name,
VariantMap const &params) override {
if (name == "get_cutoff") {
return m_ia_si.get()->cut;
}
return InteractionPotentialInterface<CoreInteraction>::do_call_method(
name, params);
}
};
#endif // LJCOS2

class NonBondedInteractionHandle
: public AutoParameters<NonBondedInteractionHandle> {
std::array<int, 2> m_types = {-1, -1};
Expand All @@ -235,6 +270,9 @@ class NonBondedInteractionHandle
#ifdef LJCOS
std::shared_ptr<InteractionLJcos> m_ljcos;
#endif
#ifdef LJCOS2
std::shared_ptr<InteractionLJcos2> m_ljcos2;
#endif

template <class T>
auto make_autoparameter(std::shared_ptr<T> &member, const char *key) const {
Expand All @@ -261,6 +299,9 @@ class NonBondedInteractionHandle
#endif
#ifdef LJCOS
make_autoparameter(m_ljcos, "lennard_jones_cos"),
#endif
#ifdef LJCOS2
make_autoparameter(m_ljcos2, "lennard_jones_cos2"),
#endif
});
}
Expand Down Expand Up @@ -311,6 +352,10 @@ class NonBondedInteractionHandle
#ifdef LJCOS
set_member<InteractionLJcos>(m_ljcos, "lennard_jones_cos",
"Interactions::InteractionLJcos", params);
#endif
#ifdef LJCOS2
set_member<InteractionLJcos2>(m_ljcos2, "lennard_jones_cos2",
"Interactions::InteractionLJcos2", params);
#endif
}

Expand Down
3 changes: 3 additions & 0 deletions src/script_interface/interactions/initialize.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,9 @@ void initialize(Utils::Factory<ObjectHandle> *om) {
#ifdef LJCOS
om->register_new<InteractionLJcos>("Interactions::InteractionLJcos");
#endif
#ifdef LJCOS2
om->register_new<InteractionLJcos2>("Interactions::InteractionLJcos2");
#endif
#ifdef WCA
om->register_new<InteractionWCA>("Interactions::InteractionWCA");
#endif
Expand Down
7 changes: 7 additions & 0 deletions testsuite/python/interactions_non-bonded_interface.py
Original file line number Diff line number Diff line change
Expand Up @@ -254,6 +254,13 @@ def test_lj_cos_exceptions(self):
{"epsilon": 1., "sigma": 1., "cutoff": 1.5, "offset": 0.2},
("epsilon", "sigma"))

@utx.skipIfMissingFeatures("LJCOS2")
def test_lj_cos2_exceptions(self):
self.check_potential_exceptions(
espressomd.interactions.LennardJonesCos2Interaction,
{"epsilon": 1., "sigma": 1., "width": 1.5, "offset": 0.2},
("epsilon", "sigma"))


if __name__ == "__main__":
ut.main()

0 comments on commit a9f5a24

Please sign in to comment.