From f4aabd3cde47e67af4d0754e89bdfae4c3bffa67 Mon Sep 17 00:00:00 2001 From: mayukh33 Date: Thu, 7 May 2026 17:45:36 +0200 Subject: [PATCH 01/17] Add two-dimensional linear interpolation --- src/LinearInterpolator2D.h | 134 +++++++++++++++++++++++++++++++++++++ 1 file changed, 134 insertions(+) create mode 100644 src/LinearInterpolator2D.h diff --git a/src/LinearInterpolator2D.h b/src/LinearInterpolator2D.h new file mode 100644 index 00000000..2906f883 --- /dev/null +++ b/src/LinearInterpolator2D.h @@ -0,0 +1,134 @@ +#ifndef AZPLUGINS_LINEAR_INTERPOLATOR_2D_H_ +#define AZPLUGINS_LINEAR_INTERPOLATOR_2D_H_ + +#include +#include +#include + +#include "hoomd/HOOMDMath.h" +#include "hoomd/Index1D.h" + +#if defined(__HIPCC__) || defined(__CUDACC__) +#define AZPLUGINS_HOSTDEVICE __host__ __device__ +#define AZPLUGINS_FORCEINLINE __forceinline__ +#else +#define AZPLUGINS_HOSTDEVICE +#define AZPLUGINS_FORCEINLINE inline +#endif + +namespace hoomd + { +namespace azplugins + { + +template class LinearInterpolator2D + { + public: + AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE LinearInterpolator2D() : m_data(nullptr), m_indexer() + { + m_n[0] = 0; + m_n[1] = 0; + for (int d = 0; d < 2; ++d) + { + m_lo[d] = Scalar(0); + m_hi[d] = Scalar(0); + m_dx[d] = Scalar(0); + } + } + + //! Constructor + AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE + LinearInterpolator2D(const T* data, + const unsigned int* n, + const Scalar* lo, + const Scalar* hi) + : m_data(data), m_indexer(n[0], n[1]) + { + for (int d = 0; d < 2; ++d) + { + const unsigned int nd = n[d]; + assert(nd >= 2); + + m_n[d] = nd; + m_lo[d] = lo[d]; + m_hi[d] = hi[d]; + m_dx[d] = (m_hi[d] - m_lo[d]) / Scalar(nd - 1); + } + } + + //! Destructor + ~LinearInterpolator2D() + { + delete[] m_data; + } + + //! Interpolate at (x0, x1) + AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar operator()(Scalar x0, Scalar x1) const + { + const Scalar x[2] = {x0, x1}; + + int bin[2]; + Scalar frac[2]; + + for (int d = 0; d < 2; ++d) + { + const unsigned int nd = m_n[d]; + const Scalar f = (x[d] - m_lo[d]) / m_dx[d]; + int b = static_cast(std::floor(f)); + + // If exactly at the top boundary, shift into the last valid cell + if (f == Scalar(nd - 1) && x[d] == m_hi[d]) + { + --b; + } + + assert(b >= 0); + assert(b < static_cast(nd - 1)); + bin[d] = b; + frac[d] = f - Scalar(b); + } + + const unsigned int b0 = static_cast(bin[0]); + const unsigned int b1 = static_cast(bin[1]); + const Scalar c00 = m_data[m_indexer(b0, b1)]; + const Scalar c01 = m_data[m_indexer(b0, b1 + 1)]; + const Scalar c10 = m_data[m_indexer(b0 + 1, b1)]; + const Scalar c11 = m_data[m_indexer(b0 + 1, b1 + 1)]; + + // Bilinear interpolation + const Scalar t0 = frac[0]; + const Scalar omt0 = Scalar(1) - t0; + const Scalar c0 = c00 * omt0 + c10 * t0; + const Scalar c1 = c01 * omt0 + c11 * t0; + + const Scalar t1 = frac[1]; + return c0 * (Scalar(1) - t1) + c1 * t1; + } + + //! Lower bound for a given dimension + AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar getLo(int dim) const + { + return m_lo[dim]; + } + //! Upper bound for a given dimension + AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar getHi(int dim) const + { + return m_hi[dim]; + } + + private: + const T* m_data; + Scalar m_lo[2]; + Scalar m_hi[2]; + Scalar m_dx[2]; + unsigned int m_n[2]; + Index2D m_indexer; + }; + + } // namespace azplugins + } // namespace hoomd + +#undef AZPLUGINS_HOSTDEVICE +#undef AZPLUGINS_FORCEINLINE + +#endif // AZPLUGINS_LINEAR_INTERPOLATOR_2D_H_ From 3fc3f51217dbd813e87e09b4901efde3dda7db4d Mon Sep 17 00:00:00 2001 From: mayukh33 Date: Thu, 7 May 2026 17:46:57 +0200 Subject: [PATCH 02/17] Add header file evaporation dependent PerturbedLJ potential --- src/PerturbedLennardJonesEvap.h | 101 ++++++++++++++++++++++++++++++++ 1 file changed, 101 insertions(+) create mode 100644 src/PerturbedLennardJonesEvap.h diff --git a/src/PerturbedLennardJonesEvap.h b/src/PerturbedLennardJonesEvap.h new file mode 100644 index 00000000..39db8f2a --- /dev/null +++ b/src/PerturbedLennardJonesEvap.h @@ -0,0 +1,101 @@ +// 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_PERTURBED_LENNARD_JONES_EVAP_H_ +#define AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_H_ + +#include +#include +#include + +#include "hoomd/BoxDim.h" +#include "hoomd/ForceCompute.h" +#include "hoomd/GPUArray.h" +#include "hoomd/HOOMDMath.h" +#include "hoomd/Index1D.h" +#include "hoomd/VectorMath.h" +#include "hoomd/md/NeighborList.h" + +#include "LinearInterpolator2D.h" + +namespace hoomd + { +namespace azplugins + { + +class PYBIND11_EXPORT PerturbedLennardJonesEvap : public ForceCompute + { + public: + PerturbedLennardJonesEvap(std::shared_ptr sysdef, + std::shared_ptr nlist, + Scalar r_cut, + Scalar epsilon, + Scalar sigma, + Scalar dt, + Scalar initial_delta, + Scalar solventD, + const Scalar* lambda_data, + const unsigned int* lambda_shape, + const Scalar* domain); + + virtual ~PerturbedLennardJonesEvap(); + + Scalar getRCut() const + { + return m_r_cut; + } + + Scalar getEpsilon() const + { + return m_epsilon; + } + + Scalar getSigma() const + { + return m_sigma; + } + + Scalar getInitial_delta() const + { + return m_initial_delta; + } + + Scalar getSolventD() const + { + return m_solventD; + } + + const GPUArray& getLambdaDomain() const + { + return m_domain; + } + + protected: + std::shared_ptr m_nlist; //!< Neighbor list + + // Lennard-Jones parameters + Scalar m_r_cut; + Scalar m_r_cutsq; + Scalar m_epsilon; + Scalar m_sigma; + Scalar m_dt; + Scalar m_initial_delta; + Scalar m_solventD; + + GPUArray m_domain; //!< Scaled coordinate domain [y_lo, y_hi, t_lo, t_hi] + GPUArray m_lambda_data; //!< Flattened scaled position scaled time data + GPUArray m_lambda_shape; //!< [ny, nt] + + void computeForces(uint64_t timestep) override; + }; + +namespace detail + { +void export_PerturbedLennardJonesEvap(pybind11::module& m); + } // end namespace detail + + } // end namespace azplugins + } // end namespace hoomd + +#endif // AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_H_ From c019b9f8757b5a1c441bfb458a949562e991d0c2 Mon Sep 17 00:00:00 2001 From: mayukh33 Date: Thu, 7 May 2026 18:11:59 +0200 Subject: [PATCH 03/17] Fix pre-commit error --- src/LinearInterpolator2D.h | 41 ++++++++++++------------ src/PerturbedLennardJonesEvap.h | 56 ++++++++++++++++----------------- 2 files changed, 49 insertions(+), 48 deletions(-) diff --git a/src/LinearInterpolator2D.h b/src/LinearInterpolator2D.h index 2906f883..1d1e8532 100644 --- a/src/LinearInterpolator2D.h +++ b/src/LinearInterpolator2D.h @@ -1,3 +1,7 @@ +// 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_LINEAR_INTERPOLATOR_2D_H_ #define AZPLUGINS_LINEAR_INTERPOLATOR_2D_H_ @@ -38,10 +42,7 @@ template class LinearInterpolator2D //! Constructor AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE - LinearInterpolator2D(const T* data, - const unsigned int* n, - const Scalar* lo, - const Scalar* hi) + LinearInterpolator2D(const T* data, const unsigned int* n, const Scalar* lo, const Scalar* hi) : m_data(data), m_indexer(n[0], n[1]) { for (int d = 0; d < 2; ++d) @@ -49,7 +50,7 @@ template class LinearInterpolator2D const unsigned int nd = n[d]; assert(nd >= 2); - m_n[d] = nd; + m_n[d] = nd; m_lo[d] = lo[d]; m_hi[d] = hi[d]; m_dx[d] = (m_hi[d] - m_lo[d]) / Scalar(nd - 1); @@ -58,9 +59,9 @@ template class LinearInterpolator2D //! Destructor ~LinearInterpolator2D() - { - delete[] m_data; - } + { + delete[] m_data; + } //! Interpolate at (x0, x1) AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar operator()(Scalar x0, Scalar x1) const @@ -81,22 +82,22 @@ template class LinearInterpolator2D { --b; } - - assert(b >= 0); - assert(b < static_cast(nd - 1)); - bin[d] = b; + + assert(b >= 0); + assert(b < static_cast(nd - 1)); + bin[d] = b; frac[d] = f - Scalar(b); } const unsigned int b0 = static_cast(bin[0]); const unsigned int b1 = static_cast(bin[1]); - const Scalar c00 = m_data[m_indexer(b0, b1)]; - const Scalar c01 = m_data[m_indexer(b0, b1 + 1)]; + const Scalar c00 = m_data[m_indexer(b0, b1)]; + const Scalar c01 = m_data[m_indexer(b0, b1 + 1)]; const Scalar c10 = m_data[m_indexer(b0 + 1, b1)]; const Scalar c11 = m_data[m_indexer(b0 + 1, b1 + 1)]; // Bilinear interpolation - const Scalar t0 = frac[0]; + const Scalar t0 = frac[0]; const Scalar omt0 = Scalar(1) - t0; const Scalar c0 = c00 * omt0 + c10 * t0; const Scalar c1 = c01 * omt0 + c11 * t0; @@ -117,12 +118,12 @@ template class LinearInterpolator2D } private: - const T* m_data; - Scalar m_lo[2]; - Scalar m_hi[2]; - Scalar m_dx[2]; + const T* m_data; + Scalar m_lo[2]; + Scalar m_hi[2]; + Scalar m_dx[2]; unsigned int m_n[2]; - Index2D m_indexer; + Index2D m_indexer; }; } // namespace azplugins diff --git a/src/PerturbedLennardJonesEvap.h b/src/PerturbedLennardJonesEvap.h index 39db8f2a..ad34461a 100644 --- a/src/PerturbedLennardJonesEvap.h +++ b/src/PerturbedLennardJonesEvap.h @@ -6,8 +6,8 @@ #define AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_H_ #include -#include #include +#include #include "hoomd/BoxDim.h" #include "hoomd/ForceCompute.h" @@ -32,9 +32,9 @@ class PYBIND11_EXPORT PerturbedLennardJonesEvap : public ForceCompute Scalar r_cut, Scalar epsilon, Scalar sigma, - Scalar dt, - Scalar initial_delta, - Scalar solventD, + Scalar dt, + Scalar initial_delta, + Scalar solventD, const Scalar* lambda_data, const unsigned int* lambda_shape, const Scalar* domain); @@ -42,34 +42,34 @@ class PYBIND11_EXPORT PerturbedLennardJonesEvap : public ForceCompute virtual ~PerturbedLennardJonesEvap(); Scalar getRCut() const - { - return m_r_cut; - } + { + return m_r_cut; + } Scalar getEpsilon() const - { - return m_epsilon; - } - + { + return m_epsilon; + } + Scalar getSigma() const - { - return m_sigma; - } + { + return m_sigma; + } Scalar getInitial_delta() const - { - return m_initial_delta; - } + { + return m_initial_delta; + } Scalar getSolventD() const - { - return m_solventD; - } - - const GPUArray& getLambdaDomain() const - { - return m_domain; - } + { + return m_solventD; + } + + const GPUArray& getLambdaDomain() const + { + return m_domain; + } protected: std::shared_ptr m_nlist; //!< Neighbor list @@ -83,9 +83,9 @@ class PYBIND11_EXPORT PerturbedLennardJonesEvap : public ForceCompute Scalar m_initial_delta; Scalar m_solventD; - GPUArray m_domain; //!< Scaled coordinate domain [y_lo, y_hi, t_lo, t_hi] - GPUArray m_lambda_data; //!< Flattened scaled position scaled time data - GPUArray m_lambda_shape; //!< [ny, nt] + GPUArray m_domain; //!< Scaled coordinate domain [y_lo, y_hi, t_lo, t_hi] + GPUArray m_lambda_data; //!< Flattened scaled position scaled time data + GPUArray m_lambda_shape; //!< [ny, nt] void computeForces(uint64_t timestep) override; }; From 04abad5fd76d9619218bdf31a6b73df0f587111d Mon Sep 17 00:00:00 2001 From: mayukh33 Date: Wed, 13 May 2026 14:32:15 +0200 Subject: [PATCH 04/17] Implement perturbedLJ for evaporation with interpolated variant --- src/CMakeLists.txt | 3 + src/LinearInterpolator1D.h | 92 +++++++ src/LinearInterpolator2D.h | 11 +- src/PerturbedLennardJonesEvap.h | 311 +++++++++++++++++++++--- src/VariantInterpolated.h | 111 +++++++++ src/__init__.py | 2 +- src/export_PerturbedLennardJonesEvap.cc | 71 ++++++ src/export_VariantInterpolated.cc | 51 ++++ src/module.cc | 7 + src/pair.py | 64 +++++ src/variant.py | 28 +++ 11 files changed, 706 insertions(+), 45 deletions(-) create mode 100644 src/LinearInterpolator1D.h create mode 100644 src/VariantInterpolated.h create mode 100644 src/export_PerturbedLennardJonesEvap.cc create mode 100644 src/export_VariantInterpolated.cc create mode 100644 src/variant.py diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a5529d98..7f0ea7eb 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -5,6 +5,8 @@ set(COMPONENT_NAME azplugins) set(_${COMPONENT_NAME}_sources ConstantFlow.cc export_ImagePotentialBondHarmonic.cc + export_PerturbedLennardJonesEvap.cc + export_VariantInterpolated.cc module.cc ParabolicFlow.cc VelocityCompute.cc @@ -32,6 +34,7 @@ set(python_files external.py flow.py pair.py + variant.py wall.py ) diff --git a/src/LinearInterpolator1D.h b/src/LinearInterpolator1D.h new file mode 100644 index 00000000..5ed7dc7d --- /dev/null +++ b/src/LinearInterpolator1D.h @@ -0,0 +1,92 @@ +// 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_LINEAR_INTERPOLATOR_1D_H_ +#define AZPLUGINS_LINEAR_INTERPOLATOR_1D_H_ + +#include +#include +#include + +#include "hoomd/HOOMDMath.h" + +#if defined(__HIPCC__) || defined(__CUDACC__) +#define AZPLUGINS_HOSTDEVICE __host__ __device__ +#define AZPLUGINS_FORCEINLINE __forceinline__ +#else +#define AZPLUGINS_HOSTDEVICE +#define AZPLUGINS_FORCEINLINE inline +#endif + +namespace hoomd + { +namespace azplugins + { + +template class LinearInterpolator1D + { + public: + AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE LinearInterpolator1D() + : m_data(nullptr), m_lo(Scalar(0)), m_hi(Scalar(0)), m_dx(Scalar(0)), m_n(0) + { + } + + //! Constructor + AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE + LinearInterpolator1D(const T* data, unsigned int n, Scalar lo, Scalar hi) + : m_data(data), m_lo(lo), m_hi(hi), m_n(n) + { + assert(n >= 2); + m_dx = (m_hi - m_lo) / Scalar(n - 1); + } + + //! Implement piecewise linear interpolation + AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar operator()(Scalar x) const + { + const Scalar f = (x - m_lo) / m_dx; + int b = static_cast(std::floor(f)); + + // If exactly at the top boundary, shift into the last valid cell + if (f == Scalar(m_n - 1) && x == m_hi) + { + --b; + } + + assert(b >= 0); + assert(b < static_cast(m_n - 1)); + const Scalar frac = f - Scalar(b); + + const Scalar c0 = m_data[b]; + const Scalar c1 = m_data[b + 1]; + + // Linear interpolation + return c0 * (Scalar(1) - frac) + c1 * frac; + } + + //! Lower bound + AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar getLo() const + { + return m_lo; + } + //! Upper bound + AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar getHi() const + { + return m_hi; + } + + private: + const T* m_data; + Scalar m_lo; + Scalar m_hi; + Scalar m_dx; + unsigned int m_n; + }; + + } // namespace azplugins + } // namespace hoomd + +#undef AZPLUGINS_HOSTDEVICE +#undef AZPLUGINS_FORCEINLINE + +#endif // AZPLUGINS_LINEAR_INTERPOLATOR_1D_H_ diff --git a/src/LinearInterpolator2D.h b/src/LinearInterpolator2D.h index 1d1e8532..e8de081a 100644 --- a/src/LinearInterpolator2D.h +++ b/src/LinearInterpolator2D.h @@ -57,12 +57,6 @@ template class LinearInterpolator2D } } - //! Destructor - ~LinearInterpolator2D() - { - delete[] m_data; - } - //! Interpolate at (x0, x1) AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar operator()(Scalar x0, Scalar x1) const { @@ -106,12 +100,13 @@ template class LinearInterpolator2D return c0 * (Scalar(1) - t1) + c1 * t1; } - //! Lower bound for a given dimension + //! Lower bound AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar getLo(int dim) const { return m_lo[dim]; } - //! Upper bound for a given dimension + + //! Upper bound AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar getHi(int dim) const { return m_hi[dim]; diff --git a/src/PerturbedLennardJonesEvap.h b/src/PerturbedLennardJonesEvap.h index ad34461a..b3e1ed0b 100644 --- a/src/PerturbedLennardJonesEvap.h +++ b/src/PerturbedLennardJonesEvap.h @@ -5,6 +5,7 @@ #ifndef AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_H_ #define AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_H_ +#include #include #include #include @@ -12,83 +13,321 @@ #include "hoomd/BoxDim.h" #include "hoomd/ForceCompute.h" #include "hoomd/GPUArray.h" -#include "hoomd/HOOMDMath.h" -#include "hoomd/Index1D.h" +#include "hoomd/Variant.h" #include "hoomd/VectorMath.h" #include "hoomd/md/NeighborList.h" #include "LinearInterpolator2D.h" +#include "VariantInterpolated.h" namespace hoomd { namespace azplugins { -class PYBIND11_EXPORT PerturbedLennardJonesEvap : public ForceCompute +struct PairParametersPerturbedLennardJonesEvap + { + Scalar epsilon_x_4; + Scalar sigma_6; + Scalar rwcasq; + +#ifndef __HIPCC__ + + //! Default constructor + PairParametersPerturbedLennardJonesEvap() : epsilon_x_4(0), sigma_6(0), rwcasq(0) { } + + PairParametersPerturbedLennardJonesEvap(Scalar epsilon, Scalar sigma) + { + const Scalar sigma_2 = sigma * sigma; + const Scalar sigma_4 = sigma_2 * sigma_2; + sigma_6 = sigma_2 * sigma_4; + epsilon_x_4 = Scalar(4.0) * epsilon; + rwcasq = std::pow(Scalar(2.0), Scalar(1.0) / Scalar(3.0)) * sigma_2; + } + + PairParametersPerturbedLennardJonesEvap(pybind11::dict v, bool managed = false) + { + auto sigma = v["sigma"].cast(); + auto epsilon = v["epsilon"].cast(); + + const Scalar sigma_2 = sigma * sigma; + const Scalar sigma_4 = sigma_2 * sigma_2; + sigma_6 = sigma_2 * sigma_4; + epsilon_x_4 = Scalar(4.0) * epsilon; + rwcasq = std::pow(Scalar(2.0), Scalar(1.0) / Scalar(3.0)) * sigma_2; + } + + pybind11::dict asDict() + { + pybind11::dict v; + v["sigma"] = std::pow(sigma_6, Scalar(1.0) / Scalar(6.0)); + v["epsilon"] = epsilon_x_4 / Scalar(4.0); + return v; + } +#endif // __HIPCC__ + }; + +class PerturbedLennardJonesEvap : public ForceCompute { public: + typedef PairParametersPerturbedLennardJonesEvap param_type; + PerturbedLennardJonesEvap(std::shared_ptr sysdef, std::shared_ptr nlist, Scalar r_cut, - Scalar epsilon, - Scalar sigma, - Scalar dt, - Scalar initial_delta, - Scalar solventD, - const Scalar* lambda_data, - const unsigned int* lambda_shape, - const Scalar* domain); + Scalar scale_factor, + const param_type& params, + bool energy_shift, + const Scalar* attraction_scale_factor_data, + const unsigned int* attraction_scale_factor_shape, + const Scalar* domain, + std::shared_ptr variant); + + ~PerturbedLennardJonesEvap() { } - virtual ~PerturbedLennardJonesEvap(); + Scalar scaleTime(uint64_t timestep) const + { + return Scalar(static_cast(timestep) / m_scale_factor); + } Scalar getRCut() const { - return m_r_cut; + return m_rcut; } Scalar getEpsilon() const { - return m_epsilon; + return epsilon_x_4 / Scalar(4.0); } Scalar getSigma() const { - return m_sigma; + return std::pow(sigma_6, Scalar(1.0) / Scalar(6.0)); } - Scalar getInitial_delta() const + protected: + std::shared_ptr m_nlist; //!< Neighbor list + Scalar epsilon_x_4; + Scalar m_rcut; + Scalar m_scale_factor; //!< Time scaling factor + Scalar lj1; + Scalar lj2; + Scalar rcutsq; + Scalar rwcasq; + Scalar sigma_6; + bool m_energy_shift; + GPUArray m_domain; //!< [y_lo, y_hi, t_lo, t_hi] + GPUArray m_attraction_scale_factor_data; //!< Flattened (y, t) data + GPUArray m_attraction_scale_factor_shape; //!< [ny, nt] + std::shared_ptr m_variant; + + void computeForces(uint64_t timestep) override; + }; + +PerturbedLennardJonesEvap::PerturbedLennardJonesEvap( + std::shared_ptr sysdef, + std::shared_ptr nlist, + Scalar rcut, + Scalar scale_factor, + const param_type& params, + bool energy_shift, + const Scalar* attraction_scale_factor_data, + const unsigned int* attraction_scale_factor_shape, + const Scalar* domain, + std::shared_ptr variant) + : ForceCompute(sysdef), m_nlist(nlist), epsilon_x_4(params.epsilon_x_4), m_rcut(rcut), + m_scale_factor(scale_factor), lj1(params.epsilon_x_4 * params.sigma_6 * params.sigma_6), + lj2(params.epsilon_x_4 * params.sigma_6), rcutsq(rcut * rcut), rwcasq(params.rwcasq), + sigma_6(params.sigma_6), m_energy_shift(energy_shift), m_variant(variant) + { + // Allocate and fill the (y, t) domain: [y_lo, y_hi, t_lo, t_hi] { - return m_initial_delta; + GPUArray domain_arr(4, m_exec_conf); + m_domain.swap(domain_arr); + + ArrayHandle h_domain(m_domain, access_location::host, access_mode::overwrite); + std::copy(domain, domain + 4, h_domain.data); } - Scalar getSolventD() const + // Allocate and fill the table shape (ny, nt) { - return m_solventD; + GPUArray shape_arr(2, m_exec_conf); + m_attraction_scale_factor_shape.swap(shape_arr); + + ArrayHandle h_shape(m_attraction_scale_factor_shape, + access_location::host, + access_mode::overwrite); + std::copy(attraction_scale_factor_shape, attraction_scale_factor_shape + 2, h_shape.data); } - const GPUArray& getLambdaDomain() const + // Allocate and fill the tabulated data { - return m_domain; + const unsigned int n_data + = attraction_scale_factor_shape[0] * attraction_scale_factor_shape[1]; + + GPUArray attraction_scale_factor_arr(n_data, m_exec_conf); + m_attraction_scale_factor_data.swap(attraction_scale_factor_arr); + + ArrayHandle h_data(m_attraction_scale_factor_data, + access_location::host, + access_mode::overwrite); + + std::copy(attraction_scale_factor_data, attraction_scale_factor_data + n_data, h_data.data); } + } - protected: - std::shared_ptr m_nlist; //!< Neighbor list +void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) + { + m_nlist->compute(timestep); - // Lennard-Jones parameters - Scalar m_r_cut; - Scalar m_r_cutsq; - Scalar m_epsilon; - Scalar m_sigma; - Scalar m_dt; - Scalar m_initial_delta; - Scalar m_solventD; + const bool third_law + = (m_nlist->getStorageMode() == hoomd::md::NeighborList::storageMode::half); - GPUArray m_domain; //!< Scaled coordinate domain [y_lo, y_hi, t_lo, t_hi] - GPUArray m_lambda_data; //!< Flattened scaled position scaled time data - GPUArray m_lambda_shape; //!< [ny, nt] + Scalar scaled_t = scaleTime(timestep); + Scalar interface_height = (*m_variant)(timestep); - void computeForces(uint64_t timestep) override; - }; + ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::read); + m_force.zeroFill(); + + ArrayHandle h_force(m_force, access_location::host, access_mode::overwrite); + ArrayHandle h_data(m_attraction_scale_factor_data, + access_location::host, + access_mode::read); + ArrayHandle h_shape(m_attraction_scale_factor_shape, + access_location::host, + access_mode::read); + ArrayHandle h_domain(m_domain, access_location::host, access_mode::read); + + // Neighbor-list arrays + ArrayHandle h_n_neigh(m_nlist->getNNeighArray(), + access_location::host, + access_mode::read); + ArrayHandle h_nlist(m_nlist->getNListArray(), + access_location::host, + access_mode::read); + ArrayHandle h_head_list(m_nlist->getHeadList(), + access_location::host, + access_mode::read); + + // Build the interpolator: lo = {y_lo, t_lo}, hi = {y_hi, t_hi} + const Scalar lo[2] = {h_domain.data[0], h_domain.data[2]}; + const Scalar hi[2] = {h_domain.data[1], h_domain.data[3]}; + LinearInterpolator2D interp(h_data.data, h_shape.data, lo, hi); + + const BoxDim box = m_pdata->getGlobalBox(); + const unsigned int N = m_pdata->getN(); + + for (unsigned int i = 0; i < N; ++i) + { + const Scalar3 pos_i = make_scalar3(h_pos.data[i].x, h_pos.data[i].y, h_pos.data[i].z); + + Scalar scaled_y_i = Scalar(h_pos.data[i].y / interface_height); + Scalar3 fi = make_scalar3(0, 0, 0); + Scalar pei = 0; + + const Scalar attraction_scale_factor_i = interp(scaled_y_i, scaled_t); + + const unsigned int size = (unsigned int)h_n_neigh.data[i]; + const size_t head = h_head_list.data[i]; + + for (unsigned int k = 0; k < size; ++k) + { + const unsigned int j = h_nlist.data[head + k]; + if (j == i) + continue; + Scalar scaled_y_j = Scalar(h_pos.data[j].y / interface_height); + Scalar3 pos_j = make_scalar3(h_pos.data[j].x, h_pos.data[j].y, h_pos.data[j].z); + + const Scalar attraction_scale_factor_j = interp(scaled_y_j, scaled_t); + + // Minimum-image + Scalar3 dx = pos_i - pos_j; + dx = box.minImage(dx); + + const Scalar rsq = dot(dx, dx); + const Scalar attraction_scale_factor_avg + = Scalar(0.5) * (attraction_scale_factor_i + attraction_scale_factor_j); + const Scalar wca_shift + = epsilon_x_4 * (Scalar(1.0) - attraction_scale_factor_avg) / Scalar(4.); + + if (rsq < rcutsq && lj1 != 0) + { + const Scalar r2inv = Scalar(1) / rsq; + const Scalar r6inv = r2inv * r2inv * r2inv; + + const Scalar pair_eng = r6inv * (lj1 * r6inv - lj2); + const Scalar force_divr + = r6inv * r2inv * (Scalar(12.0) * lj1 * r6inv - Scalar(6.0) * lj2); + + if (rsq < rwcasq) + { + pei += wca_shift * pair_eng; + if (third_law) + { + h_force.data[j].w += wca_shift * pair_eng; + } + } + else + { + fi.x += attraction_scale_factor_avg * force_divr * dx.x; + fi.y += attraction_scale_factor_avg * force_divr * dx.y; + fi.z += attraction_scale_factor_avg * force_divr * dx.z; + pei += attraction_scale_factor_avg * pair_eng; + + if (third_law) + { + h_force.data[j].x -= attraction_scale_factor_avg * force_divr * dx.x; + h_force.data[j].y -= attraction_scale_factor_avg * force_divr * dx.y; + h_force.data[j].z -= attraction_scale_factor_avg * force_divr * dx.z; + h_force.data[j].w += Scalar(0.5) * attraction_scale_factor_avg * pair_eng; + } + } + + if (m_energy_shift) + { + const Scalar rcut2inv = Scalar(1.0) / rcutsq; + const Scalar rcut6inv = rcut2inv * rcut2inv * rcut2inv; + Scalar pair_eng_shift = rcut6inv * (lj1 * rcut6inv - lj2); + const Scalar force_divr_shift + = rcut6inv * rcut2inv * (Scalar(12.0) * lj1 * rcut6inv - Scalar(6.0) * lj2); + + if (rsq < rwcasq) + { + pei += wca_shift * pair_eng_shift; + if (third_law) + { + h_force.data[j].w += Scalar(0.5) * wca_shift * pair_eng_shift; + } + } + else + { + fi.x += attraction_scale_factor_avg * force_divr_shift * dx.x; + fi.y += attraction_scale_factor_avg * force_divr_shift * dx.y; + fi.z += attraction_scale_factor_avg * force_divr_shift * dx.z; + pei += attraction_scale_factor_avg * pair_eng_shift; + + if (third_law) + { + h_force.data[j].x + -= attraction_scale_factor_avg * force_divr_shift * dx.x; + h_force.data[j].y + -= attraction_scale_factor_avg * force_divr_shift * dx.y; + h_force.data[j].z + -= attraction_scale_factor_avg * force_divr_shift * dx.z; + h_force.data[j].w + += Scalar(0.5) * attraction_scale_factor_avg * pair_eng_shift; + } + } + } + } + } + + h_force.data[i].x += fi.x; + h_force.data[i].y += fi.y; + h_force.data[i].z += fi.z; + h_force.data[i].w += Scalar(0.5) * pei; + } + } namespace detail { diff --git a/src/VariantInterpolated.h b/src/VariantInterpolated.h new file mode 100644 index 00000000..9c3455ec --- /dev/null +++ b/src/VariantInterpolated.h @@ -0,0 +1,111 @@ +// Copyright (c) 2018-2020, Michael P. Howard +// Copyright (c) 2021-2025, Auburn University +// Part of azplugins, released under the BSD 3-Clause License. + +/*! + * \file VariantInterpolated.h + * \brief Piecewise-linear variant on a uniform grid in scaled time, + * backed by LinearInterpolator1D. + */ + +#ifndef AZPLUGINS_VARIANT_INTERPOLATED_H_ +#define AZPLUGINS_VARIANT_INTERPOLATED_H_ + +#ifdef __HIPCC__ +#error This header cannot be compiled by nvcc / hipcc +#endif +#include +#include +#include +#include +#include + +#include "hoomd/GPUArray.h" +#include "hoomd/Variant.h" +#include "hoomd/VectorMath.h" + +#include "LinearInterpolator1D.h" + +namespace hoomd + { +namespace azplugins + { + +class PYBIND11_EXPORT VariantInterpolated : public Variant + { + public: + //! Construct from a 1D values array and the time domain. + VariantInterpolated(const Scalar* data, // Input interface height values + unsigned int n, // Size of the array + Scalar t_lo, // Low of the unscaled time + Scalar t_hi) // Hi of the unscaled time + + : m_n(n), m_t_lo(t_lo), m_t_hi(t_hi) + { + GPUArray data_arr(n, m_exec_conf); + m_data.swap(data_arr); + ArrayHandle h_data(m_data, access_location::host, access_mode::overwrite); + std::copy(data, data + n, h_data.data); + } + + Scalar operator()(uint64_t timestep) + { + ArrayHandle h_data(m_data, access_location::host, access_mode::read); + + // LinearInterpolator1D + LinearInterpolator1D interp(h_data.data, m_n, m_t_lo, m_t_hi); + return interp(static_cast(timestep)); + } + + Scalar getTLo() const + { + return m_t_lo; + } + + Scalar getTHi() const + { + return m_t_hi; + } + + Scalar min() + { + ArrayHandle h_data(m_data, hoomd::access_location::host, hoomd::access_mode::read); + m_min = h_data.data[0]; + for (unsigned int i = 1; i < m_n; ++i) + { + if (h_data.data[i] < m_min) + m_min = h_data.data[i]; + } + return m_min; + } + + Scalar max() + { + ArrayHandle h_data(m_data, hoomd::access_location::host, hoomd::access_mode::read); + m_max = h_data.data[0]; + for (unsigned int i = 1; i < m_n; ++i) + { + if (h_data.data[i] > m_max) + m_max = h_data.data[i]; + } + return m_max; + } + + protected: + std::shared_ptr m_exec_conf; + GPUArray m_data; + unsigned int m_n; + Scalar m_t_lo; + Scalar m_t_hi; + Scalar m_min; + Scalar m_max; + }; + +namespace detail + { +void export_VariantInterpolated(pybind11::module& m); + } // namespace detail + } // namespace azplugins + } // namespace hoomd + +#endif // AZPLUGINS_VARIANT_INTERPOLATED_H_ diff --git a/src/__init__.py b/src/__init__.py index 3e71abf0..70dbb57b 100644 --- a/src/__init__.py +++ b/src/__init__.py @@ -4,6 +4,6 @@ """azplugins.""" -from hoomd.azplugins import bond, compute, external, flow, pair, wall +from hoomd.azplugins import bond, compute, external, flow, pair, variant, wall __version__ = "1.2.0" diff --git a/src/export_PerturbedLennardJonesEvap.cc b/src/export_PerturbedLennardJonesEvap.cc new file mode 100644 index 00000000..e47d84c7 --- /dev/null +++ b/src/export_PerturbedLennardJonesEvap.cc @@ -0,0 +1,71 @@ +// Copyright (c) 2018-2020, Michael P. Howard +// Copyright (c) 2021-2025, Auburn University +// Part of azplugins, released under the BSD 3-Clause License. + +#include "PerturbedLennardJonesEvap.h" + +#include + +namespace hoomd + { +namespace azplugins + { +namespace detail + { +namespace py = pybind11; + +void export_PerturbedLennardJonesEvap(py::module& m) + { + py::class_>( + m, + "PerturbedLennardJonesEvap") + .def(py::init( + [](std::shared_ptr sysdef, + std::shared_ptr nlist, + Scalar rcut, + Scalar epsilon, + Scalar sigma, + Scalar scale_factor, + bool energy_shift, + py::array_t + attraction_scale_factor_data, + py::array_t + attraction_scale_factor_shape, + py::array_t domain, + std::shared_ptr variant) + { + if (attraction_scale_factor_shape.size() != 2) + throw std::runtime_error("lambda_shape must have 2 elements"); + if (domain.size() != 4) + throw std::runtime_error( + "domain must have 4 elements [y_lo, y_hi, t_lo, t_hi]"); + + const unsigned int* shape_ptr = attraction_scale_factor_shape.data(); + const Scalar* data_ptr = attraction_scale_factor_data.data(); + const Scalar* dom_ptr = domain.data(); + + if (attraction_scale_factor_data.size() + != static_cast(shape_ptr[0] * shape_ptr[1])) + throw std::runtime_error( + "attraction_scale_factor_data size does not match lambda_shape"); + + PairParametersPerturbedLennardJonesEvap params(epsilon, sigma); + return std::make_shared(sysdef, + nlist, + rcut, + scale_factor, + params, + energy_shift, + data_ptr, + shape_ptr, + dom_ptr, + variant); + })) + .def_property_readonly("rcut", &PerturbedLennardJonesEvap::getRCut) + .def_property_readonly("epsilon", &PerturbedLennardJonesEvap::getEpsilon) + .def_property_readonly("sigma", &PerturbedLennardJonesEvap::getSigma); + } + + } // end namespace detail + } // end namespace azplugins + } // end namespace hoomd diff --git a/src/export_VariantInterpolated.cc b/src/export_VariantInterpolated.cc new file mode 100644 index 00000000..350705a6 --- /dev/null +++ b/src/export_VariantInterpolated.cc @@ -0,0 +1,51 @@ +// Copyright (c) 2018-2020, Michael P. Howard +// Copyright (c) 2021-2025, Auburn University +// Part of azplugins, released under the BSD 3-Clause License. + +#include "VariantInterpolated.h" + +#include +#include + +#include + +namespace hoomd + { +namespace azplugins + { +namespace detail + { +namespace py = pybind11; + +void export_VariantInterpolated(py::module& m) + { + py::class_>( + m, + "VariantInterpolated") + .def(py::init( + [](py::array_t data, + unsigned int n, + Scalar t_lo, + Scalar t_hi) + + { + if (data.ndim() < 1 || data.ndim() > 2) + throw std::runtime_error("data must be 1D or 2D (N, 1)."); + if (data.ndim() == 2 && data.shape(1) != 1) + throw std::runtime_error("data must have shape (N,) or (N, 1)."); + + const py::ssize_t N_ss = data.shape(0); + if (N_ss < 2) + throw std::runtime_error("data must contain at least 2 rows."); + + const unsigned int N = static_cast(N_ss); + const Scalar* data_ptr = data.data(); + + return std::make_shared(data_ptr, N, t_lo, t_hi); + })) + .def_property_readonly("t_lo", &VariantInterpolated::getTLo) + .def_property_readonly("t_hi", &VariantInterpolated::getTHi); + } + } // namespace detail + } // namespace azplugins + } // namespace hoomd diff --git a/src/module.cc b/src/module.cc index 1ab81c7f..4ffc9223 100644 --- a/src/module.cc +++ b/src/module.cc @@ -81,6 +81,9 @@ void export_PotentialPairDPDThermoGeneralWeight(pybind11::module&); void export_PotentialExternalWallColloid(pybind11::module&); void export_PotentialExternalWallLJ93(pybind11::module&); +void export_PerturbedLennardJonesEvap(pybind11::module&); +void export_VariantInterpolated(pybind11::module&); + #ifdef ENABLE_HIP // bond void export_ImagePotentialBondHarmonicGPU(pybind11::module&); @@ -153,6 +156,10 @@ PYBIND11_MODULE(_azplugins, m) export_PotentialExternalWallColloid(m); export_PotentialExternalWallLJ93(m); + // evaporation + export_PerturbedLennardJonesEvap(m); + export_VariantInterpolated(m); + #ifdef ENABLE_HIP // bond export_ImagePotentialBondHarmonicGPU(m); diff --git a/src/pair.py b/src/pair.py index fbefb07d..48c5b2e5 100644 --- a/src/pair.py +++ b/src/pair.py @@ -4,10 +4,13 @@ """Pair potentials.""" +import numpy + from hoomd.azplugins import _azplugins from hoomd.data.parameterdicts import ParameterDict, TypeParameterDict from hoomd.data.typeparam import TypeParameter from hoomd.md import pair +from hoomd.md.force import Force from hoomd.variant import Variant @@ -426,6 +429,67 @@ def __init__(self, nlist, default_r_cut=None, default_r_on=0, mode="none"): self._add_typeparam(params) +class PerturbedLennardJonesEvap(Force): + """Perturbed Lennard Jones potential for evaporation""" + + _ext_module = _azplugins + _cpp_class_name = "PerturbedLennardJonesEvap" + + def __init__( + self, + nlist, + rcut, + epsilon, + sigma, + scale_factor, + energy_shift, + attraction_scale_factor_data, + attraction_scale_factor_shape, + domain, + variant, + ): + super().__init__() + + self._nlist = nlist + self._variant = variant + + param_dict = ParameterDict( + epsilon=float, + sigma=float, + ) + + param_dict["epsilon"] = float(epsilon) + param_dict["sigma"] = float(sigma) + + self._param_dict.update(param_dict) + + self._domain = numpy.asarray(domain, dtype=numpy.float64) + self._attraction_scale_factor_data = numpy.asarray( + attraction_scale_factor_data, dtype=numpy.float64 + ) + self._attraction_scale_factor_shape = numpy.asarray( + attraction_scale_factor_shape, dtype=numpy.uint32 + ) + + def _attach_hook(self): + self._nlist._attach(self._simulation) + + cls = getattr(self._ext_module, self._cpp_class_name) + self._cpp_obj = cls( + self._simulation.state._cpp_sys_def, + self._nlist._cpp_obj, + self.rcut, + self.epsilon, + self.sigma, + self.scale_factor, + self.energy_shift, + self._attraction_scale_factor_data, + self._attraction_scale_factor_shape, + self._domain, + self._variant, + ) + + class TwoPatchMorse(pair.aniso.AnisotropicPair): r"""Two-patch Morse potential. diff --git a/src/variant.py b/src/variant.py new file mode 100644 index 00000000..9b80bd9b --- /dev/null +++ b/src/variant.py @@ -0,0 +1,28 @@ +# Copyright (c) 2018-2020, Michael P. Howard +# Copyright (c) 2021-2025, Auburn University +# Part of azplugins, released under the BSD 3-Clause License. + +"""Variant.""" + +import numpy +import hoomd +from hoomd.azplugins import _azplugins + + +class VariantInterpolated(hoomd.variant.Variant): + """Piecewise-linear variant on a uniform grid of time""" + + def __init__(self, data, n, t_lo, t_hi): + super().__init__() + + data = numpy.asarray(data, dtype=float) + if data.size < 2: + raise ValueError("data must contain at least 2 values.") + if t_hi <= t_lo: + raise ValueError("t_hi must be greater than t_lo.") + + self._data = data + + self._cpp_obj = _azplugins.VariantInterpolated( + self._data, int(n), float(t_lo), float(t_hi) + ) From 784cbbfa1ef9a5aef782cae872592888606b2beb Mon Sep 17 00:00:00 2001 From: mayukhkundu14 <79770572+mayukh33@users.noreply.github.com> Date: Mon, 8 Jun 2026 17:51:09 +0200 Subject: [PATCH 05/17] Add test for interpolated variant and other minor fixes --- src/LinearInterpolator1D.h | 2 +- src/PerturbedLennardJonesEvap.h | 2 +- src/VariantInterpolated.h | 17 ++++++--- src/export_VariantInterpolated.cc | 4 +- src/pytest/CMakeLists.txt | 1 + src/pytest/test_variant.py | 61 +++++++++++++++++++++++++++++++ src/variant.py | 13 +++---- 7 files changed, 84 insertions(+), 16 deletions(-) create mode 100644 src/pytest/test_variant.py diff --git a/src/LinearInterpolator1D.h b/src/LinearInterpolator1D.h index 5ed7dc7d..5445a2dc 100644 --- a/src/LinearInterpolator1D.h +++ b/src/LinearInterpolator1D.h @@ -41,7 +41,7 @@ template class LinearInterpolator1D m_dx = (m_hi - m_lo) / Scalar(n - 1); } - //! Implement piecewise linear interpolation + //! Implement piecewise linear interpolation on a uniform grid AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar operator()(Scalar x) const { const Scalar f = (x - m_lo) / m_dx; diff --git a/src/PerturbedLennardJonesEvap.h b/src/PerturbedLennardJonesEvap.h index b3e1ed0b..a2bc62ed 100644 --- a/src/PerturbedLennardJonesEvap.h +++ b/src/PerturbedLennardJonesEvap.h @@ -264,7 +264,7 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) pei += wca_shift * pair_eng; if (third_law) { - h_force.data[j].w += wca_shift * pair_eng; + h_force.data[j].w += Scalar(0.5) * wca_shift * pair_eng; } } else diff --git a/src/VariantInterpolated.h b/src/VariantInterpolated.h index 9c3455ec..be3bb5b1 100644 --- a/src/VariantInterpolated.h +++ b/src/VariantInterpolated.h @@ -2,11 +2,6 @@ // Copyright (c) 2021-2025, Auburn University // Part of azplugins, released under the BSD 3-Clause License. -/*! - * \file VariantInterpolated.h - * \brief Piecewise-linear variant on a uniform grid in scaled time, - * backed by LinearInterpolator1D. - */ #ifndef AZPLUGINS_VARIANT_INTERPOLATED_H_ #define AZPLUGINS_VARIANT_INTERPOLATED_H_ @@ -66,6 +61,18 @@ class PYBIND11_EXPORT VariantInterpolated : public Variant { return m_t_hi; } + + Scalar setTLo(Scalar t_lo) + { + m_t_lo = t_lo; + return m_t_lo; + } + + Scalar setTHi(Scalar t_hi) + { + m_t_hi = t_hi; + return m_t_hi; + } Scalar min() { diff --git a/src/export_VariantInterpolated.cc b/src/export_VariantInterpolated.cc index 350705a6..d792ebba 100644 --- a/src/export_VariantInterpolated.cc +++ b/src/export_VariantInterpolated.cc @@ -43,8 +43,8 @@ void export_VariantInterpolated(py::module& m) return std::make_shared(data_ptr, N, t_lo, t_hi); })) - .def_property_readonly("t_lo", &VariantInterpolated::getTLo) - .def_property_readonly("t_hi", &VariantInterpolated::getTHi); + .def_property("t_lo", &VariantInterpolated::getTLo, &VariantInterpolated::setTLo) + .def_property("t_hi", &VariantInterpolated::getTHi, &VariantInterpolated::setTHi); } } // namespace detail } // namespace azplugins diff --git a/src/pytest/CMakeLists.txt b/src/pytest/CMakeLists.txt index 157b326b..2002932c 100644 --- a/src/pytest/CMakeLists.txt +++ b/src/pytest/CMakeLists.txt @@ -9,6 +9,7 @@ set(test_files test_pair_aniso.py test_pair_dpd.py test_wall.py + test_variant.py ) # Copy tests to the install directory. diff --git a/src/pytest/test_variant.py b/src/pytest/test_variant.py new file mode 100644 index 00000000..b4a85971 --- /dev/null +++ b/src/pytest/test_variant.py @@ -0,0 +1,61 @@ +# Copyright (c) 2018-2020, Michael P. Howard +# Copyright (c) 2021-2025, Auburn University +# Part of azplugins, released under the BSD 3-Clause License. + +"""Pair potential unit tests.""" + +import collections + +import hoomd +import hoomd.azplugins +import numpy + +import pytest + +def interpolated_eval(data, t_lo, t_hi): + data = numpy.asarray(data, dtype=float) + return hoomd.azplugins.variant.VariantInterpolated( + data=data, t_lo=t_lo, t_hi=t_hi + ) + +variant_cases = [ + ([2.0, 4.0, 6.0, 8.0], 0, 300), + ([1.0, -2.0, 5.0, 0.5, 3.0], 10, 90), + ([-1.0, 1.0], 0, 100), +] + +@pytest.mark.parametrize("data, t_lo, t_hi", variant_cases) +def test_construction(data, t_lo, t_hi): + """t_lo and t_hi are stored and exposed as read-only properties.""" + variant = interpolated_eval(data, t_lo, t_hi) + assert variant.t_lo == pytest.approx(t_lo) + assert variant.t_hi == pytest.approx(t_hi) + + variant.t_lo = 5 + variant.t_hi = 5 + + assert variant.t_lo == pytest.approx(5) + assert variant.t_hi == pytest.approx(5) + +@pytest.mark.parametrize("data, t_lo, t_hi", variant_cases) +def test_interpolated_values(data, t_lo, t_hi): + time = numpy.linspace(t_lo, t_hi, len(data)) + for timestep in range(int(t_lo), int(t_hi) + 1): + expected = numpy.interp(timestep, time, data) + assert interpolated_eval(data, t_lo, t_hi)(timestep) == pytest.approx(expected) + + +@pytest.mark.parametrize("data, t_lo, t_hi", variant_cases) +def test_boundary_values(data, t_lo, t_hi): + """At t_lo the variant returns data[0]; at t_hi it returns data[-1].""" + variant = interpolated_eval(data, t_lo, t_hi) + assert variant(int(t_lo)) == pytest.approx(data[0]) + assert variant(int(t_hi)) == pytest.approx(data[-1]) + + +@pytest.mark.parametrize("data, t_lo, t_hi", variant_cases) +def test_min_max(data, t_lo, t_hi): + """min and max report the extrema of the data array.""" + variant = interpolated_eval(data, t_lo, t_hi) + assert variant.min == pytest.approx(numpy.min(data)) + assert variant.max == pytest.approx(numpy.max(data)) \ No newline at end of file diff --git a/src/variant.py b/src/variant.py index 9b80bd9b..f85e2869 100644 --- a/src/variant.py +++ b/src/variant.py @@ -9,11 +9,11 @@ from hoomd.azplugins import _azplugins -class VariantInterpolated(hoomd.variant.Variant): - """Piecewise-linear variant on a uniform grid of time""" +class VariantInterpolated(_azplugins.VariantInterpolated, hoomd.variant.Variant): + """Piecewise-linear variant on a uniform grid of time.""" - def __init__(self, data, n, t_lo, t_hi): - super().__init__() + def __init__(self, data, t_lo, t_hi): + hoomd.variant.Variant.__init__(self) data = numpy.asarray(data, dtype=float) if data.size < 2: @@ -22,7 +22,6 @@ def __init__(self, data, n, t_lo, t_hi): raise ValueError("t_hi must be greater than t_lo.") self._data = data - - self._cpp_obj = _azplugins.VariantInterpolated( - self._data, int(n), float(t_lo), float(t_hi) + _azplugins.VariantInterpolated.__init__( + self, data, data.size, float(t_lo), float(t_hi) ) From 855ee06008346b1c71d3f4b770993046bc13cc82 Mon Sep 17 00:00:00 2001 From: mayukhkundu14 <79770572+mayukh33@users.noreply.github.com> Date: Wed, 10 Jun 2026 11:41:37 +0200 Subject: [PATCH 06/17] Fix compute forces in PerturbedLennardJonesEvap and other minor fixes --- src/PerturbedLennardJonesEvap.h | 115 ++++++++++++++++++-------------- src/VariantInterpolated.h | 5 +- src/pair.py | 8 ++- src/pytest/test_variant.py | 54 +++++++-------- 4 files changed, 101 insertions(+), 81 deletions(-) diff --git a/src/PerturbedLennardJonesEvap.h b/src/PerturbedLennardJonesEvap.h index a2bc62ed..9c2c71e0 100644 --- a/src/PerturbedLennardJonesEvap.h +++ b/src/PerturbedLennardJonesEvap.h @@ -83,7 +83,12 @@ class PerturbedLennardJonesEvap : public ForceCompute const Scalar* domain, std::shared_ptr variant); - ~PerturbedLennardJonesEvap() { } + //! Destructor + ~PerturbedLennardJonesEvap() + { + if (m_r_cut_nlist) + m_nlist->removeRCutMatrix(m_r_cut_nlist); + } Scalar scaleTime(uint64_t timestep) const { @@ -121,6 +126,9 @@ class PerturbedLennardJonesEvap : public ForceCompute GPUArray m_attraction_scale_factor_shape; //!< [ny, nt] std::shared_ptr m_variant; + std::shared_ptr> + m_r_cut_nlist; //!< Cutoff matrix shared with the neighbor list + void computeForces(uint64_t timestep) override; }; @@ -160,7 +168,6 @@ PerturbedLennardJonesEvap::PerturbedLennardJonesEvap( std::copy(attraction_scale_factor_shape, attraction_scale_factor_shape + 2, h_shape.data); } - // Allocate and fill the tabulated data { const unsigned int n_data = attraction_scale_factor_shape[0] * attraction_scale_factor_shape[1]; @@ -174,6 +181,22 @@ PerturbedLennardJonesEvap::PerturbedLennardJonesEvap( std::copy(attraction_scale_factor_data, attraction_scale_factor_data + n_data, h_data.data); } + + { + const Index2D typpair_idx(m_pdata->getNTypes()); + m_r_cut_nlist + = std::make_shared>(typpair_idx.getNumElements(), m_exec_conf); + { + ArrayHandle h_r_cut_nlist(*m_r_cut_nlist, + access_location::host, + access_mode::overwrite); + + for (unsigned int i = 0; i < m_r_cut_nlist->getNumElements(); ++i) + h_r_cut_nlist.data[i] = m_rcut; + } + m_nlist->addRCutMatrix(m_r_cut_nlist); + m_nlist->notifyRCutMatrixChange(); + } } void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) @@ -189,7 +212,7 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::read); m_force.zeroFill(); - ArrayHandle h_force(m_force, access_location::host, access_mode::overwrite); + ArrayHandle h_force(m_force, access_location::host, access_mode::readwrite); ArrayHandle h_data(m_attraction_scale_factor_data, access_location::host, access_mode::read); @@ -214,6 +237,16 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) const Scalar hi[2] = {h_domain.data[1], h_domain.data[3]}; LinearInterpolator2D interp(h_data.data, h_shape.data, lo, hi); + auto clamp_scaled_y = [](Scalar y, Scalar height) + { + const Scalar s = y / height; + if (!(s > Scalar(0.0))) + return Scalar(0.0); + if (s > Scalar(1.0)) + return Scalar(1.0); + return s; + }; + const BoxDim box = m_pdata->getGlobalBox(); const unsigned int N = m_pdata->getN(); @@ -221,7 +254,7 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) { const Scalar3 pos_i = make_scalar3(h_pos.data[i].x, h_pos.data[i].y, h_pos.data[i].z); - Scalar scaled_y_i = Scalar(h_pos.data[i].y / interface_height); + Scalar scaled_y_i = clamp_scaled_y(h_pos.data[i].y, interface_height); Scalar3 fi = make_scalar3(0, 0, 0); Scalar pei = 0; @@ -235,7 +268,7 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) const unsigned int j = h_nlist.data[head + k]; if (j == i) continue; - Scalar scaled_y_j = Scalar(h_pos.data[j].y / interface_height); + Scalar scaled_y_j = clamp_scaled_y(h_pos.data[j].y, interface_height); Scalar3 pos_j = make_scalar3(h_pos.data[j].x, h_pos.data[j].y, h_pos.data[j].z); const Scalar attraction_scale_factor_j = interp(scaled_y_j, scaled_t); @@ -247,77 +280,61 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) const Scalar rsq = dot(dx, dx); const Scalar attraction_scale_factor_avg = Scalar(0.5) * (attraction_scale_factor_i + attraction_scale_factor_j); + const Scalar wca_shift - = epsilon_x_4 * (Scalar(1.0) - attraction_scale_factor_avg) / Scalar(4.); + = epsilon_x_4 * (Scalar(1.0) - attraction_scale_factor_avg) / Scalar(4.0); if (rsq < rcutsq && lj1 != 0) { const Scalar r2inv = Scalar(1) / rsq; const Scalar r6inv = r2inv * r2inv * r2inv; - const Scalar pair_eng = r6inv * (lj1 * r6inv - lj2); - const Scalar force_divr + Scalar pair_eng = r6inv * (lj1 * r6inv - lj2); + Scalar force_divr = r6inv * r2inv * (Scalar(12.0) * lj1 * r6inv - Scalar(6.0) * lj2); if (rsq < rwcasq) { - pei += wca_shift * pair_eng; - if (third_law) - { - h_force.data[j].w += Scalar(0.5) * wca_shift * pair_eng; - } + pair_eng += wca_shift; } else { - fi.x += attraction_scale_factor_avg * force_divr * dx.x; - fi.y += attraction_scale_factor_avg * force_divr * dx.y; - fi.z += attraction_scale_factor_avg * force_divr * dx.z; - pei += attraction_scale_factor_avg * pair_eng; - - if (third_law) - { - h_force.data[j].x -= attraction_scale_factor_avg * force_divr * dx.x; - h_force.data[j].y -= attraction_scale_factor_avg * force_divr * dx.y; - h_force.data[j].z -= attraction_scale_factor_avg * force_divr * dx.z; - h_force.data[j].w += Scalar(0.5) * attraction_scale_factor_avg * pair_eng; - } + pair_eng *= attraction_scale_factor_avg; + force_divr *= attraction_scale_factor_avg; } if (m_energy_shift) { const Scalar rcut2inv = Scalar(1.0) / rcutsq; const Scalar rcut6inv = rcut2inv * rcut2inv * rcut2inv; + Scalar pair_eng_shift = rcut6inv * (lj1 * rcut6inv - lj2); - const Scalar force_divr_shift - = rcut6inv * rcut2inv * (Scalar(12.0) * lj1 * rcut6inv - Scalar(6.0) * lj2); - if (rsq < rwcasq) + if (rcutsq < rwcasq) { - pei += wca_shift * pair_eng_shift; - if (third_law) - { - h_force.data[j].w += Scalar(0.5) * wca_shift * pair_eng_shift; - } + pair_eng_shift += wca_shift; } else { - fi.x += attraction_scale_factor_avg * force_divr_shift * dx.x; - fi.y += attraction_scale_factor_avg * force_divr_shift * dx.y; - fi.z += attraction_scale_factor_avg * force_divr_shift * dx.z; - pei += attraction_scale_factor_avg * pair_eng_shift; - - if (third_law) - { - h_force.data[j].x - -= attraction_scale_factor_avg * force_divr_shift * dx.x; - h_force.data[j].y - -= attraction_scale_factor_avg * force_divr_shift * dx.y; - h_force.data[j].z - -= attraction_scale_factor_avg * force_divr_shift * dx.z; - h_force.data[j].w - += Scalar(0.5) * attraction_scale_factor_avg * pair_eng_shift; - } + pair_eng_shift *= attraction_scale_factor_avg; } + + // Apply the shift to the pair energy + pair_eng -= pair_eng_shift; + } + + fi.x += force_divr * dx.x; + fi.y += force_divr * dx.y; + fi.z += force_divr * dx.z; + + pei += pair_eng; + + if (third_law) + { + h_force.data[j].x -= force_divr * dx.x; + h_force.data[j].y -= force_divr * dx.y; + h_force.data[j].z -= force_divr * dx.z; + h_force.data[j].w += Scalar(0.5) * pair_eng; } } } diff --git a/src/VariantInterpolated.h b/src/VariantInterpolated.h index be3bb5b1..32962b38 100644 --- a/src/VariantInterpolated.h +++ b/src/VariantInterpolated.h @@ -2,7 +2,6 @@ // Copyright (c) 2021-2025, Auburn University // Part of azplugins, released under the BSD 3-Clause License. - #ifndef AZPLUGINS_VARIANT_INTERPOLATED_H_ #define AZPLUGINS_VARIANT_INTERPOLATED_H_ @@ -61,13 +60,13 @@ class PYBIND11_EXPORT VariantInterpolated : public Variant { return m_t_hi; } - + Scalar setTLo(Scalar t_lo) { m_t_lo = t_lo; return m_t_lo; } - + Scalar setTHi(Scalar t_hi) { m_t_hi = t_hi; diff --git a/src/pair.py b/src/pair.py index 48c5b2e5..5942b7a1 100644 --- a/src/pair.py +++ b/src/pair.py @@ -454,13 +454,17 @@ def __init__( self._variant = variant param_dict = ParameterDict( + rcut=float, epsilon=float, sigma=float, + scale_factor=float, + energy_shift=bool, ) - + param_dict["rcut"] = float(rcut) param_dict["epsilon"] = float(epsilon) param_dict["sigma"] = float(sigma) - + param_dict["scale_factor"] = float(scale_factor) + param_dict["energy_shift"] = bool(energy_shift) self._param_dict.update(param_dict) self._domain = numpy.asarray(domain, dtype=numpy.float64) diff --git a/src/pytest/test_variant.py b/src/pytest/test_variant.py index b4a85971..0963811d 100644 --- a/src/pytest/test_variant.py +++ b/src/pytest/test_variant.py @@ -2,21 +2,17 @@ # Copyright (c) 2021-2025, Auburn University # Part of azplugins, released under the BSD 3-Clause License. -"""Pair potential unit tests.""" - -import collections - import hoomd import hoomd.azplugins import numpy import pytest + def interpolated_eval(data, t_lo, t_hi): data = numpy.asarray(data, dtype=float) - return hoomd.azplugins.variant.VariantInterpolated( - data=data, t_lo=t_lo, t_hi=t_hi - ) + return hoomd.azplugins.variant.VariantInterpolated(data=data, t_lo=t_lo, t_hi=t_hi) + variant_cases = [ ([2.0, 4.0, 6.0, 8.0], 0, 300), @@ -24,38 +20,42 @@ def interpolated_eval(data, t_lo, t_hi): ([-1.0, 1.0], 0, 100), ] -@pytest.mark.parametrize("data, t_lo, t_hi", variant_cases) -def test_construction(data, t_lo, t_hi): - """t_lo and t_hi are stored and exposed as read-only properties.""" + +@pytest.mark.parametrize("variant_case", variant_cases) +def test_construction(variant_case): + data, t_lo, t_hi = variant_case + variant = interpolated_eval(data, t_lo, t_hi) assert variant.t_lo == pytest.approx(t_lo) assert variant.t_hi == pytest.approx(t_hi) - + variant.t_lo = 5 variant.t_hi = 5 - + assert variant.t_lo == pytest.approx(5) - assert variant.t_hi == pytest.approx(5) + assert variant.t_hi == pytest.approx(5) + -@pytest.mark.parametrize("data, t_lo, t_hi", variant_cases) -def test_interpolated_values(data, t_lo, t_hi): +@pytest.mark.parametrize("variant_case", variant_cases) +def test_interpolated_values(variant_case): + data, t_lo, t_hi = variant_case time = numpy.linspace(t_lo, t_hi, len(data)) for timestep in range(int(t_lo), int(t_hi) + 1): expected = numpy.interp(timestep, time, data) assert interpolated_eval(data, t_lo, t_hi)(timestep) == pytest.approx(expected) - -@pytest.mark.parametrize("data, t_lo, t_hi", variant_cases) -def test_boundary_values(data, t_lo, t_hi): - """At t_lo the variant returns data[0]; at t_hi it returns data[-1].""" + +@pytest.mark.parametrize("variant_case", variant_cases) +def test_boundary_values(variant_case): + data, t_lo, t_hi = variant_case variant = interpolated_eval(data, t_lo, t_hi) - assert variant(int(t_lo)) == pytest.approx(data[0]) - assert variant(int(t_hi)) == pytest.approx(data[-1]) - - -@pytest.mark.parametrize("data, t_lo, t_hi", variant_cases) -def test_min_max(data, t_lo, t_hi): - """min and max report the extrema of the data array.""" + assert variant(t_lo) == pytest.approx(data[0]) + assert variant(t_hi) == pytest.approx(data[-1]) + + +@pytest.mark.parametrize("variant_case", variant_cases) +def test_min_max(variant_case): + data, t_lo, t_hi = variant_case variant = interpolated_eval(data, t_lo, t_hi) assert variant.min == pytest.approx(numpy.min(data)) - assert variant.max == pytest.approx(numpy.max(data)) \ No newline at end of file + assert variant.max == pytest.approx(numpy.max(data)) From 6b1fa442893d5b102025491a3b87c65c90ddb29f Mon Sep 17 00:00:00 2001 From: mayukhkundu14 <79770572+mayukh33@users.noreply.github.com> Date: Fri, 12 Jun 2026 12:17:49 +0200 Subject: [PATCH 07/17] Add initial unit tests for PerturbedLennardJonesEvap and other minor fixes --- src/LinearInterpolator1D.h | 1 + src/LinearInterpolator2D.h | 1 + src/PerturbedLennardJonesEvap.h | 15 +- src/export_PerturbedLennardJonesEvap.cc | 4 +- src/pair.py | 8 +- src/pytest/CMakeLists.txt | 1 + src/pytest/test_perturbedljevap.py | 244 ++++++++++++++++++++++++ 7 files changed, 261 insertions(+), 13 deletions(-) create mode 100644 src/pytest/test_perturbedljevap.py diff --git a/src/LinearInterpolator1D.h b/src/LinearInterpolator1D.h index 5445a2dc..64f059b4 100644 --- a/src/LinearInterpolator1D.h +++ b/src/LinearInterpolator1D.h @@ -55,6 +55,7 @@ template class LinearInterpolator1D assert(b >= 0); assert(b < static_cast(m_n - 1)); + const Scalar frac = f - Scalar(b); const Scalar c0 = m_data[b]; diff --git a/src/LinearInterpolator2D.h b/src/LinearInterpolator2D.h index e8de081a..937f452a 100644 --- a/src/LinearInterpolator2D.h +++ b/src/LinearInterpolator2D.h @@ -79,6 +79,7 @@ template class LinearInterpolator2D assert(b >= 0); assert(b < static_cast(nd - 1)); + bin[d] = b; frac[d] = f - Scalar(b); } diff --git a/src/PerturbedLennardJonesEvap.h b/src/PerturbedLennardJonesEvap.h index 9c2c71e0..76c2e631 100644 --- a/src/PerturbedLennardJonesEvap.h +++ b/src/PerturbedLennardJonesEvap.h @@ -74,8 +74,8 @@ class PerturbedLennardJonesEvap : public ForceCompute PerturbedLennardJonesEvap(std::shared_ptr sysdef, std::shared_ptr nlist, - Scalar r_cut, - Scalar scale_factor, + const Scalar r_cut, + const Scalar time_scale_factor, const param_type& params, bool energy_shift, const Scalar* attraction_scale_factor_data, @@ -92,7 +92,7 @@ class PerturbedLennardJonesEvap : public ForceCompute Scalar scaleTime(uint64_t timestep) const { - return Scalar(static_cast(timestep) / m_scale_factor); + return Scalar(static_cast(timestep) / m_time_scale_factor); } Scalar getRCut() const @@ -114,7 +114,7 @@ class PerturbedLennardJonesEvap : public ForceCompute std::shared_ptr m_nlist; //!< Neighbor list Scalar epsilon_x_4; Scalar m_rcut; - Scalar m_scale_factor; //!< Time scaling factor + Scalar m_time_scale_factor; //!< Time scaling factor Scalar lj1; Scalar lj2; Scalar rcutsq; @@ -135,8 +135,8 @@ class PerturbedLennardJonesEvap : public ForceCompute PerturbedLennardJonesEvap::PerturbedLennardJonesEvap( std::shared_ptr sysdef, std::shared_ptr nlist, - Scalar rcut, - Scalar scale_factor, + const Scalar rcut, + const Scalar time_scale_factor, const param_type& params, bool energy_shift, const Scalar* attraction_scale_factor_data, @@ -144,7 +144,8 @@ PerturbedLennardJonesEvap::PerturbedLennardJonesEvap( const Scalar* domain, std::shared_ptr variant) : ForceCompute(sysdef), m_nlist(nlist), epsilon_x_4(params.epsilon_x_4), m_rcut(rcut), - m_scale_factor(scale_factor), lj1(params.epsilon_x_4 * params.sigma_6 * params.sigma_6), + m_time_scale_factor(time_scale_factor), + lj1(params.epsilon_x_4 * params.sigma_6 * params.sigma_6), lj2(params.epsilon_x_4 * params.sigma_6), rcutsq(rcut * rcut), rwcasq(params.rwcasq), sigma_6(params.sigma_6), m_energy_shift(energy_shift), m_variant(variant) { diff --git a/src/export_PerturbedLennardJonesEvap.cc b/src/export_PerturbedLennardJonesEvap.cc index e47d84c7..6160187b 100644 --- a/src/export_PerturbedLennardJonesEvap.cc +++ b/src/export_PerturbedLennardJonesEvap.cc @@ -25,7 +25,7 @@ void export_PerturbedLennardJonesEvap(py::module& m) Scalar rcut, Scalar epsilon, Scalar sigma, - Scalar scale_factor, + Scalar time_scale_factor, bool energy_shift, py::array_t attraction_scale_factor_data, @@ -53,7 +53,7 @@ void export_PerturbedLennardJonesEvap(py::module& m) return std::make_shared(sysdef, nlist, rcut, - scale_factor, + time_scale_factor, params, energy_shift, data_ptr, diff --git a/src/pair.py b/src/pair.py index 5942b7a1..1fce2885 100644 --- a/src/pair.py +++ b/src/pair.py @@ -441,7 +441,7 @@ def __init__( rcut, epsilon, sigma, - scale_factor, + time_scale_factor, energy_shift, attraction_scale_factor_data, attraction_scale_factor_shape, @@ -457,13 +457,13 @@ def __init__( rcut=float, epsilon=float, sigma=float, - scale_factor=float, + time_scale_factor=float, energy_shift=bool, ) param_dict["rcut"] = float(rcut) param_dict["epsilon"] = float(epsilon) param_dict["sigma"] = float(sigma) - param_dict["scale_factor"] = float(scale_factor) + param_dict["time_scale_factor"] = float(time_scale_factor) param_dict["energy_shift"] = bool(energy_shift) self._param_dict.update(param_dict) @@ -485,7 +485,7 @@ def _attach_hook(self): self.rcut, self.epsilon, self.sigma, - self.scale_factor, + self.time_scale_factor, self.energy_shift, self._attraction_scale_factor_data, self._attraction_scale_factor_shape, diff --git a/src/pytest/CMakeLists.txt b/src/pytest/CMakeLists.txt index 2002932c..dc8630c3 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_perturbedljevap.py test_wall.py test_variant.py ) diff --git a/src/pytest/test_perturbedljevap.py b/src/pytest/test_perturbedljevap.py new file mode 100644 index 00000000..7cacc9ee --- /dev/null +++ b/src/pytest/test_perturbedljevap.py @@ -0,0 +1,244 @@ +# Copyright (c) 2018-2020, Michael P. Howard +# Copyright (c) 2021-2025, Auburn University +# Part of azplugins, released under the BSD 3-Clause License. + +import hoomd +import hoomd.azplugins +import numpy + +import pytest + + +@pytest.fixture +def valid_args_const(): + return { + "nlist": hoomd.md.nlist.Cell(buffer=0.4), + "rcut": 3.0, + "epsilon": 1.0, + "sigma": 1.0, + "time_scale_factor": 1.0, + "energy_shift": False, + "attraction_scale_factor_data": numpy.array([[0.6, 0.6], [0.6, 0.6]]), + "attraction_scale_factor_shape": [2, 2], + "domain": [0.0, 1.0, 0.0, 100.0], + "variant": hoomd.azplugins.variant.VariantInterpolated( + [5.0, 4.0, 2.0, 1.0], 0, 300 + ), + } + + +@pytest.fixture +def two_particle_snapshot_factory(): + snap = hoomd.Snapshot() + if snap.communicator.rank == 0: + snap.configuration.box = [20, 20, 20, 0, 0, 0] + snap.particles.N = 2 + snap.particles.types = ["A"] + snap.particles.position[:] = [[0, 0, 0], [0, 1.2, 0]] + + return snap + + +def test_constructor(valid_args_const): + evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**valid_args_const) + + assert evap.epsilon == 1.0 + assert evap.sigma == 1.0 + + assert evap._nlist is valid_args_const["nlist"] + assert evap._variant is valid_args_const["variant"] + + # Check numpy array casting + numpy.testing.assert_allclose(evap._domain, valid_args_const["domain"]) + numpy.testing.assert_allclose( + evap._attraction_scale_factor_data, + valid_args_const["attraction_scale_factor_data"], + ) + numpy.testing.assert_array_equal( + evap._attraction_scale_factor_shape, + valid_args_const["attraction_scale_factor_shape"], + ) + + assert evap.rcut == 3.0 + assert evap.time_scale_factor == 1.0 + + +def test_domain_mismatch(valid_args_const, two_particle_snapshot_factory): + snap = two_particle_snapshot_factory + bad_args = valid_args_const.copy() + bad_args["domain"] = [0.0, 10.0] + + sim = hoomd.Simulation(device=hoomd.device.CPU(), seed=42) + sim.create_state_from_snapshot(snap) + + integrator = hoomd.md.Integrator(dt=0.005) + sim.operations.integrator = integrator + + evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**bad_args) + sim.operations.integrator.forces = [evap] + + with pytest.raises(RuntimeError): + sim.run(0) + + +def test_shape_mismatch(valid_args_const, two_particle_snapshot_factory): + snap = two_particle_snapshot_factory + bad_args = valid_args_const.copy() + bad_args["attraction_scale_factor_shape"] = [2, 2, 2] + + sim = hoomd.Simulation(device=hoomd.device.CPU(), seed=42) + sim.create_state_from_snapshot(snap) + + integrator = hoomd.md.Integrator(dt=0.005) + sim.operations.integrator = integrator + + evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**bad_args) + sim.operations.integrator.forces = [evap] + + with pytest.raises(RuntimeError): + sim.run(0) + + +def test_data_size_mismatch(valid_args_const, two_particle_snapshot_factory): + snap = two_particle_snapshot_factory + bad_args = valid_args_const.copy() + bad_args["attraction_scale_factor_data"] = [1.0, 1.0] + + sim = hoomd.Simulation(device=hoomd.device.CPU(), seed=42) + sim.create_state_from_snapshot(snap) + + integrator = hoomd.md.Integrator(dt=0.005) + sim.operations.integrator = integrator + + evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**bad_args) + sim.operations.integrator.forces = [evap] + + with pytest.raises(RuntimeError): + sim.run(0) + + +def test_variant_mismatch(valid_args_const, two_particle_snapshot_factory): + snap = two_particle_snapshot_factory + bad_args = valid_args_const.copy() + bad_args["variant"] = hoomd.variant.Constant( + 1.0 + ) # Invalid: Should be VariantInterpolated + + sim = hoomd.Simulation(device=hoomd.device.CPU(), seed=42) + sim.create_state_from_snapshot(snap) + + integrator = hoomd.md.Integrator(dt=0.005) + sim.operations.integrator = integrator + + evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**bad_args) + sim.operations.integrator.forces = [evap] + + with pytest.raises(TypeError): + sim.run(0) + + +"""Test energy and force calculation for single particle type +with constant attraction scalefactor. +""" + + +def test_energy_and_force_calculation_const( + valid_args_const, two_particle_snapshot_factory +): + snap = two_particle_snapshot_factory + evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**valid_args_const) + + sim = hoomd.Simulation(device=hoomd.device.CPU(), seed=42) + sim.create_state_from_snapshot(snap) + + integrator = hoomd.md.Integrator(dt=0.005) + integrator.forces = [evap] + + sim.operations.integrator = integrator + + sim.run(0) + + expected_forces = [[0.0, 1.32701601, 0.0], [0.0, -1.32701601, 0.0]] + expected_energies = [-0.26728958627492283, -0.26728958627492283] + + if sim.device.communicator.rank == 0: + numpy.testing.assert_allclose(evap.forces, expected_forces) + numpy.testing.assert_allclose(evap.energies, expected_energies) + + +@pytest.fixture( + params=[ + # (time_scale_factor, domain) + (1.0, [0.0, 1.0, 0.0, 1000.0]), + (2.0, [0.0, 1.0, 0.0, 500.0]), + ], + ids=["unscaled_time", "scaled_time"], +) +def valid_args_vary(request): + + time_scale_factor, domain = request.param + attraction_factor_table = numpy.array( + [ + [0.7, 0.6, 0.5, 0.4, 0.3, 0.2], + [0.6, 0.5, 0.4, 0.3, 0.2, 0.1], + [0.5, 0.4, 0.3, 0.2, 0.1, 0.0], + [0.4, 0.3, 0.2, 0.1, 0.0, 0.0], + ] + ) + + return { + "nlist": hoomd.md.nlist.Cell(buffer=0.4), + "rcut": 3.0, + "epsilon": 1.0, + "sigma": 1.0, + "time_scale_factor": time_scale_factor, + "energy_shift": False, + "attraction_scale_factor_data": attraction_factor_table, + "attraction_scale_factor_shape": attraction_factor_table.shape, + "domain": domain, + "variant": hoomd.azplugins.variant.VariantInterpolated( + [12, 10, 8, 6, 4, 2], + 0.0, + 1000.0, # Unscaled coordinates for the interpolated variant + ), + } + + +"""Test energy and force calculation for single particle type +with varying attraction scale factor. +""" + + +def test_energy_and_force_calculation_vary( + valid_args_vary, two_particle_snapshot_factory +): + """Energies/forces at t=0 and t=2.5, with and without time scaling.""" + snap = two_particle_snapshot_factory + evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**valid_args_vary) + + sim = hoomd.Simulation(device=hoomd.device.CPU(), seed=42) + sim.create_state_from_snapshot(snap) + + integrator = hoomd.md.Integrator(dt=0.005) + integrator.forces = [evap] + sim.operations.integrator = integrator + sim.run(0) + + expected_forces = [[0, 1.51500993942280826, 0], [0, -1.51500993942280826, 0]] + + expected_energies = [-0.305155610997203564, -0.305155610997203564] + + if sim.device.communicator.rank == 0: + numpy.testing.assert_allclose(evap.forces, expected_forces) + numpy.testing.assert_allclose(evap.energies, expected_energies) + + """Test if the potential energy and forces change as expected after a certain time. + """ + sim.run(500) # Run for 500 steps, which corresponds to t = 2.5 for dt = 0.005 + + expected_energies = [-0.18901192172298115, -0.18901192172298115] + expected_forces = [[0.0, 0.9383898894860773, 0.0], [0.0, -0.9383898894860773, 0.0]] + + if sim.device.communicator.rank == 0: + numpy.testing.assert_allclose(evap.forces, expected_forces) + numpy.testing.assert_allclose(evap.energies, expected_energies) From a25feb1b7ae564e6d8efb74769830afa1ca43784 Mon Sep 17 00:00:00 2001 From: mayukhkundu14 <79770572+mayukh33@users.noreply.github.com> Date: Mon, 15 Jun 2026 10:25:04 +0200 Subject: [PATCH 08/17] Debug perturbedljevap unit test --- src/pytest/test_perturbedljevap.py | 121 +++++++++++++++++------------ 1 file changed, 71 insertions(+), 50 deletions(-) diff --git a/src/pytest/test_perturbedljevap.py b/src/pytest/test_perturbedljevap.py index 7cacc9ee..7c6550bb 100644 --- a/src/pytest/test_perturbedljevap.py +++ b/src/pytest/test_perturbedljevap.py @@ -8,6 +8,44 @@ import pytest +_DEVICE_PARAMS = ["cpu"] + +if hoomd.version.gpu_enabled: + try: + if len(hoomd.device.GPU.get_available_devices()) > 0: + _DEVICE_PARAMS.append("gpu") + except Exception: + pass + + +@pytest.fixture(params=_DEVICE_PARAMS) +def simulation_factory(request): + """Create a Simulation on CPU, and on GPU when available.""" + + def make_simulation(snapshot): + if request.param == "cpu": + device = hoomd.device.CPU() + else: + device = hoomd.device.GPU() + + sim = hoomd.Simulation(device=device, seed=1) + sim.create_state_from_snapshot(snapshot) + return sim + + return make_simulation + + +@pytest.fixture +def two_particle_snapshot_factory(): + snap = hoomd.Snapshot() + if snap.communicator.rank == 0: + snap.configuration.box = [20, 20, 20, 0, 0, 0] + snap.particles.N = 2 + snap.particles.types = ["A"] + snap.particles.position[:] = [[0, 0, 0], [0, 1.2, 0]] + + return snap + @pytest.fixture def valid_args_const(): @@ -27,18 +65,6 @@ def valid_args_const(): } -@pytest.fixture -def two_particle_snapshot_factory(): - snap = hoomd.Snapshot() - if snap.communicator.rank == 0: - snap.configuration.box = [20, 20, 20, 0, 0, 0] - snap.particles.N = 2 - snap.particles.types = ["A"] - snap.particles.position[:] = [[0, 0, 0], [0, 1.2, 0]] - - return snap - - def test_constructor(valid_args_const): evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**valid_args_const) @@ -63,14 +89,15 @@ def test_constructor(valid_args_const): assert evap.time_scale_factor == 1.0 -def test_domain_mismatch(valid_args_const, two_particle_snapshot_factory): +def test_domain_mismatch( + valid_args_const, two_particle_snapshot_factory, simulation_factory +): snap = two_particle_snapshot_factory + sim = simulation_factory(snap) + bad_args = valid_args_const.copy() bad_args["domain"] = [0.0, 10.0] - sim = hoomd.Simulation(device=hoomd.device.CPU(), seed=42) - sim.create_state_from_snapshot(snap) - integrator = hoomd.md.Integrator(dt=0.005) sim.operations.integrator = integrator @@ -81,13 +108,14 @@ def test_domain_mismatch(valid_args_const, two_particle_snapshot_factory): sim.run(0) -def test_shape_mismatch(valid_args_const, two_particle_snapshot_factory): +def test_shape_mismatch( + valid_args_const, two_particle_snapshot_factory, simulation_factory +): snap = two_particle_snapshot_factory bad_args = valid_args_const.copy() bad_args["attraction_scale_factor_shape"] = [2, 2, 2] - sim = hoomd.Simulation(device=hoomd.device.CPU(), seed=42) - sim.create_state_from_snapshot(snap) + sim = simulation_factory(snap) integrator = hoomd.md.Integrator(dt=0.005) sim.operations.integrator = integrator @@ -99,13 +127,14 @@ def test_shape_mismatch(valid_args_const, two_particle_snapshot_factory): sim.run(0) -def test_data_size_mismatch(valid_args_const, two_particle_snapshot_factory): +def test_data_size_mismatch( + valid_args_const, two_particle_snapshot_factory, simulation_factory +): snap = two_particle_snapshot_factory bad_args = valid_args_const.copy() bad_args["attraction_scale_factor_data"] = [1.0, 1.0] - sim = hoomd.Simulation(device=hoomd.device.CPU(), seed=42) - sim.create_state_from_snapshot(snap) + sim = simulation_factory(snap) integrator = hoomd.md.Integrator(dt=0.005) sim.operations.integrator = integrator @@ -117,15 +146,16 @@ def test_data_size_mismatch(valid_args_const, two_particle_snapshot_factory): sim.run(0) -def test_variant_mismatch(valid_args_const, two_particle_snapshot_factory): +def test_variant_mismatch( + valid_args_const, two_particle_snapshot_factory, simulation_factory +): snap = two_particle_snapshot_factory bad_args = valid_args_const.copy() bad_args["variant"] = hoomd.variant.Constant( 1.0 ) # Invalid: Should be VariantInterpolated - sim = hoomd.Simulation(device=hoomd.device.CPU(), seed=42) - sim.create_state_from_snapshot(snap) + sim = simulation_factory(snap) integrator = hoomd.md.Integrator(dt=0.005) sim.operations.integrator = integrator @@ -143,13 +173,12 @@ def test_variant_mismatch(valid_args_const, two_particle_snapshot_factory): def test_energy_and_force_calculation_const( - valid_args_const, two_particle_snapshot_factory + valid_args_const, two_particle_snapshot_factory, simulation_factory ): snap = two_particle_snapshot_factory evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**valid_args_const) - sim = hoomd.Simulation(device=hoomd.device.CPU(), seed=42) - sim.create_state_from_snapshot(snap) + sim = simulation_factory(snap) integrator = hoomd.md.Integrator(dt=0.005) integrator.forces = [evap] @@ -161,16 +190,19 @@ def test_energy_and_force_calculation_const( expected_forces = [[0.0, 1.32701601, 0.0], [0.0, -1.32701601, 0.0]] expected_energies = [-0.26728958627492283, -0.26728958627492283] + forces = evap.forces + energies = evap.energies + if sim.device.communicator.rank == 0: - numpy.testing.assert_allclose(evap.forces, expected_forces) - numpy.testing.assert_allclose(evap.energies, expected_energies) + numpy.testing.assert_allclose(forces, expected_forces) + numpy.testing.assert_allclose(energies, expected_energies) @pytest.fixture( params=[ # (time_scale_factor, domain) - (1.0, [0.0, 1.0, 0.0, 1000.0]), - (2.0, [0.0, 1.0, 0.0, 500.0]), + (1.0, [0.0, 1.0, 0.0, 100.0]), + (2.0, [0.0, 1.0, 0.0, 50.0]), ], ids=["unscaled_time", "scaled_time"], ) @@ -199,7 +231,7 @@ def valid_args_vary(request): "variant": hoomd.azplugins.variant.VariantInterpolated( [12, 10, 8, 6, 4, 2], 0.0, - 1000.0, # Unscaled coordinates for the interpolated variant + 100.0, # Run for 100 timesteps, which corresponds to t = 0.5 for dt = 0.005 ), } @@ -210,14 +242,12 @@ def valid_args_vary(request): def test_energy_and_force_calculation_vary( - valid_args_vary, two_particle_snapshot_factory + valid_args_vary, two_particle_snapshot_factory, simulation_factory ): """Energies/forces at t=0 and t=2.5, with and without time scaling.""" snap = two_particle_snapshot_factory evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**valid_args_vary) - - sim = hoomd.Simulation(device=hoomd.device.CPU(), seed=42) - sim.create_state_from_snapshot(snap) + sim = simulation_factory(snap) integrator = hoomd.md.Integrator(dt=0.005) integrator.forces = [evap] @@ -228,17 +258,8 @@ def test_energy_and_force_calculation_vary( expected_energies = [-0.305155610997203564, -0.305155610997203564] + forces = evap.forces + energies = evap.energies if sim.device.communicator.rank == 0: - numpy.testing.assert_allclose(evap.forces, expected_forces) - numpy.testing.assert_allclose(evap.energies, expected_energies) - - """Test if the potential energy and forces change as expected after a certain time. - """ - sim.run(500) # Run for 500 steps, which corresponds to t = 2.5 for dt = 0.005 - - expected_energies = [-0.18901192172298115, -0.18901192172298115] - expected_forces = [[0.0, 0.9383898894860773, 0.0], [0.0, -0.9383898894860773, 0.0]] - - if sim.device.communicator.rank == 0: - numpy.testing.assert_allclose(evap.forces, expected_forces) - numpy.testing.assert_allclose(evap.energies, expected_energies) + numpy.testing.assert_allclose(forces, expected_forces) + numpy.testing.assert_allclose(energies, expected_energies) From 7b9a8491025f0d495e368241ec508ca7e8e20c86 Mon Sep 17 00:00:00 2001 From: mayukhkundu14 <79770572+mayukh33@users.noreply.github.com> Date: Mon, 15 Jun 2026 10:55:13 +0200 Subject: [PATCH 09/17] Add back time dependent unit test --- src/pytest/test_perturbedljevap.py | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/src/pytest/test_perturbedljevap.py b/src/pytest/test_perturbedljevap.py index 7c6550bb..41b3a6d9 100644 --- a/src/pytest/test_perturbedljevap.py +++ b/src/pytest/test_perturbedljevap.py @@ -2,6 +2,8 @@ # Copyright (c) 2021-2025, Auburn University # Part of azplugins, released under the BSD 3-Clause License. +"""Test Perturbed Lennard-Jones for Evaporation pair potential.""" + import hoomd import hoomd.azplugins import numpy @@ -244,7 +246,7 @@ def valid_args_vary(request): def test_energy_and_force_calculation_vary( valid_args_vary, two_particle_snapshot_factory, simulation_factory ): - """Energies/forces at t=0 and t=2.5, with and without time scaling.""" + """Energies/forces at t=0 and t=0.5, with and without time scaling.""" snap = two_particle_snapshot_factory evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**valid_args_vary) sim = simulation_factory(snap) @@ -263,3 +265,26 @@ def test_energy_and_force_calculation_vary( if sim.device.communicator.rank == 0: numpy.testing.assert_allclose(forces, expected_forces) numpy.testing.assert_allclose(energies, expected_energies) + + """Test if the potential energy and forces change as expected after a certain time. + """ + + sim.run(100) + + snapshot = sim.state.get_snapshot() + positions = snapshot.particles.position + + print(f"Particle 1 position: {positions[1]}") + + expected_energies = [-0.04900309081706919, -0.04900309081706919] + expected_forces = [ + [0.0, 0.24328626764453856, 0.0], + [0.0, -0.24328626764453856, 0.0], + ] + + forces = evap.forces + energies = evap.energies + + if sim.device.communicator.rank == 0: + numpy.testing.assert_allclose(forces, expected_forces) + numpy.testing.assert_allclose(energies, expected_energies) From dd6924936706766fe5ed4ddf5a5aa350e3656b5e Mon Sep 17 00:00:00 2001 From: mayukhkundu14 <79770572+mayukh33@users.noreply.github.com> Date: Mon, 15 Jun 2026 14:08:27 +0200 Subject: [PATCH 10/17] Fix perturbedlj unit test --- src/pytest/test_perturbedljevap.py | 56 +++++++----------------------- 1 file changed, 13 insertions(+), 43 deletions(-) diff --git a/src/pytest/test_perturbedljevap.py b/src/pytest/test_perturbedljevap.py index 41b3a6d9..d0daffe5 100644 --- a/src/pytest/test_perturbedljevap.py +++ b/src/pytest/test_perturbedljevap.py @@ -10,41 +10,16 @@ import pytest -_DEVICE_PARAMS = ["cpu"] - -if hoomd.version.gpu_enabled: - try: - if len(hoomd.device.GPU.get_available_devices()) > 0: - _DEVICE_PARAMS.append("gpu") - except Exception: - pass - - -@pytest.fixture(params=_DEVICE_PARAMS) -def simulation_factory(request): - """Create a Simulation on CPU, and on GPU when available.""" - - def make_simulation(snapshot): - if request.param == "cpu": - device = hoomd.device.CPU() - else: - device = hoomd.device.GPU() - - sim = hoomd.Simulation(device=device, seed=1) - sim.create_state_from_snapshot(snapshot) - return sim - - return make_simulation - @pytest.fixture def two_particle_snapshot_factory(): snap = hoomd.Snapshot() + d = 1 if snap.communicator.rank == 0: snap.configuration.box = [20, 20, 20, 0, 0, 0] snap.particles.N = 2 snap.particles.types = ["A"] - snap.particles.position[:] = [[0, 0, 0], [0, 1.2, 0]] + snap.particles.position[:] = [[-d/2, 0, 0], [d/2, 1.2, 0]] return snap @@ -189,8 +164,8 @@ def test_energy_and_force_calculation_const( sim.run(0) - expected_forces = [[0.0, 1.32701601, 0.0], [0.0, -1.32701601, 0.0]] - expected_energies = [-0.26728958627492283, -0.26728958627492283] + expected_forces = [[0.35032675, 0.4203921, 0.], [-0.35032675, -0.4203921, 0.]] + expected_energies = [-0.07691956918309711, -0.07691956918309711] forces = evap.forces energies = evap.energies @@ -256,9 +231,10 @@ def test_energy_and_force_calculation_vary( sim.operations.integrator = integrator sim.run(0) - expected_forces = [[0, 1.51500993942280826, 0], [0, -1.51500993942280826, 0]] - - expected_energies = [-0.305155610997203564, -0.305155610997203564] + expected_forces = [[0.39995636970934184, 0.47994764365121023, 0.0], + [-0.39995636970934184, -0.47994764365121023, 0.0]] + + expected_energies = [-0.08781650815070254, -0.08781650815070254] forces = evap.forces energies = evap.energies @@ -270,17 +246,10 @@ def test_energy_and_force_calculation_vary( """ sim.run(100) - - snapshot = sim.state.get_snapshot() - positions = snapshot.particles.position - - print(f"Particle 1 position: {positions[1]}") - - expected_energies = [-0.04900309081706919, -0.04900309081706919] - expected_forces = [ - [0.0, 0.24328626764453856, 0.0], - [0.0, -0.24328626764453856, 0.0], - ] + + expected_energies = [-0.014101921016901138, -0.014101921016901138] + expected_forces = [[0.06422657031828848, 0.07707188438194618, 0.0], + [-0.06422657031828848, -0.07707188438194618, 0.0]] forces = evap.forces energies = evap.energies @@ -288,3 +257,4 @@ def test_energy_and_force_calculation_vary( if sim.device.communicator.rank == 0: numpy.testing.assert_allclose(forces, expected_forces) numpy.testing.assert_allclose(energies, expected_energies) + From 02a6d9dcf3760ff1afdb59e620854a098e384f6a Mon Sep 17 00:00:00 2001 From: mayukhkundu14 <79770572+mayukh33@users.noreply.github.com> Date: Mon, 15 Jun 2026 15:19:49 +0200 Subject: [PATCH 11/17] Initial GPU implementation --- src/CMakeLists.txt | 1 + src/PerturbedLennardJonesEvapGPU.cuh | 254 +++++++++++++++++++++ src/PerturbedLennardJonesEvapGPU.h | 124 ++++++++++ src/export_PerturbedLennardJonesEvapGPU.cc | 70 ++++++ src/module.cc | 2 + src/pytest/test_perturbedljevap.py | 23 +- src/pytest/test_variant.py | 2 + 7 files changed, 466 insertions(+), 10 deletions(-) create mode 100644 src/PerturbedLennardJonesEvapGPU.cuh create mode 100644 src/PerturbedLennardJonesEvapGPU.h create mode 100644 src/export_PerturbedLennardJonesEvapGPU.cc diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 7f0ea7eb..b201a7c7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -14,6 +14,7 @@ set(_${COMPONENT_NAME}_sources if (ENABLE_HIP) list(APPEND _${COMPONENT_NAME}_sources export_ImagePotentialBondHarmonicGPU.cc + export_PerturbedLennardJonesEvapGPU.cc VelocityComputeGPU.cc ) endif() diff --git a/src/PerturbedLennardJonesEvapGPU.cuh b/src/PerturbedLennardJonesEvapGPU.cuh new file mode 100644 index 00000000..ffcd182e --- /dev/null +++ b/src/PerturbedLennardJonesEvapGPU.cuh @@ -0,0 +1,254 @@ +// 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_PERTURBED_LENNARD_JONES_EVAP_GPU_CUH_ +#define AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_GPU_CUH_ + +#include "hip/hip_runtime.h" +#include "hoomd/BoxDim.h" +#include "hoomd/HOOMDMath.h" +#include "hoomd/Index1D.h" +#include "hoomd/ParticleData.cuh" +#include "hoomd/TextureTools.h" + +#include "LinearInterpolator2D.h" + +#include + +namespace hoomd + { +namespace azplugins + { +namespace gpu + { + +//! Wraps arguments to kernel driver +struct compute_perturbed_lennard_jones_evap_args_t + { + //! Construct a compute_perturbed_lennard_jones_evap_args_t + compute_perturbed_lennard_jones_evap_args_t(Scalar4* _d_force, + const Scalar4* _d_pos, + const BoxDim& _box, + const unsigned int _N, + const unsigned int* _d_n_neigh, + const unsigned int* _d_nlist, + const size_t* _d_head_list, + const Scalar* _d_attraction_scale_factor_data, + const uint2 _attraction_scale_factor_shape, + const Scalar4 _domain, + const Scalar _scaled_t, + const Scalar _interface_height, + const Scalar _lj1, + const Scalar _lj2, + const Scalar _epsilon_x_4, + const Scalar _rcutsq, + const Scalar _rwcasq, + const bool _energy_shift, + const unsigned int _block_size) + : d_force(_d_force), d_pos(_d_pos), box(_box), d_n_neigh(_d_n_neigh), d_nlist(_d_nlist), + d_head_list(_d_head_list), + d_attraction_scale_factor_data(_d_attraction_scale_factor_data), + attraction_scale_factor_shape(_attraction_scale_factor_shape), domain(_domain), + scaled_t(_scaled_t), interface_height(_interface_height), lj1(_lj1), lj2(_lj2), + epsilon_x_4(_epsilon_x_4), rcutsq(_rcutsq), rwcasq(_rwcasq), energy_shift(_energy_shift), + block_size(_block_size) { }; + + Scalar4* d_force; //!< Force to write out + const Scalar4* d_pos; //!< Particle positions + const BoxDim box; //!< Box dimensions + const unsigned int N; //!< Number of particles + const unsigned int* d_n_neigh; //!< Number of neighbors per particle + const unsigned int* d_nlist; //!< Neighbor list + const size_t* d_head_list; //!< Head list indexing into the neighbor list + const Scalar* d_attraction_scale_factor_data; //!< Flattened (y, t) table data + const uint2 attraction_scale_factor_shape; //!< Table shape (ny, nt) + const Scalar4 domain; //!< Table domain [y_lo, y_hi, t_lo, t_hi] + const Scalar scaled_t; //!< Scaled simulation time + const Scalar interface_height; //!< Interface height at this timestep + const Scalar lj1; //!< 4 epsilon sigma^12 + const Scalar lj2; //!< 4 epsilon sigma^6 + const Scalar epsilon_x_4; //!< 4 epsilon + const Scalar rcutsq; //!< Squared cutoff radius + const Scalar rwcasq; //!< Squared WCA (potential minimum) radius + const bool energy_shift; //!< Whether to shift the energy at the cutoff + const unsigned int block_size; //!< Block size to execute + }; + +hipError_t compute_perturbed_lennard_jones_evap_forces( + const compute_perturbed_lennard_jones_evap_args_t& args); + +#ifdef __HIPCC__ + +namespace kernel + { + +__device__ __forceinline__ Scalar clamp_scaled_y(Scalar y, Scalar height) + { + const Scalar s = y / height; + if (!(s > Scalar(0.0))) + return Scalar(0.0); + if (s > Scalar(1.0)) + return Scalar(1.0); + return s; + } + +__global__ void +compute_perturbed_lennard_jones_evap_forces(Scalar4* d_force, + const Scalar4* d_pos, + const BoxDim& box, + const unsigned int N, + const unsigned int* d_n_neigh, + const unsigned int* d_nlist, + const size_t* d_head_list, + const LinearInterpolator2D& interp, + const Scalar interface_height, + const Scalar scaled_t, + const Scalar lj1, + const Scalar lj2, + const Scalar epsilon_x_4, + const Scalar rcutsq, + const Scalar rwcasq, + const bool energy_shift) + { + // start by identifying which particle we are to handle + const unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx >= N) + return; + + const Scalar4 postype_i = __ldg(d_pos + idx); + const Scalar3 pos_i = make_scalar3(postype_i.x, postype_i.y, postype_i.z); + + const Scalar scaled_y_i = clamp_scaled_y(postype_i.y, interface_height); + const Scalar attraction_scale_factor_i = interp(scaled_y_i, scaled_t); + + // initialize the force and energy to 0 + Scalar3 fi = make_scalar3(0, 0, 0); + Scalar pei(0); + + // load in the length of the neighbor list for this thread + const unsigned int n_neigh = d_n_neigh[idx]; + const size_t head = d_head_list[idx]; + + // loop over neighbors + for (unsigned int k = 0; k < n_neigh; ++k) + { + const unsigned int j = d_nlist[head + k]; + if (j == idx) + continue; + + const Scalar4 postype_j = __ldg(d_pos + j); + const Scalar3 pos_j = make_scalar3(postype_j.x, postype_j.y, postype_j.z); + + const Scalar scaled_y_j = clamp_scaled_y(postype_j.y, interface_height); + const Scalar attraction_scale_factor_j = interp(scaled_y_j, scaled_t); + + // minimum-image + Scalar3 dx = pos_i - pos_j; + dx = box.minImage(dx); + + const Scalar rsq = dot(dx, dx); + const Scalar attraction_scale_factor_avg + = Scalar(0.5) * (attraction_scale_factor_i + attraction_scale_factor_j); + + const Scalar wca_shift + = epsilon_x_4 * (Scalar(1.0) - attraction_scale_factor_avg) / Scalar(4.0); + + if (rsq < rcutsq && lj1 != 0) + { + const Scalar r2inv = Scalar(1) / rsq; + const Scalar r6inv = r2inv * r2inv * r2inv; + + Scalar pair_eng = r6inv * (lj1 * r6inv - lj2); + Scalar force_divr = r6inv * r2inv * (Scalar(12.0) * lj1 * r6inv - Scalar(6.0) * lj2); + + if (rsq < rwcasq) + { + pair_eng += wca_shift; + } + else + { + pair_eng *= attraction_scale_factor_avg; + force_divr *= attraction_scale_factor_avg; + } + + if (energy_shift) + { + const Scalar rcut2inv = Scalar(1.0) / rcutsq; + const Scalar rcut6inv = rcut2inv * rcut2inv * rcut2inv; + + Scalar pair_eng_shift = rcut6inv * (lj1 * rcut6inv - lj2); + + if (rcutsq < rwcasq) + { + pair_eng_shift += wca_shift; + } + else + { + pair_eng_shift *= attraction_scale_factor_avg; + } + + // apply the shift to the pair energy + pair_eng -= pair_eng_shift; + } + + fi.x += force_divr * dx.x; + fi.y += force_divr * dx.y; + fi.z += force_divr * dx.z; + + pei += pair_eng; + } + } + + d_force[idx] = make_scalar4(fi.x, fi.y, fi.z, Scalar(0.5) * pei); + } + + } // end namespace kernel + +hipError_t +compute_perturbed_lennard_jones_evap_forces(const compute_perturbed_lennard_jones_evap_args_t& args) + { + unsigned int max_block_size; + hipFuncAttributes attr; + hipFuncGetAttributes( + &attr, + reinterpret_cast(&kernel::compute_perturbed_lennard_jones_evap_forces)); + max_block_size = attr.maxThreadsPerBlock; + + const unsigned int run_block_size = min(args.block_size, max_block_size); + + // setup the grid to run the kernel + dim3 grid(args.N / run_block_size + 1, 1, 1); + dim3 threads(run_block_size, 1, 1); + + // run the kernel + hipLaunchKernelGGL((kernel::compute_perturbed_lennard_jones_evap_forces), + grid, + threads, + args.d_force, + args.d_pos, + args.N, + args.d_n_neigh, + args.d_nlist, + args.d_head_list, + args.interp, + args.scaled_t, + args.interface_height, + args.lj1, + args.lj2, + args.epsilon_x_4, + args.rcutsq, + args.rwcasq, + args.energy_shift); + + return hipSuccess; + } +#endif // __HIPCC__ + + } // end namespace gpu + } // end namespace azplugins + } // end namespace hoomd + +#endif // AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_GPU_CUH_ diff --git a/src/PerturbedLennardJonesEvapGPU.h b/src/PerturbedLennardJonesEvapGPU.h new file mode 100644 index 00000000..bf6c335c --- /dev/null +++ b/src/PerturbedLennardJonesEvapGPU.h @@ -0,0 +1,124 @@ +// Copyright (c) 2018-2020, Michael P. Howard +// Copyright (c) 2021-2025, Auburn University +// Part of azplugins, released under the BSD 3-Clause License. + +#include "PerturbedLennardJonesEvapGPU.cuh" +#include "PerturbedLennardJonesEvapGPU.h" + +namespace hoomd + { +namespace azplugins + { +class PYBIND11_EXPORT PerturbedLennardJonesEvapGPU : public PerturbedLennardJonesEvap + { + public: + //! Constructor + PerturbedLennardJonesEvapGPU(std::shared_ptr sysdef, + std::shared_ptr nlist, + const Scalar r_cut, + const Scalar scale_factor, + const param_type& params, + bool energy_shift, + const Scalar* attraction_scale_factor_data, + const unsigned int* attraction_scale_factor_shape, + const Scalar* domain, + std::shared_ptr variant); + + protected: + std::shared_ptr> m_tuner; + void computeForces(uint64_t timestep) override; + }; +PerturbedLennardJonesEvapGPU::PerturbedLennardJonesEvapGPU( + std::shared_ptr sysdef, + std::shared_ptr nlist, + const Scalar r_cut, + const Scalar scale_factor, + const param_type& params, + bool energy_shift, + const Scalar* attraction_scale_factor_data, + const unsigned int* attraction_scale_factor_shape, + const Scalar* domain, + std::shared_ptr variant) + : PerturbedLennardJonesEvap(sysdef, + nlist, + r_cut, + scale_factor, + params, + energy_shift, + attraction_scale_factor_data, + attraction_scale_factor_shape, + domain, + variant) + { + m_tuner.reset(new Autotuner<1>({AutotunerBase::makeBlockSizeRange(m_exec_conf)}, + m_exec_conf, + "perturbed_lennard_jones_evap")); + m_autotuners.push_back(m_tuner); + } + +void PerturbedLennardJonesEvapGPU::computeForces(uint64_t timestep) + { + m_nlist->compute(timestep); + + const Scalar interface_height = (*m_variant)(timestep); + const Scalar scaled_t = scaleTime(timestep); + + ArrayHandle d_pos(m_pdata->getPositions(), access_location::device, access_mode::read); + + // neighbor-list arrays + ArrayHandle d_n_neigh(m_nlist->getNNeighArray(), + access_location::device, + access_mode::read); + ArrayHandle d_nlist(m_nlist->getNListArray(), + access_location::device, + access_mode::read); + ArrayHandle d_head_list(m_nlist->getHeadList(), + access_location::device, + access_mode::read); + + ArrayHandle d_data(m_attraction_scale_factor_data, + access_location::device, + access_mode::read); + ArrayHandle h_shape(m_attraction_scale_factor_shape, + access_location::host, + access_mode::read); + ArrayHandle h_domain(m_domain, access_location::host, access_mode::read); + + // build the interpolator: lo = {y_lo, t_lo}, hi = {y_hi, t_hi} + const Scalar lo[2] = {h_domain.data[0], h_domain.data[2]}; + const Scalar hi[2] = {h_domain.data[1], h_domain.data[3]}; + const LinearInterpolator2D interp(d_data.data, h_shape.data, lo, hi); + + ArrayHandle d_force(m_force, access_location::device, access_mode::overwrite); + + m_tuner->begin(); + gpu::compute_perturbed_lennard_jones_evap(d_force.data, + d_pos.data, + this->m_pdata->getGlobalBox(), + this->m_pdata->getN(), + d_n_neigh.data, + d_nlist.data, + d_head_list.data, + interp, + interface_height, + scaled_t, + lj1, + lj2, + rcutsq, + rwcasq, + epsilon_x_4, + m_energy_shift); + if (this->m_exec_conf->isCUDAErrorCheckingEnabled()) + CHECK_CUDA_ERROR(); + m_tuner->end(); + } + +namespace detail + { +namespace py = pybind11; + +void export_PerturbedLennardJonesEvapGPU(py::module& m); + } // end namespace detail + + } // end namespace azplugins + } // end namespace hoomd diff --git a/src/export_PerturbedLennardJonesEvapGPU.cc b/src/export_PerturbedLennardJonesEvapGPU.cc new file mode 100644 index 00000000..e038dd26 --- /dev/null +++ b/src/export_PerturbedLennardJonesEvapGPU.cc @@ -0,0 +1,70 @@ +// Copyright (c) 2018-2020, Michael P. Howard +// Copyright (c) 2021-2025, Auburn University +// Part of azplugins, released under the BSD 3-Clause License. + +#include "PerturbedLennardJonesEvapGPU.h" + +#include + +namespace hoomd + { +namespace azplugins + { +namespace detail + { +namespace py = pybind11; + +void export_PerturbedLennardJonesEvapGPU(py::module& m) + { + py::class_>(m, "PerturbedLennardJonesEvapGPU") + .def(py::init( + [](std::shared_ptr sysdef, + std::shared_ptr nlist, + Scalar rcut, + Scalar epsilon, + Scalar sigma, + Scalar scale_factor, + bool energy_shift, + py::array_t + attraction_scale_factor_data, + py::array_t + attraction_scale_factor_shape, + py::array_t domain, + std::shared_ptr variant) + { + if (attraction_scale_factor_shape.size() != 2) + throw std::runtime_error("lambda_shape must have 2 elements"); + if (domain.size() != 4) + throw std::runtime_error( + "domain must have 4 elements [y_lo, y_hi, t_lo, t_hi]"); + + const unsigned int* shape_ptr = attraction_scale_factor_shape.data(); + const Scalar* data_ptr = attraction_scale_factor_data.data(); + const Scalar* dom_ptr = domain.data(); + + if (attraction_scale_factor_data.size() + != static_cast(shape_ptr[0] * shape_ptr[1])) + throw std::runtime_error( + "attraction_scale_factor_data size does not match lambda_shape"); + + PairParametersPerturbedLennardJonesEvap params(epsilon, sigma); + return std::make_shared(sysdef, + nlist, + rcut, + scale_factor, + params, + energy_shift, + data_ptr, + shape_ptr, + dom_ptr, + variant); + })) + .def_property_readonly("rcut", &PerturbedLennardJonesEvapGPU::getRCut) + .def_property_readonly("epsilon", &PerturbedLennardJonesEvapGPU::getEpsilon) + .def_property_readonly("sigma", &PerturbedLennardJonesEvapGPU::getSigma); + } + } // end namespace detail + } // end namespace azplugins + } // end namespace hoomd diff --git a/src/module.cc b/src/module.cc index 4ffc9223..b79ace94 100644 --- a/src/module.cc +++ b/src/module.cc @@ -105,6 +105,7 @@ void export_PotentialPairColloidGPU(pybind11::module&); void export_PotentialPairExpandedYukawaGPU(pybind11::module&); void export_PotentialPairHertzGPU(pybind11::module&); void export_PotentialPairPerturbedLennardJonesGPU(pybind11::module&); +void export_PerturbedLennardJonesEvapGPU(pybind11::module&); // dpd void export_PotentialPairDPDThermoGeneralWeightGPU(pybind11::module&); @@ -181,6 +182,7 @@ PYBIND11_MODULE(_azplugins, m) export_PotentialPairExpandedYukawaGPU(m); export_PotentialPairHertzGPU(m); export_PotentialPairPerturbedLennardJonesGPU(m); + export_PerturbedLennardJonesEvapGPU(m); // dpd pair export_PotentialPairDPDThermoGeneralWeightGPU(m); diff --git a/src/pytest/test_perturbedljevap.py b/src/pytest/test_perturbedljevap.py index d0daffe5..be870668 100644 --- a/src/pytest/test_perturbedljevap.py +++ b/src/pytest/test_perturbedljevap.py @@ -19,7 +19,7 @@ def two_particle_snapshot_factory(): snap.configuration.box = [20, 20, 20, 0, 0, 0] snap.particles.N = 2 snap.particles.types = ["A"] - snap.particles.position[:] = [[-d/2, 0, 0], [d/2, 1.2, 0]] + snap.particles.position[:] = [[-d / 2, 0, 0], [d / 2, 1.2, 0]] return snap @@ -164,7 +164,7 @@ def test_energy_and_force_calculation_const( sim.run(0) - expected_forces = [[0.35032675, 0.4203921, 0.], [-0.35032675, -0.4203921, 0.]] + expected_forces = [[0.35032675, 0.4203921, 0.0], [-0.35032675, -0.4203921, 0.0]] expected_energies = [-0.07691956918309711, -0.07691956918309711] forces = evap.forces @@ -231,13 +231,15 @@ def test_energy_and_force_calculation_vary( sim.operations.integrator = integrator sim.run(0) - expected_forces = [[0.39995636970934184, 0.47994764365121023, 0.0], - [-0.39995636970934184, -0.47994764365121023, 0.0]] - + expected_forces = [ + [0.39995636970934184, 0.47994764365121023, 0.0], + [-0.39995636970934184, -0.47994764365121023, 0.0], + ] expected_energies = [-0.08781650815070254, -0.08781650815070254] forces = evap.forces energies = evap.energies + if sim.device.communicator.rank == 0: numpy.testing.assert_allclose(forces, expected_forces) numpy.testing.assert_allclose(energies, expected_energies) @@ -246,10 +248,12 @@ def test_energy_and_force_calculation_vary( """ sim.run(100) - - expected_energies = [-0.014101921016901138, -0.014101921016901138] - expected_forces = [[0.06422657031828848, 0.07707188438194618, 0.0], - [-0.06422657031828848, -0.07707188438194618, 0.0]] + + expected_forces = [ + [0.06422657031828848, 0.07707188438194618, 0.0], + [-0.06422657031828848, -0.07707188438194618, 0.0], + ] + expected_energies = [-0.014101921016901138, -0.014101921016901138] forces = evap.forces energies = evap.energies @@ -257,4 +261,3 @@ def test_energy_and_force_calculation_vary( if sim.device.communicator.rank == 0: numpy.testing.assert_allclose(forces, expected_forces) numpy.testing.assert_allclose(energies, expected_energies) - diff --git a/src/pytest/test_variant.py b/src/pytest/test_variant.py index 0963811d..a28f93fc 100644 --- a/src/pytest/test_variant.py +++ b/src/pytest/test_variant.py @@ -2,6 +2,8 @@ # Copyright (c) 2021-2025, Auburn University # Part of azplugins, released under the BSD 3-Clause License. +"""Test Variant Interpolation.""" + import hoomd import hoomd.azplugins import numpy From a36afcb4d1b5b9836c4c3df427e7f25d60747897 Mon Sep 17 00:00:00 2001 From: mayukh33 Date: Wed, 17 Jun 2026 10:49:01 +0200 Subject: [PATCH 12/17] Fix linking issue and continue GPU implementation --- src/CMakeLists.txt | 3 + src/PerturbedLennardJonesEvap.cc | 230 +++++++++++++++++++++ src/PerturbedLennardJonesEvap.h | 217 +------------------ src/PerturbedLennardJonesEvapGPU.cc | 105 ++++++++++ src/PerturbedLennardJonesEvapGPU.cu | 180 ++++++++++++++++ src/PerturbedLennardJonesEvapGPU.cuh | 229 +++----------------- src/PerturbedLennardJonesEvapGPU.h | 95 +-------- src/export_PerturbedLennardJonesEvapGPU.cc | 2 +- 8 files changed, 554 insertions(+), 507 deletions(-) create mode 100644 src/PerturbedLennardJonesEvap.cc create mode 100644 src/PerturbedLennardJonesEvapGPU.cc create mode 100644 src/PerturbedLennardJonesEvapGPU.cu diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b201a7c7..b372e150 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -5,6 +5,7 @@ set(COMPONENT_NAME azplugins) set(_${COMPONENT_NAME}_sources ConstantFlow.cc export_ImagePotentialBondHarmonic.cc + PerturbedLennardJonesEvap.cc export_PerturbedLennardJonesEvap.cc export_VariantInterpolated.cc module.cc @@ -14,6 +15,7 @@ set(_${COMPONENT_NAME}_sources if (ENABLE_HIP) list(APPEND _${COMPONENT_NAME}_sources export_ImagePotentialBondHarmonicGPU.cc + PerturbedLennardJonesEvapGPU.cc export_PerturbedLennardJonesEvapGPU.cc VelocityComputeGPU.cc ) @@ -22,6 +24,7 @@ endif() # TODO: List all GPU C++ source code files in _${COMPONENT_NAME}_cu_sources. set(_${COMPONENT_NAME}_cu_sources ImagePotentialBondGPUKernel.cu + PerturbedLennardJonesEvapGPU.cu VelocityComputeGPU.cu VelocityFieldComputeGPU.cu ) diff --git a/src/PerturbedLennardJonesEvap.cc b/src/PerturbedLennardJonesEvap.cc new file mode 100644 index 00000000..ae2ef6f9 --- /dev/null +++ b/src/PerturbedLennardJonesEvap.cc @@ -0,0 +1,230 @@ +// Copyright (c) 2018-2020, Michael P. Howard +// Copyright (c) 2021-2025, Auburn University +// Part of azplugins, released under the BSD 3-Clause License. + +#include "PerturbedLennardJonesEvap.h" + +#include + +namespace hoomd + { +namespace azplugins + { + +PerturbedLennardJonesEvap::PerturbedLennardJonesEvap( + std::shared_ptr sysdef, + std::shared_ptr nlist, + const Scalar rcut, + const Scalar time_scale_factor, + const param_type& params, + bool energy_shift, + const Scalar* attraction_scale_factor_data, + const unsigned int* attraction_scale_factor_shape, + const Scalar* domain, + std::shared_ptr variant) + : ForceCompute(sysdef), m_nlist(nlist), epsilon_x_4(params.epsilon_x_4), m_rcut(rcut), + m_time_scale_factor(time_scale_factor), + lj1(params.epsilon_x_4 * params.sigma_6 * params.sigma_6), + lj2(params.epsilon_x_4 * params.sigma_6), rcutsq(rcut * rcut), rwcasq(params.rwcasq), + sigma_6(params.sigma_6), m_energy_shift(energy_shift), m_variant(variant) + { + // Allocate and fill the (y, t) domain: [y_lo, y_hi, t_lo, t_hi] + { + GPUArray domain_arr(4, m_exec_conf); + m_domain.swap(domain_arr); + + ArrayHandle h_domain(m_domain, access_location::host, access_mode::overwrite); + std::copy(domain, domain + 4, h_domain.data); + } + + // Allocate and fill the table shape (ny, nt) + { + GPUArray shape_arr(2, m_exec_conf); + m_attraction_scale_factor_shape.swap(shape_arr); + + ArrayHandle h_shape(m_attraction_scale_factor_shape, + access_location::host, + access_mode::overwrite); + std::copy(attraction_scale_factor_shape, attraction_scale_factor_shape + 2, h_shape.data); + } + + { + const unsigned int n_data + = attraction_scale_factor_shape[0] * attraction_scale_factor_shape[1]; + + GPUArray attraction_scale_factor_arr(n_data, m_exec_conf); + m_attraction_scale_factor_data.swap(attraction_scale_factor_arr); + + ArrayHandle h_data(m_attraction_scale_factor_data, + access_location::host, + access_mode::overwrite); + + std::copy(attraction_scale_factor_data, attraction_scale_factor_data + n_data, h_data.data); + } + + { + const Index2D typpair_idx(m_pdata->getNTypes()); + m_r_cut_nlist + = std::make_shared>(typpair_idx.getNumElements(), m_exec_conf); + { + ArrayHandle h_r_cut_nlist(*m_r_cut_nlist, + access_location::host, + access_mode::overwrite); + + for (unsigned int i = 0; i < m_r_cut_nlist->getNumElements(); ++i) + h_r_cut_nlist.data[i] = m_rcut; + } + m_nlist->addRCutMatrix(m_r_cut_nlist); + m_nlist->notifyRCutMatrixChange(); + } + } + +void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) + { + m_nlist->compute(timestep); + + const bool third_law + = (m_nlist->getStorageMode() == hoomd::md::NeighborList::storageMode::half); + + Scalar scaled_t = scaleTime(timestep); + Scalar interface_height = (*m_variant)(timestep); + + ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::read); + m_force.zeroFill(); + + ArrayHandle h_force(m_force, access_location::host, access_mode::readwrite); + ArrayHandle h_data(m_attraction_scale_factor_data, + access_location::host, + access_mode::read); + ArrayHandle h_shape(m_attraction_scale_factor_shape, + access_location::host, + access_mode::read); + ArrayHandle h_domain(m_domain, access_location::host, access_mode::read); + + // Neighbor-list arrays + ArrayHandle h_n_neigh(m_nlist->getNNeighArray(), + access_location::host, + access_mode::read); + ArrayHandle h_nlist(m_nlist->getNListArray(), + access_location::host, + access_mode::read); + ArrayHandle h_head_list(m_nlist->getHeadList(), + access_location::host, + access_mode::read); + + // Build the interpolator: lo = {y_lo, t_lo}, hi = {y_hi, t_hi} + const Scalar lo[2] = {h_domain.data[0], h_domain.data[2]}; + const Scalar hi[2] = {h_domain.data[1], h_domain.data[3]}; + LinearInterpolator2D interp(h_data.data, h_shape.data, lo, hi); + + auto clamp_scaled_y = [](Scalar y, Scalar height) + { + const Scalar s = y / height; + if (!(s > Scalar(0.0))) + return Scalar(0.0); + if (s > Scalar(1.0)) + return Scalar(1.0); + return s; + }; + + const BoxDim box = m_pdata->getGlobalBox(); + const unsigned int N = m_pdata->getN(); + + for (unsigned int i = 0; i < N; ++i) + { + const Scalar3 pos_i = make_scalar3(h_pos.data[i].x, h_pos.data[i].y, h_pos.data[i].z); + + Scalar scaled_y_i = clamp_scaled_y(h_pos.data[i].y, interface_height); + Scalar3 fi = make_scalar3(0, 0, 0); + Scalar pei = 0; + + const Scalar attraction_scale_factor_i = interp(scaled_y_i, scaled_t); + + const unsigned int size = (unsigned int)h_n_neigh.data[i]; + const size_t head = h_head_list.data[i]; + + for (unsigned int k = 0; k < size; ++k) + { + const unsigned int j = h_nlist.data[head + k]; + if (j == i) + continue; + Scalar scaled_y_j = clamp_scaled_y(h_pos.data[j].y, interface_height); + Scalar3 pos_j = make_scalar3(h_pos.data[j].x, h_pos.data[j].y, h_pos.data[j].z); + + const Scalar attraction_scale_factor_j = interp(scaled_y_j, scaled_t); + + // Minimum-image + Scalar3 dx = pos_i - pos_j; + dx = box.minImage(dx); + + const Scalar rsq = dot(dx, dx); + const Scalar attraction_scale_factor_avg + = Scalar(0.5) * (attraction_scale_factor_i + attraction_scale_factor_j); + + const Scalar wca_shift + = epsilon_x_4 * (Scalar(1.0) - attraction_scale_factor_avg) / Scalar(4.0); + + if (rsq < rcutsq && lj1 != 0) + { + const Scalar r2inv = Scalar(1) / rsq; + const Scalar r6inv = r2inv * r2inv * r2inv; + + Scalar pair_eng = r6inv * (lj1 * r6inv - lj2); + Scalar force_divr + = r6inv * r2inv * (Scalar(12.0) * lj1 * r6inv - Scalar(6.0) * lj2); + + if (rsq < rwcasq) + { + pair_eng += wca_shift; + } + else + { + pair_eng *= attraction_scale_factor_avg; + force_divr *= attraction_scale_factor_avg; + } + + if (m_energy_shift) + { + const Scalar rcut2inv = Scalar(1.0) / rcutsq; + const Scalar rcut6inv = rcut2inv * rcut2inv * rcut2inv; + + Scalar pair_eng_shift = rcut6inv * (lj1 * rcut6inv - lj2); + + if (rcutsq < rwcasq) + { + pair_eng_shift += wca_shift; + } + else + { + pair_eng_shift *= attraction_scale_factor_avg; + } + + // Apply the shift to the pair energy + pair_eng -= pair_eng_shift; + } + + fi.x += force_divr * dx.x; + fi.y += force_divr * dx.y; + fi.z += force_divr * dx.z; + + pei += pair_eng; + + if (third_law) + { + h_force.data[j].x -= force_divr * dx.x; + h_force.data[j].y -= force_divr * dx.y; + h_force.data[j].z -= force_divr * dx.z; + h_force.data[j].w += Scalar(0.5) * pair_eng; + } + } + } + + h_force.data[i].x += fi.x; + h_force.data[i].y += fi.y; + h_force.data[i].z += fi.z; + h_force.data[i].w += Scalar(0.5) * pei; + } + } + + } // end namespace azplugins + } // end namespace hoomd diff --git a/src/PerturbedLennardJonesEvap.h b/src/PerturbedLennardJonesEvap.h index 76c2e631..1304b673 100644 --- a/src/PerturbedLennardJonesEvap.h +++ b/src/PerturbedLennardJonesEvap.h @@ -50,7 +50,7 @@ struct PairParametersPerturbedLennardJonesEvap auto sigma = v["sigma"].cast(); auto epsilon = v["epsilon"].cast(); - const Scalar sigma_2 = sigma * sigma; + const Scalar sigma_2 = sigma * sigma; const Scalar sigma_4 = sigma_2 * sigma_2; sigma_6 = sigma_2 * sigma_4; epsilon_x_4 = Scalar(4.0) * epsilon; @@ -132,221 +132,6 @@ class PerturbedLennardJonesEvap : public ForceCompute void computeForces(uint64_t timestep) override; }; -PerturbedLennardJonesEvap::PerturbedLennardJonesEvap( - std::shared_ptr sysdef, - std::shared_ptr nlist, - const Scalar rcut, - const Scalar time_scale_factor, - const param_type& params, - bool energy_shift, - const Scalar* attraction_scale_factor_data, - const unsigned int* attraction_scale_factor_shape, - const Scalar* domain, - std::shared_ptr variant) - : ForceCompute(sysdef), m_nlist(nlist), epsilon_x_4(params.epsilon_x_4), m_rcut(rcut), - m_time_scale_factor(time_scale_factor), - lj1(params.epsilon_x_4 * params.sigma_6 * params.sigma_6), - lj2(params.epsilon_x_4 * params.sigma_6), rcutsq(rcut * rcut), rwcasq(params.rwcasq), - sigma_6(params.sigma_6), m_energy_shift(energy_shift), m_variant(variant) - { - // Allocate and fill the (y, t) domain: [y_lo, y_hi, t_lo, t_hi] - { - GPUArray domain_arr(4, m_exec_conf); - m_domain.swap(domain_arr); - - ArrayHandle h_domain(m_domain, access_location::host, access_mode::overwrite); - std::copy(domain, domain + 4, h_domain.data); - } - - // Allocate and fill the table shape (ny, nt) - { - GPUArray shape_arr(2, m_exec_conf); - m_attraction_scale_factor_shape.swap(shape_arr); - - ArrayHandle h_shape(m_attraction_scale_factor_shape, - access_location::host, - access_mode::overwrite); - std::copy(attraction_scale_factor_shape, attraction_scale_factor_shape + 2, h_shape.data); - } - - { - const unsigned int n_data - = attraction_scale_factor_shape[0] * attraction_scale_factor_shape[1]; - - GPUArray attraction_scale_factor_arr(n_data, m_exec_conf); - m_attraction_scale_factor_data.swap(attraction_scale_factor_arr); - - ArrayHandle h_data(m_attraction_scale_factor_data, - access_location::host, - access_mode::overwrite); - - std::copy(attraction_scale_factor_data, attraction_scale_factor_data + n_data, h_data.data); - } - - { - const Index2D typpair_idx(m_pdata->getNTypes()); - m_r_cut_nlist - = std::make_shared>(typpair_idx.getNumElements(), m_exec_conf); - { - ArrayHandle h_r_cut_nlist(*m_r_cut_nlist, - access_location::host, - access_mode::overwrite); - - for (unsigned int i = 0; i < m_r_cut_nlist->getNumElements(); ++i) - h_r_cut_nlist.data[i] = m_rcut; - } - m_nlist->addRCutMatrix(m_r_cut_nlist); - m_nlist->notifyRCutMatrixChange(); - } - } - -void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) - { - m_nlist->compute(timestep); - - const bool third_law - = (m_nlist->getStorageMode() == hoomd::md::NeighborList::storageMode::half); - - Scalar scaled_t = scaleTime(timestep); - Scalar interface_height = (*m_variant)(timestep); - - ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::read); - m_force.zeroFill(); - - ArrayHandle h_force(m_force, access_location::host, access_mode::readwrite); - ArrayHandle h_data(m_attraction_scale_factor_data, - access_location::host, - access_mode::read); - ArrayHandle h_shape(m_attraction_scale_factor_shape, - access_location::host, - access_mode::read); - ArrayHandle h_domain(m_domain, access_location::host, access_mode::read); - - // Neighbor-list arrays - ArrayHandle h_n_neigh(m_nlist->getNNeighArray(), - access_location::host, - access_mode::read); - ArrayHandle h_nlist(m_nlist->getNListArray(), - access_location::host, - access_mode::read); - ArrayHandle h_head_list(m_nlist->getHeadList(), - access_location::host, - access_mode::read); - - // Build the interpolator: lo = {y_lo, t_lo}, hi = {y_hi, t_hi} - const Scalar lo[2] = {h_domain.data[0], h_domain.data[2]}; - const Scalar hi[2] = {h_domain.data[1], h_domain.data[3]}; - LinearInterpolator2D interp(h_data.data, h_shape.data, lo, hi); - - auto clamp_scaled_y = [](Scalar y, Scalar height) - { - const Scalar s = y / height; - if (!(s > Scalar(0.0))) - return Scalar(0.0); - if (s > Scalar(1.0)) - return Scalar(1.0); - return s; - }; - - const BoxDim box = m_pdata->getGlobalBox(); - const unsigned int N = m_pdata->getN(); - - for (unsigned int i = 0; i < N; ++i) - { - const Scalar3 pos_i = make_scalar3(h_pos.data[i].x, h_pos.data[i].y, h_pos.data[i].z); - - Scalar scaled_y_i = clamp_scaled_y(h_pos.data[i].y, interface_height); - Scalar3 fi = make_scalar3(0, 0, 0); - Scalar pei = 0; - - const Scalar attraction_scale_factor_i = interp(scaled_y_i, scaled_t); - - const unsigned int size = (unsigned int)h_n_neigh.data[i]; - const size_t head = h_head_list.data[i]; - - for (unsigned int k = 0; k < size; ++k) - { - const unsigned int j = h_nlist.data[head + k]; - if (j == i) - continue; - Scalar scaled_y_j = clamp_scaled_y(h_pos.data[j].y, interface_height); - Scalar3 pos_j = make_scalar3(h_pos.data[j].x, h_pos.data[j].y, h_pos.data[j].z); - - const Scalar attraction_scale_factor_j = interp(scaled_y_j, scaled_t); - - // Minimum-image - Scalar3 dx = pos_i - pos_j; - dx = box.minImage(dx); - - const Scalar rsq = dot(dx, dx); - const Scalar attraction_scale_factor_avg - = Scalar(0.5) * (attraction_scale_factor_i + attraction_scale_factor_j); - - const Scalar wca_shift - = epsilon_x_4 * (Scalar(1.0) - attraction_scale_factor_avg) / Scalar(4.0); - - if (rsq < rcutsq && lj1 != 0) - { - const Scalar r2inv = Scalar(1) / rsq; - const Scalar r6inv = r2inv * r2inv * r2inv; - - Scalar pair_eng = r6inv * (lj1 * r6inv - lj2); - Scalar force_divr - = r6inv * r2inv * (Scalar(12.0) * lj1 * r6inv - Scalar(6.0) * lj2); - - if (rsq < rwcasq) - { - pair_eng += wca_shift; - } - else - { - pair_eng *= attraction_scale_factor_avg; - force_divr *= attraction_scale_factor_avg; - } - - if (m_energy_shift) - { - const Scalar rcut2inv = Scalar(1.0) / rcutsq; - const Scalar rcut6inv = rcut2inv * rcut2inv * rcut2inv; - - Scalar pair_eng_shift = rcut6inv * (lj1 * rcut6inv - lj2); - - if (rcutsq < rwcasq) - { - pair_eng_shift += wca_shift; - } - else - { - pair_eng_shift *= attraction_scale_factor_avg; - } - - // Apply the shift to the pair energy - pair_eng -= pair_eng_shift; - } - - fi.x += force_divr * dx.x; - fi.y += force_divr * dx.y; - fi.z += force_divr * dx.z; - - pei += pair_eng; - - if (third_law) - { - h_force.data[j].x -= force_divr * dx.x; - h_force.data[j].y -= force_divr * dx.y; - h_force.data[j].z -= force_divr * dx.z; - h_force.data[j].w += Scalar(0.5) * pair_eng; - } - } - } - - h_force.data[i].x += fi.x; - h_force.data[i].y += fi.y; - h_force.data[i].z += fi.z; - h_force.data[i].w += Scalar(0.5) * pei; - } - } - namespace detail { void export_PerturbedLennardJonesEvap(pybind11::module& m); diff --git a/src/PerturbedLennardJonesEvapGPU.cc b/src/PerturbedLennardJonesEvapGPU.cc new file mode 100644 index 00000000..1d5956d2 --- /dev/null +++ b/src/PerturbedLennardJonesEvapGPU.cc @@ -0,0 +1,105 @@ +// Copyright (c) 2018-2020, Michael P. Howard +// Copyright (c) 2021-2025, Auburn University +// Part of azplugins, released under the BSD 3-Clause License. + +#include "PerturbedLennardJonesEvapGPU.h" + +#include + +namespace hoomd + { +namespace azplugins + { + +PerturbedLennardJonesEvapGPU::PerturbedLennardJonesEvapGPU( + std::shared_ptr sysdef, + std::shared_ptr nlist, + const Scalar r_cut, + const Scalar scale_factor, + const param_type& params, + bool energy_shift, + const Scalar* attraction_scale_factor_data, + const unsigned int* attraction_scale_factor_shape, + const Scalar* domain, + std::shared_ptr variant) + : PerturbedLennardJonesEvap(sysdef, + nlist, + r_cut, + scale_factor, + params, + energy_shift, + attraction_scale_factor_data, + attraction_scale_factor_shape, + domain, + variant) + { + m_tuner.reset(new Autotuner<1>({AutotunerBase::makeBlockSizeRange(m_exec_conf)}, + m_exec_conf, + "perturbed_lennard_jones_evap")); + m_autotuners.push_back(m_tuner); + } + +void PerturbedLennardJonesEvapGPU::computeForces(uint64_t timestep) + { + m_nlist->compute(timestep); + + const Scalar interface_height = (*m_variant)(timestep); + const Scalar scaled_t = scaleTime(timestep); + + ArrayHandle d_pos(m_pdata->getPositions(), access_location::device, access_mode::read); + + // neighbor-list arrays + ArrayHandle d_n_neigh(m_nlist->getNNeighArray(), + access_location::device, + access_mode::read); + ArrayHandle d_nlist(m_nlist->getNListArray(), + access_location::device, + access_mode::read); + ArrayHandle d_head_list(m_nlist->getHeadList(), + access_location::device, + access_mode::read); + + ArrayHandle d_data(m_attraction_scale_factor_data, + access_location::device, + access_mode::read); + ArrayHandle h_shape(m_attraction_scale_factor_shape, + access_location::host, + access_mode::read); + ArrayHandle h_domain(m_domain, access_location::host, access_mode::read); + + // build the interpolator: lo = {y_lo, t_lo}, hi = {y_hi, t_hi} + const Scalar lo[2] = {h_domain.data[0], h_domain.data[2]}; + const Scalar hi[2] = {h_domain.data[1], h_domain.data[3]}; + const LinearInterpolator2D interp(d_data.data, h_shape.data, lo, hi); + + ArrayHandle d_force(m_force, access_location::device, access_mode::overwrite); + + m_tuner->begin(); + gpu::perturbed_lennard_jones_evap_args_t args( + d_force.data, + d_pos.data, + this->m_pdata->getGlobalBox(), + this->m_pdata->getN(), + d_n_neigh.data, + d_nlist.data, + d_head_list.data, + interp, + scaled_t, + interface_height, + lj1, + lj2, + epsilon_x_4, + rcutsq, + rwcasq, + m_energy_shift, + m_tuner->getParam()[0]); + + gpu::compute_perturbed_lennard_jones_evap_forces(args); + if (this->m_exec_conf->isCUDAErrorCheckingEnabled()) + CHECK_CUDA_ERROR(); + m_tuner->end(); + + } + + } // end namespace azplugins + } // end namespace hoomd diff --git a/src/PerturbedLennardJonesEvapGPU.cu b/src/PerturbedLennardJonesEvapGPU.cu new file mode 100644 index 00000000..1dcd96ce --- /dev/null +++ b/src/PerturbedLennardJonesEvapGPU.cu @@ -0,0 +1,180 @@ +// Copyright (c) 2018-2020, Michael P. Howard +// Copyright (c) 2021-2025, Auburn University +// Part of azplugins, released under the BSD 3-Clause License. + +#include "PerturbedLennardJonesEvapGPU.cuh" + +namespace hoomd + { +namespace azplugins + { +namespace gpu + { +namespace kernel + { + +__device__ __inline__ Scalar clamp_scaled_y(Scalar y, Scalar height) + { + const Scalar s = y / height; + if (!(s > Scalar(0.0))) + return Scalar(0.0); + if (s > Scalar(1.0)) + return Scalar(1.0); + return s; + } + +__global__ void +compute_perturbed_lennard_jones_evap_forces(Scalar4* d_force, + const Scalar4* d_pos, + const BoxDim& box, + const unsigned int N, + const unsigned int* d_n_neigh, + const unsigned int* d_nlist, + const size_t* d_head_list, + const LinearInterpolator2D& interp, + const Scalar interface_height, + const Scalar scaled_t, + const Scalar lj1, + const Scalar lj2, + const Scalar epsilon_x_4, + const Scalar rcutsq, + const Scalar rwcasq, + const bool energy_shift) + { + // start by identifying which particle we are to handle + const unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx >= N) + return; + + const Scalar4 postype_i = __ldg(d_pos + idx); + const Scalar3 pos_i = make_scalar3(postype_i.x, postype_i.y, postype_i.z); + + const Scalar scaled_y_i = clamp_scaled_y(postype_i.y, interface_height); + const Scalar attraction_scale_factor_i = interp(scaled_y_i, scaled_t); + + // initialize the force and energy to 0 + Scalar3 fi = make_scalar3(0, 0, 0); + Scalar pei = 0; + + // load in the length of the neighbor list for this thread + const unsigned int n_neigh = d_n_neigh[idx]; + const size_t head = d_head_list[idx]; + + // loop over neighbors + for (unsigned int k = 0; k < n_neigh; ++k) + { + const unsigned int j = d_nlist[head + k]; + if (j == idx) + continue; + + const Scalar4 postype_j = __ldg(d_pos + j); + const Scalar3 pos_j = make_scalar3(postype_j.x, postype_j.y, postype_j.z); + + const Scalar scaled_y_j = clamp_scaled_y(postype_j.y, interface_height); + const Scalar attraction_scale_factor_j = interp(scaled_y_j, scaled_t); + + // minimum-image + Scalar3 dx = pos_i - pos_j; + dx = box.minImage(dx); + + const Scalar rsq = dot(dx, dx); + const Scalar attraction_scale_factor_avg + = Scalar(0.5) * (attraction_scale_factor_i + attraction_scale_factor_j); + + const Scalar wca_shift + = epsilon_x_4 * (Scalar(1.0) - attraction_scale_factor_avg) / Scalar(4.0); + + if (rsq < rcutsq && lj1 != 0) + { + const Scalar r2inv = Scalar(1) / rsq; + const Scalar r6inv = r2inv * r2inv * r2inv; + + Scalar pair_eng = r6inv * (lj1 * r6inv - lj2); + Scalar force_divr = r6inv * r2inv * (Scalar(12.0) * lj1 * r6inv - Scalar(6.0) * lj2); + + if (rsq < rwcasq) + { + pair_eng += wca_shift; + } + else + { + pair_eng *= attraction_scale_factor_avg; + force_divr *= attraction_scale_factor_avg; + } + + if (energy_shift) + { + const Scalar rcut2inv = Scalar(1.0) / rcutsq; + const Scalar rcut6inv = rcut2inv * rcut2inv * rcut2inv; + + Scalar pair_eng_shift = rcut6inv * (lj1 * rcut6inv - lj2); + + if (rcutsq < rwcasq) + { + pair_eng_shift += wca_shift; + } + else + { + pair_eng_shift *= attraction_scale_factor_avg; + } + + // apply the shift to the pair energy + pair_eng -= pair_eng_shift; + } + + fi.x += force_divr * dx.x; + fi.y += force_divr * dx.y; + fi.z += force_divr * dx.z; + + pei += pair_eng; + } + } + + d_force[idx] = make_scalar4(fi.x, fi.y, fi.z, Scalar(0.5) * pei); + } + + } // end namespace kernel + +hipError_t +compute_perturbed_lennard_jones_evap_forces(const perturbed_lennard_jones_evap_args_t& args) + { + unsigned int max_block_size; + hipFuncAttributes attr; + hipFuncGetAttributes( + &attr, + reinterpret_cast(&kernel::compute_perturbed_lennard_jones_evap_forces)); + max_block_size = attr.maxThreadsPerBlock; + + const unsigned int run_block_size = min(args.block_size, max_block_size); + + dim3 grid(args.N / run_block_size + 1, 1, 1); + dim3 threads(run_block_size, 1, 1); + + hipLaunchKernelGGL((azplugins::gpu::kernel::compute_perturbed_lennard_jones_evap_forces), + grid, + threads, + 0, + 0, + args.d_force, + args.d_pos, + args.box, + args.N, + args.d_n_neigh, + args.d_nlist, + args.d_head_list, + args.interp, + args.interface_height, + args.scaled_t, + args.lj1, + args.lj2, + args.epsilon_x_4, + args.rcutsq, + args.rwcasq, + args.energy_shift); + + return hipSuccess; + } + + } // end namespace gpu + } // end namespace azplugins + } // end namespace hoomd diff --git a/src/PerturbedLennardJonesEvapGPU.cuh b/src/PerturbedLennardJonesEvapGPU.cuh index ffcd182e..02db8f81 100644 --- a/src/PerturbedLennardJonesEvapGPU.cuh +++ b/src/PerturbedLennardJonesEvapGPU.cuh @@ -24,21 +24,17 @@ namespace azplugins { namespace gpu { - //! Wraps arguments to kernel driver -struct compute_perturbed_lennard_jones_evap_args_t +struct perturbed_lennard_jones_evap_args_t { - //! Construct a compute_perturbed_lennard_jones_evap_args_t - compute_perturbed_lennard_jones_evap_args_t(Scalar4* _d_force, + perturbed_lennard_jones_evap_args_t(Scalar4* _d_force, const Scalar4* _d_pos, const BoxDim& _box, const unsigned int _N, const unsigned int* _d_n_neigh, const unsigned int* _d_nlist, const size_t* _d_head_list, - const Scalar* _d_attraction_scale_factor_data, - const uint2 _attraction_scale_factor_shape, - const Scalar4 _domain, + const LinearInterpolator2D& _interp, const Scalar _scaled_t, const Scalar _interface_height, const Scalar _lj1, @@ -48,204 +44,33 @@ struct compute_perturbed_lennard_jones_evap_args_t const Scalar _rwcasq, const bool _energy_shift, const unsigned int _block_size) - : d_force(_d_force), d_pos(_d_pos), box(_box), d_n_neigh(_d_n_neigh), d_nlist(_d_nlist), - d_head_list(_d_head_list), - d_attraction_scale_factor_data(_d_attraction_scale_factor_data), - attraction_scale_factor_shape(_attraction_scale_factor_shape), domain(_domain), - scaled_t(_scaled_t), interface_height(_interface_height), lj1(_lj1), lj2(_lj2), - epsilon_x_4(_epsilon_x_4), rcutsq(_rcutsq), rwcasq(_rwcasq), energy_shift(_energy_shift), - block_size(_block_size) { }; - - Scalar4* d_force; //!< Force to write out - const Scalar4* d_pos; //!< Particle positions - const BoxDim box; //!< Box dimensions - const unsigned int N; //!< Number of particles - const unsigned int* d_n_neigh; //!< Number of neighbors per particle - const unsigned int* d_nlist; //!< Neighbor list - const size_t* d_head_list; //!< Head list indexing into the neighbor list - const Scalar* d_attraction_scale_factor_data; //!< Flattened (y, t) table data - const uint2 attraction_scale_factor_shape; //!< Table shape (ny, nt) - const Scalar4 domain; //!< Table domain [y_lo, y_hi, t_lo, t_hi] - const Scalar scaled_t; //!< Scaled simulation time - const Scalar interface_height; //!< Interface height at this timestep - const Scalar lj1; //!< 4 epsilon sigma^12 - const Scalar lj2; //!< 4 epsilon sigma^6 - const Scalar epsilon_x_4; //!< 4 epsilon - const Scalar rcutsq; //!< Squared cutoff radius - const Scalar rwcasq; //!< Squared WCA (potential minimum) radius - const bool energy_shift; //!< Whether to shift the energy at the cutoff - const unsigned int block_size; //!< Block size to execute + : d_force(_d_force), d_pos(_d_pos), box(_box), N(_N), d_n_neigh(_d_n_neigh), + d_nlist(_d_nlist), d_head_list(_d_head_list), interp(_interp), scaled_t(_scaled_t), + interface_height(_interface_height), lj1(_lj1), lj2(_lj2), epsilon_x_4(_epsilon_x_4), + rcutsq(_rcutsq), rwcasq(_rwcasq), energy_shift(_energy_shift), block_size(_block_size) { }; + + Scalar4* d_force; + const Scalar4* d_pos; + const BoxDim box; + const unsigned int N; + const unsigned int* d_n_neigh; + const unsigned int* d_nlist; + const size_t* d_head_list; + const LinearInterpolator2D interp; + const Scalar scaled_t; + const Scalar interface_height; + const Scalar lj1; + const Scalar lj2; + const Scalar epsilon_x_4; + const Scalar rcutsq; + const Scalar rwcasq; + const bool energy_shift; + const unsigned int block_size; }; -hipError_t compute_perturbed_lennard_jones_evap_forces( - const compute_perturbed_lennard_jones_evap_args_t& args); - -#ifdef __HIPCC__ - -namespace kernel - { - -__device__ __forceinline__ Scalar clamp_scaled_y(Scalar y, Scalar height) - { - const Scalar s = y / height; - if (!(s > Scalar(0.0))) - return Scalar(0.0); - if (s > Scalar(1.0)) - return Scalar(1.0); - return s; - } - -__global__ void -compute_perturbed_lennard_jones_evap_forces(Scalar4* d_force, - const Scalar4* d_pos, - const BoxDim& box, - const unsigned int N, - const unsigned int* d_n_neigh, - const unsigned int* d_nlist, - const size_t* d_head_list, - const LinearInterpolator2D& interp, - const Scalar interface_height, - const Scalar scaled_t, - const Scalar lj1, - const Scalar lj2, - const Scalar epsilon_x_4, - const Scalar rcutsq, - const Scalar rwcasq, - const bool energy_shift) - { - // start by identifying which particle we are to handle - const unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x; - if (idx >= N) - return; - - const Scalar4 postype_i = __ldg(d_pos + idx); - const Scalar3 pos_i = make_scalar3(postype_i.x, postype_i.y, postype_i.z); - - const Scalar scaled_y_i = clamp_scaled_y(postype_i.y, interface_height); - const Scalar attraction_scale_factor_i = interp(scaled_y_i, scaled_t); - - // initialize the force and energy to 0 - Scalar3 fi = make_scalar3(0, 0, 0); - Scalar pei(0); - - // load in the length of the neighbor list for this thread - const unsigned int n_neigh = d_n_neigh[idx]; - const size_t head = d_head_list[idx]; - - // loop over neighbors - for (unsigned int k = 0; k < n_neigh; ++k) - { - const unsigned int j = d_nlist[head + k]; - if (j == idx) - continue; - - const Scalar4 postype_j = __ldg(d_pos + j); - const Scalar3 pos_j = make_scalar3(postype_j.x, postype_j.y, postype_j.z); - - const Scalar scaled_y_j = clamp_scaled_y(postype_j.y, interface_height); - const Scalar attraction_scale_factor_j = interp(scaled_y_j, scaled_t); - - // minimum-image - Scalar3 dx = pos_i - pos_j; - dx = box.minImage(dx); - - const Scalar rsq = dot(dx, dx); - const Scalar attraction_scale_factor_avg - = Scalar(0.5) * (attraction_scale_factor_i + attraction_scale_factor_j); - const Scalar wca_shift - = epsilon_x_4 * (Scalar(1.0) - attraction_scale_factor_avg) / Scalar(4.0); - - if (rsq < rcutsq && lj1 != 0) - { - const Scalar r2inv = Scalar(1) / rsq; - const Scalar r6inv = r2inv * r2inv * r2inv; - - Scalar pair_eng = r6inv * (lj1 * r6inv - lj2); - Scalar force_divr = r6inv * r2inv * (Scalar(12.0) * lj1 * r6inv - Scalar(6.0) * lj2); - - if (rsq < rwcasq) - { - pair_eng += wca_shift; - } - else - { - pair_eng *= attraction_scale_factor_avg; - force_divr *= attraction_scale_factor_avg; - } - - if (energy_shift) - { - const Scalar rcut2inv = Scalar(1.0) / rcutsq; - const Scalar rcut6inv = rcut2inv * rcut2inv * rcut2inv; - - Scalar pair_eng_shift = rcut6inv * (lj1 * rcut6inv - lj2); - - if (rcutsq < rwcasq) - { - pair_eng_shift += wca_shift; - } - else - { - pair_eng_shift *= attraction_scale_factor_avg; - } - - // apply the shift to the pair energy - pair_eng -= pair_eng_shift; - } - - fi.x += force_divr * dx.x; - fi.y += force_divr * dx.y; - fi.z += force_divr * dx.z; - - pei += pair_eng; - } - } - - d_force[idx] = make_scalar4(fi.x, fi.y, fi.z, Scalar(0.5) * pei); - } - - } // end namespace kernel - -hipError_t -compute_perturbed_lennard_jones_evap_forces(const compute_perturbed_lennard_jones_evap_args_t& args) - { - unsigned int max_block_size; - hipFuncAttributes attr; - hipFuncGetAttributes( - &attr, - reinterpret_cast(&kernel::compute_perturbed_lennard_jones_evap_forces)); - max_block_size = attr.maxThreadsPerBlock; - - const unsigned int run_block_size = min(args.block_size, max_block_size); - - // setup the grid to run the kernel - dim3 grid(args.N / run_block_size + 1, 1, 1); - dim3 threads(run_block_size, 1, 1); - - // run the kernel - hipLaunchKernelGGL((kernel::compute_perturbed_lennard_jones_evap_forces), - grid, - threads, - args.d_force, - args.d_pos, - args.N, - args.d_n_neigh, - args.d_nlist, - args.d_head_list, - args.interp, - args.scaled_t, - args.interface_height, - args.lj1, - args.lj2, - args.epsilon_x_4, - args.rcutsq, - args.rwcasq, - args.energy_shift); - - return hipSuccess; - } -#endif // __HIPCC__ +hipError_t compute_perturbed_lennard_jones_evap_forces( + const perturbed_lennard_jones_evap_args_t& args); } // end namespace gpu } // end namespace azplugins diff --git a/src/PerturbedLennardJonesEvapGPU.h b/src/PerturbedLennardJonesEvapGPU.h index bf6c335c..8672df7f 100644 --- a/src/PerturbedLennardJonesEvapGPU.h +++ b/src/PerturbedLennardJonesEvapGPU.h @@ -2,8 +2,11 @@ // Copyright (c) 2021-2025, Auburn University // Part of azplugins, released under the BSD 3-Clause License. +#ifndef AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_GPU_H_ +#define AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_GPU_H_ + #include "PerturbedLennardJonesEvapGPU.cuh" -#include "PerturbedLennardJonesEvapGPU.h" +#include "PerturbedLennardJonesEvap.h" namespace hoomd { @@ -28,97 +31,13 @@ class PYBIND11_EXPORT PerturbedLennardJonesEvapGPU : public PerturbedLennardJone std::shared_ptr> m_tuner; void computeForces(uint64_t timestep) override; }; -PerturbedLennardJonesEvapGPU::PerturbedLennardJonesEvapGPU( - std::shared_ptr sysdef, - std::shared_ptr nlist, - const Scalar r_cut, - const Scalar scale_factor, - const param_type& params, - bool energy_shift, - const Scalar* attraction_scale_factor_data, - const unsigned int* attraction_scale_factor_shape, - const Scalar* domain, - std::shared_ptr variant) - : PerturbedLennardJonesEvap(sysdef, - nlist, - r_cut, - scale_factor, - params, - energy_shift, - attraction_scale_factor_data, - attraction_scale_factor_shape, - domain, - variant) - { - m_tuner.reset(new Autotuner<1>({AutotunerBase::makeBlockSizeRange(m_exec_conf)}, - m_exec_conf, - "perturbed_lennard_jones_evap")); - m_autotuners.push_back(m_tuner); - } - -void PerturbedLennardJonesEvapGPU::computeForces(uint64_t timestep) - { - m_nlist->compute(timestep); - - const Scalar interface_height = (*m_variant)(timestep); - const Scalar scaled_t = scaleTime(timestep); - - ArrayHandle d_pos(m_pdata->getPositions(), access_location::device, access_mode::read); - - // neighbor-list arrays - ArrayHandle d_n_neigh(m_nlist->getNNeighArray(), - access_location::device, - access_mode::read); - ArrayHandle d_nlist(m_nlist->getNListArray(), - access_location::device, - access_mode::read); - ArrayHandle d_head_list(m_nlist->getHeadList(), - access_location::device, - access_mode::read); - - ArrayHandle d_data(m_attraction_scale_factor_data, - access_location::device, - access_mode::read); - ArrayHandle h_shape(m_attraction_scale_factor_shape, - access_location::host, - access_mode::read); - ArrayHandle h_domain(m_domain, access_location::host, access_mode::read); - - // build the interpolator: lo = {y_lo, t_lo}, hi = {y_hi, t_hi} - const Scalar lo[2] = {h_domain.data[0], h_domain.data[2]}; - const Scalar hi[2] = {h_domain.data[1], h_domain.data[3]}; - const LinearInterpolator2D interp(d_data.data, h_shape.data, lo, hi); - - ArrayHandle d_force(m_force, access_location::device, access_mode::overwrite); - - m_tuner->begin(); - gpu::compute_perturbed_lennard_jones_evap(d_force.data, - d_pos.data, - this->m_pdata->getGlobalBox(), - this->m_pdata->getN(), - d_n_neigh.data, - d_nlist.data, - d_head_list.data, - interp, - interface_height, - scaled_t, - lj1, - lj2, - rcutsq, - rwcasq, - epsilon_x_4, - m_energy_shift); - if (this->m_exec_conf->isCUDAErrorCheckingEnabled()) - CHECK_CUDA_ERROR(); - m_tuner->end(); - } namespace detail { -namespace py = pybind11; - -void export_PerturbedLennardJonesEvapGPU(py::module& m); +void export_PerturbedLennardJonesEvapGPU(pybind11::module& m); } // end namespace detail } // end namespace azplugins } // end namespace hoomd + +#endif // AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_GPU_H_ diff --git a/src/export_PerturbedLennardJonesEvapGPU.cc b/src/export_PerturbedLennardJonesEvapGPU.cc index e038dd26..32ecf21d 100644 --- a/src/export_PerturbedLennardJonesEvapGPU.cc +++ b/src/export_PerturbedLennardJonesEvapGPU.cc @@ -17,7 +17,7 @@ namespace py = pybind11; void export_PerturbedLennardJonesEvapGPU(py::module& m) { py::class_>(m, "PerturbedLennardJonesEvapGPU") .def(py::init( [](std::shared_ptr sysdef, From 3537a131108343ee6f976bf2d0663672f1d5b566 Mon Sep 17 00:00:00 2001 From: mayukh33 Date: Sun, 21 Jun 2026 16:55:18 +0200 Subject: [PATCH 13/17] Fix scaled y coordinate for the 2D linear interpolant --- src/PerturbedLennardJonesEvap.cc | 18 ++--- src/PerturbedLennardJonesEvap.h | 2 +- src/PerturbedLennardJonesEvapGPU.cc | 44 ++++++------- src/PerturbedLennardJonesEvapGPU.cu | 8 +-- src/PerturbedLennardJonesEvapGPU.cuh | 76 +++++++++++----------- src/PerturbedLennardJonesEvapGPU.h | 7 +- src/VariantInterpolated.h | 4 +- src/export_PerturbedLennardJonesEvap.cc | 5 +- src/export_PerturbedLennardJonesEvapGPU.cc | 6 +- src/pytest/test_perturbedljevap.py | 33 +++++----- src/pytest/test_variant.py | 1 + 11 files changed, 103 insertions(+), 101 deletions(-) diff --git a/src/PerturbedLennardJonesEvap.cc b/src/PerturbedLennardJonesEvap.cc index ae2ef6f9..7481d23b 100644 --- a/src/PerturbedLennardJonesEvap.cc +++ b/src/PerturbedLennardJonesEvap.cc @@ -28,13 +28,13 @@ PerturbedLennardJonesEvap::PerturbedLennardJonesEvap( lj2(params.epsilon_x_4 * params.sigma_6), rcutsq(rcut * rcut), rwcasq(params.rwcasq), sigma_6(params.sigma_6), m_energy_shift(energy_shift), m_variant(variant) { - // Allocate and fill the (y, t) domain: [y_lo, y_hi, t_lo, t_hi] + // Allocate and fill the domain for time: [t_lo, t_hi] { - GPUArray domain_arr(4, m_exec_conf); + GPUArray domain_arr(2, m_exec_conf); m_domain.swap(domain_arr); ArrayHandle h_domain(m_domain, access_location::host, access_mode::overwrite); - std::copy(domain, domain + 4, h_domain.data); + std::copy(domain, domain + 2, h_domain.data); } // Allocate and fill the table shape (ny, nt) @@ -113,13 +113,13 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) access_mode::read); // Build the interpolator: lo = {y_lo, t_lo}, hi = {y_hi, t_hi} - const Scalar lo[2] = {h_domain.data[0], h_domain.data[2]}; - const Scalar hi[2] = {h_domain.data[1], h_domain.data[3]}; + const Scalar lo[2] = {Scalar(0.0), h_domain.data[0]}; + const Scalar hi[2] = {Scalar(1.0), h_domain.data[1]}; // y is always scaled between 0 and 1 LinearInterpolator2D interp(h_data.data, h_shape.data, lo, hi); - auto clamp_scaled_y = [](Scalar y, Scalar height) + auto clamp_scaled_y = [](Scalar y, Scalar height, Scalar ylo) { - const Scalar s = y / height; + const Scalar s = (y - ylo) / (height - ylo); if (!(s > Scalar(0.0))) return Scalar(0.0); if (s > Scalar(1.0)) @@ -134,7 +134,7 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) { const Scalar3 pos_i = make_scalar3(h_pos.data[i].x, h_pos.data[i].y, h_pos.data[i].z); - Scalar scaled_y_i = clamp_scaled_y(h_pos.data[i].y, interface_height); + Scalar scaled_y_i = clamp_scaled_y(h_pos.data[i].y, interface_height, box.getLo().y); Scalar3 fi = make_scalar3(0, 0, 0); Scalar pei = 0; @@ -148,7 +148,7 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) const unsigned int j = h_nlist.data[head + k]; if (j == i) continue; - Scalar scaled_y_j = clamp_scaled_y(h_pos.data[j].y, interface_height); + Scalar scaled_y_j = clamp_scaled_y(h_pos.data[j].y, interface_height, box.getLo().y); Scalar3 pos_j = make_scalar3(h_pos.data[j].x, h_pos.data[j].y, h_pos.data[j].z); const Scalar attraction_scale_factor_j = interp(scaled_y_j, scaled_t); diff --git a/src/PerturbedLennardJonesEvap.h b/src/PerturbedLennardJonesEvap.h index 1304b673..67d72ee1 100644 --- a/src/PerturbedLennardJonesEvap.h +++ b/src/PerturbedLennardJonesEvap.h @@ -50,7 +50,7 @@ struct PairParametersPerturbedLennardJonesEvap auto sigma = v["sigma"].cast(); auto epsilon = v["epsilon"].cast(); - const Scalar sigma_2 = sigma * sigma; + const Scalar sigma_2 = sigma * sigma; const Scalar sigma_4 = sigma_2 * sigma_2; sigma_6 = sigma_2 * sigma_4; epsilon_x_4 = Scalar(4.0) * epsilon; diff --git a/src/PerturbedLennardJonesEvapGPU.cc b/src/PerturbedLennardJonesEvapGPU.cc index 1d5956d2..adf4de51 100644 --- a/src/PerturbedLennardJonesEvapGPU.cc +++ b/src/PerturbedLennardJonesEvapGPU.cc @@ -36,7 +36,7 @@ PerturbedLennardJonesEvapGPU::PerturbedLennardJonesEvapGPU( m_tuner.reset(new Autotuner<1>({AutotunerBase::makeBlockSizeRange(m_exec_conf)}, m_exec_conf, "perturbed_lennard_jones_evap")); - m_autotuners.push_back(m_tuner); + this->m_autotuners.push_back(m_tuner); } void PerturbedLennardJonesEvapGPU::computeForces(uint64_t timestep) @@ -68,37 +68,35 @@ void PerturbedLennardJonesEvapGPU::computeForces(uint64_t timestep) ArrayHandle h_domain(m_domain, access_location::host, access_mode::read); // build the interpolator: lo = {y_lo, t_lo}, hi = {y_hi, t_hi} - const Scalar lo[2] = {h_domain.data[0], h_domain.data[2]}; - const Scalar hi[2] = {h_domain.data[1], h_domain.data[3]}; + const Scalar lo[2] = {Scalar(0.0), h_domain.data[0]}; + const Scalar hi[2] = {Scalar(1.0), h_domain.data[1]}; const LinearInterpolator2D interp(d_data.data, h_shape.data, lo, hi); ArrayHandle d_force(m_force, access_location::device, access_mode::overwrite); - + m_tuner->begin(); - gpu::perturbed_lennard_jones_evap_args_t args( - d_force.data, - d_pos.data, - this->m_pdata->getGlobalBox(), - this->m_pdata->getN(), - d_n_neigh.data, - d_nlist.data, - d_head_list.data, - interp, - scaled_t, - interface_height, - lj1, - lj2, - epsilon_x_4, - rcutsq, - rwcasq, - m_energy_shift, - m_tuner->getParam()[0]); + gpu::perturbed_lennard_jones_evap_args_t args(d_force.data, + d_pos.data, + this->m_pdata->getGlobalBox(), + this->m_pdata->getN(), + d_n_neigh.data, + d_nlist.data, + d_head_list.data, + interp, + scaled_t, + interface_height, + lj1, + lj2, + epsilon_x_4, + rcutsq, + rwcasq, + m_energy_shift, + m_tuner->getParam()[0]); gpu::compute_perturbed_lennard_jones_evap_forces(args); if (this->m_exec_conf->isCUDAErrorCheckingEnabled()) CHECK_CUDA_ERROR(); m_tuner->end(); - } } // end namespace azplugins diff --git a/src/PerturbedLennardJonesEvapGPU.cu b/src/PerturbedLennardJonesEvapGPU.cu index 1dcd96ce..6de6cf6b 100644 --- a/src/PerturbedLennardJonesEvapGPU.cu +++ b/src/PerturbedLennardJonesEvapGPU.cu @@ -13,9 +13,9 @@ namespace gpu namespace kernel { -__device__ __inline__ Scalar clamp_scaled_y(Scalar y, Scalar height) +__device__ __inline__ Scalar clamp_scaled_y(Scalar y, Scalar height, Scalar ylo) { - const Scalar s = y / height; + const Scalar s = (y - ylo) / (height - ylo); if (!(s > Scalar(0.0))) return Scalar(0.0); if (s > Scalar(1.0)) @@ -49,7 +49,7 @@ compute_perturbed_lennard_jones_evap_forces(Scalar4* d_force, const Scalar4 postype_i = __ldg(d_pos + idx); const Scalar3 pos_i = make_scalar3(postype_i.x, postype_i.y, postype_i.z); - const Scalar scaled_y_i = clamp_scaled_y(postype_i.y, interface_height); + const Scalar scaled_y_i = clamp_scaled_y(postype_i.y, interface_height, box.getLo().y); const Scalar attraction_scale_factor_i = interp(scaled_y_i, scaled_t); // initialize the force and energy to 0 @@ -70,7 +70,7 @@ compute_perturbed_lennard_jones_evap_forces(Scalar4* d_force, const Scalar4 postype_j = __ldg(d_pos + j); const Scalar3 pos_j = make_scalar3(postype_j.x, postype_j.y, postype_j.z); - const Scalar scaled_y_j = clamp_scaled_y(postype_j.y, interface_height); + const Scalar scaled_y_j = clamp_scaled_y(postype_j.y, interface_height, box.getLo().y); const Scalar attraction_scale_factor_j = interp(scaled_y_j, scaled_t); // minimum-image diff --git a/src/PerturbedLennardJonesEvapGPU.cuh b/src/PerturbedLennardJonesEvapGPU.cuh index 02db8f81..8b39dfc4 100644 --- a/src/PerturbedLennardJonesEvapGPU.cuh +++ b/src/PerturbedLennardJonesEvapGPU.cuh @@ -2,8 +2,6 @@ // Copyright (c) 2021-2025, Auburn University // Part of azplugins, released under the BSD 3-Clause License. - - #ifndef AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_GPU_CUH_ #define AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_GPU_CUH_ @@ -28,49 +26,49 @@ namespace gpu struct perturbed_lennard_jones_evap_args_t { perturbed_lennard_jones_evap_args_t(Scalar4* _d_force, - const Scalar4* _d_pos, - const BoxDim& _box, - const unsigned int _N, - const unsigned int* _d_n_neigh, - const unsigned int* _d_nlist, - const size_t* _d_head_list, - const LinearInterpolator2D& _interp, - const Scalar _scaled_t, - const Scalar _interface_height, - const Scalar _lj1, - const Scalar _lj2, - const Scalar _epsilon_x_4, - const Scalar _rcutsq, - const Scalar _rwcasq, - const bool _energy_shift, - const unsigned int _block_size) + const Scalar4* _d_pos, + const BoxDim& _box, + const unsigned int _N, + const unsigned int* _d_n_neigh, + const unsigned int* _d_nlist, + const size_t* _d_head_list, + const LinearInterpolator2D& _interp, + const Scalar _scaled_t, + const Scalar _interface_height, + const Scalar _lj1, + const Scalar _lj2, + const Scalar _epsilon_x_4, + const Scalar _rcutsq, + const Scalar _rwcasq, + const bool _energy_shift, + const unsigned int _block_size) : d_force(_d_force), d_pos(_d_pos), box(_box), N(_N), d_n_neigh(_d_n_neigh), d_nlist(_d_nlist), d_head_list(_d_head_list), interp(_interp), scaled_t(_scaled_t), interface_height(_interface_height), lj1(_lj1), lj2(_lj2), epsilon_x_4(_epsilon_x_4), - rcutsq(_rcutsq), rwcasq(_rwcasq), energy_shift(_energy_shift), block_size(_block_size) { }; + rcutsq(_rcutsq), rwcasq(_rwcasq), energy_shift(_energy_shift), + block_size(_block_size) { }; - Scalar4* d_force; - const Scalar4* d_pos; - const BoxDim box; - const unsigned int N; - const unsigned int* d_n_neigh; - const unsigned int* d_nlist; - const size_t* d_head_list; - const LinearInterpolator2D interp; - const Scalar scaled_t; - const Scalar interface_height; - const Scalar lj1; - const Scalar lj2; - const Scalar epsilon_x_4; - const Scalar rcutsq; - const Scalar rwcasq; - const bool energy_shift; - const unsigned int block_size; + Scalar4* d_force; + const Scalar4* d_pos; + const BoxDim box; + const unsigned int N; + const unsigned int* d_n_neigh; + const unsigned int* d_nlist; + const size_t* d_head_list; + const LinearInterpolator2D interp; + const Scalar scaled_t; + const Scalar interface_height; + const Scalar lj1; + const Scalar lj2; + const Scalar epsilon_x_4; + const Scalar rcutsq; + const Scalar rwcasq; + const bool energy_shift; + const unsigned int block_size; }; - -hipError_t compute_perturbed_lennard_jones_evap_forces( - const perturbed_lennard_jones_evap_args_t& args); +hipError_t +compute_perturbed_lennard_jones_evap_forces(const perturbed_lennard_jones_evap_args_t& args); } // end namespace gpu } // end namespace azplugins diff --git a/src/PerturbedLennardJonesEvapGPU.h b/src/PerturbedLennardJonesEvapGPU.h index 8672df7f..8b75d5e6 100644 --- a/src/PerturbedLennardJonesEvapGPU.h +++ b/src/PerturbedLennardJonesEvapGPU.h @@ -5,8 +5,13 @@ #ifndef AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_GPU_H_ #define AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_GPU_H_ -#include "PerturbedLennardJonesEvapGPU.cuh" +#ifdef __HIPCC__ +#error This header cannot be compiled by nvcc +#endif + #include "PerturbedLennardJonesEvap.h" +#include "PerturbedLennardJonesEvapGPU.cuh" +#include "hoomd/Autotuner.h" namespace hoomd { diff --git a/src/VariantInterpolated.h b/src/VariantInterpolated.h index 32962b38..468d9be4 100644 --- a/src/VariantInterpolated.h +++ b/src/VariantInterpolated.h @@ -45,10 +45,10 @@ class PYBIND11_EXPORT VariantInterpolated : public Variant Scalar operator()(uint64_t timestep) { ArrayHandle h_data(m_data, access_location::host, access_mode::read); - + const Scalar t = static_cast(timestep); // LinearInterpolator1D LinearInterpolator1D interp(h_data.data, m_n, m_t_lo, m_t_hi); - return interp(static_cast(timestep)); + return interp(t); } Scalar getTLo() const diff --git a/src/export_PerturbedLennardJonesEvap.cc b/src/export_PerturbedLennardJonesEvap.cc index 6160187b..95c2d870 100644 --- a/src/export_PerturbedLennardJonesEvap.cc +++ b/src/export_PerturbedLennardJonesEvap.cc @@ -36,9 +36,8 @@ void export_PerturbedLennardJonesEvap(py::module& m) { if (attraction_scale_factor_shape.size() != 2) throw std::runtime_error("lambda_shape must have 2 elements"); - if (domain.size() != 4) - throw std::runtime_error( - "domain must have 4 elements [y_lo, y_hi, t_lo, t_hi]"); + if (domain.size() != 2) + throw std::runtime_error("domain must have 2 elements [t_lo, t_hi]"); const unsigned int* shape_ptr = attraction_scale_factor_shape.data(); const Scalar* data_ptr = attraction_scale_factor_data.data(); diff --git a/src/export_PerturbedLennardJonesEvapGPU.cc b/src/export_PerturbedLennardJonesEvapGPU.cc index 32ecf21d..87699e29 100644 --- a/src/export_PerturbedLennardJonesEvapGPU.cc +++ b/src/export_PerturbedLennardJonesEvapGPU.cc @@ -17,7 +17,7 @@ namespace py = pybind11; void export_PerturbedLennardJonesEvapGPU(py::module& m) { py::class_>(m, "PerturbedLennardJonesEvapGPU") .def(py::init( [](std::shared_ptr sysdef, @@ -36,9 +36,9 @@ void export_PerturbedLennardJonesEvapGPU(py::module& m) { if (attraction_scale_factor_shape.size() != 2) throw std::runtime_error("lambda_shape must have 2 elements"); - if (domain.size() != 4) + if (domain.size() != 2) throw std::runtime_error( - "domain must have 4 elements [y_lo, y_hi, t_lo, t_hi]"); + "domain must have 2 elements [y_lo, y_hi, t_lo, t_hi]"); const unsigned int* shape_ptr = attraction_scale_factor_shape.data(); const Scalar* data_ptr = attraction_scale_factor_data.data(); diff --git a/src/pytest/test_perturbedljevap.py b/src/pytest/test_perturbedljevap.py index be870668..f447b43c 100644 --- a/src/pytest/test_perturbedljevap.py +++ b/src/pytest/test_perturbedljevap.py @@ -14,13 +14,14 @@ @pytest.fixture def two_particle_snapshot_factory(): snap = hoomd.Snapshot() - d = 1 if snap.communicator.rank == 0: snap.configuration.box = [20, 20, 20, 0, 0, 0] snap.particles.N = 2 snap.particles.types = ["A"] - snap.particles.position[:] = [[-d / 2, 0, 0], [d / 2, 1.2, 0]] - + snap.particles.position[:] = [ + [-10, -10, -10], + [-10, -8.8, -10], + ] return snap @@ -35,7 +36,7 @@ def valid_args_const(): "energy_shift": False, "attraction_scale_factor_data": numpy.array([[0.6, 0.6], [0.6, 0.6]]), "attraction_scale_factor_shape": [2, 2], - "domain": [0.0, 1.0, 0.0, 100.0], + "domain": [0.0, 100.0], "variant": hoomd.azplugins.variant.VariantInterpolated( [5.0, 4.0, 2.0, 1.0], 0, 300 ), @@ -73,7 +74,7 @@ def test_domain_mismatch( sim = simulation_factory(snap) bad_args = valid_args_const.copy() - bad_args["domain"] = [0.0, 10.0] + bad_args["domain"] = [10.0] integrator = hoomd.md.Integrator(dt=0.005) sim.operations.integrator = integrator @@ -164,8 +165,8 @@ def test_energy_and_force_calculation_const( sim.run(0) - expected_forces = [[0.35032675, 0.4203921, 0.0], [-0.35032675, -0.4203921, 0.0]] - expected_energies = [-0.07691956918309711, -0.07691956918309711] + expected_forces = [[0.0, 1.32701601, 0.0], [0.0, -1.32701601, 0.0]] + expected_energies = [-0.267289586, -0.267289586] forces = evap.forces energies = evap.energies @@ -178,8 +179,8 @@ def test_energy_and_force_calculation_const( @pytest.fixture( params=[ # (time_scale_factor, domain) - (1.0, [0.0, 1.0, 0.0, 100.0]), - (2.0, [0.0, 1.0, 0.0, 50.0]), + (1.0, [0.0, 100.0]), + (2.0, [0.0, 50.0]), ], ids=["unscaled_time", "scaled_time"], ) @@ -206,7 +207,7 @@ def valid_args_vary(request): "attraction_scale_factor_shape": attraction_factor_table.shape, "domain": domain, "variant": hoomd.azplugins.variant.VariantInterpolated( - [12, 10, 8, 6, 4, 2], + [2, 0, -2, -4, -6, -8], 0.0, 100.0, # Run for 100 timesteps, which corresponds to t = 0.5 for dt = 0.005 ), @@ -232,10 +233,10 @@ def test_energy_and_force_calculation_vary( sim.run(0) expected_forces = [ - [0.39995636970934184, 0.47994764365121023, 0.0], - [-0.39995636970934184, -0.47994764365121023, 0.0], + [0, 1.5150099394228083, 0.0], + [0, -1.5150099394228083, 0.0], ] - expected_energies = [-0.08781650815070254, -0.08781650815070254] + expected_energies = [-0.30515561, -0.30515561] forces = evap.forces energies = evap.energies @@ -250,10 +251,10 @@ def test_energy_and_force_calculation_vary( sim.run(100) expected_forces = [ - [0.06422657031828848, 0.07707188438194618, 0.0], - [-0.06422657031828848, -0.07707188438194618, 0.0], + [0.0, 0.24328626, 0.0], + [0.0, -0.24328626, 0.0], ] - expected_energies = [-0.014101921016901138, -0.014101921016901138] + expected_energies = [-0.04900309, -0.04900309] forces = evap.forces energies = evap.energies diff --git a/src/pytest/test_variant.py b/src/pytest/test_variant.py index a28f93fc..b9772809 100644 --- a/src/pytest/test_variant.py +++ b/src/pytest/test_variant.py @@ -20,6 +20,7 @@ def interpolated_eval(data, t_lo, t_hi): ([2.0, 4.0, 6.0, 8.0], 0, 300), ([1.0, -2.0, 5.0, 0.5, 3.0], 10, 90), ([-1.0, 1.0], 0, 100), + ([2, 0, -2, -4, -6, -8], 0, 100), ] From 75bcb2b26136025d893d6ce8407e7d732837a4e4 Mon Sep 17 00:00:00 2001 From: mayukh33 Date: Wed, 24 Jun 2026 19:57:25 +0200 Subject: [PATCH 14/17] Move interpolant out of the double for loops and other fixes --- src/PerturbedLennardJonesEvap.cc | 17 ++++-- src/PerturbedLennardJonesEvapGPU.cc | 13 ++++- src/PerturbedLennardJonesEvapGPU.cu | 77 ++++++++++++++++++---------- src/PerturbedLennardJonesEvapGPU.cuh | 14 +++-- src/VariantInterpolated.h | 17 +++++- 5 files changed, 99 insertions(+), 39 deletions(-) diff --git a/src/PerturbedLennardJonesEvap.cc b/src/PerturbedLennardJonesEvap.cc index 7481d23b..d031d493 100644 --- a/src/PerturbedLennardJonesEvap.cc +++ b/src/PerturbedLennardJonesEvap.cc @@ -129,16 +129,26 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) const BoxDim box = m_pdata->getGlobalBox(); const unsigned int N = m_pdata->getN(); + const unsigned int N_tot = N + m_pdata->getNGhosts(); + + GPUArray scale_factor(N_tot, m_exec_conf); + ArrayHandle h_scale_factor(scale_factor, access_location::host, access_mode::readwrite); + + for (unsigned int k = 0; k < N_tot; k++) + { + const Scalar scaled_pos_y + = clamp_scaled_y(h_pos.data[k].y, interface_height, box.getLo().y); + h_scale_factor.data[k] = interp(scaled_pos_y, scaled_t); + } for (unsigned int i = 0; i < N; ++i) { const Scalar3 pos_i = make_scalar3(h_pos.data[i].x, h_pos.data[i].y, h_pos.data[i].z); - Scalar scaled_y_i = clamp_scaled_y(h_pos.data[i].y, interface_height, box.getLo().y); Scalar3 fi = make_scalar3(0, 0, 0); Scalar pei = 0; - const Scalar attraction_scale_factor_i = interp(scaled_y_i, scaled_t); + const Scalar attraction_scale_factor_i = h_scale_factor.data[i]; const unsigned int size = (unsigned int)h_n_neigh.data[i]; const size_t head = h_head_list.data[i]; @@ -148,10 +158,9 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) const unsigned int j = h_nlist.data[head + k]; if (j == i) continue; - Scalar scaled_y_j = clamp_scaled_y(h_pos.data[j].y, interface_height, box.getLo().y); Scalar3 pos_j = make_scalar3(h_pos.data[j].x, h_pos.data[j].y, h_pos.data[j].z); - const Scalar attraction_scale_factor_j = interp(scaled_y_j, scaled_t); + const Scalar attraction_scale_factor_j = h_scale_factor.data[j]; // Minimum-image Scalar3 dx = pos_i - pos_j; diff --git a/src/PerturbedLennardJonesEvapGPU.cc b/src/PerturbedLennardJonesEvapGPU.cc index adf4de51..bd78a4a3 100644 --- a/src/PerturbedLennardJonesEvapGPU.cc +++ b/src/PerturbedLennardJonesEvapGPU.cc @@ -70,15 +70,26 @@ void PerturbedLennardJonesEvapGPU::computeForces(uint64_t timestep) // build the interpolator: lo = {y_lo, t_lo}, hi = {y_hi, t_hi} const Scalar lo[2] = {Scalar(0.0), h_domain.data[0]}; const Scalar hi[2] = {Scalar(1.0), h_domain.data[1]}; + const LinearInterpolator2D interp(d_data.data, h_shape.data, lo, hi); + const unsigned int N = this->m_pdata->getN(); + const unsigned int n_ghost = this->m_pdata->getNGhosts(); + + GPUArray scale_factor(N + n_ghost, m_exec_conf); + ArrayHandle d_scale_factor(scale_factor, + access_location::device, + access_mode::overwrite); + ArrayHandle d_force(m_force, access_location::device, access_mode::overwrite); m_tuner->begin(); gpu::perturbed_lennard_jones_evap_args_t args(d_force.data, d_pos.data, + d_scale_factor.data, this->m_pdata->getGlobalBox(), - this->m_pdata->getN(), + N, + n_ghost, d_n_neigh.data, d_nlist.data, d_head_list.data, diff --git a/src/PerturbedLennardJonesEvapGPU.cu b/src/PerturbedLennardJonesEvapGPU.cu index 6de6cf6b..1fcf8532 100644 --- a/src/PerturbedLennardJonesEvapGPU.cu +++ b/src/PerturbedLennardJonesEvapGPU.cu @@ -12,34 +12,47 @@ namespace gpu { namespace kernel { - __device__ __inline__ Scalar clamp_scaled_y(Scalar y, Scalar height, Scalar ylo) { const Scalar s = (y - ylo) / (height - ylo); - if (!(s > Scalar(0.0))) + if (s < Scalar(0.0)) return Scalar(0.0); if (s > Scalar(1.0)) return Scalar(1.0); return s; } -__global__ void -compute_perturbed_lennard_jones_evap_forces(Scalar4* d_force, - const Scalar4* d_pos, - const BoxDim& box, - const unsigned int N, - const unsigned int* d_n_neigh, - const unsigned int* d_nlist, - const size_t* d_head_list, - const LinearInterpolator2D& interp, - const Scalar interface_height, - const Scalar scaled_t, - const Scalar lj1, - const Scalar lj2, - const Scalar epsilon_x_4, - const Scalar rcutsq, - const Scalar rwcasq, - const bool energy_shift) +__global__ void compute_attraction_scale_factor(Scalar* d_scale_factor, + const Scalar4* d_pos, + const unsigned int Ntot, + const LinearInterpolator2D interp, + const Scalar interface_height, + const Scalar scaled_t, + const Scalar y_lo) + { + const unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i >= Ntot) + return; + + const Scalar y = __ldg(d_pos + i).y; + const Scalar scaled_pos_y = clamp_scaled_y(y, interface_height, y_lo); + d_scale_factor[i] = interp(scaled_pos_y, scaled_t); + } + +__global__ void compute_perturbed_lennard_jones_evap_forces(Scalar4* d_force, + const Scalar4* d_pos, + const BoxDim& box, + const unsigned int N, + const unsigned int* d_n_neigh, + const unsigned int* d_nlist, + const size_t* d_head_list, + const Scalar* d_scale_factor, + const Scalar lj1, + const Scalar lj2, + const Scalar epsilon_x_4, + const Scalar rcutsq, + const Scalar rwcasq, + const bool energy_shift) { // start by identifying which particle we are to handle const unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x; @@ -49,9 +62,7 @@ compute_perturbed_lennard_jones_evap_forces(Scalar4* d_force, const Scalar4 postype_i = __ldg(d_pos + idx); const Scalar3 pos_i = make_scalar3(postype_i.x, postype_i.y, postype_i.z); - const Scalar scaled_y_i = clamp_scaled_y(postype_i.y, interface_height, box.getLo().y); - const Scalar attraction_scale_factor_i = interp(scaled_y_i, scaled_t); - + const Scalar attraction_scale_factor_i = __ldg(d_scale_factor + idx); // initialize the force and energy to 0 Scalar3 fi = make_scalar3(0, 0, 0); Scalar pei = 0; @@ -70,8 +81,7 @@ compute_perturbed_lennard_jones_evap_forces(Scalar4* d_force, const Scalar4 postype_j = __ldg(d_pos + j); const Scalar3 pos_j = make_scalar3(postype_j.x, postype_j.y, postype_j.z); - const Scalar scaled_y_j = clamp_scaled_y(postype_j.y, interface_height, box.getLo().y); - const Scalar attraction_scale_factor_j = interp(scaled_y_j, scaled_t); + const Scalar attraction_scale_factor_j = __ldg(d_scale_factor + j); // minimum-image Scalar3 dx = pos_i - pos_j; @@ -147,9 +157,24 @@ compute_perturbed_lennard_jones_evap_forces(const perturbed_lennard_jones_evap_a const unsigned int run_block_size = min(args.block_size, max_block_size); + const unsigned int Ntot = args.N + args.N_ghost; dim3 grid(args.N / run_block_size + 1, 1, 1); dim3 threads(run_block_size, 1, 1); + const Scalar ylo = args.box.getLo().y; + hipLaunchKernelGGL((azplugins::gpu::kernel::compute_attraction_scale_factor), + grid, + threads, + 0, + 0, + args.d_scale_factor, + args.d_pos, + Ntot, + args.interp, + args.interface_height, + args.scaled_t, + ylo); + hipLaunchKernelGGL((azplugins::gpu::kernel::compute_perturbed_lennard_jones_evap_forces), grid, threads, @@ -162,9 +187,7 @@ compute_perturbed_lennard_jones_evap_forces(const perturbed_lennard_jones_evap_a args.d_n_neigh, args.d_nlist, args.d_head_list, - args.interp, - args.interface_height, - args.scaled_t, + args.d_scale_factor, args.lj1, args.lj2, args.epsilon_x_4, diff --git a/src/PerturbedLennardJonesEvapGPU.cuh b/src/PerturbedLennardJonesEvapGPU.cuh index 8b39dfc4..4a5ed95a 100644 --- a/src/PerturbedLennardJonesEvapGPU.cuh +++ b/src/PerturbedLennardJonesEvapGPU.cuh @@ -27,8 +27,10 @@ struct perturbed_lennard_jones_evap_args_t { perturbed_lennard_jones_evap_args_t(Scalar4* _d_force, const Scalar4* _d_pos, + Scalar* _d_scale_factor, const BoxDim& _box, const unsigned int _N, + const unsigned int _n_ghost, const unsigned int* _d_n_neigh, const unsigned int* _d_nlist, const size_t* _d_head_list, @@ -42,16 +44,18 @@ struct perturbed_lennard_jones_evap_args_t const Scalar _rwcasq, const bool _energy_shift, const unsigned int _block_size) - : d_force(_d_force), d_pos(_d_pos), box(_box), N(_N), d_n_neigh(_d_n_neigh), - d_nlist(_d_nlist), d_head_list(_d_head_list), interp(_interp), scaled_t(_scaled_t), - interface_height(_interface_height), lj1(_lj1), lj2(_lj2), epsilon_x_4(_epsilon_x_4), - rcutsq(_rcutsq), rwcasq(_rwcasq), energy_shift(_energy_shift), - block_size(_block_size) { }; + : d_force(_d_force), d_pos(_d_pos), d_scale_factor(_d_scale_factor), box(_box), N(_N), + N_ghost(_n_ghost), d_n_neigh(_d_n_neigh), d_nlist(_d_nlist), d_head_list(_d_head_list), + interp(_interp), scaled_t(_scaled_t), interface_height(_interface_height), lj1(_lj1), + lj2(_lj2), epsilon_x_4(_epsilon_x_4), rcutsq(_rcutsq), rwcasq(_rwcasq), + energy_shift(_energy_shift), block_size(_block_size) { }; Scalar4* d_force; const Scalar4* d_pos; + Scalar* d_scale_factor; const BoxDim box; const unsigned int N; + const unsigned int N_ghost; const unsigned int* d_n_neigh; const unsigned int* d_nlist; const size_t* d_head_list; diff --git a/src/VariantInterpolated.h b/src/VariantInterpolated.h index 468d9be4..901c13ab 100644 --- a/src/VariantInterpolated.h +++ b/src/VariantInterpolated.h @@ -47,8 +47,21 @@ class PYBIND11_EXPORT VariantInterpolated : public Variant ArrayHandle h_data(m_data, access_location::host, access_mode::read); const Scalar t = static_cast(timestep); // LinearInterpolator1D - LinearInterpolator1D interp(h_data.data, m_n, m_t_lo, m_t_hi); - return interp(t); + if (t < m_t_lo) + { + return h_data.data[0]; + } + + if (t > m_t_hi) + { + return h_data.data[m_n]; + } + + else + { + LinearInterpolator1D interp(h_data.data, m_n, m_t_lo, m_t_hi); + return interp(t); + } } Scalar getTLo() const From f0868318742dca13d76116676028b47f3b852c8b Mon Sep 17 00:00:00 2001 From: mayukh33 Date: Tue, 30 Jun 2026 22:56:21 +0200 Subject: [PATCH 15/17] Apply suggestions and fix GPU code implementation --- src/CMakeLists.txt | 2 - src/PerturbedLennardJonesEvap.cc | 150 ++++++++++++--------- src/PerturbedLennardJonesEvap.h | 72 +--------- src/PerturbedLennardJonesEvapGPU.cc | 127 ++++++++++++----- src/PerturbedLennardJonesEvapGPU.cu | 129 ++++++------------ src/PerturbedLennardJonesEvapGPU.cuh | 27 ++-- src/PerturbedLennardJonesEvapGPU.h | 18 ++- src/export_PerturbedLennardJonesEvap.cc | 70 ---------- src/export_PerturbedLennardJonesEvapGPU.cc | 70 ---------- src/pair.py | 13 +- src/pytest/test_perturbedljevap.py | 31 +---- 11 files changed, 270 insertions(+), 439 deletions(-) delete mode 100644 src/export_PerturbedLennardJonesEvap.cc delete mode 100644 src/export_PerturbedLennardJonesEvapGPU.cc diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b372e150..a3c1b1b3 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -6,7 +6,6 @@ set(_${COMPONENT_NAME}_sources ConstantFlow.cc export_ImagePotentialBondHarmonic.cc PerturbedLennardJonesEvap.cc - export_PerturbedLennardJonesEvap.cc export_VariantInterpolated.cc module.cc ParabolicFlow.cc @@ -16,7 +15,6 @@ if (ENABLE_HIP) list(APPEND _${COMPONENT_NAME}_sources export_ImagePotentialBondHarmonicGPU.cc PerturbedLennardJonesEvapGPU.cc - export_PerturbedLennardJonesEvapGPU.cc VelocityComputeGPU.cc ) endif() diff --git a/src/PerturbedLennardJonesEvap.cc b/src/PerturbedLennardJonesEvap.cc index d031d493..99150076 100644 --- a/src/PerturbedLennardJonesEvap.cc +++ b/src/PerturbedLennardJonesEvap.cc @@ -10,7 +10,8 @@ namespace hoomd { namespace azplugins { - +namespace detail + { PerturbedLennardJonesEvap::PerturbedLennardJonesEvap( std::shared_ptr sysdef, std::shared_ptr nlist, @@ -22,11 +23,8 @@ PerturbedLennardJonesEvap::PerturbedLennardJonesEvap( const unsigned int* attraction_scale_factor_shape, const Scalar* domain, std::shared_ptr variant) - : ForceCompute(sysdef), m_nlist(nlist), epsilon_x_4(params.epsilon_x_4), m_rcut(rcut), - m_time_scale_factor(time_scale_factor), - lj1(params.epsilon_x_4 * params.sigma_6 * params.sigma_6), - lj2(params.epsilon_x_4 * params.sigma_6), rcutsq(rcut * rcut), rwcasq(params.rwcasq), - sigma_6(params.sigma_6), m_energy_shift(energy_shift), m_variant(variant) + : ForceCompute(sysdef), m_nlist(nlist), m_rcut(rcut), m_time_scale_factor(time_scale_factor), + m_params(params), m_energy_shift(energy_shift), m_variant(variant) { // Allocate and fill the domain for time: [t_lo, t_hi] { @@ -117,30 +115,35 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) const Scalar hi[2] = {Scalar(1.0), h_domain.data[1]}; // y is always scaled between 0 and 1 LinearInterpolator2D interp(h_data.data, h_shape.data, lo, hi); - auto clamp_scaled_y = [](Scalar y, Scalar height, Scalar ylo) - { - const Scalar s = (y - ylo) / (height - ylo); - if (!(s > Scalar(0.0))) - return Scalar(0.0); - if (s > Scalar(1.0)) - return Scalar(1.0); - return s; - }; - const BoxDim box = m_pdata->getGlobalBox(); const unsigned int N = m_pdata->getN(); - const unsigned int N_tot = N + m_pdata->getNGhosts(); + const unsigned int Ntot = N + m_pdata->getNGhosts(); - GPUArray scale_factor(N_tot, m_exec_conf); - ArrayHandle h_scale_factor(scale_factor, access_location::host, access_mode::readwrite); + if (m_scale_factor.getNumElements() < Ntot) + { + GPUArray tmp(Ntot, m_exec_conf); + m_scale_factor.swap(tmp); + } + ArrayHandle h_scale_factor(m_scale_factor, + access_location::host, + access_mode::readwrite); - for (unsigned int k = 0; k < N_tot; k++) + for (unsigned int k = 0; k < Ntot; k++) { - const Scalar scaled_pos_y - = clamp_scaled_y(h_pos.data[k].y, interface_height, box.getLo().y); + Scalar scaled_pos_y + = (h_pos.data[k].y - box.getLo().y) / (interface_height - box.getLo().y); + if (scaled_pos_y < 0) + { + scaled_pos_y = Scalar(0.0); + } + if (scaled_pos_y > 1) + { + scaled_pos_y = Scalar(1.0); + } h_scale_factor.data[k] = interp(scaled_pos_y, scaled_t); } + Scalar rcutsq = m_rcut * m_rcut; for (unsigned int i = 0; i < N; ++i) { const Scalar3 pos_i = make_scalar3(h_pos.data[i].x, h_pos.data[i].y, h_pos.data[i].z); @@ -170,48 +173,15 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) const Scalar attraction_scale_factor_avg = Scalar(0.5) * (attraction_scale_factor_i + attraction_scale_factor_j); - const Scalar wca_shift - = epsilon_x_4 * (Scalar(1.0) - attraction_scale_factor_avg) / Scalar(4.0); - - if (rsq < rcutsq && lj1 != 0) - { - const Scalar r2inv = Scalar(1) / rsq; - const Scalar r6inv = r2inv * r2inv * r2inv; - - Scalar pair_eng = r6inv * (lj1 * r6inv - lj2); - Scalar force_divr - = r6inv * r2inv * (Scalar(12.0) * lj1 * r6inv - Scalar(6.0) * lj2); + m_params.attraction_scale_factor = attraction_scale_factor_avg; - if (rsq < rwcasq) - { - pair_eng += wca_shift; - } - else - { - pair_eng *= attraction_scale_factor_avg; - force_divr *= attraction_scale_factor_avg; - } - - if (m_energy_shift) - { - const Scalar rcut2inv = Scalar(1.0) / rcutsq; - const Scalar rcut6inv = rcut2inv * rcut2inv * rcut2inv; - - Scalar pair_eng_shift = rcut6inv * (lj1 * rcut6inv - lj2); - - if (rcutsq < rwcasq) - { - pair_eng_shift += wca_shift; - } - else - { - pair_eng_shift *= attraction_scale_factor_avg; - } - - // Apply the shift to the pair energy - pair_eng -= pair_eng_shift; - } + Scalar force_divr = Scalar(0.0); + Scalar pair_eng = Scalar(0.0); + PairEvaluatorPerturbedLennardJones eval(rsq, rcutsq, m_params); + bool evaluated = eval.evalForceAndEnergy(force_divr, pair_eng, m_energy_shift); + if (evaluated) + { fi.x += force_divr * dx.x; fi.y += force_divr * dx.y; fi.z += force_divr * dx.z; @@ -235,5 +205,61 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) } } +namespace py = pybind11; + +void export_PerturbedLennardJonesEvap(py::module& m) + { + py::class_>( + m, + "PerturbedLennardJonesEvap") + .def(py::init( + [](std::shared_ptr sysdef, + std::shared_ptr nlist, + Scalar rcut, + Scalar epsilon, + Scalar sigma, + Scalar time_scale_factor, + bool energy_shift, + py::array_t + attraction_scale_factor_data, + py::array_t + attraction_scale_factor_shape, + py::array_t domain, + std::shared_ptr variant) + { + if (attraction_scale_factor_shape.size() != 2) + throw std::runtime_error("lambda_shape must have 2 elements"); + if (domain.size() != 2) + throw std::runtime_error("domain must have 2 elements [t_lo, t_hi]"); + + const unsigned int* shape_ptr = attraction_scale_factor_shape.data(); + const Scalar* data_ptr = attraction_scale_factor_data.data(); + const Scalar* dom_ptr = domain.data(); + + if (attraction_scale_factor_data.size() + != static_cast(shape_ptr[0] * shape_ptr[1])) + throw std::runtime_error( + "attraction_scale_factor_data size does not match lambda_shape"); + + PairParametersPerturbedLennardJones params; + const Scalar sigma_2 = sigma * sigma; + params.sigma_6 = sigma_2 * sigma_2 * sigma_2; + params.epsilon_x_4 = Scalar(4.0) * epsilon; + params.rwcasq = std::pow(Scalar(2.0), Scalar(1.0) / Scalar(3.0)) * sigma_2; + + return std::make_shared(sysdef, + nlist, + rcut, + time_scale_factor, + params, + energy_shift, + data_ptr, + shape_ptr, + dom_ptr, + variant); + })); + } + + } // end namespace detail } // end namespace azplugins } // end namespace hoomd diff --git a/src/PerturbedLennardJonesEvap.h b/src/PerturbedLennardJonesEvap.h index 67d72ee1..ffca7365 100644 --- a/src/PerturbedLennardJonesEvap.h +++ b/src/PerturbedLennardJonesEvap.h @@ -18,59 +18,19 @@ #include "hoomd/md/NeighborList.h" #include "LinearInterpolator2D.h" +#include "PairEvaluatorPerturbedLennardJones.h" #include "VariantInterpolated.h" namespace hoomd { namespace azplugins { - -struct PairParametersPerturbedLennardJonesEvap +namespace detail { - Scalar epsilon_x_4; - Scalar sigma_6; - Scalar rwcasq; - -#ifndef __HIPCC__ - - //! Default constructor - PairParametersPerturbedLennardJonesEvap() : epsilon_x_4(0), sigma_6(0), rwcasq(0) { } - - PairParametersPerturbedLennardJonesEvap(Scalar epsilon, Scalar sigma) - { - const Scalar sigma_2 = sigma * sigma; - const Scalar sigma_4 = sigma_2 * sigma_2; - sigma_6 = sigma_2 * sigma_4; - epsilon_x_4 = Scalar(4.0) * epsilon; - rwcasq = std::pow(Scalar(2.0), Scalar(1.0) / Scalar(3.0)) * sigma_2; - } - - PairParametersPerturbedLennardJonesEvap(pybind11::dict v, bool managed = false) - { - auto sigma = v["sigma"].cast(); - auto epsilon = v["epsilon"].cast(); - - const Scalar sigma_2 = sigma * sigma; - const Scalar sigma_4 = sigma_2 * sigma_2; - sigma_6 = sigma_2 * sigma_4; - epsilon_x_4 = Scalar(4.0) * epsilon; - rwcasq = std::pow(Scalar(2.0), Scalar(1.0) / Scalar(3.0)) * sigma_2; - } - - pybind11::dict asDict() - { - pybind11::dict v; - v["sigma"] = std::pow(sigma_6, Scalar(1.0) / Scalar(6.0)); - v["epsilon"] = epsilon_x_4 / Scalar(4.0); - return v; - } -#endif // __HIPCC__ - }; - class PerturbedLennardJonesEvap : public ForceCompute { public: - typedef PairParametersPerturbedLennardJonesEvap param_type; + typedef detail::PairParametersPerturbedLennardJones param_type; PerturbedLennardJonesEvap(std::shared_ptr sysdef, std::shared_ptr nlist, @@ -95,35 +55,17 @@ class PerturbedLennardJonesEvap : public ForceCompute return Scalar(static_cast(timestep) / m_time_scale_factor); } - Scalar getRCut() const - { - return m_rcut; - } - - Scalar getEpsilon() const - { - return epsilon_x_4 / Scalar(4.0); - } - - Scalar getSigma() const - { - return std::pow(sigma_6, Scalar(1.0) / Scalar(6.0)); - } - protected: std::shared_ptr m_nlist; //!< Neighbor list Scalar epsilon_x_4; Scalar m_rcut; Scalar m_time_scale_factor; //!< Time scaling factor - Scalar lj1; - Scalar lj2; - Scalar rcutsq; - Scalar rwcasq; - Scalar sigma_6; + param_type m_params; bool m_energy_shift; - GPUArray m_domain; //!< [y_lo, y_hi, t_lo, t_hi] + GPUArray m_domain; //!< [t_lo, t_hi] GPUArray m_attraction_scale_factor_data; //!< Flattened (y, t) data GPUArray m_attraction_scale_factor_shape; //!< [ny, nt] + GPUArray m_scale_factor; std::shared_ptr m_variant; std::shared_ptr> @@ -132,8 +74,6 @@ class PerturbedLennardJonesEvap : public ForceCompute void computeForces(uint64_t timestep) override; }; -namespace detail - { void export_PerturbedLennardJonesEvap(pybind11::module& m); } // end namespace detail diff --git a/src/PerturbedLennardJonesEvapGPU.cc b/src/PerturbedLennardJonesEvapGPU.cc index bd78a4a3..c370f5ea 100644 --- a/src/PerturbedLennardJonesEvapGPU.cc +++ b/src/PerturbedLennardJonesEvapGPU.cc @@ -10,7 +10,8 @@ namespace hoomd { namespace azplugins { - +namespace detail + { PerturbedLennardJonesEvapGPU::PerturbedLennardJonesEvapGPU( std::shared_ptr sysdef, std::shared_ptr nlist, @@ -33,6 +34,14 @@ PerturbedLennardJonesEvapGPU::PerturbedLennardJonesEvapGPU( domain, variant) { + if (!this->m_exec_conf->isCUDAEnabled()) + { + this->m_exec_conf->msg->error() + << "Creating a PerturbedLennardJonesPotentialGPU with no GPU in the " + << "execution configuration" << std::endl; + throw std::runtime_error("Error initializing PerturbedLennardJonesEvapPotentialGPU"); + } + m_tuner.reset(new Autotuner<1>({AutotunerBase::makeBlockSizeRange(m_exec_conf)}, m_exec_conf, "perturbed_lennard_jones_evap")); @@ -41,23 +50,21 @@ PerturbedLennardJonesEvapGPU::PerturbedLennardJonesEvapGPU( void PerturbedLennardJonesEvapGPU::computeForces(uint64_t timestep) { - m_nlist->compute(timestep); - + this->m_nlist->compute(timestep); + this->m_force.zeroFill(); const Scalar interface_height = (*m_variant)(timestep); const Scalar scaled_t = scaleTime(timestep); - ArrayHandle d_pos(m_pdata->getPositions(), access_location::device, access_mode::read); + const BoxDim box = this->m_pdata->getGlobalBox(); + const unsigned int N = this->m_pdata->getN(); + const unsigned int n_ghost = this->m_pdata->getNGhosts(); + const unsigned int Ntot = N + n_ghost; - // neighbor-list arrays - ArrayHandle d_n_neigh(m_nlist->getNNeighArray(), - access_location::device, - access_mode::read); - ArrayHandle d_nlist(m_nlist->getNListArray(), - access_location::device, - access_mode::read); - ArrayHandle d_head_list(m_nlist->getHeadList(), - access_location::device, - access_mode::read); + if (m_scale_factor.getNumElements() < Ntot) + { + GPUArray tmp(Ntot, m_exec_conf); + m_scale_factor.swap(tmp); + } ArrayHandle d_data(m_attraction_scale_factor_data, access_location::device, @@ -67,27 +74,32 @@ void PerturbedLennardJonesEvapGPU::computeForces(uint64_t timestep) access_mode::read); ArrayHandle h_domain(m_domain, access_location::host, access_mode::read); - // build the interpolator: lo = {y_lo, t_lo}, hi = {y_hi, t_hi} const Scalar lo[2] = {Scalar(0.0), h_domain.data[0]}; const Scalar hi[2] = {Scalar(1.0), h_domain.data[1]}; - const LinearInterpolator2D interp(d_data.data, h_shape.data, lo, hi); + LinearInterpolator2D interp(d_data.data, h_shape.data, lo, hi); - const unsigned int N = this->m_pdata->getN(); - const unsigned int n_ghost = this->m_pdata->getNGhosts(); - - GPUArray scale_factor(N + n_ghost, m_exec_conf); - ArrayHandle d_scale_factor(scale_factor, + ArrayHandle d_pos(m_pdata->getPositions(), access_location::device, access_mode::read); + ArrayHandle d_n_neigh(m_nlist->getNNeighArray(), + access_location::device, + access_mode::read); + ArrayHandle d_nlist(m_nlist->getNListArray(), + access_location::device, + access_mode::read); + ArrayHandle d_head_list(m_nlist->getHeadList(), + access_location::device, + access_mode::read); + ArrayHandle d_scale_factor(m_scale_factor, access_location::device, - access_mode::overwrite); - + access_mode::readwrite); ArrayHandle d_force(m_force, access_location::device, access_mode::overwrite); - m_tuner->begin(); + const Scalar rcutsq = m_rcut * m_rcut; + this->m_tuner->begin(); gpu::perturbed_lennard_jones_evap_args_t args(d_force.data, d_pos.data, d_scale_factor.data, - this->m_pdata->getGlobalBox(), + box, N, n_ghost, d_n_neigh.data, @@ -96,19 +108,74 @@ void PerturbedLennardJonesEvapGPU::computeForces(uint64_t timestep) interp, scaled_t, interface_height, - lj1, - lj2, - epsilon_x_4, rcutsq, - rwcasq, + m_params, m_energy_shift, m_tuner->getParam()[0]); gpu::compute_perturbed_lennard_jones_evap_forces(args); if (this->m_exec_conf->isCUDAErrorCheckingEnabled()) CHECK_CUDA_ERROR(); - m_tuner->end(); + this->m_tuner->end(); + } +namespace py = pybind11; + +void export_PerturbedLennardJonesEvapGPU(py::module& m) + { + py::class_>(m, "PerturbedLennardJonesEvapGPU") + .def(py::init( + [](std::shared_ptr sysdef, + std::shared_ptr nlist, + Scalar rcut, + Scalar epsilon, + Scalar sigma, + Scalar scale_factor, + bool energy_shift, + py::array_t + attraction_scale_factor_data, + py::array_t + attraction_scale_factor_shape, + py::array_t domain, + std::shared_ptr variant) + { + if (attraction_scale_factor_shape.size() != 2) + throw std::runtime_error("lambda_shape must have 2 elements"); + if (domain.size() != 2) + throw std::runtime_error( + "domain must have 2 elements [y_lo, y_hi, t_lo, t_hi]"); + + const unsigned int* shape_ptr = attraction_scale_factor_shape.data(); + const Scalar* data_ptr = attraction_scale_factor_data.data(); + const Scalar* dom_ptr = domain.data(); + + if (attraction_scale_factor_data.size() + != static_cast(shape_ptr[0] * shape_ptr[1])) + throw std::runtime_error( + "attraction_scale_factor_data size does not match lambda_shape"); + + PairParametersPerturbedLennardJones params; + const Scalar sigma_2 = sigma * sigma; + const Scalar sigma_4 = sigma_2 * sigma_2; + params.sigma_6 = sigma_2 * sigma_4; + params.epsilon_x_4 = Scalar(4.0) * epsilon; + params.rwcasq = std::pow(Scalar(2.0), Scalar(1.0) / Scalar(3.0)) * sigma_2; + params.attraction_scale_factor = Scalar(0.0); + + return std::make_shared(sysdef, + nlist, + rcut, + scale_factor, + params, + energy_shift, + data_ptr, + shape_ptr, + dom_ptr, + variant); + })); } + } // end namespace detail } // end namespace azplugins } // end namespace hoomd diff --git a/src/PerturbedLennardJonesEvapGPU.cu b/src/PerturbedLennardJonesEvapGPU.cu index 1fcf8532..c033c3ef 100644 --- a/src/PerturbedLennardJonesEvapGPU.cu +++ b/src/PerturbedLennardJonesEvapGPU.cu @@ -2,25 +2,20 @@ // Copyright (c) 2021-2025, Auburn University // Part of azplugins, released under the BSD 3-Clause License. +#include "PairEvaluatorPerturbedLennardJones.h" #include "PerturbedLennardJonesEvapGPU.cuh" +#include "hoomd/WarpTools.cuh" namespace hoomd { namespace azplugins { +namespace detail + { namespace gpu { namespace kernel { -__device__ __inline__ Scalar clamp_scaled_y(Scalar y, Scalar height, Scalar ylo) - { - const Scalar s = (y - ylo) / (height - ylo); - if (s < Scalar(0.0)) - return Scalar(0.0); - if (s > Scalar(1.0)) - return Scalar(1.0); - return s; - } __global__ void compute_attraction_scale_factor(Scalar* d_scale_factor, const Scalar4* d_pos, @@ -35,26 +30,30 @@ __global__ void compute_attraction_scale_factor(Scalar* d_scale_factor, return; const Scalar y = __ldg(d_pos + i).y; - const Scalar scaled_pos_y = clamp_scaled_y(y, interface_height, y_lo); + Scalar scaled_pos_y = (y - y_lo) / (interface_height - y_lo); + + if (!(scaled_pos_y > Scalar(0.0))) + scaled_pos_y = Scalar(0.0); + if (scaled_pos_y > Scalar(1.0)) + scaled_pos_y = Scalar(1.0); + d_scale_factor[i] = interp(scaled_pos_y, scaled_t); } -__global__ void compute_perturbed_lennard_jones_evap_forces(Scalar4* d_force, - const Scalar4* d_pos, - const BoxDim& box, - const unsigned int N, - const unsigned int* d_n_neigh, - const unsigned int* d_nlist, - const size_t* d_head_list, - const Scalar* d_scale_factor, - const Scalar lj1, - const Scalar lj2, - const Scalar epsilon_x_4, - const Scalar rcutsq, - const Scalar rwcasq, - const bool energy_shift) +__global__ void +compute_perturbed_lennard_jones_evap_forces(Scalar4* d_force, + const Scalar4* d_pos, + Scalar* d_scale_factor, + const BoxDim box, + const unsigned int N, + const unsigned int n_ghost, + const unsigned int* d_n_neigh, + const unsigned int* d_nlist, + const size_t* d_head_list, + const Scalar rcutsq, + PairParametersPerturbedLennardJones params, + const bool energy_shift) { - // start by identifying which particle we are to handle const unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x; if (idx >= N) return; @@ -67,75 +66,36 @@ __global__ void compute_perturbed_lennard_jones_evap_forces(Scalar4* d_force, Scalar3 fi = make_scalar3(0, 0, 0); Scalar pei = 0; - // load in the length of the neighbor list for this thread const unsigned int n_neigh = d_n_neigh[idx]; const size_t head = d_head_list[idx]; - // loop over neighbors + // loop over this particle's neighbors for (unsigned int k = 0; k < n_neigh; ++k) { - const unsigned int j = d_nlist[head + k]; + const unsigned int j = __ldg(d_nlist + head + k); if (j == idx) continue; const Scalar4 postype_j = __ldg(d_pos + j); const Scalar3 pos_j = make_scalar3(postype_j.x, postype_j.y, postype_j.z); - const Scalar attraction_scale_factor_j = __ldg(d_scale_factor + j); // minimum-image Scalar3 dx = pos_i - pos_j; dx = box.minImage(dx); - const Scalar rsq = dot(dx, dx); - const Scalar attraction_scale_factor_avg - = Scalar(0.5) * (attraction_scale_factor_i + attraction_scale_factor_j); - const Scalar wca_shift - = epsilon_x_4 * (Scalar(1.0) - attraction_scale_factor_avg) / Scalar(4.0); + params.attraction_scale_factor + = Scalar(0.5) * (attraction_scale_factor_i + attraction_scale_factor_j); - if (rsq < rcutsq && lj1 != 0) + Scalar force_divr = Scalar(0.0); + Scalar pair_eng = Scalar(0.0); + PairEvaluatorPerturbedLennardJones eval(rsq, rcutsq, params); + if (eval.evalForceAndEnergy(force_divr, pair_eng, energy_shift)) { - const Scalar r2inv = Scalar(1) / rsq; - const Scalar r6inv = r2inv * r2inv * r2inv; - - Scalar pair_eng = r6inv * (lj1 * r6inv - lj2); - Scalar force_divr = r6inv * r2inv * (Scalar(12.0) * lj1 * r6inv - Scalar(6.0) * lj2); - - if (rsq < rwcasq) - { - pair_eng += wca_shift; - } - else - { - pair_eng *= attraction_scale_factor_avg; - force_divr *= attraction_scale_factor_avg; - } - - if (energy_shift) - { - const Scalar rcut2inv = Scalar(1.0) / rcutsq; - const Scalar rcut6inv = rcut2inv * rcut2inv * rcut2inv; - - Scalar pair_eng_shift = rcut6inv * (lj1 * rcut6inv - lj2); - - if (rcutsq < rwcasq) - { - pair_eng_shift += wca_shift; - } - else - { - pair_eng_shift *= attraction_scale_factor_avg; - } - - // apply the shift to the pair energy - pair_eng -= pair_eng_shift; - } - fi.x += force_divr * dx.x; fi.y += force_divr * dx.y; fi.z += force_divr * dx.z; - pei += pair_eng; } } @@ -148,21 +108,21 @@ __global__ void compute_perturbed_lennard_jones_evap_forces(Scalar4* d_force, hipError_t compute_perturbed_lennard_jones_evap_forces(const perturbed_lennard_jones_evap_args_t& args) { - unsigned int max_block_size; + assert(args.block_size != 0); + hipFuncAttributes attr; hipFuncGetAttributes( &attr, reinterpret_cast(&kernel::compute_perturbed_lennard_jones_evap_forces)); - max_block_size = attr.maxThreadsPerBlock; - + const unsigned max_block_size = attr.maxThreadsPerBlock; const unsigned int run_block_size = min(args.block_size, max_block_size); - const unsigned int Ntot = args.N + args.N_ghost; - dim3 grid(args.N / run_block_size + 1, 1, 1); - dim3 threads(run_block_size, 1, 1); + const unsigned int Ntot = args.N + args.n_ghost; + dim3 threads(run_block_size, 1, 1); + dim3 grid((Ntot + run_block_size - 1) / run_block_size, 1, 1); const Scalar ylo = args.box.getLo().y; - hipLaunchKernelGGL((azplugins::gpu::kernel::compute_attraction_scale_factor), + hipLaunchKernelGGL((gpu::kernel::compute_attraction_scale_factor), grid, threads, 0, @@ -175,29 +135,28 @@ compute_perturbed_lennard_jones_evap_forces(const perturbed_lennard_jones_evap_a args.scaled_t, ylo); - hipLaunchKernelGGL((azplugins::gpu::kernel::compute_perturbed_lennard_jones_evap_forces), + hipLaunchKernelGGL((gpu::kernel::compute_perturbed_lennard_jones_evap_forces), grid, threads, 0, 0, args.d_force, args.d_pos, + args.d_scale_factor, args.box, args.N, + args.n_ghost, args.d_n_neigh, args.d_nlist, args.d_head_list, - args.d_scale_factor, - args.lj1, - args.lj2, - args.epsilon_x_4, args.rcutsq, - args.rwcasq, + args.params, args.energy_shift); return hipSuccess; } } // end namespace gpu + } // end namespace detail } // end namespace azplugins } // end namespace hoomd diff --git a/src/PerturbedLennardJonesEvapGPU.cuh b/src/PerturbedLennardJonesEvapGPU.cuh index 4a5ed95a..d1327b99 100644 --- a/src/PerturbedLennardJonesEvapGPU.cuh +++ b/src/PerturbedLennardJonesEvapGPU.cuh @@ -13,13 +13,15 @@ #include "hoomd/TextureTools.h" #include "LinearInterpolator2D.h" - +#include "PairEvaluatorPerturbedLennardJones.h" #include namespace hoomd { namespace azplugins { +namespace detail + { namespace gpu { //! Wraps arguments to kernel driver @@ -28,7 +30,7 @@ struct perturbed_lennard_jones_evap_args_t perturbed_lennard_jones_evap_args_t(Scalar4* _d_force, const Scalar4* _d_pos, Scalar* _d_scale_factor, - const BoxDim& _box, + const BoxDim _box, const unsigned int _N, const unsigned int _n_ghost, const unsigned int* _d_n_neigh, @@ -37,36 +39,30 @@ struct perturbed_lennard_jones_evap_args_t const LinearInterpolator2D& _interp, const Scalar _scaled_t, const Scalar _interface_height, - const Scalar _lj1, - const Scalar _lj2, - const Scalar _epsilon_x_4, const Scalar _rcutsq, - const Scalar _rwcasq, + PairParametersPerturbedLennardJones _params, const bool _energy_shift, const unsigned int _block_size) : d_force(_d_force), d_pos(_d_pos), d_scale_factor(_d_scale_factor), box(_box), N(_N), - N_ghost(_n_ghost), d_n_neigh(_d_n_neigh), d_nlist(_d_nlist), d_head_list(_d_head_list), - interp(_interp), scaled_t(_scaled_t), interface_height(_interface_height), lj1(_lj1), - lj2(_lj2), epsilon_x_4(_epsilon_x_4), rcutsq(_rcutsq), rwcasq(_rwcasq), - energy_shift(_energy_shift), block_size(_block_size) { }; + n_ghost(_n_ghost), d_n_neigh(_d_n_neigh), d_nlist(_d_nlist), d_head_list(_d_head_list), + interp(_interp), scaled_t(_scaled_t), interface_height(_interface_height), + rcutsq(_rcutsq), params(_params), energy_shift(_energy_shift), + block_size(_block_size) { }; Scalar4* d_force; const Scalar4* d_pos; Scalar* d_scale_factor; const BoxDim box; const unsigned int N; - const unsigned int N_ghost; + const unsigned int n_ghost; const unsigned int* d_n_neigh; const unsigned int* d_nlist; const size_t* d_head_list; const LinearInterpolator2D interp; const Scalar scaled_t; const Scalar interface_height; - const Scalar lj1; - const Scalar lj2; - const Scalar epsilon_x_4; const Scalar rcutsq; - const Scalar rwcasq; + PairParametersPerturbedLennardJones params; const bool energy_shift; const unsigned int block_size; }; @@ -75,6 +71,7 @@ hipError_t compute_perturbed_lennard_jones_evap_forces(const perturbed_lennard_jones_evap_args_t& args); } // end namespace gpu + } // end namespace detail } // end namespace azplugins } // end namespace hoomd diff --git a/src/PerturbedLennardJonesEvapGPU.h b/src/PerturbedLennardJonesEvapGPU.h index 8b75d5e6..45a027c2 100644 --- a/src/PerturbedLennardJonesEvapGPU.h +++ b/src/PerturbedLennardJonesEvapGPU.h @@ -5,22 +5,29 @@ #ifndef AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_GPU_H_ #define AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_GPU_H_ -#ifdef __HIPCC__ -#error This header cannot be compiled by nvcc -#endif +#ifdef ENABLE_HIP #include "PerturbedLennardJonesEvap.h" #include "PerturbedLennardJonesEvapGPU.cuh" #include "hoomd/Autotuner.h" +#ifdef __HIPCC__ +#error This header cannot be compiled by nvcc +#endif + namespace hoomd { namespace azplugins { -class PYBIND11_EXPORT PerturbedLennardJonesEvapGPU : public PerturbedLennardJonesEvap +namespace detail + { +class PerturbedLennardJonesEvapGPU : public PerturbedLennardJonesEvap { public: //! Constructor + + typedef detail::PairParametersPerturbedLennardJones param_type; + PerturbedLennardJonesEvapGPU(std::shared_ptr sysdef, std::shared_ptr nlist, const Scalar r_cut, @@ -37,12 +44,11 @@ class PYBIND11_EXPORT PerturbedLennardJonesEvapGPU : public PerturbedLennardJone void computeForces(uint64_t timestep) override; }; -namespace detail - { void export_PerturbedLennardJonesEvapGPU(pybind11::module& m); } // end namespace detail } // end namespace azplugins } // end namespace hoomd +#endif // ENABLE_HIP #endif // AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_GPU_H_ diff --git a/src/export_PerturbedLennardJonesEvap.cc b/src/export_PerturbedLennardJonesEvap.cc deleted file mode 100644 index 95c2d870..00000000 --- a/src/export_PerturbedLennardJonesEvap.cc +++ /dev/null @@ -1,70 +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. - -#include "PerturbedLennardJonesEvap.h" - -#include - -namespace hoomd - { -namespace azplugins - { -namespace detail - { -namespace py = pybind11; - -void export_PerturbedLennardJonesEvap(py::module& m) - { - py::class_>( - m, - "PerturbedLennardJonesEvap") - .def(py::init( - [](std::shared_ptr sysdef, - std::shared_ptr nlist, - Scalar rcut, - Scalar epsilon, - Scalar sigma, - Scalar time_scale_factor, - bool energy_shift, - py::array_t - attraction_scale_factor_data, - py::array_t - attraction_scale_factor_shape, - py::array_t domain, - std::shared_ptr variant) - { - if (attraction_scale_factor_shape.size() != 2) - throw std::runtime_error("lambda_shape must have 2 elements"); - if (domain.size() != 2) - throw std::runtime_error("domain must have 2 elements [t_lo, t_hi]"); - - const unsigned int* shape_ptr = attraction_scale_factor_shape.data(); - const Scalar* data_ptr = attraction_scale_factor_data.data(); - const Scalar* dom_ptr = domain.data(); - - if (attraction_scale_factor_data.size() - != static_cast(shape_ptr[0] * shape_ptr[1])) - throw std::runtime_error( - "attraction_scale_factor_data size does not match lambda_shape"); - - PairParametersPerturbedLennardJonesEvap params(epsilon, sigma); - return std::make_shared(sysdef, - nlist, - rcut, - time_scale_factor, - params, - energy_shift, - data_ptr, - shape_ptr, - dom_ptr, - variant); - })) - .def_property_readonly("rcut", &PerturbedLennardJonesEvap::getRCut) - .def_property_readonly("epsilon", &PerturbedLennardJonesEvap::getEpsilon) - .def_property_readonly("sigma", &PerturbedLennardJonesEvap::getSigma); - } - - } // end namespace detail - } // end namespace azplugins - } // end namespace hoomd diff --git a/src/export_PerturbedLennardJonesEvapGPU.cc b/src/export_PerturbedLennardJonesEvapGPU.cc deleted file mode 100644 index 87699e29..00000000 --- a/src/export_PerturbedLennardJonesEvapGPU.cc +++ /dev/null @@ -1,70 +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. - -#include "PerturbedLennardJonesEvapGPU.h" - -#include - -namespace hoomd - { -namespace azplugins - { -namespace detail - { -namespace py = pybind11; - -void export_PerturbedLennardJonesEvapGPU(py::module& m) - { - py::class_>(m, "PerturbedLennardJonesEvapGPU") - .def(py::init( - [](std::shared_ptr sysdef, - std::shared_ptr nlist, - Scalar rcut, - Scalar epsilon, - Scalar sigma, - Scalar scale_factor, - bool energy_shift, - py::array_t - attraction_scale_factor_data, - py::array_t - attraction_scale_factor_shape, - py::array_t domain, - std::shared_ptr variant) - { - if (attraction_scale_factor_shape.size() != 2) - throw std::runtime_error("lambda_shape must have 2 elements"); - if (domain.size() != 2) - throw std::runtime_error( - "domain must have 2 elements [y_lo, y_hi, t_lo, t_hi]"); - - const unsigned int* shape_ptr = attraction_scale_factor_shape.data(); - const Scalar* data_ptr = attraction_scale_factor_data.data(); - const Scalar* dom_ptr = domain.data(); - - if (attraction_scale_factor_data.size() - != static_cast(shape_ptr[0] * shape_ptr[1])) - throw std::runtime_error( - "attraction_scale_factor_data size does not match lambda_shape"); - - PairParametersPerturbedLennardJonesEvap params(epsilon, sigma); - return std::make_shared(sysdef, - nlist, - rcut, - scale_factor, - params, - energy_shift, - data_ptr, - shape_ptr, - dom_ptr, - variant); - })) - .def_property_readonly("rcut", &PerturbedLennardJonesEvapGPU::getRCut) - .def_property_readonly("epsilon", &PerturbedLennardJonesEvapGPU::getEpsilon) - .def_property_readonly("sigma", &PerturbedLennardJonesEvapGPU::getSigma); - } - } // end namespace detail - } // end namespace azplugins - } // end namespace hoomd diff --git a/src/pair.py b/src/pair.py index 1fce2885..67d4a28c 100644 --- a/src/pair.py +++ b/src/pair.py @@ -5,7 +5,7 @@ """Pair potentials.""" import numpy - +import hoomd from hoomd.azplugins import _azplugins from hoomd.data.parameterdicts import ParameterDict, TypeParameterDict from hoomd.data.typeparam import TypeParameter @@ -444,7 +444,6 @@ def __init__( time_scale_factor, energy_shift, attraction_scale_factor_data, - attraction_scale_factor_shape, domain, variant, ): @@ -472,13 +471,17 @@ def __init__( attraction_scale_factor_data, dtype=numpy.float64 ) self._attraction_scale_factor_shape = numpy.asarray( - attraction_scale_factor_shape, dtype=numpy.uint32 - ) + attraction_scale_factor_data, dtype=numpy.uint64 + ).shape def _attach_hook(self): self._nlist._attach(self._simulation) - cls = getattr(self._ext_module, self._cpp_class_name) + if isinstance(self._simulation.device, hoomd.device.CPU): + cls = getattr(self._ext_module, self._cpp_class_name) + else: + cls = getattr(self._ext_module, self._cpp_class_name + "GPU") + self._cpp_obj = cls( self._simulation.state._cpp_sys_def, self._nlist._cpp_obj, diff --git a/src/pytest/test_perturbedljevap.py b/src/pytest/test_perturbedljevap.py index f447b43c..f198584a 100644 --- a/src/pytest/test_perturbedljevap.py +++ b/src/pytest/test_perturbedljevap.py @@ -35,7 +35,6 @@ def valid_args_const(): "time_scale_factor": 1.0, "energy_shift": False, "attraction_scale_factor_data": numpy.array([[0.6, 0.6], [0.6, 0.6]]), - "attraction_scale_factor_shape": [2, 2], "domain": [0.0, 100.0], "variant": hoomd.azplugins.variant.VariantInterpolated( [5.0, 4.0, 2.0, 1.0], 0, 300 @@ -58,10 +57,6 @@ def test_constructor(valid_args_const): evap._attraction_scale_factor_data, valid_args_const["attraction_scale_factor_data"], ) - numpy.testing.assert_array_equal( - evap._attraction_scale_factor_shape, - valid_args_const["attraction_scale_factor_shape"], - ) assert evap.rcut == 3.0 assert evap.time_scale_factor == 1.0 @@ -86,25 +81,6 @@ def test_domain_mismatch( sim.run(0) -def test_shape_mismatch( - valid_args_const, two_particle_snapshot_factory, simulation_factory -): - snap = two_particle_snapshot_factory - bad_args = valid_args_const.copy() - bad_args["attraction_scale_factor_shape"] = [2, 2, 2] - - sim = simulation_factory(snap) - - integrator = hoomd.md.Integrator(dt=0.005) - sim.operations.integrator = integrator - - evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**bad_args) - sim.operations.integrator.forces = [evap] - - with pytest.raises(RuntimeError): - sim.run(0) - - def test_data_size_mismatch( valid_args_const, two_particle_snapshot_factory, simulation_factory ): @@ -202,9 +178,8 @@ def valid_args_vary(request): "epsilon": 1.0, "sigma": 1.0, "time_scale_factor": time_scale_factor, - "energy_shift": False, + "energy_shift": True, "attraction_scale_factor_data": attraction_factor_table, - "attraction_scale_factor_shape": attraction_factor_table.shape, "domain": domain, "variant": hoomd.azplugins.variant.VariantInterpolated( [2, 0, -2, -4, -6, -8], @@ -236,7 +211,7 @@ def test_energy_and_force_calculation_vary( [0, 1.5150099394228083, 0.0], [0, -1.5150099394228083, 0.0], ] - expected_energies = [-0.30515561, -0.30515561] + expected_energies = [-0.3032789, -0.3032789] forces = evap.forces energies = evap.energies @@ -254,7 +229,7 @@ def test_energy_and_force_calculation_vary( [0.0, 0.24328626, 0.0], [0.0, -0.24328626, 0.0], ] - expected_energies = [-0.04900309, -0.04900309] + expected_energies = [-0.048701721, -0.048701721] forces = evap.forces energies = evap.energies From cb41127bf54b0c5036efec069479aaf736d495d1 Mon Sep 17 00:00:00 2001 From: mayukh33 Date: Thu, 2 Jul 2026 15:01:10 +0200 Subject: [PATCH 16/17] Add support for particle types to set epsilon and sigma and other fixes --- src/PerturbedLennardJonesEvap.cc | 94 ++++++++++++++++++++++------ src/PerturbedLennardJonesEvap.h | 22 ++++--- src/PerturbedLennardJonesEvapGPU.cc | 20 ++---- src/PerturbedLennardJonesEvapGPU.cu | 19 ++++-- src/PerturbedLennardJonesEvapGPU.cuh | 8 ++- src/PerturbedLennardJonesEvapGPU.h | 3 +- src/pair.py | 71 +++++++++++++++++---- src/pytest/test_perturbedljevap.py | 15 ++--- 8 files changed, 180 insertions(+), 72 deletions(-) diff --git a/src/PerturbedLennardJonesEvap.cc b/src/PerturbedLennardJonesEvap.cc index 99150076..881aa842 100644 --- a/src/PerturbedLennardJonesEvap.cc +++ b/src/PerturbedLennardJonesEvap.cc @@ -5,6 +5,7 @@ #include "PerturbedLennardJonesEvap.h" #include +#include namespace hoomd { @@ -17,14 +18,13 @@ PerturbedLennardJonesEvap::PerturbedLennardJonesEvap( std::shared_ptr nlist, const Scalar rcut, const Scalar time_scale_factor, - const param_type& params, bool energy_shift, const Scalar* attraction_scale_factor_data, const unsigned int* attraction_scale_factor_shape, const Scalar* domain, std::shared_ptr variant) : ForceCompute(sysdef), m_nlist(nlist), m_rcut(rcut), m_time_scale_factor(time_scale_factor), - m_params(params), m_energy_shift(energy_shift), m_variant(variant) + m_typpair_idx(m_pdata->getNTypes()), m_energy_shift(energy_shift), m_variant(variant) { // Allocate and fill the domain for time: [t_lo, t_hi] { @@ -61,9 +61,14 @@ PerturbedLennardJonesEvap::PerturbedLennardJonesEvap( } { - const Index2D typpair_idx(m_pdata->getNTypes()); + GPUArray params_arr(m_typpair_idx.getNumElements(), m_exec_conf); + m_params.swap(params_arr); + ArrayHandle h_params(m_params, access_location::host, access_mode::overwrite); + } + + { m_r_cut_nlist - = std::make_shared>(typpair_idx.getNumElements(), m_exec_conf); + = std::make_shared>(m_typpair_idx.getNumElements(), m_exec_conf); { ArrayHandle h_r_cut_nlist(*m_r_cut_nlist, access_location::host, @@ -77,6 +82,59 @@ PerturbedLennardJonesEvap::PerturbedLennardJonesEvap( } } +void PerturbedLennardJonesEvap::validateTypes(unsigned int typ1, + unsigned int typ2, + std::string action) const + { + const unsigned int n_types = m_pdata->getNTypes(); + if (typ1 >= n_types || typ2 >= n_types) + throw std::runtime_error("Invalid type encountered when " + action); + } + +void PerturbedLennardJonesEvap::setParams(unsigned int typ1, + unsigned int typ2, + const param_type& param) + { + validateTypes(typ1, typ2, "setting params"); + ArrayHandle h_params(m_params, access_location::host, access_mode::readwrite); + + h_params.data[m_typpair_idx(typ1, typ2)] = param; + h_params.data[m_typpair_idx(typ2, typ1)] = param; + } + +void PerturbedLennardJonesEvap::setParamsPython(pybind11::tuple typ, pybind11::dict params) + { + const auto typ1 = m_pdata->getTypeByName(typ[0].cast()); + const auto typ2 = m_pdata->getTypeByName(typ[1].cast()); + + const Scalar epsilon = params["epsilon"].cast(); + const Scalar sigma = params["sigma"].cast(); + const Scalar sigma_2 = sigma * sigma; + + param_type p; + p.sigma_6 = sigma_2 * sigma_2 * sigma_2; + p.epsilon_x_4 = Scalar(4.0) * epsilon; + p.rwcasq = std::pow(Scalar(2.0), Scalar(1.0) / Scalar(3.0)) * sigma_2; + p.attraction_scale_factor = Scalar(0.0); + + setParams(typ1, typ2, p); + } + +pybind11::dict PerturbedLennardJonesEvap::getParams(pybind11::tuple typ) + { + const auto typ1 = m_pdata->getTypeByName(typ[0].cast()); + const auto typ2 = m_pdata->getTypeByName(typ[1].cast()); + validateTypes(typ1, typ2, "getting params"); + + ArrayHandle h_params(m_params, access_location::host, access_mode::read); + const param_type& p = h_params.data[m_typpair_idx(typ1, typ2)]; + + pybind11::dict v; + v["epsilon"] = p.epsilon_x_4 / Scalar(4.0); + v["sigma"] = std::pow(p.sigma_6, Scalar(1.0) / Scalar(6.0)); + return v; + } + void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) { m_nlist->compute(timestep); @@ -84,7 +142,8 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) const bool third_law = (m_nlist->getStorageMode() == hoomd::md::NeighborList::storageMode::half); - Scalar scaled_t = scaleTime(timestep); + Scalar scaled_t = Scalar(static_cast(timestep) / m_time_scale_factor); + Scalar interface_height = (*m_variant)(timestep); ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::read); @@ -98,6 +157,7 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) access_location::host, access_mode::read); ArrayHandle h_domain(m_domain, access_location::host, access_mode::read); + ArrayHandle h_params(m_params, access_location::host, access_mode::read); // Neighbor-list arrays ArrayHandle h_n_neigh(m_nlist->getNNeighArray(), @@ -147,6 +207,7 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) for (unsigned int i = 0; i < N; ++i) { const Scalar3 pos_i = make_scalar3(h_pos.data[i].x, h_pos.data[i].y, h_pos.data[i].z); + const unsigned int type_i = __scalar_as_int(h_pos.data[i].w); Scalar3 fi = make_scalar3(0, 0, 0); Scalar pei = 0; @@ -162,6 +223,7 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) if (j == i) continue; Scalar3 pos_j = make_scalar3(h_pos.data[j].x, h_pos.data[j].y, h_pos.data[j].z); + const unsigned int type_j = __scalar_as_int(h_pos.data[j].w); const Scalar attraction_scale_factor_j = h_scale_factor.data[j]; @@ -170,14 +232,15 @@ void PerturbedLennardJonesEvap::computeForces(uint64_t timestep) dx = box.minImage(dx); const Scalar rsq = dot(dx, dx); - const Scalar attraction_scale_factor_avg - = Scalar(0.5) * (attraction_scale_factor_i + attraction_scale_factor_j); + param_type param = h_params.data[m_typpair_idx(type_i, type_j)]; - m_params.attraction_scale_factor = attraction_scale_factor_avg; + // Setting averaged attraction scale factor + param.attraction_scale_factor + = Scalar(0.5) * (attraction_scale_factor_i + attraction_scale_factor_j); Scalar force_divr = Scalar(0.0); Scalar pair_eng = Scalar(0.0); - PairEvaluatorPerturbedLennardJones eval(rsq, rcutsq, m_params); + PairEvaluatorPerturbedLennardJones eval(rsq, rcutsq, param); bool evaluated = eval.evalForceAndEnergy(force_divr, pair_eng, m_energy_shift); if (evaluated) @@ -216,8 +279,6 @@ void export_PerturbedLennardJonesEvap(py::module& m) [](std::shared_ptr sysdef, std::shared_ptr nlist, Scalar rcut, - Scalar epsilon, - Scalar sigma, Scalar time_scale_factor, bool energy_shift, py::array_t @@ -241,23 +302,18 @@ void export_PerturbedLennardJonesEvap(py::module& m) throw std::runtime_error( "attraction_scale_factor_data size does not match lambda_shape"); - PairParametersPerturbedLennardJones params; - const Scalar sigma_2 = sigma * sigma; - params.sigma_6 = sigma_2 * sigma_2 * sigma_2; - params.epsilon_x_4 = Scalar(4.0) * epsilon; - params.rwcasq = std::pow(Scalar(2.0), Scalar(1.0) / Scalar(3.0)) * sigma_2; - return std::make_shared(sysdef, nlist, rcut, time_scale_factor, - params, energy_shift, data_ptr, shape_ptr, dom_ptr, variant); - })); + })) + .def("setParams", &PerturbedLennardJonesEvap::setParamsPython) + .def("getParams", &PerturbedLennardJonesEvap::getParams); } } // end namespace detail diff --git a/src/PerturbedLennardJonesEvap.h b/src/PerturbedLennardJonesEvap.h index ffca7365..501cd01a 100644 --- a/src/PerturbedLennardJonesEvap.h +++ b/src/PerturbedLennardJonesEvap.h @@ -13,6 +13,7 @@ #include "hoomd/BoxDim.h" #include "hoomd/ForceCompute.h" #include "hoomd/GPUArray.h" +#include "hoomd/Index1D.h" #include "hoomd/Variant.h" #include "hoomd/VectorMath.h" #include "hoomd/md/NeighborList.h" @@ -30,13 +31,12 @@ namespace detail class PerturbedLennardJonesEvap : public ForceCompute { public: - typedef detail::PairParametersPerturbedLennardJones param_type; - + typedef PairParametersPerturbedLennardJones param_type; + //! Constructor PerturbedLennardJonesEvap(std::shared_ptr sysdef, std::shared_ptr nlist, - const Scalar r_cut, + const Scalar rcut, const Scalar time_scale_factor, - const param_type& params, bool energy_shift, const Scalar* attraction_scale_factor_data, const unsigned int* attraction_scale_factor_shape, @@ -50,17 +50,17 @@ class PerturbedLennardJonesEvap : public ForceCompute m_nlist->removeRCutMatrix(m_r_cut_nlist); } - Scalar scaleTime(uint64_t timestep) const - { - return Scalar(static_cast(timestep) / m_time_scale_factor); - } + //! Set and get parameters for sigma and epsilon for different types + void setParams(unsigned int typ1, unsigned int typ2, const param_type& param); + void setParamsPython(pybind11::tuple typ, pybind11::dict params); + pybind11::dict getParams(pybind11::tuple typ); protected: std::shared_ptr m_nlist; //!< Neighbor list - Scalar epsilon_x_4; Scalar m_rcut; Scalar m_time_scale_factor; //!< Time scaling factor - param_type m_params; + Index2D m_typpair_idx; //!< Indexer for the type-pair parameter table + GPUArray m_params; bool m_energy_shift; GPUArray m_domain; //!< [t_lo, t_hi] GPUArray m_attraction_scale_factor_data; //!< Flattened (y, t) data @@ -71,6 +71,8 @@ class PerturbedLennardJonesEvap : public ForceCompute std::shared_ptr> m_r_cut_nlist; //!< Cutoff matrix shared with the neighbor list + void validateTypes(unsigned int typ1, unsigned int typ2, std::string action) const; + void computeForces(uint64_t timestep) override; }; diff --git a/src/PerturbedLennardJonesEvapGPU.cc b/src/PerturbedLennardJonesEvapGPU.cc index c370f5ea..ea021f05 100644 --- a/src/PerturbedLennardJonesEvapGPU.cc +++ b/src/PerturbedLennardJonesEvapGPU.cc @@ -17,7 +17,6 @@ PerturbedLennardJonesEvapGPU::PerturbedLennardJonesEvapGPU( std::shared_ptr nlist, const Scalar r_cut, const Scalar scale_factor, - const param_type& params, bool energy_shift, const Scalar* attraction_scale_factor_data, const unsigned int* attraction_scale_factor_shape, @@ -27,7 +26,6 @@ PerturbedLennardJonesEvapGPU::PerturbedLennardJonesEvapGPU( nlist, r_cut, scale_factor, - params, energy_shift, attraction_scale_factor_data, attraction_scale_factor_shape, @@ -53,12 +51,13 @@ void PerturbedLennardJonesEvapGPU::computeForces(uint64_t timestep) this->m_nlist->compute(timestep); this->m_force.zeroFill(); const Scalar interface_height = (*m_variant)(timestep); - const Scalar scaled_t = scaleTime(timestep); + const Scalar scaled_t = Scalar(static_cast(timestep) / m_time_scale_factor); const BoxDim box = this->m_pdata->getGlobalBox(); const unsigned int N = this->m_pdata->getN(); const unsigned int n_ghost = this->m_pdata->getNGhosts(); const unsigned int Ntot = N + n_ghost; + const unsigned int ntypes = this->m_pdata->getNTypes(); if (m_scale_factor.getNumElements() < Ntot) { @@ -92,6 +91,7 @@ void PerturbedLennardJonesEvapGPU::computeForces(uint64_t timestep) ArrayHandle d_scale_factor(m_scale_factor, access_location::device, access_mode::readwrite); + ArrayHandle d_params(m_params, access_location::device, access_mode::read); ArrayHandle d_force(m_force, access_location::device, access_mode::overwrite); const Scalar rcutsq = m_rcut * m_rcut; @@ -109,7 +109,8 @@ void PerturbedLennardJonesEvapGPU::computeForces(uint64_t timestep) scaled_t, interface_height, rcutsq, - m_params, + d_params.data, + ntypes, m_energy_shift, m_tuner->getParam()[0]); @@ -129,8 +130,6 @@ void export_PerturbedLennardJonesEvapGPU(py::module& m) [](std::shared_ptr sysdef, std::shared_ptr nlist, Scalar rcut, - Scalar epsilon, - Scalar sigma, Scalar scale_factor, bool energy_shift, py::array_t @@ -155,19 +154,10 @@ void export_PerturbedLennardJonesEvapGPU(py::module& m) throw std::runtime_error( "attraction_scale_factor_data size does not match lambda_shape"); - PairParametersPerturbedLennardJones params; - const Scalar sigma_2 = sigma * sigma; - const Scalar sigma_4 = sigma_2 * sigma_2; - params.sigma_6 = sigma_2 * sigma_4; - params.epsilon_x_4 = Scalar(4.0) * epsilon; - params.rwcasq = std::pow(Scalar(2.0), Scalar(1.0) / Scalar(3.0)) * sigma_2; - params.attraction_scale_factor = Scalar(0.0); - return std::make_shared(sysdef, nlist, rcut, scale_factor, - params, energy_shift, data_ptr, shape_ptr, diff --git a/src/PerturbedLennardJonesEvapGPU.cu b/src/PerturbedLennardJonesEvapGPU.cu index c033c3ef..14b2fd97 100644 --- a/src/PerturbedLennardJonesEvapGPU.cu +++ b/src/PerturbedLennardJonesEvapGPU.cu @@ -4,7 +4,6 @@ #include "PairEvaluatorPerturbedLennardJones.h" #include "PerturbedLennardJonesEvapGPU.cuh" -#include "hoomd/WarpTools.cuh" namespace hoomd { @@ -32,7 +31,7 @@ __global__ void compute_attraction_scale_factor(Scalar* d_scale_factor, const Scalar y = __ldg(d_pos + i).y; Scalar scaled_pos_y = (y - y_lo) / (interface_height - y_lo); - if (!(scaled_pos_y > Scalar(0.0))) + if (scaled_pos_y < Scalar(0.0)) scaled_pos_y = Scalar(0.0); if (scaled_pos_y > Scalar(1.0)) scaled_pos_y = Scalar(1.0); @@ -51,15 +50,19 @@ compute_perturbed_lennard_jones_evap_forces(Scalar4* d_force, const unsigned int* d_nlist, const size_t* d_head_list, const Scalar rcutsq, - PairParametersPerturbedLennardJones params, + const PairParametersPerturbedLennardJones* d_params, + const unsigned int ntypes, const bool energy_shift) { const unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x; if (idx >= N) return; + const Index2D typpair_idx(ntypes); + const Scalar4 postype_i = __ldg(d_pos + idx); const Scalar3 pos_i = make_scalar3(postype_i.x, postype_i.y, postype_i.z); + const unsigned int type_i = __scalar_as_int(postype_i.w); const Scalar attraction_scale_factor_i = __ldg(d_scale_factor + idx); // initialize the force and energy to 0 @@ -78,6 +81,7 @@ compute_perturbed_lennard_jones_evap_forces(Scalar4* d_force, const Scalar4 postype_j = __ldg(d_pos + j); const Scalar3 pos_j = make_scalar3(postype_j.x, postype_j.y, postype_j.z); + const unsigned int type_j = __scalar_as_int(postype_j.w); const Scalar attraction_scale_factor_j = __ldg(d_scale_factor + j); // minimum-image @@ -85,13 +89,17 @@ compute_perturbed_lennard_jones_evap_forces(Scalar4* d_force, dx = box.minImage(dx); const Scalar rsq = dot(dx, dx); + PairParametersPerturbedLennardJones params = d_params[typpair_idx(type_i, type_j)]; + // Setting averaged attraction scale factor params.attraction_scale_factor = Scalar(0.5) * (attraction_scale_factor_i + attraction_scale_factor_j); Scalar force_divr = Scalar(0.0); Scalar pair_eng = Scalar(0.0); PairEvaluatorPerturbedLennardJones eval(rsq, rcutsq, params); - if (eval.evalForceAndEnergy(force_divr, pair_eng, energy_shift)) + bool evaluated = eval.evalForceAndEnergy(force_divr, pair_eng, energy_shift); + + if (evaluated) { fi.x += force_divr * dx.x; fi.y += force_divr * dx.y; @@ -150,7 +158,8 @@ compute_perturbed_lennard_jones_evap_forces(const perturbed_lennard_jones_evap_a args.d_nlist, args.d_head_list, args.rcutsq, - args.params, + args.d_params, + args.ntypes, args.energy_shift); return hipSuccess; diff --git a/src/PerturbedLennardJonesEvapGPU.cuh b/src/PerturbedLennardJonesEvapGPU.cuh index d1327b99..57f3563b 100644 --- a/src/PerturbedLennardJonesEvapGPU.cuh +++ b/src/PerturbedLennardJonesEvapGPU.cuh @@ -40,13 +40,14 @@ struct perturbed_lennard_jones_evap_args_t const Scalar _scaled_t, const Scalar _interface_height, const Scalar _rcutsq, - PairParametersPerturbedLennardJones _params, + const PairParametersPerturbedLennardJones* _d_params, + const unsigned int _ntypes, const bool _energy_shift, const unsigned int _block_size) : d_force(_d_force), d_pos(_d_pos), d_scale_factor(_d_scale_factor), box(_box), N(_N), n_ghost(_n_ghost), d_n_neigh(_d_n_neigh), d_nlist(_d_nlist), d_head_list(_d_head_list), interp(_interp), scaled_t(_scaled_t), interface_height(_interface_height), - rcutsq(_rcutsq), params(_params), energy_shift(_energy_shift), + rcutsq(_rcutsq), d_params(_d_params), ntypes(_ntypes), energy_shift(_energy_shift), block_size(_block_size) { }; Scalar4* d_force; @@ -62,7 +63,8 @@ struct perturbed_lennard_jones_evap_args_t const Scalar scaled_t; const Scalar interface_height; const Scalar rcutsq; - PairParametersPerturbedLennardJones params; + const PairParametersPerturbedLennardJones* d_params; //!< Per type-pair parameter table + const unsigned int ntypes; //!< Number of particle types const bool energy_shift; const unsigned int block_size; }; diff --git a/src/PerturbedLennardJonesEvapGPU.h b/src/PerturbedLennardJonesEvapGPU.h index 45a027c2..5b777b0c 100644 --- a/src/PerturbedLennardJonesEvapGPU.h +++ b/src/PerturbedLennardJonesEvapGPU.h @@ -26,13 +26,12 @@ class PerturbedLennardJonesEvapGPU : public PerturbedLennardJonesEvap public: //! Constructor - typedef detail::PairParametersPerturbedLennardJones param_type; + typedef PairParametersPerturbedLennardJones param_type; PerturbedLennardJonesEvapGPU(std::shared_ptr sysdef, std::shared_ptr nlist, const Scalar r_cut, const Scalar scale_factor, - const param_type& params, bool energy_shift, const Scalar* attraction_scale_factor_data, const unsigned int* attraction_scale_factor_shape, diff --git a/src/pair.py b/src/pair.py index 67d4a28c..06e648ca 100644 --- a/src/pair.py +++ b/src/pair.py @@ -430,7 +430,57 @@ def __init__(self, nlist, default_r_cut=None, default_r_on=0, mode="none"): class PerturbedLennardJonesEvap(Force): - """Perturbed Lennard Jones potential for evaporation""" + r"""Perturbed Lennard-Jones potential for evaporation. + + Args: + nlist (hoomd.md.nlist.NeighborList): Neighbor list. + rcut (float): Cutoff radius :math:`[\mathrm{length}]`. + time_scale_factor (float): Factor used to rescale the timestep when + sampling the attraction scale-factor table. + energy_shift (bool): If ``True``, shift the potential to zero at ``rcut``. + attraction_scale_factor_data (numpy.ndarray): 2D ``(ny, nt)`` table of + the attraction scale factor :math:`\lambda` as a function of scaled + y coordinate and scaled time. + domain (numpy.ndarray): ``[t_lo, t_hi]`` time domain of the table. + variant (hoomd.azplugins.variant.VariantInterpolated): Interface height + as a function of timestep. + + This works similar to the standard perturbed Lennard Jones pair potential, + except the attraction scale factor :math:`\lambda` is not constant. + It is computed for every time step and depends on the particle's y coordinate. + This emulates the evaporation of the bad solvent and its concentration + gradient as the air-solvent interface recedes. For now, types are added so that + epsilon can be set to zero for other particle types for which this potential is + not valid. + + Example:: + + nl = hoomd.md.nlist.Cell(buffer=0.4) + plj_evap = azplugins.pair.PerturbedLennardJonesEvap( + nlist=nl, + rcut=3.0, + time_scale_factor=1.0, + energy_shift=True, + attraction_scale_factor_data=lambda_table, + domain=[0.0, t_max], + variant=interface_variant, + ) + plj_evap.params[("A", "A")] = dict(epsilon=1.0, sigma=1.0) + plj_evap.params[("A", "B")] = dict(epsilon=0.0, sigma=0.0) + plj_evap.params[("B", "B")] = dict(epsilon=0.0, sigma=0.0) + + .. py:attribute:: params + + The potential parameters. The dictionary has the following keys: + + * ``epsilon`` (`float`, **required**) - energy parameter + :math:`\varepsilon` :math:`[\mathrm{energy}]` + * ``sigma`` (`float`, **required**) - particle size :math:`\sigma` + :math:`[\mathrm{length}]` + + Type: :class:`~hoomd.data.typeparam.TypeParameter` [`tuple` + [``particle_type``, ``particle_type``], `dict`] + """ _ext_module = _azplugins _cpp_class_name = "PerturbedLennardJonesEvap" @@ -439,8 +489,6 @@ def __init__( self, nlist, rcut, - epsilon, - sigma, time_scale_factor, energy_shift, attraction_scale_factor_data, @@ -452,16 +500,19 @@ def __init__( self._nlist = nlist self._variant = variant + params = TypeParameter( + "params", + "particle_types", + TypeParameterDict(epsilon=float, sigma=float, len_keys=2), + ) + self._add_typeparam(params) + param_dict = ParameterDict( rcut=float, - epsilon=float, - sigma=float, time_scale_factor=float, energy_shift=bool, ) param_dict["rcut"] = float(rcut) - param_dict["epsilon"] = float(epsilon) - param_dict["sigma"] = float(sigma) param_dict["time_scale_factor"] = float(time_scale_factor) param_dict["energy_shift"] = bool(energy_shift) self._param_dict.update(param_dict) @@ -471,8 +522,8 @@ def __init__( attraction_scale_factor_data, dtype=numpy.float64 ) self._attraction_scale_factor_shape = numpy.asarray( - attraction_scale_factor_data, dtype=numpy.uint64 - ).shape + self._attraction_scale_factor_data.shape, dtype=numpy.uint64 + ) def _attach_hook(self): self._nlist._attach(self._simulation) @@ -486,8 +537,6 @@ def _attach_hook(self): self._simulation.state._cpp_sys_def, self._nlist._cpp_obj, self.rcut, - self.epsilon, - self.sigma, self.time_scale_factor, self.energy_shift, self._attraction_scale_factor_data, diff --git a/src/pytest/test_perturbedljevap.py b/src/pytest/test_perturbedljevap.py index f198584a..dbc8810f 100644 --- a/src/pytest/test_perturbedljevap.py +++ b/src/pytest/test_perturbedljevap.py @@ -30,8 +30,6 @@ def valid_args_const(): return { "nlist": hoomd.md.nlist.Cell(buffer=0.4), "rcut": 3.0, - "epsilon": 1.0, - "sigma": 1.0, "time_scale_factor": 1.0, "energy_shift": False, "attraction_scale_factor_data": numpy.array([[0.6, 0.6], [0.6, 0.6]]), @@ -44,9 +42,10 @@ def valid_args_const(): def test_constructor(valid_args_const): evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**valid_args_const) + evap.params[("A", "A")] = dict(epsilon=1.0, sigma=1.0) - assert evap.epsilon == 1.0 - assert evap.sigma == 1.0 + assert evap.params[("A", "A")]["epsilon"] == 1.0 + assert evap.params[("A", "A")]["sigma"] == 1.0 assert evap._nlist is valid_args_const["nlist"] assert evap._variant is valid_args_const["variant"] @@ -75,6 +74,7 @@ def test_domain_mismatch( sim.operations.integrator = integrator evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**bad_args) + evap.params[("A", "A")] = dict(epsilon=1.0, sigma=1.0) sim.operations.integrator.forces = [evap] with pytest.raises(RuntimeError): @@ -94,6 +94,7 @@ def test_data_size_mismatch( sim.operations.integrator = integrator evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**bad_args) + evap.params[("A", "A")] = dict(epsilon=1.0, sigma=1.0) sim.operations.integrator.forces = [evap] with pytest.raises(RuntimeError): @@ -115,6 +116,7 @@ def test_variant_mismatch( sim.operations.integrator = integrator evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**bad_args) + evap.params[("A", "A")] = dict(epsilon=1.0, sigma=1.0) sim.operations.integrator.forces = [evap] with pytest.raises(TypeError): @@ -131,7 +133,7 @@ def test_energy_and_force_calculation_const( ): snap = two_particle_snapshot_factory evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**valid_args_const) - + evap.params[("A", "A")] = dict(epsilon=1.0, sigma=1.0) sim = simulation_factory(snap) integrator = hoomd.md.Integrator(dt=0.005) @@ -175,8 +177,6 @@ def valid_args_vary(request): return { "nlist": hoomd.md.nlist.Cell(buffer=0.4), "rcut": 3.0, - "epsilon": 1.0, - "sigma": 1.0, "time_scale_factor": time_scale_factor, "energy_shift": True, "attraction_scale_factor_data": attraction_factor_table, @@ -200,6 +200,7 @@ def test_energy_and_force_calculation_vary( """Energies/forces at t=0 and t=0.5, with and without time scaling.""" snap = two_particle_snapshot_factory evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**valid_args_vary) + evap.params[("A", "A")] = dict(epsilon=1.0, sigma=1.0) sim = simulation_factory(snap) integrator = hoomd.md.Integrator(dt=0.005) From fdf388c7819a2791be980913de508a447e9983eb Mon Sep 17 00:00:00 2001 From: mayukh33 Date: Mon, 6 Jul 2026 18:54:29 +0200 Subject: [PATCH 17/17] Add more tests and documentation --- src/pair.py | 14 ++-- src/pytest/test_perturbedljevap.py | 125 +++++++++++++++++++---------- 2 files changed, 90 insertions(+), 49 deletions(-) diff --git a/src/pair.py b/src/pair.py index 06e648ca..af0254a5 100644 --- a/src/pair.py +++ b/src/pair.py @@ -436,7 +436,7 @@ class PerturbedLennardJonesEvap(Force): nlist (hoomd.md.nlist.NeighborList): Neighbor list. rcut (float): Cutoff radius :math:`[\mathrm{length}]`. time_scale_factor (float): Factor used to rescale the timestep when - sampling the attraction scale-factor table. + interpolating the attraction scale-factor table. energy_shift (bool): If ``True``, shift the potential to zero at ``rcut``. attraction_scale_factor_data (numpy.ndarray): 2D ``(ny, nt)`` table of the attraction scale factor :math:`\lambda` as a function of scaled @@ -447,11 +447,13 @@ class PerturbedLennardJonesEvap(Force): This works similar to the standard perturbed Lennard Jones pair potential, except the attraction scale factor :math:`\lambda` is not constant. - It is computed for every time step and depends on the particle's y coordinate. - This emulates the evaporation of the bad solvent and its concentration - gradient as the air-solvent interface recedes. For now, types are added so that - epsilon can be set to zero for other particle types for which this potential is - not valid. + The attraction scale factor is determined by interpolating a regular + two-dimesional array where the rows define the y coordinates and the columns + define time coordinates. It is computed for every time step and depends on + the particle's y coordinate. This emulates the evaporation of the bad solvent + and its concentration gradient as the air-solvent interface recedes. For now, + types are added so that epsilon can be set to zero for other particle types + for which this potential is not valid. Example:: diff --git a/src/pytest/test_perturbedljevap.py b/src/pytest/test_perturbedljevap.py index dbc8810f..ca4a405d 100644 --- a/src/pytest/test_perturbedljevap.py +++ b/src/pytest/test_perturbedljevap.py @@ -13,16 +13,17 @@ @pytest.fixture def two_particle_snapshot_factory(): - snap = hoomd.Snapshot() - if snap.communicator.rank == 0: - snap.configuration.box = [20, 20, 20, 0, 0, 0] - snap.particles.N = 2 - snap.particles.types = ["A"] - snap.particles.position[:] = [ - [-10, -10, -10], - [-10, -8.8, -10], - ] - return snap + + def make_snapshot(positions): + snap = hoomd.Snapshot() + if snap.communicator.rank == 0: + snap.configuration.box = [20, 20, 20, 0, 0, 0] + snap.particles.N = 2 + snap.particles.types = ["A"] + snap.particles.position[:] = positions + return snap + + return make_snapshot @pytest.fixture @@ -64,7 +65,7 @@ def test_constructor(valid_args_const): def test_domain_mismatch( valid_args_const, two_particle_snapshot_factory, simulation_factory ): - snap = two_particle_snapshot_factory + snap = two_particle_snapshot_factory([[0, 0, 0.1], [0, 0, 0.1]]) sim = simulation_factory(snap) bad_args = valid_args_const.copy() @@ -84,7 +85,8 @@ def test_domain_mismatch( def test_data_size_mismatch( valid_args_const, two_particle_snapshot_factory, simulation_factory ): - snap = two_particle_snapshot_factory + snap = two_particle_snapshot_factory([[0, 0, 0.1], [0, 0, 0.1]]) + bad_args = valid_args_const.copy() bad_args["attraction_scale_factor_data"] = [1.0, 1.0] @@ -104,7 +106,8 @@ def test_data_size_mismatch( def test_variant_mismatch( valid_args_const, two_particle_snapshot_factory, simulation_factory ): - snap = two_particle_snapshot_factory + snap = two_particle_snapshot_factory([[0, 0, 0.1], [0, 0, 0.1]]) + bad_args = valid_args_const.copy() bad_args["variant"] = hoomd.variant.Constant( 1.0 @@ -128,10 +131,27 @@ def test_variant_mismatch( """ +@pytest.mark.parametrize( + "positions, epsilon, expected_forces, expected_energies", + [ + ( + [[-10, -10, -10], [-10, -8.8, -10]], + 1.0, + [[0.0, 1.32701601, 0.0], [0.0, -1.32701601, 0.0]], + [-0.267289586, -0.267289586], + ), + ], +) def test_energy_and_force_calculation_const( - valid_args_const, two_particle_snapshot_factory, simulation_factory + valid_args_const, + two_particle_snapshot_factory, + simulation_factory, + positions, + epsilon, + expected_forces, + expected_energies, ): - snap = two_particle_snapshot_factory + snap = two_particle_snapshot_factory(positions) evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**valid_args_const) evap.params[("A", "A")] = dict(epsilon=1.0, sigma=1.0) sim = simulation_factory(snap) @@ -143,9 +163,6 @@ def test_energy_and_force_calculation_const( sim.run(0) - expected_forces = [[0.0, 1.32701601, 0.0], [0.0, -1.32701601, 0.0]] - expected_energies = [-0.267289586, -0.267289586] - forces = evap.forces energies = evap.energies @@ -178,7 +195,7 @@ def valid_args_vary(request): "nlist": hoomd.md.nlist.Cell(buffer=0.4), "rcut": 3.0, "time_scale_factor": time_scale_factor, - "energy_shift": True, + "energy_shift": False, "attraction_scale_factor_data": attraction_factor_table, "domain": domain, "variant": hoomd.azplugins.variant.VariantInterpolated( @@ -190,17 +207,52 @@ def valid_args_vary(request): """Test energy and force calculation for single particle type -with varying attraction scale factor. -""" - - +with varying attraction scale factor.""" + + +@pytest.mark.parametrize( + ( + "positions", + "epsilon", + "expected_forces_t0", + "expected_energies_t0", + "expected_forces_t100", + "expected_energies_t100", + ), + [ + ( + [[-10, -10, -10], [-10, -8.8, -10]], + 1.0, + [[0, 1.5150099394228083, 0.0], [0, -1.5150099394228083, 0.0]], + [-0.3051556, -0.3051556], + [[0.0, 0.24328626, 0.0], [0.0, -0.24328626, 0.0]], + [-0.04900309, -0.04900309], + ), + ( + [[-0.5, -0.5, 0.1], [0.5, 0.5, 0.1]], + 2, + [[1.0125, 1.0125, 0.0], [-1.0125, -1.0125, 0.0]], + [-0.196875, -0.196875], + [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], + [0.0, 0.0], + ), + ], +) def test_energy_and_force_calculation_vary( - valid_args_vary, two_particle_snapshot_factory, simulation_factory + valid_args_vary, + two_particle_snapshot_factory, + simulation_factory, + positions, + epsilon, + expected_forces_t0, + expected_energies_t0, + expected_forces_t100, + expected_energies_t100, ): """Energies/forces at t=0 and t=0.5, with and without time scaling.""" - snap = two_particle_snapshot_factory + snap = two_particle_snapshot_factory(positions) evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**valid_args_vary) - evap.params[("A", "A")] = dict(epsilon=1.0, sigma=1.0) + evap.params[("A", "A")] = dict(epsilon=epsilon, sigma=1.0) sim = simulation_factory(snap) integrator = hoomd.md.Integrator(dt=0.005) @@ -208,33 +260,20 @@ def test_energy_and_force_calculation_vary( sim.operations.integrator = integrator sim.run(0) - expected_forces = [ - [0, 1.5150099394228083, 0.0], - [0, -1.5150099394228083, 0.0], - ] - expected_energies = [-0.3032789, -0.3032789] - forces = evap.forces energies = evap.energies if sim.device.communicator.rank == 0: - numpy.testing.assert_allclose(forces, expected_forces) - numpy.testing.assert_allclose(energies, expected_energies) + numpy.testing.assert_allclose(forces, expected_forces_t0) + numpy.testing.assert_allclose(energies, expected_energies_t0) """Test if the potential energy and forces change as expected after a certain time. """ sim.run(100) - - expected_forces = [ - [0.0, 0.24328626, 0.0], - [0.0, -0.24328626, 0.0], - ] - expected_energies = [-0.048701721, -0.048701721] - forces = evap.forces energies = evap.energies if sim.device.communicator.rank == 0: - numpy.testing.assert_allclose(forces, expected_forces) - numpy.testing.assert_allclose(energies, expected_energies) + numpy.testing.assert_allclose(forces, expected_forces_t100) + numpy.testing.assert_allclose(energies, expected_energies_t100)