Misc updates and style fixes

doc/src/bond_bpm_spring.rst

@@ -64,7 +64,7 @@ go to zero, avoiding discontinuities, as bonds approach the critical strain
 
 .. math::
 
-   w = 1.0 - \left( \frac{r - r_0}{r_0 \epsilon_c} \right^4 .
+   w = 1.0 - \left( \frac{r - r_0}{r_0 \epsilon_c} \right)^8 .
 
 Finally, an additional damping force is applied to the bonded particles.
 This forces is proportional to the difference in the

doc/src/pair_bpm_spring.rst

@@ -38,8 +38,15 @@ normal velocity of particles
    F_D = - \gamma w (\hat{r} \bullet \vec{v})
 
 where :math:`\gamma` is the damping strength, :math:`\hat{r}` is the
-radial normal vector, and :math:`\vec{v}` is the velocity difference
-between the two particles.
+radial normal vector, :math:`\vec{v}` is the velocity difference
+between the two particles, and :math:`w` is a smoothing factor.
+This smoothing factor is constructed such that damping forces go to zero
+as particles come out of contact to avoid discontinuities. It is
+given by
+
+.. math::
+
+   w = 1.0 - \left( \frac{r}{r_c} \right)^8 .
 
 This pair style is designed for use in a spring-based bonded particle model.
 It mirrors the construction of the :doc:`bpm/spring <bond_bpm_spring>` bond style.

src/BPM/bond_bpm.cpp

@@ -73,9 +73,12 @@ void BondBPM::init_style()
   ifix = modify->find_fix_by_style("update/special/bonds");
   if (ifix >= 0)
     fix_update_special_bonds = (FixUpdateSpecialBonds *) modify->fix[ifix];
-  else
+  else {
     fix_update_special_bonds = nullptr;
-
+    if (force->special_lj[1] != 1.0)
+      error->all(FLERR, "Without fix update/special/bonds, BPM bond styles "
+                        "require special_bonds weight of 1.0 for first neighbors");
+  }
 
   if (force->angle || force->dihedral || force->improper)
     error->all(FLERR,

src/BPM/bond_bpm_rotational.cpp

@@ -25,11 +25,8 @@
 #include "memory.h"
 #include "modify.h"
 #include "neighbor.h"
-#include "pair.h"
 
-#include <mpi.h>
 #include <cmath>
-#include <cstring>
 
 #define EPSILON 1e-10
 
@@ -147,8 +144,7 @@ void BondBPMRotational::store_data()
       type = bond_type[i][m];
 
       //Skip if bond was turned off
-      if(type < 0)
-          continue;
+      if(type < 0) continue;
 
       // map to find index n for tag
       j = atom->map(atom->bond_atom[i][m]);
@@ -167,9 +163,9 @@ void BondBPMRotational::store_data()
 
       // Get closest image in case bonded with ghost
       domain->minimum_image(delx, dely, delz);
-
       r = sqrt(delx*delx + dely*dely + delz*delz);
       rinv = 1.0/r;
+
       fix_bond_history->update_atom_value(i, m, 0, r);
       fix_bond_history->update_atom_value(i, m, 1, delx*rinv);
       fix_bond_history->update_atom_value(i, m, 2, dely*rinv);
@@ -209,7 +205,6 @@ void BondBPMRotational::compute(int eflag, int vflag)
 
   ev_init(eflag,vflag);
 
-  double **cutsq = force->pair->cutsq;
   double **x = atom->x;
   double **v = atom->v;
   double **omega = atom->omega;

src/BPM/bond_bpm_spring.cpp

@@ -20,21 +20,13 @@
 #include "error.h"
 #include "fix_bond_history.h"
 #include "force.h"
-#include "math_const.h"
-#include "math_extra.h"
 #include "memory.h"
 #include "modify.h"
 #include "neighbor.h"
-#include "pair.h"
-
-#include <mpi.h>
-#include <cmath>
-#include <cstring>
 
 #define EPSILON 1e-10
 
 using namespace LAMMPS_NS;
-using namespace MathExtra;
 
 /* ---------------------------------------------------------------------- */
 
@@ -111,8 +103,7 @@ void BondBPMSpring::store_data()
       type = bond_type[i][m];
 
       //Skip if bond was turned off
-      if (type < 0)
-          continue;
+      if (type < 0) continue;
 
       // map to find index n
       j = atom->map(atom->bond_atom[i][m]);
@@ -125,6 +116,7 @@ void BondBPMSpring::store_data()
       // Get closest image in case bonded with ghost
       domain->minimum_image(delx, dely, delz);
       r = sqrt(delx*delx + dely*dely + delz*delz);
+
       fix_bond_history->update_atom_value(i, m, 0, r);
     }
   }
