From cb4549e0f28aab7ad1a5d6d884f44cf92b55ec18 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Wed, 4 Nov 2020 20:31:38 -0500
Subject: [PATCH 01/19] bond/react:bugfix for delete atoms feature

---
 src/USER-REACTION/fix_bond_react.cpp | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index 83a3033cc2..1d5a713f9c 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -3131,6 +3131,12 @@ void FixBondReact::update_everything()
   int Tdelta_imprp;
   MPI_Allreduce(&delta_imprp,&Tdelta_imprp,1,MPI_INT,MPI_SUM,world);
   atom->nimpropers += Tdelta_imprp;
+
+  if (ndel && (atom->map_style != Atom::MAP_NONE)) {
+    atom->nghost = 0; 
+    atom->map_init();
+    atom->map_set();
+  }
 }
 
 /* ----------------------------------------------------------------------

From 124feffafae1c16c80160925a06bf7e7f7536b97 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Wed, 4 Nov 2020 21:35:25 -0500
Subject: [PATCH 02/19] bond/react: add helpful hint for recurring mistake

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

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index 1d5a713f9c..b6c1e91fca 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -1304,7 +1304,9 @@ void FixBondReact::superimpose_algorithm()
         // let's go ahead and catch the simplest of hangs
         //if (hang_catch > onemol->natoms*4)
         if (hang_catch > atom->nlocal*30) {
-          error->one(FLERR,"Bond/react: Excessive iteration of superimpose algorithm");
+          error->one(FLERR,"Bond/react: Excessive iteration of superimpose algorithm. "
+              "Please check that all pre-reaction template atoms are linked to a bonding atom, "
+              "via at least one path that does not involve edge atoms.");
         }
       }
     }
@@ -3133,7 +3135,7 @@ void FixBondReact::update_everything()
   atom->nimpropers += Tdelta_imprp;
 
   if (ndel && (atom->map_style != Atom::MAP_NONE)) {
-    atom->nghost = 0; 
+    atom->nghost = 0;
     atom->map_init();
     atom->map_set();
   }

From badbb411eb2e48c38c67021fea666e3023e627e7 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sun, 8 Nov 2020 22:59:45 -0500
Subject: [PATCH 03/19] ringed structures bugfix

previously,  atoms at 'end' of symmetric ring could behave like edge atoms
---
 src/USER-REACTION/fix_bond_react.cpp | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index b6c1e91fca..7f32d579ff 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -1720,6 +1720,17 @@ void FixBondReact::ring_check()
   // ring_check can be made more efficient by re-introducing 'frozen' atoms
   // 'frozen' atoms have been assigned and also are no longer pioneers
 
+  // double check the number of neighbors match for all non-edge atoms
+  // otherwise, atoms at 'end' of symmetric ring can behave like edge atoms
+  for (int i = 0; i < onemol->natoms; i++) {
+    if (edge[i][rxnID] == 0) {
+      if (onemol_nxspecial[i][0] != nxspecial[atom->map(glove[i][1])][0]) {
+        status = GUESSFAIL;
+        return;
+      }
+    }
+  }
+
   for (int i = 0; i < onemol->natoms; i++) {
     for (int j = 0; j < onemol_nxspecial[i][0]; j++) {
       int ring_fail = 1;

From 4d493fd082610657792fe7ba6c051e7b6b940ac2 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Mon, 9 Nov 2020 00:05:20 -0500
Subject: [PATCH 04/19] relax cutoff checks

also fixes some issues when using hybrid pair styles
---
 src/USER-REACTION/fix_bond_react.cpp | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index 7f32d579ff..994c8cba91 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -459,6 +459,14 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
     }
   }
 
+  for (int i = 0; i < nreacts; i++) {
+    if (closeneigh[i] == -1) { // indicates will search non-bonded neighbors
+      if (cutsq[i][1] > neighbor->cutneighsq[iatomtype[i]][jatomtype[i]]) {
+        error->all(FLERR,"Bond/react: Fix bond/react cutoff is longer than pairwise cutoff");
+      }
+    }
+  }
+
   // initialize Marsaglia RNG with processor-unique seed ('prob' keyword)
 
   random = new class RanMars*[nreacts];
@@ -750,12 +758,6 @@ void FixBondReact::init()
   if (strstr(update->integrate_style,"respa"))
     nlevels_respa = ((Respa *) update->integrate)->nlevels;
 
-  // check cutoff for iatomtype,jatomtype
-  for (int i = 0; i < nreacts; i++) {
-    if (force->pair == nullptr || cutsq[i][1] > force->pair->cutsq[iatomtype[i]][jatomtype[i]])
-      error->all(FLERR,"Bond/react: Fix bond/react cutoff is longer than pairwise cutoff");
-  }
-
   // need a half neighbor list, built every Nevery steps
   int irequest = neighbor->request(this,instance_me);
   neighbor->requests[irequest]->pair = 0;

From 91d9cf97f3b438bb384eb65600f5702ace65baca Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Mon, 9 Nov 2020 13:47:29 -0500
Subject: [PATCH 05/19] bond/react:add molecule keyword

---
 src/USER-REACTION/fix_bond_react.cpp | 19 +++++++++++++++++++
 src/USER-REACTION/fix_bond_react.h   |  1 +
 2 files changed, 20 insertions(+)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index 994c8cba91..c383c709dd 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -85,6 +85,9 @@ enum{DISTANCE,ANGLE,DIHEDRAL,ARRHENIUS,RMSD};
 // keyword values that accept variables as input
 enum{NEVERY,RMIN,RMAX,PROB};
 
+// values for molecule_keyword
+enum{OFF,INTER,INTRA};
+
 /* ---------------------------------------------------------------------- */
 
 FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
