diff --git a/doc/src/fix_bond_react.txt b/doc/src/fix_bond_react.txt
index 8f71cd14ec7c4a405cd6f752db30ba7b37100e14..9c82e2f0ff2d4af88051c3a5483d2d9fd7d5bb14 100644
--- a/doc/src/fix_bond_react.txt
+++ b/doc/src/fix_bond_react.txt
@@ -36,10 +36,12 @@ react = mandatory argument indicating new reaction specification :l
template-ID(post-reacted) = ID of a molecule template containing post-reaction topology :l
map_file = name of file specifying corresponding atom-IDs in the pre- and post-reacted templates :l
zero or more individual keyword/value pairs may be appended to each react argument :l
- individual_keyword = {prob} or {stabilize_steps} or {update_edges} :l
+ individual_keyword = {prob} or {max_rxn} or {stabilize_steps} or {update_edges} :l
{prob} values = fraction seed
fraction = initiate reaction with this probability if otherwise eligible
seed = random number seed (positive integer)
+ {max_rxn} value = N
+ N = maximum number of reactions allowed to occur
{stabilize_steps} value = timesteps
timesteps = number of timesteps to apply the internally-created "nve/limit"_fix_nve_limit.html fix to reacting atoms
{update_edges} value = {none} or {charges} or {custom}
@@ -285,7 +287,8 @@ The {prob} keyword can affect whether 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
-fraction.
+fraction. Up to N reactions are permitted to occur, as optionally
+specified by the {max_rxn} keyword.
The {stabilize_steps} keyword allows for the specification of how many
timesteps a reaction site is stabilized before being returned to the
diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp
index e7dc816b8b73c3e4589dae93e5da04f67d32eecd..4d4642f102395a87598e315676f7043a3fbff86e 100644
--- a/src/USER-MISC/fix_bond_react.cpp
+++ b/src/USER-MISC/fix_bond_react.cpp
@@ -161,6 +161,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
memory->create(unreacted_mol,nreacts,"bond/react:unreacted_mol");
memory->create(reacted_mol,nreacts,"bond/react:reacted_mol");
memory->create(fraction,nreacts,"bond/react:fraction");
+ memory->create(max_rxn,nreacts,"bond/react:max_rxn");
memory->create(seed,nreacts,"bond/react:seed");
memory->create(limit_duration,nreacts,"bond/react:limit_duration");
memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag");
@@ -179,6 +180,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
for (int i = 0; i < nreacts; i++) {
fraction[i] = 1;
seed[i] = 12345;
+ max_rxn[i] = BIG;
stabilize_steps_flag[i] = 0;
update_edges_flag[i] = 0;
// set default limit duration to 60 timesteps
@@ -244,6 +246,13 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
if (seed[rxn] <= 0) error->all(FLERR,"Illegal fix bond/react command: "
"probability seed must be positive");
iarg += 3;
+ } 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[rxn] = force->inumeric(FLERR,arg[iarg+1]);
+ if (max_rxn[rxn] < 0) error->all(FLERR,"Illegal fix bond/react command: "
+ "'max_rxn' cannot be negative");
+ iarg += 2;
} else if (strcmp(arg[iarg],"stabilize_steps") == 0) {
if (stabilization_flag == 0) error->all(FLERR,"Stabilize_steps keyword "
"used without stabilization keyword");
@@ -379,9 +388,6 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
FixBondReact::~FixBondReact()
{
- // unregister callbacks to this fix from Atom class
- atom->delete_callback(id,0);
-
for (int i = 0; i < nreacts; i++) {
delete random[i];
}
@@ -396,6 +402,7 @@ FixBondReact::~FixBondReact()
memory->destroy(edge);
memory->destroy(equivalences);
memory->destroy(reverse_equiv);
+ memory->destroy(landlocked_atoms);
memory->destroy(custom_edges);
memory->destroy(delete_atoms);
@@ -405,6 +412,7 @@ FixBondReact::~FixBondReact()
memory->destroy(reacted_mol);
memory->destroy(fraction);
memory->destroy(seed);
+ memory->destroy(max_rxn);
memory->destroy(limit_duration);
memory->destroy(stabilize_steps_flag);
memory->destroy(update_edges_flag);
@@ -434,7 +442,6 @@ FixBondReact::~FixBondReact()
memory->destroy(restore);
memory->destroy(glove);
memory->destroy(pioneers);
- memory->destroy(landlocked_atoms);
memory->destroy(local_mega_glove);
memory->destroy(ghostly_mega_glove);
}
@@ -452,7 +459,7 @@ FixBondReact::~FixBondReact()
delete [] id_fix3;
}
- if (id_fix2 == NULL && modify->nfix) modify->delete_fix(id_fix2);
+ if (id_fix2 && modify->nfix) modify->delete_fix(id_fix2);
delete [] id_fix2;
delete [] statted_id;
@@ -748,6 +755,7 @@ void FixBondReact::post_integrate()
int j;
for (rxnID = 0; rxnID < nreacts; rxnID++) {
+ if (max_rxn[rxnID] <= reaction_count_total[rxnID]) continue;
for (int ii = 0; ii < nall; ii++) {
partner[ii] = 0;
finalpartner[ii] = 0;
@@ -1148,12 +1156,13 @@ void FixBondReact::superimpose_algorithm()
for (int i = 0; i < nreacts; i++)
reaction_count_total[i] += reaction_count[i];
- // this assumes compute_vector is called from process 0
- // ...so doesn't bother to bcast ghostly_rxn_count
if (me == 0)
for (int i = 0; i < nreacts; i++)
reaction_count_total[i] += ghostly_rxn_count[i];
+ // bcast ghostly_rxn_count
+ MPI_Bcast(&reaction_count_total[0], nreacts, MPI_INT, 0, world);
+
// this updates topology next step
next_reneighbor = update->ntimestep;
diff --git a/src/USER-MISC/fix_bond_react.h b/src/USER-MISC/fix_bond_react.h
index d6e7b785e7f8662d8d4eca52d5171d3a7a77be3c..2b85bbd28185e5f3dc2591b6b6744a1d5de15ab3 100644
--- a/src/USER-MISC/fix_bond_react.h
+++ b/src/USER-MISC/fix_bond_react.h
@@ -54,7 +54,7 @@ class FixBondReact : public Fix {
FILE *fp;
int *iatomtype,*jatomtype;
int *seed;
- double **cutsq,*fraction;
+ double **cutsq,*fraction,*max_rxn;
tagint lastcheck;
int stabilization_flag;
int custom_exclude_flag;