Merge branch 'master' into collected-small-changes

doc/src/comm_modify.rst

@@ -84,13 +84,15 @@ information is available, then also a heuristic based on that bond length
 is computed. It is used as communication cutoff, if there is no pair
 style present and no *comm_modify cutoff* command used. Otherwise a
 warning is printed, if this bond based estimate is larger than the
-communication cutoff used. A
+communication cutoff used. 
 
 The *cutoff/multi* option is equivalent to *cutoff*\ , but applies to
 communication mode *multi* instead. Since in this case the communication
 cutoffs are determined per atom type, a type specifier is needed and
 cutoff for one or multiple types can be extended. Also ranges of types
-using the usual asterisk notation can be given.
+using the usual asterisk notation can be given. For granular pairstyles, 
+the default cutoff is set to the sum of the current maximum atomic radii
+for each type. 
 
 These are simulation scenarios in which it may be useful or even
 necessary to set a ghost cutoff > neighbor cutoff:

doc/src/neighbor.rst

@@ -49,7 +49,9 @@ sometimes be faster.  Either style should give the same answers.
 
 The *multi* style is a modified binning algorithm that is useful for
 systems with a wide range of cutoff distances, e.g. due to different
-size particles.  For the *bin* style, the bin size is set to 1/2 of
+size particles. For granular pairstyles, cutoffs are set to the 
+sum of the maximum atomic radii for each atom type.
+For the *bin* style, the bin size is set to 1/2 of
 the largest cutoff distance between any pair of atom types and a
 single set of bins is defined to search over for all atom types.  This
 can be inefficient if one pair of types has a very long cutoff, but

src/USER-MISC/dihedral_quadratic.cpp

@@ -195,6 +195,8 @@ void DihedralQuadratic::compute(int eflag, int vflag)
     siinv = 1.0/si;
 
     double dphi = phi-phi0[type];
+    if (dphi > MY_PI) dphi -= 2*MY_PI;
+    else if (dphi < -MY_PI) dphi += 2*MY_PI;
     p = k[type]*dphi;
     pd = - 2.0 * p * siinv;
     p = p * dphi;

src/USER-OMP/dihedral_quadratic_omp.cpp

@@ -23,11 +23,13 @@
 #include "neighbor.h"
 #include "force.h"
 #include "update.h"
+#include "math_const.h"
 #include "error.h"
 
 
 #include "suffix.h"
 using namespace LAMMPS_NS;
+using namespace MathConst;
 
 #define TOLERANCE 0.05
 #define SMALL     0.001
@@ -218,6 +220,8 @@ void DihedralQuadraticOMP::eval(int nfrom, int nto, ThrData * const thr)
     siinv = 1.0/si;
 
     double dphi = phi-phi0[type];
+    if (dphi > MY_PI) dphi -= 2*MY_PI;
+    else if (dphi < -MY_PI) dphi += 2*MY_PI;
     p = k[type]*dphi;
     pd = - 2.0 * p * siinv;
     p = p * dphi;

src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp

@@ -0,0 +1,126 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#include "omp_compat.h"
+#include "npair_half_size_multi_newtoff_omp.h"
+#include "npair_omp.h"
+#include "neigh_list.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "my_page.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+NPairHalfSizeMultiNewtoffOmp::NPairHalfSizeMultiNewtoffOmp(LAMMPS *lmp) :
+  NPair(lmp) {}
+
+/* ----------------------------------------------------------------------
+   size particles
+   binned neighbor list construction with partial Newton's 3rd law
+   each owned atom i checks own bin and other bins in stencil
+   multi-type stencil is itype dependent and is distance checked
+   pair stored once if i,j are both owned and i < j
+   pair stored by me if j is ghost (also stored by proc owning j)
+------------------------------------------------------------------------- */
+
+void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list)
+{
+  const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
+  const int history = list->history;
+  const int mask_history = 3 << SBBITS;
+
+  NPAIR_OMP_INIT;
+#if defined(_OPENMP)
+#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(list)
+#endif
+  NPAIR_OMP_SETUP(nlocal);
+
+  int i,j,k,m,n,itype,jtype,ibin,ns;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  double radi,radsum,cutdistsq;
+  int *neighptr,*s;
+  double *cutsq,*distsq;
+
+  double **x = atom->x;
+  double *radius = atom->radius;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  tagint *molecule = atom->molecule;
+
+  int *ilist = list->ilist;
+  int *numneigh = list->numneigh;
+  int **firstneigh = list->firstneigh;
+
+  // each thread has its own page allocator
+  MyPage<int> &ipage = list->ipage[tid];
+  ipage.reset();
+
+  for (i = ifrom; i < ito; i++) {
+
+    n = 0;
+    neighptr = ipage.vget();
+
+    itype = type[i];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    radi = radius[i];
+
+    // loop over all atoms in other bins in stencil including self
+    // only store pair if i < j
+    // skip if i,j neighbor cutoff is less than bin distance
+    // stores own/own pairs only once
+    // stores own/ghost pairs on both procs
+
+    ibin = atom2bin[i];
+    s = stencil_multi[itype];
+    distsq = distsq_multi[itype];
+    cutsq = cutneighsq[itype];
+    ns = nstencil_multi[itype];
+    for (k = 0; k < ns; k++) {
+      for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
+        if (j <= i) continue;
+        jtype = type[j];
+        if (cutsq[jtype] < distsq[k]) continue;
+
+        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;
+        radsum = radi + radius[j];
+        cutdistsq = (radsum+skin) * (radsum+skin);
+
+        if (rsq <= cutdistsq) {
+          if (history && rsq < radsum*radsum)
+            neighptr[n++] = j ^ mask_history;
+          else
+            neighptr[n++] = j;
+        }
+      }
+    }
+
+    ilist[i] = i;
+    firstneigh[i] = neighptr;
+    numneigh[i] = n;
+    ipage.vgot(n);
+    if (ipage.status())
+      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
+  }
+  NPAIR_OMP_CLOSE;
+  list->inum = nlocal;
+}

src/USER-OMP/npair_half_size_multi_newtoff_omp.h

@@ -0,0 +1,44 @@
+/* -*- 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/size/multi/newtoff/omp,
+           NPairHalfSizeMultiNewtoffOmp,
+           NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTOFF | NP_OMP |
+           NP_ORTHO | NP_TRI)
+
+#else
+
+#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTOFF_OMP_H
+#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTOFF_OMP_H
+
+#include "npair.h"
+
+namespace LAMMPS_NS {
+
+class NPairHalfSizeMultiNewtoffOmp : public NPair {
+ public:
+  NPairHalfSizeMultiNewtoffOmp(class LAMMPS *);
+  ~NPairHalfSizeMultiNewtoffOmp() {}
+  void build(class NeighList *);
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+*/

src/USER-OMP/npair_half_size_multi_newton_omp.cpp