@@ -202,6 +205,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
   memory->create(limit_duration,nreacts,"bond/react:limit_duration");
   memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag");
   memory->create(custom_charges_fragid,nreacts,"bond/react:custom_charges_fragid");
+  memory->create(molecule_keyword,nreacts,"bond/react:molecule_keyword");
   memory->create(constraints,1,MAXCONARGS,"bond/react:constraints");
   memory->create(var_flag,NUMVARVALS,nreacts,"bond/react:var_flag");
   memory->create(var_id,NUMVARVALS,nreacts,"bond/react:var_id");
@@ -222,6 +226,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
     max_rxn[i] = INT_MAX;
     stabilize_steps_flag[i] = 0;
     custom_charges_fragid[i] = -1;
+    molecule_keyword[i] = OFF;
     // set default limit duration to 60 timesteps
     limit_duration[i] = 60;
     reaction_count[i] = 0;
@@ -377,6 +382,14 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
                                                          "'custom_charges' keyword does not exist");
         }
         iarg += 2;
+      } else if (strcmp(arg[iarg],"molecule") == 0) {
+        if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: "
+                                      "'molecule' has too few arguments");
+        if (strcmp(arg[iarg+1],"off") == 0) molecule_keyword[rxn] = OFF; //default
+        else if (strcmp(arg[iarg+1],"inter") == 0) molecule_keyword[rxn] = INTER;
+        else if (strcmp(arg[iarg+1],"intra") == 0) molecule_keyword[rxn] = INTRA;
+        else error->one(FLERR,"Bond/react: Illegal option for 'molecule' keyword");
+        iarg += 2;
       } else error->all(FLERR,"Illegal fix bond/react command: unknown keyword");
     }
   }
@@ -556,6 +569,7 @@ FixBondReact::~FixBondReact()
   memory->destroy(var_id);
   memory->destroy(stabilize_steps_flag);
   memory->destroy(custom_charges_fragid);
+  memory->destroy(molecule_keyword);
 
   memory->destroy(iatomtype);
   memory->destroy(jatomtype);
@@ -1055,6 +1069,11 @@ void FixBondReact::far_partner()
         continue;
       }
 
