Review fixes for upstream PR #5056 (dispersion/d3 OPENMP and KOKKOS variants)

doc/src/Commands_pair.rst

@@ -84,7 +84,7 @@ parenthesis: g = GPU, i = INTEL, k = KOKKOS, o = OPENMP, t = OPT.
    * :doc:`coul/tt <pair_coul_tt>`
    * :doc:`coul/wolf (ko) <pair_coul>`
    * :doc:`coul/wolf/cs (o) <pair_cs>`
-   * :doc:`dispersion/d3 <pair_dispersion_d3>`
+   * :doc:`dispersion/d3 (ko) <pair_dispersion_d3>`
    * :doc:`dpd (giko) <pair_dpd>`
    * :doc:`dpd/coul/slater/long (g) <pair_dpd_coul_slater_long>`
    * :doc:`dpd/ext (ko) <pair_dpd_ext>`

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
 

doc/src/pair_dispersion_d3.rst

@@ -1,8 +1,12 @@
 .. index:: pair_style dispersion/d3
+.. index:: pair_style dispersion/d3/kk
+.. index:: pair_style dispersion/d3/omp
 
 pair_style dispersion/d3 command
 ================================
 
+Accelerator Variants: *dispersion/d3/kk*, *dispersion/d3/omp*
+
 Syntax
 """"""
 
@@ -105,6 +109,11 @@ provide are the chemical symbols of the atoms.  The number of chemical
 symbols given must be equal to the number of atom types used and must
 match their ordering as atom types.
 
+----------
+
+.. include:: accel_styles.rst
+
+----------
 
 Mixing, shift, table, tail correction, restart, rRESPA info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

src/EXTRA-PAIR/pair_dispersion_d3.cpp

@@ -27,39 +27,26 @@
 #include "error.h"
 #include "force.h"
 #include "info.h"
+#include "math_special.h"
 #include "memory.h"
 #include "neigh_list.h"
 #include "neighbor.h"
 #include "update.h"
 
 #include <algorithm>
+#include <cctype>
 #include <cmath>
 #include <cstring>
-#include <cctype>
 #include <unordered_map>
 #include <utility>
 
 using namespace LAMMPS_NS;
-
-// global ad hoc parameters
-static constexpr double K1 = 16.0;
-static constexpr double K3 = -4.0;
-
-/*  reasonable choices for k3 are between 3 and 5 :
-    this gives smoth curves with maxima around the integer values
-    k3=3 give for CN=0 a slightly smaller value than computed
-    for the free atom. This also yields to larger CN for atoms
-    in larger molecules but with the same chemical environment
-    which is physically not right.
-    values >5 might lead to bumps in the potential.
-*/
+using namespace LAMMPS_NS::DispersionD3;
 
 static constexpr int NUM_ELEMENTS = 94;      // maximum element number
 static constexpr int N_PARS_COLS = 5;        // number of columns in C6 table
 static constexpr int N_PARS_ROWS = 32385;    // number of rows in C6 table
-
-static constexpr double autoang = 0.52917725;    // atomic units (Bohr) to Angstrom
-static constexpr double autoev = 27.21140795;    // atomic units (Hartree) to eV
+static constexpr int N_C6AB_GRID = 5;        // reference CN grid size of the C6 table
 
 #include "d3_parameters.h"
 
@@ -82,6 +69,7 @@ PairDispersionD3::PairDispersionD3(LAMMPS *lmp) :
 
   dampingCode = 0;
   s6 = s8 = s18 = rs6 = rs8 = rs18 = a1 = a2 = alpha = alpha6 = alpha8 = 0.0;
