From 57f639c0e535b781738bd394abcae45f2a0d8064 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 18 Jul 2020 12:42:47 -0600
Subject: [PATCH 01/21] bond/react:reset_mol_ids keyword

---
 src/USER-REACTION/fix_bond_react.cpp | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index 2c6c20c844..c908526067 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -144,7 +144,8 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
 
   int iarg = 3;
   stabilization_flag = 0;
-  int num_common_keywords = 1;
+  reset_mol_ids_flag = 1;
+  int num_common_keywords = 2;
   for (int m = 0; m < num_common_keywords; m++) {
     if (strcmp(arg[iarg],"stabilization") == 0) {
       if (strcmp(arg[iarg+1],"no") == 0) {
@@ -162,6 +163,14 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
         nve_limit_xmax = arg[iarg+3];
         iarg += 4;
       }
+    } else if (strcmp(arg[iarg],"reset_mol_ids") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: "
+                                    "'reset_mol_ids' keyword has too few arguments");
+      if (strcmp(arg[iarg+1],"yes") == 0) iarg += 2; //default
+      if (strcmp(arg[iarg+1],"no") == 0) {
+        reset_mol_ids_flag = 0;
+        iarg += 2;
+      }
     } else if (strcmp(arg[iarg],"react") == 0) {
       break;
     } else error->all(FLERR,"Illegal fix bond/react command: unknown keyword");
@@ -2503,11 +2512,14 @@ void FixBondReact::ghost_glovecast()
 }
 
 /* ----------------------------------------------------------------------
-update charges, types, special lists and all topology
+update molecule IDs, charges, types, special lists and all topology
 ------------------------------------------------------------------------- */
 
 void FixBondReact::update_everything()
 {
+  if (reset_mol_ids_flag)
+    input->one("reset_mol_ids " + std::string(group->names[igroup]));
+
   int *type = atom->type;
 
   int nlocal = atom->nlocal;

From 371a5c5b6156de97ba992b46304b84be2888d47b Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 18 Jul 2020 12:44:34 -0600
Subject: [PATCH 02/21] bond/react: reset_mol_ids docs

---
 doc/src/fix_bond_react.rst | 24 +++++++++++++++++-------
 1 file changed, 17 insertions(+), 7 deletions(-)

diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst
index d5f917bfdb..d1671d7c78 100644
--- a/doc/src/fix_bond_react.rst
+++ b/doc/src/fix_bond_react.rst
@@ -14,19 +14,22 @@ Syntax
      react react-ID react-group-ID Nevery Rmin Rmax template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ...
      ...
 
-* ID, group-ID are documented in :doc:`fix <fix>` command. Group-ID is ignored.
+* ID, group-ID are documented in :doc:`fix <fix>` command.
 * bond/react = style name of this fix command
 * the common keyword/values may be appended directly after 'bond/react'
 * this applies to all reaction specifications (below)
-* common_keyword = *stabilization*
+* common_keyword = *stabilization* or *reset_mol_ids*
 
   .. parsed-literal::
 
        *stabilization* values = *no* or *yes* *group-ID* *xmax*
-         *no* = no reaction site stabilization
+         *no* = no reaction site stabilization (default)
          *yes* = perform reaction site stabilization
            *group-ID* = user-assigned prefix for the dynamic group of atoms not currently involved in a reaction
            *xmax* = xmax value that is used by an internally-created :doc:`nve/limit <fix_nve_limit>` integrator
+       *reset_mol_ids* values = *yes* or *no*
+         *yes* = update molecule IDs based on new global topology (default)
+         *no* = do not update molecule IDs
 
 * react = mandatory argument indicating new reaction specification
 * react-ID = user-assigned name for the reaction
@@ -50,9 +53,9 @@ Syntax
          *stabilize_steps* value = timesteps
            timesteps = number of timesteps to apply the internally-created :doc:`nve/limit <fix_nve_limit>` fix to reacting atoms
          *update_edges* value = *none* or *charges* or *custom*
-           none = do not update topology near the edges of reaction templates
-           charges = update atomic charges of all atoms in reaction templates
-           custom = force the update of user-specified atomic charges
+           *none* = do not update topology near the edges of reaction templates
+           *charges* = update atomic charges of all atoms in reaction templates
+           *custom* = force the update of user-specified atomic charges
 
 Examples
 """"""""
@@ -154,6 +157,13 @@ due to the internal dynamic grouping performed by fix bond/react.
    If the group-ID is an existing static group, react-group-IDs
    should also be specified as this static group, or a subset.
 
+The *reset_mol_ids* keyword invokes the :doc:`reset_mol_ids <reset_mol_ids>`
+command after a reaction occurs, to ensure that molecule IDs are
+consistent with the new bond topology. The group-ID used for
+:doc:`reset_mol_ids <reset_mol_ids>` is the group-ID for this fix.
+Resetting molecule IDs is necessarily a global operation, and so can
+be slow for very large systems.
+
 The following comments pertain to each *react* argument (in other
 words, can be customized for each reaction, or reaction step):
 
@@ -553,7 +563,7 @@ Default
 """""""
 
 The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60,
-update_edges = none
+reset_mol_ids = yes, update_edges = none
 
 ----------
 

From da91f81d4062b4acad75c25af8523ea41edb91de Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 18 Jul 2020 13:42:47 -0600
Subject: [PATCH 03/21] bond/react:doc clarification

---
 doc/src/fix_bond_react.rst | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst
index d1671d7c78..82ab6beaa3 100644
--- a/doc/src/fix_bond_react.rst
+++ b/doc/src/fix_bond_react.rst
@@ -213,9 +213,9 @@ surrounding topology. As described below, the bonding atom pairs of
 the pre-reacted template are specified by atom ID in the map file. The
 pre-reacted molecule template should contain as few atoms as possible
 while still completely describing the topology of all atoms affected
-by the reaction. For example, if the force field contains dihedrals,
-the pre-reacted template should contain any atom within three bonds of
-reacting atoms.
+by the reaction (which includes all atoms that change atom type). For
+example, if the force field contains dihedrals, the pre-reacted
+template should contain any atom within three bonds of reacting atoms.
 
 Some atoms in the pre-reacted template that are not reacting may have
 missing topology with respect to the simulation. For example, the

From 6272b7d2bf55e9b8637ae33591d9a0b19d29c5e7 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 18 Jul 2020 13:52:13 -0600
Subject: [PATCH 04/21] add (undocumented) verbosity option to reset_mol_ids

---
 src/reset_mol_ids.cpp | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp
index e3cc877548..7f71aee0b4 100644
--- a/src/reset_mol_ids.cpp
+++ b/src/reset_mol_ids.cpp
@@ -56,6 +56,7 @@ void ResetMolIDs::command(int narg, char **arg)
   int compressflag = 1;
   int singleflag = 0;
   tagint offset = -1;
+  int verbose = 1;
 
   int iarg = 1;
   while (iarg < narg) {
@@ -76,6 +77,10 @@ void ResetMolIDs::command(int narg, char **arg)
       offset = utils::tnumeric(FLERR,arg[iarg+1],1,lmp);
       if (offset < -1) error->all(FLERR,"Illegal reset_mol_ids command");
       iarg += 2;
+    } else if (strcmp(arg[iarg],"verbose") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command");
+      verbose = utils::tnumeric(FLERR,arg[iarg+1],1,lmp);
+      iarg += 2;
     } else error->all(FLERR,"Illegal reset_mol_ids command");
   }
 
