Merge pull request #1507 from athomps/yarray

Back-porting of Zombie SNAP improvements

doc/src/compute_sna_atom.txt

@@ -24,12 +24,7 @@ twojmax = band limit for bispectrum components (non-negative integer) :l
 R_1, R_2,... = list of cutoff radii, one for each type (distance units) :l
 w_1, w_2,... = list of neighbor weights, one for each type  :l
 zero or more keyword/value pairs may be appended :l
-keyword = {diagonal} or {rmin0} or {switchflag} or {bzeroflag} or {quadraticflag} :l
+keyword = {rmin0} or {switchflag} or {bzeroflag} or {quadraticflag} :l
-  {diagonal} value = {0} or {1} or {2} or {3}
-     {0} = all j1, j2, j <= twojmax, j2 <= j1
-     {1} = subset satisfying j1 == j2
-     {2} = subset satisfying j1 == j2 == j3
-     {3} = subset satisfying j2 <= j1 <= j
   {rmin0} value = parameter in distance to angle conversion (distance units)
   {switchflag} value = {0} or {1}
      {0} = do not use switching function
@@ -44,7 +39,7 @@ keyword = {diagonal} or {rmin0} or {switchflag} or {bzeroflag} or {quadraticflag
 
 [Examples:]
 
-compute b all sna/atom 1.4 0.99363 6 2.0 2.4 0.75 1.0 diagonal 3 rmin0 0.0
+compute b all sna/atom 1.4 0.99363 6 2.0 2.4 0.75 1.0 rmin0 0.0
 compute db all sna/atom 1.4 0.95 6 2.0 1.0
 compute vb all sna/atom 1.4 0.95 6 2.0 1.0 :pre
 
@@ -151,7 +146,7 @@ The argument {rfac0} and the optional keyword {rmin0} define the
 linear mapping from radial distance to polar angle {theta0} on the
 3-sphere.
 
-The argument {twojmax} and the keyword {diagonal} define which
+The argument {twojmax} defines which
 bispectrum components are generated. See section below on output for a
 detailed explanation of the number of bispectrum components and the
 ordered in which they are listed.
@@ -192,23 +187,18 @@ command that includes all pairs in the neighbor list.
 Compute {sna/atom} calculates a per-atom array, each column
 corresponding to a particular bispectrum component.  The total number
 of columns and the identity of the bispectrum component contained in
-each column depend on the values of {twojmax} and {diagonal}, as
+each column depend of the value of {twojmax}, as
 described by the following piece of python code:
 
 for j1 in range(0,twojmax+1):
-    if(diagonal==2):
+    for j2 in range(0,j1+1):
-        print j1/2.,j1/2.,j1/2.
+        for j in range(j1-j2,min(twojmax,j1+j2)+1,2):
-    elif(diagonal==1):
+            if (j>=j1): print j1/2.,j2/2.,j/2. :pre
-        for j in range(0,min(twojmax,2*j1)+1,2):
+
-            print j1/2.,j1/2.,j/2.
+NOTE: the {diagonal} keyword allowing other possible choices
-    elif(diagonal==0):
+for the number of bispectrum components was removed in 2019, 
-        for j2 in range(0,j1+1):
+since all potentials use the value of 3, corresponding to the
-            for j in range(j1-j2,min(twojmax,j1+j2)+1,2):
+above set of bispectrum components.
-                print j1/2.,j2/2.,j/2.
-    elif(diagonal==3):
-        for j2 in range(0,j1+1):
-            for j in range(j1-j2,min(twojmax,j1+j2)+1,2):
-                if (j>=j1): print j1/2.,j2/2.,j/2. :pre
 
 Compute {snad/atom} evaluates a per-atom array. The columns are
 arranged into {ntypes} blocks, listed in order of atom type {I}.  Each
@@ -259,7 +249,7 @@ package"_Build_package.html doc page for more info.
 
 [Default:]
 
-The optional keyword defaults are {diagonal} = 0, {rmin0} = 0,
+The optional keyword defaults are {rmin0} = 0,
 {switchflag} = 1, {bzeroflag} = 1, {quadraticflag} = 0,
 
 :line

doc/src/pair_snap.txt

@@ -38,7 +38,7 @@ where {B_k^i} is the {k}-th bispectrum component of atom {i},
 and {beta_k^alpha_i} is the corresponding linear coefficient
 that depends on {alpha_i}, the SNAP element of atom {i}. The
 number of bispectrum components used and their definitions
-depend on the values of {twojmax} and {diagonalstyle}
+depend on the value of {twojmax}
 defined in the SNAP parameter file described below.
 The bispectrum calculation is described in more detail
 in "compute sna/atom"_compute_sna_atom.html.
@@ -125,14 +125,13 @@ This line is followed by {ncoeff} coefficients, one per line.
 The SNAP parameter file can contain blank and comment lines (start
 with #) anywhere. Each non-blank non-comment line must contain one
 keyword/value pair. The required keywords are {rcutfac} and
-{twojmax}. Optional keywords are {rfac0}, {rmin0}, {diagonalstyle},
+{twojmax}. Optional keywords are {rfac0}, {rmin0},
 {switchflag}, and {bzeroflag}.
 
 The default values for these keywords are
 
 {rfac0} = 0.99363
 {rmin0} = 0.0
-{diagonalstyle} = 3
 {switchflag} = 0
 {bzeroflag} = 1
 {quadraticflag} = 1 :ul
@@ -144,6 +143,9 @@ If {quadraticflag} is set to 1, then the SNAP energy expression includes the qua
 The SNAP element file should contain {K}({K}+1)/2 additional coefficients
 for each element, the upper-triangular elements of alpha.
 
+NOTE: The previously used {diagonalstyle} keyword was removed in 2019,
+since all known SNAP potentials use the default value of 3.
+
 :line
 
 [Mixing, shift, table, tail correction, restart, rRESPA info]:

potentials/Ta06A.snapparam

@@ -10,6 +10,5 @@ twojmax 6
 
 rfac0 0.99363
 rmin0 0
-diagonalstyle 3
 bzeroflag 0
 quadraticflag 0

potentials/W_2940_2017_2.snapparam

@@ -8,6 +8,5 @@ twojmax 8
 
 rfac0 0.99363
 rmin0 0
-diagonalstyle 3
 bzeroflag 0
 quadraticflag 0

src/KOKKOS/pair_snap_kokkos_impl.h

@@ -85,9 +85,6 @@ void PairSNAPKokkos<DeviceType>::init_style()
   if (force->newton_pair == 0)
     error->all(FLERR,"Pair style SNAP requires newton pair on");
 
-  if (diagonalstyle != 3)
-    error->all(FLERR,"Must use diagonal style = 3 with pair snap/kk");
-
   // irequest = neigh request made by parent class
 
   neighflag = lmp->kokkos->neighflag;
@@ -343,23 +340,12 @@ void PairSNAPKokkos<DeviceType>::coeff(int narg, char **arg)
   Kokkos::deep_copy(d_coeffelem,h_coeffelem);
   Kokkos::deep_copy(d_map,h_map);
 
-  // deallocate non-kokkos sna
-
-  if (sna) {
-    for (int tid = 0; tid<nthreads; tid++)
-      delete sna[tid];
-    delete [] sna;
-    sna = NULL;
-  }
-
   // allocate memory for per OpenMP thread data which
   // is wrapped into the sna class
 
   snaKK = SNAKokkos<DeviceType>(rfac0,twojmax,
-                  diagonalstyle,use_shared_arrays,
                   rmin0,switchflag,bzeroflag);
-    //if (!use_shared_arrays)
+  snaKK.grow_rij(0);
-  snaKK.grow_rij(nmax);
   snaKK.init();
 }
 
@@ -667,8 +653,6 @@ double PairSNAPKokkos<DeviceType>::memory_usage()
   int n = atom->ntypes+1;
   bytes += n*n*sizeof(int);
   bytes += n*n*sizeof(double);
-  bytes += 3*nmax*sizeof(double);
-  bytes += nmax*sizeof(int);
   bytes += (2*ncoeffall)*sizeof(double);
   bytes += (ncoeff*3)*sizeof(double);
   bytes += snaKK.memory_usage();

src/KOKKOS/sna_kokkos.h

@@ -48,7 +48,7 @@ inline
   SNAKokkos(const SNAKokkos<DeviceType>& sna, const typename Kokkos::TeamPolicy<DeviceType>::member_type& team);
 
 inline
-  SNAKokkos(double, int, int, int, double, int, int);
+  SNAKokkos(double, int, double, int, int);
 
   KOKKOS_INLINE_FUNCTION
   ~SNAKokkos();
@@ -178,12 +178,6 @@ inline
                        double, double, double, // compute_duidrj
                        double, double, double, double, double);
 
