# [lammps-users] centro-symmetry parameter and self-neighbor issue

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 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;
}
}

Craig,

Looking at your code, I think that each atom should not be in its own neighbor list and your assert should not fail. So, no, I don’t think it is safe to assume each atom will always appear in its own neighbor list as the first entry. My guess is that there’s still some subtle bug in your modified version of the code.

Paul

No neighbor list that LAMMPS produces should have itself in it.
If the cutoff is long and the box is small, it might have an image
of itself in it, but that would be as a ghost atom.

Steve