@@ -203,7 +208,7 @@ void ResetMolIDs::command(int narg, char **arg)
 
   MPI_Barrier(world);
 
-  if (comm->me == 0) {
+  if (verbose == 1 && comm->me == 0) {
     if (nchunk < 0)
       utils::logmesg(lmp,fmt::format("  number of new molecule IDs = unknown\n"));
     else

From a2547701e6bec294d362685442b6252b033ddd5f Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 18 Jul 2020 13:59:30 -0600
Subject: [PATCH 05/21] fix verbose reset_mol_ids

---
 src/reset_mol_ids.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp
index 7f71aee0b4..f62d1bf402 100644
--- a/src/reset_mol_ids.cpp
+++ b/src/reset_mol_ids.cpp
@@ -84,7 +84,7 @@ void ResetMolIDs::command(int narg, char **arg)
     } else error->all(FLERR,"Illegal reset_mol_ids command");
   }
 
-  if (comm->me == 0) utils::logmesg(lmp,"Resetting molecule IDs ...\n");
+  if (verbose == 1 && comm->me == 0) utils::logmesg(lmp,"Resetting molecule IDs ...\n");
 
   // record wall time for resetting molecule IDs
 

From e00b0e96f6a3162518ee05f48716dba3a3afe89e Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 18 Jul 2020 14:00:46 -0600
Subject: [PATCH 06/21] bond/react: prevent reset_mol_ids printing

---
 src/USER-REACTION/fix_bond_react.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index c908526067..ea64fb5447 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -2518,7 +2518,7 @@ update molecule IDs, charges, types, special lists and all topology
 void FixBondReact::update_everything()
 {
   if (reset_mol_ids_flag)
-    input->one("reset_mol_ids " + std::string(group->names[igroup]));
+    input->one("reset_mol_ids " + std::string(group->names[igroup]) + " verbose 0");
 
   int *type = atom->type;
 

From 9ec5708f2fb878a2bfc8cc7b1cd02e2997a96430 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 18 Jul 2020 14:21:10 -0600
Subject: [PATCH 07/21] Update reset_mol_ids.cpp

---
 src/reset_mol_ids.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp
index f62d1bf402..ff663ddf24 100644
--- a/src/reset_mol_ids.cpp
+++ b/src/reset_mol_ids.cpp
@@ -79,7 +79,7 @@ void ResetMolIDs::command(int narg, char **arg)
       iarg += 2;
     } else if (strcmp(arg[iarg],"verbose") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command");
-      verbose = utils::tnumeric(FLERR,arg[iarg+1],1,lmp);
+      verbose = utils::inumeric(FLERR,arg[iarg+1],1,lmp);
       iarg += 2;
     } else error->all(FLERR,"Illegal reset_mol_ids command");
   }

From fe6efe8861e287510931835b64dba97bc010946e Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 18 Jul 2020 14:29:39 -0600
Subject: [PATCH 08/21] need header file!

---
 src/USER-REACTION/fix_bond_react.h | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/USER-REACTION/fix_bond_react.h b/src/USER-REACTION/fix_bond_react.h
index e8a4253ce7..01e8a5ee96 100644
--- a/src/USER-REACTION/fix_bond_react.h
+++ b/src/USER-REACTION/fix_bond_react.h
@@ -61,6 +61,7 @@ class FixBondReact : public Fix {
   int *max_rxn,*nlocalskips,*nghostlyskips;
   tagint lastcheck;
   int stabilization_flag;
+  int reset_mol_ids_flag;
   int custom_exclude_flag;
   int *stabilize_steps_flag;
   int *update_edges_flag;

From edd3fb7108682a853dee3335f57bccf712af4c24 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 18 Jul 2020 20:51:14 -0600
Subject: [PATCH 09/21] reset_mol_ids: documented verbose option

---
 doc/src/reset_mol_ids.rst            | 9 +++++++--
 src/USER-REACTION/fix_bond_react.cpp | 2 +-
 src/reset_mol_ids.cpp                | 8 +++++---
 3 files changed, 13 insertions(+), 6 deletions(-)

diff --git a/doc/src/reset_mol_ids.rst b/doc/src/reset_mol_ids.rst
index 0d6063b3ef..9554ac86c2 100644
--- a/doc/src/reset_mol_ids.rst
+++ b/doc/src/reset_mol_ids.rst
@@ -19,6 +19,7 @@ Syntax
        *compress* value = *yes* or *no*
        *offset* value = *Noffset* >= -1
        *single* value = *yes* or *no* to treat single atoms (no bonds) as molecules
+       *verbose* value = *yes* or *no*
 
 Examples
 """"""""
@@ -95,6 +96,10 @@ is not *all* there may be collisions with the molecule IDs of other atoms.
    does not perform a full update of the bond topology data structures
    within LAMMPS.
 
+The *verbose* keyword determines if this command outputs run-time
+information to the screen and log file, including the number of new
+molecule IDs and CPU time used by the command.
+
 Restrictions
 """"""""""""
 none
@@ -112,5 +117,5 @@ Related commands
 Default
 """""""
 
-The default keyword settings are compress = yes, single = no, and
-offset = -1.
+The default keyword settings are compress = yes, single = no,
+verbose = yes, and offset = -1.
diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index ea64fb5447..54e9720ea7 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -2518,7 +2518,7 @@ update molecule IDs, charges, types, special lists and all topology
 void FixBondReact::update_everything()
 {
   if (reset_mol_ids_flag)
-    input->one("reset_mol_ids " + std::string(group->names[igroup]) + " verbose 0");
+    input->one("reset_mol_ids " + std::string(group->names[igroup]) + " verbose no");
 
   int *type = atom->type;
 
diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp
index ff663ddf24..7caa8f1e77 100644
--- a/src/reset_mol_ids.cpp
+++ b/src/reset_mol_ids.cpp
@@ -79,12 +79,14 @@ void ResetMolIDs::command(int narg, char **arg)
       iarg += 2;
     } else if (strcmp(arg[iarg],"verbose") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command");
-      verbose = utils::inumeric(FLERR,arg[iarg+1],1,lmp);
+      if (strcmp(arg[iarg+1],"yes") == 0) verbose = 1;
+      else if (strcmp(arg[iarg+1],"no") == 0) verbose = 0;
+      else error->all(FLERR,"Illegal reset_mol_ids command");
       iarg += 2;
     } else error->all(FLERR,"Illegal reset_mol_ids command");
   }
 
