Add Zetterling potential to hpmc pair potential module

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

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<SystemDefinition> sysdef)
+    : PairPotential(sysdef), m_params(m_type_param_index.getNumElements())
+    {
+    }
+
+LongReal PairPotentialZetterling::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];
+
+    // 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<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 PairPotentialZetterling::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 exportPairPotentialZetterling(pybind11::module& m)
+    {
+    pybind11::class_<PairPotentialZetterling,
+                     PairPotential,
+                     std::shared_ptr<PairPotentialZetterling>>(m, "PairPotentialZetterling")
+        .def(pybind11::init<std::shared_ptr<SystemDefinition>>())
+        .def("setParams", &PairPotentialZetterling::setParamsPython)
+        .def("getParams", &PairPotentialZetterling::getParamsPython)
+        .def_property("mode", &PairPotentialZetterling::getMode, &PairPotentialZetterling::setMode);
+    }
+    } // end namespace detail
+    } // end namespace hpmc
+    } // end namespace hoomd

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<SystemDefinition> sysdef);
+    virtual ~PairPotentialZetterling() { }
+
+    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
+        };
+
+    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<LongReal>());
+            auto r_on(v["r_on"].cast<LongReal>());
+
+            A = v["A"].cast<LongReal>();
+            alpha = v["alpha"].cast<LongReal>();
+            kf = v["kf"].cast<LongReal>();
+            B = v["B"].cast<LongReal>();
+            sigma = v["sigma"].cast<LongReal>();
+            n = v["n"].cast<LongReal>();
+            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<ParamType> m_params;
+
+    EnergyShiftMode m_mode = no_shift;
+    };
+
+    } // end namespace hpmc
+    } // end namespace hoomd

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

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

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",