@@ -142,17 +134,15 @@ void BondBPMSpring::compute(int eflag, int vflag)
     store_data();
   }
 
-  int i1,i2,itmp,m,n,type,itype,jtype;
+  int i1,i2,m,n,type,itype,jtype;
   double delx, dely, delz, delvx, delvy, delvz;
   double e, rsq, r, r0, rinv, smooth, fbond, ebond, dot;
 
   ev_init(eflag,vflag);
 
-  double **cutsq = force->pair->cutsq;
   double **x = atom->x;
   double **v = atom->v;
   double **f = atom->f;
-  tagint *tag = atom->tag;
   int **bondlist = neighbor->bondlist;
   int nbondlist = neighbor->nbondlist;
   int nlocal = atom->nlocal;
@@ -170,6 +160,10 @@ void BondBPMSpring::compute(int eflag, int vflag)
     type = bondlist[n][2];
     r0 = bondstore[n][0];
 
+    // If bond hasn't been set - should be initialized to zero
+    if (r0 < EPSILON || std::isnan(r0))
+      r0 = store_bond(n,i1,i2);
+
     delx = x[i1][0] - x[i2][0];
     dely = x[i1][1] - x[i2][1];
     delz = x[i1][2] - x[i2][2];
@@ -197,6 +191,7 @@ void BondBPMSpring::compute(int eflag, int vflag)
     smooth = (r-r0)/(r0*ecrit[type]);
     smooth *= smooth;
     smooth *= smooth;
+    smooth *= smooth;
     smooth = 1 - smooth;
     
     fbond *= rinv*smooth;

src/BPM/fix_bond_history.cpp

@@ -171,7 +171,7 @@ void FixBondHistory::pre_exchange()
       for (m = 0; m < num_bond[i2]; m++) {
         if (bond_atom[i2][m] == tag[i1]) {
           for (idata = 0; idata < ndata; idata++) {
-            stored[i1][m*ndata+idata] = bondstore[n][idata];
+            stored[i2][m*ndata+idata] = bondstore[n][idata];
           }
         }
       }

src/BPM/pair_bpm_spring.cpp

@@ -14,17 +14,13 @@
 
 #include "pair_bpm_spring.h"
 
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
 #include "atom.h"
-#include "atom_vec.h"
 #include "comm.h"
+#include "error.h"
 #include "force.h"
+#include "memory.h"
 #include "neighbor.h"
 #include "neigh_list.h"
-#include "memory.h"
-#include "error.h"
 
 using namespace LAMMPS_NS;
 
@@ -116,17 +112,17 @@ void PairBPMSpring::compute(int eflag, int vflag)
 
         rinv = 1.0/r;
         fpair = k[itype][jtype]*(cut[itype][jtype]-r);
-        if (eflag)
-          evdwl = -0.5*k[itype][jtype]*(cut[itype][jtype]-r)*(cut[itype][jtype]-r);
+        if (eflag) evdwl = -0.5 * fpair * (cut[itype][jtype]-r) * factor_lj;
 
         smooth = rsq/cutsq[itype][jtype];
         smooth *= smooth;
+        smooth *= smooth;
         smooth = 1.0 - smooth;
         delvx = vxtmp - v[j][0];
         delvy = vytmp - v[j][1];
         delvz = vztmp - v[j][2];
         dot = delx*delvx + dely*delvy + delz*delvz;
-        fpair -= gamma[itype][jtype]*dot*smooth;
+        fpair -= gamma[itype][jtype]*dot*smooth*rinv;
 
         fpair *= factor_lj*rinv;
 
@@ -190,8 +186,8 @@ void PairBPMSpring::coeff(int narg, char **arg)
   if (!allocated) allocate();
 
   int ilo,ihi,jlo,jhi;
-  utils::bounds(FLERR,arg[0],1,atom->nbondtypes,ilo,ihi,error);
-  utils::bounds(FLERR,arg[1],1,atom->nbondtypes,jlo,jhi,error);
+  utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
+  utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
 
   double k_one = utils::numeric(FLERR,arg[2],false,lmp);
   double cut_one = utils::numeric(FLERR,arg[3],false,lmp);

src/BPM/pair_bpm_spring.h

@@ -12,9 +12,9 @@
 ------------------------------------------------------------------------- */
 
 #ifdef PAIR_CLASS
-
+// clang-format off
 PairStyle(bpm/spring,PairBPMSpring)
-
+// clang-format on
 #else
 
 #ifndef LMP_PAIR_BPM_SPRING_H