Merge pull request #3696 from jtclemm/BPM

Minor updates/patches for BPM, rigid, and multi

doc/src/Modify_gran_sub_mod.rst

@@ -25,7 +25,7 @@ specific implementation. For instance, from the GranSubModNormal class the
 GranSubModNormalHooke, GranSubModNormalHertz, and GranSubModNormalJKR classes
 are derived which calculate Hookean, Hertzian, or JKR normal forces,
 respectively.  This modular structure simplifies the addition of new granular
-contact models as as one only needs to create a new GranSubMod class without
+contact models as one only needs to create a new GranSubMod class without
 having to modify the more complex PairGranular, FixGranWall, and GranularModel
 classes. Most GranSubMod methods are also already defined by the parent classes
 so new contact models typically only require edits to a few relevant methods

doc/src/bond_bpm_rotational.rst

@@ -10,7 +10,7 @@ Syntax
 
    bond_style bpm/rotational keyword value attribute1 attribute2 ...
 
-* optional keyword = *overlay/pair* or *store/local* or *smooth*
+* optional keyword = *overlay/pair* or *store/local* or *smooth* or *break/no*
 
   .. parsed-literal::
 
@@ -30,6 +30,9 @@ Syntax
        *smooth* value = *yes* or *no*
           smooths bond forces near the breaking point
 
+       *break/no*
+          indicates that bonds should not break during a run
+
 Examples
 """"""""
 
@@ -140,6 +143,12 @@ the *overlay/pair* keyword. These settings require specific
 restrictions.  Further details can be found in the `:doc: how to
 <Howto_BPM>` page on BPMs.
 
+.. versionadded:: TBD
+
+If the *break/no* keyword is used, then LAMMPS assumes bonds should not break
+during a simulation run. This will prevent some unnecessary calculation.
+However, if a bond does break, it will trigger an error.
+
 If the *store/local* keyword is used, an internal fix will track bonds that
 break during the simulation. Whenever a bond breaks, data is processed
 and transferred to an internal fix labeled *fix_ID*. This allows the

doc/src/bond_bpm_spring.rst

@@ -10,7 +10,7 @@ Syntax
 
    bond_style bpm/spring keyword value attribute1 attribute2 ...
 
-* optional keyword = *overlay/pair* or *store/local* or *smooth*
+* optional keyword = *overlay/pair* or *store/local* or *smooth* or *break/no*
 
   .. parsed-literal::
 
@@ -30,6 +30,9 @@ Syntax
        *smooth* value = *yes* or *no*
           smooths bond forces near the breaking point
 
+       *break/no*
+          indicates that bonds should not break during a run
+
 Examples
 """"""""
 
@@ -47,7 +50,7 @@ Description
 
 .. versionadded:: 4May2022
 
-The *bpm/spring* bond style computes forces and torques based on
+The *bpm/spring* bond style computes forces based on
 deviations from the initial reference state of the two atoms.  The
 reference state is stored by each bond when it is first computed in
 the setup of a run. Data is then preserved across run commands and is
@@ -56,7 +59,8 @@ the system will not reset the reference state of a bond.
 
 This bond style only applies central-body forces which conserve the
 translational and rotational degrees of freedom of a bonded set of
-particles. The force has a magnitude of
+particles based on a model described by Clemmer and Robbins
+:ref:`(Clemmer) <fragment-Clemmer>`. The force has a magnitude of
 
 .. math::
 
@@ -105,6 +109,12 @@ the *overlay/pair* keyword. These settings require specific
 restrictions.  Further details can be found in the `:doc: how to
 <Howto_BPM>` page on BPMs.
 
+.. versionadded:: TBD
+
+If the *break/no* keyword is used, then LAMMPS assumes bonds should not break
+during a simulation run. This will prevent some unnecessary calculation.
+However, if a bond does break, it will trigger an error.
+
 If the *store/local* keyword is used, an internal fix will track bonds that
 break during the simulation. Whenever a bond breaks, data is processed
 and transferred to an internal fix labeled *fix_ID*. This allows the
@@ -200,6 +210,10 @@ The option defaults are *smooth* = *yes*
 
 ----------
 
+.. _fragment-Clemmer:
+
+**(Clemmer)** Clemmer and Robbins, Phys. Rev. Lett. (2022).
+
 .. _Groot4:
 
 **(Groot)** Groot and Warren, J Chem Phys, 107, 4423-35 (1997).

doc/src/comm_modify.rst

@@ -69,6 +69,7 @@ For many systems this is an efficient algorithm, but for systems with
 widely varying cutoffs for different type pairs, the *multi* or *multi/old* mode can
 be faster.  In *multi*, each atom is assigned to a collection which should
 correspond to a set of atoms with similar interaction cutoffs.
+See the :doc:`neighbor <neighbor>` command for a detailed description of collections.
 In this case, each atom collection is assigned its own distance
 cutoff for communication purposes, and fewer atoms will be
 communicated. in *multi/old*, a similar technique is used but atoms

