diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index dbd6b58ce7..07f5100ea4 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -115,6 +115,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`property/grid ` * :doc:`property/local ` * :doc:`ptm/atom ` + * :doc:`rattlers/atom ` * :doc:`rdf ` * :doc:`reduce ` * :doc:`reduce/chunk ` diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index 7301d1345e..3283561c41 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -122,6 +122,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.rst b/doc/src/compute.rst index 6737203618..a99d7d7c9c 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -279,6 +279,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/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_contact_atom.rst b/doc/src/compute_contact_atom.rst index 31aa24aa60..b7ed062ff6 100644 --- a/doc/src/compute_contact_atom.rst +++ b/doc/src/compute_contact_atom.rst @@ -36,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 """"""""""" @@ -47,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 """""""""""" @@ -69,6 +69,3 @@ Default """"""" *group2-ID* = all - - -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 diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 0889fe281f..69a7212487 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -287,6 +287,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 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 new file mode 100644 index 0000000000..363b0a747a --- /dev/null +++ b/doc/src/fix_nonaffine_displacement.rst @@ -0,0 +1,133 @@ +.. index:: fix nonaffine/displacement + +fix nonaffine/displacement command +================================== + +Syntax +"""""" + +.. parsed-literal:: + + fix ID group nonaffine/displacement style args reference/style nstep + +* ID, group are documented in :doc:`fix ` command +* nonaffine/displacement = style name of this fix command +* nevery = calculate nonaffine displacement every this many timesteps +* style = *d2min* or *integrated* + + .. parsed-literal:: + + *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* (distance units) + *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 nonaffine displacement + +Examples +"""""""" + +.. code-block:: LAMMPS + + 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 +""""""""""" + +.. 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) `. +For each atom, the fix computes the two tensors + +.. math:: + + X = \sum_{\mathrm{neighbors}} \vec{r} \left(\vec{r}_{0} \right)^T + +and + +.. math:: + + Y = \sum_{\mathrm{neighbors}} \vec{r}_0 \left(\vec{r}_{0} \right)^T + +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:: + + 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. 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 +used in conjunction with another fix that deforms the box and displaces +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 +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 +update the reference state *nstep* timesteps before a multiple of *nevery* timesteps. + + +---------- + +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +The reference state is saved to :doc:`binary restart files `. + +None of the :doc:`fix_modify ` options are relevant to this +fix. + +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. + +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:`\sqrt{D^2_\mathrm{min}}`, the +volumetric strain, and the deviatoric strain. + +Restrictions +"""""""""""" + +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 +"""""""""""""""" + +none + +Default +""""""" + +none + +---------- + +.. _d2min-Falk: + +**(Falk)** Falk and Langer PRE, 57, 7192 (1998). diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 10ef06f89f..c31c6a0a70 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -1074,6 +1074,7 @@ facesets factorizable factorizations Fahrenberger +Falk Faken Farago Fasolino @@ -1822,6 +1823,7 @@ Lanczos Lande Landron Landsgesell +Langer langevin Langevin Langston @@ -2499,6 +2501,7 @@ noforce noguess Noid nolib +nonaffine nonequilibrium nongauss nonGaussian @@ -2571,6 +2574,7 @@ nthreads ntimestep Ntptask Ntriples +ntries ntris Ntype ntypes diff --git a/src/.gitignore b/src/.gitignore index 3ee771e139..e2c9fa6d5b 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -629,6 +629,8 @@ /compute_pressure_grem.h /compute_ptm_atom.cpp /compute_ptm_atom.h +/compute_rattlers_atom.cpp +/compute_rattlers_atom.h /compute_rigid_local.cpp /compute_rigid_local.h /compute_smd_triangle_vertices.cpp @@ -912,6 +914,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 diff --git a/src/EXTRA-COMPUTE/compute_rattlers_atom.cpp b/src/EXTRA-COMPUTE/compute_rattlers_atom.cpp new file mode 100644 index 0000000000..602923b58a --- /dev/null +++ b/src/EXTRA-COMPUTE/compute_rattlers_atom.cpp @@ -0,0 +1,312 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + 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 + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Joel Clemmer (SNL), Ishan Srivastava (LBNL) +------------------------------------------------------------------------- */ + +#include "compute_rattlers_atom.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 +#include + +using namespace LAMMPS_NS; + +enum { TYPE, RADIUS }; + +/* ---------------------------------------------------------------------- */ + +ComputeRattlersAtom::ComputeRattlersAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), ncontacts(nullptr), rattler(nullptr) +{ + 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 rattlers/atom command"); + + if (cutstyle == RADIUS && !atom->radius_flag) + 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); + + nmax = 0; + invoked_peratom = -1; + + scalar_flag = 1; + extscalar = 1; + peratom_flag = 1; + size_peratom_cols = 0; + comm_forward = 1; + comm_reverse = 1; +} + +/* ---------------------------------------------------------------------- */ + +ComputeRattlersAtom::~ComputeRattlersAtom() +{ + memory->destroy(ncontacts); + memory->destroy(rattler); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRattlersAtom::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 ComputeRattlersAtom::init_list(int /*id*/, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRattlersAtom::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 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) { + 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; + if (change_flag == 0) break; + + ntry += 1; + } + + if (change_flag == 1) + error->warning(FLERR, "Rattler calculation failed to converge within max tries"); +} + +/* ---------------------------------------------------------------------- */ + +double ComputeRattlersAtom::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 ComputeRattlersAtom::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 ComputeRattlersAtom::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 ComputeRattlersAtom::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 ComputeRattlersAtom::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/EXTRA-COMPUTE/compute_rattlers_atom.h b/src/EXTRA-COMPUTE/compute_rattlers_atom.h new file mode 100644 index 0000000000..257bae8374 --- /dev/null +++ b/src/EXTRA-COMPUTE/compute_rattlers_atom.h @@ -0,0 +1,52 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + 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 + 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/atom,ComputeRattlersAtom); +// clang-format on +#else + +#ifndef LMP_COMPUTE_RATTLERS_ATOM_H +#define LMP_COMPUTE_RATTLERS_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeRattlersAtom : public Compute { + public: + ComputeRattlersAtom(class LAMMPS *, int, char **); + ~ComputeRattlersAtom() 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/EXTRA-FIX/fix_nonaffine_displacement.cpp b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp new file mode 100644 index 0000000000..ef5481601f --- /dev/null +++ b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp @@ -0,0 +1,736 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + 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 + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Joel Clemmer (SNL), Ishan Srivastava (LBNL) +------------------------------------------------------------------------- */ + +#include "fix_nonaffine_displacement.h" + +#include "atom.h" +#include "citeme.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "fix_store_atom.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), id_fix(nullptr), X(nullptr), Y(nullptr), F(nullptr), norm(nullptr) +{ + if (narg < 4) error->all(FLERR,"Illegal fix nonaffine/displacement command"); + + 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 + 1 > narg) error->all(FLERR,"Illegal fix nonaffine/displacement command"); + nad_style = D2MIN; + if (strcmp(arg[iarg + 1], "type") == 0) { + cut_style = TYPE; + } else if (strcmp(arg[iarg + 1], "radius") == 0) { + cut_style = RADIUS; + } 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 + 2], false, lmp); + cutsq_custom = cutoff_custom * cutoff_custom; + if (cutoff_custom <= 0) + 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 += 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) { + 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) || (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]); + + if (nad_style == D2MIN) + 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; + 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 (id_fix && modify->nfix) modify->delete_fix(id_fix); + delete[] id_fix; + + 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 = "0"; + if (nad_style == D2MIN) ghost_status = "1"; + + 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) + array_atom = fix->astore; + + 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"); + + // 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) { + 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); + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +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_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(); + + 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(); +} + +/* ---------------------------------------------------------------------- */ + +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; + + for (int m = 0; m < 3; m++) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + array_atom[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; + int nall = nlocal + atom->nghost; + + if (nad_style == D2MIN) { + for (int m = 0; m < 3; m++) { + for (int i = 0; i < nall; i++) { + 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) array_atom[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() +{ + // 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 = array_atom; + 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; + int dim = domain->dimension; + + 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, 18); + + // 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 (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]; + 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; + array_atom[i][0] = sqrt(array_atom[i][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 < dim; j++) E[j][j] -= 1.0; + + evol = (E[0][0] + E[1][1] + E[2][2]) / dim; + + // Calculate deviatoric strain + 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++) + 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(X); + memory->destroy(Y); + memory->destroy(F); + memory->destroy(norm); + 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/EXTRA-FIX/fix_nonaffine_displacement.h b/src/EXTRA-FIX/fix_nonaffine_displacement.h new file mode 100644 index 0000000000..0a195dc08e --- /dev/null +++ b/src/EXTRA-FIX/fix_nonaffine_displacement.h @@ -0,0 +1,71 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + 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 + 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 *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; + 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 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"