diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst
index 16bd4ecb71..38f061e05f 100644
--- a/doc/src/fix_bond_react.rst
+++ b/doc/src/fix_bond_react.rst
@@ -457,6 +457,23 @@ example, the molecule fragment could consist of only the backbone
 atoms of a polymer chain. This constraint can be used to enforce a
 specific relative position and orientation between reacting molecules.
 
+By default, all constraints must be satisfied for the reaction to
+occur. In other words, constraints are evaluated as a series of
+logical values using the logical AND operator "&&". More complex logic
+can be achieved by explicitly adding the logical AND operator "&&" or
+the logical OR operator "||" after a given constraint command. If a
+logical operator is specified after a constraint, it must be placed
+after all constraint parameters, on the same line as the constraint
+(one per line). Similarly, parentheses can be used to group
+constraints. The expression that results from concatenating all
+constraints should be a valid logical expression that can be read by
+the :doc:`variable <variable>` command after converting each
+constraint to a logical value. Because exactly one constraint is
+allowed per line, having a valid logical expression implies that left
+parentheses "(" should only appear before a constraint, and right
+parentheses ")" should only appear after a constraint and before any
+logical operator.
+
 Once a reaction site has been successfully identified, data structures
 within LAMMPS that store bond topology are updated to reflect the
 post-reacted molecule template. All force fields with fixed bonds,
@@ -599,8 +616,8 @@ reset_mol_ids = yes, custom_charges = no, molecule = off
 
 .. _Gissinger:
 
-**(Gissinger)** Gissinger, Jensen and Wise, Polymer, 128, 211 (2017).
+**(Gissinger)** Gissinger, Jensen and Wise, Polymer, 128, 211-217 (2017).
 
 .. _Gissinger2020:
 
-**(Gissinger)** Gissinger, Jensen and Wise, Macromolecules (2020, in press).
+**(Gissinger)** Gissinger, Jensen and Wise, Macromolecules, 53, 22, 9953–9961 (2020).
diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp
index efadec1b91..2bfc7714be 100644
--- a/src/USER-REACTION/fix_bond_react.cpp
+++ b/src/USER-REACTION/fix_bond_react.cpp
@@ -82,6 +82,9 @@ enum{ACCEPT,REJECT,PROCEED,CONTINUE,GUESSFAIL,RESTORE};
 // types of available reaction constraints
 enum{DISTANCE,ANGLE,DIHEDRAL,ARRHENIUS,RMSD};
 
+// ID type used by constraint
+enum{ATOM,FRAG};
+
 // keyword values that accept variables as input
 enum{NEVERY,RMIN,RMAX,PROB};
 
@@ -115,7 +118,6 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
   global_freq = 1;
   extvector = 0;
   rxnID = 0;
-  nconstraints = 0;
   narrhenius = 0;
   status = PROCEED;
 
@@ -206,7 +208,9 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
   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(nconstraints,nreacts,"bond/react:nconstraints");
+  memory->create(constraintstr,nreacts,MAXLINE,"bond/react:constraintstr");
+  memory->create(constraints,0,nreacts,"bond/react:constraints");
   memory->create(var_flag,NUMVARVALS,nreacts,"bond/react:var_flag");
   memory->create(var_id,NUMVARVALS,nreacts,"bond/react:var_id");
   memory->create(iatomtype,nreacts,"bond/react:iatomtype");
@@ -227,6 +231,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
     stabilize_steps_flag[i] = 0;
     custom_charges_fragid[i] = -1;
     molecule_keyword[i] = OFF;
+    nconstraints[i] = 0;
     // set default limit duration to 60 timesteps
     limit_duration[i] = 60;
     reaction_count[i] = 0;
@@ -440,9 +445,11 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
 
   rrhandom = new class RanMars*[narrhenius];
   int tmp = 0;
