SNAP changes by Aidan

doc/src/compute_sna_atom.txt

@@ -231,11 +231,12 @@ 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
-a matrix of quantities that are products of two bispectrum components. If the
+the products of all distinct pairs of  bispectrum components. If the
-number of bispectrum components is {K}, then the number of matrix elements
+number of bispectrum components is {K}, then the number of distinct pairs
-is {K}^2. These are output in subblocks of {K}^2 columns, using the same
+is  {K}({K}+1)/2. These are output in subblocks of  {K}({K}+1)/2 columns, using the same
-ordering of columns and sub-blocks as was used for the bispectrum
+ordering of sub-blocks as was used for the bispectrum
-components.
+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})
 
 These values can be accessed by any command that uses per-atom values
 from a compute as input.  See "Section

doc/src/pair_snap.txt

@@ -139,10 +139,15 @@ The default values for these keywords are
 {rmin0} = 0.0
 {diagonalstyle} = 3
 {switchflag} = 0
-{bzeroflag} = 1 :ul
+{bzeroflag} = 1
+{quadraticflag} = 1 :ul
 
-Detailed definitions of these keywords are given on the "compute
+Detailed definitions for all the keywords are given on the "compute
 sna/atom"_compute_sna_atom.html doc page.
+If {quadraticflag} is set to 1, then the SNAP energy expression includes the quadratic term,
+0.5*B^t.alpha.B, where alpha is a symmetric {K} by {K} matrix.
+The SNAP element file should contain {K}({K}+1)/2 additional coefficients
+for each element, the upper-triangular elements of alpha.
 
 :line
 

examples/COUPLE/simple/simple.cpp

@@ -153,7 +153,7 @@ int main(int narg, char **arg)
     for (int i = 0; i < natoms; i++) type[i] = 1;
 
     lmp->input->one("delete_atoms group all");
-    lammps_create_atoms(lmp,natoms,NULL,type,x,v);
+    lammps_create_atoms(lmp,natoms,NULL,type,x,v,NULL,0);
     lmp->input->one("run 10");
   }
 

examples/snap/Ta06A.snapparam

@@ -8,8 +8,8 @@ twojmax 6
 
 # optional
 
-gamma 1
 rfac0 0.99363
 rmin0 0
 diagonalstyle 3
 bzeroflag 0
+quadraticflag 0

examples/snap/W_2940_2017_2.snap

@@ -5,7 +5,7 @@ variable zblcutinner equal 4
 variable zblcutouter equal 4.8
 variable zblz equal 74
 
-# Specify hybrid with SNAP, ZBL, and long-range Coulomb
+# Specify hybrid with SNAP and ZBL
 
 pair_style hybrid/overlay &
 zbl ${zblcutinner} ${zblcutouter} &

examples/snap/W_2940_2017_2.snapparam

@@ -6,8 +6,8 @@ twojmax 8
 
 # optional
 
-gamma 1
 rfac0 0.99363
 rmin0 0
 diagonalstyle 3
 bzeroflag 0
+quadraticflag 0

examples/snap/W_2940_2017_2_He_JW2013.snap

@@ -5,7 +5,7 @@ variable zblcutinner equal 4
 variable zblcutouter equal 4.8
 variable zblz equal 74
 
-# Specify hybrid with SNAP, ZBL, and long-range Coulomb
+# Specify hybrid with SNAP and ZBL
 
 pair_style hybrid/overlay zbl ${zblcutinner} ${zblcutouter} snap table spline 10000 table spline 10000
 pair_coeff 	1 1 zbl ${zblz} ${zblz}

src/SNAP/compute_sna_atom.cpp

@@ -129,7 +129,7 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
 
   ncoeff = snaptr[0]->ncoeff;
   size_peratom_cols = ncoeff;
-  if (quadraticflag) size_peratom_cols += ncoeff*ncoeff;
+  if (quadraticflag) size_peratom_cols += (ncoeff*(ncoeff+1))/2;
   peratom_flag = 1;
 
   nmax = 0;
