neighbor.cpp

/* ----------------------------------------------------------------------
   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 "limits.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "atom.h"
#include "atom_vec.h"
#include "comm.h"
#include "force.h"
#include "pair.h"
#include "domain.h"
#include "group.h"
#include "modify.h"
#include "fix.h"
#include "compute.h"
#include "update.h"
#include "respa.h"
#include "output.h"
#include "memory.h"
#include "error.h"

using namespace LAMMPS_NS;

#define RQDELTA 1
#define EXDELTA 1

#define LB_FACTOR 1.5
#define SMALL 1.0e-6
#define BIG 1.0e20
#define CUT2BIN_RATIO 100

#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))

enum{NSQ,BIN,MULTI};     // also in neigh_list.cpp

//#define NEIGH_LIST_DEBUG 1

/* ---------------------------------------------------------------------- */

Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp)
{
  MPI_Comm_rank(world,&me);
  MPI_Comm_size(world,&nprocs);

  style = BIN;
  every = 1;
  delay = 10;
  dist_check = 1;
  pgsize = 100000;
  oneatom = 2000;
  binsizeflag = 0;

  cutneighsq = NULL;
  cuttype = NULL;
  cuttypesq = NULL;
  fixchecklist = NULL;

  // coords at last neighboring

  maxhold = 0;
  xhold = NULL;

  // binning

  maxhead = 0;
  binhead = NULL;
  maxbin = 0;
  bins = NULL;

  // pair exclusion list info

  include_group = 0;

  nex_type = maxex_type = 0;
  ex1_type = ex2_type = NULL;
  ex_type = NULL;

  nex_group = maxex_group = 0;
  ex1_group = ex2_group = ex1_bit = ex2_bit = NULL;

  nex_mol = maxex_mol = 0;
  ex_mol_group = ex_mol_bit = NULL;

  // pair lists

  maxlocal = 0;
  nblist = nglist = nslist = 0;

  nlist = 0;
  lists = NULL;
  pair_build = NULL;
  stencil_create = NULL;
  blist = glist = slist = NULL;

  nrequest = maxrequest = 0;
  requests = NULL;

  old_style = BIN;
  old_nrequest = 0;
  old_requests = NULL;

  // bond lists

  maxbond = 0;
  bondlist = NULL;
  maxangle = 0;
  anglelist = NULL;
  maxdihedral = 0;
  dihedrallist = NULL;
  maximproper = 0;
  improperlist = NULL;
}

/* ---------------------------------------------------------------------- */

Neighbor::~Neighbor()
{
  memory->destroy_2d_double_array(cutneighsq);
  delete [] cuttype;
  delete [] cuttypesq;
  delete [] fixchecklist;

  memory->destroy_2d_double_array(xhold);

  memory->sfree(binhead);
  memory->sfree(bins);

  memory->sfree(ex1_type);
  memory->sfree(ex2_type);
  memory->destroy_2d_int_array(ex_type);

  memory->sfree(ex1_group);
  memory->sfree(ex2_group);
  delete [] ex1_bit;
  delete [] ex2_bit;

  memory->sfree(ex_mol_group);
  delete [] ex_mol_bit;

  for (int i = 0; i < nlist; i++) delete lists[i];
  delete [] lists;
  delete [] pair_build;
  delete [] stencil_create;
  delete [] blist;
  delete [] glist;
  delete [] slist;

  for (int i = 0; i < nrequest; i++) delete requests[i];
  memory->sfree(requests);
  for (int i = 0; i < old_nrequest; i++) delete old_requests[i];
  memory->sfree(old_requests);

  memory->destroy_2d_int_array(bondlist);
  memory->destroy_2d_int_array(anglelist);
  memory->destroy_2d_int_array(dihedrallist);
  memory->destroy_2d_int_array(improperlist);
}

/* ---------------------------------------------------------------------- */

