diff --git a/src/AMOEBA/amoeba_induce.cpp b/src/AMOEBA/amoeba_induce.cpp
index 43688cef94..c1c2cb1d51 100644
--- a/src/AMOEBA/amoeba_induce.cpp
+++ b/src/AMOEBA/amoeba_induce.cpp
@@ -20,7 +20,7 @@
 #include "domain.h"
 #include "error.h"
 #include "fft3d_wrap.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "math_const.h"
 #include "math_special.h"
 #include "memory.h"
diff --git a/src/AMOEBA/amoeba_utils.cpp b/src/AMOEBA/amoeba_utils.cpp
index 332b62708e..ac0e89a175 100644
--- a/src/AMOEBA/amoeba_utils.cpp
+++ b/src/AMOEBA/amoeba_utils.cpp
@@ -18,7 +18,7 @@
 #include "comm.h"
 #include "domain.h"
 #include "error.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "neigh_list.h"
 
 #include <cmath>
diff --git a/src/AMOEBA/pair_amoeba.cpp b/src/AMOEBA/pair_amoeba.cpp
index e8b7b71753..bf52524567 100644
--- a/src/AMOEBA/pair_amoeba.cpp
+++ b/src/AMOEBA/pair_amoeba.cpp
@@ -21,7 +21,7 @@
 #include "error.h"
 #include "fft3d_wrap.h"
 #include "fix.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "force.h"
 #include "gridcomm.h"
 #include "group.h"
@@ -783,8 +783,8 @@ 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]));
-    fixpole = dynamic_cast<FixStore *>(myfix);
+    myfix = modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 13",id_pole,group->names[0]));
+    fixpole = dynamic_cast<FixStorePeratom *>(myfix);
   }
 
   // creation of per-atom storage
@@ -795,14 +795,14 @@ 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 1 {} 3",
                                         id_udalt, group->names[0], maxualt));
-    fixudalt = dynamic_cast<FixStore *>(myfix);
+    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 1 {} 3",
                                         id_upalt, group->names[0], maxualt));
-    fixupalt = dynamic_cast<FixStore *>(myfix);
+    fixupalt = dynamic_cast<FixStorePeratom *>(myfix);
   }
 
   // create pages for storing pairwise data:
@@ -916,19 +916,22 @@ void PairAmoeba::init_style()
 
   if (id_pole) {
     myfix = modify->get_fix_by_id(id_pole);
-    if (!myfix) error->all(FLERR,"Could not find internal pair amoeba fix STORE id {}", id_pole);
-    fixpole = dynamic_cast<FixStore *>(myfix);
+    if (!myfix)
+      error->all(FLERR,"Could not find internal pair amoeba fix STORE/PERATOM id {}", id_pole);
+    fixpole = dynamic_cast<FixStorePeratom *>(myfix);
 
   }
 
   if (id_udalt) {
     myfix = modify->get_fix_by_id(id_udalt);
-    if (!myfix) error->all(FLERR,"Could not find internal pair amoeba fix STORE id {}", id_udalt);
-    fixudalt = dynamic_cast<FixStore *>(myfix);
+    if (!myfix)
+      error->all(FLERR,"Could not find internal pair amoeba fix STORE/PERATOM id {}", id_udalt);
+    fixudalt = dynamic_cast<FixStorePeratom *>(myfix);
 
     myfix = modify->get_fix_by_id(id_upalt);
-    if (!myfix) error->all(FLERR,"Could not find internal pair amoeba fix STORE id {}", id_upalt);
-    fixupalt = dynamic_cast<FixStore *>(myfix);
+    if (!myfix)
+      error->all(FLERR,"Could not find internal pair amoeba fix STORE/PERATOM id {}", id_upalt);
+    fixupalt = dynamic_cast<FixStorePeratom *>(myfix);
   }
 
   // assign hydrogen neighbors (redID) to each owned atom
diff --git a/src/AMOEBA/pair_amoeba.h b/src/AMOEBA/pair_amoeba.h
index 7b0719bca0..84bc480062 100644
--- a/src/AMOEBA/pair_amoeba.h
+++ b/src/AMOEBA/pair_amoeba.h
@@ -157,9 +157,9 @@ class PairAmoeba : public Pair {
   int *amgroup;    // AMOEBA polarization group, 1 to Ngroup
 
   char *id_pole, *id_udalt, *id_upalt;
-  class FixStore *fixpole;     // stores pole = multipole components
-  class FixStore *fixudalt;    // stores udalt = induced dipole history
-  class FixStore *fixupalt;    // stores upalt = induced dipole history
+  class FixStorePeratom *fixpole;     // stores pole = multipole components
+  class FixStorePeratom *fixudalt;    // stores udalt = induced dipole history
+  class FixStorePeratom *fixupalt;    // stores upalt = induced dipole history
 
   // static per-type properties defined in force-field file
 
diff --git a/src/CORESHELL/compute_temp_cs.cpp b/src/CORESHELL/compute_temp_cs.cpp
index 7f6fba5f6a..e9778489fd 100644
--- a/src/CORESHELL/compute_temp_cs.cpp
+++ b/src/CORESHELL/compute_temp_cs.cpp
@@ -24,7 +24,7 @@
 #include "comm.h"
 #include "domain.h"
 #include "error.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "force.h"
 #include "group.h"
 #include "memory.h"
@@ -69,11 +69,9 @@ ComputeTempCS::ComputeTempCS(LAMMPS *lmp, int narg, char **arg) :
   // create a new fix STORE style
   // id = compute-ID + COMPUTE_STORE, fix group = compute group
 
-  std::string fixcmd = id + std::string("_COMPUTE_STORE");
-  id_fix = utils::strdup(fixcmd);
-
-  fixcmd += fmt::format(" {} STORE peratom 0 1", group->names[igroup]);
-  fix = dynamic_cast<FixStore *>(modify->add_fix(fixcmd));
+  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])));
 
   // set fix store values = 0 for now
   // fill them in via setup() once Comm::borders() has been called
@@ -108,8 +106,8 @@ ComputeTempCS::~ComputeTempCS()
 
   if (modify->nfix) modify->delete_fix(id_fix);
 
-  delete [] id_fix;
-  delete [] vector;
+  delete[] id_fix;
+  delete[] vector;
   memory->destroy(vint);
 }
 
