Hi all,
I start to read the codes of LAMMPS recently and I have a question when I’m reading pair_airebo.cpp. As is known, a new neighbor list is built from the main neighbor list to calculate the REBO potential in pair_airebo.cpp. This neighbor list is built in REBO_neigh(). To understand this subroutine, I build two layers of graphene with AB stacking as the initial configuration to do some tests. I set periodic boundary condition both in x and y directions. What I want to know is how many neighbors for each atom. To do this, I add a line in pair_airebo.cpp (Bold) to output the results. The subroutine REBO_neigh() is like follows:
/* ----------------------------------------------------------------------
create REBO neighbor list from main neighbor list
REBO neighbor list stores neighbors of ghost atoms
------------------------------------------------------------------------- */
void PairAIREBO::REBO_neigh()
{
int i,j,ii,jj,n,allnum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq,dS;
int *ilist,*jlist,*numneigh,**firstneigh;
int *neighptr;
double **x = atom->x;
int *type = atom->type;
if (atom->nmax > maxlocal) {
maxlocal = atom->nmax;
memory->destroy(REBO_numneigh);
memory->sfree(REBO_firstneigh);
memory->destroy(nC);
memory->destroy(nH);
memory->create(REBO_numneigh,maxlocal,“AIREBO:numneigh”);
REBO_firstneigh = (int **) memory->smalloc(maxlocal*sizeof(int *),
“AIREBO:firstneigh”);
memory->create(nC,maxlocal,“AIREBO:nC”);
memory->create(nH,maxlocal,“AIREBO:nH”);
}
allnum = list->inum + list->gnum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// store all REBO neighs of owned and ghost atoms
// scan full neighbor list of I
ipage->reset();
for (ii = 0; ii < allnum; ii++) {
i = ilist[ii];
n = 0;
neighptr = ipage->vget();
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = map[type[i]];
nC[i] = nH[i] = 0.0;
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = map[type[j]];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delxdelx + delydely + delz*delz;
if (rsq < rcmaxsq[itype][jtype]) {
neighptr[n++] = j;
if (jtype == 0)
nC[i] += Sp(sqrt(rsq),rcmin[itype][jtype],rcmax[itype][jtype],dS);
else
nH[i] += Sp(sqrt(rsq),rcmin[itype][jtype],rcmax[itype][jtype],dS);
}
}
REBO_firstneigh[i] = neighptr;
REBO_numneigh[i] = n;
fprintf (screen,"%d %d %d %d %d %d %d %d %d %f %f %f %f\n",allnum,list->inum,n,i,j,type[i],type[j],itype,jtype,rsq,delx,dely,delz);
ipage->vgot(n);
if (ipage->status())
error->one(FLERR,“Neighbor list overflow, boost neigh_modify one”);
}
The total number of atoms in my configuration is 1760 and the number of neighbors for each atom is n. If I understand correctly, allnum is the number of atoms on each proc plus its ghost atoms.
I run lammps with one core so inum=1760. I found that when i&j < 1760, n = 3. While for i,j > 1760, the typical output is like follows:
allnum,list->inum,n,i,j,type[i],type[j],itype,jtype,rsq,delx,dely,delz
4086 1760 2 1976 2028 2 2 0 0 5.861532 0.000000 -2.421060 0.000000
4086 1760 1 2016 2033 1 2 0 0 13.310768 0.698900 -1.210540 -3.370000
4086 1760 1 2021 2081 1 1 0 0 5.861532 0.000000 -2.421060 0.000000
From these output, it seems that the number of neighbors for the gohst atom is less than 3. But this result confuses me since I set PBC in x and y directions and the number of neighbors for each atom should be always 3.
Maybe I misunderstand the meaning of ghost atom? Could anybody tell me what’s the reason?
Any input in greatly appreciated. Thank you very much!
Best,
Huang