remove USER-DPD files that should not be where they are

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

   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 authors: 

   James Larentzos and Timothy I. Mattox (Engility Corporation)

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

#include "npair_half_bin_newton_ssa.h"

#include "neighbor.h"

#include "neigh_list.h"

#include "atom.h"

#include "atom_vec.h"

#include "molecule.h"

#include "domain.h"

#include "group.h"

#include "memory.h"

#include "my_page.h"

#include "error.h"

using namespace LAMMPS_NS;

// allocate space for static class variable

// prototype for non-class function

static int *ssaAIRptr;

static int cmp_ssaAIR(const void *, const void *);

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

NPairHalfBinNewtonSSA::NPairHalfBinNewtonSSA(LAMMPS *lmp) : NPair(lmp) {}

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

   binned neighbor list construction with full Newton's 3rd law

   for use by Shardlow Spliting Algorithm

   each owned atom i checks its own bin and other bins in Newton stencil

   every pair stored exactly once by some processor

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

void NPairHalfBinNewtonSSA::build(NeighList *list)

{

  int i,j,k,n,itype,jtype,ibin,which,imol,iatom,moltemplate;

  tagint tagprev;

  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;

  int *neighptr;

  double **x = atom->x;

  int *type = atom->type;

  int *mask = atom->mask;

  tagint *tag = atom->tag;

  tagint *molecule = atom->molecule;

  tagint **special = atom->special;

  int **nspecial = atom->nspecial;

  int nlocal = atom->nlocal;

  int nall = nlocal + atom->nghost;

  if (includegroup) nlocal = atom->nfirst;

  int *ssaAIR = atom->ssaAIR;

  int *molindex = atom->molindex;

  int *molatom = atom->molatom;

  Molecule **onemols = atom->avec->onemols;

  int molecular = atom->molecular;

  if (molecular == 2) moltemplate = 1;

  else moltemplate = 0;

  int *ilist = list->ilist;

  int *numneigh = list->numneigh;

  int **firstneigh = list->firstneigh;

  MyPage<int> *ipage = list->ipage;

  int inum = 0;

  // bin owned and ghost atoms for use by Shardlow Splitting Algorithm

  // exclude ghost atoms that are not in the Active Interaction Regions (AIR)

  // NOTE to Tim: this binatomflag no longer exists

  //   the logic up higher assures that binning has been done

  //     before this build() method is called

  //   maybe this code below needs to be in a new NBinShardlow class?

  // this class also inherits NPair::nb from its parent

  //   which points to the NBin class that did the binning

  //   there are last_step variables stored there which indicate

  //   the last time binning was done

  // the basic question is what data is created/stored by SSA binning

  //   and in what class should it live?

  //   if it is created by the binning operation, then I think

  //     it should be in a new NBinShardlow class

  if (true /* binatomflag */) { // only false in Neighbor::build_one

    if (mbins > list->maxhead_ssa) {

      list->maxhead_ssa = mbins;

      memory->destroy(list->gbinhead_ssa);

      memory->destroy(list->binhead_ssa);

      memory->create(list->binhead_ssa,list->maxhead_ssa,"binhead_ssa");

      memory->create(list->gbinhead_ssa,list->maxhead_ssa,"gbinhead_ssa");

    }

    for (i = 0; i < mbins; i++) {

      list->gbinhead_ssa[i] = -1;

      list->binhead_ssa[i] = -1;

    }

    if (nall > list->maxbin_ssa) {

      list->maxbin_ssa = nall;

      memory->destroy(list->bins_ssa);

      memory->create(list->bins_ssa,list->maxbin_ssa,"bins_ssa");

    }

    // bin in reverse order so linked list will be in forward order

    if (includegroup) {

      int bitmask = group->bitmask[includegroup];

      int nowned = atom->nlocal; // NOTE: nlocal was set to atom->nfirst above

      for (i = nall-1; i >= nowned; i--) {

        if (ssaAIR[i] < 2) continue; // skip ghost atoms not in AIR

        if (mask[i] & bitmask) {

          ibin = coord2bin(x[i]);

          list->bins_ssa[i] = list->gbinhead_ssa[ibin];

          list->gbinhead_ssa[ibin] = i;

        }

      }

    } else {

      for (i = nall-1; i >= nlocal; i--) {

        if (ssaAIR[i] < 2) continue; // skip ghost atoms not in AIR

        ibin = coord2bin(x[i]);

        list->bins_ssa[i] = list->gbinhead_ssa[ibin];

        list->gbinhead_ssa[ibin] = i;

      }

    }

    for (i = nlocal-1; i >= 0; i--) {

      ibin = coord2bin(x[i]);

      list->bins_ssa[i] = list->binhead_ssa[ibin];

      list->binhead_ssa[ibin] = i;

    }

  }  // else reuse previous binning. See Neighbor::build_one comment

  ipage->reset();

  // loop over owned atoms, storing half of the neighbors

  for (i = 0; i < nlocal; i++) {

    int AIRct[8] = { 0 };

    n = 0;

    neighptr = ipage->vget();

    itype = type[i];

    xtmp = x[i][0];

    ytmp = x[i][1];

    ztmp = x[i][2];

    if (moltemplate) {

      imol = molindex[i];

      iatom = molatom[i];

      tagprev = tag[i] - iatom - 1;

    }

    // loop over rest of local atoms in i's bin

    // just store them, since j is beyond i in linked list

    for (j = list->bins_ssa[i]; j >= 0; j = list->bins_ssa[j]) {

      jtype = type[j];

      if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;

      delx = xtmp - x[j][0];

      dely = ytmp - x[j][1];

      delz = ztmp - x[j][2];

      rsq = delx*delx + dely*dely + delz*delz;

      if (rsq <= cutneighsq[itype][jtype]) {

        if (molecular) {

          if (!moltemplate)

            which = find_special(special[i],nspecial[i],tag[j]);

          else if (imol >= 0)

            which = find_special(onemols[imol]->special[iatom],

                                 onemols[imol]->nspecial[iatom],

                                 tag[j]-tagprev);

          else which = 0;

          if (which == 0) neighptr[n++] = j;

          else if (domain->minimum_image_check(delx,dely,delz))

            neighptr[n++] = j;

          else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);

        } else neighptr[n++] = j;

      }

    }

    ibin = coord2bin(x[i]);

    // loop over all local atoms in other bins in "half" stencil

    for (k = 0; k < nstencil; k++) {

      for (j = list->binhead_ssa[ibin+stencil[k]]; j >= 0; 

           j = list->bins_ssa[j]) {

        jtype = type[j];

        if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;

        delx = xtmp - x[j][0];

        dely = ytmp - x[j][1];

        delz = ztmp - x[j][2];

        rsq = delx*delx + dely*dely + delz*delz;

        if (rsq <= cutneighsq[itype][jtype]) {

          if (molecular) {

            if (!moltemplate)

              which = find_special(special[i],nspecial[i],tag[j]);

            else if (imol >= 0)

              which = find_special(onemols[imol]->special[iatom],

                                   onemols[imol]->nspecial[iatom],

                                   tag[j]-tagprev);

            else which = 0;

            if (which == 0) neighptr[n++] = j;

            else if (domain->minimum_image_check(delx,dely,delz))

              neighptr[n++] = j;

            else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);

          } else neighptr[n++] = j;

        }

      }

    }

    AIRct[0] = n;

    // loop over AIR ghost atoms in all bins in "full" stencil

    // Note: the non-AIR ghost atoms have already been filtered out

    // That is a significant time savings because of the "full" stencil

    // Note2: only non-pure locals can have ghosts as neighbors

    if (ssaAIR[i] == 1) for (k = 0; k < nstencil_ssa; k++) {

      for (j = list->gbinhead_ssa[ibin+stencil[k]]; j >= 0; 

           j = list->bins_ssa[j]) {

        jtype = type[j];

        if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;

        delx = xtmp - x[j][0];

        dely = ytmp - x[j][1];

        delz = ztmp - x[j][2];

        rsq = delx*delx + dely*dely + delz*delz;

        if (rsq <= cutneighsq[itype][jtype]) {

          if (molecular) {

            if (!moltemplate)

              which = find_special(special[i],nspecial[i],tag[j]);

            else if (imol >= 0)

              which = find_special(onemols[imol]->special[iatom],

                                   onemols[imol]->nspecial[iatom],

                                   tag[j]-tagprev);

            else which = 0;

            if (which == 0) {

              neighptr[n++] = j;

              ++(AIRct[ssaAIR[j] - 1]);

            } else if (domain->minimum_image_check(delx,dely,delz)) {

              neighptr[n++] = j;

              ++(AIRct[ssaAIR[j] - 1]);

            } else if (which > 0) {

              neighptr[n++] = j ^ (which << SBBITS);

              ++(AIRct[ssaAIR[j] - 1]);

            }

          } else {

            neighptr[n++] = j;

            ++(AIRct[ssaAIR[j] - 1]);

          }

        }

      }

    }

    ilist[inum++] = i;

    firstneigh[i] = neighptr;

    numneigh[i] = n;

    ipage->vgot(n);

    if (ipage->status())

      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");

    // sort the ghosts in the neighbor list by their ssaAIR number

    ssaAIRptr = atom->ssaAIR;

    qsort(&(neighptr[AIRct[0]]), n - AIRct[0], sizeof(int), cmp_ssaAIR);

    // do a prefix sum on the counts to turn them into indexes

    list->ndxAIR_ssa[i][0] = AIRct[0];

    for (int ndx = 1; ndx < 8; ++ndx) {

      list->ndxAIR_ssa[i][ndx] = AIRct[ndx] + list->ndxAIR_ssa[i][ndx - 1];

    }

  }

  list->inum = inum;

}

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

   comparison function invoked by qsort()

   accesses static class member ssaAIRptr, set before call to qsort()

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

