Merge remote-tracking branch 'github/develop' into distributed-grids

doc/src/Developer_updating.rst

@@ -24,6 +24,7 @@ Available topics in mostly chronological order are:
 - `Use of "override" instead of "virtual"`_
 - `Simplified and more compact neighbor list requests`_
 - `Split of fix STORE into fix STORE/GLOBAL and fix STORE/PERATOM`_
+- `Use Output::get_dump_by_id() instead of Output::find_dump()`_
 
 ----
 
@@ -382,3 +383,43 @@ New:
    FixStoreGlobal *fix = dynamic_cast<FixStoreGlobal *>(modify->get_fix_by_id(id_fix));
 
 This change is **required** or else the code will not compile.
+
+Use Output::get_dump_by_id() instead of Output::find_dump()
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. versionchanged:: TBD
+
+The accessor function to individual dump style instances has been changed
+from ``Output::find_dump()`` returning the index of the dump instance in
+the list of dumps to ``Output::get_dump_by_id()`` returning a pointer to
+the dump directly.  Example:
+
+Old:
+
+.. code-block:: C++
+
+   int idump = output->find_dump(arg[iarg+1]);
+   if (idump < 0)
+     error->all(FLERR,"Dump ID in hyper command does not exist");
+   memory->grow(dumplist,ndump+1,"hyper:dumplist");
+   dumplist[ndump++] = idump;
+
+   [...]
+
+   if (dumpflag)
+     for (int idump = 0; idump < ndump; idump++)
+       output->dump[dumplist[idump]]->write();
+
+New:
+
+.. code-block:: C++
+
+   auto idump = output->get_dump_by_id(arg[iarg+1]);
+   if (!idump) error->all(FLERR,"Dump ID {} in hyper command does not exist", arg[iarg+1]);
+   dumplist.emplace_back(idump);
+
+   [...]
+
+   if (dumpflag) for (auto idump : dumplist) idump->write();
+
+This change is **required** or else the code will not compile.

doc/src/fix_shake.rst

@@ -33,12 +33,14 @@ Syntax
        *m* value = one or more mass values
 
 * zero or more keyword/value pairs may be appended
-* keyword = *mol*
+* keyword = *mol* or *kbond*
 
   .. parsed-literal::
 
        *mol* value = template-ID
          template-ID = ID of molecule template specified in a separate :doc:`molecule <molecule>` command
+       *kbond* value = force constant
+         force constant = force constant used to apply a restraint force when used during minimization
 
 Examples
 """"""""
@@ -56,7 +58,10 @@ Description
 
 Apply bond and angle constraints to specified bonds and angles in the
 simulation by either the SHAKE or RATTLE algorithms.  This typically
-enables a longer timestep.
+enables a longer timestep.  The SHAKE or RATTLE algorithms, however, can
+*only* be applied during molecular dynamics runs.  When this fix is used
+during a minimization, the constraints are *approximated* by strong
+harmonic restraints.
 
 **SHAKE vs RATTLE:**
 
@@ -90,7 +95,7 @@ The SHAKE algorithm satisfies the first condition, i.e. the sites at
 time *n+1* will have the desired separations Dij immediately after the
 coordinates are integrated.  If we also enforce the second condition,
 the velocity components along the bonds will vanish.  RATTLE satisfies
-both conditions.  As implemented in LAMMPS, fix rattle uses fix shake
+both conditions.  As implemented in LAMMPS, *fix rattle* uses fix shake
 for satisfying the coordinate constraints. Therefore the settings and
 optional keywords are the same for both fixes, and all the information
 below about SHAKE is also relevant for RATTLE.
@@ -139,30 +144,47 @@ for.
 
 .. note::
 
-   This command works by using the current forces on atoms to
+   This command works by using the current forces on atoms to calculate
-   calculate an additional constraint force which when added will leave
+   an additional constraint force which when added will leave the atoms
-   the atoms in positions that satisfy the SHAKE constraints (e.g. bond
+   in positions that satisfy the SHAKE constraints (e.g. bond length)
-   length) after the next time integration step.  If you define fixes
+   after the next time integration step.  If you define fixes
-   (e.g. :doc:`fix efield <fix_efield>`) that add additional force to the
+   (e.g. :doc:`fix efield <fix_efield>`) that add additional force to
-   atoms after fix shake operates, then this fix will not take them into
+   the atoms after *fix shake* operates, then this fix will not take them
-   account and the time integration will typically not satisfy the SHAKE
+   into account and the time integration will typically not satisfy the
-   constraints.  The solution for this is to make sure that fix shake is
+   SHAKE constraints.  The solution for this is to make sure that fix
-   defined in your input script after any other fixes which add or change
+   shake is defined in your input script after any other fixes which add
-   forces (to atoms that fix shake operates on).
+   or change forces (to atoms that *fix shake* operates on).
 
 ----------
 
-The *mol* keyword should be used when other commands, such as :doc:`fix deposit <fix_deposit>` or :doc:`fix pour <fix_pour>`, add molecules
+The *mol* keyword should be used when other commands, such as :doc:`fix
+deposit <fix_deposit>` or :doc:`fix pour <fix_pour>`, add molecules
 on-the-fly during a simulation, and you wish to constrain the new
 molecules via SHAKE.  You specify a *template-ID* previously defined
 using the :doc:`molecule <molecule>` command, which reads a file that
 defines the molecule.  You must use the same *template-ID* that the
 command adding molecules uses.  The coordinates, atom types, special
-bond restrictions, and SHAKE info can be specified in the molecule
+bond restrictions, and SHAKE info can be specified in the molecule file.
-file.  See the :doc:`molecule <molecule>` command for details.  The only
+See the :doc:`molecule <molecule>` command for details.  The only
 settings required to be in this file (by this command) are the SHAKE
 info of atoms in the molecule.
 
+The *kbond* keyword sets the restraint force constant when *fix shake*
+or fix rattle are used during minimization.  In that case the constraint
+algorithms are *not* applied and restraint forces are used instead to
+maintain the geometries similar to the constraints.  How well the
+geometries are maintained and how quickly a minimization converges,
+depends largely on the force constant *kbond*: larger values will reduce
+the deviation from the desired geometry, but can also lead to slower
+convergence of the minimization or lead to instabilities depending on
+the minimization algorithm requiring to reduce the value of
+:doc:`timestep <timestep>`.  Even though the restraints will not
+preserve the bond lengths and angles as closely as the constraints
+during the MD, they are generally close enough so that the constraints
+will be fulfilled to the desired accuracy within a few MD steps
+following the minimization. The default value for *kbond* depends on the
+:doc:`units <units>` setting and is 1.0e6*k_B.
+
 ----------
 
 .. include:: accel_styles.rst
@@ -177,10 +199,10 @@ LAMMPS closely follows (:ref:`Andersen (1983) <Andersen3>`).
 
 .. note::
 
-   The fix rattle command modifies forces and velocities and thus
+   The *fix rattle* command modifies forces and velocities and thus
    should be defined after all other integration fixes in your input
    script.  If you define other fixes that modify velocities or forces
-   after fix rattle operates, then fix rattle will not take them into
+   after *fix rattle* operates, then *fix rattle* will not take them into
    account and the overall time integration will typically not satisfy
    the RATTLE constraints.  You can check whether the constraints work
    correctly by setting the value of RATTLE_DEBUG in src/fix_rattle.cpp
@@ -194,19 +216,32 @@ Restart, fix_modify, output, run start/stop, minimize info
 No information about these fixes is written to :doc:`binary restart
 files <restart>`.
 
-The :doc:`fix_modify <fix_modify>` *virial* option is supported by
+Both fix *shake* and fix *rattle* behave differently during a minimization
-these fixes to add the contribution due to the added forces on atoms
+in comparison to a molecular dynamics run:
-to both the global pressure and per-atom stress of the system via the
+
-:doc:`compute pressure <compute_pressure>` and :doc:`compute
+- When used during a minimization, the SHAKE or RATTLE constraint
-stress/atom <compute_stress_atom>` commands.  The former can be
+  algorithms themselves are **not** applied.  Instead the constraints
-accessed by :doc:`thermodynamic output <thermo_style>`.  The default
+  are replaced by harmonic restraint forces.  The energy and virial
-setting for this fix is :doc:`fix_modify virial yes <fix_modify>`.
+  contributions due to the restraint forces are tallied into global and
+  per-atom accumulators. The total restraint energy is also accessible
+  as a global scalar property of the fix.
+
+- During molecular dynamics runs, however, the fixes do apply the
+  requested SHAKE or RATTLE constraint algorithms.
+
+  The :doc:`fix_modify <fix_modify>` *virial* option is supported by
+  these fixes to add the contribution due to the added constraint forces
+  on atoms to both the global pressure and per-atom stress of the system
+  via the :doc:`compute pressure <compute_pressure>` and :doc:`compute
+  stress/atom <compute_stress_atom>` commands.  The former can be
+  accessed by :doc:`thermodynamic output <thermo_style>`.
+
+  The  default  setting  for  this  fix  is  :doc:`fix_modify  virial  yes
+  <fix_modify>`.  No global  or per-atom  quantities are  stored by  these
+  fixes for access by various :doc:`output commands <Howto_output>` during
+  an  MD  run.   No  parameter  of  these  fixes  can  be  used  with  the
+  *start/stop* keywords of the :doc:`run <run>` command.
 
-No global or per-atom quantities are stored by these fixes for access
-by various :doc:`output commands <Howto_output>`.  No parameter of
-these fixes can be used with the *start/stop* keywords of the
-:doc:`run <run>` command.  These fixes are not invoked during
-:doc:`energy minimization <minimize>`.
 
 Restrictions
 """"""""""""
@@ -224,21 +259,27 @@ which can lead to poor energy conservation.  You can test for this in
 your system by running a constant NVE simulation with a particular set
 of SHAKE parameters and monitoring the energy versus time.
 
