modify fix STORE/PERATOM and callers to allow ghost comm

src/AMOEBA/pair_amoeba.cpp

@@ -861,7 +861,7 @@ void PairAmoeba::init_style()
   Fix *myfix;
   if (first_flag) {
     id_pole = utils::strdup("AMOEBA_pole");
-    myfix = modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 13",id_pole,group->names[0]));
+    myfix = modify->add_fix(fmt::format("{} {} STORE/PERATOM 13 0 0 1",id_pole,group->names[0]));
     fixpole = dynamic_cast<FixStorePeratom *>(myfix);
   }
 
@@ -873,12 +873,12 @@ void PairAmoeba::init_style()
 
   if (first_flag && use_pred) {
     id_udalt = utils::strdup("AMOEBA_udalt");
-    myfix = modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 {} 3",
+    myfix = modify->add_fix(fmt::format("{} {} STORE/PERATOM {} 3 0 1",
                                         id_udalt, group->names[0], maxualt));
     fixudalt = dynamic_cast<FixStorePeratom *>(myfix);
 
     id_upalt = utils::strdup("AMOEBA_upalt");
-    myfix = modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 {} 3",
+    myfix = modify->add_fix(fmt::format("{} {} STORE/PERATOM {} 3 0 1",
                                         id_upalt, group->names[0], maxualt));
     fixupalt = dynamic_cast<FixStorePeratom *>(myfix);
   }

src/CORESHELL/compute_temp_cs.cpp

@@ -68,7 +68,7 @@ ComputeTempCS::ComputeTempCS(LAMMPS *lmp, int narg, char **arg) :
 
   id_fix = utils::strdup(id + std::string("_COMPUTE_STORE"));
   fix = dynamic_cast<FixStorePeratom *>(
-    modify->add_fix(fmt::format("{} {} STORE/PERATOM 0 1", id_fix, group->names[igroup])));
+    modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 0", id_fix, group->names[igroup])));
 
   // set fix store values = 0 for now
   // fill them in via setup() once Comm::borders() has been called

src/EXTRA-COMPUTE/compute_hma.cpp

@@ -91,7 +91,7 @@ ComputeHMA::ComputeHMA(LAMMPS *lmp, int narg, char **arg) :
 
   id_fix = utils::strdup(std::string(id)+"_COMPUTE_STORE");
   fix = dynamic_cast<FixStorePeratom *>(
-    modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 3", id_fix, group->names[igroup])));
+    modify->add_fix(fmt::format("{} {} STORE/PERATOM 3 0 0 1", id_fix, group->names[igroup])));
 
   // calculate xu,yu,zu for fix store array
   // skip if reset from restart file

src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp

@@ -17,7 +17,7 @@
 #include "comm.h"
 #include "domain.h"
 #include "error.h"
-#include "fix_property_atom.h"
+#include "fix_store_peratom.h"
 #include "force.h"
 #include "memory.h"
 #include "modify.h"
@@ -59,7 +59,7 @@ void BondHarmonicRestrain::compute(int eflag, int vflag)
   ev_init(eflag, vflag);
 
   double **x = atom->x;
-  double **x0 = (double **) atom->extract("d2_BOND_RESTRAIN_X0");
+  double **x0 = initial->astore;
   double **f = atom->f;
   int **bondlist = neighbor->bondlist;
   int nbondlist = neighbor->nbondlist;
