diff --git a/hoomd/hpmc/CMakeLists.txt b/hoomd/hpmc/CMakeLists.txt index 75301d00de..cc9a35abc0 100644 --- a/hoomd/hpmc/CMakeLists.txt +++ b/hoomd/hpmc/CMakeLists.txt @@ -43,6 +43,7 @@ set(_hpmc_sources module.cc PairPotentialExpandedGaussian.cc PairPotentialLJGauss.cc PairPotentialOPP.cc + PairPotentialZetterling.cc PairPotentialStep.cc PairPotentialUnion.cc PairPotentialAngularStep.cc @@ -83,6 +84,7 @@ set(_hpmc_headers PairPotentialExpandedGaussian.h PairPotentialLJGauss.h PairPotentialOPP.h + PairPotentialZetterling.h PairPotentialStep.h PairPotentialUnion.h PairPotentialAngularStep.h diff --git a/hoomd/hpmc/PairPotentialZetterling.cc b/hoomd/hpmc/PairPotentialZetterling.cc new file mode 100644 index 0000000000..b5a929ea1e --- /dev/null +++ b/hoomd/hpmc/PairPotentialZetterling.cc @@ -0,0 +1,107 @@ +// Copyright (c) 2009-2025 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +#include "PairPotentialZetterling.h" + +namespace hoomd + { +namespace hpmc + { + +PairPotentialZetterling::PairPotentialZetterling(std::shared_ptr sysdef) + : PairPotential(sysdef), m_params(m_type_param_index.getNumElements()) + { + } + +LongReal PairPotentialZetterling::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]; + + // Get quantities need for both energy calculation + LongReal r = fast::sqrt(r_squared); + LongReal invr = 1.0 / r; + LongReal eval_cos = fast::cos(2.0 * param.kf * r); + LongReal eval_exp = fast::exp(param.alpha * r); + LongReal eval_invr_3 = invr * invr * invr; + LongReal eval_sigma_over_r = param.sigma * invr; + + // Compute energy + LongReal term1 = param.A * eval_exp * eval_invr_3 * eval_cos; + LongReal term2 = param.B * fast::pow(eval_sigma_over_r, param.n); + LongReal energy = term1 + term2; + + 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 inv_r_cut = 1.0 / r_cut; + LongReal r_cut_eval_cos = fast::cos(2.0 * param.kf * r_cut); + LongReal r_cut_eval_exp = fast::exp(param.alpha * r_cut); + LongReal r_cut_eval_invr_3 = inv_r_cut * inv_r_cut * inv_r_cut; + LongReal r_cut_eval_sigma_over_r = param.sigma * inv_r_cut; + + // Compute energy + LongReal r_cut_term1 = param.A * r_cut_eval_exp * r_cut_eval_invr_3 * r_cut_eval_cos; + LongReal r_cut_term2 = param.B * fast::pow(r_cut_eval_sigma_over_r, param.n); + energy -= r_cut_term1 + r_cut_term2; + } + + 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 PairPotentialZetterling::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 PairPotentialZetterling::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 exportPairPotentialZetterling(pybind11::module& m) + { + pybind11::class_>(m, "PairPotentialZetterling") + .def(pybind11::init>()) + .def("setParams", &PairPotentialZetterling::setParamsPython) + .def("getParams", &PairPotentialZetterling::getParamsPython) + .def_property("mode", &PairPotentialZetterling::getMode, &PairPotentialZetterling::setMode); + } + } // end namespace detail + } // end namespace hpmc + } // end namespace hoomd diff --git a/hoomd/hpmc/PairPotentialZetterling.h b/hoomd/hpmc/PairPotentialZetterling.h new file mode 100644 index 0000000000..a2d46f9a2a --- /dev/null +++ b/hoomd/hpmc/PairPotentialZetterling.h @@ -0,0 +1,149 @@ +// Copyright (c) 2009-2025 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 Zetterling potential energy between two particles. + +For use with HPMC simulations. +*/ +class PairPotentialZetterling : public hpmc::PairPotential + { + public: + PairPotentialZetterling(std::shared_ptr sysdef); + virtual ~PairPotentialZetterling() { } + + 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 + }; + + struct ParamType + { + ParamType() + { + A = 0; + alpha = 0; + kf = 0; + B = 0; + sigma = 0; + n = 0; + } + + ParamType(pybind11::dict v) + { + auto r_cut(v["r_cut"].cast()); + auto r_on(v["r_on"].cast()); + + A = v["A"].cast(); + alpha = v["alpha"].cast(); + kf = v["kf"].cast(); + B = v["B"].cast(); + sigma = v["sigma"].cast(); + n = v["n"].cast(); + r_cut_squared = r_cut * r_cut; + r_on_squared = r_on * r_on; + } + + pybind11::dict asDict() + { + pybind11::dict result; + + result["A"] = A; + result["alpha"] = alpha; + result["kf"] = kf; + result["B"] = B; + result["sigma"] = sigma; + result["n"] = n; + result["r_cut"] = slow::sqrt(r_cut_squared); + result["r_on"] = slow::sqrt(r_on_squared); + + return result; + } + + LongReal A; + LongReal alpha; + LongReal kf; + LongReal B; + LongReal sigma; + LongReal n; + 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 9973bd143c..e0e9c2eb24 100644 --- a/hoomd/hpmc/module.cc +++ b/hoomd/hpmc/module.cc @@ -54,6 +54,8 @@ void exportPairPotentialLJGauss(pybind11::module& m); void exportPairPotentialOPP(pybind11::module& m); +void exportPairPotentialZetterling(pybind11::module& m); + void exportPairPotentialAngularStep(pybind11::module& m); void exportPairPotentialStep(pybind11::module& m); @@ -160,6 +162,7 @@ PYBIND11_MODULE(_hpmc, m) exportPairPotentialExpandedGaussian(m); exportPairPotentialLJGauss(m); exportPairPotentialOPP(m); + exportPairPotentialZetterling(m); exportPairPotentialAngularStep(m); exportPairPotentialStep(m); exportPairPotentialUnion(m); diff --git a/hoomd/hpmc/pair/CMakeLists.txt b/hoomd/hpmc/pair/CMakeLists.txt index 1ae841f254..01cfdf03d7 100644 --- a/hoomd/hpmc/pair/CMakeLists.txt +++ b/hoomd/hpmc/pair/CMakeLists.txt @@ -3,6 +3,7 @@ set(files __init__.py expanded_gaussian.py lj_gauss.py opp.py + zetterling.py pair.py step.py union.py diff --git a/hoomd/hpmc/pair/__init__.py b/hoomd/hpmc/pair/__init__.py index ecd4a3d092..c8d74ab926 100644 --- a/hoomd/hpmc/pair/__init__.py +++ b/hoomd/hpmc/pair/__init__.py @@ -30,11 +30,13 @@ from .expanded_gaussian import ExpandedGaussian from .lj_gauss import LJGauss from .opp import OPP +from .zetterling import Zetterling from .union import Union from .angular_step import AngularStep from .step import Step __all__ = [ + "Zetterling", "OPP", "AngularStep", "ExpandedGaussian", diff --git a/hoomd/hpmc/pair/zetterling.py b/hoomd/hpmc/pair/zetterling.py new file mode 100644 index 0000000000..50fcf3929b --- /dev/null +++ b/hoomd/hpmc/pair/zetterling.py @@ -0,0 +1,118 @@ +# Copyright (c) 2009-2025 The Regents of the University of Michigan. +# Part of HOOMD-blue, released under the BSD 3-Clause License. + +"""Zetterling 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 + + +@hoomd.logging.modify_namespace(("hpmc", "pair", "Zetterling")) +class Zetterling(Pair): + """Zetterling 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. + + `Zetterling` computes the oscillating pair potential between every pair + of particles in the simulation state. The functional form of the potential, + including its behavior under shifting modes. + + .. math:: + U(r) = A \\frac{\\exp{(\\alpha r)\\cos{(2 k_F r)}}}{r^3} + + B \\left( \\frac{\\sigma}{r} \\right)^n + + The potential was introduced in `F. H. M. Zetterling, M. Dzugutov, and S. Lidin + 2001`_. + + .. _F. H. M. Zetterling, M. Dzugutov, and S. Lidin 2001: + https://doi.org/10.1557/PROC-643-K9.5 + + Example:: + + zetterling = hoomd.hpmc.pair.Zetterling(mode="shift") + opp.params[("A", "A")] = { + "A": 1.58, + "alpha": -0.22, + "kf": 4.12, + "B": 0.95533, + "sigma": 1.0, + "n": 18.0, + } + opp.r_cut[("A", "A")] = 2.649 + + {inherited} + + ---------- + + **Members defined in** `Zetterling`: + + .. py:attribute:: params + + The Zetterling potential parameters. The dictionary has the following keys: + + * ``A`` (`float`, **required**) - + Energy scale of the first term :math:`A` + :math:`[\\mathrm{energy}]` + * ``alpha`` (`float`, **required**) - + Screening factor :math:`\alpha` + :math:`[\\mathrm{length}^{-1}]` + * ``kf`` (`float`, **required**) - + Wave number to mimic the Friedel oscillations effect :math:`k_F` + :math:`k_F` :math:`[\\mathrm{length}^{-1}]`. + * ``B`` (`float`, **required**) - + Energy scale of the second term :math:`B` + :math:`B` :math:`[\\mathrm{energy}]`. + * ``sigma`` (`float`, **required**) - + Repulsive core size :math:`\\sigma` :math:`[\\mathrm{length}]` + * ``n`` (`float`, **required**) - + The power to take \\sigma/r in the second term :math:`n` + :math:`[\\mathrm{dimensionless}]` + + Type: `TypeParameter` [`tuple` [``particle_type``, ``particle_type``], + `dict`] + """ + + _cpp_class_name = "PairPotentialZetterling" + + 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( + A=float, + alpha=float, + kf=float, + B=float, + sigma=float, + n=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 abecd740b0..6771878ad9 100644 --- a/hoomd/hpmc/pytest/CMakeLists.txt +++ b/hoomd/hpmc/pytest/CMakeLists.txt @@ -18,6 +18,7 @@ set(files __init__.py test_pair_expanded_gaussian.py test_pair_lj_gauss.py test_pair_opp.py + test_pair_zetterling.py test_pair_step.py test_pair_union.py test_pair_angular_step.py diff --git a/hoomd/hpmc/pytest/test_pair_zetterling.py b/hoomd/hpmc/pytest/test_pair_zetterling.py new file mode 100644 index 0000000000..359a509ddc --- /dev/null +++ b/hoomd/hpmc/pytest/test_pair_zetterling.py @@ -0,0 +1,325 @@ +# Copyright (c) 2009-2025 The Regents of the University of Michigan. +# Part of HOOMD-blue, released under the BSD 3-Clause License. + +"""Test hoomd.hpmc.pair.Zetterling and HPMC pair infrastructure.""" + +import hoomd +import pytest +import numpy as np + +valid_constructor_args = [ + {}, + dict(default_r_cut=3.0), + 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 Zetterling can be constructed with valid arguments.""" + hoomd.hpmc.pair.Zetterling(**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 Zetterling attaches.""" + zetterling = hoomd.hpmc.pair.Zetterling() + zetterling.params[("A", "A")] = dict( + A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0, r_cut=2.649 + ) + + simulation = mc_simulation_factory() + simulation.operations.integrator.pair_potentials = [zetterling] + simulation.run(0) + + assert simulation.operations.integrator._attached + assert zetterling._attached + + simulation.operations.integrator.pair_potentials.remove(zetterling) + assert not zetterling._attached + + +invalid_parameters = [ + {}, + dict(A=1.58), + dict(A=1.58, alpha=-0.22), + dict(A=1.58, alpha=-0.22, kf=4.12), + dict(A=1.58, alpha=-0.22, kf=4.12, B=0.95533), + dict(A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0), + dict(A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0), + dict(A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0, r_cut="invalid"), + dict( + A=1.58, + alpha=-0.22, + kf=4.12, + B=0.95533, + sigma=1.0, + n=18.0, + r_cut=2.649, + r_on="invalid", + ), + dict( + A=1.58, + alpha=-0.22, + kf=4.12, + B=0.95533, + sigma=1.0, + n=18.0, + r_cut=2.649, + r_on=1.0, + invalid=10, + ), +] + + +@pytest.mark.parametrize("parameters", invalid_parameters) +@pytest.mark.cpu +def test_invalid_params_on_attach(mc_simulation_factory, parameters): + """Test that Zetterling validates parameters.""" + zetterling = hoomd.hpmc.pair.Zetterling() + zetterling.params[("A", "A")] = dict( + A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0, r_cut=2.649 + ) + + # Some parameters are validated only after attaching. + simulation = mc_simulation_factory() + simulation.operations.integrator.pair_potentials = [zetterling] + simulation.run(0) + + with pytest.raises( + ( + RuntimeError, + hoomd.error.TypeConversionError, + KeyError, + ) + ): + zetterling.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 vzetterling(r, A, alpha, kf, B, sigma, n): + """Compute Zetterling energy.""" + return A * np.exp(alpha * r) / r**3 * np.cos(2 * kf * r) + B * (sigma / r) ** n + + +# (pair params, +# distance between particles, +# expected energy) +zetterling_test_parameters = [ + ( + dict(A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0, r_cut=2.649), + "none", + 1.0, + vzetterling(r=1.0, A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0), + ), + ( + dict(A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0, r_cut=2.649), + "none", + 1.5, + vzetterling(r=1.5, A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0), + ), + ( + dict(A=1.58, alpha=-0.5, kf=4.12, B=0.95533, sigma=1.0, n=18.0, r_cut=2.649), + "none", + 1.0, + vzetterling(r=1.0, A=1.58, alpha=-0.5, kf=4.12, B=0.95533, sigma=1.0, n=18.0), + ), + ( + dict(A=1.0, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0, r_cut=2.649), + "none", + 1.0, + vzetterling(r=1.0, A=1.0, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0), + ), + ( + dict(A=1.58, alpha=-0.22, kf=3.5, B=0.95533, sigma=1.0, n=18.0, r_cut=2.649), + "none", + 1.0, + vzetterling(r=1.0, A=1.58, alpha=-0.22, kf=3.5, B=0.95533, sigma=1.0, n=18.0), + ), + ( + dict(A=1.58, alpha=-0.22, kf=4.12, B=2, sigma=1.0, n=18.0, r_cut=2.649), + "none", + 1.0, + vzetterling(r=1.0, A=1.58, alpha=-0.22, kf=4.12, B=2, sigma=1.0, n=18.0), + ), + ( + dict(A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=2.0, n=18.0, r_cut=2.649), + "none", + 1.0, + vzetterling(r=1.0, A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=2.0, n=18.0), + ), + ( + dict(A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=15.0, r_cut=2.649), + "none", + 1.0, + vzetterling(r=1.0, A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=15.0), + ), + ( + dict(A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0, r_cut=2.649), + "shift", + 1.0, + vzetterling(r=1.0, A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0) + - vzetterling( + r=2.649, A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0 + ), + ), + ( + dict(A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0, r_cut=2.649), + "shift", + 3.0, + 0.0, + ), + ( + dict( + A=1.58, + alpha=-0.22, + kf=4.12, + B=0.95533, + sigma=1.0, + n=18.0, + r_cut=2.649, + r_on=0.8, + ), + "xplor", + 1.0, + vzetterling(r=1.0, A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0) + * xplor_factor(1.0, 0.8, 2.649), + ), + ( + dict( + A=1.58, + alpha=-0.22, + kf=4.12, + B=0.95533, + sigma=1.0, + n=18.0, + r_cut=2.649, + r_on=0.8, + ), + "xplor", + 1.5, + vzetterling(r=1.5, A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0) + * xplor_factor(1.5, 0.8, 2.649), + ), + ( + dict( + A=1.58, + alpha=-0.22, + kf=4.12, + B=0.95533, + sigma=1.0, + n=18.0, + r_cut=2.649, + r_on=3.0, + ), + "xplor", + 1.0, + vzetterling(r=1.0, A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0) + - vzetterling( + r=2.649, A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0 + ), + ), + ( + dict( + A=1.58, + alpha=-0.22, + kf=4.12, + B=0.95533, + sigma=1.0, + n=18.0, + r_cut=2.649, + r_on=3.0, + ), + "xplor", + 2.8, + 0, + ), +] + + +@pytest.mark.parametrize("params, mode, d, expected_energy", zetterling_test_parameters) +@pytest.mark.cpu +def test_energy(mc_simulation_factory, params, mode, d, expected_energy): + """Test that Zetterling computes the correct energies for 1 pair.""" + zetterling = hoomd.hpmc.pair.Zetterling(mode=mode) + zetterling.params[("A", "A")] = params + + simulation = mc_simulation_factory(d=d) + simulation.operations.integrator.pair_potentials = [zetterling] + simulation.run(0) + + assert zetterling.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.""" + zetterling_1 = hoomd.hpmc.pair.Zetterling() + zetterling_1.params[("A", "A")] = dict( + A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0, r_cut=2.649 + ) + + zetterling_2 = hoomd.hpmc.pair.Zetterling() + zetterling_2.params[("A", "A")] = dict( + A=1.58, alpha=-0.5, kf=4.12, B=0.95533, sigma=1.0, n=18.0, r_cut=2.649 + ) + + expected_1 = vzetterling( + 1.0, A=1.58, alpha=-0.22, kf=4.12, B=0.95533, sigma=1.0, n=18.0 + ) + expected_2 = vzetterling( + 1.0, A=1.58, alpha=-0.5, kf=4.12, B=0.95533, sigma=1.0, n=18.0 + ) + + # Some parameters are validated only after attaching. + simulation = mc_simulation_factory(1.0) + simulation.operations.integrator.pair_potentials = [zetterling_1, zetterling_2] + simulation.run(0) + + assert zetterling_1.energy == pytest.approx(expected=expected_1, rel=1e-5) + assert zetterling_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.Zetterling, + ("hpmc", "pair"), + { + "energy": { + "category": hoomd.logging.LoggerCategories.scalar, + "default": True, + } + }, + ) diff --git a/sphinx-doc/hoomd/hpmc/module-pair.rst b/sphinx-doc/hoomd/hpmc/module-pair.rst index 70487e332c..4a85e76423 100644 --- a/sphinx-doc/hoomd/hpmc/module-pair.rst +++ b/sphinx-doc/hoomd/hpmc/module-pair.rst @@ -3,7 +3,7 @@ pair .. automodule:: hoomd.hpmc.pair :members: - :exclude-members: AngularStep,ExpandedGaussian,LJGauss,LennardJones,OPP,Pair,Step,Union + :exclude-members: AngularStep,ExpandedGaussian,LJGauss,LennardJones,OPP,Pair,Step,Union,Zetterling .. rubric:: Classes @@ -18,3 +18,4 @@ pair pair/pair pair/step pair/union + pair/zetterling