static int cmp_ssaAIR(const void *iptr, const void *jptr)

{

  int i = *((int *) iptr);

  int j = *((int *) jptr);

  if (ssaAIRptr[i] < ssaAIRptr[j]) return -1;

  if (ssaAIRptr[i] > ssaAIRptr[j]) return 1;

  return 0;

}

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

   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 authors: 

   James Larentzos and Timothy I. Mattox (Engility Corporation)

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

#include "npair_half_bin_newton_ssa.h"

#include "neighbor.h"

#include "neigh_list.h"

#include "atom.h"

#include "atom_vec.h"

#include "molecule.h"

#include "domain.h"

#include "group.h"

#include "memory.h"

#include "my_page.h"

#include "error.h"

using namespace LAMMPS_NS;

// allocate space for static class variable

// prototype for non-class function

static int *ssaAIRptr;

static int cmp_ssaAIR(const void *, const void *);

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

NPairHalfBinNewtonSSA::NPairHalfBinNewtonSSA(LAMMPS *lmp) : NPair(lmp) {}

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

   binned neighbor list construction with full Newton's 3rd law

   for use by Shardlow Spliting Algorithm

   each owned atom i checks its own bin and other bins in Newton stencil

   every pair stored exactly once by some processor

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

void NPairHalfBinNewtonSSA::build(NeighList *list)