doc/src/neighbor.rst

@@ -59,9 +59,21 @@ long cutoff, but other type pairs have a much shorter cutoff. The
 sized particles, where "size" may mean the physical size of the particle
 or its cutoff distance for interacting with other particles. Different
 sets of bins are then used to construct the neighbor lists as as further
-described by Shire, Hanley, and Stratford :ref:`(Shire) <bytype-Shire>`.
-This imposes some extra setup overhead, but the searches themselves may
-be much faster. By default, each atom type defines a separate collection
+described by Shire, Hanley, and Stratford :ref:`(Shire) <multi-Shire>`
+and Monti et al. :ref:`(Monti) <multi-Monti>`. This imposes some extra
+setup overhead, but the searches themselves may be much faster.
+
+For instance in a dense binary system in d-dimensions with a ratio of the size
+of the largest to smallest collection bin :math:`\lambda`, the computational
+costs of building a default neighbor list grows as :math:`\lambda^{2d}` while
+the costs for *multi* grows as :math:`\lambda^d`, equivalent to the cost
+of force evaluations, as argued in Monti et al. :ref:`(Monti) <multi-Monti>`.
+In other words, the neighboring costs of *multi* are expected to scale the
+same as force calculations, such that its relative cost is independent of
+the particle size ratio. This is not the case for the default style which
+becomes substantially more expensive with increasing size ratios.
+
+By default in *multi*, each atom type defines a separate collection
 of particles. For systems where two or more atom types have the same
 size (either physical size or cutoff distance), the definition of
 collections can be customized, which can result in less overhead and
@@ -75,8 +87,11 @@ An alternate style, *multi/old*, sets the bin size to 1/2 of the shortest
 cutoff distance and multiple sets of bins are defined to search over for
 different atom types. This algorithm used to be the default *multi*
 algorithm in LAMMPS but was found to be significantly slower than the new
-approach. For now we are keeping the old option in case there are use cases
-where multi/old outperforms the new multi style.
+approach. For the dense binary system, computational costs still grew as
+:math:`\lambda^{2d}` at large enough :math:`\lambda`. This is equivalent
+to the default style, albeit with a smaller prefactor. For now we are
+keeping the old option in case there are use cases where multi/old
+outperforms the new multi style.
 
 .. note::
 
@@ -118,6 +133,10 @@ Default
 
 ----------
 
-.. _bytype-Shire:
+.. _multi-Shire:
 
-**(Shire)** Shire, Hanley and Stratford, Comp Part Mech, (2020).
+**(Shire)** Shire, Hanley and Stratford, Comp. Part. Mech., (2020).
+
+.. _multi-Monti:
+
+**(Monti)** Monti, Clemmer, Srivastava, Silbert, Grest, and Lechman, Phys. Rev. E, (2022).

src/BPM/bond_bpm.cpp

@@ -40,6 +40,7 @@ BondBPM::BondBPM(LAMMPS *_lmp) :
 {
   overlay_flag = 0;
   prop_atom_flag = 0;
+  break_flag = 1;
   nvalues = 0;
 
   r0_max_estimate = 0.0;
@@ -104,21 +105,21 @@ void BondBPM::init_style()
     }
   } else {
     // Require atoms know about all of their bonds and if they break
-    if (force->newton_bond)
-      error->all(FLERR, "Without overlay/pair, BPM bond styles require Newton bond off");
+    if (force->newton_bond && break_flag)
+      error->all(FLERR, "Without overlay/pair or break/no, BPM bond styles require Newton bond off");
 
     // special lj must be 0 1 1 to censor pair forces between bonded particles
     // special coulomb must be 1 1 1 to ensure all pairs are included in the
     //   neighbor list and 1-3 and 1-4 special bond lists are skipped
     if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0)
       error->all(FLERR,
-                 "Without overlay/pair, BPM bond sytles requires special LJ weights = 0,1,1");
+                 "Without overlay/pair, BPM bond styles requires special LJ weights = 0,1,1");
     if (force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 ||
         force->special_coul[3] != 1.0)
       error->all(FLERR,
-                 "Without overlay/pair, BPM bond sytles requires special Coulomb weights = 1,1,1");
+                 "Without overlay/pair, BPM bond styles requires special Coulomb weights = 1,1,1");
 
