From 75cffc78c5562fbc9bdce2e41c2a75968349173e Mon Sep 17 00:00:00 2001
From: "Steven J. Plimpton" <sjplimp@singsing.sandia.gov>
Date: Thu, 3 May 2018 10:23:50 -0600
Subject: [PATCH] updates to quadratic form of SNAP potential

---
 doc/src/compute_sna_atom.txt   |  21 ++++---
 src/SNAP/compute_sna_atom.cpp  |   6 +-
 src/SNAP/compute_snad_atom.cpp |  74 +++++++++++++---------
 src/SNAP/compute_snad_atom.h   |   2 +-
 src/SNAP/compute_snav_atom.cpp | 108 ++++++++++++++++++---------------
 src/SNAP/compute_snav_atom.h   |   3 +-
 src/SNAP/pair_snap.cpp         |  15 ++---
 src/SNAP/pair_snap.h           |   9 +--
 8 files changed, 135 insertions(+), 103 deletions(-)

diff --git a/doc/src/compute_sna_atom.txt b/doc/src/compute_sna_atom.txt
index 1c3787e696..268e23ac28 100644
--- a/doc/src/compute_sna_atom.txt
+++ b/doc/src/compute_sna_atom.txt
@@ -161,9 +161,9 @@ function.
 
 The keyword {bzeroflag} determines whether or not {B0}, the bispectrum
 components of an atom with no neighbors, are subtracted from
-the calculated bispectrum components. This optional keyword is only
-available for compute {sna/atom}, as {snad/atom} and {snav/atom}
-are unaffected by the removal of constant terms.
+the calculated bispectrum components. This optional keyword 
+normally only affects compute {sna/atom}. However, when
+{quadraticflag} is on, it also affects {snad/atom} and {snav/atom}.
 
 The keyword {quadraticflag} determines whether or not the
 quadratic analogs to the bispectrum quantities are generated.
@@ -230,13 +230,18 @@ are 30, 90, and 180, respectively. With {quadratic} value=1,
 the numbers of columns are 930, 2790, and 5580, respectively.
 
 If the {quadratic} keyword value is set to 1, then additional
-columns are appended to each per-atom array, corresponding to
+columns are generated, corresponding to
 the products of all distinct pairs of  bispectrum components. If the
 number of bispectrum components is {K}, then the number of distinct pairs
-is  {K}({K}+1)/2. These are output in subblocks of  {K}({K}+1)/2 columns, using the same
-ordering of sub-blocks as was used for the bispectrum
-components. Within each sub-block, the ordering is upper-triangular,
-(1,1),(1,2)...(1,{K}),(2,1)...({K}-1,{K}-1),({K}-1,{K}),({K},{K})
+is  {K}({K}+1)/2.
+For compute {sna/atom} these columns are appended to existing {K} columns.
+The ordering of quadratic terms is upper-triangular,
+(1,1),(1,2)...(1,{K}),(2,1)...({K}-1,{K}-1),({K}-1,{K}),({K},{K}).
+For computes {snad/atom} and {snav/atom} each set of {K}({K}+1)/2
+additional columns is inserted directly after each of sub-block
+of linear terms i.e. linear and quadratic terms are contiguous.
+So the nesting order from inside to outside is bispectrum component,
+linear then quadratic, vector/tensor component, type.
 
 These values can be accessed by any command that uses per-atom values
 from a compute as input.  See "Section
diff --git a/src/SNAP/compute_sna_atom.cpp b/src/SNAP/compute_sna_atom.cpp
index 4b114d9ce7..9f56cee3f1 100644
--- a/src/SNAP/compute_sna_atom.cpp
+++ b/src/SNAP/compute_sna_atom.cpp
@@ -276,9 +276,13 @@ void ComputeSNAAtom::compute_peratom()
         for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
           double bi = snaptr[tid]->bvec[icoeff];
 
+          // diagonal element of quadratic matrix
+
+          sna[i][ncount++] = 0.5*bi*bi;
+
           // upper-triangular elements of quadratic matrix
 
-          for (int jcoeff = icoeff; jcoeff < ncoeff; jcoeff++)
+          for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++)
             sna[i][ncount++] = bi*snaptr[tid]->bvec[jcoeff];
         }
       }
