diff --git a/doc/src/comm_modify.rst b/doc/src/comm_modify.rst index 5f084561a1..a3f7dba4eb 100644 --- a/doc/src/comm_modify.rst +++ b/doc/src/comm_modify.rst @@ -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 ` command for more information. diff --git a/doc/src/neigh_modify.rst b/doc/src/neigh_modify.rst index 33f24af96e..c6f573ef08 100644 --- a/doc/src/neigh_modify.rst +++ b/doc/src/neigh_modify.rst @@ -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 """""""""""" diff --git a/examples/multi/in.granular b/examples/multi/in.granular index 00dd4ad2bb..8f1743f90d 100644 --- a/examples/multi/in.granular +++ b/examples/multi/in.granular @@ -1,4 +1,4 @@ -# Big colloid particles and small LJ particles +# Bidisperse set of grains units lj atom_style sphere diff --git a/src/USER-OMP/npair_half_multi_newton_omp.cpp b/src/USER-OMP/npair_half_multi_newton_omp.cpp index a7a3e6d6f7..11a84a3040 100755 --- a/src/USER-OMP/npair_half_multi_newton_omp.cpp +++ b/src/USER-OMP/npair_half_multi_newton_omp.cpp @@ -100,96 +100,55 @@ 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 same group, implement with: - // loop over rest of atoms in i's bin, ghosts are at end of linked list - // if j is owned atom, store it, since j is beyond i in linked list - // if j is ghost, only store if j coords are "above and to the right" of i - - js = bins[i]; + if(igroup == jgroup) js = bins[i]; + else js = binhead_multi[jgroup][jbin]; - 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 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 + // 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 = binhead_multi[jgroup][jbin]; + // 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 - for (j = js; j >= 0; j = bins[j]) { - if(j < i) continue; + for (j = js; j >= 0; j = bins[j]) { + if(igroup != jgroup and j < i) continue; - if (j >= nlocal) { - if (x[j][2] < ztmp) continue; - if (x[j][2] == ztmp) { - if (x[j][1] < ytmp) continue; - if (x[j][1] == ytmp && x[j][0] < xtmp) continue; - } + if (j >= nlocal) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) continue; } - - jtype = type[j]; - if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - - if (rsq <= cutneighsq[itype][jtype]) { - if (molecular) { - if (!moltemplate) - which = find_special(special[i],nspecial[i],tag[j]); - else if (imol >= 0) - which = find_special(onemols[imol]->special[iatom], - onemols[imol]->nspecial[iatom], - tag[j]-tagprev); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (domain->minimum_image_check(delx,dely,delz)) - neighptr[n++] = j; - else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); - } else neighptr[n++] = j; - } - } - } + } + + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq <= cutneighsq[itype][jtype]) { + if (molecular) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } } // for all groups, loop over all atoms in other bins in stencil, store every pair diff --git a/src/USER-OMP/npair_half_size_multi_newton_omp.cpp b/src/USER-OMP/npair_half_size_multi_newton_omp.cpp index 0137873600..cec3edd02b 100755 --- a/src/USER-OMP/npair_half_size_multi_newton_omp.cpp +++ b/src/USER-OMP/npair_half_size_multi_newton_omp.cpp @@ -89,80 +89,46 @@ 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 same group, implement with: - // loop over rest of atoms in i's bin, ghosts are at end of linked list - // if j is owned atom, store it, since j is beyond i in linked list - // if j is ghost, only store if j coords are "above and to the right" of i - - 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(igroup == jgroup) js = bins[i]; + else js = binhead_multi[jgroup][jbin]; + + // 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 + + // 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 - if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + for (j = js; j >= 0; j = bins[j]) { + if(igroup != jgroup and j < i) continue; + + 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; } } - } else { - - // if different groups, implement with: - // loop over all atoms in jgroup bin - // 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 + + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; - js = binhead_multi[jgroup][jbin]; + 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); - for (j = js; j >= 0; j = bins[j]) { - if(j < i) continue; - - 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; - } + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; } } } diff --git a/src/comm.cpp b/src/comm.cpp index 542ddfdf10..867eab6256 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -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"); diff --git a/src/comm_brick.cpp b/src/comm_brick.cpp index 74865476f5..939669fb75 100644 --- a/src/comm_brick.cpp +++ b/src/comm_brick.cpp @@ -167,54 +167,55 @@ void CommBrick::setup() error->warning(FLERR,"Communication cutoff is 0.0. No ghost atoms " "will be generated. Atoms may get lost."); + + 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 + // only consider smaller jtype groups + int igroup, jgroup; + double **cutmultisq = neighbor->cutmultisq; + int *map_type_multi = neighbor->map_type_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] = MAX(cutghostmulti[i][0],sqrt(cutmultisq[igroup][jgroup])); + cutghostmulti[i][1] = MAX(cutghostmulti[i][1],sqrt(cutmultisq[igroup][jgroup])); + cutghostmulti[i][2] = MAX(cutghostmulti[i][2],sqrt(cutmultisq[igroup][jgroup])); + } + } + } else { + // otherwise, communicate a distance equal to the maximum interaction distance for each type + double *cuttype = neighbor->cuttype; + for (i = 1; i <= ntypes; 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; - - 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; - 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] = MAX(cutghostmulti[i][0],sqrt(cutmultisq[igroup][jgroup])); - cutghostmulti[i][1] = MAX(cutghostmulti[i][1],sqrt(cutmultisq[igroup][jgroup])); - cutghostmulti[i][2] = 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] = MAX(cut,cuttype[i]); - cutghostmulti[i][1] = MAX(cut,cuttype[i]); - cutghostmulti[i][2] = MAX(cut,cuttype[i]); - } - } - } - } 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; - 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]); - } + if (mode == Comm::MULTI){ + for (i = 1; i <= ntypes; i++) { + cutghostmulti[i][0] *= length0; + cutghostmulti[i][1] *= length1; + cutghostmulti[i][2] *= length2; } } } diff --git a/src/comm_tiled.cpp b/src/comm_tiled.cpp index 64d1207a64..06edbf7416 100644 --- a/src/comm_tiled.cpp +++ b/src/comm_tiled.cpp @@ -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]); } } } diff --git a/src/nbin_multi.cpp b/src/nbin_multi.cpp index b8832bdfea..04756ff803 100644 --- a/src/nbin_multi.cpp +++ b/src/nbin_multi.cpp @@ -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 ((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); - } + bytes += 2*maxatom*sizeof(int); return bytes; } diff --git a/src/nbin_multi.h b/src/nbin_multi.h index 84ec0348bc..490af10dbb 100644 --- a/src/nbin_multi.h +++ b/src/nbin_multi.h @@ -35,10 +35,6 @@ class NBinMulti : public NBin { void setup_bins(int); void bin_atoms(); double memory_usage(); - - private: - - int igroup_min(); }; } diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 94add1f324..5c8d8ec594 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -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; diff --git a/src/npair.cpp b/src/npair.cpp index facd9086d1..95a61fbefa 100644 --- a/src/npair.cpp +++ b/src/npair.cpp @@ -304,9 +304,9 @@ int NPair::coord2bin(double *x, int ig) } else iz = static_cast ((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; -} \ No newline at end of file +} diff --git a/src/npair.h b/src/npair.h index d1fcf2e6b2..c857ec0924 100644 --- a/src/npair.h +++ b/src/npair.h @@ -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 diff --git a/src/npair_half_multi_newton.cpp b/src/npair_half_multi_newton.cpp index 777aab5826..aa8f606242 100755 --- a/src/npair_half_multi_newton.cpp +++ b/src/npair_half_multi_newton.cpp @@ -88,96 +88,55 @@ 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 same group, implement with: - // loop over rest of atoms in i's bin, ghosts are at end of linked list - // if j is owned atom, store it, since j is beyond i in linked list - // if j is ghost, only store if j coords are "above and to the right" of i - - js = bins[i]; + if(igroup == jgroup) js = bins[i]; + else js = binhead_multi[jgroup][jbin]; - 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 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 + // 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 = binhead_multi[jgroup][jbin]; + // 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 - for (j = js; j >= 0; j = bins[j]) { - if(j < i) continue; + for (j = js; j >= 0; j = bins[j]) { + if(igroup != jgroup and j < i) continue; - if (j >= nlocal) { - if (x[j][2] < ztmp) continue; - if (x[j][2] == ztmp) { - if (x[j][1] < ytmp) continue; - if (x[j][1] == ytmp && x[j][0] < xtmp) continue; - } + if (j >= nlocal) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) continue; } - - jtype = type[j]; - if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + } - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - - if (rsq <= cutneighsq[itype][jtype]) { - if (molecular) { - if (!moltemplate) - which = find_special(special[i],nspecial[i],tag[j]); - else if (imol >= 0) - which = find_special(onemols[imol]->special[iatom], - onemols[imol]->nspecial[iatom], - tag[j]-tagprev); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (domain->minimum_image_check(delx,dely,delz)) - neighptr[n++] = j; - else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); - } else neighptr[n++] = j; - } - } - } + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq <= cutneighsq[itype][jtype]) { + if (molecular) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } } // for all groups, loop over all atoms in other bins in stencil, store every pair diff --git a/src/npair_half_size_multi_newton.cpp b/src/npair_half_size_multi_newton.cpp index 500e06bc35..84571a3b02 100755 --- a/src/npair_half_size_multi_newton.cpp +++ b/src/npair_half_size_multi_newton.cpp @@ -77,80 +77,46 @@ 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 same group, implement with: - // loop over rest of atoms in i's bin, ghosts are at end of linked list - // if j is owned atom, store it, since j is beyond i in linked list - // if j is ghost, only store if j coords are "above and to the right" of i - - 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(igroup == jgroup) js = bins[i]; + else js = binhead_multi[jgroup][jbin]; + + // 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 + + // 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 - if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + for (j = js; j >= 0; j = bins[j]) { + if(igroup != jgroup and j < i) continue; + + 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; } } - } else { - - // if different groups, implement with: - // loop over all atoms in jgroup bin - // 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 + + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; - js = binhead_multi[jgroup][jbin]; + 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); - for (j = js; j >= 0; j = bins[j]) { - if(j < i) continue; - - 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; - } + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + neighptr[n++] = j ^ mask_history; + else + neighptr[n++] = j; } } } diff --git a/src/npair_half_size_multi_newton_tri.cpp b/src/npair_half_size_multi_newton_tri.cpp index 0000f7b4d3..2000927c63 100755 --- a/src/npair_half_size_multi_newton_tri.cpp +++ b/src/npair_half_size_multi_newton_tri.cpp @@ -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 diff --git a/src/nstencil.cpp b/src/nstencil.cpp index acd7272245..3d86b57b9e 100644 --- a/src/nstencil.cpp +++ b/src/nstencil.cpp @@ -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 @@ -76,7 +76,7 @@ NStencil::NStencil(LAMMPS *lmp) : Pointers(lmp) distsq_multi_old = nullptr; nstencil_multi = nullptr; - stencil_multi = nullptr; + stencil_multi = nullptr; maxstencil_multi = nullptr; flag_half_multi = nullptr; @@ -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; } diff --git a/src/nstencil.h b/src/nstencil.h index b06f340d6f..b38276b45f 100644 --- a/src/nstencil.h +++ b/src/nstencil.h @@ -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 *);