ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

updating args and D2min calculation in nonaffine fix

Browse Source
This commit is contained in:
jtclemm
2023-06-15 15:51:59 -06:00
parent ae18e1e01c
commit 25f5e74e9a
11 changed files with 1430 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
doc/src/compute.rst
View File

@ -254,6 +254,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
* :doc:`property/chunk <compute_property_chunk>` - extract various per-chunk attributes
* :doc:`property/local <compute_property_local>` - convert local attributes to localvectors/arrays
* :doc:`ptm/atom <compute_ptm_atom>` - determines the local lattice structure based on the Polyhedral Template Matching method
* :doc:`rattlers <compute_rattlers>` - identify undercoordinated rattler atoms
* :doc:`rdf <compute_rdf>` - radial distribution function g(r) histogram of group of atoms
* :doc:`reduce <compute_reduce>` - combine per-atom quantities into a single global value
* :doc:`reduce/chunk <compute_reduce_chunk>` - reduce per-atom quantities within each chunk

10
doc/src/compute_contact_atom.rst
View File

@ -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 <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

74
doc/src/compute_rattlers.rst Normal file
View File

@ -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 <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 <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 <Build_package>` page for more info.
This compute requires that atoms store a radius as defined by the
:doc:`atom_style sphere <atom_style>` command.
Related commands
""""""""""""""""
:doc:`compute coord/atom <compute_coord_atom>`
Default
"""""""
*group2-ID* = all
none

1
doc/src/fix.rst
View File

@ -259,6 +259,7 @@ accelerated styles exist.
* :doc:`mvv/tdpd <fix_mvv_dpd>` - constant temperature DPD using the modified velocity-Verlet algorithm
* :doc:`neb <fix_neb>` - nudged elastic band (NEB) spring forces
* :doc:`neb/spin <fix_neb_spin>` - nudged elastic band (NEB) spring forces for spins
* :doc:`nonaffine/displacement <fix_nonaffine_displacement>` - calculate nonaffined displacement of atoms
* :doc:`nph <fix_nh>` - constant NPH time integration via Nose/Hoover
* :doc:`nph/asphere <fix_nph_asphere>` - NPH for aspherical particles
* :doc:`nph/body <fix_nph_body>` - NPH for body particles

140
doc/src/fix_nonaffine_displacement.rst Normal file
View File

@ -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 <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 <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 <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 <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 <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 <restart>`.
The :doc:`fix_modify <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 <thermo_style>`. The default setting for
this fix is :doc:`fix_modify energy no <fix_modify>`.
The :doc:`fix_modify <fix_modify>` *respa* option is supported by this
fix. This allows to set at which level of the :doc:`r-RESPA
<run_style>` 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 <Howto_output>`. 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 <run>` command. This fix is not invoked during
:doc:`energy minimization <minimize>`.
Restrictions
""""""""""""
none
Related commands
""""""""""""""""
:doc:`atom_style sphere <atom_style>`, :doc:`fix addforce <fix_addforce>`
Default
"""""""
none

52
src/GRANULAR/compute_contact_atom.cpp
View File

@ -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;
}
}
}

4
src/GRANULAR/compute_contact_atom.h
View File

@ -37,6 +37,10 @@ class ComputeContactAtom : public Compute {
private:
int nmax;
char *group2;
int jgroup, jgroupbit;
class NeighList *list;
double *contact;
};

308
src/compute_rattlers.cpp Executable file
View File

@ -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 <cmath>
#include <cstring>
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++];
}
}

52
src/compute_rattlers.h Executable file
View File

@ -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

740
src/fix_nonaffine_displacement.cpp Executable file
View File

@ -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 <cstring>
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");
}

70
src/fix_nonaffine_displacement.h Executable file
View File

@ -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