diff --git a/src/SNAP/compute_snad_atom.cpp b/src/SNAP/compute_snad_atom.cpp
index 54a6ce7612..586fb3d210 100644
--- a/src/SNAP/compute_snad_atom.cpp
+++ b/src/SNAP/compute_snad_atom.cpp
@@ -95,6 +95,11 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
         error->all(FLERR,"Illegal compute snad/atom command");
       rmin0 = atof(arg[iarg+1]);
       iarg += 2;
+    } else if (strcmp(arg[iarg],"bzeroflag") == 0) {
+      if (iarg+2 > narg)
+        error->all(FLERR,"Illegal compute snad/atom command");
+      bzeroflag = atoi(arg[iarg+1]);
+      iarg += 2;
     } else if (strcmp(arg[iarg],"switchflag") == 0) {
       if (iarg+2 > narg)
         error->all(FLERR,"Illegal compute snad/atom command");
@@ -121,16 +126,11 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
   }
 
   ncoeff = snaptr[0]->ncoeff;
-  twoncoeff = 2*ncoeff;
-  threencoeff = 3*ncoeff;
-  size_peratom_cols = threencoeff*atom->ntypes;
-  if (quadraticflag) {
-    ncoeffq = (ncoeff*(ncoeff+1))/2;
-    twoncoeffq = 2*ncoeffq;
-    threencoeffq = 3*ncoeffq;
-    size_peratom_cols +=
-      threencoeffq*atom->ntypes;
-  }
+  nperdim = ncoeff;  
+  if (quadraticflag) nperdim += (ncoeff*(ncoeff+1))/2;
+  yoffset = nperdim;
+  zoffset = 2*nperdim;
+  size_peratom_cols = 3*nperdim*atom->ntypes;
   comm_reverse = size_peratom_cols;
   peratom_flag = 1;
 
@@ -248,9 +248,10 @@ void ComputeSNADAtom::compute_peratom()
       const int* const jlist = firstneigh[i];
       const int jnum = numneigh[i];
 
-      const int typeoffset = threencoeff*(atom->type[i]-1);
-      const int quadraticoffset = threencoeff*atom->ntypes +
-        threencoeffq*(atom->type[i]-1);
+      // const int typeoffset = threencoeff*(atom->type[i]-1);
+      // const int quadraticoffset = threencoeff*atom->ntypes +
+      //   threencoeffq*(atom->type[i]-1);
+      const int typeoffset = 3*nperdim*(atom->type[i]-1);
 
       // insure rij, inside, and typej  are of size jnum
 
@@ -304,16 +305,17 @@ void ComputeSNADAtom::compute_peratom()
 
         for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
           snadi[icoeff] += snaptr[tid]->dbvec[icoeff][0];
-          snadi[icoeff+ncoeff] += snaptr[tid]->dbvec[icoeff][1];
-          snadi[icoeff+twoncoeff] += snaptr[tid]->dbvec[icoeff][2];
+          snadi[icoeff+yoffset] += snaptr[tid]->dbvec[icoeff][1];
+          snadi[icoeff+zoffset] += snaptr[tid]->dbvec[icoeff][2];
           snadj[icoeff] -= snaptr[tid]->dbvec[icoeff][0];
