From e98af882422074770d696c0fa14e389da3ce6fa4 Mon Sep 17 00:00:00 2001
From: rohskopf <drew.rohskopf@gmail.com>
Date: Fri, 1 Jul 2022 20:56:11 -0600
Subject: [PATCH] Clean up unnecessary code.

---
 examples/snap/compute_snap_dgrad.py |   2 +-
 src/ML-SNAP/compute_snap.cpp        | 447 ++++++----------------------
 src/ML-SNAP/compute_snap.h          |   9 +-
 3 files changed, 95 insertions(+), 363 deletions(-)

diff --git a/examples/snap/compute_snap_dgrad.py b/examples/snap/compute_snap_dgrad.py
index 95f2e1fa55..0919af8598 100644
--- a/examples/snap/compute_snap_dgrad.py
+++ b/examples/snap/compute_snap_dgrad.py
@@ -6,7 +6,7 @@ Purpose: Demonstrate extraction of descriptor gradient (dB/dR) array from comput
 Serial syntax:
     python compute_snap_dgrad.py
 Parallel syntax:
-    mpirun -np 2 python compute_snap_dgrad.py
+    mpirun -np 4 python compute_snap_dgrad.py
 """
 
 from __future__ import print_function
diff --git a/src/ML-SNAP/compute_snap.cpp b/src/ML-SNAP/compute_snap.cpp
index 52c1b32fbc..25c5c488fa 100644
--- a/src/ML-SNAP/compute_snap.cpp
+++ b/src/ML-SNAP/compute_snap.cpp
@@ -208,7 +208,11 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
   if (bikflag) bik_rows = natoms;
   dgrad_rows = ndims_force*natoms;
   size_array_rows = bik_rows+dgrad_rows + ndims_virial;
-  if (dgradflag) size_array_cols = nperdim + 3;
+  if (dgradflag){
+    size_array_rows = bik_rows + 3*natoms*natoms + 1;
+    size_array_cols = nperdim + 3;
+    error->warning(FLERR,"dgradflag=1 creates a N^2 array, beware of large systems.");
+  }
   else size_array_cols = nperdim*atom->ntypes + 1;
   lastcol = size_array_cols-1;
 
@@ -236,13 +240,7 @@ ComputeSnap::~ComputeSnap()
     memory->destroy(sinnerelem);
     memory->destroy(dinnerelem);
   }
-  if (dgradflag){
-    memory->destroy(dgrad);
-    memory->destroy(nneighs);
-    memory->destroy(neighsum);
-    memory->destroy(icounter);
-    memory->destroy(dbiri);
-  }
+
 }
 
 /* ---------------------------------------------------------------------- */
@@ -306,9 +304,6 @@ void ComputeSnap::init_list(int /*id*/, NeighList *ptr)
 void ComputeSnap::compute_array()
 {
 
-  if (dgradflag) get_dgrad_length();
-  //  if (dgradflag) get_dgrad_length2();
-
   int ntotal = atom->nlocal + atom->nghost;
 
   invoked_array = update->ntimestep;
@@ -448,195 +443,130 @@ void ComputeSnap::compute_array()
       for (int jj = 0; jj < ninside; jj++) {
         const int j = snaptr->inside[jj];
 
-        if (dgradflag){
-          dgrad_row_indx = 3*neighsum[atom->tag[j]-1] + 3*icounter[atom->tag[j]-1] ;
-          icounter[atom->tag[j]-1] += 1;
-        }
-
         snaptr->compute_duidrj(jj);
         snaptr->compute_dbidrj();
 
         // accumulate dBi/dRi, -dBi/dRj
 
-	if (!dgradflag) {
+	      if (!dgradflag) {
 
-	  double *snadi = snap_peratom[i]+typeoffset_local;
-	  double *snadj = snap_peratom[j]+typeoffset_local;
+	         double *snadi = snap_peratom[i]+typeoffset_local;
+	         double *snadj = snap_peratom[j]+typeoffset_local;
 
-	  for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+	         for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
 
-	    snadi[icoeff] += snaptr->dblist[icoeff][0];
-	    snadi[icoeff+yoffset] += snaptr->dblist[icoeff][1];
-	    snadi[icoeff+zoffset] += snaptr->dblist[icoeff][2];
+	           snadi[icoeff] += snaptr->dblist[icoeff][0];
+	           snadi[icoeff+yoffset] += snaptr->dblist[icoeff][1];
+	           snadi[icoeff+zoffset] += snaptr->dblist[icoeff][2];
 
-	    snadj[icoeff] -= snaptr->dblist[icoeff][0];
-	    snadj[icoeff+yoffset] -= snaptr->dblist[icoeff][1];
-	    snadj[icoeff+zoffset] -= snaptr->dblist[icoeff][2];
-	  }
+	           snadj[icoeff] -= snaptr->dblist[icoeff][0];
+	           snadj[icoeff+yoffset] -= snaptr->dblist[icoeff][1];
+	           snadj[icoeff+zoffset] -= snaptr->dblist[icoeff][2];
+	         }
 
-	  if (quadraticflag) {
-	    const int quadraticoffset = ncoeff;
-	    snadi += quadraticoffset;
-	    snadj += quadraticoffset;
-	    int ncount = 0;
-	    for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
-	      double bi = snaptr->blist[icoeff];
-	      double bix = snaptr->dblist[icoeff][0];
-	      double biy = snaptr->dblist[icoeff][1];
-	      double biz = snaptr->dblist[icoeff][2];
+	         if (quadraticflag) {
+	           const int quadraticoffset = ncoeff;
+	           snadi += quadraticoffset;
+	           snadj += quadraticoffset;
+	           int ncount = 0;
+	           for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+	             double bi = snaptr->blist[icoeff];
+	             double bix = snaptr->dblist[icoeff][0];
+	             double biy = snaptr->dblist[icoeff][1];
+	             double biz = snaptr->dblist[icoeff][2];
 
-	      // diagonal elements of quadratic matrix
+	             // diagonal elements of quadratic matrix
 
-	      double dbxtmp = bi*bix;
-	      double dbytmp = bi*biy;
-	      double dbztmp = bi*biz;
+	             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;
+	             snadi[ncount] +=         dbxtmp;
+	             snadi[ncount+yoffset] += dbytmp;
+	             snadi[ncount+zoffset] += dbztmp;
+	             snadj[ncount] -=         dbxtmp;
+	             snadj[ncount+yoffset] -= dbytmp;
+	             snadj[ncount+zoffset] -= dbztmp;
 
-	      ncount++;
+	             ncount++;
 
-	      // upper-triangular elements of quadratic matrix
+	             // upper-triangular elements of quadratic matrix
 
-	      for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
-		double dbxtmp = bi*snaptr->dblist[jcoeff][0]
-		  + bix*snaptr->blist[jcoeff];
-		double dbytmp = bi*snaptr->dblist[jcoeff][1]
-		  + biy*snaptr->blist[jcoeff];
-		double dbztmp = bi*snaptr->dblist[jcoeff][2]
-		  + biz*snaptr->blist[jcoeff];
+	             for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
+		             double dbxtmp = bi*snaptr->dblist[jcoeff][0]
+		                           + bix*snaptr->blist[jcoeff];
+		             double dbytmp = bi*snaptr->dblist[jcoeff][1]
+		                           + biy*snaptr->blist[jcoeff];
+		             double dbztmp = bi*snaptr->dblist[jcoeff][2]
+		                           + biz*snaptr->blist[jcoeff];
 
-		snadi[ncount] +=         dbxtmp;
-		snadi[ncount+yoffset] += dbytmp;
-		snadi[ncount+zoffset] += dbztmp;
-		snadj[ncount] -=         dbxtmp;
-		snadj[ncount+yoffset] -= dbytmp;
-		snadj[ncount+zoffset] -= dbztmp;
+		             snadi[ncount] +=         dbxtmp;
+		             snadi[ncount+yoffset] += dbytmp;
+		             snadi[ncount+zoffset] += dbztmp;
+		             snadj[ncount] -=         dbxtmp;
+		             snadj[ncount+yoffset] -= dbytmp;
+		             snadj[ncount+zoffset] -= dbztmp;
 
-		ncount++;
-	      }
+		             ncount++;
+	             }
+	           }
+	         }
+	       } else {
 
-	    }
+	         for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
 
-	  }
+             // add to snap array for this proc
 
-	} else {
+             // dBi/dRj
 
-	  for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+             snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 0][icoeff+3] -= snaptr->dblist[icoeff][0];
+             snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 1][icoeff+3] -= snaptr->dblist[icoeff][1];
+             snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 2][icoeff+3] -= snaptr->dblist[icoeff][2];
 
-	    // sign convention same as compute snad
-      /*
-	    dgrad[dgrad_row_indx+0][icoeff] = -snaptr->dblist[icoeff][0];
-	    dgrad[dgrad_row_indx+1][icoeff] = -snaptr->dblist[icoeff][1];
-	    dgrad[dgrad_row_indx+2][icoeff] = -snaptr->dblist[icoeff][2];
+             // dBi/dRi
 
-	    // accumulate dBi/dRi = sum (-dBi/dRj) for neighbors j of if i
+             snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 0][icoeff+3] += snaptr->dblist[icoeff][0];
+             snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 1][icoeff+3] += snaptr->dblist[icoeff][1];
+             snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 2][icoeff+3] += snaptr->dblist[icoeff][2];
+	         }
 
-	    dbiri[3*(atom->tag[i]-1)+0][icoeff] += snaptr->dblist[icoeff][0];
-	    dbiri[3*(atom->tag[i]-1)+1][icoeff] += snaptr->dblist[icoeff][1];
-	    dbiri[3*(atom->tag[i]-1)+2][icoeff] += snaptr->dblist[icoeff][2];
-      */
+	       }
 
-      // add to snap array for this proc
-      // dBi/dRj
-      snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 0][icoeff+3] -= snaptr->dblist[icoeff][0];
-      snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 1][icoeff+3] -= snaptr->dblist[icoeff][1];
-      snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 2][icoeff+3] -= snaptr->dblist[icoeff][2];
-      // dBi/dRi
-      snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 0][icoeff+3] += snaptr->dblist[icoeff][0];
-      snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 1][icoeff+3] += snaptr->dblist[icoeff][1];
-      snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 2][icoeff+3] += snaptr->dblist[icoeff][2];
-	  }
-
-    /*
-	  dgrad[dgrad_row_indx+0][ncoeff] = atom->tag[i]-1;
-	  dgrad[dgrad_row_indx+0][ncoeff+1] = atom->tag[j]-1;
-	  dgrad[dgrad_row_indx+0][ncoeff+2] = 0;
-	  dgrad[dgrad_row_indx+1][ncoeff] = atom->tag[i]-1;
-	  dgrad[dgrad_row_indx+1][ncoeff+1] = atom->tag[j]-1;
-	  dgrad[dgrad_row_indx+1][ncoeff+2] = 1;
-	  dgrad[dgrad_row_indx+2][ncoeff] = atom->tag[i]-1;
-	  dgrad[dgrad_row_indx+2][ncoeff+1] = atom->tag[j]-1;
-	  dgrad[dgrad_row_indx+2][ncoeff+2] = 2;
-
-	  dbiri[3*(atom->tag[i]-1)+0][ncoeff] = atom->tag[i]-1;
-	  dbiri[3*(atom->tag[i]-1)+0][ncoeff+1] = atom->tag[i]-1;
-	  dbiri[3*(atom->tag[i]-1)+0][ncoeff+2] = 0;
-
-	  dbiri[3*(atom->tag[i]-1)+1][ncoeff] = atom->tag[i]-1;
-	  dbiri[3*(atom->tag[i]-1)+1][ncoeff+1] = atom->tag[i]-1;
-	  dbiri[3*(atom->tag[i]-1)+1][ncoeff+2] = 1;
-
-	  dbiri[3*(atom->tag[i]-1)+2][ncoeff] = atom->tag[i]-1;
-	  dbiri[3*(atom->tag[i]-1)+2][ncoeff+1] = atom->tag[i]-1;
-	  dbiri[3*(atom->tag[i]-1)+2][ncoeff+2] = 2;
-    */
-
-    // add to snap array for this proc
-    // dBi/dRj tags
-    /*
-    snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 0][0] = atom->tag[i]-1;
-    snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 0][1] = atom->tag[j]-1;
-    snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 0][2] = 0;
-    snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 1][0] = atom->tag[i]-1;
-    snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 1][1] = atom->tag[j]-1;
-    snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 1][2] = 1;
-    snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 2][0] = atom->tag[i]-1;
-    snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 2][1] = atom->tag[j]-1;
-    snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 2][2] = 2;
-    */
-
-    // dBi/dRi tags
-    /*
-    snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 0][0] = atom->tag[i]-1;
-    snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 0][1] = atom->tag[i]-1;
-    snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 0][2] = 0;
-    snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 1][0] = atom->tag[i]-1;
-    snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 1][1] = atom->tag[i]-1;
-    snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 1][2] = 1;
-    snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 2][0] = atom->tag[i]-1;
-    snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 2][1] = atom->tag[i]-1;
-    snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 2][2] = 2;
-    */
-
-	}
-} // loop over jj inside
+      } // loop over jj inside
 
       // accumulate Bi
 
       if (!dgradflag) {
 
-	// linear contributions
+	      // linear contributions
 
         int k = typeoffset_global;
         for (int icoeff = 0; icoeff < ncoeff; icoeff++)
           snap[irow][k++] += snaptr->blist[icoeff];
 
-	// quadratic contributions
+	      // quadratic contributions
 
-	if (quadraticflag) {
-	  for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
-	    double bveci = snaptr->blist[icoeff];
-	    snap[irow][k++] += 0.5*bveci*bveci;
-	    for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
-	      double bvecj = snaptr->blist[jcoeff];
-	      snap[irow][k++] += bveci*bvecj;
-	    }
-	  }
-	}
+	      if (quadraticflag) {
+  	      for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+  	        double bveci = snaptr->blist[icoeff];
+  	        snap[irow][k++] += 0.5*bveci*bveci;
+  	        for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
+  	          double bvecj = snaptr->blist[jcoeff];
+  	          snap[irow][k++] += bveci*bvecj;
+  	        }
+  	      }
+	      }
 
       } else {
-        int k = 3;
-        for (int icoeff = 0; icoeff < ncoeff; icoeff++)
-          snap[irow][k++] += snaptr->blist[icoeff];
-	numneigh_sum += ninside;
-      }
-    }
-  }
+          int k = 3;
+          for (int icoeff = 0; icoeff < ncoeff; icoeff++)
+            snap[irow][k++] += snaptr->blist[icoeff];
+  	      numneigh_sum += ninside;
+        }
+
+    } // if (mask[i] & groupbit)
+
+  } // for (int ii = 0; ii < inum; ii++) {
 
   // accumulate bispectrum force contributions to global array
 
@@ -657,8 +587,6 @@ void ComputeSnap::compute_array()
       }
     }
 
-  } else {
-
   }
 
  // accumulate forces to global array
@@ -702,9 +630,8 @@ void ComputeSnap::compute_array()
 
   } else {
 
-    // Assign reference energy right after the dgrad rows, first column.
-    // Add 3N for the dBi/dRi rows.
-    //int irow = bik_rows + dgrad_rows + 3*natoms;
+    // assign reference energy right after the dgrad rows, first column
+
     int irow = bik_rows + 3*natoms*natoms;
     double reference_energy = c_pe->compute_scalar();
     snapall[irow][0] = reference_energy;
@@ -766,194 +693,6 @@ void ComputeSnap::dbdotr_compute()
     }
 }
 
-/* ----------------------------------------------------------------------
-   compute dgrad length
-------------------------------------------------------------------------- */
-
-void ComputeSnap::get_dgrad_length()
-{
-
-  rank = universe->me; // for MPI debugging
-
-  memory->destroy(snap);
-  memory->destroy(snapall);
-
-  // invoke full neighbor list
-
-  neighbor->build_one(list);
-  dgrad_rows = 0;
-  const int inum = list->inum;
-  const int* const ilist = list->ilist;
-  const int* const numneigh = list->numneigh;
-  int** const firstneigh = list->firstneigh;
-  int * const type = atom->type;
-  const int* const mask = atom->mask;
-  double** const x = atom->x;
-
-  memory->create(neighsum, natoms, "snap:neighsum");
-  memory->create(nneighs, natoms, "snap:nneighs");
-  memory->create(icounter, natoms, "snap:icounter");
-  memory->create(dbiri, 3*natoms,ncoeff+3, "snap:dbiri");
-  //if (atom->nlocal != natoms)
-  //  error->all(FLERR,"Compute snap dgradflag=1 does not support parallelism yet.");
-
-  for (int ii = 0; ii < 3 * natoms; ii++)
-    for (int icoeff = 0; icoeff < ncoeff; icoeff++)
-      dbiri[ii][icoeff] = 0.0;
-
-  for (int ii = 0; ii < inum; ii++) {
-    const int i = ilist[ii];
-    if (mask[i] & groupbit) {
-      icounter[i] = 0;
-      nneighs[i] = 0;
-      const double xtmp = x[i][0];
-      const double ytmp = x[i][1];
-      const double ztmp = x[i][2];
-      const int itype = type[i];
-      const int* const jlist = firstneigh[i];
-      const int jnum = numneigh[i];
-      for (int jj = 0; jj < jnum; jj++) {
-        int j = jlist[jj];
-        j &= NEIGHMASK;
-
-        const double delx = x[j][0] - xtmp;
-        const double dely = x[j][1] - ytmp;
-        const double delz = x[j][2] - ztmp;
-        const double rsq = delx * delx + dely * dely + delz * delz;
-        int jtype = type[j];
-
-        if (rsq < cutsq[itype][jtype] && rsq>1e-20) {
-          dgrad_rows++;
-          nneighs[i]++;
-        }
-      }
-    }
-  }
-
-  dgrad_rows *= ndims_force;
-
-  neighsum[0] = 0;
-  for (int ii = 1; ii < inum; ii++) {
-    const int i = ilist[ii];
-    if (mask[i] & groupbit)
-      neighsum[i] = neighsum[i-1] + nneighs[i-1];
-  }
-
-  memory->create(dgrad, dgrad_rows, ncoeff+3, "snap:dgrad");
-  for (int i = 0; i < dgrad_rows; i++)
-    for (int j = 0; j < ncoeff+3; j++)
-      dgrad[i][j] = 0.0;
-
-  // set size array rows which now depends on dgrad_rows.
-
-  //size_array_rows = bik_rows + dgrad_rows + 3*atom->nlocal + 1; // Add 3*N for dBi/dRi. and add 1 for reference energy
-  size_array_rows = bik_rows + 3*natoms*natoms + 1;
-  //printf("----- dgrad_rows, 3*natoms*natoms: %d %d\n", dgrad_rows, 3*natoms*natoms);
-
-  memory->create(snap,size_array_rows,size_array_cols, "snap:snap");
-  memory->create(snapall,size_array_rows,size_array_cols, "snap:snapall");
-  array = snapall;
-
-}
-
-/* ----------------------------------------------------------------------
-   compute dgrad length
-------------------------------------------------------------------------- */
-
-void ComputeSnap::get_dgrad_length2()
-{
-  memory->destroy(snap);
-  memory->destroy(snapall);
-
-  // invoke full neighbor list
-
-  neighbor->build_one(list);
-  dgrad_rows = 0;
-  const int inum = list->inum;
-  const int* const ilist = list->ilist;
-  const int* const numneigh = list->numneigh;
-  int** const firstneigh = list->firstneigh;
-  int * const type = atom->type;
-  const int* const mask = atom->mask;
-  double** const x = atom->x;
-
-  memory->create(neighsum, natoms, "snap:neighsum");
-  memory->create(nneighs, natoms, "snap:nneighs");
-  memory->create(icounter, natoms, "snap:icounter");
-  memory->create(dbiri, 3*natoms,ncoeff+3, "snap:dbiri");
-  if (atom->nlocal != natoms)
-    error->all(FLERR,"Compute snap dgradflag=1 does not support parallelism yet.");
-
-  for (int ii = 0; ii < 3 * natoms; ii++)
-    for (int icoeff = 0; icoeff < ncoeff; icoeff++)
-      dbiri[ii][icoeff] = 0.0;
-
-  for (int ii = 0; ii < inum; ii++) {
-    const int i = ilist[ii];
-    if (mask[i] & groupbit) {
-      icounter[i] = 0;
-      nneighs[i] = 0;
-      const double xtmp = x[i][0];
-      const double ytmp = x[i][1];
-      const double ztmp = x[i][2];
-      const int itype = type[i];
-      const int* const jlist = firstneigh[i];
-      const int jnum = numneigh[i];
-      for (int jj = 0; jj < jnum; jj++) {
-        int j = jlist[jj];
-        j &= NEIGHMASK;
-
-        const double delx = x[j][0] - xtmp;
-        const double dely = x[j][1] - ytmp;
-        const double delz = x[j][2] - ztmp;
-        const double rsq = delx * delx + dely * dely + delz * delz;
-        int jtype = type[j];
-
-        if (rsq < cutsq[itype][jtype] && rsq>1e-20) {
-          dgrad_rows++;
-          nneighs[i]++;
-        }
-      }
-    }
-  }
-
-  dgrad_rows *= ndims_force;
-
-  // loop over all atoms again to calculate neighsum
-
-  // for (int ii = 0; ii < inum; ii++) {
-  //   const int i = ilist[ii];
-  //   if (mask[i] & groupbit) {
-  //     for (int jj = 0; jj < ii; jj++) {
-  // 	const int j = ilist[jj];
-  // 	if (mask[j] & groupbit)
-  // 	  neighsum[i] += nneighs[j];
-  //     }
-  //   }
-  // }
-
-  neighsum[0] = 0;
-  for (int ii = 1; ii < inum; ii++) {
-    const int i = ilist[ii];
-    if (mask[i] & groupbit)
-      neighsum[i] = neighsum[i-1] + nneighs[i-1];
-  }
-
-  memory->create(dgrad, dgrad_rows, ncoeff+3, "snap:dgrad");
-  for (int i = 0; i < dgrad_rows; i++)
-    for (int j = 0; j < ncoeff+3; j++)
-      dgrad[i][j] = 0.0;
-
-  // set size array rows which now depends on dgrad_rows.
-
-  size_array_rows = bik_rows + dgrad_rows + 3*atom->nlocal + 1; // Add 3*N for dBi/dRi. and add 1 for reference energy
-
-  memory->create(snap,size_array_rows,size_array_cols, "snap:snap");
-  memory->create(snapall,size_array_rows,size_array_cols, "snap:snapall");
-  array = snapall;
-
-}
-
 /* ----------------------------------------------------------------------
    memory usage
 ------------------------------------------------------------------------- */
diff --git a/src/ML-SNAP/compute_snap.h b/src/ML-SNAP/compute_snap.h
index e1fb7cee49..a89a870fe8 100644
--- a/src/ML-SNAP/compute_snap.h
+++ b/src/ML-SNAP/compute_snap.h
@@ -53,20 +53,13 @@ class ComputeSnap : public Compute {
   double cutmax;
   int quadraticflag;
   int bikflag, bik_rows, dgradflag, dgrad_rows;
-  double **dgrad; // First ncoeff columns are descriptor derivatives.
-                  // Last 3 columns are indices i,j,a
-  double **dbiri; // dBi/dRi = sum(-dBi/dRj) over neighbors j
-  int *nneighs; // number of neighs inside the snap cutoff.
-  int *neighsum;
-  int *icounter; // counting atoms i for each j.
   int rank;
 
   Compute *c_pe;
   Compute *c_virial;
 
   void dbdotr_compute();
-  void get_dgrad_length();
-  void get_dgrad_length2();
+
 };
 
 }    // namespace LAMMPS_NS