@@ -0,0 +1,151 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#include "omp_compat.h"
+#include "npair_half_size_multi_newton_omp.h"
+#include "npair_omp.h"
+#include "neigh_list.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "my_page.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+NPairHalfSizeMultiNewtonOmp::NPairHalfSizeMultiNewtonOmp(LAMMPS *lmp) :
+  NPair(lmp) {}
+
+/* ----------------------------------------------------------------------
+   size particles
+   binned neighbor list construction with full Newton's 3rd law
+   each owned atom i checks its own bin and other bins in Newton stencil
+   multi-type stencil is itype dependent and is distance checked
+   every pair stored exactly once by some processor
+------------------------------------------------------------------------- */
+
+void NPairHalfSizeMultiNewtonOmp::build(NeighList *list)
+{
+  const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
+  const int history = list->history;
+  const int mask_history = 3 << SBBITS;
+
+  NPAIR_OMP_INIT;
+#if defined(_OPENMP)
+#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(list)
+#endif
+  NPAIR_OMP_SETUP(nlocal);
+
+  int i,j,k,m,n,itype,jtype,ibin,ns;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  double radi,radsum,cutdistsq;
+  int *neighptr,*s;
+  double *cutsq,*distsq;
+
+  double **x = atom->x;
+  double *radius = atom->radius;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  tagint *molecule = atom->molecule;
+
+  int *ilist = list->ilist;
+  int *numneigh = list->numneigh;
+  int **firstneigh = list->firstneigh;
+
+  // each thread has its own page allocator
+  MyPage<int> &ipage = list->ipage[tid];
+  ipage.reset();
+
+  for (i = ifrom; i < ito; i++) {
+
+    n = 0;
+    neighptr = ipage.vget();
+
+    itype = type[i];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    radi = radius[i];
+
+    // loop over rest of atoms in i's bin, ghosts are at end of linked list
+    // if j is owned atom, store it, since j is beyond i in linked list
+    // if j is ghost, only store if j coords are "above and to the right" of i
+
+    for (j = bins[i]; j >= 0; j = bins[j]) {
+      if (j >= nlocal) {
+        if (x[j][2] < ztmp) continue;
+        if (x[j][2] == ztmp) {
+          if (x[j][1] < ytmp) continue;
+          if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
+        }
+      }
+
+      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;
+      radsum = radi + radius[j];
+      cutdistsq = (radsum+skin) * (radsum+skin);
+
+      if (rsq <= cutdistsq) {
+        if (history && rsq < radsum*radsum)
+          neighptr[n++] = j ^ mask_history;
+        else
+          neighptr[n++] = j;
+      }
+    }
+
+    // loop over all atoms in other bins in stencil, store every pair
+    // skip if i,j neighbor cutoff is less than bin distance
+
+    ibin = atom2bin[i];
+    s = stencil_multi[itype];
+    distsq = distsq_multi[itype];
+    cutsq = cutneighsq[itype];
+    ns = nstencil_multi[itype];
+    for (k = 0; k < ns; k++) {
+      for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
+        jtype = type[j];
+        if (cutsq[jtype] < distsq[k]) continue;
+
+        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;
+        radsum = radi + radius[j];
+        cutdistsq = (radsum+skin) * (radsum+skin);
+
+        if (rsq <= cutdistsq) {
+          if (history && rsq < radsum*radsum)
+            neighptr[n++] = j ^ mask_history;
+          else
+            neighptr[n++] = j;
+        }
+      }
+    }
+
+    ilist[i] = i;
+    firstneigh[i] = neighptr;
+    numneigh[i] = n;
+    ipage.vgot(n);
+    if (ipage.status())
+      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
+  }
+  NPAIR_OMP_CLOSE;
+  list->inum = nlocal;
+}

src/USER-OMP/npair_half_size_multi_newton_omp.h

@@ -0,0 +1,43 @@
+/* -*- 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/size/multi/newton/omp,
+           NPairHalfSizeMultiNewtonOmp,
+           NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTON | NP_OMP | NP_ORTHO)
+
+#else
+
+#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_OMP_H
+#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_OMP_H
+
+#include "npair.h"
+
+namespace LAMMPS_NS {
+
+class NPairHalfSizeMultiNewtonOmp : public NPair {
+ public:
+  NPairHalfSizeMultiNewtonOmp(class LAMMPS *);
+  ~NPairHalfSizeMultiNewtonOmp() {}
+  void build(class NeighList *);
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+*/

src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp

@@ -0,0 +1,134 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#include "omp_compat.h"
+#include "npair_half_size_multi_newton_tri_omp.h"
+#include "npair_omp.h"
+#include "neigh_list.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "my_page.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+NPairHalfSizeMultiNewtonTriOmp::NPairHalfSizeMultiNewtonTriOmp(LAMMPS *lmp) :
+  NPair(lmp) {}
+
+/* ----------------------------------------------------------------------
+   size particles
+   binned neighbor list construction with Newton's 3rd law for triclinic
+   each owned atom i checks its own bin and other bins in triclinic stencil
+   multi-type stencil is itype dependent and is distance checked
+   every pair stored exactly once by some processor
+------------------------------------------------------------------------- */
+
+void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list)
+{
+  const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal;
+  const int history = list->history;
+  const int mask_history = 3 << SBBITS;
+
+  NPAIR_OMP_INIT;
+#if defined(_OPENMP)
+#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(list)
+#endif
+  NPAIR_OMP_SETUP(nlocal);
+
+  int i,j,k,m,n,itype,jtype,ibin,ns;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  double radi,radsum,cutdistsq;
+  int *neighptr,*s;
+  double *cutsq,*distsq;
+
+  double **x = atom->x;
+  double *radius = atom->radius;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  tagint *molecule = atom->molecule;
+
+  int *ilist = list->ilist;
+  int *numneigh = list->numneigh;
+  int **firstneigh = list->firstneigh;
+
+  // each thread has its own page allocator
+  MyPage<int> &ipage = list->ipage[tid];
+  ipage.reset();
+
+  for (i = ifrom; i < ito; i++) {
+
+    n = 0;
+    neighptr = ipage.vget();
+
+    itype = type[i];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    radi = radius[i];
+
+    // loop over all atoms in bins, including self, in stencil
+    // skip if i,j neighbor cutoff is less than bin distance
+    // bins below self are excluded from stencil
+    // pairs for atoms j "below" i are excluded
+    // below = lower z or (equal z and lower y) or (equal zy and lower x)
+    //         (equal zyx and j <= i)
+    // latter excludes self-self interaction but allows superposed atoms
+
+    ibin = atom2bin[i];
+    s = stencil_multi[itype];
+    distsq = distsq_multi[itype];
+    cutsq = cutneighsq[itype];
+    ns = nstencil_multi[itype];
+    for (k = 0; k < ns; k++) {
+      for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
+        jtype = type[j];
+        if (cutsq[jtype] < distsq[k]) continue;
+        if (x[j][2] < ztmp) continue;
+        if (x[j][2] == ztmp) {
+          if (x[j][1] < ytmp) continue;
+          if (x[j][1] == ytmp) {
+            if (x[j][0] < xtmp) continue;
+            if (x[j][0] == xtmp && j <= i) continue;
+          }
+        }
+
+        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;
+        radsum = radi + radius[j];
+        cutdistsq = (radsum+skin) * (radsum+skin);
+
+        if (rsq <= cutdistsq) {
+          if (history && rsq < radsum*radsum)
+            neighptr[n++] = j ^ mask_history;
+          else
+            neighptr[n++] = j;
+        }
+      }
+    }
+
+    ilist[i] = i;
+    firstneigh[i] = neighptr;
+    numneigh[i] = n;
+    ipage.vgot(n);
+    if (ipage.status())
+      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
+  }
+  NPAIR_OMP_CLOSE;
+  list->inum = nlocal;
+}

src/USER-OMP/npair_half_size_multi_newton_tri_omp.h

@@ -0,0 +1,43 @@
+/* -*- 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/size/multi/newton/tri/omp,
+           NPairHalfSizeMultiNewtonTriOmp,
+           NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTON | NP_TRI | NP_OMP)
+
+#else
+
+#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_TRI_OMP_H
+#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_TRI_OMP_H
+
+#include "npair.h"
+
+namespace LAMMPS_NS {
+
+class NPairHalfSizeMultiNewtonTriOmp : public NPair {
+ public:
+  NPairHalfSizeMultiNewtonTriOmp(class LAMMPS *);
+  ~NPairHalfSizeMultiNewtonTriOmp() {}
+  void build(class NeighList *);
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+*/