@@ -275,7 +275,10 @@ void ComputeSNAAtom::compute_peratom()
         int ncount = ncoeff;
         for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
           double bi = snaptr[tid]->bvec[icoeff];
-          for (int jcoeff = 0; jcoeff < ncoeff; jcoeff++)
+
+          // upper-triangular elements of quadratic matrix
+          
+          for (int jcoeff = icoeff; jcoeff < ncoeff; jcoeff++)
             sna[i][ncount++] = bi*snaptr[tid]->bvec[jcoeff];
         }
       }

src/SNAP/compute_snad_atom.cpp

@@ -125,11 +125,11 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
   threencoeff = 3*ncoeff;
   size_peratom_cols = threencoeff*atom->ntypes;
   if (quadraticflag) {
-    ncoeffsq = ncoeff*ncoeff;
+    ncoeffq = (ncoeff*(ncoeff+1))/2;
-    twoncoeffsq = 2*ncoeffsq;
+    twoncoeffq = 2*ncoeffq;
-    threencoeffsq = 3*ncoeffsq;
+    threencoeffq = 3*ncoeffq;
     size_peratom_cols +=
-      threencoeffsq*atom->ntypes;
+      threencoeffq*atom->ntypes;
   }
   comm_reverse = size_peratom_cols;
   peratom_flag = 1;
@@ -250,7 +250,7 @@ void ComputeSNADAtom::compute_peratom()
 
       const int typeoffset = threencoeff*(atom->type[i]-1);
       const int quadraticoffset = threencoeff*atom->ntypes +
-        threencoeffsq*(atom->type[i]-1);
+        threencoeffq*(atom->type[i]-1);
 
       // insure rij, inside, and typej  are of size jnum
 
@@ -320,7 +320,10 @@ void ComputeSNADAtom::compute_peratom()
             double bix = snaptr[tid]->dbvec[icoeff][0];
             double biy = snaptr[tid]->dbvec[icoeff][1];
             double biz = snaptr[tid]->dbvec[icoeff][2];
-            for (int jcoeff = 0; jcoeff < ncoeff; jcoeff++) {
+
+            // upper-triangular elements of quadratic matrix
+          
+            for (int jcoeff = icoeff; jcoeff < ncoeff; jcoeff++) {
               double dbxtmp = bi*snaptr[tid]->dbvec[jcoeff][0]
                 + bix*snaptr[tid]->bvec[jcoeff];
               double dbytmp = bi*snaptr[tid]->dbvec[jcoeff][1]
@@ -328,11 +331,11 @@ void ComputeSNADAtom::compute_peratom()
               double dbztmp = bi*snaptr[tid]->dbvec[jcoeff][2]
                 + biz*snaptr[tid]->bvec[jcoeff];
               snadi[ncount] += dbxtmp;
-              snadi[ncount+ncoeffsq] += dbytmp;
+              snadi[ncount+ncoeffq] += dbytmp;
-              snadi[ncount+twoncoeffsq] += dbztmp;
+              snadi[ncount+twoncoeffq] += dbztmp;
               snadj[ncount] -= dbxtmp;
-              snadj[ncount+ncoeffsq] -= dbytmp;
+              snadj[ncount+ncoeffq] -= dbytmp;
-              snadj[ncount+twoncoeffsq] -= dbztmp;
+              snadj[ncount+twoncoeffq] -= dbztmp;
               ncount++;
             }
           }
@@ -385,7 +388,7 @@ double ComputeSNADAtom::memory_usage()
   bytes += 3*njmax*sizeof(double);
   bytes += njmax*sizeof(int);
   bytes += threencoeff*atom->ntypes;
-  if (quadraticflag) bytes += threencoeffsq*atom->ntypes;
+  if (quadraticflag) bytes += threencoeffq*atom->ntypes;
   bytes += snaptr[0]->memory_usage()*comm->nthreads;
   return bytes;
 }

src/SNAP/compute_snad_atom.h

