Merge pull request #1817 from glotzerlab/Feature/Expanded-Gaussian-Fo…

…r-HPMC

Add expanded gaussian pair for hpmc

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

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<SystemDefinition> sysdef)
+    : PairPotential(sysdef), m_params(m_type_param_index.getNumElements())
+    {
+    }
+
+LongReal PairPotentialExpandedGaussian::energy(const LongReal r_squared,
+                                               const vec3<LongReal>& r_ij,
+                                               const unsigned int type_i,
+                                               const quat<LongReal>& q_i,
+                                               const LongReal charge_i,
+                                               const unsigned int type_j,
+                                               const quat<LongReal>& 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<std::string>());
+    auto type_j = pdata->getTypeByName(typ[1].cast<std::string>());
+    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<std::string>());
+    auto type_j = pdata->getTypeByName(typ[1].cast<std::string>());
+    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_<PairPotentialExpandedGaussian,
+                     PairPotential,
+                     std::shared_ptr<PairPotentialExpandedGaussian>>(
+        m,
+        "PairPotentialExpandedGaussian")
+        .def(pybind11::init<std::shared_ptr<SystemDefinition>>())
+        .def("setParams", &PairPotentialExpandedGaussian::setParamsPython)
+        .def("getParams", &PairPotentialExpandedGaussian::getParamsPython)
+        .def_property("mode",
+                      &PairPotentialExpandedGaussian::getMode,
+                      &PairPotentialExpandedGaussian::setMode);
+    }
+    } // end namespace detail
+    } // end namespace hpmc
+    } // end namespace hoomd

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<SystemDefinition> sysdef);
+    virtual ~PairPotentialExpandedGaussian() { }
+
+    virtual LongReal energy(const LongReal r_squared,
+                            const vec3<LongReal>& r_ij,
+                            const unsigned int type_i,
+                            const quat<LongReal>& q_i,
+                            const LongReal charge_i,
+                            const unsigned int type_j,
+                            const quat<LongReal>& 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<LongReal>());
+            auto r_cut(v["r_cut"].cast<LongReal>());
+            auto r_on(v["r_on"].cast<LongReal>());
+
+            sigma_2 = sigma * sigma;
+            epsilon = v["epsilon"].cast<LongReal>();
+            delta = v["delta"].cast<LongReal>();
+            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<ParamType> m_params;
+
+    EnergyShiftMode m_mode = no_shift;
+    };
+
+    } // end namespace hpmc
+    } // end namespace hoomd

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);

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

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

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