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

Merge pull request #2481 from lammps/energy_virial_constants

Browse Source 
Introduce enums for energy and virial flags
This commit is contained in:
Richard Berger
2020-11-23 10:35:20 -05:00
committed by GitHub
parent d1af1aa12d e1835250c7
commit b8d821e1f9
107 changed files with 576 additions and 397 deletions
Download Patch File Download Diff File Expand all files Collapse all files

236
doc/src/compute_stress_atom.rst
View File

@ -33,32 +33,31 @@ Examples
Description
"""""""""""
Define a computation that computes per-atom stress
tensor for each atom in a group. In case of compute *stress/atom*,
the tensor for each atom is symmetric with 6
components and is stored as a 6-element vector in the following order:
:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`.
In case of compute *centroid/stress/atom*,
the tensor for each atom is asymmetric with 9 components
and is stored as a 9-element vector in the following order:
:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`,
:math:`yx`, :math:`zx`, :math:`zy`.
See the :doc:`compute pressure <compute_pressure>` command if you want the stress tensor
Define a computation that computes per-atom stress tensor for each
atom in a group. In case of compute *stress/atom*, the tensor for
each atom is symmetric with 6 components and is stored as a 6-element
vector in the following order: :math:`xx`, :math:`yy`, :math:`zz`,
:math:`xy`, :math:`xz`, :math:`yz`. In case of compute
*centroid/stress/atom*, the tensor for each atom is asymmetric with 9
components and is stored as a 9-element vector in the following order:
:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`,
:math:`yz`, :math:`yx`, :math:`zx`, :math:`zy`. See the :doc:`compute
pressure <compute_pressure>` command if you want the stress tensor
(pressure) of the entire system.
The stress tensor for atom :math:`I` is given by the following formula,
where :math:`a` and :math:`b` take on values :math:`x`, :math:`y`, :math:`z`
to generate the components of the tensor:
The stress tensor for atom :math:`I` is given by the following
formula, where :math:`a` and :math:`b` take on values :math:`x`,
:math:`y`, :math:`z` to generate the components of the tensor:
.. math::
S_{ab} = - m v_a v_b - W_{ab}
The first term is a kinetic energy contribution for atom :math:`I`. See
details below on how the specified *temp-ID* can affect the velocities
used in this calculation. The second term is the virial
contribution due to intra and intermolecular interactions,
where the exact computation details are determined by the compute style.
The first term is a kinetic energy contribution for atom :math:`I`.
See details below on how the specified *temp-ID* can affect the
velocities used in this calculation. The second term is the virial
contribution due to intra and intermolecular interactions, where the
exact computation details are determined by the compute style.
In case of compute *stress/atom*, the virial contribution is:
@ -68,29 +67,26 @@ In case of compute *stress/atom*, the virial contribution is:
& + \frac{1}{3} \sum_{n = 1}^{N_a} (r_{1_a} F_{1_b} + r_{2_a} F_{2_b} + r_{3_a} F_{3_b}) + \frac{1}{4} \sum_{n = 1}^{N_d} (r_{1_a} F_{1_b} + r_{2_a} F_{2_b} + r_{3_a} F_{3_b} + r_{4_a} F_{4_b}) \\
& + \frac{1}{4} \sum_{n = 1}^{N_i} (r_{1_a} F_{1_b} + r_{2_a} F_{2_b} + r_{3_a} F_{3_b} + r_{4_a} F_{4_b}) + {\rm Kspace}(r_{i_a},F_{i_b}) + \sum_{n = 1}^{N_f} r_{i_a} F_{i_b}
The first term is a pairwise energy
contribution where :math:`n` loops over the :math:`N_p`
neighbors of atom :math:`I`, :math:`\mathbf{r}_1` and :math:`\mathbf{r}_2`
are the positions of the 2 atoms in the pairwise interaction,
and :math:`\mathbf{F}_1` and :math:`\mathbf{F}_2` are the forces
on the 2 atoms resulting from the pairwise interaction.
The second term is a bond contribution of
similar form for the :math:`N_b` bonds which atom :math:`I` is part of.
There are similar terms for the :math:`N_a` angle, :math:`N_d` dihedral,
and :math:`N_i` improper interactions atom :math:`I` is part of.
There is also a term for the KSpace
contribution from long-range Coulombic interactions, if defined.
Finally, there is a term for the :math:`N_f` :doc:`fixes <fix>` that apply
internal constraint forces to atom :math:`I`. Currently, only the
:doc:`fix shake <fix_shake>` and :doc:`fix rigid <fix_rigid>` commands
contribute to this term.
As the coefficients in the formula imply, a virial contribution
produced by a small set of atoms (e.g. 4 atoms in a dihedral or 3
atoms in a Tersoff 3-body interaction) is assigned in equal portions
to each atom in the set. E.g. 1/4 of the dihedral virial to each of
the 4 atoms, or 1/3 of the fix virial due to SHAKE constraints applied
to atoms in a water molecule via the :doc:`fix shake <fix_shake>`
command.
The first term is a pairwise energy contribution where :math:`n` loops
over the :math:`N_p` neighbors of atom :math:`I`, :math:`\mathbf{r}_1`
and :math:`\mathbf{r}_2` are the positions of the 2 atoms in the
pairwise interaction, and :math:`\mathbf{F}_1` and
:math:`\mathbf{F}_2` are the forces on the 2 atoms resulting from the
pairwise interaction. The second term is a bond contribution of
similar form for the :math:`N_b` bonds which atom :math:`I` is part
of. There are similar terms for the :math:`N_a` angle, :math:`N_d`
dihedral, and :math:`N_i` improper interactions atom :math:`I` is part
of. There is also a term for the KSpace contribution from long-range
Coulombic interactions, if defined. Finally, there is a term for the
:math:`N_f` :doc:`fixes <fix>` that apply internal constraint forces
to atom :math:`I`. Currently, only the :doc:`fix shake <fix_shake>`
and :doc:`fix rigid <fix_rigid>` commands contribute to this term. As
the coefficients in the formula imply, a virial contribution produced
by a small set of atoms (e.g. 4 atoms in a dihedral or 3 atoms in a
Tersoff 3-body interaction) is assigned in equal portions to each atom
in the set. E.g. 1/4 of the dihedral virial to each of the 4 atoms,
or 1/3 of the fix virial due to SHAKE constraints applied to atoms in
a water molecule via the :doc:`fix shake <fix_shake>` command.
In case of compute *centroid/stress/atom*, the virial contribution is:
@ -99,71 +95,69 @@ In case of compute *centroid/stress/atom*, the virial contribution is:
W_{ab} & = \sum_{n = 1}^{N_p} r_{I0_a} F_{I_b} + \sum_{n = 1}^{N_b} r_{I0_a} F_{I_b} + \sum_{n = 1}^{N_a} r_{I0_a} F_{I_b} + \sum_{n = 1}^{N_d} r_{I0_a} F_{I_b} + \sum_{n = 1}^{N_i} r_{I0_a} F_{I_b} \\
& + {\rm Kspace}(r_{i_a},F_{i_b}) + \sum_{n = 1}^{N_f} r_{i_a} F_{i_b}
As with compute *stress/atom*, the first, second, third, fourth and fifth terms
are pairwise, bond, angle, dihedral and improper contributions,
but instead of assigning the virial contribution equally to each atom,
only the force :math:`\mathbf{F}_I` acting on atom :math:`I`
due to the interaction and the relative
position :math:`\mathbf{r}_{I0}` of the atom :math:`I` to the geometric center
of the interacting atoms, i.e. centroid, is used.
As the geometric center is different
for each interaction, the :math:`\mathbf{r}_{I0}` also differs.
The sixth and seventh terms, Kspace and :doc:`fix <fix>` contribution
respectively, are computed identical to compute *stress/atom*.
Although the total system virial is the same as compute *stress/atom*,
compute *centroid/stress/atom* is know to result in more consistent
heat flux values for angle, dihedrals and improper contributions
when computed via :doc:`compute heat/flux <compute_heat_flux>`.
As with compute *stress/atom*, the first, second, third, fourth and
fifth terms are pairwise, bond, angle, dihedral and improper
contributions, but instead of assigning the virial contribution
equally to each atom, only the force :math:`\mathbf{F}_I` acting on
atom :math:`I` due to the interaction and the relative position
:math:`\mathbf{r}_{I0}` of the atom :math:`I` to the geometric center
of the interacting atoms, i.e. centroid, is used. As the geometric
center is different for each interaction, the :math:`\mathbf{r}_{I0}`
also differs. The sixth and seventh terms, Kspace and :doc:`fix
<fix>` contribution respectively, are computed identical to compute
*stress/atom*. Although the total system virial is the same as
compute *stress/atom*, compute *centroid/stress/atom* is know to
result in more consistent heat flux values for angle, dihedrals and
improper contributions when computed via :doc:`compute heat/flux
<compute_heat_flux>`.
If no extra keywords are listed, the kinetic contribution
all of the virial contribution terms are
included in the per-atom stress tensor. If any extra keywords are
listed, only those terms are summed to compute the tensor. The
*virial* keyword means include all terms except the kinetic energy
*ke*\ .
If no extra keywords are listed, the kinetic contribution all of the
virial contribution terms are included in the per-atom stress tensor.
If any extra keywords are listed, only those terms are summed to
compute the tensor. The *virial* keyword means include all terms
except the kinetic energy *ke*\ .
Note that the stress for each atom is due to its interaction with all
other atoms in the simulation, not just with other atoms in the group.
Details of how compute *stress/atom* obtains the virial for individual atoms for
either pairwise or many-body potentials, and including the effects of
periodic boundary conditions is discussed in :ref:`(Thompson) <Thompson2>`.
The basic idea for many-body potentials is to treat each component of
the force computation between a small cluster of atoms in the same
manner as in the formula above for bond, angle, dihedral, etc
interactions. Namely the quantity :math:`\mathbf{r} \cdot \mathbf{F}`
is summed over the atoms in
the interaction, with the :math:`r` vectors unwrapped by periodic boundaries
so that the cluster of atoms is close together. The total
Details of how compute *stress/atom* obtains the virial for individual
atoms for either pairwise or many-body potentials, and including the
effects of periodic boundary conditions is discussed in
:ref:`(Thompson) <Thompson2>`. The basic idea for many-body
potentials is to treat each component of the force computation between
a small cluster of atoms in the same manner as in the formula above
for bond, angle, dihedral, etc interactions. Namely the quantity
:math:`\mathbf{r} \cdot \mathbf{F}` is summed over the atoms in the
interaction, with the :math:`r` vectors unwrapped by periodic
boundaries so that the cluster of atoms is close together. The total
contribution for the cluster interaction is divided evenly among those
atoms. Details of how compute *centroid/stress/atom* obtains
the virial for individual atoms
is given in :ref:`(Surblys) <Surblys1>`,
where the idea is that the virial of the atom :math:`I`
is the result of only the force :math:`\mathbf{F}_I` on the atom due
to the interaction
and its positional vector :math:`\mathbf{r}_{I0}`,
relative to the geometric center of the
interacting atoms, regardless of the number of participating atoms.
The periodic boundary treatment is identical to
atoms.
Details of how compute *centroid/stress/atom* obtains the virial for
individual atoms is given in :ref:`(Surblys) <Surblys1>`, where the
idea is that the virial of the atom :math:`I` is the result of only
the force :math:`\mathbf{F}_I` on the atom due to the interaction and
its positional vector :math:`\mathbf{r}_{I0}`, relative to the
geometric center of the interacting atoms, regardless of the number of
participating atoms. The periodic boundary treatment is identical to
that of compute *stress/atom*, and both of them reduce to identical
expressions for two-body interactions,
i.e. computed values for contributions from bonds and two-body pair styles,
such as :doc:`Lennard-Jones <pair_lj>`, will be the same,
while contributions from angles, dihedrals and impropers will be different.
expressions for two-body interactions, i.e. computed values for
contributions from bonds and two-body pair styles, such as
:doc:`Lennard-Jones <pair_lj>`, will be the same, while contributions
from angles, dihedrals and impropers will be different.
The :doc:`dihedral_style charmm <dihedral_charmm>` style calculates
pairwise interactions between 1-4 atoms. The virial contribution of
these terms is included in the pair virial, not the dihedral virial.
The KSpace contribution is calculated using the method in
:ref:`(Heyes) <Heyes2>` for the Ewald method and by the methodology described
in :ref:`(Sirk) <Sirk1>` for PPPM. The choice of KSpace solver is specified
by the :doc:`kspace_style pppm <kspace_style>` command. Note that for
PPPM, the calculation requires 6 extra FFTs each timestep that
per-atom stress is calculated. Thus it can significantly increase the
cost of the PPPM calculation if it is needed on a large fraction of
the simulation timesteps.
:ref:`(Heyes) <Heyes2>` for the Ewald method and by the methodology
described in :ref:`(Sirk) <Sirk1>` for PPPM. The choice of KSpace
solver is specified by the :doc:`kspace_style pppm <kspace_style>`
command. Note that for PPPM, the calculation requires 6 extra FFTs
each timestep that per-atom stress is calculated. Thus it can
significantly increase the cost of the PPPM calculation if it is
needed on a large fraction of the simulation timesteps.
The *temp-ID* argument can be used to affect the per-atom velocities
used in the kinetic energy contribution to the total stress. If the
@ -189,10 +183,10 @@ See the :doc:`compute voronoi/atom <compute_voronoi_atom>` command for
one possible way to estimate a per-atom volume.
Thus, if the diagonal components of the per-atom stress tensor are
summed for all atoms in the system and the sum is divided by :math:`dV`, where
:math:`d` = dimension and :math:`V` is the volume of the system,
the result should be :math:`-P`, where :math:`P`
is the total pressure of the system.
summed for all atoms in the system and the sum is divided by
:math:`dV`, where :math:`d` = dimension and :math:`V` is the volume of
the system, the result should be :math:`-P`, where :math:`P` is the
total pressure of the system.
These lines in an input script for a 3d system should yield that
result. I.e. the last 2 columns of thermo output will be the same:
@ -207,33 +201,43 @@ result. I.e. the last 2 columns of thermo output will be the same:
.. note::
The per-atom stress does not include any Lennard-Jones tail
corrections to the pressure added by the :doc:`pair_modify tail yes <pair_modify>` command, since those are contributions to the
global system pressure.
corrections to the pressure added by the :doc:`pair_modify tail yes
<pair_modify>` command, since those are contributions to the global
system pressure.
Output info
"""""""""""
This compute *stress/atom* calculates a per-atom array with 6 columns, which can be
accessed by indices 1-6 by any command that uses per-atom values from
a compute as input.
The compute *centroid/stress/atom* produces a per-atom array with 9 columns,
but otherwise can be used in an identical manner to compute *stress/atom*.
See the :doc:`Howto output <Howto_output>` doc page
for an overview of LAMMPS output options.
Compute *stress/atom* calculates a per-atom array with 6 columns,
which can be accessed by indices 1-6 by any command that uses per-atom
values from a compute as input. Compute *centroid/stress/atom*
produces a per-atom array with 9 columns, but otherwise can be used in
an identical manner to compute *stress/atom*. See the :doc:`Howto
output <Howto_output>` doc page for an overview of LAMMPS output
options.
The per-atom array values will be in pressure\*volume
:doc:`units <units>` as discussed above.
The per-atom array values will be in pressure\*volume :doc:`units
<units>` as discussed above.
Restrictions
""""""""""""
Currently (Spring 2020), compute *centroid/stress/atom* does not support
pair styles with many-body interactions, such as :doc:`Tersoff
<pair_tersoff>`, or pair styles with long-range Coulomb interactions.
LAMMPS will generate an error in such cases. In principal, equivalent
formulation to that of angle, dihedral and improper contributions in the
virial :math:`W_{ab}` formula can also be applied to the many-body pair
styles, and is planned in the future.
Currently, compute *centroid/stress/atom* does not support pair styles
with many-body interactions (:doc:`EAM <pair_eam>` is an exception,
since its computations are performed pairwise), nor granular pair
styles with pairwise forces which are not aligned with the vector
between the pair of particles. All bond styles are supported. All
angle, dihedral, improper styles are supported with the exception of
USER-INTEL and KOKKOS variants of specific styles. It also does not
support models with long-range Coulombic or dispersion forces,
i.e. the kspace_style command in LAMMPS. It also does not support the
following fixes which add rigid-body constraints: :doc:`fix shake
<fix_shake>`, :doc:`fix rattle <fix_rattle>`, :doc:`fix rigid
<fix_rigid>`, :doc:`fix rigid/small <fix_rigid>`.
LAMMPS will generate an error if one of these options is included in
your model. Extension of centroid stress calculations to these force
and fix styles is planned for the futre.
Related commands
""""""""""""""""