-  // if number of atoms are small use per atom arrays
-  // for twojmax arrays, rij, inside, bvec
-  // this will increase the memory footprint considerably,
-  // but allows parallel filling and reuse of these arrays
-  int use_shared_arrays;
-
   // Sets the style for the switching function
   // 0 = none
   // 1 = cosine

src/KOKKOS/sna_kokkos_impl.h

@@ -27,19 +27,17 @@ static const double MY_PI  = 3.14159265358979323846; // pi
 template<class DeviceType>
 inline
 SNAKokkos<DeviceType>::SNAKokkos(double rfac0_in,
-         int twojmax_in, int diagonalstyle_in, int use_shared_arrays_in,
+         int twojmax_in,
          double rmin0_in, int switch_flag_in, int bzero_flag_in)
 {
   wself = 1.0;
 
-  use_shared_arrays = use_shared_arrays_in;
   rfac0 = rfac0_in;
   rmin0 = rmin0_in;
   switch_flag = switch_flag_in;
   bzero_flag = bzero_flag_in;
 
   twojmax = twojmax_in;
-  diagonalstyle = diagonalstyle_in;
 
   ncoeff = compute_ncoeff();
 
@@ -70,14 +68,12 @@ KOKKOS_INLINE_FUNCTION
 SNAKokkos<DeviceType>::SNAKokkos(const SNAKokkos<DeviceType>& sna, const typename Kokkos::TeamPolicy<DeviceType>::member_type& team) {
   wself = sna.wself;
 
-  use_shared_arrays = sna.use_shared_arrays;
   rfac0 = sna.rfac0;
   rmin0 = sna.rmin0;
   switch_flag = sna.switch_flag;
   bzero_flag = sna.bzero_flag;
 
   twojmax = sna.twojmax;
-  diagonalstyle = sna.diagonalstyle;
 
   ncoeff = sna.ncoeff;
   nmax = sna.nmax;
@@ -104,48 +100,45 @@ template<class DeviceType>
 inline
 void SNAKokkos<DeviceType>::build_indexlist()
 {
-  if(diagonalstyle == 3) {
+  int idxj_count = 0;
-    int idxj_count = 0;
+  int idxj_full_count = 0;
-    int idxj_full_count = 0;
 
-    for(int j1 = 0; j1 <= twojmax; j1++)
+  for(int j1 = 0; j1 <= twojmax; j1++)
-      for(int j2 = 0; j2 <= j1; j2++)
+    for(int j2 = 0; j2 <= j1; j2++)
-        for(int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2) {
+      for(int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2) {
-          if (j >= j1) idxj_count++;
+	if (j >= j1) idxj_count++;
-          idxj_full_count++;
+	idxj_full_count++;
-        }
+      }
 
-    // indexList can be changed here
+  // indexList can be changed here
 
-    idxj = Kokkos::View<SNAKK_LOOPINDICES*, DeviceType>("SNAKokkos::idxj",idxj_count);
+  idxj = Kokkos::View<SNAKK_LOOPINDICES*, DeviceType>("SNAKokkos::idxj",idxj_count);
-    idxj_full = Kokkos::View<SNAKK_LOOPINDICES*, DeviceType>("SNAKokkos::idxj_full",idxj_full_count);
+  idxj_full = Kokkos::View<SNAKK_LOOPINDICES*, DeviceType>("SNAKokkos::idxj_full",idxj_full_count);
-    auto h_idxj = Kokkos::create_mirror_view(idxj);
+  auto h_idxj = Kokkos::create_mirror_view(idxj);
-    auto h_idxj_full = Kokkos::create_mirror_view(idxj_full);
+  auto h_idxj_full = Kokkos::create_mirror_view(idxj_full);
 
-    idxj_max = idxj_count;
+  idxj_max = idxj_count;
-    idxj_full_max = idxj_full_count;
+  idxj_full_max = idxj_full_count;
 
-    idxj_count = 0;
+  idxj_count = 0;
-    idxj_full_count = 0;
+  idxj_full_count = 0;
 
-    for(int j1 = 0; j1 <= twojmax; j1++)
+  for(int j1 = 0; j1 <= twojmax; j1++)
-      for(int j2 = 0; j2 <= j1; j2++)
+    for(int j2 = 0; j2 <= j1; j2++)
-        for(int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2) {
+      for(int j = abs(j1 - j2); j <= MIN(twojmax, j1 + j2); j += 2) {
-          if (j >= j1) {
+	if (j >= j1) {
-            h_idxj[idxj_count].j1 = j1;
+	  h_idxj[idxj_count].j1 = j1;
-            h_idxj[idxj_count].j2 = j2;
+	  h_idxj[idxj_count].j2 = j2;
-            h_idxj[idxj_count].j = j;
+	  h_idxj[idxj_count].j = j;
-            idxj_count++;
+	  idxj_count++;
-          }
+	}
-          h_idxj_full[idxj_full_count].j1 = j1;
+	h_idxj_full[idxj_full_count].j1 = j1;
-          h_idxj_full[idxj_full_count].j2 = j2;
+	h_idxj_full[idxj_full_count].j2 = j2;
-          h_idxj_full[idxj_full_count].j = j;
+	h_idxj_full[idxj_full_count].j = j;
-          idxj_full_count++;
+	idxj_full_count++;
-        }
+      }
-    Kokkos::deep_copy(idxj,h_idxj);
+  Kokkos::deep_copy(idxj,h_idxj);
-    Kokkos::deep_copy(idxj_full,h_idxj_full);
+  Kokkos::deep_copy(idxj_full,h_idxj_full);
-
-  }
 
 }
 /* ---------------------------------------------------------------------- */
@@ -1223,26 +1216,10 @@ int SNAKokkos<DeviceType>::compute_ncoeff()
   ncount = 0;
 
   for (int j1 = 0; j1 <= twojmax; j1++)
-    if(diagonalstyle == 0) {
+    for (int j2 = 0; j2 <= j1; j2++)
-      for (int j2 = 0; j2 <= j1; j2++)
-        for (int j = abs(j1 - j2);
-            j <= MIN(twojmax, j1 + j2); j += 2)
-          ncount++;
-    } else if(diagonalstyle == 1) {
-      int j2 = j1;
-
       for (int j = abs(j1 - j2);
-          j <= MIN(twojmax, j1 + j2); j += 2)
+	   j <= MIN(twojmax, j1 + j2); j += 2)
-        ncount++;
+	if (j >= j1) ncount++;
-    } else if(diagonalstyle == 2) {
-      ncount++;
-    } else if(diagonalstyle == 3) {
-      for (int j2 = 0; j2 <= j1; j2++)
-        for (int j = abs(j1 - j2);
-            j <= MIN(twojmax, j1 + j2); j += 2)
-          if (j >= j1) ncount++;
-    }
-
   return ncount;
 }
 

