|
| 1 | +// Copyright (c) 2009-2025 The Regents of the University of Michigan. |
| 2 | +// Part of HOOMD-blue, released under the BSD 3-Clause License. |
| 3 | + |
| 4 | +#ifndef __PAIR_EVALUATOR_ZETTERLING_H__ |
| 5 | +#define __PAIR_EVALUATOR_ZETTERLING_H__ |
| 6 | + |
| 7 | +#ifndef __HIPCC__ |
| 8 | +#include <string> |
| 9 | +#endif |
| 10 | + |
| 11 | +#include "hoomd/HOOMDMath.h" |
| 12 | + |
| 13 | +/*! \file EvaluatorPairZetterling.h |
| 14 | + \brief Defines the pair evaluator class for Zetterling potential |
| 15 | +*/ |
| 16 | + |
| 17 | +// need to declare these class methods with __device__ qualifiers when building |
| 18 | +// in nvcc DEVICE is __host__ __device__ when included in nvcc and blank when |
| 19 | +// included into the host compiler |
| 20 | +#ifdef __HIPCC__ |
| 21 | +#define DEVICE __device__ |
| 22 | +#define HOSTDEVICE __host__ __device__ |
| 23 | +#else |
| 24 | +#define DEVICE |
| 25 | +#define HOSTDEVICE |
| 26 | +#endif |
| 27 | + |
| 28 | +namespace hoomd |
| 29 | + { |
| 30 | +namespace md |
| 31 | + { |
| 32 | +//! Class for evaluating the Zetterling pair potential |
| 33 | +/*! <b>General Overview</b> |
| 34 | +
|
| 35 | + <b>OPP specifics</b> |
| 36 | +
|
| 37 | + EvaluatorPairZetterling evaluates the function: |
| 38 | + \f{equation*} |
| 39 | + V_{\mathrm{Zetterling}}(r) = |
| 40 | + A \frac{\exp{(\alpha r)\cos{(2 k_F r)}}}{r^3} |
| 41 | + + B \left( \frac{\sigma}{r} \right)^n) |
| 42 | + \f} |
| 43 | +
|
| 44 | +*/ |
| 45 | +class EvaluatorPairZetterling |
| 46 | + { |
| 47 | + public: |
| 48 | + //! Define the parameter type used by this pair potential evaluator |
| 49 | + struct param_type |
| 50 | + { |
| 51 | + Scalar A; |
| 52 | + Scalar alpha; |
| 53 | + Scalar kf; |
| 54 | + Scalar B; |
| 55 | + Scalar sigma; |
| 56 | + Scalar n; |
| 57 | + |
| 58 | + DEVICE void load_shared(char*& ptr, unsigned int& available_bytes) { } |
| 59 | + |
| 60 | + HOSTDEVICE void allocate_shared(char*& ptr, unsigned int& available_bytes) const { } |
| 61 | + |
| 62 | +#ifdef ENABLE_HIP |
| 63 | + //! Set CUDA memory hints |
| 64 | + void set_memory_hint() const |
| 65 | + { |
| 66 | + // default implementation does nothing |
| 67 | + } |
| 68 | +#endif |
| 69 | + |
| 70 | +#ifndef __HIPCC__ |
| 71 | + param_type() : A(0), alpha(0), kf(0), B(0), sigma(0), n(0) { } |
| 72 | + |
| 73 | + param_type(pybind11::dict v, bool managed = false) |
| 74 | + { |
| 75 | + A = v["A"].cast<Scalar>(); |
| 76 | + alpha = v["alpha"].cast<Scalar>(); |
| 77 | + kf = v["kf"].cast<Scalar>(); |
| 78 | + B = v["B"].cast<Scalar>(); |
| 79 | + sigma = v["sigma"].cast<Scalar>(); |
| 80 | + n = v["n"].cast<Scalar>(); |
| 81 | + } |
| 82 | + |
| 83 | + pybind11::dict asDict() |
| 84 | + { |
| 85 | + pybind11::dict v; |
| 86 | + v["A"] = A; |
| 87 | + v["alpha"] = alpha; |
| 88 | + v["kf"] = kf; |
| 89 | + v["B"] = B; |
| 90 | + v["sigma"] = sigma; |
| 91 | + v["n"] = n; |
| 92 | + return v; |
| 93 | + } |
| 94 | +#endif |
| 95 | + }; |
| 96 | + |
| 97 | + //! Constructs the pair potential evaluator |
| 98 | + /*! \param _rsq Squared distance between the particles |
| 99 | + \param _rcutsq Squared distance at which the potential goes to 0 |
| 100 | + \param _params Per type pair parameters of this potential |
| 101 | + */ |
| 102 | + DEVICE EvaluatorPairZetterling(Scalar _rsq, Scalar _rcutsq, const param_type& _params) |
| 103 | + : rsq(_rsq), rcutsq(_rcutsq), params(_params) |
| 104 | + { |
| 105 | + } |
| 106 | + |
| 107 | + //! Zetterling doesn't use charge |
| 108 | + DEVICE static bool needsCharge() |
| 109 | + { |
| 110 | + return false; |
| 111 | + } |
| 112 | + |
| 113 | + //! Accept the optional charge values. |
| 114 | + /*! \param qi Charge of particle i |
| 115 | + \param qj Charge of particle j |
| 116 | + */ |
| 117 | + DEVICE void setCharge(Scalar qi, Scalar qj) { } |
| 118 | + |
| 119 | + //! Evaluate the force and energy |
| 120 | + /*! \param force_divr Output parameter to write the computed force |
| 121 | + * divided by r. |
| 122 | + * \param pair_eng Output parameter to write the computed pair energy |
| 123 | + * \param energy_shift If true, the potential must be shifted so that |
| 124 | + * V(r) is continuous at the cutoff |
| 125 | +
|
| 126 | + * \return True if they are evaluated or false if they are not because |
| 127 | + * we are beyond the cutoff |
| 128 | + */ |
| 129 | + DEVICE bool evalForceAndEnergy(Scalar& force_divr, Scalar& pair_eng, bool energy_shift) |
| 130 | + { |
| 131 | + if (rsq < rcutsq) |
| 132 | + { |
| 133 | + // Get quantities need for both energy and force calculation |
| 134 | + Scalar r(fast::sqrt(rsq)); |
| 135 | + Scalar eval_sin, eval_cos; |
| 136 | + fast::sincos(Scalar(2.0) * params.kf * r, eval_sin, eval_cos); |
| 137 | + Scalar screening(fast::exp(params.alpha * r)); |
| 138 | + Scalar inv_r(Scalar(1.0) / r); |
| 139 | + Scalar inv_r_2(Scalar(1.0) / rsq); |
| 140 | + Scalar inv_r_3(inv_r_2 * inv_r); |
| 141 | + Scalar inv_r_5(inv_r_3 * inv_r_2); |
| 142 | + Scalar power_sigma_over_r(fast::pow(params.sigma * inv_r, params.n)); |
| 143 | + |
| 144 | + // Compute energy |
| 145 | + Scalar term1(params.A * screening * eval_cos * inv_r_3); |
| 146 | + Scalar term2(params.B * power_sigma_over_r); |
| 147 | + pair_eng = term1 + term2; |
| 148 | + |
| 149 | + // Compute force |
| 150 | + Scalar deriv_term1( |
| 151 | + -params.A * screening |
| 152 | + * ((params.alpha * r - Scalar(3.0)) * eval_cos - 2 * params.kf * r * eval_sin) |
| 153 | + * inv_r_5); |
| 154 | + Scalar deriv_term2(params.B * params.n * power_sigma_over_r * inv_r_2); |
| 155 | + force_divr = deriv_term1 + deriv_term2; |
| 156 | + |
| 157 | + if (energy_shift) |
| 158 | + { |
| 159 | + Scalar r_cut(fast::sqrt(rcutsq)); |
| 160 | + Scalar screening_r_cut(fast::exp(params.alpha * r_cut)); |
| 161 | + Scalar inv_rcut(Scalar(1.0) / r_cut); |
| 162 | + Scalar inv_rcut_3(inv_rcut * inv_rcut * inv_rcut); |
| 163 | + Scalar term1_rcut(params.A * screening_r_cut |
| 164 | + * fast::cos(Scalar(2.0) * params.kf * r_cut) * inv_rcut_3); |
| 165 | + Scalar term2_rcut(params.B * fast::pow(params.sigma * inv_rcut, params.n)); |
| 166 | + pair_eng -= term1_rcut + term2_rcut; |
| 167 | + } |
| 168 | + |
| 169 | + return true; |
| 170 | + } |
| 171 | + else |
| 172 | + { |
| 173 | + return false; |
| 174 | + } |
| 175 | + } |
| 176 | + |
| 177 | + DEVICE Scalar evalPressureLRCIntegral() |
| 178 | + { |
| 179 | + return 0; |
| 180 | + } |
| 181 | + |
| 182 | + DEVICE Scalar evalEnergyLRCIntegral() |
| 183 | + { |
| 184 | + return 0; |
| 185 | + } |
| 186 | + |
| 187 | +#ifndef __HIPCC__ |
| 188 | + //! Get the name of this potential |
| 189 | + /*! \returns The potential name. |
| 190 | + */ |
| 191 | + static std::string getName() |
| 192 | + { |
| 193 | + return std::string("zetterling"); |
| 194 | + } |
| 195 | + |
| 196 | + std::string getShapeSpec() const |
| 197 | + { |
| 198 | + throw std::runtime_error("Shape definition not supported for this pair potential."); |
| 199 | + } |
| 200 | +#endif |
| 201 | + |
| 202 | + protected: |
| 203 | + Scalar rsq; /// Stored rsq from the constructor |
| 204 | + Scalar rcutsq; /// Stored rcutsq from the constructor |
| 205 | + param_type params; /// Stored pair parameters for a given type pair |
| 206 | + }; |
| 207 | + |
| 208 | + } // end namespace md |
| 209 | + } // end namespace hoomd |
| 210 | + |
| 211 | +#endif // __PAIR_EVALUATOR_ZETTERLING_H__ |
0 commit comments