+  max_mxci = 0;
 }
 
 /* ----------------------------------------------------------------------
@@ -90,6 +78,8 @@ PairDispersionD3::PairDispersionD3(LAMMPS *lmp) :
 
 PairDispersionD3::~PairDispersionD3()
 {
+  if (copymode) return;
+
   if (allocated) {
     memory->destroy(setflag);
     memory->destroy(cutsq);
@@ -126,7 +116,7 @@ void PairDispersionD3::allocate()
   memory->create(rcov, n + 1, "pair:rcov");
 
   memory->create(r0ab, n + 1, n + 1, "pair:r0ab");
-  memory->create(c6ab, n + 1, n + 1, 5, 5, 3, "pair:c6ab");
+  memory->create(c6ab, n + 1, n + 1, N_C6AB_GRID, N_C6AB_GRID, 3, "pair:c6ab");
 }
 
 /* ----------------------------------------------------------------------
@@ -250,11 +240,15 @@ void PairDispersionD3::read_c6ab(int *atomic_numbers, int ntypes)
   for (int i = 0; i < N_PARS_ROWS; i++) {
     const double ref_c6 = c6ab_table[i][0];
 
-    int atom_number_1 = (int)std::round(c6ab_table[i][1]);
-    int atom_number_2 = (int)std::round(c6ab_table[i][2]);
+    int atom_number_1 = (int) std::round(c6ab_table[i][1]);
+    int atom_number_2 = (int) std::round(c6ab_table[i][2]);
 
     set_limit_in_pars_array(atom_number_1, atom_number_2, grid_i, grid_j);
 
+    if ((grid_i >= N_C6AB_GRID) || (grid_j >= N_C6AB_GRID))
+      error->all(FLERR, Error::NOLASTLINE,
+                 "Reference C6 grid index out of range in the compiled in D3 parameters");
+
     std::vector<int> idx_atoms_1 = is_int_in_array(atomic_numbers, ntypes, atom_number_1);
     if (idx_atoms_1.empty()) continue;
 
@@ -293,7 +287,7 @@ void PairDispersionD3::coeff(int narg, char **arg)
   if (!allocated) allocate();
 
   std::string element;
-  int *atomic_numbers = (int *) malloc(sizeof(int) * ntypes);
+  std::vector<int> atomic_numbers(ntypes);
   for (int i = 0; i < ntypes; i++) {
     element = arg[i + 2];
     atomic_numbers[i] = find_atomic_number(element);
@@ -317,12 +311,14 @@ void PairDispersionD3::coeff(int narg, char **arg)
   }
 
   // set r0ab
-  read_r0ab(atomic_numbers, ntypes);
+  read_r0ab(atomic_numbers.data(), ntypes);
 
-  // read c6ab
-  read_c6ab(atomic_numbers, ntypes);
+  // read c6ab and record the largest reference CN grid index that is in use.
+  // The accelerated variants size their coefficient tables from it.
 
-  free(atomic_numbers);
+  read_c6ab(atomic_numbers.data(), ntypes);
+  max_mxci = 0;
+  for (int i = 1; i <= ntypes; i++) max_mxci = std::max(max_mxci, mxci[i]);
 }
 
 /* ----------------------------------------------------------------------
@@ -377,7 +373,7 @@ void PairDispersionD3::calc_coordination_number()
       if (rsq > cn_thr) continue;
 
       double rr = sqrt(rsq);
-      double rcov_ij = (rcov[itype] + rcov[jtype]) * autoang;
+      double rcov_ij = (rcov[itype] + rcov[jtype]) * AUTOANG;
       double cn_ij = 1.0 / (1.0 + exp(-K1 * ((rcov_ij / rr) - 1.0)));
 
       // update coordination number
@@ -396,10 +392,8 @@ void PairDispersionD3::calc_coordination_number()
    Get derivative of C6
 ------------------------------------------------------------------------- */
 
