diff --git a/cmake/Modules/Packages/MACHDYN.cmake b/cmake/Modules/Packages/MACHDYN.cmake
index 6d941f9798..fab532541e 100644
--- a/cmake/Modules/Packages/MACHDYN.cmake
+++ b/cmake/Modules/Packages/MACHDYN.cmake
@@ -7,8 +7,8 @@ endif()
 option(DOWNLOAD_EIGEN3 "Download Eigen3 instead of using an already installed one)" ${DOWNLOAD_EIGEN3_DEFAULT})
 if(DOWNLOAD_EIGEN3)
   message(STATUS "Eigen3 download requested - we will build our own")
-  set(EIGEN3_URL "https://gitlab.com/libeigen/eigen/-/archive/3.3.7/eigen-3.3.7.tar.gz" CACHE STRING "URL for Eigen3 tarball")
-  set(EIGEN3_MD5 "9e30f67e8531477de4117506fe44669b" CACHE STRING "MD5 checksum of Eigen3 tarball")
+  set(EIGEN3_URL "https://gitlab.com/libeigen/eigen/-/archive/3.3.9/eigen-3.3.9.tar.gz" CACHE STRING "URL for Eigen3 tarball")
+  set(EIGEN3_MD5 "609286804b0f79be622ccf7f9ff2b660" CACHE STRING "MD5 checksum of Eigen3 tarball")
   mark_as_advanced(EIGEN3_URL)
   mark_as_advanced(EIGEN3_MD5)
   include(ExternalProject)
@@ -30,3 +30,8 @@ else()
   endif()
   target_link_libraries(lammps PRIVATE Eigen3::Eigen)
 endif()
+
+# PGI/Nvidia compiler internals collide with vector intrinsics support in Eigen3
+if((CMAKE_CXX_COMPILER_ID STREQUAL "PGI") OR (CMAKE_CXX_COMPILER_ID STREQUAL "NVHPC"))
+  target_compile_definitions(lammps PRIVATE -DEIGEN_DONT_VECTORIZE)
+endif()
diff --git a/lib/machdyn/Install.py b/lib/machdyn/Install.py
index 16dd3038a8..2e90c9ca0f 100644
--- a/lib/machdyn/Install.py
+++ b/lib/machdyn/Install.py
@@ -17,11 +17,12 @@ parser = ArgumentParser(prog='Install.py',
 
 # settings
 
-version = '3.3.7'
+version = '3.3.9'
 tarball = "eigen.tar.gz"
 
 # known checksums for different Eigen versions. used to validate the download.
 checksums = { \
+              '3.3.9' : '609286804b0f79be622ccf7f9ff2b660', \
               '3.3.7' : '9e30f67e8531477de4117506fe44669b' \
 }
 
diff --git a/lib/machdyn/Makefile.lammps b/lib/machdyn/Makefile.lammps
index 6c914708dd..71db4fd4bc 100644
--- a/lib/machdyn/Makefile.lammps
+++ b/lib/machdyn/Makefile.lammps
@@ -1,5 +1,5 @@
 # Settings that the LAMMPS build will import when this package library is used
 
-machdyn_SYSINC = -I../../lib/includelink/eigen3
-machdyn_SYSLIB = 
-machdyn_SYSPATH = 
+machdyn_SYSINC =
+machdyn_SYSLIB =
+machdyn_SYSPATH =
diff --git a/src/CORESHELL/compute_temp_cs.cpp b/src/CORESHELL/compute_temp_cs.cpp
index ed6634f4b0..04c7a743d8 100644
--- a/src/CORESHELL/compute_temp_cs.cpp
+++ b/src/CORESHELL/compute_temp_cs.cpp
@@ -19,19 +19,19 @@
 
 #include "compute_temp_cs.h"
 
-#include <cstring>
-
 #include "atom.h"
 #include "atom_vec.h"
+#include "comm.h"
 #include "domain.h"
-#include "update.h"
+#include "error.h"
+#include "fix_store.h"
 #include "force.h"
 #include "group.h"
-#include "modify.h"
-#include "fix_store.h"
-#include "comm.h"
 #include "memory.h"
-#include "error.h"
+#include "modify.h"
+#include "update.h"
+
+#include <cstring>
 
 
 using namespace LAMMPS_NS;
@@ -74,8 +74,7 @@ ComputeTempCS::ComputeTempCS(LAMMPS *lmp, int narg, char **arg) :
   strcpy(id_fix,fixcmd.c_str());
 
   fixcmd += fmt::format(" {} STORE peratom 0 1", group->names[igroup]);
-  modify->add_fix(fixcmd);
-  fix = (FixStore *) modify->fix[modify->nfix-1];
+  fix = (FixStore *)modify->add_fix(fixcmd);
 
   // set fix store values = 0 for now
   // fill them in via setup() once Comm::borders() has been called
diff --git a/src/DPD-REACT/fix_rx.cpp b/src/DPD-REACT/fix_rx.cpp
index 52dcaaf2d5..e50e145853 100644
--- a/src/DPD-REACT/fix_rx.cpp
+++ b/src/DPD-REACT/fix_rx.cpp
@@ -18,6 +18,7 @@
 #include "comm.h"
 #include "domain.h"
 #include "error.h"
+#include "fix_property_atom.h"
 #include "force.h"
 #include "group.h"
 #include "math_special.h"
@@ -344,12 +345,9 @@ void FixRX::post_constructor()
   newcmd1 += " ghost yes";
   newcmd2 += " ghost yes";
 
-  modify->add_fix(newcmd1);
-  fix_species = (FixPropertyAtom *) modify->fix[modify->nfix-1];
-  restartFlag = modify->fix[modify->nfix-1]->restart_reset;
-
-  modify->add_fix(newcmd2);
-  fix_species_old = (FixPropertyAtom *) modify->fix[modify->nfix-1];
+  fix_species = (FixPropertyAtom *) modify->add_fix(newcmd1);
+  restartFlag = fix_species->restart_reset;
+  fix_species_old = (FixPropertyAtom *) modify->add_fix(newcmd2);
 
   if (nspecies==0) error->all(FLERR,"There are no rx species specified.");
 
diff --git a/src/DPD-REACT/fix_rx.h b/src/DPD-REACT/fix_rx.h
index debf4d911e..1d65c4c09e 100644
--- a/src/DPD-REACT/fix_rx.h
+++ b/src/DPD-REACT/fix_rx.h
@@ -136,7 +136,7 @@ class FixRX : public Fix {
   //!< ODE Solver diagnostics.
   void odeDiagnostics(void);
 
- private:
+ protected:
   char *kineticsFile;
   char *id_fix_species, *id_fix_species_old;
   class FixPropertyAtom *fix_species, *fix_species_old;
diff --git a/src/EXTRA-COMPUTE/compute_hma.cpp b/src/EXTRA-COMPUTE/compute_hma.cpp
index 09a2840906..043dcec6cf 100644
--- a/src/EXTRA-COMPUTE/compute_hma.cpp
+++ b/src/EXTRA-COMPUTE/compute_hma.cpp
@@ -91,9 +91,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");
-  modify->add_fix(fmt::format("{} {} STORE peratom 1 3",
-                              id_fix, group->names[igroup]));
-  fix = (FixStore *) modify->fix[modify->nfix-1];
+  fix = (FixStore *)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
diff --git a/src/FEP/fix_adapt_fep.cpp b/src/FEP/fix_adapt_fep.cpp
index b5aad03887..3b4290f2c9 100644
--- a/src/FEP/fix_adapt_fep.cpp
+++ b/src/FEP/fix_adapt_fep.cpp
@@ -216,10 +216,8 @@ void FixAdaptFEP::post_constructor()
   id_fix_chg = nullptr;
 
   if (diamflag) {
-    auto cmd = fmt::format("{}_FIX_STORE_DIAM {} STORE peratom 1 1",
-                           group->names[igroup]);
-    modify->add_fix(cmd);
-    fix_diam = (FixStore *) modify->fix[modify->nfix-1];
+    auto cmd = fmt::format("{}_FIX_STORE_DIAM {} STORE peratom 1 1", group->names[igroup]);
+    fix_diam = (FixStore *) modify->add_fix(cmd);
 
     if (fix_diam->restart_reset) fix_diam->restart_reset = 0;
     else {
@@ -236,10 +234,8 @@ void FixAdaptFEP::post_constructor()
   }
 
   if (chgflag) {
-    auto cmd = fmt::format("{}_FIX_STORE_CHG {} STORE peratom 1 1",
-                           group->names[igroup]);
-    modify->add_fix(cmd);
-    fix_chg = (FixStore *) modify->fix[modify->nfix-1];
+    auto cmd = fmt::format("{}_FIX_STORE_CHG {} STORE peratom 1 1", group->names[igroup]);
+    fix_chg = (FixStore *) modify->add_fix(cmd);
 
     if (fix_chg->restart_reset) fix_chg->restart_reset = 0;
     else {
diff --git a/src/GPU/Install.sh b/src/GPU/Install.sh
index 2e2b32d96a..a87d2165d9 100755
--- a/src/GPU/Install.sh
+++ b/src/GPU/Install.sh
@@ -113,8 +113,8 @@ action pair_lj_cut_coul_msm_gpu.cpp pair_lj_cut_coul_msm.cpp
 action pair_lj_cut_coul_msm_gpu.h pair_lj_cut_coul_msm.h
 action pair_lj_cut_gpu.cpp
 action pair_lj_cut_gpu.h
-action pair_lj_smooth_gpu.cpp
-action pair_lj_smooth_gpu.h
+action pair_lj_smooth_gpu.cpp pair_lj_smooth.cpp
+action pair_lj_smooth_gpu.h pair_lj_smooth.h
 action pair_lj_expand_gpu.cpp
 action pair_lj_expand_gpu.h
 action pair_lj_expand_coul_long_gpu.cpp pair_lj_expand_coul_long.cpp
diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp
index 7137b911f1..e6013b9940 100644
--- a/src/GRANULAR/pair_gran_hooke_history.cpp
+++ b/src/GRANULAR/pair_gran_hooke_history.cpp
@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -68,8 +67,8 @@ PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp)
   // this is so final order of Modify:fix will conform to input script
 
   fix_history = nullptr;
-  modify->add_fix("NEIGH_HISTORY_HH_DUMMY"+std::to_string(instance_me)+" all DUMMY");
-  fix_dummy = (FixDummy *) modify->fix[modify->nfix-1];
+  fix_dummy = (FixDummy *) modify->add_fix("NEIGH_HISTORY_HH_DUMMY"
+                                           + std::to_string(instance_me) + " all DUMMY");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -437,20 +436,19 @@ void PairGranHookeHistory::init_style()
 
   if (history && (fix_history == nullptr)) {
     auto cmd = fmt::format("NEIGH_HISTORY_HH{} all NEIGH_HISTORY {}", instance_me, size_history);
-    modify->replace_fix("NEIGH_HISTORY_HH_DUMMY"+std::to_string(instance_me),cmd,1);
-    int ifix = modify->find_fix("NEIGH_HISTORY_HH"+std::to_string(instance_me));
-    fix_history = (FixNeighHistory *) modify->fix[ifix];
+    fix_history = (FixNeighHistory *) modify->replace_fix("NEIGH_HISTORY_HH_DUMMY"
+                                                          + std::to_string(instance_me),cmd,1);
     fix_history->pair = this;
   }
 
   // check for FixFreeze and set freeze_group_bit
 
-  for (i = 0; i < modify->nfix; i++)
-    if (strcmp(modify->fix[i]->style,"freeze") == 0) break;
-  if (i < modify->nfix) freeze_group_bit = modify->fix[i]->groupbit;
-  else freeze_group_bit = 0;
+  int ifreeze = modify->find_fix_by_style("^freeze");
+  if (ifreeze < 0) freeze_group_bit = 0;
+  else freeze_group_bit = modify->fix[ifreeze]->groupbit;
 
   // check for FixRigid so can extract rigid body masses
+  // FIXME: this only catches the first rigid fix, there may be multiple.
 
   fix_rigid = nullptr;
   for (i = 0; i < modify->nfix; i++)
@@ -459,15 +457,8 @@ void PairGranHookeHistory::init_style()
 
   // check for FixPour and FixDeposit so can extract particle radii
 
-  int ipour;
-  for (ipour = 0; ipour < modify->nfix; ipour++)
-    if (strcmp(modify->fix[ipour]->style,"pour") == 0) break;
-  if (ipour == modify->nfix) ipour = -1;
-
-  int idep;
-  for (idep = 0; idep < modify->nfix; idep++)
-    if (strcmp(modify->fix[idep]->style,"deposit") == 0) break;
-  if (idep == modify->nfix) idep = -1;
+  int ipour = modify->find_fix_by_style("^pour");
+  int idep  = modify->find_fix_by_style("^deposit");
 
   // set maxrad_dynamic and maxrad_frozen for each type
   // include future FixPour and FixDeposit particles as dynamic
@@ -498,10 +489,8 @@ void PairGranHookeHistory::init_style()
     else
       onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]],radius[i]);
 
-  MPI_Allreduce(&onerad_dynamic[1],&maxrad_dynamic[1],atom->ntypes,
-                MPI_DOUBLE,MPI_MAX,world);
-  MPI_Allreduce(&onerad_frozen[1],&maxrad_frozen[1],atom->ntypes,
-                MPI_DOUBLE,MPI_MAX,world);
+  MPI_Allreduce(&onerad_dynamic[1],&maxrad_dynamic[1],atom->ntypes,MPI_DOUBLE,MPI_MAX,world);
+  MPI_Allreduce(&onerad_frozen[1],&maxrad_frozen[1],atom->ntypes,MPI_DOUBLE,MPI_MAX,world);
 
   // set fix which stores history info
 
diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp
index cca349b819..d0e75d4812 100644
--- a/src/GRANULAR/pair_granular.cpp
+++ b/src/GRANULAR/pair_granular.cpp
@@ -20,25 +20,24 @@
 
 #include "pair_granular.h"
 
-#include <cmath>
-#include <cstring>
-
 #include "atom.h"
-#include "force.h"
-#include "update.h"
-#include "modify.h"
+#include "comm.h"
+#include "error.h"
 #include "fix.h"
 #include "fix_dummy.h"
 #include "fix_neigh_history.h"
-#include "comm.h"
-#include "neighbor.h"
-#include "neigh_list.h"
-#include "neigh_request.h"
-#include "memory.h"
-#include "error.h"
+#include "force.h"
 #include "math_const.h"
 #include "math_special.h"
+#include "memory.h"
+#include "modify.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "neighbor.h"
+#include "update.h"
 
+#include <cmath>
+#include <cstring>
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
@@ -103,8 +102,7 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp)
   // this is so final order of Modify:fix will conform to input script
 
   fix_history = nullptr;
-  modify->add_fix("NEIGH_HISTORY_GRANULAR_DUMMY all DUMMY");
-  fix_dummy = (FixDummy *) modify->fix[modify->nfix-1];
+  fix_dummy = (FixDummy *) modify->add_fix("NEIGH_HISTORY_GRANULAR_DUMMY all DUMMY");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1122,21 +1120,21 @@ void PairGranular::init_style()
   // this is so its order in the fix list is preserved
 
   if (use_history && fix_history == nullptr) {
-    modify->replace_fix("NEIGH_HISTORY_GRANULAR_DUMMY","NEIGH_HISTORY_GRANULAR"
-                        " all NEIGH_HISTORY " + std::to_string(size_history),1);
-    int ifix = modify->find_fix("NEIGH_HISTORY_GRANULAR");
-    fix_history = (FixNeighHistory *) modify->fix[ifix];
+    fix_history = (FixNeighHistory *) modify->replace_fix("NEIGH_HISTORY_GRANULAR_DUMMY",
+                                                          "NEIGH_HISTORY_GRANULAR"
+                                                          " all NEIGH_HISTORY "
+                                                          + std::to_string(size_history),1);
     fix_history->pair = this;
   }
 
   // check for FixFreeze and set freeze_group_bit
 
