Merge branch 'reset-molecules' of github.com:lammps/lammps into reset-molecules

.. 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.

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.

are N different molecules in the group.  Only molecule IDs for atoms

in the specified group are reset.

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.

The *offset* keyword can be used to change the range of new molecule

IDs assigned.  If *Noffset* = 0 (the default) and the specified group

is *all*, then new molecule IDs will be from 1 to N.  If *Noffset* = 0

and the group is not all, then new molecule IDs will be from M+1 to

M+N, where M is the largest molecule ID for any atom not in the group.

This insures new molecule IDs do not collide with existing molecule

IDs that are not changed by this command.

If *Noffset* is set to a value greater than 0, then new molecule IDs

will be from Noffset+1 to Noffset+N.  If the group is not all, it is

up to you to insure the choice of *Noffset* does not produce

collisions between existing and new molecule IDs.

.. note::

   The same as explained for the :doc:`compute fragment/atom

   <compute_cluster_atom>` command, molecules are identified using the

   current bond topology within each fragment.  This will not account

   for bonds broken by the :doc:`bond_style quartic <bond_quartic>`

   command because it does not perform a full update of the bond

   topology data structures within LAMMPS.

Restrictions

""""""""""""

:doc:`fix evaporate <fix_evaporate>`,

:doc:`delete_atoms <delete_atoms>`

**Default:** none

**Default:**

The default keyword setting is offset = 0.

  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];

  tagint 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");

  if (comm->me == 0) utils::logmesg(lmp,"Resetting molecule IDs ...\n");

  // record wall time for resetting molecule IDs

  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]));

  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

  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++)

    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 a 32-bit int

  icompute = modify->find_compute(idchunk);

  ComputeChunkAtom *cca = (ComputeChunkAtom *) modify->compute[icompute];

  ids = cca->vector_atom;

  int nchunk = cca->nchunk;

  // copy chunkID to molecule ID

  // if offset = 0 and group != all,

  // reset offset to largest molID of non-group atoms

  if (offset == 0 && groupbit != 1) {

    tagint mymol = 0;

    for (int i = 0; i < nlocal; i++)

    molecule[i] = static_cast<tagint> (ids[i]);

      if (!(mask[i] & groupbit))

	mymol = MAX(mymol,molecule[i]);

    MPI_Allreduce(&mymol,&offset,1,MPI_LMP_TAGINT,MPI_MAX,world);

  }

  // copy chunkID to molecule ID + offset

  // only for atoms in the group

  for (int i = 0; i < nlocal; i++)

    if (mask[i] & groupbit)

      molecule[i] = static_cast<tagint> (ids[i]) + offset;

  // clean up

  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));

  }

}