From 3dbfe26b6df145d369dda65b02a9007f5fe5b74e Mon Sep 17 00:00:00 2001 From: jtclemm Date: Tue, 23 Apr 2024 20:58:20 -0600 Subject: [PATCH] Extra D2min options/checks for undercoord particles --- doc/src/fix_nonaffine_displacement.rst | 15 ++++- src/EXTRA-FIX/fix_nonaffine_displacement.cpp | 59 ++++++++++++++++---- src/EXTRA-FIX/fix_nonaffine_displacement.h | 4 +- 3 files changed, 64 insertions(+), 14 deletions(-) diff --git a/doc/src/fix_nonaffine_displacement.rst b/doc/src/fix_nonaffine_displacement.rst index 0a271ebc32..becc47eec8 100644 --- a/doc/src/fix_nonaffine_displacement.rst +++ b/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 ` 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 diff --git a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp index eaf45f4e59..846f434dce 100644 --- a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp +++ b/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,14 +487,29 @@ void FixNonaffineDisplacement::calculate_D2Min() } if (dim == 3) { - invert3(Y_tmp, Y_inv); + 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; - 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; + 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) - D2min[i] /= norm[i]; - else - D2min[i] = 0.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]; 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"); } } diff --git a/src/EXTRA-FIX/fix_nonaffine_displacement.h b/src/EXTRA-FIX/fix_nonaffine_displacement.h index c7177bd3d9..b0e9c464ca 100644 --- a/src/EXTRA-FIX/fix_nonaffine_displacement.h +++ b/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