diff --git a/doc/module-azplugins-pair.rst b/doc/module-azplugins-pair.rst index 56ecbf27..3fd43a58 100644 --- a/doc/module-azplugins-pair.rst +++ b/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, diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 030a8d17..68e51d2d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -32,6 +32,7 @@ set(_pair_evaluators Colloid Hertz PerturbedLennardJones + ACSW ) # TODO: Add names of all dpd evaluators diff --git a/src/PairEvaluatorACSW.h b/src/PairEvaluatorACSW.h new file mode 100644 index 00000000..e6bcbda7 --- /dev/null +++ b/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 + +#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(); + sigma = v["sigma"].cast(); + a = v["a"].cast(); + q = v["q"].cast(); + } + + 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_ diff --git a/src/module.cc b/src/module.cc index f3e183d5..0e4bf5a7 100644 --- a/src/module.cc +++ b/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); diff --git a/src/pair.py b/src/pair.py index a0a49539..4fb1cb3a 100644 --- a/src/pair.py +++ b/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) diff --git a/src/pytest/test_pair.py b/src/pytest/test_pair.py index 006d14ed..a172e8b0 100644 --- a/src/pytest/test_pair.py +++ b/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