-  for (i = 0; i < modify->nfix; i++)
-    if (strcmp(modify->fix[i]->style,"freeze") == 0) break;
-  if (i < modify->nfix) freeze_group_bit = modify->fix[i]->groupbit;
-  else freeze_group_bit = 0;
+  int ifix = modify->find_fix_by_style("^freeze");
+  if (ifix < 0) freeze_group_bit = 0;
+  else freeze_group_bit = modify->fix[ifix]->groupbit;
 
   // check for FixRigid so can extract rigid body masses
+  // FIXME: this only catches the first rigid fix, there may be multiple.
 
   fix_rigid = nullptr;
   for (i = 0; i < modify->nfix; i++)
@@ -1145,15 +1143,8 @@ void PairGranular::init_style()
 
   // check for FixPour and FixDeposit so can extract particle radii
 
-  int ipour;
-  for (ipour = 0; ipour < modify->nfix; ipour++)
-    if (strcmp(modify->fix[ipour]->style,"pour") == 0) break;
-  if (ipour == modify->nfix) ipour = -1;
-
-  int idep;
-  for (idep = 0; idep < modify->nfix; idep++)
-    if (strcmp(modify->fix[idep]->style,"deposit") == 0) break;
-  if (idep == modify->nfix) idep = -1;
+  int ipour = modify->find_fix_by_style("^pour");
+  int idep  = modify->find_fix_by_style("^deposit");
 
   // set maxrad_dynamic and maxrad_frozen for each type
   // include future FixPour and FixDeposit particles as dynamic
@@ -1187,10 +1178,8 @@ void PairGranular::init_style()
     }
   }
 
-  MPI_Allreduce(&onerad_dynamic[1],&maxrad_dynamic[1],atom->ntypes,
-                MPI_DOUBLE,MPI_MAX,world);
-  MPI_Allreduce(&onerad_frozen[1],&maxrad_frozen[1],atom->ntypes,
-                MPI_DOUBLE,MPI_MAX,world);
+  MPI_Allreduce(&onerad_dynamic[1],&maxrad_dynamic[1],atom->ntypes,MPI_DOUBLE,MPI_MAX,world);
+  MPI_Allreduce(&onerad_frozen[1],&maxrad_frozen[1],atom->ntypes,MPI_DOUBLE,MPI_MAX,world);
 
   // set fix which stores history info
 
diff --git a/src/KOKKOS/fix_rx_kokkos.cpp b/src/KOKKOS/fix_rx_kokkos.cpp
index 7bc7bf95c3..15b4a39849 100644
--- a/src/KOKKOS/fix_rx_kokkos.cpp
+++ b/src/KOKKOS/fix_rx_kokkos.cpp
@@ -13,24 +13,25 @@
 ------------------------------------------------------------------------- */
 
 #include "fix_rx_kokkos.h"
-#include <cstring>
-#include "atom_masks.h"
+
 #include "atom_kokkos.h"
-#include "force.h"
-#include "memory_kokkos.h"
-#include "update.h"
-#include "modify.h"
-#include "neighbor.h"
-#include "neigh_list_kokkos.h"
-#include "neigh_request.h"
-#include "error.h"
-#include "math_special_kokkos.h"
+#include "atom_masks.h"
 #include "comm.h"
 #include "domain.h"
+#include "error.h"
+#include "fix_property_atom.h"
+#include "force.h"
 #include "kokkos.h"
-
+#include "math_special_kokkos.h"
+#include "memory_kokkos.h"
+#include "modify.h"
+#include "neigh_list_kokkos.h"
+#include "neigh_request.h"
+#include "neighbor.h"
+#include "update.h"
 
 #include <cfloat> // DBL_EPSILON
+#include <cstring>
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
@@ -104,7 +105,7 @@ void FixRxKokkos<DeviceType>::post_constructor()
   FixRX::post_constructor();
 
   // Need a copy of this
-  this->my_restartFlag = modify->fix[modify->nfix-1]->restart_reset;
+  this->my_restartFlag = fix_species->restart_reset;
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp b/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp
index 93f8e37006..fab33d0ec7 100644
--- a/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp
+++ b/src/KOKKOS/pair_gran_hooke_history_kokkos.cpp
@@ -13,17 +13,18 @@
 ------------------------------------------------------------------------- */
 
 #include "pair_gran_hooke_history_kokkos.h"
-#include "kokkos.h"
+
 #include "atom_kokkos.h"
 #include "atom_masks.h"
-#include "memory_kokkos.h"
+#include "error.h"
+#include "fix_neigh_history_kokkos.h"
 #include "force.h"
-#include "neighbor.h"
+#include "kokkos.h"
+#include "memory_kokkos.h"
+#include "modify.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
-#include "error.h"
-#include "modify.h"
-#include "fix_neigh_history_kokkos.h"
+#include "neighbor.h"
 #include "update.h"
 
 using namespace LAMMPS_NS;
@@ -67,20 +68,13 @@ void PairGranHookeHistoryKokkos<DeviceType>::init_style()
   // this is so its order in the fix list is preserved
 
   if (history && fix_history == nullptr) {
-    char dnumstr[16];
-    sprintf(dnumstr,"%d",3);
-    char **fixarg = new char*[4];
-    fixarg[0] = (char *) "NEIGH_HISTORY_HH";
-    fixarg[1] = (char *) "all";
+    auto cmd = std::string("NEIGH_HISTORY_HH") + std::to_string(instance_me) + " all ";
     if (execution_space == Device)
-      fixarg[2] = (char *) "NEIGH_HISTORY/KK/DEVICE";
+      cmd += "NEIGH_HISTORY/KK/DEVICE 3";
     else
-      fixarg[2] = (char *) "NEIGH_HISTORY/KK/HOST";
-    fixarg[3] = dnumstr;
-    modify->replace_fix("NEIGH_HISTORY_HH_DUMMY",4,fixarg,1);
-    delete [] fixarg;
-    int ifix = modify->find_fix("NEIGH_HISTORY_HH");
-    fix_history = (FixNeighHistory *) modify->fix[ifix];
+      cmd += "NEIGH_HISTORY/KK/HOST 3";
+    fix_history = (FixNeighHistory *)
+      modify->replace_fix("NEIGH_HISTORY_HH_DUMMY"+std::to_string(instance_me),cmd,1);
     fix_history->pair = this;
     fix_historyKK = (FixNeighHistoryKokkos<DeviceType> *)fix_history;
   }
diff --git a/src/MISC/pair_srp.cpp b/src/MISC/pair_srp.cpp
index 332e103b12..0fef989c88 100644
--- a/src/MISC/pair_srp.cpp
+++ b/src/MISC/pair_srp.cpp
@@ -84,8 +84,7 @@ PairSRP::PairSRP(LAMMPS *lmp) : Pair(lmp), fix_id(nullptr)
   //   will be invoked before other fixes that migrate atoms
   //   this is checked for in FixSRP
 
-  modify->add_fix(fmt::format("{:02d}_FIX_SRP all SRP",srp_instance));
-  f_srp = (FixSRP *) modify->fix[modify->nfix-1];
+  f_srp = (FixSRP *) modify->add_fix(fmt::format("{:02d}_FIX_SRP all SRP",srp_instance));
   ++srp_instance;
 }
 
