--- fix_bond_react.h.bak	2019-01-11 09:28:43.241582631 +0100
+++ fix_bond_react.h	2019-01-11 09:32:42.823094556 +0100
@@ -60,6 +60,8 @@
   int custom_exclude_flag;
   int *stabilize_steps_flag;
   int *update_edges_flag;
+  int *max_rxn_flag,*max_rxn;
+  int *var_fraction_flag, *var_fraction_id;
   int status;
   int *groupbits;
 
--- fix_bond_react.cpp.bak	2019-01-11 09:28:34.661528484 +0100
+++ fix_bond_react.cpp	2019-01-11 11:08:53.071508816 +0100
@@ -38,6 +38,8 @@
 #include "citeme.h"
 #include "memory.h"
 #include "error.h"
+#include "input.h"
+#include "variable.h"
 
 #include <algorithm>
 
@@ -165,6 +167,10 @@
   memory->create(limit_duration,nreacts,"bond/react:limit_duration");
   memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag");
   memory->create(update_edges_flag,nreacts,"bond/react:update_edges_flag");
+  memory->create(var_fraction_flag,nreacts,"bond/react:var_fraction_flag");
+  memory->create(var_fraction_id,nreacts,"bond/react:var_fraction_id");
+  memory->create(max_rxn_flag,nreacts,"bond/react:max_rxn_flag");
+  memory->create(max_rxn,nreacts,"bond/react:max_rxn");
   memory->create(iatomtype,nreacts,"bond/react:iatomtype");
   memory->create(jatomtype,nreacts,"bond/react:jatomtype");
   memory->create(ibonding,nreacts,"bond/react:ibonding");
@@ -178,9 +184,13 @@
 
   for (int i = 0; i < nreacts; i++) {
     fraction[i] = 1;
+    var_fraction_flag[i] = 0;
+    var_fraction_id[i] = 0;
     seed[i] = 12345;
     stabilize_steps_flag[i] = 0;
     update_edges_flag[i] = 0;
+    max_rxn_flag[i] = 0;
+    max_rxn[i] = 0;
     // set default limit duration to 60 timesteps
     limit_duration[i] = 60;
     reaction_count[i] = 0;
@@ -236,7 +246,19 @@
       if (strcmp(arg[iarg],"prob") == 0) {
         if (iarg+3 > narg) error->all(FLERR,"Illegal fix bond/react command: "
                                       "'prob' keyword has too few arguments");
-        fraction[rxn] = force->numeric(FLERR,arg[iarg+1]);
+        if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
+	  var_fraction_flag[rxn] = 1;
+	  int n = strlen(&arg[iarg+1][2]) + 1;
+	  char *fracstr = new char[n];
+	  strcpy(fracstr,&arg[iarg+1][2]);
+	  var_fraction_id[rxn] = input->variable->find(fracstr);
+	  if (var_fraction_id[rxn] < 0)
+	    error->all(FLERR,"variable name for fix bond/react does not exist");
+	  if (! input->variable->equalstyle(var_fraction_id[rxn]))
+	    error->all(FLERR,"variable in bond/react is not equal style");
+	  fraction[rxn] = input->variable->compute_equal(var_fraction_id[rxnID]);
+	  delete [] fracstr;
+	} else fraction[rxn] = force->numeric(FLERR,arg[iarg+1]);
         seed[rxn] = force->inumeric(FLERR,arg[iarg+2]);
         if (fraction[rxn] < 0.0 || fraction[rxn] > 1.0)
           error->all(FLERR,"Illegal fix bond/react command: "
@@ -260,6 +282,14 @@
         else if (strcmp(arg[iarg+1],"custom") == 0) update_edges_flag[rxn] = 2;
         else error->all(FLERR,"Illegal value for 'update_edges' keyword'");
         iarg += 2;
+      } else if (strcmp(arg[iarg],"max_rxn") == 0) {
+	if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: "
+				      "'max_rxn' has too few arguments");
+	max_rxn_flag[rxn] = 1;
+	max_rxn[rxn] = force->inumeric(FLERR,arg[iarg+1]);
+	if (max_rxn[rxn] <= 0) error->all(FLERR,"Illegal fix bond/react command "
+					  "'max_rxn' must be positive");
+	iarg += 2;
       } else error->all(FLERR,"Illegal fix bond/react command: unknown keyword");
     }
   }
@@ -404,10 +434,14 @@
   memory->destroy(unreacted_mol);
   memory->destroy(reacted_mol);
   memory->destroy(fraction);
+  memory->destroy(var_fraction_flag);
+  memory->destroy(var_fraction_id);
   memory->destroy(seed);
   memory->destroy(limit_duration);
   memory->destroy(stabilize_steps_flag);
   memory->destroy(update_edges_flag);
+  memory->destroy(max_rxn_flag);
+  memory->destroy(max_rxn);
 
   memory->destroy(iatomtype);
   memory->destroy(jatomtype);
@@ -748,6 +782,7 @@
 
   int j;
   for (rxnID = 0; rxnID < nreacts; rxnID++) {
+    if (max_rxn_flag[rxnID] && (max_rxn[rxnID] <= reaction_count_total[rxnID])) continue;
     for (int ii = 0; ii < nall; ii++) {
       partner[ii] = 0;
       finalpartner[ii] = 0;
@@ -771,7 +806,9 @@
     // each atom now knows its winning partner
     // for prob check, generate random value for each atom with a bond partner
     // forward comm of partner and random value, so ghosts have it
-
+    if (var_fraction_flag[rxnID] == 1) {
+      fraction[rxnID] = input->variable->compute_equal(var_fraction_id[rxnID]);
+    }
     if (fraction[rxnID] < 1.0) {
       for (int i = 0; i < nlocal; i++)
       if (partner[i]) probability[i] = random[rxnID]->uniform();