void Neighbor::init()
{
  int i,j,m,n;

  ncalls = ndanger = 0;
  dimension = domain->dimension;
  triclinic = domain->triclinic;
  newton_pair = force->newton_pair;

  // error check

  if (delay > 0 && (delay % every) != 0)
    error->all("Neighbor delay must be 0 or multiple of every setting");

  if (pgsize < 10*oneatom)
    error->all("Neighbor page size must be >= 10x the one atom setting");

  // ------------------------------------------------------------------
  // settings

  // set neighbor cutoffs (force cutoff + skin)
  // trigger determines when atoms migrate and neighbor lists are rebuilt
  // cutneigh and cutneighsq determine what pairs go into neighbor list
  //   set to 0 if cutforce = 0
  // cutneighmin/max used for neighbor bin sizes
  // cutghost determines comm distance = max of cutneigh & skin
  //   skin is only > cutneighmax when cutneighmax = 0.0 due to pair = NULL
  //   in this case, are running a bond-only simulation
  //   still need ghost atom communication for bonds and atom exchange
  //     triggered by atoms moving 1/2 skin distance, but no neigh lists formed
  //   user should set skin big enough to find all bonds

  triggersq = 0.25*skin*skin;

  n = atom->ntypes;
  if (cutneighsq == NULL) {
    cutneighsq = memory->create_2d_double_array(n+1,n+1,"neigh:cutneighsq");
    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);
    }
  }
  cutghost = MAX(cutneighmax,skin);
  cutneighmaxsq = cutneighmax * cutneighmax;

  // check other classes that can induce reneighboring in decide()

  restart_check = 0;
  if (output->restart_every) 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

  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) special_flag[1] = special_flag[2] = special_flag[3] = 2;

  // rRESPA cutoffs

  int respa = 0;
  if (update->whichflag == 0 && strcmp(update->integrate_style,"respa") == 0) {
    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);
  }

  // ------------------------------------------------------------------
  // xhold, bins, exclusion lists

  // free xhold and bins if not needed for this run

  if (dist_check == 0) {
    memory->destroy_2d_double_array(xhold);
    maxhold = 0;
    xhold = NULL;
  }

  if (style == NSQ) {
    memory->sfree(bins);
    memory->sfree(binhead);
    maxbin = maxhead = 0;
    binhead = NULL;
    bins = NULL;
  }

  // 1st time allocation of xhold and bins

  if (dist_check) {
    if (maxhold == 0) {
      maxhold = atom->nmax;
      xhold = memory->create_2d_double_array(maxhold,3,"neigh:xhold");
    }
  }

  if (style != NSQ) {
    if (maxbin == 0) {
      maxbin = atom->nmax;
      bins = (int *) memory->smalloc(maxbin*sizeof(int),"bins");
    }
  }
    
  // exclusion lists for type, group, molecule settings from neigh_modify

  n = atom->ntypes;

  if (nex_type == 0 && nex_group == 0 && nex_mol == 0) exclude = 0;
  else exclude = 1;

  if (nex_type) {
    memory->destroy_2d_int_array(ex_type);
    ex_type = (int **) memory->create_2d_int_array(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("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) {
    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) {
    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]];
  }

  // ------------------------------------------------------------------
  // pairwise lists

  // test if pairwise lists need to be re-created
  // no need to re-create if:
  // neigh style has not changed and current requests = old requests

  int same = 1;
  if (style != old_style) same = 0;
  if (nrequest != old_nrequest) same = 0;
  else
    for (i = 0; i < nrequest; i++)
      if (requests[i]->identical(old_requests[i]) == 0) same = 0;

#ifdef NEIGH_LIST_DEBUG
  if (comm->me == 0) printf("SAME flag %d\n",same);
