From f557dfdebbd413441f372b64ea6f96471884ac97 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Thu, 11 Apr 2024 23:36:11 -0400
Subject: [PATCH 01/17] mention typelabel paper

---
 doc/src/Howto_type_labels.rst               | 25 +++++++++++++--------
 doc/src/labelmap.rst                        |  9 +++++++-
 doc/utils/sphinx-config/false_positives.txt |  1 +
 3 files changed, 25 insertions(+), 10 deletions(-)

diff --git a/doc/src/Howto_type_labels.rst b/doc/src/Howto_type_labels.rst
index 8f03f45ef9..d7460d50a7 100644
--- a/doc/src/Howto_type_labels.rst
+++ b/doc/src/Howto_type_labels.rst
@@ -14,16 +14,17 @@ wherever they appear in LAMMPS input or output files.  The total number
 Ntypes for each interaction is "locked in" when the simulation box
 is created.
 
-A recent addition to LAMMPS is the option to use strings - referred
-to as type labels - as an alternative.  Using type labels instead of
+A recent addition to LAMMPS is the option to use strings - referred to
+as type labels - as an alternative.  Using type labels instead of
 numeric types can be advantageous in various scenarios.  For example,
-type labels can make inputs more readable and generic (i.e. usable through
-the :doc:`include command <include>` for different systems with different
-numerical values assigned to types.  This generality also applies to
-other inputs like data files read by :doc:`read_data <read_data>` or
-molecule template files read by the :doc:`molecule <molecule>`
-command.  See below for a list of other commands that can use
-type labels in different ways.
+type labels can make inputs more readable and generic (i.e. usable
+through the :doc:`include command <include>` for different systems with
+different numerical values assigned to types.  This generality also
+applies to other inputs like data files read by :doc:`read_data
+<read_data>` or molecule template files read by the :doc:`molecule
+<molecule>` command.  A discussion of the current type label support can
+be found in :ref:`(Gissinger) <Typelabel24>`.  See below for a list of
+other commands that can use type labels in different ways.
 
 LAMMPS will *internally* continue to use numeric types, which means
 that many previous restrictions still apply.  For example, the total
@@ -124,3 +125,9 @@ between the files.  The creation of simulation-ready reaction templates
 for :doc:`fix bond/react <fix_bond_react>` is much simpler when using
 type labels, and results in templates that can be used without
 modification in multiple simulations or different systems.
+
+-----------
+
+.. _Typelabel24:
+
+**(Gissinger)** J. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M. Choi, D. Karls, A. Stukowski, W, Im, H. Heinz, A. Kohlmeyer, and E. Tadmor, J Phys Chem B, 128, 3282-3297 (2024)
diff --git a/doc/src/labelmap.rst b/doc/src/labelmap.rst
index 9e3d705101..5e0d1e92e9 100644
--- a/doc/src/labelmap.rst
+++ b/doc/src/labelmap.rst
@@ -43,7 +43,8 @@ The label map can also be defined by the :doc:`read_data <read_data>`
 command when it reads these sections in a data file: Atom Type Labels,
 Bond Type Labels, etc.  See the :doc:`Howto type labels
 <Howto_type_labels>` doc page for a general discussion of how type
-labels can be used.
+labels can be used.  See :ref:`(Gissinger) <Typelabel>` for a discussion
+of the type label implementation in LAMMPS and its uses.
 
 Valid type labels can contain any alphanumeric character, but must not
 start with a number, a '#', or a '*' character.  They can contain other
@@ -98,3 +99,9 @@ Default
 """""""
 
 none
+
+-----------
+
+.. _Typelabel:
+
+**(Gissinger)** J. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M. Choi, D. Karls, A. Stukowski, W, Im, H. Heinz, A. Kohlmeyer, and E. Tadmor, J Phys Chem B, 128, 3282-3297 (2024)
diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt
index 88ac0e3140..cf4e48feff 100644
--- a/doc/utils/sphinx-config/false_positives.txt
+++ b/doc/utils/sphinx-config/false_positives.txt
@@ -476,6 +476,7 @@ ChiralIDs
 chirality
 Cho
 Chodera
+Choi
 ChooseOffset
 chris
 Christoph

From 441a521d8bacf91c0e0091ed437c5ca5879ad896 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Thu, 11 Apr 2024 23:52:34 -0400
Subject: [PATCH 02/17] initialize all class pointers to null

---
 src/EXTRA-PAIR/pair_pedone.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/EXTRA-PAIR/pair_pedone.cpp b/src/EXTRA-PAIR/pair_pedone.cpp
index c5f7f49117..9b8ce451d9 100644
--- a/src/EXTRA-PAIR/pair_pedone.cpp
+++ b/src/EXTRA-PAIR/pair_pedone.cpp
@@ -32,8 +32,8 @@ using namespace LAMMPS_NS;
 /* ---------------------------------------------------------------------- */
 
 PairPedone::PairPedone(LAMMPS *lmp) :
-    Pair(lmp), d0(nullptr), alpha(nullptr), r0(nullptr), c0(nullptr), pedone1(nullptr),
-    pedone2(nullptr)
+    Pair(lmp), cut(nullptr), d0(nullptr), alpha(nullptr), r0(nullptr), c0(nullptr),
+    pedone1(nullptr), pedone2(nullptr), offset(nullptr)
 {
   writedata = 1;
 }

From 4bd983ce6aa656418ad4baf79f5eb245ab72ffc9 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Fri, 12 Apr 2024 00:24:21 -0400
Subject: [PATCH 03/17] make warnings in scatter/gather into errors

---
 src/library.cpp | 430 +++++++++++++++++++++++++-----------------------
 1 file changed, 222 insertions(+), 208 deletions(-)

diff --git a/src/library.cpp b/src/library.cpp
index f81f52305e..75d74f4bf8 100644
--- a/src/library.cpp
+++ b/src/library.cpp
@@ -527,7 +527,7 @@ void lammps_file(void *handle, const char *filename)
   BEGIN_CAPTURE
   {
     if (lmp->update->whichflag != 0)
-      lmp->error->all(FLERR, "Library error: issuing LAMMPS commands during a run is not allowed");
+      lmp->error->all(FLERR, "Issuing LAMMPS commands during a run is not allowed");
     else
       lmp->input->file(filename);
   }
@@ -564,8 +564,7 @@ char *lammps_command(void *handle, const char *cmd)
   BEGIN_CAPTURE
   {
     if (lmp->update->whichflag != 0)
-      lmp->error->all(FLERR,"Library error: issuing LAMMPS commands "
-                      "during a run is not allowed.");
+      lmp->error->all(FLERR, "Issuing LAMMPS command during a run is not allowed.");
     else
       result = lmp->input->one(cmd);
   }
@@ -641,7 +640,7 @@ void lammps_commands_string(void *handle, const char *str)
   BEGIN_CAPTURE
   {
     if (lmp->update->whichflag != 0) {
-      lmp->error->all(FLERR,"Library error: issuing LAMMPS command during run");
+      lmp->error->all(FLERR, "Issuing LAMMPS commands during a run is not allowed");
     }
 
     std::size_t cursor = 0;
@@ -949,9 +948,9 @@ void lammps_extract_box(void *handle, double *boxlo, double *boxhi,
   BEGIN_CAPTURE
   {
     // do nothing if box does not yet exist
-    if ((lmp->domain->box_exist == 0)
-        && (lmp->comm->me == 0)) {
-      lmp->error->warning(FLERR,"Calling lammps_extract_box without a box");
+    if (lmp->domain->box_exist == 0) {
+      if (lmp->comm->me == 0)
+        lmp->error->warning(FLERR, "Call to lammps_extract_box() without a box ignored");
       return;
     }
 
@@ -1011,12 +1010,12 @@ void lammps_reset_box(void *handle, double *boxlo, double *boxhi,
   BEGIN_CAPTURE
   {
     if (lmp->atom->natoms > 0)
-      lmp->error->all(FLERR,"Calling lammps_reset_box not supported when atoms exist");
+      lmp->error->all(FLERR, "Calling lammps_reset_box() not supported when atoms exist");
 
     // warn and do nothing if no box exists
     if (lmp->domain->box_exist == 0) {
       if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"Ignoring call to lammps_reset_box without a box");
+        lmp->error->warning(FLERR,"Call to lammps_reset_box() without a box ignored");
       return;
     }
 
@@ -2629,7 +2628,14 @@ x[0][2], x[1][0], x[1][1], x[1][2], x[2][0], :math:`\dots`);
 *natoms*), as queried by :cpp:func:`lammps_get_natoms`,
 :cpp:func:`lammps_extract_global`, or :cpp:func:`lammps_extract_setting`.
 
-This function is not compatible with ``-DLAMMPS_BIGBIG``.
+.. admonition:: Restrictions
+    :class: warning
+
+    This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
+    Atom IDs must be defined and consecutive.
+
+    The total number of atoms must not be more than 2147483647 (max 32-bit signed int).
 
 \endverbatim
  *
@@ -2650,8 +2656,7 @@ This function is not compatible with ``-DLAMMPS_BIGBIG``.
      Allreduce to sum vector into data across all procs
 ------------------------------------------------------------------------- */
 
-void lammps_gather_atoms(void *handle, const char *name, int type, int count,
-                         void *data)
+void lammps_gather_atoms(void *handle, const char *name, int type, int count, void *data)
 {
   auto lmp = (LAMMPS *) handle;
 
@@ -2671,8 +2676,7 @@ void lammps_gather_atoms(void *handle, const char *name, int type, int count,
       flag = 1;
     if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
     if (flag) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"Library error in lammps_gather_atoms");
+      lmp->error->all(FLERR,"lammps_gather_atoms(): Atom-IDs must exist and be consecutive");
       return;
     }
 
@@ -2680,8 +2684,7 @@ void lammps_gather_atoms(void *handle, const char *name, int type, int count,
 
     void *vptr = lmp->atom->extract(name);
     if (vptr == nullptr) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"lammps_gather_atoms: unknown property name");
+      lmp->error->all(FLERR, "lammps_gather_atoms(): unknown property {}", name);
       return;
     }
 
@@ -2756,8 +2759,7 @@ void lammps_gather_atoms(void *handle, const char *name, int type, int count,
       MPI_Allreduce(copy,data,count*natoms,MPI_DOUBLE,MPI_SUM,lmp->world);
       lmp->memory->destroy(copy);
     } else {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"lammps_gather_atoms: unsupported data type");
+      lmp->error->all(FLERR,"lammps_gather_atoms(): unsupported data type");
       return;
     }
 #endif
@@ -2785,10 +2787,17 @@ groups total, but not in order by atom ID (e.g., if *name* is *x* and *count*
 is 3, then *data* might be something like x[10][0], x[10][1], x[10][2],
 x[2][0], x[2][1], x[2][2], x[4][0], :math:`\dots`); *data* must be
 pre-allocated by the caller to length (*count* :math:`\times` *natoms*), as
-queried by :cpp:func:`lammps_get_natoms`,
-:cpp:func:`lammps_extract_global`, or :cpp:func:`lammps_extract_setting`.
+queried by :cpp:func:`lammps_get_natoms`, :cpp:func:`lammps_extract_global`,
+or :cpp:func:`lammps_extract_setting`.
 
-This function is not compatible with ``-DLAMMPS_BIGBIG``.
+.. admonition:: Restrictions
+    :class: warning
+
+    This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
+    Atom IDs must be defined.
+
+    The total number of atoms must not be more than 2147483647 (max 32-bit signed int).
 
 \endverbatim
  *
@@ -2828,8 +2837,7 @@ void lammps_gather_atoms_concat(void *handle, const char *name, int type,
     if (lmp->atom->tag_enable == 0) flag = 1;
     if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
     if (flag) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"Library error in lammps_gather_atoms");
+      lmp->error->all(FLERR,"lammps_gather_atoms_concat(): Atom-IDs must exist");
       return;
     }
 
@@ -2837,8 +2845,7 @@ void lammps_gather_atoms_concat(void *handle, const char *name, int type,
 
     void *vptr = lmp->atom->extract(name);
     if (vptr == nullptr) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"lammps_gather_atoms: unknown property name");
+      lmp->error->all(FLERR,"lammps_gather_atoms_concat(): unknown property {}", name);
       return;
     }
 
@@ -2947,7 +2954,14 @@ x[100][2], x[57][0], x[57][1], x[57][2], x[210][0], :math:`\dots`);
 *data* must be pre-allocated by the caller to length
 (*count* :math:`\times` *ndata*).
 