+      if (molecule_keyword[rxnID] == INTER)
+        if (atom->molecule[i] == atom->molecule[j]) continue;
+      else if (molecule_keyword[rxnID] == INTRA)
+        if (atom->molecule[i] != atom->molecule[j]) continue;
+
       jtype = type[j];
       possible = 0;
       if (itype == iatomtype[rxnID] && jtype == jatomtype[rxnID]) {
diff --git a/src/USER-REACTION/fix_bond_react.h b/src/USER-REACTION/fix_bond_react.h
index 895a63b7ae..6105859136 100644
--- a/src/USER-REACTION/fix_bond_react.h
+++ b/src/USER-REACTION/fix_bond_react.h
@@ -65,6 +65,7 @@ class FixBondReact : public Fix {
   int custom_exclude_flag;
   int *stabilize_steps_flag;
   int *custom_charges_fragid;
+  int *molecule_keyword;
   int nconstraints;
   int narrhenius;
   double **constraints;

From 958ab461b36c5f84c5bf627e737342ee17b4c349 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Mon, 9 Nov 2020 14:06:17 -0500
Subject: [PATCH 06/19] bond/react:molecule keyword docs

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

diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst
index 8a23149aa9..1a26f7dea4 100644
--- a/doc/src/fix_bond_react.rst
+++ b/doc/src/fix_bond_react.rst
@@ -55,6 +55,10 @@ Syntax
          *custom_charges* value = *no* or *fragmentID*
            no = update all atomic charges (default)
            fragmentID = ID of molecule fragment whose charges are updated
+         *molecule* value = *off* or *inter* or *intra*
+           off = allow both inter- and intramolecular reactions (default)
+           inter = search for reactions between molecules with different IDs
+           intra = search for reactions within the same molecule
 
 Examples
 """"""""
@@ -491,12 +495,20 @@ situations, decreasing rather than increasing this parameter will
 result in an increase in stability.
 
 The *custom_charges* keyword can be used to specify which atoms'
-atomic charges are updated. When the value is set to 'no,' all atomic
+atomic charges are updated. When the value is set to 'no', all atomic
 charges are updated to those specified by the post-reaction template
 (default). Otherwise, the value should be the name of a molecule
 fragment defined in the pre-reaction molecule template. In this case,
 only the atomic charges of atoms in the molecule fragment are updated.
 
+The *molecule* keyword can be used to force the reaction to be
+intermolecular, intramolecular or either. When the value is set to
+'off', molecule IDs are not considered when searching for reactions
+(default). When the value is set to 'inter', the bonding atoms must
+have different molecule IDs in order to be considered for the
+reaction. When the value is set to 'intra', only atoms with the same
+molecule ID are considered for the reaction.
+
 A few other considerations:
 
 Many reactions result in one or more atoms that are considered
@@ -558,7 +570,7 @@ These is 1 quantity for each react argument:
 No parameter of this fix can be used with the *start/stop* keywords
 of the :doc:`run <run>` command.  This fix is not invoked during :doc:`energy minimization <minimize>`.
 
-When fix bond/react is 'unfixed,' all internally-created groups are
+When fix bond/react is 'unfixed', all internally-created groups are
 deleted. Therefore, fix bond/react can only be unfixed after unfixing
 all other fixes that use any group created by fix bond/react.
 
@@ -581,7 +593,7 @@ Default
 """""""
 
 The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60,
-reset_mol_ids = yes, custom_charges = no
+reset_mol_ids = yes, custom_charges = no, molecule = off
 
 ----------
 

From f6975bf4eb5a7ac5347aad117a0164532ad2b422 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Mon, 9 Nov 2020 14:10:41 -0500
Subject: [PATCH 07/19] also add molecule_keyword to close_partner, for
 consistency

---
 src/USER-REACTION/fix_bond_react.cpp | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index c383c709dd..12e44c6a9e 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -1153,6 +1153,11 @@ void FixBondReact::close_partner()
       if (i_limit_tags[i2] != 0) continue;
       if (itype != iatomtype[rxnID] || jtype != jatomtype[rxnID]) continue;
 
+      if (molecule_keyword[rxnID] == INTER)
+        if (atom->molecule[i1] == atom->molecule[i2]) continue;
+      else if (molecule_keyword[rxnID] == INTRA)
+        if (atom->molecule[i1] != atom->molecule[i2]) continue;
+
       delx = x[i1][0] - x[i2][0];
       dely = x[i1][1] - x[i2][1];
       delz = x[i1][2] - x[i2][2];

From 6ece19f91940d4de51566b48e7cd59a64d28cba0 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Mon, 9 Nov 2020 16:29:52 -0500
Subject: [PATCH 08/19] Update fix_bond_react.rst

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

diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst
index 1a26f7dea4..1139d4e9fa 100644
--- a/doc/src/fix_bond_react.rst
+++ b/doc/src/fix_bond_react.rst
@@ -506,8 +506,8 @@ intermolecular, intramolecular or either. When the value is set to
 'off', molecule IDs are not considered when searching for reactions
 (default). When the value is set to 'inter', the bonding atoms must
 have different molecule IDs in order to be considered for the
-reaction. When the value is set to 'intra', only atoms with the same
-molecule ID are considered for the reaction.
+reaction. When the value is set to 'intra', only bonding atoms with
+the same molecule ID are considered for the reaction.
 
 A few other considerations:
 

From 64cc0adb9e7072fb302ba011b3635e1640809718 Mon Sep 17 00:00:00 2001
From: jrgissing <jrgiss05@gmail.com>
Date: Wed, 11 Nov 2020 13:27:34 -0500
Subject: [PATCH 09/19] bond/react: add new reference

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

diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst
index 1139d4e9fa..92005bd2b2 100644
--- a/doc/src/fix_bond_react.rst
+++ b/doc/src/fix_bond_react.rst
@@ -600,3 +600,7 @@ reset_mol_ids = yes, custom_charges = no, molecule = off
 .. _Gissinger:
 
 **(Gissinger)** Gissinger, Jensen and Wise, Polymer, 128, 211 (2017).
+
+.. _Gissinger:
+
+**(Gissinger)** Gissinger, Jensen and Wise, Macromolecules (2020, in press).

From c013e6db10b21b8596e294c344e949bd6b0338a7 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Fri, 13 Nov 2020 00:50:01 -0500
Subject: [PATCH 10/19] revert cutoff check

add temporary fix for hybrid pair_style
---
 src/USER-REACTION/fix_bond_react.cpp | 15 +++++++--------
 1 file changed, 7 insertions(+), 8 deletions(-)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index 12e44c6a9e..bda2b1c148 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -472,14 +472,6 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
     }
   }
 
-  for (int i = 0; i < nreacts; i++) {
-    if (closeneigh[i] == -1) { // indicates will search non-bonded neighbors
-      if (cutsq[i][1] > neighbor->cutneighsq[iatomtype[i]][jatomtype[i]]) {
-        error->all(FLERR,"Bond/react: Fix bond/react cutoff is longer than pairwise cutoff");
-      }
-    }
-  }
-
   // initialize Marsaglia RNG with processor-unique seed ('prob' keyword)
 
   random = new class RanMars*[nreacts];
