Added adjusted GCMS potential

doc/module-azplugins-pair.rst

@@ -12,6 +12,7 @@ azplugins.pair
 .. autosummary::
     :nosignatures:
 
+    ACSW
     Colloid
     DPDGeneralWeight
     Hertz
@@ -22,7 +23,8 @@ azplugins.pair
 
 .. automodule:: hoomd.azplugins.pair
     :synopsis: Pair potentials.
-    :members: Colloid,
+    :members: ACSW,
+        Colloid,
         DPDGeneralWeight,
         Hertz,
         PerturbedLennardJones,

src/CMakeLists.txt

@@ -32,6 +32,7 @@ set(_pair_evaluators
     Colloid
     Hertz
     PerturbedLennardJones
+    ACSW
     )
 
 # TODO: Add names of all dpd evaluators

src/PairEvaluatorACSW.h

@@ -0,0 +1,144 @@
+// Copyright (c) 2018-2020, Michael P. Howard
+// Copyright (c) 2021-2024, Auburn University
+// Part of azplugins, released under the BSD 3-Clause License.
+
+#ifndef AZPLUGINS_PAIR_EVALUATOR_ACSW_H_
+#define AZPLUGINS_PAIR_EVALUATOR_ACSW_H_
+
+#include "PairEvaluator.h"
+#include <stdio.h>
+
+#ifdef __HIPCC__
+#define DEVICE __device__
+#else
+#define DEVICE
+#endif
+
+namespace hoomd
+    {
+namespace azplugins
+    {
+namespace detail
+    {
+
+struct PairParametersACSW : public PairParameters
+    {
+#ifndef __HIPCC__
+    PairParametersACSW() : w(0), sigma(0), a(0), q(0) { }
+
+    PairParametersACSW(pybind11::dict v, bool managed = false)
+        {
+        w = v["w"].cast<Scalar>();
+        sigma = v["sigma"].cast<Scalar>();
+        a = v["a"].cast<Scalar>();
+        q = v["q"].cast<Scalar>();
+        }
+
+    pybind11::dict asDict()
+        {
+        pybind11::dict v;
+        v["w"] = w;
+        v["sigma"] = sigma;
+        v["a"] = a;
+        v["q"] = q;
+        return v;
+        }
+#endif // __HIPCC__
+
+    Scalar w;     //!< well width [length]
+    Scalar sigma; //!< Hard core repulsion distance [length]
+    Scalar a;     //!< Well depth [energy]
+    Scalar q;     //!< steepness
+    }
+#if HOOMD_LONGREAL_SIZE == 32
+    __attribute__((aligned(16)));
+#else
+    __attribute__((aligned(32)));
+#endif
+
+//! Class for evaluating the adjusted generalized continuous multiple step pair potential
+/*!
+ * This is the generalized continuous multiple step pair potential
+ * modified such that there is only one sigma and one minimum present
+ * with any combination of well depth and width causing the potential
+ * to have a value of zero at sigma.
+ */
+class PairEvaluatorACSW : public PairEvaluator
+    {
+    public:
+    typedef PairParametersACSW param_type;
+
+    //! Constructor
+    /*!
+     * \param _rsq Squared distance between particles
+     * \param _rcutsq Cutoff radius squared
+     * \param _params Pair potential parameters, given by typedef above
+     *
+     * The functor initializes its members from \a _params.
+     */
+    DEVICE PairEvaluatorACSW(Scalar _rsq, Scalar _rcutsq, const param_type& _params)
+        : PairEvaluator(_rsq, _rcutsq)
+        {
+        w = _params.w;
+        sigma = _params.sigma;
+        a = _params.a;
+        q = _params.q;
+        }
+
+    //! Evaluate the force and energy
+    /*!
+     * \param force_divr Holds the computed force divided by r
+     * \param pair_eng Holds the computed pair energy
+     * \param energy_shift If true, the potential is shifted to zero at the cutoff
+     *
+     * \returns True if the energy calculation occurs
+     *
+     * The calculation does not occur if the pair distance is greater than the cutoff
+     * or if the potential is scaled to zero.
+     */
+
+    DEVICE bool evalForceAndEnergy(Scalar& force_divr, Scalar& pair_eng, bool energy_shift)
+        {
+        // compute the force divided by r in force_divr
+        if (rsq < rcutsq && a != Scalar(0))
+            {
+            Scalar r = fast::sqrt(rsq);
+            Scalar rinv = 1 / r;
+            Scalar w_rinv_shifted = w / (r - sigma + w);
+            Scalar core_repuls = pow(w_rinv_shifted, q);
+            Scalar exponent = exp(q * (r - sigma - w) / w);
+            force_divr = -a * rinv
+                         * (q * core_repuls * w_rinv_shifted
+                            - q * exponent / (w * pow(1 + exponent, 2.0)));
+            pair_eng = a * (1 / (1 + exponent) - core_repuls);
+            return true;
+            }
+        else
+            return false;
+        }
+
+#ifndef __HIPCC__
+    //! Get the name of this potential
+    /*! \returns The potential name. Must be short and all lowercase, as this is the name energies
+       will be logged as via analyze.log.
+    */
+    static std::string getName()
+        {
+        return std::string("acsw");
+        }
+#endif
+
+    protected:
+    Scalar w;     //!< well width
+    Scalar sigma; //!< Hard core repulsion distance
+    Scalar a;     //!< well depth
+    Scalar q;     //!< Steepness
+    };
+
+    } // end namespace detail
+    } // end namespace azplugins
+    } // end namespace hoomd
+
+#undef DEVICE
+
+#endif // AZPLUGINS_PAIR_EVALUATOR_ACSW_H_

