Added compute snap descriptor gradient example.

examples/snap/README.md

@@ -0,0 +1,9 @@
+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.

examples/snap/compute_snap_dgrad.py

@@ -0,0 +1,130 @@
+"""
+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.
+Serial syntax:
+    python compute_snap_dgrad.py
+Parallel syntax:
+    mpirun -np 2 python compute_snap_dgrad.py
+"""
+
+from __future__ import print_function
+import sys
+import ctypes
+import numpy as np
+
+# uncomment this if running in parallel via mpi4py
+#me = 0
+#from mpi4py import MPI
+#me = MPI.COMM_WORLD.Get_rank()
+#nprocs = MPI.COMM_WORLD.Get_size()
+
+from lammps import lammps, LMP_TYPE_ARRAY, LMP_STYLE_GLOBAL
+cmds = ["-screen", "none", "-log", "none"]
+lmp = lammps(cmdargs=cmds)
+
+def run_lammps(dbirjflag):
+    lmp.command("clear")
+    lmp.command("units metal")
+    lmp.command("boundary	p p p")
+    lmp.command("atom_modify	map hash")
+    lmp.command(f"lattice         bcc {latparam}")
+    lmp.command(f"region		box block 0 {nx} 0 {ny} 0 {nz}")
+    lmp.command(f"create_box	{ntypes} box")
+    lmp.command(f"create_atoms	{ntypes} box")
+    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}'
+    lmp.command(f"pair_style 	zero {rcutfac}")
+    lmp.command(f"pair_coeff 	* *")
+    lmp.command(f"pair_style 	zbl {zblcutinner} {zblcutouter}")
+    lmp.command(f"pair_coeff 	* * {zblz} {zblz}")
+    # set up compute snap generating global array
+    lmp.command(f"compute 	snap all snap {snap_options}")
+    # Run
+    lmp.command(f"thermo 		100")
+    lmp.command(f"run {nsteps}")
+
+# Declare simulation/structure variables
+nsteps=0
+nrep=2
+latparam=2.0
+ntypes=2
+nx=nrep
+ny=nrep
+nz=nrep
+
+# Declare compute snap variables
+twojmax=8
+m = (twojmax/2)+1
+rcutfac=1.0
+rfac0=0.99363
+rmin0=0
+radelem1=2.3
+radelem2=2.0
+wj1=1.0
+wj2=0.96
+quadratic=0
+bzero=0
+switch=0
+bikflag=1
+dbirjflag=1
+
+# Declare reference potential variables
+zblcutinner=4.0
+zblcutouter=4.8
+zblz=73
+
+# Number of descriptors
+if (twojmax % 2 == 0):
+    nd = int(m*(m+1)*(2*m+1)/6)
+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(1)
+
+# Get global snap array
+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
+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.
+array1 = np.zeros((3*natoms,nd))
+a = 0
+for k in range(0,nd):
+    for l in range(0,dbdr_length):
+        j = force_indices[l,0]
+        i = force_indices[l,1]
+        #array1[3*(i-1)+a,k] += dBdR[l,k]
+        array1[3*(i)+a,k] += dBdR[l,k]
+        a = a+1
+        if (a>2):
+            a=0
+
+# Run lammps with dbirjflag off
+run_lammps(0)
+
+# Get global snap array
+lmp_snap = lmp.numpy.extract_compute("snap",0, 2)
+natoms = lmp.get_natoms()
+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.
+summ = array1 + array2
+#np.savetxt("sum.dat", summ)
+
+print(f"Maximum difference in descriptor sums: {np.max(summ)}")
+print(f"Maximum difference in reference forces: {np.max(fref1-fref2)}")
+print(f"Difference in reference energy: {np.max(eref1-eref2)}")

python/examples/in.lj

@@ -0,0 +1,25 @@
+# 3d Lennard-Jones melt
+
+units		lj
+atom_style	atomic
+atom_modify	map array
+
+lattice		fcc 0.8442
+region		box block 0 4 0 4 0 4
+create_box	1 box
+create_atoms	1 box
+mass		1 1.0
+
+velocity	all create 1.44 87287 loop geom
+
+pair_style	lj/cut 2.5
+pair_coeff	1 1 1.0 1.0 2.5
+
+neighbor	0.3 bin
+neigh_modify	delay 0 every 20 check no
+
+fix		1 all nve
+
+variable        fx atom fx
+
+run		10