@@ -772,6 +764,13 @@ void FixBondReact::init()
   if (strstr(update->integrate_style,"respa"))
     nlevels_respa = ((Respa *) update->integrate)->nlevels;
 
+    // check cutoff for iatomtype,jatomtype
+    for (int i = 0; i < nreacts; i++) {
+      if (!utils::strmatch(force->pair_style,"^hybrid"))
+        if (force->pair == nullptr || cutsq[i][1] > force->pair->cutsq[iatomtype[i]][jatomtype[i]])
+          error->all(FLERR,"Bond/react: Fix bond/react cutoff is longer than pairwise cutoff");
+    }
+
   // need a half neighbor list, built every Nevery steps
   int irequest = neighbor->request(this,instance_me);
   neighbor->requests[irequest]->pair = 0;

From 348afb68676a22d7bdfdb09c088c799df74cbb18 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 14 Nov 2020 17:27:49 -0500
Subject: [PATCH 11/19] correct cutsq check in proper place

---
 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 bda2b1c148..9f880b10dd 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -766,7 +766,7 @@ void FixBondReact::init()
 
     // check cutoff for iatomtype,jatomtype
     for (int i = 0; i < nreacts; i++) {
-      if (!utils::strmatch(force->pair_style,"^hybrid"))
+      if (closeneigh[i] == -1) // indicates will search for non-bonded bonding atoms
         if (force->pair == nullptr || cutsq[i][1] > force->pair->cutsq[iatomtype[i]][jatomtype[i]])
           error->all(FLERR,"Bond/react: Fix bond/react cutoff is longer than pairwise cutoff");
     }

From aa129bc2188fa84ab61e5a1a28637d053c83ea86 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Thu, 19 Nov 2020 13:29:55 -0500
Subject: [PATCH 12/19] bond/react: update doc terminology

change 'bonding atoms' to 'initiator atoms'
(bonds need not form or break between initiator atoms)
---
 doc/src/fix_bond_react.rst | 100 ++++++++++++++++++-------------------
 1 file changed, 50 insertions(+), 50 deletions(-)

diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst
index 92005bd2b2..3b98f44713 100644
--- a/doc/src/fix_bond_react.rst
+++ b/doc/src/fix_bond_react.rst
@@ -35,8 +35,8 @@ Syntax
 * react-ID = user-assigned name for the reaction
 * react-group-ID = only atoms in this group are considered for the reaction
 * Nevery = attempt reaction every this many steps
-* Rmin = bonding pair atoms must be separated by more than Rmin to initiate reaction (distance units)
-* Rmax = bonding pair atoms must be separated by less than Rmax to initiate reaction (distance units)
+* Rmin = initiator atoms must be separated by more than Rmin to initiate reaction (distance units)
+* Rmax = initiator atoms must be separated by less than Rmax to initiate reaction (distance units)
 * template-ID(pre-reacted) = ID of a molecule template containing pre-reaction topology
 * template-ID(post-reacted) = ID of a molecule template containing post-reaction topology
 * map_file = name of file specifying corresponding atom-IDs in the pre- and post-reacted templates
@@ -176,8 +176,8 @@ timesteps. *Nevery* can be specified with an equal-style
 integer.
 
 Three physical conditions must be met for a reaction to occur. First,
-a bonding atom pair must be identified within the reaction distance
-cutoffs. Second, the topology surrounding the bonding atom pair must
+an initiator atom pair must be identified within the reaction distance
+cutoffs. Second, the topology surrounding the initiator atom pair must
 match the topology of the pre-reaction template. Only atom types and
 bond connectivity are used to identify a valid reaction site (not bond
 types, etc.). Finally, any reaction constraints listed in the map file
@@ -185,10 +185,10 @@ types, etc.). Finally, any reaction constraints listed in the map file
 reaction site is eligible to be modified to match the post-reaction
 template.
 
-A bonding atom pair will be identified if several conditions are met.
-First, a pair of atoms I,J within the specified react-group-ID of type
-itype and jtype must be separated by a distance between *Rmin* and
-*Rmax*\ . *Rmin* and *Rmax* can be specified with equal-style
+An initiator atom pair will be identified if several conditions are
+met. First, a pair of atoms I,J within the specified react-group-ID of
+type itype and jtype must be separated by a distance between *Rmin*
+and *Rmax*\ . *Rmin* and *Rmax* can be specified with equal-style
 :doc:`variables <variable>`. For example, these reaction cutoffs can
 be a function of the reaction conversion using the following commands:
 
@@ -198,20 +198,20 @@ be a function of the reaction conversion using the following commands:
    fix myrxn all bond/react react myrxn1 all 1 0 v_rmax mol1 mol2 map_file.txt
    variable rmax equal 3+f_myrxn[1]/100 # arbitrary function of reaction count
 
-The following criteria are used if multiple candidate bonding atom
-pairs are identified within the cutoff distance: 1) If the bonding
+The following criteria are used if multiple candidate initiator atom
+pairs are identified within the cutoff distance: 1) If the initiator
 atoms in the pre-reaction template are not 1-2 neighbors (i.e. not
 directly bonded) the closest potential partner is chosen. 2)
-Otherwise, if the bonding atoms in the pre-reaction template are 1-2
+Otherwise, if the initiator atoms in the pre-reaction template are 1-2
 neighbors (i.e. directly bonded) the farthest potential partner is
 chosen. 3) Then, if both an atom I and atom J have each other as their
-bonding partners, these two atoms are identified as the bonding atom
-pair of the reaction site. Note that it can be helpful to select
-unique atom types for the bonding atoms: if a bonding atom pair is
-identified, as described in the previous steps, but does not
+initiator partners, these two atoms are identified as the initiator
+atom pair of the reaction site. Note that it can be helpful to select
+unique atom types for the initiator atoms: if an initiator atom pair
+is identified, as described in the previous steps, but does not
 correspond to the same pair specified in the pre-reaction template, an
 otherwise eligible reaction could be prevented from occurring. Once
-this unique bonding atom pair is identified for each reaction, there
+this unique initiator atom pair is identified for each reaction, there
 could be two or more reactions that involve the same atom on the same
 timestep. If this is the case, only one such reaction is permitted to
 occur. This reaction is chosen randomly from all potential reactions
@@ -221,7 +221,7 @@ with user-specified probabilities.
 
 The pre-reacted molecule template is specified by a molecule command.
 This molecule template file contains a sample reaction site and its
-surrounding topology. As described below, the bonding atom pairs of
+surrounding topology. As described below, the initiator 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
@@ -235,18 +235,18 @@ missing topology with respect to the simulation. For example, the
 pre-reacted template may contain an atom that, in the simulation, is
 currently connected to the rest of a long polymer chain. These are
 referred to as edge atoms, and are also specified in the map file. All
-pre-reaction template atoms should be linked to a bonding atom, via at
-least one path that does not involve edge atoms. When the pre-reaction
-template contains edge atoms, not all atoms, bonds, charges, etc.
-specified in the reaction templates will be updated. Specifically,
-topology that involves only atoms that are 'too near' to template
-edges will not be updated. The definition of 'too near the edge'
-depends on which interactions are defined in the simulation. If the
-simulation has defined dihedrals, atoms within two bonds of edge atoms
-are considered 'too near the edge.' If the simulation defines angles,
-but not dihedrals, atoms within one bond of edge atoms are considered
-'too near the edge.' If just bonds are defined, only edge atoms are
-considered 'too near the edge.'
+pre-reaction template atoms should be linked to an initiator atom, via
+at least one path that does not involve edge atoms. When the
+pre-reaction template contains edge atoms, not all atoms, bonds,
+charges, etc. specified in the reaction templates will be updated.
+Specifically, topology that involves only atoms that are 'too near' to
+template edges will not be updated. The definition of 'too near the
+edge' depends on which interactions are defined in the simulation. If
+the simulation has defined dihedrals, atoms within two bonds of edge
+atoms are considered 'too near the edge.' If the simulation defines
+angles, but not dihedrals, atoms within one bond of edge atoms are
+considered 'too near the edge.' If just bonds are defined, only edge
+atoms are considered 'too near the edge.'
 
 .. note::
 
@@ -302,23 +302,23 @@ The optional keywords are 'edgeIDs', 'deleteIDs', 'chiralIDs' and
 
 The body of the map file contains two mandatory sections and four
 optional sections. The first mandatory section begins with the keyword
-'BondingIDs' and lists the atom IDs of the bonding atom pair in the
-pre-reacted molecule template. The second mandatory section begins
-with the keyword 'Equivalences' and lists a one-to-one correspondence
-between atom IDs of the pre- and post-reacted templates. The first
-column is an atom ID of the pre-reacted molecule template, and the
-second column is the corresponding atom ID of the post-reacted
-molecule template. The first optional section begins with the keyword
-'EdgeIDs' and lists the atom IDs of edge atoms in the pre-reacted
-molecule template. The second optional section begins with the keyword
-'DeleteIDs' and lists the atom IDs of pre-reaction template atoms to
-delete. The third optional section begins with the keyword 'ChiralIDs'
-lists the atom IDs of chiral atoms whose handedness should be
-enforced. The fourth optional section begins with the keyword
-'Constraints' and lists additional criteria that must be satisfied in
-order for the reaction to occur. Currently, there are five types of
-constraints available, as discussed below: 'distance', 'angle',
-'dihedral', 'arrhenius', and 'rmsd'.
+'InitiatorIDs' and lists the two atom IDs of the initiator atom pair
+in the pre-reacted molecule template. The second mandatory section
+begins with the keyword 'Equivalences' and lists a one-to-one
+correspondence between atom IDs of the pre- and post-reacted
+templates. The first column is an atom ID of the pre-reacted molecule
+template, and the second column is the corresponding atom ID of the
+post-reacted molecule template. The first optional section begins with
+the keyword 'EdgeIDs' and lists the atom IDs of edge atoms in the
+pre-reacted molecule template. The second optional section begins with
+the keyword 'DeleteIDs' and lists the atom IDs of pre-reaction
+template atoms to delete. The third optional section begins with the
+keyword 'ChiralIDs' lists the atom IDs of chiral atoms whose
+handedness should be enforced. The fourth optional section begins with
+the keyword 'Constraints' and lists additional criteria that must be
+satisfied in order for the reaction to occur. Currently, there are
+five types of constraints available, as discussed below: 'distance',
+'angle', 'dihedral', 'arrhenius', and 'rmsd'.
 
 A sample map file is given below:
 