-  for (int i = 0; i < nconstraints; i++) {
-    if (constraints[i][1] == ARRHENIUS) {
-      rrhandom[tmp++] = new RanMars(lmp,(int) constraints[i][6] + me);
+  for (int i = 0; i < nreacts; i++) {
+    for (int j = 0; j < nconstraints[i]; j++) {
+      if (constraints[j][i].type == ARRHENIUS) {
+        rrhandom[tmp++] = new RanMars(lmp,(int) constraints[j][i].par[4] + me);
+      }
     }
   }
 
@@ -562,6 +569,9 @@ FixBondReact::~FixBondReact()
   memory->destroy(stabilize_steps_flag);
   memory->destroy(custom_charges_fragid);
   memory->destroy(molecule_keyword);
+  memory->destroy(constraints);
+  memory->destroy(nconstraints);
+  // need to delete rxn_name and constraintstr
 
   memory->destroy(iatomtype);
   memory->destroy(jatomtype);
@@ -1795,163 +1805,178 @@ int FixBondReact::check_constraints()
   tagint atom1,atom2;
   double **x = atom->x;
 
-  for (int i = 0; i < nconstraints; i++) {
-    if (constraints[i][0] == rxnID) {
-      if (constraints[i][1] == DISTANCE) {
-        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][6] || rsq > constraints[i][7]) return 0;
-      } else if (constraints[i][1] == ANGLE) {
-        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);
+  int *satisfied;
+  memory->create(satisfied,nconstraints[rxnID],"bond/react:satisfied");
+  for (int i = 0; i < nconstraints[rxnID]; i++)
+    satisfied[i] = 1;
 
-        // 1st bond
-        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);
+  for (int i = 0; i < nconstraints[rxnID]; i++) {
+    if (constraints[i][rxnID].type == DISTANCE) {
+      get_IDcoords(constraints[i][rxnID].idtype[0], constraints[i][rxnID].id[0], x1);
+      get_IDcoords(constraints[i][rxnID].idtype[1], constraints[i][rxnID].id[1], 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][rxnID].par[0] || rsq > constraints[i][rxnID].par[1]) satisfied[i] = 0;
+    } else if (constraints[i][rxnID].type == ANGLE) {
+      get_IDcoords(constraints[i][rxnID].idtype[0], constraints[i][rxnID].id[0], x1);
+      get_IDcoords(constraints[i][rxnID].idtype[1], constraints[i][rxnID].id[1], x2);
+      get_IDcoords(constraints[i][rxnID].idtype[2], constraints[i][rxnID].id[2], x3);
 
-        // 2nd bond
-        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);
+      // 1st bond
+      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);
 
-        // angle (cos and sin)
-        c = delx1*delx2 + dely1*dely2 + delz1*delz2;
-        c /= r1*r2;
-        if (c > 1.0) c = 1.0;
-        if (c < -1.0) c = -1.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
-        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);
+      // 2nd bond
+      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);
 
