Merge pull request #4359 from jtclemm/BPM

Adding new force options to BPM package

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 <bond_bpm_spring>`, only applies
 pairwise, central body forces. Point particles must have :doc:`bond atom
 style <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
 <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 <fix_nve_asphere>`.
 
 In addition to bond styles, a new pair style :doc:`pair bpm/spring
 <pair_bpm_spring>` was added to accompany the bpm/spring bond
-style. This pair style is simply a hookean repulsion with similar
+style. By default, this pair style is simply a hookean repulsion with
-velocity damping as its sister bond style.
+similar velocity damping as its sister bond style, but optional
+arguments can be used to modify the force.
 
 ----------
 

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 <bond_coeff>` command as in the example above, or in
-the data file or restart files read by the :doc:`read_data
-<read_data>` or :doc:`read_restart <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) <multibody-Clemmer>`,
+
+.. 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 <bond_coeff>` command as in the example above, or in
+the data file or restart files read by the :doc:`read_data
+<read_data>` or :doc:`read_restart <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).

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.
+where :math:`k` is a stiffness, :math:`r_c` is the cutoff
-An additional damping force is also applied to interacting
+length, and :math:`k_a` is an optional anharmonic cubic prefactor
-particles. The force is proportional to the difference in the
+that can be enabled using the *anharmonic* keyword. The anharmonic
-normal velocity of particles
+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*
+

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            2 all local 100 brokenDump f_brkbond[1] f_brkbond[2] f_brkbond[3]
-dump_modify     2 header no
+#dump_modify     2 header no
 
 run             7500

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            2 all local 100 brokenDump f_brkbond[1] f_brkbond[2] f_brkbond[3]
-dump_modify     2 header no
+#dump_modify     2 header no
 
 run             7500

src/BPM/bond_bpm.cpp

@@ -48,16 +48,17 @@ static const char cite_bpm[] =
 /* ---------------------------------------------------------------------- */
 
 BondBPM::BondBPM(LAMMPS *_lmp) :
-    Bond(_lmp), id_fix_dummy(nullptr), id_fix_dummy2(nullptr), id_fix_update(nullptr),
+    Bond(_lmp), id_fix_dummy_special(nullptr), id_fix_dummy_history(nullptr),
-    id_fix_bond_history(nullptr), id_fix_store_local(nullptr), id_fix_prop_atom(nullptr),
+    id_fix_update_special_bonds(nullptr), id_fix_bond_history(nullptr), id_fix_store_local(nullptr),
-    fix_store_local(nullptr), fix_bond_history(nullptr), fix_update_special_bonds(nullptr),
+    id_fix_property_atom(nullptr), fix_store_local(nullptr), fix_bond_history(nullptr),
-    pack_choice(nullptr), output_data(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));
+  id_fix_dummy_special = utils::strdup(fmt::format("BPM_DUMMY_SPECIAL_{}", instance_total));
-  modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy));
+  modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy_special));
 
-  id_fix_dummy2 = utils::strdup(fmt::format("BPM_DUMMY2_{}", instance_total));
+  id_fix_dummy_history = utils::strdup(fmt::format("BPM_DUMMY_HISTORY_{}", instance_total));
-  modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy2));
+  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_dummy_special) modify->delete_fix(id_fix_dummy_special);
-  if (id_fix_dummy2) modify->delete_fix(id_fix_dummy2);
+  if (id_fix_dummy_history) modify->delete_fix(id_fix_dummy_history);
-  if (id_fix_update) modify->delete_fix(id_fix_update);
+  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_update_special_bonds;
-  delete[] id_fix_dummy2;
-  delete[] id_fix_update;
   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 ||
+      if (force->special_lj[1] != 1.0 || force->special_lj[2] != 1.0 ||
-          force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 || force->special_coul[3] != 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");
+                   "With overlay/pair yes, BPM bond styles require a value of 1.0 for all "
-      if (id_fix_update) {
+                   "special_bonds weights");
-        modify->delete_fix(id_fix_update);
+      if (id_fix_update_special_bonds) {
-        delete[] id_fix_update;
+        modify->delete_fix(id_fix_update_special_bonds);
-        id_fix_update = nullptr;
+        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 ||
+      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) {
+      if (id_fix_dummy_special && break_flag) {
-        id_fix_update = utils::strdup("BPM_UPDATE_SPECIAL_BONDS");
+        id_fix_update_special_bonds = utils::strdup("BPM_UPDATE_SPECIAL_BONDS");
-        fix_update_special_bonds = dynamic_cast<FixUpdateSpecialBonds *>(modify->replace_fix(
+        auto newfix = modify->replace_fix(
-            id_fix_dummy, fmt::format("{} all UPDATE_SPECIAL_BONDS", id_fix_update), 1));
+            id_fix_dummy_special,
-        delete[] id_fix_dummy;
+            fmt::format("{} all UPDATE_SPECIAL_BONDS", id_fix_update_special_bonds), 1);
-        id_fix_dummy = nullptr;
+        fix_update_special_bonds = dynamic_cast<FixUpdateSpecialBonds *>(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 ||
+    if (force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 ||
-        force->special_coul[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<FixBondHistory *>(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<FixBondHistory *>(modify->replace_fix(
+    auto newfix = modify->replace_fix(
-        id_fix_dummy2, fmt::format("{} all BOND_HISTORY {} {}", id_fix_bond_history, update_flag, nhistory), 1));
+        id_fix_dummy_history,
-    delete[] id_fix_dummy2;
+        fmt::format("{} all BOND_HISTORY {} {}", id_fix_bond_history, update_flag, nhistory), 1);
-    id_fix_dummy2 = nullptr;
+    fix_bond_history = dynamic_cast<FixBondHistory *>(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)
+  if (!break_flag) error->one(FLERR, "BPM bond broke with break no option");
-    error->one(FLERR, "BPM bond broke with break no option");
 
   int nlocal = atom->nlocal;
   if (fix_store_local) {

src/BPM/bond_bpm.h

@@ -42,9 +42,9 @@ class BondBPM : public Bond {
 
   std::vector<int> leftover_iarg;
 
-  char *id_fix_dummy, *id_fix_dummy2;
+  char *id_fix_dummy_special, *id_fix_dummy_history;
-  char *id_fix_update, *id_fix_bond_history;
+  char *id_fix_update_special_bonds, *id_fix_bond_history;
-  char *id_fix_store_local, *id_fix_prop_atom;
+  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;

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");

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)
+  if (volume_flag) {
-    fix_bond_history->compress_history();
+    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)
+  // Update vol0 to account for any broken bonds
-    fix_bond_history->uncompress_history();
+  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++];
+}

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;
+  double *k, *av, *ecrit, *gamma;
-  int smooth_flag, normalize_flag;
+  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

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,23 +297,51 @@ 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)
 {
+  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]);
   }
+}
 
 /* ----------------------------------------------------------------------
    proc 0 writes all pairs to data file
@@ -298,17 +349,23 @@ void PairBPMSpring::write_data(FILE *fp)
 
 void PairBPMSpring::write_data_all(FILE *fp)
 {
+  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]);
   }
+}
 
 /* ---------------------------------------------------------------------- */
 
 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;

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();
 };

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",
+    {"GPU", "OMP", "INTEL", "property/atom", "cmap", "cmap3", "rx", "deprecated", "STORE/KIM",
-     "deprecated", "STORE/KIM", "amoeba/pitorsion", "amoeba/bitorsion",
+     "amoeba/pitorsion", "amoeba/bitorsion", "DUMMY", nullptr};
-     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