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

Misc comments, typos, and cleanups

Browse Source
This commit is contained in:
Joel Clemmer
2021-01-12 23:04:14 -07:00
parent d5f34f6296
commit a25c77e512
18 changed files with 285 additions and 476 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
doc/src/comm_modify.rst
View File

@ -11,7 +11,7 @@ Syntax
comm_modify keyword value ...
* zero or more keyword/value pairs may be appended
* keyword = *mode* or *cutoff* or *cutoff/multi* or *cutoff/bytype* or *group* or *vel*
* keyword = *mode* or *cutoff* or *cutoff/multi* or *multi/reduce* or *group* or *vel*
.. parsed-literal::
@ -20,7 +20,7 @@ Syntax
*cutoff/multi* type value
type = atom type or type range (supports asterisk notation)
value = Rcut (distance units) = communicate atoms for selected types from this far away
*cutoff/byptype* arg = none = communicate atoms for each type by a distance equal to the largest interaction distance for that type
*multi/reduce* arg = none = reduce number of communicated ghost atoms for multi style
*group* value = group-ID = only communicate atoms in the group
*vel* value = *yes* or *no* = do or do not communicate velocity info with ghost atoms
@ -96,10 +96,10 @@ using the usual asterisk notation can be given. For granular pair styles,
the default cutoff is set to the sum of the current maximum atomic radii
for each type.
The *multi/reduce* option applies to *multi* and automatically sets communication
The *multi/reduce* option applies to *multi* and sets communication
cutoffs for different sized particles based on the largest interaction distance
between two particles in the same multi grouping. This reduces the number of
ghost that need to be communicated This method is only compatible with the
ghost atoms that need to be communicated. This method is only compatible with the
*multi* neighbor style and requires only half neighbor lists and Newton on.
See the :doc:`neighbor multi <neighbor>` command for more information.

12
doc/src/neigh_modify.rst
View File

@ -47,9 +47,9 @@ Syntax
N = max number of neighbors of one atom
*binsize* value = size
size = bin size for neighbor list construction (distance units)
*multi/custom* values = N types
*multi/custom* values = N arg1 ... argN
N = number of custom groups
types = N separate types or groups of types (see below)
arg = N separate types or ranges of types (see below)
Examples
""""""""
@ -208,12 +208,12 @@ overhead. You must first specify the number of custom groups N to be
defined followed by N ranges of types. The range can be specified as a
single numeric value, or a wildcard asterisk can be used to specify a range
of values. This takes the form "\*" or "\*n" or "n\*" or "m\*n". For
example, if N = the number of atom types, then an asterisk with no numeric
values means all types from 1 to N. A leading asterisk means all types
from 1 to n (inclusive). A trailing asterisk means all types from n to N
example, if M = the number of atom types, then an asterisk with no numeric
values means all types from 1 to M. A leading asterisk means all types
from 1 to n (inclusive). A trailing asterisk means all types from n to M
(inclusive). A middle asterisk means all types from m to n (inclusive).
Note that any atom types not included in a custom group will be automatically
placed within a new, separate group.
placed within a separate group.
Restrictions
""""""""""""

2
examples/multi/in.granular
View File

@ -1,4 +1,4 @@
# Big colloid particles and small LJ particles
# Bidisperse set of grains
units lj
atom_style sphere

53
src/USER-OMP/npair_half_multi_newton_omp.cpp
View File

