ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

add a compute aggregate/atom, that combines the rules for compute cluster/atom and fragment/atom

Browse Source
This commit is contained in:
Axel Kohlmeyer
2017-08-21 13:12:43 -04:00
parent c895df73d6
commit 5a0c3aea8a
3 changed files with 363 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

33
doc/src/compute_cluster_atom.txt
View File

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

270
src/compute_aggregate_atom.cpp Normal file
View File

@ -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 <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;
}

70
src/compute_aggregate_atom.h Normal file
View File

@ -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.
*/