allow wildcards with "cutoff" keyword to fix reaxff/species

this also switched to using fmtlib for column aligned output formatting
and re-applies clang-format.

doc/src/fix_reaxff_species.rst

@@ -25,7 +25,7 @@ Syntax
   .. parsed-literal::
 
        *cutoff* value = I J Cutoff
-         I, J = atom types
+         I, J = atom types (see asterisk form below)
          Cutoff = Bond-order cutoff value for this pair of atom types
        *element* value = Element1, Element2, ...
        *position* value = posfreq filepos
@@ -49,7 +49,7 @@ Examples
 .. code-block:: LAMMPS
 
    fix 1 all reaxff/species 10 10 100 species.out
-   fix 1 all reaxff/species 1 2 20 species.out cutoff 1 1 0.40 cutoff 1 2 0.55
+   fix 1 all reaxff/species 1 2 20 species.out cutoff 1 1 0.40 cutoff 1 2*3 0.55
    fix 1 all reaxff/species 1 100 100 species.out element Au O H position 1000 AuOH.pos
    fix 1 all reaxff/species 1 100 100 species.out delete species.del masslimit 0 50
 
@@ -88,13 +88,24 @@ If the filename ends with ".gz", the output file is written in gzipped
 format.  A gzipped dump file will be about 3x smaller than the text version,
 but will also take longer to write.
 
+.. versionadded:: TBD
+
+   Support for wildcards added
+
 Optional keyword *cutoff* can be assigned to change the minimum
 bond-order values used in identifying chemical bonds between pairs of
 atoms.  Bond-order cutoffs should be carefully chosen, as bond-order
-cutoffs that are too small may include too many bonds (which will
+cutoffs that are too small may include too many bonds (which will result
-result in an error), while cutoffs that are too large will result in
+in an error), while cutoffs that are too large will result in fragmented
-fragmented molecules.  The default cutoff of 0.3 usually gives good
+molecules.  The default cutoff of 0.3 usually gives good results.  A
-results.
+wildcard asterisk can be used in place of or in conjunction with the I,J
+arguments to set the bond-order cutoff for multiple pairs of atom types.
+This takes the form "\*" or "\*n" or "n\*" or "m\*n".  If :math:`N` is
+the number of atom types, then an asterisk with no numeric values means
+all types from 1 to :math:`N`.  A leading asterisk means all types from
+1 to n (inclusive).  A trailing asterisk means all types from n to
+:math:`N` (inclusive).  A middle asterisk means all types from m to n
+(inclusive).
 
 The optional keyword *element* can be used to specify the chemical
 symbol printed for each LAMMPS atom type. The number of symbols must

src/REAXFF/fix_reaxff_species.cpp

@@ -48,15 +48,16 @@ using namespace LAMMPS_NS;
 using namespace FixConst;
 
 static const char cite_reaxff_species_delete[] =
-  "fix reaxff/species, 'delete' keyword: https://doi.org/10.1016/j.carbon.2022.11.002\n\n"
+    "fix reaxff/species, 'delete' keyword: https://doi.org/10.1016/j.carbon.2022.11.002\n\n"
-  "@Article{Gissinger23,\n"
+    "@Article{Gissinger23,\n"
-  " author = {J. R. Gissinger, S. R. Zavada, J. G. Smith, J. Kemppainen, I. Gallegos, G. M. Odegard, E. J. Siochi, K. E. Wise},\n"
+    " author = {J. R. Gissinger, S. R. Zavada, J. G. Smith, J. Kemppainen, I. Gallegos, G. M. "
-  " title = {Predicting char yield of high-temperature resins},\n"
+    "Odegard, E. J. Siochi, K. E. Wise},\n"
-  " journal = {Carbon},\n"
+    " title = {Predicting char yield of high-temperature resins},\n"
-  " year =    2023,\n"
+    " journal = {Carbon},\n"
-  " volume =  202,\n"
+    " year =    2023,\n"
-  " pages =   {336-347}\n"
+    " volume =  202,\n"
-  "}\n\n";
+    " pages =   {336-347}\n"
+    "}\n\n";
 
 /* ---------------------------------------------------------------------- */
 
@@ -148,13 +149,11 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) :
   setupflag = 0;
 
   // set default bond order cutoff
-  int itype, jtype;
+  double bo_cut = 0.30;
-  double bo_cut;
-  bg_cut = 0.30;
   int np1 = ntypes + 1;
   memory->create(BOCut, np1, np1, "reaxff/species:BOCut");
   for (int i = 1; i < np1; i++)
-    for (int j = 1; j < np1; j++) BOCut[i][j] = bg_cut;
+    for (int j = 1; j < np1; j++) BOCut[i][j] = bo_cut;
 
   // optional args
   eletype = nullptr;
@@ -172,16 +171,19 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) :
     // set BO cutoff
     if (strcmp(arg[iarg], "cutoff") == 0) {
       if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix reaxff/species cutoff", error);
-      itype = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
+      int ilo, ihi, jlo, jhi;
-      jtype = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
+      utils::bounds(FLERR, arg[iarg + 1], 1, atom->ntypes, ilo, ihi, error);
+      utils::bounds(FLERR, arg[iarg + 2], 1, atom->ntypes, jlo, jhi, error);
       bo_cut = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
-      if ((itype <= 0) || (jtype <= 0) || (itype > ntypes) || (jtype > ntypes))
-        error->all(FLERR, "Fix reaxff/species cutoff atom type(s) out of range");
       if ((bo_cut > 1.0) || (bo_cut < 0.0))
         error->all(FLERR, "Fix reaxff/species invalid cutoff value: {}", bo_cut);
 
-      BOCut[itype][jtype] = bo_cut;
+      for (int i = ilo; i <= ihi; ++i) {
-      BOCut[jtype][itype] = bo_cut;
+        for (int j = MAX(jlo, i); j <= jhi; ++j) {
+          BOCut[i][j] = bo_cut;
+          BOCut[j][i] = bo_cut;
+        }
+      }
       iarg += 4;
 
       // modify element type names
@@ -240,17 +242,21 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) :
         error->all(FLERR, "Unknown fix reaxff/species delete option: {}", arg[iarg]);
       // rate limit when deleting molecules
     } else if (strcmp(arg[iarg], "delete_rate_limit") == 0) {
-      if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix reaxff/species delete_rate_limit", error);
+      if (iarg + 3 > narg)
+        utils::missing_cmd_args(FLERR, "fix reaxff/species delete_rate_limit", error);
       delete_Nlimit_varid = -1;
-      if (strncmp(arg[iarg+1],"v_",2) == 0) {
+      if (strncmp(arg[iarg + 1], "v_", 2) == 0) {
-        delete_Nlimit_varname = &arg[iarg+1][2];
+        delete_Nlimit_varname = &arg[iarg + 1][2];
         delete_Nlimit_varid = input->variable->find(delete_Nlimit_varname.c_str());
         if (delete_Nlimit_varid < 0)
-          error->all(FLERR,"Fix reaxff/species: Variable name {} does not exist",delete_Nlimit_varname);
+          error->all(FLERR, "Fix reaxff/species: Variable name {} does not exist",
+                     delete_Nlimit_varname);
         if (!input->variable->equalstyle(delete_Nlimit_varid))
-          error->all(FLERR,"Fix reaxff/species: Variable {} is not equal-style",delete_Nlimit_varname);
+          error->all(FLERR, "Fix reaxff/species: Variable {} is not equal-style",
-      } else delete_Nlimit = utils::numeric(FLERR, arg[iarg+1], false, lmp);
+                     delete_Nlimit_varname);
-      delete_Nsteps = utils::numeric(FLERR, arg[iarg+2], false, lmp);
+      } else
+        delete_Nlimit = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
+      delete_Nsteps = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
       iarg += 3;
       // position of molecules
     } else if (strcmp(arg[iarg], "position") == 0) {
@@ -292,10 +298,9 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) :
 
   if (delete_Nsteps > 0) {
     if (lmp->citeme) lmp->citeme->add(cite_reaxff_species_delete);
-    memory->create(delete_Tcount,delete_Nsteps,"reaxff/species:delete_Tcount");
+    memory->create(delete_Tcount, delete_Nsteps, "reaxff/species:delete_Tcount");
 
-    for (int i = 0; i < delete_Nsteps; i++)
+    for (int i = 0; i < delete_Nsteps; i++) delete_Tcount[i] = -1;
-      delete_Tcount[i] = -1;
     delete_Tcount[0] = 0;
   }
 
@@ -393,9 +398,11 @@ void FixReaxFFSpecies::init()
   if (delete_Nsteps > 0) {
     delete_Nlimit_varid = input->variable->find(delete_Nlimit_varname.c_str());
     if (delete_Nlimit_varid < 0)
-      error->all(FLERR,"Fix reaxff/species: Variable name {} does not exist",delete_Nlimit_varname);
+      error->all(FLERR, "Fix reaxff/species: Variable name {} does not exist",
+                 delete_Nlimit_varname);
     if (!input->variable->equalstyle(delete_Nlimit_varid))
-      error->all(FLERR,"Fix reaxff/species: Variable {} is not equal-style",delete_Nlimit_varname);
+      error->all(FLERR, "Fix reaxff/species: Variable {} is not equal-style",
+                 delete_Nlimit_varname);
   }
 }
 
@@ -427,8 +434,7 @@ void FixReaxFFSpecies::Output_ReaxFF_Bonds(bigint ntimestep, FILE * /*fp*/)
   if (ntimestep != nvalid) {
     // push back delete_Tcount on every step
     if (delete_Nsteps > 0)
-      for (int i = delete_Nsteps-1; i > 0; i--)
+      for (int i = delete_Nsteps - 1; i > 0; i--) delete_Tcount[i] = delete_Tcount[i - 1];
-        delete_Tcount[i] = delete_Tcount[i-1];
     return;
   }
 
@@ -732,31 +738,31 @@ void FixReaxFFSpecies::WriteFormulas(int Nmole, int Nspec)
   int i, j, itemp;
   bigint ntimestep = update->ntimestep;
 
-  fprintf(fp, "# Timestep     No_Moles     No_Specs     ");
+  fprintf(fp, "#  Timestep    No_Moles    No_Specs");
 
   Nmoltype = 0;
 
   for (i = 0; i < Nspec; i++) nd[i] = CheckExistence(i, ntypes);
 
   for (i = 0; i < Nmoltype; i++) {
+    std::string molname;
     for (j = 0; j < ntypes; j++) {
       itemp = MolType[ntypes * i + j];
       if (itemp != 0) {
         if (eletype)
-          fprintf(fp, "%s", eletype[j]);
+          molname += eletype[j];
         else
-          fprintf(fp, "%c", ele[j]);
+          molname += ele[j];
-        if (itemp != 1) fprintf(fp, "%d", itemp);
+        if (itemp != 1) molname += std::to_string(itemp);
       }
     }
-    fputs("\t", fp);
+    fmt::print(fp, " {:>11}", molname);
   }
   fputs("\n", fp);
 
-  fmt::print(fp, "{} {:11} {:11}\t", ntimestep, Nmole, Nspec);
+  fmt::print(fp, "{:>11} {:>11} {:>11}", ntimestep, Nmole, Nspec);
-
+  for (i = 0; i < Nmoltype; i++) fmt::print(fp, " {:>11}", NMol[i]);
-  for (i = 0; i < Nmoltype; i++) fprintf(fp, " %d\t", NMol[i]);
+  fputs("\n", fp);
-  fprintf(fp, "\n");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -884,8 +890,8 @@ void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec)
   int ndeletions;
   int headroom = -1;
   if (delete_Nsteps > 0) {
-    if (delete_Tcount[delete_Nsteps-1] == -1) return;
+    if (delete_Tcount[delete_Nsteps - 1] == -1) return;
-    ndeletions = delete_Tcount[0] - delete_Tcount[delete_Nsteps-1];
+    ndeletions = delete_Tcount[0] - delete_Tcount[delete_Nsteps - 1];
     if (delete_Nlimit_varid > -1)
       delete_Nlimit = input->variable->compute_equal(delete_Nlimit_varid);
     headroom = MAX(0, delete_Nlimit - ndeletions);
@@ -925,13 +931,11 @@ void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec)
   std::random_device rnd;
   std::minstd_rand park_rng(rnd());
   int *molrange;
-  memory->create(molrange,Nmole,"reaxff/species:molrange");
+  memory->create(molrange, Nmole, "reaxff/species:molrange");
-  for (m = 0; m < Nmole; m++)
+  for (m = 0; m < Nmole; m++) molrange[m] = m + 1;
-    molrange[m] = m + 1;
   if (delete_Nsteps > 0) {
     // shuffle index when using rate_limit, in case order is biased
-    if (comm->me == 0)
+    if (comm->me == 0) std::shuffle(&molrange[0], &molrange[Nmole], park_rng);
-      std::shuffle(&molrange[0],&molrange[Nmole], park_rng);
     MPI_Bcast(&molrange[0], Nmole, MPI_INT, 0, world);
   }
 
@@ -1060,11 +1064,9 @@ void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec)
     }
   }
 
-
   // push back delete_Tcount on every step
   if (delete_Nsteps > 0) {
-    for (i = delete_Nsteps-1; i > 0; i--)
+    for (i = delete_Nsteps - 1; i > 0; i--) delete_Tcount[i] = delete_Tcount[i - 1];
-      delete_Tcount[i] = delete_Tcount[i-1];
     delete_Tcount[0] += this_delete_Tcount;
   }
 

src/REAXFF/fix_reaxff_species.h

@@ -51,8 +51,6 @@ class FixReaxFFSpecies : public Fix {
   int *Mol2Spec;
   double *clusterID;
   AtomCoord *x0;
-
-  double bg_cut;
   double **BOCut;
 
   std::vector<std::string> del_species;