-          snadj[icoeff+ncoeff] -= snaptr[tid]->dbvec[icoeff][1];
-          snadj[icoeff+twoncoeff] -= snaptr[tid]->dbvec[icoeff][2];
+          snadj[icoeff+yoffset] -= snaptr[tid]->dbvec[icoeff][1];
+          snadj[icoeff+zoffset] -= snaptr[tid]->dbvec[icoeff][2];
         }
 
         if (quadraticflag) {
-          double *snadi = snad[i]+quadraticoffset;
-          double *snadj = snad[j]+quadraticoffset;
+          const int quadraticoffset = ncoeff;
+          snadi += quadraticoffset;
+          snadj += quadraticoffset;
           int ncount = 0;
           for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
             double bi = snaptr[tid]->bvec[icoeff];
@@ -321,21 +323,36 @@ void ComputeSNADAtom::compute_peratom()
             double biy = snaptr[tid]->dbvec[icoeff][1];
             double biz = snaptr[tid]->dbvec[icoeff][2];
 
+            // diagonal elements of quadratic matrix
+
+            double dbxtmp = bi*bix;
+            double dbytmp = bi*biy;
+            double dbztmp = bi*biz;
+
+            snadi[ncount] +=         dbxtmp;
+            snadi[ncount+yoffset] += dbytmp;
+            snadi[ncount+zoffset] += dbztmp;
+            snadj[ncount] -=         dbxtmp;
+            snadj[ncount+yoffset] -= dbytmp;
+            snadj[ncount+zoffset] -= dbztmp;
+            ncount++;
+
             // upper-triangular elements of quadratic matrix
 
-            for (int jcoeff = icoeff; jcoeff < ncoeff; jcoeff++) {
+            for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
               double dbxtmp = bi*snaptr[tid]->dbvec[jcoeff][0]
                 + bix*snaptr[tid]->bvec[jcoeff];
               double dbytmp = bi*snaptr[tid]->dbvec[jcoeff][1]
                 + biy*snaptr[tid]->bvec[jcoeff];
               double dbztmp = bi*snaptr[tid]->dbvec[jcoeff][2]
                 + biz*snaptr[tid]->bvec[jcoeff];
-              snadi[ncount] += dbxtmp;
-              snadi[ncount+ncoeffq] += dbytmp;
-              snadi[ncount+twoncoeffq] += dbztmp;
-              snadj[ncount] -= dbxtmp;
-              snadj[ncount+ncoeffq] -= dbytmp;
-              snadj[ncount+twoncoeffq] -= dbztmp;
+
+              snadi[ncount] +=         dbxtmp;
+              snadi[ncount+yoffset] += dbytmp;
+              snadi[ncount+zoffset] += dbztmp;
+              snadj[ncount] -=         dbxtmp;
+              snadj[ncount+yoffset] -= dbytmp;
+              snadj[ncount+zoffset] -= dbztmp;
               ncount++;
             }
           }
@@ -361,7 +378,7 @@ int ComputeSNADAtom::pack_reverse_comm(int n, int first, double *buf)
   for (i = first; i < last; i++)
     for (icoeff = 0; icoeff < size_peratom_cols; icoeff++)
       buf[m++] = snad[i][icoeff];
-  return comm_reverse;
+  return m;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -387,8 +404,7 @@ double ComputeSNADAtom::memory_usage()
   double bytes = nmax*size_peratom_cols * sizeof(double);
   bytes += 3*njmax*sizeof(double);
   bytes += njmax*sizeof(int);
-  bytes += threencoeff*atom->ntypes;
-  if (quadraticflag) bytes += threencoeffq*atom->ntypes;
+  bytes += 3*nperdim*atom->ntypes;
   bytes += snaptr[0]->memory_usage()*comm->nthreads;
   return bytes;
 }
diff --git a/src/SNAP/compute_snad_atom.h b/src/SNAP/compute_snad_atom.h
index 5e40632d8c..a33e6047c2 100644
--- a/src/SNAP/compute_snad_atom.h
+++ b/src/SNAP/compute_snad_atom.h
@@ -37,7 +37,7 @@ class ComputeSNADAtom : public Compute {
 
  private:
   int nmax, njmax, diagonalstyle;
-  int ncoeff, twoncoeff, threencoeff, ncoeffq, twoncoeffq, threencoeffq;
+  int ncoeff, nperdim, yoffset, zoffset;
   double **cutsq;
   class NeighList *list;
   double **snad;
diff --git a/src/SNAP/compute_snav_atom.cpp b/src/SNAP/compute_snav_atom.cpp
index 0278be5a97..1248a223c3 100644
--- a/src/SNAP/compute_snav_atom.cpp
+++ b/src/SNAP/compute_snav_atom.cpp
@@ -96,6 +96,11 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
         error->all(FLERR,"Illegal compute snav/atom command");
       switchflag = atoi(arg[iarg+1]);
       iarg += 2;
+    } else if (strcmp(arg[iarg],"bzeroflag") == 0) {
+      if (iarg+2 > narg)
+        error->all(FLERR,"Illegal compute snav/atom command");
+      bzeroflag = atoi(arg[iarg+1]);
+      iarg += 2;
     } else if (strcmp(arg[iarg],"quadraticflag") == 0) {
       if (iarg+2 > narg)
         error->all(FLERR,"Illegal compute snav/atom command");
@@ -117,22 +122,9 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
   }
 
   ncoeff = snaptr[0]->ncoeff;
