From c57ed87e9a9f73dfbc4e6e62a5f98ac46e083f48 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 23 Nov 2023 15:36:27 +0200 Subject: [PATCH 01/18] Addition of conical indenter --- src/fix_indent.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/fix_indent.h b/src/fix_indent.h index 527e9ec277..6f33f6fbb1 100644 --- a/src/fix_indent.h +++ b/src/fix_indent.h @@ -49,6 +49,10 @@ class FixIndent : public Fix { int cdim, varflag; int ilevel_respa; + char *rlostr, *rhistr, *lostr, *histr; + int rlovar, rhivar, lovar, hivar; + double rlovalue, rhivalue, lovalue, hivalue; + void options(int, char **); }; From a90c7b42f99bf61d79a7cdacc69394b56de829b0 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 23 Nov 2023 15:40:34 +0200 Subject: [PATCH 02/18] Include code for conical indenter in fix_indent.cpp --- src/fix_indent.cpp | 466 ++++++++++++++++++++++++++++++++++++++------- 1 file changed, 398 insertions(+), 68 deletions(-) diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp index 9adb337dd6..d3e80ecd79 100644 --- a/src/fix_indent.cpp +++ b/src/fix_indent.cpp @@ -23,6 +23,7 @@ #include "error.h" #include "input.h" #include "lattice.h" +#include "math_extra.h" #include "modify.h" #include "respa.h" #include "update.h" @@ -34,14 +35,30 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum{NONE,SPHERE,CYLINDER,PLANE}; -enum{INSIDE,OUTSIDE}; +enum{NONE, SPHERE, CYLINDER, PLANE, CONE}; +enum{INSIDE, OUTSIDE}; + +static bool PointInsideCone(int dir, double *center, double lo, + double hi, double rlo, double rhi, double *point); + +static void DistanceExteriorPoint(int dir, double *center, + double lo, double hi, double rlo, double rhi, double &x, + double &y, double &z); + +static void DistanceInteriorPoint(int dir, double *center, + double lo, double hi, double rlo, double rhi, double &x, + double &y, double &z); + +static void point_on_line_segment(double *a, double *b, double *c, double *d); + +static double closest(double *x, double *near, double *nearest, double dsq); /* ---------------------------------------------------------------------- */ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), - xstr(nullptr), ystr(nullptr), zstr(nullptr), rstr(nullptr), pstr(nullptr) + xstr(nullptr), ystr(nullptr), zstr(nullptr), rstr(nullptr), pstr(nullptr), + rlostr(nullptr), rhistr(nullptr), lostr(nullptr), histr(nullptr) { if (narg < 4) utils::missing_cmd_args(FLERR, "fix indent", error); @@ -56,6 +73,7 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : ilevel_respa = 0; k = utils::numeric(FLERR,arg[3],false,lmp); + if (k < 0.0) error->all(FLERR, "Illegal fix indent force constant: {}", k); k3 = k/3.0; // read options from end of input line @@ -79,6 +97,30 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : if (!ystr) yvalue *= yscale; if (!zstr) zvalue *= zscale; if (!rstr) rvalue *= xscale; + } else if (istyle == CONE) { + + if (!xstr) xvalue *= xscale; + if (!ystr) yvalue *= yscale; + if (!zstr) zvalue *= zscale; + + double scaling_factor = 1.0; + switch (cdim) { + case 0: + scaling_factor = xscale; + break; + case 1: + scaling_factor = yscale; + break; + case 2: + scaling_factor = zscale; + break; + } + + if (!rlostr) rlovalue *= scaling_factor; + if (!rhistr) rhivalue *= scaling_factor; + if (!lostr) lovalue *= scaling_factor; + if (!histr) hivalue *= scaling_factor; + } else if (istyle == PLANE) { if (cdim == 0 && !pstr) pvalue *= xscale; else if (cdim == 1 && !pstr) pvalue *= yscale; @@ -86,7 +128,7 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : } else error->all(FLERR,"Unknown fix indent keyword: {}", istyle); varflag = 0; - if (xstr || ystr || zstr || rstr || pstr) varflag = 1; + if (xstr || ystr || zstr || rstr || pstr || rlostr || rhistr || lostr || histr) varflag = 1; indenter_flag = 0; indenter[0] = indenter[1] = indenter[2] = indenter[3] = 0.0; @@ -101,6 +143,10 @@ FixIndent::~FixIndent() delete [] zstr; delete [] rstr; delete [] pstr; + delete [] rlostr; + delete [] rhistr; + delete [] lostr; + delete [] histr; } /* ---------------------------------------------------------------------- */ @@ -154,6 +200,38 @@ void FixIndent::init() error->all(FLERR,"Variable {} for fix indent is invalid style", pstr); } + if (rlostr) { + rlovar = input->variable->find(rlostr); + if (rlovar < 0) + error->all(FLERR,"Variable {} for fix indent does not exist", rlostr); + if (!input->variable->equalstyle(rlovar)) + error->all(FLERR,"Variable {} for fix indent is invalid style", rlostr); + } + + if (rhistr) { + rhivar = input->variable->find(rhistr); + if (rhivar < 0) + error->all(FLERR,"Variable {} for fix indent does not exist", rhistr); + if (!input->variable->equalstyle(rhivar)) + error->all(FLERR,"Variable {} for fix indent is invalid style", rhistr); + } + + if (lostr) { + lovar = input->variable->find(lostr); + if (lovar < 0) + error->all(FLERR,"Variable {} for fix indent does not exist", lostr); + if (!input->variable->equalstyle(lovar)) + error->all(FLERR,"Variable {} for fix indent is invalid style", lostr); + } + + if (histr) { + hivar = input->variable->find(histr); + if (hivar < 0) + error->all(FLERR,"Variable {} for fix indent does not exist", histr); + if (!input->variable->equalstyle(hivar)) + error->all(FLERR,"Variable {} for fix indent is invalid style", histr); + } + if (utils::strmatch(update->integrate_style,"^respa")) { ilevel_respa = (dynamic_cast(update->integrate))->nlevels-1; if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa); @@ -192,32 +270,28 @@ void FixIndent::post_force(int /*vflag*/) indenter_flag = 0; indenter[0] = indenter[1] = indenter[2] = indenter[3] = 0.0; + // ctr = current indenter centerz + double ctr[3] {xvalue, yvalue, zvalue}; + if (xstr) ctr[0] = input->variable->compute_equal(xvar); + if (ystr) ctr[1] = input->variable->compute_equal(yvar); + if (zstr) ctr[2] = input->variable->compute_equal(zvar); + + double **x = atom->x; + double **f = atom->f; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double delx, dely, delz, r, dr, fmag, fx, fy, fz; + // spherical indenter if (istyle == SPHERE) { - // ctr = current indenter center - // remap into periodic box - - double ctr[3]; - if (xstr) ctr[0] = input->variable->compute_equal(xvar); - else ctr[0] = xvalue; - if (ystr) ctr[1] = input->variable->compute_equal(yvar); - else ctr[1] = yvalue; - if (zstr) ctr[2] = input->variable->compute_equal(zvar); - else ctr[2] = zvalue; + // remap indenter center into periodic box domain->remap(ctr); - double radius; - if (rstr) radius = input->variable->compute_equal(rvar); - else radius = rvalue; - - double **x = atom->x; - double **f = atom->f; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - double delx,dely,delz,r,dr,fmag,fx,fy,fz; + double radius { rstr ? input->variable->compute_equal(rvar) : rvalue}; + if (radius < 0.0) error->all(FLERR, "Illegal fix indent sphere radius: {}", radius); for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { @@ -254,38 +328,11 @@ void FixIndent::post_force(int /*vflag*/) // remap into periodic box // 3rd coord is just near box for remap(), since isn't used - double ctr[3]; - if (cdim == 0) { - ctr[0] = domain->boxlo[0]; - if (ystr) ctr[1] = input->variable->compute_equal(yvar); - else ctr[1] = yvalue; - if (zstr) ctr[2] = input->variable->compute_equal(zvar); - else ctr[2] = zvalue; - } else if (cdim == 1) { - if (xstr) ctr[0] = input->variable->compute_equal(xvar); - else ctr[0] = xvalue; - ctr[1] = domain->boxlo[1]; - if (zstr) ctr[2] = input->variable->compute_equal(zvar); - else ctr[2] = zvalue; - } else { - if (xstr) ctr[0] = input->variable->compute_equal(xvar); - else ctr[0] = xvalue; - if (ystr) ctr[1] = input->variable->compute_equal(yvar); - else ctr[1] = yvalue; - ctr[2] = domain->boxlo[2]; - } + ctr[cdim] = domain->boxlo[cdim]; domain->remap(ctr); - double radius; - if (rstr) radius = input->variable->compute_equal(rvar); - else radius = rvalue; - - double **x = atom->x; - double **f = atom->f; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - double delx,dely,delz,r,dr,fmag,fx,fy,fz; + double radius { rstr ? input->variable->compute_equal(rvar) : rvalue}; + if (radius < 0.0) error->all(FLERR, "Illegal fix indent cylinder radius: {}", radius); for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { @@ -324,31 +371,85 @@ void FixIndent::post_force(int /*vflag*/) indenter[3] -= fz; } + // conical indenter + + } else if (istyle == CONE) { + + double radiuslo { rlostr ? input->variable->compute_equal(rlovar) : rlovalue }; + if (radiuslo < 0.0) error->all(FLERR, "Illegal fix indent cone lower radius: {}", radiuslo); + + double radiushi { rhistr ? input->variable->compute_equal(rhivar) : rhivalue }; + if (radiushi < 0.0) error->all(FLERR, "Illegal fix indent cone high radius: {}", radiushi); + + double initial_lo { lostr ? input->variable->compute_equal(lovar) : lovalue }; + double initial_hi { histr ? input->variable->compute_equal(hivar) : hivalue }; + + ctr[cdim] = 0.5 * (initial_hi + initial_lo); + + domain->remap(ctr); + + double hi = ctr[cdim] + 0.5 * (initial_hi - initial_lo); + double lo = ctr[cdim] - 0.5 * (initial_hi - initial_lo); + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + + delx = x[i][0] - ctr[0]; + dely = x[i][1] - ctr[1]; + delz = x[i][2] - ctr[2]; + domain->minimum_image(delx, dely, delz); + + double x0[3] {delx + ctr[0], dely + ctr[1], delz + ctr[2]}; + r = sqrt(delx * delx + dely * dely + delz * delz); + + // find if the particle is inside or outside the cone + bool point_inside_cone = PointInsideCone(cdim, ctr, lo, hi, radiuslo, radiushi, x0); + + if (side == INSIDE && point_inside_cone) continue; + if (side == OUTSIDE && !point_inside_cone) continue; + + // find the distance between the point and the cone + if (point_inside_cone) { + DistanceInteriorPoint(cdim, ctr, lo, hi, radiuslo, radiushi, x0[0], x0[1], x0[2]); + } else { + DistanceExteriorPoint(cdim, ctr, lo, hi, radiuslo, radiushi, x0[0], x0[1], x0[2]); + } + + // compute the force from the center of the cone - it is different from the approach of fix wall/region + dr = sqrt(x0[0] * x0[0] + x0[1] * x0[1] + x0[2] * x0[2]); + + int force_sign = { point_inside_cone ? 1 : -1 }; + fmag = force_sign * k * dr * dr; + + fx = delx*fmag/r; + fy = dely*fmag/r; + fz = delz*fmag/r; + f[i][0] += fx; + f[i][1] += fy; + f[i][2] += fz; + indenter[0] -= k3 * dr * dr * dr; + indenter[1] -= fx; + indenter[2] -= fy; + indenter[3] -= fz; + } + } + // planar indenter } else { // plane = current plane position - double plane; - if (pstr) plane = input->variable->compute_equal(pvar); - else plane = pvalue; - - double **x = atom->x; - double **f = atom->f; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - double dr,fatom; + double plane { pstr ? input->variable->compute_equal(pvar) : pvalue}; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { dr = planeside * (plane - x[i][cdim]); if (dr >= 0.0) continue; - fatom = -planeside * k*dr*dr; - f[i][cdim] += fatom; - indenter[0] -= k3 * dr*dr*dr; - indenter[cdim+1] -= fatom; + fmag = -planeside * k * dr * dr; + f[i][cdim] += fmag; + indenter[0] -= k3 * dr * dr * dr; + indenter[cdim+1] -= fmag; } } @@ -487,6 +588,64 @@ void FixIndent::options(int narg, char **arg) istyle = PLANE; iarg += 4; + } else if (strcmp(arg[iarg],"cone") == 0) { + + if (iarg+8 > narg) utils::missing_cmd_args(FLERR, "fix indent cone", error); + + if (strcmp(arg[iarg+1],"x") == 0) { + cdim = 0; + + if (utils::strmatch(arg[iarg+2],"^v_")) { + ystr = utils::strdup(arg[iarg+2]+2); + } else yvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); + + if (utils::strmatch(arg[iarg+3],"^v_")) { + zstr = utils::strdup(arg[iarg+3]+2); + } else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); + + } else if (strcmp(arg[iarg+1],"y") == 0) { + cdim = 1; + + if (utils::strmatch(arg[iarg+2],"^v_")) { + xstr = utils::strdup(arg[iarg+2]+2); + } else xvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); + + if (utils::strmatch(arg[iarg+3],"^v_")) { + zstr = utils::strdup(arg[iarg+3]+2); + } else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); + + } else if (strcmp(arg[iarg+1],"z") == 0) { + cdim = 2; + + if (utils::strmatch(arg[iarg+2],"^v_")) { + xstr = utils::strdup(arg[iarg+2]+2); + } else xvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); + + if (utils::strmatch(arg[iarg+3],"^v_")) { + ystr = utils::strdup(arg[iarg+3]+2); + } else yvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); + + } else error->all(FLERR,"Unknown fix indent cone argument: {}", arg[iarg+1]); + + if (utils::strmatch(arg[iarg+4],"^v_")) { + rlostr = utils::strdup(arg[iarg+4]+2); + } else rlovalue = utils::numeric(FLERR,arg[iarg+4],false,lmp); + + if (utils::strmatch(arg[iarg+5],"^v_")) { + rhistr = utils::strdup(arg[iarg+5]+2); + } else rhivalue = utils::numeric(FLERR,arg[iarg+5],false,lmp); + + if (utils::strmatch(arg[iarg+6],"^v_")) { + lostr = utils::strdup(arg[iarg+6]+2); + } else lovalue = utils::numeric(FLERR,arg[iarg+6],false,lmp); + + if (utils::strmatch(arg[iarg+7],"^v_")) { + histr = utils::strdup(arg[iarg+7]+2); + } else hivalue = utils::numeric(FLERR,arg[iarg+7],false,lmp); + + istyle = CONE; + iarg += 8; + } else if (strcmp(arg[iarg],"units") == 0) { if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix indent units", error); if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; @@ -503,3 +662,174 @@ void FixIndent::options(int narg, char **arg) } else error->all(FLERR,"Unknown fix indent argument: {}", arg[iarg]); } } + +/* ---------------------------------------------------------------------- + determines if a point is inside (true) or outside (false) of a cone +------------------------------------------------------------------------- */ + +bool PointInsideCone(int dir, double *center, double lo, + double hi, double rlo, double rhi, double *x) +{ + if ((x[dir] > hi) || (x[dir] < lo)) return false; + + double del[3] {x[0] - center[0], x[1] - center[1], x[2] - center[2]}; + del[dir] = 0.0; + + double dist = sqrt(del[0] * del[0] + del[1] * del[1] + del[2] * del[2]); + double currentradius = rlo + (x[dir] - lo) * (rhi - rlo) / (hi - lo); + + if (dist > currentradius) return false; + + return true; +} + +/* ---------------------------------------------------------------------- + distance between an exterior point and a cone +------------------------------------------------------------------------- */ +void DistanceExteriorPoint(int dir, double *center, double lo, double hi, + double rlo, double rhi, double &x, double &y, double &z) +{ + double xp[3], nearest[3], corner1[3], corner2[3]; + + double point[3] {x, y, z}; + + double del[3] {x - center[0], y - center[1], z - center[2]}; + del[dir] = 0.0; + + double r = sqrt(del[0] * del[0] + del[1] * del[1] + del[2] * del[2]); + + corner1[0] = center[0] + del[0] * rlo / r; + corner1[1] = center[1] + del[1] * rlo / r; + corner1[2] = center[2] + del[2] * rlo / r; + corner1[dir] = lo; + + corner2[0] = center[0] + del[0] * rhi / r; + corner2[1] = center[1] + del[1] * rhi / r; + corner2[2] = center[2] + del[2] * rhi / r; + corner2[dir] = hi; + + double corner3[3] {center[0], center[1], center[2]}; + corner3[dir] = lo; + + double corner4[3] {center[0], center[1], center[2]}; + corner4[dir] = hi; + + // initialize distance to a big number + double distsq = 1.0e20; + + // check the first triangle + point_on_line_segment(corner1, corner2, point, xp); + distsq = closest(point, xp, nearest, distsq); + + // check the second triangle + point_on_line_segment(corner1, corner3, point, xp); + distsq = closest(point, xp, nearest, distsq); + + // check the third triangle + point_on_line_segment(corner2, corner4, point, xp); + distsq = closest(point, xp, nearest, distsq); + + x -= nearest[0]; + y -= nearest[1]; + z -= nearest[2]; + + return; +} + +/* ---------------------------------------------------------------------- + distance between an interior point and a cone +------------------------------------------------------------------------- */ + +void DistanceInteriorPoint(int dir, double *center, + double lo, double hi, double rlo, double rhi, double &x, + double &y, double &z) +{ + double r, dist_disk, dist_surf; + double surflo[3], surfhi[3], xs[3]; + double initial_point[3] {x, y, z}; + double point[3] {0.0, 0.0, 0.0}; + + // initial check with the two disks + if ( (initial_point[dir] - lo) < (hi - initial_point[dir]) ) { + dist_disk = (initial_point[dir] - lo) * (initial_point[dir] - lo); + point[dir] = initial_point[dir] - lo; + } else { + dist_disk = (hi - initial_point[dir]) * (hi - initial_point[dir]); + point[dir] = initial_point[dir] - hi; + } + + // check with the points in the conical surface + double del[3] {x - center[0], y - center[1], z - center[2]}; + del[dir] = 0.0; + r = sqrt(del[0] * del[0] + del[1] * del[1] + del[2] * del[2]); + + surflo[0] = center[0] + del[0] * rlo / r; + surflo[1] = center[1] + del[1] * rlo / r; + surflo[2] = center[2] + del[2] * rlo / r; + surflo[dir] = lo; + + surfhi[0] = center[0] + del[0] * rhi / r; + surfhi[1] = center[1] + del[1] * rhi / r; + surfhi[2] = center[2] + del[2] * rhi / r; + surfhi[dir] = hi; + + point_on_line_segment(surflo, surfhi, initial_point, xs); + + double dx[3] {initial_point[0] - xs[0], initial_point[1] - xs[1], initial_point[2] - xs[2]}; + dist_surf = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; + if (dist_surf < dist_disk) { + x = dx[0]; + y = dx[1]; + z = dx[2]; + } else { + x = point[0]; + y = point[1]; + z = point[2]; + } + + return; +} + +/* ---------------------------------------------------------------------- + helper function extracted from region.cpp +------------------------------------------------------------------------- */ + +void point_on_line_segment(double *a, double *b, double *c, double *d) +{ + double ba[3], ca[3]; + + MathExtra::sub3(b, a, ba); + MathExtra::sub3(c, a, ca); + double t = MathExtra::dot3(ca, ba) / MathExtra::dot3(ba, ba); + if (t <= 0.0) { + d[0] = a[0]; + d[1] = a[1]; + d[2] = a[2]; + } else if (t >= 1.0) { + d[0] = b[0]; + d[1] = b[1]; + d[2] = b[2]; + } else { + d[0] = a[0] + t * ba[0]; + d[1] = a[1] + t * ba[1]; + d[2] = a[2] + t * ba[2]; + } +} + +/* ---------------------------------------------------------------------- + helper function extracted from region_cone.cpp +------------------------------------------------------------------------- */ + +double closest(double *x, double *near, double *nearest, double dsq) +{ + double dx = x[0] - near[0]; + double dy = x[1] - near[1]; + double dz = x[2] - near[2]; + double rsq = dx * dx + dy * dy + dz * dz; + if (rsq >= dsq) return dsq; + + nearest[0] = near[0]; + nearest[1] = near[1]; + nearest[2] = near[2]; + return rsq; +} From e1c028c785fface11d8eba7d6ab57e42d278e654 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 23 Nov 2023 16:02:22 +0200 Subject: [PATCH 03/18] Update documentation Describe the arguments for cone indenter style --- doc/src/fix_indent.rst | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/doc/src/fix_indent.rst b/doc/src/fix_indent.rst index 15790e15d0..d478d7dc50 100644 --- a/doc/src/fix_indent.rst +++ b/doc/src/fix_indent.rst @@ -14,17 +14,23 @@ Syntax * indent = style name of this fix command * K = force constant for indenter surface (force/distance\^2 units) * one or more keyword/value pairs may be appended -* keyword = *sphere* or *cylinder* or *plane* or *side* or *units* +* keyword = *sphere* or *cone* or *cylinder* or *plane* or *side* or *units* .. parsed-literal:: *sphere* args = x y z R - x,y,z = position of center of indenter (distance units) + x, y, z = position of center of indenter (distance units) R = sphere radius of indenter (distance units) - any of x,y,z,R can be a variable (see below) + any of x, y, z, R can be a variable (see below) + *cone* args = dim c1 c2 radlo radhi lo hi + dim = *x* or *y* or *z* = axis of cone + c1, c2 = coords of cone axis in other 2 dimensions (distance units) + radlo,radhi = cone radii at lo and hi end (distance units) + lo,hi = bounds of cone in dim (distance units) + any of c1, c2, radlo, radhi, lo, hi can be a variable (see below) *cylinder* args = dim c1 c2 R dim = *x* or *y* or *z* = axis of cylinder - c1,c2 = coords of cylinder axis in other 2 dimensions (distance units) + c1, c2 = coords of cylinder axis in other 2 dimensions (distance units) R = cylinder radius of indenter (distance units) any of c1,c2,R can be a variable (see below) *plane* args = dim pos side @@ -57,7 +63,7 @@ material or as an obstacle in a flow. Or it can be used as a constraining wall around a simulation; see the discussion of the *side* keyword below. -The indenter can either be spherical or cylindrical or planar. You +The indenter can either be spherical or conical or cylindrical or planar. You must set one of those 3 keywords. A spherical indenter exerts a force of magnitude @@ -75,9 +81,9 @@ A cylindrical indenter exerts the same force, except that *r* is the distance from the atom to the center axis of the cylinder. The cylinder extends infinitely along its axis. -Spherical and cylindrical indenters account for periodic boundaries in +Spherical, conical and cylindrical indenters account for periodic boundaries in two ways. First, the center point of a spherical indenter (x,y,z) or -axis of a cylindrical indenter (c1,c2) is remapped back into the +axis of a conical/cylindrical indenter (c1,c2) is remapped back into the simulation box, if the box is periodic in a particular dimension. This occurs every timestep if the indenter geometry is specified with a variable (see below), e.g. it is moving over time. Second, the From ac9afb26dd7ea27b40bbcd2bd44fcdf1ae3932ae Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 4 Jan 2024 17:38:14 +0200 Subject: [PATCH 04/18] aesthetic optimization for xscale, yscale and zscale --- src/fix_indent.cpp | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp index d3e80ecd79..28dcac08d0 100644 --- a/src/fix_indent.cpp +++ b/src/fix_indent.cpp @@ -82,13 +82,9 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : // setup scaling - double xscale,yscale,zscale; - if (scaleflag) { - xscale = domain->lattice->xlattice; - yscale = domain->lattice->ylattice; - zscale = domain->lattice->zlattice; - } - else xscale = yscale = zscale = 1.0; + const double xscale { scaleflag ? domain->lattice->xlattice : 1.0}; + const double yscale { scaleflag ? domain->lattice->ylattice : 1.0}; + const double zscale { scaleflag ? domain->lattice->zlattice : 1.0}; // apply scaling factors to geometry From 921ddffda234de84c9696a17ed895b05d934ec15 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sun, 14 Jan 2024 17:18:33 +0100 Subject: [PATCH 05/18] small code simplification for the cylinder indenter --- src/fix_indent.cpp | 25 +++++++------------------ 1 file changed, 7 insertions(+), 18 deletions(-) diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp index 28dcac08d0..65cf91c2d2 100644 --- a/src/fix_indent.cpp +++ b/src/fix_indent.cpp @@ -332,21 +332,10 @@ void FixIndent::post_force(int /*vflag*/) for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (cdim == 0) { - delx = 0; - dely = x[i][1] - ctr[1]; - delz = x[i][2] - ctr[2]; - } else if (cdim == 1) { - delx = x[i][0] - ctr[0]; - dely = 0; - delz = x[i][2] - ctr[2]; - } else { - delx = x[i][0] - ctr[0]; - dely = x[i][1] - ctr[1]; - delz = 0; - } - domain->minimum_image(delx,dely,delz); - r = sqrt(delx*delx + dely*dely + delz*delz); + double del[3] {x[i][0] - ctr[0], x[i][1] - ctr[1], x[i][2] - ctr[2]}; + del[cdim] = 0; + domain->minimum_image(del[0], del[1], del[2]); + r = sqrt(del[0]*del[0] + del[1]*del[1] + del[2]*del[2]); if (side == OUTSIDE) { dr = r - radius; fmag = k*dr*dr; @@ -355,9 +344,9 @@ void FixIndent::post_force(int /*vflag*/) fmag = -k*dr*dr; } if (dr >= 0.0) continue; - fx = delx*fmag/r; - fy = dely*fmag/r; - fz = delz*fmag/r; + fx = del[0]*fmag/r; + fy = del[1]*fmag/r; + fz = del[2]*fmag/r; f[i][0] += fx; f[i][1] += fy; f[i][2] += fz; From 4e77556610b7955ab10b7ba363dd2ff065976e3d Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 26 Feb 2024 15:59:03 -0700 Subject: [PATCH 06/18] minor changes to source and doc files --- doc/src/fix_indent.rst | 82 +++++---- src/fix_indent.cpp | 371 ++++++++++++++++++++++------------------- src/fix_indent.h | 13 ++ 3 files changed, 261 insertions(+), 205 deletions(-) diff --git a/doc/src/fix_indent.rst b/doc/src/fix_indent.rst index d478d7dc50..5658c06373 100644 --- a/doc/src/fix_indent.rst +++ b/doc/src/fix_indent.rst @@ -8,13 +8,12 @@ Syntax .. code-block:: LAMMPS - fix ID group-ID indent K keyword values ... + fix ID group-ID indent K gstyle args keyword value ... * ID, group-ID are documented in :doc:`fix ` command * indent = style name of this fix command * K = force constant for indenter surface (force/distance\^2 units) -* one or more keyword/value pairs may be appended -* keyword = *sphere* or *cone* or *cylinder* or *plane* or *side* or *units* +* gstyle = *sphere* or *cylinder* or *cone* or *plane* .. parsed-literal:: @@ -22,22 +21,28 @@ Syntax x, y, z = position of center of indenter (distance units) R = sphere radius of indenter (distance units) any of x, y, z, R can be a variable (see below) + *cylinder* args = dim c1 c2 R + dim = *x* or *y* or *z* = axis of cylinder + c1, c2 = coords of cylinder axis in other 2 dimensions (distance units) + R = cylinder radius of indenter (distance units) + any of c1,c2,R can be a variable (see below) *cone* args = dim c1 c2 radlo radhi lo hi dim = *x* or *y* or *z* = axis of cone c1, c2 = coords of cone axis in other 2 dimensions (distance units) radlo,radhi = cone radii at lo and hi end (distance units) lo,hi = bounds of cone in dim (distance units) any of c1, c2, radlo, radhi, lo, hi can be a variable (see below) - *cylinder* args = dim c1 c2 R - dim = *x* or *y* or *z* = axis of cylinder - c1, c2 = coords of cylinder axis in other 2 dimensions (distance units) - R = cylinder radius of indenter (distance units) - any of c1,c2,R can be a variable (see below) *plane* args = dim pos side dim = *x* or *y* or *z* = plane perpendicular to this dimension pos = position of plane in dimension x, y, or z (distance units) pos can be a variable (see below) side = *lo* or *hi* + +* zero or more keyword/value pairs may be appended +* keyword = *side* or *units* + + .. parsed-literal:: + *side* value = *in* or *out* *in* = the indenter acts on particles inside the sphere or cylinder *out* = the indenter acts on particles outside the sphere or cylinder @@ -63,8 +68,8 @@ material or as an obstacle in a flow. Or it can be used as a constraining wall around a simulation; see the discussion of the *side* keyword below. -The indenter can either be spherical or conical or cylindrical or planar. You -must set one of those 3 keywords. +The *gstyle* geometry of the indenter can either be a sphere, a +cylinder, a cone, or a plane. A spherical indenter exerts a force of magnitude @@ -81,15 +86,20 @@ A cylindrical indenter exerts the same force, except that *r* is the distance from the atom to the center axis of the cylinder. The cylinder extends infinitely along its axis. -Spherical, conical and cylindrical indenters account for periodic boundaries in -two ways. First, the center point of a spherical indenter (x,y,z) or -axis of a conical/cylindrical indenter (c1,c2) is remapped back into the -simulation box, if the box is periodic in a particular dimension. -This occurs every timestep if the indenter geometry is specified with -a variable (see below), e.g. it is moving over time. Second, the -calculation of distance to the indenter center or axis accounts for -periodic boundaries. Both of these mean that an indenter can -effectively move through and straddle one or more periodic boundaries. +A conical indenter is similar to a cylindrical indenter except that it +has a finite length (between *lo* and *hi*), and that two different +radii (one at each end, *radlo* and *radhi*) can be defined. + +Spherical, cylindrical, and conical indenters account for periodic +boundaries in two ways. First, the center point of a spherical +indenter (x,y,z) or axis of a cylindrical/conical indenter (c1,c2) is +remapped back into the simulation box, if the box is periodic in a +particular dimension. This occurs every timestep if the indenter +geometry is specified with a variable (see below), e.g. it is moving +over time. Second, the calculation of distance to the indenter center +or axis accounts for periodic boundaries. Both of these mean that an +indenter can effectively move through and straddle one or more +periodic boundaries. A planar indenter is really an axis-aligned infinite-extent wall exerting the same force on atoms in the system, where *R* is the @@ -103,9 +113,13 @@ is specified as *hi*\ . Any of the 4 quantities defining a spherical indenter's geometry can be specified as an equal-style :doc:`variable `, namely *x*, -*y*, *z*, or *R*\ . Similarly, for a cylindrical indenter, any of *c1*, -*c2*, or *R*, can be a variable. For a planar indenter, *pos* can be -a variable. If the value is a variable, it should be specified as +*y*, *z*, or *R*\ . For a cylindrical indenter, any of the 3 +quantities *c1*, *c2*, or *R*, can be a variable. For a conical +indenter, any of the 6 quantities *c1*, *c2*, *radlo*, *radhi*, *lo*, +or *hi* can be a variable. For a planar indenter, the single value +*pos* can be a variable. + +If any of these values is a variable, it should be specified as v_name, where name is the variable name. In this case, the variable will be evaluated each timestep, and its value used to define the indenter geometry. @@ -116,7 +130,8 @@ command keywords for the simulation box parameters and timestep and elapsed time. Thus it is easy to specify indenter properties that change as a function of time or span consecutive runs in a continuous fashion. For the latter, see the *start* and *stop* keywords of the -:doc:`run ` command and the *elaplong* keyword of :doc:`thermo_style custom ` for details. +:doc:`run ` command and the *elaplong* keyword of +:doc:`thermo_style custom ` for details. For example, if a spherical indenter's x-position is specified as v_x, then this variable definition will keep it's center at a relative @@ -147,12 +162,13 @@ rate. If the *side* keyword is specified as *out*, which is the default, then particles outside the indenter are pushed away from its outer -surface, as described above. This only applies to spherical or -cylindrical indenters. If the *side* keyword is specified as *in*, -the action of the indenter is reversed. Particles inside the indenter -are pushed away from its inner surface. In other words, the indenter -is now a containing wall that traps the particles inside it. If the -radius shrinks over time, it will squeeze the particles. +surface, as described above. This only applies to spherical, +cylindrical, and conical indenters. If the *side* keyword is +specified as *in*, the action of the indenter is reversed. Particles +inside the indenter are pushed away from its inner surface. In other +words, the indenter is now a containing wall that traps the particles +inside it. If the radius shrinks over time, it will squeeze the +particles. The *units* keyword determines the meaning of the distance units used to define the indenter geometry. A *box* value selects standard @@ -172,10 +188,10 @@ lattice spacings in a variable formula. The force constant *K* is not affected by the *units* keyword. It is always in force/distance\^2 units where force and distance are defined -by the :doc:`units ` command. If you wish K to be scaled by the -lattice spacing, you can define K with a variable whose formula -contains *xlat*, *ylat*, *zlat* keywords of the -:doc:`thermo_style ` command, e.g. +by the :doc:`units ` command. If you wish K to be scaled by +the lattice spacing, you can define K with a variable whose formula +contains *xlat*, *ylat*, *zlat* keywords of the :doc:`thermo_style +` command, e.g. .. code-block:: LAMMPS diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp index 65cf91c2d2..8d450bee75 100644 --- a/src/fix_indent.cpp +++ b/src/fix_indent.cpp @@ -38,21 +38,6 @@ using namespace FixConst; enum{NONE, SPHERE, CYLINDER, PLANE, CONE}; enum{INSIDE, OUTSIDE}; -static bool PointInsideCone(int dir, double *center, double lo, - double hi, double rlo, double rhi, double *point); - -static void DistanceExteriorPoint(int dir, double *center, - double lo, double hi, double rlo, double rhi, double &x, - double &y, double &z); - -static void DistanceInteriorPoint(int dir, double *center, - double lo, double hi, double rlo, double rhi, double &x, - double &y, double &z); - -static void point_on_line_segment(double *a, double *b, double *c, double *d); - -static double closest(double *x, double *near, double *nearest, double dsq); - /* ---------------------------------------------------------------------- */ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : @@ -76,10 +61,11 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : if (k < 0.0) error->all(FLERR, "Illegal fix indent force constant: {}", k); k3 = k/3.0; - // read options from end of input line - - options(narg-4,&arg[4]); + // read geometry of indenter and optional args + int iarg = geometry(narg-4,&arg[4]) + 4; + options(narg-iarg,&arg[iarg]); + // setup scaling const double xscale { scaleflag ? domain->lattice->xlattice : 1.0}; @@ -93,8 +79,8 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : if (!ystr) yvalue *= yscale; if (!zstr) zvalue *= zscale; if (!rstr) rvalue *= xscale; + } else if (istyle == CONE) { - if (!xstr) xvalue *= xscale; if (!ystr) yvalue *= yscale; if (!zstr) zvalue *= zscale; @@ -121,6 +107,7 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : if (cdim == 0 && !pstr) pvalue *= xscale; else if (cdim == 1 && !pstr) pvalue *= yscale; else if (cdim == 2 && !pstr) pvalue *= zscale; + } else error->all(FLERR,"Unknown fix indent keyword: {}", istyle); varflag = 0; @@ -195,7 +182,6 @@ void FixIndent::init() if (!input->variable->equalstyle(pvar)) error->all(FLERR,"Variable {} for fix indent is invalid style", pstr); } - if (rlostr) { rlovar = input->variable->find(rlostr); if (rlovar < 0) @@ -203,7 +189,6 @@ void FixIndent::init() if (!input->variable->equalstyle(rlovar)) error->all(FLERR,"Variable {} for fix indent is invalid style", rlostr); } - if (rhistr) { rhivar = input->variable->find(rhistr); if (rhivar < 0) @@ -211,7 +196,6 @@ void FixIndent::init() if (!input->variable->equalstyle(rhivar)) error->all(FLERR,"Variable {} for fix indent is invalid style", rhistr); } - if (lostr) { lovar = input->variable->find(lostr); if (lovar < 0) @@ -219,7 +203,6 @@ void FixIndent::init() if (!input->variable->equalstyle(lovar)) error->all(FLERR,"Variable {} for fix indent is invalid style", lostr); } - if (histr) { hivar = input->variable->find(histr); if (hivar < 0) @@ -267,6 +250,7 @@ void FixIndent::post_force(int /*vflag*/) indenter[0] = indenter[1] = indenter[2] = indenter[3] = 0.0; // ctr = current indenter centerz + double ctr[3] {xvalue, yvalue, zvalue}; if (xstr) ctr[0] = input->variable->compute_equal(xvar); if (ystr) ctr[1] = input->variable->compute_equal(yvar); @@ -284,6 +268,7 @@ void FixIndent::post_force(int /*vflag*/) if (istyle == SPHERE) { // remap indenter center into periodic box + domain->remap(ctr); double radius { rstr ? input->variable->compute_equal(rvar) : rvalue}; @@ -387,20 +372,24 @@ void FixIndent::post_force(int /*vflag*/) double x0[3] {delx + ctr[0], dely + ctr[1], delz + ctr[2]}; r = sqrt(delx * delx + dely * dely + delz * delz); - // find if the particle is inside or outside the cone + // check if particle is inside or outside the cone + bool point_inside_cone = PointInsideCone(cdim, ctr, lo, hi, radiuslo, radiushi, x0); if (side == INSIDE && point_inside_cone) continue; if (side == OUTSIDE && !point_inside_cone) continue; // find the distance between the point and the cone + if (point_inside_cone) { DistanceInteriorPoint(cdim, ctr, lo, hi, radiuslo, radiushi, x0[0], x0[1], x0[2]); } else { DistanceExteriorPoint(cdim, ctr, lo, hi, radiuslo, radiushi, x0[0], x0[1], x0[2]); } - // compute the force from the center of the cone - it is different from the approach of fix wall/region + // compute the force from the center of the cone + // this is different from how it is done in fix wall/region + dr = sqrt(x0[0] * x0[0] + x0[1] * x0[1] + x0[2] * x0[2]); int force_sign = { point_inside_cone ? 1 : -1 }; @@ -486,10 +475,10 @@ double FixIndent::compute_vector(int n) } /* ---------------------------------------------------------------------- - parse optional parameters at end of input line + parse input args for geometry of indenter ------------------------------------------------------------------------- */ -void FixIndent::options(int narg, char **arg) +int FixIndent::geometry(int narg, char **arg) { if (narg < 0) utils::missing_cmd_args(FLERR, "fix indent", error); @@ -499,139 +488,168 @@ void FixIndent::options(int narg, char **arg) scaleflag = 1; side = OUTSIDE; + // sphere + + if (strcmp(arg[0],"sphere") == 0) { + if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword"); + if (5 > narg) utils::missing_cmd_args(FLERR, "fix indent sphere", error); + + if (utils::strmatch(arg[1],"^v_")) { + xstr = utils::strdup(arg[1]+2); + } else xvalue = utils::numeric(FLERR,arg[1],false,lmp); + if (utils::strmatch(arg[2],"^v_")) { + ystr = utils::strdup(arg[2]+2); + } else yvalue = utils::numeric(FLERR,arg[2],false,lmp); + if (utils::strmatch(arg[3],"^v_")) { + zstr = utils::strdup(arg[3]+2); + } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); + if (utils::strmatch(arg[4],"^v_")) { + rstr = utils::strdup(arg[4]+2); + } else rvalue = utils::numeric(FLERR,arg[4],false,lmp); + + istyle = SPHERE; + return 5; + } + + // cylinder + + if (strcmp(arg[0],"cylinder") == 0) { + if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword"); + if (5 > narg) utils::missing_cmd_args(FLERR, "fix indent cylinder", error); + + if (strcmp(arg[1],"x") == 0) { + cdim = 0; + if (utils::strmatch(arg[2],"^v_")) { + ystr = utils::strdup(arg[2]+2); + } else yvalue = utils::numeric(FLERR,arg[2],false,lmp); + if (utils::strmatch(arg[3],"^v_")) { + zstr = utils::strdup(arg[3]+2); + } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); + } else if (strcmp(arg[1],"y") == 0) { + cdim = 1; + if (utils::strmatch(arg[2],"^v_")) { + xstr = utils::strdup(arg[2]+2); + } else xvalue = utils::numeric(FLERR,arg[2],false,lmp); + if (utils::strmatch(arg[3],"^v_")) { + zstr = utils::strdup(arg[3]+2); + } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); + } else if (strcmp(arg[1],"z") == 0) { + cdim = 2; + if (utils::strmatch(arg[2],"^v_")) { + xstr = utils::strdup(arg[2]+2); + } else xvalue = utils::numeric(FLERR,arg[2],false,lmp); + if (utils::strmatch(arg[3],"^v_")) { + ystr = utils::strdup(arg[3]+2); + } else yvalue = utils::numeric(FLERR,arg[3],false,lmp); + } else error->all(FLERR,"Unknown fix indent cylinder argument: {}", arg[1]); + + if (utils::strmatch(arg[4],"^v_")) { + rstr = utils::strdup(arg[4]+2); + } else rvalue = utils::numeric(FLERR,arg[4],false,lmp); + + istyle = CYLINDER; + return 5; + } + + // cone + + if (strcmp(arg[0],"cone") == 0) { + if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword"); + if (8 > narg) utils::missing_cmd_args(FLERR, "fix indent cone", error); + + if (strcmp(arg[1],"x") == 0) { + cdim = 0; + + if (utils::strmatch(arg[2],"^v_")) { + ystr = utils::strdup(arg[2]+2); + } else yvalue = utils::numeric(FLERR,arg[2],false,lmp); + + if (utils::strmatch(arg[3],"^v_")) { + zstr = utils::strdup(arg[3]+2); + } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); + + } else if (strcmp(arg[1],"y") == 0) { + cdim = 1; + + if (utils::strmatch(arg[2],"^v_")) { + xstr = utils::strdup(arg[2]+2); + } else xvalue = utils::numeric(FLERR,arg[2],false,lmp); + + if (utils::strmatch(arg[3],"^v_")) { + zstr = utils::strdup(arg[3]+2); + } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); + + } else if (strcmp(arg[1],"z") == 0) { + cdim = 2; + + if (utils::strmatch(arg[2],"^v_")) { + xstr = utils::strdup(arg[2]+2); + } else xvalue = utils::numeric(FLERR,arg[2],false,lmp); + + if (utils::strmatch(arg[3],"^v_")) { + ystr = utils::strdup(arg[3]+2); + } else yvalue = utils::numeric(FLERR,arg[3],false,lmp); + + } else error->all(FLERR,"Unknown fix indent cone argument: {}", arg[1]); + + if (utils::strmatch(arg[4],"^v_")) { + rlostr = utils::strdup(arg[4]+2); + } else rlovalue = utils::numeric(FLERR,arg[4],false,lmp); + + if (utils::strmatch(arg[5],"^v_")) { + rhistr = utils::strdup(arg[5]+2); + } else rhivalue = utils::numeric(FLERR,arg[5],false,lmp); + + if (utils::strmatch(arg[6],"^v_")) { + lostr = utils::strdup(arg[6]+2); + } else lovalue = utils::numeric(FLERR,arg[6],false,lmp); + + if (utils::strmatch(arg[7],"^v_")) { + histr = utils::strdup(arg[7]+2); + } else hivalue = utils::numeric(FLERR,arg[7],false,lmp); + + istyle = CONE; + return 8; + } + + // plane + + if (strcmp(arg[0],"plane") == 0) { + if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword"); + if (4 > narg) utils::missing_cmd_args(FLERR, "fix indent plane", error); + if (strcmp(arg[1],"x") == 0) cdim = 0; + else if (strcmp(arg[1],"y") == 0) cdim = 1; + else if (strcmp(arg[1],"z") == 0) cdim = 2; + else error->all(FLERR,"Unknown fix indent plane argument: {}", arg[1]); + + if (utils::strmatch(arg[2],"^v_")) { + pstr = utils::strdup(arg[2]+2); + } else pvalue = utils::numeric(FLERR,arg[2],false,lmp); + + if (strcmp(arg[3],"lo") == 0) planeside = -1; + else if (strcmp(arg[3],"hi") == 0) planeside = 1; + else error->all(FLERR,"Unknown fix indent plane argument: {}", arg[3]); + istyle = PLANE; + return 4; + } + + // invalid istyle arg + + error->all(FLERR,"Unknown fix indent argument: {}", arg[0]); + + return 0; +} + +/* ---------------------------------------------------------------------- + parse optional input args +------------------------------------------------------------------------- */ + +void FixIndent::options(int narg, char **arg) +{ int iarg = 0; + while (iarg < narg) { - if (strcmp(arg[iarg],"sphere") == 0) { - if (iarg+5 > narg) utils::missing_cmd_args(FLERR, "fix indent sphere", error); - - if (utils::strmatch(arg[iarg+1],"^v_")) { - xstr = utils::strdup(arg[iarg+1]+2); - } else xvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (utils::strmatch(arg[iarg+2],"^v_")) { - ystr = utils::strdup(arg[iarg+2]+2); - } else yvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (utils::strmatch(arg[iarg+3],"^v_")) { - zstr = utils::strdup(arg[iarg+3]+2); - } else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - if (utils::strmatch(arg[iarg+4],"^v_")) { - rstr = utils::strdup(arg[iarg+4]+2); - } else rvalue = utils::numeric(FLERR,arg[iarg+4],false,lmp); - - istyle = SPHERE; - iarg += 5; - - } else if (strcmp(arg[iarg],"cylinder") == 0) { - if (iarg+5 > narg) utils::missing_cmd_args(FLERR, "fix indent cylinder", error); - - if (strcmp(arg[iarg+1],"x") == 0) { - cdim = 0; - if (utils::strmatch(arg[iarg+2],"^v_")) { - ystr = utils::strdup(arg[iarg+2]+2); - } else yvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (utils::strmatch(arg[iarg+3],"^v_")) { - zstr = utils::strdup(arg[iarg+3]+2); - } else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - } else if (strcmp(arg[iarg+1],"y") == 0) { - cdim = 1; - if (utils::strmatch(arg[iarg+2],"^v_")) { - xstr = utils::strdup(arg[iarg+2]+2); - } else xvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (utils::strmatch(arg[iarg+3],"^v_")) { - zstr = utils::strdup(arg[iarg+3]+2); - } else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - } else if (strcmp(arg[iarg+1],"z") == 0) { - cdim = 2; - if (utils::strmatch(arg[iarg+2],"^v_")) { - xstr = utils::strdup(arg[iarg+2]+2); - } else xvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (utils::strmatch(arg[iarg+3],"^v_")) { - ystr = utils::strdup(arg[iarg+3]+2); - } else yvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - } else error->all(FLERR,"Unknown fix indent cylinder argument: {}", arg[iarg+1]); - - if (utils::strmatch(arg[iarg+4],"^v_")) { - rstr = utils::strdup(arg[iarg+4]+2); - } else rvalue = utils::numeric(FLERR,arg[iarg+4],false,lmp); - - istyle = CYLINDER; - iarg += 5; - - } else if (strcmp(arg[iarg],"plane") == 0) { - if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "fix indent plane", error); - if (strcmp(arg[iarg+1],"x") == 0) cdim = 0; - else if (strcmp(arg[iarg+1],"y") == 0) cdim = 1; - else if (strcmp(arg[iarg+1],"z") == 0) cdim = 2; - else error->all(FLERR,"Unknown fix indent plane argument: {}", arg[iarg+1]); - - if (utils::strmatch(arg[iarg+2],"^v_")) { - pstr = utils::strdup(arg[iarg+2]+2); - } else pvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - - if (strcmp(arg[iarg+3],"lo") == 0) planeside = -1; - else if (strcmp(arg[iarg+3],"hi") == 0) planeside = 1; - else error->all(FLERR,"Unknown fix indent plane argument: {}", arg[iarg+3]); - istyle = PLANE; - iarg += 4; - - } else if (strcmp(arg[iarg],"cone") == 0) { - - if (iarg+8 > narg) utils::missing_cmd_args(FLERR, "fix indent cone", error); - - if (strcmp(arg[iarg+1],"x") == 0) { - cdim = 0; - - if (utils::strmatch(arg[iarg+2],"^v_")) { - ystr = utils::strdup(arg[iarg+2]+2); - } else yvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - - if (utils::strmatch(arg[iarg+3],"^v_")) { - zstr = utils::strdup(arg[iarg+3]+2); - } else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - - } else if (strcmp(arg[iarg+1],"y") == 0) { - cdim = 1; - - if (utils::strmatch(arg[iarg+2],"^v_")) { - xstr = utils::strdup(arg[iarg+2]+2); - } else xvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - - if (utils::strmatch(arg[iarg+3],"^v_")) { - zstr = utils::strdup(arg[iarg+3]+2); - } else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - - } else if (strcmp(arg[iarg+1],"z") == 0) { - cdim = 2; - - if (utils::strmatch(arg[iarg+2],"^v_")) { - xstr = utils::strdup(arg[iarg+2]+2); - } else xvalue = utils::numeric(FLERR,arg[iarg+2],false,lmp); - - if (utils::strmatch(arg[iarg+3],"^v_")) { - ystr = utils::strdup(arg[iarg+3]+2); - } else yvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); - - } else error->all(FLERR,"Unknown fix indent cone argument: {}", arg[iarg+1]); - - if (utils::strmatch(arg[iarg+4],"^v_")) { - rlostr = utils::strdup(arg[iarg+4]+2); - } else rlovalue = utils::numeric(FLERR,arg[iarg+4],false,lmp); - - if (utils::strmatch(arg[iarg+5],"^v_")) { - rhistr = utils::strdup(arg[iarg+5]+2); - } else rhivalue = utils::numeric(FLERR,arg[iarg+5],false,lmp); - - if (utils::strmatch(arg[iarg+6],"^v_")) { - lostr = utils::strdup(arg[iarg+6]+2); - } else lovalue = utils::numeric(FLERR,arg[iarg+6],false,lmp); - - if (utils::strmatch(arg[iarg+7],"^v_")) { - histr = utils::strdup(arg[iarg+7]+2); - } else hivalue = utils::numeric(FLERR,arg[iarg+7],false,lmp); - - istyle = CONE; - iarg += 8; - - } else if (strcmp(arg[iarg],"units") == 0) { + if (strcmp(arg[iarg],"units") == 0) { if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix indent units", error); if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; @@ -644,16 +662,17 @@ void FixIndent::options(int narg, char **arg) else if (strcmp(arg[iarg+1],"out") == 0) side = OUTSIDE; else error->all(FLERR,"Unknown fix indent side argument: {}", arg[iarg+1]); iarg += 2; + } else error->all(FLERR,"Unknown fix indent argument: {}", arg[iarg]); } } /* ---------------------------------------------------------------------- - determines if a point is inside (true) or outside (false) of a cone + determines if a point is inside (true) or outside (false) of a cone ------------------------------------------------------------------------- */ -bool PointInsideCone(int dir, double *center, double lo, - double hi, double rlo, double rhi, double *x) +bool FixIndent::PointInsideCone(int dir, double *center, double lo, + double hi, double rlo, double rhi, double *x) { if ((x[dir] > hi) || (x[dir] < lo)) return false; @@ -669,10 +688,12 @@ bool PointInsideCone(int dir, double *center, double lo, } /* ---------------------------------------------------------------------- - distance between an exterior point and a cone + distance between an exterior point and a cone ------------------------------------------------------------------------- */ -void DistanceExteriorPoint(int dir, double *center, double lo, double hi, - double rlo, double rhi, double &x, double &y, double &z) + +void FixIndent::DistanceExteriorPoint(int dir, double *center, double lo, double hi, + double rlo, double rhi, + double &x, double &y, double &z) { double xp[3], nearest[3], corner1[3], corner2[3]; @@ -700,17 +721,21 @@ void DistanceExteriorPoint(int dir, double *center, double lo, double hi, corner4[dir] = hi; // initialize distance to a big number + double distsq = 1.0e20; // check the first triangle + point_on_line_segment(corner1, corner2, point, xp); distsq = closest(point, xp, nearest, distsq); // check the second triangle + point_on_line_segment(corner1, corner3, point, xp); distsq = closest(point, xp, nearest, distsq); // check the third triangle + point_on_line_segment(corner2, corner4, point, xp); distsq = closest(point, xp, nearest, distsq); @@ -722,12 +747,12 @@ void DistanceExteriorPoint(int dir, double *center, double lo, double hi, } /* ---------------------------------------------------------------------- - distance between an interior point and a cone + distance between an interior point and a cone ------------------------------------------------------------------------- */ -void DistanceInteriorPoint(int dir, double *center, - double lo, double hi, double rlo, double rhi, double &x, - double &y, double &z) +void FixIndent::DistanceInteriorPoint(int dir, double *center, + double lo, double hi, double rlo, double rhi, double &x, + double &y, double &z) { double r, dist_disk, dist_surf; double surflo[3], surfhi[3], xs[3]; @@ -735,6 +760,7 @@ void DistanceInteriorPoint(int dir, double *center, double point[3] {0.0, 0.0, 0.0}; // initial check with the two disks + if ( (initial_point[dir] - lo) < (hi - initial_point[dir]) ) { dist_disk = (initial_point[dir] - lo) * (initial_point[dir] - lo); point[dir] = initial_point[dir] - lo; @@ -744,6 +770,7 @@ void DistanceInteriorPoint(int dir, double *center, } // check with the points in the conical surface + double del[3] {x - center[0], y - center[1], z - center[2]}; del[dir] = 0.0; r = sqrt(del[0] * del[0] + del[1] * del[1] + del[2] * del[2]); @@ -776,10 +803,10 @@ void DistanceInteriorPoint(int dir, double *center, } /* ---------------------------------------------------------------------- - helper function extracted from region.cpp + helper function extracted from region.cpp ------------------------------------------------------------------------- */ -void point_on_line_segment(double *a, double *b, double *c, double *d) +void FixIndent::point_on_line_segment(double *a, double *b, double *c, double *d) { double ba[3], ca[3]; @@ -802,10 +829,10 @@ void point_on_line_segment(double *a, double *b, double *c, double *d) } /* ---------------------------------------------------------------------- - helper function extracted from region_cone.cpp + helper function extracted from region_cone.cpp ------------------------------------------------------------------------- */ -double closest(double *x, double *near, double *nearest, double dsq) +double FixIndent::closest(double *x, double *near, double *nearest, double dsq) { double dx = x[0] - near[0]; double dy = x[1] - near[1]; diff --git a/src/fix_indent.h b/src/fix_indent.h index 6f33f6fbb1..202a138729 100644 --- a/src/fix_indent.h +++ b/src/fix_indent.h @@ -53,7 +53,20 @@ class FixIndent : public Fix { int rlovar, rhivar, lovar, hivar; double rlovalue, rhivalue, lovalue, hivalue; + // methods for argument + + int geometry(int, char **); void options(int, char **); + + // methods for conical indenter + + bool PointInsideCone(int, double *, double, double, double, double, double *); + void DistanceExteriorPoint(int, double *, double, double, double, double, + double &, double &, double &); + void DistanceInteriorPoint(int, double *, double, double, double, double, + double &, double &, double &); + void point_on_line_segment(double *, double *, double *, double *); + double closest(double *, double *, double *, double); }; } // namespace LAMMPS_NS From 363db81be1f03351353d87f5cf459eaf52730d9f Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 26 Feb 2024 16:05:19 -0700 Subject: [PATCH 07/18] tweak a comment --- src/fix_indent.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fix_indent.h b/src/fix_indent.h index 202a138729..7224908390 100644 --- a/src/fix_indent.h +++ b/src/fix_indent.h @@ -53,7 +53,7 @@ class FixIndent : public Fix { int rlovar, rhivar, lovar, hivar; double rlovalue, rhivalue, lovalue, hivalue; - // methods for argument + // methods for argument parsing int geometry(int, char **); void options(int, char **); From 5ad4545273a66c36b4da032287b79f19fdcf0068 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 26 Feb 2024 16:07:14 -0700 Subject: [PATCH 08/18] fix where initialization of options is done --- src/fix_indent.cpp | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp index 8d450bee75..3599c4ddaa 100644 --- a/src/fix_indent.cpp +++ b/src/fix_indent.cpp @@ -347,7 +347,6 @@ void FixIndent::post_force(int /*vflag*/) double radiuslo { rlostr ? input->variable->compute_equal(rlovar) : rlovalue }; if (radiuslo < 0.0) error->all(FLERR, "Illegal fix indent cone lower radius: {}", radiuslo); - double radiushi { rhistr ? input->variable->compute_equal(rhivar) : rhivalue }; if (radiushi < 0.0) error->all(FLERR, "Illegal fix indent cone high radius: {}", radiushi); @@ -485,8 +484,6 @@ int FixIndent::geometry(int narg, char **arg) istyle = NONE; xstr = ystr = zstr = rstr = pstr = nullptr; xvalue = yvalue = zvalue = rvalue = pvalue = 0.0; - scaleflag = 1; - side = OUTSIDE; // sphere @@ -559,33 +556,27 @@ int FixIndent::geometry(int narg, char **arg) if (strcmp(arg[1],"x") == 0) { cdim = 0; - if (utils::strmatch(arg[2],"^v_")) { ystr = utils::strdup(arg[2]+2); } else yvalue = utils::numeric(FLERR,arg[2],false,lmp); - if (utils::strmatch(arg[3],"^v_")) { zstr = utils::strdup(arg[3]+2); } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); } else if (strcmp(arg[1],"y") == 0) { cdim = 1; - if (utils::strmatch(arg[2],"^v_")) { xstr = utils::strdup(arg[2]+2); } else xvalue = utils::numeric(FLERR,arg[2],false,lmp); - if (utils::strmatch(arg[3],"^v_")) { zstr = utils::strdup(arg[3]+2); } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); } else if (strcmp(arg[1],"z") == 0) { cdim = 2; - if (utils::strmatch(arg[2],"^v_")) { xstr = utils::strdup(arg[2]+2); } else xvalue = utils::numeric(FLERR,arg[2],false,lmp); - if (utils::strmatch(arg[3],"^v_")) { ystr = utils::strdup(arg[3]+2); } else yvalue = utils::numeric(FLERR,arg[3],false,lmp); @@ -595,15 +586,12 @@ int FixIndent::geometry(int narg, char **arg) if (utils::strmatch(arg[4],"^v_")) { rlostr = utils::strdup(arg[4]+2); } else rlovalue = utils::numeric(FLERR,arg[4],false,lmp); - if (utils::strmatch(arg[5],"^v_")) { rhistr = utils::strdup(arg[5]+2); } else rhivalue = utils::numeric(FLERR,arg[5],false,lmp); - if (utils::strmatch(arg[6],"^v_")) { lostr = utils::strdup(arg[6]+2); } else lovalue = utils::numeric(FLERR,arg[6],false,lmp); - if (utils::strmatch(arg[7],"^v_")) { histr = utils::strdup(arg[7]+2); } else hivalue = utils::numeric(FLERR,arg[7],false,lmp); @@ -646,6 +634,9 @@ int FixIndent::geometry(int narg, char **arg) void FixIndent::options(int narg, char **arg) { + scaleflag = 1; + side = OUTSIDE; + int iarg = 0; while (iarg < narg) { From ba5c1a4ac39f84e38d21c9d25b24b6ce8c8e14e4 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Fri, 8 Mar 2024 15:35:48 +0100 Subject: [PATCH 09/18] minor modifications in fix_indent.rst --- doc/src/fix_indent.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/src/fix_indent.rst b/doc/src/fix_indent.rst index 5658c06373..e041f9f29b 100644 --- a/doc/src/fix_indent.rst +++ b/doc/src/fix_indent.rst @@ -44,8 +44,8 @@ Syntax .. parsed-literal:: *side* value = *in* or *out* - *in* = the indenter acts on particles inside the sphere or cylinder - *out* = the indenter acts on particles outside the sphere or cylinder + *in* = the indenter acts on particles inside the sphere or cylinder or cone + *out* = the indenter acts on particles outside the sphere or cylinder or cone *units* value = *lattice* or *box* lattice = the geometry is defined in lattice units box = the geometry is defined in simulation box units @@ -64,7 +64,7 @@ Description Insert an indenter within a simulation box. The indenter repels all atoms in the group that touch it, so it can be used to push into a -material or as an obstacle in a flow. Or it can be used as a +material or as an obstacle in a flow. Alternatively, it can be used as a constraining wall around a simulation; see the discussion of the *side* keyword below. From ebe57ce9eb329018b64ff503651c81b721091adc Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Fri, 8 Mar 2024 15:45:46 +0100 Subject: [PATCH 10/18] removing whitespaces - fix_indent.cpp --- src/fix_indent.cpp | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp index 3599c4ddaa..87ed5091bd 100644 --- a/src/fix_indent.cpp +++ b/src/fix_indent.cpp @@ -65,7 +65,7 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : int iarg = geometry(narg-4,&arg[4]) + 4; options(narg-iarg,&arg[iarg]); - + // setup scaling const double xscale { scaleflag ? domain->lattice->xlattice : 1.0}; @@ -268,7 +268,7 @@ void FixIndent::post_force(int /*vflag*/) if (istyle == SPHERE) { // remap indenter center into periodic box - + domain->remap(ctr); double radius { rstr ? input->variable->compute_equal(rvar) : rvalue}; @@ -388,7 +388,7 @@ void FixIndent::post_force(int /*vflag*/) // compute the force from the center of the cone // this is different from how it is done in fix wall/region - + dr = sqrt(x0[0] * x0[0] + x0[1] * x0[1] + x0[2] * x0[2]); int force_sign = { point_inside_cone ? 1 : -1 }; @@ -509,7 +509,7 @@ int FixIndent::geometry(int narg, char **arg) } // cylinder - + if (strcmp(arg[0],"cylinder") == 0) { if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword"); if (5 > narg) utils::missing_cmd_args(FLERR, "fix indent cylinder", error); @@ -539,21 +539,21 @@ int FixIndent::geometry(int narg, char **arg) ystr = utils::strdup(arg[3]+2); } else yvalue = utils::numeric(FLERR,arg[3],false,lmp); } else error->all(FLERR,"Unknown fix indent cylinder argument: {}", arg[1]); - + if (utils::strmatch(arg[4],"^v_")) { rstr = utils::strdup(arg[4]+2); } else rvalue = utils::numeric(FLERR,arg[4],false,lmp); - + istyle = CYLINDER; return 5; } // cone - + if (strcmp(arg[0],"cone") == 0) { if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword"); if (8 > narg) utils::missing_cmd_args(FLERR, "fix indent cone", error); - + if (strcmp(arg[1],"x") == 0) { cdim = 0; if (utils::strmatch(arg[2],"^v_")) { @@ -562,7 +562,7 @@ int FixIndent::geometry(int narg, char **arg) if (utils::strmatch(arg[3],"^v_")) { zstr = utils::strdup(arg[3]+2); } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); - + } else if (strcmp(arg[1],"y") == 0) { cdim = 1; if (utils::strmatch(arg[2],"^v_")) { @@ -571,7 +571,7 @@ int FixIndent::geometry(int narg, char **arg) if (utils::strmatch(arg[3],"^v_")) { zstr = utils::strdup(arg[3]+2); } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); - + } else if (strcmp(arg[1],"z") == 0) { cdim = 2; if (utils::strmatch(arg[2],"^v_")) { @@ -580,9 +580,9 @@ int FixIndent::geometry(int narg, char **arg) if (utils::strmatch(arg[3],"^v_")) { ystr = utils::strdup(arg[3]+2); } else yvalue = utils::numeric(FLERR,arg[3],false,lmp); - + } else error->all(FLERR,"Unknown fix indent cone argument: {}", arg[1]); - + if (utils::strmatch(arg[4],"^v_")) { rlostr = utils::strdup(arg[4]+2); } else rlovalue = utils::numeric(FLERR,arg[4],false,lmp); @@ -638,7 +638,7 @@ void FixIndent::options(int narg, char **arg) side = OUTSIDE; int iarg = 0; - + while (iarg < narg) { if (strcmp(arg[iarg],"units") == 0) { if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix indent units", error); @@ -653,7 +653,7 @@ void FixIndent::options(int narg, char **arg) else if (strcmp(arg[iarg+1],"out") == 0) side = OUTSIDE; else error->all(FLERR,"Unknown fix indent side argument: {}", arg[iarg+1]); iarg += 2; - + } else error->all(FLERR,"Unknown fix indent argument: {}", arg[iarg]); } } @@ -712,21 +712,21 @@ void FixIndent::DistanceExteriorPoint(int dir, double *center, double lo, double corner4[dir] = hi; // initialize distance to a big number - + double distsq = 1.0e20; // check the first triangle - + point_on_line_segment(corner1, corner2, point, xp); distsq = closest(point, xp, nearest, distsq); // check the second triangle - + point_on_line_segment(corner1, corner3, point, xp); distsq = closest(point, xp, nearest, distsq); // check the third triangle - + point_on_line_segment(corner2, corner4, point, xp); distsq = closest(point, xp, nearest, distsq); @@ -751,7 +751,7 @@ void FixIndent::DistanceInteriorPoint(int dir, double *center, double point[3] {0.0, 0.0, 0.0}; // initial check with the two disks - + if ( (initial_point[dir] - lo) < (hi - initial_point[dir]) ) { dist_disk = (initial_point[dir] - lo) * (initial_point[dir] - lo); point[dir] = initial_point[dir] - lo; @@ -761,7 +761,7 @@ void FixIndent::DistanceInteriorPoint(int dir, double *center, } // check with the points in the conical surface - + double del[3] {x - center[0], y - center[1], z - center[2]}; del[dir] = 0.0; r = sqrt(del[0] * del[0] + del[1] * del[1] + del[2] * del[2]); From 01628dfc614799adb2efb1c5414de0e40d0cf7d9 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Fri, 8 Mar 2024 15:46:27 +0100 Subject: [PATCH 11/18] removing whitespaces in fix_indent.h --- src/fix_indent.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/fix_indent.h b/src/fix_indent.h index 7224908390..37e1623df9 100644 --- a/src/fix_indent.h +++ b/src/fix_indent.h @@ -54,12 +54,12 @@ class FixIndent : public Fix { double rlovalue, rhivalue, lovalue, hivalue; // methods for argument parsing - + int geometry(int, char **); void options(int, char **); // methods for conical indenter - + bool PointInsideCone(int, double *, double, double, double, double, double *); void DistanceExteriorPoint(int, double *, double, double, double, double, double &, double &, double &); From 94e9fe5df38db71b01080a970c2ce2d9175ad74e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 8 Mar 2024 11:08:06 -0500 Subject: [PATCH 12/18] flag error with explanation when removed keyword "reax/c" is used instead of a file not found --- src/REAXFF/fix_qeq_reaxff.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/REAXFF/fix_qeq_reaxff.cpp b/src/REAXFF/fix_qeq_reaxff.cpp index ab561de0a7..2c3089b5e8 100644 --- a/src/REAXFF/fix_qeq_reaxff.cpp +++ b/src/REAXFF/fix_qeq_reaxff.cpp @@ -141,7 +141,7 @@ FixQEqReaxFF::FixQEqReaxFF(LAMMPS *lmp, int narg, char **arg) : // perform initial allocation of atom-based arrays // register with Atom class - reaxff = dynamic_cast(force->pair_match("^reax..",0)); + reaxff = dynamic_cast(force->pair_match("^reaxff",0)); s_hist = t_hist = nullptr; atom->add_callback(Atom::GROW); @@ -217,6 +217,8 @@ void FixQEqReaxFF::pertype_parameters(char *arg) if (chi == nullptr || eta == nullptr || gamma == nullptr) error->all(FLERR, "Fix qeq/reaxff could not extract params from pair reaxff"); return; + } else if (utils::strmatch(arg,"^reax/c")) { + error->all(FLERR, "Fix qeq/reaxff keyword 'reax/c' is obsolete; please use 'reaxff'"); } reaxflag = 0; From bd0eb1ec843b82525dfee31e6a904b8e7bb32a68 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 8 Mar 2024 11:13:15 -0500 Subject: [PATCH 13/18] silence compiler warning, cosmetic --- src/EXTRA-FIX/fix_nonaffine_displacement.cpp | 6 ++++-- src/EXTRA-FIX/fix_nonaffine_displacement.h | 1 - 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp index ee14822e98..eaf45f4e59 100644 --- a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp +++ b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp @@ -65,7 +65,8 @@ static const char cite_nonaffine_d2min[] = /* ---------------------------------------------------------------------- */ FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), id_fix(nullptr), X(nullptr), Y(nullptr), F(nullptr), norm(nullptr), D2min(nullptr) + Fix(lmp, narg, arg), id_fix(nullptr), fix(nullptr), D2min(nullptr), X(nullptr), Y(nullptr), + F(nullptr), norm(nullptr) { if (narg < 4) utils::missing_cmd_args(FLERR,"fix nonaffine/displacement", error); @@ -85,7 +86,8 @@ FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char * } else if (strcmp(arg[iarg + 1], "radius") == 0) { cut_style = RADIUS; } else if (strcmp(arg[iarg + 1], "custom") == 0) { - if (iarg + 2 > narg) utils::missing_cmd_args(FLERR,"fix nonaffine/displacement custom", error); + if (iarg + 2 > narg) + utils::missing_cmd_args(FLERR,"fix nonaffine/displacement custom", error); if ((neighbor->style == Neighbor::MULTI) || (neighbor->style == Neighbor::MULTI_OLD)) error->all(FLERR, "Fix nonaffine/displacement with custom cutoff requires neighbor style 'bin' or 'nsq'"); cut_style = CUSTOM; diff --git a/src/EXTRA-FIX/fix_nonaffine_displacement.h b/src/EXTRA-FIX/fix_nonaffine_displacement.h index 79dbdabf49..c7177bd3d9 100644 --- a/src/EXTRA-FIX/fix_nonaffine_displacement.h +++ b/src/EXTRA-FIX/fix_nonaffine_displacement.h @@ -57,7 +57,6 @@ class FixNonaffineDisplacement : public Fix { class NeighList *list; // half neighbor list - void integrate_velocity(); void calculate_D2Min(); void save_reference_state(); From e7d77b62445d9d85c19c3cd3954efdf64fcafcee Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 8 Mar 2024 12:23:02 -0500 Subject: [PATCH 14/18] enable and apply clang-format, revert to some older code constructs where equivalent --- src/fix_indent.cpp | 475 ++++++++++++++++++++++++--------------------- 1 file changed, 251 insertions(+), 224 deletions(-) diff --git a/src/fix_indent.cpp b/src/fix_indent.cpp index 87ed5091bd..a8b14940f4 100644 --- a/src/fix_indent.cpp +++ b/src/fix_indent.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -35,15 +34,14 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum{NONE, SPHERE, CYLINDER, PLANE, CONE}; -enum{INSIDE, OUTSIDE}; +enum { NONE, SPHERE, CYLINDER, PLANE, CONE }; +enum { INSIDE, OUTSIDE }; /* ---------------------------------------------------------------------- */ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), - xstr(nullptr), ystr(nullptr), zstr(nullptr), rstr(nullptr), pstr(nullptr), - rlostr(nullptr), rhistr(nullptr), lostr(nullptr), histr(nullptr) + Fix(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr), rstr(nullptr), pstr(nullptr), + rlostr(nullptr), rhistr(nullptr), lostr(nullptr), histr(nullptr) { if (narg < 4) utils::missing_cmd_args(FLERR, "fix indent", error); @@ -57,20 +55,20 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : respa_level_support = 1; ilevel_respa = 0; - k = utils::numeric(FLERR,arg[3],false,lmp); + k = utils::numeric(FLERR, arg[3], false, lmp); if (k < 0.0) error->all(FLERR, "Illegal fix indent force constant: {}", k); - k3 = k/3.0; + k3 = k / 3.0; // read geometry of indenter and optional args - int iarg = geometry(narg-4,&arg[4]) + 4; - options(narg-iarg,&arg[iarg]); + int iarg = geometry(narg - 4, &arg[4]) + 4; + options(narg - iarg, &arg[iarg]); // setup scaling - const double xscale { scaleflag ? domain->lattice->xlattice : 1.0}; - const double yscale { scaleflag ? domain->lattice->ylattice : 1.0}; - const double zscale { scaleflag ? domain->lattice->zlattice : 1.0}; + const double xscale{scaleflag ? domain->lattice->xlattice : 1.0}; + const double yscale{scaleflag ? domain->lattice->ylattice : 1.0}; + const double zscale{scaleflag ? domain->lattice->zlattice : 1.0}; // apply scaling factors to geometry @@ -79,7 +77,7 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : if (!ystr) yvalue *= yscale; if (!zstr) zvalue *= zscale; if (!rstr) rvalue *= xscale; - + } else if (istyle == CONE) { if (!xstr) xvalue *= xscale; if (!ystr) yvalue *= yscale; @@ -104,11 +102,15 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : if (!histr) hivalue *= scaling_factor; } else if (istyle == PLANE) { - if (cdim == 0 && !pstr) pvalue *= xscale; - else if (cdim == 1 && !pstr) pvalue *= yscale; - else if (cdim == 2 && !pstr) pvalue *= zscale; + if (cdim == 0 && !pstr) + pvalue *= xscale; + else if (cdim == 1 && !pstr) + pvalue *= yscale; + else if (cdim == 2 && !pstr) + pvalue *= zscale; - } else error->all(FLERR,"Unknown fix indent keyword: {}", istyle); + } else + error->all(FLERR, "Unknown fix indent keyword: {}", istyle); varflag = 0; if (xstr || ystr || zstr || rstr || pstr || rlostr || rhistr || lostr || histr) varflag = 1; @@ -121,15 +123,15 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) : FixIndent::~FixIndent() { - delete [] xstr; - delete [] ystr; - delete [] zstr; - delete [] rstr; - delete [] pstr; - delete [] rlostr; - delete [] rhistr; - delete [] lostr; - delete [] histr; + delete[] xstr; + delete[] ystr; + delete[] zstr; + delete[] rstr; + delete[] pstr; + delete[] rlostr; + delete[] rhistr; + delete[] lostr; + delete[] histr; } /* ---------------------------------------------------------------------- */ @@ -149,71 +151,62 @@ void FixIndent::init() { if (xstr) { xvar = input->variable->find(xstr); - if (xvar < 0) - error->all(FLERR,"Variable {} for fix indent does not exist", xstr); + if (xvar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", xstr); if (!input->variable->equalstyle(xvar)) - error->all(FLERR,"Variable {} for fix indent is invalid style", xstr); + error->all(FLERR, "Variable {} for fix indent is invalid style", xstr); } if (ystr) { yvar = input->variable->find(ystr); - if (yvar < 0) - error->all(FLERR,"Variable {} for fix indent does not exist", ystr); + if (yvar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", ystr); if (!input->variable->equalstyle(yvar)) - error->all(FLERR,"Variable {} for fix indent is invalid style", ystr); + error->all(FLERR, "Variable {} for fix indent is invalid style", ystr); } if (zstr) { zvar = input->variable->find(zstr); - if (zvar < 0) - error->all(FLERR,"Variable {} for fix indent does not exist", zstr); + if (zvar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", zstr); if (!input->variable->equalstyle(zvar)) - error->all(FLERR,"Variable {} for fix indent is invalid style", zstr); + error->all(FLERR, "Variable {} for fix indent is invalid style", zstr); } if (rstr) { rvar = input->variable->find(rstr); - if (rvar < 0) - error->all(FLERR,"Variable {} for fix indent does not exist", rstr); + if (rvar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", rstr); if (!input->variable->equalstyle(rvar)) - error->all(FLERR,"Variable {} for fix indent is invalid style", rstr); + error->all(FLERR, "Variable {} for fix indent is invalid style", rstr); } if (pstr) { pvar = input->variable->find(pstr); - if (pvar < 0) - error->all(FLERR,"Variable {} for fix indent does not exist", pstr); + if (pvar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", pstr); if (!input->variable->equalstyle(pvar)) - error->all(FLERR,"Variable {} for fix indent is invalid style", pstr); + error->all(FLERR, "Variable {} for fix indent is invalid style", pstr); } if (rlostr) { rlovar = input->variable->find(rlostr); - if (rlovar < 0) - error->all(FLERR,"Variable {} for fix indent does not exist", rlostr); + if (rlovar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", rlostr); if (!input->variable->equalstyle(rlovar)) - error->all(FLERR,"Variable {} for fix indent is invalid style", rlostr); + error->all(FLERR, "Variable {} for fix indent is invalid style", rlostr); } if (rhistr) { rhivar = input->variable->find(rhistr); - if (rhivar < 0) - error->all(FLERR,"Variable {} for fix indent does not exist", rhistr); + if (rhivar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", rhistr); if (!input->variable->equalstyle(rhivar)) - error->all(FLERR,"Variable {} for fix indent is invalid style", rhistr); + error->all(FLERR, "Variable {} for fix indent is invalid style", rhistr); } if (lostr) { lovar = input->variable->find(lostr); - if (lovar < 0) - error->all(FLERR,"Variable {} for fix indent does not exist", lostr); + if (lovar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", lostr); if (!input->variable->equalstyle(lovar)) - error->all(FLERR,"Variable {} for fix indent is invalid style", lostr); + error->all(FLERR, "Variable {} for fix indent is invalid style", lostr); } if (histr) { hivar = input->variable->find(histr); - if (hivar < 0) - error->all(FLERR,"Variable {} for fix indent does not exist", histr); + if (hivar < 0) error->all(FLERR, "Variable {} for fix indent does not exist", histr); if (!input->variable->equalstyle(hivar)) - error->all(FLERR,"Variable {} for fix indent is invalid style", histr); + error->all(FLERR, "Variable {} for fix indent is invalid style", histr); } - if (utils::strmatch(update->integrate_style,"^respa")) { - ilevel_respa = (dynamic_cast(update->integrate))->nlevels-1; - if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa); + if (utils::strmatch(update->integrate_style, "^respa")) { + ilevel_respa = (dynamic_cast(update->integrate))->nlevels - 1; + if (respa_level >= 0) ilevel_respa = MIN(respa_level, ilevel_respa); } } @@ -221,11 +214,11 @@ void FixIndent::init() void FixIndent::setup(int vflag) { - if (utils::strmatch(update->integrate_style,"^verlet")) + if (utils::strmatch(update->integrate_style, "^verlet")) post_force(vflag); else { (dynamic_cast(update->integrate))->copy_flevel_f(ilevel_respa); - post_force_respa(vflag,ilevel_respa,0); + post_force_respa(vflag, ilevel_respa, 0); (dynamic_cast(update->integrate))->copy_f_flevel(ilevel_respa); } } @@ -250,8 +243,8 @@ void FixIndent::post_force(int /*vflag*/) indenter[0] = indenter[1] = indenter[2] = indenter[3] = 0.0; // ctr = current indenter centerz - - double ctr[3] {xvalue, yvalue, zvalue}; + + double ctr[3] = {xvalue, yvalue, zvalue}; if (xstr) ctr[0] = input->variable->compute_equal(xvar); if (ystr) ctr[1] = input->variable->compute_equal(yvar); if (zstr) ctr[2] = input->variable->compute_equal(zvar); @@ -271,7 +264,7 @@ void FixIndent::post_force(int /*vflag*/) domain->remap(ctr); - double radius { rstr ? input->variable->compute_equal(rvar) : rvalue}; + double radius = rstr ? input->variable->compute_equal(rvar) : rvalue; if (radius < 0.0) error->all(FLERR, "Illegal fix indent sphere radius: {}", radius); for (int i = 0; i < nlocal; i++) @@ -279,29 +272,29 @@ void FixIndent::post_force(int /*vflag*/) delx = x[i][0] - ctr[0]; dely = x[i][1] - ctr[1]; delz = x[i][2] - ctr[2]; - domain->minimum_image(delx,dely,delz); - r = sqrt(delx*delx + dely*dely + delz*delz); + domain->minimum_image(delx, dely, delz); + r = sqrt(delx * delx + dely * dely + delz * delz); if (side == OUTSIDE) { dr = r - radius; - fmag = k*dr*dr; + fmag = k * dr * dr; } else { dr = radius - r; - fmag = -k*dr*dr; + fmag = -k * dr * dr; } if (dr >= 0.0) continue; - fx = delx*fmag/r; - fy = dely*fmag/r; - fz = delz*fmag/r; + fx = delx * fmag / r; + fy = dely * fmag / r; + fz = delz * fmag / r; f[i][0] += fx; f[i][1] += fy; f[i][2] += fz; - indenter[0] -= k3 * dr*dr*dr; + indenter[0] -= k3 * dr * dr * dr; indenter[1] -= fx; indenter[2] -= fy; indenter[3] -= fz; } - // cylindrical indenter + // cylindrical indenter } else if (istyle == CYLINDER) { @@ -312,46 +305,46 @@ void FixIndent::post_force(int /*vflag*/) ctr[cdim] = domain->boxlo[cdim]; domain->remap(ctr); - double radius { rstr ? input->variable->compute_equal(rvar) : rvalue}; + double radius{rstr ? input->variable->compute_equal(rvar) : rvalue}; if (radius < 0.0) error->all(FLERR, "Illegal fix indent cylinder radius: {}", radius); for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - double del[3] {x[i][0] - ctr[0], x[i][1] - ctr[1], x[i][2] - ctr[2]}; + double del[3] = {x[i][0] - ctr[0], x[i][1] - ctr[1], x[i][2] - ctr[2]}; del[cdim] = 0; domain->minimum_image(del[0], del[1], del[2]); - r = sqrt(del[0]*del[0] + del[1]*del[1] + del[2]*del[2]); + r = sqrt(del[0] * del[0] + del[1] * del[1] + del[2] * del[2]); if (side == OUTSIDE) { dr = r - radius; - fmag = k*dr*dr; + fmag = k * dr * dr; } else { dr = radius - r; - fmag = -k*dr*dr; + fmag = -k * dr * dr; } if (dr >= 0.0) continue; - fx = del[0]*fmag/r; - fy = del[1]*fmag/r; - fz = del[2]*fmag/r; + fx = del[0] * fmag / r; + fy = del[1] * fmag / r; + fz = del[2] * fmag / r; f[i][0] += fx; f[i][1] += fy; f[i][2] += fz; - indenter[0] -= k3 * dr*dr*dr; + indenter[0] -= k3 * dr * dr * dr; indenter[1] -= fx; indenter[2] -= fy; indenter[3] -= fz; } - // conical indenter + // conical indenter } else if (istyle == CONE) { - double radiuslo { rlostr ? input->variable->compute_equal(rlovar) : rlovalue }; + double radiuslo{rlostr ? input->variable->compute_equal(rlovar) : rlovalue}; if (radiuslo < 0.0) error->all(FLERR, "Illegal fix indent cone lower radius: {}", radiuslo); - double radiushi { rhistr ? input->variable->compute_equal(rhivar) : rhivalue }; + double radiushi{rhistr ? input->variable->compute_equal(rhivar) : rhivalue}; if (radiushi < 0.0) error->all(FLERR, "Illegal fix indent cone high radius: {}", radiushi); - double initial_lo { lostr ? input->variable->compute_equal(lovar) : lovalue }; - double initial_hi { histr ? input->variable->compute_equal(hivar) : hivalue }; + double initial_lo{lostr ? input->variable->compute_equal(lovar) : lovalue}; + double initial_hi{histr ? input->variable->compute_equal(hivar) : hivalue}; ctr[cdim] = 0.5 * (initial_hi + initial_lo); @@ -368,18 +361,18 @@ void FixIndent::post_force(int /*vflag*/) delz = x[i][2] - ctr[2]; domain->minimum_image(delx, dely, delz); - double x0[3] {delx + ctr[0], dely + ctr[1], delz + ctr[2]}; + double x0[3] = {delx + ctr[0], dely + ctr[1], delz + ctr[2]}; r = sqrt(delx * delx + dely * dely + delz * delz); // check if particle is inside or outside the cone - + bool point_inside_cone = PointInsideCone(cdim, ctr, lo, hi, radiuslo, radiushi, x0); if (side == INSIDE && point_inside_cone) continue; if (side == OUTSIDE && !point_inside_cone) continue; // find the distance between the point and the cone - + if (point_inside_cone) { DistanceInteriorPoint(cdim, ctr, lo, hi, radiuslo, radiushi, x0[0], x0[1], x0[2]); } else { @@ -391,12 +384,12 @@ void FixIndent::post_force(int /*vflag*/) dr = sqrt(x0[0] * x0[0] + x0[1] * x0[1] + x0[2] * x0[2]); - int force_sign = { point_inside_cone ? 1 : -1 }; + int force_sign = {point_inside_cone ? 1 : -1}; fmag = force_sign * k * dr * dr; - fx = delx*fmag/r; - fy = dely*fmag/r; - fz = delz*fmag/r; + fx = delx * fmag / r; + fy = dely * fmag / r; + fz = delz * fmag / r; f[i][0] += fx; f[i][1] += fy; f[i][2] += fz; @@ -407,13 +400,13 @@ void FixIndent::post_force(int /*vflag*/) } } - // planar indenter + // planar indenter } else { // plane = current plane position - double plane { pstr ? input->variable->compute_equal(pvar) : pvalue}; + double plane{pstr ? input->variable->compute_equal(pvar) : pvalue}; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { @@ -422,7 +415,7 @@ void FixIndent::post_force(int /*vflag*/) fmag = -planeside * k * dr * dr; f[i][cdim] += fmag; indenter[0] -= k3 * dr * dr * dr; - indenter[cdim+1] -= fmag; + indenter[cdim + 1] -= fmag; } } @@ -452,7 +445,7 @@ double FixIndent::compute_scalar() // only sum across procs one time if (indenter_flag == 0) { - MPI_Allreduce(indenter,indenter_all,4,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(indenter, indenter_all, 4, MPI_DOUBLE, MPI_SUM, world); indenter_flag = 1; } return indenter_all[0]; @@ -467,10 +460,10 @@ double FixIndent::compute_vector(int n) // only sum across procs one time if (indenter_flag == 0) { - MPI_Allreduce(indenter,indenter_all,4,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(indenter, indenter_all, 4, MPI_DOUBLE, MPI_SUM, world); indenter_flag = 1; } - return indenter_all[n+1]; + return indenter_all[n + 1]; } /* ---------------------------------------------------------------------- @@ -486,63 +479,75 @@ int FixIndent::geometry(int narg, char **arg) xvalue = yvalue = zvalue = rvalue = pvalue = 0.0; // sphere - - if (strcmp(arg[0],"sphere") == 0) { + + if (strcmp(arg[0], "sphere") == 0) { if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword"); if (5 > narg) utils::missing_cmd_args(FLERR, "fix indent sphere", error); - if (utils::strmatch(arg[1],"^v_")) { - xstr = utils::strdup(arg[1]+2); - } else xvalue = utils::numeric(FLERR,arg[1],false,lmp); - if (utils::strmatch(arg[2],"^v_")) { - ystr = utils::strdup(arg[2]+2); - } else yvalue = utils::numeric(FLERR,arg[2],false,lmp); - if (utils::strmatch(arg[3],"^v_")) { - zstr = utils::strdup(arg[3]+2); - } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); - if (utils::strmatch(arg[4],"^v_")) { - rstr = utils::strdup(arg[4]+2); - } else rvalue = utils::numeric(FLERR,arg[4],false,lmp); - + if (utils::strmatch(arg[1], "^v_")) { + xstr = utils::strdup(arg[1] + 2); + } else + xvalue = utils::numeric(FLERR, arg[1], false, lmp); + if (utils::strmatch(arg[2], "^v_")) { + ystr = utils::strdup(arg[2] + 2); + } else + yvalue = utils::numeric(FLERR, arg[2], false, lmp); + if (utils::strmatch(arg[3], "^v_")) { + zstr = utils::strdup(arg[3] + 2); + } else + zvalue = utils::numeric(FLERR, arg[3], false, lmp); + if (utils::strmatch(arg[4], "^v_")) { + rstr = utils::strdup(arg[4] + 2); + } else + rvalue = utils::numeric(FLERR, arg[4], false, lmp); + istyle = SPHERE; return 5; } // cylinder - if (strcmp(arg[0],"cylinder") == 0) { + if (strcmp(arg[0], "cylinder") == 0) { if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword"); if (5 > narg) utils::missing_cmd_args(FLERR, "fix indent cylinder", error); - if (strcmp(arg[1],"x") == 0) { + if (strcmp(arg[1], "x") == 0) { cdim = 0; - if (utils::strmatch(arg[2],"^v_")) { - ystr = utils::strdup(arg[2]+2); - } else yvalue = utils::numeric(FLERR,arg[2],false,lmp); - if (utils::strmatch(arg[3],"^v_")) { - zstr = utils::strdup(arg[3]+2); - } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); - } else if (strcmp(arg[1],"y") == 0) { + if (utils::strmatch(arg[2], "^v_")) { + ystr = utils::strdup(arg[2] + 2); + } else + yvalue = utils::numeric(FLERR, arg[2], false, lmp); + if (utils::strmatch(arg[3], "^v_")) { + zstr = utils::strdup(arg[3] + 2); + } else + zvalue = utils::numeric(FLERR, arg[3], false, lmp); + } else if (strcmp(arg[1], "y") == 0) { cdim = 1; - if (utils::strmatch(arg[2],"^v_")) { - xstr = utils::strdup(arg[2]+2); - } else xvalue = utils::numeric(FLERR,arg[2],false,lmp); - if (utils::strmatch(arg[3],"^v_")) { - zstr = utils::strdup(arg[3]+2); - } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); - } else if (strcmp(arg[1],"z") == 0) { + if (utils::strmatch(arg[2], "^v_")) { + xstr = utils::strdup(arg[2] + 2); + } else + xvalue = utils::numeric(FLERR, arg[2], false, lmp); + if (utils::strmatch(arg[3], "^v_")) { + zstr = utils::strdup(arg[3] + 2); + } else + zvalue = utils::numeric(FLERR, arg[3], false, lmp); + } else if (strcmp(arg[1], "z") == 0) { cdim = 2; - if (utils::strmatch(arg[2],"^v_")) { - xstr = utils::strdup(arg[2]+2); - } else xvalue = utils::numeric(FLERR,arg[2],false,lmp); - if (utils::strmatch(arg[3],"^v_")) { - ystr = utils::strdup(arg[3]+2); - } else yvalue = utils::numeric(FLERR,arg[3],false,lmp); - } else error->all(FLERR,"Unknown fix indent cylinder argument: {}", arg[1]); + if (utils::strmatch(arg[2], "^v_")) { + xstr = utils::strdup(arg[2] + 2); + } else + xvalue = utils::numeric(FLERR, arg[2], false, lmp); + if (utils::strmatch(arg[3], "^v_")) { + ystr = utils::strdup(arg[3] + 2); + } else + yvalue = utils::numeric(FLERR, arg[3], false, lmp); + } else + error->all(FLERR, "Unknown fix indent cylinder argument: {}", arg[1]); - if (utils::strmatch(arg[4],"^v_")) { - rstr = utils::strdup(arg[4]+2); - } else rvalue = utils::numeric(FLERR,arg[4],false,lmp); + if (utils::strmatch(arg[4], "^v_")) { + rstr = utils::strdup(arg[4] + 2); + } else + rvalue = utils::numeric(FLERR, arg[4], false, lmp); istyle = CYLINDER; return 5; @@ -550,81 +555,100 @@ int FixIndent::geometry(int narg, char **arg) // cone - if (strcmp(arg[0],"cone") == 0) { + if (strcmp(arg[0], "cone") == 0) { if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword"); if (8 > narg) utils::missing_cmd_args(FLERR, "fix indent cone", error); - if (strcmp(arg[1],"x") == 0) { + if (strcmp(arg[1], "x") == 0) { cdim = 0; - if (utils::strmatch(arg[2],"^v_")) { - ystr = utils::strdup(arg[2]+2); - } else yvalue = utils::numeric(FLERR,arg[2],false,lmp); - if (utils::strmatch(arg[3],"^v_")) { - zstr = utils::strdup(arg[3]+2); - } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); + if (utils::strmatch(arg[2], "^v_")) { + ystr = utils::strdup(arg[2] + 2); + } else + yvalue = utils::numeric(FLERR, arg[2], false, lmp); + if (utils::strmatch(arg[3], "^v_")) { + zstr = utils::strdup(arg[3] + 2); + } else + zvalue = utils::numeric(FLERR, arg[3], false, lmp); - } else if (strcmp(arg[1],"y") == 0) { + } else if (strcmp(arg[1], "y") == 0) { cdim = 1; - if (utils::strmatch(arg[2],"^v_")) { - xstr = utils::strdup(arg[2]+2); - } else xvalue = utils::numeric(FLERR,arg[2],false,lmp); - if (utils::strmatch(arg[3],"^v_")) { - zstr = utils::strdup(arg[3]+2); - } else zvalue = utils::numeric(FLERR,arg[3],false,lmp); + if (utils::strmatch(arg[2], "^v_")) { + xstr = utils::strdup(arg[2] + 2); + } else + xvalue = utils::numeric(FLERR, arg[2], false, lmp); + if (utils::strmatch(arg[3], "^v_")) { + zstr = utils::strdup(arg[3] + 2); + } else + zvalue = utils::numeric(FLERR, arg[3], false, lmp); - } else if (strcmp(arg[1],"z") == 0) { + } else if (strcmp(arg[1], "z") == 0) { cdim = 2; - if (utils::strmatch(arg[2],"^v_")) { - xstr = utils::strdup(arg[2]+2); - } else xvalue = utils::numeric(FLERR,arg[2],false,lmp); - if (utils::strmatch(arg[3],"^v_")) { - ystr = utils::strdup(arg[3]+2); - } else yvalue = utils::numeric(FLERR,arg[3],false,lmp); + if (utils::strmatch(arg[2], "^v_")) { + xstr = utils::strdup(arg[2] + 2); + } else + xvalue = utils::numeric(FLERR, arg[2], false, lmp); + if (utils::strmatch(arg[3], "^v_")) { + ystr = utils::strdup(arg[3] + 2); + } else + yvalue = utils::numeric(FLERR, arg[3], false, lmp); - } else error->all(FLERR,"Unknown fix indent cone argument: {}", arg[1]); + } else + error->all(FLERR, "Unknown fix indent cone argument: {}", arg[1]); + + if (utils::strmatch(arg[4], "^v_")) { + rlostr = utils::strdup(arg[4] + 2); + } else + rlovalue = utils::numeric(FLERR, arg[4], false, lmp); + if (utils::strmatch(arg[5], "^v_")) { + rhistr = utils::strdup(arg[5] + 2); + } else + rhivalue = utils::numeric(FLERR, arg[5], false, lmp); + if (utils::strmatch(arg[6], "^v_")) { + lostr = utils::strdup(arg[6] + 2); + } else + lovalue = utils::numeric(FLERR, arg[6], false, lmp); + if (utils::strmatch(arg[7], "^v_")) { + histr = utils::strdup(arg[7] + 2); + } else + hivalue = utils::numeric(FLERR, arg[7], false, lmp); - if (utils::strmatch(arg[4],"^v_")) { - rlostr = utils::strdup(arg[4]+2); - } else rlovalue = utils::numeric(FLERR,arg[4],false,lmp); - if (utils::strmatch(arg[5],"^v_")) { - rhistr = utils::strdup(arg[5]+2); - } else rhivalue = utils::numeric(FLERR,arg[5],false,lmp); - if (utils::strmatch(arg[6],"^v_")) { - lostr = utils::strdup(arg[6]+2); - } else lovalue = utils::numeric(FLERR,arg[6],false,lmp); - if (utils::strmatch(arg[7],"^v_")) { - histr = utils::strdup(arg[7]+2); - } else hivalue = utils::numeric(FLERR,arg[7],false,lmp); - istyle = CONE; return 8; } // plane - - if (strcmp(arg[0],"plane") == 0) { + + if (strcmp(arg[0], "plane") == 0) { if (istyle != NONE) error->all(FLERR, "Fix indent requires a single geometry keyword"); if (4 > narg) utils::missing_cmd_args(FLERR, "fix indent plane", error); - if (strcmp(arg[1],"x") == 0) cdim = 0; - else if (strcmp(arg[1],"y") == 0) cdim = 1; - else if (strcmp(arg[1],"z") == 0) cdim = 2; - else error->all(FLERR,"Unknown fix indent plane argument: {}", arg[1]); - - if (utils::strmatch(arg[2],"^v_")) { - pstr = utils::strdup(arg[2]+2); - } else pvalue = utils::numeric(FLERR,arg[2],false,lmp); - - if (strcmp(arg[3],"lo") == 0) planeside = -1; - else if (strcmp(arg[3],"hi") == 0) planeside = 1; - else error->all(FLERR,"Unknown fix indent plane argument: {}", arg[3]); + if (strcmp(arg[1], "x") == 0) + cdim = 0; + else if (strcmp(arg[1], "y") == 0) + cdim = 1; + else if (strcmp(arg[1], "z") == 0) + cdim = 2; + else + error->all(FLERR, "Unknown fix indent plane argument: {}", arg[1]); + + if (utils::strmatch(arg[2], "^v_")) { + pstr = utils::strdup(arg[2] + 2); + } else + pvalue = utils::numeric(FLERR, arg[2], false, lmp); + + if (strcmp(arg[3], "lo") == 0) + planeside = -1; + else if (strcmp(arg[3], "hi") == 0) + planeside = 1; + else + error->all(FLERR, "Unknown fix indent plane argument: {}", arg[3]); istyle = PLANE; return 4; } - + // invalid istyle arg - - error->all(FLERR,"Unknown fix indent argument: {}", arg[0]); - + + error->all(FLERR, "Unknown fix indent argument: {}", arg[0]); + return 0; } @@ -640,21 +664,28 @@ void FixIndent::options(int narg, char **arg) int iarg = 0; while (iarg < narg) { - if (strcmp(arg[iarg],"units") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix indent units", error); - if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; - else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; - else error->all(FLERR,"Unknown fix indent units argument: {}", arg[iarg+1]); + if (strcmp(arg[iarg], "units") == 0) { + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix indent units", error); + if (strcmp(arg[iarg + 1], "box") == 0) + scaleflag = 0; + else if (strcmp(arg[iarg + 1], "lattice") == 0) + scaleflag = 1; + else + error->all(FLERR, "Unknown fix indent units argument: {}", arg[iarg + 1]); iarg += 2; - } else if (strcmp(arg[iarg],"side") == 0) { - if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix indent side", error); - if (strcmp(arg[iarg+1],"in") == 0) side = INSIDE; - else if (strcmp(arg[iarg+1],"out") == 0) side = OUTSIDE; - else error->all(FLERR,"Unknown fix indent side argument: {}", arg[iarg+1]); + } else if (strcmp(arg[iarg], "side") == 0) { + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix indent side", error); + if (strcmp(arg[iarg + 1], "in") == 0) + side = INSIDE; + else if (strcmp(arg[iarg + 1], "out") == 0) + side = OUTSIDE; + else + error->all(FLERR, "Unknown fix indent side argument: {}", arg[iarg + 1]); iarg += 2; - } else error->all(FLERR,"Unknown fix indent argument: {}", arg[iarg]); + } else + error->all(FLERR, "Unknown fix indent argument: {}", arg[iarg]); } } @@ -662,12 +693,12 @@ void FixIndent::options(int narg, char **arg) determines if a point is inside (true) or outside (false) of a cone ------------------------------------------------------------------------- */ -bool FixIndent::PointInsideCone(int dir, double *center, double lo, - double hi, double rlo, double rhi, double *x) +bool FixIndent::PointInsideCone(int dir, double *center, double lo, double hi, double rlo, + double rhi, double *x) { if ((x[dir] > hi) || (x[dir] < lo)) return false; - double del[3] {x[0] - center[0], x[1] - center[1], x[2] - center[2]}; + double del[3] = {x[0] - center[0], x[1] - center[1], x[2] - center[2]}; del[dir] = 0.0; double dist = sqrt(del[0] * del[0] + del[1] * del[1] + del[2] * del[2]); @@ -682,17 +713,14 @@ bool FixIndent::PointInsideCone(int dir, double *center, double lo, distance between an exterior point and a cone ------------------------------------------------------------------------- */ -void FixIndent::DistanceExteriorPoint(int dir, double *center, double lo, double hi, - double rlo, double rhi, - double &x, double &y, double &z) +void FixIndent::DistanceExteriorPoint(int dir, double *center, double lo, double hi, double rlo, + double rhi, double &x, double &y, double &z) { double xp[3], nearest[3], corner1[3], corner2[3]; + double point[3] = {x, y, z}; + double del[3] = {x - center[0], y - center[1], z - center[2]}; - double point[3] {x, y, z}; - - double del[3] {x - center[0], y - center[1], z - center[2]}; del[dir] = 0.0; - double r = sqrt(del[0] * del[0] + del[1] * del[1] + del[2] * del[2]); corner1[0] = center[0] + del[0] * rlo / r; @@ -705,10 +733,10 @@ void FixIndent::DistanceExteriorPoint(int dir, double *center, double lo, double corner2[2] = center[2] + del[2] * rhi / r; corner2[dir] = hi; - double corner3[3] {center[0], center[1], center[2]}; + double corner3[3] = {center[0], center[1], center[2]}; corner3[dir] = lo; - double corner4[3] {center[0], center[1], center[2]}; + double corner4[3] = {center[0], center[1], center[2]}; corner4[dir] = hi; // initialize distance to a big number @@ -741,18 +769,17 @@ void FixIndent::DistanceExteriorPoint(int dir, double *center, double lo, double distance between an interior point and a cone ------------------------------------------------------------------------- */ -void FixIndent::DistanceInteriorPoint(int dir, double *center, - double lo, double hi, double rlo, double rhi, double &x, - double &y, double &z) +void FixIndent::DistanceInteriorPoint(int dir, double *center, double lo, double hi, double rlo, + double rhi, double &x, double &y, double &z) { double r, dist_disk, dist_surf; double surflo[3], surfhi[3], xs[3]; - double initial_point[3] {x, y, z}; - double point[3] {0.0, 0.0, 0.0}; + double initial_point[3] = {x, y, z}; + double point[3] = {0.0, 0.0, 0.0}; // initial check with the two disks - if ( (initial_point[dir] - lo) < (hi - initial_point[dir]) ) { + if ((initial_point[dir] - lo) < (hi - initial_point[dir])) { dist_disk = (initial_point[dir] - lo) * (initial_point[dir] - lo); point[dir] = initial_point[dir] - lo; } else { @@ -762,7 +789,7 @@ void FixIndent::DistanceInteriorPoint(int dir, double *center, // check with the points in the conical surface - double del[3] {x - center[0], y - center[1], z - center[2]}; + double del[3] = {x - center[0], y - center[1], z - center[2]}; del[dir] = 0.0; r = sqrt(del[0] * del[0] + del[1] * del[1] + del[2] * del[2]); @@ -778,7 +805,7 @@ void FixIndent::DistanceInteriorPoint(int dir, double *center, point_on_line_segment(surflo, surfhi, initial_point, xs); - double dx[3] {initial_point[0] - xs[0], initial_point[1] - xs[1], initial_point[2] - xs[2]}; + double dx[3] = {initial_point[0] - xs[0], initial_point[1] - xs[1], initial_point[2] - xs[2]}; dist_surf = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; if (dist_surf < dist_disk) { x = dx[0]; From 5b16cf97736efccff080be19768f688f7d237362 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 10 Mar 2024 16:19:22 -0400 Subject: [PATCH 15/18] use std::move() to avoid extra copy of temporaries --- src/ELECTRODE/fix_electrode_conp.cpp | 8 ++++---- src/GRANULAR/granular_model.cpp | 5 +++-- src/reader_native.cpp | 2 +- 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/ELECTRODE/fix_electrode_conp.cpp b/src/ELECTRODE/fix_electrode_conp.cpp index 94c085de5c..f3c17062d2 100644 --- a/src/ELECTRODE/fix_electrode_conp.cpp +++ b/src/ELECTRODE/fix_electrode_conp.cpp @@ -915,20 +915,20 @@ void FixElectrodeConp::update_charges() a = ele_ele_interaction(q_local); r = add_nlocalele(b, a); } else { - r = add_nlocalele(r, scale_vector(alpha, y)); + r = add_nlocalele(r, scale_vector(alpha, std::move(y))); } auto p = constraint_projection(r); double dot_new = dot_nlocalele(r, p); - d = add_nlocalele(p, scale_vector(dot_new / dot_old, d)); + d = add_nlocalele(std::move(p), scale_vector(dot_new / dot_old, d)); delta = dot_nlocalele(r, d); dot_old = dot_new; } - recompute_potential(b, q_local); + recompute_potential(std::move(b), q_local); if (delta > cg_threshold && comm->me == 0) error->warning(FLERR, "CG threshold not reached"); } else { error->all(FLERR, "This algorithm is not implemented, yet"); } - set_charges(q_local); + set_charges(std::move(q_local)); update_time += MPI_Wtime() - start; } diff --git a/src/GRANULAR/granular_model.cpp b/src/GRANULAR/granular_model.cpp index 6de147b34a..14431f41b4 100644 --- a/src/GRANULAR/granular_model.cpp +++ b/src/GRANULAR/granular_model.cpp @@ -32,6 +32,7 @@ #include #include +#include using namespace LAMMPS_NS; using namespace Granular_NS; @@ -333,11 +334,11 @@ void GranularModel::read_restart(FILE *fp) utils::sfread(FLERR, &num_char, sizeof(int), 1, fp, nullptr, error); MPI_Bcast(&num_char, 1, MPI_INT, 0, world); - std::string model_name (num_char, ' '); + std::string model_name(num_char, ' '); if (comm->me == 0) utils::sfread(FLERR, const_cast(model_name.data()), sizeof(char),num_char, fp, nullptr, error); MPI_Bcast(const_cast(model_name.data()), num_char, MPI_CHAR, 0, world); - construct_sub_model(model_name, (SubModelType) i); + construct_sub_model(std::move(model_name), (SubModelType) i); if (comm->me == 0) utils::sfread(FLERR, &num_coeff, sizeof(int), 1, fp, nullptr, error); diff --git a/src/reader_native.cpp b/src/reader_native.cpp index ae59ca6805..4dac65e3cb 100644 --- a/src/reader_native.cpp +++ b/src/reader_native.cpp @@ -289,7 +289,7 @@ bigint ReaderNative::read_header(double box[3][3], int &boxinfo, int &triclinic, labelline = line + strlen("ITEM: ATOMS "); } - Tokenizer tokens(labelline); + Tokenizer tokens(std::move(labelline)); std::map labels; nwords = 0; From c3c61a368d1471a836561523daa657ee1e4cfef2 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 11 Mar 2024 15:20:42 -0600 Subject: [PATCH 16/18] Bugfix: both pair hybrid and hybrid/overlay cannot fuse force zeroing --- src/KOKKOS/pair_kokkos.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/KOKKOS/pair_kokkos.h b/src/KOKKOS/pair_kokkos.h index 87324b49b9..54502f290c 100644 --- a/src/KOKKOS/pair_kokkos.h +++ b/src/KOKKOS/pair_kokkos.h @@ -1020,7 +1020,7 @@ template EV_FLOAT pair_compute (PairStyle* fpair, NeighListKokkos* list) { EV_FLOAT ev; if (fpair->neighflag == FULL) { - if (utils::strmatch(fpair->lmp->force->pair_style,"^hybrid/overlay")) { + if (utils::strmatch(fpair->lmp->force->pair_style,"^hybrid")) { fpair->fuse_force_clear_flag = 0; ev = pair_compute_neighlist (fpair,list); } else { From d0d4cf9ad02dcbcfbe5e7684824d6604c67f495a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 11 Mar 2024 22:40:21 -0400 Subject: [PATCH 17/18] add DOIs for recent stable releases --- doc/src/Intro_citing.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/src/Intro_citing.rst b/doc/src/Intro_citing.rst index 69ccab6162..8b906bb725 100644 --- a/doc/src/Intro_citing.rst +++ b/doc/src/Intro_citing.rst @@ -47,6 +47,8 @@ In addition there are DOIs generated for individual stable releases: - 3 March 2020 version: `DOI:10.5281/zenodo.3726417 `_ - 29 October 2020 version: `DOI:10.5281/zenodo.4157471 `_ - 29 September 2021 version: `DOI:10.5281/zenodo.6386596 `_ +- 23 June 2022 version: `DOI:10.5281/zenodo.10806836 `_ +- 2 August 2023 version: `DOI:10.5281/zenodo.10806852 `_ Home page ^^^^^^^^^ From aa7e2da33dac69fc5b65d5d88d130733ca5f5746 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Wed, 13 Mar 2024 15:04:19 -0600 Subject: [PATCH 18/18] Misc small patches --- doc/src/fix_deform_pressure.rst | 2 +- src/EXTRA-FIX/fix_deform_pressure.cpp | 4 +--- src/EXTRA-FIX/fix_deform_pressure.h | 1 - src/neighbor.cpp | 3 +++ 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/src/fix_deform_pressure.rst b/doc/src/fix_deform_pressure.rst index f85ad37238..c814aa892f 100644 --- a/doc/src/fix_deform_pressure.rst +++ b/doc/src/fix_deform_pressure.rst @@ -91,7 +91,7 @@ corresponding component of the pressure tensor. This option attempts to maintain a specified target pressure using a linear controller where the box length :math:`L` evolves according to the equation -.. parsed-literal:: +.. math:: \frac{d L(t)}{dt} = L(t) k (P_t - P) diff --git a/src/EXTRA-FIX/fix_deform_pressure.cpp b/src/EXTRA-FIX/fix_deform_pressure.cpp index ffa3f11d92..51ea75cfed 100644 --- a/src/EXTRA-FIX/fix_deform_pressure.cpp +++ b/src/EXTRA-FIX/fix_deform_pressure.cpp @@ -767,8 +767,7 @@ void FixDeformPressure::apply_box() if (fabs(v_rate) > max_h_rate) v_rate = max_h_rate * v_rate / fabs(v_rate); - set_extra[6].cumulative_strain += update->dt * v_rate; - scale = (1.0 + set_extra[6].cumulative_strain); + scale = (1.0 + update->dt * v_rate); for (i = 0; i < 3; i++) { shift = 0.5 * (set[i].hi_target - set[i].lo_target) * scale; set[i].lo_target = 0.5 * (set[i].lo_start + set[i].hi_start) - shift; @@ -843,7 +842,6 @@ void FixDeformPressure::restart(char *buf) set_extra[i].saved = set_extra_restart[i].saved; set_extra[i].prior_rate = set_extra_restart[i].prior_rate; set_extra[i].prior_pressure = set_extra_restart[i].prior_pressure; - set_extra[i].cumulative_strain = set_extra_restart[i].cumulative_strain; } } diff --git a/src/EXTRA-FIX/fix_deform_pressure.h b/src/EXTRA-FIX/fix_deform_pressure.h index 5a0d844bad..10af1e5ba3 100644 --- a/src/EXTRA-FIX/fix_deform_pressure.h +++ b/src/EXTRA-FIX/fix_deform_pressure.h @@ -51,7 +51,6 @@ class FixDeformPressure : public FixDeform { struct SetExtra { double ptarget, pgain; double prior_pressure, prior_rate; - double cumulative_strain; int saved; char *pstr; int pvar, pvar_flag; diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 662494ea7b..c5cbe0e885 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -427,6 +427,9 @@ void Neighbor::init() } } } else { + if (!force->pair) + error->all(FLERR, "Cannot use collection/interval command without defining a pairstyle"); + if (force->pair->finitecutflag) { finite_cut_flag = 1; // If cutoffs depend on finite atom sizes, use radii of intervals to find cutoffs