diff --git a/src/neigh_full.cpp b/src/neigh_full.cpp index 98812d077f..970e3c5725 100644 --- a/src/neigh_full.cpp +++ b/src/neigh_full.cpp @@ -163,6 +163,7 @@ void Neighbor::full_nsq_ghost(NeighList *list) // loop over all atoms, owned and ghost // skip i = j + // no molecular test when i = ghost atom if (i < nlocal) { for (j = 0; j < nall; j++) { @@ -195,8 +196,7 @@ void Neighbor::full_nsq_ghost(NeighList *list) delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; - if (rsq <= cutneighghostsq[itype][jtype]) - neighptr[n++] = j; + if (rsq <= cutneighghostsq[itype][jtype]) neighptr[n++] = j; } } @@ -371,6 +371,7 @@ void Neighbor::full_bin_ghost(NeighList *list) // loop over all atoms in surrounding bins in stencil including self // when i is a ghost atom, must check if stencil bin is out of bounds // skip i = j + // no molecular test when i = ghost atom if (i < nlocal) { ibin = coord2bin(x[i]); @@ -418,8 +419,7 @@ void Neighbor::full_bin_ghost(NeighList *list) delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; - if (rsq <= cutneighghostsq[itype][jtype]) - neighptr[n++] = j; + if (rsq <= cutneighghostsq[itype][jtype]) neighptr[n++] = j; } } } diff --git a/src/neigh_half_bin.cpp b/src/neigh_half_bin.cpp index dfe8456aea..b74ab51189 100644 --- a/src/neigh_half_bin.cpp +++ b/src/neigh_half_bin.cpp @@ -187,6 +187,7 @@ void Neighbor::half_bin_no_newton_ghost(NeighList *list) // stores own/own pairs only once // stores own/ghost pairs with owned atom only, on both procs // stores ghost/ghost pairs only once + // no molecular test when i = ghost atom if (i < nlocal) { ibin = coord2bin(x[i]); @@ -235,8 +236,7 @@ void Neighbor::half_bin_no_newton_ghost(NeighList *list) delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; - if (rsq <= cutneighghostsq[itype][jtype]) - neighptr[n++] = j; + if (rsq <= cutneighghostsq[itype][jtype]) neighptr[n++] = j; } } } diff --git a/src/neigh_half_nsq.cpp b/src/neigh_half_nsq.cpp index 98a8893c64..d18555fc6c 100644 --- a/src/neigh_half_nsq.cpp +++ b/src/neigh_half_nsq.cpp @@ -32,8 +32,6 @@ void Neighbor::half_nsq_no_newton(NeighList *list) double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; - // loop over each atom, storing neighbors - int **special = atom->special; int **nspecial = atom->nspecial; int *tag = atom->tag; @@ -59,6 +57,8 @@ void Neighbor::half_nsq_no_newton(NeighList *list) int npage = 0; int npnt = 0; + // loop over owned atoms, storing neighbors + for (i = 0; i < nlocal; i++) { if (pgsize - npnt < oneatom) { @@ -76,6 +76,7 @@ void Neighbor::half_nsq_no_newton(NeighList *list) ztmp = x[i][2]; // loop over remaining atoms, owned and ghost + // only store pair if i < j for (j = i+1; j < nall; j++) { if (includegroup && !(mask[j] & bitmask)) continue; @@ -109,6 +110,119 @@ void Neighbor::half_nsq_no_newton(NeighList *list) list->inum = inum; } +/* ---------------------------------------------------------------------- + N^2 / 2 search for neighbor pairs with partial Newton's 3rd law + include neighbors of ghost atoms, but no "special neighbors" for ghosts + pair stored once if i,j are both owned and i < j + pair stored by me if i owned and j ghost (also stored by proc owning j) + pair stored once if i,j are both ghost and i < j +------------------------------------------------------------------------- */ + +void Neighbor::half_nsq_no_newton_ghost(NeighList *list) +{ + int i,j,n,itype,jtype,which,bitmask; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr; + + int **special = atom->special; + int **nspecial = atom->nspecial; + int *tag = atom->tag; + + double **x = atom->x; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + int molecular = atom->molecular; + if (includegroup) { + nlocal = atom->nfirst; + bitmask = group->bitmask[includegroup]; + } + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int **pages = list->pages; + + int inum = 0; + int npage = 0; + int npnt = 0; + + // loop over owned & ghost atoms, storing neighbors + + for (i = 0; i < nall; i++) { + + if (pgsize - npnt < oneatom) { + npnt = 0; + npage++; + if (npage == list->maxpage) pages = list->add_pages(); + } + + neighptr = &pages[npage][npnt]; + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over remaining atoms, owned and ghost + // only store pair if i < j + // stores own/own pairs only once + // stores own/ghost pairs with owned atom only, on both procs + // stores ghost/ghost pairs only once + // no molecular test when i = ghost atom + + if (i < nlocal) { + for (j = i+1; j < nall; j++) { + if (includegroup && !(mask[j] & bitmask)) continue; + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq <= cutneighsq[itype][jtype]) { + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + + } else { + for (j = i+1; j < nall; j++) { + if (includegroup && !(mask[j] & bitmask)) continue; + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq <= cutneighsq[itype][jtype]) neighptr[n++] = j; + } + } + + ilist[inum++] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + } + + list->inum = atom->nlocal; + list->gnum = inum - atom->nlocal; +} + /* ---------------------------------------------------------------------- N^2 / 2 search for neighbor pairs with full Newton's 3rd law every pair stored exactly once by some processor diff --git a/src/neigh_stencil.cpp b/src/neigh_stencil.cpp index 53c71286da..0474e3348b 100644 --- a/src/neigh_stencil.cpp +++ b/src/neigh_stencil.cpp @@ -61,6 +61,28 @@ void Neighbor::stencil_half_bin_2d_no_newton(NeighList *list, /* ---------------------------------------------------------------------- */ +void Neighbor::stencil_half_ghost_bin_2d_no_newton(NeighList *list, + int sx, int sy, int sz) +{ + int i,j,k; + int *stencil = list->stencil; + int **stencilxyz = list->stencilxyz; + int nstencil = 0; + + for (j = -sy; j <= sy; j++) + for (i = -sx; i <= sx; i++) + if (bin_distance(i,j,0) < cutneighmaxsq) { + stencilxyz[nstencil][0] = i; + stencilxyz[nstencil][1] = j; + stencilxyz[nstencil][2] = 0; + stencil[nstencil++] = j*mbinx + i; + } + + list->nstencil = nstencil; +} + +/* ---------------------------------------------------------------------- */ + void Neighbor::stencil_half_bin_3d_no_newton(NeighList *list, int sx, int sy, int sz) { diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 52c249196e..386d962694 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -826,24 +826,48 @@ void Neighbor::choose_build(int index, NeighRequest *rq) } else if (rq->half) { if (style == NSQ) { if (rq->newton == 0) { - if (newton_pair == 0) pb = &Neighbor::half_nsq_no_newton; - else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton; + if (newton_pair == 0) { + if (rq->ghost == 0) pb = &Neighbor::half_nsq_no_newton; + else if (includegroup) + error->all(FLERR,"Neighbor include group not allowed " + "with ghost neighbors"); + else pb = &Neighbor::half_nsq_no_newton_ghost; + } else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton; } else if (rq->newton == 1) { pb = &Neighbor::half_nsq_newton; - } else if (rq->newton == 2) pb = &Neighbor::half_nsq_no_newton; + } else if (rq->newton == 2) { + if (rq->ghost == 0) pb = &Neighbor::half_nsq_no_newton; + else if (includegroup) + error->all(FLERR,"Neighbor include group not allowed " + "with ghost neighbors"); + else pb = &Neighbor::half_nsq_no_newton_ghost; + } } else if (style == BIN) { if (rq->newton == 0) { - if (newton_pair == 0) pb = &Neighbor::half_bin_no_newton; - else if (triclinic == 0) pb = &Neighbor::half_bin_newton; - else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri; + if (newton_pair == 0) { + if (rq->ghost == 0) pb = &Neighbor::half_bin_no_newton; + else if (includegroup) + error->all(FLERR,"Neighbor include group not allowed " + "with ghost neighbors"); + else pb = &Neighbor::half_bin_no_newton_ghost; + } else if (triclinic == 0) { + pb = &Neighbor::half_bin_newton; + } else if (triclinic == 1) + pb = &Neighbor::half_bin_newton_tri; } else if (rq->newton == 1) { if (triclinic == 0) pb = &Neighbor::half_bin_newton; else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri; } else if (rq->newton == 2) { - if (rq->ghost) pb = &Neighbor::half_bin_no_newton_ghost; - else pb = &Neighbor::half_bin_no_newton; + if (rq->ghost == 0) pb = &Neighbor::half_bin_no_newton; + else if (includegroup) + error->all(FLERR,"Neighbor include group not allowed " + "with ghost neighbors"); + else pb = &Neighbor::half_bin_no_newton_ghost; } } else if (style == MULTI) { + if (rq->ghost == 1) + error->all(FLERR, + "Neighbor multi not yet enabled for ghost neighbors"); if (rq->newton == 0) { if (newton_pair == 0) pb = &Neighbor::half_multi_no_newton; else if (triclinic == 0) pb = &Neighbor::half_multi_newton; @@ -860,17 +884,18 @@ void Neighbor::choose_build(int index, NeighRequest *rq) else if (includegroup) error->all(FLERR, "Neighbor include group not allowed with ghost neighbors"); - else if (rq->ghost == 1) pb = &Neighbor::full_nsq_ghost; + else pb = &Neighbor::full_nsq_ghost; } else if (style == BIN) { if (rq->ghost == 0) pb = &Neighbor::full_bin; else if (includegroup) error->all(FLERR, "Neighbor include group not allowed with ghost neighbors"); - else if (rq->ghost == 1) pb = &Neighbor::full_bin_ghost; + else pb = &Neighbor::full_bin_ghost; } else if (style == MULTI) { - if (rq->ghost == 0) pb = &Neighbor::full_multi; - else error->all(FLERR, - "Neighbor multi not yet enabled for ghost neighbors"); + if (rq->ghost == 1) + error->all(FLERR, + "Neighbor multi not yet enabled for ghost neighbors"); + pb = &Neighbor::full_multi; } } else if (rq->gran) { @@ -895,6 +920,9 @@ void Neighbor::choose_build(int index, NeighRequest *rq) } else if (style == MULTI) error->all(FLERR,"Neighbor multi not yet enabled for rRESPA"); } + + // OMP versions of build methods + } else { if (rq->copy) pb = &Neighbor::copy_from; @@ -912,23 +940,52 @@ void Neighbor::choose_build(int index, NeighRequest *rq) } else if (rq->half) { if (style == NSQ) { if (rq->newton == 0) { - if (newton_pair == 0) pb = &Neighbor::half_nsq_no_newton_omp; - else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton_omp; + if (newton_pair == 0) { + if (rq->ghost == 0) pb = &Neighbor::half_nsq_no_newton_omp; + else if (includegroup) + error->all(FLERR,"Neighbor include group not allowed " + "with ghost neighbors"); + // NOTE: missing OMP version of this one + else pb = &Neighbor::half_nsq_no_newton_ghost; + } else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton_omp; } else if (rq->newton == 1) { pb = &Neighbor::half_nsq_newton_omp; } else if (rq->newton == 2) { - pb = &Neighbor::half_nsq_no_newton_omp; + if (rq->ghost == 0) pb = &Neighbor::half_nsq_no_newton_omp; + else if (includegroup) + error->all(FLERR,"Neighbor include group not allowed " + "with ghost neighbors"); + // NOTE: missing OMP version of this one + else pb = &Neighbor::half_nsq_no_newton_ghost; } } else if (style == BIN) { if (rq->newton == 0) { - if (newton_pair == 0) pb = &Neighbor::half_bin_no_newton_omp; - else if (triclinic == 0) pb = &Neighbor::half_bin_newton_omp; - else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri_omp; + if (newton_pair == 0) { + if (rq->ghost == 0) pb = &Neighbor::half_bin_no_newton_omp; + else if (includegroup) + error->all(FLERR,"Neighbor include group not allowed " + "with ghost neighbors"); + // NOTE: missing OMP version of this one + else pb = &Neighbor::half_bin_no_newton_ghost; + } else if (triclinic == 0) { + pb = &Neighbor::half_bin_newton_omp; + } else if (triclinic == 1) + pb = &Neighbor::half_bin_newton_tri_omp; } else if (rq->newton == 1) { if (triclinic == 0) pb = &Neighbor::half_bin_newton_omp; else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri_omp; - } else if (rq->newton == 2) pb = &Neighbor::half_bin_no_newton_omp; + } else if (rq->newton == 2) { + if (rq->ghost == 0) pb = &Neighbor::half_bin_no_newton_omp; + else if (includegroup) + error->all(FLERR,"Neighbor include group not allowed " + "with ghost neighbors"); + // NOTE: missing OMP version of this one + else pb = &Neighbor::half_bin_no_newton_ghost; + } } else if (style == MULTI) { + if (rq->ghost == 1) + error->all(FLERR, + "Neighbor multi not yet enabled for ghost neighbors"); if (rq->newton == 0) { if (newton_pair == 0) pb = &Neighbor::half_multi_no_newton_omp; else if (triclinic == 0) pb = &Neighbor::half_multi_newton_omp; @@ -945,17 +1002,18 @@ void Neighbor::choose_build(int index, NeighRequest *rq) else if (includegroup) error->all(FLERR, "Neighbor include group not allowed with ghost neighbors"); - else if (rq->ghost == 1) pb = &Neighbor::full_nsq_ghost_omp; + else pb = &Neighbor::full_nsq_ghost_omp; } else if (style == BIN) { if (rq->ghost == 0) pb = &Neighbor::full_bin_omp; else if (includegroup) error->all(FLERR, "Neighbor include group not allowed with ghost neighbors"); - else if (rq->ghost == 1) pb = &Neighbor::full_bin_ghost_omp; + else pb = &Neighbor::full_bin_ghost_omp; } else if (style == MULTI) { - if (rq->ghost == 0) pb = &Neighbor::full_multi_omp; - else error->all(FLERR, - "Neighbor multi not yet enabled for ghost neighbors"); + if (rq->ghost == 1) + error->all(FLERR, + "Neighbor multi not yet enabled for ghost neighbors"); + pb = &Neighbor::full_multi_omp; } } else if (rq->gran) { @@ -1004,10 +1062,13 @@ void Neighbor::choose_stencil(int index, NeighRequest *rq) if (style == BIN) { if (rq->newton == 0) { if (newton_pair == 0) { - if (dimension == 2) - sc = &Neighbor::stencil_half_bin_2d_no_newton; - else if (dimension == 3) - sc = &Neighbor::stencil_half_bin_3d_no_newton; + if (dimension == 2) { + if (rq->ghost) sc = &Neighbor::stencil_half_ghost_bin_2d_no_newton; + else sc = &Neighbor::stencil_half_bin_2d_no_newton; + } else if (dimension == 3) { + if (rq->ghost) sc = &Neighbor::stencil_half_ghost_bin_3d_no_newton; + else sc = &Neighbor::stencil_half_bin_3d_no_newton; + } } else if (triclinic == 0) { if (dimension == 2) sc = &Neighbor::stencil_half_bin_2d_newton; @@ -1033,7 +1094,8 @@ void Neighbor::choose_stencil(int index, NeighRequest *rq) } } else if (rq->newton == 2) { if (dimension == 2) - sc = &Neighbor::stencil_half_bin_2d_no_newton; + if (rq->ghost) sc = &Neighbor::stencil_half_ghost_bin_2d_no_newton; + else sc = &Neighbor::stencil_half_bin_2d_no_newton; else if (dimension == 3) { if (rq->ghost) sc = &Neighbor::stencil_half_ghost_bin_3d_no_newton; else sc = &Neighbor::stencil_half_bin_3d_no_newton; diff --git a/src/neighbor.h b/src/neighbor.h index b93e2fab9a..96b8ff9a9c 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -177,6 +177,7 @@ class Neighbor : protected Pointers { PairPtr *pair_build; void half_nsq_no_newton(class NeighList *); + void half_nsq_no_newton_ghost(class NeighList *); void half_nsq_newton(class NeighList *); void half_bin_no_newton(class NeighList *); @@ -223,6 +224,7 @@ class Neighbor : protected Pointers { StencilPtr *stencil_create; void stencil_half_bin_2d_no_newton(class NeighList *, int, int, int); + void stencil_half_ghost_bin_2d_no_newton(class NeighList *, int, int, int); void stencil_half_bin_3d_no_newton(class NeighList *, int, int, int); void stencil_half_ghost_bin_3d_no_newton(class NeighList *, int, int, int); void stencil_half_bin_2d_newton(class NeighList *, int, int, int);