python/examples/simple.py

@@ -5,10 +5,10 @@
 # Purpose: mimic operation of examples/COUPLE/simple/simple.cpp via Python
 
 # Serial syntax: simple.py in.lammps
-#                in.lammps = LAMMPS input script
+#                in.simple = LAMMPS input script
 
-# Parallel syntax: mpirun -np 4 simple.py in.lammps
+# Parallel syntax: mpirun -np 4 python simple.py in.simple
-#                  in.lammps = LAMMPS input script
+#                  in.simple = LAMMPS input script
 # also need to uncomment mpi4py sections below
 
 from __future__ import print_function
@@ -27,9 +27,9 @@ infile = sys.argv[1]
 me = 0
 
 # uncomment this if running in parallel via mpi4py
-#from mpi4py import MPI
+from mpi4py import MPI
-#me = MPI.COMM_WORLD.Get_rank()
+me = MPI.COMM_WORLD.Get_rank()
-#nprocs = MPI.COMM_WORLD.Get_size()
+nprocs = MPI.COMM_WORLD.Get_size()
 
 from lammps import lammps
 lmp = lammps()
@@ -122,10 +122,10 @@ if me == 0: print("Gather post scatter subset:",
 boxlo,boxhi,xy,yz,xz,periodicity,box_change = lmp.extract_box()
 if me == 0: print("Box info",boxlo,boxhi,xy,yz,xz,periodicity,box_change)
 
-lmp.reset_box([0,0,0],[10,10,8],0,0,0)
+#lmp.reset_box([0,0,0],[10,10,8],0,0,0)
 
-boxlo,boxhi,xy,yz,xz,periodicity,box_change = lmp.extract_box()
+#boxlo,boxhi,xy,yz,xz,periodicity,box_change = lmp.extract_box()
-if me == 0: print("Box info",boxlo,boxhi,xy,yz,xz,periodicity,box_change)
+#if me == 0: print("Box info",boxlo,boxhi,xy,yz,xz,periodicity,box_change)
 
 # uncomment if running in parallel via mpi4py
-#print("Proc %d out of %d procs has" % (me,nprocs), lmp)
+print("Proc %d out of %d procs has" % (me,nprocs), lmp)

src/ML-SNAP/compute_snap.cpp

@@ -13,7 +13,6 @@
 ------------------------------------------------------------------------- */
 
 #include "compute_snap.h"
-
 #include "sna.h"
 #include "atom.h"
 #include "update.h"
@@ -25,6 +24,7 @@
 #include "comm.h"
 #include "memory.h"
 #include "error.h"
+#include "universe.h" // For MPI
 
 #include <cstring>
 
@@ -82,7 +82,6 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
   memory->create(cutsq,ntypes+1,ntypes+1,"snap: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++) {
@@ -201,10 +200,9 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
   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;
-  if (dbirjflag) size_array_cols = nperdim+3; // plus 3 for tag[i], tag[2], and cartesian index
+  if (dbirjflag) 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;
 
@@ -220,41 +218,25 @@ ComputeSnap::~ComputeSnap()
 {
 
   memory->destroy(snap);
-  //printf("---------------- 1\n");
   memory->destroy(snapall);
-  //printf("---------------- 2\n");
   memory->destroy(snap_peratom);
-  //printf("---------------- 3\n");
   memory->destroy(radelem);
-  //printf("---------------- 4\n");
   memory->destroy(wjelem);
-  //printf("---------------- 5\n");
   memory->destroy(cutsq);
-  //printf("---------------- 6\n");
   delete snaptr;
-  //printf("---------------- 7\n");
   if (chemflag) memory->destroy(map);
 
   if (switchinnerflag) {
     memory->destroy(sinnerelem);
     memory->destroy(dinnerelem);
   }
-  //printf("---------------- 8\n");
   if (dbirjflag){
-    //printf("dbirj_rows: %d\n", dbirj_rows);
-    //printf("----- destroy dbirj\n");
     memory->destroy(dbirj);
-    //printf("----- 1-1-1-1-1-1\n");
     memory->destroy(nneighs);
-    //printf("----- 2-1-1-1-1-1\n");
     memory->destroy(neighsum);
-    //printf("----- 3-1-1-1-1-1\n");
     memory->destroy(icounter);
-    //printf("----- 4-1-1-1-1-1\n");
     memory->destroy(dbiri);
-    //printf("----- 5-1-1-1-1-1\n");
   }
-  //printf("---------------- 9\n");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -265,7 +247,6 @@ void ComputeSnap::init()
     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");
   }
 
@@ -279,7 +260,6 @@ void ComputeSnap::init()
 
   // allocate memory for global array
 
-  //printf("----- dbirjflag: %d\n", dbirjflag);
   memory->create(snap,size_array_rows,size_array_cols,
                  "snap:snap");
   memory->create(snapall,size_array_rows,size_array_cols,
@@ -319,18 +299,10 @@ void ComputeSnap::init_list(int /*id*/, NeighList *ptr)
 void ComputeSnap::compute_array()
 {
 
-  //printf(" -----2 dbirjflag: %d\n", dbirjflag);
   if (dbirjflag){
-    //printf("----- dbirjflag true.\n");
     get_dbirj_length();
-    //printf("----- got dbirj_length\n");
   }
-  //else{
-  //  printf("----- dbirjflag false.\n");
-  //}
 
-  //printf("----- cutmax cutforce: %f %f\n", cutmax, force->pair->cutforce);
-  //else{
   int ntotal = atom->nlocal + atom->nghost;
 
   invoked_array = update->ntimestep;
@@ -343,11 +315,10 @@ void ComputeSnap::compute_array()
     memory->create(snap_peratom,nmax,size_peratom,
                    "snap:snap_peratom");
   }
+
   // clear global array
-  //printf("size_array_rows: %d\n", size_array_rows);
   for (int irow = 0; irow < size_array_rows; irow++){
     for (int icoeff = 0; icoeff < size_array_cols; icoeff++){
-      //printf("%d %d\n", irow, icoeff);
       snap[irow][icoeff] = 0.0;
     }
   }
@@ -372,18 +343,13 @@ void ComputeSnap::compute_array()
 
   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;
   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];
@@ -408,12 +374,8 @@ void ComputeSnap::compute_array()
       // typej = types of neighbors of I within cutoff
       // note Rij sign convention => dU/dRij = dU/dRj = -dU/dRi
 
-      /*
+      // assign quantities in snaptr
-      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];
@@ -428,7 +390,6 @@ void ComputeSnap::compute_array()
         if (chemflag)
           jelem = map[jtype];
         if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
-          //printf("cutsq: %f\n", cutsq[itype][jtype]);
           snaptr->rij[ninside][0] = delx;
           snaptr->rij[ninside][1] = dely;
           snaptr->rij[ninside][2] = delz;
@@ -444,88 +405,36 @@ void ComputeSnap::compute_array()
         }
       }
 
-      /*
+      // compute bispectrum for atom i
-      Now that we have assigned neighbor quantities with previous loop, we are ready to compute things.
+
-      Here we compute the wigner functions (U), Z is some other quantity, and bi is bispectrum.
-      */
       snaptr->compute_ui(ninside, ielem);
       snaptr->compute_zi();
       snaptr->compute_bi(ielem);
 
-      /*
+      // loop over neighbors for descriptors derivatives
-      Looks like this loop computes derivatives.
+
-      How does snaptr know what atom I we're dealing with?
-      I think it only needs neighbor info, and then it goes from there.
-      */
-      //printf("----- Derivative loop - looping over neighbors j.\n");
-      //printf("-----   ninside: %d\n", ninside); // numneighs of I within cutoff
       for (int jj = 0; jj < ninside; jj++) {
-        //printf("-----   jj: %d\n", jj);
         const int j = snaptr->inside[jj];
-        //printf("-----   jj, j, jtag: %d %d %d\n", jj, j, atom->tag[j]);
+
-        //int dbirj_row_indx = 3*neighsum[i] + 3*jj ; // THIS IS WRONG, SEE NEXT LINE.
-        //int dbirj_row_indx = 3*neighsum[j] + 3*i_indx; // need to get i_indx.
-        // How to get i_indx?
-        /*
-        i_indx must start at zero and end at (nneighs[j]-1).
-        We can start a counter for each atom j.
-        Maybe this icounter can serve as an index for i as a neighbor of j.
-        icounter starts at zero and ends at (nneighs[j]-1).
-        */
-        //icounter[atom->tag[j]-1] += 1;
         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[j]-1] + 3*icounter[atom->tag[j]-1] ;
-          //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;
         }
-        //int dbirj_row_indx = 3*neighsum[atom->tag[j]-1] + 3*icounter[atom->tag[j]-1] ; // THIS IS WRONG, SEE NEXT VAR.
+
-        //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;
-        /*
-        j is an atom index starting from 0.
-        Use atom->tag[j] to get the atom in the box (index starts at 1).
-        Need to make sure that the order of these ij pairs is the same when multiplying by dE/dD later.
-        */
-        //printf("-----   jj, j, jtag: %d %d %d\n", jj, j, atom->tag[j]);
         snaptr->compute_duidrj(jj);
         snaptr->compute_dbidrj();
 
         // Accumulate dBi/dRi, -dBi/dRj
 
-        /*
-        snap_peratom[i] has type double * because each atom index has indices for descriptors.
-        */
         double *snadi = snap_peratom[i]+typeoffset_local;
         double *snadj = snap_peratom[j]+typeoffset_local;
-        //printf("-----   ncoeff: %d\n", ncoeff);
-        //printf("snadi: %f %f %f %f %f\n", snadi[0], snadi[1], snadi[2], snadi[3], snadi[4]);
-        //printf("----- typeoffset_local: %d\n", typeoffset_local);
-        //printf("snadi: ");
-        //for (int s=0; s<(ncoeff*3); s++){
-        //  printf("%f ", snadi[s]);
-        //}
-        /*
-        printf("snadj: ");
-        for (int s=0; s<(ncoeff*3); s++){
-          printf("%f ", snadj[s]);
-        }
-        */
-        //printf("\n");
-        for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
-          //printf("-----     dblist[icoeff]: %f %f %f\n", snaptr->dblist[icoeff][0], snaptr->dblist[icoeff][1], snaptr->dblist[icoeff][2]);
-          /*
-          I think these are the descriptor derivatives.
-          Desriptor derivatives wrt atom i.
-          What exactly is being summed here?
-          This is a loop over descriptors or coeff k.
 
-          */
+        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];
-          /*
+
-          Descriptor derivatives wrt atom j
-          */
           snadj[icoeff] -= snaptr->dblist[icoeff][0];
           snadj[icoeff+yoffset] -= snaptr->dblist[icoeff][1];
           snadj[icoeff+zoffset] -= snaptr->dblist[icoeff][2];
