diff --git a/doc/src/Eqs/fix_mvv_dpd.jpg b/doc/src/Eqs/fix_mvv_dpd.jpg index f62ae28bc..b5f39bb47 100644 Binary files a/doc/src/Eqs/fix_mvv_dpd.jpg and b/doc/src/Eqs/fix_mvv_dpd.jpg differ diff --git a/doc/src/Eqs/fix_mvv_dpd.tex b/doc/src/Eqs/fix_mvv_dpd.tex index 4652d54b7..685d52776 100644 --- a/doc/src/Eqs/fix_mvv_dpd.tex +++ b/doc/src/Eqs/fix_mvv_dpd.tex @@ -1,21 +1,21 @@ \documentclass[12pt]{article} \begin{document} $$ v(t+\frac{\Delta t}{2}) = v(t) + \frac{\Delta t}{2}\cdot a(t), $$ $$ r(t+\Delta t) = r(t) + \Delta t\cdot v(t+\frac{\Delta t}{2}), $$ $$ a(t+\Delta t) = \frac{1}{m}\cdot F\left[ r(t+\Delta t), v(t) +\lambda \cdot \Delta t\cdot a(t)\right], $$ $$ - v(t+\Delta t) = v(t+\frac{\Delta t}{2}) + \frac{\Delta t}{2}\cdot a(t++\Delta t), + v(t+\Delta t) = v(t+\frac{\Delta t}{2}) + \frac{\Delta t}{2}\cdot a(t+\Delta t) $$ \end{document} diff --git a/src/neighbor.cpp b/src/neighbor.cpp index ef150902e..3943d876f 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -1,2462 +1,2466 @@ /* ---------------------------------------------------------------------- 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 #include #include #include #include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" +#include "nbin.h" +#include "nstencil.h" +#include "npair.h" +#include "ntopo.h" #include "style_nbin.h" #include "style_nstencil.h" #include "style_npair.h" #include "style_ntopo.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" #include using namespace LAMMPS_NS; using namespace NeighConst; #define RQDELTA 1 #define EXDELTA 1 #define BIG 1.0e20 enum{NSQ,BIN,MULTI}; // also in NBin, NeighList, NStencil enum{NONE,ALL,PARTIAL,TEMPLATE}; 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), pairclass(NULL), pairnames(NULL), pairmasks(NULL) { MPI_Comm_rank(world,&me); MPI_Comm_size(world,&nprocs); firsttime = 1; style = 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 = NULL; cutneighghostsq = NULL; cuttype = NULL; cuttypesq = NULL; fixchecklist = NULL; // pairwise neighbor lists and associated data structs nlist = 0; lists = NULL; nbin = 0; neigh_bin = NULL; nstencil = 0; neigh_stencil = NULL; neigh_pair = NULL; nstencil_perpetual = 0; slist = NULL; npair_perpetual = 0; plist = NULL; nrequest = maxrequest = 0; requests = NULL; old_nrequest = 0; old_requests = NULL; old_style = style; old_triclinic = 0; old_pgsize = pgsize; old_oneatom = oneatom; zeroes = NULL; binclass = NULL; binnames = NULL; binmasks = NULL; stencilclass = NULL; stencilnames = NULL; stencilmasks = NULL; // topology lists bondwhich = anglewhich = dihedralwhich = improperwhich = NONE; neigh_bond = NULL; neigh_angle = NULL; neigh_dihedral = NULL; neigh_improper = NULL; // coords at last neighboring maxhold = 0; xhold = NULL; lastcall = -1; last_setup_bins = -1; // 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 = ex_mol_intra = NULL; // 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 [] zeroes; 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); } /* ---------------------------------------------------------------------- */ 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 = 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; // 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; } // fixchecklist = 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; // ------------------------------------------------------------------ // xhold array // free if not needed for this run if (dist_check == 0) { memory->destroy(xhold); maxhold = 0; xhold = NULL; } // 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"); // ------------------------------------------------------------------ // 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 (int i = 0; i < nrequest; i++) { delete requests[i]; requests[i] = NULL; } 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" #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] = &bin_creator; \ binmasks[nbclass++] = bitmasks; #include "style_nbin.h" #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" #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] = &stencil_creator; \ stencilmasks[nsclass++] = bitmasks; #include "style_nstencil.h" #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" #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] = &pair_creator; \ pairmasks[npclass++] = bitmasks; #include "style_npair.h" #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 == 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) other = point history and rRESPA lists at their partner lists // (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: // (1) before (2), b/c (2) needs to know history partner pairings // (2) after (1), b/c (2) may also need to create new history lists // (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_other(); 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; 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; for (j = 0; j < nbin; j++) if (neigh_bin[j]->istyle == flag) break; if (j < nbin && !requests[i]->unique) { 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; for (j = 0; j < nstencil; j++) if (neigh_stencil[j]->istyle == flag) break; if (j < nstencil && !requests[i]->unique) { 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 == NULL) 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] = NULL; continue; } PairCreator pair_creator = pairclass[flag-1]; 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 == NULL) 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 == NULL) 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 // allocate dnum vector of zeroes if set int dnummax = 0; for (i = 0; i < nlist; i++) { if (lists[i]->copy) continue; lists[i]->setup_pages(pgsize,oneatom); dnummax = MAX(dnummax,lists[i]->dnum); } if (dnummax) { delete [] zeroes; zeroes = new double[dnummax]; for (i = 0; i < dnummax; i++) zeroes[i] = 0.0; } // first-time allocation of per-atom data for lists that are built and store // lists that are not built: granhistory, respa inner/middle (no neigh_pair) // 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]->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 = NULL; 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 == NULL) 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 only for history, respaouter, custom cutoff lists ------------------------------------------------------------------------- */ void Neighbor::morph_other() { NeighRequest *irq; for (int i = 0; i < nrequest; i++) { irq = requests[i]; // if history, point this list and partner list at each other if (irq->history) { irq->historylist = i-1; requests[i-1]->history_partner = 1; requests[i-1]->historylist = i; } // if respaouter, point all associated rRESPA lists at each other if (irq->respaouter) { if (requests[i-1]->respainner) { irq->respainnerlist = i-1; requests[i-1]->respaouterlist = i; } else { irq->respamiddlelist = i-1; requests[i-1]->respaouterlist = i; requests[i-1]->respainnerlist = i-1; irq->respainnerlist = i-2; requests[i-2]->respaouterlist = i; requests[i-2]->respamiddlelist = i-1; } } // 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; } } /* ---------------------------------------------------------------------- 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->history) continue; if (irq->respainner || irq->respamiddle) continue; 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 // no need to check respaouter b/c it stores same pairs // no need to check dnum b/c only set for history // NOTE: need check for 2 Kokkos flags? if (irq->ghost != jrq->ghost) continue; if (irq->size != jrq->size) 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 list // b/c child skip lists store their own history info 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; // Kokkos doesn't yet support half from full if (irq->kokkos_host) continue; if (irq->kokkos_device) continue; // these lists are created other ways, no need for halffull // do want to process skip lists if (irq->history) continue; if (irq->respainner || irq->respamiddle) continue; 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 // no need to check respaouter b/c it stores same pairs // no need to check dnum b/c only set for history if (irq->ghost != jrq->ghost) continue; if (irq->size != jrq->size) 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; // these lists are created other ways, no need to copy // skip lists are eligible to become a copy list if (irq->history) continue; if (irq->respainner || irq->respamiddle) 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 history list or a respa middle/inner list if (jrq->history) continue; 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 respaouter b/c it stores same pairs // no need to check omp b/c it stores same pairs // no need to check dnum b/c only set for history // NOTE: need check for 2 Kokkos flags? if (irq->size != jrq->size) 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) 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 = 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 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 == 2) 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 == 2) 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 == 2) 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 == 2) 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,m; char str[128]; NeighRequest *rq; FILE *out; 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]; 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++; } for (m = 0; m < 2; m++) { if (m == 0) out = screen; else out = logfile; if (out) { fprintf(out,"Neighbor list info ...\n"); fprintf(out," update every %d steps, delay %d steps, check %s\n", every,delay,dist_check ? "yes" : "no"); fprintf(out," max neighbors/atom: %d, page size: %d\n", oneatom, pgsize); fprintf(out," master list distance cutoff = %g\n",cutneighmax); fprintf(out," ghost atom cutoff = %g\n",cutghost); if (style != NSQ) fprintf(out," binsize = %g, bins = %g %g %g\n",binsize, ceil(bbox[0]/binsize), ceil(bbox[1]/binsize), ceil(bbox[2]/binsize)); fprintf(out," %d neighbor lists, " "perpetual/occasional/extra = %d %d %d\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); sprintf(str," (%d) pair %s",i+1,pname); } else if (rq->fix) { sprintf(str," (%d) fix %s",i+1,((Fix *) rq->requestor)->style); } else if (rq->compute) { sprintf(str," (%d) compute %s",i+1, ((Compute *) rq->requestor)->style); } else if (rq->command) { sprintf(str," (%d) command %s",i+1,rq->command_style); } else if (rq->neigh) { sprintf(str," (%d) neighbor class addition",i+1); } fprintf(out,"%s",str); if (rq->occasional) fprintf(out,", occasional"); else fprintf(out,", perpetual"); // order these to get single output of most relevant if (rq->history) fprintf(out,", history for (%d)",rq->historylist+1); else if (rq->copy) fprintf(out,", copy from (%d)",rq->copylist+1); else if (rq->halffull) fprintf(out,", half/full from (%d)",rq->halffulllist+1); else if (rq->skip) fprintf(out,", skip from (%d)",rq->skiplist+1); fprintf(out,"\n"); // list of neigh list attributes fprintf(out," attributes: "); if (rq->half) fprintf(out,"half"); else if (rq->full) fprintf(out,"full"); if (rq->newton == 0) { if (force->newton_pair) fprintf(out,", newton on"); else fprintf(out,", newton off"); } else if (rq->newton == 1) fprintf(out,", newton on"); else if (rq->newton == 2) fprintf(out,", newton off"); if (rq->ghost) fprintf(out,", ghost"); if (rq->size) fprintf(out,", size"); if (rq->history) fprintf(out,", history"); if (rq->granonesided) fprintf(out,", onesided"); if (rq->respainner) fprintf(out,", respa outer"); if (rq->respamiddle) fprintf(out,", respa middle"); if (rq->respaouter) fprintf(out,", respa inner"); if (rq->bond) fprintf(out,", bond"); if (rq->omp) fprintf(out,", omp"); if (rq->intel) fprintf(out,", intel"); if (rq->kokkos_device) fprintf(out,", kokkos_device"); if (rq->kokkos_host) fprintf(out,", kokkos_host"); if (rq->ssa) fprintf(out,", ssa"); if (rq->cut) fprintf(out,", cut %g",rq->cutoff); if (rq->off2on) fprintf(out,", off2on"); fprintf(out,"\n"); fprintf(out," "); if (lists[i]->pair_method == 0) fprintf(out,"pair build: none\n"); else fprintf(out,"pair build: %s\n",pairnames[lists[i]->pair_method-1]); fprintf(out," "); if (lists[i]->stencil_method == 0) fprintf(out,"stencil: none\n"); else fprintf(out,"stencil: %s\n", stencilnames[lists[i]->stencil_method-1]); fprintf(out," "); if (lists[i]->bin_method == 0) fprintf(out,"bin: none\n"); else fprintf(out,"bin: %s\n",binnames[lists[i]->bin_method-1]); } /* fprintf(out," %d stencil methods\n",nstencil); for (i = 0; i < nstencil; i++) fprintf(out," (%d) %s\n", i+1,stencilnames[neigh_stencil[i]->istyle-1]); fprintf(out," %d bin methods\n",nbin); for (i = 0; i < nbin; i++) fprintf(out," (%d) %s\n",i+1,binnames[neigh_bin[i]->istyle-1]); */ } } } /* ---------------------------------------------------------------------- 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(lmp); old_requests[i]->copy_request(requests[i],1); } old_style = style; old_triclinic = triclinic; old_pgsize = pgsize; old_oneatom = oneatom; } /* ---------------------------------------------------------------------- find and return request made by classptr if not found or classpt = NULL, return NULL ------------------------------------------------------------------------- */ NeighRequest *Neighbor::find_request(void *classptr) { if (classptr == NULL) return NULL; for (int i = 0; i < nrequest; i++) if (requests[i]->requestor == classptr) return requests[i]; return NULL; } /* ---------------------------------------------------------------------- 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 == NSQ) return 0; if (rq->skip || rq->copy || rq->halffull) return 0; if (rq->history) return 0; if (rq->respainner || rq->respamiddle) 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; 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 == NSQ) return 0; if (rq->skip || rq->copy || rq->halffull) return 0; if (rq->history) return 0; if (rq->respainner || rq->respamiddle) return 0; // convert newton request to newtflag = on or off int newtflag; 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; //printf("STENCIL RQ FLAGS: hff %d %d n %d g %d s %d newtflag %d\n", // rq->half,rq->full,rq->newton,rq->ghost,rq->ssa, // newtflag); // 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 newton %d newtoff %d ghost %d ssa %d\n", // i,mask & NS_HALF,mask & NS_FULL,mask & NS_NEWTON, // mask & NS_NEWTOFF,mask & NS_GHOST,mask & NS_SSA); // exactly one of half or full is set and must match if (rq->half) { if (!(mask & NS_HALF)) continue; } else if (rq->full) { if (!(mask & NS_FULL)) continue; } // newtflag is on or off and must match if (newtflag) { if (!(mask & NS_NEWTON)) continue; } else if (!newtflag) { if (!(mask & NS_NEWTOFF)) 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 BIN or MULTI and must match if (style == BIN) { if (!(mask & NS_BIN)) continue; } else if (style == 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) { // no neighbor list build performed if (rq->history) return 0; if (rq->respainner || rq->respamiddle) return 0; // 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 int newtflag; 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; 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 != !(mask & NP_KOKKOS_DEVICE)) continue; if (!rq->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) { if (mask & NP_ATOMONLY) continue; } else if (!molecular) { 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 and must match if (style == NSQ) { if (!(mask & NP_NSQ)) continue; } else if (style == BIN) { if (!(mask & NP_BIN)) continue; } else if (style == 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); requests[nrequest]->index = nrequest; requests[nrequest]->requestor = requestor; requests[nrequest]->requestor_instance = instance; nrequest++; return nrequest-1; } /* ---------------------------------------------------------------------- one instance per entry in style_neigh_bin.h ------------------------------------------------------------------------- */ template NBin *Neighbor::bin_creator(LAMMPS *lmp) { return new T(lmp); } /* ---------------------------------------------------------------------- one instance per entry in style_neigh_stencil.h ------------------------------------------------------------------------- */ template NStencil *Neighbor::stencil_creator(LAMMPS *lmp) { return new T(lmp); } /* ---------------------------------------------------------------------- one instance per entry in style_neigh_pair.h ------------------------------------------------------------------------- */ template NPair *Neighbor::pair_creator(LAMMPS *lmp) { return new T(lmp); } /* ---------------------------------------------------------------------- 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 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 (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; // 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 != NSQ) { for (int 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]->grow(nlocal,nall); neigh_pair[m]->build_setup(); neigh_pair[m]->build(lists[m]); } // build topology lists for bonds/angles/etc if (atom->molecular && 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 == NULL) error->all(FLERR,"Trying to build an occasional neighbor list " "before initialization completed"); // 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, // 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); // 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->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 = 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 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 = 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/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 error->all(FLERR,"Illegal neigh_modify command"); } } /* ---------------------------------------------------------------------- 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; } /* ---------------------------------------------------------------------- return # of bytes of allocated memory ------------------------------------------------------------------------- */ bigint Neighbor::memory_usage() { bigint 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; }