2792 lines
88 KiB
C++
2792 lines
88 KiB
C++
// clang-format off
|
|
/* ----------------------------------------------------------------------
|
|
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
|
https://www.lammps.org/, 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.
|
|
------------------------------------------------------------------------- */
|
|
|
|
/* ----------------------------------------------------------------------
|
|
Contributing author (triclinic and multi-neigh) : Pieter in 't Veld (SNL)
|
|
Contributing author (improved multi-neigh) : Joel Clemmer (SNL)
|
|
------------------------------------------------------------------------- */
|
|
|
|
#include "neighbor.h"
|
|
|
|
#include "atom.h"
|
|
#include "atom_vec.h"
|
|
#include "citeme.h"
|
|
#include "comm.h"
|
|
#include "compute.h"
|
|
#include "domain.h"
|
|
#include "error.h"
|
|
#include "fix.h"
|
|
#include "force.h"
|
|
#include "group.h"
|
|
#include "memory.h"
|
|
#include "modify.h"
|
|
#include "nbin.h"
|
|
#include "neigh_list.h"
|
|
#include "neigh_request.h"
|
|
#include "npair.h"
|
|
#include "nstencil.h"
|
|
#include "ntopo.h"
|
|
#include "output.h"
|
|
#include "pair.h"
|
|
#include "pair_hybrid.h"
|
|
#include "respa.h"
|
|
#include "style_nbin.h" // IWYU pragma: keep
|
|
#include "style_npair.h" // IWYU pragma: keep
|
|
#include "style_nstencil.h" // IWYU pragma: keep
|
|
#include "style_ntopo.h" // IWYU pragma: keep
|
|
#include "tokenizer.h"
|
|
#include "update.h"
|
|
|
|
#include <cmath>
|
|
#include <cstring>
|
|
|
|
using namespace LAMMPS_NS;
|
|
using namespace NeighConst;
|
|
|
|
#define RQDELTA 1
|
|
#define EXDELTA 1
|
|
#define DELTA_PERATOM 64
|
|
|
|
#define BIG 1.0e20
|
|
|
|
enum{NONE,ALL,PARTIAL,TEMPLATE};
|
|
|
|
static const char cite_neigh_multi_old[] =
|
|
"neighbor multi/old command: doi:10.1016/j.cpc.2008.03.005\n\n"
|
|
"@Article{Intveld08,\n"
|
|
" author = {P.{\\,}J.~in{\\,}'t~Veld and S.{\\,}J.~Plimpton"
|
|
" and G.{\\,}S.~Grest},\n"
|
|
" title = {Accurate and Efficient Methods for Modeling Colloidal\n"
|
|
" Mixtures in an Explicit Solvent using Molecular Dynamics},\n"
|
|
" journal = {Comp.~Phys.~Comm.},\n"
|
|
" year = 2008,\n"
|
|
" volume = 179,\n"
|
|
" pages = {320--329}\n"
|
|
"}\n\n";
|
|
|
|
static const char cite_neigh_multi[] =
|
|
"neighbor multi command: doi:10.1016/j.cpc.2008.03.005, doi:10.1007/s40571-020-00361-2\n\n"
|
|
"@Article{Intveld08,\n"
|
|
" author = {P.{\\,}J.~in{\\,}'t~Veld and S.{\\,}J.~Plimpton"
|
|
" and G.{\\,}S.~Grest},\n"
|
|
" title = {Accurate and Efficient Methods for Modeling Colloidal\n"
|
|
" Mixtures in an Explicit Solvent using Molecular Dynamics},\n"
|
|
" journal = {Comp.~Phys.~Comm.},\n"
|
|
" year = 2008,\n"
|
|
" 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"
|
|
" detection applied to investigate the quasi-static limit},\n"
|
|
" journal = {Computational Particle Mechanics},\n"
|
|
" year = {2020}\n"
|
|
"}\n\n";
|
|
|
|
// template for factory functions:
|
|
// there will be one instance for each style keyword in the respective style_xxx.h files
|
|
|
|
template <typename S, typename T> static S *style_creator(LAMMPS *lmp)
|
|
{
|
|
return new T(lmp);
|
|
}
|
|
|
|
//#define NEIGH_LIST_DEBUG 1
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp),
|
|
pairclass(nullptr), pairnames(nullptr), pairmasks(nullptr)
|
|
{
|
|
MPI_Comm_rank(world,&me);
|
|
MPI_Comm_size(world,&nprocs);
|
|
|
|
firsttime = 1;
|
|
|
|
style = Neighbor::BIN;
|
|
every = 1;
|
|
delay = 10;
|
|
dist_check = 1;
|
|
pgsize = 100000;
|
|
oneatom = 2000;
|
|
binsizeflag = 0;
|
|
build_once = 0;
|
|
cluster_check = 0;
|
|
ago = -1;
|
|
|
|
cutneighmax = 0.0;
|
|
cutneighsq = nullptr;
|
|
cutneighghostsq = nullptr;
|
|
cuttype = nullptr;
|
|
cuttypesq = nullptr;
|
|
fixchecklist = nullptr;
|
|
|
|
// pairwise neighbor lists and associated data structs
|
|
|
|
nlist = 0;
|
|
lists = nullptr;
|
|
|
|
nbin = 0;
|
|
neigh_bin = nullptr;
|
|
|
|
nstencil = 0;
|
|
neigh_stencil = nullptr;
|
|
|
|
neigh_pair = nullptr;
|
|
|
|
nstencil_perpetual = 0;
|
|
slist = nullptr;
|
|
|
|
npair_perpetual = 0;
|
|
plist = nullptr;
|
|
|
|
nrequest = maxrequest = 0;
|
|
requests = nullptr;
|
|
|
|
old_nrequest = 0;
|
|
old_requests = nullptr;
|
|
|
|
old_style = style;
|
|
old_triclinic = 0;
|
|
old_pgsize = pgsize;
|
|
old_oneatom = oneatom;
|
|
|
|
binclass = nullptr;
|
|
binnames = nullptr;
|
|
binmasks = nullptr;
|
|
stencilclass = nullptr;
|
|
stencilnames = nullptr;
|
|
stencilmasks = nullptr;
|
|
|
|
// topology lists
|
|
|
|
bondwhich = anglewhich = dihedralwhich = improperwhich = NONE;
|
|
|
|
neigh_bond = nullptr;
|
|
neigh_angle = nullptr;
|
|
neigh_dihedral = nullptr;
|
|
neigh_improper = nullptr;
|
|
|
|
// coords at last neighboring
|
|
|
|
maxhold = 0;
|
|
xhold = nullptr;
|
|
lastcall = -1;
|
|
last_setup_bins = -1;
|
|
|
|
// pair exclusion list info
|
|
|
|
includegroup = 0;
|
|
|
|
nex_type = maxex_type = 0;
|
|
ex1_type = ex2_type = nullptr;
|
|
ex_type = nullptr;
|
|
|
|
nex_group = maxex_group = 0;
|
|
ex1_group = ex2_group = ex1_bit = ex2_bit = nullptr;
|
|
|
|
nex_mol = maxex_mol = 0;
|
|
ex_mol_group = ex_mol_bit = ex_mol_intra = nullptr;
|
|
|
|
// Multi data
|
|
|
|
type2collection = nullptr;
|
|
collection2cut = nullptr;
|
|
collection = nullptr;
|
|
cutcollectionsq = nullptr;
|
|
custom_collection_flag = 0;
|
|
interval_collection_flag = 0;
|
|
nmax_collection = 0;
|
|
|
|
// Kokkos setting
|
|
|
|
copymode = 0;
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
Neighbor::~Neighbor()
|
|
{
|
|
if (copymode) return;
|
|
|
|
memory->destroy(cutneighsq);
|
|
memory->destroy(cutneighghostsq);
|
|
delete[] cuttype;
|
|
delete[] cuttypesq;
|
|
delete[] fixchecklist;
|
|
|
|
for (int i = 0; i < nlist; i++) delete lists[i];
|
|
for (int i = 0; i < nbin; i++) delete neigh_bin[i];
|
|
for (int i = 0; i < nstencil; i++) delete neigh_stencil[i];
|
|
for (int i = 0; i < nlist; i++) delete neigh_pair[i];
|
|
delete[] lists;
|
|
delete[] neigh_bin;
|
|
delete[] neigh_stencil;
|
|
delete[] neigh_pair;
|
|
|
|
delete[] slist;
|
|
delete[] plist;
|
|
|
|
for (int i = 0; i < nrequest; i++)
|
|
if (requests[i]) delete requests[i];
|
|
memory->sfree(requests);
|
|
for (int i = 0; i < old_nrequest; i++)
|
|
if (old_requests[i]) delete old_requests[i];
|
|
memory->sfree(old_requests);
|
|
|
|
delete[] binclass;
|
|
delete[] binnames;
|
|
delete[] binmasks;
|
|
delete[] stencilclass;
|
|
delete[] stencilnames;
|
|
delete[] stencilmasks;
|
|
delete[] pairclass;
|
|
delete[] pairnames;
|
|
delete[] pairmasks;
|
|
|
|
delete neigh_bond;
|
|
delete neigh_angle;
|
|
delete neigh_dihedral;
|
|
delete neigh_improper;
|
|
|
|
memory->destroy(xhold);
|
|
|
|
memory->destroy(ex1_type);
|
|
memory->destroy(ex2_type);
|
|
memory->destroy(ex_type);
|
|
|
|
memory->destroy(ex1_group);
|
|
memory->destroy(ex2_group);
|
|
delete[] ex1_bit;
|
|
delete[] ex2_bit;
|
|
|
|
memory->destroy(ex_mol_group);
|
|
delete[] ex_mol_bit;
|
|
memory->destroy(ex_mol_intra);
|
|
|
|
memory->destroy(type2collection);
|
|
memory->destroy(collection2cut);
|
|
memory->destroy(collection);
|
|
memory->destroy(cutcollectionsq);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
void Neighbor::init()
|
|
{
|
|
int i,j,n;
|
|
|
|
ncalls = ndanger = 0;
|
|
dimension = domain->dimension;
|
|
triclinic = domain->triclinic;
|
|
newton_pair = force->newton_pair;
|
|
|
|
// error check
|
|
|
|
if (delay > 0 && (delay % every) != 0)
|
|
error->all(FLERR,"Neighbor delay must be 0 or multiple of every setting");
|
|
|
|
if (pgsize < 10*oneatom)
|
|
error->all(FLERR,"Neighbor page size must be >= 10x the one atom setting");
|
|
|
|
// ------------------------------------------------------------------
|
|
// settings
|
|
|
|
// bbox lo/hi ptrs = bounding box of entire domain, stored by Domain
|
|
|
|
if (triclinic == 0) {
|
|
bboxlo = domain->boxlo;
|
|
bboxhi = domain->boxhi;
|
|
} else {
|
|
bboxlo = domain->boxlo_bound;
|
|
bboxhi = domain->boxhi_bound;
|
|
}
|
|
|
|
// set neighbor cutoffs (force cutoff + skin)
|
|
// trigger determines when atoms migrate and neighbor lists are rebuilt
|
|
// needs to be non-zero for migration distance check
|
|
// even if pair = nullptr and no neighbor lists are used
|
|
// cutneigh = force cutoff + skin if cutforce > 0, else cutneigh = 0
|
|
// cutneighghost = pair cutghost if it requests it, else same as cutneigh
|
|
|
|
triggersq = 0.25*skin*skin;
|
|
boxcheck = 0;
|
|
if (domain->box_change && (domain->xperiodic || domain->yperiodic ||
|
|
(dimension == 3 && domain->zperiodic)))
|
|
boxcheck = 1;
|
|
|
|
n = atom->ntypes;
|
|
if (cutneighsq == nullptr) {
|
|
if (lmp->kokkos) init_cutneighsq_kokkos(n);
|
|
else memory->create(cutneighsq,n+1,n+1,"neigh:cutneighsq");
|
|
memory->create(cutneighghostsq,n+1,n+1,"neigh:cutneighghostsq");
|
|
cuttype = new double[n+1];
|
|
cuttypesq = new double[n+1];
|
|
}
|
|
|
|
double cutoff,delta,cut;
|
|
cutneighmin = BIG;
|
|
cutneighmax = 0.0;
|
|
|
|
for (i = 1; i <= n; i++) {
|
|
cuttype[i] = cuttypesq[i] = 0.0;
|
|
for (j = 1; j <= n; j++) {
|
|
if (force->pair) cutoff = sqrt(force->pair->cutsq[i][j]);
|
|
else cutoff = 0.0;
|
|
if (cutoff > 0.0) delta = skin;
|
|
else delta = 0.0;
|
|
cut = cutoff + delta;
|
|
|
|
cutneighsq[i][j] = cut*cut;
|
|
cuttype[i] = MAX(cuttype[i],cut);
|
|
cuttypesq[i] = MAX(cuttypesq[i],cut*cut);
|
|
cutneighmin = MIN(cutneighmin,cut);
|
|
cutneighmax = MAX(cutneighmax,cut);
|
|
|
|
if (force->pair && force->pair->ghostneigh) {
|
|
cut = force->pair->cutghost[i][j] + skin;
|
|
cutneighghostsq[i][j] = cut*cut;
|
|
} else cutneighghostsq[i][j] = cut*cut;
|
|
}
|
|
}
|
|
cutneighmaxsq = cutneighmax * cutneighmax;
|
|
|
|
// Define cutoffs for multi
|
|
if (style == Neighbor::MULTI) {
|
|
int icollection, jcollection;
|
|
|
|
// If collections not yet defined, create default map using types
|
|
if (!custom_collection_flag) {
|
|
ncollections = n;
|
|
interval_collection_flag = 0;
|
|
if (!type2collection)
|
|
memory->create(type2collection,n+1,"neigh:type2collection");
|
|
for (i = 1; i <= n; i++)
|
|
type2collection[i] = i-1;
|
|
}
|
|
|
|
memory->grow(cutcollectionsq, ncollections, ncollections, "neigh:cutcollectionsq");
|
|
|
|
// 3 possible ways of defining collections
|
|
// 1) Types are used to define collections
|
|
// Each collection loops through its owned types, and uses cutneighsq to calculate its cutoff
|
|
// 2) Collections are defined by intervals, point particles
|
|
// Types are first sorted into collections based on cutneighsq[i][i]
|
|
// Each collection loops through its owned types, and uses cutneighsq to calculate its cutoff
|
|
// 3) Collections are defined by intervals, finite particles
|
|
//
|
|
|
|
// Define collection cutoffs
|
|
for (i = 0; i < ncollections; i++)
|
|
for (j = 0; j < ncollections; j++)
|
|
cutcollectionsq[i][j] = 0.0;
|
|
|
|
if (!interval_collection_flag) {
|
|
finite_cut_flag = 0;
|
|
for (i = 1; i <= n; i++){
|
|
icollection = type2collection[i];
|
|
for (j = 1; j <= n; j++){
|
|
jcollection = type2collection[j];
|
|
if (cutneighsq[i][j] > cutcollectionsq[icollection][jcollection]) {
|
|
cutcollectionsq[icollection][jcollection] = cutneighsq[i][j];
|
|
cutcollectionsq[jcollection][icollection] = cutneighsq[i][j];
|
|
}
|
|
}
|
|
}
|
|
} else {
|
|
if (force->pair->finitecutflag) {
|
|
finite_cut_flag = 1;
|
|
// If cutoffs depend on finite atom sizes, use radii of intervals to find cutoffs
|
|
double ri, rj, tmp;
|
|
for (i = 0; i < ncollections; i++){
|
|
ri = collection2cut[i]*0.5;
|
|
for (j = 0; j < ncollections; j++){
|
|
rj = collection2cut[j]*0.5;
|
|
tmp = force->pair->radii2cut(ri, rj) + skin;
|
|
cutcollectionsq[i][j] = tmp*tmp;
|
|
}
|
|
}
|
|
} else {
|
|
finite_cut_flag = 0;
|
|
|
|
// Map types to collections
|
|
if (!type2collection)
|
|
memory->create(type2collection,n+1,"neigh:type2collection");
|
|
|
|
for (i = 1; i <= n; i++)
|
|
type2collection[i] = -1;
|
|
|
|
double cuttmp;
|
|
for (i = 1; i <= n; i++){
|
|
// Remove skin added to cutneighsq
|
|
cuttmp = sqrt(cutneighsq[i][i]) - skin;
|
|
for (icollection = 0; icollection < ncollections; icollection ++){
|
|
if (collection2cut[icollection] >= cuttmp) {
|
|
type2collection[i] = icollection;
|
|
break;
|
|
}
|
|
}
|
|
|
|
if (type2collection[i] == -1)
|
|
error->all(FLERR, "Pair cutoff exceeds interval cutoffs for multi");
|
|
}
|
|
|
|
// Define cutoffs
|
|
for (i = 1; i <= n; i++){
|
|
icollection = type2collection[i];
|
|
for (j = 1; j <= n; j++){
|
|
jcollection = type2collection[j];
|
|
if (cutneighsq[i][j] > cutcollectionsq[icollection][jcollection]) {
|
|
cutcollectionsq[icollection][jcollection] = cutneighsq[i][j];
|
|
cutcollectionsq[jcollection][icollection] = cutneighsq[i][j];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// rRESPA cutoffs
|
|
|
|
int respa = 0;
|
|
if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) {
|
|
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
|
|
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
|
|
}
|
|
|
|
if (respa) {
|
|
double *cut_respa = ((Respa *) update->integrate)->cutoff;
|
|
cut_inner_sq = (cut_respa[1] + skin) * (cut_respa[1] + skin);
|
|
cut_middle_sq = (cut_respa[3] + skin) * (cut_respa[3] + skin);
|
|
cut_middle_inside_sq = (cut_respa[0] - skin) * (cut_respa[0] - skin);
|
|
if (cut_respa[0]-skin < 0) cut_middle_inside_sq = 0.0;
|
|
}
|
|
|
|
// fixchecklist = other classes that can induce reneighboring in decide()
|
|
|
|
restart_check = 0;
|
|
if (output->restart_flag) restart_check = 1;
|
|
|
|
delete[] fixchecklist;
|
|
fixchecklist = nullptr;
|
|
fixchecklist = new int[modify->nfix];
|
|
|
|
fix_check = 0;
|
|
for (i = 0; i < modify->nfix; i++)
|
|
if (modify->fix[i]->force_reneighbor)
|
|
fixchecklist[fix_check++] = i;
|
|
|
|
must_check = 0;
|
|
if (restart_check || fix_check) must_check = 1;
|
|
|
|
// set special_flag for 1-2, 1-3, 1-4 neighbors
|
|
// flag[0] is not used, flag[1] = 1-2, flag[2] = 1-3, flag[3] = 1-4
|
|
// flag = 0 if both LJ/Coulomb special values are 0.0
|
|
// flag = 1 if both LJ/Coulomb special values are 1.0
|
|
// flag = 2 otherwise or if KSpace solver is enabled
|
|
// pairwise portion of KSpace solver uses all 1-2,1-3,1-4 neighbors
|
|
// or selected Coulomb-approixmation pair styles require it
|
|
|
|
if (force->special_lj[1] == 0.0 && force->special_coul[1] == 0.0)
|
|
special_flag[1] = 0;
|
|
else if (force->special_lj[1] == 1.0 && force->special_coul[1] == 1.0)
|
|
special_flag[1] = 1;
|
|
else special_flag[1] = 2;
|
|
|
|
if (force->special_lj[2] == 0.0 && force->special_coul[2] == 0.0)
|
|
special_flag[2] = 0;
|
|
else if (force->special_lj[2] == 1.0 && force->special_coul[2] == 1.0)
|
|
special_flag[2] = 1;
|
|
else special_flag[2] = 2;
|
|
|
|
if (force->special_lj[3] == 0.0 && force->special_coul[3] == 0.0)
|
|
special_flag[3] = 0;
|
|
else if (force->special_lj[3] == 1.0 && force->special_coul[3] == 1.0)
|
|
special_flag[3] = 1;
|
|
else special_flag[3] = 2;
|
|
|
|
// We cannot remove special neighbors with kspace or kspace-like pair styles
|
|
// as the exclusion needs to remove the full coulomb and not the damped interaction.
|
|
// Special treatment is required for hybrid pair styles since Force::pair_match()
|
|
// will only return a non-null pointer if there is only one substyle of the kind.
|
|
|
|
if (force->kspace) {
|
|
special_flag[1] = special_flag[2] = special_flag[3] = 2;
|
|
} else {
|
|
PairHybrid *ph = reinterpret_cast<PairHybrid *>(force->pair_match("^hybrid",0));
|
|
if (ph) {
|
|
int flag=0;
|
|
for (int isub=0; isub < ph->nstyles; ++isub) {
|
|
if (force->pair_match("coul/wolf",0,isub)
|
|
|| force->pair_match("coul/dsf",0,isub)
|
|
|| force->pair_match("coul/exclude",0)
|
|
|| force->pair_match("thole",0,isub))
|
|
++flag;
|
|
}
|
|
if (flag)
|
|
special_flag[1] = special_flag[2] = special_flag[3] = 2;
|
|
} else {
|
|
if (force->pair_match("coul/wolf",0)
|
|
|| force->pair_match("coul/dsf",0)
|
|
|| force->pair_match("coul/exclude",0)
|
|
|| force->pair_match("thole",0))
|
|
special_flag[1] = special_flag[2] = special_flag[3] = 2;
|
|
}
|
|
}
|
|
|
|
// ------------------------------------------------------------------
|
|
// xhold array
|
|
|
|
// free if not needed for this run
|
|
|
|
if (dist_check == 0) {
|
|
memory->destroy(xhold);
|
|
maxhold = 0;
|
|
xhold = nullptr;
|
|
}
|
|
|
|
// first time allocation
|
|
|
|
if (dist_check) {
|
|
if (maxhold == 0) {
|
|
maxhold = atom->nmax;
|
|
memory->create(xhold,maxhold,3,"neigh:xhold");
|
|
}
|
|
}
|
|
|
|
// ------------------------------------------------------------------
|
|
// exclusion lists
|
|
|
|
// depend on type, group, molecule settings from neigh_modify
|
|
// warn if exclusions used with KSpace solver
|
|
|
|
n = atom->ntypes;
|
|
|
|
if (nex_type == 0 && nex_group == 0 && nex_mol == 0) exclude = 0;
|
|
else exclude = 1;
|
|
|
|
if (nex_type) {
|
|
if (lmp->kokkos)
|
|
init_ex_type_kokkos(n);
|
|
else {
|
|
memory->destroy(ex_type);
|
|
memory->create(ex_type,n+1,n+1,"neigh:ex_type");
|
|
}
|
|
|
|
for (i = 1; i <= n; i++)
|
|
for (j = 1; j <= n; j++)
|
|
ex_type[i][j] = 0;
|
|
|
|
for (i = 0; i < nex_type; i++) {
|
|
if (ex1_type[i] <= 0 || ex1_type[i] > n ||
|
|
ex2_type[i] <= 0 || ex2_type[i] > n)
|
|
error->all(FLERR,"Invalid atom type in neighbor exclusion list");
|
|
ex_type[ex1_type[i]][ex2_type[i]] = 1;
|
|
ex_type[ex2_type[i]][ex1_type[i]] = 1;
|
|
}
|
|
}
|
|
|
|
if (nex_group) {
|
|
if (lmp->kokkos)
|
|
init_ex_bit_kokkos();
|
|
else {
|
|
delete[] ex1_bit;
|
|
delete[] ex2_bit;
|
|
ex1_bit = new int[nex_group];
|
|
ex2_bit = new int[nex_group];
|
|
}
|
|
|
|
for (i = 0; i < nex_group; i++) {
|
|
ex1_bit[i] = group->bitmask[ex1_group[i]];
|
|
ex2_bit[i] = group->bitmask[ex2_group[i]];
|
|
}
|
|
}
|
|
|
|
if (nex_mol) {
|
|
if (lmp->kokkos)
|
|
init_ex_mol_bit_kokkos();
|
|
else {
|
|
delete[] ex_mol_bit;
|
|
ex_mol_bit = new int[nex_mol];
|
|
}
|
|
|
|
for (i = 0; i < nex_mol; i++)
|
|
ex_mol_bit[i] = group->bitmask[ex_mol_group[i]];
|
|
}
|
|
|
|
if (exclude && force->kspace && me == 0)
|
|
error->warning(FLERR,"Neighbor exclusions used with KSpace solver "
|
|
"may give inconsistent Coulombic energies");
|
|
|
|
if (lmp->kokkos)
|
|
set_binsize_kokkos();
|
|
|
|
// ------------------------------------------------------------------
|
|
// create pairwise lists
|
|
// one-time call to init_styles() to scan style files and setup
|
|
// init_pair() creates auxiliary classes: NBin, NStencil, NPair
|
|
|
|
if (firsttime) init_styles();
|
|
firsttime = 0;
|
|
|
|
int same = init_pair();
|
|
|
|
// invoke copy_neighbor_info() in Bin,Stencil,Pair classes
|
|
// copied once per run in case any cutoff, exclusion, special info changed
|
|
|
|
for (i = 0; i < nbin; i++) neigh_bin[i]->copy_neighbor_info();
|
|
for (i = 0; i < nstencil; i++) neigh_stencil[i]->copy_neighbor_info();
|
|
for (i = 0; i < nlist; i++)
|
|
if (neigh_pair[i]) neigh_pair[i]->copy_neighbor_info();
|
|
|
|
if (!same && comm->me == 0) print_pairwise_info();
|
|
|
|
// can now delete requests so next run can make new ones
|
|
// print_pairwise_info() made use of requests
|
|
// set of NeighLists now stores all needed info
|
|
|
|
for (i = 0; i < nrequest; i++) {
|
|
delete requests[i];
|
|
requests[i] = nullptr;
|
|
}
|
|
nrequest = 0;
|
|
|
|
// ------------------------------------------------------------------
|
|
// create topology lists
|
|
// instantiated topo styles can change from run to run
|
|
|
|
init_topology();
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
create and initialize lists of Nbin, Nstencil, NPair classes
|
|
lists have info on all classes in 3 style*.h files
|
|
cannot do this in constructor, b/c too early to instantiate classes
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::init_styles()
|
|
{
|
|
// extract info from NBin classes listed in style_nbin.h
|
|
|
|
nbclass = 0;
|
|
|
|
#define NBIN_CLASS
|
|
#define NBinStyle(key,Class,bitmasks) nbclass++;
|
|
#include "style_nbin.h" // IWYU pragma: keep
|
|
#undef NBinStyle
|
|
#undef NBIN_CLASS
|
|
|
|
binclass = new BinCreator[nbclass];
|
|
binnames = new char*[nbclass];
|
|
binmasks = new int[nbclass];
|
|
nbclass = 0;
|
|
|
|
#define NBIN_CLASS
|
|
#define NBinStyle(key,Class,bitmasks) \
|
|
binnames[nbclass] = (char *) #key; \
|
|
binclass[nbclass] = &style_creator<NBin, Class>; \
|
|
binmasks[nbclass++] = bitmasks;
|
|
#include "style_nbin.h" // IWYU pragma: keep
|
|
#undef NBinStyle
|
|
#undef NBIN_CLASS
|
|
|
|
// extract info from NStencil classes listed in style_nstencil.h
|
|
|
|
nsclass = 0;
|
|
|
|
#define NSTENCIL_CLASS
|
|
#define NStencilStyle(key,Class,bitmasks) nsclass++;
|
|
#include "style_nstencil.h" // IWYU pragma: keep
|
|
#undef NStencilStyle
|
|
#undef NSTENCIL_CLASS
|
|
|
|
stencilclass = new StencilCreator[nsclass];
|
|
stencilnames = new char*[nsclass];
|
|
stencilmasks = new int[nsclass];
|
|
nsclass = 0;
|
|
|
|
#define NSTENCIL_CLASS
|
|
#define NStencilStyle(key,Class,bitmasks) \
|
|
stencilnames[nsclass] = (char *) #key; \
|
|
stencilclass[nsclass] = &style_creator<NStencil, Class>; \
|
|
stencilmasks[nsclass++] = bitmasks;
|
|
#include "style_nstencil.h" // IWYU pragma: keep
|
|
#undef NStencilStyle
|
|
#undef NSTENCIL_CLASS
|
|
|
|
// extract info from NPair classes listed in style_npair.h
|
|
|
|
npclass = 0;
|
|
|
|
#define NPAIR_CLASS
|
|
#define NPairStyle(key,Class,bitmasks) npclass++;
|
|
#include "style_npair.h" // IWYU pragma: keep
|
|
#undef NPairStyle
|
|
#undef NPAIR_CLASS
|
|
|
|
pairclass = new PairCreator[npclass];
|
|
pairnames = new char*[npclass];
|
|
pairmasks = new int[npclass];
|
|
npclass = 0;
|
|
|
|
#define NPAIR_CLASS
|
|
#define NPairStyle(key,Class,bitmasks) \
|
|
pairnames[npclass] = (char *) #key; \
|
|
pairclass[npclass] = &style_creator<NPair, Class>; \
|
|
pairmasks[npclass++] = bitmasks;
|
|
#include "style_npair.h" // IWYU pragma: keep
|
|
#undef NPairStyle
|
|
#undef NPAIR_CLASS
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
create and initialize NPair classes
|
|
------------------------------------------------------------------------- */
|
|
|
|
int Neighbor::init_pair()
|
|
{
|
|
int i,j,k,m;
|
|
|
|
// test if pairwise lists need to be re-created
|
|
// no need to re-create if:
|
|
// neigh style, triclinic, pgsize, oneatom have not changed
|
|
// current requests = old requests
|
|
// so just return:
|
|
// delete requests so next run can make new ones
|
|
// current set of NeighLists already stores all needed info
|
|
// requests are compared via identical() before:
|
|
// any requests are morphed using logic below
|
|
// any requests are added below, e.g. as parents of pair hybrid skip lists
|
|
// copy them via requests_new2old() BEFORE any changes made to requests
|
|
// necessary b/c morphs can change requestor settings (see comment below)
|
|
|
|
int same = 1;
|
|
if (style != old_style) same = 0;
|
|
if (triclinic != old_triclinic) same = 0;
|
|
if (pgsize != old_pgsize) same = 0;
|
|
if (oneatom != old_oneatom) same = 0;
|
|
|
|
if (nrequest != old_nrequest) same = 0;
|
|
else
|
|
for (i = 0; i < nrequest; i++)
|
|
if (requests[i]->identical(old_requests[i]) == 0) same = 0;
|
|
|
|
#ifdef NEIGH_LIST_DEBUG
|
|
if (comm->me == 0) printf("SAME flag %d\n",same);
|
|
#endif
|
|
|
|
if (same) return same;
|
|
requests_new2old();
|
|
|
|
// delete old lists since creating new ones
|
|
|
|
for (i = 0; i < nlist; i++) delete lists[i];
|
|
for (i = 0; i < nbin; i++) delete neigh_bin[i];
|
|
for (i = 0; i < nstencil; i++) delete neigh_stencil[i];
|
|
for (i = 0; i < nlist; i++) delete neigh_pair[i];
|
|
delete[] lists;
|
|
delete[] neigh_bin;
|
|
delete[] neigh_stencil;
|
|
delete[] neigh_pair;
|
|
|
|
// error check on requests
|
|
// do not allow occasional, ghost, bin list
|
|
// b/c it still uses variant of coord2bin() in NPair() method
|
|
// instead of atom2bin, this could cause error b/c stoms have
|
|
// moved out of proc domain by time occasional list is built
|
|
// solution would be to use a different NBin variant
|
|
// that used Npair::coord2bin(x,ix,iy,iz) (then delete it from NPair)
|
|
// and stored the ix,iy,iz values for all atoms (including ghosts)
|
|
// at time of binning when neighbor lists are rebuilt,
|
|
// similar to what vanilla Nbin::coord2atom() does now in atom2bin
|
|
|
|
if (style == Neighbor::BIN) {
|
|
for (i = 0; i < nrequest; i++)
|
|
if (requests[i]->occasional && requests[i]->ghost)
|
|
error->all(FLERR,"Cannot request an occasional binned neighbor list "
|
|
"with ghost info");
|
|
}
|
|
|
|
// morph requests in various ways
|
|
// purpose is to avoid duplicate or inefficient builds
|
|
// may add new requests if a needed request to derive from does not exist
|
|
// methods:
|
|
// (1) unique = create unique lists if cutoff is explicitly set
|
|
// (2) skip = create any new non-skip lists needed by pair hybrid skip lists
|
|
// (3) granular = adjust parent and skip lists for granular onesided usage
|
|
// (4) h/f = pair up any matching half/full lists
|
|
// (5) copy = convert as many lists as possible to copy lists
|
|
// order of morph methods matters:
|
|
// (3) after (2), b/c it adjusts lists created by (2)
|
|
// (4) after (2) and (3),
|
|
// b/c (2) may create new full lists, (3) may change them
|
|
// (5) last, after all lists are finalized, so all possible copies found
|
|
|
|
int nrequest_original = nrequest;
|
|
|
|
morph_unique();
|
|
morph_skip();
|
|
morph_granular(); // this method can change flags set by requestor
|
|
morph_halffull();
|
|
morph_copy();
|
|
|
|
// create new lists, one per request including added requests
|
|
// wait to allocate initial pages until copy lists are detected
|
|
// NOTE: can I allocate now, instead of down below?
|
|
|
|
nlist = nrequest;
|
|
|
|
lists = new NeighList*[nrequest];
|
|
neigh_bin = new NBin*[nrequest];
|
|
neigh_stencil = new NStencil*[nrequest];
|
|
neigh_pair = new NPair*[nrequest];
|
|
|
|
// allocate new lists
|
|
// pass list ptr back to requestor (except for Command class)
|
|
// only for original requests, not ones added by Neighbor class
|
|
|
|
for (i = 0; i < nrequest; i++) {
|
|
if (requests[i]->kokkos_host || requests[i]->kokkos_device)
|
|
create_kokkos_list(i);
|
|
else lists[i] = new NeighList(lmp);
|
|
lists[i]->index = i;
|
|
lists[i]->requestor = requests[i]->requestor;
|
|
|
|
if (requests[i]->pair) {
|
|
lists[i]->requestor_type = NeighList::PAIR;
|
|
} else if (requests[i]->fix) {
|
|
lists[i]->requestor_type = NeighList::FIX;
|
|
} else if (requests[i]->compute) {
|
|
lists[i]->requestor_type = NeighList::COMPUTE;
|
|
}
|
|
|
|
if (requests[i]->pair && i < nrequest_original) {
|
|
Pair *pair = (Pair *) requests[i]->requestor;
|
|
pair->init_list(requests[i]->id,lists[i]);
|
|
} else if (requests[i]->fix && i < nrequest_original) {
|
|
Fix *fix = (Fix *) requests[i]->requestor;
|
|
fix->init_list(requests[i]->id,lists[i]);
|
|
} else if (requests[i]->compute && i < nrequest_original) {
|
|
Compute *compute = (Compute *) requests[i]->requestor;
|
|
compute->init_list(requests[i]->id,lists[i]);
|
|
}
|
|
}
|
|
|
|
// invoke post_constructor() for all lists
|
|
// copies info from requests to lists, sets ptrs to related lists
|
|
|
|
for (i = 0; i < nrequest; i++)
|
|
lists[i]->post_constructor(requests[i]);
|
|
|
|
// assign Bin,Stencil,Pair style to each list
|
|
|
|
int flag;
|
|
for (i = 0; i < nrequest; i++) {
|
|
flag = choose_bin(requests[i]);
|
|
lists[i]->bin_method = flag;
|
|
if (flag < 0)
|
|
error->all(FLERR,"Requested neighbor bin option does not exist");
|
|
|
|
flag = choose_stencil(requests[i]);
|
|
lists[i]->stencil_method = flag;
|
|
if (flag < 0)
|
|
error->all(FLERR,"Requested neighbor stencil method does not exist");
|
|
|
|
flag = choose_pair(requests[i]);
|
|
lists[i]->pair_method = flag;
|
|
if (flag < 0)
|
|
error->all(FLERR,"Requested neighbor pair method does not exist");
|
|
}
|
|
|
|
// instantiate unique Bin,Stencil classes in neigh_bin & neigh_stencil vecs
|
|
// unique = only one of its style, or request unique flag set (custom cutoff)
|
|
|
|
nbin = 0;
|
|
for (i = 0; i < nrequest; i++) {
|
|
requests[i]->index_bin = -1;
|
|
flag = lists[i]->bin_method;
|
|
if (flag == 0) continue;
|
|
if (!requests[i]->unique) {
|
|
for (j = 0; j < nbin; j++)
|
|
if (neigh_bin[j]->istyle == flag &&
|
|
neigh_bin[j]->cutoff_custom == 0.0) break;
|
|
if (j < nbin) {
|
|
requests[i]->index_bin = j;
|
|
continue;
|
|
}
|
|
}
|
|
|
|
BinCreator &bin_creator = binclass[flag-1];
|
|
neigh_bin[nbin] = bin_creator(lmp);
|
|
neigh_bin[nbin]->post_constructor(requests[i]);
|
|
neigh_bin[nbin]->istyle = flag;
|
|
|
|
requests[i]->index_bin = nbin;
|
|
nbin++;
|
|
}
|
|
|
|
nstencil = 0;
|
|
for (i = 0; i < nrequest; i++) {
|
|
requests[i]->index_stencil = -1;
|
|
flag = lists[i]->stencil_method;
|
|
if (flag == 0) continue;
|
|
if (!requests[i]->unique) {
|
|
for (j = 0; j < nstencil; j++)
|
|
if (neigh_stencil[j]->istyle == flag &&
|
|
neigh_stencil[j]->cutoff_custom == 0.0) break;
|
|
if (j < nstencil) {
|
|
requests[i]->index_stencil = j;
|
|
continue;
|
|
}
|
|
}
|
|
|
|
StencilCreator &stencil_creator = stencilclass[flag-1];
|
|
neigh_stencil[nstencil] = stencil_creator(lmp);
|
|
neigh_stencil[nstencil]->post_constructor(requests[i]);
|
|
neigh_stencil[nstencil]->istyle = flag;
|
|
|
|
if (lists[i]->bin_method > 0) {
|
|
neigh_stencil[nstencil]->nb = neigh_bin[requests[i]->index_bin];
|
|
if (neigh_stencil[nstencil]->nb == nullptr)
|
|
error->all(FLERR,"Could not assign bin method to neighbor stencil");
|
|
}
|
|
|
|
requests[i]->index_stencil = nstencil;
|
|
nstencil++;
|
|
}
|
|
|
|
// instantiate one Pair class per list in neigh_pair vec
|
|
|
|
for (i = 0; i < nrequest; i++) {
|
|
requests[i]->index_pair = -1;
|
|
flag = lists[i]->pair_method;
|
|
if (flag == 0) {
|
|
neigh_pair[i] = nullptr;
|
|
continue;
|
|
}
|
|
|
|
PairCreator &pair_creator = pairclass[flag-1];
|
|
lists[i]->np = neigh_pair[i] = pair_creator(lmp);
|
|
neigh_pair[i]->post_constructor(requests[i]);
|
|
neigh_pair[i]->istyle = flag;
|
|
|
|
if (lists[i]->bin_method > 0) {
|
|
neigh_pair[i]->nb = neigh_bin[requests[i]->index_bin];
|
|
if (neigh_pair[i]->nb == nullptr)
|
|
error->all(FLERR,"Could not assign bin method to neighbor pair");
|
|
}
|
|
if (lists[i]->stencil_method > 0) {
|
|
neigh_pair[i]->ns = neigh_stencil[requests[i]->index_stencil];
|
|
if (neigh_pair[i]->ns == nullptr)
|
|
error->all(FLERR,"Could not assign stencil method to neighbor pair");
|
|
}
|
|
|
|
requests[i]->index_pair = i;
|
|
}
|
|
|
|
// allocate initial pages for each list, except if copy flag set
|
|
|
|
for (i = 0; i < nlist; i++) {
|
|
if (lists[i]->copy && !lists[i]->kk2cpu)
|
|
continue;
|
|
lists[i]->setup_pages(pgsize,oneatom);
|
|
}
|
|
|
|
// first-time allocation of per-atom data for lists that are built and store
|
|
// lists that do not store: copy
|
|
// use atom->nmax for both grow() args
|
|
// i.e. grow first time to expanded size to avoid future reallocs
|
|
// also Kokkos list initialization
|
|
|
|
int maxatom = atom->nmax;
|
|
for (i = 0; i < nlist; i++) {
|
|
if (neigh_pair[i] && (!lists[i]->copy || lists[i]->kk2cpu))
|
|
lists[i]->grow(maxatom,maxatom);
|
|
}
|
|
|
|
// plist = indices of perpetual NPair classes
|
|
// perpetual = non-occasional, re-built at every reneighboring
|
|
// slist = indices of perpetual NStencil classes
|
|
// perpetual = used by any perpetual NPair class
|
|
|
|
delete[] slist;
|
|
delete[] plist;
|
|
nstencil_perpetual = npair_perpetual = 0;
|
|
slist = new int[nstencil];
|
|
plist = new int[nlist];
|
|
|
|
for (i = 0; i < nlist; i++) {
|
|
if (lists[i]->occasional == 0 && lists[i]->pair_method)
|
|
plist[npair_perpetual++] = i;
|
|
}
|
|
|
|
for (i = 0; i < nstencil; i++) {
|
|
flag = 0;
|
|
for (j = 0; j < npair_perpetual; j++)
|
|
if (lists[plist[j]]->stencil_method == neigh_stencil[i]->istyle)
|
|
flag = 1;
|
|
if (flag) slist[nstencil_perpetual++] = i;
|
|
}
|
|
|
|
// reorder plist vector if necessary
|
|
// relevant for lists that are derived from a parent list:
|
|
// half-full,copy,skip
|
|
// the child index must appear in plist after the parent index
|
|
// swap two indices within plist when dependency is mis-ordered
|
|
// start double loop check again whenever a swap is made
|
|
// done when entire double loop test results in no swaps
|
|
|
|
NeighList *ptr;
|
|
|
|
int done = 0;
|
|
while (!done) {
|
|
done = 1;
|
|
for (i = 0; i < npair_perpetual; i++) {
|
|
for (k = 0; k < 3; k++) {
|
|
ptr = nullptr;
|
|
if (k == 0) ptr = lists[plist[i]]->listcopy;
|
|
if (k == 1) ptr = lists[plist[i]]->listskip;
|
|
if (k == 2) ptr = lists[plist[i]]->listfull;
|
|
if (ptr == nullptr) continue;
|
|
for (m = 0; m < nrequest; m++)
|
|
if (ptr == lists[m]) break;
|
|
for (j = 0; j < npair_perpetual; j++)
|
|
if (m == plist[j]) break;
|
|
if (j < i) continue;
|
|
int tmp = plist[i]; // swap I,J indices
|
|
plist[i] = plist[j];
|
|
plist[j] = tmp;
|
|
done = 0;
|
|
break;
|
|
}
|
|
if (!done) break;
|
|
}
|
|
}
|
|
|
|
// debug output
|
|
|
|
#ifdef NEIGH_LIST_DEBUG
|
|
for (i = 0; i < nrequest; i++) lists[i]->print_attributes();
|
|
#endif
|
|
|
|
return same;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
scan NeighRequests to set additional flags:
|
|
custom cutoff lists and accelerator lists
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::morph_unique()
|
|
{
|
|
NeighRequest *irq;
|
|
|
|
for (int i = 0; i < nrequest; i++) {
|
|
irq = requests[i];
|
|
|
|
// if cut flag set by requestor, set unique flag
|
|
// this forces Pair,Stencil,Bin styles to be instantiated separately
|
|
|
|
if (irq->cut) irq->unique = 1;
|
|
|
|
// avoid flagging a neighbor list as both INTEL and OPENMP
|
|
|
|
if (irq->intel) irq->omp = 0;
|
|
|
|
// avoid flagging a neighbor list as both KOKKOS and INTEL or OPENMP
|
|
|
|
if (irq->kokkos_host || irq->kokkos_device) irq->omp = irq->intel = 0;
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
scan NeighRequests to process all skip lists
|
|
look for a matching non-skip list
|
|
if one exists, point at it via skiplist
|
|
else make new parent via copy_request() and point at it
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::morph_skip()
|
|
{
|
|
int i,j,inewton,jnewton;
|
|
NeighRequest *irq,*jrq,*nrq;
|
|
|
|
for (i = 0; i < nrequest; i++) {
|
|
irq = requests[i];
|
|
|
|
// only processing skip lists
|
|
|
|
if (!irq->skip) continue;
|
|
|
|
// these lists are created other ways, no need for skipping
|
|
// halffull list and its full parent may both skip,
|
|
// but are checked to insure matching skip info
|
|
|
|
if (irq->halffull) continue;
|
|
if (irq->copy) continue;
|
|
|
|
// check all other lists
|
|
|
|
for (j = 0; j < nrequest; j++) {
|
|
if (i == j) continue;
|
|
jrq = requests[j];
|
|
|
|
// can only skip from a perpetual non-skip list
|
|
|
|
if (jrq->occasional) continue;
|
|
if (jrq->skip) continue;
|
|
|
|
// both lists must be half, or both full
|
|
|
|
if (irq->half != jrq->half) continue;
|
|
if (irq->full != jrq->full) continue;
|
|
|
|
// both lists must be newton on, or both newton off
|
|
// IJ newton = 1 for newton on, 2 for newton off
|
|
|
|
inewton = irq->newton;
|
|
if (inewton == 0) inewton = force->newton_pair ? 1 : 2;
|
|
jnewton = jrq->newton;
|
|
if (jnewton == 0) jnewton = force->newton_pair ? 1 : 2;
|
|
if (inewton != jnewton) continue;
|
|
|
|
// these flags must be same,
|
|
// else 2 lists do not store same pairs
|
|
// or their data structures are different
|
|
// this includes custom cutoff set by requestor
|
|
// NOTE: need check for 2 Kokkos flags?
|
|
|
|
if (irq->ghost != jrq->ghost) continue;
|
|
if (irq->size != jrq->size) continue;
|
|
if (irq->history != jrq->history) continue;
|
|
if (irq->bond != jrq->bond) continue;
|
|
if (irq->omp != jrq->omp) continue;
|
|
if (irq->intel != jrq->intel) continue;
|
|
if (irq->kokkos_host != jrq->kokkos_host) continue;
|
|
if (irq->kokkos_device != jrq->kokkos_device) continue;
|
|
if (irq->ssa != jrq->ssa) continue;
|
|
if (irq->cut != jrq->cut) continue;
|
|
if (irq->cutoff != jrq->cutoff) continue;
|
|
|
|
// 2 lists are a match
|
|
|
|
break;
|
|
}
|
|
|
|
// if matching list exists, point to it
|
|
// else create a new identical list except non-skip
|
|
// for new list, set neigh = 1, skip = 0, no skip vec/array,
|
|
// copy unique flag (since copy_request() will not do it)
|
|
// note: parents of skip lists do not have associated history
|
|
// b/c child skip lists have the associated history
|
|
|
|
if (j < nrequest) irq->skiplist = j;
|
|
else {
|
|
int newrequest = request(this,-1);
|
|
irq->skiplist = newrequest;
|
|
|
|
nrq = requests[newrequest];
|
|
nrq->copy_request(irq,0);
|
|
nrq->pair = nrq->fix = nrq->compute = nrq->command = 0;
|
|
nrq->neigh = 1;
|
|
nrq->skip = 0;
|
|
if (irq->unique) nrq->unique = 1;
|
|
}
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
scan NeighRequests just added by morph_skip for hybrid granular
|
|
adjust newton/oneside parent settings if children require onesided skipping
|
|
also set children off2on flag if parent becomes a newton off list
|
|
this is needed because line/gran and tri/gran pair styles
|
|
require onesided neigh lists and system newton on,
|
|
but parent list must be newton off to enable the onesided skipping
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::morph_granular()
|
|
{
|
|
int i,j;
|
|
NeighRequest *irq,*jrq;
|
|
|
|
for (i = 0; i < nrequest; i++) {
|
|
irq = requests[i];
|
|
|
|
// only examine NeighRequests added by morph_skip()
|
|
// only those with size attribute for granular systems
|
|
|
|
if (!irq->neigh) continue;
|
|
if (!irq->size) continue;
|
|
|
|
// check children of this list
|
|
|
|
int onesided = -1;
|
|
for (j = 0; j < nrequest; j++) {
|
|
jrq = requests[j];
|
|
|
|
// only consider JRQ pair, size lists that skip from Irq list
|
|
|
|
if (!jrq->pair) continue;
|
|
if (!jrq->size) continue;
|
|
if (!jrq->skip || jrq->skiplist != i) continue;
|
|
|
|
// onesided = -1 if no children
|
|
// onesided = 0/1 = child granonesided value if same for all children
|
|
// onesided = 2 if children have different granonesided values
|
|
|
|
if (onesided < 0) onesided = jrq->granonesided;
|
|
else if (onesided != jrq->granonesided) onesided = 2;
|
|
if (onesided == 2) break;
|
|
}
|
|
|
|
// if onesided = 2, parent has children with both granonesided = 0/1
|
|
// force parent newton off (newton = 2) to enable onesided skip by child
|
|
// set parent granonesided = 0, so it stores all neighs in usual manner
|
|
// set off2on = 1 for all children, since they expect newton on lists
|
|
// this is b/c granonesided only set by line/gran and tri/gran which
|
|
// both require system newton on
|
|
|
|
if (onesided == 2) {
|
|
irq->newton = 2;
|
|
irq->granonesided = 0;
|
|
|
|
for (j = 0; j < nrequest; j++) {
|
|
jrq = requests[j];
|
|
|
|
// only consider JRQ pair, size lists that skip from Irq list
|
|
|
|
if (!jrq->pair) continue;
|
|
if (!jrq->size) continue;
|
|
if (!jrq->skip || jrq->skiplist != i) continue;
|
|
|
|
jrq->off2on = 1;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
scan NeighRequests for possible half lists to derive from full lists
|
|
if 2 requests match, set half list to derive from full list
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::morph_halffull()
|
|
{
|
|
int i,j;
|
|
NeighRequest *irq,*jrq;
|
|
|
|
for (i = 0; i < nrequest; i++) {
|
|
irq = requests[i];
|
|
|
|
// only processing half lists
|
|
|
|
if (!irq->half) continue;
|
|
|
|
// these lists are created other ways, no need for halffull
|
|
// do want to process skip lists
|
|
|
|
if (irq->copy) continue;
|
|
|
|
// check all other lists
|
|
|
|
for (j = 0; j < nrequest; j++) {
|
|
if (i == j) continue;
|
|
jrq = requests[j];
|
|
|
|
// can only derive from a perpetual full list
|
|
// newton setting of derived list does not matter
|
|
|
|
if (jrq->occasional) continue;
|
|
if (!jrq->full) continue;
|
|
|
|
// these flags must be same,
|
|
// else 2 lists do not store same pairs
|
|
// or their data structures are different
|
|
// this includes custom cutoff set by requestor
|
|
|
|
if (irq->ghost != jrq->ghost) continue;
|
|
if (irq->size != jrq->size) continue;
|
|
if (irq->history != jrq->history) continue;
|
|
if (irq->bond != jrq->bond) continue;
|
|
if (irq->omp != jrq->omp) continue;
|
|
if (irq->intel != jrq->intel) continue;
|
|
if (irq->kokkos_host != jrq->kokkos_host) continue;
|
|
if (irq->kokkos_device != jrq->kokkos_device) continue;
|
|
if (irq->ssa != jrq->ssa) continue;
|
|
if (irq->cut != jrq->cut) continue;
|
|
if (irq->cutoff != jrq->cutoff) continue;
|
|
|
|
// skip flag must be same
|
|
// if both are skip lists, skip info must match
|
|
|
|
if (irq->skip != jrq->skip) continue;
|
|
if (irq->skip && irq->same_skip(jrq) == 0) continue;
|
|
|
|
// 2 lists are a match
|
|
|
|
break;
|
|
}
|
|
|
|
// if matching list exists, point to it
|
|
|
|
if (j < nrequest) {
|
|
irq->halffull = 1;
|
|
irq->halffulllist = j;
|
|
}
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
scan NeighRequests for possible copies
|
|
if 2 requests match, turn one into a copy of the other
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::morph_copy()
|
|
{
|
|
int i,j,inewton,jnewton;
|
|
NeighRequest *irq,*jrq;
|
|
|
|
for (i = 0; i < nrequest; i++) {
|
|
irq = requests[i];
|
|
|
|
// this list is already a copy list due to another morph method
|
|
|
|
if (irq->copy) continue;
|
|
|
|
// check all other lists
|
|
|
|
for (j = 0; j < nrequest; j++) {
|
|
if (i == j) continue;
|
|
jrq = requests[j];
|
|
|
|
// other list is already copied from this one
|
|
|
|
if (jrq->copy && jrq->copylist == i) continue;
|
|
|
|
// other list (jrq) to copy from must be perpetual
|
|
// list that becomes a copy list (irq) can be perpetual or occasional
|
|
// if both lists are perpetual, require j < i
|
|
// to prevent circular dependence with 3 or more copies of a list
|
|
|
|
if (jrq->occasional) continue;
|
|
if (!irq->occasional && j > i) continue;
|
|
|
|
// both lists must be half, or both full
|
|
|
|
if (irq->half != jrq->half) continue;
|
|
if (irq->full != jrq->full) continue;
|
|
|
|
// both lists must be newton on, or both newton off
|
|
// IJ newton = 1 for newton on, 2 for newton off
|
|
|
|
inewton = irq->newton;
|
|
if (inewton == 0) inewton = force->newton_pair ? 1 : 2;
|
|
jnewton = jrq->newton;
|
|
if (jnewton == 0) jnewton = force->newton_pair ? 1 : 2;
|
|
if (inewton != jnewton) continue;
|
|
|
|
// ok for non-ghost list to copy from ghost list, but not vice versa
|
|
|
|
if (irq->ghost && !jrq->ghost) continue;
|
|
|
|
// do not copy from a list with respa middle/inner
|
|
// b/c its outer list will not be complete
|
|
|
|
if (jrq->respamiddle) continue;
|
|
if (jrq->respainner) continue;
|
|
|
|
// these flags must be same,
|
|
// else 2 lists do not store same pairs
|
|
// or their data structures are different
|
|
// this includes custom cutoff set by requestor
|
|
// no need to check omp b/c it stores same pairs
|
|
// NOTE: need check for 2 Kokkos flags?
|
|
|
|
if (irq->size != jrq->size) continue;
|
|
if (irq->history != jrq->history) continue;
|
|
if (irq->bond != jrq->bond) continue;
|
|
if (irq->intel != jrq->intel) continue;
|
|
if (irq->kokkos_host && !jrq->kokkos_host) continue;
|
|
if (irq->kokkos_device && !jrq->kokkos_device) continue;
|
|
if (irq->ssa != jrq->ssa) continue;
|
|
if (irq->cut != jrq->cut) continue;
|
|
if (irq->cutoff != jrq->cutoff) continue;
|
|
|
|
// skip flag must be same
|
|
// if both are skip lists, skip info must match
|
|
|
|
if (irq->skip != jrq->skip) continue;
|
|
if (irq->skip && irq->same_skip(jrq) == 0) continue;
|
|
|
|
// 2 lists are a match
|
|
|
|
break;
|
|
}
|
|
|
|
// turn list I into a copy of list J
|
|
// do not copy a list from another copy list, but from its parent list
|
|
|
|
if (j < nrequest) {
|
|
irq->copy = 1;
|
|
if (jrq->copy) irq->copylist = jrq->copylist;
|
|
else irq->copylist = j;
|
|
}
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
create and initialize NTopo classes
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::init_topology()
|
|
{
|
|
int i,m;
|
|
|
|
if (atom->molecular == Atom::ATOMIC) return;
|
|
|
|
// set flags that determine which topology neighbor classes to use
|
|
// these settings could change from run to run, depending on fixes defined
|
|
// bonds,etc can only be broken for atom->molecular = Atom::MOLECULAR, not Atom::TEMPLATE
|
|
// SHAKE sets bonds and angles negative
|
|
// gcmc sets all bonds, angles, etc negative
|
|
// bond_quartic sets bonds to 0
|
|
// delete_bonds sets all interactions negative
|
|
|
|
int bond_off = 0;
|
|
int angle_off = 0;
|
|
for (i = 0; i < modify->nfix; i++)
|
|
if (utils::strmatch(modify->fix[i]->style,"^shake")
|
|
|| utils::strmatch(modify->fix[i]->style,"^rattle"))
|
|
bond_off = angle_off = 1;
|
|
if (force->bond && force->bond_match("quartic")) bond_off = 1;
|
|
|
|
if (atom->avec->bonds_allow && atom->molecular == Atom::MOLECULAR) {
|
|
for (i = 0; i < atom->nlocal; i++) {
|
|
if (bond_off) break;
|
|
for (m = 0; m < atom->num_bond[i]; m++)
|
|
if (atom->bond_type[i][m] <= 0) bond_off = 1;
|
|
}
|
|
}
|
|
|
|
if (atom->avec->angles_allow && atom->molecular == Atom::MOLECULAR) {
|
|
for (i = 0; i < atom->nlocal; i++) {
|
|
if (angle_off) break;
|
|
for (m = 0; m < atom->num_angle[i]; m++)
|
|
if (atom->angle_type[i][m] <= 0) angle_off = 1;
|
|
}
|
|
}
|
|
|
|
int dihedral_off = 0;
|
|
if (atom->avec->dihedrals_allow && atom->molecular == Atom::MOLECULAR) {
|
|
for (i = 0; i < atom->nlocal; i++) {
|
|
if (dihedral_off) break;
|
|
for (m = 0; m < atom->num_dihedral[i]; m++)
|
|
if (atom->dihedral_type[i][m] <= 0) dihedral_off = 1;
|
|
}
|
|
}
|
|
|
|
int improper_off = 0;
|
|
if (atom->avec->impropers_allow && atom->molecular == Atom::MOLECULAR) {
|
|
for (i = 0; i < atom->nlocal; i++) {
|
|
if (improper_off) break;
|
|
for (m = 0; m < atom->num_improper[i]; m++)
|
|
if (atom->improper_type[i][m] <= 0) improper_off = 1;
|
|
}
|
|
}
|
|
|
|
for (i = 0; i < modify->nfix; i++)
|
|
if ((strcmp(modify->fix[i]->style,"gcmc") == 0))
|
|
bond_off = angle_off = dihedral_off = improper_off = 1;
|
|
|
|
// sync on/off settings across all procs
|
|
|
|
int onoff = bond_off;
|
|
MPI_Allreduce(&onoff,&bond_off,1,MPI_INT,MPI_MAX,world);
|
|
onoff = angle_off;
|
|
MPI_Allreduce(&onoff,&angle_off,1,MPI_INT,MPI_MAX,world);
|
|
onoff = dihedral_off;
|
|
MPI_Allreduce(&onoff,&dihedral_off,1,MPI_INT,MPI_MAX,world);
|
|
onoff = improper_off;
|
|
MPI_Allreduce(&onoff,&improper_off,1,MPI_INT,MPI_MAX,world);
|
|
|
|
// instantiate NTopo classes
|
|
|
|
if (atom->avec->bonds_allow) {
|
|
int old_bondwhich = bondwhich;
|
|
if (atom->molecular == Atom::TEMPLATE) bondwhich = TEMPLATE;
|
|
else if (bond_off) bondwhich = PARTIAL;
|
|
else bondwhich = ALL;
|
|
if (!neigh_bond || bondwhich != old_bondwhich) {
|
|
delete neigh_bond;
|
|
if (bondwhich == ALL)
|
|
neigh_bond = new NTopoBondAll(lmp);
|
|
else if (bondwhich == PARTIAL)
|
|
neigh_bond = new NTopoBondPartial(lmp);
|
|
else if (bondwhich == TEMPLATE)
|
|
neigh_bond = new NTopoBondTemplate(lmp);
|
|
}
|
|
}
|
|
|
|
if (atom->avec->angles_allow) {
|
|
int old_anglewhich = anglewhich;
|
|
if (atom->molecular == Atom::TEMPLATE) anglewhich = TEMPLATE;
|
|
else if (angle_off) anglewhich = PARTIAL;
|
|
else anglewhich = ALL;
|
|
if (!neigh_angle || anglewhich != old_anglewhich) {
|
|
delete neigh_angle;
|
|
if (anglewhich == ALL)
|
|
neigh_angle = new NTopoAngleAll(lmp);
|
|
else if (anglewhich == PARTIAL)
|
|
neigh_angle = new NTopoAnglePartial(lmp);
|
|
else if (anglewhich == TEMPLATE)
|
|
neigh_angle = new NTopoAngleTemplate(lmp);
|
|
}
|
|
}
|
|
|
|
if (atom->avec->dihedrals_allow) {
|
|
int old_dihedralwhich = dihedralwhich;
|
|
if (atom->molecular == Atom::TEMPLATE) dihedralwhich = TEMPLATE;
|
|
else if (dihedral_off) dihedralwhich = PARTIAL;
|
|
else dihedralwhich = ALL;
|
|
if (!neigh_dihedral || dihedralwhich != old_dihedralwhich) {
|
|
delete neigh_dihedral;
|
|
if (dihedralwhich == ALL)
|
|
neigh_dihedral = new NTopoDihedralAll(lmp);
|
|
else if (dihedralwhich == PARTIAL)
|
|
neigh_dihedral = new NTopoDihedralPartial(lmp);
|
|
else if (dihedralwhich == TEMPLATE)
|
|
neigh_dihedral = new NTopoDihedralTemplate(lmp);
|
|
}
|
|
}
|
|
|
|
if (atom->avec->impropers_allow) {
|
|
int old_improperwhich = improperwhich;
|
|
if (atom->molecular == Atom::TEMPLATE) improperwhich = TEMPLATE;
|
|
else if (improper_off) improperwhich = PARTIAL;
|
|
else improperwhich = ALL;
|
|
if (!neigh_improper || improperwhich != old_improperwhich) {
|
|
delete neigh_improper;
|
|
if (improperwhich == ALL)
|
|
neigh_improper = new NTopoImproperAll(lmp);
|
|
else if (improperwhich == PARTIAL)
|
|
neigh_improper = new NTopoImproperPartial(lmp);
|
|
else if (improperwhich == TEMPLATE)
|
|
neigh_improper = new NTopoImproperTemplate(lmp);
|
|
}
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
output summary of pairwise neighbor list info
|
|
only called by proc 0
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::print_pairwise_info()
|
|
{
|
|
int i;
|
|
NeighRequest *rq;
|
|
|
|
const double cutghost = MAX(cutneighmax,comm->cutghostuser);
|
|
|
|
double binsize, bbox[3];
|
|
bbox[0] = bboxhi[0]-bboxlo[0];
|
|
bbox[1] = bboxhi[1]-bboxlo[1];
|
|
bbox[2] = bboxhi[2]-bboxlo[2];
|
|
if (binsizeflag) binsize = binsize_user;
|
|
else if (style == Neighbor::BIN) binsize = 0.5*cutneighmax;
|
|
else binsize = 0.5*cutneighmin;
|
|
if (binsize == 0.0) binsize = bbox[0];
|
|
|
|
int nperpetual = 0;
|
|
int noccasional = 0;
|
|
int nextra = 0;
|
|
for (i = 0; i < nlist; i++) {
|
|
if (lists[i]->pair_method == 0) nextra++;
|
|
else if (lists[i]->occasional) noccasional++;
|
|
else nperpetual++;
|
|
}
|
|
|
|
std::string out = "Neighbor list info ...\n";
|
|
out += fmt::format(" update every {} steps, delay {} steps, check {}\n",
|
|
every,delay,dist_check ? "yes" : "no");
|
|
out += fmt::format(" max neighbors/atom: {}, page size: {}\n",
|
|
oneatom, pgsize);
|
|
out += fmt::format(" master list distance cutoff = {:.8g}\n",cutneighmax);
|
|
out += fmt::format(" ghost atom cutoff = {:.8g}\n",cutghost);
|
|
if (style != Neighbor::NSQ)
|
|
out += fmt::format(" binsize = {:.8g}, bins = {:g} {:g} {:g}\n",binsize,
|
|
ceil(bbox[0]/binsize), ceil(bbox[1]/binsize),
|
|
ceil(bbox[2]/binsize));
|
|
|
|
out += fmt::format(" {} neighbor lists, perpetual/occasional/extra = {} {} {}\n",
|
|
nlist,nperpetual,noccasional,nextra);
|
|
|
|
for (i = 0; i < nlist; i++) {
|
|
rq = requests[i];
|
|
if (rq->pair) {
|
|
char *pname = force->pair_match_ptr((Pair *) rq->requestor);
|
|
if (pname) out += fmt::format(" ({}) pair {}",i+1,pname);
|
|
else out += fmt::format(" ({}) pair (none)",i+1);
|
|
} else if (rq->fix) {
|
|
out += fmt::format(" ({}) fix {}",i+1,((Fix *) rq->requestor)->style);
|
|
} else if (rq->compute) {
|
|
out += fmt::format(" ({}) compute {}",i+1,((Compute *) rq->requestor)->style);
|
|
} else if (rq->command) {
|
|
out += fmt::format(" ({}) command {}",i+1,rq->command_style);
|
|
} else if (rq->neigh) {
|
|
out += fmt::format(" ({}) neighbor class addition",i+1);
|
|
}
|
|
|
|
if (rq->occasional) out += ", occasional";
|
|
else out += ", perpetual";
|
|
|
|
// order these to get single output of most relevant
|
|
|
|
if (rq->copy)
|
|
out += fmt::format(", copy from ({})",rq->copylist+1);
|
|
else if (rq->halffull)
|
|
out += fmt::format(", half/full from ({})",rq->halffulllist+1);
|
|
else if (rq->skip)
|
|
out += fmt::format(", skip from ({})",rq->skiplist+1);
|
|
out += "\n";
|
|
|
|
// list of neigh list attributes
|
|
|
|
out += " attributes: ";
|
|
if (rq->half) out += "half";
|
|
else if (rq->full) out += "full";
|
|
|
|
if (rq->newton == 0) {
|
|
if (force->newton_pair) out += ", newton on";
|
|
else out += ", newton off";
|
|
} else if (rq->newton == 1) out += ", newton on";
|
|
else if (rq->newton == 2) out += ", newton off";
|
|
|
|
if (rq->ghost) out += ", ghost";
|
|
if (rq->size) out += ", size";
|
|
if (rq->history) out += ", history";
|
|
if (rq->granonesided) out += ", onesided";
|
|
if (rq->respamiddle) out += ", respa outer/middle/inner";
|
|
else if (rq->respainner) out += ", respa outer/inner";
|
|
if (rq->bond) out += ", bond";
|
|
if (rq->omp) out += ", omp";
|
|
if (rq->intel) out += ", intel";
|
|
if (rq->kokkos_device) out += ", kokkos_device";
|
|
if (rq->kokkos_host) out += ", kokkos_host";
|
|
if (rq->ssa) out += ", ssa";
|
|
if (rq->cut) out += fmt::format(", cut {}",rq->cutoff);
|
|
if (rq->off2on) out += ", off2on";
|
|
out += "\n";
|
|
|
|
out += " ";
|
|
if (lists[i]->pair_method == 0) out += "pair build: none\n";
|
|
else out += fmt::format("pair build: {}\n",pairnames[lists[i]->pair_method-1]);
|
|
|
|
out += " ";
|
|
if (lists[i]->stencil_method == 0) out += "stencil: none\n";
|
|
else out += fmt::format("stencil: {}\n",stencilnames[lists[i]->stencil_method-1]);
|
|
|
|
out += " ";
|
|
if (lists[i]->bin_method == 0) out += "bin: none\n";
|
|
else out += fmt::format("bin: {}\n",binnames[lists[i]->bin_method-1]);
|
|
|
|
}
|
|
utils::logmesg(lmp,out);
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
make copy of current requests and Neighbor params
|
|
used to compare to when next run occurs
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::requests_new2old()
|
|
{
|
|
for (int i = 0; i < old_nrequest; i++) delete old_requests[i];
|
|
memory->sfree(old_requests);
|
|
|
|
old_nrequest = nrequest;
|
|
old_requests = (NeighRequest **)
|
|
memory->smalloc(old_nrequest*sizeof(NeighRequest *),"neighbor:old_requests");
|
|
|
|
for (int i = 0; i < old_nrequest; i++)
|
|
old_requests[i] = new NeighRequest(requests[i]);
|
|
|
|
old_style = style;
|
|
old_triclinic = triclinic;
|
|
old_pgsize = pgsize;
|
|
old_oneatom = oneatom;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
find and return request made by classptr
|
|
if not found or classptr = nullptr, return nullptr
|
|
TODO: should have optional argument "id" to match ID if multiple requests
|
|
------------------------------------------------------------------------- */
|
|
|
|
NeighRequest *Neighbor::find_request(void *classptr)
|
|
{
|
|
if (classptr == nullptr) return nullptr;
|
|
|
|
for (int i = 0; i < nrequest; i++)
|
|
if (requests[i]->requestor == classptr) return requests[i];
|
|
|
|
return nullptr;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
return vector with neighbor list requests from pair styles
|
|
------------------------------------------------------------------------- */
|
|
|
|
const std::vector<NeighRequest *> Neighbor::get_pair_requests() const
|
|
{
|
|
std::vector<NeighRequest *> matches;
|
|
for (int i=0; i < nrequest; ++i)
|
|
if (requests[i]->pair) matches.push_back(requests[i]);
|
|
return matches;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
find and return list requested by classptr
|
|
if not found or classptr = nullptr, return nullptr
|
|
TODO: should have optional argument "id" to match ID if multiple requests
|
|
------------------------------------------------------------------------- */
|
|
|
|
NeighList *Neighbor::find_list(void *classptr)
|
|
{
|
|
if (classptr == nullptr) return nullptr;
|
|
|
|
for (int i = 0; i < nrequest; i++)
|
|
if (lists[i]->requestor == classptr) return lists[i];
|
|
|
|
return nullptr;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
assign NBin class to a NeighList
|
|
use neigh request settings to build mask
|
|
match mask to list of masks of known Nbin classes
|
|
return index+1 of match in list of masks
|
|
return 0 for no binning
|
|
return -1 if no match
|
|
------------------------------------------------------------------------- */
|
|
|
|
int Neighbor::choose_bin(NeighRequest *rq)
|
|
{
|
|
// no binning needed
|
|
|
|
if (style == Neighbor::NSQ) return 0;
|
|
if (rq->skip || rq->copy || rq->halffull) return 0;
|
|
|
|
// use request settings to match exactly one NBin class mask
|
|
// checks are bitwise using NeighConst bit masks
|
|
|
|
int mask;
|
|
|
|
for (int i = 0; i < nbclass; i++) {
|
|
mask = binmasks[i];
|
|
|
|
// require match of these request flags and mask bits
|
|
// (!A != !B) is effectively a logical xor
|
|
|
|
if (!rq->intel != !(mask & NB_INTEL)) continue;
|
|
if (!rq->ssa != !(mask & NB_SSA)) continue;
|
|
if (!rq->kokkos_device != !(mask & NB_KOKKOS_DEVICE)) continue;
|
|
if (!rq->kokkos_host != !(mask & NB_KOKKOS_HOST)) continue;
|
|
|
|
// 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;
|
|
}
|
|
|
|
// error return if matched none
|
|
|
|
return -1;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
assign NStencil class to a NeighList
|
|
use neigh request settings to build mask
|
|
match mask to list of masks of known NStencil classes
|
|
return index+1 of match in list of masks
|
|
return 0 for no binning
|
|
return -1 if no match
|
|
------------------------------------------------------------------------- */
|
|
|
|
int Neighbor::choose_stencil(NeighRequest *rq)
|
|
{
|
|
// no stencil creation needed
|
|
|
|
if (style == Neighbor::NSQ) return 0;
|
|
if (rq->skip || rq->copy || rq->halffull) return 0;
|
|
|
|
// convert newton request to newtflag = on or off
|
|
|
|
int newtflag = 1;
|
|
if (rq->newton == 0 && newton_pair) newtflag = 1;
|
|
else if (rq->newton == 0 && !newton_pair) newtflag = 0;
|
|
else if (rq->newton == 1) newtflag = 1;
|
|
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;
|
|
|
|
//printf("STENCIL RQ FLAGS: hff %d %d n %d g %d s %d newtflag %d fullflag %d\n",
|
|
// rq->half,rq->full,rq->newton,rq->ghost,rq->ssa,
|
|
// newtflag, fullflag);
|
|
|
|
// use request and system settings to match exactly one NStencil class mask
|
|
// checks are bitwise using NeighConst bit masks
|
|
|
|
int mask;
|
|
|
|
for (int i = 0; i < nsclass; i++) {
|
|
mask = stencilmasks[i];
|
|
|
|
//printf("III %d: half %d full %d ghost %d ssa %d\n",
|
|
// i,mask & NS_HALF,mask & NS_FULL,mask & NS_GHOST,mask & NS_SSA);
|
|
|
|
// exactly one of half or full is set and must match
|
|
|
|
if (fullflag) {
|
|
if (!(mask & NS_FULL)) continue;
|
|
} else {
|
|
if (!(mask & NS_HALF)) continue;
|
|
}
|
|
|
|
// require match of these request flags and mask bits
|
|
// (!A != !B) is effectively a logical xor
|
|
|
|
if (!rq->ghost != !(mask & NS_GHOST)) continue;
|
|
if (!rq->ssa != !(mask & NS_SSA)) continue;
|
|
|
|
// neighbor style is one of BIN, MULTI_OLD, or MULTI and must match
|
|
|
|
if (style == Neighbor::BIN) {
|
|
if (!(mask & NS_BIN)) continue;
|
|
} else if (style == Neighbor::MULTI_OLD) {
|
|
if (!(mask & NS_MULTI_OLD)) continue;
|
|
} else if (style == Neighbor::MULTI) {
|
|
if (!(mask & NS_MULTI)) continue;
|
|
}
|
|
|
|
// dimension is 2 or 3 and must match
|
|
|
|
if (dimension == 2) {
|
|
if (!(mask & NS_2D)) continue;
|
|
} else if (dimension == 3) {
|
|
if (!(mask & NS_3D)) continue;
|
|
}
|
|
|
|
// domain triclinic flag is on or off and must match
|
|
|
|
if (triclinic) {
|
|
if (!(mask & NS_TRI)) continue;
|
|
} else if (!triclinic) {
|
|
if (!(mask & NS_ORTHO)) continue;
|
|
}
|
|
|
|
return i+1;
|
|
}
|
|
|
|
// error return if matched none
|
|
|
|
return -1;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
assign NPair class to a NeighList
|
|
use neigh request settings to build mask
|
|
match mask to list of masks of known NPair classes
|
|
return index+1 of match in list of masks
|
|
return 0 for no binning
|
|
return -1 if no match
|
|
------------------------------------------------------------------------- */
|
|
|
|
int Neighbor::choose_pair(NeighRequest *rq)
|
|
{
|
|
// error check for includegroup with ghost neighbor request
|
|
|
|
if (includegroup && rq->ghost)
|
|
error->all(FLERR,"Neighbor include group not allowed with ghost neighbors");
|
|
|
|
// convert newton request to newtflag = on or off
|
|
|
|
bool newtflag;
|
|
if (rq->newton == 0 && newton_pair) newtflag = true;
|
|
else if (rq->newton == 0 && !newton_pair) newtflag = false;
|
|
else if (rq->newton == 1) newtflag = true;
|
|
else if (rq->newton == 2) newtflag = false;
|
|
else error->all(FLERR,"Illegal 'newton' flag in neighbor list request");
|
|
|
|
int molecular = atom->molecular;
|
|
|
|
//printf("PAIR RQ FLAGS: hf %d %d n %d g %d sz %d gos %d r %d b %d o %d i %d "
|
|
// "kk %d %d ss %d dn %d sk %d cp %d hf %d oo %d\n",
|
|
// rq->half,rq->full,rq->newton,rq->ghost,rq->size,
|
|
// rq->granonesided,rq->respaouter,rq->bond,rq->omp,rq->intel,
|
|
// rq->kokkos_host,rq->kokkos_device,rq->ssa,rq->dnum,
|
|
// rq->skip,rq->copy,rq->halffull,rq->off2on);
|
|
|
|
// use request and system settings to match exactly one NPair class mask
|
|
// checks are bitwise using NeighConst bit masks
|
|
|
|
int mask;
|
|
|
|
for (int i = 0; i < npclass; i++) {
|
|
mask = pairmasks[i];
|
|
|
|
//printf(" PAIR NAMES i %d %d name %s mask %d\n",i,nrequest,
|
|
// pairnames[i],pairmasks[i]);
|
|
|
|
// if copy request, no further checks needed, just return or continue
|
|
// Kokkos device/host flags must also match in order to copy
|
|
|
|
if (rq->copy) {
|
|
if (!(mask & NP_COPY)) continue;
|
|
if (rq->kokkos_device || rq->kokkos_host) {
|
|
if (!rq->kokkos_device != !(mask & NP_KOKKOS_DEVICE)) continue;
|
|
if (!rq->kokkos_host != !(mask & NP_KOKKOS_HOST)) continue;
|
|
}
|
|
if (!requests[rq->copylist]->kokkos_device != !(mask & NP_KOKKOS_DEVICE)) continue;
|
|
if (!requests[rq->copylist]->kokkos_host != !(mask & NP_KOKKOS_HOST)) continue;
|
|
return i+1;
|
|
}
|
|
|
|
// exactly one of half or full is set and must match
|
|
|
|
if (rq->half) {
|
|
if (!(mask & NP_HALF)) continue;
|
|
} else if (rq->full) {
|
|
if (!(mask & NP_FULL)) continue;
|
|
}
|
|
|
|
// newtflag is on or off and must match
|
|
|
|
if (newtflag) {
|
|
if (!(mask & NP_NEWTON)) continue;
|
|
} else if (!newtflag) {
|
|
if (!(mask & NP_NEWTOFF)) continue;
|
|
}
|
|
|
|
// if molecular on, do not match ATOMONLY (b/c a MOLONLY Npair exists)
|
|
// if molecular off, do not match MOLONLY (b/c an ATOMONLY Npair exists)
|
|
|
|
if (molecular != Atom::ATOMIC) {
|
|
if (mask & NP_ATOMONLY) continue;
|
|
} else if (molecular == Atom::ATOMIC) {
|
|
if (mask & NP_MOLONLY) continue;
|
|
}
|
|
|
|
// require match of these request flags and mask bits
|
|
// (!A != !B) is effectively a logical xor
|
|
|
|
if (!rq->ghost != !(mask & NP_GHOST)) continue;
|
|
if (!rq->size != !(mask & NP_SIZE)) continue;
|
|
if (!rq->respaouter != !(mask & NP_RESPA)) continue;
|
|
if (!rq->granonesided != !(mask & NP_ONESIDE)) continue;
|
|
if (!rq->respaouter != !(mask & NP_RESPA)) continue;
|
|
if (!rq->bond != !(mask & NP_BOND)) continue;
|
|
if (!rq->omp != !(mask & NP_OMP)) continue;
|
|
if (!rq->intel != !(mask & NP_INTEL)) continue;
|
|
if (!rq->kokkos_device != !(mask & NP_KOKKOS_DEVICE)) continue;
|
|
if (!rq->kokkos_host != !(mask & NP_KOKKOS_HOST)) continue;
|
|
if (!rq->ssa != !(mask & NP_SSA)) continue;
|
|
|
|
if (!rq->skip != !(mask & NP_SKIP)) continue;
|
|
|
|
if (!rq->halffull != !(mask & NP_HALF_FULL)) continue;
|
|
if (!rq->off2on != !(mask & NP_OFF2ON)) continue;
|
|
|
|
// neighbor style is one of NSQ, BIN, MULTI_OLD, or MULTI and must match
|
|
|
|
if (style == Neighbor::NSQ) {
|
|
if (!(mask & NP_NSQ)) continue;
|
|
} else if (style == Neighbor::BIN) {
|
|
if (!(mask & NP_BIN)) continue;
|
|
} else if (style == Neighbor::MULTI_OLD) {
|
|
if (!(mask & NP_MULTI_OLD)) continue;
|
|
} else if (style == Neighbor::MULTI) {
|
|
if (!(mask & NP_MULTI)) continue;
|
|
}
|
|
|
|
// domain triclinic flag is on or off and must match
|
|
|
|
if (triclinic) {
|
|
if (!(mask & NP_TRI)) continue;
|
|
} else if (!triclinic) {
|
|
if (!(mask & NP_ORTHO)) continue;
|
|
}
|
|
|
|
return i+1;
|
|
}
|
|
|
|
// error return if matched none
|
|
|
|
return -1;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
called by other classes to request a pairwise neighbor list
|
|
------------------------------------------------------------------------- */
|
|
|
|
int Neighbor::request(void *requestor, int instance)
|
|
{
|
|
if (nrequest == maxrequest) {
|
|
maxrequest += RQDELTA;
|
|
requests = (NeighRequest **)
|
|
memory->srealloc(requests,maxrequest*sizeof(NeighRequest *), "neighbor:requests");
|
|
}
|
|
|
|
requests[nrequest] = new NeighRequest(lmp, nrequest, requestor, instance);
|
|
nrequest++;
|
|
return nrequest-1;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
called by other classes to request a pairwise neighbor list
|
|
------------------------------------------------------------------------- */
|
|
|
|
NeighRequest *Neighbor::add_request(Pair *requestor, int flags)
|
|
{
|
|
int irequest = request(requestor, requestor->instance_me);
|
|
auto req = requests[irequest];
|
|
req->apply_flags(flags);
|
|
return req;
|
|
}
|
|
|
|
NeighRequest *Neighbor::add_request(Fix *requestor, int flags)
|
|
{
|
|
int irequest = request(requestor, requestor->instance_me);
|
|
auto req = requests[irequest];
|
|
req->pair = 0;
|
|
req->fix = 1;
|
|
req->apply_flags(flags);
|
|
return req;
|
|
}
|
|
|
|
NeighRequest *Neighbor::add_request(Compute *requestor, int flags)
|
|
{
|
|
int irequest = request(requestor, requestor->instance_me);
|
|
auto req = requests[irequest];
|
|
req->pair = 0;
|
|
req->compute = 1;
|
|
req->apply_flags(flags);
|
|
return req;
|
|
}
|
|
|
|
NeighRequest *Neighbor::add_request(Command *requestor, const char *style, int flags)
|
|
{
|
|
int irequest = request(requestor, 0);
|
|
auto req = requests[irequest];
|
|
req->pair = 0;
|
|
req->command = 1;
|
|
req->occasional = 1;
|
|
req->command_style = style;
|
|
req->apply_flags(flags);
|
|
return req;
|
|
}
|
|
|
|
// set neighbor list request OpenMP flag
|
|
|
|
void Neighbor::set_omp_neighbor(int flag)
|
|
{
|
|
// flag *all* neighbor list requests as OPENMP threaded,
|
|
// but skip lists already flagged as INTEL threaded
|
|
for (int i = 0; i < nrequest; ++i)
|
|
if (!requests[i]->intel) requests[i]->omp = flag;
|
|
}
|
|
|
|
/* report if there is a neighbor list with the intel flag set */
|
|
|
|
bool Neighbor::has_intel_request() const
|
|
{
|
|
return (((nrequest > 0) && (requests[0]->intel > 0))
|
|
|| ((old_nrequest > 0) && (old_requests[0]->intel > 0)));
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
setup neighbor binning and neighbor stencils
|
|
called before run and every reneighbor if box size/shape changes
|
|
only operates on perpetual lists
|
|
build_one() operates on occasional lists
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::setup_bins()
|
|
{
|
|
// invoke setup_bins() for all NBin
|
|
// actual binning is performed in build()
|
|
|
|
for (int i = 0; i < nbin; i++)
|
|
neigh_bin[i]->setup_bins(style);
|
|
|
|
// invoke create_setup() and create() for all perpetual NStencil
|
|
// same ops performed for occasional lists in build_one()
|
|
|
|
for (int i = 0; i < nstencil_perpetual; i++) {
|
|
neigh_stencil[slist[i]]->create_setup();
|
|
neigh_stencil[slist[i]]->create();
|
|
}
|
|
|
|
last_setup_bins = update->ntimestep;
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
int Neighbor::decide()
|
|
{
|
|
if (must_check) {
|
|
bigint n = update->ntimestep;
|
|
if (restart_check && n == output->next_restart) return 1;
|
|
for (int i = 0; i < fix_check; i++)
|
|
if (n == modify->fix[fixchecklist[i]]->next_reneighbor) return 1;
|
|
}
|
|
|
|
ago++;
|
|
if (ago >= delay && ago % every == 0) {
|
|
if (build_once) return 0;
|
|
if (dist_check == 0) return 1;
|
|
return check_distance();
|
|
} else return 0;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
if any atom moved trigger distance (half of neighbor skin) return 1
|
|
shrink trigger distance if box size has changed
|
|
conservative shrink procedure:
|
|
compute distance each of 8 corners of box has moved since last reneighbor
|
|
reduce skin distance by sum of 2 largest of the 8 values
|
|
if reduced skin distance is negative, set to zero
|
|
new trigger = 1/2 of reduced skin distance
|
|
for orthogonal box, only need 2 lo/hi corners
|
|
for triclinic, need all 8 corners since deformations can displace all 8
|
|
------------------------------------------------------------------------- */
|
|
|
|
int Neighbor::check_distance()
|
|
{
|
|
double delx,dely,delz,rsq;
|
|
double delta,deltasq,delta1,delta2;
|
|
|
|
if (boxcheck) {
|
|
if (triclinic == 0) {
|
|
delx = bboxlo[0] - boxlo_hold[0];
|
|
dely = bboxlo[1] - boxlo_hold[1];
|
|
delz = bboxlo[2] - boxlo_hold[2];
|
|
delta1 = sqrt(delx*delx + dely*dely + delz*delz);
|
|
delx = bboxhi[0] - boxhi_hold[0];
|
|
dely = bboxhi[1] - boxhi_hold[1];
|
|
delz = bboxhi[2] - boxhi_hold[2];
|
|
delta2 = sqrt(delx*delx + dely*dely + delz*delz);
|
|
delta = 0.5 * (skin - (delta1+delta2));
|
|
if (delta < 0.0) delta = 0.0;
|
|
deltasq = delta*delta;
|
|
} else {
|
|
domain->box_corners();
|
|
delta1 = delta2 = 0.0;
|
|
for (int i = 0; i < 8; i++) {
|
|
delx = corners[i][0] - corners_hold[i][0];
|
|
dely = corners[i][1] - corners_hold[i][1];
|
|
delz = corners[i][2] - corners_hold[i][2];
|
|
delta = sqrt(delx*delx + dely*dely + delz*delz);
|
|
if (delta > delta1) delta1 = delta;
|
|
else if (delta > delta2) delta2 = delta;
|
|
}
|
|
delta = 0.5 * (skin - (delta1+delta2));
|
|
if (delta < 0.0) delta = 0.0;
|
|
deltasq = delta*delta;
|
|
}
|
|
} else deltasq = triggersq;
|
|
|
|
double **x = atom->x;
|
|
int nlocal = atom->nlocal;
|
|
if (includegroup) nlocal = atom->nfirst;
|
|
|
|
int flag = 0;
|
|
for (int i = 0; i < nlocal; i++) {
|
|
delx = x[i][0] - xhold[i][0];
|
|
dely = x[i][1] - xhold[i][1];
|
|
delz = x[i][2] - xhold[i][2];
|
|
rsq = delx*delx + dely*dely + delz*delz;
|
|
if (rsq > deltasq) flag = 1;
|
|
}
|
|
|
|
int flagall;
|
|
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
|
|
if (flagall && ago == MAX(every,delay)) ndanger++;
|
|
return flagall;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
build perpetual neighbor lists
|
|
called at setup and every few timesteps during run or minimization
|
|
topology lists also built if topoflag = 1 (Kokkos calls with topoflag=0)
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::build(int topoflag)
|
|
{
|
|
int i,m;
|
|
|
|
ago = 0;
|
|
ncalls++;
|
|
lastcall = update->ntimestep;
|
|
|
|
int nlocal = atom->nlocal;
|
|
int nall = nlocal + atom->nghost;
|
|
// rebuild collection array from scratch
|
|
if (style == Neighbor::MULTI) build_collection(0);
|
|
|
|
// check that using special bond flags will not overflow neigh lists
|
|
|
|
if (nall > NEIGHMASK)
|
|
error->one(FLERR,"Too many local+ghost atoms for neighbor list");
|
|
|
|
// store current atom positions and box size if needed
|
|
|
|
if (dist_check) {
|
|
double **x = atom->x;
|
|
if (includegroup) nlocal = atom->nfirst;
|
|
if (atom->nmax > maxhold) {
|
|
maxhold = atom->nmax;
|
|
memory->destroy(xhold);
|
|
memory->create(xhold,maxhold,3,"neigh:xhold");
|
|
}
|
|
for (i = 0; i < nlocal; i++) {
|
|
xhold[i][0] = x[i][0];
|
|
xhold[i][1] = x[i][1];
|
|
xhold[i][2] = x[i][2];
|
|
}
|
|
if (boxcheck) {
|
|
if (triclinic == 0) {
|
|
boxlo_hold[0] = bboxlo[0];
|
|
boxlo_hold[1] = bboxlo[1];
|
|
boxlo_hold[2] = bboxlo[2];
|
|
boxhi_hold[0] = bboxhi[0];
|
|
boxhi_hold[1] = bboxhi[1];
|
|
boxhi_hold[2] = bboxhi[2];
|
|
} else {
|
|
domain->box_corners();
|
|
corners = domain->corners;
|
|
for (i = 0; i < 8; i++) {
|
|
corners_hold[i][0] = corners[i][0];
|
|
corners_hold[i][1] = corners[i][1];
|
|
corners_hold[i][2] = corners[i][2];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// bin atoms for all NBin instances
|
|
// not just NBin associated with perpetual lists, also occasional lists
|
|
// b/c cannot wait to bin occasional lists in build_one() call
|
|
// if bin then, atoms may have moved outside of proc domain & bin extent,
|
|
// leading to errors or even a crash
|
|
|
|
if (style != Neighbor::NSQ) {
|
|
if (last_setup_bins < 0) setup_bins();
|
|
for (i = 0; i < nbin; i++) {
|
|
neigh_bin[i]->bin_atoms_setup(nall);
|
|
neigh_bin[i]->bin_atoms();
|
|
}
|
|
}
|
|
|
|
// build pairwise lists for all perpetual NPair/NeighList
|
|
// grow() with nlocal/nall args so that only realloc if have to
|
|
|
|
for (i = 0; i < npair_perpetual; i++) {
|
|
m = plist[i];
|
|
if (!lists[m]->copy || lists[m]->kk2cpu)
|
|
lists[m]->grow(nlocal,nall);
|
|
neigh_pair[m]->build_setup();
|
|
neigh_pair[m]->build(lists[m]);
|
|
}
|
|
|
|
// build topology lists for bonds/angles/etc
|
|
|
|
if ((atom->molecular != Atom::ATOMIC) && topoflag) build_topology();
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
build topology neighbor lists: bond, angle, dihedral, improper
|
|
copy their list info back to Neighbor for access by bond/angle/etc classes
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::build_topology()
|
|
{
|
|
if (force->bond) {
|
|
neigh_bond->build();
|
|
nbondlist = neigh_bond->nbondlist;
|
|
bondlist = neigh_bond->bondlist;
|
|
}
|
|
if (force->angle) {
|
|
neigh_angle->build();
|
|
nanglelist = neigh_angle->nanglelist;
|
|
anglelist = neigh_angle->anglelist;
|
|
}
|
|
if (force->dihedral) {
|
|
neigh_dihedral->build();
|
|
ndihedrallist = neigh_dihedral->ndihedrallist;
|
|
dihedrallist = neigh_dihedral->dihedrallist;
|
|
}
|
|
if (force->improper) {
|
|
neigh_improper->build();
|
|
nimproperlist = neigh_improper->nimproperlist;
|
|
improperlist = neigh_improper->improperlist;
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
build a single occasional pairwise neighbor list indexed by I
|
|
called by other classes
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::build_one(class NeighList *mylist, int preflag)
|
|
{
|
|
// check if list structure is initialized
|
|
|
|
if (mylist == nullptr)
|
|
error->all(FLERR,"Trying to build an occasional neighbor list before initialization complete");
|
|
|
|
// build_one() should never be invoked on a perpetual list
|
|
|
|
if (!mylist->occasional) error->all(FLERR,"Neighbor::build_one() invoked on perpetual list");
|
|
|
|
// no need to build if already built since last re-neighbor
|
|
// preflag is set by fix bond/create and fix bond/swap
|
|
// b/c they invoke build_one() on same step neigh list is re-built,
|
|
// but before re-build, so need to use ">" instead of ">="
|
|
|
|
NPair *np = neigh_pair[mylist->index];
|
|
|
|
if (preflag) {
|
|
if (np->last_build > lastcall) return;
|
|
} else {
|
|
if (np->last_build >= lastcall) return;
|
|
}
|
|
|
|
// if this is copy list and parent is occasional list,
|
|
// or this is halffull and parent is occasional list,
|
|
// or this is skip list and parent is occasional list,
|
|
// insure parent is current
|
|
|
|
if (mylist->listcopy && mylist->listcopy->occasional)
|
|
build_one(mylist->listcopy,preflag);
|
|
if (mylist->listfull && mylist->listfull->occasional)
|
|
build_one(mylist->listfull,preflag);
|
|
if (mylist->listskip && mylist->listskip->occasional)
|
|
build_one(mylist->listskip,preflag);
|
|
|
|
// create stencil if hasn't been created since last setup_bins() call
|
|
|
|
NStencil *ns = np->ns;
|
|
if (ns && ns->last_stencil < last_setup_bins) {
|
|
ns->create_setup();
|
|
ns->create();
|
|
}
|
|
|
|
// build the list
|
|
|
|
if (!mylist->copy || mylist->kk2cpu)
|
|
mylist->grow(atom->nlocal,atom->nlocal+atom->nghost);
|
|
np->build_setup();
|
|
np->build(mylist);
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
set neighbor style and skin distance
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::set(int narg, char **arg)
|
|
{
|
|
if (narg != 2) error->all(FLERR,"Illegal neighbor command");
|
|
|
|
skin = utils::numeric(FLERR,arg[0],false,lmp);
|
|
if (skin < 0.0) error->all(FLERR,"Illegal neighbor command");
|
|
|
|
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;
|
|
ncollections = atom->ntypes;
|
|
} else if (strcmp(arg[1],"multi/old") == 0) style = Neighbor::MULTI_OLD;
|
|
else error->all(FLERR,"Illegal neighbor command");
|
|
|
|
if (style == Neighbor::MULTI_OLD && lmp->citeme) lmp->citeme->add(cite_neigh_multi_old);
|
|
if (style == Neighbor::MULTI && lmp->citeme) lmp->citeme->add(cite_neigh_multi);
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
reset timestamps in all NeignBin, NStencil, NPair classes
|
|
so that neighbor lists will rebuild properly with timestep change
|
|
ditto for lastcall and last_setup_bins
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::reset_timestep(bigint /*ntimestep*/)
|
|
{
|
|
for (int i = 0; i < nbin; i++)
|
|
neigh_bin[i]->last_bin = -1;
|
|
for (int i = 0; i < nstencil; i++)
|
|
neigh_stencil[i]->last_stencil = -1;
|
|
for (int i = 0; i < nlist; i++) {
|
|
if (!neigh_pair[i]) continue;
|
|
neigh_pair[i]->last_build = -1;
|
|
}
|
|
|
|
lastcall = -1;
|
|
last_setup_bins = -1;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
modify parameters of the pair-wise neighbor build
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::modify_params(int narg, char **arg)
|
|
{
|
|
int iarg = 0;
|
|
while (iarg < narg) {
|
|
if (strcmp(arg[iarg],"every") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
|
|
every = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
|
|
if (every <= 0) error->all(FLERR,"Illegal neigh_modify command");
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"delay") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
|
|
delay = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
|
|
if (delay < 0) error->all(FLERR,"Illegal neigh_modify command");
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"check") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
|
|
dist_check = utils::logical(FLERR,arg[iarg+1],false,lmp);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"once") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
|
|
build_once = utils::logical(FLERR,arg[iarg+1],false,lmp);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"page") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
|
|
old_pgsize = pgsize;
|
|
pgsize = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"one") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
|
|
old_oneatom = oneatom;
|
|
oneatom = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"binsize") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
|
|
binsize_user = utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
|
if (binsize_user <= 0.0) binsizeflag = 0;
|
|
else binsizeflag = 1;
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"cluster") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
|
|
cluster_check = utils::logical(FLERR,arg[iarg+1],false,lmp);
|
|
iarg += 2;
|
|
|
|
} else if (strcmp(arg[iarg],"include") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
|
|
includegroup = group->find(arg[iarg+1]);
|
|
if (includegroup < 0)
|
|
error->all(FLERR,"Invalid group ID in neigh_modify command");
|
|
if (includegroup && (atom->firstgroupname == nullptr ||
|
|
strcmp(arg[iarg+1],atom->firstgroupname) != 0))
|
|
error->all(FLERR,
|
|
"Neigh_modify include group != atom_modify first group");
|
|
iarg += 2;
|
|
|
|
} else if (strcmp(arg[iarg],"exclude") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
|
|
|
|
if (strcmp(arg[iarg+1],"type") == 0) {
|
|
if (iarg+4 > narg) error->all(FLERR,"Illegal neigh_modify command");
|
|
if (nex_type == maxex_type) {
|
|
maxex_type += EXDELTA;
|
|
memory->grow(ex1_type,maxex_type,"neigh:ex1_type");
|
|
memory->grow(ex2_type,maxex_type,"neigh:ex2_type");
|
|
}
|
|
ex1_type[nex_type] = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
|
|
ex2_type[nex_type] = utils::inumeric(FLERR,arg[iarg+3],false,lmp);
|
|
nex_type++;
|
|
iarg += 4;
|
|
|
|
} else if (strcmp(arg[iarg+1],"group") == 0) {
|
|
if (iarg+4 > narg) error->all(FLERR,"Illegal neigh_modify command");
|
|
if (nex_group == maxex_group) {
|
|
maxex_group += EXDELTA;
|
|
memory->grow(ex1_group,maxex_group,"neigh:ex1_group");
|
|
memory->grow(ex2_group,maxex_group,"neigh:ex2_group");
|
|
}
|
|
ex1_group[nex_group] = group->find(arg[iarg+2]);
|
|
ex2_group[nex_group] = group->find(arg[iarg+3]);
|
|
if (ex1_group[nex_group] == -1 || ex2_group[nex_group] == -1)
|
|
error->all(FLERR,"Invalid group ID in neigh_modify command");
|
|
nex_group++;
|
|
iarg += 4;
|
|
|
|
} else if (strcmp(arg[iarg+1],"molecule/inter") == 0 ||
|
|
strcmp(arg[iarg+1],"molecule/intra") == 0) {
|
|
if (iarg+3 > narg) error->all(FLERR,"Illegal neigh_modify command");
|
|
if (atom->molecule_flag == 0)
|
|
error->all(FLERR,"Neigh_modify exclude molecule "
|
|
"requires atom attribute molecule");
|
|
if (nex_mol == maxex_mol) {
|
|
maxex_mol += EXDELTA;
|
|
memory->grow(ex_mol_group,maxex_mol,"neigh:ex_mol_group");
|
|
if (lmp->kokkos)
|
|
grow_ex_mol_intra_kokkos();
|
|
else
|
|
memory->grow(ex_mol_intra,maxex_mol,"neigh:ex_mol_intra");
|
|
}
|
|
ex_mol_group[nex_mol] = group->find(arg[iarg+2]);
|
|
if (ex_mol_group[nex_mol] == -1)
|
|
error->all(FLERR,"Invalid group ID in neigh_modify command");
|
|
if (strcmp(arg[iarg+1],"molecule/intra") == 0)
|
|
ex_mol_intra[nex_mol] = 1;
|
|
else
|
|
ex_mol_intra[nex_mol] = 0;
|
|
nex_mol++;
|
|
iarg += 3;
|
|
|
|
} else if (strcmp(arg[iarg+1],"none") == 0) {
|
|
nex_type = nex_group = nex_mol = 0;
|
|
iarg += 2;
|
|
|
|
} else error->all(FLERR,"Illegal neigh_modify command");
|
|
} else if (strcmp(arg[iarg],"collection/interval") == 0) {
|
|
if (style != Neighbor::MULTI)
|
|
error->all(FLERR,"Cannot use collection/interval command without multi setting");
|
|
|
|
if (iarg+2 > narg)
|
|
error->all(FLERR,"Invalid collection/interval command");
|
|
ncollections = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
|
|
if (ncollections < 1)
|
|
error->all(FLERR,"Invalid collection/interval command");
|
|
if (iarg+1+ncollections > narg)
|
|
error->all(FLERR,"Invalid collection/interval command");
|
|
|
|
int i;
|
|
|
|
// Invalidate old user cutoffs
|
|
|
|
comm->ncollections_cutoff = 0;
|
|
interval_collection_flag = 1;
|
|
custom_collection_flag = 1;
|
|
memory->grow(collection2cut,ncollections,"neigh:collection2cut");
|
|
|
|
// Set upper cutoff for each collection
|
|
|
|
double cut_interval;
|
|
for (i = 0; i < ncollections; i++){
|
|
cut_interval = utils::numeric(FLERR,arg[iarg+2+i],false,lmp);
|
|
collection2cut[i] = cut_interval;
|
|
|
|
if (i != 0)
|
|
if (collection2cut[i-1] >= collection2cut[i])
|
|
error->all(FLERR,"Nonsequential interval cutoffs in collection/interval setting");
|
|
}
|
|
|
|
iarg += 2 + ncollections;
|
|
} else if (strcmp(arg[iarg],"collection/type") == 0) {
|
|
if (style != Neighbor::MULTI)
|
|
error->all(FLERR,"Cannot use collection/type command without multi setting");
|
|
|
|
if (iarg+2 > narg)
|
|
error->all(FLERR,"Invalid collection/type command");
|
|
ncollections = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
|
|
if (ncollections < 1)
|
|
error->all(FLERR,"Invalid collection/interval command");
|
|
if (iarg+1+ncollections > narg)
|
|
error->all(FLERR,"Invalid collection/type command");
|
|
|
|
int ntypes = atom->ntypes;
|
|
int nlo, nhi, i, k;
|
|
|
|
// Invalidate old user cutoffs
|
|
|
|
comm->ncollections_cutoff = 0;
|
|
interval_collection_flag = 0;
|
|
custom_collection_flag = 1;
|
|
if (!type2collection)
|
|
memory->create(type2collection,ntypes+1,"neigh:type2collection");
|
|
|
|
// Erase previous mapping
|
|
|
|
for (i = 1; i <= ntypes; i++)
|
|
type2collection[i] = -1;
|
|
|
|
// For each custom range, define mapping for types in interval
|
|
|
|
for (i = 0; i < ncollections; i++){
|
|
std::vector<std::string> words = Tokenizer(arg[iarg+2+i], ",").as_vector();
|
|
for (const auto &word : words) {
|
|
utils::bounds(FLERR,word,1,ntypes,nlo,nhi,error);
|
|
for (k = nlo; k <= nhi; k++) {
|
|
if (type2collection[k] != -1)
|
|
error->all(FLERR,"Type specified more than once in collection/type commnd");
|
|
type2collection[k] = i;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Check for undefined atom type
|
|
|
|
for (i = 1; i <= ntypes; i++){
|
|
if (type2collection[i] == -1) {
|
|
error->all(FLERR,"Type missing in collection/type commnd");
|
|
}
|
|
}
|
|
|
|
iarg += 2 + ncollections;
|
|
} else error->all(FLERR,"Illegal neigh_modify command");
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
convenience function to allow modifying parameters from a single string
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::modify_params(const std::string &modcmd)
|
|
{
|
|
auto args = utils::split_words(modcmd);
|
|
char **newarg = new char*[args.size()];
|
|
int i=0;
|
|
for (const auto &arg : args) {
|
|
newarg[i++] = (char *)arg.c_str();
|
|
}
|
|
modify_params(args.size(),newarg);
|
|
delete[] newarg;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
remove the first group-group exclusion matching group1, group2
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::exclusion_group_group_delete(int group1, int group2)
|
|
{
|
|
int m, mlast;
|
|
for (m = 0; m < nex_group; m++)
|
|
if (ex1_group[m] == group1 && ex2_group[m] == group2 )
|
|
break;
|
|
|
|
mlast = m;
|
|
if (mlast == nex_group)
|
|
error->all(FLERR,"Unable to find group-group exclusion");
|
|
|
|
for (m = mlast+1; m < nex_group; m++) {
|
|
ex1_group[m-1] = ex1_group[m];
|
|
ex2_group[m-1] = ex2_group[m];
|
|
ex1_bit[m-1] = ex1_bit[m];
|
|
ex2_bit[m-1] = ex2_bit[m];
|
|
}
|
|
nex_group--;
|
|
}
|
|
|
|
|
|
/* ----------------------------------------------------------------------
|
|
return the value of exclude - used to check compatibility with GPU
|
|
------------------------------------------------------------------------- */
|
|
|
|
int Neighbor::exclude_setting()
|
|
{
|
|
return exclude;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
check if any of the old requested neighbor lists are full
|
|
------------------------------------------------------------------------- */
|
|
|
|
int Neighbor::any_full()
|
|
{
|
|
int any_full = 0;
|
|
for (int i = 0; i < old_nrequest; i++) {
|
|
if (old_requests[i]->full) any_full = 1;
|
|
}
|
|
return any_full;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
populate collection array for multi starting at the index istart
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Neighbor::build_collection(int istart)
|
|
{
|
|
if (style != Neighbor::MULTI)
|
|
error->all(FLERR, "Cannot define atom collections without neighbor style multi");
|
|
|
|
int nmax = atom->nlocal+atom->nghost;
|
|
if (nmax > nmax_collection) {
|
|
nmax_collection = nmax+DELTA_PERATOM;
|
|
memory->grow(collection, nmax_collection, "neigh:collection");
|
|
}
|
|
|
|
if (finite_cut_flag) {
|
|
double cut;
|
|
int icollection;
|
|
for (int i = istart; i < nmax; i++){
|
|
cut = force->pair->atom2cut(i);
|
|
collection[i] = -1;
|
|
|
|
for (icollection = 0; icollection < ncollections; icollection++){
|
|
if (collection2cut[icollection] >= cut) {
|
|
collection[i] = icollection;
|
|
break;
|
|
}
|
|
}
|
|
|
|
if (collection[i] == -1)
|
|
error->one(FLERR, "Atom cutoff exceeds interval cutoffs for multi");
|
|
}
|
|
} else {
|
|
int *type = atom->type;
|
|
for (int i = istart; i < nmax; i++){
|
|
collection[i] = type2collection[type[i]];
|
|
}
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
return # of bytes of allocated memory
|
|
------------------------------------------------------------------------- */
|
|
|
|
double Neighbor::memory_usage()
|
|
{
|
|
double bytes = 0;
|
|
bytes += memory->usage(xhold,maxhold,3);
|
|
|
|
for (int i = 0; i < nlist; i++)
|
|
if (lists[i]) bytes += lists[i]->memory_usage();
|
|
for (int i = 0; i < nstencil; i++)
|
|
bytes += neigh_stencil[i]->memory_usage();
|
|
for (int i = 0; i < nbin; i++)
|
|
bytes += neigh_bin[i]->memory_usage();
|
|
|
|
if (neigh_bond) bytes += neigh_bond->memory_usage();
|
|
if (neigh_angle) bytes += neigh_angle->memory_usage();
|
|
if (neigh_dihedral) bytes += neigh_dihedral->memory_usage();
|
|
if (neigh_improper) bytes += neigh_improper->memory_usage();
|
|
|
|
return bytes;
|
|
}
|