-        vb1x = x1[0] - x2[0];
-        vb1y = x1[1] - x2[1];
-        vb1z = x1[2] - x2[2];
-        domain->minimum_image(vb1x,vb1y,vb1z);
+      // angle (cos and sin)
+      c = delx1*delx2 + dely1*dely2 + delz1*delz2;
+      c /= r1*r2;
+      if (c > 1.0) c = 1.0;
+      if (c < -1.0) c = -1.0;
+      if (acos(c) < constraints[i][rxnID].par[0] || acos(c) > constraints[i][rxnID].par[1]) satisfied[i] = 0;
+    } else if (constraints[i][rxnID].type == DIHEDRAL) {
+      // phi calculation from dihedral style harmonic
+      get_IDcoords(constraints[i][rxnID].idtype[0], constraints[i][rxnID].id[0], x1);
+      get_IDcoords(constraints[i][rxnID].idtype[1], constraints[i][rxnID].id[1], x2);
+      get_IDcoords(constraints[i][rxnID].idtype[2], constraints[i][rxnID].id[2], x3);
+      get_IDcoords(constraints[i][rxnID].idtype[3], constraints[i][rxnID].id[3], x4);
 
-        vb2x = x3[0] - x2[0];
-        vb2y = x3[1] - x2[1];
-        vb2z = x3[2] - x2[2];
-        domain->minimum_image(vb2x,vb2y,vb2z);
+      vb1x = x1[0] - x2[0];
+      vb1y = x1[1] - x2[1];
+      vb1z = x1[2] - x2[2];
+      domain->minimum_image(vb1x,vb1y,vb1z);
 
-        vb2xm = -vb2x;
-        vb2ym = -vb2y;
-        vb2zm = -vb2z;
-        domain->minimum_image(vb2xm,vb2ym,vb2zm);
+      vb2x = x3[0] - x2[0];
+      vb2y = x3[1] - x2[1];
+      vb2z = x3[2] - x2[2];
+      domain->minimum_image(vb2x,vb2y,vb2z);
 
-        vb3x = x4[0] - x3[0];
-        vb3y = x4[1] - x3[1];
-        vb3z = x4[2] - x3[2];
-        domain->minimum_image(vb3x,vb3y,vb3z);
+      vb2xm = -vb2x;
+      vb2ym = -vb2y;
+      vb2zm = -vb2z;
+      domain->minimum_image(vb2xm,vb2ym,vb2zm);
 
-        ax = vb1y*vb2zm - vb1z*vb2ym;
-        ay = vb1z*vb2xm - vb1x*vb2zm;
-        az = vb1x*vb2ym - vb1y*vb2xm;
-        bx = vb3y*vb2zm - vb3z*vb2ym;
-        by = vb3z*vb2xm - vb3x*vb2zm;
-        bz = vb3x*vb2ym - vb3y*vb2xm;
+      vb3x = x4[0] - x3[0];
+      vb3y = x4[1] - x3[1];
+      vb3z = x4[2] - x3[2];
+      domain->minimum_image(vb3x,vb3y,vb3z);
 
-        rasq = ax*ax + ay*ay + az*az;
-        rbsq = bx*bx + by*by + bz*bz;
-        rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm;
-        rg = sqrt(rgsq);
+      ax = vb1y*vb2zm - vb1z*vb2ym;
+      ay = vb1z*vb2xm - vb1x*vb2zm;
+      az = vb1x*vb2ym - vb1y*vb2xm;
+      bx = vb3y*vb2zm - vb3z*vb2ym;
+      by = vb3z*vb2xm - vb3x*vb2zm;
+      bz = vb3x*vb2ym - vb3y*vb2xm;
 
-        ra2inv = rb2inv = 0.0;
-        if (rasq > 0) ra2inv = 1.0/rasq;
-        if (rbsq > 0) rb2inv = 1.0/rbsq;
-        rabinv = sqrt(ra2inv*rb2inv);
+      rasq = ax*ax + ay*ay + az*az;
+      rbsq = bx*bx + by*by + bz*bz;
+      rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm;
+      rg = sqrt(rgsq);
 
-        c = (ax*bx + ay*by + az*bz)*rabinv;
-        s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z);
+      ra2inv = rb2inv = 0.0;
+      if (rasq > 0) ra2inv = 1.0/rasq;
+      if (rbsq > 0) rb2inv = 1.0/rbsq;
+      rabinv = sqrt(ra2inv*rb2inv);
 
-        if (c > 1.0) c = 1.0;
-        if (c < -1.0) c = -1.0;
-        phi = atan2(s,c);
+      c = (ax*bx + ay*by + az*bz)*rabinv;
+      s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z);
 