-double *PairDispersionD3::get_dC6(int iat, int jat, double cni, double cnj)
+void PairDispersionD3::get_dC6(int iat, int jat, double cni, double cnj, double *c6_res)
 {
-
-  static double c6_res[3] = {};
   double c6_ref, cni_ref, cnj_ref;
   double c6mem, r_save, r;
   double expterm, term;
@@ -417,7 +411,7 @@ double *PairDispersionD3::get_dC6(int iat, int jat, double cni, double cnj)
     for (int cj = 0; cj <= mxci[jat]; cj++) {
 
       c6_ref = c6ab[iat][jat][ci][cj][0];
-      c6_ref *= autoev * pow(autoang, 6);
+      c6_ref *= AUTOEV * AUTOANG6;
 
       if (c6_ref > 0) {
         cni_ref = c6ab[iat][jat][ci][cj][1];
@@ -457,7 +451,6 @@ double *PairDispersionD3::get_dC6(int iat, int jat, double cni, double cnj)
     c6_res[1] = 0;
     c6_res[2] = 0;
   }
-  return c6_res;
 }
 
 /* ----------------------------------------------------------------------
@@ -510,16 +503,17 @@ void PairDispersionD3::compute(int eflag, int vflag)
 
       if (rsq < cutsq[type[i]][type[j]]) {
 
-        double r = sqrt(rsq);
         double r2inv = 1.0 / rsq;
         double r6inv = r2inv * r2inv * r2inv;
         double r8inv = r2inv * r2inv * r2inv * r2inv;
         double r10inv = r2inv * r2inv * r2inv * r2inv * r2inv;
 
-        double *c6_res = get_dC6(type[i], type[j], cn[i], cn[j]);
+        // get_dC6 writes {C6, dC6/dCN_i, dC6/dCN_j}
+        double c6_res[3] = {};
+        get_dC6(type[i], type[j], cn[i], cn[j], c6_res);
 
         double C6 = c6_res[0];
-        double C8 = 3.0 * C6 * r2r4[type[i]] * r2r4[type[j]] * autoang * autoang;
+        double C8 = 3.0 * C6 * r2r4[type[i]] * r2r4[type[j]] * AUTOANG * AUTOANG;
 
         double alpha6 = alpha;
         double alpha8 = alpha + 2;
@@ -530,13 +524,18 @@ void PairDispersionD3::compute(int eflag, int vflag)
 
         switch (dampingCode) {
 
-          case 1: {    // original
+          // Written to avoid using sqrt and pow()
+          case 1: /* Original damping */
+          {
+            double ip6 = rs6 * r0ab[type[i]][type[j]];
+            double ip8 = rs8 * r0ab[type[i]][type[j]];
 
-            double r0 = r / r0ab[type[i]][type[j]];
+            double half_alpha6 = 0.5 * alpha6;
+            double half_alpha8 = 0.5 * alpha8;
 
-            t6 = pow(rs6 / r0, alpha6);
+            t6 = MathSpecial::powauto(ip6, alpha6) * MathSpecial::powauto(rsq, -half_alpha6);
             damp6 = 1.0 / (1.0 + 6.0 * t6);
-            t8 = pow(rs8 / r0, alpha8);
+            t8 = MathSpecial::powauto(ip8, alpha8) * MathSpecial::powauto(rsq, -half_alpha8);
             damp8 = 1.0 / (1.0 + 6.0 * t8);
 
             e6 = C6 * damp6 * r6inv;
@@ -551,14 +550,15 @@ void PairDispersionD3::compute(int eflag, int vflag)
             fpair = fpair1 + fpair2;
             fpair *= factor_lj;
           } break;
-
+          // Written to avoid pow
           case 2: {    // zerom
 
+            double r = sqrt(rsq);
             double r0 = r0ab[type[i]][type[j]];
 
-            t6 = pow((r / (rs6 * r0)) + rs8 * r0, -alpha6);
+            t6 = MathSpecial::powauto((r / (rs6 * r0)) + rs8 * r0, -alpha6);
             damp6 = 1.0 / (1.0 + 6.0 * t6);
-            t8 = pow((r / r0) + rs8 * r0, -alpha8);
+            t8 = MathSpecial::powauto((r / r0) + rs8 * r0, -alpha8);
             damp8 = 1.0 / (1.0 + 6.0 * t8);
 
             e6 = C6 * damp6 * r6inv;
@@ -578,37 +578,21 @@ void PairDispersionD3::compute(int eflag, int vflag)
             fpair *= factor_lj;
           } break;
 
-          case 3: {    // bj
+          case 3:      // bj
+          case 4: {    // bjm, same functional form as bj, different parameters
 
             double r0 = sqrt(C8 / C6);
 
             double r4 = rsq * rsq;
             double r6 = rsq * rsq * rsq;
             double r8 = rsq * rsq * rsq * rsq;
 
-            t6 = r6 + pow((a1 * r0 + a2), 6);
-            t8 = r8 + pow((a1 * r0 + a2), 8);
-
-            e6 = C6 / t6;
-            e8 = C8 / t8;
-
-            tmp6 = 6.0 * s6 * C6 * r4 / (t6 * t6);
-            tmp8 = 8.0 * s8 * C8 * r6 / (t8 * t8);
-
-            fpair = -(tmp6 + tmp8);
-            fpair *= factor_lj;
-          } break;
-
-          case 4: {    // bjm
+            double d = a1 * r0 + a2;
+            double d2 = d * d;
+            double d4 = d2 * d2;
 
-            double r0 = sqrt(C8 / C6);
-
-            double r4 = rsq * rsq;
-            double r6 = rsq * rsq * rsq;
-            double r8 = rsq * rsq * rsq * rsq;
-
-            t6 = r6 + pow((a1 * r0 + a2), 6);
-            t8 = r8 + pow((a1 * r0 + a2), 8);
+            t6 = r6 + MathSpecial::cube(d2);
+            t8 = r8 + MathSpecial::square(d4);
 
             e6 = C6 / t6;
             e8 = C8 / t8;
@@ -685,7 +669,7 @@ void PairDispersionD3::compute(int eflag, int vflag)
 
         // here we calculate dcn = dCNi/dr = dCNj/dr
         if (rsq < cn_thr) {
-          double rcovij = (rcov[type[i]] + rcov[type[j]]) * autoang;
+          double rcovij = (rcov[type[i]] + rcov[type[j]]) * AUTOANG;
           double expterm = exp(-K1 * (rcovij / r - 1.0));
           dcn = -K1 * rcovij * expterm / (rsq * (expterm + 1.0) * (expterm + 1.0));
         } else {
@@ -1039,7 +1023,7 @@ void PairDispersionD3::set_funcpar(std::string &functional_name)
           break;
       }
 
-      rs8 = rs8 / autoang;
+      rs8 = rs8 / AUTOANG;
     } break;
 
     case 3: {    // bj
@@ -1356,7 +1340,7 @@ void PairDispersionD3::set_funcpar(std::string &functional_name)
           break;
       }
 
-      a2 = a2 * autoang;
+      a2 = a2 * AUTOANG;
     } break;
 
     case 4: {    // bjm
@@ -1417,8 +1401,7 @@ void PairDispersionD3::set_funcpar(std::string &functional_name)
           break;
       }
 
-      a2 = a2 * autoang;
-
+      a2 = a2 * AUTOANG;
     } break;
     default:
       // this should not happen with the error check in the init_style function
