Merge pull request #624 from akohlmey/compute-fragment

"erotate/sphere"_compute_erotate_sphere.html,

"erotate/sphere/atom"_compute_erotate_sphere_atom.html,

"event/displace"_compute_event_displace.html,

"fragment/atom"_compute_cluster_atom.html,

"global/atom"_compute_global_atom.html,

"group/group"_compute_group_group.html,

"gyration"_compute_gyration.html,

:line

compute cluster/atom command :h3

compute fragment/atom command :h3

compute aggregate/atom command :h3

[Syntax:]

compute ID group-ID cluster/atom cutoff :pre

compute ID group-ID cluster/atom cutoff

compute ID group-ID fragment/atom

compute ID group-ID aggregate/atom cutoff :pre

ID, group-ID are documented in "compute"_compute.html command

cluster/atom = style name of this 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) :ul

[Examples:]

compute 1 all cluster/atom 1.0 :pre

compute 1 all cluster/atom 3.5

compute 1 all fragment/atom :pre

compute 1 all aggregate/atom 3.5 :pre

[Description:]

Define a computation that assigns each atom a cluster ID.

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

or aggregate ID.

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

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

cluster.  The ID of every atom in the cluster will be the smallest

atom ID of any atom in the cluster.

cluster.

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

which has an explicit bond (i.e. defined via a "data file"_read_data.html,

the "create_bonds"_create_bonds.html command, or through fixes like

"fix bond/create"_fix_bond_create.html, "fix bond/swap"_fix_bond_swap.html,

or "fix bond/break"_fix_bond_break.html).  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

distance from one or more atoms within a fragment that is part of

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.

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

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

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

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.

NOTE: If you have a bonded system, then the settings of

"special_bonds"_special_bonds.html command can remove pairwise

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

   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator

   http://lammps.sandia.gov, Sandia National Laboratories

   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract

   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains

   certain rights in this software.  This software is distributed under

   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.

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

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

   Contributing author: Axel Kohlmeyer (Temple U)

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

#include <string.h>

#include "compute_aggregate_atom.h"

#include "atom.h"

#include "atom_vec.h"

#include "update.h"

#include "modify.h"

#include "neighbor.h"

#include "neigh_list.h"

#include "neigh_request.h"

#include "pair.h"

#include "force.h"

#include "comm.h"

#include "memory.h"

#include "error.h"

#include "group.h"

using namespace LAMMPS_NS;

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

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

  Compute(lmp, narg, arg),

  aggregateID(NULL)

{

  if (narg != 4) error->all(FLERR,"Illegal compute aggregate/atom command");

  double cutoff = force->numeric(FLERR,arg[3]);

  cutsq = cutoff*cutoff;

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

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

  peratom_flag = 1;

  size_peratom_cols = 0;

  comm_forward = 1;

  nmax = 0;

}

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

ComputeAggregateAtom::~ComputeAggregateAtom()

{

  memory->destroy(aggregateID);

}

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

void ComputeAggregateAtom::init()

{

  if (atom->tag_enable == 0)

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

  if (force->bond == NULL)

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

  if (force->pair == NULL)

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

  if (sqrt(cutsq) > force->pair->cutforce)

    error->all(FLERR,

               "Compute cluster/atom cutoff is longer than pairwise cutoff");

  // need an occasional full neighbor list

  // full required so that pair of atoms on 2 procs both set their clusterID

  int irequest = neighbor->request(this,instance_me);

  neighbor->requests[irequest]->pair = 0;

  neighbor->requests[irequest]->compute = 1;

  neighbor->requests[irequest]->half = 0;

  neighbor->requests[irequest]->full = 1;

  neighbor->requests[irequest]->occasional = 1;

  int count = 0;

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

    if (strcmp(modify->compute[i]->style,"aggregate/atom") == 0) count++;

  if (count > 1 && comm->me == 0)

    error->warning(FLERR,"More than one compute aggregate/atom");

}

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

void ComputeAggregateAtom::init_list(int id, NeighList *ptr)

{

  list = ptr;

}

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

void ComputeAggregateAtom::compute_peratom()