-  if (verbose == 1 && comm->me == 0) utils::logmesg(lmp,"Resetting molecule IDs ...\n");
+  if (verbose && (comm->me == 0)) utils::logmesg(lmp,"Resetting molecule IDs ...\n");
 
   // record wall time for resetting molecule IDs
 
@@ -208,7 +210,7 @@ void ResetMolIDs::command(int narg, char **arg)
 
   MPI_Barrier(world);
 
-  if (verbose == 1 && comm->me == 0) {
+  if (verbose && (comm->me == 0)) {
     if (nchunk < 0)
       utils::logmesg(lmp,fmt::format("  number of new molecule IDs = unknown\n"));
     else

From 845b9185010b12b8a655c487f9d134d5acc45f52 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 18 Jul 2020 20:59:17 -0600
Subject: [PATCH 10/21] probably better reset_mol_id doc version

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

diff --git a/doc/src/reset_mol_ids.rst b/doc/src/reset_mol_ids.rst
index 9554ac86c2..d0e8bc179a 100644
--- a/doc/src/reset_mol_ids.rst
+++ b/doc/src/reset_mol_ids.rst
@@ -97,8 +97,10 @@ is not *all* there may be collisions with the molecule IDs of other atoms.
    within LAMMPS.
 
 The *verbose* keyword determines if this command outputs run-time
-information to the screen and log file, including the number of new
-molecule IDs and CPU time used by the command.
+information to the screen and log file. If the setting is *yes*
+(the default), the command prints out information including the number
+of new molecule IDs and CPU time used by the command. If the setting
+is *no*, no information is printed.
 
 Restrictions
 """"""""""""

From 295d75f230d231f5a696fc5d9400a90a4f9f3883 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 18 Jul 2020 21:14:36 -0600
Subject: [PATCH 11/21] guess should reset_mol_IDs *after* updating bonds

---
 src/USER-REACTION/fix_bond_react.cpp | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index 54e9720ea7..6ce19bf22a 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -2517,9 +2517,6 @@ update molecule IDs, charges, types, special lists and all topology
 
 void FixBondReact::update_everything()
 {
-  if (reset_mol_ids_flag)
-    input->one("reset_mol_ids " + std::string(group->names[igroup]) + " verbose no");
-
   int *type = atom->type;
 
   int nlocal = atom->nlocal;
@@ -3054,6 +3051,10 @@ void FixBondReact::update_everything()
   atom->natoms -= ndel;
   // done deleting atoms
 
+  // reset mol ids
+  if (reset_mol_ids_flag)
+    input->one("reset_mol_ids " + std::string(group->names[igroup]) + " verbose no");
+
   // something to think about: this could done much more concisely if
   // all atom-level info (bond,angles, etc...) were kinda inherited from a common data struct --JG
 

From cf83ce674548f4424166d095816f66f47ed17c30 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 25 Jul 2020 00:52:39 -0600
Subject: [PATCH 12/21] reset_mol_ids->reset() version

---
 doc/src/reset_mol_ids.rst            | 11 +---
 src/USER-REACTION/fix_bond_react.cpp | 13 ++++-
 src/USER-REACTION/fix_bond_react.h   |  1 +
 src/reset_mol_ids.cpp                | 84 +++++++++++++++-------------
 src/reset_mol_ids.h                  |  7 +++
 5 files changed, 66 insertions(+), 50 deletions(-)

diff --git a/doc/src/reset_mol_ids.rst b/doc/src/reset_mol_ids.rst
index d0e8bc179a..0d6063b3ef 100644
--- a/doc/src/reset_mol_ids.rst
+++ b/doc/src/reset_mol_ids.rst
@@ -19,7 +19,6 @@ Syntax
        *compress* value = *yes* or *no*
        *offset* value = *Noffset* >= -1
        *single* value = *yes* or *no* to treat single atoms (no bonds) as molecules
-       *verbose* value = *yes* or *no*
 
 Examples
 """"""""
@@ -96,12 +95,6 @@ is not *all* there may be collisions with the molecule IDs of other atoms.
    does not perform a full update of the bond topology data structures
    within LAMMPS.
 
-The *verbose* keyword determines if this command outputs run-time
-information to the screen and log file. If the setting is *yes*
-(the default), the command prints out information including the number
-of new molecule IDs and CPU time used by the command. If the setting
-is *no*, no information is printed.
-
 Restrictions
 """"""""""""
 none
@@ -119,5 +112,5 @@ Related commands
 Default
 """""""
 
-The default keyword settings are compress = yes, single = no,
-verbose = yes, and offset = -1.
+The default keyword settings are compress = yes, single = no, and
+offset = -1.
diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index 6ce19bf22a..16a9ca29f1 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -33,6 +33,7 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu)
 #include "neigh_list.h"
 #include "neigh_request.h"
 #include "random_mars.h"
+#include "reset_mol_ids.h"
 #include "molecule.h"
 #include "group.h"
 #include "citeme.h"
@@ -92,6 +93,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
   fix1 = NULL;
   fix2 = NULL;
   fix3 = NULL;
+  reset_mol_ids = NULL;
 
   if (narg < 8) error->all(FLERR,"Illegal fix bond/react command: "
                            "too few arguments");
@@ -166,7 +168,11 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
     } else if (strcmp(arg[iarg],"reset_mol_ids") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: "
                                     "'reset_mol_ids' keyword has too few arguments");
-      if (strcmp(arg[iarg+1],"yes") == 0) iarg += 2; //default
+      if (strcmp(arg[iarg+1],"yes") == 0) { // default
+        delete reset_mol_ids;
+        reset_mol_ids = new ResetMolIDs(lmp);
+        iarg += 2;
+      }
       if (strcmp(arg[iarg+1],"no") == 0) {
         reset_mol_ids_flag = 0;
         iarg += 2;
@@ -508,6 +514,8 @@ FixBondReact::~FixBondReact()
   }
   delete [] random;
 
+  delete reset_mol_ids;
+
   memory->destroy(partner);
   memory->destroy(finalpartner);
   memory->destroy(ncreate);
@@ -3052,8 +3060,7 @@ void FixBondReact::update_everything()
   // done deleting atoms
 
   // reset mol ids
-  if (reset_mol_ids_flag)
-    input->one("reset_mol_ids " + std::string(group->names[igroup]) + " verbose no");
+  if (reset_mol_ids_flag) reset_mol_ids->reset(group->names[igroup]);
 
   // something to think about: this could done much more concisely if
   // all atom-level info (bond,angles, etc...) were kinda inherited from a common data struct --JG
diff --git a/src/USER-REACTION/fix_bond_react.h b/src/USER-REACTION/fix_bond_react.h
index 01e8a5ee96..1b98614064 100644
--- a/src/USER-REACTION/fix_bond_react.h
+++ b/src/USER-REACTION/fix_bond_react.h
@@ -94,6 +94,7 @@ class FixBondReact : public Fix {
   class RanMars **random; // random number for 'prob' keyword
   class RanMars **rrhandom; // random number for Arrhenius constraint
   class NeighList *list;
+  class ResetMolIDs *reset_mol_ids; // class for resetting mol IDs
 
   int *reacted_mol,*unreacted_mol;
   int *limit_duration; // indicates how long to relax
diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp
index 7caa8f1e77..e6add7df9b 100644
--- a/src/reset_mol_ids.cpp
+++ b/src/reset_mol_ids.cpp
@@ -33,9 +33,15 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-ResetMolIDs::ResetMolIDs(LAMMPS *lmp) : Pointers(lmp) {}
+ResetMolIDs::ResetMolIDs(LAMMPS *lmp) : Pointers(lmp) {
+  compressflag = 1;
+  singleflag = 0;
+  offset = -1;
+}
 
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   called as reset_mol_ids command in input script
+------------------------------------------------------------------------- */
 
 void ResetMolIDs::command(int narg, char **arg)
 {
@@ -49,14 +55,7 @@ void ResetMolIDs::command(int narg, char **arg)
   // process args
 
   if (narg < 1) error->all(FLERR,"Illegal reset_mol_ids command");
-  int igroup = group->find(arg[0]);
-  if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID");
-  int groupbit = group->bitmask[igroup];
-
-  int compressflag = 1;
-  int singleflag = 0;
-  tagint offset = -1;
-  int verbose = 1;
+  char *groupid = arg[0];
 
   int iarg = 1;
   while (iarg < narg) {
@@ -77,39 +76,61 @@ void ResetMolIDs::command(int narg, char **arg)
       offset = utils::tnumeric(FLERR,arg[iarg+1],1,lmp);
       if (offset < -1) error->all(FLERR,"Illegal reset_mol_ids command");
       iarg += 2;
-    } else if (strcmp(arg[iarg],"verbose") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command");
-      if (strcmp(arg[iarg+1],"yes") == 0) verbose = 1;
-      else if (strcmp(arg[iarg+1],"no") == 0) verbose = 0;
-      else error->all(FLERR,"Illegal reset_mol_ids command");
-      iarg += 2;
     } else error->all(FLERR,"Illegal reset_mol_ids command");
   }
 
-  if (verbose && (comm->me == 0)) utils::logmesg(lmp,"Resetting molecule IDs ...\n");
+  if (comm->me == 0) utils::logmesg(lmp,"Resetting molecule IDs ...\n");
 
   // record wall time for resetting molecule IDs
 
   MPI_Barrier(world);
   double time1 = MPI_Wtime();
 
+  // initialize system since comm->borders() will be invoked
+
+  lmp->init();
+
+  // reset molecule IDs
+
+  reset(groupid);
+
+  // total time
+
+  MPI_Barrier(world);
+
+  if (comm->me == 0) {
+    if (nchunk < 0)
+      utils::logmesg(lmp,fmt::format("  number of new molecule IDs = unknown\n"));
+    else
+      utils::logmesg(lmp,fmt::format("  number of new molecule IDs = {}\n",nchunk));
+    utils::logmesg(lmp,fmt::format("  reset_mol_ids CPU = {:.3f} seconds\n",
+                                   MPI_Wtime()-time1));
+  }
+}
+
+/* ----------------------------------------------------------------------
+   called from command() and directly from fixes that update molecules
+------------------------------------------------------------------------- */
+
+void ResetMolIDs::reset(char *groupid)
+{
+  int igroup = group->find(groupid);
+  if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID");
+  int groupbit = group->bitmask[igroup];
+
   // create instances of compute fragment/atom, compute reduce (if needed),
   // and compute chunk/atom.  all use the group-ID for this command
 
   const std::string idfrag = "reset_mol_ids_FRAGMENT_ATOM";
   if (singleflag)
-    modify->add_compute(fmt::format("{} {} fragment/atom single yes",idfrag,arg[0]));
+    modify->add_compute(fmt::format("{} {} fragment/atom single yes",idfrag,groupid));
   else
-    modify->add_compute(fmt::format("{} {} fragment/atom single no",idfrag,arg[0]));
+    modify->add_compute(fmt::format("{} {} fragment/atom single no",idfrag,groupid));
 
   const std::string idchunk = "reset_mol_ids_CHUNK_ATOM";
   if (compressflag)
     modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes",
-                                    idchunk,arg[0]));
-
-  // initialize system since comm->borders() will be invoked
-
-  lmp->init();
+                                    idchunk,groupid));
 
   // setup domain, communication
   // exchange will clear map, borders will reset
@@ -145,7 +166,7 @@ void ResetMolIDs::command(int narg, char **arg)
   // if compressflag = 0, done
   // set nchunk = -1 since cannot easily determine # of new molecule IDs
 
-  int nchunk = -1;
+  nchunk = -1;
 
   // if compressflag = 1, invoke peratom method of compute chunk/atom
   // will compress new molecule IDs to be contiguous 1 to Nmol
@@ -205,17 +226,4 @@ void ResetMolIDs::command(int narg, char **arg)
 
   modify->delete_compute(idfrag);
   if (compressflag) modify->delete_compute(idchunk);
-
-  // total time
-
-  MPI_Barrier(world);
-
-  if (verbose && (comm->me == 0)) {
-    if (nchunk < 0)
-      utils::logmesg(lmp,fmt::format("  number of new molecule IDs = unknown\n"));
-    else
-      utils::logmesg(lmp,fmt::format("  number of new molecule IDs = {}\n",nchunk));
-    utils::logmesg(lmp,fmt::format("  reset_mol_ids CPU = {:.3f} seconds\n",
-                                   MPI_Wtime()-time1));
-  }
 }
diff --git a/src/reset_mol_ids.h b/src/reset_mol_ids.h
index 83c0b5d80f..cb7c472724 100644
--- a/src/reset_mol_ids.h
+++ b/src/reset_mol_ids.h
@@ -28,6 +28,13 @@ class ResetMolIDs : protected Pointers {
  public:
   ResetMolIDs(class LAMMPS *);
   void command(int, char **);
+  void reset(char *);
+
+private:
+  int nchunk;
+  int compressflag; // 1 = contiguous values for new IDs
+  int singleflag; // 0 = mol IDs of single atoms set to 0
+  tagint offset; // offset for contiguous mol ID values
 };
 
 }

