From ab6e3568083a2e29c19e6f9bd0a3fe4771b02eaf Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 2 Dec 2011 18:25:56 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7290 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-OMP/neigh_derive_omp.cpp | 184 +++++ src/USER-OMP/neigh_full_omp.cpp | 564 +++++++++++++++ src/USER-OMP/neigh_gran_omp.cpp | 657 +++++++++++++++++ src/USER-OMP/neigh_half_bin_omp.cpp | 364 ++++++++++ src/USER-OMP/neigh_half_multi_omp.cpp | 392 ++++++++++ src/USER-OMP/neigh_half_nsq_omp.cpp | 222 ++++++ src/USER-OMP/neigh_respa_omp.cpp | 985 ++++++++++++++++++++++++++ src/USER-OMP/neighbor_omp.h | 60 ++ 8 files changed, 3428 insertions(+) create mode 100644 src/USER-OMP/neigh_derive_omp.cpp create mode 100644 src/USER-OMP/neigh_full_omp.cpp create mode 100644 src/USER-OMP/neigh_gran_omp.cpp create mode 100644 src/USER-OMP/neigh_half_bin_omp.cpp create mode 100644 src/USER-OMP/neigh_half_multi_omp.cpp create mode 100644 src/USER-OMP/neigh_half_nsq_omp.cpp create mode 100644 src/USER-OMP/neigh_respa_omp.cpp create mode 100644 src/USER-OMP/neighbor_omp.h diff --git a/src/USER-OMP/neigh_derive_omp.cpp b/src/USER-OMP/neigh_derive_omp.cpp new file mode 100644 index 0000000000..3b8fd79a20 --- /dev/null +++ b/src/USER-OMP/neigh_derive_omp.cpp @@ -0,0 +1,184 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "neighbor.h" +#include "neighbor_omp.h" +#include "neigh_list.h" +#include "atom.h" +#include "comm.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- + build half list from full list + pair stored once if i,j are both owned and i < j + pair stored by me if j is ghost (also stored by proc owning j) + works if full list is a skip list +------------------------------------------------------------------------- */ + +void Neighbor::half_from_full_no_newton_omp(NeighList *list) +{ + const int inum_full = list->listfull->inum; + + NEIGH_OMP_INIT; + +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(inum_full); + + int i,j,ii,jj,n,jnum,joriginal; + int *neighptr,*jlist; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int *ilist_full = list->listfull->ilist; + int *numneigh_full = list->listfull->numneigh; + int **firstneigh_full = list->listfull->firstneigh; + + // each thread works on its own page + int npage = tid; + int npnt = 0; + + // loop over atoms in full list + + for (ii = ifrom; ii < ito; ii++) { + + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + // only one thread at a time may check whether we + // need new neighbor list pages and then add to them. +#if defined(_OPENMP) +#pragma omp critical +#endif + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + neighptr = &(list->pages[npage][npnt]); + n = 0; + + // loop over parent full list + + i = ilist_full[ii]; + jlist = firstneigh_full[i]; + jnum = numneigh_full[i]; + + for (jj = 0; jj < jnum; jj++) { + joriginal = jlist[jj]; + j = joriginal & NEIGHMASK; + if (j > i) neighptr[n++] = joriginal; + } + + ilist[ii] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = inum_full; +} + +/* ---------------------------------------------------------------------- + build half list from full list + pair stored once if i,j are both owned and i < j + if j is ghost, only store if j coords are "above and to the right" of i + works if full list is a skip list +------------------------------------------------------------------------- */ + +void Neighbor::half_from_full_newton_omp(NeighList *list) +{ + const int inum_full = list->listfull->inum; + + NEIGH_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(inum_full); + + int i,j,ii,jj,n,jnum,joriginal; + int *neighptr,*jlist; + double xtmp,ytmp,ztmp; + + double **x = atom->x; + int nlocal = atom->nlocal; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int *ilist_full = list->listfull->ilist; + int *numneigh_full = list->listfull->numneigh; + int **firstneigh_full = list->listfull->firstneigh; + + // each thread works on its own page + int npage = tid; + int npnt = 0; + + // loop over parent full list + + for (ii = ifrom; ii < ito; ii++) { + + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + // only one thread at a time may check whether we + // need new neighbor list pages and then add to them. +#if defined(_OPENMP) +#pragma omp critical +#endif + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + neighptr = &(list->pages[npage][npnt]); + n = 0; + + i = ilist_full[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over full neighbor list + + jlist = firstneigh_full[i]; + jnum = numneigh_full[i]; + + for (jj = 0; jj < jnum; jj++) { + joriginal = jlist[jj]; + j = joriginal & NEIGHMASK; + if (j < nlocal) { + if (i > j) continue; + } else { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + } + neighptr[n++] = joriginal; + } + + ilist[ii] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = inum_full; +} + diff --git a/src/USER-OMP/neigh_full_omp.cpp b/src/USER-OMP/neigh_full_omp.cpp new file mode 100644 index 0000000000..5f3ebd4d96 --- /dev/null +++ b/src/USER-OMP/neigh_full_omp.cpp @@ -0,0 +1,564 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "neighbor.h" +#include "neighbor_omp.h" +#include "neigh_list.h" +#include "atom.h" +#include "comm.h" +#include "group.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- + N^2 search for all neighbors + every neighbor pair appears in list of both atoms i and j +------------------------------------------------------------------------- */ + +void Neighbor::full_nsq_omp(NeighList *list) +{ + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0; + + NEIGH_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,n,itype,jtype,which; + 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 nall = atom->nlocal + atom->nghost; + int molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + + int npage = tid; + int npnt = 0; + + // loop over owned atoms, storing neighbors + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + neighptr = &(list->pages[npage][npnt]); + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over all atoms, owned and ghost + // skip i = j + + for (j = 0; j < nall; j++) { + if (includegroup && !(mask[j] & bitmask)) continue; + if (i == j) 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; + list->gnum = 0; +} + +/* ---------------------------------------------------------------------- + N^2 search for all neighbors + include neighbors of ghost atoms (no "special neighbors" for ghosts) + every neighbor pair appears in list of both atoms i and j +------------------------------------------------------------------------- */ + +void Neighbor::full_nsq_ghost_omp(NeighList *list) +{ + const int nlocal = atom->nlocal; + const int nall = nlocal + atom->nghost; + + NEIGH_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(nall); + + int i,j,n,itype,jtype,which; + 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 molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + + int npage = tid; + int npnt = 0; + + // loop over owned & ghost atoms, storing neighbors + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + neighptr = &(list->pages[npage][npnt]); + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over all atoms, owned and ghost + // skip i = j + + if (i < nlocal) { + for (j = 0; j < nall; j++) { + if (i == j) 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + } else { + for (j = 0; j < nall; j++) { + if (i == j) 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 <= cutneighghostsq[itype][jtype]) + neighptr[n++] = j; + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; + list->gnum = nall - nlocal; +} + +/* ---------------------------------------------------------------------- + binned neighbor list construction for all neighbors + every neighbor pair appears in list of both atoms i and j +------------------------------------------------------------------------- */ + +void Neighbor::full_bin_omp(NeighList *list) +{ + // bin owned & ghost atoms + + bin_atoms(); + + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + + NEIGH_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,k,n,itype,jtype,ibin,which; + 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 molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int nstencil = list->nstencil; + int *stencil = list->stencil; + + int npage = tid; + int npnt = 0; + + // loop over owned atoms, storing neighbors + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + neighptr = &(list->pages[npage][npnt]); + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over all atoms in surrounding bins in stencil including self + // skip i = j + + ibin = coord2bin(x[i]); + + for (k = 0; k < nstencil; k++) { + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (i == j) 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; + list->gnum = 0; +} + +/* ---------------------------------------------------------------------- + binned neighbor list construction for all neighbors + include neighbors of ghost atoms (no "special neighbors" for ghosts) + every neighbor pair appears in list of both atoms i and j +------------------------------------------------------------------------- */ + +void Neighbor::full_bin_ghost_omp(NeighList *list) +{ + // bin owned & ghost atoms + + bin_atoms(); + + const int nlocal = atom->nlocal; + const int nall = nlocal + atom->nghost; + + NEIGH_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(nall); + + int i,j,k,n,itype,jtype,ibin,which; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int xbin,ybin,zbin,xbin2,ybin2,zbin2; + 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 molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int nstencil = list->nstencil; + int *stencil = list->stencil; + int **stencilxyz = list->stencilxyz; + + int npage = tid; + int npnt = 0; + + // loop over owned & ghost atoms, storing neighbors + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + neighptr = &(list->pages[npage][npnt]); + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // 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 + + if (i < nlocal) { + ibin = coord2bin(x[i]); + for (k = 0; k < nstencil; k++) { + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (i == j) 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + } + + } else { + ibin = coord2bin(x[i],xbin,ybin,zbin); + for (k = 0; k < nstencil; k++) { + xbin2 = xbin + stencilxyz[k][0]; + ybin2 = ybin + stencilxyz[k][1]; + zbin2 = zbin + stencilxyz[k][2]; + if (xbin2 < 0 || xbin2 >= mbinx || + ybin2 < 0 || ybin2 >= mbiny || + zbin2 < 0 || zbin2 >= mbinz) continue; + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (i == j) 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 <= cutneighghostsq[itype][jtype]) + neighptr[n++] = j; + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; + list->gnum = nall - nlocal; +} + +/* ---------------------------------------------------------------------- + binned neighbor list construction for all neighbors + multi-type stencil is itype dependent and is distance checked + every neighbor pair appears in list of both atoms i and j +------------------------------------------------------------------------- */ + +void Neighbor::full_multi_omp(NeighList *list) +{ + // bin local & ghost atoms + + bin_atoms(); + + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + + NEIGH_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,k,n,itype,jtype,ibin,which,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr,*s; + double *cutsq,*distsq; + + // loop over each atom, storing neighbors + + 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 molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int *nstencil_multi = list->nstencil_multi; + int **stencil_multi = list->stencil_multi; + double **distsq_multi = list->distsq_multi; + + int npage = tid; + int npnt = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + neighptr = &(list->pages[npage][npnt]); + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over all atoms in other bins in stencil, including self + // skip if i,j neighbor cutoff is less than bin distance + // skip i = j + + ibin = coord2bin(x[i]); + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + if (i == j) continue; + + 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; + list->gnum = 0; +} diff --git a/src/USER-OMP/neigh_gran_omp.cpp b/src/USER-OMP/neigh_gran_omp.cpp new file mode 100644 index 0000000000..eca58ddccd --- /dev/null +++ b/src/USER-OMP/neigh_gran_omp.cpp @@ -0,0 +1,657 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "neighbor.h" +#include "neighbor_omp.h" +#include "neigh_list.h" +#include "atom.h" +#include "comm.h" +#include "group.h" +#include "fix_shear_history.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- + granular particles + N^2 / 2 search for neighbor pairs with partial Newton's 3rd law + shear history must be accounted for when a neighbor pair is added + pair added to list if atoms i and j are both owned and i < j + pair added if j is ghost (also stored by proc owning j) +------------------------------------------------------------------------- */ + +void Neighbor::granular_nsq_no_newton_omp(NeighList *list) +{ + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0; + + FixShearHistory * const fix_history = list->fix_history; + NeighList * listgranhistory = list->listgranhistory; + + NEIGH_OMP_INIT; + + if (fix_history) + if (nthreads > listgranhistory->maxpage) + listgranhistory->add_pages(nthreads - listgranhistory->maxpage); + +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list,listgranhistory) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,m,n,nn; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutsq; + int *neighptr,*touchptr; + double *shearptr; + + int *npartner,**partner; + double ***shearpartner; + int **firsttouch; + double **firstshear; + + double **x = atom->x; + double *radius = atom->radius; + int *tag = atom->tag; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int nall = atom->nlocal + atom->nghost; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + + if (fix_history) { + npartner = fix_history->npartner; + partner = fix_history->partner; + shearpartner = fix_history->shearpartner; + firsttouch = listgranhistory->firstneigh; + firstshear = listgranhistory->firstdouble; + } + + int npage = tid; + int npnt = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + { + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) { + list->add_pages(nthreads); + if (fix_history) + listgranhistory->add_pages(nthreads); + } + } + + n = nn = 0; + neighptr = &(list->pages[npage][npnt]); + if (fix_history) { + touchptr = &(listgranhistory->pages[npage][npnt]); + shearptr = &(listgranhistory->dpages[npage][3*npnt]); + } + } + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + + // loop over remaining atoms, owned and ghost + + for (j = i+1; j < nall; j++) { + if (includegroup && !(mask[j] & bitmask)) continue; + if (exclude && exclusion(i,j,type[i],type[j],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; + radsum = radi + radius[j]; + cutsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutsq) { + neighptr[n] = j; + + if (fix_history) { + if (rsq < radsum*radsum) { + for (m = 0; m < npartner[i]; m++) + if (partner[i][m] == tag[j]) break; + if (m < npartner[i]) { + touchptr[n] = 1; + shearptr[nn++] = shearpartner[i][m][0]; + shearptr[nn++] = shearpartner[i][m][1]; + shearptr[nn++] = shearpartner[i][m][2]; + } else { + touchptr[n] = 0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + } + } else { + touchptr[n] = 0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + } + } + + n++; + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + if (fix_history) { + firsttouch[i] = touchptr; + firstshear[i] = shearptr; + } + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; +} + +/* ---------------------------------------------------------------------- + granular particles + N^2 / 2 search for neighbor pairs with full Newton's 3rd law + no shear history is allowed for this option + pair added to list if atoms i and j are both owned and i < j + if j is ghost only me or other proc adds pair + decision based on itag,jtag tests +------------------------------------------------------------------------- */ + +void Neighbor::granular_nsq_newton_omp(NeighList *list) +{ + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0; + + NEIGH_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,n,itag,jtag; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutsq; + int *neighptr; + + double **x = atom->x; + double *radius = atom->radius; + int *tag = atom->tag; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int nall = atom->nlocal + atom->nghost; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + + int npage = tid; + int npnt = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + neighptr = &(list->pages[npage][npnt]); + n = 0; + + itag = tag[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + + // loop over remaining atoms, owned and ghost + + for (j = i+1; j < nall; j++) { + if (includegroup && !(mask[j] & bitmask)) continue; + + if (j >= nlocal) { + 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) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + } + } + + if (exclude && exclusion(i,j,type[i],type[j],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; + radsum = radi + radius[j]; + cutsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutsq) neighptr[n++] = j; + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; +} + +/* ---------------------------------------------------------------------- + granular particles + binned neighbor list construction with partial Newton's 3rd law + shear history must be accounted for when a neighbor pair is added + each owned atom i checks own bin and surrounding bins in non-Newton stencil + pair stored once if i,j are both owned and i < j + pair stored by me if j is ghost (also stored by proc owning j) +------------------------------------------------------------------------- */ + +void Neighbor::granular_bin_no_newton_omp(NeighList *list) +{ + // bin local & ghost atoms + + bin_atoms(); + + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + + FixShearHistory * const fix_history = list->fix_history; + NeighList * listgranhistory = list->listgranhistory; + + NEIGH_OMP_INIT; + + if (fix_history) + if (nthreads > listgranhistory->maxpage) + listgranhistory->add_pages(nthreads - listgranhistory->maxpage); + +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list,listgranhistory) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,k,m,n,nn,ibin; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutsq; + int *neighptr,*touchptr; + double *shearptr; + + int *npartner,**partner; + double ***shearpartner; + int **firsttouch; + double **firstshear; + + // loop over each atom, storing neighbors + + double **x = atom->x; + double *radius = atom->radius; + int *tag = atom->tag; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int nstencil = list->nstencil; + int *stencil = list->stencil; + + if (fix_history) { + npartner = fix_history->npartner; + partner = fix_history->partner; + shearpartner = fix_history->shearpartner; + firsttouch = listgranhistory->firstneigh; + firstshear = listgranhistory->firstdouble; + } + + int npage = tid; + int npnt = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + { + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) { + list->add_pages(nthreads); + if (fix_history) + listgranhistory->add_pages(nthreads); + } + } + + n = nn = 0; + neighptr = &(list->pages[npage][npnt]); + if (fix_history) { + touchptr = &(listgranhistory->pages[npage][npnt]); + shearptr = &(listgranhistory->dpages[npage][3*npnt]); + } + } + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + ibin = coord2bin(x[i]); + + // loop over all atoms in surrounding bins in stencil including self + // only store pair if i < j + // stores own/own pairs only once + // stores own/ghost pairs on both procs + + for (k = 0; k < nstencil; k++) { + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (j <= i) continue; + if (exclude && exclusion(i,j,type[i],type[j],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; + radsum = radi + radius[j]; + cutsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutsq) { + neighptr[n] = j; + + if (fix_history) { + if (rsq < radsum*radsum) { + for (m = 0; m < npartner[i]; m++) + if (partner[i][m] == tag[j]) break; + if (m < npartner[i]) { + touchptr[n] = 1; + shearptr[nn++] = shearpartner[i][m][0]; + shearptr[nn++] = shearpartner[i][m][1]; + shearptr[nn++] = shearpartner[i][m][2]; + } else { + touchptr[n] = 0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + } + } else { + touchptr[n] = 0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + shearptr[nn++] = 0.0; + } + } + + n++; + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + if (fix_history) { + firsttouch[i] = touchptr; + firstshear[i] = shearptr; + } + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; +} + +/* ---------------------------------------------------------------------- + granular particles + binned neighbor list construction with full Newton's 3rd law + no shear history is allowed for this option + each owned atom i checks its own bin and other bins in Newton stencil + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void Neighbor::granular_bin_newton_omp(NeighList *list) +{ + // bin local & ghost atoms + + bin_atoms(); + + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + + NEIGH_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,k,n,ibin; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutsq; + int *neighptr; + + // loop over each atom, storing neighbors + + double **x = atom->x; + double *radius = atom->radius; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int nstencil = list->nstencil; + int *stencil = list->stencil; + + int npage = tid; + int npnt = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + n = 0; + neighptr = &(list->pages[npage][npnt]); + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + + // loop over rest of atoms in i's bin, ghosts are at end of linked list + // if j is owned atom, store it, since j is beyond i in linked list + // if j is ghost, only store if j coords are "above and to the right" of i + + for (j = bins[i]; j >= 0; j = bins[j]) { + if (j >= nlocal) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + } + + if (exclude && exclusion(i,j,type[i],type[j],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; + radsum = radi + radius[j]; + cutsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutsq) neighptr[n++] = j; + } + + // loop over all atoms in other bins in stencil, store every pair + + ibin = coord2bin(x[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; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutsq) neighptr[n++] = j; + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; +} + +/* ---------------------------------------------------------------------- + granular particles + binned neighbor list construction with Newton's 3rd law for triclinic + no shear history is allowed for this option + each owned atom i checks its own bin and other bins in triclinic stencil + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void Neighbor::granular_bin_newton_tri_omp(NeighList *list) +{ + // bin local & ghost atoms + + bin_atoms(); + + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + + NEIGH_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,k,n,ibin; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutsq; + int *neighptr; + + // loop over each atom, storing neighbors + + double **x = atom->x; + double *radius = atom->radius; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int nstencil = list->nstencil; + int *stencil = list->stencil; + + int npage = tid; + int npnt = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + n = 0; + neighptr = &(list->pages[npage][npnt]); + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + + // loop over all atoms in bins in stencil + // pairs for atoms j "below" i are excluded + // below = lower z or (equal z and lower y) or (equal zy and lower x) + // (equal zyx and j <= i) + // latter excludes self-self interaction but allows superposed atoms + + ibin = coord2bin(x[i]); + for (k = 0; k < nstencil; k++) { + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp) { + if (x[j][0] < xtmp) continue; + if (x[j][0] == xtmp && j <= i) continue; + } + } + + if (exclude && exclusion(i,j,type[i],type[j],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; + radsum = radi + radius[j]; + cutsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutsq) neighptr[n++] = j; + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; +} diff --git a/src/USER-OMP/neigh_half_bin_omp.cpp b/src/USER-OMP/neigh_half_bin_omp.cpp new file mode 100644 index 0000000000..3df11d38ac --- /dev/null +++ b/src/USER-OMP/neigh_half_bin_omp.cpp @@ -0,0 +1,364 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "neighbor.h" +#include "neighbor_omp.h" +#include "neigh_list.h" +#include "atom.h" +#include "comm.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- + binned neighbor list construction with partial Newton's 3rd law + each owned atom i checks own bin and other bins in stencil + pair stored once if i,j are both owned and i < j + pair stored by me if j is ghost (also stored by proc owning j) +------------------------------------------------------------------------- */ + +void Neighbor::half_bin_no_newton_omp(NeighList *list) +{ + // bin local & ghost atoms + + bin_atoms(); + + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + + NEIGH_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,k,n,itype,jtype,ibin,which; + 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; + + double **x = atom->x; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int nstencil = list->nstencil; + int *stencil = list->stencil; + + // each thread works on its own page + int npage = tid; + int npnt = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + neighptr = &(list->pages[npage][npnt]); + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over all atoms in other bins in stencil including self + // only store pair if i < j + // stores own/own pairs only once + // stores own/ghost pairs on both procs + + ibin = coord2bin(x[i]); + + for (k = 0; k < nstencil; k++) { + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (j <= i) 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; +} + +/* ---------------------------------------------------------------------- + binned neighbor list construction with full Newton's 3rd law + each owned atom i checks its own bin and other bins in Newton stencil + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void Neighbor::half_bin_newton_omp(NeighList *list) +{ + // bin local & ghost atoms + + bin_atoms(); + + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + + NEIGH_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,k,n,itype,jtype,ibin,which; + 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; + + double **x = atom->x; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int nstencil = list->nstencil; + int *stencil = list->stencil; + + // each thread works on its own page + int npage = tid; + int npnt = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + neighptr = &(list->pages[npage][npnt]); + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over rest of atoms in i's bin, ghosts are at end of linked list + // if j is owned atom, store it, since j is beyond i in linked list + // if j is ghost, only store if j coords are "above and to the right" of i + + for (j = bins[i]; j >= 0; j = bins[j]) { + if (j >= nlocal) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + + // loop over all atoms in other bins in stencil, store every pair + + ibin = coord2bin(x[i]); + for (k = 0; k < nstencil; k++) { + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; +} + +/* ---------------------------------------------------------------------- + binned neighbor list construction with Newton's 3rd law for triclinic + each owned atom i checks its own bin and other bins in triclinic stencil + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void Neighbor::half_bin_newton_tri_omp(NeighList *list) +{ + // bin local & ghost atoms + + bin_atoms(); + + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + + NEIGH_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,k,n,itype,jtype,ibin,which; + 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; + + double **x = atom->x; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int nstencil = list->nstencil; + int *stencil = list->stencil; + + // each thread works on its own page + int npage = tid; + int npnt = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + neighptr = &(list->pages[npage][npnt]); + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over all atoms in bins in stencil + // pairs for atoms j "below" i are excluded + // below = lower z or (equal z and lower y) or (equal zy and lower x) + // (equal zyx and j <= i) + // latter excludes self-self interaction but allows superposed atoms + + ibin = coord2bin(x[i]); + for (k = 0; k < nstencil; k++) { + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp) { + if (x[j][0] < xtmp) continue; + if (x[j][0] == xtmp && j <= i) 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; +} diff --git a/src/USER-OMP/neigh_half_multi_omp.cpp b/src/USER-OMP/neigh_half_multi_omp.cpp new file mode 100644 index 0000000000..dcc2661310 --- /dev/null +++ b/src/USER-OMP/neigh_half_multi_omp.cpp @@ -0,0 +1,392 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "neighbor.h" +#include "neighbor_omp.h" +#include "neigh_list.h" +#include "atom.h" +#include "comm.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- + binned neighbor list construction with partial Newton's 3rd law + each owned atom i checks own bin and other bins in stencil + multi-type stencil is itype dependent and is distance checked + pair stored once if i,j are both owned and i < j + pair stored by me if j is ghost (also stored by proc owning j) +------------------------------------------------------------------------- */ + +void Neighbor::half_multi_no_newton_omp(NeighList *list) +{ + // bin local & ghost atoms + + bin_atoms(); + + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + + NEIGH_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,k,n,itype,jtype,ibin,which,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr,*s; + double *cutsq,*distsq; + + // loop over each atom, storing neighbors + + 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 molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int *nstencil_multi = list->nstencil_multi; + int **stencil_multi = list->stencil_multi; + double **distsq_multi = list->distsq_multi; + + // each thread works on its own page + int npage = tid; + int npnt = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + neighptr = &(list->pages[npage][npnt]); + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over all atoms in other bins in stencil including self + // only store pair if i < j + // skip if i,j neighbor cutoff is less than bin distance + // stores own/own pairs only once + // stores own/ghost pairs on both procs + + ibin = coord2bin(x[i]); + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + if (j <= i) continue; + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + + 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; +} + +/* ---------------------------------------------------------------------- + binned neighbor list construction with full Newton's 3rd law + each owned atom i checks its own bin and other bins in Newton stencil + multi-type stencil is itype dependent and is distance checked + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void Neighbor::half_multi_newton_omp(NeighList *list) +{ + // bin local & ghost atoms + + bin_atoms(); + + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + + NEIGH_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,k,n,itype,jtype,ibin,which,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr,*s; + double *cutsq,*distsq; + + // loop over each atom, storing neighbors + + 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 molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int *nstencil_multi = list->nstencil_multi; + int **stencil_multi = list->stencil_multi; + double **distsq_multi = list->distsq_multi; + + // each thread works on its own page + int npage = tid; + int npnt = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + neighptr = &(list->pages[npage][npnt]); + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over rest of atoms in i's bin, ghosts are at end of linked list + // if j is owned atom, store it, since j is beyond i in linked list + // if j is ghost, only store if j coords are "above and to the right" of i + + for (j = bins[i]; j >= 0; j = bins[j]) { + if (j >= nlocal) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + + // 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]); + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + + 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; +} + +/* ---------------------------------------------------------------------- + binned neighbor list construction with Newton's 3rd law for triclinic + each owned atom i checks its own bin and other bins in triclinic stencil + multi-type stencil is itype dependent and is distance checked + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void Neighbor::half_multi_newton_tri_omp(NeighList *list) +{ + // bin local & ghost atoms + + bin_atoms(); + + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + + NEIGH_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,k,n,itype,jtype,ibin,which,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr,*s; + double *cutsq,*distsq; + + // loop over each atom, storing neighbors + + 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 molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int *nstencil_multi = list->nstencil_multi; + int **stencil_multi = list->stencil_multi; + double **distsq_multi = list->distsq_multi; + + // each thread works on its own page + int npage = tid; + int npnt = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + neighptr = &(list->pages[npage][npnt]); + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over all atoms in bins, including self, in stencil + // skip if i,j neighbor cutoff is less than bin distance + // bins below self are excluded from stencil + // pairs for atoms j "below" i are excluded + // below = lower z or (equal z and lower y) or (equal zy and lower x) + // (equal zyx and j <= i) + // latter excludes self-self interaction but allows superposed atoms + + ibin = coord2bin(x[i]); + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp) { + if (x[j][0] < xtmp) continue; + if (x[j][0] == xtmp && j <= i) continue; + } + } + + 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; +} diff --git a/src/USER-OMP/neigh_half_nsq_omp.cpp b/src/USER-OMP/neigh_half_nsq_omp.cpp new file mode 100644 index 0000000000..0d389b45f6 --- /dev/null +++ b/src/USER-OMP/neigh_half_nsq_omp.cpp @@ -0,0 +1,222 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "neighbor.h" +#include "neighbor_omp.h" +#include "neigh_list.h" +#include "atom.h" +#include "comm.h" +#include "group.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- + N^2 / 2 search for neighbor pairs with partial Newton's 3rd law + pair stored once if i,j are both owned and i < j + pair stored by me if j is ghost (also stored by proc owning j) +------------------------------------------------------------------------- */ + +void Neighbor::half_nsq_no_newton_omp(NeighList *list) +{ + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0; + + NEIGH_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,n,itype,jtype,which; + 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; + + double **x = atom->x; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int nall = atom->nlocal + atom->nghost; + int molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + + int npage = tid; + int npnt = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + neighptr = &(list->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 + + 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; +} + +/* ---------------------------------------------------------------------- + N^2 / 2 search for neighbor pairs with full Newton's 3rd law + every pair stored exactly once by some processor + decision on ghost atoms based on itag,jtag tests +------------------------------------------------------------------------- */ + +void Neighbor::half_nsq_newton_omp(NeighList *list) +{ + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0; + + NEIGH_OMP_INIT; +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,n,itype,jtype,itag,jtag,which; + 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; + + double **x = atom->x; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int nall = atom->nlocal + atom->nghost; + int molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + + int npage = tid; + int npnt = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage >= list->maxpage) list->add_pages(nthreads); + } + + neighptr = &(list->pages[npage][npnt]); + n = 0; + + itag = tag[i]; + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over remaining atoms, owned and ghost + // itag = jtag is possible for long cutoffs that include images of self + + for (j = i+1; j < nall; j++) { + if (includegroup && !(mask[j] & bitmask)) continue; + + if (j >= nlocal) { + 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) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; +} diff --git a/src/USER-OMP/neigh_respa_omp.cpp b/src/USER-OMP/neigh_respa_omp.cpp new file mode 100644 index 0000000000..f8487244a5 --- /dev/null +++ b/src/USER-OMP/neigh_respa_omp.cpp @@ -0,0 +1,985 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "neighbor.h" +#include "neighbor_omp.h" +#include "neigh_list.h" +#include "atom.h" +#include "comm.h" +#include "group.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- + multiple respa lists + N^2 / 2 search for neighbor pairs with partial Newton's 3rd law + pair added to list if atoms i and j are both owned and i < j + pair added if j is ghost (also stored by proc owning j) +------------------------------------------------------------------------- */ + +void Neighbor::respa_nsq_no_newton_omp(NeighList *list) +{ + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0; + + NEIGH_OMP_INIT; + + NeighList *listinner = list->listinner; + if (nthreads > listinner->maxpage) + listinner->add_pages(nthreads - listinner->maxpage); + + NeighList *listmiddle; + const int respamiddle = list->respamiddle; + if (respamiddle) { + listmiddle = list->listmiddle; + if (nthreads > listmiddle->maxpage) + listmiddle->add_pages(nthreads - listmiddle->maxpage); + } + +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list,listinner,listmiddle) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,n,itype,jtype,which,n_inner,n_middle; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr,*neighptr_inner,*neighptr_middle; + + // loop over each atom, storing neighbors + + 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 nall = atom->nlocal + atom->nghost; + int molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + + int *ilist_inner = listinner->ilist; + int *numneigh_inner = listinner->numneigh; + int **firstneigh_inner = listinner->firstneigh; + + int *ilist_middle,*numneigh_middle,**firstneigh_middle; + if (respamiddle) { + ilist_middle = listmiddle->ilist; + numneigh_middle = listmiddle->numneigh; + firstneigh_middle = listmiddle->firstneigh; + } + + int npage = tid; + int npnt = 0; + int npage_inner = tid; + int npnt_inner = 0; + int npage_middle = tid; + int npnt_middle = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage == list->maxpage) list->add_pages(nthreads); + } + neighptr = &(list->pages[npage][npnt]); + n = 0; + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt_inner < oneatom) { + npnt_inner = 0; + npage_inner += nthreads; + if (npage_inner == listinner->maxpage) listinner->add_pages(nthreads); + } + neighptr_inner = &(listinner->pages[npage_inner][npnt_inner]); + n_inner = 0; + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (respamiddle) { + if (pgsize - npnt_middle < oneatom) { + npnt_middle = 0; + npage_middle += nthreads; + if (npage_middle == listmiddle->maxpage) listmiddle->add_pages(nthreads); + } + neighptr_middle = &(listmiddle->pages[npage_middle][npnt_middle]); + n_middle = 0; + } + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over remaining atoms, owned and ghost + + 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + + if (rsq < cut_inner_sq) { + if (which == 0) neighptr_inner[n_inner++] = j; + else if (which > 0) neighptr_inner[n_inner++] = j ^ (which << SBBITS); + } + + if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { + if (which == 0) neighptr_middle[n_middle++] = j; + else if (which > 0) + neighptr_middle[n_middle++] = j ^ (which << SBBITS); + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + + ilist_inner[i] = i; + firstneigh_inner[i] = neighptr_inner; + numneigh_inner[i] = n_inner; + npnt_inner += n_inner; + if (npnt_inner >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + + if (respamiddle) { + ilist_middle[i] = i; + firstneigh_middle[i] = neighptr_middle; + numneigh_middle[i] = n_middle; + npnt_middle += n_middle; + if (npnt_middle >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; + listinner->inum = nlocal; + if (respamiddle) listmiddle->inum = nlocal; +} + +/* ---------------------------------------------------------------------- + multiple respa lists + N^2 / 2 search for neighbor pairs with full Newton's 3rd law + pair added to list if atoms i and j are both owned and i < j + if j is ghost only me or other proc adds pair + decision based on itag,jtag tests +------------------------------------------------------------------------- */ + +void Neighbor::respa_nsq_newton_omp(NeighList *list) +{ + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0; + + NEIGH_OMP_INIT; + + NeighList *listinner = list->listinner; + if (nthreads > listinner->maxpage) + listinner->add_pages(nthreads - listinner->maxpage); + + NeighList *listmiddle; + const int respamiddle = list->respamiddle; + if (respamiddle) { + listmiddle = list->listmiddle; + if (nthreads > listmiddle->maxpage) + listmiddle->add_pages(nthreads - listmiddle->maxpage); + } + +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list,listinner,listmiddle) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,n,itype,jtype,itag,jtag,which,n_inner,n_middle; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr,*neighptr_inner,*neighptr_middle; + + // loop over each atom, storing neighbors + + 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 nall = atom->nlocal + atom->nghost; + int molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + + int *ilist_inner = listinner->ilist; + int *numneigh_inner = listinner->numneigh; + int **firstneigh_inner = listinner->firstneigh; + + int *ilist_middle,*numneigh_middle,**firstneigh_middle; + if (respamiddle) { + ilist_middle = listmiddle->ilist; + numneigh_middle = listmiddle->numneigh; + firstneigh_middle = listmiddle->firstneigh; + } + + int npage = tid; + int npnt = 0; + int npage_inner = tid; + int npnt_inner = 0; + int npage_middle = tid; + int npnt_middle = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage == list->maxpage) list->add_pages(nthreads); + } + neighptr = &(list->pages[npage][npnt]); + n = 0; + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt_inner < oneatom) { + npnt_inner = 0; + npage_inner += nthreads; + if (npage_inner == listinner->maxpage) listinner->add_pages(nthreads); + } + neighptr_inner = &(listinner->pages[npage_inner][npnt_inner]); + n_inner = 0; + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (respamiddle) { + if (pgsize - npnt_middle < oneatom) { + npnt_middle = 0; + npage_middle += nthreads; + if (npage_middle == listmiddle->maxpage) listmiddle->add_pages(nthreads); + } + neighptr_middle = &(listmiddle->pages[npage_middle][npnt_middle]); + n_middle = 0; + } + + itag = tag[i]; + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over remaining atoms, owned and ghost + + for (j = i+1; j < nall; j++) { + if (includegroup && !(mask[j] & bitmask)) continue; + + if (j >= nlocal) { + 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) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + + if (rsq < cut_inner_sq) { + if (which == 0) neighptr_inner[n_inner++] = j; + else if (which > 0) neighptr_inner[n_inner++] = j ^ (which << SBBITS); + } + + if (respamiddle && + rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { + if (which == 0) neighptr_middle[n_middle++] = j; + else if (which > 0) + neighptr_middle[n_middle++] = j ^ (which << SBBITS); + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + + ilist_inner[i] = i; + firstneigh_inner[i] = neighptr_inner; + numneigh_inner[i] = n_inner; + npnt_inner += n_inner; + if (npnt_inner >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + + if (respamiddle) { + ilist_middle[i] = i; + firstneigh_middle[i] = neighptr_middle; + numneigh_middle[i] = n_middle; + npnt_middle += n_middle; + if (npnt_middle >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; + listinner->inum = nlocal; + if (respamiddle) listmiddle->inum = nlocal; +} + +/* ---------------------------------------------------------------------- + multiple respa lists + binned neighbor list construction with partial Newton's 3rd law + each owned atom i checks own bin and surrounding bins in non-Newton stencil + pair stored once if i,j are both owned and i < j + pair stored by me if j is ghost (also stored by proc owning j) +------------------------------------------------------------------------- */ + +void Neighbor::respa_bin_no_newton_omp(NeighList *list) +{ + // bin local & ghost atoms + + bin_atoms(); + + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + + NEIGH_OMP_INIT; + + NeighList *listinner = list->listinner; + if (nthreads > listinner->maxpage) + listinner->add_pages(nthreads - listinner->maxpage); + + NeighList *listmiddle; + const int respamiddle = list->respamiddle; + if (respamiddle) { + listmiddle = list->listmiddle; + if (nthreads > listmiddle->maxpage) + listmiddle->add_pages(nthreads - listmiddle->maxpage); + } + +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list,listinner,listmiddle) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,k,n,itype,jtype,ibin,which,n_inner,n_middle; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr,*neighptr_inner,*neighptr_middle; + + // loop over each atom, storing neighbors + + 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 molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int nstencil = list->nstencil; + int *stencil = list->stencil; + + int *ilist_inner = listinner->ilist; + int *numneigh_inner = listinner->numneigh; + int **firstneigh_inner = listinner->firstneigh; + + int *ilist_middle,*numneigh_middle,**firstneigh_middle; + if (respamiddle) { + ilist_middle = listmiddle->ilist; + numneigh_middle = listmiddle->numneigh; + firstneigh_middle = listmiddle->firstneigh; + } + + int npage = tid; + int npnt = 0; + int npage_inner = tid; + int npnt_inner = 0; + int npage_middle = tid; + int npnt_middle = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage == list->maxpage) list->add_pages(nthreads); + } + neighptr = &(list->pages[npage][npnt]); + n = 0; + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt_inner < oneatom) { + npnt_inner = 0; + npage_inner += nthreads; + if (npage_inner == listinner->maxpage) listinner->add_pages(nthreads); + } + neighptr_inner = &(listinner->pages[npage_inner][npnt_inner]); + n_inner = 0; + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (respamiddle) { + if (pgsize - npnt_middle < oneatom) { + npnt_middle = 0; + npage_middle += nthreads; + if (npage_middle == listmiddle->maxpage) listmiddle->add_pages(nthreads); + } + neighptr_middle = &(listmiddle->pages[npage_middle][npnt_middle]); + n_middle = 0; + } + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + ibin = coord2bin(x[i]); + + // loop over all atoms in surrounding bins in stencil including self + // only store pair if i < j + // stores own/own pairs only once + // stores own/ghost pairs on both procs + + for (k = 0; k < nstencil; k++) { + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (j <= i) 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + + if (rsq < cut_inner_sq) { + if (which == 0) neighptr_inner[n_inner++] = j; + else if (which > 0) + neighptr_inner[n_inner++] = j ^ (which << SBBITS); + } + + if (respamiddle && + rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { + if (which == 0) neighptr_middle[n_middle++] = j; + else if (which > 0) + neighptr_middle[n_middle++] = j ^ (which << SBBITS); + } + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + + ilist_inner[i] = i; + firstneigh_inner[i] = neighptr_inner; + numneigh_inner[i] = n_inner; + npnt_inner += n_inner; + if (npnt_inner >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + + if (respamiddle) { + ilist_middle[i] = i; + firstneigh_middle[i] = neighptr_middle; + numneigh_middle[i] = n_middle; + npnt_middle += n_middle; + if (npnt_middle >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; + listinner->inum = nlocal; + if (respamiddle) listmiddle->inum = nlocal; +} + +/* ---------------------------------------------------------------------- + multiple respa lists + binned neighbor list construction with full Newton's 3rd law + each owned atom i checks its own bin and other bins in Newton stencil + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void Neighbor::respa_bin_newton_omp(NeighList *list) +{ + // bin local & ghost atoms + + bin_atoms(); + + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + + NEIGH_OMP_INIT; + + NeighList *listinner = list->listinner; + if (nthreads > listinner->maxpage) + listinner->add_pages(nthreads - listinner->maxpage); + + NeighList *listmiddle; + const int respamiddle = list->respamiddle; + if (respamiddle) { + listmiddle = list->listmiddle; + if (nthreads > listmiddle->maxpage) + listmiddle->add_pages(nthreads - listmiddle->maxpage); + } + +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list,listinner,listmiddle) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,k,n,itype,jtype,ibin,which,n_inner,n_middle; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr,*neighptr_inner,*neighptr_middle; + + // loop over each atom, storing neighbors + + 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 molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int nstencil = list->nstencil; + int *stencil = list->stencil; + + int *ilist_inner = listinner->ilist; + int *numneigh_inner = listinner->numneigh; + int **firstneigh_inner = listinner->firstneigh; + + int *ilist_middle,*numneigh_middle,**firstneigh_middle; + if (respamiddle) { + ilist_middle = listmiddle->ilist; + numneigh_middle = listmiddle->numneigh; + firstneigh_middle = listmiddle->firstneigh; + } + + int npage = tid; + int npnt = 0; + int npage_inner = tid; + int npnt_inner = 0; + int npage_middle = tid; + int npnt_middle = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage == list->maxpage) list->add_pages(nthreads); + } + neighptr = &(list->pages[npage][npnt]); + n = 0; + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt_inner < oneatom) { + npnt_inner = 0; + npage_inner += nthreads; + if (npage_inner == listinner->maxpage) listinner->add_pages(nthreads); + } + neighptr_inner = &(listinner->pages[npage_inner][npnt_inner]); + n_inner = 0; + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (respamiddle) { + if (pgsize - npnt_middle < oneatom) { + npnt_middle = 0; + npage_middle += nthreads; + if (npage_middle == listmiddle->maxpage) listmiddle->add_pages(nthreads); + } + neighptr_middle = &(listmiddle->pages[npage_middle][npnt_middle]); + n_middle = 0; + } + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over rest of atoms in i's bin, ghosts are at end of linked list + // if j is owned atom, store it, since j is beyond i in linked list + // if j is ghost, only store if j coords are "above and to the right" of i + + for (j = bins[i]; j >= 0; j = bins[j]) { + if (j >= nlocal) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + + if (rsq < cut_inner_sq) { + if (which == 0) neighptr_inner[n_inner++] = j; + else if (which > 0) neighptr_inner[n_inner++] = j ^ (which << SBBITS); + } + + if (respamiddle && + rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { + if (which == 0) neighptr_middle[n_middle++] = j; + else if (which > 0) + neighptr_middle[n_middle++] = j ^ (which << SBBITS); + } + } + } + + // loop over all atoms in other bins in stencil, store every pair + + ibin = coord2bin(x[i]); + for (k = 0; k < nstencil; k++) { + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + + if (rsq < cut_inner_sq) { + if (which == 0) neighptr_inner[n_inner++] = j; + else if (which > 0) + neighptr_inner[n_inner++] = j ^ (which << SBBITS); + } + + if (respamiddle && + rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { + if (which == 0) neighptr_middle[n_middle++] = j; + else if (which > 0) + neighptr_middle[n_middle++] = j ^ (which << SBBITS); + } + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + + ilist_inner[i] = i; + firstneigh_inner[i] = neighptr_inner; + numneigh_inner[i] = n_inner; + npnt_inner += n_inner; + if (npnt_inner >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + + if (respamiddle) { + ilist_middle[i] = i; + firstneigh_middle[i] = neighptr_middle; + numneigh_middle[i] = n_middle; + npnt_middle += n_middle; + if (npnt_middle >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; + listinner->inum = nlocal; + if (respamiddle) listmiddle->inum = nlocal; +} + +/* ---------------------------------------------------------------------- + multiple respa lists + binned neighbor list construction with Newton's 3rd law for triclinic + each owned atom i checks its own bin and other bins in triclinic stencil + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void Neighbor::respa_bin_newton_tri_omp(NeighList *list) +{ + // bin local & ghost atoms + + bin_atoms(); + + const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + + NEIGH_OMP_INIT; + + NeighList *listinner = list->listinner; + if (nthreads > listinner->maxpage) + listinner->add_pages(nthreads - listinner->maxpage); + + NeighList *listmiddle; + const int respamiddle = list->respamiddle; + if (respamiddle) { + listmiddle = list->listmiddle; + if (nthreads > listmiddle->maxpage) + listmiddle->add_pages(nthreads - listmiddle->maxpage); + } + +#if defined(_OPENMP) +#pragma omp parallel default(none) shared(list,listinner,listmiddle) +#endif + NEIGH_OMP_SETUP(nlocal); + + int i,j,k,n,itype,jtype,ibin,which,n_inner,n_middle; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr,*neighptr_inner,*neighptr_middle; + + // loop over each atom, storing neighbors + + 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 molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int nstencil = list->nstencil; + int *stencil = list->stencil; + + int *ilist_inner = listinner->ilist; + int *numneigh_inner = listinner->numneigh; + int **firstneigh_inner = listinner->firstneigh; + + int *ilist_middle,*numneigh_middle,**firstneigh_middle; + if (respamiddle) { + ilist_middle = listmiddle->ilist; + numneigh_middle = listmiddle->numneigh; + firstneigh_middle = listmiddle->firstneigh; + } + + int npage = tid; + int npnt = 0; + int npage_inner = tid; + int npnt_inner = 0; + int npage_middle = tid; + int npnt_middle = 0; + + for (i = ifrom; i < ito; i++) { + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt < oneatom) { + npnt = 0; + npage += nthreads; + if (npage == list->maxpage) list->add_pages(nthreads); + } + neighptr = &(list->pages[npage][npnt]); + n = 0; + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (pgsize - npnt_inner < oneatom) { + npnt_inner = 0; + npage_inner += nthreads; + if (npage_inner == listinner->maxpage) listinner->add_pages(nthreads); + } + neighptr_inner = &(listinner->pages[npage_inner][npnt_inner]); + n_inner = 0; + +#if defined(_OPENMP) +#pragma omp critical +#endif + if (respamiddle) { + if (pgsize - npnt_middle < oneatom) { + npnt_middle = 0; + npage_middle += nthreads; + if (npage_middle == listmiddle->maxpage) listmiddle->add_pages(nthreads); + } + neighptr_middle = &(listmiddle->pages[npage_middle][npnt_middle]); + n_middle = 0; + } + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over all atoms in bins in stencil + // pairs for atoms j "below" i are excluded + // below = lower z or (equal z and lower y) or (equal zy and lower x) + // (equal zyx and j <= i) + // latter excludes self-self interaction but allows superposed atoms + + ibin = coord2bin(x[i]); + for (k = 0; k < nstencil; k++) { + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp) { + if (x[j][0] < xtmp) continue; + if (x[j][0] == xtmp && j <= i) 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 ^ (which << SBBITS); + } else neighptr[n++] = j; + + if (rsq < cut_inner_sq) { + if (which == 0) neighptr_inner[n_inner++] = j; + else if (which > 0) + neighptr_inner[n_inner++] = j ^ (which << SBBITS); + } + + if (respamiddle && + rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { + if (which == 0) neighptr_middle[n_middle++] = j; + else if (which > 0) + neighptr_middle[n_middle++] = j ^ (which << SBBITS); + } + } + } + } + + ilist[i] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom || npnt >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + + ilist_inner[i] = i; + firstneigh_inner[i] = neighptr_inner; + numneigh_inner[i] = n_inner; + npnt_inner += n_inner; + if (npnt_inner >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + + if (respamiddle) { + ilist_middle[i] = i; + firstneigh_middle[i] = neighptr_middle; + numneigh_middle[i] = n_middle; + npnt_middle += n_middle; + if (npnt_middle >= pgsize) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + } + } + NEIGH_OMP_CLOSE; + list->inum = nlocal; + listinner->inum = nlocal; + if (respamiddle) listmiddle->inum = nlocal; +} diff --git a/src/USER-OMP/neighbor_omp.h b/src/USER-OMP/neighbor_omp.h new file mode 100644 index 0000000000..8cdd61d5b1 --- /dev/null +++ b/src/USER-OMP/neighbor_omp.h @@ -0,0 +1,60 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_NEIGHBOR_OMP_H +#define LMP_NEIGHBOR_OMP_H + +#if defined(_OPENMP) +#include +#endif + +namespace LAMMPS_NS { + +// these macros hide some ugly and redundant OpenMP related stuff +#if defined(_OPENMP) + +// make sure we have at least one page for each thread +#define NEIGH_OMP_INIT \ + const int nthreads = comm->nthreads; \ + if (nthreads > list->maxpage) \ + list->add_pages(nthreads - list->maxpage) + +// get thread id and then assign each thread a fixed chunk of atoms +#define NEIGH_OMP_SETUP(num) \ + { \ + const int tid = omp_get_thread_num(); \ + const int idelta = 1 + num/nthreads; \ + const int ifrom = tid*idelta; \ + int ito = ifrom + idelta; \ + if (ito > num) \ + ito = num + +#define NEIGH_OMP_CLOSE } + +#else /* !defined(_OPENMP) */ + +#define NEIGH_OMP_INIT \ + const int nthreads = comm->nthreads; + +#define NEIGH_OMP_SETUP(num) \ + const int tid = 0; \ + const int ifrom = 0; \ + const int ito = num + +#define NEIGH_OMP_CLOSE + +#endif + +} + +#endif