@@ -546,11 +455,11 @@ void ComputeSnap::compute_array()
               dbirj[dbirj_row_indx+2][ncoeff+1] = atom->tag[j]-1;
               dbirj[dbirj_row_indx+2][ncoeff+2] = 2;
             }
-            // Accumulate dBi/dRi = sum (-dBi/dRj) for neighbors j of if i.
+            // Accumulate dBi/dRi = sum (-dBi/dRj) for neighbors j of if i
             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];
-            // Get last columns
+            // Get last columns which are i, j, and Cartesian index
             if (icoeff==(ncoeff-1)){
               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;
@@ -617,18 +526,13 @@ void ComputeSnap::compute_array()
           }
 
         }
-      } // for (int jj = 0; jj < ninside; jj++)
+      }
-      //printf("---- irow after jj loop: %d\n", irow);
-
-      // Accumulate Bi
-      //printf("----- Accumulate Bi.\n");
 
       // linear contributions
       int k;
       if (dbirjflag) k = 0;
       else k = typeoffset_global;
       for (int icoeff = 0; icoeff < ncoeff; icoeff++){
-        //printf("----- %d %f\n", icoeff, snaptr->blist[icoeff]);
         snap[irow][k++] += snaptr->blist[icoeff];
       }
 
@@ -649,85 +553,25 @@ void ComputeSnap::compute_array()
     numneigh_sum += ninside;
   } // for (int ii = 0; ii < inum; ii++)
 
