git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@2498 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
366
src/compute_cna_atom.cpp
Normal file
366
src/compute_cna_atom.cpp
Normal file
@ -0,0 +1,366 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
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 MAX(a,b) ((a) > (b) ? (a) : (b))
|
||||
#define MIN(a,b) ((a) < (b) ? (a) : (b))
|
||||
|
||||
#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)
|
||||
{
|
||||
if (narg != 4) error->all("Illegal compute cna/atom command");
|
||||
|
||||
peratom_flag = 1;
|
||||
size_peratom = 0;
|
||||
|
||||
double cutoff = atof(arg[3]);
|
||||
if (cutoff < 0.0) error->all("Illegal compute cna/atom command");
|
||||
cutsq = cutoff*cutoff;
|
||||
|
||||
nmax = 0;
|
||||
nearest = NULL;
|
||||
nnearest = NULL;
|
||||
pattern = NULL;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeCNAAtom::~ComputeCNAAtom()
|
||||
{
|
||||
memory->destroy_2d_int_array(nearest);
|
||||
memory->sfree(nnearest);
|
||||
memory->sfree(pattern);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeCNAAtom::init()
|
||||
{
|
||||
if (force->pair == NULL)
|
||||
error->all("Compute cna/atom requires a pair style be defined");
|
||||
if (sqrt(cutsq) > force->pair->cutforce)
|
||||
error->all("Compute cna/atom cutoff is longer than pairwise cutoff");
|
||||
|
||||
double cutneighbor = force->pair->cutforce + neighbor->skin;
|
||||
if (2.0*sqrt(cutsq) > cutneighbor && comm->me == 0)
|
||||
error->warning("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("More than one compute cna/atom defined");
|
||||
|
||||
// need an occasional full neighbor list
|
||||
|
||||
int irequest = neighbor->request((void *) this);
|
||||
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->nlocal > nmax) {
|
||||
memory->destroy_2d_int_array(nearest);
|
||||
memory->sfree(nnearest);
|
||||
memory->sfree(pattern);
|
||||
nmax = atom->nmax;
|
||||
|
||||
nearest = memory->create_2d_int_array(nmax,MAXNEAR,"cna:nearest");
|
||||
nnearest = (int *) memory->smalloc(nmax*sizeof(int),"cna:nnearest");
|
||||
pattern = (double *) memory->smalloc(nmax*sizeof(double),
|
||||
"cna:cna_pattern");
|
||||
scalar_atom = pattern;
|
||||
}
|
||||
|
||||
// invoke full 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;
|
||||
|
||||
// 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];
|
||||
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(str);
|
||||
}
|
||||
|
||||
// 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];
|
||||
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; 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(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;
|
||||
}
|
||||
41
src/compute_cna_atom.h
Normal file
41
src/compute_cna_atom.h
Normal file
@ -0,0 +1,41 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifndef COMPUTE_CNA_ATOM_H
|
||||
#define COMPUTE_CNA_ATOM_H
|
||||
|
||||
#include "compute.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class ComputeCNAAtom : public Compute {
|
||||
public:
|
||||
ComputeCNAAtom(class LAMMPS *, int, char **);
|
||||
~ComputeCNAAtom();
|
||||
void init();
|
||||
void init_list(int, class NeighList *);
|
||||
void compute_peratom();
|
||||
double memory_usage();
|
||||
|
||||
private:
|
||||
int nmax;
|
||||
double cutsq;
|
||||
class NeighList *list;
|
||||
int **nearest;
|
||||
int *nnearest;
|
||||
double *pattern;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
Reference in New Issue
Block a user