Merge pull request #3905 from jrgissing/bond_react_fixes_aug23

Assorted bond/react fixes

examples/PACKAGES/reaction/create_atoms_polystyrene/in.grow_styrene

@@ -40,7 +40,7 @@ fix 1 statted_grp_REACT nvt temp $T $T 100
 
 fix 4 bond_react_MASTER_group temp/rescale 1 $T $T 1 1
 
-thermo_style custom step temp press density f_myrxns[1]
+thermo_style custom step temp press density f_myrxns[*]
 
 thermo 100
 

examples/PACKAGES/reaction/nylon,6-6_melt/in.large_nylon_melt

@@ -50,7 +50,7 @@ fix 1 statted_grp_REACT nvt temp 800 800 100
 # you can use the internally created 'bond_react_MASTER_group', like so:
 # fix 2 bond_react_MASTER_group temp/rescale 1 800 800 10 1
 
-thermo_style custom step temp press density f_myrxns[1] f_myrxns[2] # cumulative reaction counts
+thermo_style custom step temp press density f_myrxns[*] # cumulative reaction counts
 
 # restart 100 restart1 restart2
 

examples/PACKAGES/reaction/tiny_epoxy/in.tiny_epoxy.stabilized

@@ -20,7 +20,8 @@ improper_style class2
 special_bonds lj/coul 0 0 1
 pair_modify tail yes mix sixthpower
 
-read_data tiny_epoxy.data
+read_data tiny_epoxy.data &
+  extra/special/per/atom 25
 
 velocity all create 300.0 4928459 dist gaussian
 
@@ -44,7 +45,7 @@ fix rxns all bond/react stabilization yes statted_grp .03 &
 
 fix 1 statted_grp_REACT nvt temp 300 300 100
 
-thermo_style custom step temp f_rxns[1] f_rxns[2] f_rxns[3] f_rxns[4]
+thermo_style custom step temp f_rxns[*]
 
 run 2000
 

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

@@ -50,7 +50,7 @@ fix 1 statted_grp_REACT nvt temp 300 300 100
 # by using the internally-created 'bond_react_MASTER_group', like so:
 fix 4 bond_react_MASTER_group temp/rescale 1 300 300 10 1
 
-thermo_style custom step temp press density f_myrxns[1] f_myrxns[2]
+thermo_style custom step temp press density f_myrxns[*]
 
 # restart 100 restart1 restart2
 

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

@@ -54,7 +54,7 @@ fix 1 statted_grp_REACT nvt temp 300 300 100
 # by using the internally-created 'bond_react_MASTER_group', like so:
 fix 4 bond_react_MASTER_group temp/rescale 1 300 300 10 1
 
-thermo_style custom step temp press density v_prob1 v_prob2 f_myrxns[1] f_myrxns[2]
+thermo_style custom step temp press density v_prob1 v_prob2 f_myrxns[*]
 
 # restart 100 restart1 restart2
 

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

@@ -47,7 +47,7 @@ fix myrxns all bond/react stabilization no &
 
 fix 1 all nve/limit .03
 
-thermo_style custom step temp press density f_myrxns[1] f_myrxns[2]
+thermo_style custom step temp press density f_myrxns[*]
 
 # restart 100 restart1 restart2
 

examples/PACKAGES/reaction/tiny_polystyrene/in.tiny_polystyrene.stabilized

@@ -51,7 +51,7 @@ fix 1 statted_grp_REACT nvt temp $T $T 100
 
 fix 4 bond_react_MASTER_group temp/rescale 1 $T $T 1 1
 
-thermo_style custom step temp press density f_rxn1[1] f_rxn1[2] f_rxn1[3]
+thermo_style custom step temp press density f_rxn1[*]
 
 run 10000
 

src/REACTION/README

@@ -25,4 +25,5 @@ The REACTER methodology is detailed in:
     https://doi.org/10.1021/acs.macromol.0c02012
 
 This package was created by Jacob Gissinger
-(jacob.r.gissinger@gmail.com) at the NASA Langley Research Center.
+(jgissing@stevens.edu) while at the NASA Langley Research Center
+and Stevens Institute of Technology.

src/REACTION/fix_bond_react.cpp

@@ -13,7 +13,7 @@ See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
-Contributing Author: Jacob Gissinger (jacob.r.gissinger@gmail.com)
+Contributing Author: Jacob Gissinger (jgissing@stevens.edu)
 ------------------------------------------------------------------------- */
 
 #include "fix_bond_react.h"
@@ -670,15 +670,6 @@ FixBondReact::~FixBondReact()
   memory->destroy(ghostly_rxn_count);
   memory->destroy(reaction_count_total);
 
-  if (newton_bond == 0) {
-    memory->destroy(xspecial);
-    memory->destroy(nxspecial);
-    memory->destroy(onemol_xspecial);
-    memory->destroy(onemol_nxspecial);
-    memory->destroy(twomol_xspecial);
-    memory->destroy(twomol_nxspecial);
-  }
-
   if (attempted_rxn == 1) {
     memory->destroy(restore_pt);
     memory->destroy(restore);
@@ -827,11 +818,10 @@ void FixBondReact::init()
     nlevels_respa = (dynamic_cast<Respa *>(update->integrate))->nlevels;
 
   // check cutoff for iatomtype,jatomtype
-  for (int i = 0; i < nreacts; i++) {
   if (!utils::strmatch(force->pair_style,"^hybrid"))
-      if (force->pair == nullptr || cutsq[i][1] > force->pair->cutsq[iatomtype[i]][jatomtype[i]])
+    for (int i = 0; i < nreacts; i++)
+      if (force->pair == nullptr || (closeneigh[i] < 0 && cutsq[i][1] > force->pair->cutsq[iatomtype[i]][jatomtype[i]]))
         error->all(FLERR,"Fix bond/react: Fix bond/react cutoff is longer than pairwise cutoff");
-  }
 
   // need a half neighbor list, built every Nevery steps
   neighbor->add_request(this, NeighConst::REQ_OCCASIONAL);
@@ -931,29 +921,10 @@ void FixBondReact::post_integrate()
 
   neighbor->build_one(list,1);
 
-  // here we define a full special list, independent of Newton setting
+  // here we define a full special list
-  if (newton_bond == 1) {
+  // may need correction for unusual special bond settings
   nxspecial = atom->nspecial;
   xspecial = atom->special;
-  } else {
-    int nall = atom->nlocal + atom->nghost;
-    memory->destroy(nxspecial);
-    memory->destroy(xspecial);
-    memory->create(nxspecial,nall,3,"bond/react:nxspecial");
-    memory->create(xspecial,nall,atom->maxspecial,"bond/react:xspecial");
-    for (int i = 0; i < atom->nlocal; i++) {
-      nxspecial[i][0] = atom->num_bond[i];
-      for (int j = 0; j < nxspecial[i][0]; j++) {
-        xspecial[i][j] = atom->bond_atom[i][j];
-      }
-      nxspecial[i][1] = atom->nspecial[i][1];
-      nxspecial[i][2] = atom->nspecial[i][2];
-      int joffset = nxspecial[i][0] - atom->nspecial[i][0];
-      for (int j = nxspecial[i][0]; j < nxspecial[i][2]; j++) {
-        xspecial[i][j+joffset] = atom->special[i][j];
-      }
-    }
-  }
 
   int j;
   for (rxnID = 0; rxnID < nreacts; rxnID++) {
@@ -2541,49 +2512,15 @@ int FixBondReact::get_chirality(double four_coords[12])
 
 /* ----------------------------------------------------------------------
   Get xspecials for current molecule templates
+  may need correction when specials defined explicitly in molecule templates
 ------------------------------------------------------------------------- */
 
 void FixBondReact::get_molxspecials()
 {
-  if (newton_bond == 1) {
   onemol_nxspecial = onemol->nspecial;
   onemol_xspecial = onemol->special;
   twomol_nxspecial = twomol->nspecial;
   twomol_xspecial = twomol->special;
-  } else {
-    memory->destroy(onemol_nxspecial);
-    memory->destroy(onemol_xspecial);
-    memory->create(onemol_nxspecial,onemol->natoms,3,"bond/react:onemol_nxspecial");
-    memory->create(onemol_xspecial,onemol->natoms,atom->maxspecial,"bond/react:onemol_xspecial");
-    for (int i = 0; i < onemol->natoms; i++) {
-      onemol_nxspecial[i][0] = onemol->num_bond[i];
-      for (int j = 0; j < onemol_nxspecial[i][0]; j++) {
-        onemol_xspecial[i][j] = onemol->bond_atom[i][j];
-      }
-      onemol_nxspecial[i][1] = onemol->nspecial[i][1];
-      onemol_nxspecial[i][2] = onemol->nspecial[i][2];
-      int joffset = onemol_nxspecial[i][0] - onemol->nspecial[i][0];
-      for (int j = onemol_nxspecial[i][0]; j < onemol_nxspecial[i][2]; j++) {
-        onemol_xspecial[i][j+joffset] = onemol->special[i][j];
-      }
-    }
-    memory->destroy(twomol_nxspecial);
-    memory->destroy(twomol_xspecial);
-    memory->create(twomol_nxspecial,twomol->natoms,3,"bond/react:twomol_nxspecial");
-    memory->create(twomol_xspecial,twomol->natoms,atom->maxspecial,"bond/react:twomol_xspecial");
-    for (int i = 0; i < twomol->natoms; i++) {
-      twomol_nxspecial[i][0] = twomol->num_bond[i];
-      for (int j = 0; j < twomol_nxspecial[i][0]; j++) {
-        twomol_xspecial[i][j] = twomol->bond_atom[i][j];
-      }
-      twomol_nxspecial[i][1] = twomol->nspecial[i][1];
-      twomol_nxspecial[i][2] = twomol->nspecial[i][2];
-      int joffset = twomol_nxspecial[i][0] - twomol->nspecial[i][0];
-      for (int j = twomol_nxspecial[i][0]; j < twomol_nxspecial[i][2]; j++) {
-        twomol_xspecial[i][j+joffset] = twomol->special[i][j];
-      }
-    }
-  }
 }
 
 /* ----------------------------------------------------------------------
@@ -2682,15 +2619,42 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
   }
 
   // also, if atoms change number of bonds, but aren't landlocked, that could be bad
+  int warnflag = 0;
   if (comm->me == 0)
     for (int i = 0; i < twomol->natoms; i++) {
       if ((create_atoms[i][myrxn] == 0) &&
           (twomol_nxspecial[i][0] != onemol_nxspecial[equivalences[i][1][myrxn]-1][0]) &&
-          (landlocked_atoms[i][myrxn] == 0))
+          (landlocked_atoms[i][myrxn] == 0)) {
-        error->warning(FLERR, "Fix bond/react: Atom affected by reaction {} is too close "
+        warnflag = 1;
-                       "to template edge",rxn_name[myrxn]);
         break;
       }
+    }
+
+  // also, if an atom changes any of its bonds, but is not landlocked, that could be bad
+  int thereflag;
+  if (comm->me == 0)
+    for (int i = 0; i < twomol->natoms; i++) {
+      if (landlocked_atoms[i][myrxn] == 1) continue;
+      for (int j = 0; j < twomol_nxspecial[i][0]; j++) {
+        int oneneighID = equivalences[twomol_xspecial[i][j]-1][1][myrxn];
+        int ii = equivalences[i][1][myrxn] - 1;
+        thereflag = 0;
+        for (int k = 0; k < onemol_nxspecial[ii][0]; k++) {
+          if (oneneighID == onemol_xspecial[ii][k]) {
+            thereflag = 1;
+            break;
+          }
+        }
+        if (thereflag == 0) {
+          warnflag = 1;
+          break;
+        }
+      }
+      if (warnflag == 1) break;
+    }
+
+  if (comm->me == 0 && warnflag == 1) error->warning(FLERR, "Fix bond/react: Atom affected "
+                       "by reaction {} is too close to template edge",rxn_name[myrxn]);
 
   // finally, if a created atom is not landlocked, bad!
   for (int i = 0; i < twomol->natoms; i++) {
@@ -3349,7 +3313,7 @@ void FixBondReact::update_everything()
         dynamic_cast<FixBondHistory *>(ihistory)->clear_cache();
 
     // Angles! First let's delete all angle info:
-    if (force->angle && twomol->angleflag) {
+    if (force->angle) {
       int *num_angle = atom->num_angle;
       int **angle_type = atom->angle_type;
       tagint **angle_atom1 = atom->angle_atom1;
@@ -3390,6 +3354,7 @@ void FixBondReact::update_everything()
           }
         }
         // now let's add the new angle info.
+        if (twomol->angleflag) {
           for (int j = 0; j < twomol->natoms; j++) {
             int jj = equivalences[j][1][rxnID]-1;
             if (atom->map(update_mega_glove[jj+1][i]) < nlocal && atom->map(update_mega_glove[jj+1][i]) >= 0) {
@@ -3424,9 +3389,10 @@ void FixBondReact::update_everything()
           }
         }
       }
+    }
 
     // Dihedrals! first let's delete all dihedral info for landlocked atoms
-    if (force->dihedral && twomol->dihedralflag) {
+    if (force->dihedral) {
       int *num_dihedral = atom->num_dihedral;
       int **dihedral_type = atom->dihedral_type;
       tagint **dihedral_atom1 = atom->dihedral_atom1;
@@ -3470,6 +3436,7 @@ void FixBondReact::update_everything()
           }
         }
         // now let's add new dihedral info
+        if (twomol->dihedralflag) {
           for (int j = 0; j < twomol->natoms; j++) {
             int jj = equivalences[j][1][rxnID]-1;
             if (atom->map(update_mega_glove[jj+1][i]) < nlocal && atom->map(update_mega_glove[jj+1][i]) >= 0) {
@@ -3507,9 +3474,10 @@ void FixBondReact::update_everything()
           }
         }
       }
+    }
 
     // finally IMPROPERS!!!! first let's delete all improper info for landlocked atoms
-    if (force->improper && twomol->improperflag) {
+    if (force->improper) {
       int *num_improper = atom->num_improper;
       int **improper_type = atom->improper_type;
       tagint **improper_atom1 = atom->improper_atom1;
@@ -3553,6 +3521,7 @@ void FixBondReact::update_everything()
           }
         }
         // now let's add new improper info
+        if (twomol->improperflag) {
           for (int j = 0; j < twomol->natoms; j++) {
             int jj = equivalences[j][1][rxnID]-1;
             if (atom->map(update_mega_glove[jj+1][i]) < nlocal && atom->map(update_mega_glove[jj+1][i]) >= 0) {
@@ -3590,6 +3559,7 @@ void FixBondReact::update_everything()
           }
         }
       }
+    }
 
   }
 
@@ -3895,7 +3865,8 @@ int FixBondReact::insert_atoms(tagint **my_update_mega_glove, int iupdate)
         // 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 vtnorm = sqrt(t / (force->mvv2e / (dimension * force->boltz)) / atom->mass[twomol->type[m]]);
+        double mymass = atom->rmass ? atom->rmass[n] : atom->mass[twomol->type[m]];
+        double vtnorm = sqrt(t / (force->mvv2e / (dimension * force->boltz)) / mymass);
         v[n][0] = random[rxnID]->uniform();
         v[n][1] = random[rxnID]->uniform();
         v[n][2] = random[rxnID]->uniform();

src/REACTION/fix_bond_react.h

@@ -12,7 +12,7 @@
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
-   Contributing Author: Jacob Gissinger (jacob.r.gissinger@gmail.com)
+   Contributing Author: Jacob Gissinger (jgissing@stevens.edu)
 ------------------------------------------------------------------------- */
 
 #ifdef FIX_CLASS
@@ -139,7 +139,7 @@ class FixBondReact : public Fix {
   int avail_guesses;     // num of restore points available
   int *guess_branch;     // used when there is more than two choices when guessing
   int **restore_pt;      // contains info about restore points
-  tagint **restore;      // contaings info about restore points
+  tagint **restore;      // contains info about restore points
   int *pioneer_count;    // counts pioneers
 
   int **edge;                // atoms in molecule templates with incorrect valences