From 061229093cb7c28058fdbdcff2fd13c9226d78de Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Tue, 10 Nov 2020 22:44:55 -0700 Subject: [PATCH] Adding npair classes --- src/npair.cpp | 65 +++++ src/npair.h | 14 ++ ...multi_tiered.cpp => npair_full_multi2.cpp} | 83 ++++--- ...ull_multi_tiered.h => npair_full_multi2.h} | 16 +- src/npair_half_multi2_newtoff.cpp | 143 +++++++++++ src/npair_half_multi2_newtoff.h | 47 ++++ src/npair_half_multi2_newton.cpp | 225 +++++++++--------- src/npair_half_multi2_newton.h | 16 +- src/npair_half_multi2_newton_tri.cpp | 172 +++++++------ src/npair_half_multi2_newton_tri.h | 16 +- src/npair_half_size_multi2_newtoff.cpp | 68 +++--- src/npair_half_size_multi2_newtoff.h | 16 +- src/npair_half_size_multi2_newton.cpp | 183 +++++++------- src/npair_half_size_multi2_newton.h | 16 +- src/npair_half_size_multi2_newton_tri.cpp | 145 ++++++----- src/npair_half_size_multi2_newton_tri.h | 16 +- src/nstencil.h | 4 +- 17 files changed, 738 insertions(+), 507 deletions(-) rename src/{npair_full_multi_tiered.cpp => npair_full_multi2.cpp} (63%) rename src/{npair_full_multi_tiered.h => npair_full_multi2.h} (75%) create mode 100755 src/npair_half_multi2_newtoff.cpp create mode 100755 src/npair_half_multi2_newtoff.h diff --git a/src/npair.cpp b/src/npair.cpp index 698ff28067..920728a692 100644 --- a/src/npair.cpp +++ b/src/npair.cpp @@ -134,6 +134,25 @@ void NPair::copy_bin_info() atom2bin = nb->atom2bin; bins = nb->bins; binhead = nb->binhead; + + nbinx_multi2 = nb->nbinx_multi2; + nbiny_multi2 = nb->nbiny_multi2; + nbinz_multi2 = nb->nbinz_multi2; + mbins_multi2 = nb->mbins_multi2; + mbinx_multi2 = nb->mbinx_multi2; + mbiny_multi2 = nb->mbiny_multi2; + mbinz_multi2 = nb->mbinz_multi2; + mbinxlo_multi2 = nb->mbinxlo_multi2; + mbinylo_multi2 = nb->mbinylo_multi2; + mbinzlo_multi2 = nb->mbinzlo_multi2; + + bininvx_multi2 = nb->bininvx_multi2; + bininvy_multi2 = nb->bininvy_multi2; + bininvz_multi2 = nb->bininvz_multi2; + + atom2bin_multi2 = nb->atom2bin_multi2; + bins_multi2 = nb->bins_multi2; + binhead_multi2 = nb->binhead_multi2; } /* ---------------------------------------------------------------------- @@ -148,6 +167,9 @@ void NPair::copy_stencil_info() nstencil_multi = ns->nstencil_multi; stencil_multi = ns->stencil_multi; distsq_multi = ns->distsq_multi; + + nstencil_multi2 = ns->nstencil_multi2; + stencil_multi2 = ns->stencil_multi2; } /* ---------------------------------------------------------------------- @@ -241,3 +263,46 @@ int NPair::coord2bin(double *x, int &ix, int &iy, int &iz) return iz*mbiny*mbinx + iy*mbinx + ix; } + +/* ---------------------------------------------------------------------- + same as coord2bin in NbinMulti2 +------------------------------------------------------------------------- */ + +int NPair::coord2bin(double *x, int it) +{ + int ix,iy,iz; + int ibin; + + if (!std::isfinite(x[0]) || !std::isfinite(x[1]) || !std::isfinite(x[2])) + error->one(FLERR,"Non-numeric positions - simulation unstable"); + + if (x[0] >= bboxhi[0]) + ix = static_cast ((x[0]-bboxhi[0])*bininvx_multi2[it]) + nbinx_multi2[it]; + else if (x[0] >= bboxlo[0]) { + ix = static_cast ((x[0]-bboxlo[0])*bininvx_multi2[it]); + ix = MIN(ix,nbinx_multi2[it]-1); + } else + ix = static_cast ((x[0]-bboxlo[0])*bininvx_multi2[it]) - 1; + + if (x[1] >= bboxhi[1]) + iy = static_cast ((x[1]-bboxhi[1])*bininvy_multi2[it]) + nbiny_multi2[it]; + else if (x[1] >= bboxlo[1]) { + iy = static_cast ((x[1]-bboxlo[1])*bininvy_multi2[it]); + iy = MIN(iy,nbiny_multi2[it]-1); + } else + iy = static_cast ((x[1]-bboxlo[1])*bininvy_multi2[it]) - 1; + + if (x[2] >= bboxhi[2]) + iz = static_cast ((x[2]-bboxhi[2])*bininvz_multi2[it]) + nbinz_multi2[it]; + else if (x[2] >= bboxlo[2]) { + iz = static_cast ((x[2]-bboxlo[2])*bininvz_multi2[it]); + iz = MIN(iz,nbinz_multi2[it]-1); + } else + iz = static_cast ((x[2]-bboxlo[2])*bininvz_multi2[it]) - 1; + + + ibin = (iz-mbinzlo_multi2[it])*mbiny_multi2[it]*mbinx_multi2[it] + + (iy-mbinylo_multi2[it])*mbinx_multi2[it] + + (ix-mbinxlo_multi2[it]); + return ibin; +} \ No newline at end of file diff --git a/src/npair.h b/src/npair.h index ec75042302..87b2584521 100644 --- a/src/npair.h +++ b/src/npair.h @@ -79,6 +79,14 @@ class NPair : protected Pointers { int *atom2bin,*bins; int *binhead; + int *nbinx_multi2, *nbiny_multi2, *nbinz_multi2; + int *mbins_multi2; + int *mbinx_multi2, *mbiny_multi2, *mbinz_multi2; + int *mbinxlo_multi2, *mbinylo_multi2, *mbinzlo_multi2; + double *bininvx_multi2, *bininvy_multi2, *bininvz_multi2; + int **binhead_multi2,**bins_multi2; + int **atom2bin_multi2; + // data from NStencil class int nstencil; @@ -88,6 +96,9 @@ class NPair : protected Pointers { int **stencil_multi; double **distsq_multi; + int ** nstencil_multi2; + int *** stencil_multi2; + // data common to all NPair variants int molecular; @@ -102,6 +113,9 @@ class NPair : protected Pointers { int coord2bin(double *); // mapping atom coord to a bin int coord2bin(double *, int &, int &, int&); // ditto + int coord2bin(double *, int); // mapping atom coord to type bin + + // find_special: determine if atom j is in special list of atom i // if it is not, return 0 // if it is and special flag is 0 (both coeffs are 0.0), return -1 diff --git a/src/npair_full_multi_tiered.cpp b/src/npair_full_multi2.cpp similarity index 63% rename from src/npair_full_multi_tiered.cpp rename to src/npair_full_multi2.cpp index 3b9a89bb56..4fbdc6148b 100644 --- a/src/npair_full_multi_tiered.cpp +++ b/src/npair_full_multi2.cpp @@ -11,7 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include "npair_full_bytype.h" +#include "npair_full_multi2.h" #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" @@ -19,28 +19,27 @@ #include "domain.h" #include "my_page.h" #include "error.h" -#include "nbin.h" -#include "nstencil.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NPairFullBytype::NPairFullBytype(LAMMPS *lmp) : NPair(lmp) {} +NPairFullMulti2::NPairFullMulti2(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- + REWRITE 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 - KS ADJUST ------------------------------------------------------------------------- */ -void NPairFullBytype::build(NeighList *list) +void NPairFullMulti2::build(NeighList *list) { - int i,j,k,n,itype,jtype,ibin,which,ns,imol,iatom,moltemplate; + int i,j,k,n,itype,jtype,ktype,ibin,kbin,which,ns,imol,iatom,moltemplate; tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; + int js; double **x = atom->x; int *type = atom->type; @@ -84,48 +83,46 @@ void NPairFullBytype::build(NeighList *list) // skip if i,j neighbor cutoff is less than bin distance // skip i = j - int kbin; - ibin = nb->atom2bin_type[itype][i]; - for (int ktype = 1; ktype <= atom->ntypes; ktype++) { + ibin = atom2bin_multi2[itype][i]; + for (ktype = 1; ktype <= atom->ntypes; ktype++) { if (itype == ktype) { - kbin = ibin; - } - else { - // Locate i in ktype bin - kbin = nb->coord2bin(x[i], ktype); + kbin = ibin; + } else { + // Locate i in ktype bin + kbin = coord2bin(x[i], ktype); } - s = this->ns->stencil_type[itype][ktype]; - ns = this->ns->nstencil_type[itype][ktype]; + s = stencil_multi2[itype][ktype]; + ns = nstencil_multi2[itype][ktype]; for (k = 0; k < ns; k++) { - int js = this->nb->binhead_type[ktype][kbin + s[k]]; - for (j = js; j >= 0; j = this->nb->bins_type[ktype][j]) { - if (i == j) continue; + js = binhead_multi2[ktype][kbin + s[k]]; + for (j = js; j >= 0; j = bins_multi2[ktype][j]) { + if (i == j) continue; - jtype = type[j]; - if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) 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 != Atom::ATOMIC) { - if (!moltemplate) - which = find_special(special[i],nspecial[i],tag[j]); - else if (imol >= 0) - which = find_special(onemols[imol]->special[iatom], - onemols[imol]->nspecial[iatom], - tag[j]-tagprev); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (domain->minimum_image_check(delx,dely,delz)) - neighptr[n++] = j; - else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); - } else neighptr[n++] = j; - } - } + 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 != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } } } diff --git a/src/npair_full_multi_tiered.h b/src/npair_full_multi2.h similarity index 75% rename from src/npair_full_multi_tiered.h rename to src/npair_full_multi2.h index fdcd6b5fb1..f552e5bf47 100644 --- a/src/npair_full_multi_tiered.h +++ b/src/npair_full_multi2.h @@ -13,23 +13,23 @@ #ifdef NPAIR_CLASS -NPairStyle(full/bytype, - NPairFullBytype, - NP_FULL | NP_BYTYPE | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI) +NPairStyle(full/multi2, + NPairFullMulti2, + NP_FULL | NP_MULTI2 | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI) #else -#ifndef LMP_NPAIR_FULL_BYTYPE_H -#define LMP_NPAIR_FULL_BYTYPE_H +#ifndef LMP_NPAIR_FULL_MULTI2_H +#define LMP_NPAIR_FULL_MULTI2_H #include "npair.h" namespace LAMMPS_NS { -class NPairFullBytype : public NPair { +class NPairFullMulti2 : public NPair { public: - NPairFullBytype(class LAMMPS *); - ~NPairFullBytype() {} + NPairFullMulti2(class LAMMPS *); + ~NPairFullMulti2() {} void build(class NeighList *); }; diff --git a/src/npair_half_multi2_newtoff.cpp b/src/npair_half_multi2_newtoff.cpp new file mode 100755 index 0000000000..83a8054247 --- /dev/null +++ b/src/npair_half_multi2_newtoff.cpp @@ -0,0 +1,143 @@ +/* ---------------------------------------------------------------------- + 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 "npair_half_multi2_newtoff.h" +#include "neigh_list.h" +#include "atom.h" +#include "atom_vec.h" +#include "molecule.h" +#include "domain.h" +#include "my_page.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NPairHalfMulti2Newtoff::NPairHalfMulti2Newtoff(LAMMPS *lmp) : NPair(lmp) {} + +/* ---------------------------------------------------------------------- + REWRITE + 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 NPairHalfMulti2Newtoff::build(NeighList *list) +{ + int i,j,k,n,itype,jtype,ktype,ibin,kbin,which,ns,imol,iatom,moltemplate; + tagint tagprev; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr,*s; + int js; + + double **x = atom->x; + int *type = atom->type; + int *mask = atom->mask; + tagint *tag = atom->tag; + tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + int nlocal = atom->nlocal; + if (includegroup) nlocal = atom->nfirst; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + if (molecular == Atom::TEMPLATE) moltemplate = 1; + else moltemplate = 0; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + MyPage *ipage = list->ipage; + + int inum = 0; + ipage->reset(); + + for (i = 0; i < nlocal; i++) { + n = 0; + neighptr = ipage->vget(); + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } + + // 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 = nb->atom2bin_type[itype][i]; + for (ktype = 1; ktype <= atom->ntypes; ktype++) { + if (itype == ktype) { + kbin = ibin; + } else { + // Locate i in ktype bin + kbin = coord2bin(x[i], ktype); + } + + s = stencil_multi2[itype][ktype]; + ns = nstencil_multi2[itype][ktype]; + for (k = 0; k < ns; k++) { + js = binhead_multi2[ktype][kbin + s[k]]; + for (j = js; j >=0; j = nb->bins_multi2[ktype][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 != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + } + } + + ilist[inum++] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + ipage->vgot(n); + if (ipage->status()) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + } + + list->inum = inum; +} diff --git a/src/npair_half_multi2_newtoff.h b/src/npair_half_multi2_newtoff.h new file mode 100755 index 0000000000..30a6d3164d --- /dev/null +++ b/src/npair_half_multi2_newtoff.h @@ -0,0 +1,47 @@ +/* -*- c++ -*- ---------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef NPAIR_CLASS + +NPairStyle(half/multi2/newtoff, + NPairHalfMulti2Newtoff, + NP_HALF | NP_MULTI2 | NP_NEWTOFF | NP_ORTHO | NP_TRI) + +#else + +#ifndef LMP_NPAIR_HALF_MULTI2_NEWTOFF_H +#define LMP_NPAIR_HALF_MULTI2_NEWTOFF_H + +#include "npair.h" + +namespace LAMMPS_NS { + +class NPairHalfMulti2Newtoff : public NPair { + public: + NPairHalfMulti2Newtoff(class LAMMPS *); + ~NPairHalfMulti2Newtoff() {} + void build(class NeighList *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Neighbor list overflow, boost neigh_modify one + +UNDOCUMENTED + +*/ diff --git a/src/npair_half_multi2_newton.cpp b/src/npair_half_multi2_newton.cpp index 66899564be..58397f645e 100755 --- a/src/npair_half_multi2_newton.cpp +++ b/src/npair_half_multi2_newton.cpp @@ -11,9 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include -#include "npair_half_bytype_newton.h" -#include "neighbor.h" +#include "npair_half_multi2_newton.h" #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" @@ -22,14 +20,11 @@ #include "my_page.h" #include "error.h" -#include "nbin.h" -#include "nstencil.h" - using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NPairHalfBytypeNewton::NPairHalfBytypeNewton(LAMMPS *lmp) : NPair(lmp) {} +NPairHalfMulti2Newton::NPairHalfMulti2Newton(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- KS REWRTIE @@ -39,12 +34,13 @@ NPairHalfBytypeNewton::NPairHalfBytypeNewton(LAMMPS *lmp) : NPair(lmp) {} every pair stored exactly once by some processor ------------------------------------------------------------------------- */ -void NPairHalfBytypeNewton::build(NeighList *list) +void NPairHalfMulti2Newton::build(NeighList *list) { - int i,j,k,n,itype,jtype,ibin,which,ns,imol,iatom,moltemplate; + int i,j,k,n,itype,jtype,ktype,ibin,kbin,which,ns,imol,iatom,moltemplate; tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; + int js; double **x = atom->x; int *type = atom->type; @@ -84,129 +80,124 @@ void NPairHalfBytypeNewton::build(NeighList *list) tagprev = tag[i] - iatom - 1; } - int js; - int kbin; - // own type: loop over atoms ahead in bin, including ghosts at end of list // if j is owned atom, store by virtue of being ahead of i in list // if j is ghost, store if x[j] "above and to right of" x[i] - ibin = nb->atom2bin_type[itype][i]; - - for (int ktype = 1; ktype <= atom->ntypes; ktype++) { + ibin = atom2bin_multi2[itype][i]; + for (ktype = 1; ktype <= atom->ntypes; ktype++) { if (itype == ktype) { - // own bin ... - js = nb->bins_type[itype][i]; - for (j = js; j >= 0; j = nb->bins_type[itype][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; + // own bin ... + js = bins_multi2[itype][i]; + for (j = js; j >= 0; j = bins_multi2[itype][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) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; + } } - } - 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) { - if (!moltemplate) - which = find_special(special[i],nspecial[i],tag[j]); - else if (imol >= 0) - which = find_special(onemols[imol]->special[iatom], - onemols[imol]->nspecial[iatom], - tag[j]-tagprev); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (domain->minimum_image_check(delx,dely,delz)) - neighptr[n++] = j; - else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); - } else neighptr[n++] = j; - } - } - - // loop over all atoms in other bins in stencil, store every pair - // skip if i,j neighbor cutoff is less than bin distance - - s = this->ns->stencil_type[itype][itype]; - ns = this->ns->nstencil_type[itype][itype]; - for (k = 0; k < ns; k++) { - js = nb->binhead_type[itype][ibin + s[k]]; - for (j = js; j >= 0; j = nb->bins_type[itype][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) { - if (!moltemplate) - which = find_special(special[i],nspecial[i],tag[j]); - else if (imol >= 0) - which = find_special(onemols[imol]->special[iatom], - onemols[imol]->nspecial[iatom], - tag[j]-tagprev); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (domain->minimum_image_check(delx,dely,delz)) - neighptr[n++] = j; - else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); - } else neighptr[n++] = j; + // loop over all atoms in other bins in stencil, store every pair + // skip if i,j neighbor cutoff is less than bin distance + + s = stencil_multi2[itype][itype]; + ns = nstencil_multi2[itype][itype]; + for (k = 0; k < ns; k++) { + js = binhead_multi2[itype][ibin + s[k]]; + for (j = js; j >= 0; j = bins_multi2[itype][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) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } } - } - } - - } - else { + } else { // smaller -> larger: locate i in the ktype bin structure - kbin = nb->coord2bin(x[i], ktype); - s = this->ns->stencil_type[itype][ktype]; - ns = this->ns->nstencil_type[itype][ktype]; - - for (k = 0; k < ns; k++) { - js = nb->binhead_type[ktype][kbin + s[k]]; - for (j = js; j >= 0; j = nb->bins_type[ktype][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) { - if (!moltemplate) - which = find_special(special[i],nspecial[i],tag[j]); - else if (imol >= 0) - which = find_special(onemols[imol]->special[iatom], - onemols[imol]->nspecial[iatom], - tag[j]-tagprev); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (domain->minimum_image_check(delx,dely,delz)) - neighptr[n++] = j; - else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); - } else neighptr[n++] = j; + kbin = coord2bin(x[i], ktype); + s = stencil_multi2[itype][ktype]; + ns = nstencil_multi2[itype][ktype]; + + for (k = 0; k < ns; k++) { + js = binhead_multi2[ktype][kbin + s[k]]; + for (j = js; j >= 0; j = bins_multi2[ktype][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) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } } - } - } } } + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; diff --git a/src/npair_half_multi2_newton.h b/src/npair_half_multi2_newton.h index 8be7292219..8037d2e172 100755 --- a/src/npair_half_multi2_newton.h +++ b/src/npair_half_multi2_newton.h @@ -13,23 +13,23 @@ #ifdef NPAIR_CLASS -NPairStyle(half/bytype/newton, - NPairHalfBytypeNewton, - NP_HALF | NP_BYTYPE | NP_NEWTON | NP_ORTHO) +NPairStyle(half/multi2/newton, + NPairHalfMulti2Newton, + NP_HALF | NP_MULTI2 | NP_NEWTON | NP_ORTHO) #else -#ifndef LMP_NPAIR_HALF_BYTYPE_NEWTON_H -#define LMP_NPAIR_HALF_BYTYPE_NEWTON_H +#ifndef LMP_NPAIR_HALF_MULTI2_NEWTON_H +#define LMP_NPAIR_HALF_MULTI2_NEWTON_H #include "npair.h" namespace LAMMPS_NS { -class NPairHalfBytypeNewton : public NPair { +class NPairHalfMulti2Newton : public NPair { public: - NPairHalfBytypeNewton(class LAMMPS *); - ~NPairHalfBytypeNewton() {} + NPairHalfMulti2Newton(class LAMMPS *); + ~NPairHalfMulti2Newton() {} void build(class NeighList *); }; diff --git a/src/npair_half_multi2_newton_tri.cpp b/src/npair_half_multi2_newton_tri.cpp index 6b28e76789..c741a860f1 100755 --- a/src/npair_half_multi2_newton_tri.cpp +++ b/src/npair_half_multi2_newton_tri.cpp @@ -11,9 +11,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include -#include "npair_half_bytype_newton_tri.h" -#include "neighbor.h" +#include "npair_half_multi2_newton_tri.h" #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" @@ -22,14 +20,11 @@ #include "my_page.h" #include "error.h" -#include "nbin.h" -#include "nstencil.h" - using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NPairHalfBytypeNewtonTri::NPairHalfBytypeNewtonTri(LAMMPS *lmp) : NPair(lmp) {} +NPairHalfMulti2NewtonTri::NPairHalfMulti2NewtonTri(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- KS REWRTIE @@ -39,12 +34,13 @@ NPairHalfBytypeNewtonTri::NPairHalfBytypeNewtonTri(LAMMPS *lmp) : NPair(lmp) {} every pair stored exactly once by some processor ------------------------------------------------------------------------- */ -void NPairHalfBytypeNewtonTri::build(NeighList *list) +void NPairHalfMulti2NewtonTri::build(NeighList *list) { - int i,j,k,n,itype,jtype,ibin,which,ns,imol,iatom,moltemplate; + int i,j,k,n,itype,jtype,ktype,ibin,kbin,which,ns,imol,iatom,moltemplate; tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; + int js; double **x = atom->x; int *type = atom->type; @@ -83,101 +79,97 @@ void NPairHalfBytypeNewtonTri::build(NeighList *list) iatom = molatom[i]; tagprev = tag[i] - iatom - 1; } - - int js; - int kbin; - + // own type: loop over atoms ahead in bin, including ghosts at end of list // if j is owned atom, store by virtue of being ahead of i in list // if j is ghost, store if x[j] "above and to right of" x[i] - ibin = nb->atom2bin_type[itype][i]; + ibin = atom2bin_multi2[itype][i]; - for (int ktype = 1; ktype <= atom->ntypes; ktype++) { + for (ktype = 1; ktype <= atom->ntypes; ktype++) { if (itype == ktype) { - // loop over all atoms in other bins in stencil, store every pair - // skip if i,j neighbor cutoff is less than bin distance - - s = this->ns->stencil_type[itype][itype]; - ns = this->ns->nstencil_type[itype][itype]; - for (k = 0; k < ns; k++) { - js = nb->binhead_type[itype][ibin + s[k]]; - for (j = js; j >= 0; j = nb->bins_type[itype][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) { - if (!moltemplate) - which = find_special(special[i],nspecial[i],tag[j]); - else if (imol >= 0) - which = find_special(onemols[imol]->special[iatom], - onemols[imol]->nspecial[iatom], - tag[j]-tagprev); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (domain->minimum_image_check(delx,dely,delz)) - neighptr[n++] = j; - else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); - } else neighptr[n++] = j; + // loop over all atoms in other bins in stencil, store every pair + // skip if i,j neighbor cutoff is less than bin distance + + s = stencil_multi2[itype][itype]; + ns = nstencil_multi2[itype][itype]; + for (k = 0; k < ns; k++) { + js = binhead_multi2[itype][ibin + s[k]]; + for (j = js; j >= 0; j = bins_multi2[itype][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) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } } - } - } - - } - else { + } else { // smaller -> larger: locate i in the ktype bin structure - kbin = nb->coord2bin(x[i], ktype); - s = this->ns->stencil_type[itype][ktype]; - ns = this->ns->nstencil_type[itype][ktype]; - - for (k = 0; k < ns; k++) { - js = nb->binhead_type[ktype][kbin + s[k]]; - for (j = js; j >= 0; j = nb->bins_type[ktype][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) { - if (!moltemplate) - which = find_special(special[i],nspecial[i],tag[j]); - else if (imol >= 0) - which = find_special(onemols[imol]->special[iatom], - onemols[imol]->nspecial[iatom], - tag[j]-tagprev); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (domain->minimum_image_check(delx,dely,delz)) - neighptr[n++] = j; - else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); - } else neighptr[n++] = j; + kbin = coord2bin(x[i], ktype); + s = stencil_multi2[itype][ktype]; + ns = nstencil_multi2[itype][ktype]; + + for (k = 0; k < ns; k++) { + js = binhead_multi2[ktype][kbin + s[k]]; + for (j = js; j >= 0; j = bins_multi2[ktype][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) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } } - } - } } } + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; diff --git a/src/npair_half_multi2_newton_tri.h b/src/npair_half_multi2_newton_tri.h index 0582801028..0787860c52 100755 --- a/src/npair_half_multi2_newton_tri.h +++ b/src/npair_half_multi2_newton_tri.h @@ -13,23 +13,23 @@ #ifdef NPAIR_CLASS -NPairStyle(half/bytype/newton/tri, - NPairHalfBytypeNewtonTri, - NP_HALF | NP_BYTYPE | NP_NEWTON | NP_TRI) +NPairStyle(half/multi2/newton/tri, + NPairHalfMulti2NewtonTri, + NP_HALF | NP_MULTI2 | NP_NEWTON | NP_TRI) #else -#ifndef LMP_NPAIR_HALF_BYTYPE_NEWTON_TRI_H -#define LMP_NPAIR_HALF_BYTYPE_NEWTON_TRI_H +#ifndef LMP_NPAIR_HALF_MULTI2_NEWTON_TRI_H +#define LMP_NPAIR_HALF_MULTI2_NEWTON_TRI_H #include "npair.h" namespace LAMMPS_NS { -class NPairHalfBytypeNewtonTri : public NPair { +class NPairHalfMulti2NewtonTri : public NPair { public: - NPairHalfBytypeNewtonTri(class LAMMPS *); - ~NPairHalfBytypeNewtonTri() {} + NPairHalfMulti2NewtonTri(class LAMMPS *); + ~NPairHalfMulti2NewtonTri() {} void build(class NeighList *); }; diff --git a/src/npair_half_size_multi2_newtoff.cpp b/src/npair_half_size_multi2_newtoff.cpp index 220ae747a7..0ea9e493ce 100644 --- a/src/npair_half_size_multi2_newtoff.cpp +++ b/src/npair_half_size_multi2_newtoff.cpp @@ -12,20 +12,18 @@ es certain rights in this software. This software is distributed under ------------------------------------------------------------------------- */ #include -#include "npair_half_size_bytype_newtoff.h" +#include "npair_half_size_multi2_newtoff.h" #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" #include "my_page.h" #include "error.h" -#include "nbin.h" -#include "nstencil.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NPairHalfSizeBytypeNewtoff::NPairHalfSizeBytypeNewtoff(LAMMPS *lmp) : NPair(lmp) {} +NPairHalfSizeMulti2Newtoff::NPairHalfSizeMulti2Newtoff(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- REWRITE @@ -36,12 +34,13 @@ NPairHalfSizeBytypeNewtoff::NPairHalfSizeBytypeNewtoff(LAMMPS *lmp) : NPair(lmp) pair stored by me if j is ghost (also stored by proc owning j) ------------------------------------------------------------------------- */ -void NPairHalfSizeBytypeNewtoff::build(NeighList *list) +void NPairHalfSizeMulti2Newtoff::build(NeighList *list) { - int i,j,k,n,itype,jtype,ibin,ns; + int i,j,k,n,itype,jtype,ktype,ibin,kbin,ns; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; + int js; double **x = atom->x; double *radius = atom->radius; @@ -78,41 +77,40 @@ void NPairHalfSizeBytypeNewtoff::build(NeighList *list) // stores own/own pairs only once // stores own/ghost pairs on both procs - int kbin; - ibin = nb->atom2bin_type[itype][i]; - for (int ktype = 1; ktype <= atom->ntypes; ktype++) { + for (ktype = 1; ktype <= atom->ntypes; ktype++) { if (itype == ktype) { - kbin = ibin; + kbin = ibin; + } else { + // Locate i in ktype bin + kbin = coord2bin(x[i], ktype); } - else { - // Locate i in ktype bin - kbin = nb->coord2bin(x[i], ktype); - } - s = this->ns->stencil_type[itype][ktype]; - ns = this->ns->nstencil_type[itype][ktype]; + + s = stencil_multi2[itype][ktype]; + ns = nstencil_multi2[itype][ktype]; for (k = 0; k < ns; k++) { - int js = nb->binhead_type[ktype][kbin + s[k]]; - for (j = js; j >=0; j = nb->bins_type[ktype][j]) { - if (j <= i) continue; - jtype = type[j]; + js = binhead_multi2[ktype][kbin + s[k]]; + for (j = js; j >=0; j = nb->bins_multi2[ktype][j]) { + if (j <= i) continue; + + jtype = type[j]; - if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) 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; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); - 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]; - cutdistsq = (radsum+skin) * (radsum+skin); - - if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; - } - } + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } } } diff --git a/src/npair_half_size_multi2_newtoff.h b/src/npair_half_size_multi2_newtoff.h index 183bd39fb3..15540666c3 100644 --- a/src/npair_half_size_multi2_newtoff.h +++ b/src/npair_half_size_multi2_newtoff.h @@ -13,23 +13,23 @@ #ifdef NPAIR_CLASS -NPairStyle(half/size/bytype/newtoff, - NPairHalfSizeBytypeNewtoff, - NP_HALF | NP_SIZE | NP_BYTYPE | NP_NEWTOFF | NP_ORTHO | NP_TRI) +NPairStyle(half/size/multi2/newtoff, + NPairHalfSizeMulti2Newtoff, + NP_HALF | NP_SIZE | NP_MULTI2 | NP_NEWTOFF | NP_ORTHO | NP_TRI) #else -#ifndef LMP_NPAIR_HALF_SIZE_BYTYPE_NEWTOFF_H -#define LMP_NPAIR_HALF_SIZE_BYTYPE_NEWTOFF_H +#ifndef LMP_NPAIR_HALF_SIZE_MULTI2_NEWTOFF_H +#define LMP_NPAIR_HALF_SIZE_MULTI2_NEWTOFF_H #include "npair.h" namespace LAMMPS_NS { -class NPairHalfSizeBytypeNewtoff : public NPair { +class NPairHalfSizeMulti2Newtoff : public NPair { public: - NPairHalfSizeBytypeNewtoff(class LAMMPS *); - ~NPairHalfSizeBytypeNewtoff() {} + NPairHalfSizeMulti2Newtoff(class LAMMPS *); + ~NPairHalfSizeMulti2Newtoff() {} void build(class NeighList *); }; diff --git a/src/npair_half_size_multi2_newton.cpp b/src/npair_half_size_multi2_newton.cpp index 339c4859a6..3b89c259ec 100755 --- a/src/npair_half_size_multi2_newton.cpp +++ b/src/npair_half_size_multi2_newton.cpp @@ -11,21 +11,18 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include -#include "npair_half_size_bytype_newton.h" +#include "npair_half_size_multi2_newton.h" #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" #include "my_page.h" #include "error.h" -#include "nbin.h" -#include "nstencil.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NPairHalfSizeBytypeNewton::NPairHalfSizeBytypeNewton(LAMMPS *lmp) : NPair(lmp) {} +NPairHalfSizeMulti2Newton::NPairHalfSizeMulti2Newton(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- KS REWRTIE @@ -35,9 +32,9 @@ NPairHalfSizeBytypeNewton::NPairHalfSizeBytypeNewton(LAMMPS *lmp) : NPair(lmp) { every pair stored exactly once by some processor ------------------------------------------------------------------------- */ -void NPairHalfSizeBytypeNewton::build(NeighList *list) +void NPairHalfSizeMulti2Newton::build(NeighList *list) { - int i,j,k,n,itype,jtype,ibin,ns; + int i,j,k,n,itype,jtype,ktype,ibin,kbin,ns,js; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -71,110 +68,104 @@ void NPairHalfSizeBytypeNewton::build(NeighList *list) ztmp = x[i][2]; radi = radius[i]; - int js; - int kbin; - // own type: loop over atoms ahead in bin, including ghosts at end of list // if j is owned atom, store by virtue of being ahead of i in list // if j is ghost, store if x[j] "above and to right of" x[i] - ibin = nb->atom2bin_type[itype][i]; + ibin = nb->atom2bin_multi2[itype][i]; - for (int ktype = 1; ktype <= atom->ntypes; ktype++) { + for (ktype = 1; ktype <= atom->ntypes; ktype++) { if (itype == ktype) { - // own bin ... - js = nb->bins_type[itype][i]; - for (j = js; j >= 0; j = nb->bins_type[itype][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; + // own bin ... + js = bins_multi2[itype][i]; + for (j = js; j >= 0; j = bins_multi2[itype][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; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = 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; - radsum = radi + radius[j]; - cutdistsq = (radsum+skin) * (radsum+skin); - - if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - 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 - - s = this->ns->stencil_type[itype][itype]; - ns = this->ns->nstencil_type[itype][itype]; - for (k = 0; k < ns; k++) { - js = nb->binhead_type[itype][ibin + s[k]]; - for (j = js; j >= 0; j = nb->bins_type[itype][j]) { - jtype = type[j]; - // KS. CHECK ME 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; - radsum = radi + radius[j]; - cutdistsq = (radsum+skin) * (radsum+skin); - - if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - 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 + + s = stencil_multi2[itype][itype]; + ns = nstencil_multi2[itype][itype]; + for (k = 0; k < ns; k++) { + js = binhead_multi2[itype][ibin + s[k]]; + for (j = js; j >= 0; j = bins_multi2[itype][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; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } } - } - } - - } - else { - // KS + } else { // smaller -> larger: locate i in the ktype bin structure - kbin = nb->coord2bin(x[i], ktype); - - s = this->ns->stencil_type[itype][ktype]; - ns = this->ns->nstencil_type[itype][ktype]; - for (k = 0; k < ns; k++) { - js = nb->binhead_type[ktype][kbin + s[k]]; - for (j = js; j >= 0; j = nb->bins_type[ktype][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; - radsum = radi + radius[j]; - cutdistsq = (radsum+skin) * (radsum+skin); - - if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + kbin = coord2bin(x[i], ktype); + + s = stencil_multi2[itype][ktype]; + ns = nstencil_multi2[itype][ktype]; + for (k = 0; k < ns; k++) { + js = binhead_multi2[ktype][kbin + s[k]]; + for (j = js; j >= 0; j = bins_multi2[ktype][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; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } } - } - } } } + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; diff --git a/src/npair_half_size_multi2_newton.h b/src/npair_half_size_multi2_newton.h index 3cd1ee05d1..8c7bc1cc9e 100755 --- a/src/npair_half_size_multi2_newton.h +++ b/src/npair_half_size_multi2_newton.h @@ -13,23 +13,23 @@ #ifdef NPAIR_CLASS -NPairStyle(half/size/bytype/newton, - NPairHalfSizeBytypeNewton, - NP_HALF | NP_SIZE | NP_BYTYPE | NP_NEWTON | NP_ORTHO) +NPairStyle(half/size/multi2/newton, + NPairHalfSizeMulti2Newton, + NP_HALF | NP_SIZE | NP_MULTI2 | NP_NEWTON | NP_ORTHO) #else -#ifndef LMP_NPAIR_HALF_SIZE_BYTYPE_NEWTON_H -#define LMP_NPAIR_HALF_SIZE_BYTYPE_NEWTON_H +#ifndef LMP_NPAIR_HALF_SIZE_MULTI2_NEWTON_H +#define LMP_NPAIR_HALF_SIZE_MULTI2_NEWTON_H #include "npair.h" namespace LAMMPS_NS { -class NPairHalfSizeBytypeNewton : public NPair { +class NPairHalfSizeMulti2Newton : public NPair { public: - NPairHalfSizeBytypeNewton(class LAMMPS *); - ~NPairHalfSizeBytypeNewton() {} + NPairHalfSizeMulti2Newton(class LAMMPS *); + ~NPairHalfSizeMulti2Newton() {} void build(class NeighList *); }; diff --git a/src/npair_half_size_multi2_newton_tri.cpp b/src/npair_half_size_multi2_newton_tri.cpp index 62d499b10f..1c0c1d514f 100755 --- a/src/npair_half_size_multi2_newton_tri.cpp +++ b/src/npair_half_size_multi2_newton_tri.cpp @@ -11,21 +11,18 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include -#include "npair_half_size_bytype_newton_tri.h" +#include "npair_half_size_multi2_newton_tri.h" #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" #include "my_page.h" #include "error.h" -#include "nbin.h" -#include "nstencil.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NPairHalfSizeBytypeNewtonTri::NPairHalfSizeBytypeNewtonTri(LAMMPS *lmp) : NPair(lmp) {} +NPairHalfSizeMulti2NewtonTri::NPairHalfSizeMulti2NewtonTri(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- KS REWRTIE @@ -35,9 +32,9 @@ NPairHalfSizeBytypeNewtonTri::NPairHalfSizeBytypeNewtonTri(LAMMPS *lmp) : NPair( every pair stored exactly once by some processor ------------------------------------------------------------------------- */ -void NPairHalfSizeBytypeNewtonTri::build(NeighList *list) +void NPairHalfSizeMulti2NewtonTri::build(NeighList *list) { - int i,j,k,n,itype,jtype,ibin,ns; + int i,j,k,n,itype,jtype,ktype,ibin,kbin,ns,js; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -71,89 +68,83 @@ void NPairHalfSizeBytypeNewtonTri::build(NeighList *list) ztmp = x[i][2]; radi = radius[i]; - int js; - int kbin; - // own type: loop over atoms ahead in bin, including ghosts at end of list // if j is owned atom, store by virtue of being ahead of i in list // if j is ghost, store if x[j] "above and to right of" x[i] - ibin = nb->atom2bin_type[itype][i]; + ibin = atom2bin_multi2[itype][i]; - for (int ktype = 1; ktype <= atom->ntypes; ktype++) { + for (ktype = 1; ktype <= atom->ntypes; ktype++) { if (itype == ktype) { - // loop over all atoms in other bins in stencil, store every pair - // skip if i,j neighbor cutoff is less than bin distance - - s = this->ns->stencil_type[itype][itype]; - ns = this->ns->nstencil_type[itype][itype]; - for (k = 0; k < ns; k++) { - js = nb->binhead_type[itype][ibin + s[k]]; - for (j = js; j >= 0; j = nb->bins_type[itype][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]; - // KS. CHECK ME 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; - radsum = radi + radius[j]; - cutdistsq = (radsum+skin) * (radsum+skin); - - if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - 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 + + s = stencil_multi2[itype][itype]; + ns = nstencil_multi2[itype][itype]; + for (k = 0; k < ns; k++) { + js = binhead_multi2[itype][ibin + s[k]]; + for (j = js; j >= 0; j = bins_multi2[itype][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; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } } - } - } - - } - else { - // KS + } else { // smaller -> larger: locate i in the ktype bin structure - kbin = nb->coord2bin(x[i], ktype); - - s = this->ns->stencil_type[itype][ktype]; - ns = this->ns->nstencil_type[itype][ktype]; - for (k = 0; k < ns; k++) { - js = nb->binhead_type[ktype][kbin + s[k]]; - for (j = js; j >= 0; j = nb->bins_type[ktype][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; - radsum = radi + radius[j]; - cutdistsq = (radsum+skin) * (radsum+skin); - - if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + kbin = coord2bin(x[i], ktype); + + s = stencil_multi2[itype][ktype]; + ns = nstencil_multi2[itype][ktype]; + for (k = 0; k < ns; k++) { + js = binhead_multi2[ktype][kbin + s[k]]; + for (j = js; j >= 0; j = bins_multi2[ktype][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; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); + + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; + } + } } - } - } } } + ilist[inum++] = i; firstneigh[i] = neighptr; numneigh[i] = n; diff --git a/src/npair_half_size_multi2_newton_tri.h b/src/npair_half_size_multi2_newton_tri.h index 0f770e780d..6ce4f6d2fb 100755 --- a/src/npair_half_size_multi2_newton_tri.h +++ b/src/npair_half_size_multi2_newton_tri.h @@ -13,23 +13,23 @@ #ifdef NPAIR_CLASS -NPairStyle(half/size/bytype/newton/tri, - NPairHalfSizeBytypeNewtonTri, - NP_HALF | NP_SIZE | NP_BYTYPE | NP_NEWTON | NP_TRI) +NPairStyle(half/size/multi2/newton/tri, + NPairHalfSizeMulti2NewtonTri, + NP_HALF | NP_SIZE | NP_MULTI2 | NP_NEWTON | NP_TRI) #else -#ifndef LMP_NPAIR_HALF_SIZE_BYTYPE_NEWTON_TRI_H -#define LMP_NPAIR_HALF_SIZE_BYTYPE_NEWTON_TRI_H +#ifndef LMP_NPAIR_HALF_SIZE_MULTI2_NEWTON_TRI_H +#define LMP_NPAIR_HALF_SIZE_MULTI2_NEWTON_TRI_H #include "npair.h" namespace LAMMPS_NS { -class NPairHalfSizeBytypeNewtonTri : public NPair { +class NPairHalfSizeMulti2NewtonTri : public NPair { public: - NPairHalfSizeBytypeNewtonTri(class LAMMPS *); - ~NPairHalfSizeBytypeNewtonTri() {} + NPairHalfSizeMulti2NewtonTri(class LAMMPS *); + ~NPairHalfSizeMulti2NewtonTri() {} void build(class NeighList *); }; diff --git a/src/nstencil.h b/src/nstencil.h index 7c6a04a56a..faf4ca0872 100644 --- a/src/nstencil.h +++ b/src/nstencil.h @@ -34,7 +34,9 @@ class NStencil : protected Pointers { int *** stencil_multi2; // list of bin offsets in each multi2 stencil int ** maxstencil_multi2; // max stencil size for each multi2 stencil - ^do i need a multi 2 analog to distsq_multi? + // Note distsq_multi is used in multi to quickly skip bins beyond interaction cutoff + // Not quite sure why bins are beyond this distance? Have to think + // Probably not needed for multi2 since bins are more efficiently chosen int sx,sy,sz; // extent of stencil in each dim int **sx_multi2; // analogs for multi tiered