-This function is not compatible with ``-DLAMMPS_BIGBIG``.
+.. admonition:: Restrictions
+    :class: warning
+
+    This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
+    Atom IDs must be defined and an :doc:`atom map must be enabled <atom_modify>`
+
+    The total number of atoms must not be more than 2147483647 (max 32-bit signed int).
 
 \endverbatim
  *
@@ -2993,16 +3007,13 @@ void lammps_gather_atoms_subset(void *handle, const char *name, int type,
     if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
     if (lmp->atom->map_style == Atom::MAP_NONE) flag = 1;
     if (flag) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"Library error in lammps_gather_atoms_subset: atoms must have mappable ids");
+      lmp->error->all(FLERR,"lammps_gather_atoms_subset(): Atom-IDs must exist and be mapped");
       return;
     }
 
     void *vptr = lmp->atom->extract(name);
     if (vptr == nullptr) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"lammps_gather_atoms_subset: "
-                            "unknown property name");
+      lmp->error->all(FLERR,"lammps_gather_atoms_subset(): unknown property {}", name);
       return;
     }
 
@@ -3110,7 +3121,15 @@ atom ID (e.g., if *name* is *x* and *count* = 3, then
 *data* = {x[0][0], x[0][1], x[0][2], x[1][0], x[1][1], x[1][2], x[2][0],
 :math:`\dots`}); *data* must be of length (*count* :math:`\times` *natoms*).
 
-This function is not compatible with ``-DLAMMPS_BIGBIG``.
+.. admonition:: Restrictions
+    :class: warning
+
+    This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
+    Atom IDs must be defined, must be consecutive, and an
+    :doc:`atom map must be enabled <atom_modify>`
+
+    The total number of atoms must not be more than 2147483647 (max 32-bit signed int).
 
 \endverbatim
  *
@@ -3152,8 +3171,8 @@ void lammps_scatter_atoms(void *handle, const char *name, int type, int count,
     if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
     if (lmp->atom->map_style == Atom::MAP_NONE) flag = 1;
     if (flag) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"Library error in lammps_scatter_atoms: ids must exist, be consecutive, and be mapped");
+      lmp->error->all(FLERR,"lammps_scatter_atoms(): "
+                      "Atom-IDs must exist, be consecutive, and be mapped");
       return;
     }
 
@@ -3161,9 +3180,7 @@ void lammps_scatter_atoms(void *handle, const char *name, int type, int count,
 
     void *vptr = lmp->atom->extract(name);
     if (vptr == nullptr) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,
-                            "lammps_scatter_atoms: unknown property name");
+      lmp->error->all(FLERR, "lammps_scatter_atoms(): unknown property {}", name);
       return;
     }
 
@@ -3244,7 +3261,14 @@ to be the array {x[1][0], x[1][1], x[1][2], x[100][0], x[100][1], x[100][2],
 x[57][0], x[57][1], x[57][2]}, then *count* = 3, *ndata* = 3, and *ids* would
 be {1, 100, 57}.
 
-This function is not compatible with ``-DLAMMPS_BIGBIG``.
+.. admonition:: Restrictions
+    :class: warning
+
+    This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
+    Atom IDs must be defined and an :doc:`atom map must be enabled <atom_modify>`
+
+    The total number of atoms must not be more than 2147483647 (max 32-bit signed int).
 
 \endverbatim
  *
@@ -3301,16 +3325,13 @@ void lammps_scatter_atoms_subset(void *handle, const char *name, int type,
     if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
     if (lmp->atom->map_style == Atom::MAP_NONE) flag = 1;
     if (flag) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"Library error in lammps_scatter_atoms_subset: atoms must have mapped ids");
+      lmp->error->all(FLERR,"lammps_scatter_atoms_subset(): Atom-IDs must exist and be mapped");
       return;
     }
 
     void *vptr = lmp->atom->extract(name);
     if (vptr == nullptr) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,
-                            "lammps_scatter_atoms_subset: unknown property name");
+      lmp->error->all(FLERR, "lammps_scatter_atoms_subset(): unknown property {}", name);
       return;
     }
 
@@ -3849,7 +3870,14 @@ The *data* array will be ordered in groups of *count* values, sorted by atom ID
 This function will return an error if fix or compute data are requested and the fix or compute ID
 given does not have per-atom data.
 
-This function is not compatible with ``-DLAMMPS_BIGBIG``.
+.. admonition:: Restrictions
+    :class: warning
+
+    This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
+    Atom IDs must be defined and must be consecutive.
+
+    The total number of atoms must not be more than 2147483647 (max 32-bit signed int).
 
 \endverbatim
  *
@@ -3899,7 +3927,7 @@ void lammps_gather(void *handle, const char *name, int type, int count, void *da
   BEGIN_CAPTURE
   {
 #if defined(LAMMPS_BIGBIG)
-    lmp->error->all(FLERR,"Library function lammps_gather not compatible with -DLAMMPS_BIGBIG");
+    lmp->error->all(FLERR, "Library function lammps_gather() is not compatible with -DLAMMPS_BIGBIG");
 #else
     int i,j,offset,ltype;
 
@@ -3910,8 +3938,7 @@ void lammps_gather(void *handle, const char *name, int type, int count, void *da
       flag = 1;
     if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
     if (flag) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"Library error in lammps_gather");
+      lmp->error->all(FLERR,"lammps_gather(): Atom-IDs must exist, and be consecutive");
       return;
     }
 
@@ -3921,28 +3948,25 @@ void lammps_gather(void *handle, const char *name, int type, int count, void *da
     // fix
 
     if (vptr==nullptr && utils::strmatch(name,"^f_")) {
+      const char *fixid = name+2;
 
-      auto fix = lmp->modify->get_fix_by_id(&name[2]);
+      auto fix = lmp->modify->get_fix_by_id(fixid);
       if (!fix) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather: unknown fix id");
+        lmp->error->all(FLERR,"lammps_gather(): unknown fix id {}", fixid);
         return;
       }
 
       if (fix->peratom_flag == 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather: fix does not return peratom data");
+        lmp->error->all(FLERR,"lammps_gather(): fix {} does not return peratom data", fixid);
         return;
       }
       if ((count > 1) && (fix->size_peratom_cols != count)) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather: count != values peratom for fix");
+        lmp->error->all(FLERR,"lammps_gather: count != values peratom for fix {}", fixid);
         return;
       }
 
       if (lmp->update->ntimestep % fix->peratom_freq) {
-        if (lmp->comm->me == 0)
-          lmp->error->all(FLERR,"lammps_gather: fix not computed at compatible time");
+        lmp->error->all(FLERR,"lammps_gather: fix {} not computed at compatible time", fixid);
         return;
       }
 
@@ -3953,22 +3977,19 @@ void lammps_gather(void *handle, const char *name, int type, int count, void *da
     // compute
 
     if (vptr==nullptr && utils::strmatch(name,"^c_")) {
-
-      auto compute = lmp->modify->get_compute_by_id(&name[2]);
+      const char *compid = name+2;
+      auto compute = lmp->modify->get_compute_by_id(compid);
       if (!compute) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather: unknown compute id");
+        lmp->error->all(FLERR,"lammps_gather(): unknown compute id {}", compid);
         return;
       }
 
       if (compute->peratom_flag == 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather: compute does not return peratom data");
+        lmp->error->all(FLERR,"lammps_gather(): compute {} does not return peratom data", compid);
         return;
       }
       if ((count > 1) && (compute->size_peratom_cols != count)) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather: count != values peratom for compute");
+        lmp->error->all(FLERR,"lammps_gather(): count != values peratom for compute {}", compid);
         return;
       }
 
@@ -3984,28 +4005,26 @@ void lammps_gather(void *handle, const char *name, int type, int count, void *da
     if ((vptr == nullptr) && utils::strmatch(name,"^[id]2?_")) {
 
       int idx,icol;
-      if (utils::strmatch(name,"^[id]_")) idx = lmp->atom->find_custom(&name[2],ltype,icol);
-      else idx = lmp->atom->find_custom(&name[3],ltype,icol);
+      const char *propid;
+      if (utils::strmatch(name,"^[id]_")) propid = name+2;
+      else propid = name+3;
+      idx = lmp->atom->find_custom(propid,ltype,icol);
 
       if (idx < 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather: unknown property/atom id");
+        lmp->error->all(FLERR,"lammps_gather(): unknown property/atom id {}", propid);
         return;
       }
 
       if (ltype != type) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather: mismatch property/atom type");
+        lmp->error->all(FLERR,"lammps_gather(): mismatch property/atom type for {}", propid);
         return;
       }
       if ((count == 1) && (icol != 0)) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather: mismatch property/atom count");
+        lmp->error->all(FLERR,"lammps_gather(): mismatch property/atom count for {}", propid);
         return;
       }
       if ((count > 1) && (icol != count)) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather: mismatch property/atom count");
+        lmp->error->all(FLERR,"lammps_gather(): mismatch property/atom count for {}", propid);
         return;
       }
 
@@ -4021,8 +4040,7 @@ void lammps_gather(void *handle, const char *name, int type, int count, void *da
     // no match
 
     if (vptr == nullptr) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"lammps_gather: undefined property name");
+      lmp->error->all(FLERR,"lammps_gather(): undefined property {}", name);
       return;
     }
 
