From 25f5e74e9a7b240fe123f59c7749714d7304af46 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 15 Jun 2023 15:51:59 -0600 Subject: [PATCH 01/18] updating args and D2min calculation in nonaffine fix --- doc/src/compute.rst | 1 + doc/src/compute_contact_atom.rst | 10 +- doc/src/compute_rattlers.rst | 74 +++ doc/src/fix.rst | 1 + doc/src/fix_nonaffine_displacement.rst | 140 +++++ src/GRANULAR/compute_contact_atom.cpp | 52 +- src/GRANULAR/compute_contact_atom.h | 4 + src/compute_rattlers.cpp | 308 ++++++++++ src/compute_rattlers.h | 52 ++ src/fix_nonaffine_displacement.cpp | 740 +++++++++++++++++++++++++ src/fix_nonaffine_displacement.h | 70 +++ 11 files changed, 1430 insertions(+), 22 deletions(-) create mode 100644 doc/src/compute_rattlers.rst create mode 100644 doc/src/fix_nonaffine_displacement.rst create mode 100755 src/compute_rattlers.cpp create mode 100755 src/compute_rattlers.h create mode 100755 src/fix_nonaffine_displacement.cpp create mode 100755 src/fix_nonaffine_displacement.h diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 6fdedbbb95..0f6bfa75d9 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -254,6 +254,7 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`property/chunk ` - extract various per-chunk attributes * :doc:`property/local ` - convert local attributes to localvectors/arrays * :doc:`ptm/atom ` - determines the local lattice structure based on the Polyhedral Template Matching method +* :doc:`rattlers ` - identify undercoordinated rattler atoms * :doc:`rdf ` - radial distribution function g(r) histogram of group of atoms * :doc:`reduce ` - combine per-atom quantities into a single global value * :doc:`reduce/chunk ` - reduce per-atom quantities within each chunk diff --git a/doc/src/compute_contact_atom.rst b/doc/src/compute_contact_atom.rst index 20dcbfae29..c26f2b4252 100644 --- a/doc/src/compute_contact_atom.rst +++ b/doc/src/compute_contact_atom.rst @@ -8,10 +8,11 @@ Syntax .. parsed-literal:: - compute ID group-ID contact/atom + compute ID group-ID contact/atom group2-ID * ID, group-ID are documented in :doc:`compute ` command * contact/atom = style name of this compute command +* group2-ID = optional argument select group-ID to restrict which atoms to consider for contacts (see below) Examples """""""" @@ -19,6 +20,7 @@ Examples .. code-block:: LAMMPS compute 1 all contact/atom + compute 1 all contact/atom mygroup Description """"""""""" @@ -34,6 +36,9 @@ sum of the radii of the two particles. The value of the contact number will be 0.0 for atoms not in the specified compute group. +The optional *group2-ID* argument allows to specify from which group atoms +contribute to the coordination number. Default setting is group 'all'. + Output info """"""""""" @@ -63,4 +68,7 @@ Related commands Default """"""" +*group2-ID* = all + + none diff --git a/doc/src/compute_rattlers.rst b/doc/src/compute_rattlers.rst new file mode 100644 index 0000000000..c26f2b4252 --- /dev/null +++ b/doc/src/compute_rattlers.rst @@ -0,0 +1,74 @@ +.. index:: compute contact/atom + +compute contact/atom command +============================ + +Syntax +"""""" + +.. parsed-literal:: + + compute ID group-ID contact/atom group2-ID + +* ID, group-ID are documented in :doc:`compute ` command +* contact/atom = style name of this compute command +* group2-ID = optional argument select group-ID to restrict which atoms to consider for contacts (see below) + +Examples +"""""""" + +.. code-block:: LAMMPS + + compute 1 all contact/atom + compute 1 all contact/atom mygroup + +Description +""""""""""" + +Define a computation that calculates the number of contacts +for each atom in a group. + +The contact number is defined for finite-size spherical particles as +the number of neighbor atoms which overlap the central particle, +meaning that their distance of separation is less than or equal to the +sum of the radii of the two particles. + +The value of the contact number will be 0.0 for atoms not in the +specified compute group. + +The optional *group2-ID* argument allows to specify from which group atoms +contribute to the coordination number. Default setting is group 'all'. + +Output info +""""""""""" + +This compute calculates a per-atom vector, whose values can be +accessed by any command that uses per-atom values from a compute as +input. See the :doc:`Howto output ` page for an +overview of LAMMPS output options. + +The per-atom vector values will be a number >= 0.0, as explained +above. + +Restrictions +"""""""""""" + +This compute is part of the GRANULAR package. It is only enabled if +LAMMPS was built with that package. See the +:doc:`Build package ` page for more info. + +This compute requires that atoms store a radius as defined by the +:doc:`atom_style sphere ` command. + +Related commands +"""""""""""""""" + +:doc:`compute coord/atom ` + +Default +""""""" + +*group2-ID* = all + + +none diff --git a/doc/src/fix.rst b/doc/src/fix.rst index b0ec47fbe6..b07aa2770e 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -259,6 +259,7 @@ accelerated styles exist. * :doc:`mvv/tdpd ` - constant temperature DPD using the modified velocity-Verlet algorithm * :doc:`neb ` - nudged elastic band (NEB) spring forces * :doc:`neb/spin ` - nudged elastic band (NEB) spring forces for spins +* :doc:`nonaffine/displacement ` - calculate nonaffined displacement of atoms * :doc:`nph ` - constant NPH time integration via Nose/Hoover * :doc:`nph/asphere ` - NPH for aspherical particles * :doc:`nph/body ` - NPH for body particles diff --git a/doc/src/fix_nonaffine_displacement.rst b/doc/src/fix_nonaffine_displacement.rst new file mode 100644 index 0000000000..7b0f7669ed --- /dev/null +++ b/doc/src/fix_nonaffine_displacement.rst @@ -0,0 +1,140 @@ +.. index:: fix gravity +.. index:: fix gravity/omp +.. index:: fix gravity/kk + +fix gravity command +=================== + +Accelerator Variants: *gravity/omp*, *gravity/kk* + +Syntax +"""""" + +.. parsed-literal:: + + fix ID group gravity magnitude style args + +* ID, group are documented in :doc:`fix ` command +* gravity = style name of this fix command +* magnitude = size of acceleration (force/mass units) +* magnitude can be a variable (see below) +* style = *chute* or *spherical* or *gradient* or *vector* + + .. parsed-literal:: + + *chute* args = angle + angle = angle in +x away from -z or -y axis in 3d/2d (in degrees) + angle can be a variable (see below) + *spherical* args = phi theta + phi = azimuthal angle from +x axis (in degrees) + theta = angle from +z or +y axis in 3d/2d (in degrees) + phi or theta can be a variable (see below) + *vector* args = x y z + x y z = vector direction to apply the acceleration + x or y or z can be a variable (see below) + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix 1 all gravity 1.0 chute 24.0 + fix 1 all gravity v_increase chute 24.0 + fix 1 all gravity 1.0 spherical 0.0 -180.0 + fix 1 all gravity 10.0 spherical v_phi v_theta + fix 1 all gravity 100.0 vector 1 1 0 + +Description +""""""""""" + +Impose an additional acceleration on each particle in the group. This +fix is typically used with granular systems to include a "gravity" +term acting on the macroscopic particles. More generally, it can +represent any kind of driving field, e.g. a pressure gradient inducing +a Poiseuille flow in a fluid. Note that this fix operates differently +than the :doc:`fix addforce ` command. The addforce fix +adds the same force to each atom, independent of its mass. This +command imparts the same acceleration to each atom (force/mass). + +The *magnitude* of the acceleration is specified in force/mass units. +For granular systems (LJ units) this is typically 1.0. See the +:doc:`units ` command for details. + +Style *chute* is typically used for simulations of chute flow where +the specified *angle* is the chute angle, with flow occurring in the +x +direction. For 3d systems, the tilt is away from the z axis; for 2d +systems, the tilt is away from the y axis. + +Style *spherical* allows an arbitrary 3d direction to be specified for +the acceleration vector. *Phi* and *theta* are defined in the usual +spherical coordinates. Thus for acceleration acting in the -z +direction, *theta* would be 180.0 (or -180.0). *Theta* = 90.0 and +*phi* = -90.0 would mean acceleration acts in the -y direction. For +2d systems, *phi* is ignored and *theta* is an angle in the xy plane +where *theta* = 0.0 is the y-axis. + +Style *vector* imposes an acceleration in the vector direction given +by (x,y,z). Only the direction of the vector is important; it's +length is ignored. For 2d systems, the *z* component is ignored. + +Any of the quantities *magnitude*, *angle*, *phi*, *theta*, *x*, *y*, +*z* which define the gravitational magnitude and direction, can be +specified as an equal-style :doc:`variable `. If the value 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 determine the quantity. You should +insure that the variable calculates a result in the appropriate units, +e.g. force/mass or degrees. + +Equal-style variables can specify formulas with various mathematical +functions, and include :doc:`thermo_style ` command +keywords for the simulation box parameters and timestep and elapsed +time. Thus it is easy to specify a time-dependent gravitational +field. + +---------- + +.. include:: accel_styles.rst + +---------- + +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +No information about this fix is written to :doc:`binary restart files `. + +The :doc:`fix_modify ` *energy* option is supported by +this fix to add the gravitational potential energy of the system to +the global potential energy of the system as part of +:doc:`thermodynamic output `. The default setting for +this fix is :doc:`fix_modify energy no `. + +The :doc:`fix_modify ` *respa* option is supported by this +fix. This allows to set at which level of the :doc:`r-RESPA +` integrator the fix is adding its forces. Default is the +outermost level. + +This fix computes a global scalar which can be accessed by various +:doc:`output commands `. This scalar is the +gravitational potential energy of the particles in the defined field, +namely mass \* (g dot x) for each particles, where x and mass are the +particles position and mass, and g is the gravitational field. The +scalar value calculated by this fix is "extensive". + +No parameter of this fix can be used with the *start/stop* keywords of +the :doc:`run ` command. This fix is not invoked during +:doc:`energy minimization `. + +Restrictions +"""""""""""" + none + +Related commands +"""""""""""""""" + +:doc:`atom_style sphere `, :doc:`fix addforce ` + +Default +""""""" + +none diff --git a/src/GRANULAR/compute_contact_atom.cpp b/src/GRANULAR/compute_contact_atom.cpp index 5a07d14eb8..26018f7bfb 100644 --- a/src/GRANULAR/compute_contact_atom.cpp +++ b/src/GRANULAR/compute_contact_atom.cpp @@ -18,6 +18,7 @@ #include "comm.h" #include "error.h" #include "force.h" +#include "group.h" #include "memory.h" #include "modify.h" #include "neigh_list.h" @@ -34,7 +35,16 @@ ComputeContactAtom::ComputeContactAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), contact(nullptr) { - if (narg != 3) error->all(FLERR,"Illegal compute contact/atom command"); + if (narg != 3 && narg != 4) error->all(FLERR,"Illegal compute contact/atom command"); + + jgroup = group->find("all"); + jgroupbit = group->bitmask[jgroup]; + if (narg == 4) { + group2 = utils::strdup(arg[3]); + jgroup = group->find(group2); + if (jgroup == -1) error->all(FLERR, "Compute contact/atom group2 ID does not exist"); + jgroupbit = group->bitmask[jgroup]; + } peratom_flag = 1; size_peratom_cols = 0; @@ -120,28 +130,28 @@ void ComputeContactAtom::compute_peratom() for (ii = 0; ii < inum; ii++) { i = ilist[ii]; - if (mask[i] & groupbit) { - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - radi = radius[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + radi = radius[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - radsum = radi + radius[j]; - radsumsq = radsum*radsum; - if (rsq <= radsumsq) { - contact[i] += 1.0; - contact[j] += 1.0; - } + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + radsum = radi + radius[j]; + radsumsq = radsum * radsum; + if (rsq <= radsumsq) { + // Only tally for atoms in compute group (groupbit) if neighbor is in group2 (jgroupbit) + if ((mask[i] & groupbit) && (mask[j] & jgroupbit)) contact[i] += 1.0; + if ((mask[j] & groupbit) && (mask[i] & jgroupbit)) contact[j] += 1.0; } } } diff --git a/src/GRANULAR/compute_contact_atom.h b/src/GRANULAR/compute_contact_atom.h index f01815e666..5262c1b29b 100644 --- a/src/GRANULAR/compute_contact_atom.h +++ b/src/GRANULAR/compute_contact_atom.h @@ -37,6 +37,10 @@ class ComputeContactAtom : public Compute { private: int nmax; + + char *group2; + int jgroup, jgroupbit; + class NeighList *list; double *contact; }; diff --git a/src/compute_rattlers.cpp b/src/compute_rattlers.cpp new file mode 100755 index 0000000000..5a6d2b022d --- /dev/null +++ b/src/compute_rattlers.cpp @@ -0,0 +1,308 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "compute_rattlers.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "neighbor.h" +#include "pair.h" +#include "update.h" +#include "utils.h" + +#include +#include + +using namespace LAMMPS_NS; + +enum { TYPE, RADIUS }; + +/* ---------------------------------------------------------------------- */ + +ComputeRattlers::ComputeRattlers(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), ncontacts(nullptr), rattler(nullptr) +{ + if (narg != 6) error->all(FLERR, "Illegal compute fabric command"); + + if (strcmp(arg[3], "type") == 0) + cutstyle = TYPE; + else if (strcmp(arg[3], "radius") == 0) + cutstyle = RADIUS; + else + error->all(FLERR, "Illegal compute fabric command"); + + if (cutstyle == RADIUS && !atom->radius_flag) + error->all(FLERR, "Compute fabric radius style requires atom attribute radius"); + + ncontacts_rattler = utils::inumeric(FLERR, arg[4], false, lmp); + max_tries = utils::inumeric(FLERR, arg[5], false, lmp); + + nmax = 0; + invoked_peratom = -1; + + scalar_flag = 1; + extscalar = 1; + peratom_flag = 1; + size_peratom_cols = 0; + comm_forward = 1; + comm_reverse = 1; +} + +/* ---------------------------------------------------------------------- */ + +ComputeRattlers::~ComputeRattlers() +{ + memory->destroy(ncontacts); + memory->destroy(rattler); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRattlers::init() +{ + if (force->pair == nullptr) error->all(FLERR, "No pair style is defined for compute rattlers"); + + // Cannot calculate distance from radii for JKR/DMT + if (force->pair->beyond_contact) + error->all(FLERR, "Compute rattlers does not currently support pair styles that extend beyond contact"); + + // need an occasional half neighbor list + // set size to same value as request made by force->pair + // this should enable it to always be a copy list (e.g. for granular pstyle) + + auto pairrequest = neighbor->find_request(force->pair); + if (pairrequest && pairrequest->get_size()) + neighbor->add_request(this, NeighConst::REQ_SIZE | NeighConst::REQ_OCCASIONAL); + else + neighbor->add_request(this, NeighConst::REQ_OCCASIONAL); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRattlers::init_list(int /*id*/, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRattlers::compute_peratom() +{ + if (invoked_peratom == update->ntimestep) return; + invoked_peratom = update->ntimestep; + + int i, j, ii, jj, inum, jnum, itype, jtype, tmp_flag; + tagint itag, jtag; + double xtmp, ytmp, ztmp, delx, dely, delz; + double r, rinv, rsq, radsum; + + if (nmax < atom->nmax) { + nmax = atom->nmax; + memory->destroy(ncontacts); + memory->destroy(rattler); + memory->create(ncontacts, nmax, "rattlers:ncontacts"); + memory->create(rattler, nmax, "rattlers:rattler"); + vector_atom = rattler; + } + + for (i = 0; i < nmax; i++) rattler[i] = 0; + + int *ilist, *jlist, *numneigh, **firstneigh; + + double **x = atom->x; + double *radius = atom->radius; + tagint *tag = atom->tag; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + // invoke half neighbor list (will copy or build if necessary) + neighbor->build_one(list); + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + Pair *pair = force->pair; + double **cutsq = force->pair->cutsq; + + int change_flag = 1; + int ntry = 0; + while (ntry < max_tries) { + // Quit when answer stops evolving + if (change_flag == 0) break; + change_flag = 0; + + for (i = 0; i < nmax; i++) ncontacts[i] = 0; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (!(mask[i] & groupbit)) continue; + if (rattler[i] == 1) continue; + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itag = tag[i]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + if (!(mask[j] & groupbit)) continue; + if (rattler[j] == 1) continue; + + // itag = jtag is possible for long cutoffs that include images of self + + if (newton_pair == 0 && j >= nlocal) { + jtag = tag[j]; + if (itag > jtag) { + if ((itag + jtag) % 2 == 0) continue; + } else if (itag < jtag) { + if ((itag + jtag) % 2 == 1) continue; + } else { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + } + } + + jtype = type[j]; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + + if (cutstyle == TYPE) { + if (rsq >= cutsq[itype][jtype]) continue; + } else { + radsum = radius[i] + radius[j]; + if (rsq >= radsum * radsum) continue; + } + ncontacts[i] += 1; + if (newton_pair || j < nlocal) + ncontacts[j] += 1; + } + } + + // add contributions from ghosts + if (force->newton_pair) comm->reverse_comm(this); + + // Set flags for rattlers + for (i = 0; i < atom->nlocal; i++) { + if (ncontacts[i] < ncontacts_rattler && rattler[i] == 0) { + rattler[i] = 1; + change_flag = 1; + } + } + + comm->forward_comm(this); + + MPI_Allreduce(&change_flag, &tmp_flag, 1, MPI_INT, MPI_MAX, world); + change_flag = tmp_flag; + + ntry += 1; + } + + if (change_flag == 1) + error->warning(FLERR, "Rattler calculation failed to converge within max tries"); +} + +/* ---------------------------------------------------------------------- */ + +double ComputeRattlers::compute_scalar() { + if (invoked_peratom != update->ntimestep) + compute_peratom(); + + invoked_scalar = update->ntimestep; + + double total_rattlers = 0; + for (int i = 0; i < atom->nlocal; i++) { + if (rattler[i] == 1) { + total_rattlers += 1; + } + } + + //Total across processors + MPI_Allreduce(&total_rattlers, &scalar, 1, MPI_DOUBLE, MPI_SUM, world); + return scalar; +} + +/* ---------------------------------------------------------------------- */ + +int ComputeRattlers::pack_reverse_comm(int n, int first, double *buf) +{ + int i, m, last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + buf[m++] = ubuf(ncontacts[i]).d; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRattlers::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i, j, m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + ncontacts[j] += (int) ubuf(buf[m++]).i; + } +} + +/* ---------------------------------------------------------------------- */ + +int ComputeRattlers::pack_forward_comm(int n, int *list, double *buf, + int /*pbc_flag*/, int * /*pbc*/) +{ + int i, j, m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = rattler[j]; + } + + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRattlers::unpack_forward_comm(int n, int first, double *buf) +{ + int i, m, last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + rattler[i] = buf[m++]; + } +} diff --git a/src/compute_rattlers.h b/src/compute_rattlers.h new file mode 100755 index 0000000000..08f8d7e71d --- /dev/null +++ b/src/compute_rattlers.h @@ -0,0 +1,52 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS +// clang-format off +ComputeStyle(rattlers,ComputeRattlers); +// clang-format on +#else + +#ifndef LMP_COMPUTE_RATTLERS_H +#define LMP_COMPUTE_RATTLERS_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeRattlers : public Compute { + public: + ComputeRattlers(class LAMMPS *, int, char **); + ~ComputeRattlers() override; + void init() override; + void init_list(int, class NeighList *) override; + void compute_peratom() override; + double compute_scalar() override; + int pack_forward_comm(int, int *, double *, int, int *) override; + void unpack_forward_comm(int, int, double *) override; + int pack_reverse_comm(int, int, double *) override; + void unpack_reverse_comm(int, int *, double *) override; + + private: + int pstyle, cutstyle; + int ncontacts_rattler, max_tries, nmax, invoked_peratom; + int *ncontacts; + double *rattler; + class NeighList *list; + +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/fix_nonaffine_displacement.cpp b/src/fix_nonaffine_displacement.cpp new file mode 100755 index 0000000000..15c678f83d --- /dev/null +++ b/src/fix_nonaffine_displacement.cpp @@ -0,0 +1,740 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "fix_nonaffine_displacement.h" + +#include "atom.h" +#include "citeme.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "group.h" +#include "math_extra.h" +#include "memory.h" +#include "modify.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "neighbor.h" +#include "pair.h" +#include "update.h" + +#include + +using namespace LAMMPS_NS; +using namespace FixConst; +using namespace MathExtra; + +enum { TYPE, RADIUS, CUSTOM }; +enum { INTEGRATED, D2MIN }; +enum { FIXED, OFFSET, UPDATE }; + +static const char cite_nonaffine_d2min[] = + "@article{PhysRevE.57.7192,\n" + " title = {Dynamics of viscoplastic deformation in amorphous solids},\n" + " author = {Falk, M. L. and Langer, J. S.},\n" + " journal = {Phys. Rev. E},\n" + " volume = {57},\n" + " issue = {6},\n" + " pages = {7192--7205},\n" + " numpages = {0},\n" + " year = {1998},\n" + " month = {Jun},\n" + " publisher = {American Physical Society},\n" + " doi = {10.1103/PhysRevE.57.7192},\n" + "url = {https://link.aps.org/doi/10.1103/PhysRevE.57.7192}\n" + "}\n\n"; + +/* ---------------------------------------------------------------------- */ + +FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), new_fix_id(nullptr), X(nullptr), Y(nullptr), F(nullptr), norm (nullptr) +{ + if (narg < 4) error->all(FLERR,"Illegal fix nonaffine/displacement command"); + + int iarg = 3; + if (strcmp(arg[iarg], "integrated") == 0) { + nad_style = INTEGRATED; + nevery = 1; + iarg += 1; + } else if (strcmp(arg[iarg], "d2min") == 0) { + if (iarg + 2 > narg) error->all(FLERR,"Illegal fix nonaffine/displacement command"); + + nad_style = D2MIN; + nevery = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + if (nevery <= 0) error->all(FLERR,"Illegal nevery value {} in fix nonaffine/displacement", nevery); + + if (strcmp(arg[iarg + 2], "type") == 0) { + cut_style = TYPE; + } else if (strcmp(arg[iarg + 2], "radius") == 0) { + cut_style = RADIUS; + } else if (strcmp(arg[iarg + 2], "custom") == 0) { + if (iarg + 3 > narg) error->all(FLERR,"Illegal fix nonaffine/displacement command"); + cut_style = CUSTOM; + cutoff_custom = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + cutsq_custom = cutoff_custom * cutoff_custom; + if (cutoff_custom <= 0) + error->all(FLERR, "Illegal custom cutoff length {}", arg[iarg + 3]); + iarg += 1; + } else error->all(FLERR,"Illegal cutoff style {} in fix nonaffine/displacement", arg[iarg + 2]); + iarg += 3; + } + + if (iarg + 2 > narg) error->all(FLERR,"Illegal fix nonaffine/displacement command"); + if (strcmp(arg[iarg], "fixed") == 0) { + reference_style = FIXED; + reference_timestep = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + if (update_timestep < 0) + error->all(FLERR, "Illegal reference timestep {} in fix nonaffine/displacement", arg[iarg + 1]); + } else if (strcmp(arg[iarg], "update") == 0) { + reference_style = UPDATE; + update_timestep = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + if (update_timestep < 0) + error->all(FLERR, "Illegal update timestep {} in fix nonaffine/displacement", arg[iarg + 1]); + } else if (strcmp(arg[iarg], "offset") == 0) { + reference_style = OFFSET; + offset_timestep = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + if (offset_timestep < 0) + error->all(FLERR, "Illegal offset timestep {} in fix nonaffine/displacement", arg[iarg + 1]); + } else error->all(FLERR,"Illegal reference style {} in fix nonaffine/displacement", arg[iarg]); + + if (cut_style == RADIUS && (!atom->radius_flag)) + error->all(FLERR, "Fix nonaffine/displacement radius style requires atom attribute radius"); + + peratom_flag = 1; + peratom_freq = nevery; + nmax = -1; + reference_saved = 0; + restart_global = 1; + + size_peratom_cols = 3; + comm_reverse = 0; + comm_forward = 0; + if (nad_style == D2MIN) { + comm_reverse = 18; + comm_forward = 9; + } + + if (nad_style == D2MIN && lmp->citeme) lmp->citeme->add(cite_nonaffine_d2min); +} + +/* ---------------------------------------------------------------------- */ + +FixNonaffineDisplacement::~FixNonaffineDisplacement() +{ + if (new_fix_id && modify->nfix) modify->delete_fix(new_fix_id); + delete[] new_fix_id; + + if (nad_style == D2MIN) { + memory->destroy(X); + memory->destroy(Y); + memory->destroy(F); + memory->destroy(norm); + memory->destroy(array_atom); + } +} + +/* ---------------------------------------------------------------------- */ + +int FixNonaffineDisplacement::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::post_constructor() +{ + // Create persistent peratom storage for either an integrated velocity or reference position + // Ghost atoms need reference coordinates for D2min + std::string ghost_status = "no"; + if (nad_style == D2MIN) ghost_status = "yes"; + + new_fix_id = utils::strdup(id + std::string("_FIX_PA")); + modify->add_fix(fmt::format("{} {} property/atom d2_nad 3 ghost {}", new_fix_id, group->names[igroup], ghost_status)); + int tmp1, tmp2; + nad_index = atom->find_custom("nad",tmp1,tmp2); + + if (nad_style == INTEGRATED) { + double **nad = atom->darray[nad_index]; + array_atom = nad; + } + + if (nad_style == D2MIN) { + grow_arrays(atom->nmax); + for (int i = 0; i < atom->nlocal; i++) + for (int j = 0; j < 3; j++) array_atom[i][j] = 0.0; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::init() +{ + dtv = update->dt; + + if (!reference_saved && (reference_style == FIXED) && (update->ntimestep > reference_timestep)) + error->all(FLERR, "Initial timestep exceeds that of the reference state in fix nonaffine/displacement"); + + if (nad_style == D2MIN) { + if (!force->pair && (cut_style == TYPE)) + error->all(FLERR,"Fix nonaffine/displacement D2Min option requires a pair style be defined " + "or cutoff specified"); + + if (cut_style == CUSTOM) { + double skin = neighbor->skin; + mycutneigh = cutoff_custom + skin; + + double cutghost; // as computed by Neighbor and Comm + if (force->pair) + cutghost = MAX(force->pair->cutforce + skin, comm->cutghostuser); + else + cutghost = comm->cutghostuser; + + if (mycutneigh > cutghost) + error->all(FLERR,"Fix nonaffine/displacement D2Min option cutoff exceeds ghost atom range - use comm_modify cutoff command"); + } + + // need an occasional half neighbor list + + if (cut_style == RADIUS) { + auto req = neighbor->add_request(this, NeighConst::REQ_SIZE | NeighConst::REQ_OCCASIONAL); + } else { + auto req = neighbor->add_request(this, NeighConst::REQ_OCCASIONAL); + if (cut_style == CUSTOM) req->set_cutoff(mycutneigh); + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::init_list(int /*id*/, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::setup(int vflag) +{ + post_force(0); // Save state if needed before starting the 1st timestep +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::post_force(int /*vflag*/) +{ + if (reference_style == FIXED) + if (update->ntimestep == reference_timestep) + save_reference_state(); + + if (reference_style == UPDATE) + if ((update->ntimestep % update_timestep) == 0) + save_reference_state(); + + if (reference_style == OFFSET) + if (((update->ntimestep - offset_timestep) % nevery) == 0) + save_reference_state(); + + if (update->setupflag) return; + + if (reference_style == FIXED && (update->ntimestep < reference_timestep)) return; + + if (nad_style == INTEGRATED) { + integrate_velocity(); + } else { + if ((update->ntimestep % nevery) == 0) calculate_D2Min(); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::write_restart(FILE *fp) +{ + if (comm->me == 0) { + int size = sizeof(int); + fwrite(&size, sizeof(int), 1, fp); + fwrite(&reference_saved, sizeof(int), 1, fp); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::restart(char *buf) +{ + reference_saved = (int) ubuf(buf[0]).i; +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::integrate_velocity() +{ + int i,n; + dtv = update->dt; + + double **v = atom->v; + + int *mask = atom->mask; + int nlocal = atom->nlocal; + double **nad = atom->darray[nad_index]; + + for (int m = 0; m < 3; m++) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + nad[i][m] += dtv * v[i][m]; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::save_reference_state() +{ + int i, n; + double **x = atom->x; + + int *mask = atom->mask; + int nlocal = atom->nlocal; + double **nad = atom->darray[nad_index]; + + if (nad_style == D2MIN) { + for (int m = 0; m < 3; m++) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) nad[i][m] = x[i][m]; + } + } + } else { + for (int m = 0; m < 3; m++) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) nad[i][m] = 0.0; + } + } + } + + if (nad_style == D2MIN) { + xprd0 = domain->xprd; + yprd0 = domain->yprd; + zprd0 = domain->zprd; + xprd0_half = domain->xprd_half; + yprd0_half = domain->yprd_half; + zprd0_half = domain->zprd_half; + xy0 = domain->xy; + xz0 = domain->xz; + yz0 = domain->yz; + } + + reference_saved = 1; +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::calculate_D2Min() +{ + if (!reference_saved) { + error->warning(FLERR, "Calculating D2Min without a saved reference state"); + return; + } + + // invoke half neighbor list (will copy or build if necessary) + neighbor->build_one(list); + + if (atom->nmax > nmax) + grow_arrays(atom->nmax); + + int i, j, k, l, ii, jj, inum, jnum, itype, jtype; + double evol, j2, edev; + double r[3], r0[3], rsq, rsq0, radsum, temp[3]; + double X_tmp[3][3], Y_tmp[3][3], F_tmp[3][3], E[3][3]; + double Y_inv[3][3] = {0.0}; // Zero for 2d since not all entries used + int *ilist, *jlist, *numneigh, **firstneigh; + + double **x = atom->x; + double **x0 = atom->darray[nad_index]; + double *radius = atom->radius; + tagint *tag = atom->tag; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + Pair *pair = force->pair; + double **cutsq; + if (pair) cutsq = force->pair->cutsq; + + for (i = 0; i < nmax; i++) { + for (k = 0; k < 3; k++) { + for (l = 0; l < 3; l++) { + X[i][k][l] = 0.0; + Y[i][k][l] = 0.0; + } + } + norm[i] = 0; + array_atom[i][0] = 0; + } + + // First loop through neighbors + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (!(mask[i] & groupbit)) continue; + + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + if (!(mask[j] & groupbit)) continue; + + jtype = type[j]; + r[0] = x[i][0] - x[j][0]; + r[1] = x[i][1] - x[j][1]; + r[2] = x[i][2] - x[j][2]; + rsq = lensq3(r); + + // Only include contributions from atoms that are CURRENTLY neighbors + if (cut_style == TYPE) { + if (rsq > cutsq[itype][jtype]) continue; + } else if (cut_style == CUSTOM) { + if (rsq > cutsq_custom) continue; + } else { + radsum = radius[i] + radius[j]; + if (rsq > radsum * radsum) continue; + } + + r0[0] = x0[i][0] - x0[j][0]; + r0[1] = x0[i][1] - x0[j][1]; + r0[2] = x0[i][2] - x0[j][2]; + minimum_image0(r0); + + // Using notation from Falk & Langer 1998 + outer3(r, r0, X_tmp); + outer3(r0, r0, Y_tmp); + + for (k = 0; k < 3; k++) { + for (l = 0; l < 3; l++) { + X[i][k][l] += X_tmp[k][l]; + Y[i][k][l] += Y_tmp[k][l]; + } + } + + if (newton_pair || j < nlocal) { + for (k = 0; k < 3; k++) { + for (l = 0; l < 3; l++) { + X[j][k][l] += X_tmp[k][l]; + Y[j][k][l] += Y_tmp[k][l]; + } + } + } + } + } + + comm_flag = 0; + if (newton_pair) comm->reverse_comm(this); + + // Calculate contributions to strain tensor + double denom; + for (i = 0; i < nlocal; i++) { + if (!(mask[i] & groupbit)) continue; + for (j = 0; j < 3; j++) { + for (k = 0; k < 3; k++) { + Y_tmp[j][k] = Y[i][j][k]; + X_tmp[j][k] = X[i][j][k]; + } + } + + if (domain->dimension == 3) { + invert3(Y_tmp, Y_inv); + } else { + denom = Y_tmp[0][0] * Y_tmp[1][1] - Y_tmp[0][1] * Y_tmp[1][0]; + if (denom != 0.0) denom = 1.0 / denom; + Y_inv[0][0] = Y_tmp[1][1] * denom; + Y_inv[0][1] = -Y_tmp[0][1] * denom; + Y_inv[1][0] = -Y_tmp[1][0] * denom; + Y_inv[1][1] = Y_tmp[0][0] * denom; + } + + times3(X_tmp, Y_inv, F_tmp); + + for (j = 0; j < 3; j++) { + for (k = 0; k < 3; k++) { + F[i][j][k] = F_tmp[j][k]; + } + } + } + + comm->forward_comm(this); + + // Second loop through neighbors + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (!(mask[i] & groupbit)) continue; + + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + if (!(mask[j] & groupbit)) continue; + + jtype = type[j]; + r[0] = x[i][0] - x[j][0]; + r[1] = x[i][1] - x[j][1]; + r[2] = x[i][2] - x[j][2]; + rsq = lensq3(r); + + // Only include contributions from atoms that are CURRENTLY neighbors + if (cut_style == TYPE) { + if (rsq >= cutsq[itype][jtype]) continue; + } else if (cut_style == CUSTOM) { + if (rsq >= cutsq_custom) continue; + } else { + radsum = radius[i] + radius[j]; + if (rsq >= radsum * radsum) continue; + } + + r0[0] = x0[i][0] - x0[j][0]; + r0[1] = x0[i][1] - x0[j][1]; + r0[2] = x0[i][2] - x0[j][2]; + minimum_image0(r0); + + // E * r0 + for (k = 0; k < 3; k++) { + temp[k] = 0.0; + for (l = 0; l < 3; l++) + temp[k] += F[i][k][l] * r0[l]; + } + + sub3(r, temp, temp); + array_atom[i][0] += lensq3(temp); + norm[i] += 1; + + if (newton_pair || j < nlocal) { + for (k = 0; k < 3; k++) { + temp[k] = 0.0; + for (l = 0; l < 3; l++) + temp[k] += F[j][k][l] * r0[l]; + } + + sub3(r, temp, temp); + array_atom[j][0] += lensq3(temp); + norm[j] += 1; + } + } + } + + comm_flag = 1; + if (newton_pair) comm->reverse_comm(this, 2); + + for (i = 0; i < nlocal; i++) { + if (!(mask[i] & groupbit)) continue; + + if (norm[i] != 0) + array_atom[i][0] /= norm[i]; + else + array_atom[i][0] = 0.0; + + for (j = 0; j < 3; j++) + for (k = 0; k < 3; k++) + F_tmp[j][k] = F[i][j][k]; + + transpose_times3(F_tmp, F_tmp, E); + for (j = 0; j < 3; j++) E[j][j] -= 1.0; + + evol = (E[0][0] + E[1][1] + E[2][2]) / domain->dimension; + + // Calculate deviatoric strain + for (j = 0; j < 3; j++) E[j][j] -= evol; + j2 = 0.0; + for (j = 0; j < 3; j++) + for (k = 0; k < 3; k++) + j2 += E[j][k] * E[j][k]; + + edev = sqrt(0.5 * j2); + + array_atom[i][1] = evol; + array_atom[i][2] = edev; + } +} + +/* ---------------------------------------------------------------------- */ + +int FixNonaffineDisplacement::pack_reverse_comm(int n, int first, double *buf) +{ + int i, m, last, k, l; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (comm_flag == 0) { + for (k = 0; k < 3; k++) { + for (l = 0; l < 3; l++) { + buf[m++] = X[i][k][l]; + buf[m++] = Y[i][k][l]; + } + } + } else { + buf[m++] = array_atom[i][0]; + buf[m++] = ubuf(norm[i]).d; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i, j, m, k, l; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + if (comm_flag == 0) { + for (k = 0; k < 3; k++) { + for (l = 0; l < 3; l++) { + X[j][k][l] += buf[m++]; + Y[j][k][l] += buf[m++]; + } + } + } else { + array_atom[j][0] += buf[m++]; + norm[j] += (int) ubuf(buf[m++]).i; + } + } +} + +/* ---------------------------------------------------------------------- */ + +int FixNonaffineDisplacement::pack_forward_comm(int n, int *list, double *buf, + int /*pbc_flag*/, int * /*pbc*/) +{ + int i, j, m, k, l; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + for (k = 0; k < 3; k++) { + for (l = 0; l < 3; l ++) { + buf[m++] = F[j][k][l]; + } + } + } + + return m; +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::unpack_forward_comm(int n, int first, double *buf) +{ + int i, m, last, k, l; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + for (k = 0; k < 3; k++) { + for (l = 0; l < 3; l ++) { + F[i][k][l] = buf[m++]; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::minimum_image0(double *delta) +{ + if (domain->triclinic == 0) { + if (domain->xperiodic) { + while (fabs(delta[0]) > xprd0_half) { + if (delta[0] < 0.0) delta[0] += xprd0; + else delta[0] -= xprd0; + } + } + if (domain->yperiodic) { + while (fabs(delta[1]) > yprd0_half) { + if (delta[1] < 0.0) delta[1] += yprd0; + else delta[1] -= yprd0; + } + } + if (domain->zperiodic) { + while (fabs(delta[2]) > zprd0_half) { + if (delta[2] < 0.0) delta[2] += zprd0; + else delta[2] -= zprd0; + } + } + + } else { + if (domain->zperiodic) { + while (fabs(delta[2]) > zprd0_half) { + if (delta[2] < 0.0) { + delta[2] += zprd0; + delta[1] += yz0; + delta[0] += xz0; + } else { + delta[2] -= zprd0; + delta[1] -= yz0; + delta[0] -= xz0; + } + } + } + if (domain->yperiodic) { + while (fabs(delta[1]) > yprd0_half) { + if (delta[1] < 0.0) { + delta[1] += yprd0; + delta[0] += xy0; + } else { + delta[1] -= yprd0; + delta[0] -= xy0; + } + } + } + if (domain->xperiodic) { + while (fabs(delta[0]) > xprd0_half) { + if (delta[0] < 0.0) delta[0] += xprd0; + else delta[0] -= xprd0; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::grow_arrays(int nmax_new) +{ + nmax = nmax_new; + memory->destroy(array_atom); + memory->destroy(X); + memory->destroy(Y); + memory->destroy(F); + memory->destroy(norm); + memory->create(array_atom, nmax, 3, "fix_nonaffine_displacement:array_atom"); + memory->create(X, nmax, 3, 3, "fix_nonaffine_displacement:X"); + memory->create(Y, nmax, 3, 3, "fix_nonaffine_displacement:Y"); + memory->create(F, nmax, 3, 3, "fix_nonaffine_displacement:F"); + memory->create(norm, nmax, "fix_nonaffine_displacement:norm"); +} diff --git a/src/fix_nonaffine_displacement.h b/src/fix_nonaffine_displacement.h new file mode 100755 index 0000000000..9a914cf5df --- /dev/null +++ b/src/fix_nonaffine_displacement.h @@ -0,0 +1,70 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS +// clang-format off +FixStyle(nonaffine/displacement,FixNonaffineDisplacement) +// clang-format on +#else + +#ifndef LMP_FIX_NONAFFINE_DISPLACEMENT_H +#define LMP_FIX_NONAFFINE_DISPLACEMENT_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixNonaffineDisplacement : public Fix { + public: + FixNonaffineDisplacement(class LAMMPS *, int, char **); + ~FixNonaffineDisplacement() override; + int setmask() override; + void post_constructor() override; + void init() override; + void init_list(int, class NeighList *) override; + void setup(int); + void post_force(int) override; + void write_restart(FILE *fp) override; + void restart(char *buf) override; + int pack_forward_comm(int, int *, double *, int, int *) override; + void unpack_forward_comm(int, int, double *) override; + int pack_reverse_comm(int, int, double *) override; + void unpack_reverse_comm(int, int *, double *) override; + + private: + double dtv; + char *new_fix_id; + int nad_index, nmax, comm_flag; + int nad_style, cut_style; + int reference_style, offset_timestep, reference_timestep, update_timestep; + int reference_saved; + double cutoff_custom, cutsq_custom, mycutneigh; + double xprd0, yprd0, zprd0, xprd0_half, yprd0_half, zprd0_half, xy0, xz0, yz0; + + double ***X, ***Y, ***F; + int *norm; + + class NeighList *list; // half neighbor list + + + void integrate_velocity(); + void calculate_D2Min(); + void save_reference_state(); + void minimum_image0(double *); + void grow_arrays(int); +}; + +} // namespace LAMMPS_NS + +#endif +#endif From 2fbaeb4fc7033484f5498c541e3c41430b9873a8 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 15 Jun 2023 16:03:30 -0600 Subject: [PATCH 02/18] Fixing merge conflict --- doc/src/compute_contact_atom.rst | 3 --- 1 file changed, 3 deletions(-) diff --git a/doc/src/compute_contact_atom.rst b/doc/src/compute_contact_atom.rst index 8e94460877..384ede290f 100644 --- a/doc/src/compute_contact_atom.rst +++ b/doc/src/compute_contact_atom.rst @@ -50,9 +50,6 @@ overview of LAMMPS output options. The per-atom vector values will be a number :math:`\ge 0.0`, as explained above. -The optional *group2-ID* argument allows to specify from which group atoms -contribute to the coordination number. Default setting is group 'all.' - Restrictions """""""""""" From 7012e6ddd466684806832b5f2bbbd253d491048b Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 15 Jun 2023 19:06:14 -0600 Subject: [PATCH 03/18] Debugging and cleaning up D2min calculation --- src/fix_nonaffine_displacement.cpp | 38 +++++++++++++++--------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/fix_nonaffine_displacement.cpp b/src/fix_nonaffine_displacement.cpp index 15c678f83d..e5c5eba4e2 100755 --- a/src/fix_nonaffine_displacement.cpp +++ b/src/fix_nonaffine_displacement.cpp @@ -105,7 +105,7 @@ FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char * } else if (strcmp(arg[iarg], "offset") == 0) { reference_style = OFFSET; offset_timestep = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); - if (offset_timestep < 0) + if ((offset_timestep <= 0) || (offset_timestep > nevery)) error->all(FLERR, "Illegal offset timestep {} in fix nonaffine/displacement", arg[iarg + 1]); } else error->all(FLERR,"Illegal reference style {} in fix nonaffine/displacement", arg[iarg]); @@ -237,6 +237,14 @@ void FixNonaffineDisplacement::setup(int vflag) void FixNonaffineDisplacement::post_force(int /*vflag*/) { + if (reference_saved && !update->setupflag) { + if (nad_style == INTEGRATED) { + integrate_velocity(); + } else { + if ((update->ntimestep % nevery) == 0) calculate_D2Min(); + } + } + if (reference_style == FIXED) if (update->ntimestep == reference_timestep) save_reference_state(); @@ -246,18 +254,8 @@ void FixNonaffineDisplacement::post_force(int /*vflag*/) save_reference_state(); if (reference_style == OFFSET) - if (((update->ntimestep - offset_timestep) % nevery) == 0) + if (((update->ntimestep + offset_timestep) % nevery) == 0) save_reference_state(); - - if (update->setupflag) return; - - if (reference_style == FIXED && (update->ntimestep < reference_timestep)) return; - - if (nad_style == INTEGRATED) { - integrate_velocity(); - } else { - if ((update->ntimestep % nevery) == 0) calculate_D2Min(); - } } /* ---------------------------------------------------------------------- */ @@ -309,17 +307,18 @@ void FixNonaffineDisplacement::save_reference_state() int *mask = atom->mask; int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; double **nad = atom->darray[nad_index]; if (nad_style == D2MIN) { for (int m = 0; m < 3; m++) { - for (int i = 0; i < nlocal; i++) { + for (int i = 0; i < nall; i++) { if (mask[i] & groupbit) nad[i][m] = x[i][m]; } } } else { for (int m = 0; m < 3; m++) { - for (int i = 0; i < nlocal; i++) { + for (int i = 0; i < nall; i++) { if (mask[i] & groupbit) nad[i][m] = 0.0; } } @@ -370,6 +369,7 @@ void FixNonaffineDisplacement::calculate_D2Min() int *mask = atom->mask; int nlocal = atom->nlocal; int newton_pair = force->newton_pair; + int dim = domain->dimension; inum = list->inum; ilist = list->ilist; @@ -450,7 +450,7 @@ void FixNonaffineDisplacement::calculate_D2Min() } comm_flag = 0; - if (newton_pair) comm->reverse_comm(this); + if (newton_pair) comm->reverse_comm(this, 18); // Calculate contributions to strain tensor double denom; @@ -463,7 +463,7 @@ void FixNonaffineDisplacement::calculate_D2Min() } } - if (domain->dimension == 3) { + if (dim == 3) { invert3(Y_tmp, Y_inv); } else { denom = Y_tmp[0][0] * Y_tmp[1][1] - Y_tmp[0][1] * Y_tmp[1][0]; @@ -562,12 +562,12 @@ void FixNonaffineDisplacement::calculate_D2Min() F_tmp[j][k] = F[i][j][k]; transpose_times3(F_tmp, F_tmp, E); - for (j = 0; j < 3; j++) E[j][j] -= 1.0; + for (j = 0; j < dim; j++) E[j][j] -= 1.0; - evol = (E[0][0] + E[1][1] + E[2][2]) / domain->dimension; + evol = (E[0][0] + E[1][1] + E[2][2]) / dim; // Calculate deviatoric strain - for (j = 0; j < 3; j++) E[j][j] -= evol; + for (j = 0; j < dim; j++) E[j][j] -= evol; j2 = 0.0; for (j = 0; j < 3; j++) for (k = 0; k < 3; k++) From a1513a7d3b04c760fff222a275c7e18d09afd540 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 15 Jun 2023 23:06:50 -0600 Subject: [PATCH 04/18] rattler documentation --- doc/src/compute_contact_atom.rst | 3 -- doc/src/compute_rattlers.rst | 65 +++++++++++++------------- doc/src/fix_nonaffine_displacement.rst | 10 ++-- 3 files changed, 36 insertions(+), 42 deletions(-) diff --git a/doc/src/compute_contact_atom.rst b/doc/src/compute_contact_atom.rst index 384ede290f..b7ed062ff6 100644 --- a/doc/src/compute_contact_atom.rst +++ b/doc/src/compute_contact_atom.rst @@ -69,6 +69,3 @@ Default """"""" *group2-ID* = all - - -none diff --git a/doc/src/compute_rattlers.rst b/doc/src/compute_rattlers.rst index c26f2b4252..a4ae2d3d16 100644 --- a/doc/src/compute_rattlers.rst +++ b/doc/src/compute_rattlers.rst @@ -1,6 +1,6 @@ -.. index:: compute contact/atom +.. index:: compute rattlers -compute contact/atom command +compute rattlers command ============================ Syntax @@ -8,67 +8,68 @@ Syntax .. parsed-literal:: - compute ID group-ID contact/atom group2-ID + compute ID group-ID rattlers cutoff zmin ntries * ID, group-ID are documented in :doc:`compute ` command -* contact/atom = style name of this compute command -* group2-ID = optional argument select group-ID to restrict which atoms to consider for contacts (see below) +* rattlers = style name of this compute command +* cutoff = *type* or *radius* + + .. parsed-literal:: + + *type* = cutoffs determined based on atom types + *radius* = cutoffs determined based on atom diameters (atom style sphere) + +* zmin = minimum coordination for a non-rattler particle +* ntries = maximum number of interations to remove rattlers Examples """""""" .. code-block:: LAMMPS - compute 1 all contact/atom - compute 1 all contact/atom mygroup + compute 1 all rattlers type 4 10 Description """"""""""" -Define a computation that calculates the number of contacts -for each atom in a group. - -The contact number is defined for finite-size spherical particles as -the number of neighbor atoms which overlap the central particle, -meaning that their distance of separation is less than or equal to the -sum of the radii of the two particles. - -The value of the contact number will be 0.0 for atoms not in the -specified compute group. - -The optional *group2-ID* argument allows to specify from which group atoms -contribute to the coordination number. Default setting is group 'all'. +Define a compute that identifies rattlers in a system. Rattlers are +identified using an interative approach. The coordination number of +all atoms is first calculated. The *type* and *radius* settings are +used to select whether interactions cutoffs are determined by atom +types or by the sum of atomic radii (atom style sphere), respectively. +Rattlers are then identified as particles with a coordination number +less than *zmin* and are removed from consideration. A coordination +number is then recalaculated, excluding previously identified rattlers, +to identify new rattlers. This process is looped, up to a maximum +of *ntries*, until no new rattlers are identified and the remaining +atoms form a stable network of contacts. Output info """"""""""" -This compute calculates a per-atom vector, whose values can be -accessed by any command that uses per-atom values from a compute as -input. See the :doc:`Howto output ` page for an +This compute calculates a per-atom vector and a global scalar. The vector +designates which atoms are rattlers, indicated by a value 1. Non-rattlers +have a value of 0. The global scalar returns the total number of rattlers +in the system. See the :doc:`Howto output ` page for an overview of LAMMPS output options. -The per-atom vector values will be a number >= 0.0, as explained -above. - Restrictions """""""""""" -This compute is part of the GRANULAR package. It is only enabled if +This compute is part of the EXTRA-COMPUTE package. It is only enabled if LAMMPS was built with that package. See the :doc:`Build package ` page for more info. -This compute requires that atoms store a radius as defined by the -:doc:`atom_style sphere ` command. +The *radius* cutoff option requires that atoms store a radius as defined by the +:doc:`atom_style sphere ` or similar commands. Related commands """""""""""""""" :doc:`compute coord/atom ` +:doc:`compute contact/atom ` Default """"""" -*group2-ID* = all - - none diff --git a/doc/src/fix_nonaffine_displacement.rst b/doc/src/fix_nonaffine_displacement.rst index 7b0f7669ed..51084c9501 100644 --- a/doc/src/fix_nonaffine_displacement.rst +++ b/doc/src/fix_nonaffine_displacement.rst @@ -1,18 +1,14 @@ -.. index:: fix gravity -.. index:: fix gravity/omp -.. index:: fix gravity/kk +.. index:: fix nonaffine/displacement -fix gravity command +fix nonaffine/displacement command =================== -Accelerator Variants: *gravity/omp*, *gravity/kk* - Syntax """""" .. parsed-literal:: - fix ID group gravity magnitude style args + fix ID group nonaffine/displacement style args reference/style args * ID, group are documented in :doc:`fix ` command * gravity = style name of this fix command From ca636ffa7d44de7f218c2c439ed3133c861aaccc Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 16 Jun 2023 16:28:48 -0600 Subject: [PATCH 05/18] Adding d2min doc --- doc/src/fix_nonaffine_displacement.rst | 149 +++++++++++-------------- src/GRANULAR/compute_fabric.cpp | 4 + src/compute_rattlers.cpp | 4 + src/fix_nonaffine_displacement.cpp | 29 ++--- 4 files changed, 91 insertions(+), 95 deletions(-) diff --git a/doc/src/fix_nonaffine_displacement.rst b/doc/src/fix_nonaffine_displacement.rst index 51084c9501..fc3a51b1b6 100644 --- a/doc/src/fix_nonaffine_displacement.rst +++ b/doc/src/fix_nonaffine_displacement.rst @@ -8,129 +8,114 @@ Syntax .. parsed-literal:: - fix ID group nonaffine/displacement style args reference/style args + fix ID group nonaffine/displacement style args reference/style nstep * ID, group are documented in :doc:`fix ` command -* gravity = style name of this fix command -* magnitude = size of acceleration (force/mass units) -* magnitude can be a variable (see below) -* style = *chute* or *spherical* or *gradient* or *vector* +* nonaffine/displacement = style name of this fix command +* nevery = calculate nonaffine displacement every this many timesteps +* style = *d2min* or *integrated* .. parsed-literal:: - *chute* args = angle - angle = angle in +x away from -z or -y axis in 3d/2d (in degrees) - angle can be a variable (see below) - *spherical* args = phi theta - phi = azimuthal angle from +x axis (in degrees) - theta = angle from +z or +y axis in 3d/2d (in degrees) - phi or theta can be a variable (see below) - *vector* args = x y z - x y z = vector direction to apply the acceleration - x or y or z can be a variable (see below) + *d2min* args = cutoff args + cutoff = *type* or *radius* or *custom* + *type* args = none, cutoffs determined by atom types + *radius* args = none, cutoffs determined based on atom diameters (atom style sphere) + *custom* args = *rmax*, cutoff set by a constant numeric value *rmax* + *integrated* args = none + +* reference/style = *fixed* or *update* or *offset* + + .. parsed-literal:: + + *fixed* = use a fixed reference frame at *nstep* + *update* = update the reference frame every *nstep* timesteps + *offset* = update the reference frame *nstep* timesteps before calculating the non/affine displacement Examples """""""" .. code-block:: LAMMPS - fix 1 all gravity 1.0 chute 24.0 - fix 1 all gravity v_increase chute 24.0 - fix 1 all gravity 1.0 spherical 0.0 -180.0 - fix 1 all gravity 10.0 spherical v_phi v_theta - fix 1 all gravity 100.0 vector 1 1 0 + fix 1 all nonaffine/displacement 100 integrated update 100 + fix 1 all nonaffine/displacement 1000 d2min type fixed 0 + fix 1 all nonaffine/displacement 1000 d2min custom 2.0 offset 100 Description """"""""""" -Impose an additional acceleration on each particle in the group. This -fix is typically used with granular systems to include a "gravity" -term acting on the macroscopic particles. More generally, it can -represent any kind of driving field, e.g. a pressure gradient inducing -a Poiseuille flow in a fluid. Note that this fix operates differently -than the :doc:`fix addforce ` command. The addforce fix -adds the same force to each atom, independent of its mass. This -command imparts the same acceleration to each atom (force/mass). +This fix computes different metrics of the nonaffine displacement of +particles. The first metric, *d2min* calculates the :math:`D^2_\mathrm{min}` +nonaffine displacement by Falk and Langer in :ref:`(Falk) `. +For each atom, the fix computes the two tensors -The *magnitude* of the acceleration is specified in force/mass units. -For granular systems (LJ units) this is typically 1.0. See the -:doc:`units ` command for details. +.. math:: -Style *chute* is typically used for simulations of chute flow where -the specified *angle* is the chute angle, with flow occurring in the +x -direction. For 3d systems, the tilt is away from the z axis; for 2d -systems, the tilt is away from the y axis. + X = \sum_{\mathrm{neighbors}} \vec{r} \left(\vec{r}_{0} \right)^T -Style *spherical* allows an arbitrary 3d direction to be specified for -the acceleration vector. *Phi* and *theta* are defined in the usual -spherical coordinates. Thus for acceleration acting in the -z -direction, *theta* would be 180.0 (or -180.0). *Theta* = 90.0 and -*phi* = -90.0 would mean acceleration acts in the -y direction. For -2d systems, *phi* is ignored and *theta* is an angle in the xy plane -where *theta* = 0.0 is the y-axis. +and -Style *vector* imposes an acceleration in the vector direction given -by (x,y,z). Only the direction of the vector is important; it's -length is ignored. For 2d systems, the *z* component is ignored. +.. math:: -Any of the quantities *magnitude*, *angle*, *phi*, *theta*, *x*, *y*, -*z* which define the gravitational magnitude and direction, can be -specified as an equal-style :doc:`variable `. If the value 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 determine the quantity. You should -insure that the variable calculates a result in the appropriate units, -e.g. force/mass or degrees. + Y = \sum_{\mathrm{neighbors}} \vec{r}_0 \left(\vec{r}_{0} \right)^T -Equal-style variables can specify formulas with various mathematical -functions, and include :doc:`thermo_style ` command -keywords for the simulation box parameters and timestep and elapsed -time. Thus it is easy to specify a time-dependent gravitational -field. +where the neighbors include all other atoms within the distance criterion +set by the cutoff option, discussed below, :math:`\vec{r}` is the current +displacement between particles, and :math:`\vec{r}_0` is the reference +displacement. A deformation gradient tensor is then calculated as +:math:`F = X Y^{-1}` from which ----------- +.. math:: -.. include:: accel_styles.rst + D^2_\mathrm{min} = \sum_{\mathrm{neighbors}} \left| \vec{r} - F \vec{r}_0 \right|^2 + +and a strain tensor is calculated :math:`E = F F^{T} - I` where :math:`I` +is the identity tensor. + +The *integrated* style simply integrates the velocity of particles +every timestep to calculate a displacement. This style only works if +used in conjunction with another fix that deforms the box and displaces +atom positions such as :doc:`the remap x option of fix deform `, +:doc:`fix press/berendsen `, or :doc:`fix nh `. ---------- Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" -No information about this fix is written to :doc:`binary restart files `. +The reference state is saved to :doc:`binary restart files `. -The :doc:`fix_modify ` *energy* option is supported by -this fix to add the gravitational potential energy of the system to -the global potential energy of the system as part of -:doc:`thermodynamic output `. The default setting for -this fix is :doc:`fix_modify energy no `. +None of the :doc:`fix_modify ` options are relevant to this +fix. -The :doc:`fix_modify ` *respa* option is supported by this -fix. This allows to set at which level of the :doc:`r-RESPA -` integrator the fix is adding its forces. Default is the -outermost level. +This fix computes a peratom array with 3 columns, which can be accessed +by indices 1-3 using any command that uses per-atom values from a fix +as input. -This fix computes a global scalar which can be accessed by various -:doc:`output commands `. This scalar is the -gravitational potential energy of the particles in the defined field, -namely mass \* (g dot x) for each particles, where x and mass are the -particles position and mass, and g is the gravitational field. The -scalar value calculated by this fix is "extensive". - -No parameter of this fix can be used with the *start/stop* keywords of -the :doc:`run ` command. This fix is not invoked during -:doc:`energy minimization `. +For the *integrated* style, the three columns are the nonaffine +displacements in the x, y, and z directions. For the *d2min* style, +the three columns are the calculated :math:`D^2_\mathrm{min}`, the +volumetric strain, and the deviatoric strain. Restrictions """""""""""" - none + +This compute is part of the EXTRA-FIX package. It is only enabled if +LAMMPS was built with that package. See the +:doc:`Build package ` page for more info. Related commands """""""""""""""" -:doc:`atom_style sphere `, :doc:`fix addforce ` +none Default """"""" none + +---------- + +.. _nh-Martyna: + +**(Martyna)** Martyna, Tobias and Klein, J Chem Phys, 101, 4177 (1994). \ No newline at end of file diff --git a/src/GRANULAR/compute_fabric.cpp b/src/GRANULAR/compute_fabric.cpp index fb95a8b446..adaf242c92 100644 --- a/src/GRANULAR/compute_fabric.cpp +++ b/src/GRANULAR/compute_fabric.cpp @@ -11,6 +11,10 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Joel Clemmer (SNL), Ishan Srivastava (LBNL) +------------------------------------------------------------------------- */ + #include "compute_fabric.h" #include "atom.h" diff --git a/src/compute_rattlers.cpp b/src/compute_rattlers.cpp index 5a6d2b022d..e7e99a7471 100755 --- a/src/compute_rattlers.cpp +++ b/src/compute_rattlers.cpp @@ -11,6 +11,10 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Joel Clemmer (SNL), Ishan Srivastava (LBNL) +------------------------------------------------------------------------- */ + #include "compute_rattlers.h" #include "atom.h" diff --git a/src/fix_nonaffine_displacement.cpp b/src/fix_nonaffine_displacement.cpp index e5c5eba4e2..425800fb45 100755 --- a/src/fix_nonaffine_displacement.cpp +++ b/src/fix_nonaffine_displacement.cpp @@ -12,6 +12,10 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Joel Clemmer (SNL), Ishan Srivastava (LBNL) +------------------------------------------------------------------------- */ + #include "fix_nonaffine_displacement.h" #include "atom.h" @@ -63,31 +67,30 @@ FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char * { if (narg < 4) error->all(FLERR,"Illegal fix nonaffine/displacement command"); - int iarg = 3; + nevery = utils::inumeric(FLERR, arg[3], false, lmp); + if (nevery <= 0) error->all(FLERR,"Illegal nevery value {} in fix nonaffine/displacement", nevery); + + int iarg = 4; if (strcmp(arg[iarg], "integrated") == 0) { nad_style = INTEGRATED; nevery = 1; iarg += 1; } else if (strcmp(arg[iarg], "d2min") == 0) { - if (iarg + 2 > narg) error->all(FLERR,"Illegal fix nonaffine/displacement command"); - + if (iarg + 1 > narg) error->all(FLERR,"Illegal fix nonaffine/displacement command"); nad_style = D2MIN; - nevery = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); - if (nevery <= 0) error->all(FLERR,"Illegal nevery value {} in fix nonaffine/displacement", nevery); - - if (strcmp(arg[iarg + 2], "type") == 0) { + if (strcmp(arg[iarg + 1], "type") == 0) { cut_style = TYPE; - } else if (strcmp(arg[iarg + 2], "radius") == 0) { + } else if (strcmp(arg[iarg + 1], "radius") == 0) { cut_style = RADIUS; - } else if (strcmp(arg[iarg + 2], "custom") == 0) { - if (iarg + 3 > narg) error->all(FLERR,"Illegal fix nonaffine/displacement command"); + } else if (strcmp(arg[iarg + 1], "custom") == 0) { + if (iarg + 2 > narg) error->all(FLERR,"Illegal fix nonaffine/displacement command"); cut_style = CUSTOM; - cutoff_custom = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + cutoff_custom = utils::numeric(FLERR, arg[iarg + 2], false, lmp); cutsq_custom = cutoff_custom * cutoff_custom; if (cutoff_custom <= 0) - error->all(FLERR, "Illegal custom cutoff length {}", arg[iarg + 3]); + error->all(FLERR, "Illegal custom cutoff length {}", arg[iarg + 2]); iarg += 1; - } else error->all(FLERR,"Illegal cutoff style {} in fix nonaffine/displacement", arg[iarg + 2]); + } else error->all(FLERR,"Illegal cutoff style {} in fix nonaffine/displacement", arg[iarg + 1]); iarg += 3; } From 2a432bdaf1957df008e548f18744c787ee8ffe71 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 16 Jun 2023 16:32:15 -0600 Subject: [PATCH 06/18] Clarifying output as sqrt --- doc/src/fix_nonaffine_displacement.rst | 4 ++-- src/fix_nonaffine_displacement.cpp | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/doc/src/fix_nonaffine_displacement.rst b/doc/src/fix_nonaffine_displacement.rst index fc3a51b1b6..892a33bb64 100644 --- a/doc/src/fix_nonaffine_displacement.rst +++ b/doc/src/fix_nonaffine_displacement.rst @@ -21,7 +21,7 @@ Syntax cutoff = *type* or *radius* or *custom* *type* args = none, cutoffs determined by atom types *radius* args = none, cutoffs determined based on atom diameters (atom style sphere) - *custom* args = *rmax*, cutoff set by a constant numeric value *rmax* + *custom* args = *rmax*, cutoff set by a constant numeric value *rmax* (distance units) *integrated* args = none * reference/style = *fixed* or *update* or *offset* @@ -94,7 +94,7 @@ as input. For the *integrated* style, the three columns are the nonaffine displacements in the x, y, and z directions. For the *d2min* style, -the three columns are the calculated :math:`D^2_\mathrm{min}`, the +the three columns are the calculated :math:`\sqrt{D^2_\mathrm{min}}`, the volumetric strain, and the deviatoric strain. Restrictions diff --git a/src/fix_nonaffine_displacement.cpp b/src/fix_nonaffine_displacement.cpp index 425800fb45..24f7b2e70b 100755 --- a/src/fix_nonaffine_displacement.cpp +++ b/src/fix_nonaffine_displacement.cpp @@ -559,6 +559,7 @@ void FixNonaffineDisplacement::calculate_D2Min() array_atom[i][0] /= norm[i]; else array_atom[i][0] = 0.0; + array_atom[i][0] = sqrt(array_atom[i][0]); for (j = 0; j < 3; j++) for (k = 0; k < 3; k++) From 7fa1f4b3b4e5d8b908526a4246ce540413e93272 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Sat, 15 Jul 2023 22:10:19 -0600 Subject: [PATCH 07/18] Various clean ups, moving files --- src/{ => EXTRA-COMPUTE}/compute_rattlers.cpp | 8 +-- src/{ => EXTRA-COMPUTE}/compute_rattlers.h | 0 .../fix_nonaffine_displacement.cpp | 51 ++++++++++--------- .../fix_nonaffine_displacement.h | 0 4 files changed, 30 insertions(+), 29 deletions(-) rename src/{ => EXTRA-COMPUTE}/compute_rattlers.cpp (98%) rename src/{ => EXTRA-COMPUTE}/compute_rattlers.h (100%) rename src/{ => EXTRA-FIX}/fix_nonaffine_displacement.cpp (94%) rename src/{ => EXTRA-FIX}/fix_nonaffine_displacement.h (100%) diff --git a/src/compute_rattlers.cpp b/src/EXTRA-COMPUTE/compute_rattlers.cpp similarity index 98% rename from src/compute_rattlers.cpp rename to src/EXTRA-COMPUTE/compute_rattlers.cpp index e7e99a7471..5613afb707 100755 --- a/src/compute_rattlers.cpp +++ b/src/EXTRA-COMPUTE/compute_rattlers.cpp @@ -113,7 +113,7 @@ void ComputeRattlers::compute_peratom() int i, j, ii, jj, inum, jnum, itype, jtype, tmp_flag; tagint itag, jtag; double xtmp, ytmp, ztmp, delx, dely, delz; - double r, rinv, rsq, radsum; + double rsq, radsum; if (nmax < atom->nmax) { nmax = atom->nmax; @@ -150,8 +150,6 @@ void ComputeRattlers::compute_peratom() int change_flag = 1; int ntry = 0; while (ntry < max_tries) { - // Quit when answer stops evolving - if (change_flag == 0) break; change_flag = 0; for (i = 0; i < nmax; i++) ncontacts[i] = 0; @@ -227,6 +225,7 @@ void ComputeRattlers::compute_peratom() MPI_Allreduce(&change_flag, &tmp_flag, 1, MPI_INT, MPI_MAX, world); change_flag = tmp_flag; + if (change_flag == 0) break; ntry += 1; } @@ -237,7 +236,8 @@ void ComputeRattlers::compute_peratom() /* ---------------------------------------------------------------------- */ -double ComputeRattlers::compute_scalar() { +double ComputeRattlers::compute_scalar() +{ if (invoked_peratom != update->ntimestep) compute_peratom(); diff --git a/src/compute_rattlers.h b/src/EXTRA-COMPUTE/compute_rattlers.h similarity index 100% rename from src/compute_rattlers.h rename to src/EXTRA-COMPUTE/compute_rattlers.h diff --git a/src/fix_nonaffine_displacement.cpp b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp similarity index 94% rename from src/fix_nonaffine_displacement.cpp rename to src/EXTRA-FIX/fix_nonaffine_displacement.cpp index 24f7b2e70b..f18fdd2c4a 100755 --- a/src/fix_nonaffine_displacement.cpp +++ b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp @@ -63,7 +63,7 @@ static const char cite_nonaffine_d2min[] = /* ---------------------------------------------------------------------- */ FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), new_fix_id(nullptr), X(nullptr), Y(nullptr), F(nullptr), norm (nullptr) + Fix(lmp, narg, arg), new_fix_id(nullptr), X(nullptr), Y(nullptr), F(nullptr), norm(nullptr) { if (narg < 4) error->all(FLERR,"Illegal fix nonaffine/displacement command"); @@ -91,8 +91,8 @@ FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char * error->all(FLERR, "Illegal custom cutoff length {}", arg[iarg + 2]); iarg += 1; } else error->all(FLERR,"Illegal cutoff style {} in fix nonaffine/displacement", arg[iarg + 1]); - iarg += 3; - } + iarg += 2; + } else error->all(FLERR,"Illegal nonaffine displacement style {} in fix nonaffine/displacement", arg[iarg]); if (iarg + 2 > narg) error->all(FLERR,"Illegal fix nonaffine/displacement command"); if (strcmp(arg[iarg], "fixed") == 0) { @@ -112,8 +112,9 @@ FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char * error->all(FLERR, "Illegal offset timestep {} in fix nonaffine/displacement", arg[iarg + 1]); } else error->all(FLERR,"Illegal reference style {} in fix nonaffine/displacement", arg[iarg]); - if (cut_style == RADIUS && (!atom->radius_flag)) - error->all(FLERR, "Fix nonaffine/displacement radius style requires atom attribute radius"); + if (nad_style == D2MIN) + if (cut_style == RADIUS && (!atom->radius_flag)) + error->all(FLERR, "Fix nonaffine/displacement radius style requires atom attribute radius"); peratom_flag = 1; peratom_freq = nevery; @@ -169,7 +170,7 @@ void FixNonaffineDisplacement::post_constructor() new_fix_id = utils::strdup(id + std::string("_FIX_PA")); modify->add_fix(fmt::format("{} {} property/atom d2_nad 3 ghost {}", new_fix_id, group->names[igroup], ghost_status)); int tmp1, tmp2; - nad_index = atom->find_custom("nad",tmp1,tmp2); + nad_index = atom->find_custom("nad", tmp1, tmp2); if (nad_style == INTEGRATED) { double **nad = atom->darray[nad_index]; @@ -189,35 +190,35 @@ void FixNonaffineDisplacement::init() { dtv = update->dt; - if (!reference_saved && (reference_style == FIXED) && (update->ntimestep > reference_timestep)) + if ((!reference_saved) && (reference_style == FIXED) && (update->ntimestep > reference_timestep)) error->all(FLERR, "Initial timestep exceeds that of the reference state in fix nonaffine/displacement"); if (nad_style == D2MIN) { - if (!force->pair && (cut_style == TYPE)) + if ((!force->pair) && (cut_style == TYPE)) error->all(FLERR,"Fix nonaffine/displacement D2Min option requires a pair style be defined " "or cutoff specified"); - if (cut_style == CUSTOM) { - double skin = neighbor->skin; - mycutneigh = cutoff_custom + skin; - - double cutghost; // as computed by Neighbor and Comm - if (force->pair) - cutghost = MAX(force->pair->cutforce + skin, comm->cutghostuser); - else - cutghost = comm->cutghostuser; - - if (mycutneigh > cutghost) - error->all(FLERR,"Fix nonaffine/displacement D2Min option cutoff exceeds ghost atom range - use comm_modify cutoff command"); - } - // need an occasional half neighbor list if (cut_style == RADIUS) { auto req = neighbor->add_request(this, NeighConst::REQ_SIZE | NeighConst::REQ_OCCASIONAL); } else { auto req = neighbor->add_request(this, NeighConst::REQ_OCCASIONAL); - if (cut_style == CUSTOM) req->set_cutoff(mycutneigh); + if (cut_style == CUSTOM) { + double skin = neighbor->skin; + mycutneigh = cutoff_custom + skin; + + double cutghost; // as computed by Neighbor and Comm + if (force->pair) + cutghost = MAX(force->pair->cutforce + skin, comm->cutghostuser); + else + cutghost = comm->cutghostuser; + + if (mycutneigh > cutghost) + error->all(FLERR,"Fix nonaffine/displacement D2Min option cutoff exceeds ghost atom range - use comm_modify cutoff command"); + + req->set_cutoff(mycutneigh); + } } } } @@ -240,7 +241,7 @@ void FixNonaffineDisplacement::setup(int vflag) void FixNonaffineDisplacement::post_force(int /*vflag*/) { - if (reference_saved && !update->setupflag) { + if (reference_saved && (!update->setupflag)) { if (nad_style == INTEGRATED) { integrate_velocity(); } else { @@ -422,7 +423,7 @@ void FixNonaffineDisplacement::calculate_D2Min() if (rsq > cutsq_custom) continue; } else { radsum = radius[i] + radius[j]; - if (rsq > radsum * radsum) continue; + if (rsq > (radsum * radsum)) continue; } r0[0] = x0[i][0] - x0[j][0]; diff --git a/src/fix_nonaffine_displacement.h b/src/EXTRA-FIX/fix_nonaffine_displacement.h similarity index 100% rename from src/fix_nonaffine_displacement.h rename to src/EXTRA-FIX/fix_nonaffine_displacement.h From a80739c537f9711c88900bae784ab1ad8e01b8fc Mon Sep 17 00:00:00 2001 From: jtclemm Date: Sat, 15 Jul 2023 22:12:10 -0600 Subject: [PATCH 08/18] Doc and gitignore updates --- doc/src/compute_rattlers.rst | 10 +++++----- doc/src/fix_nonaffine_displacement.rst | 2 +- src/.gitignore | 4 ++++ 3 files changed, 10 insertions(+), 6 deletions(-) diff --git a/doc/src/compute_rattlers.rst b/doc/src/compute_rattlers.rst index a4ae2d3d16..a0f0d1b54a 100644 --- a/doc/src/compute_rattlers.rst +++ b/doc/src/compute_rattlers.rst @@ -20,7 +20,7 @@ Syntax *radius* = cutoffs determined based on atom diameters (atom style sphere) * zmin = minimum coordination for a non-rattler particle -* ntries = maximum number of interations to remove rattlers +* ntries = maximum number of iterations to remove rattlers Examples """""""" @@ -35,12 +35,12 @@ Description Define a compute that identifies rattlers in a system. Rattlers are identified using an interative approach. The coordination number of all atoms is first calculated. The *type* and *radius* settings are -used to select whether interactions cutoffs are determined by atom +used to select whether interaction cutoffs are determined by atom types or by the sum of atomic radii (atom style sphere), respectively. Rattlers are then identified as particles with a coordination number -less than *zmin* and are removed from consideration. A coordination -number is then recalaculated, excluding previously identified rattlers, -to identify new rattlers. This process is looped, up to a maximum +less than *zmin* and are removed from consideration. Atomic coordination +numbers are then recalculated, excluding previously identified rattlers, +to identify a new set of rattlers. This process is iterated, up to a maximum of *ntries*, until no new rattlers are identified and the remaining atoms form a stable network of contacts. diff --git a/doc/src/fix_nonaffine_displacement.rst b/doc/src/fix_nonaffine_displacement.rst index 892a33bb64..466a7b1ab0 100644 --- a/doc/src/fix_nonaffine_displacement.rst +++ b/doc/src/fix_nonaffine_displacement.rst @@ -30,7 +30,7 @@ Syntax *fixed* = use a fixed reference frame at *nstep* *update* = update the reference frame every *nstep* timesteps - *offset* = update the reference frame *nstep* timesteps before calculating the non/affine displacement + *offset* = update the reference frame *nstep* timesteps before calculating the nonaffine displacement Examples """""""" diff --git a/src/.gitignore b/src/.gitignore index 2cb2fd315b..006b468c08 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -598,6 +598,8 @@ /compute_pressure_grem.h /compute_ptm_atom.cpp /compute_ptm_atom.h +/compute_rattlers.cpp +/compute_rattlers.h /compute_rigid_local.cpp /compute_rigid_local.h /compute_smd_triangle_vertices.cpp @@ -889,6 +891,8 @@ /fix_nvt_sllod_eff.h /fix_nve_tri.cpp /fix_nve_tri.h +/fix_nonaffine_displacement.cpp +/fix_nonaffine_displacement.h /fix_oneway.cpp /fix_oneway.h /fix_orient_bcc.cpp From f4000efd8a43e8131a3669e445da62699e253c6e Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 17 Jul 2023 11:52:06 -0600 Subject: [PATCH 09/18] Adding reference section to doc pages, update misc error messages --- doc/src/fix_nonaffine_displacement.rst | 12 +++++++++++- src/EXTRA-FIX/fix_nonaffine_displacement.cpp | 8 +++----- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/doc/src/fix_nonaffine_displacement.rst b/doc/src/fix_nonaffine_displacement.rst index 466a7b1ab0..99d98b451b 100644 --- a/doc/src/fix_nonaffine_displacement.rst +++ b/doc/src/fix_nonaffine_displacement.rst @@ -70,7 +70,9 @@ displacement. A deformation gradient tensor is then calculated as D^2_\mathrm{min} = \sum_{\mathrm{neighbors}} \left| \vec{r} - F \vec{r}_0 \right|^2 and a strain tensor is calculated :math:`E = F F^{T} - I` where :math:`I` -is the identity tensor. +is the identity tensor. This calculation is only performed on timesteps that +are a multiple of *nevery* (including timestep zero). Data accessed before +this occurs will simply be zeroed. The *integrated* style simply integrates the velocity of particles every timestep to calculate a displacement. This style only works if @@ -78,6 +80,14 @@ used in conjunction with another fix that deforms the box and displaces atom positions such as :doc:`the remap x option of fix deform `, :doc:`fix press/berendsen `, or :doc:`fix nh `. +Both of these methods require defining a reference state. With the *fixed* reference +style, the user picks a specific timestep *nstep* from which particle positions are saved. +If peratom data is accessed from this compute prior to this timestep, it will simply be +zeroed. The *update* reference style implies the reference state will be updated every +*nstep* timesteps. The *offset* reference only applies to the *d2min* metric and will +update the reference state *nstep* timesteps before a multiple of *nevery* timesteps. + + ---------- Restart, fix_modify, output, run start/stop, minimize info diff --git a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp index f18fdd2c4a..bb5c2ca0c3 100755 --- a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp +++ b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp @@ -116,6 +116,9 @@ FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char * if (cut_style == RADIUS && (!atom->radius_flag)) error->all(FLERR, "Fix nonaffine/displacement radius style requires atom attribute radius"); + if (nad_style == INTEGRATED && reference_style == OFFSET) + error->all(FLERR, "Fix nonaffine/displacement cannot use the integrated style with an offset reference state"); + peratom_flag = 1; peratom_freq = nevery; nmax = -1; @@ -347,11 +350,6 @@ void FixNonaffineDisplacement::save_reference_state() void FixNonaffineDisplacement::calculate_D2Min() { - if (!reference_saved) { - error->warning(FLERR, "Calculating D2Min without a saved reference state"); - return; - } - // invoke half neighbor list (will copy or build if necessary) neighbor->build_one(list); From 1ffb2289be15a40757f39855f9d8a0301bbcb04b Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 17 Jul 2023 11:59:50 -0600 Subject: [PATCH 10/18] explicitly zeroing arrays --- doc/src/fix_nonaffine_displacement.rst | 2 +- src/EXTRA-FIX/fix_nonaffine_displacement.cpp | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/src/fix_nonaffine_displacement.rst b/doc/src/fix_nonaffine_displacement.rst index 99d98b451b..85f51b2419 100644 --- a/doc/src/fix_nonaffine_displacement.rst +++ b/doc/src/fix_nonaffine_displacement.rst @@ -81,7 +81,7 @@ atom positions such as :doc:`the remap x option of fix deform `, :doc:`fix press/berendsen `, or :doc:`fix nh `. Both of these methods require defining a reference state. With the *fixed* reference -style, the user picks a specific timestep *nstep* from which particle positions are saved. +style, the user picks a specific timestep *nstep* at which particle positions are saved. If peratom data is accessed from this compute prior to this timestep, it will simply be zeroed. The *update* reference style implies the reference state will be updated every *nstep* timesteps. The *offset* reference only applies to the *d2min* metric and will diff --git a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp index bb5c2ca0c3..5572aab1bf 100755 --- a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp +++ b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp @@ -180,11 +180,11 @@ void FixNonaffineDisplacement::post_constructor() array_atom = nad; } - if (nad_style == D2MIN) { + if (nad_style == D2MIN) grow_arrays(atom->nmax); - for (int i = 0; i < atom->nlocal; i++) - for (int j = 0; j < 3; j++) array_atom[i][j] = 0.0; - } + + for (int i = 0; i < atom->nlocal; i++) + for (int j = 0; j < 3; j++) array_atom[i][j] = 0.0; } /* ---------------------------------------------------------------------- */ From 2d897ac8f3eda1352210721c83d3640e3dff3686 Mon Sep 17 00:00:00 2001 From: Joel Thomas Clemmer Date: Mon, 17 Jul 2023 14:22:22 -0600 Subject: [PATCH 11/18] Fixing doc build errors --- doc/src/Commands_compute.rst | 1 + doc/src/Commands_fix.rst | 1 + doc/src/compute_rattlers.rst | 2 +- doc/src/fix_nonaffine_displacement.rst | 8 ++++---- 4 files changed, 7 insertions(+), 5 deletions(-) diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 755000c976..712391c574 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -114,6 +114,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`property/grid ` * :doc:`property/local ` * :doc:`ptm/atom ` + * :doc:`rattlers ` * :doc:`rdf ` * :doc:`reduce ` * :doc:`reduce/chunk ` diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index 6fe321e3c9..3ab4533c80 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -123,6 +123,7 @@ OPT. * :doc:`mvv/tdpd ` * :doc:`neb ` * :doc:`neb/spin ` + * :doc:`nonaffine/displacement ` * :doc:`nph (ko) ` * :doc:`nph/asphere (o) ` * :doc:`nph/body ` diff --git a/doc/src/compute_rattlers.rst b/doc/src/compute_rattlers.rst index a0f0d1b54a..f873a65dbf 100644 --- a/doc/src/compute_rattlers.rst +++ b/doc/src/compute_rattlers.rst @@ -1,7 +1,7 @@ .. index:: compute rattlers compute rattlers command -============================ +======================== Syntax """""" diff --git a/doc/src/fix_nonaffine_displacement.rst b/doc/src/fix_nonaffine_displacement.rst index 85f51b2419..7764b64066 100644 --- a/doc/src/fix_nonaffine_displacement.rst +++ b/doc/src/fix_nonaffine_displacement.rst @@ -1,7 +1,7 @@ .. index:: fix nonaffine/displacement fix nonaffine/displacement command -=================== +================================== Syntax """""" @@ -77,7 +77,7 @@ this occurs will simply be zeroed. The *integrated* style simply integrates the velocity of particles every timestep to calculate a displacement. This style only works if used in conjunction with another fix that deforms the box and displaces -atom positions such as :doc:`the remap x option of fix deform `, +atom positions such as :doc:`fix deform ` with remap x, :doc:`fix press/berendsen `, or :doc:`fix nh `. Both of these methods require defining a reference state. With the *fixed* reference @@ -126,6 +126,6 @@ none ---------- -.. _nh-Martyna: +.. _d2min-Falk: -**(Martyna)** Martyna, Tobias and Klein, J Chem Phys, 101, 4177 (1994). \ No newline at end of file +**(Falk)** Falk and Langer PRE, 57, 7192 (1998). \ No newline at end of file From e9223fc5af8bb9b0b5dd8fc945c40dfa73417f5c Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 17 Jul 2023 14:34:46 -0600 Subject: [PATCH 12/18] Fixing LAMMPS headers --- src/EXTRA-COMPUTE/compute_rattlers.cpp | 3 ++- src/EXTRA-COMPUTE/compute_rattlers.h | 2 +- src/EXTRA-FIX/fix_nonaffine_displacement.cpp | 4 ++-- src/EXTRA-FIX/fix_nonaffine_displacement.h | 4 ++-- 4 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_rattlers.cpp b/src/EXTRA-COMPUTE/compute_rattlers.cpp index 5613afb707..5de246f990 100755 --- a/src/EXTRA-COMPUTE/compute_rattlers.cpp +++ b/src/EXTRA-COMPUTE/compute_rattlers.cpp @@ -1,7 +1,8 @@ +// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov + LAMMPS development team: developers@lammps.org Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains diff --git a/src/EXTRA-COMPUTE/compute_rattlers.h b/src/EXTRA-COMPUTE/compute_rattlers.h index 08f8d7e71d..5ed0cd5882 100755 --- a/src/EXTRA-COMPUTE/compute_rattlers.h +++ b/src/EXTRA-COMPUTE/compute_rattlers.h @@ -1,7 +1,7 @@ /* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov + LAMMPS development team: developers@lammps.org Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains diff --git a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp index 5572aab1bf..704ac93382 100755 --- a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp +++ b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp @@ -1,8 +1,8 @@ // clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains diff --git a/src/EXTRA-FIX/fix_nonaffine_displacement.h b/src/EXTRA-FIX/fix_nonaffine_displacement.h index 9a914cf5df..5bbc335892 100755 --- a/src/EXTRA-FIX/fix_nonaffine_displacement.h +++ b/src/EXTRA-FIX/fix_nonaffine_displacement.h @@ -1,7 +1,7 @@ /* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains From bb2c286f2728391f8d8e35d45592a54c47f51173 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 18 Aug 2023 06:46:10 -0400 Subject: [PATCH 13/18] fix whitespace --- doc/src/fix_nonaffine_displacement.rst | 2 +- src/GRANULAR/gran_sub_mod_normal.h | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/src/fix_nonaffine_displacement.rst b/doc/src/fix_nonaffine_displacement.rst index 7764b64066..1936f44bde 100644 --- a/doc/src/fix_nonaffine_displacement.rst +++ b/doc/src/fix_nonaffine_displacement.rst @@ -128,4 +128,4 @@ none .. _d2min-Falk: -**(Falk)** Falk and Langer PRE, 57, 7192 (1998). \ No newline at end of file +**(Falk)** Falk and Langer PRE, 57, 7192 (1998). diff --git a/src/GRANULAR/gran_sub_mod_normal.h b/src/GRANULAR/gran_sub_mod_normal.h index c1511fdfa5..6f2f3aabbe 100644 --- a/src/GRANULAR/gran_sub_mod_normal.h +++ b/src/GRANULAR/gran_sub_mod_normal.h @@ -94,7 +94,7 @@ namespace Granular_NS { void coeffs_to_local() override; void mix_coeffs(double *, double *) override; private: - int mixed_coefficients; + int mixed_coefficients; }; /* ---------------------------------------------------------------------- */ @@ -110,7 +110,7 @@ namespace Granular_NS { protected: double k, cohesion; double F_pulloff, Fne; - int mixed_coefficients; + int mixed_coefficients; }; /* ---------------------------------------------------------------------- */ @@ -129,7 +129,7 @@ namespace Granular_NS { protected: double k, cohesion; double Emix, F_pulloff, Fne; - int mixed_coefficients; + int mixed_coefficients; }; } // namespace Granular_NS From 67f42fa84f723dd6cc08f6269f3ccd0537ff2588 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 18 Aug 2023 06:52:21 -0400 Subject: [PATCH 14/18] add version tags, sort out spelling issues --- doc/src/compute.rst | 2 +- doc/src/compute_rattlers.rst | 4 +++- doc/src/fix.rst | 2 +- doc/src/fix_nonaffine_displacement.rst | 2 ++ doc/utils/sphinx-config/false_positives.txt | 4 ++++ 5 files changed, 11 insertions(+), 3 deletions(-) diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 5dcef4dd0e..6667f9c2b4 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -268,7 +268,7 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`property/grid ` - convert per-grid attributes to per-grid vectors/arrays * :doc:`property/local ` - convert local attributes to local vectors/arrays * :doc:`ptm/atom ` - determines the local lattice structure based on the Polyhedral Template Matching method -* :doc:`rattlers ` - identify undercoordinated rattler atoms +* :doc:`rattlers ` - identify under-coordinated rattler atoms * :doc:`rdf ` - radial distribution function :math:`g(r)` histogram of group of atoms * :doc:`reduce ` - combine per-atom quantities into a single global value * :doc:`reduce/chunk ` - reduce per-atom quantities within each chunk diff --git a/doc/src/compute_rattlers.rst b/doc/src/compute_rattlers.rst index f873a65dbf..58dba3a126 100644 --- a/doc/src/compute_rattlers.rst +++ b/doc/src/compute_rattlers.rst @@ -32,8 +32,10 @@ Examples Description """"""""""" +.. versionadded:: TBD + Define a compute that identifies rattlers in a system. Rattlers are -identified using an interative approach. The coordination number of +identified using an interactive approach. The coordination number of all atoms is first calculated. The *type* and *radius* settings are used to select whether interaction cutoffs are determined by atom types or by the sum of atomic radii (atom style sphere), respectively. diff --git a/doc/src/fix.rst b/doc/src/fix.rst index f05c7927c1..23346ccf67 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -275,7 +275,7 @@ accelerated styles exist. * :doc:`mvv/tdpd ` - constant temperature DPD using the modified velocity-Verlet algorithm * :doc:`neb ` - nudged elastic band (NEB) spring forces * :doc:`neb/spin ` - nudged elastic band (NEB) spring forces for spins -* :doc:`nonaffine/displacement ` - calculate nonaffined displacement of atoms +* :doc:`nonaffine/displacement ` - calculate nonaffine displacement of atoms * :doc:`nph ` - constant NPH time integration via Nose/Hoover * :doc:`nph/asphere ` - NPH for aspherical particles * :doc:`nph/body ` - NPH for body particles diff --git a/doc/src/fix_nonaffine_displacement.rst b/doc/src/fix_nonaffine_displacement.rst index 1936f44bde..363b0a747a 100644 --- a/doc/src/fix_nonaffine_displacement.rst +++ b/doc/src/fix_nonaffine_displacement.rst @@ -44,6 +44,8 @@ Examples Description """"""""""" +.. versionadded:: TBD + This fix computes different metrics of the nonaffine displacement of particles. The first metric, *d2min* calculates the :math:`D^2_\mathrm{min}` nonaffine displacement by Falk and Langer in :ref:`(Falk) `. diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index c7b73c9f14..82901ad649 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -1073,6 +1073,7 @@ facesets factorizable factorizations Fahrenberger +Falk Faken Farago Fasolino @@ -1816,6 +1817,7 @@ Lanczos Lande Landron Landsgesell +Langer langevin Langevin Langston @@ -2491,6 +2493,7 @@ noforce noguess Noid nolib +nonaffine nonequilibrium nongauss nonGaussian @@ -2563,6 +2566,7 @@ nthreads ntimestep Ntptask Ntriples +ntries ntris Ntype ntypes From 9702a7a9d437d1cf59ef5e506fb157c3609b6b4c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 18 Aug 2023 06:56:14 -0400 Subject: [PATCH 15/18] don't need to list utils.h as include as it is always included for styles --- src/EXTRA-COMPUTE/compute_rattlers.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_rattlers.cpp b/src/EXTRA-COMPUTE/compute_rattlers.cpp index 5de246f990..7a3771cdf5 100755 --- a/src/EXTRA-COMPUTE/compute_rattlers.cpp +++ b/src/EXTRA-COMPUTE/compute_rattlers.cpp @@ -28,7 +28,6 @@ #include "neighbor.h" #include "pair.h" #include "update.h" -#include "utils.h" #include #include From 9ab79c745b631403351de442be523fc1a4cf765c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 18 Aug 2023 07:03:02 -0400 Subject: [PATCH 16/18] no executable permissions on source code --- src/EXTRA-COMPUTE/compute_rattlers.cpp | 0 src/EXTRA-COMPUTE/compute_rattlers.h | 0 src/EXTRA-FIX/fix_nonaffine_displacement.cpp | 0 src/EXTRA-FIX/fix_nonaffine_displacement.h | 0 4 files changed, 0 insertions(+), 0 deletions(-) mode change 100755 => 100644 src/EXTRA-COMPUTE/compute_rattlers.cpp mode change 100755 => 100644 src/EXTRA-COMPUTE/compute_rattlers.h mode change 100755 => 100644 src/EXTRA-FIX/fix_nonaffine_displacement.cpp mode change 100755 => 100644 src/EXTRA-FIX/fix_nonaffine_displacement.h diff --git a/src/EXTRA-COMPUTE/compute_rattlers.cpp b/src/EXTRA-COMPUTE/compute_rattlers.cpp old mode 100755 new mode 100644 diff --git a/src/EXTRA-COMPUTE/compute_rattlers.h b/src/EXTRA-COMPUTE/compute_rattlers.h old mode 100755 new mode 100644 diff --git a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp old mode 100755 new mode 100644 diff --git a/src/EXTRA-FIX/fix_nonaffine_displacement.h b/src/EXTRA-FIX/fix_nonaffine_displacement.h old mode 100755 new mode 100644 From c5deb581c28c46778c241772a12412914bc40780 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 21 Sep 2023 22:12:58 +0200 Subject: [PATCH 17/18] Updates to address PR comments --- src/.gitignore | 4 +-- ...rattlers.cpp => compute_rattlers_atom.cpp} | 28 +++++++-------- ...ute_rattlers.h => compute_rattlers_atom.h} | 12 +++---- src/EXTRA-FIX/fix_nonaffine_displacement.cpp | 35 ++++++++----------- src/EXTRA-FIX/fix_nonaffine_displacement.h | 5 +-- 5 files changed, 39 insertions(+), 45 deletions(-) rename src/EXTRA-COMPUTE/{compute_rattlers.cpp => compute_rattlers_atom.cpp} (89%) rename src/EXTRA-COMPUTE/{compute_rattlers.h => compute_rattlers_atom.h} (84%) diff --git a/src/.gitignore b/src/.gitignore index a170061ea7..955e740511 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -602,8 +602,8 @@ /compute_pressure_grem.h /compute_ptm_atom.cpp /compute_ptm_atom.h -/compute_rattlers.cpp -/compute_rattlers.h +/compute_rattlers_atom.cpp +/compute_rattlers_atom.h /compute_rigid_local.cpp /compute_rigid_local.h /compute_smd_triangle_vertices.cpp diff --git a/src/EXTRA-COMPUTE/compute_rattlers.cpp b/src/EXTRA-COMPUTE/compute_rattlers_atom.cpp similarity index 89% rename from src/EXTRA-COMPUTE/compute_rattlers.cpp rename to src/EXTRA-COMPUTE/compute_rattlers_atom.cpp index 7a3771cdf5..602923b58a 100644 --- a/src/EXTRA-COMPUTE/compute_rattlers.cpp +++ b/src/EXTRA-COMPUTE/compute_rattlers_atom.cpp @@ -16,7 +16,7 @@ Contributing authors: Joel Clemmer (SNL), Ishan Srivastava (LBNL) ------------------------------------------------------------------------- */ -#include "compute_rattlers.h" +#include "compute_rattlers_atom.h" #include "atom.h" #include "comm.h" @@ -38,20 +38,20 @@ enum { TYPE, RADIUS }; /* ---------------------------------------------------------------------- */ -ComputeRattlers::ComputeRattlers(LAMMPS *lmp, int narg, char **arg) : +ComputeRattlersAtom::ComputeRattlersAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), ncontacts(nullptr), rattler(nullptr) { - if (narg != 6) error->all(FLERR, "Illegal compute fabric command"); + if (narg != 6) error->all(FLERR, "Illegal compute rattlers/atom command"); if (strcmp(arg[3], "type") == 0) cutstyle = TYPE; else if (strcmp(arg[3], "radius") == 0) cutstyle = RADIUS; else - error->all(FLERR, "Illegal compute fabric command"); + error->all(FLERR, "Illegal compute rattlers/atom command"); if (cutstyle == RADIUS && !atom->radius_flag) - error->all(FLERR, "Compute fabric radius style requires atom attribute radius"); + error->all(FLERR, "Compute rattlers/atom radius style requires atom attribute radius"); ncontacts_rattler = utils::inumeric(FLERR, arg[4], false, lmp); max_tries = utils::inumeric(FLERR, arg[5], false, lmp); @@ -69,7 +69,7 @@ ComputeRattlers::ComputeRattlers(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ -ComputeRattlers::~ComputeRattlers() +ComputeRattlersAtom::~ComputeRattlersAtom() { memory->destroy(ncontacts); memory->destroy(rattler); @@ -77,7 +77,7 @@ ComputeRattlers::~ComputeRattlers() /* ---------------------------------------------------------------------- */ -void ComputeRattlers::init() +void ComputeRattlersAtom::init() { if (force->pair == nullptr) error->all(FLERR, "No pair style is defined for compute rattlers"); @@ -98,14 +98,14 @@ void ComputeRattlers::init() /* ---------------------------------------------------------------------- */ -void ComputeRattlers::init_list(int /*id*/, NeighList *ptr) +void ComputeRattlersAtom::init_list(int /*id*/, NeighList *ptr) { list = ptr; } /* ---------------------------------------------------------------------- */ -void ComputeRattlers::compute_peratom() +void ComputeRattlersAtom::compute_peratom() { if (invoked_peratom == update->ntimestep) return; invoked_peratom = update->ntimestep; @@ -236,7 +236,7 @@ void ComputeRattlers::compute_peratom() /* ---------------------------------------------------------------------- */ -double ComputeRattlers::compute_scalar() +double ComputeRattlersAtom::compute_scalar() { if (invoked_peratom != update->ntimestep) compute_peratom(); @@ -257,7 +257,7 @@ double ComputeRattlers::compute_scalar() /* ---------------------------------------------------------------------- */ -int ComputeRattlers::pack_reverse_comm(int n, int first, double *buf) +int ComputeRattlersAtom::pack_reverse_comm(int n, int first, double *buf) { int i, m, last; @@ -271,7 +271,7 @@ int ComputeRattlers::pack_reverse_comm(int n, int first, double *buf) /* ---------------------------------------------------------------------- */ -void ComputeRattlers::unpack_reverse_comm(int n, int *list, double *buf) +void ComputeRattlersAtom::unpack_reverse_comm(int n, int *list, double *buf) { int i, j, m; @@ -284,7 +284,7 @@ void ComputeRattlers::unpack_reverse_comm(int n, int *list, double *buf) /* ---------------------------------------------------------------------- */ -int ComputeRattlers::pack_forward_comm(int n, int *list, double *buf, +int ComputeRattlersAtom::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) { int i, j, m; @@ -300,7 +300,7 @@ int ComputeRattlers::pack_forward_comm(int n, int *list, double *buf, /* ---------------------------------------------------------------------- */ -void ComputeRattlers::unpack_forward_comm(int n, int first, double *buf) +void ComputeRattlersAtom::unpack_forward_comm(int n, int first, double *buf) { int i, m, last; diff --git a/src/EXTRA-COMPUTE/compute_rattlers.h b/src/EXTRA-COMPUTE/compute_rattlers_atom.h similarity index 84% rename from src/EXTRA-COMPUTE/compute_rattlers.h rename to src/EXTRA-COMPUTE/compute_rattlers_atom.h index 5ed0cd5882..257bae8374 100644 --- a/src/EXTRA-COMPUTE/compute_rattlers.h +++ b/src/EXTRA-COMPUTE/compute_rattlers_atom.h @@ -13,21 +13,21 @@ #ifdef COMPUTE_CLASS // clang-format off -ComputeStyle(rattlers,ComputeRattlers); +ComputeStyle(rattlers/atom,ComputeRattlersAtom); // clang-format on #else -#ifndef LMP_COMPUTE_RATTLERS_H -#define LMP_COMPUTE_RATTLERS_H +#ifndef LMP_COMPUTE_RATTLERS_ATOM_H +#define LMP_COMPUTE_RATTLERS_ATOM_H #include "compute.h" namespace LAMMPS_NS { -class ComputeRattlers : public Compute { +class ComputeRattlersAtom : public Compute { public: - ComputeRattlers(class LAMMPS *, int, char **); - ~ComputeRattlers() override; + ComputeRattlersAtom(class LAMMPS *, int, char **); + ~ComputeRattlersAtom() override; void init() override; void init_list(int, class NeighList *) override; void compute_peratom() override; diff --git a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp index 704ac93382..ef5481601f 100644 --- a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp +++ b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp @@ -23,6 +23,7 @@ #include "comm.h" #include "domain.h" #include "error.h" +#include "fix_store_atom.h" #include "force.h" #include "group.h" #include "math_extra.h" @@ -63,7 +64,7 @@ static const char cite_nonaffine_d2min[] = /* ---------------------------------------------------------------------- */ FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), new_fix_id(nullptr), X(nullptr), Y(nullptr), F(nullptr), norm(nullptr) + Fix(lmp, narg, arg), id_fix(nullptr), X(nullptr), Y(nullptr), F(nullptr), norm(nullptr) { if (narg < 4) error->all(FLERR,"Illegal fix nonaffine/displacement command"); @@ -140,8 +141,8 @@ FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char * FixNonaffineDisplacement::~FixNonaffineDisplacement() { - if (new_fix_id && modify->nfix) modify->delete_fix(new_fix_id); - delete[] new_fix_id; + if (id_fix && modify->nfix) modify->delete_fix(id_fix); + delete[] id_fix; if (nad_style == D2MIN) { memory->destroy(X); @@ -167,18 +168,14 @@ void FixNonaffineDisplacement::post_constructor() { // Create persistent peratom storage for either an integrated velocity or reference position // Ghost atoms need reference coordinates for D2min - std::string ghost_status = "no"; - if (nad_style == D2MIN) ghost_status = "yes"; + std::string ghost_status = "0"; + if (nad_style == D2MIN) ghost_status = "1"; - new_fix_id = utils::strdup(id + std::string("_FIX_PA")); - modify->add_fix(fmt::format("{} {} property/atom d2_nad 3 ghost {}", new_fix_id, group->names[igroup], ghost_status)); - int tmp1, tmp2; - nad_index = atom->find_custom("nad", tmp1, tmp2); + id_fix = utils::strdup(id + std::string("_FIX_PA")); + fix = dynamic_cast(modify->add_fix(fmt::format("{} {} STORE/ATOM 3 0 {} 1", id_fix, group->names[igroup], ghost_status))); - if (nad_style == INTEGRATED) { - double **nad = atom->darray[nad_index]; - array_atom = nad; - } + if (nad_style == INTEGRATED) + array_atom = fix->astore; if (nad_style == D2MIN) grow_arrays(atom->nmax); @@ -294,12 +291,11 @@ void FixNonaffineDisplacement::integrate_velocity() int *mask = atom->mask; int nlocal = atom->nlocal; - double **nad = atom->darray[nad_index]; for (int m = 0; m < 3; m++) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - nad[i][m] += dtv * v[i][m]; + array_atom[i][m] += dtv * v[i][m]; } } } @@ -315,18 +311,17 @@ void FixNonaffineDisplacement::save_reference_state() int *mask = atom->mask; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; - double **nad = atom->darray[nad_index]; if (nad_style == D2MIN) { for (int m = 0; m < 3; m++) { for (int i = 0; i < nall; i++) { - if (mask[i] & groupbit) nad[i][m] = x[i][m]; + if (mask[i] & groupbit) array_atom[i][m] = x[i][m]; } } } else { for (int m = 0; m < 3; m++) { for (int i = 0; i < nall; i++) { - if (mask[i] & groupbit) nad[i][m] = 0.0; + if (mask[i] & groupbit) array_atom[i][m] = 0.0; } } } @@ -364,7 +359,7 @@ void FixNonaffineDisplacement::calculate_D2Min() int *ilist, *jlist, *numneigh, **firstneigh; double **x = atom->x; - double **x0 = atom->darray[nad_index]; + double **x0 = array_atom; double *radius = atom->radius; tagint *tag = atom->tag; int *type = atom->type; @@ -730,12 +725,10 @@ void FixNonaffineDisplacement::minimum_image0(double *delta) void FixNonaffineDisplacement::grow_arrays(int nmax_new) { nmax = nmax_new; - memory->destroy(array_atom); memory->destroy(X); memory->destroy(Y); memory->destroy(F); memory->destroy(norm); - memory->create(array_atom, nmax, 3, "fix_nonaffine_displacement:array_atom"); memory->create(X, nmax, 3, 3, "fix_nonaffine_displacement:X"); memory->create(Y, nmax, 3, 3, "fix_nonaffine_displacement:Y"); memory->create(F, nmax, 3, 3, "fix_nonaffine_displacement:F"); diff --git a/src/EXTRA-FIX/fix_nonaffine_displacement.h b/src/EXTRA-FIX/fix_nonaffine_displacement.h index 5bbc335892..0a195dc08e 100644 --- a/src/EXTRA-FIX/fix_nonaffine_displacement.h +++ b/src/EXTRA-FIX/fix_nonaffine_displacement.h @@ -43,8 +43,9 @@ class FixNonaffineDisplacement : public Fix { private: double dtv; - char *new_fix_id; - int nad_index, nmax, comm_flag; + char *id_fix; + class FixStoreAtom *fix; + int nmax, comm_flag; int nad_style, cut_style; int reference_style, offset_timestep, reference_timestep, update_timestep; int reference_saved; From c22d6a9e4e49a1a632194fafc6df46d1d04ee5d8 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 21 Sep 2023 22:13:16 +0200 Subject: [PATCH 18/18] update to doc pages --- doc/src/Commands_compute.rst | 2 +- doc/src/compute.rst | 2 +- doc/src/compute_rattlers.rst | 77 -------------------------- doc/src/compute_rattlers_atom.rst | 92 +++++++++++++++++++++++++++++++ 4 files changed, 94 insertions(+), 79 deletions(-) delete mode 100644 doc/src/compute_rattlers.rst create mode 100644 doc/src/compute_rattlers_atom.rst diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 01d3c5c161..07f5100ea4 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -115,7 +115,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`property/grid ` * :doc:`property/local ` * :doc:`ptm/atom ` - * :doc:`rattlers ` + * :doc:`rattlers/atom ` * :doc:`rdf ` * :doc:`reduce ` * :doc:`reduce/chunk ` diff --git a/doc/src/compute.rst b/doc/src/compute.rst index d63bb1abb1..8471f71fa8 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -269,7 +269,7 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`property/grid ` - convert per-grid attributes to per-grid vectors/arrays * :doc:`property/local ` - convert local attributes to local vectors/arrays * :doc:`ptm/atom ` - determines the local lattice structure based on the Polyhedral Template Matching method -* :doc:`rattlers ` - identify under-coordinated rattler atoms +* :doc:`rattlers/atom ` - identify under-coordinated rattler atoms * :doc:`rdf ` - radial distribution function :math:`g(r)` histogram of group of atoms * :doc:`reduce ` - combine per-atom quantities into a single global value * :doc:`reduce/chunk ` - reduce per-atom quantities within each chunk diff --git a/doc/src/compute_rattlers.rst b/doc/src/compute_rattlers.rst deleted file mode 100644 index 58dba3a126..0000000000 --- a/doc/src/compute_rattlers.rst +++ /dev/null @@ -1,77 +0,0 @@ -.. index:: compute rattlers - -compute rattlers command -======================== - -Syntax -"""""" - -.. parsed-literal:: - - compute ID group-ID rattlers cutoff zmin ntries - -* ID, group-ID are documented in :doc:`compute ` command -* rattlers = style name of this compute command -* cutoff = *type* or *radius* - - .. parsed-literal:: - - *type* = cutoffs determined based on atom types - *radius* = cutoffs determined based on atom diameters (atom style sphere) - -* zmin = minimum coordination for a non-rattler particle -* ntries = maximum number of iterations to remove rattlers - -Examples -"""""""" - -.. code-block:: LAMMPS - - compute 1 all rattlers type 4 10 - -Description -""""""""""" - -.. versionadded:: TBD - -Define a compute that identifies rattlers in a system. Rattlers are -identified using an interactive approach. The coordination number of -all atoms is first calculated. The *type* and *radius* settings are -used to select whether interaction cutoffs are determined by atom -types or by the sum of atomic radii (atom style sphere), respectively. -Rattlers are then identified as particles with a coordination number -less than *zmin* and are removed from consideration. Atomic coordination -numbers are then recalculated, excluding previously identified rattlers, -to identify a new set of rattlers. This process is iterated, up to a maximum -of *ntries*, until no new rattlers are identified and the remaining -atoms form a stable network of contacts. - -Output info -""""""""""" - -This compute calculates a per-atom vector and a global scalar. The vector -designates which atoms are rattlers, indicated by a value 1. Non-rattlers -have a value of 0. The global scalar returns the total number of rattlers -in the system. See the :doc:`Howto output ` page for an -overview of LAMMPS output options. - -Restrictions -"""""""""""" - -This compute is part of the EXTRA-COMPUTE package. It is only enabled if -LAMMPS was built with that package. See the -:doc:`Build package ` page for more info. - -The *radius* cutoff option requires that atoms store a radius as defined by the -:doc:`atom_style sphere ` or similar commands. - -Related commands -"""""""""""""""" - -:doc:`compute coord/atom ` -:doc:`compute contact/atom ` - -Default -""""""" - -none diff --git a/doc/src/compute_rattlers_atom.rst b/doc/src/compute_rattlers_atom.rst new file mode 100644 index 0000000000..a69d091466 --- /dev/null +++ b/doc/src/compute_rattlers_atom.rst @@ -0,0 +1,92 @@ +.. index:: compute rattlers/atom + +compute rattlers/atom command +======================== + +Syntax +"""""" + +.. parsed-literal:: + + compute ID group-ID rattlers/atom cutoff zmin ntries + +* ID, group-ID are documented in :doc:`compute ` command +* rattlers/atom = style name of this compute command +* cutoff = *type* or *radius* + + .. parsed-literal:: + + *type* = cutoffs determined based on atom types + *radius* = cutoffs determined based on atom diameters (atom style sphere) + +* zmin = minimum coordination for a non-rattler atom +* ntries = maximum number of iterations to remove rattlers + +Examples +"""""""" + +.. code-block:: LAMMPS + + compute 1 all rattlers/atom type 4 10 + +Description +""""""""""" + +.. versionadded:: TBD + +Define a compute that identifies rattlers in a system. Rattlers are often +identified in granular or glassy packings as undercoordinated atoms that +do not have the required number of contacts to constrain their translational +degrees of freedom. Such atoms are not considered rigid and can often freely +rattle around in the system. This compute identifies rattlers which can be +helpful for excluding them from analysis or providing extra damping forces +to accelerate relaxation processes. + +Rattlers are identified using an interactive approach. The coordination +number of all atoms is first calculated. The *type* and *radius* settings +are used to select whether interaction cutoffs are determined by atom +types or by the sum of atomic radii (atom style sphere), respectively. +Rattlers are then identified as atoms with a coordination number less +than *zmin* and are removed from consideration. Atomic coordination +numbers are then recalculated, excluding previously identified rattlers, +to identify a new set of rattlers. This process is iterated up to a maximum +of *ntries* or until no new rattlers are identified and the remaining +atoms form a stable network of contacts. + +In dense homogeneous systems where the average atom coordination number +is expected to be larger than *zmin*, this process usually only takes a few +iterations and a value of *ntries* around ten may be sufficient. In systems +with significant heterogeneity or average coordination numbers less than +*zmin*, an appropriate value of *ntries* depends heavily on the specific +system. For instance, a linear chain of N rattler atoms with a *zmin* of 2 +would take N/2 iterations to identify that all the atoms are rattlers. + +Output info +""""""""""" + +This compute calculates a per-atom vector and a global scalar. The vector +designates which atoms are rattlers, indicated by a value 1. Non-rattlers +have a value of 0. The global scalar returns the total number of rattlers +in the system. See the :doc:`Howto output ` page for an +overview of LAMMPS output options. + +Restrictions +"""""""""""" + +This compute is part of the EXTRA-COMPUTE package. It is only enabled if +LAMMPS was built with that package. See the +:doc:`Build package ` page for more info. + +The *radius* cutoff option requires that atoms store a radius as defined by the +:doc:`atom_style sphere ` or similar commands. + +Related commands +"""""""""""""""" + +:doc:`compute coord/atom ` +:doc:`compute contact/atom ` + +Default +""""""" + +none