diff --git a/src/OPENMP/pair_reaxff_omp.cpp b/src/OPENMP/pair_reaxff_omp.cpp
index 69d68e67e8..6cf3e90575 100644
--- a/src/OPENMP/pair_reaxff_omp.cpp
+++ b/src/OPENMP/pair_reaxff_omp.cpp
@@ -136,10 +136,8 @@ void PairReaxFFOMP::init_style()
     error->warning(FLERR,"Total cutoff < 2*bond cutoff. May need to use an "
                    "increased neighbor list skin.");
 
-  if (fix_reaxff == nullptr) {
-    modify->add_fix(fmt::format("{} all REAXFF",fix_id));
-    fix_reaxff = (FixReaxFF *) modify->fix[modify->nfix-1];
-  }
+  if (fix_reaxff == nullptr)
+    fix_reaxff = (FixReaxFF *) modify->add_fix(fmt::format("{} all REAXFF",fix_id));
 
   api->control->nthreads = comm->nthreads;
 }
diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp
index 46d7efdf8d..e52347727f 100644
--- a/src/REACTION/fix_bond_react.cpp
+++ b/src/REACTION/fix_bond_react.cpp
@@ -701,37 +701,25 @@ it will have the name 'i_limit_tags' and will be intitialized to 0 (not in group
 void FixBondReact::post_constructor()
 {
   // let's add the limit_tags per-atom property fix
-  std::string cmd = std::string("bond_react_props_internal");
-  id_fix2 = new char[cmd.size()+1];
-  strcpy(id_fix2,cmd.c_str());
-
-  int ifix = modify->find_fix(id_fix2);
-  if (ifix == -1) {
-    cmd += std::string(" all property/atom i_limit_tags i_react_tags ghost yes");
-    modify->add_fix(cmd);
-  }
+  id_fix2 = utils::strdup("bond_react_props_internal");
+  if (modify->find_fix(id_fix2) < 0)
+    modify->add_fix(std::string(id_fix2)+" all property/atom i_limit_tags i_react_tags ghost yes");
 
   // create master_group if not already existing
   // NOTE: limit_tags and react_tags automaticaly intitialized to zero (unless read from restart)
   group->find_or_create(master_group);
-  cmd = fmt::format("{} dynamic all property limit_tags",master_group);
+  std::string cmd = fmt::format("{} dynamic all property limit_tags",master_group);
   group->assign(cmd);
 
   if (stabilization_flag == 1) {
     int groupid = group->find(exclude_group);
     // create exclude_group if not already existing, or use as parent group if static
     if (groupid == -1 || group->dynamic[groupid] == 0) {
-      // create stabilization per-atom property
-      cmd = std::string("bond_react_stabilization_internal");
-      id_fix3 = new char[cmd.size()+1];
-      strcpy(id_fix3,cmd.c_str());
 
-      ifix = modify->find_fix(id_fix3);
-      if (ifix == -1) {
-        cmd += std::string(" all property/atom i_statted_tags ghost yes");
-        modify->add_fix(cmd);
-        fix3 = modify->fix[modify->nfix-1];
-      }
+      // create stabilization per-atom property
+      id_fix3 = utils::strdup("bond_react_stabilization_internal");
+      if (modify->find_fix(id_fix3) < 0)
+        fix3 = modify->add_fix(std::string(id_fix3) + " all property/atom i_statted_tags ghost yes");
 
       statted_id = utils::strdup("statted_tags");
 
@@ -770,7 +758,7 @@ void FixBondReact::post_constructor()
 
       // this returns names of corresponding property
       int unused;
-      char * idprop;
+      char *idprop;
       idprop = (char *) fix->extract("property",unused);
       if (idprop == nullptr)
         error->all(FLERR,"Exclude group must be a per-atom property group");
@@ -787,17 +775,10 @@ void FixBondReact::post_constructor()
 
 
     // let's create a new nve/limit fix to limit newly reacted atoms
-    cmd = std::string("bond_react_MASTER_nve_limit");
-    id_fix1 = new char[cmd.size()+1];
-    strcpy(id_fix1,cmd.c_str());
-
-    ifix = modify->find_fix(id_fix1);
-
-    if (ifix == -1) {
-      cmd += fmt::format(" {} nve/limit  {}",master_group,nve_limit_xmax);
-      modify->add_fix(cmd);
-      fix1 = modify->fix[modify->nfix-1];
-    }
+    id_fix1 = utils::strdup("bond_react_MASTER_nve_limit");
+    if (modify->find_fix(id_fix1) < 0)
+      fix1 = modify->add_fix(fmt::format("{} {} nve/limit  {}",
+                                         id_fix1,master_group,nve_limit_xmax));
   }
 }
 
diff --git a/src/REAXFF/fix_reaxff_species.cpp b/src/REAXFF/fix_reaxff_species.cpp
index 53746b968b..4fc08975a8 100644
--- a/src/REAXFF/fix_reaxff_species.cpp
+++ b/src/REAXFF/fix_reaxff_species.cpp
@@ -309,113 +309,16 @@ void FixReaxFFSpecies::init()
 
   if (!setupflag) {
     // create a compute to store properties
-    create_compute();
+    modify->add_compute("SPECATOM all SPEC/ATOM q x y z vx vy vz abo01 abo02 abo03 abo04 "
+                        "abo05 abo06 abo07 abo08 abo09 abo10 abo11 abo12 abo13 abo14 "
+                        "abo15 abo16 abo17 abo18 abo19 abo20 abo21 abo22 abo23 abo24");
 
     // create a fix to point to fix_ave_atom for averaging stored properties
-    create_fix();
-
+    auto fixcmd = fmt::format("SPECBOND all ave/atom {} {} {}",tmparg[0],tmparg[1],tmparg[2]);
+    for (int i = 1; i < 32; ++i) fixcmd += " c_SPECATOM[" + std::to_string(i) + "]";
+    f_SPECBOND = (FixAveAtom *) modify->add_fix(fixcmd);
     setupflag = 1;
   }
-
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixReaxFFSpecies::create_compute()
-{
-  int narg;
-  char **args;
-
-  narg = 34;
-  args = new char*[narg];
-  args[0]  = (char *) "SPECATOM";
-  args[1]  = (char *) "all";
-  args[2]  = (char *) "SPEC/ATOM";
-  args[3]  = (char *) "q";
-  args[4]  = (char *) "x";
-  args[5]  = (char *) "y";
-  args[6]  = (char *) "z";
-  args[7]  = (char *) "vx";
-  args[8]  = (char *) "vy";
-  args[9]  = (char *) "vz";
-  args[10] = (char *) "abo01";
-  args[11] = (char *) "abo02";
-  args[12] = (char *) "abo03";
-  args[13] = (char *) "abo04";
-  args[14] = (char *) "abo05";
-  args[15] = (char *) "abo06";
-  args[16] = (char *) "abo07";
-  args[17] = (char *) "abo08";
-  args[18] = (char *) "abo09";
-  args[19] = (char *) "abo10";
-  args[20] = (char *) "abo11";
-  args[21] = (char *) "abo12";
-  args[22] = (char *) "abo13";
-  args[23] = (char *) "abo14";
-  args[24] = (char *) "abo15";
-  args[25] = (char *) "abo16";
-  args[26] = (char *) "abo17";
-  args[27] = (char *) "abo18";
-  args[28] = (char *) "abo19";
-  args[29] = (char *) "abo20";
-  args[30] = (char *) "abo21";
-  args[31] = (char *) "abo22";
-  args[32] = (char *) "abo23";
-  args[33] = (char *) "abo24";
-  modify->add_compute(narg,args);
-  delete [] args;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixReaxFFSpecies::create_fix()
-{
-  int narg;
-  char **args;
-
-  narg = 37;
-  args = new char*[narg];
-  args[0]  = (char *) "SPECBOND";
-  args[1]  = (char *) "all";
-  args[2]  = (char *) "ave/atom";
-  args[3]  = tmparg[0];
-  args[4]  = tmparg[1];
-  args[5]  = tmparg[2];
-  args[6]  = (char *) "c_SPECATOM[1]";   // q, array_atoms[i][0]
-  args[7]  = (char *) "c_SPECATOM[2]";   // x, 1
-  args[8]  = (char *) "c_SPECATOM[3]";   // y, 2
-  args[9]  = (char *) "c_SPECATOM[4]";   // z, 3
-  args[10] = (char *) "c_SPECATOM[5]";   // vx, 4
-  args[11] = (char *) "c_SPECATOM[6]";   // vy, 5
-  args[12] = (char *) "c_SPECATOM[7]";   // vz, 6
-  args[13] = (char *) "c_SPECATOM[8]";   // abo01, 7
-  args[14] = (char *) "c_SPECATOM[9]";
-  args[15] = (char *) "c_SPECATOM[10]";
-  args[16] = (char *) "c_SPECATOM[11]";
-  args[17] = (char *) "c_SPECATOM[12]";
-  args[18] = (char *) "c_SPECATOM[13]";
-  args[19] = (char *) "c_SPECATOM[14]";
-  args[20] = (char *) "c_SPECATOM[15]";
-  args[21] = (char *) "c_SPECATOM[16]";
-  args[22] = (char *) "c_SPECATOM[17]";
-  args[23] = (char *) "c_SPECATOM[18]";
-  args[24] = (char *) "c_SPECATOM[19]"; // abo12, 18
-  args[25] = (char *) "c_SPECATOM[20]";
-  args[26] = (char *) "c_SPECATOM[21]";
-  args[27] = (char *) "c_SPECATOM[22]";
-  args[28] = (char *) "c_SPECATOM[23]";
-  args[29] = (char *) "c_SPECATOM[24]";
-  args[30] = (char *) "c_SPECATOM[25]";
-  args[31] = (char *) "c_SPECATOM[26]";
-  args[32] = (char *) "c_SPECATOM[27]";
-  args[33] = (char *) "c_SPECATOM[28]";
-  args[34] = (char *) "c_SPECATOM[29]";
-  args[35] = (char *) "c_SPECATOM[30]";
-  args[36] = (char *) "c_SPECATOM[31]";
-  modify->add_fix(narg,args);
-  f_SPECBOND = (FixAveAtom *) modify->fix[modify->nfix-1];
-  delete [] args;
-
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/REAXFF/fix_reaxff_species.h b/src/REAXFF/fix_reaxff_species.h
index 1957c1f8b1..b65ea49476 100644
--- a/src/REAXFF/fix_reaxff_species.h
+++ b/src/REAXFF/fix_reaxff_species.h
@@ -61,8 +61,6 @@ class FixReaxFFSpecies : public Fix {
   char *ele, **eletype, *filepos;
 
   void Output_ReaxFF_Bonds(bigint, FILE *);
-  void create_compute();
-  void create_fix();
   AtomCoord chAnchor(AtomCoord, AtomCoord);
   virtual void FindMolecule();
   void SortMolecule(int &);
@@ -82,7 +80,6 @@ class FixReaxFFSpecies : public Fix {
   class NeighList *list;
   class FixAveAtom *f_SPECBOND;
   class PairReaxFF *reaxff;
-
 };
 }
 
diff --git a/src/REAXFF/pair_reaxff.cpp b/src/REAXFF/pair_reaxff.cpp
index 0aea45b9cd..758ee70ab7 100644
--- a/src/REAXFF/pair_reaxff.cpp
+++ b/src/REAXFF/pair_reaxff.cpp
@@ -368,10 +368,8 @@ void PairReaxFF::init_style()
     error->warning(FLERR,"Total cutoff < 2*bond cutoff. May need to use an "
                    "increased neighbor list skin.");
 
-  if (fix_reaxff == nullptr) {
-    modify->add_fix(fmt::format("{} all REAXFF",fix_id));
-    fix_reaxff = (FixReaxFF *) modify->fix[modify->nfix-1];
-  }
+  if (fix_reaxff == nullptr)
+    fix_reaxff = (FixReaxFF *) modify->add_fix(fmt::format("{} all REAXFF",fix_id));
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/REPLICA/hyper.cpp b/src/REPLICA/hyper.cpp
index 84ff43ec85..d391604530 100644
--- a/src/REPLICA/hyper.cpp
+++ b/src/REPLICA/hyper.cpp
@@ -101,8 +101,7 @@ void Hyper::command(int narg, char **arg)
 
   // create FixEventHyper class to store event and pre-quench states
 
-  modify->add_fix("hyper_event all EVENT/HYPER");
-  fix_event = (FixEventHyper *) modify->fix[modify->nfix-1];
+  fix_event = (FixEventHyper *) modify->add_fix("hyper_event all EVENT/HYPER");
 
   // create Finish for timing output
 
diff --git a/src/REPLICA/prd.cpp b/src/REPLICA/prd.cpp
index 23e572e88b..501e7675a8 100644
--- a/src/REPLICA/prd.cpp
+++ b/src/REPLICA/prd.cpp
@@ -145,8 +145,7 @@ void PRD::command(int narg, char **arg)
 
   // create ComputeTemp class to monitor temperature
 
-  modify->add_compute("prd_temp all temp");
-  temperature = modify->compute[modify->ncompute-1];
+  temperature = modify->add_compute("prd_temp all temp");
 
   // create Velocity class for velocity creation in dephasing
   // pass it temperature compute, loop_setting, dist_setting settings
@@ -168,8 +167,7 @@ void PRD::command(int narg, char **arg)
 
   // create FixEventPRD class to store event and pre-quench states
 
-  modify->add_fix("prd_event all EVENT/PRD");
-  fix_event = (FixEventPRD *) modify->fix[modify->nfix-1];
+  fix_event = (FixEventPRD *) modify->add_fix("prd_event all EVENT/PRD");
 
   // create Finish for timing output
 
diff --git a/src/REPLICA/tad.cpp b/src/REPLICA/tad.cpp
index 5631944dbf..9b6ecf0b0c 100644
--- a/src/REPLICA/tad.cpp
+++ b/src/REPLICA/tad.cpp
@@ -131,28 +131,11 @@ void TAD::command(int narg, char **arg)
 
   // create FixEventTAD object to store last event
 
-  int narg2 = 3;
-  char **args = new char*[narg2];
-  args[0] = (char *) "tad_event";
-  args[1] = (char *) "all";
-  args[2] = (char *) "EVENT/TAD";
-  modify->add_fix(narg2,args);
-  fix_event = (FixEventTAD *) modify->fix[modify->nfix-1];
-  delete [] args;
+  fix_event = (FixEventTAD *) modify->add_fix("tad_event all EVENT/TAD");
 
   // create FixStore object to store revert state
 
-  narg2 = 6;
-  args = new char*[narg2];
-  args[0] = (char *) "tad_revert";
-  args[1] = (char *) "all";
-  args[2] = (char *) "STORE";
-  args[3] = (char *) "peratom";
-  args[4] = (char *) "0";
-  args[5] = (char *) "7";
-  modify->add_fix(narg2,args);
-  fix_revert = (FixStore *) modify->fix[modify->nfix-1];
-  delete [] args;
+  fix_revert = (FixStore *) modify->add_fix("tad_revert all STORE peratom 0 7");
 
   // create Finish for timing output
 
@@ -195,13 +178,10 @@ void TAD::command(int narg, char **arg)
 
   // set minimize style for quench
 
-  narg2 = 1;
-  args = new char*[narg2];
+  char *args[1];
   args[0] = min_style;
 
-  update->create_minimize(narg2,args,1);
-
-  delete [] args;
+  update->create_minimize(1,args,1);
 
   // init minimizer settings and minimizer itself
 
@@ -691,25 +671,13 @@ void TAD::perform_neb(int ievent)
 
   // create FixNEB object to support NEB
 
-  int narg2 = 4;
-  char **args = new char*[narg2];
-  args[0] = (char *) "neb";
-  args[1] = (char *) "all";
-  args[2] = (char *) "neb";
-  args[3] = (char *) "1.0";
-  modify->add_fix(narg2,args);
-  fix_neb = (Fix *) modify->fix[modify->nfix-1];
-  delete [] args;
+  fix_neb = (Fix *) modify->add_fix("neb all neb 1.0");
 
   // switch minimize style to quickmin for NEB
 
-  narg2 = 1;
-  args = new char*[narg2];
+  char *args[1];
   args[0] = min_style_neb;
-
-  update->create_minimize(narg2,args,1);
-
-  delete [] args;
+  update->create_minimize(1,args,1);
 
   // create NEB object
 
@@ -769,17 +737,12 @@ void TAD::perform_neb(int ievent)
 
   // switch minimize style back for quench
 
-  narg2 = 1;
-  args = new char*[narg2];
   args[0] = min_style;
-
-  update->create_minimize(narg2,args,1);
+  update->create_minimize(1,args,1);
 
   update->etol = etol;
   update->ftol = ftol;
 
-  delete [] args;
-
   // clean up
 
   modify->delete_fix("neb");
@@ -895,25 +858,14 @@ void TAD::delete_event_list() {
 
 void TAD::add_event()
 {
+  if (n_event_list == nmax_event_list)
+    grow_event_list(nmax_event_list+nmin_event_list);
 
   // create FixEventTAD object to store possible event
 
-  int narg = 3;
-  char **args = new char*[narg];
-
-  char str[128];
-  sprintf(str,"tad_event_%d",n_event_list);
-
-  args[0] = str;
-  args[1] = (char *) "all";
-  args[2] = (char *) "EVENT/TAD";
-  modify->add_fix(narg,args);
-
-  if (n_event_list == nmax_event_list)
-    grow_event_list(nmax_event_list+nmin_event_list);
-  n_event_list += 1;
-  int ievent = n_event_list-1;
-  fix_event_list[ievent] = (FixEventTAD *) modify->fix[modify->nfix-1];
+  int ievent = n_event_list++;
+  fix_event_list[ievent]
+    = (FixEventTAD *) modify->add_fix(fmt::format("tad_event_{} all EVENT/TAD", ievent));
 
   // store quenched state for new event
 
@@ -923,11 +875,6 @@ void TAD::add_event()
 
   fix_event->restore_state_quench();
   fix_event_list[ievent]->store_state_quench();
-
-  // string clean-up
-
-  delete [] args;
-
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/balance.cpp b/src/balance.cpp
index af142f404d..feb5eae708 100644
--- a/src/balance.cpp
+++ b/src/balance.cpp
@@ -499,8 +499,7 @@ void Balance::weight_storage(char *prefix)
   int ifix = modify->find_fix(cmd);
   if (ifix < 1) {
     cmd += " all STORE peratom 0 1";
-    modify->add_fix(cmd);
-    fixstore = (FixStore *) modify->fix[modify->nfix-1];
+    fixstore = (FixStore *) modify->add_fix(cmd);
   } else fixstore = (FixStore *) modify->fix[ifix];
 
   // do not carry weights with atoms during normal atom migration
diff --git a/src/compute_chunk_atom.cpp b/src/compute_chunk_atom.cpp
index 81eafcf8d6..5d789f722a 100644
--- a/src/compute_chunk_atom.cpp
+++ b/src/compute_chunk_atom.cpp
@@ -562,13 +562,9 @@ void ComputeChunkAtom::init()
   // fixstore initializes all values to 0.0
 
   if ((idsflag == ONCE || lockcount) && !fixstore) {
-    std::string cmd = id + std::string("_COMPUTE_STORE");
-    id_fix = new char[cmd.size()+1];
-    strcpy(id_fix,cmd.c_str());
-
-    cmd += fmt::format(" {} STORE peratom 1 1", group->names[igroup]);
-    modify->add_fix(cmd);
-    fixstore = (FixStore *) modify->fix[modify->nfix-1];
+    id_fix = utils::strdup(id + std::string("_COMPUTE_STORE"));
+    fixstore = (FixStore *) modify->add_fix(fmt::format("{} {} STORE peratom 1 1",
+                                                        id_fix, group->names[igroup]));
   }
 
   if ((idsflag != ONCE && !lockcount) && fixstore) {
diff --git a/src/compute_displace_atom.cpp b/src/compute_displace_atom.cpp
index 68776c1936..821525c53a 100644
--- a/src/compute_displace_atom.cpp
+++ b/src/compute_displace_atom.cpp
@@ -74,10 +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");
-  std::string cmd = id_fix + fmt::format(" {} STORE peratom 1 3",
-                                         group->names[igroup]);
-  modify->add_fix(cmd);
-  fix = (FixStore *) modify->fix[modify->nfix-1];
+  fix = (FixStore *) 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
diff --git a/src/compute_msd.cpp b/src/compute_msd.cpp
index ce9e10304e..209eec2810 100644
--- a/src/compute_msd.cpp
+++ b/src/compute_msd.cpp
@@ -14,16 +14,15 @@
 
 #include "compute_msd.h"
 
-#include <cstring>
-
 #include "atom.h"
-#include "update.h"
-#include "group.h"
 #include "domain.h"
-#include "modify.h"
-#include "fix_store.h"
 #include "error.h"
+#include "fix_store.h"
+#include "group.h"
+#include "modify.h"
+#include "update.h"
 
+#include <cstring>
 
 using namespace LAMMPS_NS;
 
@@ -66,13 +65,9 @@ ComputeMSD::ComputeMSD(LAMMPS *lmp, int narg, char **arg) :
   // create a new fix STORE style for reference positions
   // id = compute-ID + COMPUTE_STORE, fix group = compute group
 
-  std::string fixcmd = id + std::string("_COMPUTE_STORE");
-  id_fix = new char[fixcmd.size()+1];
-  strcpy(id_fix,fixcmd.c_str());
-
-  fixcmd += fmt::format(" {} STORE peratom 1 3",group->names[igroup]);
-  modify->add_fix(fixcmd);
-  fix = (FixStore *) modify->fix[modify->nfix-1];
+  id_fix = utils::strdup(id + std::string("_COMPUTE_STORE"));
+  fix = (FixStore *) 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
diff --git a/src/compute_msd_chunk.cpp b/src/compute_msd_chunk.cpp
index 74528da2f6..0d5b146338 100644
--- a/src/compute_msd_chunk.cpp
+++ b/src/compute_msd_chunk.cpp
@@ -14,18 +14,17 @@
 
 #include "compute_msd_chunk.h"
 
-#include <cstring>
-
 #include "atom.h"
-#include "group.h"
-#include "update.h"
-#include "modify.h"
 #include "compute_chunk_atom.h"
 #include "domain.h"
-#include "fix_store.h"
-#include "memory.h"
 #include "error.h"
+#include "fix_store.h"
+#include "group.h"
+#include "memory.h"
+#include "modify.h"
+#include "update.h"
 
+#include <cstring>
 
 using namespace LAMMPS_NS;
 
@@ -52,16 +51,14 @@ ComputeMSDChunk::ComputeMSDChunk(LAMMPS *lmp, int narg, char **arg) :
 
   // create a new fix STORE style for reference positions
   // id = compute-ID + COMPUTE_STORE, fix group = compute group
-  // do not know size of array at this point, just allocate 1x3 array
+  // do not know size of array at this point, just allocate 1x1 array
   // fix creation must be done now so that a restart run can
   //   potentially re-populate the fix array (and change it to correct size)
   // otherwise size reset and init will be done in setup()
 
   id_fix = utils::strdup(std::string(id) + "_COMPUTE_STORE");
-  std::string fixcmd = id_fix
-    + fmt::format(" {} STORE global 1 1",group->names[igroup]);
-  modify->add_fix(fixcmd);
-  fix = (FixStore *) modify->fix[modify->nfix-1];
+  fix = (FixStore *) modify->add_fix(fmt::format("{} {} STORE global 1 1",
+                                                 id_fix,group->names[igroup]));
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/compute_vacf.cpp b/src/compute_vacf.cpp
index ba88c6d2c1..90dcbf8d4f 100644
--- a/src/compute_vacf.cpp
+++ b/src/compute_vacf.cpp
@@ -14,8 +14,6 @@
 
 #include "compute_vacf.h"
 
-#include <cstring>
-
 #include "atom.h"
 #include "update.h"
 #include "group.h"
@@ -23,6 +21,7 @@
 #include "fix_store.h"
 #include "error.h"
 
+#include <cstring>
 
 using namespace LAMMPS_NS;
 
@@ -42,12 +41,9 @@ ComputeVACF::ComputeVACF(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 = new char[fixcmd.size()+1];
-  strcpy(id_fix,fixcmd.c_str());
-  fixcmd += fmt::format(" {} STORE peratom 1 3", group->names[igroup]);
-  modify->add_fix(fixcmd);
-  fix = (FixStore *) modify->fix[modify->nfix-1];
+  id_fix = utils::strdup(id + std::string("_COMPUTE_STORE"));
+  fix = (FixStore *) 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
diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp
index 6c417026ca..6bb7653e3b 100644
--- a/src/dump_custom.cpp
+++ b/src/dump_custom.cpp
@@ -1960,18 +1960,15 @@ int DumpCustom::modify_param(int narg, char **arg)
       thresh_last[nthresh] = -1;
     } else {
       thresh_fix = (FixStore **)
-        memory->srealloc(thresh_fix,(nthreshlast+1)*sizeof(FixStore *),
-                         "dump:thresh_fix");
+        memory->srealloc(thresh_fix,(nthreshlast+1)*sizeof(FixStore *),"dump:thresh_fix");
       thresh_fixID = (char **)
-        memory->srealloc(thresh_fixID,(nthreshlast+1)*sizeof(char *),
-                         "dump:thresh_fixID");
+        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);
-      modify->add_fix(fmt::format("{} {} STORE peratom 1 1",threshid,
-                                  group->names[igroup]));
-      thresh_fix[nthreshlast] = (FixStore *) modify->fix[modify->nfix-1];
+      threshid += fmt::format(" {} STORE peratom 1 1", group->names[igroup]);
+      thresh_fix[nthreshlast] = (FixStore *) modify->add_fix(threshid);
 
       thresh_last[nthreshlast] = nthreshlast;
       thresh_first[nthreshlast] = 1;
diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp
index 670c0dc4be..e5c2c3dece 100644
--- a/src/fix_adapt.cpp
+++ b/src/fix_adapt.cpp
@@ -252,12 +252,9 @@ void FixAdapt::post_constructor()
   id_fix_chg = nullptr;
 
   if (diamflag && atom->radius_flag) {
-    std::string fixcmd = id + std::string("_FIX_STORE_DIAM");
-    id_fix_diam = utils::strdup(fixcmd);
-    fixcmd += fmt::format(" {} STORE peratom 1 1",group->names[igroup]);
-    modify->add_fix(fixcmd);
-    fix_diam = (FixStore *) modify->fix[modify->nfix-1];
-
+    id_fix_diam = utils::strdup(id + std::string("_FIX_STORE_DIAM"));
+    fix_diam = (FixStore *) 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;
@@ -273,12 +270,9 @@ void FixAdapt::post_constructor()
   }
 
   if (chgflag && atom->q_flag) {
-    std::string fixcmd = id + std::string("_FIX_STORE_CHG");
-    id_fix_chg = utils::strdup(fixcmd);
-    fixcmd += fmt::format(" {} STORE peratom 1 1",group->names[igroup]);
-    modify->add_fix(fixcmd);
-    fix_chg = (FixStore *) modify->fix[modify->nfix-1];
-
+    id_fix_chg = utils::strdup(id + std::string("_FIX_STORE_CHG"));
+    fix_chg = (FixStore *) 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;
diff --git a/src/min.cpp b/src/min.cpp
index 8fe73dc9e1..69e24a541f 100644
--- a/src/min.cpp
+++ b/src/min.cpp
@@ -121,8 +121,7 @@ void Min::init()
   // create fix needed for storing atom-based quantities
   // will delete it at end of run
 
-  modify->add_fix("MINIMIZE all MINIMIZE");
-  fix_minimize = (FixMinimize *) modify->fix[modify->nfix-1];
+  fix_minimize = (FixMinimize *) modify->add_fix("MINIMIZE all MINIMIZE");
 
   // clear out extra global and per-atom dof
   // will receive requests for new per-atom dof during pair init()
diff --git a/src/modify.cpp b/src/modify.cpp
index b2670fe5f9..995b3b82ac 100644
--- a/src/modify.cpp
+++ b/src/modify.cpp
@@ -37,7 +37,6 @@ using namespace FixConst;
 
 #define DELTA 4
 #define BIG 1.0e20
-#define NEXCEPT 7       // change when add to exceptions in add_fix()
 
 /* ---------------------------------------------------------------------- */
 
@@ -803,7 +802,7 @@ int Modify::min_reset_ref()
    add a new fix or replace one with same ID
 ------------------------------------------------------------------------- */
 
-void Modify::add_fix(int narg, char **arg, int trysuffix)
+Fix *Modify::add_fix(int narg, char **arg, int trysuffix)
 {
   if (narg < 3) error->all(FLERR,"Illegal fix command");
 
@@ -956,13 +955,14 @@ void Modify::add_fix(int narg, char **arg, int trysuffix)
   if (newflag) nfix++;
   fmask[ifix] = fix[ifix]->setmask();
   fix[ifix]->post_constructor();
+  return fix[ifix];
 }
 
 /* ----------------------------------------------------------------------
    convenience function to allow adding a fix from a single string
 ------------------------------------------------------------------------- */
 
-void Modify::add_fix(const std::string &fixcmd, int trysuffix)
+Fix *Modify::add_fix(const std::string &fixcmd, int trysuffix)
 {
   auto args = utils::split_words(fixcmd);
   std::vector<char *> newarg(args.size());
@@ -970,7 +970,7 @@ void Modify::add_fix(const std::string &fixcmd, int trysuffix)
   for (const auto &arg : args) {
     newarg[i++] = (char *)arg.c_str();
   }
-  add_fix(args.size(),newarg.data(),trysuffix);
+  return add_fix(args.size(),newarg.data(),trysuffix);
 }
 
 
@@ -981,8 +981,7 @@ void Modify::add_fix(const std::string &fixcmd, int trysuffix)
         replace it later with the desired Fix instance
 ------------------------------------------------------------------------- */
 
-void Modify::replace_fix(const char *replaceID,
-                         int narg, char **arg, int trysuffix)
+Fix *Modify::replace_fix(const char *replaceID, int narg, char **arg, int trysuffix)
 {
   int ifix = find_fix(replaceID);
   if (ifix < 0) error->all(FLERR,"Modify replace_fix ID {} could not be found", replaceID);
@@ -1007,24 +1006,22 @@ void Modify::replace_fix(const char *replaceID,
   // invoke add_fix
   // it will find and overwrite the replaceID fix
 
-  add_fix(narg,arg,trysuffix);
+  return add_fix(narg,arg,trysuffix);
 }
 
 /* ----------------------------------------------------------------------
    convenience function to allow replacing a fix from a single string
 ------------------------------------------------------------------------- */
 
-void Modify::replace_fix(const std::string &oldfix,
-                         const std::string &fixcmd, int trysuffix)
+Fix *Modify::replace_fix(const std::string &oldfix, const std::string &fixcmd, int trysuffix)
 {
   auto args = utils::split_words(fixcmd);
-  char **newarg = new char*[args.size()];
-  int i=0;
+  std::vector<char *> newarg(args.size());
+  int i = 0;
   for (const auto &arg : args) {
     newarg[i++] = (char *)arg.c_str();
   }
-  replace_fix(oldfix.c_str(),args.size(),newarg,trysuffix);
-  delete[] newarg;
+  return replace_fix(oldfix.c_str(),args.size(),newarg.data(),trysuffix);
 }
 
 /* ----------------------------------------------------------------------
@@ -1214,7 +1211,7 @@ int Modify::check_rigid_list_overlap(int *select)
    add a new compute
 ------------------------------------------------------------------------- */
 
-void Modify::add_compute(int narg, char **arg, int trysuffix)
+Compute *Modify::add_compute(int narg, char **arg, int trysuffix)
 {
   if (narg < 3) error->all(FLERR,"Illegal compute command");
 
@@ -1269,23 +1266,22 @@ void Modify::add_compute(int narg, char **arg, int trysuffix)
   if (compute[ncompute] == nullptr)
     error->all(FLERR,utils::check_packages_for_style("compute",arg[2],lmp));
 
-  ncompute++;
+  return compute[ncompute++];
 }
 
 /* ----------------------------------------------------------------------
    convenience function to allow adding a compute from a single string
 ------------------------------------------------------------------------- */
 
-void Modify::add_compute(const std::string &computecmd, int trysuffix)
+Compute *Modify::add_compute(const std::string &computecmd, int trysuffix)
 {
   auto args = utils::split_words(computecmd);
-  char **newarg = new char*[args.size()];
+  std::vector<char *>newarg(args.size());
   int i=0;
   for (const auto &arg : args) {
     newarg[i++] = (char *)arg.c_str();
   }
-  add_compute(args.size(),newarg,trysuffix);
-  delete[] newarg;
+  return add_compute(args.size(),newarg.data(),trysuffix);
 }
 
 
diff --git a/src/modify.h b/src/modify.h
index 5111ef43ab..9446f285e7 100644
--- a/src/modify.h
+++ b/src/modify.h
@@ -100,18 +100,18 @@ class Modify : protected Pointers {
   virtual int min_dof();
   virtual int min_reset_ref();
 
-  void add_fix(int, char **, int trysuffix = 1);
-  void add_fix(const std::string &, int trysuffix = 1);
-  void replace_fix(const char *, int, char **, int trysuffix = 1);
-  void replace_fix(const std::string &, const std::string &, int trysuffix = 1);
+  Fix *add_fix(int, char **, int trysuffix = 1);
+  Fix *add_fix(const std::string &, int trysuffix = 1);
+  Fix *replace_fix(const char *, int, char **, int trysuffix = 1);
+  Fix *replace_fix(const std::string &, const std::string &, int trysuffix = 1);
   void modify_fix(int, char **);
   void delete_fix(const std::string &);
   void delete_fix(int);
   int find_fix(const std::string &);
   int find_fix_by_style(const char *);
 
-  void add_compute(int, char **, int trysuffix = 1);
-  void add_compute(const std::string &, int trysuffix = 1);
+  Compute *add_compute(int, char **, int trysuffix = 1);
+  Compute *add_compute(const std::string &, int trysuffix = 1);
   void modify_compute(int, char **);
   void delete_compute(const std::string &);
   void delete_compute(int);
diff --git a/src/respa.cpp b/src/respa.cpp
index c51ec04131..5190441994 100644
--- a/src/respa.cpp
+++ b/src/respa.cpp
@@ -302,9 +302,8 @@ void Respa::init()
   // if supported, we also store torques on a per-level basis
 
   std::string cmd = fmt::format("RESPA all RESPA {}",nlevels);
-  if (atom->torque_flag) modify->add_fix(cmd + " torque");
-  else modify->add_fix(cmd);
-  fix_respa = (FixRespa *) modify->fix[modify->nfix-1];
+  if (atom->torque_flag) cmd += " torque";
+  fix_respa = (FixRespa *) modify->add_fix(cmd);
 
   // insure respa inner/middle/outer is using Pair class that supports it
 
diff --git a/src/variable.cpp b/src/variable.cpp
index 2ef1957a09..a7b2955b33 100644
--- a/src/variable.cpp
+++ b/src/variable.cpp
@@ -5065,16 +5065,10 @@ VarReader::VarReader(LAMMPS *lmp, char *name, char *file, int flag) :
 
   if (style == Variable::ATOMFILE) {
     if (atom->map_style == Atom::MAP_NONE)
-      error->all(FLERR,"Cannot use atomfile-style "
-                 "variable unless an atom map exists");
-
-    std::string cmd = name + std::string("_VARIABLE_STORE");
-    id_fix = utils::strdup(cmd);
-
-    cmd += " all STORE peratom 0 1";
-    modify->add_fix(cmd);
-    fixstore = (FixStore *) modify->fix[modify->nfix-1];
+      error->all(FLERR,"Cannot use atomfile-style variable unless an atom map exists");
 
+    id_fix = utils::strdup(std::string(name) + "_VARIABLE_STORE");
+    fixstore = (FixStore *) modify->add_fix(std::string(id_fix) + " all STORE peratom 0 1");
     buffer = new char[CHUNK*MAXLINE];
   }
 }
diff --git a/src/velocity.cpp b/src/velocity.cpp
index acc374fa4a..fa09f451df 100644
--- a/src/velocity.cpp
+++ b/src/velocity.cpp
@@ -187,11 +187,11 @@ void Velocity::create(double t_desired, int seed)
   Compute *temperature_nobias = nullptr;
 
   if (temperature == nullptr || bias_flag) {
-    modify->add_compute(fmt::format("velocity_temp {} temp",group->names[igroup]));
+    auto newcompute = modify->add_compute(fmt::format("velocity_temp {} temp",group->names[igroup]));
     if (temperature == nullptr) {
-      temperature = modify->compute[modify->ncompute-1];
+      temperature = newcompute;
       tcreate_flag = 1;
-    } else temperature_nobias = modify->compute[modify->ncompute-1];
+    } else temperature_nobias = newcompute;
   }
 
   // initialize temperature computation(s)
@@ -575,8 +575,7 @@ void Velocity::scale(int /*narg*/, char **arg)
 
   int tflag = 0;
   if (temperature == nullptr) {
-    modify->add_compute(fmt::format("velocity_temp {} temp",group->names[igroup]));
-    temperature = modify->compute[modify->ncompute-1];
+    temperature = modify->add_compute(fmt::format("velocity_temp {} temp",group->names[igroup]));
     tflag = 1;
   }
 
@@ -849,15 +848,11 @@ void Velocity::options(int narg, char **arg)
       iarg += 2;
     } else if (strcmp(arg[iarg],"temp") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
-      int icompute;
-      for (icompute = 0; icompute < modify->ncompute; icompute++)
-        if (strcmp(arg[iarg+1],modify->compute[icompute]->id) == 0) break;
-      if (icompute == modify->ncompute)
-        error->all(FLERR,"Could not find velocity temperature ID");
+      int icompute = modify->find_compute(arg[iarg+1]);
+      if (icompute < 0) error->all(FLERR,"Could not find velocity temperature ID");
       temperature = modify->compute[icompute];
       if (temperature->tempflag == 0)
-        error->all(FLERR,
-                   "Velocity temperature ID does not compute temperature");
+        error->all(FLERR,"Velocity temperature ID does not compute temperature");
       iarg += 2;
     } else if (strcmp(arg[iarg],"bias") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");