Merge pull request #4293 from jrgissing/bond_react_fixes_summ24

Bond/react: fixes and error checking

doc/src/fix_bond_react.rst

@@ -785,3 +785,7 @@ reset_mol_ids = yes, custom_charges = no, molecule = off, modify_create = *fit a
 .. _Gissinger2020:
 
 **(Gissinger2020)** Gissinger, Jensen and Wise, Macromolecules, 53, 22, 9953-9961 (2020).
+
+.. _Gissinger2024:
+
+**(Gissinger2024)** Gissinger, Jensen and Wise, Computer Physics Communications, 304, 109287 (2024).

examples/PACKAGES/reaction/tiny_nylon/in.tiny_nylon.stabilized

@@ -44,6 +44,7 @@ thermo 50
 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
+  react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map rescale_charges yes
 
 fix 1 statted_grp_REACT nvt temp 300 300 100
 

examples/PACKAGES/reaction/tiny_nylon/in.tiny_nylon.stabilized_variable_probability

@@ -47,7 +47,7 @@ thermo 50
 
 fix myrxns all bond/react stabilization yes statted_grp .03 &
   react rxn1 all 1 0.0 5.0 mol1 mol2 rxn1_stp1_map prob v_prob1 1234 &
-  react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map prob v_prob2 1234
+  react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map prob v_prob2 1234 rescale_charges yes
 
 fix 1 statted_grp_REACT nvt temp 300 300 100
 

examples/PACKAGES/reaction/tiny_nylon/in.tiny_nylon.unstabilized

@@ -44,7 +44,7 @@ thermo 50
 
 fix myrxns all bond/react stabilization no &
   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
+  react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map rescale_charges yes
 
 fix 1 all nve/limit .03
 

examples/PACKAGES/reaction/tiny_nylon/rxn1_stp1_reacted.molecule_template

@@ -48,27 +48,6 @@ Types
      17 hc
      18 hc
 
-Charges
-
-      1 -0.300000
-      2  0.000000
-      3  0.000000
-      4  0.000000
-      5  0.000000
-      6  0.000000
-      7  0.000000
-      8  0.000000
-      9  0.000000
-     10  0.300000
-     11  0.000000
-     12  0.000000
-     13  0.000000
-     14  0.000000
-     15  0.000000
-     16  0.000000
-     17  0.000000
-     18  0.000000
-
 Molecules
 
       1      1

examples/PACKAGES/reaction/tiny_nylon/rxn1_stp2_reacted.molecule_template

@@ -44,21 +44,21 @@ Types
 
 Charges
 
-      1 -0.300000
+      1 -0.60533
-      2  0.000000
+      2 -0.01149
-      3  0.000000
+      3 -0.76306
-      4  0.410000
+      4 0.38
-      5  0.000000
+      5 0.29346
-      6  0.000000
+      6 0.18360
-      7  0.000000
+      7 0.15396
-      8  0.000000
+      8 -0.72636
-      9  0.000000
+      9 -0.27437
-     10  0.300000
+     10 0.40603
-     11  0.000000
+     11 -0.65530
-     12 -0.820000
+     12 -0.76
-     13  0.000000
+     13 0.21423
-     14  0.000000
+     14 0.18949
-     15  0.410000
+     15 0.38
 
 Molecules
 

src/REACTION/fix_bond_react.cpp

@@ -58,7 +58,7 @@ using namespace MathConst;
 
 static const char cite_fix_bond_react[] =
     "fix bond/react: reacter.org doi:10.1016/j.polymer.2017.09.038, "
-    "doi:10.1021/acs.macromol.0c02012\n\n"
+    "doi:10.1021/acs.macromol.0c02012, doi:10.1016/j.cpc.2024.109287\n\n"
     "@Article{Gissinger17,\n"
     " author = {J. R. Gissinger and B. D. Jensen and K. E. Wise},\n"
     " title = {Modeling Chemical Reactions in Classical Molecular Dynamics Simulations},\n"
@@ -75,6 +75,14 @@ static const char cite_fix_bond_react[] =
     " volume =  53,\n"
     " number =  22,\n"
     " pages =   {9953--9961}\n"
