small modification to fix bond/react to allow equal style variables as probability o

On branch master
Your branch is up to date with 'origin/master'.

Changes to be committed:
	modified:   doc/src/fix_bond_react.rst
	new file:   examples/USER/reaction/nylon,6-6_melt/in.large_nylon_melt_variable_probability
	modified:   src/USER-REACTION/fix_bond_react.cpp
	modified:   src/USER-REACTION/fix_bond_react.h

doc/src/fix_bond_react.rst

@@ -419,8 +419,10 @@ it occurs:
 
 The *prob* keyword can affect whether or not an eligible reaction
 actually occurs. The fraction setting must be a value between 0.0 and
-1.0. A uniform random number between 0.0 and 1.0 is generated and the
-eligible reaction only occurs if the random number is less than the
+1.0 or can be an equal style variable. In the later case the variable
+is evaluated during runtime and adjusted to be between 0.0 and 1.0 if 
+necessary. A uniform random number between 0.0 and 1.0 is generated and
+the eligible reaction only occurs if the random number is less than the
 fraction. Up to N reactions are permitted to occur, as optionally
 specified by the *max_rxn* keyword.
 

examples/USER/reaction/nylon,6-6_melt/in.large_nylon_melt_variable_probability

@@ -0,0 +1,55 @@
+# 35,000 atom nylon melt example
+
+units real
+
+boundary p p p
+
+atom_style full
+
+kspace_style pppm 1.0e-4
+
+pair_style lj/class2/coul/long 8.5
+
+angle_style class2
+
+bond_style class2
+
+dihedral_style class2
+
+improper_style class2
+
+read_data large_nylon_melt.data.gz
+
+variable runsteps equal 200
+varaible prob equal step/v_runsteps
+
+velocity all create 800.0 4928459 dist gaussian
+
+molecule mol1 rxn1_stp1_unreacted.data_template
+molecule mol2 rxn1_stp1_reacted.data_template
+molecule mol3 rxn1_stp2_unreacted.data_template
+molecule mol4 rxn1_stp2_reacted.data_template
+
+thermo 50
+
+# dump 1 all xyz 100 test_vis.xyz
+
+fix myrxns all bond/react stabilization yes statted_grp .03 &
+  react rxn1 all 1 0.0 2.9 mol1 mol2 rxn1_stp1_map &
+  react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map
+
+# stable at 800K
+fix 1 statted_grp_REACT nvt temp 800 800 100
+
+# in order to customize behavior of reacting atoms,
+# you can use the internally created 'bond_react_MASTER_group', like so:
+# fix 2 bond_react_MASTER_group temp/rescale 1 800 800 10 1
+
+thermo_style custom step temp press density f_myrxns[1] f_myrxns[2] # cumulative reaction counts
+
+# restart 100 restart1 restart2
+
+run 200
+
+# write_restart restart_longrun
+# write_data restart_longrun.data

src/USER-REACTION/fix_bond_react.cpp

@@ -39,6 +39,8 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu)
 #include "math_extra.h"
 #include "memory.h"
 #include "error.h"
+#include "input.h"
+#include "variable.h"
 
 #include <algorithm>
 
@@ -175,6 +177,8 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
   memory->create(nghostlyskips,nreacts,"bond/react:nghostlyskips");
   memory->create(seed,nreacts,"bond/react:seed");
   memory->create(limit_duration,nreacts,"bond/react:limit_duration");
+  memory->create(var_fraction_flag,nreacts,"bond/react:var_fraction_flag");
+  memory->create(var_fraction_id,nreacts,"bond/react:var_fraction_id");
   memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag");
   memory->create(update_edges_flag,nreacts,"bond/react:update_edges_flag");
   memory->create(constraints,1,MAXCONARGS,"bond/react:constraints");
@@ -193,6 +197,8 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
     fraction[i] = 1;
     seed[i] = 12345;
     max_rxn[i] = INT_MAX;
+    var_fraction_flag[i] = 0;
+    var_fraction_id[i] = 0;
     stabilize_steps_flag[i] = 0;
     update_edges_flag[i] = 0;
     // set default limit duration to 60 timesteps
@@ -251,7 +257,23 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
       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]);
+	// check if probability is a variable
+        if (strstr(arg[iarg+1],"v_") == arg[iarg+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]);
+	   var_fraction_flag[rxn] = 1.0;
+ 	   delete [] fracstr;
+ 	} else {
+	  // otherwise probability should be a number
+	  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: "
@@ -447,6 +469,8 @@ FixBondReact::~FixBondReact()
   memory->destroy(nlocalskips);
   memory->destroy(nghostlyskips);
   memory->destroy(limit_duration);
+  memory->destroy(var_fraction_flag);
+  memory->destroy(var_fraction_id);
   memory->destroy(stabilize_steps_flag);
   memory->destroy(update_edges_flag);
 
@@ -824,10 +848,17 @@ void FixBondReact::post_integrate()
       comm->reverse_comm_fix(this);
     }
 
+    // update reaction probability
+    if (var_fraction_flag[rxnID]) {
+      fraction[rxnID] = input->variable->compute_equal(var_fraction_id[rxnID]);
+      if (fraction[rxnID] < 0.0) fraction[rxnID] = 0.0;
+      if (fraction[rxnID] > 1.0) fraction[rxnID] = 1.0;
+    }
+
     // 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 (fraction[rxnID] < 1.0) {
       for (int i = 0; i < nlocal; i++)
       if (partner[i]) probability[i] = random[rxnID]->uniform();

src/USER-REACTION/fix_bond_react.h

@@ -64,6 +64,7 @@ class FixBondReact : public Fix {
   int custom_exclude_flag;
   int *stabilize_steps_flag;
   int *update_edges_flag;
+  int *var_fraction_flag, *var_fraction_id;
   int nconstraints;
   int narrhenius;
   double **constraints;