src/module.cc

@@ -61,6 +61,7 @@ void export_AnisoPotentialPairTwoPatchMorse(pybind11::module&);
 void export_PotentialPairColloid(pybind11::module&);
 void export_PotentialPairHertz(pybind11::module&);
 void export_PotentialPairPerturbedLennardJones(pybind11::module&);
+void export_PotentialPairACSW(pybind11::module&);
 
 // dpd
 void export_PotentialPairDPDThermoGeneralWeight(pybind11::module&);
@@ -74,6 +75,7 @@ void export_AnisoPotentialPairTwoPatchMorseGPU(pybind11::module&);
 void export_PotentialPairColloidGPU(pybind11::module&);
 void export_PotentialPairHertzGPU(pybind11::module&);
 void export_PotentialPairPerturbedLennardJonesGPU(pybind11::module&);
+void export_PotentialPairACSWGPU(pybind11::module&);
 
 // dpd
 void export_PotentialPairDPDThermoGeneralWeightGPU(pybind11::module&);
@@ -101,6 +103,7 @@ PYBIND11_MODULE(_azplugins, m)
     export_PotentialPairColloid(m);
     export_PotentialPairHertz(m);
     export_PotentialPairPerturbedLennardJones(m);
+    export_PotentialPairACSW(m);
 
     // dpd pair
     export_PotentialPairDPDThermoGeneralWeight(m);
@@ -114,6 +117,7 @@ PYBIND11_MODULE(_azplugins, m)
     export_PotentialPairColloidGPU(m);
     export_PotentialPairHertzGPU(m);
     export_PotentialPairPerturbedLennardJonesGPU(m);
+    export_PotentialPairACSWGPU(m);
 
     // dpd pair
     export_PotentialPairDPDThermoGeneralWeightGPU(m);

src/pair.py

@@ -459,3 +459,107 @@ def __init__(self, nlist, default_r_cut=None, mode='none'):
             ),
         )
         self._add_typeparam(params)
+
+
+class ACSW(pair.Pair):
+    r"""Adjusted Continuous Square Well potential.
+
+    Args:
+        nlist (hoomd.md.nlist.NeighborList): Neighbor list.
+        r_cut (float): Default cutoff radius :math:`[\mathrm{length}]`.
+        mode (str): Energy shifting/smoothing mode.
+
+    `ACSW` is a Generalized Continuous Multiple Step potential adjusted such
+    that it is a single continuous square well where all :math:`a` and :math:`w`
+    values at a distance of :math:`\sigma` will have a potential value of zero.
+
+
+    The original GCMS potential from |Munguia-Valadez et al.|_ is:
+
+    .. math::
+
+        U(r) = U_{\rm{core}}(r, \sigma_0, w_0, q_0) + \sum^{n}_{k=1}U_{\rm{step}}(r,
+        \sigma_k, w_k, q_k)
+
+    where :math:`U_{\rm{core}}` is
+
+    .. math::
+
+        U_{\rm{core}} = \left(\frac{w_0}{r - \sigma_0 + w_0}\right)^{q_0}
+
+
+    and :math:`U_{\rm{step}}` is
+
+    .. math::
+
+        U_{\rm{step}} = \left(\frac{a_k}{1 + \exp[q_k(r-\sigma_k-w_k)/w_k]}\right)
+
+    The `ACSW` potential is
+
+    .. math::
+
+        U(r) = a\left[U_{\rm{step}}(r, \sigma, w, q) -
+        U_{\rm{core}}(r, \sigma, w, q)\right]
+
+    where
+
+    .. math::
+
+        U_{\rm{core}} = \left(\frac{w}{r - \sigma + w}\right)^q
+
+    and :math:`U_{\rm{step}}` is
+
+    .. math::
+
+        U_{\rm{step}} = \left(\frac{1}{1 + \exp[q(r-\sigma-w)/w]}\right)
+
+    Here, :math:`w` is the well width, :math:`a` is the signed well depth,
+    :math:`\sigma` is the hard core repulsion diameter, and :math:`q` is the
+    steepness of the potential.
+
+    Example::
+
+        nl = nlist.Cell()
+        agcms = pair.AGCMS(default_r_cut=4.0, nlist=nl)
+        # particle types A interacting
+        agcms.params[('A', 'A')] = dict(w=1.0, sigma=2.0, a=2.0, q=16.0)
+
+    .. py:attribute:: params
+
+        The `ACSW` potential parameters. The dictionary has the following
+        keys:
+
+        * ``w`` (`float`, **required**) - Well width :math:`w`
+          :math:`[\mathrm{length}]`
+        * ``sigma`` (`float`, **required**) - Hard core repulsion diameter
+          :math:`\sigma` :math:`[\mathrm{length}]`
+        * ``a`` (`float`, **required**) - Depth of well :math:`a`
+          :math:`[\mathrm{energy}]`
+        * ``q`` (`float`, **required**) - Steepness parameter for well
+
+        Type: :class:`~hoomd.data.typeparam.TypeParameter` [`tuple`
+        [``particle_type``, ``particle_type``], `dict`]
+
+    .. py:attribute:: mode
+
+        Energy shifting/smoothing mode: ``"none"``.
+
+        Type: `str`
+
+    .. |acutei| unicode:: U+00ED .. acute i
+    .. |Munguia-Valadez et al.| replace:: Mungu\ |acutei|\ a-Valadez et al.
+    .. _Munguia-Valadez et al.: https://doi.org/10.1088/1361-648X/ac4fe8
+
+    """
+
+    _ext_module = _azplugins
+    _cpp_class_name = 'PotentialPairACSW'
+
+    def __init__(self, nlist, default_r_cut=None, default_r_on=0, mode='none'):
+        super().__init__(nlist, default_r_cut, default_r_on, mode)
+        params = TypeParameter(
+            'params',
+            'particle_types',
+            TypeParameterDict(w=float, sigma=float, a=float, q=float, len_keys=2),
+        )
+        self._add_typeparam(params)

src/pytest/test_pair.py

@@ -238,6 +238,63 @@
     ),
 ]
 
+# ACSW
+potential_tests += [
+    # test the calculation of force and potential
+    # test goes to zero outside cutoff
+    PotentialTestCase(
+        hoomd.azplugins.pair.ACSW,
+        {'w': 1.0, 'sigma': 1.0, 'a': -1.0, 'q': 16.0},
+        3.0,
+        False,
+        4.0,
+        0.0,
+        0.0,
+    ),
+    # change w to check for changing
+    # well width
+    PotentialTestCase(
+        hoomd.azplugins.pair.ACSW,
+        {'w': 1.25, 'sigma': 1.0, 'a': -1.0, 'q': 16.0},
+        3.0,
+        False,
+        1.75,
+        -0.9977990978,
+        -0.015776422214703146,
+    ),
+    # change sigma so that now the potential will
+    # be shifted to the right
+    PotentialTestCase(
+        hoomd.azplugins.pair.ACSW,
+        {'w': 1.0, 'sigma': 2.0, 'a': -1.0, 'q': 16.0},
+        3.0,
+        False,
+        2.5,
+        -0.998142211,
+        0.010875544898269139,
+    ),
+    # change well depth to increase the attractiveness
+    PotentialTestCase(
+        hoomd.azplugins.pair.ACSW,
+        {'w': 1.0, 'sigma': 1.0, 'a': -5.0, 'q': 16.0},
+        3.0,
+        False,
+        1.50,
+        -4.990711055,
+        0.0543777724491345696,
+    ),
+    # change q to increase stiffness
+    PotentialTestCase(
+        hoomd.azplugins.pair.ACSW,
+        {'w': 1.0, 'sigma': 1.0, 'a': -1.0, 'q': 50.0},
+        3.0,
+        False,
+        1.50,
+        -0.9999999984,
+        5.1583220989569564 * 10 ** (-8.0),
+    ),
+]
+
 
 @pytest.mark.parametrize(
     'potential_test', potential_tests, ids=lambda x: x.potential.__name__
@@ -284,12 +341,14 @@ def test_energy_and_force(
     # test that the energies match reference values, half goes to each particle
     energies = potential.energies
     e = potential_test.energy
+    print(energies, e)
     if sim.device.communicator.rank == 0:
         numpy.testing.assert_array_almost_equal(energies, [0.5 * e, 0.5 * e], decimal=4)
 
     # test that the forces match reference values, should be directed along x
     forces = potential.forces
     f = potential_test.force
+    print(forces, f)
     if sim.device.communicator.rank == 0:
         numpy.testing.assert_array_almost_equal(
             forces, [[-f, 0, 0], [f, 0, 0]], decimal=4