From 25f5e74e9a7b240fe123f59c7749714d7304af46 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 15 Jun 2023 15:51:59 -0600 Subject: [PATCH 01/29] 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/29] 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/29] 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/29] 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/29] 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/29] 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/29] 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/29] 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/29] 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/29] 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/29] 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/29] 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/29] 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/29] 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/29] 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/29] 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/29] 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/29] 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 From 77549414aecc88e58aaf6a51b0299af70e7dd7e5 Mon Sep 17 00:00:00 2001 From: Yury Lysogorskiy Date: Wed, 22 Nov 2023 10:19:18 +0100 Subject: [PATCH 19/29] - implementation of ZBL-core repulsion with smooth transition from ACE to ZBL + Kokkos implementation - automatic BBasis(.yaml) to CTildeBasis(.yace) conversion for pair_pace --- cmake/Modules/Packages/ML-PACE.cmake | 53 +++++----- lib/pace/Install.py | 4 +- src/KOKKOS/pair_pace_extrapolation_kokkos.cpp | 87 ++++++++++++++++- src/KOKKOS/pair_pace_extrapolation_kokkos.h | 7 ++ src/KOKKOS/pair_pace_kokkos.cpp | 96 +++++++++++++++++-- src/KOKKOS/pair_pace_kokkos.h | 7 ++ src/ML-PACE/pair_pace.cpp | 12 ++- 7 files changed, 228 insertions(+), 38 deletions(-) diff --git a/cmake/Modules/Packages/ML-PACE.cmake b/cmake/Modules/Packages/ML-PACE.cmake index ce8f02f5f4..8a92dee395 100644 --- a/cmake/Modules/Packages/ML-PACE.cmake +++ b/cmake/Modules/Packages/ML-PACE.cmake @@ -1,33 +1,40 @@ -set(PACELIB_URL "https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/v.2023.10.04.tar.gz" CACHE STRING "URL for PACE evaluator library sources") +set(PACELIB_URL "https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/v.2023.11.22.tar.gz" CACHE STRING "URL for PACE evaluator library sources") -set(PACELIB_MD5 "70ff79f4e59af175e55d24f3243ad1ff" CACHE STRING "MD5 checksum of PACE evaluator library tarball") +set(PACELIB_MD5 "c8e811f96d761ef8863f5b88a3fd36f4" CACHE STRING "MD5 checksum of PACE evaluator library tarball") mark_as_advanced(PACELIB_URL) mark_as_advanced(PACELIB_MD5) GetFallbackURL(PACELIB_URL PACELIB_FALLBACK) -# download library sources to build folder -if(EXISTS ${CMAKE_BINARY_DIR}/libpace.tar.gz) - file(MD5 ${CMAKE_BINARY_DIR}/libpace.tar.gz DL_MD5) -endif() -if(NOT "${DL_MD5}" STREQUAL "${PACELIB_MD5}") - message(STATUS "Downloading ${PACELIB_URL}") - file(DOWNLOAD ${PACELIB_URL} ${CMAKE_BINARY_DIR}/libpace.tar.gz STATUS DL_STATUS SHOW_PROGRESS) - file(MD5 ${CMAKE_BINARY_DIR}/libpace.tar.gz DL_MD5) - if((NOT DL_STATUS EQUAL 0) OR (NOT "${DL_MD5}" STREQUAL "${PACELIB_MD5}")) - message(WARNING "Download from primary URL ${PACELIB_URL} failed\nTrying fallback URL ${PACELIB_FALLBACK}") - file(DOWNLOAD ${PACELIB_FALLBACK} ${CMAKE_BINARY_DIR}/libpace.tar.gz EXPECTED_HASH MD5=${PACELIB_MD5} SHOW_PROGRESS) - endif() +# LOCAL_ML-PACE points to top-level dir with local lammps-user-pace repo, +# to make it easier to check local build without going through the public github releases +if(LOCAL_ML-PACE) + set(lib-pace "${LOCAL_ML-PACE}") else() - message(STATUS "Using already downloaded archive ${CMAKE_BINARY_DIR}/libpace.tar.gz") -endif() + # download library sources to build folder + if(EXISTS ${CMAKE_BINARY_DIR}/libpace.tar.gz) + file(MD5 ${CMAKE_BINARY_DIR}/libpace.tar.gz DL_MD5) + endif() + if(NOT "${DL_MD5}" STREQUAL "${PACELIB_MD5}") + message(STATUS "Downloading ${PACELIB_URL}") + #file(DOWNLOAD ${PACELIB_URL} ${CMAKE_BINARY_DIR}/libpace.tar.gz STATUS DL_STATUS SHOW_PROGRESS) + #file(MD5 ${CMAKE_BINARY_DIR}/libpace.tar.gz DL_MD5) + if((NOT DL_STATUS EQUAL 0) OR (NOT "${DL_MD5}" STREQUAL "${PACELIB_MD5}")) + message(WARNING "Download from primary URL ${PACELIB_URL} failed\nTrying fallback URL ${PACELIB_FALLBACK}") + file(DOWNLOAD ${PACELIB_FALLBACK} ${CMAKE_BINARY_DIR}/libpace.tar.gz EXPECTED_HASH MD5=${PACELIB_MD5} SHOW_PROGRESS) + endif() + else() + message(STATUS "Using already downloaded archive ${CMAKE_BINARY_DIR}/libpace.tar.gz") + endif() -# uncompress downloaded sources -execute_process( - COMMAND ${CMAKE_COMMAND} -E remove_directory lammps-user-pace* - COMMAND ${CMAKE_COMMAND} -E tar xzf libpace.tar.gz - WORKING_DIRECTORY ${CMAKE_BINARY_DIR} -) -get_newest_file(${CMAKE_BINARY_DIR}/lammps-user-pace-* lib-pace) + + # uncompress downloaded sources + execute_process( + COMMAND ${CMAKE_COMMAND} -E remove_directory lammps-user-pace* + COMMAND ${CMAKE_COMMAND} -E tar xzf libpace.tar.gz + WORKING_DIRECTORY ${CMAKE_BINARY_DIR} + ) + get_newest_file(${CMAKE_BINARY_DIR}/lammps-user-pace-* lib-pace) +endif() add_subdirectory(${lib-pace} build-pace) set_target_properties(pace PROPERTIES CXX_EXTENSIONS ON OUTPUT_NAME lammps_pace${LAMMPS_MACHINE}) diff --git a/lib/pace/Install.py b/lib/pace/Install.py index 8d31852e44..8b89b3dec1 100644 --- a/lib/pace/Install.py +++ b/lib/pace/Install.py @@ -18,11 +18,11 @@ from install_helpers import fullpath, geturl, checkmd5sum, getfallback # settings thisdir = fullpath('.') -version ='v.2023.10.04' +version ='v.2023.11.22' # known checksums for different PACE versions. used to validate the download. checksums = { \ - 'v.2023.10.04': '70ff79f4e59af175e55d24f3243ad1ff' + 'v.2023.11.22': 'c8e811f96d761ef8863f5b88a3fd36f4' } parser = ArgumentParser(prog='Install.py', description="LAMMPS library build wrapper script") diff --git a/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp b/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp index 61722bf62d..b361cfe486 100644 --- a/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp +++ b/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp @@ -123,6 +123,9 @@ void PairPACEExtrapolationKokkos::grow(int natom, int maxneigh) // hard-core repulsion MemKK::realloc_kokkos(rho_core, "pace:rho_core", natom); MemKK::realloc_kokkos(dF_drho_core, "pace:dF_drho_core", natom); + MemKK::realloc_kokkos(dF_dfcut, "pace:dF_dfcut", natom); + MemKK::realloc_kokkos(d_d_min, "pace:r_min_pair", natom); + MemKK::realloc_kokkos(d_jj_min, "pace:j_min_pair", natom); MemKK::realloc_kokkos(dB_flatten, "pace:dB_flatten", natom, idx_ms_combs_max, basis_set->rankmax); @@ -219,6 +222,24 @@ void PairPACEExtrapolationKokkos::copy_pertype() Kokkos::deep_copy(d_wpre, h_wpre); Kokkos::deep_copy(d_mexp, h_mexp); + + + // ZBL core-rep + MemKK::realloc_kokkos(d_cut_in, "pace:d_cut_in", nelements, nelements); + MemKK::realloc_kokkos(d_dcut_in, "pace:d_dcut_in", nelements, nelements); + auto h_cut_in = Kokkos::create_mirror_view(d_cut_in); + auto h_dcut_in = Kokkos::create_mirror_view(d_dcut_in); + + for (int mu_i = 0; mu_i < nelements; ++mu_i) { + for (int mu_j = 0; mu_j < nelements; ++mu_j) { + h_cut_in(mu_i,mu_j) = basis_set->map_bond_specifications.at({mu_i,mu_j}).rcut_in; + h_dcut_in(mu_i,mu_j) = basis_set->map_bond_specifications.at({mu_i,mu_j}).dcut_in; + } + } + Kokkos::deep_copy(d_cut_in, h_cut_in); + Kokkos::deep_copy(d_dcut_in, h_dcut_in); + + is_zbl = basis_set->radial_functions->inner_cutoff_type == "zbl"; } /* ---------------------------------------------------------------------- */ @@ -631,6 +652,10 @@ void PairPACEExtrapolationKokkos::compute(int eflag_in, int vflag_in Kokkos::deep_copy(A_rank1, 0.0); Kokkos::deep_copy(rhos, 0.0); + Kokkos::deep_copy(rho_core, 0.0); + Kokkos::deep_copy(d_d_min, PairPACEExtrapolation::aceimpl->basis_set->cutoffmax); + Kokkos::deep_copy(d_jj_min, -1); + Kokkos::deep_copy(projections, 0.0); Kokkos::deep_copy(d_gamma, 0.0); @@ -799,6 +824,7 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeNeig const X_FLOAT ytmp = x(i,1); const X_FLOAT ztmp = x(i,2); const int jnum = d_numneigh[i]; + const int mu_i = d_map(type(i)); // get a pointer to scratch memory // This is used to cache whether or not an atom is within the cutoff @@ -858,6 +884,32 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeNeig } offset++; }); + + if(is_zbl) { + //adapted from https://www.osti.gov/servlets/purl/1429450 + using minloc_value_type=Kokkos::MinLoc::value_type; + minloc_value_type djjmin; + djjmin.val=1e20; + djjmin.loc=-1; + Kokkos::MinLoc reducer_scalar(djjmin); + // loop over ncount (actual neighbours withing cutoff) rather than jnum (total number of neigh in cutoff+skin) + Kokkos::parallel_reduce(Kokkos::TeamThreadRange(team, ncount), + [&](const int offset, minloc_value_type &min_d_dist) { + int j = d_nearest(ii,offset); + j &= NEIGHMASK; + const int jtype = type(j); + auto r = d_rnorms(ii,offset); + const int mu_j = d_map(type(j)); + const F_FLOAT d = r - (d_cut_in(mu_i, mu_j) - d_dcut_in(mu_i, mu_j)); + if (d < min_d_dist.val) { + min_d_dist.val = d; + min_d_dist.loc = offset; + } + + }, reducer_scalar); + d_d_min(ii) = djjmin.val; + d_jj_min(ii) = djjmin.loc;// d_jj_min should be NOT in 0..jnum range, but in 0..d_ncount(<=jnum) + } } /* ---------------------------------------------------------------------- */ @@ -1056,19 +1108,38 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeFS, const int ndensity = d_ndensity(mu_i); double evdwl, fcut, dfcut; + double evdwl_cut; evdwl = fcut = dfcut = 0.0; inner_cutoff(rho_core(ii), rho_cut, drho_cut, fcut, dfcut); FS_values_and_derivatives(ii, evdwl, mu_i); - dF_drho_core(ii) = evdwl * dfcut + 1; + if(is_zbl) { + if (d_jj_min(ii) != -1) { + const int mu_jmin = d_mu(ii,d_jj_min(ii)); + F_FLOAT dcutin = d_dcut_in(mu_i, mu_jmin); + F_FLOAT transition_coordinate = dcutin - d_d_min(ii); // == cutin - r_min + cutoff_func_poly(transition_coordinate, dcutin, dcutin, fcut, dfcut); + dfcut = -dfcut; // invert, because rho_core = cutin - r_min + } else { + // no neighbours + fcut = 1; + dfcut = 0; + } + evdwl_cut = evdwl * fcut + rho_core(ii) * (1 - fcut); // evdwl * fcut + rho_core_uncut - rho_core_uncut* fcut + dF_drho_core(ii) = 1 - fcut; + dF_dfcut(ii) = evdwl * dfcut - rho_core(ii) * dfcut; + } else { + inner_cutoff(rho_core(ii), rho_cut, drho_cut, fcut, dfcut); + dF_drho_core(ii) = evdwl * dfcut + 1; + evdwl_cut = evdwl * fcut + rho_core(ii); + } for (int p = 0; p < ndensity; ++p) dF_drho(ii, p) *= fcut; // tally energy contribution if (eflag) { - double evdwl_cut = evdwl * fcut + rho_core(ii); // E0 shift evdwl_cut += d_E0vals(mu_i); e_atom(ii) = evdwl_cut; @@ -1240,6 +1311,15 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeDeri f_ij(ii, jj, 0) = scale * f_ji[0] + fpair * r_hat[0]; f_ij(ii, jj, 1) = scale * f_ji[1] + fpair * r_hat[1]; f_ij(ii, jj, 2) = scale * f_ji[2] + fpair * r_hat[2]; + + if(is_zbl) { + if(jj==d_jj_min(ii)) { + // DCRU = 1.0 + f_ij(ii, jj, 0) += dF_dfcut(ii) * r_hat[0]; + f_ij(ii, jj, 1) += dF_dfcut(ii) * r_hat[1]; + f_ij(ii, jj, 2) += dF_dfcut(ii) * r_hat[2]; + } + } } /* ---------------------------------------------------------------------- */ @@ -1777,6 +1857,7 @@ double PairPACEExtrapolationKokkos::memory_usage() bytes += MemKK::memory_usage(weights_rank1); bytes += MemKK::memory_usage(rho_core); bytes += MemKK::memory_usage(dF_drho_core); + bytes += MemKK::memory_usage(dF_dfcut); bytes += MemKK::memory_usage(dB_flatten); bytes += MemKK::memory_usage(fr); bytes += MemKK::memory_usage(dfr); @@ -1794,6 +1875,8 @@ double PairPACEExtrapolationKokkos::memory_usage() bytes += MemKK::memory_usage(d_mu); bytes += MemKK::memory_usage(d_rhats); bytes += MemKK::memory_usage(d_rnorms); + bytes += MemKK::memory_usage(d_d_min); + bytes += MemKK::memory_usage(d_jj_min); bytes += MemKK::memory_usage(d_nearest); bytes += MemKK::memory_usage(f_ij); bytes += MemKK::memory_usage(d_rho_core_cutoff); diff --git a/src/KOKKOS/pair_pace_extrapolation_kokkos.h b/src/KOKKOS/pair_pace_extrapolation_kokkos.h index 55bcf4fead..092171cfb8 100644 --- a/src/KOKKOS/pair_pace_extrapolation_kokkos.h +++ b/src/KOKKOS/pair_pace_extrapolation_kokkos.h @@ -130,6 +130,7 @@ class PairPACEExtrapolationKokkos : public PairPACEExtrapolation { tdual_fparams k_cutsq, k_scale; typedef Kokkos::View t_fparams; t_fparams d_cutsq, d_scale; + t_fparams d_cut_in, d_dcut_in; // inner cutoff typename AT::t_int_1d d_map; @@ -240,6 +241,7 @@ class PairPACEExtrapolationKokkos : public PairPACEExtrapolation { t_ace_2d cr; t_ace_2d dcr; t_ace_1d dF_drho_core; + t_ace_1d dF_dfcut; // radial functions t_ace_4d fr; @@ -282,6 +284,11 @@ class PairPACEExtrapolationKokkos : public PairPACEExtrapolation { t_ace_3d3 d_rhats; t_ace_2i d_nearest; + // for ZBL core-rep implementation + t_ace_1d d_d_min; // [i] -> min-d for atom ii, d=d = r - (cut_in(mu_i, mu_j) - dcut_in(mu_i, mu_j)) + t_ace_1i d_jj_min; // [i] -> jj-index of nearest neigh (by r-(cut_in-dcut_in) criterion) + bool is_zbl; + // per-type t_ace_1i d_ndensity; t_ace_1i d_npoti; diff --git a/src/KOKKOS/pair_pace_kokkos.cpp b/src/KOKKOS/pair_pace_kokkos.cpp index 153a6d0333..62b7bf5d01 100644 --- a/src/KOKKOS/pair_pace_kokkos.cpp +++ b/src/KOKKOS/pair_pace_kokkos.cpp @@ -121,6 +121,9 @@ void PairPACEKokkos::grow(int natom, int maxneigh) // hard-core repulsion MemKK::realloc_kokkos(rho_core, "pace:rho_core", natom); MemKK::realloc_kokkos(dF_drho_core, "pace:dF_drho_core", natom); + MemKK::realloc_kokkos(dF_dfcut, "pace:dF_dfcut", natom); + MemKK::realloc_kokkos(d_d_min, "pace:r_min_pair", natom); + MemKK::realloc_kokkos(d_jj_min, "pace:j_min_pair", natom); MemKK::realloc_kokkos(dB_flatten, "pace:dB_flatten", natom, idx_rho_max, basis_set->rankmax); } @@ -212,6 +215,23 @@ void PairPACEKokkos::copy_pertype() Kokkos::deep_copy(d_wpre, h_wpre); Kokkos::deep_copy(d_mexp, h_mexp); + + // ZBL core-rep + MemKK::realloc_kokkos(d_cut_in, "pace:d_cut_in", nelements, nelements); + MemKK::realloc_kokkos(d_dcut_in, "pace:d_dcut_in", nelements, nelements); + auto h_cut_in = Kokkos::create_mirror_view(d_cut_in); + auto h_dcut_in = Kokkos::create_mirror_view(d_dcut_in); + + for (int mu_i = 0; mu_i < nelements; ++mu_i) { + for (int mu_j = 0; mu_j < nelements; ++mu_j) { + h_cut_in(mu_i,mu_j) = basis_set->map_bond_specifications.at({mu_i,mu_j}).rcut_in; + h_dcut_in(mu_i,mu_j) = basis_set->map_bond_specifications.at({mu_i,mu_j}).dcut_in; + } + } + Kokkos::deep_copy(d_cut_in, h_cut_in); + Kokkos::deep_copy(d_dcut_in, h_dcut_in); + + is_zbl = basis_set->radial_functions->inner_cutoff_type == "zbl"; } /* ---------------------------------------------------------------------- */ @@ -588,6 +608,8 @@ void PairPACEKokkos::compute(int eflag_in, int vflag_in) Kokkos::deep_copy(A_rank1, 0.0); Kokkos::deep_copy(rhos, 0.0); Kokkos::deep_copy(rho_core, 0.0); + Kokkos::deep_copy(d_d_min, PairPACE::aceimpl->basis_set->cutoffmax); + Kokkos::deep_copy(d_jj_min, -1); EV_FLOAT ev_tmp; @@ -665,6 +687,7 @@ void PairPACEKokkos::compute(int eflag_in, int vflag_in) Kokkos::parallel_for("ComputeDerivative",policy_derivative,*this); } + //ComputeForce { if (evflag) { @@ -741,6 +764,7 @@ void PairPACEKokkos::operator() (TagPairPACEComputeNeigh,const typen const X_FLOAT ytmp = x(i,1); const X_FLOAT ztmp = x(i,2); const int jnum = d_numneigh[i]; + const int mu_i = d_map(type(i)); // get a pointer to scratch memory // This is used to cache whether or not an atom is within the cutoff @@ -800,6 +824,32 @@ void PairPACEKokkos::operator() (TagPairPACEComputeNeigh,const typen } offset++; }); + + if(is_zbl) { + //adapted from https://www.osti.gov/servlets/purl/1429450 + using minloc_value_type=Kokkos::MinLoc::value_type; + minloc_value_type djjmin; + djjmin.val=1e20; + djjmin.loc=-1; + Kokkos::MinLoc reducer_scalar(djjmin); + // loop over ncount (actual neighbours withing cutoff) rather than jnum (total number of neigh in cutoff+skin) + Kokkos::parallel_reduce(Kokkos::TeamThreadRange(team, ncount), + [&](const int offset, minloc_value_type &min_d_dist) { + int j = d_nearest(ii,offset); + j &= NEIGHMASK; + const int jtype = type(j); + auto r = d_rnorms(ii,offset); + const int mu_j = d_map(type(j)); + const F_FLOAT d = r - (d_cut_in(mu_i, mu_j) - d_dcut_in(mu_i, mu_j)); + if (d < min_d_dist.val) { + min_d_dist.val = d; + min_d_dist.loc = offset; + } + + }, reducer_scalar); + d_d_min(ii) = djjmin.val; + d_jj_min(ii) = djjmin.loc;// d_jj_min should be NOT in 0..jnum range, but in 0..d_ncount(<=jnum) + } } /* ---------------------------------------------------------------------- */ @@ -990,22 +1040,38 @@ void PairPACEKokkos::operator() (TagPairPACEComputeFS, const int& ii const int ndensity = d_ndensity(mu_i); double evdwl, fcut, dfcut; + double evdwl_cut; evdwl = fcut = dfcut = 0.0; - inner_cutoff(rho_core(ii), rho_cut, drho_cut, fcut, dfcut); FS_values_and_derivatives(ii, evdwl, mu_i); - - dF_drho_core(ii) = evdwl * dfcut + 1; + if(is_zbl) { + if (d_jj_min(ii) != -1) { + const int mu_jmin = d_mu(ii,d_jj_min(ii)); + F_FLOAT dcutin = d_dcut_in(mu_i, mu_jmin); + F_FLOAT transition_coordinate = dcutin - d_d_min(ii); // == cutin - r_min + cutoff_func_poly(transition_coordinate, dcutin, dcutin, fcut, dfcut); + dfcut = -dfcut; // invert, because rho_core = cutin - r_min + } else { + // no neighbours + fcut = 1; + dfcut = 0; + } + evdwl_cut = evdwl * fcut + rho_core(ii) * (1 - fcut); // evdwl * fcut + rho_core_uncut - rho_core_uncut* fcut + dF_drho_core(ii) = 1 - fcut; + dF_dfcut(ii) = evdwl * dfcut - rho_core(ii) * dfcut; + } else { + inner_cutoff(rho_core(ii), rho_cut, drho_cut, fcut, dfcut); + dF_drho_core(ii) = evdwl * dfcut + 1; + evdwl_cut = evdwl * fcut + rho_core(ii); + } for (int p = 0; p < ndensity; ++p) - dF_drho(ii, p) *= fcut; - + dF_drho(ii, p) *= fcut; // tally energy contribution if (eflag) { - double evdwl_cut = evdwl * fcut + rho_core(ii); - // E0 shift - evdwl_cut += d_E0vals(mu_i); - e_atom(ii) = evdwl_cut; + // E0 shift + evdwl_cut += d_E0vals(mu_i); + e_atom(ii) = evdwl_cut; } } @@ -1146,6 +1212,15 @@ void PairPACEKokkos::operator() (TagPairPACEComputeDerivative, const f_ij(ii, jj, 0) = scale * f_ji[0] + fpair * r_hat[0]; f_ij(ii, jj, 1) = scale * f_ji[1] + fpair * r_hat[1]; f_ij(ii, jj, 2) = scale * f_ji[2] + fpair * r_hat[2]; + + if(is_zbl) { + if(jj==d_jj_min(ii)) { + // DCRU = 1.0 + f_ij(ii, jj, 0) += dF_dfcut(ii) * r_hat[0]; + f_ij(ii, jj, 1) += dF_dfcut(ii) * r_hat[1]; + f_ij(ii, jj, 2) += dF_dfcut(ii) * r_hat[2]; + } + } } /* ---------------------------------------------------------------------- */ @@ -1683,6 +1758,7 @@ double PairPACEKokkos::memory_usage() bytes += MemKK::memory_usage(weights_rank1); bytes += MemKK::memory_usage(rho_core); bytes += MemKK::memory_usage(dF_drho_core); + bytes += MemKK::memory_usage(dF_dfcut); bytes += MemKK::memory_usage(dB_flatten); bytes += MemKK::memory_usage(fr); bytes += MemKK::memory_usage(dfr); @@ -1700,6 +1776,8 @@ double PairPACEKokkos::memory_usage() bytes += MemKK::memory_usage(d_mu); bytes += MemKK::memory_usage(d_rhats); bytes += MemKK::memory_usage(d_rnorms); + bytes += MemKK::memory_usage(d_d_min); + bytes += MemKK::memory_usage(d_jj_min); bytes += MemKK::memory_usage(d_nearest); bytes += MemKK::memory_usage(f_ij); bytes += MemKK::memory_usage(d_rho_core_cutoff); diff --git a/src/KOKKOS/pair_pace_kokkos.h b/src/KOKKOS/pair_pace_kokkos.h index 39cfd100f8..4a69e8e258 100644 --- a/src/KOKKOS/pair_pace_kokkos.h +++ b/src/KOKKOS/pair_pace_kokkos.h @@ -121,6 +121,7 @@ class PairPACEKokkos : public PairPACE { tdual_fparams k_cutsq, k_scale; typedef Kokkos::View t_fparams; t_fparams d_cutsq, d_scale; + t_fparams d_cut_in, d_dcut_in; // inner cutoff typename AT::t_int_1d d_map; @@ -228,6 +229,7 @@ class PairPACEKokkos : public PairPACE { t_ace_2d cr; t_ace_2d dcr; t_ace_1d dF_drho_core; + t_ace_1d dF_dfcut; // radial functions t_ace_4d fr; @@ -265,6 +267,11 @@ class PairPACEKokkos : public PairPACE { t_ace_3d3 d_rhats; t_ace_2i d_nearest; + // for ZBL core-rep implementation + t_ace_1d d_d_min; // [i] -> min-d for atom ii, d=d = r - (cut_in(mu_i, mu_j) - dcut_in(mu_i, mu_j)) + t_ace_1i d_jj_min; // [i] -> jj-index of nearest neigh (by r-(cut_in-dcut_in) criterion) + bool is_zbl; + // per-type t_ace_1i d_ndensity; t_ace_1i d_npoti; diff --git a/src/ML-PACE/pair_pace.cpp b/src/ML-PACE/pair_pace.cpp index 57f12597d1..7ef6789598 100644 --- a/src/ML-PACE/pair_pace.cpp +++ b/src/ML-PACE/pair_pace.cpp @@ -45,7 +45,8 @@ Copyright 2021 Yury Lysogorskiy^1, Cas van der Oord^2, Anton Bochkarev^1, #include "ace-evaluator/ace_evaluator.h" #include "ace-evaluator/ace_recursive.h" #include "ace-evaluator/ace_version.h" - +#include "ace/ace_b_basis.h" +#include namespace LAMMPS_NS { struct ACEImpl { ACEImpl() : basis_set(nullptr), ace(nullptr) {} @@ -287,7 +288,14 @@ void PairPACE::coeff(int narg, char **arg) //load potential file delete aceimpl->basis_set; if (comm->me == 0) utils::logmesg(lmp, "Loading {}\n", potential_file_name); - aceimpl->basis_set = new ACECTildeBasisSet(potential_file_name); + // if potential is in ACEBBasisSet (YAML) format, then convert to ACECTildeBasisSet automatically + if (utils::strmatch(potential_file_name,".*\\.yaml$")) { + ACEBBasisSet bBasisSet = ACEBBasisSet(potential_file_name); + ACECTildeBasisSet cTildeBasisSet = bBasisSet.to_ACECTildeBasisSet(); + aceimpl->basis_set = new ACECTildeBasisSet(cTildeBasisSet); + } else { + aceimpl->basis_set = new ACECTildeBasisSet(potential_file_name); + } if (comm->me == 0) { utils::logmesg(lmp, "Total number of basis functions\n"); From 2e421b2eace20f6c42c185325f81073d69b86d3f Mon Sep 17 00:00:00 2001 From: Yury Lysogorskiy Date: Wed, 22 Nov 2023 10:36:39 +0100 Subject: [PATCH 20/29] fix ML_PACE.cmake --- cmake/Modules/Packages/ML-PACE.cmake | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cmake/Modules/Packages/ML-PACE.cmake b/cmake/Modules/Packages/ML-PACE.cmake index 8a92dee395..50072f2d28 100644 --- a/cmake/Modules/Packages/ML-PACE.cmake +++ b/cmake/Modules/Packages/ML-PACE.cmake @@ -16,8 +16,8 @@ else() endif() if(NOT "${DL_MD5}" STREQUAL "${PACELIB_MD5}") message(STATUS "Downloading ${PACELIB_URL}") - #file(DOWNLOAD ${PACELIB_URL} ${CMAKE_BINARY_DIR}/libpace.tar.gz STATUS DL_STATUS SHOW_PROGRESS) - #file(MD5 ${CMAKE_BINARY_DIR}/libpace.tar.gz DL_MD5) + file(DOWNLOAD ${PACELIB_URL} ${CMAKE_BINARY_DIR}/libpace.tar.gz STATUS DL_STATUS SHOW_PROGRESS) + file(MD5 ${CMAKE_BINARY_DIR}/libpace.tar.gz DL_MD5) if((NOT DL_STATUS EQUAL 0) OR (NOT "${DL_MD5}" STREQUAL "${PACELIB_MD5}")) message(WARNING "Download from primary URL ${PACELIB_URL} failed\nTrying fallback URL ${PACELIB_FALLBACK}") file(DOWNLOAD ${PACELIB_FALLBACK} ${CMAKE_BINARY_DIR}/libpace.tar.gz EXPECTED_HASH MD5=${PACELIB_MD5} SHOW_PROGRESS) From bf498022cc339d05f8d1097f51dc84bbc014e7a6 Mon Sep 17 00:00:00 2001 From: Yury Lysogorskiy Date: Wed, 22 Nov 2023 14:44:49 +0100 Subject: [PATCH 21/29] use inum instead of list->inum --- src/ML-PACE/pair_pace.cpp | 4 ++-- src/ML-PACE/pair_pace_extrapolation.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ML-PACE/pair_pace.cpp b/src/ML-PACE/pair_pace.cpp index 7ef6789598..5ceb426b7d 100644 --- a/src/ML-PACE/pair_pace.cpp +++ b/src/ML-PACE/pair_pace.cpp @@ -147,7 +147,7 @@ void PairPACE::compute(int eflag, int vflag) //determine the maximum number of neighbours int max_jnum = 0; int nei = 0; - for (ii = 0; ii < list->inum; ii++) { + for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jnum = numneigh[i]; nei = nei + jnum; @@ -157,7 +157,7 @@ void PairPACE::compute(int eflag, int vflag) aceimpl->ace->resize_neighbours_cache(max_jnum); //loop over atoms - for (ii = 0; ii < list->inum; ii++) { + for (ii = 0; ii < inum; ii++) { i = list->ilist[ii]; const int itype = type[i]; diff --git a/src/ML-PACE/pair_pace_extrapolation.cpp b/src/ML-PACE/pair_pace_extrapolation.cpp index dc0fb1848b..ebd4018fc7 100644 --- a/src/ML-PACE/pair_pace_extrapolation.cpp +++ b/src/ML-PACE/pair_pace_extrapolation.cpp @@ -170,7 +170,7 @@ void PairPACEExtrapolation::compute(int eflag, int vflag) //determine the maximum number of neighbours int max_jnum = 0; int nei = 0; - for (ii = 0; ii < list->inum; ii++) { + for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jnum = numneigh[i]; nei = nei + jnum; @@ -183,7 +183,7 @@ void PairPACEExtrapolation::compute(int eflag, int vflag) aceimpl->rec_ace->resize_neighbours_cache(max_jnum); //loop over atoms - for (ii = 0; ii < list->inum; ii++) { + for (ii = 0; ii < inum; ii++) { i = list->ilist[ii]; const int itype = type[i]; From c0631c9bd2a9f52ab42ac16cd83992b9d9d51c0f Mon Sep 17 00:00:00 2001 From: Yury Lysogorskiy Date: Sat, 25 Nov 2023 23:04:18 +0100 Subject: [PATCH 22/29] add per-atom array corerep_factor to all four pair styles: pace, pace/extrapolation, pace/kk, pace/extrapolation/kk upd PACELIB_URL to v.2023.11.25 (+md5sum) in ML-PACE.cmake and Install.py --- cmake/Modules/Packages/ML-PACE.cmake | 4 +- lib/pace/Install.py | 4 +- src/KOKKOS/pair_pace_extrapolation_kokkos.cpp | 39 ++++++++++--- src/KOKKOS/pair_pace_extrapolation_kokkos.h | 7 ++- src/KOKKOS/pair_pace_kokkos.cpp | 57 +++++++++++++++++++ src/KOKKOS/pair_pace_kokkos.h | 10 +++- src/ML-PACE/pair_pace.cpp | 38 +++++++++++++ src/ML-PACE/pair_pace.h | 4 ++ src/ML-PACE/pair_pace_extrapolation.cpp | 25 +++++++- src/ML-PACE/pair_pace_extrapolation.h | 8 ++- 10 files changed, 177 insertions(+), 19 deletions(-) diff --git a/cmake/Modules/Packages/ML-PACE.cmake b/cmake/Modules/Packages/ML-PACE.cmake index 50072f2d28..1d76012803 100644 --- a/cmake/Modules/Packages/ML-PACE.cmake +++ b/cmake/Modules/Packages/ML-PACE.cmake @@ -1,6 +1,6 @@ -set(PACELIB_URL "https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/v.2023.11.22.tar.gz" CACHE STRING "URL for PACE evaluator library sources") +set(PACELIB_URL "https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/v.2023.11.25.tar.gz" CACHE STRING "URL for PACE evaluator library sources") -set(PACELIB_MD5 "c8e811f96d761ef8863f5b88a3fd36f4" CACHE STRING "MD5 checksum of PACE evaluator library tarball") +set(PACELIB_MD5 "fda61c88a6c6a335d5fc10a86a5aa710" CACHE STRING "MD5 checksum of PACE evaluator library tarball") mark_as_advanced(PACELIB_URL) mark_as_advanced(PACELIB_MD5) GetFallbackURL(PACELIB_URL PACELIB_FALLBACK) diff --git a/lib/pace/Install.py b/lib/pace/Install.py index 8b89b3dec1..0d5294e00b 100644 --- a/lib/pace/Install.py +++ b/lib/pace/Install.py @@ -18,11 +18,11 @@ from install_helpers import fullpath, geturl, checkmd5sum, getfallback # settings thisdir = fullpath('.') -version ='v.2023.11.22' +version ='v.2023.11.25' # known checksums for different PACE versions. used to validate the download. checksums = { \ - 'v.2023.11.22': 'c8e811f96d761ef8863f5b88a3fd36f4' + 'v.2023.11.25': 'fda61c88a6c6a335d5fc10a86a5aa710' } parser = ArgumentParser(prog='Install.py', description="LAMMPS library build wrapper script") diff --git a/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp b/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp index b361cfe486..d47c8a432a 100644 --- a/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp +++ b/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp @@ -126,6 +126,7 @@ void PairPACEExtrapolationKokkos::grow(int natom, int maxneigh) MemKK::realloc_kokkos(dF_dfcut, "pace:dF_dfcut", natom); MemKK::realloc_kokkos(d_d_min, "pace:r_min_pair", natom); MemKK::realloc_kokkos(d_jj_min, "pace:j_min_pair", natom); + MemKK::realloc_kokkos(d_corerep, "pace:corerep", natom); // per-atom corerep MemKK::realloc_kokkos(dB_flatten, "pace:dB_flatten", natom, idx_ms_combs_max, basis_set->rankmax); @@ -594,13 +595,20 @@ void PairPACEExtrapolationKokkos::compute(int eflag_in, int vflag_in d_vatom = k_vatom.view(); } - if (gamma_flag && atom->nlocal > nmax) { + if (flag_compute_extrapolation_grade && atom->nlocal > nmax) { memory->destroy(extrapolation_grade_gamma); nmax = atom->nlocal; memory->create(extrapolation_grade_gamma, nmax, "pace/atom:gamma"); //zeroify array memset(extrapolation_grade_gamma, 0, nmax * sizeof(*extrapolation_grade_gamma)); } + if (flag_corerep_factor && atom->nlocal > nmax_corerep) { + memory->destroy(corerep_factor); + nmax_corerep = atom->nlocal; + memory->create(corerep_factor, nmax_corerep, "pace/atom:corerep"); + //zeroify array + memset(corerep_factor, 0, nmax_corerep * sizeof(*corerep_factor)); + } copymode = 1; if (!force->newton_pair) @@ -658,6 +666,7 @@ void PairPACEExtrapolationKokkos::compute(int eflag_in, int vflag_in Kokkos::deep_copy(projections, 0.0); Kokkos::deep_copy(d_gamma, 0.0); + Kokkos::deep_copy(d_corerep, 0.0); EV_FLOAT ev_tmp; @@ -721,7 +730,7 @@ void PairPACEExtrapolationKokkos::compute(int eflag_in, int vflag_in } //ComputeGamma - if (gamma_flag) { + if (flag_compute_extrapolation_grade) { typename Kokkos::RangePolicy policy_gamma(0,chunk_size); Kokkos::parallel_for("ComputeGamma",policy_gamma,*this); } @@ -763,12 +772,17 @@ void PairPACEExtrapolationKokkos::compute(int eflag_in, int vflag_in } ev += ev_tmp; - //if gamma_flag - copy current d_gamma to extrapolation_grade_gamma - if (gamma_flag){ + //if flag_compute_extrapolation_grade - copy current d_gamma to extrapolation_grade_gamma + if (flag_compute_extrapolation_grade){ h_gamma = Kokkos::create_mirror_view(d_gamma); Kokkos::deep_copy(h_gamma, d_gamma); memcpy(extrapolation_grade_gamma+chunk_offset, (void *) h_gamma.data(), sizeof(double)*chunk_size); } + if (flag_corerep_factor) { + h_corerep = Kokkos::create_mirror_view(d_corerep); + Kokkos::deep_copy(h_corerep,d_corerep); + memcpy(corerep_factor+chunk_offset, (void *) h_corerep.data(), sizeof(double)*chunk_size); + } chunk_offset += chunk_size; } // end while @@ -1050,7 +1064,7 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeRho, //gamma_i - if (gamma_flag) + if (flag_compute_extrapolation_grade) Kokkos::atomic_add(&projections(ii, func_ind), d_gen_cgs(mu_i, idx_ms_comb) * A_cur); } else { // rank > 1 @@ -1089,7 +1103,7 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeRho, Kokkos::atomic_add(&rhos(ii, p), B.real_part_product(d_coeffs(mu_i, func_ind, p) * d_gen_cgs(mu_i, idx_ms_comb))); } //gamma_i - if (gamma_flag) + if (flag_compute_extrapolation_grade) Kokkos::atomic_add(&projections(ii, func_ind), B.real_part_product(d_gen_cgs(mu_i, idx_ms_comb))); } } @@ -1144,6 +1158,9 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeFS, evdwl_cut += d_E0vals(mu_i); e_atom(ii) = evdwl_cut; } + + if (flag_corerep_factor) + d_corerep(ii) = 1-fcut; } /* ---------------------------------------------------------------------- */ @@ -1858,6 +1875,7 @@ double PairPACEExtrapolationKokkos::memory_usage() bytes += MemKK::memory_usage(rho_core); bytes += MemKK::memory_usage(dF_drho_core); bytes += MemKK::memory_usage(dF_dfcut); + bytes += MemKK::memory_usage(d_corerep); bytes += MemKK::memory_usage(dB_flatten); bytes += MemKK::memory_usage(fr); bytes += MemKK::memory_usage(dfr); @@ -1925,9 +1943,10 @@ double PairPACEExtrapolationKokkos::memory_usage() template void *PairPACEExtrapolationKokkos::extract(const char *str, int &dim) { - //check if str=="gamma_flag" then compute extrapolation grades on this iteration dim = 0; - if (strcmp(str, "gamma_flag") == 0) return (void *) &gamma_flag; + //check if str=="flag_compute_extrapolation_grade" then compute extrapolation grades on this iteration + if (strcmp(str, "gamma_flag") == 0) return (void *) &flag_compute_extrapolation_grade; + if (strcmp(str, "corerep_flag") == 0) return (void *) &flag_corerep_factor; dim = 2; if (strcmp(str, "scale") == 0) return (void *) scale; @@ -1950,6 +1969,10 @@ void *PairPACEExtrapolationKokkos::extract_peratom(const char *str, ncol = 0; return (void *) extrapolation_grade_gamma; } + if (strcmp(str, "corerep") == 0) { + ncol = 0; + return (void *) corerep_factor; + } return nullptr; } diff --git a/src/KOKKOS/pair_pace_extrapolation_kokkos.h b/src/KOKKOS/pair_pace_extrapolation_kokkos.h index 092171cfb8..d8c0038f3b 100644 --- a/src/KOKKOS/pair_pace_extrapolation_kokkos.h +++ b/src/KOKKOS/pair_pace_extrapolation_kokkos.h @@ -106,7 +106,7 @@ class PairPACEExtrapolationKokkos : public PairPACEExtrapolation { protected: int inum, maxneigh, chunk_size, chunk_offset, idx_ms_combs_max, total_num_functions_max; int host_flag; - int gamma_flag; + //int gamma_flag; int eflag, vflag; @@ -235,13 +235,16 @@ class PairPACEExtrapolationKokkos : public PairPACEExtrapolation { t_ace_2d rhos; t_ace_2d dF_drho; + t_ace_3c dB_flatten; + // hard-core repulsion t_ace_1d rho_core; - t_ace_3c dB_flatten; t_ace_2d cr; t_ace_2d dcr; t_ace_1d dF_drho_core; t_ace_1d dF_dfcut; + t_ace_1d d_corerep; + th_ace_1d h_corerep; // radial functions t_ace_4d fr; diff --git a/src/KOKKOS/pair_pace_kokkos.cpp b/src/KOKKOS/pair_pace_kokkos.cpp index 62b7bf5d01..84431a5908 100644 --- a/src/KOKKOS/pair_pace_kokkos.cpp +++ b/src/KOKKOS/pair_pace_kokkos.cpp @@ -124,6 +124,8 @@ void PairPACEKokkos::grow(int natom, int maxneigh) MemKK::realloc_kokkos(dF_dfcut, "pace:dF_dfcut", natom); MemKK::realloc_kokkos(d_d_min, "pace:r_min_pair", natom); MemKK::realloc_kokkos(d_jj_min, "pace:j_min_pair", natom); + MemKK::realloc_kokkos(d_corerep, "pace:corerep", natom); // per-atom corerep + MemKK::realloc_kokkos(dB_flatten, "pace:dB_flatten", natom, idx_rho_max, basis_set->rankmax); } @@ -555,6 +557,13 @@ void PairPACEKokkos::compute(int eflag_in, int vflag_in) memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom"); d_vatom = k_vatom.view(); } + if (flag_corerep_factor && atom->nlocal > nmax_corerep) { + memory->destroy(corerep_factor); + nmax_corerep = atom->nlocal; + memory->create(corerep_factor, nmax_corerep, "pace/atom:corerep"); + //zeroify array + memset(corerep_factor, 0, nmax_corerep * sizeof(*corerep_factor)); + } copymode = 1; if (!force->newton_pair) @@ -610,6 +619,7 @@ void PairPACEKokkos::compute(int eflag_in, int vflag_in) Kokkos::deep_copy(rho_core, 0.0); Kokkos::deep_copy(d_d_min, PairPACE::aceimpl->basis_set->cutoffmax); Kokkos::deep_copy(d_jj_min, -1); + Kokkos::deep_copy(d_corerep, 0.0); EV_FLOAT ev_tmp; @@ -709,6 +719,13 @@ void PairPACEKokkos::compute(int eflag_in, int vflag_in) } } ev += ev_tmp; + + if (flag_corerep_factor) { + h_corerep = Kokkos::create_mirror_view(d_corerep); + Kokkos::deep_copy(h_corerep,d_corerep); + memcpy(corerep_factor+chunk_offset, (void *) h_corerep.data(), sizeof(double)*chunk_size); + } + chunk_offset += chunk_size; } // end while @@ -1073,6 +1090,9 @@ void PairPACEKokkos::operator() (TagPairPACEComputeFS, const int& ii evdwl_cut += d_E0vals(mu_i); e_atom(ii) = evdwl_cut; } + + if (flag_corerep_factor) + d_corerep(ii) = 1-fcut; } /* ---------------------------------------------------------------------- */ @@ -1759,6 +1779,7 @@ double PairPACEKokkos::memory_usage() bytes += MemKK::memory_usage(rho_core); bytes += MemKK::memory_usage(dF_drho_core); bytes += MemKK::memory_usage(dF_dfcut); + bytes += MemKK::memory_usage(d_corerep); bytes += MemKK::memory_usage(dB_flatten); bytes += MemKK::memory_usage(fr); bytes += MemKK::memory_usage(dfr); @@ -1814,6 +1835,42 @@ double PairPACEKokkos::memory_usage() return bytes; } + +/* ---------------------------------------------------------------------- + extract method for extracting value of scale variable + ---------------------------------------------------------------------- */ + +template +void *PairPACEKokkos::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str, "corerep_flag") == 0) return (void *) &flag_corerep_factor; + + dim = 2; + if (strcmp(str, "scale") == 0) return (void *) scale; + return nullptr; +} + +/* ---------------------------------------------------------------------- + peratom requests from FixPair + return ptr to requested data + also return ncol = # of quantites per atom + 0 = per-atom vector + 1 or more = # of columns in per-atom array + return NULL if str is not recognized +---------------------------------------------------------------------- */ + +template +void *PairPACEKokkos::extract_peratom(const char *str, int &ncol) +{ + if (strcmp(str, "corerep") == 0) { + ncol = 0; + return (void *) corerep_factor; + } + + return nullptr; +} + /* ---------------------------------------------------------------------- */ namespace LAMMPS_NS { diff --git a/src/KOKKOS/pair_pace_kokkos.h b/src/KOKKOS/pair_pace_kokkos.h index 4a69e8e258..2c5aaa846e 100644 --- a/src/KOKKOS/pair_pace_kokkos.h +++ b/src/KOKKOS/pair_pace_kokkos.h @@ -95,6 +95,9 @@ class PairPACEKokkos : public PairPACE { KOKKOS_INLINE_FUNCTION void operator() (TagPairPACEComputeForce,const int& ii, EV_FLOAT&) const; + void *extract(const char *str, int &dim) override; + void *extract_peratom(const char *str, int &ncol) override; + protected: int inum, maxneigh, chunk_size, chunk_offset, idx_rho_max; int host_flag; @@ -210,6 +213,8 @@ class PairPACEKokkos : public PairPACE { typedef Kokkos::View t_ace_4c; typedef Kokkos::View t_ace_4c3; + typedef typename Kokkos::View::HostMirror th_ace_1d; + t_ace_3d A_rank1; t_ace_4c A; @@ -223,13 +228,16 @@ class PairPACEKokkos : public PairPACE { t_ace_2d rhos; t_ace_2d dF_drho; + t_ace_3c dB_flatten; + // hard-core repulsion t_ace_1d rho_core; - t_ace_3c dB_flatten; t_ace_2d cr; t_ace_2d dcr; t_ace_1d dF_drho_core; t_ace_1d dF_dfcut; + t_ace_1d d_corerep; + th_ace_1d h_corerep; // radial functions t_ace_4d fr; diff --git a/src/ML-PACE/pair_pace.cpp b/src/ML-PACE/pair_pace.cpp index 5ceb426b7d..933d4abecb 100644 --- a/src/ML-PACE/pair_pace.cpp +++ b/src/ML-PACE/pair_pace.cpp @@ -88,6 +88,10 @@ PairPACE::PairPACE(LAMMPS *lmp) : Pair(lmp) one_coeff = 1; manybody_flag = 1; + nmax_corerep = 0; + flag_corerep_factor = 0; + corerep_factor = nullptr; + aceimpl = new ACEImpl; recursive = false; @@ -110,6 +114,7 @@ PairPACE::~PairPACE() memory->destroy(setflag); memory->destroy(cutsq); memory->destroy(scale); + memory->destroy(corerep_factor); } } @@ -144,6 +149,14 @@ void PairPACE::compute(int eflag, int vflag) // the pointer to the list of neighbors of "i" firstneigh = list->firstneigh; + if (flag_corerep_factor && atom->nlocal > nmax_corerep) { + memory->destroy(corerep_factor); + nmax_corerep = atom->nlocal; + memory->create(corerep_factor, nmax_corerep, "pace/atom:corerep_factor"); + //zeroify array + memset(corerep_factor, 0, nmax_corerep * sizeof(*corerep_factor)); + } + //determine the maximum number of neighbours int max_jnum = 0; int nei = 0; @@ -182,6 +195,9 @@ void PairPACE::compute(int eflag, int vflag) error->one(FLERR, e.what()); } + if (flag_corerep_factor) + corerep_factor[i] = 1 - aceimpl->ace->ace_fcut; + // 'compute_atom' will update the `aceimpl->ace->e_atom` and `aceimpl->ace->neighbours_forces(jj, alpha)` arrays for (jj = 0; jj < jnum; jj++) { @@ -382,7 +398,29 @@ double PairPACE::init_one(int i, int j) ---------------------------------------------------------------------- */ void *PairPACE::extract(const char *str, int &dim) { + dim = 0; + //check if str=="corerep_flag" then compute extrapolation grades on this iteration + if (strcmp(str, "corerep_flag") == 0) return (void *) &flag_corerep_factor; + dim = 2; if (strcmp(str, "scale") == 0) return (void *) scale; return nullptr; } + +/* ---------------------------------------------------------------------- + peratom requests from FixPair + return ptr to requested data + also return ncol = # of quantites per atom + 0 = per-atom vector + 1 or more = # of columns in per-atom array + return NULL if str is not recognized +---------------------------------------------------------------------- */ +void *PairPACE::extract_peratom(const char *str, int &ncol) +{ + if (strcmp(str, "corerep") == 0) { + ncol = 0; + return (void *) corerep_factor; + } + + return nullptr; +} \ No newline at end of file diff --git a/src/ML-PACE/pair_pace.h b/src/ML-PACE/pair_pace.h index 94649ecaab..9b5d2c5480 100644 --- a/src/ML-PACE/pair_pace.h +++ b/src/ML-PACE/pair_pace.h @@ -48,11 +48,15 @@ class PairPACE : public Pair { double init_one(int, int) override; void *extract(const char *, int &) override; + void *extract_peratom(const char *, int &) override; protected: struct ACEImpl *aceimpl; + int nmax_corerep = 0; virtual void allocate(); + double *corerep_factor; //per-atom core-rep factor (= 1 - fcut) + int flag_corerep_factor; double **scale; bool recursive; // "recursive" option for ACERecursiveEvaluator diff --git a/src/ML-PACE/pair_pace_extrapolation.cpp b/src/ML-PACE/pair_pace_extrapolation.cpp index ebd4018fc7..d9b8d3588a 100644 --- a/src/ML-PACE/pair_pace_extrapolation.cpp +++ b/src/ML-PACE/pair_pace_extrapolation.cpp @@ -93,11 +93,14 @@ PairPACEExtrapolation::PairPACEExtrapolation(LAMMPS *lmp) : Pair(lmp) manybody_flag = 1; nmax = 0; + nmax_corerep = 0; aceimpl = new ACEALImpl; scale = nullptr; flag_compute_extrapolation_grade = 0; extrapolation_grade_gamma = nullptr; + flag_corerep_factor = 0; + corerep_factor = nullptr; chunksize = 4096; } @@ -118,6 +121,7 @@ PairPACEExtrapolation::~PairPACEExtrapolation() memory->destroy(scale); memory->destroy(map); memory->destroy(extrapolation_grade_gamma); + memory->destroy(corerep_factor); } } @@ -166,6 +170,13 @@ void PairPACEExtrapolation::compute(int eflag, int vflag) //zeroify array memset(extrapolation_grade_gamma, 0, nmax * sizeof(*extrapolation_grade_gamma)); } + if (flag_corerep_factor && atom->nlocal > nmax_corerep) { + memory->destroy(corerep_factor); + nmax_corerep = atom->nlocal; + memory->create(corerep_factor, nmax_corerep, "pace/atom:corerep_factor"); + //zeroify array + memset(corerep_factor, 0, nmax_corerep * sizeof(*corerep_factor)); + } //determine the maximum number of neighbours int max_jnum = 0; @@ -216,6 +227,11 @@ void PairPACEExtrapolation::compute(int eflag, int vflag) if (flag_compute_extrapolation_grade) extrapolation_grade_gamma[i] = aceimpl->ace->max_gamma_grade; + if (flag_corerep_factor) { + corerep_factor[i] = 1 - (flag_compute_extrapolation_grade ? aceimpl->ace->ace_fcut + : aceimpl->rec_ace->ace_fcut); + } + Array2D &neighbours_forces = (flag_compute_extrapolation_grade ? aceimpl->ace->neighbours_forces : aceimpl->rec_ace->neighbours_forces); @@ -437,9 +453,11 @@ double PairPACEExtrapolation::init_one(int i, int j) ---------------------------------------------------------------------- */ void *PairPACEExtrapolation::extract(const char *str, int &dim) { - //check if str=="gamma_flag" then compute extrapolation grades on this iteration dim = 0; + //check if str=="gamma_flag" then compute extrapolation grades on this iteration if (strcmp(str, "gamma_flag") == 0) return (void *) &flag_compute_extrapolation_grade; + //check if str=="corerep_flag" then compute extrapolation grades on this iteration + if (strcmp(str, "corerep_flag") == 0) return (void *) &flag_corerep_factor; dim = 2; if (strcmp(str, "scale") == 0) return (void *) scale; @@ -461,5 +479,10 @@ void *PairPACEExtrapolation::extract_peratom(const char *str, int &ncol) return (void *) extrapolation_grade_gamma; } + if (strcmp(str, "corerep") == 0) { + ncol = 0; + return (void *) corerep_factor; + } + return nullptr; } diff --git a/src/ML-PACE/pair_pace_extrapolation.h b/src/ML-PACE/pair_pace_extrapolation.h index 6f7eeb279e..2dcec04d4b 100644 --- a/src/ML-PACE/pair_pace_extrapolation.h +++ b/src/ML-PACE/pair_pace_extrapolation.h @@ -47,13 +47,15 @@ class PairPACEExtrapolation : public Pair { protected: struct ACEALImpl *aceimpl; - int nmax; + int nmax = 0, nmax_corerep = 0; virtual void allocate(); std::vector element_names; // list of elements (used by dump pace/extrapolation) - double *extrapolation_grade_gamma; //per-atom gamma value + double *extrapolation_grade_gamma = nullptr; //per-atom gamma value + double *corerep_factor = nullptr; //per-atom core-rep factor (= 1 - fcut) - int flag_compute_extrapolation_grade; + int flag_compute_extrapolation_grade = 0; + int flag_corerep_factor = 0; double **scale; From d4fe21f34df7e02f9df359547b381924926b36ba Mon Sep 17 00:00:00 2001 From: Yury Lysogorskiy Date: Sat, 25 Nov 2023 23:37:01 +0100 Subject: [PATCH 23/29] update pair_pace.rst documentation page --- doc/src/pair_pace.rst | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/doc/src/pair_pace.rst b/doc/src/pair_pace.rst index d815f663fe..001214370c 100644 --- a/doc/src/pair_pace.rst +++ b/doc/src/pair_pace.rst @@ -40,6 +40,9 @@ Examples pair_style pace product chunksize 2048 pair_coeff * * Cu-PBE-core-rep.ace Cu + pair_style pace + pair_coeff * * Cu.yaml Cu + pair_style pace/extrapolation pair_coeff * * Cu.yaml Cu.asi Cu @@ -64,7 +67,7 @@ specifies an ACE coefficient file followed by N additional arguments specifying the mapping of ACE elements to LAMMPS atom types, where N is the number of LAMMPS atom types: -* ACE coefficient file +* ACE coefficient file (.yaml or .yace/.ace format) * N element names = mapping of ACE elements to atom types Only a single pair_coeff command is used with the *pace* style which @@ -136,6 +139,22 @@ product B-basis evaluator is always used and only *linear* ASI is supported. See the :doc:`pair_coeff ` page for alternate ways to specify the path for the ACE coefficient file. +Core repulsion +""""""""""""""""""" +The ACE potential can be configured to initiate core-repulsion from an inner cutoff, +seamlessly transitioning from ACE to ZBL. The core repulsion factor can be accessed +as a per-atom quantity, as demonstrated in the example below: + +.. code-block:: LAMMPS + + pair_style pace + pair_coeff * * CuNi.yaml Cu Ni + + fix pace_corerep all pair 1 pace corerep 1 + +In this case, per-atom `f_pace_corerep` quantities represent the fraction of ZBL +core-repulsion for each atom. + Mixing, shift, table, tail correction, restart, rRESPA info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" From 99f0a7aa8ec8c6bab20009dcda221019b2bc730a Mon Sep 17 00:00:00 2001 From: Yury Lysogorskiy Date: Mon, 27 Nov 2023 22:38:35 +0100 Subject: [PATCH 24/29] upd version string in ML-PACE package --- cmake/Modules/Packages/ML-PACE.cmake | 4 ++-- lib/pace/Install.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/cmake/Modules/Packages/ML-PACE.cmake b/cmake/Modules/Packages/ML-PACE.cmake index 1d76012803..248b8eea76 100644 --- a/cmake/Modules/Packages/ML-PACE.cmake +++ b/cmake/Modules/Packages/ML-PACE.cmake @@ -1,6 +1,6 @@ -set(PACELIB_URL "https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/v.2023.11.25.tar.gz" CACHE STRING "URL for PACE evaluator library sources") +set(PACELIB_URL "https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/v.2023.11.25.fix.tar.gz" CACHE STRING "URL for PACE evaluator library sources") -set(PACELIB_MD5 "fda61c88a6c6a335d5fc10a86a5aa710" CACHE STRING "MD5 checksum of PACE evaluator library tarball") +set(PACELIB_MD5 "b45de9a633f42ed65422567e3ce56f9f" CACHE STRING "MD5 checksum of PACE evaluator library tarball") mark_as_advanced(PACELIB_URL) mark_as_advanced(PACELIB_MD5) GetFallbackURL(PACELIB_URL PACELIB_FALLBACK) diff --git a/lib/pace/Install.py b/lib/pace/Install.py index 0d5294e00b..fcd9497937 100644 --- a/lib/pace/Install.py +++ b/lib/pace/Install.py @@ -18,11 +18,11 @@ from install_helpers import fullpath, geturl, checkmd5sum, getfallback # settings thisdir = fullpath('.') -version ='v.2023.11.25' +version ='v.2023.11.25.fix' # known checksums for different PACE versions. used to validate the download. checksums = { \ - 'v.2023.11.25': 'fda61c88a6c6a335d5fc10a86a5aa710' + 'v.2023.11.25.fix': 'b45de9a633f42ed65422567e3ce56f9f' } parser = ArgumentParser(prog='Install.py', description="LAMMPS library build wrapper script") From 631dcc36dee037c4216f0854599e7dd7424c769f Mon Sep 17 00:00:00 2001 From: Yury Lysogorskiy Date: Mon, 27 Nov 2023 23:51:51 +0100 Subject: [PATCH 25/29] BUGFIX: TagPairPACEComputeNeigh - explicitly check if atom has neighbours (ncount>0), only them run parallel_reduce, otherwise set d_d_min and d_jj_min to default values --- src/KOKKOS/pair_pace_extrapolation_kokkos.cpp | 47 ++++++++++--------- src/KOKKOS/pair_pace_kokkos.cpp | 43 +++++++++-------- 2 files changed, 50 insertions(+), 40 deletions(-) diff --git a/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp b/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp index d47c8a432a..9d47a4a278 100644 --- a/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp +++ b/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp @@ -901,28 +901,33 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeNeig if(is_zbl) { //adapted from https://www.osti.gov/servlets/purl/1429450 - using minloc_value_type=Kokkos::MinLoc::value_type; - minloc_value_type djjmin; - djjmin.val=1e20; - djjmin.loc=-1; - Kokkos::MinLoc reducer_scalar(djjmin); - // loop over ncount (actual neighbours withing cutoff) rather than jnum (total number of neigh in cutoff+skin) - Kokkos::parallel_reduce(Kokkos::TeamThreadRange(team, ncount), - [&](const int offset, minloc_value_type &min_d_dist) { - int j = d_nearest(ii,offset); - j &= NEIGHMASK; - const int jtype = type(j); - auto r = d_rnorms(ii,offset); - const int mu_j = d_map(type(j)); - const F_FLOAT d = r - (d_cut_in(mu_i, mu_j) - d_dcut_in(mu_i, mu_j)); - if (d < min_d_dist.val) { - min_d_dist.val = d; - min_d_dist.loc = offset; - } + if(ncount>0) { + using minloc_value_type=Kokkos::MinLoc::value_type; + minloc_value_type djjmin; + djjmin.val=1e20; + djjmin.loc=-1; + Kokkos::MinLoc reducer_scalar(djjmin); + // loop over ncount (actual neighbours withing cutoff) rather than jnum (total number of neigh in cutoff+skin) + Kokkos::parallel_reduce(Kokkos::TeamThreadRange(team, ncount), + [&](const int offset, minloc_value_type &min_d_dist) { + int j = d_nearest(ii,offset); + j &= NEIGHMASK; + const int jtype = type(j); + auto r = d_rnorms(ii,offset); + const int mu_j = d_map(type(j)); + const F_FLOAT d = r - (d_cut_in(mu_i, mu_j) - d_dcut_in(mu_i, mu_j)); + if (d < min_d_dist.val) { + min_d_dist.val = d; + min_d_dist.loc = offset; + } - }, reducer_scalar); - d_d_min(ii) = djjmin.val; - d_jj_min(ii) = djjmin.loc;// d_jj_min should be NOT in 0..jnum range, but in 0..d_ncount(<=jnum) + }, reducer_scalar); + d_d_min(ii) = djjmin.val; + d_jj_min(ii) = djjmin.loc;// d_jj_min should be NOT in 0..jnum range, but in 0..d_ncount(<=jnum) + } else { + d_d_min(ii) = 1e20; + d_jj_min(ii) = -1; + } } } diff --git a/src/KOKKOS/pair_pace_kokkos.cpp b/src/KOKKOS/pair_pace_kokkos.cpp index 84431a5908..4e297ec7ed 100644 --- a/src/KOKKOS/pair_pace_kokkos.cpp +++ b/src/KOKKOS/pair_pace_kokkos.cpp @@ -844,28 +844,33 @@ void PairPACEKokkos::operator() (TagPairPACEComputeNeigh,const typen if(is_zbl) { //adapted from https://www.osti.gov/servlets/purl/1429450 - using minloc_value_type=Kokkos::MinLoc::value_type; - minloc_value_type djjmin; - djjmin.val=1e20; - djjmin.loc=-1; - Kokkos::MinLoc reducer_scalar(djjmin); - // loop over ncount (actual neighbours withing cutoff) rather than jnum (total number of neigh in cutoff+skin) - Kokkos::parallel_reduce(Kokkos::TeamThreadRange(team, ncount), - [&](const int offset, minloc_value_type &min_d_dist) { - int j = d_nearest(ii,offset); - j &= NEIGHMASK; - const int jtype = type(j); - auto r = d_rnorms(ii,offset); - const int mu_j = d_map(type(j)); - const F_FLOAT d = r - (d_cut_in(mu_i, mu_j) - d_dcut_in(mu_i, mu_j)); - if (d < min_d_dist.val) { - min_d_dist.val = d; - min_d_dist.loc = offset; - } + if(ncount>0) { + using minloc_value_type=Kokkos::MinLoc::value_type; + minloc_value_type djjmin; + djjmin.val=1e20; + djjmin.loc=-1; + Kokkos::MinLoc reducer_scalar(djjmin); + // loop over ncount (actual neighbours withing cutoff) rather than jnum (total number of neigh in cutoff+skin) + Kokkos::parallel_reduce(Kokkos::TeamThreadRange(team, ncount), + [&](const int offset, minloc_value_type &min_d_dist) { + int j = d_nearest(ii,offset); + j &= NEIGHMASK; + const int jtype = type(j); + auto r = d_rnorms(ii,offset); + const int mu_j = d_map(type(j)); + const F_FLOAT d = r - (d_cut_in(mu_i, mu_j) - d_dcut_in(mu_i, mu_j)); + if (d < min_d_dist.val) { + min_d_dist.val = d; + min_d_dist.loc = offset; + } - }, reducer_scalar); + }, reducer_scalar); d_d_min(ii) = djjmin.val; d_jj_min(ii) = djjmin.loc;// d_jj_min should be NOT in 0..jnum range, but in 0..d_ncount(<=jnum) + } else { + d_d_min(ii) = 1e20; + d_jj_min(ii) = -1; + } } } From 1f7d262bd73ba018912231cad672ad17a4ec63e1 Mon Sep 17 00:00:00 2001 From: jbcouli Date: Mon, 4 Dec 2023 16:01:39 -0700 Subject: [PATCH 26/29] add mol keyword to fix rigid doc --- doc/src/fix_rigid.rst | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/doc/src/fix_rigid.rst b/doc/src/fix_rigid.rst index a50e215681..3174a0929c 100644 --- a/doc/src/fix_rigid.rst +++ b/doc/src/fix_rigid.rst @@ -80,7 +80,7 @@ Syntax groupID1, groupID2, ... = list of N group IDs * zero or more keyword/value pairs may be appended -* keyword = *langevin* or *reinit* or *temp* or *iso* or *aniso* or *x* or *y* or *z* or *couple* or *tparam* or *pchain* or *dilate* or *force* or *torque* or *infile* or *gravity* +* keyword = *langevin* or *reinit* or *temp* or *mol* or *iso* or *aniso* or *x* or *y* or *z* or *couple* or *tparam* or *pchain* or *dilate* or *force* or *torque* or *infile* or *gravity* .. parsed-literal:: @@ -92,6 +92,8 @@ Syntax *temp* values = Tstart Tstop Tdamp Tstart,Tstop = desired temperature at start/stop of run (temperature units) Tdamp = temperature damping parameter (time units) + *mol* value = template-ID + template-ID = ID of molecule template specified in a separate :doc:`molecule ` command *iso* or *aniso* values = Pstart Pstop Pdamp Pstart,Pstop = scalar external pressure at start/end of run (pressure units) Pdamp = pressure damping parameter (time units) From 69fb814b33980f807c24c8438316afc8f04d9825 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 4 Dec 2023 17:31:38 -0700 Subject: [PATCH 27/29] Remove duplicated code, small whitespace tweaks --- src/KOKKOS/pair_pace_extrapolation_kokkos.cpp | 11 ++--- src/KOKKOS/pair_pace_kokkos.cpp | 46 ++----------------- src/KOKKOS/pair_pace_kokkos.h | 3 -- 3 files changed, 9 insertions(+), 51 deletions(-) diff --git a/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp b/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp index 9d47a4a278..0980ad776d 100644 --- a/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp +++ b/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp @@ -899,9 +899,9 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeNeig offset++; }); - if(is_zbl) { + if (is_zbl) { //adapted from https://www.osti.gov/servlets/purl/1429450 - if(ncount>0) { + if (ncount > 0) { using minloc_value_type=Kokkos::MinLoc::value_type; minloc_value_type djjmin; djjmin.val=1e20; @@ -920,7 +920,6 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeNeig min_d_dist.val = d; min_d_dist.loc = offset; } - }, reducer_scalar); d_d_min(ii) = djjmin.val; d_jj_min(ii) = djjmin.loc;// d_jj_min should be NOT in 0..jnum range, but in 0..d_ncount(<=jnum) @@ -1133,7 +1132,7 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeFS, inner_cutoff(rho_core(ii), rho_cut, drho_cut, fcut, dfcut); FS_values_and_derivatives(ii, evdwl, mu_i); - if(is_zbl) { + if (is_zbl) { if (d_jj_min(ii) != -1) { const int mu_jmin = d_mu(ii,d_jj_min(ii)); F_FLOAT dcutin = d_dcut_in(mu_i, mu_jmin); @@ -1334,8 +1333,8 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeDeri f_ij(ii, jj, 1) = scale * f_ji[1] + fpair * r_hat[1]; f_ij(ii, jj, 2) = scale * f_ji[2] + fpair * r_hat[2]; - if(is_zbl) { - if(jj==d_jj_min(ii)) { + if (is_zbl) { + if (jj==d_jj_min(ii)) { // DCRU = 1.0 f_ij(ii, jj, 0) += dF_dfcut(ii) * r_hat[0]; f_ij(ii, jj, 1) += dF_dfcut(ii) * r_hat[1]; diff --git a/src/KOKKOS/pair_pace_kokkos.cpp b/src/KOKKOS/pair_pace_kokkos.cpp index 4e297ec7ed..805d7f68bb 100644 --- a/src/KOKKOS/pair_pace_kokkos.cpp +++ b/src/KOKKOS/pair_pace_kokkos.cpp @@ -697,7 +697,6 @@ void PairPACEKokkos::compute(int eflag_in, int vflag_in) Kokkos::parallel_for("ComputeDerivative",policy_derivative,*this); } - //ComputeForce { if (evflag) { @@ -842,7 +841,7 @@ void PairPACEKokkos::operator() (TagPairPACEComputeNeigh,const typen offset++; }); - if(is_zbl) { + if (is_zbl) { //adapted from https://www.osti.gov/servlets/purl/1429450 if(ncount>0) { using minloc_value_type=Kokkos::MinLoc::value_type; @@ -863,7 +862,6 @@ void PairPACEKokkos::operator() (TagPairPACEComputeNeigh,const typen min_d_dist.val = d; min_d_dist.loc = offset; } - }, reducer_scalar); d_d_min(ii) = djjmin.val; d_jj_min(ii) = djjmin.loc;// d_jj_min should be NOT in 0..jnum range, but in 0..d_ncount(<=jnum) @@ -1066,7 +1064,7 @@ void PairPACEKokkos::operator() (TagPairPACEComputeFS, const int& ii evdwl = fcut = dfcut = 0.0; FS_values_and_derivatives(ii, evdwl, mu_i); - if(is_zbl) { + if (is_zbl) { if (d_jj_min(ii) != -1) { const int mu_jmin = d_mu(ii,d_jj_min(ii)); F_FLOAT dcutin = d_dcut_in(mu_i, mu_jmin); @@ -1238,8 +1236,8 @@ void PairPACEKokkos::operator() (TagPairPACEComputeDerivative, const f_ij(ii, jj, 1) = scale * f_ji[1] + fpair * r_hat[1]; f_ij(ii, jj, 2) = scale * f_ji[2] + fpair * r_hat[2]; - if(is_zbl) { - if(jj==d_jj_min(ii)) { + if (is_zbl) { + if (jj==d_jj_min(ii)) { // DCRU = 1.0 f_ij(ii, jj, 0) += dF_dfcut(ii) * r_hat[0]; f_ij(ii, jj, 1) += dF_dfcut(ii) * r_hat[1]; @@ -1840,42 +1838,6 @@ double PairPACEKokkos::memory_usage() return bytes; } - -/* ---------------------------------------------------------------------- - extract method for extracting value of scale variable - ---------------------------------------------------------------------- */ - -template -void *PairPACEKokkos::extract(const char *str, int &dim) -{ - dim = 0; - if (strcmp(str, "corerep_flag") == 0) return (void *) &flag_corerep_factor; - - dim = 2; - if (strcmp(str, "scale") == 0) return (void *) scale; - return nullptr; -} - -/* ---------------------------------------------------------------------- - peratom requests from FixPair - return ptr to requested data - also return ncol = # of quantites per atom - 0 = per-atom vector - 1 or more = # of columns in per-atom array - return NULL if str is not recognized ----------------------------------------------------------------------- */ - -template -void *PairPACEKokkos::extract_peratom(const char *str, int &ncol) -{ - if (strcmp(str, "corerep") == 0) { - ncol = 0; - return (void *) corerep_factor; - } - - return nullptr; -} - /* ---------------------------------------------------------------------- */ namespace LAMMPS_NS { diff --git a/src/KOKKOS/pair_pace_kokkos.h b/src/KOKKOS/pair_pace_kokkos.h index 2c5aaa846e..36486f8628 100644 --- a/src/KOKKOS/pair_pace_kokkos.h +++ b/src/KOKKOS/pair_pace_kokkos.h @@ -95,9 +95,6 @@ class PairPACEKokkos : public PairPACE { KOKKOS_INLINE_FUNCTION void operator() (TagPairPACEComputeForce,const int& ii, EV_FLOAT&) const; - void *extract(const char *str, int &dim) override; - void *extract_peratom(const char *str, int &ncol) override; - protected: int inum, maxneigh, chunk_size, chunk_offset, idx_rho_max; int host_flag; From 250307f7c612d7b9a54fa8803ea2db9508d23fed Mon Sep 17 00:00:00 2001 From: jbcouli Date: Mon, 4 Dec 2023 21:13:50 -0700 Subject: [PATCH 28/29] add description for undocumented Body entries in molecule command --- doc/src/Howto_body.rst | 2 +- doc/src/molecule.rst | 72 ++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 70 insertions(+), 4 deletions(-) diff --git a/doc/src/Howto_body.rst b/doc/src/Howto_body.rst index 115b7797c8..968e10edd8 100644 --- a/doc/src/Howto_body.rst +++ b/doc/src/Howto_body.rst @@ -335,7 +335,7 @@ faces are listed, so that M = 6 + 3\*N + 1. The integer line has three values: number of vertices (N), number of edges (E) and number of faces (F). The floating point line(s) list 6 moments of inertia followed by the coordinates of the N vertices (x1 -to zN) as 3N values, followed by 2N vertex indices corresponding to +to zN) as 3N values, followed by 2E vertex indices corresponding to the end points of the E edges, then 4\*F vertex indices defining F faces. The last value is the diameter value = the rounded diameter of the sphere that surrounds each vertex. The diameter value can be diff --git a/doc/src/molecule.rst b/doc/src/molecule.rst index 480e175e7b..1de905886c 100644 --- a/doc/src/molecule.rst +++ b/doc/src/molecule.rst @@ -154,21 +154,25 @@ These are the recognized header keywords. Header lines can come in any order. The numeric value(s) are read from the beginning of the line. The keyword should appear at the end of the line. All these settings have default values, as explained below. A line need only -appear if the value(s) are different than the default. +appear if the value(s) are different than the default, except when +defining a *body* particle, which requires setting the number of +*atoms* to 1, and setting the *inertia* in a specific section (see below). * N *atoms* = # of atoms N in molecule, default = 0 * Nb *bonds* = # of bonds Nb in molecule, default = 0 * Na *angles* = # of angles Na in molecule, default = 0 * Nd *dihedrals* = # of dihedrals Nd in molecule, default = 0 * Ni *impropers* = # of impropers Ni in molecule, default = 0 -* Nf *fragments* = # of fragments in molecule, default = 0 +* Nf *fragments* = # of fragments Nf in molecule, default = 0 +* Ninteger Ndouble *body* = # of integer and floating-point values in body +particle, default = 0 * Mtotal *mass* = total mass of molecule * Xc Yc Zc *com* = coordinates of center-of-mass of molecule * Ixx Iyy Izz Ixy Ixz Iyz *inertia* = 6 components of inertia tensor of molecule For *mass*, *com*, and *inertia*, the default is for LAMMPS to calculate this quantity itself if needed, assuming the molecules -consists of a set of point particles or finite-size particles (with a +consist of a set of point particles or finite-size particles (with a non-zero diameter) that do not overlap. If finite-size particles in the molecule do overlap, LAMMPS will not account for the overlap effects when calculating any of these 3 quantities, so you should @@ -188,6 +192,7 @@ These are the allowed section keywords for the body of the file. * *Bonds, Angles, Dihedrals, Impropers* = molecular topology sections * *Special Bond Counts, Special Bonds* = special neighbor info * *Shake Flags, Shake Atoms, Shake Bond Types* = SHAKE info +* *Body Integers, Body Doubles* = body-property sections For the Types, Bonds, Angles, Dihedrals, and Impropers sections, each atom/bond/angle/etc type can be specified either as a number (numeric @@ -515,6 +520,67 @@ of SHAKE clusters. ---------- +*Body Integers* section: + +* one line +* line syntax: N E F +* N = number of sub-particles or number or vertices +* E,F = number of edges and faces + +This section is only needed when the molecule is a body particle. the other +Body section must also appear in the file. + +The total number of values that must appear is determined by the body style, and +must be equal to the Ninteger value given in the *body* header. + +For *nparticle* and *rounded/polygon*, only the number of sub-particles or +vertices N is required, and Ninteger should have a value of 1. + +For *rounded/polyhedron*, the number of edges E and faces F is required, and +Ninteger should have a value of 3. + +See the :doc:`Howto body ` page for a further description of +the file format. + +---------- + +*Body Doubles* section: + +* first line +* line syntax: Ixx Iyy Izz Ixy Ixz Iyz +* Ixx Iyy Izz Ixy Ixz Iyz = 6 components of inertia tensor of body particle +* one line per sub-particle or vertex +* line syntax: x y z +* x, y, z = coordinates of sub-particle or vertex +* one line per edge +* line syntax: N1 N2 +* N1, N2 = vertex indices +* one line per face +* line syntax: N1 N2 N3 N4 +* N1, N2, N3, N4 = vertex indices +* last line +* line syntax: diam +* diam = rounded diameter that surrounds each vertex + +This section is only needed when the molecule is a body particle. the other +Body section must also appear in the file. + +The total number of values that must appear is determined by the body style, and +must be equal to the Ndouble value given in the *body* header. The 6 moments of +inertia and the 3N coordinates of the sub-particles or vertices are required +for all body styles. + +For *rounded/polygon*, the E = 6 + 3*N + 1 edges are automatically determined +from the vertices. + +For *rounded/polyhedron*, the 2E vertex indices for the end points of the edges +and 4F vertex indices defining the faces are required. + +See the :doc:`Howto body ` page for a further description of +the file format. + +---------- + Restrictions """""""""""" From c90999a083204d74d8e329ec92cad46ddfdd1a45 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 5 Dec 2023 21:56:29 -0500 Subject: [PATCH 29/29] cosmetic --- src/KOKKOS/pair_pace_extrapolation_kokkos.h | 1 - src/ML-PACE/pair_pace.cpp | 4 ++-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/KOKKOS/pair_pace_extrapolation_kokkos.h b/src/KOKKOS/pair_pace_extrapolation_kokkos.h index d8c0038f3b..aa6c49c36d 100644 --- a/src/KOKKOS/pair_pace_extrapolation_kokkos.h +++ b/src/KOKKOS/pair_pace_extrapolation_kokkos.h @@ -106,7 +106,6 @@ class PairPACEExtrapolationKokkos : public PairPACEExtrapolation { protected: int inum, maxneigh, chunk_size, chunk_offset, idx_ms_combs_max, total_num_functions_max; int host_flag; - //int gamma_flag; int eflag, vflag; diff --git a/src/ML-PACE/pair_pace.cpp b/src/ML-PACE/pair_pace.cpp index 933d4abecb..e9bd25f9d7 100644 --- a/src/ML-PACE/pair_pace.cpp +++ b/src/ML-PACE/pair_pace.cpp @@ -46,7 +46,7 @@ Copyright 2021 Yury Lysogorskiy^1, Cas van der Oord^2, Anton Bochkarev^1, #include "ace-evaluator/ace_recursive.h" #include "ace-evaluator/ace_version.h" #include "ace/ace_b_basis.h" -#include + namespace LAMMPS_NS { struct ACEImpl { ACEImpl() : basis_set(nullptr), ace(nullptr) {} @@ -423,4 +423,4 @@ void *PairPACE::extract_peratom(const char *str, int &ncol) } return nullptr; -} \ No newline at end of file +}