Extra D2min options/checks for undercoord particles

doc/src/fix_nonaffine_displacement.rst

@@ -8,7 +8,7 @@ Syntax
 
 .. parsed-literal::
 
-   fix ID group nonaffine/displacement style args reference/style nstep
+   fix ID group nonaffine/displacement style args reference/style nstep keyword values
 
 * ID, group are documented in :doc:`fix <fix>` command
 * nonaffine/displacement = style name of this fix command
@@ -32,6 +32,13 @@ Syntax
        *update* = update the reference frame every *nstep* timesteps
        *offset* = update the reference frame *nstep* timesteps before calculating the nonaffine displacement
 
+* zero or more keyword/value pairs may be appended
+
+  .. parsed-literal::
+
+       *z/min* values = zmin
+         zmin = minimum coordination number to calculate D2min
+
 Examples
 """"""""
 
@@ -76,6 +83,12 @@ is the identity tensor. This calculation is only performed on timesteps that
 are a multiple of *nevery* (including timestep zero). Data accessed before
 this occurs will simply be zeroed.
 
+For particles with low coordination numbers, calculations of :math:`D^2_\mathrm{min}`
+may have poor accuracy. An optional minimum coordination number can be defined using
+the *z/min* keyword. If any particles have fewer than the specified number of particles
+in the cutoff distance or in contact, the above calculations will be skipped and the
+peratom array entries will be zero.
+
 The *integrated* style simply integrates the velocity of particles
 every timestep to calculate a displacement. This style only works if
 used in conjunction with another fix that deforms the box and displaces

src/EXTRA-FIX/fix_nonaffine_displacement.cpp

@@ -46,6 +46,8 @@ enum { TYPE, RADIUS, CUSTOM };
 enum { INTEGRATED, D2MIN };
 enum { FIXED, OFFSET, UPDATE };
 
+static constexpr double EPSILON = 1.0e-15;
+
 static const char cite_nonaffine_d2min[] =
   "@article{PhysRevE.57.7192,\n"
   " title = {Dynamics of viscoplastic deformation in amorphous solids},\n"
@@ -66,7 +68,7 @@ static const char cite_nonaffine_d2min[] =
 
 FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg), id_fix(nullptr), fix(nullptr), D2min(nullptr), X(nullptr), Y(nullptr),
-  F(nullptr), norm(nullptr)
+  F(nullptr), norm(nullptr), singular(nullptr)
 {
   if (narg < 4) utils::missing_cmd_args(FLERR,"fix nonaffine/displacement", error);
 
@@ -74,6 +76,8 @@ FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char *
   if (nevery <= 0) error->all(FLERR,"Illegal nevery value {} in fix nonaffine/displacement", nevery);
 
   reference_timestep = update_timestep = offset_timestep = -1;
+  z_min = 0;
+
   int iarg = 4;
   if (strcmp(arg[iarg], "integrated") == 0) {
     nad_style = INTEGRATED;
@@ -117,6 +121,16 @@ FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char *
     if ((offset_timestep <= 0) || (offset_timestep > nevery))
       error->all(FLERR, "Illegal offset timestep {} in fix nonaffine/displacement", arg[iarg + 1]);
   } else error->all(FLERR,"Illegal reference style {} in fix nonaffine/displacement", arg[iarg]);
+  iarg += 2;
+
+  while (iarg < narg) {
+    if (strcmp(arg[iarg], "z/min") == 0) {
+      if (iarg + 2 > narg) utils::missing_cmd_args(FLERR,"fix nonaffine/displacement", error);
+      z_min = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
+      if (z_min < 0) error->all(FLERR, "Minimum coordination must be positive");
+      iarg += 2;
+    } else error->all(FLERR,"Illegal keyword {} in fix nonaffine/displacement", arg[iarg]);
+  }
 
   if (nad_style == D2MIN)
     if (cut_style == RADIUS && (!atom->radius_flag))
@@ -151,6 +165,7 @@ FixNonaffineDisplacement::~FixNonaffineDisplacement()
     memory->destroy(Y);
     memory->destroy(F);
     memory->destroy(norm);
+    memory->destroy(singular);
     memory->destroy(D2min);
   }
 