-        ANDgate = 0;
-        if (constraints[i][10] < constraints[i][11]) {
-          if (phi > constraints[i][10] && phi < constraints[i][11]) ANDgate = 1;
-        } else {
-          if (phi > constraints[i][10] || phi < constraints[i][11]) 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][12] || phi < constraints[i][13]) ANDgate = 1;
-        }
-        if (ANDgate != 1) return 0;
-      } else if (constraints[i][1] == ARRHENIUS) {
-        t = get_temperature();
-        prrhob = constraints[i][3]*pow(t,constraints[i][4])*
-          exp(-constraints[i][5]/(force->boltz*t));
-        if (prrhob < rrhandom[(int) constraints[i][2]]->uniform()) return 0;
-      } else if (constraints[i][1] == RMSD) {
-        // call superpose
-        int iatom;
-        int iref = -1; // choose first atom as reference
-        int n2superpose = 0;
-        double **xfrozen; // coordinates for the "frozen" target molecule
-        double **xmobile; // coordinates for the "mobile" molecule
-        int ifragment = constraints[i][3];
-        if (ifragment >= 0) {
-          for (int j = 0; j < onemol->natoms; j++)
-            if (onemol->fragmentmask[ifragment][j]) n2superpose++;
-          memory->create(xfrozen,n2superpose,3,"bond/react:xfrozen");
-          memory->create(xmobile,n2superpose,3,"bond/react:xmobile");
-          int myincr = 0;
-          for (int j = 0; j < onemol->natoms; j++) {
-            if (onemol->fragmentmask[ifragment][j]) {
-              iatom = atom->map(glove[j][1]);
-              if (iref == -1) iref = iatom;
-              iatom = domain->closest_image(iref,iatom);
-              for (int k = 0; k < 3; k++) {
-                xfrozen[myincr][k] = x[iatom][k];
-                xmobile[myincr][k] = onemol->x[j][k];
-              }
-              myincr++;
-            }
-          }
-        } else {
-          int iatom;
-          int iref = -1; // choose first atom as reference
-          n2superpose = onemol->natoms;
-          memory->create(xfrozen,n2superpose,3,"bond/react:xfrozen");
-          memory->create(xmobile,n2superpose,3,"bond/react:xmobile");
-          for (int j = 0; j < n2superpose; j++) {
+      if (c > 1.0) c = 1.0;
+      if (c < -1.0) c = -1.0;
+      phi = atan2(s,c);
+
+      ANDgate = 0;
+      if (constraints[i][rxnID].par[0] < constraints[i][rxnID].par[1]) {
+        if (phi > constraints[i][rxnID].par[0] && phi < constraints[i][rxnID].par[1]) ANDgate = 1;
+      } else {
+        if (phi > constraints[i][rxnID].par[0] || phi < constraints[i][rxnID].par[1]) ANDgate = 1;
+      }
+      if (constraints[i][rxnID].par[2] < constraints[i][rxnID].par[3]) {
+        if (phi > constraints[i][rxnID].par[2] && phi < constraints[i][rxnID].par[3]) ANDgate = 1;
+      } else {
+        if (phi > constraints[i][rxnID].par[2] || phi < constraints[i][rxnID].par[3]) ANDgate = 1;
+      }
+      if (ANDgate != 1) satisfied[i] = 0;
+    } else if (constraints[i][rxnID].type == ARRHENIUS) {
+      t = get_temperature();
+      prrhob = constraints[i][rxnID].par[1]*pow(t,constraints[i][rxnID].par[2])*
+        exp(-constraints[i][rxnID].par[3]/(force->boltz*t));
+      if (prrhob < rrhandom[(int) constraints[i][rxnID].par[0]]->uniform()) satisfied[i] = 0;
+    } else if (constraints[i][rxnID].type == RMSD) {
+      // call superpose
+      int iatom;
+      int iref = -1; // choose first atom as reference
+      int n2superpose = 0;
+      double **xfrozen; // coordinates for the "frozen" target molecule
+      double **xmobile; // coordinates for the "mobile" molecule
+      int ifragment = constraints[i][rxnID].id[0];
+      if (ifragment >= 0) {
+        for (int j = 0; j < onemol->natoms; j++)
+          if (onemol->fragmentmask[ifragment][j]) n2superpose++;
+        memory->create(xfrozen,n2superpose,3,"bond/react:xfrozen");
+        memory->create(xmobile,n2superpose,3,"bond/react:xmobile");
+        int myincr = 0;
+        for (int j = 0; j < onemol->natoms; j++) {
+          if (onemol->fragmentmask[ifragment][j]) {
             iatom = atom->map(glove[j][1]);
             if (iref == -1) iref = iatom;
             iatom = domain->closest_image(iref,iatom);
             for (int k = 0; k < 3; k++) {
-              xfrozen[j][k] = x[iatom][k];
-              xmobile[j][k] = onemol->x[j][k];
+              xfrozen[myincr][k] = x[iatom][k];
+              xmobile[myincr][k] = onemol->x[j][k];
             }
+            myincr++;
+          }
+        }
+      } else {
+        int iatom;
+        int iref = -1; // choose first atom as reference
+        n2superpose = onemol->natoms;
+        memory->create(xfrozen,n2superpose,3,"bond/react:xfrozen");
+        memory->create(xmobile,n2superpose,3,"bond/react:xmobile");
+        for (int j = 0; j < n2superpose; j++) {
+          iatom = atom->map(glove[j][1]);
+          if (iref == -1) iref = iatom;
+          iatom = domain->closest_image(iref,iatom);
+          for (int k = 0; k < 3; k++) {
+            xfrozen[j][k] = x[iatom][k];
+            xmobile[j][k] = onemol->x[j][k];
           }
         }
-        Superpose3D<double, double **> superposer(n2superpose);
-        double rmsd = superposer.Superpose(xfrozen, xmobile);
-        if (rmsd > constraints[i][2]) return 0;
-        memory->destroy(xfrozen);
-        memory->destroy(xmobile);
       }
+      Superpose3D<double, double **> superposer(n2superpose);
+      double rmsd = superposer.Superpose(xfrozen, xmobile);
+      memory->destroy(xfrozen);
+      memory->destroy(xmobile);
+      if (rmsd > constraints[i][rxnID].par[0]) satisfied[i] = 0;
     }
   }
 
+  if (nconstraints[rxnID] > 0) {
+    char evalstr[MAXLINE],*ptr,valstr;
+    strcpy(evalstr,constraintstr[rxnID]);
+    for (int i = 0; i < nconstraints[rxnID]; i++) {
+      sprintf(&valstr,"%d", satisfied[i]);
+      ptr = strchr(evalstr,'C');
+      *ptr = valstr;
+    }
+    double verdict = input->variable->evaluate_boolean(evalstr);
+    if (verdict == 0.0) return 0;
+  }
+
   // let's also check chirality within 'check_constraint'
   for (int i = 0; i < onemol->natoms; i++) {
     if (chiral_atoms[i][0][rxnID] == 1) {
@@ -1975,6 +2000,7 @@ int FixBondReact::check_constraints()
     }
   }
 
+  memory->destroy(satisfied);
   return 1;
 }
 
