From f4aabd3cde47e67af4d0754e89bdfae4c3bffa67 Mon Sep 17 00:00:00 2001 From: mayukh33 Date: Thu, 7 May 2026 17:45:36 +0200 Subject: [PATCH 1/4] 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 0000000..2906f88 --- /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 2/4] 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 0000000..39db8f2 --- /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 3/4] 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 2906f88..1d1e853 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 39db8f2..ad34461 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 4/4] 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 a5529d9..7f0ea7e 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 0000000..5ed7dc7 --- /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 1d1e853..e8de081 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 ad34461..b3e1ed0 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 0000000..9c3455e --- /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 3e71abf..70dbb57 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 0000000..e47d84c --- /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 0000000..350705a --- /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 1ab81c7..4ffc922 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 fbefb07..48c5b2e 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 0000000..9b80bd9 --- /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) + )