#endif

  // if old and new are not the same, create new pairwise lists

  if (!same) {

    // delete old lists and create new ones

    for (i = 0; i < nlist; i++) delete lists[i];
    delete [] lists;
    delete [] pair_build;
    delete [] stencil_create;

    nlist = nrequest;
    lists = new NeighList*[nlist];
    pair_build = new PairPtr[nlist];
    stencil_create = new StencilPtr[nlist];

    // create individual lists, one per request
    // copy dnum setting from request to list
    // pass list ptr back to requestor (except for Command class)

    for (i = 0; i < nlist; i++) {
      lists[i] = new NeighList(lmp,pgsize);
      lists[i]->index = i;
      lists[i]->dnum = requests[i]->dnum;

      if (requests[i]->pair) {
	Pair *pair = (Pair *) requests[i]->requestor;
	pair->init_list(requests[i]->id,lists[i]);
      } else if (requests[i]->fix) {
	Fix *fix = (Fix *) requests[i]->requestor;
	fix->init_list(requests[i]->id,lists[i]);
      } else if (requests[i]->compute) {
	Compute *compute = (Compute *) requests[i]->requestor;
	compute->init_list(requests[i]->id,lists[i]);
      }
    }

    // detect lists that are connected to other lists
    // if-then-else sequence is important
    //   since don't want to re-process skip or copy lists further down
    // copy: point this list at request->otherlist, could be a skip list
    // skip: point this list at request->otherlist, copy skip info from request
    // half_from_full: point this list at preceeding full list
    // granhistory: set preceeding list's listgranhistory to this list
    //   also set precedding list's ptr to FixShearHistory
    // respaouter: point this list at preceeding 1/2 inner/middle lists
    // pair and half: if there is a full non-occasional non-skip list
    //   change this list to half_from_full and point at the full list
    //   parent could be copy list or pair or fix
    // fix/compute requests:
    //   kind of request = half or full, occasional or not doesn't matter
    //   if request = half and non-skip pair half/respaouter exists,
    //     become copy of that list
    //   if request = full and non-skip pair full exists,
    //     become copy of that list
    //   if request = half and non-skip pair full exists,
    //     become half_from_full of that list
    //   if no matches, do nothing, fix/compute list will be built directly
    //   ok if parent is copy list

    for (i = 0; i < nlist; i++) {
      if (requests[i]->copy)
	lists[i]->listcopy = lists[requests[i]->otherlist];

      else if (requests[i]->skip) {
	lists[i]->listskip = lists[requests[i]->otherlist];
	lists[i]->copy_skip_info(requests[i]->iskip,requests[i]->ijskip);

      } else if (requests[i]->half_from_full)
	lists[i]->listfull = lists[i-1];

      else if (requests[i]->granhistory) {
	lists[i-1]->listgranhistory = lists[i];
	for (int ifix = 0; ifix < modify->nfix; ifix++)
	  if (strcmp(modify->fix[ifix]->style,"SHEAR_HISTORY") == 0) 
	    lists[i-1]->fix_history = (FixShearHistory *) modify->fix[ifix];
 
      } else if (requests[i]->respaouter) {
	if (requests[i-1]->respainner) {
	  lists[i]->respamiddle = 0;
	  lists[i]->listinner = lists[i-1];
	} else {
	  lists[i]->respamiddle = 1;
	  lists[i]->listmiddle = lists[i-1];
	  lists[i]->listinner = lists[i-2];
	}

      } else if (requests[i]->pair && requests[i]->half) {
	for (j = 0; j < nlist; j++)
	  if (requests[j]->full && requests[j]->occasional == 0 &&
	      requests[j]->skip == 0) break;
	if (j < nlist) {
	  requests[i]->half = 0;
	  requests[i]->half_from_full = 1;
	  lists[i]->listfull = lists[j];
	}

      } else if (requests[i]->fix || requests[i]->compute) {
	for (j = 0; j < nlist; j++) {
	  if (requests[i]->half && requests[j]->pair &&
	      requests[j]->skip == 0 && requests[j]->half) break;
	  if (requests[i]->full && requests[j]->pair &&
	      requests[j]->skip == 0 && requests[j]->full) break;
	  if (requests[i]->half && requests[j]->pair &&
	      requests[j]->skip == 0 && requests[j]->respaouter) break;
	}
	if (j < nlist) {
	  requests[i]->copy = 1;
	  lists[i]->listcopy = lists[j];
	} else {
	  for (j = 0; j < nlist; j++) {
	    if (requests[i]->half && requests[j]->pair &&
		requests[j]->skip == 0 && requests[j]->full) break;
	  }
	  if (j < nlist) {
	    requests[i]->half = 0;
	    requests[i]->half_from_full = 1;
	    lists[i]->listfull = lists[j];
	  }
	}
      }
    }

    // set ptrs to pair_build and stencil_create functions for each list
    // ptrs set to NULL if not set explicitly

    for (i = 0; i < nlist; i++) {
      choose_build(i,requests[i]);
      if (style != NSQ) choose_stencil(i,requests[i]);
      else stencil_create[i] = NULL;
    }

    // set each list's build/grow/stencil flags based on neigh request
    // buildflag = 1 if its pair_build() invoked every reneighbor
    // growflag = 1 if it stores atom-based arrays and pages
    // stencilflag = 1 if it stores stencil arrays

    for (i = 0; i < nlist; i++) {
      lists[i]->buildflag = 1;
      if (pair_build[i] == NULL) lists[i]->buildflag = 0;
      if (requests[i]->occasional) lists[i]->buildflag = 0;

      lists[i]->growflag = 1;
      if (requests[i]->copy) lists[i]->growflag = 0;

      lists[i]->stencilflag = 1;
      if (style == NSQ) lists[i]->stencilflag = 0;
      if (stencil_create[i] == NULL) lists[i]->stencilflag = 0;
    }

#ifdef NEIGH_LIST_DEBUG
    for (i = 0; i < nlist; i++) lists[i]->print_attributes();
#endif

    // allocate atom arrays and 1st pages of lists that store them

    maxlocal = atom->nmax;
    for (i = 0; i < nlist; i++)
      if (lists[i]->growflag) {
	lists[i]->grow(maxlocal);
	lists[i]->add_pages();
      }

    // setup 3 vectors of pairwise neighbor lists
    // blist = lists whose pair_build() is invoked every reneighbor
    // glist = lists who store atom arrays which are used every reneighbor
    // slist = lists who store stencil arrays which are used every reneighbor
    // blist and glist vectors are used by neighbor::build()
    // slist vector is used by neighbor::setup_bins()

    nblist = nglist = nslist = 0;
    delete [] blist;
    delete [] glist;
    delete [] slist;
    blist = new int[nlist];
    glist = new int[nlist];
    slist = new int[nlist];

    for (i = 0; i < nlist; i++) {
      if (lists[i]->buildflag) blist[nblist++] = i;
      if (lists[i]->growflag && requests[i]->occasional == 0)
	glist[nglist++] = i;
      if (lists[i]->stencilflag && requests[i]->occasional == 0)
	slist[nslist++] = i;
    }

#ifdef NEIGH_LIST_DEBUG
    print_lists_of_lists();
#endif

    // reorder build vector if necessary
    // relevant for lists that copy/skip/half-full from parent
    // the derived list must appear in blist after the parent list
    // no occasional lists are in build vector
    // swap two lists within blist when dependency is mis-ordered
    // done when entire pass thru blist results in no swaps

    int done = 0;
    while (!done) {
      done = 1;
      for (i = 0; i < nblist; i++) {
	NeighList *ptr = NULL;
	if (lists[blist[i]]->listfull) ptr = lists[blist[i]]->listfull;
	if (lists[blist[i]]->listcopy) ptr = lists[blist[i]]->listcopy;
	if (lists[blist[i]]->listskip) ptr = lists[blist[i]]->listskip;
	if (ptr == NULL) continue;
	for (m = 0; m < nlist; m++)
	  if (ptr == lists[m]) break;
	for (j = 0; j < nblist; j++)
	  if (m == blist[j]) break;
	if (j < i) continue;
	int tmp = blist[i];
	blist[i] = blist[j];
	blist[j] = tmp;
	done = 0;
	break;
      }
    }

#ifdef NEIGH_LIST_DEBUG
    print_lists_of_lists();
#endif
  }

  // delete old requests
  // copy current requests and style to old for next run
  
  for (i = 0; i < old_nrequest; i++) delete old_requests[i];
  memory->sfree(old_requests);
  old_nrequest = nrequest;
  old_requests = requests;
  nrequest = maxrequest = 0;
  requests = NULL;
  old_style = style;

  // ------------------------------------------------------------------
  // topology lists

  // 1st time allocation of topology lists

  if (atom->molecular && atom->nbonds && maxbond == 0) {
    if (nprocs == 1) maxbond = atom->nbonds;
    else maxbond = static_cast<int> (LB_FACTOR * atom->nbonds / nprocs);
    bondlist = memory->create_2d_int_array(maxbond,3,"neigh:bondlist");
  }

  if (atom->molecular && atom->nangles && maxangle == 0) {
    if (nprocs == 1) maxangle = atom->nangles;
    else maxangle = static_cast<int> (LB_FACTOR * atom->nangles / nprocs);
    anglelist =  memory->create_2d_int_array(maxangle,4,"neigh:anglelist");
  }

  if (atom->molecular && atom->ndihedrals && maxdihedral == 0) {
    if (nprocs == 1) maxdihedral = atom->ndihedrals;
    else maxdihedral = static_cast<int> 
	   (LB_FACTOR * atom->ndihedrals / nprocs);
    dihedrallist = 
      memory->create_2d_int_array(maxdihedral,5,"neigh:dihedrallist");
  }

  if (atom->molecular && atom->nimpropers && maximproper == 0) {
    if (nprocs == 1) maximproper = atom->nimpropers;
    else maximproper = static_cast<int>
	   (LB_FACTOR * atom->nimpropers / nprocs);
    improperlist = 
      memory->create_2d_int_array(maximproper,5,"neigh:improperlist");
  }

  // set flags that determine which topology neighboring routines to use
  // SHAKE sets bonds and angles 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)
      bond_off = angle_off = 1;
  if (force->bond && force->bond_match("quartic")) bond_off = 1;

  if (atom->avec->bonds_allow) {
    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) {
    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) {
    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) {
    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;
    }
  }

  // set ptrs to topology build functions

  if (bond_off) bond_build = &Neighbor::bond_partial;
  else bond_build = &Neighbor::bond_all;

  if (angle_off) angle_build = &Neighbor::angle_partial;
  else angle_build = &Neighbor::angle_all;

  if (dihedral_off) dihedral_build = &Neighbor::dihedral_partial;
  else dihedral_build = &Neighbor::dihedral_all;

  if (improper_off) improper_build = &Neighbor::improper_partial;
  else improper_build = &Neighbor::improper_all;

  // set topology neighbor list counts to 0
  // in case all are turned off but potential is still defined

  nbondlist = nanglelist = ndihedrallist = nimproperlist = 0;
}