{

  int i,j,k,n,itype,jtype,ibin,which,imol,iatom,moltemplate;

  tagint tagprev;

  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;

  int *neighptr;

  double **x = atom->x;

  int *type = atom->type;

  int *mask = atom->mask;

  tagint *tag = atom->tag;

  tagint *molecule = atom->molecule;

  tagint **special = atom->special;

  int **nspecial = atom->nspecial;

  int nlocal = atom->nlocal;

  int nall = nlocal + atom->nghost;

  if (includegroup) nlocal = atom->nfirst;

  int *ssaAIR = atom->ssaAIR;

  int *molindex = atom->molindex;

  int *molatom = atom->molatom;

  Molecule **onemols = atom->avec->onemols;

  int molecular = atom->molecular;

  if (molecular == 2) moltemplate = 1;

  else moltemplate = 0;

  int *ilist = list->ilist;

  int *numneigh = list->numneigh;

  int **firstneigh = list->firstneigh;

  MyPage<int> *ipage = list->ipage;

  int inum = 0;

  // bin owned and ghost atoms for use by Shardlow Splitting Algorithm

  // exclude ghost atoms that are not in the Active Interaction Regions (AIR)

  // NOTE to Tim: this binatomflag no longer exists

  //   the logic up higher assures that binning has been done

  //     before this build() method is called

  //   maybe this code below needs to be in a new NBinShardlow class?

  // this class also inherits NPair::nb from its parent

  //   which points to the NBin class that did the binning

  //   there are last_step variables stored there which indicate

  //   the last time binning was done

  // the basic question is what data is created/stored by SSA binning

  //   and in what class should it live?

  //   if it is created by the binning operation, then I think

  //     it should be in a new NBinShardlow class

  if (true /* binatomflag */) { // only false in Neighbor::build_one

    if (mbins > list->maxhead_ssa) {

      list->maxhead_ssa = mbins;

      memory->destroy(list->gbinhead_ssa);

      memory->destroy(list->binhead_ssa);

      memory->create(list->binhead_ssa,list->maxhead_ssa,"binhead_ssa");

      memory->create(list->gbinhead_ssa,list->maxhead_ssa,"gbinhead_ssa");

    }

    for (i = 0; i < mbins; i++) {

      list->gbinhead_ssa[i] = -1;

      list->binhead_ssa[i] = -1;

    }

    if (nall > list->maxbin_ssa) {

      list->maxbin_ssa = nall;

      memory->destroy(list->bins_ssa);

      memory->create(list->bins_ssa,list->maxbin_ssa,"bins_ssa");

    }

    // bin in reverse order so linked list will be in forward order

    if (includegroup) {

      int bitmask = group->bitmask[includegroup];

      int nowned = atom->nlocal; // NOTE: nlocal was set to atom->nfirst above

      for (i = nall-1; i >= nowned; i--) {

        if (ssaAIR[i] < 2) continue; // skip ghost atoms not in AIR

        if (mask[i] & bitmask) {

          ibin = coord2bin(x[i]);

          list->bins_ssa[i] = list->gbinhead_ssa[ibin];

          list->gbinhead_ssa[ibin] = i;

        }

      }

    } else {

      for (i = nall-1; i >= nlocal; i--) {

        if (ssaAIR[i] < 2) continue; // skip ghost atoms not in AIR

        ibin = coord2bin(x[i]);

        list->bins_ssa[i] = list->gbinhead_ssa[ibin];

        list->gbinhead_ssa[ibin] = i;

      }

    }

    for (i = nlocal-1; i >= 0; i--) {

      ibin = coord2bin(x[i]);

      list->bins_ssa[i] = list->binhead_ssa[ibin];

      list->binhead_ssa[ibin] = i;

    }

  }  // else reuse previous binning. See Neighbor::build_one comment

  ipage->reset();

  // loop over owned atoms, storing half of the neighbors

  for (i = 0; i < nlocal; i++) {

    int AIRct[8] = { 0 };

    n = 0;

    neighptr = ipage->vget();

    itype = type[i];

    xtmp = x[i][0];

    ytmp = x[i][1];

    ztmp = x[i][2];

    if (moltemplate) {

      imol = molindex[i];

      iatom = molatom[i];

      tagprev = tag[i] - iatom - 1;

    }

    // loop over rest of local atoms in i's bin

    // just store them, since j is beyond i in linked list

    for (j = list->bins_ssa[i]; j >= 0; j = list->bins_ssa[j]) {

      jtype = type[j];

      if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;

      delx = xtmp - x[j][0];

      dely = ytmp - x[j][1];

      delz = ztmp - x[j][2];

      rsq = delx*delx + dely*dely + delz*delz;

      if (rsq <= cutneighsq[itype][jtype]) {

        if (molecular) {

          if (!moltemplate)

            which = find_special(special[i],nspecial[i],tag[j]);

          else if (imol >= 0)

            which = find_special(onemols[imol]->special[iatom],

                                 onemols[imol]->nspecial[iatom],

                                 tag[j]-tagprev);

          else which = 0;

          if (which == 0) neighptr[n++] = j;

          else if (domain->minimum_image_check(delx,dely,delz))

            neighptr[n++] = j;

          else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);

        } else neighptr[n++] = j;

      }

    }

    ibin = coord2bin(x[i]);

    // loop over all local atoms in other bins in "half" stencil

    for (k = 0; k < nstencil; k++) {

      for (j = list->binhead_ssa[ibin+stencil[k]]; j >= 0; 

           j = list->bins_ssa[j]) {

        jtype = type[j];

        if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;

        delx = xtmp - x[j][0];

        dely = ytmp - x[j][1];

        delz = ztmp - x[j][2];

        rsq = delx*delx + dely*dely + delz*delz;

        if (rsq <= cutneighsq[itype][jtype]) {

          if (molecular) {

            if (!moltemplate)

              which = find_special(special[i],nspecial[i],tag[j]);

            else if (imol >= 0)

              which = find_special(onemols[imol]->special[iatom],

                                   onemols[imol]->nspecial[iatom],

                                   tag[j]-tagprev);

            else which = 0;

            if (which == 0) neighptr[n++] = j;

            else if (domain->minimum_image_check(delx,dely,delz))

              neighptr[n++] = j;

            else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);

          } else neighptr[n++] = j;

        }

      }

    }

    AIRct[0] = n;

    // loop over AIR ghost atoms in all bins in "full" stencil

    // Note: the non-AIR ghost atoms have already been filtered out

    // That is a significant time savings because of the "full" stencil

    // Note2: only non-pure locals can have ghosts as neighbors

    if (ssaAIR[i] == 1) for (k = 0; k < nstencil_ssa; k++) {

      for (j = list->gbinhead_ssa[ibin+stencil[k]]; j >= 0; 

           j = list->bins_ssa[j]) {

        jtype = type[j];

        if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;

        delx = xtmp - x[j][0];

        dely = ytmp - x[j][1];

        delz = ztmp - x[j][2];

        rsq = delx*delx + dely*dely + delz*delz;

        if (rsq <= cutneighsq[itype][jtype]) {

          if (molecular) {

            if (!moltemplate)

              which = find_special(special[i],nspecial[i],tag[j]);

            else if (imol >= 0)

              which = find_special(onemols[imol]->special[iatom],

                                   onemols[imol]->nspecial[iatom],

                                   tag[j]-tagprev);

            else which = 0;

            if (which == 0) {

              neighptr[n++] = j;

              ++(AIRct[ssaAIR[j] - 1]);

            } else if (domain->minimum_image_check(delx,dely,delz)) {

              neighptr[n++] = j;

              ++(AIRct[ssaAIR[j] - 1]);

            } else if (which > 0) {

              neighptr[n++] = j ^ (which << SBBITS);

              ++(AIRct[ssaAIR[j] - 1]);

            }

          } else {

            neighptr[n++] = j;

            ++(AIRct[ssaAIR[j] - 1]);

          }

        }

      }

    }

    ilist[inum++] = i;

    firstneigh[i] = neighptr;

    numneigh[i] = n;

    ipage->vgot(n);

    if (ipage->status())

      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");

    // sort the ghosts in the neighbor list by their ssaAIR number

    ssaAIRptr = atom->ssaAIR;

    qsort(&(neighptr[AIRct[0]]), n - AIRct[0], sizeof(int), cmp_ssaAIR);

    // do a prefix sum on the counts to turn them into indexes

    list->ndxAIR_ssa[i][0] = AIRct[0];

    for (int ndx = 1; ndx < 8; ++ndx) {

      list->ndxAIR_ssa[i][ndx] = AIRct[ndx] + list->ndxAIR_ssa[i][ndx - 1];

    }

  }

  list->inum = inum;

}

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

   comparison function invoked by qsort()

   accesses static class member ssaAIRptr, set before call to qsort()

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