@@ -162,12 +162,12 @@ void BondHarmonicRestrain::init_style()
   if (natoms < 0) {
 
     // create internal fix to store initial positions
-    initial = dynamic_cast<FixPropertyAtom *>(
+    initial = dynamic_cast<FixStorePeratom *>(
-        modify->add_fix("BOND_RESTRAIN_X0 all property/atom d2_BOND_RESTRAIN_X0 3 ghost yes"));
+      modify->add_fix("BOND_RESTRAIN_X0 all STORE/PERATOM 3 0 1 1"));
     if (!initial) error->all(FLERR, "Failure to create internal per-atom storage");
 
     natoms = atom->natoms;
-    double *const *const x0 = (double **) atom->extract("d2_BOND_RESTRAIN_X0");
+    double **x0 = initial->astore;
     const double *const *const x = atom->x;
     for (int i = 0; i < atom->nlocal; ++i)
       for (int j = 0; j < 3; ++j) x0[i][j] = x[i][j];
@@ -176,10 +176,10 @@ void BondHarmonicRestrain::init_style()
 
     // after a restart, natoms is set but initial is a null pointer.
     // we add the fix, but do not initialize it.  It will pull the data from the restart.
-
+    
     if (!initial) {
-      initial = dynamic_cast<FixPropertyAtom *>(
+      initial = dynamic_cast<FixStorePeratom *>(
-          modify->add_fix("BOND_RESTRAIN_X0 all property/atom d2_BOND_RESTRAIN_X0 3 ghost yes"));
+	modify->add_fix("BOND_RESTRAIN_X0 all STORE/PERATOM 3 0 1 1"));
       if (!initial) error->all(FLERR, "Failure to create internal per-atom storage");
     }
   }
