From 120f5cf7f120f4aa30c677ddaef1cee77129dc9d Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 16 Mar 2023 15:20:18 -0600 Subject: [PATCH 1/7] Minor patches to BPM, multi, and rigid --- doc/src/bond_bpm_rotational.rst | 9 ++++++++- doc/src/bond_bpm_spring.rst | 28 ++++++++++++++++++-------- doc/src/comm_modify.rst | 1 + doc/src/neighbor.rst | 21 +++++++++++++++----- src/BPM/bond_bpm.cpp | 16 ++++++++++----- src/BPM/bond_bpm.h | 2 +- src/RIGID/fix_rigid.cpp | 34 +++++++++++++++++--------------- src/fix_update_special_bonds.cpp | 3 +++ src/neighbor.cpp | 10 ++++++++++ 9 files changed, 88 insertions(+), 36 deletions(-) diff --git a/doc/src/bond_bpm_rotational.rst b/doc/src/bond_bpm_rotational.rst index bc2f383828..0975fb35f6 100644 --- a/doc/src/bond_bpm_rotational.rst +++ b/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 """""""" @@ -138,6 +141,10 @@ the *overlay/pair* keyword. These settings require specific restrictions. Further details can be found in the `:doc: how to ` page on BPMs. +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 diff --git a/doc/src/bond_bpm_spring.rst b/doc/src/bond_bpm_spring.rst index de31357afe..fe5c63216d 100644 --- a/doc/src/bond_bpm_spring.rst +++ b/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 """""""" @@ -45,16 +48,17 @@ Examples Description """"""""""" -The *bpm/spring* bond style computes forces and torques 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 -written to :doc:`binary restart files ` such that restarting -the system will not reset the reference state of a bond. +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 written to +:doc:`binary restart files ` such that restarting 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) `. The force has a magnitude of .. math:: @@ -103,6 +107,10 @@ the *overlay/pair* keyword. These settings require specific restrictions. Further details can be found in the `:doc: how to ` page on BPMs. +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 @@ -198,6 +206,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). diff --git a/doc/src/comm_modify.rst b/doc/src/comm_modify.rst index d0526b792b..14f55378c2 100644 --- a/doc/src/comm_modify.rst +++ b/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 ` 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 diff --git a/doc/src/neighbor.rst b/doc/src/neighbor.rst index 663170ef47..bf9a13923c 100644 --- a/doc/src/neighbor.rst +++ b/doc/src/neighbor.rst @@ -59,9 +59,16 @@ 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) `. -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) ` +and Monti et al. :ref:`(Monti) `. This imposes some extra +setup overhead, but the searches themselves may be much faster. For +instance in a dense binary system with a ratio of the size of the largest +to smallest collection bin :math:`\lamda`, the computational costs of +building a default neighbor list grows as :math:`\lamda^6` while the costs +for *multi* grows as :math:`\lamda^3`, equivalent to the cost of force +evaluations, as identified in Monti et al. :ref:`(Monti) `. + +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 @@ -118,6 +125,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). diff --git a/src/BPM/bond_bpm.cpp b/src/BPM/bond_bpm.cpp index 5a20059860..724ec0a7bc 100644 --- a/src/BPM/bond_bpm.cpp +++ b/src/BPM/bond_bpm.cpp @@ -39,6 +39,7 @@ BondBPM::BondBPM(LAMMPS *_lmp) : { overlay_flag = 0; prop_atom_flag = 0; + break_flag = 1; nvalues = 0; r0_max_estimate = 0.0; @@ -103,21 +104,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(modify->replace_fix( id_fix_dummy, fmt::format("{} all UPDATE_SPECIAL_BONDS", id_fix_update), 1)); @@ -187,6 +188,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); diff --git a/src/BPM/bond_bpm.h b/src/BPM/bond_bpm.h index 7a32ec2d0d..470307add1 100644 --- a/src/BPM/bond_bpm.h +++ b/src/BPM/bond_bpm.h @@ -52,7 +52,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); diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index f0bcf630b7..caae812468 100644 --- a/src/RIGID/fix_rigid.cpp +++ b/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 (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; } diff --git a/src/fix_update_special_bonds.cpp b/src/fix_update_special_bonds.cpp index 581918dbd3..9ffd84c7de 100644 --- a/src/fix_update_special_bonds.cpp +++ b/src/fix_update_special_bonds.cpp @@ -73,6 +73,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"); } /* ---------------------------------------------------------------------- diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 298ec29994..6696971e7d 100644 --- a/src/neighbor.cpp +++ b/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: From ac9389f5cbdaa2d0e16bb805c432498fb391fcf7 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 20 Mar 2023 14:24:21 -0600 Subject: [PATCH 2/7] Slight rewrite --- doc/src/neighbor.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/neighbor.rst b/doc/src/neighbor.rst index bf9a13923c..b5ca4d15a8 100644 --- a/doc/src/neighbor.rst +++ b/doc/src/neighbor.rst @@ -66,7 +66,7 @@ instance in a dense binary system with a ratio of the size of the largest to smallest collection bin :math:`\lamda`, the computational costs of building a default neighbor list grows as :math:`\lamda^6` while the costs for *multi* grows as :math:`\lamda^3`, equivalent to the cost of force -evaluations, as identified in Monti et al. :ref:`(Monti) `. +evaluations, as argued in Monti et al. :ref:`(Monti) `. By default in *multi*, each atom type defines a separate collection of particles. For systems where two or more atom types have the same From d15e13d4756270b42fb7762282028cd69916fddf Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 20 Mar 2023 14:56:39 -0600 Subject: [PATCH 3/7] Reverting mistakenly deleted line, fixing duplicated text in granular doc --- doc/src/Modify_gran_sub_mod.rst | 2 +- doc/src/bond_bpm_spring.rst | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/doc/src/Modify_gran_sub_mod.rst b/doc/src/Modify_gran_sub_mod.rst index e1d559e6bf..b3b4a13cac 100644 --- a/doc/src/Modify_gran_sub_mod.rst +++ b/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 diff --git a/doc/src/bond_bpm_spring.rst b/doc/src/bond_bpm_spring.rst index 1de814b488..93908710d1 100644 --- a/doc/src/bond_bpm_spring.rst +++ b/doc/src/bond_bpm_spring.rst @@ -48,6 +48,8 @@ Examples Description """"""""""" +.. versionadded:: 4May2022 + 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 From badfd0bc40072d2da1a4de14eb0c098682b360bb Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 20 Mar 2023 15:34:52 -0600 Subject: [PATCH 4/7] Specifying dimensions, lamda->lambda --- doc/src/neighbor.rst | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/src/neighbor.rst b/doc/src/neighbor.rst index b5ca4d15a8..c4a2caac81 100644 --- a/doc/src/neighbor.rst +++ b/doc/src/neighbor.rst @@ -62,11 +62,11 @@ sets of bins are then used to construct the neighbor lists as as further described by Shire, Hanley, and Stratford :ref:`(Shire) ` and Monti et al. :ref:`(Monti) `. This imposes some extra setup overhead, but the searches themselves may be much faster. For -instance in a dense binary system with a ratio of the size of the largest -to smallest collection bin :math:`\lamda`, the computational costs of -building a default neighbor list grows as :math:`\lamda^6` while the costs -for *multi* grows as :math:`\lamda^3`, equivalent to the cost of force -evaluations, as argued in Monti et al. :ref:`(Monti) `. +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^d` 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) `. By default in *multi*, each atom type defines a separate collection of particles. For systems where two or more atom types have the same From 1101383b51b966fc59b354e6c6f28cef84ee1db5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 20 Mar 2023 22:07:02 -0400 Subject: [PATCH 5/7] add versionadded tags --- doc/src/bond_bpm_rotational.rst | 2 ++ doc/src/bond_bpm_spring.rst | 2 ++ 2 files changed, 4 insertions(+) diff --git a/doc/src/bond_bpm_rotational.rst b/doc/src/bond_bpm_rotational.rst index 19fd607169..0b3fa4442c 100644 --- a/doc/src/bond_bpm_rotational.rst +++ b/doc/src/bond_bpm_rotational.rst @@ -143,6 +143,8 @@ the *overlay/pair* keyword. These settings require specific restrictions. Further details can be found in the `:doc: how to ` 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. diff --git a/doc/src/bond_bpm_spring.rst b/doc/src/bond_bpm_spring.rst index 93908710d1..5391d81420 100644 --- a/doc/src/bond_bpm_spring.rst +++ b/doc/src/bond_bpm_spring.rst @@ -109,6 +109,8 @@ the *overlay/pair* keyword. These settings require specific restrictions. Further details can be found in the `:doc: how to ` 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. From 669397b0929024a5f82088502f0c6ce4cdddaae4 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Wed, 22 Mar 2023 09:36:29 -0600 Subject: [PATCH 6/7] fixing exponent in multi documentation --- doc/src/neighbor.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/neighbor.rst b/doc/src/neighbor.rst index c4a2caac81..4ee3d859f8 100644 --- a/doc/src/neighbor.rst +++ b/doc/src/neighbor.rst @@ -64,7 +64,7 @@ and Monti et al. :ref:`(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^d` while +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) `. From 1370f0571486880179310483b1c6977dd517403d Mon Sep 17 00:00:00 2001 From: jtclemm Date: Wed, 22 Mar 2023 11:33:17 -0600 Subject: [PATCH 7/7] Elaborating on the scaling of multi --- doc/src/neighbor.rst | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/doc/src/neighbor.rst b/doc/src/neighbor.rst index 4ee3d859f8..9dcdd1daf2 100644 --- a/doc/src/neighbor.rst +++ b/doc/src/neighbor.rst @@ -61,12 +61,17 @@ 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) ` and Monti et al. :ref:`(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 +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) `. +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 @@ -82,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::