static int cmp_ssaAIR(const void *iptr, const void *jptr)

{

  int i = *((int *) iptr);

  int j = *((int *) jptr);

  if (ssaAIRptr[i] < ssaAIRptr[j]) return -1;

  if (ssaAIRptr[i] > ssaAIRptr[j]) return 1;

  return 0;

}

/* -*- c++ -*- ----------------------------------------------------------

   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.

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

#ifdef NPAIR_CLASS

NPairStyle(half/bin/newton/ssa,

           NPairHalfBinNewtonSSA,

           NP_BIN | NP_NEWTON | NP_ORTHO | NP_SSA)

#else

#ifndef LMP_NPAIR_HALF_BIN_NEWTON_SSA_H

#define LMP_NPAIR_HALF_BIN_NEWTON_SSA_H

#include "npair.h"

namespace LAMMPS_NS {

class NPairHalfBinNewtonSSA : public NPair {

 public:

  NPairHalfBinNewtonSSA(class LAMMPS *);

  ~NPairHalfBinNewtonSSA() {}

  void build(class NeighList *);

};

}

#endif

#endif

/* ERROR/WARNING messages:

*/

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

   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 authors: 

   James Larentzos and Timothy I. Mattox (Engility Corporation)

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

#include "nstencil_half_bin_2d_newton_ssa.h"

