minor adjustments to new reset_mol_ids command

.. parsed-literal::

   compute ID group-ID cluster/atom cutoff

   compute ID group-ID fragment/atom [singlezero]

   compute ID group-ID fragment/atom keyword value ...

   compute ID group-ID aggregate/atom cutoff

* ID, group-ID are documented in :doc:`compute <compute>` command

* *cluster/atom* or *fragment/atom* or *aggregate/atom* = style name of this compute command

* cutoff = distance within which to label atoms as part of same cluster (distance units)

* *singlezero* = keyword to trigger assigning an ID of 0 to fragments with a single atom (optional)

* zero or more keyword/value pairs may be appended to *fragment/atom*

* keyword = *single*

    .. parsed-literal::

       *single* value = *yes* or *no* to treat single atoms (no bonds) as fragments

Examples

""""""""

   compute 1 all cluster/atom 3.5

   compute 1 all fragment/atom

   compute 1 all fragment/atom single no

   compute 1 all aggregate/atom 3.5

Description

"""""""""""

Define a computation that assigns each atom a cluster, fragment,

or aggregate ID.

Define a computation that assigns each atom a cluster, fragment, or

aggregate ID.  Only atoms in the compute group are clustered and

assigned cluster IDs. Atoms not in the compute group are assigned an

ID = 0.

A cluster is defined as a set of atoms, each of which is within the

cutoff distance from one or more other atoms in the cluster.  If an

cluster.

A fragment is similarly defined as a set of atoms, each of which has a

bond to another atom in the fragment, as defined initially via the

:doc:`data file <read_data>` or :doc:`create_bonds <create_bonds>`

commands, or by fixes which dynamically create or break bonds like

:doc:`fix bond/react <fix_bond_react>`, :doc:`fix bond/create

<fix_bond_create>`, :doc:`fix bond/swap <fix_bond_swap>`, or :doc:`fix

bond/break <fix_bond_break>`.  The cluster ID or fragment ID of every

atom in the cluster will be set to the smallest atom ID of any atom in

the cluster or fragment, respectively.

The *singlezero* keyword turns on a special treatment for fragments,

where all fragments within the compute group that contain only a single

atom will have a fragment ID of 0.  This can be useful in cases where

the fragment IDs are used as input for other commands in LAMMPS that

treat such single atoms different.

bond to another atom in the fragment.  Bonds can be defined initially

via the :doc:`data file <read_data>` or :doc:`create_bonds

<create_bonds>` commands, or dynamically by fixes which create or

break bonds like :doc:`fix bond/react <fix_bond_react>`, :doc:`fix

bond/create <fix_bond_create>`, :doc:`fix bond/swap <fix_bond_swap>`,

or :doc:`fix bond/break <fix_bond_break>`.  The cluster ID or fragment

ID of every atom in the cluster will be set to the smallest atom ID of

any atom in the cluster or fragment, respectively.

For the *fragment/atom* style, the *single* keyword determines whether

single atoms (not bonded to another atom) are treated as one-atom

fragments or not, based on the *yes* or *no* setting.  If the setting

is *no* (the default), their fragment IDs are set to 0.

An aggregate is defined by combining the rules for clusters and

fragments, i.e. a set of atoms, where each of it is within the cutoff

the same cluster. This measure can be used to track molecular assemblies

like micelles.

Only atoms in the compute group are clustered and assigned cluster

IDs. Atoms not in the compute group are assigned a cluster ID = 0.

For fragments, only bonds where **both** atoms of the bond are included

in the compute group are assigned to fragments, so that only fragments

are detected where **all** atoms are in the compute group. Thus atoms

may be included in the compute group, yes still have a fragment ID of 0.

For computes *cluster/atom* and *aggregate/atom* the neighbor list needed

to compute this quantity is constructed each time the calculation is

performed (i.e. each time a snapshot of atoms is dumped).  Thus it can be

inefficient to compute/dump this quantity too frequently or to have

multiple compute/dump commands, each of a *cluster/atom* or

*aggregate/atom* style.

For computes *cluster/atom* and *aggregate/atom* a neighbor list

needed to compute cluster IDs is constructed each time the compute is

invoked.  Thus it can be inefficient to compute/dump this quantity too

frequently or to have multiple *cluster/atom* or *aggregate/atom*

style computes.

.. note::

.. note::

   For the compute fragment/atom style, each fragment is 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.

   using the current bond topology.  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.

**Output info:**

:doc:`compute coord/atom <compute_coord_atom>`

**Default:** none

**Default:**

The default for fragment/atom is single no.

* group-ID = ID of group of atoms whose molecule IDs will be reset

* zero or more keyword/value pairs may be appended

* keyword = *offset* or *singlezero*

* keyword = *compress* or *offset* or *single*

  .. parsed-literal::

       *offset* value = *Noffset* or auto

       *singlezero* = assign molecule ID 0 to atoms without bonds

       *compress* value = *yes* or *no*

       *offset* value = *Noffset* >= -1

       *single* value = *yes* or *no* to treat single atoms (no bonds) as molecules

Examples

""""""""

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

Reset molecule IDs for a group of atoms based on current bond

connectivity.  This will typically create a new set of molecule IDs

for atoms in the group.  Only molecule IDs for atoms in the specified

group are reset; molecule IDs for atoms not in the group are not

changed.

For purposes of this operation, molecules are identified by the

current bond connectivity in the system, which may or may not be

consistent with current molecule IDs.  A molecule is a set of atoms,

each of which is bonded to one or more atoms in the set.  Once new

molecules are identified and a molecule ID assigned to each one, this

command will update the current molecule ID for each atom in the group

with a (potentially) new ID.  Note that if the group is set so as to

exclude atoms within molecules, one molecule may become several.  For

example if the group excludes atoms in the midddle of a linear chain,

then each end of the chain becomes an independent molecules.

This can be a useful operation to perform after running reactive

molecular dynamics run with :doc:`fix bond/react <fix_bond_react>`,

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* = auto (the default) and the specified group

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

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 number instead of 'auto', 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.

The *singlezero* keyword turns on a special treatment for single atoms

without bonds.  Instead of assigning each atom a different molecule ID

those atoms will all be given the molecule ID 0.  This can be useful

when the molecule ID is used to group atoms where atoms with a molecule

ID of 0 will be considered as individual atoms; for example when using

:doc:`fix rigid <fix_rigid>` with the *molecule* option.

The *compress* keyword determines how new molecule IDs are assigned.

If the setting is *no* (the default), the molecule ID of every atom in

the molecule will be set to the smallest atom ID of any atom in the

molecule.  If the setting is *yes*, and there are N molecules in the

group, the new molecule IDs will be a set of N contiguous values.  See

the *offset* keyword for more details.

The *single* keyword determines whether single atoms (not bonded to

another atom) are treated as one-atom molecules or not, based on the

*yes* or *no* setting.  If the setting is *no* (the default), their

molecule IDs are set to 0.  This setting can be important if the new

molecule IDs will be used as input to other commands such as

:doc:`compute chunk/atom molecule <compute_chunk_atom>` or :doc:`fix

rigid molecule <fix_rigid>`.

The *offset* keyword is only used if the *compress* setting is *yes*.

Its default value is *Noffset* = -1.  In that case, if the specified

group is *all*, then the new compressed molecule IDs will range from 1

to N.  If the specified group is not *all* and the largest molecule ID

in the non-group atoms is M, then the new compressed molecule IDs will

range from M+1 to M+N, so as to not collide with existing molecule

IDs.  If an *Noffset* >= 0 is specified, then the new compressed

molecule IDs will range from *Noffset*+1 to *Noffset*+N.  If the group

is not *all* it is up to you to insure there are no collisions with

the molecule IDs of non-group atoms.

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

   current bond topology.  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

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

**Default:**

The default keyword setting is offset = 0.

The default keyword settings are compress = no, single = no, and

offset = -1.

  Compute(lmp, narg, arg),

  fragmentID(NULL)

{

  singleflag = 0;

  if ((narg < 3) || (narg > 4))

    error->all(FLERR,"Illegal compute fragment/atom command");

  if ((narg == 4) && (strcmp(arg[3],"singlezero") == 0))

    singleflag = 1;

  if (atom->avec->bonds_allow == 0)

    error->all(FLERR,"Compute fragment/atom used when bonds are not allowed");

  size_peratom_cols = 0;

  comm_forward = 1;

  // process optional args

  singleflag = 0;

  int iarg = 3;

  while (iarg < narg) {

    if (strcmp(arg[iarg],"single") == 0) {

      if (iarg+2 > narg) error->all(FLERR,"Illegal compute fragment/atom command");

      if (strcmp(arg[iarg+1],"yes") == 0) singleflag = 1;

      else if (strcmp(arg[iarg+1],"no") == 0) singleflag = 0;

      else error->all(FLERR,"Illegal compute fragment/atom command");

      iarg += 2;

    } else error->all(FLERR,"Illegal compute fragment/atom command");

  }

  nmax = 0;

  stack = NULL;

  clist = NULL;

  // owned + ghost atoms start with fragmentID = atomID

  // atoms not in group have fragmentID = 0

  // if singleflag is set atoms without bonds have fragmentID 0 as well.

  tagint *tag = atom->tag;

  int *mask = atom->mask;

  // acquire fragmentIDs of ghost atoms

  // loop over clusters of atoms, which include ghost atoms

  // set fragmentIDs for each cluster to min framentID in the clusters

  // if singleflag = 0 atoms without bonds are assigned fragmentID = 0

  // iterate until no changes to ghost atom fragmentIDs

  commflag = 1;

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

      // skip atom I if not in group or already marked

      // also skip and set fragment ID to zero if singleflag is set

      // and the atom is an isolated atom without bonds

      // if singleflag = 0 and atom has no bond partners, fragID = 0 and done

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

      if (markflag[i]) continue;

      if (singleflag && (nspecial[i][0] == 0)) {

        fragmentID[i] = 0;

        markflag[i] = 1;

      if (!singleflag && (nspecial[i][0] == 0)) {

	fragmentID[i] = 0.0;

	continue;

      }

#include "modify.h"

#include "compute_fragment_atom.h"

#include "compute_chunk_atom.h"

#include "compute_reduce.h"

#include "utils.h"

#include "error.h"

#include "fmt/format.h"

  if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID");

  int groupbit = group->bitmask[igroup];

  tagint offset = -1;

  int compressflag = 0;

  int singleflag = 0;

  tagint offset = -1;

  int iarg = 1;

  while (iarg < narg) {

    if (strcmp(arg[iarg],"offset") == 0) {

    if (strcmp(arg[iarg],"compress") == 0) {

      if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command");

      if (strcmp(arg[iarg+1],"yes") == 0) compressflag = 1;

      else if (strcmp(arg[iarg+1],"no") == 0) compressflag = 0;

      else error->all(FLERR,"Illegal reset_mol_ids command");

      iarg += 2;

    } else if (strcmp(arg[iarg],"single") == 0) {

      if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command");

      if (strcmp(arg[iarg+1],"yes") == 0) singleflag = 1;

      else if (strcmp(arg[iarg+1],"no") == 0) singleflag = 0;

      else error->all(FLERR,"Illegal reset_mol_ids command");

      iarg += 2;

    } else if (strcmp(arg[iarg],"offset") == 0) {

      if (iarg+2 > narg) error->all(FLERR,"Illegal reset_mol_ids command");

      if (strcmp(arg[iarg+1],"auto") == 0) {

        offset = -1;

      } else {

      offset = utils::tnumeric(FLERR,arg[iarg+1],1,lmp);

        if (offset < 0) error->all(FLERR,"Illegal reset_mol_ids command");

      }

      if (offset < -1) error->all(FLERR,"Illegal reset_mol_ids command");

      iarg += 2;

    } else if (strcmp(arg[iarg],"singlezero") == 0) {

      singleflag = 1;

      ++iarg;

    } else error->all(FLERR,"Illegal reset_mol_ids command");

  }

  const std::string idfrag = "reset_mol_ids_FRAGMENT_ATOM";

  if (singleflag)

    modify->add_compute(fmt::format("{} {} fragment/atom singlezero",idfrag,arg[0]));

    modify->add_compute(fmt::format("{} {} fragment/atom single yes",idfrag,arg[0]));

  else

    modify->add_compute(fmt::format("{} {} fragment/atom",idfrag,arg[0]));

  const std::string idmin = "reset_mol_ids_COMPUTE_MINFRAG";

  if (singleflag)

    modify->add_compute(fmt::format("{} {} reduce min c_{}",idmin,arg[0],idfrag));

    modify->add_compute(fmt::format("{} {} fragment/atom single no",idfrag,arg[0]));

  const std::string idchunk = "reset_mol_ids_CHUNK_ATOM";

  if (compressflag)

    modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes",

				    idchunk,arg[0]));

  // invoke peratom method of compute fragment/atom

  // walks bond connectivity and assigns each atom a fragment ID

  // if singleflag = 0, atoms w/out bonds will be assigned fragID = 0

  int icompute = modify->find_compute(idfrag);

  ComputeFragmentAtom *cfa = (ComputeFragmentAtom *) modify->compute[icompute];

  cfa->compute_peratom();

  double *ids = cfa->vector_atom;

  double *fragIDs = cfa->vector_atom;

  // when the lowest fragment ID is 0, in case the singlezero option is used

  // we need to remove adjust the chunk ID, so that the resuling molecule ID is correct.

  int adjid = 0;

  if (singleflag) {

    icompute = modify->find_compute(idmin);

    ComputeReduce *crd = (ComputeReduce *) modify->compute[icompute];

    if (crd->compute_scalar() == 0.0) adjid = 1;

  }

  for (int i = 0; i < atom->nlocal; i++)

    printf("FragID %d %g\n",atom->tag[i],fragIDs[i]);

  // copy fragmentID to molecule ID

  // only for atoms in the group

  // copy fragID to molecule ID for atoms in group

  tagint *molecule = atom->molecule;

  int *mask = atom->mask;

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

    if (mask[i] & groupbit)

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

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

  // invoke peratom method of compute chunk/atom

  // compress new molecule IDs to be contiguous 1 to Nmol

  // NOTE: use of compute chunk/atom limits # of molecules to a 32-bit int

  // if compressflag = 0, done

  // set nchunk = -1 since cannot easily determine # of new molecule IDs

  int nchunk = -1;

  // if compressflag = 1, invoke peratom method of compute chunk/atom

  // will compress new molecule IDs to be contiguous 1 to Nmol

  // NOTE: use of compute chunk/atom limits Nmol to a 32-bit int

  if (compressflag) {

    icompute = modify->find_compute(idchunk);

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

    cca->compute_peratom();

  ids = cca->vector_atom;

  int nchunk = cca->nchunk;

    double *chunkIDs = cca->vector_atom;

    nchunk = cca->nchunk;

    for (int i = 0; i < atom->nlocal; i++)

      printf("ChunkID %d %g\n",atom->tag[i],chunkIDs[i]);

    // if singleflag = 0, check if any single (no-bond) atoms exist

    // if so, they have fragID = 0, and compression just set their chunkID = 1

    // singleexist = 0/1 if no/yes single atoms exist with chunkID = 1

    int singleexist = 0;

    if (!singleflag) {

      int mysingle = 0;

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

	if (mask[i] & groupbit)

	  if (fragIDs[i] == 0.0) mysingle = 1;

      MPI_Allreduce(&mysingle,&singleexist,1,MPI_INT,MPI_MAX,world);

      if (singleexist) nchunk--;

    }

  // if offset = -1 (i.e. "auto") and group != all,

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

  // otherwise set to 0.

    // if offset < 0 (default), reset it

    // if group = all, offset = 0

    // else offset = largest molID of non-group atoms

  if (offset == -1) {

    if (offset < 0) {

      if (groupbit != 1) {

	tagint mymol = 0;

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

      } else offset = 0;

    }

  // copy chunkID to molecule ID + offset

  // only for atoms in the group

    // reset molecule ID for all atoms in group

    // newID = chunkID + offset

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

      if (mask[i] & groupbit) {

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

      if (singleflag && adjid) {

	tagint newid =  static_cast<tagint>(chunkIDs[i]);

	if (singleexist) {

	  if (newid == 1) newid = 0;

	  else newid += offset - 1;

	} else newid += offset;

	molecule[i] = newid;

      } else molecule[i] = newid + offset;

      }

    }

  }

  // clean up

  modify->delete_compute(idchunk);

  modify->delete_compute(idfrag);

  if (singleflag) modify->delete_compute(idmin);

  if (compressflag) modify->delete_compute(idchunk);

  // total time

  MPI_Barrier(world);

  if (comm->me == 0) {

    if (nchunk < 0)

      utils::logmesg(lmp,fmt::format("  number of new molecule IDs = unknown\n"));

    else

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