diff --git a/doc/src/bond_bpm_spring.rst b/doc/src/bond_bpm_spring.rst index 1a67c09a70..e149a45d2e 100644 --- a/doc/src/bond_bpm_spring.rst +++ b/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 diff --git a/doc/src/pair_bpm_spring.rst b/doc/src/pair_bpm_spring.rst index 19bb987ae3..33883f05ef 100755 --- a/doc/src/pair_bpm_spring.rst +++ b/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 style. diff --git a/src/BPM/bond_bpm.cpp b/src/BPM/bond_bpm.cpp index a6e5d4583c..eacfd064e5 100644 --- a/src/BPM/bond_bpm.cpp +++ b/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, diff --git a/src/BPM/bond_bpm_rotational.cpp b/src/BPM/bond_bpm_rotational.cpp index aa7de10c13..0f62d88966 100644 --- a/src/BPM/bond_bpm_rotational.cpp +++ b/src/BPM/bond_bpm_rotational.cpp @@ -25,11 +25,8 @@ #include "memory.h" #include "modify.h" #include "neighbor.h" -#include "pair.h" -#include #include -#include #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; diff --git a/src/BPM/bond_bpm_spring.cpp b/src/BPM/bond_bpm_spring.cpp index 028c972500..b418a19583 100644 --- a/src/BPM/bond_bpm_spring.cpp +++ b/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 -#include -#include #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; diff --git a/src/BPM/fix_bond_history.cpp b/src/BPM/fix_bond_history.cpp index f60661ef2f..5184882961 100644 --- a/src/BPM/fix_bond_history.cpp +++ b/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]; } } } diff --git a/src/BPM/pair_bpm_spring.cpp b/src/BPM/pair_bpm_spring.cpp index 29ea1c0f66..abaf868313 100644 --- a/src/BPM/pair_bpm_spring.cpp +++ b/src/BPM/pair_bpm_spring.cpp @@ -14,17 +14,13 @@ #include "pair_bpm_spring.h" -#include -#include -#include #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); diff --git a/src/BPM/pair_bpm_spring.h b/src/BPM/pair_bpm_spring.h index 86be2a9867..216fd0a9dc 100644 --- a/src/BPM/pair_bpm_spring.h +++ b/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