Homec4science
Diffusion lammps (master) Browse

lammps/src/USER-MISC/compute_cnp_atom.cpp
9a7207e34cb4master

Last Change
Raw File

compute_cnp_atom.cpp
View Options

/* ----------------------------------------------------------------------
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.
Common Neighbor Parameter as proposed in:
Tsuzuki, Branicio, Rino, Comput Phys Comm, 177, 518 (2007)
Cite: http://dx.doi.org/10.1063/1.2197987
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Paulo Branicio (University of Southern California)
branicio@usc.edu
------------------------------------------------------------------------- */
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include "compute_cnp_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;
//define maximum values
#define MAXNEAR 24
#define MAXCOMMON 12
enum{NCOMMON};
/* ---------------------------------------------------------------------- */
ComputeCNPAtom::ComputeCNPAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
nearest(NULL), nnearest(NULL), cnpv(NULL)
{
if (narg != 4) error->all(FLERR,"Illegal compute cnp/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 cnp/atom command");
cutsq = cutoff*cutoff;
// apply check for single type atoms in compute group
int lasttype = -1;
int n = -1;
for (int i=0; i < atom->nlocal; ++i) {
if (atom->mask[i] & groupbit) {
if (lasttype != atom->type[i]) {
lasttype = atom->type[i];
++n;
}
}
}
int all_n = 0;
MPI_Allreduce(&n,&all_n,1,MPI_INT,MPI_MAX,world);
if (all_n > 0)
error->warning(FLERR,"Compute cnp/atom requested on multi-type system");
nmax = 0;
}
/* ---------------------------------------------------------------------- */
ComputeCNPAtom::~ComputeCNPAtom()
{
memory->destroy(nearest);
memory->destroy(nnearest);
memory->destroy(cnpv);
}
/* ---------------------------------------------------------------------- */
void ComputeCNPAtom::init()
{
if (force->pair == NULL)
error->all(FLERR,"Compute cnp/atom requires a pair style be defined");
if (sqrt(cutsq) > force->pair->cutforce)
error->all(FLERR,"Compute cnp/atom cutoff is longer than pairwise cutoff");
if (2.0*sqrt(cutsq) > force->pair->cutforce + neighbor->skin &&
comm->me == 0)
error->warning(FLERR,"Compute cnp/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,"cnp/atom") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute cnp/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 ComputeCNPAtom::init_list(int id, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
void ComputeCNPAtom::compute_peratom()
{
int i,j,k,ii,jj,kk,m,n,inum,jnum,inear,jnear;
int firstflag,ncommon;
int *ilist,*jlist,*numneigh,**firstneigh;
int onenearest[MAXNEAR];
int common[MAXCOMMON];
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
double xjtmp,yjtmp,zjtmp,rjkx,rjky,rjkz;
invoked_peratom = update->ntimestep;
// grow arrays if necessary
if (atom->nmax > nmax) {
memory->destroy(nearest);
memory->destroy(nnearest);
memory->destroy(cnpv);
nmax = atom->nmax;
memory->create(nearest,nmax,MAXNEAR,"cnp:nearest");
memory->create(nnearest,nmax,"cnp:nnearest");
memory->create(cnpv,nmax,"cnp:cnp_cnpv");
vector_atom = cnpv;
}
// 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 neigbors 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 CNP 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 CNP for %d atoms",nerrorall);
error->warning(FLERR,str,0);
}
// compute CNP value 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];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
// reset cnpv
cnpv[i] = 0.0;
// skip computation of cnpv for atoms outside the compute group
if (!(mask[i] & groupbit)) continue;
// loop over nearest neighbors of I to build cnp data structure
// cnp[k][NCOMMON] = # of common neighbors of I with each of its neighbors
for (m = 0; m < nnearest[i]; m++) {
j = nearest[i][m];
xjtmp = x[j][0];
yjtmp = x[j][1];
zjtmp = x[j][2];
// 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
// find common neighbors of i and j using near 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;
}
}
// find common neighbors of i and j using full neighbor list
} else {
jlist = firstneigh[i];
jnum = numneigh[i];
n = 0;
for (kk = 0; kk < jnum; kk++) {
k = jlist[kk];
k &= NEIGHMASK;
if (k == j) continue;
delx = xjtmp - x[k][0];
dely = yjtmp - x[k][1];
delz = zjtmp - 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;
}
}
}
// Calculate and update sum |Rik+Rjk|ˆ2
rjkx = 0.0;
rjky = 0.0;
rjkz = 0.0;
for (kk = 0; kk < ncommon; kk++) {
k = common[kk];
rjkx += 2.0*x[k][0] - xjtmp - xtmp;
rjky += 2.0*x[k][1] - yjtmp - ytmp;
rjkz += 2.0*x[k][2] - zjtmp - ztmp;
}
// update cnpv with summed (valuejk)
cnpv[i] += rjkx*rjkx + rjky*rjky + rjkz*rjkz;
// end of loop over j atoms
}
// normalize cnp by the number of nearest neighbors
cnpv[i] = cnpv[i] / nnearest[i];
// end of loop over i atoms
}
// 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 CNP %d times",nerrorall);
error->warning(FLERR,str);
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeCNPAtom::memory_usage()
{
double bytes = nmax * sizeof(int);
bytes += nmax * MAXNEAR * sizeof(int);
bytes += nmax * sizeof(double);
return bytes;
}
c4science · Help