diff --git a/doc/src/reset_mol_ids.rst b/doc/src/reset_mol_ids.rst index 8d2d6eeab5..0bfcccb219 100644 --- a/doc/src/reset_mol_ids.rst +++ b/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 `, :doc:`fix -bond/create `, or :doc:`fix bond/break -`. It can also be useful after molecules have been -deleted with :doc:`delete_atoms ` or after a simulation -which has lost molecules, e.g. via the :doc:`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 `, +:doc:`fix bond/create `, or :doc:`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 ` command or after a simulation which +has lost molecules, e.g. via the :doc:`fix evaporate ` +command. Restrictions """""""""""" diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp index 81ecf12a42..0b51dd08c8 100644 --- a/src/reset_mol_ids.cpp +++ b/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 (ids[i]); + if (mask[i] & groupbit) + molecule[i] = static_cast (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 (ids[i]); - + if (mask[i] & groupbit) + molecule[i] = static_cast (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)); + } }