-SHAKE or RATTLE should not be used to constrain an angle at 180
+SHAKE or RATTLE should not be used to constrain an angle at 180 degrees
-degrees (e.g. linear CO2 molecule).  This causes numeric difficulties.
+(e.g. a linear CO2 molecule).  This causes a divergence when solving the
-You can use :doc:`fix rigid or fix rigid/small <fix_rigid>` instead to
+constraint equations numerically.  You can use :doc:`fix rigid or fix
-make a linear molecule rigid.
+rigid/small <fix_rigid>` instead to make a linear molecule rigid.
+
+When used during minimization choosing a too large value of the *kbond*
+can make minimization very inefficient and also cause stability problems
+with some minimization algorithms.  Sometimes those can be avoided by
+reducing the :doc:`timestep <timestep>`.
 
 Related commands
 """"""""""""""""
 
-none
+`fix rigid <fix_rigid>`, `fix ehex <fix_ehex>`,
+`fix nve/manifold/rattle <fix_nve_manifold_rattle>`
 
 
 Default
 """""""
 
-none
+kbond = 1.0e9*k_B
 
 ----------
 

doc/src/minimize.rst

@@ -32,30 +32,51 @@ Description
 Perform an energy minimization of the system, by iteratively adjusting
 atom coordinates.  Iterations are terminated when one of the stopping
 criteria is satisfied.  At that point the configuration will hopefully
-be in local potential energy minimum.  More precisely, the
+be in a local potential energy minimum.  More precisely, the
 configuration should approximate a critical point for the objective
 function (see below), which may or may not be a local minimum.
 
-The minimization algorithm used is set by the
+The minimization algorithm used is set by the :doc:`min_style
-:doc:`min_style <min_style>` command.  Other options are set by the
+<min_style>` command.  Other options are set by the :doc:`min_modify
-:doc:`min_modify <min_modify>` command.  Minimize commands can be
+<min_modify>` command.  Minimize commands can be interspersed with
-interspersed with :doc:`run <run>` commands to alternate between
+:doc:`run <run>` commands to alternate between relaxation and dynamics.
-relaxation and dynamics.  The minimizers bound the distance atoms move
+The minimizers bound the distance atoms may move in one iteration, so
-in one iteration, so that you can relax systems with highly overlapped
+that you can relax systems with highly overlapped atoms (large energies
-atoms (large energies and forces) by pushing the atoms off of each
+and forces) by pushing the atoms off of each other.
-other.
 
-Alternate means of relaxing a system are to run dynamics with a small
+.. admonition:: Neighbor list update settings
-or :doc:`limited timestep <fix_nve_limit>`.  Or dynamics can be run
+   :class: note
-using :doc:`fix viscous <fix_viscous>` to impose a damping force that
+
-slowly drains all kinetic energy from the system.  The :doc:`pair_style soft <pair_soft>` potential can be used to un-overlap atoms while
+   The distance that atoms can move during individual minimization steps
-running dynamics.
+   can be quite large, especially at the beginning of a minimization.
-un-overlap atoms while running dynamics.
+   Thus `neighbor list settings <neigh_modify>` of *every = 1* and
+   *delay = 0* are **required**.  This may be combined with either
+   *check = no* (always update the neighbor list) or *check = yes* (only
+   update the neighbor list if at least one atom has moved more than
+   half the `neighbor list skin <neighbor>` distance since the last
+   reneighboring).  Using *check = yes* is recommended since it avoids
+   unneeded reneighboring steps when the system is closer to the minimum
+   and thus atoms move only small distances.  Using *check = no* may
+   be required for debugging or when coupling LAMMPS with external
+   codes that require a predictable sequence of neighbor list updates.
+
+   If the settings are **not** *every = 1* and *delay = 0*, LAMMPS
+   will temporarily apply a `neigh_modify every 1 delay 0 check yes
+   <neigh_modify>` setting during the minimization and restore the
+   original setting at the end of the minimization.  A corresponding
+   message will be printed to the screen and log file, if this happens.
+
+Alternate means of relaxing a system are to run dynamics with a small or
+:doc:`limited timestep <fix_nve_limit>`.  Or dynamics can be run using
+:doc:`fix viscous <fix_viscous>` to impose a damping force that slowly
+drains all kinetic energy from the system.  The :doc:`pair_style soft
+<pair_soft>` potential can be used to un-overlap atoms while running
+dynamics.
 
 Note that you can minimize some atoms in the system while holding the
-coordinates of other atoms fixed by applying :doc:`fix setforce
+coordinates of other atoms fixed by applying :doc:`fix setforce 0.0 0.0
-<fix_setforce>` to the other atoms.  See a fuller discussion of using
+0.0 <fix_setforce>` to the other atoms.  See a more detailed discussion
-fixes while minimizing below.
+of :ref:`using fixes while minimizing below <fix_minimize>`.
 
 The :doc:`minimization styles <min_style>` *cg*, *sd*, and *hftn*
 involves an outer iteration loop which sets the search direction along
@@ -74,10 +95,11 @@ they require a :doc:`timestep <timestep>` be defined.
 
 .. note::
 
-   The damped dynamic minimizers use whatever timestep you have
+   The damped dynamic minimizer algorithms will use the timestep you
-   defined via the :doc:`timestep <timestep>` command.  Often they
+   have defined via the :doc:`timestep <timestep>` command or its
-   will converge more quickly if you use a timestep about 10x larger
+   default value.  Often they will converge more quickly if you use a
-   than you would normally use for dynamics simulations.
+   timestep about 10x larger than you would normally use for regular
+   molecular dynamics simulations.
 
 ----------
 
@@ -210,34 +232,42 @@ The iterations and force evaluation values are what is checked by the
 
 .. note::
 
-   There are several force fields in LAMMPS which have
+   There are several force fields in LAMMPS which have discontinuities
-   discontinuities or other approximations which may prevent you from
+   or other approximations which may prevent you from performing an
-   performing an energy minimization to high tolerances.  For example,
+   energy minimization to tight tolerances.  For example, you should use
-   you should use a :doc:`pair style <pair_style>` that goes to 0.0 at the
+   a :doc:`pair style <pair_style>` that goes to 0.0 at the cutoff
-   cutoff distance when performing minimization (even if you later change
+   distance when performing minimization (even if you later change it
-   it when running dynamics).  If you do not do this, the total energy of
+   when running dynamics).  If you do not do this, the total energy of
    the system will have discontinuities when the relative distance
-   between any pair of atoms changes from cutoff+epsilon to
+   between any pair of atoms changes from cutoff *plus* epsilon to
-   cutoff-epsilon and the minimizer may behave poorly.  Some of the
+   cutoff *minus* epsilon and the minimizer may thus behave poorly.
-   many-body potentials use splines and other internal cutoffs that
+   Some of the many-body potentials use splines and other internal
-   inherently have this problem.  The :doc:`long-range Coulombic styles <kspace_style>` (PPPM, Ewald) are approximate to within the
+   cutoffs that inherently have this problem.  The :doc:`long-range
-   user-specified tolerance, which means their energy and forces may not
+   Coulombic styles <kspace_style>` (PPPM, Ewald) are approximate to
-   agree to a higher precision than the Kspace-specified tolerance.  In
+   within the user-specified tolerance, which means their energy and
-   all these cases, the minimizer may give up and stop before finding a
+   forces may not agree to a higher precision than the Kspace-specified
-   minimum to the specified energy or force tolerance.
+   tolerance.  This agreement is further reduced when using tabulation
+   to speed up the computation of the real-space part of the Coulomb
+   interactions, which is enabled by default.  In all these cases, the
+   minimizer may give up and stop before finding a minimum to the
+   specified energy or force tolerance.
 
 Note that a cutoff Lennard-Jones potential (and others) can be shifted
-so that its energy is 0.0 at the cutoff via the
+so that its energy is 0.0 at the cutoff via the :doc:`pair_modify
-:doc:`pair_modify <pair_modify>` command.  See the doc pages for
+<pair_modify>` command.  See the doc pages for individual :doc:`pair
-individual :doc:`pair styles <pair_style>` for details.  Note that
+styles <pair_style>` for details.  Note that most Coulombic potentials
-Coulombic potentials always have a cutoff, unless versions with a
+have a cutoff, unless versions with a long-range component are used
-long-range component are used (e.g. :doc:`pair_style lj/cut/coul/long <pair_lj_cut_coul>`).  The CHARMM potentials go to 0.0 at
+(e.g. :doc:`pair_style lj/cut/coul/long <pair_lj_cut_coul>`) or some
-the cutoff (e.g. :doc:`pair_style lj/charmm/coul/charmm <pair_charmm>`),
+other damping/smoothing schemes are used.  The CHARMM potentials go to
-as do the GROMACS potentials (e.g. :doc:`pair_style lj/gromacs <pair_gromacs>`).
+0.0 at the cutoff (e.g. :doc:`pair_style lj/charmm/coul/charmm
+<pair_charmm>`), as do the GROMACS potentials (e.g. :doc:`pair_style
+lj/gromacs <pair_gromacs>`).
 
 If a soft potential (:doc:`pair_style soft <pair_soft>`) is used the
 Astop value is used for the prefactor (no time dependence).
 
+.. _fix_minimize:
+
 The :doc:`fix box/relax <fix_box_relax>` command can be used to apply an
 external pressure to the simulation box and allow it to shrink/expand
 during the minimization.

doc/utils/sphinx-config/false_positives.txt

@@ -13,6 +13,7 @@ Accelrys
 acceptor
 Acceptor
 acceptors
+accessor
 accomx
 accomy
 accomz
@@ -1623,6 +1624,7 @@ kb
 kB
 kbit
 kbits
+kbond
 kcal
 kcl
 Kd

src/AMOEBA/amoeba_induce.cpp

@@ -53,8 +53,7 @@ enum{GORDON1,GORDON2};
 void PairAmoeba::induce()
 {
   bool done;
-  int i,j,m,itype;
+  int i,j,m,itype,iter;
-  int iter,maxiter;
   double polmin;
   double eps,epsold;
   double epsd,epsp;
@@ -159,7 +158,6 @@ void PairAmoeba::induce()
 
   if (poltyp == MUTUAL) {
     done = false;
-    maxiter = 100;
     iter = 0;
     polmin = 0.00000001;
     eps = 100.0;

src/DIFFRACTION/compute_saed.cpp

@@ -37,7 +37,7 @@ using namespace LAMMPS_NS;
 using namespace MathConst;
 
 static const char cite_compute_saed_c[] =
-  "compute_saed command: doi:10.1088/0965-0393/21/5/055020\n\n"
+  "compute saed command: doi:10.1088/0965-0393/21/5/055020\n\n"
   "@Article{Coleman13,\n"
   " author = {S. P. Coleman and D. E. Spearot and L. Capolungo},\n"
   " title = {Virtual Diffraction Analysis of {Ni} [010] Symmetric Tilt Grain Boundaries},\n"

src/DIFFRACTION/compute_xrd.cpp

@@ -38,7 +38,7 @@ using namespace LAMMPS_NS;
 using namespace MathConst;
 
 static const char cite_compute_xrd_c[] =
-  "compute_xrd command: doi:10.1088/0965-0393/21/5/055020\n\n"
+  "compute xrd command: doi:10.1088/0965-0393/21/5/055020\n\n"
   "@Article{Coleman13,\n"
   " author = {S. P. Coleman and D. E. Spearot and L. Capolungo},\n"
   " title = {Virtual Diffraction Analysis of {Ni} [010] Symmetric Tilt Grain Boundaries},\n"

src/H5MD/dump_h5md.cpp

@@ -73,7 +73,7 @@ DumpH5MD::DumpH5MD(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg)
   format_default = nullptr;
   flush_flag = 0;
   unwrap_flag = 0;
-  datafile_from_dump = -1;
+  other_dump = nullptr;
   author_name = nullptr;
 
   every_dump = utils::inumeric(FLERR,arg[3],false,lmp);
@@ -141,11 +141,8 @@ DumpH5MD::DumpH5MD(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg)
       }
       if (box_is_set||create_group_is_set)
         error->all(FLERR, "Cannot set file_from in dump h5md after box or create_group");
-      int idump;
+      other_dump = dynamic_cast<DumpH5MD *>(output->get_dump_by_id(arg[iarg+1]));
-      for (idump = 0; idump < output->ndump; idump++)
+      if (!other_dump) error->all(FLERR,"Cound not find dump_modify H5MD dump ID {}", arg[iarg+1]);
-        if (strcmp(arg[iarg+1],output->dump[idump]->id) == 0) break;
-      if (idump == output->ndump) error->all(FLERR,"Cound not find dump_modify ID");
-      datafile_from_dump = idump;
       do_box=false;
       create_group=false;
       iarg+=2;
@@ -262,7 +259,7 @@ void DumpH5MD::openfile()
   }
 
   if (me == 0) {
-    if (datafile_from_dump<0) {
+    if (!other_dump) {
       if (author_name==nullptr) {
         datafile = h5md_create_file(filename, "N/A", nullptr, "lammps", LAMMPS_VERSION);
       } else {
@@ -296,8 +293,6 @@ void DumpH5MD::openfile()
         h5md_write_string_attribute(particles_data.group, "charge", "type", "effective");
       }
     } else {
-      DumpH5MD* other_dump;
-      other_dump=(DumpH5MD*)output->dump[datafile_from_dump];
       datafile = other_dump->datafile;
       group_name_length = strlen(group->names[igroup]);
       group_name = new char[group_name_length];

src/H5MD/dump_h5md.h

@@ -35,9 +35,9 @@ class DumpH5MD : public Dump {
   int natoms, ntotal;
   int unwrap_flag;    // 1 if atom coords are unwrapped, 0 if no
   h5md_file datafile;
-  int datafile_from_dump;
   h5md_particles_group particles_data;
   char *author_name;
+  DumpH5MD *other_dump;
 
   bool do_box;
   bool create_group;

src/KIM/kim_command.cpp

@@ -70,7 +70,7 @@
 using namespace LAMMPS_NS;
 
 static constexpr const char *const cite_openkim =
-  "OpenKIM: https://doi.org/10.1007/s11837-011-0102-6\n\n"
+  "OpenKIM Project: doi:10.1007/s11837-011-0102-6\n\n"
   "@Article{tadmor:elliott:2011,\n"
   " author = {E. B. Tadmor and R. S. Elliott and J. P. Sethna and R. E. Miller "
   "and C. A. Becker},\n"
@@ -85,7 +85,7 @@ static constexpr const char *const cite_openkim =
   "}\n\n";
 
 static constexpr const char *const cite_openkim_query =
-  "OpenKIM query: https://doi.org/10.1063/5.0014267\n\n"
+  "OpenKIM query: doi:10.1063/5.0014267\n\n"
   "@Article{karls:bierbaum:2020,\n"
   " author = {D. S. Karls and M. Bierbaum and A. A. Alemi and R. S. Elliott "
   "and J. P. Sethna and E. B. Tadmor},\n"

src/KIM/pair_kim.cpp

@@ -77,7 +77,7 @@
 using namespace LAMMPS_NS;
 
 static constexpr const char *const cite_openkim =
-  "OpenKIM: https://doi.org/10.1007/s11837-011-0102-6\n\n"
+  "OpenKIM Project: doi:10.1007/s11837-011-0102-6\n\n"
   "@Article{tadmor:elliott:2011,\n"
   " author = {E. B. Tadmor and R. S. Elliott and J. P. Sethna and R. E. Miller "
   "and C. A. Becker},\n"

src/KOKKOS/comm_kokkos.cpp

@@ -125,19 +125,19 @@ void CommKokkos::init()
   if (force->pair && (force->pair->execution_space == Host))
     check_reverse += force->pair->comm_reverse;
 
-  for (int i = 0; i < modify->nfix; i++) {
+  for (const auto &fix : modify->get_fix_list()) {
-    check_forward += modify->fix[i]->comm_forward;
+    check_forward += fix->comm_forward;
-    check_reverse += modify->fix[i]->comm_reverse;
+    check_reverse += fix->comm_reverse;
   }
 
-  for (int i = 0; i < modify->ncompute; i++) {
+  for (const auto &compute : modify->get_compute_list()) {
-    check_forward += modify->compute[i]->comm_forward;
+    check_forward += compute->comm_forward;
-    check_reverse += modify->compute[i]->comm_reverse;
+    check_reverse += compute->comm_reverse;
   }
 
-  for (int i = 0; i < output->ndump; i++) {
+  for (const auto &dump : output->get_dump_list()) {
-    check_forward += output->dump[i]->comm_forward;
+    check_forward += dump->comm_forward;
-    check_reverse += output->dump[i]->comm_reverse;
+    check_reverse += dump->comm_reverse;
   }
 
   if (force->newton == 0) check_reverse = 0;

src/KOKKOS/fix_shake_kokkos.cpp

@@ -12,12 +12,8 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
-#include <mpi.h>
-#include <cmath>
-#include <cstdlib>
-#include <cstring>
-#include <cstdio>
 #include "fix_shake_kokkos.h"
+
 #include "fix_rattle.h"
 #include "atom_kokkos.h"
 #include "atom_vec.h"
@@ -38,6 +34,9 @@
 #include "kokkos.h"
 #include "atom_masks.h"
 
+#include <cmath>
+#include <cstring>
+
 using namespace LAMMPS_NS;
 using namespace FixConst;
 using namespace MathConst;
@@ -181,6 +180,16 @@ void FixShakeKokkos<DeviceType>::init()
 }
 
 
+/* ----------------------------------------------------------------------
+   run setup for minimization.
+------------------------------------------------------------------------- */
+
+template<class DeviceType>
+void FixShakeKokkos<DeviceType>::min_setup(int /*vflag*/)
+{
+  error->all(FLERR, "Cannot yet use fix {} during minimization with Kokkos", style);
+}
+
 /* ----------------------------------------------------------------------
    build list of SHAKE clusters to constrain
    if one or more atoms in cluster are on this proc,
@@ -193,6 +202,7 @@ void FixShakeKokkos<DeviceType>::pre_neighbor()
   // local copies of atom quantities
   // used by SHAKE until next re-neighboring
 
+  ebond = 0.0;
   x = atom->x;
   v = atom->v;
   f = atom->f;
@@ -292,7 +302,7 @@ void FixShakeKokkos<DeviceType>::pre_neighbor()
 
   if (h_error_flag() == 1) {
     error->one(FLERR,"Shake atoms missing on proc "
-                                 "{} at step {}",me,update->ntimestep);
+                                 "{} at step {}",comm->me,update->ntimestep);
   }
 }
 
@@ -303,6 +313,7 @@ void FixShakeKokkos<DeviceType>::pre_neighbor()
 template<class DeviceType>
 void FixShakeKokkos<DeviceType>::post_force(int vflag)
 {
+  ebond = 0.0;
   copymode = 1;
 
   d_x = atomKK->k_x.view<DeviceType>();
@@ -341,7 +352,7 @@ void FixShakeKokkos<DeviceType>::post_force(int vflag)
   // communicate results if necessary
 
   unconstrained_update();
-  if (nprocs > 1) comm->forward_comm(this);
+  if (comm->nprocs > 1) comm->forward_comm(this);
   k_xshake.sync<DeviceType>();
 
   // virial setup
@@ -1702,7 +1713,7 @@ void FixShakeKokkos<DeviceType>::correct_coordinates(int vflag) {
 
   double **xtmp = xshake;
   xshake = x;
-  if (nprocs > 1) {
+  if (comm->nprocs > 1) {
     forward_comm_device = 0;
     comm->forward_comm(this);
     forward_comm_device = 1;

src/KOKKOS/fix_shake_kokkos.h

@@ -51,6 +51,7 @@ class FixShakeKokkos : public FixShake, public KokkosBase {
   FixShakeKokkos(class LAMMPS *, int, char **);
   ~FixShakeKokkos() override;
   void init() override;
+  void min_setup(int) override;
   void pre_neighbor() override;
   void post_force(int) override;
 

src/MISC/pair_srp_react.h

@@ -34,7 +34,7 @@ class PairSRPREACT : public PairSRP {
  private:
   char *idbreak;
   char *idcreate;
-  bool bond_create, bond_break;
+  bool bond_break, bond_create;
 };
 }    // namespace LAMMPS_NS
 #endif

src/ML-SNAP/compute_snap.cpp

@@ -539,18 +539,18 @@ void ComputeSnap::compute_array()
         for (int icoeff = 0; icoeff < ncoeff; icoeff++)
           snap[irow][k++] += snaptr->blist[icoeff];
 
-          // quadratic contributions
+        // quadratic contributions
 
-          if (quadraticflag) {
+        if (quadraticflag) {
-            for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+          for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
-              double bveci = snaptr->blist[icoeff];
+            double bveci = snaptr->blist[icoeff];
-              snap[irow][k++] += 0.5*bveci*bveci;
+            snap[irow][k++] += 0.5*bveci*bveci;
-              for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
+            for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
-                double bvecj = snaptr->blist[jcoeff];
+              double bvecj = snaptr->blist[jcoeff];
-                snap[irow][k++] += bveci*bvecj;
+              snap[irow][k++] += bveci*bvecj;
-              }
             }
           }
+        }
 
       } else {
           int k = 3;

src/MPIIO/dump_custom_mpiio.cpp

@@ -231,10 +231,12 @@ void DumpCustomMPIIO::init_style()
   // lo priority = line, medium priority = int/float, hi priority = column
 
   auto words = utils::split_words(format);
-  if ((int) words.size() < nfield) error->all(FLERR, "Dump_modify format line is too short");
+  if ((int) words.size() < nfield)
+    error->all(FLERR, "Dump_modify format line is too short: {}", format);
 
   int i = 0;
   for (const auto &word : words) {
+    if (i >= nfield) break;
     delete[] vformat[i];
 
     if (format_column_user[i])

src/REAXFF/fix_reaxff_species.cpp

@@ -827,12 +827,11 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec)
 
 void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec)
 {
-  int i, j, m, n, k, itype, cid;
+  int i, j, m, n, itype, cid;
   int ndel, ndelone, count, count_tmp;
   int *Nameall;
   int *mask = atom->mask;
   double localmass, totalmass;
-  double **spec_atom = f_SPECBOND->array_atom;
   std::string species_str;
 
   AtomVec *avec = atom->avec;

src/REPLICA/hyper.cpp

@@ -39,7 +39,7 @@ enum{NOHYPER,GLOBAL,LOCAL};
 
 /* ---------------------------------------------------------------------- */
 
-Hyper::Hyper(LAMMPS *lmp) : Command(lmp), dumplist(nullptr) {}
+Hyper::Hyper(LAMMPS *_lmp) : Command(_lmp) {}
 
 /* ----------------------------------------------------------------------
    perform hyperdynamics simulation
@@ -177,9 +177,7 @@ void Hyper::command(int narg, char **arg)
 
   fix_event->store_state_quench();
   quench(1);
-  if (dumpflag)
+  if (dumpflag) for (auto idump : dumplist) idump->write();
-    for (int idump = 0; idump < ndump; idump++)
-      output->dump[dumplist[idump]]->write();
   fix_event->store_event();
   if (hyperenable) fix_hyper->build_bond_list(0);
   fix_event->restore_state_quench();
@@ -208,9 +206,7 @@ void Hyper::command(int narg, char **arg)
       nevent++;
       nevent_atoms += ecount;
 
-      if (dumpflag)
+      if (dumpflag) for (auto idump : dumplist) idump->write();
-        for (int idump = 0; idump < ndump; idump++)
-          output->dump[dumplist[idump]]->write();
       fix_event->store_event();
       if (hyperenable) fix_hyper->build_bond_list(ecount);
 
@@ -345,7 +341,6 @@ void Hyper::command(int narg, char **arg)
 
   delete [] id_fix;
   delete [] id_compute;
-  memory->destroy(dumplist);
   delete finish;
   modify->delete_fix("hyper_event");
 
@@ -444,8 +439,6 @@ void Hyper::options(int narg, char **arg)
   maxiter = 40;
   maxeval = 50;
   dumpflag = 0;
-  ndump = 0;
-  dumplist = nullptr;
   rebond = 0;
 
   int iarg = 0;
@@ -462,11 +455,9 @@ void Hyper::options(int narg, char **arg)
     } else if (strcmp(arg[iarg],"dump") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal hyper command");
       dumpflag = 1;
-      int idump = output->find_dump(arg[iarg+1]);
+      auto idump = output->get_dump_by_id(arg[iarg+1]);
-      if (idump < 0)
+      if (!idump) error->all(FLERR,"Dump ID {} in hyper command does not exist", arg[iarg+1]);
-        error->all(FLERR,"Dump ID in hyper command does not exist");
+      dumplist.emplace_back(idump);
-      memory->grow(dumplist,ndump+1,"hyper:dumplist");
-      dumplist[ndump++] = idump;
       iarg += 2;
 
     } else if (strcmp(arg[iarg],"rebond") == 0) {

src/REPLICA/hyper.h

@@ -35,8 +35,8 @@ class Hyper : public Command {
   int t_event;
   double etol, ftol;
   int maxiter, maxeval;
-  int stepmode, dumpflag, ndump, rebond;
+  int stepmode, dumpflag, rebond;
-  int *dumplist;
+  std::vector<class Dump *> dumplist;
 
   int neigh_every, neigh_delay, neigh_dist_check;
   int quench_reneighbor;

src/RIGID/fix_rattle.cpp

@@ -81,7 +81,7 @@ FixRattle::~FixRattle()
 {
   memory->destroy(vp);
 
-  if (RATTLE_DEBUG) {
+#if RATTLE_DEBUG
 
     // communicate maximum distance error
 
@@ -91,13 +91,11 @@ FixRattle::~FixRattle()
     MPI_Reduce(&derr_max, &global_derr_max, 1 , MPI_DOUBLE, MPI_MAX, 0, world);
     MPI_Reduce(&verr_max, &global_verr_max, 1 , MPI_DOUBLE, MPI_MAX, 0, world);
 
-    MPI_Comm_rank (world, &npid); // Find out process rank
+    if (comm->me == 0 && screen) {
-
-    if (npid == 0 && screen) {
       fprintf(screen, "RATTLE: Maximum overall relative position error ( (r_ij-d_ij)/d_ij ): %.10g\n", global_derr_max);
       fprintf(screen, "RATTLE: Maximum overall absolute velocity error (r_ij * v_ij): %.10g\n", global_verr_max);
     }
-  }
+#endif
 }
 
 /* ---------------------------------------------------------------------- */
@@ -110,7 +108,10 @@ int FixRattle::setmask()
   mask |= POST_FORCE_RESPA;
   mask |= FINAL_INTEGRATE;
   mask |= FINAL_INTEGRATE_RESPA;
-  if (RATTLE_DEBUG) mask |= END_OF_STEP;
+  mask |= MIN_POST_FORCE;
+#if RATTLE_DEBUG
+  mask |= END_OF_STEP;
+#endif
   return mask;
 }
 
@@ -156,7 +157,7 @@ void FixRattle::post_force(int vflag)
 
   // communicate the unconstrained velocities
 
-  if (nprocs > 1) {
+  if (comm->nprocs > 1) {
     comm_mode = VP;
     comm->forward_comm(this);
   }
@@ -188,7 +189,7 @@ void FixRattle::post_force_respa(int vflag, int ilevel, int /*iloop*/)
 
   // communicate the unconstrained velocities
 
-  if (nprocs > 1) {
+  if (comm->nprocs > 1) {
     comm_mode = VP;
     comm->forward_comm(this);
   }
@@ -718,7 +719,7 @@ void FixRattle::unpack_forward_comm(int n, int first, double *buf)
 
 void FixRattle::shake_end_of_step(int vflag) {
 
-  if (nprocs > 1) {
+  if (comm->nprocs > 1) {
     comm_mode = V;
     comm->forward_comm(this);
   }
@@ -738,7 +739,6 @@ void FixRattle::correct_coordinates(int vflag) {
   FixShake::correct_coordinates(vflag);
 }
 
-
 /* ----------------------------------------------------------------------
    Remove the velocity component along any bond.
 ------------------------------------------------------------------------- */
@@ -759,7 +759,7 @@ void FixRattle::correct_velocities() {
 
   // communicate the unconstrained velocities
 
-  if (nprocs > 1) {
+  if (comm->nprocs > 1) {
     comm_mode = VP;
     comm->forward_comm(this);
   }
@@ -769,14 +769,13 @@ void FixRattle::correct_velocities() {
   int m;
   for (int i = 0; i < nlist; i++) {
     m = list[i];
-    if      (shake_flag[m] == 2)        vrattle2(m);
+    if      (shake_flag[m] == 2)  vrattle2(m);
-    else if (shake_flag[m] == 3)        vrattle3(m);
+    else if (shake_flag[m] == 3)  vrattle3(m);
-    else if (shake_flag[m] == 4)        vrattle4(m);
+    else if (shake_flag[m] == 4)  vrattle4(m);
-    else                                vrattle3angle(m);
+    else                          vrattle3angle(m);
   }
 }
 
-
 /* ----------------------------------------------------------------------
    DEBUGGING methods
    The functions below allow you to check whether the
@@ -788,16 +787,16 @@ void FixRattle::correct_velocities() {
 
 void FixRattle::end_of_step()
 {
-  if (nprocs > 1) {
+  if (comm->nprocs > 1) {
-       comm_mode = V;
+    comm_mode = V;
-       comm->forward_comm(this);
+    comm->forward_comm(this);
   }
-  if (!check_constraints(v, RATTLE_TEST_POS, RATTLE_TEST_VEL) && RATTLE_RAISE_ERROR) {
+#if RATTLE_RAISE_ERROR
+  if (!check_constraints(v, RATTLE_TEST_POS, RATTLE_TEST_VEL))
     error->one(FLERR, "Rattle failed ");
-  }
+#endif
 }
 
-
 /* ---------------------------------------------------------------------- */
 
 bool FixRattle::check_constraints(double **v, bool checkr, bool checkv)
@@ -807,12 +806,11 @@ bool FixRattle::check_constraints(double **v, bool checkr, bool checkv)
   int i=0;
   while (i < nlist && ret) {
     m = list[i];
-    if      (shake_flag[m] == 2)     ret =   check2(v,m,checkr,checkv);
+    if      (shake_flag[m] == 2)     ret = check2(v,m,checkr,checkv);
-    else if (shake_flag[m] == 3)     ret =   check3(v,m,checkr,checkv);
+    else if (shake_flag[m] == 3)     ret = check3(v,m,checkr,checkv);
-    else if (shake_flag[m] == 4)     ret =   check4(v,m,checkr,checkv);
+    else if (shake_flag[m] == 4)     ret = check4(v,m,checkr,checkv);
-    else                             ret =   check3angle(v,m,checkr,checkv);
+    else                             ret = check3angle(v,m,checkr,checkv);
     i++;
-    if (!RATTLE_RAISE_ERROR)         ret = true;
   }
   return ret;
 }
@@ -834,14 +832,10 @@ bool FixRattle::check2(double **v, int m, bool checkr, bool checkv)
   MathExtra::sub3(v[i1],v[i0],v01);
 
   stat = !(checkr && (fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1) > tol));
-  if (!stat)
+  if (!stat) error->one(FLERR,"Coordinate constraints are not satisfied up to desired tolerance ");
-     error->one(FLERR,"Coordinate constraints are not satisfied "
-                "up to desired tolerance ");
 
   stat = !(checkv && (fabs(MathExtra::dot3(r01,v01)) > tol));
-  if (!stat)
+  if (!stat) error->one(FLERR,"Velocity constraints are not satisfied up to desired tolerance ");
-     error->one(FLERR,"Velocity constraints are not satisfied "
-                "up to desired tolerance ");
   return stat;
 }
 
@@ -871,15 +865,11 @@ bool FixRattle::check3(double **v, int m, bool checkr, bool checkv)
 
   stat = !(checkr && (fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1) > tol ||
                       fabs(sqrt(MathExtra::dot3(r02,r02))-bond2) > tol));
-  if (!stat)
+  if (!stat) error->one(FLERR,"Coordinate constraints are not satisfied up to desired tolerance ");
-     error->one(FLERR,"Coordinate constraints are not satisfied "
-                "up to desired tolerance ");
 
   stat = !(checkv && (fabs(MathExtra::dot3(r01,v01)) > tol ||
                       fabs(MathExtra::dot3(r02,v02)) > tol));
-  if (!stat)
+  if (!stat) error->one(FLERR,"Velocity constraints are not satisfied up to desired tolerance ");
-     error->one(FLERR,"Velocity constraints are not satisfied "
-                "up to desired tolerance ");
   return stat;
 }
 
@@ -914,16 +904,12 @@ bool FixRattle::check4(double **v, int m, bool checkr, bool checkv)
   stat = !(checkr && (fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1) > tol ||
                       fabs(sqrt(MathExtra::dot3(r02,r02))-bond2) > tol ||
                       fabs(sqrt(MathExtra::dot3(r03,r03))-bond3) > tol));
-  if (!stat)
+  if (!stat) error->one(FLERR,"Coordinate constraints are not satisfied up to desired tolerance ");
-     error->one(FLERR,"Coordinate constraints are not satisfied "
-                "up to desired tolerance ");
 
   stat = !(checkv && (fabs(MathExtra::dot3(r01,v01)) > tol ||
                       fabs(MathExtra::dot3(r02,v02)) > tol ||
                       fabs(MathExtra::dot3(r03,v03)) > tol));
-  if (!stat)
+  if (!stat) error->one(FLERR,"Velocity constraints are not satisfied up to desired tolerance ");
-     error->one(FLERR,"Velocity constraints are not satisfied "
-                "up to desired tolerance ");
   return stat;
 }
 
@@ -954,25 +940,19 @@ bool FixRattle::check3angle(double **v, int m, bool checkr, bool checkv)
   MathExtra::sub3(v[i2],v[i0],v02);
   MathExtra::sub3(v[i2],v[i1],v12);
 
-
-
   double db1 = fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1);
   double db2 = fabs(sqrt(MathExtra::dot3(r02,r02))-bond2);
   double db12 = fabs(sqrt(MathExtra::dot3(r12,r12))-bond12);
 
-
+  stat = !(checkr && (db1 > tol || db2 > tol || db12 > tol));
-  stat = !(checkr && (db1 > tol ||
-                      db2 > tol ||
-                      db12 > tol));
 
   if (derr_max < db1/bond1)    derr_max = db1/bond1;
   if (derr_max < db2/bond2)    derr_max = db2/bond2;
   if (derr_max < db12/bond12)  derr_max = db12/bond12;
 
-
+#if RATTLE_RAISE_ERROR
-  if (!stat && RATTLE_RAISE_ERROR)
+  if (!stat) error->one(FLERR,"Coordinate constraints are not satisfied up to desired tolerance ");
-     error->one(FLERR,"Coordinate constraints are not satisfied "
+#endif
-                "up to desired tolerance ");
 
   double dv1 = fabs(MathExtra::dot3(r01,v01));
   double dv2 = fabs(MathExtra::dot3(r02,v02));
@@ -982,16 +962,10 @@ bool FixRattle::check3angle(double **v, int m, bool checkr, bool checkv)
   if (verr_max < dv2)    verr_max = dv2;
   if (verr_max < dv12)   verr_max = dv12;
 
+  stat = !(checkv && (dv1 > tol || dv2 > tol || dv12> tol));
 
-  stat = !(checkv && (dv1 > tol ||
+#if RATTLE_RAISE_ERROR
-                      dv2 > tol ||
+  if (!stat) error->one(FLERR,"Velocity constraints are not satisfied up to desired tolerance!");
-                      dv12> tol));
+#endif
-
-
-  if (!stat && RATTLE_RAISE_ERROR)
-     error->one(FLERR,"Velocity constraints are not satisfied "
-                "up to desired tolerance!");
-
-
   return stat;
 }

src/RIGID/fix_shake.cpp

@@ -41,8 +41,8 @@ using namespace MathConst;
 
 #define RVOUS 1   // 0 for irregular, 1 for all2all
 
-#define BIG 1.0e20
+static constexpr double BIG = 1.0e20;
-#define MASSDELTA 0.1
+static constexpr double MASSDELTA = 0.1;
 
 /* ---------------------------------------------------------------------- */
 
@@ -52,26 +52,31 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
   loop_respa(nullptr), step_respa(nullptr), x(nullptr), v(nullptr), f(nullptr), ftmp(nullptr),
   vtmp(nullptr), mass(nullptr), rmass(nullptr), type(nullptr), shake_flag(nullptr),
   shake_atom(nullptr), shake_type(nullptr), xshake(nullptr), nshake(nullptr), list(nullptr),
-  b_count(nullptr), b_count_all(nullptr), b_atom(nullptr), b_atom_all(nullptr), b_ave(nullptr), b_max(nullptr),
+  b_count(nullptr), b_count_all(nullptr), b_ave(nullptr), b_max(nullptr), b_min(nullptr),
-  b_min(nullptr), b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), a_count(nullptr),
+  b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), a_count(nullptr),
   a_count_all(nullptr), a_ave(nullptr), a_max(nullptr), a_min(nullptr), a_ave_all(nullptr),
   a_max_all(nullptr), a_min_all(nullptr), atommols(nullptr), onemols(nullptr)
 {
-  MPI_Comm_rank(world,&me);
+  energy_global_flag = energy_peratom_flag = 1;
-  MPI_Comm_size(world,&nprocs);
-
   virial_global_flag = virial_peratom_flag = 1;
-  thermo_virial = 1;
+  thermo_energy = thermo_virial = 1;
   create_attribute = 1;
   dof_flag = 1;
+  scalar_flag = 1;
   stores_ids = 1;
   centroidstressflag = CENTROID_AVAIL;
+  next_output = -1;
+
+  // to avoid uninitialized access
+  vflag_post_force = 0;
+  eflag_pre_reverse = 0;
+  ebond = 0.0;
 
   // error check
 
   molecular = atom->molecular;
   if (molecular == Atom::ATOMIC)
-    error->all(FLERR,"Cannot use fix shake with non-molecular system");
+    error->all(FLERR,"Cannot use fix {} with non-molecular system", style);
 
   // perform initial allocation of atom-based arrays
   // register with Atom class
@@ -92,8 +97,9 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
   comm_forward = 3;
 
   // parse SHAKE args
+  auto mystyle = fmt::format("fix {}",style);
 
-  if (narg < 8) error->all(FLERR,"Illegal fix shake command");
+  if (narg < 8) utils::missing_cmd_args(FLERR,mystyle, error);
 
   tolerance = utils::numeric(FLERR,arg[3],false,lmp);
   max_iter = utils::inumeric(FLERR,arg[4],false,lmp);
@@ -133,49 +139,55 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
     else if (mode == 'b') {
       int i = utils::inumeric(FLERR,arg[next],false,lmp);
       if (i < 1 || i > atom->nbondtypes)
-        error->all(FLERR,"Invalid bond type index for fix shake");
+        error->all(FLERR,"Invalid bond type index for {}", mystyle);
       bond_flag[i] = 1;
 
     } else if (mode == 'a') {
       int i = utils::inumeric(FLERR,arg[next],false,lmp);
       if (i < 1 || i > atom->nangletypes)
-        error->all(FLERR,"Invalid angle type index for fix shake");
+        error->all(FLERR,"Invalid angle type index for {}", mystyle);
       angle_flag[i] = 1;
 
     } else if (mode == 't') {
       int i = utils::inumeric(FLERR,arg[next],false,lmp);
       if (i < 1 || i > atom->ntypes)
-        error->all(FLERR,"Invalid atom type index for fix shake");
+        error->all(FLERR,"Invalid atom type index for {}", mystyle);
       type_flag[i] = 1;
 
     } else if (mode == 'm') {
       double massone = utils::numeric(FLERR,arg[next],false,lmp);
-      if (massone == 0.0) error->all(FLERR,"Invalid atom mass for fix shake");
+      if (massone == 0.0) error->all(FLERR,"Invalid atom mass for {}", mystyle);
       if (nmass == atom->ntypes)
-        error->all(FLERR,"Too many masses for fix shake");
+        error->all(FLERR,"Too many masses for {}", mystyle);
       mass_list[nmass++] = massone;
 
-    } else error->all(FLERR,"Illegal fix shake command");
+    } else error->all(FLERR,"Unknown {} command option: {}", mystyle, arg[next]);
     next++;
   }
 
   // parse optional args
 
   onemols = nullptr;
+  kbond = 1.0e6*force->boltz;
 
   int iarg = next;
   while (iarg < narg) {
-    if (strcmp(arg[next],"mol") == 0) {
+    if (strcmp(arg[iarg],"mol") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal fix shake command");
+      if (iarg+2 > narg) utils::missing_cmd_args(FLERR,mystyle+" mol",error);
       int imol = atom->find_molecule(arg[iarg+1]);
       if (imol == -1)
-        error->all(FLERR,"Molecule template ID for fix shake does not exist");
+        error->all(FLERR,"Molecule template ID {} for {} does not exist", mystyle, arg[iarg+1]);
-      if (atom->molecules[imol]->nset > 1 && comm->me == 0)
+      if ((atom->molecules[imol]->nset > 1) && (comm->me == 0))
-        error->warning(FLERR,"Molecule template for fix shake has multiple molecules");
+        error->warning(FLERR,"Molecule template for {} has multiple molecules", mystyle);
       onemols = &atom->molecules[imol];
       nmol = onemols[0]->nset;
       iarg += 2;
-    } else error->all(FLERR,"Illegal fix shake command");
+    } else if (strcmp(arg[iarg],"kbond") == 0) {
+      if (iarg+2 > narg) utils::missing_cmd_args(FLERR,mystyle+" kbond",error);
+      kbond = utils::numeric(FLERR, arg[iarg+1], false, lmp);
+      if (kbond < 0) error->all(FLERR,"Illegal {} kbond value {}. Must be >= 0.0", mystyle, kbond);
+      iarg += 2;
+    } else error->all(FLERR,"Unknown {} command option: {}", mystyle, arg[iarg]);
   }
 
   // error check for Molecule template
@@ -183,7 +195,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
   if (onemols) {
     for (int i = 0; i < nmol; i++)
       if (onemols[i]->shakeflag == 0)
-        error->all(FLERR,"Fix shake molecule template must have shake info");
+        error->all(FLERR,"Fix {} molecule template must have shake info", style);
   }
 
   // allocate bond and angle distance arrays, indexed from 1 to n
@@ -197,8 +209,6 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
     int nb = atom->nbondtypes + 1;
     b_count = new int[nb];
     b_count_all = new int[nb];
-    b_atom = new int[nb];
-    b_atom_all = new int[nb];
     b_ave = new double[nb];
     b_ave_all = new double[nb];
     b_max = new double[nb];
@@ -291,8 +301,6 @@ FixShake::~FixShake()
   if (output_every) {
     delete[] b_count;
     delete[] b_count_all;
-    delete[] b_atom;
-    delete[] b_atom_all;
     delete[] b_ave;
     delete[] b_ave_all;
     delete[] b_max;
@@ -321,6 +329,9 @@ int FixShake::setmask()
   mask |= PRE_NEIGHBOR;
   mask |= POST_FORCE;
   mask |= POST_FORCE_RESPA;
+  mask |= MIN_PRE_REVERSE;
+  mask |= MIN_POST_FORCE;
+  mask |= POST_RUN;
   return mask;
 }
 
@@ -335,28 +346,24 @@ void FixShake::init()
   double rsq,angle;
 
   // error if more than one shake fix
+  auto pattern = fmt::format("^{}",style);
 
-  int count = 0;
+  if (modify->get_fix_by_style(pattern).size() > 1)
-  for (i = 0; i < modify->nfix; i++)
+    error->all(FLERR,"More than one fix {} instance",style);
-    if (strcmp(modify->fix[i]->style,"shake") == 0) count++;
-  if (count > 1) error->all(FLERR,"More than one fix shake");
 
   // cannot use with minimization since SHAKE turns off bonds
   // that should contribute to potential energy
 
-  if (update->whichflag == 2)
+  if ((comm->me == 0) && (update->whichflag == 2))
-    error->all(FLERR,"Fix shake cannot be used with minimization");
+    error->warning(FLERR,"Using fix {} with minimization.\n  Substituting constraints with "
+                   "harmonic restraint forces using kbond={:.4g}", style, kbond);
 
-  // error if npt,nph fix comes before shake fix
+  // error if a fix changing the box comes before shake fix
-
+  bool boxflag = false;
-  for (i = 0; i < modify->nfix; i++) {
+  for (auto ifix : modify->get_fix_list()) {
-    if (strcmp(modify->fix[i]->style,"npt") == 0) break;
+   if (boxflag && utils::strmatch(ifix->style,pattern))
-    if (strcmp(modify->fix[i]->style,"nph") == 0) break;
+     error->all(FLERR,"Fix {} must come before any box changing fix", style);
-  }
+    if (ifix->box_change) boxflag = true;
-  if (i < modify->nfix) {
-    for (int j = i; j < modify->nfix; j++)
-      if (strcmp(modify->fix[j]->style,"shake") == 0)
-        error->all(FLERR,"Shake fix must come before NPT/NPH fix");
   }
 
   // if rRESPA, find associated fix that must exist
@@ -379,7 +386,7 @@ void FixShake::init()
   // set equilibrium bond distances
 
   if (force->bond == nullptr)
-    error->all(FLERR,"Bond potential must be defined for SHAKE");
+    error->all(FLERR,"Bond style must be defined for fix {}",style);
   for (i = 1; i <= atom->nbondtypes; i++)
     bond_distance[i] = force->bond->equilibrium_distance(i);
 
@@ -390,7 +397,7 @@ void FixShake::init()
   for (i = 1; i <= atom->nangletypes; i++) {
     if (angle_flag[i] == 0) continue;
     if (force->angle == nullptr)
-      error->all(FLERR,"Angle potential must be defined for SHAKE");
+      error->all(FLERR,"Angle style must be defined for fix {}",style);
 
     // scan all atoms for a SHAKE angle cluster
     // extract bond types for the 2 bonds in the cluster
@@ -417,7 +424,7 @@ void FixShake::init()
     // error check for any bond types that are not the same
 
     MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_MAX,world);
-    if (flag_all) error->all(FLERR,"Shake angles have different bond types");
+    if (flag_all) error->all(FLERR,"Fix {} angles have different bond types", style);
 
     // insure all procs have bond types
 
@@ -494,6 +501,23 @@ void FixShake::setup(int vflag)
   shake_end_of_step(vflag);
 }
 
+/* ----------------------------------------------------------------------
+   during minimization fix SHAKE adds strong bond forces
+------------------------------------------------------------------------- */
+
+void FixShake::min_setup(int vflag)
+{
+  pre_neighbor();
+  min_post_force(vflag);
+}
+
+/* --------------------------------------------------------------------- */
+
+void FixShake::setup_pre_reverse(int eflag, int vflag)
+{
+  min_pre_reverse(eflag,vflag);
+}
+
 /* ----------------------------------------------------------------------
    build list of SHAKE clusters to constrain
    if one or more atoms in cluster are on this proc,
@@ -502,6 +526,7 @@ void FixShake::setup(int vflag)
 
 void FixShake::pre_neighbor()
 {
+  ebond = 0.0;
   int atom1,atom2,atom3,atom4;
 
   // local copies of atom quantities
@@ -533,19 +558,17 @@ void FixShake::pre_neighbor()
         atom1 = atom->map(shake_atom[i][0]);
         atom2 = atom->map(shake_atom[i][1]);
         if (atom1 == -1 || atom2 == -1)
-          error->one(FLERR,"Shake atoms {} {} missing on proc "
+          error->one(FLERR,"Shake atoms {} {} missing on proc {} at step {}",shake_atom[i][0],
-                                       "{} at step {}",shake_atom[i][0],
+                     shake_atom[i][1],comm->me,update->ntimestep);
-                                       shake_atom[i][1],me,update->ntimestep);
         if (i <= atom1 && i <= atom2) list[nlist++] = i;
       } else if (shake_flag[i] % 2 == 1) {
         atom1 = atom->map(shake_atom[i][0]);
         atom2 = atom->map(shake_atom[i][1]);
         atom3 = atom->map(shake_atom[i][2]);
         if (atom1 == -1 || atom2 == -1 || atom3 == -1)
-          error->one(FLERR,"Shake atoms {} {} {} missing on proc "
+          error->one(FLERR,"Shake atoms {} {} {} missing on proc {} at step {}",shake_atom[i][0],
-                                       "{} at step {}",shake_atom[i][0],
                                        shake_atom[i][1],shake_atom[i][2],
-                                       me,update->ntimestep);
+                                       comm->me,update->ntimestep);
         if (i <= atom1 && i <= atom2 && i <= atom3) list[nlist++] = i;
       } else {
         atom1 = atom->map(shake_atom[i][0]);
@@ -553,10 +576,9 @@ void FixShake::pre_neighbor()
         atom3 = atom->map(shake_atom[i][2]);
         atom4 = atom->map(shake_atom[i][3]);
         if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1)
-          error->one(FLERR,"Shake atoms {} {} {} {} missing on "
+          error->one(FLERR,"Shake atoms {} {} {} {} missing on proc {} at step {}",shake_atom[i][0],
-                                       "proc {} at step {}",shake_atom[i][0],
                                        shake_atom[i][1],shake_atom[i][2],
-                                       shake_atom[i][3],me,update->ntimestep);
+                                       shake_atom[i][3],comm->me,update->ntimestep);
         if (i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4)
           list[nlist++] = i;
       }
@@ -575,11 +597,13 @@ void FixShake::post_force(int vflag)
   // communicate results if necessary
 
   unconstrained_update();
-  if (nprocs > 1) comm->forward_comm(this);
+  if (comm->nprocs > 1) comm->forward_comm(this);
 
   // virial setup
 
-  v_init(vflag);
+  int eflag = eflag_pre_reverse;
+  ev_init(eflag, vflag);
+  ebond = 0.0;
 
   // loop over clusters to add constraint forces
 
@@ -619,7 +643,7 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop)
   // communicate results if necessary
 
   unconstrained_update_respa(ilevel);
-  if (nprocs > 1) comm->forward_comm(this);
+  if (comm->nprocs > 1) comm->forward_comm(this);
 
   // virial setup only needed on last iteration of innermost level
   //   and if pressure is requested
@@ -644,6 +668,59 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop)
   vflag_post_force = vflag;
 }
 
+/* ----------------------------------------------------------------------
+   store eflag so it can be used in min_post_force
+------------------------------------------------------------------------- */
+
+void FixShake::min_pre_reverse(int eflag, int /*vflag*/)
+{
+  eflag_pre_reverse = eflag;
+}
+
+/* ----------------------------------------------------------------------
+   substitute shake constraints with very strong bonds
+------------------------------------------------------------------------- */
+
+void FixShake::min_post_force(int vflag)
+{
+  if (output_every) {
+    bigint ntimestep = update->ntimestep;
+    if (next_output == ntimestep) stats();
+
+    next_output = ntimestep + output_every;
+    if (ntimestep % output_every != 0)
+      next_output = (ntimestep/output_every)*output_every + output_every;
+  } else next_output = -1;
+
+  int eflag = eflag_pre_reverse;
+  ev_init(eflag, vflag);
+
+  x = atom->x;
+  f = atom->f;
+  nlocal = atom->nlocal;
+  ebond = 0.0;
+
+  // loop over shake clusters to add restraint forces
+
+  for (int i = 0; i < nlist; i++) {
+    int m = list[i];
+    if (shake_flag[m] == 2) {
+      bond_force(shake_atom[m][0], shake_atom[m][1], bond_distance[shake_type[m][0]]);
+    } else if (shake_flag[m] == 3) {
+      bond_force(shake_atom[m][0], shake_atom[m][1], bond_distance[shake_type[m][0]]);
+      bond_force(shake_atom[m][0], shake_atom[m][2], bond_distance[shake_type[m][1]]);
+    } else if (shake_flag[m] == 4) {
+      bond_force(shake_atom[m][0], shake_atom[m][1], bond_distance[shake_type[m][0]]);
+      bond_force(shake_atom[m][0], shake_atom[m][2], bond_distance[shake_type[m][1]]);
+      bond_force(shake_atom[m][0], shake_atom[m][3], bond_distance[shake_type[m][2]]);
+    } else {
+      bond_force(shake_atom[m][0], shake_atom[m][1], bond_distance[shake_type[m][0]]);
+      bond_force(shake_atom[m][0], shake_atom[m][2], bond_distance[shake_type[m][1]]);
+      bond_force(shake_atom[m][1], shake_atom[m][2], angle_distance[shake_type[m][2]]);
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    count # of degrees-of-freedom removed by SHAKE for atoms in igroup
 ------------------------------------------------------------------------- */
@@ -690,10 +767,7 @@ void FixShake::find_clusters()
   tagint tagprev;
   double massone;
 
-  if ((me == 0) && screen) {
+  if (comm->me == 0) utils::logmesg(lmp, "Finding {} clusters ...\n",utils::uppercase(style));
-    if (!rattle) fputs("Finding SHAKE clusters ...\n",screen);
-    else fputs("Finding RATTLE clusters ...\n",screen);
-  }
   atommols = atom->avec->onemols;
 
   tagint *tag = atom->tag;
@@ -805,7 +879,7 @@ void FixShake::find_clusters()
     }
 
   MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world);
-  if (flag_all) error->all(FLERR,"Did not find fix shake partner info");
+  if (flag_all) error->all(FLERR,"Did not find fix {} partner info", style);
 
   // -----------------------------------------------------
   // identify SHAKEable bonds
@@ -1015,7 +1089,7 @@ void FixShake::find_clusters()
   tmp = count4;
   MPI_Allreduce(&tmp,&count4,1,MPI_INT,MPI_SUM,world);
 
-  if (me == 0) {
+  if (comm->me == 0) {
     utils::logmesg(lmp,"{:>8} = # of size 2 clusters\n"
                    "{:>8} = # of size 3 clusters\n"
                    "{:>8} = # of size 4 clusters\n"
@@ -1043,8 +1117,8 @@ void FixShake::atom_owners()
   // one datum for each owned atom: datum = owning proc, atomID
 
   for (int i = 0; i < nlocal; i++) {
-    proclist[i] = tag[i] % nprocs;
+    proclist[i] = tag[i] % comm->nprocs;
-    idbuf[i].me = me;
+    idbuf[i].me = comm->me;
     idbuf[i].atomID = tag[i];
   }
 
@@ -1089,7 +1163,7 @@ void FixShake::partner_info(int *npartner, tagint **partner_tag,
   // set values in 4 partner arrays for all partner atoms I own
   // also setup input buf to rendezvous comm
   // input datums = pair of bonded atoms where I do not own partner
-  // owning proc for each datum = partner_tag % nprocs
+  // owning proc for each datum = partner_tag % comm->nprocs
   // datum: atomID = partner_tag (off-proc), partnerID = tag (on-proc)
   //        4 values for my owned atom
 
@@ -1127,7 +1201,7 @@ void FixShake::partner_info(int *npartner, tagint **partner_tag,
         }
 
       } else {
-        proclist[nsend] = partner_tag[i][j] % nprocs;
+        proclist[nsend] = partner_tag[i][j] % comm->nprocs;
         inbuf[nsend].atomID = partner_tag[i][j];
         inbuf[nsend].partnerID = tag[i];
         inbuf[nsend].mask = mask[i];
@@ -1217,7 +1291,7 @@ void FixShake::nshake_info(int *npartner, tagint **partner_tag,
   // set partner_nshake for all partner atoms I own
   // also setup input buf to rendezvous comm
   // input datums = pair of bonded atoms where I do not own partner
-  // owning proc for each datum = partner_tag % nprocs
+  // owning proc for each datum = partner_tag % comm->nprocs
   // datum: atomID = partner_tag (off-proc), partnerID = tag (on-proc)
   //        nshake value for my owned atom
 
@@ -1231,7 +1305,7 @@ void FixShake::nshake_info(int *npartner, tagint **partner_tag,
       if (m >= 0 && m < nlocal) {
         partner_nshake[i][j] = nshake[m];
       } else {
-        proclist[nsend] = partner_tag[i][j] % nprocs;
+        proclist[nsend] = partner_tag[i][j] % comm->nprocs;
         inbuf[nsend].atomID = partner_tag[i][j];
         inbuf[nsend].partnerID = tag[i];
         inbuf[nsend].nshake = nshake[i];
@@ -1295,7 +1369,7 @@ void FixShake::shake_info(int *npartner, tagint **partner_tag,
   // set 3 shake arrays for all partner atoms I own
   // also setup input buf to rendezvous comm
   // input datums = partner atom where I do not own partner
-  // owning proc for each datum = partner_tag % nprocs
+  // owning proc for each datum = partner_tag % comm->nprocs
   // datum: atomID = partner_tag (off-proc)
   //        values in 3 shake arrays
 
@@ -1317,7 +1391,7 @@ void FixShake::shake_info(int *npartner, tagint **partner_tag,
         shake_type[m][2] = shake_type[i][2];
 
       } else {
-        proclist[nsend] = partner_tag[i][j] % nprocs;
+        proclist[nsend] = partner_tag[i][j] % comm->nprocs;
         inbuf[nsend].atomID = partner_tag[i][j];
         inbuf[nsend].shake_flag = shake_flag[i];
         inbuf[nsend].shake_atom[0] = shake_atom[i][0];
@@ -2451,13 +2525,66 @@ void FixShake::shake3angle(int m)
   }
 }
 
+/* ----------------------------------------------------------------------
+   apply bond force for minimization
+------------------------------------------------------------------------- */
+
+void FixShake::bond_force(tagint id1, tagint id2, double length)
+{
+  int i1 = atom->map(id1);
+  int i2 = atom->map(id2);
+
+  if ((i1 < 0) || (i2 < 0)) return;
+
+  // distance vec between atoms, with PBC
+
+  double delx = x[i1][0] - x[i2][0];
+  double dely = x[i1][1] - x[i2][1];
+  double delz = x[i1][2] - x[i2][2];
+  domain->minimum_image(delx, dely, delz);
+
+  // compute and apply force
+
+  const double r = sqrt(delx * delx + dely * dely + delz * delz);
+  const double dr = r - length;
+  const double rk = kbond * dr;
+  const double fbond = (r > 0.0) ? -2.0 * rk / r : 0.0;
+  const double eb = rk*dr;
+  int list[2];
+  int nlist = 0;
+
+  if (i1 < nlocal) {
+    f[i1][0] += delx * fbond;
+    f[i1][1] += dely * fbond;
+    f[i1][2] += delz * fbond;
+    list[nlist++] = i1;
+    ebond += 0.5*eb;
+  }
+  if (i2 < nlocal) {
+    f[i2][0] -= delx * fbond;
+    f[i2][1] -= dely * fbond;
+    f[i2][2] -= delz * fbond;
+    list[nlist++] = i2;
+    ebond += 0.5*eb;
+  }
+  if (evflag) {
+    double v[6];
+    v[0] = 0.5 * delx * delx * fbond;
+    v[1] = 0.5 * dely * dely * fbond;
+    v[2] = 0.5 * delz * delz * fbond;
+    v[3] = 0.5 * delx * dely * fbond;
+    v[4] = 0.5 * delx * delz * fbond;
+    v[5] = 0.5 * dely * delz * fbond;
+    ev_tally(nlist, list, 2.0, eb, v);
+  }
+}
+
 /* ----------------------------------------------------------------------
    print-out bond & angle statistics
 ------------------------------------------------------------------------- */
 
 void FixShake::stats()
 {
-  int i,j,m,n,iatom,jatom,katom;
   double delx,dely,delz;
   double r,r1,r2,r3,angle;
 
@@ -2466,13 +2593,12 @@ void FixShake::stats()
   int nb = atom->nbondtypes + 1;
   int na = atom->nangletypes + 1;
 
-  for (i = 0; i < nb; i++) {
+  for (int i = 0; i < nb; i++) {
     b_count[i] = 0;
     b_ave[i] = b_max[i] = 0.0;
     b_min[i] = BIG;
-    b_atom[i] = -1;
   }
-  for (i = 0; i < na; i++) {
+  for (int i = 0; i < na; i++) {
     a_count[i] = 0;
     a_ave[i] = a_max[i] = 0.0;
     a_min[i] = BIG;
@@ -2484,25 +2610,26 @@ void FixShake::stats()
   double **x = atom->x;
   int nlocal = atom->nlocal;
 
-  for (i = 0; i < nlocal; i++) {
+  for (int ii = 0; ii < nlist; ++ii) {
-    if (shake_flag[i] == 0) continue;
+    int i = list[ii];
+    int n = shake_flag[i];
+    if (n == 0) continue;
 
     // bond stats
 
-    n = shake_flag[i];
     if (n == 1) n = 3;
-    iatom = atom->map(shake_atom[i][0]);
+    int iatom = atom->map(shake_atom[i][0]);
-    for (j = 1; j < n; j++) {
+    for (int j = 1; j < n; j++) {
-      jatom = atom->map(shake_atom[i][j]);
+      int jatom = atom->map(shake_atom[i][j]);
+      if (jatom >= nlocal) continue;
       delx = x[iatom][0] - x[jatom][0];
       dely = x[iatom][1] - x[jatom][1];
       delz = x[iatom][2] - x[jatom][2];
       domain->minimum_image(delx,dely,delz);
-      r = sqrt(delx*delx + dely*dely + delz*delz);
 
-      m = shake_type[i][j-1];
+      r = sqrt(delx*delx + dely*dely + delz*delz);
+      int m = shake_type[i][j-1];
       b_count[m]++;
-      b_atom[m] = n;
       b_ave[m] += r;
       b_max[m] = MAX(b_max[m],r);
       b_min[m] = MIN(b_min[m],r);
@@ -2511,9 +2638,13 @@ void FixShake::stats()
     // angle stats
 
     if (shake_flag[i] == 1) {
-      iatom = atom->map(shake_atom[i][0]);
+      int iatom = atom->map(shake_atom[i][0]);
-      jatom = atom->map(shake_atom[i][1]);
+      int jatom = atom->map(shake_atom[i][1]);
-      katom = atom->map(shake_atom[i][2]);
+      int katom = atom->map(shake_atom[i][2]);
+      int n = 0;
+      if (iatom < nlocal) ++n;
+      if (jatom < nlocal) ++n;
+      if (katom < nlocal) ++n;
 
       delx = x[iatom][0] - x[jatom][0];
       dely = x[iatom][1] - x[jatom][1];
@@ -2535,9 +2666,9 @@ void FixShake::stats()
 
       angle = acos((r1*r1 + r2*r2 - r3*r3) / (2.0*r1*r2));
       angle *= 180.0/MY_PI;
-      m = shake_type[i][2];
+      int m = shake_type[i][2];
-      a_count[m]++;
+      a_count[m] += n;
-      a_ave[m] += angle;
+      a_ave[m] += n*angle;
       a_max[m] = MAX(a_max[m],angle);
       a_min[m] = MIN(a_min[m],angle);
     }
@@ -2546,7 +2677,6 @@ void FixShake::stats()
   // sum across all procs
 
   MPI_Allreduce(b_count,b_count_all,nb,MPI_INT,MPI_SUM,world);
-  MPI_Allreduce(b_atom,b_atom_all,nb,MPI_INT,MPI_MAX,world);
   MPI_Allreduce(b_ave,b_ave_all,nb,MPI_DOUBLE,MPI_SUM,world);
   MPI_Allreduce(b_max,b_max_all,nb,MPI_DOUBLE,MPI_MAX,world);
   MPI_Allreduce(b_min,b_min_all,nb,MPI_DOUBLE,MPI_MIN,world);
@@ -2558,16 +2688,17 @@ void FixShake::stats()
 
   // print stats only for non-zero counts
 
-  if (me == 0) {
+  if (comm->me == 0) {
     const int width = log10((MAX(MAX(1,nb),na)))+2;
-    auto mesg = fmt::format("SHAKE stats (type/ave/delta/count) on step {}\n", update->ntimestep);
+    auto mesg = fmt::format("{} stats (type/ave/delta/count) on step {}\n",
-    for (i = 1; i < nb; i++) {
+                            utils::uppercase(style), update->ntimestep);
+    for (int i = 1; i < nb; i++) {
       const auto bcnt = b_count_all[i];
       if (bcnt)
         mesg += fmt::format("Bond:  {:>{}d}   {:<9.6} {:<11.6} {:>8d}\n",i,width,
-                            b_ave_all[i]/bcnt,b_max_all[i]-b_min_all[i],bcnt/b_atom_all[i]);
+                            b_ave_all[i]/bcnt,b_max_all[i]-b_min_all[i],bcnt);
     }
-    for (i = 1; i < na; i++) {
+    for (int i = 1; i < na; i++) {
       const auto acnt = a_count_all[i];
       if (acnt)
         mesg += fmt::format("Angle: {:>{}d}   {:<9.6} {:<11.6} {:>8d}\n",i,width,
@@ -3005,6 +3136,26 @@ void *FixShake::extract(const char *str, int &dim)
   return nullptr;
 }
 
+/* ----------------------------------------------------------------------
+   energy due to restraint forces
+------------------------------------------------------------------------- */
+
+double FixShake::compute_scalar()
+{
+  double all;
+  MPI_Allreduce(&ebond, &all, 1, MPI_DOUBLE, MPI_SUM, world);
+  return all;
+}
+
+/*  ----------------------------------------------------------------------
+   print shake stats at the end of a minimization
+------------------------------------------------------------------------- */
+void FixShake::post_run()
+{
+  if ((update->whichflag == 2) && (output_every > 0)) stats();
+}
+
+
 /* ----------------------------------------------------------------------
    add coordinate constraining forces
    this method is called at the end of a timestep
@@ -3109,7 +3260,7 @@ void FixShake::correct_coordinates(int vflag) {
 
   double **xtmp = xshake;
   xshake = x;
-  if (nprocs > 1) {
+  if (comm->nprocs > 1) {
     comm->forward_comm(this);
   }
   xshake = xtmp;

src/RIGID/fix_shake.h

@@ -34,9 +34,14 @@ class FixShake : public Fix {
   int setmask() override;
   void init() override;
   void setup(int) override;
+  void setup_pre_reverse(int, int) override;
+  void min_setup(int) override;
   void pre_neighbor() override;
   void post_force(int) override;
   void post_force_respa(int, int, int) override;
+  void min_pre_reverse(int, int) override;
+  void min_post_force(int) override;
+  void post_run() override;
 
   double memory_usage() override;
   void grow_arrays(int) override;
@@ -57,16 +62,17 @@ class FixShake : public Fix {
   int dof(int) override;
   void reset_dt() override;
   void *extract(const char *, int &) override;
+  double compute_scalar() override;
 
  protected:
   int vflag_post_force;    // store the vflag of last post_force call
+  int eflag_pre_reverse;   // store the eflag of last pre_reverse call
   int respa;               // 0 = vel. Verlet, 1 = respa
-  int me, nprocs;
+  int rattle;              // 0 = SHAKE, 1 = RATTLE
-  int rattle;            // 0 = SHAKE, 1 = RATTLE
+  double tolerance;        // SHAKE tolerance
-  double tolerance;      // SHAKE tolerance
+  int max_iter;            // max # of SHAKE iterations
-  int max_iter;          // max # of SHAKE iterations
+  int output_every;        // SHAKE stat output every so often
-  int output_every;      // SHAKE stat output every so often
+  bigint next_output;      // timestep for next output
-  bigint next_output;    // timestep for next output
 
   // settings from input command
   int *bond_flag, *angle_flag;    // bond/angle types to constrain
@@ -76,6 +82,8 @@ class FixShake : public Fix {
 
   int molecular;                             // copy of atom->molecular
   double *bond_distance, *angle_distance;    // constraint distances
+  double kbond;                              // force constant for restraint
+  double ebond;                              // energy of bond restraints
 
   class FixRespa *fix_respa;    // rRESPA fix needed by SHAKE
   int nlevels_respa;            // copies of needed rRESPA variables
@@ -108,8 +116,7 @@ class FixShake : public Fix {
   int nlist, maxlist;    // size and max-size of list
 
   // stat quantities
-  int *b_count, *b_count_all, *b_atom,
+  int *b_count, *b_count_all;                   // counts for each bond type, atoms in bond cluster
-      *b_atom_all;                              // counts for each bond type, atoms in bond cluster
   double *b_ave, *b_max, *b_min;                // ave/max/min dist for each bond type
   double *b_ave_all, *b_max_all, *b_min_all;    // MPI summing arrays
   int *a_count, *a_count_all;                   // ditto for angle types
@@ -133,6 +140,7 @@ class FixShake : public Fix {
   void shake3(int);
   void shake4(int);
   void shake3angle(int);
+  void bond_force(tagint, tagint, double);
   void stats();
   int bondtype_findset(int, tagint, tagint, int);
   int angletype_findset(int, tagint, tagint, int);

src/SMTBQ/pair_smtbq.cpp

@@ -1232,8 +1232,6 @@ void PairSMTBQ::tabqeq()
   double aCoeff,bCoeff,rcoupe,nang;
 
   int n = atom->ntypes;
-  int nlocal = atom->nlocal;
-  int nghost = atom->nghost;
   nmax = atom->nmax;
   nntype = int((n+1)*n/2);
 
@@ -1242,7 +1240,7 @@ void PairSMTBQ::tabqeq()
 
 #if VERBOSE
   printf ("kmax %d, ds %f, nmax %d\n",kmax,ds,nmax);
-  printf ("nlocal = %d, nghost = %d\n",nlocal,nghost);
+  printf ("nlocal = %d, nghost = %d\n",atom->nlocal,atom->nghost);
   printf ("nntypes %d, kmax %d, rc %f, n %d\n",nntype,kmax,rc,n);
 #endif
 

src/atom.cpp

@@ -1891,11 +1891,10 @@ void Atom::add_molecule(int narg, char **arg)
    return -1 if does not exist
 ------------------------------------------------------------------------- */
 
-int Atom::find_molecule(char *id)
+int Atom::find_molecule(const char *id)
 {
   if (id == nullptr) return -1;
-  int imol;
+  for (int imol = 0; imol < nmolecule; imol++)
-  for (imol = 0; imol < nmolecule; imol++)
     if (strcmp(id,molecules[imol]->id) == 0) return imol;
   return -1;
 }

src/atom.h

@@ -344,7 +344,7 @@ class Atom : protected Pointers {
   int shape_consistency(int, double &, double &, double &);
 
   void add_molecule(int, char **);
-  int find_molecule(char *);
+  int find_molecule(const char *);
   std::vector<Molecule *>get_molecule_by_id(const std::string &);
   void add_molecule_atom(class Molecule *, int, int, tagint);
 

src/comm.cpp

@@ -222,18 +222,18 @@ void Comm::init()
   if (force->bond) maxreverse = MAX(maxreverse,force->bond->comm_reverse);
 
   for (const auto &fix : fix_list) {
-    maxforward = MAX(maxforward,fix->comm_forward);
+    maxforward = MAX(maxforward, fix->comm_forward);
-    maxreverse = MAX(maxreverse,fix->comm_reverse);
+    maxreverse = MAX(maxreverse, fix->comm_reverse);
   }
 
-  for (int i = 0; i < modify->ncompute; i++) {
+  for (const auto &compute : modify->get_compute_list()) {
-    maxforward = MAX(maxforward,modify->compute[i]->comm_forward);
+    maxforward = MAX(maxforward,compute->comm_forward);
-    maxreverse = MAX(maxreverse,modify->compute[i]->comm_reverse);
+    maxreverse = MAX(maxreverse,compute->comm_reverse);
   }
 
-  for (int i = 0; i < output->ndump; i++) {
+  for (const auto &dump: output->get_dump_list()) {
-    maxforward = MAX(maxforward,output->dump[i]->comm_forward);
+    maxforward = MAX(maxforward,dump->comm_forward);
-    maxreverse = MAX(maxreverse,output->dump[i]->comm_reverse);
+    maxreverse = MAX(maxreverse,dump->comm_reverse);
   }
 
   if (force->newton == 0) maxreverse = 0;
@@ -247,8 +247,7 @@ void Comm::init()
   maxexchange_atom = atom->avec->maxexchange;
 
   maxexchange_fix_dynamic = 0;
-  for (const auto &fix : fix_list)
+  for (const auto &fix : fix_list) if (fix->maxexchange_dynamic) maxexchange_fix_dynamic = 1;
-    if (fix->maxexchange_dynamic) maxexchange_fix_dynamic = 1;
 
   if ((mode == Comm::MULTI) && (neighbor->style != Neighbor::MULTI))
     error->all(FLERR,"Cannot use comm mode multi without multi-style neighbor lists");
@@ -270,8 +269,7 @@ void Comm::init()
 void Comm::init_exchange()
 {
   maxexchange_fix = 0;
-  for (const auto &fix : modify->get_fix_list())
+  for (const auto &fix : modify->get_fix_list()) maxexchange_fix += fix->maxexchange;
-    maxexchange_fix += fix->maxexchange;
 
   maxexchange = maxexchange_atom + maxexchange_fix;
   bufextra = maxexchange + BUFEXTRA;

src/compute_centroid_stress_atom.cpp

@@ -37,7 +37,7 @@ enum { NOBIAS, BIAS };
 
 static const char cite_centroid_angle_improper_dihedral[] =
     "compute centroid/stress/atom for angles, impropers and dihedrals: doi:10.1103/PhysRevE.99.051301\n\n"
-    "@article{PhysRevE.99.051301,\n"
+    "@article{Surblys2019,\n"
     " title = {Application of Atomic Stress to Compute Heat Flux via Molecular\n"
     "    Dynamics for Systems With Many-Body Interactions},\n"
     " author = {Surblys, Donatas and Matsubara, Hiroki and Kikugawa, Gota and Ohara, Taku},\n"
@@ -52,7 +52,7 @@ static const char cite_centroid_angle_improper_dihedral[] =
 
 static const char cite_centroid_shake_rigid[] =
     "compute centroid/stress/atom for constrained dynamics: doi:10.1063/5.0070930\n\n"
-    "@article{doi:10.1063/5.0070930,\n"
+    "@article{Surblys2021,\n"
     " author = {Surblys, Donatas and Matsubara, Hiroki and Kikugawa, Gota and Ohara, Taku},\n"
     " journal = {Journal of Applied Physics},\n"
     " title = {Methodology and Meaning of Computing Heat Flux via Atomic Stress in Systems with\n"

src/compute_pe_atom.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -13,30 +12,31 @@
 ------------------------------------------------------------------------- */
 
 #include "compute_pe_atom.h"
-#include <cstring>
+
-#include "atom.h"
-#include "update.h"
-#include "comm.h"
-#include "force.h"
-#include "pair.h"
-#include "bond.h"
 #include "angle.h"
+#include "atom.h"
+#include "bond.h"
+#include "comm.h"
 #include "dihedral.h"
+#include "error.h"
+#include "force.h"
 #include "improper.h"
 #include "kspace.h"
-#include "modify.h"
 #include "memory.h"
-#include "error.h"
+#include "modify.h"
+#include "pair.h"
+#include "update.h"
+
+#include <cstring>
 
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
 ComputePEAtom::ComputePEAtom(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg),
+    Compute(lmp, narg, arg), energy(nullptr)
-  energy(nullptr)
 {
-  if (narg < 3) error->all(FLERR,"Illegal compute pe/atom command");
+  if (narg < 3) error->all(FLERR, "Illegal compute pe/atom command");
 
   peratom_flag = 1;
   size_peratom_cols = 0;
@@ -56,14 +56,22 @@ ComputePEAtom::ComputePEAtom(LAMMPS *lmp, int narg, char **arg) :
     fixflag = 0;
     int iarg = 3;
     while (iarg < narg) {
-      if (strcmp(arg[iarg],"pair") == 0) pairflag = 1;
+      if (strcmp(arg[iarg], "pair") == 0)
-      else if (strcmp(arg[iarg],"bond") == 0) bondflag = 1;
+        pairflag = 1;
-      else if (strcmp(arg[iarg],"angle") == 0) angleflag = 1;
+      else if (strcmp(arg[iarg], "bond") == 0)
-      else if (strcmp(arg[iarg],"dihedral") == 0) dihedralflag = 1;
+        bondflag = 1;
-      else if (strcmp(arg[iarg],"improper") == 0) improperflag = 1;
+      else if (strcmp(arg[iarg], "angle") == 0)
-      else if (strcmp(arg[iarg],"kspace") == 0) kspaceflag = 1;
+        angleflag = 1;
-      else if (strcmp(arg[iarg],"fix") == 0) fixflag = 1;
+      else if (strcmp(arg[iarg], "dihedral") == 0)
-      else error->all(FLERR,"Illegal compute pe/atom command");
+        dihedralflag = 1;
+      else if (strcmp(arg[iarg], "improper") == 0)
+        improperflag = 1;
+      else if (strcmp(arg[iarg], "kspace") == 0)
+        kspaceflag = 1;
+      else if (strcmp(arg[iarg], "fix") == 0)
+        fixflag = 1;
+      else
+        error->all(FLERR, "Illegal compute pe/atom command");
       iarg++;
     }
   }
@@ -86,7 +94,7 @@ void ComputePEAtom::compute_peratom()
 
   invoked_peratom = update->ntimestep;
   if (update->eflag_atom != invoked_peratom)
-    error->all(FLERR,"Per-atom energy was not tallied on needed timestep");
+    error->all(FLERR, "Per-atom energy was not tallied on needed timestep");
 
   // grow local energy array if necessary
   // needs to be atom->nmax in length
@@ -94,7 +102,7 @@ void ComputePEAtom::compute_peratom()
   if (atom->nmax > nmax) {
     memory->destroy(energy);
     nmax = atom->nmax;
-    memory->create(energy,nmax,"pe/atom:energy");
+    memory->create(energy, nmax, "pe/atom:energy");
     vector_atom = energy;
   }
 
@@ -153,13 +161,11 @@ void ComputePEAtom::compute_peratom()
   // add in per-atom contributions from relevant fixes
   // always only for owned atoms, not ghost
 
-  if (fixflag && modify->n_energy_atom)
+  if (fixflag && modify->n_energy_atom) modify->energy_atom(nlocal, energy);
-    modify->energy_atom(nlocal,energy);
 
   // communicate ghost energy between neighbor procs
 
-  if (force->newton || (force->kspace && force->kspace->tip4pflag))
+  if (force->newton || (force->kspace && force->kspace->tip4pflag)) comm->reverse_comm(this);
-    comm->reverse_comm(this);
 
   // zero energy of atoms not in group
   // only do this after comm since ghost contributions must be included
@@ -174,7 +180,7 @@ void ComputePEAtom::compute_peratom()
 
 int ComputePEAtom::pack_reverse_comm(int n, int first, double *buf)
 {
-  int i,m,last;
+  int i, m, last;
 
   m = 0;
   last = first + n;
@@ -186,7 +192,7 @@ int ComputePEAtom::pack_reverse_comm(int n, int first, double *buf)
 
 void ComputePEAtom::unpack_reverse_comm(int n, int *list, double *buf)
 {
-  int i,j,m;
+  int i, j, m;
 
   m = 0;
   for (i = 0; i < n; i++) {
@@ -201,6 +207,6 @@ void ComputePEAtom::unpack_reverse_comm(int n, int *list, double *buf)
 
 double ComputePEAtom::memory_usage()
 {
-  double bytes = (double)nmax * sizeof(double);
+  double bytes = (double) nmax * sizeof(double);
   return bytes;
 }

src/dump_custom.cpp

@@ -262,10 +262,11 @@ void DumpCustom::init_style()
 
   auto words = utils::split_words(format);
   if ((int) words.size() < nfield)
-    error->all(FLERR,"Dump_modify format line is too short");
+    error->all(FLERR,"Dump_modify format line is too short: {}", format);
 
   int i=0;
   for (const auto &word : words) {
+    if (i >= nfield) break;
     delete[] vformat[i];
 
     if (format_column_user[i])

src/dump_local.cpp

@@ -185,10 +185,11 @@ void DumpLocal::init_style()
 
   auto words = utils::split_words(format);
   if ((int) words.size() <  size_one)
-    error->all(FLERR,"Dump_modify format line is too short");
+    error->all(FLERR,"Dump_modify format line is too short: {}", format);
 
   int i=0;
   for (const auto &word : words) {
+    if (i >= size_one) break;
     delete[] vformat[i];
 
     if (format_column_user[i])

src/fix_langevin.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -44,21 +43,20 @@
 using namespace LAMMPS_NS;
 using namespace FixConst;
 
-enum{NOBIAS,BIAS};
+enum { NOBIAS, BIAS };
-enum{CONSTANT,EQUAL,ATOM};
+enum { CONSTANT, EQUAL, ATOM };
 
-#define SINERTIA 0.4          // moment of inertia prefactor for sphere
+#define SINERTIA 0.4    // moment of inertia prefactor for sphere
-#define EINERTIA 0.2          // moment of inertia prefactor for ellipsoid
+#define EINERTIA 0.2    // moment of inertia prefactor for ellipsoid
 
 /* ---------------------------------------------------------------------- */
 
 FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg),
+    Fix(lmp, narg, arg), gjfflag(0), gfactor1(nullptr), gfactor2(nullptr), ratio(nullptr),
-  gjfflag(0), gfactor1(nullptr), gfactor2(nullptr), ratio(nullptr), tstr(nullptr),
+    tstr(nullptr), flangevin(nullptr), tforce(nullptr), franprev(nullptr), lv(nullptr),
-  flangevin(nullptr), tforce(nullptr), franprev(nullptr),
+    id_temp(nullptr), random(nullptr)
-  lv(nullptr), id_temp(nullptr), random(nullptr)
 {
-  if (narg < 7) error->all(FLERR,"Illegal fix langevin command");
+  if (narg < 7) error->all(FLERR, "Illegal fix langevin command");
 
   dynamic_group_allow = 1;
   scalar_flag = 1;
@@ -67,30 +65,30 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
   ecouple_flag = 1;
   nevery = 1;
 
-  if (utils::strmatch(arg[3],"^v_")) {
+  if (utils::strmatch(arg[3], "^v_")) {
-    tstr = utils::strdup(arg[3]+2);
+    tstr = utils::strdup(arg[3] + 2);
   } else {
-    t_start = utils::numeric(FLERR,arg[3],false,lmp);
+    t_start = utils::numeric(FLERR, arg[3], false, lmp);
     t_target = t_start;
     tstyle = CONSTANT;
   }
 
-  t_stop = utils::numeric(FLERR,arg[4],false,lmp);
+  t_stop = utils::numeric(FLERR, arg[4], false, lmp);
-  t_period = utils::numeric(FLERR,arg[5],false,lmp);
+  t_period = utils::numeric(FLERR, arg[5], false, lmp);
-  seed = utils::inumeric(FLERR,arg[6],false,lmp);
+  seed = utils::inumeric(FLERR, arg[6], false, lmp);
 
-  if (t_period <= 0.0) error->all(FLERR,"Fix langevin period must be > 0.0");
+  if (t_period <= 0.0) error->all(FLERR, "Fix langevin period must be > 0.0");
-  if (seed <= 0) error->all(FLERR,"Illegal fix langevin command");
+  if (seed <= 0) error->all(FLERR, "Illegal fix langevin command");
 
   // initialize Marsaglia RNG with processor-unique seed
 
-  random = new RanMars(lmp,seed + comm->me);
+  random = new RanMars(lmp, seed + comm->me);
 
   // allocate per-type arrays for force prefactors
 
-  gfactor1 = new double[atom->ntypes+1];
+  gfactor1 = new double[atom->ntypes + 1];
-  gfactor2 = new double[atom->ntypes+1];
+  gfactor2 = new double[atom->ntypes + 1];
-  ratio = new double[atom->ntypes+1];
+  ratio = new double[atom->ntypes + 1];
 
   // optional args
 
@@ -105,48 +103,48 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
 
   int iarg = 7;
   while (iarg < narg) {
-    if (strcmp(arg[iarg],"angmom") == 0) {
+    if (strcmp(arg[iarg], "angmom") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
+      if (iarg + 2 > narg) error->all(FLERR, "Illegal fix langevin command");
-      if (strcmp(arg[iarg+1],"no") == 0) ascale = 0.0;
+      if (strcmp(arg[iarg + 1], "no") == 0)
-      else ascale = utils::numeric(FLERR,arg[iarg+1],false,lmp);
+        ascale = 0.0;
+      else
+        ascale = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
       iarg += 2;
-    } else if (strcmp(arg[iarg],"gjf") == 0) {
+    } else if (strcmp(arg[iarg], "gjf") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
+      if (iarg + 2 > narg) error->all(FLERR, "Illegal fix langevin command");
-      if (strcmp(arg[iarg+1],"no") == 0) {
+      if (strcmp(arg[iarg + 1], "no") == 0) {
         gjfflag = 0;
         osflag = 0;
-      }
+      } else if (strcmp(arg[iarg + 1], "vfull") == 0) {
-      else if (strcmp(arg[iarg+1],"vfull") == 0) {
         gjfflag = 1;
         osflag = 1;
-      }
+      } else if (strcmp(arg[iarg + 1], "vhalf") == 0) {
-      else if (strcmp(arg[iarg+1],"vhalf") == 0) {
         gjfflag = 1;
         osflag = 0;
-      }
+      } else
-      else error->all(FLERR,"Illegal fix langevin command");
+        error->all(FLERR, "Illegal fix langevin command");
       iarg += 2;
-    } else if (strcmp(arg[iarg],"omega") == 0) {
+    } else if (strcmp(arg[iarg], "omega") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
+      if (iarg + 2 > narg) error->all(FLERR, "Illegal fix langevin command");
-      oflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
+      oflag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
       iarg += 2;
-    } else if (strcmp(arg[iarg],"scale") == 0) {
+    } else if (strcmp(arg[iarg], "scale") == 0) {
-      if (iarg+3 > narg) error->all(FLERR,"Illegal fix langevin command");
+      if (iarg + 3 > narg) error->all(FLERR, "Illegal fix langevin command");
-      int itype = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
+      int itype = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
-      double scale = utils::numeric(FLERR,arg[iarg+2],false,lmp);
+      double scale = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
-      if (itype <= 0 || itype > atom->ntypes)
+      if (itype <= 0 || itype > atom->ntypes) error->all(FLERR, "Illegal fix langevin command");
-        error->all(FLERR,"Illegal fix langevin command");
       ratio[itype] = scale;
       iarg += 3;
-    } else if (strcmp(arg[iarg],"tally") == 0) {
+    } else if (strcmp(arg[iarg], "tally") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
+      if (iarg + 2 > narg) error->all(FLERR, "Illegal fix langevin command");
-      tallyflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
+      tallyflag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
       iarg += 2;
-    } else if (strcmp(arg[iarg],"zero") == 0) {
+    } else if (strcmp(arg[iarg], "zero") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
+      if (iarg + 2 > narg) error->all(FLERR, "Illegal fix langevin command");
-      zeroflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
+      zeroflag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
       iarg += 2;
-    } else error->all(FLERR,"Illegal fix langevin command");
+    } else
+      error->all(FLERR, "Illegal fix langevin command");
   }
 
   // set temperature = nullptr, user can override via fix_modify if wants bias
@@ -175,7 +173,7 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
     FixLangevin::grow_arrays(atom->nmax);
     atom->add_callback(Atom::GROW);
 
-  // initialize franprev to zero
+    // initialize franprev to zero
 
     int nlocal = atom->nlocal;
     for (int i = 0; i < nlocal; i++) {
@@ -194,18 +192,18 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
 FixLangevin::~FixLangevin()
 {
   delete random;
-  delete [] tstr;
+  delete[] tstr;
-  delete [] gfactor1;
+  delete[] gfactor1;
-  delete [] gfactor2;
+  delete[] gfactor2;
-  delete [] ratio;
+  delete[] ratio;
-  delete [] id_temp;
+  delete[] id_temp;
   memory->destroy(flangevin);
   memory->destroy(tforce);
 
   if (gjfflag) {
     memory->destroy(franprev);
     memory->destroy(lv);
-    atom->delete_callback(id,Atom::GROW);
+    atom->delete_callback(id, Atom::GROW);
   }
 }
 
@@ -226,34 +224,38 @@ int FixLangevin::setmask()
 void FixLangevin::init()
 {
   if (gjfflag) {
-    if (t_period*2 == update->dt)
+    if (t_period * 2 == update->dt)
-      error->all(FLERR,"Fix langevin gjf cannot have t_period equal to dt/2");
+      error->all(FLERR, "Fix langevin gjf cannot have t_period equal to dt/2");
 
     // warn if any integrate fix comes after this one
     int before = 1;
     int flag = 0;
     for (int i = 0; i < modify->nfix; i++) {
-      if (strcmp(id,modify->fix[i]->id) == 0) before = 0;
+      auto ifix = modify->get_fix_by_index(i);
-      else if ((modify->fmask[i] && utils::strmatch(modify->fix[i]->style,"^nve")) && before) flag = 1;
+      if (strcmp(id, ifix->id) == 0)
+        before = 0;
+      else if ((modify->fmask[i] && utils::strmatch(ifix->style, "^nve")) && before)
+        flag = 1;
     }
-    if (flag)
+    if (flag) error->all(FLERR, "Fix langevin gjf should come before fix nve");
-      error->all(FLERR,"Fix langevin gjf should come before fix nve");
   }
 
   if (oflag && !atom->sphere_flag)
-    error->all(FLERR,"Fix langevin omega requires atom style sphere");
+    error->all(FLERR, "Fix langevin omega requires atom style sphere");
   if (ascale && !atom->ellipsoid_flag)
-    error->all(FLERR,"Fix langevin angmom requires atom style ellipsoid");
+    error->all(FLERR, "Fix langevin angmom requires atom style ellipsoid");
 
   // check variable
 
   if (tstr) {
     tvar = input->variable->find(tstr);
-    if (tvar < 0)
+    if (tvar < 0) error->all(FLERR, "Variable name {} for fix langevin does not exist", tstr);
-      error->all(FLERR,"Variable name for fix langevin does not exist");
+    if (input->variable->equalstyle(tvar))
-    if (input->variable->equalstyle(tvar)) tstyle = EQUAL;
+      tstyle = EQUAL;
-    else if (input->variable->atomstyle(tvar)) tstyle = ATOM;
+    else if (input->variable->atomstyle(tvar))
-    else error->all(FLERR,"Variable for fix langevin is invalid style");
+      tstyle = ATOM;
+    else
+      error->all(FLERR, "Variable {} for fix langevin is invalid style", tstr);
   }
 
   // if oflag or ascale set, check that all group particles are finite-size
@@ -265,14 +267,12 @@ void FixLangevin::init()
 
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & groupbit)
-        if (radius[i] == 0.0)
+        if (radius[i] == 0.0) error->one(FLERR, "Fix langevin omega requires extended particles");
-          error->one(FLERR,"Fix langevin omega requires extended particles");
   }
 
   if (ascale) {
-    avec = dynamic_cast<AtomVecEllipsoid *>( atom->style_match("ellipsoid"));
+    avec = dynamic_cast<AtomVecEllipsoid *>(atom->style_match("ellipsoid"));
-    if (!avec)
+    if (!avec) error->all(FLERR, "Fix langevin angmom requires atom style ellipsoid");
-      error->all(FLERR,"Fix langevin angmom requires atom style ellipsoid");
 
     int *ellipsoid = atom->ellipsoid;
     int *mask = atom->mask;
@@ -280,8 +280,7 @@ void FixLangevin::init()
 
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & groupbit)
-        if (ellipsoid[i] < 0)
+        if (ellipsoid[i] < 0) error->one(FLERR, "Fix langevin angmom requires extended particles");
-          error->one(FLERR,"Fix langevin angmom requires extended particles");
   }
 
   // set force prefactors
@@ -289,30 +288,30 @@ void FixLangevin::init()
   if (!atom->rmass) {
     for (int i = 1; i <= atom->ntypes; i++) {
       gfactor1[i] = -atom->mass[i] / t_period / force->ftm2v;
+      gfactor2[i] = sqrt(atom->mass[i]) / force->ftm2v;
       if (gjfflag)
-        gfactor2[i] = sqrt(atom->mass[i]) *
+        gfactor2[i] *= sqrt(2.0 * force->boltz / t_period / update->dt / force->mvv2e);
-          sqrt(2.0*force->boltz/t_period/update->dt/force->mvv2e) /
-          force->ftm2v;
       else
-        gfactor2[i] = sqrt(atom->mass[i]) *
+        gfactor2[i] *= sqrt(24.0 * force->boltz / t_period / update->dt / force->mvv2e);
-          sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) /
+      gfactor1[i] *= 1.0 / ratio[i];
-          force->ftm2v;
+      gfactor2[i] *= 1.0 / sqrt(ratio[i]);
-      gfactor1[i] *= 1.0/ratio[i];
-      gfactor2[i] *= 1.0/sqrt(ratio[i]);
     }
   }
 
-  if (temperature && temperature->tempbias) tbiasflag = BIAS;
+  if (temperature && temperature->tempbias)
-  else tbiasflag = NOBIAS;
+    tbiasflag = BIAS;
+  else
+    tbiasflag = NOBIAS;
 
-  if (utils::strmatch(update->integrate_style,"^respa"))
+  if (utils::strmatch(update->integrate_style, "^respa")) {
-    nlevels_respa = (dynamic_cast<Respa *>( update->integrate))->nlevels;
+    nlevels_respa = (static_cast<Respa *>(update->integrate))->nlevels;
+    if (gjfflag) error->all(FLERR, "Fix langevin gjf and run style respa are not compatible");
+  }
 
-  if (utils::strmatch(update->integrate_style,"^respa") && gjfflag)
+  if (gjfflag) {
-    error->all(FLERR,"Fix langevin gjf and respa are not compatible");
+    gjfa = (1.0 - update->dt / 2.0 / t_period) / (1.0 + update->dt / 2.0 / t_period);
-
+    gjfsib = sqrt(1.0 + update->dt / 2.0 / t_period);
-  if (gjfflag) gjfa = (1.0-update->dt/2.0/t_period)/(1.0+update->dt/2.0/t_period);
+  }
-  if (gjfflag) gjfsib = sqrt(1.0+update->dt/2.0/t_period);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -336,13 +335,11 @@ void FixLangevin::setup(int vflag)
           v[i][0] -= dtfm * f[i][0];
           v[i][1] -= dtfm * f[i][1];
           v[i][2] -= dtfm * f[i][2];
-          if (tbiasflag)
+          if (tbiasflag) temperature->remove_bias(i, v[i]);
-            temperature->remove_bias(i,v[i]);
+          v[i][0] /= gjfa * gjfsib * gjfsib;
-          v[i][0] /= gjfa*gjfsib*gjfsib;
+          v[i][1] /= gjfa * gjfsib * gjfsib;
-          v[i][1] /= gjfa*gjfsib*gjfsib;
+          v[i][2] /= gjfa * gjfsib * gjfsib;
-          v[i][2] /= gjfa*gjfsib*gjfsib;
+          if (tbiasflag) temperature->restore_bias(i, v[i]);
-          if (tbiasflag)
-            temperature->restore_bias(i,v[i]);
         }
 
     } else {
@@ -352,22 +349,21 @@ void FixLangevin::setup(int vflag)
           v[i][0] -= dtfm * f[i][0];
           v[i][1] -= dtfm * f[i][1];
           v[i][2] -= dtfm * f[i][2];
-          if (tbiasflag)
+          if (tbiasflag) temperature->remove_bias(i, v[i]);
-            temperature->remove_bias(i,v[i]);
+          v[i][0] /= gjfa * gjfsib * gjfsib;
-          v[i][0] /= gjfa*gjfsib*gjfsib;
+          v[i][1] /= gjfa * gjfsib * gjfsib;
-          v[i][1] /= gjfa*gjfsib*gjfsib;
+          v[i][2] /= gjfa * gjfsib * gjfsib;
-          v[i][2] /= gjfa*gjfsib*gjfsib;
+          if (tbiasflag) temperature->restore_bias(i, v[i]);
-          if (tbiasflag)
-            temperature->restore_bias(i,v[i]);
         }
     }
   }
-  if (utils::strmatch(update->integrate_style,"^verlet"))
+  if (utils::strmatch(update->integrate_style, "^verlet"))
     post_force(vflag);
   else {
-    (dynamic_cast<Respa *>( update->integrate))->copy_flevel_f(nlevels_respa-1);
+    auto respa = static_cast<Respa *>(update->integrate);
-    post_force_respa(vflag,nlevels_respa-1,0);
+    respa->copy_flevel_f(nlevels_respa - 1);
-    (dynamic_cast<Respa *>( update->integrate))->copy_f_flevel(nlevels_respa-1);
+    post_force_respa(vflag, nlevels_respa - 1, 0);
+    respa->copy_f_flevel(nlevels_respa - 1);
   }
   if (gjfflag) {
     double dtfm;
@@ -390,7 +386,7 @@ void FixLangevin::setup(int vflag)
           lv[i][1] = v[i][1];
           lv[i][2] = v[i][2];
         }
-//
+      //
     } else {
       for (int i = 0; i < nlocal; i++)
         if (mask[i] & groupbit) {
@@ -427,6 +423,7 @@ void FixLangevin::initial_integrate(int /* vflag */)
 }
 
 /* ---------------------------------------------------------------------- */
+// clang-format off
 
 void FixLangevin::post_force(int /*vflag*/)
 {
@@ -575,8 +572,7 @@ void FixLangevin::post_force_respa(int vflag, int ilevel, int /*iloop*/)
    modify forces using one of the many Langevin styles
 ------------------------------------------------------------------------- */
 
-template < int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY,
+template<int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY, int Tp_BIAS, int Tp_RMASS, int Tp_ZERO>
-           int Tp_BIAS, int Tp_RMASS, int Tp_ZERO >
 void FixLangevin::post_force_templated()
 {
   double gamma1,gamma2;
@@ -988,6 +984,7 @@ void FixLangevin::end_of_step()
   energy += energy_onestep*update->dt;
 }
 
+// clang-format on
 /* ---------------------------------------------------------------------- */
 
 void FixLangevin::reset_target(double t_new)
@@ -1001,20 +998,17 @@ void FixLangevin::reset_dt()
 {
   if (atom->mass) {
     for (int i = 1; i <= atom->ntypes; i++) {
+      gfactor2[i] = sqrt(atom->mass[i]) / force->ftm2v;
       if (gjfflag)
-        gfactor2[i] = sqrt(atom->mass[i]) *
+        gfactor2[i] *= sqrt(2.0 * force->boltz / t_period / update->dt / force->mvv2e);
-          sqrt(2.0*force->boltz/t_period/update->dt/force->mvv2e) /
-          force->ftm2v;
       else
-        gfactor2[i] = sqrt(atom->mass[i]) *
+        gfactor2[i] *= sqrt(24.0 * force->boltz / t_period / update->dt / force->mvv2e);
-          sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) /
+      gfactor2[i] *= 1.0 / sqrt(ratio[i]);
-          force->ftm2v;
-      gfactor2[i] *= 1.0/sqrt(ratio[i]);
     }
   }
   if (gjfflag) {
-    gjfa = (1.0-update->dt/2.0/t_period)/(1.0+update->dt/2.0/t_period);
+    gjfa = (1.0 - update->dt / 2.0 / t_period) / (1.0 + update->dt / 2.0 / t_period);
-    gjfsib = sqrt(1.0+update->dt/2.0/t_period);
+    gjfsib = sqrt(1.0 + update->dt / 2.0 / t_period);
   }
 }
 
@@ -1022,21 +1016,19 @@ void FixLangevin::reset_dt()
 
 int FixLangevin::modify_param(int narg, char **arg)
 {
-  if (strcmp(arg[0],"temp") == 0) {
+  if (strcmp(arg[0], "temp") == 0) {
-    if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
+    if (narg < 2) utils::missing_cmd_args(FLERR, "fix_modify", error);
-    delete [] id_temp;
+    delete[] id_temp;
     id_temp = utils::strdup(arg[1]);
-
+    temperature = modify->get_compute_by_id(id_temp);
-    int icompute = modify->find_compute(id_temp);
+    if (!temperature)
-    if (icompute < 0)
+      error->all(FLERR, "Could not find fix_modify temperature compute ID: {}", id_temp);
-      error->all(FLERR,"Could not find fix_modify temperature ID");
-    temperature = modify->compute[icompute];
 
     if (temperature->tempflag == 0)
-      error->all(FLERR,
+      error->all(FLERR, "Fix_modify temperature compute {} does not compute temperature", id_temp);
-                 "Fix_modify temperature ID does not compute temperature");
     if (temperature->igroup != igroup && comm->me == 0)
-      error->warning(FLERR,"Group for fix_modify temp != fix group");
+      error->warning(FLERR, "Group for fix_modify temp != fix group: {} vs {}",
+                     group->names[igroup], group->names[temperature->igroup]);
     return 2;
   }
   return 0;
@@ -1059,29 +1051,27 @@ double FixLangevin::compute_scalar()
     if (!gjfflag) {
       for (int i = 0; i < nlocal; i++)
         if (mask[i] & groupbit)
-          energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
+          energy_onestep +=
-                            flangevin[i][2]*v[i][2];
+              flangevin[i][0] * v[i][0] + flangevin[i][1] * v[i][1] + flangevin[i][2] * v[i][2];
-      energy = 0.5*energy_onestep*update->dt;
+      energy = 0.5 * energy_onestep * update->dt;
     } else {
       for (int i = 0; i < nlocal; i++)
         if (mask[i] & groupbit) {
-          if (tbiasflag)
+          if (tbiasflag) temperature->remove_bias(i, lv[i]);
-            temperature->remove_bias(i, lv[i]);
+          energy_onestep +=
-          energy_onestep += flangevin[i][0]*lv[i][0] + flangevin[i][1]*lv[i][1] +
+              flangevin[i][0] * lv[i][0] + flangevin[i][1] * lv[i][1] + flangevin[i][2] * lv[i][2];
-                            flangevin[i][2]*lv[i][2];
+          if (tbiasflag) temperature->restore_bias(i, lv[i]);
-          if (tbiasflag)
-            temperature->restore_bias(i, lv[i]);
         }
-      energy = -0.5*energy_onestep*update->dt;
+      energy = -0.5 * energy_onestep * update->dt;
     }
   }
 
   // convert midstep energy back to previous fullstep energy
 
-  double energy_me = energy - 0.5*energy_onestep*update->dt;
+  double energy_me = energy - 0.5 * energy_onestep * update->dt;
 
   double energy_all;
-  MPI_Allreduce(&energy_me,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
+  MPI_Allreduce(&energy_me, &energy_all, 1, MPI_DOUBLE, MPI_SUM, world);
   return -energy_all;
 }
 
@@ -1092,9 +1082,7 @@ double FixLangevin::compute_scalar()
 void *FixLangevin::extract(const char *str, int &dim)
 {
   dim = 0;
-  if (strcmp(str,"t_target") == 0) {
+  if (strcmp(str, "t_target") == 0) { return &t_target; }
-    return &t_target;
-  }
   return nullptr;
 }
 
@@ -1105,9 +1093,9 @@ void *FixLangevin::extract(const char *str, int &dim)
 double FixLangevin::memory_usage()
 {
   double bytes = 0.0;
-  if (gjfflag) bytes += (double)atom->nmax*6 * sizeof(double);
+  if (gjfflag) bytes += (double) atom->nmax * 6 * sizeof(double);
-  if (tallyflag || osflag) bytes += (double)atom->nmax*3 * sizeof(double);
+  if (tallyflag || osflag) bytes += (double) atom->nmax * 3 * sizeof(double);
-  if (tforce) bytes += (double)atom->nmax * sizeof(double);
+  if (tforce) bytes += (double) atom->nmax * sizeof(double);
   return bytes;
 }
 
@@ -1117,8 +1105,8 @@ double FixLangevin::memory_usage()
 
 void FixLangevin::grow_arrays(int nmax)
 {
-  memory->grow(franprev,nmax,3,"fix_langevin:franprev");
+  memory->grow(franprev, nmax, 3, "fix_langevin:franprev");
-  memory->grow(lv,nmax,3,"fix_langevin:lv");
+  memory->grow(lv, nmax, 3, "fix_langevin:lv");
 }
 
 /* ----------------------------------------------------------------------

src/fix_press_berendsen.cpp

@@ -21,6 +21,7 @@
 #include "error.h"
 #include "fix_deform.h"
 #include "force.h"
+#include "group.h"
 #include "kspace.h"
 #include "modify.h"
 #include "update.h"
@@ -240,14 +241,14 @@ FixPressBerendsen::FixPressBerendsen(LAMMPS *lmp, int narg, char **arg) :
 
 FixPressBerendsen::~FixPressBerendsen()
 {
-  delete [] rfix;
+  delete[] rfix;
 
   // delete temperature and pressure if fix created them
 
   if (tflag) modify->delete_compute(id_temp);
   if (pflag) modify->delete_compute(id_press);
-  delete [] id_temp;
+  delete[] id_temp;
-  delete [] id_press;
+  delete[] id_press;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -268,9 +269,9 @@ void FixPressBerendsen::init()
 
   // insure no conflict with fix deform
 
-  for (int i = 0; i < modify->nfix; i++)
+  for (const auto &ifix : modify->get_fix_list())
-    if (strcmp(modify->fix[i]->style,"deform") == 0) {
+    if (strcmp(ifix->style, "^deform") == 0) {
-      int *dimflag = (dynamic_cast<FixDeform *>( modify->fix[i]))->dimflag;
+      int *dimflag = static_cast<FixDeform *>(ifix)->dimflag;
       if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) ||
           (p_flag[2] && dimflag[2]))
         error->all(FLERR,"Cannot use fix press/berendsen and "
@@ -279,18 +280,16 @@ void FixPressBerendsen::init()
 
   // set temperature and pressure ptrs
 
-  int icompute = modify->find_compute(id_temp);
+  temperature = modify->get_compute_by_id(id_temp);
-  if (icompute < 0)
+  if (!temperature)
-    error->all(FLERR,"Temperature ID for fix press/berendsen does not exist");
+    error->all(FLERR, "Temperature compute ID {} for fix press/berendsen does not exist", id_temp);
-  temperature = modify->compute[icompute];
 
   if (temperature->tempbias) which = BIAS;
   else which = NOBIAS;
 
-  icompute = modify->find_compute(id_press);
+  pressure = modify->get_compute_by_id(id_press);
-  if (icompute < 0)
+  if (!pressure)
-    error->all(FLERR,"Pressure ID for fix press/berendsen does not exist");
+    error->all(FLERR, "Pressure compute ID {} for fix press/berendsen does not exist", id_press);
-  pressure = modify->compute[icompute];
 
   // Kspace setting
 
@@ -300,13 +299,13 @@ void FixPressBerendsen::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;
 
   for (int i = 0; i < modify->nfix; i++)
     if (modify->fix[i]->rigid_flag) nrigid++;
-  if (nrigid) {
+  if (nrigid > 0) {
     rfix = new int[nrigid];
     nrigid = 0;
     for (int i = 0; i < modify->nfix; i++)
@@ -460,24 +459,25 @@ int FixPressBerendsen::modify_param(int narg, char **arg)
       modify->delete_compute(id_temp);
       tflag = 0;
     }
-    delete [] id_temp;
+    delete[] id_temp;
     id_temp = utils::strdup(arg[1]);
 
-    int icompute = modify->find_compute(arg[1]);
+    temperature = modify->get_compute_by_id(arg[1]);
-    if (icompute < 0) error->all(FLERR,"Could not find fix_modify temperature ID");
+    if (!temperature)
-    temperature = modify->compute[icompute];
+      error->all(FLERR,"Could not find fix_modify temperature compute ID: ", arg[1]);
 
     if (temperature->tempflag == 0)
-      error->all(FLERR,"Fix_modify temperature ID does not compute temperature");
+      error->all(FLERR,"Fix_modify temperature compute {} does not compute temperature", arg[1]);
     if (temperature->igroup != 0 && comm->me == 0)
-      error->warning(FLERR,"Temperature for NPT is not for group all");
+      error->warning(FLERR,"Temperature compute {} for fix {} is not for group all: {}",
+                     arg[1], style, group->names[temperature->igroup]);
 
     // reset id_temp of pressure to new temperature ID
 
-    icompute = modify->find_compute(id_press);
+    auto icompute = modify->get_compute_by_id(id_press);
-    if (icompute < 0)
+    if (!icompute)
-      error->all(FLERR,"Pressure ID for fix press/berendsen does not exist");
+      error->all(FLERR,"Pressure compute ID {} for fix {} does not exist", id_press, style);
-    modify->compute[icompute]->reset_extra_compute_fix(id_temp);
+    icompute->reset_extra_compute_fix(id_temp);
 
     return 2;
 
@@ -487,15 +487,13 @@ int FixPressBerendsen::modify_param(int narg, char **arg)
       modify->delete_compute(id_press);
       pflag = 0;
     }
-    delete [] id_press;
+    delete[] id_press;
     id_press = utils::strdup(arg[1]);
 
-    int icompute = modify->find_compute(arg[1]);
+    pressure = modify->get_compute_by_id(arg[1]);
-    if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID");
+    if (pressure) error->all(FLERR,"Could not find fix_modify pressure compute ID: {}", arg[1]);
-    pressure = modify->compute[icompute];
-
     if (pressure->pressflag == 0)
-      error->all(FLERR,"Fix_modify pressure ID does not compute pressure");
+      error->all(FLERR,"Fix_modify pressure compute {} does not compute pressure", arg[1]);
     return 2;
   }
   return 0;

src/fix_restrain.cpp

@@ -19,17 +19,18 @@
 
 #include "fix_restrain.h"
 
-#include <cmath>
-#include <cstring>
 #include "atom.h"
-#include "force.h"
-#include "update.h"
-#include "domain.h"
 #include "comm.h"
-#include "respa.h"
+#include "domain.h"
+#include "error.h"
+#include "force.h"
 #include "math_const.h"
 #include "memory.h"
-#include "error.h"
+#include "respa.h"
+#include "update.h"
+
+#include <cmath>
+#include <cstring>
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
@@ -271,14 +272,12 @@ void FixRestrain::restrain_bond(int m)
   if (newton_bond) {
     if (i2 == -1 || i2 >= nlocal) return;
     if (i1 == -1)
-      error->one(FLERR,"Restrain atoms {} {} missing on "
+      error->one(FLERR,"Restrain atoms {} {} missing on proc {} at step {}", ids[m][0],ids[m][1],
-                                   "proc {} at step {}", ids[m][0],ids[m][1],
                                    comm->me,update->ntimestep);
   } else {
     if ((i1 == -1 || i1 >= nlocal) && (i2 == -1 || i2 >= nlocal)) return;
     if (i1 == -1 || i2 == -1)
-      error->one(FLERR,"Restrain atoms {} {} missing on "
+      error->one(FLERR,"Restrain atoms {} {} missing on proc {} at step {}", ids[m][0],ids[m][1],
-                                   "proc {} at step {}", ids[m][0],ids[m][1],
                                    comm->me,update->ntimestep);
   }
 

src/fix_temp_berendsen.cpp

@@ -40,7 +40,8 @@ FixTempBerendsen::FixTempBerendsen(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg),
   tstr(nullptr), id_temp(nullptr), tflag(0)
 {
-  if (narg != 6) error->all(FLERR,"Illegal fix temp/berendsen command");
+  if (narg != 6)
+    error->all(FLERR,"Illegal fix {} command: expected 6 arguments but found {}", style, narg);
 
   // Berendsen thermostat should be applied every step
 
@@ -68,7 +69,7 @@ FixTempBerendsen::FixTempBerendsen(LAMMPS *lmp, int narg, char **arg) :
   // error checks
 
   if (t_period <= 0.0)
-    error->all(FLERR,"Fix temp/berendsen period must be > 0.0");
+    error->all(FLERR,"Fix temp/berendsen Tdamp period must be > 0.0");
 
   // create a new compute temp style
   // id = fix-ID + temp, compute group = fix group
@@ -84,12 +85,12 @@ FixTempBerendsen::FixTempBerendsen(LAMMPS *lmp, int narg, char **arg) :
 
 FixTempBerendsen::~FixTempBerendsen()
 {
-  delete [] tstr;
+  delete[] tstr;
 
   // delete temperature if fix created it
 
   if (tflag) modify->delete_compute(id_temp);
-  delete [] id_temp;
+  delete[] id_temp;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -110,18 +111,17 @@ void FixTempBerendsen::init()
   if (tstr) {
     tvar = input->variable->find(tstr);
     if (tvar < 0)
-      error->all(FLERR,"Variable name for fix temp/berendsen does not exist");
+      error->all(FLERR,"Variable name {} for fix temp/berendsen does not exist", tstr);
     if (input->variable->equalstyle(tvar)) tstyle = EQUAL;
-    else error->all(FLERR,"Variable for fix temp/berendsen is invalid style");
+    else error->all(FLERR,"Variable {} for fix temp/berendsen is invalid style", tstr);
   }
 
-  int icompute = modify->find_compute(id_temp);
+  temperature = modify->get_compute_by_id(id_temp);
-  if (icompute < 0)
+  if (!temperature)
-    error->all(FLERR,"Temperature ID for fix temp/berendsen does not exist");
+    error->all(FLERR,"Temperature compute ID {} for fix {} does not exist", id_temp, style);
-  temperature = modify->compute[icompute];
 
   if (modify->check_rigid_group_overlap(groupbit))
-    error->warning(FLERR,"Cannot thermostat atoms in rigid bodies");
+    error->warning(FLERR,"Cannot thermostat atoms in rigid bodies with fix {}", style);
 
   if (temperature->tempbias) which = BIAS;
   else which = NOBIAS;
@@ -139,8 +139,7 @@ void FixTempBerendsen::end_of_step()
   if (tdof < 1) return;
 
   if (t_current == 0.0)
-    error->all(FLERR,
+    error->all(FLERR, "Computed current temperature for fix temp/berendsen must not be 0.0");
-               "Computed temperature for fix temp/berendsen cannot be 0.0");
 
   double delta = update->ntimestep - update->beginstep;
   if (delta != 0.0) delta /= update->endstep - update->beginstep;
@@ -154,8 +153,8 @@ void FixTempBerendsen::end_of_step()
     modify->clearstep_compute();
     t_target = input->variable->compute_equal(tvar);
     if (t_target < 0.0)
-      error->one(FLERR,
+      error->one(FLERR, "Fix temp/berendsen variable {} returned negative temperature",
-                 "Fix temp/berendsen variable returned negative temperature");
+                 input->variable->names[tvar]);
     modify->addstep_compute(update->ntimestep + nevery);
   }
 
@@ -198,24 +197,23 @@ void FixTempBerendsen::end_of_step()
 int FixTempBerendsen::modify_param(int narg, char **arg)
 {
   if (strcmp(arg[0],"temp") == 0) {
-    if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
+    if (narg < 2) utils::missing_cmd_args(FLERR, "fix_modify", error);
     if (tflag) {
       modify->delete_compute(id_temp);
       tflag = 0;
     }
-    delete [] id_temp;
+    delete[] id_temp;
     id_temp = utils::strdup(arg[1]);
 
-    int icompute = modify->find_compute(id_temp);
+    temperature = modify->get_compute_by_id(id_temp);
-    if (icompute < 0)
+    if (!temperature)
-      error->all(FLERR,"Could not find fix_modify temperature ID");
+      error->all(FLERR,"Could not find fix_modify temperature compute {}", id_temp);
-    temperature = modify->compute[icompute];
 
     if (temperature->tempflag == 0)
-      error->all(FLERR,
+      error->all(FLERR, "Fix_modify temperature compute {} does not compute temperature", id_temp);
-                 "Fix_modify temperature ID does not compute temperature");
     if (temperature->igroup != igroup && comm->me == 0)
-      error->warning(FLERR,"Group for fix_modify temp != fix group");
+      error->warning(FLERR, "Group for fix_modify temp != fix group: {} vs {}",
+                     group->names[igroup], group->names[temperature->igroup]);
     return 2;
   }
   return 0;

src/fix_temp_rescale.cpp

@@ -43,7 +43,7 @@ FixTempRescale::FixTempRescale(LAMMPS *lmp, int narg, char **arg) :
   if (narg < 8) utils::missing_cmd_args(FLERR, "fix temp/rescale", error);
 
   nevery = utils::inumeric(FLERR,arg[3],false,lmp);
-  if (nevery <= 0) error->all(FLERR, "Invalid fix temp/rescale argument: {}", nevery);
+  if (nevery <= 0) error->all(FLERR, "Invalid fix temp/rescale every argument: {}", nevery);
 
   restart_global = 1;
   scalar_flag = 1;
@@ -66,6 +66,10 @@ FixTempRescale::FixTempRescale(LAMMPS *lmp, int narg, char **arg) :
   t_window = utils::numeric(FLERR,arg[6],false,lmp);
   fraction = utils::numeric(FLERR,arg[7],false,lmp);
 
+  if (t_stop < 0) error->all(FLERR, "Invalid fix temp/rescale Tstop argument: {}", t_stop);
+  if (t_window < 0) error->all(FLERR, "Invalid fix temp/rescale window argument: {}", t_window);
+  if (fraction <= 0) error->all(FLERR, "Invalid fix temp/rescale fraction argument: {}", fraction);
+
   // create a new compute temp
   // id = fix-ID + temp, compute group = fix group
 
@@ -80,12 +84,12 @@ FixTempRescale::FixTempRescale(LAMMPS *lmp, int narg, char **arg) :
 
 FixTempRescale::~FixTempRescale()
 {
-  delete [] tstr;
+  delete[] tstr;
 
   // delete temperature if fix created it
 
   if (tflag) modify->delete_compute(id_temp);
-  delete [] id_temp;
+  delete[] id_temp;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -106,15 +110,14 @@ void FixTempRescale::init()
   if (tstr) {
     tvar = input->variable->find(tstr);
     if (tvar < 0)
-      error->all(FLERR,"Variable {} for fix temp/rescale does not exist", tstr);
+      error->all(FLERR,"Variable name {} for fix temp/rescale does not exist", tstr);
     if (input->variable->equalstyle(tvar)) tstyle = EQUAL;
     else error->all(FLERR,"Variable {} for fix temp/rescale is invalid style", tstr);
   }
 
-  int icompute = modify->find_compute(id_temp);
+  temperature = modify->get_compute_by_id(id_temp);
-  if (icompute < 0)
+  if (!temperature)
-    error->all(FLERR,"Temperature ID for fix temp/rescale does not exist");
+    error->all(FLERR,"Temperature ID {} for fix temp/rescale does not exist", id_temp);
-  temperature = modify->compute[icompute];
 
   if (temperature->tempbias) which = BIAS;
   else which = NOBIAS;
@@ -147,8 +150,7 @@ void FixTempRescale::end_of_step()
     modify->clearstep_compute();
     t_target = input->variable->compute_equal(tvar);
     if (t_target < 0.0)
-      error->one(FLERR,
+      error->one(FLERR, "Fix temp/rescale variable returned negative temperature");
-                 "Fix temp/rescale variable returned negative temperature");
     modify->addstep_compute(update->ntimestep + nevery);
   }
 
@@ -200,19 +202,18 @@ int FixTempRescale::modify_param(int narg, char **arg)
       modify->delete_compute(id_temp);
       tflag = 0;
     }
-    delete [] id_temp;
+    delete[] id_temp;
     id_temp = utils::strdup(arg[1]);
 
-    int icompute = modify->find_compute(id_temp);
+    temperature = modify->get_compute_by_id(id_temp);
-    if (icompute < 0)
+    if (!temperature)
-      error->all(FLERR,"Could not find fix_modify temperature ID {}", arg[1]);
+      error->all(FLERR,"Could not find fix_modify temperature compute {}", id_temp);
-    temperature = modify->compute[icompute];
 
     if (temperature->tempflag == 0)
-      error->all(FLERR,
+      error->all(FLERR, "Fix_modify temperature compute {} does not compute temperature", id_temp);
-                 "Fix_modify temperature ID does not compute temperature");
     if (temperature->igroup != igroup && comm->me == 0)
-      error->warning(FLERR,"Group for fix_modify temp != fix group");
+      error->warning(FLERR, "Group for fix_modify temp != fix group: {} vs {}",
+                     group->names[igroup], group->names[temperature->igroup]);
     return 2;
   }
   return 0;

src/fix_vector.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/, Sandia National Laboratories
@@ -26,41 +25,34 @@
 using namespace LAMMPS_NS;
 using namespace FixConst;
 
-enum{ONE,RUNNING,WINDOW};
+enum { ONE, RUNNING, WINDOW };
-enum{SCALAR,VECTOR};
+enum { SCALAR, VECTOR };
-
 
 /* ---------------------------------------------------------------------- */
 
 FixVector::FixVector(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg),
+    Fix(lmp, narg, arg), vector(nullptr), array(nullptr)
-  nvalues(0), which(nullptr), argindex(nullptr), value2index(nullptr), ids(nullptr), vector(nullptr), array(nullptr)
 {
-  if (narg < 5) error->all(FLERR,"Illegal fix vector command");
+  if (narg < 5) utils::missing_cmd_args(FLERR, "fix vector", error);
 
-  nevery = utils::inumeric(FLERR,arg[3],false,lmp);
+  nevery = utils::inumeric(FLERR, arg[3], false, lmp);
-  if (nevery <= 0) error->all(FLERR,"Illegal fix vector command");
+  if (nevery <= 0) error->all(FLERR, "Invalid fix vector every argument: {}", nevery);
 
-  nvalues = narg-4;
+  values.clear();
-  which = new int[nvalues];
-  argindex = new int[nvalues];
-  value2index = new int[nvalues];
-  ids = new char*[nvalues];
-
-  nvalues = 0;
   for (int iarg = 4; iarg < narg; iarg++) {
     ArgInfo argi(arg[iarg]);
 
-    which[nvalues] = argi.get_type();
+    value_t val;
-    argindex[nvalues] = argi.get_index1();
+    val.which = argi.get_type();
-    ids[nvalues] = argi.copy_name();
+    val.argindex = argi.get_index1();
+    val.id = argi.get_name();
+    val.val.c = nullptr;
 
-    if ((argi.get_type() == ArgInfo::UNKNOWN)
+    if ((argi.get_type() == ArgInfo::UNKNOWN) || (argi.get_type() == ArgInfo::NONE) ||
-        || (argi.get_type() == ArgInfo::NONE)
+        (argi.get_dim() > 1))
-        || (argi.get_dim() > 1))
+      error->all(FLERR, "Invalid fix vector argument: {}", arg[iarg]);
-      error->all(FLERR,"Illegal fix vector command");
 
-    nvalues++;
+    values.push_back(val);
   }
 
   // setup and error check
@@ -68,59 +60,70 @@ FixVector::FixVector(LAMMPS *lmp, int narg, char **arg) :
   // this fix produces either a global vector or array
   // intensive/extensive flags set by compute,fix,variable that produces value
 
-  int value,finalvalue;
+  int value, finalvalue;
-  for (int i = 0; i < nvalues; i++) {
+  bool first = true;
-    if (which[i] == ArgInfo::COMPUTE) {
+  for (auto &val : values) {
-      auto icompute = modify->get_compute_by_id(ids[i]);
+    if (val.which == ArgInfo::COMPUTE) {
-      if (!icompute) error->all(FLERR,"Compute ID {} for fix vector does not exist",ids[i]);
+      auto icompute = modify->get_compute_by_id(val.id);
-      if (argindex[i] == 0 && icompute->scalar_flag == 0)
+      if (!icompute) error->all(FLERR, "Compute ID {} for fix vector does not exist", val.id);
-        error->all(FLERR,"Fix vector compute {} does not calculate a scalar",ids[i]);
+      if (val.argindex == 0 && icompute->scalar_flag == 0)
-      if (argindex[i] && icompute->vector_flag == 0)
+        error->all(FLERR, "Fix vector compute {} does not calculate a scalar", val.id);
-        error->all(FLERR,"Fix vector compute {} does not calculate a vector",ids[i]);
+      if (val.argindex && icompute->vector_flag == 0)
-      if (argindex[i] && argindex[i] > icompute->size_vector)
+        error->all(FLERR, "Fix vector compute {} does not calculate a vector", val.id);
-        error->all(FLERR,"Fix vector compute {} vector is accessed out-of-range",ids[i]);
+      if (val.argindex && (val.argindex > icompute->size_vector))
+        error->all(FLERR, "Fix vector compute {} vector is accessed out-of-range", val.id);
 
-      if (argindex[i] == 0) value = icompute->extscalar;
+      if (val.argindex == 0)
-      else if (icompute->extvector >= 0) value = icompute->extvector;
+        value = icompute->extscalar;
-      else value = icompute->extlist[argindex[i]-1];
+      else if (icompute->extvector >= 0)
+        value = icompute->extvector;
+      else
+        value = icompute->extlist[val.argindex - 1];
+      val.val.c = icompute;
 
-    } else if (which[i] == ArgInfo::FIX) {
+    } else if (val.which == ArgInfo::FIX) {
-      auto ifix = modify->get_fix_by_id(ids[i]);
+      auto ifix = modify->get_fix_by_id(val.id);
-      if (!ifix) error->all(FLERR,"Fix ID {} for fix vector does not exist",ids[i]);
+      if (!ifix) error->all(FLERR, "Fix ID {} for fix vector does not exist", val.id);
-      if (argindex[i] == 0 && ifix->scalar_flag == 0)
+      if (val.argindex == 0 && ifix->scalar_flag == 0)
-        error->all(FLERR,"Fix vector fix {} does not calculate a scalar",ids[i]);
+        error->all(FLERR, "Fix vector fix {} does not calculate a scalar", val.id);
-      if (argindex[i] && ifix->vector_flag == 0)
+      if (val.argindex && ifix->vector_flag == 0)
-        error->all(FLERR,"Fix vector fix {} does not calculate a vector",ids[i]);
+        error->all(FLERR, "Fix vector fix {} does not calculate a vector", val.id);
-      if (argindex[i] && argindex[i] > ifix->size_vector)
+      if (val.argindex && val.argindex > ifix->size_vector)
-        error->all(FLERR,"Fix vector fix {} vector is accessed out-of-range",ids[i]);
+        error->all(FLERR, "Fix vector fix {} vector is accessed out-of-range", val.id);
       if (nevery % ifix->global_freq)
-        error->all(FLERR,"Fix for fix {} vector not computed at compatible time",ids[i]);
+        error->all(FLERR, "Fix for fix {} vector not computed at compatible time", val.id);
 
-      if (argindex[i] == 0) value = ifix->extvector;
+      if (val.argindex == 0)
-      else value = ifix->extarray;
+        value = ifix->extvector;
+      else
+        value = ifix->extarray;
+      val.val.f = ifix;
 
-    } else if (which[i] == ArgInfo::VARIABLE) {
+    } else if (val.which == ArgInfo::VARIABLE) {
-      int ivariable = input->variable->find(ids[i]);
+      int ivariable = input->variable->find(val.id.c_str());
       if (ivariable < 0)
-        error->all(FLERR,"Variable name {} for fix vector does not exist",ids[i]);
+        error->all(FLERR, "Variable name {} for fix vector does not exist", val.id);
-      if (argindex[i] == 0 && input->variable->equalstyle(ivariable) == 0)
+      if (val.argindex == 0 && input->variable->equalstyle(ivariable) == 0)
-        error->all(FLERR,"Fix vector variable {} is not equal-style variable",ids[i]);
+        error->all(FLERR, "Fix vector variable {} is not equal-style variable", val.id);
-      if (argindex[i] && input->variable->vectorstyle(ivariable) == 0)
+      if (val.argindex && input->variable->vectorstyle(ivariable) == 0)
-        error->all(FLERR,"Fix vector variable {} is not vector-style variable",ids[i]);
+        error->all(FLERR, "Fix vector variable {} is not vector-style variable", val.id);
       value = 0;
+      val.val.v = ivariable;
     }
 
-    if (i == 0) finalvalue = value;
+    if (first) {
-    else if (value != finalvalue)
+      finalvalue = value;
-      error->all(FLERR,"Fix vector cannot set output array intensive/extensive from these inputs");
+      first = false;
+    } else if (value != finalvalue)
+      error->all(FLERR, "Fix vector cannot set output array intensive/extensive from these inputs");
   }
 
-  if (nvalues == 1) {
+  if (values.size() == 1) {
     vector_flag = 1;
     extvector = finalvalue;
   } else {
     array_flag = 1;
-    size_array_cols = nvalues;
+    size_array_cols = values.size();
     extarray = finalvalue;
   }
 
@@ -132,15 +135,17 @@ FixVector::FixVector(LAMMPS *lmp, int narg, char **arg) :
   vector = nullptr;
   array = nullptr;
   ncount = ncountmax = 0;
-  if (nvalues == 1) size_vector = 0;
+  if (values.size() == 1)
-  else size_array_rows = 0;
+    size_vector = 0;
+  else
+    size_array_rows = 0;
 
   // nextstep = next step on which end_of_step does something
   // add nextstep to all computes that store invocation times
   // since don't know a priori which are invoked by this fix
   // once in end_of_step() can set timestep for ones actually invoked
 
-  nextstep = (update->ntimestep/nevery)*nevery;
+  nextstep = (update->ntimestep / nevery) * nevery;
   if (nextstep < update->ntimestep) nextstep += nevery;
   modify->addstep_compute_all(nextstep);
 
@@ -153,12 +158,7 @@ FixVector::FixVector(LAMMPS *lmp, int narg, char **arg) :
 
 FixVector::~FixVector()
 {
-  delete [] which;
+  values.clear();
-  delete [] argindex;
-  delete [] value2index;
-  for (int i = 0; i < nvalues; i++) delete [] ids[i];
-  delete [] ids;
-
   memory->destroy(vector);
   memory->destroy(array);
 }
@@ -178,22 +178,19 @@ void FixVector::init()
 {
   // set current indices for all computes,fixes,variables
 
-  for (int i = 0; i < nvalues; i++) {
+  for (auto &val : values) {
-    if (which[i] == ArgInfo::COMPUTE) {
+    if (val.which == ArgInfo::COMPUTE) {
-      int icompute = modify->find_compute(ids[i]);
+      val.val.c = modify->get_compute_by_id(val.id);
-      if (icompute < 0) error->all(FLERR,"Compute ID {} for fix vector does not exist",id[i]);
+      if (!val.val.c) error->all(FLERR, "Compute ID {} for fix vector does not exist", val.id);
-      value2index[i] = icompute;
 
-    } else if (which[i] == ArgInfo::FIX) {
+    } else if (val.which == ArgInfo::FIX) {
-      int ifix = modify->find_fix(ids[i]);
+      val.val.f = modify->get_fix_by_id(val.id);
-      if (ifix < 0) error->all(FLERR,"Fix ID {} for fix vector does not exist",id[i]);
+      if (!val.val.f) error->all(FLERR, "Fix ID {} for fix vector does not exist", val.id);
-      value2index[i] = ifix;
 
-    } else if (which[i] == ArgInfo::VARIABLE) {
+    } else if (val.which == ArgInfo::VARIABLE) {
-      int ivariable = input->variable->find(ids[i]);
+      int ivariable = input->variable->find(val.id.c_str());
-      if (ivariable < 0)
+      if (ivariable < 0) error->all(FLERR, "Variable name for fix vector does not exist");
-        error->all(FLERR,"Variable name for fix vector does not exist");
+      val.val.v = ivariable;
-      value2index[i] = ivariable;
     }
   }
 
@@ -201,11 +198,13 @@ void FixVector::init()
   // use endstep to allow for subsequent runs with "pre no"
   // nsize = # of entries from initialstep to finalstep
 
-  bigint finalstep = update->endstep/nevery * nevery;
+  bigint finalstep = update->endstep / nevery * nevery;
   if (finalstep > update->endstep) finalstep -= nevery;
-  ncountmax = (finalstep-initialstep)/nevery + 1;
+  ncountmax = (finalstep - initialstep) / nevery + 1;
-  if (nvalues == 1) memory->grow(vector,ncountmax,"vector:vector");
+  if (values.size() == 1)
-  else memory->grow(array,ncountmax,nvalues,"vector:array");
+    memory->grow(vector, ncountmax, "vector:vector");
+  else
+    memory->grow(array, ncountmax, values.size(), "vector:array");
 }
 
 /* ----------------------------------------------------------------------
@@ -224,61 +223,64 @@ void FixVector::end_of_step()
   // skip if not step which requires doing something
 
   if (update->ntimestep != nextstep) return;
-  if (ncount == ncountmax)
+  if (ncount == ncountmax) error->all(FLERR, "Overflow of allocated fix vector storage");
-    error->all(FLERR,"Overflow of allocated fix vector storage");
 
   // accumulate results of computes,fixes,variables to local copy
   // compute/fix/variable may invoke computes so wrap with clear/add
 
   double *result;
-  if (nvalues == 1) result = &vector[ncount];
+  if (values.size() == 1)
-  else result = array[ncount];
+    result = &vector[ncount];
+  else
+    result = array[ncount];
 
   modify->clearstep_compute();
 
-  for (int i = 0; i < nvalues; i++) {
+  int i = 0;
-    int m = value2index[i];
+  for (auto &val : values) {
 
     // invoke compute if not previously invoked
 
-    if (which[i] == ArgInfo::COMPUTE) {
+    if (val.which == ArgInfo::COMPUTE) {
-      auto compute = modify->get_compute_by_index(m);
 
-      if (argindex[i] == 0) {
+      if (val.argindex == 0) {
-        if (!(compute->invoked_flag & Compute::INVOKED_SCALAR)) {
+        if (!(val.val.c->invoked_flag & Compute::INVOKED_SCALAR)) {
-          compute->compute_scalar();
+          val.val.c->compute_scalar();
-          compute->invoked_flag |= Compute::INVOKED_SCALAR;
+          val.val.c->invoked_flag |= Compute::INVOKED_SCALAR;
         }
-        result[i] = compute->scalar;
+        result[i] = val.val.c->scalar;
       } else {
-        if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) {
+        if (!(val.val.c->invoked_flag & Compute::INVOKED_VECTOR)) {
-          compute->compute_vector();
+          val.val.c->compute_vector();
-          compute->invoked_flag |= Compute::INVOKED_VECTOR;
+          val.val.c->invoked_flag |= Compute::INVOKED_VECTOR;
         }
-        result[i] = compute->vector[argindex[i]-1];
+        result[i] = val.val.c->vector[val.argindex - 1];
       }
 
-    // access fix fields, guaranteed to be ready
+      // access fix fields, guaranteed to be ready
 
-    } else if (which[i] == ArgInfo::FIX) {
+    } else if (val.which == ArgInfo::FIX) {
-      if (argindex[i] == 0)
+      if (val.argindex == 0)
-        result[i] = modify->get_fix_by_index(m)->compute_scalar();
+        result[i] = val.val.f->compute_scalar();
       else
-        result[i] = modify->get_fix_by_index(m)->compute_vector(argindex[i]-1);
+        result[i] = val.val.f->compute_vector(val.argindex - 1);
 
-    // evaluate equal-style or vector-style variable
+      // evaluate equal-style or vector-style variable
 
-    } else if (which[i] == ArgInfo::VARIABLE) {
+    } else if (val.which == ArgInfo::VARIABLE) {
-      if (argindex[i] == 0)
+      if (val.argindex == 0)
-        result[i] = input->variable->compute_equal(m);
+        result[i] = input->variable->compute_equal(val.val.v);
       else {
         double *varvec;
-        int nvec = input->variable->compute_vector(m,&varvec);
+        int nvec = input->variable->compute_vector(val.val.v, &varvec);
-        int index = argindex[i];
+        int index = val.argindex;
-        if (nvec < index) result[i] = 0.0;
+        if (nvec < index)
-        else result[i] = varvec[index-1];
+          result[i] = 0.0;
+        else
+          result[i] = varvec[index - 1];
       }
     }
+    ++i;
   }
 
   // trigger computes on next needed step
@@ -289,8 +291,10 @@ void FixVector::end_of_step()
   // update size of vector or array
 
   ncount++;
-  if (nvalues == 1) size_vector++;
+  if (values.size() == 1)
-  else size_array_rows++;
+    size_vector++;
+  else
+    size_array_rows++;
 }
 
 /* ----------------------------------------------------------------------

src/fix_vector.h

@@ -36,9 +36,17 @@ class FixVector : public Fix {
   double compute_array(int, int) override;
 
  private:
-  int nvalues;
+  struct value_t {
-  int *which, *argindex, *value2index;
+    int which;
-  char **ids;
+    int argindex;
+    std::string id;
+    union {
+      class Compute *c;
+      class Fix *f;
+      int v;
+    } val;
+  };
+  std::vector<value_t> values;
 
   bigint nextstep, initialstep;
 
@@ -47,8 +55,6 @@ class FixVector : public Fix {
   double *vector;
   double **array;
 };
-
 }    // namespace LAMMPS_NS
-
 #endif
 #endif

src/group.cpp

@@ -101,22 +101,21 @@ void Group::assign(int narg, char **arg)
   // clear mask of each atom assigned to this group
 
   if (strcmp(arg[1],"delete") == 0) {
-    if (narg != 2) error->all(FLERR,"Illegal group command");
+    if (narg != 2) error->all(FLERR,"Illegal group command: too many arguments");
     int igroup = find(arg[0]);
-    if (igroup == -1) error->all(FLERR,"Could not find group delete group ID");
+    if (igroup == -1) error->all(FLERR,"Could not find group delete group ID {}",arg[0]);
     if (igroup == 0) error->all(FLERR,"Cannot delete group all");
-    for (const auto &fix : modify->get_fix_list())
+    for (const auto &i : modify->get_fix_list())
-      if (fix->igroup == igroup)
+      if (i->igroup == igroup)
-        error->all(FLERR,"Cannot delete group currently used by a fix");
+        error->all(FLERR,"Cannot delete group {} currently used by fix ID {}",arg[0],i->id);
-    for (i = 0; i < modify->ncompute; i++)
+    for (const auto &i : modify->get_compute_list())
-      if (modify->compute[i]->igroup == igroup)
+      if (i->igroup == igroup)
-        error->all(FLERR,"Cannot delete group currently used by a compute");
+        error->all(FLERR,"Cannot delete group {} currently used by compute ID {}",arg[0],i->id);
-    for (i = 0; i < output->ndump; i++)
+    for (const auto &i : output->get_dump_list())
-      if (output->dump[i]->igroup == igroup)
+      if (i->igroup == igroup)
-        error->all(FLERR,"Cannot delete group currently used by a dump");
+        error->all(FLERR,"Cannot delete group {} currently used by dump ID {}",arg[0],i->id);
     if (atom->firstgroupname && strcmp(arg[0],atom->firstgroupname) == 0)
-      error->all(FLERR,
+      error->all(FLERR,"Cannot delete group {} currently used by atom_modify first",arg[0]);
-                 "Cannot delete group currently used by atom_modify first");
 
     int *mask = atom->mask;
     int nlocal = atom->nlocal;

src/info.cpp

@@ -887,42 +887,17 @@ bool Info::is_defined(const char *category, const char *name)
   if ((category == nullptr) || (name == nullptr)) return false;
 
   if (strcmp(category,"compute") == 0) {
-    int ncompute = modify->ncompute;
+    if (modify->get_compute_by_id(name)) return true;
-    Compute **compute = modify->compute;
-    for (int i=0; i < ncompute; ++i) {
-      if (strcmp(compute[i]->id,name) == 0)
-        return true;
-    }
   } else if (strcmp(category,"dump") == 0) {
-    int ndump = output->ndump;
+    if (output->get_dump_by_id(name)) return true;
-    Dump **dump = output->dump;
-    for (int i=0; i < ndump; ++i) {
-      if (strcmp(dump[i]->id,name) == 0)
-        return true;
-    }
   } else if (strcmp(category,"fix") == 0) {
-    for (const auto &fix : modify->get_fix_list()) {
+    if (modify->get_fix_by_id(name)) return true;
-      if (strcmp(fix->id,name) == 0)
-        return true;
-    }
   } else if (strcmp(category,"group") == 0) {
-    int ngroup = group->ngroup;
+    if (group->find(name) >= 0) return true;
-    char **names = group->names;
-    for (int i=0; i < ngroup; ++i) {
-      if (strcmp(names[i],name) == 0)
-        return true;
-    }
   } else if (strcmp(category,"region") == 0) {
-    for (auto &reg : domain->get_region_list())
+    if (domain->get_region_by_id(name)) return true;
-      if (strcmp(reg->id,name) == 0) return true;
   } else if (strcmp(category,"variable") == 0) {
-    int nvar = input->variable->nvar;
+    if (input->variable->find(name) >= 0) return true;
-    char **names = input->variable->names;
-
-    for (int i=0; i < nvar; ++i) {
-      if (strcmp(names[i],name) == 0)
-        return true;
-    }
   } else error->all(FLERR,"Unknown category for info is_defined(): {}", category);
 
   return false;

src/library.cpp

@@ -4761,43 +4761,19 @@ int lammps_has_id(void *handle, const char *category, const char *name) {
   auto lmp = (LAMMPS *) handle;
 
   if (strcmp(category,"compute") == 0) {
-    int ncompute = lmp->modify->ncompute;
+    if (lmp->modify->get_compute_by_id(name)) return 1;
-    Compute **compute = lmp->modify->compute;
-    for (int i=0; i < ncompute; ++i) {
-      if (strcmp(name,compute[i]->id) == 0) return 1;
-    }
   } else if (strcmp(category,"dump") == 0) {
-    int ndump = lmp->output->ndump;
+    if (lmp->output->get_dump_by_id(name)) return 1;
-    Dump **dump = lmp->output->dump;
-    for (int i=0; i < ndump; ++i) {
-      if (strcmp(name,dump[i]->id) == 0) return 1;
-    }
   } else if (strcmp(category,"fix") == 0) {
-    for (auto &ifix : lmp->modify->get_fix_list()) {
+    if (lmp->modify->get_fix_by_id(name)) return 1;
-      if (strcmp(name,ifix->id) == 0) return 1;
-    }
   } else if (strcmp(category,"group") == 0) {
-    int ngroup = lmp->group->ngroup;
+    if (lmp->group->find(name) >= 0) return 1;
-    char **groups = lmp->group->names;
-    for (int i=0; i < ngroup; ++i) {
-      if (strcmp(groups[i],name) == 0) return 1;
-    }
   } else if (strcmp(category,"molecule") == 0) {
-    int nmolecule = lmp->atom->nmolecule;
+    if (lmp->atom->find_molecule(name) >= 0) return 1;
-    Molecule **molecule = lmp->atom->molecules;
-    for (int i=0; i < nmolecule; ++i) {
-      if (strcmp(name,molecule[i]->id) == 0) return 1;
-    }
   } else if (strcmp(category,"region") == 0) {
-    for (auto &reg : lmp->domain->get_region_list()) {
+    if (lmp->domain->get_region_by_id(name)) return 1;
-      if (strcmp(name,reg->id) == 0) return 1;
-    }
   } else if (strcmp(category,"variable") == 0) {
-    int nvariable = lmp->input->variable->nvar;
+    if (lmp->input->variable->find(name) >= 0) return 1;
-    char **varnames = lmp->input->variable->names;
-    for (int i=0; i < nvariable; ++i) {
-      if (strcmp(name,varnames[i]) == 0) return 1;
-    }
   }
   return 0;
 }
@@ -4823,11 +4799,11 @@ categories.
 int lammps_id_count(void *handle, const char *category) {
   auto lmp = (LAMMPS *) handle;
   if (strcmp(category,"compute") == 0) {
-    return lmp->modify->ncompute;
+    return lmp->modify->get_compute_list().size();
   } else if (strcmp(category,"dump") == 0) {
-    return lmp->output->ndump;
+    return lmp->output->get_dump_list().size();
   } else if (strcmp(category,"fix") == 0) {
-    return lmp->modify->nfix;
+    return lmp->modify->get_fix_list().size();
   } else if (strcmp(category,"group") == 0) {
     return lmp->group->ngroup;
   } else if (strcmp(category,"molecule") == 0) {
@@ -4865,6 +4841,7 @@ set to an empty string, otherwise 1.
  */
 int lammps_id_name(void *handle, const char *category, int idx, char *buffer, int buf_size) {
   auto lmp = (LAMMPS *) handle;
+  if (idx < 0) return 0;
 
   if (strcmp(category,"compute") == 0) {
     auto icompute = lmp->modify->get_compute_by_index(idx);
@@ -4873,8 +4850,9 @@ int lammps_id_name(void *handle, const char *category, int idx, char *buffer, in
       return 1;
     }
   } else if (strcmp(category,"dump") == 0) {
-    if ((idx >=0) && (idx < lmp->output->ndump)) {
+    auto idump = lmp->output->get_dump_by_index(idx);
-      strncpy(buffer, lmp->output->dump[idx]->id, buf_size);
+    if (idump) {
+      strncpy(buffer, idump->id, buf_size);
       return 1;
     }
   } else if (strcmp(category,"fix") == 0) {

src/min.cpp

@@ -94,13 +94,13 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp)
 
 Min::~Min()
 {
-  delete [] elist_global;
+  delete[] elist_global;
-  delete [] elist_atom;
+  delete[] elist_atom;
-  delete [] vlist_global;
+  delete[] vlist_global;
-  delete [] vlist_atom;
+  delete[] vlist_atom;
-  delete [] cvlist_atom;
+  delete[] cvlist_atom;
 
-  delete [] fextra;
+  delete[] fextra;
 
   memory->sfree(xextra_atom);
   memory->sfree(fextra_atom);
@@ -128,7 +128,7 @@ void Min::init()
   // can then add vectors to fix_minimize in setup()
 
   nextra_global = 0;
-  delete [] fextra;
+  delete[] fextra;
   fextra = nullptr;
 
   nextra_atom = 0;
@@ -182,16 +182,15 @@ void Min::init()
   neigh_delay = neighbor->delay;
   neigh_dist_check = neighbor->dist_check;
 
-  if (neigh_every != 1 || neigh_delay != 0 || neigh_dist_check != 1) {
+  if ((neigh_every != 1) || (neigh_delay != 0)) {
     if (comm->me == 0)
-      error->warning(FLERR, "Using 'neigh_modify every 1 delay 0 check"
+      utils::logmesg(lmp, "Switching to 'neigh_modify every 1 delay 0 check yes' "
-                     " yes' setting during minimization");
+                     "setting during minimization\n");
+    neighbor->every = 1;
+    neighbor->delay = 0;
+    neighbor->dist_check = 1;
   }
 
-  neighbor->every = 1;
-  neighbor->delay = 0;
-  neighbor->dist_check = 1;
-
   niter = neval = 0;
 
   // store timestep size (important for variable timestep minimizer)
@@ -751,11 +750,11 @@ void Min::modify_params(int narg, char **arg)
 
 void Min::ev_setup()
 {
-  delete [] elist_global;
+  delete[] elist_global;
-  delete [] elist_atom;
+  delete[] elist_atom;
-  delete [] vlist_global;
+  delete[] vlist_global;
-  delete [] vlist_atom;
+  delete[] vlist_atom;
-  delete [] cvlist_atom;
+  delete[] cvlist_atom;
   elist_global = elist_atom = nullptr;
   vlist_global = vlist_atom = cvlist_atom = nullptr;
 

src/modify.cpp

@@ -810,7 +810,7 @@ void Modify::reset_grid()
 
 Fix *Modify::add_fix(int narg, char **arg, int trysuffix)
 {
-  if (narg < 3) error->all(FLERR, "Illegal fix command");
+  if (narg < 3) utils::missing_cmd_args(FLERR, "fix", error);
 
   // cannot define fix before box exists unless style is in exception list
   // don't like this way of checking for exceptions by adding fixes to list,
@@ -1001,7 +1001,7 @@ Fix *Modify::replace_fix(const char *replaceID, int narg, char **arg, int trysuf
   // change ID, igroup, style of fix being replaced to match new fix
   // requires some error checking on arguments for new fix
 
-  if (narg < 3) error->all(FLERR, "Illegal replace_fix invocation");
+  if (narg < 3) error->all(FLERR, "Not enough arguments for replace_fix invocation");
   if (get_fix_by_id(arg[0])) error->all(FLERR, "Replace_fix ID {} is already in use", arg[0]);
 
   delete[] oldfix->id;
@@ -1039,16 +1039,11 @@ Fix *Modify::replace_fix(const std::string &oldfix, const std::string &fixcmd, i
 
 void Modify::modify_fix(int narg, char **arg)
 {
-  if (narg < 2) error->all(FLERR, "Illegal fix_modify command");
+  if (narg < 2) utils::missing_cmd_args(FLERR, "fix_modify", error);
 
-  // lookup Fix ID
+  auto ifix = get_fix_by_id(arg[0]);
-
+  if (!ifix) error->all(FLERR, "Could not find fix_modify ID {}", arg[0]);
-  int ifix;
+  ifix->modify_params(narg - 1, &arg[1]);
-  for (ifix = 0; ifix < nfix; ifix++)
-    if (strcmp(arg[0], fix[ifix]->id) == 0) break;
-  if (ifix == nfix) error->all(FLERR, "Could not find fix_modify ID {}", arg[0]);
-
-  fix[ifix]->modify_params(narg - 1, &arg[1]);
 }
 
 /* ----------------------------------------------------------------------
@@ -1239,13 +1234,11 @@ int Modify::check_rigid_list_overlap(int *select)
 
 Compute *Modify::add_compute(int narg, char **arg, int trysuffix)
 {
-  if (narg < 3) error->all(FLERR, "Illegal compute command");
+  if (narg < 3) utils::missing_cmd_args(FLERR, "compute", error);
 
   // error check
 
-  for (int icompute = 0; icompute < ncompute; icompute++)
+  if (get_compute_by_id(arg[0])) error->all(FLERR, "Reuse of compute ID '{}'", arg[0]);
-    if (strcmp(arg[0], compute[icompute]->id) == 0)
-      error->all(FLERR, "Reuse of compute ID '{}'", arg[0]);
 
   // extend Compute list if necessary
 
@@ -1312,16 +1305,13 @@ Compute *Modify::add_compute(const std::string &computecmd, int trysuffix)
 
 void Modify::modify_compute(int narg, char **arg)
 {
-  if (narg < 2) error->all(FLERR, "Illegal compute_modify command");
+  if (narg < 2) utils::missing_cmd_args(FLERR, "compute_modify",error);
 
   // lookup Compute ID
 
-  int icompute;
+  auto icompute = get_compute_by_id(arg[0]);
-  for (icompute = 0; icompute < ncompute; icompute++)
+  if (!icompute) error->all(FLERR, "Could not find compute_modify ID {}", arg[0]);
-    if (strcmp(arg[0], compute[icompute]->id) == 0) break;
+  icompute->modify_params(narg - 1, &arg[1]);
-  if (icompute == ncompute) error->all(FLERR, "Could not find compute_modify ID {}", arg[0]);
-
-  compute[icompute]->modify_params(narg - 1, &arg[1]);
 }
 
 /* ----------------------------------------------------------------------

src/neighbor.cpp

@@ -89,11 +89,6 @@ static const char cite_neigh_multi[] =
   " volume =  179,\n"
   " pages =   {320--329}\n"
   "}\n\n"
-  "@article{Stratford2018,\n"
-  " author = {Stratford, Kevin and Shire, Tom and Hanley, Kevin},\n"
-  " title = {Implementation of multi-level contact detection in LAMMPS},\n"
-  " year = {2018}\n"
-  "}\n\n"
   "@article{Shire2020,\n"
   " author = {Shire, Tom and Hanley, Kevin J. and Stratford, Kevin},\n"
   " title = {{DEM} Simulations of Polydisperse Media: Efficient Contact\n"
@@ -1718,7 +1713,7 @@ void Neighbor::print_pairwise_info()
   }
 
   std::string out = "Neighbor list info ...\n";
-  out += fmt::format("  update every {} steps, delay {} steps, check {}\n",
+  out += fmt::format("  update: every = {} steps, delay = {} steps, check = {}\n",
                      every,delay,dist_check ? "yes" : "no");
   out += fmt::format("  max neighbors/atom: {}, page size: {}\n",
                      oneatom, pgsize);

src/output.cpp

@@ -783,6 +783,7 @@ Dump *Output::add_dump(int narg, char **arg)
   next_dump[idump] = 0;
 
   ndump++;
+  dump_list = std::vector<Dump *>(dump, dump + ndump);
   return dump[idump];
 }
 
@@ -792,16 +793,13 @@ Dump *Output::add_dump(int narg, char **arg)
 
 void Output::modify_dump(int narg, char **arg)
 {
-  if (narg < 1) utils::missing_cmd_args(FLERR, "dump_modify",error);
+  if (narg < 2) utils::missing_cmd_args(FLERR, "dump_modify",error);
 
   // find which dump it is
 
-  int idump;
+  auto idump = get_dump_by_id(arg[0]);
-  for (idump = 0; idump < ndump; idump++)
+  if (!idump) error->all(FLERR,"Could not find dump_modify ID: {}", arg[0]);
-    if (strcmp(arg[0],dump[idump]->id) == 0) break;
+  idump->modify_params(narg-1,&arg[1]);
-  if (idump == ndump) error->all(FLERR,"Could not find dump_modify ID: {}", arg[0]);
-
-  dump[idump]->modify_params(narg-1,&arg[1]);
 }
 
 /* ----------------------------------------------------------------------
@@ -833,21 +831,7 @@ void Output::delete_dump(const std::string &id)
     ivar_dump[i-1] = ivar_dump[i];
   }
   ndump--;
-}
+  dump_list = std::vector<Dump *>(dump, dump + ndump);
-
-/* ----------------------------------------------------------------------
-   find a dump by ID
-   return index of dump or -1 if not found
-------------------------------------------------------------------------- */
-
-int Output::find_dump(const char *id)
-{
-  if (id == nullptr) return -1;
-  int idump;
-  for (idump = 0; idump < ndump; idump++)
-    if (strcmp(id,dump[idump]->id) == 0) break;
-  if (idump == ndump) return -1;
-  return idump;
 }
 
 /* ----------------------------------------------------------------------
@@ -855,13 +839,23 @@ int Output::find_dump(const char *id)
    return pointer to dump
 ------------------------------------------------------------------------- */
 
-Dump *Output::get_dump_by_id(const std::string &id)
+Dump *Output::get_dump_by_id(const std::string &id) const
 {
   if (id.empty()) return nullptr;
   for (int idump = 0; idump < ndump; idump++) if (id == dump[idump]->id) return dump[idump];
   return nullptr;
 }
 
+/* ----------------------------------------------------------------------
+   return list of dumps as vector
+------------------------------------------------------------------------- */
+
+const std::vector<Dump *> &Output::get_dump_list()
+{
+  dump_list = std::vector<Dump *>(dump, dump + ndump);
+  return dump_list;
+}
+
 /* ----------------------------------------------------------------------
    set thermo output frequency from input script
 ------------------------------------------------------------------------- */

src/output.h

@@ -80,12 +80,17 @@ class Output : protected Pointers {
   void reset_timestep(bigint);    // reset output which depends on timestep
   void reset_dt();                // reset output which depends on timestep size
 
-  Dump *add_dump(int, char **);                 // add a Dump to Dump list
+  Dump *add_dump(int, char **);                       // add a Dump to Dump list
-  void modify_dump(int, char **);               // modify a Dump
+  void modify_dump(int, char **);                     // modify a Dump
-  void delete_dump(const std::string &);        // delete a Dump from Dump list
+  void delete_dump(const std::string &);              // delete a Dump from Dump list
-  int find_dump(const char *);                  // find a Dump ID
+  Dump *get_dump_by_id(const std::string &) const;    // find a Dump by ID
-  Dump *get_dump_by_id(const std::string &);    // find a Dump by ID
+  Dump *get_dump_by_index(int idx) const              // find a Dump by index in Dump list
-  int check_time_dumps(bigint);                 // check if any time dump is output now
+  {
+    return ((idx >= 0) && (idx < ndump)) ? dump[idx] : nullptr;
+  }
+
+  const std::vector<Dump *> &get_dump_list();    // get vector with all dumps
+  int check_time_dumps(bigint);                  // check if any time dump is output now
 
   void set_thermo(int, char **);        // set thermo output freqquency
   void create_thermo(int, char **);     // create a thermo style
@@ -94,6 +99,7 @@ class Output : protected Pointers {
   void memory_usage();    // print out memory usage
 
  private:
+  std::vector<Dump *> dump_list;
   void calculate_next_dump(int, int, bigint);
 };
 

src/read_data.cpp

@@ -17,6 +17,7 @@
 #include "angle.h"
 #include "atom.h"
 #include "atom_vec.h"
+#include "atom_vec_body.h"
 #include "atom_vec_ellipsoid.h"
 #include "atom_vec_line.h"
 #include "atom_vec_tri.h"
@@ -99,30 +100,30 @@ ReadData::ReadData(LAMMPS *lmp) : Command(lmp)
   // pointers to atom styles that store bonus info
 
   nellipsoids = 0;
-  avec_ellipsoid = dynamic_cast<AtomVecEllipsoid *>( atom->style_match("ellipsoid"));
+  avec_ellipsoid = dynamic_cast<AtomVecEllipsoid *>(atom->style_match("ellipsoid"));
   nlines = 0;
-  avec_line = dynamic_cast<AtomVecLine *>( atom->style_match("line"));
+  avec_line = dynamic_cast<AtomVecLine *>(atom->style_match("line"));
   ntris = 0;
-  avec_tri = dynamic_cast<AtomVecTri *>( atom->style_match("tri"));
+  avec_tri = dynamic_cast<AtomVecTri *>(atom->style_match("tri"));
   nbodies = 0;
-  avec_body = (AtomVecBody *) atom->style_match("body");
+  avec_body = dynamic_cast<AtomVecBody *>(atom->style_match("body"));
 }
 
 /* ---------------------------------------------------------------------- */
 
 ReadData::~ReadData()
 {
-  delete [] line;
+  delete[] line;
-  delete [] keyword;
+  delete[] keyword;
-  delete [] style;
+  delete[] style;
-  delete [] buffer;
+  delete[] buffer;
   memory->sfree(coeffarg);
 
   for (int i = 0; i < nfix; i++) {
-    delete [] fix_header[i];
+    delete[] fix_header[i];
-    delete [] fix_section[i];
+    delete[] fix_section[i];
   }
-  memory->destroy(fix_index);
+  memory->sfree(fix_index);
   memory->sfree(fix_header);
   memory->sfree(fix_section);
 }
@@ -285,16 +286,13 @@ void ReadData::command(int narg, char **arg)
     } else if (strcmp(arg[iarg],"fix") == 0) {
       if (iarg+4 > narg)
         error->all(FLERR,"Illegal read_data command");
-      memory->grow(fix_index,nfix+1,"read_data:fix_index");
+      fix_index = (Fix **)memory->srealloc(fix_index,(nfix+1)*sizeof(Fix *),"read_data:fix_index");
-      fix_header = (char **) memory->srealloc(fix_header,(nfix+1)*sizeof(char *),
+      fix_header = (char **) memory->srealloc(fix_header,(nfix+1)*sizeof(char *),"read_data:fix_header");
-                                              "read_data:fix_header");
+      fix_section = (char **) memory->srealloc(fix_section,(nfix+1)*sizeof(char *),"read_data:fix_section");
-      fix_section = (char **) memory->srealloc(fix_section,(nfix+1)*sizeof(char *),
-                                               "read_data:fix_section");
       if (is_data_section(arg[iarg+3]))
-        error->all(FLERR,"Custom data section name {} for fix {} collides with existing "
+        error->all(FLERR,"Custom data section name {} for fix {} collides with existing data section",arg[iarg+3],arg[iarg+1]);
-                   "data section",arg[iarg+3],arg[iarg+1]);
+      fix_index[nfix] = modify->get_fix_by_id(arg[iarg+1]);
-      fix_index[nfix] = modify->find_fix(arg[iarg+1]);
+      if (!fix_index[nfix]) error->all(FLERR,"Fix ID {} for read_data does not exist",arg[iarg+1]);
-      if (fix_index[nfix] < 0) error->all(FLERR,"Fix ID for read_data does not exist");
       if (strcmp(arg[iarg+2],"NULL") == 0) fix_header[nfix] = nullptr;
       else fix_header[nfix] = utils::strdup(arg[iarg+2]);
       if (strcmp(arg[iarg+3],"NULL") == 0) fix_section[nfix] = utils::strdup(arg[iarg+1]);
@@ -743,8 +741,7 @@ void ReadData::command(int narg, char **arg)
         for (i = 0; i < nfix; i++)
           if (strcmp(keyword,fix_section[i]) == 0) {
             if (firstpass) fix(fix_index[i],keyword);
-            else skip_lines(modify->fix[fix_index[i]]->
+            else skip_lines(fix_index[i]->read_data_skip_lines(keyword));
-                            read_data_skip_lines(keyword));
             break;
           }
         if (i == nfix)
@@ -1011,7 +1008,7 @@ void ReadData::header(int firstpass)
       for (n = 0; n < nfix; n++) {
         if (!fix_header[n]) continue;
         if (strstr(line,fix_header[n])) {
-          modify->fix[fix_index[n]]->read_data_header(line);
+          fix_index[n]->read_data_header(line);
           break;
         }
       }
@@ -1757,7 +1754,7 @@ void ReadData::mass()
     atom->set_mass(FLERR,buf,toffset);
     buf = next + 1;
   }
-  delete [] original;
+  delete[] original;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1780,7 +1777,7 @@ void ReadData::paircoeffs()
     force->pair->coeff(ncoeffarg,coeffarg);
     buf = next + 1;
   }
-  delete [] original;
+  delete[] original;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1807,7 +1804,7 @@ void ReadData::pairIJcoeffs()
       force->pair->coeff(ncoeffarg,coeffarg);
       buf = next + 1;
     }
-  delete [] original;
+  delete[] original;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1832,7 +1829,7 @@ void ReadData::bondcoeffs()
     force->bond->coeff(ncoeffarg,coeffarg);
     buf = next + 1;
   }
-  delete [] original;
+  delete[] original;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1859,7 +1856,7 @@ void ReadData::anglecoeffs(int which)
     force->angle->coeff(ncoeffarg,coeffarg);
     buf = next + 1;
   }
-  delete [] original;
+  delete[] original;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1889,7 +1886,7 @@ void ReadData::dihedralcoeffs(int which)
     force->dihedral->coeff(ncoeffarg,coeffarg);
     buf = next + 1;
   }
-  delete [] original;
+  delete[] original;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1914,7 +1911,7 @@ void ReadData::impropercoeffs(int which)
     force->improper->coeff(ncoeffarg,coeffarg);
     buf = next + 1;
   }
-  delete [] original;
+  delete[] original;
 }
 
 /* ----------------------------------------------------------------------
@@ -1922,18 +1919,18 @@ void ReadData::impropercoeffs(int which)
    n = index of fix
 ------------------------------------------------------------------------- */
 
-void ReadData::fix(int ifix, char *keyword)
+void ReadData::fix(Fix *ifix, char *keyword)
 {
   int nchunk,eof;
 
-  bigint nline = modify->fix[ifix]->read_data_skip_lines(keyword);
+  bigint nline = ifix->read_data_skip_lines(keyword);
 
   bigint nread = 0;
   while (nread < nline) {
     nchunk = MIN(nline-nread,CHUNK);
     eof = utils::read_lines_from_file(fp,nchunk,MAXLINE,buffer,me,world);
-    if (eof) error->all(FLERR,"Unexpected end of data file");
+    if (eof) error->all(FLERR,"Unexpected end of data file while reading section {}",keyword);
-    modify->fix[ifix]->read_data_section(keyword,nchunk,buffer,id_offset);
+    ifix->read_data_section(keyword,nchunk,buffer,id_offset);
     nread += nchunk;
   }
 }

src/read_data.h

@@ -23,7 +23,7 @@ CommandStyle(read_data,ReadData);
 #include "command.h"
 
 namespace LAMMPS_NS {
-
+class Fix;
 class ReadData : public Command {
  public:
   ReadData(class LAMMPS *);
@@ -73,7 +73,7 @@ class ReadData : public Command {
   int groupbit;
 
   int nfix;
-  int *fix_index;
+  Fix **fix_index;
   char **fix_header;
   char **fix_section;
 
@@ -108,7 +108,7 @@ class ReadData : public Command {
   void dihedralcoeffs(int);
   void impropercoeffs(int);
 
-  void fix(int, char *);
+  void fix(Fix *, char *);
 };
 
 }    // namespace LAMMPS_NS

src/region.cpp

@@ -325,7 +325,7 @@ void Region::options(int narg, char **arg)
         error->all(FLERR, "Illegal region units: {}", arg[iarg+1]);
       iarg += 2;
     } else if (strcmp(arg[iarg], "side") == 0) {
-      if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "region side", error); 
+      if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "region side", error);
       if (strcmp(arg[iarg + 1], "in") == 0)
         interior = 1;
       else if (strcmp(arg[iarg + 1], "out") == 0)

unittest/commands/test_groups.cpp

@@ -122,13 +122,13 @@ TEST_F(GroupTest, EmptyDelete)
     TEST_FAILURE(".*ERROR: Illegal group command.*", command("group new2 delete xxx"););
     TEST_FAILURE(".*ERROR: Cannot delete group all.*", command("group all delete"););
     TEST_FAILURE(".*ERROR: Could not find group delete.*", command("group new0 delete"););
-    TEST_FAILURE(".*ERROR: Cannot delete group currently used by a fix.*",
+    TEST_FAILURE(".*ERROR: Cannot delete group new2 currently used by fix.*",
                  command("group new2 delete"););
-    TEST_FAILURE(".*ERROR: Cannot delete group currently used by a compute.*",
+    TEST_FAILURE(".*ERROR: Cannot delete group new3 currently used by compute.*",
                  command("group new3 delete"););
-    TEST_FAILURE(".*ERROR: Cannot delete group currently used by a dump.*",
+    TEST_FAILURE(".*ERROR: Cannot delete group new4 currently used by dump.*",
                  command("group new4 delete"););
-    TEST_FAILURE(".*ERROR: Cannot delete group currently used by atom_modify.*",
+    TEST_FAILURE(".*ERROR: Cannot delete group new5 currently used by atom_modify.*",
                  command("group new5 delete"););
 }