Hi everyone,

I am trying to implement a new pair style which requires me to keep track of a number of 1D per-atom arrays. I am using an old-ish release, 7Aug19. but I think that should not be an issue given the simplicity of what I’m trying to do.

If it is of interest to anyone, I am trying to code a double-sum implementation of the 3-body Stillinger-Weber potential (pair_sw)

I am using the pair_eam style as a guide to work out how to communicate and update the per-atom arrays that I have defined. I have also gone through the mail list archive but not exhaustively; forgive me if I’ve missed something obvious.

My issue is that atoms of a given “tag” show up multiple times as ghost atom neighbors of a given atom. I need to understand how to filter out/choose and include only the “correct” pairs from these.

These ghost atom neighbors have different distances with respect to the reference (ilist[i]) atom and I find my code computing the interaction between a pair (itag,jtag) multiple times.

I have a self-written code of the pair style where I calculate the interaction due to any pair once using the nearest image distance in a cubic box with PBC. I am using the numbers from that code as a reference for testing.

Details on test:

N=4 (need more than 2 because the potential is a 3-body potential).

run 1 step and observe total energy and forces on each particle

box-length ~2.2163 with a cut-off of ~3.75

Results compared with those from pair_sw on an identical configuration. Checking if total energy and forces on each particle are the same.

I’ve attached a code snippet with (I think) the main details needed. I have a number of per-atom arrays but I’ve written out only two to get the point across

I’ve only included the compute (with allocation of the new arrays), pack/unpack_reverse_comm routines. There is a debug line where I print out for a given i,j, pair (i, j taken from list->ilist and jlist), the tags and the distance rij.

I’ve omitted details specific to how the potential is calculated, as well as other debug lines, expecting that they are not relevant to the issue.

void PairSWDS::compute(int eflag, int vflag)

{

int i,j,k,ii,jj,kk,inum,jnum,jnumm1;

int itype,jtype,ktype,ijparam,ikparam,ijkparam;

tagint itag,jtag;

double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;

double rsq,rsq1,rsq2,rp,rq;

double delr1[3],delr2[3],fj[3],fk[3];

int *ilist,*jlist,*numneigh,**firstneigh;

evdwl = 0.0;

if (eflag || vflag) ev_setup(eflag,vflag);

else evflag = vflag_fdotr = 0;

//printf(“computing DS SW\n”);

double **x = atom->x;

double **f = atom->f;

tagint *tag = atom->tag;

int *type = atom->type;

nlocal = atom->nlocal;

nall = nlocal + atom->nghost;

newton_pair = force->newton_pair;

printf(“newton pair is %d\n”,newton_pair);

inum = list->inum;

ilist = list->ilist;

numneigh = list->numneigh;

firstneigh = list->firstneigh;

memory->destroy(hi);

memory->destroy(si);

memory->create(hi,nall); // oversized if newton_pair=0

memory->create(si,nall);

int looplimit;

if (newton_pair) looplimit=nall;

else looplimit = nlocal;

for(i=0;i<looplimit;i++){

hi[i]=0.0;

si[i]=0.0;

}

for (ii = 0; ii < inum; ii++) {

i = ilist[ii];

//printf(“hi i, i %f %d\n”,hi[i],i);

itag = tag[i];

itype = map[type[i]];

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;

jtag = tag[j];

if (itag > jtag) {

if ((itag+jtag) 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) 2 == 1) continue;

} else {

if (x[j][2] < ztmp) continue;

if (x[j][2] == ztmp && x[j][1] < ytmp) continue;

if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;

//if (itag==jtag) continue;

}

jtype = map[type[j]];

delx = xtmp - x[j][0];

dely = ytmp - x[j][1];

delz = ztmp - x[j][2];

rsq = delx*delx + dely*dely + delz*delz;

ijparam = elem2param[itype][jtype][jtype];

if (rsq > params[ijparam].cutsq ) continue; // apt cut-off for given pair of atoms, their coeffs

rij=sqrt(rsq);

// Note : debug print below

printf("%d %d i j %d %d itag jtag %f rij \n",i,j,itag,jtag,rij);

hi[i] += exp(rsq);

si[i] += 2*rsq;
// when do we update these?
hi[j] += exp(rsq);
si[j] += 2*rsq;

}

}

if (!newton_pair) comm->reverse_comm_pair(this); // comm needed if supplied half-neighbor list with newton_pair=0

for (ii = 0; ii < inum; ii++) {

i = ilist[ii];

//printf(“hi i, i %f %d\n”,hi[i],i);

itag = tag[i];

itype = map[type[i]];

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;

jtag = tag[j];

if (itag > jtag) {

if ((itag+jtag) 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) 2 == 1) continue;

} else {

if (x[j][2] < ztmp) continue;

if (x[j][2] == ztmp && x[j][1] < ytmp) continue;

if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;

//if (itag==jtag) continue;

}

jtype = map[type[j]];

delx = xtmp - x[j][0];

dely = ytmp - x[j][1];

delz = ztmp - x[j][2];

rsq = delx*delx + dely*dely + delz*delz;

ijparam = elem2param[itype][jtype][jtype];

if (rsq > params[ijparam].cutsq ) continue; // apt cut-off for given pair of atoms, their coeffs

evdwl = 1/rsq + hi[i] - hi[j] + si[i] + si[j];

fijx = (hi[i] - hi[j] + si[i] + si[j])*delx/sqrt(rsq);

fijy = (hi[i] - hi[j] + si[i] + si[j])*dely/sqrt(rsq);

fijz = (hi[i] - hi[j] + si[i] + si[j])*delz/sqrt(rsq);

f[i][0] += fijx;

f[i][1] += fijy;

f[i][2] += fijz;

f[j][0] -= fijx;

f[j][1] -= fijy;

f[j][2] -= fijz;

if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,evdwl,0.0,

fijx,fijy,fijz,delx,dely,delz);

}

}

}

int PairSWDS::pack_reverse_comm(int n, int first, double *buf)
{
int ii,i,m,last;
int count=0;
m = 0; // buffer index counter - buffer entries not indexed with i directly
last = first + n;
for (ii = first; ii < last; ii++) // loops over indices split to the different procs
{
i=ii-first;
buf[m++] = hi[i];
buf[m++] = si[i];
}
// for(i=first;i<last;i++)
// {
// printf(“packing buffer %d %d %f %f\n”,i,i-first,hi[i-first],buf[13*(i-first)]);

// }

return m; // size of the full buffer that is communicated

}

/* ---------------------------------------------------------------------- */

void PairSWDS::unpack_reverse_comm(int n, int *list, double *buf)

{

int ii,i,j,m;

m = 0; // final value is size of communicated buffer - entries in buffer are not indexed with i

printf(“unpacking %d\n”,n);

for (ii = 0; ii < n; ii++) // loops over indices provided from the other processors

{

j = list[ii];

//printf("%d %d %d %f %f\n",i,j,m,buf[m],hi[j]);

hi[j] += buf[m++];

si[j] += buf[m++];

}

}

// output of debug line - N=4 atoms - box size 2.2163 (all lengths Angstrom)

0 1 i j 1 2 itag jtag 2.494940 rij

0 2 i j 1 3 itag jtag 3.534513 rij <<

0 13 i j 1 2 itag jtag 2.435488 rij <<

0 14 i j 1 3 itag jtag 3.590037 rij

0 44 i j 1 3 itag jtag 3.545202 rij

0 45 i j 1 4 itag jtag 2.278151 rij <<

0 49 i j 1 4 itag jtag 2.572599 rij

0 56 i j 1 3 itag jtag 3.600561 rij

1 24 i j 2 1 itag jtag 2.435488 rij <<

1 31 i j 2 4 itag jtag 3.720288 rij

1 0 i j 2 1 itag jtag 2.494940 rij

1 2 i j 2 3 itag jtag 2.589012 rij

1 6 i j 2 3 itag jtag 2.509788 rij <<

1 73 i j 2 4 itag jtag 3.359584 rij <<

1 49 i j 2 4 itag jtag 3.556079 rij

2 108 i j 3 1 itag jtag 3.600561 rij

2 84 i j 3 1 itag jtag 3.545202 rij

2 24 i j 3 1 itag jtag 3.590037 rij

2 27 i j 3 4 itag jtag 2.492355 rij <<

2 9 i j 3 2 itag jtag 2.509788 rij

2 0 i j 3 1 itag jtag 3.534513 rij <<

2 1 i j 3 2 itag jtag 2.589012 rij

2 3 i j 3 4 itag jtag 2.623361 rij

3 92 i j 4 1 itag jtag 2.572599 rij

3 93 i j 4 2 itag jtag 3.556079 rij

3 84 i j 4 1 itag jtag 2.278151 rij <<

3 105 i j 4 2 itag jtag 3.359584 rij <<

3 2 i j 4 3 itag jtag 2.623361 rij

3 21 i j 4 2 itag jtag 3.720288 rij

3 14 i j 4 3 itag jtag 2.492355 rij <<

The lines marked with “<<” are the “correct” computations, i.e., I expect that the others will be filtered out with the use of appropriate “if” conditions to give me the correct result.

Apologies for the long post, I will be grateful for any help.

Best regards,

Yagyik