@@ -1447,6 +1430,14 @@ void PairDispersionD3::init_style()
 {
   if (atom->tag_enable == 0) error->all(FLERR, "Pair style D3 requires atom IDs");
 
+  // settings() only accepts known damping keywords, so this can not trigger.
+  // It is checked once here so that the threaded and device variants of the
+  // compute loop below do not have to report an error from inside a parallel
+  // region or a device kernel.
+
+  if ((dampingCode < 1) || (dampingCode > 4))
+    error->all(FLERR, Error::NOLASTLINE, "Damping code {} unknown", dampingCode);
+
   // need an half neighbor list
   neighbor->add_request(this);
 }
@@ -1531,8 +1522,8 @@ double PairDispersionD3::memory_usage()
 {
   double bytes = Pair::memory_usage();
   int n = atom->ntypes;
-  // c6ab[n+1][n+1][5][5][3] coefficient table
-  bytes += (double)(n+1)*(n+1)*5*5*3 * sizeof(double);
+  // c6ab[n+1][n+1][N_C6AB_GRID][N_C6AB_GRID][3] coefficient table
+  bytes += (double) (n + 1) * (n + 1) * N_C6AB_GRID * N_C6AB_GRID * 3 * sizeof(double);
   // per-atom coordination number and C6 derivative arrays
   bytes += (double) nmax * 2 * sizeof(double);    // cn[nmax] + dc6[nmax]
   return bytes;