src/npair_half_size_multi_newtoff.cpp

@@ -0,0 +1,117 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#include "npair_half_size_multi_newtoff.h"
+#include "neigh_list.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "my_page.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+NPairHalfSizeMultiNewtoff::NPairHalfSizeMultiNewtoff(LAMMPS *lmp) : NPair(lmp) {}
+
+/* ----------------------------------------------------------------------
+   size particles
+   binned neighbor list construction with partial Newton's 3rd law
+   each owned atom i checks own bin and other bins in stencil
+   multi-type stencil is itype dependent and is distance checked
+   pair stored once if i,j are both owned and i < j
+   pair stored by me if j is ghost (also stored by proc owning j)
+------------------------------------------------------------------------- */
+
+void NPairHalfSizeMultiNewtoff::build(NeighList *list)
+{
+  int i,j,k,m,n,itype,jtype,ibin,ns;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  double radi,radsum,cutdistsq;
+  int *neighptr,*s;
+  double *cutsq,*distsq;
+
+  double **x = atom->x;
+  double *radius = atom->radius;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  tagint *molecule = atom->molecule;
+  int nlocal = atom->nlocal;
+  if (includegroup) nlocal = atom->nfirst;
+
+  int history = list->history;
+  int *ilist = list->ilist;
+  int *numneigh = list->numneigh;
+  int **firstneigh = list->firstneigh;
+  MyPage<int> *ipage = list->ipage;
+
+  int mask_history = 3 << SBBITS;
+
+  int inum = 0;
+  ipage->reset();
+
+  for (i = 0; i < nlocal; i++) {
+    n = 0;
+    neighptr = ipage->vget();
+
+    itype = type[i];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    radi = radius[i];
+
+    // loop over all atoms in other bins in stencil including self
+    // only store pair if i < j
+    // skip if i,j neighbor cutoff is less than bin distance
+    // stores own/own pairs only once
+    // stores own/ghost pairs on both procs
+
+    ibin = atom2bin[i];
+    s = stencil_multi[itype];
+    distsq = distsq_multi[itype];
+    cutsq = cutneighsq[itype];
+    ns = nstencil_multi[itype];
+    for (k = 0; k < ns; k++) {
+      for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
+        if (j <= i) continue;
+        jtype = type[j];
+        if (cutsq[jtype] < distsq[k]) continue;
+
+        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;
+        radsum = radi + radius[j];
+        cutdistsq = (radsum+skin) * (radsum+skin);
+
+        if (rsq <= cutdistsq) {
+          if (history && rsq < radsum*radsum)
+            neighptr[n++] = j ^ mask_history;
+          else
+            neighptr[n++] = j;
+        }
+      }
+    }
+
+    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");
+  }
+
+  list->inum = inum;
+}

src/npair_half_size_multi_newtoff.h

@@ -0,0 +1,46 @@
+/* -*- 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/size/multi/newtoff,
+           NPairHalfSizeMultiNewtoff,
+           NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTOFF | NP_ORTHO | NP_TRI)
+
+#else
+
+#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTOFF_H
+#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTOFF_H
+
+#include "npair.h"
+
+namespace LAMMPS_NS {
+
+class NPairHalfSizeMultiNewtoff : public NPair {
+ public:
+  NPairHalfSizeMultiNewtoff(class LAMMPS *);
+  ~NPairHalfSizeMultiNewtoff() {}
+  void build(class NeighList *);
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Neighbor list overflow, boost neigh_modify one
+
+UNDOCUMENTED
+*/

src/npair_half_size_multi_newton.cpp

