bond/react: Arrhenius constraint

doc/src/Eqs/fix_bond_react.tex

@@ -0,0 +1,9 @@
+\documentstyle[12pt]{article}
+\pagestyle{empty}
+\begin{document}
+
+\begin{eqnarray*}
+   k = AT^{n}e^{\frac{-E_{a}}{k_{B}T}}
+\end{eqnarray*}                           
+
+\end{document}

doc/src/fix_bond_react.txt

@@ -266,7 +266,7 @@ either 'none' or 'charges.' Further details are provided in the
 discussion of the 'update_edges' keyword. 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 two types of constraints available, as discussed below.
+there are three types of constraints available, as discussed below.
 
 A sample map file is given below:
 
@@ -320,6 +320,27 @@ 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 'arrhenius' imposes an additional reaction
+probability according to the temperature-dependent Arrhenius equation:
+
+:c,image(Eqs/fix_bond_react.jpg)
+
+The Arrhenius constraint has the following syntax:
+
+arrhenius {A} {n} {E_a} {seed} :pre
+
+where 'arrhenius' is the required keyword, {A} is the pre-exponential
+factor, {n} is the exponent of the temperature dependence, {E_a} is
+the activation energy ("units"_units.html of energy), and {seed} is a
+random number seed. The temperature is defined as the instantaneous
+temperature averaged over all atoms in the reaction site, and is
+calculated in the same manner as for example
+"compute_temp_chunk"_compute_temp_chunk.html. Currently, there are no
+options for additional temperature averaging or velocity-biased
+temperature calculations. A uniform random number between 0 and 1 is
+generated using {seed}; if this number is less than the result of the
+Arrhenius equation above, the reaction is permitted occur.
+
 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,

src/USER-MISC/fix_bond_react.cpp

@@ -71,7 +71,7 @@ static const char cite_fix_bond_react[] =
 enum{ACCEPT,REJECT,PROCEED,CONTINUE,GUESSFAIL,RESTORE};
 
 // types of available reaction constraints
-enum{DISTANCE,ANGLE};
+enum{DISTANCE,ANGLE,ARRHENIUS};
 
 /* ---------------------------------------------------------------------- */
 
@@ -100,6 +100,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
   extvector = 0;
   rxnID = 0;
   nconstraints = 0;
+  narrhenius = 0;
   status = PROCEED;
 
   nxspecial = NULL;
@@ -322,6 +323,16 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
     find_landlocked_atoms(i);
   }
 
+  // initialize Marsaglia RNG with processor-unique seed (Arrhenius prob)
+
+  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++) {
     delete [] files[i];
   }
@@ -348,7 +359,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
     }
   }
 
-  // initialize Marsaglia RNG with processor-unique seed
+  // initialize Marsaglia RNG with processor-unique seed ('prob' keyword)
 
   random = new class RanMars*[nreacts];
   for (int i = 0; i < nreacts; i++) {
@@ -1640,7 +1651,7 @@ int FixBondReact::check_constraints()
   tagint atom1,atom2,atom3;
   double delx,dely,delz,rsq;
   double delx1,dely1,delz1,delx2,dely2,delz2;
-  double rsq1,rsq2,r1,r2,c;
+  double rsq1,rsq2,r1,r2,c,t,prrhob;
 
   double **x = atom->x;
 
@@ -1680,12 +1691,54 @@ int FixBondReact::check_constraints()
         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;
+      } 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;
       }
     }
   }
   return 1;
 }
 
+/* ----------------------------------------------------------------------
+compute local temperature: average over all atoms in reaction template
+------------------------------------------------------------------------- */
+
+double FixBondReact::get_temperature()
+{
+  int i,ilocal;
+  double adof = domain->dimension;
+
+  double **v = atom->v;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int *type = atom->type;
+
+  double t = 0.0;
+
+  if (rmass) {
+    for (i = 0; i < onemol->natoms; i++) {
+      ilocal = atom->map(glove[i][1]);
+      t += (v[ilocal][0]*v[ilocal][0] + v[ilocal][1]*v[ilocal][1] +
+        v[ilocal][2]*v[ilocal][2]) * rmass[ilocal];
+    }
+  } else {
+    for (i = 0; i < onemol->natoms; i++) {
+      ilocal = atom->map(glove[i][1]);
+      t += (v[ilocal][0]*v[ilocal][0] + v[ilocal][1]*v[ilocal][1] +
+        v[ilocal][2]*v[ilocal][2]) * mass[type[ilocal]];
+    }
+  }
+
+  // final temperature
+  double dof = adof*onemol->natoms;
+  double tfactor = force->mvv2e / (dof * force->boltz);
+  t *= tfactor;
+  return t;
+}
+
 /* ----------------------------------------------------------------------
   Get xspecials for current molecule templates
 ------------------------------------------------------------------------- */
@@ -2947,6 +3000,14 @@ void FixBondReact::Constraints(char *line, int myrxn)
       constraints[nconstraints][4] = tmp[2];
       constraints[nconstraints][5] = tmp[3]/180.0 * MY_PI;
       constraints[nconstraints][6] = tmp[4]/180.0 * MY_PI;
+    } else if (strcmp(constraint_type,"arrhenius") == 0) {
+      constraints[nconstraints][1] = ARRHENIUS;
+      constraints[nconstraints][2] = 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];
     } else
       error->one(FLERR,"Bond/react: Illegal constraint type in 'Constraints' section of map file");
     nconstraints++;

src/USER-MISC/fix_bond_react.h

@@ -65,6 +65,7 @@ class FixBondReact : public Fix {
   int *stabilize_steps_flag;
   int *update_edges_flag;
   int nconstraints;
+  int narrhenius;
   double **constraints;
   int status;
   int *groupbits;
@@ -88,7 +89,8 @@ class FixBondReact : public Fix {
   Fix *fix2;              // properties/atom used to indicate 1) relaxing atoms
                           //                                  2) to which 'react' atom belongs
   Fix *fix3;              // property/atom used for system-wide thermostat
-  class RanMars **random;
+  class RanMars **random; // random number for 'prob' keyword
+  class RanMars **rrhandom; // random number for Arrhenius constraint
   class NeighList *list;
 
   int *reacted_mol,*unreacted_mol;
@@ -156,6 +158,7 @@ class FixBondReact : public Fix {
   void inner_crosscheck_loop();
   void ring_check();
   int check_constraints();
+  double get_temperature();
 
   void open(char *);
   void readline(char *);