@@ -37,7 +37,7 @@ class ComputeSNADAtom : public Compute {
 
  private:
   int nmax, njmax, diagonalstyle;
-  int ncoeff, twoncoeff, threencoeff, ncoeffsq, twoncoeffsq, threencoeffsq;
+  int ncoeff, twoncoeff, threencoeff, ncoeffq, twoncoeffq, threencoeffq;
   double **cutsq;
   class NeighList *list;
   double **snad;

src/SNAP/compute_snav_atom.cpp

@@ -124,14 +124,14 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
   sixncoeff = 6*ncoeff;
   size_peratom_cols = sixncoeff*atom->ntypes;
   if (quadraticflag) {
-    ncoeffsq = ncoeff*ncoeff;
+    ncoeffq = ncoeff*ncoeff;
-    twoncoeffsq = 2*ncoeffsq;
+    twoncoeffq = 2*ncoeffq;
-    threencoeffsq = 3*ncoeffsq;
+    threencoeffq = 3*ncoeffq;
-    fourncoeffsq = 4*ncoeffsq;
+    fourncoeffq = 4*ncoeffq;
-    fivencoeffsq = 5*ncoeffsq;
+    fivencoeffq = 5*ncoeffq;
-    sixncoeffsq = 6*ncoeffsq;
+    sixncoeffq = 6*ncoeffq;
     size_peratom_cols +=
-      sixncoeffsq*atom->ntypes;
+      sixncoeffq*atom->ntypes;
   }
   comm_reverse = size_peratom_cols;
   peratom_flag = 1;
@@ -253,7 +253,7 @@ void ComputeSNAVAtom::compute_peratom()
 
       const int typeoffset = sixncoeff*(atom->type[i]-1);
       const int quadraticoffset = sixncoeff*atom->ntypes +
-        sixncoeffsq*(atom->type[i]-1);
+        sixncoeffq*(atom->type[i]-1);
 
       // insure rij, inside, and typej  are of size jnum
 
@@ -330,7 +330,10 @@ void ComputeSNAVAtom::compute_peratom()
             double bix = snaptr[tid]->dbvec[icoeff][0];
             double biy = snaptr[tid]->dbvec[icoeff][1];
             double biz = snaptr[tid]->dbvec[icoeff][2];
-            for (int jcoeff = 0; jcoeff < ncoeff; jcoeff++) {
+
+            // upper-triangular elements of quadratic matrix
+          
+            for (int jcoeff = icoeff; jcoeff < ncoeff; jcoeff++) {
               double dbxtmp = bi*snaptr[tid]->dbvec[jcoeff][0]
                 + bix*snaptr[tid]->bvec[jcoeff];
               double dbytmp = bi*snaptr[tid]->dbvec[jcoeff][1]
@@ -338,17 +341,17 @@ void ComputeSNAVAtom::compute_peratom()
               double dbztmp = bi*snaptr[tid]->dbvec[jcoeff][2]
                 + biz*snaptr[tid]->bvec[jcoeff];
               snavi[ncount] +=               dbxtmp*xtmp;
-              snavi[ncount+ncoeffsq] +=      dbytmp*ytmp;
+              snavi[ncount+ncoeffq] +=      dbytmp*ytmp;
-              snavi[ncount+twoncoeffsq] +=   dbztmp*ztmp;
+              snavi[ncount+twoncoeffq] +=   dbztmp*ztmp;
-              snavi[ncount+threencoeffsq] += dbytmp*ztmp;
+              snavi[ncount+threencoeffq] += dbytmp*ztmp;
-              snavi[ncount+fourncoeffsq] +=  dbxtmp*ztmp;
+              snavi[ncount+fourncoeffq] +=  dbxtmp*ztmp;
-              snavi[ncount+fivencoeffsq] +=  dbxtmp*ytmp;
+              snavi[ncount+fivencoeffq] +=  dbxtmp*ytmp;
               snavj[ncount] -=               dbxtmp*x[j][0];
-              snavj[ncount+ncoeffsq] -=      dbytmp*x[j][1];
+              snavj[ncount+ncoeffq] -=      dbytmp*x[j][1];
-              snavj[ncount+twoncoeffsq] -=   dbztmp*x[j][2];
+              snavj[ncount+twoncoeffq] -=   dbztmp*x[j][2];
-              snavj[ncount+threencoeffsq] -= dbytmp*x[j][2];
+              snavj[ncount+threencoeffq] -= dbytmp*x[j][2];
-              snavj[ncount+fourncoeffsq] -=  dbxtmp*x[j][2];
+              snavj[ncount+fourncoeffq] -=  dbxtmp*x[j][2];
-              snavj[ncount+fivencoeffsq] -=  dbxtmp*x[j][1];
+              snavj[ncount+fivencoeffq] -=  dbxtmp*x[j][1];
               ncount++;
             }
           }
@@ -401,7 +404,7 @@ double ComputeSNAVAtom::memory_usage()
   bytes += 3*njmax*sizeof(double);
   bytes += njmax*sizeof(int);
   bytes += sixncoeff*atom->ntypes;
-  if (quadraticflag) bytes += sixncoeffsq*atom->ntypes;
+  if (quadraticflag) bytes += sixncoeffq*atom->ntypes;
   bytes += snaptr[0]->memory_usage()*comm->nthreads;
   return bytes;
 }