diff --git a/src/CORESHELL/compute_temp_cs.h b/src/CORESHELL/compute_temp_cs.h
index 96ad6ca91a..d8368827d5 100644
--- a/src/CORESHELL/compute_temp_cs.h
+++ b/src/CORESHELL/compute_temp_cs.h
@@ -54,7 +54,7 @@ class ComputeTempCS : public Compute {
   double **vint;
 
   char *id_fix;
-  class FixStore *fix;
+  class FixStorePeratom *fix;
 
   void dof_compute();
   void vcm_pairs();
diff --git a/src/EXTRA-COMPUTE/compute_hma.cpp b/src/EXTRA-COMPUTE/compute_hma.cpp
index bb182de7af..12d8325957 100644
--- a/src/EXTRA-COMPUTE/compute_hma.cpp
+++ b/src/EXTRA-COMPUTE/compute_hma.cpp
@@ -54,7 +54,7 @@ https://doi.org/10.1103/PhysRevE.92.043303
 #include "domain.h"
 #include "error.h"
 #include "fix.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "force.h"
 #include "group.h"
 #include "improper.h"
@@ -90,8 +90,8 @@ ComputeHMA::ComputeHMA(LAMMPS *lmp, int narg, char **arg) :
   // our new fix's group = same as compute group
 
   id_fix = utils::strdup(std::string(id)+"_COMPUTE_STORE");
-  fix = dynamic_cast<FixStore *>(modify->add_fix(fmt::format("{} {} STORE peratom 1 3",
-                                                id_fix, group->names[igroup])));
+  fix = dynamic_cast<FixStorePeratom *>(
+    modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 3", id_fix, group->names[igroup])));
 
   // calculate xu,yu,zu for fix store array
   // skip if reset from restart file
@@ -196,9 +196,8 @@ void ComputeHMA::setup()
 
   // set fix which stores original atom coords
 
-  int ifix2 = modify->find_fix(id_fix);
-  if (ifix2 < 0) error->all(FLERR,"Could not find hma store fix ID");
-  fix = dynamic_cast<FixStore *>( modify->fix[ifix2]);
+  fix = dynamic_cast<FixStorePeratom *>(modify->get_fix_by_id(id_fix));
+  if (!fix) error->all(FLERR,"Could not find hma per-atom store fix ID {}", id_fix);
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/EXTRA-COMPUTE/compute_hma.h b/src/EXTRA-COMPUTE/compute_hma.h
index d15c6f4317..12af0e809e 100644
--- a/src/EXTRA-COMPUTE/compute_hma.h
+++ b/src/EXTRA-COMPUTE/compute_hma.h
@@ -43,7 +43,7 @@ class ComputeHMA : public Compute {
   char *id_fix;
   char *id_temp;
   double finaltemp;
-  class FixStore *fix;
+  class FixStorePeratom *fix;
   double boltz, nktv2p, inv_volume;
   double deltaPcap;
   double virial_compute(int);
diff --git a/src/EXTRA-COMPUTE/compute_msd_nongauss.cpp b/src/EXTRA-COMPUTE/compute_msd_nongauss.cpp
index 0f68e03ca7..24920e50cf 100644
--- a/src/EXTRA-COMPUTE/compute_msd_nongauss.cpp
+++ b/src/EXTRA-COMPUTE/compute_msd_nongauss.cpp
@@ -19,7 +19,7 @@
 
 #include "atom.h"
 #include "domain.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "group.h"
 #include "update.h"
 
diff --git a/src/EXTRA-FIX/fix_npt_cauchy.cpp b/src/EXTRA-FIX/fix_npt_cauchy.cpp
index 4a1e11d963..c066d471d6 100644
--- a/src/EXTRA-FIX/fix_npt_cauchy.cpp
+++ b/src/EXTRA-FIX/fix_npt_cauchy.cpp
@@ -24,7 +24,7 @@
 #include "domain.h"
 #include "error.h"
 #include "fix_deform.h"
-#include "fix_store.h"
+#include "fix_store_global.h"
 #include "force.h"
 #include "group.h"
 #include "irregular.h"
@@ -59,7 +59,7 @@ FixNPTCauchy::FixNPTCauchy(LAMMPS *lmp, int narg, char **arg) :
   id_temp(nullptr), id_press(nullptr),
   eta(nullptr), eta_dot(nullptr), eta_dotdot(nullptr),
   eta_mass(nullptr), etap(nullptr), etap_dot(nullptr), etap_dotdot(nullptr),
-  etap_mass(nullptr), id_store(nullptr),init_store(nullptr)
+  etap_mass(nullptr), id_store(nullptr), init_store(nullptr)
 {
   if (narg < 4) error->all(FLERR,"Illegal fix npt/cauchy command");
 
@@ -78,7 +78,6 @@ FixNPTCauchy::FixNPTCauchy(LAMMPS *lmp, int narg, char **arg) :
   initRUN = 0;
   restartPK = 0;
   restart_global = 1;
-  restart_stored = 0;
 
   // default values
 
@@ -275,7 +274,7 @@ FixNPTCauchy::FixNPTCauchy(LAMMPS *lmp, int narg, char **arg) :
       if (strcmp(arg[iarg+1],"all") == 0) allremap = 1;
       else {
         allremap = 0;
-        delete [] id_dilate;
+        delete[] id_dilate;
         id_dilate = utils::strdup(arg[iarg+1]);
         int idilate = group->find(id_dilate);
         if (idilate == -1)
@@ -619,32 +618,32 @@ FixNPTCauchy::~FixNPTCauchy()
 {
   if (copymode) return;
 
-  delete [] id_dilate;
-  delete [] rfix;
+  delete[] id_dilate;
+  delete[] rfix;
 
-  delete [] id_store;
+  delete[] id_store;
   delete irregular;
 
   // delete temperature and pressure if fix created them
 
   if (tcomputeflag) modify->delete_compute(id_temp);
-  delete [] id_temp;
+  delete[] id_temp;
 
   if (tstat_flag) {
-    delete [] eta;
-    delete [] eta_dot;
-    delete [] eta_dotdot;
-    delete [] eta_mass;
+    delete[] eta;
+    delete[] eta_dot;
+    delete[] eta_dotdot;
+    delete[] eta_mass;
   }
 
   if (pstat_flag) {
     if (pcomputeflag) modify->delete_compute(id_press);
-    delete [] id_press;
+    delete[] id_press;
     if (mpchain) {
-      delete [] etap;
-      delete [] etap_dot;
-      delete [] etap_dotdot;
-      delete [] etap_mass;
+      delete[] etap;
+      delete[] etap_dot;
+      delete[] etap_dotdot;
+      delete[] etap_mass;
     }
   }
 }
@@ -762,7 +761,7 @@ void FixNPTCauchy::init()
   // detect if any rigid fixes exist so rigid bodies move when box is remapped
   // rfix[] = indices to each fix rigid
 
-  delete [] rfix;
+  delete[] rfix;
   nrigid = 0;
   rfix = nullptr;
 
@@ -1430,7 +1429,7 @@ int FixNPTCauchy::modify_param(int narg, char **arg)
       modify->delete_compute(id_temp);
       tcomputeflag = 0;
     }
-    delete [] id_temp;
+    delete[] id_temp;
     id_temp = utils::strdup(arg[1]);
 
     int icompute = modify->find_compute(arg[1]);
@@ -1439,8 +1438,7 @@ int FixNPTCauchy::modify_param(int narg, char **arg)
     temperature = modify->compute[icompute];
 
     if (temperature->tempflag == 0)
-      error->all(FLERR,
-                 "Fix_modify temperature ID does not compute temperature");
+      error->all(FLERR,"Fix_modify temperature ID does not compute temperature");
     if (temperature->igroup != 0 && comm->me == 0)
       error->warning(FLERR,"Temperature for fix modify is not for group all");
 
@@ -1462,7 +1460,7 @@ int FixNPTCauchy::modify_param(int narg, char **arg)
       modify->delete_compute(id_press);
       pcomputeflag = 0;
     }
-    delete [] id_press;
+    delete[] id_press;
     id_press = utils::strdup(arg[1]);
 
     int icompute = modify->find_compute(arg[1]);
@@ -2451,22 +2449,19 @@ void FixNPTCauchy::CauchyStat_init()
     utils::logmesg(lmp, mesg);
   }
 
-  if (!id_store)
-    id_store = utils::strdup(std::string(id) + "_FIX_NH_STORE");
-  restart_stored = modify->find_fix(id_store);
+  if (!id_store) id_store = utils::strdup(std::string(id) + "_FIX_NH_STORE");
+  init_store = dynamic_cast<FixStoreGlobal *>(modify->get_fix_by_id(id_store));
 
-  if (restartPK==1 && restart_stored < 0)
-    error->all(FLERR,"Illegal npt/cauchy command.  Continuation run"
+  if ((restartPK == 1) && !init_store)
+    error->all(FLERR,"Illegal fix npt/cauchy command.  Continuation run"
                " must follow a previously equilibrated npt/cauchy run");
 
   if (alpha<=0.0)
     error->all(FLERR,"Illegal fix npt/cauchy command: Alpha cannot be zero or negative.");
 
-  if (restart_stored < 0) {
-    modify->add_fix(std::string(id_store) + " all STORE global 1 6");
-    restart_stored = modify->find_fix(id_store);
-  }
-  init_store = dynamic_cast<FixStore *>(modify->fix[restart_stored]);
+  if (!init_store)
+    init_store = dynamic_cast<FixStoreGlobal *>(
+      modify->add_fix(std::string(id_store) + " all STORE global 1 6"));
 
   initRUN = 0;
   initPK = 1;
diff --git a/src/EXTRA-FIX/fix_npt_cauchy.h b/src/EXTRA-FIX/fix_npt_cauchy.h
index 593520cd41..b17d039724 100644
--- a/src/EXTRA-FIX/fix_npt_cauchy.h
+++ b/src/EXTRA-FIX/fix_npt_cauchy.h
@@ -149,24 +149,23 @@ class FixNPTCauchy : public Fix {
   void nh_omega_dot();
 
   // Implementation of CauchyStat
-  char *id_store;                // fix id of the STORE fix for retaining data
-  class FixStore *init_store;    // fix pointer to STORE fix
-  double H0[3][3];               // shape matrix for the undeformed cell
-  double h_old[6];               // previous time step shape matrix for
-                                 // the undeformed cell
-  double invH0[3][3];            // inverse of H0;
-  double CSvol0;                 // volume of undeformed cell
-  double setPK[3][3];            // current set values of the PK stress
-                                 // (this is modified until the cauchy
-                                 // stress converges)
-  double alpha;                  // integration parameter for the cauchystat
-  int initPK;                    // 1 if setPK needs to be initialized either
-                                 // from cauchy or restart, else 0
-  int restartPK;                 // Read PK stress from the previous run
-  int restart_stored;            // values of PK stress from the previous step stored
-  int initRUN;                   // 0 if run not initialized
-                                 // (pressure->vector not computed yet),
-                                 // else 1 (pressure->vector available)
+  char *id_store;                      // fix id of the STORE fix for retaining data
+  class FixStoreGlobal *init_store;    // fix pointer to STORE fix
+  double H0[3][3];                     // shape matrix for the undeformed cell
+  double h_old[6];                     // previous time step shape matrix for
+                                       // the undeformed cell
+  double invH0[3][3];                  // inverse of H0;
+  double CSvol0;                       // volume of undeformed cell
+  double setPK[3][3];                  // current set values of the PK stress
+                                       // (this is modified until the cauchy
+                                       // stress converges)
+  double alpha;                        // integration parameter for the cauchystat
+  int initPK;                          // 1 if setPK needs to be initialized either
+                                       // from cauchy or restart, else 0
+  int restartPK;                       // Read PK stress from the previous run
+  int initRUN;                         // 0 if run not initialized
+                                       // (pressure->vector not computed yet),
+                                       // else 1 (pressure->vector available)
 
   void CauchyStat_init();
   void CauchyStat_cleanup();
diff --git a/src/FEP/fix_adapt_fep.cpp b/src/FEP/fix_adapt_fep.cpp
index 090931ae31..0fad4f362d 100644
--- a/src/FEP/fix_adapt_fep.cpp
+++ b/src/FEP/fix_adapt_fep.cpp
@@ -20,7 +20,7 @@
 
 #include "atom.h"
 #include "error.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "force.h"
 #include "group.h"
 #include "input.h"
@@ -210,9 +210,9 @@ void FixAdaptFEP::post_constructor()
   id_fix_chg = nullptr;
 
   if (diamflag) {
-    auto cmd = fmt::format("{}_FIX_STORE_DIAM {} STORE peratom 1 1", group->names[igroup]);
-    fix_diam = dynamic_cast<FixStore *>( modify->add_fix(cmd));
-
+    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])));
     if (fix_diam->restart_reset) fix_diam->restart_reset = 0;
     else {
       double *vec = fix_diam->vstore;
@@ -228,9 +228,9 @@ void FixAdaptFEP::post_constructor()
   }
 
   if (chgflag) {
-    auto cmd = fmt::format("{}_FIX_STORE_CHG {} STORE peratom 1 1", group->names[igroup]);
-    fix_chg = dynamic_cast<FixStore *>( modify->add_fix(cmd));
-
+    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])));
     if (fix_chg->restart_reset) fix_chg->restart_reset = 0;
     else {
       double *vec = fix_chg->vstore;
@@ -333,18 +333,16 @@ void FixAdaptFEP::init()
   // fixes that store initial per-atom values
 
   if (id_fix_diam) {
-    int ifix = modify->find_fix(id_fix_diam);
-    if (ifix < 0) error->all(FLERR,"Could not find fix adapt storage fix ID");
-    fix_diam = dynamic_cast<FixStore *>( modify->fix[ifix]);
+    fix_diam = dynamic_cast<FixStorePeratom *>(modify->get_fix_by_id(id_fix_diam));
+    if (!fix_diam) error->all(FLERR,"Could not find fix adapt/fep storage fix ID {}", id_fix_diam);
   }
   if (id_fix_chg) {
-    int ifix = modify->find_fix(id_fix_chg);
-    if (ifix < 0) error->all(FLERR,"Could not find fix adapt storage fix ID");
-    fix_chg = dynamic_cast<FixStore *>( modify->fix[ifix]);
+    fix_chg = dynamic_cast<FixStorePeratom *>(modify->get_fix_by_id(id_fix_chg));
+    if (!fix_chg) error->all(FLERR,"Could not find fix adapt/fep storage fix ID {}", id_fix_chg);
   }
 
   if (utils::strmatch(update->integrate_style,"^respa"))
-    nlevels_respa = (dynamic_cast<Respa *>( update->integrate))->nlevels;
+    nlevels_respa = (dynamic_cast<Respa *>(update->integrate))->nlevels;
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/FEP/fix_adapt_fep.h b/src/FEP/fix_adapt_fep.h
index ca0255962c..86b6a93733 100644
--- a/src/FEP/fix_adapt_fep.h
+++ b/src/FEP/fix_adapt_fep.h
@@ -46,7 +46,7 @@ class FixAdaptFEP : public Fix {
   int anypair;
   int nlevels_respa;
   char *id_fix_diam, *id_fix_chg;
-  class FixStore *fix_diam, *fix_chg;
+  class FixStorePeratom *fix_diam, *fix_chg;
 
   struct Adapt {
     int which, ivar;
diff --git a/src/REPLICA/tad.cpp b/src/REPLICA/tad.cpp
index 91ea2a0a29..490e98f168 100644
--- a/src/REPLICA/tad.cpp
+++ b/src/REPLICA/tad.cpp
@@ -24,7 +24,7 @@
 #include "error.h"
 #include "finish.h"
 #include "fix_event_tad.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "force.h"
 #include "integrate.h"
 #include "memory.h"
@@ -54,9 +54,9 @@ TAD::TAD(LAMMPS *lmp) : Command(lmp)
 TAD::~TAD()
 {
   memory->sfree(fix_event_list);
-  if (neb_logfilename != nullptr) delete [] neb_logfilename;
-  delete [] min_style;
-  delete [] min_style_neb;
+  if (neb_logfilename != nullptr) delete[] neb_logfilename;
+  delete[] min_style;
+  delete[] min_style_neb;
 }
 
 /* ----------------------------------------------------------------------
@@ -131,11 +131,11 @@ void TAD::command(int narg, char **arg)
 
   // create FixEventTAD object to store last event
 
-  fix_event = dynamic_cast<FixEventTAD *>( modify->add_fix("tad_event all EVENT/TAD"));
+  fix_event = dynamic_cast<FixEventTAD *>(modify->add_fix("tad_event all EVENT/TAD"));
 
-  // create FixStore object to store revert state
+  // create FixStorePeratom object to store revert state
 
-  fix_revert = dynamic_cast<FixStore *>( modify->add_fix("tad_revert all STORE peratom 0 7"));
+  fix_revert = dynamic_cast<FixStorePeratom *>(modify->add_fix("tad_revert all STORE/PERATOM 0 7"));
 
   // create Finish for timing output
 
@@ -388,7 +388,7 @@ void TAD::command(int narg, char **arg)
   neighbor->dist_check = neigh_dist_check;
 
 
-  delete [] id_compute;
+  delete[] id_compute;
   delete finish;
   modify->delete_fix("tad_event");
   modify->delete_fix("tad_revert");
@@ -578,7 +578,7 @@ void TAD::options(int narg, char **arg)
 
     } else if (strcmp(arg[iarg],"neb_style") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal tad command");
-      delete [] min_style_neb;
+      delete[] min_style_neb;
       min_style_neb = utils::strdup(arg[iarg+1]);
       iarg += 2;
 
@@ -589,7 +589,7 @@ void TAD::options(int narg, char **arg)
       iarg += 2;
 
     } else if (strcmp(arg[iarg],"neb_log") == 0) {
-      delete [] neb_logfilename;
+      delete[] neb_logfilename;
       if (iarg+2 > narg) error->all(FLERR,"Illegal tad command");
       if (strcmp(arg[iarg+1],"none") == 0) neb_logfilename = nullptr;
       else {
diff --git a/src/REPLICA/tad.h b/src/REPLICA/tad.h
index b16477096c..07c8764686 100644
--- a/src/REPLICA/tad.h
+++ b/src/REPLICA/tad.h
@@ -48,15 +48,15 @@ class TAD : public Command {
   double time_dynamics, time_quench, time_neb, time_comm, time_output;
   double time_start;
 
-  class NEB *neb;                  // NEB object
-  class Fix *fix_neb;              // FixNEB object
-  class Compute *compute_event;    // compute to detect event
-  class FixEventTAD *fix_event;    // current event/state
-  class FixStore *fix_revert;      // revert state
-  FixEventTAD **fix_event_list;    // list of possible events
-  int n_event_list;                // number of events
-  int nmax_event_list;             // allocated events
-  int nmin_event_list;             // minimum allocation
+  class NEB *neb;                       // NEB object
+  class Fix *fix_neb;                   // FixNEB object
+  class Compute *compute_event;         // compute to detect event
+  class FixEventTAD *fix_event;         // current event/state
+  class FixStorePeratom *fix_revert;    // revert state
+  FixEventTAD **fix_event_list;         // list of possible events
+  int n_event_list;                     // number of events
+  int nmax_event_list;                  // allocated events
+  int nmin_event_list;                  // minimum allocation
 
   char *neb_logfilename;    // filename for ulogfile_neb
   FILE *uscreen_neb;        // neb universe screen output
diff --git a/src/VTK/dump_vtk.cpp b/src/VTK/dump_vtk.cpp
index f346badfb5..61897735c7 100644
--- a/src/VTK/dump_vtk.cpp
+++ b/src/VTK/dump_vtk.cpp
@@ -31,7 +31,7 @@
 #include "domain.h"
 #include "error.h"
 #include "fix.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "force.h"
 #include "group.h"
 #include "input.h"
@@ -2357,16 +2357,16 @@ int DumpVTK::modify_param(int narg, char **arg)
       thresh_value[nthresh] = utils::numeric(FLERR,arg[3],false,lmp);
       thresh_last[nthresh] = -1;
     } else {
-      thresh_fix = (FixStore **)
-        memory->srealloc(thresh_fix,(nthreshlast+1)*sizeof(FixStore *),"dump:thresh_fix");
+      thresh_fix = (FixStorePeratom **)
+        memory->srealloc(thresh_fix,(nthreshlast+1)*sizeof(FixStorePeratom *),"dump:thresh_fix");
       thresh_fixID = (char **)
         memory->srealloc(thresh_fixID,(nthreshlast+1)*sizeof(char *),"dump:thresh_fixID");
       memory->grow(thresh_first,(nthreshlast+1),"dump:thresh_first");
 
       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]);
-      thresh_fix[nthreshlast] = (FixStore *) modify->add_fix(threshid);
+      threshid += fmt::format(" {} STORE/PERATOM 1 1", group->names[igroup]);
+      thresh_fix[nthreshlast] = dynamic_cast<FixStorePeratom *>(modify->add_fix(threshid));
 
       thresh_last[nthreshlast] = nthreshlast;
       thresh_first[nthreshlast] = 1;
diff --git a/src/balance.cpp b/src/balance.cpp
index 0e75419433..c67e561738 100644
--- a/src/balance.cpp
+++ b/src/balance.cpp
@@ -26,7 +26,7 @@
 #include "neighbor.h"
 #include "comm.h"
 #include "domain.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "imbalance.h"
 #include "imbalance_group.h"
 #include "imbalance_neigh.h"
@@ -78,26 +78,26 @@ Balance::~Balance()
   memory->destroy(proccost);
   memory->destroy(allproccost);
 
-  delete [] user_xsplit;
-  delete [] user_ysplit;
-  delete [] user_zsplit;
+  delete[] user_xsplit;
+  delete[] user_ysplit;
+  delete[] user_zsplit;
 
   if (shift_allocate) {
-    delete [] bdim;
-    delete [] onecost;
-    delete [] allcost;
-    delete [] sum;
-    delete [] target;
-    delete [] lo;
-    delete [] hi;
-    delete [] losum;
-    delete [] hisum;
+    delete[] bdim;
+    delete[] onecost;
+    delete[] allcost;
+    delete[] sum;
+    delete[] target;
+    delete[] lo;
+    delete[] hi;
+    delete[] losum;
+    delete[] hisum;
   }
 
   delete rcb;
 
   for (int i = 0; i < nimbalance; i++) delete imbalances[i];
-  delete [] imbalances;
+  delete[] imbalances;
 
   // check nfix in case all fixes have already been deleted
 
@@ -143,7 +143,7 @@ void Balance::command(int narg, char **arg)
         if (1 + procgrid[0]-1 > narg)
           error->all(FLERR,"Illegal balance command");
         xflag = USER;
-        delete [] user_xsplit;
+        delete[] user_xsplit;
         user_xsplit = new double[procgrid[0]+1];
         user_xsplit[0] = 0.0;
         iarg++;
@@ -163,7 +163,7 @@ void Balance::command(int narg, char **arg)
         if (1 + procgrid[1]-1 > narg)
           error->all(FLERR,"Illegal balance command");
         yflag = USER;
-        delete [] user_ysplit;
+        delete[] user_ysplit;
         user_ysplit = new double[procgrid[1]+1];
         user_ysplit[0] = 0.0;
         iarg++;
@@ -183,7 +183,7 @@ void Balance::command(int narg, char **arg)
         if (1 + procgrid[2]-1 > narg)
           error->all(FLERR,"Illegal balance command");
         zflag = USER;
-        delete [] user_zsplit;
+        delete[] user_zsplit;
         user_zsplit = new double[procgrid[2]+1];
         user_zsplit[0] = 0.0;
         iarg++;
@@ -496,8 +496,9 @@ void Balance::weight_storage(char *prefix)
   if (prefix) cmd = prefix;
   cmd += "IMBALANCE_WEIGHTS";
 
-  fixstore = dynamic_cast<FixStore *>( modify->get_fix_by_id(cmd));
-  if (!fixstore) fixstore = dynamic_cast<FixStore *>( modify->add_fix(cmd + " all STORE peratom 0 1"));
+  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"));
 
   // do not carry weights with atoms during normal atom migration
 
diff --git a/src/balance.h b/src/balance.h
index 252960743b..29870ea52a 100644
--- a/src/balance.h
+++ b/src/balance.h
@@ -27,10 +27,10 @@ namespace LAMMPS_NS {
 class Balance : public Command {
  public:
   class RCB *rcb;
-  class FixStore *fixstore;    // per-atom weights stored in FixStore
-  int wtflag;                  // 1 if particle weighting is used
-  int varflag;                 // 1 if weight style var(iable) is used
-  int outflag;                 // 1 for output of balance results to file
+  class FixStorePeratom *fixstore;    // per-atom weights stored in FixStorePeratom
+  int wtflag;                         // 1 if particle weighting is used
+  int varflag;                        // 1 if weight style var(iable) is used
+  int outflag;                        // 1 for output of balance results to file
 
   Balance(class LAMMPS *);
   ~Balance() override;
@@ -90,8 +90,6 @@ class Balance : public Command {
   void debug_shift_output(int, int, int, double *);
 #endif
 };
-
 }    // namespace LAMMPS_NS
-
 #endif
 #endif
diff --git a/src/compute_chunk_atom.cpp b/src/compute_chunk_atom.cpp
index 4bc92ea476..f16ab05598 100644
--- a/src/compute_chunk_atom.cpp
+++ b/src/compute_chunk_atom.cpp
@@ -21,7 +21,7 @@
 #include "domain.h"
 #include "error.h"
 #include "fix.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "group.h"
 #include "input.h"
 #include "lattice.h"
@@ -571,8 +571,8 @@ void ComputeChunkAtom::init()
 
   if ((idsflag == ONCE || lockcount) && !fixstore) {
     id_fix = utils::strdup(id + std::string("_COMPUTE_STORE"));
-    fixstore = dynamic_cast<FixStore *>(
-        modify->add_fix(fmt::format("{} {} STORE peratom 1 1", id_fix, group->names[igroup])));
+    fixstore = dynamic_cast<FixStorePeratom *>(
+      modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 1", id_fix, group->names[igroup])));
   }
 
   if ((idsflag != ONCE && !lockcount) && fixstore) {
diff --git a/src/compute_chunk_atom.h b/src/compute_chunk_atom.h
index a52d1c5445..2c1f3ef000 100644
--- a/src/compute_chunk_atom.h
+++ b/src/compute_chunk_atom.h
@@ -93,7 +93,7 @@ class ComputeChunkAtom : public Compute {
   double *varatom;
 
   char *id_fix;
-  class FixStore *fixstore;
+  class FixStorePeratom *fixstore;
 
   class Fix *lockfix;            // ptr to FixAveChunk that is locking out setups
                                  // null pointer if no lock currently in place
diff --git a/src/compute_displace_atom.cpp b/src/compute_displace_atom.cpp
index ba9b5c997f..6322bc1c00 100644
--- a/src/compute_displace_atom.cpp
+++ b/src/compute_displace_atom.cpp
@@ -17,7 +17,7 @@
 #include "atom.h"
 #include "domain.h"
 #include "error.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "group.h"
 #include "input.h"
 #include "memory.h"
@@ -74,8 +74,8 @@ ComputeDisplaceAtom::ComputeDisplaceAtom(LAMMPS *lmp, int narg, char **arg) :
   // id = compute-ID + COMPUTE_STORE, fix group = compute group
 
   id_fix = utils::strdup(std::string(id) + "_COMPUTE_STORE");
-  fix = dynamic_cast<FixStore *>( modify->add_fix(fmt::format("{} {} STORE peratom 1 3",
-                                                 id_fix, group->names[igroup])));
+  fix = dynamic_cast<FixStorePeratom *>(
+    modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 3", id_fix, group->names[igroup])));
 
   // calculate xu,yu,zu for fix store array
   // skip if reset from restart file
@@ -120,7 +120,7 @@ void ComputeDisplaceAtom::init()
 {
   // set fix which stores original atom coords
 
-  fix = dynamic_cast<FixStore *>( modify->get_fix_by_id(id_fix));
+  fix = dynamic_cast<FixStorePeratom *>(modify->get_fix_by_id(id_fix));
   if (!fix) error->all(FLERR,"Could not find compute displace/atom fix with ID {}", id_fix);
 
   if (refreshflag) {
diff --git a/src/compute_displace_atom.h b/src/compute_displace_atom.h
index 6ba0fb3d9f..b0e6710d8d 100644
--- a/src/compute_displace_atom.h
+++ b/src/compute_displace_atom.h
@@ -38,7 +38,7 @@ class ComputeDisplaceAtom : public Compute {
   int nmax;
   double **displace;
   char *id_fix;
-  class FixStore *fix;
+  class FixStorePeratom *fix;
 
   int refreshflag, ivar, nvmax;    // refresh option is enabled
   char *rvar;                      // for incremental dumps
diff --git a/src/compute_msd.cpp b/src/compute_msd.cpp
index 9fbaa603ec..7ecb89b4c0 100644
--- a/src/compute_msd.cpp
+++ b/src/compute_msd.cpp
@@ -16,7 +16,7 @@
 #include "atom.h"
 #include "domain.h"
 #include "error.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "group.h"
 #include "modify.h"
 #include "update.h"
@@ -27,7 +27,8 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-ComputeMSD::ComputeMSD(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), id_fix(nullptr)
+ComputeMSD::ComputeMSD(LAMMPS *lmp, int narg, char **arg)
+  : Compute(lmp, narg, arg), id_fix(nullptr)
 {
   if (narg < 3) error->all(FLERR, "Illegal compute msd command");
 
@@ -63,8 +64,8 @@ ComputeMSD::ComputeMSD(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, a
   // id = compute-ID + COMPUTE_STORE, fix group = compute group
 
   id_fix = utils::strdup(id + std::string("_COMPUTE_STORE"));
-  fix = dynamic_cast<FixStore *>(
-      modify->add_fix(fmt::format("{} {} STORE peratom 1 3", id_fix, group->names[igroup])));
+  fix = dynamic_cast<FixStorePeratom *>(
+      modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 3", id_fix, group->names[igroup])));
 
   // calculate xu,yu,zu for fix store array
   // skip if reset from restart file
@@ -127,7 +128,7 @@ void ComputeMSD::init()
 {
   // set fix which stores reference atom coords
 
-  fix = dynamic_cast<FixStore *>(modify->get_fix_by_id(id_fix));
+  fix = dynamic_cast<FixStorePeratom *>(modify->get_fix_by_id(id_fix));
   if (!fix) error->all(FLERR, "Could not find compute msd fix with ID {}", id_fix);
 
   // nmsd = # of atoms in group
diff --git a/src/compute_msd.h b/src/compute_msd.h
index c7ab7fa083..11ac804259 100644
--- a/src/compute_msd.h
+++ b/src/compute_msd.h
@@ -39,7 +39,7 @@ class ComputeMSD : public Compute {
   bigint nmsd;
   double masstotal;
   char *id_fix;
-  class FixStore *fix;
+  class FixStorePeratom *fix;
 };
 
 }    // namespace LAMMPS_NS
diff --git a/src/compute_msd_chunk.cpp b/src/compute_msd_chunk.cpp
index 8624520cf7..c2ec394db9 100644
--- a/src/compute_msd_chunk.cpp
+++ b/src/compute_msd_chunk.cpp
@@ -18,7 +18,7 @@
 #include "compute_chunk_atom.h"
 #include "domain.h"
 #include "error.h"
-#include "fix_store.h"
+#include "fix_store_global.h"
 #include "group.h"
 #include "memory.h"
 #include "modify.h"
@@ -31,8 +31,8 @@ using namespace LAMMPS_NS;
 /* ---------------------------------------------------------------------- */
 
 ComputeMSDChunk::ComputeMSDChunk(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg),
-  idchunk(nullptr), id_fix(nullptr), massproc(nullptr), masstotal(nullptr), com(nullptr), comall(nullptr), msd(nullptr)
+  Compute(lmp, narg, arg), idchunk(nullptr), id_fix(nullptr), massproc(nullptr),
+  masstotal(nullptr), com(nullptr), comall(nullptr), msd(nullptr)
 {
   if (narg != 4) error->all(FLERR,"Illegal compute msd/chunk command");
 
@@ -57,8 +57,8 @@ ComputeMSDChunk::ComputeMSDChunk(LAMMPS *lmp, int narg, char **arg) :
   // otherwise size reset and init will be done in setup()
 
   id_fix = utils::strdup(std::string(id) + "_COMPUTE_STORE");
-  fix = dynamic_cast<FixStore *>( modify->add_fix(fmt::format("{} {} STORE global 1 1",
-                                                 id_fix,group->names[igroup])));
+  fix = dynamic_cast<FixStoreGlobal *>(
+    modify->add_fix(fmt::format("{} {} STORE/GLOBAL 1 1", id_fix,group->names[igroup])));
 }
 
 /* ---------------------------------------------------------------------- */
@@ -93,7 +93,7 @@ void ComputeMSDChunk::init()
   // if firstflag, will be created in setup()
 
   if (!firstflag) {
-    fix = dynamic_cast<FixStore *>( modify->get_fix_by_id(id_fix));
+    fix = dynamic_cast<FixStoreGlobal *>(modify->get_fix_by_id(id_fix));
     if (!fix) error->all(FLERR,"Could not find compute msd/chunk fix with ID {}", id_fix);
   }
 }
diff --git a/src/compute_msd_chunk.h b/src/compute_msd_chunk.h
index 86fabd5efd..5e745bb39a 100644
--- a/src/compute_msd_chunk.h
+++ b/src/compute_msd_chunk.h
@@ -45,7 +45,7 @@ class ComputeMSDChunk : public Compute {
   char *idchunk;
   class ComputeChunkAtom *cchunk;
   char *id_fix;
-  class FixStore *fix;
+  class FixStoreGlobal *fix;
   int firstflag;
 
   double *massproc, *masstotal;
@@ -54,8 +54,6 @@ class ComputeMSDChunk : public Compute {
 
   void allocate();
 };
-
 }    // namespace LAMMPS_NS
-
 #endif
 #endif
diff --git a/src/compute_vacf.cpp b/src/compute_vacf.cpp
index ce3ba557aa..07669ec84f 100644
--- a/src/compute_vacf.cpp
+++ b/src/compute_vacf.cpp
@@ -18,7 +18,7 @@
 #include "update.h"
 #include "group.h"
 #include "modify.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
@@ -26,8 +26,7 @@ using namespace LAMMPS_NS;
 /* ---------------------------------------------------------------------- */
 
 ComputeVACF::ComputeVACF(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg),
-  id_fix(nullptr)
+  Compute(lmp, narg, arg), id_fix(nullptr)
 {
   if (narg < 3) error->all(FLERR,"Illegal compute vacf command");
 
@@ -40,7 +39,8 @@ ComputeVACF::ComputeVACF(LAMMPS *lmp, int narg, char **arg) :
   // id = compute-ID + COMPUTE_STORE, fix group = compute group
 
   id_fix = utils::strdup(id + std::string("_COMPUTE_STORE"));
-  fix = dynamic_cast<FixStore *>( modify->add_fix(fmt::format("{} {} STORE peratom 1 3", id_fix, group->names[igroup])));
+  fix = dynamic_cast<FixStorePeratom *>(
+    modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 3", id_fix, group->names[igroup])));
 
   // store current velocities in fix store array
   // skip if reset from restart file
@@ -84,7 +84,7 @@ void ComputeVACF::init()
 {
   // set fix which stores original atom velocities
 
-  fix = dynamic_cast<FixStore *>( modify->get_fix_by_id(id_fix));
+  fix = dynamic_cast<FixStorePeratom *>(modify->get_fix_by_id(id_fix));
   if (!fix) error->all(FLERR,"Could not find compute vacf fix ID {}", id_fix);
 
   // nvacf = # of atoms in group
diff --git a/src/compute_vacf.h b/src/compute_vacf.h
index 00ddb8be0f..b48bfcf475 100644
--- a/src/compute_vacf.h
+++ b/src/compute_vacf.h
@@ -35,7 +35,7 @@ class ComputeVACF : public Compute {
  protected:
   bigint nvacf;
   char *id_fix;
-  class FixStore *fix;
+  class FixStorePeratom *fix;
 };
 
 }    // namespace LAMMPS_NS
diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp
index 87bd04b2a9..9389839aca 100644
--- a/src/dump_custom.cpp
+++ b/src/dump_custom.cpp
@@ -20,7 +20,7 @@
 #include "domain.h"
 #include "error.h"
 #include "fix.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "group.h"
 #include "input.h"
 #include "memory.h"
@@ -2000,16 +2000,16 @@ int DumpCustom::modify_param(int narg, char **arg)
       thresh_value[nthresh] = utils::numeric(FLERR,arg[3],false,lmp);
       thresh_last[nthresh] = -1;
     } else {
-      thresh_fix = (FixStore **)
-        memory->srealloc(thresh_fix,(nthreshlast+1)*sizeof(FixStore *),"dump:thresh_fix");
+      thresh_fix = (FixStorePeratom **)
+        memory->srealloc(thresh_fix,(nthreshlast+1)*sizeof(FixStorePeratom *),"dump:thresh_fix");
       thresh_fixID = (char **)
         memory->srealloc(thresh_fixID,(nthreshlast+1)*sizeof(char *),"dump:thresh_fixID");
       memory->grow(thresh_first,(nthreshlast+1),"dump:thresh_first");
 
       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]);
-      thresh_fix[nthreshlast] = dynamic_cast<FixStore *>( modify->add_fix(threshid));
+      threshid += fmt::format(" {} STORE/PERATOM 1 1", group->names[igroup]);
+      thresh_fix[nthreshlast] = dynamic_cast<FixStorePeratom *>(modify->add_fix(threshid));
 
       thresh_last[nthreshlast] = nthreshlast;
       thresh_first[nthreshlast] = 1;
diff --git a/src/dump_custom.h b/src/dump_custom.h
index dd653d5e98..92c50fc173 100644
--- a/src/dump_custom.h
+++ b/src/dump_custom.h
@@ -37,19 +37,19 @@ class DumpCustom : public Dump {
   int nevery;        // dump frequency for output
   char *idregion;    // region ID, nullptr if no region
 
-  int nthresh;                    // # of defined thresholds
-  int nthreshlast;                // # of defined thresholds with value = LAST
-                                  //
-  int *thresh_array;              // array to threshold on for each nthresh
-  int *thresh_op;                 // threshold operation for each nthresh
-  double *thresh_value;           // threshold value for each nthresh
-  int *thresh_last;               // for threshold value = LAST,
-                                  // index into thresh_fix
-                                  // -1 if not LAST, value is numeric
-                                  //
-  class FixStore **thresh_fix;    // stores values for each threshold LAST
-  char **thresh_fixID;            // IDs of thresh_fixes
-  int *thresh_first;              // 1 the first time a FixStore values accessed
+  int nthresh;                           // # of defined thresholds
+  int nthreshlast;                       // # of defined thresholds with value = LAST
+                                         //
+  int *thresh_array;                     // array to threshold on for each nthresh
+  int *thresh_op;                        // threshold operation for each nthresh
+  double *thresh_value;                  // threshold value for each nthresh
+  int *thresh_last;                      // for threshold value = LAST,
+                                         // index into thresh_fix
+                                         // -1 if not LAST, value is numeric
+                                         //
+  class FixStorePeratom **thresh_fix;    // stores values for each threshold LAST
+  char **thresh_fixID;                   // IDs of thresh_fixes
+  int *thresh_first;                     // 1 the first time a FixStore values accessed
 
   int expand;        // flag for whether field args were expanded
   char **earg;       // field names with wildcard expansion
diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp
index c54e569c3a..3fa6f4b8df 100644
--- a/src/fix_adapt.cpp
+++ b/src/fix_adapt.cpp
@@ -19,7 +19,7 @@
 #include "bond.h"
 #include "domain.h"
 #include "error.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "force.h"
 #include "group.h"
 #include "input.h"
@@ -278,8 +278,8 @@ void FixAdapt::post_constructor()
 
   if (diamflag && atom->radius_flag) {
     id_fix_diam = utils::strdup(id + std::string("_FIX_STORE_DIAM"));
-    fix_diam = dynamic_cast<FixStore *>( modify->add_fix(fmt::format("{} {} STORE peratom 1 1",
-                                                        id_fix_diam,group->names[igroup])));
+    fix_diam = dynamic_cast<FixStorePeratom *>(
+      modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 1", id_fix_diam,group->names[igroup])));
     if (fix_diam->restart_reset) fix_diam->restart_reset = 0;
     else {
       double *vec = fix_diam->vstore;
@@ -296,8 +296,8 @@ void FixAdapt::post_constructor()
 
   if (chgflag && atom->q_flag) {
     id_fix_chg = utils::strdup(id + std::string("_FIX_STORE_CHG"));
-    fix_chg = dynamic_cast<FixStore *>( modify->add_fix(fmt::format("{} {} STORE peratom 1 1",
-                                                       id_fix_chg,group->names[igroup])));
+    fix_chg = dynamic_cast<FixStorePeratom *>(
+      modify->add_fix(fmt::format("{} {} STORE/PERATOM 1 1",id_fix_chg,group->names[igroup])));
     if (fix_chg->restart_reset) fix_chg->restart_reset = 0;
     else {
       double *vec = fix_chg->vstore;
@@ -492,18 +492,16 @@ void FixAdapt::init()
   // fixes that store initial per-atom values
 
   if (id_fix_diam) {
-    int ifix = modify->find_fix(id_fix_diam);
-    if (ifix < 0) error->all(FLERR,"Could not find fix adapt storage fix ID");
-    fix_diam = dynamic_cast<FixStore *>( modify->fix[ifix]);
+    fix_diam = dynamic_cast<FixStorePeratom *>(modify->get_fix_by_id(id_fix_diam));
+    if (!fix_diam) error->all(FLERR,"Could not find fix adapt storage fix ID {}", id_fix_diam);
   }
   if (id_fix_chg) {
-    int ifix = modify->find_fix(id_fix_chg);
-    if (ifix < 0) error->all(FLERR,"Could not find fix adapt storage fix ID");
-    fix_chg = dynamic_cast<FixStore *>( modify->fix[ifix]);
+    fix_chg = dynamic_cast<FixStorePeratom *>(modify->get_fix_by_id(id_fix_chg));
+    if (!fix_chg) error->all(FLERR,"Could not find fix adapt storage fix ID {}", id_fix_chg);
   }
 
   if (utils::strmatch(update->integrate_style,"^respa"))
-    nlevels_respa = (dynamic_cast<Respa *>( update->integrate))->nlevels;
+    nlevels_respa = (dynamic_cast<Respa *>(update->integrate))->nlevels;
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/fix_adapt.h b/src/fix_adapt.h
index cdd90fc225..d4db692895 100644
--- a/src/fix_adapt.h
+++ b/src/fix_adapt.h
@@ -48,7 +48,7 @@ class FixAdapt : public Fix {
   int anypair, anybond, anyangle;
   int nlevels_respa;
   char *id_fix_diam, *id_fix_chg;
-  class FixStore *fix_diam, *fix_chg;
+  class FixStorePeratom *fix_diam, *fix_chg;
   double previous_diam_scale, previous_chg_scale;
   int discflag;
 
diff --git a/src/fix_balance.cpp b/src/fix_balance.cpp
index 6b5c4b6ab1..1ce3c212d4 100644
--- a/src/fix_balance.cpp
+++ b/src/fix_balance.cpp
@@ -19,7 +19,7 @@
 #include "comm.h"
 #include "domain.h"
 #include "error.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "force.h"
 #include "irregular.h"
 #include "kspace.h"
diff --git a/src/variable.cpp b/src/variable.cpp
index 84ab0d69ae..168ccfa1cd 100644
--- a/src/variable.cpp
+++ b/src/variable.cpp
@@ -20,7 +20,7 @@
 #include "domain.h"
 #include "error.h"
 #include "fix.h"
-#include "fix_store.h"
+#include "fix_store_peratom.h"
 #include "group.h"
 #include "info.h"
 #include "input.h"
@@ -4911,7 +4911,8 @@ VarReader::VarReader(LAMMPS *lmp, char *name, char *file, int flag) :
       error->all(FLERR,"Cannot use atomfile-style variable unless an atom map exists");
 
     id_fix = utils::strdup(std::string(name) + "_VARIABLE_STORE");
-    fixstore = dynamic_cast<FixStore *>( modify->add_fix(std::string(id_fix) + " all STORE peratom 0 1"));
+    fixstore = dynamic_cast<FixStorePeratom *>(
+      modify->add_fix(std::string(id_fix) + " all STORE/PERATOM 0 1"));
     buffer = new char[CHUNK*MAXLINE];
   }
 }
diff --git a/src/variable.h b/src/variable.h
index 52632ba6ce..1e05fd3904 100644
--- a/src/variable.h
+++ b/src/variable.h
@@ -151,7 +151,7 @@ class Variable : protected Pointers {
 
 class VarReader : protected Pointers {
  public:
-  class FixStore *fixstore;
+  class FixStorePeratom *fixstore;
   char *id_fix;
 
   VarReader(class LAMMPS *, char *, char *, int);