@@ -331,7 +331,7 @@ A sample map file is given below:
    7 equivalences
    2 edgeIDs
 
-   BondingIDs
+   InitiatorIDs
 
    3
    5
@@ -466,7 +466,7 @@ A few capabilities to note: 1) You may specify as many *react*
 arguments as desired. For example, you could break down a complicated
 reaction mechanism into several reaction steps, each defined by its
 own *react* argument. 2) While typically a bond is formed or removed
-between the bonding atom pairs specified in the pre-reacted molecule
+between the initiator atoms specified in the pre-reacted molecule
 template, this is not required. 3) By reversing the order of the pre-
 and post- reacted molecule templates in another *react* argument, you
 can allow for the possibility of one or more reverse reactions.
@@ -504,9 +504,9 @@ only the atomic charges of atoms in the molecule fragment are updated.
 The *molecule* keyword can be used to force the reaction to be
 intermolecular, intramolecular or either. When the value is set to
 'off', molecule IDs are not considered when searching for reactions
-(default). When the value is set to 'inter', the bonding atoms must
+(default). When the value is set to 'inter', the initiator atoms must
 have different molecule IDs in order to be considered for the
-reaction. When the value is set to 'intra', only bonding atoms with
+reaction. When the value is set to 'intra', only initiator atoms with
 the same molecule ID are considered for the reaction.
 
 A few other considerations:

From 5aff81946b305db2422942b48132356e5a359d9a Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Thu, 19 Nov 2020 13:48:11 -0500
Subject: [PATCH 13/19] bond/react: update examples with new keyword

---
 examples/USER/reaction/nylon,6-6_melt/rxn1_stp1_map            | 2 +-
 examples/USER/reaction/nylon,6-6_melt/rxn1_stp2_map            | 2 +-
 examples/USER/reaction/tiny_epoxy/rxn1_stp1.map                | 2 +-
 examples/USER/reaction/tiny_epoxy/rxn1_stp2.map                | 2 +-
 examples/USER/reaction/tiny_epoxy/rxn2_stp1.map                | 2 +-
 examples/USER/reaction/tiny_epoxy/rxn2_stp2.map                | 2 +-
 examples/USER/reaction/tiny_nylon/rxn1_stp1_map                | 2 +-
 examples/USER/reaction/tiny_nylon/rxn1_stp2_map                | 2 +-
 examples/USER/reaction/tiny_polystyrene/2styrene_map           | 2 +-
 examples/USER/reaction/tiny_polystyrene/chain_chain_map        | 2 +-
 examples/USER/reaction/tiny_polystyrene/chain_plus_styrene_map | 2 +-
 11 files changed, 11 insertions(+), 11 deletions(-)

diff --git a/examples/USER/reaction/nylon,6-6_melt/rxn1_stp1_map b/examples/USER/reaction/nylon,6-6_melt/rxn1_stp1_map
index 44f7ad8137..0de209fdfd 100644
--- a/examples/USER/reaction/nylon,6-6_melt/rxn1_stp1_map
+++ b/examples/USER/reaction/nylon,6-6_melt/rxn1_stp1_map
@@ -3,7 +3,7 @@ this is a nominal superimpose file
 2 edgeIDs
 18 equivalences
 
-BondingIDs
+InitiatorIDs
 
 10
 1
diff --git a/examples/USER/reaction/nylon,6-6_melt/rxn1_stp2_map b/examples/USER/reaction/nylon,6-6_melt/rxn1_stp2_map
index 35fe47fdb3..8389688dd5 100644
--- a/examples/USER/reaction/nylon,6-6_melt/rxn1_stp2_map
+++ b/examples/USER/reaction/nylon,6-6_melt/rxn1_stp2_map
@@ -3,7 +3,7 @@ this is a nominal superimpose file
 2 edgeIDs
 15 equivalences
 
-BondingIDs
+InitiatorIDs
 
 4
 12
diff --git a/examples/USER/reaction/tiny_epoxy/rxn1_stp1.map b/examples/USER/reaction/tiny_epoxy/rxn1_stp1.map
index fb8720be63..9901104f03 100644
--- a/examples/USER/reaction/tiny_epoxy/rxn1_stp1.map
+++ b/examples/USER/reaction/tiny_epoxy/rxn1_stp1.map
@@ -3,7 +3,7 @@ this is a map file
 1 edgeIDs
 31 equivalences
 
-BondingIDs
+InitiatorIDs
 
 15
 1
