From 1ebd4ccb303d9340ce2c3f76563da87086d8770b Mon Sep 17 00:00:00 2001 From: Marc Descoteaux <87877404+MarcD3@users.noreply.github.com> Date: Thu, 25 Jun 2026 11:13:05 -0400 Subject: [PATCH 01/17] Update test for new naming scheme --- unittest/force-styles/tests/manybody-pair-dispersion_d3.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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: ! "" From 991ce331e0c15998f0770c4ee0231ffdc0d2d5b9 Mon Sep 17 00:00:00 2001 From: Marc Descoteaux <87877404+MarcD3@users.noreply.github.com> Date: Thu, 25 Jun 2026 11:16:06 -0400 Subject: [PATCH 02/17] Add Accelerated Dispersion Variants --- src/EXTRA-PAIR/pair_dispersion_d3.cpp | 73 +- src/EXTRA-PAIR/pair_dispersion_d3.h | 55 +- src/KOKKOS/pair_dispersion_d3_kokkos.cpp | 821 +++++++++++++++ src/KOKKOS/pair_dispersion_d3_kokkos.h | 1190 ++++++++++++++++++++++ src/OPENMP/pair_dispersion_d3_omp.cpp | 639 ++++++++++++ src/OPENMP/pair_dispersion_d3_omp.h | 55 + 6 files changed, 2802 insertions(+), 31 deletions(-) create mode 100644 src/KOKKOS/pair_dispersion_d3_kokkos.cpp create mode 100644 src/KOKKOS/pair_dispersion_d3_kokkos.h create mode 100644 src/OPENMP/pair_dispersion_d3_omp.cpp create mode 100644 src/OPENMP/pair_dispersion_d3_omp.h diff --git a/src/EXTRA-PAIR/pair_dispersion_d3.cpp b/src/EXTRA-PAIR/pair_dispersion_d3.cpp index 1c86bb074f8..5919e0e7c30 100644 --- a/src/EXTRA-PAIR/pair_dispersion_d3.cpp +++ b/src/EXTRA-PAIR/pair_dispersion_d3.cpp @@ -33,9 +33,9 @@ #include "update.h" #include +#include #include #include -#include #include #include @@ -58,8 +58,8 @@ 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 double AUTOANG = 0.52917725; // atomic units (Bohr) to Angstrom +static constexpr double AUTOEV = 27.21140795; // atomic units (Hartree) to eV #include "d3_parameters.h" @@ -82,6 +82,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 +91,8 @@ PairDispersionD3::PairDispersionD3(LAMMPS *lmp) : PairDispersionD3::~PairDispersionD3() { + if (copymode) return; + if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); @@ -250,8 +253,8 @@ 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); @@ -319,8 +322,18 @@ void PairDispersionD3::coeff(int narg, char **arg) // set r0ab read_r0ab(atomic_numbers, ntypes); - // read c6ab + // read c6ab and determine max grid size read_c6ab(atomic_numbers, ntypes); + max_mxci = 0; + for (int i = 1; i <= ntypes; i++) { + if (mxci[i] > max_mxci) max_mxci = mxci[i]; + } + if (max_mxci < 4) max_mxci = 4; + if (max_mxci > 4) { + memory->destroy(c6ab); + memory->create(c6ab, ntypes + 1, ntypes + 1, max_mxci + 1, max_mxci + 1, 3, "pair:c6ab"); + read_c6ab(atomic_numbers, ntypes); + } free(atomic_numbers); } @@ -377,7 +390,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 @@ -417,7 +430,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 * pow6(AUTOANG); if (c6_ref > 0) { cni_ref = c6ab[iat][jat][ci][cj][1]; @@ -510,7 +523,6 @@ 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; @@ -519,7 +531,7 @@ void PairDispersionD3::compute(int eflag, int vflag) double *c6_res = get_dC6(type[i], type[j], cn[i], cn[j]); 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 +542,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 = pow_general(ip6, alpha6) * pow_general(rsq, -half_alpha6); damp6 = 1.0 / (1.0 + 6.0 * t6); - t8 = pow(rs8 / r0, alpha8); + t8 = pow_general(ip8, alpha8) * pow_general(rsq, -half_alpha8); damp8 = 1.0 / (1.0 + 6.0 * t8); e6 = C6 * damp6 * r6inv; @@ -551,14 +568,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 = pow_general((r / (rs6 * r0)) + rs8 * r0, -alpha6); damp6 = 1.0 / (1.0 + 6.0 * t6); - t8 = pow((r / r0) + rs8 * r0, -alpha8); + t8 = pow_general((r / r0) + rs8 * r0, -alpha8); damp8 = 1.0 / (1.0 + 6.0 * t8); e6 = C6 * damp6 * r6inv; @@ -586,8 +604,8 @@ 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); + t6 = r6 + pow6(a1 * r0 + a2); + t8 = r8 + pow8(a1 * r0 + a2); e6 = C6 / t6; e8 = C8 / t8; @@ -607,8 +625,8 @@ 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); + t6 = r6 + pow6(a1 * r0 + a2); + t8 = r8 + pow8(a1 * r0 + a2); e6 = C6 / t6; e8 = C8 / t8; @@ -685,7 +703,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 +1057,7 @@ void PairDispersionD3::set_funcpar(std::string &functional_name) break; } - rs8 = rs8 / autoang; + rs8 = rs8 / AUTOANG; } break; case 3: { // bj @@ -1061,7 +1079,7 @@ void PairDispersionD3::set_funcpar(std::string &functional_name) {"hf/mixed", 41}, {"hf/sv", 42}, {"hf/minis", 43}, {"b3-lyp/6-31gd", 44}, {"hcth120", 45}, {"pw1pw", 46}, {"pwgga", 47}, {"hsesol", 48}, {"hf3c", 49}, {"hf3cv", 50}, {"pbeh3c", 51}, {"pbeh-3c", 52}, - {"mn15", 53}, {"r2scan", 54}}; + {"mn15", 53}}; int functionalCode = functionalMap[functional_name]; switch (functionalCode) { @@ -1356,7 +1374,7 @@ void PairDispersionD3::set_funcpar(std::string &functional_name) break; } - a2 = a2 * autoang; + a2 = a2 * AUTOANG; } break; case 4: { // bjm @@ -1417,8 +1435,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 @@ -1532,7 +1549,7 @@ 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); + bytes += (double) (n + 1) * (n + 1) * 5 * 5 * 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..d7251eec680 100644 --- a/src/EXTRA-PAIR/pair_dispersion_d3.h +++ b/src/EXTRA-PAIR/pair_dispersion_d3.h @@ -21,8 +21,56 @@ PairStyle(dispersion/d3,PairDispersionD3); #define LMP_PAIR_DISPERSION_D3_H #include "pair.h" +#include +#include namespace LAMMPS_NS { +template inline T pow2(T x) +{ return x * x; } + +template inline T pow4(T x) +{ + const T x2 = x * x; + return x2 * x2; +} + +template inline T pow6(T x) +{ return pow2(x) * pow4(x); } + +template inline T pow8(T x) +{ + const T x4 = pow4(x); + return x4 * x4; +} + +template inline T pow_int(T base, int exponent) +{ + if (exponent == 0) return static_cast(1); + + if (exponent < 0) return static_cast(1) / pow_int(base, -exponent); + + T result = static_cast(1); + while (exponent > 0) { + if (exponent & 1) result *= base; + base *= base; + exponent >>= 1; + } + return result; +} + +template inline bool is_integer_value(T x) +{ + const T nearest_int = std::round(x); + return x == nearest_int; +} + +template inline T pow_general(T base, T exponent) +{ + if (is_integer_value(exponent)) return pow_int(base, static_cast(std::round(exponent))); + + using std::pow; + return pow(base, exponent); +} class PairDispersionD3 : public Pair { @@ -48,15 +96,16 @@ 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) diff --git a/src/KOKKOS/pair_dispersion_d3_kokkos.cpp b/src/KOKKOS/pair_dispersion_d3_kokkos.cpp new file mode 100644 index 00000000000..75f013a3345 --- /dev/null +++ b/src/KOKKOS/pair_dispersion_d3_kokkos.cpp @@ -0,0 +1,821 @@ +// 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; + +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 + +#include "d3_parameters.h" + +/* ---------------------------------------------------------------------- */ + +template +PairDispersionD3Kokkos::PairDispersionD3Kokkos(LAMMPS *lmp) : PairDispersionD3(lmp) +{ + respa_enable = 0; + + nmax = 0; + comm_forward = 2; + comm_reverse = 2; + + restartinfo = 0; + one_coeff = 1; + single_enable = 0; + + dampingCode = 0; + s6 = s8 = s18 = rs6 = rs8 = rs18 = a1 = a2 = alpha = alpha6 = alpha8 = 0.0; + + kokkosable = 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 + if (newton_pair){ + Kokkos::parallel_for( + Kokkos::RangePolicy(0, nall), + PairDispersionD3InitializeFunctor{d_cn, d_dc6}); + } + else { + Kokkos::parallel_for( + Kokkos::RangePolicy(0, nlocal), + PairDispersionD3InitializeFunctor{d_cn, d_dc6}); + } + + // calculate coordination number + if (newton_pair) { + if (neighflag == FULL) { + // Initialize coordination number kernel + PairDispersionD3CoordinationNumberKernel cnkernel( + x, type, d_rcov, d_cn, d_ilist, d_numneigh, d_neighbors, nlocal, nall, cn_thr); + + // Do parallel computation + Kokkos::parallel_for(Kokkos::RangePolicy(0, inum),cnkernel); + + // Contribute to coordination number + cnkernel.contribute(); + + } else if (neighflag == HALFTHREAD) { + PairDispersionD3CoordinationNumberKernel cnkernel( + x, type, d_rcov, d_cn, d_ilist, d_numneigh, d_neighbors, nlocal, nall, cn_thr); + Kokkos::parallel_for(Kokkos::RangePolicy(0, inum),cnkernel); + cnkernel.contribute(); + } else if (neighflag == HALF) { + PairDispersionD3CoordinationNumberKernel cnkernel( + x, type, d_rcov, d_cn, d_ilist, d_numneigh, d_neighbors, nlocal, nall, cn_thr); + Kokkos::parallel_for(Kokkos::RangePolicy(0, inum),cnkernel); + cnkernel.contribute(); + } else { + error->all(FLERR, "Invalid neighflag in PairDispersionD3Kokkos"); + } + + } else { + if (neighflag == FULL) { + PairDispersionD3CoordinationNumberKernel cnkernel( + x, type, d_rcov, d_cn, d_ilist, d_numneigh, d_neighbors, nlocal, nall, cn_thr); + Kokkos::parallel_for(Kokkos::RangePolicy(0, inum),cnkernel); + cnkernel.contribute(); + } else if (neighflag == HALFTHREAD) { + PairDispersionD3CoordinationNumberKernel cnkernel( + x, type, d_rcov, d_cn, d_ilist, d_numneigh, d_neighbors, nlocal, nall, cn_thr); + Kokkos::parallel_for(Kokkos::RangePolicy(0, inum),cnkernel); + cnkernel.contribute(); + } else if (neighflag == HALF) { + PairDispersionD3CoordinationNumberKernel cnkernel( + x, type, d_rcov, d_cn, d_ilist, d_numneigh, d_neighbors, nlocal, nall, cn_thr); + Kokkos::parallel_for(Kokkos::RangePolicy(0, inum),cnkernel); + cnkernel.contribute(); + } else { + error->all(FLERR, "Invalid neighflag in PairDispersionD3Kokkos"); + } + } + + // sync to host before host-side communication to avoid dual modification + k_cn.sync_host(); + k_dc6.sync_host(); + + // communicate coordination number + communicationStage = 1; + if (newton_pair) comm->reverse_comm(this); + comm->forward_comm(this); +} + +/* ---------------------------------------------------------------------- + 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.template modify(); + return cut; +} + + +/* ---------------------------------------------------------------------- + init pair style +------------------------------------------------------------------------- */ + +template +void PairDispersionD3Kokkos::init_style() +{ + PairDispersionD3::init_style(); + + // adjust neighbor list request for KOKKOS + + neighflag = lmp->kokkos->neighflag; + 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); + k_c6ab = decltype(k_c6ab)("pair:c6ab", n+1, n+1, this->max_mxci + 1, this->max_mxci + 1, 3); + + d_r2r4 = k_r2r4.template view(); + d_rcov = k_rcov.template view(); + d_mxci = k_mxci.template view(); + d_r0ab = k_r0ab.template view(); + d_c6ab = k_c6ab.template view(); +} + + +/* ---------------------------------------------------------------------- + Coeff: read from pair_coeff + pair_coeff * * path_r0ab.csv path_c6ab.csv functional element1 element2 ... +------------------------------------------------------------------------- */ + +template +void PairDispersionD3Kokkos::coeff(int narg, char **arg) +{ + int ntypes = atom->ntypes; + if (narg != ntypes + 2) error->all(FLERR, "Pair_coeff * * needs: element1 element2 ..."); + + if (!allocated) allocate(); + std::string element; + int *atomic_numbers = (int *) malloc(sizeof(int) * ntypes); + for (int i = 0; i < ntypes; i++) { + element = arg[i + 2]; + atomic_numbers[i] = find_atomic_number(element); + if (atomic_numbers[i] < 0) + error->all(FLERR, Error::NOLASTLINE, "Element {} not supported", element); + } + + int count = 0; + for (int i = 1; i <= ntypes; i++) { + for (int j = 1; j <= ntypes; j++) { + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients" + utils::errorurl(21)); + + for (int i = 1; i <= ntypes; i++) { + r2r4[i] = r2r4_ref[atomic_numbers[i - 1]]; + rcov[i] = rcov_ref[atomic_numbers[i - 1]]; + } + + // set r0ab + read_r0ab(atomic_numbers, ntypes); + + // read c6ab and determine max grid size + read_c6ab(atomic_numbers, ntypes); + max_mxci = 0; + for (int i = 1; i <= ntypes; i++) { + if (mxci[i] > max_mxci) max_mxci = mxci[i]; + } + if (max_mxci < 4) max_mxci = 4; + if (max_mxci > 4) { + memory->destroy(c6ab); + memory->create(c6ab, ntypes + 1, ntypes + 1, max_mxci + 1, max_mxci + 1, 3, "pair:c6ab"); + read_c6ab(atomic_numbers, ntypes); + } + + const int desired_c6ab = max_mxci + 1; + if (k_c6ab.extent_int(2) != desired_c6ab || k_c6ab.extent_int(3) != desired_c6ab) { + k_c6ab = decltype(k_c6ab)("pair:c6ab", ntypes + 1, ntypes + 1, desired_c6ab, desired_c6ab, 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 <= this->max_mxci; ci++) { + for (int cj = 0; cj <= this->max_mxci; 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(); + + free(atomic_numbers); +} + + +/* ---------------------------------------------------------------------- + Compute : energy, force, and stress (Required) +------------------------------------------------------------------------- */ + + +template +void PairDispersionD3Kokkos::compute(int eflag_in, int vflag_in) +{ + eflag = eflag_in; + vflag = vflag_in; + + if (neighflag == FULL) no_virial_fdotr_compute = 1; + + 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; + + + if constexpr (std::is_same_v) { + need_dup = false; + } else 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, "Invalid neighflag in PairDispersionD3Kokkos"); + } + + 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, "Invalid neighflag in PairDispersionD3Kokkos"); + } + + 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, f, type, + d_cutsq, d_cn, d_dc6, + d_r2r4, d_r0ab, d_rcov, 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, nall, + eflag, vflag_either, + eflag_global, eflag_atom, + vflag_global, vflag_atom, + dampingCode, + s6, s8, rs6, rs8, + a1, a2, alpha, + 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_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, f, type, + d_cutsq, d_dc6, d_rcov, + d_numneigh, d_neighbors, d_ilist, + nall, + 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..9a98c31ca0e --- /dev/null +++ b/src/KOKKOS/pair_dispersion_d3_kokkos.h @@ -0,0 +1,1190 @@ +/* -*- 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 "pair_dispersion_d3.h" +#include "pair_kokkos.h" + +// 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. +*/ + +static constexpr double AUTOANG = 0.52917725; // atomic units (Bohr) to Angstrom +static constexpr double AUTOEV = 27.21140795; // atomic units (Hartree) to eV + + +namespace LAMMPS_NS { + +template +KOKKOS_INLINE_FUNCTION +T kk_pow2(T x) +{ + return x * x; +} + +template +KOKKOS_INLINE_FUNCTION +T kk_pow4(T x) +{ + const T x2 = x * x; + return x2 * x2; +} + +template +KOKKOS_INLINE_FUNCTION +T kk_pow6(T x) +{ + return kk_pow2(x) * kk_pow4(x); +} + +template +KOKKOS_INLINE_FUNCTION +T kk_pow8(T x) +{ + const T x4 = kk_pow4(x); + return x4 * x4; +} + +template +KOKKOS_INLINE_FUNCTION +T kk_pow_int(T base, int exponent) +{ + if (exponent == 0) + return static_cast(1); + + if (exponent < 0) + return static_cast(1) / kk_pow_int(base, -exponent); + + T result = static_cast(1); + while (exponent > 0) + { + if (exponent & 1) + result *= base; + base *= base; + exponent >>= 1; + } + return result; +} + +template +KOKKOS_INLINE_FUNCTION +bool kk_is_integer(T x) +{ + const T nearest_int = Kokkos::round(x); + return x == nearest_int; +} + +template +KOKKOS_INLINE_FUNCTION +T kk_pow_general(T base, T exponent) +{ + if (kk_is_integer(exponent)) + return kk_pow_int(base, static_cast(Kokkos::round(exponent))); + + return Kokkos::pow(base, exponent); +} + +/* ---------------------------------------------------------------------- */ +// 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 = std::conditional_t, + Kokkos::Experimental::ScatterNonDuplicated, + 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; + int nall; + 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, int nall_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), nall(nall_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 + 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 = std::conditional_t, + Kokkos::Experimental::ScatterNonDuplicated, + 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_kkacc_1d_3 f; + 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; + typename AT::t_kkfloat_1d d_rcov; + 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 nall; + 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; + KK_FLOAT cn_thr; + + PairDispersionD3KernelA( + const typename AT::t_kkfloat_1d_3_lr_randomread &x_in, + const typename AT::t_kkacc_1d_3 &f_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 typename AT::t_kkfloat_1d &d_rcov_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 nall_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, KK_FLOAT cn_thr_in) + : x(x_in), f(f_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_rcov(d_rcov_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), nall(nall_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), 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; + } + + // Init per-thread reduction variable + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void init(value_type &ev) const { + ev.evdwl = 0.0; + for (int i = 0; i < 6; i++) ev.v[i] = 0.0; + } + + // Combine per-thread reduction variables + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void join(value_type &dst, const value_type &src) const { + dst.evdwl += src.evdwl; + for (int i = 0; i < 6; i++) dst.v[i] += src.v[i]; + } + + /* ---------------------------------------------------------------------- + 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; + + const KK_FLOAT autoang6 = kk_pow6(AUTOANG); + + 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 = 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; + } + } + } + + if (den > 1.0E-99) { + 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 = kk_pow_general(ip6, alpha6) * kk_pow_general(rsq, -half_alpha6); + t8 = kk_pow_general(ip8, alpha8) * kk_pow_general(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 = kk_pow_general((r / (rs6 * r0)) + rs8 * r0, -alpha6); + damp6 = 1.0f / (1.0f + 6.0f * t6); + t8 = kk_pow_general((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 + + 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; + + t6 = r6 + kk_pow6(a1 * r0 + a2); + t8 = r8 + kk_pow8(a1 * r0 + a2); + + 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 + + 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; + + t6 = r6 + kk_pow6(a1 * r0 + a2); + t8 = r8 + kk_pow8(a1 * r0 + a2); + + 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; + + default: { + Kokkos::abort("PairDispersionD3Kokkos: invalid dampingCode"); + } break; + } + + 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 = std::conditional_t, + Kokkos::Experimental::ScatterNonDuplicated, + 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_kkacc_1d_3 f; + 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; + int nall; + + 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_kkacc_1d_3 &f_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, + int nall_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), f(f_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), nall(nall_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; + } + + // Init per-thread reduction variable + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void init(value_type &ev) const { + ev.evdwl = 0.0; + for (int i = 0; i < 6; i++) ev.v[i] = 0.0; + } + + // Combine thread reduction variables + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void join(value_type &dst, const value_type &src) const { + dst.evdwl += src.evdwl; + for (int i = 0; i < 6; i++) dst.v[i] += src.v[i]; + } + + + // 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 = 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(); + void compute(int, int) override; + double init_one(int, int) override; + void init_style() override; + void allocate(); + 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; + // NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + static int sbmask(const int &j) { + return j >> SBBITS & 3; + } + + + 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 + + typename AT::tdual_kkfloat_2d k_cutsq; // cutoff^2 table (dual view) + 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_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..59fae8d364b --- /dev/null +++ b/src/OPENMP/pair_dispersion_d3_omp.cpp @@ -0,0 +1,639 @@ +// 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 "domain.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "neigh_list.h" +#include "neighbor.h" +#include "suffix.h" +#include "update.h" // needed for Update::ntimestep + +#include "omp_compat.h" +#include +#include +#include // for setenv +#include +#include + + +using namespace LAMMPS_NS; + +// global ad hoc parameters - copied from pair_dispersion_d3.cpp +static constexpr double K1 = 16.0; +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 + +/* ---------------------------------------------------------------------- */ + +PairDispersionD3OMP::PairDispersionD3OMP(LAMMPS *lmp) : + PairDispersionD3(lmp), ThrOMP(lmp, THR_PAIR) +{ + suffix_flag |= Suffix::OMP; +} + + + + +void PairDispersionD3OMP::calc_coordination_number() +{ + int nlocal = atom->nlocal; + const int nthreads = comm->nthreads; + int nall = nlocal + atom->nghost; + + + int newton_pair = force->newton_pair; + + if (atom->nmax > nmax) { + nmax = atom->nmax; + memory->grow(cn, nmax, "pair:cn"); + memory->grow(dc6, nmax, "pair:dc6"); + } + + // zero out coordination number + memset(cn, 0, sizeof(double) * (newton_pair ? nall : nlocal)); + memset(dc6, 0, sizeof(double) * (newton_pair ? nall : nlocal)); + + int inum = list->inum; + + // Begin parallel region, the central atoms indexed by ii are assigned to different threads. + #if defined(_OPENMP) + #pragma omp parallel default(none) \ + shared(force,atom,list,fix,cn,dc6,cutsq,r0ab,r2r4,rcov,s6,s8,rs6,rs8,alpha,a1,a2,dampingCode,cn_thr,eatom,vatom,stderr,comm) \ + 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 * _noalias const ilist = list->ilist; + const int * _noalias const numneigh = list->numneigh; + const int * const * const firstneigh = list->firstneigh; + + // Thread-local cn array to avoid race conditions + double *thr_cn = new double[atom->nmax](); // Initialize to zero + + 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; } + } + } + + // Contribute thread-local cn to global cn + for (int i = 0; i < atom->nmax; i++) { + // No need for atomic if newton_pair is false + #pragma omp atomic + cn[i] += thr_cn[i]; + } + + delete[] thr_cn; + + +} + + +void PairDispersionD3OMP::compute(int eflag, int vflag) +{ + ev_init(eflag, vflag); + // 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; + + // Zero out dc6 values before OpenMP section + memset(dc6, 0, sizeof(double) * nall); + +// Parallel direct force computation and some other quantities calculation. +#if defined(_OPENMP) +#pragma omp parallel default(none) \ +shared(eflag,vflag,force,atom,list,fix,cn,dc6,cutsq,r0ab,r2r4,rcov,s6,s8,rs6,rs8,alpha,a1,a2,dampingCode,cn_thr,eatom,vatom,stderr,comm) \ +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 + + + // 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 default(none) \ + shared(eflag,vflag,force,atom,list,fix,cn,dc6,cutsq,r0ab,r2r4,rcov,s6,s8,rs6,rs8,alpha,a1,a2,dampingCode,cn_thr,eatom,vatom,stderr,comm) \ + 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_either, thr); + } //end of omp parallel region + + + if (vflag_fdotr) virial_fdotr_compute(); +} + +// Modified from original code to avoid race conditions +void PairDispersionD3OMP::get_dC6(int iat, int jat, double cni, double cnj, double c6_res[3]) +{ + + + double c6_ref, cni_ref, cnj_ref; + double c6mem, r_save, r; + double expterm, term; + double num, den, d_num_i, d_num_j, d_den_i, d_den_j; + + c6mem = -1.0e20, r_save = 1.0e20; + num = 0; + den = 0; + d_num_i = 0; + d_num_j = 0; + d_den_i = 0; + d_den_j = 0; + + for (int ci = 0; ci <= mxci[iat]; ci++) { + for (int cj = 0; cj <= mxci[jat]; cj++) { + + c6_ref = c6ab[iat][jat][ci][cj][0]; + c6_ref *= AUTOEV * pow6(AUTOANG); + + if (c6_ref > 0) { + cni_ref = c6ab[iat][jat][ci][cj][1]; + cnj_ref = 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 = exp(static_cast(K3) * static_cast(r)); + + num += c6_ref * expterm; + den += expterm; + + expterm = expterm * 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; + } + } + } + + if (den > 1.0E-99) { + c6_res[0] = num / den; + c6_res[1] = ((d_num_i * den) - (d_den_i * num)) / (den * den); + c6_res[2] = ((d_num_j * den) - (d_den_j * num)) / (den * den); + } else { + c6_res[0] = c6mem; + c6_res[1] = 0; + c6_res[2] = 0; + } + +} + + +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; + double evdwl = 0.0; + + // Thread-local dc6 array to avoid race conditions + double *thr_dc6 = new double[atom->nmax](); // Initialize to zero. Consider if we can dodge this. + + // 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; + + // Modified from original code to avoid race conditions + 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 = pow_general(ip6, alpha6) * pow_general(rsq, -half_alpha6); + t8 = pow_general(ip8, alpha8) * pow_general(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 = pow_general((r / (rs6 * r0)) + rs8 * r0, -alpha6); + damp6 = 1.0 / (1.0 + 6.0 * t6); + t8 = pow_general((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 + double r0 = sqrt(C8 / C6); + + double r4 = rsq * rsq; + double r6 = rsq * rsq * rsq; + double r8 = rsq * rsq * rsq * rsq; + + t6 = r6 + pow6(a1 * r0 + a2); + t8 = r8 + pow8(a1 * r0 + a2); + + 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 r0 = sqrt(C8 / C6); + + double r4 = rsq * rsq; + double r6 = rsq * rsq * rsq; + double r8 = rsq * rsq * rsq * rsq; + + t6 = r6 + pow6(a1 * r0 + a2); + t8 = r8 + pow8(a1 * r0 + a2); + + 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; + + default: { + error->all(FLERR, Error::NOLASTLINE, "Damping code {} unknown", dampingCode); + } 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); + } + } + } + + // Contribute thread-local dc6 to global dc6 + for (int i = 0; i < atom->nmax; i++) { + // Possibly no need for atomic if newton_pair is false + #pragma omp atomic + dc6[i] += thr_dc6[i]; + } + + delete[] thr_dc6; +} + +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 dc6tmp,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(); + 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..c874b822142 --- /dev/null +++ b/src/OPENMP/pair_dispersion_d3_omp.h @@ -0,0 +1,55 @@ +/* -*- 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(); + + void get_dC6(int iat, int jat, double cni, double cnj, double c6_res[3]); + + template + void eval_second_phase(int ifrom, int ito, ThrData *const thr); + + int thr_evflag; // Flag for ev_tally thread usage +}; + +} // namespace LAMMPS_NS + +#endif +#endif From fe058ca659b41da37497de591adf26befe3c9c6c Mon Sep 17 00:00:00 2001 From: Marc Descoteaux <87877404+MarcD3@users.noreply.github.com> Date: Thu, 25 Jun 2026 11:21:53 -0400 Subject: [PATCH 03/17] replace missing r2scan case --- src/EXTRA-PAIR/pair_dispersion_d3.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/EXTRA-PAIR/pair_dispersion_d3.cpp b/src/EXTRA-PAIR/pair_dispersion_d3.cpp index 5919e0e7c30..31a631a89cd 100644 --- a/src/EXTRA-PAIR/pair_dispersion_d3.cpp +++ b/src/EXTRA-PAIR/pair_dispersion_d3.cpp @@ -1079,7 +1079,7 @@ void PairDispersionD3::set_funcpar(std::string &functional_name) {"hf/mixed", 41}, {"hf/sv", 42}, {"hf/minis", 43}, {"b3-lyp/6-31gd", 44}, {"hcth120", 45}, {"pw1pw", 46}, {"pwgga", 47}, {"hsesol", 48}, {"hf3c", 49}, {"hf3cv", 50}, {"pbeh3c", 51}, {"pbeh-3c", 52}, - {"mn15", 53}}; + {"mn15", 53}, {"r2scan", 54}}; int functionalCode = functionalMap[functional_name]; switch (functionalCode) { From 6bd43a684ae19c7743992874b35116ef7d6c12f0 Mon Sep 17 00:00:00 2001 From: Marc Descoteaux <87877404+MarcD3@users.noreply.github.com> Date: Thu, 25 Jun 2026 12:04:53 -0400 Subject: [PATCH 04/17] update documentation for accelerated variants --- doc/src/pair_dispersion_d3.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/doc/src/pair_dispersion_d3.rst b/doc/src/pair_dispersion_d3.rst index eec12157f26..e303d2496d1 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 """""" From b42efb556ac5301b83f66d0e048ea8db79d245f3 Mon Sep 17 00:00:00 2001 From: Marc Descoteaux <87877404+MarcD3@users.noreply.github.com> Date: Thu, 25 Jun 2026 15:08:25 -0400 Subject: [PATCH 05/17] apply whitespace corrections --- src/KOKKOS/pair_dispersion_d3_kokkos.cpp | 12 ++--- src/KOKKOS/pair_dispersion_d3_kokkos.h | 6 +-- src/OPENMP/pair_dispersion_d3_omp.cpp | 68 ++++++++++++------------ 3 files changed, 43 insertions(+), 43 deletions(-) diff --git a/src/KOKKOS/pair_dispersion_d3_kokkos.cpp b/src/KOKKOS/pair_dispersion_d3_kokkos.cpp index 75f013a3345..0addda6356f 100644 --- a/src/KOKKOS/pair_dispersion_d3_kokkos.cpp +++ b/src/KOKKOS/pair_dispersion_d3_kokkos.cpp @@ -14,8 +14,8 @@ /* ---------------------------------------------------------------------- Contributing authors: - Yizhong R. Hu - Marc L. Descoteaux + Yizhong R. Hu + Marc L. Descoteaux Ulrik Unneberg William C. Witt Affiliation: Harvard University @@ -49,7 +49,7 @@ template PairDispersionD3Kokkos::PairDispersionD3Kokkos(LAMMPS *lmp) : PairDispersionD3(lmp) { respa_enable = 0; - + nmax = 0; comm_forward = 2; comm_reverse = 2; @@ -371,9 +371,9 @@ void PairDispersionD3Kokkos::compute(int eflag_in, int vflag_in) k_cutsq.template sync(); k_mxci.template sync(); k_c6ab.template sync(); - if (eflag || vflag) + if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else + else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); @@ -393,7 +393,7 @@ void PairDispersionD3Kokkos::compute(int eflag_in, int vflag_in) d_neighbors = k_list->d_neighbors; d_ilist = k_list->d_ilist; inum = list->inum; - + if constexpr (std::is_same_v) { need_dup = false; diff --git a/src/KOKKOS/pair_dispersion_d3_kokkos.h b/src/KOKKOS/pair_dispersion_d3_kokkos.h index 9a98c31ca0e..1e2454ef0a0 100644 --- a/src/KOKKOS/pair_dispersion_d3_kokkos.h +++ b/src/KOKKOS/pair_dispersion_d3_kokkos.h @@ -655,7 +655,7 @@ struct PairDispersionD3KernelA { t6 = kk_pow_general(ip6, alpha6) * kk_pow_general(rsq, -half_alpha6); t8 = kk_pow_general(ip8, alpha8) * kk_pow_general(rsq, -half_alpha8); - + damp6 = 1.0f / (1.0f + 6.0f * t6); damp8 = 1.0f / (1.0f + 6.0f * t8); @@ -893,7 +893,7 @@ struct PairDispersionD3KernelB { dst.evdwl += src.evdwl; for (int i = 0; i < 6; i++) dst.v[i] += src.v[i]; } - + // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION @@ -1079,7 +1079,7 @@ class PairDispersionD3Kokkos : public PairDispersionD3, public KokkosBase { typedef ArrayTypes HAT; typedef EV_FLOAT value_type; PairDispersionD3Kokkos(class LAMMPS *); - ~PairDispersionD3Kokkos() override; + ~PairDispersionD3Kokkos() override; void calc_coordination_number(); void compute(int, int) override; diff --git a/src/OPENMP/pair_dispersion_d3_omp.cpp b/src/OPENMP/pair_dispersion_d3_omp.cpp index 59fae8d364b..9142abc768d 100644 --- a/src/OPENMP/pair_dispersion_d3_omp.cpp +++ b/src/OPENMP/pair_dispersion_d3_omp.cpp @@ -64,8 +64,8 @@ void PairDispersionD3OMP::calc_coordination_number() { int nlocal = atom->nlocal; const int nthreads = comm->nthreads; - int nall = nlocal + atom->nghost; - + int nall = nlocal + atom->nghost; + int newton_pair = force->newton_pair; @@ -107,7 +107,7 @@ void PairDispersionD3OMP::calc_coordination_number() thr->timer(Timer::PAIR); } - + // communicate coordination number communicationStage = 1; if (newton_pair) comm->reverse_comm(this); @@ -118,7 +118,7 @@ template void PairDispersionD3OMP::eval_coordination(int iifrom, int iito, ThrData * const thr) { - const auto * _noalias const x = (dbl3_t *) atom->x[0]; + const auto * _noalias const x = (dbl3_t *) atom->x[0]; const int * _noalias const type = atom->type; const int nlocal = atom->nlocal; const int * _noalias const ilist = list->ilist; @@ -160,7 +160,7 @@ void PairDispersionD3OMP::eval_coordination(int iifrom, int iito, ThrData * cons if (NEWTON_PAIR || j < nlocal) { thr_cn[j] += cn_ij; } } } - + // Contribute thread-local cn to global cn for (int i = 0; i < atom->nmax; i++) { // No need for atomic if newton_pair is false @@ -221,19 +221,19 @@ firstprivate(inum,nthreads,nall) } thr->timer(Timer::PAIR); } // end of omp parallel region - + // 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 default(none) \ shared(eflag,vflag,force,atom,list,fix,cn,dc6,cutsq,r0ab,r2r4,rcov,s6,s8,rs6,rs8,alpha,a1,a2,dampingCode,cn_thr,eatom,vatom,stderr,comm) \ @@ -269,7 +269,7 @@ firstprivate(inum,nthreads,nall) reduce_thr(this, eflag, vflag_either, thr); } //end of omp parallel region - + if (vflag_fdotr) virial_fdotr_compute(); } @@ -355,7 +355,7 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const // Thread-local dc6 array to avoid race conditions double *thr_dc6 = new double[atom->nmax](); // Initialize to zero. Consider if we can dodge this. - + // Loop over assigned atoms for (int ii = iifrom; ii < iito; ++ii) { int i = ilist[ii]; @@ -368,48 +368,48 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const 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; // Modified from original code to avoid race conditions - double c6_res[3] = {}; + 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 */ + // 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]]; @@ -419,7 +419,7 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const t6 = pow_general(ip6, alpha6) * pow_general(rsq, -half_alpha6); t8 = pow_general(ip8, alpha8) * pow_general(rsq, -half_alpha8); - + damp6 = 1.0 / (1.0 + 6.0 * t6); damp8 = 1.0 / (1.0 + 6.0 * t8); @@ -435,7 +435,7 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const fpair = fpair1 + fpair2; fpair *= factor_lj; } break; - // Written to avoid pow + // Written to avoid pow case 2: { // zerom double r = sqrt(rsq); @@ -502,27 +502,27 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const fpair = -(tmp6 + tmp8); fpair *= factor_lj; } break; - + default: { error->all(FLERR, Error::NOLASTLINE, "Damping code {} unknown", dampingCode); } 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) { + + if (NEWTON_PAIR || j < nlocal) { double dc6_contrib_j = rest * c6_res[2]; - thr_dc6[j] += dc6_contrib_j; + thr_dc6[j] += dc6_contrib_j; } - + f[i].x += delx * fpair; f[i].y += dely * fpair; f[i].z += delz * fpair; @@ -552,7 +552,7 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const 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; @@ -573,7 +573,7 @@ void PairDispersionD3OMP::eval_second_phase(int iifrom, int iito, ThrData * cons 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; @@ -616,7 +616,7 @@ void PairDispersionD3OMP::eval_second_phase(int iifrom, int iito, ThrData * cons } // 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); + if (EVFLAG) ev_tally_thr(this, i, j, nlocal, NEWTON_PAIR, 0.0, 0.0, fpair, delx, dely, delz, thr); } } From 46577252c350944b346f08180c1c0cee7f1a3c3d Mon Sep 17 00:00:00 2001 From: Marc Descoteaux <87877404+MarcD3@users.noreply.github.com> Date: Thu, 25 Jun 2026 15:20:27 -0400 Subject: [PATCH 06/17] use the LMP macros for omp compatibility --- src/OPENMP/pair_dispersion_d3_omp.cpp | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/OPENMP/pair_dispersion_d3_omp.cpp b/src/OPENMP/pair_dispersion_d3_omp.cpp index 9142abc768d..c781162e650 100644 --- a/src/OPENMP/pair_dispersion_d3_omp.cpp +++ b/src/OPENMP/pair_dispersion_d3_omp.cpp @@ -83,8 +83,7 @@ void PairDispersionD3OMP::calc_coordination_number() // Begin parallel region, the central atoms indexed by ii are assigned to different threads. #if defined(_OPENMP) - #pragma omp parallel default(none) \ - shared(force,atom,list,fix,cn,dc6,cutsq,r0ab,r2r4,rcov,s6,s8,rs6,rs8,alpha,a1,a2,dampingCode,cn_thr,eatom,vatom,stderr,comm) \ + #pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED() \ firstprivate(inum,nthreads) #endif { @@ -189,8 +188,7 @@ void PairDispersionD3OMP::compute(int eflag, int vflag) // Parallel direct force computation and some other quantities calculation. #if defined(_OPENMP) -#pragma omp parallel default(none) \ -shared(eflag,vflag,force,atom,list,fix,cn,dc6,cutsq,r0ab,r2r4,rcov,s6,s8,rs6,rs8,alpha,a1,a2,dampingCode,cn_thr,eatom,vatom,stderr,comm) \ +#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag) \ firstprivate(inum,nthreads,nall) #endif { @@ -235,8 +233,7 @@ firstprivate(inum,nthreads,nall) // Process the second phase with the combined dc6 values #if defined(_OPENMP) - #pragma omp parallel default(none) \ - shared(eflag,vflag,force,atom,list,fix,cn,dc6,cutsq,r0ab,r2r4,rcov,s6,s8,rs6,rs8,alpha,a1,a2,dampingCode,cn_thr,eatom,vatom,stderr,comm) \ + #pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag) \ firstprivate(inum,nthreads,nall) #endif { @@ -636,4 +633,3 @@ double PairDispersionD3OMP::memory_usage() bytes += PairDispersionD3::memory_usage(); return bytes; } - From cdc4be7ae3cffaf0f45c41c104de434fc827a5bd Mon Sep 17 00:00:00 2001 From: Marc Descoteaux <87877404+MarcD3@users.noreply.github.com> Date: Thu, 9 Jul 2026 18:36:19 -0400 Subject: [PATCH 07/17] improve whitespace and formatting --- src/OPENMP/pair_dispersion_d3_omp.cpp | 54 ++++++++------------------- 1 file changed, 16 insertions(+), 38 deletions(-) diff --git a/src/OPENMP/pair_dispersion_d3_omp.cpp b/src/OPENMP/pair_dispersion_d3_omp.cpp index c781162e650..5f2c89b836f 100644 --- a/src/OPENMP/pair_dispersion_d3_omp.cpp +++ b/src/OPENMP/pair_dispersion_d3_omp.cpp @@ -18,7 +18,6 @@ Affiliation: Harvard University ------------------------------------------------------------------------- */ - #include "pair_dispersion_d3_omp.h" #include "atom.h" @@ -32,13 +31,13 @@ #include "suffix.h" #include "update.h" // needed for Update::ntimestep -#include "omp_compat.h" #include #include #include // for setenv #include #include +#include "omp_compat.h" using namespace LAMMPS_NS; @@ -57,8 +56,7 @@ PairDispersionD3OMP::PairDispersionD3OMP(LAMMPS *lmp) : suffix_flag |= Suffix::OMP; } - - +/* ---------------------------------------------------------------------- */ void PairDispersionD3OMP::calc_coordination_number() { @@ -66,7 +64,6 @@ void PairDispersionD3OMP::calc_coordination_number() const int nthreads = comm->nthreads; int nall = nlocal + atom->nghost; - int newton_pair = force->newton_pair; if (atom->nmax > nmax) { @@ -102,9 +99,7 @@ void PairDispersionD3OMP::calc_coordination_number() eval_coordination<0>(ifrom,ito,thr); } - thr->timer(Timer::PAIR); - } // communicate coordination number @@ -113,6 +108,8 @@ void PairDispersionD3OMP::calc_coordination_number() comm->forward_comm(this); } +/* ---------------------------------------------------------------------- */ + template void PairDispersionD3OMP::eval_coordination(int iifrom, int iito, ThrData * const thr) { @@ -166,12 +163,10 @@ void PairDispersionD3OMP::eval_coordination(int iifrom, int iito, ThrData * cons #pragma omp atomic cn[i] += thr_cn[i]; } - delete[] thr_cn; - - } +/* ---------------------------------------------------------------------- */ void PairDispersionD3OMP::compute(int eflag, int vflag) { @@ -220,7 +215,6 @@ firstprivate(inum,nthreads,nall) thr->timer(Timer::PAIR); } // end of omp parallel region - // Communication stage 2 for dc6 values in preparation for calculation of indirect forces in the second phase communicationStage = 2; if (force->newton_pair) { @@ -229,9 +223,7 @@ firstprivate(inum,nthreads,nall) 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) @@ -243,7 +235,6 @@ firstprivate(inum,nthreads,nall) 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) { @@ -266,15 +257,15 @@ firstprivate(inum,nthreads,nall) reduce_thr(this, eflag, vflag_either, thr); } //end of omp parallel region - if (vflag_fdotr) virial_fdotr_compute(); } -// Modified from original code to avoid race conditions +/* ---------------------------------------------------------------------- + Modified from serial code to avoid race conditions + ---------------------------------------------------------------------- */ + void PairDispersionD3OMP::get_dC6(int iat, int jat, double cni, double cnj, double c6_res[3]) { - - double c6_ref, cni_ref, cnj_ref; double c6mem, r_save, r; double expterm, term; @@ -332,14 +323,13 @@ void PairDispersionD3OMP::get_dC6(int iat, int jat, double cni, double cnj, doub c6_res[1] = 0; c6_res[2] = 0; } - } +/* ---------------------------------------------------------------------- */ 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; @@ -357,7 +347,6 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const 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; @@ -365,13 +354,11 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const 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; @@ -382,17 +369,15 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const 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; // Modified from original code to avoid race conditions - double c6_res[3] = {}; + 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; @@ -405,7 +390,7 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const // Damping code selection - now using the passed dampingCode parameter switch (dampingCode) { - // Written to avoid using sqrt and pow() + // Written to avoid using sqrt and pow() case 1: /* Original damping */ { double ip6 = rs6 * r0ab[type[i]][type[j]]; @@ -432,7 +417,7 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const fpair = fpair1 + fpair2; fpair *= factor_lj; } break; - // Written to avoid pow + // Written to avoid pow case 2: { // zerom double r = sqrt(rsq); @@ -518,8 +503,6 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const thr_dc6[j] += dc6_contrib_j; } - - f[i].x += delx * fpair; f[i].y += dely * fpair; f[i].z += delz * fpair; @@ -546,10 +529,11 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const delete[] thr_dc6; } +/* ---------------------------------------------------------------------- */ + 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; @@ -561,8 +545,6 @@ void PairDispersionD3OMP::eval_second_phase(int iifrom, int iito, ThrData * cons double dc6tmp,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]; @@ -590,7 +572,6 @@ void PairDispersionD3OMP::eval_second_phase(int iifrom, int iito, ThrData * cons 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)); @@ -613,18 +594,15 @@ void PairDispersionD3OMP::eval_second_phase(int iifrom, int iito, ThrData * cons } // 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); - + 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() From 80e916b773b6f49f615a225620f41e6ef6244e22 Mon Sep 17 00:00:00 2001 From: Marc Descoteaux <87877404+MarcD3@users.noreply.github.com> Date: Thu, 9 Jul 2026 18:44:34 -0400 Subject: [PATCH 08/17] use vector --- src/OPENMP/pair_dispersion_d3_omp.cpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/OPENMP/pair_dispersion_d3_omp.cpp b/src/OPENMP/pair_dispersion_d3_omp.cpp index 5f2c89b836f..4ba1cb176bf 100644 --- a/src/OPENMP/pair_dispersion_d3_omp.cpp +++ b/src/OPENMP/pair_dispersion_d3_omp.cpp @@ -36,6 +36,7 @@ #include // for setenv #include #include +#include #include "omp_compat.h" @@ -122,7 +123,7 @@ void PairDispersionD3OMP::eval_coordination(int iifrom, int iito, ThrData * cons const int * const * const firstneigh = list->firstneigh; // Thread-local cn array to avoid race conditions - double *thr_cn = new double[atom->nmax](); // Initialize to zero + auto thr_cn = std::vector(atom->nmax); // Initialize to zero for (int ii = iifrom; ii < iito; ii++) { @@ -163,7 +164,6 @@ void PairDispersionD3OMP::eval_coordination(int iifrom, int iito, ThrData * cons #pragma omp atomic cn[i] += thr_cn[i]; } - delete[] thr_cn; } /* ---------------------------------------------------------------------- */ @@ -341,7 +341,7 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const double evdwl = 0.0; // Thread-local dc6 array to avoid race conditions - double *thr_dc6 = new double[atom->nmax](); // Initialize to zero. Consider if we can dodge this. + auto thr_dc6 = std::vector(atom->nmax); // Initialize to zero. // Loop over assigned atoms for (int ii = iifrom; ii < iito; ++ii) { @@ -525,8 +525,6 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const #pragma omp atomic dc6[i] += thr_dc6[i]; } - - delete[] thr_dc6; } /* ---------------------------------------------------------------------- */ From 1cc39699acfc500bde00821e4a0e3f2854e3db68 Mon Sep 17 00:00:00 2001 From: Marc Descoteaux <87877404+MarcD3@users.noreply.github.com> Date: Thu, 9 Jul 2026 18:52:54 -0400 Subject: [PATCH 09/17] work on formatting and whitespace --- src/KOKKOS/pair_dispersion_d3_kokkos.cpp | 28 +++++++++++++++--------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/src/KOKKOS/pair_dispersion_d3_kokkos.cpp b/src/KOKKOS/pair_dispersion_d3_kokkos.cpp index 0addda6356f..65ef49ccb72 100644 --- a/src/KOKKOS/pair_dispersion_d3_kokkos.cpp +++ b/src/KOKKOS/pair_dispersion_d3_kokkos.cpp @@ -68,7 +68,6 @@ PairDispersionD3Kokkos::PairDispersionD3Kokkos(LAMMPS *lmp) : PairDi datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK; } - /* ---------------------------------------------------------------------- */ template @@ -190,7 +189,6 @@ double PairDispersionD3Kokkos::init_one(int i, int j) return cut; } - /* ---------------------------------------------------------------------- init pair style ------------------------------------------------------------------------- */ @@ -210,7 +208,6 @@ void PairDispersionD3Kokkos::init_style() if (neighflag == FULL) request->enable_full(); } - /* ---------------------------------------------------------------------- allocate req. arrays ------------------------------------------------------------------------- */ @@ -238,7 +235,6 @@ void PairDispersionD3Kokkos::allocate() d_c6ab = k_c6ab.template view(); } - /* ---------------------------------------------------------------------- Coeff: read from pair_coeff pair_coeff * * path_r0ab.csv path_c6ab.csv functional element1 element2 ... @@ -338,12 +334,10 @@ void PairDispersionD3Kokkos::coeff(int narg, char **arg) free(atomic_numbers); } - /* ---------------------------------------------------------------------- Compute : energy, force, and stress (Required) ------------------------------------------------------------------------- */ - template void PairDispersionD3Kokkos::compute(int eflag_in, int vflag_in) { @@ -387,14 +381,12 @@ void PairDispersionD3Kokkos::compute(int eflag_in, int vflag_in) 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; - if constexpr (std::is_same_v) { need_dup = false; } else if (neighflag == FULL) { @@ -569,6 +561,8 @@ void PairDispersionD3Kokkos::launch_kernel_A(EV_FLOAT &ev) } } +/* ---------------------------------------------------------------------- */ + template template void PairDispersionD3Kokkos::dispatch_kernel_A(EV_FLOAT &ev) @@ -582,6 +576,8 @@ void PairDispersionD3Kokkos::dispatch_kernel_A(EV_FLOAT &ev) } } +/* ---------------------------------------------------------------------- */ + template template void PairDispersionD3Kokkos::launch_kernel_B(EV_FLOAT &ev) @@ -613,6 +609,8 @@ void PairDispersionD3Kokkos::launch_kernel_B(EV_FLOAT &ev) } } +/* ---------------------------------------------------------------------- */ + template template void PairDispersionD3Kokkos::dispatch_kernel_B(EV_FLOAT &ev) @@ -626,12 +624,10 @@ void PairDispersionD3Kokkos::dispatch_kernel_B(EV_FLOAT &ev) } } - /* ---------------------------------------------------------------------- Communication section ------------------------------------------------------------------------- */ - template int PairDispersionD3Kokkos::pack_forward_comm_kokkos(int n, DAT::tdual_int_1d k_sendlist, DAT::tdual_double_1d &buf, @@ -654,6 +650,8 @@ int PairDispersionD3Kokkos::pack_forward_comm_kokkos(int n, DAT::tdu return n; } +/* ---------------------------------------------------------------------- */ + template void PairDispersionD3Kokkos::unpack_forward_comm_kokkos(int n, int first, DAT::tdual_double_1d &buf) { @@ -674,6 +672,8 @@ void PairDispersionD3Kokkos::unpack_forward_comm_kokkos(int n, int f } } +/* ---------------------------------------------------------------------- */ + template int PairDispersionD3Kokkos::pack_reverse_comm_kokkos(int n, int first, DAT::tdual_double_1d &buf) { @@ -693,6 +693,8 @@ int PairDispersionD3Kokkos::pack_reverse_comm_kokkos(int n, int firs return n; } +/* ---------------------------------------------------------------------- */ + template void PairDispersionD3Kokkos::unpack_reverse_comm_kokkos(int n, DAT::tdual_int_1d k_recvlist, DAT::tdual_double_1d &buf) { @@ -745,6 +747,8 @@ int PairDispersionD3Kokkos::pack_forward_comm(int n, int *list, doub return m; } +/* ---------------------------------------------------------------------- */ + template void PairDispersionD3Kokkos::unpack_forward_comm(int n, int first, double *buf) { @@ -769,6 +773,8 @@ void PairDispersionD3Kokkos::unpack_forward_comm(int n, int first, d } } +/* ---------------------------------------------------------------------- */ + template int PairDispersionD3Kokkos::pack_reverse_comm(int n, int first, double *buf) { @@ -787,6 +793,8 @@ int PairDispersionD3Kokkos::pack_reverse_comm(int n, int first, doub return m; } +/* ---------------------------------------------------------------------- */ + template void PairDispersionD3Kokkos::unpack_reverse_comm(int n, int *list, double *buf) { From 6bc529d1cb94dda5611049b68e9f2c974a7e8bbe Mon Sep 17 00:00:00 2001 From: Marc Descoteaux <87877404+MarcD3@users.noreply.github.com> Date: Thu, 9 Jul 2026 18:53:54 -0400 Subject: [PATCH 10/17] comment box format --- src/OPENMP/pair_dispersion_d3_omp.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OPENMP/pair_dispersion_d3_omp.cpp b/src/OPENMP/pair_dispersion_d3_omp.cpp index 4ba1cb176bf..49a5f20a58f 100644 --- a/src/OPENMP/pair_dispersion_d3_omp.cpp +++ b/src/OPENMP/pair_dispersion_d3_omp.cpp @@ -262,7 +262,7 @@ firstprivate(inum,nthreads,nall) /* ---------------------------------------------------------------------- Modified from serial code to avoid race conditions - ---------------------------------------------------------------------- */ +------------------------------------------------------------------------- */ void PairDispersionD3OMP::get_dC6(int iat, int jat, double cni, double cnj, double c6_res[3]) { From 648f37fc7e907f1065037e91e953df92ef101bd4 Mon Sep 17 00:00:00 2001 From: Marc Descoteaux <87877404+MarcD3@users.noreply.github.com> Date: Thu, 9 Jul 2026 19:24:03 -0400 Subject: [PATCH 11/17] remove empty omp shared clause --- src/OPENMP/pair_dispersion_d3_omp.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OPENMP/pair_dispersion_d3_omp.cpp b/src/OPENMP/pair_dispersion_d3_omp.cpp index 49a5f20a58f..3181c52b748 100644 --- a/src/OPENMP/pair_dispersion_d3_omp.cpp +++ b/src/OPENMP/pair_dispersion_d3_omp.cpp @@ -81,7 +81,7 @@ void PairDispersionD3OMP::calc_coordination_number() // Begin parallel region, the central atoms indexed by ii are assigned to different threads. #if defined(_OPENMP) - #pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED() \ + #pragma omp parallel LMP_DEFAULT_NONE \ firstprivate(inum,nthreads) #endif { From f80d8a7df4c16402d5506ede21e2e00d9f7170dc Mon Sep 17 00:00:00 2001 From: Marc Descoteaux <87877404+MarcD3@users.noreply.github.com> Date: Mon, 20 Jul 2026 20:07:38 -0600 Subject: [PATCH 12/17] Add powauto for cases where the integer-nature of the exponent may not be known --- src/KOKKOS/math_special_kokkos.h | 15 +++++++++++++++ src/math_special.h | 11 +++++++++++ 2 files changed, 26 insertions(+) diff --git a/src/KOKKOS/math_special_kokkos.h b/src/KOKKOS/math_special_kokkos.h index 5ecb23a891e..9e84d0f9e36 100644 --- a/src/KOKKOS/math_special_kokkos.h +++ b/src/KOKKOS/math_special_kokkos.h @@ -234,6 +234,21 @@ namespace LAMMPS_NS::MathSpecialKokkos { return (n > 0) ? yy : static_cast(1) / yy; } + template +// NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + static T powauto(const T x, const T n) + { + if (x == static_cast(0)) return static_cast(0); + const T rounded = round(n); + + if (n == rounded) { + return powint(x, static_cast(rounded)); + } + + return 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/math_special.h b/src/math_special.h index 3cce45538fa..032aad2c105 100644 --- a/src/math_special.h +++ b/src/math_special.h @@ -164,6 +164,17 @@ namespace LAMMPS_NS::MathSpecial { return (n > 0) ? yy : 1.0 / yy; } + static inline double powauto(const double x, const double n) + { + if (n == 0.0) return 1.0; + if (x == 0.0) return 0.0; // Use MathSpecial powint convention (different from std::pow convention) + 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. From b454b455685cd0762d4ebf511a34f5d35a82fe5a Mon Sep 17 00:00:00 2001 From: Marc Descoteaux <87877404+MarcD3@users.noreply.github.com> Date: Mon, 20 Jul 2026 20:07:50 -0600 Subject: [PATCH 13/17] simplify pow()-family usage --- src/EXTRA-PAIR/pair_dispersion_d3.cpp | 29 +++++--- src/EXTRA-PAIR/pair_dispersion_d3.h | 47 ------------- src/KOKKOS/pair_dispersion_d3_kokkos.h | 97 +++++--------------------- src/OPENMP/pair_dispersion_d3_omp.cpp | 29 +++++--- 4 files changed, 59 insertions(+), 143 deletions(-) diff --git a/src/EXTRA-PAIR/pair_dispersion_d3.cpp b/src/EXTRA-PAIR/pair_dispersion_d3.cpp index 31a631a89cd..dbc445f6182 100644 --- a/src/EXTRA-PAIR/pair_dispersion_d3.cpp +++ b/src/EXTRA-PAIR/pair_dispersion_d3.cpp @@ -27,6 +27,7 @@ #include "error.h" #include "force.h" #include "info.h" +#include "math_special.h" #include "memory.h" #include "neigh_list.h" #include "neighbor.h" @@ -430,7 +431,9 @@ 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 * pow6(AUTOANG); + double autoang6 = AUTOANG * AUTOANG * AUTOANG; + autoang6 = MathSpecial::square(autoang6); + c6_ref *= AUTOEV * autoang6; if (c6_ref > 0) { cni_ref = c6ab[iat][jat][ci][cj][1]; @@ -551,9 +554,9 @@ void PairDispersionD3::compute(int eflag, int vflag) double half_alpha6 = 0.5 * alpha6; double half_alpha8 = 0.5 * alpha8; - t6 = pow_general(ip6, alpha6) * pow_general(rsq, -half_alpha6); + t6 = MathSpecial::powauto(ip6, alpha6) * MathSpecial::powauto(rsq, -half_alpha6); damp6 = 1.0 / (1.0 + 6.0 * t6); - t8 = pow_general(ip8, alpha8) * pow_general(rsq, -half_alpha8); + t8 = MathSpecial::powauto(ip8, alpha8) * MathSpecial::powauto(rsq, -half_alpha8); damp8 = 1.0 / (1.0 + 6.0 * t8); e6 = C6 * damp6 * r6inv; @@ -574,9 +577,9 @@ void PairDispersionD3::compute(int eflag, int vflag) double r = sqrt(rsq); double r0 = r0ab[type[i]][type[j]]; - t6 = pow_general((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_general((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; @@ -604,8 +607,12 @@ void PairDispersionD3::compute(int eflag, int vflag) double r6 = rsq * rsq * rsq; double r8 = rsq * rsq * rsq * rsq; - t6 = r6 + pow6(a1 * r0 + a2); - t8 = r8 + pow8(a1 * r0 + a2); + 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; @@ -625,8 +632,12 @@ void PairDispersionD3::compute(int eflag, int vflag) double r6 = rsq * rsq * rsq; double r8 = rsq * rsq * rsq * rsq; - t6 = r6 + pow6(a1 * r0 + a2); - t8 = r8 + pow8(a1 * r0 + a2); + 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; diff --git a/src/EXTRA-PAIR/pair_dispersion_d3.h b/src/EXTRA-PAIR/pair_dispersion_d3.h index d7251eec680..f9c375060fb 100644 --- a/src/EXTRA-PAIR/pair_dispersion_d3.h +++ b/src/EXTRA-PAIR/pair_dispersion_d3.h @@ -25,53 +25,6 @@ PairStyle(dispersion/d3,PairDispersionD3); #include namespace LAMMPS_NS { -template inline T pow2(T x) -{ return x * x; } - -template inline T pow4(T x) -{ - const T x2 = x * x; - return x2 * x2; -} - -template inline T pow6(T x) -{ return pow2(x) * pow4(x); } - -template inline T pow8(T x) -{ - const T x4 = pow4(x); - return x4 * x4; -} - -template inline T pow_int(T base, int exponent) -{ - if (exponent == 0) return static_cast(1); - - if (exponent < 0) return static_cast(1) / pow_int(base, -exponent); - - T result = static_cast(1); - while (exponent > 0) { - if (exponent & 1) result *= base; - base *= base; - exponent >>= 1; - } - return result; -} - -template inline bool is_integer_value(T x) -{ - const T nearest_int = std::round(x); - return x == nearest_int; -} - -template inline T pow_general(T base, T exponent) -{ - if (is_integer_value(exponent)) return pow_int(base, static_cast(std::round(exponent))); - - using std::pow; - return pow(base, exponent); -} - class PairDispersionD3 : public Pair { public: diff --git a/src/KOKKOS/pair_dispersion_d3_kokkos.h b/src/KOKKOS/pair_dispersion_d3_kokkos.h index 1e2454ef0a0..421f6c0c778 100644 --- a/src/KOKKOS/pair_dispersion_d3_kokkos.h +++ b/src/KOKKOS/pair_dispersion_d3_kokkos.h @@ -24,6 +24,7 @@ PairStyle(dispersion/d3/kk/host,PairDispersionD3Kokkos); #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" @@ -46,75 +47,6 @@ static constexpr double AUTOEV = 27.21140795; // atomic units (Hartree) to eV namespace LAMMPS_NS { -template -KOKKOS_INLINE_FUNCTION -T kk_pow2(T x) -{ - return x * x; -} - -template -KOKKOS_INLINE_FUNCTION -T kk_pow4(T x) -{ - const T x2 = x * x; - return x2 * x2; -} - -template -KOKKOS_INLINE_FUNCTION -T kk_pow6(T x) -{ - return kk_pow2(x) * kk_pow4(x); -} - -template -KOKKOS_INLINE_FUNCTION -T kk_pow8(T x) -{ - const T x4 = kk_pow4(x); - return x4 * x4; -} - -template -KOKKOS_INLINE_FUNCTION -T kk_pow_int(T base, int exponent) -{ - if (exponent == 0) - return static_cast(1); - - if (exponent < 0) - return static_cast(1) / kk_pow_int(base, -exponent); - - T result = static_cast(1); - while (exponent > 0) - { - if (exponent & 1) - result *= base; - base *= base; - exponent >>= 1; - } - return result; -} - -template -KOKKOS_INLINE_FUNCTION -bool kk_is_integer(T x) -{ - const T nearest_int = Kokkos::round(x); - return x == nearest_int; -} - -template -KOKKOS_INLINE_FUNCTION -T kk_pow_general(T base, T exponent) -{ - if (kk_is_integer(exponent)) - return kk_pow_int(base, static_cast(Kokkos::round(exponent))); - - return Kokkos::pow(base, exponent); -} - /* ---------------------------------------------------------------------- */ // Functor to initialize cn and dc6 arrays /* ---------------------------------------------------------------------- */ @@ -438,7 +370,8 @@ struct PairDispersionD3KernelA { d_den_i = 0.0; d_den_j = 0.0; - const KK_FLOAT autoang6 = kk_pow6(AUTOANG); + KK_FLOAT autoang6 = AUTOANG * AUTOANG * AUTOANG; + autoang6 = MathSpecialKokkos::square(autoang6); int maxci = d_mxci(iat); int maxcj = d_mxci(jat); @@ -653,8 +586,8 @@ struct PairDispersionD3KernelA { const KK_FLOAT half_alpha6 = 0.5 * alpha6; const KK_FLOAT half_alpha8 = 0.5 * alpha8; - t6 = kk_pow_general(ip6, alpha6) * kk_pow_general(rsq, -half_alpha6); - t8 = kk_pow_general(ip8, alpha8) * kk_pow_general(rsq, -half_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); @@ -678,9 +611,9 @@ struct PairDispersionD3KernelA { const KK_FLOAT r0 = d_r0ab(itype, jtype); const KK_FLOAT r = Kokkos::sqrt(rsq); - t6 = kk_pow_general((r / (rs6 * r0)) + rs8 * r0, -alpha6); + t6 = MathSpecialKokkos::powauto((r / (rs6 * r0)) + rs8 * r0, -alpha6); damp6 = 1.0f / (1.0f + 6.0f * t6); - t8 = kk_pow_general((r / r0) + rs8 * r0, -alpha8); + t8 = MathSpecialKokkos::powauto((r / r0) + rs8 * r0, -alpha8); damp8 = 1.0f / (1.0f + 6.0f * t8); e6 = C6 * damp6 * r6inv; @@ -708,8 +641,12 @@ struct PairDispersionD3KernelA { KK_FLOAT r6 = rsq * rsq * rsq; KK_FLOAT r8 = r6 * rsq; - t6 = r6 + kk_pow6(a1 * r0 + a2); - t8 = r8 + kk_pow8(a1 * r0 + a2); + 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; @@ -729,8 +666,12 @@ struct PairDispersionD3KernelA { KK_FLOAT r6 = rsq * rsq * rsq; KK_FLOAT r8 = r6 * rsq; - t6 = r6 + kk_pow6(a1 * r0 + a2); - t8 = r8 + kk_pow8(a1 * r0 + a2); + 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; diff --git a/src/OPENMP/pair_dispersion_d3_omp.cpp b/src/OPENMP/pair_dispersion_d3_omp.cpp index 3181c52b748..a6294ea3c98 100644 --- a/src/OPENMP/pair_dispersion_d3_omp.cpp +++ b/src/OPENMP/pair_dispersion_d3_omp.cpp @@ -25,6 +25,7 @@ #include "domain.h" #include "error.h" #include "force.h" +#include "math_special.h" #include "memory.h" #include "neigh_list.h" #include "neighbor.h" @@ -283,7 +284,9 @@ void PairDispersionD3OMP::get_dC6(int iat, int jat, double cni, double cnj, doub for (int cj = 0; cj <= mxci[jat]; cj++) { c6_ref = c6ab[iat][jat][ci][cj][0]; - c6_ref *= AUTOEV * pow6(AUTOANG); + double autoang6 = AUTOANG * AUTOANG * AUTOANG; + autoang6 = MathSpecial::square(autoang6); + c6_ref *= AUTOEV * autoang6; if (c6_ref > 0) { cni_ref = c6ab[iat][jat][ci][cj][1]; @@ -399,8 +402,8 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const double half_alpha6 = 0.5 * alpha6; double half_alpha8 = 0.5 * alpha8; - t6 = pow_general(ip6, alpha6) * pow_general(rsq, -half_alpha6); - t8 = pow_general(ip8, alpha8) * pow_general(rsq, -half_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); @@ -423,9 +426,9 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const double r = sqrt(rsq); double r0 = r0ab[type[i]][type[j]]; - t6 = pow_general((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_general((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; @@ -452,8 +455,12 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const double r6 = rsq * rsq * rsq; double r8 = rsq * rsq * rsq * rsq; - t6 = r6 + pow6(a1 * r0 + a2); - t8 = r8 + pow8(a1 * r0 + a2); + 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; @@ -472,8 +479,12 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const double r6 = rsq * rsq * rsq; double r8 = rsq * rsq * rsq * rsq; - t6 = r6 + pow6(a1 * r0 + a2); - t8 = r8 + pow8(a1 * r0 + a2); + 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; From 754211e0376679e9326f53724fc3f3b748ce75fd Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 25 Aug 2026 16:00:20 +0000 Subject: [PATCH 14/17] fix CI blockers for dispersion/d3 accelerator variants - remove trailing whitespace in MathSpecial::powauto() - register pair_dispersion_d3_kokkos.{cpp,h} in src/KOKKOS/Install.sh so the traditional make build installs the style (the base style lives in the optional EXTRA-PAIR package, hence the dependency argument) - add the accelerator code letters to the dispersion/d3 entry in doc/src/Commands_pair.rst - add the missing ".. include:: accel_styles.rst" block to doc/src/pair_dispersion_d3.rst, which doc/utils/check-styles.py requires once a page lists accelerated style index entries --- doc/src/Commands_pair.rst | 2 +- doc/src/pair_dispersion_d3.rst | 5 +++++ src/KOKKOS/Install.sh | 2 ++ src/math_special.h | 2 +- 4 files changed, 9 insertions(+), 2 deletions(-) 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/pair_dispersion_d3.rst b/doc/src/pair_dispersion_d3.rst index e303d2496d1..1e8aa449e2f 100644 --- a/doc/src/pair_dispersion_d3.rst +++ b/doc/src/pair_dispersion_d3.rst @@ -109,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/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/math_special.h b/src/math_special.h index 032aad2c105..2f3f35cf3f3 100644 --- a/src/math_special.h +++ b/src/math_special.h @@ -170,7 +170,7 @@ namespace LAMMPS_NS::MathSpecial { if (x == 0.0) return 0.0; // Use MathSpecial powint convention (different from std::pow convention) const double rounded = std::round(n); if (n == rounded) { - return powint(x, static_cast(rounded)); + return powint(x, static_cast(rounded)); } return std::pow(x, n); } From 759fcaf5a56f97710b1a1bfcd2214e8abc1ff8ac Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 25 Aug 2026 16:33:51 +0000 Subject: [PATCH 15/17] clean up, document and test MathSpecial::powauto() The KOKKOS version was missing the "n == 0" shortcut, so powauto(0, 0) returned 0 on the device and 1 on the host, and it called round() and pow() unqualified, which resolves to the global double overloads in device code and promotes for a KK_FLOAT argument. Use Kokkos::round() and Kokkos::pow() as required for device kernels, and take the arguments by const reference like the neighboring powint(). Add the doxygen comment block that every other function in these headers has, the matching entry in the programmer guide, and unit tests covering the integer, fractional and zero cases against powint() and std::pow(). --- doc/src/Developer_utils.rst | 3 +++ src/KOKKOS/math_special_kokkos.h | 24 +++++++++++------ src/math_special.h | 18 ++++++++++--- unittest/utils/test_math_special.cpp | 40 ++++++++++++++++++++++++++++ 4 files changed, 73 insertions(+), 12 deletions(-) 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/src/KOKKOS/math_special_kokkos.h b/src/KOKKOS/math_special_kokkos.h index 9e84d0f9e36..98885338d20 100644 --- a/src/KOKKOS/math_special_kokkos.h +++ b/src/KOKKOS/math_special_kokkos.h @@ -234,19 +234,27 @@ 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) + 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 = round(n); - - if (n == rounded) { - return powint(x, static_cast(rounded)); - } - - return pow(x, n); + 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 diff --git a/src/math_special.h b/src/math_special.h index 2f3f35cf3f3..2296a2cbf58 100644 --- a/src/math_special.h +++ b/src/math_special.h @@ -164,14 +164,24 @@ 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; // Use MathSpecial powint convention (different from std::pow convention) + if (x == 0.0) return 0.0; const double rounded = std::round(n); - if (n == rounded) { - return powint(x, static_cast(rounded)); - } + if (n == rounded) return powint(x, static_cast(rounded)); return std::pow(x, n); } 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() // ------------------------------------------------------------------------- From 4bf355bb677801783f0d2858da2846ffedf6d246 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 25 Aug 2026 16:34:14 +0000 Subject: [PATCH 16/17] dispersion/d3: fix the /omp virial and share code between the variants The /omp variant tallied the virial incorrectly. compute() passed vflag_either to reduce_thr(), which bit-tests its argument against the raw VIRIAL_* flags, while ev_setup() reduces vflag_either to 0/1. So "1 & (VIRIAL_ATOM|VIRIAL_CENTROID)" was always 0 and the per-atom virial was never reduced across threads, and "1 & (VIRIAL_PAIR|VIRIAL_FDOTR)" was 1, so the per-thread fdotr virial was added on top of the trailing virial_fdotr_compute() call. With "compute stress/atom" active the global pressure came out 3.4x too large and the per-atom stress lost everything computed by threads 1..N-1. With only the fdotr virial requested the two errors cancelled, which is why the force styles test did not catch it. Pass vflag and drop the trailing virial_fdotr_compute(), as in every other OPENMP pair style. Energy, pressure tensor and summed per-atom stress are now bit identical to the plain style for 1, 2, 4 and 8 threads. The thread private cn and dc6 accumulators were std::vector of atom->nmax, allocated and zeroed per thread on every timestep and reduced with an "omp atomic" loop over all of nmax, which also read the elements past nall that no memset had initialized. Switch to the scheme pair_adp/omp uses: grow cn and dc6 to nthreads*nmax, let each thread accumulate into its own nall long slice, and reduce with data_reduce_thr(). PairDispersionD3::get_dC6() returned a pointer to a function local "static double c6_res[3]", so it was not usable from a threaded or a device loop. Both accelerated variants worked around that by carrying a verbatim ~60 line copy of it. Give the base version an output parameter instead and delete both copies, so a future fix to the interpolation cannot miss a variant. The KOKKOS variant also carried a copy of coeff() with the device transfer appended, which pulled the 2.9 MB d3_parameters.h - and with it a second copy of the 32385 x 5 reference C6 table - into a second translation unit. Call the base version and keep only the device transfer. The constants K1, K3, AUTOANG and AUTOEV existed in three places, and in pair_dispersion_d3_kokkos.h they sat at global namespace scope in a style header, which is what akohlmey objected to for the pow() helpers. Define them once in namespace LAMMPS_NS::DispersionD3 in the base header, together with AUTOANG6, which all three variants recomputed inside the innermost loop of get_dC6(). Declare allocate() and calc_coordination_number() virtual in the base header and mark the overrides, as required for a KOKKOS style; they were hidden, not overridden, and only happened to work because coeff() and compute() are overridden as well. The "unknown damping code" check moves to PairDispersionD3::init_style(), so neither the threaded loop nor the device kernel has to report an error from a context that cannot. Further fixes found while reviewing: - coeff() reallocated c6ab and re-read the reference table whenever the largest CN grid index exceeded 4. The largest encoded index in d3_parameters.h is 4, so the branch was unreachable - and had it ever fired, the first read_c6ab() call would already have written out of bounds on the fixed 5x5 buffer. Drop the branch, keep the max_mxci value the KOKKOS view sizing needs, and bounds check the grid indices in read_c6ab(), where they are produced. - coeff() leaked the malloc()ed element list when an unknown element name triggered error->all(); use std::vector. - k_cutsq was a plain dual view of KK_FLOAT while cutsq is double **, which only compiles when KK_FLOAT is double. Use DAT::ttransform_kkfloat_2d, as in pair_lj_cut_kokkos. - The kernel init()/join() reduction hooks zeroed and summed only evdwl and v[6] of EV_FLOAT, leaving ecoul indeterminate. Drop them and use the default constructor and operator+= of s_EV_FLOAT. - The "den > 1.0E-99" fallback threshold in the device version of get_dC6() underflows to zero in a single precision build; use the smallest normalized KK_FLOAT instead. - reverse_comm_device was never set, so the device pack/unpack methods were dead code and every newton step round-tripped cn and dc6 through the host. - exp() was called unqualified in three device kernels. - Drop the unused f, d_rcov, nall and cn_thr members shipped with the kernel functors, the unused sbmask() of the pair class, the dead grid dimension constants, the unused includes, the unused thr_evflag member, and the constructor assignments the base class already makes. - Collapse the six way copy of the coordination number kernel launch into a template, and the byte identical bj/bjm damping cases into one. --- src/EXTRA-PAIR/pair_dispersion_d3.cpp | 99 ++++------- src/EXTRA-PAIR/pair_dispersion_d3.h | 37 +++- src/KOKKOS/pair_dispersion_d3_kokkos.cpp | 205 +++++++---------------- src/KOKKOS/pair_dispersion_d3_kokkos.h | 149 ++++------------ src/OPENMP/pair_dispersion_d3_omp.cpp | 192 ++++++--------------- src/OPENMP/pair_dispersion_d3_omp.h | 6 +- 6 files changed, 208 insertions(+), 480 deletions(-) diff --git a/src/EXTRA-PAIR/pair_dispersion_d3.cpp b/src/EXTRA-PAIR/pair_dispersion_d3.cpp index dbc445f6182..911b94799c1 100644 --- a/src/EXTRA-PAIR/pair_dispersion_d3.cpp +++ b/src/EXTRA-PAIR/pair_dispersion_d3.cpp @@ -41,26 +41,12 @@ #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" @@ -130,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"); } /* ---------------------------------------------------------------------- @@ -259,6 +245,10 @@ void PairDispersionD3::read_c6ab(int *atomic_numbers, int ntypes) 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; @@ -297,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); @@ -321,22 +311,14 @@ void PairDispersionD3::coeff(int narg, char **arg) } // set r0ab - read_r0ab(atomic_numbers, ntypes); + read_r0ab(atomic_numbers.data(), ntypes); - // read c6ab and determine max grid size - read_c6ab(atomic_numbers, ntypes); - max_mxci = 0; - for (int i = 1; i <= ntypes; i++) { - if (mxci[i] > max_mxci) max_mxci = mxci[i]; - } - if (max_mxci < 4) max_mxci = 4; - if (max_mxci > 4) { - memory->destroy(c6ab); - memory->create(c6ab, ntypes + 1, ntypes + 1, max_mxci + 1, max_mxci + 1, 3, "pair: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]); } /* ---------------------------------------------------------------------- @@ -410,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; @@ -431,9 +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]; - double autoang6 = AUTOANG * AUTOANG * AUTOANG; - autoang6 = MathSpecial::square(autoang6); - c6_ref *= AUTOEV * autoang6; + c6_ref *= AUTOEV * AUTOANG6; if (c6_ref > 0) { cni_ref = c6ab[iat][jat][ci][cj][1]; @@ -473,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; } /* ---------------------------------------------------------------------- @@ -531,7 +508,9 @@ void PairDispersionD3::compute(int eflag, int vflag) 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; @@ -599,32 +578,8 @@ void PairDispersionD3::compute(int eflag, int vflag) fpair *= factor_lj; } break; - case 3: { // bj - - 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; - - case 4: { // bjm + case 3: // bj + case 4: { // bjm, same functional form as bj, different parameters double r0 = sqrt(C8 / C6); @@ -1475,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); } @@ -1559,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 f9c375060fb..e9298345935 100644 --- a/src/EXTRA-PAIR/pair_dispersion_d3.h +++ b/src/EXTRA-PAIR/pair_dispersion_d3.h @@ -21,8 +21,33 @@ PairStyle(dispersion/d3,PairDispersionD3); #define LMP_PAIR_DISPERSION_D3_H #include "pair.h" -#include -#include + +// 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 { @@ -65,10 +90,10 @@ class PairDispersionD3 : public Pair { 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); @@ -77,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/pair_dispersion_d3_kokkos.cpp b/src/KOKKOS/pair_dispersion_d3_kokkos.cpp index 65ef49ccb72..c99cb4e295c 100644 --- a/src/KOKKOS/pair_dispersion_d3_kokkos.cpp +++ b/src/KOKKOS/pair_dispersion_d3_kokkos.cpp @@ -37,31 +37,17 @@ using namespace LAMMPS_NS; -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 - -#include "d3_parameters.h" - /* ---------------------------------------------------------------------- */ template PairDispersionD3Kokkos::PairDispersionD3Kokkos(LAMMPS *lmp) : PairDispersionD3(lmp) { - respa_enable = 0; - - nmax = 0; - comm_forward = 2; - comm_reverse = 2; - - restartinfo = 0; - one_coeff = 1; - single_enable = 0; + kokkosable = 1; - dampingCode = 0; - s6 = s8 = s18 = rs6 = rs8 = rs18 = a1 = a2 = alpha = alpha6 = alpha8 = 0.0; + // 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; - kokkosable = 1; atomKK = (AtomKokkos *) atom; execution_space = ExecutionSpaceFromDevice::space; datamask_read = X_MASK | F_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK; @@ -106,76 +92,49 @@ void PairDispersionD3Kokkos::calc_coordination_number() k_cn.template modify(); k_dc6.template modify(); - // zero out coordination number - if (newton_pair){ - Kokkos::parallel_for( - Kokkos::RangePolicy(0, nall), - PairDispersionD3InitializeFunctor{d_cn, d_dc6}); - } - else { - Kokkos::parallel_for( - Kokkos::RangePolicy(0, nlocal), - PairDispersionD3InitializeFunctor{d_cn, d_dc6}); - } + // 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 (newton_pair) { - if (neighflag == FULL) { - // Initialize coordination number kernel - PairDispersionD3CoordinationNumberKernel cnkernel( - x, type, d_rcov, d_cn, d_ilist, d_numneigh, d_neighbors, nlocal, nall, cn_thr); - - // Do parallel computation - Kokkos::parallel_for(Kokkos::RangePolicy(0, inum),cnkernel); - - // Contribute to coordination number - cnkernel.contribute(); - - } else if (neighflag == HALFTHREAD) { - PairDispersionD3CoordinationNumberKernel cnkernel( - x, type, d_rcov, d_cn, d_ilist, d_numneigh, d_neighbors, nlocal, nall, cn_thr); - Kokkos::parallel_for(Kokkos::RangePolicy(0, inum),cnkernel); - cnkernel.contribute(); - } else if (neighflag == HALF) { - PairDispersionD3CoordinationNumberKernel cnkernel( - x, type, d_rcov, d_cn, d_ilist, d_numneigh, d_neighbors, nlocal, nall, cn_thr); - Kokkos::parallel_for(Kokkos::RangePolicy(0, inum),cnkernel); - cnkernel.contribute(); - } else { - error->all(FLERR, "Invalid neighflag in PairDispersionD3Kokkos"); - } + if (neighflag == FULL) { + dispatch_coordination_kernel(); + } else if (neighflag == HALFTHREAD) { + dispatch_coordination_kernel(); + } else if (neighflag == HALF) { + dispatch_coordination_kernel(); } else { - if (neighflag == FULL) { - PairDispersionD3CoordinationNumberKernel cnkernel( - x, type, d_rcov, d_cn, d_ilist, d_numneigh, d_neighbors, nlocal, nall, cn_thr); - Kokkos::parallel_for(Kokkos::RangePolicy(0, inum),cnkernel); - cnkernel.contribute(); - } else if (neighflag == HALFTHREAD) { - PairDispersionD3CoordinationNumberKernel cnkernel( - x, type, d_rcov, d_cn, d_ilist, d_numneigh, d_neighbors, nlocal, nall, cn_thr); - Kokkos::parallel_for(Kokkos::RangePolicy(0, inum),cnkernel); - cnkernel.contribute(); - } else if (neighflag == HALF) { - PairDispersionD3CoordinationNumberKernel cnkernel( - x, type, d_rcov, d_cn, d_ilist, d_numneigh, d_neighbors, nlocal, nall, cn_thr); - Kokkos::parallel_for(Kokkos::RangePolicy(0, inum),cnkernel); - cnkernel.contribute(); - } else { - error->all(FLERR, "Invalid neighflag in PairDispersionD3Kokkos"); - } + error->all(FLERR, "Must use half or full neighbor list style with pair dispersion/d3/kk"); } - // sync to host before host-side communication to avoid dual modification - k_cn.sync_host(); - k_dc6.sync_host(); - // 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 ------------------------------------------------------------------------- */ @@ -185,7 +144,7 @@ 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.template modify(); + k_cutsq.modify_host(); return cut; } @@ -201,6 +160,12 @@ void PairDispersionD3Kokkos::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); @@ -226,13 +191,13 @@ void PairDispersionD3Kokkos::allocate() 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); - k_c6ab = decltype(k_c6ab)("pair:c6ab", n+1, n+1, this->max_mxci + 1, this->max_mxci + 1, 3); d_r2r4 = k_r2r4.template view(); d_rcov = k_rcov.template view(); d_mxci = k_mxci.template view(); d_r0ab = k_r0ab.template view(); - d_c6ab = k_c6ab.template view(); + + // k_c6ab is created in coeff(), where the reference CN grid size is known } /* ---------------------------------------------------------------------- @@ -243,57 +208,21 @@ void PairDispersionD3Kokkos::allocate() template void PairDispersionD3Kokkos::coeff(int narg, char **arg) { - int ntypes = atom->ntypes; - if (narg != ntypes + 2) error->all(FLERR, "Pair_coeff * * needs: element1 element2 ..."); + // the base class parses the arguments and fills the host side tables + PairDispersionD3::coeff(narg, arg); - if (!allocated) allocate(); - std::string element; - int *atomic_numbers = (int *) malloc(sizeof(int) * ntypes); - for (int i = 0; i < ntypes; i++) { - element = arg[i + 2]; - atomic_numbers[i] = find_atomic_number(element); - if (atomic_numbers[i] < 0) - error->all(FLERR, Error::NOLASTLINE, "Element {} not supported", element); - } + const int ntypes = atom->ntypes; - int count = 0; - for (int i = 1; i <= ntypes; i++) { - for (int j = 1; j <= ntypes; j++) { - setflag[i][j] = 1; - count++; - } - } + // now that max_mxci is known, resize the C6 table to the grid actually used - if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients" + utils::errorurl(21)); - - for (int i = 1; i <= ntypes; i++) { - r2r4[i] = r2r4_ref[atomic_numbers[i - 1]]; - rcov[i] = rcov_ref[atomic_numbers[i - 1]]; - } - - // set r0ab - read_r0ab(atomic_numbers, ntypes); - - // read c6ab and determine max grid size - read_c6ab(atomic_numbers, ntypes); - max_mxci = 0; - for (int i = 1; i <= ntypes; i++) { - if (mxci[i] > max_mxci) max_mxci = mxci[i]; - } - if (max_mxci < 4) max_mxci = 4; - if (max_mxci > 4) { - memory->destroy(c6ab); - memory->create(c6ab, ntypes + 1, ntypes + 1, max_mxci + 1, max_mxci + 1, 3, "pair:c6ab"); - read_c6ab(atomic_numbers, ntypes); - } - - const int desired_c6ab = max_mxci + 1; - if (k_c6ab.extent_int(2) != desired_c6ab || k_c6ab.extent_int(3) != desired_c6ab) { - k_c6ab = decltype(k_c6ab)("pair:c6ab", ntypes + 1, ntypes + 1, desired_c6ab, desired_c6ab, 3); + 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(); @@ -309,13 +238,9 @@ void PairDispersionD3Kokkos::coeff(int narg, char **arg) 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 <= this->max_mxci; ci++) { - for (int cj = 0; cj <= this->max_mxci; cj++) { - for (int k = 0; k < 3; k++) { - h_c6ab(i, j, ci, cj, k) = c6ab[i][j][ci][cj][k]; - } - } - } + 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]; } } @@ -330,8 +255,6 @@ void PairDispersionD3Kokkos::coeff(int narg, char **arg) k_mxci.template sync(); k_r0ab.template sync(); k_c6ab.template sync(); - - free(atomic_numbers); } /* ---------------------------------------------------------------------- @@ -344,8 +267,6 @@ void PairDispersionD3Kokkos::compute(int eflag_in, int vflag_in) eflag = eflag_in; vflag = vflag_in; - if (neighflag == FULL) no_virial_fdotr_compute = 1; - ev_init(eflag,vflag,0); // reallocate per-atom arrays if necessary @@ -444,7 +365,7 @@ void PairDispersionD3Kokkos::compute(int eflag_in, int vflag_in) } else if (neighflag == FULL) { dispatch_kernel_A(ev); } else { - error->all(FLERR, "Invalid neighflag in PairDispersionD3Kokkos"); + error->all(FLERR, "Must use half or full neighbor list style with pair dispersion/d3/kk"); } if (evflag) ev_all += ev; @@ -473,7 +394,7 @@ void PairDispersionD3Kokkos::compute(int eflag_in, int vflag_in) } else if (neighflag == FULL) { dispatch_kernel_B(ev); } else { - error->all(FLERR, "Invalid neighflag in PairDispersionD3Kokkos"); + error->all(FLERR, "Must use half or full neighbor list style with pair dispersion/d3/kk"); } if (evflag) ev_all += ev; @@ -531,23 +452,22 @@ template void PairDispersionD3Kokkos::launch_kernel_A(EV_FLOAT &ev) { auto functor = PairDispersionD3KernelA( - x, f, type, + x, type, d_cutsq, d_cn, d_dc6, - d_r2r4, d_r0ab, d_rcov, d_c6ab, d_mxci, + 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, nall, + nlocal, eflag, vflag_either, eflag_global, eflag_atom, vflag_global, vflag_atom, dampingCode, s6, s8, rs6, rs8, - a1, a2, alpha, - cn_thr); + a1, a2, alpha); if constexpr (EVFLAG) { Kokkos::parallel_reduce( @@ -583,10 +503,9 @@ template void PairDispersionD3Kokkos::launch_kernel_B(EV_FLOAT &ev) { auto functor = PairDispersionD3KernelB( - x, f, type, + x, type, d_cutsq, d_dc6, d_rcov, d_numneigh, d_neighbors, d_ilist, - nall, dup_f, ndup_f, dup_eatom, ndup_eatom, dup_vatom, ndup_vatom, diff --git a/src/KOKKOS/pair_dispersion_d3_kokkos.h b/src/KOKKOS/pair_dispersion_d3_kokkos.h index 421f6c0c778..ba9fd4c6c55 100644 --- a/src/KOKKOS/pair_dispersion_d3_kokkos.h +++ b/src/KOKKOS/pair_dispersion_d3_kokkos.h @@ -28,25 +28,16 @@ PairStyle(dispersion/d3/kk/host,PairDispersionD3Kokkos); #include "pair_dispersion_d3.h" #include "pair_kokkos.h" -// 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. -*/ - -static constexpr double AUTOANG = 0.52917725; // atomic units (Bohr) to Angstrom -static constexpr double AUTOEV = 27.21140795; // atomic units (Hartree) to eV - +#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 /* ---------------------------------------------------------------------- */ @@ -147,7 +138,6 @@ struct PairDispersionD3CoordinationNumberKernel { typename AT::t_int_1d d_numneigh; typename AT::t_neighbors_2d d_neighbors; int nlocal; - int nall; KK_FLOAT cn_thr; KKScatterView(d_cn); } @@ -202,7 +192,7 @@ struct PairDispersionD3CoordinationNumberKernel { 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 + exp(-K1 * ((rcov_ij / rr) - 1.0))); + 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; @@ -236,14 +226,12 @@ struct PairDispersionD3KernelA { using NonDupScatterView = KKScatterView; typename AT::t_kkfloat_1d_3_lr_randomread x; - typename AT::t_kkacc_1d_3 f; 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; - typename AT::t_kkfloat_1d d_rcov; Kokkos::View d_c6ab; typename AT::t_int_1d d_mxci; @@ -262,7 +250,6 @@ struct PairDispersionD3KernelA { KK_FLOAT special_lj[4]; int nlocal; - int nall; int eflag; int vflag_either; int eflag_global; @@ -272,18 +259,15 @@ struct PairDispersionD3KernelA { int dampingCode; KK_FLOAT s6, s8, rs6, rs8, a1, a2, alpha; - KK_FLOAT cn_thr; PairDispersionD3KernelA( const typename AT::t_kkfloat_1d_3_lr_randomread &x_in, - const typename AT::t_kkacc_1d_3 &f_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 typename AT::t_kkfloat_1d &d_rcov_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, @@ -298,25 +282,25 @@ struct PairDispersionD3KernelA { const DupScatterView &dup_dc6_in, const NonDupScatterView &ndup_dc6_in, const KK_FLOAT *special_lj_in, - int nlocal_in, int nall_in, int eflag_in, int vflag_either_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, KK_FLOAT cn_thr_in) - : x(x_in), f(f_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_rcov(d_rcov_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), nall(nall_in), eflag(eflag_in), vflag_either(vflag_either_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), cn_thr(cn_thr_in) + alpha(alpha_in) { special_lj[0] = special_lj_in[0]; special_lj[1] = special_lj_in[1]; @@ -331,22 +315,6 @@ struct PairDispersionD3KernelA { return j >> SBBITS & 3; } - // Init per-thread reduction variable - // NOLINTNEXTLINE - KOKKOS_INLINE_FUNCTION - void init(value_type &ev) const { - ev.evdwl = 0.0; - for (int i = 0; i < 6; i++) ev.v[i] = 0.0; - } - - // Combine per-thread reduction variables - // NOLINTNEXTLINE - KOKKOS_INLINE_FUNCTION - void join(value_type &dst, const value_type &src) const { - dst.evdwl += src.evdwl; - for (int i = 0; i < 6; i++) dst.v[i] += src.v[i]; - } - /* ---------------------------------------------------------------------- Get derivative of C6 on device ------------------------------------------------------------------------- */ @@ -370,15 +338,12 @@ struct PairDispersionD3KernelA { d_den_i = 0.0; d_den_j = 0.0; - KK_FLOAT autoang6 = AUTOANG * AUTOANG * AUTOANG; - autoang6 = MathSpecialKokkos::square(autoang6); - 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; + c6_ref *= AUTOEV * AUTOANG6; if (c6_ref > 0) { cni_ref = d_c6ab(iat, jat, ci, cj, 1); @@ -391,7 +356,7 @@ struct PairDispersionD3KernelA { c6mem = c6_ref; } - expterm = exp(static_cast(K3) * r); + expterm = Kokkos::exp(static_cast(K3) * r); num += c6_ref * expterm; den += expterm; @@ -409,7 +374,10 @@ struct PairDispersionD3KernelA { } } - if (den > 1.0E-99) { + // 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); @@ -633,32 +601,8 @@ struct PairDispersionD3KernelA { fpair *= factor_lj; } break; - case 3: { // bj - - 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; - - case 4: { // bjm + case 3: // bj + case 4: { // bjm, same functional form as bj, different parameters const KK_FLOAT r0 = Kokkos::sqrt(C8 / C6); @@ -683,9 +627,7 @@ struct PairDispersionD3KernelA { fpair *= factor_lj; } break; - default: { - Kokkos::abort("PairDispersionD3Kokkos: invalid dampingCode"); - } break; + // no default case: dampingCode is validated in init_style() } if (EVFLAG) evdwl = -(s6 * e6 + s8 * e8) * factor_lj; @@ -748,7 +690,6 @@ struct PairDispersionD3KernelB { using NonDupScatterView = KKScatterView; typename AT::t_kkfloat_1d_3_lr_randomread x; - typename AT::t_kkacc_1d_3 f; typename AT::t_int_1d_randomread type; typename AT::t_kkfloat_2d d_cutsq; typename AT::t_kkfloat_1d d_dc6; @@ -756,7 +697,6 @@ struct PairDispersionD3KernelB { typename AT::t_int_1d d_numneigh; typename AT::t_neighbors_2d d_neighbors; typename AT::t_int_1d d_ilist; - int nall; DupScatterView dup_f; NonDupScatterView ndup_f; @@ -777,7 +717,6 @@ struct PairDispersionD3KernelB { PairDispersionD3KernelB( const typename AT::t_kkfloat_1d_3_lr_randomread &x_in, - const typename AT::t_kkacc_1d_3 &f_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, @@ -785,7 +724,6 @@ struct PairDispersionD3KernelB { 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, - int nall_in, const DupScatterView &dup_f_in, const NonDupScatterView &ndup_f_in, const DupScatterView &dup_eatom_in, @@ -796,9 +734,9 @@ struct PairDispersionD3KernelB { 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), f(f_in), type(type_in), d_cutsq(d_cutsq_in), d_dc6(d_dc6_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), nall(nall_in), dup_f(dup_f_in), ndup_f(ndup_f_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), @@ -819,22 +757,6 @@ struct PairDispersionD3KernelB { return j >> SBBITS & 3; } - // Init per-thread reduction variable - // NOLINTNEXTLINE - KOKKOS_INLINE_FUNCTION - void init(value_type &ev) const { - ev.evdwl = 0.0; - for (int i = 0; i < 6; i++) ev.v[i] = 0.0; - } - - // Combine thread reduction variables - // NOLINTNEXTLINE - KOKKOS_INLINE_FUNCTION - void join(value_type &dst, const value_type &src) const { - dst.evdwl += src.evdwl; - for (int i = 0; i < 6; i++) dst.v[i] += src.v[i]; - } - // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION @@ -971,7 +893,7 @@ struct PairDispersionD3KernelB { if (rsq < cn_thr) { const KK_FLOAT rcovij = (d_rcov(itype) + d_rcov(jtype)) * AUTOANG; - const KK_FLOAT expterm = exp(-K1 * (rcovij / r - 1.0)); + 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; @@ -1022,11 +944,11 @@ class PairDispersionD3Kokkos : public PairDispersionD3, public KokkosBase { PairDispersionD3Kokkos(class LAMMPS *); ~PairDispersionD3Kokkos() override; - void calc_coordination_number(); + void calc_coordination_number() override; void compute(int, int) override; double init_one(int, int) override; void init_style() override; - void allocate(); + void allocate() override; void coeff(int, char **) override; int pack_forward_comm_kokkos(int, DAT::tdual_int_1d, DAT::tdual_double_1d&, @@ -1038,12 +960,6 @@ class PairDispersionD3Kokkos : public PairDispersionD3, public KokkosBase { void unpack_forward_comm(int, int, double *) override; int pack_reverse_comm(int, int, double *) override; void unpack_reverse_comm(int, int *, double *) override; - // NOLINTNEXTLINE - KOKKOS_INLINE_FUNCTION - static int sbmask(const int &j) { - return j >> SBBITS & 3; - } - protected: typename AT::t_kkfloat_1d_3_lr_randomread x; // atom positions @@ -1095,7 +1011,7 @@ class PairDispersionD3Kokkos : public PairDispersionD3, public KokkosBase { typename AT::t_kkfloat_2d d_r0ab; // device R0 table Kokkos::View d_c6ab; // device C6 table - typename AT::tdual_kkfloat_2d k_cutsq; // cutoff^2 table (dual view) + 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 @@ -1111,6 +1027,9 @@ class PairDispersionD3Kokkos : public PairDispersionD3, public KokkosBase { 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); diff --git a/src/OPENMP/pair_dispersion_d3_omp.cpp b/src/OPENMP/pair_dispersion_d3_omp.cpp index a6294ea3c98..fe9611f145c 100644 --- a/src/OPENMP/pair_dispersion_d3_omp.cpp +++ b/src/OPENMP/pair_dispersion_d3_omp.cpp @@ -22,33 +22,19 @@ #include "atom.h" #include "comm.h" -#include "domain.h" -#include "error.h" #include "force.h" #include "math_special.h" #include "memory.h" #include "neigh_list.h" -#include "neighbor.h" #include "suffix.h" -#include "update.h" // needed for Update::ntimestep -#include #include -#include // for setenv #include -#include -#include #include "omp_compat.h" using namespace LAMMPS_NS; - -// global ad hoc parameters - copied from pair_dispersion_d3.cpp -static constexpr double K1 = 16.0; -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 +using namespace LAMMPS_NS::DispersionD3; /* ---------------------------------------------------------------------- */ @@ -62,23 +48,21 @@ PairDispersionD3OMP::PairDispersionD3OMP(LAMMPS *lmp) : void PairDispersionD3OMP::calc_coordination_number() { - int nlocal = atom->nlocal; const int nthreads = comm->nthreads; - int nall = nlocal + atom->nghost; + const int nall = atom->nlocal + atom->nghost; + const int newton_pair = force->newton_pair; - 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, nmax, "pair:cn"); - memory->grow(dc6, nmax, "pair:dc6"); + memory->grow(cn, nthreads * nmax, "pair:cn"); + memory->grow(dc6, nthreads * nmax, "pair:dc6"); } - // zero out coordination number - memset(cn, 0, sizeof(double) * (newton_pair ? nall : nlocal)); - memset(dc6, 0, sizeof(double) * (newton_pair ? nall : nlocal)); - - int inum = list->inum; + const int inum = list->inum; // Begin parallel region, the central atoms indexed by ii are assigned to different threads. #if defined(_OPENMP) @@ -119,12 +103,16 @@ void PairDispersionD3OMP::eval_coordination(int iifrom, int iito, ThrData * cons 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; - // Thread-local cn array to avoid race conditions - auto thr_cn = std::vector(atom->nmax); // Initialize to zero + // 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++) { @@ -159,12 +147,11 @@ void PairDispersionD3OMP::eval_coordination(int iifrom, int iito, ThrData * cons } } - // Contribute thread-local cn to global cn - for (int i = 0; i < atom->nmax; i++) { - // No need for atomic if newton_pair is false - #pragma omp atomic - cn[i] += thr_cn[i]; - } + // sum the per thread copies into cn[0] ... cn[nall-1] + + sync_threads(); + data_reduce_thr(cn, nall, nthreads, 1, tid); + sync_threads(); } /* ---------------------------------------------------------------------- */ @@ -172,6 +159,11 @@ void PairDispersionD3OMP::eval_coordination(int iifrom, int iito, ThrData * cons 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(); @@ -179,9 +171,6 @@ void PairDispersionD3OMP::compute(int eflag, int vflag) const int nthreads = comm->nthreads; const int inum = list->inum; - // Zero out dc6 values before OpenMP section - memset(dc6, 0, sizeof(double) * nall); - // Parallel direct force computation and some other quantities calculation. #if defined(_OPENMP) #pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag) \ @@ -216,6 +205,12 @@ firstprivate(inum,nthreads,nall) 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) { @@ -255,77 +250,8 @@ firstprivate(inum,nthreads,nall) else eval_second_phase<0,0,0>(ifrom, ito, thr); } thr->timer(Timer::PAIR); - reduce_thr(this, eflag, vflag_either, thr); + reduce_thr(this, eflag, vflag, thr); } //end of omp parallel region - - if (vflag_fdotr) virial_fdotr_compute(); -} - -/* ---------------------------------------------------------------------- - Modified from serial code to avoid race conditions -------------------------------------------------------------------------- */ - -void PairDispersionD3OMP::get_dC6(int iat, int jat, double cni, double cnj, double c6_res[3]) -{ - double c6_ref, cni_ref, cnj_ref; - double c6mem, r_save, r; - double expterm, term; - double num, den, d_num_i, d_num_j, d_den_i, d_den_j; - - c6mem = -1.0e20, r_save = 1.0e20; - num = 0; - den = 0; - d_num_i = 0; - d_num_j = 0; - d_den_i = 0; - d_den_j = 0; - - for (int ci = 0; ci <= mxci[iat]; ci++) { - for (int cj = 0; cj <= mxci[jat]; cj++) { - - c6_ref = c6ab[iat][jat][ci][cj][0]; - double autoang6 = AUTOANG * AUTOANG * AUTOANG; - autoang6 = MathSpecial::square(autoang6); - c6_ref *= AUTOEV * autoang6; - - if (c6_ref > 0) { - cni_ref = c6ab[iat][jat][ci][cj][1]; - cnj_ref = 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 = exp(static_cast(K3) * static_cast(r)); - - num += c6_ref * expterm; - den += expterm; - - expterm = expterm * 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; - } - } - } - - if (den > 1.0E-99) { - c6_res[0] = num / den; - c6_res[1] = ((d_num_i * den) - (d_den_i * num)) / (den * den); - c6_res[2] = ((d_num_j * den) - (d_den_j * num)) / (den * den); - } else { - c6_res[0] = c6mem; - c6_res[1] = 0; - c6_res[2] = 0; - } } /* ---------------------------------------------------------------------- */ @@ -341,10 +267,14 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const 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; - // Thread-local dc6 array to avoid race conditions - auto thr_dc6 = std::vector(atom->nmax); // Initialize to zero. + // 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) { @@ -377,7 +307,7 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const double r8inv = r2inv * r2inv * r2inv * r2inv; double r10inv = r2inv * r2inv * r2inv * r2inv * r2inv; - // Modified from original code to avoid race conditions + // get_dC6 writes {C6, dC6/dCN_i, dC6/dCN_j} double c6_res[3] = {}; get_dC6(itype, jtype, cn[i], cn[j], c6_res); @@ -448,7 +378,8 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const 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); double r4 = rsq * rsq; @@ -471,34 +402,6 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const fpair = -(tmp6 + tmp8); fpair *= factor_lj; } break; - - case 4: { // bjm - 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; - - default: { - error->all(FLERR, Error::NOLASTLINE, "Damping code {} unknown", dampingCode); - } break; } if (EFLAG) evdwl = -(s6 * e6 + s8 * e8) * factor_lj; @@ -530,12 +433,11 @@ void PairDispersionD3OMP::eval_first_phase(int iifrom, int iito, ThrData * const } } - // Contribute thread-local dc6 to global dc6 - for (int i = 0; i < atom->nmax; i++) { - // Possibly no need for atomic if newton_pair is false - #pragma omp atomic - dc6[i] += thr_dc6[i]; - } + // sum the per thread copies into dc6[0] ... dc6[nall-1] + + sync_threads(); + data_reduce_thr(dc6, nall, nthreads, 1, tid); + sync_threads(); } /* ---------------------------------------------------------------------- */ @@ -618,5 +520,7 @@ 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 index c874b822142..ef0cd80917b 100644 --- a/src/OPENMP/pair_dispersion_d3_omp.h +++ b/src/OPENMP/pair_dispersion_d3_omp.h @@ -39,14 +39,10 @@ class PairDispersionD3OMP : public PairDispersionD3, public ThrOMP { template void eval_coordination(int iifrom, int iito, ThrData *const thr); - void calc_coordination_number(); - - void get_dC6(int iat, int jat, double cni, double cnj, double c6_res[3]); + void calc_coordination_number() override; template void eval_second_phase(int ifrom, int ito, ThrData *const thr); - - int thr_evflag; // Flag for ev_tally thread usage }; } // namespace LAMMPS_NS From edfbfb3500b998f1f8535de9ab23d07783191bd0 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 25 Aug 2026 17:18:58 +0000 Subject: [PATCH 17/17] KOKKOS: use the standard scatter view duplication rule in dispersion/d3/kk The kernels forced ScatterNonDuplicated whenever DeviceType was LMPDeviceType. In a GPU build that matches NeedDup_v anyway, but in a CPU only build LMPDeviceType and LMPHostType are the same type, so the "half" neighbor style with more than one thread - which maps to HALFTHREAD - fell back to atomics instead of the duplicated scatter views every other pair style uses there. Use NeedDup_v directly, as PairComputeFunctor does, and drop the matching special case from the need_dup computation in compute() so the two stay consistent (they must agree, or ScatterViewHelper hands the functor the wrong half of the dup_/ndup_ pair). Also drop two variables left unused by the previous commit, so both changed OPENMP and EXTRA-PAIR files compile clean under -Wall -Wextra. Checked against the plain style with "compute stress/atom" active: energy, pressure tensor and summed per-atom stress are bit identical for neigh half and neigh full, newton on and off, at 1 and 4 threads. --- src/KOKKOS/pair_dispersion_d3_kokkos.cpp | 7 ++++--- src/KOKKOS/pair_dispersion_d3_kokkos.h | 12 +++--------- src/OPENMP/pair_dispersion_d3_omp.cpp | 3 +-- 3 files changed, 8 insertions(+), 14 deletions(-) diff --git a/src/KOKKOS/pair_dispersion_d3_kokkos.cpp b/src/KOKKOS/pair_dispersion_d3_kokkos.cpp index c99cb4e295c..33e8db392e9 100644 --- a/src/KOKKOS/pair_dispersion_d3_kokkos.cpp +++ b/src/KOKKOS/pair_dispersion_d3_kokkos.cpp @@ -308,9 +308,10 @@ void PairDispersionD3Kokkos::compute(int eflag_in, int vflag_in) d_ilist = k_list->d_ilist; inum = list->inum; - if constexpr (std::is_same_v) { - need_dup = false; - } else if (neighflag == FULL) { + // 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>; diff --git a/src/KOKKOS/pair_dispersion_d3_kokkos.h b/src/KOKKOS/pair_dispersion_d3_kokkos.h index ba9fd4c6c55..5d7027939bd 100644 --- a/src/KOKKOS/pair_dispersion_d3_kokkos.h +++ b/src/KOKKOS/pair_dispersion_d3_kokkos.h @@ -122,9 +122,7 @@ struct PairDispersionD3UnpackReverseCommFunctor { template struct PairDispersionD3CoordinationNumberKernel { typedef ArrayTypes AT; - using DUP = std::conditional_t, - Kokkos::Experimental::ScatterNonDuplicated, - NeedDup_v>; + using DUP = NeedDup_v; using ScatterAccess = std::conditional_t< std::is_same_v, Kokkos::Experimental::ScatterNonAtomic, @@ -211,9 +209,7 @@ template struct PairDispersionD3KernelA { typedef ArrayTypes AT; using value_type = EV_FLOAT; - using DUP = std::conditional_t, - Kokkos::Experimental::ScatterNonDuplicated, - NeedDup_v>; + using DUP = NeedDup_v; using ScatterAccess = std::conditional_t< std::is_same_v, Kokkos::Experimental::ScatterNonAtomic, @@ -675,9 +671,7 @@ template struct PairDispersionD3KernelB { typedef ArrayTypes AT; using value_type = EV_FLOAT; - using DUP = std::conditional_t, - Kokkos::Experimental::ScatterNonDuplicated, - NeedDup_v>; + using DUP = NeedDup_v; using ScatterAccess = std::conditional_t< std::is_same_v, Kokkos::Experimental::ScatterNonAtomic, diff --git a/src/OPENMP/pair_dispersion_d3_omp.cpp b/src/OPENMP/pair_dispersion_d3_omp.cpp index fe9611f145c..27bb4cda2e2 100644 --- a/src/OPENMP/pair_dispersion_d3_omp.cpp +++ b/src/OPENMP/pair_dispersion_d3_omp.cpp @@ -49,7 +49,6 @@ PairDispersionD3OMP::PairDispersionD3OMP(LAMMPS *lmp) : void PairDispersionD3OMP::calc_coordination_number() { const int nthreads = comm->nthreads; - const int nall = atom->nlocal + atom->nghost; const int newton_pair = force->newton_pair; // cn and dc6 hold one copy per thread, each nall long. The threads only @@ -454,7 +453,7 @@ void PairDispersionD3OMP::eval_second_phase(int iifrom, int iito, ThrData * cons const int * _noalias const numneigh = list->numneigh; const int * const * const firstneigh = list->firstneigh; - double dc6tmp,xtmp,ytmp,ztmp,delx,dely,delz,rsq,factor_lj,dcn,rcovij,expterm,fpair,fxtmp,fytmp,fztmp,r; + 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++) {