diff --git a/doc/src/compute_sna_atom.rst b/doc/src/compute_sna_atom.rst
index 54a6df02a2..e3255f201d 100644
--- a/doc/src/compute_sna_atom.rst
+++ b/doc/src/compute_sna_atom.rst
@@ -33,7 +33,7 @@ Syntax
 * R_1, R_2,... = list of cutoff radii, one for each type (distance units)
 * w_1, w_2,... = list of neighbor weights, one for each type
 * zero or more keyword/value pairs may be appended
-* keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag* or *bikflag* or *switchinnerflag* or *sinner* or *dinner*
+* keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag* or *bikflag* or *switchinnerflag* or *sinner* or *dinner* or *dgradflag*
 
   .. parsed-literal::
 
@@ -66,6 +66,9 @@ Syntax
           *sinnerlist* = *ntypes* values of *Sinner* (distance units)
        *dinner* values = *dinnerlist*
           *dinnerlist* = *ntypes* values of *Dinner* (distance units)
+       *dgradflag* value = *0* or *1*
+          *0* = bispectrum descriptor gradients are summed over neighbors
+          *1* = bispectrum descriptor gradients are not summed over neighbors
 
 Examples
 """"""""
@@ -340,6 +343,14 @@ When the central atom and the neighbor atom have different types,
 the values of :math:`S_{inner}` and :math:`D_{inner}` are
 the arithmetic means of the values for both types.
 