@@ -4122,7 +4140,14 @@ pre-allocated by the caller to length (*count* :math:`\times` *natoms*), as quer
 :cpp:func:`lammps_get_natoms`, :cpp:func:`lammps_extract_global`, or
 :cpp:func:`lammps_extract_setting`.
 
-This function is not compatible with ``-DLAMMPS_BIGBIG``.
+.. admonition:: Restrictions
+    :class: warning
+
+    This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
+    Atom IDs must be defined.
+
+    The total number of atoms must be less than 2147483647 (max 32-bit signed int).
 
 \endverbatim
  *
@@ -4173,8 +4198,8 @@ void lammps_gather_concat(void *handle, const char *name, int type, int count,
   BEGIN_CAPTURE
   {
 #if defined(LAMMPS_BIGBIG)
-    lmp->error->all(FLERR,"Library function lammps_gather_concat"
-                          " not compatible with -DLAMMPS_BIGBIG");
+    lmp->error->all(FLERR,"Library function lammps_gather_concat()"
+                          " is not compatible with -DLAMMPS_BIGBIG");
 #else
     int i,offset,ltype;
 
@@ -4184,8 +4209,7 @@ void lammps_gather_concat(void *handle, const char *name, int type, int count,
     if (lmp->atom->tag_enable == 0) flag = 1;
     if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
     if (flag) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"Library error in lammps_gather_concat");
+      lmp->error->all(FLERR,"lammps_gather_concat(): atom-IDs must exist");
       return;
     }
 
@@ -4195,27 +4219,24 @@ void lammps_gather_concat(void *handle, const char *name, int type, int count,
     // fix
 
     if (vptr==nullptr && utils::strmatch(name,"^f_")) {
-
-      auto fix = lmp->modify->get_fix_by_id(&name[2]);
+      const char *fixid = name+2;
+      auto fix = lmp->modify->get_fix_by_id(fixid);
       if (!fix) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_concat: unknown fix id");
+        lmp->error->all(FLERR,"lammps_gather_concat(): unknown fix id {}", fixid);
         return;
       }
 
       if (fix->peratom_flag == 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_concat: fix does not return peratom data");
+        lmp->error->all(FLERR,"lammps_gather_concat(): fix {} does not return peratom data", fixid);
         return;
       }
       if ((count > 1) && (fix->size_peratom_cols != count)) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_concat: count != values peratom for fix");
+        lmp->error->all(FLERR,"lammps_gather_concat(): count != values peratom for fix {}", fixid);
         return;
       }
       if (lmp->update->ntimestep % fix->peratom_freq) {
-        if (lmp->comm->me == 0)
-          lmp->error->all(FLERR,"lammps_gather_concat: fix not computed at compatible time");
+        lmp->error->all(FLERR,"lammps_gather_concat(): fix {} not computed at compatible time",
+                        fixid);
         return;
       }
 
@@ -4227,21 +4248,21 @@ void lammps_gather_concat(void *handle, const char *name, int type, int count,
 
     if (vptr==nullptr && utils::strmatch(name,"^c_")) {
 
-      auto compute = lmp->modify->get_compute_by_id(&name[2]);
+      const char *compid = name + 2;
+      auto compute = lmp->modify->get_compute_by_id(compid);
       if (!compute) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_concat: unknown compute id");
+        lmp->error->all(FLERR,"lammps_gather_concat(): unknown compute id {}", compid);
         return;
       }
 
       if (compute->peratom_flag == 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_concat: compute does not return peratom data");
+        lmp->error->all(FLERR,"lammps_gather_concat(): compute {} does not return peratom data",
+                        compid);
         return;
       }
       if ((count > 1) && (compute->size_peratom_cols != count)) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_concat: count != values peratom for compute");
+        lmp->error->all(FLERR,"lammps_gather_concat(): count != values peratom for compute {}",
+          compid);
         return;
       }
 
@@ -4257,28 +4278,26 @@ void lammps_gather_concat(void *handle, const char *name, int type, int count,
     if ((vptr==nullptr) && utils::strmatch(name,"^[id]2?_")) {
 
       int idx,icol;
-      if (utils::strmatch(name,"^[id]_")) idx = lmp->atom->find_custom(&name[2],ltype,icol);
-      else idx = lmp->atom->find_custom(&name[3],ltype,icol);
+      const char *propid;
+      if (utils::strmatch(name,"^[id]_")) propid = name + 2;
+      else propid = name + 3;
+      idx = lmp->atom->find_custom(propid,ltype,icol);
 
       if (idx < 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_concat: unknown property/atom id");
+        lmp->error->all(FLERR,"lammps_gather_concat(): unknown property/atom id {}", propid);
         return;
       }
 
       if (ltype != type) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_concat: mismatch property/atom type");
+        lmp->error->all(FLERR,"lammps_gather_concat(): mismatch property/atom {} type", propid);
         return;
       }
       if ((count == 1) && (icol != 0)) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_concat: mismatch property/atom count");
+        lmp->error->all(FLERR,"lammps_gather_concat(): mismatch property/atom {} count", propid);
         return;
       }
       if ((count > 1) && (icol != count)) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_concat: mismatch property/atom count");
+        lmp->error->all(FLERR,"lammps_gather_concat(): mismatch property/atom {} count", propid);
         return;
       }
 
@@ -4294,8 +4313,7 @@ void lammps_gather_concat(void *handle, const char *name, int type, int count,
     // no match
 
     if (vptr == nullptr) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"lammps_gather_concat: undefined property name");
+      lmp->error->all(FLERR,"lammps_gather_concat(): undefined property {}", name);
       return;
     }
 
@@ -4405,7 +4423,14 @@ look like {x[100][0], x[100][1], x[100][2], x[57][0], x[57][1], x[57][2], x[210]
 :math:`\dots`}); *ids* must be provided by the user with length *ndata*, and *data* must be
 pre-allocated by the caller to length (*count*\ :math:`{}\times{}`\ *ndata*).
 
-This function is not compatible with ``-DLAMMPS_BIGBIG``.
+.. admonition:: Restrictions
+    :class: warning
+
+    This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
+    Atom IDs must be defined and an :doc:`atom map must be enabled <atom_modify>`
+
+    The total number of atoms must not be more than 2147483647 (max 32-bit signed int).
 
 \endverbatim
  *
@@ -4471,8 +4496,7 @@ void lammps_gather_subset(void *handle, const char *name, int type, int count,
     if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
     if (lmp->atom->map_style == Atom::MAP_NONE) flag = 1;
     if (flag) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"Library error in lammps_gather_subset");
+      lmp->error->all (FLERR,"lammps_gather_subset(): atom IDs must be enabled and mapped");
       return;
     }
 
@@ -4481,26 +4505,24 @@ void lammps_gather_subset(void *handle, const char *name, int type, int count,
     // fix
 
     if (vptr==nullptr && utils::strmatch(name,"^f_")) {
+      const char *fixid = name + 2;
 
-      auto fix = lmp->modify->get_fix_by_id(&name[2]);
+      auto fix = lmp->modify->get_fix_by_id(fixid);
       if (!fix) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_subset: unknown fix id");
+        lmp->error->all(FLERR,"lammps_gather_subset(): unknown fix id {}", fixid);
         return;
       }
 
       if (fix->peratom_flag == 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_subset: fix does not return peratom data");
+        lmp->error->all(FLERR,"lammps_gather_subset(): fix {} does not return peratom data", fixid);
         return;
       }
       if ((count > 1) && (fix->size_peratom_cols != count)) {
-        lmp->error->warning(FLERR,"lammps_gather_subset: count != values peratom for fix");
+        lmp->error->all(FLERR,"lammps_gather_subset(): count != values peratom for fix {}", fixid);
         return;
       }
       if (lmp->update->ntimestep % fix->peratom_freq) {
-        if (lmp->comm->me == 0)
-          lmp->error->all(FLERR,"lammps_gather_subset: fix not computed at compatible time");
+        lmp->error->all(FLERR,"lammps_gather_subset(): fix {} not computed at compatible time", fixid);
         return;
       }
 
@@ -4511,22 +4533,21 @@ void lammps_gather_subset(void *handle, const char *name, int type, int count,
     // compute
 
     if (vptr==nullptr && utils::strmatch(name,"^c_")) {
-
-      auto compute = lmp->modify->get_compute_by_id(&name[2]);
+      const char *compid = name + 2;
+      auto compute = lmp->modify->get_compute_by_id(compid);
       if (!compute) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_subset: unknown compute id");
+        lmp->error->all(FLERR,"lammps_gather_subset(): unknown compute id {}", compid);
         return;
       }
 
       if (compute->peratom_flag == 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_subset: compute does not return peratom data");
+        lmp->error->all(FLERR,"lammps_gather_subset(): compute {} does not return peratom data",
+                        compid);
         return;
       }
       if ((count > 1) && (compute->size_peratom_cols != count)) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_subset: count != values peratom for compute");
+        lmp->error->all(FLERR,"lammps_gather_subset(): count != values peratom for compute {}",
+                        compid);
         return;
       }
 
@@ -4542,28 +4563,26 @@ void lammps_gather_subset(void *handle, const char *name, int type, int count,
     if ((vptr == nullptr) && utils::strmatch(name,"^[id]2?_")) {
 
       int idx,icol;
-      if (utils::strmatch(name,"^[id]_")) idx = lmp->atom->find_custom(&name[2],ltype,icol);
-      else idx = lmp->atom->find_custom(&name[3],ltype,icol);
+      const char *propid;
+      if (utils::strmatch(name,"^[id]_")) propid = name + 2;
+      else propid = name + 3;
+      idx = lmp->atom->find_custom(propid,ltype,icol);
 
       if (idx < 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_subset: unknown property/atom id");
+        lmp->error->all(FLERR,"lammps_gather_subset(): unknown property/atom id {}", propid);
         return;
       }
 
       if (ltype != type) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_subset: mismatch property/atom type");
+        lmp->error->all(FLERR,"lammps_gather_subset(): mismatch property/atom {} type", propid);
         return;
       }
       if (count == 1 && icol != 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_subset: mismatch property/atom count");
+        lmp->error->all(FLERR,"lammps_gather_subset(): mismatch property/atom {} count", propid);
         return;
       }
       if (count > 1 && icol != count) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather_subset: mismatch property/atom count");
+        lmp->error->all(FLERR,"lammps_gather_subset(): mismatch property/atom {} count", propid);
         return;
       }
 
@@ -4579,8 +4598,7 @@ void lammps_gather_subset(void *handle, const char *name, int type, int count,
     // no match
 
     if (vptr == nullptr) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"lammps_gather_subset: undefined property name");
+      lmp->error->all(FLERR,"lammps_gather_subset(): undefined property {}", name);
       return;
     }
 
@@ -4688,7 +4706,15 @@ The *data* array needs to be ordered in groups of *count* values, sorted by atom
 *name* is *x* and *count* = 3, then *data* = {x[0][0], x[0][1], x[0][2], x[1][0], x[1][1],
 x[1][2], x[2][0], :math:`\dots`}); *data* must be of length (*count* :math:`\times` *natoms*).
 