src/SNAP/compute_sna_atom.cpp

@@ -25,7 +25,6 @@
 #include "comm.h"
 #include "memory.h"
 #include "error.h"
-#include "openmp_snap.h"
 
 using namespace LAMMPS_NS;
 
@@ -45,7 +44,6 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
 
   // default values
 
-  diagonalstyle = 0;
   rmin0 = 0.0;
   switchflag = 1;
   bzeroflag = 1;
@@ -85,14 +83,7 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
   int iarg = nargmin;
 
   while (iarg < narg) {
-    if (strcmp(arg[iarg],"diagonal") == 0) {
+    if (strcmp(arg[iarg],"rmin0") == 0) {
-      if (iarg+2 > narg)
-        error->all(FLERR,"Illegal compute sna/atom command");
-      diagonalstyle = atoi(arg[iarg+1]);
-      if (diagonalstyle < 0 || diagonalstyle > 3)
-        error->all(FLERR,"Illegal compute sna/atom command");
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"rmin0") == 0) {
       if (iarg+2 > narg)
         error->all(FLERR,"Illegal compute sna/atom command");
       rmin0 = atof(arg[iarg+1]);
@@ -115,28 +106,16 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
     } else error->all(FLERR,"Illegal compute sna/atom command");
   }
 
-  nthreads = comm->nthreads;
+  snaptr = new SNA(lmp,rfac0,twojmax,
-  snaptr = new SNA*[nthreads];
+                   rmin0,switchflag,bzeroflag);
-#if defined(_OPENMP)
-#pragma omp parallel default(none) shared(lmp,rfac0,twojmax,rmin0,switchflag,bzeroflag)
-#endif
-  {
-    int tid = omp_get_thread_num();
 
-    // always unset use_shared_arrays since it does not work with computes
+  ncoeff = snaptr->ncoeff;
-    snaptr[tid] = new SNA(lmp,rfac0,twojmax,diagonalstyle,
-                          0 /*use_shared_arrays*/, rmin0,switchflag,bzeroflag);
-  }
-
-  ncoeff = snaptr[0]->ncoeff;
   size_peratom_cols = ncoeff;
   if (quadraticflag) size_peratom_cols += (ncoeff*(ncoeff+1))/2;
   peratom_flag = 1;
 
   nmax = 0;
-  njmax = 0;
   sna = NULL;
-
 }
 
 /* ---------------------------------------------------------------------- */
@@ -147,9 +126,7 @@ ComputeSNAAtom::~ComputeSNAAtom()
   memory->destroy(radelem);
   memory->destroy(wjelem);
   memory->destroy(cutsq);
-  for (int tid = 0; tid<nthreads; tid++)
+  delete snaptr;
-    delete snaptr[tid];
-  delete [] snaptr;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -176,13 +153,7 @@ void ComputeSNAAtom::init()
     if (strcmp(modify->compute[i]->style,"sna/atom") == 0) count++;
   if (count > 1 && comm->me == 0)
     error->warning(FLERR,"More than one compute sna/atom");