/* ---------------------------------------------------------------------- */

int Neighbor::request(void *requestor)
{
  if (nrequest == maxrequest) {
    maxrequest += RQDELTA;
    requests = (NeighRequest **) 
      memory->srealloc(requests,maxrequest*sizeof(NeighRequest *),
		       "neighbor:requests");
  }

  requests[nrequest] = new NeighRequest(lmp);
  requests[nrequest]->requestor = requestor;
  nrequest++;
  return nrequest-1;
}

/* ----------------------------------------------------------------------
   determine which pair_build function each neigh list needs
   based on settings of neigh request
   copy -> copy_from function
   skip -> granular function if gran with granhistory
           respa function if respaouter
	   skip_from function for everything else
   half_from_full, half, full, gran, respaouter -> choose by newton and tri
     style NSQ options = newton off, newton on
     style BIN options = newton off, newton on and not tri, newton on and tri
     stlye MULTI options = same options as BIN
   if none of these, ptr = NULL since pair_build is not invoked for this list
   use "else if" b/c skip,copy can be set in addition to half,full,etc
------------------------------------------------------------------------- */

void Neighbor::choose_build(int index, NeighRequest *rq)
{
  PairPtr pb = NULL;

  if (rq->copy) pb = &Neighbor::copy_from;

  else if (rq->skip) {
    if (rq->gran && lists[index]->listgranhistory)
      pb = &Neighbor::skip_from_granular;
    else if (rq->respaouter) pb = &Neighbor::skip_from_respa;
    else pb = &Neighbor::skip_from;

  } else if (rq->half_from_full) {
    if (newton_pair == 0) pb = &Neighbor::half_from_full_no_newton;
    else if (newton_pair == 1) pb = &Neighbor::half_from_full_newton;

  } else if (rq->half) {
    if (style == NSQ) {
      if (newton_pair == 0) pb = &Neighbor::half_nsq_no_newton;
      else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton;
    } else if (style == BIN) {
      if (newton_pair == 0) pb = &Neighbor::half_bin_no_newton;
      else if (triclinic == 0) pb = &Neighbor::half_bin_newton;
      else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri;
    } else if (style == MULTI) {
      if (newton_pair == 0) pb = &Neighbor::half_multi_no_newton;
      else if (triclinic == 0) pb = &Neighbor::half_multi_newton;
      else if (triclinic == 1) pb = &Neighbor::half_multi_newton_tri;
    }

  } else if (rq->full) {
    if (style == NSQ) pb = &Neighbor::full_nsq;
    else if (style == BIN) pb = &Neighbor::full_bin;
    else if (style == MULTI) pb = &Neighbor::full_multi;

  } else if (rq->gran) {
    if (style == NSQ) {
      if (newton_pair == 0) pb = &Neighbor::granular_nsq_no_newton;
      else if (newton_pair == 1) pb = &Neighbor::granular_nsq_newton;
    } else if (style == BIN) {
      if (newton_pair == 0) pb = &Neighbor::granular_bin_no_newton;
      else if (triclinic == 0) pb = &Neighbor::granular_bin_newton;
      else if (triclinic == 1) pb = &Neighbor::granular_bin_newton_tri;
    } else if (style == MULTI)
      error->all("Neighbor multi not yet enabled for granular");

  } else if (rq->respaouter) {
    if (style == NSQ) {
      if (newton_pair == 0) pb = &Neighbor::respa_nsq_no_newton;
      else if (newton_pair == 1) pb = &Neighbor::respa_nsq_newton;
    } else if (style == BIN) {
      if (newton_pair == 0) pb = &Neighbor::respa_bin_no_newton;
      else if (triclinic == 0) pb = &Neighbor::respa_bin_newton;
      else if (triclinic == 1) pb = &Neighbor::respa_bin_newton_tri;
    } else if (style == MULTI)
      error->all("Neighbor multi not yet enabled for rRESPA");
  }

  pair_build[index] = pb;
}

