Homec4science
Diffusion lammps (efficient_neuronet) Browse

lammps/src/neighbor.cpp
2e7404c96c42efficient_neuronet

Last Change
Raw File

neighbor.cpp
View Options

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author (triclinic and multi-neigh) : Pieter in 't Veld (SNL)
------------------------------------------------------------------------- */
#include <mpi.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "atom.h"
#include "atom_vec.h"
#include "comm.h"
#include "force.h"
#include "pair.h"
#include "domain.h"
#include "group.h"
#include "modify.h"
#include "fix.h"
#include "compute.h"
#include "update.h"
#include "respa.h"
#include "output.h"
#include "citeme.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define RQDELTA 1
#define EXDELTA 1
#define LB_FACTOR 1.5
#define SMALL 1.0e-6
#define BIG 1.0e20
#define CUT2BIN_RATIO 100
enum{NSQ,BIN,MULTI}; // also in neigh_list.cpp
static const char cite_neigh_multi[] =
"neighbor multi command:\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";
//#define NEIGH_LIST_DEBUG 1
/* ---------------------------------------------------------------------- */
Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp)
{
MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs);
style = BIN;
every = 1;
delay = 10;
dist_check = 1;
pgsize = 100000;
oneatom = 2000;
binsizeflag = 0;
build_once = 0;
cluster_check = 0;
binatomflag = 1;
ago = -1;
cutneighmax = 0.0;
cutneighsq = NULL;
cutneighghostsq = NULL;
cuttype = NULL;
cuttypesq = NULL;
fixchecklist = NULL;
// coords at last neighboring
maxhold = 0;
xhold = NULL;
lastcall = -1;
// binning
maxhead = 0;
binhead = NULL;
maxbin = 0;
bins = NULL;
// pair exclusion list info
includegroup = 0;
nex_type = maxex_type = 0;
ex1_type = ex2_type = NULL;
ex_type = NULL;
nex_group = maxex_group = 0;
ex1_group = ex2_group = ex1_bit = ex2_bit = NULL;
nex_mol = maxex_mol = 0;
ex_mol_group = ex_mol_bit = NULL;
// pair lists
maxatom = 0;
nblist = nglist = nslist = 0;
nlist = 0;
lists = NULL;
pair_build = NULL;
stencil_create = NULL;
blist = glist = slist = NULL;
anyghostlist = 0;
nrequest = maxrequest = 0;
requests = NULL;
old_nrequest = 0;
old_requests = NULL;
old_style = style;
old_triclinic = 0;
old_pgsize = pgsize;
old_oneatom = oneatom;
old_every = every;
old_delay = delay;
old_check = dist_check;
old_cutoff = cutneighmax;
zeroes = NULL;
// bond lists
maxbond = 0;
bondlist = NULL;
maxangle = 0;
anglelist = NULL;
maxdihedral = 0;
dihedrallist = NULL;
maximproper = 0;
improperlist = NULL;
copymode = 0;
}
/* ---------------------------------------------------------------------- */
Neighbor::~Neighbor()
{
if (copymode) return;
memory->destroy(cutneighsq);
memory->destroy(cutneighghostsq);
delete [] cuttype;
delete [] cuttypesq;
delete [] fixchecklist;
memory->destroy(xhold);
memory->destroy(binhead);
memory->destroy(bins);
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;
for (int i = 0; i < nlist; i++) delete lists[i];
delete [] lists;
delete [] pair_build;
delete [] stencil_create;
delete [] blist;
delete [] glist;
delete [] slist;
for (int i = 0; i < nrequest; i++) delete requests[i];
memory->sfree(requests);
for (int i = 0; i < old_nrequest; i++) delete old_requests[i];
memory->sfree(old_requests);
delete [] zeroes;
memory->destroy(bondlist);
memory->destroy(anglelist);
memory->destroy(dihedrallist);
memory->destroy(improperlist);
}
/* ---------------------------------------------------------------------- */
void Neighbor::init()
{
int i,j,m,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 = 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 = NULL 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 == NULL) {
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;
// check other classes that can induce reneighboring in decide()
restart_check = 0;
if (output->restart_flag) restart_check = 1;
delete [] fixchecklist;
fixchecklist = NULL;
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;
if (force->kspace || force->pair_match("coul/wolf",0) ||
force->pair_match("coul/dsf",0) || force->pair_match("thole",0))
special_flag[1] = special_flag[2] = special_flag[3] = 2;
// maxwt = max multiplicative factor on atom indices stored in neigh list
maxwt = 0;
if (special_flag[1] == 2) maxwt = 2;
if (special_flag[2] == 2) maxwt = 3;
if (special_flag[3] == 2) maxwt = 4;
// rRESPA cutoffs
int respa = 0;
if (update->whichflag == 1 && strstr(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;
}
// ------------------------------------------------------------------
// xhold, bins, exclusion lists
// free xhold and bins if not needed for this run
if (dist_check == 0) {
memory->destroy(xhold);
maxhold = 0;
xhold = NULL;
}
if (style == NSQ) {
memory->destroy(bins);
memory->destroy(binhead);
maxbin = maxhead = 0;
binhead = NULL;
bins = NULL;
}
// 1st time allocation of xhold and bins
if (dist_check) {
if (maxhold == 0) {
maxhold = atom->nmax;
memory->create(xhold,maxhold,3,"neigh:xhold");
}
}
if (style != NSQ) {
if (maxbin == 0) {
maxbin = atom->nmax;
memory->create(bins,maxbin,"bins");
}
}
// exclusion lists for 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");
// ------------------------------------------------------------------
// pairwise lists
// 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
// first archive request params for current requests
// before Neighbor possibly changes them below
for (i = 0; i < nrequest; i++) requests[i]->archive();
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 old and new are not the same, create new pairwise lists
if (!same) {
// delete old lists and create new ones
for (i = 0; i < nlist; i++) delete lists[i];
delete [] lists;
delete [] pair_build;
delete [] stencil_create;
if (lmp->kokkos) nlist = init_lists_kokkos();
else nlist = nrequest;
lists = new NeighList*[nrequest];
pair_build = new PairPtr[nrequest];
stencil_create = new StencilPtr[nrequest];
// initialize to NULL since some may be Kokkos lists
for (i = 0; i < nrequest; i++) {
lists[i] = NULL;
pair_build[i] = NULL;
stencil_create[i] = NULL;
}
// create individual lists, one per request
// pass list ptr back to requestor (except for Command class)
// wait to allocate initial pages until copy lists are detected
for (i = 0; i < nrequest; i++) {
if (requests[i]->kokkos_host || requests[i]->kokkos_device) continue;
lists[i] = new NeighList(lmp);
lists[i]->index = i;
if (requests[i]->pair) {
Pair *pair = (Pair *) requests[i]->requestor;
pair->init_list(requests[i]->id,lists[i]);
} else if (requests[i]->fix) {
Fix *fix = (Fix *) requests[i]->requestor;
fix->init_list(requests[i]->id,lists[i]);
} else if (requests[i]->compute) {
Compute *compute = (Compute *) requests[i]->requestor;
compute->init_list(requests[i]->id,lists[i]);
}
}
// detect lists that are connected to other lists
// if-then-else sequence and processed flag is important
// since don't want to re-process skip or copy lists further down
int processed;
for (i = 0; i < nrequest; i++) {
if (!lists[i]) continue;
processed = 0;
// copy: point this list at request->otherlist, could be a skip list
if (requests[i]->copy) {
lists[i]->listcopy = lists[requests[i]->otherlist];
processed = 1;
// skip: point this list at request->otherlist,
// copy skip info from request
// skip list still needs to have granhistory or respa info added below
} else if (requests[i]->skip) {
lists[i]->listskip = lists[requests[i]->otherlist];
lists[i]->copy_skip_info(requests[i]->iskip,requests[i]->ijskip);
processed = 1;
// half_from_full: point this list at full list that comes right before
// will only be case if pair style requested one after other
} else if (requests[i]->half_from_full) {
lists[i]->listfull = lists[i-1];
processed = 1;
}
// granhistory: set preceeding list's listgranhistory to this list
// also set FH ptr in preceeding list to FSH class created by pair
if (requests[i]->granhistory) {
lists[i-1]->listgranhistory = lists[i];
lists[i-1]->fix_history = requests[i]->fix_history;
processed = 1;
// respaouter: point this list at preceeding 1/2 inner/middle lists
} else if (requests[i]->respaouter) {
if (requests[i-1]->respainner) {
lists[i]->respamiddle = 0;
lists[i]->listinner = lists[i-1];
} else {
lists[i]->respamiddle = 1;
lists[i]->listmiddle = lists[i-1];
lists[i]->listinner = lists[i-2];
}
processed = 1;
}
if (processed) continue;
// pair and half and newton != 2:
// if there is a full non-occasional non-skip list
// change this list to half_from_full and point at the full list
// parent could be copy list or pair or fix
// could remove newton != 2 check if added half_from_full_no_newton_ghost
// option in neigh_derive.cpp and below in choose_build()
// this would require full list had ghost info
// would be useful when reax/c used in hybrid mode, e.g. with airebo
if (requests[i]->pair && requests[i]->half && requests[i]->newton != 2) {
for (j = 0; j < nrequest; j++) {
if (!lists[j]) continue;
if (requests[j]->full && requests[j]->occasional == 0 &&
requests[j]->skip == 0) break;
}
if (j < nrequest) {
requests[i]->half = 0;
requests[i]->half_from_full = 1;
lists[i]->listfull = lists[j];
}
// fix/compute requests:
// whether request is occasional or not doesn't matter
// if request = half and non-skip pair half/respaouter exists,
// if request = full and non-skip pair full exists,
// if request = half and non-skip pair full exists,
// if no matches, do nothing
// fix/compute list will be built independently as needed
// ok if parent is itself a copy list
} else if (requests[i]->fix || requests[i]->compute) {
for (j = 0; j < nrequest; j++) {
if (!lists[j]) continue;
if (requests[i]->half && requests[j]->pair &&
requests[j]->skip == 0 && requests[j]->half) break;
if (requests[i]->full && requests[j]->pair &&
requests[j]->skip == 0 && requests[j]->full) break;
if (requests[i]->gran && requests[j]->pair &&
requests[j]->skip == 0 && requests[j]->gran) break;
if (requests[i]->half && requests[j]->pair &&
requests[j]->skip == 0 && requests[j]->respaouter) break;
}
if (j < nrequest) {
requests[i]->copy = 1;
requests[i]->otherlist = j;
lists[i]->listcopy = lists[j];
} else {
for (j = 0; j < nrequest; j++) {
if (!lists[j]) continue;
if (requests[i]->half && requests[j]->pair &&
requests[j]->skip == 0 && requests[j]->full) break;
}
if (j < nrequest) {
requests[i]->half = 0;
requests[i]->half_from_full = 1;
lists[i]->listfull = lists[j];
}
}
}
}
// allocate initial pages for each list, except if listcopy set
// allocate dnum vector of zeroes if set
int dnummax = 0;
for (i = 0; i < nrequest; i++) {
if (!lists[i]) continue;
if (!lists[i]->listcopy) {
lists[i]->setup_pages(pgsize,oneatom,requests[i]->dnum);
dnummax = MAX(dnummax,requests[i]->dnum);
}
}
if (dnummax) {
delete [] zeroes;
zeroes = new double[dnummax];
for (i = 0; i < dnummax; i++) zeroes[i] = 0.0;
}
// set ptrs to pair_build and stencil_create functions for each list
// ptrs set to NULL if not set explicitly
for (i = 0; i < nrequest; i++) {
choose_build(i,requests[i]);
if (style != NSQ) choose_stencil(i,requests[i]);
else stencil_create[i] = NULL;
}
// set each list's build/grow/stencil/ghost flags based on neigh request
// buildflag = 1 if its pair_build() invoked every reneighbor
// growflag = 1 if it stores atom-based arrays and pages
// stencilflag = 1 if it stores stencil arrays
// ghostflag = 1 if it stores neighbors of ghosts
// anyghostlist = 1 if any non-occasional list stores neighbors of ghosts
anyghostlist = 0;
int anybuild = 0;
for (i = 0; i < nrequest; i++) {
if (lists[i]) {
lists[i]->buildflag = 1;
if (pair_build[i] == NULL) lists[i]->buildflag = 0;
if (requests[i]->occasional) lists[i]->buildflag = 0;
if (lists[i]->buildflag) anybuild = 1;
lists[i]->growflag = 1;
if (requests[i]->copy) lists[i]->growflag = 0;
lists[i]->stencilflag = 1;
if (style == NSQ) lists[i]->stencilflag = 0;
if (stencil_create[i] == NULL) lists[i]->stencilflag = 0;
lists[i]->ghostflag = 0;
if (requests[i]->ghost) lists[i]->ghostflag = 1;
if (requests[i]->ghost && !requests[i]->occasional) anyghostlist = 1;
lists[i]->ssaflag = 0;
if (requests[i]->ssa) lists[i]->ssaflag = 1;
} else init_list_flags1_kokkos(i);
}
// no request has the buildflag set, so set it for the first request only
// this insure binning is done for any occasional neighbor lists
if (!anybuild) {
for (i = 0; i < nrequest; i++) {
if (lists[i]) {
lists[i]->buildflag = 1;
break;
}
}
}
#ifdef NEIGH_LIST_DEBUG
for (i = 0; i < nrequest; i++)
if (comm->me == 0) printf("Build/stencil methods: %d: %p %p\n",
i,pair_build[i],stencil_create[i]);
for (i = 0; i < nrequest; i++) lists[i]->print_attributes();
#endif
// allocate atom arrays for neighbor lists that store them
maxatom = atom->nmax;
for (i = 0; i < nrequest; i++) {
if (lists[i]) {
if (lists[i]->growflag) lists[i]->grow(maxatom);
} else init_list_grow_kokkos(i);
}
// setup 3 vectors of pairwise neighbor lists
// blist = lists whose pair_build() is invoked every reneighbor
// glist = lists who store atom arrays which are used every reneighbor
// slist = lists who store stencil arrays which are used every reneighbor
// blist and glist vectors are used by neighbor::build()
// slist vector is used by neighbor::setup_bins()
nblist = nglist = nslist = 0;
delete [] blist;
delete [] glist;
delete [] slist;
blist = new int[nrequest];
glist = new int[nrequest];
slist = new int[nrequest];
for (i = 0; i < nrequest; i++) {
if (lists[i]) {
if (lists[i]->buildflag) blist[nblist++] = i;
if (lists[i]->growflag && requests[i]->occasional == 0)
glist[nglist++] = i;
if (lists[i]->stencilflag && requests[i]->occasional == 0)
slist[nslist++] = i;
} else init_list_flags2_kokkos(i);
}
// no request had buildflag set, so we set it for the first request
// we also need to set other occasional neighbor list properties
if (!anybuild)
for (i = 0; i < nrequest; i++)
if (lists[i] && lists[i]->buildflag) {
if (lists[i]->growflag) glist[nglist++] = i;
if (lists[i]->stencilflag) slist[nslist++] = i;
}
#ifdef NEIGH_LIST_DEBUG
print_lists_of_lists();
#endif
// reorder build vector if necessary
// relevant for lists that copy/skip/half-full from parent
// the derived list must appear in blist after the parent list
// no occasional lists are in build vector
// swap two lists within blist when dependency is mis-ordered
// done when entire pass thru blist results in no swaps
int done = 0;
while (!done) {
done = 1;
for (i = 0; i < nblist; i++) {
if (!lists[blist[i]]) continue;
NeighList *ptr = NULL;
if (lists[blist[i]]->listfull) ptr = lists[blist[i]]->listfull;
if (lists[blist[i]]->listcopy) ptr = lists[blist[i]]->listcopy;
if (lists[blist[i]]->listskip) ptr = lists[blist[i]]->listskip;
if (ptr == NULL) continue;
for (m = 0; m < nrequest; m++)
if (ptr == lists[m]) break;
for (j = 0; j < nblist; j++)
if (m == blist[j]) break;
if (j < i) continue;
int tmp = blist[i];
blist[i] = blist[j];
blist[j] = tmp;
done = 0;
break;
}
}
#ifdef NEIGH_LIST_DEBUG
print_lists_of_lists();
#endif
}
// output neighbor list info, only first time or when info changes
if (!same || every != old_every || delay != old_delay ||
old_check != dist_check || old_cutoff != cutneighmax) {
if (me == 0) {
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 == BIN) binsize = 0.5*cutneighmax;
else binsize = 0.5*cutneighmin;
if (binsize == 0.0) binsize = bbox[0];
if (logfile) {
fprintf(logfile,"Neighbor list info ...\n");
fprintf(logfile," %d neighbor list requests\n",nrequest);
fprintf(logfile," update every %d steps, delay %d steps, check %s\n",
every,delay,dist_check ? "yes" : "no");
fprintf(logfile," max neighbors/atom: %d, page size: %d\n",
oneatom, pgsize);
fprintf(logfile," master list distance cutoff = %g\n",cutneighmax);
fprintf(logfile," ghost atom cutoff = %g\n",cutghost);
if (style != NSQ)
fprintf(logfile," binsize = %g -> bins = %g %g %g\n",binsize,
ceil(bbox[0]/binsize), ceil(bbox[1]/binsize),
ceil(bbox[2]/binsize));
}
if (screen) {
fprintf(screen,"Neighbor list info ...\n");
fprintf(screen," %d neighbor list requests\n",nrequest);
fprintf(screen," update every %d steps, delay %d steps, check %s\n",
every,delay,dist_check ? "yes" : "no");
fprintf(screen," max neighbors/atom: %d, page size: %d\n",
oneatom, pgsize);
fprintf(screen," master list distance cutoff = %g\n",cutneighmax);
fprintf(screen," ghost atom cutoff = %g\n",cutghost);
if (style != NSQ)
fprintf(screen," binsize = %g, bins = %g %g %g\n",binsize,
ceil(bbox[0]/binsize), ceil(bbox[1]/binsize),
ceil(bbox[2]/binsize));
}
}
}
// mark all current requests as processed
// delete old requests
// copy current requests and style to old for next run
for (i = 0; i < nrequest; i++) requests[i]->unprocessed = 0;
for (i = 0; i < old_nrequest; i++) delete old_requests[i];
memory->sfree(old_requests);
old_nrequest = nrequest;
old_requests = requests;
nrequest = maxrequest = 0;
requests = NULL;
old_style = style;
old_triclinic = triclinic;
old_pgsize = pgsize;
old_oneatom = oneatom;
old_every = every;
old_delay = delay;
old_check = dist_check;
old_cutoff = cutneighmax;
// ------------------------------------------------------------------
// topology lists
// 1st time allocation of topology lists
if (lmp->kokkos) {
init_topology_kokkos();
return;
}
if (atom->molecular && atom->nbonds && maxbond == 0) {
if (nprocs == 1) maxbond = atom->nbonds;
else maxbond = static_cast<int> (LB_FACTOR * atom->nbonds / nprocs);
memory->create(bondlist,maxbond,3,"neigh:bondlist");
}
if (atom->molecular && atom->nangles && maxangle == 0) {
if (nprocs == 1) maxangle = atom->nangles;
else maxangle = static_cast<int> (LB_FACTOR * atom->nangles / nprocs);
memory->create(anglelist,maxangle,4,"neigh:anglelist");
}
if (atom->molecular && atom->ndihedrals && maxdihedral == 0) {
if (nprocs == 1) maxdihedral = atom->ndihedrals;
else maxdihedral = static_cast<int>
(LB_FACTOR * atom->ndihedrals / nprocs);
memory->create(dihedrallist,maxdihedral,5,"neigh:dihedrallist");
}
if (atom->molecular && atom->nimpropers && maximproper == 0) {
if (nprocs == 1) maximproper = atom->nimpropers;
else maximproper = static_cast<int>
(LB_FACTOR * atom->nimpropers / nprocs);
memory->create(improperlist,maximproper,5,"neigh:improperlist");
}
// set flags that determine which topology neighboring routines to use
// bonds,etc can only be broken for atom->molecular = 1, not 2
// 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 ((strcmp(modify->fix[i]->style,"shake") == 0)
|| (strcmp(modify->fix[i]->style,"rattle") == 0))
bond_off = angle_off = 1;
if (force->bond && force->bond_match("quartic")) bond_off = 1;
if (atom->avec->bonds_allow && atom->molecular == 1) {
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 == 1) {
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 == 1) {
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 == 1) {
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 on_or_off = bond_off;
MPI_Allreduce(&on_or_off,&bond_off,1,MPI_INT,MPI_MAX,world);
on_or_off = angle_off;
MPI_Allreduce(&on_or_off,&angle_off,1,MPI_INT,MPI_MAX,world);
on_or_off = dihedral_off;
MPI_Allreduce(&on_or_off,&dihedral_off,1,MPI_INT,MPI_MAX,world);
on_or_off = improper_off;
MPI_Allreduce(&on_or_off,&improper_off,1,MPI_INT,MPI_MAX,world);
// set ptrs to topology build functions
if (atom->molecular == 2) bond_build = &Neighbor::bond_template;
else if (bond_off) bond_build = &Neighbor::bond_partial;
else bond_build = &Neighbor::bond_all;
if (atom->molecular == 2) angle_build = &Neighbor::angle_template;
else if (angle_off) angle_build = &Neighbor::angle_partial;
else angle_build = &Neighbor::angle_all;
if (atom->molecular == 2) dihedral_build = &Neighbor::dihedral_template;
else if (dihedral_off) dihedral_build = &Neighbor::dihedral_partial;
else dihedral_build = &Neighbor::dihedral_all;
if (atom->molecular == 2) improper_build = &Neighbor::improper_template;
else if (improper_off) improper_build = &Neighbor::improper_partial;
else improper_build = &Neighbor::improper_all;
// set topology neighbor list counts to 0
// in case all are turned off but potential is still defined
nbondlist = nanglelist = ndihedrallist = nimproperlist = 0;
}
/* ---------------------------------------------------------------------- */
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);
requests[nrequest]->requestor = requestor;
requests[nrequest]->requestor_instance = instance;
nrequest++;
return nrequest-1;
}
/* ----------------------------------------------------------------------
determine which pair_build function each neigh list needs
based on settings of neigh request
copy -> copy_from function
skip -> granular function if gran, several options
respa function if respaouter,
skip_from function for everything else
ssa -> special case for USER-DPD pair styles
half_from_full, half, full, gran, respaouter ->
choose by newton and rq->newton and triclinic settings
style NSQ options = newton off, newton on
style BIN options = newton off, newton on and not tri, newton on and tri
stlye MULTI options = same options as BIN
if none of these, ptr = NULL since pair_build is not invoked for this list
use "else if" logic b/c skip,copy can be set in addition to half,full,etc
------------------------------------------------------------------------- */
void Neighbor::choose_build(int index, NeighRequest *rq)
{
PairPtr pb = NULL;
if (rq->omp == 0 && rq->intel == 0) {
if (rq->copy) pb = &Neighbor::copy_from;
else if (rq->skip) {
if (rq->gran) {
NeighRequest *otherrq = requests[rq->otherlist];
if (otherrq->newton == 0) {
pb = &Neighbor::skip_from_granular;
} else if (otherrq->newton == 1) {
error->all(FLERR,"Neighbor build method not supported");
} else if (otherrq->newton == 2) {
if (rq->granonesided == 0)
pb = &Neighbor::skip_from_granular_off2on;
else if (rq->granonesided == 1)
pb = &Neighbor::skip_from_granular_off2on_onesided;
}
}
else if (rq->respaouter) pb = &Neighbor::skip_from_respa;
else pb = &Neighbor::skip_from;
} else if (rq->ssa) {
if (rq->half_from_full) pb = &Neighbor::half_from_full_newton_ssa;
else pb = &Neighbor::half_bin_newton_ssa;
} else if (rq->half_from_full) {
if (rq->newton == 0) {
if (newton_pair == 0) pb = &Neighbor::half_from_full_no_newton;
else if (newton_pair == 1) pb = &Neighbor::half_from_full_newton;
} else if (rq->newton == 1) {
pb = &Neighbor::half_from_full_newton;
} else if (rq->newton == 2) {
pb = &Neighbor::half_from_full_no_newton;
}
} else if (rq->half) {
if (style == NSQ) {
if (rq->newton == 0) {
if (newton_pair == 0) {
if (rq->ghost == 0) pb = &Neighbor::half_nsq_no_newton;
else if (includegroup)
error->all(FLERR,"Neighbor include group not allowed "
"with ghost neighbors");
else pb = &Neighbor::half_nsq_no_newton_ghost;
} else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton;
} else if (rq->newton == 1) {
pb = &Neighbor::half_nsq_newton;
} else if (rq->newton == 2) {
if (rq->ghost == 0) pb = &Neighbor::half_nsq_no_newton;
else if (includegroup)
error->all(FLERR,"Neighbor include group not allowed "
"with ghost neighbors");
else pb = &Neighbor::half_nsq_no_newton_ghost;
}
} else if (style == BIN) {
if (rq->newton == 0) {
if (newton_pair == 0) {
if (rq->ghost == 0) pb = &Neighbor::half_bin_no_newton;
else if (includegroup)
error->all(FLERR,"Neighbor include group not allowed "
"with ghost neighbors");
else pb = &Neighbor::half_bin_no_newton_ghost;
} else if (triclinic == 0) {
pb = &Neighbor::half_bin_newton;
} else if (triclinic == 1)
pb = &Neighbor::half_bin_newton_tri;
} else if (rq->newton == 1) {
if (triclinic == 0) pb = &Neighbor::half_bin_newton;
else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri;
} else if (rq->newton == 2) {
if (rq->ghost == 0) pb = &Neighbor::half_bin_no_newton;
else if (includegroup)
error->all(FLERR,"Neighbor include group not allowed "
"with ghost neighbors");
else pb = &Neighbor::half_bin_no_newton_ghost;
}
} else if (style == MULTI) {
if (rq->ghost == 1)
error->all(FLERR,
"Neighbor multi not yet enabled for ghost neighbors");
if (rq->newton == 0) {
if (newton_pair == 0) pb = &Neighbor::half_multi_no_newton;
else if (triclinic == 0) pb = &Neighbor::half_multi_newton;
else if (triclinic == 1) pb = &Neighbor::half_multi_newton_tri;
} else if (rq->newton == 1) {
if (triclinic == 0) pb = &Neighbor::half_multi_newton;
else if (triclinic == 1) pb = &Neighbor::half_multi_newton_tri;
} else if (rq->newton == 2) pb = &Neighbor::half_multi_no_newton;
}
} else if (rq->full) {
if (style == NSQ) {
if (rq->ghost == 0) pb = &Neighbor::full_nsq;
else if (includegroup)
error->all(FLERR,
"Neighbor include group not allowed with ghost neighbors");
else pb = &Neighbor::full_nsq_ghost;
} else if (style == BIN) {
if (rq->ghost == 0) pb = &Neighbor::full_bin;
else if (includegroup)
error->all(FLERR,
"Neighbor include group not allowed with ghost neighbors");
else pb = &Neighbor::full_bin_ghost;
} else if (style == MULTI) {
if (rq->ghost == 1)
error->all(FLERR,
"Neighbor multi not yet enabled for ghost neighbors");
pb = &Neighbor::full_multi;
}
} else if (rq->gran) {
if (rq->newton == 0) {
if (style == NSQ) {
if (newton_pair == 0) pb = &Neighbor::granular_nsq_no_newton;
else if (newton_pair == 1) {
if (rq->granonesided == 0) pb = &Neighbor::granular_nsq_newton;
else pb = &Neighbor::granular_nsq_newton_onesided;
}
} else if (style == BIN) {
if (newton_pair == 0) pb = &Neighbor::granular_bin_no_newton;
else if (newton_pair == 1) {
if (triclinic == 0) {
if (rq->granonesided == 0) pb = &Neighbor::granular_bin_newton;
else pb = &Neighbor::granular_bin_newton_onesided;
} else if (triclinic == 1) {
if (rq->granonesided == 0)
pb = &Neighbor::granular_bin_newton_tri;
else error->all(FLERR,"Neighbor build method not supported");
}
}
} else if (style == MULTI)
error->all(FLERR,"Neighbor multi not yet enabled for granular");
} else if (rq->newton == 1) {
error->all(FLERR,"Neighbor build method not yet supported");
} else if (rq->newton == 2) {
if (style == NSQ) pb = &Neighbor::granular_nsq_no_newton;
else if (style == BIN) {
if (triclinic == 0) pb = &Neighbor::granular_bin_no_newton;
else if (triclinic == 1)
error->all(FLERR,"Neighbor build method not yet supported");
} else if (style == MULTI)
error->all(FLERR,"Neighbor multi not yet enabled for granular");
}
} else if (rq->respaouter) {
if (style == NSQ) {
if (newton_pair == 0) pb = &Neighbor::respa_nsq_no_newton;
else if (newton_pair == 1) pb = &Neighbor::respa_nsq_newton;
} else if (style == BIN) {
if (newton_pair == 0) pb = &Neighbor::respa_bin_no_newton;
else if (triclinic == 0) pb = &Neighbor::respa_bin_newton;
else if (triclinic == 1) pb = &Neighbor::respa_bin_newton_tri;
} else if (style == MULTI)
error->all(FLERR,"Neighbor multi not yet enabled for rRESPA");
}
// OMP versions of build methods
} else {
if (rq->copy) pb = &Neighbor::copy_from;
else if (rq->skip) {
if (rq->gran && lists[index]->listgranhistory)
pb = &Neighbor::skip_from_granular;
else if (rq->respaouter) pb = &Neighbor::skip_from_respa;
else pb = &Neighbor::skip_from;
} else if (rq->half_from_full) {
if (newton_pair == 0) pb = &Neighbor::half_from_full_no_newton_omp;
else if (newton_pair == 1) pb = &Neighbor::half_from_full_newton_omp;
} else if (rq->half) {
if (style == NSQ) {
if (rq->newton == 0) {
if (newton_pair == 0) {
if (rq->ghost == 0) pb = &Neighbor::half_nsq_no_newton_omp;
else if (includegroup)
error->all(FLERR,"Neighbor include group not allowed "
"with ghost neighbors");
else pb = &Neighbor::half_nsq_no_newton_ghost_omp;
} else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton_omp;
} else if (rq->newton == 1) {
pb = &Neighbor::half_nsq_newton_omp;
} else if (rq->newton == 2) {
if (rq->ghost == 0) pb = &Neighbor::half_nsq_no_newton_omp;
else if (includegroup)
error->all(FLERR,"Neighbor include group not allowed "
"with ghost neighbors");
else pb = &Neighbor::half_nsq_no_newton_ghost_omp;
}
} else if (style == BIN) {
if (rq->newton == 0) {
if (newton_pair == 0) {
if (rq->ghost == 0) {
if (rq->intel) pb = &Neighbor::half_bin_no_newton_intel;
else pb = &Neighbor::half_bin_no_newton_omp;
} else if (includegroup)
error->all(FLERR,"Neighbor include group not allowed "
"with ghost neighbors");
else pb = &Neighbor::half_bin_no_newton_ghost_omp;
} else if (triclinic == 0) {
if (rq->intel) pb = &Neighbor::half_bin_newton_intel;
else pb = &Neighbor::half_bin_newton_omp;
} else if (triclinic == 1) {
if (rq->intel) pb = &Neighbor::half_bin_newton_tri_intel;
else pb = &Neighbor::half_bin_newton_tri_omp;
}
} else if (rq->newton == 1) {
if (triclinic == 0) {
if (rq->intel) pb = &Neighbor::half_bin_newton_intel;
else pb = &Neighbor::half_bin_newton_omp;
} else if (triclinic == 1) {
if (rq->intel) pb = &Neighbor::half_bin_newton_tri_intel;
else pb = &Neighbor::half_bin_newton_tri_omp;
}
} else if (rq->newton == 2) {
if (rq->ghost == 0) {
if (rq->intel) pb = &Neighbor::half_bin_no_newton_intel;
else pb = &Neighbor::half_bin_no_newton_omp;
} else if (includegroup)
error->all(FLERR,"Neighbor include group not allowed "
"with ghost neighbors");
else pb = &Neighbor::half_bin_no_newton_ghost_omp;
}
} else if (style == MULTI) {
if (rq->ghost == 1)
error->all(FLERR,
"Neighbor multi not yet enabled for ghost neighbors");
if (rq->newton == 0) {
if (newton_pair == 0) pb = &Neighbor::half_multi_no_newton_omp;
else if (triclinic == 0) pb = &Neighbor::half_multi_newton_omp;
else if (triclinic == 1) pb = &Neighbor::half_multi_newton_tri_omp;
} else if (rq->newton == 1) {
if (triclinic == 0) pb = &Neighbor::half_multi_newton_omp;
else if (triclinic == 1) pb = &Neighbor::half_multi_newton_tri_omp;
} else if (rq->newton == 2) pb = &Neighbor::half_multi_no_newton_omp;
}
} else if (rq->full) {
if (style == NSQ) {
if (rq->ghost == 0) pb = &Neighbor::full_nsq_omp;
else if (includegroup)
error->all(FLERR,
"Neighbor include group not allowed with ghost neighbors");
else pb = &Neighbor::full_nsq_ghost_omp;
} else if (style == BIN) {
if (rq->ghost == 0) {
if (rq->intel) pb = &Neighbor::full_bin_intel;
else pb = &Neighbor::full_bin_omp;
} else if (includegroup)
error->all(FLERR,
"Neighbor include group not allowed with ghost neighbors");
else pb = &Neighbor::full_bin_ghost_omp;
} else if (style == MULTI) {
if (rq->ghost == 1)
error->all(FLERR,
"Neighbor multi not yet enabled for ghost neighbors");
pb = &Neighbor::full_multi_omp;
}
} else if (rq->gran) {
if (style == NSQ) {
if (newton_pair == 0) pb = &Neighbor::granular_nsq_no_newton_omp;
else if (newton_pair == 1) pb = &Neighbor::granular_nsq_newton_omp;
} else if (style == BIN) {
if (newton_pair == 0) pb = &Neighbor::granular_bin_no_newton_omp;
else if (triclinic == 0) pb = &Neighbor::granular_bin_newton_omp;
else if (triclinic == 1) pb = &Neighbor::granular_bin_newton_tri_omp;
} else if (style == MULTI)
error->all(FLERR,"Neighbor multi not yet enabled for granular");
} else if (rq->respaouter) {
if (style == NSQ) {
if (newton_pair == 0) pb = &Neighbor::respa_nsq_no_newton_omp;
else if (newton_pair == 1) pb = &Neighbor::respa_nsq_newton_omp;
} else if (style == BIN) {
if (newton_pair == 0) pb = &Neighbor::respa_bin_no_newton_omp;
else if (triclinic == 0) pb = &Neighbor::respa_bin_newton_omp;
else if (triclinic == 1) pb = &Neighbor::respa_bin_newton_tri_omp;
} else if (style == MULTI)
error->all(FLERR,"Neighbor multi not yet enabled for rRESPA");
}
}
pair_build[index] = pb;
}
/* ----------------------------------------------------------------------
determine which stencil_create function each neigh list needs
based on settings of neigh request, only called if style != NSQ
skip or copy or half_from_full -> no stencil
ssa = special case for USER-DPD pair styles
half, gran, respaouter, full -> choose by newton and tri and dimension
if none of these, ptr = NULL since this list needs no stencils
use "else if" b/c skip,copy can be set in addition to half,full,etc
------------------------------------------------------------------------- */
void Neighbor::choose_stencil(int index, NeighRequest *rq)
{
StencilPtr sc = NULL;
if (rq->skip || rq->copy || rq->half_from_full) sc = NULL;
else if (rq->ssa) {
if (dimension == 2) sc = &Neighbor::stencil_half_bin_2d_ssa;
else if (dimension == 3) sc = &Neighbor::stencil_half_bin_3d_ssa;
} else if (rq->half || rq->gran || rq->respaouter) {
if (style == BIN) {
if (rq->newton == 0) {
if (newton_pair == 0) {
if (dimension == 2) {
if (rq->ghost) sc = &Neighbor::stencil_half_ghost_bin_2d_no_newton;
else sc = &Neighbor::stencil_half_bin_2d_no_newton;
} else if (dimension == 3) {
if (rq->ghost) sc = &Neighbor::stencil_half_ghost_bin_3d_no_newton;
else sc = &Neighbor::stencil_half_bin_3d_no_newton;
}
} else if (triclinic == 0) {
if (dimension == 2)
sc = &Neighbor::stencil_half_bin_2d_newton;
else if (dimension == 3)
sc = &Neighbor::stencil_half_bin_3d_newton;
} else if (triclinic == 1) {
if (dimension == 2)
sc = &Neighbor::stencil_half_bin_2d_newton_tri;
else if (dimension == 3)
sc = &Neighbor::stencil_half_bin_3d_newton_tri;
}
} else if (rq->newton == 1) {
if (triclinic == 0) {
if (dimension == 2)
sc = &Neighbor::stencil_half_bin_2d_newton;
else if (dimension == 3)
sc = &Neighbor::stencil_half_bin_3d_newton;
} else if (triclinic == 1) {
if (dimension == 2)
sc = &Neighbor::stencil_half_bin_2d_newton_tri;
else if (dimension == 3)
sc = &Neighbor::stencil_half_bin_3d_newton_tri;
}
} else if (rq->newton == 2) {
if (dimension == 2)
if (rq->ghost) sc = &Neighbor::stencil_half_ghost_bin_2d_no_newton;
else sc = &Neighbor::stencil_half_bin_2d_no_newton;
else if (dimension == 3) {
if (rq->ghost) sc = &Neighbor::stencil_half_ghost_bin_3d_no_newton;
else sc = &Neighbor::stencil_half_bin_3d_no_newton;
}
}
} else if (style == MULTI) {
if (rq->newton == 0) {
if (newton_pair == 0) {
if (dimension == 2)
sc = &Neighbor::stencil_half_multi_2d_no_newton;
else if (dimension == 3)
sc = &Neighbor::stencil_half_multi_3d_no_newton;
} else if (triclinic == 0) {
if (dimension == 2)
sc = &Neighbor::stencil_half_multi_2d_newton;
else if (dimension == 3)
sc = &Neighbor::stencil_half_multi_3d_newton;
} else if (triclinic == 1) {
if (dimension == 2)
sc = &Neighbor::stencil_half_multi_2d_newton_tri;
else if (dimension == 3)
sc = &Neighbor::stencil_half_multi_3d_newton_tri;
}
} else if (rq->newton == 1) {
if (triclinic == 0) {
if (dimension == 2)
sc = &Neighbor::stencil_half_multi_2d_newton;
else if (dimension == 3)
sc = &Neighbor::stencil_half_multi_3d_newton;
} else if (triclinic == 1) {
if (dimension == 2)
sc = &Neighbor::stencil_half_multi_2d_newton_tri;
else if (dimension == 3)
sc = &Neighbor::stencil_half_multi_3d_newton_tri;
}
} else if (rq->newton == 2) {
if (dimension == 2)
sc = &Neighbor::stencil_half_multi_2d_no_newton;
else if (dimension == 3)
sc = &Neighbor::stencil_half_multi_3d_no_newton;
}
}
} else if (rq->full) {
if (style == BIN) {
if (dimension == 2) {
if (rq->ghost) sc = &Neighbor::stencil_full_ghost_bin_2d;
else sc = &Neighbor::stencil_full_bin_2d;
}
else if (dimension == 3) {
if (rq->ghost) sc = &Neighbor::stencil_full_ghost_bin_3d;
else sc = &Neighbor::stencil_full_bin_3d;
}
} else if (style == MULTI) {
if (dimension == 2) sc = &Neighbor::stencil_full_multi_2d;
else if (dimension == 3) sc = &Neighbor::stencil_full_multi_3d;
}
}
stencil_create[index] = sc;
}
/* ---------------------------------------------------------------------- */
void Neighbor::print_lists_of_lists()
{
if (comm->me == 0) {
printf("Build lists = %d: ",nblist);
for (int i = 0; i < nblist; i++) printf("%d ",blist[i]);
printf("\n");
printf("Grow lists = %d: ",nglist);
for (int i = 0; i < nglist; i++) printf("%d ",glist[i]);
printf("\n");
printf("Stencil lists = %d: ",nslist);
for (int i = 0; i < nslist; i++) printf("%d ",slist[i]);
printf("\n");
}
}
/* ---------------------------------------------------------------------- */
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
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));
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));
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, USER-CUDA calls with topoflag = 0
------------------------------------------------------------------------- */
void Neighbor::build(int topoflag)
{
int i;
ago = 0;
ncalls++;
lastcall = update->ntimestep;
// store current atom positions and box size if needed
if (dist_check) {
double **x = atom->x;
int nlocal = atom->nlocal;
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];
}
}
}
}
// if any lists store neighbors of ghosts:
// invoke grow() if nlocal+nghost exceeds previous list size
// else only invoke grow() if nlocal exceeds previous list size
// only for lists with growflag set and which are perpetual (glist)
if (anyghostlist && atom->nmax > maxatom) {
maxatom = atom->nmax;
for (i = 0; i < nglist; i++) lists[glist[i]]->grow(maxatom);
} else if (atom->nmax > maxatom) {
maxatom = atom->nmax;
for (i = 0; i < nglist; i++) lists[glist[i]]->grow(maxatom);
}
// extend atom bin list if necessary
if (style != NSQ && atom->nmax > maxbin) {
maxbin = atom->nmax;
memory->destroy(bins);
memory->create(bins,maxbin,"bins");
}
// check that using special bond flags will not overflow neigh lists
if (atom->nlocal+atom->nghost > NEIGHMASK)
error->one(FLERR,"Too many local+ghost atoms for neighbor list");
// invoke building of pair and molecular topology neighbor lists
// only for pairwise lists with buildflag set
// blist is for standard neigh lists, otherwise is a Kokkos list
for (i = 0; i < nblist; i++)
(this->*pair_build[blist[i]])(lists[blist[i]]);
if (atom->molecular && topoflag) build_topology();
}
/* ----------------------------------------------------------------------
build all topology neighbor lists every few timesteps
normally built with pair lists, but USER-CUDA separates them
------------------------------------------------------------------------- */
void Neighbor::build_topology()
{
if (force->bond) (this->*bond_build)();
if (force->angle) (this->*angle_build)();
if (force->dihedral) (this->*dihedral_build)();
if (force->improper) (this->*improper_build)();
}
/* ----------------------------------------------------------------------
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 == NULL)
error->all(FLERR,"Trying to build an occasional neighbor list "
"before initialization completed");
// 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 ">="
if (preflag) {
if (mylist->last_build > lastcall) return;
} else {
if (mylist->last_build >= lastcall) return;
}
mylist->last_build = update->ntimestep;
// update stencils and grow atom arrays as needed
// only for relevant settings of stencilflag and growflag
// grow atom array for this list to current size of perpetual lists
if (mylist->stencilflag) {
mylist->stencil_allocate(smax,style);
(this->*stencil_create[mylist->index])(mylist,sx,sy,sz);
}
if (mylist->growflag) mylist->grow(maxatom);
// build list I, turning off atom binning
// binning results from last re-neighbor should be used instead
// if re-bin now, atoms may have moved outside of proc domain & bin extent,
// leading to errors or even a crash
binatomflag = 0;
(this->*pair_build[mylist->index])(mylist);
binatomflag = 1;
}
/* ----------------------------------------------------------------------
setup neighbor binning parameters
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 = lowest global bin any of my ghost atoms could fall into
mbinhi = highest global bin any of my ghost atoms could fall into
mbin = number of bins I need in a dimension
------------------------------------------------------------------------- */
void Neighbor::setup_bins()
{
// bbox = size of bbox of entire domain
// bsubbox lo/hi = bounding box of my subdomain extended by comm->cutghost
// for triclinic:
// bbox bounds all 8 corners of tilted box
// subdomain is in lamda coords
// include dimension-dependent extension via comm->cutghost
// domain->bbox() converts lamda extent to box coords and computes bbox
double bbox[3],bsubboxlo[3],bsubboxhi[3];
double *cutghost = comm->cutghost;
if (triclinic == 0) {
bsubboxlo[0] = domain->sublo[0] - cutghost[0];
bsubboxlo[1] = domain->sublo[1] - cutghost[1];
bsubboxlo[2] = domain->sublo[2] - cutghost[2];
bsubboxhi[0] = domain->subhi[0] + cutghost[0];
bsubboxhi[1] = domain->subhi[1] + cutghost[1];
bsubboxhi[2] = domain->subhi[2] + cutghost[2];
} else {
double lo[3],hi[3];
lo[0] = domain->sublo_lamda[0] - cutghost[0];
lo[1] = domain->sublo_lamda[1] - cutghost[1];
lo[2] = domain->sublo_lamda[2] - cutghost[2];
hi[0] = domain->subhi_lamda[0] + cutghost[0];
hi[1] = domain->subhi_lamda[1] + cutghost[1];
hi[2] = domain->subhi_lamda[2] + cutghost[2];
domain->bbox(lo,hi,bsubboxlo,bsubboxhi);
}
bbox[0] = bboxhi[0] - bboxlo[0];
bbox[1] = bboxhi[1] - bboxlo[1];
bbox[2] = bboxhi[2] - bboxlo[2];
// optimal bin size is roughly 1/2 the cutoff
// for BIN style, binsize = 1/2 of max neighbor cutoff
// for MULTI style, binsize = 1/2 of min neighbor cutoff
// special case of all cutoffs = 0.0, binsize = box size
double binsize_optimal;
if (binsizeflag) binsize_optimal = binsize_user;
else if (style == BIN) binsize_optimal = 0.5*cutneighmax;
else binsize_optimal = 0.5*cutneighmin;
if (binsize_optimal == 0.0) binsize_optimal = bbox[0];
double binsizeinv = 1.0/binsize_optimal;
// test for too many global bins in any dimension due to huge global domain
if (bbox[0]*binsizeinv > MAXSMALLINT || bbox[1]*binsizeinv > MAXSMALLINT ||
bbox[2]*binsizeinv > MAXSMALLINT)
error->all(FLERR,"Domain too large for neighbor bins");
// create actual bins
// always have one bin even if cutoff > bbox
// for 2d, nbinz = 1
nbinx = static_cast<int> (bbox[0]*binsizeinv);
nbiny = static_cast<int> (bbox[1]*binsizeinv);
if (dimension == 3) nbinz = static_cast<int> (bbox[2]*binsizeinv);
else nbinz = 1;
if (nbinx == 0) nbinx = 1;
if (nbiny == 0) nbiny = 1;
if (nbinz == 0) nbinz = 1;
// compute actual bin size for nbins to fit into box exactly
// error if actual bin size << cutoff, since will create a zillion bins
// this happens when nbin = 1 and box size << cutoff
// typically due to non-periodic, flat system in a particular dim
// in that extreme case, should use NSQ not BIN neighbor style
binsizex = bbox[0]/nbinx;
binsizey = bbox[1]/nbiny;
binsizez = bbox[2]/nbinz;
bininvx = 1.0 / binsizex;
bininvy = 1.0 / binsizey;
bininvz = 1.0 / binsizez;
if (binsize_optimal*bininvx > CUT2BIN_RATIO ||
binsize_optimal*bininvy > CUT2BIN_RATIO ||
binsize_optimal*bininvz > CUT2BIN_RATIO)
error->all(FLERR,"Cannot use neighbor bins - box size << cutoff");
// mbinlo/hi = lowest and highest global bins my ghost atoms could be in
// coord = lowest and highest values of coords for my ghost atoms
// static_cast(-1.5) = -1, so subract additional -1
// add in SMALL for round-off safety
int mbinxhi,mbinyhi,mbinzhi;
double coord;
coord = bsubboxlo[0] - SMALL*bbox[0];
mbinxlo = static_cast<int> ((coord-bboxlo[0])*bininvx);
if (coord < bboxlo[0]) mbinxlo = mbinxlo - 1;
coord = bsubboxhi[0] + SMALL*bbox[0];
mbinxhi = static_cast<int> ((coord-bboxlo[0])*bininvx);
coord = bsubboxlo[1] - SMALL*bbox[1];
mbinylo = static_cast<int> ((coord-bboxlo[1])*bininvy);
if (coord < bboxlo[1]) mbinylo = mbinylo - 1;
coord = bsubboxhi[1] + SMALL*bbox[1];
mbinyhi = static_cast<int> ((coord-bboxlo[1])*bininvy);
if (dimension == 3) {
coord = bsubboxlo[2] - SMALL*bbox[2];
mbinzlo = static_cast<int> ((coord-bboxlo[2])*bininvz);
if (coord < bboxlo[2]) mbinzlo = mbinzlo - 1;
coord = bsubboxhi[2] + SMALL*bbox[2];
mbinzhi = static_cast<int> ((coord-bboxlo[2])*bininvz);
}
// extend bins by 1 to insure stencil extent is included
// if 2d, only 1 bin in z
mbinxlo = mbinxlo - 1;
mbinxhi = mbinxhi + 1;
mbinx = mbinxhi - mbinxlo + 1;
mbinylo = mbinylo - 1;
mbinyhi = mbinyhi + 1;
mbiny = mbinyhi - mbinylo + 1;
if (dimension == 3) {
mbinzlo = mbinzlo - 1;
mbinzhi = mbinzhi + 1;
} else mbinzlo = mbinzhi = 0;
mbinz = mbinzhi - mbinzlo + 1;
// memory for bin ptrs
bigint bbin = ((bigint) mbinx) * ((bigint) mbiny) * ((bigint) mbinz);
if (bbin > MAXSMALLINT) error->one(FLERR,"Too many neighbor bins");
mbins = bbin;
if (mbins > maxhead) {
maxhead = mbins;
memory->destroy(binhead);
// USER-INTEL package requires one additional element
#if defined(LMP_USER_INTEL)
memory->create(binhead,maxhead + 1,"neigh:binhead");
#else
memory->create(binhead,maxhead,"neigh:binhead");
#endif
}
// create stencil of bins to search over in neighbor list construction
// sx,sy,sz = max range of stencil in each dim
// smax = max possible size of entire 3d stencil
// stencil is empty if cutneighmax = 0.0
sx = static_cast<int> (cutneighmax*bininvx);
if (sx*binsizex < cutneighmax) sx++;
sy = static_cast<int> (cutneighmax*bininvy);
if (sy*binsizey < cutneighmax) sy++;
sz = static_cast<int> (cutneighmax*bininvz);
if (sz*binsizez < cutneighmax) sz++;
if (dimension == 2) sz = 0;
smax = (2*sx+1) * (2*sy+1) * (2*sz+1);
// create stencils for pairwise neighbor lists
// only done for lists with stencilflag and buildflag set
for (int i = 0; i < nslist; i++) {
if (lists[slist[i]]) {
lists[slist[i]]->stencil_allocate(smax,style);
(this->*stencil_create[slist[i]])(lists[slist[i]],sx,sy,sz);
} else setup_bins_kokkos(i);
}
}
/* ----------------------------------------------------------------------
compute closest distance between central bin (0,0,0) and bin (i,j,k)
------------------------------------------------------------------------- */
double Neighbor::bin_distance(int i, int j, int k)
{
double delx,dely,delz;
if (i > 0) delx = (i-1)*binsizex;
else if (i == 0) delx = 0.0;
else delx = (i+1)*binsizex;
if (j > 0) dely = (j-1)*binsizey;
else if (j == 0) dely = 0.0;
else dely = (j+1)*binsizey;
if (k > 0) delz = (k-1)*binsizez;
else if (k == 0) delz = 0.0;
else delz = (k+1)*binsizez;
return (delx*delx + dely*dely + delz*delz);
}
/* ----------------------------------------------------------------------
set neighbor style and skin distance
------------------------------------------------------------------------- */
void Neighbor::set(int narg, char **arg)
{
if (narg != 2) error->all(FLERR,"Illegal neighbor command");
skin = force->numeric(FLERR,arg[0]);
if (skin < 0.0) error->all(FLERR,"Illegal neighbor command");
if (strcmp(arg[1],"nsq") == 0) style = NSQ;
else if (strcmp(arg[1],"bin") == 0) style = BIN;
else if (strcmp(arg[1],"multi") == 0) style = MULTI;
else error->all(FLERR,"Illegal neighbor command");
if (style == MULTI && lmp->citeme) lmp->citeme->add(cite_neigh_multi);
}
/* ----------------------------------------------------------------------
reset timestamps in all NeighList classes
so that neighbor lists will rebuild properly with timestep change
------------------------------------------------------------------------- */
void Neighbor::reset_timestep(bigint ntimestep)
{
for (int i = 0; i < nlist; i++)
lists[i]->last_build = -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 = force->inumeric(FLERR,arg[iarg+1]);
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 = force->inumeric(FLERR,arg[iarg+1]);
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");
if (strcmp(arg[iarg+1],"yes") == 0) dist_check = 1;
else if (strcmp(arg[iarg+1],"no") == 0) dist_check = 0;
else error->all(FLERR,"Illegal neigh_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"once") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) build_once = 1;
else if (strcmp(arg[iarg+1],"no") == 0) build_once = 0;
else error->all(FLERR,"Illegal neigh_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"page") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
old_pgsize = pgsize;
pgsize = force->inumeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"one") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
old_oneatom = oneatom;
oneatom = force->inumeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"binsize") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
binsize_user = force->numeric(FLERR,arg[iarg+1]);
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");
if (strcmp(arg[iarg+1],"yes") == 0) cluster_check = 1;
else if (strcmp(arg[iarg+1],"no") == 0) cluster_check = 0;
else error->all(FLERR,"Illegal neigh_modify command");
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 == NULL ||
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] = force->inumeric(FLERR,arg[iarg+2]);
ex2_type[nex_type] = force->inumeric(FLERR,arg[iarg+3]);
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") == 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");
}
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");
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 error->all(FLERR,"Illegal neigh_modify command");
}
}
/* ----------------------------------------------------------------------
bin owned and ghost atoms
------------------------------------------------------------------------- */
void Neighbor::bin_atoms()
{
int i,ibin;
for (i = 0; i < mbins; i++) binhead[i] = -1;
// bin in reverse order so linked list will be in forward order
// also puts ghost atoms at end of list, which is necessary
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
if (includegroup) {
int bitmask = group->bitmask[includegroup];
for (i = nall-1; i >= nlocal; i--) {
if (mask[i] & bitmask) {
ibin = coord2bin(x[i]);
bins[i] = binhead[ibin];
binhead[ibin] = i;
}
}
for (i = atom->nfirst-1; i >= 0; i--) {
ibin = coord2bin(x[i]);
bins[i] = binhead[ibin];
binhead[ibin] = i;
}
} else {
for (i = nall-1; i >= 0; i--) {
ibin = coord2bin(x[i]);
bins[i] = binhead[ibin];
binhead[ibin] = i;
}
}
}
/* ----------------------------------------------------------------------
convert atom coords into local bin #
for orthogonal, only ghost atoms will have coord >= bboxhi or coord < bboxlo
take special care to insure ghosts are in correct bins even w/ roundoff
hi ghost atoms = nbin,nbin+1,etc
owned atoms = 0 to nbin-1
lo ghost atoms = -1,-2,etc
this is necessary so that both procs on either side of PBC
treat a pair of atoms straddling the PBC in a consistent way
for triclinic, doesn't matter since stencil & neigh list built differently
------------------------------------------------------------------------- */
int Neighbor::coord2bin(double *x)
{
int ix,iy,iz;
if (!ISFINITE(x[0]) || !ISFINITE(x[1]) || !ISFINITE(x[2]))
error->one(FLERR,"Non-numeric positions - simulation unstable");
if (x[0] >= bboxhi[0])
ix = static_cast<int> ((x[0]-bboxhi[0])*bininvx) + nbinx;
else if (x[0] >= bboxlo[0]) {
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx);
ix = MIN(ix,nbinx-1);
} else
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx) - 1;
if (x[1] >= bboxhi[1])
iy = static_cast<int> ((x[1]-bboxhi[1])*bininvy) + nbiny;
else if (x[1] >= bboxlo[1]) {
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy);
iy = MIN(iy,nbiny-1);
} else
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy) - 1;
if (x[2] >= bboxhi[2])
iz = static_cast<int> ((x[2]-bboxhi[2])*bininvz) + nbinz;
else if (x[2] >= bboxlo[2]) {
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz);
iz = MIN(iz,nbinz-1);
} else
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz) - 1;
return (iz-mbinzlo)*mbiny*mbinx + (iy-mbinylo)*mbinx + (ix-mbinxlo);
}
/* ----------------------------------------------------------------------
same as coord2bin, but also return ix,iy,iz offsets in each dim
------------------------------------------------------------------------- */
int Neighbor::coord2bin(double *x, int &ix, int &iy, int &iz)
{
if (!ISFINITE(x[0]) || !ISFINITE(x[1]) || !ISFINITE(x[2]))
error->one(FLERR,"Non-numeric positions - simulation unstable");
if (x[0] >= bboxhi[0])
ix = static_cast<int> ((x[0]-bboxhi[0])*bininvx) + nbinx;
else if (x[0] >= bboxlo[0]) {
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx);
ix = MIN(ix,nbinx-1);
} else
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx) - 1;
if (x[1] >= bboxhi[1])
iy = static_cast<int> ((x[1]-bboxhi[1])*bininvy) + nbiny;
else if (x[1] >= bboxlo[1]) {
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy);
iy = MIN(iy,nbiny-1);
} else
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy) - 1;
if (x[2] >= bboxhi[2])
iz = static_cast<int> ((x[2]-bboxhi[2])*bininvz) + nbinz;
else if (x[2] >= bboxlo[2]) {
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz);
iz = MIN(iz,nbinz-1);
} else
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz) - 1;
ix -= mbinxlo;
iy -= mbinylo;
iz -= mbinzlo;
return iz*mbiny*mbinx + iy*mbinx + ix;
}
/* ----------------------------------------------------------------------
test if atom pair i,j is excluded from neighbor list
due to type, group, molecule settings from neigh_modify command
return 1 if should be excluded, 0 if included
------------------------------------------------------------------------- */
int Neighbor::exclusion(int i, int j, int itype, int jtype,
int *mask, tagint *molecule) const {
int m;
if (nex_type && ex_type[itype][jtype]) return 1;
if (nex_group) {
for (m = 0; m < nex_group; m++) {
if (mask[i] & ex1_bit[m] && mask[j] & ex2_bit[m]) return 1;
if (mask[i] & ex2_bit[m] && mask[j] & ex1_bit[m]) return 1;
}
}
if (nex_mol) {
for (m = 0; m < nex_mol; m++)
if (mask[i] & ex_mol_bit[m] && mask[j] & ex_mol_bit[m] &&
molecule[i] == molecule[j]) return 1;
}
return 0;
}
/* ----------------------------------------------------------------------
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 # of bytes of allocated memory
------------------------------------------------------------------------- */
bigint Neighbor::memory_usage()
{
bigint bytes = 0;
bytes += memory->usage(xhold,maxhold,3);
if (style != NSQ) {
bytes += memory->usage(bins,maxbin);
bytes += memory->usage(binhead,maxhead);
}
for (int i = 0; i < nrequest; i++)
if (lists[i]) bytes += lists[i]->memory_usage();
bytes += memory->usage(bondlist,maxbond,3);
bytes += memory->usage(anglelist,maxangle,4);
bytes += memory->usage(dihedrallist,maxdihedral,5);
bytes += memory->usage(improperlist,maximproper,5);
return bytes;
}
/* ----------------------------------------------------------------------
return the value of exclude - used to check compatibility with GPU
------------------------------------------------------------------------- */
int Neighbor::exclude_setting()
{
return exclude;
}
c4science · Help