From a64a9c12535fd2d4456d0a19bab5c5c231a0a26f Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 23 May 2020 22:13:20 -0600
Subject: [PATCH 01/11] distance constraint: fragment support

---
 src/USER-REACTION/fix_bond_react.cpp | 83 +++++++++++++++++++++++-----
 src/USER-REACTION/fix_bond_react.h   | 18 ++++--
 2 files changed, 80 insertions(+), 21 deletions(-)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index 39750dd7da..6f57a6182f 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -1757,6 +1757,7 @@ evaluate constraints: return 0 if any aren't satisfied
 int FixBondReact::check_constraints()
 {
   tagint atom1,atom2,atom3,atom4;
+  double x1[3],x2[3];
   double delx,dely,delz,rsq;
   double delx1,dely1,delz1,delx2,dely2,delz2;
   double rsq1,rsq2,r1,r2,c,t,prrhob;
@@ -1771,14 +1772,14 @@ int FixBondReact::check_constraints()
   for (int i = 0; i < nconstraints; i++) {
     if (constraints[i][0] == rxnID) {
       if (constraints[i][1] == DISTANCE) {
-        atom1 = atom->map(glove[(int) constraints[i][2]-1][1]);
-        atom2 = atom->map(glove[(int) constraints[i][3]-1][1]);
-        delx = x[atom1][0] - x[atom2][0];
-        dely = x[atom1][1] - x[atom2][1];
-        delz = x[atom1][2] - x[atom2][2];
+        get_IDcoords((int) constraints[i][2], (int) constraints[i][3], x1);
+        get_IDcoords((int) constraints[i][4], (int) constraints[i][5], x2);
+        delx = x1[0] - x2[0];
+        dely = x1[1] - x2[1];
+        delz = x1[2] - x2[2];
         domain->minimum_image(delx,dely,delz); // ghost location fix
         rsq = delx*delx + dely*dely + delz*delz;
-        if (rsq < constraints[i][4] || rsq > constraints[i][5]) return 0;
+        if (rsq < constraints[i][6] || rsq > constraints[i][7]) return 0;
       } else if (constraints[i][1] == ANGLE) {
         atom1 = atom->map(glove[(int) constraints[i][2]-1][1]);
         atom2 = atom->map(glove[(int) constraints[i][3]-1][1]);
@@ -1904,6 +1905,37 @@ int FixBondReact::check_constraints()
   return 1;
 }
 
+/* ----------------------------------------------------------------------
+return pre-reaction atom or fragment location
+fragment: given pre-reacted molID (onemol) and fragID,
+          return geometric center (of mapped simulation atoms)
+------------------------------------------------------------------------- */
+
+void FixBondReact::get_IDcoords(int mode, int myID, double *center)
+{
+  double **x = atom->x;
+  if (mode == 1) {
+    int iatom = atom->map(glove[myID-1][1]);
+    for (int i = 0; i < 3; i++)
+      center[i] = x[iatom][i];
+  } else {
+    int nfragatoms = 0;
+    for (int i = 0; i < 3; i++)
+      center[i] = 0;
+
+    for (int i = 0; i < onemol->natoms; i++) {
+      if (onemol->fragmentmask[myID][i]) {
+        for (int j = 0; j < 3; j++) {
+          center[j] += x[atom->map(glove[i][1])][j];
+        }
+        nfragatoms++;
+      }
+    }
+    for (int i = 0; i < 3; i++)
+      center[i] /= nfragatoms;
+  }
+}
+
 /* ----------------------------------------------------------------------
 compute local temperature: average over all atoms in reaction template
 ------------------------------------------------------------------------- */
@@ -3255,21 +3287,21 @@ void FixBondReact::ChiralCenters(char *line, int myrxn)
 void FixBondReact::Constraints(char *line, int myrxn)
 {
   double tmp[MAXCONARGS];
-  int n = strlen("distance") + 1;
-  char *constraint_type = new char[n];
+  char **strargs;
+  memory->create(strargs,MAXCONARGS,MAXLINE,"bond/react:strargs");
+  char *constraint_type = new char[MAXLINE];
   for (int i = 0; i < nconstr; i++) {
     readline(line);
     sscanf(line,"%s",constraint_type);
     constraints[nconstraints][0] = myrxn;
     if (strcmp(constraint_type,"distance") == 0) {
       constraints[nconstraints][1] = DISTANCE;
-      sscanf(line,"%*s %lg %lg %lg %lg",&tmp[0],&tmp[1],&tmp[2],&tmp[3]);
-      if (tmp[0] > onemol->natoms || tmp[1] > onemol->natoms)
-        error->one(FLERR,"Bond/react: Invalid template atom ID in map file");
-      constraints[nconstraints][2] = tmp[0];
-      constraints[nconstraints][3] = tmp[1];
-      constraints[nconstraints][4] = tmp[2]*tmp[2]; // using square of distance
-      constraints[nconstraints][5] = tmp[3]*tmp[3];
+      sscanf(line,"%*s %s %s %lg %lg",strargs[0],strargs[1],&tmp[0],&tmp[1]);
+      readID(strargs[0], nconstraints, 2, 3);
+      readID(strargs[1], nconstraints, 4, 5);
+      // cutoffs
+      constraints[nconstraints][6] = tmp[0]*tmp[0]; // using square of distance
+      constraints[nconstraints][7] = tmp[1]*tmp[1];
     } else if (strcmp(constraint_type,"angle") == 0) {
       constraints[nconstraints][1] = ANGLE;
       sscanf(line,"%*s %lg %lg %lg %lg %lg",&tmp[0],&tmp[1],&tmp[2],&tmp[3],&tmp[4]);
@@ -3310,6 +3342,27 @@ void FixBondReact::Constraints(char *line, int myrxn)
     nconstraints++;
   }
   delete [] constraint_type;
+  memory->destroy(strargs);
+}
+
+/* ----------------------------------------------------------------------
+if ID starts with character, assume it is a pre-reaction molecule fragment ID
+otherwise, it is a pre-reaction atom ID
+---------------------------------------------------------------------- */
+
+void FixBondReact::readID(char *strarg, int iconstr, int mode, int myID)
+{
+  if (isalpha(strarg[0])) {
+    constraints[iconstr][mode] = 0; // fragment vs. atom ID flag
+    int ifragment = onemol->findfragment(strarg);
+    if (ifragment < 0) error->one(FLERR,"Bond/react: Molecule fragment does not exist");
+    constraints[iconstr][myID] = ifragment;
+  } else {
+    constraints[iconstr][mode] = 1; // fragment vs. atom ID flag
+    int iatom = atoi(strarg);
+    if (iatom > onemol->natoms) error->one(FLERR,"Bond/react: Invalid template atom ID in map file");
+    constraints[iconstr][myID] = iatom;
+  }
 }
 
 void FixBondReact::open(char *file)
diff --git a/src/USER-REACTION/fix_bond_react.h b/src/USER-REACTION/fix_bond_react.h
index 722a533bbf..e8a4253ce7 100644
--- a/src/USER-REACTION/fix_bond_react.h
+++ b/src/USER-REACTION/fix_bond_react.h
@@ -147,12 +147,13 @@ class FixBondReact : public Fix {
   int glove_counter; // used to determine when to terminate Superimpose Algorithm
 
   void read(int);
-  void EdgeIDs(char *,int);
-  void Equivalences(char *,int);
-  void CustomEdges(char *,int);
-  void DeleteAtoms(char *,int);
-  void ChiralCenters(char *,int);
-  void Constraints(char *,int);
+  void EdgeIDs(char *, int);
+  void Equivalences(char *, int);
+  void CustomEdges(char *, int);
+  void DeleteAtoms(char *, int);
+  void ChiralCenters(char *, int);
+  void Constraints(char *, int);
+  void readID(char *, int, int, int);
 
   void make_a_guess ();
   void neighbor_loop();
@@ -161,6 +162,7 @@ class FixBondReact : public Fix {
   void inner_crosscheck_loop();
   void ring_check();
   int check_constraints();
+  void get_IDcoords(int, int, double *);
   double get_temperature();
   int get_chirality(double[12]); // get handedness given an ordered set of coordinates
 
@@ -288,4 +290,8 @@ E: Bond/react: Variable is not equal-style
 
 Self-explanatory.
 
+E: Bond/react: Molecule fragment does not exist
+
+Self-explanatory.
+
 */

From 74d58778b489760cc1b74b7cbf4eaa912bcb3ed9 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 23 May 2020 22:29:14 -0600
Subject: [PATCH 02/11] angle constraint: fragment support

---
 src/USER-REACTION/fix_bond_react.cpp | 36 +++++++++++++---------------
 1 file changed, 17 insertions(+), 19 deletions(-)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index 6f57a6182f..02b06313d4 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -1757,7 +1757,7 @@ evaluate constraints: return 0 if any aren't satisfied
 int FixBondReact::check_constraints()
 {
   tagint atom1,atom2,atom3,atom4;
-  double x1[3],x2[3];
+  double x1[3],x2[3],x3[3];
   double delx,dely,delz,rsq;
   double delx1,dely1,delz1,delx2,dely2,delz2;
   double rsq1,rsq2,r1,r2,c,t,prrhob;
@@ -1781,22 +1781,22 @@ int FixBondReact::check_constraints()
         rsq = delx*delx + dely*dely + delz*delz;
         if (rsq < constraints[i][6] || rsq > constraints[i][7]) return 0;
       } else if (constraints[i][1] == ANGLE) {
-        atom1 = atom->map(glove[(int) constraints[i][2]-1][1]);
-        atom2 = atom->map(glove[(int) constraints[i][3]-1][1]);
-        atom3 = atom->map(glove[(int) constraints[i][4]-1][1]);
+        get_IDcoords((int) constraints[i][2], (int) constraints[i][3], x1);
+        get_IDcoords((int) constraints[i][4], (int) constraints[i][5], x2);
+        get_IDcoords((int) constraints[i][6], (int) constraints[i][7], x3);
 
         // 1st bond
-        delx1 = x[atom1][0] - x[atom2][0];
-        dely1 = x[atom1][1] - x[atom2][1];
-        delz1 = x[atom1][2] - x[atom2][2];
+        delx1 = x1[0] - x2[0];
+        dely1 = x1[1] - x2[1];
+        delz1 = x1[2] - x2[2];
         domain->minimum_image(delx1,dely1,delz1); // ghost location fix
         rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
         r1 = sqrt(rsq1);
 
         // 2nd bond
-        delx2 = x[atom3][0] - x[atom2][0];
-        dely2 = x[atom3][1] - x[atom2][1];
-        delz2 = x[atom3][2] - x[atom2][2];
+        delx2 = x3[0] - x2[0];
+        dely2 = x3[1] - x2[1];
+        delz2 = x3[2] - x2[2];
         domain->minimum_image(delx2,dely2,delz2); // ghost location fix
         rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
         r2 = sqrt(rsq2);
@@ -1806,7 +1806,7 @@ int FixBondReact::check_constraints()
         c /= r1*r2;
         if (c > 1.0) c = 1.0;
         if (c < -1.0) c = -1.0;
-        if (acos(c) < constraints[i][5] || acos(c) > constraints[i][6]) return 0;
+        if (acos(c) < constraints[i][8] || acos(c) > constraints[i][9]) return 0;
       } else if (constraints[i][1] == DIHEDRAL) {
         // phi calculation from dihedral style harmonic
         atom1 = atom->map(glove[(int) constraints[i][2]-1][1]);
@@ -3304,14 +3304,12 @@ void FixBondReact::Constraints(char *line, int myrxn)
       constraints[nconstraints][7] = tmp[1]*tmp[1];
     } else if (strcmp(constraint_type,"angle") == 0) {
       constraints[nconstraints][1] = ANGLE;
-      sscanf(line,"%*s %lg %lg %lg %lg %lg",&tmp[0],&tmp[1],&tmp[2],&tmp[3],&tmp[4]);
-      if (tmp[0] > onemol->natoms || tmp[1] > onemol->natoms || tmp[2] > onemol->natoms)
-        error->one(FLERR,"Bond/react: Invalid template atom ID in map file");
-      constraints[nconstraints][2] = tmp[0];
-      constraints[nconstraints][3] = tmp[1];
-      constraints[nconstraints][4] = tmp[2];
-      constraints[nconstraints][5] = tmp[3]/180.0 * MY_PI;
-      constraints[nconstraints][6] = tmp[4]/180.0 * MY_PI;
+      sscanf(line,"%*s %s %s %s %lg %lg",strargs[0],strargs[1],strargs[2],&tmp[0],&tmp[1]);
+      readID(strargs[0], nconstraints, 2, 3);
+      readID(strargs[1], nconstraints, 4, 5);
+      readID(strargs[2], nconstraints, 6, 7);
+      constraints[nconstraints][8] = tmp[0]/180.0 * MY_PI;
+      constraints[nconstraints][9] = tmp[1]/180.0 * MY_PI;
     } else if (strcmp(constraint_type,"dihedral") == 0) {
       constraints[nconstraints][1] = DIHEDRAL;
       tmp[6] = 181.0; // impossible range

From 4250def29a0cacc9f0c236d8de0dd3710f36eea2 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 23 May 2020 22:58:37 -0600
Subject: [PATCH 03/11] dihedral constraint: fragment support

---
 src/USER-REACTION/fix_bond_react.cpp | 59 +++++++++++++---------------
 1 file changed, 28 insertions(+), 31 deletions(-)

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index 02b06313d4..df93ae982a 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -62,7 +62,7 @@ static const char cite_fix_bond_react[] =
 #define BIG 1.0e20
 #define DELTA 16
 #define MAXGUESS 20 // max # of guesses allowed by superimpose algorithm
-#define MAXCONARGS 10 // max # of arguments for any type of constraint + rxnID
+#define MAXCONARGS 14 // max # of arguments for any type of constraint + rxnID
 #define NUMVARVALS 4 // max # of keyword values that have variables as input
 
 // various statuses of superimpose algorithm:
@@ -1756,8 +1756,7 @@ evaluate constraints: return 0 if any aren't satisfied
 
 int FixBondReact::check_constraints()
 {
-  tagint atom1,atom2,atom3,atom4;
-  double x1[3],x2[3],x3[3];
+  double x1[3],x2[3],x3[3],x4[3];
   double delx,dely,delz,rsq;
   double delx1,dely1,delz1,delx2,dely2,delz2;
   double rsq1,rsq2,r1,r2,c,t,prrhob;
@@ -1767,6 +1766,7 @@ int FixBondReact::check_constraints()
   double s,phi;
   int ANDgate;
 
+  tagint atom1,atom2;
   double **x = atom->x;
 
   for (int i = 0; i < nconstraints; i++) {
@@ -1809,19 +1809,19 @@ int FixBondReact::check_constraints()
         if (acos(c) < constraints[i][8] || acos(c) > constraints[i][9]) return 0;
       } else if (constraints[i][1] == DIHEDRAL) {
         // phi calculation from dihedral style harmonic
-        atom1 = atom->map(glove[(int) constraints[i][2]-1][1]);
-        atom2 = atom->map(glove[(int) constraints[i][3]-1][1]);
-        atom3 = atom->map(glove[(int) constraints[i][4]-1][1]);
-        atom4 = atom->map(glove[(int) constraints[i][5]-1][1]);
+        get_IDcoords((int) constraints[i][2], (int) constraints[i][3], x1);
+        get_IDcoords((int) constraints[i][4], (int) constraints[i][5], x2);
+        get_IDcoords((int) constraints[i][6], (int) constraints[i][7], x3);
+        get_IDcoords((int) constraints[i][8], (int) constraints[i][9], x4);
 
-        vb1x = x[atom1][0] - x[atom2][0];
-        vb1y = x[atom1][1] - x[atom2][1];
-        vb1z = x[atom1][2] - x[atom2][2];
+        vb1x = x1[0] - x2[0];
+        vb1y = x1[1] - x2[1];
+        vb1z = x1[2] - x2[2];
         domain->minimum_image(vb1x,vb1y,vb1z);
 
-        vb2x = x[atom3][0] - x[atom2][0];
-        vb2y = x[atom3][1] - x[atom2][1];
-        vb2z = x[atom3][2] - x[atom2][2];
+        vb2x = x3[0] - x2[0];
+        vb2y = x3[1] - x2[1];
+        vb2z = x3[2] - x2[2];
         domain->minimum_image(vb2x,vb2y,vb2z);
 
         vb2xm = -vb2x;
@@ -1829,9 +1829,9 @@ int FixBondReact::check_constraints()
         vb2zm = -vb2z;
         domain->minimum_image(vb2xm,vb2ym,vb2zm);
 
-        vb3x = x[atom4][0] - x[atom3][0];
-        vb3y = x[atom4][1] - x[atom3][1];
-        vb3z = x[atom4][2] - x[atom3][2];
+        vb3x = x4[0] - x3[0];
+        vb3y = x4[1] - x3[1];
+        vb3z = x4[2] - x3[2];
         domain->minimum_image(vb3x,vb3y,vb3z);
 
         ax = vb1y*vb2zm - vb1z*vb2ym;
@@ -3312,21 +3312,18 @@ void FixBondReact::Constraints(char *line, int myrxn)
       constraints[nconstraints][9] = tmp[1]/180.0 * MY_PI;
     } else if (strcmp(constraint_type,"dihedral") == 0) {
       constraints[nconstraints][1] = DIHEDRAL;
-      tmp[6] = 181.0; // impossible range
-      tmp[7] = 182.0;
-      sscanf(line,"%*s %lg %lg %lg %lg %lg %lg %lg %lg",&tmp[0],&tmp[1],
-             &tmp[2],&tmp[3],&tmp[4],&tmp[5],&tmp[6],&tmp[7]);
-      if (tmp[0] > onemol->natoms || tmp[1] > onemol->natoms ||
-          tmp[2] > onemol->natoms || tmp[3] > onemol->natoms)
-        error->one(FLERR,"Bond/react: Invalid template atom ID in map file");
-      constraints[nconstraints][2] = tmp[0];
-      constraints[nconstraints][3] = tmp[1];
-      constraints[nconstraints][4] = tmp[2];
-      constraints[nconstraints][5] = tmp[3];
-      constraints[nconstraints][6] = tmp[4]/180.0 * MY_PI;
-      constraints[nconstraints][7] = tmp[5]/180.0 * MY_PI;
-      constraints[nconstraints][8] = tmp[6]/180.0 * MY_PI;
-      constraints[nconstraints][9] = tmp[7]/180.0 * MY_PI;
+      tmp[2] = 181.0; // impossible range
+      tmp[3] = 182.0;
+      sscanf(line,"%*s %s %s %s %s %lg %lg %lg %lg",strargs[0],strargs[1],
+             strargs[2],strargs[3],&tmp[0],&tmp[1],&tmp[2],&tmp[3]);
+      readID(strargs[0], nconstraints, 2, 3);
+      readID(strargs[1], nconstraints, 4, 5);
+      readID(strargs[2], nconstraints, 6, 7);
+      readID(strargs[3], nconstraints, 8, 9);
+      constraints[nconstraints][10] = tmp[0]/180.0 * MY_PI;
+      constraints[nconstraints][11] = tmp[1]/180.0 * MY_PI;
+      constraints[nconstraints][12] = tmp[2]/180.0 * MY_PI;
+      constraints[nconstraints][13] = tmp[3]/180.0 * MY_PI;
     } else if (strcmp(constraint_type,"arrhenius") == 0) {
       constraints[nconstraints][1] = ARRHENIUS;
       constraints[nconstraints][2] = narrhenius++;

From 60e0a8a6a80220a8f90af34a664e7dfb64eabb7e Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 23 May 2020 23:12:16 -0600
Subject: [PATCH 04/11] Update fix_bond_react.rst

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

diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst
index fc260de324..453b31e12f 100755
--- a/doc/src/fix_bond_react.rst
+++ b/doc/src/fix_bond_react.rst
@@ -300,7 +300,8 @@ either 'none' or 'charges.' Further details are provided in the
 discussion of the 'update_edges' keyword. The fifth 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 four types of constraints available, as discussed below.
+there are four types of constraints available, as discussed below:
+'distance', 'angle', 'dihedral', and 'arrhenius'.
 
 A sample map file is given below:
 
@@ -353,8 +354,9 @@ has syntax as follows:
    distance *ID1* *ID2* *rmin* *rmax*
 
 where 'distance' is the required keyword, *ID1* and *ID2* are
-pre-reaction atom IDs, and these two atoms must be separated by a
-distance between *rmin* and *rmax* for the reaction to occur.
+pre-reaction atom IDs (or molecule-fragment IDs, see below), and these
+two atoms must be separated by a distance between *rmin* and *rmax*
+for the reaction to occur.
 
 The constraint of type 'angle' has the following syntax:
 
@@ -363,11 +365,11 @@ The constraint of type 'angle' has the following syntax:
    angle *ID1* *ID2* *ID3* *amin* *amax*
 
 where 'angle' is the required keyword, *ID1*\ , *ID2* and *ID3* are
-pre-reaction atom IDs, and these three atoms must form an angle
-between *amin* and *amax* for the reaction to occur (where *ID2* is
-the central atom). Angles must be specified in degrees. This
-constraint can be used to enforce a certain orientation between
-reacting molecules.
+pre-reaction atom IDs (or molecule-fragment IDs, see below), and these
+three atoms must form an angle between *amin* and *amax* for the
+reaction to occur (where *ID2* is the central atom). Angles must be
+specified in degrees. This constraint can be used to enforce a certain
+orientation between reacting molecules.
 
 The constraint of type 'dihedral' has the following syntax:
 
@@ -376,15 +378,23 @@ The constraint of type 'dihedral' has the following syntax:
    dihedral *ID1* *ID2* *ID3* *ID4* *amin* *amax* *amin2* *amax2*
 
 where 'dihedral' is the required keyword, and *ID1*\ , *ID2*\ , *ID3*
-and *ID4* are pre-reaction atom IDs. Dihedral angles are calculated in
-the interval (-180,180]. Refer to the :doc:`dihedral style <dihedral_style>`
-documentation for further details on convention. If *amin* is less
-than *amax*, these four atoms must form a dihedral angle greater than
-*amin* **and** less than *amax* for the reaction to occur. If *amin*
-is greater than *amax*, these four atoms must form a dihedral angle
-greater than *amin* **or** less than *amax* for the reaction to occur.
-Angles must be specified in degrees. Optionally, a second range of
-permissible angles *amin2*-*amax2* can be specified.
+and *ID4* are pre-reaction atom IDs (or molecule-fragment IDs, see
+below). Dihedral angles are calculated in the interval (-180,180].
+Refer to the :doc:`dihedral style <dihedral_style>` documentation for
+further details on convention. If *amin* is less than *amax*, these
+four atoms must form a dihedral angle greater than *amin* **and** less
+than *amax* for the reaction to occur. If *amin* is greater than
+*amax*, these four atoms must form a dihedral angle greater than
+*amin* **or** less than *amax* for the reaction to occur. Angles must
+be specified in degrees. Optionally, a second range of permissible
+angles *amin2*-*amax2* can be specified.
+
+For the 'distance', 'angle', and 'dihedral' constraints (explained
+above), atom IDs can be replaced by pre-reaction molecule-fragment
+IDs. The molecule-fragment ID must begin with a letter. The location
+of the ID is the average of all atom positions in the fragment. The
+molecule fragment must have been defined in the :doc:`molecule <molecule>`
+command for the pre-reaction template.
 
 The constraint of type 'arrhenius' imposes an additional reaction
 probability according to the temperature-dependent Arrhenius equation:

From faec8ac2be1f5924aa9b655c875c63f779be7ace Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 23 May 2020 23:22:56 -0600
Subject: [PATCH 05/11] correctly update dihedral constraint

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

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index df93ae982a..c9ffd45774 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -1859,15 +1859,15 @@ int FixBondReact::check_constraints()
         phi = atan2(s,c);
 
         ANDgate = 0;
-        if (constraints[i][6] < constraints[i][7]) {
-          if (phi > constraints[i][6] && phi < constraints[i][7]) ANDgate = 1;
+        if (constraints[i][10] < constraints[i][11]) {
+          if (phi > constraints[i][10] && phi < constraints[i][11]) ANDgate = 1;
         } else {
-          if (phi > constraints[i][6] || phi < constraints[i][7]) ANDgate = 1;
+          if (phi > constraints[i][10] || phi < constraints[i][11]) ANDgate = 1;
         }
-        if (constraints[i][8] < constraints[i][9]) {
-          if (phi > constraints[i][8] && phi < constraints[i][9]) ANDgate = 1;
+        if (constraints[i][12] < constraints[i][13]) {
+          if (phi > constraints[i][12] && phi < constraints[i][13]) ANDgate = 1;
         } else {
-          if (phi > constraints[i][8] || phi < constraints[i][9]) ANDgate = 1;
+          if (phi > constraints[i][12] || phi < constraints[i][13]) ANDgate = 1;
         }
         if (ANDgate != 1) return 0;
       } else if (constraints[i][1] == ARRHENIUS) {

From 375fb4b3145717734b708bfae99843d04c722acb Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sat, 23 May 2020 23:32:49 -0600
Subject: [PATCH 06/11] bond/react: change 'general name' for method

and update contact info
---
 doc/src/Packages_details.rst                | 2 +-
 doc/utils/sphinx-config/false_positives.txt | 2 +-
 src/USER-REACTION/README                    | 8 ++++----
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/doc/src/Packages_details.rst b/doc/src/Packages_details.rst
index 8251c5301e..2353df777f 100644
--- a/doc/src/Packages_details.rst
+++ b/doc/src/Packages_details.rst
@@ -2064,7 +2064,7 @@ molecules, and chiral-sensitive reactions.
 * examples/USER/reaction
 * `2017 LAMMPS Workshop <https://lammps.sandia.gov/workshops/Aug17/pdf/gissinger.pdf>`_
 * `2019 LAMMPS Workshop <https://lammps.sandia.gov/workshops/Aug19/talk_gissinger.pdf>`_
-* disarmmd.org
+* reacter.org
 
 ----------
 
diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt
index 3d9f3ddcc9..25dc652321 100644
--- a/doc/utils/sphinx-config/false_positives.txt
+++ b/doc/utils/sphinx-config/false_positives.txt
@@ -629,7 +629,6 @@ dipolar
 dir
 Direc
 dirname
-disarmmd
 discoverable
 discretization
 discretized
@@ -2456,6 +2455,7 @@ rdc
 rdf
 RDideal
 rdx
+reacter
 README
 realtime
 reamin
diff --git a/src/USER-REACTION/README b/src/USER-REACTION/README
index f099ace103..ed2e760e64 100644
--- a/src/USER-REACTION/README
+++ b/src/USER-REACTION/README
@@ -1,4 +1,4 @@
-This package implements the DisARMMD protocol (disarmmd.org) as
+This package implements the REACTER protocol (reacter.org) as
 "fix bond/react." This fix allows for complex topology changes
 during a running MD simulation, when using classical force fields.
 Topology changes are defined in pre- and post-reaction molecule
@@ -6,7 +6,7 @@ templates and can include creation and deletion of bonds, angles,
 dihedrals, impropers, atom types, bond types, angle types,
 dihedral types, improper types, and/or atomic charges.
 
-The DisARMMD protocol is a method for modeling chemical reactions in
+The REACTER protocol is a method for modeling chemical reactions in
 classical molecular dynamics simulations. It was developed to build
 physically-realistic initial configurations for amorphous or
 crosslinked materials. Any number of competing or reversible reaction
@@ -15,9 +15,9 @@ advanced options currently available include reaction constraints
 (e.g. angle and Arrhenius constraints), deletion of reaction
 byproducts or other small molecules, and chiral-sensitive reactions.
 
-The DisARMMD methodology is detailed in:
+The REACTER methodology is detailed in:
     Gissinger et al., Polymer 128, 211-217 (2017)
     https://doi.org/10.1016/j.polymer.2017.09.038
 
-This package was created by Jacob Gissinger (info@disarmmd.org),
+This package was created by Jacob Gissinger (jrgiss05@gmail.com),
 while at the NASA Langley Research Center.

From 18320ded892b9c80aaad70559d75bacc2efb4190 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Sun, 24 May 2020 14:47:50 -0600
Subject: [PATCH 07/11] image correction for fragment location

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

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index c9ffd45774..d8bd61855d 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -1919,14 +1919,20 @@ void FixBondReact::get_IDcoords(int mode, int myID, double *center)
     for (int i = 0; i < 3; i++)
       center[i] = x[iatom][i];
   } else {
+    int iref = -1; // choose first atom as reference
+    int iatom;
     int nfragatoms = 0;
     for (int i = 0; i < 3; i++)
       center[i] = 0;
 
     for (int i = 0; i < onemol->natoms; i++) {
       if (onemol->fragmentmask[myID][i]) {
+        if (iref == -1)
+          iref = atom->map(glove[i][1]);
         for (int j = 0; j < 3; j++) {
-          center[j] += x[atom->map(glove[i][1])][j];
+          iatom = atom->map(glove[i][1]);
+          iatom = domain->closest_image(iref,iatom);
+          center[j] += x[iatom][j];
         }
         nfragatoms++;
       }

From db13dff4996342d23cffa6aa4e2611ac8f0a50a7 Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Wed, 27 May 2020 23:32:55 -0600
Subject: [PATCH 08/11] minor efficiency correction

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

diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index d8bd61855d..6ba5d1ce49 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -1929,11 +1929,10 @@ void FixBondReact::get_IDcoords(int mode, int myID, double *center)
       if (onemol->fragmentmask[myID][i]) {
         if (iref == -1)
           iref = atom->map(glove[i][1]);
-        for (int j = 0; j < 3; j++) {
-          iatom = atom->map(glove[i][1]);
-          iatom = domain->closest_image(iref,iatom);
+        iatom = atom->map(glove[i][1]);
+        iatom = domain->closest_image(iref,iatom);
+        for (int j = 0; j < 3; j++)
           center[j] += x[iatom][j];
-        }
         nfragatoms++;
       }
     }

From b0f6eafac8fcbc036df74051a877cd42f585e5eb Mon Sep 17 00:00:00 2001
From: Jacob Gissinger <jrgiss05@gmail.com>
Date: Wed, 27 May 2020 23:37:35 -0600
Subject: [PATCH 09/11] minor 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 453b31e12f..023d0cf4d5 100755
--- a/doc/src/fix_bond_react.rst
+++ b/doc/src/fix_bond_react.rst
@@ -392,9 +392,9 @@ angles *amin2*-*amax2* can be specified.
 For the 'distance', 'angle', and 'dihedral' constraints (explained
 above), atom IDs can be replaced by pre-reaction molecule-fragment
 IDs. The molecule-fragment ID must begin with a letter. The location
-of the ID is the average of all atom positions in the fragment. The
-molecule fragment must have been defined in the :doc:`molecule <molecule>`
-command for the pre-reaction template.
+of the ID is the unweighted average of all atom positions in the
+fragment. The molecule fragment must have been defined in the
+:doc:`molecule <molecule>` command for the pre-reaction template.
 
 The constraint of type 'arrhenius' imposes an additional reaction
 probability according to the temperature-dependent Arrhenius equation:

From acf280811a6e4d48395c8f189cfb2160a2927841 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Sun, 31 May 2020 11:14:50 -0400
Subject: [PATCH 10/11] step version strings for next patch release

---
 doc/lammps.1  | 2 +-
 src/version.h | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/doc/lammps.1 b/doc/lammps.1
index 80d07cbf0a..d9af339e0a 100644
--- a/doc/lammps.1
+++ b/doc/lammps.1
@@ -1,4 +1,4 @@
-.TH LAMMPS "5 May 2020" "2020-05-5"
+.TH LAMMPS "2 June 2020" "2020-06-02"
 .SH NAME
 .B LAMMPS
 \- Molecular Dynamics Simulator.
diff --git a/src/version.h b/src/version.h
index 6907e1c0a3..301c39539e 100644
--- a/src/version.h
+++ b/src/version.h
@@ -1 +1 @@
-#define LAMMPS_VERSION "5 May 2020"
+#define LAMMPS_VERSION "2 Jun 2020"

From 09ce0d1198a6c517568df0159998d4479f24ed42 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Sun, 31 May 2020 16:33:37 -0400
Subject: [PATCH 11/11] rephrase to use more common term

---
 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 023d0cf4d5..d5f917bfdb 100755
--- a/doc/src/fix_bond_react.rst
+++ b/doc/src/fix_bond_react.rst
@@ -392,7 +392,7 @@ angles *amin2*-*amax2* can be specified.
 For the 'distance', 'angle', and 'dihedral' constraints (explained
 above), atom IDs can be replaced by pre-reaction molecule-fragment
 IDs. The molecule-fragment ID must begin with a letter. The location
-of the ID is the unweighted average of all atom positions in the
+of the ID is the geometric center of all atom positions in the
 fragment. The molecule fragment must have been defined in the
 :doc:`molecule <molecule>` command for the pre-reaction template.