src/SNAP/compute_snav_atom.h

@@ -38,7 +38,7 @@ class ComputeSNAVAtom : public Compute {
  private:
   int nmax, njmax, diagonalstyle;
   int ncoeff, twoncoeff, threencoeff, fourncoeff, fivencoeff, sixncoeff;
-  int ncoeffsq, twoncoeffsq, threencoeffsq, fourncoeffsq, fivencoeffsq, sixncoeffsq;
+  int ncoeffq, twoncoeffq, threencoeffq, fourncoeffq, fivencoeffq, sixncoeffq;
   double **cutsq;
   class NeighList *list;
   double **snav;

src/SNAP/pair_snap.cpp

@@ -35,6 +35,10 @@ using namespace LAMMPS_NS;
 #define MAXLINE 1024
 #define MAXWORD 3
 
+// Outstanding issues with quadratic term
+// 1. there seems to a problem with compute_optimized energy calc
+// it does not match compute_regular, even when quadratic coeffs = 0
+
 /* ---------------------------------------------------------------------- */
 
 PairSNAP::PairSNAP(LAMMPS *lmp) : Pair(lmp)
@@ -232,10 +236,6 @@ void PairSNAP::compute_regular(int eflag, int vflag)
 
     snaptr->compute_ui(ninside);
     snaptr->compute_zi();
-    if (!gammaoneflag) {
-      snaptr->compute_bi();
-      snaptr->copy_bi2bvec();
-    }
 
     // for neighbors of I within cutoff:
     // compute dUi/drj and dBi/drj
@@ -255,17 +255,41 @@ void PairSNAP::compute_regular(int eflag, int vflag)
       fij[1] = 0.0;
       fij[2] = 0.0;
 
+      // linear contributions
+      
       for (int k = 1; k <= ncoeff; k++) {
-	double bgb;
+	double bgb = coeffi[k];
-	if (gammaoneflag)
-	  bgb = coeffi[k];
-	else bgb = coeffi[k]*
-	       gamma*pow(snaptr->bvec[k-1],gamma-1.0);
 	fij[0] += bgb*snaptr->dbvec[k-1][0];
 	fij[1] += bgb*snaptr->dbvec[k-1][1];
 	fij[2] += bgb*snaptr->dbvec[k-1][2];
       }
 
+      // quadratic contributions
+      
+      if (quadraticflag) {
+        snaptr->compute_bi();
+        snaptr->copy_bi2bvec();
+        int k = ncoeff+1;
+        for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+          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];
+          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];
+            k++;
+          }
+        }
+      }
+      
       f[i][0] += fij[0];
       f[i][1] += fij[1];
       f[i][2] += fij[2];
@@ -285,14 +309,33 @@ void PairSNAP::compute_regular(int eflag, int vflag)
       // evdwl = energy of atom I, sum over coeffs_k * Bi_k
 
       evdwl = coeffi[0];
