diff --git a/hoomd/hpmc/CMakeLists.txt b/hoomd/hpmc/CMakeLists.txt index b3b65a9de5..c5af5e08cd 100644 --- a/hoomd/hpmc/CMakeLists.txt +++ b/hoomd/hpmc/CMakeLists.txt @@ -41,6 +41,7 @@ set(_hpmc_sources module.cc ExternalPotentialLinear.cc PairPotential.cc PairPotentialLennardJones.cc + PairPotentialExpandedGaussian.cc PairPotentialStep.cc PairPotentialUnion.cc PairPotentialAngularStep.cc @@ -85,6 +86,7 @@ set(_hpmc_headers OBBTree.h PairPotential.h PairPotentialLennardJones.h + PairPotentialExpandedGaussian.h PairPotentialStep.h PairPotentialUnion.h PairPotentialAngularStep.h diff --git a/hoomd/hpmc/PairPotentialExpandedGaussian.cc b/hoomd/hpmc/PairPotentialExpandedGaussian.cc new file mode 100644 index 0000000000..9d8ee2063c --- /dev/null +++ b/hoomd/hpmc/PairPotentialExpandedGaussian.cc @@ -0,0 +1,98 @@ +// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +#include "PairPotentialExpandedGaussian.h" + +namespace hoomd + { +namespace hpmc + { + +PairPotentialExpandedGaussian::PairPotentialExpandedGaussian( + std::shared_ptr sysdef) + : PairPotential(sysdef), m_params(m_type_param_index.getNumElements()) + { + } + +LongReal PairPotentialExpandedGaussian::energy(const LongReal r_squared, + const vec3& r_ij, + const unsigned int type_i, + const quat& q_i, + const LongReal charge_i, + const unsigned int type_j, + const quat& q_j, + const LongReal charge_j) const + { + unsigned int param_index = m_type_param_index(type_i, type_j); + const auto& param = m_params[param_index]; + + LongReal r = fast::sqrt(r_squared); + LongReal rmd_2 = (r - param.delta) * (r - param.delta); + LongReal rmd_over_sigma_2 = rmd_2 / param.sigma_2; + LongReal exp_val = fast::exp(-LongReal(1.0) / LongReal(2.0) * rmd_over_sigma_2); + LongReal energy = param.epsilon * exp_val; + + if (m_mode == shift || (m_mode == xplor && param.r_on_squared >= param.r_cut_squared)) + { + LongReal r_cut = fast::sqrt(param.r_cut_squared); + LongReal rcutmd_2 = (r_cut - param.delta) * (r_cut - param.delta); + LongReal rcutmd_over_sigma_2 = rcutmd_2 / param.sigma_2; + energy -= param.epsilon * fast::exp(-LongReal(1.0) / LongReal(2.0) * rcutmd_over_sigma_2); + } + + if (m_mode == xplor && r_squared > param.r_on_squared) + { + LongReal a = param.r_cut_squared - param.r_on_squared; + LongReal denominator = a * a * a; + + LongReal b = param.r_cut_squared - r_squared; + LongReal numerator = b * b + * (param.r_cut_squared + LongReal(2.0) * r_squared + - LongReal(3.0) * param.r_on_squared); + energy *= numerator / denominator; + } + + return energy; + } + +void PairPotentialExpandedGaussian::setParamsPython(pybind11::tuple typ, pybind11::dict params) + { + auto pdata = m_sysdef->getParticleData(); + auto type_i = pdata->getTypeByName(typ[0].cast()); + auto type_j = pdata->getTypeByName(typ[1].cast()); + unsigned int param_index_1 = m_type_param_index(type_i, type_j); + m_params[param_index_1] = ParamType(params); + unsigned int param_index_2 = m_type_param_index(type_j, type_i); + m_params[param_index_2] = ParamType(params); + + notifyRCutChanged(); + } + +pybind11::dict PairPotentialExpandedGaussian::getParamsPython(pybind11::tuple typ) + { + auto pdata = m_sysdef->getParticleData(); + auto type_i = pdata->getTypeByName(typ[0].cast()); + auto type_j = pdata->getTypeByName(typ[1].cast()); + unsigned int param_index = m_type_param_index(type_i, type_j); + return m_params[param_index].asDict(); + } + +namespace detail + { +void exportPairPotentialExpandedGaussian(pybind11::module& m) + { + pybind11::class_>( + m, + "PairPotentialExpandedGaussian") + .def(pybind11::init>()) + .def("setParams", &PairPotentialExpandedGaussian::setParamsPython) + .def("getParams", &PairPotentialExpandedGaussian::getParamsPython) + .def_property("mode", + &PairPotentialExpandedGaussian::getMode, + &PairPotentialExpandedGaussian::setMode); + } + } // end namespace detail + } // end namespace hpmc + } // end namespace hoomd diff --git a/hoomd/hpmc/PairPotentialExpandedGaussian.h b/hoomd/hpmc/PairPotentialExpandedGaussian.h new file mode 100644 index 0000000000..dd237d8826 --- /dev/null +++ b/hoomd/hpmc/PairPotentialExpandedGaussian.h @@ -0,0 +1,141 @@ +// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +#pragma once + +#include "PairPotential.h" + +namespace hoomd + { +namespace hpmc + { + +/*** Compute Expanded Gaussian energy between two particles. + +For use with HPMC simulations. +*/ +class PairPotentialExpandedGaussian : public hpmc::PairPotential + { + public: + PairPotentialExpandedGaussian(std::shared_ptr sysdef); + virtual ~PairPotentialExpandedGaussian() { } + + virtual LongReal energy(const LongReal r_squared, + const vec3& r_ij, + const unsigned int type_i, + const quat& q_i, + const LongReal charge_i, + const unsigned int type_j, + const quat& q_j, + const LongReal charge_j) const; + + /// Compute the non-additive cuttoff radius + virtual LongReal computeRCutNonAdditive(unsigned int type_i, unsigned int type_j) const + { + unsigned int param_index = m_type_param_index(type_i, type_j); + return slow::sqrt(m_params[param_index].r_cut_squared); + } + + /// Set type pair dependent parameters to the potential. + void setParamsPython(pybind11::tuple typ, pybind11::dict params); + + /// Get type pair dependent parameters. + pybind11::dict getParamsPython(pybind11::tuple typ); + + void setMode(const std::string& mode_str) + { + if (mode_str == "none") + { + m_mode = no_shift; + } + else if (mode_str == "shift") + { + m_mode = shift; + } + else if (mode_str == "xplor") + { + m_mode = xplor; + } + else + { + throw std::domain_error("Invalid mode " + mode_str); + } + } + + std::string getMode() + { + std::string result = "none"; + + if (m_mode == shift) + { + result = "shift"; + } + if (m_mode == xplor) + { + result = "xplor"; + } + + return result; + } + + protected: + /// Shifting modes that can be applied to the energy + enum EnergyShiftMode + { + no_shift = 0, + shift, + xplor + }; + + /// Type pair parameters of EG potential + struct ParamType + { + ParamType() + { + sigma_2 = 0; + epsilon = 0; + delta = 0; + r_cut_squared = 0; + r_on_squared = 0; + } + + ParamType(pybind11::dict v) + { + auto sigma(v["sigma"].cast()); + auto r_cut(v["r_cut"].cast()); + auto r_on(v["r_on"].cast()); + + sigma_2 = sigma * sigma; + epsilon = v["epsilon"].cast(); + delta = v["delta"].cast(); + r_cut_squared = r_cut * r_cut; + r_on_squared = r_on * r_on; + } + + pybind11::dict asDict() + { + pybind11::dict result; + + result["sigma"] = pow(sigma_2, 1. / 2.); + result["epsilon"] = epsilon; + result["delta"] = delta; + result["r_cut"] = slow::sqrt(r_cut_squared); + result["r_on"] = slow::sqrt(r_on_squared); + + return result; + } + + LongReal sigma_2; + LongReal epsilon; + LongReal delta; + LongReal r_cut_squared; + LongReal r_on_squared; + }; + + std::vector m_params; + + EnergyShiftMode m_mode = no_shift; + }; + + } // end namespace hpmc + } // end namespace hoomd diff --git a/hoomd/hpmc/module.cc b/hoomd/hpmc/module.cc index cf29b53800..0aaf99effd 100644 --- a/hoomd/hpmc/module.cc +++ b/hoomd/hpmc/module.cc @@ -49,6 +49,8 @@ void exportPairPotential(pybind11::module& m); void exportPairPotentialLennardJones(pybind11::module& m); +void exportPairPotentialExpandedGaussian(pybind11::module& m); + void exportPairPotentialAngularStep(pybind11::module& m); void exportPairPotentialStep(pybind11::module& m); @@ -154,6 +156,7 @@ PYBIND11_MODULE(_hpmc, m) exportPairPotential(m); exportPairPotentialLennardJones(m); + exportPairPotentialExpandedGaussian(m); exportPairPotentialAngularStep(m); exportPairPotentialStep(m); exportPairPotentialUnion(m); diff --git a/hoomd/hpmc/pair/CMakeLists.txt b/hoomd/hpmc/pair/CMakeLists.txt index 5d1601d8f6..4170d377f3 100644 --- a/hoomd/hpmc/pair/CMakeLists.txt +++ b/hoomd/hpmc/pair/CMakeLists.txt @@ -1,5 +1,6 @@ set(files __init__.py lennard_jones.py + expanded_gaussian.py pair.py step.py union.py diff --git a/hoomd/hpmc/pair/__init__.py b/hoomd/hpmc/pair/__init__.py index db1fe0bbe2..ffb72ccbb1 100644 --- a/hoomd/hpmc/pair/__init__.py +++ b/hoomd/hpmc/pair/__init__.py @@ -30,6 +30,7 @@ from . import user from .pair import Pair from .lennard_jones import LennardJones +from .expanded_gaussian import ExpandedGaussian from .union import Union from .angular_step import AngularStep from .step import Step diff --git a/hoomd/hpmc/pair/expanded_gaussian.py b/hoomd/hpmc/pair/expanded_gaussian.py new file mode 100644 index 0000000000..8f611259f7 --- /dev/null +++ b/hoomd/hpmc/pair/expanded_gaussian.py @@ -0,0 +1,107 @@ +# Copyright (c) 2009-2024 The Regents of the University of Michigan. +# Part of HOOMD-blue, released under the BSD 3-Clause License. + +"""Expanded Gaussian pair potential. + +.. invisible-code-block: python + + simulation = hoomd.util.make_example_simulation() + sphere = hoomd.hpmc.integrate.Sphere() + sphere.shape['A'] = dict(diameter=0.0) + simulation.operations.integrator = sphere +""" + +import hoomd + +from .pair import Pair +from hoomd.data.typeconverter import positive_real + + +@hoomd.logging.modify_namespace(('hpmc', 'pair', 'ExpandedGaussian')) +class ExpandedGaussian(Pair): + """Expanded Gaussian pair potential (HPMC). + + Args: + default_r_cut (float): Default cutoff radius :math:`[\\mathrm{length}]`. + default_r_on (float): Default XPLOR on radius + :math:`[\\mathrm{length}]`. + mode (str): Energy shifting/smoothing mode. + + `ExpandedGaussian` computes the Expanded Gaussian pair potential between + every pair of particles in the simulation state. The functional form of the + potential, including its behavior under shifting modes, is identical to that + in the MD pair potential `hoomd.md.pair.ExpandedGaussian`. + + See Also: + `hoomd.md.pair.ExpandedGaussian` + + `hoomd.md.pair` + + .. rubric:: Example + + .. code-block:: python + + expanded_gaussian = hoomd.hpmc.pair.ExpandedGaussian() + expanded_gaussian.params[('A', 'A')] = dict(epsilon=1.0, + sigma=1.0, + delta=1.0, + r_cut=2.5) + simulation.operations.integrator.pair_potentials = [expanded_gaussian] + + .. py:attribute:: params + + The potential parameters. The dictionary has the following keys: + + * ``epsilon`` (`float`, **required**) - + Energy well depth :math:`\\varepsilon` :math:`[\\mathrm{energy}]`. + * ``sigma`` (`float`, **required**) - + Characteristic length scale :math:`\\sigma` + :math:`[\\mathrm{length}]`. + * ``delta`` (`float`, **required**) - + Characteristic length scale :math:`\\delta` + :math:`[\\mathrm{length}]`. + * ``r_cut`` (`float`): Cutoff radius :math:`[\\mathrm{length}]`. + Defaults to the value given in ``default_r_cut`` on construction. + * ``r_on`` (`float`): XPLOR on radius :math:`[\\mathrm{length}]`. + Defaults to the value given in ``default_r_on`` on construction. + + Type: `TypeParameter` [`tuple` [``particle_type``, ``particle_type``], + `dict`] + + .. py:attribute:: mode + + The energy shifting/smoothing mode: Possible values are: + ``"none"``, ``"shift"``, and ``"xplor"``. + + .. rubric:: Example + + .. code-block:: python + + expanded_gaussian.mode = 'shift' + + Type: `str` + """ + _cpp_class_name = "PairPotentialExpandedGaussian" + + def __init__(self, default_r_cut=None, default_r_on=0.0, mode='none'): + if default_r_cut is None: + default_r_cut = float + else: + default_r_cut = float(default_r_cut) + + params = hoomd.data.typeparam.TypeParameter( + 'params', 'particle_types', + hoomd.data.parameterdicts.TypeParameterDict( + epsilon=float, + sigma=positive_real, + delta=float, + r_cut=default_r_cut, + r_on=float(default_r_on), + len_keys=2)) + self._add_typeparam(params) + + self._param_dict.update( + hoomd.data.parameterdicts.ParameterDict( + mode=hoomd.data.typeconverter.OnlyFrom(("none", "shift", + "xplor")))) + self.mode = mode diff --git a/hoomd/hpmc/pytest/CMakeLists.txt b/hoomd/hpmc/pytest/CMakeLists.txt index e085c9603d..2bfff5d32c 100644 --- a/hoomd/hpmc/pytest/CMakeLists.txt +++ b/hoomd/hpmc/pytest/CMakeLists.txt @@ -14,6 +14,7 @@ set(files __init__.py test_shape_utils.py test_move_size_tuner.py test_pair_lennard_jones.py + test_pair_expanded_gaussian.py test_pair_step.py test_pair_user.py test_pair_union.py diff --git a/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py new file mode 100644 index 0000000000..e3040e8b5d --- /dev/null +++ b/hoomd/hpmc/pytest/test_pair_expanded_gaussian.py @@ -0,0 +1,222 @@ +# Copyright (c) 2009-2024 The Regents of the University of Michigan. +# Part of HOOMD-blue, released under the BSD 3-Clause License. + +"""Test hoomd.hpmc.pair.ExpandedGaussian and HPMC pair infrastructure.""" + +import hoomd +import pytest +import math + +valid_constructor_args = [ + {}, + dict(default_r_cut=2.5), + dict(default_r_on=2.0), + dict(mode='shift'), +] + + +@pytest.mark.parametrize("constructor_args", valid_constructor_args) +def test_valid_construction(device, constructor_args): + """Test that ExpandedGaussian can be constructed with valid arguments.""" + hoomd.hpmc.pair.ExpandedGaussian(**constructor_args) + + +@pytest.fixture(scope='session') +def mc_simulation_factory(simulation_factory, two_particle_snapshot_factory): + """Make a MC simulation with two particles separate dy by a distance d.""" + + def make_simulation(d=1): + snapshot = two_particle_snapshot_factory(d=d) + simulation = simulation_factory(snapshot) + + sphere = hoomd.hpmc.integrate.Sphere() + sphere.shape['A'] = dict(diameter=0) + simulation.operations.integrator = sphere + + return simulation + + return make_simulation + + +@pytest.mark.cpu +def test_attaching(mc_simulation_factory): + """Test that ExpandedGaussian attaches.""" + expanded_gauss = hoomd.hpmc.pair.ExpandedGaussian() + expanded_gauss.params[('A', 'A')] = dict(epsilon=1.0, + sigma=1.0, + delta=1.0, + r_cut=2.5) + + simulation = mc_simulation_factory() + simulation.operations.integrator.pair_potentials = [expanded_gauss] + simulation.run(0) + + assert simulation.operations.integrator._attached + assert expanded_gauss._attached + + simulation.operations.integrator.pair_potentials.remove(expanded_gauss) + assert not expanded_gauss._attached + + +invalid_parameters = [ + {}, + dict(epsilon=1.0), + dict(epsilon=1.0, sigma=1.0), + dict(epsilon=1.0, sigma=1.0, delta=1.0, r_cut='invalid'), + dict(epsilon=1.0, sigma=1.0, delta=1.0, r_cut=2.5, r_on='invalid'), + dict(epsilon=1.0, sigma=1.0, delta=1.0, r_cut=2.5, r_on=2.0, invalid=10), +] + + +@pytest.mark.parametrize("parameters", invalid_parameters) +@pytest.mark.cpu +def test_invalid_params_on_attach(mc_simulation_factory, parameters): + """Test that ExpandedGaussian validates parameters.""" + expanded_gauss = hoomd.hpmc.pair.ExpandedGaussian() + expanded_gauss.params[('A', 'A')] = dict(epsilon=1.0, + sigma=1.0, + delta=1.0, + r_cut=2.5) + + # Some parameters are validated only after attaching. + simulation = mc_simulation_factory() + simulation.operations.integrator.pair_potentials = [expanded_gauss] + simulation.run(0) + + with pytest.raises(( + RuntimeError, + hoomd.error.TypeConversionError, + KeyError, + )): + expanded_gauss.params[('A', 'A')] = parameters + + +def xplor_factor(r, r_on, r_cut): + """Compute the XPLOR smoothing factor.""" + if r < r_on: + return 1 + if r < r_cut: + denominator = (r_cut**2 - r_on**2)**3 + numerator = (r_cut**2 - r**2)**2 * (r_cut**2 + 2 * r**2 - 3 * r_on**2) + return numerator / denominator + + return 0 + + +def eg(r, epsilon, sigma, delta): + """Compute the eg energy.""" + return epsilon * math.exp(-0.5 * (((r - delta) / sigma)**2)) + + +# (pair params, +# mode, +# distance between particles, +# expected energy) +expanded_gauss_test_parameters = [ + ( + dict(epsilon=2.0, sigma=1.5, delta=1.0, r_cut=2.5), + 'none', + 3.0, + 0.0, + ), + ( + dict(epsilon=2.0, sigma=1.5, delta=1.0, r_cut=2.5), + 'none', + 1.0, + 2.0, + ), + ( + dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=2.0), + 'none', + 1.5, + eg(1.5, 5, 1.1, 1.0), + ), + ( + dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=2.0), + 'shift', + 1.5, + eg(1.5, 5, 1.1, 1.0) - eg(2.0, 5, 1.1, 1.0), + ), + ( + dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=2.5, r_on=2.0), + 'xplor', + 1.5, + eg(1.5, 5.0, 1.1, 1.0) * xplor_factor(1.5, 2.0, 2.5), + ), + ( + dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=2.5, r_on=2.0), + 'xplor', + 2.3, + eg(2.3, 5, 1.1, 1.0) * xplor_factor(2.3, 2.0, 2.5), + ), + ( + dict(epsilon=5.0, sigma=1.1, delta=1.0, r_cut=2.0, r_on=3), + 'xplor', + 1.5, + eg(1.5, 5, 1.1, 1.0) - eg(2, 5, 1.1, 1.0), + ), + ( + dict(epsilon=1.0, sigma=1, delta=1.0, r_cut=3.0, r_on=4), + 'xplor', + 3.2, + 0, + ), +] + + +@pytest.mark.parametrize('params, mode, d, expected_energy', + expanded_gauss_test_parameters) +@pytest.mark.cpu +def test_energy(mc_simulation_factory, params, mode, d, expected_energy): + """Test that ExpandedGaussian computes the correct energies for 1 pair.""" + expanded_gauss = hoomd.hpmc.pair.ExpandedGaussian(mode=mode) + expanded_gauss.params[('A', 'A')] = params + + simulation = mc_simulation_factory(d=d) + simulation.operations.integrator.pair_potentials = [expanded_gauss] + simulation.run(0) + + assert expanded_gauss.energy == pytest.approx(expected=expected_energy, + rel=1e-5) + + +@pytest.mark.cpu +def test_multiple_pair_potentials(mc_simulation_factory): + """Test that energy operates correctly with multiple pair potentials.""" + expanded_gauss_1 = hoomd.hpmc.pair.ExpandedGaussian() + expanded_gauss_1.params[('A', 'A')] = dict(epsilon=1.0, + sigma=1.0, + delta=1.0, + r_cut=2.5) + + expanded_gauss_2 = hoomd.hpmc.pair.ExpandedGaussian() + expanded_gauss_2.params[('A', 'A')] = dict(epsilon=2.0, + sigma=1.0, + delta=1.0, + r_cut=2.5) + expected_1 = eg(1.5, 1.0, 1.0, 1.0) + expected_2 = eg(1.5, 2.0, 1.0, 1.0) + + # Some parameters are validated only after attaching. + simulation = mc_simulation_factory(1.5) + simulation.operations.integrator.pair_potentials = [ + expanded_gauss_1, expanded_gauss_2 + ] + simulation.run(0) + + assert expanded_gauss_1.energy == pytest.approx(expected=expected_1, + rel=1e-5) + assert expanded_gauss_2.energy == pytest.approx(expected=expected_2, + rel=1e-5) + assert simulation.operations.integrator.pair_energy == pytest.approx( + expected=expected_1 + expected_2, rel=1e-5) + + +def test_logging(): + hoomd.conftest.logging_check( + hoomd.hpmc.pair.ExpandedGaussian, ('hpmc', 'pair'), { + 'energy': { + 'category': hoomd.logging.LoggerCategories.scalar, + 'default': True + } + }) diff --git a/sphinx-doc/module-hpmc-pair.rst b/sphinx-doc/module-hpmc-pair.rst index af00e3ce67..05e643ffc8 100644 --- a/sphinx-doc/module-hpmc-pair.rst +++ b/sphinx-doc/module-hpmc-pair.rst @@ -12,6 +12,7 @@ hoomd.hpmc.pair :nosignatures: AngularStep + ExpandedGaussian LennardJones Pair Step @@ -23,6 +24,7 @@ hoomd.hpmc.pair :synopsis: Pair potentials for HPMC. :show-inheritance: :members: AngularStep, + ExpandedGaussian, LennardJones, Pair, Step,