@@ -395,6 +410,7 @@ void FixNonaffineDisplacement::calculate_D2Min()
       }
     }
     norm[i] = 0;
+    singular[i] = 0;
     D2min[i] = 0;
   }
 
@@ -471,15 +487,30 @@ void FixNonaffineDisplacement::calculate_D2Min()
     }
 
     if (dim == 3) {
+      denom = det3(Y_tmp);
+      if (fabs(denom) < EPSILON) {
+        singular[i] = 1;
+        for (j = 0; j < 3; j++)
+          for (k = 0; k < 3; k++)
+            Y_inv[j][k] = 0.0;
+      } else {
         invert3(Y_tmp, Y_inv);
+      }
     } else {
       denom = Y_tmp[0][0] * Y_tmp[1][1] - Y_tmp[0][1] * Y_tmp[1][0];
-      if (denom != 0.0) denom = 1.0 / denom;
+      if (fabs(denom) < EPSILON) {
+        singular[i] = 1;
+        for (j = 0; j < 2; j++)
+          for (k = 0; k < 2; k++)
+            Y_inv[j][k] = 0.0;
+      } else {
+        denom = 1.0 / denom;
         Y_inv[0][0] = Y_tmp[1][1] * denom;
         Y_inv[0][1] = -Y_tmp[0][1] * denom;
         Y_inv[1][0] = -Y_tmp[1][0] * denom;
         Y_inv[1][1] = Y_tmp[0][0] * denom;
       }
+    }
 
     times3(X_tmp, Y_inv, F_tmp);
 
@@ -559,10 +590,14 @@ void FixNonaffineDisplacement::calculate_D2Min()
   for (i = 0; i < nlocal; i++) {
     if (!(mask[i] & groupbit)) continue;
 
-    if (norm[i] != 0)
+    if (norm[i] < z_min || singular[i] == 1) {
+      array_atom[i][0] = 0.0;
+      array_atom[i][1] = 0.0;
+      array_atom[i][2] = 0.0;
+      continue;
+    }
+
     D2min[i] /= norm[i];
-    else
-      D2min[i] = 0.0;
 
     for (j = 0; j < 3; j++)
       for (k = 0; k < 3; k++)
@@ -743,10 +778,12 @@ void FixNonaffineDisplacement::grow_arrays(int nmax_new)
     memory->destroy(F);
     memory->destroy(D2min);
     memory->destroy(norm);
+    memory->destroy(singular);
     memory->create(X, nmax, 3, 3, "fix_nonaffine_displacement:X");
     memory->create(Y, nmax, 3, 3, "fix_nonaffine_displacement:Y");
     memory->create(F, nmax, 3, 3, "fix_nonaffine_displacement:F");
     memory->create(D2min, nmax, "fix_nonaffine_displacement:D2min");
     memory->create(norm, nmax, "fix_nonaffine_displacement:norm");
+    memory->create(singular, nmax, "fix_nonaffine_displacement:singular");
   }
 }

src/EXTRA-FIX/fix_nonaffine_displacement.h

@@ -48,12 +48,12 @@ class FixNonaffineDisplacement : public Fix {
   int nmax, comm_flag;
   int nad_style, cut_style;
   int reference_style, offset_timestep, reference_timestep, update_timestep;
-  int reference_saved;
+  int reference_saved, z_min;
   double cutoff_custom, cutsq_custom, mycutneigh;
   double xprd0, yprd0, zprd0, xprd0_half, yprd0_half, zprd0_half, xy0, xz0, yz0;
 
   double *D2min, ***X, ***Y, ***F;
-  int *norm;
+  int *norm, *singular;
 
   class NeighList *list;    // half neighbor list