Unverified Commit fd95fc98 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

add support for auto offset and singlezero option

parent d3853af4
Loading
Loading
Loading
Loading
+19 −9
Original line number Diff line number Diff line
@@ -15,11 +15,12 @@ Syntax

* group-ID = ID of group of atoms whose molecule IDs will be reset
* zero or more keyword/value pairs may be appended
* keyword = *offset*
* keyword = *offset* or *singlezero*

  .. parsed-literal::

       *offset* value = *Noffset*
       *offset* value = *Noffset* or auto
       *singlezero* = assign molecule ID 0 to atoms without bonds

Examples
""""""""
@@ -27,7 +28,9 @@ Examples
.. code-block:: LAMMPS

   reset_mol_ids all
   reset_mol_ids all offset 10 singlezero
   reset_mol_ids solvent offset 1000
   reset_mol_ids solvent offset auto

Description
"""""""""""
@@ -39,7 +42,7 @@ 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
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.
@@ -54,18 +57,25 @@ 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.
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 value greater than 0, then new molecule IDs
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.

.. note::

   The same as explained for the :doc:`compute fragment/atom
+1 −0
Original line number Diff line number Diff line
@@ -2058,6 +2058,7 @@ nodeless
nodeset
nodesets
Noehring
Noffset
noforce
Noid
nolib
+48 −11
Original line number Diff line number Diff line
@@ -25,6 +25,7 @@
#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"
@@ -53,15 +54,23 @@ void ResetMolIDs::command(int narg, char **arg)
  if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID");
  int groupbit = group->bitmask[igroup];

  tagint offset = 0;
  tagint offset = -1;
  int singleflag = 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");
      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");
      }
      iarg += 2;
    } else if (strcmp(arg[iarg],"singlezero") == 0) {
      singleflag = 1;
      ++iarg;
    } else error->all(FLERR,"Illegal reset_mol_ids command");
  }

@@ -72,12 +81,19 @@ void ResetMolIDs::command(int narg, char **arg)
  MPI_Barrier(world);
  double time1 = MPI_Wtime();

  // create instances of compute fragment/atom and compute chunk/atom
  // both use the group-ID for this command
  // create instances of compute fragment/atom, compute reduce (if needed),
  // and compute chunk/atom.  all use the group-ID for this command

  const std::string idfrag = "reset_mol_ids_FRAGMENT_ATOM";
  if (singleflag)
    modify->add_compute(fmt::format("{} {} fragment/atom singlezero",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));

  const std::string idchunk = "reset_mol_ids_CHUNK_ATOM";
  modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes",
                                  idchunk,arg[0]));
@@ -106,6 +122,16 @@ void ResetMolIDs::command(int narg, char **arg)
  cfa->compute_peratom();
  double *ids = 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;
  }

  // copy fragmentID to molecule ID
  // only for atoms in the group

@@ -127,28 +153,39 @@ void ResetMolIDs::command(int narg, char **arg)
  ids = cca->vector_atom;
  int nchunk = cca->nchunk;

  // if offset = 0 and group != all,
  // 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 && groupbit != 1) {
  if (offset == -1) {
    if (groupbit != 1) {
    tagint mymol = 0;
    for (int i = 0; i < nlocal; i++)
      if (!(mask[i] & groupbit))
        mymol = MAX(mymol,molecule[i]);
    MPI_Allreduce(&mymol,&offset,1,MPI_LMP_TAGINT,MPI_MAX,world);
    } else offset = 0;
  }

  // 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;
  for (int i = 0; i < nlocal; i++) {
    if (mask[i] & groupbit) {
      tagint newid =  static_cast<tagint>(ids[i]);
      if (singleflag && adjid) {
        if (newid == 1) newid = 0;
        else newid += offset - 1;
        molecule[i] = newid;
      } else molecule[i] = newid + offset;
    }
  }

  // clean up

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

  // total time