From f6d91b3b2c92d4f8c820265a3bcea5a6649c9daa Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Tue, 11 Aug 2020 15:02:37 -0400
Subject: [PATCH 13/21] move domain/comm commands

---
 src/reset_mol_ids.cpp | 24 ++++++++++++------------
 1 file changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp
index e6add7df9b..aa0fe37ce7 100644
--- a/src/reset_mol_ids.cpp
+++ b/src/reset_mol_ids.cpp
@@ -90,6 +90,18 @@ void ResetMolIDs::command(int narg, char **arg)
 
   lmp->init();
 
+  // setup domain, communication
+  // exchange will clear map, borders will reset
+  // this is the map needed to lookup current global IDs for bond topology
+
+  if (domain->triclinic) domain->x2lamda(atom->nlocal);
+  domain->pbc();
+  domain->reset_box();
+  comm->setup();
+  comm->exchange();
+  comm->borders();
+  if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
+
   // reset molecule IDs
 
   reset(groupid);
@@ -132,18 +144,6 @@ void ResetMolIDs::reset(char *groupid)
     modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes",
                                     idchunk,groupid));
 
-  // setup domain, communication
-  // exchange will clear map, borders will reset
-  // this is the map needed to lookup current global IDs for bond topology
-
-  if (domain->triclinic) domain->x2lamda(atom->nlocal);
-  domain->pbc();
-  domain->reset_box();
-  comm->setup();
-  comm->exchange();
-  comm->borders();
-  if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
-
   // invoke peratom method of compute fragment/atom
   // walks bond connectivity and assigns each atom a fragment ID
   // if singleflag = 0, atoms w/out bonds will be assigned fragID = 0

From 3c69ebc6696f83f81f5b6f28bce35193d5bf93c6 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Tue, 11 Aug 2020 17:12:36 -0400
Subject: [PATCH 14/21] reset_mold_ids: add create_computes

---
 src/reset_mol_ids.cpp | 55 ++++++++++++++++++++++++++++++-------------
 src/reset_mol_ids.h   | 10 +++++++-
 2 files changed, 48 insertions(+), 17 deletions(-)

diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp
index aa0fe37ce7..efeb6a672a 100644
--- a/src/reset_mol_ids.cpp
+++ b/src/reset_mol_ids.cpp
@@ -17,7 +17,6 @@
 
 #include "reset_mol_ids.h"
 #include <mpi.h>
-#include <string>
 #include "atom.h"
 #include "domain.h"
 #include "comm.h"
@@ -34,11 +33,24 @@ using namespace LAMMPS_NS;
 /* ---------------------------------------------------------------------- */
 
 ResetMolIDs::ResetMolIDs(LAMMPS *lmp) : Pointers(lmp) {
+  cfa = NULL;
+  cca = NULL;
+
+  // default settings
+
   compressflag = 1;
   singleflag = 0;
   offset = -1;
 }
 
+/* ---------------------------------------------------------------------- */
+
+ResetMolIDs::~ResetMolIDs()
+{
+  modify->delete_compute(idfrag);
+  if (compressflag) modify->delete_compute(idchunk);
+}
+
 /* ----------------------------------------------------------------------
    called as reset_mol_ids command in input script
 ------------------------------------------------------------------------- */
@@ -102,9 +114,13 @@ void ResetMolIDs::command(int narg, char **arg)
   comm->borders();
   if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
 
+  // create computes
+
+  create_computes(groupid);
+
   // reset molecule IDs
 
-  reset(groupid);
+  reset();
 
   // total time
 
