diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 8c50d5ff49..ee3bf865aa 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -7,6 +7,11 @@ Change Log 6.0.0 (not released) ^^^^^^^^^^^^^^^^^^^^ +*Added* + +* Add thermodynamically consistent frictional contact forces: ``hoomd.md.pair.friction`` + (`#2116 `__). + *Fixed* * Use the provided alpha parameter in ``make_pppm_coulomb_forces`` diff --git a/hoomd/RNGIdentifiers.h b/hoomd/RNGIdentifiers.h index bdcab2fa07..730dbad997 100644 --- a/hoomd/RNGIdentifiers.h +++ b/hoomd/RNGIdentifiers.h @@ -38,7 +38,7 @@ struct RNGIdentifier static const uint8_t Unused3 = 12; static const uint8_t Unused4 = 13; static const uint8_t Unused5 = 14; - static const uint8_t Unused6 = 15; + static const uint8_t EvaluatorPairFrictionLJBase = 15; static const uint8_t Unused7 = 16; static const uint8_t Unused8 = 17; static const uint8_t Unused9 = 18; diff --git a/hoomd/md/CMakeLists.txt b/hoomd/md/CMakeLists.txt index c22524bb19..fdec5294bd 100644 --- a/hoomd/md/CMakeLists.txt +++ b/hoomd/md/CMakeLists.txt @@ -23,6 +23,7 @@ set(_md_sources module-md.cc FIREEnergyMinimizer.cc ForceComposite.cc ForceDistanceConstraint.cc + FrictionPairFrictionLJBase.cc HalfStepHook.cc HarmonicAngleForceCompute.cc HarmonicDihedralForceCompute.cc @@ -109,6 +110,8 @@ set(_md_headers ActiveForceComputeGPU.h EvaluatorPairALJ.h EvaluatorPairBuckingham.h EvaluatorPairDipole.h + EvaluatorPairFrictionLJBase.h + EvaluatorPairFrictionLJVariants.h EvaluatorPairYLZ.h EvaluatorPairDLVO.h EvaluatorPairDPDThermoLJ.h @@ -216,6 +219,9 @@ set(_md_headers ActiveForceComputeGPU.h PPPMForceComputeGPU.h PPPMForceComputeGPU.cuh PPPMForceCompute.h + FrictionPairGPU.cuh + FrictionPairGPU.h + FrictionPair.h TableAngleForceComputeGPU.h TableAngleForceCompute.h TableAngleForceGPU.cuh @@ -286,6 +292,7 @@ list(APPEND _md_sources ActiveForceComputeGPU.cc OPLSDihedralForceComputeGPU.cc PeriodicImproperForceComputeGPU.cc PPPMForceComputeGPU.cc + FrictionPairFrictionLJBaseGPU.cc TableAngleForceComputeGPU.cc TableDihedralForceComputeGPU.cc TriangleAreaConservationMeshForceComputeGPU.cc @@ -327,6 +334,7 @@ set(_md_cu_sources ActiveForceComputeGPU.cu OPLSDihedralForceGPU.cu PeriodicImproperForceGPU.cu PPPMForceComputeGPU.cu + FrictionPairFrictionLJBaseGPUKernel.cu TableAngleForceGPU.cu TableDihedralForceGPU.cu TriangleAreaConservationMeshForceComputeGPU.cu diff --git a/hoomd/md/EvaluatorPairFrictionLJBase.h b/hoomd/md/EvaluatorPairFrictionLJBase.h new file mode 100644 index 0000000000..166fce3e7d --- /dev/null +++ b/hoomd/md/EvaluatorPairFrictionLJBase.h @@ -0,0 +1,371 @@ +// Copyright (c) 2009-2025 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +#ifndef __PAIR_EVALUATOR_FRICTIONLJBASE_H__ +#define __PAIR_EVALUATOR_FRICTIONLJBASE_H__ + +#ifndef __HIPCC__ +#include +#endif + +#ifdef ENABLE_HIP +#include +#endif +#include "hoomd/RNGIdentifiers.h" +#include "hoomd/RandomNumbers.h" +#include "hoomd/VectorMath.h" + +// need to declare these class methods with __device__ qualifiers when building +// in nvcc. HOSTDEVICE is __host__ __device__ when included in nvcc and blank +// when included into the host compiler +#ifdef __HIPCC__ +#define HOSTDEVICE __host__ __device__ +#define DEVICE __device__ +#else +#define HOSTDEVICE +#define DEVICE +#endif + +#define TOLERANCE_DU 1e-6 + +namespace hoomd + { +namespace md + { +template class EvaluatorPairFrictionLJBase + { + public: + // Add the eval_factors() method declaration + HOSTDEVICE void eval_factors(Scalar& factor_f, Scalar& factor_r, Scalar w, Scalar du) + { + static_cast(this)->eval_factors(factor_f, factor_r, w, du); + } + + struct param_type + { + Scalar sigma_6; // Sigma^6 parameter of the LJ-Potential + Scalar epsilon_x_4; // 4*Epsilon parameter of the LJ-Potential + Scalar + gamma; // Gamma parameter of the frictional contacts (optional depends on friction type) + Scalar + kappa; // kappa parameter of the frictional contacts (optional depends on friction type) + Scalar pair_temp; // Temperature of the pairwise thermostat + +#ifdef ENABLE_HIP + //! Set CUDA memory hints + void set_memory_hint() const + { + // default implementation does nothing + } +#endif + + //! Load dynamic data members into shared memory and increase pointer + /*! \param ptr Pointer to load data to (will be incremented) + \param available_bytes Size of remaining shared memory + allocation + */ + DEVICE void load_shared(char*& ptr, unsigned int& available_bytes) { } + + HOSTDEVICE void allocate_shared(char*& ptr, unsigned int& available_bytes) const { } + + HOSTDEVICE param_type() : sigma_6(0), epsilon_x_4(0), gamma(0), kappa(0), pair_temp(0) { } + +#ifndef __HIPCC__ + + param_type(pybind11::dict v, bool managed) + { + auto sigma(v["sigma"].cast()); + auto epsilon(v["epsilon"].cast()); + Scalar gamma_f = v.contains("gamma_f") ? v["gamma_f"].cast() : Scalar(0.0); + Scalar kappa_f = v.contains("kappa_f") ? v["kappa_f"].cast() : Scalar(0.0); + auto kT(v["kT"].cast()); + + sigma_6 = sigma * sigma * sigma * sigma * sigma * sigma; + epsilon_x_4 = Scalar(4.0) * epsilon; + gamma = gamma_f; + kappa = kappa_f; + pair_temp = kT; + } + + pybind11::object toPython() + { + pybind11::dict v; + v["sigma"] = pow(sigma_6, 1. / 6.); + v["epsilon"] = epsilon_x_4 / 4.0; + if (gamma != Scalar(0.0)) + v["gamma_f"] = gamma; + if (kappa != Scalar(0.0)) + v["kappa_f"] = kappa; + v["kT"] = Scalar(1.0) * pair_temp; + return v; + } + +#endif + } +#if HOOMD_LONGREAL_SIZE == 32 + __attribute__((aligned(8))); +#else + __attribute__((aligned(16))); +#endif + + //! Constructs the pair potential evaluator + /*! \param _dr Displacement vector between particle centers of mass + \param _rcutsq Squared distance at which the potential goes to 0 + \param _angvel_i Angular veloctiy ot i^{th} particle + \param _angvel_j Angular veloctiy ot j^{th} particle + \param _dia_i Diameter ot i^{th} particle + \param _dia_j Diameter ot i^{th} particle + \param _epsilon_x_4 Epsilon of the LJ Potential + \param _sigma_6 Sigma of the LJ Potential + \param gamma Gamma parameter of the frictional contact + \param kappa Kappa parameter of the frictional contact + \param pair_temp Pair temperature of the frictional contact + \param _params Per type pair parameters of this potential + */ + HOSTDEVICE EvaluatorPairFrictionLJBase(Scalar3& _dr, + Scalar3& _angvel_i, + Scalar3& _angvel_j, + Scalar3& _dv, + Scalar _dia_i, + Scalar _dia_j, + Scalar _rcutsq, + const param_type& _params) + : dr(_dr), rcutsq(_rcutsq), angvel_i(_angvel_i), angvel_j(_angvel_j), dv(_dv), + dia_i(_dia_i), dia_j(_dia_j), + lj1(_params.epsilon_x_4 * _params.sigma_6 * _params.sigma_6), + lj2(_params.epsilon_x_4 * _params.sigma_6), gamma(_params.gamma), kappa(_params.kappa), + pair_temp(_params.pair_temp) + { + } + + //! Whether the pair potential needs particle tags. + HOSTDEVICE static bool needsTags() + { + return false; + } + + //! whether pair potential requires charges + HOSTDEVICE static bool needsCharge() + { + return true; + } + + //! whether pair potential requires nu_ito + HOSTDEVICE static bool needsNu() + { + return true; + } + + // Seed, Timestep, and the particle ids are necessary for the correlation of the pair noise + // (equation 26 and 27 of manuscript) + HOSTDEVICE void + set_seed_ij_timestep(uint16_t seed, unsigned int i, unsigned int j, uint64_t timestep) + { + m_seed = seed; + m_i = i; + m_j = j; + m_timestep = timestep; + } + + //! Set the timestep size + HOSTDEVICE void setDeltaT(Scalar dt) + { + m_deltaT = dt; + } + + //! Accept the optional tags + /*! \param tag_i Tag of particle i + \param tag_j Tag of particle j + */ + HOSTDEVICE void setTags(unsigned int tagi, unsigned int tagj) { } + + //! Accept the optional charge values + /*! \param qi Charge of particle i + \param qj Charge of particle j + */ + HOSTDEVICE void setCharge(Scalar qi, Scalar qj) { } + + //! Accept the optional nu value + /*! \param nu_ito nu value for the ito formalism + */ + HOSTDEVICE void setNu(Scalar nu_ito) + { + nu_ij = nu_ito; + } + //! Set the third law value + /*! \param third_law + */ + HOSTDEVICE void setThirdLaw(bool third_law) + { + m_third_law = third_law; + } + + //! Evaluate the force and energy + /*! \param force Output parameter to write the computed force. + \param pair_eng Output parameter to write the computed pair energy. + \param torque_i The torque exerted on the i^th particle. + \param torque_j The torque exerted on the j^th particle. + \return True if they are evaluated or false if they are not because + we are beyond the cutoff. + */ + HOSTDEVICE bool evaluate(Scalar3& force, Scalar& pair_eng, Scalar3& torque_i, Scalar3& torque_j) + { + vec3 rvec(dr); + vec3 w_i(angvel_i); + vec3 w_j(angvel_j); + vec3 v_ij(dv); + Scalar d_i(dia_i); + Scalar d_j(dia_j); + + Scalar rsq = dot(rvec, rvec); + + if (rsq > rcutsq) + return false; + + //! Define the force and torques which get applied to the particles + vec3 f; + vec3 t_j(0.0, 0.0, 0.0); + vec3 t_i(0.0, 0.0, 0.0); + + //! Calculation of the repulsive LJ-Force + Scalar rinv = fast::rsqrt(rsq); + Scalar r2inv = Scalar(1.0) / rsq; + Scalar r6inv = r2inv * r2inv * r2inv; + + //! Calculation of the repulsive LJ-Force + Scalar force_divr = r2inv * r6inv * (Scalar(12.0) * lj1 * r6inv - Scalar(6.0) * lj2); + f = force_divr * rvec; + + //! Calculation of the rotational friction Force + + //! e_ij: unit vector between center of masses + vec3 e_ij = Scalar(-1.0) * rvec * rinv; + + //! Project v_ij onto perpendicular to e_ij + vec3 P_e_v = v_ij - (dot(v_ij, e_ij) * e_ij); + + //! Calculate (w_i*R_i+w_j*R_j) + vec3 wiRiwjRj = Scalar(0.5) * (d_i * w_i + d_j * w_j); + + //! Calculate the relative tangential velocity u_ij at the contact point + vec3 u_ij = P_e_v - cross(wiRiwjRj, e_ij); + Scalar du = fast::sqrt(dot(u_ij, u_ij)); + + //! du should not be too small to avoid division by zero + if (TOLERANCE_DU < du) + { + unsigned int m_oi, m_oj; + int sign_xi = 1; + // initialize the RNG + if (m_i > m_j) + { + m_oi = m_j; + m_oj = m_i; + } + else + { + m_oi = m_i; + m_oj = m_j; + if (!m_third_law) + sign_xi + = -1; // If the neighborlist is only halfe the antisymmetric correlations + // have to be handled here! (Maybe there is a better way?) + } + + //! Init Random number Generator + hoomd::RandomGenerator rng( + hoomd::Seed(hoomd::RNGIdentifier::EvaluatorPairFrictionLJBase, m_timestep, m_seed), + hoomd::Counter(m_oi, m_oj)); + + //! Distant dependend factor + Scalar w = fast::sqrt(dot(f, f)); + + Scalar factor_f = Scalar(0.0); + Scalar factor_r = Scalar(0.0); + + //! How the factors are calculated depends on the friction type (currently implemented + //! in EvaluatorPairFrictionLJVariants.h is linear, constant and coulombNewton) + eval_factors(factor_f, factor_r, w, du); + + //! Calculation of the Rotational Friction force and torque + vec3 f_f = factor_f * (P_e_v + cross(e_ij, wiRiwjRj)); + vec3 exff = cross(e_ij, f_f); + + //! Noise for rotational friction + Scalar sigma_f = fast::sqrt(pair_temp / m_deltaT); + + Scalar xi_x = sign_xi * hoomd::NormalDistribution(sigma_f)(rng); + Scalar xi_y = sign_xi * hoomd::NormalDistribution(sigma_f)(rng); + Scalar xi_z = sign_xi * hoomd::NormalDistribution(sigma_f)(rng); + + Scalar N_x = hoomd::NormalDistribution(sigma_f)(rng); + Scalar N_y = hoomd::NormalDistribution(sigma_f)(rng); + Scalar N_z = hoomd::NormalDistribution(sigma_f)(rng); + + vec3 xi = vec3(xi_x, xi_y, xi_z); + vec3 N = vec3(N_x, N_y, N_z); + + vec3 f_r = factor_r * (xi - dot(xi, e_ij) * e_ij - cross(e_ij, N)); + vec3 t_r = factor_r * (cross(e_ij, xi) + N - dot(N, e_ij) * e_ij); + + //! Calculate the torque on the particles + t_i = (d_i / Scalar(2.0)) * (exff + t_r); + t_j = (d_j / Scalar(2.0)) * (exff + t_r); + + //! Add all forces + f = f + f_f + f_r; + } + + force = vec_to_scalar3(f); + torque_i = vec_to_scalar3(t_i); + torque_j = vec_to_scalar3(t_j); + // pair_eng = e; + + return true; + } + + DEVICE Scalar evalPressureLRCIntegral() + { + return 0; + } + + DEVICE Scalar evalEnergyLRCIntegral() + { + return 0; + } + +#ifndef __HIPCC__ + //! Get the name of the potential + /*! \returns The potential name. + */ + static std::string getName() + { + return "FrictionLJBase"; + } +#endif + + protected: + Scalar3 dr; //!< Stored vector pointing between particle centers of mass + Scalar rcutsq; //!< Stored rcutsq from the constructor + Scalar3 angvel_i, + angvel_j; //!< Stored angular momentum of ith and jth particle from the constructor + Scalar3 dv; //!< Stored velocity difference vij between the ith and jth particle + Scalar dia_i, dia_j; //!< Stored diameter of ith and jth particle from the constructor + Scalar nu_ij; //!< Factor which depends if Ito or the Stratonovich case. + bool m_third_law; //!< Boolean storing if only a half neighborlist is used + Scalar lj1, lj2; //!< lj1 and lj2 parameter extracted from the params passed to the constructor + Scalar gamma; //!< Optional gamma parameter from the constructor + Scalar kappa; //!< Optional kappa parameter from the constructor + Scalar pair_temp; //!< User set temperature for the DPD like PRNG + Scalar m_deltaT; //!< Timestep size stored from constructor + uint16_t m_seed; //!< User set seed for thermostat PRNG + unsigned int m_i; //!< Index of first particle. For use in PRNG + unsigned int m_j; //!< Index of second particle. For use in PRNG + uint64_t m_timestep; //!< timestep for use in PRNG + // const param_type ¶ms; //!< The pair potential parameters + }; + + } // end namespace md + } // end namespace hoomd + +#endif // __PAIR_EVALUATOR_FRICTIONLJBASE_H__ diff --git a/hoomd/md/EvaluatorPairFrictionLJVariants.h b/hoomd/md/EvaluatorPairFrictionLJVariants.h new file mode 100644 index 0000000000..9c0bd4bd0a --- /dev/null +++ b/hoomd/md/EvaluatorPairFrictionLJVariants.h @@ -0,0 +1,123 @@ +// Copyright (c) 2009-2025 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +#ifndef __PAIR_EVALUATOR_FRICTIONLJVARIANTS_H__ +#define __PAIR_EVALUATOR_FRICTIONLJVARIANTS_H__ + +/*! \file EvaluatorPairFrictionLJVariants.h + \brief Defines the different variants of frictional contact interactions +*/ + +// need to declare these class methods with __device__ qualifiers when building +// in nvcc. HOSTDEVICE is __host__ __device__ when included in nvcc and blank +// when included into the host compiler + +#ifdef __HIPCC__ +#define HOSTDEVICE __host__ __device__ +#define DEVICE __device__ +#else +#define HOSTDEVICE +#define DEVICE +#endif + +namespace hoomd + { +namespace md + { + +class EvaluatorPairFrictionLJLinear + : public EvaluatorPairFrictionLJBase + { + public: + using EvaluatorPairFrictionLJBase::EvaluatorPairFrictionLJBase; + + HOSTDEVICE void eval_factors(Scalar& factor_f, Scalar& factor_r, Scalar w, Scalar du) + { + Scalar f = gamma; + + factor_f = w * f; + factor_r = fast::sqrt(factor_f); + } + +#ifndef __HIPCC__ + static std::string getName() + { + return "FrictionLJLinear"; + } +#endif + }; + +class EvaluatorPairFrictionLJCoulomb + : public EvaluatorPairFrictionLJBase + { + public: + using EvaluatorPairFrictionLJBase::EvaluatorPairFrictionLJBase; + + HOSTDEVICE void eval_factors(Scalar& factor_f, Scalar& factor_r, Scalar w, Scalar du) + { + Scalar f = kappa / du; + Scalar D = Scalar(0.0); + if (pair_temp > 0.0) + D = kappa * fast::sqrt(M_PI / (Scalar(2.0) * pair_temp * nu_ij)) + * exp(du * du / (2 * pair_temp * nu_ij)) + * erfc(du / fast::sqrt(Scalar(2.0) * pair_temp * nu_ij)); + + factor_f = w * f; + factor_r = fast::sqrt(w * D); + } + +#ifndef __HIPCC__ + static std::string getName() + { + return "FrictionLJCoulomb"; + } +#endif + }; + +class EvaluatorPairFrictionLJCoulombNewton + : public EvaluatorPairFrictionLJBase + { + public: + using EvaluatorPairFrictionLJBase< + EvaluatorPairFrictionLJCoulombNewton>::EvaluatorPairFrictionLJBase; + + HOSTDEVICE void eval_factors(Scalar& factor_f, Scalar& factor_r, Scalar w, Scalar du) + { + Scalar f = Scalar(0.0); + Scalar D = Scalar(0.0); + + if (du < (w * kappa / gamma)) + { + if (pair_temp > 0.0) + D = w * kappa * fast::sqrt(M_PI / (Scalar(2.0) * pair_temp * nu_ij)) + * exp(du * du / (Scalar(2.0) * pair_temp * nu_ij)) + * erfc(w * kappa / (gamma * fast::sqrt(Scalar(2.0) * pair_temp * nu_ij))) + + gamma + * (Scalar(1.0) + - exp((du * du - (w * kappa / gamma) * (w * kappa / gamma)) + / (Scalar(2.0) * pair_temp * nu_ij))); + f = gamma; + } + else + { + if (pair_temp > 0.0) + D = w * kappa * fast::sqrt(M_PI / (Scalar(2.0) * pair_temp * nu_ij)) + * exp(du * du / (Scalar(2.0) * pair_temp * nu_ij)) + * erfc(du / fast::sqrt(Scalar(2.0) * pair_temp * nu_ij)); + f = w * kappa / du; + } + + factor_f = f; + factor_r = fast::sqrt(D); + } + +#ifndef __HIPCC__ + static std::string getName() + { + return "FrictionLJCoulombNewton"; + } +#endif + }; + } // namespace md + } // namespace hoomd +#endif // __PAIR_EVALUATOR_FRICTIONLJVARIANTS_H__ diff --git a/hoomd/md/FrictionPair.h b/hoomd/md/FrictionPair.h new file mode 100644 index 0000000000..18f39bb6d6 --- /dev/null +++ b/hoomd/md/FrictionPair.h @@ -0,0 +1,592 @@ +// Copyright (c) 2009-2025 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +#ifndef __FRICTION_POTENTIAL_PAIR_H__ +#define __FRICTION_POTENTIAL_PAIR_H__ + +#include +#include +#include +#include + +#ifdef ENABLE_HIP +#include +#endif + +#include "NeighborList.h" +#include "hoomd/ForceCompute.h" + +#include "hoomd/ManagedArray.h" +#include "hoomd/VectorMath.h" + +/*! \file FrictionPair.h + \brief Defines the template class for friction pair potentials + \details The heart of the code that computes friction pair potentials is in this file. + \note This header cannot be compiled by nvcc +*/ + +#ifdef __HIPCC__ +#error This header cannot be compiled by nvcc +#endif + +#include + +namespace hoomd + { +namespace md + { + +//! Template class for computing frictional contacts +/*! Overview: + FrictionPair is modeled after AnisoPotentialPair and computes frictional pair forces (and + torques) between all particles in the simulation. It employs the use of a neighbor list to limit + the number of computations done to only those particles with the cutoff radius of each other. The + computation of the force (and torque) is not performed directly by this class, but by an + friction_evaluator class (e.g. EvaluatorPairFrictionLJBase) which is passed in as a template + parameter so the computations are performed as efficiently as possible. + + Implementation details + FrictionPair is accessing the velocity, angular momentum, diameter, and moment of inertia of + all particles. It is used to calculate the relative velocity between the particles dv and the + angular velocity angvel_i/j of each particle. This information is passed to the + friction_evaluator class for the calculation of frictional pair force (and torque). +*/ + +template class FrictionPair : public ForceCompute + { + public: + //! Param type from friction_evaluator + typedef typename friction_evaluator::param_type param_type; + + //! Construct the pair potential + FrictionPair(std::shared_ptr sysdef, std::shared_ptr nlist); + //! Destructor + virtual ~FrictionPair(); + + //! Set the pair parameters for a single type pair + virtual void setParams(unsigned int typ1, unsigned int typ2, const param_type& param); + + virtual void setParamsPython(pybind11::tuple typ, pybind11::object params); + + /// Get params for a single type pair using a tuple of strings + virtual pybind11::object getParamsPython(pybind11::tuple typ); + + //! Set the rcut for a single type pair + virtual void setRcut(unsigned int typ1, unsigned int typ2, Scalar rcut); + + /// Get the r_cut for a single type pair + Scalar getRCut(pybind11::tuple types); + + /// Set the rcut for a single type pair using a tuple of strings + virtual void setRCutPython(pybind11::tuple types, Scalar r_cut); + + /// Validate that types are within Ntypes + virtual void validateTypes(unsigned int typ1, unsigned int typ2, std::string action); + + void setShiftModePython(std::string mode) + { + if (mode != "none") + { + throw std::runtime_error("Invalid energy shift mode."); + } + } + + /// Get the mode used for the energy shifting + std::string getShiftMode() + { + return "none"; + } + + virtual void notifyDetach() + { + if (m_attached) + { + m_nlist->removeRCutMatrix(m_r_cut_nlist); + } + m_attached = false; + } + +#ifdef ENABLE_MPI + //! Get ghost particle fields requested by this pair potential + virtual CommFlags getRequestedCommFlags(uint64_t timestep); +#endif + + /// Start autotuning kernel launch parameters + virtual void startAutotuning() + { + ForceCompute::startAutotuning(); + + // Start autotuning the neighbor list. + m_nlist->startAutotuning(); + } + + /// Check if autotuning is complete. + virtual bool isAutotuningComplete() + { + bool result = ForceCompute::isAutotuningComplete(); + result = result && m_nlist->isAutotuningComplete(); + return result; + } + + protected: + std::shared_ptr m_nlist; //!< The neighborlist to use for the computation + Index2D m_typpair_idx; //!< Helper class for indexing per type pair arrays + GPUArray m_rcutsq; //!< Cutoff radius squared per type pair + std::vector> + m_params; //!< Pair parameters per type pair + + /// Track whether we have attached to the Simulation object + bool m_attached = true; + + /// r_cut (not squared) given to the neighbor list + std::shared_ptr> m_r_cut_nlist; + + //! Actually compute the forces + virtual void computeForces(uint64_t timestep); + }; + +/*! \param sysdef System to compute forces on + \param nlist Neighborlist to use for computing the forces +*/ +template +FrictionPair::FrictionPair(std::shared_ptr sysdef, + std::shared_ptr nlist) + : ForceCompute(sysdef), m_nlist(nlist), m_typpair_idx(m_pdata->getNTypes()) + { + m_exec_conf->msg->notice(5) << "Constructing FrictionPair<" << friction_evaluator::getName() + << ">" << std::endl; + assert(m_pdata); + assert(m_nlist); + + GPUArray rcutsq(m_typpair_idx.getNumElements(), m_exec_conf); + m_rcutsq.swap(rcutsq); + GPUArray ronsq(m_typpair_idx.getNumElements(), m_exec_conf); + std::vector> params( + static_cast(m_typpair_idx.getNumElements()), + param_type(), + hoomd::detail::managed_allocator(m_exec_conf->isCUDAEnabled())); + m_params.swap(params); + m_r_cut_nlist = std::make_shared>(m_typpair_idx.getNumElements(), m_exec_conf); + nlist->addRCutMatrix(m_r_cut_nlist); + } + +template FrictionPair::~FrictionPair() + { + m_exec_conf->msg->notice(5) << "Destroying FrictionPair<" << friction_evaluator::getName() + << ">" << std::endl; + + if (m_attached) + { + m_nlist->removeRCutMatrix(m_r_cut_nlist); + } + } + +/*! \param typ1 First type index in the pair + \param typ2 Second type index in the pair + \param param Parameter to set + \note When setting the value for (\a typ1, \a typ2), the parameter for (\a typ2, \a typ1) is + automatically set. +*/ +template +void FrictionPair::setParams(unsigned int typ1, + unsigned int typ2, + const param_type& param) + { + validateTypes(typ1, typ2, "setting params"); + m_params[m_typpair_idx(typ1, typ2)] = param; + m_params[m_typpair_idx(typ2, typ1)] = param; + } + +template +void FrictionPair::setParamsPython(pybind11::tuple typ, pybind11::object params) + { + auto typ1 = m_pdata->getTypeByName(typ[0].cast()); + auto typ2 = m_pdata->getTypeByName(typ[1].cast()); + setParams(typ1, typ2, param_type(params, m_exec_conf->isCUDAEnabled())); + } + +template +pybind11::object FrictionPair::getParamsPython(pybind11::tuple typ) + { + auto typ1 = m_pdata->getTypeByName(typ[0].cast()); + auto typ2 = m_pdata->getTypeByName(typ[1].cast()); + validateTypes(typ1, typ2, "getting params"); + + return m_params[m_typpair_idx(typ1, typ2)].toPython(); + } + +template +void FrictionPair::validateTypes(unsigned int typ1, + unsigned int typ2, + std::string action) + { + auto n_types = this->m_pdata->getNTypes(); + if (typ1 >= n_types || typ2 >= n_types) + { + throw std::runtime_error("Error in" + action + " for pair potential. Invalid type"); + } + } + +/*! \param typ1 First type index in the pair + \param typ2 Second type index in the pair + \param rcut Cutoff radius to set + \note When setting the value for (\a typ1, \a typ2), the parameter for (\a typ2, \a typ1) is + automatically set. +*/ +template +void FrictionPair::setRcut(unsigned int typ1, unsigned int typ2, Scalar rcut) + { + validateTypes(typ1, typ2, "setting r_cut"); + { + // store r_cut**2 for use internally + ArrayHandle h_rcutsq(m_rcutsq, access_location::host, access_mode::readwrite); + h_rcutsq.data[m_typpair_idx(typ1, typ2)] = rcut * rcut; + h_rcutsq.data[m_typpair_idx(typ2, typ1)] = rcut * rcut; + + // store r_cut unmodified for so the neighbor list knows what particles to include + ArrayHandle h_r_cut_nlist(*m_r_cut_nlist, + access_location::host, + access_mode::readwrite); + h_r_cut_nlist.data[m_typpair_idx(typ1, typ2)] = rcut; + h_r_cut_nlist.data[m_typpair_idx(typ2, typ1)] = rcut; + } + + // notify the neighbor list that we have changed r_cut values + m_nlist->notifyRCutMatrixChange(); + } + +template +void FrictionPair::setRCutPython(pybind11::tuple types, Scalar r_cut) + { + auto typ1 = m_pdata->getTypeByName(types[0].cast()); + auto typ2 = m_pdata->getTypeByName(types[1].cast()); + setRcut(typ1, typ2, r_cut); + } + +template +Scalar FrictionPair::getRCut(pybind11::tuple types) + { + auto typ1 = m_pdata->getTypeByName(types[0].cast()); + auto typ2 = m_pdata->getTypeByName(types[1].cast()); + validateTypes(typ1, typ2, "getting r_cut."); + ArrayHandle h_rcutsq(m_rcutsq, access_location::host, access_mode::read); + return sqrt(h_rcutsq.data[m_typpair_idx(typ1, typ2)]); + } + +/*! \post The pair forces are computed for the given timestep. The neighborlist's compute method is + called to ensure that it is up to date before proceeding. + + \param timestep specifies the current time step of the simulation +*/ +template +void FrictionPair::computeForces(uint64_t timestep) + { + // start by updating the neighborlist + m_nlist->compute(timestep); + + // depending on the neighborlist settings, we can take advantage of newton's third law + // to reduce computations at the cost of memory access complexity: set that flag now + bool third_law = m_nlist->getStorageMode() == NeighborList::half; + + // access the neighbor list, particle data, and system box + 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); + + ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::read); + ArrayHandle h_vel(m_pdata->getVelocities(), access_location::host, access_mode::read); + ArrayHandle h_charge(m_pdata->getCharges(), access_location::host, access_mode::read); + ArrayHandle h_orientation(m_pdata->getOrientationArray(), + access_location::host, + access_mode::read); + ArrayHandle h_angmom(m_pdata->getAngularMomentumArray(), + access_location::host, + access_mode::read); + ArrayHandle h_diameter(m_pdata->getDiameters(), + access_location::host, + access_mode::read); + ArrayHandle h_moment_inertia(m_pdata->getMomentsOfInertiaArray(), + access_location::host, + access_mode::read); + + ArrayHandle h_tag(m_pdata->getTags(), access_location::host, access_mode::read); + + // force arrays + ArrayHandle h_force(m_force, access_location::host, access_mode::overwrite); + ArrayHandle h_torque(m_torque, access_location::host, access_mode::overwrite); + ArrayHandle h_virial(m_virial, access_location::host, access_mode::overwrite); + + const BoxDim box = m_pdata->getBox(); + ArrayHandle h_rcutsq(m_rcutsq, access_location::host, access_mode::read); + { + // need to start from a zero force, energy and virial + m_force.zeroFill(); + m_torque.zeroFill(); + m_virial.zeroFill(); + + PDataFlags flags = this->m_pdata->getFlags(); + bool compute_virial = flags[pdata_flag::pressure_tensor]; + + uint16_t seed = this->m_sysdef->getSeed(); + + // for each particle + for (int i = 0; i < (int)m_pdata->getN(); i++) + { + // access the particle's position and type + Scalar3 pi = make_scalar3(h_pos.data[i].x, h_pos.data[i].y, h_pos.data[i].z); + + quat quat_i(h_orientation.data[i]); + + // get the velocitie of i-th particle and transform it from body to principle frame + Scalar3 vel_i = make_scalar3(h_vel.data[i].x, h_vel.data[i].y, h_vel.data[i].z); + + unsigned int typei = __scalar_as_int(h_pos.data[i].w); + Scalar dia_i = h_diameter.data[i]; + + // calculate angmom_i (angular momentum of th particle in the body frame) + quat angmom_i(h_angmom.data[i]); + quat qxP_i = conj(quat_i) * angmom_i; + vec3 bf_vel_i = Scalar(0.5) + * vec3(qxP_i.v.x / h_moment_inertia.data[i].x, + qxP_i.v.y / h_moment_inertia.data[i].y, + qxP_i.v.z / h_moment_inertia.data[i].z); + + // Rotate angular velocity into global frame + vec3 gf_angvel_i = rotate(quat_i, bf_vel_i); + Scalar3 angvel_i = make_scalar3(gf_angvel_i.x, gf_angvel_i.y, gf_angvel_i.z); + + // sanity check + assert(typei < m_pdata->getNTypes()); + + // access charge (if needed) + Scalar qi = Scalar(0.0); + Scalar mass_i = Scalar(0.0); + if (friction_evaluator::needsCharge()) + qi = h_charge.data[i]; + + if (friction_evaluator::needsNu()) + mass_i = h_vel.data[i].w; + + // initialize current particle force, torque, potential energy, and virial to 0 + Scalar fxi = Scalar(0.0); + Scalar fyi = Scalar(0.0); + Scalar fzi = Scalar(0.0); + Scalar txi = Scalar(0.0); + Scalar tyi = Scalar(0.0); + Scalar tzi = Scalar(0.0); + Scalar pei = Scalar(0.0); + Scalar virialxxi = 0.0; + Scalar virialxyi = 0.0; + Scalar virialxzi = 0.0; + Scalar virialyyi = 0.0; + Scalar virialyzi = 0.0; + Scalar virialzzi = 0.0; + + // loop over all of the neighbors of this particle + const size_t myHead = h_head_list.data[i]; + const unsigned int size = (unsigned int)h_n_neigh.data[i]; + for (unsigned int k = 0; k < size; k++) + { + // access the index of this neighbor + unsigned int j = h_nlist.data[myHead + k]; + assert(j < m_pdata->getN() + m_pdata->getNGhosts()); + + // calculate dr_ji + Scalar3 pj = make_scalar3(h_pos.data[j].x, h_pos.data[j].y, h_pos.data[j].z); + Scalar3 dx = pi - pj; + quat quat_j(h_orientation.data[j]); + + // get the velocitie of jth particle and transform it from body to principle frame + Scalar3 vel_j = make_scalar3(h_vel.data[j].x, h_vel.data[j].y, h_vel.data[j].z); + + // calculate dv_ij + Scalar3 dv = vel_j - vel_i; + + // calculate angmom_j (angular momentum of jth particle in the global frame) + quat angmom_j(h_angmom.data[j]); + quat qxP_j = conj(quat_j) * angmom_j; + vec3 bf_vel_j = Scalar(0.5) + * vec3(qxP_j.v.x / h_moment_inertia.data[j].x, + qxP_j.v.y / h_moment_inertia.data[j].y, + qxP_j.v.z / h_moment_inertia.data[j].z); + + // Rotate angular velocity into global frame + vec3 gf_angvel_j = rotate(quat_j, bf_vel_j); + Scalar3 angvel_j = make_scalar3(gf_angvel_j.x, gf_angvel_j.y, gf_angvel_j.z); + + // get the diameter of jth particle + Scalar dia_j = h_diameter.data[j]; + + // access the type of the neighbor particle + unsigned int typej = __scalar_as_int(h_pos.data[j].w); + assert(typej < m_pdata->getNTypes()); + + // access charge (if needed) + Scalar qj = Scalar(0.0); + Scalar mass_j = Scalar(0.0); + if (friction_evaluator::needsCharge()) + qj = h_charge.data[j]; + if (friction_evaluator::needsNu()) + mass_j = h_vel.data[j].w; + + // apply periodic boundary conditions + dx = box.minImage(dx); + + // get parameters for this type pair + unsigned int typpair_idx = m_typpair_idx(typei, typej); + const param_type& param = m_params[typpair_idx]; + Scalar rcutsq = h_rcutsq.data[typpair_idx]; + + // compute the force and potential energy + Scalar3 force = make_scalar3(0.0, 0.0, 0.0); + Scalar3 torque_i = make_scalar3(0.0, 0.0, 0.0); + Scalar3 torque_j = make_scalar3(0.0, 0.0, 0.0); + + Scalar pair_eng = Scalar(0.0); + + friction_evaluator eval(dx, angvel_i, angvel_j, dv, dia_i, dia_j, rcutsq, param); + + // set seed using global tags + unsigned int tagi = h_tag.data[i]; + unsigned int tagj = h_tag.data[j]; + + eval.set_seed_ij_timestep(seed, tagi, tagj, timestep); + eval.setDeltaT(this->m_deltaT); + eval.setThirdLaw(third_law); + + if (friction_evaluator::needsCharge()) + eval.setCharge(qi, qj); + if (friction_evaluator::needsTags()) + eval.setTags(h_tag.data[i], h_tag.data[j]); + if (friction_evaluator::needsNu()) + { + // Calculate nu for the Ito formalism + Scalar nu_ito + = ((Scalar(1.0) / mass_i) + (Scalar(1.0) / mass_j)) + + (((dia_j * dia_j / Scalar(4.0)) / h_moment_inertia.data[j].z) + + ((dia_i * dia_i / Scalar(4.0)) / h_moment_inertia.data[i].z)); + eval.setNu(nu_ito); + } + + bool evaluated = eval.evaluate(force, pair_eng, torque_i, torque_j); + + if (evaluated) + { + Scalar3 force2 = Scalar(0.5) * force; + + // add the force, potential energy and virial to the particle i + fxi += force.x; + fyi += force.y; + fzi += force.z; + txi += torque_i.x; + tyi += torque_i.y; + tzi += torque_i.z; + pei += pair_eng * Scalar(0.5); + + if (compute_virial) + { + virialxxi += dx.x * force2.x; + virialxyi += dx.y * force2.x; + virialxzi += dx.z * force2.x; + virialyyi += dx.y * force2.y; + virialyzi += dx.z * force2.y; + virialzzi += dx.z * force2.z; + } + + // add the force to particle j if we are using the third law + if (third_law && j < m_pdata->getN()) + { + h_force.data[j].x -= force.x; + h_force.data[j].y -= force.y; + h_force.data[j].z -= force.z; + h_torque.data[j].x += torque_j.x; + h_torque.data[j].y += torque_j.y; + h_torque.data[j].z += torque_j.z; + h_force.data[j].w += pair_eng * Scalar(0.5); + if (compute_virial) + { + h_virial.data[0 * m_virial_pitch + j] += dx.x * force2.x; + h_virial.data[1 * m_virial_pitch + j] += dx.y * force2.x; + h_virial.data[2 * m_virial_pitch + j] += dx.z * force2.x; + h_virial.data[3 * m_virial_pitch + j] += dx.y * force2.y; + h_virial.data[4 * m_virial_pitch + j] += dx.z * force2.y; + h_virial.data[5 * m_virial_pitch + j] += dx.z * force2.z; + } + } + } + } + + // finally, increment the force, potential energy and virial for particle i + h_force.data[i].x += fxi; + h_force.data[i].y += fyi; + h_force.data[i].z += fzi; + h_torque.data[i].x += txi; + h_torque.data[i].y += tyi; + h_torque.data[i].z += tzi; + h_force.data[i].w += pei; + if (compute_virial) + { + h_virial.data[0 * m_virial_pitch + i] += virialxxi; + h_virial.data[1 * m_virial_pitch + i] += virialxyi; + h_virial.data[2 * m_virial_pitch + i] += virialxzi; + h_virial.data[3 * m_virial_pitch + i] += virialyyi; + h_virial.data[4 * m_virial_pitch + i] += virialyzi; + h_virial.data[5 * m_virial_pitch + i] += virialzzi; + } + } + } + } + +#ifdef ENABLE_MPI +/*! \param timestep Current time step + */ +template +CommFlags FrictionPair::getRequestedCommFlags(uint64_t timestep) + { + CommFlags flags = CommFlags(0); + + // we need orientations for frictiontropic ptls + flags[comm_flag::orientation] = 1; + + if (friction_evaluator::needsCharge()) + flags[comm_flag::charge] = 1; + + // with rigid bodies, include net torque + flags[comm_flag::net_torque] = 1; + + flags |= ForceCompute::getRequestedCommFlags(timestep); + + return flags; + } +#endif + +namespace detail + { +//! Export this pair potential to python +/*! \param name Name of the class in the exported python module + \tparam T Evaluator type to export. +*/ +template void export_FrictionPair(pybind11::module& m, const std::string& name) + { + pybind11::class_, ForceCompute, std::shared_ptr>> + frictionpotentialpair(m, name.c_str()); + frictionpotentialpair + .def(pybind11::init, std::shared_ptr>()) + .def("setParams", &FrictionPair::setParamsPython) + .def("getParams", &FrictionPair::getParamsPython) + .def("setRCut", &FrictionPair::setRCutPython) + .def("getRCut", &FrictionPair::getRCut) + .def_property("mode", &FrictionPair::getShiftMode, &FrictionPair::setShiftModePython); + } + + } // end namespace detail + } // end namespace md + } // end namespace hoomd + +#endif // __FRICTION_POTENTIAL_PAIR_H__ diff --git a/hoomd/md/FrictionPairFrictionLJBase.cc b/hoomd/md/FrictionPairFrictionLJBase.cc new file mode 100644 index 0000000000..ffbe89b122 --- /dev/null +++ b/hoomd/md/FrictionPairFrictionLJBase.cc @@ -0,0 +1,34 @@ +// Copyright (c) 2009-2025 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +#include "EvaluatorPairFrictionLJBase.h" +#include "EvaluatorPairFrictionLJVariants.h" +#include "FrictionPair.h" + +namespace hoomd + { +namespace md + { +namespace detail + { + +void export_FrictionPairFrictionLJLinear(pybind11::module& m) + { + export_FrictionPair(m, "FrictionPairFrictionLJLinear"); + } + +void export_FrictionPairFrictionLJCoulomb(pybind11::module& m) + { + export_FrictionPair(m, "FrictionPairFrictionLJCoulomb"); + } + +void export_FrictionPairFrictionLJCoulombNewton(pybind11::module& m) + { + export_FrictionPair( + m, + "FrictionPairFrictionLJCoulombNewton"); + } + + } // end namespace detail + } // end namespace md + } // end namespace hoomd diff --git a/hoomd/md/FrictionPairFrictionLJBaseGPU.cc b/hoomd/md/FrictionPairFrictionLJBaseGPU.cc new file mode 100644 index 0000000000..6ae705aaa8 --- /dev/null +++ b/hoomd/md/FrictionPairFrictionLJBaseGPU.cc @@ -0,0 +1,33 @@ +// Copyright (c) 2009-2025 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +#include "EvaluatorPairFrictionLJBase.h" +#include "EvaluatorPairFrictionLJVariants.h" +#include "FrictionPairGPU.h" + +namespace hoomd + { +namespace md + { +namespace detail + { +void export_FrictionPairFrictionLJLinearGPU(pybind11::module& m) + { + export_FrictionPairGPU(m, "FrictionPairFrictionLJLinearGPU"); + } + +void export_FrictionPairFrictionLJCoulombGPU(pybind11::module& m) + { + export_FrictionPairGPU(m, "FrictionPairFrictionLJCoulombGPU"); + } + +void export_FrictionPairFrictionLJCoulombNewtonGPU(pybind11::module& m) + { + export_FrictionPairGPU( + m, + "FrictionPairFrictionLJCoulombNewtonGPU"); + } + + } // end namespace detail + } // end namespace md + } // end namespace hoomd diff --git a/hoomd/md/FrictionPairFrictionLJBaseGPUKernel.cu b/hoomd/md/FrictionPairFrictionLJBaseGPUKernel.cu new file mode 100644 index 0000000000..3a2f302cc7 --- /dev/null +++ b/hoomd/md/FrictionPairFrictionLJBaseGPUKernel.cu @@ -0,0 +1,30 @@ +// Copyright (c) 2009-2025 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +#include "EvaluatorPairFrictionLJBase.h" +#include "EvaluatorPairFrictionLJVariants.h" +#include "FrictionPairGPU.cuh" + +namespace hoomd + { +namespace md + { +namespace kernel + { +template hipError_t __attribute__((visibility("default"))) +gpu_compute_pair_friction_forces( + const a_pair_args_t& pair_args, + const EvaluatorPairFrictionLJLinear::param_type* d_param); + +template hipError_t __attribute__((visibility("default"))) +gpu_compute_pair_friction_forces( + const a_pair_args_t& pair_args, + const EvaluatorPairFrictionLJCoulomb::param_type* d_param); + +template hipError_t __attribute__((visibility("default"))) +gpu_compute_pair_friction_forces( + const a_pair_args_t& pair_args, + const EvaluatorPairFrictionLJCoulombNewton::param_type* d_param); + } // end namespace kernel + } // end namespace md + } // end namespace hoomd diff --git a/hoomd/md/FrictionPairGPU.cuh b/hoomd/md/FrictionPairGPU.cuh new file mode 100644 index 0000000000..3932542ae0 --- /dev/null +++ b/hoomd/md/FrictionPairGPU.cuh @@ -0,0 +1,624 @@ +// Copyright (c) 2009-2025 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +#include "hip/hip_runtime.h" +#include "hoomd/HOOMDMath.h" +#include "hoomd/Index1D.h" +#include "hoomd/ParticleData.cuh" +#include "hoomd/TextureTools.h" +#include "hoomd/VectorMath.h" + +#ifdef __HIPCC__ +#include "hoomd/WarpTools.cuh" +#endif // __HIPCC__ + +/*! \file FrictionPairGPU.cuh + \brief Defines templated GPU kernel code for calculating the friction ptl pair forces and + torques +*/ + +#ifndef __FRICTION_PAIR_GPU_CUH__ +#define __FRICTION_PAIR_GPU_CUH__ + +//! Maximum number of threads (width of a warp) +// currently this is hardcoded, we should set it to the max of platforms +#if defined(__HIP_PLATFORM_NVCC__) +const int gpu_friction_pair_force_max_tpp = 32; +#elif defined(__HIP_PLATFORM_HCC__) +const int gpu_friction_pair_force_max_tpp = 64; +#endif + +namespace hoomd + { +namespace md + { +namespace kernel + { +//! Wraps arguments to gpu_cgpf +struct a_pair_args_t + { + //! Construct a pair_args_t + a_pair_args_t(Scalar4* _d_force, + Scalar4* _d_torque, + Scalar* _d_virial, + size_t _virial_pitch, + const unsigned int _N, + const unsigned int _n_max, + const Scalar4* _d_pos, + const Scalar4* _d_vel, + const Scalar* _d_charge, + const Scalar4* _d_orientation, + const Scalar4* _d_angmom, + const Scalar* _d_diameter, + const Scalar3* _d_moment_inertia, + const unsigned int* _d_tag, + const BoxDim& _box, + const bool _third_law, + const unsigned int _dim, + uint16_t _seed, + uint64_t _timestep, + const Scalar _deltaT, + const unsigned int* _d_n_neigh, + const unsigned int* _d_nlist, + const size_t* _d_head_list, + const Scalar* _d_rcutsq, + const unsigned int _ntypes, + const unsigned int _block_size, + const unsigned int _compute_virial, + const unsigned int _threads_per_particle, + const hipDeviceProp_t& _devprop) + : d_force(_d_force), d_torque(_d_torque), d_virial(_d_virial), virial_pitch(_virial_pitch), + N(_N), n_max(_n_max), d_pos(_d_pos), d_vel(_d_vel), d_charge(_d_charge), + d_orientation(_d_orientation), d_angmom(_d_angmom), d_diameter(_d_diameter), + d_moment_inertia(_d_moment_inertia), d_tag(_d_tag), box(_box), third_law(_third_law), + dim(_dim), seed(_seed), timestep(_timestep), deltaT(_deltaT), d_n_neigh(_d_n_neigh), + d_nlist(_d_nlist), d_head_list(_d_head_list), d_rcutsq(_d_rcutsq), ntypes(_ntypes), + block_size(_block_size), compute_virial(_compute_virial), + threads_per_particle(_threads_per_particle), devprop(_devprop) { }; + + Scalar4* d_force; //!< Force to write out + Scalar4* d_torque; //!< Torque to write out + Scalar* d_virial; //!< Virial to write out + const size_t virial_pitch; //!< The pitch of the 2D array of virial matrix elements + const unsigned int N; //!< number of particles + const unsigned int n_max; //!< maximum size of particle data arrays + const Scalar4* d_pos; //!< particle positions + const Scalar4* d_vel; //!< particle velocity + const Scalar* d_charge; //!< particle charges + const Scalar4* d_orientation; //!< particle orientation to compute forces over + const Scalar4* d_angmom; //!< particle angular momentum + const Scalar* d_diameter; //!< particle diameter + const Scalar3* d_moment_inertia; //!< particle moment of inertia + const unsigned int* d_tag; //!< particle tags to compute forces over + const BoxDim box; //!< Simulation box in GPU format + const bool third_law; //!< Boolean storing if only a half neighborlist is used + const unsigned int dim; //!< Dimension of the simulation + uint16_t seed; //!< Seed of the simulation + uint64_t timestep; //!< Current timestep + const Scalar deltaT; //!< Timestep dt size of the simulation + const unsigned int* + d_n_neigh; //!< Device array listing the number of neighbors on each particle + const unsigned int* d_nlist; //!< Device array listing the neighbors of each particle + const size_t* d_head_list; //!< Device array listing beginning of each particle's neighbors + const Scalar* d_rcutsq; //!< Device array listing r_cut squared per particle type pair + const unsigned int ntypes; //!< Number of particle types in the simulation + const unsigned int block_size; //!< Block size to execute + const unsigned int compute_virial; //!< Flag to indicate if virials should be computed + const unsigned int threads_per_particle; //!< Number of threads to launch per particle + const hipDeviceProp_t& devprop; //!< CUDA device properties + }; + +#ifdef __HIPCC__ + +//! Kernel for calculating pair forces +/*! This kernel is called to calculate the pair forces on all N particles. Actual evaluation of the + potentials and forces for each pair is handled via the template class \a evaluator. + + \param d_force Device memory to write computed forces + \param d_torque Device memory to write computed torques + \param d_virial Device memory to write computed virials + \param virial_pitch pitch of 2D virial array + \param N number of particles in system + \param d_pos particle positions + \param d_vel particle velocity + \param d_charge particle charges + \param d_orientation Quaternion data on the GPU to calculate forces on + \param d_angmom particle angular momentum + \param d_diameter particle diameter + \param d_moment_inertia particle moment of inertia + \param d_tag Tag data on the GPU to calculate forces on + \param box Box dimensions used to implement periodic boundary conditions + \param third_law Boolean storing if only a half neighborlist is used + \param dim Dimension of the simulation + \param d_seed Seed of the simulation + \param d_timestep Current timestep + \param d_deltaT Timestep size of the simulation + \param d_n_neigh Device memory array listing the number of neighbors for each particle + \param d_nlist Device memory array containing the neighbor list contents + \param d_head_list Device memory array listing beginning of each particle's neighbors + \param d_params Parameters for the potential, stored per type pair + \param d_rcutsq rcut squared, stored per type pair + \param ntypes Number of types in the simulation + \param tpp Number of threads per particle + + \a d_params and \a d_rcutsq must be indexed with an Index2DUpperTriangular(typei, typej) to + access the unique value for that type pair. These values are all cached into shared memory for + quick access, so a dynamic amount of shared memory must be allocated for this kernel launch. The + amount is (2*sizeof(Scalar) + sizeof(typename evaluator::param_type)) * + typpair_idx.getNumElements() + + Certain options are controlled via template parameters to avoid the performance hit when they + are not enabled. \tparam evaluator EvaluatorPair class to evaluate V(r) and -delta V(r)/r \tparam + No energy shifting is done. \tparam compute_virial When + non-zero, the virial tensor is computed. When zero, the virial tensor is not computed. + + Implementation details + Each block will calculate the forces on a block of particles. + Each thread will calculate the total force on one particle. + The neighborlist is arranged in columns so that reads are fully coalesced when doing this. +*/ +template +__global__ void +gpu_compute_pair_friction_forces_kernel(Scalar4* d_force, + Scalar4* d_torque, + Scalar* d_virial, + const size_t virial_pitch, + const unsigned int N, + const Scalar4* d_pos, + const Scalar4* d_vel, + const Scalar* d_charge, + const Scalar4* d_orientation, + const Scalar4* d_angmom, + const Scalar* d_diameter, + const Scalar3* d_moment_inertia, + const unsigned int* d_tag, + const BoxDim box, + const bool third_law, + const unsigned int dim, + const uint16_t d_seed, + const uint64_t d_timestep, + const Scalar d_deltaT, + const unsigned int* d_n_neigh, + const unsigned int* d_nlist, + const size_t* d_head_list, + const typename evaluator::param_type* d_params, + const Scalar* d_rcutsq, + const unsigned int ntypes, + unsigned int max_extra_bytes) + { + Index2D typpair_idx(ntypes); + const unsigned int num_typ_parameters = typpair_idx.getNumElements(); + + // shared arrays for per type pair parameters + HIP_DYNAMIC_SHARED(char, s_data) + typename evaluator::param_type* s_params = (typename evaluator::param_type*)(&s_data[0]); + Scalar* s_rcutsq + = (Scalar*)(&s_data[num_typ_parameters * sizeof(typename evaluator::param_type)]); + + // load in the per type pair parameters + for (unsigned int cur_offset = 0; cur_offset < num_typ_parameters; cur_offset += blockDim.x) + { + if (cur_offset + threadIdx.x < num_typ_parameters) + { + s_rcutsq[cur_offset + threadIdx.x] = d_rcutsq[cur_offset + threadIdx.x]; + } + } + + unsigned int param_size + = num_typ_parameters * sizeof(typename evaluator::param_type) / sizeof(int); + for (unsigned int cur_offset = 0; cur_offset < param_size; cur_offset += blockDim.x) + { + if (cur_offset + threadIdx.x < param_size) + { + ((int*)s_params)[cur_offset + threadIdx.x] = ((int*)d_params)[cur_offset + threadIdx.x]; + } + } + + __syncthreads(); + + // initialize extra shared mem + char* s_extra = (char*)(s_params + num_typ_parameters); + + unsigned int available_bytes = max_extra_bytes; + for (unsigned int cur_pair = 0; cur_pair < typpair_idx.getNumElements(); ++cur_pair) + s_params[cur_pair].load_shared(s_extra, available_bytes); + + __syncthreads(); + + // start by identifying which particle we are to handle + unsigned int idx; + idx = blockIdx.x * (blockDim.x / tpp) + threadIdx.x / tpp; + bool active = true; + if (idx >= N) + { + // need to mask this thread, but still participate in warp-level reduction + active = false; + } + + // initialize the force to 0 + Scalar4 force = make_scalar4(Scalar(0), Scalar(0), Scalar(0), Scalar(0)); + Scalar4 torque = make_scalar4(Scalar(0), Scalar(0), Scalar(0), Scalar(0)); + Scalar virialxx = Scalar(0); + Scalar virialxy = Scalar(0); + Scalar virialxz = Scalar(0); + Scalar virialyy = Scalar(0); + Scalar virialyz = Scalar(0); + Scalar virialzz = Scalar(0); + + if (active) + { + // load in the length of the neighbor list + unsigned int n_neigh = d_n_neigh[idx]; + + // read in the particle data + Scalar4 postypei = __ldg(d_pos + idx); + Scalar4 veltypei = __ldg(d_vel + idx); + Scalar3 posi = make_scalar3(postypei.x, postypei.y, postypei.z); + Scalar3 veli = make_scalar3(veltypei.x, veltypei.y, veltypei.z); + quat quati(__ldg(d_orientation + idx)); + quat angmomi(__ldg(d_angmom + idx)); + Scalar diai = __ldg(d_diameter + idx); + Scalar3 moment_inertia_i = d_moment_inertia[idx]; + + Scalar qi = Scalar(0); + Scalar massi = Scalar(0); + if (evaluator::needsCharge()) + qi = __ldg(d_charge + idx); + if (evaluator::needsNu()) + massi = veltypei.w; + + // calculate angular momentum of i-th particle in the body frame + quat qxP_i = conj(quati) * angmomi; + vec3 bf_vel_i = Scalar(0.5) + * vec3(qxP_i.v.x / moment_inertia_i.x, + qxP_i.v.y / moment_inertia_i.y, + qxP_i.v.z / moment_inertia_i.z); + + // Rotate angular velocity into global frame + vec3 gf_angvel_i = rotate(quati, bf_vel_i); + Scalar3 angvel_i = make_scalar3(gf_angvel_i.x, gf_angvel_i.y, gf_angvel_i.z); + + size_t my_head = d_head_list[idx]; + unsigned int cur_j = 0; + + unsigned int next_j(0); + next_j = threadIdx.x % tpp < n_neigh ? __ldg(d_nlist + my_head + threadIdx.x % tpp) : 0; + + // loop over neighbors + for (int neigh_idx = threadIdx.x % tpp; neigh_idx < n_neigh; neigh_idx += tpp) + { + { + // read the current neighbor index + // prefetch the next value and set the current one + cur_j = next_j; + if (neigh_idx + tpp < n_neigh) + { + next_j = __ldg(d_nlist + my_head + neigh_idx + tpp); + } + + // get the neighbor's particle data + Scalar4 postypej = __ldg(d_pos + cur_j); + Scalar4 veltypej = __ldg(d_vel + cur_j); + Scalar3 posj = make_scalar3(postypej.x, postypej.y, postypej.z); + Scalar3 velj = make_scalar3(veltypej.x, veltypej.y, veltypej.z); + + quat quatj(__ldg(d_orientation + cur_j)); + quat angmomj(__ldg(d_angmom + cur_j)); + Scalar diaj = __ldg(d_diameter + cur_j); + Scalar3 moment_inertia_j = d_moment_inertia[cur_j]; + + Scalar qj = Scalar(0); + Scalar massj = Scalar(0); + if (evaluator::needsCharge()) + qj = __ldg(d_charge + cur_j); + if (evaluator::needsNu()) + massj = veltypej.w; + + // calculate dv_ij + Scalar3 dv = velj - veli; + + // calculate angular momentum of i-th particle in the body frame + quat qxP_j = conj(quatj) * angmomj; + vec3 bf_vel_j = Scalar(0.5) + * vec3(qxP_j.v.x / moment_inertia_j.x, + qxP_j.v.y / moment_inertia_j.y, + qxP_j.v.z / moment_inertia_j.z); + + // Rotate angular velocity into global frame + vec3 gf_angvel_j = rotate(quatj, bf_vel_j); + Scalar3 angvel_j = make_scalar3(gf_angvel_j.x, gf_angvel_j.y, gf_angvel_j.z); + + // calculate dr (with periodic boundary conditions) + Scalar3 dx = posi - posj; + + // apply periodic boundary conditions + dx = box.minImage(dx); + + // calculate r squared + Scalar rsq = dot(dx, dx); + + // access the per type pair parameters + unsigned int typpair + = typpair_idx(__scalar_as_int(postypei.w), __scalar_as_int(postypej.w)); + Scalar rcutsq = s_rcutsq[typpair]; + const typename evaluator::param_type& param = s_params[typpair]; + + // evaluate the potential + Scalar3 jforce = {Scalar(0), Scalar(0), Scalar(0)}; + Scalar3 torquei = {Scalar(0), Scalar(0), Scalar(0)}; + Scalar3 torquej = {Scalar(0), Scalar(0), Scalar(0)}; + Scalar pair_eng = Scalar(0); + + // constructor call + evaluator eval(dx, angvel_i, angvel_j, dv, diai, diaj, rcutsq, param); + + // Special Potential Pair DPD like Requirements + + // set seed using global tags + + unsigned int tagi = __ldg(d_tag + idx); + unsigned int tagj = __ldg(d_tag + cur_j); + + eval.set_seed_ij_timestep(d_seed, tagi, tagj, d_timestep); + eval.setDeltaT(d_deltaT); + eval.setThirdLaw(third_law); + + if (evaluator::needsCharge()) + eval.setCharge(qi, qj); + if (evaluator::needsTags()) + eval.setTags(__ldg(d_tag + idx), __ldg(d_tag + cur_j)); + + if (evaluator::needsNu()) + { + // Calculate nu for the Ito formalism + Scalar nu_ito = ((Scalar(1.0) / massi) + (Scalar(1.0) / massj)) + + (((diaj * diaj / Scalar(4.0)) / moment_inertia_j.x) + + ((diai * diai / Scalar(4.0)) / moment_inertia_i.x)); + eval.setNu(nu_ito); + } + + // call evaluator + eval.evaluate(jforce, pair_eng, torquei, torquej); + + if (dim == 2) + { + jforce.z = 0.0; + torquei.x = 0.0; + torquei.y = 0.0; + torquej.x = 0.0; + torquej.y = 0.0; + } + + // calculate the virial + if (compute_virial) + { + Scalar3 jforce2 = Scalar(0.5) * jforce; + virialxx += dx.x * jforce2.x; + virialxy += dx.y * jforce2.x; + virialxz += dx.z * jforce2.x; + virialyy += dx.y * jforce2.y; + virialyz += dx.z * jforce2.y; + virialzz += dx.z * jforce2.z; + } + + // add up the force vector components + force.x += jforce.x; + force.y += jforce.y; + force.z += jforce.z; + torque.x += torquei.x; + torque.y += torquei.y; + torque.z += torquei.z; + + force.w += pair_eng; + } + } + + // potential energy per particle must be halved + force.w *= Scalar(0.5); + } + + // reduce force over threads in cta + hoomd::detail::WarpReduce reducer; + force.x = reducer.Sum(force.x); + force.y = reducer.Sum(force.y); + force.z = reducer.Sum(force.z); + force.w = reducer.Sum(force.w); + + torque.x = reducer.Sum(torque.x); + torque.y = reducer.Sum(torque.y); + torque.z = reducer.Sum(torque.z); + + // now that the force calculation is complete, write out the result + if (active && threadIdx.x % tpp == 0) + { + d_force[idx] = force; + d_torque[idx] = torque; + } + + if (compute_virial) + { + virialxx = reducer.Sum(virialxx); + virialxy = reducer.Sum(virialxy); + virialxz = reducer.Sum(virialxz); + virialyy = reducer.Sum(virialyy); + virialyz = reducer.Sum(virialyz); + virialzz = reducer.Sum(virialzz); + + // if we are the first thread in the cta, write out virial to global mem + if (active && threadIdx.x % tpp == 0) + { + d_virial[0 * virial_pitch + idx] = virialxx; + d_virial[1 * virial_pitch + idx] = virialxy; + d_virial[2 * virial_pitch + idx] = virialxz; + d_virial[3 * virial_pitch + idx] = virialyy; + d_virial[4 * virial_pitch + idx] = virialyz; + d_virial[5 * virial_pitch + idx] = virialzz; + } + } + } + +//! Friction pair force compute kernel launcher +/*! + * \tparam evaluator EvaluatorPair class to evaluate V(r) and -delta V(r)/r + * \tparam compute_virial When non-zero, the virial tensor is computed. When zero, the virial tensor + * is not computed. \tparam tpp Number of threads to use per particle, must be power of 2 and + * smaller than warp size + * + * Partial function template specialization is not allowed in C++, so instead we have to wrap this + * with a struct that we are allowed to partially specialize. + */ +template +struct FrictionPairForceComputeKernel + { + //! Launcher for the pair force kernel + /*! + * \param pair_args Other arguments to pass onto the kernel + * \param N Number of particles to compute. + * \param params Parameters for the potential, stored per type pair + */ + + static void launch(const a_pair_args_t& pair_args, + unsigned int N, + const typename evaluator::param_type* params) + { + if (tpp == pair_args.threads_per_particle) + { + unsigned int block_size = pair_args.block_size; + + Index2D typpair_idx(pair_args.ntypes); + size_t shared_bytes = (2 * sizeof(Scalar) + sizeof(typename evaluator::param_type)) + * typpair_idx.getNumElements(); + + unsigned int max_block_size; + hipFuncAttributes attr; + hipFuncGetAttributes( + &attr, + reinterpret_cast( + &gpu_compute_pair_friction_forces_kernel)); + int max_threads = attr.maxThreadsPerBlock; + // number of threads has to be multiple of warp size + max_block_size = max_threads - max_threads % gpu_friction_pair_force_max_tpp; + + unsigned int base_shared_bytes; + base_shared_bytes = (unsigned int)(shared_bytes + attr.sharedSizeBytes); + + if (base_shared_bytes > pair_args.devprop.sharedMemPerBlock) + { + throw std::runtime_error("Pair potential parameters exceed the available shared " + "memory per block."); + } + + unsigned int max_extra_bytes + = (unsigned int)(pair_args.devprop.sharedMemPerBlock - base_shared_bytes); + unsigned int extra_bytes; + // determine dynamically requested shared memory + char* ptr = (char*)nullptr; + unsigned int available_bytes = max_extra_bytes; + for (unsigned int i = 0; i < typpair_idx.getNumElements(); ++i) + { + params[i].allocate_shared(ptr, available_bytes); + } + extra_bytes = max_extra_bytes - available_bytes; + + shared_bytes += extra_bytes; + + block_size = block_size < max_block_size ? block_size : max_block_size; + dim3 grid(N / (block_size / tpp) + 1, 1, 1); + + hipLaunchKernelGGL( + (gpu_compute_pair_friction_forces_kernel), + dim3(grid), + dim3(block_size), + shared_bytes, + 0, + pair_args.d_force, + pair_args.d_torque, + pair_args.d_virial, + pair_args.virial_pitch, + N, + pair_args.d_pos, + pair_args.d_vel, + pair_args.d_charge, + pair_args.d_orientation, + pair_args.d_angmom, + pair_args.d_diameter, + pair_args.d_moment_inertia, + pair_args.d_tag, + pair_args.box, + pair_args.third_law, + pair_args.dim, + pair_args.seed, + pair_args.timestep, + pair_args.deltaT, + pair_args.d_n_neigh, + pair_args.d_nlist, + pair_args.d_head_list, + params, + pair_args.d_rcutsq, + pair_args.ntypes, + max_extra_bytes); + } + else + { + FrictionPairForceComputeKernel::launch(pair_args, + N, + params); + } + } + }; + +//! Template specialization to do nothing for the tpp = 0 case +template +struct FrictionPairForceComputeKernel + { + static void launch(const a_pair_args_t& pair_args, + unsigned int N, + const typename evaluator::param_type* d_params) + { + // do nothing + } + }; + +//! Kernel driver that computes lj forces on the GPU for FrictionPairGPU +/*! \param pair_args Other arguments to pass onto the kernel + \param d_params Parameters for the potential, stored per type pair + + This is just a driver function for gpu_compute_pair_friction_forces_kernel(), see it for + details. +*/ +template +hipError_t gpu_compute_pair_friction_forces(const a_pair_args_t& pair_args, + const typename evaluator::param_type* d_params) + { + assert(d_params); + assert(pair_args.d_rcutsq); + assert(pair_args.ntypes > 0); + + // run the kernel + if (pair_args.compute_virial) + { + FrictionPairForceComputeKernel::launch( + pair_args, + pair_args.N, + d_params); + } + else + { + FrictionPairForceComputeKernel::launch( + pair_args, + pair_args.N, + d_params); + } + return hipSuccess; + } +#else +template +hipError_t gpu_compute_pair_friction_forces(const a_pair_args_t& pair_args, + const typename evaluator::param_type* d_params); +#endif + + } // end namespace kernel + } // end namespace md + } // end namespace hoomd + +#endif // __FRICTION_PAIR_GPU_CUH__ diff --git a/hoomd/md/FrictionPairGPU.h b/hoomd/md/FrictionPairGPU.h new file mode 100644 index 0000000000..ce26e9b1e9 --- /dev/null +++ b/hoomd/md/FrictionPairGPU.h @@ -0,0 +1,226 @@ +// Copyright (c) 2009-2025 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +#ifndef __FRICTION_PAIR_GPU_H__ +#define __FRICTION_PAIR_GPU_H__ + +#ifdef ENABLE_HIP + +#include "FrictionPair.h" +#include "FrictionPairGPU.cuh" +#include "hoomd/Autotuner.h" + +/*! \file FrictionPairGPU.h + \brief Defines the template class for frictional pair interactions on the GPU + \note This header cannot be compiled by nvcc +*/ + +#ifdef __HIPCC__ +#error This header cannot be compiled by nvcc +#endif + +#include + +namespace hoomd + { +namespace md + { +//! Template class for computing friction pair potentials on the GPU +/*! Derived from FrictionPair, this class provides exactly the same interface for computing + friction pair interactions, forces and torques. In the same way as FrictionPair, this class + serves as a shell dealing with all the details common to every frictional contact calculation + while te \a evaluator calculates \f$V(\vec r,\vec e_i, \vec e_j)\f$ in a generic way. + + \tparam evaluator EvaluatorPair class used to evaluate potential, force and torque. + \sa export_FrictionPairGPU() +*/ +template class FrictionPairGPU : public FrictionPair + { + public: + //! Construct the pair potential + FrictionPairGPU(std::shared_ptr sysdef, std::shared_ptr nlist); + //! Destructor + virtual ~FrictionPairGPU() { }; + + virtual void + setParams(unsigned int typ1, unsigned int typ2, const typename evaluator::param_type& param); + + protected: + std::shared_ptr> m_tuner; //!< Autotuner for block size and threads per particle + + //! Actually compute the forces + virtual void computeForces(uint64_t timestep); + }; + +template +FrictionPairGPU::FrictionPairGPU(std::shared_ptr sysdef, + std::shared_ptr nlist) + : FrictionPair(sysdef, nlist) + { + // can't run on the GPU if there aren't any GPUs in the execution configuration + if (!this->m_exec_conf->isCUDAEnabled()) + { + this->m_exec_conf->msg->error() + << "ai_pair." << evaluator::getName() + << ": Creating a FrictionPairGPU with no GPU in the execution configuration" + << std::endl + << std::endl; + throw std::runtime_error("Error initializing FrictionPairGPU"); + } + + // Initialize autotuner that tunes block sizes and threads per particle. + m_tuner.reset(new Autotuner<2>({AutotunerBase::makeBlockSizeRange(this->m_exec_conf), + AutotunerBase::getTppListPow2(this->m_exec_conf)}, + this->m_exec_conf, + "friction_pair_" + evaluator::getName())); + this->m_autotuners.push_back(m_tuner); + +#ifdef ENABLE_MPI + // synchronize autotuner results across ranks + m_tuner->setSync(bool(this->m_pdata->getDomainDecomposition())); +#endif + } + +template void FrictionPairGPU::computeForces(uint64_t timestep) + { + this->m_nlist->compute(timestep); + + // The GPU implementation CANNOT handle a half neighborlist, error out now + bool third_law = this->m_nlist->getStorageMode() == NeighborList::half; + if (third_law) + { + this->m_exec_conf->msg->error() + << "ai_pair." << evaluator::getName() + << ": FrictionPairGPU cannot handle a half neighborlist" << std::endl + << std::endl; + throw std::runtime_error("Error computing forces in FrictionPairGPU"); + } + + // access the neighbor list + ArrayHandle d_n_neigh(this->m_nlist->getNNeighArray(), + access_location::device, + access_mode::read); + ArrayHandle d_nlist(this->m_nlist->getNListArray(), + access_location::device, + access_mode::read); + ArrayHandle d_head_list(this->m_nlist->getHeadList(), + access_location::device, + access_mode::read); + + // access the particle data + ArrayHandle d_pos(this->m_pdata->getPositions(), + access_location::device, + access_mode::read); + ArrayHandle d_vel(this->m_pdata->getVelocities(), + access_location::device, + access_mode::read); + ArrayHandle d_charge(this->m_pdata->getCharges(), + access_location::device, + access_mode::read); + ArrayHandle d_orientation(this->m_pdata->getOrientationArray(), + access_location::device, + access_mode::read); + ArrayHandle d_angmom(this->m_pdata->getAngularMomentumArray(), + access_location::device, + access_mode::read); + ArrayHandle d_diameter(this->m_pdata->getDiameters(), + access_location::device, + access_mode::read); + ArrayHandle d_moment_inertia(this->m_pdata->getMomentsOfInertiaArray(), + access_location::device, + access_mode::read); + ArrayHandle d_tag(this->m_pdata->getTags(), + access_location::device, + access_mode::read); + + BoxDim box = this->m_pdata->getBox(); + unsigned int dim = this->m_sysdef->getNDimensions(); + + // uint16_t seed = this->m_sysdef->getSeed(); + // Scalar deltaT = this->m_deltaT; + + // access parameters + ArrayHandle d_rcutsq(this->m_rcutsq, access_location::device, access_mode::read); + ArrayHandle d_force(this->m_force, access_location::device, access_mode::overwrite); + ArrayHandle d_torque(this->m_torque, access_location::device, access_mode::overwrite); + ArrayHandle d_virial(this->m_virial, access_location::device, access_mode::overwrite); + + // access flags + PDataFlags flags = this->m_pdata->getFlags(); + + this->m_exec_conf->setDevice(); + + this->m_tuner->begin(); + unsigned int block_size = this->m_tuner->getParam()[0]; + unsigned int threads_per_particle = this->m_tuner->getParam()[1]; + + kernel::gpu_compute_pair_friction_forces( + kernel::a_pair_args_t(d_force.data, + d_torque.data, + d_virial.data, + this->m_virial.getPitch(), + this->m_pdata->getN(), + this->m_pdata->getMaxN(), + d_pos.data, + d_vel.data, + d_charge.data, + d_orientation.data, + d_angmom.data, + d_diameter.data, + d_moment_inertia.data, + d_tag.data, + box, + third_law, + dim, + this->m_sysdef->getSeed(), + timestep, + this->m_deltaT, + d_n_neigh.data, + d_nlist.data, + d_head_list.data, + d_rcutsq.data, + this->m_pdata->getNTypes(), + block_size, + flags[pdata_flag::pressure_tensor], + threads_per_particle, + this->m_exec_conf->dev_prop), + this->m_params.data()); + + this->m_tuner->end(); + + if (this->m_exec_conf->isCUDAErrorCheckingEnabled()) + CHECK_CUDA_ERROR(); + } + +template +void FrictionPairGPU::setParams(unsigned int typ1, + unsigned int typ2, + const typename evaluator::param_type& param) + { + FrictionPair::setParams(typ1, typ2, param); + this->m_params[this->m_typpair_idx(typ1, typ2)].set_memory_hint(); + this->m_params[this->m_typpair_idx(typ2, typ1)].set_memory_hint(); + } + +namespace detail + { +//! Export this pair potential to python +/*! \param name Name of the class in the exported python module + \tparam T Class type to export. \b Must be an instantiated FrictionPairGPU class template. + \tparam Base Base class of \a T. \b Must be Pair with the same evaluator as + used in \a T. +*/ +template void export_FrictionPairGPU(pybind11::module& m, const std::string& name) + { + pybind11::class_, FrictionPair, std::shared_ptr>>( + m, + name.c_str()) + .def(pybind11::init, std::shared_ptr>()); + } + + } // end namespace detail + } // end namespace md + } // end namespace hoomd + +#endif // ENABLE_HIP +#endif // __FRICTION_PAIR_GPU_H__ diff --git a/hoomd/md/module-md.cc b/hoomd/md/module-md.cc index db2f64fe7b..ea168671f1 100644 --- a/hoomd/md/module-md.cc +++ b/hoomd/md/module-md.cc @@ -84,6 +84,9 @@ void export_AnisoPotentialPairPatchyMie(pybind11::module& m); void export_AnisoPotentialPairPatchyYukawa(pybind11::module& m); void export_AnisoPotentialPairPatchyTable(pybind11::module& m); +void export_FrictionPairFrictionLJLinear(pybind11::module& m); +void export_FrictionPairFrictionLJCoulomb(pybind11::module& m); +void export_FrictionPairFrictionLJCoulombNewton(pybind11::module& m); void export_MeshForceCompute(pybind11::module& m); void export_PotentialBondHarmonic(pybind11::module& m); @@ -250,6 +253,10 @@ void export_AnisoPotentialPairPatchyMieGPU(pybind11::module& m); void export_AnisoPotentialPairPatchyYukawaGPU(pybind11::module& m); void export_AnisoPotentialPairPatchyTableGPU(pybind11::module& m); +void export_FrictionPairFrictionLJLinearGPU(pybind11::module& m); +void export_FrictionPairFrictionLJCoulombGPU(pybind11::module& m); +void export_FrictionPairFrictionLJCoulombNewtonGPU(pybind11::module& m); + void export_PotentialBondHarmonicGPU(pybind11::module& m); void export_PotentialBondFENEGPU(pybind11::module& m); void export_PotentialBondTetherGPU(pybind11::module& m); @@ -399,6 +406,10 @@ PYBIND11_MODULE(_md, m) export_AnisoPotentialPairPatchyYukawa(m); export_AnisoPotentialPairPatchyTable(m); + export_FrictionPairFrictionLJLinear(m); + export_FrictionPairFrictionLJCoulomb(m); + export_FrictionPairFrictionLJCoulombNewton(m); + export_PotentialPairDPDThermoDPD(m); export_PotentialPairDPDThermoLJ(m); @@ -493,6 +504,10 @@ PYBIND11_MODULE(_md, m) export_AnisoPotentialPairYLZGPU(m); export_AnisoPotentialPairGBGPU(m); + export_FrictionPairFrictionLJLinearGPU(m); + export_FrictionPairFrictionLJCoulombGPU(m); + export_FrictionPairFrictionLJCoulombNewtonGPU(m); + export_AnisoPotentialPairPatchyExpandedGaussianGPU(m); ; export_AnisoPotentialPairPatchyExpandedLJGPU(m); diff --git a/hoomd/md/pair/CMakeLists.txt b/hoomd/md/pair/CMakeLists.txt index 38e79960de..c46dedb93e 100644 --- a/hoomd/md/pair/CMakeLists.txt +++ b/hoomd/md/pair/CMakeLists.txt @@ -1,6 +1,7 @@ set(files __init__.py pair.py aniso.py + friction.py ) install(FILES ${files} diff --git a/hoomd/md/pair/__init__.py b/hoomd/md/pair/__init__.py index ee09159bce..3a0356728b 100644 --- a/hoomd/md/pair/__init__.py +++ b/hoomd/md/pair/__init__.py @@ -131,6 +131,7 @@ """ from . import aniso +from . import friction from .pair import ( Pair, LJ, @@ -192,4 +193,5 @@ "Yukawa", "Zetterling", "aniso", + "friction", ] diff --git a/hoomd/md/pair/friction.py b/hoomd/md/pair/friction.py new file mode 100644 index 0000000000..dc846d156e --- /dev/null +++ b/hoomd/md/pair/friction.py @@ -0,0 +1,425 @@ +# Copyright (c) 2009-2025 The Regents of the University of Michigan. +# Part of HOOMD-blue, released under the BSD 3-Clause License. + +r"""Frictional pair force classes apply a force, and torque on every +particle in the simulation state. The following general expression +for Markovian tangential friction forces is implemented for +interactions between two spherical particles and is discussed in detail in +`Hofmann et. al. 2025`_. For two particles :math:`i` and +:math:`j` with radii :math:`R_{i,j}`, center positions :math:`\mathbf{r}_{i,j}`, +angular velocities :math:`\mathbf{\omega}_{i,j}`, and translational velocities +:math:`\mathbf{v}_{i,j}`, their surface velocities at the contact point are given by + +.. _Hofmann et. al. 2025: https://doi.org/10.48550/arXiv.2507.16388 + +.. math:: + + \begin{align*} + \mathbf{u}_i &= \mathbf{v}_i+\mathbf{\omega}_i \times \mathbf{\hat{r}}_{ij}R_i \\ + \mathbf{u}_j &= \mathbf{v}_j-\mathbf{\omega}_j \times \mathbf{\hat{r}}_{ij}R_j \, , + \end{align*} + +where :math:`\mathbf{\hat{r}}_{ij}=\mathbf{r}_{ij}/r_{i,j}`. With these expressions, +we calculate the relative tangential velocity :math:`\mathbf{u}^\perp_{i,j}` at the +contact point + +.. math:: + + \mathbf{u}^\perp_{i,j} = \mathbf{P}(\mathbf{\hat{r}}_{ij})(\mathbf{v}_j + -\mathbf{v}_i) -(\mathbf{\omega}_iR_i+\mathbf{\omega}_jR_j) + \times \mathbf{\hat{r}}_{ij}\, , + +wit the projection operator :math:`\mathbf{P}(\mathbf{\hat{r}}_{ij})=1 +-\mathbf{\hat{r}}_{ij}\mathbf{\hat{r}}_{ij}`. We model the tangential +friction force at the contact point very general as + +.. math:: + + \mathbf{F}^\mathrm{f,contact}_i = -\mathbf{F}^\mathrm{f,contact}_j = f(u^\perp_{i,j} + ,r_{i,j})\mathbf{\hat{u}}^\perp_{i,j} + +with :math:`\mathbf{\hat{u}}^\perp_{i,j}=\mathbf{u}^\perp_{i,j}/u^\perp_{i,j}`, where +:math:`f(u^\perp_{i,j},r_{i,j})` is an arbitrary function. The surface force +:math:`\mathbf{F}^\mathrm{f,contact}_i` generates a center-of-mass force and torque +acting on particle :math:`i`, + +.. math:: + + \begin{align*} + \mathbf{F}^\mathrm{f}_{ij} &= \mathbf{F}^\mathrm{f,contact}_i \\ + \mathbf{\tau}^\mathrm{f}_{ij} &= R_i\mathbf{\hat{r}}_{ij} \times + \mathbf{F}^\mathrm{f,contact}_i\, , + \end{align*} + +which is a pair friction force and torque resulting from the friction with the particle +:math:`j`. + +The functional form of :math:`f(u^\perp_{i,j},r_{i,j})` specifies the frictional model. + +`FrictionalPair` does not support any shifting modes. + +.. invisible-code-block: python + + neighbor_list = hoomd.md.nlist.Cell(buffer = 0.4) + simulation = hoomd.util.make_example_simulation() + simulation.operations.integrator = hoomd.md.Integrator( + dt=0.001, + integrate_rotational_dof = True) +""" + +import inspect + +from hoomd.md.pair.pair import Pair +from hoomd.data.parameterdicts import TypeParameterDict +from hoomd.data.typeparam import TypeParameter + + +class FrictionalPair(Pair): + r"""Base class friction pair force. + + `FrictionalPair` is the base class for all frictional pair forces. + + Warning: + This class should not be instantiated by users. The class can be used + for `isinstance` or `issubclass` checks. + """ + + __doc__ = inspect.cleandoc(__doc__) + "\n\n" + inspect.cleandoc(Pair._doc_inherited) + + _accepted_modes = ("none",) + + def __init__(self, nlist, default_r_cut=None, mode="none"): + super().__init__(nlist, default_r_cut, 0.0, mode) + + +class FrictionLJCoulomb(FrictionalPair): + r"""Coulomb frictional model pair force with the LJ conservative force. + + Args: + nlist (hoomd.md.nlist.NeighborList): Neighbor list + default_r_cut (float): Default cutoff radius :math:`[\mathrm{length}]`. + + `FrictionLJCoulomb` computes the frictional interaction + between pairs of particles with a Coulomb friction model + as described in `Hofmann et. al. 2025`_. + + .. _Hofmann et. al. 2025: https://doi.org/10.48550/arXiv.2507.16388 + + The Coulomb friction model is defined by the function + + .. math:: + + f_\mathrm{C}(u,r) = w(r)\kappa_\mathrm{f}\, . + + And a repulsive Weeks-Chandler-Anderson (WCA) potential + + .. math:: + + \begin{align*} + U_\mathrm{WCA}(r) &= 4 \varepsilon \left[ \left( + \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} + \right)^{6} \right] \\ + w(r)&=-\frac{\mathrm{d}}{\mathrm{d}r}U_\mathrm{WCA}(r) + \end{align*} + + Since the frictional force results from coarse-graining the microscopic equations + of motion, a corresponding pairwise noise that satisfies a fluctuation-dissipation + relation must be applied. For this a stochastic force and torque is added + + .. math:: + + \begin{align*} + D_\mathrm{C}(u,r) &= \frac{w(r)\kappa_\mathrm{f}\sqrt{\pi}}{\sqrt{2k_\mathrm{B} + T\nu}}\mathrm{e}^{\frac{u^2}{2k_\mathrm{B}T\nu}} + \mathrm{Erfc}\big(\frac{u} + {\sqrt{2k_\mathrm{B}T\nu}} \big) \\ + \mathbf{F}^\mathrm{R}_{ij} = -\mathbf{F}^\mathrm{R}_{ji} &= \sqrt{D_\mathrm{C} + (u^\perp_{ij},r_{ij})}\Big[\mathbf{P} + (\mathbf{\hat{r}}_{ij})\mathbf{\xi}^\mathrm{f}_{ij} + - \mathbf{\hat{r}}_{ij} \times \mathbf{N}^\mathrm{f} + _{ij}\Big] \\ + \frac{\mathbf{\tau}^\mathrm{R}_{ij}}{R_i} = + \frac{\mathbf{\tau}^\mathrm{R}_{ji}}{R_j} &= + \sqrt{D_\mathrm{C}(u^\perp_{ij},r_{ij})}\Big[ + \mathbf{\hat{r}}_{ij} \times + \mathbf{\xi}^\mathrm{f}_{ij}+\mathbf{P} + (\mathbf{e}_{ij})\mathbf{N}^\mathrm{f}_{ij}\Big]\, , + \end{align*} + + where :math:`\nu=(1/m_i+1/m_j)+(R_i^2/\mathbf{I}_i+R_j^2/\mathbf{I}_j))`, + and :math:`\mathbf{\xi}^\mathrm{f}_{ij}` and :math:`\mathbf{N}^\mathrm{f}_{ij}` + are three dimensional Gaussian white noise vectors with correlations + + .. math:: + + \begin{align*} + \langle \mathbf{\xi}^\mathrm{f}_{ij}(t) \mathbf{\xi}^\mathrm{f}_{kl}(t') + \rangle &= \mathbf{1}k_\mathrm{B}T + (\delta_{ik}\delta_{jl}-\delta_{il} + \delta_{jk})\delta(t-t') \\ + \langle \mathbf{N}^\mathrm{f}_{ij}(t) \mathbf{N}^\mathrm{f}_{kl}(t') \rangle &= + \mathbf{1}k_\mathrm{B}T(\delta_{ik} + \delta_{jl}+\delta_{il}\delta_{jk}) + \delta(t-t')\, . + \end{align*} + + .. rubric:: Example: + + .. code-block:: python + + coulomb_lj = hoomd.md.pair.friction.FrictionLJCoulomb(nlist=neighbor_list, + default_r_cut=3) + + coulomb_lj_params = {'epsilon':1, 'sigma':1, 'kappa_f':3, 'kT':1} + + coulomb_lj.params.default = coulomb_lj_params + simulation.operations.integrator.forces = [coulomb_lj] + + {inherited} + """ + + _cpp_class_name = "FrictionPairFrictionLJCoulomb" + __doc__ = inspect.cleandoc(__doc__).replace( + "{inherited}", inspect.cleandoc(FrictionalPair._doc_inherited) + ) + + def __init__(self, nlist, default_r_cut=None): + super().__init__(nlist, default_r_cut, "none") + params = TypeParameter( + "params", + "particle_types", + TypeParameterDict( + epsilon=float, sigma=float, kappa_f=float, kT=float, len_keys=2 + ), + ) + self._add_typeparam(params) + + +class FrictionLJCoulombNewton(FrictionalPair): + r"""CoulombNewton frictional model pair force with the LJ conservative force. + + Args: + nlist (hoomd.md.nlist.NeighborList): Neighbor list + default_r_cut (float): Default cutoff radius :math:`[\mathrm{length}]`. + + `FrictionLJCoulombNewton` computes the frictional interaction + between pairs of particles with a Coulomb-Newton friction model as described in + `Hofmann et. al. 2025`_. + + .. _Hofmann et. al. 2025: https://doi.org/10.48550/arXiv.2507.16388 + + The Coulomb-Newton friction model is defined by the function + + .. math:: + + f_\mathrm{CN}(u,r) = \mathrm{min}[\gamma_\mathrm{f}u,w(r)\kappa_\mathrm{f}]\, . + + And a repulsive Weeks-Chandler-Anderson (WCA) potential + + .. math:: + + \begin{align*} + U_\mathrm{WCA}(r) &= 4 \varepsilon \left[ \left( + \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} + \right)^{6} \right] \\ + w(r)&=-\frac{\mathrm{d}}{\mathrm{d}r}U_\mathrm{WCA}(r) + \end{align*} + + Since the frictional force results from coarse-graining the microscopic equations + of motion, a corresponding pairwise noise that satisfies a fluctuation-dissipation + relation must be applied. For this a stochastic force and torque is added + + .. math:: + + \begin{align*} + D_\mathrm{CN}(u,r) &= \begin{cases} + \frac{w(r)\kappa_\mathrm{f}\sqrt{\pi}} + {\sqrt{2k_\mathrm{B}T\nu}}\mathrm{e}^{\frac{u^2} + {2k_\mathrm{B}T\nu}}\mathrm{Erfc}\big(\frac{u} + {\sqrt{2k_\mathrm{B}T\nu}} \big) & ,\, u\ge + \frac{w(r)\kappa_\mathrm{f}}{\gamma_\mathrm{f}} \\ + \gamma_\mathrm{f}\left(1-\mathrm{e}^{\left(u^2 + -(\frac{w(r)\kappa_\mathrm{f}}{\gamma_\mathrm{f}})^2 + \right)/2k_\mathrm{B}T\nu} \right) + \frac{w(r) + \kappa_\mathrm{f}\sqrt{\pi}}{\sqrt{2k_\mathrm{B}T\nu}} + \mathrm{e}^\frac{u^2}{2k_\mathrm{B}T\nu}\mathrm{Erfc} + \left(\frac{w(r)\kappa_\mathrm{f}/\gamma_\mathrm{f}} + {\sqrt{2k_\mathrm{B}T\nu}}\right) &,\, u < \frac{w(r) + \kappa_\mathrm{f}}{\gamma_\mathrm{f}} + \end{cases} \\ + \mathbf{F}^\mathrm{R}_{ij} = -\mathbf{F}^\mathrm{R}_{ji} &= \sqrt{D_\mathrm{CN} + (u^\perp_{ij},r_{ij})}\Big[\mathbf{P} + (\mathbf{\hat{r}}_{ij}) + \mathbf{\xi}^\mathrm{f}_{ij} - + \mathbf{\hat{r}}_{ij} \times + \mathbf{N}^\mathrm{f}_{ij}\Big] \\ + \frac{\mathbf{\tau}^\mathrm{R}_{ij}}{R_i} = \frac{\mathbf{\tau}^\mathrm{R}_{ji}} + {R_j} &= \sqrt{D_\mathrm{CN} + (u^\perp_{ij},r_{ij})}\Big[\mathbf{\hat{r}}_{ij} + \times \mathbf{\xi}^\mathrm{f}_{ij}+\mathbf{P} + (\mathbf{e}_{ij})\mathbf{N}^\mathrm{f}_{ij} + \Big]\, , + \end{align*} + + where :math:`\nu=(1/m_i+1/m_j)+(R_i^2/\mathbf{I}_i+R_j^2/\mathbf{I}_j))`, and + :math:`\mathbf{\xi}^\mathrm{f}_{ij}` and :math:`\mathbf{N}^\mathrm{f}_{ij}` are + three dimensional Gaussian white noise vectors with correlations + + .. math:: + + \begin{align*} + \langle \mathbf{\xi}^\mathrm{f}_{ij}(t) \mathbf{\xi}^\mathrm{f}_{kl}(t') + \rangle &= \mathbf{1}k_\mathrm{B}T + (\delta_{ik}\delta_{jl}-\delta_{il} + \delta_{jk})\delta(t-t') \\ + \langle \mathbf{N}^\mathrm{f}_{ij}(t) \mathbf{N}^\mathrm{f}_{kl}(t') \rangle &= + \mathbf{1}k_\mathrm{B}T + (\delta_{ik}\delta_{jl}+\delta_{il} + \delta_{jk})\delta(t-t')\, . + \end{align*} + + .. rubric:: Example: + + .. code-block:: python + + coulombNewton_lj = hoomd.md.pair.friction.FrictionLJCoulombNewton( + nlist=neighbor_list, + default_r_cut=3) + + coulombNewton_lj_params = { 'epsilon':1, 'sigma':1, 'gamma_f':1, 'kappa_f':3, + 'kT':1} + + coulombNewton_lj.params.default = coulombNewton_lj_params + simulation.operations.integrator.forces = [coulombNewton_lj] + + {inherited} + """ + + _cpp_class_name = "FrictionPairFrictionLJCoulombNewton" + __doc__ = inspect.cleandoc(__doc__).replace( + "{inherited}", inspect.cleandoc(FrictionalPair._doc_inherited) + ) + + def __init__(self, nlist, default_r_cut=None): + super().__init__(nlist, default_r_cut, "none") + params = TypeParameter( + "params", + "particle_types", + TypeParameterDict( + epsilon=float, + sigma=float, + gamma_f=float, + kappa_f=float, + kT=float, + len_keys=2, + ), + ) + self._add_typeparam(params) + + +class FrictionLJLinear(FrictionalPair): + r"""Linear frictional model pair force with the LJ conservative force. + + Args: + nlist (hoomd.md.nlist.NeighborList): Neighbor list + default_r_cut (float): Default cutoff radius :math:`[\mathrm{length}]`. + + `FrictionLJLinear` computes the frictional interaction + between pairs of particles with a linear friction model as described in + `Hofmann et. al. 2025`_. + + .. _Hofmann et. al. 2025: https://doi.org/10.48550/arXiv.2507.16388 + + The linear friction model is defined by the function + + .. math:: + + f_\mathrm{l}(u,r) = w(r)\gamma_\mathrm{f}u\, . + + And a repulsive Weeks-Chandler-Anderson (WCA) potential + + .. math:: + + \begin{align*} + U_\mathrm{WCA}(r) &= 4 \varepsilon \left[ \left( + \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} + \right)^{6} \right] \\ + w(r)&=-\frac{\mathrm{d}}{\mathrm{d}r}U_\mathrm{WCA}(r) + \end{align*} + + Since the frictional force results from coarse-graining the microscopic equations + of motion, a corresponding pairwise noise that satisfies a fluctuation-dissipation + relation must be applied. For this a stochastic force and torque is added + + .. math:: + + \begin{align*} + D_\mathrm{l}(u,r) &= w(r)\gamma_\mathrm{f} \\ + \mathbf{F}^\mathrm{R}_{ij} = -\mathbf{F}^\mathrm{R}_{ji} &= \sqrt{D_\mathrm{l} + (u^\perp_{ij},r_{ij})}\Big[\mathbf{P} + (\mathbf{\hat{r}}_{ij})\mathbf{\xi}^\mathrm{f}_{ij} + - \mathbf{\hat{r}}_{ij} \times + \mathbf{N}^\mathrm{f}_{ij}\Big] \\ + \frac{\mathbf{\tau}^\mathrm{R}_{ij}}{R_i} = \frac{\mathbf{\tau}^\mathrm{R}_{ji}} + {R_j} &= \sqrt{D_\mathrm{l} + (u^\perp_{ij},r_{ij})}\Big[\mathbf{\hat{r}}_{ij} + \times \mathbf{\xi}^\mathrm{f}_{ij}+\mathbf{P} + (\mathbf{e}_{ij})\mathbf{N}^\mathrm{f}_{ij} + \Big]\, , + \end{align*} + + where :math:`\nu=(1/m_i+1/m_j)+(R_i^2/\mathbf{I}_i+R_j^2/\mathbf{I}_j))`, and + :math:`\mathbf{\xi}^\mathrm{f}_{ij}` and :math:`\mathbf{N}^\mathrm{f}_{ij}` are + three dimensional Gaussian white noise vectors with correlations + + .. math:: + + \begin{align*} + \langle \mathbf{\xi}^\mathrm{f}_{ij}(t) \mathbf{\xi}^\mathrm{f}_{kl}(t') + \rangle &= \mathbf{1}k_\mathrm{B}T + (\delta_{ik}\delta_{jl}-\delta_{il} + \delta_{jk})\delta(t-t') \\ + \langle \mathbf{N}^\mathrm{f}_{ij}(t) \mathbf{N}^\mathrm{f}_{kl}(t') \rangle &= + \mathbf{1}k_\mathrm{B}T(\delta_{ik} + \delta_{jl}+\delta_{il}\delta_{jk}) + \delta(t-t')\, . + \end{align*} + + .. rubric:: Example: + + .. code-block:: python + + linear_lj = hoomd.md.pair.friction.FrictionLJLinear(nlist=neighbor_list, + default_r_cut=3) + + linear_lj_params = {'epsilon':1, 'sigma':1, 'gamma_f':1, 'kT':1} + + linear_lj.params.default = linear_lj_params + simulation.operations.integrator.forces = [linear_lj] + + {inherited} + """ + + _cpp_class_name = "FrictionPairFrictionLJLinear" + __doc__ = inspect.cleandoc(__doc__).replace( + "{inherited}", inspect.cleandoc(FrictionalPair._doc_inherited) + ) + + def __init__(self, nlist, default_r_cut=None): + super().__init__(nlist, default_r_cut, "none") + params = TypeParameter( + "params", + "particle_types", + TypeParameterDict( + epsilon=float, sigma=float, gamma_f=float, kT=float, len_keys=2 + ), + ) + self._add_typeparam(params) + + +__all__ = [ + "FrictionLJCoulomb", + "FrictionLJCoulombNewton", + "FrictionLJLinear", + "FrictionalPair", +] diff --git a/hoomd/md/pytest/test_friction.py b/hoomd/md/pytest/test_friction.py new file mode 100644 index 0000000000..f742be1e7e --- /dev/null +++ b/hoomd/md/pytest/test_friction.py @@ -0,0 +1,250 @@ +# Copyright (c) 2009-2025 The Regents of the University of Michigan. +# Part of HOOMD-blue, released under the BSD 3-Clause License. + +import hoomd +import pytest +import numpy + +TOLERANCES = {"rtol": 1e-2, "atol": 1e-5} + +friction_test_parameters = [ + # test for LJLinear friction model + ( + hoomd.md.pair.friction.FrictionLJLinear, + {}, + { + "epsilon": 1, + "sigma": 1, + "gamma_f": 6, + "kT": 1, + }, + [[-0.35, -0.35, -0.1], [0.35, 0.35, 0.1]], + [[0.8, 0.0, 0.0], [-0.8, 0.0, 0.0]], + [[0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 0.0, 1.0]], + [ + [0.19503148, -0.97515738, 0.0390063, 0.09751574], + [0.81649658, 0.0, -0.40824829, -0.40824829], + ], + [38.39953386, -152.42732975, 299.06283406], + [118.73327022, -99.83859521, -66.13136255], + ), + # test for the LJCoulomb friction model + ( + hoomd.md.pair.friction.FrictionLJCoulomb, + {}, + { + "epsilon": 1, + "sigma": 1, + "kappa_f": 5, + "kT": 1, + }, + [[-0.35, -0.35, -0.1], [0.35, 0.35, 0.1]], + [[0.8, 0.0, 0.0], [-0.8, 0.0, 0.0]], + [[0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 0.0, 1.0]], + [ + [0.19503148, -0.97515738, 0.0390063, 0.09751574], + [0.81649658, 0.0, -0.40824829, -0.40824829], + ], + [1.57890627, -54.46188002, 85.05595653], + [34.8688109, -29.31994622, -19.4210264], + ), + # test for the LJCoulombNewton friction model + ( + hoomd.md.pair.friction.FrictionLJCoulombNewton, + {}, + { + "epsilon": 1, + "sigma": 1, + "gamma_f": 6, + "kappa_f": 5, + "kT": 1, + }, + [[-0.35, -0.35, -0.1], [0.35, 0.35, 0.1]], + [[0.8, 0.0, 0.0], [-0.8, 0.0, 0.0]], + [[0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 0.0, 1.0]], + [ + [0.19503148, -0.97515738, 0.0390063, 0.09751574], + [0.81649658, 0.0, -0.40824829, -0.40824829], + ], + [-11.09870197, -20.73166752, 11.37184162], + [5.99367133, -5.03986561, -3.33832001], + ), +] + + +@pytest.fixture(scope="session") +def friction_snapshot_factory(device): + def make_snapshot( + position_i=numpy.array([[]]), + position_j=numpy.array([]), + velocity_i=numpy.array([]), + velocity_j=numpy.array([]), + angmom_i=numpy.array([]), + angmom_j=numpy.array([]), + orientation_i=(), + orientation_j=(), + dimensions=3, + L=20, + ): + snapshot = hoomd.Snapshot(device.communicator) + if snapshot.communicator.rank == 0: + N = 2 + box = [L, L, L, 0, 0, 0] + if dimensions == 2: + box[2] = 0 + snapshot.configuration.box = box + snapshot.particles.N = N + snapshot.particles.diameter[:] = [1] * N + snapshot.particles.mass[:] = [1] * N + snapshot.particles.position[:] = [position_i, position_j] + snapshot.particles.velocity[:] = [velocity_i, velocity_j] + snapshot.particles.angmom[:] = [angmom_i, angmom_j] + snapshot.particles.orientation[:] = [orientation_i, orientation_j] + snapshot.particles.types = ["A", "B"] + snapshot.particles.typeid[:] = [0, 1] + snapshot.particles.moment_inertia[:] = [(0.1, 0.1, 0.1)] * N + + return snapshot + + return make_snapshot + + +@pytest.mark.parametrize( + "friction_cls, friction_args, params, positions," + "velocities, angmom, orientations, force, torque", + friction_test_parameters, +) +def test_params( + friction_cls, + friction_args, + params, + positions, + velocities, + angmom, + orientations, + force, + torque, +): + rcut = 2.0 ** (1.0 / 6.0) * params["sigma"] + friction = friction_cls( + nlist=hoomd.md.nlist.Cell(buffer=0.4), default_r_cut=rcut, **friction_args + ) + + friction.params.default = params + for key in params: + assert friction.params[("A", "A")][key] == pytest.approx(params[key]) + assert friction.params[("A", "B")][key] == pytest.approx(params[key]) + assert friction.params[("B", "B")][key] == pytest.approx(params[key]) + + +# Test the only the deterministic forces and torques (zero temperature) +@pytest.mark.parametrize( + "friction_cls, friction_args, params, positions," + "velocities, angmom, orientations, force, torque", + friction_test_parameters, +) +def test_deterministic_forces_torques( + friction_snapshot_factory, + simulation_factory, + friction_cls, + friction_args, + params, + positions, + velocities, + angmom, + orientations, + force, + torque, +): + snapshot = friction_snapshot_factory( + position_i=positions[0], + position_j=positions[1], + velocity_i=velocities[0], + velocity_j=velocities[1], + angmom_i=angmom[0], + angmom_j=angmom[1], + orientation_i=orientations[0], + orientation_j=orientations[1], + ) + sim = simulation_factory(snapshot) + rcut = 2.0 ** (1.0 / 6.0) * params["sigma"] + friction = friction_cls( + nlist=hoomd.md.nlist.Cell(buffer=0.4), default_r_cut=rcut, **friction_args + ) + # Set temperature to zero + params["kT"] = 0 + friction.params.default = params + + sim.operations.integrator = hoomd.md.Integrator( + dt=0.001, forces=[friction], integrate_rotational_dof=True + ) + sim.run(0) + + sim_forces = friction.forces + sim_torques = friction.torques + if sim.device.communicator.rank == 0: + sim.orientations = snapshot.particles.orientation + + numpy.testing.assert_allclose(sim_forces[0], force, **TOLERANCES) + + numpy.testing.assert_allclose( + sim_forces[1], [-force[0], -force[1], -force[2]], **TOLERANCES + ) + + numpy.testing.assert_allclose(sim_torques[0], torque, **TOLERANCES) + + numpy.testing.assert_allclose(sim_torques[1], torque, **TOLERANCES) + + +# Test the frictional model for some timesteps (nonzero temperature) +@pytest.mark.parametrize( + "friction_cls, friction_args, params, positions," + "velocities, angmom, orientations, force, torque", + friction_test_parameters, +) +def test_forces_torques( + friction_snapshot_factory, + simulation_factory, + friction_cls, + friction_args, + params, + positions, + velocities, + angmom, + orientations, + force, + torque, +): + snapshot = friction_snapshot_factory( + position_i=positions[0], + position_j=positions[1], + velocity_i=velocities[0], + velocity_j=velocities[1], + angmom_i=angmom[0], + angmom_j=angmom[1], + orientation_i=orientations[0], + orientation_j=orientations[1], + ) + sim = simulation_factory(snapshot) + rcut = 2.0 ** (1.0 / 6.0) * params["sigma"] + friction = friction_cls( + nlist=hoomd.md.nlist.Cell(buffer=0.4), default_r_cut=rcut, **friction_args + ) + friction.params.default = params + sim.operations.integrator = hoomd.md.Integrator( + dt=0.001, forces=[friction], integrate_rotational_dof=True + ) + # Run some timesteps + sim.run(5) + sim_forces = friction.forces + sim_torques = friction.torques + if sim.device.communicator.rank == 0: + # check nonzero force and torques + for x in sim_forces[0]: + assert x != 0 + for x in sim_forces[1]: + assert x != 0 + for x in sim_torques[0]: + assert x != 0 + for x in sim_torques[1]: + assert x != 0 diff --git a/sphinx-doc/credits.rst b/sphinx-doc/credits.rst index 0af76c6a6c..e671996fce 100644 --- a/sphinx-doc/credits.rst +++ b/sphinx-doc/credits.rst @@ -134,4 +134,5 @@ The following people have contributed to HOOMD-blue: * Nathan Barrett, Pritzker School of Molecular Engineering * Domagoj Fijan, University of Michigan * Cristina Butu, Columbia University +* Kay Hofmann, Johannes Gutenberg University Mainz * Shuo-Lin Weng, Texas A&M University diff --git a/sphinx-doc/hoomd/md/module-pair.rst b/sphinx-doc/hoomd/md/module-pair.rst index d83e3f62bb..44e2ea81b2 100644 --- a/sphinx-doc/hoomd/md/module-pair.rst +++ b/sphinx-doc/hoomd/md/module-pair.rst @@ -11,6 +11,7 @@ pair :maxdepth: 1 pair/module-aniso + pair/module-friction .. rubric:: Classes diff --git a/sphinx-doc/hoomd/md/pair/friction/frictionalpair.rst b/sphinx-doc/hoomd/md/pair/friction/frictionalpair.rst new file mode 100644 index 0000000000..6b2787892f --- /dev/null +++ b/sphinx-doc/hoomd/md/pair/friction/frictionalpair.rst @@ -0,0 +1,8 @@ +FrictionalPair +============== + +.. py:currentmodule:: hoomd.md.pair.friction + +.. autoclass:: FrictionalPair + :members: + :show-inheritance: diff --git a/sphinx-doc/hoomd/md/pair/friction/frictionljcoulomb.rst b/sphinx-doc/hoomd/md/pair/friction/frictionljcoulomb.rst new file mode 100644 index 0000000000..14c3fb311e --- /dev/null +++ b/sphinx-doc/hoomd/md/pair/friction/frictionljcoulomb.rst @@ -0,0 +1,8 @@ +FrictionLJCoulomb +================= + +.. py:currentmodule:: hoomd.md.pair.friction + +.. autoclass:: FrictionLJCoulomb + :members: + :show-inheritance: diff --git a/sphinx-doc/hoomd/md/pair/friction/frictionljcoulombnewton.rst b/sphinx-doc/hoomd/md/pair/friction/frictionljcoulombnewton.rst new file mode 100644 index 0000000000..f46fa1eec4 --- /dev/null +++ b/sphinx-doc/hoomd/md/pair/friction/frictionljcoulombnewton.rst @@ -0,0 +1,8 @@ +FrictionLJCoulombNewton +======================= + +.. py:currentmodule:: hoomd.md.pair.friction + +.. autoclass:: FrictionLJCoulombNewton + :members: + :show-inheritance: diff --git a/sphinx-doc/hoomd/md/pair/friction/frictionljlinear.rst b/sphinx-doc/hoomd/md/pair/friction/frictionljlinear.rst new file mode 100644 index 0000000000..ae9cc270c2 --- /dev/null +++ b/sphinx-doc/hoomd/md/pair/friction/frictionljlinear.rst @@ -0,0 +1,8 @@ +FrictionLJLinear +================ + +.. py:currentmodule:: hoomd.md.pair.friction + +.. autoclass:: FrictionLJLinear + :members: + :show-inheritance: diff --git a/sphinx-doc/hoomd/md/pair/module-friction.rst b/sphinx-doc/hoomd/md/pair/module-friction.rst new file mode 100644 index 0000000000..a375b23dcc --- /dev/null +++ b/sphinx-doc/hoomd/md/pair/module-friction.rst @@ -0,0 +1,16 @@ +friction +======== + +.. automodule:: hoomd.md.pair.friction + :members: + :exclude-members: FrictionLJCoulomb,FrictionLJCoulombNewton,FrictionLJLinear,FrictionalPair + +.. rubric:: Classes + +.. toctree:: + :maxdepth: 1 + + friction/frictionljcoulomb + friction/frictionljcoulombnewton + friction/frictionljlinear + friction/frictionalpair