neighbor.cpp
7dfc6b7eab71lammm-devel

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 "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 <map>

			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;

			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 [] 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.25skinskin;
			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<Class>; \
			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<Class>; \
			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<Class>; \
			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) unique = create unique lists if cutoff is explicitly set
			// (2) skip = create any new non-skip lists needed by pair hybrid skip lists
			// (3) granular = adjust parent and skip lists for granular onesided usage
			// (4) h/f = pair up any matching half/full lists
			// (5) copy = convert as many lists as possible to copy lists
			// order of morph methods matters:
			// (3) after (2), b/c it adjusts lists created by (2)
			// (4) after (2) and (3),
			// b/c (2) may create new full lists, (3) may change them
			// (5) last, after all lists are finalized, so all possible copies found

			int nrequest_original = nrequest;

			morph_unique();
			morph_skip();
			morph_granular(); // this method can change flags set by requestor
			morph_halffull();
			morph_copy();

			// create new lists, one per request including added requests
			// wait to allocate initial pages until copy lists are detected
			// NOTE: can I allocate now, instead of down below?

			nlist = nrequest;

			lists = new NeighList*[nrequest];
			neigh_bin = new NBin*[nrequest];
			neigh_stencil = new NStencil*[nrequest];
			neigh_pair = new NPair*[nrequest];

			// allocate new lists
			// pass list ptr back to requestor (except for Command class)
			// only for original requests, not ones added by Neighbor class

			for (i = 0; i < nrequest; i++) {
			if (requests[i]->kokkos_host \|\| requests[i]->kokkos_device)
			create_kokkos_list(i);
			else lists[i] = new NeighList(lmp);
			lists[i]->index = i;

			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];
			lists[i]->np = neigh_pair[i] = pair_creator(lmp);
			neigh_pair[i]->post_constructor(requests[i]);
			neigh_pair[i]->istyle = flag;

			if (lists[i]->bin_method > 0) {
			neigh_pair[i]->nb = neigh_bin[requests[i]->index_bin];
			if (neigh_pair[i]->nb == 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

			for (i = 0; i < nlist; i++) {
			if (lists[i]->copy) continue;
			lists[i]->setup_pages(pgsize,oneatom);
			}

			// first-time allocation of per-atom data for lists that are built and store
			// lists that do not store: copy
			// use atom->nmax for both grow() args
			// i.e. grow first time to expanded size to avoid future reallocs
			// also Kokkos list initialization

			int maxatom = atom->nmax;
			for (i = 0; i < nlist; i++)
			if (neigh_pair[i] && !lists[i]->copy) lists[i]->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 custom cutoff lists
			------------------------------------------------------------------------- */

			void Neighbor::morph_unique()
			{
			NeighRequest *irq;

			for (int i = 0; i < nrequest; i++) {
			irq = requests[i];

			// if cut flag set by requestor, set unique flag
			// this forces Pair,Stencil,Bin styles to be instantiated separately

			if (irq->cut) irq->unique = 1;
			}
			}

			/* ----------------------------------------------------------------------
			scan NeighRequests to process all skip lists
			look for a matching non-skip list
			if one exists, point at it via skiplist
			else make new parent via copy_request() and point at it
			------------------------------------------------------------------------- */

			void Neighbor::morph_skip()
			{
			int i,j,inewton,jnewton;
			NeighRequest irq,jrq,*nrq;

			for (i = 0; i < nrequest; i++) {
			irq = requests[i];

			// only processing skip lists

			if (!irq->skip) continue;

			// these lists are created other ways, no need for skipping
			// halffull list and its full parent may both skip,
			// but are checked to insure matching skip info

			if (irq->halffull) continue;
			if (irq->copy) continue;

			// check all other lists

			for (j = 0; j < nrequest; j++) {
			if (i == j) continue;
			jrq = requests[j];

			// can only skip from a perpetual non-skip list

			if (jrq->occasional) continue;
			if (jrq->skip) continue;

			// both lists must be half, or both full

			if (irq->half != jrq->half) continue;
			if (irq->full != jrq->full) continue;

			// both lists must be newton on, or both newton off
			// IJ newton = 1 for newton on, 2 for newton off

			inewton = irq->newton;
			if (inewton == 0) inewton = force->newton_pair ? 1 : 2;
			jnewton = jrq->newton;
			if (jnewton == 0) jnewton = force->newton_pair ? 1 : 2;
			if (inewton != jnewton) continue;

			// these flags must be same,
			// else 2 lists do not store same pairs
			// or their data structures are different
			// this includes custom cutoff set by requestor
			// NOTE: need check for 2 Kokkos flags?

			if (irq->ghost != jrq->ghost) continue;
			if (irq->size != jrq->size) continue;
			if (irq->history != jrq->history) continue;
			if (irq->bond != jrq->bond) continue;
			if (irq->omp != jrq->omp) continue;
			if (irq->intel != jrq->intel) continue;
			if (irq->kokkos_host != jrq->kokkos_host) continue;
			if (irq->kokkos_device != jrq->kokkos_device) continue;
			if (irq->ssa != jrq->ssa) continue;
			if (irq->cut != jrq->cut) continue;
			if (irq->cutoff != jrq->cutoff) continue;

			// 2 lists are a match

			break;
			}

			// if matching list exists, point to it
			// else create a new identical list except non-skip
			// for new list, set neigh = 1, skip = 0, no skip vec/array,
			// copy unique flag (since copy_request() will not do it)
			// note: parents of skip lists do not have associated history
			// b/c child skip lists have the associated history

			if (j < nrequest) irq->skiplist = j;
			else {
			int newrequest = request(this,-1);
			irq->skiplist = newrequest;

			nrq = requests[newrequest];
			nrq->copy_request(irq,0);
			nrq->pair = nrq->fix = nrq->compute = nrq->command = 0;
			nrq->neigh = 1;
			nrq->skip = 0;
			if (irq->unique) nrq->unique = 1;
			}
			}
			}

			/* ----------------------------------------------------------------------
			scan NeighRequests just added by morph_skip for hybrid granular
			adjust newton/oneside parent settings if children require onesided skipping
			also set children off2on flag if parent becomes a newton off list
			this is needed because line/gran and tri/gran pair styles
			require onesided neigh lists and system newton on,
			but parent list must be newton off to enable the onesided skipping
			------------------------------------------------------------------------- */

			void Neighbor::morph_granular()
			{
			int i,j;
			NeighRequest irq,jrq;

			for (i = 0; i < nrequest; i++) {
			irq = requests[i];

			// only examine NeighRequests added by morph_skip()
			// only those with size attribute for granular systems

			if (!irq->neigh) continue;
			if (!irq->size) continue;

			// check children of this list

			int onesided = -1;
			for (j = 0; j < nrequest; j++) {
			jrq = requests[j];

			// only consider JRQ pair, size lists that skip from Irq list

			if (!jrq->pair) continue;
			if (!jrq->size) continue;
			if (!jrq->skip \|\| jrq->skiplist != i) continue;

			// onesided = -1 if no children
			// onesided = 0/1 = child granonesided value if same for all children
			// onesided = 2 if children have different granonesided values

			if (onesided < 0) onesided = jrq->granonesided;
			else if (onesided != jrq->granonesided) onesided = 2;
			if (onesided == 2) break;
			}

			// if onesided = 2, parent has children with both granonesided = 0/1
			// force parent newton off (newton = 2) to enable onesided skip by child
			// set parent granonesided = 0, so it stores all neighs in usual manner
			// set off2on = 1 for all children, since they expect newton on lists
			// this is b/c granonesided only set by line/gran and tri/gran which
			// both require system newton on

			if (onesided == 2) {
			irq->newton = 2;
			irq->granonesided = 0;

			for (j = 0; j < nrequest; j++) {
			jrq = requests[j];

			// only consider JRQ pair, size lists that skip from Irq list

			if (!jrq->pair) continue;
			if (!jrq->size) continue;
			if (!jrq->skip \|\| jrq->skiplist != i) continue;

			jrq->off2on = 1;
			}
			}
			}
			}

			/* ----------------------------------------------------------------------
			scan NeighRequests for possible half lists to derive from full lists
			if 2 requests match, set half list to derive from full list
			------------------------------------------------------------------------- */

			void Neighbor::morph_halffull()
			{
			int i,j;
			NeighRequest irq,jrq;

			for (i = 0; i < nrequest; i++) {
			irq = requests[i];

			// only processing half lists

			if (!irq->half) continue;

			// 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->copy) continue;

			// check all other lists

			for (j = 0; j < nrequest; j++) {
			if (i == j) continue;
			jrq = requests[j];

			// can only derive from a perpetual full list
			// newton setting of derived list does not matter

			if (jrq->occasional) continue;
			if (!jrq->full) continue;

			// these flags must be same,
			// else 2 lists do not store same pairs
			// or their data structures are different
			// this includes custom cutoff set by requestor

			if (irq->ghost != jrq->ghost) continue;
			if (irq->size != jrq->size) continue;
			if (irq->history != jrq->history) continue;
			if (irq->bond != jrq->bond) continue;
			if (irq->omp != jrq->omp) continue;
			if (irq->intel != jrq->intel) continue;
			if (irq->kokkos_host != jrq->kokkos_host) continue;
			if (irq->kokkos_device != jrq->kokkos_device) continue;
			if (irq->ssa != jrq->ssa) continue;
			if (irq->cut != jrq->cut) continue;
			if (irq->cutoff != jrq->cutoff) continue;

			// skip flag must be same
			// if both are skip lists, skip info must match

			if (irq->skip != jrq->skip) continue;
			if (irq->skip && irq->same_skip(jrq) == 0) continue;

			// 2 lists are a match

			break;
			}

			// if matching list exists, point to it

			if (j < nrequest) {
			irq->halffull = 1;
			irq->halffulllist = j;
			}
			}
			}

			/* ----------------------------------------------------------------------
			scan NeighRequests for possible copies
			if 2 requests match, turn one into a copy of the other
			------------------------------------------------------------------------- */

			void Neighbor::morph_copy()
			{
			int i,j,inewton,jnewton;
			NeighRequest irq,jrq;

			for (i = 0; i < nrequest; i++) {
			irq = requests[i];

			// this list is already a copy list due to another morph method

			if (irq->copy) continue;

			// check all other lists

			for (j = 0; j < nrequest; j++) {
			if (i == j) continue;
			jrq = requests[j];

			// other list is already copied from this one

			if (jrq->copy && jrq->copylist == i) continue;

			// other list (jrq) to copy from must be perpetual
			// list that becomes a copy list (irq) can be perpetual or occasional
			// if both lists are perpetual, require j < i
			// to prevent circular dependence with 3 or more copies of a list

			if (jrq->occasional) continue;
			if (!irq->occasional && j > i) continue;

			// both lists must be half, or both full

			if (irq->half != jrq->half) continue;
			if (irq->full != jrq->full) continue;

			// both lists must be newton on, or both newton off
			// IJ newton = 1 for newton on, 2 for newton off

			inewton = irq->newton;
			if (inewton == 0) inewton = force->newton_pair ? 1 : 2;
			jnewton = jrq->newton;
			if (jnewton == 0) jnewton = force->newton_pair ? 1 : 2;
			if (inewton != jnewton) continue;

			// ok for non-ghost list to copy from ghost list, but not vice versa

			if (irq->ghost && !jrq->ghost) continue;

			// do not copy from a list with respa middle/inner
			// b/c its outer list will not be complete

			if (jrq->respamiddle) continue;
			if (jrq->respainner) continue;

			// these flags must be same,
			// else 2 lists do not store same pairs
			// or their data structures are different
			// this includes custom cutoff set by requestor
			// no need to check omp b/c it stores same pairs
			// NOTE: need check for 2 Kokkos flags?

			if (irq->size != jrq->size) continue;
			if (irq->history != jrq->history) continue;
			if (irq->bond != jrq->bond) continue;
			if (irq->intel != jrq->intel) continue;
			if (irq->kokkos_host != jrq->kokkos_host) continue;
			if (irq->kokkos_device != jrq->kokkos_device) continue;
			if (irq->ssa != jrq->ssa) continue;
			if (irq->cut != jrq->cut) continue;
			if (irq->cutoff != jrq->cutoff) continue;

			// skip flag must be same
			// if both are skip lists, skip info must match

			if (irq->skip != jrq->skip) continue;
			if (irq->skip && irq->same_skip(jrq) == 0) continue;

			// 2 lists are a match

			break;
			}

			// turn list I into a copy of list J
			// do not copy a list from another copy list, but from its parent list

			if (j < nrequest) {
			irq->copy = 1;
			if (jrq->copy) irq->copylist = jrq->copylist;
			else irq->copylist = j;
			}
			}
			}

			/* ----------------------------------------------------------------------
			create and initialize NTopo classes
			------------------------------------------------------------------------- */

			void Neighbor::init_topology()
			{
			int i,m;

			if (!atom->molecular) 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->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->respamiddle) fprintf(out,", respa outer/middle/inner");
			else if (rq->respainner) fprintf(out,", respa outer/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_nrequestsizeof(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;

			// 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;

			// 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)
			{
			// 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,maxrequestsizeof(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 <typename T>
			NBin Neighbor::bin_creator(LAMMPS lmp)
			{
			return new T(lmp);
			}

			/* ----------------------------------------------------------------------
			one instance per entry in style_neigh_stencil.h
			------------------------------------------------------------------------- */

			template <typename T>
			NStencil Neighbor::stencil_creator(LAMMPS lmp)
			{
			return new T(lmp);
			}

			/* ----------------------------------------------------------------------
			one instance per entry in style_neigh_pair.h
			------------------------------------------------------------------------- */

			template <typename T>
			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(delxdelx + delydely + delz*delz);
			delx = bboxhi[0] - boxhi_hold[0];
			dely = bboxhi[1] - boxhi_hold[1];
			delz = bboxhi[2] - boxhi_hold[2];
			delta2 = sqrt(delxdelx + delydely + 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(delxdelx + delydely + 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 = delxdelx + delydely + 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;
			}

lammps/src/neighbor.cpp7dfc6b7eab71lammm-devel

neighbor.cppView Options

lammps/src/neighbor.cpp
7dfc6b7eab71lammm-devel

neighbor.cpp
View Options