diff --git a/src/EXTRA-FIX/fix_ttm.cpp b/src/EXTRA-FIX/fix_ttm.cpp
index eaf095d329..798ab97dca 100644
--- a/src/EXTRA-FIX/fix_ttm.cpp
+++ b/src/EXTRA-FIX/fix_ttm.cpp
@@ -37,8 +37,11 @@ using namespace LAMMPS_NS;
 using namespace FixConst;
 
 #define MAXLINE 1024
-#define OFFSET 16384    // to avoid outside-of-box atoms being rounded incorrectly
-#define SHIFT 0.5       // 0.5 for nearest grid point, 0.0 for lower-left grid point
+//#define OFFSET 16384    // to avoid outside-of-box atoms being rounded incorrectly
+//#define SHIFT 0.5       // 0.5 for nearest grid point, 0.0 for lower-left grid point
+
+#define OFFSET 0    // to avoid outside-of-box atoms being rounded incorrectly
+#define SHIFT 0.0       // 0.5 for nearest grid point, 0.0 for lower-left grid point
 
 /* ---------------------------------------------------------------------- */
 
@@ -114,7 +117,6 @@ FixTTM::FixTTM(LAMMPS *lmp, int narg, char **arg) :
   ngridtotal = totalgrid;
 
   // allocate per-atom flangevin and zero it
-  // NOTE: is init to zero necessary?
 
   flangevin = nullptr;
   grow_arrays(atom->nmax);
@@ -130,33 +132,13 @@ FixTTM::FixTTM(LAMMPS *lmp, int narg, char **arg) :
   atom->add_callback(Atom::GROW);
   atom->add_callback(Atom::RESTART);
 
+  // determines which class deallocate_grid() is called from
+
   deallocate_flag = 0;
 }
 
 /* ---------------------------------------------------------------------- */
 
