diff --git a/doc/module-azplugins-wall.rst b/doc/module-azplugins-wall.rst index 86cdab8..337aa70 100644 --- a/doc/module-azplugins-wall.rst +++ b/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: diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 64a36dd..a5529d9 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -51,6 +51,7 @@ set(_pair_evaluators # TODO: Add names of all wall evaluators set(_wall_evaluators + Colloid LJ93 ) diff --git a/src/WallEvaluatorColloid.h b/src/WallEvaluatorColloid.h index 00b3711..29cc729 100644 --- a/src/WallEvaluatorColloid.h +++ b/src/WallEvaluatorColloid.h @@ -11,23 +11,62 @@ #ifndef AZPLUGINS_WALL_EVALUATOR_COLLOID_H_ #define AZPLUGINS_WALL_EVALUATOR_COLLOID_H_ -#ifndef NVCC -#include -#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()); + auto sigma(v["sigma"].cast()); + + 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(); + } + + 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(force_divr, rsq); if (energy_shift) @@ -176,7 +185,7 @@ class WallEvaluatorColloid return false; } -#ifndef NVCC +#ifndef __HIPCC__ //! Return the name of this potential static std::string getName() { @@ -184,18 +193,15 @@ class WallEvaluatorColloid } #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_ diff --git a/src/WallEvaluatorLJ93.h b/src/WallEvaluatorLJ93.h index 662d629..3b67d65 100644 --- a/src/WallEvaluatorLJ93.h +++ b/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. */ diff --git a/src/module.cc b/src/module.cc index ec89000..1ab81c7 100644 --- a/src/module.cc +++ b/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 diff --git a/src/pytest/test_wall.py b/src/pytest/test_wall.py index 98add73..d2a6e05 100644 --- a/src/pytest/test_wall.py +++ b/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 diff --git a/src/wall.py b/src/wall.py index 71524f1..d9a8639 100644 --- a/src/wall.py +++ b/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.