From b0af54ac374f90adb28c0b5a22f2f92495c2a21e Mon Sep 17 00:00:00 2001 From: jrgissing Date: Sun, 13 Jan 2019 00:07:03 -0700 Subject: [PATCH] bond/react:limit reaction occurrences --- doc/src/fix_bond_react.txt | 7 +++++-- src/USER-MISC/fix_bond_react.cpp | 23 ++++++++++++++++------- src/USER-MISC/fix_bond_react.h | 2 +- 3 files changed, 22 insertions(+), 10 deletions(-) diff --git a/doc/src/fix_bond_react.txt b/doc/src/fix_bond_react.txt index 8f71cd14ec..9c82e2f0ff 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 e7dc816b8b..4d4642f102 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 d6e7b785e7..2b85bbd281 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;