-void FixTTM::post_constructor()
-{
-  // allocate 3d grid variables
-
-  allocate_grid();
-
-  // zero net_energy_transfer
-  // in case compute_vector accesses it on timestep 0
-
-  outflag = 0;
-  for (int ix = 0; ix < nxgrid; ix++)
-    for (int iy = 0; iy < nygrid; iy++)
-      for (int iz = 0; iz < nzgrid; iz++)
-        net_energy_transfer_all[ix][iy][iz] = 0;
-
-  // set initial electron temperatures from user input file
-
-  read_electron_temperatures(infile);
-}
-
-/* ---------------------------------------------------------------------- */
-
 FixTTM::~FixTTM()
 {
   delete [] infile;
@@ -173,6 +155,26 @@ FixTTM::~FixTTM()
 
 /* ---------------------------------------------------------------------- */
 
+void FixTTM::post_constructor()
+{
+  // allocate global grid on each proc
+  // needs to be done in post_contructor() beccause is virtual method
+
+  allocate_grid();
+
+  // zero net_energy_transfer
+  // in case compute_vector accesses it on timestep 0
+
+  outflag = 0;
+  memset(&net_energy_transfer[0][0][0],0,ngridtotal*sizeof(double));
+
+  // set initial electron temperatures from user input file
+
+  read_electron_temperatures(infile);
+}
+
+/* ---------------------------------------------------------------------- */
+
 int FixTTM::setmask()
 {
   int mask = 0;
@@ -269,10 +271,10 @@ void FixTTM::post_force(int /*vflag*/)
       if (iy >= nygrid) iy -= nygrid;
       if (iz >= nzgrid) iz -= nzgrid;
 
-      if (T_electron[ix][iy][iz] < 0)
+      if (T_electron[iz][iy][ix] < 0)
         error->all(FLERR,"Electronic temperature dropped below zero");
 
-      double tsqrt = sqrt(T_electron[ix][iy][iz]);
+      double tsqrt = sqrt(T_electron[iz][iy][ix]);
 
       gamma1 = gfactor1[type[i]];
       double vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
@@ -320,10 +322,10 @@ void FixTTM::end_of_step()
   int nlocal = atom->nlocal;
 
   outflag = 0;
-  for (ix = 0; ix < nxgrid; ix++)
+  for (iz = 0; iz < nzgrid; iz++)
     for (iy = 0; iy < nygrid; iy++)
-      for (iz = 0; iz < nzgrid; iz++)
-        net_energy_transfer[ix][iy][iz] = 0;
+      for (ix = 0; ix < nxgrid; ix++)
+        net_energy_transfer[iz][iy][ix] = 0;
 
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
@@ -339,7 +341,7 @@ void FixTTM::end_of_step()
       if (ix >= nxgrid) ix -= nxgrid;
       if (iy >= nygrid) iy -= nygrid;
       if (iz >= nzgrid) iz -= nzgrid;
-      net_energy_transfer[ix][iy][iz] +=
+      net_energy_transfer[iz][iy][ix] +=
         (flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
          flangevin[i][2]*v[i][2]);
     }
@@ -375,18 +377,18 @@ void FixTTM::end_of_step()
   // finite difference iterations to update T_electron
 
   for (int istep = 0; istep < num_inner_timesteps; istep++) {
-
-    for (ix = 0; ix < nxgrid; ix++)
+    
+    for (iz = 0; iz < nzgrid; iz++)
       for (iy = 0; iy < nygrid; iy++)
-        for (iz = 0; iz < nzgrid; iz++)
-          T_electron_old[ix][iy][iz] =
-            T_electron[ix][iy][iz];
+        for (ix = 0; ix < nxgrid; ix++)
+          T_electron_old[iz][iy][ix] =
+            T_electron[iz][iy][ix];
 
     // compute new electron T profile
 
-    for (ix = 0; ix < nxgrid; ix++)
+    for (iz = 0; iz < nzgrid; iz++)
       for (iy = 0; iy < nygrid; iy++)
-        for (iz = 0; iz < nzgrid; iz++) {
+        for (ix = 0; ix < nxgrid; ix++) {
           int right_xnode = ix + 1;
           int right_ynode = iy + 1;
           int right_znode = iz + 1;
@@ -400,22 +402,22 @@ void FixTTM::end_of_step()
           if (left_ynode == -1) left_ynode = nygrid - 1;
           if (left_znode == -1) left_znode = nzgrid - 1;
 
-          T_electron[ix][iy][iz] =
-            T_electron_old[ix][iy][iz] +
+          T_electron[iz][iy][ix] =
+            T_electron_old[iz][iy][ix] +
             inner_dt/(electronic_specific_heat*electronic_density) *
             (electronic_thermal_conductivity *
 
-             ((T_electron_old[right_xnode][iy][iz] +
-               T_electron_old[left_xnode][iy][iz] -
-               2*T_electron_old[ix][iy][iz])/dx/dx +
-              (T_electron_old[ix][right_ynode][iz] +
-               T_electron_old[ix][left_ynode][iz] -
-               2*T_electron_old[ix][iy][iz])/dy/dy +
-              (T_electron_old[ix][iy][right_znode] +
-               T_electron_old[ix][iy][left_znode] -
-               2*T_electron_old[ix][iy][iz])/dz/dz) -
-
-             (net_energy_transfer_all[ix][iy][iz])/del_vol);
+             ((T_electron_old[iz][iy][right_xnode] +
+               T_electron_old[iz][iy][left_xnode] -
+               2*T_electron_old[iz][iy][ix])/dx/dx +
+              (T_electron_old[iz][right_ynode][ix] +
+               T_electron_old[iz][left_ynode][ix] -
+               2*T_electron_old[iz][iy][ix])/dy/dy +
+              (T_electron_old[right_znode][iy][iz] +
+               T_electron_old[left_znode][iy][ix] -
+               2*T_electron_old[iz][iy][ix])/dz/dz) -
+             
+             (net_energy_transfer_all[iz][iy][ix])/del_vol);
         }
   }
 }