-    if (id_fix_dummy) {
+    if (id_fix_dummy && break_flag) {
       id_fix_update = utils::strdup("BPM_UPDATE_SPECIAL_BONDS");
       fix_update_special_bonds = dynamic_cast<FixUpdateSpecialBonds *>(modify->replace_fix(
           id_fix_dummy, fmt::format("{} all UPDATE_SPECIAL_BONDS", id_fix_update), 1));
@@ -188,6 +189,9 @@ void BondBPM::settings(int narg, char **arg)
     } else if (strcmp(arg[iarg], "overlay/pair") == 0) {
       overlay_flag = 1;
       iarg++;
+    } else if (strcmp(arg[iarg], "break/no") == 0) {
+      break_flag = 0;
+      iarg++;
     } else {
       leftover_iarg.push_back(iarg);
       iarg++;
@@ -328,6 +332,8 @@ void BondBPM::read_restart(FILE *fp)
 
 void BondBPM::process_broken(int i, int j)
 {
+  if (break_flag)
+    error->one(FLERR, "BPM bond broke with break/no option");
   if (fix_store_local) {
     for (int n = 0; n < nvalues; n++) (this->*pack_choice[n])(n, i, j);
 

src/BPM/bond_bpm.h

@@ -50,7 +50,7 @@ class BondBPM : public Bond {
   FnPtrPack *pack_choice;    // ptrs to pack functions
   double *output_data;
 
-  int prop_atom_flag, nvalues, overlay_flag;
+  int prop_atom_flag, nvalues, overlay_flag, break_flag;
   int index_x_ref, index_y_ref, index_z_ref;
 
   void pack_id1(int, int, int);

src/RIGID/fix_rigid.cpp

@@ -2523,9 +2523,17 @@ int FixRigid::pack_exchange(int i, double *buf)
   buf[2] = displace[i][0];
   buf[3] = displace[i][1];
   buf[4] = displace[i][2];
-  if (!extended) return 5;
+
+  // must also pack vatom if per-atom virial calculated on this timestep
+  // since vatom is calculated before and after atom migration
 
   int m = 5;
+  if (vflag_atom)
+    for (int k = 0; k < 6; k++)
+      buf[m++] = vatom[i][k];
+
+  if (!extended) return m;
+
   buf[m++] = eflags[i];
   for (int j = 0; j < orientflag; j++)
     buf[m++] = orient[i][j];
@@ -2535,13 +2543,6 @@ int FixRigid::pack_exchange(int i, double *buf)
     buf[m++] = dorient[i][2];
   }
 
-  // must also pack vatom if per-atom virial calculated on this timestep
-  // since vatom is calculated before and after atom migration
-
-  if (vflag_atom)
-    for (int k = 0; k < 6; k++)
-      buf[m++] = vatom[i][k];
-
   return m;
 }
 
@@ -2556,9 +2557,17 @@ int FixRigid::unpack_exchange(int nlocal, double *buf)
   displace[nlocal][0] = buf[2];
   displace[nlocal][1] = buf[3];
   displace[nlocal][2] = buf[4];
-  if (!extended) return 5;
+
+  // must also unpack vatom if per-atom virial calculated on this timestep
+  // since vatom is calculated before and after atom migration
 
   int m = 5;
+  if (vflag_atom)
+    for (int k = 0; k < 6; k++)
+      vatom[nlocal][k] = buf[m++];
+
+  if (!extended) return m;
+
   eflags[nlocal] = static_cast<int> (buf[m++]);
   for (int j = 0; j < orientflag; j++)
     orient[nlocal][j] = buf[m++];
@@ -2568,13 +2577,6 @@ int FixRigid::unpack_exchange(int nlocal, double *buf)
     dorient[nlocal][2] = buf[m++];
   }
 
-  // must also unpack vatom if per-atom virial calculated on this timestep
-  // since vatom is calculated before and after atom migration
-
-  if (vflag_atom)
-    for (int k = 0; k < 6; k++)
-      vatom[nlocal][k] = buf[m++];
-
   return m;
 }
 

src/fix_update_special_bonds.cpp

@@ -72,6 +72,9 @@ void FixUpdateSpecialBonds::setup(int /*vflag*/)
       force->special_coul[3] != 1.0)
     error->all(FLERR, "Fix update/special/bonds requires special Coulomb weights = 1,1,1");
   // Implies neighbor->special_flag = [X, 2, 1, 1]
+
+  if (utils::strmatch(force->pair_style, "^hybrid"))
+    error->all(FLERR, "Cannot use fix update/special/bonds with hybrid pair styles");
 }
 
 /* ----------------------------------------------------------------------

src/neighbor.cpp

@@ -95,6 +95,16 @@ static const char cite_neigh_multi[] =
   "          Detection Applied to Investigate the Quasi-Static Limit},\n"
   " journal = {Computational Particle Mechanics},\n"
   " year = {2020}\n"
+  "@article{Monti2022,\n"
+  " author = {Monti, Joseph M. and Clemmer, Joel T. and Srivastava, \n"
+  "           Ishan and Silbert, Leonardo E. and Grest, Gary S. \n"
+  "           and Lechman, Jeremy B.},\n"
+  " title = {Large-scale frictionless jamming with power-law particle \n"
+  "          size distributions},\n"
+  " journal = {Phys. Rev. E},\n"
+  " volume = {106}\n"
+  " issue = {3}\n"
+  " year = {2022}\n"
   "}\n\n";
 
 // template for factory functions: