diff --git a/src/compute_cna_atom.cpp b/src/compute_cna_atom.cpp index 46c6726b0..024a95d9d 100644 --- a/src/compute_cna_atom.cpp +++ b/src/compute_cna_atom.cpp @@ -1,367 +1,365 @@ /* ---------------------------------------------------------------------- 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: Wan Liang (Chinese Academy of Sciences) ------------------------------------------------------------------------- */ #include <string.h> #include <stdlib.h> #include "compute_cna_atom.h" #include "atom.h" #include "update.h" #include "force.h" #include "pair.h" #include "modify.h" #include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" #include "comm.h" #include "memory.h" #include "error.h" #include <math.h> using namespace LAMMPS_NS; #define MAXNEAR 16 #define MAXCOMMON 8 enum{UNKNOWN,FCC,HCP,BCC,ICOS,OTHER}; enum{NCOMMON,NBOND,MAXBOND,MINBOND}; /* ---------------------------------------------------------------------- */ ComputeCNAAtom::ComputeCNAAtom(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) + Compute(lmp, narg, arg), + nearest(NULL), nnearest(NULL), pattern(NULL) { if (narg != 4) error->all(FLERR,"Illegal compute cna/atom command"); peratom_flag = 1; size_peratom_cols = 0; double cutoff = force->numeric(FLERR,arg[3]); if (cutoff < 0.0) error->all(FLERR,"Illegal compute cna/atom command"); cutsq = cutoff*cutoff; nmax = 0; - nearest = NULL; - nnearest = NULL; - pattern = NULL; } /* ---------------------------------------------------------------------- */ ComputeCNAAtom::~ComputeCNAAtom() { memory->destroy(nearest); memory->destroy(nnearest); memory->destroy(pattern); } /* ---------------------------------------------------------------------- */ void ComputeCNAAtom::init() { if (force->pair == NULL) error->all(FLERR,"Compute cna/atom requires a pair style be defined"); if (sqrt(cutsq) > force->pair->cutforce) error->all(FLERR,"Compute cna/atom cutoff is longer than pairwise cutoff"); // cannot use neighbor->cutneighmax b/c neighbor has not yet been init if (2.0*sqrt(cutsq) > force->pair->cutforce + neighbor->skin && comm->me == 0) error->warning(FLERR,"Compute cna/atom cutoff may be too large to find " "ghost atom neighbors"); int count = 0; for (int i = 0; i < modify->ncompute; i++) if (strcmp(modify->compute[i]->style,"cna/atom") == 0) count++; if (count > 1 && comm->me == 0) error->warning(FLERR,"More than one compute cna/atom defined"); // need an occasional full neighbor list 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; } /* ---------------------------------------------------------------------- */ void ComputeCNAAtom::init_list(int id, NeighList *ptr) { list = ptr; } /* ---------------------------------------------------------------------- */ void ComputeCNAAtom::compute_peratom() { int i,j,k,ii,jj,kk,m,n,inum,jnum,inear,jnear; int firstflag,ncommon,nbonds,maxbonds,minbonds; int nfcc,nhcp,nbcc4,nbcc6,nico,cj,ck,cl,cm; int *ilist,*jlist,*numneigh,**firstneigh; int cna[MAXNEAR][4],onenearest[MAXNEAR]; int common[MAXCOMMON],bonds[MAXCOMMON]; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; invoked_peratom = update->ntimestep; // grow arrays if necessary if (atom->nmax > nmax) { memory->destroy(nearest); memory->destroy(nnearest); memory->destroy(pattern); nmax = atom->nmax; memory->create(nearest,nmax,MAXNEAR,"cna:nearest"); memory->create(nnearest,nmax,"cna:nnearest"); memory->create(pattern,nmax,"cna:cna_pattern"); vector_atom = pattern; } // invoke full neighbor list (will copy or build if necessary) neighbor->build_one(list); inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; // find the neigbours of each atom within cutoff using full neighbor list // nearest[] = atom indices of nearest neighbors, up to MAXNEAR // do this for all atoms, not just compute group // since CNA calculation requires neighbors of neighbors double **x = atom->x; int *mask = atom->mask; int nlocal = atom->nlocal; int nerror = 0; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; jlist = firstneigh[i]; jnum = numneigh[i]; n = 0; 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; if (rsq < cutsq) { if (n < MAXNEAR) nearest[i][n++] = j; else { nerror++; break; } } } nnearest[i] = n; } // warning message int nerrorall; MPI_Allreduce(&nerror,&nerrorall,1,MPI_INT,MPI_SUM,world); if (nerrorall && comm->me == 0) { char str[128]; sprintf(str,"Too many neighbors in CNA for %d atoms",nerrorall); error->warning(FLERR,str,0); } // compute CNA for each atom in group // only performed if # of nearest neighbors = 12 or 14 (fcc,hcp) nerror = 0; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; if (!(mask[i] & groupbit)) { pattern[i] = UNKNOWN; continue; } if (nnearest[i] != 12 && nnearest[i] != 14) { pattern[i] = OTHER; continue; } // loop over near neighbors of I to build cna data structure // cna[k][NCOMMON] = # of common neighbors of I with each of its neighs // cna[k][NBONDS] = # of bonds between those common neighbors // cna[k][MAXBOND] = max # of bonds of any common neighbor // cna[k][MINBOND] = min # of bonds of any common neighbor for (m = 0; m < nnearest[i]; m++) { j = nearest[i][m]; // common = list of neighbors common to atom I and atom J // if J is an owned atom, use its near neighbor list to find them // if J is a ghost atom, use full neighbor list of I to find them // in latter case, must exclude J from I's neighbor list if (j < nlocal) { firstflag = 1; ncommon = 0; for (inear = 0; inear < nnearest[i]; inear++) for (jnear = 0; jnear < nnearest[j]; jnear++) if (nearest[i][inear] == nearest[j][jnear]) { if (ncommon < MAXCOMMON) common[ncommon++] = nearest[i][inear]; else if (firstflag) { nerror++; firstflag = 0; } } } else { xtmp = x[j][0]; ytmp = x[j][1]; ztmp = x[j][2]; jlist = firstneigh[i]; jnum = numneigh[i]; n = 0; for (kk = 0; kk < jnum; kk++) { k = jlist[kk]; k &= NEIGHMASK; if (k == j) continue; delx = xtmp - x[k][0]; dely = ytmp - x[k][1]; delz = ztmp - x[k][2]; rsq = delx*delx + dely*dely + delz*delz; if (rsq < cutsq) { if (n < MAXNEAR) onenearest[n++] = k; else break; } } firstflag = 1; ncommon = 0; for (inear = 0; inear < nnearest[i]; inear++) for (jnear = 0; (jnear < n) && (n < MAXNEAR); jnear++) if (nearest[i][inear] == onenearest[jnear]) { if (ncommon < MAXCOMMON) common[ncommon++] = nearest[i][inear]; else if (firstflag) { nerror++; firstflag = 0; } } } cna[m][NCOMMON] = ncommon; // calculate total # of bonds between common neighbor atoms // also max and min # of common atoms any common atom is bonded to // bond = pair of atoms within cutoff for (n = 0; n < ncommon; n++) bonds[n] = 0; nbonds = 0; for (jj = 0; jj < ncommon; jj++) { j = common[jj]; xtmp = x[j][0]; ytmp = x[j][1]; ztmp = x[j][2]; for (kk = jj+1; kk < ncommon; kk++) { k = common[kk]; delx = xtmp - x[k][0]; dely = ytmp - x[k][1]; delz = ztmp - x[k][2]; rsq = delx*delx + dely*dely + delz*delz; if (rsq < cutsq) { nbonds++; bonds[jj]++; bonds[kk]++; } } } cna[m][NBOND] = nbonds; maxbonds = 0; minbonds = MAXCOMMON; for (n = 0; n < ncommon; n++) { maxbonds = MAX(bonds[n],maxbonds); minbonds = MIN(bonds[n],minbonds); } cna[m][MAXBOND] = maxbonds; cna[m][MINBOND] = minbonds; } // detect CNA pattern of the atom nfcc = nhcp = nbcc4 = nbcc6 = nico = 0; pattern[i] = OTHER; if (nnearest[i] == 12) { for (inear = 0; inear < 12; inear++) { cj = cna[inear][NCOMMON]; ck = cna[inear][NBOND]; cl = cna[inear][MAXBOND]; cm = cna[inear][MINBOND]; if (cj == 4 && ck == 2 && cl == 1 && cm == 1) nfcc++; else if (cj == 4 && ck == 2 && cl == 2 && cm == 0) nhcp++; else if (cj == 5 && ck == 5 && cl == 2 && cm == 2) nico++; } if (nfcc == 12) pattern[i] = FCC; else if (nfcc == 6 && nhcp == 6) pattern[i] = HCP; else if (nico == 12) pattern[i] = ICOS; } else if (nnearest[i] == 14) { for (inear = 0; inear < 14; inear++) { cj = cna[inear][NCOMMON]; ck = cna[inear][NBOND]; cl = cna[inear][MAXBOND]; cm = cna[inear][MINBOND]; if (cj == 4 && ck == 4 && cl == 2 && cm == 2) nbcc4++; else if (cj == 6 && ck == 6 && cl == 2 && cm == 2) nbcc6++; } if (nbcc4 == 6 && nbcc6 == 8) pattern[i] = BCC; } } // warning message MPI_Allreduce(&nerror,&nerrorall,1,MPI_INT,MPI_SUM,world); if (nerrorall && comm->me == 0) { char str[128]; sprintf(str,"Too many common neighbors in CNA %d times",nerrorall); error->warning(FLERR,str); } } /* ---------------------------------------------------------------------- memory usage of local atom-based array ------------------------------------------------------------------------- */ double ComputeCNAAtom::memory_usage() { double bytes = nmax * sizeof(int); bytes += nmax * MAXNEAR * sizeof(int); bytes += nmax * sizeof(double); return bytes; } diff --git a/src/compute_com_chunk.cpp b/src/compute_com_chunk.cpp index ebf2b30f6..3eb686783 100644 --- a/src/compute_com_chunk.cpp +++ b/src/compute_com_chunk.cpp @@ -1,243 +1,242 @@ /* ---------------------------------------------------------------------- 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 <string.h> #include "compute_com_chunk.h" #include "atom.h" #include "update.h" #include "modify.h" #include "compute_chunk_atom.h" #include "domain.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; enum{ONCE,NFREQ,EVERY}; /* ---------------------------------------------------------------------- */ ComputeCOMChunk::ComputeCOMChunk(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) + Compute(lmp, narg, arg), + idchunk(NULL), masstotal(NULL), massproc(NULL), com(NULL), comall(NULL) { if (narg != 4) error->all(FLERR,"Illegal compute com/chunk command"); array_flag = 1; size_array_cols = 3; size_array_rows = 0; size_array_rows_variable = 1; extarray = 0; // ID of compute chunk/atom int n = strlen(arg[3]) + 1; idchunk = new char[n]; strcpy(idchunk,arg[3]); init(); // chunk-based data nchunk = 1; maxchunk = 0; - massproc = masstotal = NULL; - com = comall = NULL; allocate(); firstflag = massneed = 1; } /* ---------------------------------------------------------------------- */ ComputeCOMChunk::~ComputeCOMChunk() { delete [] idchunk; memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(com); memory->destroy(comall); } /* ---------------------------------------------------------------------- */ void ComputeCOMChunk::init() { int icompute = modify->find_compute(idchunk); if (icompute < 0) error->all(FLERR,"Chunk/atom compute does not exist for compute com/chunk"); cchunk = (ComputeChunkAtom *) modify->compute[icompute]; if (strcmp(cchunk->style,"chunk/atom") != 0) error->all(FLERR,"Compute com/chunk does not use chunk/atom compute"); } /* ---------------------------------------------------------------------- */ void ComputeCOMChunk::setup() { // one-time calculation of per-chunk mass // done in setup, so that ComputeChunkAtom::setup() is already called if (firstflag && cchunk->idsflag == ONCE) { compute_array(); firstflag = massneed = 0; } } /* ---------------------------------------------------------------------- */ void ComputeCOMChunk::compute_array() { int index; double massone; double unwrap[3]; invoked_array = update->ntimestep; // compute chunk/atom assigns atoms to chunk IDs // extract ichunk index vector from compute // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms nchunk = cchunk->setup_chunks(); cchunk->compute_ichunk(); int *ichunk = cchunk->ichunk; if (nchunk > maxchunk) allocate(); size_array_rows = nchunk; // zero local per-chunk values for (int i = 0; i < nchunk; i++) com[i][0] = com[i][1] = com[i][2] = 0.0; if (massneed) for (int i = 0; i < nchunk; i++) massproc[i] = 0.0; // compute COM for each chunk double **x = atom->x; int *mask = atom->mask; int *type = atom->type; imageint *image = atom->image; double *mass = atom->mass; double *rmass = atom->rmass; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { index = ichunk[i]-1; if (index < 0) continue; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; domain->unmap(x[i],image[i],unwrap); com[index][0] += unwrap[0] * massone; com[index][1] += unwrap[1] * massone; com[index][2] += unwrap[2] * massone; if (massneed) massproc[index] += massone; } MPI_Allreduce(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); if (massneed) MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world); for (int i = 0; i < nchunk; i++) { if (masstotal[i] > 0.0) { comall[i][0] /= masstotal[i]; comall[i][1] /= masstotal[i]; comall[i][2] /= masstotal[i]; } else comall[i][0] = comall[i][1] = comall[i][2] = 0.0; } } /* ---------------------------------------------------------------------- lock methods: called by fix ave/time these methods insure vector/array size is locked for Nfreq epoch by passing lock info along to compute chunk/atom ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- increment lock counter ------------------------------------------------------------------------- */ void ComputeCOMChunk::lock_enable() { cchunk->lockcount++; } /* ---------------------------------------------------------------------- decrement lock counter in compute chunk/atom, it if still exists ------------------------------------------------------------------------- */ void ComputeCOMChunk::lock_disable() { int icompute = modify->find_compute(idchunk); if (icompute >= 0) { cchunk = (ComputeChunkAtom *) modify->compute[icompute]; cchunk->lockcount--; } } /* ---------------------------------------------------------------------- calculate and return # of chunks = length of vector/array ------------------------------------------------------------------------- */ int ComputeCOMChunk::lock_length() { nchunk = cchunk->setup_chunks(); return nchunk; } /* ---------------------------------------------------------------------- set the lock from startstep to stopstep ------------------------------------------------------------------------- */ void ComputeCOMChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) { cchunk->lock(fixptr,startstep,stopstep); } /* ---------------------------------------------------------------------- unset the lock ------------------------------------------------------------------------- */ void ComputeCOMChunk::unlock(Fix *fixptr) { cchunk->unlock(fixptr); } /* ---------------------------------------------------------------------- free and reallocate per-chunk arrays ------------------------------------------------------------------------- */ void ComputeCOMChunk::allocate() { memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(com); memory->destroy(comall); maxchunk = nchunk; memory->create(massproc,maxchunk,"com/chunk:massproc"); memory->create(masstotal,maxchunk,"com/chunk:masstotal"); memory->create(com,maxchunk,3,"com/chunk:com"); memory->create(comall,maxchunk,3,"com/chunk:comall"); array = comall; } /* ---------------------------------------------------------------------- memory usage of local data ------------------------------------------------------------------------- */ double ComputeCOMChunk::memory_usage() { double bytes = (bigint) maxchunk * 2 * sizeof(double); bytes += (bigint) maxchunk * 2*3 * sizeof(double); return bytes; } diff --git a/src/compute_contact_atom.cpp b/src/compute_contact_atom.cpp index 99c4766b6..c86829bc1 100644 --- a/src/compute_contact_atom.cpp +++ b/src/compute_contact_atom.cpp @@ -1,197 +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) + Compute(lmp, narg, arg), + contact(NULL) { 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(this,instance_me); 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); 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 m; } /* ---------------------------------------------------------------------- */ 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_coord_atom.cpp b/src/compute_coord_atom.cpp index c6aa561eb..5747c869d 100644 --- a/src/compute_coord_atom.cpp +++ b/src/compute_coord_atom.cpp @@ -1,227 +1,226 @@ /* ---------------------------------------------------------------------- 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_coord_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; /* ---------------------------------------------------------------------- */ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) + Compute(lmp, narg, arg), + typelo(NULL), typehi(NULL), cvec(NULL), carray(NULL) { if (narg < 4) error->all(FLERR,"Illegal compute coord/atom command"); double cutoff = force->numeric(FLERR,arg[3]); cutsq = cutoff*cutoff; ncol = narg-4 + 1; int ntypes = atom->ntypes; typelo = new int[ncol]; typehi = new int[ncol]; if (narg == 4) { ncol = 1; typelo[0] = 1; typehi[0] = ntypes; } else { ncol = 0; int iarg = 4; while (iarg < narg) { force->bounds(arg[iarg],ntypes,typelo[ncol],typehi[ncol]); if (typelo[ncol] > typehi[ncol]) error->all(FLERR,"Illegal compute coord/atom command"); ncol++; iarg++; } } peratom_flag = 1; if (ncol == 1) size_peratom_cols = 0; else size_peratom_cols = ncol; nmax = 0; - cvec = NULL; - carray = NULL; } /* ---------------------------------------------------------------------- */ ComputeCoordAtom::~ComputeCoordAtom() { delete [] typelo; delete [] typehi; memory->destroy(cvec); memory->destroy(carray); } /* ---------------------------------------------------------------------- */ void ComputeCoordAtom::init() { if (force->pair == NULL) error->all(FLERR,"Compute coord/atom requires a pair style be defined"); if (sqrt(cutsq) > force->pair->cutforce) error->all(FLERR, "Compute coord/atom cutoff is longer than pairwise cutoff"); // need an occasional full neighbor list 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,"coord/atom") == 0) count++; if (count > 1 && comm->me == 0) error->warning(FLERR,"More than one compute coord/atom"); } /* ---------------------------------------------------------------------- */ void ComputeCoordAtom::init_list(int id, NeighList *ptr) { list = ptr; } /* ---------------------------------------------------------------------- */ void ComputeCoordAtom::compute_peratom() { int i,j,m,ii,jj,inum,jnum,jtype,n; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *ilist,*jlist,*numneigh,**firstneigh; double *count; invoked_peratom = update->ntimestep; // grow coordination array if necessary if (atom->nmax > nmax) { if (ncol == 1) { memory->destroy(cvec); nmax = atom->nmax; memory->create(cvec,nmax,"coord/atom:cvec"); vector_atom = cvec; } else { memory->destroy(carray); nmax = atom->nmax; memory->create(carray,nmax,ncol,"coord/atom:carray"); array_atom = carray; } } // invoke full neighbor list (will copy or build if necessary) neighbor->build_one(list); inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; // compute coordination number(s) for each atom in group // use full neighbor list to count atoms less than cutoff double **x = atom->x; int *type = atom->type; int *mask = atom->mask; if (ncol == 1) { 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]; jlist = firstneigh[i]; jnum = numneigh[i]; n = 0; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; j &= NEIGHMASK; jtype = type[j]; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; if (rsq < cutsq && jtype >= typelo[0] && jtype <= typehi[0]) n++; } cvec[i] = n; } else cvec[i] = 0.0; } } else { for (ii = 0; ii < inum; ii++) { i = ilist[ii]; count = carray[i]; for (m = 0; m < ncol; m++) count[m] = 0.0; if (mask[i] & groupbit) { xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; jlist = firstneigh[i]; jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; j &= NEIGHMASK; jtype = type[j]; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; if (rsq < cutsq) { for (m = 0; m < ncol; m++) if (jtype >= typelo[m] && jtype <= typehi[m]) count[m] += 1.0; } } } } } } /* ---------------------------------------------------------------------- memory usage of local atom-based array ------------------------------------------------------------------------- */ double ComputeCoordAtom::memory_usage() { double bytes = ncol*nmax * sizeof(double); return bytes; }