Adding compatability with MC fixes, set_array to property/atom, faster update/special/bonds, single methods, and misc small changes

doc/src/Howto_bpm.rst

@@ -109,13 +109,14 @@ velocity damping as its sister bond style.
 
 ----------
 
-While LAMMPS has many untilites to create and delete bonds, only a
+While LAMMPS has many untilites to create and delete bonds, the
-select few are compatible with BPM bond styles. They include:
+following are currently compatible with BPM bond styles:
 
 * :doc:`create_bonds <create_bonds>`
 * :doc:`delete_bonds <delete_bonds>`
 * :doc:`fix bond/create <fix_bond_create>`
 * :doc:`fix bond/break <fix_bond_break>`
+* :doc:`fix bond/swap <fix_bond_swap>`
 
 Note :doc:`bond_create <bond_create>` requires certain special_bonds settings.
 To subtract pair interactions, one will need to switch between different

doc/src/Howto_broken_bonds.rst

@@ -31,7 +31,9 @@ page.
 
 Bonds can also be broken by fixes which change bond topology, including
 :doc:`fix bond/break <fix_bond_break>` and
-:doc:`fix bond/react <fix_bond_react>`.
+:doc:`fix bond/react <fix_bond_react>`. These fixes will automatically
+trigger a rebuild of the neighbor list and update special bond data structures
+when bonds are broken.
 
 Note that when bonds are dumped to a file via the :doc:`dump local <dump>` command, bonds with type 0 are not included.  The
 :doc:`delete_bonds <delete_bonds>` command can also be used to query the

doc/src/bond_bpm_rotational.rst

@@ -150,13 +150,18 @@ status of broken bonds or permanently delete them, e.g.:
 
 ----------
 
-Restart
+Restart and other info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
 This bond style writes the reference state of each bond to
 :doc:`binary restart files <restart>`. Loading a restart file will
 properly resume bonds.
 
+The single() function of these pair styles returns 0.0 for the energy
+of a pairwise interaction, since energy is not conserved in these
+dissipative potentials.  It also returns only the normal component of
+the pairwise interaction force.
+
 Restrictions
 """"""""""""
 

doc/src/bond_bpm_spring.rst

@@ -113,13 +113,17 @@ query the status of broken bonds or permanently delete them, e.g.:
 
 ----------
 
-Restart
+Restart and other info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
 This bond style writes the reference state of each bond to
 :doc:`binary restart files <restart>`. Loading a restart
 file will properly resume bonds.
 
+The single() function of these pair styles returns 0.0 for the energy
+of a pairwise interaction, since energy is not conserved in these
+dissipative potentials.
+
 Restrictions
 """"""""""""
 

examples/bpm/impact/in.bpm.impact.rotational

@@ -21,10 +21,7 @@ displace_atoms  all random 0.1 0.1 0.1 134598738
 
 neighbor        1.0 bin
 pair_style      gran/hooke/history 1.0 NULL 0.5 NULL 0.1 1
-bond_style      bpm/rotational store/local 2 time id1 id2
 pair_coeff      1 1
-bond_coeff      1 1.0 0.2 0.02 0.02 0.05 0.01 0.01 0.01 0.1 0.02 0.002 0.002
-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
@@ -35,6 +32,10 @@ create_bonds    many projectile projectile 2 0.0 1.5
 neighbor        0.3 bin
 special_bonds   lj 0.0 1.0 1.0 coul 1.0 1.0 1.0
 
+bond_style      bpm/rotational store/local 2 time id1 id2
+bond_coeff      1 1.0 0.2 0.02 0.02 0.05 0.01 0.01 0.01 0.1 0.02 0.002 0.002
+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
+
 velocity        projectile set -0.05 0.0 0.0
 compute         nbond all nbond/atom
 compute         tbond all reduce sum c_nbond

examples/bpm/impact/in.bpm.impact.spring

@@ -23,10 +23,7 @@ mass            1 1.0
 
 neighbor        1.0 bin
 pair_style      bpm/spring
-bond_style      bpm/spring store/local 2 time id1 id2
 pair_coeff      1 1 1.0 1.0 1.0
-bond_coeff      1 1.0 0.04 1.0
-bond_coeff      2 1.0 0.20 1.0
 
 fix             1 all nve
 fix             2 all store/local 100 3
@@ -37,6 +34,10 @@ create_bonds    many projectile projectile 2 0.0 1.5
 neighbor        0.3 bin
 special_bonds   lj 0.0 1.0 1.0 coul 1.0 1.0 1.0
 
+bond_style      bpm/spring store/local 2 time id1 id2
+bond_coeff      1 1.0 0.04 1.0
+bond_coeff      2 1.0 0.20 1.0
+
 velocity        projectile set -0.05 0.0 0.0
 compute         nbond all nbond/atom
 compute         tbond all reduce sum c_nbond

src/.gitignore

@@ -265,8 +265,6 @@
 /compute_nbond_atom.h
 /fix_nve_sphere_bpm.cpp
 /fix_nve_sphere_bpm.h
-/fix_update_special_bonds.cpp
-/fix_update_special_bonds.h
 /pair_bpm_spring.cpp
 /pair_bpm_spring.h
 

src/BPM/bond_bpm.cpp

@@ -25,6 +25,9 @@
 #include "modify.h"
 #include "update.h"
 
+#include <string>
+#include <vector>
+
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
@@ -49,7 +52,7 @@ BondBPM::BondBPM(LAMMPS *lmp) : Bond(lmp)
   // 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");
+  id_fix_dummy = utils::strdup("BPM_DUMMY" + std::to_string(instance_me));
   modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy));
 }
 
@@ -58,14 +61,16 @@ BondBPM::BondBPM(LAMMPS *lmp) : Bond(lmp)
 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);
