diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index c4f539339ce..5134055a7ad 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -84,7 +84,7 @@ parenthesis: g = GPU, i = INTEL, k = KOKKOS, o = OPENMP, t = OPT. * :doc:`coul/tt ` * :doc:`coul/wolf (ko) ` * :doc:`coul/wolf/cs (o) ` - * :doc:`dispersion/d3 ` + * :doc:`dispersion/d3 (ko) ` * :doc:`dpd (giko) ` * :doc:`dpd/coul/slater/long (g) ` * :doc:`dpd/ext (ko) ` diff --git a/doc/src/Developer_utils.rst b/doc/src/Developer_utils.rst index 0ce2bd31f62..cff1385daa3 100644 --- a/doc/src/Developer_utils.rst +++ b/doc/src/Developer_utils.rst @@ -332,6 +332,9 @@ mathematical functions for a variety of applications. .. doxygenfunction:: powint :project: progguide +.. doxygenfunction:: powauto + :project: progguide + .. doxygenfunction:: powsinxx :project: progguide diff --git a/doc/src/pair_dispersion_d3.rst b/doc/src/pair_dispersion_d3.rst index eec12157f26..1e8aa449e2f 100644 --- a/doc/src/pair_dispersion_d3.rst +++ b/doc/src/pair_dispersion_d3.rst @@ -1,8 +1,12 @@ .. index:: pair_style dispersion/d3 +.. index:: pair_style dispersion/d3/kk +.. index:: pair_style dispersion/d3/omp pair_style dispersion/d3 command ================================ +Accelerator Variants: *dispersion/d3/kk*, *dispersion/d3/omp* + Syntax """""" @@ -105,6 +109,11 @@ provide are the chemical symbols of the atoms. The number of chemical symbols given must be equal to the number of atom types used and must match their ordering as atom types. +---------- + +.. include:: accel_styles.rst + +---------- Mixing, shift, table, tail correction, restart, rRESPA info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/src/EXTRA-PAIR/pair_dispersion_d3.cpp b/src/EXTRA-PAIR/pair_dispersion_d3.cpp index 1c86bb074f8..911b94799c1 100644 --- a/src/EXTRA-PAIR/pair_dispersion_d3.cpp +++ b/src/EXTRA-PAIR/pair_dispersion_d3.cpp @@ -27,39 +27,26 @@ #include "error.h" #include "force.h" #include "info.h" +#include "math_special.h" #include "memory.h" #include "neigh_list.h" #include "neighbor.h" #include "update.h" #include +#include #include #include -#include #include #include using namespace LAMMPS_NS; - -// global ad hoc parameters -static constexpr double K1 = 16.0; -static constexpr double K3 = -4.0; - -/* reasonable choices for k3 are between 3 and 5 : - this gives smoth curves with maxima around the integer values - k3=3 give for CN=0 a slightly smaller value than computed - for the free atom. This also yields to larger CN for atoms - in larger molecules but with the same chemical environment - which is physically not right. - values >5 might lead to bumps in the potential. -*/ +using namespace LAMMPS_NS::DispersionD3; static constexpr int NUM_ELEMENTS = 94; // maximum element number static constexpr int N_PARS_COLS = 5; // number of columns in C6 table static constexpr int N_PARS_ROWS = 32385; // number of rows in C6 table - -static constexpr double autoang = 0.52917725; // atomic units (Bohr) to Angstrom -static constexpr double autoev = 27.21140795; // atomic units (Hartree) to eV +static constexpr int N_C6AB_GRID = 5; // reference CN grid size of the C6 table #include "d3_parameters.h" @@ -82,6 +69,7 @@ PairDispersionD3::PairDispersionD3(LAMMPS *lmp) : dampingCode = 0; s6 = s8 = s18 = rs6 = rs8 = rs18 = a1 = a2 = alpha = alpha6 = alpha8 = 0.0; + max_mxci = 0; } /* ---------------------------------------------------------------------- @@ -90,6 +78,8 @@ PairDispersionD3::PairDispersionD3(LAMMPS *lmp) : PairDispersionD3::~PairDispersionD3() { + if (copymode) return; + if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); @@ -126,7 +116,7 @@ void PairDispersionD3::allocate() memory->create(rcov, n + 1, "pair:rcov"); memory->create(r0ab, n + 1, n + 1, "pair:r0ab"); - memory->create(c6ab, n + 1, n + 1, 5, 5, 3, "pair:c6ab"); + memory->create(c6ab, n + 1, n + 1, N_C6AB_GRID, N_C6AB_GRID, 3, "pair:c6ab"); } /* ---------------------------------------------------------------------- @@ -250,11 +240,15 @@ void PairDispersionD3::read_c6ab(int *atomic_numbers, int ntypes) for (int i = 0; i < N_PARS_ROWS; i++) { const double ref_c6 = c6ab_table[i][0]; - int atom_number_1 = (int)std::round(c6ab_table[i][1]); - int atom_number_2 = (int)std::round(c6ab_table[i][2]); + int atom_number_1 = (int) std::round(c6ab_table[i][1]); + int atom_number_2 = (int) std::round(c6ab_table[i][2]); set_limit_in_pars_array(atom_number_1, atom_number_2, grid_i, grid_j); + if ((grid_i >= N_C6AB_GRID) || (grid_j >= N_C6AB_GRID)) + error->all(FLERR, Error::NOLASTLINE, + "Reference C6 grid index out of range in the compiled in D3 parameters"); + std::vector idx_atoms_1 = is_int_in_array(atomic_numbers, ntypes, atom_number_1); if (idx_atoms_1.empty()) continue; @@ -293,7 +287,7 @@ void PairDispersionD3::coeff(int narg, char **arg) if (!allocated) allocate(); std::string element; - int *atomic_numbers = (int *) malloc(sizeof(int) * ntypes); + std::vector atomic_numbers(ntypes); for (int i = 0; i < ntypes; i++) { element = arg[i + 2]; atomic_numbers[i] = find_atomic_number(element); @@ -317,12 +311,14 @@ void PairDispersionD3::coeff(int narg, char **arg) } // set r0ab - read_r0ab(atomic_numbers, ntypes); + read_r0ab(atomic_numbers.data(), ntypes); - // read c6ab - read_c6ab(atomic_numbers, ntypes); + // read c6ab and record the largest reference CN grid index that is in use. + // The accelerated variants size their coefficient tables from it. - free(atomic_numbers); + read_c6ab(atomic_numbers.data(), ntypes); + max_mxci = 0; + for (int i = 1; i <= ntypes; i++) max_mxci = std::max(max_mxci, mxci[i]); } /* ---------------------------------------------------------------------- @@ -377,7 +373,7 @@ void PairDispersionD3::calc_coordination_number() if (rsq > cn_thr) continue; double rr = sqrt(rsq); - double rcov_ij = (rcov[itype] + rcov[jtype]) * autoang; + double rcov_ij = (rcov[itype] + rcov[jtype]) * AUTOANG; double cn_ij = 1.0 / (1.0 + exp(-K1 * ((rcov_ij / rr) - 1.0))); // update coordination number @@ -396,10 +392,8 @@ void PairDispersionD3::calc_coordination_number() Get derivative of C6 ------------------------------------------------------------------------- */ -double *PairDispersionD3::get_dC6(int iat, int jat, double cni, double cnj) +void PairDispersionD3::get_dC6(int iat, int jat, double cni, double cnj, double *c6_res) { - - static double c6_res[3] = {}; double c6_ref, cni_ref, cnj_ref; double c6mem, r_save, r; double expterm, term; @@ -417,7 +411,7 @@ double *PairDispersionD3::get_dC6(int iat, int jat, double cni, double cnj) for (int cj = 0; cj <= mxci[jat]; cj++) { c6_ref = c6ab[iat][jat][ci][cj][0]; - c6_ref *= autoev * pow(autoang, 6); + c6_ref *= AUTOEV * AUTOANG6; if (c6_ref > 0) { cni_ref = c6ab[iat][jat][ci][cj][1]; @@ -457,7 +451,6 @@ double *PairDispersionD3::get_dC6(int iat, int jat, double cni, double cnj) c6_res[1] = 0; c6_res[2] = 0; } - return c6_res; } /* ---------------------------------------------------------------------- @@ -510,16 +503,17 @@ void PairDispersionD3::compute(int eflag, int vflag) if (rsq < cutsq[type[i]][type[j]]) { - double r = sqrt(rsq); double r2inv = 1.0 / rsq; double r6inv = r2inv * r2inv * r2inv; double r8inv = r2inv * r2inv * r2inv * r2inv; double r10inv = r2inv * r2inv * r2inv * r2inv * r2inv; - double *c6_res = get_dC6(type[i], type[j], cn[i], cn[j]); + // get_dC6 writes {C6, dC6/dCN_i, dC6/dCN_j} + double c6_res[3] = {}; + get_dC6(type[i], type[j], cn[i], cn[j], c6_res); double C6 = c6_res[0]; - double C8 = 3.0 * C6 * r2r4[type[i]] * r2r4[type[j]] * autoang * autoang; + double C8 = 3.0 * C6 * r2r4[type[i]] * r2r4[type[j]] * AUTOANG * AUTOANG; double alpha6 = alpha; double alpha8 = alpha + 2; @@ -530,13 +524,18 @@ void PairDispersionD3::compute(int eflag, int vflag) switch (dampingCode) { - case 1: { // original + // Written to avoid using sqrt and pow() + case 1: /* Original damping */ + { + double ip6 = rs6 * r0ab[type[i]][type[j]]; + double ip8 = rs8 * r0ab[type[i]][type[j]]; - double r0 = r / r0ab[type[i]][type[j]]; + double half_alpha6 = 0.5 * alpha6; + double half_alpha8 = 0.5 * alpha8; - t6 = pow(rs6 / r0, alpha6); + t6 = MathSpecial::powauto(ip6, alpha6) * MathSpecial::powauto(rsq, -half_alpha6); damp6 = 1.0 / (1.0 + 6.0 * t6); - t8 = pow(rs8 / r0, alpha8); + t8 = MathSpecial::powauto(ip8, alpha8) * MathSpecial::powauto(rsq, -half_alpha8); damp8 = 1.0 / (1.0 + 6.0 * t8); e6 = C6 * damp6 * r6inv; @@ -551,14 +550,15 @@ void PairDispersionD3::compute(int eflag, int vflag) fpair = fpair1 + fpair2; fpair *= factor_lj; } break; - + // Written to avoid pow case 2: { // zerom + double r = sqrt(rsq); double r0 = r0ab[type[i]][type[j]]; - t6 = pow((r / (rs6 * r0)) + rs8 * r0, -alpha6); + t6 = MathSpecial::powauto((r / (rs6 * r0)) + rs8 * r0, -alpha6); damp6 = 1.0 / (1.0 + 6.0 * t6); - t8 = pow((r / r0) + rs8 * r0, -alpha8); + t8 = MathSpecial::powauto((r / r0) + rs8 * r0, -alpha8); damp8 = 1.0 / (1.0 + 6.0 * t8); e6 = C6 * damp6 * r6inv; @@ -578,7 +578,8 @@ void PairDispersionD3::compute(int eflag, int vflag) fpair *= factor_lj; } break; - case 3: { // bj + case 3: // bj + case 4: { // bjm, same functional form as bj, different parameters double r0 = sqrt(C8 / C6); @@ -586,29 +587,12 @@ void PairDispersionD3::compute(int eflag, int vflag) double r6 = rsq * rsq * rsq; double r8 = rsq * rsq * rsq * rsq; - t6 = r6 + pow((a1 * r0 + a2), 6); - t8 = r8 + pow((a1 * r0 + a2), 8); - - e6 = C6 / t6; - e8 = C8 / t8; - - tmp6 = 6.0 * s6 * C6 * r4 / (t6 * t6); - tmp8 = 8.0 * s8 * C8 * r6 / (t8 * t8); - - fpair = -(tmp6 + tmp8); - fpair *= factor_lj; - } break; - - case 4: { // bjm + double d = a1 * r0 + a2; + double d2 = d * d; + double d4 = d2 * d2; - double r0 = sqrt(C8 / C6); - - double r4 = rsq * rsq; - double r6 = rsq * rsq * rsq; - double r8 = rsq * rsq * rsq * rsq; - - t6 = r6 + pow((a1 * r0 + a2), 6); - t8 = r8 + pow((a1 * r0 + a2), 8); + t6 = r6 + MathSpecial::cube(d2); + t8 = r8 + MathSpecial::square(d4); e6 = C6 / t6; e8 = C8 / t8; @@ -685,7 +669,7 @@ void PairDispersionD3::compute(int eflag, int vflag) // here we calculate dcn = dCNi/dr = dCNj/dr if (rsq < cn_thr) { - double rcovij = (rcov[type[i]] + rcov[type[j]]) * autoang; + double rcovij = (rcov[type[i]] + rcov[type[j]]) * AUTOANG; double expterm = exp(-K1 * (rcovij / r - 1.0)); dcn = -K1 * rcovij * expterm / (rsq * (expterm + 1.0) * (expterm + 1.0)); } else { @@ -1039,7 +1023,7 @@ void PairDispersionD3::set_funcpar(std::string &functional_name) break; } - rs8 = rs8 / autoang; + rs8 = rs8 / AUTOANG; } break; case 3: { // bj @@ -1356,7 +1340,7 @@ void PairDispersionD3::set_funcpar(std::string &functional_name) break; } - a2 = a2 * autoang; + a2 = a2 * AUTOANG; } break; case 4: { // bjm @@ -1417,8 +1401,7 @@ void PairDispersionD3::set_funcpar(std::string &functional_name) break; } - a2 = a2 * autoang; - + a2 = a2 * AUTOANG; } break; default: // this should not happen with the error check in the init_style function @@ -1447,6 +1430,14 @@ void PairDispersionD3::init_style() { if (atom->tag_enable == 0) error->all(FLERR, "Pair style D3 requires atom IDs"); + // settings() only accepts known damping keywords, so this can not trigger. + // It is checked once here so that the threaded and device variants of the + // compute loop below do not have to report an error from inside a parallel + // region or a device kernel. + + if ((dampingCode < 1) || (dampingCode > 4)) + error->all(FLERR, Error::NOLASTLINE, "Damping code {} unknown", dampingCode); + // need an half neighbor list neighbor->add_request(this); } @@ -1531,8 +1522,8 @@ double PairDispersionD3::memory_usage() { double bytes = Pair::memory_usage(); int n = atom->ntypes; - // c6ab[n+1][n+1][5][5][3] coefficient table - bytes += (double)(n+1)*(n+1)*5*5*3 * sizeof(double); + // c6ab[n+1][n+1][N_C6AB_GRID][N_C6AB_GRID][3] coefficient table + bytes += (double) (n + 1) * (n + 1) * N_C6AB_GRID * N_C6AB_GRID * 3 * sizeof(double); // per-atom coordination number and C6 derivative arrays bytes += (double) nmax * 2 * sizeof(double); // cn[nmax] + dc6[nmax] return bytes; diff --git a/src/EXTRA-PAIR/pair_dispersion_d3.h b/src/EXTRA-PAIR/pair_dispersion_d3.h index a9d450f3d57..e9298345935 100644 --- a/src/EXTRA-PAIR/pair_dispersion_d3.h +++ b/src/EXTRA-PAIR/pair_dispersion_d3.h @@ -22,8 +22,34 @@ PairStyle(dispersion/d3,PairDispersionD3); #include "pair.h" -namespace LAMMPS_NS { +// global ad hoc parameters of the D3 model, shared by the plain pair style and +// its accelerated variants. They live in a named namespace, not at global +// scope, because this is a style header and gets included all over the place. + +namespace LAMMPS_NS::DispersionD3 { + + static constexpr double K1 = 16.0; + + /* reasonable choices for k3 are between 3 and 5 : + this gives smooth curves with maxima around the integer values + k3=3 give for CN=0 a slightly smaller value than computed + for the free atom. This also yields to larger CN for atoms + in larger molecules but with the same chemical environment + which is physically not right. + values >5 might lead to bumps in the potential. + */ + static constexpr double K3 = -4.0; + + static constexpr double AUTOANG = 0.52917725; // atomic units (Bohr) to Angstrom + static constexpr double AUTOEV = 27.21140795; // atomic units (Hartree) to eV + + // conversion factor for the tabulated C6 reference values + static constexpr double AUTOANG3 = AUTOANG * AUTOANG * AUTOANG; + static constexpr double AUTOANG6 = AUTOANG3 * AUTOANG3; +} // namespace LAMMPS_NS::DispersionD3 + +namespace LAMMPS_NS { class PairDispersionD3 : public Pair { public: @@ -48,25 +74,26 @@ class PairDispersionD3 : public Pair { double rthr; // R^2 distance to cutoff for D3_calculation double cn_thr; // R^2 distance to cutoff for CN_calculation - int dampingCode; - double s6, s8, s18, rs6, rs8, rs18; // XC parameters - double a1, a2, alpha, alpha6, alpha8; + int dampingCode; // Which damping function to use + double s6, s8, s18, rs6, rs8, rs18; // XC parameters + double a1, a2, alpha, alpha6, alpha8; // XC parameters double *r2r4; // scale r4/r2 values of the atoms by sqrt(Z) double *rcov; // covalent radii int *mxci; // How large the grid for c6 interpolation double **r0ab; // cut-off radii for all element pairs double *****c6ab; // C6 for all element pairs + int max_mxci; // Maximum grid size of the C_i coefficient double *cn; // Coordination numbers double *dc6; // dC6i(iat) saves dE_dsp/dCN(iat) int communicationStage; // communication stage double memory_usage() override; - void allocate(); + virtual void allocate(); virtual void set_funcpar(std::string &); - void calc_coordination_number(); + virtual void calc_coordination_number(); int find_atomic_number(std::string &); std::vector is_int_in_array(int *, int, int); @@ -75,7 +102,9 @@ class PairDispersionD3 : public Pair { void set_limit_in_pars_array(int &, int &, int &, int &); void read_c6ab(int *, int); - double *get_dC6(int, int, double, double); + // writes {C6, dC6/dCN_i, dC6/dCN_j} to c6_res; must not use static storage, + // it is called concurrently from the threaded and device variants + void get_dC6(int, int, double, double, double *); }; } // namespace LAMMPS_NS #endif diff --git a/src/KOKKOS/Install.sh b/src/KOKKOS/Install.sh index 4e1429da7d7..4f717213087 100755 --- a/src/KOKKOS/Install.sh +++ b/src/KOKKOS/Install.sh @@ -394,6 +394,8 @@ action pair_coul_shield_kokkos.cpp pair_coul_shield.cpp action pair_coul_shield_kokkos.h pair_coul_shield.h action pair_cosine_squared_kokkos.cpp pair_cosine_squared.cpp action pair_cosine_squared_kokkos.h pair_cosine_squared.h +action pair_dispersion_d3_kokkos.cpp pair_dispersion_d3.cpp +action pair_dispersion_d3_kokkos.h pair_dispersion_d3.h action pair_dpd_ext_kokkos.cpp pair_dpd_ext.cpp action pair_dpd_ext_kokkos.h pair_dpd_ext.h action pair_dpd_ext_tstat_kokkos.cpp pair_dpd_ext_tstat.cpp diff --git a/src/KOKKOS/math_special_kokkos.h b/src/KOKKOS/math_special_kokkos.h index 5ecb23a891e..98885338d20 100644 --- a/src/KOKKOS/math_special_kokkos.h +++ b/src/KOKKOS/math_special_kokkos.h @@ -234,6 +234,29 @@ namespace LAMMPS_NS::MathSpecialKokkos { return (n > 0) ? yy : static_cast(1) / yy; } + /* Version of pow(x,n) for an exponent that is only known at run time + * + * Device capable equivalent of MathSpecial::powauto(). Dispatches to + * powint() when the exponent happens to have an integer value, and to + * Kokkos::pow() otherwise. Follows the powint() convention of returning + * 0 for x == 0, which differs from pow() for a negative exponent. + * + * \param x base + * \param n exponent + * \return value of x^n */ + + template +// NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + static T powauto(const T &x, const T &n) + { + if (n == static_cast(0)) return static_cast(1); + if (x == static_cast(0)) return static_cast(0); + const T rounded = Kokkos::round(n); + if (n == rounded) return powint(x, static_cast(rounded)); + return Kokkos::pow(x, n); + } + /* Fast inline version of (sin(x)/x)^n as used by PPPM kspace styles * * This is an optimized function to compute (sin(x)/x)^n as frequently used by PPPM. diff --git a/src/KOKKOS/pair_dispersion_d3_kokkos.cpp b/src/KOKKOS/pair_dispersion_d3_kokkos.cpp new file mode 100644 index 00000000000..33e8db392e9 --- /dev/null +++ b/src/KOKKOS/pair_dispersion_d3_kokkos.cpp @@ -0,0 +1,749 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: + Yizhong R. Hu + Marc L. Descoteaux + Ulrik Unneberg + William C. Witt + Affiliation: Harvard University +------------------------------------------------------------------------- */ + +#include "pair_dispersion_d3_kokkos.h" + +#include "atom_kokkos.h" +#include "atom_masks.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "kokkos.h" +#include "memory_kokkos.h" +#include "neigh_list_kokkos.h" +#include "neigh_request.h" +#include "neighbor.h" +#include "pair_kokkos.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +template +PairDispersionD3Kokkos::PairDispersionD3Kokkos(LAMMPS *lmp) : PairDispersionD3(lmp) +{ + kokkosable = 1; + + // cn and dc6 are reduced with reverse_comm() and are device resident, so + // let CommKokkos use the device pack/unpack methods below + reverse_comm_device = 1; + + atomKK = (AtomKokkos *) atom; + execution_space = ExecutionSpaceFromDevice::space; + datamask_read = X_MASK | F_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK; + datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK; +} + +/* ---------------------------------------------------------------------- */ + +template +PairDispersionD3Kokkos::~PairDispersionD3Kokkos() +{ + if (copymode) return; + + if (allocated) { + memoryKK->destroy_kokkos(k_eatom,eatom); + memoryKK->destroy_kokkos(k_vatom,vatom); + memoryKK->destroy_kokkos(k_cutsq,cutsq); + cutsq = nullptr; // prevent base destructor from double-freeing + } +} + +/* ---------------------------------------------------------------------- + Calculate coordination number of atoms +------------------------------------------------------------------------- */ + +template +void PairDispersionD3Kokkos::calc_coordination_number() +{ + if (atom->nmax > nmax) { + nmax = atom->nmax; + k_cn = DAT::tdual_kkfloat_1d("pair:cn",nmax); + k_dc6 = DAT::tdual_kkfloat_1d("pair:dc6",nmax); + d_cn = k_cn.template view(); + d_dc6 = k_dc6.template view(); + h_cn = k_cn.view_host(); + h_dc6 = k_dc6.view_host(); + } + + // cn/dc6 are recomputed from scratch; discard prior sync state + k_cn.clear_sync_state(); + k_dc6.clear_sync_state(); + k_cn.template modify(); + k_dc6.template modify(); + + // zero out coordination number and dC6 + + Kokkos::parallel_for( + Kokkos::RangePolicy(0, newton_pair ? nall : nlocal), + PairDispersionD3InitializeFunctor{d_cn, d_dc6}); + + // calculate coordination number + + if (neighflag == FULL) { + dispatch_coordination_kernel(); + } else if (neighflag == HALFTHREAD) { + dispatch_coordination_kernel(); + } else if (neighflag == HALF) { + dispatch_coordination_kernel(); + } else { + error->all(FLERR, "Must use half or full neighbor list style with pair dispersion/d3/kk"); + } + + // communicate coordination number + communicationStage = 1; + if (newton_pair) comm->reverse_comm(this); + comm->forward_comm(this); +} + +/* ---------------------------------------------------------------------- */ + +template +template +void PairDispersionD3Kokkos::dispatch_coordination_kernel() +{ + if (newton_pair) { + PairDispersionD3CoordinationNumberKernel cnkernel( + x, type, d_rcov, d_cn, d_ilist, d_numneigh, d_neighbors, nlocal, cn_thr); + Kokkos::parallel_for(Kokkos::RangePolicy(0, inum), cnkernel); + cnkernel.contribute(); + } else { + PairDispersionD3CoordinationNumberKernel cnkernel( + x, type, d_rcov, d_cn, d_ilist, d_numneigh, d_neighbors, nlocal, cn_thr); + Kokkos::parallel_for(Kokkos::RangePolicy(0, inum), cnkernel); + cnkernel.contribute(); + } +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +template +double PairDispersionD3Kokkos::init_one(int i, int j) +{ + const double cut = PairDispersionD3::init_one(i, j); + // Since cutsq is written on host by base class we mark the host view modified + k_cutsq.modify_host(); + return cut; +} + +/* ---------------------------------------------------------------------- + init pair style +------------------------------------------------------------------------- */ + +template +void PairDispersionD3Kokkos::init_style() +{ + PairDispersionD3::init_style(); + + // adjust neighbor list request for KOKKOS + + neighflag = lmp->kokkos->neighflag; + + // a full neighbor list visits every pair twice, so the pairwise forces do + // not add up to the fdotr virial and it has to be tallied explicitly + + if (neighflag == FULL) no_virial_fdotr_compute = 1; + + auto request = neighbor->find_request(this); + request->set_kokkos_host(std::is_same_v && + !std::is_same_v); + request->set_kokkos_device(std::is_same_v); + if (neighflag == FULL) request->enable_full(); +} + +/* ---------------------------------------------------------------------- + allocate req. arrays +------------------------------------------------------------------------- */ + +template +void PairDispersionD3Kokkos::allocate() +{ + PairDispersionD3::allocate(); + + int n = atom->ntypes; + memory->destroy(cutsq); + memoryKK->create_kokkos(k_cutsq,cutsq,n+1,n+1,"pair:cutsq"); + d_cutsq = k_cutsq.template view(); + + k_r2r4 = DAT::tdual_kkfloat_1d("pair:r2r4", n+1); + k_rcov = DAT::tdual_kkfloat_1d("pair:rcov", n+1); + k_mxci = DAT::tdual_int_1d("pair:mxci", n+1); + k_r0ab = DAT::tdual_kkfloat_2d("pair:r0ab", n+1, n+1); + + d_r2r4 = k_r2r4.template view(); + d_rcov = k_rcov.template view(); + d_mxci = k_mxci.template view(); + d_r0ab = k_r0ab.template view(); + + // k_c6ab is created in coeff(), where the reference CN grid size is known +} + +/* ---------------------------------------------------------------------- + Coeff: read from pair_coeff + pair_coeff * * path_r0ab.csv path_c6ab.csv functional element1 element2 ... +------------------------------------------------------------------------- */ + +template +void PairDispersionD3Kokkos::coeff(int narg, char **arg) +{ + // the base class parses the arguments and fills the host side tables + PairDispersionD3::coeff(narg, arg); + + const int ntypes = atom->ntypes; + + // now that max_mxci is known, resize the C6 table to the grid actually used + + const int ngrid = max_mxci + 1; + if ((k_c6ab.extent_int(2) != ngrid) || (k_c6ab.extent_int(3) != ngrid)) { + k_c6ab = decltype(k_c6ab)("pair:c6ab", ntypes + 1, ntypes + 1, ngrid, ngrid, 3); + d_c6ab = k_c6ab.template view(); + } + + // copy coefficients to device + + auto h_r2r4 = k_r2r4.view_host(); + auto h_rcov = k_rcov.view_host(); + auto h_mxci = k_mxci.view_host(); + auto h_r0ab = k_r0ab.view_host(); + auto h_c6ab = k_c6ab.view_host(); + + for (int i = 1; i <= ntypes; i++) { + h_r2r4(i) = r2r4[i]; + h_rcov(i) = rcov[i]; + h_mxci(i) = mxci[i]; + } + + for (int i = 1; i <= ntypes; i++) { + for (int j = 1; j <= ntypes; j++) { + h_r0ab(i, j) = r0ab[i][j]; + for (int ci = 0; ci < ngrid; ci++) + for (int cj = 0; cj < ngrid; cj++) + for (int k = 0; k < 3; k++) h_c6ab(i, j, ci, cj, k) = c6ab[i][j][ci][cj][k]; + } + } + + k_r2r4.template modify(); + k_rcov.template modify(); + k_mxci.template modify(); + k_r0ab.template modify(); + k_c6ab.template modify(); + + k_r2r4.template sync(); + k_rcov.template sync(); + k_mxci.template sync(); + k_r0ab.template sync(); + k_c6ab.template sync(); +} + +/* ---------------------------------------------------------------------- + Compute : energy, force, and stress (Required) +------------------------------------------------------------------------- */ + +template +void PairDispersionD3Kokkos::compute(int eflag_in, int vflag_in) +{ + eflag = eflag_in; + vflag = vflag_in; + + ev_init(eflag,vflag,0); + + // reallocate per-atom arrays if necessary + + if (eflag_atom) { + memoryKK->destroy_kokkos(k_eatom,eatom); + memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom"); + d_eatom = k_eatom.view(); + } + if (vflag_atom) { + memoryKK->destroy_kokkos(k_vatom,vatom); + memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom"); + d_vatom = k_vatom.view(); + } + + atomKK->sync(execution_space,datamask_read); + k_cutsq.template sync(); + k_mxci.template sync(); + k_c6ab.template sync(); + if (eflag || vflag) + atomKK->modified(execution_space,datamask_modify); + else + atomKK->modified(execution_space,F_MASK); + + x = atomKK->k_x.view(); + f = atomKK->k_f.view(); + type = atomKK->k_type.view(); + nlocal = atom->nlocal; + nall = atom->nlocal + atom->nghost; + newton_pair = force->newton_pair; + special_lj[0] = force->special_lj[0]; + special_lj[1] = force->special_lj[1]; + special_lj[2] = force->special_lj[2]; + special_lj[3] = force->special_lj[3]; + + NeighListKokkos* k_list = static_cast*>(list); + d_numneigh = k_list->d_numneigh; + d_neighbors = k_list->d_neighbors; + d_ilist = k_list->d_ilist; + inum = list->inum; + + // must stay in sync with the DUP alias of the kernels below, or the wrong + // one of the dup_/ndup_ scatter view pair is handed to the functor + + if (neighflag == FULL) { + need_dup = std::is_same_v, Kokkos::Experimental::ScatterDuplicated>; + } else if (neighflag == HALFTHREAD) { + need_dup = std::is_same_v, Kokkos::Experimental::ScatterDuplicated>; + } else { + need_dup = std::is_same_v, Kokkos::Experimental::ScatterDuplicated>; + } + + EV_FLOAT ev_all = {}; + + // Calculate the coordination number for each atom + calc_coordination_number(); + + // Since communication impacts host view we refresh device view before the device kernels + k_cn.template sync(); + + // Could move these scatter view manipulations into a separate function + dup_eatom = {}; + dup_vatom = {}; + ndup_eatom = {}; + ndup_vatom = {}; + dup_dc6 = {}; + ndup_dc6 = {}; + if (need_dup) { + dup_f = Kokkos::Experimental::create_scatter_view(f); + } + ndup_f = Kokkos::Experimental::create_scatter_view(f); + if (need_dup) { + dup_dc6 = Kokkos::Experimental::create_scatter_view(d_dc6); + } + ndup_dc6 = Kokkos::Experimental::create_scatter_view(d_dc6); + + if (eflag_atom) { + if (need_dup) { + dup_eatom = Kokkos::Experimental::create_scatter_view(d_eatom); + } + ndup_eatom = Kokkos::Experimental::create_scatter_view(d_eatom); + } + if (vflag_atom) { + if (need_dup) { + dup_vatom = Kokkos::Experimental::create_scatter_view(d_vatom); + } + ndup_vatom = Kokkos::Experimental::create_scatter_view(d_vatom); + } + + // first: compute forces, energies, and dC6 + EV_FLOAT ev = {}; + + if (neighflag == HALF) { + dispatch_kernel_A(ev); + } else if (neighflag == HALFTHREAD) { + dispatch_kernel_A(ev); + } else if (neighflag == FULL) { + dispatch_kernel_A(ev); + } else { + error->all(FLERR, "Must use half or full neighbor list style with pair dispersion/d3/kk"); + } + + if (evflag) ev_all += ev; + + if (need_dup) + Kokkos::Experimental::contribute(d_dc6, dup_dc6); + else + Kokkos::Experimental::contribute(d_dc6, ndup_dc6); + k_dc6.template modify(); + + // communicate derivatives of C6 + communicationStage = 2; + if (newton_pair) comm->reverse_comm(this); + comm->forward_comm(this); + + // Communication updates host view so we must refresh device view before kernel B + k_dc6.template sync(); + + // second: compute and apply force contribution from dC6 + ev = {}; + + if (neighflag == HALF) { + dispatch_kernel_B(ev); + } else if (neighflag == HALFTHREAD) { + dispatch_kernel_B(ev); + } else if (neighflag == FULL) { + dispatch_kernel_B(ev); + } else { + error->all(FLERR, "Must use half or full neighbor list style with pair dispersion/d3/kk"); + } + + if (evflag) ev_all += ev; + + if (need_dup) + Kokkos::Experimental::contribute(f, dup_f); + else + Kokkos::Experimental::contribute(f, ndup_f); + + if (eflag_global) eng_vdwl += ev_all.evdwl; + if (vflag_global) { + virial[0] += ev_all.v[0]; + virial[1] += ev_all.v[1]; + virial[2] += ev_all.v[2]; + virial[3] += ev_all.v[3]; + virial[4] += ev_all.v[4]; + virial[5] += ev_all.v[5]; + } + + if (eflag_atom) { + if (need_dup) + Kokkos::Experimental::contribute(d_eatom, dup_eatom); + else + Kokkos::Experimental::contribute(d_eatom, ndup_eatom); + k_eatom.template modify(); + k_eatom.sync_host(); + } + + if (vflag_atom) { + if (need_dup) + Kokkos::Experimental::contribute(d_vatom, dup_vatom); + else + Kokkos::Experimental::contribute(d_vatom, ndup_vatom); + k_vatom.template modify(); + k_vatom.sync_host(); + } + + if (vflag_fdotr) pair_virial_fdotr_compute(this); + + // free duplicated memory + if (need_dup) { + dup_f = {}; + dup_eatom = {}; + dup_vatom = {}; + dup_dc6 = {}; + } +} + +/* ---------------------------------------------------------------------- + Templated kernel methods +------------------------------------------------------------------------- */ + +template +template +void PairDispersionD3Kokkos::launch_kernel_A(EV_FLOAT &ev) +{ + auto functor = PairDispersionD3KernelA( + x, type, + d_cutsq, d_cn, d_dc6, + d_r2r4, d_r0ab, d_c6ab, d_mxci, + d_numneigh, d_neighbors, d_ilist, + dup_f, ndup_f, + dup_eatom, ndup_eatom, + dup_vatom, ndup_vatom, + dup_dc6, ndup_dc6, + special_lj, + nlocal, + eflag, vflag_either, + eflag_global, eflag_atom, + vflag_global, vflag_atom, + dampingCode, + s6, s8, rs6, rs8, + a1, a2, alpha); + + if constexpr (EVFLAG) { + Kokkos::parallel_reduce( + Kokkos::RangePolicy(0, inum), + functor, + ev); + } else { + Kokkos::parallel_for( + Kokkos::RangePolicy(0, inum), + functor); + } +} + +/* ---------------------------------------------------------------------- */ + +template +template +void PairDispersionD3Kokkos::dispatch_kernel_A(EV_FLOAT &ev) +{ + if (newton_pair) { + if (evflag) launch_kernel_A(ev); + else launch_kernel_A(ev); + } else { + if (evflag) launch_kernel_A(ev); + else launch_kernel_A(ev); + } +} + +/* ---------------------------------------------------------------------- */ + +template +template +void PairDispersionD3Kokkos::launch_kernel_B(EV_FLOAT &ev) +{ + auto functor = PairDispersionD3KernelB( + x, type, + d_cutsq, d_dc6, d_rcov, + d_numneigh, d_neighbors, d_ilist, + dup_f, ndup_f, + dup_eatom, ndup_eatom, + dup_vatom, ndup_vatom, + special_lj, + nlocal, + eflag, vflag_either, + eflag_global, eflag_atom, + vflag_global, vflag_atom, + cn_thr); + + if constexpr (EVFLAG) { + Kokkos::parallel_reduce( + Kokkos::RangePolicy(0, inum), + functor, + ev); + } else { + Kokkos::parallel_for( + Kokkos::RangePolicy(0, inum), + functor); + } +} + +/* ---------------------------------------------------------------------- */ + +template +template +void PairDispersionD3Kokkos::dispatch_kernel_B(EV_FLOAT &ev) +{ + if (newton_pair) { + if (evflag) launch_kernel_B(ev); + else launch_kernel_B(ev); + } else { + if (evflag) launch_kernel_B(ev); + else launch_kernel_B(ev); + } +} + +/* ---------------------------------------------------------------------- + Communication section +------------------------------------------------------------------------- */ + +template +int PairDispersionD3Kokkos::pack_forward_comm_kokkos(int n, DAT::tdual_int_1d k_sendlist, + DAT::tdual_double_1d &buf, + int /*pbc_flag*/, int * /*pbc*/) +{ + d_sendlist = k_sendlist.view(); + v_buf = buf.view(); + if (communicationStage == 1) { + k_cn.template sync(); + Kokkos::parallel_for( + Kokkos::RangePolicy(0, n), + PairDispersionD3PackForwardCommFunctor{d_cn, d_sendlist, v_buf}); + } + if (communicationStage == 2) { + k_dc6.template sync(); + Kokkos::parallel_for( + Kokkos::RangePolicy(0, n), + PairDispersionD3PackForwardCommFunctor{d_dc6, d_sendlist, v_buf}); + } + return n; +} + +/* ---------------------------------------------------------------------- */ + +template +void PairDispersionD3Kokkos::unpack_forward_comm_kokkos(int n, int first, DAT::tdual_double_1d &buf) +{ + v_buf = buf.view(); + if (communicationStage == 1) { + k_cn.template sync(); + Kokkos::parallel_for( + Kokkos::RangePolicy(0, n), + PairDispersionD3UnpackForwardCommFunctor{d_cn, first, v_buf}); + k_cn.template modify(); + } + if (communicationStage == 2) { + k_dc6.template sync(); + Kokkos::parallel_for( + Kokkos::RangePolicy(0, n), + PairDispersionD3UnpackForwardCommFunctor{d_dc6, first, v_buf}); + k_dc6.template modify(); + } +} + +/* ---------------------------------------------------------------------- */ + +template +int PairDispersionD3Kokkos::pack_reverse_comm_kokkos(int n, int first, DAT::tdual_double_1d &buf) +{ + v_buf = buf.view(); + if (communicationStage == 1) { + k_cn.template sync(); + Kokkos::parallel_for( + Kokkos::RangePolicy(0, n), + PairDispersionD3PackReverseCommFunctor{d_cn, first, v_buf}); + } + if (communicationStage == 2) { + k_dc6.template sync(); + Kokkos::parallel_for( + Kokkos::RangePolicy(0, n), + PairDispersionD3PackReverseCommFunctor{d_dc6, first, v_buf}); + } + return n; +} + +/* ---------------------------------------------------------------------- */ + +template +void PairDispersionD3Kokkos::unpack_reverse_comm_kokkos(int n, DAT::tdual_int_1d k_recvlist, DAT::tdual_double_1d &buf) +{ + d_recvlist = k_recvlist.view(); + v_buf = buf.view(); + if (communicationStage == 1) { + k_cn.template sync(); + Kokkos::parallel_for( + Kokkos::RangePolicy(0, n), + PairDispersionD3UnpackReverseCommFunctor{d_cn, d_recvlist, v_buf}); + k_cn.template modify(); + } + if (communicationStage == 2) { + k_dc6.template sync(); + Kokkos::parallel_for( + Kokkos::RangePolicy(0, n), + PairDispersionD3UnpackReverseCommFunctor{d_dc6, d_recvlist, v_buf}); + k_dc6.template modify(); + } +} + +/* ---------------------------------------------------------------------- */ + +template +int PairDispersionD3Kokkos::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, + int * /*pbc*/) +{ + int i, j, m; + + m = 0; + if (communicationStage == 1) { + + k_cn.sync_host(); + + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = h_cn(j); + } + } + if (communicationStage == 2) { + + k_dc6.sync_host(); + + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = h_dc6(j); + } + } + + return m; +} + +/* ---------------------------------------------------------------------- */ + +template +void PairDispersionD3Kokkos::unpack_forward_comm(int n, int first, double *buf) +{ + int i, m, last; + + m = 0; + last = first + n; + if (communicationStage == 1) { + k_cn.sync_host(); + + for (i = first; i < last; i++) { h_cn(i) = buf[m++]; } + + k_cn.modify_host(); + } + if (communicationStage == 2) { + + k_dc6.sync_host(); + + for (i = first; i < last; i++) { h_dc6(i) = buf[m++]; } + + k_dc6.modify_host(); + } +} + +/* ---------------------------------------------------------------------- */ + +template +int PairDispersionD3Kokkos::pack_reverse_comm(int n, int first, double *buf) +{ + int i, m, last; + + m = 0; + last = first + n; + if (communicationStage == 1) { + k_cn.sync_host(); + for (i = first; i < last; i++) { buf[m++] = h_cn(i); } + } + if (communicationStage == 2) { + k_dc6.sync_host(); + for (i = first; i < last; i++) { buf[m++] = h_dc6(i); } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +template +void PairDispersionD3Kokkos::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i, j, m; + + m = 0; + if (communicationStage == 1) { + k_cn.sync_host(); + for (i = 0; i < n; i++) { + j = list[i]; + h_cn(j) += buf[m++]; + } + k_cn.modify_host(); + } + if (communicationStage == 2) { + k_dc6.sync_host(); + for (i = 0; i < n; i++) { + j = list[i]; + h_dc6(j) += buf[m++]; + } + k_dc6.modify_host(); + } +} + +/* ---------------------------------------------------------------------- */ + +namespace LAMMPS_NS { +template class PairDispersionD3Kokkos; +#ifdef LMP_KOKKOS_GPU +template class PairDispersionD3Kokkos; +#endif +} diff --git a/src/KOKKOS/pair_dispersion_d3_kokkos.h b/src/KOKKOS/pair_dispersion_d3_kokkos.h new file mode 100644 index 00000000000..5d7027939bd --- /dev/null +++ b/src/KOKKOS/pair_dispersion_d3_kokkos.h @@ -0,0 +1,1044 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(dispersion/d3/kk,PairDispersionD3Kokkos); +PairStyle(dispersion/d3/kk/device,PairDispersionD3Kokkos); +PairStyle(dispersion/d3/kk/host,PairDispersionD3Kokkos); +// clang-format on +#else + +// clang-format off +#ifndef LMP_PAIR_DISPERSION_D3_KOKKOS_H +#define LMP_PAIR_DISPERSION_D3_KOKKOS_H + +#include "kokkos_base.h" +#include "math_special_kokkos.h" +#include "pair_dispersion_d3.h" +#include "pair_kokkos.h" + +#include + +namespace LAMMPS_NS { + +using DispersionD3::AUTOANG; +using DispersionD3::AUTOANG6; +using DispersionD3::AUTOEV; +using DispersionD3::K1; +using DispersionD3::K3; + +/* ---------------------------------------------------------------------- */ +// Functor to initialize cn and dc6 arrays +/* ---------------------------------------------------------------------- */ + +template +struct PairDispersionD3InitializeFunctor { + typedef ArrayTypes AT; + + typename AT::t_kkfloat_1d d_cn; + typename AT::t_kkfloat_1d d_dc6; + + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void operator()(const int &i) const { + d_cn(i) = 0.0; + d_dc6(i) = 0.0; + } +}; + +template +struct PairDispersionD3PackForwardCommFunctor { + typedef ArrayTypes AT; + + typename AT::t_kkfloat_1d d_arr; + typename AT::t_int_1d d_sendlist; + typename AT::t_double_1d_um v_buf; + + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void operator()(const int &i) const { + int j = d_sendlist(i); + v_buf(i) = d_arr(j); + } +}; + +template +struct PairDispersionD3UnpackForwardCommFunctor { + typedef ArrayTypes AT; + + typename AT::t_kkfloat_1d d_arr; + int first; + typename AT::t_double_1d_um v_buf; + + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void operator()(const int &i) const { + d_arr(i + first) = v_buf(i); + } +}; + +template +struct PairDispersionD3PackReverseCommFunctor { + typedef ArrayTypes AT; + + typename AT::t_kkfloat_1d d_arr; + int first; + typename AT::t_double_1d_um v_buf; + + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void operator()(const int &i) const { + v_buf(i) = d_arr(i + first); + } +}; + +template +struct PairDispersionD3UnpackReverseCommFunctor { + typedef ArrayTypes AT; + + typename AT::t_kkfloat_1d d_arr; + typename AT::t_int_1d d_recvlist; + typename AT::t_double_1d_um v_buf; + + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void operator()(const int &i) const { + int j = d_recvlist(i); + d_arr(j) += v_buf(i); + } +}; + +template +struct PairDispersionD3CoordinationNumberKernel { + typedef ArrayTypes AT; + using DUP = NeedDup_v; + using ScatterAccess = std::conditional_t< + std::is_same_v, + Kokkos::Experimental::ScatterNonAtomic, + Kokkos::Experimental::ScatterAtomic>; + + typename AT::t_kkfloat_1d_3_lr_randomread x; + typename AT::t_int_1d_randomread type; + typename AT::t_kkfloat_1d d_rcov; + typename AT::t_kkfloat_1d d_cn; + typename AT::t_int_1d d_ilist; + typename AT::t_int_1d d_numneigh; + typename AT::t_neighbors_2d d_neighbors; + int nlocal; + KK_FLOAT cn_thr; + + KKScatterView::value, KKScatterSum, DUP> dup_cn; + + PairDispersionD3CoordinationNumberKernel( + const typename AT::t_kkfloat_1d_3_lr_randomread &x_in, + const typename AT::t_int_1d_randomread &type_in, + const typename AT::t_kkfloat_1d &d_rcov_in, + const typename AT::t_kkfloat_1d &d_cn_in, + const typename AT::t_int_1d &d_ilist_in, + const typename AT::t_int_1d &d_numneigh_in, + const typename AT::t_neighbors_2d &d_neighbors_in, + int nlocal_in, KK_FLOAT cn_thr_in) + : x(x_in), type(type_in), d_rcov(d_rcov_in), d_cn(d_cn_in), + d_ilist(d_ilist_in), d_numneigh(d_numneigh_in), + d_neighbors(d_neighbors_in), nlocal(nlocal_in), cn_thr(cn_thr_in) + { + dup_cn = Kokkos::Experimental::create_scatter_view(d_cn); + } + + void contribute() { + Kokkos::Experimental::contribute(d_cn, dup_cn); + } + + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void operator()(const int &ii) const { + auto a_cn = dup_cn.template access(); + + int i = d_ilist(ii); + int itype = type(i); + int jnum = d_numneigh(i); + + const KK_FLOAT xtmp = x(i, 0); + const KK_FLOAT ytmp = x(i, 1); + const KK_FLOAT ztmp = x(i, 2); + + KK_FLOAT cn_i = 0.0; + + for (int jj = 0; jj < jnum; jj++) { + int j = d_neighbors(i,jj); + j &= NEIGHMASK; + int jtype = type(j); + + const KK_FLOAT delx = xtmp - x(j, 0); + const KK_FLOAT dely = ytmp - x(j, 1); + const KK_FLOAT delz = ztmp - x(j, 2); + const KK_FLOAT rsq = delx * delx + dely * dely + delz * delz; + + if (rsq > cn_thr) continue; + + const KK_FLOAT rr = Kokkos::sqrt(rsq); + const KK_FLOAT rcov_ij = (d_rcov(itype) + d_rcov(jtype)) * AUTOANG; + const KK_FLOAT cn_ij = 1.0 / (1.0 + Kokkos::exp(-K1 * ((rcov_ij / rr) - 1.0))); + + cn_i += cn_ij; + if (NEIGHFLAG != FULL && (NEWTON_PAIR || j < nlocal)) a_cn(j) += cn_ij; + } + + a_cn(i) += cn_i; + } +}; + + +/* ---------------------------------------------------------------------- + Kernel A: compute energy/force and dC6 +------------------------------------------------------------------------- */ + +template +struct PairDispersionD3KernelA { + typedef ArrayTypes AT; + using value_type = EV_FLOAT; + using DUP = NeedDup_v; + using ScatterAccess = std::conditional_t< + std::is_same_v, + Kokkos::Experimental::ScatterNonAtomic, + Kokkos::Experimental::ScatterAtomic>; + + typedef typename KKDevice::value KKDeviceType; + template + using DupScatterView = KKScatterView; + template + using NonDupScatterView = KKScatterView; + + typename AT::t_kkfloat_1d_3_lr_randomread x; + typename AT::t_int_1d_randomread type; + typename AT::t_kkfloat_2d d_cutsq; + typename AT::t_kkfloat_1d d_cn; + typename AT::t_kkfloat_1d d_dc6; + typename AT::t_kkfloat_1d d_r2r4; + typename AT::t_kkfloat_2d d_r0ab; + Kokkos::View d_c6ab; + typename AT::t_int_1d d_mxci; + + typename AT::t_int_1d d_numneigh; + typename AT::t_neighbors_2d d_neighbors; + typename AT::t_int_1d d_ilist; + + DupScatterView dup_f; + NonDupScatterView ndup_f; + DupScatterView dup_eatom; + NonDupScatterView ndup_eatom; + DupScatterView dup_vatom; + NonDupScatterView ndup_vatom; + DupScatterView dup_dc6; + NonDupScatterView ndup_dc6; + + KK_FLOAT special_lj[4]; + int nlocal; + int eflag; + int vflag_either; + int eflag_global; + int eflag_atom; + int vflag_global; + int vflag_atom; + + int dampingCode; + KK_FLOAT s6, s8, rs6, rs8, a1, a2, alpha; + + PairDispersionD3KernelA( + const typename AT::t_kkfloat_1d_3_lr_randomread &x_in, + const typename AT::t_int_1d_randomread &type_in, + const typename AT::t_kkfloat_2d &d_cutsq_in, + const typename AT::t_kkfloat_1d &d_cn_in, + const typename AT::t_kkfloat_1d &d_dc6_in, + const typename AT::t_kkfloat_1d &d_r2r4_in, + const typename AT::t_kkfloat_2d &d_r0ab_in, + const Kokkos::View &d_c6ab_in, + const typename AT::t_int_1d &d_mxci_in, + const typename AT::t_int_1d &d_numneigh_in, + const typename AT::t_neighbors_2d &d_neighbors_in, + const typename AT::t_int_1d &d_ilist_in, + const DupScatterView &dup_f_in, + const NonDupScatterView &ndup_f_in, + const DupScatterView &dup_eatom_in, + const NonDupScatterView &ndup_eatom_in, + const DupScatterView &dup_vatom_in, + const NonDupScatterView &ndup_vatom_in, + const DupScatterView &dup_dc6_in, + const NonDupScatterView &ndup_dc6_in, + const KK_FLOAT *special_lj_in, + int nlocal_in, int eflag_in, int vflag_either_in, + int eflag_global_in, int eflag_atom_in, int vflag_global_in, int vflag_atom_in, + int dampingCode_in, KK_FLOAT s6_in, KK_FLOAT s8_in, + KK_FLOAT rs6_in, KK_FLOAT rs8_in, KK_FLOAT a1_in, KK_FLOAT a2_in, + KK_FLOAT alpha_in) + : x(x_in), type(type_in), d_cutsq(d_cutsq_in), d_cn(d_cn_in), d_dc6(d_dc6_in), + d_r2r4(d_r2r4_in), d_r0ab(d_r0ab_in), + d_c6ab(d_c6ab_in), d_mxci(d_mxci_in), + d_numneigh(d_numneigh_in), d_neighbors(d_neighbors_in), d_ilist(d_ilist_in), + dup_f(dup_f_in), ndup_f(ndup_f_in), + dup_eatom(dup_eatom_in), ndup_eatom(ndup_eatom_in), + dup_vatom(dup_vatom_in), ndup_vatom(ndup_vatom_in), + dup_dc6(dup_dc6_in), ndup_dc6(ndup_dc6_in), + nlocal(nlocal_in), eflag(eflag_in), vflag_either(vflag_either_in), + eflag_global(eflag_global_in), eflag_atom(eflag_atom_in), + vflag_global(vflag_global_in), vflag_atom(vflag_atom_in), + dampingCode(dampingCode_in), s6(s6_in), s8(s8_in), + rs6(rs6_in), rs8(rs8_in), a1(a1_in), a2(a2_in), + alpha(alpha_in) + { + special_lj[0] = special_lj_in[0]; + special_lj[1] = special_lj_in[1]; + special_lj[2] = special_lj_in[2]; + special_lj[3] = special_lj_in[3]; + } + + // Extract special bond mask + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + static int sbmask(const int &j) { + return j >> SBBITS & 3; + } + + /* ---------------------------------------------------------------------- + Get derivative of C6 on device + ------------------------------------------------------------------------- */ + + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void get_dC6_kokkos(KK_FLOAT &c6, KK_FLOAT &dc6i, KK_FLOAT &dc6j, + const int &iat, const int &jat, const KK_FLOAT &cni, const KK_FLOAT &cnj) const + { + KK_FLOAT c6_ref, cni_ref, cnj_ref; + KK_FLOAT c6mem, r_save, r; + KK_FLOAT expterm, term; + KK_FLOAT num, den, d_num_i, d_num_j, d_den_i, d_den_j; + + c6mem = -1.0e20; + r_save = 1.0e20; + num = 0.0; + den = 0.0; + d_num_i = 0.0; + d_num_j = 0.0; + d_den_i = 0.0; + d_den_j = 0.0; + + int maxci = d_mxci(iat); + int maxcj = d_mxci(jat); + for (int ci = 0; ci <= maxci; ci++) { + for (int cj = 0; cj <= maxcj; cj++) { + c6_ref = d_c6ab(iat, jat, ci, cj, 0); + c6_ref *= AUTOEV * AUTOANG6; + + if (c6_ref > 0) { + cni_ref = d_c6ab(iat, jat, ci, cj, 1); + cnj_ref = d_c6ab(iat, jat, ci, cj, 2); + + r = (cni - cni_ref) * (cni - cni_ref) + (cnj - cnj_ref) * (cnj - cnj_ref); + + if (r < r_save) { + r_save = r; + c6mem = c6_ref; + } + + expterm = Kokkos::exp(static_cast(K3) * r); + + num += c6_ref * expterm; + den += expterm; + + expterm = expterm * static_cast(2.0 * K3); + + term = expterm * (cni - cni_ref); + d_num_i += c6_ref * term; + d_den_i += term; + + term = expterm * (cnj - cnj_ref); + d_num_j += c6_ref * term; + d_den_j += term; + } + } + } + + // the reference threshold of 1.0e-99 underflows to zero in a single + // precision build, so use the smallest normalized value of KK_FLOAT + + if (den > std::numeric_limits::min()) { + c6 = num / den; + dc6i = ((d_num_i * den) - (d_den_i * num)) / (den * den); + dc6j = ((d_num_j * den) - (d_den_j * num)) / (den * den); + } else { + c6 = c6mem; + dc6i = 0; + dc6j = 0; + } + } + + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void ev_tally(value_type &ev, const int &i, const int &j, + const KK_FLOAT &epair, const KK_FLOAT &fpair, const KK_FLOAT &delx, + const KK_FLOAT &dely, const KK_FLOAT &delz) const + { + const int EFLAG = eflag; + const int VFLAG = vflag_either; + + if (EFLAG) { + const KK_FLOAT epairhalf = 0.5 * epair; + if (eflag_global) { + if (NEIGHFLAG!=FULL && (NEWTON_PAIR || j < nlocal)) { + ev.evdwl += epair; + } else { + ev.evdwl += epairhalf; + } + } + if (eflag_atom) { + auto v_eatom = ScatterViewHelper::get(dup_eatom,ndup_eatom); + auto a_eatom = v_eatom.template access(); + if (NEIGHFLAG!=FULL) { + if (NEWTON_PAIR || i < nlocal) a_eatom[i] += epairhalf; + if (NEWTON_PAIR || j < nlocal) a_eatom[j] += epairhalf; + } else { + a_eatom[i] += epairhalf; + } + } + } + + if (VFLAG) { + const KK_FLOAT v0 = delx*delx*fpair; + const KK_FLOAT v1 = dely*dely*fpair; + const KK_FLOAT v2 = delz*delz*fpair; + const KK_FLOAT v3 = delx*dely*fpair; + const KK_FLOAT v4 = delx*delz*fpair; + const KK_FLOAT v5 = dely*delz*fpair; + + if (vflag_global) { + if (NEIGHFLAG!=FULL) { + if (NEWTON_PAIR || i < nlocal) { + ev.v[0] += 0.5*v0; + ev.v[1] += 0.5*v1; + ev.v[2] += 0.5*v2; + ev.v[3] += 0.5*v3; + ev.v[4] += 0.5*v4; + ev.v[5] += 0.5*v5; + } + if (NEWTON_PAIR || j < nlocal) { + ev.v[0] += 0.5*v0; + ev.v[1] += 0.5*v1; + ev.v[2] += 0.5*v2; + ev.v[3] += 0.5*v3; + ev.v[4] += 0.5*v4; + ev.v[5] += 0.5*v5; + } + } else { + ev.v[0] += 0.5*v0; + ev.v[1] += 0.5*v1; + ev.v[2] += 0.5*v2; + ev.v[3] += 0.5*v3; + ev.v[4] += 0.5*v4; + ev.v[5] += 0.5*v5; + } + } + + if (vflag_atom) { + auto v_vatom = ScatterViewHelper::get(dup_vatom,ndup_vatom); + auto a_vatom = v_vatom.template access(); + if (NEIGHFLAG!=FULL) { + if (NEWTON_PAIR || i < nlocal) { + a_vatom(i,0) += 0.5*v0; + a_vatom(i,1) += 0.5*v1; + a_vatom(i,2) += 0.5*v2; + a_vatom(i,3) += 0.5*v3; + a_vatom(i,4) += 0.5*v4; + a_vatom(i,5) += 0.5*v5; + } + if (NEWTON_PAIR || j < nlocal) { + a_vatom(j,0) += 0.5*v0; + a_vatom(j,1) += 0.5*v1; + a_vatom(j,2) += 0.5*v2; + a_vatom(j,3) += 0.5*v3; + a_vatom(j,4) += 0.5*v4; + a_vatom(j,5) += 0.5*v5; + } + } else { + a_vatom(i,0) += 0.5*v0; + a_vatom(i,1) += 0.5*v1; + a_vatom(i,2) += 0.5*v2; + a_vatom(i,3) += 0.5*v3; + a_vatom(i,4) += 0.5*v4; + a_vatom(i,5) += 0.5*v5; + } + } + } + } + + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void operator()(const int &ii, value_type &ev) const { + auto v_f = ScatterViewHelper::get(dup_f,ndup_f); + auto a_f = v_f.template access(); + + auto v_dc6 = ScatterViewHelper::get(dup_dc6,ndup_dc6); + auto a_dc6 = v_dc6.template access(); + + const int i = d_ilist(ii); + const int itype = type(i); + int jnum = d_numneigh(i); + + const KK_FLOAT xtmp = x(i, 0); + const KK_FLOAT ytmp = x(i, 1); + const KK_FLOAT ztmp = x(i, 2); + + KK_FLOAT fxtmp = 0.0; + KK_FLOAT fytmp = 0.0; + KK_FLOAT fztmp = 0.0; + KK_FLOAT dc6_i = 0.0; + + for (int jj = 0; jj < jnum; jj++) { + int j = d_neighbors(i, jj); + KK_FLOAT factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + const KK_FLOAT delx = xtmp - x(j, 0); + const KK_FLOAT dely = ytmp - x(j, 1); + const KK_FLOAT delz = ztmp - x(j, 2); + + const KK_FLOAT rsq = delx * delx + dely * dely + delz * delz; + + const int jtype = type(j); + + if (rsq < d_cutsq(itype, jtype)) { + + const KK_FLOAT r2inv = 1.0f / rsq; + const KK_FLOAT r6inv = r2inv * r2inv * r2inv; + const KK_FLOAT r8inv = r6inv * r2inv; + const KK_FLOAT r10inv = r8inv * r2inv; + + KK_FLOAT c6 = 0.0; + KK_FLOAT dc6i = 0.0; + KK_FLOAT dc6j = 0.0; + get_dC6_kokkos(c6, dc6i, dc6j, itype, jtype, d_cn(i), d_cn(j)); + + const KK_FLOAT C6 = c6; + const KK_FLOAT C8 = 3.0 * C6 * d_r2r4(itype) * d_r2r4(jtype) * AUTOANG * AUTOANG; + + const KK_FLOAT alpha6 = alpha; + const KK_FLOAT alpha8 = alpha + 2; + + KK_FLOAT t6, t8, damp6, damp8, e6, e8; + KK_FLOAT tmp6, tmp8, fpair1, fpair2, fpair; + KK_FLOAT evdwl = 0.0; + t6 = t8 = e6 = e8 = fpair = fpair1 = fpair2 = 0.0; + + switch (dampingCode) { + + case 1: { // original + + const KK_FLOAT ip6 = rs6 * d_r0ab(itype, jtype); + const KK_FLOAT ip8 = rs8 * d_r0ab(itype, jtype); + + const KK_FLOAT half_alpha6 = 0.5 * alpha6; + const KK_FLOAT half_alpha8 = 0.5 * alpha8; + + t6 = MathSpecialKokkos::powauto(ip6, alpha6) * MathSpecialKokkos::powauto(rsq, -half_alpha6); + t8 = MathSpecialKokkos::powauto(ip8, alpha8) * MathSpecialKokkos::powauto(rsq, -half_alpha8); + + damp6 = 1.0f / (1.0f + 6.0f * t6); + damp8 = 1.0f / (1.0f + 6.0f * t8); + + e6 = C6 * damp6 * r6inv; + e8 = C8 * damp8 * r8inv; + + tmp6 = 6 * s6 * C6 * r8inv * damp6; + tmp8 = 8 * s8 * C8 * r10inv * damp8; + + fpair1 = -tmp6 - tmp8; + fpair2 = tmp6 * alpha6 * t6 * damp6 + (3.0f / 4) * tmp8 * alpha8 * t8 * damp8; + + fpair = fpair1 + fpair2; + fpair *= factor_lj; + + } break; + + case 2: { // zerom + + const KK_FLOAT r0 = d_r0ab(itype, jtype); + const KK_FLOAT r = Kokkos::sqrt(rsq); + + t6 = MathSpecialKokkos::powauto((r / (rs6 * r0)) + rs8 * r0, -alpha6); + damp6 = 1.0f / (1.0f + 6.0f * t6); + t8 = MathSpecialKokkos::powauto((r / r0) + rs8 * r0, -alpha8); + damp8 = 1.0f / (1.0f + 6.0f * t8); + + e6 = C6 * damp6 * r6inv; + e8 = C8 * damp8 * r8inv; + + tmp6 = 6 * s6 * C6 * r8inv * damp6; + tmp8 = 8 * s8 * C8 * r10inv * damp8; + + fpair1 = -tmp6 - tmp8; + + const KK_FLOAT fp26 = tmp6 * alpha6 * t6 * damp6 * r / (r + rs6 * rs8 * r0 * r0); + const KK_FLOAT fp28 = tmp8 * alpha8 * t8 * damp8 * r / (r + rs8 * r0 * r0); + + fpair2 = fp26 + (3.0f / 4) * fp28; + + fpair = fpair1 + fpair2; + fpair *= factor_lj; + } break; + + case 3: // bj + case 4: { // bjm, same functional form as bj, different parameters + + const KK_FLOAT r0 = Kokkos::sqrt(C8 / C6); + + const KK_FLOAT r4 = rsq * rsq; + KK_FLOAT r6 = rsq * rsq * rsq; + KK_FLOAT r8 = r6 * rsq; + + const KK_FLOAT d = a1 * r0 + a2; + const KK_FLOAT d2 = d * d; + const KK_FLOAT d4 = d2 * d2; + + t6 = r6 + MathSpecialKokkos::cube(d2); + t8 = r8 + MathSpecialKokkos::square(d4); + + e6 = C6 / t6; + e8 = C8 / t8; + + tmp6 = 6.0 * s6 * C6 * r4 / (t6 * t6); + tmp8 = 8.0 * s8 * C8 * r6 / (t8 * t8); + + fpair = -(tmp6 + tmp8); + fpair *= factor_lj; + } break; + + // no default case: dampingCode is validated in init_style() + } + + if (EVFLAG) evdwl = -(s6 * e6 + s8 * e8) * factor_lj; + + const KK_FLOAT rest = (s6 * e6 + s8 * e8) / C6; + + dc6_i += rest * dc6i; + if (NEIGHFLAG != FULL && (NEWTON_PAIR || j < nlocal)) { a_dc6(j) += rest * dc6j; } + + fxtmp += delx * fpair; + fytmp += dely * fpair; + fztmp += delz * fpair; + + if (NEIGHFLAG != FULL && (NEWTON_PAIR || j < nlocal)) { + a_f(j, 0) -= delx * fpair; + a_f(j, 1) -= dely * fpair; + a_f(j, 2) -= delz * fpair; + } + + if (EVFLAG) { + ev_tally(ev, i, j, evdwl, fpair, delx, dely, delz); + } + } + } + + a_f(i, 0) += fxtmp; + a_f(i, 1) += fytmp; + a_f(i, 2) += fztmp; + a_dc6(i) += dc6_i; + } + + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void operator()(const int &ii) const { + value_type ev; + operator()(ii, ev); + } +}; + +/* ---------------------------------------------------------------------- + Kernel B: compute force contribution from dC6 +------------------------------------------------------------------------- */ + +template +struct PairDispersionD3KernelB { + typedef ArrayTypes AT; + using value_type = EV_FLOAT; + using DUP = NeedDup_v; + using ScatterAccess = std::conditional_t< + std::is_same_v, + Kokkos::Experimental::ScatterNonAtomic, + Kokkos::Experimental::ScatterAtomic>; + + typedef typename KKDevice::value KKDeviceType; + template + using DupScatterView = KKScatterView; + template + using NonDupScatterView = KKScatterView; + + typename AT::t_kkfloat_1d_3_lr_randomread x; + typename AT::t_int_1d_randomread type; + typename AT::t_kkfloat_2d d_cutsq; + typename AT::t_kkfloat_1d d_dc6; + typename AT::t_kkfloat_1d d_rcov; + typename AT::t_int_1d d_numneigh; + typename AT::t_neighbors_2d d_neighbors; + typename AT::t_int_1d d_ilist; + + DupScatterView dup_f; + NonDupScatterView ndup_f; + DupScatterView dup_eatom; + NonDupScatterView ndup_eatom; + DupScatterView dup_vatom; + NonDupScatterView ndup_vatom; + + KK_FLOAT special_lj[4]; + int nlocal; + int eflag; + int vflag_either; + int eflag_global; + int eflag_atom; + int vflag_global; + int vflag_atom; + KK_FLOAT cn_thr; + + PairDispersionD3KernelB( + const typename AT::t_kkfloat_1d_3_lr_randomread &x_in, + const typename AT::t_int_1d_randomread &type_in, + const typename AT::t_kkfloat_2d &d_cutsq_in, + const typename AT::t_kkfloat_1d &d_dc6_in, + const typename AT::t_kkfloat_1d &d_rcov_in, + const typename AT::t_int_1d &d_numneigh_in, + const typename AT::t_neighbors_2d &d_neighbors_in, + const typename AT::t_int_1d &d_ilist_in, + const DupScatterView &dup_f_in, + const NonDupScatterView &ndup_f_in, + const DupScatterView &dup_eatom_in, + const NonDupScatterView &ndup_eatom_in, + const DupScatterView &dup_vatom_in, + const NonDupScatterView &ndup_vatom_in, + const KK_FLOAT *special_lj_in, + int nlocal_in, int eflag_in, int vflag_either_in, + int eflag_global_in, int eflag_atom_in, int vflag_global_in, int vflag_atom_in, + KK_FLOAT cn_thr_in) + : x(x_in), type(type_in), d_cutsq(d_cutsq_in), d_dc6(d_dc6_in), + d_rcov(d_rcov_in), d_numneigh(d_numneigh_in), d_neighbors(d_neighbors_in), + d_ilist(d_ilist_in), dup_f(dup_f_in), ndup_f(ndup_f_in), + dup_eatom(dup_eatom_in), ndup_eatom(ndup_eatom_in), + dup_vatom(dup_vatom_in), ndup_vatom(ndup_vatom_in), + nlocal(nlocal_in), eflag(eflag_in), vflag_either(vflag_either_in), + eflag_global(eflag_global_in), eflag_atom(eflag_atom_in), + vflag_global(vflag_global_in), vflag_atom(vflag_atom_in), + cn_thr(cn_thr_in) + { + special_lj[0] = special_lj_in[0]; + special_lj[1] = special_lj_in[1]; + special_lj[2] = special_lj_in[2]; + special_lj[3] = special_lj_in[3]; + } + + // Extract special bond mask + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + static int sbmask(const int &j) { + return j >> SBBITS & 3; + } + + + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void ev_tally(value_type &ev, const int &i, const int &j, + const KK_FLOAT &epair, const KK_FLOAT &fpair, const KK_FLOAT &delx, + const KK_FLOAT &dely, const KK_FLOAT &delz) const + { + const int EFLAG = eflag; + const int VFLAG = vflag_either; + + if (EFLAG) { + const KK_FLOAT epairhalf = 0.5 * epair; + if (eflag_global) { + if (NEIGHFLAG!=FULL && (NEWTON_PAIR || j < nlocal)) { + ev.evdwl += epair; + } else { + ev.evdwl += epairhalf; + } + } + if (eflag_atom) { + auto v_eatom = ScatterViewHelper::get(dup_eatom,ndup_eatom); + auto a_eatom = v_eatom.template access(); + if (NEIGHFLAG!=FULL) { + if (NEWTON_PAIR || i < nlocal) a_eatom[i] += epairhalf; + if (NEWTON_PAIR || j < nlocal) a_eatom[j] += epairhalf; + } else { + a_eatom[i] += epairhalf; + } + } + } + + if (VFLAG) { + const KK_FLOAT v0 = delx*delx*fpair; + const KK_FLOAT v1 = dely*dely*fpair; + const KK_FLOAT v2 = delz*delz*fpair; + const KK_FLOAT v3 = delx*dely*fpair; + const KK_FLOAT v4 = delx*delz*fpair; + const KK_FLOAT v5 = dely*delz*fpair; + + if (vflag_global) { + if (NEIGHFLAG!=FULL) { + if (NEWTON_PAIR || i < nlocal) { + ev.v[0] += 0.5*v0; + ev.v[1] += 0.5*v1; + ev.v[2] += 0.5*v2; + ev.v[3] += 0.5*v3; + ev.v[4] += 0.5*v4; + ev.v[5] += 0.5*v5; + } + if (NEWTON_PAIR || j < nlocal) { + ev.v[0] += 0.5*v0; + ev.v[1] += 0.5*v1; + ev.v[2] += 0.5*v2; + ev.v[3] += 0.5*v3; + ev.v[4] += 0.5*v4; + ev.v[5] += 0.5*v5; + } + } else { + ev.v[0] += 0.5*v0; + ev.v[1] += 0.5*v1; + ev.v[2] += 0.5*v2; + ev.v[3] += 0.5*v3; + ev.v[4] += 0.5*v4; + ev.v[5] += 0.5*v5; + } + } + + if (vflag_atom) { + auto v_vatom = ScatterViewHelper::get(dup_vatom,ndup_vatom); + auto a_vatom = v_vatom.template access(); + if (NEIGHFLAG!=FULL) { + if (NEWTON_PAIR || i < nlocal) { + a_vatom(i,0) += 0.5*v0; + a_vatom(i,1) += 0.5*v1; + a_vatom(i,2) += 0.5*v2; + a_vatom(i,3) += 0.5*v3; + a_vatom(i,4) += 0.5*v4; + a_vatom(i,5) += 0.5*v5; + } + if (NEWTON_PAIR || j < nlocal) { + a_vatom(j,0) += 0.5*v0; + a_vatom(j,1) += 0.5*v1; + a_vatom(j,2) += 0.5*v2; + a_vatom(j,3) += 0.5*v3; + a_vatom(j,4) += 0.5*v4; + a_vatom(j,5) += 0.5*v5; + } + } else { + a_vatom(i,0) += 0.5*v0; + a_vatom(i,1) += 0.5*v1; + a_vatom(i,2) += 0.5*v2; + a_vatom(i,3) += 0.5*v3; + a_vatom(i,4) += 0.5*v4; + a_vatom(i,5) += 0.5*v5; + } + } + } + } + + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void operator()(const int &ii, value_type &ev) const { + auto v_f = ScatterViewHelper::get(dup_f,ndup_f); + auto a_f = v_f.template access(); + + const int i = d_ilist(ii); + const int itype = type(i); + int jnum = d_numneigh(i); + + const KK_FLOAT xtmp = x(i, 0); + const KK_FLOAT ytmp = x(i, 1); + const KK_FLOAT ztmp = x(i, 2); + + KK_FLOAT fxtmp = 0.0; + KK_FLOAT fytmp = 0.0; + KK_FLOAT fztmp = 0.0; + + for (int jj = 0; jj < jnum; jj++) { + int j = d_neighbors(i, jj); + KK_FLOAT factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + const KK_FLOAT delx = xtmp - x(j, 0); + const KK_FLOAT dely = ytmp - x(j, 1); + const KK_FLOAT delz = ztmp - x(j, 2); + const KK_FLOAT rsq = delx * delx + dely * dely + delz * delz; + + const int jtype = type(j); + + if (rsq < d_cutsq(itype, jtype)) { + + const KK_FLOAT r = Kokkos::sqrt(rsq); + KK_FLOAT dcn; + + if (rsq < cn_thr) { + const KK_FLOAT rcovij = (d_rcov(itype) + d_rcov(jtype)) * AUTOANG; + const KK_FLOAT expterm = Kokkos::exp(-K1 * (rcovij / r - 1.0)); + dcn = -K1 * rcovij * expterm / (rsq * (expterm + 1.0) * (expterm + 1.0)); + } else { + dcn = 0.0; + } + + KK_FLOAT fpair = dcn * (d_dc6(i) + d_dc6(j)) / r; + fpair *= factor_lj; + + fxtmp += delx * fpair; + fytmp += dely * fpair; + fztmp += delz * fpair; + + if (NEIGHFLAG != FULL && (NEWTON_PAIR || j < nlocal)) { + a_f(j, 0) -= delx * fpair; + a_f(j, 1) -= dely * fpair; + a_f(j, 2) -= delz * fpair; + } + + if (EVFLAG) { + const KK_FLOAT epair = 0.0; + ev_tally(ev, i, j, epair, fpair, delx, dely, delz); + } + } + } + + a_f(i, 0) += fxtmp; + a_f(i, 1) += fytmp; + a_f(i, 2) += fztmp; + } + + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void operator()(const int &ii) const { + value_type ev; + operator()(ii, ev); + } +}; + +template +class PairDispersionD3Kokkos : public PairDispersionD3, public KokkosBase { + public: + enum {EnabledNeighFlags=FULL|HALFTHREAD|HALF}; + enum {COUL_FLAG=0}; + typedef DeviceType device_type; + typedef ArrayTypes AT; + typedef ArrayTypes HAT; + typedef EV_FLOAT value_type; + PairDispersionD3Kokkos(class LAMMPS *); + ~PairDispersionD3Kokkos() override; + + void calc_coordination_number() override; + void compute(int, int) override; + double init_one(int, int) override; + void init_style() override; + void allocate() override; + void coeff(int, char **) override; + + int pack_forward_comm_kokkos(int, DAT::tdual_int_1d, DAT::tdual_double_1d&, + int, int *) override; + void unpack_forward_comm_kokkos(int, int, DAT::tdual_double_1d&) override; + int pack_reverse_comm_kokkos(int, int, DAT::tdual_double_1d&) override; + void unpack_reverse_comm_kokkos(int, DAT::tdual_int_1d, DAT::tdual_double_1d&) override; + int pack_forward_comm(int, int *, double *, int, int *) override; + void unpack_forward_comm(int, int, double *) override; + int pack_reverse_comm(int, int, double *) override; + void unpack_reverse_comm(int, int *, double *) override; + + protected: + typename AT::t_kkfloat_1d_3_lr_randomread x; // atom positions + typename AT::t_kkacc_1d_3 f; // atom forces + typename AT::t_int_1d_randomread type; // atom types + + DAT::ttransform_kkacc_1d k_eatom; // per-atom energy (dual view) + DAT::ttransform_kkacc_1d_6 k_vatom; // per-atom virial (dual view) + typename AT::t_kkacc_1d d_eatom; // device view of per-atom energy + typename AT::t_kkacc_1d_6 d_vatom; // device view of per-atom virial + + KK_FLOAT special_lj[4]; // special-bond scaling + int inum; // number of neighbor list atoms + bool need_dup; // whether duplicated scatter is required + + typedef typename KKDevice::value KKDeviceType; + template + using DupScatterView = KKScatterView; + + template + using NonDupScatterView = KKScatterView; + + DupScatterView dup_f; // duplicated force + DupScatterView dup_eatom; // duplicated energy + DupScatterView dup_vatom; // duplicated virial + DupScatterView dup_dc6; // duplicated dC6 + NonDupScatterView ndup_f; // non-dup force + NonDupScatterView ndup_eatom; // non-dup energy + NonDupScatterView ndup_vatom; // non-dup virial + NonDupScatterView ndup_dc6; // non-dup dC6 + + + DAT::tdual_kkfloat_1d k_cn; // coordination numbers (dual view) + DAT::tdual_kkfloat_1d k_dc6; // dC6 values (dual view) + typename AT::t_kkfloat_1d d_cn; // device CN + typename AT::t_kkfloat_1d d_dc6; // device dC6 + HAT::t_kkfloat_1d h_cn; // host CN + HAT::t_kkfloat_1d h_dc6; // host dC6 + + DAT::tdual_kkfloat_1d k_r2r4; // r2r4 table (dual view) + DAT::tdual_kkfloat_1d k_rcov; // covalent radii (dual view) + DAT::tdual_int_1d k_mxci; // max C6 grid index (dual view) + DAT::tdual_kkfloat_2d k_r0ab; // R0 table (dual view) + Kokkos::DualView k_c6ab; // C6 table (dual view) + + typename AT::t_kkfloat_1d d_r2r4; // device r2r4 + typename AT::t_kkfloat_1d d_rcov; // device covalent radii + typename AT::t_int_1d d_mxci; // device max C6 grid index + typename AT::t_kkfloat_2d d_r0ab; // device R0 table + Kokkos::View d_c6ab; // device C6 table + + DAT::ttransform_kkfloat_2d k_cutsq; // cutoff^2 table (host double, device KK_FLOAT) + typename AT::t_kkfloat_2d d_cutsq; // device cutoff^2 table + + typename AT::t_neighbors_2d d_neighbors; // neighbor list + typename AT::t_int_1d d_ilist; // neighbor list indices + typename AT::t_int_1d d_numneigh; // neighbor counts + + typename AT::t_int_1d d_sendlist, d_recvlist; // comm lists + typename AT::t_double_1d_um v_buf; // comm buffer + + int neighflag, newton_pair; // neighbor/newton settings + int nlocal, nall, eflag, vflag; // local/total counts and flags + + friend void pair_virial_fdotr_compute>(PairDispersionD3Kokkos*); + + // To make the compute() function cleaner: + template + void dispatch_coordination_kernel(); + + template + void dispatch_kernel_A(EV_FLOAT &ev); + + template + void launch_kernel_A(EV_FLOAT &ev); + + template + void dispatch_kernel_B(EV_FLOAT &ev); + + template + void launch_kernel_B(EV_FLOAT &ev); + +}; + +} + +#endif +#endif diff --git a/src/OPENMP/pair_dispersion_d3_omp.cpp b/src/OPENMP/pair_dispersion_d3_omp.cpp new file mode 100644 index 00000000000..27bb4cda2e2 --- /dev/null +++ b/src/OPENMP/pair_dispersion_d3_omp.cpp @@ -0,0 +1,525 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + This software is distributed under the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: + Ulrik Unneberg + Marc L. Descoteaux + Yizhong R. Hu + William C. Witt + Affiliation: Harvard University +------------------------------------------------------------------------- */ + +#include "pair_dispersion_d3_omp.h" + +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "math_special.h" +#include "memory.h" +#include "neigh_list.h" +#include "suffix.h" + +#include +#include + +#include "omp_compat.h" + +using namespace LAMMPS_NS; +using namespace LAMMPS_NS::DispersionD3; + +/* ---------------------------------------------------------------------- */ + +PairDispersionD3OMP::PairDispersionD3OMP(LAMMPS *lmp) : + PairDispersionD3(lmp), ThrOMP(lmp, THR_PAIR) +{ + suffix_flag |= Suffix::OMP; +} + +/* ---------------------------------------------------------------------- */ + +void PairDispersionD3OMP::calc_coordination_number() +{ + const int nthreads = comm->nthreads; + const int newton_pair = force->newton_pair; + + // cn and dc6 hold one copy per thread, each nall long. The threads only + // ever touch their own copy; data_reduce_thr() sums the copies into the + // first nall elements at the end of the respective loops. + + if (atom->nmax > nmax) { + nmax = atom->nmax; + memory->grow(cn, nthreads * nmax, "pair:cn"); + memory->grow(dc6, nthreads * nmax, "pair:dc6"); + } + + const int inum = list->inum; + + // Begin parallel region, the central atoms indexed by ii are assigned to different threads. + #if defined(_OPENMP) + #pragma omp parallel LMP_DEFAULT_NONE \ + firstprivate(inum,nthreads) + #endif + { + int ifrom, ito, tid; + + // Set up the starting and ending indices for each thread + loop_setup_thr(ifrom, ito, tid, inum, nthreads); + ThrData *thr = fix->get_thr(tid); + thr->timer(Timer::START); + + // Calculate coordination number with the helper functions + // The flags need to be constants for the template instantiation + if (force->newton_pair) { + eval_coordination<1>(ifrom,ito,thr); + } else{ + eval_coordination<0>(ifrom,ito,thr); + } + + thr->timer(Timer::PAIR); + } + + // communicate coordination number + communicationStage = 1; + if (newton_pair) comm->reverse_comm(this); + comm->forward_comm(this); +} + +/* ---------------------------------------------------------------------- */ + +template +void PairDispersionD3OMP::eval_coordination(int iifrom, int iito, ThrData * const thr) +{ + + const auto * _noalias const x = (dbl3_t *) atom->x[0]; + const int * _noalias const type = atom->type; + const int nlocal = atom->nlocal; + const int nall = nlocal + atom->nghost; + const int nthreads = comm->nthreads; + const int tid = thr->get_tid(); + const int * _noalias const ilist = list->ilist; + const int * _noalias const numneigh = list->numneigh; + const int * const * const firstneigh = list->firstneigh; + + // this thread's private copy of the coordination number accumulator + double * _noalias const thr_cn = cn + tid * nall; + memset(thr_cn, 0, sizeof(double) * nall); + + for (int ii = iifrom; ii < iito; ii++) { + + int i = ilist[ii]; + int itype = type[i]; + const int * _noalias const jlist = firstneigh[i]; + int jnum = numneigh[i]; + + for (int jj = 0; jj < jnum; jj++) { + + int j = jlist[jj]; + j &= NEIGHMASK; + int jtype = type[j]; + + double delrj[3]; + delrj[0] = x[i].x - x[j].x; + delrj[1] = x[i].y - x[j].y; + delrj[2] = x[i].z - x[j].z; + + double rsq = delrj[0] * delrj[0] + delrj[1] * delrj[1] + delrj[2] * delrj[2]; + + // if the atoms are too far away don't consider the contribution + if (rsq > cn_thr) continue; + + double rr = sqrt(rsq); + double rcov_ij = (rcov[itype] + rcov[jtype]) * AUTOANG; + double cn_ij = 1.0 / (1.0 + exp(-K1 * ((rcov_ij / rr) - 1.0))); + + // update coordination number on a thread-local array + thr_cn[i] += cn_ij; + if (NEWTON_PAIR || j < nlocal) { thr_cn[j] += cn_ij; } + } + } + + // sum the per thread copies into cn[0] ... cn[nall-1] + + sync_threads(); + data_reduce_thr(cn, nall, nthreads, 1, tid); + sync_threads(); +} + +/* ---------------------------------------------------------------------- */ + +void PairDispersionD3OMP::compute(int eflag, int vflag) +{ + ev_init(eflag, vflag); + + // dampingCode is validated in PairDispersionD3::init_style(), so the loops + // below need no default case (error->all() must not be called from inside a + // parallel region) + + // First call coordination number calculation + calc_coordination_number(); + + const int nall = atom->nlocal + atom->nghost; + const int nthreads = comm->nthreads; + const int inum = list->inum; + +// Parallel direct force computation and some other quantities calculation. +#if defined(_OPENMP) +#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag) \ +firstprivate(inum,nthreads,nall) +#endif + { + int ifrom, ito, tid; + + loop_setup_thr(ifrom, ito, tid, inum, nthreads); + ThrData *thr = fix->get_thr(tid); + thr->timer(Timer::START); + ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr); + + // Call the helper eval function with the appropriate flags for the first phase of the computation + // Again, the flags need to be constants for the template instantiation + if (evflag) { + if (eflag) { + if (force->newton_pair) { + eval_first_phase<1,1,1>(ifrom, ito, thr); + } + else { + eval_first_phase<1,1,0>(ifrom, ito, thr); + } + } else { + if (force->newton_pair) eval_first_phase<1,0,1>(ifrom, ito, thr); + else eval_first_phase<1,0,0>(ifrom, ito, thr); + } + } else { + if (force->newton_pair) eval_first_phase<0,0,1>(ifrom, ito, thr); + else eval_first_phase<0,0,0>(ifrom, ito, thr); + } + thr->timer(Timer::PAIR); + } // end of omp parallel region + + // Both phases tally into the same ThrData, so ev_setup_thr() is called only + // in the first region (it zeroes the per thread accumulators) and + // reduce_thr() only at the end of the second one. The per thread force + // arrays therefore stay unreduced across the communication below, which is + // safe because it only exchanges dc6. + + // Communication stage 2 for dc6 values in preparation for calculation of indirect forces in the second phase + communicationStage = 2; + if (force->newton_pair) { + comm->reverse_comm(this); + } + + comm->forward_comm(this); + + // Process the second phase with the combined dc6 values + #if defined(_OPENMP) + #pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag) \ + firstprivate(inum,nthreads,nall) + #endif + { + int ifrom, ito, tid; + + loop_setup_thr(ifrom, ito, tid, inum, nthreads); + ThrData *thr = fix->get_thr(tid); + thr->timer(Timer::START); + + // Call the helper eval function with the appropriate flags for the second phase of the computation + // Again, the flags need to be constants for the template instantiation + if (evflag) { + if (eflag) { + if (force->newton_pair) { + eval_second_phase<1,1,1>(ifrom, ito, thr); + } + else { + eval_second_phase<1,1,0>(ifrom, ito, thr); + } + } else { + if (force->newton_pair) eval_second_phase<1,0,1>(ifrom, ito, thr); + else eval_second_phase<1,0,0>(ifrom, ito, thr); + } + } else { + if (force->newton_pair) eval_second_phase<0,0,1>(ifrom, ito, thr); + else eval_second_phase<0,0,0>(ifrom, ito, thr); + } + thr->timer(Timer::PAIR); + reduce_thr(this, eflag, vflag, thr); + } //end of omp parallel region +} + +/* ---------------------------------------------------------------------- */ + +template +void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const thr) +{ + const auto * _noalias const x = (dbl3_t *) atom->x[0]; + auto * _noalias const f = (dbl3_t *) thr->get_f()[0]; + const int * _noalias const type = atom->type; + const int nlocal = atom->nlocal; + const double * _noalias const special_lj = force->special_lj; + const int * _noalias const ilist = list->ilist; + const int * _noalias const numneigh = list->numneigh; + const int * const * const firstneigh = list->firstneigh; + const int nall = nlocal + atom->nghost; + const int nthreads = comm->nthreads; + const int tid = thr->get_tid(); + double evdwl = 0.0; + + // this thread's private copy of the dE/dC6 accumulator + double * _noalias const thr_dc6 = dc6 + tid * nall; + memset(thr_dc6, 0, sizeof(double) * nall); + + // Loop over assigned atoms + for (int ii = iifrom; ii < iito; ++ii) { + int i = ilist[ii]; + + double xtmp = x[i].x; + double ytmp = x[i].y; + double ztmp = x[i].z; + int itype = type[i]; + int jnum = numneigh[i]; + const int * _noalias const jlist = firstneigh[i]; + + for (int jj = 0; jj < jnum; jj++) { + int j = jlist[jj]; + double factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + double delx = xtmp - x[j].x; + double dely = ytmp - x[j].y; + double delz = ztmp - x[j].z; + + double rsq = delx * delx + dely * dely + delz * delz; + + int jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + + double r2inv = 1.0 / rsq; + double r6inv = r2inv * r2inv * r2inv; + double r8inv = r2inv * r2inv * r2inv * r2inv; + double r10inv = r2inv * r2inv * r2inv * r2inv * r2inv; + + // get_dC6 writes {C6, dC6/dCN_i, dC6/dCN_j} + double c6_res[3] = {}; + get_dC6(itype, jtype, cn[i], cn[j], c6_res); + + double C6 = c6_res[0]; + double C8 = 3.0 * C6 * r2r4[itype] * r2r4[jtype] * AUTOANG * AUTOANG; + + double alpha6 = alpha; + double alpha8 = alpha + 2; + + double t6, t8, damp6, damp8, e6, e8; + double tmp6, tmp8, fpair1, fpair2, fpair; + t6 = t8 = e6 = e8 = evdwl = fpair = fpair1 = fpair2 = 0.0; + + // Damping code selection - now using the passed dampingCode parameter + switch (dampingCode) { + // Written to avoid using sqrt and pow() + case 1: /* Original damping */ + { + double ip6 = rs6 * r0ab[type[i]][type[j]]; + double ip8 = rs8 * r0ab[type[i]][type[j]]; + + double half_alpha6 = 0.5 * alpha6; + double half_alpha8 = 0.5 * alpha8; + + t6 = MathSpecial::powauto(ip6, alpha6) * MathSpecial::powauto(rsq, -half_alpha6); + t8 = MathSpecial::powauto(ip8, alpha8) * MathSpecial::powauto(rsq, -half_alpha8); + + damp6 = 1.0 / (1.0 + 6.0 * t6); + damp8 = 1.0 / (1.0 + 6.0 * t8); + + e6 = C6 * damp6 * r6inv; + e8 = C8 * damp8 * r8inv; + + tmp6 = 6 * s6 * C6 * r8inv * damp6; + tmp8 = 8 * s8 * C8 * r10inv * damp8; + + fpair1 = -tmp6 - tmp8; + fpair2 = tmp6 * alpha6 * t6 * damp6 + (3.0 / 4.0) * tmp8 * alpha8 * t8 * damp8; + + fpair = fpair1 + fpair2; + fpair *= factor_lj; + } break; + // Written to avoid pow + case 2: { // zerom + + double r = sqrt(rsq); + double r0 = r0ab[type[i]][type[j]]; + + t6 = MathSpecial::powauto((r / (rs6 * r0)) + rs8 * r0, -alpha6); + damp6 = 1.0 / (1.0 + 6.0 * t6); + t8 = MathSpecial::powauto((r / r0) + rs8 * r0, -alpha8); + damp8 = 1.0 / (1.0 + 6.0 * t8); + + e6 = C6 * damp6 * r6inv; + e8 = C8 * damp8 * r8inv; + + tmp6 = 6 * s6 * C6 * r8inv * damp6; + tmp8 = 8 * s8 * C8 * r10inv * damp8; + + fpair1 = -tmp6 - tmp8; + + double fp26 = tmp6 * alpha6 * t6 * damp6 * r / (r + rs6 * rs8 * r0 * r0); + double fp28 = tmp8 * alpha8 * t8 * damp8 * r / (r + rs8 * r0 * r0); + + fpair2 = fp26 + (3.0 / 4.0) * fp28; + + fpair = fpair1 + fpair2; + fpair *= factor_lj; + } break; + + case 3: // bj + case 4: { // bjm, same functional form as bj, different parameters + double r0 = sqrt(C8 / C6); + + double r4 = rsq * rsq; + double r6 = rsq * rsq * rsq; + double r8 = rsq * rsq * rsq * rsq; + + double d = a1 * r0 + a2; + double d2 = d * d; + double d4 = d2 * d2; + + t6 = r6 + MathSpecial::cube(d2); + t8 = r8 + MathSpecial::square(d4); + + e6 = C6 / t6; + e8 = C8 / t8; + + tmp6 = 6.0 * s6 * C6 * r4 / (t6 * t6); + tmp8 = 8.0 * s8 * C8 * r6 / (t8 * t8); + + fpair = -(tmp6 + tmp8); + fpair *= factor_lj; + } break; + } + + if (EFLAG) evdwl = -(s6 * e6 + s8 * e8) * factor_lj; + + double rest = (s6 * e6 + s8 * e8) / C6; + + // Update thread-local dc6 + double dc6_contrib_i = rest * c6_res[1]; + thr_dc6[i] += dc6_contrib_i; + + if (NEWTON_PAIR || j < nlocal) { + double dc6_contrib_j = rest * c6_res[2]; + thr_dc6[j] += dc6_contrib_j; + } + + f[i].x += delx * fpair; + f[i].y += dely * fpair; + f[i].z += delz * fpair; + + if (NEWTON_PAIR || j < nlocal) { + f[j].x -= delx * fpair; + f[j].y -= dely * fpair; + f[j].z -= delz * fpair; + } + + // Update energy and virial + if (EVFLAG) ev_tally_thr(this, i, j, nlocal, NEWTON_PAIR, evdwl, 0.0, fpair, delx, dely, delz, thr); + } + } + } + + // sum the per thread copies into dc6[0] ... dc6[nall-1] + + sync_threads(); + data_reduce_thr(dc6, nall, nthreads, 1, tid); + sync_threads(); +} + +/* ---------------------------------------------------------------------- */ + +template +void PairDispersionD3OMP::eval_second_phase(int iifrom, int iito, ThrData * const thr) +{ + const auto * _noalias const x = (dbl3_t *) atom->x[0]; + auto * _noalias const f = (dbl3_t *) thr->get_f()[0]; + const int * _noalias const type = atom->type; + const int nlocal = atom->nlocal; + const double * _noalias const special_lj = force->special_lj; + const int * _noalias const ilist = list->ilist; + const int * _noalias const numneigh = list->numneigh; + const int * const * const firstneigh = list->firstneigh; + + double xtmp,ytmp,ztmp,delx,dely,delz,rsq,factor_lj,dcn,rcovij,expterm,fpair,fxtmp,fytmp,fztmp,r; + + // Loop over assigned center atoms + for (int ii = iifrom; ii < iito; ii++) { + int i = ilist[ii]; + int itype = type[i]; + xtmp = x[i].x; + ytmp = x[i].y; + ztmp = x[i].z; + + int jnum = numneigh[i]; + const int * _noalias const jlist = firstneigh[i]; + fxtmp=fytmp=fztmp=0.0; + // Neighbor atom + for (int jj = 0; jj < jnum; jj++) { + int j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j].x; + dely = ytmp - x[j].y; + delz = ztmp - x[j].z; + + rsq = delx * delx + dely * dely + delz * delz; + int jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + + if (rsq < cn_thr) { + rcovij = (rcov[type[i]] + rcov[type[j]]) * AUTOANG; + expterm = exp(-K1 * (rcovij / r - 1.0)); + dcn = -K1 * rcovij * expterm / (rsq * (expterm + 1.0) * (expterm + 1.0)); + + } else { + dcn = 0.0; + } + + fpair = dcn * (dc6[i] + dc6[j]) / r; + fpair *= factor_lj; + + fxtmp += delx * fpair; + fytmp += dely * fpair; + fztmp += delz * fpair; + if (NEWTON_PAIR || j < nlocal) { + f[j].x -= delx * fpair; + f[j].y -= dely * fpair; + f[j].z -= delz * fpair; + } + + // Update virial (no energy contributions in this phase) + if (EVFLAG) ev_tally_thr(this, i, j, nlocal, NEWTON_PAIR, 0.0, 0.0, fpair, delx, dely, delz, thr); + } + } + f[i].x += fxtmp; + f[i].y += fytmp; + f[i].z += fztmp; + } +} + +/* ---------------------------------------------------------------------- */ + +double PairDispersionD3OMP::memory_usage() +{ + double bytes = memory_usage_thr(); + bytes += PairDispersionD3::memory_usage(); + // cn and dc6 hold comm->nthreads copies here, the base class counts one each + bytes += (double) (comm->nthreads - 1) * nmax * 2 * sizeof(double); + return bytes; +} diff --git a/src/OPENMP/pair_dispersion_d3_omp.h b/src/OPENMP/pair_dispersion_d3_omp.h new file mode 100644 index 00000000000..ef0cd80917b --- /dev/null +++ b/src/OPENMP/pair_dispersion_d3_omp.h @@ -0,0 +1,51 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(dispersion/d3/omp,PairDispersionD3OMP); +// clang-format on +#else + +#ifndef LMP_PAIR_DISPERSION_D3_OMP_H +#define LMP_PAIR_DISPERSION_D3_OMP_H + +#include "pair_dispersion_d3.h" +#include "thr_omp.h" + +namespace LAMMPS_NS { + +class PairDispersionD3OMP : public PairDispersionD3, public ThrOMP { + + public: + PairDispersionD3OMP(class LAMMPS *); + + void compute(int, int) override; + double memory_usage() override; + + private: + template + void eval_first_phase(int iifrom, int iito, ThrData *const thr); + + template void eval_coordination(int iifrom, int iito, ThrData *const thr); + + void calc_coordination_number() override; + + template + void eval_second_phase(int ifrom, int ito, ThrData *const thr); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/math_special.h b/src/math_special.h index 3cce45538fa..2296a2cbf58 100644 --- a/src/math_special.h +++ b/src/math_special.h @@ -164,6 +164,27 @@ namespace LAMMPS_NS::MathSpecial { return (n > 0) ? yy : 1.0 / yy; } + /* Version of pow(x,n) for an exponent that is only known at run time + * + * Dispatches to powint() when the exponent happens to have an integer + * value, and to std::pow() otherwise. Useful where a force field + * parameter is a floating point number that is an integer for most + * parameter sets. Follows the powint() convention of returning 0 for + * x == 0, which differs from std::pow() for a negative exponent. + * + * \param x base + * \param n exponent + * \return value of x^n */ + + static inline double powauto(const double x, const double n) + { + if (n == 0.0) return 1.0; + if (x == 0.0) return 0.0; + const double rounded = std::round(n); + if (n == rounded) return powint(x, static_cast(rounded)); + return std::pow(x, n); + } + /* Fast inline version of (sin(x)/x)^n as used by PPPM kspace styles * * This is an optimized function to compute (sin(x)/x)^n as frequently used by PPPM. diff --git a/unittest/force-styles/tests/manybody-pair-dispersion_d3.yaml b/unittest/force-styles/tests/manybody-pair-dispersion_d3.yaml index 46297ca6d64..e3728b77c33 100644 --- a/unittest/force-styles/tests/manybody-pair-dispersion_d3.yaml +++ b/unittest/force-styles/tests/manybody-pair-dispersion_d3.yaml @@ -8,7 +8,7 @@ prerequisites: ! | pre_commands: ! "" post_commands: ! "" input_file: in.manybody -pair_style: dispersion/d3 zero pbe 12.0 12.0 +pair_style: dispersion/d3 original pbe 12.0 12.0 pair_coeff: ! | * * Si Si Si Si Si Si Si Si extract: ! "" diff --git a/unittest/utils/test_math_special.cpp b/unittest/utils/test_math_special.cpp index 47dbe2b5ae9..892aad63a92 100644 --- a/unittest/utils/test_math_special.cpp +++ b/unittest/utils/test_math_special.cpp @@ -156,6 +156,46 @@ TEST(MathSpecial, powint_consistency_with_pow) } } +// ------------------------------------------------------------------------- +// powauto() +// ------------------------------------------------------------------------- + +TEST(MathSpecial, powauto_integer_exponent) +{ + // must agree exactly with powint() where the exponent has an integer value + EXPECT_DOUBLE_EQ(MathSpecial::powauto(2.0, 10.0), MathSpecial::powint(2.0, 10)); + EXPECT_DOUBLE_EQ(MathSpecial::powauto(1.7, 14.0), MathSpecial::powint(1.7, 14)); + EXPECT_DOUBLE_EQ(MathSpecial::powauto(1.7, -7.0), MathSpecial::powint(1.7, -7)); + EXPECT_DOUBLE_EQ(MathSpecial::powauto(-2.0, 3.0), -8.0); +} + +TEST(MathSpecial, powauto_fractional_exponent) +{ + EXPECT_DOUBLE_EQ(MathSpecial::powauto(4.0, 0.5), 2.0); + EXPECT_DOUBLE_EQ(MathSpecial::powauto(2.0, -6.5), std::pow(2.0, -6.5)); + EXPECT_DOUBLE_EQ(MathSpecial::powauto(1.7, 3.25), std::pow(1.7, 3.25)); +} + +TEST(MathSpecial, powauto_edge_cases) +{ + // zero exponent wins over zero base, as in powint() and std::pow() + EXPECT_DOUBLE_EQ(MathSpecial::powauto(0.0, 0.0), 1.0); + EXPECT_DOUBLE_EQ(MathSpecial::powauto(5.0, 0.0), 1.0); + // powint() convention for a zero base, which differs from std::pow() + EXPECT_DOUBLE_EQ(MathSpecial::powauto(0.0, 3.0), 0.0); + EXPECT_DOUBLE_EQ(MathSpecial::powauto(0.0, -3.0), 0.0); +} + +TEST(MathSpecial, powauto_consistency_with_pow) +{ + const double base = 1.7; + for (int n = -5; n <= 10; n++) { + EXPECT_NEAR(MathSpecial::powauto(base, (double) n), std::pow(base, (double) n), + std::abs(std::pow(base, (double) n)) * 1e-14) + << "powauto(" << base << ", " << n << ")"; + } +} + // ------------------------------------------------------------------------- // powsinxx() // -------------------------------------------------------------------------