diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index 14cabc6c27..b7b5089aed 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -243,7 +243,6 @@ OPT. * :doc:`ttm/grid ` * :doc:`ttm/mod ` * :doc:`tune/kspace ` - * :doc:`update/special/bonds ` * :doc:`vector ` * :doc:`viscosity ` * :doc:`viscous ` diff --git a/doc/src/Howto_bpm.rst b/doc/src/Howto_bpm.rst index 020b57821a..4474e17058 100644 --- a/doc/src/Howto_bpm.rst +++ b/doc/src/Howto_bpm.rst @@ -26,37 +26,37 @@ preserved across run commands and is written to :doc:`binary restart files ` -work differently for BPM bond styles. There are two possible special -bond settings which determine how pair interactions work between bonded -particles. First, one can simply overlay pair interactions such that all -bonded particles also feel pair interactions. This can be accomplished by -simply turning off all special bonds by setting - -.. code-block:: LAMMPS - - special_bonds lj/coul 1 1 1 - -Alternatively, one can censor all pair interactions between bonded particles. +work differently for BPM bond styles. There are two possible settings +which determine how pair interactions work between bonded +particles. +First, one can censor all pair interactions between bonded particles. Unlike :doc:`bond quartic `, this is not done by subtracting pair forces during the bond computation but rather by dynamically updating -the special bond list. To do this, one must both define an instance of -:doc:`fix update/special/bonds ` and have the special bond -settings +the special bond list. This is the default behavior of BPM bond styles +and is done by updating the 1-2 special bond list as bonds break. +To do this, LAMMPS requires :doc:`newton ` bond off such that all +processors containing an atom know when a bond breaks. Additionally, +one must use the following special bond settings .. code-block:: LAMMPS special_bonds lj 0 1 1 coul 1 1 1 -This fix ensures the 1-2 special bond list remains updated as bonds break. The fix -also requires :doc:`newton ` bond off such that whena bond breaks between -atoms across multiple processors, all processors are aware of the event. -The special bond settings then accomplish two tasks. First, they turns off 1-3 and +These settings accomplish two goals. First, they turns off 1-3 and 1-4 special bond lists, which are not currently supported for BPMs. As BPMs often have dense bond networks, generating 1-3 and 1-4 special bond lists is expensive. By setting the lj weight for 1-2 bonds to zero, this censors pairwise interactions. -However, setting a nonzero coul weight for 1-2 bonds ensures all bonded -neighbors are included in the neighbor list. All bonded neighbors must be included -in neighbor lists as they could become unbonded at any timestep. +By setting a nonzero coul weight for 1-2 bonds ensures all bonded neighbors are +still included in the neighbor list in case bonds break between neighbor list builds. + +Alternatively, one can simply overlay pair interactions such that all +bonded particles also feel pair interactions. This can be accomplished by +using the *overlay/pair* keyword in the bond settings and by +using the following special bond settings + +.. code-block:: LAMMPS + + special_bonds lj/coul 1 1 1 Currently there are two types of bonds included in this package. The first bond style, :doc:`bond bpm/spring `, only applies pairwise, diff --git a/doc/src/Packages_details.rst b/doc/src/Packages_details.rst index 281c853eb0..e5e548a341 100644 --- a/doc/src/Packages_details.rst +++ b/doc/src/Packages_details.rst @@ -307,7 +307,6 @@ models for mesoscale simulations of solids and fracture. See the * :doc:`bond_style bpm/spring ` * :doc:`compute nbonds/atom ` * :doc:`fix nve/sphere/bpm ` -* :doc:`fix update/special/bonds ` * :doc:`pair_style bpm/spring ` * examples/bpm diff --git a/doc/src/bond_bpm_rotational.rst b/doc/src/bond_bpm_rotational.rst index 3fcc3d49ec..bffcf0d142 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 = *store/local* +* optional keyword = *overlay/pair* or *store/local* .. parsed-literal:: @@ -121,10 +121,11 @@ or :doc:`read_restart ` commands: * :math:`\gamma_r` (distance*force/seconds/radians units) * :math:`\gamma_t` (distance*force/seconds/radians units) -As bonds can be broken between neighbor list builds, particular -:doc:`special_bonds ` are required. See the `:doc: how to ` -page on BPMs or `:doc: fix update/special/bonds ` -for details. +By default, pair forces are not calculated between bonded particles. +Pair forces can alternatively be overlaid on top of bond forces +using the *overlay/pair* keyword. These settings require specific +:doc:`special_bonds ` settings described in the restrictions. +Further details can be found in the `:doc: how to ` page on BPMs. This bond style tracks broken bonds and can record them using an instance of :doc:`fix store/local ` if the *store/local* keyword is @@ -157,8 +158,18 @@ This bond style can only be used if LAMMPS was built with the BPM package. See the :doc:`Build package ` doc page for more info. -The *bpm/rotational* style requires 1-3 and 1-4 :doc:`special_bonds ` -be turned off using the :doc:`special_bonds ` command. +By default if pair interactions are censored, this bond style requires setting + +.. code-block:: LAMMPS + + special_bonds lj 0 1 1 coul 1 1 1 + +and :doc:`newton ` must be set to bond off. +If the *overlay/pair* option is used, this bond style alternatively requires setting + +.. code-block:: LAMMPS + + special_bonds lj/coul 1 1 1 The *bpm/rotational* style requires :doc:`atom style sphere/bpm `. diff --git a/doc/src/bond_bpm_spring.rst b/doc/src/bond_bpm_spring.rst index e149a45d2e..0235fb52bc 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 = *store/local* +* optional keyword = *overlay/pair* or *store/local* .. parsed-literal:: @@ -88,10 +88,11 @@ or :doc:`read_restart ` commands: * :math:`\epsilon_c` (unit less) * :math:`\gamma` (force/velocity units) -As bonds can be broken between neighbor list builds, particular -:doc:`special_bonds ` are required. See the `:doc: how to ` -page on BPMs or `:doc: fix update/special/bonds ` -for details. +By default, pair forces are not calculated between bonded particles. +Pair forces can alternatively be overlaid on top of bond forces +using the *overlay/pair* keyword. These settings require specific +:doc:`special_bonds ` settings described in the restrictions. +Further details can be found in the `:doc: how to ` page on BPMs. This bond style tracks broken bonds and can record them using an instance of :doc:`fix store/local ` if the *store/local* keyword is @@ -124,8 +125,18 @@ This bond style can only be used if LAMMPS was built with the BPM package. See the :doc:`Build package ` doc page for more info. -The *bpm/spring* style requires 1-3 and 1-4 :doc:`special_bonds ` -be turned off using the :doc:`special_bonds ` command. +By default if pair interactions are censored, this bond style requires setting + +.. code-block:: LAMMPS + + special_bonds lj 0 1 1 coul 1 1 1 + +and :doc:`newton ` must be set to bond off. +If the *overlay/pair* option is used, this bond style alternatively requires setting + +.. code-block:: LAMMPS + + special_bonds lj/coul 1 1 1 Related commands """""""""""""""" diff --git a/doc/src/fix_update_special_bonds.rst b/doc/src/fix_update_special_bonds.rst deleted file mode 100644 index be693118cd..0000000000 --- a/doc/src/fix_update_special_bonds.rst +++ /dev/null @@ -1,55 +0,0 @@ -.. index:: fix update/special/bonds - -fix update/special/bonds command -====================== - -Syntax -"""""" - -.. parsed-literal:: - - fix ID group-ID update/special/bonds - -* ID, group-ID are documented in :doc:`fix ` command -* update/special/bonds = style name of this fix command - -Examples -"""""""" - -.. code-block:: LAMMPS - - fix 1 all update/special/bonds - -Description -""""""""""" - -This fix is used to update the 1-2 special bond list for BPM bond styles. -This feature is used to censor pair forces between bonded particles. -See the :doc:`BPM how to ` for more information. - ----------- - -Restart, fix_modify, output, run start/stop, minimize info -""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" - -No information about this fix is written to :doc:`binary restart files `. None of the :doc:`fix_modify ` options -are relevant to this fix. No global or per-atom quantities are stored -by this fix for access by various :doc:`output commands `. -No parameter of this fix can be used with the *start/stop* keywords of -the :doc:`run ` command. This fix is not invoked during :doc:`energy minimization `. - -Restrictions -"""""""""""" - -This fix requires :doc:`newton ` bond off. - -Related commands -"""""""""""""""" - -:doc:`bond bpm/rotational ` - -Default -""""""" - -none - diff --git a/examples/bpm/impact/impact_rotational_bpm.lmp b/examples/bpm/impact/impact_rotational_bpm.lmp index f79a72fd9f..3c4667ca1b 100644 --- a/examples/bpm/impact/impact_rotational_bpm.lmp +++ b/examples/bpm/impact/impact_rotational_bpm.lmp @@ -28,7 +28,6 @@ bond_coeff 2 1.0 0.2 0.02 0.02 0.20 0.04 0.04 0.04 0.1 0.02 0.002 0.002 fix 1 all nve/sphere/bpm fix 2 all store/local 100 3 -fix 3 all update/special/bonds create_bonds many plate plate 1 0.0 1.5 create_bonds many projectile projectile 2 0.0 1.5 diff --git a/examples/bpm/impact/impact_spring_bpm.lmp b/examples/bpm/impact/impact_spring_bpm.lmp index bbb968bee2..ef5506218c 100644 --- a/examples/bpm/impact/impact_spring_bpm.lmp +++ b/examples/bpm/impact/impact_spring_bpm.lmp @@ -30,7 +30,6 @@ bond_coeff 2 1.0 0.20 1.0 fix 1 all nve fix 2 all store/local 100 3 -fix 3 all update/special/bonds create_bonds many plate plate 1 0.0 1.5 create_bonds many projectile projectile 2 0.0 1.5 diff --git a/examples/bpm/pour/pour_bpm.lmp b/examples/bpm/pour/pour_bpm.lmp index bb582413c3..98dca9d3a7 100644 --- a/examples/bpm/pour/pour_bpm.lmp +++ b/examples/bpm/pour/pour_bpm.lmp @@ -26,7 +26,6 @@ fix 2 all wall/gran/region hertz/history 1.0 NULL 0.5 NULL 0.1 1 reg fix 3 all gravity 1e-4 vector 0 0 -1 fix 4 all deposit 40 0 1500 712511343 mol my_mol region dropzone near 2.0 vz -0.05 -0.05 fix 5 all nve/sphere/bpm -fix 6 all update/special/bonds timestep 0.05 thermo_style custom step ke pe pxx pyy pzz c_tbond diff --git a/src/BPM/bond_bpm.cpp b/src/BPM/bond_bpm.cpp index eacfd064e5..df5bad4cda 100644 --- a/src/BPM/bond_bpm.cpp +++ b/src/BPM/bond_bpm.cpp @@ -33,10 +33,11 @@ BondBPM::BondBPM(LAMMPS *lmp) : Bond(lmp) { id_fix_store_local = nullptr; id_fix_prop_atom = nullptr; + id_fix_update = nullptr; fix_store_local = nullptr; - fix_update_special_bonds = nullptr; fix_bond_history = nullptr; + overlay_flag = 0; prop_atom_flag = 0; nvalues = 0; output_data = nullptr; @@ -44,6 +45,12 @@ BondBPM::BondBPM(LAMMPS *lmp) : Bond(lmp) r0_max_estimate = 0.0; max_stretch = 1.0; + + // create dummy fix as placeholder for FixUpdateSpecialBonds + // this is so final order of Modify:fix will conform to input script + + id_fix_dummy = utils::strdup("BPM_DUMMY"); + modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy)); } /* ---------------------------------------------------------------------- */ @@ -53,6 +60,13 @@ BondBPM::~BondBPM() delete [] pack_choice; delete [] id_fix_store_local; delete [] id_fix_prop_atom; + + if (id_fix_dummy) modify->delete_fix(id_fix_dummy); + if (id_fix_update) modify->delete_fix(id_fix_update); + + delete [] id_fix_dummy; + delete [] id_fix_update; + memory->destroy(output_data); } @@ -70,14 +84,38 @@ void BondBPM::init_style() fix_store_local->nvalues = nvalues; } - ifix = modify->find_fix_by_style("update/special/bonds"); - if (ifix >= 0) - fix_update_special_bonds = (FixUpdateSpecialBonds *) modify->fix[ifix]; - else { - fix_update_special_bonds = nullptr; + + if (overlay_flag) { if (force->special_lj[1] != 1.0) - error->all(FLERR, "Without fix update/special/bonds, BPM bond styles " + error->all(FLERR, "With overlay/pair, BPM bond styles " "require special_bonds weight of 1.0 for first neighbors"); + if (id_fix_update) { + modify->delete_fix(id_fix_update); + delete [] id_fix_update; + id_fix_update = nullptr; + } + } 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 sytles 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"); + 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"); + + if (id_fix_dummy) { + id_fix_update = utils::strdup("BPM_update_special_bonds"); + fix_update_special_bonds = (FixUpdateSpecialBonds *) modify->replace_fix(id_fix_dummy, + fmt::format("{} all update/special/bonds", id_fix_update),1); + delete [] id_fix_dummy; + id_fix_dummy = nullptr; + } } if (force->angle || force->dihedral || force->improper) @@ -105,38 +143,44 @@ void BondBPM::settings(int narg, char **arg) { int iarg = 0; while (iarg < narg) { - if (id_fix_store_local) { - if (strcmp(arg[iarg], "id1") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_id1; - } else if (strcmp(arg[iarg], "id2") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_id2; - } else if (strcmp(arg[iarg], "time") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_time; - } else if (strcmp(arg[iarg], "x") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_x; - } else if (strcmp(arg[iarg], "y") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_y; - } else if (strcmp(arg[iarg], "z") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_z; - } else if (strcmp(arg[iarg], "x/ref") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_x_ref; - prop_atom_flag = 1; - } else if (strcmp(arg[iarg], "y/ref") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_y_ref; - prop_atom_flag = 1; - } else if (strcmp(arg[iarg], "z/ref") == 0) { - pack_choice[nvalues++] = &BondBPM::pack_z_ref; - prop_atom_flag = 1; - } else { - error->all(FLERR, "Illegal bond_style command"); - } - } else if (strcmp(arg[iarg], "store/local") == 0) { + if (strcmp(arg[iarg], "store/local") == 0) { id_fix_store_local = utils::strdup(arg[iarg+1]); - iarg ++; nvalues = 0; pack_choice = new FnPtrPack[narg - iarg - 1]; + iarg += 2; + while (iarg < narg) { + if (strcmp(arg[iarg], "id1") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_id1; + } else if (strcmp(arg[iarg], "id2") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_id2; + } else if (strcmp(arg[iarg], "time") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_time; + } else if (strcmp(arg[iarg], "x") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_x; + } else if (strcmp(arg[iarg], "y") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_y; + } else if (strcmp(arg[iarg], "z") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_z; + } else if (strcmp(arg[iarg], "x/ref") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_x_ref; + prop_atom_flag = 1; + } else if (strcmp(arg[iarg], "y/ref") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_y_ref; + prop_atom_flag = 1; + } else if (strcmp(arg[iarg], "z/ref") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_z_ref; + prop_atom_flag = 1; + } else { + break; + } + iarg ++; + } + } else if (strcmp(arg[iarg], "overlay/pair") == 0) { + overlay_flag = 1; + iarg ++; + } else { + error->all(FLERR, "Illegal pair_style command"); } - iarg ++; } if (id_fix_store_local) { @@ -147,7 +191,7 @@ void BondBPM::settings(int narg, char **arg) // Use store property to save reference positions as it can transfer to ghost atoms if (prop_atom_flag == 1) { - id_fix_prop_atom = utils::strdup("BPM_PROPERTY_ATOM"); + id_fix_prop_atom = utils::strdup("BPM_property_atom"); int ifix = modify->find_fix(id_fix_prop_atom); if (ifix < 0) { modify->add_fix(fmt::format("{} all property/atom " diff --git a/src/BPM/bond_bpm.h b/src/BPM/bond_bpm.h index 53a202802d..4836a08a64 100644 --- a/src/BPM/bond_bpm.h +++ b/src/BPM/bond_bpm.h @@ -25,7 +25,7 @@ class BondBPM : public Bond { virtual void compute(int, int) = 0; virtual void coeff(int, char **) = 0; virtual void init_style(); - void settings(int, char **); + void settings(int, char **); double equilibrium_distance(int); void write_restart(FILE *){}; void read_restart(FILE *){}; @@ -35,18 +35,19 @@ class BondBPM : public Bond { protected: double r0_max_estimate; double max_stretch; - + + char *id_fix_dummy, *id_fix_update; char *id_fix_store_local, *id_fix_prop_atom; class FixStoreLocal *fix_store_local; - class FixUpdateSpecialBonds *fix_update_special_bonds; class FixBondHistory *fix_bond_history; + class FixUpdateSpecialBonds *fix_update_special_bonds; - void process_broken(int, int); + void process_broken(int, int); typedef void (BondBPM::*FnPtrPack)(int,int,int); FnPtrPack *pack_choice; // ptrs to pack functions double *output_data; - int prop_atom_flag, nvalues; + int prop_atom_flag, nvalues, overlay_flag; int index_x_ref, index_y_ref, index_z_ref; void pack_id1(int,int,int); @@ -57,7 +58,7 @@ class BondBPM : public Bond { void pack_z(int,int,int); void pack_x_ref(int,int,int); void pack_y_ref(int,int,int); - void pack_z_ref(int,int,int); + void pack_z_ref(int,int,int); }; } // namespace LAMMPS_NS diff --git a/src/MISC/pair_tracker.cpp b/src/MISC/pair_tracker.cpp index 0297cbc938..cd6ed67615 100644 --- a/src/MISC/pair_tracker.cpp +++ b/src/MISC/pair_tracker.cpp @@ -385,33 +385,35 @@ void PairTracker::init_style() "Pair tracker incompatible with granular pairstyles that extend beyond contact"); // check for FixFreeze and set freeze_group_bit - int ifix = modify->find_fix_by_style("^freeze"); - if (ifix < 0) freeze_group_bit = 0; - else freeze_group_bit = modify->fix[ifix]->groupbit; + + auto fixlist = modify->get_fix_by_style("^freeze"); + if (fixlist.size() == 0) + freeze_group_bit = 0; + else if (fixlist.size() > 1) + error->all(FLERR, "Only one fix freeze command at a time allowed"); + else + freeze_group_bit = fixlist.front()->groupbit; // check for FixPour and FixDeposit so can extract particle radii - int ipour; - for (ipour = 0; ipour < modify->nfix; ipour++) - if (strcmp(modify->fix[ipour]->style, "pour") == 0) break; - if (ipour == modify->nfix) ipour = -1; - int idep; - for (idep = 0; idep < modify->nfix; idep++) - if (strcmp(modify->fix[idep]->style, "deposit") == 0) break; - if (idep == modify->nfix) idep = -1; + auto pours = modify->get_fix_by_style("^pour"); + auto deps = modify->get_fix_by_style("^deposit"); // set maxrad_dynamic and maxrad_frozen for each type // include future FixPour and FixDeposit particles as dynamic + int itype; for (i = 1; i <= atom->ntypes; i++) { onerad_dynamic[i] = onerad_frozen[i] = 0.0; - if (ipour >= 0) { + for (auto &ipour : pours) { itype = i; - onerad_dynamic[i] = *((double *) modify->fix[ipour]->extract("radius", itype)); + double maxrad = *((double *) ipour->extract("radius", itype)); + if (maxrad > 0.0) onerad_dynamic[i] = maxrad; } - if (idep >= 0) { + for (auto &idep : deps) { itype = i; - onerad_dynamic[i] = *((double *) modify->fix[idep]->extract("radius", itype)); + double maxrad = *((double *) idep->extract("radius", itype)); + if (maxrad > 0.0) onerad_dynamic[i] = maxrad; } }