add offset option to reset_mol_ids command

doc/src/reset_mol_ids.rst

@@ -11,33 +11,47 @@ Syntax
 
 .. parsed-literal::
 
-   reset_mol_ids group-ID
+   reset_mol_ids group-ID keyword value ...
 
 * group-ID = ID of group of atoms whose molecule IDs will be reset
+* zero or more keyword/value pairs may be appended
+* keyword = *offset*
 
+  .. parsed-literal::
+
+       *offset* value = *Noffset*
+       
 Examples
 """"""""
 
 .. code-block:: LAMMPS
 
    reset_mol_ids all
-   reset_mol_ids solvent
+   reset_mol_ids solvent offset 1000
 
 Description
 """""""""""
 
 Reset molecule IDs for a group of atoms.  This will create a new set
 of molecule IDs that are numbered contiguously from 1 to N, if there
-are N different molecule IDs used by the group.  Only molecule IDs for
-atoms in the specified group are reset.
+are N different molecules in the group.  Only molecule IDs for atoms
+in the specified group are reset.
 
-This can be useful to invoke after performing a reactive molecular
-dynamics run with :doc:`fix bond/react <fix_bond_react>`, :doc:`fix
-bond/create <fix_bond_create>`, or :doc:`fix bond/break
-<fix_bond_break>`. It can also be useful after molecules have been
-deleted with :doc:`delete_atoms <delete_atoms>` or after a simulation
-which has lost molecules, e.g. via the :doc:`fix evaporate
-<fix_evaporate>` command.
+For purposes of this operation, molecules are identified by the bond
+topology of the system, not by the current molecule IDs.  A molecule
+is a set of atoms, each is which is bonded to one or more atoms in the
+set.  If an atom is not bonded to any other atom, it becomes its own
+1-atom molecule.  Once new molecules are identified, this command will
+overwrite the current molecule ID for each atom with a new ID.
+
+This can be a useful operation to perform after running reactive
+molecular dynamics run with :doc:`fix bond/react <fix_bond_react>`,
+:doc:`fix bond/create <fix_bond_create>`, or :doc:`fix bond/break
+<fix_bond_break>`, all of which can change molecule topologies. It can
+also be useful after molecules have been deleted with the
+:doc:`delete_atoms <delete_atoms>` command or after a simulation which
+has lost molecules, e.g. via the :doc:`fix evaporate <fix_evaporate>`
+command.
 
 Restrictions
 """"""""""""

src/reset_mol_ids.cpp

@@ -46,12 +46,26 @@ void ResetMolIDs::command(int narg, char **arg)
   if (atom->molecular != 1)
     error->all(FLERR,"Can only use reset_mol_ids on molecular systems");
 
-  if (narg != 1) error->all(FLERR,"Illegal reset_mol_ids command");
+  // process args
+  
+  if (narg < 1) error->all(FLERR,"Illegal reset_mol_ids command");
   int igroup = group->find(arg[0]);
   if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID");
+  int groupbit = group->bitmask[igroup];
 
-  if (comm->me == 0)
-    utils::logmesg(lmp,"Resetting molecule IDs ...\n");
+  bigint offset = 0;
+  
+  int iarg = 1;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"offset") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command");
+      offset = utils::tnumeric(FLERR,arg[iarg+1],1,lmp);
+      if (offset < 0) error->all(FLERR,"Illegal reset_mol_ids command");
+      iarg += 2;
+    } else error->all(FLERR,"Illegal reset_mol_ids command");
+  }
+
+  if (comm->me == 0) utils::logmesg(lmp,"Resetting molecule IDs ...\n");
 
   // record wall time for resetting molecule IDs
 
@@ -59,13 +73,14 @@ void ResetMolIDs::command(int narg, char **arg)
   double time1 = MPI_Wtime();
 
   // create instances of compute fragment/atom and compute chunk/atom
-
+  // both use the group-ID for this command
+  
   const std::string idfrag = "reset_mol_ids_FRAGMENT_ATOM";
-  modify->add_compute(fmt::format("{} {} fragment/atom",idfrag, arg[0]));
+  modify->add_compute(fmt::format("{} {} fragment/atom",idfrag,arg[0]));
 
   const std::string idchunk = "reset_mol_ids_CHUNK_ATOM";
-  modify->add_compute(fmt::format("{} all chunk/atom molecule compress yes",
-                                  idchunk));
+  modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes",
+                                  idchunk,arg[0]));
 
   // initialize system since comm->borders() will be invoked
 
@@ -92,27 +107,33 @@ void ResetMolIDs::command(int narg, char **arg)
   double *ids = cfa->vector_atom;
 
   // copy fragmentID to molecule ID
-
+  // only for atoms in the group
+  
   tagint *molecule = atom->molecule;
+  int *mask = atom->mask;
   int nlocal = atom->nlocal;
   
   for (int i = 0; i < nlocal; i++)
-    molecule[i] = static_cast<tagint> (ids[i]);
+    if (mask[i] & groupbit)
+      molecule[i] = static_cast<tagint> (ids[i]);
 
   // invoke peratom method of compute chunk/atom
-  // compresses new molecule IDs to be contiguous 1 to Nmol
-
+  // compress new molecule IDs to be contiguous 1 to Nmol
+  // NOTE: use of compute chunk/atom limits # of molecules to 32-bit int
+  
   icompute = modify->find_compute(idchunk);
   ComputeChunkAtom *cca = (ComputeChunkAtom *) modify->compute[icompute];
   cca->compute_peratom();
   ids = cca->vector_atom;
   int nchunk = cca->nchunk;
   
-  // copy chunkID to molecule ID
-
+  // copy chunkID to molecule ID + offset
+  // only for atoms in the group
+  
   for (int i = 0; i < nlocal; i++)
-    molecule[i] = static_cast<tagint> (ids[i]);
-
+    if (mask[i] & groupbit)
+      molecule[i] = static_cast<tagint> (ids[i]) + offset;
+  
   // clean up
 
   modify->delete_compute(idchunk);
@@ -122,7 +143,9 @@ void ResetMolIDs::command(int narg, char **arg)
 
   MPI_Barrier(world);
   
-  if (comm->me == 0)
+  if (comm->me == 0) {
+    utils::logmesg(lmp,fmt::format("  number of new molecule IDs = {}\n",nchunk));
     utils::logmesg(lmp,fmt::format("  reset_mol_ids CPU = {:.3f} seconds\n",
                                    MPI_Wtime()-time1));
+  }
 }