diff --git a/doc/src/Howto_bpm.rst b/doc/src/Howto_bpm.rst index 86e091fcef..3dee00fc15 100644 --- a/doc/src/Howto_bpm.rst +++ b/doc/src/Howto_bpm.rst @@ -42,7 +42,8 @@ Currently, there are two types of bonds included in the BPM package. The first bond style, :doc:`bond bpm/spring `, only applies pairwise, central body forces. Point particles must have :doc:`bond atom style ` and may be thought of as nodes in a spring -network. Alternatively, the second bond style, :doc:`bond bpm/rotational +network. An optional multibody term can be used to adjust the network's +Poisson's ratio. Alternatively, the second bond style, :doc:`bond bpm/rotational `, resolves tangential forces and torques arising with the shearing, bending, and twisting of the bond due to rotation or displacement of particles. Particles are similar to those used in the @@ -55,8 +56,9 @@ orientation similar to :doc:`fix nve/asphere `. In addition to bond styles, a new pair style :doc:`pair bpm/spring ` was added to accompany the bpm/spring bond -style. This pair style is simply a hookean repulsion with similar -velocity damping as its sister bond style. +style. By default, this pair style is simply a hookean repulsion with +similar velocity damping as its sister bond style, but optional +arguments can be used to modify the force. ---------- diff --git a/doc/src/bond_bpm_spring.rst b/doc/src/bond_bpm_spring.rst index 81582bd5ea..4bc5cd3463 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* or *break* +* optional keyword = *overlay/pair* or *store/local* or *smooth* or *break* or *volume/factor* .. parsed-literal:: @@ -36,6 +36,9 @@ Syntax *break* value = *yes* or *no* indicates whether bonds break during a run + *volume/factor* value = *yes* or *no* + indicates whether forces include the volumetric contribution + Examples """""""" @@ -44,6 +47,9 @@ Examples bond_style bpm/spring bond_coeff 1 1.0 0.05 0.1 + bond_style bpm/spring volume/factor yes + bond_coeff 1 1.0 0.05 0.1 0.5 + bond_style bpm/spring myfix 1000 time id1 id2 dump 1 all local 1000 dump.broken f_myfix[1] f_myfix[2] f_myfix[3] dump_modify 1 write_header no @@ -97,15 +103,6 @@ approach the critical strain w = 1.0 - \left( \frac{r - r_0}{r_0 \epsilon_c} \right)^8 . -The following coefficients must be defined for each bond type via the -:doc:`bond_coeff ` command as in the example above, or in -the data file or restart files read by the :doc:`read_data -` or :doc:`read_restart ` commands: - -* :math:`k` (force/distance units) -* :math:`\epsilon_c` (unit less) -* :math:`\gamma` (force/velocity units) - If the *normalize* keyword is set to *yes*, the elastic bond force will be normalized by :math:`r_0` such that :math:`k` must be given in force units. @@ -123,6 +120,43 @@ during a simulation run. This will prevent some unnecessary calculation. However, if a bond reaches a strain greater than :math:`\epsilon_c`, it will trigger an error. +.. versionadded:: TBD + +The *volume/factor* keyword toggles whether an additional multibody +contribution is added to he force using the formulation in +:ref:`(Clemmer2) `, + +.. math:: + + \alpha_v \left(\left[\frac{V_i + V_j}{V_{0,i} + V_{0,j}}\right]^{1/3} - \frac{r_{ij}}{r_{0,ij}}\right) + +where :math:`\alpha_v` is a user specified coefficient and :math:`V_i` +and :math:`V_{0,i}` are estimates of the current and local volume +of atom :math:`i`. These volumes are calculated as the sum of current +or initial bond lengths cubed. In 2D, the volume is replaced with an area +calculated using bond lengths squared and the cube root in the above equation +is accordingly replaced with a square root. This approximation assumes bonds +are evenly distributed on a spherical surface and neglects constant prefactors +which are irrelevant since only the ratio of volumes matters. This term may be +used to adjust the Poisson's ratio. + +If a bond is broken (or created), :math:`V_{0,i}` is updated by subtracting +(or adding) that bond's contribution. + +The following coefficients must be defined for each bond type via the +:doc:`bond_coeff ` command as in the example above, or in +the data file or restart files read by the :doc:`read_data +` or :doc:`read_restart ` commands: + +* :math:`k` (force/distance units) +* :math:`\epsilon_c` (unit less) +* :math:`\gamma` (force/velocity units) + +Additionally, if *volume/factor* is set to *yes*, a fourth coefficient +must be provided: + +* :math:`a_v` (force units) + 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 @@ -213,7 +247,7 @@ Related commands Default """"""" -The option defaults are *overlay/pair* = *no*, *smooth* = *yes*, *normalize* = *no*, and *break* = *yes* +The option defaults are *overlay/pair* = *no*, *smooth* = *yes*, *normalize* = *no*, *break* = *yes*, and *volume/factor* = *no* ---------- @@ -224,3 +258,7 @@ The option defaults are *overlay/pair* = *no*, *smooth* = *yes*, *normalize* = * .. _Groot4: **(Groot)** Groot and Warren, J Chem Phys, 107, 4423-35 (1997). + +.. _multibody-Clemmer: + +**(Clemmer2)** Clemmer, Monti, Lechman, Soft Matter, 20, 1702 (2024). diff --git a/doc/src/pair_bpm_spring.rst b/doc/src/pair_bpm_spring.rst index 7dfa9bc12c..068efff577 100644 --- a/doc/src/pair_bpm_spring.rst +++ b/doc/src/pair_bpm_spring.rst @@ -8,7 +8,14 @@ Syntax .. code-block:: LAMMPS - pair_style bpm/spring + pair_style bpm/spring keyword value ... + +* optional keyword = *anharmonic* + + .. parsed-literal:: + + *anharmonic* value = *yes* or *no* + whether forces include the anharmonic term Examples """""""" @@ -17,7 +24,8 @@ Examples pair_style bpm/spring pair_coeff * * 1.0 1.0 1.0 - pair_coeff 1 1 1.0 1.0 1.0 + pair_style bpm/spring anharmonic yes + pair_coeff 1 1 1.0 1.0 1.0 50.0 Description """"""""""" @@ -28,12 +36,16 @@ Style *bpm/spring* computes pairwise forces with the formula .. math:: - F = k (r - r_c) + F = k (r - r_c) + k_a (r - r_c)^3 -where :math:`k` is a stiffness and :math:`r_c` is the cutoff length. -An additional damping force is also applied to interacting -particles. The force is proportional to the difference in the -normal velocity of particles +where :math:`k` is a stiffness, :math:`r_c` is the cutoff +length, and :math:`k_a` is an optional anharmonic cubic prefactor +that can be enabled using the *anharmonic* keyword. The anharmonic +term may be useful in scenarios that need to prevent large particle overlap. + +An additional damping force is also applied to interacting particles. +The force is proportional to the difference in the normal velocity of +particles .. math:: @@ -73,6 +85,12 @@ commands, or by mixing as described below: * :math:`r_c` (distance units) * :math:`\gamma` (force/velocity units) +.. versionadded:: TBD + +Additionally, if *anharmonic* is set to *yes*, a fourth coefficient +must be provided: + +* :math:`k_a` (force/distance\^3 units) ---------- @@ -117,4 +135,5 @@ Related commands Default """"""" -none +The option defaults are *anharmonic* = *no* + diff --git a/examples/bpm/impact/in.bpm.impact.rotational b/examples/bpm/impact/in.bpm.impact.rotational index 698fa22fb8..caca37ddc0 100644 --- a/examples/bpm/impact/in.bpm.impact.rotational +++ b/examples/bpm/impact/in.bpm.impact.rotational @@ -45,7 +45,7 @@ thermo 100 thermo_modify lost ignore lost/bond ignore #dump 1 all custom 100 atomDump id radius x y z c_nbond -dump 2 all local 100 brokenDump f_brkbond[1] f_brkbond[2] f_brkbond[3] -dump_modify 2 header no +#dump 2 all local 100 brokenDump f_brkbond[1] f_brkbond[2] f_brkbond[3] +#dump_modify 2 header no run 7500 diff --git a/examples/bpm/impact/in.bpm.impact.spring b/examples/bpm/impact/in.bpm.impact.spring index 7c5c56841b..d135de6f9e 100644 --- a/examples/bpm/impact/in.bpm.impact.spring +++ b/examples/bpm/impact/in.bpm.impact.spring @@ -47,7 +47,7 @@ thermo 100 thermo_modify lost ignore lost/bond ignore #dump 1 all custom 100 atomDump id x y z c_nbond -dump 2 all local 100 brokenDump f_brkbond[1] f_brkbond[2] f_brkbond[3] -dump_modify 2 header no +#dump 2 all local 100 brokenDump f_brkbond[1] f_brkbond[2] f_brkbond[3] +#dump_modify 2 header no run 7500 diff --git a/src/BPM/bond_bpm.cpp b/src/BPM/bond_bpm.cpp index 130080b349..3ed4224a94 100644 --- a/src/BPM/bond_bpm.cpp +++ b/src/BPM/bond_bpm.cpp @@ -33,31 +33,32 @@ using namespace LAMMPS_NS; static const char cite_bpm[] = - "BPM bond style: doi:10.1039/D3SM01373A\n\n" - "@Article{Clemmer2024,\n" - " author = {Clemmer, Joel T. and Monti, Joseph M. and Lechman, Jeremy B.},\n" - " title = {A soft departure from jamming: the compaction of deformable\n" - " granular matter under high pressures},\n" - " journal = {Soft Matter},\n" - " year = 2024,\n" - " volume = 20,\n" - " number = 8,\n" - " pages = {1702--1718}\n" - "}\n\n"; + "BPM bond style: doi:10.1039/D3SM01373A\n\n" + "@Article{Clemmer2024,\n" + " author = {Clemmer, Joel T. and Monti, Joseph M. and Lechman, Jeremy B.},\n" + " title = {A soft departure from jamming: the compaction of deformable\n" + " granular matter under high pressures},\n" + " journal = {Soft Matter},\n" + " year = 2024,\n" + " volume = 20,\n" + " number = 8,\n" + " pages = {1702--1718}\n" + "}\n\n"; /* ---------------------------------------------------------------------- */ BondBPM::BondBPM(LAMMPS *_lmp) : - Bond(_lmp), id_fix_dummy(nullptr), id_fix_dummy2(nullptr), id_fix_update(nullptr), - id_fix_bond_history(nullptr), id_fix_store_local(nullptr), id_fix_prop_atom(nullptr), - fix_store_local(nullptr), fix_bond_history(nullptr), fix_update_special_bonds(nullptr), - pack_choice(nullptr), output_data(nullptr) + Bond(_lmp), id_fix_dummy_special(nullptr), id_fix_dummy_history(nullptr), + id_fix_update_special_bonds(nullptr), id_fix_bond_history(nullptr), id_fix_store_local(nullptr), + id_fix_property_atom(nullptr), fix_store_local(nullptr), fix_bond_history(nullptr), + fix_update_special_bonds(nullptr), pack_choice(nullptr), output_data(nullptr) { overlay_flag = 0; ignore_special_flag = 0; - prop_atom_flag = 0; + property_atom_flag = 0; break_flag = 1; nvalues = 0; + writedata = 0; nhistory = 0; update_flag = 0; @@ -69,11 +70,11 @@ BondBPM::BondBPM(LAMMPS *_lmp) : // this is so final order of Modify:fix will conform to input script // BondHistory technically only needs this if updateflag = 1 - id_fix_dummy = utils::strdup(fmt::format("BPM_DUMMY_{}", instance_total)); - modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy)); + id_fix_dummy_special = utils::strdup(fmt::format("BPM_DUMMY_SPECIAL_{}", instance_total)); + modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy_special)); - id_fix_dummy2 = utils::strdup(fmt::format("BPM_DUMMY2_{}", instance_total)); - modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy2)); + id_fix_dummy_history = utils::strdup(fmt::format("BPM_DUMMY_HISTORY_{}", instance_total)); + modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy_history)); if (lmp->citeme) lmp->citeme->add(cite_bpm); } @@ -84,19 +85,17 @@ BondBPM::~BondBPM() { delete[] pack_choice; - if (id_fix_dummy) modify->delete_fix(id_fix_dummy); - if (id_fix_dummy2) modify->delete_fix(id_fix_dummy2); - if (id_fix_update) modify->delete_fix(id_fix_update); + if (id_fix_dummy_special) modify->delete_fix(id_fix_dummy_special); + if (id_fix_dummy_history) modify->delete_fix(id_fix_dummy_history); + if (id_fix_update_special_bonds) modify->delete_fix(id_fix_update_special_bonds); if (fix_bond_history) modify->delete_fix(id_fix_bond_history); if (id_fix_store_local) modify->delete_fix(id_fix_store_local); - if (id_fix_prop_atom) modify->delete_fix(id_fix_prop_atom); + if (id_fix_property_atom) modify->delete_fix(id_fix_property_atom); - delete[] id_fix_dummy; - delete[] id_fix_dummy2; - delete[] id_fix_update; + delete[] id_fix_update_special_bonds; delete[] id_fix_bond_history; delete[] id_fix_store_local; - delete[] id_fix_prop_atom; + delete[] id_fix_property_atom; memory->destroy(output_data); } @@ -116,19 +115,22 @@ void BondBPM::init_style() if (!ignore_special_flag) { if (overlay_flag) { - if (force->special_lj[1] != 1.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 || - force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 || force->special_coul[3] != 1.0) + if (force->special_lj[1] != 1.0 || force->special_lj[2] != 1.0 || + force->special_lj[3] != 1.0 || force->special_coul[1] != 1.0 || + force->special_coul[2] != 1.0 || force->special_coul[3] != 1.0) error->all(FLERR, - "With overlay/pair yes, BPM bond styles require a value of 1.0 for all special_bonds weights"); - if (id_fix_update) { - modify->delete_fix(id_fix_update); - delete[] id_fix_update; - id_fix_update = nullptr; + "With overlay/pair yes, BPM bond styles require a value of 1.0 for all " + "special_bonds weights"); + if (id_fix_update_special_bonds) { + modify->delete_fix(id_fix_update_special_bonds); + delete[] id_fix_update_special_bonds; + id_fix_update_special_bonds = nullptr; } } else { // Require atoms know about all of their bonds and if they break if (force->newton_bond && break_flag) - error->all(FLERR, "With overlay/pair no, or break yes, BPM bond styles require Newton bond off"); + error->all(FLERR, + "With overlay/pair no, or break yes, BPM bond styles require Newton bond off"); // special lj must be 0 1 1 to censor pair forces between bonded particles if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0) @@ -136,23 +138,27 @@ void BondBPM::init_style() "With overlay/pair no, BPM bond styles require special LJ weights = 0,1,1"); // if bonds can break, 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 (break_flag && (force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 || - force->special_coul[3] != 1.0)) + if (break_flag && + (force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 || + force->special_coul[3] != 1.0)) error->all(FLERR, - "With overlay/pair no, and break yes, BPM bond styles requires special Coulomb weights = 1,1,1"); + "With overlay/pair no, and break yes, BPM bond styles requires special Coulomb " + "weights = 1,1,1"); - 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)); - delete[] id_fix_dummy; - id_fix_dummy = nullptr; + if (id_fix_dummy_special && break_flag) { + id_fix_update_special_bonds = utils::strdup("BPM_UPDATE_SPECIAL_BONDS"); + auto newfix = modify->replace_fix( + id_fix_dummy_special, + fmt::format("{} all UPDATE_SPECIAL_BONDS", id_fix_update_special_bonds), 1); + fix_update_special_bonds = dynamic_cast(newfix); + delete[] id_fix_dummy_special; + id_fix_dummy_special = nullptr; } } // special 1-3 and 1-4 weights must be 1 to prevent building 1-3 and 1-4 special bond lists - if (force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 || force->special_coul[2] != 1.0 || - force->special_coul[3] != 1.0) + if (force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 || + force->special_coul[2] != 1.0 || force->special_coul[3] != 1.0) error->all(FLERR, "Bond style bpm requires 1-3 and 1-4 special weights of 1.0"); } @@ -170,6 +176,17 @@ void BondBPM::init_style() hybrid_flag = 0; for (int i = 1; i <= atom->nbondtypes; i++) if (!setflag[i]) hybrid_flag = 1; + + // Set up necessary history fix + if (!fix_bond_history) { + auto newfix = modify->replace_fix( + id_fix_dummy_history, + fmt::format("{} all BOND_HISTORY {} {}", id_fix_bond_history, update_flag, nhistory), 1); + fix_bond_history = dynamic_cast(newfix); + delete[] id_fix_dummy_history; + id_fix_dummy_history = nullptr; + } + fix_bond_history->setflag = setflag; } @@ -208,20 +225,21 @@ void BondBPM::settings(int narg, char **arg) 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; + property_atom_flag = 1; } else if (strcmp(arg[iarg], "y/ref") == 0) { pack_choice[nvalues++] = &BondBPM::pack_y_ref; - prop_atom_flag = 1; + property_atom_flag = 1; } else if (strcmp(arg[iarg], "z/ref") == 0) { pack_choice[nvalues++] = &BondBPM::pack_z_ref; - prop_atom_flag = 1; + property_atom_flag = 1; } else { break; } iarg++; } } else if (strcmp(arg[iarg], "overlay/pair") == 0) { - if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command, missing option for overlay/pair"); + if (iarg + 1 > narg) + error->all(FLERR, "Illegal bond bpm command, missing option for overlay/pair"); overlay_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg += 2; } else if (strcmp(arg[iarg], "break") == 0) { @@ -249,17 +267,17 @@ void BondBPM::settings(int narg, char **arg) // Use property/atom to save reference positions as it can transfer to ghost atoms // This won't work for instances where bonds are added (e.g. fix pour) but in those cases // a reference state isn't well defined - if (prop_atom_flag == 1) { + if (property_atom_flag == 1) { - id_fix_prop_atom = utils::strdup("BPM_property_atom"); + id_fix_property_atom = utils::strdup("BPM_property_atom"); char *x_ref_id = utils::strdup("BPM_X_REF"); char *y_ref_id = utils::strdup("BPM_Y_REF"); char *z_ref_id = utils::strdup("BPM_Z_REF"); - ifix = modify->get_fix_by_id(id_fix_prop_atom); + ifix = modify->get_fix_by_id(id_fix_property_atom); if (!ifix) ifix = modify->add_fix(fmt::format("{} all property/atom d_{} d_{} d_{} ghost yes", - id_fix_prop_atom, x_ref_id, y_ref_id, z_ref_id)); + id_fix_property_atom, x_ref_id, y_ref_id, z_ref_id)); int type_flag; int col_flag; @@ -290,10 +308,12 @@ void BondBPM::settings(int narg, char **arg) // Set up necessary history fix if (!fix_bond_history) { - fix_bond_history = dynamic_cast(modify->replace_fix( - id_fix_dummy2, fmt::format("{} all BOND_HISTORY {} {}", id_fix_bond_history, update_flag, nhistory), 1)); - delete[] id_fix_dummy2; - id_fix_dummy2 = nullptr; + auto newfix = modify->replace_fix( + id_fix_dummy_history, + fmt::format("{} all BOND_HISTORY {} {}", id_fix_bond_history, update_flag, nhistory), 1); + fix_bond_history = dynamic_cast(newfix); + delete[] id_fix_dummy_history; + id_fix_dummy_history = nullptr; } } @@ -376,8 +396,7 @@ 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 (!break_flag) error->one(FLERR, "BPM bond broke with break no option"); int nlocal = atom->nlocal; if (fix_store_local) { @@ -415,8 +434,8 @@ void BondBPM::process_broken(int i, int j) n = num_bond[i]; bond_type[i][m] = bond_type[i][n - 1]; bond_atom[i][m] = bond_atom[i][n - 1]; - for (auto &ihistory: histories) { - auto fix_bond_history2 = dynamic_cast (ihistory); + for (auto &ihistory : histories) { + auto fix_bond_history2 = dynamic_cast(ihistory); fix_bond_history2->shift_history(i, m, n - 1); fix_bond_history2->delete_history(i, n - 1); } @@ -432,8 +451,8 @@ void BondBPM::process_broken(int i, int j) n = num_bond[j]; bond_type[j][m] = bond_type[j][n - 1]; bond_atom[j][m] = bond_atom[j][n - 1]; - for (auto &ihistory: histories) { - auto fix_bond_history2 = dynamic_cast (ihistory); + for (auto &ihistory : histories) { + auto fix_bond_history2 = dynamic_cast(ihistory); fix_bond_history2->shift_history(j, m, n - 1); fix_bond_history2->delete_history(j, n - 1); } diff --git a/src/BPM/bond_bpm.h b/src/BPM/bond_bpm.h index 28e4e7187e..42c44d61ce 100644 --- a/src/BPM/bond_bpm.h +++ b/src/BPM/bond_bpm.h @@ -42,9 +42,9 @@ class BondBPM : public Bond { std::vector leftover_iarg; - char *id_fix_dummy, *id_fix_dummy2; - char *id_fix_update, *id_fix_bond_history; - char *id_fix_store_local, *id_fix_prop_atom; + char *id_fix_dummy_special, *id_fix_dummy_history; + char *id_fix_update_special_bonds, *id_fix_bond_history; + char *id_fix_store_local, *id_fix_property_atom; class FixStoreLocal *fix_store_local; class FixBondHistory *fix_bond_history; class FixUpdateSpecialBonds *fix_update_special_bonds; @@ -54,7 +54,7 @@ class BondBPM : public Bond { FnPtrPack *pack_choice; // ptrs to pack functions double *output_data; - int prop_atom_flag, nvalues, overlay_flag, break_flag, ignore_special_flag; + int property_atom_flag, nvalues, overlay_flag, break_flag, ignore_special_flag; int index_x_ref, index_y_ref, index_z_ref; int n_histories; diff --git a/src/BPM/bond_bpm_rotational.cpp b/src/BPM/bond_bpm_rotational.cpp index 630f6ac906..2d5cb0652b 100644 --- a/src/BPM/bond_bpm_rotational.cpp +++ b/src/BPM/bond_bpm_rotational.cpp @@ -51,6 +51,7 @@ BondBPMRotational::BondBPMRotational(LAMMPS *_lmp) : partial_flag = 1; smooth_flag = 1; normalize_flag = 0; + writedata = 0; nhistory = 4; id_fix_bond_history = utils::strdup("HISTORY_BPM_ROTATIONAL"); diff --git a/src/BPM/bond_bpm_spring.cpp b/src/BPM/bond_bpm_spring.cpp index 2863bbf317..37437d5b60 100644 --- a/src/BPM/bond_bpm_spring.cpp +++ b/src/BPM/bond_bpm_spring.cpp @@ -33,17 +33,25 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ BondBPMSpring::BondBPMSpring(LAMMPS *_lmp) : - BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr) + BondBPM(_lmp), k(nullptr), av(nullptr), ecrit(nullptr), gamma(nullptr), + id_fix_property_bond(nullptr), vol_current(nullptr), dvol0(nullptr) { partial_flag = 1; smooth_flag = 1; normalize_flag = 0; + volume_flag = 0; + writedata = 0; nhistory = 1; id_fix_bond_history = utils::strdup("HISTORY_BPM_SPRING"); single_extra = 1; svector = new double[1]; + + nmax = 0; + + comm_forward = 0; + comm_reverse = 0; } /* ---------------------------------------------------------------------- */ @@ -51,13 +59,20 @@ BondBPMSpring::BondBPMSpring(LAMMPS *_lmp) : BondBPMSpring::~BondBPMSpring() { delete[] svector; + if (id_fix_property_bond && modify->nfix) { + modify->delete_fix(id_fix_property_bond); + delete[] id_fix_property_bond; + } if (allocated) { memory->destroy(setflag); memory->destroy(k); memory->destroy(ecrit); memory->destroy(gamma); + memory->destroy(av); } + + memory->destroy(dvol0); } /* ---------------------------------------------------------------------- @@ -134,18 +149,56 @@ void BondBPMSpring::store_data() void BondBPMSpring::compute(int eflag, int vflag) { + if (hybrid_flag) fix_bond_history->compress_history(); + + int i, bond_change_flag; + double *vol0, *vol; + + if (volume_flag) { + vol0 = atom->dvector[index_vol0]; + vol = atom->dvector[index_vol]; + + // grow + initialize dvol0 as necessary + if (nmax < atom->nmax) { + nmax = atom->nmax; + memory->create(dvol0, nmax, "bond/bpm/spring:dvol0"); + for (i = 0; i < nmax; i++) dvol0[i] = 0.0; + } + } if (!fix_bond_history->stored_flag) { fix_bond_history->stored_flag = true; store_data(); + + if (volume_flag) { + vol_current = vol0; + bond_change_flag = calculate_vol(); + + // zero dvol0, not needed since vol0 just calculated + for (i = 0; i < nmax; i++) dvol0[i] = 0.0; + } } - if (hybrid_flag) - fix_bond_history->compress_history(); + if (volume_flag) { + vol_current = vol; + bond_change_flag = calculate_vol(); + + // Update vol0 to account for any new bonds + if (bond_change_flag) { + update_vol0(); + + // forward new vol0 to ghosts before force calculation + vol_current = vol0; + comm->forward_comm(this); + + bond_change_flag = 0; + } + } int i1, i2, itmp, n, type; double delx, dely, delz, delvx, delvy, delvz; double e, rsq, r, r0, rinv, smooth, fbond, dot; + double vol_sum, vol0_sum, vol_temp; ev_init(eflag, vflag); @@ -157,6 +210,8 @@ void BondBPMSpring::compute(int eflag, int vflag) int nbondlist = neighbor->nbondlist; int nlocal = atom->nlocal; int newton_bond = force->newton_bond; + double dim = domain->dimension; + double invdim = 1.0 / dim; double **bondstore = fix_bond_history->bondstore; @@ -193,6 +248,15 @@ void BondBPMSpring::compute(int eflag, int vflag) if (fabs(e) > ecrit[type]) { bondlist[n][2] = 0; process_broken(i1, i2); + + if (volume_flag) { + bond_change_flag = 1; + vol_temp = r0 * r0; + if (dim == 3) vol_temp *= r0; + if (newton_bond || i1 < nlocal) dvol0[i1] -= vol_temp; + if (newton_bond || i2 < nlocal) dvol0[i2] -= vol_temp; + } + continue; } @@ -202,6 +266,12 @@ void BondBPMSpring::compute(int eflag, int vflag) else fbond = k[type] * (r0 - r); + if (volume_flag) { + vol_sum = vol[i1] + vol[i2]; + vol0_sum = vol0[i1] + vol0[i2]; + fbond += av[type] * (pow(vol_sum / vol0_sum, invdim) - 1.0 - e); + } + delvx = v[i1][0] - v[i2][0]; delvy = v[i1][1] - v[i2][1]; delvz = v[i1][2] - v[i2][2]; @@ -233,8 +303,77 @@ void BondBPMSpring::compute(int eflag, int vflag) if (evflag) ev_tally(i1, i2, nlocal, newton_bond, 0.0, fbond, delx, dely, delz); } - if (hybrid_flag) - fix_bond_history->uncompress_history(); + // Update vol0 to account for any broken bonds + if (volume_flag && bond_change_flag) update_vol0(); + + if (hybrid_flag) fix_bond_history->uncompress_history(); +} + +/* ---------------------------------------------------------------------- */ + +int BondBPMSpring::calculate_vol() +{ + int n, i1, i2; + double r0, delx, dely, delz, rsq, vol_temp; + + int nlocal = atom->nlocal; + int ntotal = nlocal + atom->nghost; + int newton_bond = force->newton_bond; + int dim = domain->dimension; + + double **x = atom->x; + int **bondlist = neighbor->bondlist; + int nbondlist = neighbor->nbondlist; + double **bondstore = fix_bond_history->bondstore; + + for (n = 0; n < ntotal; n++) vol_current[n] = 0.0; + + int bond_change_flag = 0; + + for (n = 0; n < nbondlist; n++) { + if (bondlist[n][2] <= 0) continue; + i1 = bondlist[n][0]; + i2 = bondlist[n][1]; + r0 = bondstore[n][0]; + + delx = x[i1][0] - x[i2][0]; + dely = x[i1][1] - x[i2][1]; + delz = x[i1][2] - x[i2][2]; + rsq = delx * delx + dely * dely + delz * delz; + + vol_temp = rsq; + if (dim == 3) vol_temp *= sqrt(rsq); + + if (newton_bond || i1 < nlocal) vol_current[i1] += vol_temp; + if (newton_bond || i2 < nlocal) vol_current[i2] += vol_temp; + + // If bond hasn't been set - increment dvol0 too to update vol0 + if (r0 < EPSILON || std::isnan(r0)) { + bond_change_flag = 1; + if (newton_bond || i1 < nlocal) dvol0[i1] += vol_temp; + if (newton_bond || i2 < nlocal) dvol0[i2] += vol_temp; + } + } + + if (newton_bond) comm->reverse_comm(this); + comm->forward_comm(this); + + return bond_change_flag; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPMSpring::update_vol0() +{ + // accumulate changes in vol0 from ghosts + vol_current = dvol0; + if (force->newton_bond) comm->reverse_comm(this); + + double *vol0 = atom->dvector[index_vol0]; + for (int i = 0; i < atom->nlocal; i++) vol0[i] += dvol0[i]; + + // zero dvol0 for next change + for (int i = 0; i < nmax; i++) dvol0[i] = 0.0; } /* ---------------------------------------------------------------------- */ @@ -247,6 +386,7 @@ void BondBPMSpring::allocate() memory->create(k, np1, "bond:k"); memory->create(ecrit, np1, "bond:ecrit"); memory->create(gamma, np1, "bond:gamma"); + memory->create(av, np1, "bond:av"); memory->create(setflag, np1, "bond:setflag"); for (int i = 1; i < np1; i++) setflag[i] = 0; @@ -258,7 +398,8 @@ void BondBPMSpring::allocate() void BondBPMSpring::coeff(int narg, char **arg) { - if (narg != 4) error->all(FLERR, "Incorrect args for bond coefficients"); + if ((!volume_flag && narg != 4) || (volume_flag && narg != 5)) + error->all(FLERR, "Incorrect args for bond coefficients"); if (!allocated) allocate(); int ilo, ihi; @@ -268,11 +409,15 @@ void BondBPMSpring::coeff(int narg, char **arg) double ecrit_one = utils::numeric(FLERR, arg[2], false, lmp); double gamma_one = utils::numeric(FLERR, arg[3], false, lmp); + double av_one = 0.0; + if (volume_flag) av_one = utils::numeric(FLERR, arg[4], false, lmp); + int count = 0; for (int i = ilo; i <= ihi; i++) { k[i] = k_one; ecrit[i] = ecrit_one; gamma[i] = gamma_one; + av[i] = av_one; setflag[i] = 1; count++; @@ -292,6 +437,16 @@ void BondBPMSpring::init_style() if (comm->ghost_velocity == 0) error->all(FLERR, "Bond bpm/spring requires ghost atoms store velocity"); + + if (volume_flag && !id_fix_property_bond) { + id_fix_property_bond = utils::strdup("BOND_BPM_SPRING_FIX_PROPERTY_ATOM"); + modify->add_fix(fmt::format("{} all property/atom d_vol d_vol0 ghost yes writedata no", + id_fix_property_bond)); + + int tmp1 = 0, tmp2 = 0; + index_vol = atom->find_custom("vol", tmp1, tmp2); + index_vol0 = atom->find_custom("vol0", tmp1, tmp2); + } } /* ---------------------------------------------------------------------- */ @@ -308,8 +463,23 @@ void BondBPMSpring::settings(int narg, char **arg) smooth_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp); i += 1; } else if (strcmp(arg[iarg], "normalize") == 0) { - if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command, missing option for normalize"); + if (iarg + 1 > narg) + error->all(FLERR, "Illegal bond bpm command, missing option for normalize"); normalize_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp); + i += 1; + } else if (strcmp(arg[iarg], "volume/factor") == 0) { + if (iarg + 1 > narg) + error->all(FLERR, "Illegal bond bpm command, missing option for volume/factor"); + volume_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp); + + if (volume_flag) { + comm_forward = 1; + comm_reverse = 1; + } else { + comm_forward = 0; + comm_reverse = 0; + } + i += 1; } else { error->all(FLERR, "Illegal bond bpm command, invalid argument {}", arg[iarg]); @@ -329,6 +499,7 @@ void BondBPMSpring::write_restart(FILE *fp) fwrite(&k[1], sizeof(double), atom->nbondtypes, fp); fwrite(&ecrit[1], sizeof(double), atom->nbondtypes, fp); fwrite(&gamma[1], sizeof(double), atom->nbondtypes, fp); + fwrite(&av[1], sizeof(double), atom->nbondtypes, fp); } /* ---------------------------------------------------------------------- @@ -345,10 +516,12 @@ void BondBPMSpring::read_restart(FILE *fp) utils::sfread(FLERR, &k[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); utils::sfread(FLERR, &ecrit[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); utils::sfread(FLERR, &gamma[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); + utils::sfread(FLERR, &av[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); } MPI_Bcast(&k[1], atom->nbondtypes, MPI_DOUBLE, 0, world); MPI_Bcast(&ecrit[1], atom->nbondtypes, MPI_DOUBLE, 0, world); MPI_Bcast(&gamma[1], atom->nbondtypes, MPI_DOUBLE, 0, world); + MPI_Bcast(&av[1], atom->nbondtypes, MPI_DOUBLE, 0, world); for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1; } @@ -360,6 +533,8 @@ void BondBPMSpring::read_restart(FILE *fp) void BondBPMSpring::write_restart_settings(FILE *fp) { fwrite(&smooth_flag, sizeof(int), 1, fp); + fwrite(&normalize_flag, sizeof(int), 1, fp); + fwrite(&volume_flag, sizeof(int), 1, fp); } /* ---------------------------------------------------------------------- @@ -368,8 +543,14 @@ void BondBPMSpring::write_restart_settings(FILE *fp) void BondBPMSpring::read_restart_settings(FILE *fp) { - if (comm->me == 0) utils::sfread(FLERR, &smooth_flag, sizeof(int), 1, fp, nullptr, error); + if (comm->me == 0) { + utils::sfread(FLERR, &smooth_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &normalize_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &volume_flag, sizeof(int), 1, fp, nullptr, error); + } MPI_Bcast(&smooth_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&normalize_flag, 1, MPI_INT, 0, world); + MPI_Bcast(&volume_flag, 1, MPI_INT, 0, world); } /* ---------------------------------------------------------------------- */ @@ -385,12 +566,22 @@ double BondBPMSpring::single(int type, double rsq, int i, int j, double &fforce) double r = sqrt(rsq); double rinv = 1.0 / r; + double e = (r - r0) / r0; if (normalize_flag) - fforce = k[type] * (r0 - r) / r0; + fforce = -k[type] * e; else fforce = k[type] * (r0 - r); + if (volume_flag) { + double invdim = 1.0 / domain->dimension; + double *vol0 = atom->dvector[index_vol0]; + double *vol = atom->dvector[index_vol]; + double vol_sum = vol[i] + vol[j]; + double vol0_sum = vol0[i] + vol0[j]; + fforce += av[type] * (pow(vol_sum / vol0_sum, invdim) - 1.0 - e); + } + double **x = atom->x; double **v = atom->v; double delx = x[i][0] - x[j][0]; @@ -418,3 +609,49 @@ double BondBPMSpring::single(int type, double rsq, int i, int j, double &fforce) return 0.0; } + +/* ---------------------------------------------------------------------- */ + +int BondBPMSpring::pack_reverse_comm(int n, int first, double *buf) +{ + int i, m, last; + m = 0; + last = first + n; + for (i = first; i < last; i++) buf[m++] = vol_current[i]; + return m; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPMSpring::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i, j, m; + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + vol_current[j] += buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int BondBPMSpring::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) +{ + int i, j, m; + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = vol_current[j]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPMSpring::unpack_forward_comm(int n, int first, double *buf) +{ + int i, m, last; + m = 0; + last = first + n; + for (i = first; i < last; i++) vol_current[i] = buf[m++]; +} diff --git a/src/BPM/bond_bpm_spring.h b/src/BPM/bond_bpm_spring.h index 93f4b49a26..8e5a80229d 100644 --- a/src/BPM/bond_bpm_spring.h +++ b/src/BPM/bond_bpm_spring.h @@ -37,14 +37,24 @@ class BondBPMSpring : public BondBPM { void write_restart_settings(FILE *) override; void read_restart_settings(FILE *) override; double single(int, double, int, int, double &) override; + int pack_forward_comm(int, int *, double *, int, int *) override; + void unpack_forward_comm(int, int, double *) override; + int pack_reverse_comm(int, int, double *) override; + void unpack_reverse_comm(int, int *, double *) override; protected: - double *k, *ecrit, *gamma; - int smooth_flag, normalize_flag; + double *k, *av, *ecrit, *gamma; + int smooth_flag, normalize_flag, volume_flag; + + int index_vol, index_vol0, nmax; + char *id_fix_property_bond; + double *vol_current, *dvol0; void allocate(); void store_data(); double store_bond(int, int, int); + int calculate_vol(); + void update_vol0(); }; } // namespace LAMMPS_NS diff --git a/src/BPM/pair_bpm_spring.cpp b/src/BPM/pair_bpm_spring.cpp index 3407e0274b..bf58b9e515 100644 --- a/src/BPM/pair_bpm_spring.cpp +++ b/src/BPM/pair_bpm_spring.cpp @@ -27,9 +27,10 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -PairBPMSpring::PairBPMSpring(LAMMPS *_lmp) : Pair(_lmp), k(nullptr), cut(nullptr), gamma(nullptr) +PairBPMSpring::PairBPMSpring(LAMMPS *_lmp) : Pair(_lmp), k(nullptr), ka(nullptr), cut(nullptr), gamma(nullptr) { writedata = 1; + anharmonic_flag = 0; } /* ---------------------------------------------------------------------- */ @@ -41,6 +42,7 @@ PairBPMSpring::~PairBPMSpring() memory->destroy(cutsq); memory->destroy(k); + memory->destroy(ka); memory->destroy(cut); memory->destroy(gamma); } @@ -51,7 +53,7 @@ PairBPMSpring::~PairBPMSpring() void PairBPMSpring::compute(int eflag, int vflag) { int i, j, ii, jj, inum, jnum, itype, jtype; - double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair; + double xtmp, ytmp, ztmp, delx, dely, delz, dr, evdwl, fpair; double r, rsq, rinv, factor_lj; int *ilist, *jlist, *numneigh, **firstneigh; double vxtmp, vytmp, vztmp, delvx, delvy, delvz, dot, smooth; @@ -107,7 +109,11 @@ void PairBPMSpring::compute(int eflag, int vflag) r = sqrt(rsq); rinv = 1.0 / r; - fpair = k[itype][jtype] * (cut[itype][jtype] - r); + dr = r - cut[itype][jtype]; + + fpair = -k[itype][jtype] * dr; + if (anharmonic_flag) + fpair += -ka[itype][jtype] * dr * dr * dr; smooth = rsq / cutsq[itype][jtype]; smooth *= smooth; @@ -156,6 +162,7 @@ void PairBPMSpring::allocate() memory->create(cutsq, np1, np1, "pair:cutsq"); memory->create(k, np1, np1, "pair:k"); + memory->create(ka, np1, np1, "pair:ka"); memory->create(cut, np1, np1, "pair:cut"); memory->create(gamma, np1, np1, "pair:gamma"); } @@ -164,9 +171,17 @@ void PairBPMSpring::allocate() global settings ------------------------------------------------------------------------- */ -void PairBPMSpring::settings(int narg, char ** /*arg*/) +void PairBPMSpring::settings(int narg, char ** arg) { - if (narg != 0) error->all(FLERR, "Illegal pair_style command"); + int iarg = 0; + while (iarg < narg) { + if (strcmp(arg[iarg], "anharmonic") == 0) { + if (iarg + 1 >= narg) + utils::missing_cmd_args(FLERR, "pair_coeff bpm/spring anharmonic", error); + anharmonic_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else error->all(FLERR, "Illegal pair_style command {}", arg[iarg]); + } } /* ---------------------------------------------------------------------- @@ -175,7 +190,8 @@ void PairBPMSpring::settings(int narg, char ** /*arg*/) void PairBPMSpring::coeff(int narg, char **arg) { - if (narg != 5) error->all(FLERR, "Incorrect args for pair coefficients"); + if ((!anharmonic_flag && narg != 5) || (anharmonic_flag && narg != 6)) + error->all(FLERR, "Incorrect args for pair coefficients"); if (!allocated) allocate(); int ilo, ihi, jlo, jhi; @@ -186,6 +202,10 @@ void PairBPMSpring::coeff(int narg, char **arg) double cut_one = utils::numeric(FLERR, arg[3], false, lmp); double gamma_one = utils::numeric(FLERR, arg[4], false, lmp); + double ka_one = 0.0; + if (anharmonic_flag) + ka_one = utils::numeric(FLERR, arg[5], false, lmp); + if (cut_one <= 0.0) error->all(FLERR, "Incorrect args for pair coefficients"); int count = 0; @@ -194,6 +214,7 @@ void PairBPMSpring::coeff(int narg, char **arg) k[i][j] = k_one; cut[i][j] = cut_one; gamma[i][j] = gamma_one; + ka[i][j] = ka_one; setflag[i][j] = 1; count++; @@ -230,6 +251,7 @@ double PairBPMSpring::init_one(int i, int j) cut[j][i] = cut[i][j]; k[j][i] = k[i][j]; gamma[j][i] = gamma[i][j]; + ka[j][i] = ka[i][j]; return cut[i][j]; } @@ -250,6 +272,7 @@ void PairBPMSpring::write_restart(FILE *fp) fwrite(&k[i][j], sizeof(double), 1, fp); fwrite(&cut[i][j], sizeof(double), 1, fp); fwrite(&gamma[i][j], sizeof(double), 1, fp); + fwrite(&ka[i][j], sizeof(double), 1, fp); } } } @@ -274,22 +297,50 @@ void PairBPMSpring::read_restart(FILE *fp) utils::sfread(FLERR, &k[i][j], sizeof(double), 1, fp, nullptr, error); utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error); utils::sfread(FLERR, &gamma[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &ka[i][j], sizeof(double), 1, fp, nullptr, error); } MPI_Bcast(&k[i][j], 1, MPI_DOUBLE, 0, world); MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world); MPI_Bcast(&gamma[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&ka[i][j], 1, MPI_DOUBLE, 0, world); } } } + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairBPMSpring::write_restart_settings(FILE *fp) +{ + fwrite(&anharmonic_flag, sizeof(int), 1, fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairBPMSpring::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) + utils::sfread(FLERR, &anharmonic_flag, sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&anharmonic_flag, 1, MPI_INT, 0, world); +} + /* ---------------------------------------------------------------------- proc 0 writes to data file ------------------------------------------------------------------------- */ void PairBPMSpring::write_data(FILE *fp) { - for (int i = 1; i <= atom->ntypes; i++) - fprintf(fp, "%d %g %g %g\n", i, k[i][i], cut[i][i], gamma[i][i]); + if (anharmonic_flag) { + for (int i = 1; i <= atom->ntypes; i++) + fprintf(fp, "%d %g %g %g %g\n", i, k[i][i], cut[i][i], gamma[i][i], ka[i][i]); + } else { + for (int i = 1; i <= atom->ntypes; i++) + fprintf(fp, "%d %g %g %g\n", i, k[i][i], cut[i][i], gamma[i][i]); + } } /* ---------------------------------------------------------------------- @@ -298,9 +349,15 @@ void PairBPMSpring::write_data(FILE *fp) void PairBPMSpring::write_data_all(FILE *fp) { - for (int i = 1; i <= atom->ntypes; i++) - for (int j = i; j <= atom->ntypes; j++) - fprintf(fp, "%d %d %g %g %g\n", i, j, k[i][j], cut[i][j], gamma[i][j]); + if (anharmonic_flag) { + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i; j <= atom->ntypes; j++) + fprintf(fp, "%d %d %g %g %g %g\n", i, j, k[i][j], cut[i][j], gamma[i][j], ka[i][j]); + } else { + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i; j <= atom->ntypes; j++) + fprintf(fp, "%d %d %g %g %g\n", i, j, k[i][j], cut[i][j], gamma[i][j]); + } } /* ---------------------------------------------------------------------- */ @@ -308,7 +365,7 @@ void PairBPMSpring::write_data_all(FILE *fp) double PairBPMSpring::single(int i, int j, int itype, int jtype, double rsq, double /*factor_coul*/, double factor_lj, double &fforce) { - double fpair, r, rinv; + double fpair, r, rinv, dr; double delx, dely, delz, delvx, delvy, delvz, dot, smooth; if (rsq > cutsq[itype][jtype]) return 0.0; @@ -319,7 +376,10 @@ double PairBPMSpring::single(int i, int j, int itype, int jtype, double rsq, dou r = sqrt(rsq); rinv = 1.0 / r; - fpair = k[itype][jtype] * (cut[itype][jtype] - r); + dr = r - cut[itype][jtype]; + fpair = -k[itype][jtype] * dr; + if (anharmonic_flag) + fpair += -ka[itype][jtype] * dr * dr * dr; smooth = rsq / cutsq[itype][jtype]; smooth *= smooth; diff --git a/src/BPM/pair_bpm_spring.h b/src/BPM/pair_bpm_spring.h index c10e4a3400..515317e045 100644 --- a/src/BPM/pair_bpm_spring.h +++ b/src/BPM/pair_bpm_spring.h @@ -29,18 +29,21 @@ class PairBPMSpring : public Pair { PairBPMSpring(class LAMMPS *); ~PairBPMSpring() override; void compute(int, int) override; - void settings(int, char **) override; void coeff(int, char **) override; + void settings(int, char **) override; void init_style() override; double init_one(int, int) override; void write_restart(FILE *) override; void read_restart(FILE *) override; + void write_restart_settings(FILE *) override; + void read_restart_settings(FILE *) override; void write_data(FILE *) override; void write_data_all(FILE *) override; double single(int, int, int, int, double, double, double, double &) override; protected: - double **k, **cut, **gamma; + int anharmonic_flag; + double **k, **ka, **cut, **gamma; void allocate(); }; diff --git a/src/modify.cpp b/src/modify.cpp index 1a9c056027..89cc1d0643 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -825,16 +825,16 @@ Fix *Modify::add_fix(int narg, char **arg, int trysuffix) // clang-format off const char *exceptions[] = - {"GPU", "OMP", "INTEL", "property/atom", "cmap", "cmap3", "rx", - "deprecated", "STORE/KIM", "amoeba/pitorsion", "amoeba/bitorsion", - nullptr}; + {"GPU", "OMP", "INTEL", "property/atom", "cmap", "cmap3", "rx", "deprecated", "STORE/KIM", + "amoeba/pitorsion", "amoeba/bitorsion", "DUMMY", nullptr}; // clang-format on if (domain->box_exist == 0) { int m; for (m = 0; exceptions[m] != nullptr; m++) if (strcmp(arg[2], exceptions[m]) == 0) break; - if (exceptions[m] == nullptr) error->all(FLERR, "Fix command before simulation box is defined"); + if (exceptions[m] == nullptr) + error->all(FLERR, "Fix {} command before simulation box is defined", arg[2]); } // check group ID