-  twoncoeff = 2*ncoeff;
-  threencoeff = 3*ncoeff;
-  fourncoeff = 4*ncoeff;
-  fivencoeff = 5*ncoeff;
-  sixncoeff = 6*ncoeff;
-  size_peratom_cols = sixncoeff*atom->ntypes;
-  if (quadraticflag) {
-    ncoeffq = ncoeff*ncoeff;
-    twoncoeffq = 2*ncoeffq;
-    threencoeffq = 3*ncoeffq;
-    fourncoeffq = 4*ncoeffq;
-    fivencoeffq = 5*ncoeffq;
-    sixncoeffq = 6*ncoeffq;
-    size_peratom_cols +=
-      sixncoeffq*atom->ntypes;
-  }
+  nperdim = ncoeff;  
+  if (quadraticflag) nperdim += (ncoeff*(ncoeff+1))/2;
+  size_peratom_cols = 6*nperdim*atom->ntypes;
   comm_reverse = size_peratom_cols;
   peratom_flag = 1;
 
@@ -251,9 +243,7 @@ void ComputeSNAVAtom::compute_peratom()
       const int* const jlist = firstneigh[i];
       const int jnum = numneigh[i];
 
-      const int typeoffset = sixncoeff*(atom->type[i]-1);
-      const int quadraticoffset = sixncoeff*atom->ntypes +
-        sixncoeffq*(atom->type[i]-1);
+      const int typeoffset = 6*nperdim*(atom->type[i]-1);
 
       // insure rij, inside, and typej  are of size jnum
 