@@ -0,0 +1,143 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#include "npair_half_size_multi_newton.h"
+#include "neigh_list.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "my_page.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+NPairHalfSizeMultiNewton::NPairHalfSizeMultiNewton(LAMMPS *lmp) : NPair(lmp) {}
+
+/* ----------------------------------------------------------------------
+   size particles
+   binned neighbor list construction with full Newton's 3rd law
+   each owned atom i checks its own bin and other bins in Newton stencil
+   multi-type stencil is itype dependent and is distance checked
+   every pair stored exactly once by some processor
+------------------------------------------------------------------------- */
+
+void NPairHalfSizeMultiNewton::build(NeighList *list)
+{
+  int i,j,k,m,n,itype,jtype,ibin,ns;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  double radi,radsum,cutdistsq;
+  int *neighptr,*s;
+  double *cutsq,*distsq;
+
+  double **x = atom->x;
+  double *radius = atom->radius;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  tagint *molecule = atom->molecule;
+  int nlocal = atom->nlocal;
+  if (includegroup) nlocal = atom->nfirst;
+
+  int history = list->history;
+  int *ilist = list->ilist;
+  int *numneigh = list->numneigh;
+  int **firstneigh = list->firstneigh;
+  MyPage<int> *ipage = list->ipage;
+
+  int mask_history = 3 << SBBITS;
+
+  int inum = 0;
+  ipage->reset();
+
+  for (i = 0; i < nlocal; i++) {
+    n = 0;
+    neighptr = ipage->vget();
+
+    itype = type[i];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    radi = radius[i];
+
+    // loop over rest of atoms in i's bin, ghosts are at end of linked list
+    // if j is owned atom, store it, since j is beyond i in linked list
+    // if j is ghost, only store if j coords are "above and to the right" of i
+
+    for (j = bins[i]; j >= 0; j = bins[j]) {
+      if (j >= nlocal) {
+        if (x[j][2] < ztmp) continue;
+        if (x[j][2] == ztmp) {
+          if (x[j][1] < ytmp) continue;
+          if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
+        }
+      }
+
+      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;
+      radsum = radi + radius[j];
+      cutdistsq = (radsum+skin) * (radsum+skin);
+
+      if (rsq <= cutdistsq) {
+        if (history && rsq < radsum*radsum) 
+          neighptr[n++] = j ^ mask_history;
+        else 
+          neighptr[n++] = j;
+      }
+    }
+
+    // loop over all atoms in other bins in stencil, store every pair
+    // skip if i,j neighbor cutoff is less than bin distance
+
+    ibin = atom2bin[i];
+    s = stencil_multi[itype];
+    distsq = distsq_multi[itype];
+    cutsq = cutneighsq[itype];
+    ns = nstencil_multi[itype];
+    for (k = 0; k < ns; k++) {
+      for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
+        jtype = type[j];
+        if (cutsq[jtype] < distsq[k]) continue;
+
+        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;
+        radsum = radi + radius[j];
+        cutdistsq = (radsum+skin) * (radsum+skin);
+
+        if (rsq <= cutdistsq) {
+          if (history && rsq < radsum*radsum) 
+            neighptr[n++] = j ^ mask_history;
+          else
+            neighptr[n++] = j;
+        }
+      }
+    }
+
+    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");
+  }
+
+  list->inum = inum;
+}

src/npair_half_size_multi_newton.h

@@ -0,0 +1,46 @@
+/* -*- 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/size/multi/newton,
+           NPairHalfSizeMultiNewton,
+           NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTON | NP_ORTHO)
+
+#else
+
+#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_H
+#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_H
+
+#include "npair.h"
+
+namespace LAMMPS_NS {
+
+class NPairHalfSizeMultiNewton : public NPair {
+ public:
+  NPairHalfSizeMultiNewton(class LAMMPS *);
+  ~NPairHalfSizeMultiNewton() {}
+  void build(class NeighList *);
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Neighbor list overflow, boost neigh_modify one
+
+UNDOCUMENTED
+*/

src/npair_half_size_multi_newton_tri.cpp