/* ----------------------------------------------------------------------
   determine which stencil_create function each neigh list needs
   based on settings of neigh request, only called if style != NSQ
   skip or copy or half_from_full -> no stencil
   half, gran, respaouter, full -> choose by newton and tri and dimension
   if none of these, ptr = NULL since this list needs no stencils
   use "else if" b/c skip,copy can be set in addition to half,full,etc
------------------------------------------------------------------------- */

void Neighbor::choose_stencil(int index, NeighRequest *rq)
{
  StencilPtr sc = NULL;

  if (rq->skip || rq->copy || rq->half_from_full) sc = NULL;

  else if (rq->half || rq->gran || rq->respaouter) {
    if (style == BIN) {
      if (newton_pair == 0) {
	if (dimension == 2) sc = &Neighbor::stencil_half_bin_2d_no_newton;
	else if (dimension == 3) sc = &Neighbor::stencil_half_bin_3d_no_newton;
      } else if (triclinic == 0) {
	if (dimension == 2) sc = &Neighbor::stencil_half_bin_2d_newton;
	else if (dimension == 3) sc = &Neighbor::stencil_half_bin_3d_newton;
      } else if (triclinic == 1) {
	if (dimension == 2) sc = &Neighbor::stencil_half_bin_2d_newton_tri;
	else if (dimension == 3) 
	  sc = &Neighbor::stencil_half_bin_3d_newton_tri;
      }
    } else if (style == MULTI) {
      if (newton_pair == 0) {
	if (dimension == 2) sc = &Neighbor::stencil_half_multi_2d_no_newton;
	else if (dimension == 3)
	  sc = &Neighbor::stencil_half_multi_3d_no_newton;
      } else if (triclinic == 0) {
	if (dimension == 2) sc = &Neighbor::stencil_half_multi_2d_newton;
	else if (dimension == 3) sc = &Neighbor::stencil_half_multi_3d_newton;
      } else if (triclinic == 1) {
	if (dimension == 2) sc = &Neighbor::stencil_half_multi_2d_newton_tri;
	else if (dimension == 3) 
	  sc = &Neighbor::stencil_half_multi_3d_newton_tri;
      }
    }

  } else if (rq->full) {
    if (style == BIN) {
      if (dimension == 2) sc = &Neighbor::stencil_full_bin_2d;
      else if (dimension == 3) sc = &Neighbor::stencil_full_bin_3d;
    } else if (style == MULTI) {
      if (dimension == 2) sc = &Neighbor::stencil_full_multi_2d;
      else if (dimension == 3) sc = &Neighbor::stencil_full_multi_3d;
    }
  }

  stencil_create[index] = sc;
}