@@ -1987,7 +2013,7 @@ fragment: given pre-reacted molID (onemol) and fragID,
 void FixBondReact::get_IDcoords(int mode, int myID, double *center)
 {
   double **x = atom->x;
-  if (mode == 1) {
+  if (mode == ATOM) {
     int iatom = atom->map(glove[myID-1][1]);
     for (int i = 0; i < 3; i++)
       center[i] = x[iatom][i];
@@ -3216,8 +3242,8 @@ void FixBondReact::read(int myrxn)
     else if (strstr(line,"deleteIDs")) sscanf(line,"%d",&ndelete);
     else if (strstr(line,"chiralIDs")) sscanf(line,"%d",&nchiral);
     else if (strstr(line,"constraints")) {
-      sscanf(line,"%d",&nconstr);
-      memory->grow(constraints,nconstraints+nconstr,MAXCONARGS,"bond/react:constraints");
+      sscanf(line,"%d",&nconstraints[myrxn]);
+      memory->grow(constraints,nconstraints[myrxn],nreacts,"bond/react:constraints");
     } else break;
   }
 
@@ -3252,7 +3278,7 @@ void FixBondReact::read(int myrxn)
     } else if (strcmp(keyword,"ChiralIDs") == 0) {
       ChiralCenters(line, myrxn);
     } else if (strcmp(keyword,"Constraints") == 0) {
-      Constraints(line, myrxn);
+      ReadConstraints(line, myrxn);
     } else error->one(FLERR,"Bond/react: Unknown section in map file");
 
     parse_keyword(1,line,keyword);
@@ -3350,69 +3376,89 @@ void FixBondReact::ChiralCenters(char *line, int myrxn)
   }
 }
 