-This function is not compatible with ``-DLAMMPS_BIGBIG``.
+.. admonition:: Restrictions
+    :class: warning
+
+    This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
+    Atom IDs must be defined, must be consecutive, and an
+    :doc:`atom map must be enabled <atom_modify>`
+
+    The total number of atoms must not be more than 2147483647 (max 32-bit signed int).
 
 \endverbatim
  *
@@ -4751,8 +4777,7 @@ void lammps_scatter(void *handle, const char *name, int type, int count,
     if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
     if (lmp->atom->map_style == Atom::MAP_NONE) flag = 1;
     if (flag) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"Library error in lammps_scatter");
+      lmp->error->all(FLERR,"lammps_scatter(): atom IDs must be defined, consecutive, and mapped");
       return;
     }
 
@@ -4762,22 +4787,19 @@ void lammps_scatter(void *handle, const char *name, int type, int count,
     // fix
 
     if (vptr==nullptr && utils::strmatch(name,"^f_")) {
-
-      auto fix = lmp->modify->get_fix_by_id(&name[2]);
+      const char *fixid = name + 2;
+      auto fix = lmp->modify->get_fix_by_id(fixid);
       if (!fix) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter: unknown fix id");
+        lmp->error->all(FLERR,"lammps_scatter(): unknown fix id {}", fixid);
         return;
       }
 
       if (fix->peratom_flag == 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter: fix does not return peratom data");
+        lmp->error->all(FLERR,"lammps_scatter(): fix {} does not return peratom data", fixid);
         return;
       }
       if ((count > 1) && (fix->size_peratom_cols != count)) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter: count != values peratom for fix");
+        lmp->error->all(FLERR,"lammps_scatter(): count != values peratom for fix {}", fixid);
         return;
       }
 
@@ -4788,22 +4810,19 @@ void lammps_scatter(void *handle, const char *name, int type, int count,
     // compute
 
     if (vptr==nullptr && utils::strmatch(name,"^c_")) {
-
-      auto compute = lmp->modify->get_compute_by_id(&name[2]);
+      const char *compid = name + 2;
+      auto compute = lmp->modify->get_compute_by_id(compid);
       if (!compute) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter: unknown compute id");
+        lmp->error->all(FLERR,"lammps_scatter(): unknown compute id {}",compid);
         return;
       }
 
       if (compute->peratom_flag == 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter: compute does not return peratom data");
+        lmp->error->all(FLERR,"lammps_scatter(): compute {} does not return peratom data", compid);
         return;
       }
       if ((count > 1) && (compute->size_peratom_cols != count)) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter: count != values peratom for compute");
+        lmp->error->all(FLERR,"lammps_scatter(): count != values peratom for compute {}", compid);
         return;
       }
 
@@ -4819,28 +4838,26 @@ void lammps_scatter(void *handle, const char *name, int type, int count,
     if ((vptr == nullptr) && utils::strmatch(name,"^[id]2?_")) {
 
       int idx,icol;
-      if (utils::strmatch(name,"^[id]_")) idx = lmp->atom->find_custom(&name[2],ltype,icol);
-      else idx = lmp->atom->find_custom(&name[3],ltype,icol);
+      const char *propid;
+      if (utils::strmatch(name,"^[id]_")) propid = name + 2;
+      else propid = name + 3;
+      idx = lmp->atom->find_custom(propid,ltype,icol);
 
       if (idx < 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter: unknown property/atom id");
+        lmp->error->all(FLERR,"lammps_scatter(): unknown property/atom id {}", propid);
         return;
       }
 
       if (ltype != type) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter: mismatch property/atom type");
+        lmp->error->all(FLERR,"lammps_scatter(): mismatch property/atom {} type", propid);
         return;
       }
       if (count == 1 && icol != 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter: mismatch property/atom count");
+        lmp->error->all(FLERR,"lammps_scatter(): mismatch property/atom {} count", propid);
         return;
       }
       if (count > 1 && icol != count) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter: mismatch property/atom count");
+        lmp->error->all(FLERR,"lammps_scatter(): mismatch property/atom {} count", propid);
         return;
       }
 
@@ -4856,8 +4873,7 @@ void lammps_scatter(void *handle, const char *name, int type, int count,
     // no match
 
     if (vptr == nullptr) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"lammps_scatter: unknown property name");
+      lmp->error->all(FLERR,"lammps_scatter(): unknown property {}", name);
       return;
     }
 
@@ -4938,7 +4954,14 @@ to be the array {x[1][0], x[1][1], x[1][2], x[100][0], x[100][1], x[100][2],
 x[57][0], x[57][1], x[57][2]}, then *count* = 3, *ndata* = 3, and *ids* would
 be {1, 100, 57}.
 
-This function is not compatible with ``-DLAMMPS_BIGBIG``.
+.. admonition:: Restrictions
+    :class: warning
+
+    This function is not compatible with ``-DLAMMPS_BIGBIG``.
+
+    Atom IDs must be defined and an :doc:`atom map must be enabled <atom_modify>`
+
+    The total number of atoms must not be more than 2147483647 (max 32-bit signed int).
 
 \endverbatim
  *
@@ -4979,8 +5002,8 @@ This function is not compatible with ``-DLAMMPS_BIGBIG``.
      loop over Ndata, if I own atom ID, set its values from data
 ------------------------------------------------------------------------- */
 
-void lammps_scatter_subset(void *handle, const char *name,int type, int count,
-                                 int ndata, int *ids, void *data)
+void lammps_scatter_subset(void *handle, const char *name,int type, int count, int ndata,
+                           int *ids, void *data)
 {
   auto lmp = (LAMMPS *) handle;
 
@@ -5001,8 +5024,7 @@ void lammps_scatter_subset(void *handle, const char *name,int type, int count,
     if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
     if (lmp->atom->map_style == Atom::MAP_NONE) flag = 1;
     if (flag) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"Library error in lammps_scatter_atoms_subset");
+      lmp->error->all(FLERR,"lammps_scatter_subset(): atom-IDs must be defined and mapped");
       return;
     }
 
@@ -5011,22 +5033,19 @@ void lammps_scatter_subset(void *handle, const char *name,int type, int count,
     // fix
 
     if (vptr==nullptr && utils::strmatch(name,"^f_")) {
-
-      auto fix = lmp->modify->get_fix_by_id(&name[2]);
+      const char *fixid = name + 2;
+      auto fix = lmp->modify->get_fix_by_id(fixid);
       if (!fix) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter_subset: unknown fix id");
+        lmp->error->all(FLERR,"lammps_scatter_subset(): unknown fix id {}", fixid);
         return;
       }
 
       if (fix->peratom_flag == 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter_subset: fix does not return peratom data");
+        lmp->error->all(FLERR,"lammps_scatter_subset(): fix {} does not return peratom data", fixid);
         return;
       }
       if ((count > 1) && (fix->size_peratom_cols != count)) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter_subset: count != values peratom for fix");
+        lmp->error->all(FLERR,"lammps_scatter_subset(): count != values peratom for fix {}", fixid);
         return;
       }
 
@@ -5037,22 +5056,21 @@ void lammps_scatter_subset(void *handle, const char *name,int type, int count,
     // compute
 
     if (vptr==nullptr && utils::strmatch(name,"^c_")) {
-
-      auto compute = lmp->modify->get_compute_by_id(&name[2]);
+      const char *compid = name + 2;
+      auto compute = lmp->modify->get_compute_by_id(compid);
       if (!compute) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter_subset: unknown compute id");
+        lmp->error->all(FLERR,"lammps_scatter_subset(): unknown compute id {}", compid);
         return;
       }
 
       if (compute->peratom_flag == 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter_subset: compute does not return peratom data");
+        lmp->error->all(FLERR,"lammps_scatter_subset(): compute {} does not return peratom data",
+                        compid);
         return;
       }
       if ((count > 1) && (compute->size_peratom_cols != count)) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter_subset: count != values peratom for compute");
+        lmp->error->all(FLERR,"lammps_scatter_subset(): count != values peratom for compute {}",
+                        compid);
         return;
       }
 
@@ -5068,28 +5086,26 @@ void lammps_scatter_subset(void *handle, const char *name,int type, int count,
     if ((vptr == nullptr) && utils::strmatch(name,"^[id]2?_")) {
 
       int idx,icol;
-      if (utils::strmatch(name,"^[id]_")) idx = lmp->atom->find_custom(&name[2],ltype,icol);
-      else idx = lmp->atom->find_custom(&name[3],ltype,icol);
+      const char *propid;
+      if (utils::strmatch(name,"^[id]_")) propid = name + 2;
+      else propid = name + 3;
+      idx = lmp->atom->find_custom(propid,ltype,icol);
 
       if (idx < 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter_subset: unknown property/atom id");
+        lmp->error->all(FLERR,"lammps_scatter_subset(): unknown property/atom id {}", propid);
         return;
       }
 
       if (ltype != type) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_scatter_subset: mismatch property/atom type");
+        lmp->error->all(FLERR,"lammps_scatter_subset(): mismatch property/atom {} type", propid);
         return;
       }
       if (count == 1 && icol != 0) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather: mismatch property/atom count");
+        lmp->error->all(FLERR,"lammps_scatter_subset(): mismatch property/atom {} count", propid);
         return;
       }
       if (count > 1 && icol != count) {
-        if (lmp->comm->me == 0)
-          lmp->error->warning(FLERR,"lammps_gather: mismatch property/atom count");
+        lmp->error->all(FLERR,"lammps_scatter_subset(): mismatch property/atom {} count", propid);
         return;
       }
 