/* ---------------------------------------------------------------------- */

void Neighbor::print_lists_of_lists()
{
  if (comm->me == 0) {
    printf("Build lists = %d: ",nblist);
    for (int i = 0; i < nblist; i++) printf("%d ",blist[i]);
    printf("\n");
    printf("Grow lists = %d: ",nglist);
    for (int i = 0; i < nglist; i++) printf("%d ",glist[i]);
    printf("\n");
    printf("Stencil lists = %d: ",nslist);
    for (int i = 0; i < nslist; i++) printf("%d ",slist[i]);
    printf("\n");
  }
}

/* ---------------------------------------------------------------------- */

int Neighbor::decide()
{
  if (must_check) {
    int 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 (dist_check == 0) return 1;
    else return check_distance();
  } else return 0;
}

/* ---------------------------------------------------------------------- */

int Neighbor::check_distance()
{
  double delx,dely,delz,rsq;

  double **x = atom->x;
  int nlocal = atom->nlocal;
  if (include_group) 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 > triggersq) flag = 1;
  }

  int flagall;
  MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
  if (flagall && ago == MAX(every,delay)) ndanger++;
  return flagall;
}

/* ----------------------------------------------------------------------
   build all perpetual neighbor lists every few timesteps
   pairwise & topology lists are created as needed
------------------------------------------------------------------------- */

