diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a5529d98..a3c1b1b3 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 + PerturbedLennardJonesEvap.cc + export_VariantInterpolated.cc module.cc ParabolicFlow.cc VelocityCompute.cc @@ -12,6 +14,7 @@ set(_${COMPONENT_NAME}_sources if (ENABLE_HIP) list(APPEND _${COMPONENT_NAME}_sources export_ImagePotentialBondHarmonicGPU.cc + PerturbedLennardJonesEvapGPU.cc VelocityComputeGPU.cc ) endif() @@ -19,6 +22,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 ) @@ -32,6 +36,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..64f059b4 --- /dev/null +++ b/src/LinearInterpolator1D.h @@ -0,0 +1,93 @@ +// 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 on a uniform grid + 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 new file mode 100644 index 00000000..937f452a --- /dev/null +++ b/src/LinearInterpolator2D.h @@ -0,0 +1,131 @@ +// 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_ + +#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); + } + } + + //! 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 + AZPLUGINS_HOSTDEVICE AZPLUGINS_FORCEINLINE Scalar getLo(int dim) const + { + return m_lo[dim]; + } + + //! Upper bound + 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_ diff --git a/src/PerturbedLennardJonesEvap.cc b/src/PerturbedLennardJonesEvap.cc new file mode 100644 index 00000000..881aa842 --- /dev/null +++ b/src/PerturbedLennardJonesEvap.cc @@ -0,0 +1,321 @@ +// 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 +#include + +namespace hoomd + { +namespace azplugins + { +namespace detail + { +PerturbedLennardJonesEvap::PerturbedLennardJonesEvap( + std::shared_ptr sysdef, + std::shared_ptr nlist, + const Scalar rcut, + const Scalar time_scale_factor, + 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_typpair_idx(m_pdata->getNTypes()), m_energy_shift(energy_shift), m_variant(variant) + { + // Allocate and fill the domain for time: [t_lo, t_hi] + { + 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 + 2, 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); + } + + { + 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>(m_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::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); + + const bool third_law + = (m_nlist->getStorageMode() == hoomd::md::NeighborList::storageMode::half); + + 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); + 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); + ArrayHandle h_params(m_params, 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] = {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); + + const BoxDim box = m_pdata->getGlobalBox(); + const unsigned int N = m_pdata->getN(); + const unsigned int Ntot = N + m_pdata->getNGhosts(); + + 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 < Ntot; k++) + { + 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); + const unsigned int type_i = __scalar_as_int(h_pos.data[i].w); + + Scalar3 fi = make_scalar3(0, 0, 0); + Scalar pei = 0; + + 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]; + + for (unsigned int k = 0; k < size; ++k) + { + const unsigned int j = h_nlist.data[head + k]; + 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]; + + // Minimum-image + Scalar3 dx = pos_i - pos_j; + dx = box.minImage(dx); + + const Scalar rsq = dot(dx, dx); + param_type param = h_params.data[m_typpair_idx(type_i, type_j)]; + + // 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, param); + 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; + + 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 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 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"); + + return std::make_shared(sysdef, + nlist, + rcut, + time_scale_factor, + energy_shift, + data_ptr, + shape_ptr, + dom_ptr, + variant); + })) + .def("setParams", &PerturbedLennardJonesEvap::setParamsPython) + .def("getParams", &PerturbedLennardJonesEvap::getParams); + } + + } // end namespace detail + } // end namespace azplugins + } // end namespace hoomd diff --git a/src/PerturbedLennardJonesEvap.h b/src/PerturbedLennardJonesEvap.h new file mode 100644 index 00000000..501cd01a --- /dev/null +++ b/src/PerturbedLennardJonesEvap.h @@ -0,0 +1,85 @@ +// 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 + +#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" + +#include "LinearInterpolator2D.h" +#include "PairEvaluatorPerturbedLennardJones.h" +#include "VariantInterpolated.h" + +namespace hoomd + { +namespace azplugins + { +namespace detail + { +class PerturbedLennardJonesEvap : public ForceCompute + { + public: + typedef PairParametersPerturbedLennardJones param_type; + //! Constructor + PerturbedLennardJonesEvap(std::shared_ptr sysdef, + std::shared_ptr nlist, + const Scalar rcut, + const Scalar time_scale_factor, + bool energy_shift, + const Scalar* attraction_scale_factor_data, + const unsigned int* attraction_scale_factor_shape, + const Scalar* domain, + std::shared_ptr variant); + + //! Destructor + ~PerturbedLennardJonesEvap() + { + if (m_r_cut_nlist) + m_nlist->removeRCutMatrix(m_r_cut_nlist); + } + + //! 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 m_rcut; + Scalar m_time_scale_factor; //!< Time scaling factor + 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 + GPUArray m_attraction_scale_factor_shape; //!< [ny, nt] + GPUArray m_scale_factor; + std::shared_ptr m_variant; + + 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; + }; + +void export_PerturbedLennardJonesEvap(pybind11::module& m); + } // end namespace detail + + } // end namespace azplugins + } // end namespace hoomd + +#endif // AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_H_ diff --git a/src/PerturbedLennardJonesEvapGPU.cc b/src/PerturbedLennardJonesEvapGPU.cc new file mode 100644 index 00000000..ea021f05 --- /dev/null +++ b/src/PerturbedLennardJonesEvapGPU.cc @@ -0,0 +1,171 @@ +// 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 + { +PerturbedLennardJonesEvapGPU::PerturbedLennardJonesEvapGPU( + std::shared_ptr sysdef, + std::shared_ptr nlist, + const Scalar r_cut, + const Scalar scale_factor, + 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, + energy_shift, + attraction_scale_factor_data, + attraction_scale_factor_shape, + 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")); + this->m_autotuners.push_back(m_tuner); + } + +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 = 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) + { + GPUArray tmp(Ntot, m_exec_conf); + m_scale_factor.swap(tmp); + } + + 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); + + const Scalar lo[2] = {Scalar(0.0), h_domain.data[0]}; + const Scalar hi[2] = {Scalar(1.0), h_domain.data[1]}; + + LinearInterpolator2D interp(d_data.data, h_shape.data, lo, hi); + + 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::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; + this->m_tuner->begin(); + gpu::perturbed_lennard_jones_evap_args_t args(d_force.data, + d_pos.data, + d_scale_factor.data, + box, + N, + n_ghost, + d_n_neigh.data, + d_nlist.data, + d_head_list.data, + interp, + scaled_t, + interface_height, + rcutsq, + d_params.data, + ntypes, + m_energy_shift, + m_tuner->getParam()[0]); + + gpu::compute_perturbed_lennard_jones_evap_forces(args); + if (this->m_exec_conf->isCUDAErrorCheckingEnabled()) + CHECK_CUDA_ERROR(); + 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 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"); + + return std::make_shared(sysdef, + nlist, + rcut, + scale_factor, + 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 new file mode 100644 index 00000000..14b2fd97 --- /dev/null +++ b/src/PerturbedLennardJonesEvapGPU.cu @@ -0,0 +1,171 @@ +// Copyright (c) 2018-2020, Michael P. Howard +// Copyright (c) 2021-2025, Auburn University +// Part of azplugins, released under the BSD 3-Clause License. + +#include "PairEvaluatorPerturbedLennardJones.h" +#include "PerturbedLennardJonesEvapGPU.cuh" + +namespace hoomd + { +namespace azplugins + { +namespace detail + { +namespace gpu + { +namespace kernel + { + +__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; + 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, + 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, + 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 + Scalar3 fi = make_scalar3(0, 0, 0); + Scalar pei = 0; + + const unsigned int n_neigh = d_n_neigh[idx]; + const size_t head = d_head_list[idx]; + + // loop over this particle's neighbors + for (unsigned int k = 0; k < n_neigh; ++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 unsigned int type_j = __scalar_as_int(postype_j.w); + 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); + + 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); + bool evaluated = eval.evalForceAndEnergy(force_divr, pair_eng, energy_shift); + + if (evaluated) + { + 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) + { + assert(args.block_size != 0); + + hipFuncAttributes attr; + hipFuncGetAttributes( + &attr, + reinterpret_cast(&kernel::compute_perturbed_lennard_jones_evap_forces)); + 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 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((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((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.rcutsq, + args.d_params, + args.ntypes, + 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 new file mode 100644 index 00000000..57f3563b --- /dev/null +++ b/src/PerturbedLennardJonesEvapGPU.cuh @@ -0,0 +1,80 @@ +// 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 "PairEvaluatorPerturbedLennardJones.h" +#include + +namespace hoomd + { +namespace azplugins + { +namespace detail + { +namespace gpu + { +//! Wraps arguments to kernel driver +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, + const LinearInterpolator2D& _interp, + const Scalar _scaled_t, + const Scalar _interface_height, + const Scalar _rcutsq, + 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), d_params(_d_params), ntypes(_ntypes), 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; + const LinearInterpolator2D interp; + const Scalar scaled_t; + const Scalar interface_height; + const Scalar rcutsq; + 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; + }; + +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 + +#endif // AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_GPU_CUH_ diff --git a/src/PerturbedLennardJonesEvapGPU.h b/src/PerturbedLennardJonesEvapGPU.h new file mode 100644 index 00000000..5b777b0c --- /dev/null +++ b/src/PerturbedLennardJonesEvapGPU.h @@ -0,0 +1,53 @@ +// 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_H_ +#define AZPLUGINS_PERTURBED_LENNARD_JONES_EVAP_GPU_H_ + +#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 + { +namespace detail + { +class PerturbedLennardJonesEvapGPU : public PerturbedLennardJonesEvap + { + public: + //! Constructor + + typedef PairParametersPerturbedLennardJones param_type; + + PerturbedLennardJonesEvapGPU(std::shared_ptr sysdef, + std::shared_ptr nlist, + const Scalar r_cut, + const Scalar scale_factor, + 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; + }; + +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/VariantInterpolated.h b/src/VariantInterpolated.h new file mode 100644 index 00000000..901c13ab --- /dev/null +++ b/src/VariantInterpolated.h @@ -0,0 +1,130 @@ +// 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_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); + const Scalar t = static_cast(timestep); + // LinearInterpolator1D + 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 + { + return m_t_lo; + } + + Scalar getTHi() const + { + 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() + { + 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 052d4c13..d892805a 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.3.0" diff --git a/src/export_VariantInterpolated.cc b/src/export_VariantInterpolated.cc new file mode 100644 index 00000000..d792ebba --- /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("t_lo", &VariantInterpolated::getTLo, &VariantInterpolated::setTLo) + .def_property("t_hi", &VariantInterpolated::getTHi, &VariantInterpolated::setTHi); + } + } // namespace detail + } // namespace azplugins + } // namespace hoomd diff --git a/src/module.cc b/src/module.cc index d8658225..bca26e9a 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&); @@ -102,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&); @@ -153,6 +157,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); @@ -174,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/pair.py b/src/pair.py index 9c4a4e68..51b29aca 100644 --- a/src/pair.py +++ b/src/pair.py @@ -4,10 +4,13 @@ """Pair potentials.""" +import numpy +import hoomd 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,125 @@ def __init__(self, nlist, default_r_cut=None, default_r_on=0, mode="none"): self._add_typeparam(params) +class PerturbedLennardJonesEvap(Force): + 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 + 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 + 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. + 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:: + + 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" + + def __init__( + self, + nlist, + rcut, + time_scale_factor, + energy_shift, + attraction_scale_factor_data, + domain, + variant, + ): + super().__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, + time_scale_factor=float, + energy_shift=bool, + ) + param_dict["rcut"] = float(rcut) + param_dict["time_scale_factor"] = float(time_scale_factor) + param_dict["energy_shift"] = bool(energy_shift) + 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( + self._attraction_scale_factor_data.shape, dtype=numpy.uint64 + ) + + def _attach_hook(self): + self._nlist._attach(self._simulation) + + 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, + self.rcut, + self.time_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/pytest/CMakeLists.txt b/src/pytest/CMakeLists.txt index 157b326b..dc8630c3 100644 --- a/src/pytest/CMakeLists.txt +++ b/src/pytest/CMakeLists.txt @@ -8,7 +8,9 @@ set(test_files test_pair.py test_pair_aniso.py test_pair_dpd.py + test_perturbedljevap.py test_wall.py + test_variant.py ) # Copy tests to the install directory. diff --git a/src/pytest/test_perturbedljevap.py b/src/pytest/test_perturbedljevap.py new file mode 100644 index 00000000..ca4a405d --- /dev/null +++ b/src/pytest/test_perturbedljevap.py @@ -0,0 +1,279 @@ +# Copyright (c) 2018-2020, Michael P. Howard +# 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 + +import pytest + + +@pytest.fixture +def two_particle_snapshot_factory(): + + 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 +def valid_args_const(): + return { + "nlist": hoomd.md.nlist.Cell(buffer=0.4), + "rcut": 3.0, + "time_scale_factor": 1.0, + "energy_shift": False, + "attraction_scale_factor_data": numpy.array([[0.6, 0.6], [0.6, 0.6]]), + "domain": [0.0, 100.0], + "variant": hoomd.azplugins.variant.VariantInterpolated( + [5.0, 4.0, 2.0, 1.0], 0, 300 + ), + } + + +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.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"] + + # 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"], + ) + + assert evap.rcut == 3.0 + assert evap.time_scale_factor == 1.0 + + +def test_domain_mismatch( + valid_args_const, two_particle_snapshot_factory, simulation_factory +): + snap = two_particle_snapshot_factory([[0, 0, 0.1], [0, 0, 0.1]]) + sim = simulation_factory(snap) + + bad_args = valid_args_const.copy() + bad_args["domain"] = [10.0] + + integrator = hoomd.md.Integrator(dt=0.005) + 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): + sim.run(0) + + +def test_data_size_mismatch( + valid_args_const, two_particle_snapshot_factory, simulation_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] + + sim = simulation_factory(snap) + + integrator = hoomd.md.Integrator(dt=0.005) + 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): + sim.run(0) + + +def test_variant_mismatch( + valid_args_const, two_particle_snapshot_factory, simulation_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 + ) # Invalid: Should be VariantInterpolated + + sim = simulation_factory(snap) + + integrator = hoomd.md.Integrator(dt=0.005) + 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): + sim.run(0) + + +"""Test energy and force calculation for single particle type +with constant attraction scalefactor. +""" + + +@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, + positions, + epsilon, + expected_forces, + expected_energies, +): + 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) + + integrator = hoomd.md.Integrator(dt=0.005) + integrator.forces = [evap] + + sim.operations.integrator = integrator + + sim.run(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) + + +@pytest.fixture( + params=[ + # (time_scale_factor, domain) + (1.0, [0.0, 100.0]), + (2.0, [0.0, 50.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, + "time_scale_factor": time_scale_factor, + "energy_shift": False, + "attraction_scale_factor_data": attraction_factor_table, + "domain": domain, + "variant": hoomd.azplugins.variant.VariantInterpolated( + [2, 0, -2, -4, -6, -8], + 0.0, + 100.0, # Run for 100 timesteps, which corresponds to t = 0.5 for dt = 0.005 + ), + } + + +"""Test energy and force calculation for single particle type +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, + 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(positions) + evap = hoomd.azplugins.pair.PerturbedLennardJonesEvap(**valid_args_vary) + evap.params[("A", "A")] = dict(epsilon=epsilon, sigma=1.0) + sim = simulation_factory(snap) + + integrator = hoomd.md.Integrator(dt=0.005) + integrator.forces = [evap] + sim.operations.integrator = integrator + sim.run(0) + + forces = evap.forces + energies = evap.energies + + if sim.device.communicator.rank == 0: + 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) + forces = evap.forces + energies = evap.energies + + if sim.device.communicator.rank == 0: + numpy.testing.assert_allclose(forces, expected_forces_t100) + numpy.testing.assert_allclose(energies, expected_energies_t100) diff --git a/src/pytest/test_variant.py b/src/pytest/test_variant.py new file mode 100644 index 00000000..b9772809 --- /dev/null +++ b/src/pytest/test_variant.py @@ -0,0 +1,64 @@ +# Copyright (c) 2018-2020, Michael P. Howard +# 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 + +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), + ([2, 0, -2, -4, -6, -8], 0, 100), +] + + +@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) + + +@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("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(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)) diff --git a/src/variant.py b/src/variant.py new file mode 100644 index 00000000..f85e2869 --- /dev/null +++ b/src/variant.py @@ -0,0 +1,27 @@ +# 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(_azplugins.VariantInterpolated, hoomd.variant.Variant): + """Piecewise-linear variant on a uniform grid of time.""" + + def __init__(self, data, t_lo, t_hi): + hoomd.variant.Variant.__init__(self) + + 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 + _azplugins.VariantInterpolated.__init__( + self, data, data.size, float(t_lo), float(t_hi) + )