+The keyword *dgradflag* determines whether or not to sum the bispectrum descriptor gradients over neighboring atoms *i'*
+as explained with *snad/atom* above. If *dgradflag* is set to 1 then the descriptor gradient rows of the global snap array
+are not summed over atoms *i'*. Instead, each row corresponds to a single term :math:`\frac{\partial {B_{i,k}  }}{\partial {r}^a_j}`
+where :math:`a` is the Cartesian direction for the gradient. This also changes
+the number of columns to be equal to the number of bispectrum components, with 3
+additional columns representing the indices :math:`i`, :math:`j`, and :math:`a`,
+as explained more in the Output info section below. The option *dgradflag=1* must be used with *bikflag=1*.
+
 .. note::
 
    If you have a bonded system, then the settings of :doc:`special_bonds
@@ -435,6 +446,36 @@ components. For the purposes of handling contributions to force, virial,
 and quadratic combinations, these :math:`N_{elem}^3` sub-blocks are
 treated as a single block of :math:`K N_{elem}^3` columns.
 
+If the *bik* keyword is set to 1, then the first :math:`N` rows of the snap array
+correspond to :math:`B_{i,k}` instead of the sum over atoms :math:`i`. In this case, the entries in the final column for these rows
+are set to zero.
+
+If the *dgradflag* keyword is set to 1, this changes the structure of the snap array completely.
+Here the *snad/atom* quantities are replaced with rows corresponding to descriptor
+gradient components
+
+.. math::
+
+  \frac{\partial {B_{i,k}  }}{\partial {r}^a_j}
+
+where :math:`a` is the Cartesian direction for the gradient. The rows are organized in chunks, where each chunk corresponds to
+an atom :math:`j` in the system of :math:`N` atoms. The rows in an atom :math:`j` chunk correspond to neighbors :math:`i` of :math:`j`.
+The number of rows in the atom :math:`j` chunk is therefore equal to the number of neighbors :math:`N_{neighs}[j]` within the SNAP
+potential cutoff radius of atom :math:`j`, times 3 for each Cartesian direction.
+The total number of rows for these descriptor gradients is therefore
+
+.. math::
+
+  3 \sum_j^{N} N_{neighs}[j].
+
+For *dgradflag=1*, the number of columns is equal to the number of bispectrum components,
+plus 3 additional columns representing the indices :math:`i`, :math:`j`, and :math:`a` which
+identify the atoms :math:`i` and :math:`j`, and Cartesian direction :math:`a` for which
+a particular gradient :math:`\frac{\partial {B_{i,k}  }}{\partial {r}^a_j}` belongs to. The reference energy and forces are also located in different parts of the array.
+The last 3 columns of the first :math:`N` rows belong to the reference potential force components.
+The first column of the last row, after the first :math:`N + 3 \sum_j^{N} N_{neighs}[j]` rows,
+contains the reference potential energy. The virial components are not used with this option.
+
 These values can be accessed by any command that uses per-atom values
 from a compute as input.  See the :doc:`Howto output <Howto_output>` doc
 page for an overview of LAMMPS output options. To see how this command
diff --git a/examples/snap/README.md b/examples/snap/README.md
index 7c34fe7021..e1bd54ff37 100644
--- a/examples/snap/README.md
+++ b/examples/snap/README.md
@@ -1,9 +1 @@
-See `compute_snap_dgrad.py` for a test that compares the dBi/dRj from compute snap (`dbirjflag=1`) to the sum of dBi/dRj from usual compute snap (`dbirjflag=0`).
-
-The format of the global array from `dbirjflag=1` is as follows.
-
-The first N rows belong to bispectrum components for each atom, if we use `bikflag=1`. The first `K` columns correspond to each bispectrum coefficient. The final 3 columns contain reference force components for each atom.
-
-The rows after the first N rows contain dBi/dRj values for all pairs. These values are arranged in row chunks for each atom `j`, where all the rows in a chunk are associated with the neighbors `i` of `j`, as well as the self-terms where `i=j`. So for atom `j`, the number of rows is equal to the number of atoms within the SNAP cutoff, plus 1 for the `i=j` terms, times 3 for each Cartesian component. The total number of dBi/dRj rows is therefore equal to `N*(Nneigh+1)*3`, and `Nneigh` may be different for each atom. To facilitate with determining which row belong to which atom pair `ij`, the last 3 columns contain indices; the 3rd to last column contains global indices of atoms `i` (the neighbors), the 2nd to last column contains global indices of atoms `j`, and the last column contains an index 0,1,2 for the Cartesian component. Like the `bik` rows, the first `K` columns correspond to each bispectrum coefficient.
-
-Finally, the first column of the last row contains the reference energy.
+See `compute_snap_dgrad.py` for a test that compares the dBi/dRj from compute snap (`dgradflag=1`) to the sum of dBi/dRj from usual compute snap (`dgradflag=0`).
diff --git a/examples/snap/compute_snap_dgrad.py b/examples/snap/compute_snap_dgrad.py
index dc06ae1576..84d0b8113a 100644
--- a/examples/snap/compute_snap_dgrad.py
+++ b/examples/snap/compute_snap_dgrad.py
@@ -1,8 +1,8 @@
 """
 compute_snap_dgrad.py
 Purpose: Demonstrate extraction of descriptor gradient (dB/dR) array from compute snap.
-         Show that dBi/dRj components summed over neighbors i yields same output as regular compute snap with dbirjflag=0.
-         This shows that the dBi/dRj components extracted with dbirjflag=1 are correct.
+         Show that dBi/dRj components summed over neighbors i yields same output as regular compute snap with dgradflag=0.
+         This shows that the dBi/dRj components extracted with dgradflag=1 are correct.
 Serial syntax:
     python compute_snap_dgrad.py
 Parallel syntax:
@@ -24,7 +24,7 @@ from lammps import lammps, LMP_TYPE_ARRAY, LMP_STYLE_GLOBAL
 cmds = ["-screen", "none", "-log", "none"]
 lmp = lammps(cmdargs=cmds)
 
-def run_lammps(dbirjflag):
+def run_lammps(dgradflag):
     lmp.command("clear")
     lmp.command("units metal")
     lmp.command("boundary	p p p")
@@ -36,7 +36,7 @@ def run_lammps(dbirjflag):
     lmp.command("mass 		* 180.88")
     lmp.command("displace_atoms 	all random 0.01 0.01 0.01 123456")
     # Pair style
-    snap_options=f'{rcutfac} {rfac0} {twojmax} {radelem1} {radelem2} {wj1} {wj2} rmin0 {rmin0} quadraticflag {quadratic} bzeroflag {bzero} switchflag {switch} bikflag {bikflag} dbirjflag {dbirjflag}'
+    snap_options=f'{rcutfac} {rfac0} {twojmax} {radelem1} {radelem2} {wj1} {wj2} rmin0 {rmin0} quadraticflag {quadratic} bzeroflag {bzero} switchflag {switch} bikflag {bikflag} dgradflag {dgradflag}'
     lmp.command(f"pair_style 	zero {rcutfac}")
     lmp.command(f"pair_coeff 	* *")
     lmp.command(f"pair_style 	zbl {zblcutinner} {zblcutouter}")
@@ -70,7 +70,7 @@ quadratic=0
 bzero=0
 switch=0
 bikflag=1
-dbirjflag=1
+dgradflag=1
 
 # Declare reference potential variables
 zblcutinner=4.0
@@ -84,7 +84,8 @@ else:
     nd = int(m*(m+1)*(m+2)/3)
 print(f"Number of descriptors based on twojmax : {nd}")
 
-# Run lammps with dbirjflag on
+# Run lammps with dgradflag on
+print("Running with dgradflag on")
 run_lammps(1)
 
 # Get global snap array
@@ -93,12 +94,12 @@ lmp_snap = lmp.numpy.extract_compute("snap",0, 2)
 # Extract dBj/dRi (includes dBi/dRi)
 natoms = lmp.get_natoms()
 fref1 = lmp_snap[0:natoms,-3:].flatten()
-eref1 = lmp_snap[-6,0]
-dbdr_length = np.shape(lmp_snap)[0]-(natoms)-6 # Length of neighborlist pruned dbdr array
+eref1 = lmp_snap[-1,0] #lmp_snap[-6,0]
+dbdr_length = np.shape(lmp_snap)[0]-(natoms) - 1 #-6 # Length of neighborlist pruned dbdr array
 dBdR = lmp_snap[natoms:(natoms+dbdr_length),:]
 force_indices = lmp_snap[natoms:(natoms+dbdr_length),-3:].astype(np.int32)
 
-# Sum over neighbors j for each atom i, like dbirjflag=0 does.
+# Sum over neighbors j for each atom i, like dgradflag=0 does.
 array1 = np.zeros((3*natoms,nd))
 a = 0
 for k in range(0,nd):
@@ -111,7 +112,8 @@ for k in range(0,nd):
         if (a>2):
             a=0
 
-# Run lammps with dbirjflag off
+# Run lammps with dgradflag off
+print("Running with dgradflag off")
 run_lammps(0)
 
 # Get global snap array
@@ -121,7 +123,7 @@ fref2 = lmp_snap[natoms:(natoms+3*natoms),-1]
 eref2 = lmp_snap[0,-1]
 array2 = lmp_snap[natoms:natoms+(3*natoms), nd:-1]
 
-# Sum the arrays obtained from dbirjflag on and off.
+# Sum the arrays obtained from dgradflag on and off.
 summ = array1 + array2
 #np.savetxt("sum.dat", summ)
 
diff --git a/src/ML-SNAP/compute_snap.cpp b/src/ML-SNAP/compute_snap.cpp
index d1b0980dfb..88360a5005 100644
--- a/src/ML-SNAP/compute_snap.cpp
+++ b/src/ML-SNAP/compute_snap.cpp
@@ -56,7 +56,7 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
   bzeroflag = 1;
   quadraticflag = 0;
   bikflag = 0;
-  dbirjflag = 0;
+  dgradflag = 0;
   chemflag = 0;
   bnormflag = 0;
   wselfallflag = 0;
@@ -148,10 +148,10 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
         error->all(FLERR,"Illegal compute snap command");
       bikflag = atoi(arg[iarg+1]);
       iarg += 2;
-    } else if (strcmp(arg[iarg],"dbirjflag") == 0) {
+    } else if (strcmp(arg[iarg],"dgradflag") == 0) {
       if (iarg+2 > narg)
         error->all(FLERR,"Illegal compute snap command");
-      dbirjflag = atoi(arg[iarg+1]);
+      dgradflag = atoi(arg[iarg+1]);
       iarg += 2;
     } else if (strcmp(arg[iarg],"switchinnerflag") == 0) {
       if (iarg+2 > narg)
@@ -185,6 +185,9 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
   if (!switchinnerflag && (sinnerflag || dinnerflag))
     error->all(FLERR,"Illegal compute snap command: switchinnerflag = 0, unexpected sinner/dinner keyword");
 
+  if (dgradflag && !bikflag)
+    error->all(FLERR,"Illegal compute snap command: dgradflag=1 requires bikflag=1");
+
   snaptr = new SNA(lmp, rfac0, twojmax,
                    rmin0, switchflag, bzeroflag,
                    chemflag, bnormflag, wselfallflag,
@@ -200,9 +203,9 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
   natoms = atom->natoms;
   bik_rows = 1;
   if (bikflag) bik_rows = natoms;
-  dbirj_rows = ndims_force*natoms;
-  size_array_rows = bik_rows+dbirj_rows+ndims_virial;
-  if (dbirjflag) size_array_cols = nperdim+3; // plus 3 for tag[i], tag[j], and cartesian index
+  dgrad_rows = ndims_force*natoms;
+  size_array_rows = bik_rows+dgrad_rows+ndims_virial;
+  if (dgradflag) size_array_cols = nperdim+3; // plus 3 for tag[i], tag[j], and cartesian index
   else size_array_cols = nperdim*atom->ntypes+1;
   lastcol = size_array_cols-1;
 
@@ -230,8 +233,8 @@ ComputeSnap::~ComputeSnap()
     memory->destroy(sinnerelem);
     memory->destroy(dinnerelem);
   }
-  if (dbirjflag){
-    memory->destroy(dbirj);
+  if (dgradflag){
+    memory->destroy(dgrad);
     memory->destroy(nneighs);
     memory->destroy(neighsum);
     memory->destroy(icounter);
@@ -260,10 +263,19 @@ void ComputeSnap::init()
 
   // allocate memory for global array
 
-  memory->create(snap,size_array_rows,size_array_cols,
-                 "snap:snap");
-  memory->create(snapall,size_array_rows,size_array_cols,
-                 "snap:snapall");
+  if (dgradflag){
+    // Initially allocate natoms^2 rows, will prune with neighborlist later
+    memory->create(snap,natoms*natoms,ncoeff+3,
+                   "snap:snap");
+    memory->create(snapall,natoms*natoms,ncoeff+3,
+                   "snap:snapall");
+  }
+  else{
+    memory->create(snap,size_array_rows,size_array_cols,
+                   "snap:snap");
+    memory->create(snapall,size_array_rows,size_array_cols,
+                   "snap:snapall");
+  }
   array = snapall;
 
   // find compute for reference energy
@@ -299,8 +311,8 @@ void ComputeSnap::init_list(int /*id*/, NeighList *ptr)
 void ComputeSnap::compute_array()
 {
 
-  if (dbirjflag){
-    get_dbirj_length();
+  if (dgradflag){
+    get_dgrad_length();
   }
 
   int ntotal = atom->nlocal + atom->nghost;
@@ -345,7 +357,7 @@ void ComputeSnap::compute_array()
   const int* const mask = atom->mask;
   int ninside;
   int numneigh_sum = 0;
-  int dbirj_row_indx;
+  int dgrad_row_indx;
   for (int ii = 0; ii < inum; ii++) {
     int irow = 0;
     if (bikflag) irow = atom->tag[ilist[ii] & NEIGHMASK]-1;
@@ -416,8 +428,8 @@ void ComputeSnap::compute_array()
       for (int jj = 0; jj < ninside; jj++) {
         const int j = snaptr->inside[jj];
 
-        if (dbirjflag){
-          dbirj_row_indx = 3*neighsum[atom->tag[j]-1] + 3*icounter[atom->tag[j]-1] ;
+        if (dgradflag){
+          dgrad_row_indx = 3*neighsum[atom->tag[j]-1] + 3*icounter[atom->tag[j]-1] ;
           icounter[atom->tag[j]-1] += 1;
         }
 
@@ -440,20 +452,20 @@ void ComputeSnap::compute_array()
           snadj[icoeff+zoffset] -= snaptr->dblist[icoeff][2];
 
 
-          if (dbirjflag){
-            dbirj[dbirj_row_indx+0][icoeff] = snaptr->dblist[icoeff][0];
-            dbirj[dbirj_row_indx+1][icoeff] = snaptr->dblist[icoeff][1];
-            dbirj[dbirj_row_indx+2][icoeff] = snaptr->dblist[icoeff][2];
+          if (dgradflag){
+            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];
             if (icoeff==(ncoeff-1)){
-              dbirj[dbirj_row_indx+0][ncoeff] = atom->tag[i]-1;
-              dbirj[dbirj_row_indx+0][ncoeff+1] = atom->tag[j]-1;
-              dbirj[dbirj_row_indx+0][ncoeff+2] = 0;
-              dbirj[dbirj_row_indx+1][ncoeff] = atom->tag[i]-1;
-              dbirj[dbirj_row_indx+1][ncoeff+1] = atom->tag[j]-1;
-              dbirj[dbirj_row_indx+1][ncoeff+2] = 1;
-              dbirj[dbirj_row_indx+2][ncoeff] = atom->tag[i]-1;
-              dbirj[dbirj_row_indx+2][ncoeff+1] = atom->tag[j]-1;
-              dbirj[dbirj_row_indx+2][ncoeff+2] = 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;
             }
             // Accumulate dBi/dRi = sum (-dBi/dRj) for neighbors j of if i
             dbiri[3*(atom->tag[i]-1)+0][icoeff] -= snaptr->dblist[icoeff][0];
@@ -530,7 +542,7 @@ void ComputeSnap::compute_array()
 
       // linear contributions
       int k;
-      if (dbirjflag) k = 0;
+      if (dgradflag) k = 0;
       else k = typeoffset_global;
       for (int icoeff = 0; icoeff < ncoeff; icoeff++){
         snap[irow][k++] += snaptr->blist[icoeff];
@@ -555,7 +567,7 @@ void ComputeSnap::compute_array()
 
   // Accumulate contributions to global array
 
-  if (dbirjflag){
+  if (dgradflag){
 
     int dbiri_indx;
     int irow;
@@ -568,34 +580,34 @@ void ComputeSnap::compute_array()
 
             for (int jj=0; jj<nneighs[i]; jj++){
 
-              int dbirj_row_indx = 3*neighsum[atom->tag[i]-1] + 3*jj;
+              int dgrad_row_indx = 3*neighsum[atom->tag[i]-1] + 3*jj;
               int snap_row_indx = 3*neighsum[atom->tag[i]-1] + 3*(atom->tag[i]-1) + 3*jj;
               irow = snap_row_indx + bik_rows;
 
               // x-coordinate
-              snap[irow][icoeff+typeoffset_global] += dbirj[dbirj_row_indx+0][icoeff];
+              snap[irow][icoeff+typeoffset_global] += dgrad[dgrad_row_indx+0][icoeff];
               if (icoeff==(ncoeff-1)){
-                snap[irow][ncoeff] += dbirj[dbirj_row_indx+0][ncoeff];
-                snap[irow][ncoeff+1] += dbirj[dbirj_row_indx+0][ncoeff+1];
-                snap[irow][ncoeff+2] += dbirj[dbirj_row_indx+0][ncoeff+2];
+                snap[irow][ncoeff] += dgrad[dgrad_row_indx+0][ncoeff];
+                snap[irow][ncoeff+1] += dgrad[dgrad_row_indx+0][ncoeff+1];
+                snap[irow][ncoeff+2] += dgrad[dgrad_row_indx+0][ncoeff+2];
               }
               irow++;
 
               // y-coordinate
-              snap[irow][icoeff+typeoffset_global] += dbirj[dbirj_row_indx+1][icoeff];
+              snap[irow][icoeff+typeoffset_global] += dgrad[dgrad_row_indx+1][icoeff];
               if (icoeff==(ncoeff-1)){
-                snap[irow][ncoeff] += dbirj[dbirj_row_indx+1][ncoeff];
-                snap[irow][ncoeff+1] += dbirj[dbirj_row_indx+1][ncoeff+1];
-                snap[irow][ncoeff+2] += dbirj[dbirj_row_indx+1][ncoeff+2];
+                snap[irow][ncoeff] += dgrad[dgrad_row_indx+1][ncoeff];
+                snap[irow][ncoeff+1] += dgrad[dgrad_row_indx+1][ncoeff+1];
+                snap[irow][ncoeff+2] += dgrad[dgrad_row_indx+1][ncoeff+2];
               }
               irow++;
 
               // z-coordinate
-              snap[irow][icoeff+typeoffset_global]   += dbirj[dbirj_row_indx+2][icoeff];
+              snap[irow][icoeff+typeoffset_global]   += dgrad[dgrad_row_indx+2][icoeff];
               if (icoeff==(ncoeff-1)){
-                snap[irow][ncoeff] += dbirj[dbirj_row_indx+2][ncoeff];
-                snap[irow][ncoeff+1] += dbirj[dbirj_row_indx+2][ncoeff+1];
-                snap[irow][ncoeff+2] += dbirj[dbirj_row_indx+2][ncoeff+2];
+                snap[irow][ncoeff] += dgrad[dgrad_row_indx+2][ncoeff];
+                snap[irow][ncoeff+1] += dgrad[dgrad_row_indx+2][ncoeff+1];
+                snap[irow][ncoeff+2] += dgrad[dgrad_row_indx+2][ncoeff+2];
               }
               dbiri_indx = dbiri_indx+3;
             }
@@ -659,8 +671,8 @@ void ComputeSnap::compute_array()
 
  // accumulate forces to global array
 
- if (dbirjflag){
-   // for dbirjflag=1, put forces at last 3 columns of bik rows
+ if (dgradflag){
+   // for dgradflag=1, put forces at last 3 columns of bik rows
    for (int i=0; i<atom->nlocal; i++){
      int iglobal = atom->tag[i];
      snap[iglobal-1][ncoeff+0] = atom->f[i][0];
@@ -681,8 +693,8 @@ void ComputeSnap::compute_array()
 
   // accumulate bispectrum virial contributions to global array
 
-  if (dbirjflag){
-    // no virial terms for dbirj yet
+  if (dgradflag){
+    // no virial terms for dgrad yet
   }
   else{
     dbdotr_compute();
@@ -691,10 +703,10 @@ void ComputeSnap::compute_array()
   // sum up over all processes
   MPI_Allreduce(&snap[0][0],&snapall[0][0],size_array_rows*size_array_cols,MPI_DOUBLE,MPI_SUM,world);
   // assign energy to last column
-  if (dbirjflag){
-    // Assign reference energy right after the dbirj rows, first column.
+  if (dgradflag){
+    // Assign reference energy right after the dgrad rows, first column.
     // Add 3N for the dBi/dRi rows.
-    int irow = bik_rows + dbirj_rows + 3*natoms;
+    int irow = bik_rows + dgrad_rows + 3*natoms;
     double reference_energy = c_pe->compute_scalar();
     snapall[irow][0] = reference_energy;
   }
@@ -709,10 +721,11 @@ void ComputeSnap::compute_array()
   // switch to Voigt notation
 
   c_virial->compute_vector();
-  if (dbirjflag){
-    // no virial terms for dbirj yet
+  if (dgradflag){
+    // no virial terms for dgrad yet
   }
   else{
+    c_virial->compute_vector();
     int irow = 3*natoms+bik_rows;
     snapall[irow++][lastcol] = c_virial->vector[0];
     snapall[irow++][lastcol] = c_virial->vector[1];
@@ -760,10 +773,10 @@ void ComputeSnap::dbdotr_compute()
 }
 
 /* ----------------------------------------------------------------------
-   compute dbirj length
+   compute dgrad length
 ------------------------------------------------------------------------- */
 
-void ComputeSnap::get_dbirj_length()
+void ComputeSnap::get_dgrad_length()
 {
 
   int rank = universe->me; // for MPI debugging
@@ -772,7 +785,7 @@ void ComputeSnap::get_dbirj_length()
   memory->destroy(snapall);
   // invoke full neighbor list
   neighbor->build_one(list);
-  dbirj_rows = 0;
+  dgrad_rows = 0;
   const int inum = list->inum;
   const int* const ilist = list->ilist;
   const int* const numneigh = list->numneigh;
@@ -791,7 +804,7 @@ void ComputeSnap::get_dbirj_length()
   memory->create(icounter, natoms, "snap:icounter");
   memory->create(dbiri, 3*natoms,ncoeff+3, "snap:dbiri");
   if (atom->nlocal != natoms){
-    error->all(FLERR,"Compute snap dbirjflag=1 does not support parallelism.");
+    error->all(FLERR,"Compute snap dgradflag=1 does not support parallelism.");
   }
   for (int ii=0; ii<3*natoms; ii++){
     for (int icoeff=0; icoeff<ncoeff; icoeff++){
@@ -822,7 +835,7 @@ void ComputeSnap::get_dbirj_length()
         int jtype = type[j];
 
         if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
-          dbirj_rows += 1; //jnum + 1;
+          dgrad_rows += 1; //jnum + 1;
           jnum_cutoff += 1;
           nneighs[i]+=1;
         }
@@ -830,7 +843,7 @@ void ComputeSnap::get_dbirj_length()
     }
   }
 
-  dbirj_rows *= ndims_force;
+  dgrad_rows *= ndims_force;
 
   // Loop over all atoms again to calculate neighsum.
   for (int ii = 0; ii < inum; ii++) {
@@ -850,14 +863,15 @@ void ComputeSnap::get_dbirj_length()
     }
   }
 
-  memory->create(dbirj, dbirj_rows, ncoeff+3, "snap:dbirj");
-  for (int i=0; i<dbirj_rows; i++){
+  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++){
-      dbirj[i][j]=0.0;
+      dgrad[i][j]=0.0;
     }
   }
-  // Set size array rows which now depends on dbirj_rows.
-  size_array_rows = bik_rows+dbirj_rows+ndims_virial+3*atom->nlocal; // Add 3*N for dBi/dRi
+  // Set size array rows which now depends on dgrad_rows.
+  //size_array_rows = bik_rows+dgrad_rows+ndims_virial+3*atom->nlocal; // Add 3*N for dBi/dRi
+  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");
@@ -871,7 +885,7 @@ void ComputeSnap::get_dbirj_length()
 
 double ComputeSnap::compute_scalar()
 {
-  if (dbirjflag) get_dbirj_length();
+  if (dgradflag) get_dgrad_length();
   return size_array_rows;
 }
 
diff --git a/src/ML-SNAP/compute_snap.h b/src/ML-SNAP/compute_snap.h
index 3b03512efc..447f6d286d 100644
--- a/src/ML-SNAP/compute_snap.h
+++ b/src/ML-SNAP/compute_snap.h
@@ -56,8 +56,8 @@ class ComputeSnap : public Compute {
   int quadraticflag;
   //int bikflag;
   //int bik_rows;
-  int bikflag, bik_rows, dbirjflag, dbirj_rows;
-  double **dbirj;
+  int bikflag, bik_rows, dgradflag, dgrad_rows;
+  double **dgrad;
   double **dbiri; // dBi/dRi = sum(-dBi/dRj) over neighbors j
   int *nneighs; // number of neighs inside the snap cutoff.
   int *neighsum;
@@ -67,7 +67,7 @@ class ComputeSnap : public Compute {
   Compute *c_virial;
 
   void dbdotr_compute();
-  void get_dbirj_length();
+  void get_dgrad_length();
 };
 
 }    // namespace LAMMPS_NS
diff --git a/src/ML-SNAP/compute_snapneigh.cpp b/src/ML-SNAP/compute_snapneigh.cpp
deleted file mode 100644
index eef18a93f4..0000000000
--- a/src/ML-SNAP/compute_snapneigh.cpp
+++ /dev/null
@@ -1,563 +0,0 @@
-// clang-format off
-/* ----------------------------------------------------------------------
-   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
-   https://www.lammps.org/, Sandia National Laboratories
-   Steve Plimpton, sjplimp@sandia.gov
-
-   Copyright (2003) Sandia Corporation.  Under the terms of Contract
-   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
-   certain rights in this software.  This software is distributed under
-   the GNU General Public License.
-
-   See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-#include "compute_snapneigh.h"
-
-#include "sna.h"
-#include "atom.h"
-#include "update.h"
-#include "modify.h"
-#include "neighbor.h"
-#include "neigh_list.h"
-#include "force.h"
-#include "pair.h"
-#include "comm.h"
-#include "memory.h"
-#include "error.h"
-
-#include <cstring>
-
-using namespace LAMMPS_NS;
-
-enum{SCALAR,VECTOR,ARRAY};
-
-ComputeSnapneigh::ComputeSnapneigh(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), cutsq(nullptr), list(nullptr), radelem(nullptr), wjelem(nullptr),
-  sinnerelem(nullptr), dinnerelem(nullptr)
-{
-
-  array_flag = 1;
-  //vector_flag = 1;
-  extarray = 0;
-
-  double rfac0, rmin0;
-  int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag;
-
-  int ntypes = atom->ntypes;
-  int nargmin = 6+2*ntypes;
-
-  if (narg < nargmin) error->all(FLERR,"Illegal compute snap command");
-
-  // default values
-
-  rmin0 = 0.0;
-  switchflag = 1;
-  bzeroflag = 1;
-  quadraticflag = 0;
-  bikflag = 0;
-  dbirjflag = 0;
-  chemflag = 0;
-  bnormflag = 0;
-  wselfallflag = 0;
-  switchinnerflag = 0;
-  nelements = 1;
-
-  // process required arguments
-
-  memory->create(radelem,ntypes+1,"snapneigh:radelem"); // offset by 1 to match up with types
-  memory->create(wjelem,ntypes+1,"snapneigh:wjelem");
-  rcutfac = atof(arg[3]);
-  rfac0 = atof(arg[4]);
-  twojmax = atoi(arg[5]);
-  for (int i = 0; i < ntypes; i++)
-    radelem[i+1] = atof(arg[6+i]);
-  for (int i = 0; i < ntypes; i++)
-    wjelem[i+1] = atof(arg[6+ntypes+i]);
-
-  // construct cutsq
-
-  double cut;
-  cutmax = 0.0;
-  memory->create(cutsq,ntypes+1,ntypes+1,"snapneigh:cutsq");
-  for (int i = 1; i <= ntypes; i++) {
-    cut = 2.0*radelem[i]*rcutfac;
-    //printf("cut: %f\n", cut);
-    if (cut > cutmax) cutmax = cut;
-    cutsq[i][i] = cut*cut;
-    for (int j = i+1; j <= ntypes; j++) {
-      cut = (radelem[i]+radelem[j])*rcutfac;
-      cutsq[i][j] = cutsq[j][i] = cut*cut;
-    }
-  }
-
-  // set local input checks
-
-  int sinnerflag = 0;
-  int dinnerflag = 0;
-
-  // process optional args
-
-  int iarg = nargmin;
-
-  while (iarg < narg) {
-    if (strcmp(arg[iarg],"rmin0") == 0) {
-      if (iarg+2 > narg)
-        error->all(FLERR,"Illegal compute snap command");
-      rmin0 = atof(arg[iarg+1]);
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"bzeroflag") == 0) {
-      if (iarg+2 > narg)
-        error->all(FLERR,"Illegal compute snap command");
-      bzeroflag = atoi(arg[iarg+1]);
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"switchflag") == 0) {
-      if (iarg+2 > narg)
-        error->all(FLERR,"Illegal compute snap command");
-      switchflag = atoi(arg[iarg+1]);
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"quadraticflag") == 0) {
-      if (iarg+2 > narg)
-        error->all(FLERR,"Illegal compute snap command");
-      quadraticflag = atoi(arg[iarg+1]);
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"chem") == 0) {
-      if (iarg+2+ntypes > narg)
-        error->all(FLERR,"Illegal compute snap command");
-      chemflag = 1;
-      memory->create(map,ntypes+1,"compute_snapneigh:map");
-      nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
-      for (int i = 0; i < ntypes; i++) {
-        int jelem = utils::inumeric(FLERR,arg[iarg+2+i],false,lmp);
-        if (jelem < 0 || jelem >= nelements)
-          error->all(FLERR,"Illegal compute snap command");
-        map[i+1] = jelem;
-      }
-      iarg += 2+ntypes;
-    } else if (strcmp(arg[iarg],"bnormflag") == 0) {
-      if (iarg+2 > narg)
-        error->all(FLERR,"Illegal compute snap command");
-      bnormflag = atoi(arg[iarg+1]);
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"wselfallflag") == 0) {
-      if (iarg+2 > narg)
-        error->all(FLERR,"Illegal compute snap command");
-      wselfallflag = atoi(arg[iarg+1]);
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"bikflag") == 0) {
-      if (iarg+2 > narg)
-        error->all(FLERR,"Illegal compute snap command");
-      bikflag = atoi(arg[iarg+1]);
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"dbirjflag") == 0) {
-      if (iarg+2 > narg)
-        error->all(FLERR,"Illegal compute snap command");
-      dbirjflag = atoi(arg[iarg+1]);
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"switchinnerflag") == 0) {
-      if (iarg+2 > narg)
-        error->all(FLERR,"Illegal compute snap command");
-      switchinnerflag = atoi(arg[iarg+1]);
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"sinner") == 0) {
-      iarg++;
-      if (iarg+ntypes > narg)
-        error->all(FLERR,"Illegal compute snap command");
-      memory->create(sinnerelem,ntypes+1,"snapneigh:sinnerelem");
-      for (int i = 0; i < ntypes; i++)
-        sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp);
-      sinnerflag = 1;
-      iarg += ntypes;
-    } else if (strcmp(arg[iarg],"dinner") == 0) {
-      iarg++;
-      if (iarg+ntypes > narg)
-        error->all(FLERR,"Illegal compute snap command");
-      memory->create(dinnerelem,ntypes+1,"snapneigh:dinnerelem");
-      for (int i = 0; i < ntypes; i++)
-        dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp);
-      dinnerflag = 1;
-      iarg += ntypes;
-    } else error->all(FLERR,"Illegal compute snap command");
-  }
-
-  if (switchinnerflag && !(sinnerflag && dinnerflag))
-    error->all(FLERR,"Illegal compute snap command: switchinnerflag = 1, missing sinner/dinner keyword");
-
-  if (!switchinnerflag && (sinnerflag || dinnerflag))
-    error->all(FLERR,"Illegal compute snap command: switchinnerflag = 0, unexpected sinner/dinner keyword");
-
-  /*
-  snaptr = new SNA(lmp, rfac0, twojmax,
-                   rmin0, switchflag, bzeroflag,
-                   chemflag, bnormflag, wselfallflag,
-                   nelements, switchinnerflag);
-  */
-
-  //ncoeff = snaptr->ncoeff;
-  nperdim = ncoeff;
-  if (quadraticflag) nperdim += (ncoeff*(ncoeff+1))/2;
-  ndims_force = 3;
-  ndims_virial = 6;
-  yoffset = nperdim;
-  zoffset = 2*nperdim;
-  natoms = atom->natoms;
-  bik_rows = 1;
-  if (bikflag) bik_rows = natoms;
-  //size_array_rows = bik_rows+ndims_force*natoms+ndims_virial;
-  dbirj_rows = ndims_force*natoms;
-  size_array_rows = bik_rows+dbirj_rows+ndims_virial;
-  size_array_cols = nperdim*atom->ntypes+1;
-  lastcol = size_array_cols-1;
-
-  ndims_peratom = ndims_force;
-  size_peratom = ndims_peratom*nperdim*atom->ntypes;
-
-  nmax = 0;
-
-}
-
-/* ---------------------------------------------------------------------- */
-
-ComputeSnapneigh::~ComputeSnapneigh()
-{
-
-  memory->destroy(neighs);
-  memory->destroy(radelem);
-  memory->destroy(wjelem);
-  memory->destroy(cutsq);
-  //delete snaptr;
-
-  if (chemflag) memory->destroy(map);
-
-  if (switchinnerflag) {
-    memory->destroy(sinnerelem);
-    memory->destroy(dinnerelem);
-  }
-
-  if (dbirjflag){
-    //printf("dbirj_rows: %d\n", dbirj_rows);
-    //memory->destroy(dbirj);
-    memory->destroy(nneighs);
-    memory->destroy(neighsum);
-    memory->destroy(icounter);
-    //memory->destroy(dbiri);
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void ComputeSnapneigh::init()
-{
-  if (force->pair == nullptr)
-    error->all(FLERR,"Compute snap requires a pair style be defined");
-
-  if (cutmax > force->pair->cutforce){
-    //printf("----- cutmax cutforce: %f %f\n", cutmax, force->pair->cutforce);
-    error->all(FLERR,"Compute snap cutoff is longer than pairwise cutoff");
-  }
-
-  // need an occasional full neighbor list
-
-  neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL);
-
-}
-
-
-/* ---------------------------------------------------------------------- */
-
-void ComputeSnapneigh::init_list(int /*id*/, NeighList *ptr)
-{
-  list = ptr;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void ComputeSnapneigh::compute_array()
-{
-
-  if (dbirjflag){
-    //printf("----- dbirjflag true.\n");
-    get_dbirj_length();
-    //printf("----- got dbirj_length\n");
-  }
-
-  //printf("----- cutmax cutforce: %f %f\n", cutmax, force->pair->cutforce);
-  //else{
-  int ntotal = atom->nlocal + atom->nghost;
-
-  invoked_array = update->ntimestep;
-
-  // invoke full neighbor list (will copy or build if necessary)
-  neighbor->build_one(list);
-
-  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;
-
-  // compute sna derivatives for each atom in group
-  // use full neighbor list to count atoms less than cutoff
-
-  double** const x = atom->x;
-  const int* const mask = atom->mask;
-  //printf("----- inum: %d\n", inum);
-  //printf("----- NEIGHMASK: %d\n", NEIGHMASK);
-  int ninside;
-  int numneigh_sum = 0;
-  int dbirj_row_indx;
-  int dbiri_indx=0;
-  for (int ii = 0; ii < inum; ii++) {
-    int irow = 0;
-    if (bikflag) irow = atom->tag[ilist[ii] & NEIGHMASK]-1;
-    //printf("----- i, itag: %d %d\n", ilist[ii] & NEIGHMASK, atom->tag[ilist[ii]]);
-    const int i = ilist[ii];
-    //printf("----- ii, i: %d %d\n", ii, i);
-    //printf("----- mask[i] groupbit: %d %d\n", mask[i], groupbit);
-    if (mask[i] & groupbit) {
-
-      const double xtmp = x[i][0];
-      const double ytmp = x[i][1];
-      const double ztmp = x[i][2];
-      const int itype = type[i];
-      int ielem = 0;
-      if (chemflag)
-        ielem = map[itype];
-      const double radi = radelem[itype];
-      const int* const jlist = firstneigh[i];
-      const int jnum = numneigh[i];
-      const int typeoffset_local = ndims_peratom*nperdim*(itype-1);
-      const int typeoffset_global = nperdim*(itype-1);
-
-      // insure rij, inside, and typej  are of size jnum
-
-      //snaptr->grow_rij(jnum);
-
-      // rij[][3] = displacements between atom I and those neighbors
-      // inside = indices of neighbors of I within cutoff
-      // typej = types of neighbors of I within cutoff
-      // note Rij sign convention => dU/dRij = dU/dRj = -dU/dRi
-
-      /*
-      This loop assigns quantities in snaptr.
-      snaptr is a SNA class instance, see sna.h
-
-      */
-      //int ninside = 0;
-      ninside=0;
-      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];
-        int jelem = 0;
-        if (chemflag)
-          jelem = map[jtype];
-        if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
-          if (dbirjflag){
-            //dbirj_row_indx = 3*neighsum[atom->tag[j]-1] + 3*icounter[atom->tag[j]-1] ; // THIS IS WRONG, SEE NEXT VAR.
-            //dbirj_row_indx = 3*neighsum[atom->tag[i]-1] + 3*(atom->tag[i]-1) + 3*icounter[atom->tag[j]-1]; // 3*tagi is to leave space for dBi/dRi
-            dbirj_row_indx = 3*neighsum[atom->tag[j]-1] + 3*icounter[atom->tag[j]-1] + 3*(atom->tag[j]-1); // THIS IS WRONG, SEE NEXT VAR.
-            //printf("--- %d %d %d %d\n", dbirj_row_indx, 3*neighsum[atom->tag[i]-1], 3*(atom->tag[i]-1), jj);
-            //printf("jtag, icounter, dbirj_row_indx: %d, %d, %d %d %d\n", atom->tag[j], icounter[atom->tag[j]-1], dbirj_row_indx+0, dbirj_row_indx+1, dbirj_row_indx+2);
-            icounter[atom->tag[j]-1] += 1;
-
-            neighs[dbirj_row_indx+0][0] = atom->tag[i];
-            neighs[dbirj_row_indx+1][0] = atom->tag[i];
-            neighs[dbirj_row_indx+2][0] = atom->tag[i];
-
-            neighs[dbirj_row_indx+0][1] = atom->tag[j];
-            neighs[dbirj_row_indx+1][1] = atom->tag[j];
-            neighs[dbirj_row_indx+2][1] = atom->tag[j];
-
-            neighs[dbirj_row_indx+0][2] = 0;
-            neighs[dbirj_row_indx+1][2] = 1;
-            neighs[dbirj_row_indx+2][2] = 2;
-
-            dbiri_indx = dbiri_indx+3;
-          }
-        }
-      }
-
-      //printf("--- dbiri_indx: %d\n", dbiri_indx);
-      // Put dBi/dRi in
-
-      neighs[dbiri_indx+0][0] = atom->tag[i];
-      neighs[dbiri_indx+1][0] = atom->tag[i];
-      neighs[dbiri_indx+2][0] = atom->tag[i];
-
-      neighs[dbiri_indx+0][1] = atom->tag[i];
-      neighs[dbiri_indx+1][1] = atom->tag[i];
-      neighs[dbiri_indx+2][1] = atom->tag[i];
-
-      neighs[dbiri_indx+0][2] = 0;
-      neighs[dbiri_indx+1][2] = 1;
-      neighs[dbiri_indx+2][2] = 2;
-
-      dbiri_indx = dbiri_indx+3;
-    }
-
-    numneigh_sum += ninside;
-  } // for (int ii = 0; ii < inum; ii++)
-
-
-  // Check icounter.
-  /*
-  for (int ii = 0; ii < inum; ii++) {
-    const int i = ilist[ii];
-    if (mask[i] & groupbit) {
-      printf("icounter[i]: %d\n", icounter[i]);
-    }
-  }
-  */
-
-
-  // sum up over all processes
-  // I'll need to do something like this...
-  /*
-  MPI_Allreduce(&snap[0][0],&snapall[0][0],size_array_rows*size_array_cols,MPI_DOUBLE,MPI_SUM,world);
-  // assign energy to last column
-  for (int i = 0; i < bik_rows; i++) snapall[i][lastcol] = 0;
-  int irow = 0;
-  double reference_energy = c_pe->compute_scalar();
-  snapall[irow][lastcol] = reference_energy;
-  */
-
-  /*
-  for (int i=0; i<dbirj_rows; i++){
-    printf("----- %d: %f %f\n", i, array[i][0], array[i][1]);
-  }
-
-  //printf("vector[0]: %d\n", vector[0]);
-  printf("----- End of compute snapneigh.\n");
-  */
-}
-
-/* ----------------------------------------------------------------------
-   compute dbirj length
-------------------------------------------------------------------------- */
-
-void ComputeSnapneigh::get_dbirj_length()
-{
-  // invoke full neighbor list (will copy or build if necessary)
-  neighbor->build_one(list);
-  dbirj_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;
-  //printf("----- inum: %d\n", inum);
-  memory->create(neighsum, atom->nlocal, "snapneigh:neighsum");
-  memory->create(nneighs, atom->nlocal, "snapneigh:nneighs");
-  memory->create(icounter, atom->nlocal, "snapneigh:icounter");
-  //memory->create(dbiri, 3*inum,ncoeff, "snapneigh:dbiri");
-  /*
-  for (int ii=0; ii<inum; 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;
-      neighsum[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];
-      //printf("-----   jnum: %d\n", jnum);
-      int jnum_cutoff = 0;
-      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) {
-          dbirj_rows += 1; //jnum + 1;
-          jnum_cutoff += 1;
-          nneighs[i]+=1;
-        }
-      }
-      //printf("----- jnum_cutoff: %d\n", jnum_cutoff);
-    }
-  }
-
-  dbirj_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) {
-      //printf("nneighs[i]: %d\n", nneighs[i]);
-      //neighsum[i] = 0;
-      //printf("i nneighs[i]: %d %d\n", i, nneighs[i]);
-      if (i==0){
-        neighsum[i]=0;
-      }
-      else{
-        for (int jj=0; jj < ii; jj++){
-          const int j = ilist[jj];
-          if (mask[j] & groupbit) {
-            //printf("  j nneighs[j-1]: %d %d\n", j, nneighs[j]);
-            neighsum[i] += nneighs[j];
-          }
-        }
-        //neighsum[i] += 1; // Add 1 for the self term dBi/dRi
-      }
-    }
-    //printf("%d\n", neighsum[i]);
-  }
-
-  size_array_rows = dbirj_rows+(3*atom->nlocal);
-  size_array_cols = 3;
-  //memory->create(dbirj, dbirj_rows, ncoeff, "snapneigh:dbirj");
-  memory->create(neighs, size_array_rows, size_array_cols, "snapneigh:neighs");
-
-  //vector = neighs;
-  array = neighs;
-  // Set size array rows which now depends on dbirj_rows.
-  //size_array_rows = bik_rows+dbirj_rows+ndims_virial;
-  //printf("----- dbirj_rows: %d\n", dbirj_rows);
-  //printf("----- end of dbirj length.\n");
-
-}
-
-/* ----------------------------------------------------------------------
-   compute array length
-------------------------------------------------------------------------- */
-
-double ComputeSnapneigh::compute_scalar()
-{
-  if (dbirjflag) get_dbirj_length();
-  return size_array_rows;
-}
-
-/* ----------------------------------------------------------------------
-   memory usage
-------------------------------------------------------------------------- */
-
-double ComputeSnapneigh::memory_usage()
-{
-
-  double bytes = (double)size_array_rows*size_array_cols *
-    sizeof(double);                                     // array
-
-  return bytes;
-}
diff --git a/src/ML-SNAP/compute_snapneigh.h b/src/ML-SNAP/compute_snapneigh.h
deleted file mode 100644
index 1ac712a101..0000000000
--- a/src/ML-SNAP/compute_snapneigh.h
+++ /dev/null
@@ -1,73 +0,0 @@
-/* -*- c++ -*- ----------------------------------------------------------
-   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
-   https://www.lammps.org/, Sandia National Laboratories
-   Steve Plimpton, sjplimp@sandia.gov
-
-   Copyright (2003) Sandia Corporation.  Under the terms of Contract
-   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
-   certain rights in this software.  This software is distributed under
-   the GNU General Public License.
-
-   See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
-
-#ifdef COMPUTE_CLASS
-// clang-format off
-ComputeStyle(snapneigh,ComputeSnapneigh);
-// clang-format on
-#else
-
-#ifndef LMP_COMPUTE_SNAPNEIGH_H
-#define LMP_COMPUTE_SNAPNEIGH_H
-
-#include "compute.h"
-
-namespace LAMMPS_NS {
-
-class ComputeSnapneigh : public Compute {
- public:
-  ComputeSnapneigh(class LAMMPS *, int, char **);
-  ~ComputeSnapneigh() override;
-  void init() override;
-  void init_list(int, class NeighList *) override;
-  void compute_array() override;
-  double compute_scalar() override;
-  double memory_usage() override;
-
- private:
-  FILE * fh_d;
-  int natoms, nmax, size_peratom, lastcol;
-  int ncoeff, nperdim, yoffset, zoffset;
-  int ndims_peratom, ndims_force, ndims_virial;
-  double **cutsq;
-  class NeighList *list;
-  double **snap, **snapall;
-  double **snap_peratom;
-  double rcutfac;
-  double *radelem;
-  double *wjelem;
-  int *map;    // map types to [0,nelements)
-  int nelements, chemflag;
-  int switchinnerflag;
-  double *sinnerelem;
-  double *dinnerelem;
-  //class SNA *snaptr;
-  double cutmax;
-  int quadraticflag;
-  //int bikflag;
-  //int bik_rows;
-  int bikflag, bik_rows, dbirjflag, dbirj_rows;
-  double **dbirj;
-  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.
-  double **neighs; // neighborlist for neural networks
-
-  void get_dbirj_length();
-};
-
-}    // namespace LAMMPS_NS
-
-#endif
-#endif