@@ -5105,9 +5121,7 @@ void lammps_scatter_subset(void *handle, const char *name,int type, int count,
     // no match
 
     if (vptr == nullptr) {
-      if (lmp->comm->me == 0)
-        lmp->error->warning(FLERR,"lammps_scatter_atoms_subset: "
-                                  "unknown property name");
+      lmp->error->all(FLERR,"lammps_scatter_subset(): unknown property {}", name);
       return;
     }
 
@@ -5246,7 +5260,7 @@ int lammps_create_atoms(void *handle, int n, const tagint *id, const int *type,
     // error if box does not exist or tags not defined
 
     int flag = 0;
-    std::string msg("Failure in lammps_create_atoms: ");
+    std::string msg("Failure in lammps_create_atoms(): ");
     if (lmp->domain->box_exist == 0) {
       flag = 1;
       msg += "trying to create atoms before before simulation box is defined";
@@ -5257,7 +5271,7 @@ int lammps_create_atoms(void *handle, int n, const tagint *id, const int *type,
     }
 
     if (flag) {
-      if (lmp->comm->me == 0) lmp->error->warning(FLERR,msg);
+      lmp->error->all(FLERR, msg);
       return -1;
     }
 

From 7f68620fa97342434b1a801e6d4b5856742318d8 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Fri, 12 Apr 2024 23:46:35 -0400
Subject: [PATCH 04/17] add example with more compact initialization for
 scatter

---
 doc/src/Python_scatter.rst | 17 +++++++++++++++--
 1 file changed, 15 insertions(+), 2 deletions(-)

diff --git a/doc/src/Python_scatter.rst b/doc/src/Python_scatter.rst
index 002a36cc98..e4444fabe7 100644
--- a/doc/src/Python_scatter.rst
+++ b/doc/src/Python_scatter.rst
@@ -54,8 +54,21 @@ like this:
    x[3] = x coord of atom with ID 2
    ...
    x[n3-1] = z coord of atom with ID natoms
-   lmp.scatter_atoms("x",1,3,x)
+   lmp.scatter_atoms("x", 1, 3, x)
+
+The coordinates can also be provided as arguments to the initializer of x:
+
+.. code-block:: python
+
+   from ctypes import c_double
+   natoms = 2
+   n3 = 3*natoms
+   # init in constructor
+   x = (n3*c_double)(0.0, 0.0, 0.0, 1.0, 1.0, 1.0)
+   lmp.scatter_atoms("x", 1, 3, x)
+   # or using a list
+   coords = [1.0, 2.0, 3.0, -3.0, -2.0, -1.0]
+   x = (c_double*len(coords))(*coords)
 
 Alternatively, you can just change values in the vector returned by
 the gather methods, since they are also ctypes vectors.
-

From 18d45d1ff0bbf5491d62db627cca0c549837227c Mon Sep 17 00:00:00 2001
From: yuhldr <yuhldr@qq.com>
Date: Sun, 14 Apr 2024 17:02:23 +0800
Subject: [PATCH 05/17] pylammps: fix get atom.mass by  atom.type

---
 python/lammps/pylammps.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/lammps/pylammps.py b/python/lammps/pylammps.py
index c693b7d565..05d6d577b9 100644
--- a/python/lammps/pylammps.py
+++ b/python/lammps/pylammps.py
@@ -196,7 +196,7 @@ class Atom(object):
 
     :type: float
     """
-    return self.get("mass", self.index)
+    return self.get("mass", self.type)
 
   @property
   def radius(self):

From 962219a446b3e78c636a694f40f2938bbdf076ff Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Sun, 14 Apr 2024 18:20:19 -0400
Subject: [PATCH 06/17] make PyLammps mass property compatible with per-atom
 masses.

---
 python/lammps/pylammps.py | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/python/lammps/pylammps.py b/python/lammps/pylammps.py
index 05d6d577b9..b959cf7d65 100644
--- a/python/lammps/pylammps.py
+++ b/python/lammps/pylammps.py
@@ -192,11 +192,23 @@ class Atom(object):
   @property
   def mass(self):
     """
-    Return the atom mass
+    Return the atom mass as a per-atom property.
+    This returns either the per-type mass or the per-atom
+    mass (AKA 'rmass') depending on what is available with
+    preference for the per-atom mass.
+
+    .. versionchanged:: TBD
+
+       Support both per-type and per-atom masses. With
+       per-type return "mass[type[i]]" else return "rmass[i]".
+       Per-atom mass is preferred if available.
 
     :type: float
     """
-    return self.get("mass", self.type)
+    if self._pylmp.lmp.extract_setting('rmass_flag'):
+      return self.get("rmass", self.index)
+    else:
+      return self.get("mass", self.type)
 
   @property
   def radius(self):

From a9b9f7f2c7835e97d4c75b8f3706878aad94e2c2 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Mon, 15 Apr 2024 01:45:05 -0400
Subject: [PATCH 07/17] correct fix ttm/mod example input

---
 examples/ttm/Si.ttm_mod                       |  4 +-
 ...tm.mod.g++.1 => log.15Apr24.ttm.mod.g++.1} | 43 ++++++++++---------
 ...tm.mod.g++.4 => log.15Apr24.ttm.mod.g++.4} | 42 +++++++++---------
 3 files changed, 45 insertions(+), 44 deletions(-)
 rename examples/ttm/{log.18May23.ttm.mod.g++.1 => log.15Apr24.ttm.mod.g++.1} (66%)
 rename examples/ttm/{log.18May23.ttm.mod.g++.4 => log.15Apr24.ttm.mod.g++.4} (65%)

diff --git a/examples/ttm/Si.ttm_mod b/examples/ttm/Si.ttm_mod
index 6e047857ff..ad407352e7 100644
--- a/examples/ttm/Si.ttm_mod
+++ b/examples/ttm/Si.ttm_mod
@@ -17,9 +17,9 @@ rho_e, electrons/volume units
 D_e, length^2/time units
 20000
 gamma_p, mass/time units
-39.235
-gamma_s, mass/time units
 24.443
+gamma_s, mass/time units
+39.235
 v_0, length/time units
 79.76
 I_0, energy/(time*length^2) units
diff --git a/examples/ttm/log.18May23.ttm.mod.g++.1 b/examples/ttm/log.15Apr24.ttm.mod.g++.1
similarity index 66%
rename from examples/ttm/log.18May23.ttm.mod.g++.1
rename to examples/ttm/log.15Apr24.ttm.mod.g++.1
index b97e8ab0ea..4b4f88b2a2 100644
--- a/examples/ttm/log.18May23.ttm.mod.g++.1
+++ b/examples/ttm/log.15Apr24.ttm.mod.g++.1
@@ -1,4 +1,5 @@
-LAMMPS (28 Mar 2023 - Development)
+LAMMPS (7 Feb 2024 - Development - patch_7Feb2024_update1-394-g03ab36a37d)
+  using 1 OpenMP thread(s) per MPI task
 units           metal
 atom_style      atomic
 boundary        p p p
@@ -13,7 +14,7 @@ mass 1 28.0855
 create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 basis 4 1 basis 5 1 basis 6 1 basis 7 1 basis 8 1
 Created 8000 atoms
   using lattice units in orthogonal box = (0 0 0) to (54.309 54.309 54.309)
-  create_atoms CPU = 0.002 seconds
+  create_atoms CPU = 0.001 seconds
 
 pair_style      sw
 pair_coeff      * * Si.sw Si
@@ -79,30 +80,30 @@ Neighbor list info ...
 Per MPI rank memory allocation (min/avg/max) = 4.433 | 4.433 | 4.433 Mbytes
    Step          Temp          TotEng      f_twotemp[1]   f_twotemp[2] 
          0                     0   -34692.79996100604   -52.79390940511979                    0
-       100     2.004897156140836   -34690.27961013186    -55.3499730543189  0.01301140393178352
-       200     2.837118035232607   -34687.74741132015   -57.93445748841876  0.02696025968760173
-       300     4.263087164947482   -34684.98084093686   -60.75945453846793  0.02175636603841567
-       400     5.568003854939066   -34682.25271040963   -63.56896518300499  0.03000618483472749
-       500     6.225602451570786   -34679.49948952029   -66.40897551884574  0.02768827702656703
-       600     7.608847536264781   -34676.69728436362   -69.32060611557282  0.05579466731854091
-       700     9.049471241531297   -34674.00093206036   -72.10055094219462 0.004335980559879032
-       800     9.826796099683211   -34671.27720242751   -74.95010610862134  0.02371649678091515
-       900      11.8609224958918   -34668.35091308811   -77.98544170794545 0.004658649791374908
-      1000     13.88037467640968   -34665.35025858006   -81.16445160194111  0.07684078334464743
-Loop time of 2.48942 on 1 procs for 1000 steps with 8000 atoms
+       100     1.255921182965094   -34691.22889627319   -54.38067722556279 0.004868249873095404
+       200     1.858362347834853    -34689.5405389424   -56.09419523244324  0.01649190747838086
+       300     2.581575104085017    -34687.9650112138   -57.69350558275053  0.01683584513983131
+       400      3.47533128765632    -34686.2796683925   -59.40465113478642 0.005727647825729662
+       500     4.080137293185865   -34684.25857873315   -61.46449138661911 0.005828121949923951
+       600     4.816104423494803   -34682.51412688349   -63.25804498666959  0.02397283419020746
+       700     5.937291156573137   -34680.64941595491   -65.17152689673857  0.02604017750117964
+       800     6.487028971399661   -34678.87151939966   -66.99420300650799 0.009720189851817886
+       900     7.461479797687167   -34677.29259652842   -68.63442522233655  0.02576822683306545
+      1000     8.696444335455215   -34675.39247806347   -70.59264558122587   0.0147252863003017
+Loop time of 5.11497 on 1 procs for 1000 steps with 8000 atoms
 
-Performance: 3.471 ns/day, 6.915 hours/ns, 401.700 timesteps/s, 3.214 Matom-step/s
-100.0% CPU use with 1 MPI tasks x no OpenMP threads
+Performance: 1.689 ns/day, 14.208 hours/ns, 195.505 timesteps/s, 1.564 Matom-step/s
+99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
 
 MPI task timing breakdown:
 Section |  min time  |  avg time  |  max time  |%varavg| %total
 ---------------------------------------------------------------
-Pair    | 2.126      | 2.126      | 2.126      |   0.0 | 85.40
+Pair    | 4.3498     | 4.3498     | 4.3498     |   0.0 | 85.04
 Neigh   | 0          | 0          | 0          |   0.0 |  0.00
-Comm    | 0.016147   | 0.016147   | 0.016147   |   0.0 |  0.65
-Output  | 0.0013116  | 0.0013116  | 0.0013116  |   0.0 |  0.05
-Modify  | 0.33864    | 0.33864    | 0.33864    |   0.0 | 13.60
-Other   |            | 0.007318   |            |       |  0.29
+Comm    | 0.037981   | 0.037981   | 0.037981   |   0.0 |  0.74
+Output  | 0.0025641  | 0.0025641  | 0.0025641  |   0.0 |  0.05
+Modify  | 0.71279    | 0.71279    | 0.71279    |   0.0 | 13.94
+Other   |            | 0.01179    |            |       |  0.23
 
 Nlocal:           8000 ave        8000 max        8000 min
 Histogram: 1 0 0 0 0 0 0 0 0 0
@@ -117,4 +118,4 @@ Total # of neighbors = 272000
 Ave neighs/atom = 34
 Neighbor list builds = 0
 Dangerous builds = 0
-Total wall time: 0:00:02
+Total wall time: 0:00:05
diff --git a/examples/ttm/log.18May23.ttm.mod.g++.4 b/examples/ttm/log.15Apr24.ttm.mod.g++.4
similarity index 65%
rename from examples/ttm/log.18May23.ttm.mod.g++.4
rename to examples/ttm/log.15Apr24.ttm.mod.g++.4
index ea675c8594..231a9af2c9 100644
--- a/examples/ttm/log.18May23.ttm.mod.g++.4
+++ b/examples/ttm/log.15Apr24.ttm.mod.g++.4
@@ -1,5 +1,5 @@
-LAMMPS (28 Mar 2023 - Development)
-WARNING: Using I/O redirection is unreliable with parallel runs. Better to use the -in switch to read input files. (../lammps.cpp:531)
+LAMMPS (7 Feb 2024 - Development - patch_7Feb2024_update1-394-g03ab36a37d)
+  using 1 OpenMP thread(s) per MPI task
 units           metal
 atom_style      atomic
 boundary        p p p
@@ -14,7 +14,7 @@ mass 1 28.0855
 create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 basis 4 1 basis 5 1 basis 6 1 basis 7 1 basis 8 1
 Created 8000 atoms
   using lattice units in orthogonal box = (0 0 0) to (54.309 54.309 54.309)
-  create_atoms CPU = 0.001 seconds
+  create_atoms CPU = 0.002 seconds
 
 pair_style      sw
 pair_coeff      * * Si.sw Si
@@ -80,30 +80,30 @@ Neighbor list info ...
 Per MPI rank memory allocation (min/avg/max) = 3.436 | 3.436 | 3.436 Mbytes
    Step          Temp          TotEng      f_twotemp[1]   f_twotemp[2] 
          0                     0   -34692.79996100361   -52.79390940511979                    0
-       100     1.852689977101411   -34690.49204900486   -55.14271612882064  0.02726188676577098
-       200     2.735750477179192   -34688.11139028054   -57.57110998717798  0.03387986355513584
-       300     3.931848271449558   -34685.54667417785   -60.18684521127231  0.02261256315262403
-       400     5.462009198576365   -34682.74455105668   -63.05420336037233 0.002402241637719578
-       500     6.267811692893873   -34679.96493887379   -65.93304222280051  0.02448378880218699
-       600      7.21148216150661   -34677.41455784726   -68.58391420045926  0.04114045759945374
-       700      8.84660534187052   -34674.40610468235   -71.68798344296859  0.02372984027434538
-       800      10.1748456457686   -34671.08749605772   -75.11943618276216 0.007538225788030307
-       900     11.27479036162859    -34668.4118066423   -77.92921692176756  0.02537529314475071
-      1000     13.26881394868076   -34665.56617589539   -80.91544540266317  0.03112665440209921
-Loop time of 0.995347 on 4 procs for 1000 steps with 8000 atoms
+       100      1.20337355884597   -34691.30677367127   -54.30747356568817  0.01557346850238741
+       200     1.709631732825883   -34689.83859944795    -55.7982356998371  0.02508386983502213
+       300     2.488524478071323   -34688.26307995134    -57.3977272154369  0.02664346353990833
+       400      3.38535890366476   -34686.51395648598   -59.17547816947624  0.02164200191836632
+       500     3.838163353802383   -34684.79466673204   -60.92228950760077 0.005860499116196545
+       600     4.675913079756001   -34683.03448988724   -62.72423959707044   0.0106700119158327
+       700     5.637185532827328   -34681.25888274477    -64.5491928842093  0.01568536325219336
+       800     6.316986413957468   -34679.29231578312   -66.57005328290739  0.02035373879569394
+       900     7.211479047111087   -34677.61236020172   -68.30976417874265  0.03239086895076279
+      1000     8.431725106300505   -34675.81097854214     -70.161139196977  0.01219385884660358
+Loop time of 1.73439 on 4 procs for 1000 steps with 8000 atoms
 
-Performance: 8.680 ns/day, 2.765 hours/ns, 1004.675 timesteps/s, 8.037 Matom-step/s
-97.9% CPU use with 4 MPI tasks x no OpenMP threads
+Performance: 4.982 ns/day, 4.818 hours/ns, 576.572 timesteps/s, 4.613 Matom-step/s
+99.7% CPU use with 4 MPI tasks x 1 OpenMP threads
 
 MPI task timing breakdown:
 Section |  min time  |  avg time  |  max time  |%varavg| %total
 ---------------------------------------------------------------
-Pair    | 0.65351    | 0.6616     | 0.66783    |   0.8 | 66.47
+Pair    | 1.127      | 1.1392     | 1.1511     |   1.0 | 65.68
 Neigh   | 0          | 0          | 0          |   0.0 |  0.00
-Comm    | 0.041606   | 0.048314   | 0.056589   |   2.9 |  4.85
-Output  | 0.0014609  | 0.0014742  | 0.0014968  |   0.0 |  0.15
-Modify  | 0.27934    | 0.28016    | 0.28089    |   0.1 | 28.15
-Other   |            | 0.003798   |            |       |  0.38
+Comm    | 0.068488   | 0.082304   | 0.094797   |   4.1 |  4.75
+Output  | 0.0024745  | 0.0025221  | 0.0025705  |   0.1 |  0.15
+Modify  | 0.50194    | 0.50329    | 0.50522    |   0.2 | 29.02
+Other   |            | 0.007117   |            |       |  0.41
 
 Nlocal:           2000 ave        2000 max        2000 min
 Histogram: 4 0 0 0 0 0 0 0 0 0

From 3a94e4df2d2e981a1e4fd957947d16d8d5a85075 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Mon, 15 Apr 2024 02:50:16 -0400
Subject: [PATCH 08/17] print single warning when some rho[i] exceeded rhomax
 of the current EAM potential

---
 doc/src/pair_eam.rst         | 14 ++++++++++++++
 src/INTEL/pair_eam_intel.cpp | 18 +++++++++++++++---
 src/MANYBODY/pair_eam.cpp    | 18 +++++++++++++++++-
 src/MANYBODY/pair_eam.h      |  3 ++-
 src/OPENMP/pair_eam_omp.cpp  | 17 ++++++++++++++++-
 src/OPT/pair_eam_opt.cpp     | 17 ++++++++++++++++-
 6 files changed, 80 insertions(+), 7 deletions(-)

diff --git a/doc/src/pair_eam.rst b/doc/src/pair_eam.rst
index 654f49c166..bb6a01eabf 100644
--- a/doc/src/pair_eam.rst
+++ b/doc/src/pair_eam.rst
@@ -140,6 +140,20 @@ The OpenKIM Project at
 provides EAM potentials that can be used directly in LAMMPS with the
 :doc:`kim command <kim_commands>` interface.
 
+.. warning::
+
+   The EAM potential files tabulate the embedding energy as a function
+   of the local electron density :math:`\rho`.  When atoms get too
+   close, this electron density may exceed the range for which the
+   embedding energy was tabulated for.  For simplicity and to avoid
+   errors during equilibration of randomized geometries, LAMMPS will
+   assume a linearly increasing embedding energy for electron densities
+   beyond the maximum tabulated value.  This usually means that the EAM
+   model is not a good model for the kind of system under investigation.
+   LAMMPS will print a single warning when this happens.  It may be
+   harmless at the beginning of an equilibration but would be a big
+   concern for accuracy if it happens during production runs.
+
 ----------
 
 For style *eam*, potential values are read from a file that is in the
diff --git a/src/INTEL/pair_eam_intel.cpp b/src/INTEL/pair_eam_intel.cpp
index bd78c3239d..4fd527858e 100644
--- a/src/INTEL/pair_eam_intel.cpp
+++ b/src/INTEL/pair_eam_intel.cpp
@@ -234,7 +234,7 @@ void PairEAMIntel::eval(const int offload, const int vflag,
   const int istride = fc.rhor_istride();
   const int jstride = fc.rhor_jstride();
   const int fstride = fc.frho_stride();
-
+  int beyond_rhomax = 0;
   {
     #if defined(__MIC__) && defined(_LMP_INTEL_OFFLOAD)
     *timer_compute = MIC_Wtime();
@@ -453,7 +453,10 @@ void PairEAMIntel::eval(const int offload, const int vflag,
         if (EFLAG) {
           flt_t phi = ((frho_spline_e[ioff].a*p + frho_spline_e[ioff].b)*p +
                        frho_spline_e[ioff].c)*p + frho_spline_e[ioff].d;
-          if (rho[i] > frhomax) phi += fp_f[i] * (rho[i]-frhomax);
+          if (rho[i] > frhomax) {
+            phi += fp_f[i] * (rho[i]-frhomax);
+            beyond_rhomax = 1;
+          }
           if (!ONETYPE) {
             const int ptr_off=itype*ntypes + itype;
             oscale = scale_f[ptr_off];
@@ -653,6 +656,16 @@ void PairEAMIntel::eval(const int offload, const int vflag,
   else
     fix->stop_watch(TIME_HOST_PAIR);
 
+  if (EFLAG && (exceeded_rhomax >= 0)) {
+    MPI_Allreduce(&beyond_rhomax, &exceeded_rhomax, 1, MPI_INT, MPI_MAX, world);
+    if (exceeded_rhomax > 0) {
+      if (comm->me == 0)
+        error->warning(FLERR, "Local rho[i] exceeded rhomax of EAM potential table. "
+                       "Computed embedding term is unreliable.");
+      exceeded_rhomax = -1;
+    }
+  }
+
   if (EFLAG || vflag)
     fix->add_result_array(f_start, ev_global, offload, eatom, 0, vflag);
   else
@@ -849,4 +862,3 @@ void PairEAMIntel::unpack_forward_comm(int n, int first, double *buf,
   last = first + n;
   for (i = first; i < last; i++) fp_f[i] = buf[m++];
 }
-
diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp
index 669a5cadbb..78690d38bf 100644
--- a/src/MANYBODY/pair_eam.cpp
+++ b/src/MANYBODY/pair_eam.cpp
@@ -43,6 +43,7 @@ PairEAM::PairEAM(LAMMPS *lmp) : Pair(lmp)
   unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
 
   nmax = 0;
+  exceeded_rhomax = 0;
   rho = nullptr;
   fp = nullptr;
   numforce = nullptr;
@@ -145,6 +146,7 @@ void PairEAM::compute(int eflag, int vflag)
   double *coeff;
   int *ilist,*jlist,*numneigh,**firstneigh;
 
+  int beyond_rhomax = 0;
   evdwl = 0.0;
   ev_init(eflag,vflag);
 
@@ -237,7 +239,9 @@ void PairEAM::compute(int eflag, int vflag)
     fp[i] = (coeff[0]*p + coeff[1])*p + coeff[2];
     if (eflag) {
       phi = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
-      if (rho[i] > rhomax) phi += fp[i] * (rho[i]-rhomax);
+      if (rho[i] > rhomax) {
+        phi += fp[i] * (rho[i]-rhomax);
+      }
       phi *= scale[type[i]][type[i]];
       if (eflag_global) eng_vdwl += phi;
       if (eflag_atom) eatom[i] += phi;
@@ -322,6 +326,16 @@ void PairEAM::compute(int eflag, int vflag)
     }
   }
 
+  if (eflag && (exceeded_rhomax >= 0)) {
+    MPI_Allreduce(&beyond_rhomax, &exceeded_rhomax, 1, MPI_INT, MPI_MAX, world);
+    if (exceeded_rhomax > 0) {
+      if (comm->me == 0)
+        error->warning(FLERR, "Local rho[i] exceeded rhomax of EAM potential table. "
+                       "Computed embedding term is unreliable.");
+      exceeded_rhomax = -1;
+    }
+  }
+
   if (vflag_fdotr) virial_fdotr_compute();
 }
 
@@ -424,6 +438,8 @@ void PairEAM::init_style()
 
   neighbor->add_request(this);
   embedstep = -1;
+
+  exceeded_rhomax = 0;
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/MANYBODY/pair_eam.h b/src/MANYBODY/pair_eam.h
index 67bbde570d..7ecc2ce0df 100644
--- a/src/MANYBODY/pair_eam.h
+++ b/src/MANYBODY/pair_eam.h
@@ -63,7 +63,8 @@ class PairEAM : public Pair {
   void swap_eam(double *, double **) override;
 
  protected:
-  int nmax;    // allocated size of per-atom arrays
+  int nmax;               // allocated size of per-atom arrays
+  int exceeded_rhomax;    // 0 or 1 when rho[i] exceeded rhomax, -1 when not to check
   double cutforcesq;
   double **scale;
   bigint embedstep;    // timestep, the embedding term was computed
diff --git a/src/OPENMP/pair_eam_omp.cpp b/src/OPENMP/pair_eam_omp.cpp
index e99fbedbb7..57aff9ea3a 100644
--- a/src/OPENMP/pair_eam_omp.cpp
+++ b/src/OPENMP/pair_eam_omp.cpp
@@ -17,6 +17,7 @@
 
 #include "atom.h"
 #include "comm.h"
+#include "error.h"
 #include "force.h"
 #include "memory.h"
 #include "neigh_list.h"
@@ -102,6 +103,7 @@ void PairEAMOMP::eval(int iifrom, int iito, ThrData * const thr)
   double *coeff;
   int *ilist,*jlist,*numneigh,**firstneigh;
 
+  int beyond_rhomax = 0;
   evdwl = 0.0;
 
   const auto * _noalias const x = (dbl3_t *) atom->x[0];
@@ -203,7 +205,10 @@ void PairEAMOMP::eval(int iifrom, int iito, ThrData * const thr)
     fp[i] = (coeff[0]*p + coeff[1])*p + coeff[2];
     if (EFLAG) {
       phi = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
-      if (rho[i] > rhomax) phi += fp[i] * (rho[i]-rhomax);
+      if (rho[i] > rhomax) {
+        phi += fp[i] * (rho[i]-rhomax);
+        beyond_rhomax = 1;
+      }
       e_tally_thr(this, i, i, nlocal, NEWTON_PAIR, scale[type[i]][type[i]]*phi, 0.0, thr);
     }
   }
@@ -298,6 +303,16 @@ void PairEAMOMP::eval(int iifrom, int iito, ThrData * const thr)
     f[i].y += fytmp;
     f[i].z += fztmp;
   }
+
+  if (EFLAG && (exceeded_rhomax >= 0)) {
+    MPI_Allreduce(&beyond_rhomax, &exceeded_rhomax, 1, MPI_INT, MPI_MAX, world);
+    if (exceeded_rhomax > 0) {
+      if (comm->me == 0)
+        error->warning(FLERR, "Local rho[i] exceeded rhomax of EAM potential table. "
+                       "Computed embedding term is unreliable.");
+      exceeded_rhomax = -1;
+    }
+  }
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/OPT/pair_eam_opt.cpp b/src/OPT/pair_eam_opt.cpp
index 0560b0693a..6b4f1063e2 100644
--- a/src/OPT/pair_eam_opt.cpp
+++ b/src/OPT/pair_eam_opt.cpp
@@ -23,6 +23,7 @@
 
 #include "atom.h"
 #include "comm.h"
+#include "error.h"
 #include "force.h"
 #include "memory.h"
 #include "neigh_list.h"
@@ -118,6 +119,7 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairEAMOpt::eval()
 
   int ntypes = atom->ntypes;
   int ntypes2 = ntypes * ntypes;
+  int beyond_rhomax = 0;
 
   auto *_noalias fast_alpha =
       (fast_alpha_t *) malloc((size_t) ntypes2 * (nr + 1) * sizeof(fast_alpha_t));
@@ -251,7 +253,10 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairEAMOpt::eval()
     fp[i] = (coeff[0] * p + coeff[1]) * p + coeff[2];
     if (EFLAG) {
       double phi = ((coeff[3] * p + coeff[4]) * p + coeff[5]) * p + coeff[6];
-      if (rho[i] > rhomax) phi += fp[i] * (rho[i] - rhomax);
+      if (rho[i] > rhomax) {
+        phi += fp[i] * (rho[i] - rhomax);
+        beyond_rhomax = 1;
+      }
       phi *= scale[type[i]][type[i]];
       if (eflag_global) eng_vdwl += phi;
       if (eflag_atom) eatom[i] += phi;
@@ -361,5 +366,15 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairEAMOpt::eval()
   free(fast_gamma);
   fast_gamma = nullptr;
 
+  if (EFLAG && (exceeded_rhomax >= 0)) {
+    MPI_Allreduce(&beyond_rhomax, &exceeded_rhomax, 1, MPI_INT, MPI_MAX, world);
+    if (exceeded_rhomax > 0) {
+      if (comm->me == 0)
+        error->warning(FLERR, "Local rho[i] exceeded rhomax of EAM potential table. "
+                       "Computed embedding term is unreliable.");
+      exceeded_rhomax = -1;
+    }
+  }
+
   if (vflag_fdotr) virial_fdotr_compute();
 }

From b2fded39332a1f009dfdac487a06ecf6fe6de19c Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Mon, 15 Apr 2024 11:24:54 -0400
Subject: [PATCH 09/17] restore missing line

---
 src/MANYBODY/pair_eam.cpp | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp
index 78690d38bf..fdbedf6804 100644
--- a/src/MANYBODY/pair_eam.cpp
+++ b/src/MANYBODY/pair_eam.cpp
@@ -241,6 +241,7 @@ void PairEAM::compute(int eflag, int vflag)
       phi = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
       if (rho[i] > rhomax) {
         phi += fp[i] * (rho[i]-rhomax);
+        beyond_rhomax = 1;
       }
       phi *= scale[type[i]][type[i]];
       if (eflag_global) eng_vdwl += phi;

From 16e5fce31e2dca5fc9a0a63f199206a36a421b9d Mon Sep 17 00:00:00 2001
From: Steve Plimpton <sjplimp@gmail.com>
Date: Mon, 15 Apr 2024 17:24:47 -0600
Subject: [PATCH 10/17] adjust error messaging for rho[i] > rhomax in EAM

---
 src/MANYBODY/pair_eam.cpp | 14 +++++++-------
 src/MANYBODY/pair_eam.h   |  6 ++++--
 2 files changed, 11 insertions(+), 9 deletions(-)

diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp
index fdbedf6804..8ea4efe868 100644
--- a/src/MANYBODY/pair_eam.cpp
+++ b/src/MANYBODY/pair_eam.cpp
@@ -43,7 +43,6 @@ PairEAM::PairEAM(LAMMPS *lmp) : Pair(lmp)
   unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
 
   nmax = 0;
-  exceeded_rhomax = 0;
   rho = nullptr;
   fp = nullptr;
   numforce = nullptr;
@@ -146,10 +145,11 @@ void PairEAM::compute(int eflag, int vflag)
   double *coeff;
   int *ilist,*jlist,*numneigh,**firstneigh;
 
-  int beyond_rhomax = 0;
   evdwl = 0.0;
   ev_init(eflag,vflag);
 
+  int beyond_rhomax = 0;
+  
   // grow energy and fp arrays if necessary
   // need to be atom->nmax in length
 
@@ -327,13 +327,13 @@ void PairEAM::compute(int eflag, int vflag)
     }
   }
 
-  if (eflag && (exceeded_rhomax >= 0)) {
+  if (eflag && (!exceeded_rhomax)) {
     MPI_Allreduce(&beyond_rhomax, &exceeded_rhomax, 1, MPI_INT, MPI_MAX, world);
-    if (exceeded_rhomax > 0) {
+    if (exceeded_rhomax) {
       if (comm->me == 0)
-        error->warning(FLERR, "Local rho[i] exceeded rhomax of EAM potential table. "
-                       "Computed embedding term is unreliable.");
-      exceeded_rhomax = -1;
+        error->warning(FLERR,
+                       "A per-atom density exceeded rhomax of EAM potential table - "
+                       "a linear extrapolation to the energy was made");
     }
   }
 
diff --git a/src/MANYBODY/pair_eam.h b/src/MANYBODY/pair_eam.h
index 7ecc2ce0df..23c417d1a5 100644
--- a/src/MANYBODY/pair_eam.h
+++ b/src/MANYBODY/pair_eam.h
@@ -64,10 +64,12 @@ class PairEAM : public Pair {
 
  protected:
   int nmax;               // allocated size of per-atom arrays
-  int exceeded_rhomax;    // 0 or 1 when rho[i] exceeded rhomax, -1 when not to check
   double cutforcesq;
   double **scale;
-  bigint embedstep;    // timestep, the embedding term was computed
+  bigint embedstep;       // timestep, the embedding term was computed
+
+  int exceeded_rhomax;    // global flag for whether rho[i] has exceeded rhomax
+                          // on a step energy is computed - 0 = no, 1 = yes
 
   // per-atom arrays
 

From 8dd1d12a73d23e213548ec1c6d027ee1d287acc4 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Mon, 15 Apr 2024 21:04:02 -0400
Subject: [PATCH 11/17] apply clang-format and whitespace cleanup

---
 src/MANYBODY/pair_eam.cpp | 2 +-
 src/MANYBODY/pair_eam.h   | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp
index 8ea4efe868..d1b56a9eae 100644
--- a/src/MANYBODY/pair_eam.cpp
+++ b/src/MANYBODY/pair_eam.cpp
@@ -149,7 +149,7 @@ void PairEAM::compute(int eflag, int vflag)
   ev_init(eflag,vflag);
 
   int beyond_rhomax = 0;
-  
+
   // grow energy and fp arrays if necessary
   // need to be atom->nmax in length
 
diff --git a/src/MANYBODY/pair_eam.h b/src/MANYBODY/pair_eam.h
index 23c417d1a5..24221a07ce 100644
--- a/src/MANYBODY/pair_eam.h
+++ b/src/MANYBODY/pair_eam.h
@@ -63,10 +63,10 @@ class PairEAM : public Pair {
   void swap_eam(double *, double **) override;
 
  protected:
-  int nmax;               // allocated size of per-atom arrays
+  int nmax;    // allocated size of per-atom arrays
   double cutforcesq;
   double **scale;
-  bigint embedstep;       // timestep, the embedding term was computed
+  bigint embedstep;    // timestep, the embedding term was computed
 
   int exceeded_rhomax;    // global flag for whether rho[i] has exceeded rhomax
                           // on a step energy is computed - 0 = no, 1 = yes

From 28620a66316b89e3cac342f32e295a576ad4bef2 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Mon, 15 Apr 2024 21:04:57 -0400
Subject: [PATCH 12/17] propagate changes to other EAM implementations

---
 src/INTEL/pair_eam_intel.cpp | 10 +++++-----
 src/OPENMP/pair_eam_omp.cpp  | 10 +++++-----
 src/OPT/pair_eam_opt.cpp     | 10 +++++-----
 3 files changed, 15 insertions(+), 15 deletions(-)

diff --git a/src/INTEL/pair_eam_intel.cpp b/src/INTEL/pair_eam_intel.cpp
index 4fd527858e..ed92a41a11 100644
--- a/src/INTEL/pair_eam_intel.cpp
+++ b/src/INTEL/pair_eam_intel.cpp
@@ -656,13 +656,13 @@ void PairEAMIntel::eval(const int offload, const int vflag,
   else
     fix->stop_watch(TIME_HOST_PAIR);
 
-  if (EFLAG && (exceeded_rhomax >= 0)) {
+  if (EFLAG && (!exceeded_rhomax)) {
     MPI_Allreduce(&beyond_rhomax, &exceeded_rhomax, 1, MPI_INT, MPI_MAX, world);
-    if (exceeded_rhomax > 0) {
+    if (exceeded_rhomax) {
       if (comm->me == 0)
-        error->warning(FLERR, "Local rho[i] exceeded rhomax of EAM potential table. "
-                       "Computed embedding term is unreliable.");
-      exceeded_rhomax = -1;
+        error->warning(FLERR,
+                       "A per-atom density exceeded rhomax of EAM potential table - "
+                       "a linear extrapolation to the energy was made");
     }
   }
 
diff --git a/src/OPENMP/pair_eam_omp.cpp b/src/OPENMP/pair_eam_omp.cpp
index 57aff9ea3a..1d3c7ded3a 100644
--- a/src/OPENMP/pair_eam_omp.cpp
+++ b/src/OPENMP/pair_eam_omp.cpp
@@ -304,13 +304,13 @@ void PairEAMOMP::eval(int iifrom, int iito, ThrData * const thr)
     f[i].z += fztmp;
   }
 
-  if (EFLAG && (exceeded_rhomax >= 0)) {
+  if (EFLAG && (!exceeded_rhomax)) {
     MPI_Allreduce(&beyond_rhomax, &exceeded_rhomax, 1, MPI_INT, MPI_MAX, world);
-    if (exceeded_rhomax > 0) {
+    if (exceeded_rhomax) {
       if (comm->me == 0)
-        error->warning(FLERR, "Local rho[i] exceeded rhomax of EAM potential table. "
-                       "Computed embedding term is unreliable.");
-      exceeded_rhomax = -1;
+        error->warning(FLERR,
+                       "A per-atom density exceeded rhomax of EAM potential table - "
+                       "a linear extrapolation to the energy was made");
     }
   }
 }
diff --git a/src/OPT/pair_eam_opt.cpp b/src/OPT/pair_eam_opt.cpp
index 6b4f1063e2..20463515da 100644
--- a/src/OPT/pair_eam_opt.cpp
+++ b/src/OPT/pair_eam_opt.cpp
@@ -366,13 +366,13 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairEAMOpt::eval()
   free(fast_gamma);
   fast_gamma = nullptr;
 
-  if (EFLAG && (exceeded_rhomax >= 0)) {
+  if (EFLAG && (!exceeded_rhomax)) {
     MPI_Allreduce(&beyond_rhomax, &exceeded_rhomax, 1, MPI_INT, MPI_MAX, world);
-    if (exceeded_rhomax > 0) {
+    if (exceeded_rhomax) {
       if (comm->me == 0)
-        error->warning(FLERR, "Local rho[i] exceeded rhomax of EAM potential table. "
-                       "Computed embedding term is unreliable.");
-      exceeded_rhomax = -1;
+        error->warning(FLERR,
+                       "A per-atom density exceeded rhomax of EAM potential table - "
+                       "a linear extrapolation to the energy was made");
     }
   }
 

From a810f1eca83df21039846e6a97b7a47e54c08b50 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jacob.r.gissinger@gmail.com>
Date: Tue, 16 Apr 2024 00:20:18 -0400
Subject: [PATCH 13/17] Update Howto_type_labels.rst

---
 doc/src/Howto_type_labels.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/src/Howto_type_labels.rst b/doc/src/Howto_type_labels.rst
index d7460d50a7..6b734f5c8d 100644
--- a/doc/src/Howto_type_labels.rst
+++ b/doc/src/Howto_type_labels.rst
@@ -130,4 +130,4 @@ modification in multiple simulations or different systems.
 
 .. _Typelabel24:
 
-**(Gissinger)** J. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M. Choi, D. Karls, A. Stukowski, W, Im, H. Heinz, A. Kohlmeyer, and E. Tadmor, J Phys Chem B, 128, 3282-3297 (2024)
+**(Gissinger)** J. R. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M. Choi, D. S. Karls, A. Stukowski, W. Im, H. Heinz, A. Kohlmeyer, and E. B. Tadmor, J Phys Chem B, 128, 3282-3297 (2024).

From 49b5361088d3c961e74b987b6975635aafd8083b Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jacob.r.gissinger@gmail.com>
Date: Tue, 16 Apr 2024 00:24:25 -0400
Subject: [PATCH 14/17] Update labelmap.rst

---
 doc/src/labelmap.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/src/labelmap.rst b/doc/src/labelmap.rst
index 5e0d1e92e9..70ca9262da 100644
--- a/doc/src/labelmap.rst
+++ b/doc/src/labelmap.rst
@@ -104,4 +104,4 @@ none
 
 .. _Typelabel:
 
-**(Gissinger)** J. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M. Choi, D. Karls, A. Stukowski, W, Im, H. Heinz, A. Kohlmeyer, and E. Tadmor, J Phys Chem B, 128, 3282-3297 (2024)
+**(Gissinger)** J. R. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M. Choi, D. S. Karls, A. Stukowski, W. Im, H. Heinz, A. Kohlmeyer, and E. B. Tadmor, J Phys Chem B, 128, 3282-3297 (2024).

From 95e43644002928b1be0fec1b501e2b0931af6d4b Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Tue, 16 Apr 2024 19:51:01 -0400
Subject: [PATCH 15/17] breathe is currently not compatible with sphinx 7.3

---
 doc/utils/requirements.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/utils/requirements.txt b/doc/utils/requirements.txt
index 5bb8e3911d..816d52bf54 100644
--- a/doc/utils/requirements.txt
+++ b/doc/utils/requirements.txt
@@ -1,4 +1,4 @@
-Sphinx >= 5.3.0, <8.0
+Sphinx >= 5.3.0, <7.3
 sphinxcontrib-spelling
 sphinxcontrib-jquery
 git+https://github.com/akohlmey/sphinx-fortran@parallel-read

From 98eb34561554f3a7e428f44e818ce0d6f3b5d043 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Tue, 16 Apr 2024 20:07:30 -0400
Subject: [PATCH 16/17] small tweak

---
 doc/src/Intro_features.rst | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/doc/src/Intro_features.rst b/doc/src/Intro_features.rst
index f7b4dd319b..98e9d981ed 100644
--- a/doc/src/Intro_features.rst
+++ b/doc/src/Intro_features.rst
@@ -29,7 +29,7 @@ General features
 * spatial decomposition of simulation domain for MPI parallelism
 * particle decomposition inside spatial decomposition for OpenMP and GPU parallelism
 * GPLv2 licensed open-source distribution
-* highly portable C++-11
+* highly portable C++-11 (optional packages may require C++17)
 * modular code with most functionality in optional packages
 * only depends on MPI library for basic parallel functionality, MPI stub for serial compilation
 * other libraries are optional and only required for specific packages
@@ -81,7 +81,7 @@ commands)
 * pairwise potentials: Lennard-Jones, Buckingham, Morse, Born-Mayer-Huggins, Yukawa, soft, Class II (COMPASS), hydrogen bond, harmonic, gaussian, tabulated, scripted
 * charged pairwise potentials: Coulombic, point-dipole
 * many-body potentials: EAM, Finnis/Sinclair, MEAM, MEAM+SW, EIM, EDIP, ADP, Stillinger-Weber, Tersoff, REBO, AIREBO, ReaxFF, COMB, Streitz-Mintmire, 3-body polymorphic, BOP, Vashishta
-* machine learning potentials: ACE, AGNI, GAP, Behler-Parrinello (N2P2), POD, RANN
+* machine learning potentials: ACE, AGNI, GAP, Behler-Parrinello (N2P2), POD, RANN, SNAP
 * interfaces to ML potentials distributed by external groups: ANI, ChIMES, DeepPot, HIPNN, MTP
 * long-range interactions for charge, point-dipoles, and LJ dispersion:  Ewald, Wolf, PPPM (similar to particle-mesh Ewald), MSM, ScaFaCoS
 * polarization models: :doc:`QEq <fix_qeq>`, :doc:`core/shell model <Howto_coreshell>`, :doc:`Drude dipole model <Howto_drude>`

From 101a937870fc0fa9844068e0db27a0ec75e3bdb7 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Tue, 16 Apr 2024 23:02:29 -0400
Subject: [PATCH 17/17] rephrase warning about EAM deansity exceeding tabulated
 embedding term

---
 doc/src/pair_eam.rst | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/doc/src/pair_eam.rst b/doc/src/pair_eam.rst
index bb6a01eabf..cdb6487981 100644
--- a/doc/src/pair_eam.rst
+++ b/doc/src/pair_eam.rst
@@ -145,14 +145,15 @@ provides EAM potentials that can be used directly in LAMMPS with the
    The EAM potential files tabulate the embedding energy as a function
    of the local electron density :math:`\rho`.  When atoms get too
    close, this electron density may exceed the range for which the
-   embedding energy was tabulated for.  For simplicity and to avoid
-   errors during equilibration of randomized geometries, LAMMPS will
+   embedding energy was tabulated for.  To avoid crashes, LAMMPS will
    assume a linearly increasing embedding energy for electron densities
-   beyond the maximum tabulated value.  This usually means that the EAM
-   model is not a good model for the kind of system under investigation.
-   LAMMPS will print a single warning when this happens.  It may be
-   harmless at the beginning of an equilibration but would be a big
-   concern for accuracy if it happens during production runs.
+   beyond the maximum tabulated value.  LAMMPS will print a warning when
+   this happens.  It may be acceptable at the beginning of an
+   equilibration (e.g. when using randomized coordinates) but would be a
+   big concern for accuracy if it happens during production runs.  The
+   EAM potential file triggering the warning during production is thus
+   not a good choice, and the EAM model in general not likely a good
+   model for the kind of system under investigation.
 
 ----------