-#if defined(_OPENMP)
+  snaptr->init();
-#pragma omp parallel default(none)
-#endif
-  {
-    int tid = omp_get_thread_num();
-    snaptr[tid]->init();
-  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -223,11 +194,7 @@ void ComputeSNAAtom::compute_peratom()
   double** const x = atom->x;
   const int* const mask = atom->mask;
 
-#if defined(_OPENMP)
-#pragma omp parallel for default(none)
-#endif
   for (int ii = 0; ii < inum; ii++) {
-    const int tid = omp_get_thread_num();
     const int i = ilist[ii];
     if (mask[i] & groupbit) {
 
@@ -241,7 +208,7 @@ void ComputeSNAAtom::compute_peratom()
 
       // insure rij, inside, and typej  are of size jnum
 
-      snaptr[tid]->grow_rij(jnum);
+      snaptr->grow_rij(jnum);
 
       // rij[][3] = displacements between atom I and those neighbors
       // inside = indices of neighbors of I within cutoff
@@ -258,26 +225,25 @@ void ComputeSNAAtom::compute_peratom()
         const double rsq = delx*delx + dely*dely + delz*delz;
         int jtype = type[j];
         if (rsq < cutsq[itype][jtype] && rsq>1e-20) {
-          snaptr[tid]->rij[ninside][0] = delx;
+          snaptr->rij[ninside][0] = delx;
-          snaptr[tid]->rij[ninside][1] = dely;
+          snaptr->rij[ninside][1] = dely;
-          snaptr[tid]->rij[ninside][2] = delz;
+          snaptr->rij[ninside][2] = delz;
-          snaptr[tid]->inside[ninside] = j;
+          snaptr->inside[ninside] = j;
-          snaptr[tid]->wj[ninside] = wjelem[jtype];
+          snaptr->wj[ninside] = wjelem[jtype];
-          snaptr[tid]->rcutij[ninside] = (radi+radelem[jtype])*rcutfac;
+          snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac;
           ninside++;
         }
       }
 
-      snaptr[tid]->compute_ui(ninside);
+      snaptr->compute_ui(ninside);
-      snaptr[tid]->compute_zi();
+      snaptr->compute_zi();
-      snaptr[tid]->compute_bi();
+      snaptr->compute_bi();
-      snaptr[tid]->copy_bi2bvec();
       for (int icoeff = 0; icoeff < ncoeff; icoeff++)
-        sna[i][icoeff] = snaptr[tid]->bvec[icoeff];
+        sna[i][icoeff] = snaptr->blist[icoeff];
       if (quadraticflag) {
         int ncount = ncoeff;
         for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
-          double bi = snaptr[tid]->bvec[icoeff];
+          double bi = snaptr->blist[icoeff];
 
           // diagonal element of quadratic matrix
 
@@ -286,7 +252,7 @@ void ComputeSNAAtom::compute_peratom()
           // upper-triangular elements of quadratic matrix
 
           for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++)
-            sna[i][ncount++] = bi*snaptr[tid]->bvec[jcoeff];
+            sna[i][ncount++] = bi*snaptr->blist[jcoeff];
         }
       }
     } else {
@@ -302,10 +268,9 @@ void ComputeSNAAtom::compute_peratom()
 
 double ComputeSNAAtom::memory_usage()
 {
-  double bytes = nmax*size_peratom_cols * sizeof(double);
+  double bytes = nmax*size_peratom_cols * sizeof(double); // sna
-  bytes += 3*njmax*sizeof(double);
+  bytes += snaptr->memory_usage();                        // SNA object
-  bytes += njmax*sizeof(int);
+
-  bytes += snaptr[0]->memory_usage()*comm->nthreads;
   return bytes;
 }
 

src/SNAP/compute_sna_atom.h

@@ -34,7 +34,7 @@ class ComputeSNAAtom : public Compute {
   double memory_usage();
 
  private:
-  int nmax, njmax, diagonalstyle;
+  int nmax;
   int ncoeff;
   double **cutsq;
   class NeighList *list;
@@ -42,10 +42,9 @@ class ComputeSNAAtom : public Compute {
   double rcutfac;
   double *radelem;
   double *wjelem;
-  class SNA** snaptr;
+  class SNA* snaptr;
   double cutmax;
   int quadraticflag;
-  int nthreads;
 };
 
 }

src/SNAP/compute_snad_atom.cpp

@@ -25,7 +25,6 @@
 #include "comm.h"
 #include "memory.h"
 #include "error.h"
-#include "openmp_snap.h"
 
 using namespace LAMMPS_NS;
 
@@ -45,7 +44,6 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
 
   // default values
 
-  diagonalstyle = 0;
   rmin0 = 0.0;
   switchflag = 1;
   bzeroflag = 1;
@@ -83,14 +81,7 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
   int iarg = nargmin;
 
   while (iarg < narg) {
-    if (strcmp(arg[iarg],"diagonal") == 0) {
+    if (strcmp(arg[iarg],"rmin0") == 0) {
-      if (iarg+2 > narg)
-        error->all(FLERR,"Illegal compute snad/atom command");
-      diagonalstyle = atof(arg[iarg+1]);
-      if (diagonalstyle < 0 || diagonalstyle > 3)
-        error->all(FLERR,"Illegal compute snad/atom command");
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"rmin0") == 0) {
       if (iarg+2 > narg)
         error->all(FLERR,"Illegal compute snad/atom command");
       rmin0 = atof(arg[iarg+1]);
@@ -113,20 +104,10 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
     } else error->all(FLERR,"Illegal compute snad/atom command");
   }
 
-  nthreads = comm->nthreads;
+  snaptr = new SNA(lmp,rfac0,twojmax,
-  snaptr = new SNA*[nthreads];
+                   rmin0,switchflag,bzeroflag);
-#if defined(_OPENMP)
-#pragma omp parallel default(none) shared(lmp,rfac0,twojmax,rmin0,switchflag,bzeroflag)
-#endif
-  {
-    int tid = omp_get_thread_num();
 
-    // always unset use_shared_arrays since it does not work with computes
+  ncoeff = snaptr->ncoeff;
-    snaptr[tid] = new SNA(lmp,rfac0,twojmax,diagonalstyle,
-                          0 /*use_shared_arrays*/, rmin0,switchflag,bzeroflag);
-  }
-
-  ncoeff = snaptr[0]->ncoeff;
   nperdim = ncoeff;
   if (quadraticflag) nperdim += (ncoeff*(ncoeff+1))/2;
   yoffset = nperdim;
@@ -136,9 +117,7 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
   peratom_flag = 1;
 
   nmax = 0;
-  njmax = 0;
   snad = NULL;
-
 }
 
 /* ---------------------------------------------------------------------- */
@@ -149,9 +128,7 @@ ComputeSNADAtom::~ComputeSNADAtom()
   memory->destroy(radelem);
   memory->destroy(wjelem);
   memory->destroy(cutsq);
-  for (int tid = 0; tid<nthreads; tid++)
+  delete snaptr;
-    delete snaptr[tid];
-  delete [] snaptr;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -178,13 +155,7 @@ void ComputeSNADAtom::init()
     if (strcmp(modify->compute[i]->style,"snad/atom") == 0) count++;
   if (count > 1 && comm->me == 0)
     error->warning(FLERR,"More than one compute snad/atom");
