From b1b014aed3e15e2d8ca02b68a0cf567edcff67e8 Mon Sep 17 00:00:00 2001 From: Joel Clemmer Date: Fri, 2 Oct 2020 15:09:47 -0600 Subject: [PATCH] Fixing merge conflicts --- src/comm.cpp | 6 + src/comm.h | 1 + src/comm_brick.cpp | 6 + src/nbin.cpp | 9 + src/nbin.h | 14 + src/nbin_bytype.cpp | 466 ++++++++++++++++++++++++ src/nbin_bytype.h | 56 +++ src/nbin_standard.h | 2 +- src/neighbor.cpp | 17 +- src/neighbor.h | 6 +- src/npair_full_bytype.cpp | 144 ++++++++ src/npair_full_bytype.h | 43 +++ src/npair_half_bytype_newton.cpp | 219 +++++++++++ src/npair_half_bytype_newton.h | 43 +++ src/npair_half_size_bytype_newtoff.cpp | 130 +++++++ src/npair_half_size_bytype_newtoff.h | 43 +++ src/npair_half_size_bytype_newton.cpp | 189 ++++++++++ src/npair_half_size_bytype_newton.h | 43 +++ src/nstencil.h | 4 + src/nstencil_full_bytype_3d.cpp | 192 ++++++++++ src/nstencil_full_bytype_3d.h | 53 +++ src/nstencil_half_bytype_2d_newton.cpp | 192 ++++++++++ src/nstencil_half_bytype_2d_newton.h | 52 +++ src/nstencil_half_bytype_3d_newtoff.cpp | 190 ++++++++++ src/nstencil_half_bytype_3d_newtoff.h | 51 +++ src/nstencil_half_bytype_3d_newton.cpp | 194 ++++++++++ src/nstencil_half_bytype_3d_newton.h | 52 +++ 27 files changed, 2413 insertions(+), 4 deletions(-) create mode 100644 src/nbin_bytype.cpp create mode 100644 src/nbin_bytype.h create mode 100644 src/npair_full_bytype.cpp create mode 100644 src/npair_full_bytype.h create mode 100755 src/npair_half_bytype_newton.cpp create mode 100755 src/npair_half_bytype_newton.h create mode 100644 src/npair_half_size_bytype_newtoff.cpp create mode 100644 src/npair_half_size_bytype_newtoff.h create mode 100755 src/npair_half_size_bytype_newton.cpp create mode 100755 src/npair_half_size_bytype_newton.h create mode 100644 src/nstencil_full_bytype_3d.cpp create mode 100644 src/nstencil_full_bytype_3d.h create mode 100644 src/nstencil_half_bytype_2d_newton.cpp create mode 100644 src/nstencil_half_bytype_2d_newton.h create mode 100644 src/nstencil_half_bytype_3d_newtoff.cpp create mode 100644 src/nstencil_half_bytype_3d_newtoff.h create mode 100644 src/nstencil_half_bytype_3d_newton.cpp create mode 100644 src/nstencil_half_bytype_3d_newton.h diff --git a/src/comm.cpp b/src/comm.cpp index 32a4152294..1efe966459 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -76,6 +76,7 @@ Comm::Comm(LAMMPS *lmp) : Pointers(lmp) grid2proc = nullptr; xsplit = ysplit = zsplit = nullptr; rcbnew = 0; + multi_bytype = 0; // use of OpenMP threads // query OpenMP for number of threads/process set by user at run-time @@ -326,6 +327,11 @@ void Comm::modify_params(int narg, char **arg) for (i=nlo; i<=nhi; ++i) cutusermulti[i] = cut; iarg += 3; + } else if (strcmp(arg[iarg],"cutoff/bytype") == 0) { + if (mode == Comm::SINGLE) + error->all(FLERR,"Use cutoff/bytype in mode multi only"); + multi_bytype = 1; + iarg += 1; } else if (strcmp(arg[iarg],"vel") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal comm_modify command"); if (strcmp(arg[iarg+1],"yes") == 0) ghost_velocity = 1; diff --git a/src/comm.h b/src/comm.h index c8d7add79c..e5fc1bdee4 100644 --- a/src/comm.h +++ b/src/comm.h @@ -161,6 +161,7 @@ class Comm : protected Pointers { int (*)(int, char *, int &, int *&, char *&, void *), int, char *&, int, void *, int); void rendezvous_stats(int, int, int, int, int, int, bigint); + int multi_bytype; // 1 if multi cutoff is intra-type cutoff public: enum{MULTIPLE}; diff --git a/src/comm_brick.cpp b/src/comm_brick.cpp index 374c394b0b..57703b578c 100644 --- a/src/comm_brick.cpp +++ b/src/comm_brick.cpp @@ -181,6 +181,12 @@ void CommBrick::setup() cutghostmulti[i][0] = MAX(cut,cuttype[i]); cutghostmulti[i][1] = MAX(cut,cuttype[i]); cutghostmulti[i][2] = MAX(cut,cuttype[i]); + if (multi_bytype == 1) { + // Set the BYTYPE cutoff + cutghostmulti[i][0] = sqrt(neighbor->cutneighsq[i][i]); + cutghostmulti[i][1] = sqrt(neighbor->cutneighsq[i][i]); + cutghostmulti[i][2] = sqrt(neighbor->cutneighsq[i][i]); + } } } diff --git a/src/nbin.cpp b/src/nbin.cpp index 1fe011e9f2..64ae7e7b83 100644 --- a/src/nbin.cpp +++ b/src/nbin.cpp @@ -150,6 +150,15 @@ int NBin::coord2bin(double *x) return (iz-mbinzlo)*mbiny*mbinx + (iy-mbinylo)*mbinx + (ix-mbinxlo); } +/* ---------------------------------------------------------------------- + to be overridden by NBinType + ------------------------------------------------------------------------- */ + +int NBin::coord2bin(double * x, int itype) +{ + error->all(FLERR,"coord2bin(x, itype) not available.\n"); + return -1; +} /* ---------------------------------------------------------------------- */ diff --git a/src/nbin.h b/src/nbin.h index 4bfe579514..96cc7eac64 100644 --- a/src/nbin.h +++ b/src/nbin.h @@ -37,6 +37,18 @@ class NBin : protected Pointers { double cutoff_custom; // cutoff set by requestor + // Analogues for NBinType + int * nbinx_type, * nbiny_type, * nbinz_type; + int * mbins_type; + int * mbinx_type, * mbiny_type, * mbinz_type; + int * mbinxlo_type, * mbinylo_type, * mbinzlo_type; + double * binsizex_type, * binsizey_type, * binsizez_type; + double * bininvx_type, * bininvy_type, * bininvz_type; + + int ** binhead_type; + int ** bins_type; + int ** atom2bin_type; + NBin(class LAMMPS *); ~NBin(); void post_constructor(class NeighRequest *); @@ -50,6 +62,8 @@ class NBin : protected Pointers { // Kokkos package int kokkos; // 1 if class stores Kokkos data + // For NBinType + virtual int coord2bin(double *, int); protected: diff --git a/src/nbin_bytype.cpp b/src/nbin_bytype.cpp new file mode 100644 index 0000000000..f719ed116f --- /dev/null +++ b/src/nbin_bytype.cpp @@ -0,0 +1,466 @@ +/* ---------------------------------------------------------------------- + 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 "nbin_bytype.h" +#include "neighbor.h" +#include "atom.h" +#include "group.h" +#include "domain.h" +#include "comm.h" +#include "update.h" +#include "error.h" + +#include +#include "memory.h" + +using namespace LAMMPS_NS; + +#define SMALL 1.0e-6 // Duplicated from NBinStandard +#define CUT2BIN_RATIO 100 // Duplicated from NBinStandard + +/* ---------------------------------------------------------------------- */ + +NBinBytype::NBinBytype(LAMMPS *lmp) : NBin(lmp) { + + nbinx_type = NULL; nbiny_type = NULL; nbinz_type = NULL; + mbins_type = NULL; + mbinx_type = NULL; mbiny_type = NULL, mbinz_type = NULL; + mbinxlo_type = NULL; + mbinylo_type = NULL; + mbinzlo_type = NULL; + binsizex_type = NULL; binsizey_type = NULL; binsizez_type = NULL; + bininvx_type = NULL; bininvy_type = NULL; bininvz_type = NULL; + binhead_type = NULL; + bins_type = NULL; + atom2bin_type = NULL; + maxtypes = 0; + maxbins_type = NULL; +} + +NBinBytype::~NBinBytype() { + + memory->destroy(nbinx_type); + memory->destroy(nbiny_type); + memory->destroy(nbinz_type); + memory->destroy(mbins_type); + memory->destroy(mbinx_type); + memory->destroy(mbiny_type); + memory->destroy(mbinz_type); + memory->destroy(mbinxlo_type); + memory->destroy(mbinylo_type); + memory->destroy(mbinzlo_type); + + memory->destroy(binsizex_type); + memory->destroy(binsizey_type); + memory->destroy(binsizez_type); + memory->destroy(bininvx_type); + memory->destroy(bininvy_type); + memory->destroy(bininvz_type); + + for (int n = 1; n <= maxtypes; n++) { + memory->destroy(binhead_type[n]); + memory->destroy(bins_type[n]); + memory->destroy(atom2bin_type[n]); + } + delete [] binhead_type; + delete [] bins_type; + delete [] atom2bin_type; + + memory->destroy(maxbins_type); +} + +/* ---------------------------------------------------------------------- + arrange storage for types + allows ntypes to change, but not currently expected after initialisation + ------------------------------------------------------------------------- */ + +void NBinBytype::setup_types() { + + int n; + + if (atom->ntypes > maxtypes) { + + // Clear any/all memory for existing types + + for (n = 1; n <= maxtypes; n++) { + memory->destroy(atom2bin_type[n]); + memory->destroy(binhead_type[n]); + memory->destroy(bins_type[n]); + } + delete [] atom2bin_type; + delete [] binhead_type; + delete [] bins_type; + + // Realloacte at updated maxtypes + + maxtypes = atom->ntypes; + + atom2bin_type = new int*[maxtypes+1](); + binhead_type = new int*[maxtypes+1](); + bins_type = new int*[maxtypes+1](); + + memory->destroy(nbinx_type); + memory->destroy(nbiny_type); + memory->destroy(nbinz_type); + memory->create(nbinx_type, maxtypes+1, "nBinType:nbinx_type"); + memory->create(nbiny_type, maxtypes+1, "nBinType:nbiny_type"); + memory->create(nbinz_type, maxtypes+1, "nBinType:nbinz_type"); + + memory->destroy(mbins_type); + memory->destroy(mbinx_type); + memory->destroy(mbiny_type); + memory->destroy(mbinz_type); + memory->create(mbins_type, maxtypes+1, "nBinType:mbins_type"); + memory->create(mbinx_type, maxtypes+1, "nBinType:mbinx_type"); + memory->create(mbiny_type, maxtypes+1, "nBinType:mbiny_type"); + memory->create(mbinz_type, maxtypes+1, "nBinType:mbinz_type"); + + memory->destroy(mbinxlo_type); + memory->destroy(mbinylo_type); + memory->destroy(mbinzlo_type); + memory->create(mbinxlo_type, maxtypes+1, "nBinType:mbinxlo_type"); + memory->create(mbinylo_type, maxtypes+1, "nBinType:mbinylo_type"); + memory->create(mbinzlo_type, maxtypes+1, "nBinType:mbinzlo_type"); + + memory->destroy(binsizex_type); + memory->destroy(binsizey_type); + memory->destroy(binsizez_type); + memory->create(binsizex_type, maxtypes+1, "nBinType:binsizex_type"); + memory->create(binsizey_type, maxtypes+1, "nBinType:binsizey_type"); + memory->create(binsizez_type, maxtypes+1, "nBinType:binsizez_type"); + + memory->destroy(bininvx_type); + memory->destroy(bininvy_type); + memory->destroy(bininvz_type); + memory->create(bininvx_type, maxtypes+1, "nBinType:bininvx_type"); + memory->create(bininvy_type, maxtypes+1, "nBinType:bininvy_type"); + memory->create(bininvz_type, maxtypes+1, "nBinType:bininvz_type"); + + memory->destroy(maxbins_type); + memory->create(maxbins_type, maxtypes+1, "nBinType:maxbins_type"); + // make sure reallocation occurs in bin_atoms_setup() + for (n = 1; n <= maxtypes; n++) { + maxbins_type[n] = 0; + } + maxatom = 0; + } + +} + +/* --------------------------------------------------------------------- + identify index of type with smallest cutoff + --------------------------------------------------------------------- */ + +int NBinBytype::itype_min() { + + int itypemin = 1; + double ** cutneighsq; + + cutneighsq = neighbor->cutneighsq; + + for (int n = 1; n <= atom->ntypes; n++) { + if (cutneighsq[n][n] < cutneighsq[itypemin][itypemin]) { + itypemin = n; + } + } + + return itypemin; +} + +/* ---------------------------------------------------------------------- + setup neighbor binning geometry + ---------------------------------------------------------------------- */ + +void NBinBytype::setup_bins(int style) +{ + int n; + int itypemin; + + setup_types(); + itypemin = itype_min(); + + // bbox = size of bbox of entire domain + // bsubbox lo/hi = bounding box of my subdomain extended by comm->cutghost + // for triclinic: + // bbox bounds all 8 corners of tilted box + // subdomain is in lamda coords + // include dimension-dependent extension via comm->cutghost + // domain->bbox() converts lamda extent to box coords and computes bbox + + double bbox[3],bsubboxlo[3],bsubboxhi[3]; + double *cutghost = comm->cutghost; + + if (triclinic == 0) { + bsubboxlo[0] = domain->sublo[0] - cutghost[0]; + bsubboxlo[1] = domain->sublo[1] - cutghost[1]; + bsubboxlo[2] = domain->sublo[2] - cutghost[2]; + bsubboxhi[0] = domain->subhi[0] + cutghost[0]; + bsubboxhi[1] = domain->subhi[1] + cutghost[1]; + bsubboxhi[2] = domain->subhi[2] + cutghost[2]; + } else { + double lo[3],hi[3]; + lo[0] = domain->sublo_lamda[0] - cutghost[0]; + lo[1] = domain->sublo_lamda[1] - cutghost[1]; + lo[2] = domain->sublo_lamda[2] - cutghost[2]; + hi[0] = domain->subhi_lamda[0] + cutghost[0]; + hi[1] = domain->subhi_lamda[1] + cutghost[1]; + hi[2] = domain->subhi_lamda[2] + cutghost[2]; + domain->bbox(lo,hi,bsubboxlo,bsubboxhi); + } + + bbox[0] = bboxhi[0] - bboxlo[0]; + bbox[1] = bboxhi[1] - bboxlo[1]; + bbox[2] = bboxhi[2] - bboxlo[2]; + + // For each type... + + for (n = 1; n <= atom->ntypes; n++) { + // binsize_user only relates to smallest type + // optimal bin size is roughly 1/2 the type-type cutoff + // special case of all cutoffs = 0.0, binsize = box size + + double binsize_optimal; + if (n == itypemin && binsizeflag) binsize_optimal = binsize_user; + else binsize_optimal = 0.5*sqrt(neighbor->cutneighsq[n][n]); + if (binsize_optimal == 0.0) binsize_optimal = bbox[0]; + double binsizeinv = 1.0/binsize_optimal; + + // test for too many global bins in any dimension due to huge global domain + + if (bbox[0]*binsizeinv > MAXSMALLINT || bbox[1]*binsizeinv > MAXSMALLINT || + bbox[2]*binsizeinv > MAXSMALLINT) + error->all(FLERR,"Domain too large for neighbor bins"); + + // create actual bins + // always have one bin even if cutoff > bbox + // for 2d, nbinz_type[n] = 1 + + nbinx_type[n] = static_cast (bbox[0]*binsizeinv); + nbiny_type[n] = static_cast (bbox[1]*binsizeinv); + if (dimension == 3) nbinz_type[n] = static_cast (bbox[2]*binsizeinv); + else nbinz_type[n] = 1; + + if (nbinx_type[n] == 0) nbinx_type[n] = 1; + if (nbiny_type[n] == 0) nbiny_type[n] = 1; + if (nbinz_type[n] == 0) nbinz_type[n] = 1; + + // compute actual bin size for nbins to fit into box exactly + // error if actual bin size << cutoff, since will create a zillion bins + // this happens when nbin = 1 and box size << cutoff + // typically due to non-periodic, flat system in a particular dim + // in that extreme case, should use NSQ not BIN neighbor style + + binsizex_type[n] = bbox[0]/nbinx_type[n]; + binsizey_type[n] = bbox[1]/nbiny_type[n]; + binsizez_type[n] = bbox[2]/nbinz_type[n]; + + bininvx_type[n] = 1.0 / binsizex_type[n]; + bininvy_type[n] = 1.0 / binsizey_type[n]; + bininvz_type[n] = 1.0 / binsizez_type[n]; + + if (binsize_optimal*bininvx_type[n] > CUT2BIN_RATIO || + binsize_optimal*bininvy_type[n] > CUT2BIN_RATIO || + binsize_optimal*bininvz_type[n] > CUT2BIN_RATIO) + error->all(FLERR,"Cannot use neighbor bins - box size << cutoff"); + + // mbinlo/hi = lowest and highest global bins my ghost atoms could be in + // coord = lowest and highest values of coords for my ghost atoms + // static_cast(-1.5) = -1, so subract additional -1 + // add in SMALL for round-off safety + + int mbinxhi,mbinyhi,mbinzhi; + double coord; + + coord = bsubboxlo[0] - SMALL*bbox[0]; + mbinxlo_type[n] = static_cast ((coord-bboxlo[0])*bininvx_type[n]); + if (coord < bboxlo[0]) mbinxlo_type[n] = mbinxlo_type[n] - 1; + coord = bsubboxhi[0] + SMALL*bbox[0]; + mbinxhi = static_cast ((coord-bboxlo[0])*bininvx_type[n]); + + coord = bsubboxlo[1] - SMALL*bbox[1]; + mbinylo_type[n] = static_cast ((coord-bboxlo[1])*bininvy_type[n]); + if (coord < bboxlo[1]) mbinylo_type[n] = mbinylo_type[n] - 1; + coord = bsubboxhi[1] + SMALL*bbox[1]; + mbinyhi = static_cast ((coord-bboxlo[1])*bininvy_type[n]); + + if (dimension == 3) { + coord = bsubboxlo[2] - SMALL*bbox[2]; + mbinzlo_type[n] = static_cast ((coord-bboxlo[2])*bininvz_type[n]); + if (coord < bboxlo[2]) mbinzlo_type[n] = mbinzlo_type[n] - 1; + coord = bsubboxhi[2] + SMALL*bbox[2]; + mbinzhi = static_cast ((coord-bboxlo[2])*bininvz_type[n]); + } + + // extend bins by 1 to insure stencil extent is included + // for 2d, only 1 bin in z + + mbinxlo_type[n] = mbinxlo_type[n] - 1; + mbinxhi = mbinxhi + 1; + mbinx_type[n] = mbinxhi - mbinxlo_type[n] + 1; + + mbinylo_type[n] = mbinylo_type[n] - 1; + mbinyhi = mbinyhi + 1; + mbiny_type[n] = mbinyhi - mbinylo_type[n] + 1; + + if (dimension == 3) { + mbinzlo_type[n] = mbinzlo_type[n] - 1; + mbinzhi = mbinzhi + 1; + } else mbinzlo_type[n] = mbinzhi = 0; + mbinz_type[n] = mbinzhi - mbinzlo_type[n] + 1; + + bigint bbin = ((bigint) mbinx_type[n]) + * ((bigint) mbiny_type[n]) * ((bigint) mbinz_type[n]) + 1; + if (bbin > MAXSMALLINT) error->one(FLERR,"Too many neighbor bins"); + mbins_type[n] = bbin; + } + +} + +/* ---------------------------------------------------------------------- + bin owned and ghost atoms by type +------------------------------------------------------------------------- */ + +void NBinBytype::bin_atoms() +{ + int i,ibin,n; + + last_bin = update->ntimestep; + for (n = 1; n <= maxtypes; n++) { + for (i = 0; i < mbins_type[n]; i++) binhead_type[n][i] = -1; + } + + // bin in reverse order so linked list will be in forward order + // also puts ghost atoms at end of list, which is necessary + + double **x = atom->x; + int *mask = atom->mask; + int *type = atom->type; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + + if (includegroup) { + int bitmask = group->bitmask[includegroup]; + for (i = nall-1; i >= nlocal; i--) { + if (mask[i] & bitmask) { + n = type[i]; + ibin = coord2bin(x[i], n); + atom2bin_type[n][i] = ibin; + bins_type[n][i] = binhead_type[n][ibin]; + binhead_type[n][ibin] = i; + } + } + for (i = atom->nfirst-1; i >= 0; i--) { + n = type[i]; + ibin = coord2bin(x[i], n); + atom2bin_type[n][i] = ibin; + bins_type[n][i] = binhead_type[n][ibin]; + binhead_type[n][ibin] = i; + } + } else { + for (i = nall-1; i >= 0; i--) { + n = type[i]; + ibin = coord2bin(x[i], n); + atom2bin_type[n][i] = ibin; + bins_type[n][i] = binhead_type[n][ibin]; + binhead_type[n][ibin] = i; + } + } +} + + +/* ---------------------------------------------------------------------- + allocate/reallocate vectors +------------------------------------------------------------------------- */ + +void NBinBytype::bin_atoms_setup(int nall) { + + // all atom2bin, bins must be of length nall + if (nall > maxatom) { + for (int n = 1; n <= maxtypes; n++) { + memory->destroy(bins_type[n]); + memory->destroy(atom2bin_type[n]); + memory->create(bins_type[n], nall, "NBinBytype:bin_type"); + memory->create(atom2bin_type[n], nall, "NBinBytype:atom2bin_type"); + } + maxatom = nall; + } + + for (int n = 1; n <= maxtypes; n++) { + if (mbins_type[n] > maxbins_type[n]) { + maxbins_type[n] = mbins_type[n]; + memory->destroy(binhead_type[n]); + memory->create(binhead_type[n], mbins_type[n], "NBinBytype:mbins_type"); + } + } + +} + + +/* ---------------------------------------------------------------------- + convert atom coords into local bin # for bin type it +------------------------------------------------------------------------- */ + +int NBinBytype::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_type[it]) + nbinx_type[it]; + else if (x[0] >= bboxlo[0]) { + ix = static_cast ((x[0]-bboxlo[0])*bininvx_type[it]); + ix = MIN(ix,nbinx_type[it]-1); + } else + ix = static_cast ((x[0]-bboxlo[0])*bininvx_type[it]) - 1; + + if (x[1] >= bboxhi[1]) + iy = static_cast ((x[1]-bboxhi[1])*bininvy_type[it]) + nbiny_type[it]; + else if (x[1] >= bboxlo[1]) { + iy = static_cast ((x[1]-bboxlo[1])*bininvy_type[it]); + iy = MIN(iy,nbiny_type[it]-1); + } else + iy = static_cast ((x[1]-bboxlo[1])*bininvy_type[it]) - 1; + + if (x[2] >= bboxhi[2]) + iz = static_cast ((x[2]-bboxhi[2])*bininvz_type[it]) + nbinz_type[it]; + else if (x[2] >= bboxlo[2]) { + iz = static_cast ((x[2]-bboxlo[2])*bininvz_type[it]); + iz = MIN(iz,nbinz_type[it]-1); + } else + iz = static_cast ((x[2]-bboxlo[2])*bininvz_type[it]) - 1; + + + ibin = (iz-mbinzlo_type[it])*mbiny_type[it]*mbinx_type[it] + + (iy-mbinylo_type[it])*mbinx_type[it] + + (ix-mbinxlo_type[it]); + return ibin; +} + + +/* ---------------------------------------------------------------------- + tot up for types + ---------------------------------------------------------------------- */ + +bigint NBinBytype::memory_usage() +{ + bigint bytes = 0; + + for (int m = 1; m < maxtypes; m++) { + bytes += 2*maxatom*sizeof(int); + bytes += maxbins_type[m]*sizeof(int); + } + return bytes; +} diff --git a/src/nbin_bytype.h b/src/nbin_bytype.h new file mode 100644 index 0000000000..89e6ca3650 --- /dev/null +++ b/src/nbin_bytype.h @@ -0,0 +1,56 @@ +/* -*- 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 NBIN_CLASS + +NBinStyle(bytype, + NBinBytype, + NB_BYTYPE) + +#else + +#ifndef LMP_NBIN_BYTYPE_H +#define LMP_NBIN_BYTYPE_H + +#include "nbin.h" + +namespace LAMMPS_NS { + +class NBinBytype : public NBin { + public: + + NBinBytype(class LAMMPS *); + ~NBinBytype(); + void bin_atoms_setup(int); + void setup_bins(int); + void bin_atoms(); + + int coord2bin(double *x, int itype); + bigint memory_usage(); + + private: + int maxtypes; + int * maxbins_type; + + void setup_types(); + int itype_min(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/nbin_standard.h b/src/nbin_standard.h index 14b587c24c..06defdb101 100644 --- a/src/nbin_standard.h +++ b/src/nbin_standard.h @@ -15,7 +15,7 @@ NBinStyle(standard, NBinStandard, - 0) + NB_STANDARD) #else diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 58737b815d..6e861b5a2b 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -1620,6 +1620,14 @@ int Neighbor::choose_bin(NeighRequest *rq) if (!rq->kokkos_device != !(mask & NB_KOKKOS_DEVICE)) continue; if (!rq->kokkos_host != !(mask & NB_KOKKOS_HOST)) continue; + // neighbor style is BIN or MULTI or BYTYPE and must match + + if (style == Neighbor::BIN || style == Neighbor::MULTI) { + if (!(mask & NB_STANDARD)) continue; + } else if (style == Neighbor::BYTYPE) { + if (!(mask & NB_BYTYPE)) continue; + } + return i+1; } @@ -1690,12 +1698,14 @@ int Neighbor::choose_stencil(NeighRequest *rq) if (!rq->ghost != !(mask & NS_GHOST)) continue; if (!rq->ssa != !(mask & NS_SSA)) continue; - // neighbor style is BIN or MULTI and must match + // neighbor style is one of BIN, MULTI or BYTYPE and must match if (style == Neighbor::BIN) { if (!(mask & NS_BIN)) continue; } else if (style == Neighbor::MULTI) { if (!(mask & NS_MULTI)) continue; + } else if (style == Neighbor::BYTYPE) { + if (!(mask & NS_BYTYPE)) continue; } // dimension is 2 or 3 and must match @@ -1825,7 +1835,7 @@ int Neighbor::choose_pair(NeighRequest *rq) if (!rq->halffull != !(mask & NP_HALF_FULL)) continue; if (!rq->off2on != !(mask & NP_OFF2ON)) continue; - // neighbor style is one of NSQ,BIN,MULTI and must match + // neighbor style is one of NSQ,BIN,MULTI or BYTYPE and must match if (style == Neighbor::NSQ) { if (!(mask & NP_NSQ)) continue; @@ -1833,6 +1843,8 @@ int Neighbor::choose_pair(NeighRequest *rq) if (!(mask & NP_BIN)) continue; } else if (style == Neighbor::MULTI) { if (!(mask & NP_MULTI)) continue; + } else if (style == Neighbor::BYTYPE) { + if (!(mask & NP_BYTYPE)) continue; } // domain triclinic flag is on or off and must match @@ -2199,6 +2211,7 @@ void Neighbor::set(int narg, char **arg) if (strcmp(arg[1],"nsq") == 0) style = Neighbor::NSQ; else if (strcmp(arg[1],"bin") == 0) style = Neighbor::BIN; else if (strcmp(arg[1],"multi") == 0) style = Neighbor::MULTI; + else if (strcmp(arg[1],"bytype") == 0) style = Neighbor::BYTYPE; else error->all(FLERR,"Illegal neighbor command"); if (style == Neighbor::MULTI && lmp->citeme) lmp->citeme->add(cite_neigh_multi); diff --git a/src/neighbor.h b/src/neighbor.h index 9ee2af9c75..e5b76dd41b 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -20,7 +20,7 @@ namespace LAMMPS_NS { class Neighbor : protected Pointers { public: - enum{NSQ,BIN,MULTI}; + enum{NSQ,BIN,MULTI,BYTYPE}; int style; // 0,1,2 = nsq, bin, multi int every; // build every this many steps int delay; // delay build for this many steps @@ -239,6 +239,8 @@ namespace NeighConst { static const int NB_KOKKOS_DEVICE = 1<<1; static const int NB_KOKKOS_HOST = 1<<2; static const int NB_SSA = 1<<3; + static const int NB_BYTYPE = 1<<4; + static const int NB_STANDARD = 1<<5; static const int NS_BIN = 1<<0; static const int NS_MULTI = 1<<1; @@ -252,6 +254,7 @@ namespace NeighConst { static const int NS_TRI = 1<<9; static const int NS_GHOST = 1<<10; static const int NS_SSA = 1<<11; + static const int NS_BYTYPE = 1<<12; static const int NP_NSQ = 1<<0; static const int NP_BIN = 1<<1; @@ -278,6 +281,7 @@ namespace NeighConst { static const int NP_SKIP = 1<<22; static const int NP_HALF_FULL = 1<<23; static const int NP_OFF2ON = 1<<24; + static const int NP_BYTYPE = 1<<25; } } diff --git a/src/npair_full_bytype.cpp b/src/npair_full_bytype.cpp new file mode 100644 index 0000000000..d244c263ed --- /dev/null +++ b/src/npair_full_bytype.cpp @@ -0,0 +1,144 @@ +/* ---------------------------------------------------------------------- + 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_full_bytype.h" +#include "neighbor.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" + +#include "nbin.h" +#include "nstencil.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NPairFullBytype::NPairFullBytype(LAMMPS *lmp) : NPair(lmp) {} + +/* ---------------------------------------------------------------------- + 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) +{ + int i,j,k,n,itype,jtype,ibin,which,ns,imol,iatom,moltemplate; + tagint tagprev; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr,*s; + + 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 == 2) 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 + // 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++) { + if (itype == ktype) { + kbin = ibin; + } + 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]; + 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]) { + jtype = type[j]; + 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) { + 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; + list->gnum = 0; +} diff --git a/src/npair_full_bytype.h b/src/npair_full_bytype.h new file mode 100644 index 0000000000..9ce221ff22 --- /dev/null +++ b/src/npair_full_bytype.h @@ -0,0 +1,43 @@ +/* -*- 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(full/bytype, + NPairFullBytype, + NP_FULL | NP_BYTYPE | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI) + +#else + +#ifndef LMP_NPAIR_FULL_BYTYPE_H +#define LMP_NPAIR_FULL_BYTYPE_H + +#include "npair.h" + +namespace LAMMPS_NS { + +class NPairFullBytype : public NPair { + public: + NPairFullBytype(class LAMMPS *); + ~NPairFullBytype() {} + void build(class NeighList *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/npair_half_bytype_newton.cpp b/src/npair_half_bytype_newton.cpp new file mode 100755 index 0000000000..66899564be --- /dev/null +++ b/src/npair_half_bytype_newton.cpp @@ -0,0 +1,219 @@ +/* ---------------------------------------------------------------------- + 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 +#include "npair_half_bytype_newton.h" +#include "neighbor.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" + +#include "nbin.h" +#include "nstencil.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NPairHalfBytypeNewton::NPairHalfBytypeNewton(LAMMPS *lmp) : NPair(lmp) {} + +/* ---------------------------------------------------------------------- + KS REWRTIE + 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 NPairHalfBytypeNewton::build(NeighList *list) +{ + int i,j,k,n,itype,jtype,ibin,which,ns,imol,iatom,moltemplate; + tagint tagprev; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr,*s; + + 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 == 2) 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; + } + + 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++) { + + 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; + } + } + + 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; + } + } + } + + } + 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; + } + } + } + } + } + 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_bytype_newton.h b/src/npair_half_bytype_newton.h new file mode 100755 index 0000000000..8be7292219 --- /dev/null +++ b/src/npair_half_bytype_newton.h @@ -0,0 +1,43 @@ +/* -*- 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/bytype/newton, + NPairHalfBytypeNewton, + NP_HALF | NP_BYTYPE | NP_NEWTON | NP_ORTHO) + +#else + +#ifndef LMP_NPAIR_HALF_BYTYPE_NEWTON_H +#define LMP_NPAIR_HALF_BYTYPE_NEWTON_H + +#include "npair.h" + +namespace LAMMPS_NS { + +class NPairHalfBytypeNewton : public NPair { + public: + NPairHalfBytypeNewton(class LAMMPS *); + ~NPairHalfBytypeNewton() {} + void build(class NeighList *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/npair_half_size_bytype_newtoff.cpp b/src/npair_half_size_bytype_newtoff.cpp new file mode 100644 index 0000000000..f065c60736 --- /dev/null +++ b/src/npair_half_size_bytype_newtoff.cpp @@ -0,0 +1,130 @@ +/* ---------------------------------------------------------------------- + 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 +es 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 +#include "npair_half_size_bytype_newtoff.h" +#include "neighbor.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) {} + +/* ---------------------------------------------------------------------- + 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 NPairHalfSizeBytypeNewtoff::build(NeighList *list) +{ + int i,j,k,n,itype,jtype,ibin,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutdistsq; + int *neighptr,*s; + + double **x = atom->x; + double *radius = atom->radius; + int *type = atom->type; + int *mask = atom->mask; + tagint *molecule = atom->molecule; + int nlocal = atom->nlocal; + if (includegroup) nlocal = atom->nfirst; + + int history = list->history; + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + MyPage *ipage = list->ipage; + + int mask_history = 3 << SBBITS; + + 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]; + radi = radius[i]; + + // 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 + + int kbin; + + ibin = nb->atom2bin_type[itype][i]; + for (int ktype = 1; ktype <= atom->ntypes; ktype++) { + if (itype == ktype) { + kbin = ibin; + } + 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]; + 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]; + + 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; + 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_size_bytype_newtoff.h b/src/npair_half_size_bytype_newtoff.h new file mode 100644 index 0000000000..e131a12f4b --- /dev/null +++ b/src/npair_half_size_bytype_newtoff.h @@ -0,0 +1,43 @@ +/* -*- 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/size/bytype/newtoff, + NPairHalfSizeBytypeNewtoff, + NP_HALF | NP_SIZE | NP_BYTYPE | NP_NEWTOFF | NP_ORTHO | NP_TRI) + +#else + +#ifndef LMP_NPAIR_HALF_SIZE_BYTYPE_NEWTOFF_H +#define LMP_NPAIR_HALF_SIZE_BYTYPE_NEWTOFF_H + +#include "npair.h" + +namespace LAMMPS_NS { + +class NPairHalfSizeBytypeNewtoff : public NPair { + public: + NPairHalfSizeBytypeNewtoff(class LAMMPS *); + ~NPairHalfSizeBytypeNewtoff() {} + void build(class NeighList *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/npair_half_size_bytype_newton.cpp b/src/npair_half_size_bytype_newton.cpp new file mode 100755 index 0000000000..3981635439 --- /dev/null +++ b/src/npair_half_size_bytype_newton.cpp @@ -0,0 +1,189 @@ +/* ---------------------------------------------------------------------- + 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 +#include "npair_half_size_bytype_newton.h" +#include "neighbor.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) {} + +/* ---------------------------------------------------------------------- + KS REWRTIE + 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 NPairHalfSizeBytypeNewton::build(NeighList *list) +{ + int i,j,k,n,itype,jtype,ibin,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + double radi,radsum,cutdistsq; + int *neighptr,*s; + + double **x = atom->x; + double *radius = atom->radius; + int *type = atom->type; + int *mask = atom->mask; + tagint *molecule = atom->molecule; + int nlocal = atom->nlocal; + if (includegroup) nlocal = atom->nfirst; + + int history = list->history; + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + MyPage *ipage = list->ipage; + + int mask_history = 3 << SBBITS; + + 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]; + 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]; + + for (int 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; + } + } + + 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; + } + } + } + + } + else { + // KS + // 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; + } + } + } + } + } + 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_size_bytype_newton.h b/src/npair_half_size_bytype_newton.h new file mode 100755 index 0000000000..a6bb7a00a8 --- /dev/null +++ b/src/npair_half_size_bytype_newton.h @@ -0,0 +1,43 @@ +/* -*- 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/size/bytype/newton, + NPairHalfSizeBytypeNewton, + NP_HALF | NP_SIZE | NP_BYTYPE | NP_NEWTON | NP_ORTHO) + +#else + +#ifndef LMP_NPAIR_HALF_SIZE_BYTYPE_NEWTON_H +#define LMP_NPAIR_HALF_SIZE_BYTYPE_NEWTON_H + +#include "npair.h" + +namespace LAMMPS_NS { + +class NPairHalfSizeBytypeNewton : public NPair { + public: + NPairHalfSizeBytypeNewton(class LAMMPS *); + ~NPairHalfSizeBytypeNewton() {} + void build(class NeighList *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/nstencil.h b/src/nstencil.h index 75e678b483..bdd656b937 100644 --- a/src/nstencil.h +++ b/src/nstencil.h @@ -34,6 +34,10 @@ class NStencil : protected Pointers { double cutoff_custom; // cutoff set by requestor + // BYTYPE stencils + int ** nstencil_type; // No. of bins for stencil[itype][jtype] + int *** stencil_type; // Stencil for [itype][jtype] + NStencil(class LAMMPS *); virtual ~NStencil(); void post_constructor(class NeighRequest *); diff --git a/src/nstencil_full_bytype_3d.cpp b/src/nstencil_full_bytype_3d.cpp new file mode 100644 index 0000000000..3936ec2483 --- /dev/null +++ b/src/nstencil_full_bytype_3d.cpp @@ -0,0 +1,192 @@ +/* ---------------------------------------------------------------------- + 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 "nstencil_full_bytype_3d.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "nbin.h" +#include "memory.h" +#include "atom.h" + +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NStencilFullBytype3d::NStencilFullBytype3d(LAMMPS *lmp) : + NStencil(lmp) +{ + maxstencil_type = NULL; +} + +NStencilFullBytype3d::~NStencilFullBytype3d() { + + if (maxstencil_type) { + memory->destroy(nstencil_type); + for (int i = 1; i <= atom->ntypes; i++) { + for (int j = 0; j <= atom->ntypes; j++) { + if (maxstencil_type[i][j] > 0) memory->destroy(stencil_type[i][j]); + } + delete [] stencil_type[i]; + } + delete [] stencil_type; + memory->destroy(maxstencil_type); + } +} + +void NStencilFullBytype3d::copy_bin_info_bytype(int itype) { + + mbinx = nb->mbinx_type[itype]; + mbiny = nb->mbiny_type[itype]; + mbinz = nb->mbinz_type[itype]; + binsizex = nb->binsizex_type[itype]; + binsizey = nb->binsizey_type[itype]; + binsizez = nb->binsizez_type[itype]; + bininvx = nb->bininvx_type[itype]; + bininvy = nb->bininvy_type[itype]; + bininvz = nb->bininvz_type[itype]; +} + +int NStencilFullBytype3d::copy_neigh_info_bytype(int itype) { + + cutneighmaxsq = neighbor->cutneighsq[itype][itype]; + cutneighmax = sqrt(cutneighmaxsq); + cuttypesq = neighbor->cuttypesq; + + // sx,sy,sz = max range of stencil in each dim + // smax = max possible size of entire 3d stencil + // stencil will be empty if cutneighmax = 0.0 + + sx = static_cast (cutneighmax*bininvx); + if (sx*binsizex < cutneighmax) sx++; + sy = static_cast (cutneighmax*bininvy); + if (sy*binsizey < cutneighmax) sy++; + sz = static_cast (cutneighmax*bininvz); + if (sz*binsizez < cutneighmax) sz++; + + return ((2*sx+1) * (2*sy+1) * (2*sz+1)); +} + +void NStencilFullBytype3d::create_setup() +{ + + int itype, jtype; + int maxtypes; + int smax; + + //printf("NStencilFullBytype3d::create_steup()\n"); + maxtypes = atom->ntypes; + + if (maxstencil_type == NULL) { + memory->create(maxstencil_type, maxtypes+1, maxtypes+1, "BAD A"); + memory->create(nstencil_type, maxtypes+1, maxtypes+1, "BAD B"); + stencil_type = new int**[maxtypes+1](); + for (itype = 1; itype <= maxtypes; ++itype) { + stencil_type[itype] = new int*[maxtypes+1](); + for (jtype = 1; jtype <= maxtypes; ++jtype) { + maxstencil_type[itype][jtype] = 0; + } + } + } + + // like -> like => use standard newtoff stencil at bin + + for (itype = 1; itype <= maxtypes; ++itype) { + copy_bin_info_bytype(itype); + smax = copy_neigh_info_bytype(itype); + if (smax > maxstencil_type[itype][itype]) { + maxstencil_type[itype][itype] = smax; + memory->destroy(stencil_type[itype][itype]); + memory->create(stencil_type[itype][itype], smax, + "NStencilFullBytype::create_steup() stencil"); + } + create_newtoff(itype, itype, cutneighmaxsq); + } + + // smaller -> larger => use existing newtoff stencil in larger bin + // larger -> smaller => use multi-like stencil for small-large in smaller bin + // If types are same cutoff, use existing like-like stencil. + + for (itype = 1; itype <= maxtypes; ++itype) { + for (jtype = 1; jtype <= maxtypes; ++jtype) { + if (itype == jtype) continue; + if (cuttypesq[itype] <= cuttypesq[jtype]) { + // Potential destroy/create problem? + nstencil_type[itype][jtype] = nstencil_type[jtype][jtype]; + stencil_type[itype][jtype] = stencil_type[jtype][jtype]; + } + else { + copy_bin_info_bytype(jtype); + // smax = copy_neigh_info_bytype(jtype); + + cutneighmaxsq = cuttypesq[jtype]; + cutneighmax = sqrt(cutneighmaxsq); + sx = static_cast (cutneighmax*bininvx); + if (sx*binsizex < cutneighmax) sx++; + sy = static_cast (cutneighmax*bininvy); + if (sy*binsizey < cutneighmax) sy++; + sz = static_cast (cutneighmax*bininvz); + if (sz*binsizez < cutneighmax) sz++; + + smax = (2*sx+1) * (2*sy+1) * (2*sz+1); + if (smax > maxstencil_type[itype][jtype]) { + maxstencil_type[itype][jtype] = smax; + memory->destroy(stencil_type[itype][jtype]); + memory->create(stencil_type[itype][jtype], smax, "Bad C"); + } + create_newtoff(itype, jtype, cuttypesq[jtype]); + } + } + } + + //for (itype = 1; itype <= maxtypes; itype++) { + // for (jtype = 1; jtype <= maxtypes; jtype++) { + // printf("i j n %d %d %d\n", itype, jtype, nstencil_type[itype][jtype]); + // printf("i j n %d %d %d\n", itype, jtype, maxstencil_type[itype][jtype]); + // } + // } + +} + +void NStencilFullBytype3d::create_newtoff(int itype, int jtype, double cutsq) { + + int i, j, k, ns; + + ns = 0; + + for (k = -sz; k <= sz; k++) + for (j = -sy; j <= sy; j++) + for (i = -sx; i <= sx; i++) + if (bin_distance(i,j,k) < cutsq) + stencil_type[itype][jtype][ns++] = k*mbiny*mbinx + j*mbinx + i; + + nstencil_type[itype][jtype] = ns; +} + +/* ---------------------------------------------------------------------- + create stencil based on bin geometry and cutoff +------------------------------------------------------------------------- */ + +void NStencilFullBytype3d::create() +{ + //int i,j,k; + + //nstencil = 0; + + //for (k = -sz; k <= sz; k++) + // for (j = -sy; j <= sy; j++) + // for (i = -sx; i <= sx; i++) + // if (bin_distance(i,j,k) < cutneighmaxsq) + // stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i; +} diff --git a/src/nstencil_full_bytype_3d.h b/src/nstencil_full_bytype_3d.h new file mode 100644 index 0000000000..6e61365d17 --- /dev/null +++ b/src/nstencil_full_bytype_3d.h @@ -0,0 +1,53 @@ +/* -*- 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 NSTENCIL_CLASS + +NStencilStyle(full/bytype/3d, + NStencilFullBytype3d, + NS_FULL | NS_BYTYPE | NS_3D | + NS_NEWTON | NS_NEWTOFF | NS_ORTHO | NS_TRI) + +#else + +#ifndef LMP_NSTENCIL_FULL_BYTYPE_3D_H +#define LMP_NSTENCIL_FULL_BYTYPE_3D_H + +#include "nstencil.h" + +namespace LAMMPS_NS { + +class NStencilFullBytype3d : public NStencil { + public: + NStencilFullBytype3d(class LAMMPS *); + ~NStencilFullBytype3d(); + void create(); + void create_setup(); + +private: + int ** maxstencil_type; + + void copy_bin_info_bytype(int); + int copy_neigh_info_bytype(int); + void create_newtoff(int, int, double); + +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/nstencil_half_bytype_2d_newton.cpp b/src/nstencil_half_bytype_2d_newton.cpp new file mode 100644 index 0000000000..18ed0972ff --- /dev/null +++ b/src/nstencil_half_bytype_2d_newton.cpp @@ -0,0 +1,192 @@ +/* ---------------------------------------------------------------------- + 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 "nstencil_half_bytype_2d_newton.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "nbin.h" +#include "memory.h" +#include "atom.h" + +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NStencilHalfBytype2dNewton::NStencilHalfBytype2dNewton(LAMMPS *lmp) : + NStencil(lmp) +{ + maxstencil_type = NULL; +} + +NStencilHalfBytype2dNewton::~NStencilHalfBytype2dNewton() { + + memory->destroy(nstencil_type); + for (int i = 1; i <= atom->ntypes; i++) { + for (int j = 0; j <= atom->ntypes; j++) { + if (maxstencil_type[i][j] > 0) memory->destroy(stencil_type[i][j]); + } + delete [] stencil_type[i]; + } + delete [] stencil_type; + memory->destroy(maxstencil_type); +} + +// KS To superclass +void NStencilHalfBytype2dNewton::copy_bin_info_bytype(int itype) { + + mbinx = nb->mbinx_type[itype]; + mbiny = nb->mbiny_type[itype]; + mbinz = nb->mbinz_type[itype]; + binsizex = nb->binsizex_type[itype]; + binsizey = nb->binsizey_type[itype]; + binsizez = nb->binsizez_type[itype]; + bininvx = nb->bininvx_type[itype]; + bininvy = nb->bininvy_type[itype]; + bininvz = nb->bininvz_type[itype]; +} + +// KS To superclass? +int NStencilHalfBytype2dNewton::copy_neigh_info_bytype(int itype) { + + cutneighmaxsq = neighbor->cutneighsq[itype][itype]; + cutneighmax = sqrt(cutneighmaxsq); + cuttypesq = neighbor->cuttypesq; + + // sx,sy,sz = max range of stencil in each dim + // smax = max possible size of entire 2d stencil + // stencil will be empty if cutneighmax = 0.0 + + sx = static_cast (cutneighmax*bininvx); + if (sx*binsizex < cutneighmax) sx++; + sy = static_cast (cutneighmax*bininvy); + if (sy*binsizey < cutneighmax) sy++; + sz = static_cast (cutneighmax*bininvz); + if (sz*binsizez < cutneighmax) sz++; + + return ((2*sx+1) * (2*sy+1) * (2*sz+1)); +} + +void NStencilHalfBytype2dNewton::create_setup() +{ + + int itype, jtype; + int maxtypes; + int smax; + + // maxstencil_type to superclass? + maxtypes = atom->ntypes; + + if (maxstencil_type == NULL) { + memory->create(maxstencil_type, maxtypes+1, maxtypes+1, "maxstencil_type"); + memory->create(nstencil_type, maxtypes+1, maxtypes+1, "nstencil_type"); + stencil_type = new int**[maxtypes+1](); + for (itype = 1; itype <= maxtypes; ++itype) { + stencil_type[itype] = new int*[maxtypes+1](); + for (jtype = 1; jtype <= maxtypes; ++jtype) { + maxstencil_type[itype][jtype] = 0; + nstencil_type[itype][jtype] = 0; + } + } + } + + // like -> like => use standard Newton stencil at bin + + for (itype = 1; itype <= maxtypes; ++itype) { + copy_bin_info_bytype(itype); + smax = copy_neigh_info_bytype(itype); + if (smax > maxstencil_type[itype][itype]) { + maxstencil_type[itype][itype] = smax; + memory->destroy(stencil_type[itype][itype]); + memory->create(stencil_type[itype][itype], smax, + "NStencilHalfBytypeNewton::create_steup() stencil"); + } + create_newton(itype, itype, cutneighmaxsq); + } + + // Cross types: "Newton on" reached by using Newton off stencil and + // looking one way through hierarchy + // smaller -> larger => use Newton off stencil in larger bin + // larger -> smaller => no nstecil required + // If cut offs are same, use existing type-type stencil + + for (itype = 1; itype <= maxtypes; ++itype) { + for (jtype = 1; jtype <= maxtypes; ++jtype) { + if (itype == jtype) continue; + if (cuttypesq[itype] == cuttypesq[jtype]) { + nstencil_type[itype][jtype] = nstencil_type[jtype][jtype]; + stencil_type[itype][jtype] = stencil_type[jtype][jtype]; + } + else if (cuttypesq[itype] < cuttypesq[jtype]) { + copy_bin_info_bytype(jtype); + + cutneighmaxsq = cuttypesq[jtype]; + cutneighmax = sqrt(cutneighmaxsq); + sx = static_cast (cutneighmax*bininvx); + if (sx*binsizex < cutneighmax) sx++; + sy = static_cast (cutneighmax*bininvy); + if (sy*binsizey < cutneighmax) sy++; + sz = static_cast (cutneighmax*bininvz); + if (sz*binsizez < cutneighmax) sz++; + + smax = (2*sx+1) * (2*sy+1) * (2*sz+1); + if (smax > maxstencil_type[itype][jtype]) { + maxstencil_type[itype][jtype] = smax; + memory->destroy(stencil_type[itype][jtype]); + memory->create(stencil_type[itype][jtype], smax, "stencil_type[]"); + } + create_newtoff(itype, jtype, cuttypesq[jtype]); + } + } + } + +} + +void NStencilHalfBytype2dNewton::create_newton(int itype, int jtype, double cutsq) { + + int i, j, ns; + + ns = 0; + + for (j = 0; j <= sy; j++) + for (i = -sx; i <= sx; i++) + if (j > 0 || (j == 0 && i > 0)) { + if (bin_distance(i,j,0) < cutsq) + stencil_type[itype][jtype][ns++] = j*mbinx + i; + } + nstencil_type[itype][jtype] = ns; +} + +void NStencilHalfBytype2dNewton::create_newtoff(int itype, int jtype, double cutsq) { + + int i, j, ns; + + ns = 0; + + for (j = -sy; j <= sy; j++) + for (i = -sx; i <= sx; i++) + if (bin_distance(i,j,0) < cutsq) + stencil_type[itype][jtype][ns++] = j*mbinx + i; + + nstencil_type[itype][jtype] = ns; +} + +/* ---------------------------------------------------------------------- + create stencil based on bin geometry and cutoff +------------------------------------------------------------------------- */ + +void NStencilHalfBytype2dNewton::create() +{ + // KS. Move "creation" here. +} diff --git a/src/nstencil_half_bytype_2d_newton.h b/src/nstencil_half_bytype_2d_newton.h new file mode 100644 index 0000000000..4d33c26a71 --- /dev/null +++ b/src/nstencil_half_bytype_2d_newton.h @@ -0,0 +1,52 @@ +/* -*- 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 NSTENCIL_CLASS + +NStencilStyle(half/bytype/2d/newton, + NStencilHalfBytype2dNewton, + NS_HALF | NS_BYTYPE | NS_2D | NS_NEWTON | NS_ORTHO) + +#else + +#ifndef LMP_NSTENCIL_HALF_BYTYPE_2D_NEWTON_H +#define LMP_NSTENCIL_HALF_BYTYPE_2D_NEWTON_H + +#include "nstencil.h" + +namespace LAMMPS_NS { + +class NStencilHalfBytype2dNewton : public NStencil { + public: + NStencilHalfBytype2dNewton(class LAMMPS *); + ~NStencilHalfBytype2dNewton(); + void create_setup(); + void create(); + + private: + int ** maxstencil_type; + + void copy_bin_info_bytype(int); + int copy_neigh_info_bytype(int); + void create_newton(int, int, double); + void create_newtoff(int, int, double); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/nstencil_half_bytype_3d_newtoff.cpp b/src/nstencil_half_bytype_3d_newtoff.cpp new file mode 100644 index 0000000000..f111715e5d --- /dev/null +++ b/src/nstencil_half_bytype_3d_newtoff.cpp @@ -0,0 +1,190 @@ +/* ---------------------------------------------------------------------- + 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 "nstencil_half_bytype_3d_newtoff.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "nbin.h" +#include "memory.h" +#include "atom.h" + +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NStencilHalfBytype3dNewtoff::NStencilHalfBytype3dNewtoff(LAMMPS *lmp) : + NStencil(lmp) +{ + maxstencil_type = NULL; +} + +NStencilHalfBytype3dNewtoff::~NStencilHalfBytype3dNewtoff() { + + memory->destroy(nstencil_type); + for (int i = 1; i <= atom->ntypes; i++) { + for (int j = 0; j <= atom->ntypes; j++) { + if (maxstencil_type[i][j] > 0) memory->destroy(stencil_type[i][j]); + } + delete [] stencil_type[i]; + } + delete [] stencil_type; + memory->destroy(maxstencil_type); +} + +void NStencilHalfBytype3dNewtoff::copy_bin_info_bytype(int itype) { + + mbinx = nb->mbinx_type[itype]; + mbiny = nb->mbiny_type[itype]; + mbinz = nb->mbinz_type[itype]; + binsizex = nb->binsizex_type[itype]; + binsizey = nb->binsizey_type[itype]; + binsizez = nb->binsizez_type[itype]; + bininvx = nb->bininvx_type[itype]; + bininvy = nb->bininvy_type[itype]; + bininvz = nb->bininvz_type[itype]; +} + +int NStencilHalfBytype3dNewtoff::copy_neigh_info_bytype(int itype) { + + cutneighmaxsq = neighbor->cutneighsq[itype][itype]; + cutneighmax = sqrt(cutneighmaxsq); + cuttypesq = neighbor->cuttypesq; + + // sx,sy,sz = max range of stencil in each dim + // smax = max possible size of entire 3d stencil + // stencil will be empty if cutneighmax = 0.0 + + sx = static_cast (cutneighmax*bininvx); + if (sx*binsizex < cutneighmax) sx++; + sy = static_cast (cutneighmax*bininvy); + if (sy*binsizey < cutneighmax) sy++; + sz = static_cast (cutneighmax*bininvz); + if (sz*binsizez < cutneighmax) sz++; + + return ((2*sx+1) * (2*sy+1) * (2*sz+1)); +} + +void NStencilHalfBytype3dNewtoff::create_setup() +{ + + int itype, jtype; + int maxtypes; + int smax; + + //printf("NStencilHalfBytype3dNewtoff::create_steup()\n"); + maxtypes = atom->ntypes; + + if (maxstencil_type == NULL) { + memory->create(maxstencil_type, maxtypes+1, maxtypes+1, "BAD A"); + memory->create(nstencil_type, maxtypes+1, maxtypes+1, "BAD B"); + stencil_type = new int**[maxtypes+1](); + for (itype = 1; itype <= maxtypes; ++itype) { + stencil_type[itype] = new int*[maxtypes+1](); + for (jtype = 1; jtype <= maxtypes; ++jtype) { + maxstencil_type[itype][jtype] = 0; + } + } + } + + // like -> like => use standard newtoff stencil at bin + + for (itype = 1; itype <= maxtypes; ++itype) { + copy_bin_info_bytype(itype); + smax = copy_neigh_info_bytype(itype); + if (smax > maxstencil_type[itype][itype]) { + maxstencil_type[itype][itype] = smax; + memory->destroy(stencil_type[itype][itype]); + memory->create(stencil_type[itype][itype], smax, + "NStencilHalfBytypeNewtoff::create_steup() stencil"); + } + create_newtoff(itype, itype, cutneighmaxsq); + } + + // smaller -> larger => use existing newtoff stencil in larger bin + // larger -> smaller => use multi-like stencil for small-large in smaller bin + // If types are same cutoff, use existing like-like stencil. + + for (itype = 1; itype <= maxtypes; ++itype) { + for (jtype = 1; jtype <= maxtypes; ++jtype) { + if (itype == jtype) continue; + if (cuttypesq[itype] <= cuttypesq[jtype]) { + // Potential destroy/create problem? + nstencil_type[itype][jtype] = nstencil_type[jtype][jtype]; + stencil_type[itype][jtype] = stencil_type[jtype][jtype]; + } + else { + copy_bin_info_bytype(jtype); + // smax = copy_neigh_info_bytype(jtype); + + cutneighmaxsq = cuttypesq[jtype]; + cutneighmax = sqrt(cutneighmaxsq); + sx = static_cast (cutneighmax*bininvx); + if (sx*binsizex < cutneighmax) sx++; + sy = static_cast (cutneighmax*bininvy); + if (sy*binsizey < cutneighmax) sy++; + sz = static_cast (cutneighmax*bininvz); + if (sz*binsizez < cutneighmax) sz++; + + smax = (2*sx+1) * (2*sy+1) * (2*sz+1); + if (smax > maxstencil_type[itype][jtype]) { + maxstencil_type[itype][jtype] = smax; + memory->destroy(stencil_type[itype][jtype]); + memory->create(stencil_type[itype][jtype], smax, "Bad C"); + } + create_newtoff(itype, jtype, cuttypesq[jtype]); + } + } + } + + //for (itype = 1; itype <= maxtypes; itype++) { + // for (jtype = 1; jtype <= maxtypes; jtype++) { + // printf("i j n %d %d %d\n", itype, jtype, nstencil_type[itype][jtype]); + // printf("i j n %d %d %d\n", itype, jtype, maxstencil_type[itype][jtype]); + // } + // } + +} + +void NStencilHalfBytype3dNewtoff::create_newtoff(int itype, int jtype, double cutsq) { + + int i, j, k, ns; + + ns = 0; + + for (k = -sz; k <= sz; k++) + for (j = -sy; j <= sy; j++) + for (i = -sx; i <= sx; i++) + if (bin_distance(i,j,k) < cutsq) + stencil_type[itype][jtype][ns++] = k*mbiny*mbinx + j*mbinx + i; + + nstencil_type[itype][jtype] = ns; +} + +/* ---------------------------------------------------------------------- + create stencil based on bin geometry and cutoff +------------------------------------------------------------------------- */ + +void NStencilHalfBytype3dNewtoff::create() +{ + //int i,j,k; + + //nstencil = 0; + + //for (k = -sz; k <= sz; k++) + // for (j = -sy; j <= sy; j++) + // for (i = -sx; i <= sx; i++) + // if (bin_distance(i,j,k) < cutneighmaxsq) + // stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i; +} diff --git a/src/nstencil_half_bytype_3d_newtoff.h b/src/nstencil_half_bytype_3d_newtoff.h new file mode 100644 index 0000000000..bf83a31e9a --- /dev/null +++ b/src/nstencil_half_bytype_3d_newtoff.h @@ -0,0 +1,51 @@ +/* -*- 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 NSTENCIL_CLASS + +NStencilStyle(half/bytype/3d/newtoff, + NStencilHalfBytype3dNewtoff, + NS_HALF | NS_BYTYPE | NS_3D | NS_NEWTOFF | NS_ORTHO | NS_TRI) + +#else + +#ifndef LMP_NSTENCIL_HALF_BYTYPE_3D_NEWTOFF_H +#define LMP_NSTENCIL_HALF_BYTYPE_3D_NEWTOFF_H + +#include "nstencil.h" + +namespace LAMMPS_NS { + +class NStencilHalfBytype3dNewtoff : public NStencil { + public: + NStencilHalfBytype3dNewtoff(class LAMMPS *); + ~NStencilHalfBytype3dNewtoff(); + void create_setup(); + void create(); + + private: + int ** maxstencil_type; + + void copy_bin_info_bytype(int); + int copy_neigh_info_bytype(int); + void create_newtoff(int, int, double); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/ diff --git a/src/nstencil_half_bytype_3d_newton.cpp b/src/nstencil_half_bytype_3d_newton.cpp new file mode 100644 index 0000000000..94c1ecd94d --- /dev/null +++ b/src/nstencil_half_bytype_3d_newton.cpp @@ -0,0 +1,194 @@ +/* ---------------------------------------------------------------------- + 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 "nstencil_half_bytype_3d_newton.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "nbin.h" +#include "memory.h" +#include "atom.h" + +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NStencilHalfBytype3dNewton::NStencilHalfBytype3dNewton(LAMMPS *lmp) : + NStencil(lmp) +{ + maxstencil_type = NULL; +} + +NStencilHalfBytype3dNewton::~NStencilHalfBytype3dNewton() { + + memory->destroy(nstencil_type); + for (int i = 1; i <= atom->ntypes; i++) { + for (int j = 0; j <= atom->ntypes; j++) { + if (maxstencil_type[i][j] > 0) memory->destroy(stencil_type[i][j]); + } + delete [] stencil_type[i]; + } + delete [] stencil_type; + memory->destroy(maxstencil_type); +} + +// KS To superclass +void NStencilHalfBytype3dNewton::copy_bin_info_bytype(int itype) { + + mbinx = nb->mbinx_type[itype]; + mbiny = nb->mbiny_type[itype]; + mbinz = nb->mbinz_type[itype]; + binsizex = nb->binsizex_type[itype]; + binsizey = nb->binsizey_type[itype]; + binsizez = nb->binsizez_type[itype]; + bininvx = nb->bininvx_type[itype]; + bininvy = nb->bininvy_type[itype]; + bininvz = nb->bininvz_type[itype]; +} + +// KS To superclass? +int NStencilHalfBytype3dNewton::copy_neigh_info_bytype(int itype) { + + cutneighmaxsq = neighbor->cutneighsq[itype][itype]; + cutneighmax = sqrt(cutneighmaxsq); + cuttypesq = neighbor->cuttypesq; + + // sx,sy,sz = max range of stencil in each dim + // smax = max possible size of entire 3d stencil + // stencil will be empty if cutneighmax = 0.0 + + sx = static_cast (cutneighmax*bininvx); + if (sx*binsizex < cutneighmax) sx++; + sy = static_cast (cutneighmax*bininvy); + if (sy*binsizey < cutneighmax) sy++; + sz = static_cast (cutneighmax*bininvz); + if (sz*binsizez < cutneighmax) sz++; + + return ((2*sx+1) * (2*sy+1) * (2*sz+1)); +} + +void NStencilHalfBytype3dNewton::create_setup() +{ + + int itype, jtype; + int maxtypes; + int smax; + + // maxstencil_type to superclass? + maxtypes = atom->ntypes; + + if (maxstencil_type == NULL) { + memory->create(maxstencil_type, maxtypes+1, maxtypes+1, "maxstencil_type"); + memory->create(nstencil_type, maxtypes+1, maxtypes+1, "nstencil_type"); + stencil_type = new int**[maxtypes+1](); + for (itype = 1; itype <= maxtypes; ++itype) { + stencil_type[itype] = new int*[maxtypes+1](); + for (jtype = 1; jtype <= maxtypes; ++jtype) { + maxstencil_type[itype][jtype] = 0; + nstencil_type[itype][jtype] = 0; + } + } + } + + // like -> like => use standard Newton stencil at bin + + for (itype = 1; itype <= maxtypes; ++itype) { + copy_bin_info_bytype(itype); + smax = copy_neigh_info_bytype(itype); + if (smax > maxstencil_type[itype][itype]) { + maxstencil_type[itype][itype] = smax; + memory->destroy(stencil_type[itype][itype]); + memory->create(stencil_type[itype][itype], smax, + "NStencilHalfBytypeNewton::create_steup() stencil"); + } + create_newton(itype, itype, cutneighmaxsq); + } + + // Cross types: "Newton on" reached by using Newton off stencil and + // looking one way through hierarchy + // smaller -> larger => use Newton off stencil in larger bin + // larger -> smaller => no nstecil required + // If cut offs are same, use existing type-type stencil + + for (itype = 1; itype <= maxtypes; ++itype) { + for (jtype = 1; jtype <= maxtypes; ++jtype) { + if (itype == jtype) continue; + if (cuttypesq[itype] == cuttypesq[jtype]) { + nstencil_type[itype][jtype] = nstencil_type[jtype][jtype]; + stencil_type[itype][jtype] = stencil_type[jtype][jtype]; + } + else if (cuttypesq[itype] < cuttypesq[jtype]) { + copy_bin_info_bytype(jtype); + + cutneighmaxsq = cuttypesq[jtype]; + cutneighmax = sqrt(cutneighmaxsq); + sx = static_cast (cutneighmax*bininvx); + if (sx*binsizex < cutneighmax) sx++; + sy = static_cast (cutneighmax*bininvy); + if (sy*binsizey < cutneighmax) sy++; + sz = static_cast (cutneighmax*bininvz); + if (sz*binsizez < cutneighmax) sz++; + + smax = (2*sx+1) * (2*sy+1) * (2*sz+1); + if (smax > maxstencil_type[itype][jtype]) { + maxstencil_type[itype][jtype] = smax; + memory->destroy(stencil_type[itype][jtype]); + memory->create(stencil_type[itype][jtype], smax, "stencil_type[]"); + } + create_newtoff(itype, jtype, cuttypesq[jtype]); + } + } + } + +} + +void NStencilHalfBytype3dNewton::create_newton(int itype, int jtype, double cutsq) { + + int i, j, k, ns; + + ns = 0; + + for (k = 0; k <= sz; k++) + for (j = -sy; j <= sy; j++) + for (i = -sx; i <= sx; i++) + if (k > 0 || j > 0 || (j == 0 && i > 0)) { + if (bin_distance(i,j,k) < cutsq) + stencil_type[itype][jtype][ns++] = k*mbiny*mbinx + j*mbinx + i; + } + nstencil_type[itype][jtype] = ns; +} + +void NStencilHalfBytype3dNewton::create_newtoff(int itype, int jtype, double cutsq) { + + int i, j, k, ns; + + ns = 0; + + for (k = -sz; k <= sz; k++) + for (j = -sy; j <= sy; j++) + for (i = -sx; i <= sx; i++) + if (bin_distance(i,j,k) < cutsq) + stencil_type[itype][jtype][ns++] = k*mbiny*mbinx + j*mbinx + i; + + nstencil_type[itype][jtype] = ns; +} + +/* ---------------------------------------------------------------------- + create stencil based on bin geometry and cutoff +------------------------------------------------------------------------- */ + +void NStencilHalfBytype3dNewton::create() +{ + // KS. Move "creation" here. +} diff --git a/src/nstencil_half_bytype_3d_newton.h b/src/nstencil_half_bytype_3d_newton.h new file mode 100644 index 0000000000..1245e2cd08 --- /dev/null +++ b/src/nstencil_half_bytype_3d_newton.h @@ -0,0 +1,52 @@ +/* -*- 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 NSTENCIL_CLASS + +NStencilStyle(half/bytype/3d/newton, + NStencilHalfBytype3dNewton, + NS_HALF | NS_BYTYPE | NS_3D | NS_NEWTON | NS_ORTHO | NS_TRI) + +#else + +#ifndef LMP_NSTENCIL_HALF_BYTYPE_3D_NEWTON_H +#define LMP_NSTENCIL_HALF_BYTYPE_3D_NEWTON_H + +#include "nstencil.h" + +namespace LAMMPS_NS { + +class NStencilHalfBytype3dNewton : public NStencil { + public: + NStencilHalfBytype3dNewton(class LAMMPS *); + ~NStencilHalfBytype3dNewton(); + void create_setup(); + void create(); + + private: + int ** maxstencil_type; + + void copy_bin_info_bytype(int); + int copy_neigh_info_bytype(int); + void create_newton(int, int, double); + void create_newtoff(int, int, double); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +*/