From 6f66e6c45445bfbd459ac2beeeb7a8ea1fa1c083 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 18 Aug 2017 19:07:32 -0400 Subject: [PATCH 1/7] add new compute fragment/atom --- src/compute_cluster_atom.h | 2 +- src/compute_fragment_atom.cpp | 210 ++++++++++++++++++++++++++++++++++ src/compute_fragment_atom.h | 67 +++++++++++ 3 files changed, 278 insertions(+), 1 deletion(-) create mode 100644 src/compute_fragment_atom.cpp create mode 100644 src/compute_fragment_atom.h diff --git a/src/compute_cluster_atom.h b/src/compute_cluster_atom.h index aae4747781..253396b40e 100644 --- a/src/compute_cluster_atom.h +++ b/src/compute_cluster_atom.h @@ -59,7 +59,7 @@ E: Cannot use compute cluster/atom unless atoms have IDs Atom IDs are used to identify clusters. -E: Compute cluster/atom requires a pair style be defined +E: Compute cluster/atom requires a pair style to be defined This is so that the pair style defines a cutoff distance which is used to find clusters. diff --git a/src/compute_fragment_atom.cpp b/src/compute_fragment_atom.cpp new file mode 100644 index 0000000000..02a0b8e7dc --- /dev/null +++ b/src/compute_fragment_atom.cpp @@ -0,0 +1,210 @@ +/* ---------------------------------------------------------------------- + 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 +#include "compute_fragment_atom.h" +#include "atom.h" +#include "atom_vec.h" +#include "update.h" +#include "modify.h" +#include "neighbor.h" +#include "force.h" +#include "pair.h" +#include "comm.h" +#include "memory.h" +#include "error.h" + +#include "group.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeFragmentAtom::ComputeFragmentAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), + fragmentID(NULL) +{ + if (narg != 3) error->all(FLERR,"Illegal compute fragment/atom command"); + + if (atom->avec->bonds_allow == 0) + error->all(FLERR,"Compute fragment/atom used when bonds are not allowed"); + + peratom_flag = 1; + size_peratom_cols = 0; + comm_forward = 1; + + nmax = 0; +} + +/* ---------------------------------------------------------------------- */ + +ComputeFragmentAtom::~ComputeFragmentAtom() +{ + memory->destroy(fragmentID); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeFragmentAtom::init() +{ + 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"); + + int count = 0; + for (int i = 0; i < modify->ncompute; i++) + if (strcmp(modify->compute[i]->style,"fragment/atom") == 0) count++; + if (count > 1 && comm->me == 0) + error->warning(FLERR,"More than one compute fragment/atom"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeFragmentAtom::compute_peratom() +{ + int i,j,k; + + invoked_peratom = update->ntimestep; + + // grow fragmentID array if necessary + + if (atom->nmax > nmax) { + memory->destroy(fragmentID); + nmax = atom->nmax; + memory->create(fragmentID,nmax,"fragment/atom:fragmentID"); + vector_atom = fragmentID; + } + + int nbondlist = neighbor->nbondlist; + int **bondlist = neighbor->bondlist; + + // if group is dynamic, insure ghost atom masks are current + + if (group->dynamic[igroup]) { + commflag = 0; + comm->forward_comm_compute(this); + } + + // every bond starts in its own fragment, + // with fragmentID = MIN(b1atomID,b2atomID) + // only bonds wholly contained in the group are considered + + tagint *tag = atom->tag; + int *mask = atom->mask; + + for (i = 0; i < nbondlist; i++) { + const int b1 = bondlist[i][0]; + const int b2 = bondlist[i][1]; + + if ((mask[b1] & groupbit) && (mask[b2] & groupbit)) + fragmentID[b1] = fragmentID[b2] = MIN(tag[b1],tag[b2]); + else fragmentID[b1] = fragmentID[b2] = 0; + } + + // loop until no more changes on any proc: + // 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 + + commflag = 1; + int nlocal = atom->nlocal; + int *num_bond = atom->num_bond; + tagint **bond_atom = atom->bond_atom; + + 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++) { + k = bond_atom[i][j]; + if (!(mask[k] & groupbit)) continue; + if (fragmentID[i] == fragmentID[k]) continue; + + fragmentID[i] = fragmentID[k] = MIN(fragmentID[i],fragmentID[k]); + 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 ComputeFragmentAtom::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++] = fragmentID[j]; + } + } else { + int *mask = atom->mask; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = ubuf(mask[j]).d; + } + } + + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeFragmentAtom::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++) fragmentID[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 ComputeFragmentAtom::memory_usage() +{ + double bytes = nmax * sizeof(double); + return bytes; +} diff --git a/src/compute_fragment_atom.h b/src/compute_fragment_atom.h new file mode 100644 index 0000000000..a56239dbda --- /dev/null +++ b/src/compute_fragment_atom.h @@ -0,0 +1,67 @@ +/* -*- 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(fragment/atom,ComputeFragmentAtom) + +#else + +#ifndef LMP_COMPUTE_FRAGMENT_ATOM_H +#define LMP_COMPUTE_FRAGMENT_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeFragmentAtom : public Compute { + public: + ComputeFragmentAtom(class LAMMPS *, int, char **); + ~ComputeFragmentAtom(); + void init(); + 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 *fragmentID; +}; + +} + +#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 fragment/atom unless atoms have IDs + +Atom IDs are used to identify fragments. + +E: Compute fragment/atom requires a bond style to be defined + +This is so that a bond list is generated which is used to find fragments. + +W: More than one compute fragment/atom + +It is not efficient to use compute fragment/atom more than once. + +*/ From c5ce3ffe60e40f186cb8b63b4233ccfed7d229dd Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 20 Aug 2017 09:18:04 -0400 Subject: [PATCH 2/7] use list of bonds per atom instead of bondlist, as that will work with shake as well --- src/compute_fragment_atom.cpp | 29 ++++++++++------------------- 1 file changed, 10 insertions(+), 19 deletions(-) diff --git a/src/compute_fragment_atom.cpp b/src/compute_fragment_atom.cpp index 02a0b8e7dc..5024ada4f4 100644 --- a/src/compute_fragment_atom.cpp +++ b/src/compute_fragment_atom.cpp @@ -90,9 +90,6 @@ void ComputeFragmentAtom::compute_peratom() vector_atom = fragmentID; } - int nbondlist = neighbor->nbondlist; - int **bondlist = neighbor->bondlist; - // if group is dynamic, insure ghost atom masks are current if (group->dynamic[igroup]) { @@ -100,21 +97,17 @@ void ComputeFragmentAtom::compute_peratom() comm->forward_comm_compute(this); } - // every bond starts in its own fragment, - // with fragmentID = MIN(b1atomID,b2atomID) - // only bonds wholly contained in the group are considered + // each atom starts in its own fragment, + int nlocal = atom->nlocal; tagint *tag = atom->tag; int *mask = atom->mask; + int *num_bond = atom->num_bond; + tagint **bond_atom = atom->bond_atom; - for (i = 0; i < nbondlist; i++) { - const int b1 = bondlist[i][0]; - const int b2 = bondlist[i][1]; - - if ((mask[b1] & groupbit) && (mask[b2] & groupbit)) - fragmentID[b1] = fragmentID[b2] = MIN(tag[b1],tag[b2]); - else fragmentID[b1] = fragmentID[b2] = 0; - } + for (i = 0; i < nlocal + atom->nghost; i++) + if (mask[i] & groupbit) fragmentID[i] = tag[i]; + else fragmentID[i] = 0; // loop until no more changes on any proc: // acquire fragmentIDs of ghost atoms @@ -124,9 +117,6 @@ void ComputeFragmentAtom::compute_peratom() // then check if any proc made changes commflag = 1; - int nlocal = atom->nlocal; - int *num_bond = atom->num_bond; - tagint **bond_atom = atom->bond_atom; int change,done,anychange; @@ -138,9 +128,10 @@ void ComputeFragmentAtom::compute_peratom() done = 1; for (i = 0; i < nlocal; i++) { if (!(mask[i] & groupbit)) continue; - + for (j = 0; j < num_bond[i]; j++) { - k = bond_atom[i][j]; + k = atom->map(bond_atom[i][j]); + if (k < 0) continue; if (!(mask[k] & groupbit)) continue; if (fragmentID[i] == fragmentID[k]) continue; From 35fd82b602828463dfa151325086b2e2d67309b2 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 20 Aug 2017 09:19:04 -0400 Subject: [PATCH 3/7] trim unneeded includes --- src/compute_fragment_atom.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/compute_fragment_atom.cpp b/src/compute_fragment_atom.cpp index 5024ada4f4..7417c26831 100644 --- a/src/compute_fragment_atom.cpp +++ b/src/compute_fragment_atom.cpp @@ -21,9 +21,7 @@ #include "atom_vec.h" #include "update.h" #include "modify.h" -#include "neighbor.h" #include "force.h" -#include "pair.h" #include "comm.h" #include "memory.h" #include "error.h" From 4dc1195cd8511cb08b203eef1da0ac3c5a10fd49 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 20 Aug 2017 09:41:49 -0400 Subject: [PATCH 4/7] add docs for compute fragment/atom --- doc/src/Section_commands.txt | 1 + doc/src/compute_cluster_atom.txt | 38 +++++++++++++++++++++----------- 2 files changed, 26 insertions(+), 13 deletions(-) diff --git a/doc/src/Section_commands.txt b/doc/src/Section_commands.txt index 571c6c4920..0211bbe055 100644 --- a/doc/src/Section_commands.txt +++ b/doc/src/Section_commands.txt @@ -775,6 +775,7 @@ KOKKOS, o = USER-OMP, t = OPT. "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, diff --git a/doc/src/compute_cluster_atom.txt b/doc/src/compute_cluster_atom.txt index 147d06c2a8..c01c2010c5 100644 --- a/doc/src/compute_cluster_atom.txt +++ b/doc/src/compute_cluster_atom.txt @@ -7,37 +7,49 @@ :line compute cluster/atom command :h3 +compute fragment/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 :pre ID, group-ID are documented in "compute"_compute.html command -cluster/atom = style name of this compute command -cutoff = distance within which to label atoms as part of same cluster (distance units) :ul +{cluster/atom} or {fragment/atom} = style name of this compute command +cutoff = distance within which to label atoms as part of same cluster (distance units, cluster/atom only) :ul [Examples:] -compute 1 all cluster/atom 1.0 :pre +compute 1 all cluster/atom 3.5 +compute 1 all fragment/atom :pre [Description:] -Define a computation that assigns each atom a cluster ID. +Define a computation that assigns each atom a cluster ID or fragment 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 define 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. Only atoms in the compute group are clustered and assigned cluster -IDs. Atoms not in the compute group are assigned a cluster ID = 0. +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 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. -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 compute {cluster/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} style. NOTE: If you have a bonded system, then the settings of "special_bonds"_special_bonds.html command can remove pairwise From c895df73d624570703d2c56256c87abe92421ccf Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 20 Aug 2017 09:47:51 -0400 Subject: [PATCH 5/7] skip over disabled bonds --- src/compute_fragment_atom.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/compute_fragment_atom.cpp b/src/compute_fragment_atom.cpp index 7417c26831..2dfb20a570 100644 --- a/src/compute_fragment_atom.cpp +++ b/src/compute_fragment_atom.cpp @@ -101,6 +101,7 @@ void ComputeFragmentAtom::compute_peratom() 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; for (i = 0; i < nlocal + atom->nghost; i++) @@ -128,6 +129,7 @@ void ComputeFragmentAtom::compute_peratom() 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; From 5a0c3aea8ac1413ae026e011171a70b37378ec28 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 21 Aug 2017 13:12:43 -0400 Subject: [PATCH 6/7] add a compute aggregate/atom, that combines the rules for compute cluster/atom and fragment/atom --- doc/src/compute_cluster_atom.txt | 33 ++-- src/compute_aggregate_atom.cpp | 270 +++++++++++++++++++++++++++++++ src/compute_aggregate_atom.h | 70 ++++++++ 3 files changed, 363 insertions(+), 10 deletions(-) create mode 100644 src/compute_aggregate_atom.cpp create mode 100644 src/compute_aggregate_atom.h diff --git a/doc/src/compute_cluster_atom.txt b/doc/src/compute_cluster_atom.txt index c01c2010c5..0aa38ae590 100644 --- a/doc/src/compute_cluster_atom.txt +++ b/doc/src/compute_cluster_atom.txt @@ -8,29 +8,35 @@ compute cluster/atom command :h3 compute fragment/atom command :h3 +compute aggregate/atom command :h3 [Syntax:] compute ID group-ID cluster/atom cutoff -compute ID group-ID fragment/atom :pre +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} or {fragment/atom} = style name of this compute command -cutoff = distance within which to label atoms as part of same cluster (distance units, cluster/atom only) :ul +{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 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 or fragment 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. A fragment is similarly define as a set of atoms, each of +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, @@ -38,6 +44,12 @@ 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. For fragments, only bonds where [both] atoms of the bond are included @@ -45,11 +57,12 @@ 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 compute {cluster/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} style. +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 diff --git a/src/compute_aggregate_atom.cpp b/src/compute_aggregate_atom.cpp new file mode 100644 index 0000000000..1155ac437a --- /dev/null +++ b/src/compute_aggregate_atom.cpp @@ -0,0 +1,270 @@ +/* ---------------------------------------------------------------------- + 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 +#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; +} diff --git a/src/compute_aggregate_atom.h b/src/compute_aggregate_atom.h new file mode 100644 index 0000000000..8170aabc7f --- /dev/null +++ b/src/compute_aggregate_atom.h @@ -0,0 +1,70 @@ +/* -*- 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. + +*/ From 24c00b1f7a202b45100eb59f94843b484bcc5c3a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 21 Aug 2017 13:12:48 -0400 Subject: [PATCH 7/7] fix typo --- src/compute_cluster_atom.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/compute_cluster_atom.cpp b/src/compute_cluster_atom.cpp index 80c002d39f..5ee6368504 100644 --- a/src/compute_cluster_atom.cpp +++ b/src/compute_cluster_atom.cpp @@ -63,7 +63,7 @@ void ComputeClusterAtom::init() 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");