void Neighbor::build()
{
  int i;

  ago = 0;
  ncalls++;

  // store current atom positions if needed

  if (dist_check) {
    double **x = atom->x;
    int nlocal = atom->nlocal;
    if (nlocal > maxhold) {
      maxhold = atom->nmax;
      memory->destroy_2d_double_array(xhold);
      xhold = memory->create_2d_double_array(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 necessary, extend atom arrays in pairwise lists
  // only for lists with growflag set and which are used every reneighbor

  if (atom->nlocal > maxlocal) {
    maxlocal = atom->nmax;
    for (i = 0; i < nglist; i++) lists[glist[i]]->grow(maxlocal);
  }

  // extend atom bin list if necessary

  if (style != NSQ && atom->nmax > maxbin) {
    maxbin = atom->nmax;
    memory->sfree(bins);
    bins = (int *) memory->smalloc(maxbin*sizeof(int),"bins");
  }

  // invoke building of pair and molecular neighbor lists
  // only for pairwise lists with buildflag set

  for (i = 0; i < nblist; i++)
    (this->*pair_build[blist[i]])(lists[blist[i]]);

  if (atom->molecular) {
    if (force->bond) (this->*bond_build)();
    if (force->angle) (this->*angle_build)();
    if (force->dihedral) (this->*dihedral_build)();
    if (force->improper) (this->*improper_build)();