@@ -307,23 +297,24 @@ 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+ncoeff]      += snaptr[tid]->dbvec[icoeff][1]*ytmp;
-          snavi[icoeff+twoncoeff]   += snaptr[tid]->dbvec[icoeff][2]*ztmp;
-          snavi[icoeff+threencoeff] += snaptr[tid]->dbvec[icoeff][1]*ztmp;
-          snavi[icoeff+fourncoeff]  += snaptr[tid]->dbvec[icoeff][0]*ztmp;
-          snavi[icoeff+fivencoeff]  += snaptr[tid]->dbvec[icoeff][0]*ytmp;
-          snavj[icoeff]             -= snaptr[tid]->dbvec[icoeff][0]*x[j][0];
-          snavj[icoeff+ncoeff]      -= snaptr[tid]->dbvec[icoeff][1]*x[j][1];
-          snavj[icoeff+twoncoeff]   -= snaptr[tid]->dbvec[icoeff][2]*x[j][2];
-          snavj[icoeff+threencoeff] -= snaptr[tid]->dbvec[icoeff][1]*x[j][2];
-          snavj[icoeff+fourncoeff]  -= snaptr[tid]->dbvec[icoeff][0]*x[j][2];
-          snavj[icoeff+fivencoeff]  -= snaptr[tid]->dbvec[icoeff][0]*x[j][1];
+          snavi[icoeff]           += snaptr[tid]->dbvec[icoeff][0]*xtmp;
+          snavi[icoeff+nperdim]   += snaptr[tid]->dbvec[icoeff][1]*ytmp;
+          snavi[icoeff+2*nperdim] += snaptr[tid]->dbvec[icoeff][2]*ztmp;
+          snavi[icoeff+3*nperdim] += snaptr[tid]->dbvec[icoeff][1]*ztmp;
+          snavi[icoeff+4*nperdim] += snaptr[tid]->dbvec[icoeff][0]*ztmp;
+          snavi[icoeff+5*nperdim] += snaptr[tid]->dbvec[icoeff][0]*ytmp;
+          snavj[icoeff]           -= snaptr[tid]->dbvec[icoeff][0]*x[j][0];
+          snavj[icoeff+nperdim]   -= snaptr[tid]->dbvec[icoeff][1]*x[j][1];
+          snavj[icoeff+2*nperdim] -= snaptr[tid]->dbvec[icoeff][2]*x[j][2];
+          snavj[icoeff+3*nperdim] -= snaptr[tid]->dbvec[icoeff][1]*x[j][2];
+          snavj[icoeff+4*nperdim] -= snaptr[tid]->dbvec[icoeff][0]*x[j][2];
+          snavj[icoeff+5*nperdim] -= snaptr[tid]->dbvec[icoeff][0]*x[j][1];
         }
 
         if (quadraticflag) {
-          double *snavi = snav[i]+quadraticoffset;
-          double *snavj = snav[j]+quadraticoffset;
+          const int quadraticoffset = ncoeff;
+          snavi += quadraticoffset;
+          snavj += quadraticoffset;
           int ncount = 0;
           for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
             double bi = snaptr[tid]->bvec[icoeff];
@@ -331,27 +322,46 @@ void ComputeSNAVAtom::compute_peratom()
             double biy = snaptr[tid]->dbvec[icoeff][1];
             double biz = snaptr[tid]->dbvec[icoeff][2];
 
+            // diagonal element of quadratic matrix
+            
+            double dbxtmp = bi*bix;
+            double dbytmp = bi*biy;
+            double dbztmp = bi*biz;
+            snavi[ncount] +=           dbxtmp*xtmp;
+            snavi[ncount+nperdim] +=   dbytmp*ytmp;
+            snavi[ncount+2*nperdim] += dbztmp*ztmp;
+            snavi[ncount+3*nperdim] += dbytmp*ztmp;
+            snavi[ncount+4*nperdim] += dbxtmp*ztmp;
+            snavi[ncount+5*nperdim] += dbxtmp*ytmp;
+            snavj[ncount] -=            dbxtmp*x[j][0];
+            snavj[ncount+nperdim] -=    dbytmp*x[j][1];
+            snavj[ncount+2*nperdim] -=  dbztmp*x[j][2];
+            snavj[ncount+3*nperdim] -=  dbytmp*x[j][2];
+            snavj[ncount+4*nperdim] -=  dbxtmp*x[j][2];
+            snavj[ncount+5*nperdim] -=  dbxtmp*x[j][1];
+            ncount++;
+
             // upper-triangular elements of quadratic matrix
 
-            for (int jcoeff = icoeff; jcoeff < ncoeff; jcoeff++) {
+            for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
               double dbxtmp = bi*snaptr[tid]->dbvec[jcoeff][0]
                 + bix*snaptr[tid]->bvec[jcoeff];
               double dbytmp = bi*snaptr[tid]->dbvec[jcoeff][1]
                 + biy*snaptr[tid]->bvec[jcoeff];
               double dbztmp = bi*snaptr[tid]->dbvec[jcoeff][2]
                 + biz*snaptr[tid]->bvec[jcoeff];
-              snavi[ncount] +=               dbxtmp*xtmp;
-              snavi[ncount+ncoeffq] +=      dbytmp*ytmp;
-              snavi[ncount+twoncoeffq] +=   dbztmp*ztmp;
-              snavi[ncount+threencoeffq] += dbytmp*ztmp;
-              snavi[ncount+fourncoeffq] +=  dbxtmp*ztmp;
-              snavi[ncount+fivencoeffq] +=  dbxtmp*ytmp;
-              snavj[ncount] -=               dbxtmp*x[j][0];
-              snavj[ncount+ncoeffq] -=      dbytmp*x[j][1];
-              snavj[ncount+twoncoeffq] -=   dbztmp*x[j][2];
-              snavj[ncount+threencoeffq] -= dbytmp*x[j][2];
-              snavj[ncount+fourncoeffq] -=  dbxtmp*x[j][2];
-              snavj[ncount+fivencoeffq] -=  dbxtmp*x[j][1];
+              snavi[ncount] +=           dbxtmp*xtmp;
+              snavi[ncount+nperdim] +=   dbytmp*ytmp;
+              snavi[ncount+2*nperdim] += dbztmp*ztmp;
+              snavi[ncount+3*nperdim] += dbytmp*ztmp;
+              snavi[ncount+4*nperdim] += dbxtmp*ztmp;
+              snavi[ncount+5*nperdim] += dbxtmp*ytmp;
+              snavj[ncount] -=           dbxtmp*x[j][0];
+              snavj[ncount+nperdim] -=   dbytmp*x[j][1];
+              snavj[ncount+2*nperdim] -= dbztmp*x[j][2];
+              snavj[ncount+3*nperdim] -= dbytmp*x[j][2];
+              snavj[ncount+4*nperdim] -= dbxtmp*x[j][2];
+              snavj[ncount+5*nperdim] -= dbxtmp*x[j][1];
               ncount++;
             }
           }
@@ -377,7 +387,7 @@ int ComputeSNAVAtom::pack_reverse_comm(int n, int first, double *buf)
   for (i = first; i < last; i++)
     for (icoeff = 0; icoeff < size_peratom_cols; icoeff++)
       buf[m++] = snav[i][icoeff];
-  return comm_reverse;
+  return m;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -403,8 +413,8 @@ double ComputeSNAVAtom::memory_usage()
   double bytes = nmax*size_peratom_cols * sizeof(double);
   bytes += 3*njmax*sizeof(double);
   bytes += njmax*sizeof(int);
