Skip to content

Commit

Permalink
Rewrite generic Lennard-Jones interaction interface
Browse files Browse the repository at this point in the history
  • Loading branch information
@jngrad
jngrad committed Sep 5, 2022
1 parent a42a586 commit 8576bf8
Show file tree
Hide file tree
Showing 10 changed files with 191 additions and 195 deletions.
55 changes: 24 additions & 31 deletions src/core/nonbonded_interactions/ljgen.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,43 +25,36 @@
#include "ljgen.hpp"

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

#include <utils/constants.hpp>

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 <stdexcept>

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
23 changes: 7 additions & 16 deletions src/core/nonbonded_interactions/ljgen.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,24 +47,15 @@

#include <cmath>

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
Expand All @@ -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 *
Expand All @@ -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
Expand All @@ -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) +
Expand All @@ -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
3 changes: 1 addition & 2 deletions src/core/nonbonded_interactions/nonbonded_interaction_data.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
20 changes: 18 additions & 2 deletions src/core/nonbonded_interactions/nonbonded_interaction_data.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@

#include <algorithm>
#include <cassert>
#include <cmath>
#include <memory>
#include <string>
#include <vector>
Expand Down Expand Up @@ -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 */
Expand Down
29 changes: 0 additions & 29 deletions src/python/espressomd/interactions.pxd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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,
Expand Down
140 changes: 34 additions & 106 deletions src/python/espressomd/interactions.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
----------
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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:
Expand Down
Loading

0 comments on commit 8576bf8

Please sign in to comment.