-#if defined(_OPENMP)
+  snaptr->init();
-#pragma omp parallel default(none)
-#endif
-  {
-    int tid = omp_get_thread_num();
-    snaptr[tid]->init();
-  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -235,11 +206,7 @@ void ComputeSNADAtom::compute_peratom()
   double** const x = atom->x;
   const int* const mask = atom->mask;
 
-#if defined(_OPENMP)
-#pragma omp parallel for default(none)
-#endif
   for (int ii = 0; ii < inum; ii++) {
-    const int tid = omp_get_thread_num();
     const int i = ilist[ii];
     if (mask[i] & groupbit) {
 
@@ -258,7 +225,7 @@ void ComputeSNADAtom::compute_peratom()
 
       // insure rij, inside, and typej  are of size jnum
 
-      snaptr[tid]->grow_rij(jnum);
+      snaptr->grow_rij(jnum);
 
       // rij[][3] = displacements between atom I and those neighbors
       // inside = indices of neighbors of I within cutoff
@@ -276,30 +243,28 @@ void ComputeSNADAtom::compute_peratom()
         const double rsq = delx*delx + dely*dely + delz*delz;
         int jtype = type[j];
         if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
-          snaptr[tid]->rij[ninside][0] = delx;
+          snaptr->rij[ninside][0] = delx;
-          snaptr[tid]->rij[ninside][1] = dely;
+          snaptr->rij[ninside][1] = dely;
-          snaptr[tid]->rij[ninside][2] = delz;
+          snaptr->rij[ninside][2] = delz;
-          snaptr[tid]->inside[ninside] = j;
+          snaptr->inside[ninside] = j;
-          snaptr[tid]->wj[ninside] = wjelem[jtype];
+          snaptr->wj[ninside] = wjelem[jtype];
-          snaptr[tid]->rcutij[ninside] = (radi+radelem[jtype])*rcutfac;
+          snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac;
           ninside++;
         }
       }
 
-      snaptr[tid]->compute_ui(ninside);
+      snaptr->compute_ui(ninside);
-      snaptr[tid]->compute_zi();
+      snaptr->compute_zi();
       if (quadraticflag) {
-        snaptr[tid]->compute_bi();
+        snaptr->compute_bi();
-        snaptr[tid]->copy_bi2bvec();
       }
 
       for (int jj = 0; jj < ninside; jj++) {
-        const int j = snaptr[tid]->inside[jj];
+        const int j = snaptr->inside[jj];
-        snaptr[tid]->compute_duidrj(snaptr[tid]->rij[jj],
+        snaptr->compute_duidrj(snaptr->rij[jj],
-                                    snaptr[tid]->wj[jj],
+                                    snaptr->wj[jj],
-                                    snaptr[tid]->rcutij[jj]);
+                                    snaptr->rcutij[jj]);
-        snaptr[tid]->compute_dbidrj();
+        snaptr->compute_dbidrj();
-        snaptr[tid]->copy_dbi2dbvec();
 
         // Accumulate -dBi/dRi, -dBi/dRj
 
@@ -307,12 +272,12 @@ void ComputeSNADAtom::compute_peratom()
         double *snadj = snad[j]+typeoffset;
 
         for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
-          snadi[icoeff] += snaptr[tid]->dbvec[icoeff][0];
+          snadi[icoeff] += snaptr->dblist[icoeff][0];
-          snadi[icoeff+yoffset] += snaptr[tid]->dbvec[icoeff][1];
+          snadi[icoeff+yoffset] += snaptr->dblist[icoeff][1];
-          snadi[icoeff+zoffset] += snaptr[tid]->dbvec[icoeff][2];
+          snadi[icoeff+zoffset] += snaptr->dblist[icoeff][2];
-          snadj[icoeff] -= snaptr[tid]->dbvec[icoeff][0];
+          snadj[icoeff] -= snaptr->dblist[icoeff][0];
-          snadj[icoeff+yoffset] -= snaptr[tid]->dbvec[icoeff][1];
+          snadj[icoeff+yoffset] -= snaptr->dblist[icoeff][1];
-          snadj[icoeff+zoffset] -= snaptr[tid]->dbvec[icoeff][2];
+          snadj[icoeff+zoffset] -= snaptr->dblist[icoeff][2];
         }
 
         if (quadraticflag) {
@@ -321,10 +286,10 @@ void ComputeSNADAtom::compute_peratom()
           snadj += quadraticoffset;
           int ncount = 0;
           for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
-            double bi = snaptr[tid]->bvec[icoeff];
+            double bi = snaptr->blist[icoeff];
-            double bix = snaptr[tid]->dbvec[icoeff][0];
+            double bix = snaptr->dblist[icoeff][0];
-            double biy = snaptr[tid]->dbvec[icoeff][1];
+            double biy = snaptr->dblist[icoeff][1];
-            double biz = snaptr[tid]->dbvec[icoeff][2];
+            double biz = snaptr->dblist[icoeff][2];
 
             // diagonal elements of quadratic matrix
 
@@ -343,12 +308,12 @@ void ComputeSNADAtom::compute_peratom()
             // upper-triangular elements of quadratic matrix
 
             for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
-              double dbxtmp = bi*snaptr[tid]->dbvec[jcoeff][0]
+              double dbxtmp = bi*snaptr->dblist[jcoeff][0]
-                + bix*snaptr[tid]->bvec[jcoeff];
+                + bix*snaptr->blist[jcoeff];
-              double dbytmp = bi*snaptr[tid]->dbvec[jcoeff][1]
+              double dbytmp = bi*snaptr->dblist[jcoeff][1]
-                + biy*snaptr[tid]->bvec[jcoeff];
+                + biy*snaptr->blist[jcoeff];
-              double dbztmp = bi*snaptr[tid]->dbvec[jcoeff][2]
+              double dbztmp = bi*snaptr->dblist[jcoeff][2]
-                + biz*snaptr[tid]->bvec[jcoeff];
+                + biz*snaptr->blist[jcoeff];
 
               snadi[ncount] +=         dbxtmp;
               snadi[ncount+yoffset] += dbytmp;
@@ -404,10 +369,9 @@ void ComputeSNADAtom::unpack_reverse_comm(int n, int *list, double *buf)
 
 double ComputeSNADAtom::memory_usage()
 {
-  double bytes = nmax*size_peratom_cols * sizeof(double);
+
-  bytes += 3*njmax*sizeof(double);
+  double bytes = nmax*size_peratom_cols * sizeof(double); // snad
-  bytes += njmax*sizeof(int);
+  bytes += snaptr->memory_usage();                        // SNA object
-  bytes += 3*nperdim*atom->ntypes;
+
-  bytes += snaptr[0]->memory_usage()*comm->nthreads;
   return bytes;
 }

src/SNAP/compute_snad_atom.h

@@ -36,7 +36,7 @@ class ComputeSNADAtom : public Compute {
   double memory_usage();
 
  private:
-  int nmax, njmax, diagonalstyle;
+  int nmax;
   int ncoeff, nperdim, yoffset, zoffset;
   double **cutsq;
   class NeighList *list;
@@ -44,10 +44,9 @@ class ComputeSNADAtom : public Compute {
   double rcutfac;
   double *radelem;
   double *wjelem;
-  class SNA** snaptr;
+  class SNA* snaptr;
   double cutmax;
   int quadraticflag;
-  int nthreads;
 };
 
 }

src/SNAP/compute_snav_atom.cpp

@@ -25,7 +25,6 @@
 #include "comm.h"
 #include "memory.h"
 #include "error.h"
-#include "openmp_snap.h"
 
 using namespace LAMMPS_NS;
 
@@ -45,7 +44,6 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
 
   // default values
 
-  diagonalstyle = 0;
   rmin0 = 0.0;
   switchflag = 1;
   bzeroflag = 1;
@@ -79,14 +77,7 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
   int iarg = nargmin;
 
   while (iarg < narg) {
-    if (strcmp(arg[iarg],"diagonal") == 0) {
+    if (strcmp(arg[iarg],"rmin0") == 0) {
-      if (iarg+2 > narg)
-        error->all(FLERR,"Illegal compute snav/atom command");
-      diagonalstyle = atof(arg[iarg+1]);
-      if (diagonalstyle < 0 || diagonalstyle > 3)
-        error->all(FLERR,"Illegal compute snav/atom command");
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"rmin0") == 0) {
       if (iarg+2 > narg)
         error->all(FLERR,"Illegal compute snav/atom command");
       rmin0 = atof(arg[iarg+1]);
@@ -109,20 +100,10 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
     } else error->all(FLERR,"Illegal compute snav/atom command");
   }
 
-  nthreads = comm->nthreads;
+  snaptr = new SNA(lmp,rfac0,twojmax,
-  snaptr = new SNA*[nthreads];
+                   rmin0,switchflag,bzeroflag);
-#if defined(_OPENMP)
-#pragma omp parallel default(none) shared(lmp,rfac0,twojmax,rmin0,switchflag,bzeroflag)
-#endif
-  {
-    int tid = omp_get_thread_num();
 
-    // always unset use_shared_arrays since it does not work with computes
+  ncoeff = snaptr->ncoeff;
-    snaptr[tid] = new SNA(lmp,rfac0,twojmax,diagonalstyle,
-                          0 /*use_shared_arrays*/, rmin0,switchflag,bzeroflag);
-  }
-
-  ncoeff = snaptr[0]->ncoeff;
   nperdim = ncoeff;
   if (quadraticflag) nperdim += (ncoeff*(ncoeff+1))/2;
   size_peratom_cols = 6*nperdim*atom->ntypes;
@@ -130,7 +111,6 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
   peratom_flag = 1;
 
   nmax = 0;
-  njmax = 0;
   snav = NULL;
 
 }
@@ -144,9 +124,7 @@ ComputeSNAVAtom::~ComputeSNAVAtom()
   memory->destroy(wjelem);
   memory->destroy(cutsq);
 
-  for (int tid = 0; tid<nthreads; tid++)
+  delete snaptr;
-    delete snaptr[tid];
-  delete [] snaptr;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -174,13 +152,7 @@ void ComputeSNAVAtom::init()
     if (strcmp(modify->compute[i]->style,"snav/atom") == 0) count++;
   if (count > 1 && comm->me == 0)
     error->warning(FLERR,"More than one compute snav/atom");
-#if defined(_OPENMP)
+  snaptr->init();
-#pragma omp parallel default(none)
-#endif
-  {
-    int tid = omp_get_thread_num();
-    snaptr[tid]->init();
-  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -230,11 +202,7 @@ void ComputeSNAVAtom::compute_peratom()
   double** const x = atom->x;
   const int* const mask = atom->mask;
 
-#if defined(_OPENMP)
-#pragma omp parallel for default(none)
-#endif
   for (int ii = 0; ii < inum; ii++) {
-    const int tid = omp_get_thread_num();
     const int i = ilist[ii];
     if (mask[i] & groupbit) {
 
@@ -251,7 +219,7 @@ void ComputeSNAVAtom::compute_peratom()
 
       // insure rij, inside, and typej  are of size jnum
 
-          snaptr[tid]->grow_rij(jnum);
+      snaptr->grow_rij(jnum);
 
       // rij[][3] = displacements between atom I and those neighbors
       // inside = indices of neighbors of I within cutoff
@@ -269,31 +237,29 @@ void ComputeSNAVAtom::compute_peratom()
         const double rsq = delx*delx + dely*dely + delz*delz;
         int jtype = type[j];
         if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
-          snaptr[tid]->rij[ninside][0] = delx;
+          snaptr->rij[ninside][0] = delx;
-          snaptr[tid]->rij[ninside][1] = dely;
+          snaptr->rij[ninside][1] = dely;
-          snaptr[tid]->rij[ninside][2] = delz;
+          snaptr->rij[ninside][2] = delz;
-          snaptr[tid]->inside[ninside] = j;
+          snaptr->inside[ninside] = j;
-          snaptr[tid]->wj[ninside] = wjelem[jtype];
+          snaptr->wj[ninside] = wjelem[jtype];
-          snaptr[tid]->rcutij[ninside] = (radi+radelem[jtype])*rcutfac;
+          snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac;
           ninside++;
         }
       }
 
-      snaptr[tid]->compute_ui(ninside);
+      snaptr->compute_ui(ninside);
-      snaptr[tid]->compute_zi();
+      snaptr->compute_zi();
       if (quadraticflag) {
-        snaptr[tid]->compute_bi();
+        snaptr->compute_bi();
-        snaptr[tid]->copy_bi2bvec();
       }
 
       for (int jj = 0; jj < ninside; jj++) {
-        const int j = snaptr[tid]->inside[jj];
+        const int j = snaptr->inside[jj];
 
-        snaptr[tid]->compute_duidrj(snaptr[tid]->rij[jj],
+        snaptr->compute_duidrj(snaptr->rij[jj],
-                                    snaptr[tid]->wj[jj],
+                                    snaptr->wj[jj],
-                                    snaptr[tid]->rcutij[jj]);
+                                    snaptr->rcutij[jj]);
-        snaptr[tid]->compute_dbidrj();
+        snaptr->compute_dbidrj();
-        snaptr[tid]->copy_dbi2dbvec();
 
         // Accumulate -dBi/dRi*Ri, -dBi/dRj*Rj
 
@@ -301,18 +267,18 @@ void ComputeSNAVAtom::compute_peratom()
         double *snavj = snav[j]+typeoffset;
 
         for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
-          snavi[icoeff]           += snaptr[tid]->dbvec[icoeff][0]*xtmp;
+          snavi[icoeff]           += snaptr->dblist[icoeff][0]*xtmp;
-          snavi[icoeff+nperdim]   += snaptr[tid]->dbvec[icoeff][1]*ytmp;
+          snavi[icoeff+nperdim]   += snaptr->dblist[icoeff][1]*ytmp;
-          snavi[icoeff+2*nperdim] += snaptr[tid]->dbvec[icoeff][2]*ztmp;
+          snavi[icoeff+2*nperdim] += snaptr->dblist[icoeff][2]*ztmp;
-          snavi[icoeff+3*nperdim] += snaptr[tid]->dbvec[icoeff][1]*ztmp;
+          snavi[icoeff+3*nperdim] += snaptr->dblist[icoeff][1]*ztmp;
-          snavi[icoeff+4*nperdim] += snaptr[tid]->dbvec[icoeff][0]*ztmp;
+          snavi[icoeff+4*nperdim] += snaptr->dblist[icoeff][0]*ztmp;
-          snavi[icoeff+5*nperdim] += snaptr[tid]->dbvec[icoeff][0]*ytmp;
+          snavi[icoeff+5*nperdim] += snaptr->dblist[icoeff][0]*ytmp;
-          snavj[icoeff]           -= snaptr[tid]->dbvec[icoeff][0]*x[j][0];
+          snavj[icoeff]           -= snaptr->dblist[icoeff][0]*x[j][0];
-          snavj[icoeff+nperdim]   -= snaptr[tid]->dbvec[icoeff][1]*x[j][1];
+          snavj[icoeff+nperdim]   -= snaptr->dblist[icoeff][1]*x[j][1];
-          snavj[icoeff+2*nperdim] -= snaptr[tid]->dbvec[icoeff][2]*x[j][2];
+          snavj[icoeff+2*nperdim] -= snaptr->dblist[icoeff][2]*x[j][2];
-          snavj[icoeff+3*nperdim] -= snaptr[tid]->dbvec[icoeff][1]*x[j][2];
+          snavj[icoeff+3*nperdim] -= snaptr->dblist[icoeff][1]*x[j][2];
-          snavj[icoeff+4*nperdim] -= snaptr[tid]->dbvec[icoeff][0]*x[j][2];
+          snavj[icoeff+4*nperdim] -= snaptr->dblist[icoeff][0]*x[j][2];
-          snavj[icoeff+5*nperdim] -= snaptr[tid]->dbvec[icoeff][0]*x[j][1];
+          snavj[icoeff+5*nperdim] -= snaptr->dblist[icoeff][0]*x[j][1];
         }
 
         if (quadraticflag) {
@@ -321,10 +287,10 @@ void ComputeSNAVAtom::compute_peratom()
           snavj += quadraticoffset;
           int ncount = 0;
           for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
-            double bi = snaptr[tid]->bvec[icoeff];
+            double bi = snaptr->blist[icoeff];
-            double bix = snaptr[tid]->dbvec[icoeff][0];
+            double bix = snaptr->dblist[icoeff][0];
-            double biy = snaptr[tid]->dbvec[icoeff][1];
+            double biy = snaptr->dblist[icoeff][1];
-            double biz = snaptr[tid]->dbvec[icoeff][2];
+            double biz = snaptr->dblist[icoeff][2];
 
             // diagonal element of quadratic matrix
 
@@ -348,12 +314,12 @@ void ComputeSNAVAtom::compute_peratom()
             // upper-triangular elements of quadratic matrix
 
             for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
-              double dbxtmp = bi*snaptr[tid]->dbvec[jcoeff][0]
+              double dbxtmp = bi*snaptr->dblist[jcoeff][0]
-                + bix*snaptr[tid]->bvec[jcoeff];
+                + bix*snaptr->blist[jcoeff];
-              double dbytmp = bi*snaptr[tid]->dbvec[jcoeff][1]
+              double dbytmp = bi*snaptr->dblist[jcoeff][1]
-                + biy*snaptr[tid]->bvec[jcoeff];
+                + biy*snaptr->blist[jcoeff];
-              double dbztmp = bi*snaptr[tid]->dbvec[jcoeff][2]
+              double dbztmp = bi*snaptr->dblist[jcoeff][2]
-                + biz*snaptr[tid]->bvec[jcoeff];
+                + biz*snaptr->blist[jcoeff];
               snavi[ncount] +=           dbxtmp*xtmp;
               snavi[ncount+nperdim] +=   dbytmp*ytmp;
               snavi[ncount+2*nperdim] += dbztmp*ztmp;
@@ -414,11 +380,8 @@ void ComputeSNAVAtom::unpack_reverse_comm(int n, int *list, double *buf)
 
 double ComputeSNAVAtom::memory_usage()
 {
-  double bytes = nmax*size_peratom_cols * sizeof(double);
+  double bytes = nmax*size_peratom_cols * sizeof(double); // snav
-  bytes += 3*njmax*sizeof(double);
+  bytes += snaptr->memory_usage();                        // SNA object
-  bytes += njmax*sizeof(int);
+
-  bytes += 6*nperdim*atom->ntypes;
-  if (quadraticflag) bytes += 6*nperdim*atom->ntypes;
-  bytes += snaptr[0]->memory_usage()*comm->nthreads;
   return bytes;
 }

src/SNAP/compute_snav_atom.h

@@ -36,7 +36,7 @@ class ComputeSNAVAtom : public Compute {
   double memory_usage();
 
  private:
-  int nmax, njmax, diagonalstyle;
+  int nmax;
   int ncoeff, nperdim;
   double **cutsq;
   class NeighList *list;
@@ -44,9 +44,8 @@ class ComputeSNAVAtom : public Compute {
   double rcutfac;
   double *radelem;
   double *wjelem;
-  class SNA** snaptr;
+  class SNA* snaptr;
   int quadraticflag;
-  int nthreads;
 };
 
 }

src/SNAP/openmp_snap.h

@@ -1,16 +0,0 @@
-
-#ifndef LMP_OPENMP_SNAP_H
-#define LMP_OPENMP_SNAP_H
-
-#if defined(_OPENMP)
-#include <omp.h>
-#else
-enum omp_sched_t {omp_sched_static, omp_sched_dynamic, omp_sched_guided, omp_sched_auto};
-inline int omp_get_thread_num() { return 0;}
-inline int omp_set_num_threads(int num_threads) {return 1;}
-/* inline int __sync_fetch_and_add(int* ptr, int value) {int tmp = *ptr; ptr[0]+=value; return tmp;} */
-inline void omp_set_schedule(omp_sched_t schedule,int modifier=1) {}
-inline int omp_in_parallel() {return 0;}
-#endif
-
-#endif

src/SNAP/pair_snap.h

@@ -29,8 +29,6 @@ public:
   PairSNAP(class LAMMPS *);
   ~PairSNAP();
   virtual void compute(int, int);
-  void compute_regular(int, int);
-  void compute_optimized(int, int);
   void settings(int, char **);
   virtual void coeff(int, char **);
   virtual void init_style();
@@ -42,56 +40,14 @@ public:
 
 protected:
   int ncoeffq, ncoeffall;
-  double **bvec, ***dbvec;
+  class SNA* snaptr;
-  class SNA** sna;
-  int nmax;
-  int nthreads;
   virtual void allocate();
   void read_files(char *, char *);
   inline int equal(double* x,double* y);
   inline double dist2(double* x,double* y);
-  double extra_cutoff();
-  void load_balance();
-  void set_sna_to_shared(int snaid,int i);
-  void build_per_atom_arrays();
 
-  int schedule_user;
+  void compute_beta();
-  double schedule_time_guided;
+  void compute_bispectrum();
-  double schedule_time_dynamic;
-
-  int ncalls_neigh;
-  int do_load_balance;
-  int ilistmask_max;
-  int* ilistmask;
-  int ghostinum;
-  int ghostilist_max;
-  int* ghostilist;
-  int ghostnumneigh_max;
-  int* ghostnumneigh;
-  int* ghostneighs;
-  int* ghostfirstneigh;
-  int ghostneighs_total;
-  int ghostneighs_max;
-
-  int use_optimized;
-  int use_shared_arrays;
-
-  int i_max;
-  int i_neighmax;
-  int i_numpairs;
-  int **i_pairs;
-  double ***i_rij;
-  int **i_inside;
-  double **i_wj;
-  double **i_rcutij;
-  int *i_ninside;
-  double ****i_uarraytot_r, ****i_uarraytot_i;
-  double ******i_zarray_r, ******i_zarray_i;
-
-#ifdef TIMING_INFO
-  //  timespec starttime, endtime;
-  double timers[4];
-#endif
 
   double rcutmax;               // max cutoff for all elements
   int nelements;                // # of unique elements
@@ -99,10 +55,13 @@ protected:
   double *radelem;              // element radii
   double *wjelem;               // elements weights
   double **coeffelem;           // element bispectrum coefficients
+  double** beta;                // betas for all atoms in list
+  double** bispectrum;          // bispectrum components for all atoms in list
   int *map;                     // mapping from atom types to elements
-  int twojmax, diagonalstyle, switchflag, bzeroflag;
+  int twojmax, switchflag, bzeroflag;
   double rfac0, rmin0, wj1, wj2;
   int rcutfacflag, twojmaxflag; // flags for required parameters
+  int beta_max;                 // length of beta
 };
 
 }
@@ -124,15 +83,6 @@ Self-explanatory.  Check the input script syntax and compare to the
 documentation for the command.  You can use -echo screen as a
 command-line option when running LAMMPS to see the offending line.
 
-E: Must set number of threads via package omp command
-
-Because you are using the USER-OMP package, set the number of threads
-via its settings, not by the pair_style snap nthreads setting.
-
-W: Communication cutoff is too small for SNAP micro load balancing, increased to %lf
-
-Self-explanatory.
-
 E: Incorrect args for pair coefficients
 
 Self-explanatory.  Check the input script or data file.

src/SNAP/sna.h

@@ -18,20 +18,22 @@
 #ifndef LMP_SNA_H
 #define LMP_SNA_H
 
-#include <complex>
 #include "pointers.h"
-#include <ctime>
 
 namespace LAMMPS_NS {
 
-struct SNA_LOOPINDICES {
+struct SNA_ZINDICES {
+  int j1, j2, j, ma1min, ma2max, mb1min, mb2max, na, nb, jju;
+};
+
+struct SNA_BINDICES {
   int j1, j2, j;
 };
 
 class SNA : protected Pointers {
 
 public:
-  SNA(LAMMPS*, double, int, int, int, double, int, int);
+  SNA(LAMMPS*, double, int, double, int, int);
 
   SNA(LAMMPS* lmp) : Pointers(lmp) {};
   ~SNA();
@@ -44,31 +46,21 @@ public:
   // functions for bispectrum coefficients
 
   void compute_ui(int);
-  void compute_ui_omp(int, int);
   void compute_zi();
-  void compute_zi_omp(int);
+  void compute_yi(const double*);
+  void compute_yterm(int, int, int, const double*);
   void compute_bi();
-  void copy_bi2bvec();
 
   // functions for derivatives
 
   void compute_duidrj(double*, double, double);
   void compute_dbidrj();
-  void compute_dbidrj_nonsymm();
+  void compute_deidrj(double*);
-  void copy_dbi2dbvec();
   double compute_sfac(double, double);
   double compute_dsfac(double, double);
 
-#ifdef TIMING_INFO
+  double* blist;
-  double* timers;
+  double** dblist;
-  timespec starttime, endtime;
-  int print;
-  int counter;
-#endif
-
-  //per sna class instance for OMP use
-
-  double* bvec, ** dbvec;
   double** rij;
   int* inside;
   double* wj;
@@ -77,29 +69,32 @@ public:
 
   void grow_rij(int);
 
-  int twojmax, diagonalstyle;
+  int twojmax;
-  double*** uarraytot_r, *** uarraytot_i;
-  double***** zarray_r, ***** zarray_i;
-  double*** uarraytot_r_b, *** uarraytot_i_b;
-  double***** zarray_r_b, ***** zarray_i_b;
-  double*** uarray_r, *** uarray_i;
 
 private:
   double rmin0, rfac0;
 
-  //use indexlist instead of loops, constructor generates these
-  SNA_LOOPINDICES* idxj;
-  int idxj_max;
   // data for bispectrum coefficients
 
-  double***** cgarray;
+  SNA_ZINDICES* idxz;
+  SNA_BINDICES* idxb;
+  int idxcg_max, idxu_max, idxz_max, idxb_max;
+
   double** rootpqarray;
-  double*** barray;
+  double* cglist;  
+  int*** idxcg_block; 
 
-  // derivatives of data
+  double* ulisttot_r, * ulisttot_i;
+  double* ulist_r, * ulist_i;
+  int* idxu_block;
 
-  double**** duarray_r, **** duarray_i;
+  double* zlist_r, * zlist_i;
-  double**** dbarray;
+  int*** idxz_block;
+
+  int*** idxb_block;
+
+  double** dulist_r, ** dulist_i;
+  double* ylist_r, * ylist_i;
 
   static const int nmaxfactorial = 167;
   static const double nfac_table[];
@@ -109,28 +104,16 @@ private:
   void destroy_twojmax_arrays();
   void init_clebsch_gordan();
   void init_rootpqarray();
-  void jtostr(char*, int);
-  void mtostr(char*, int, int);
-  void print_clebsch_gordan(FILE*);
   void zero_uarraytot();
   void addself_uarraytot(double);
   void add_uarraytot(double, double, double);
-  void add_uarraytot_omp(double, double, double);
   void compute_uarray(double, double, double,
                       double, double);
-  void compute_uarray_omp(double, double, double,
-                          double, double, int);
   double deltacg(int, int, int);
   int compute_ncoeff();
   void compute_duarray(double, double, double,
                        double, double, double, double, double);
 
-  // if number of atoms are small use per atom arrays
-  // for twojmax arrays, rij, inside, bvec
-  // this will increase the memory footprint considerably,
-  // but allows parallel filling and reuse of these arrays
-  int use_shared_arrays;
-
   // Sets the style for the switching function
   // 0 = none
   // 1 = cosine
@@ -140,7 +123,7 @@ private:
   double wself;
 
   int bzero_flag; // 1 if bzero subtracted from barray
-  double *bzero;  // array of B values for isolated atoms
+  double* bzero;  // array of B values for isolated atoms
 };
 
 }