diff --git a/examples/USER/reaction/tiny_epoxy/rxn1_stp2.map b/examples/USER/reaction/tiny_epoxy/rxn1_stp2.map
index 9e39e57310..b83dbaf6c9 100644
--- a/examples/USER/reaction/tiny_epoxy/rxn1_stp2.map
+++ b/examples/USER/reaction/tiny_epoxy/rxn1_stp2.map
@@ -3,7 +3,7 @@ this is a map file
 1 edgeIDs
 31 equivalences
 
-BondingIDs
+InitiatorIDs
 
 4
 25
diff --git a/examples/USER/reaction/tiny_epoxy/rxn2_stp1.map b/examples/USER/reaction/tiny_epoxy/rxn2_stp1.map
index cbe2a44f36..46b916c2ee 100644
--- a/examples/USER/reaction/tiny_epoxy/rxn2_stp1.map
+++ b/examples/USER/reaction/tiny_epoxy/rxn2_stp1.map
@@ -3,7 +3,7 @@ this is a map file
 2 edgeIDs
 42 equivalences
 
-BondingIDs
+InitiatorIDs
 
 15
 32
diff --git a/examples/USER/reaction/tiny_epoxy/rxn2_stp2.map b/examples/USER/reaction/tiny_epoxy/rxn2_stp2.map
index f90915a54f..0d23a884e7 100644
--- a/examples/USER/reaction/tiny_epoxy/rxn2_stp2.map
+++ b/examples/USER/reaction/tiny_epoxy/rxn2_stp2.map
@@ -3,7 +3,7 @@ this is a map file
 2 edgeIDs
 42 equivalences
 
-BondingIDs
+InitiatorIDs
 
 35
 24
diff --git a/examples/USER/reaction/tiny_nylon/rxn1_stp1_map b/examples/USER/reaction/tiny_nylon/rxn1_stp1_map
index 44f7ad8137..0de209fdfd 100644
--- a/examples/USER/reaction/tiny_nylon/rxn1_stp1_map
+++ b/examples/USER/reaction/tiny_nylon/rxn1_stp1_map
@@ -3,7 +3,7 @@ this is a nominal superimpose file
 2 edgeIDs
 18 equivalences
 
-BondingIDs
+InitiatorIDs
 
 10
 1
diff --git a/examples/USER/reaction/tiny_nylon/rxn1_stp2_map b/examples/USER/reaction/tiny_nylon/rxn1_stp2_map
index 35fe47fdb3..8389688dd5 100644
--- a/examples/USER/reaction/tiny_nylon/rxn1_stp2_map
+++ b/examples/USER/reaction/tiny_nylon/rxn1_stp2_map
@@ -3,7 +3,7 @@ this is a nominal superimpose file
 2 edgeIDs
 15 equivalences
 
-BondingIDs
+InitiatorIDs
 
 4
 12
diff --git a/examples/USER/reaction/tiny_polystyrene/2styrene_map b/examples/USER/reaction/tiny_polystyrene/2styrene_map
index be2a789671..e48bbac724 100644
--- a/examples/USER/reaction/tiny_polystyrene/2styrene_map
+++ b/examples/USER/reaction/tiny_polystyrene/2styrene_map
@@ -3,7 +3,7 @@ this is a map file
 0 edgeIDs
 32 equivalences
 
-BondingIDs
+InitiatorIDs
 
 4
 30
diff --git a/examples/USER/reaction/tiny_polystyrene/chain_chain_map b/examples/USER/reaction/tiny_polystyrene/chain_chain_map
index 5a0fde6742..c4813238ef 100644
--- a/examples/USER/reaction/tiny_polystyrene/chain_chain_map
+++ b/examples/USER/reaction/tiny_polystyrene/chain_chain_map
@@ -8,7 +8,7 @@ EdgeIDs
 17
 34
 
-BondingIDs
+InitiatorIDs
 
 14
 38
diff --git a/examples/USER/reaction/tiny_polystyrene/chain_plus_styrene_map b/examples/USER/reaction/tiny_polystyrene/chain_plus_styrene_map
index 56808b0962..7e61e49c66 100644
--- a/examples/USER/reaction/tiny_polystyrene/chain_plus_styrene_map
+++ b/examples/USER/reaction/tiny_polystyrene/chain_plus_styrene_map
@@ -7,7 +7,7 @@ EdgeIDs
 
 30
 
-BondingIDs
+InitiatorIDs
 
 14
 34

From 2a765efc8c0f17f766dc475c64e58746db9592e9 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Thu, 19 Nov 2020 13:50:30 -0500
Subject: [PATCH 14/19] bond/react: update map file section