+    "}\n\n"
+    "@Article{Gissinger24,\n"
+    " author = {J. R. Gissinger, B. D. Jensen, K. E. Wise},\n"
+    " title = {Molecular Modeling of Reactive Systems with REACTER},\n"
+    " journal = {Computer Physics Communications},\n"
+    " year =    2024,\n"
+    " volume =  304,\n"
+    " number =  109287\n"
     "}\n\n";
 
 static constexpr double BIG = 1.0e20;
@@ -225,8 +233,8 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
   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(nlocalskips,nreacts,"bond/react:nlocalskips");
+  memory->create(nlocalkeep,nreacts,"bond/react:nlocalkeep");
-  memory->create(nghostlyskips,nreacts,"bond/react:nghostlyskips");
+  memory->create(nghostlykeep,nreacts,"bond/react:nghostlykeep");
   memory->create(seed,nreacts,"bond/react:seed");
   memory->create(limit_duration,nreacts,"bond/react:limit_duration");
   memory->create(rate_limit,3,nreacts,"bond/react:rate_limit");
@@ -486,10 +494,6 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
     get_molxspecials();
     read_map_file(i);
     fclose(fp);
-    if (ncreate == 0 && onemol->natoms != twomol->natoms)
-      error->all(FLERR,"Fix bond/react: Reaction templates must contain the same number of atoms");
-    else if (ncreate > 0 && onemol->natoms + ncreate != twomol->natoms)
-      error->all(FLERR,"Fix bond/react: Incorrect number of created atoms");
     iatomtype[i] = onemol->type[ibonding[i]-1];
     jatomtype[i] = onemol->type[jbonding[i]-1];
     find_landlocked_atoms(i);
@@ -644,8 +648,8 @@ FixBondReact::~FixBondReact()
   memory->destroy(fraction);
   memory->destroy(seed);
   memory->destroy(max_rxn);
-  memory->destroy(nlocalskips);
+  memory->destroy(nlocalkeep);
-  memory->destroy(nghostlyskips);
+  memory->destroy(nghostlykeep);
   memory->destroy(limit_duration);
   memory->destroy(var_flag);
   memory->destroy(var_id);
@@ -716,6 +720,7 @@ int FixBondReact::setmask()
   int mask = 0;
   mask |= POST_INTEGRATE;
   mask |= POST_INTEGRATE_RESPA;
+  mask |= POST_FORCE;
   return mask;
 }
 
@@ -872,8 +877,8 @@ void FixBondReact::post_integrate()
     reaction_count[i] = 0;
     local_rxn_count[i] = 0;
     ghostly_rxn_count[i] = 0;
-    nlocalskips[i] = 0;
+    nlocalkeep[i] = INT_MAX;
-    nghostlyskips[i] = 0;
+    nghostlykeep[i] = INT_MAX;
     // update reaction probability
     if (var_flag[PROB][i])
       fraction[i] = input->variable->compute_equal(var_id[PROB][i]);
@@ -1424,10 +1429,13 @@ void FixBondReact::superimpose_algorithm()
   MPI_Allreduce(&local_rxn_count[0],&reaction_count[0],nreacts,MPI_INT,MPI_SUM,world);
 
   int rxnflag = 0;
+  int *delta_rxn;
+  memory->create(delta_rxn,nreacts,"bond/react:delta_rxn");
   if (comm->me == 0)
     for (int i = 0; i < nreacts; i++) {
-      reaction_count_total[i] += reaction_count[i] + ghostly_rxn_count[i];
+      delta_rxn[i] = reaction_count[i] + ghostly_rxn_count[i];
-      rxnflag += reaction_count[i] + ghostly_rxn_count[i];
+      reaction_count_total[i] += delta_rxn[i];
+      rxnflag += delta_rxn[i];
     }
 
   MPI_Bcast(&reaction_count_total[0], nreacts, MPI_INT, 0, world);
@@ -1460,42 +1468,43 @@ void FixBondReact::superimpose_algorithm()
     if (overstep > 0) {
       // let's randomly choose rxns to skip, unbiasedly from local and ghostly
       int *local_rxncounts;
-      int *all_localskips;
+      int *all_localkeep;
       memory->create(local_rxncounts,nprocs,"bond/react:local_rxncounts");
-      memory->create(all_localskips,nprocs,"bond/react:all_localskips");
+      memory->create(all_localkeep,nprocs,"bond/react:all_localkeep");
       MPI_Gather(&local_rxn_count[i],1,MPI_INT,local_rxncounts,1,MPI_INT,0,world);
       if (comm->me == 0) {
-        int delta_rxn = reaction_count[i] + ghostly_rxn_count[i];
         // when using variable input for rate_limit, rate_limit_overstep could be > delta_rxn (below)
         // we need to limit overstep to the number of reactions on this timestep
         // essentially skipping all reactions, would be more efficient to use a skip_all flag
-        if (overstep > delta_rxn) overstep = delta_rxn;
+        if (overstep > delta_rxn[i]) overstep = delta_rxn[i];
+        int nkeep = delta_rxn[i] - overstep;
         int *rxn_by_proc;
-        memory->create(rxn_by_proc,delta_rxn,"bond/react:rxn_by_proc");
+        memory->create(rxn_by_proc,delta_rxn[i],"bond/react:rxn_by_proc");
-        for (int j = 0; j < delta_rxn; j++)
+        for (int j = 0; j < delta_rxn[i]; j++)
           rxn_by_proc[j] = -1; // corresponds to ghostly
         int itemp = 0;
         for (int j = 0; j < nprocs; j++)
           for (int k = 0; k < local_rxncounts[j]; k++)
             rxn_by_proc[itemp++] = j;
-        std::shuffle(&rxn_by_proc[0],&rxn_by_proc[delta_rxn], park_rng);
+        std::shuffle(&rxn_by_proc[0],&rxn_by_proc[delta_rxn[i]], park_rng);
         for (int j = 0; j < nprocs; j++)
-          all_localskips[j] = 0;
+          all_localkeep[j] = 0;
-        nghostlyskips[i] = 0;
+        nghostlykeep[i] = 0;
-        for (int j = 0; j < overstep; j++) {
+        for (int j = 0; j < nkeep; j++) {
-          if (rxn_by_proc[j] == -1) nghostlyskips[i]++;
+          if (rxn_by_proc[j] == -1) nghostlykeep[i]++;
-          else all_localskips[rxn_by_proc[j]]++;
+          else all_localkeep[rxn_by_proc[j]]++;
         }
         memory->destroy(rxn_by_proc);
         reaction_count_total[i] -= overstep;
       }
-      MPI_Scatter(&all_localskips[0],1,MPI_INT,&nlocalskips[i],1,MPI_INT,0,world);
+      MPI_Scatter(&all_localkeep[0],1,MPI_INT,&nlocalkeep[i],1,MPI_INT,0,world);
-      MPI_Bcast(&nghostlyskips[i],1,MPI_INT,0,world);
+      MPI_Bcast(&nghostlykeep[i],1,MPI_INT,0,world);
       memory->destroy(local_rxncounts);
-      memory->destroy(all_localskips);
+      memory->destroy(all_localkeep);
     }
   }
   MPI_Bcast(&reaction_count_total[0], nreacts, MPI_INT, 0, world);
+  memory->destroy(delta_rxn);
 
   // this updates topology next step
   next_reneighbor = update->ntimestep;
@@ -2965,6 +2974,8 @@ void FixBondReact::update_everything()
   int *type = atom->type;
   int **nspecial = atom->nspecial;
   tagint **special = atom->special;
+  tagint *tag = atom->tag;
+  AtomVec *avec = atom->avec;
 
   int **bond_type = atom->bond_type;
   tagint **bond_atom = atom->bond_atom;
@@ -2977,13 +2988,16 @@ void FixBondReact::update_everything()
   memory->create(mark,nmark,"bond/react:mark");
   for (int i = 0; i < nmark; i++) mark[i] = 0;
 
+  // used when creating atoms
+  addatomtag = 0;
+  for (int i = 0; i < nlocal; i++) addatomtag = MAX(addatomtag,tag[i]);
+  MPI_Allreduce(MPI_IN_PLACE,&addatomtag,1,MPI_LMP_TAGINT,MPI_MAX,world);
+  addatoms.clear();
+
   // flag used to delete special interactions
   int *delflag;
   memory->create(delflag,atom->maxspecial,"bond/react:delflag");
 
-  tagint *tag = atom->tag;
-  AtomVec *avec = atom->avec;
-
   // used when creating atoms
   int inserted_atoms_flag = 0;
 
@@ -3026,13 +3040,14 @@ void FixBondReact::update_everything()
 
   for (int pass = 0; pass < 2; pass++) {
     update_num_mega = 0;
-    int *iskip = new int[nreacts];
+    int *noccur = new int[nreacts];
-    for (int i = 0; i < nreacts; i++) iskip[i] = 0;
+    for (int i = 0; i < nreacts; i++) noccur[i] = 0;
     if (pass == 0) {
       for (int i = 0; i < local_num_mega; i++) {
         rxnID = (int) local_mega_glove[0][i];
         // reactions already shuffled from dedup procedure, so can skip first N
-        if (iskip[rxnID]++ < nlocalskips[rxnID]) continue;
+        // wait, this check needs to be after add atoms, because they can also be 'skipped' due to overlap
+        if (noccur[rxnID] >= nlocalkeep[rxnID]) continue;
 
         // this will be overwritten if reaction skipped by create_atoms below
         update_mega_glove[0][update_num_mega] = (tagint) local_mega_glove[0][i];
@@ -3043,13 +3058,14 @@ void FixBondReact::update_everything()
         if (create_atoms_flag[rxnID] == 1) {
           onemol = atom->molecules[unreacted_mol[rxnID]];
           twomol = atom->molecules[reacted_mol[rxnID]];
-          if (insert_atoms(update_mega_glove,update_num_mega)) {
+          if (insert_atoms_setup(update_mega_glove,update_num_mega)) {
             inserted_atoms_flag = 1;
           } else { // create aborted
             reaction_count_total[rxnID]--;
             continue;
           }
         }
+        noccur[rxnID]++;
 
         if (rescale_charges_flag[rxnID]) sim_total_charges[update_num_mega] = local_mega_glove[1][i];
         update_num_mega++;
@@ -3058,7 +3074,7 @@ void FixBondReact::update_everything()
       for (int i = 0; i < global_megasize; i++) {
         rxnID = (int) global_mega_glove[0][i];
         // reactions already shuffled from dedup procedure, so can skip first N
-        if (iskip[rxnID]++ < nghostlyskips[rxnID]) continue;
+        if (noccur[rxnID] >= nghostlykeep[rxnID]) continue;
 
         // this will be overwritten if reaction skipped by create_atoms below
         update_mega_glove[0][update_num_mega] = (tagint) global_mega_glove[0][i];
@@ -3071,30 +3087,49 @@ void FixBondReact::update_everything()
         if (create_atoms_flag[rxnID] == 1) {
           onemol = atom->molecules[unreacted_mol[rxnID]];
           twomol = atom->molecules[reacted_mol[rxnID]];
-          if (insert_atoms(update_mega_glove,update_num_mega)) {
+          if (insert_atoms_setup(update_mega_glove,update_num_mega)) {
             inserted_atoms_flag = 1;
           } else { // create aborted
             reaction_count_total[rxnID]--;
             continue;
           }
         }
+        noccur[rxnID]++;
 
         if (rescale_charges_flag[rxnID]) sim_total_charges[update_num_mega] = global_mega_glove[1][i];
         update_num_mega++;
       }
     }
-    delete [] iskip;
+    delete [] noccur;
 
     if (update_num_mega == 0) continue;
 
-    // if inserted atoms and global map exists, reset map now instead
+    // insert all atoms for all rxns here
-    //   of waiting for comm since other pre-exchange fixes may use it
+    if (inserted_atoms_flag == 1) {
-    // invoke map_init() b/c atom count has grown
+      // clear to-be-overwritten ghost info
-    // do this once after all atom insertions
+      atom->nghost = 0;
-    if (inserted_atoms_flag == 1 && atom->map_style != Atom::MAP_NONE) {
+      atom->avec->clear_bonus();
+
+      for (auto & myaddatom : addatoms) {
+        atom->avec->create_atom(myaddatom.type,myaddatom.x);
+        int n = atom->nlocal - 1;
+        atom->tag[n] = myaddatom.tag;
+        atom->molecule[n] = myaddatom.molecule;
+        atom->mask[n] = myaddatom.mask;
+        atom->image[n] = myaddatom.image;
+        atom->v[n][0] = myaddatom.v[0];
+        atom->v[n][1] = myaddatom.v[1];
+        atom->v[n][2] = myaddatom.v[2];
+        if (atom->rmass) atom->rmass[n]= myaddatom.rmass;
+        modify->create_attribute(n);
+      }
+
+      // reset atom->map
+      if (atom->map_style != Atom::MAP_NONE) {
         atom->map_init();
         atom->map_set();
       }
+    }
 
     // mark to-delete atoms
     nlocal = atom->nlocal;
@@ -3620,10 +3655,6 @@ void FixBondReact::update_everything()
 
   atom->natoms -= ndel;
   // done deleting atoms
-
-  // reset mol ids
-  if (reset_mol_ids_flag) reset_mol_ids->reset();
-
   // something to think about: this could done much more concisely if
   // all atom-level info (bond,angles, etc...) were kinda inherited from a common data struct --JG
 
@@ -3651,10 +3682,11 @@ void FixBondReact::update_everything()
 }
 
 /* ----------------------------------------------------------------------
-insert created atoms
+setup for inserting created atoms
+atoms for all rxns are actually created all at once in update_everything
 ------------------------------------------------------------------------- */
 
-int FixBondReact::insert_atoms(tagint **my_update_mega_glove, int iupdate)
+int FixBondReact::insert_atoms_setup(tagint **my_update_mega_glove, int iupdate)
 {
   // inserting atoms based off fix_deposit->pre_exchange
   int flag;
@@ -3682,14 +3714,9 @@ int FixBondReact::insert_atoms(tagint **my_update_mega_glove, int iupdate)
   tagint *molecule = atom->molecule;
   int nlocal = atom->nlocal;
 
-  tagint maxtag_all,maxmol_all;
+  tagint maxmol_all;
-  tagint max = 0;
+  for (int i = 0; i < nlocal; i++) maxmol_all = MAX(maxmol_all,molecule[i]);
-  for (int i = 0; i < nlocal; i++) max = MAX(max,tag[i]);
+  MPI_Allreduce(MPI_IN_PLACE,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
-  MPI_Allreduce(&max,&maxtag_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
-
-  max = 0;
-  for (int i = 0; i < nlocal; i++) max = MAX(max,molecule[i]);
-  MPI_Allreduce(&max,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
 
   int dimension = domain->dimension;
 
@@ -3786,6 +3813,26 @@ int FixBondReact::insert_atoms(tagint **my_update_mega_glove, int iupdate)
         if (abortflag) break;
       }
     }
+    // also check against previous to-be-added atoms
+    if (!abortflag) {
+      for (auto & myaddatom : addatoms) {
+        for (int m = 0; m < twomol->natoms; m++) {
+          if (create_atoms[m][rxnID] == 1) {
+            delx = coords[m][0] - myaddatom.x[0];
+            dely = coords[m][1] - myaddatom.x[1];
+            delz = coords[m][2] - myaddatom.x[2];
+            domain->minimum_image(delx,dely,delz);
+            rsq = delx*delx + dely*dely + delz*delz;
+            if (rsq < overlapsq[rxnID]) {
+              abortflag = 1;
+              break;
+            }
+          }
+        }
+        if (abortflag) break;
+      }
+    }
+
     MPI_Allreduce(MPI_IN_PLACE,&abortflag,1,MPI_INT,MPI_MAX,world);
     if (abortflag) {
       memory->destroy(coords);
@@ -3794,12 +3841,6 @@ int FixBondReact::insert_atoms(tagint **my_update_mega_glove, int iupdate)
     }
   }
 
-  // clear ghost count and any ghost bonus data internal to AtomVec
-  // same logic as beginning of Comm::exchange()
-  // do it now b/c inserting atoms will overwrite ghost atoms
-  atom->nghost = 0;
-  atom->avec->clear_bonus();
-
   // check if new atoms are in my sub-box or above it if I am highest proc
   // if so, add atom to my list via create_atom()
   // initialize additional info about the atoms
@@ -3842,40 +3883,46 @@ int FixBondReact::insert_atoms(tagint **my_update_mega_glove, int iupdate)
       }
 
       int root = 0;
+      addatomtag++;
       if (flag) {
+        struct AddAtom myaddatom;
         root = comm->me;
 
-        atom->avec->create_atom(twomol->type[m],coords[m]);
+        myaddatom.type = twomol->type[m];
-        int n = atom->nlocal - 1;
+        myaddatom.x[0] = coords[m][0];
-        atom->tag[n] = maxtag_all + add_count;
+        myaddatom.x[1] = coords[m][1];
+        myaddatom.x[2] = coords[m][2];
+        myaddatom.tag = addatomtag;
 
         // locally update mega_glove
-        my_update_mega_glove[preID][iupdate] = atom->tag[n];
+        my_update_mega_glove[preID][iupdate] = myaddatom.tag;
 
+        // !! could do better job choosing mol ID for added atoms
         if (atom->molecule_flag) {
           if (twomol->moleculeflag) {
-            atom->molecule[n] = maxmol_all + twomol->molecule[m];
+            myaddatom.molecule = maxmol_all + twomol->molecule[m];
           } else {
-            atom->molecule[n] = maxmol_all + 1;
+            myaddatom.molecule = maxmol_all + 1;
           }
         }
 
-        atom->mask[n] = 1 | groupbit;
+        myaddatom.mask = 1 | groupbit;
-        atom->image[n] = imageflags[m];
+        myaddatom.image = imageflags[m];
 
         // guess a somewhat reasonable initial velocity based on reaction site
         // further control is possible using bond_react_MASTER_group
         // compute |velocity| corresponding to a given temperature t, using specific atom's mass
-        double mymass = atom->rmass ? atom->rmass[n] : atom->mass[twomol->type[m]];
+        myaddatom.rmass = atom->rmass ? twomol->rmass[m] : atom->mass[twomol->type[m]];
-        double vtnorm = sqrt(t / (force->mvv2e / (dimension * force->boltz)) / mymass);
+        double vtnorm = sqrt(t / (force->mvv2e / (dimension * force->boltz)) / myaddatom.rmass);
-        v[n][0] = random[rxnID]->uniform();
+        double myv[3];
-        v[n][1] = random[rxnID]->uniform();
+        myv[0] = random[rxnID]->uniform();
-        v[n][2] = random[rxnID]->uniform();
+        myv[1] = random[rxnID]->uniform();
-        double vnorm = sqrt(v[n][0]*v[n][0] + v[n][1]*v[n][1] + v[n][2]*v[n][2]);
+        myv[2] = random[rxnID]->uniform();
-        v[n][0] = v[n][0]/vnorm*vtnorm;
+        double vnorm = sqrt(myv[0]*myv[0] + myv[1]*myv[1] + myv[2]*myv[2]);
-        v[n][1] = v[n][1]/vnorm*vtnorm;
+        myaddatom.v[0] = myv[0]/vnorm*vtnorm;
-        v[n][2] = v[n][2]/vnorm*vtnorm;
+        myaddatom.v[1] = myv[1]/vnorm*vtnorm;
-        modify->create_attribute(n);
+        myaddatom.v[2] = myv[2]/vnorm*vtnorm;
+        addatoms.push_back(myaddatom);
       }
       // globally update mega_glove and equivalences
       MPI_Allreduce(MPI_IN_PLACE,&root,1,MPI_INT,MPI_SUM,world);
@@ -3888,12 +3935,11 @@ int FixBondReact::insert_atoms(tagint **my_update_mega_glove, int iupdate)
   }
 
   // reset global natoms here
-  // reset atom map elsewhere, after all calls to 'insert_atoms'
+  // reset atom map elsewhere, after all calls to 'insert_atoms_setup'
   atom->natoms += add_count;
   if (atom->natoms < 0)
     error->all(FLERR,"Too many total atoms");
-  maxtag_all += add_count;
+  if (addatomtag >= MAXTAGINT)
-  if (maxtag_all >= MAXTAGINT)
     error->all(FLERR,"New atom IDs exceed maximum allowed ID");
   // atom creation successful
   memory->destroy(coords);
@@ -3970,6 +4016,11 @@ void FixBondReact::read_map_file(int myrxn)
     } else break;
   }
 
+  if (ncreate == 0 && onemol->natoms != twomol->natoms)
+    error->all(FLERR,"Fix bond/react: Reaction templates must contain the same number of atoms");
+  else if (ncreate > 0 && onemol->natoms + ncreate != twomol->natoms)
+    error->all(FLERR,"Fix bond/react: Incorrect number of created atoms");
+
   // grab keyword and skip next line
 
   parse_keyword(0,line,keyword);
@@ -4012,6 +4063,13 @@ void FixBondReact::read_map_file(int myrxn)
 
   }
 
