Created in.grid.local, fixed some problems in src

examples/snap/in.grid.local

@@ -0,0 +1,73 @@
+# Demonstrate bispectrum computes with local grid
+
+# 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 mygridlocal all sna/grid/local 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
+thermo_modify norm yes
+
+dump mydump_b all custom 1000 dump_b id c_b[*]
+dump mydump_bgridlocal all local 1000 dump_bgridlocal index c_mygridlocal[*]
+
+# run
+
+run 0
+

src/ML-SNAP/compute_sna_grid_local.cpp

@@ -0,0 +1,308 @@
+/* ----------------------------------------------------------------------
+   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_local.h"
+#include "compute_sna_grid_local.h"
+#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 "domain.h"
+#include "comm.h"
+#include "memory.h"
+#include "error.h"
+#include "tokenizer.h"
+
+#include <cmath>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+
+ComputeSNAGridLocal::ComputeSNAGridLocal(LAMMPS *lmp, int narg, char **arg) :
+  ComputeGridLocal(lmp, narg, arg), cutsq(nullptr),
+  radelem(nullptr), wjelem(nullptr)
+{
+  double rfac0, rmin0;
+  int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag;
+
+  // 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/local command");
+
+  // default values
+
+  rmin0 = 0.0;
+  switchflag = 1;
+  bzeroflag = 1;
+  quadraticflag = 0;
+  chemflag = 0;
+  bnormflag = 0;
+  wselfallflag = 0;
+  nelements = 1;
+  
+  // process required arguments
+
+  memory->create(radelem,ntypes+1,"sna/grid/local:radelem"); // offset by 1 to match up with types
+  memory->create(wjelem,ntypes+1,"sna/grid/local: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/local: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/local 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/local 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/local 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/local command");
+      quadraticflag = atoi(arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"chem") == 0) {
+      if (iarg+2 > narg)
+        error->all(FLERR,"Illegal compute sna/grid/local command");
+      chemflag = 1;
+      memory->create(map,ntypes+1,"compute_sna_grid_local: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 sna/grid/local command");
+        map[i+1] = jelem;
+      }
+      iarg += 2+ntypes;
+    } else if (strcmp(arg[iarg],"bnormflag") == 0) {
+      if (iarg+2 > narg)
+        error->all(FLERR,"Illegal compute sna/grid/local command");
+      bnormflag = atoi(arg[iarg+1]);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"wselfallflag") == 0) {
+      if (iarg+2 > narg)
+        error->all(FLERR,"Illegal compute sna/grid/local command");
+      wselfallflag = atoi(arg[iarg+1]);
+      iarg += 2;
+    } else error->all(FLERR,"Illegal compute sna/grid/local command");
+
+  }
+
+  snaptr = new SNA(lmp, rfac0, twojmax,
+                   rmin0, switchflag, bzeroflag,
+                   chemflag, bnormflag, wselfallflag, nelements);
+
+  ncoeff = snaptr->ncoeff;
+  nvalues = ncoeff;
+  if (quadraticflag) nvalues += (ncoeff*(ncoeff+1))/2;
+  size_local_cols = size_local_cols_base + nvalues;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeSNAGridLocal::~ComputeSNAGridLocal()
+{
+  memory->destroy(radelem);
+  memory->destroy(wjelem);
+  memory->destroy(cutsq);
+  delete snaptr;
+
+  if (chemflag) memory->destroy(map);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeSNAGridLocal::init()
+{
+  if (force->pair == nullptr)
+    error->all(FLERR,"Compute sna/grid/local requires a pair style be defined");
+
+  if (cutmax > force->pair->cutforce)
+    error->all(FLERR,"Compute sna/grid/local 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/local") == 0) count++;
+  if (count > 1 && comm->me == 0)
+    error->warning(FLERR,"More than one compute sna/grid/local");
+  snaptr->init();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeSNAGridLocal::init_list(int /*id*/, NeighList *ptr)
+{
+  list = ptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeSNAGridLocal::compute_local()
+{
+  invoked_array = update->ntimestep;
+
+  // compute sna for each gridpoint
+
+  double** const x = atom->x;
+  const int* const mask = atom->mask;
+  int * const type = atom->type;
+  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++) {
+	double xgrid[3];
+	grid2x(ix, iy, iz, xgrid);
+	const double xtmp = xgrid[0];
+	const double ytmp = xgrid[1];
+	const double ztmp = xgrid[2];
+
+	// currently, all grid points are type 1
+	
+	const int itype = 1;
+	int ielem = 0;
+	if (chemflag)
+	  ielem = map[itype];
+	const double radi = radelem[itype];
+
+	// 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];
+	  int jelem = 0;
+	  if (chemflag)
+	    jelem = map[jtype];
+	  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;
+	    snaptr->element[ninside] = jelem; // element index for multi-element snap
+	    ninside++;
+	  }
+	}
+
+	snaptr->compute_ui(ninside, ielem);
+	snaptr->compute_zi();
+	snaptr->compute_bi(ielem);
+
+	// linear contributions
+
+	for (int icoeff = 0; icoeff < ncoeff; icoeff++)
+	  gridlocal[size_local_cols_base+icoeff][iz][iy][ix] = snaptr->blist[icoeff];
+	
+	// quadratic contributions
+
+	if (quadraticflag) {
+	  int ncount = ncoeff;
+	  for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+	    double bveci = snaptr->blist[icoeff];
+	    gridlocal[size_local_cols_base+ncount++][iz][iy][ix] = 0.5*bveci*bveci;
+	    for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++)
+	      gridlocal[size_local_cols_base+ncount++][iz][iy][ix] = bveci*snaptr->blist[jcoeff];
+	  }
+	}
+      }
+
+  // copy 4d array to 2d array
+
+  copy_gridlocal_to_local_array();
+}
+
+
+/* ----------------------------------------------------------------------
+   memory usage
+------------------------------------------------------------------------- */
+
+double ComputeSNAGridLocal::memory_usage()
+{
+  double nbytes = snaptr->memory_usage();    // SNA object
+  int n = atom->ntypes+1;
+  nbytes += (double)n*sizeof(int);            // map
+
+  return nbytes;
+}
+