@@ -121,35 +137,49 @@ void ResetMolIDs::command(int narg, char **arg)
 }
 
 /* ----------------------------------------------------------------------
-   called from command() and directly from fixes that update molecules
+   create computes used by reset_mol_ids
 ------------------------------------------------------------------------- */
 
-void ResetMolIDs::reset(char *groupid)
+void ResetMolIDs::create_computes(char *groupid)
 {
   int igroup = group->find(groupid);
   if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID");
-  int groupbit = group->bitmask[igroup];
+  groupbit = group->bitmask[igroup];
 
   // create instances of compute fragment/atom, compute reduce (if needed),
   // and compute chunk/atom.  all use the group-ID for this command
 
-  const std::string idfrag = "reset_mol_ids_FRAGMENT_ATOM";
+  idfrag = "reset_mol_ids_FRAGMENT_ATOM";
   if (singleflag)
     modify->add_compute(fmt::format("{} {} fragment/atom single yes",idfrag,groupid));
   else
     modify->add_compute(fmt::format("{} {} fragment/atom single no",idfrag,groupid));
 
-  const std::string idchunk = "reset_mol_ids_CHUNK_ATOM";
+  idchunk = "reset_mol_ids_CHUNK_ATOM";
   if (compressflag)
     modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes",
                                     idchunk,groupid));
 
+  int icompute = modify->find_compute(idfrag);
+  cfa = (ComputeFragmentAtom *) modify->compute[icompute];
+
+
+  if (compressflag) {
+    icompute = modify->find_compute(idchunk);
+    cca = (ComputeChunkAtom *) modify->compute[icompute];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   called from command() and directly from fixes that update molecules
+------------------------------------------------------------------------- */
+
+void ResetMolIDs::reset()
+{
   // invoke peratom method of compute fragment/atom
   // walks bond connectivity and assigns each atom a fragment ID
   // if singleflag = 0, atoms w/out bonds will be assigned fragID = 0
 
-  int icompute = modify->find_compute(idfrag);
-  ComputeFragmentAtom *cfa = (ComputeFragmentAtom *) modify->compute[icompute];
   cfa->compute_peratom();
   double *fragIDs = cfa->vector_atom;
 
@@ -173,8 +203,6 @@ void ResetMolIDs::reset(char *groupid)
   // NOTE: use of compute chunk/atom limits Nmol to a 32-bit int
 
   if (compressflag) {
-    icompute = modify->find_compute(idchunk);
-    ComputeChunkAtom *cca = (ComputeChunkAtom *) modify->compute[icompute];
     cca->compute_peratom();
     double *chunkIDs = cca->vector_atom;
     nchunk = cca->nchunk;
@@ -221,9 +249,4 @@ void ResetMolIDs::reset(char *groupid)
       }
     }
   }
-
-  // clean up
-
-  modify->delete_compute(idfrag);
-  if (compressflag) modify->delete_compute(idchunk);
 }
diff --git a/src/reset_mol_ids.h b/src/reset_mol_ids.h
index cb7c472724..62543f6eda 100644
--- a/src/reset_mol_ids.h
+++ b/src/reset_mol_ids.h
@@ -21,20 +21,28 @@ CommandStyle(reset_mol_ids,ResetMolIDs)
 #define LMP_RESET_MOL_IDS_H
 
 #include "pointers.h"
+#include <string>
 
 namespace LAMMPS_NS {
 
 class ResetMolIDs : protected Pointers {
  public:
   ResetMolIDs(class LAMMPS *);
+  ~ResetMolIDs();
   void command(int, char **);
-  void reset(char *);
+  void create_computes(char *);
+  void reset();
 
 private:
+  std::string idfrag, idchunk;
   int nchunk;
+  int groupbit;
   int compressflag; // 1 = contiguous values for new IDs
   int singleflag; // 0 = mol IDs of single atoms set to 0
   tagint offset; // offset for contiguous mol ID values
+
+  class ComputeFragmentAtom *cfa;
+  class ComputeChunkAtom *cca;
 };
 
 }

From 9d486d734b3365e835f3b28b17eaa115c2f741b1 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Tue, 11 Aug 2020 17:29:27 -0400
Subject: [PATCH 15/21] update bond/react for reset_mol_ids->create_computes

