More edits and documentation

doc/src/comm_modify.rst

@@ -11,7 +11,7 @@ Syntax
    comm_modify keyword value ...
 
 * zero or more keyword/value pairs may be appended
-* keyword = *mode* or *cutoff* or *cutoff/multi* or *group* or *vel*
+* keyword = *mode* or *cutoff* or *cutoff/multi* or *cutoff/bytype* or *group* or *vel*
 
   .. parsed-literal::
 
@@ -20,6 +20,7 @@ Syntax
        *cutoff/multi* type value
           type = atom type or type range (supports asterisk notation)
           value = Rcut (distance units) = communicate atoms for selected types from this far away
+       *cutoff/byptype* arg = none = communicate atoms for each type by a distance equal to the largest interaction distance for that type
        *group* value = group-ID = only communicate atoms in the group
        *vel* value = *yes* or *no* = do or do not communicate velocity info with ghost atoms
 
@@ -92,6 +93,12 @@ 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.
 
+The *cutoff/bytype* option applies to *multi* and automtically sets communication 
+cutoffs for each particle type based on the largest interaction distance
+between two particles of the same type. This method is only compatible 
+with Newton on and the *bytype* neighbor style. See the :doc:`neighbor bytype <neighbor>` 
+command for more information.
+
 These are simulation scenarios in which it may be useful or even
 necessary to set a ghost cutoff > neighbor cutoff:
 

doc/src/neighbor.rst

@@ -11,7 +11,7 @@ Syntax
    neighbor skin style
 
 * skin = extra distance beyond force cutoff (distance units)
-* style = *bin* or *nsq* or *multi*
+* style = *bin* or *nsq* or *multi* or *bytype*
 
 Examples
 """"""""
@@ -60,6 +60,14 @@ This imposes some extra setup overhead, but the searches themselves
 may be much faster for the short-cutoff cases.  See the :doc:`comm_modify mode multi <comm_modify>` command for a communication option
 that may also be beneficial for simulations of this kind.
 
+The *bytype* style is an extension of the *multi* style that was 
+presented by Shire, Hanley, and Stratford :ref:`(Shire) <bytype-Shire>`. 
+For style *bytype*, different bin lists are created for each different 
+type and separate bin sizes are generated. Whether *bytype* or *multi* 
+is faster may depend on the specifics of your system. See the 
+:doc:`comm_modify mode bytype <comm_modify>` command for a compatible 
+communication option.
+
 The :doc:`neigh_modify <neigh_modify>` command has additional options
 that control how often neighbor lists are built and which pairs are
 stored in the list.
@@ -86,3 +94,9 @@ Default
 | 0.001 bin for units = si, skin = 0.001 meters = 1.0 mm
 | 0.1 bin for units = cgs, skin = 0.1 cm = 1.0 mm
 |
+
+----------
+
+.. _bytype-Shire:
+
+**(Shire)** Shire, Hanley and Stratford, Comp Part Mech, (2020).

src/comm.h

@@ -199,6 +199,10 @@ E: Use cutoff keyword to set cutoff in single mode
 
 Mode is single so cutoff/multi keyword cannot be used.
 
+E: Use cutoff/bytype in mode multi only
+
+Mode is single so cutoff/bytype keyword cannot be used.
+
 E: Cannot set cutoff/multi before simulation box is defined
 
 Self-explanatory.

src/comm_tiled.cpp

@@ -183,6 +183,12 @@ void CommTiled::setup()
       cutghostmulti[i][0] = MAX(cut,cuttype[i]);
       cutghostmulti[i][1] = MAX(cut,cuttype[i]);
       cutghostmulti[i][2] = MAX(cut,cuttype[i]);
+      if (multi_bytype == 1) {
+        // Set the BYTYPE cutoff
+        cutghostmulti[i][0] = sqrt(neighbor->cutneighsq[i][i]);
+        cutghostmulti[i][1] = sqrt(neighbor->cutneighsq[i][i]);
+        cutghostmulti[i][2] = sqrt(neighbor->cutneighsq[i][i]);
+      }      
     }
   }
 

src/nbin_bytype.cpp

@@ -19,14 +19,12 @@
 #include "comm.h"
 #include "update.h"
 #include "error.h"
-
-#include <assert.h>
 #include "memory.h"
 
 using namespace LAMMPS_NS;
 
-#define SMALL 1.0e-6         // Duplicated from NBinStandard
+#define SMALL 1.0e-6
-#define CUT2BIN_RATIO 100    // Duplicated from NBinStandard
+#define CUT2BIN_RATIO 100
 
 /* ---------------------------------------------------------------------- */
 
@@ -238,7 +236,7 @@ void NBinBytype::setup_bins(int style)
     // test for too many global bins in any dimension due to huge global domain
 
     if (bbox[0]*binsizeinv > MAXSMALLINT || bbox[1]*binsizeinv > MAXSMALLINT ||
-	bbox[2]*binsizeinv > MAXSMALLINT)
+	    bbox[2]*binsizeinv > MAXSMALLINT)
       error->all(FLERR,"Domain too large for neighbor bins");
 
     // create actual bins
@@ -269,8 +267,8 @@ void NBinBytype::setup_bins(int style)
     bininvz_type[n] = 1.0 / binsizez_type[n];
 
     if (binsize_optimal*bininvx_type[n] > CUT2BIN_RATIO ||
-	binsize_optimal*bininvy_type[n] > CUT2BIN_RATIO ||
+	    binsize_optimal*bininvy_type[n] > CUT2BIN_RATIO ||
-	binsize_optimal*bininvz_type[n] > CUT2BIN_RATIO)
+	    binsize_optimal*bininvz_type[n] > CUT2BIN_RATIO)
       error->all(FLERR,"Cannot use neighbor bins - box size << cutoff");
 
     // mbinlo/hi = lowest and highest global bins my ghost atoms could be in
@@ -352,8 +350,8 @@ void NBinBytype::bin_atoms()
     int bitmask = group->bitmask[includegroup];
     for (i = nall-1; i >= nlocal; i--) {
       if (mask[i] & bitmask) {
-	n = type[i];
+        n = type[i];
-	ibin = coord2bin(x[i], n);
+        ibin = coord2bin(x[i], n);
         atom2bin_type[n][i] = ibin;
         bins_type[n][i] = binhead_type[n][ibin];
         binhead_type[n][ibin] = i;

src/nbin_bytype.h

@@ -53,4 +53,20 @@ class NBinBytype : public NBin {
 
 /* ERROR/WARNING messages:
 
+E: Domain too large for neighbor bins
+
+UNDOCUMENTED
+
+E: Cannot use neighbor bins - box size << cutoff
+
+UNDOCUMENTED
+
+E: Too many neighbor bins
+
+UNDOCUMENTED
+
+E Non-numeric positions - simulation unstable
+
+UNDOCUMENTED
+
 */

src/npair_full_bytype.cpp

@@ -12,7 +12,6 @@
 ------------------------------------------------------------------------- */
 
 #include "npair_full_bytype.h"
-#include "neighbor.h"
 #include "neigh_list.h"
 #include "atom.h"
 #include "atom_vec.h"
@@ -20,7 +19,6 @@
 #include "domain.h"
 #include "my_page.h"
 #include "error.h"
-
 #include "nbin.h"
 #include "nstencil.h"
 
@@ -102,9 +100,9 @@ void NPairFullBytype::build(NeighList *list)
       for (k = 0; k < ns; k++) {
 	int js = this->nb->binhead_type[ktype][kbin + s[k]];
 	for (j = js; j >= 0; j = this->nb->bins_type[ktype][j]) {
-	  jtype = type[j];
 	  if (i == j) continue;
-
+	  
+      jtype = type[j];
 	  if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
 
 	  delx = xtmp - x[j][0];

src/npair_full_bytype.h

@@ -40,4 +40,8 @@ class NPairFullBytype : public NPair {
 
 /* ERROR/WARNING messages:
 
+E: Neighbor list overflow, boost neigh_modify one
+
+UNDOCUMENTED
+
 */

src/npair_half_size_bytype_newtoff.cpp

@@ -13,13 +13,11 @@ es   certain rights in this software.  This software is distributed under
 
 #include <string.h>
 #include "npair_half_size_bytype_newtoff.h"
-#include "neighbor.h"
 #include "neigh_list.h"
 #include "atom.h"
 #include "atom_vec.h"
 #include "my_page.h"
 #include "error.h"
-
 #include "nbin.h"
 #include "nstencil.h"
 

src/npair_half_size_bytype_newtoff.h

@@ -40,4 +40,8 @@ class NPairHalfSizeBytypeNewtoff : public NPair {
 
 /* ERROR/WARNING messages:
 
+E: Neighbor list overflow, boost neigh_modify one
+
+UNDOCUMENTED
+
 */

src/npair_half_size_bytype_newton.cpp

@@ -13,13 +13,11 @@
 
 #include <string.h>
 #include "npair_half_size_bytype_newton.h"
-#include "neighbor.h"
 #include "neigh_list.h"
 #include "atom.h"
 #include "atom_vec.h"
 #include "my_page.h"
 #include "error.h"
-
 #include "nbin.h"
 #include "nstencil.h"
 

src/npair_half_size_bytype_newton.h

@@ -40,4 +40,8 @@ class NPairHalfSizeBytypeNewton : public NPair {
 
 /* ERROR/WARNING messages:
 
-*/
+E: Neighbor list overflow, boost neigh_modify one
+
+UNDOCUMENTED
+
+*/

src/nstencil_full_bytype_3d.cpp

@@ -17,7 +17,6 @@
 #include "nbin.h"
 #include "memory.h"
 #include "atom.h"
-
 #include <math.h>
 
 using namespace LAMMPS_NS;
@@ -30,6 +29,8 @@ NStencilFullBytype3d::NStencilFullBytype3d(LAMMPS *lmp) :
   maxstencil_type = NULL;
 }
 
+/* ---------------------------------------------------------------------- */
+
 NStencilFullBytype3d::~NStencilFullBytype3d() {
 
   if (maxstencil_type) {
@@ -45,6 +46,8 @@ NStencilFullBytype3d::~NStencilFullBytype3d() {
   }
 }
 
+/* ---------------------------------------------------------------------- */
+
 void NStencilFullBytype3d::copy_bin_info_bytype(int itype) {
 
   mbinx = nb->mbinx_type[itype];
@@ -58,6 +61,8 @@ void NStencilFullBytype3d::copy_bin_info_bytype(int itype) {
   bininvz = nb->bininvz_type[itype];
 }
 
+/* ---------------------------------------------------------------------- */
+
 int NStencilFullBytype3d::copy_neigh_info_bytype(int itype) {
 
   cutneighmaxsq = neighbor->cutneighsq[itype][itype];
@@ -78,6 +83,8 @@ int NStencilFullBytype3d::copy_neigh_info_bytype(int itype) {
   return ((2*sx+1) * (2*sy+1) * (2*sz+1));
 }
 
+/* ---------------------------------------------------------------------- */
+
 void NStencilFullBytype3d::create_setup()
 {
 
@@ -85,7 +92,6 @@ void NStencilFullBytype3d::create_setup()
   int maxtypes;
   int smax;
 
-  //printf("NStencilFullBytype3d::create_steup()\n");
   maxtypes = atom->ntypes;
 
   if (maxstencil_type == NULL) {
@@ -159,6 +165,8 @@ void NStencilFullBytype3d::create_setup()
 
 }
 
+/* ---------------------------------------------------------------------- */
+
 void NStencilFullBytype3d::create_newtoff(int itype, int jtype, double cutsq) {
 
   int i, j, k, ns;

src/nstencil_half_bytype_2d_newton.cpp

@@ -17,7 +17,6 @@
 #include "nbin.h"
 #include "memory.h"
 #include "atom.h"
-
 #include <math.h>
 
 using namespace LAMMPS_NS;
@@ -43,6 +42,8 @@ NStencilHalfBytype2dNewton::~NStencilHalfBytype2dNewton() {
   memory->destroy(maxstencil_type);
 }
 
+/* ---------------------------------------------------------------------- */
+
 // KS To superclass
 void NStencilHalfBytype2dNewton::copy_bin_info_bytype(int itype) {
 
@@ -57,6 +58,8 @@ void NStencilHalfBytype2dNewton::copy_bin_info_bytype(int itype) {
   bininvz = nb->bininvz_type[itype];
 }
 
+/* ---------------------------------------------------------------------- */
+
 // KS To superclass?
 int NStencilHalfBytype2dNewton::copy_neigh_info_bytype(int itype) {
 
@@ -78,6 +81,8 @@ int NStencilHalfBytype2dNewton::copy_neigh_info_bytype(int itype) {
   return ((2*sx+1) * (2*sy+1) * (2*sz+1));
 }
 
+/* ---------------------------------------------------------------------- */
+
 void NStencilHalfBytype2dNewton::create_setup()
 {
 
@@ -153,6 +158,8 @@ void NStencilHalfBytype2dNewton::create_setup()
 
 }
 
+/* ---------------------------------------------------------------------- */
+
 void NStencilHalfBytype2dNewton::create_newton(int itype, int jtype, double cutsq) {
 
   int i, j, ns;
@@ -168,6 +175,8 @@ void NStencilHalfBytype2dNewton::create_newton(int itype, int jtype, double cuts
   nstencil_type[itype][jtype] = ns;
 }
 
+/* ---------------------------------------------------------------------- */
+
 void NStencilHalfBytype2dNewton::create_newtoff(int itype, int jtype, double cutsq) {
 
   int i, j, ns;

src/nstencil_half_bytype_3d_newtoff.cpp

@@ -17,7 +17,6 @@
 #include "nbin.h"
 #include "memory.h"
 #include "atom.h"
-
 #include <math.h>
 
 using namespace LAMMPS_NS;
@@ -30,6 +29,8 @@ NStencilHalfBytype3dNewtoff::NStencilHalfBytype3dNewtoff(LAMMPS *lmp) :
   maxstencil_type = NULL;
 }
 
+/* ---------------------------------------------------------------------- */
+
 NStencilHalfBytype3dNewtoff::~NStencilHalfBytype3dNewtoff() {
 
   memory->destroy(nstencil_type);
@@ -43,6 +44,8 @@ NStencilHalfBytype3dNewtoff::~NStencilHalfBytype3dNewtoff() {
   memory->destroy(maxstencil_type);
 }
 
+/* ---------------------------------------------------------------------- */
+
 void NStencilHalfBytype3dNewtoff::copy_bin_info_bytype(int itype) {
 
   mbinx = nb->mbinx_type[itype];
@@ -56,6 +59,8 @@ void NStencilHalfBytype3dNewtoff::copy_bin_info_bytype(int itype) {
   bininvz = nb->bininvz_type[itype];
 }
 
+/* ---------------------------------------------------------------------- */
+
 int NStencilHalfBytype3dNewtoff::copy_neigh_info_bytype(int itype) {
 
   cutneighmaxsq = neighbor->cutneighsq[itype][itype];
@@ -76,6 +81,8 @@ int NStencilHalfBytype3dNewtoff::copy_neigh_info_bytype(int itype) {
   return ((2*sx+1) * (2*sy+1) * (2*sz+1));
 }
 
+/* ---------------------------------------------------------------------- */
+
 void NStencilHalfBytype3dNewtoff::create_setup()
 {
 
@@ -83,7 +90,6 @@ void NStencilHalfBytype3dNewtoff::create_setup()
   int maxtypes;
   int smax;
 
-  //printf("NStencilHalfBytype3dNewtoff::create_steup()\n");
   maxtypes = atom->ntypes;
 
   if (maxstencil_type == NULL) {
@@ -157,6 +163,8 @@ void NStencilHalfBytype3dNewtoff::create_setup()
 
 }
 
+/* ---------------------------------------------------------------------- */
+
 void NStencilHalfBytype3dNewtoff::create_newtoff(int itype, int jtype, double cutsq) {
 
   int i, j, k, ns;

src/nstencil_half_bytype_3d_newton.cpp

@@ -17,7 +17,6 @@
 #include "nbin.h"
 #include "memory.h"
 #include "atom.h"
-
 #include <math.h>
 
 using namespace LAMMPS_NS;
@@ -30,6 +29,8 @@ NStencilHalfBytype3dNewton::NStencilHalfBytype3dNewton(LAMMPS *lmp) :
   maxstencil_type = NULL;
 }
 
+/* ---------------------------------------------------------------------- */
+
 NStencilHalfBytype3dNewton::~NStencilHalfBytype3dNewton() {
 
   memory->destroy(nstencil_type);
@@ -43,6 +44,8 @@ NStencilHalfBytype3dNewton::~NStencilHalfBytype3dNewton() {
   memory->destroy(maxstencil_type);
 }
 
+/* ---------------------------------------------------------------------- */
+
 // KS To superclass
 void NStencilHalfBytype3dNewton::copy_bin_info_bytype(int itype) {
 
@@ -57,6 +60,8 @@ void NStencilHalfBytype3dNewton::copy_bin_info_bytype(int itype) {
   bininvz = nb->bininvz_type[itype];
 }
 
+/* ---------------------------------------------------------------------- */
+
 // KS To superclass?
 int NStencilHalfBytype3dNewton::copy_neigh_info_bytype(int itype) {
 
@@ -78,6 +83,8 @@ int NStencilHalfBytype3dNewton::copy_neigh_info_bytype(int itype) {
   return ((2*sx+1) * (2*sy+1) * (2*sz+1));
 }
 
+/* ---------------------------------------------------------------------- */
+
 void NStencilHalfBytype3dNewton::create_setup()
 {
 
@@ -153,6 +160,8 @@ void NStencilHalfBytype3dNewton::create_setup()
 
 }
 
+/* ---------------------------------------------------------------------- */
+
 void NStencilHalfBytype3dNewton::create_newton(int itype, int jtype, double cutsq) {
 
   int i, j, k, ns;
@@ -169,6 +178,8 @@ void NStencilHalfBytype3dNewton::create_newton(int itype, int jtype, double cuts
   nstencil_type[itype][jtype] = ns;
 }
 
+/* ---------------------------------------------------------------------- */
+
 void NStencilHalfBytype3dNewton::create_newtoff(int itype, int jtype, double cutsq) {
 
   int i, j, k, ns;