src/ML-SNAP/compute_sna_grid_local.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.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(sna/grid/local,ComputeSNAGridLocal)
+
+#else
+
+#ifndef LMP_COMPUTE_SNA_GRID_LOCAL_H
+#define LMP_COMPUTE_SNA_GRID_LOCAL_H
+
+#include "compute_grid_local.h"
+
+namespace LAMMPS_NS {
+
+class ComputeSNAGridLocal : public ComputeGridLocal {
+ public:
+  ComputeSNAGridLocal(class LAMMPS *, int, char **);
+  ~ComputeSNAGridLocal();
+  void init();
+  void init_list(int, class NeighList *);
+  void compute_local();
+  double memory_usage();
+
+ private:
+  int ncoeff;
+  double **cutsq;
+  class NeighList *list;
+  double rcutfac;
+  double *radelem;
+  double *wjelem;
+  int *map;    // map types to [0,nelements)
+  int nelements, chemflag;
+  class SNA *snaptr;
+  double cutmax;
+  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/local requires a pair style be defined
+
+Self-explanatory.
+
+E: Compute sna/grid/local cutoff is longer than pairwise cutoff
+
+Self-explanatory.
+
+W: More than one compute sna/grid/local
+
+Self-explanatory.
+
+*/

src/compute_grid_local.cpp

@@ -0,0 +1,253 @@
+/* ----------------------------------------------------------------------
+   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_local.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;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeGridLocal::ComputeGridLocal(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg), gridlocal(nullptr), alocal(nullptr)
+{
+  if (narg < 6) error->all(FLERR,"Illegal compute grid/local command");
+
+  local_flag = 1;
+  size_local_cols = 0;
+  size_local_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/local command");
+    nx = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
+    ny = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
+    nz = utils::inumeric(FLERR,arg[iarg+3],false,lmp);
+    if (nx <= 0 || ny <= 0 || nz <= 0)
+      error->all(FLERR,"All grid/local dimensions must be positive");
+    iarg += 4;
+  } else error->all(FLERR,"Illegal compute grid/local command");
+
+  nargbase = iarg - iarg0;
+
+  size_local_cols_base = 6;
+  gridlocal_allocated = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeGridLocal::~ComputeGridLocal()
+{
+  if (gridlocal_allocated) {
+    gridlocal_allocated = 0;
+    memory->destroy4d_offset(gridlocal,nzlo,nylo,nxlo);
+  }
+  memory->destroy(alocal);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeGridLocal::init()
+{
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeGridLocal::setup()
+{
+  set_grid_global();
+  set_grid_local();
+  allocate();
+}
+
+/* ----------------------------------------------------------------------
+   convert global array indexes to box coords
+------------------------------------------------------------------------- */
+
+void ComputeGridLocal::grid2x(int ix, int iy, int iz, double *x)
+{
+  x[0] = ix*delx;
+  x[1] = iy*dely;
+  x[2] = iz*delz;
+
+  if (triclinic) domain->lamda2x(x, x);
+}
+
+/* ----------------------------------------------------------------------
+   free and reallocate arrays
+------------------------------------------------------------------------- */
+
+void ComputeGridLocal::allocate()
+{
+  // allocate local array
+
+  if (gridlocal_allocated) {
+    gridlocal_allocated = 0;
+    // MEMORY LEAK!!
+    // can't seem to free this memory without seg-fault
+    // printf("Before allocate destroy4d, proc %d %p\n",comm->me,gridlocal);
+    // memory->destroy4d_offset(gridlocal,nzlo,nylo,nxlo);
+    // printf("After allocate destroy4d, proc %d %p\n",comm->me,gridlocal);
+  }
+
+  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");
+  }
+}
+
+/* ----------------------------------------------------------------------
+   set global grid
+------------------------------------------------------------------------- */
+
+void ComputeGridLocal::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 ComputeGridLocal::set_grid_local()
+{
+  // nx,ny,nz = extent of global grid
+  // indices into the global grid range from 0 to N-1 in each dim
+  // if grid point is inside my sub-domain I own it,
+  //   this includes sub-domain lo boundary but excludes hi boundary
+  // ixyz lo/hi = inclusive lo/hi bounds of global grid sub-brick I own
+  // if proc owns no grid cells in a dim, then ilo > ihi
+  // if 2 procs share a boundary a grid point is exactly on,
+  //   the 2 equality if tests insure a consistent decision
+  //   as to which proc owns it
+
+  double xfraclo,xfrachi,yfraclo,yfrachi,zfraclo,zfrachi;
+
+  if (comm->layout != Comm::LAYOUT_TILED) {
+    xfraclo = comm->xsplit[comm->myloc[0]];
+    xfrachi = comm->xsplit[comm->myloc[0]+1];
+    yfraclo = comm->ysplit[comm->myloc[1]];
+    yfrachi = comm->ysplit[comm->myloc[1]+1];
+    zfraclo = comm->zsplit[comm->myloc[2]];
+    zfrachi = comm->zsplit[comm->myloc[2]+1];
+  } else {
+    xfraclo = comm->mysplit[0][0];
+    xfrachi = comm->mysplit[0][1];
+    yfraclo = comm->mysplit[1][0];
+    yfrachi = comm->mysplit[1][1];
+    zfraclo = comm->mysplit[2][0];
+    zfrachi = comm->mysplit[2][1];
+  }
+
+  nxlo = static_cast<int> (xfraclo * nx);
+  if (1.0*nxlo != xfraclo*nx) nxlo++;
+  nxhi = static_cast<int> (xfrachi * nx);
+  if (1.0*nxhi == xfrachi*nx) nxhi--;
+
+  nylo = static_cast<int> (yfraclo * ny);
+  if (1.0*nylo != yfraclo*ny) nylo++;
+  nyhi = static_cast<int> (yfrachi * ny);
+  if (1.0*nyhi == yfrachi*ny) nyhi--;
+
+  nzlo = static_cast<int> (zfraclo * nz);
+  if (1.0*nzlo != zfraclo*nz) nzlo++;
+  nzhi = static_cast<int> (zfrachi * nz);
+  if (1.0*nzhi == zfrachi*nz) nzhi--;
+
+  ngridlocal = (nxhi - nxlo + 1) * (nyhi - nylo + 1) * (nzhi - nzlo + 1);
+  size_local_rows = ngridlocal;
+
+  memory->destroy(alocal);
+  memory->create(alocal, size_local_rows, size_local_cols, "compute/grid/local:alocal");
+  array_local = alocal;
+
+  int igrid = 0;
+  for (int iz = nzlo; iz <= nzhi; iz++)
+    for (int iy = nylo; iy <= nyhi; iy++)
+      for (int ix = nxlo; ix <= nxhi; ix++) {
+	alocal[igrid][0] = ix;
+	alocal[igrid][1] = iy;
+	alocal[igrid][2] = iz;
+	double xgrid[3];
+	grid2x(ix, iy, iz, xgrid);
+	alocal[igrid][3] = xgrid[0];
+	alocal[igrid][4] = xgrid[1];
+	alocal[igrid][5] = xgrid[2];
+	igrid++;
+      }
+}
+
+/* ----------------------------------------------------------------------
+   copy the 4d gridlocal array values to the 2d local array
+------------------------------------------------------------------------- */
+
+void ComputeGridLocal::copy_gridlocal_to_local_array()
+{
+  int igrid = 0;
+  for (int iz = nzlo; iz <= nzhi; iz++)
+    for (int iy = nylo; iy <= nyhi; iy++)
+      for (int ix = nxlo; ix <= nxhi; ix++) {
+	for (int icol = size_local_cols_base; icol < size_local_cols; icol++)
+	  alocal[igrid][icol] = gridlocal[icol][iz][iy][ix];
+	igrid++;
+      }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local data
+------------------------------------------------------------------------- */
+
+double ComputeGridLocal::memory_usage()
+{
+  int nbytes = size_array_cols*ngridlocal*sizeof(double); // gridlocal
+  return nbytes;
+}

src/compute_grid_local.h

@@ -0,0 +1,69 @@
+/* -*- 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_LOCAL_H
+#define LMP_COMPUTE_GRID_LOCAL_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeGridLocal : public Compute {
+ public:
+
+  ComputeGridLocal(class LAMMPS *, int, char **);
+  virtual ~ComputeGridLocal();
+  void init();
+  void setup();
+  virtual void compute_local() = 0;
+
+  double memory_usage();
+
+ protected:
+  int nx, ny, nz;                      // global grid dimensions
+  int nxlo, nxhi, nylo, nyhi, nzlo, nzhi; // local grid bounds, inclusive
+  int ngridlocal;                      // number of local grid points
+  int nvalues;                         // number of values per grid point
+  double ****gridlocal;                // local grid
+  double **alocal;                     // pointer to Compute::array_local
+  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_local_cols_base;            // number of columns used for coords, etc.
+  int gridlocal_allocated;             // shows if gridlocal allocated
+
+  void allocate();
+  void grid2x(int, int, int, double*); // convert global indices to coordinates
+  void set_grid_global();              // set global grid
+  void set_grid_local();               // set bounds for local grid
+  void copy_gridlocal_to_local_array();// copy 4d gridlocal array to 2d local array
+ 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.
+
+*/