Merge pull request #3858 from jtclemm/new-computes
Adding tools to identify rattlers and measure nonaffine displacement
This commit is contained in:
Axel Kohlmeyer
@ -115,6 +115,7 @@ KOKKOS, o = OPENMP, t = OPT.
|
||||
* :doc:`property/grid <compute_property_grid>`
|
||||
* :doc:`property/local <compute_property_local>`
|
||||
* :doc:`ptm/atom <compute_ptm_atom>`
|
||||
* :doc:`rattlers/atom <compute_rattlers_atom>`
|
||||
* :doc:`rdf <compute_rdf>`
|
||||
* :doc:`reduce <compute_reduce>`
|
||||
* :doc:`reduce/chunk <compute_reduce_chunk>`
|
||||
|
||||
@ -122,6 +122,7 @@ OPT.
|
||||
* :doc:`mvv/tdpd <fix_mvv_dpd>`
|
||||
* :doc:`neb <fix_neb>`
|
||||
* :doc:`neb/spin <fix_neb_spin>`
|
||||
* :doc:`nonaffine/displacement <fix_nonaffine_displacement>`
|
||||
* :doc:`nph (ko) <fix_nh>`
|
||||
* :doc:`nph/asphere (o) <fix_nph_asphere>`
|
||||
* :doc:`nph/body <fix_nph_body>`
|
||||
|
||||
@ -279,6 +279,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
|
||||
* :doc:`property/grid <compute_property_grid>` - convert per-grid attributes to per-grid vectors/arrays
|
||||
* :doc:`property/local <compute_property_local>` - convert local attributes to local vectors/arrays
|
||||
* :doc:`ptm/atom <compute_ptm_atom>` - determines the local lattice structure based on the Polyhedral Template Matching method
|
||||
* :doc:`rattlers/atom <compute_rattlers_atom>` - identify under-coordinated rattler atoms
|
||||
* :doc:`rdf <compute_rdf>` - radial distribution function :math:`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
|
||||
|
||||
@ -36,6 +36,9 @@ sum of the radii of the two particles.
|
||||
The value of the contact number will be 0.0 for atoms not in the
|
||||
specified compute group.
|
||||
|
||||
The optional *group2-ID* argument allows to specify from which group atoms
|
||||
contribute to the coordination number. Default setting is group 'all'.
|
||||
|
||||
Output info
|
||||
"""""""""""
|
||||
|
||||
@ -47,9 +50,6 @@ overview of LAMMPS output options.
|
||||
The per-atom vector values will be a number :math:`\ge 0.0`, as explained
|
||||
above.
|
||||
|
||||
The optional *group2-ID* argument allows to specify from which group atoms
|
||||
contribute to the coordination number. Default setting is group 'all.'
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
@ -69,6 +69,3 @@ Default
|
||||
"""""""
|
||||
|
||||
*group2-ID* = all
|
||||
|
||||
|
||||
none
|
||||
|
||||
92
doc/src/compute_rattlers_atom.rst
Normal file
92
doc/src/compute_rattlers_atom.rst
Normal file
@ -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 <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 <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 <Build_package>` page for more info.
|
||||
|
||||
The *radius* cutoff option requires that atoms store a radius as defined by the
|
||||
:doc:`atom_style sphere <atom_style>` or similar commands.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`compute coord/atom <compute_coord_atom>`
|
||||
:doc:`compute contact/atom <compute_contact_atom>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
|
||||
none
|
||||
@ -287,6 +287,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 nonaffine 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
|
||||
|
||||
133
doc/src/fix_nonaffine_displacement.rst
Normal file
133
doc/src/fix_nonaffine_displacement.rst
Normal file
@ -0,0 +1,133 @@
|
||||
.. index:: fix nonaffine/displacement
|
||||
|
||||
fix nonaffine/displacement command
|
||||
==================================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
fix ID group nonaffine/displacement style args reference/style nstep
|
||||
|
||||
* ID, group are documented in :doc:`fix <fix>` command
|
||||
* nonaffine/displacement = style name of this fix command
|
||||
* nevery = calculate nonaffine displacement every this many timesteps
|
||||
* style = *d2min* or *integrated*
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
*d2min* args = cutoff args
|
||||
cutoff = *type* or *radius* or *custom*
|
||||
*type* args = none, cutoffs determined by atom types
|
||||
*radius* args = none, cutoffs determined based on atom diameters (atom style sphere)
|
||||
*custom* args = *rmax*, cutoff set by a constant numeric value *rmax* (distance units)
|
||||
*integrated* args = none
|
||||
|
||||
* reference/style = *fixed* or *update* or *offset*
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
*fixed* = use a fixed reference frame at *nstep*
|
||||
*update* = update the reference frame every *nstep* timesteps
|
||||
*offset* = update the reference frame *nstep* timesteps before calculating the nonaffine displacement
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
fix 1 all nonaffine/displacement 100 integrated update 100
|
||||
fix 1 all nonaffine/displacement 1000 d2min type fixed 0
|
||||
fix 1 all nonaffine/displacement 1000 d2min custom 2.0 offset 100
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
.. versionadded:: TBD
|
||||
|
||||
This fix computes different metrics of the nonaffine displacement of
|
||||
particles. The first metric, *d2min* calculates the :math:`D^2_\mathrm{min}`
|
||||
nonaffine displacement by Falk and Langer in :ref:`(Falk) <d2min-Falk>`.
|
||||
For each atom, the fix computes the two tensors
|
||||
|
||||
.. math::
|
||||
|
||||
X = \sum_{\mathrm{neighbors}} \vec{r} \left(\vec{r}_{0} \right)^T
|
||||
|
||||
and
|
||||
|
||||
.. math::
|
||||
|
||||
Y = \sum_{\mathrm{neighbors}} \vec{r}_0 \left(\vec{r}_{0} \right)^T
|
||||
|
||||
where the neighbors include all other atoms within the distance criterion
|
||||
set by the cutoff option, discussed below, :math:`\vec{r}` is the current
|
||||
displacement between particles, and :math:`\vec{r}_0` is the reference
|
||||
displacement. A deformation gradient tensor is then calculated as
|
||||
:math:`F = X Y^{-1}` from which
|
||||
|
||||
.. math::
|
||||
|
||||
D^2_\mathrm{min} = \sum_{\mathrm{neighbors}} \left| \vec{r} - F \vec{r}_0 \right|^2
|
||||
|
||||
and a strain tensor is calculated :math:`E = F F^{T} - I` where :math:`I`
|
||||
is the identity tensor. This calculation is only performed on timesteps that
|
||||
are a multiple of *nevery* (including timestep zero). Data accessed before
|
||||
this occurs will simply be zeroed.
|
||||
|
||||
The *integrated* style simply integrates the velocity of particles
|
||||
every timestep to calculate a displacement. This style only works if
|
||||
used in conjunction with another fix that deforms the box and displaces
|
||||
atom positions such as :doc:`fix deform <fix_deform>` with remap x,
|
||||
:doc:`fix press/berendsen <fix_press_berendsen>`, or :doc:`fix nh <fix_nh>`.
|
||||
|
||||
Both of these methods require defining a reference state. With the *fixed* reference
|
||||
style, the user picks a specific timestep *nstep* at which particle positions are saved.
|
||||
If peratom data is accessed from this compute prior to this timestep, it will simply be
|
||||
zeroed. The *update* reference style implies the reference state will be updated every
|
||||
*nstep* timesteps. The *offset* reference only applies to the *d2min* metric and will
|
||||
update the reference state *nstep* timesteps before a multiple of *nevery* timesteps.
|
||||
|
||||
|
||||
----------
|
||||
|
||||
Restart, fix_modify, output, run start/stop, minimize info
|
||||
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
|
||||
|
||||
The reference state is saved to :doc:`binary restart files <restart>`.
|
||||
|
||||
None of the :doc:`fix_modify <fix_modify>` options are relevant to this
|
||||
fix.
|
||||
|
||||
This fix computes a peratom array with 3 columns, which can be accessed
|
||||
by indices 1-3 using any command that uses per-atom values from a fix
|
||||
as input.
|
||||
|
||||
For the *integrated* style, the three columns are the nonaffine
|
||||
displacements in the x, y, and z directions. For the *d2min* style,
|
||||
the three columns are the calculated :math:`\sqrt{D^2_\mathrm{min}}`, the
|
||||
volumetric strain, and the deviatoric strain.
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
This compute is part of the EXTRA-FIX package. It is only enabled if
|
||||
LAMMPS was built with that package. See the
|
||||
:doc:`Build package <Build_package>` page for more info.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
none
|
||||
|
||||
Default
|
||||
"""""""
|
||||
|
||||
none
|
||||
|
||||
----------
|
||||
|
||||
.. _d2min-Falk:
|
||||
|
||||
**(Falk)** Falk and Langer PRE, 57, 7192 (1998).
|
||||
@ -1074,6 +1074,7 @@ facesets
|
||||
factorizable
|
||||
factorizations
|
||||
Fahrenberger
|
||||
Falk
|
||||
Faken
|
||||
Farago
|
||||
Fasolino
|
||||
@ -1822,6 +1823,7 @@ Lanczos
|
||||
Lande
|
||||
Landron
|
||||
Landsgesell
|
||||
Langer
|
||||
langevin
|
||||
Langevin
|
||||
Langston
|
||||
@ -2499,6 +2501,7 @@ noforce
|
||||
noguess
|
||||
Noid
|
||||
nolib
|
||||
nonaffine
|
||||
nonequilibrium
|
||||
nongauss
|
||||
nonGaussian
|
||||
@ -2571,6 +2574,7 @@ nthreads
|
||||
ntimestep
|
||||
Ntptask
|
||||
Ntriples
|
||||
ntries
|
||||
ntris
|
||||
Ntype
|
||||
ntypes
|
||||
|
||||
4
src/.gitignore
vendored
4
src/.gitignore
vendored
@ -629,6 +629,8 @@
|
||||
/compute_pressure_grem.h
|
||||
/compute_ptm_atom.cpp
|
||||
/compute_ptm_atom.h
|
||||
/compute_rattlers_atom.cpp
|
||||
/compute_rattlers_atom.h
|
||||
/compute_rigid_local.cpp
|
||||
/compute_rigid_local.h
|
||||
/compute_smd_triangle_vertices.cpp
|
||||
@ -912,6 +914,8 @@
|
||||
/fix_nvt_sllod_eff.h
|
||||
/fix_nve_tri.cpp
|
||||
/fix_nve_tri.h
|
||||
/fix_nonaffine_displacement.cpp
|
||||
/fix_nonaffine_displacement.h
|
||||
/fix_oneway.cpp
|
||||
/fix_oneway.h
|
||||
/fix_orient_bcc.cpp
|
||||
|
||||
312
src/EXTRA-COMPUTE/compute_rattlers_atom.cpp
Normal file
312
src/EXTRA-COMPUTE/compute_rattlers_atom.cpp
Normal file
@ -0,0 +1,312 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing authors: Joel Clemmer (SNL), Ishan Srivastava (LBNL)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "compute_rattlers_atom.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "memory.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neighbor.h"
|
||||
#include "pair.h"
|
||||
#include "update.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
enum { TYPE, RADIUS };
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeRattlersAtom::ComputeRattlersAtom(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg), ncontacts(nullptr), rattler(nullptr)
|
||||
{
|
||||
if (narg != 6) error->all(FLERR, "Illegal compute rattlers/atom command");
|
||||
|
||||
if (strcmp(arg[3], "type") == 0)
|
||||
cutstyle = TYPE;
|
||||
else if (strcmp(arg[3], "radius") == 0)
|
||||
cutstyle = RADIUS;
|
||||
else
|
||||
error->all(FLERR, "Illegal compute rattlers/atom command");
|
||||
|
||||
if (cutstyle == RADIUS && !atom->radius_flag)
|
||||
error->all(FLERR, "Compute rattlers/atom radius style requires atom attribute radius");
|
||||
|
||||
ncontacts_rattler = utils::inumeric(FLERR, arg[4], false, lmp);
|
||||
max_tries = utils::inumeric(FLERR, arg[5], false, lmp);
|
||||
|
||||
nmax = 0;
|
||||
invoked_peratom = -1;
|
||||
|
||||
scalar_flag = 1;
|
||||
extscalar = 1;
|
||||
peratom_flag = 1;
|
||||
size_peratom_cols = 0;
|
||||
comm_forward = 1;
|
||||
comm_reverse = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeRattlersAtom::~ComputeRattlersAtom()
|
||||
{
|
||||
memory->destroy(ncontacts);
|
||||
memory->destroy(rattler);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeRattlersAtom::init()
|
||||
{
|
||||
if (force->pair == nullptr) error->all(FLERR, "No pair style is defined for compute rattlers");
|
||||
|
||||
// Cannot calculate distance from radii for JKR/DMT
|
||||
if (force->pair->beyond_contact)
|
||||
error->all(FLERR, "Compute rattlers does not currently support pair styles that extend beyond contact");
|
||||
|
||||
// need an occasional half neighbor list
|
||||
// set size to same value as request made by force->pair
|
||||
// this should enable it to always be a copy list (e.g. for granular pstyle)
|
||||
|
||||
auto pairrequest = neighbor->find_request(force->pair);
|
||||
if (pairrequest && pairrequest->get_size())
|
||||
neighbor->add_request(this, NeighConst::REQ_SIZE | NeighConst::REQ_OCCASIONAL);
|
||||
else
|
||||
neighbor->add_request(this, NeighConst::REQ_OCCASIONAL);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeRattlersAtom::init_list(int /*id*/, NeighList *ptr)
|
||||
{
|
||||
list = ptr;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeRattlersAtom::compute_peratom()
|
||||
{
|
||||
if (invoked_peratom == update->ntimestep) return;
|
||||
invoked_peratom = update->ntimestep;
|
||||
|
||||
int i, j, ii, jj, inum, jnum, itype, jtype, tmp_flag;
|
||||
tagint itag, jtag;
|
||||
double xtmp, ytmp, ztmp, delx, dely, delz;
|
||||
double rsq, radsum;
|
||||
|
||||
if (nmax < atom->nmax) {
|
||||
nmax = atom->nmax;
|
||||
memory->destroy(ncontacts);
|
||||
memory->destroy(rattler);
|
||||
memory->create(ncontacts, nmax, "rattlers:ncontacts");
|
||||
memory->create(rattler, nmax, "rattlers:rattler");
|
||||
vector_atom = rattler;
|
||||
}
|
||||
|
||||
for (i = 0; i < nmax; i++) rattler[i] = 0;
|
||||
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
|
||||
double **x = atom->x;
|
||||
double *radius = atom->radius;
|
||||
tagint *tag = atom->tag;
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
int newton_pair = force->newton_pair;
|
||||
|
||||
// invoke half neighbor list (will copy or build if necessary)
|
||||
neighbor->build_one(list);
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
Pair *pair = force->pair;
|
||||
double **cutsq = force->pair->cutsq;
|
||||
|
||||
int change_flag = 1;
|
||||
int ntry = 0;
|
||||
while (ntry < max_tries) {
|
||||
change_flag = 0;
|
||||
|
||||
for (i = 0; i < nmax; i++) ncontacts[i] = 0;
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
if (!(mask[i] & groupbit)) continue;
|
||||
if (rattler[i] == 1) continue;
|
||||
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itag = tag[i];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
if (!(mask[j] & groupbit)) continue;
|
||||
if (rattler[j] == 1) continue;
|
||||
|
||||
// itag = jtag is possible for long cutoffs that include images of self
|
||||
|
||||
if (newton_pair == 0 && j >= nlocal) {
|
||||
jtag = tag[j];
|
||||
if (itag > jtag) {
|
||||
if ((itag + jtag) % 2 == 0) continue;
|
||||
} else if (itag < jtag) {
|
||||
if ((itag + jtag) % 2 == 1) continue;
|
||||
} else {
|
||||
if (x[j][2] < ztmp) continue;
|
||||
if (x[j][2] == ztmp) {
|
||||
if (x[j][1] < ytmp) continue;
|
||||
if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
jtype = type[j];
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
|
||||
if (cutstyle == TYPE) {
|
||||
if (rsq >= cutsq[itype][jtype]) continue;
|
||||
} else {
|
||||
radsum = radius[i] + radius[j];
|
||||
if (rsq >= radsum * radsum) continue;
|
||||
}
|
||||
ncontacts[i] += 1;
|
||||
if (newton_pair || j < nlocal)
|
||||
ncontacts[j] += 1;
|
||||
}
|
||||
}
|
||||
|
||||
// add contributions from ghosts
|
||||
if (force->newton_pair) comm->reverse_comm(this);
|
||||
|
||||
// Set flags for rattlers
|
||||
for (i = 0; i < atom->nlocal; i++) {
|
||||
if (ncontacts[i] < ncontacts_rattler && rattler[i] == 0) {
|
||||
rattler[i] = 1;
|
||||
change_flag = 1;
|
||||
}
|
||||
}
|
||||
|
||||
comm->forward_comm(this);
|
||||
|
||||
MPI_Allreduce(&change_flag, &tmp_flag, 1, MPI_INT, MPI_MAX, world);
|
||||
change_flag = tmp_flag;
|
||||
if (change_flag == 0) break;
|
||||
|
||||
ntry += 1;
|
||||
}
|
||||
|
||||
if (change_flag == 1)
|
||||
error->warning(FLERR, "Rattler calculation failed to converge within max tries");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double ComputeRattlersAtom::compute_scalar()
|
||||
{
|
||||
if (invoked_peratom != update->ntimestep)
|
||||
compute_peratom();
|
||||
|
||||
invoked_scalar = update->ntimestep;
|
||||
|
||||
double total_rattlers = 0;
|
||||
for (int i = 0; i < atom->nlocal; i++) {
|
||||
if (rattler[i] == 1) {
|
||||
total_rattlers += 1;
|
||||
}
|
||||
}
|
||||
|
||||
//Total across processors
|
||||
MPI_Allreduce(&total_rattlers, &scalar, 1, MPI_DOUBLE, MPI_SUM, world);
|
||||
return scalar;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int ComputeRattlersAtom::pack_reverse_comm(int n, int first, double *buf)
|
||||
{
|
||||
int i, m, last;
|
||||
|
||||
m = 0;
|
||||
last = first + n;
|
||||
for (i = first; i < last; i++) {
|
||||
buf[m++] = ubuf(ncontacts[i]).d;
|
||||
}
|
||||
return m;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeRattlersAtom::unpack_reverse_comm(int n, int *list, double *buf)
|
||||
{
|
||||
int i, j, m;
|
||||
|
||||
m = 0;
|
||||
for (i = 0; i < n; i++) {
|
||||
j = list[i];
|
||||
ncontacts[j] += (int) ubuf(buf[m++]).i;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int ComputeRattlersAtom::pack_forward_comm(int n, int *list, double *buf,
|
||||
int /*pbc_flag*/, int * /*pbc*/)
|
||||
{
|
||||
int i, j, m;
|
||||
|
||||
m = 0;
|
||||
for (i = 0; i < n; i++) {
|
||||
j = list[i];
|
||||
buf[m++] = rattler[j];
|
||||
}
|
||||
|
||||
return m;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeRattlersAtom::unpack_forward_comm(int n, int first, double *buf)
|
||||
{
|
||||
int i, m, last;
|
||||
|
||||
m = 0;
|
||||
last = first + n;
|
||||
for (i = first; i < last; i++) {
|
||||
rattler[i] = buf[m++];
|
||||
}
|
||||
}
|
||||
52
src/EXTRA-COMPUTE/compute_rattlers_atom.h
Normal file
52
src/EXTRA-COMPUTE/compute_rattlers_atom.h
Normal file
@ -0,0 +1,52 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef COMPUTE_CLASS
|
||||
// clang-format off
|
||||
ComputeStyle(rattlers/atom,ComputeRattlersAtom);
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
#ifndef LMP_COMPUTE_RATTLERS_ATOM_H
|
||||
#define LMP_COMPUTE_RATTLERS_ATOM_H
|
||||
|
||||
#include "compute.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class ComputeRattlersAtom : public Compute {
|
||||
public:
|
||||
ComputeRattlersAtom(class LAMMPS *, int, char **);
|
||||
~ComputeRattlersAtom() override;
|
||||
void init() override;
|
||||
void init_list(int, class NeighList *) override;
|
||||
void compute_peratom() override;
|
||||
double compute_scalar() override;
|
||||
int pack_forward_comm(int, int *, double *, int, int *) override;
|
||||
void unpack_forward_comm(int, int, double *) override;
|
||||
int pack_reverse_comm(int, int, double *) override;
|
||||
void unpack_reverse_comm(int, int *, double *) override;
|
||||
|
||||
private:
|
||||
int pstyle, cutstyle;
|
||||
int ncontacts_rattler, max_tries, nmax, invoked_peratom;
|
||||
int *ncontacts;
|
||||
double *rattler;
|
||||
class NeighList *list;
|
||||
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
736
src/EXTRA-FIX/fix_nonaffine_displacement.cpp
Normal file
736
src/EXTRA-FIX/fix_nonaffine_displacement.cpp
Normal file
@ -0,0 +1,736 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing authors: Joel Clemmer (SNL), Ishan Srivastava (LBNL)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "fix_nonaffine_displacement.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "citeme.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "fix_store_atom.h"
|
||||
#include "force.h"
|
||||
#include "group.h"
|
||||
#include "math_extra.h"
|
||||
#include "memory.h"
|
||||
#include "modify.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neighbor.h"
|
||||
#include "pair.h"
|
||||
#include "update.h"
|
||||
|
||||
#include <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), id_fix(nullptr), X(nullptr), Y(nullptr), F(nullptr), norm(nullptr)
|
||||
{
|
||||
if (narg < 4) error->all(FLERR,"Illegal fix nonaffine/displacement command");
|
||||
|
||||
nevery = utils::inumeric(FLERR, arg[3], false, lmp);
|
||||
if (nevery <= 0) error->all(FLERR,"Illegal nevery value {} in fix nonaffine/displacement", nevery);
|
||||
|
||||
int iarg = 4;
|
||||
if (strcmp(arg[iarg], "integrated") == 0) {
|
||||
nad_style = INTEGRATED;
|
||||
nevery = 1;
|
||||
iarg += 1;
|
||||
} else if (strcmp(arg[iarg], "d2min") == 0) {
|
||||
if (iarg + 1 > narg) error->all(FLERR,"Illegal fix nonaffine/displacement command");
|
||||
nad_style = D2MIN;
|
||||
if (strcmp(arg[iarg + 1], "type") == 0) {
|
||||
cut_style = TYPE;
|
||||
} else if (strcmp(arg[iarg + 1], "radius") == 0) {
|
||||
cut_style = RADIUS;
|
||||
} else if (strcmp(arg[iarg + 1], "custom") == 0) {
|
||||
if (iarg + 2 > narg) error->all(FLERR,"Illegal fix nonaffine/displacement command");
|
||||
cut_style = CUSTOM;
|
||||
cutoff_custom = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
|
||||
cutsq_custom = cutoff_custom * cutoff_custom;
|
||||
if (cutoff_custom <= 0)
|
||||
error->all(FLERR, "Illegal custom cutoff length {}", arg[iarg + 2]);
|
||||
iarg += 1;
|
||||
} else error->all(FLERR,"Illegal cutoff style {} in fix nonaffine/displacement", arg[iarg + 1]);
|
||||
iarg += 2;
|
||||
} else error->all(FLERR,"Illegal nonaffine displacement style {} in fix nonaffine/displacement", arg[iarg]);
|
||||
|
||||
if (iarg + 2 > narg) error->all(FLERR,"Illegal fix nonaffine/displacement command");
|
||||
if (strcmp(arg[iarg], "fixed") == 0) {
|
||||
reference_style = FIXED;
|
||||
reference_timestep = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
|
||||
if (update_timestep < 0)
|
||||
error->all(FLERR, "Illegal reference timestep {} in fix nonaffine/displacement", arg[iarg + 1]);
|
||||
} else if (strcmp(arg[iarg], "update") == 0) {
|
||||
reference_style = UPDATE;
|
||||
update_timestep = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
|
||||
if (update_timestep < 0)
|
||||
error->all(FLERR, "Illegal update timestep {} in fix nonaffine/displacement", arg[iarg + 1]);
|
||||
} else if (strcmp(arg[iarg], "offset") == 0) {
|
||||
reference_style = OFFSET;
|
||||
offset_timestep = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
|
||||
if ((offset_timestep <= 0) || (offset_timestep > nevery))
|
||||
error->all(FLERR, "Illegal offset timestep {} in fix nonaffine/displacement", arg[iarg + 1]);
|
||||
} else error->all(FLERR,"Illegal reference style {} in fix nonaffine/displacement", arg[iarg]);
|
||||
|
||||
if (nad_style == D2MIN)
|
||||
if (cut_style == RADIUS && (!atom->radius_flag))
|
||||
error->all(FLERR, "Fix nonaffine/displacement radius style requires atom attribute radius");
|
||||
|
||||
if (nad_style == INTEGRATED && reference_style == OFFSET)
|
||||
error->all(FLERR, "Fix nonaffine/displacement cannot use the integrated style with an offset reference state");
|
||||
|
||||
peratom_flag = 1;
|
||||
peratom_freq = nevery;
|
||||
nmax = -1;
|
||||
reference_saved = 0;
|
||||
restart_global = 1;
|
||||
|
||||
size_peratom_cols = 3;
|
||||
comm_reverse = 0;
|
||||
comm_forward = 0;
|
||||
if (nad_style == D2MIN) {
|
||||
comm_reverse = 18;
|
||||
comm_forward = 9;
|
||||
}
|
||||
|
||||
if (nad_style == D2MIN && lmp->citeme) lmp->citeme->add(cite_nonaffine_d2min);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixNonaffineDisplacement::~FixNonaffineDisplacement()
|
||||
{
|
||||
if (id_fix && modify->nfix) modify->delete_fix(id_fix);
|
||||
delete[] id_fix;
|
||||
|
||||
if (nad_style == D2MIN) {
|
||||
memory->destroy(X);
|
||||
memory->destroy(Y);
|
||||
memory->destroy(F);
|
||||
memory->destroy(norm);
|
||||
memory->destroy(array_atom);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixNonaffineDisplacement::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= POST_FORCE;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNonaffineDisplacement::post_constructor()
|
||||
{
|
||||
// Create persistent peratom storage for either an integrated velocity or reference position
|
||||
// Ghost atoms need reference coordinates for D2min
|
||||
std::string ghost_status = "0";
|
||||
if (nad_style == D2MIN) ghost_status = "1";
|
||||
|
||||
id_fix = utils::strdup(id + std::string("_FIX_PA"));
|
||||
fix = dynamic_cast<FixStoreAtom *>(modify->add_fix(fmt::format("{} {} STORE/ATOM 3 0 {} 1", id_fix, group->names[igroup], ghost_status)));
|
||||
|
||||
if (nad_style == INTEGRATED)
|
||||
array_atom = fix->astore;
|
||||
|
||||
if (nad_style == D2MIN)
|
||||
grow_arrays(atom->nmax);
|
||||
|
||||
for (int i = 0; i < atom->nlocal; i++)
|
||||
for (int j = 0; j < 3; j++) array_atom[i][j] = 0.0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNonaffineDisplacement::init()
|
||||
{
|
||||
dtv = update->dt;
|
||||
|
||||
if ((!reference_saved) && (reference_style == FIXED) && (update->ntimestep > reference_timestep))
|
||||
error->all(FLERR, "Initial timestep exceeds that of the reference state in fix nonaffine/displacement");
|
||||
|
||||
if (nad_style == D2MIN) {
|
||||
if ((!force->pair) && (cut_style == TYPE))
|
||||
error->all(FLERR,"Fix nonaffine/displacement D2Min option requires a pair style be defined "
|
||||
"or cutoff specified");
|
||||
|
||||
// need an occasional half neighbor list
|
||||
|
||||
if (cut_style == RADIUS) {
|
||||
auto req = neighbor->add_request(this, NeighConst::REQ_SIZE | NeighConst::REQ_OCCASIONAL);
|
||||
} else {
|
||||
auto req = neighbor->add_request(this, NeighConst::REQ_OCCASIONAL);
|
||||
if (cut_style == CUSTOM) {
|
||||
double skin = neighbor->skin;
|
||||
mycutneigh = cutoff_custom + skin;
|
||||
|
||||
double cutghost; // as computed by Neighbor and Comm
|
||||
if (force->pair)
|
||||
cutghost = MAX(force->pair->cutforce + skin, comm->cutghostuser);
|
||||
else
|
||||
cutghost = comm->cutghostuser;
|
||||
|
||||
if (mycutneigh > cutghost)
|
||||
error->all(FLERR,"Fix nonaffine/displacement D2Min option cutoff exceeds ghost atom range - use comm_modify cutoff command");
|
||||
|
||||
req->set_cutoff(mycutneigh);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNonaffineDisplacement::init_list(int /*id*/, NeighList *ptr)
|
||||
{
|
||||
list = ptr;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNonaffineDisplacement::setup(int vflag)
|
||||
{
|
||||
post_force(0); // Save state if needed before starting the 1st timestep
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNonaffineDisplacement::post_force(int /*vflag*/)
|
||||
{
|
||||
if (reference_saved && (!update->setupflag)) {
|
||||
if (nad_style == INTEGRATED) {
|
||||
integrate_velocity();
|
||||
} else {
|
||||
if ((update->ntimestep % nevery) == 0) calculate_D2Min();
|
||||
}
|
||||
}
|
||||
|
||||
if (reference_style == FIXED)
|
||||
if (update->ntimestep == reference_timestep)
|
||||
save_reference_state();
|
||||
|
||||
if (reference_style == UPDATE)
|
||||
if ((update->ntimestep % update_timestep) == 0)
|
||||
save_reference_state();
|
||||
|
||||
if (reference_style == OFFSET)
|
||||
if (((update->ntimestep + offset_timestep) % nevery) == 0)
|
||||
save_reference_state();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNonaffineDisplacement::write_restart(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
int size = sizeof(int);
|
||||
fwrite(&size, sizeof(int), 1, fp);
|
||||
fwrite(&reference_saved, sizeof(int), 1, fp);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNonaffineDisplacement::restart(char *buf)
|
||||
{
|
||||
reference_saved = (int) ubuf(buf[0]).i;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNonaffineDisplacement::integrate_velocity()
|
||||
{
|
||||
int i,n;
|
||||
dtv = update->dt;
|
||||
|
||||
double **v = atom->v;
|
||||
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (int m = 0; m < 3; m++) {
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
array_atom[i][m] += dtv * v[i][m];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNonaffineDisplacement::save_reference_state()
|
||||
{
|
||||
int i, n;
|
||||
double **x = atom->x;
|
||||
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
int nall = nlocal + atom->nghost;
|
||||
|
||||
if (nad_style == D2MIN) {
|
||||
for (int m = 0; m < 3; m++) {
|
||||
for (int i = 0; i < nall; i++) {
|
||||
if (mask[i] & groupbit) array_atom[i][m] = x[i][m];
|
||||
}
|
||||
}
|
||||
} else {
|
||||
for (int m = 0; m < 3; m++) {
|
||||
for (int i = 0; i < nall; i++) {
|
||||
if (mask[i] & groupbit) array_atom[i][m] = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (nad_style == D2MIN) {
|
||||
xprd0 = domain->xprd;
|
||||
yprd0 = domain->yprd;
|
||||
zprd0 = domain->zprd;
|
||||
xprd0_half = domain->xprd_half;
|
||||
yprd0_half = domain->yprd_half;
|
||||
zprd0_half = domain->zprd_half;
|
||||
xy0 = domain->xy;
|
||||
xz0 = domain->xz;
|
||||
yz0 = domain->yz;
|
||||
}
|
||||
|
||||
reference_saved = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNonaffineDisplacement::calculate_D2Min()
|
||||
{
|
||||
// invoke half neighbor list (will copy or build if necessary)
|
||||
neighbor->build_one(list);
|
||||
|
||||
if (atom->nmax > nmax)
|
||||
grow_arrays(atom->nmax);
|
||||
|
||||
int i, j, k, l, ii, jj, inum, jnum, itype, jtype;
|
||||
double evol, j2, edev;
|
||||
double r[3], r0[3], rsq, rsq0, radsum, temp[3];
|
||||
double X_tmp[3][3], Y_tmp[3][3], F_tmp[3][3], E[3][3];
|
||||
double Y_inv[3][3] = {0.0}; // Zero for 2d since not all entries used
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
|
||||
double **x = atom->x;
|
||||
double **x0 = array_atom;
|
||||
double *radius = atom->radius;
|
||||
tagint *tag = atom->tag;
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
int newton_pair = force->newton_pair;
|
||||
int dim = domain->dimension;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
Pair *pair = force->pair;
|
||||
double **cutsq;
|
||||
if (pair) cutsq = force->pair->cutsq;
|
||||
|
||||
for (i = 0; i < nmax; i++) {
|
||||
for (k = 0; k < 3; k++) {
|
||||
for (l = 0; l < 3; l++) {
|
||||
X[i][k][l] = 0.0;
|
||||
Y[i][k][l] = 0.0;
|
||||
}
|
||||
}
|
||||
norm[i] = 0;
|
||||
array_atom[i][0] = 0;
|
||||
}
|
||||
|
||||
// First loop through neighbors
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
if (!(mask[i] & groupbit)) continue;
|
||||
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
if (!(mask[j] & groupbit)) continue;
|
||||
|
||||
jtype = type[j];
|
||||
r[0] = x[i][0] - x[j][0];
|
||||
r[1] = x[i][1] - x[j][1];
|
||||
r[2] = x[i][2] - x[j][2];
|
||||
rsq = lensq3(r);
|
||||
|
||||
// Only include contributions from atoms that are CURRENTLY neighbors
|
||||
if (cut_style == TYPE) {
|
||||
if (rsq > cutsq[itype][jtype]) continue;
|
||||
} else if (cut_style == CUSTOM) {
|
||||
if (rsq > cutsq_custom) continue;
|
||||
} else {
|
||||
radsum = radius[i] + radius[j];
|
||||
if (rsq > (radsum * radsum)) continue;
|
||||
}
|
||||
|
||||
r0[0] = x0[i][0] - x0[j][0];
|
||||
r0[1] = x0[i][1] - x0[j][1];
|
||||
r0[2] = x0[i][2] - x0[j][2];
|
||||
minimum_image0(r0);
|
||||
|
||||
// Using notation from Falk & Langer 1998
|
||||
outer3(r, r0, X_tmp);
|
||||
outer3(r0, r0, Y_tmp);
|
||||
|
||||
for (k = 0; k < 3; k++) {
|
||||
for (l = 0; l < 3; l++) {
|
||||
X[i][k][l] += X_tmp[k][l];
|
||||
Y[i][k][l] += Y_tmp[k][l];
|
||||
}
|
||||
}
|
||||
|
||||
if (newton_pair || j < nlocal) {
|
||||
for (k = 0; k < 3; k++) {
|
||||
for (l = 0; l < 3; l++) {
|
||||
X[j][k][l] += X_tmp[k][l];
|
||||
Y[j][k][l] += Y_tmp[k][l];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
comm_flag = 0;
|
||||
if (newton_pair) comm->reverse_comm(this, 18);
|
||||
|
||||
// Calculate contributions to strain tensor
|
||||
double denom;
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
if (!(mask[i] & groupbit)) continue;
|
||||
for (j = 0; j < 3; j++) {
|
||||
for (k = 0; k < 3; k++) {
|
||||
Y_tmp[j][k] = Y[i][j][k];
|
||||
X_tmp[j][k] = X[i][j][k];
|
||||
}
|
||||
}
|
||||
|
||||
if (dim == 3) {
|
||||
invert3(Y_tmp, Y_inv);
|
||||
} else {
|
||||
denom = Y_tmp[0][0] * Y_tmp[1][1] - Y_tmp[0][1] * Y_tmp[1][0];
|
||||
if (denom != 0.0) denom = 1.0 / denom;
|
||||
Y_inv[0][0] = Y_tmp[1][1] * denom;
|
||||
Y_inv[0][1] = -Y_tmp[0][1] * denom;
|
||||
Y_inv[1][0] = -Y_tmp[1][0] * denom;
|
||||
Y_inv[1][1] = Y_tmp[0][0] * denom;
|
||||
}
|
||||
|
||||
times3(X_tmp, Y_inv, F_tmp);
|
||||
|
||||
for (j = 0; j < 3; j++) {
|
||||
for (k = 0; k < 3; k++) {
|
||||
F[i][j][k] = F_tmp[j][k];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
comm->forward_comm(this);
|
||||
|
||||
// Second loop through neighbors
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
if (!(mask[i] & groupbit)) continue;
|
||||
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
if (!(mask[j] & groupbit)) continue;
|
||||
|
||||
jtype = type[j];
|
||||
r[0] = x[i][0] - x[j][0];
|
||||
r[1] = x[i][1] - x[j][1];
|
||||
r[2] = x[i][2] - x[j][2];
|
||||
rsq = lensq3(r);
|
||||
|
||||
// Only include contributions from atoms that are CURRENTLY neighbors
|
||||
if (cut_style == TYPE) {
|
||||
if (rsq >= cutsq[itype][jtype]) continue;
|
||||
} else if (cut_style == CUSTOM) {
|
||||
if (rsq >= cutsq_custom) continue;
|
||||
} else {
|
||||
radsum = radius[i] + radius[j];
|
||||
if (rsq >= radsum * radsum) continue;
|
||||
}
|
||||
|
||||
r0[0] = x0[i][0] - x0[j][0];
|
||||
r0[1] = x0[i][1] - x0[j][1];
|
||||
r0[2] = x0[i][2] - x0[j][2];
|
||||
minimum_image0(r0);
|
||||
|
||||
// E * r0
|
||||
for (k = 0; k < 3; k++) {
|
||||
temp[k] = 0.0;
|
||||
for (l = 0; l < 3; l++)
|
||||
temp[k] += F[i][k][l] * r0[l];
|
||||
}
|
||||
|
||||
sub3(r, temp, temp);
|
||||
array_atom[i][0] += lensq3(temp);
|
||||
norm[i] += 1;
|
||||
|
||||
if (newton_pair || j < nlocal) {
|
||||
for (k = 0; k < 3; k++) {
|
||||
temp[k] = 0.0;
|
||||
for (l = 0; l < 3; l++)
|
||||
temp[k] += F[j][k][l] * r0[l];
|
||||
}
|
||||
|
||||
sub3(r, temp, temp);
|
||||
array_atom[j][0] += lensq3(temp);
|
||||
norm[j] += 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
comm_flag = 1;
|
||||
if (newton_pair) comm->reverse_comm(this, 2);
|
||||
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
if (!(mask[i] & groupbit)) continue;
|
||||
|
||||
if (norm[i] != 0)
|
||||
array_atom[i][0] /= norm[i];
|
||||
else
|
||||
array_atom[i][0] = 0.0;
|
||||
array_atom[i][0] = sqrt(array_atom[i][0]);
|
||||
|
||||
for (j = 0; j < 3; j++)
|
||||
for (k = 0; k < 3; k++)
|
||||
F_tmp[j][k] = F[i][j][k];
|
||||
|
||||
transpose_times3(F_tmp, F_tmp, E);
|
||||
for (j = 0; j < dim; j++) E[j][j] -= 1.0;
|
||||
|
||||
evol = (E[0][0] + E[1][1] + E[2][2]) / dim;
|
||||
|
||||
// Calculate deviatoric strain
|
||||
for (j = 0; j < dim; j++) E[j][j] -= evol;
|
||||
j2 = 0.0;
|
||||
for (j = 0; j < 3; j++)
|
||||
for (k = 0; k < 3; k++)
|
||||
j2 += E[j][k] * E[j][k];
|
||||
|
||||
edev = sqrt(0.5 * j2);
|
||||
|
||||
array_atom[i][1] = evol;
|
||||
array_atom[i][2] = edev;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixNonaffineDisplacement::pack_reverse_comm(int n, int first, double *buf)
|
||||
{
|
||||
int i, m, last, k, l;
|
||||
|
||||
m = 0;
|
||||
last = first + n;
|
||||
for (i = first; i < last; i++) {
|
||||
if (comm_flag == 0) {
|
||||
for (k = 0; k < 3; k++) {
|
||||
for (l = 0; l < 3; l++) {
|
||||
buf[m++] = X[i][k][l];
|
||||
buf[m++] = Y[i][k][l];
|
||||
}
|
||||
}
|
||||
} else {
|
||||
buf[m++] = array_atom[i][0];
|
||||
buf[m++] = ubuf(norm[i]).d;
|
||||
}
|
||||
}
|
||||
return m;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNonaffineDisplacement::unpack_reverse_comm(int n, int *list, double *buf)
|
||||
{
|
||||
int i, j, m, k, l;
|
||||
|
||||
m = 0;
|
||||
for (i = 0; i < n; i++) {
|
||||
j = list[i];
|
||||
if (comm_flag == 0) {
|
||||
for (k = 0; k < 3; k++) {
|
||||
for (l = 0; l < 3; l++) {
|
||||
X[j][k][l] += buf[m++];
|
||||
Y[j][k][l] += buf[m++];
|
||||
}
|
||||
}
|
||||
} else {
|
||||
array_atom[j][0] += buf[m++];
|
||||
norm[j] += (int) ubuf(buf[m++]).i;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixNonaffineDisplacement::pack_forward_comm(int n, int *list, double *buf,
|
||||
int /*pbc_flag*/, int * /*pbc*/)
|
||||
{
|
||||
int i, j, m, k, l;
|
||||
|
||||
m = 0;
|
||||
for (i = 0; i < n; i++) {
|
||||
j = list[i];
|
||||
for (k = 0; k < 3; k++) {
|
||||
for (l = 0; l < 3; l ++) {
|
||||
buf[m++] = F[j][k][l];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return m;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNonaffineDisplacement::unpack_forward_comm(int n, int first, double *buf)
|
||||
{
|
||||
int i, m, last, k, l;
|
||||
|
||||
m = 0;
|
||||
last = first + n;
|
||||
for (i = first; i < last; i++) {
|
||||
for (k = 0; k < 3; k++) {
|
||||
for (l = 0; l < 3; l ++) {
|
||||
F[i][k][l] = buf[m++];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNonaffineDisplacement::minimum_image0(double *delta)
|
||||
{
|
||||
if (domain->triclinic == 0) {
|
||||
if (domain->xperiodic) {
|
||||
while (fabs(delta[0]) > xprd0_half) {
|
||||
if (delta[0] < 0.0) delta[0] += xprd0;
|
||||
else delta[0] -= xprd0;
|
||||
}
|
||||
}
|
||||
if (domain->yperiodic) {
|
||||
while (fabs(delta[1]) > yprd0_half) {
|
||||
if (delta[1] < 0.0) delta[1] += yprd0;
|
||||
else delta[1] -= yprd0;
|
||||
}
|
||||
}
|
||||
if (domain->zperiodic) {
|
||||
while (fabs(delta[2]) > zprd0_half) {
|
||||
if (delta[2] < 0.0) delta[2] += zprd0;
|
||||
else delta[2] -= zprd0;
|
||||
}
|
||||
}
|
||||
|
||||
} else {
|
||||
if (domain->zperiodic) {
|
||||
while (fabs(delta[2]) > zprd0_half) {
|
||||
if (delta[2] < 0.0) {
|
||||
delta[2] += zprd0;
|
||||
delta[1] += yz0;
|
||||
delta[0] += xz0;
|
||||
} else {
|
||||
delta[2] -= zprd0;
|
||||
delta[1] -= yz0;
|
||||
delta[0] -= xz0;
|
||||
}
|
||||
}
|
||||
}
|
||||
if (domain->yperiodic) {
|
||||
while (fabs(delta[1]) > yprd0_half) {
|
||||
if (delta[1] < 0.0) {
|
||||
delta[1] += yprd0;
|
||||
delta[0] += xy0;
|
||||
} else {
|
||||
delta[1] -= yprd0;
|
||||
delta[0] -= xy0;
|
||||
}
|
||||
}
|
||||
}
|
||||
if (domain->xperiodic) {
|
||||
while (fabs(delta[0]) > xprd0_half) {
|
||||
if (delta[0] < 0.0) delta[0] += xprd0;
|
||||
else delta[0] -= xprd0;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixNonaffineDisplacement::grow_arrays(int nmax_new)
|
||||
{
|
||||
nmax = nmax_new;
|
||||
memory->destroy(X);
|
||||
memory->destroy(Y);
|
||||
memory->destroy(F);
|
||||
memory->destroy(norm);
|
||||
memory->create(X, nmax, 3, 3, "fix_nonaffine_displacement:X");
|
||||
memory->create(Y, nmax, 3, 3, "fix_nonaffine_displacement:Y");
|
||||
memory->create(F, nmax, 3, 3, "fix_nonaffine_displacement:F");
|
||||
memory->create(norm, nmax, "fix_nonaffine_displacement:norm");
|
||||
}
|
||||
71
src/EXTRA-FIX/fix_nonaffine_displacement.h
Normal file
71
src/EXTRA-FIX/fix_nonaffine_displacement.h
Normal file
@ -0,0 +1,71 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
// clang-format off
|
||||
FixStyle(nonaffine/displacement,FixNonaffineDisplacement)
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_NONAFFINE_DISPLACEMENT_H
|
||||
#define LMP_FIX_NONAFFINE_DISPLACEMENT_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixNonaffineDisplacement : public Fix {
|
||||
public:
|
||||
FixNonaffineDisplacement(class LAMMPS *, int, char **);
|
||||
~FixNonaffineDisplacement() override;
|
||||
int setmask() override;
|
||||
void post_constructor() override;
|
||||
void init() override;
|
||||
void init_list(int, class NeighList *) override;
|
||||
void setup(int);
|
||||
void post_force(int) override;
|
||||
void write_restart(FILE *fp) override;
|
||||
void restart(char *buf) override;
|
||||
int pack_forward_comm(int, int *, double *, int, int *) override;
|
||||
void unpack_forward_comm(int, int, double *) override;
|
||||
int pack_reverse_comm(int, int, double *) override;
|
||||
void unpack_reverse_comm(int, int *, double *) override;
|
||||
|
||||
private:
|
||||
double dtv;
|
||||
char *id_fix;
|
||||
class FixStoreAtom *fix;
|
||||
int nmax, comm_flag;
|
||||
int nad_style, cut_style;
|
||||
int reference_style, offset_timestep, reference_timestep, update_timestep;
|
||||
int reference_saved;
|
||||
double cutoff_custom, cutsq_custom, mycutneigh;
|
||||
double xprd0, yprd0, zprd0, xprd0_half, yprd0_half, zprd0_half, xy0, xz0, yz0;
|
||||
|
||||
double ***X, ***Y, ***F;
|
||||
int *norm;
|
||||
|
||||
class NeighList *list; // half neighbor list
|
||||
|
||||
|
||||
void integrate_velocity();
|
||||
void calculate_D2Min();
|
||||
void save_reference_state();
|
||||
void minimum_image0(double *);
|
||||
void grow_arrays(int);
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
@ -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"
|
||||
|
||||
Reference in New Issue
Block a user