---
 src/USER-REACTION/fix_bond_react.cpp | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index 16a9ca29f1..9b33630b59 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -171,6 +171,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
       if (strcmp(arg[iarg+1],"yes") == 0) { // default
         delete reset_mol_ids;
         reset_mol_ids = new ResetMolIDs(lmp);
+        reset_mol_ids->create_computes(group->names[igroup]);
         iarg += 2;
       }
       if (strcmp(arg[iarg+1],"no") == 0) {
@@ -246,9 +247,9 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
     if (n > MAXLINE) error->all(FLERR,"Reaction name (react-ID) is too long (limit: 256 characters)");
     strncpy(rxn_name[rxn],arg[iarg++],n);
 
-    int igroup = group->find(arg[iarg++]);
-    if (igroup == -1) error->all(FLERR,"Could not find fix group ID");
-    groupbits[rxn] = group->bitmask[igroup];
+    int groupid = group->find(arg[iarg++]);
+    if (groupid == -1) error->all(FLERR,"Could not find fix group ID");
+    groupbits[rxn] = group->bitmask[groupid];
 
     if (strncmp(arg[iarg],"v_",2) == 0) {
       n = strlen(&arg[iarg][2]) + 1;
@@ -640,9 +641,9 @@ void FixBondReact::post_constructor()
   group->assign(cmd);
 
   if (stabilization_flag == 1) {
-    int igroup = group->find(exclude_group);
+    int groupid = group->find(exclude_group);
     // create exclude_group if not already existing, or use as parent group if static
-    if (igroup == -1 || group->dynamic[igroup] == 0) {
+    if (groupid == -1 || group->dynamic[groupid] == 0) {
       // create stabilization per-atom property
       cmd = std::string("bond_react_stabilization_internal");
       id_fix3 = new char[cmd.size()+1];
@@ -672,7 +673,7 @@ void FixBondReact::post_constructor()
       strcat(exclude_group,"_REACT");
 
       group->find_or_create(exclude_group);
-      if (igroup == -1)
+      if (groupid == -1)
         cmd = fmt::format("{} dynamic all property statted_tags",exclude_group);
       else
         cmd = fmt::format("{} dynamic {} property statted_tags",exclude_group,exclude_PARENT_group);
@@ -3060,7 +3061,7 @@ void FixBondReact::update_everything()
   // done deleting atoms
 
   // reset mol ids
-  if (reset_mol_ids_flag) reset_mol_ids->reset(group->names[igroup]);
+  if (reset_mol_ids_flag) reset_mol_ids->reset();
 
   // something to think about: this could done much more concisely if
   // all atom-level info (bond,angles, etc...) were kinda inherited from a common data struct --JG

From df497e4853ac26faf7d13d43efbcb7bfdc254aec Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Tue, 18 Aug 2020 17:53:07 -0400
Subject: [PATCH 16/21] bond/react: clarify bond-type-checking error

---
 src/USER-REACTION/fix_bond_react.cpp | 6 +++---
 src/USER-REACTION/fix_bond_react.h   | 4 +++-
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index 9b33630b59..57bc42cffd 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -2079,7 +2079,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
   for (int i = 0; i < twomol->natoms; i++) {
     if (twomol->type[i] != onemol->type[equivalences[i][1][myrxn]-1] && landlocked_atoms[i][myrxn] == 0) {
       char str[128];
-      snprintf(str,128,"Bond/react: Atom affected by reaction %s too close to template edge",rxn_name[myrxn]);
+      snprintf(str,128,"Bond/react: Atom type affected by reaction %s too close to template edge",rxn_name[myrxn]);
       error->all(FLERR,str);
     }
   }
@@ -2098,7 +2098,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
             if (onemol_batom == equivalences[twomol_atomj-1][1][myrxn]) {
               if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomi-1][m]) {
                 char str[128];
-                snprintf(str,128,"Bond/react: Atom affected by reaction %s too close to template edge",rxn_name[myrxn]);
+                snprintf(str,128,"Bond/react: Bond type affected by reaction %s too close to template edge",rxn_name[myrxn]);
                 error->all(FLERR,str);
               }
             }
@@ -2110,7 +2110,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
               if (onemol_batom == equivalences[i][1][myrxn]) {
                 if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomj-1][m]) {
                   char str[128];
-                  snprintf(str,128,"Bond/react: Atom affected by reaction %s too close to template edge",rxn_name[myrxn]);
+                  snprintf(str,128,"Bond/react: Bond type affected by reaction %s too close to template edge",rxn_name[myrxn]);
                   error->all(FLERR,str);
                 }
               }
diff --git a/src/USER-REACTION/fix_bond_react.h b/src/USER-REACTION/fix_bond_react.h
index 1b98614064..61a1bd4213 100644
--- a/src/USER-REACTION/fix_bond_react.h
+++ b/src/USER-REACTION/fix_bond_react.h
@@ -242,8 +242,10 @@ E: Bond/react: Invalid template atom ID in map file
 Atom IDs in molecule templates range from 1 to the number of atoms in the template.
 
 E or W: Bond/react: Atom affected by reaction %s too close to template edge
+        Bond/react: Atom type affected by reaction %s too close to template edge
+        Bond/react: Bond type affected by reaction %s too close to template edge
 
-This means an atom which changes type or connectivity during the
+This means an atom (or bond) that changes type or connectivity during the
 reaction is too close to an 'edge' atom defined in the map file. This
 could cause incorrect assignment of bonds, angle, etc. Generally, this
 means you must include more atoms in your templates, such that there

From e4ab49c2e5d0856b78047088f0442be2d801bb8f Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Tue, 18 Aug 2020 18:12:01 -0400
Subject: [PATCH 17/21] bond/react: bond-type-checking docs

---
 doc/src/Errors_messages.rst | 3 ++-
 doc/src/fix_bond_react.rst  | 7 ++++---
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/doc/src/Errors_messages.rst b/doc/src/Errors_messages.rst
index f3be94a239..96e21da681 100644
--- a/doc/src/Errors_messages.rst
+++ b/doc/src/Errors_messages.rst
@@ -502,7 +502,8 @@ Doc page with :doc:`WARNING messages <Errors_warnings>`
 *Bond/react: Unknown section in map file*
    Please ensure reaction map files are properly formatted.
 
-*Bond/react: Atom affected by reaction too close to template edge*
+*Bond/react: Atom type affected by reaction too close to template edge*
+*Bond/react: Bond type affected by reaction too close to template edge*
    This means an atom which changes type or connectivity during the
    reaction is too close to an 'edge' atom defined in the map
    file. This could cause incorrect assignment of bonds, angle, etc.
diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst
index 82ab6beaa3..bbc6fbf864 100644
--- a/doc/src/fix_bond_react.rst
+++ b/doc/src/fix_bond_react.rst
@@ -213,9 +213,10 @@ surrounding topology. As described below, the bonding atom pairs of
 the pre-reacted template are specified by atom ID in the map file. The
 pre-reacted molecule template should contain as few atoms as possible
 while still completely describing the topology of all atoms affected
-by the reaction (which includes all atoms that change atom type). For
-example, if the force field contains dihedrals, the pre-reacted
-template should contain any atom within three bonds of reacting atoms.
+by the reaction (which includes all atoms that change atom type or
+connectivity, and all bonds that change bond type). For example, if
+the force field contains dihedrals, the pre-reacted template should
+contain any atom within three bonds of reacting atoms.
 
 Some atoms in the pre-reacted template that are not reacting may have
 missing topology with respect to the simulation. For example, the

From 6fc2ab07ef4797372063c349247877d0d9d02cd8 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Fri, 21 Aug 2020 14:52:39 -0400
Subject: [PATCH 18/21] reset_mol_ids: unique created computes

---
 src/USER-REACTION/fix_bond_react.cpp |  2 +-
 src/reset_mol_ids.cpp                | 11 ++++++-----
 src/reset_mol_ids.h                  |  2 +-
 3 files changed, 8 insertions(+), 7 deletions(-)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index 57bc42cffd..a1154c9221 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -171,7 +171,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
       if (strcmp(arg[iarg+1],"yes") == 0) { // default
         delete reset_mol_ids;
         reset_mol_ids = new ResetMolIDs(lmp);
-        reset_mol_ids->create_computes(group->names[igroup]);
+        reset_mol_ids->create_computes(id,group->names[igroup]);
         iarg += 2;
       }
       if (strcmp(arg[iarg+1],"no") == 0) {
diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp
index efeb6a672a..1461a49f32 100644
--- a/src/reset_mol_ids.cpp
+++ b/src/reset_mol_ids.cpp
@@ -116,7 +116,7 @@ void ResetMolIDs::command(int narg, char **arg)
 
   // create computes
 
-  create_computes(groupid);
+  create_computes(NULL,groupid);
 
   // reset molecule IDs
 
@@ -140,22 +140,23 @@ void ResetMolIDs::command(int narg, char **arg)
    create computes used by reset_mol_ids
 ------------------------------------------------------------------------- */
 
-void ResetMolIDs::create_computes(char *groupid)
+void ResetMolIDs::create_computes(char *fixid, char *groupid)
 {
   int igroup = group->find(groupid);
   if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID");
   groupbit = group->bitmask[igroup];
 
   // create instances of compute fragment/atom, compute reduce (if needed),
-  // and compute chunk/atom.  all use the group-ID for this command
+  // and compute chunk/atom.  all use the group-ID for this command.
+  // 'fixid' allows for creating independent instances of the computes
 
-  idfrag = "reset_mol_ids_FRAGMENT_ATOM";
+  idfrag = fmt::format("{}_reset_mol_ids_FRAGMENT_ATOM",fixid);
   if (singleflag)
     modify->add_compute(fmt::format("{} {} fragment/atom single yes",idfrag,groupid));
   else
     modify->add_compute(fmt::format("{} {} fragment/atom single no",idfrag,groupid));
 
-  idchunk = "reset_mol_ids_CHUNK_ATOM";
+  idchunk = fmt::format("{}_reset_mol_ids_CHUNK_ATOM",fixid);
   if (compressflag)
     modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes",
                                     idchunk,groupid));
diff --git a/src/reset_mol_ids.h b/src/reset_mol_ids.h
index 62543f6eda..9d2dd24bc8 100644
--- a/src/reset_mol_ids.h
+++ b/src/reset_mol_ids.h
@@ -30,7 +30,7 @@ class ResetMolIDs : protected Pointers {
   ResetMolIDs(class LAMMPS *);
   ~ResetMolIDs();
   void command(int, char **);
-  void create_computes(char *);
+  void create_computes(char *, char *);
   void reset();
 
 private:

From 01dd80f35e7709bb15fcfeadbf6c6f9ceb960de6 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sun, 23 Aug 2020 11:21:43 -0400
Subject: [PATCH 19/21] bond/react: actually make reset_mol_ids the default

---
 src/USER-REACTION/fix_bond_react.cpp | 21 ++++++++++-----------
 1 file changed, 10 insertions(+), 11 deletions(-)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index a1154c9221..29aa476148 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -168,21 +168,20 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
     } else if (strcmp(arg[iarg],"reset_mol_ids") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: "
                                     "'reset_mol_ids' keyword has too few arguments");
-      if (strcmp(arg[iarg+1],"yes") == 0) { // default
-        delete reset_mol_ids;
-        reset_mol_ids = new ResetMolIDs(lmp);
-        reset_mol_ids->create_computes(id,group->names[igroup]);
-        iarg += 2;
-      }
-      if (strcmp(arg[iarg+1],"no") == 0) {
-        reset_mol_ids_flag = 0;
-        iarg += 2;
-      }
+      if (strcmp(arg[iarg+1],"yes") == 0) ; // default
+      if (strcmp(arg[iarg+1],"no") == 0) reset_mol_ids_flag = 0;
+      iarg += 2;
     } else if (strcmp(arg[iarg],"react") == 0) {
       break;
     } else error->all(FLERR,"Illegal fix bond/react command: unknown keyword");
   }
 
