add option to update all atoms' atomic charges

option to update all atomic charges, even when edge atoms are defined

doc/src/fix_bond_react.txt

@@ -41,7 +41,10 @@ react = mandatory argument indicating new reaction specification :l
       fraction = initiate reaction with this probability if otherwise eligible
       seed = random number seed (positive integer)
     {stabilize_steps} value = timesteps
-      timesteps = number of timesteps to apply internally created nve/limit.html :pre
+      timesteps = number of timesteps to apply internally created nve/limit.html
+    {update_edges} value = {none} or {charges} :l
+      none = do not update topology near the edges of reaction templates
+      charges = update atomic charges of all atoms in reaction templates :pre
 :ule
 
 [Examples:]
@@ -155,7 +158,17 @@ Some atoms in the pre-reacted template that are not reacting may have
 missing topology with respect to the simulation. For example, the
 pre-reacted template may contain an atom that would connect to the
 rest of a long polymer chain. These are referred to as edge atoms, and
-are also specified in the map file.
+are also specified in the map file. When the pre-reaction template
+contains edge atoms, not all atoms, bonds, charges, etc. specified in
+the reaction templates will be updated. Specifically, topology that
+involves only atoms that are 'too near' to template edges will not be
+updated. The definition of 'too near the edge' depends on which
+interactions are defined in the simulation. If the simulation has
+defined dihedrals, atoms within two bonds of edge atoms are considered
+'too near the edge.' If the simulation defines angles, but not
+dihedrals, atoms within one bond of edge atoms are considered 'too
+near the edge.' If just bonds are defined, only edge atoms are
+considered 'too near the edge.'
 
 Note that some care must be taken when a building a molecule template
 for a given simulation. All atom types in the pre-reacted template
@@ -255,6 +268,12 @@ The {stabilize_steps} keyword allows for the specification of how many
 timesteps a reaction site is stabilized before being returned to the
 overall system thermostat.
 
+The {update_edges} keyword can increase the number of atoms whose
+atomic charges are updated, when the pre-reaction template contains
+edge atoms. When the value is set to 'charges,' all atoms' atomic
+charges are updated to those specified by the post-reaction template,
+including atoms near the edge of reaction templates.
+
 In order to produce the most physical behavior, this 'reaction site
 equilibration time' should be tuned to be as small as possible while
 retaining stability for a given system or reaction step. After a
@@ -323,7 +342,7 @@ bond/break"_fix_bond_break.html, "fix bond/swap"_fix_bond_swap.html,
 
 [Default:]
 
-The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60
+The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60, update_edges = none
 
 :line
 

src/USER-MISC/fix_bond_react.cpp

@@ -163,6 +163,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
   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");
+  memory->create(update_edges_flag,nreacts,"bond/react:update_edges_flag");
   memory->create(iatomtype,nreacts,"bond/react:iatomtype");
   memory->create(jatomtype,nreacts,"bond/react:jatomtype");
   memory->create(ibonding,nreacts,"bond/react:ibonding");
@@ -178,6 +179,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
     fraction[i] = 1;
     seed[i] = 12345;
     stabilize_steps_flag[i] = 0;
+    update_edges_flag[i] = 0;
     // set default limit duration to 60 timesteps
     limit_duration[i] = 60;
     reaction_count[i] = 0;
@@ -249,6 +251,11 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
         limit_duration[rxn] = force->numeric(FLERR,arg[iarg+1]);
         stabilize_steps_flag[rxn] = 1;
         iarg += 2;
+      } else if (strcmp(arg[iarg],"update_edges") == 0) {
+        if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: "
+                                      "'update_edges' has too few arguments");
+        if (strcmp(arg[iarg+1],"charges") == 0) update_edges_flag[rxn] = 1;
+        iarg += 2;
       } else error->all(FLERR,"Illegal fix bond/react command: unknown keyword");
     }
   }
@@ -379,6 +386,7 @@ FixBondReact::~FixBondReact()
   memory->destroy(seed);
   memory->destroy(limit_duration);
   memory->destroy(stabilize_steps_flag);
+  memory->destroy(update_edges_flag);
 
   memory->destroy(iatomtype);
   memory->destroy(jatomtype);
@@ -2022,13 +2030,13 @@ void FixBondReact::update_everything()
     }
 
     // update charges and types of landlocked atoms
-    // here, add check for charge instead of requiring it
     for (int i = 0; i < update_num_mega; i++) {
       rxnID = update_mega_glove[0][i];
       twomol = atom->molecules[reacted_mol[rxnID]];
       for (int j = 0; j < twomol->natoms; j++) {
         int jj = equivalences[j][1][rxnID]-1;
-        if (landlocked_atoms[j][rxnID] == 1 && atom->map(update_mega_glove[jj+1][i]) >= 0 &&
+        if ((landlocked_atoms[j][rxnID] == 1 || update_edges_flag[rxnID] == 1) &&
+            atom->map(update_mega_glove[jj+1][i]) >= 0 &&
             atom->map(update_mega_glove[jj+1][i]) < nlocal) {
           type[atom->map(update_mega_glove[jj+1][i])] = twomol->type[j];
           if (twomol->qflag && atom->q_flag) {

src/USER-MISC/fix_bond_react.h

@@ -58,6 +58,7 @@ class FixBondReact : public Fix {
   tagint lastcheck;
   int stabilization_flag;
   int *stabilize_steps_flag;
+  int *update_edges_flag;
   int status;
   int *groupbits;