bugfix for 2 recenty reported neighbor issues, also a ReaxFF fix species issue

This commit is contained in:
Steve Plimpton 2017-05-16 15:51:41 -06:00
parent 06c151421c
commit 7f9a331c73
47 changed files with 80 additions and 98 deletions

View File

@ -35,7 +35,7 @@ cutoff.
In contrast to "pair_style yukawa"_pair_yukawa.html, this functional
form arises from the Coulombic interaction between two colloid
particles, screened due to the presence of an electrolyte, see the
book by "Safran"_#Safran for a derivation in the context of DVLO
book by "Safran"_#Safran for a derivation in the context of DLVO
theory. "Pair_style yukawa"_pair_yukawa.html is a screened Coulombic
potential between two point-charges and uses no such approximation.

View File

@ -286,9 +286,6 @@ void FixQEq::setup_pre_force(int vflag)
if (force->newton_pair == 0)
error->all(FLERR,"QEQ with 'newton pair off' not supported");
// should not be needed
// neighbor->build_one(list);
deallocate_storage();
allocate_storage();

View File

@ -76,6 +76,7 @@ void NBinSSA::bin_atoms()
if (ssaAIR[i] < 2) continue; // skip ghost atoms not in AIR
if (mask[i] & bitmask) {
ibin = coord2bin(x[i]);
atom2bin[i] = ibin;
bins_ssa[i] = gbinhead_ssa[ibin];
gbinhead_ssa[ibin] = i;
}
@ -84,12 +85,14 @@ void NBinSSA::bin_atoms()
for (i = nall-1; i >= nlocal; i--) {
if (ssaAIR[i] < 2) continue; // skip ghost atoms not in AIR
ibin = coord2bin(x[i]);
atom2bin[i] = ibin;
bins_ssa[i] = gbinhead_ssa[ibin];
gbinhead_ssa[ibin] = i;
}
}
for (i = nlocal-1; i >= 0; i--) {
ibin = coord2bin(x[i]);
atom2bin[i] = ibin;
bins_ssa[i] = binhead_ssa[ibin];
binhead_ssa[ibin] = i;
}

View File

@ -141,7 +141,7 @@ void NPairHalfBinNewtonSSA::build(NeighList *list)
}
}
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
// loop over all local atoms in other bins in "half" stencil

View File