@@ -462,18 +464,18 @@ void FixTTM::read_electron_temperatures(const char *filename)
     if (T_tmp < 0.0)
       error->one(FLERR,"Fix ttm electron temperatures must be > 0.0");
 
-    T_electron[ix][iy][iz] = T_tmp;
-    T_initial_set[ix][iy][iz] = 1;
+    T_electron[iz][iy][ix] = T_tmp;
+    T_initial_set[iz][iy][ix] = 1;
   }
 
   fclose(fp);
 
   // check completeness of input data
 
-  for (int ix = 0; ix < nxgrid; ix++)
+  for (int iz = 0; iz < nzgrid; iz++)
     for (int iy = 0; iy < nygrid; iy++)
-      for (int iz = 0; iz < nzgrid; iz++)
-        if (T_initial_set[ix][iy][iz] == 0)
+      for (int ix = 0; ix < nxgrid; ix++)
+        if (T_initial_set[iz][iy][ix] == 0)
           error->one(FLERR,"Initial temperatures not all set in fix ttm");
 
   memory->destroy(T_initial_set);
@@ -494,7 +496,7 @@ void FixTTM::reset_dt()
 
 void FixTTM::grow_arrays(int ngrow)
 {
-  memory->grow(flangevin,ngrow,3,"TTM:flangevin");
+  memory->grow(flangevin,ngrow,3,"ttm:flangevin");
 }
 
 /* ----------------------------------------------------------------------
@@ -504,15 +506,20 @@ void FixTTM::grow_arrays(int ngrow)
 void FixTTM::write_restart(FILE *fp)
 {
   double *rlist;
-  memory->create(rlist,nxgrid*nygrid*nzgrid+1,"TTM:rlist");
+  memory->create(rlist,nxgrid*nygrid*nzgrid+4,"ttm:rlist");
 
   int n = 0;
+  rlist[n++] = nxgrid;
+  rlist[n++] = nygrid;
+  rlist[n++] = nzgrid;
   rlist[n++] = seed;
 
-  for (int ix = 0; ix < nxgrid; ix++)
+  // store global grid values
+
+  for (int iz = 0; iz < nzgrid; iz++)
     for (int iy = 0; iy < nygrid; iy++)
-      for (int iz = 0; iz < nzgrid; iz++)
-        rlist[n++] =  T_electron[ix][iy][iz];
+      for (int ix = 0; ix < nxgrid; ix++)
+        rlist[n++] =  T_electron[iz][iy][ix];
 
   if (comm->me == 0) {
     int size = n * sizeof(double);
@@ -532,18 +539,28 @@ void FixTTM::restart(char *buf)
   int n = 0;
   double *rlist = (double *) buf;
 
+  // check that restart grid size is same as current grid size
+
+  int nxgrid_old = static_cast<int> (rlist[n++]);
+  int nygrid_old = static_cast<int> (rlist[n++]);
+  int nzgrid_old = static_cast<int> (rlist[n++]);
+
+  if (nxgrid_old != nxgrid || nygrid_old != nygrid || nzgrid_old != nzgrid)
+    error->all(FLERR,"Must restart fix ttm with same grid size");
+
   // change RN seed from initial seed, to avoid same Langevin factors
-  // just increment by 1, since for RanMars that is new RN streamd
+  // just increment by 1, since for RanMars that is a new RN stream
 
   seed = static_cast<int> (rlist[n++]) + 1;
-
-  for (int ix = 0; ix < nxgrid; ix++)
-    for (int iy = 0; iy < nygrid; iy++)
-      for (int iz = 0; iz < nzgrid; iz++)
-        T_electron[ix][iy][iz] = rlist[n++];
-
   delete random;
   random = new RanMars(lmp,seed+comm->me);
+
+  // restore global grid values
+
+  for (int iz = 0; iz < nzgrid; iz++)
+    for (int iy = 0; iy < nygrid; iy++)
+      for (int ix = 0; ix < nxgrid; ix++)
+        T_electron[iz][iy][ix] = rlist[n++];
 }
 
 /* ----------------------------------------------------------------------
@@ -617,14 +634,14 @@ double FixTTM::compute_vector(int n)
     double dz = domain->zprd/nzgrid;
     double del_vol = dx*dy*dz;
 
-    for (ix = 0; ix < nxgrid; ix++)
+    for (iz = 0; iz < nzgrid; iz++)
       for (iy = 0; iy < nygrid; iy++)
-        for (iz = 0; iz < nzgrid; iz++) {
+        for (ix = 0; ix < nxgrid; ix++) {
           e_energy +=
-            T_electron[ix][iy][iz]*electronic_specific_heat*
+            T_electron[iz][iy][ix]*electronic_specific_heat*
             electronic_density*del_vol;
           transfer_energy +=
-            net_energy_transfer_all[ix][iy][iz]*update->dt;
+            net_energy_transfer_all[iz][iy][ix]*update->dt;
         }
 
     outflag = 1;
@@ -653,11 +670,11 @@ double FixTTM::memory_usage()
 
 void FixTTM::allocate_grid()
 {
-  memory->create(T_electron_old,nxgrid,nygrid,nzgrid,"ttm:T_electron_old");
-  memory->create(T_electron,nxgrid,nygrid,nzgrid,"ttm:T_electron");
-  memory->create(net_energy_transfer,nxgrid,nygrid,nzgrid,
+  memory->create(T_electron_old,nzgrid,nygrid,nxgrid,"ttm:T_electron_old");
+  memory->create(T_electron,nzgrid,nygrid,nxgrid,"ttm:T_electron");
+  memory->create(net_energy_transfer,nzgrid,nygrid,nxgrid,
                  "ttm:net_energy_transfer");
-  memory->create(net_energy_transfer_all,nxgrid,nygrid,nzgrid,
+  memory->create(net_energy_transfer_all,nzgrid,nygrid,nxgrid,
                  "ttm:net_energy_transfer_all");
 }
 
diff --git a/src/EXTRA-FIX/fix_ttm.h b/src/EXTRA-FIX/fix_ttm.h
index a7e7e0ea45..0b8900977e 100644
--- a/src/EXTRA-FIX/fix_ttm.h
+++ b/src/EXTRA-FIX/fix_ttm.h
@@ -30,7 +30,7 @@ class FixTTM : public Fix {
   virtual ~FixTTM();
   virtual void post_constructor();
   int setmask();
-  void init();
+  virtual void init();
   void setup(int);
   void post_force_setup(int);
   virtual void post_force(int);
diff --git a/src/EXTRA-FIX/fix_ttm_grid.cpp b/src/EXTRA-FIX/fix_ttm_grid.cpp
index 3ca94e10b4..b9e52b1d68 100644
--- a/src/EXTRA-FIX/fix_ttm_grid.cpp
+++ b/src/EXTRA-FIX/fix_ttm_grid.cpp
@@ -43,7 +43,10 @@ static constexpr int CHUNK = 1024;
 /* ---------------------------------------------------------------------- */
 
 FixTTMGrid::FixTTMGrid(LAMMPS *lmp, int narg, char **arg) : 