#include "neighbor.h"

#include "neigh_list.h"

using namespace LAMMPS_NS;

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

NStencilHalfBin2dNewtonSSA::NStencilHalfBin2dNewtonSSA(LAMMPS *lmp) : 

  NStencil(lmp) {}

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

   create stencil based on bin geometry and cutoff

   stencil = bins whose closest corner to central bin is within cutoff

   sx,sy,sz = bin bounds = furthest the stencil could possibly extend

   3d creates xyz stencil, 2d creates xy stencil

   for half list with newton on:

     stencil is bins to the "upper right" of central bin

     stencil does not include self

   additionally, includes the bins beyond nstencil that are needed

     to locate all the Active Interaction Region (AIR) ghosts for SSA

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

void NStencilHalfBin2dNewtonSSA::create()

{

  int i,j,pos = 0;

  for (j = 0; j <= sy; j++)

    for (i = -sx; i <= sx; i++)

      if (j > 0 || (j == 0 && i > 0))

        if (bin_distance(i,j,0) < cutneighmaxsq)

          stencil[pos++] = j*mbinx + i;

  nstencil = pos; // record where normal half stencil ends

  // include additional bins for AIR ghosts only

  for (j = -sy; j <= 0; j++)

    for (i = -sx; i <= sx; i++) {

      if (j == 0 && i > 0) continue;

      if (bin_distance(i,j,0) < cutneighmaxsq)

        stencil[pos++] = j*mbinx + i;

    }

  nstencil_ssa = pos; // record where full stencil ends

}