{

  int i,j,k;

  invoked_peratom = update->ntimestep;

  // grow aggregateID array if necessary

  if (atom->nmax > nmax) {

    memory->destroy(aggregateID);

    nmax = atom->nmax;

    memory->create(aggregateID,nmax,"aggregate/atom:aggregateID");

    vector_atom = aggregateID;

  }

  // invoke full neighbor list (will copy or build if necessary)

  neighbor->build_one(list);

  // if group is dynamic, insure ghost atom masks are current

  if (group->dynamic[igroup]) {

    commflag = 0;

    comm->forward_comm_compute(this);

  }

  // each atom starts in its own aggregate,

  int nlocal = atom->nlocal;

  int inum = list->inum;

  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;

  int *ilist = list->ilist;

  int *numneigh = list->numneigh;

  int **firstneigh = list->firstneigh;

  double **x = atom->x;

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

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

    else aggregateID[i] = 0;

  // loop until no more changes on any proc:

  // acquire aggregateIDs of ghost atoms

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

  // if both atoms are in aggregate, assign lowest aggregateID to both

  // then loop over my atoms, checking distance to neighbors

  // if both atoms are in cluster, assign lowest clusterID to both

  // iterate until no changes in my atoms

  // then check if any proc made changes

  commflag = 1;

  int change,done,anychange;

  while (1) {

    comm->forward_comm_compute(this);

    change = 0;

    while (1) {

      done = 1;

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

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

        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 (aggregateID[i] == aggregateID[k]) continue;

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

          done = 0;

        }

      }

      for (int ii = 0; ii < inum; ii++) {

        i = ilist[ii];

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

        const double xtmp = x[i][0];

        const double ytmp = x[i][1];

        const double ztmp = x[i][2];

        int *jlist = firstneigh[i];

        const int jnum = numneigh[i];

        for (int jj = 0; jj < jnum; jj++) {

          j = jlist[jj];

          j &= NEIGHMASK;

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

          if (aggregateID[i] == aggregateID[j]) continue;

          const double delx = xtmp - x[j][0];

          const double dely = ytmp - x[j][1];

          const double delz = ztmp - x[j][2];

          const double rsq = delx*delx + dely*dely + delz*delz;

          if (rsq < cutsq) {

            aggregateID[i] = aggregateID[j]

              = MIN(aggregateID[i],aggregateID[j]);

            done = 0;

          }

        }

      }

      if (!done) change = 1;

      if (done) break;

    }

    // stop if all procs are done

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

    if (!anychange) break;

  }

}

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

int ComputeAggregateAtom::pack_forward_comm(int n, int *list, double *buf,

                                          int pbc_flag, int *pbc)

{

  int i,j,m;

  m = 0;

  if (commflag) {

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

      j = list[i];

      buf[m++] = aggregateID[j];

    }

  } else {

    int *mask = atom->mask;

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

      j = list[i];

      buf[m++] = ubuf(mask[j]).d;

    }

  }

  return m;

}

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

void ComputeAggregateAtom::unpack_forward_comm(int n, int first, double *buf)

{

  int i,m,last;

  m = 0;

  last = first + n;

  if (commflag)

    for (i = first; i < last; i++) aggregateID[i] = buf[m++];

  else {

    int *mask = atom->mask;

    for (i = first; i < last; i++) mask[i] = (int) ubuf(buf[m++]).i;

  }

}

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

   memory usage of local atom-based array

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

double ComputeAggregateAtom::memory_usage()

{

  double bytes = nmax * sizeof(double);

  return bytes;

}

/* -*- c++ -*- ----------------------------------------------------------

   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator

   http://lammps.sandia.gov, Sandia National Laboratories

   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract

   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains

   certain rights in this software.  This software is distributed under

   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.

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

#ifdef COMPUTE_CLASS

ComputeStyle(aggregate/atom,ComputeAggregateAtom)

#else

#ifndef LMP_COMPUTE_AGGREGATE_ATOM_H

#define LMP_COMPUTE_AGGREGATE_ATOM_H

#include "compute.h"

namespace LAMMPS_NS {

class ComputeAggregateAtom : public Compute {

 public:

  ComputeAggregateAtom(class LAMMPS *, int, char **);

  ~ComputeAggregateAtom();

  void init();

  void init_list(int, class NeighList *);

  void compute_peratom();

  int pack_forward_comm(int, int *, double *, int, int *);

  void unpack_forward_comm(int, int, double *);

  double memory_usage();

 private:

  int nmax,commflag;

  double cutsq;

  class NeighList *list;

  double *aggregateID;

};

}

#endif

#endif

/* ERROR/WARNING messages:

E: Illegal ... command

Self-explanatory.  Check the input script syntax and compare to the

documentation for the command.  You can use -echo screen as a

command-line option when running LAMMPS to see the offending line.

E: Cannot use compute aggregate/atom unless atoms have IDs

Atom IDs are used to identify aggregates.

E: Compute aggregate/atom requires a bond style to be defined

This is so that a bond list is generated which is used to find aggregates.

W: More than one compute aggregate/atom

It is not efficient to use compute aggregate/atom  more than once.

*/

  if (atom->tag_enable == 0)

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

  if (force->pair == NULL)

    error->all(FLERR,"Compute cluster/atom requires a pair style be defined");

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

  if (sqrt(cutsq) > force->pair->cutforce)

    error->all(FLERR,

               "Compute cluster/atom cutoff is longer than pairwise cutoff");