-  FixTTM(lmp, narg, arg) {}
+  FixTTM(lmp, narg, arg)
+{
+  skin_original = neighbor->skin;
+}
 
 /* ---------------------------------------------------------------------- */
 
@@ -57,13 +60,15 @@ FixTTMGrid::~FixTTMGrid()
 
 void FixTTMGrid::post_constructor()
 {
-  // allocate 3d grid variables
+  // allocate global grid on each proc
+  // needs to be done in post_contructor() beccause is virtual method
 
   allocate_grid();
 
   // zero net_energy_transfer
   // in case compute_vector accesses it on timestep 0
 
+  outflag = 0;
   memset(&net_energy_transfer[nzlo_out][nylo_out][nxlo_out],0,
          ngridout*sizeof(double));
 
@@ -74,6 +79,16 @@ void FixTTMGrid::post_constructor()
 
 /* ---------------------------------------------------------------------- */
 
+void FixTTMGrid::init()
+{
+  FixTTM::init();
+
+  if (neighbor->skin > skin_original)
+    error->all(FLERR,"Cannot extend neighbor skin after fix ttm/griddefined");
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixTTMGrid::post_force(int /*vflag*/)
 {
   int ix,iy,iz;
@@ -232,10 +247,6 @@ void FixTTMGrid::end_of_step()
 
   gc->forward_comm(GridComm::FIX,this,1,sizeof(double),0,
                    gc_buf1,gc_buf2,MPI_DOUBLE);
-
-  // assign electron temperature to each atom for fix output
-
-
 }
 
 /* ----------------------------------------------------------------------
@@ -411,8 +422,6 @@ void FixTTMGrid::allocate_grid()
   //   global grid that I own without ghost cells
   // both non-tiled and tiled proc layouts use 0-1 fractional subdomain info
 
-  // NOTE: replace with comm->partition_grid()
-
   if (comm->layout != Comm::LAYOUT_TILED) {
     nxlo_in = static_cast<int> (comm->xsplit[comm->myloc[0]] * nxgrid);
     nxhi_in = static_cast<int> (comm->xsplit[comm->myloc[0]+1] * nxgrid) - 1;
@@ -481,11 +490,11 @@ void FixTTMGrid::allocate_grid()
   memory->create(gc_buf2,ngc_buf2,"ttm/grid:gc_buf2");
 
   memory->create3d_offset(T_electron_old,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                          nxlo_out,nxhi_out,"ttm:T_electron_old");
+                          nxlo_out,nxhi_out,"ttm/grid:T_electron_old");
   memory->create3d_offset(T_electron,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                          nxlo_out,nxhi_out,"ttm:T_electron");
+                          nxlo_out,nxhi_out,"ttm/grid:T_electron");
   memory->create3d_offset(net_energy_transfer,nzlo_out,nzhi_out,nylo_out,nyhi_out,
-                          nxlo_out,nxhi_out,"ttm:net_energy_transfer");
+                          nxlo_out,nxhi_out,"ttm/grid:net_energy_transfer");
 }
 
 /* ----------------------------------------------------------------------
@@ -503,32 +512,6 @@ void FixTTMGrid::deallocate_grid()
   memory->destroy3d_offset(net_energy_transfer,nzlo_out,nylo_out,nxlo_out);
 }
 
-/* ----------------------------------------------------------------------
-   use state info from restart file to restart the Fix
-------------------------------------------------------------------------- */
-
-void FixTTMGrid::restart(char *buf)
-{
-  int ix,iy,iz;
-
-  int n = 0;
-  double *rlist = (double *) buf;
-
-  // the seed must be changed from the initial seed
-  // NOTE: 0.5 is whacky, could go to zero
-  // NOTE: maybe GridComm should provide a method to pack a grid with bounds
-
-  seed = static_cast<int> (0.5*rlist[n++]);
-
-  for (iz = nzlo_in; iz <= nzhi_in; iz++)
-    for (iy = nylo_in; iy <= nyhi_in; iy++)
-      for (ix = nxlo_in; ix <= nxhi_in; ix++)
-        T_electron[iz][iy][ix] = rlist[n++];
-
-  delete random;
-  random = new RanMars(lmp,seed+comm->me);
-}
-
 /* ----------------------------------------------------------------------
    pack entire state of Fix into one write
 ------------------------------------------------------------------------- */
@@ -538,12 +521,17 @@ void FixTTMGrid::write_restart(FILE *fp)
   int ix,iy,iz;
 
   double *rlist;
-  memory->create(rlist,nxgrid*nygrid*nzgrid+1,"TTM:rlist");
+  memory->create(rlist,nxgrid*nygrid*nzgrid+4,"ttm/grid:rlist");
 
   int n = 0;
+  rlist[n++] = nxgrid;
+  rlist[n++] = nygrid;
+  rlist[n++] = nzgrid;
   rlist[n++] = seed;
 
-  // NOTE: this is a parallel grid now, not a serial grid
+  // gather rest of rlist on proc 0 as global grid values
+
+  gc->gather(GridComm::FIX,this,1,sizeof(double),0,&rlist[4],MPI_DOUBLE);
 
   for (iz = nzlo_in; iz <= nzhi_in; iz++)
     for (iy = nylo_in; iy <= nyhi_in; iy++)
@@ -559,6 +547,86 @@ void FixTTMGrid::write_restart(FILE *fp)
   memory->destroy(rlist);
 }
 
+/* ----------------------------------------------------------------------
+   use state info from restart file to restart the Fix
+------------------------------------------------------------------------- */
+
+void FixTTMGrid::restart(char *buf)
+{
+  int ix,iy,iz;
+
+  int n = 0;
+  double *rlist = (double *) buf;
+
+  // check that restart grid size is same as current grid size
+
+  int nxgrid_old = static_cast<int> (rlist[n++]);
+  int nygrid_old = static_cast<int> (rlist[n++]);
+  int nzgrid_old = static_cast<int> (rlist[n++]);
+
+  if (nxgrid_old != nxgrid || nygrid_old != nygrid || nzgrid_old != nzgrid)
+    error->all(FLERR,"Must restart fix ttm/grid with same grid size");
+
+  // change RN seed from initial seed, to avoid same Langevin factors
+  // just increment by 1, since for RanMars that is a new RN stream
+
+  seed = static_cast<int> (rlist[n++]) + 1;
+  delete random;
+  random = new RanMars(lmp,seed+comm->me);
+
+  // extract this proc's local grid values from global grid in rlist
+
+  int iglobal;
+
+  for (iz = nzlo_in; iz <= nzhi_in; iz++)
+    for (iy = nylo_in; iy <= nyhi_in; iy++)
+      for (ix = nxlo_in; ix <= nxhi_in; ix++) {
+        iglobal = nygrid*nxgrid*iz + nxgrid*iy + ix;
+        T_electron[iz][iy][ix] = rlist[n+iglobal];
+      }
+}
+
+/* ----------------------------------------------------------------------
+   pack values from local grid into buf
+------------------------------------------------------------------------- */
+
+void FixTTMGrid::pack_gather_grid(int which, void *vbuf)
+{
+  int ix,iy,iz;
+
+  double *buf = (double *) vbuf;
+
+  int m = 0;
+  for (iz = nzlo_in; iz <= nzhi_in; iz++)
+    for (iy = nylo_in; iy <= nyhi_in; iy++)
+      for (ix = nxlo_in; ix <= nxhi_in; ix++)
+        buf[m++] = T_electron[iz][iy][ix];
+}
+
+/* ----------------------------------------------------------------------
+   unpack values from buf into global gbuf based on their indices
+------------------------------------------------------------------------- */
+
+void FixTTMGrid::unpack_gather_grid(int which, void *vbuf, void *vgbuf,
+                                    int xlo, int xhi, int ylo, int yhi, 
+                                    int zlo, int zhi)
+{
+  int ix,iy,iz;
+
+  double *buf = (double *) vbuf;
+  double *gbuf = (double *) vgbuf;
+
+  int iglobal;
+  int ilocal = 0;
+
+  for (iz = zlo; iz <= zhi; iz++)
+    for (iy = ylo; iy <= yhi; iy++)
+      for (ix = xlo; ix <= xhi; ix++) {
+        iglobal = nygrid*nxgrid*iz + nxgrid*iy + ix;
+        gbuf[iglobal] = buf[ilocal++];
+      }
+}
+
 /* ----------------------------------------------------------------------
    return the energy of the electronic subsystem 
    or the net_energy transfer between the subsystems
@@ -609,4 +677,3 @@ double FixTTMGrid::memory_usage()
   bytes += (double)3*ngridout * sizeof(double);
   return bytes;
 }
-
diff --git a/src/EXTRA-FIX/fix_ttm_grid.h b/src/EXTRA-FIX/fix_ttm_grid.h
index 617612d641..dec5ca871b 100644
--- a/src/EXTRA-FIX/fix_ttm_grid.h
+++ b/src/EXTRA-FIX/fix_ttm_grid.h
@@ -29,6 +29,7 @@ class FixTTMGrid : public FixTTM {
   FixTTMGrid(class LAMMPS *, int, char **);
   ~FixTTMGrid();
   void post_constructor();
+  void init();
   void post_force(int);
   void end_of_step();
 
@@ -38,6 +39,8 @@ class FixTTMGrid : public FixTTM {
   void unpack_forward_grid(int, void *, int, int *);
   void pack_reverse_grid(int, void *, int, int *);
   void unpack_reverse_grid(int, void *, int, int *);
+  void pack_gather_grid(int, void *);
+  void unpack_gather_grid(int, void *, void *, int, int, int, int, int, int);
 
   void write_restart(FILE *);
   void restart(char *);
@@ -49,6 +52,7 @@ class FixTTMGrid : public FixTTM {
   int nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in;
   int nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out;
   double delxinv,delyinv,delzinv;
+  double skin_original;
 
   class GridComm *gc;
   int ngc_buf1,ngc_buf2;
diff --git a/src/fix.h b/src/fix.h
index a48c4eeb6a..cacbf04a2a 100644
--- a/src/fix.h
+++ b/src/fix.h
@@ -209,6 +209,9 @@ class Fix : protected Pointers {
   virtual void unpack_forward_grid(int, void *, int, int *) {};
   virtual void pack_reverse_grid(int, void *, int, int *) {};
   virtual void unpack_reverse_grid(int, void *, int, int *) {};
+  virtual void pack_gather_grid(int, void *) {};
+  virtual void unpack_gather_grid(int, void *, void *, 
+                                  int, int, int, int, int, int) {};
 
   virtual double compute_scalar() { return 0.0; }
   virtual double compute_vector(int) { return 0.0; }
diff --git a/src/gridcomm.cpp b/src/gridcomm.cpp
index fc6d8296d8..098a410a54 100644
--- a/src/gridcomm.cpp
+++ b/src/gridcomm.cpp
@@ -1125,6 +1125,81 @@ reverse_comm_tiled(T *ptr, int nper, int nbyte, int which,
   }
 }
 
+/* ----------------------------------------------------------------------
+   gather global grid values to proc 0
+   proc 0 pings each proc for its contribution
+   pack/unpack operations are performed by caller via callbacks
+------------------------------------------------------------------------- */
+
+void GridComm::gather(int caller, void *ptr, int nper, int nbyte, int which,
+                      void *buf, MPI_Datatype datatype)
+{
+  int me = comm->me;
+  Fix *fptr = (Fix *) ptr;
+
+  // maxsize = max grid data owned by any proc
+
+  int mysize = (inxhi-inxlo+1) * (inyhi-inylo+1) * (inzhi-inzlo+1);
+  mysize *= nper;
+  int maxsize;
+  MPI_Allreduce(&mysize,&maxsize,1,MPI_INT,MPI_MAX,world);
+
+  // pack my data via callback to caller
+
+  char *mybuf;
+  if (me == 0) memory->create(mybuf,maxsize*nbyte,"GridComm:mybuf");
+  else memory->create(mybuf,mysize*nbyte,"GridComm:mybuf");
+  fptr->pack_gather_grid(which,mybuf);
+
+  // ping each proc for its data
+  // unpack into full buffer via callback to caller
+
+  int xlo,xhi,ylo,yhi,zlo,zhi,tmp;
+  int bounds[6];
+
+  if (me == 0) {
+    MPI_Status status;
+    MPI_Request request;
+
+    for (int iproc = 0; iproc < nprocs; iproc++) {
+      if (iproc) {
+        MPI_Irecv(mybuf,maxsize,datatype,iproc,0,world,&request);
+        MPI_Send(&tmp,0,MPI_INT,iproc,0,world);
+        MPI_Wait(&request,&status);
+        MPI_Recv(bounds,6,MPI_INT,iproc,0,world,&status);
+        xlo = bounds[0];
+        xhi = bounds[1];
+        ylo = bounds[2];
+        yhi = bounds[3];
+        zlo = bounds[4];
+        zhi = bounds[5];
+      } else {
+        xlo = inxlo;
+        xhi = inxhi;
+        ylo = inylo;
+        yhi = inyhi;
+        zlo = inzlo;
+        zhi = inzhi;
+      }
+
+      fptr->unpack_gather_grid(which,mybuf,buf,xlo,xhi,ylo,yhi,zlo,zhi);
+    }
+
+  } else {
+    MPI_Recv(&tmp,0,MPI_INT,0,0,world,MPI_STATUS_IGNORE);
+    MPI_Rsend(mybuf,mysize,datatype,0,0,world);
+    bounds[0] = inxlo;
+    bounds[1] = inxhi;
+    bounds[2] = inylo;
+    bounds[3] = inyhi;
+    bounds[4] = inzlo;
+    bounds[5] = inzhi;
+    MPI_Send(bounds,6,MPI_INT,0,0,world);
+  }
+
+  memory->destroy(mybuf);
+}
+
 /* ----------------------------------------------------------------------
    create swap stencil for grid own/ghost communication
    swaps covers all 3 dimensions and both directions
diff --git a/src/gridcomm.h b/src/gridcomm.h
index 86f2c30297..4b3c3de1b8 100644
--- a/src/gridcomm.h
+++ b/src/gridcomm.h
@@ -31,6 +31,7 @@ class GridComm : protected Pointers {
   int ghost_adjacent();
   void forward_comm(int, void *, int, int, int, void *, void *, MPI_Datatype);
   void reverse_comm(int, void *, int, int, int, void *, void *, MPI_Datatype);
+  void gather(int, void *, int, int, int, void *, MPI_Datatype);
 
  protected:
   int me, nprocs;
diff --git a/src/modify.cpp b/src/modify.cpp
index 995b3b82ac..06b3c6cdb9 100644
--- a/src/modify.cpp
+++ b/src/modify.cpp
@@ -915,6 +915,18 @@ Fix *Modify::add_fix(int narg, char **arg, int trysuffix)
   if (fix[ifix] == nullptr)
     error->all(FLERR,utils::check_packages_for_style("fix",arg[2],lmp));
 
+  // increment nfix (if new)
+
+  if (newflag) nfix++;
+
+  // post_constructor() can call virtual methods in parent or child
+  //   which would otherwise not yet be visible in child class
+  // post_constructor() allows new fix to create other fixes
+  // nfix increment must come first so recursive call to add_fix within
+  //   post_constructor() will see updated nfix
+
+  fix[ifix]->post_constructor();
+
   // check if Fix is in restart_global list
   // if yes, pass state info to the Fix so it can reset itself
 
@@ -946,15 +958,12 @@ Fix *Modify::add_fix(int narg, char **arg, int trysuffix)
                        fix[ifix]->style,fix[ifix]->id);
     }
 
-  // increment nfix (if new)
   // set fix mask values
-  // post_constructor() allows new fix to create other fixes
-  // nfix increment comes first so that recursive call to add_fix within
-  //   post_constructor() will see updated nfix
 
-  if (newflag) nfix++;
   fmask[ifix] = fix[ifix]->setmask();
-  fix[ifix]->post_constructor();
+
+  // return pointer to fix
+
   return fix[ifix];
 }
 
@@ -973,7 +982,6 @@ Fix *Modify::add_fix(const std::string &fixcmd, int trysuffix)
   return add_fix(args.size(),newarg.data(),trysuffix);
 }
 
-
 /* ----------------------------------------------------------------------
    replace replaceID fix with a new fix
    this is used by callers to preserve ordering of fixes