+  // error check
+  for (int i = 0; i < onemol->natoms; i++) {
+    int my_equiv = reverse_equiv[i][1][myrxn];
+    if (create_atoms[my_equiv-1][myrxn] == 1)
+      error->all(FLERR,"Fix bond/react: Created atoms cannot also be listed in Equivalences section\n");
+  }
+
   // error check
   if (bondflag == 0 || equivflag == 0)
     error->all(FLERR,"Fix bond/react: Map file missing InitiatorIDs or Equivalences section\n");
@@ -4071,6 +4129,8 @@ void FixBondReact::CreateAtoms(char *line, int myrxn)
     readline(line);
     rv = sscanf(line,"%d",&tmp);
     if (rv != 1) error->one(FLERR, "CreateIDs section is incorrectly formatted");
+    if (tmp > twomol->natoms)
+      error->one(FLERR,"Fix bond/react: Invalid atom ID in CreateIDs section of map file");
     create_atoms[tmp-1][myrxn] = 1;
   }
   if (twomol->xflag == 0)
@@ -4331,6 +4391,13 @@ void FixBondReact::post_integrate_respa(int ilevel, int /*iloop*/)
 
 /* ---------------------------------------------------------------------- */
 
+void FixBondReact::post_force(int /*vflag*/)
+{
+  if (reset_mol_ids_flag) reset_mol_ids->reset();
+}
+
+/* ---------------------------------------------------------------------- */
+
 int FixBondReact::pack_forward_comm(int n, int *list, double *buf,
                                     int /*pbc_flag*/, int * /*pbc*/)
 {

src/REACTION/fix_bond_react.h

@@ -46,6 +46,7 @@ class FixBondReact : public Fix {
   void init_list(int, class NeighList *) override;
   void post_integrate() override;
   void post_integrate_respa(int, int) override;
+  void post_force(int) override;
 
   int pack_forward_comm(int, int *, double *, int, int *) override;
   void unpack_forward_comm(int, int, double *) override;
@@ -62,7 +63,7 @@ class FixBondReact : public Fix {
   int *iatomtype, *jatomtype;
   int *seed;
   double **cutsq, *fraction;
-  int *max_rxn, *nlocalskips, *nghostlyskips;
+  int *max_rxn, *nlocalkeep, *nghostlykeep;
   tagint lastcheck;
   int stabilization_flag;
   int reset_mol_ids_flag;
@@ -215,7 +216,7 @@ class FixBondReact : public Fix {
   void glove_ghostcheck();
   void ghost_glovecast();
   void update_everything();
-  int insert_atoms(tagint **, int);
+  int insert_atoms_setup(tagint **, int);
   void unlimit_bond(); // removes atoms from stabilization, and other post-reaction every-step operations
   void dedup_mega_gloves(int);    //dedup global mega_glove
   void write_restart(FILE *) override;
@@ -245,6 +246,15 @@ class FixBondReact : public Fix {
   std::map<std::set<tagint>, int> atoms2bond;    // maps atom pair to index of local bond array
   std::vector<std::vector<Constraint>> constraints;
 
+  tagint addatomtag;
+  struct AddAtom {
+    tagint tag, molecule;
+    int type, mask;
+    imageint image;
+    double rmass, x[3], v[3];
+  };
+  std::vector<AddAtom> addatoms;
+
   // DEBUG
 
   void print_bb();