@ -100,62 +100,22 @@ void NPairHalfMultiNewtonOmp::build(NeighList *list)
if(igroup == jgroup) jbin = ibin;
else jbin = coord2bin(x[i], jgroup);
// if same size: uses half stencil so check central bin
if(cutmultisq[igroup][igroup] == cutmultisq[jgroup][jgroup]){
// if same size: use half stencil
if(igroup == jgroup){
if(igroup == jgroup) js = bins[i];
else js = binhead_multi[jgroup][jbin];
// if same group, implement with:
// loop over rest of atoms in i's bin, ghosts are at end of linked list
// if same group,
// if j is owned atom, store it, since j is beyond i in linked list
// if j is ghost, only store if j coords are "above and to the right" of i
js = bins[i];
for (j = js; j >= 0; j = bins[j]) {
if (j >= nlocal) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
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 {
// if different groups, implement with:
// loop over all atoms in jgroup bin
// if different groups,
// if j is owned atom, store it if j > i
// if j is ghost, only store if j coords are "above and to the right" of i
js = binhead_multi[jgroup][jbin];
for (j = js; j >= 0; j = bins[j]) {
if(j < i) continue;
if(igroup != jgroup and j < i) continue;
if (j >= nlocal) {
if (x[j][2] < ztmp) continue;
@ -190,7 +150,6 @@ void NPairHalfMultiNewtonOmp::build(NeighList *list)
}
}
}
}
// for all groups, loop over all atoms in other bins in stencil, store every pair
// stencil is empty if i larger than j

46
src/USER-OMP/npair_half_size_multi_newton_omp.cpp
View File

@ -89,55 +89,22 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list)
if(igroup == jgroup) jbin = ibin;
else jbin = coord2bin(x[i], jgroup);
// if same size: uses half stencil so check central bin
if(cutmultisq[igroup][igroup] == cutmultisq[jgroup][jgroup]){
// if same size: use half stencil
if(igroup == jgroup){
if(igroup == jgroup) js = bins[i];
else js = binhead_multi[jgroup][jbin];
// if same group, implement with:
// loop over rest of atoms in i's bin, ghosts are at end of linked list
// if same group,
// if j is owned atom, store it, since j is beyond i in linked list
// if j is ghost, only store if j coords are "above and to the right" of i
js = bins[i];
for (j = js; j >= 0; j = bins[j]) {
if (j >= nlocal) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
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 {
// if different groups, implement with:
// loop over all atoms in jgroup bin
// if different groups,
// if j is owned atom, store it if j > i
// if j is ghost, only store if j coords are "above and to the right" of i
js = binhead_multi[jgroup][jbin];
for (j = js; j >= 0; j = bins[j]) {
if(j < i) continue;
if(igroup != jgroup and j < i) continue;
if (j >= nlocal) {
if (x[j][2] < ztmp) continue;
@ -165,7 +132,6 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list)
}
}
}
}
// for all groups, loop over all atoms in other bins in stencil, store every pair
// stencil is empty if i larger than j

1
src/comm.cpp
View File

@ -242,7 +242,6 @@ void Comm::init()
for (int i = 0; i < nfix; i++)
if (fix[i]->maxexchange_dynamic) maxexchange_fix_dynamic = 1;
// Can't used multi/reduce communication with Newton off or full neighbor lits
if(multi_reduce){
if (force->newton == 0)
error->all(FLERR,"Cannot use multi/reduce communication with Newton off");

68
src/comm_brick.cpp
View File

@ -167,16 +167,12 @@ void CommBrick::setup()
error->warning(FLERR,"Communication cutoff is 0.0. No ghost atoms "
"will be generated. Atoms may get lost.");
if (triclinic == 0) {
prd = domain->prd;
sublo = domain->sublo;
subhi = domain->subhi;
cutghost[0] = cutghost[1] = cutghost[2] = cut;
if (mode == Comm::MULTI) {
if (multi_reduce) {
// If using multi/reduce, communicate itype particles a distance equal
// to the max of itype-jtype group interaction given smaller jtype group
// to the max of itype-jtype group interaction
// only consider smaller jtype groups
int igroup, jgroup;
double **cutmultisq = neighbor->cutmultisq;
int *map_type_multi = neighbor->map_type_multi;
@ -203,18 +199,23 @@ void CommBrick::setup()
}
}
} else {
// If using a single binlist, use the max itype-jtype interaction distance for communication
// otherwise, communicate a distance equal to the maximum interaction distance for each type
double *cuttype = neighbor->cuttype;
for (i = 1; i <= ntypes; i++) {
cut = 0.0;
if (cutusermulti) cut = cutusermulti[i];
cutghostmulti[i][0] = MAX(cut,cuttype[i]);
cutghostmulti[i][1] = MAX(cut,cuttype[i]);
cutghostmulti[i][2] = MAX(cut,cuttype[i]);
double tmp = 0.0;
if (cutusermulti) tmp = cutusermulti[i];
cutghostmulti[i][0] = MAX(tmp,cuttype[i]);
cutghostmulti[i][1] = MAX(tmp,cuttype[i]);
cutghostmulti[i][2] = MAX(tmp,cuttype[i]);
}
}
}
if (triclinic == 0) {
prd = domain->prd;
sublo = domain->sublo;
subhi = domain->subhi;
cutghost[0] = cutghost[1] = cutghost[2] = cut;
} else {
prd = domain->prd_lamda;
sublo = domain->sublo_lamda;
@ -227,46 +228,11 @@ void CommBrick::setup()
cutghost[1] = cut * length1;
length2 = h_inv[2];
cutghost[2] = cut * length2;
if (mode == Comm::MULTI) {
if (multi_reduce) {
// If using multi/reduce, communicate itype particles a distance equal
// to the max of itype-jtype group interaction given smaller jtype group
int igroup, jgroup;
double **cutmultisq = neighbor->cutmultisq;
int *map_type_multi = neighbor->map_type_multi;
if (mode == Comm::MULTI){
for (i = 1; i <= ntypes; i++) {
if (cutusermulti) {
cutghostmulti[i][0] = cutusermulti[i];
cutghostmulti[i][1] = cutusermulti[i];
cutghostmulti[i][2] = cutusermulti[i];
} else {
cutghostmulti[i][0] = 0.0;
cutghostmulti[i][1] = 0.0;
cutghostmulti[i][2] = 0.0;
}
igroup = map_type_multi[i];
for (j = 1; j <= ntypes; j++){
jgroup = map_type_multi[j];
if(cutmultisq[jgroup][jgroup] > cutmultisq[igroup][igroup]) continue;
cutghostmulti[i][0] = length0 * MAX(cutghostmulti[i][0],sqrt(cutmultisq[igroup][jgroup]));
cutghostmulti[i][1] = length1 * MAX(cutghostmulti[i][1],sqrt(cutmultisq[igroup][jgroup]));
cutghostmulti[i][2] = length2 * MAX(cutghostmulti[i][2],sqrt(cutmultisq[igroup][jgroup]));
}
}
} else {
// If using a single binlist, use the max itype-jtype interaction distance for communication
double *cuttype = neighbor->cuttype;
for (i = 1; i <= ntypes; i++) {
cut = 0.0;
if (cutusermulti) cut = cutusermulti[i];
cutghostmulti[i][0] = length0 * MAX(cut,cuttype[i]);
cutghostmulti[i][1] = length1 * MAX(cut,cuttype[i]);
cutghostmulti[i][2] = length2 * MAX(cut,cuttype[i]);
}
cutghostmulti[i][0] *= length0;
cutghostmulti[i][1] *= length1;
cutghostmulti[i][2] *= length2;
}
}
}

16
src/comm_tiled.cpp
View File

@ -175,10 +175,10 @@ void CommTiled::setup()
// check that cutoff < any periodic box length
if (mode == Comm::MULTI) {
double cut;
if (multi_reduce) {
// If using multi/reduce, communicate itype particles a distance equal
// to the max of itype-jtype group interaction given smaller jtype group
// to the max of itype-jtype group interaction
// only consider smaller jtype groups
int igroup, jgroup;
double **cutmultisq = neighbor->cutmultisq;
int *map_type_multi = neighbor->map_type_multi;
@ -205,14 +205,14 @@ void CommTiled::setup()
}
}
} else {
// If using a single binlist, use the max itype-jtype interaction distance for communication
// otherwise, communicate a distance equal to the maximum interaction distance for each type
double *cuttype = neighbor->cuttype;
for (i = 1; i <= ntypes; i++) {
cut = 0.0;
if (cutusermulti) cut = cutusermulti[i];
cutghostmulti[i][0] = MAX(cut,cuttype[i]);
cutghostmulti[i][1] = MAX(cut,cuttype[i]);
cutghostmulti[i][2] = MAX(cut,cuttype[i]);
double tmp = 0.0;
if (cutusermulti) tmp = cutusermulti[i];
cutghostmulti[i][0] = MAX(tmp,cuttype[i]);
cutghostmulti[i][1] = MAX(tmp,cuttype[i]);
cutghostmulti[i][2] = MAX(tmp,cuttype[i]);
}
}
}

67
src/nbin_multi.cpp
View File

@ -59,41 +59,41 @@ void NBinMulti::bin_atoms_setup(int nall)
}
}
/* ---------------------------------------------------------------------
Identify index of group with smallest cutoff
------------------------------------------------------------------------ */
int NBinMulti::igroup_min()
{
int imin = 0;
for (int n = 0; n < maxgroups; n++)
if (cutmultisq[n][n] < cutmultisq[imin][imin]) imin = n;
return imin;
}
/* ----------------------------------------------------------------------
setup neighbor binning geometry
---------------------------------------------------------------------- */
bin numbering in each dimension is global:
0 = 0.0 to binsize, 1 = binsize to 2*binsize, etc
nbin-1,nbin,etc = bbox-binsize to bbox, bbox to bbox+binsize, etc
-1,-2,etc = -binsize to 0.0, -2*binsize to -binsize, etc
code will work for any binsize
since next(xyz) and stencil extend as far as necessary
binsize = 1/2 of cutoff is roughly optimal
for orthogonal boxes:
a dim must be filled exactly by integer # of bins
in periodic, procs on both sides of PBC must see same bin boundary
in non-periodic, coord2bin() still assumes this by use of nbin xyz
for triclinic boxes:
tilted simulation box cannot contain integer # of bins
stencil & neigh list built differently to account for this
mbinlo_multi = lowest global bin any of my ghost atoms could fall into for each grouping
mbinhi_multi = highest global bin any of my ghost atoms could fall into for each grouping
mbin_multi = number of bins I need in a dimension for each grouping
------------------------------------------------------------------------- */
void NBinMulti::setup_bins(int style)
{
int n;
int igroupmin;
// Initialize arrays
if (n_multi_groups > maxgroups) {
// Clear any/all memory for existing types
// Clear any/all memory for existing groupings
for (n = 0; n < maxgroups; n++)
memory->destroy(binhead_multi[n]);
delete [] binhead_multi;
// Realloacte at updated maxtypes
// Realloacte at updated maxgroups
maxgroups = n_multi_groups;
binhead_multi = new int*[maxgroups]();
@ -138,14 +138,18 @@ void NBinMulti::setup_bins(int style)
memory->destroy(maxbins_multi);
memory->create(maxbins_multi, maxgroups, "neigh:maxbins_multi");
// make sure reallocation occurs in bin_atoms_setup()
// ensure reallocation occurs in bin_atoms_setup()
for (n = 0; n < maxgroups; n++) {
maxbins_multi[n] = 0;
}
maxatom = 0;
}
igroupmin = igroup_min();
// Identify smallest group
int igroupmin = 0;
for (n = 0; n < maxgroups; n++)
if (cutmultisq[n][n] < cutmultisq[igroupmin][igroupmin])
igroupmin = n;
// bbox = size of bbox of entire domain
// bsubbox lo/hi = bounding box of my subdomain extended by comm->cutghost
@ -182,16 +186,18 @@ void NBinMulti::setup_bins(int style)
// For each grouping...
double binsize_optimal, binsizeinv, coord;
int mbinxhi,mbinyhi,mbinzhi;
for (n = 0; n < maxgroups; n++) {
// binsize_user only relates to smallest type
// optimal bin size is roughly 1/2 the type-type cutoff
// binsize_user only relates to smallest group
// optimal bin size is roughly 1/2 the group-group cutoff
// special case of all cutoffs = 0.0, binsize = box size
double binsize_optimal;
if (n == igroupmin && binsizeflag) binsize_optimal = binsize_user;
else binsize_optimal = 0.5*sqrt(cutmultisq[n][n]);
if (binsize_optimal == 0.0) binsize_optimal = bbox[0];
double binsizeinv = 1.0/binsize_optimal;
binsizeinv = 1.0/binsize_optimal;
// test for too many global bins in any dimension due to huge global domain
@ -216,7 +222,7 @@ void NBinMulti::setup_bins(int style)
// 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
// in that extreme case, cannot use multi, should use NSQ not BIN neighbor style
binsizex_multi[n] = bbox[0]/nbinx_multi[n];
binsizey_multi[n] = bbox[1]/nbiny_multi[n];
@ -236,9 +242,6 @@ void NBinMulti::setup_bins(int style)
// 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_multi[n] = static_cast<int> ((coord-bboxlo[0])*bininvx_multi[n]);
if (coord < bboxlo[0]) mbinxlo_multi[n] = mbinxlo_multi[n] - 1;
@ -340,10 +343,8 @@ void NBinMulti::bin_atoms()
double NBinMulti::memory_usage()
{
double bytes = 0;
for (int m = 0; m < maxgroups; m++) {
for (int m = 0; m < maxgroups; m++)
bytes += maxbins_multi[m]*sizeof(int);
bytes += 2*maxatom*sizeof(int);
}
return bytes;
}

4
src/nbin_multi.h
View File

@ -35,10 +35,6 @@ class NBinMulti : public NBin {
void setup_bins(int);
void bin_atoms();
double memory_usage();
private:
int igroup_min();
};
}

21
src/neighbor.cpp
View File

@ -83,6 +83,11 @@ static const char cite_neigh_multi[] =
" volume = 179,\n"
" pages = {320--329}\n"
"}\n\n"
"@article{Stratford2018,\n"
" author = {Stratford, Kevin and Shire, Tom and Hanley, Kevin},\n"
" title = {Implementation of multi-level contact detection in LAMMPS},\n"
" year = {2018}\n"
"}\n\n"
"@article{Shire2020,\n"
" author = {Shire, Tom and Hanley, Kevin J. and Stratford, Kevin},\n"
" title = {DEM simulations of polydisperse media: efficient contact\n"
@ -348,7 +353,8 @@ void Neighbor::init()
// multi cutoffs
if(style == Neighbor::MULTI){
int igroup, jgroup;
// If not defined, create default mapping
// If not defined from custom grouping, create default map
if(not map_type_multi) {
n_multi_groups = n;
memory->create(map_type_multi,n+1,"neigh:map_type_multi");
@ -1678,12 +1684,11 @@ 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 MULTI_OLD and must match
if (style == Neighbor::BIN || style == Neighbor::MULTI_OLD) {
if (!(mask & NB_STANDARD)) continue;
} else if (style == Neighbor::MULTI) {
// multi neighbor style require multi bin style
if (style == Neighbor::MULTI) {
if (!(mask & NB_MULTI)) continue;
} else {
if (!(mask & NB_STANDARD)) continue;
}
return i+1;
@ -1719,7 +1724,6 @@ int Neighbor::choose_stencil(NeighRequest *rq)
else if (rq->newton == 2) newtflag = 0;
// request a full stencil if building full neighbor list or newton is off
int fullflag = 0;
if (rq->full) fullflag = 1;
if (!newtflag) fullflag = 1;
@ -1972,7 +1976,6 @@ NPair *Neighbor::pair_creator(LAMMPS *lmp)
/* ----------------------------------------------------------------------
setup neighbor binning and neighbor stencils
called before run and every reneighbor if box size/shape changes
initialize default settings for multi before run
only operates on perpetual lists
build_one() operates on occasional lists
------------------------------------------------------------------------- */
@ -2458,7 +2461,7 @@ void Neighbor::modify_params(int narg, char **arg)
n_multi_groups += 1;
}
// Create seprate group for each undefined atom type
// Create separate group for each undefined atom type
for(i = 1; i <= ntypes; i++){
if(map_type_multi[i] == -1){
map_type_multi[i] = n_multi_groups;

8
src/npair.cpp
View File

@ -304,9 +304,9 @@ int NPair::coord2bin(double *x, int ig)
} else
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz_multi[ig]) - 1;
ibin = (iz-mbinzlo_multi[ig])*mbiny_multi[ig]*mbinx_multi[ig]
+ (iy-mbinylo_multi[ig])*mbinx_multi[ig]
+ (ix-mbinxlo_multi[ig]);
ix -= mbinxlo_multi[ig];
iy -= mbinylo_multi[ig];
iz -= mbinzlo_multi[ig];
ibin = iz*mbiny_multi[ig]*mbinx_multi[ig] + iy*mbinx_multi[ig] + ix;
return ibin;
}

2
src/npair.h
View File

@ -115,7 +115,7 @@ class NPair : protected Pointers {
int coord2bin(double *); // mapping atom coord to a bin
int coord2bin(double *, int &, int &, int&); // ditto
int coord2bin(double *, int); // mapping atom coord to type bin
int coord2bin(double *, int); // mapping atom coord to group bin
// find_special: determine if atom j is in special list of atom i

53
src/npair_half_multi_newton.cpp
View File

@ -88,62 +88,22 @@ void NPairHalfMultiNewton::build(NeighList *list)
if(igroup == jgroup) jbin = ibin;
else jbin = coord2bin(x[i], jgroup);
// if same size: uses half stencil so check central bin
if(cutmultisq[igroup][igroup] == cutmultisq[jgroup][jgroup]){
// if same size: use half stencil
if(igroup == jgroup){
if(igroup == jgroup) js = bins[i];
else js = binhead_multi[jgroup][jbin];
// if same group, implement with:
// loop over rest of atoms in i's bin, ghosts are at end of linked list
// if same group,
// if j is owned atom, store it, since j is beyond i in linked list
// if j is ghost, only store if j coords are "above and to the right" of i
js = bins[i];
for (j = js; j >= 0; j = bins[j]) {
if (j >= nlocal) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
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 {
// if different groups, implement with:
// loop over all atoms in jgroup bin
// if different groups,
// if j is owned atom, store it if j > i
// if j is ghost, only store if j coords are "above and to the right" of i
js = binhead_multi[jgroup][jbin];
for (j = js; j >= 0; j = bins[j]) {
if(j < i) continue;
if(igroup != jgroup and j < i) continue;
if (j >= nlocal) {
if (x[j][2] < ztmp) continue;
@ -178,7 +138,6 @@ void NPairHalfMultiNewton::build(NeighList *list)
}
}
}
}
// for all groups, loop over all atoms in other bins in stencil, store every pair
// stencil is empty if i larger than j

46
src/npair_half_size_multi_newton.cpp
View File

@ -77,55 +77,22 @@ void NPairHalfSizeMultiNewton::build(NeighList *list)
if(igroup == jgroup) jbin = ibin;
else jbin = coord2bin(x[i], jgroup);
// if same size: uses half stencil so check central bin
if(cutmultisq[igroup][igroup] == cutmultisq[jgroup][jgroup]){
// if same size: use half stencil
if(igroup == jgroup){
if(igroup == jgroup) js = bins[i];
else js = binhead_multi[jgroup][jbin];
// if same group, implement with:
// loop over rest of atoms in i's bin, ghosts are at end of linked list
// if same group,
// if j is owned atom, store it, since j is beyond i in linked list
// if j is ghost, only store if j coords are "above and to the right" of i
js = bins[i];
for (j = js; j >= 0; j = bins[j]) {
if (j >= nlocal) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
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 {
// if different groups, implement with:
// loop over all atoms in jgroup bin
// if different groups,
// if j is owned atom, store it if j > i
// if j is ghost, only store if j coords are "above and to the right" of i
js = binhead_multi[jgroup][jbin];
for (j = js; j >= 0; j = bins[j]) {
if(j < i) continue;
if(igroup != jgroup and j < i) continue;
if (j >= nlocal) {
if (x[j][2] < ztmp) continue;
@ -153,7 +120,6 @@ void NPairHalfSizeMultiNewton::build(NeighList *list)
}
}
}
}
// for all groups, loop over all atoms in other bins in stencil, store every pair
// stencil is empty if i larger than j

1
src/npair_half_size_multi_newton_tri.cpp
View File

@ -77,7 +77,6 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
if(igroup == jgroup) jbin = ibin;
else jbin = coord2bin(x[i], jgroup);
// loop over all atoms in bins in stencil
// stencil is empty if i larger than j
// stencil is half if i same size as j

27
src/nstencil.cpp
View File

@ -55,7 +55,7 @@ using namespace LAMMPS_NS;
create one stencil for each igroup-jgroup pairing
the full/half stencil label refers to the same-group stencil
a half list with newton on has a half same-group stencil
a full list or half list with newton of has a full same-group stencil
a full list or half list with newton off has a full same-group stencil
cross group stencils are always full to allow small-to-large lookups
for orthogonal boxes, a half stencil includes bins to the "upper right" of central bin
for triclinic, a half stencil includes bins in the z (3D) or y (2D) plane of self and above
@ -230,7 +230,7 @@ void NStencil::create_setup()
int smax = (2*sx+1) * (2*sy+1) * (2*sz+1);
// reallocate stencil structs if necessary
// for BIN and MULTI styles
// for BIN and MULTI_OLD styles
if (neighstyle == Neighbor::BIN) {
if (smax > maxstencil) {
@ -315,9 +315,8 @@ void NStencil::create_setup()
set_stencil_properties();
// Allocate arrays to store stencils
if (!maxstencil_multi) {
memory->create(maxstencil_multi, n, n, "neighstencil::stencil_multi");
memory->create(maxstencil_multi, n, n, "neighstencil::maxstencil_multi");
memory->create(nstencil_multi, n, n, "neighstencil::nstencil_multi");
stencil_multi = new int**[n]();
for (i = 0; i < n; ++i) {
@ -325,6 +324,7 @@ void NStencil::create_setup()
for (j = 0; j < n; ++j) {
maxstencil_multi[i][j] = 0;
nstencil_multi[i][j] = 0;
stencil_multi[i][j] = nullptr;
}
}
}
@ -400,21 +400,21 @@ double NStencil::bin_distance(int i, int j, int k)
compute closest distance for a given atom grouping
------------------------------------------------------------------------- */
double NStencil::bin_distance_multi(int i, int j, int k, int group)
double NStencil::bin_distance_multi(int i, int j, int k, int ig)
{
double delx,dely,delz;
if (i > 0) delx = (i-1)*binsizex_multi[group];
if (i > 0) delx = (i-1)*binsizex_multi[ig];
else if (i == 0) delx = 0.0;
else delx = (i+1)*binsizex_multi[group];
else delx = (i+1)*binsizex_multi[ig];
if (j > 0) dely = (j-1)*binsizey_multi[group];
if (j > 0) dely = (j-1)*binsizey_multi[ig];
else if (j == 0) dely = 0.0;
else dely = (j+1)*binsizey_multi[group];
else dely = (j+1)*binsizey_multi[ig];
if (k > 0) delz = (k-1)*binsizez_multi[group];
if (k > 0) delz = (k-1)*binsizez_multi[ig];
else if (k == 0) delz = 0.0;
else delz = (k+1)*binsizez_multi[group];
else delz = (k+1)*binsizez_multi[ig];
return (delx*delx + dely*dely + delz*delz);
}
@ -435,13 +435,8 @@ double NStencil::memory_usage()
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
bytes += maxstencil_multi[i][j] * sizeof(int);
bytes += maxstencil_multi[i][j] * sizeof(int);
bytes += maxstencil_multi[i][j] * sizeof(double);
}
}
bytes += 2 * n * n * sizeof(bool);
bytes += 6 * n * n * sizeof(int);
bytes += 4 * n * n * sizeof(double);
}
return bytes;
}

6
src/nstencil.h
View File

@ -30,7 +30,7 @@ class NStencil : protected Pointers {
int *nstencil_multi_old; // # bins in each type-based old multi stencil
int **stencil_multi_old; // list of bin offsets in each stencil
double **distsq_multi_old; // sq distances to bins in each stencil
int ** nstencil_multi; // # bins bins in each itype-jtype multi stencil
int ** nstencil_multi; // # bins bins in each igroup-jgroup multi stencil
int *** stencil_multi; // list of bin offsets in each multi stencil
int ** maxstencil_multi; // max stencil size for each multi stencil
@ -42,8 +42,8 @@ class NStencil : protected Pointers {
double cutoff_custom; // cutoff set by requestor
// Arrays to store options for multi itype-jtype stencils
bool **flag_half_multi; // flag creation of a half stencil for itype-jtype
bool **flag_skip_multi; // skip creation of itype-jtype stencils (for newton on)
bool **flag_half_multi; // flag creation of a half stencil for igroup-jgroup
bool **flag_skip_multi; // skip creation of igroup-jgroup stencils (for newton on)
int **bin_group_multi; // what group to use for bin information
NStencil(class LAMMPS *);