diff --git a/examples/snap/grid.py b/examples/snap/grid.py
new file mode 100755
index 0000000000..9fe88fd1b4
--- /dev/null
+++ b/examples/snap/grid.py
@@ -0,0 +1,104 @@
+#!/Users/athomps/miniconda3/bin/python3.7
+#
+# An example of SNAP grid from LAMMPS Python interface
+#
+# https://lammps.sandia.gov/doc/Python_library.html
+
+
+from lammps import lammps
+import lammps_utils
+
+# define command line input variables
+
+ngridx = 2
+ngridy = 3
+ngridz = 4
+twojmax = 2
+
+lmp_cmdargs = ["-echo","screen"]
+lmp_cmdargs = lammps_utils.set_cmdlinevars(lmp_cmdargs,
+    {
+        "ngridx":ngridx,
+        "ngridy":ngridy,
+        "ngridz":ngridz,
+        "twojmax":twojmax
+        }
+    )
+
+# launch LAMMPS instance
+
+lmp = lammps(cmdargs=lmp_cmdargs)
+
+# run LAMMPS input script
+
+lmp.file("in.grid.python")
+
+# get quantities from LAMMPS 
+
+num_atoms = lmp.get_natoms() 
+
+# set things not accessible from LAMMPS
+
+# first 3 cols are x, y, z, coords
+
+ncols0 = 3 
+
+# analytical relation
+
+ncoeff = (twojmax+2)*(twojmax+3)*(twojmax+4)
+ncoeff = ncoeff // 24 # integer division
+ncols = ncols0+ncoeff
+
+# get B_0 at position (0,0,0) in 4 different ways
+
+# 1. from comute sna/atom
+
+bptr = lmp.extract_compute("b", 1, 2) # 1 = per-atom data, 2 = array
+print("b = ",bptr[0][0])
+
+# 2. from compute sna/grid
+
+bgridptr = lmp.extract_compute("bgrid", 0, 2) # 0 = style global, 2 = type array
+print("bgrid = ",bgridptr[0][3])
+
+# 3. from Numpy array pointing to sna/atom array
+ 
+bptr_np = lammps_utils.extract_compute_np(lmp,"b",1,2,(num_atoms,ncoeff))
+print("b_np = ",bptr_np[0][0])
+
+# 4. from Numpy array pointing to sna/grid array
+
+bgridptr_np = lammps_utils.extract_compute_np(lmp,"bgrid",0,2,(ngridz,ngridy,ngridx,ncols))
+print("bgrid_np = ",bgridptr_np[0][0][0][ncols0+0])
+
+# print out the LAMMPS array to a file
+
+outfile = open("bgrid.dat",'w')
+igrid = 0
+for iz in range(ngridz):
+    for iy in range(ngridy):
+        for ix in range(ngridx):
+            outfile.write("x, y, z = %g %g %g\n" % 
+                          (bgridptr[igrid][0],
+                           bgridptr[igrid][1],
+                           bgridptr[igrid][2]))
+            for icoeff in range(ncoeff):
+                outfile.write("%g " % bgridptr[igrid][ncols0+icoeff])
+            outfile.write("\n")
+            igrid += 1
+outfile.close()
+
+# print out the Numpy array to a file
+
+outfile = open("bgrid_np.dat",'w')
+for iz in range(ngridz):
+    for iy in range(ngridy):
+        for ix in range(ngridx):
+            outfile.write("x, y, z = %g %g %g\n" % 
+                          (bgridptr_np[iz][iy][ix][0],
+                           bgridptr_np[iz][iy][ix][1],
+                           bgridptr_np[iz][iy][ix][2]))
+            for icoeff in range(ncoeff):
+                outfile.write("%g " % bgridptr_np[iz][iy][ix][ncols0+icoeff])
+            outfile.write("\n")
+outfile.close()
diff --git a/examples/snap/in.grid b/examples/snap/in.grid
new file mode 100644
index 0000000000..44c42673e0
--- /dev/null
+++ b/examples/snap/in.grid
@@ -0,0 +1,83 @@
+# Demonstrate bispectrum computes
+
+# CORRECTNESS: thermo output for c_mygrid[*][1] and c_mygrid[*][8] should
+# match the values in dump_b: 108.173 3.21778 0.712238 7.06634 1.04273 
+  
+# Initialize simulation
+
+variable nsteps index 0
+variable nrep index 1
+variable a index 3.316
+variable ngrid index 2
+
+units		metal
+
+# generate the box and atom positions using a BCC lattice
+
+variable nx equal ${nrep}
+variable ny equal ${nrep}
+variable nz equal ${nrep}
+
+boundary	p p p
+
+lattice		custom $a &
+		a1 1 0 0 &
+		a2 0 1 0  &
+		a3 0 0 1 &
+		basis 0 0 0 &
+		basis 0.5 0.5 0.5 &
+#		origin 0.25 0.25 0.25
+
+region		box block 0 ${nx} 0 ${ny} 0 ${nz}
+create_box	1 box
+create_atoms	1 box
+
+mass 1 180.88
+
+# choose potential
+
+include Ta06A.snap
+
+# define grid compute and atom compute
+
+group 		snapgroup type 1
+variable 	twojmax equal 2
+variable 	rcutfac equal 4.67637
+variable 	rfac0 equal 0.99363
+variable 	rmin0 equal 0
+variable 	wj equal 1
+variable 	radelem equal 0.5
+variable 	bzero equal 0
+variable 	quad equal 0
+variable 	switch equal 1
+
+compute b all sna/atom &
+	${rcutfac} ${rfac0} ${twojmax} ${radelem} &
+	${wj} rmin0 ${rmin0} bzeroflag ${bzero} &
+	quadraticflag ${quad} switchflag ${switch}
+
+compute mygrid all sna/grid grid ${ngrid} ${ngrid} ${ngrid} &
+	${rcutfac} ${rfac0} ${twojmax} ${radelem} &
+	${wj} rmin0 ${rmin0} bzeroflag ${bzero} &
+	quadraticflag ${quad} switchflag ${switch}
+
+# define output
+
+# mygrid is ngrid by (3+nbis) = 8x8
+thermo_style	custom step temp ke pe vol &
+	c_mygrid[1][1] c_mygrid[2][1] c_mygrid[3][1] c_mygrid[4][1] c_mygrid[5][1] c_mygrid[6][1] c_mygrid[7][1] c_mygrid[8][1] &
+	c_mygrid[1][2] c_mygrid[2][2] c_mygrid[3][2] c_mygrid[4][2] c_mygrid[5][2] c_mygrid[6][2] c_mygrid[7][2] c_mygrid[8][2] &
+	c_mygrid[1][3] c_mygrid[2][3] c_mygrid[3][3] c_mygrid[4][3] c_mygrid[5][3] c_mygrid[6][3] c_mygrid[7][3] c_mygrid[8][3] &
+	c_mygrid[1][4] c_mygrid[2][4] c_mygrid[3][4] c_mygrid[4][4] c_mygrid[5][4] c_mygrid[6][4] c_mygrid[7][4] c_mygrid[8][4] &
+	c_mygrid[1][5] c_mygrid[2][5] c_mygrid[3][5] c_mygrid[4][5] c_mygrid[5][5] c_mygrid[6][5] c_mygrid[7][5] c_mygrid[8][5] &
+	c_mygrid[1][6] c_mygrid[2][6] c_mygrid[3][6] c_mygrid[4][6] c_mygrid[5][6] c_mygrid[6][6] c_mygrid[7][6] c_mygrid[8][6] &
+	c_mygrid[1][7] c_mygrid[2][7] c_mygrid[3][7] c_mygrid[4][7] c_mygrid[5][7] c_mygrid[6][7] c_mygrid[7][7] c_mygrid[8][7] &
+	c_mygrid[1][8] c_mygrid[2][8] c_mygrid[3][8] c_mygrid[4][8] c_mygrid[5][8] c_mygrid[6][8] c_mygrid[7][8] c_mygrid[8][8]
+thermo_modify norm yes
+
+dump mydump_b all custom 1000 dump_b id c_b[*]
+
+# run
+
+run 0
+
diff --git a/examples/snap/in.grid.python b/examples/snap/in.grid.python
new file mode 100644
index 0000000000..13fe6d0638
--- /dev/null
+++ b/examples/snap/in.grid.python
@@ -0,0 +1,61 @@
+# Demonstrate bispectrum per-atom and grid computes
+
+# Invoked from grid.py
+# pass in values for ngridx, ngridy, ngridz, twojmax
+
+# Initialize simulation
+
+variable nsteps equal 0
+variable nrep equal 1
+variable a equal 3.316
+
+units		metal
+
+# generate the box and atom positions using a BCC lattice
+
+variable nx equal ${nrep}
+variable ny equal ${nrep}
+variable nz equal ${nrep}
+
+boundary	p p p
+
+lattice		custom ${a} a1 1 0 0 a2 0 1 0 a3 0 0 1 basis 0 0 0 basis 0.5 0.5 0.5
+
+region		box block 0 ${nx} 0 ${ny} 0 ${nz}
+create_box	1 box
+create_atoms	1 box
+
+mass 1 180.88
+
+# define grid compute and atom compute
+
+group 		snapgroup type 1
+variable 	rcutfac equal 4.67637
+variable 	rfac0 equal 0.99363
+variable 	rmin0 equal 0
+variable 	wj equal 1
+variable 	radelem equal 0.5
+variable 	bzero equal 0
+variable 	quad equal 0
+variable 	switch equal 1
+
+compute b all sna/atom ${rcutfac} ${rfac0} ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} bzeroflag ${bzero} quadraticflag ${quad} switchflag ${switch}
+
+compute bgrid all sna/grid grid ${ngridx} ${ngridy} ${ngridz} ${rcutfac} ${rfac0} ${twojmax} ${radelem} ${wj} rmin0 ${rmin0} bzeroflag ${bzero}	quadraticflag ${quad} switchflag ${switch}
+
+# create dummy potential for neighbor list
+
+variable rcutneigh equal 2.0*${rcutfac}*${radelem}
+pair_style zero ${rcutneigh}
+pair_coeff * *
+
+# define output
+
+thermo_style	custom step temp ke pe vol c_bgrid[1][1]
+thermo_modify norm yes
+
+dump mydump_b all custom 1000 dump_b id c_b[*]
+
+# run
+
+run 0
diff --git a/examples/snap/in.grid.tri b/examples/snap/in.grid.tri
new file mode 100644
index 0000000000..f375c18c42
--- /dev/null
+++ b/examples/snap/in.grid.tri
@@ -0,0 +1,135 @@
+# Demonstrate bispectrum computes for triclinic cell
+
+# This triclinic cell has 6 times the volume of the single
+# unit cell used by in.grid
+# and contains 12 atoms.  It is a 3x2x1 supercell
+# with each unit cell containing 2 atoms and the
+# reduced lattice vectors are [1 0 0], [1 1 0], and [1 1 1].
+# The grid is listed in x-fastest order
+# CORRECTNESS: thermo output for c_mygrid[*][7] should contain
+# the same 2 values that appear in c_mygrid[*][7] for in.grid,
+# 7.0663376 and 13.808803, corresponding to atom and insterstitial sites,
+# respectively, and with occurrences 12 and 36, respectively.
+
+# Initialize simulation
+
+variable nsteps index 0
+variable nrep index 1
+variable a index 3.316
+variable ngrid index 4
+
+variable ngridx equal 3*${ngrid}
+variable ngridy equal 2*${ngrid}
+variable ngridz equal 1*${ngrid}
+
+variable nrepx equal 1*${nrep}
+variable nrepy equal 1*${nrep}
+variable nrepz equal 1*${nrep}
+
+units		metal
+
+# generate the box and atom positions using a BCC lattice
+
+variable nx equal ${nrepx}
+variable ny equal ${nrepy}
+variable nz equal ${nrepz}
+
+boundary	p p p
+
+lattice		custom $a &
+		a1 1 0 0 &
+		a2 1 1 0  &
+		a3 1 1 1 &
+		basis 0 0 0 &
+		basis 0.0 0.0 0.5 &
+#		origin 0.25 0.25 0.25
+
+box 		tilt large
+region		box prism 0 ${nx} 0 ${ny} 0 ${nz} ${ny} ${nz} ${nz}
+create_box	1 box
+create_atoms	1 box
+
+mass 1 180.88
+
+# choose potential
+
+include Ta06A.snap
+
+# define grid compute and atom compute
+
+group 		snapgroup type 1
+variable 	twojmax equal 2
+variable 	rcutfac equal 4.67637
+variable 	rfac0 equal 0.99363
+variable 	rmin0 equal 0
+variable 	wj equal 1
+variable 	radelem equal 0.5
+variable 	bzero equal 0
+variable 	quad equal 0
+variable 	switch equal 1
+
+compute b all sna/atom &
+	${rcutfac} ${rfac0} ${twojmax} ${radelem} &
+	${wj} rmin0 ${rmin0} bzeroflag ${bzero} &
+	quadraticflag ${quad} switchflag ${switch}
+
+compute mygrid all sna/grid grid ${ngridx} ${ngridy} ${ngridz} &
+	${rcutfac} ${rfac0} ${twojmax} ${radelem} &
+	${wj} rmin0 ${rmin0} bzeroflag ${bzero} &
+	quadraticflag ${quad} switchflag ${switch}
+
+# define output
+
+# mygrid is ngrid by (3+nbis) = 384x8
+
+thermo_style	custom step temp ke pe vol &
+		c_mygrid[1][7] c_mygrid[2][7] c_mygrid[3][7] c_mygrid[4][7] c_mygrid[5][7] c_mygrid[6][7] c_mygrid[7][7] c_mygrid[8][7] c_mygrid[9][7] &
+		c_mygrid[10][7] c_mygrid[11][7] c_mygrid[12][7] c_mygrid[13][7] c_mygrid[14][7] c_mygrid[15][7] c_mygrid[16][7] c_mygrid[17][7] c_mygrid[18][7] c_mygrid[19][7] &
+ 		c_mygrid[20][7] c_mygrid[21][7] c_mygrid[22][7] c_mygrid[23][7] c_mygrid[24][7] c_mygrid[25][7] c_mygrid[26][7] c_mygrid[27][7] c_mygrid[28][7] c_mygrid[29][7] &
+ 		c_mygrid[30][7] c_mygrid[31][7] c_mygrid[32][7] c_mygrid[33][7] c_mygrid[34][7] c_mygrid[35][7] c_mygrid[36][7] c_mygrid[37][7] c_mygrid[38][7] c_mygrid[39][7] &
+ 		c_mygrid[40][7] c_mygrid[41][7] c_mygrid[42][7] c_mygrid[43][7] c_mygrid[44][7] c_mygrid[45][7] c_mygrid[46][7] c_mygrid[47][7] c_mygrid[48][7] c_mygrid[49][7] &
+ 		c_mygrid[50][7] c_mygrid[51][7] c_mygrid[52][7] c_mygrid[53][7] c_mygrid[54][7] c_mygrid[55][7] c_mygrid[56][7] c_mygrid[57][7] c_mygrid[58][7] c_mygrid[59][7] &
+ 		c_mygrid[60][7] c_mygrid[61][7] c_mygrid[62][7] c_mygrid[63][7] c_mygrid[64][7] c_mygrid[65][7] c_mygrid[66][7] c_mygrid[67][7] c_mygrid[68][7] c_mygrid[69][7] &
+ 		c_mygrid[70][7] c_mygrid[71][7] c_mygrid[72][7] c_mygrid[73][7] c_mygrid[74][7] c_mygrid[75][7] c_mygrid[76][7] c_mygrid[77][7] c_mygrid[78][7] c_mygrid[79][7] &
+ 		c_mygrid[80][7] c_mygrid[81][7] c_mygrid[82][7] c_mygrid[83][7] c_mygrid[84][7] c_mygrid[85][7] c_mygrid[86][7] c_mygrid[87][7] c_mygrid[88][7] c_mygrid[89][7] &
+ 		c_mygrid[90][7] c_mygrid[91][7] c_mygrid[92][7] c_mygrid[93][7] c_mygrid[94][7] c_mygrid[95][7] c_mygrid[96][7] c_mygrid[97][7] c_mygrid[98][7] c_mygrid[99][7] &
+		c_mygrid[100][7] c_mygrid[101][7] c_mygrid[102][7] c_mygrid[103][7] c_mygrid[104][7] c_mygrid[105][7] c_mygrid[106][7] c_mygrid[107][7] c_mygrid[108][7] c_mygrid[109][7] &
+		c_mygrid[110][7] c_mygrid[111][7] c_mygrid[112][7] c_mygrid[113][7] c_mygrid[114][7] c_mygrid[115][7] c_mygrid[116][7] c_mygrid[117][7] c_mygrid[118][7] c_mygrid[119][7] &
+ 		c_mygrid[120][7] c_mygrid[121][7] c_mygrid[122][7] c_mygrid[123][7] c_mygrid[124][7] c_mygrid[125][7] c_mygrid[126][7] c_mygrid[127][7] c_mygrid[128][7] c_mygrid[129][7] &
+ 		c_mygrid[130][7] c_mygrid[131][7] c_mygrid[132][7] c_mygrid[133][7] c_mygrid[134][7] c_mygrid[135][7] c_mygrid[136][7] c_mygrid[137][7] c_mygrid[138][7] c_mygrid[139][7] &
+ 		c_mygrid[140][7] c_mygrid[141][7] c_mygrid[142][7] c_mygrid[143][7] c_mygrid[144][7] c_mygrid[145][7] c_mygrid[146][7] c_mygrid[147][7] c_mygrid[148][7] c_mygrid[149][7] &
+ 		c_mygrid[150][7] c_mygrid[151][7] c_mygrid[152][7] c_mygrid[153][7] c_mygrid[154][7] c_mygrid[155][7] c_mygrid[156][7] c_mygrid[157][7] c_mygrid[158][7] c_mygrid[159][7] &
+ 		c_mygrid[160][7] c_mygrid[161][7] c_mygrid[162][7] c_mygrid[163][7] c_mygrid[164][7] c_mygrid[165][7] c_mygrid[166][7] c_mygrid[167][7] c_mygrid[168][7] c_mygrid[169][7] &
+ 		c_mygrid[170][7] c_mygrid[171][7] c_mygrid[172][7] c_mygrid[173][7] c_mygrid[174][7] c_mygrid[175][7] c_mygrid[176][7] c_mygrid[177][7] c_mygrid[178][7] c_mygrid[179][7] &
+ 		c_mygrid[180][7] c_mygrid[181][7] c_mygrid[182][7] c_mygrid[183][7] c_mygrid[184][7] c_mygrid[185][7] c_mygrid[186][7] c_mygrid[187][7] c_mygrid[188][7] c_mygrid[189][7] &
+ 		c_mygrid[190][7] c_mygrid[191][7] c_mygrid[192][7] c_mygrid[193][7] c_mygrid[194][7] c_mygrid[195][7] c_mygrid[196][7] c_mygrid[197][7] c_mygrid[198][7] c_mygrid[199][7] &
+		c_mygrid[200][7] c_mygrid[201][7] c_mygrid[202][7] c_mygrid[203][7] c_mygrid[204][7] c_mygrid[205][7] c_mygrid[206][7] c_mygrid[207][7] c_mygrid[208][7] c_mygrid[209][7] &
+		c_mygrid[210][7] c_mygrid[211][7] c_mygrid[212][7] c_mygrid[213][7] c_mygrid[214][7] c_mygrid[215][7] c_mygrid[216][7] c_mygrid[217][7] c_mygrid[218][7] c_mygrid[219][7] &
+ 		c_mygrid[220][7] c_mygrid[221][7] c_mygrid[222][7] c_mygrid[223][7] c_mygrid[224][7] c_mygrid[225][7] c_mygrid[226][7] c_mygrid[227][7] c_mygrid[228][7] c_mygrid[229][7] &
+ 		c_mygrid[230][7] c_mygrid[231][7] c_mygrid[232][7] c_mygrid[233][7] c_mygrid[234][7] c_mygrid[235][7] c_mygrid[236][7] c_mygrid[237][7] c_mygrid[238][7] c_mygrid[239][7] &
+ 		c_mygrid[240][7] c_mygrid[241][7] c_mygrid[242][7] c_mygrid[243][7] c_mygrid[244][7] c_mygrid[245][7] c_mygrid[246][7] c_mygrid[247][7] c_mygrid[248][7] c_mygrid[249][7] &
+ 		c_mygrid[250][7] c_mygrid[251][7] c_mygrid[252][7] c_mygrid[253][7] c_mygrid[254][7] c_mygrid[255][7] c_mygrid[256][7] c_mygrid[257][7] c_mygrid[258][7] c_mygrid[259][7] &
+ 		c_mygrid[260][7] c_mygrid[261][7] c_mygrid[262][7] c_mygrid[263][7] c_mygrid[264][7] c_mygrid[265][7] c_mygrid[266][7] c_mygrid[267][7] c_mygrid[268][7] c_mygrid[269][7] &
+ 		c_mygrid[270][7] c_mygrid[271][7] c_mygrid[272][7] c_mygrid[273][7] c_mygrid[274][7] c_mygrid[275][7] c_mygrid[276][7] c_mygrid[277][7] c_mygrid[278][7] c_mygrid[279][7] &
+ 		c_mygrid[280][7] c_mygrid[281][7] c_mygrid[282][7] c_mygrid[283][7] c_mygrid[284][7] c_mygrid[285][7] c_mygrid[286][7] c_mygrid[287][7] c_mygrid[288][7] c_mygrid[289][7] &
+ 		c_mygrid[290][7] c_mygrid[291][7] c_mygrid[292][7] c_mygrid[293][7] c_mygrid[294][7] c_mygrid[295][7] c_mygrid[296][7] c_mygrid[297][7] c_mygrid[298][7] c_mygrid[299][7] &
+		c_mygrid[300][7] c_mygrid[301][7] c_mygrid[302][7] c_mygrid[303][7] c_mygrid[304][7] c_mygrid[305][7] c_mygrid[306][7] c_mygrid[307][7] c_mygrid[308][7] c_mygrid[309][7] &
+		c_mygrid[310][7] c_mygrid[311][7] c_mygrid[312][7] c_mygrid[313][7] c_mygrid[314][7] c_mygrid[315][7] c_mygrid[316][7] c_mygrid[317][7] c_mygrid[318][7] c_mygrid[319][7] &
+ 		c_mygrid[320][7] c_mygrid[321][7] c_mygrid[322][7] c_mygrid[323][7] c_mygrid[324][7] c_mygrid[325][7] c_mygrid[326][7] c_mygrid[327][7] c_mygrid[328][7] c_mygrid[329][7] &
+ 		c_mygrid[330][7] c_mygrid[331][7] c_mygrid[332][7] c_mygrid[333][7] c_mygrid[334][7] c_mygrid[335][7] c_mygrid[336][7] c_mygrid[337][7] c_mygrid[338][7] c_mygrid[339][7] &
+ 		c_mygrid[340][7] c_mygrid[341][7] c_mygrid[342][7] c_mygrid[343][7] c_mygrid[344][7] c_mygrid[345][7] c_mygrid[346][7] c_mygrid[347][7] c_mygrid[348][7] c_mygrid[349][7] &
+ 		c_mygrid[350][7] c_mygrid[351][7] c_mygrid[352][7] c_mygrid[353][7] c_mygrid[354][7] c_mygrid[355][7] c_mygrid[356][7] c_mygrid[357][7] c_mygrid[358][7] c_mygrid[359][7] &
+ 		c_mygrid[360][7] c_mygrid[361][7] c_mygrid[362][7] c_mygrid[363][7] c_mygrid[364][7] c_mygrid[365][7] c_mygrid[366][7] c_mygrid[367][7] c_mygrid[368][7] c_mygrid[369][7] &
+ 		c_mygrid[370][7] c_mygrid[371][7] c_mygrid[372][7] c_mygrid[373][7] c_mygrid[374][7] c_mygrid[375][7] c_mygrid[376][7] c_mygrid[377][7] c_mygrid[378][7] c_mygrid[379][7] &
+ 		c_mygrid[380][7] c_mygrid[381][7] c_mygrid[382][7] c_mygrid[383][7] c_mygrid[384][7]
+
+thermo_modify norm yes
+
+dump mydump_b all custom 1000 dump_b id c_b[*]
+
+# run
+
+run 0
+
+
+
diff --git a/src/SNAP/compute_sna_grid.cpp b/src/SNAP/compute_sna_grid.cpp
new file mode 100644
index 0000000000..d2c2ae74ca
--- /dev/null
+++ b/src/SNAP/compute_sna_grid.cpp
@@ -0,0 +1,264 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, 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_grid.h"
+#include "compute_sna_grid.h"
+#include <cstring>
+#include <cstdlib>
+#include "sna.h"
+#include "atom.h"
+#include "update.h"
+#include "modify.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "force.h"
+#include "pair.h"
+#include "comm.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+ComputeSNAGrid::ComputeSNAGrid(LAMMPS *lmp, int narg, char **arg) :
+  ComputeGrid(lmp, narg, arg), cutsq(NULL), sna(NULL),
+  radelem(NULL), wjelem(NULL)
+{
+  double rmin0, rfac0;
+  int twojmax, switchflag, bzeroflag;
+  radelem = NULL;
+  wjelem = NULL;
+
+  // skip over arguments used by base class
+  // so that argument positions are identical to 
+  // regular per-atom compute
+
+  arg += nargbase;
+  narg -= nargbase;
+
+  int ntypes = atom->ntypes;
+  int nargmin = 6+2*ntypes;
+
+  if (narg < nargmin) error->all(FLERR,"Illegal compute sna/grid command");
+
+  // default values
+
+  rmin0 = 0.0;
+  switchflag = 1;
+  bzeroflag = 1;
+  quadraticflag = 0;
+
+  // offset by 1 to match up with types
+
+  memory->create(radelem,ntypes+1,"sna/grid:radelem");
+  memory->create(wjelem,ntypes+1,"sna/grid: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,"sna/grid:cutsq");
+  for(int i = 1; i <= ntypes; i++) {
+    cut = 2.0*radelem[i]*rcutfac;
+    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;
+    }
+  }
+
+  // process optional args
+
+  int iarg = nargmin;
+
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"rmin0") == 0) {
+      if (iarg+2 > narg)
+        error->all(FLERR,"Illegal compute sna/grid command");
+      rmin0 = atof(arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"switchflag") == 0) {
+      if (iarg+2 > narg)
+        error->all(FLERR,"Illegal compute sna/grid command");
+      switchflag = atoi(arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"bzeroflag") == 0) {
+      if (iarg+2 > narg)
+        error->all(FLERR,"Illegal compute sna/grid command");
+      bzeroflag = atoi(arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"quadraticflag") == 0) {
+      if (iarg+2 > narg)
+        error->all(FLERR,"Illegal compute sna/grid command");
+      quadraticflag = atoi(arg[iarg+1]);
+      iarg += 2;
+    } else error->all(FLERR,"Illegal compute sna/grid command");
+
+  }
+
+  snaptr = new SNA(lmp,rfac0,twojmax,
+                   rmin0,switchflag,bzeroflag);
+
+  ncoeff = snaptr->ncoeff;
+  nvalues = ncoeff;
+  if (quadraticflag) nvalues += (ncoeff*(ncoeff+1))/2;
+  size_array_cols = size_array_cols_base + nvalues;
+  array_flag = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeSNAGrid::~ComputeSNAGrid()
+{
+  memory->destroy(sna);
+  memory->destroy(radelem);
+  memory->destroy(wjelem);
+  memory->destroy(cutsq);
+  delete snaptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeSNAGrid::init()
+{
+  if (force->pair == NULL)
+    error->all(FLERR,"Compute sna/grid requires a pair style be defined");
+
+  if (cutmax > force->pair->cutforce)
+    error->all(FLERR,"Compute sna/grid cutoff is longer than pairwise cutoff");
+
+  // need an occasional full neighbor list
+
+  int irequest = neighbor->request(this,instance_me);
+  neighbor->requests[irequest]->pair = 0;
+  neighbor->requests[irequest]->compute = 1;
+  neighbor->requests[irequest]->half = 0;
+  neighbor->requests[irequest]->full = 1;
+  neighbor->requests[irequest]->occasional = 1;
+
+  int count = 0;
+  for (int i = 0; i < modify->ncompute; i++)
+    if (strcmp(modify->compute[i]->style,"sna/grid") == 0) count++;
+  if (count > 1 && comm->me == 0)
+    error->warning(FLERR,"More than one compute sna/grid");
+  snaptr->init();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeSNAGrid::compute_array()
+{
+  invoked_array = update->ntimestep;
+
+  int * const type = atom->type;
+
+  // compute sna for each gridpoint
+
+  double** const x = atom->x;
+  const int* const mask = atom->mask;
+  const int ntotal = atom->nlocal + atom->nghost;
+
+  // insure rij, inside, and typej are of size jnum
+  
+  snaptr->grow_rij(ntotal);
+
+  for (int iz = nzlo; iz <= nzhi; iz++)
+    for (int iy = nylo; iy <= nyhi; iy++)
+      for (int ix = nxlo; ix <= nxhi; ix++) {
+	const int igrid = iz*(nx*ny) + iy*nx + ix;
+	const double xtmp = grid[igrid][0];
+	const double ytmp = grid[igrid][1];
+	const double ztmp = grid[igrid][2];
+
+	// 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
+
+	int ninside = 0;
+	for (int j = 0; j < ntotal; j++) {
+
+	  // check that j is in compute group
+
+	  if (!(mask[j] & groupbit)) continue;
+
+	  const double delx = xtmp - x[j][0];
+	  const double dely = ytmp - x[j][1];
+	  const double delz = ztmp - x[j][2];
+	  const double rsq = delx*delx + dely*dely + delz*delz;
+	  int jtype = type[j];
+	  if (rsq < cutsq[jtype][jtype] && rsq>1e-20) {
+	    snaptr->rij[ninside][0] = delx;
+	    snaptr->rij[ninside][1] = dely;
+	    snaptr->rij[ninside][2] = delz;
+	    snaptr->inside[ninside] = j;
+	    snaptr->wj[ninside] = wjelem[jtype];
+	    snaptr->rcutij[ninside] = 2.0*radelem[jtype]*rcutfac;
+	    ninside++;
+	  }
+	}
+
+	snaptr->compute_ui(ninside);
+	snaptr->compute_zi();
+	snaptr->compute_bi();
+	for (int icoeff = 0; icoeff < ncoeff; icoeff++)
+	  gridlocal[size_array_cols_base+icoeff][iz][iy][ix] = snaptr->blist[icoeff];
+	if (quadraticflag) {
+	  int ncount = ncoeff;
+	  for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+	    double bi = snaptr->blist[icoeff];
+
+	    // diagonal element of quadratic matrix
+
+	    gridlocal[size_array_cols_base+ncount++][iz][iy][ix] = 0.5*bi*bi;
+
+	    // upper-triangular elements of quadratic matrix
+
+	    for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++)
+	      gridlocal[size_array_cols_base+ncount++][iz][iy][ix] = bi*snaptr->blist[jcoeff];
+      }
+    }
+  }
+
+  for (int iz = nzlo; iz <= nzhi; iz++)
+    for (int iy = nylo; iy <= nyhi; iy++)
+      for (int ix = nxlo; ix <= nxhi; ix++) {
+	const int igrid = iz*(nx*ny) + iy*nx + ix;
+	for (int j = 0; j < nvalues; j++)
+	  grid[igrid][size_array_cols_base + j] = gridlocal[size_array_cols_base + j][iz][iy][ix];
+      }
+  MPI_Allreduce(&grid[0][0],&gridall[0][0],ngrid*size_array_cols,MPI_DOUBLE,MPI_SUM,world);
+
+}
+
+
+/* ----------------------------------------------------------------------
+   memory usage
+------------------------------------------------------------------------- */
+
+double ComputeSNAGrid::memory_usage()
+{
+  double nbytes = snaptr->memory_usage();                        // SNA object
+
+  return nbytes;
+}
+
diff --git a/src/SNAP/compute_sna_grid.h b/src/SNAP/compute_sna_grid.h
new file mode 100644
index 0000000000..0242a2962b
--- /dev/null
+++ b/src/SNAP/compute_sna_grid.h
@@ -0,0 +1,71 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, 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
+
+ComputeStyle(sna/grid,ComputeSNAGrid)
+
+#else
+
+#ifndef LMP_COMPUTE_SNA_GRID_H
+#define LMP_COMPUTE_SNA_GRID_H
+
+#include "compute_grid.h"
+
+namespace LAMMPS_NS {
+
+class ComputeSNAGrid : public ComputeGrid {
+ public:
+  ComputeSNAGrid(class LAMMPS *, int, char **);
+  ~ComputeSNAGrid();
+  void init();
+  void compute_array();
+  double memory_usage();
+
+ private:
+  int ncoeff;
+  double **cutsq;
+  double **sna;
+  double rcutfac;
+  double *radelem;
+  double *wjelem;
+  class SNA* snaptr;
+  int quadraticflag;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Compute sna/grid requires a pair style be defined
+
+Self-explanatory.
+
+E: Compute sna/grid cutoff is longer than pairwise cutoff
+
+Self-explanatory.
+
+W: More than one compute sna/grid
+
+Self-explanatory.
+
+*/
diff --git a/src/compute_grid.cpp b/src/compute_grid.cpp
new file mode 100644
index 0000000000..3cbcb2c2db
--- /dev/null
+++ b/src/compute_grid.cpp
@@ -0,0 +1,289 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, 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_grid.h"
+#include <mpi.h>
+#include <cstring>
+#include "atom.h"
+#include "update.h"
+#include "modify.h"
+#include "domain.h"
+#include "force.h"
+#include "memory.h"
+#include "error.h"
+#include "comm.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeGrid::ComputeGrid(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg), grid(NULL), local_flags(NULL), gridlocal(NULL)
+{
+  if (narg < 6) error->all(FLERR,"Illegal compute grid command");
+
+  array_flag = 1;
+  size_array_cols = 0;
+  size_array_rows = 0;
+  extarray = 0;
+
+  int iarg0 = 3;
+  int iarg = iarg0;
+  if (strcmp(arg[iarg],"grid") == 0) {
+    if (iarg+4 > narg) error->all(FLERR,"Illegal compute grid command");
+    nx = force->inumeric(FLERR,arg[iarg+1]);
+    ny = force->inumeric(FLERR,arg[iarg+2]);
+    nz = force->inumeric(FLERR,arg[iarg+3]);
+    if (nx <= 0 || ny <= 0 || nz <= 0)
+      error->all(FLERR,"All grid dimensions must be positive");
+    iarg += 4;
+  } else error->all(FLERR,"Illegal compute grid command");
+
+  nargbase = iarg - iarg0;
+
+  size_array_rows = ngrid = nx*ny*nz;
+  size_array_cols_base = 3;
+  gridlocal_allocated = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeGrid::~ComputeGrid()
+{
+  memory->destroy(grid);
+  memory->destroy(local_flags);
+  if (gridlocal_allocated)
+    memory->destroy4d_offset(gridlocal,nzlo,nylo,nxlo);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeGrid::init()
+{
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeGrid::setup()
+{
+  set_grid_global();
+  set_grid_local();
+  allocate();
+  assign_coords();
+  assign_local_flags();
+}
+
+/* ----------------------------------------------------------------------
+   convert global array index to box coords
+------------------------------------------------------------------------- */
+
+void ComputeGrid::grid2x(int igrid, double *x)
+{
+  int iz = igrid / (nx*ny);
+  igrid -= iz * (nx*ny);
+  int iy = igrid / nx;
+  igrid -= iy * nx;
+  int ix = igrid;
+
+  x[0] = ix*delx;
+  x[1] = iy*dely;
+  x[2] = iz*delz;
+
+  if (triclinic) domain->lamda2x(x, x);
+
+}
+
+/* ----------------------------------------------------------------------
+   convert global array index to box coords
+------------------------------------------------------------------------- */
+
+void ComputeGrid::grid2ix(int igrid, int& ix, int& iy, int& iz)
+{
+  iz = igrid / (nx*ny);
+  igrid -= iz * (nx*ny);
+  iy = igrid / nx;
+  igrid -= iy * nx;
+  ix = igrid;
+}
+
+// /* ----------------------------------------------------------------------
+//    check if grid point is local
+// ------------------------------------------------------------------------- */
+
+// int ComputeGrid::check_local(int igrid)
+// {
+//   double x[3];
+
+//   int iz = igrid / (nx*ny);
+//   igrid -= iz * (nx*ny);
+//   int iy = igrid / nx;
+//   igrid -= iy * nx;
+//   int ix = igrid;
+
+//   x[0] = ix*delx;
+//   x[1] = iy*dely;
+//   x[2] = iz*delz;
+
+//   int islocal = 
+//     x[0] >= sublo[0] && x[0] < subhi[0] &&
+//     x[1] >= sublo[1] && x[1] < subhi[1] &&
+//     x[2] >= sublo[2] && x[2] < subhi[2];
+
+//   return islocal;
+// }
+
+/* ----------------------------------------------------------------------
+   check if grid point is local
+------------------------------------------------------------------------- */
+
+int ComputeGrid::check_local(int igrid)
+{
+  int ix, iy, iz;
+
+  grid2ix(igrid, ix, iy, iz);
+
+  int islocal = 
+    ix >= nxlo && ix <= nxhi &&
+    iy >= nylo && iy <= nyhi &&
+    iz >= nzlo && iz <= nzhi;
+
+  return islocal;
+}
+
+/* ----------------------------------------------------------------------
+   copy coords to global array
+------------------------------------------------------------------------- */
+
+void ComputeGrid::assign_coords()
+{
+  double x[3];
+  for (int igrid = 0; igrid < ngrid; igrid++) {
+    grid2x(igrid,x);
+    grid[igrid][0] = x[0];
+    grid[igrid][1] = x[1];
+    grid[igrid][2] = x[2];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   copy coords to global array
+------------------------------------------------------------------------- */
+
+void ComputeGrid::assign_local_flags()
+{
+  for (int igrid = 0; igrid < ngrid; igrid++) {
+    if (check_local(igrid))
+      local_flags[igrid] = 1;
+    else {
+      local_flags[igrid] = 0;
+      memset(grid[igrid],0,size_array_cols * sizeof(double));
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   free and reallocate arrays
+------------------------------------------------------------------------- */
+
+void ComputeGrid::allocate()
+{
+  // allocate arrays
+
+  memory->destroy(grid);
+  memory->destroy(local_flags);
+  if (gridlocal_allocated)
+    memory->destroy4d_offset(gridlocal,nzlo,nylo,nxlo);
+  memory->create(grid,size_array_rows,size_array_cols,"grid:grid");
+  memory->create(gridall,size_array_rows,size_array_cols,"grid:gridall");
+  memory->create(local_flags,size_array_rows,"grid:local_flags");
+  if (nxlo <= nxhi && nylo <= nyhi && nzlo <= nzhi) {
+    gridlocal_allocated = 1;
+    memory->create4d_offset(gridlocal,size_array_cols,nzlo,nzhi,nylo,nyhi,
+			    nxlo,nxhi,"grid:gridlocal");
+  }
+  array = gridall;
+}
+
+
+/* ----------------------------------------------------------------------
+   set global grid
+------------------------------------------------------------------------- */
+
+void ComputeGrid::set_grid_global()
+{
+ // calculate grid layout
+
+  triclinic = domain->triclinic;
+
+  if (triclinic == 0) {
+    prd = domain->prd;
+    boxlo = domain->boxlo;
+    sublo = domain->sublo;
+    subhi = domain->subhi;
+  } else {
+    prd = domain->prd_lamda;
+    boxlo = domain->boxlo_lamda;
+    sublo = domain->sublo_lamda;
+    subhi = domain->subhi_lamda;
+  }
+
+  double xprd = prd[0];
+  double yprd = prd[1];
+  double zprd = prd[2];
+
+  delxinv = nx/xprd;
+  delyinv = ny/yprd;
+  delzinv = nz/zprd;
+  
+  delx = 1.0/delxinv;
+  dely = 1.0/delyinv;
+  delz = 1.0/delzinv;
+}
+
+/* ----------------------------------------------------------------------
+   set local subset of grid that I own
+   n xyz lo/hi = 3d brick that I own (inclusive)
+------------------------------------------------------------------------- */
+
+void ComputeGrid::set_grid_local()
+{
+  // global indices of grid range from 0 to N-1
+  // nlo,nhi = lower/upper limits of the 3d sub-brick of
+  //   global grid that I own without ghost cells
+
+  nxlo = static_cast<int> (comm->xsplit[comm->myloc[0]] * nx);
+  nxhi = static_cast<int> (comm->xsplit[comm->myloc[0]+1] * nx) - 1;
+
+  nylo = static_cast<int> (comm->ysplit[comm->myloc[1]] * ny);
+  nyhi = static_cast<int> (comm->ysplit[comm->myloc[1]+1] * ny) - 1;
+
+  nzlo = static_cast<int> (comm->zsplit[comm->myloc[2]] * nz);
+  nzhi = static_cast<int> (comm->zsplit[comm->myloc[2]+1] * nz) - 1;
+
+  ngridlocal = (nxhi - nxlo + 1) * (nyhi - nylo + 1) * (nzhi - nzlo + 1);
+
+  printf("me = %d n = %d x %d %d y %d %d z %d %d \n", comm->me, ngridlocal, nxlo, nxhi, nylo, nyhi, nzlo, nzhi);
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local data
+------------------------------------------------------------------------- */
+
+double ComputeGrid::memory_usage()
+{
+  double nbytes = size_array_rows*size_array_cols * 
+    sizeof(double);                             // grid
+  nbytes += size_array_rows*sizeof(int); // local_flags
+  nbytes += size_array_cols*ngridlocal*sizeof(double); // gridlocal
+  return nbytes;
+}
diff --git a/src/compute_grid.h b/src/compute_grid.h
new file mode 100644
index 0000000000..1b2797732d
--- /dev/null
+++ b/src/compute_grid.h
@@ -0,0 +1,75 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, 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.
+------------------------------------------------------------------------- */
+
+#ifndef LMP_COMPUTE_GRID_H
+#define LMP_COMPUTE_GRID_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeGrid : public Compute {
+ public:
+
+  ComputeGrid(class LAMMPS *, int, char **);
+  virtual ~ComputeGrid();
+  void init();
+  void setup();
+  virtual void compute_array() = 0;
+
+  double memory_usage();
+
+ protected:
+  int nx, ny, nz;                      // global grid dimensions
+  int nxlo, nxhi, nylo, nyhi, nzlo, nzhi; // local grid bounds, inclusive
+  int ngrid;                           // number of global grid points
+  int ngridlocal;                      // number of local grid points
+  int nvalues;                         // number of values per grid point
+  double **grid;                       // global grid
+  double **gridall;                    // global grid summed over procs
+  double ****gridlocal;                // local grid
+  int triclinic;                       // triclinic flag
+  double *boxlo, *prd;                 // box info (units real/ortho or reduced/tri)
+  double *sublo, *subhi;               // subdomain info (units real/ortho or reduced/tri)
+  double delxinv,delyinv,delzinv;      // inverse grid spacing
+  double delx,dely,delz;               // grid spacing
+  int nargbase;                        // number of base class args 
+  double cutmax;                       // largest cutoff distance
+  int size_array_cols_base;            // number of columns used for coords, etc.
+  int *local_flags;                    // local flag for each grid point
+  int gridlocal_allocated;             // shows if gridlocal allocated
+
+  void allocate();
+  void grid2x(int, double*);           // convert grid point to coord
+  void grid2ix(int, int&, int&, int&); // convert grid point to ix, iy, iz
+  void assign_coords();                // assign coords for grid
+  void assign_local_flags();           // set local flag for each grid point
+  int check_local(int);                // check if grid point is local
+  void set_grid_global();              // set global grid
+  void set_grid_local();               // set bounds for local grid
+ private:
+};
+
+}
+
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+*/