-void FixBondReact::Constraints(char *line, int myrxn)
+void FixBondReact::ReadConstraints(char *line, int myrxn)
 {
   double tmp[MAXCONARGS];
-  char **strargs;
+  char **strargs,*ptr;
   memory->create(strargs,MAXCONARGS,MAXLINE,"bond/react:strargs");
   char *constraint_type = new char[MAXLINE];
-  for (int i = 0; i < nconstr; i++) {
+  strcpy(constraintstr[myrxn],"("); // string for boolean constraint logic
+  for (int i = 0; i < nconstraints[myrxn]; i++) {
     readline(line);
+    if (ptr = strrchr(line,'(')) { // reverse char search
+      strncat(constraintstr[myrxn],line,ptr-line+1);
+      line = ptr + 1;
+    }
+    // 'C' indicates where to sub in next constraint
+    strcat(constraintstr[myrxn],"C");
+    if (ptr = strchr(line,')')) {
+      strncat(constraintstr[myrxn],ptr,strrchr(line,')')-ptr+1);
+    }
+    if (ptr = strstr(line,"&&")) {
+      strcat(constraintstr[myrxn],"&&");
+      *ptr = '\0';
+    } else if (ptr = strstr(line,"||")) {
+      strcat(constraintstr[myrxn],"||");
+      *ptr = '\0';
+    } else if (i+1 < nconstraints[myrxn]){
+      strcat(constraintstr[myrxn],"&&");
+    }
+    if (ptr = strchr(line,')'))
+      *ptr = '\0';
     sscanf(line,"%s",constraint_type);
-    constraints[nconstraints][0] = myrxn;
     if (strcmp(constraint_type,"distance") == 0) {
-      constraints[nconstraints][1] = DISTANCE;
+      constraints[i][myrxn].type = DISTANCE;
       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);
+      readID(strargs[0], i, myrxn, 0);
+      readID(strargs[1], i, myrxn, 1);
       // cutoffs
-      constraints[nconstraints][6] = tmp[0]*tmp[0]; // using square of distance
-      constraints[nconstraints][7] = tmp[1]*tmp[1];
+      constraints[i][myrxn].par[0] = tmp[0]*tmp[0]; // using square of distance
+      constraints[i][myrxn].par[1] = tmp[1]*tmp[1];
     } else if (strcmp(constraint_type,"angle") == 0) {
-      constraints[nconstraints][1] = ANGLE;
+      constraints[i][myrxn].type = ANGLE;
       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;
+      readID(strargs[0], i, myrxn, 0);
+      readID(strargs[1], i, myrxn, 1);
+      readID(strargs[2], i, myrxn, 2);
+      constraints[i][myrxn].par[0] = tmp[0]/180.0 * MY_PI;
+      constraints[i][myrxn].par[1] = tmp[1]/180.0 * MY_PI;
     } else if (strcmp(constraint_type,"dihedral") == 0) {
-      constraints[nconstraints][1] = DIHEDRAL;
+      constraints[i][myrxn].type = DIHEDRAL;
       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;
+      readID(strargs[0], i, myrxn, 0);
+      readID(strargs[1], i, myrxn, 1);
+      readID(strargs[2], i, myrxn, 2);
+      readID(strargs[3], i, myrxn, 3);
+      constraints[i][myrxn].par[0] = tmp[0]/180.0 * MY_PI;
+      constraints[i][myrxn].par[1] = tmp[1]/180.0 * MY_PI;
+      constraints[i][myrxn].par[2] = tmp[2]/180.0 * MY_PI;
+      constraints[i][myrxn].par[3] = tmp[3]/180.0 * MY_PI;
     } else if (strcmp(constraint_type,"arrhenius") == 0) {
-      constraints[nconstraints][1] = ARRHENIUS;
-      constraints[nconstraints][2] = narrhenius++;
+      constraints[i][myrxn].type = ARRHENIUS;
+      constraints[i][myrxn].par[0] = narrhenius++;
       sscanf(line,"%*s %lg %lg %lg %lg",&tmp[0],&tmp[1],&tmp[2],&tmp[3]);
-      constraints[nconstraints][3] = tmp[0];
-      constraints[nconstraints][4] = tmp[1];
-      constraints[nconstraints][5] = tmp[2];
-      constraints[nconstraints][6] = tmp[3];
+      constraints[i][myrxn].par[1] = tmp[0];
+      constraints[i][myrxn].par[2] = tmp[1];
+      constraints[i][myrxn].par[3] = tmp[2];
+      constraints[i][myrxn].par[4] = tmp[3];
     } else if (strcmp(constraint_type,"rmsd") == 0) {
-      constraints[nconstraints][1] = RMSD;
+      constraints[i][myrxn].type = RMSD;
       strcpy(strargs[0],"0");
       sscanf(line,"%*s %lg %s",&tmp[0],strargs[0]);
-      constraints[nconstraints][2] = tmp[0]; // RMSDmax
-      constraints[nconstraints][3] = -1; // optional molecule fragment
+      constraints[i][myrxn].par[0] = tmp[0]; // RMSDmax
+      constraints[i][myrxn].id[0] = -1; // optional molecule fragment
       if (isalpha(strargs[0][0])) {
         int ifragment = onemol->findfragment(strargs[0]);
         if (ifragment < 0) error->one(FLERR,"Bond/react: Molecule fragment does not exist");
-        else constraints[nconstraints][3] = ifragment;
+        else constraints[i][myrxn].id[0] = ifragment;
       }
     } else
       error->one(FLERR,"Bond/react: Illegal constraint type in 'Constraints' section of map file");
-    nconstraints++;
   }
+  strcat(constraintstr[myrxn],")"); // close boolean constraint logic string
   delete [] constraint_type;
   memory->destroy(strargs);
 }
