From c3f96727ca810ba4c80c8b638c81089a90409fef Mon Sep 17 00:00:00 2001 From: Jinny Cha Date: Tue, 2 Dec 2025 10:11:21 -0600 Subject: [PATCH 1/4] Refactor wall potential --- src/CMakeLists.txt | 14 +++ src/WallEvaluator.h | 134 ++++++++++++++++++++ src/WallEvaluatorLJ93.h | 110 ++++++++-------- src/WallPotentials.cu | 27 ---- src/WallPotentials.h | 91 -------------- src/__init__.py | 2 +- src/export_PotentialExternalWall.cc.inc | 33 +++++ src/module.cc | 6 + src/pytest/CMakeLists.txt | 1 + src/pytest/test_wall.py | 161 ++++++++++++++++++++++++ src/wall.py | 72 +++++++++++ 11 files changed, 478 insertions(+), 173 deletions(-) create mode 100644 src/WallEvaluator.h delete mode 100644 src/WallPotentials.cu delete mode 100644 src/WallPotentials.h create mode 100644 src/export_PotentialExternalWall.cc.inc create mode 100644 src/pytest/test_wall.py create mode 100644 src/wall.py diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ebde4dd5..96e7f145 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -32,6 +32,7 @@ set(python_files external.py flow.py pair.py + wall.py ) # TODO: Add names of all bond evaluators @@ -48,6 +49,11 @@ set(_pair_evaluators PerturbedLennardJones ) +# TODO: Add names of all wall evaluators +set(_wall_evaluators + LJ93 + ) + # TODO: Add names of all dpd evaluators set(_dpd_evaluators GeneralWeight @@ -96,6 +102,14 @@ foreach(_evaluator ${_pair_evaluators}) endif() endforeach() +# process wall potentials +foreach(_evaluator ${_wall_evaluators}) + configure_file(export_PotentialExternalWall.cc.inc + export_PotentialExternalWall${_evaluator}.cc + @ONLY) + set(_${COMPONENT_NAME}_sources ${_${COMPONENT_NAME}_sources} export_PotentialExternalWall${_evaluator}.cc) +endforeach() + # process DPD potentials foreach(_evaluator ${_dpd_evaluators}) configure_file(export_PotentialPairDPDThermo.cc.inc diff --git a/src/WallEvaluator.h b/src/WallEvaluator.h new file mode 100644 index 00000000..170295aa --- /dev/null +++ b/src/WallEvaluator.h @@ -0,0 +1,134 @@ +// Copyright (c) 2018-2020, Michael P. Howard +// Copyright (c) 2021-2025, Auburn University +// Part of azplugins, released under the BSD 3-Clause License. + +#ifndef AZPLUGINS_WALL_EVALUATOR_H_ +#define AZPLUGINS_WALL_EVALUATOR_H_ + +#ifndef __HIPCC__ +#include +#include +#endif // __HIPCC__ + +#include "hoomd/HOOMDMath.h" + +#ifdef __HIPCC__ +#define DEVICE __device__ +#define HOSTDEVICE __host__ __device__ +#else +#define DEVICE +#define HOSTDEVICE +#endif // __HIPCC__ + +namespace hoomd + { +namespace azplugins + { +namespace detail + { +struct WallParameters + { + +#ifndef __HIPCC__ + WallParameters() { } + + WallParameters(pybind11::dict v, bool managed = false) { } + + pybind11::dict asDict() + { + return pybind11::dict(); + } +#endif //__HIPCC__ + + DEVICE void load_shared(char*& ptr, unsigned int& available_bytes) { } + + HOSTDEVICE void allocate_shared(char*& ptr, unsigned int& available_bytes) const { } + +#ifdef ENABLE_HIP + void set_memory_hint() const { } +#endif + }; + +class WallEvaluator + { + public: + typedef WallParameters param_type; + + DEVICE WallEvaluator(const Scalar _rsq, const Scalar _rcutsq): rsq(_rsq), rcutsq(_rcutsq) {} + + //! Base potential does not need charge + DEVICE static bool needsCharge() + { + return false; + } + + //! Accept the optional charge values + /*! + * \param qi Charge of particle i + * \param qj Charge of particle j + */ + DEVICE void setCharge(Scalar qi, Scalar qj) { } + + //! Evaluate the force and energy + /*! + * \param force_divr Holds the computed force divided by r + * \param energy 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& energy, bool energy_shift) + { + return false; + } + + //! Evaluate the long-ranged correction to the pressure. + /*! + * \returns Default value of 0 (no correction). + */ + DEVICE Scalar evalPressureLRCIntegral() + { + return Scalar(0.0); + } + + //! Evaluate the long-ranged correction to the energy. + /*! + * \returns Default value of 0 (no correction). + */ + DEVICE Scalar evalEnergyLRCIntegral() + { + return Scalar(0.0); + } + +#ifndef __HIPCC__ + //! Get the name of this potential + /*! + * This method must be overridden by deriving classes. + */ + static std::string getName() + { + throw std::runtime_error("Name not defined for this wall potential."); + } + + std::string getShapeSpec() const + { + throw std::runtime_error("Shape definition not supported for this wall potential."); + } +#endif // __HIPCC__ + + protected: + Scalar rsq; //!< Squared distance between particles + Scalar rcutsq; //!< Squared cutoff distance + }; + + } // end namespace detail + } // end namespace azplugins + } // end namespace hoomd + +#undef DEVICE +#undef HOSTDEVICE + +#endif // AZPLUGINS_WALL_EVALUATOR_H_ \ No newline at end of file diff --git a/src/WallEvaluatorLJ93.h b/src/WallEvaluatorLJ93.h index 74fa224b..caad6d25 100644 --- a/src/WallEvaluatorLJ93.h +++ b/src/WallEvaluatorLJ93.h @@ -10,23 +10,54 @@ #ifndef AZPLUGINS_WALL_EVALUATOR_LJ_93_H_ #define AZPLUGINS_WALL_EVALUATOR_LJ_93_H_ -#ifndef NVCC -#include -#endif - -#include "hoomd/HOOMDMath.h" +#include "WallEvaluator.h" -#ifdef NVCC +#ifdef __HIPCC__ #define DEVICE __device__ #else #define DEVICE #endif +namespace hoomd + { namespace azplugins { namespace detail { -//! Evaluates the Lennard-Jones 9-3 wall force +//! Define the paramter type used by this wall potential evaluator +struct WallParametersLJ93 : public WallParameters + { +#ifndef __HIPCC__ + WallParametersLJ93(): sigma_3(0), epsilon(0) { } + + WallParametersLJ93(pybind11::dict v, bool managed = false) + { + auto sigma(v["sigma"].cast()); + epsilon = v["epsilon"].cast(); + + sigma_3 = sigma * sigma * sigma; + } + + pybind11::dict asDict() + { + pybind11::dict v; + v["sigma"] = std::cbrt(sigma_3); + v["epsilon"] = epsilon; + return v; + } +#endif // __HIPCC__ + + Scalar sigma_3; + Scalar epsilon; + } + +#if HOOMD_LONGREAL_SIZE == 32 + __attribute__((aligned(16))); +#else + __attribute__((aligned(32))); +#endif + +//! Class for evaluatring the Lennard-Jones 9-3 wall force /*! * WallEvaluatorLJ93 computes the Lennard-Jones 9-3 wall potential, which is derived from * integrating the standard Lennard-Jones potential between a point particle and a half plane: @@ -47,55 +78,26 @@ namespace detail * \f[ F(r)/r = \frac{\varepsilon}{r^2} \left ( \frac{6}{5}\left(\frac{\sigma}{r}\right)^9 - 3 * \left(\frac{\sigma}{r}\right)^3 \right) \f] */ -class WallEvaluatorLJ93 +class WallEvaluatorLJ93 : public WallEvaluator { public: - //! Define the parameter type used by this wall potential evaluator - typedef Scalar2 param_type; + typedef WallParametersLJ93 param_type; //! Constructor /*! - * \param _rsq Squared distance between particles + * \param _rsq Sqaured distance between particles * \param _rcutsq Cutoff radius squared - * \param _params Pair potential parameters, given by typedef above - * + * \param _params Wall potential paramters, given by typedef above + * * The functor initializes its members from \a _params. - */ + */ DEVICE WallEvaluatorLJ93(Scalar _rsq, Scalar _rcutsq, const param_type& _params) - : rsq(_rsq), rcutsq(_rcutsq), lj1(_params.x), lj2(_params.y) - { - } - - //! LJ 9-3 doesn't use diameter - DEVICE static bool needsDiameter() - { - return false; - } - //! 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) { } - - //! LJ 9-3 doesn't use charge - DEVICE static bool needsCharge() + : WallEvaluator(_rsq, _rcutsq) { - return false; + lj1 = (Scalar(2.0) / Scalar(15.0)) * _params.epsilon * _params.sigma_3 * _params.sigma_3 * _params.sigma_3; + lj2 = _params.epsilon * _params.sigma_3; } - //! 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) { } - + //! Evaluate the force and energy /*! * \param force_divr Holds the computed force divided by r @@ -134,8 +136,8 @@ class WallEvaluatorLJ93 else return false; } - -#ifndef NVCC + +#ifndef __HIPCC__ //! Return the name of this potential static std::string getName() { @@ -143,14 +145,14 @@ class WallEvaluatorLJ93 } #endif - protected: - Scalar rsq; //!< Stored rsq from the constructor - Scalar rcutsq; //!< Stored rcutsq from the constructor - Scalar lj1; //!< lj1 parameter extracted from the params passed to the constructor - Scalar lj2; //!< lj2 parameter extracted from the params passed to the constructor + private: + Scalar lj1; + Scalar lj2; }; + } // end namespace detail } // end namespace azplugins + } // end namespace hoomd #undef DEVICE -#endif // AZPLUGINS_WALL_EVALUATOR_LJ_93_H_ +#endif // AZPLUGINS_WALL_EVALUATOR_LJ_93_H_ \ No newline at end of file diff --git a/src/WallPotentials.cu b/src/WallPotentials.cu deleted file mode 100644 index b58296f1..00000000 --- a/src/WallPotentials.cu +++ /dev/null @@ -1,27 +0,0 @@ -// Copyright (c) 2018-2020, Michael P. Howard -// Copyright (c) 2021-2025, Auburn University -// Part of azplugins, released under the BSD 3-Clause License. - -/*! - * \file WallPotentials.cu - * \brief Defines the driver functions for computing wall forces on the GPU - * - * The driver function must have an explicit template instantiation for each - * evaluator. The template must be in the global namespace that is used by hoomd. - */ - -#include "WallPotentials.h" -#include "hoomd/md/EvaluatorWalls.h" -#include "hoomd/md/PotentialExternalGPU.cuh" - -//! Evaluator for colloid (integrated Lennard-Jones) wall potential -template cudaError_t gpu_cpef>( - const external_potential_args_t& external_potential_args, - const typename EvaluatorWalls::param_type* d_params, - const typename EvaluatorWalls::field_type* d_field); - -//! Evaluator for Lennard-Jones 9-3 wall potential -template cudaError_t gpu_cpef>( - const external_potential_args_t& external_potential_args, - const typename EvaluatorWalls::param_type* d_params, - const typename EvaluatorWalls::field_type* d_field); diff --git a/src/WallPotentials.h b/src/WallPotentials.h deleted file mode 100644 index 456741da..00000000 --- a/src/WallPotentials.h +++ /dev/null @@ -1,91 +0,0 @@ -// Copyright (c) 2018-2020, Michael P. Howard -// Copyright (c) 2021-2025, Auburn University -// Part of azplugins, released under the BSD 3-Clause License. - -/*! - * \file WallPotentials.h - * \brief Convenience inclusion of all wall potential evaluators. - * - * In HOOMD-blue, wall potentials are templated on a base class ForceCompute called - * PotentialExternal (via a templated EvaluatorWall functor). This avoids code duplication of basic - * routines. - * - * To add a new wall potential, take the following steps: - * 1. Create an evaluator functor for your potential, for example WallEvaluatorMyGreatPotential.h. - * This file should be added to the list of includes below. You can follow one of the other - * evaluator functors as an example for the details. - * - * 2. Explicitly instantiate a template for a CUDA driver for your potential in WallPotentials.cu. - * - * 3. Expose the wall potential on the python level in module.cc with export_wall_potential and add - * the mirror python object to wall.py. - * - * 4. Write a unit test for the potential in test-py. Two types of tests should be conducted: one - * that checks that all methods work on the python object, and one that validates the force and - * energy for the particle at a fixed distance from the wall. - */ - -#ifndef AZPLUGINS_WALL_POTENTIALS_H_ -#define AZPLUGINS_WALL_POTENTIALS_H_ - -// All wall potential evaluators should be included here -#include "WallEvaluatorColloid.h" -#include "WallEvaluatorLJ93.h" - -/* - * The code below handles python exports using a templated function, and so should - * not be compiled in NVCC. - */ -#ifndef NVCC -#include "hoomd/extern/pybind/include/pybind11/pybind11.h" - -#include "hoomd/md/EvaluatorWalls.h" -#include "hoomd/md/PotentialExternal.h" -#ifdef ENABLE_CUDA -#include "hoomd/md/PotentialExternalGPU.h" -#endif - -namespace azplugins - { -namespace detail - { -//! Convenience function to export wall potential to python -/*! - * \param m Python module - * \param name Name of CPU class - * \tparam evaluator Wall potential evaluator functor - * - * A CPU object called \a name is exported to python, along with an - * object for the wall parameters (name_params) and convenience method - * for creating the parameters (make_name_params()). If CUDA is enabled, - * a corresponding GPU class called nameGPU is exported. - */ -template void export_wall_potential(pybind11::module& m, const std::string& name) - { - namespace py = pybind11; - typedef ::EvaluatorWalls wall_evaluator; - typedef ::PotentialExternal wall_potential_cpu; - export_PotentialExternal(m, name); - -#ifdef ENABLE_CUDA - typedef ::PotentialExternalGPU wall_potential_gpu; - export_PotentialExternalGPU(m, name + "GPU"); -#endif // ENABLE_CUDA - - py::class_>( - m, - (wall_evaluator::getName() + "_params").c_str()) - .def(py::init<>()) - .def_readwrite("params", &wall_evaluator::param_type::params) - .def_readwrite("rextrap", &wall_evaluator::param_type::rextrap) - .def_readwrite("rcutsq", &wall_evaluator::param_type::rcutsq); - m.def(std::string("make_" + wall_evaluator::getName() + "_params").c_str(), - &make_wall_params); - } - - } // end namespace detail - } // end namespace azplugins -#endif // NVCC - -#endif // AZPLUGINS_WALL_POTENTIALS_H_ diff --git a/src/__init__.py b/src/__init__.py index a7ca8a70..3e71abf0 100644 --- a/src/__init__.py +++ b/src/__init__.py @@ -4,6 +4,6 @@ """azplugins.""" -from hoomd.azplugins import bond, compute, external, flow, pair +from hoomd.azplugins import bond, compute, external, flow, pair, wall __version__ = "1.2.0" diff --git a/src/export_PotentialExternalWall.cc.inc b/src/export_PotentialExternalWall.cc.inc new file mode 100644 index 00000000..be769757 --- /dev/null +++ b/src/export_PotentialExternalWall.cc.inc @@ -0,0 +1,33 @@ +// Copyright (c) 2018-2020, Michael P. Howard +// Copyright (c) 2021-2025, Auburn University +// Part of azplugins, released under the BSD 3-Clause License. + +// Adapted from hoomd/md/export_PotentialExternalWall.cc.inc of HOOMD-blue. +// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +// See CMakeLists.txt for the source of these variables to be processed by CMake's +// configure_file(). + +// clang-format off +#include "hoomd/md/PotentialExternal.h" +#include "hoomd/md/EvaluatorWalls.h" +#include "WallEvaluator@_evaluator@.h" + +#define EVALUATOR_CLASS hoomd::md::EvaluatorWalls +#define EXPORT_FUNCTION export_PotentialExternalWall@_evaluator@ +// clang-format on + +namespace hoomd + { +namespace azplugins + { +namespace detail + { +void EXPORT_FUNCTION(pybind11::module& m) + { + hoomd::md::detail::export_PotentialExternal(m, "WallsPotential@_evaluator@"); + } + } // end namespace detail + } // end namespace azplugins + } // end namespace hoomd \ No newline at end of file diff --git a/src/module.cc b/src/module.cc index 90e25494..54a264bf 100644 --- a/src/module.cc +++ b/src/module.cc @@ -77,6 +77,9 @@ void export_PotentialPairPerturbedLennardJones(pybind11::module&); // dpd void export_PotentialPairDPDThermoGeneralWeight(pybind11::module&); +// wall +void export_PotentialExternalWallLJ93(pybind11::module&); + #ifdef ENABLE_HIP // bond void export_ImagePotentialBondHarmonicGPU(pybind11::module&); @@ -141,6 +144,9 @@ PYBIND11_MODULE(_azplugins, m) // dpd pair export_PotentialPairDPDThermoGeneralWeight(m); + // wall + export_PotentialExternalWallLJ93(m); + #ifdef ENABLE_HIP // bond export_ImagePotentialBondHarmonicGPU(m); diff --git a/src/pytest/CMakeLists.txt b/src/pytest/CMakeLists.txt index 1c6f60ba..157b326b 100644 --- a/src/pytest/CMakeLists.txt +++ b/src/pytest/CMakeLists.txt @@ -8,6 +8,7 @@ set(test_files test_pair.py test_pair_aniso.py test_pair_dpd.py + test_wall.py ) # Copy tests to the install directory. diff --git a/src/pytest/test_wall.py b/src/pytest/test_wall.py new file mode 100644 index 00000000..7cc8b9f6 --- /dev/null +++ b/src/pytest/test_wall.py @@ -0,0 +1,161 @@ +# Copyright (c) 2018-2020, Michael P. Howard +# Copyright (c) 2021-2025, Auburn University +# Part of azplugins, released under the BSD 3-Clause License. + +"""Wall potential unit tests.""" + +import collections + +import hoomd +import hoomd.azplugins +import numpy + +import pytest + +PotentialTestCase = collections.namedtuple( + "PotentialTestCase", + ["potential", "params", "wall", "position", "energy", "force"], +) + +potential_tests = [] + +# Lennard Jones 9-3 +potential_tests += [ + # test the calculation of force and potential for plane + PotentialTestCase( + hoomd.azplugins.wall.LJ93, + {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0, "r_extrap": 0.0}, + hoomd.wall.Plane(origin=(0, 0, 0), normal=(0, 0, 1)), + (0, 0, 1.5), + -0.2558, + -0.5718, + ), + + PotentialTestCase( + hoomd.azplugins.wall.LJ93, + {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0, "r_extrap": 0.0}, + hoomd.wall.Plane(origin=(0, 0, 0), normal=(0, 0, 1)), + (0, 0, 3.5), + 0, + 0, + ), + + # test the calculation of force and potential for sphere + PotentialTestCase( + hoomd.azplugins.wall.LJ93, + {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0, "r_extrap": 0.0}, + hoomd.wall.Sphere(radius=5, origin=(0, 0, 0), inside=True), + (0, 0, 3.5), + -0.2558, + 0.5718, + ), + + PotentialTestCase( + hoomd.azplugins.wall.LJ93, + {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0, "r_extrap": 0.0}, + hoomd.wall.Sphere(radius=5, origin=(0, 0, 0), inside=True), + (0, 0, 2.0), + 0, + 0, + ), + PotentialTestCase( + hoomd.azplugins.wall.LJ93, + {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0, "r_extrap": 0.0}, + hoomd.wall.Sphere(radius=5, origin=(0, 0, 0), inside=True), + (0, 0, 5.5), + 0, + 0, + ), + + # test the calculation of force and potential for cylinder + PotentialTestCase( + hoomd.azplugins.wall.LJ93, + {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0, "r_extrap": 0.0}, + hoomd.wall.Cylinder(radius=5, origin=(0, 0, 0), axis=(1, 0, 0), inside=True), + (0, 0, 3.5), + -0.2558, + 0.5718, + ), + + PotentialTestCase( + hoomd.azplugins.wall.LJ93, + {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0, "r_extrap": 0.0}, + hoomd.wall.Cylinder(radius=5, origin=(0, 0, 0), axis=(1, 0, 0), inside=True), + (0, 0, 5.5), + 0, + 0, + ), + + PotentialTestCase( + hoomd.azplugins.wall.LJ93, + {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0, "r_extrap": 0.0}, + hoomd.wall.Sphere(radius=5, origin=(0, 0, 0), inside=True), + (0, 0, 2.0), + 0, + 0, + ), +] + +@pytest.mark.parametrize( +"potential_test", potential_tests, ids=lambda x: x.potential.__name__ +) +def test_energy_and_force( +simulation_factory, one_particle_snapshot_factory, potential_test +): + """Test energy and force evaluation.""" + # make one particle test configuration + sim = simulation_factory( + one_particle_snapshot_factory(position=potential_test.position, L=20) + ) + + # setup dummy NVE integrator + integrator = hoomd.md.Integrator(dt=0.001) + nve = hoomd.md.methods.ConstantVolume(hoomd.filter.All()) + integrator.methods = [nve] + + # setup wall potential + wall = potential_test.wall + potential = potential_test.potential( + walls=[wall] + ) + potential.params["A"] = potential_test.params + integrator.forces = [potential] + + # calculate energies and forces + sim.operations.integrator = integrator + sim.run(0) + + # test that parameters are still correct after attach runs + assert potential.params["A"] == potential_test.params + + # test that the energies match reference values, half goes to each particle + energies = potential.energies + e = potential_test.energy + if sim.device.communicator.rank == 0: + numpy.testing.assert_array_almost_equal(energies, e, decimal=4) + + # test that the forces match reference values, should be directed along x + forces = potential.forces + f = potential_test.force + if sim.device.communicator.rank == 0: + numpy.testing.assert_array_almost_equal( + forces, [[0, 0, f]], decimal=4 + ) + + # test if particle stays on correct side of the wall + for _ in range(10): + sim.run(10) + with sim.state.cpu_local_snapshot as snap: + pos = snap.particles.position + + if wall == hoomd.wall.Plane: + assert numpy.all(snap.particles.position[:, 2] > 0) + + elif wall == hoomd.wall.Sphere: + assert numpy.all(numpy.linalg.norm(snap.particles.position, axis=1) < wall.radius) + + elif wall == hoomd.wall.Cylinder: + axis = numpy.array(wall.axis, dtype=float) + axis /= numpy.linalg.norm(axis) + radial = pos - numpy.dot(pos, axis) * axis + assert numpy.linalg.norm(radial) < wall.radius \ No newline at end of file diff --git a/src/wall.py b/src/wall.py new file mode 100644 index 00000000..1015bbb0 --- /dev/null +++ b/src/wall.py @@ -0,0 +1,72 @@ +# Copyright (c) 2018-2020, Michael P. Howard +# Copyright (c) 2021-2025, Auburn University +# Part of azplugins, released under the BSD 3-Clause License. + +"""Wall potentials.""" + +from hoomd.azplugins import _azplugins +from hoomd.data.parameterdicts import ParameterDict, TypeParameterDict +from hoomd.data.typeparam import TypeParameter +from hoomd.md.external import wall +from hoomd.variant import Variant + + +class LJ93(wall.WallPotential): + r"""Lennard-Jones 9-3 wall potential. + + Args: + walls (`list` [`hoomd.wall.WallGeometry`]) : A list of wall definitions. + + `LJ93` is Lennard-Jones 9-3 wall potential, which is derived from integrating + the standard Lennard-Jones potential between a point particle and a half plane: + + .. math:: + + U(r) = \varepsilon \left(\frac{2}{15} \left(\frac{\sigma}{r} \right)^9 + - \left(\frac{\sigma}{r}\right)^3 \right)] + + where: + + * :math:`\sigma` - diameter of Lennard-Jones particles in the wall + * :math:`\varepsilon` - effective Hamaker constant + + Example:: + + walls = [hoomd.wall.Plane(origin=(0, 0, 0), normal=(0, 1,0))] + lj93 = azplugins.wall.LJ93(walls=walls) + lj93.params["A"] = dict( + sigma=1.0, + epsilon=1.0, + r_cut=3.0, + r_extrap=0.0, + ) + + .. py:attribute:: params + + The `LJ93` potential parameters. The dictionary has the following + keys: + + * ``sigma`` (`float`, **required**) - Lennard-Jones particle size + :math:`\sigma` :math:`[\mathrm{length}]` + * ``epsilon`` (`float`, **required**) - energy parameter + :math:`\varepsilon` :math:`[\mathrm{energy}]` + * ``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 = "WallsPotentialLJ93" + + def __init__(self, walls): + super().__init__(walls) + params = TypeParameter( + "params", + "particle_types", + TypeParameterDict(sigma=float, epsilon=float, r_cut=float, r_extrap=float, len_keys=1), + ) + self._add_typeparam(params) \ No newline at end of file From 8eee09fc881e6727cbd371eb7110ce52f413c61e Mon Sep 17 00:00:00 2001 From: Jinny Cha Date: Fri, 5 Dec 2025 15:05:04 -0600 Subject: [PATCH 2/4] Revise unit test and documentation --- doc/index.rst | 1 + doc/module-azplugins-wall.rst | 23 +++++++++++++ src/WallEvaluatorLJ93.h | 4 +-- src/pytest/test_wall.py | 65 ++++------------------------------- 4 files changed, 33 insertions(+), 60 deletions(-) create mode 100644 doc/module-azplugins-wall.rst diff --git a/doc/index.rst b/doc/index.rst index 0345cbb0..b840d45a 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -56,6 +56,7 @@ Contents module-azplugins-external module-azplugins-flow module-azplugins-pair + module-azplugins-wall .. toctree:: :maxdepth: 1 diff --git a/doc/module-azplugins-wall.rst b/doc/module-azplugins-wall.rst new file mode 100644 index 00000000..86cdab8b --- /dev/null +++ b/doc/module-azplugins-wall.rst @@ -0,0 +1,23 @@ +.. Copyright (c) 2018-2020, Michael P. Howard +.. Copyright (c) 2021-2025, Auburn University +.. Part of azplugins, released under the BSD 3-Clause License. + +azplugins.wall +-------------- + +.. rubric:: Overview + +.. py:currentmodule:: hoomd.azplugins.wall + +.. autosummary:: + :nosignatures: + + LJ93 + +.. rubric:: Details + +.. automodule:: hoomd.azplugins.wall + :synopsis: Wall potentials. + :members: LJ93 + :no-inherited-members: + :show-inheritance: diff --git a/src/WallEvaluatorLJ93.h b/src/WallEvaluatorLJ93.h index 5695543b..97ccf8f5 100644 --- a/src/WallEvaluatorLJ93.h +++ b/src/WallEvaluatorLJ93.h @@ -52,9 +52,9 @@ struct WallParametersLJ93 : public PairParameters } #if HOOMD_LONGREAL_SIZE == 32 - __attribute__((aligned(16))); + __attribute__((aligned(8))); #else - __attribute__((aligned(32))); + __attribute__((aligned(16))); #endif //! Class for evaluatring the Lennard-Jones 9-3 wall force diff --git a/src/pytest/test_wall.py b/src/pytest/test_wall.py index cad64ff0..71686164 100644 --- a/src/pytest/test_wall.py +++ b/src/pytest/test_wall.py @@ -14,80 +14,28 @@ PotentialTestCase = collections.namedtuple( "PotentialTestCase", - ["potential", "params", "wall", "position", "energy", "force"], + ["potential", "params", "position", "energy", "force"], ) potential_tests = [] # Lennard Jones 9-3 potential_tests += [ - # test the calculation of force and potential for plane + # test the calculation of force and potential PotentialTestCase( hoomd.azplugins.wall.LJ93, - {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0, "r_extrap": 0.0}, - hoomd.wall.Plane(origin=(0, 0, 0), normal=(0, 0, 1)), + {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0}, (0, 0, 1.5), -0.2558, -0.5718, ), PotentialTestCase( hoomd.azplugins.wall.LJ93, - {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0, "r_extrap": 0.0}, - hoomd.wall.Plane(origin=(0, 0, 0), normal=(0, 0, 1)), + {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0}, (0, 0, 3.5), 0, 0, ), - # test the calculation of force and potential for sphere - PotentialTestCase( - hoomd.azplugins.wall.LJ93, - {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0, "r_extrap": 0.0}, - hoomd.wall.Sphere(radius=5, origin=(0, 0, 0), inside=True), - (0, 0, 3.5), - -0.2558, - 0.5718, - ), - PotentialTestCase( - hoomd.azplugins.wall.LJ93, - {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0, "r_extrap": 0.0}, - hoomd.wall.Sphere(radius=5, origin=(0, 0, 0), inside=True), - (0, 0, 2.0), - 0, - 0, - ), - PotentialTestCase( - hoomd.azplugins.wall.LJ93, - {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0, "r_extrap": 0.0}, - hoomd.wall.Sphere(radius=5, origin=(0, 0, 0), inside=True), - (0, 0, 5.5), - 0, - 0, - ), - # test the calculation of force and potential for cylinder - PotentialTestCase( - hoomd.azplugins.wall.LJ93, - {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0, "r_extrap": 0.0}, - hoomd.wall.Cylinder(radius=5, origin=(0, 0, 0), axis=(1, 0, 0), inside=True), - (0, 0, 3.5), - -0.2558, - 0.5718, - ), - PotentialTestCase( - hoomd.azplugins.wall.LJ93, - {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0, "r_extrap": 0.0}, - hoomd.wall.Cylinder(radius=5, origin=(0, 0, 0), axis=(1, 0, 0), inside=True), - (0, 0, 5.5), - 0, - 0, - ), - PotentialTestCase( - hoomd.azplugins.wall.LJ93, - {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0, "r_extrap": 0.0}, - hoomd.wall.Sphere(radius=5, origin=(0, 0, 0), inside=True), - (0, 0, 2.0), - 0, - 0, - ), ] @@ -109,7 +57,7 @@ def test_energy_and_force( integrator.methods = [nve] # setup wall potential - wall = potential_test.wall + wall = hoomd.wall.Plane(origin=(0, 0, 0), normal=(0, 0, 1)) potential = potential_test.potential(walls=[wall]) potential.params["A"] = potential_test.params integrator.forces = [potential] @@ -119,7 +67,8 @@ def test_energy_and_force( sim.run(0) # test that parameters are still correct after attach runs - assert potential.params["A"] == potential_test.params + for attr in potential_test.params: + assert potential.params["A"][attr] == potential_test.params[attr] # test that the energies match reference values, half goes to each particle energies = potential.energies From 39006284b6ff68f80cb7cdbb5f7500668c33ea60 Mon Sep 17 00:00:00 2001 From: Jinny Cha Date: Mon, 8 Dec 2025 17:53:32 -0600 Subject: [PATCH 3/4] Rename parameter epsilon to A --- src/WallEvaluatorLJ93.h | 34 +++++++++++++++++----------------- src/pytest/test_wall.py | 4 ++-- src/wall.py | 18 ++++++++++-------- 3 files changed, 29 insertions(+), 27 deletions(-) diff --git a/src/WallEvaluatorLJ93.h b/src/WallEvaluatorLJ93.h index 97ccf8f5..662d629b 100644 --- a/src/WallEvaluatorLJ93.h +++ b/src/WallEvaluatorLJ93.h @@ -28,12 +28,12 @@ namespace detail struct WallParametersLJ93 : public PairParameters { #ifndef __HIPCC__ - WallParametersLJ93() : sigma_3(0), epsilon(0) { } + WallParametersLJ93() : sigma_3(0), A(0) { } WallParametersLJ93(pybind11::dict v, bool managed = false) { auto sigma(v["sigma"].cast()); - epsilon = v["epsilon"].cast(); + A = v["A"].cast(); sigma_3 = sigma * sigma * sigma; } @@ -42,13 +42,13 @@ struct WallParametersLJ93 : public PairParameters { pybind11::dict v; v["sigma"] = std::cbrt(sigma_3); - v["epsilon"] = epsilon; + v["A"] = A; return v; } #endif // __HIPCC__ - Scalar sigma_3; - Scalar epsilon; + Scalar sigma_3; //!< Lennard-Jones sigma + Scalar A; //!< Hamaker constant } #if HOOMD_LONGREAL_SIZE == 32 @@ -57,25 +57,25 @@ struct WallParametersLJ93 : public PairParameters __attribute__((aligned(16))); #endif -//! Class for evaluatring the Lennard-Jones 9-3 wall force +//! Class for evaluating the Lennard-Jones 9-3 wall force /*! * WallEvaluatorLJ93 computes the Lennard-Jones 9-3 wall potential, which is derived from * integrating the standard Lennard-Jones potential between a point particle and a half plane: * - * \f[ V(r) = \varepsilon \left( \frac{2}{15}\left(\frac{\sigma}{r}\right)^9 - - * \left(\frac{\sigma}{r}\right)^3 \right) \f] + * \f[ V(r) = A \left[ \frac{2}{15}\left(\frac{\sigma}{r}\right)^9 - + * \left(\frac{\sigma}{r}\right)^3 \right] \f] * - * where \f$\sigma\f$ is the diameter of Lennard-Jones particles in the wall, and \f$ \varepsilon - * \f$ is the effective Hamaker constant \f$ \varepsilon = (2/3) \pi \varepsilon_{\rm LJ} \rho_{\rm - * w} \sigma^3 \f$ with \f$\varepsilon_{\rm LJ}\f$ the energy of interaction and \f$\rho_{\rm w}\f$ - * the density of particles in the wall. Evaluation of this energy is simplified into the following + * where \f$\sigma\f$ is the diameter of Lennard-Jones particles in the wall, and \f$ A \f$ is the + * Hamaker constant \f$ A = (2/3) \pi \varepsilon_{\rm LJ} \rho_{\rm w} \sigma^3 \f$ with + * \f$\varepsilon_{\rm LJ}\f$ the energy of interaction and \f$\rho_{\rm w}\f$ the density of + * particles in the wall. Evaluation of this energy is simplified into the following * parameters: * - * - \verbatim lj1 = (2.0/15.0) * epsilon * pow(sigma,9.0) \endverbatim - * - \verbatim lj2 = epsilon * pow(sigma,3.0) \endverbatim + * - \verbatim lj1 = (2.0/15.0) * A * pow(sigma,9.0) \endverbatim + * - \verbatim lj2 = A * pow(sigma,3.0) \endverbatim * * The force acting on the particle is then - * \f[ F(r)/r = \frac{\varepsilon}{r^2} \left ( \frac{6}{5}\left(\frac{\sigma}{r}\right)^9 - 3 + * \f[ F(r)/r = \frac{A}{r^2} \left ( \frac{6}{5}\left(\frac{\sigma}{r}\right)^9 - 3 * \left(\frac{\sigma}{r}\right)^3 \right) \f] */ class WallEvaluatorLJ93 : public PairEvaluator @@ -94,9 +94,9 @@ class WallEvaluatorLJ93 : public PairEvaluator DEVICE WallEvaluatorLJ93(Scalar _rsq, Scalar _rcutsq, const param_type& _params) : PairEvaluator(_rsq, _rcutsq) { - lj1 = (Scalar(2.0) / Scalar(15.0)) * _params.epsilon * _params.sigma_3 * _params.sigma_3 + lj1 = (Scalar(2.0) / Scalar(15.0)) * _params.A * _params.sigma_3 * _params.sigma_3 * _params.sigma_3; - lj2 = _params.epsilon * _params.sigma_3; + lj2 = _params.A * _params.sigma_3; } //! Evaluate the force and energy diff --git a/src/pytest/test_wall.py b/src/pytest/test_wall.py index 71686164..98add731 100644 --- a/src/pytest/test_wall.py +++ b/src/pytest/test_wall.py @@ -24,14 +24,14 @@ # test the calculation of force and potential PotentialTestCase( hoomd.azplugins.wall.LJ93, - {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0}, + {"sigma": 1.0, "A": 1.0, "r_cut": 3.0}, (0, 0, 1.5), -0.2558, -0.5718, ), PotentialTestCase( hoomd.azplugins.wall.LJ93, - {"sigma": 1.0, "epsilon": 1.0, "r_cut": 3.0}, + {"sigma": 1.0, "A": 1.0, "r_cut": 3.0}, (0, 0, 3.5), 0, 0, diff --git a/src/wall.py b/src/wall.py index 9c637156..a0bce050 100644 --- a/src/wall.py +++ b/src/wall.py @@ -21,23 +21,25 @@ class LJ93(wall.WallPotential): .. math:: - U(r) = \varepsilon \left(\frac{2}{15} \left(\frac{\sigma}{r} \right)^9 - - \left(\frac{\sigma}{r}\right)^3 \right)] + U(r) = A \left[\frac{2}{15} \left(\frac{\sigma}{r} \right)^9 + - \left(\frac{\sigma}{r}\right)^3 \right] where: * :math:`\sigma` - diameter of Lennard-Jones particles in the wall - * :math:`\varepsilon` - effective Hamaker constant + * :math:`A` - Hamaker constant, related to the Lennard-Jones parameters as + :math:`A = \frac{2}{3} \pi \varepsilon \sigma^3 \rho`, where :math:`\rho` + is the number density of particles in the wall. Example:: walls = [hoomd.wall.Plane(origin=(0, 0, 0), normal=(0, 1, 0))] lj93 = azplugins.wall.LJ93(walls=walls) + # For rho = 1, epsilon = 1, sigma = 1, A = 2*pi/3 lj93.params["A"] = dict( sigma=1.0, - epsilon=1.0, + A=2 * numpy.pi / 3, r_cut=3.0, - r_extrap=0.0, ) .. py:attribute:: params @@ -47,8 +49,8 @@ class LJ93(wall.WallPotential): * ``sigma`` (`float`, **required**) - Lennard-Jones particle size :math:`\sigma` :math:`[\mathrm{length}]` - * ``epsilon`` (`float`, **required**) - energy parameter - :math:`\varepsilon` :math:`[\mathrm{energy}]` + * ``A`` (`float`, **required**) - Hamaker constant + :math:`A` :math:`[\mathrm{energy}]` * ``r_cut`` (`float`, **required**) - The cut off distance for the wall potential :math:`[\mathrm{length}]` * ``r_extrap`` (`float`, **optional**) - The distance to @@ -67,7 +69,7 @@ def __init__(self, walls): "params", "particle_types", TypeParameterDict( - sigma=float, epsilon=float, r_cut=float, r_extrap=0.0, len_keys=1 + sigma=float, A=float, r_cut=float, r_extrap=0.0, len_keys=1 ), ) self._add_typeparam(params) From 90d306a002edbf5bbe1e6e6519d313efbb1c4011 Mon Sep 17 00:00:00 2001 From: Jinny Cha Date: Sun, 21 Dec 2025 00:58:54 -0600 Subject: [PATCH 4/4] Address review comments --- src/WallEvaluatorColloid.h | 40 ++++++++++++++++++++++---------------- src/wall.py | 4 ++-- 2 files changed, 25 insertions(+), 19 deletions(-) diff --git a/src/WallEvaluatorColloid.h b/src/WallEvaluatorColloid.h index b91dffb0..29cc7294 100644 --- a/src/WallEvaluatorColloid.h +++ b/src/WallEvaluatorColloid.h @@ -29,29 +29,35 @@ namespace detail struct WallParametersColloid : public PairParameters { #ifndef __HIPCC__ - WallParametersColloid() : A(0), a(0), sigma_3(0) { } + WallParametersColloid() : c_1(0), c_2(0), a(0) { } WallParametersColloid(pybind11::dict v, bool managed = false) { - A = v["A"].cast(); - a = v["a"].cast(); + auto A(v["A"].cast()); auto sigma(v["sigma"].cast()); - sigma_3 = sigma * sigma * sigma; + + 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"] = std::cbrt(sigma_3); + v["sigma"] = pow(sigma_6, 1. / 6.); return v; } #endif // __HIPCC__ - Scalar A; //!< Hamaker constant - Scalar a; //!< particle radius - Scalar sigma_3; //!< Lennard-Jones sigma + Scalar c_1; //!< First coefficient + Scalar c_2; //!< Second coefficient + Scalar a; //!< Particle radius } #if HOOMD_LONGREAL_SIZE == 32 @@ -99,9 +105,9 @@ class WallEvaluatorColloid : public PairEvaluator DEVICE WallEvaluatorColloid(Scalar _rsq, Scalar _rcutsq, const param_type& _params) : PairEvaluator(_rsq, _rcutsq) { + c_1 = _params.c_1; + c_2 = _params.c_2; a = _params.a; - lj1 = _params.A * _params.sigma_3 * _params.sigma_3 / Scalar(7560); - lj2 = _params.A / Scalar(6); } //! Computes the colloid-wall interaction potential @@ -138,18 +144,18 @@ class WallEvaluatorColloid : public PairEvaluator if (force) { Scalar arinv8 = Scalar(8.0) * arinv; - force_divr = Scalar(6.0) * lj1 + 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 -= lj2 * (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 = lj1 + 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 -= lj2 * (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; } @@ -166,7 +172,7 @@ class WallEvaluatorColloid : public PairEvaluator */ DEVICE bool evalForceAndEnergy(Scalar& force_divr, Scalar& energy, bool energy_shift) { - if (rsq < rcutsq && lj1 != 0 && a > 0) + if (rsq < rcutsq && c_1 != 0 && a > 0) { energy = computePotential(force_divr, rsq); if (energy_shift) @@ -188,9 +194,9 @@ class WallEvaluatorColloid : public PairEvaluator #endif 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 - Scalar lj1; //!< Prefactor for first term (includes Hamaker constant) - Scalar lj2; //!< Prefactor for second term (includes Hamaker constant) }; } // end namespace detail diff --git a/src/wall.py b/src/wall.py index 3fd0c0d9..d9a86392 100644 --- a/src/wall.py +++ b/src/wall.py @@ -17,8 +17,8 @@ class Colloid(wall.WallPotential): 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 :math:`a` - and a half plane: + the standard Lennard-Jones potential between a spherical particle of radius + :math:`a` and a half plane: .. math::