+  if (id_fix_store_local) modify->delete_fix(id_fix_store_local);
+  if (id_fix_prop_atom) modify->delete_fix(id_fix_prop_atom);
 
   delete [] id_fix_dummy;
   delete [] id_fix_update;
+  delete [] id_fix_store_local;
+  delete [] id_fix_prop_atom;
 
   memory->destroy(output_data);
 }
@@ -110,7 +115,7 @@ void BondBPM::init_style()
       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");
+      id_fix_update = utils::strdup("BPM_update_special_bonds" + std::to_string(instance_me));
       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;
@@ -141,11 +146,14 @@ void BondBPM::init_style()
 
 void BondBPM::settings(int narg, char **arg)
 {
+  leftover_args.clear();
+
+  int local_freq;
   int iarg = 0;
   while (iarg < narg) {
     if (strcmp(arg[iarg], "store/local") == 0) {
-      id_fix_store_local = utils::strdup(arg[iarg+1]);
       nvalues = 0;
+      local_freq = utils::inumeric(FLERR, arg[iarg+1], false, lmp);
       pack_choice = new FnPtrPack[narg - iarg - 1];
       iarg += 2;
       while (iarg < narg) {
@@ -179,31 +187,44 @@ void BondBPM::settings(int narg, char **arg)
       overlay_flag = 1;
       iarg ++;
     } else {
-      error->all(FLERR, "Illegal pair_style command");
+      leftover_args.push_back(iarg);
+      iarg ++;
     }
   }
 
-  if (id_fix_store_local) {
+  if (nvalues != 0 && !id_fix_store_local) {
+    //Todo, assign ID and create fix id_fix_store_local = utils::strdup(arg[iarg+1]);
+
     if (nvalues == 0) error->all(FLERR,
         "Bond style bpm/rotational must include at least one value to output");
     memory->create(output_data, nvalues, "bond/bpm:output_data");
 
     // Use store property 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) {
 
-      id_fix_prop_atom = utils::strdup("BPM_property_atom");
+      id_fix_prop_atom = utils::strdup("BPM_property_atom" + std::to_string(instance_me));
       int ifix = modify->find_fix(id_fix_prop_atom);
+
+      char *x_ref_id = utils::strdup("BPM_X_REF" + std::to_string(instance_me));
+      char *y_ref_id = utils::strdup("BPM_Y_REF" + std::to_string(instance_me));
+      char *z_ref_id = utils::strdup("BPM_Z_REF" + std::to_string(instance_me));
       if (ifix < 0) {
-        modify->add_fix(fmt::format("{} all property/atom "
+        modify->add_fix(fmt::format("{} all property/atom {} {} {} ghost yes",
-            "d_BPM_X_REF d_BPM_Y_REF d_BPM_Z_REF ghost yes", id_fix_prop_atom));
+          id_fix_prop_atom, x_ref_id, y_ref_id, z_ref_id));
         ifix = modify->find_fix(id_fix_prop_atom);
       }
 
       int type_flag;
       int col_flag;
-      index_x_ref = atom->find_custom("BPM_X_REF", type_flag, col_flag);
+      index_x_ref = atom->find_custom(x_ref_id, type_flag, col_flag);
-      index_y_ref = atom->find_custom("BPM_Y_REF", type_flag, col_flag);
+      index_y_ref = atom->find_custom(y_ref_id, type_flag, col_flag);
-      index_z_ref = atom->find_custom("BPM_Z_REF", type_flag, col_flag);
+      index_z_ref = atom->find_custom(z_ref_id, type_flag, col_flag);
+
+      delete [] x_ref_id;
+      delete [] y_ref_id;
+      delete [] z_ref_id;
 
       if (modify->fix[ifix]->restart_reset) {
           modify->fix[ifix]->restart_reset = 0;

src/BPM/bond_bpm.h

@@ -16,6 +16,8 @@
 
 #include "bond.h"
 
+#include <vector>
+
 namespace LAMMPS_NS {
 
 class BondBPM : public Bond {
@@ -36,6 +38,8 @@ class BondBPM : public Bond {
   double r0_max_estimate;
   double max_stretch;
 
+  std::vector<int> leftover_args;
+
   char *id_fix_dummy, *id_fix_update;
   char *id_fix_store_local, *id_fix_prop_atom;
   class FixStoreLocal *fix_store_local;

src/BPM/bond_bpm_rotational.cpp

@@ -176,78 +176,29 @@ void BondBPMRotational::store_data()
   fix_bond_history->post_neighbor();
 }
 
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+  Calculate forces using formulation in:
+    1) Y. Wang Acta Geotechnica 2009
+    2) P. Mora & Y. Wang Advances in Geomcomputing 2009
+---------------------------------------------------------------------- */
 
-void BondBPMRotational::compute(int eflag, int vflag)
+double BondBPMRotational::elastic_forces(int i1, int i2, int type, double& Fr,
+     double r_mag, double r0_mag, double r_mag_inv, double* rhat, double* r,
+     double* r0, double* force1on2, double* torque1on2, double* torque2on1)
 {
-
+  int m;
-  if (! fix_bond_history->stored_flag) {
+  double breaking, temp, r0_dot_rb, c, gamma;
-    fix_bond_history->stored_flag = true;
+  double psi, theta, cos_phi, sin_phi;
-    store_data();
-  }
-
-  int i1,i2,itmp,m,n,type,itype,jtype;
-  double q1[4], q2[4], r[3], r0[3];
-  double rsq, r0_mag, r_mag, r_mag_inv, Fr, Fs_mag;
-  double Tt_mag, Tb_mag, breaking, smooth;
-  double force1on2[3], torque1on2[3], torque2on1[3];
-  double rhat[3], wn1[3], wn2[3], wxn1[3], wxn2[3], vroll[3];
-  double w1dotr, w2dotr, v1dotr, v2dotr;
-  double vn1[3], vn2[3], vt1[3], vt2[3], tmp[3], s1[3], s2[3], tdamp[3];
-  double tor1, tor2, tor3, fs1, fs2, fs3;
-
-  double q2inv[4], rb[3], rb_x_r0[3], s[3], t[3], Fs[3];
-  double q21[4], qp21[4], Tbp[3], Ttp[3];
-  double Tsp[3], Fsp[3], Tt[3], Tb[3], Ts[3], F_rot[3], T_rot[3];
-  double mq[4], mqinv[4], Ttmp[3], Ftmp[3], qtmp[4];
-  double r0_dot_rb, gamma, c, psi, theta, sin_phi, cos_phi, temp;
   double mag_in_plane, mag_out_plane;
+  double  Fs_mag, Tt_mag, Tb_mag;
 
-  ev_init(eflag,vflag);
+  double q1[4], q2[4];
+  double q2inv[4], mq[4], mqinv[4], qp21[4], q21[4], qtmp[4];
+  double rb[3], rb_x_r0[3], s[3], t[3];
+  double Fs[3], Fsp[3], F_rot[3], Ftmp[3];
+  double Ts[3], Tb[3], Tt[3], Tbp[3], Ttp[3], Tsp[3], T_rot[3], Ttmp[3];
 
-  double **x = atom->x;
-  double **v = atom->v;
-  double **omega = atom->omega;
-  double **f = atom->f;
-  double **torque = atom->torque;
-  double *radius = atom->radius;
   double **quat = atom->quat;
-  tagint *tag = atom->tag;
-  int **bondlist = neighbor->bondlist;
-  int nbondlist = neighbor->nbondlist;
-  int nlocal = atom->nlocal;
-  int newton_bond = force->newton_bond;
-
-  double **bondstore = fix_bond_history->bondstore;
-
-  for (n = 0; n < nbondlist; n++) {
-
-    // skip bond if already broken
-    if (bondlist[n][2] <= 0) continue;
-
-    i1 = bondlist[n][0];
-    i2 = bondlist[n][1];
-    type = bondlist[n][2];
-    r0_mag = bondstore[n][0];
-
-    // Ensure pair is always ordered such that r0 points in
-    // a consistent direction and to ensure numerical operations
-    // are identical to minimize the possibility that a bond straddling
-    // an mpi grid (newton off) doesn't break on one proc but not the other
-    if (tag[i2] < tag[i1]) {
-      itmp = i1;
-      i1 = i2;
-      i2 = itmp;
-    }
-
-    // If bond hasn't been set - should be initialized to zero
-    if (r0_mag < EPSILON || std::isnan(r0_mag))
-      r0_mag = store_bond(n,i1,i2);
-
-    r0[0] = bondstore[n][1];
-    r0[1] = bondstore[n][2];
-    r0[2] = bondstore[n][3];
-    MathExtra::scale3(r0_mag, r0);
 
   q1[0] = quat[i1][0];
   q1[1] = quat[i1][1];
@@ -259,20 +210,6 @@ void BondBPMRotational::compute(int eflag, int vflag)
   q2[2] = quat[i2][2];
   q2[3] = quat[i2][3];
 
-    // Note this is the reverse of Mora & Wang
-    MathExtra::sub3(x[i1], x[i2], r);
-
-    rsq = MathExtra::lensq3(r);
-    r_mag = sqrt(rsq);
-    r_mag_inv = 1.0/r_mag;
-    MathExtra::scale3(r_mag_inv, r, rhat);
-
-    // ------------------------------------------------------//
-    //  Calculate forces using formulation in:
-    //    1) Y. Wang Acta Geotechnica 2009
-    //    2) P. Mora & Y. Wang Advances in Geomcomputing 2009
-    // ------------------------------------------------------//
-
   // Calculate normal forces, rb = bond vector in particle 1's frame
   MathExtra::qconjugate(q2, q2inv);
   MathExtra::quatrotvec(q2inv, r, rb);
@@ -416,27 +353,29 @@ void BondBPMRotational::compute(int eflag, int vflag)
   Tt_mag = MathExtra::len3(Tt);
   Tb_mag = MathExtra::len3(Tb);
 
-    // ------------------------------------------------------//
-    //  Check if bond breaks
-    // ------------------------------------------------------//
-
   breaking = Fr/Fcr[type] + Fs_mag/Fcs[type] + Tb_mag/Tcb[type] + Tt_mag/Tct[type];
   if (breaking < 0.0) breaking = 0.0;
 
-    if (breaking >= 1.0) {
+  return breaking;
-      bondlist[n][2] = 0;
+}
-      process_broken(i1, i2);
-      continue;
-    }
 
-    smooth = breaking*breaking;
+/* ----------------------------------------------------------------------
-    smooth = 1.0 - smooth*smooth;
+  Calculate damping using formulation in
+        Y. Wang, F. Alonso-Marroquin, W. Guo 2015
+  Note: n points towards 1 vs pointing towards 2
+---------------------------------------------------------------------- */
 
-    // ------------------------------------------------------//
+void BondBPMRotational::damping_forces(int i1, int i2, int type, double& Fr,
-    //  Calculate damping using formulation in
+     double* rhat, double* r, double* force1on2, double* torque1on2,
-    //    Y. Wang, F. Alonso-Marroquin, W. Guo 2015
+     double* torque2on1)
-    // ------------------------------------------------------//
+{
-    // Note: n points towards 1 vs pointing towards 2
+  double v1dotr, v2dotr, w1dotr, w2dotr;
+  double s1[3], s2[3], tdamp[3], tmp[3];
+  double vn1[3], vn2[3], vt1[3], vt2[3], vroll[3];
+  double wxn1[3], wxn2[3], wn1[3], wn2[3];
+
+  double **v = atom->v;
+  double **omega = atom->omega;
 
   // Damp normal velocity difference
   v1dotr = MathExtra::dot3(v[i1],rhat);
@@ -447,6 +386,7 @@ void BondBPMRotational::compute(int eflag, int vflag)
 
   MathExtra::sub3(vn1, vn2, tmp);
   MathExtra::scale3(gnorm[type], tmp);
+  Fr = MathExtra::lensq3(tmp);
   MathExtra::add3(force1on2, tmp, force1on2);
 
   // Damp tangential objective velocities
@@ -498,6 +438,92 @@ void BondBPMRotational::compute(int eflag, int vflag)
   MathExtra::add3(torque1on2, tdamp, torque1on2);
   MathExtra::scale3(-1.0, tdamp);
   MathExtra::add3(torque2on1, tdamp, torque2on1);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondBPMRotational::compute(int eflag, int vflag)
+{
+
+  if (! fix_bond_history->stored_flag) {
+    fix_bond_history->stored_flag = true;
+    store_data();
+  }
+
+  int i1,i2,itmp,n,type;
+  double r[3], r0[3], rhat[3];
+  double delx, dely, delz, rsq, r0_mag, r_mag, r_mag_inv;
+  double Fr, breaking, smooth;
+  double force1on2[3], torque1on2[3], torque2on1[3];
+
+  ev_init(eflag,vflag);
+
+  double **x = atom->x;
+  double **f = atom->f;
+  double **torque = atom->torque;
+  double *radius = atom->radius;
+  tagint *tag = atom->tag;
+  int **bondlist = neighbor->bondlist;
+  int nbondlist = neighbor->nbondlist;
+  int nlocal = atom->nlocal;
+  int newton_bond = force->newton_bond;
+
+  double **bondstore = fix_bond_history->bondstore;
+
+  for (n = 0; n < nbondlist; n++) {
+
+    // skip bond if already broken
+    if (bondlist[n][2] <= 0) continue;
+
+    i1 = bondlist[n][0];
+    i2 = bondlist[n][1];
+    type = bondlist[n][2];
+    r0_mag = bondstore[n][0];
+
+    // Ensure pair is always ordered such that r0 points in
+    // a consistent direction and to ensure numerical operations
+    // are identical to minimize the possibility that a bond straddling
+    // an mpi grid (newton off) doesn't break on one proc but not the other
+    if (tag[i2] < tag[i1]) {
+      itmp = i1;
+      i1 = i2;
+      i2 = itmp;
+    }
+
+    // If bond hasn't been set - should be initialized to zero
+    if (r0_mag < EPSILON || std::isnan(r0_mag))
+      r0_mag = store_bond(n,i1,i2);
+
+    r0[0] = bondstore[n][1];
+    r0[1] = bondstore[n][2];
+    r0[2] = bondstore[n][3];
+    MathExtra::scale3(r0_mag, r0);
+
+    // Note this is the reverse of Mora & Wang
+    MathExtra::sub3(x[i1], x[i2], r);
+
+    rsq = MathExtra::lensq3(r);
+    r_mag = sqrt(rsq);
+    r_mag_inv = 1.0/r_mag;
+    MathExtra::scale3(r_mag_inv, r, rhat);
+
+    // ------------------------------------------------------//
+    //  Calculate forces, check if bond breaks
+    // ------------------------------------------------------//
+
+    breaking = elastic_forces(i1, i2, type, Fr, r_mag, r0_mag, r_mag_inv,
+            rhat, r, r0, force1on2, torque1on2, torque2on1);
+
+    if (breaking >= 1.0) {
+      bondlist[n][2] = 0;
+      process_broken(i1, i2);
+      continue;
+    }
+
+    damping_forces(i1, i2, type, Fr, rhat, r, force1on2, torque1on2, torque2on1);
+
+    smooth = breaking*breaking;
+    smooth = 1.0 - smooth*smooth;
 
     // ------------------------------------------------------//
     //  Apply forces and torques to particles
@@ -617,7 +643,18 @@ void BondBPMRotational::init_style()
 
   if (!fix_bond_history)
     fix_bond_history = (FixBondHistory *) modify->add_fix(
-         "BOND_HISTORY_BPM_ROTATIONAL all BOND_HISTORY 0 4");
+         "HISTORY_BPM_ROTATIONAL" + std::to_string(instance_me) + " all BOND_HISTORY 0 4");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondBPMRotational::settings(int narg, char **arg)
+{
+  BondBPM::settings(narg, arg);
+
+  for (int iarg : leftover_args) {
+    error->all(FLERR, "Illegal bond_style command");
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -698,10 +735,42 @@ double BondBPMRotational::single(int type, double rsq, int i, int j,
   // Not yet enabled
   if (type <= 0) return 0.0;
 
-  //double r0;
+  int itmp;
-  //for (int n = 0; n < atom->num_bond[i]; n ++) {
+  if (atom->tag[j] < atom->tag[i]) {
-  //  if (atom->bond_atom[i][n] == atom->tag[j]) {
+    itmp = i;
-  //    r0 = fix_bond_history->get_atom_value(i, n, 0);
+    i = j;
-  //  }
+    j = itmp;
-  //}
+  }
+
+  double r0_mag, r_mag, r_mag_inv;
+  double r0[3], r[3], rhat[3];
+  for (int n = 0; n < atom->num_bond[i]; n ++) {
+    if (atom->bond_atom[i][n] == atom->tag[j]) {
+      r0_mag = fix_bond_history->get_atom_value(i, n, 0);
+      r0[0] = fix_bond_history->get_atom_value(i, n, 1);
+      r0[1] = fix_bond_history->get_atom_value(i, n, 2);
+      r0[2] = fix_bond_history->get_atom_value(i, n, 3);
+    }
+  }
+
+  double **x = atom->x;
+  MathExtra::scale3(r0_mag, r0);
+  MathExtra::sub3(x[i], x[j], r);
+
+  r_mag = sqrt(rsq);
+  r_mag_inv = 1.0/r_mag;
+  MathExtra::scale3(r_mag_inv, r, rhat);
+
+  double breaking, smooth, Fr;
+  double force1on2[3], torque1on2[3], torque2on1[3];
+  breaking = elastic_forces(i, j, type, Fr, r_mag, r0_mag, r_mag_inv,
+          rhat, r, r0, force1on2, torque1on2, torque2on1);
+  fforce = Fr;
+  damping_forces(i, j, type, Fr, rhat, r, force1on2, torque1on2, torque2on1);
+  fforce += Fr;
+
+  smooth = breaking*breaking;
+  smooth = 1.0 - smooth*smooth;
+  fforce *= smooth;
+  return 0.0;
 }

src/BPM/bond_bpm_rotational.h

@@ -31,6 +31,7 @@ class BondBPMRotational : public BondBPM {
   virtual void compute(int, int);
   void coeff(int, char **);
   void init_style();
+  void settings(int, char **);
   void write_restart(FILE *);
   void read_restart(FILE *);
   void write_data(FILE *);
@@ -41,6 +42,11 @@ class BondBPMRotational : public BondBPM {
   double *Fcr, *Fcs, *Tct, *Tcb;
   double acos_limit(double);
 
+  double elastic_forces(int, int, int, double &, double, double, double,
+     double*, double*, double*, double*, double*, double*);
+  void damping_forces(int, int, int, double &, double*, double*, double*,
+     double*, double*);
+
   void allocate();
   void store_data();
   double store_bond(int, int, int);

src/BPM/bond_bpm_spring.cpp

@@ -136,7 +136,7 @@ void BondBPMSpring::compute(int eflag, int vflag)
 
   int i1,i2,itmp,m,n,type,itype,jtype;
   double delx, dely, delz, delvx, delvy, delvz;
-  double e, rsq, r, r0, rinv, smooth, fbond, ebond, dot;
+  double e, rsq, r, r0, rinv, smooth, fbond, dot;
 
   ev_init(eflag,vflag);
 
@@ -190,7 +190,6 @@ void BondBPMSpring::compute(int eflag, int vflag)
 
     rinv = 1.0/r;
     fbond = k[type]*(r0-r);
-    if (eflag) ebond = -0.5*fbond*(r0-r);
 
     delvx = v[i1][0] - v[i2][0];
     delvy = v[i1][1] - v[i2][1];
@@ -218,7 +217,7 @@ void BondBPMSpring::compute(int eflag, int vflag)
       f[i2][2] -= delz*fbond;
     }
 
-    if (evflag) ev_tally(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz);
+    if (evflag) ev_tally(i1,i2,nlocal,newton_bond,0.0,fbond,delx,dely,delz);
   }
 }
 
@@ -280,7 +279,18 @@ void BondBPMSpring::init_style()
 
   if (!fix_bond_history)
     fix_bond_history = (FixBondHistory *) modify->add_fix(
-         "BOND_HISTORY_BPM_SPRING all BOND_HISTORY 0 1");
+         "HISTORY_BPM_SPRING_" + std::to_string(instance_me) + " all BOND_HISTORY 0 1");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondBPMSpring::settings(int narg, char **arg)
+{
+  BondBPM::settings(narg, arg);
+
+  for (int iarg : leftover_args) {
+    error->all(FLERR, "Illegal bond_style command");
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -329,13 +339,37 @@ void BondBPMSpring::write_data(FILE *fp)
 double BondBPMSpring::single(int type, double rsq, int i, int j,
                            double &fforce)
 {
-  // Not yet enabled
   if (type <= 0) return 0.0;
 
-  //double r0;
+  double r0;
-  //for (int n = 0; n < atom->num_bond[i]; n ++) {
+  for (int n = 0; n < atom->num_bond[i]; n ++) {
-  //  if (atom->bond_atom[i][n] == atom->tag[j]) {
+    if (atom->bond_atom[i][n] == atom->tag[j]) {
-  //    r0 = fix_bond_history->get_atom_value(i, n, 0);
+      r0 = fix_bond_history->get_atom_value(i, n, 0);
-  //  }
+    }
-  //}
+  }
+
+  double r = sqrt(rsq);
+  double rinv = 1.0/r;
+  double e = (r - r0)/r0;
+  fforce = k[type]*(r0-r);
+
+  double **x = atom->x;
+  double **v = atom->v;
+  double delx = x[i][0] - x[j][0];
+  double dely = x[i][1] - x[j][1];
+  double delz = x[i][2] - x[j][2];
+  double delvx = v[i][0] - v[j][0];
+  double delvy = v[i][1] - v[j][1];
+  double delvz = v[i][2] - v[j][2];
+  double dot = delx*delvx + dely*delvy + delz*delvz;
+  fforce -= gamma[type]*dot*rinv;
+
+  double smooth = (r-r0)/(r0*ecrit[type]);
+  smooth *= smooth;
+  smooth *= smooth;
+  smooth *= smooth;
+  smooth = 1 - smooth;
+
+  fforce *= rinv*smooth;
+  return 0.0;
 }

src/BPM/bond_bpm_spring.h

@@ -31,6 +31,7 @@ class BondBPMSpring : public BondBPM {
   virtual void compute(int, int);
   void coeff(int, char **);
   void init_style();
+  void settings(int, char **);
   void write_restart(FILE *);
   void read_restart(FILE *);
   void write_data(FILE *);

src/BPM/pair_bpm_spring.cpp

@@ -112,7 +112,6 @@ void PairBPMSpring::compute(int eflag, int vflag)
 
         rinv = 1.0/r;
         fpair = k[itype][jtype]*(cut[itype][jtype]-r);
-        if (eflag) evdwl = -0.5 * fpair * (cut[itype][jtype]-r) * factor_lj;
 
         smooth = rsq/cutsq[itype][jtype];
         smooth *= smooth;
@@ -125,6 +124,7 @@ void PairBPMSpring::compute(int eflag, int vflag)
         fpair -= gamma[itype][jtype]*dot*smooth*rinv;
 
         fpair *= factor_lj*rinv;
+        if (eflag) evdwl = 0.0;
 
         f[i][0] += delx*fpair;
         f[i][1] += dely*fpair;
@@ -307,7 +307,7 @@ double PairBPMSpring::single(int i, int j, int itype, int jtype, double rsq,
                             double factor_coul, double factor_lj,
                             double &fforce)
 {
-  double fpair,philj,r,rinv;
+  double fpair,r,rinv;
   double delx, dely, delz, delvx, delvy, delvz, dot, smooth;
 
   if(rsq > cutsq[itype][jtype]) return 0.0;
@@ -319,7 +319,6 @@ double PairBPMSpring::single(int i, int j, int itype, int jtype, double rsq,
   rinv = 1.0/r;
 
   fpair = k[itype][jtype]*(cut[itype][jtype]-r);
-  philj = -0.5*k[itype][jtype]*(cut[itype][jtype]-r)*(cut[itype][jtype]-r);
 
   smooth = rsq/cutsq[itype][jtype];
   smooth *= smooth;
@@ -336,5 +335,5 @@ double PairBPMSpring::single(int i, int j, int itype, int jtype, double rsq,
   fpair *= factor_lj;
   fforce = fpair;
 
-  return philj;
+  return 0.0;
 }

src/bond.cpp

@@ -28,6 +28,10 @@ using namespace LAMMPS_NS;
 
 enum{NONE,LINEAR,SPLINE};
 
+// allocate space for static class instance variable and initialize it
+
+int Bond::instance_total = 0;
+
 /* -----------------------------------------------------------------------
    set bond contribution to Vdwl energy to 0.0
    a particular bond style can override this
@@ -35,6 +39,8 @@ enum{NONE,LINEAR,SPLINE};
 
 Bond::Bond(LAMMPS *lmp) : Pointers(lmp)
 {
+  instance_me = instance_total++;
+
   energy = 0.0;
   virial[0] = virial[1] = virial[2] = virial[3] = virial[4] = virial[5] = 0.0;
   writedata = 1;

src/bond.h

@@ -23,6 +23,8 @@ class Bond : protected Pointers {
   friend class FixOMP;
 
  public:
+  static int instance_total; // # of Bond classes ever instantiated
+
   int allocated;
   int *setflag;
   int partial_flag;          // 1 if bond type can be set to 0 and deleted
@@ -68,6 +70,7 @@ class Bond : protected Pointers {
   void write_file(int, char **);
 
  protected:
+  int instance_me;    // which Bond class instantiation I am
   int suffix_flag;    // suffix compatibility flag
 
   int evflag;

src/fix_property_atom.cpp

@@ -38,6 +38,7 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
 
   restart_peratom = 1;
   wd_section = 1;
+  create_attribute = 1;
 
   int iarg = 3;
   nvalue = narg-iarg;
@@ -554,25 +555,45 @@ void FixPropertyAtom::copy_arrays(int i, int j, int /*delflag*/)
       atom->q[j] = atom->q[i];
     else if (styles[nv] == RMASS)
       atom->rmass[j] = atom->rmass[i];
-    else if (styles[nv] == IVEC) {
+    else if (styles[nv] == IVEC)
       atom->ivector[index[nv]][j] = atom->ivector[index[nv]][i];
-      atom->ivector[index[nv]][i] = 0;
+    else if (styles[nv] == DVEC)
-    } else if (styles[nv] == DVEC) {
       atom->dvector[index[nv]][j] = atom->dvector[index[nv]][i];
-      atom->dvector[index[nv]][i] = 0.0;
+    else if (styles[nv] == IARRAY) {
-    } else if (styles[nv] == IARRAY) {
       ncol = cols[nv];
-      for (k = 0; k < ncol; k++) {
+      for (k = 0; k < ncol; k++)
         atom->iarray[index[nv]][j][k] = atom->iarray[index[nv]][i][k];
-        atom->iarray[index[nv]][i][k] = 0;
-      }
     } else if (styles[nv] == DARRAY) {
       ncol = cols[nv];
-      for (k = 0; k < ncol; k++) {
+      for (k = 0; k < ncol; k++)
         atom->darray[index[nv]][j][k] = atom->darray[index[nv]][i][k];
-        atom->darray[index[nv]][i][k] = 0.0;
     }
   }
+}
+
+
+/* ----------------------------------------------------------------------
+   initialize one atom's storage values, called when atom is created
+------------------------------------------------------------------------- */
+
+void FixPropertyAtom::set_arrays(int i)
+{
+  int k,ncol;
+
+  for (int nv = 0; nv < nvalue; nv++) {
+    if (styles[nv] == IVEC)
+      atom->ivector[index[nv]][i] = 0;
+    else if (styles[nv] == DVEC)
+      atom->dvector[index[nv]][i] = 0.0;
+    else if (styles[nv] == IARRAY) {
+      ncol = cols[nv];
+      for (k = 0; k < ncol; k++)
+        atom->iarray[index[nv]][i][k] = 0;
+    } else if (styles[nv] == DARRAY) {
+      ncol = cols[nv];
+      for (k = 0; k < ncol; k++)
+        atom->darray[index[nv]][i][k] = 0.0;
+    }
   }
 }
 
@@ -706,24 +727,16 @@ int FixPropertyAtom::pack_exchange(int i, double *buf)
     if (styles[nv] == MOLECULE) buf[m++] = ubuf(atom->molecule[i]).d;
     else if (styles[nv] == CHARGE) buf[m++] = atom->q[i];
     else if (styles[nv] == RMASS) buf[m++] = atom->rmass[i];
-    else if (styles[nv] == IVEC) {
+    else if (styles[nv] == IVEC) buf[m++] = ubuf(atom->ivector[index[nv]][i]).d;
-      buf[m++] = ubuf(atom->ivector[index[nv]][i]).d;
+    else if (styles[nv] == DVEC) buf[m++] = atom->dvector[index[nv]][i];
-      atom->ivector[index[nv]][i] = 0;
+    else if (styles[nv] == IARRAY) {
-    } else if (styles[nv] == DVEC) {
-      buf[m++] = atom->dvector[index[nv]][i];
-      atom->dvector[index[nv]][i] = 0.0;
-    } else if (styles[nv] == IARRAY) {
       ncol = cols[nv];
-      for (k = 0; k < ncol; k++) {
+      for (k = 0; k < ncol; k++)
         buf[m++] = ubuf(atom->iarray[index[nv]][i][k]).d;
-        atom->iarray[index[nv]][i][k] = 0;
-      }
     } else if (styles[nv] == DARRAY) {
       ncol = cols[nv];
-      for (k = 0; k < ncol; k++) {
+      for (k = 0; k < ncol; k++)
         buf[m++] = atom->darray[index[nv]][i][k];
-        atom->darray[index[nv]][i][k] = 0.0;
-      }
     }
   }
 

src/fix_property_atom.h

@@ -40,6 +40,7 @@ class FixPropertyAtom : public Fix {
 
   virtual void grow_arrays(int);
   void copy_arrays(int, int, int);
+  void set_arrays(int);
   int pack_border(int, int *, double *);
   int unpack_border(int, int, double *);
   int pack_exchange(int, double *);

src/fix_store_local.h

@@ -13,7 +13,7 @@
 
 #ifdef FIX_CLASS
 // clang-format off
-FixStyle(store/local,FixStoreLocal);
+FixStyle(STORE_LOCAL,FixStoreLocal);
 // clang-format on
 #else
 

src/fix_update_special_bonds.cpp

@@ -24,7 +24,7 @@
 #include "pair.h"
 
 #include <cstring>
-#include <set>
+#include <vector>
 #include <utility>
 
 using namespace LAMMPS_NS;
@@ -36,7 +36,6 @@ FixUpdateSpecialBonds::FixUpdateSpecialBonds(LAMMPS *lmp, int narg, char **arg)
   Fix(lmp, narg, arg)
 {
   if (narg != 3) error->all(FLERR,"Illegal fix update/special/bonds command");
-  comm_forward = 1+atom->maxspecial;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -76,6 +75,7 @@ void FixUpdateSpecialBonds::setup(int /*vflag*/)
       force->special_coul[3] != 1.0)
     error->all(FLERR,"Fix update/special/bonds requires special Coulomb weights = 1,1,1");
 
+  new_broken_pairs.clear();
   broken_pairs.clear();
 }
 
@@ -86,7 +86,7 @@ void FixUpdateSpecialBonds::setup(int /*vflag*/)
 void FixUpdateSpecialBonds::pre_exchange()
 {
   int i, j, key, m, n1, n3;
-  tagint min_tag, max_tag;
+  tagint tagi, tagj;
   int nlocal = atom->nlocal;
 
   tagint *tag = atom->tag;
@@ -94,19 +94,19 @@ void FixUpdateSpecialBonds::pre_exchange()
   int **nspecial = atom->nspecial;
   tagint **special = atom->special;
 
-  for (auto const &key : broken_pairs) {
+  for (auto const &it : broken_pairs) {
-    min_tag = key.first;
+    tagi = it.first;
-    max_tag = key.second;
+    tagj = it.second;
 
-    i = atom->map(min_tag);
+    i = atom->map(tagi);
-    j = atom->map(max_tag);
+    j = atom->map(tagj);
 
     // remove i from special bond list for atom j and vice versa
     if (i < nlocal) {
       slist = special[i];
       n1 = nspecial[i][0];
       for (m = 0; m < n1; m++)
-        if (slist[m] == max_tag) break;
+        if (slist[m] == tagj) break;
       n3 = nspecial[i][2];
       for (; m < n3-1; m++) slist[m] = slist[m+1];
       nspecial[i][0]--;
@@ -118,7 +118,7 @@ void FixUpdateSpecialBonds::pre_exchange()
       slist = special[j];
       n1 = nspecial[j][0];
       for (m = 0; m < n1; m++)
-        if (slist[m] == min_tag) break;
+        if (slist[m] == tagi) break;
       n3 = nspecial[j][2];
       for (; m < n3-1; m++) slist[m] = slist[m+1];
       nspecial[j][0]--;
@@ -127,9 +127,6 @@ void FixUpdateSpecialBonds::pre_exchange()
     }
   }
 
-  // Forward updated special bond list
-  comm->forward_comm_fix(this);
-
   broken_pairs.clear();
 }
 
@@ -139,82 +136,51 @@ void FixUpdateSpecialBonds::pre_exchange()
 
 void FixUpdateSpecialBonds::pre_force(int /*vflag*/)
 {
-  int i,j,n,m,ii,jj,inum,jnum;
+  int i1,i2,j,jj,jnum;
-  int *ilist,*jlist,*numneigh,**firstneigh;
+  int *jlist,*numneigh,**firstneigh;
-  tagint min_tag, max_tag;
+  tagint tag1, tag2;
-  std::pair <tagint, tagint> key;
 
-  int **bond_type = atom->bond_type;
-  int *num_bond = atom->num_bond;
-  tagint **bond_atom = atom->bond_atom;
   int nlocal = atom->nlocal;
 
   tagint *tag = atom->tag;
-  NeighList *list = force->pair->list;
+  NeighList *list = force->pair->list; // may need to be generalized to work with pair hybrid*
-  inum = list->inum;
-  ilist = list->ilist;
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
-  for (ii = 0; ii < inum; ii++) {
+  // In theory could communicate a list of broken bonds to neighboring processors here
-    i = ilist[ii];
+  // to remove restriction that users use Newton bond off
-    jlist = firstneigh[i];
-    jnum = numneigh[i];
 
+  for (auto const &it : new_broken_pairs) {
+    tag1 = it.first;
+    tag2 = it.second;
+    i1 = atom->map(tag1);
+    i2 = atom->map(tag2);
+    
+    // Loop through atoms of owned atoms i j
+    if (i1 < nlocal) {
+      jlist = firstneigh[i1];
+      jnum = numneigh[i1];    
       for (jj = 0; jj < jnum; jj++) {
         j = jlist[jj];
-      j &= NEIGHMASK;
+        j &= SPECIALMASK; // Clear special bond bits
-
+        if (tag[j] == tag2)
-      min_tag = tag[i];
+          jlist[jj] = j;
-      max_tag = tag[j];
-      if (max_tag < min_tag) {
-        min_tag = tag[j];
-        max_tag = tag[i];
-      }
-      key = std::make_pair(min_tag, max_tag);
-
-      if (broken_pairs.find(key) != broken_pairs.end())
-        jlist[jj] = j; // Clear special bond bits
       }
     }
-}
 
-/* ---------------------------------------------------------------------- */
+    if (i2 < nlocal) {
-
+      jlist = firstneigh[i2];
-int FixUpdateSpecialBonds::pack_forward_comm(int n, int *list, double *buf,
+      jnum = numneigh[i2];
-                                    int /*pbc_flag*/, int * /*pbc*/)
+      for (jj = 0; jj < jnum; jj++) {
-{
+        j = jlist[jj];
-  int i,j,k,m,ns;
+        j &= SPECIALMASK; // Clear special bond bits
-  int **nspecial = atom->nspecial;
+        if (tag[j] == tag1)
-  tagint **special = atom->special;
+          jlist[jj] = j;
-
-  m = 0;
-  for (i = 0; i < n; i++) {
-    j = list[i];
-    ns = nspecial[j][0];
-    buf[m++] = ubuf(ns).d;
-    for (k = 0; k < ns; k++)
-      buf[m++] = ubuf(special[j][k]).d;
       }
-  return m;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixUpdateSpecialBonds::unpack_forward_comm(int n, int first, double *buf)
-{
-  int i,j,m,ns,last;
-  int **nspecial = atom->nspecial;
-  tagint **special = atom->special;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
-    ns = (int) ubuf(buf[m++]).i;
-    nspecial[i][0] = ns;
-    for (j = 0; j < ns; j++)
-      special[i][j] = (tagint) ubuf(buf[m++]).i;
     }
+  }
+
+  new_broken_pairs.clear();
 }
 
 /* ---------------------------------------------------------------------- */
@@ -222,14 +188,8 @@ void FixUpdateSpecialBonds::unpack_forward_comm(int n, int first, double *buf)
 void FixUpdateSpecialBonds::add_broken_bond(int i, int j)
 {
   tagint *tag = atom->tag;
+  std::pair <tagint, tagint> tag_pair = std::make_pair(tag[i],tag[j]);
 
-  tagint min_tag = tag[i];
+  new_broken_pairs.push_back(tag_pair);
-  tagint max_tag = tag[j];
+  broken_pairs.push_back(tag_pair);
-  if (max_tag < min_tag) {
-    min_tag = tag[j];
-    max_tag = tag[i];
-  }
-  std::pair <tagint, tagint> key = std::make_pair(min_tag, max_tag);
-
-  broken_pairs.insert(key);
 }

src/fix_update_special_bonds.h

@@ -22,7 +22,7 @@ FixStyle(UPDATE_SPECIAL_BONDS,FixUpdateSpecialBonds)
 
 #include "fix.h"
 
-#include <set>
+#include <vector>
 #include <utility>
 
 namespace LAMMPS_NS {
@@ -35,15 +35,13 @@ class FixUpdateSpecialBonds : public Fix {
   void setup(int);
   void pre_exchange();
   void pre_force(int); 
-  int pack_forward_comm(int, int *, double *, int, int *);
-  void unpack_forward_comm(int, int, double *);  
   void add_broken_bond(int,int);
 
  protected:
-  std::set <std::pair<tagint, tagint>> broken_pairs;
+  // Create two arrays to store bonds broken this timestep (new)
-  inline int sbmask(int j) const {
+  // and since the last neighbor list build
-    return j >> SBBITS & 3;
+  std::vector <std::pair<tagint, tagint>> new_broken_pairs;
-  }
+  std::vector <std::pair<tagint, tagint>> broken_pairs;
 };
 
 }    // namespace LAMMPS_NS

src/lmptype.h

@@ -55,14 +55,15 @@
 
 namespace LAMMPS_NS {
 
-// reserve 2 hi bits in molecular system neigh list for special bonds flag
+// reserve 2 highest bits in molecular system neigh list for special bonds flag
-// reserve 3rd last bit in neigh list for fix neigh/history flag
+// reserve 3rd highest bit in neigh list for fix neigh/history flag
 // max local + ghost atoms per processor = 2^29 - 1
 
 #define SBBITS 30
 #define HISTBITS 29
 #define NEIGHMASK 0x1FFFFFFF
 #define HISTMASK 0xDFFFFFFF
+#define SPECIALMASK 0x3FFFFFFF
 
 // default to 32-bit smallint and other ints, 64-bit bigint