@@ -3422,18 +3468,18 @@ 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)
+void FixBondReact::readID(char *strarg, int iconstr, int myrxn, int i)
 {
   if (isalpha(strarg[0])) {
-    constraints[iconstr][mode] = 0; // fragment vs. atom ID flag
+    constraints[iconstr][myrxn].idtype[i] = FRAG; // 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;
+    constraints[iconstr][myrxn].id[i] = ifragment;
   } else {
-    constraints[iconstr][mode] = 1; // fragment vs. atom ID flag
+    constraints[iconstr][myrxn].idtype[i] = ATOM; // 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;
+    constraints[iconstr][myrxn].id[i] = iatom;
   }
 }
 
diff --git a/src/USER-REACTION/fix_bond_react.h b/src/USER-REACTION/fix_bond_react.h
index 6105859136..a66f57c2bc 100644
--- a/src/USER-REACTION/fix_bond_react.h
+++ b/src/USER-REACTION/fix_bond_react.h
@@ -31,7 +31,9 @@ namespace LAMMPS_NS {
 class FixBondReact : public Fix {
  public:
 
-  enum {MAXLINE=256};
+  enum {MAXLINE=256}; // max length of line read from files
+  enum {MAXCONIDS=4}; // max # of IDs used by any constraint
+  enum {MAXCONPAR=4}; // max # of constraint parameters
 
   FixBondReact(class LAMMPS *, int, char **);
   ~FixBondReact();
@@ -66,9 +68,9 @@ class FixBondReact : public Fix {
   int *stabilize_steps_flag;
   int *custom_charges_fragid;
   int *molecule_keyword;
-  int nconstraints;
+  int *nconstraints;
+  char **constraintstr;
   int narrhenius;
-  double **constraints;
   int **var_flag,**var_id; // for keyword values with variable inputs
   int status;
   int *groupbits;
@@ -114,7 +116,7 @@ class FixBondReact : public Fix {
 
   int *ibonding,*jbonding;
   int *closeneigh; // indicates if bonding atoms of a rxn are 1-2, 1-3, or 1-4 neighbors
-  int nedge,nequivalent,ndelete,nchiral,nconstr; // # edge, equivalent atoms in mapping file
+  int nedge,nequivalent,ndelete,nchiral; // # edge, equivalent atoms in mapping file
   int attempted_rxn; // there was an attempt!
   int *local_rxn_count;
   int *ghostly_rxn_count;
@@ -155,7 +157,7 @@ class FixBondReact : public Fix {
   void DeleteAtoms(char *, int);
   void CustomCharges(int, int);
   void ChiralCenters(char *, int);
-  void Constraints(char *, int);
+  void ReadConstraints(char *, int);
   void readID(char *, int, int, int);
 
   void make_a_guess ();
@@ -195,6 +197,14 @@ class FixBondReact : public Fix {
   };
   Set *set;
 
+  struct Constraint {
+    int type;
+    int id[MAXCONIDS];
+    int idtype[MAXCONIDS];
+    double par[MAXCONPAR];
+  };
+  Constraint **constraints;
+
   // DEBUG
 
   void print_bb();