+  if (reset_mol_ids_flag) {
+    delete reset_mol_ids;
+    reset_mol_ids = new ResetMolIDs(lmp);
+    reset_mol_ids->create_computes(id,group->names[igroup]);
+  }
+
   // set up common variables as vectors of length 'nreacts'
   // nevery, cutoff, onemol, twomol, superimpose file
 
@@ -245,7 +244,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
 
     int n = strlen(arg[iarg]) + 1;
     if (n > MAXLINE) error->all(FLERR,"Reaction name (react-ID) is too long (limit: 256 characters)");
-    strncpy(rxn_name[rxn],arg[iarg++],n);
+    strcpy(rxn_name[rxn],arg[iarg++]);
 
     int groupid = group->find(arg[iarg++]);
     if (groupid == -1) error->all(FLERR,"Could not find fix group ID");

From 63abb2dff919327995dc350b70f44800176d7eef Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sun, 23 Aug 2020 11:32:54 -0400
Subject: [PATCH 20/21] fix broken reset_mol_ids command

---
 src/reset_mol_ids.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp
index 1461a49f32..de76ef1732 100644
--- a/src/reset_mol_ids.cpp
+++ b/src/reset_mol_ids.cpp
@@ -116,7 +116,7 @@ void ResetMolIDs::command(int narg, char **arg)
 
   // create computes
 
-  create_computes(NULL,groupid);
+  create_computes("COMMAND",groupid);
 
   // reset molecule IDs
 

From 10080079e3da86a39cc615a236242262710b91a9 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sun, 23 Aug 2020 11:44:48 -0400
Subject: [PATCH 21/21] ISO compliance

---
 src/reset_mol_ids.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp
index de76ef1732..7ce027a657 100644
--- a/src/reset_mol_ids.cpp
+++ b/src/reset_mol_ids.cpp
@@ -116,7 +116,7 @@ void ResetMolIDs::command(int narg, char **arg)
 
   // create computes
 
-  create_computes("COMMAND",groupid);
+  create_computes((char *) "COMMAND",groupid);
 
   // reset molecule IDs