diff --git a/src/compute_contact_atom.cpp b/src/compute_contact_atom.cpp new file mode 100644 index 0000000000..ef3d4ff581 --- /dev/null +++ b/src/compute_contact_atom.cpp @@ -0,0 +1,197 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "math.h" +#include "string.h" +#include "stdlib.h" +#include "compute_contact_atom.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "force.h" +#include "pair.h" +#include "comm.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeContactAtom::ComputeContactAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 3) error->all(FLERR,"Illegal compute contact/atom command"); + + peratom_flag = 1; + size_peratom_cols = 0; + comm_reverse = 1; + + nmax = 0; + contact = NULL; + + // error checks + + if (!atom->sphere_flag) + error->all(FLERR,"Compute contact/atom requires atom style sphere"); +} + +/* ---------------------------------------------------------------------- */ + +ComputeContactAtom::~ComputeContactAtom() +{ + memory->destroy(contact); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeContactAtom::init() +{ + if (force->pair == NULL) + error->all(FLERR,"Compute contact/atom requires a pair style be defined"); + + int count = 0; + for (int i = 0; i < modify->ncompute; i++) + if (strcmp(modify->compute[i]->style,"contact/atom") == 0) count++; + if (count > 1 && comm->me == 0) + error->warning(FLERR,"More than one compute contact/atom"); + + // need an occasional neighbor list + + int irequest = neighbor->request((void *) this); + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->gran = 1; + neighbor->requests[irequest]->pair = 0; + neighbor->requests[irequest]->compute = 1; + neighbor->requests[irequest]->occasional = 1; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeContactAtom::init_list(int id, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeContactAtom::compute_peratom() +{ + int i,j,ii,jj,inum,jnum; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,radsumsq; + int *ilist,*jlist,*numneigh,**firstneigh; + + invoked_peratom = update->ntimestep; + + // grow contact array if necessary + + if (atom->nmax > nmax) { + memory->destroy(contact); + nmax = atom->nmax; + memory->create(contact,nmax,"contact/atom:contact"); + vector_atom = contact; + } + + // invoke neighbor list (will copy or build if necessary) + + neighbor->build_one(list->index); + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // compute number of contacts for each atom in group + // contact if distance <= sum of radii + // tally for both I and J + + double **x = atom->x; + double *radius = atom->radius; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + + for (i = 0; i < nall; i++) contact[i] = 0.0; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (mask[i] & groupbit) { + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + radsumsq = radsum*radsum; + if (rsq <= radsumsq) { + contact[i] += 1.0; + contact[j] += 1.0; + } + } + } + } + + // communicate ghost atom counts between neighbor procs if necessary + + if (force->newton_pair) comm->reverse_comm_compute(this); +} + +/* ---------------------------------------------------------------------- */ + +int ComputeContactAtom::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++] = contact[i]; + return 1; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeContactAtom::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + contact[j] += buf[m++]; + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double ComputeContactAtom::memory_usage() +{ + double bytes = nmax * sizeof(double); + return bytes; +} diff --git a/src/compute_contact_atom.h b/src/compute_contact_atom.h new file mode 100644 index 0000000000..c602131544 --- /dev/null +++ b/src/compute_contact_atom.h @@ -0,0 +1,65 @@ +/* ---------------------------------------------------------------------- + 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(contact/atom,ComputeContactAtom) + +#else + +#ifndef LMP_COMPUTE_CONTACT_ATOM_H +#define LMP_COMPUTE_CONTACT_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeContactAtom : public Compute { + public: + ComputeContactAtom(class LAMMPS *, int, char **); + ~ComputeContactAtom(); + void init(); + void init_list(int, class NeighList *); + void compute_peratom(); + int pack_reverse_comm(int, int, double *); + void unpack_reverse_comm(int, int *, double *); + double memory_usage(); + + private: + int nmax; + class NeighList *list; + double *contact; +}; + +} + +#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: Compute contact/atom requires a pair style be defined + +Self-explantory. + +W: More than one compute contact/atom + +It is not efficient to use compute contact/atom more than once. + +*/ diff --git a/src/neigh_request.h b/src/neigh_request.h index f5947cf729..af3e39839d 100644 --- a/src/neigh_request.h +++ b/src/neigh_request.h @@ -26,18 +26,18 @@ class NeighRequest : protected Pointers { // which class is requesting the list, one flag is 1, others are 0 - int pair; + int pair; // set by default int fix; int compute; int command; - // kind of list requested + // kind of list requested, one flag is 1, others are 0 // set by requesting class - int half; // 1 if half neigh list + int half; // 1 if half neigh list (set by default) int full; // 1 if full neigh list - int gran; // 1 if granular list + int gran; // 1 if granular list int granhistory; // 1 if granular history list int respainner; // 1 if a rRESPA inner list @@ -59,7 +59,7 @@ class NeighRequest : protected Pointers { int newton; // 0 if user of list wants no encoding of special bond flags and all neighs - // 1 if user of list wants special bond flags encoded + // 1 if user of list wants special bond flags encoded, set by default int special; @@ -89,7 +89,7 @@ class NeighRequest : protected Pointers { int *iskip; // iskip[i] if atoms of type I are not in list int **ijskip; // ijskip[i][j] if pairs of type I,J are not in list - int otherlist; // other list to copy or skip from + int otherlist; // index of other list to copy or skip from // methods diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 08a4b317b8..7dbcbb701f 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -528,6 +528,8 @@ void Neighbor::init() requests[j]->skip == 0 && requests[j]->half) break; if (requests[i]->full && requests[j]->pair && requests[j]->skip == 0 && requests[j]->full) break; + if (requests[i]->gran && requests[j]->pair && + requests[j]->skip == 0 && requests[j]->gran) break; if (requests[i]->half && requests[j]->pair && requests[j]->skip == 0 && requests[j]->respaouter) break; }