-      if (gammaoneflag) {
+      if (!quadraticflag) {
-	snaptr->compute_bi();
+        snaptr->compute_bi();
-	snaptr->copy_bi2bvec();
+        snaptr->copy_bi2bvec();
-	for (int k = 1; k <= ncoeff; k++)
+      }
-	  evdwl += coeffi[k]*snaptr->bvec[k-1];
+      
-      } else
+      // E = beta.B + 0.5*B^t.alpha.B
-      	for (int k = 1; k <= ncoeff; k++)
+      // coeff[k] = beta[k-1] or
-      	  evdwl += coeffi[k]*pow(snaptr->bvec[k-1],gamma);
+      // coeff[k] = alpha_ii or
+      // coeff[k] = alpha_ij = alpha_ji, j != i
+
+      // linear contributions
+
+      for (int k = 1; k <= ncoeff; k++)
+        evdwl += coeffi[k]*snaptr->bvec[k-1];
+
+      // quadratic contributions
+
+      if (quadraticflag) {
+        int k = ncoeff+1;
+        for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+          double bveci = snaptr->bvec[icoeff];
+          evdwl += 0.5*coeffi[k++]*bveci*bveci;
+          for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
+            evdwl += coeffi[k++]*bveci*snaptr->bvec[jcoeff];
+          }
+        }
+      }
       ev_tally_full(i,2.0*evdwl,0.0,0.0,delx,dely,delz);
     }
 
@@ -562,26 +605,46 @@ void PairSNAP::compute_optimized(int eflag, int vflag)
 
         sna[tid]->compute_dbidrj();
         sna[tid]->copy_dbi2dbvec();
-	if (!gammaoneflag) {
-	  sna[tid]->compute_bi();
-	  sna[tid]->copy_bi2bvec();
-	}
 
         fij[0] = 0.0;
         fij[1] = 0.0;
         fij[2] = 0.0;
 
+        // linear contributions
+
         for (k = 1; k <= ncoeff; k++) {
-	  double bgb;
+	  double bgb = coeffi[k];
-	  if (gammaoneflag)
-	    bgb = coeffi[k];
-	  else bgb = coeffi[k]*
-		 gamma*pow(sna[tid]->bvec[k-1],gamma-1.0);
 	  fij[0] += bgb*sna[tid]->dbvec[k-1][0];
 	  fij[1] += bgb*sna[tid]->dbvec[k-1][1];
 	  fij[2] += bgb*sna[tid]->dbvec[k-1][2];
         }
 
+      // quadratic contributions
+        
+      if (quadraticflag) {
+        sna[tid]->compute_bi();
+        sna[tid]->copy_bi2bvec();
+        int k = ncoeff+1;
+        for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+          double bveci = sna[tid]->bvec[icoeff];
+          double fack = coeffi[k]*bveci;
+          double* dbveci = sna[tid]->dbvec[icoeff];
+          fij[0] += fack*sna[tid]->dbvec[icoeff][0];
+          fij[1] += fack*sna[tid]->dbvec[icoeff][1];
+          fij[2] += fack*sna[tid]->dbvec[icoeff][2];
+          k++;
+          for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
+            double facki = coeffi[k]*bveci;
+            double fackj = coeffi[k]*sna[tid]->bvec[jcoeff];
+            double* dbvecj = sna[tid]->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];
+            k++;
+          }
+        }
+      }
+
 #if defined(_OPENMP)
 #pragma omp critical
 #endif
@@ -606,15 +669,33 @@ void PairSNAP::compute_optimized(int eflag, int vflag)
 
       if (eflag&&pairs[iijj][1] == 0) {
 	evdwl = coeffi[0];
-	if (gammaoneflag) {
-	  sna[tid]->compute_bi();
-	  sna[tid]->copy_bi2bvec();
-	  for (int k = 1; k <= ncoeff; k++)
-	    evdwl += coeffi[k]*sna[tid]->bvec[k-1];
-	} else
-	  for (int k = 1; k <= ncoeff; k++)
-	    evdwl += coeffi[k]*pow(sna[tid]->bvec[k-1],gamma);
 
+        sna[tid]->compute_bi();
+        sna[tid]->copy_bi2bvec();
+
+        // E = beta.B + 0.5*B^t.alpha.B
+        // coeff[k] = beta[k-1] or
+        // coeff[k] = alpha_ii or
+        // coeff[k] = alpha_ij = alpha_ji, j != i
+
+        // linear contributions
+        
+        for (int k = 1; k <= ncoeff; k++)
+          evdwl += coeffi[k]*sna[tid]->bvec[k-1];
+
+        // quadratic contributions
+        
+        if (quadraticflag) {
+          int k = ncoeff+1;
+          for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+            double bveci = sna[tid]->bvec[icoeff];
+            evdwl += 0.5*coeffi[k++]*bveci*bveci;
+            for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
+              evdwl += coeffi[k++]*bveci*sna[tid]->bvec[jcoeff];
+            }
+          }
+        }
+        
 #if defined(_OPENMP)
 #pragma omp critical
 #endif
