ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fixing merge conflicts

Browse Source
This commit is contained in:
Joel Clemmer
2020-10-02 15:09:47 -06:00
parent 51d55aa036
commit b1b014aed3
27 changed files with 2413 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
src/comm.cpp
View File

@ -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;

1
src/comm.h
View File

@ -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};

6
src/comm_brick.cpp
View File

@ -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]);
}
}
}

9
src/nbin.cpp
View File

@ -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;
}
/* ---------------------------------------------------------------------- */

14
src/nbin.h
View File

@ -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:

466
src/nbin_bytype.cpp Normal file
View File

@ -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 <assert.h>
#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<int> (bbox[0]*binsizeinv);
nbiny_type[n] = static_cast<int> (bbox[1]*binsizeinv);
if (dimension == 3) nbinz_type[n] = static_cast<int> (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<int> ((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<int> ((coord-bboxlo[0])*bininvx_type[n]);
coord = bsubboxlo[1] - SMALL*bbox[1];
mbinylo_type[n] = static_cast<int> ((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<int> ((coord-bboxlo[1])*bininvy_type[n]);
if (dimension == 3) {
coord = bsubboxlo[2] - SMALL*bbox[2];
mbinzlo_type[n] = static_cast<int> ((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<int> ((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<int> ((x[0]-bboxhi[0])*bininvx_type[it]) + nbinx_type[it];
else if (x[0] >= bboxlo[0]) {
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx_type[it]);
ix = MIN(ix,nbinx_type[it]-1);
} else
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx_type[it]) - 1;
if (x[1] >= bboxhi[1])
iy = static_cast<int> ((x[1]-bboxhi[1])*bininvy_type[it]) + nbiny_type[it];
else if (x[1] >= bboxlo[1]) {
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy_type[it]);
iy = MIN(iy,nbiny_type[it]-1);
} else
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy_type[it]) - 1;
if (x[2] >= bboxhi[2])
iz = static_cast<int> ((x[2]-bboxhi[2])*bininvz_type[it]) + nbinz_type[it];
else if (x[2] >= bboxlo[2]) {
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz_type[it]);
iz = MIN(iz,nbinz_type[it]-1);
} else
iz = static_cast<int> ((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;
}

56
src/nbin_bytype.h Normal file
View File

@ -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:
*/

2
src/nbin_standard.h
View File

@ -15,7 +15,7 @@
NBinStyle(standard,
NBinStandard,
0)
NB_STANDARD)
#else

17
src/neighbor.cpp
View File

@ -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);

6
src/neighbor.h
View File

@ -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;
}
}

144
src/npair_full_bytype.cpp Normal file
View File

@ -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<int> *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;
}

43
src/npair_full_bytype.h Normal file
View File

@ -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:
*/

219
src/npair_half_bytype_newton.cpp Executable file
View File

@ -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 <string.h>
#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<int> *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;
}

43
src/npair_half_bytype_newton.h Executable file
View File

@ -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:
*/

130
src/npair_half_size_bytype_newtoff.cpp Normal file
View File

@ -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 <string.h>
#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<int> *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;
}

43
src/npair_half_size_bytype_newtoff.h Normal file
View File

@ -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:
*/

189
src/npair_half_size_bytype_newton.cpp Executable file
View File

@ -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 <string.h>
#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<int> *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;
}

43
src/npair_half_size_bytype_newton.h Executable file
View File

@ -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:
*/

4
src/nstencil.h
View File

@ -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 *);

192
src/nstencil_full_bytype_3d.cpp Normal file
View File

@ -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 <math.h>
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<int> (cutneighmax*bininvx);
if (sx*binsizex < cutneighmax) sx++;
sy = static_cast<int> (cutneighmax*bininvy);
if (sy*binsizey < cutneighmax) sy++;
sz = static_cast<int> (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<int> (cutneighmax*bininvx);
if (sx*binsizex < cutneighmax) sx++;
sy = static_cast<int> (cutneighmax*bininvy);
if (sy*binsizey < cutneighmax) sy++;
sz = static_cast<int> (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;
}

53
src/nstencil_full_bytype_3d.h Normal file
View File

@ -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:
*/

192
src/nstencil_half_bytype_2d_newton.cpp Normal file
View File

@ -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 <math.h>
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<int> (cutneighmax*bininvx);
if (sx*binsizex < cutneighmax) sx++;
sy = static_cast<int> (cutneighmax*bininvy);
if (sy*binsizey < cutneighmax) sy++;
sz = static_cast<int> (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<int> (cutneighmax*bininvx);
if (sx*binsizex < cutneighmax) sx++;
sy = static_cast<int> (cutneighmax*bininvy);
if (sy*binsizey < cutneighmax) sy++;
sz = static_cast<int> (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.
}

52
src/nstencil_half_bytype_2d_newton.h Normal file
View File

@ -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:
*/

190
src/nstencil_half_bytype_3d_newtoff.cpp Normal file
View File

@ -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 <math.h>
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<int> (cutneighmax*bininvx);
if (sx*binsizex < cutneighmax) sx++;
sy = static_cast<int> (cutneighmax*bininvy);
if (sy*binsizey < cutneighmax) sy++;
sz = static_cast<int> (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<int> (cutneighmax*bininvx);
if (sx*binsizex < cutneighmax) sx++;
sy = static_cast<int> (cutneighmax*bininvy);
if (sy*binsizey < cutneighmax) sy++;
sz = static_cast<int> (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;
}

51
src/nstencil_half_bytype_3d_newtoff.h Normal file
View File

@ -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:
*/

194
src/nstencil_half_bytype_3d_newton.cpp Normal file
View File

@ -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 <math.h>
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<int> (cutneighmax*bininvx);
if (sx*binsizex < cutneighmax) sx++;
sy = static_cast<int> (cutneighmax*bininvy);
if (sy*binsizey < cutneighmax) sy++;
sz = static_cast<int> (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<int> (cutneighmax*bininvx);
if (sx*binsizex < cutneighmax) sx++;
sy = static_cast<int> (cutneighmax*bininvy);
if (sy*binsizey < cutneighmax) sy++;
sz = static_cast<int> (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.
}

52
src/nstencil_half_bytype_3d_newton.h Normal file
View File

@ -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:
*/