new reset_mol_ids command

   * :doc:`replicate <replicate>`

   * :doc:`rerun <rerun>`

   * :doc:`reset_ids <reset_ids>`

   * :doc:`reset_mol_ids <reset_mol_ids>`

   * :doc:`reset_timestep <reset_timestep>`

   * :doc:`restart <restart>`

   * :doc:`run <run>`

   replicate

   rerun

   reset_ids

   reset_mol_ids

   reset_timestep

   restart

   run

   compute 1 all cluster/atom 3.5

   compute 1 all fragment/atom

   compute 1 all aggregate/atom 3.5

Description

atom has no neighbors within the cutoff distance, then it is a 1-atom

cluster.

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

which has an explicit bond (i.e. defined via a :doc:`data file <read_data>`,

the :doc:`create_bonds <create_bonds>` command, or through fixes like

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

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.

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

   :doc:`special_bonds <special_bonds>` command that includes all pairs in

   the neighbor list.

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

**Output info:**

This compute calculates a per-atom vector, which can be accessed by

.. index:: reset_mol_ids

reset_mol_ids command

=================

Syntax

""""""

Syntax

""""""

.. parsed-literal::

   reset_mol_ids group-ID

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

Examples

""""""""

.. code-block:: LAMMPS

   reset_mol_ids all

   reset_mol_ids solvent

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/delete

<fix_bond_delete>`. It can also be useful after any simulation which

has lost molecules, e.g. via the :doc:`fix evaporate <fix_evaporate>`

command.

Restrictions

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

none

Related commands

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

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

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

:doc:`fix bond/delete <fix_bond_delete>`,

:doc:`fix evaporate <fix_evaporate>`

**Default:** none

#include "update.h"

#include "modify.h"

#include "force.h"

#include "group.h"

#include "comm.h"

#include "memory.h"

#include "error.h"

#include "group.h"

using namespace LAMMPS_NS;

#define BIG 1.0e20

/* ---------------------------------------------------------------------- */

ComputeFragmentAtom::ComputeFragmentAtom(LAMMPS *lmp, int narg, char **arg) :

  peratom_flag = 1;

  size_peratom_cols = 0;

  comm_forward = 1;

  comm_reverse = 1;

  nmax = 0;

  stack = NULL;

  clist = NULL;

  markflag = NULL;

}

/* ---------------------------------------------------------------------- */

ComputeFragmentAtom::~ComputeFragmentAtom()

{

  memory->destroy(stack);

  memory->destroy(clist);

  memory->destroy(markflag);

  memory->destroy(fragmentID);

}

{

  if (atom->tag_enable == 0)

    error->all(FLERR,"Cannot use compute fragment/atom unless atoms have IDs");

  if (force->bond == NULL)

    error->all(FLERR,"Compute fragment/atom requires a bond style to be defined");

  if (!atom->molecular)

    error->all(FLERR,"Compute fragment/atom requires a molecular system");

  int count = 0;

  for (int i = 0; i < modify->ncompute; i++)

void ComputeFragmentAtom::compute_peratom()

{

  int i,j,k;

  int i,j,k,m,n;

  int nstack,ncluster,done,alldone;

  double newID,cID;

  tagint *list;

  invoked_peratom = update->ntimestep;

  // grow fragmentID array if necessary

  // grow work and fragmentID vectors if necessary

  if (atom->nmax > nmax) {

    memory->destroy(stack);

    memory->destroy(clist);

    memory->destroy(markflag);

    memory->destroy(fragmentID);

    nmax = atom->nmax;

    memory->create(stack,nmax,"fragment/atom:stack");

    memory->create(clist,nmax,"fragment/atom:clist");

    memory->create(markflag,nmax,"fragment/atom:markflag");

    memory->create(fragmentID,nmax,"fragment/atom:fragmentID");

    vector_atom = fragmentID;

  }

    comm->forward_comm_compute(this);

  }

  // each atom starts in its own fragment,

  // owned + ghost atoms start with fragmentID = atomID

  // atoms not in group have fragmentID = 0

  int nlocal = atom->nlocal;

  tagint *tag = atom->tag;

  int *mask = atom->mask;

  int *num_bond = atom->num_bond;

  int **bond_type = atom->bond_type;

  tagint **bond_atom = atom->bond_atom;

  tagint **special = atom->special;

  int **nspecial = atom->nspecial;

  int nlocal = atom->nlocal;

  int nall = nlocal + atom->nghost;

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

  for (i = 0; i < nall; i++)

    if (mask[i] & groupbit) fragmentID[i] = tag[i];

    else fragmentID[i] = 0;

  // loop until no more changes on any proc:

  // loop until no ghost atom fragment ID is changed

  // acquire fragmentIDs of ghost atoms

  // loop over my atoms, and check atoms bound to it

  // if both atoms are in fragment, assign lowest fragmentID to both

  // iterate until no changes in my atoms

  // then check if any proc made changes

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

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

  // iterate until no changes to ghost atom fragmentIDs

  commflag = 1;

  int change,done,anychange;

  int iteration = 0;

  while (1) {

    iteration++;

    comm->forward_comm_compute(this);

    done = 1;

    // reverse communication when bonds are not stored on every processor

    // set markflag = 0 for all owned atoms, for new iteration

    if (force->newton_bond)

      comm->reverse_comm_compute(this);

    for (i = 0; i < nlocal; i++) markflag[i] = 0;

    // loop over all owned atoms

    // each unmarked atom starts a cluster search

    change = 0;

    while (1) {

      done = 1;

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

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

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

      if (markflag[i]) continue;

      // find one cluster of bond-connected atoms

      // ncluster = # of owned and ghost atoms in cluster

      // clist = vector of local indices of the ncluster atoms

      // stack is used to walk the bond topology

      ncluster = nstack = 0;

      stack[nstack++] = i;

      while (nstack) {

	j = stack[--nstack];

	clist[ncluster++] = j;

	markflag[j] = 1;

	n = nspecial[j][0];

	list = special[j];

        for (m = 0; m < n; m++) {

          k = atom->map(list[m]);

	  // skip bond neighbor K if not in group or already marked

        for (j = 0; j < num_bond[i]; j++) {

          if (bond_type[i][j] == 0) continue;

          k = atom->map(bond_atom[i][j]);

	  if (k < 0) continue;

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

          if (fragmentID[i] == fragmentID[k]) continue;

	  if (k < nlocal && markflag[k]) continue;

	  // owned bond neighbors are added to stack for further walking

	  // ghost bond neighbors are added directly w/out use of stack

          fragmentID[i] = fragmentID[k] = MIN(fragmentID[i],fragmentID[k]);

          done = 0;

	  if (k < nlocal) stack[nstack++] = k;

	  else clist[ncluster++] = k;

	}

      }

      if (!done) change = 1;

      if (done) break;

      // newID = minimum fragment ID in cluster list, including ghost atoms

      newID = BIG;

      for (m = 0; m < ncluster; m++) {

	cID = fragmentID[clist[m]];

	newID = MIN(newID,cID);

      }

      // set fragmentID = newID for all atoms in cluster, including ghost atoms

      // not done with iterations if change the fragmentID of a ghost atom

      for (m = 0; m < ncluster; m++) {

	j = clist[m];

	if (j >= nlocal && fragmentID[j] != newID) done = 0;

	fragmentID[j] = newID;

      }

    }

    // stop if all procs are done

    MPI_Allreduce(&change,&anychange,1,MPI_INT,MPI_MAX,world);

    if (!anychange) break;

    MPI_Allreduce(&done,&alldone,1,MPI_INT,MPI_MIN,world);

    if (alldone) break;

  }

}

  }

}

/* ---------------------------------------------------------------------- */

int ComputeFragmentAtom::pack_reverse_comm(int n, int first, double *buf)

{

  int i,m,last;

  m = 0;

  last = first + n;

  for (i = first; i < last; i++) {

    buf[m++] = fragmentID[i];

  }

  return m;

}

/* ---------------------------------------------------------------------- */

void ComputeFragmentAtom::unpack_reverse_comm(int n, int *list, double *buf)

{

  int i,j,m;

  m = 0;

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

    j = list[i];

    double x = buf[m++];

    // only overwrite local IDs with values lower than current ones

    fragmentID[j] = MIN(x,fragmentID[j]);

  }

}

/* ----------------------------------------------------------------------

   memory usage of local atom-based array

   memory usage of local atom-based arrays

------------------------------------------------------------------------- */

double ComputeFragmentAtom::memory_usage()

{

  double bytes = nmax * sizeof(double);

  bytes += 3*nmax * sizeof(int);

  return bytes;

}