@@ -1363,6 +1444,22 @@ void PairSNAP::coeff(int narg, char **arg)
 
   read_files(coefffilename,paramfilename);
 
+  if (!quadraticflag)
+    ncoeff = ncoeffall - 1;
+  else {
+
+    // ncoeffall should be (ncoeff+2)*(ncoeff+1)/2
+    // so, ncoeff = floor(sqrt(2*ncoeffall))-1
+    
+    ncoeff = sqrt(2*ncoeffall)-1; 
+    ncoeffq = (ncoeff*(ncoeff+1))/2;
+    int ntmp = 1+ncoeff+ncoeffq;
+    if (ntmp != ncoeffall) {
+      printf("ncoeffall = %d ntmp = %d ncoeff = %d \n",ncoeffall,ntmp,ncoeff);
+      error->all(FLERR,"Incorrect SNAP coeff file");
+    }
+  }
+
   // read args that map atom types to SNAP elements
   // map[i] = which element the Ith atom type is, -1 if not mapped
   // map[0] is not used
@@ -1416,6 +1513,7 @@ void PairSNAP::coeff(int narg, char **arg)
       sna[tid]->grow_rij(nmax);
   }
 
+  printf("ncoeff = %d snancoeff = %d \n",ncoeff,sna[0]->ncoeff);
   if (ncoeff != sna[0]->ncoeff) {
     printf("ncoeff = %d snancoeff = %d \n",ncoeff,sna[0]->ncoeff);
     error->all(FLERR,"Incorrect SNAP parameter file");
@@ -1470,7 +1568,7 @@ double PairSNAP::init_one(int i, int j)
 void PairSNAP::read_files(char *coefffilename, char *paramfilename)
 {
 
-  // open SNAP ceofficient file on proc 0
+  // open SNAP coefficient file on proc 0
 
   FILE *fpcoeff;
   if (comm->me == 0) {
@@ -1518,13 +1616,13 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
   words[iword] = strtok(NULL,"' \t\n\r\f");
 
   int nelemfile = atoi(words[0]);
-  ncoeff = atoi(words[1])-1;
+  ncoeffall = atoi(words[1]);
-
+  
   // Set up element lists
 
   memory->create(radelem,nelements,"pair:radelem");
   memory->create(wjelem,nelements,"pair:wjelem");
-  memory->create(coeffelem,nelements,ncoeff+1,"pair:coeffelem");
+  memory->create(coeffelem,nelements,ncoeffall,"pair:coeffelem");
 
   int *found = new int[nelements];
   for (int ielem = 0; ielem < nelements; ielem++)
@@ -1569,7 +1667,7 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
       if (strcmp(elemtmp,elements[ielem]) == 0) break;
     if (ielem == nelements) {
       if (comm->me == 0)
-	for (int icoeff = 0; icoeff <= ncoeff; icoeff++)
+	for (int icoeff = 0; icoeff < ncoeffall; icoeff++)
 	  ptr = fgets(line,MAXLINE,fpcoeff);
       continue;
     }
@@ -1578,7 +1676,7 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
 
     if (found[ielem]) {
       if (comm->me == 0)
-	for (int icoeff = 0; icoeff <= ncoeff; icoeff++)
+	for (int icoeff = 0; icoeff < ncoeffall; icoeff++)
 	  ptr = fgets(line,MAXLINE,fpcoeff);
       continue;
     }
@@ -1595,7 +1693,7 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
 			  elements[ielem], radelem[ielem], wjelem[ielem]);
     }
 
-    for (int icoeff = 0; icoeff <= ncoeff; icoeff++) {
+    for (int icoeff = 0; icoeff < ncoeffall; icoeff++) {
       if (comm->me == 0) {
 	ptr = fgets(line,MAXLINE,fpcoeff);
 	if (ptr == NULL) {
@@ -1629,13 +1727,12 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
 
   // Set defaults for optional keywords
 
-  gamma = 1.0;
-  gammaoneflag = 1;
   rfac0 = 0.99363;
   rmin0 = 0.0;
   diagonalstyle = 3;
   switchflag = 1;
   bzeroflag = 1;
+  quadraticflag = 0;
   
   // open SNAP parameter file on proc 0
 
@@ -1689,9 +1786,7 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
     } else if (strcmp(keywd,"twojmax") == 0) {
       twojmax = atoi(keyval);
       twojmaxflag = 1;
-    } else if (strcmp(keywd,"gamma") == 0)
+    } else if (strcmp(keywd,"rfac0") == 0)
-      gamma = atof(keyval);
-    else if (strcmp(keywd,"rfac0") == 0)
       rfac0 = atof(keyval);
     else if (strcmp(keywd,"rmin0") == 0)
       rmin0 = atof(keyval);
@@ -1701,6 +1796,8 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
       switchflag = atoi(keyval);
     else if (strcmp(keywd,"bzeroflag") == 0)
       bzeroflag = atoi(keyval);
+    else if (strcmp(keywd,"quadraticflag") == 0)
+      quadraticflag = atoi(keyval);
     else
       error->all(FLERR,"Incorrect SNAP parameter file");
   }
@@ -1708,9 +1805,6 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
   if (rcutfacflag == 0 || twojmaxflag == 0)
     error->all(FLERR,"Incorrect SNAP parameter file");
 
-  if (gamma == 1.0) gammaoneflag = 1;
-  else gammaoneflag = 0;
-
   delete[] found;
 }
 
@@ -1726,7 +1820,7 @@ double PairSNAP::memory_usage()
   bytes += n*n*sizeof(double);
   bytes += 3*nmax*sizeof(double);
   bytes += nmax*sizeof(int);
-  bytes += (2*ncoeff+1)*sizeof(double);
+  bytes += (2*ncoeffall)*sizeof(double);
   bytes += (ncoeff*3)*sizeof(double);
   bytes += sna[0]->memory_usage()*nthreads;
   return bytes;

src/SNAP/pair_snap.h

@@ -38,7 +38,7 @@ public:
   double memory_usage();
 
 protected:
-  int ncoeff;
+  int ncoeff, ncoeffq, ncoeffall;
   double **bvec, ***dbvec;
   class SNA** sna;
   int nmax;
@@ -89,7 +89,6 @@ protected:
   //  timespec starttime, endtime;
   double timers[4];
 #endif
-  double gamma;
 
   double rcutmax;               // max cutoff for all elements
   int nelements;                // # of unique elements
@@ -98,10 +97,9 @@ protected:
   double *wjelem;               // elements weights
   double **coeffelem;           // element bispectrum coefficients
   int *map;                     // mapping from atom types to elements
-  int twojmax, diagonalstyle, switchflag, bzeroflag;
+  int twojmax, diagonalstyle, switchflag, bzeroflag, quadraticflag;
   double rcutfac, rfac0, rmin0, wj1, wj2;
   int rcutfacflag, twojmaxflag; // flags for required parameters
-  int gammaoneflag;              // 1 if parameter gamma is 1
 };
 
 }

src/npair_full_bin_atomonly.h

@@ -14,7 +14,7 @@
 #ifdef NPAIR_CLASS
 
 NPairStyle(full/bin/atomonly,
-           NPairFullBin,
+           NPairFullBinAtomonly,
            NP_FULL | NP_BIN | NP_ATOMONLY |
 	   NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI)