Hi all.
I’m trying to write a compute which is very close in spirit to compute_centro_atom.cpp.
In the ComputeCentro code, for each atom, the 12 nearest neighbors are identified and then used in a subsequent calculation. The problem is that it seems like the neighbor list for the current atom includes the atom itself, which seems to me like it would make the calculation for the centro symmetry parameter wrong. e.g. the assert line in the final loop of my code (below) fails.
I could just look for the n+1 nearest neighbors and ignore the first one assuming that the atom itself will always appear in its own neighbor list. But is this a safe assumption?
Thanks,
Craig
-------------------------- code snippet -----------------------------
void ComputeQ6Atom::compute_peratom()
{
int i,j,k,ii,jj,kk,n,inum,jnum;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq,value;
int *ilist,*jlist,*numneigh,**firstneigh;
double pairs[66];
invoked |= INVOKED_PERATOM;
// grow q6 array if necessary
if (atom->nlocal > nmax) {
memory->destroy_2d_double_array(q6);
nmax = atom->nmax;
q6 = memory->create_2d_double_array(nmax,2,“compute/q6/atom:q6”);
vector_atom = q6;
}
// 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;
// compute q6 parameter for each atom in group
// use full neighbor list
double **x = atom->x;
int *mask = atom->mask;
int nall = atom->nlocal + atom->nghost;
double cutsq = force->pair->cutforce * force->pair->cutforce;
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];
// insure distsq and nearest arrays are long enough
if (jnum > maxneigh) {
memory->sfree(distsq);
memory->sfree(nearest);
maxneigh = jnum;
distsq = (double ) memory->smalloc(maxneighsizeof(double),
“compute/q6/atom:q6”);
nearest = (int ) memory->smalloc(maxneighsizeof(int),
“compute/q6/atom:q6”);
}
// loop over list of all neighbors within force cutoff
// distsq[] = distance sq to each
// nearest[] = atom indices of neighbors
n = 0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
if (j >= nall) j %= nall;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delxdelx + delydely + delz*delz;
if (rsq < cutsq) {
distsq[n] = rsq;
nearest[n++] = j;
}
}
// if not at least 6 neighbors, q6 = 0.0
if (n < 6) {
q6[i][0] = 0.0;
q6[i][1] = 0.0;
continue;
}
// store 6 nearest neighs in 1st 6 locations of distsq and nearest
select2(6,n,distsq,nearest);
// R = Ri + Rj for each of 66 i,j pairs among 12 neighbors
// pairs = squared length of each R
double q6Re=0.0;
double q6Im=0.0;
for (j=0; j<6; j++){
assert(nearest[j]!=i);
double dx=x[nearest[j]][0]-x[i][0];
double dy=x[nearest[j]][1]-x[i][1];
double r=sqrt(dxdx+dydy);
printf("%f\t%f\n",dx,dy);
q6Re+=dx/r;
q6Im+=dy/r;
}
q6[i][0]=q6Re;
q6[i][1]=q6Im;
} else q6[i][0]=q6[i][1] = 0.0;
}
}