@ -97,7 +97,7 @@ void NPairFullBinGhostOmp::build(NeighList *list)
// no molecular test when i = ghost atom
if (i < nlocal) {
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (i == j) continue;

View File

@ -90,7 +90,7 @@ void NPairFullBinOmp::build(NeighList *list)
// loop over all atoms in surrounding bins in stencil including self
// skip i = j
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {

View File

@ -94,7 +94,7 @@ void NPairFullMultiOmp::build(NeighList *list)
// skip if i,j neighbor cutoff is less than bin distance
// skip i = j
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
s = stencil_multi[itype];
distsq = distsq_multi[itype];
cutsq = cutneighsq[itype];

View File

@ -103,7 +103,7 @@ void NPairHalfBinNewtoffGhostOmp::build(NeighList *list)
// no molecular test when i = ghost atom
if (i < nlocal) {
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {

View File

@ -94,7 +94,7 @@ void NPairHalfBinNewtoffOmp::build(NeighList *list)
// stores own/own pairs only once
// stores own/ghost pairs on both procs
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {

View File

@ -130,7 +130,7 @@ void NPairHalfBinNewtonOmp::build(NeighList *list)
// loop over all atoms in other bins in stencil, store every pair
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
jtype = type[j];

View File

@ -94,7 +94,7 @@ void NPairHalfBinNewtonTriOmp::build(NeighList *list)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (x[j][2] < ztmp) continue;

View File

@ -97,7 +97,7 @@ void NPairHalfMultiNewtoffOmp::build(NeighList *list)
// stores own/own pairs only once
// stores own/ghost pairs on both procs
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
s = stencil_multi[itype];
distsq = distsq_multi[itype];
cutsq = cutneighsq[itype];

View File

@ -131,7 +131,7 @@ void NPairHalfMultiNewtonOmp::build(NeighList *list)
// loop over all atoms in other bins in stencil, store every pair
// skip if i,j neighbor cutoff is less than bin distance
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
s = stencil_multi[itype];
distsq = distsq_multi[itype];
cutsq = cutneighsq[itype];

View File

@ -99,7 +99,7 @@ void NPairHalfMultiNewtonTriOmp::build(NeighList *list)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
s = stencil_multi[itype];
distsq = distsq_multi[itype];
cutsq = cutneighsq[itype];

View File

@ -117,7 +117,7 @@ void NPairHalfRespaBinNewtoffOmp::build(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];

View File

@ -176,7 +176,7 @@ void NPairHalfRespaBinNewtonOmp::build(NeighList *list)
// loop over all atoms in other bins in stencil, store every pair
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
jtype = type[j];

View File

@ -128,7 +128,7 @@ void NPairHalfRespaBinNewtonTriOmp::build(NeighList *list)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (x[j][2] < ztmp) continue;

View File

@ -113,7 +113,7 @@ void NPairHalfSizeBinNewtoffOmp::build(NeighList *list)
ytmp = x[i][1];
ztmp = x[i][2];
radi = radius[i];
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
// loop over all atoms in surrounding bins in stencil including self
// only store pair if i < j

View File

@ -168,7 +168,7 @@ void NPairHalfSizeBinNewtonOmp::build(NeighList *list)
// loop over all atoms in other bins in stencil, store every pair
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (exclude && exclusion(i,j,type[i],type[j],mask,molecule)) continue;

View File

@ -84,7 +84,7 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (x[j][2] < ztmp) continue;

View File

@ -416,9 +416,6 @@ void FixQEqReax::init_taper()
void FixQEqReax::setup_pre_force(int vflag)
{
// should not be needed
// neighbor->build_one(list);
deallocate_storage();
allocate_storage();

View File

@ -500,7 +500,7 @@ void FixReaxCSpecies::Output_ReaxC_Bonds(bigint ntimestep, FILE *fp)
AtomCoord FixReaxCSpecies::chAnchor(AtomCoord in1, AtomCoord in2)
{
if (in1.x < in2.x)
if (in1.x <= in2.x)
return in1;
return in2;
}

View File

@ -29,6 +29,7 @@ NBin::NBin(LAMMPS *lmp) : Pointers(lmp)
maxbin = maxatom = 0;
binhead = NULL;
bins = NULL;
atom2bin = NULL;
// geometry settings
@ -42,6 +43,7 @@ NBin::~NBin()
{
memory->destroy(binhead);
memory->destroy(bins);
memory->destroy(atom2bin);
}
/* ---------------------------------------------------------------------- */
@ -87,12 +89,15 @@ void NBin::bin_atoms_setup(int nall)
memory->create(binhead,maxbin,"neigh:binhead");
}
// bins = per-atom vector
// bins and atom2bin = per-atom vectors
// for both local and ghost atoms
if (nall > maxatom) {
maxatom = nall;
memory->destroy(bins);
memory->create(bins,maxatom,"neigh:bins");
memory->destroy(atom2bin);
memory->create(atom2bin,maxatom,"neigh:atom2bin");
}
}
@ -148,6 +153,6 @@ bigint NBin::memory_usage()
{
bigint bytes = 0;
bytes += maxbin*sizeof(int);
bytes += maxatom*sizeof(int);
bytes += 2*maxatom*sizeof(int);
return bytes;
}

View File

@ -31,10 +31,11 @@ class NBin : protected Pointers {
double binsizex,binsizey,binsizez; // bin sizes and inverse sizes
double bininvx,bininvy,bininvz;
int *binhead; // index of first atom in each bin
int *bins; // index of next atom in same bin
int *binhead; // index of first atom in each bin
int *bins; // index of next atom in same bin
int *atom2bin; // bin assignment for each atom (local+ghost)
double cutoff_custom; // cutoff set by requestor
double cutoff_custom; // cutoff set by requestor
NBin(class LAMMPS *);
~NBin();

View File

@ -211,12 +211,14 @@ void NBinStandard::bin_atoms()
for (i = nall-1; i >= nlocal; i--) {
if (mask[i] & bitmask) {
ibin = coord2bin(x[i]);
atom2bin[i] = ibin;
bins[i] = binhead[ibin];
binhead[ibin] = i;
}
}
for (i = atom->nfirst-1; i >= 0; i--) {
ibin = coord2bin(x[i]);
atom2bin[i] = ibin;
bins[i] = binhead[ibin];
binhead[ibin] = i;
}
@ -224,6 +226,7 @@ void NBinStandard::bin_atoms()
} else {
for (i = nall-1; i >= 0; i--) {
ibin = coord2bin(x[i]);
atom2bin[i] = ibin;
bins[i] = binhead[ibin];
binhead[ibin] = i;
}

View File

@ -50,6 +50,7 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp)
// default is no Intel-specific neighbor list build
// default is no Kokkos neighbor list build
// default is no Shardlow Splitting Algorithm (SSA) neighbor list build
// default is no list-specific cutoff
// default is no storage of auxiliary floating point values
occasional = 0;

View File

@ -640,6 +640,24 @@ int Neighbor::init_pair()
delete [] neigh_stencil;
delete [] neigh_pair;
// error check on requests
// do not allow occasional, ghost, bin list
// b/c it still uses variant of coord2bin() in NPair() method
// instead of atom2bin, this could cause error b/c stoms have
// moved out of proc domain by time occasional list is built
// solution would be to use a different NBin variant
// that used Npair::coord2bin(x,ix,iy,iz) (then delete it from NPair)
// and stored the ix,iy,iz values for all atoms (including ghosts)
// at time of binning when neighbor lists are rebuilt,
// similar to what vanilla Nbin::coord2atom() does now in atom2bin
if (style == BIN) {
for (i = 0; i < nrequest; i++)
if (requests[i]->occasional && requests[i]->ghost)
error->all(FLERR,"Cannot request an occasional binned neighbor list "
"with ghost info");
}
// morph requests in various ways
// purpose is to avoid duplicate or inefficient builds
// may add new requests if a needed request to derive from does not exist
@ -1669,7 +1687,6 @@ int Neighbor::choose_stencil(NeighRequest *rq)
else if (rq->newton == 1) newtflag = 1;
else if (rq->newton == 2) newtflag = 0;
//printf("STENCIL RQ FLAGS: hff %d %d n %d g %d s %d newtflag %d\n",
// rq->half,rq->full,rq->newton,rq->ghost,rq->ssa,
// newtflag);
@ -2087,7 +2104,7 @@ void Neighbor::build(int topoflag)
}
// bin atoms for all NBin instances
// not just NBin associated with perpetual lists
// not just NBin associated with perpetual lists, also occasional lists
// b/c cannot wait to bin occasional lists in build_one() call
// if bin then, atoms may have moved outside of proc domain & bin extent,
// leading to errors or even a crash
@ -2193,6 +2210,7 @@ void Neighbor::build_one(class NeighList *mylist, int preflag)
// build the list
if (!mylist->copy) mylist->grow(atom->nlocal,atom->nlocal+atom->nghost);
np->build_setup();
np->build(mylist);
}

View File

@ -128,6 +128,7 @@ void NPair::copy_bin_info()
bininvy = nb->bininvy;
bininvz = nb->bininvz;
atom2bin = nb->atom2bin;
bins = nb->bins;
binhead = nb->binhead;
}
@ -198,53 +199,8 @@ int NPair::exclusion(int i, int j, int itype, int jtype,
}
/* ----------------------------------------------------------------------
convert atom coords into local bin #
for orthogonal, only ghost atoms will have coord >= bboxhi or coord < bboxlo
take special care to insure ghosts are in correct bins even w/ roundoff
hi ghost atoms = nbin,nbin+1,etc
owned atoms = 0 to nbin-1
lo ghost atoms = -1,-2,etc
this is necessary so that both procs on either side of PBC
treat a pair of atoms straddling the PBC in a consistent way
for triclinic, doesn't matter since stencil & neigh list built differently
------------------------------------------------------------------------- */
int NPair::coord2bin(double *x)
{
int ix,iy,iz;
if (!ISFINITE(x[0]) || !ISFINITE(x[1]) || !ISFINITE(x[2]))
error->one(FLERR,"Non-numeric positions - simulation unstable");
if (x[0] >= bboxhi[0])
ix = static_cast<int> ((x[0]-bboxhi[0])*bininvx) + nbinx;
else if (x[0] >= bboxlo[0]) {
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx);
ix = MIN(ix,nbinx-1);
} else
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx) - 1;
if (x[1] >= bboxhi[1])
iy = static_cast<int> ((x[1]-bboxhi[1])*bininvy) + nbiny;
else if (x[1] >= bboxlo[1]) {
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy);
iy = MIN(iy,nbiny-1);
} else
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy) - 1;
if (x[2] >= bboxhi[2])
iz = static_cast<int> ((x[2]-bboxhi[2])*bininvz) + nbinz;
else if (x[2] >= bboxlo[2]) {
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz);
iz = MIN(iz,nbinz-1);
} else
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz) - 1;
return (iz-mbinzlo)*mbiny*mbinx + (iy-mbinylo)*mbinx + (ix-mbinxlo);
}
/* ----------------------------------------------------------------------
same as coord2bin, but also return ix,iy,iz offsets in each dim
same as coord2bin in Nbin, but also return ix,iy,iz offsets in each dim
used by some of the ghost neighbor lists
------------------------------------------------------------------------- */
int NPair::coord2bin(double *x, int &ix, int &iy, int &iz)