-  bytes += sixncoeff*atom->ntypes;
-  if (quadraticflag) bytes += sixncoeffq*atom->ntypes;
+  bytes += 6*nperdim*atom->ntypes;
+  if (quadraticflag) bytes += 6*nperdim*atom->ntypes;
   bytes += snaptr[0]->memory_usage()*comm->nthreads;
   return bytes;
 }
diff --git a/src/SNAP/compute_snav_atom.h b/src/SNAP/compute_snav_atom.h
index 35f1478393..2eb9fb804f 100644
--- a/src/SNAP/compute_snav_atom.h
+++ b/src/SNAP/compute_snav_atom.h
@@ -37,8 +37,7 @@ class ComputeSNAVAtom : public Compute {
 
  private:
   int nmax, njmax, diagonalstyle;
-  int ncoeff, twoncoeff, threencoeff, fourncoeff, fivencoeff, sixncoeff;
-  int ncoeffq, twoncoeffq, threencoeffq, fourncoeffq, fivencoeffq, sixncoeffq;
+  int ncoeff, nperdim;
   double **cutsq;
   class NeighList *list;
   double **snav;
diff --git a/src/SNAP/pair_snap.cpp b/src/SNAP/pair_snap.cpp
index 377235685c..8714d55c5c 100644
--- a/src/SNAP/pair_snap.cpp
+++ b/src/SNAP/pair_snap.cpp
@@ -278,14 +278,15 @@ void PairSNAP::compute_regular(int eflag, int vflag)
           double bveci = snaptr->bvec[icoeff];
           double fack = coeffi[k]*bveci;
           double* dbveci = snaptr->dbvec[icoeff];
-          fij[0] += fack*snaptr->dbvec[icoeff][0];
-          fij[1] += fack*snaptr->dbvec[icoeff][1];
-          fij[2] += fack*snaptr->dbvec[icoeff][2];
+          fij[0] += fack*dbveci[0];
+          fij[1] += fack*dbveci[1];
+          fij[2] += fack*dbveci[2];
           k++;
           for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
             double facki = coeffi[k]*bveci;
             double fackj = coeffi[k]*snaptr->bvec[jcoeff];
             double* dbvecj = snaptr->dbvec[jcoeff];
+
             fij[0] += facki*dbvecj[0]+fackj*dbveci[0];
             fij[1] += facki*dbvecj[1]+fackj*dbveci[1];
             fij[2] += facki*dbvecj[2]+fackj*dbveci[2];
@@ -1529,10 +1530,10 @@ void PairSNAP::coeff(int narg, char **arg)
   }
 
   if (comm->me == 0)
-    printf("ncoeff = %d snancoeff = %d \n",ncoeff,sna[0]->ncoeff);
-  if (ncoeff != sna[0]->ncoeff) {
-    error->all(FLERR,"Incorrect SNAP parameter file");
-  }
+    if (ncoeff != sna[0]->ncoeff) {
+      printf("ncoeff = %d snancoeff = %d \n",ncoeff,sna[0]->ncoeff);
+      error->all(FLERR,"Incorrect SNAP parameter file");
+    }
 
   // Calculate maximum cutoff for all elements
 
diff --git a/src/SNAP/pair_snap.h b/src/SNAP/pair_snap.h
index d39cb0f8d4..b60ab3c3e4 100644
--- a/src/SNAP/pair_snap.h
+++ b/src/SNAP/pair_snap.h
@@ -37,11 +37,8 @@ public:
   virtual double init_one(int, int);
   virtual double memory_usage();
 
-  double rcutfac, quadraticflag; // declared public to workaround gcc 4.9
-  int ncoeff;                    //  compiler bug, manifest in KOKKOS package
-
 protected:
-  int ncoeffq, ncoeffall;
+  int ncoeff, ncoeffq, ncoeffall;
   double **bvec, ***dbvec;
   class SNA** sna;
   int nmax;
@@ -100,8 +97,8 @@ protected:
   double *wjelem;               // elements weights
   double **coeffelem;           // element bispectrum coefficients
   int *map;                     // mapping from atom types to elements
-  int twojmax, diagonalstyle, switchflag, bzeroflag;
-  double rfac0, rmin0, wj1, wj2;
+  int twojmax, diagonalstyle, switchflag, bzeroflag, quadraticflag;
+  double rcutfac, rfac0, rmin0, wj1, wj2;
   int rcutfacflag, twojmaxflag; // flags for required parameters
 };
 
-- 
GitLab