retain old section title for backwards compatibility
---
 src/USER-REACTION/fix_bond_react.cpp | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index 9f880b10dd..3fd76e3779 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -3228,7 +3228,9 @@ void FixBondReact::read(int myrxn)
 
   int equivflag = 0, bondflag = 0;
   while (strlen(keyword)) {
-    if (strcmp(keyword,"BondingIDs") == 0) {
+    if (strcmp(keyword,"InitiatorIDs") == 0 || strcmp(keyword,"BondingIDs") == 0) {
+      if (strcmp(keyword,"BondingIDs") == 0)
+        error->warning(FLERR,"Bond/react: The BondingIDs section title has been deprecated. Please use InitiatorIDs instead.");
       bondflag = 1;
       readline(line);
       sscanf(line,"%d",&ibonding[myrxn]);
@@ -3257,7 +3259,7 @@ void FixBondReact::read(int myrxn)
 
   // error check
   if (bondflag == 0 || equivflag == 0)
-    error->all(FLERR,"Bond/react: Map file missing BondingIDs or Equivalences section\n");
+    error->all(FLERR,"Bond/react: Map file missing InitiatorIDs or Equivalences section\n");
 }
 
 void FixBondReact::EdgeIDs(char *line, int myrxn)

From dba84be75a96874916c459991f51cd43b02957c9 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Thu, 19 Nov 2020 13:52:58 -0500
Subject: [PATCH 15/19] Update false_positives.txt

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

diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt
index 3ef0b904eb..40552f4198 100644
--- a/doc/utils/sphinx-config/false_positives.txt
+++ b/doc/utils/sphinx-config/false_positives.txt
@@ -277,7 +277,6 @@ Boltzman
 BondAngle
 BondBond
 bondchk
-BondingIDs
 bondmax
 bondtype
 Bonet
@@ -1324,6 +1323,7 @@ inhomogeneous
 init
 initialdelay
 initializations
+InitiatorIDs
 initio
 InP
 inregion

From f42e907d30596b89b3553f17451e435c5f5b287d Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Thu, 19 Nov 2020 19:06:33 -0500
Subject: [PATCH 16/19] update terminology

---
 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 3fd76e3779..810adf4b48 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -1330,7 +1330,7 @@ void FixBondReact::superimpose_algorithm()
         //if (hang_catch > onemol->natoms*4)
         if (hang_catch > atom->nlocal*30) {
           error->one(FLERR,"Bond/react: Excessive iteration of superimpose algorithm. "
-              "Please check that all pre-reaction template atoms are linked to a bonding atom, "
+              "Please check that all pre-reaction template atoms are linked to an initiator atom, "
               "via at least one path that does not involve edge atoms.");
         }
       }

From fcf5f34a4984ad14a09756808bd36f1ef21224cc Mon Sep 17 00:00:00 2001
From: jrgissing <jrgiss05@gmail.com>
Date: Fri, 27 Nov 2020 12:08:09 -0500
Subject: [PATCH 17/19] memory fix

---
 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 810adf4b48..937f22c857 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -1228,7 +1228,7 @@ void FixBondReact::superimpose_algorithm()
   memory->create(glove,max_natoms,2,"bond/react:glove");
   memory->create(restore_pt,MAXGUESS,4,"bond/react:restore_pt");
   memory->create(pioneers,max_natoms,"bond/react:pioneers");
-  memory->create(restore,max_natoms,MAXGUESS,"bond/react:restore");
+  memory->create(restore,max_natoms,MAXGUESS*4,"bond/react:restore");
   memory->create(local_mega_glove,max_natoms+1,allncreate,"bond/react:local_mega_glove");
   memory->create(ghostly_mega_glove,max_natoms+1,allncreate,"bond/react:ghostly_mega_glove");
 

From 8df0ee0dfaaf315f8b3d0942b083fed08f7f8ece Mon Sep 17 00:00:00 2001
From: jrgissing <jrgiss05@gmail.com>
Date: Fri, 27 Nov 2020 12:15:31 -0500
Subject: [PATCH 18/19] make new warning more gentle

---
 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 937f22c857..29a99db75f 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -3230,7 +3230,7 @@ void FixBondReact::read(int myrxn)
   while (strlen(keyword)) {
     if (strcmp(keyword,"InitiatorIDs") == 0 || strcmp(keyword,"BondingIDs") == 0) {
       if (strcmp(keyword,"BondingIDs") == 0)
-        error->warning(FLERR,"Bond/react: The BondingIDs section title has been deprecated. Please use InitiatorIDs instead.");
+        if (me == 0) error->warning(FLERR,"Bond/react: The BondingIDs section title has been deprecated. Please use InitiatorIDs instead.");
       bondflag = 1;
       readline(line);
       sscanf(line,"%d",&ibonding[myrxn]);

From ce3c92d52e150dd6ff323c7165e77124b3ff7de4 Mon Sep 17 00:00:00 2001
From: Richard Berger <richard.berger@temple.edu>
Date: Mon, 30 Nov 2020 14:55:35 -0500
Subject: [PATCH 19/19] Avoid duplicate doc target

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

diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst
index 3b98f44713..16bd4ecb71 100644
--- a/doc/src/fix_bond_react.rst
+++ b/doc/src/fix_bond_react.rst
@@ -601,6 +601,6 @@ reset_mol_ids = yes, custom_charges = no, molecule = off
 
 **(Gissinger)** Gissinger, Jensen and Wise, Polymer, 128, 211 (2017).
 
-.. _Gissinger:
+.. _Gissinger2020:
 
 **(Gissinger)** Gissinger, Jensen and Wise, Macromolecules (2020, in press).