Refactor wall potential, colloid potential

doc/module-azplugins-wall.rst

@@ -12,12 +12,14 @@ azplugins.wall
 .. autosummary::
     :nosignatures:
 
+    Colloid
     LJ93
 
 .. rubric:: Details
 
 .. automodule:: hoomd.azplugins.wall
     :synopsis: Wall potentials.
-    :members: LJ93
+    :members: Colloid,
+        LJ93
     :no-inherited-members:
     :show-inheritance:

src/CMakeLists.txt

@@ -51,6 +51,7 @@ set(_pair_evaluators
 
 # TODO: Add names of all wall evaluators
 set(_wall_evaluators
+    Colloid
     LJ93
     )
 

src/WallEvaluatorColloid.h

@@ -11,23 +11,62 @@
 #ifndef AZPLUGINS_WALL_EVALUATOR_COLLOID_H_
 #define AZPLUGINS_WALL_EVALUATOR_COLLOID_H_
 
-#ifndef NVCC
-#include <string>
-#endif
-
-#include "hoomd/HOOMDMath.h"
+#include "PairEvaluator.h"
 
-#ifdef NVCC
+#ifdef __HIPCC__
 #define DEVICE __device__
 #else
 #define DEVICE
 #endif
 
+namespace hoomd
+    {
 namespace azplugins
     {
 namespace detail
     {
-//! Evaluates the Lennard-Jones colloid wall force
+//! Define the parameter type used by this wall potential evaluator
+struct WallParametersColloid : public PairParameters
+    {
+#ifndef __HIPCC__
+    WallParametersColloid() : c_1(0), c_2(0), a(0) { }
+
+    WallParametersColloid(pybind11::dict v, bool managed = false)
+        {
+        auto A(v["A"].cast<Scalar>());
+        auto sigma(v["sigma"].cast<Scalar>());
+
+        const Scalar sigma_3 = sigma * sigma * sigma;
+        c_1 = A * sigma_3 * sigma_3 / Scalar(7560);
+        c_2 = A / Scalar(6);
+        a = v["a"].cast<Scalar>();
+        }
+
+    pybind11::dict asDict()
+        {
+        const Scalar A = Scalar(6) * c_2;
+        const Scalar sigma_6 = Scalar(7560) * c_1 / A;
+
+        pybind11::dict v;
+        v["A"] = A;
+        v["a"] = a;
+        v["sigma"] = pow(sigma_6, 1. / 6.);
+        return v;
+        }
+#endif // __HIPCC__
+
+    Scalar c_1; //!< First coefficient
+    Scalar c_2; //!< Second coefficient
+    Scalar a;   //!< Particle radius
+    }
+
+#if HOOMD_LONGREAL_SIZE == 32
+    __attribute__((aligned(8)));
+#else
+    __attribute__((aligned(16)));
+#endif
+
+//! Class for evaluating the Lennard-Jones colloid wall force
 /*!
  * The Lennard-Jones colloid wall potential is derived from integrating the standard Lennard-Jones
  * potential between a spherical particle of radius \f$ a \f$ and a half plane, and it takes the
@@ -49,11 +88,11 @@ namespace detail
  * where \f$ \rho_{\rm w} \f$ and \f$ \rho_{\rm c} \f$ are the wall and colloid densities and \f$
  * \epsilon \f$ is the Lennard-Jones interaction energy.
  */
-class WallEvaluatorColloid
+class WallEvaluatorColloid : public PairEvaluator
     {
     public:
     //! Define the parameter type used by this wall potential evaluator
-    typedef Scalar2 param_type;
+    typedef WallParametersColloid param_type;
 
     //! Constructor
     /*!
@@ -64,43 +103,13 @@ class WallEvaluatorColloid
      * The functor initializes its members from \a _params.
      */
     DEVICE WallEvaluatorColloid(Scalar _rsq, Scalar _rcutsq, const param_type& _params)
-        : rsq(_rsq), rcutsq(_rcutsq), A(_params.x), B(_params.y)
+        : PairEvaluator(_rsq, _rcutsq)
         {
+        c_1 = _params.c_1;
+        c_2 = _params.c_2;
+        a = _params.a;
         }
 
-    //! Colloid diameter is needed
-    DEVICE static bool needsDiameter()
-        {
-        return true;
-        }
-    //! Accept the optional diameter values
-    /*!
-     * \param di Diameter of particle
-     * \param dj Dummy diameter
-     *
-     * \note The way HOOMD computes wall forces by recycling evaluators requires that we give
-     *       a second diameter, even though this is meaningless for the potential.
-     */
-    DEVICE void setDiameter(Scalar di, Scalar dj)
-        {
-        a = Scalar(0.5) * di;
-        }
-
-    //! Colloid wall potential doesn't use charge
-    DEVICE static bool needsCharge()
-        {
-        return false;
-        }
-    //! Accept the optional charge values
-    /*!
-     * \param qi Charge of particle
-     * \param qj Dummy charge
-     *
-     * \note The way HOOMD computes wall forces by recycling evaluators requires that we give
-     *       a second charge, even though this is meaningless for the potential.
-     */
-    DEVICE void setCharge(Scalar qi, Scalar qj) { }
-
     //! Computes the colloid-wall interaction potential
     /*!
      * \param force_divr Force divided by r
@@ -135,18 +144,18 @@ class WallEvaluatorColloid
         if (force)
             {
             Scalar arinv8 = Scalar(8.0) * arinv;
-            force_divr = Scalar(6.0) * A
+            force_divr = Scalar(6.0) * c_1
                          * ((arinv8 - Scalar(1.0)) * r_minus_a_inv2 * r_minus_a_inv6
                             + (arinv8 + Scalar(1.0)) * r_plus_a_inv2 * r_plus_a_inv6);
-            force_divr -= B * (Scalar(4.0) * a * a * arinv * r2_minus_a2_inv * r2_minus_a2_inv);
+            force_divr -= c_2 * (Scalar(4.0) * a * a * arinv * r2_minus_a2_inv * r2_minus_a2_inv);
             }
 
         // energy
         Scalar a7 = Scalar(7.0) * a;
-        Scalar energy = A
+        Scalar energy = c_1
                         * ((a7 - r) * r_minus_a_inv * r_minus_a_inv6
                            + (a7 + r) * r_plus_a_inv * r_plus_a_inv6);
-        energy -= B * (Scalar(2.0) * a * r * r2_minus_a2_inv + log(r_plus_a_inv / r_minus_a_inv));
+        energy -= c_2 * (Scalar(2.0) * a * r * r2_minus_a2_inv + log(r_plus_a_inv / r_minus_a_inv));
         return energy;
         }
 
@@ -163,7 +172,7 @@ class WallEvaluatorColloid
      */
     DEVICE bool evalForceAndEnergy(Scalar& force_divr, Scalar& energy, bool energy_shift)
         {
-        if (rsq < rcutsq && A != 0 && a > 0)
+        if (rsq < rcutsq && c_1 != 0 && a > 0)
             {
             energy = computePotential<true>(force_divr, rsq);
             if (energy_shift)
@@ -176,26 +185,23 @@ class WallEvaluatorColloid
             return false;
         }
 
-#ifndef NVCC
+#ifndef __HIPCC__
     //! Return the name of this potential
     static std::string getName()
         {
         return std::string("colloid");
         }
 #endif
 
-    protected:
-    Scalar rsq;    //!< Stored rsq from the constructor
-    Scalar rcutsq; //!< Stored rcutsq from the constructor
-
-    Scalar A; //!< Prefactor for first term (includes Hammaker constant)
-    Scalar B; //!< Prefactor for second term (includes Hammaker constant)
-
-    Scalar a; //!< The particle radius
+    private:
+    Scalar c_1; //!< Prefactor for first term (includes Hamaker constant)
+    Scalar c_2; //!< Prefactor for second term (includes Hamaker constant)
+    Scalar a;   //!< The particle radius
     };
 
     } // end namespace detail
-    } // namespace azplugins
+    } // end namespace azplugins
+    } // end namespace hoomd
 
 #undef DEVICE
 #endif // AZPLUGINS_WALL_EVALUATOR_COLLOID_H_

src/WallEvaluatorLJ93.h

@@ -24,7 +24,7 @@ namespace azplugins
     {
 namespace detail
     {
-//! Define the paramter type used by this wall potential evaluator
+//! Define the parameter type used by this wall potential evaluator
 struct WallParametersLJ93 : public PairParameters
     {
 #ifndef __HIPCC__
@@ -85,9 +85,9 @@ class WallEvaluatorLJ93 : public PairEvaluator
 
     //! Constructor
     /*!
-     * \param _rsq Sqaured distance between particles
+     * \param _rsq Squared distance between particles
      * \param _rcutsq Cutoff radius squared
-     * \param _params Wall potential paramters, given by typedef above
+     * \param _params Wall potential parameters, given by typedef above
      *
      * The functor initializes its members from \a _params.
      */

src/module.cc

@@ -78,6 +78,7 @@ void export_PotentialPairPerturbedLennardJones(pybind11::module&);
 void export_PotentialPairDPDThermoGeneralWeight(pybind11::module&);
 
 // wall
+void export_PotentialExternalWallColloid(pybind11::module&);
 void export_PotentialExternalWallLJ93(pybind11::module&);
 
 #ifdef ENABLE_HIP
@@ -106,6 +107,7 @@ void export_PotentialPairPerturbedLennardJonesGPU(pybind11::module&);
 void export_PotentialPairDPDThermoGeneralWeightGPU(pybind11::module&);
 
 // wall
+void export_PotentialExternalWallColloidGPU(pybind11::module&);
 void export_PotentialExternalWallLJ93GPU(pybind11::module&);
 
 #endif // ENABLE_HIP
@@ -148,6 +150,7 @@ PYBIND11_MODULE(_azplugins, m)
     export_PotentialPairDPDThermoGeneralWeight(m);
 
     // wall
+    export_PotentialExternalWallColloid(m);
     export_PotentialExternalWallLJ93(m);
 
 #ifdef ENABLE_HIP
@@ -176,6 +179,7 @@ PYBIND11_MODULE(_azplugins, m)
     export_PotentialPairDPDThermoGeneralWeightGPU(m);
 
     // wall
+    export_PotentialExternalWallColloidGPU(m);
     export_PotentialExternalWallLJ93GPU(m);
 
 #endif // ENABLE_HIP

src/pytest/test_wall.py

@@ -19,6 +19,25 @@
 
 potential_tests = []
 
+# Colloid
+potential_tests += [
+    # test the calculation of force and potential
+    PotentialTestCase(
+        hoomd.azplugins.wall.Colloid,
+        {"A": 1.0, "a": 1.0, "sigma": 1.0, "r_cut": 3.0},
+        (0, 0, 1.5),
+        -0.0292,
+        0.8940,
+    ),
+    PotentialTestCase(
+        hoomd.azplugins.wall.Colloid,
+        {"A": 1.0, "a": 1.0, "sigma": 1.0, "r_cut": 3.0},
+        (0, 0, 3.5),
+        0,
+        0,
+    ),
+]
+
 # Lennard Jones 9-3
 potential_tests += [
     # test the calculation of force and potential

src/wall.py

@@ -10,6 +10,77 @@
 from hoomd.md.external import wall
 
 
+class Colloid(wall.WallPotential):
+    r"""Colloid wall potential.
+
+    Args:
+        walls (`list` [`hoomd.wall.WallGeometry`]) : A list of wall definitions.
+
+    `Colloid` is Lennard-Jones colloid wall potential, which is derived from integrating
+    the standard Lennard-Jones potential between a spherical particle of radius
+    :math:`a` and a half plane:
+
+    .. math::
+
+        U(r) = \frac{A \sigma^6}{7560} \left[\frac{7a-r}{(r-a)^7} +
+            \frac{7a+r}{(r+a)^7} \right] - \frac{A}{6} \left[\frac{2ar}{r^2-a^2}
+            + \ln\left(\frac{z-a}{z+a}\right) \right]
+
+    where:
+
+    * :math:`\sigma` - diameter of Lennard-Jones particles in the sphere and the wall
+    * :math:`A` - Hamaker constant, related to the Lennard-Jones parameters as
+        :math:`A = 4 \pi^2 \varepsilon \sigma^6 \rho_{\rm w} \rho_{\rm c}`, where
+        :math:`\rho_{\rm w}` is the number density of particles in the wall,
+        :math:`\rho_{\rm c}` is the number density of particles in the sphere, and
+        :math:`\varepsilon` is the Lennard-Jones energy parameter.
+
+    Example::
+
+        walls = [hoomd.wall.Plane(origin=(0, 0, 0), normal=(0, 1, 0))]
+        colloid = azplugins.wall.Colloid(walls=walls)
+        colloid.params["A"] = dict(
+            A=4 * numpy.pi ** 2,  # rho = 1, epsilon = 1, sigma = 1
+            a=1.0
+            sigma=1.0,
+            r_cut=3.0,
+        )
+
+    .. py:attribute:: params
+
+        The `Colloid` potential parameters. The dictionary has the following
+        keys:
+
+        * ``A`` (`float`, **required**) - Hamaker constant
+          :math:`A` :math:`[\mathrm{energy}]`
+        * ``a`` (`float`, **required**) - Radius of particle
+          :math:`a` :math:`[\mathrm{length}]`
+        * ``sigma`` (`float`, **required**) - Lennard-Jones particle size
+          :math:`\sigma` :math:`[\mathrm{length}]`
+        * ``r_cut`` (`float`, **required**) - The cut off distance for
+          the wall potential :math:`[\mathrm{length}]`
+        * ``r_extrap`` (`float`, **optional**) - The distance to
+          extrapolate the potential, defaults to 0.
+          :math:`[\mathrm{length}]`.
+
+        Type: :class:`~hoomd.data.typeparam.TypeParameter` [``particle_type``], `dict`]
+    """
+
+    _ext_module = _azplugins
+    _cpp_class_name = "WallsPotentialColloid"
+
+    def __init__(self, walls):
+        super().__init__(walls)
+        params = TypeParameter(
+            "params",
+            "particle_types",
+            TypeParameterDict(
+                A=float, a=float, sigma=float, r_cut=float, r_extrap=0.0, len_keys=1
+            ),
+        )
+        self._add_typeparam(params)
+
+
 class LJ93(wall.WallPotential):
     r"""Lennard-Jones 9-3 wall potential.