@@ -230,7 +230,7 @@ void BondHarmonicRestrain::write_data(FILE *fp)
 
 double BondHarmonicRestrain::single(int type, double rsq, int i, int j, double &fforce)
 {
-  double **x0 = (double **) atom->extract("d2_BOND_RESTRAIN_X0");
+  double **x0 = initial->astore;
 
   double delx = x0[i][0] - x0[j][0];
   double dely = x0[i][1] - x0[j][1];

src/EXTRA-MOLECULE/bond_harmonic_restrain.h

@@ -41,7 +41,7 @@ class BondHarmonicRestrain : public Bond {
  protected:
   double *k;
   bigint natoms;
-  class FixPropertyAtom *initial;
+  class FixStorePeratom *initial;
 
   virtual void allocate();
 };

src/FEP/fix_adapt_fep.cpp

@@ -212,7 +212,7 @@ void FixAdaptFEP::post_constructor()
   if (diamflag) {
     id_fix_diam = utils::strdup(id + std::string("_FIX_STORE_DIAM"));
     fix_diam = dynamic_cast<FixStorePeratom *>(
-      modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 1", id_fix_diam,group->names[igroup])));
+      modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 1", id_fix_diam,group->names[igroup])));
     if (fix_diam->restart_reset) fix_diam->restart_reset = 0;
     else {
       double *vec = fix_diam->vstore;
@@ -230,7 +230,7 @@ void FixAdaptFEP::post_constructor()
   if (chgflag) {
     id_fix_chg = utils::strdup(id + std::string("_FIX_STORE_CHG"));
     fix_chg = dynamic_cast<FixStorePeratom *>(
-      modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 1",id_fix_chg,group->names[igroup])));
+      modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 1",id_fix_chg,group->names[igroup])));
     if (fix_chg->restart_reset) fix_chg->restart_reset = 0;
     else {
       double *vec = fix_chg->vstore;

src/REPLICA/tad.cpp

@@ -135,7 +135,7 @@ void TAD::command(int narg, char **arg)
 
   // create FixStorePeratom object to store revert state
 
-  fix_revert = dynamic_cast<FixStorePeratom *>(modify->add_fix("tad_revert all STORE/PERATOM 0 7"));
+  fix_revert = dynamic_cast<FixStorePeratom *>(modify->add_fix("tad_revert all STORE/PERATOM 7 0 0 0"));
 
   // create Finish for timing output
 

src/VTK/dump_vtk.cpp

@@ -2365,7 +2365,7 @@ int DumpVTK::modify_param(int narg, char **arg)
 
       std::string threshid = fmt::format("{}{}_DUMP_STORE",id,nthreshlast);
       thresh_fixID[nthreshlast] = utils::strdup(threshid);
-      threshid += fmt::format(" {} STORE/PERATOM 1 1", group->names[igroup]);
+      threshid += fmt::format(" {} STORE/PERATOM 1 0 0 1", group->names[igroup]);
       thresh_fix[nthreshlast] = dynamic_cast<FixStorePeratom *>(modify->add_fix(threshid));
 
       thresh_last[nthreshlast] = nthreshlast;

src/balance.cpp

@@ -515,7 +515,7 @@ void Balance::weight_storage(char *prefix)
 
   fixstore = dynamic_cast<FixStorePeratom *>(modify->get_fix_by_id(cmd));
   if (!fixstore)
-    fixstore = dynamic_cast<FixStorePeratom *>(modify->add_fix(cmd + " all STORE/PERATOM 0 1"));
+    fixstore = dynamic_cast<FixStorePeratom *>(modify->add_fix(cmd + " all STORE/PERATOM 1 0 0 0"));
 
   // do not carry weights with atoms during normal atom migration
 

src/compute_chunk_atom.cpp

@@ -572,7 +572,7 @@ void ComputeChunkAtom::init()
   if ((idsflag == ONCE || lockcount) && !fixstore) {
     id_fix = utils::strdup(id + std::string("_COMPUTE_STORE"));
     fixstore = dynamic_cast<FixStorePeratom *>(
-        modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 1", id_fix, group->names[igroup])));
+        modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 1", id_fix, group->names[igroup])));
   }
 
   if ((idsflag != ONCE && !lockcount) && fixstore) {

src/compute_displace_atom.cpp

@@ -75,7 +75,7 @@ ComputeDisplaceAtom::ComputeDisplaceAtom(LAMMPS *lmp, int narg, char **arg) :
 
   id_fix = utils::strdup(std::string(id) + "_COMPUTE_STORE");
   fix = dynamic_cast<FixStorePeratom *>(
-    modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 3", id_fix, group->names[igroup])));
+    modify->add_fix(fmt::format("{} {} STORE/PERATOM 3 0 0 1", id_fix, group->names[igroup])));
 
   // calculate xu,yu,zu for fix store array
   // skip if reset from restart file

src/compute_msd.cpp

@@ -64,7 +64,7 @@ ComputeMSD::ComputeMSD(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, a
 
   id_fix = utils::strdup(id + std::string("_COMPUTE_STORE"));
   fix = dynamic_cast<FixStorePeratom *>(
-      modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 3", id_fix, group->names[igroup])));
+      modify->add_fix(fmt::format("{} {} STORE/PERATOM 3 0 0 1", id_fix, group->names[igroup])));
 
   // calculate xu,yu,zu for fix store array
   // skip if reset from restart file

src/compute_vacf.cpp

@@ -40,7 +40,7 @@ ComputeVACF::ComputeVACF(LAMMPS *lmp, int narg, char **arg) :
 
   id_fix = utils::strdup(id + std::string("_COMPUTE_STORE"));
   fix = dynamic_cast<FixStorePeratom *>(
-    modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 3", id_fix, group->names[igroup])));
+    modify->add_fix(fmt::format("{} {} STORE/PERATOM 3 0 0 1", id_fix, group->names[igroup])));
 
   // store current velocities in fix store array
   // skip if reset from restart file

src/dump_custom.cpp

@@ -2013,7 +2013,7 @@ int DumpCustom::modify_param(int narg, char **arg)
 
       std::string threshid = fmt::format("{}{}_DUMP_STORE",id,nthreshlast);
       thresh_fixID[nthreshlast] = utils::strdup(threshid);
-      threshid += fmt::format(" {} STORE/PERATOM 1 1", group->names[igroup]);
+      threshid += fmt::format(" {} STORE/PERATOM 1 0 0 1", group->names[igroup]);
       thresh_fix[nthreshlast] = dynamic_cast<FixStorePeratom *>(modify->add_fix(threshid));
 
       thresh_last[nthreshlast] = nthreshlast;

src/fix_adapt.cpp

@@ -279,7 +279,7 @@ void FixAdapt::post_constructor()
   if (diamflag && atom->radius_flag) {
     id_fix_diam = utils::strdup(id + std::string("_FIX_STORE_DIAM"));
     fix_diam = dynamic_cast<FixStorePeratom *>(
-      modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 1", id_fix_diam,group->names[igroup])));
+      modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 1", id_fix_diam,group->names[igroup])));
     if (fix_diam->restart_reset) fix_diam->restart_reset = 0;
     else {
       double *vec = fix_diam->vstore;
@@ -297,7 +297,7 @@ void FixAdapt::post_constructor()
   if (chgflag && atom->q_flag) {
     id_fix_chg = utils::strdup(id + std::string("_FIX_STORE_CHG"));
     fix_chg = dynamic_cast<FixStorePeratom *>(
-      modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 1",id_fix_chg,group->names[igroup])));
+      modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 0 0 1",id_fix_chg,group->names[igroup])));
     if (fix_chg->restart_reset) fix_chg->restart_reset = 0;
     else {
       double *vec = fix_chg->vstore;

src/fix_store_peratom.cpp

@@ -22,50 +22,57 @@
 using namespace LAMMPS_NS;
 using namespace FixConst;
 
+// INTERNAL fix for storing/communicating per-atom values
+// syntax: id group style n1 n2 gflag rflag 
+//   N1 = 1, N2 = 0 is per-atom vector, single value per atom
+//   N1 > 1, N2 = 0 is per-atom array, N1 values per atom
+//   N1 > 0, N2 > 0 is per-atom tensor, N1xN2 array per atom
+//   gflag = 0/1, no/yes communicate per-atom values with ghost atoms
+//   rflag = 0/1, no/yes store per-atom values in restart file
+
 /* ---------------------------------------------------------------------- */
 
 FixStorePeratom::FixStorePeratom(LAMMPS *lmp, int narg, char **arg) :
     Fix(lmp, narg, arg), vstore(nullptr), astore(nullptr)
 {
-  if (narg != 5 && narg != 6)
+  if (narg != 8)
     error->all(FLERR, "Illegal fix STORE/PERATOM command: number of args");
 
-  // syntax: id group style 0/1 n1 n2 (n3), last arg optional
-  //   0/1 flag = not-store or store peratom values in restart file
-  //   N2 = 1 and no n3 is vector, N2 > 1 and no n3 is array, N3 = tensor
-  //   nvalues = # of peratom values, N = 1 is vector, N > 1 is array
-
   disable = 0;
-  vecflag = arrayflag = tensorflag = 0;
 
-  restart_peratom = utils::inumeric(FLERR, arg[3], false, lmp);
+  n1 = utils::inumeric(FLERR, arg[4], false, lmp);
-  n2 = utils::inumeric(FLERR, arg[4], false, lmp);
+  n2 = utils::inumeric(FLERR, arg[5], false, lmp);
-  if (narg == 6)
+  ghostflag = utils::inumeric(FLERR, arg[6], false, lmp);
-    n3 = utils::inumeric(FLERR, arg[5], false, lmp);
+  restartflag = utils::inumeric(FLERR, arg[7], false, lmp);
-  else
+
-    n3 = 1;
+  vecflag = arrayflag = tensorflag = 0;
-  if (restart_peratom < 0 || restart_peratom > 1)
+  if (n1 == 1 && n2 == 0) vecflag = 1;
-    error->all(FLERR, "Illegal fix STORE/PERATOM restart flag: {}", restart_peratom);
+  else if (n1 > 1 && n2 == 0) arrayflag = 1;
-  if (n2 <= 0 || n3 <= 0)
+  else if (n1 > 0 && n2 > 0) tensorflag = 1;
-    error->all(FLERR, "Illegal fix STORE/PERATOM dimension args: must be >0: {} {}", n2, n3);
+  else error->all(FLERR, "Illegal fix STORE/PERATOM dimension args: {} {}", n1, n2);
-  if (n2 == 1 && narg == 5)
+  
-    vecflag = 1;
+  if (restartflag < 0 || restartflag > 1)
-  else if (narg == 5)
+    error->all(FLERR, "Illegal fix STORE/PERATOM restart flag: {}", restartflag);
-    arrayflag = 1;
+  if (ghostflag < 0 || ghostflag > 1)
-  else
+    error->all(FLERR, "Illegal fix STORE/PERATOM ghost flag: {}", ghostflag);
-    tensorflag = 1;
+
-  nvalues = n2 * n3;
+  if (vecflag || arrayflag) nvalues = n1;
+  else nvalues = n1 * n2;
   nbytes = nvalues * sizeof(double);
 
+  if (ghostflag) comm_border = nvalues;
+  maxexchange = nvalues;
+  
+  // allocate data structs and register with Atom class
+
   vstore = nullptr;
   astore = nullptr;
   tstore = nullptr;
 
-  // allocate data structs and register with Atom class
-
   FixStorePeratom::grow_arrays(atom->nmax);
   atom->add_callback(Atom::GROW);
-  if (restart_peratom) atom->add_callback(Atom::RESTART);
+  if (restartflag) atom->add_callback(Atom::RESTART);
+  if (ghostflag) atom->add_callback(Atom::BORDER);
 
   // zero the storage
 
@@ -74,13 +81,12 @@ FixStorePeratom::FixStorePeratom(LAMMPS *lmp, int narg, char **arg) :
     for (int i = 0; i < nlocal; i++) vstore[i] = 0.0;
   } else if (arrayflag) {
     for (int i = 0; i < nlocal; i++)
-      for (int j = 0; j < n2; j++) astore[i][j] = 0.0;
+      for (int j = 0; j < n1; j++) astore[i][j] = 0.0;
   } else if (tensorflag) {
     for (int i = 0; i < nlocal; i++)
-      for (int j = 0; j < n2; j++)
+      for (int j = 0; j < n1; j++)
-        for (int k = 0; k < n3; k++) tstore[i][j][k] = 0.0;
+        for (int k = 0; k < n2; k++) tstore[i][j][k] = 0.0;
   }
-  maxexchange = nvalues;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -90,7 +96,8 @@ FixStorePeratom::~FixStorePeratom()
   // unregister callbacks to this fix from Atom class
 
   atom->delete_callback(id, Atom::GROW);
-  if (restart_peratom) atom->delete_callback(id, Atom::RESTART);
+  if (restartflag) atom->delete_callback(id, Atom::RESTART);
+  if (ghostflag) atom->delete_callback(id, Atom::BORDER);
 
   memory->destroy(vstore);
   memory->destroy(astore);
@@ -114,9 +121,9 @@ void FixStorePeratom::grow_arrays(int nmax)
   if (vecflag)
     memory->grow(vstore, nmax, "store:vstore");
   else if (arrayflag)
-    memory->grow(astore, nmax, n2, "store:astore");
+    memory->grow(astore, nmax, n1, "store:astore");
   else if (tensorflag)
-    memory->grow(tstore, nmax, n2, n3, "store:tstore");
+    memory->grow(tstore, nmax, n1, n2, "store:tstore");
 }
 
 /* ----------------------------------------------------------------------
@@ -136,6 +143,62 @@ void FixStorePeratom::copy_arrays(int i, int j, int /*delflag*/)
   }
 }
 
+/* ----------------------------------------------------------------------
+   pack values for border communication at re-neighboring
+------------------------------------------------------------------------- */
+
+int FixStorePeratom::pack_border(int n, int *list, double *buf)
+{
+  int i, j, k, ncol;
+
+  int m = 0;
+  if (vecflag) {
+    for (i = 0; i < n; i++) {
+      j = list[i];
+      buf[m++] = vstore[j];
+    }
+  } else if (arrayflag) {
+    for (i = 0; i < n; i++) {
+      j = list[i];
+      for (k = 0; k < nvalues; k++)
+	buf[m++] = astore[j][k];
+    }
+  } else if (tensorflag) {
+    for (i = 0; i < n; i++) {
+      j = list[i];
+      memcpy(&buf[m], &tstore[j][0][0], nbytes);
+      m += nvalues;
+    }
+  }
+
+  return m;
+}
+
+/* ----------------------------------------------------------------------
+   unpack values for border communication at re-neighboring
+------------------------------------------------------------------------- */
+
+int FixStorePeratom::unpack_border(int n, int first, double *buf)
+{
+  int i, k, last, ncol;
+
+  int m = 0;
+  last = first + n;
+  if (vecflag) {
+    for (i = first; i < last; i++)
+      vstore[i] = buf[m++];
+  } else if (arrayflag) {
+    for (i = first; i < last; i++)
+      for (k = 0; k < nvalues; k++)
+	astore[i][k] = buf[m++];
+  } else if (tensorflag) {
+    for (i = first; i < last; i++) {
+      memcpy(&tstore[i][0][0], &buf[m], nbytes);
+      m += nvalues;
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    pack values in local atom-based array for exchange with another proc
 ------------------------------------------------------------------------- */
@@ -147,7 +210,8 @@ int FixStorePeratom::pack_exchange(int i, double *buf)
   if (vecflag) {
     buf[0] = vstore[i];
   } else if (arrayflag) {
-    for (int m = 0; m < nvalues; m++) buf[m] = astore[i][m];
+    for (int m = 0; m < nvalues; m++)
+      buf[m] = astore[i][m];
   } else if (tensorflag) {
     memcpy(buf, &tstore[i][0][0], nbytes);
   }
@@ -166,7 +230,8 @@ int FixStorePeratom::unpack_exchange(int nlocal, double *buf)
   if (vecflag) {
     vstore[nlocal] = buf[0];
   } else if (arrayflag) {
-    for (int m = 0; m < nvalues; m++) astore[nlocal][m] = buf[m];
+    for (int m = 0; m < nvalues; m++)
+      astore[nlocal][m] = buf[m];
   } else if (tensorflag) {
     memcpy(&tstore[nlocal][0][0], buf, nbytes);
   }
@@ -191,7 +256,8 @@ int FixStorePeratom::pack_restart(int i, double *buf)
   if (vecflag) {
     buf[1] = vstore[i];
   } else if (arrayflag) {
-    for (int m = 0; m < nvalues; m++) buf[m + 1] = astore[i][m];
+    for (int m = 0; m < nvalues; m++)
+      buf[m + 1] = astore[i][m];
   } else if (tensorflag) {
     memcpy(&buf[1], &tstore[i][0][0], nbytes);
   }
@@ -219,7 +285,8 @@ void FixStorePeratom::unpack_restart(int nlocal, int nth)
   if (vecflag) {
     vstore[nlocal] = extra[nlocal][m];
   } else if (arrayflag) {
-    for (int i = 0; i < nvalues; i++) astore[nlocal][i] = extra[nlocal][m++];
+    for (int i = 0; i < nvalues; i++)
+      astore[nlocal][i] = extra[nlocal][m++];
   } else if (tensorflag) {
     memcpy(&tstore[nlocal][0][0], &extra[nlocal][m], nbytes);
   }
@@ -251,5 +318,5 @@ int FixStorePeratom::size_restart(int /*nlocal*/)
 
 double FixStorePeratom::memory_usage()
 {
-  return (double) atom->nmax * n2 * n3 * sizeof(double);
+  return (double) atom->nmax * nvalues * sizeof(double);
 }

src/fix_store_peratom.h

@@ -37,6 +37,8 @@ class FixStorePeratom : public Fix {
 
   void grow_arrays(int) override;
   void copy_arrays(int, int, int) override;
+  int pack_border(int, int *, double *) override;
+  int unpack_border(int, int, double *) override;
   int pack_exchange(int, double *) override;
   int unpack_exchange(int, double *) override;
   int pack_restart(int, double *) override;
@@ -50,8 +52,10 @@ class FixStorePeratom : public Fix {
   int vecflag;       // 1 if ncol=1
   int arrayflag;     // 1 if a 2d array (vector per atom)
   int tensorflag;    // 1 if a 3d array (array per atom)
-
+  int ghostflag;     // 0/1 to communicate values with ghost atoms
-  int n2, n3;     // size of 3d dims of per-atom data struct
+  int restartflag;   // 0/1 to store values in restart files
+  
+  int n1, n2;     // size of 3d dims of per-atom data struct
   int nvalues;    // number of per-atom values
   int nbytes;     // number of per-atom bytes
 };

src/variable.cpp

@@ -5028,7 +5028,7 @@ VarReader::VarReader(LAMMPS *lmp, char *name, char *file, int flag) :
 
     id_fix = utils::strdup(std::string(name) + "_VARIABLE_STORE");
     fixstore = dynamic_cast<FixStorePeratom *>(
-      modify->add_fix(std::string(id_fix) + " all STORE/PERATOM 0 1"));
+      modify->add_fix(std::string(id_fix) + " all STORE/PERATOM 1 0 0 0"));
     buffer = new char[CHUNK*MAXLINE];
   }
 }