-  //printf("----- bik_rows: %d\n", bik_rows);
+  // Accumulate contributions to global array
-  //printf("----- bikflag: %d\n", bikflag);
-
-  // 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 all the derivatives we calculated to check usual compute snap output.
-  /*
-  if (dbirjflag){
-    fh_d = fopen("DEBUG", "w");
-    int row_indx=0;
-    double sum;
-    for (int ii=0; ii<inum; ii++){
-      printf("ii: %d\n", ii);
-      for (int a=0; a<3; a++){
-        for (int k=0; k<ncoeff; k++){
-          //double sumx = 0.0;
-          //double sumy = 0.0;
-          //double sumz = 0.0;
-          sum=0.0;
-          for (int jj=0; jj<nneighs[ii]; jj++){
-            row_indx=3*neighsum[ii]+3*jj;
-            //sumx+=dbirj[row_indx+a][k];
-            //sumy+=dbirj[row_indx+1][k];
-            //sumz+=dbirj[row_indx+2][k];
-            sum+=dbirj[row_indx+a][k];
-          }
-          // Add dBi/dRi
-          //sumx+=dbiri[3*ii+0][k];
-          //sumy+=dbiri[3*ii+1][k];
-          //sumz+=dbiri[3*ii+2][k];
-
-          sum+=dbiri[3*ii+a][k];
-
-          //printf("%d %f %f %f\n", ii, sumx, sumy, sumz);
-          fprintf(fh_d, "%f ", sum);
-        }
-        fprintf(fh_d, "\n");
-
-      }
-    }
-    fclose(fh_d);
-  }
-  */
 
   if (dbirjflag){
 
-    //printf(" Accumulate to global array.\n");
-    //printf("ntypes: %d\n", atom->ntypes);
-    //for (int itype = 0; itype < atom->ntypes; itype++) {
     int dbiri_indx;
     int irow;
     for (int itype=0; itype<1; itype++){
       const int typeoffset_local = ndims_peratom*nperdim*itype;
       const int typeoffset_global = nperdim*itype;
-      //printf("----- nperdim: %d\n", nperdim);
       for (int icoeff = 0; icoeff < nperdim; icoeff++) {
-        //printf("-----   icoeff: %d\n", icoeff);
         dbiri_indx=0;
         for (int i = 0; i < atom->nlocal; i++) {
-            //printf("i: %d\n", i);
 
-            //int dbiri_indx = 0;
-            //int irow;
             for (int jj=0; jj<nneighs[i]; jj++){
-              //printf("  jj: %d\n", jj);
+
               int dbirj_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;
-              //printf("snap_row_indx: %d\n", snap_row_indx);
-              //int irow = dbirj_row_indx+bik_rows;
               irow = snap_row_indx + bik_rows;
-              //printf("  row_indx, irow: %d %d\n", dbirj_row_indx, irow);
+
               // x-coordinate
               snap[irow][icoeff+typeoffset_global] += dbirj[dbirj_row_indx+0][icoeff];
               if (icoeff==(ncoeff-1)){
@@ -736,7 +580,7 @@ void ComputeSnap::compute_array()
                 snap[irow][ncoeff+2] += dbirj[dbirj_row_indx+0][ncoeff+2];
               }
               irow++;
-              //printf("    irow: %d\n", irow);
+
               // y-coordinate
               snap[irow][icoeff+typeoffset_global] += dbirj[dbirj_row_indx+1][icoeff];
               if (icoeff==(ncoeff-1)){
@@ -745,7 +589,7 @@ void ComputeSnap::compute_array()
                 snap[irow][ncoeff+2] += dbirj[dbirj_row_indx+1][ncoeff+2];
               }
               irow++;
-              //printf("    irow: %d\n", irow);
+
               // z-coordinate
               snap[irow][icoeff+typeoffset_global]   += dbirj[dbirj_row_indx+2][icoeff];
               if (icoeff==(ncoeff-1)){
@@ -755,11 +599,11 @@ void ComputeSnap::compute_array()
               }
               dbiri_indx = dbiri_indx+3;
             }
+
             // Put dBi/dRi at end of each dBj/dRi chunk.
-            //int dbiri_row_indx;
+
             irow = dbiri_indx + bik_rows;
-            //printf("dbiri_indx: %d\n", dbiri_indx);
+
-            //printf("dbiri: %f %f %f\n", dbiri[3*i+0][icoeff], dbiri[3*i+1][icoeff], dbiri[3*i+2][icoeff]);
             // x-coordinate
             snap[irow][icoeff+typeoffset_global] += dbiri[3*i+0][icoeff];
             if (icoeff==(ncoeff-1)){
@@ -768,7 +612,7 @@ void ComputeSnap::compute_array()
               snap[irow][ncoeff+2] += dbiri[3*i+0][ncoeff+2];
             }
             irow++;
-            //printf("    irow: %d\n", irow);
+
             // y-coordinate
             snap[irow][icoeff+typeoffset_global] += dbiri[3*i+1][icoeff];
             if (icoeff==(ncoeff-1)){
@@ -777,7 +621,8 @@ void ComputeSnap::compute_array()
               snap[irow][ncoeff+2] += dbiri[3*i+1][ncoeff+2];
             }
             irow++;
-            //printf("    irow: %d\n", irow);
+
+            // z-coordinate
             snap[irow][icoeff+typeoffset_global]   += dbiri[3*i+2][icoeff];
             if (icoeff==(ncoeff-1)){
               snap[irow][ncoeff] += dbiri[3*i+2][ncoeff];
@@ -789,31 +634,21 @@ void ComputeSnap::compute_array()
         }
       }
     }
-    //printf(" Accumulated to global array.\n");
 
   }
 
   else{
-    //printf("----- Accumulate bispecturm force contributions to global array.\n");
+
     // accumulate bispectrum force contributions to global array
-    //printf("----- ntotal, nmax, natoms: %d %d %d\n", ntotal, nmax, atom->natoms);
+
     for (int itype = 0; itype < atom->ntypes; itype++) {
       const int typeoffset_local = ndims_peratom*nperdim*itype;
       const int typeoffset_global = nperdim*itype;
-      //printf("----- nperdim: %d\n", nperdim);
-      /*nperdim = ncoeff set previsouly*/
       for (int icoeff = 0; icoeff < nperdim; icoeff++) {
-        //printf("-----   icoeff: %d\n", icoeff);
         for (int i = 0; i < ntotal; i++) {
           double *snadi = snap_peratom[i]+typeoffset_local;
           int iglobal = atom->tag[i];
-          if (icoeff==4){
-            if ( (snadi[icoeff] != 0.0) || (snadi[icoeff+yoffset] != 0.0) || (snadi[icoeff+zoffset] != 0.0) ){
-              //printf("%d %d %f %f %f\n", i, iglobal, snadi[icoeff], snadi[icoeff+yoffset], snadi[icoeff+zoffset]);
-            }
-          }
           int irow = 3*(iglobal-1)+bik_rows;
-          //printf("-----     snadi[icoeff]: %f\n", snadi[icoeff]);
           snap[irow++][icoeff+typeoffset_global] += snadi[icoeff];
           snap[irow++][icoeff+typeoffset_global] += snadi[icoeff+yoffset];
           snap[irow][icoeff+typeoffset_global] += snadi[icoeff+zoffset];
@@ -822,15 +657,11 @@ void ComputeSnap::compute_array()
     }
   }
 
-  //printf("----- Accumulate forces to global array.\n");
-  /*
-  These are the last columns.
-  */
  // accumulate forces to global array
+
  if (dbirjflag){
+   // for dbirjflag=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];
      snap[iglobal-1][ncoeff+1] = atom->f[i][1];
@@ -842,11 +673,8 @@ void ComputeSnap::compute_array()
     for (int i = 0; i < atom->nlocal; i++) {
       int iglobal = atom->tag[i];
       int irow = 3*(iglobal-1)+bik_rows;
-      //printf("---- irow: %d\n", irow);
       snap[irow++][lastcol] = atom->f[i][0];
-      //printf("---- irow: %d\n", irow);
       snap[irow++][lastcol] = atom->f[i][1];
-      //printf("---- irow: %d\n", irow);
       snap[irow][lastcol] = atom->f[i][2];
     }
   }
@@ -854,7 +682,7 @@ void ComputeSnap::compute_array()
   // accumulate bispectrum virial contributions to global array
 
   if (dbirjflag){
-
+    // no virial terms for dbirj yet
   }
   else{
     dbdotr_compute();
@@ -865,8 +693,8 @@ void ComputeSnap::compute_array()
   // assign energy to last column
   if (dbirjflag){
     // Assign reference energy right after the dbirj rows, first column.
-    int irow = bik_rows + dbirj_rows + 3*natoms; // add 3N for the dBi/dRi rows
+    // Add 3N for the dBi/dRi rows.
-    //printf("irow bik_rows dbirj_rows: %d %d %d\n", irow, bik_rows, dbirj_rows);
+    int irow = bik_rows + dbirj_rows + 3*natoms;
     double reference_energy = c_pe->compute_scalar();
     snapall[irow][0] = reference_energy;
   }
@@ -882,7 +710,7 @@ void ComputeSnap::compute_array()
 
   c_virial->compute_vector();
   if (dbirjflag){
-
+    // no virial terms for dbirj yet
   }
   else{
     int irow = 3*natoms+bik_rows;
@@ -894,7 +722,6 @@ void ComputeSnap::compute_array()
     snapall[irow][lastcol] = c_virial->vector[3];
   }
 
-  //}// else
 }
 
 /* ----------------------------------------------------------------------
@@ -906,14 +733,7 @@ void ComputeSnap::dbdotr_compute()
 {
   double **x = atom->x;
 
-  int irow0;
+  int irow0 = bik_rows+ndims_force*natoms;
-  if (dbirjflag){
-    irow0 = bik_rows+dbirj_rows+(3*natoms);
-  }
-  else{
-    irow0 = bik_rows+ndims_force*natoms;
-  }
-  //int irow0 = bik_rows+ndims_force*natoms;
 
   // sum over bispectrum contributions to forces
   // on all particles including ghosts
@@ -946,9 +766,11 @@ void ComputeSnap::dbdotr_compute()
 void ComputeSnap::get_dbirj_length()
 {
 
+  int rank = universe->me; // for MPI debugging
+
   memory->destroy(snap);
   memory->destroy(snapall);
-  // invoke full neighbor list (will copy or build if necessary)
+  // invoke full neighbor list
   neighbor->build_one(list);
   dbirj_rows = 0;
   const int inum = list->inum;
@@ -958,12 +780,20 @@ void ComputeSnap::get_dbirj_length()
   int * const type = atom->type;
   const int* const mask = atom->mask;
   double** const x = atom->x;
-  //printf("----- inum: %d\n", inum);
+  /*
   memory->create(neighsum, inum, "snap:neighsum");
   memory->create(nneighs, inum, "snap:nneighs");
   memory->create(icounter, inum, "snap:icounter");
   memory->create(dbiri, 3*atom->nlocal,ncoeff+3, "snap:dbiri");
-  for (int ii=0; ii<3*atom->nlocal; ii++){
+  */
+  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 dbirjflag=1 does not support parallelism.");
+  }
+  for (int ii=0; ii<3*natoms; ii++){
     for (int icoeff=0; icoeff<ncoeff; icoeff++){
       dbiri[ii][icoeff]=0.0;
     }
@@ -980,7 +810,6 @@ void ComputeSnap::get_dbirj_length()
       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];
@@ -998,7 +827,6 @@ void ComputeSnap::get_dbirj_length()
           nneighs[i]+=1;
         }
       }
-      //printf("----- jnum_cutoff: %d\n", jnum_cutoff);
     }
   }
 
@@ -1008,9 +836,6 @@ void ComputeSnap::get_dbirj_length()
   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;
       }
@@ -1018,18 +843,13 @@ void ComputeSnap::get_dbirj_length()
         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]);
   }
 
-  //printf("----- dbirj_rows: %d\n", dbirj_rows);
-
   memory->create(dbirj, dbirj_rows, ncoeff+3, "snap:dbirj");
   for (int i=0; i<dbirj_rows; i++){
     for (int j=0; j<ncoeff+3; j++){
@@ -1037,10 +857,7 @@ void ComputeSnap::get_dbirj_length()
     }
   }
   // Set size array rows which now depends on dbirj_rows.
-  //size_array_rows = bik_rows+dbirj_rows+ndims_virial;
   size_array_rows = bik_rows+dbirj_rows+ndims_virial+3*atom->nlocal; // Add 3*N for dBi/dRi
-  //printf("----- dbirj_rows: %d\n", dbirj_rows);
-  //printf("----- end of dbirj length.\n");
 
   memory->create(snap,size_array_rows,size_array_cols, "snap:snap");
   memory->create(snapall,size_array_rows,size_array_cols, "snap:snapall");