@@ -0,0 +1,125 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#include "npair_half_size_multi_newton_tri.h"
+#include "neigh_list.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "my_page.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+NPairHalfSizeMultiNewtonTri::NPairHalfSizeMultiNewtonTri(LAMMPS *lmp) : NPair(lmp) {}
+
+/* ----------------------------------------------------------------------
+   binned neighbor list construction with Newton's 3rd law for triclinic
+   each owned atom i checks its own bin and other bins in triclinic stencil
+   multi-type stencil is itype dependent and is distance checked
+   every pair stored exactly once by some processor
+------------------------------------------------------------------------- */
+
+void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
+{
+  int i,j,k,m,n,itype,jtype,ibin,ns;
+  double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
+  double radi,radsum,cutdistsq;
+  int *neighptr,*s;
+  double *cutsq,*distsq;
+
+  double **x = atom->x;
+  double *radius = atom->radius;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  tagint *molecule = atom->molecule;
+  int nlocal = atom->nlocal;
+  if (includegroup) nlocal = atom->nfirst;
+
+  int history = list->history;
+  int *ilist = list->ilist;
+  int *numneigh = list->numneigh;
+  int **firstneigh = list->firstneigh;
+  MyPage<int> *ipage = list->ipage;
+
+  int mask_history = 3 << SBBITS;
+
+  int inum = 0;
+  ipage->reset();
+
+  for (i = 0; i < nlocal; i++) {
+    n = 0;
+    neighptr = ipage->vget();
+
+    itype = type[i];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    radi = radius[i];
+
+
+    // loop over all atoms in bins, including self, in stencil
+    // skip if i,j neighbor cutoff is less than bin distance
+    // bins below self are excluded from stencil
+    // pairs for atoms j "below" i are excluded
+    // below = lower z or (equal z and lower y) or (equal zy and lower x)
+    //         (equal zyx and j <= i)
+    // latter excludes self-self interaction but allows superposed atoms
+
+    ibin = atom2bin[i];
+    s = stencil_multi[itype];
+    distsq = distsq_multi[itype];
+    cutsq = cutneighsq[itype];
+    ns = nstencil_multi[itype];
+    for (k = 0; k < ns; k++) {
+      for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
+        jtype = type[j];
+        if (cutsq[jtype] < distsq[k]) continue;
+        if (x[j][2] < ztmp) continue;
+        if (x[j][2] == ztmp) {
+          if (x[j][1] < ytmp) continue;
+          if (x[j][1] == ytmp) {
+            if (x[j][0] < xtmp) continue;
+            if (x[j][0] == xtmp && j <= i) continue;
+          }
+        }
+
+        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;
+        radsum = radi + radius[j];
+        cutdistsq = (radsum+skin) * (radsum+skin);
+
+        if (rsq <= cutdistsq) {
+          if (history && rsq < radsum*radsum) 
+            neighptr[n++] = j ^ mask_history;
+          else
+            neighptr[n++] = j;
+        }
+      }
+    }
+
+    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");
+  }
+
+  list->inum = inum;
+}

src/npair_half_size_multi_newton_tri.h

@@ -0,0 +1,46 @@
+/* -*- 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/size/multi/newton/tri,
+           NPairHalfSizeMultiNewtonTri,
+           NP_HALF | NP_SIZE | NP_MULTI | NP_NEWTON | NP_TRI)
+
+#else
+
+#ifndef LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_TRI_H
+#define LMP_NPAIR_HALF_SIZE_MULTI_NEWTON_TRI_H
+
+#include "npair.h"
+
+namespace LAMMPS_NS {
+
+class NPairHalfSizeMultiNewtonTri : public NPair {
+ public:
+  NPairHalfSizeMultiNewtonTri(class LAMMPS *);
+  ~NPairHalfSizeMultiNewtonTri() {}
+  void build(class NeighList *);
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Neighbor list overflow, boost neigh_modify one
+
+UNDOCUMENTED
+*/