View File

@ -77,7 +77,7 @@ class NPair : protected Pointers {
int mbinx,mbiny,mbinz;
int mbinxlo,mbinylo,mbinzlo;
double bininvx,bininvy,bininvz;
int *bins;
int *atom2bin,*bins;
int *binhead;
// data from NStencil class

View File

@ -80,7 +80,7 @@ void NPairFullBin::build(NeighList *list)
// loop over all atoms in surrounding bins in stencil including self
// skip i = j
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {

View File

@ -64,7 +64,7 @@ void NPairFullBinAtomonly::build(NeighList *list)
// loop over all atoms in surrounding bins in stencil including self
// skip i = j
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {

View File

@ -87,7 +87,7 @@ void NPairFullBinGhost::build(NeighList *list)
// no molecular test when i = ghost atom
if (i < nlocal) {
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (i == j) continue;

View File

@ -83,7 +83,7 @@ void NPairFullMulti::build(NeighList *list)
// skip if i,j neighbor cutoff is less than bin distance
// skip i = j
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
s = stencil_multi[itype];
distsq = distsq_multi[itype];
cutsq = cutneighsq[itype];

View File

@ -90,7 +90,8 @@ void NPairHalfBinAtomonlyNewton::build(NeighList *list)
// loop over all atoms in other bins in stencil, store every pair
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
jtype = type[j];

View File

@ -84,7 +84,7 @@ void NPairHalfBinNewtoff::build(NeighList *list)
// stores own/own pairs only once
// stores own/ghost pairs on both procs
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {

View File

@ -92,7 +92,7 @@ void NPairHalfBinNewtoffGhost::build(NeighList *list)
// no molecular test when i = ghost atom
if (i < nlocal) {
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {

View File

@ -119,7 +119,7 @@ void NPairHalfBinNewton::build(NeighList *list)
// loop over all atoms in other bins in stencil, store every pair
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
jtype = type[j];

View File

@ -84,7 +84,7 @@ void NPairHalfBinNewtonTri::build(NeighList *list)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (x[j][2] < ztmp) continue;

View File

@ -87,7 +87,7 @@ void NPairHalfMultiNewtoff::build(NeighList *list)
// stores own/own pairs only once
// stores own/ghost pairs on both procs
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
s = stencil_multi[itype];
distsq = distsq_multi[itype];
cutsq = cutneighsq[itype];

View File

@ -121,7 +121,7 @@ void NPairHalfMultiNewton::build(NeighList *list)
// loop over all atoms in other bins in stencil, store every pair
// skip if i,j neighbor cutoff is less than bin distance
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
s = stencil_multi[itype];
distsq = distsq_multi[itype];
cutsq = cutneighsq[itype];

View File

@ -88,7 +88,7 @@ void NPairHalfMultiNewtonTri::build(NeighList *list)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
s = stencil_multi[itype];
distsq = distsq_multi[itype];
cutsq = cutneighsq[itype];

View File

@ -101,7 +101,7 @@ void NPairHalfRespaBinNewtoff::build(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];

View File

@ -160,7 +160,7 @@ void NPairHalfRespaBinNewton::build(NeighList *list)
// loop over all atoms in other bins in stencil, store every pair
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
jtype = type[j];

View File

@ -113,7 +113,7 @@ void NPairHalfRespaBinNewtonTri::build(NeighList *list)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (x[j][2] < ztmp) continue;

View File

@ -105,7 +105,7 @@ void NPairHalfSizeBinNewtoff::build(NeighList *list)
ytmp = x[i][1];
ztmp = x[i][2];
radi = radius[i];
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
// loop over all atoms in surrounding bins in stencil including self
// only store pair if i < j

View File

@ -156,7 +156,7 @@ void NPairHalfSizeBinNewton::build(NeighList *list)
// loop over all atoms in other bins in stencil, store every pair
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (exclude && exclusion(i,j,type[i],type[j],mask,molecule)) continue;

View File

@ -112,7 +112,7 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
ibin = coord2bin(x[i]);
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (x[j][2] < ztmp) continue;