Merge pull request #3797 from Compizfox/develop
Include bond forces in `compute stress/cartesian`
This commit is contained in:
Axel Kohlmeyer
@ -153,11 +153,11 @@ KOKKOS, o = OPENMP, t = OPT.
|
||||
* :doc:`sph/t/atom <compute_sph_t_atom>`
|
||||
* :doc:`spin <compute_spin>`
|
||||
* :doc:`stress/atom <compute_stress_atom>`
|
||||
* :doc:`stress/cartesian <compute_stress_profile>`
|
||||
* :doc:`stress/cylinder <compute_stress_profile>`
|
||||
* :doc:`stress/cartesian <compute_stress_cartesian>`
|
||||
* :doc:`stress/cylinder <compute_stress_curvilinear>`
|
||||
* :doc:`stress/mop <compute_stress_mop>`
|
||||
* :doc:`stress/mop/profile <compute_stress_mop>`
|
||||
* :doc:`stress/spherical <compute_stress_profile>`
|
||||
* :doc:`stress/spherical <compute_stress_curvilinear>`
|
||||
* :doc:`stress/tally <compute_tally>`
|
||||
* :doc:`tdpd/cc/atom <compute_tdpd_cc_atom>`
|
||||
* :doc:`temp (k) <compute_temp>`
|
||||
|
||||
@ -307,11 +307,11 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
|
||||
* :doc:`sph/t/atom <compute_sph_t_atom>` - per-atom internal temperature of Smooth-Particle Hydrodynamics atoms
|
||||
* :doc:`spin <compute_spin>` - magnetic quantities for a system of atoms having spins
|
||||
* :doc:`stress/atom <compute_stress_atom>` - stress tensor for each atom
|
||||
* :doc:`stress/cartesian <compute_stress_profile>` - stress tensor in cartesian coordinates
|
||||
* :doc:`stress/cylinder <compute_stress_profile>` - stress tensor in cylindrical coordinates
|
||||
* :doc:`stress/cartesian <compute_stress_cartesian>` - stress tensor in cartesian coordinates
|
||||
* :doc:`stress/cylinder <compute_stress_curvilinear>` - stress tensor in cylindrical coordinates
|
||||
* :doc:`stress/mop <compute_stress_mop>` - normal components of the local stress tensor using the method of planes
|
||||
* :doc:`stress/mop/profile <compute_stress_mop>` - profile of the normal components of the local stress tensor using the method of planes
|
||||
* :doc:`stress/spherical <compute_stress_profile>` - stress tensor in spherical coordinates
|
||||
* :doc:`stress/spherical <compute_stress_curvilinear>` - stress tensor in spherical coordinates
|
||||
* :doc:`stress/tally <compute_tally>` - stress between two groups of atoms via the tally callback mechanism
|
||||
* :doc:`tdpd/cc/atom <compute_tdpd_cc_atom>` - per-atom chemical concentration of a specified species for each tDPD particle
|
||||
* :doc:`temp <compute_temp>` - temperature of group of atoms
|
||||
|
||||
122
doc/src/compute_stress_cartesian.rst
Normal file
122
doc/src/compute_stress_cartesian.rst
Normal file
@ -0,0 +1,122 @@
|
||||
.. index:: compute stress/cartesian
|
||||
|
||||
compute stress/cartesian command
|
||||
==================================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
compute ID group-ID stress/cartesian args
|
||||
|
||||
* ID, group-ID are documented in :doc:`compute <compute>` command
|
||||
* args = argument specific to the compute style
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
*stress/cartesian* args = dim1 bin_width1 dim2 bin_width2 keyword
|
||||
dim1 = *x* or *y* or *z*
|
||||
bin_width1 = width of the bin
|
||||
dim2 = *x* or *y* or *z* or *NULL*
|
||||
bin_width2 = width of the bin
|
||||
keyword = *ke* or *pair* or *bond*
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
compute 1 all stress/cartesian x 0.1 NULL 0
|
||||
compute 1 all stress/cartesian y 0.1 z 0.1
|
||||
compute 1 all stress/cartesian x 0.1 NULL 0 ke pair
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
Compute *stress/cartesian* defines computations that calculate profiles of the
|
||||
diagonal components of the local stress tensor over one or two Cartesian
|
||||
dimensions, as described in :ref:`(Ikeshoji)<Ikeshoji2>`. The stress tensor is
|
||||
split into a kinetic contribution :math:`P^k` and a virial contribution
|
||||
:math:`P^v`. The sum gives the total stress tensor :math:`P = P^k+P^v`.
|
||||
This compute obeys momentum balance through fluid interfaces. They use the
|
||||
Irving--Kirkwood contour, which is the straight line between particle pairs.
|
||||
|
||||
.. versionadded:: TBD
|
||||
|
||||
Added support for bond styles
|
||||
|
||||
This compute only supports pair and bond (no angle, dihedral, improper,
|
||||
or kspace) forces. By default, if no extra keywords are specified, all
|
||||
supported contributions to the stress are included (ke, pair, bond). If any
|
||||
keywords are specified, then only those components are summed.
|
||||
|
||||
Output info
|
||||
"""""""""""
|
||||
|
||||
The output columns for *stress/cartesian* are the position of the
|
||||
center of the local volume in the first and second dimensions, number
|
||||
density, :math:`P^k_{xx}`, :math:`P^k_{yy}`, :math:`P^k_{zz}`,
|
||||
:math:`P^v_{xx}`, :math:`P^v_{yy}`, and :math:`P^v_{zz}`. There are 8
|
||||
columns when one dimension is specified and 9 columns when two
|
||||
dimensions are specified. The number of bins (rows) is
|
||||
:math:`(L_1/b_1)(L_2/b_2)`, where :math:`L_1` and :math:`L_2` are the lengths
|
||||
of the simulation box in the specified dimensions and :math:`b_1` and
|
||||
:math:`b_2` are the specified bin widths. When only one dimension is
|
||||
specified, the number of bins (rows) is :math:`L_1/b_1`.
|
||||
|
||||
This array can be output with :doc:`fix ave/time <fix_ave_time>`,
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
compute p all stress/cartesian x 0.1
|
||||
fix 2 all ave/time 100 1 100 c_p[*] file dump_p.out mode vector
|
||||
|
||||
The values calculated by this compute are "intensive". The stress
|
||||
values will be in pressure :doc:`units <units>`. The number density
|
||||
values are in inverse volume :doc:`units <units>`.
|
||||
|
||||
NOTE 1: The local stress does not include any Lennard-Jones tail
|
||||
corrections to the stress added by the
|
||||
:doc:`pair_modify tail yes <pair_modify>`
|
||||
command, since those are contributions to the global system pressure.
|
||||
|
||||
NOTE 2: The local stress profiles generated by these computes are
|
||||
similar to those obtained by the
|
||||
:doc:`method-of-planes (MOP) <compute_stress_mop>`.
|
||||
A key difference is that compute
|
||||
:doc:`stress/mop/profile <compute_stress_mop>`
|
||||
considers particles crossing a set of planes, while
|
||||
*stress/cartesian* computes averages for a set of small volumes.
|
||||
Moreover, this compute computes the diagonal components of the stress
|
||||
tensor :math:`P_{xx}`, :math:`P_{yy}`, and :math:`P_{zz}`, while
|
||||
*stress/mop/profile* computes the components
|
||||
:math:`P_{ix}`, :math:`P_{iy}`, and :math:`P_{iz}`, where :math:`i` is the
|
||||
direction normal to the plane.
|
||||
|
||||
More information on the similarities and differences can be found in
|
||||
:ref:`(Ikeshoji)<Ikeshoji2>`.
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
These computes calculate the stress tensor contributions for pair and bond
|
||||
forces only (no angle, dihedral, improper, or kspace force).
|
||||
It requires pairwise force calculations not available for most
|
||||
many-body pair styles.
|
||||
|
||||
These computes are part of the EXTRA-COMPUTE package. They are only
|
||||
enabled if LAMMPS was built with that package. See the :doc:`Build
|
||||
package <Build_package>` doc page for more info.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`compute stress/atom <compute_stress_atom>`, :doc:`compute pressure <compute_pressure>`,
|
||||
:doc:`compute stress/mop/profile <compute_stress_mop>`, :doc:`compute stress/spherical <compute_stress_curvilinear>`,
|
||||
:doc:`compute stress/cylinder <compute_stress_curvilinear>`
|
||||
|
||||
----------
|
||||
|
||||
.. _Ikeshoji2:
|
||||
|
||||
**(Ikeshoji)** Ikeshoji, Hafskjold, Furuholt, Mol Sim, 29, 101-109, (2003).
|
||||
@ -1,11 +1,6 @@
|
||||
.. index:: compute stress/cartesian
|
||||
.. index:: compute stress/cylinder
|
||||
.. index:: compute stress/spherical
|
||||
|
||||
|
||||
compute stress/cartesian command
|
||||
==================================
|
||||
|
||||
compute stress/cylinder command
|
||||
=================================
|
||||
|
||||
@ -20,15 +15,11 @@ Syntax
|
||||
compute ID group-ID style args
|
||||
|
||||
* ID, group-ID are documented in :doc:`compute <compute>` command
|
||||
* style = stress/cartesian or stress/spherical or stress/cylinder
|
||||
* style = stress/spherical or stress/cylinder
|
||||
* args = argument specific to the compute style
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
*stress/cartesian* args = dim bin_width
|
||||
dim = *x* or *y* or *z*
|
||||
bin_width = width of the bin
|
||||
one or two dim/bin_width pairs may be appended
|
||||
*stress/cylinder* args = zlo zh Rmax bin_width keyword
|
||||
zlo = minimum z-boundary for cylinder
|
||||
zhi = maximum z-boundary for cylinder
|
||||
@ -45,9 +36,6 @@ Examples
|
||||
""""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
compute 1 all stress/cartesian x 0.1
|
||||
compute 1 all stress/cartesian y 0.25 z 0.1
|
||||
compute 1 all stress/cylinder -10.0 10.0 15.0 0.25
|
||||
compute 1 all stress/cylinder -10.0 10.0 15.0 0.25 ke no
|
||||
compute 1 all stress/spherical 0 0 0 0.1 10
|
||||
@ -55,21 +43,19 @@ Examples
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
Compute *stress/cartesian*, compute *stress/cylinder*, and compute
|
||||
Compute *stress/cylinder*, and compute
|
||||
*stress/spherical* define computations that calculate profiles of the
|
||||
diagonal components of the local stress tensor in the specified
|
||||
coordinate system. The stress tensor is split into a kinetic
|
||||
contribution :math:`P^k` and a virial contribution :math:`P^v`. The sum
|
||||
gives the total stress tensor :math:`P = P^k+P^v`. These computes can
|
||||
for example be used to calculate the diagonal components of the local
|
||||
stress tensor of interfaces with flat, cylindrical, or spherical
|
||||
stress tensor of surfaces with cylindrical or spherical
|
||||
symmetry. These computes obeys momentum balance through fluid
|
||||
interfaces. They use the Irving--Kirkwood contour, which is the straight
|
||||
line between particle pairs.
|
||||
|
||||
The *stress/cartesian* computes the stress profile along one or two
|
||||
Cartesian coordinates, as described in :ref:`(Ikeshoji)<Ikeshoji2>`. The
|
||||
compute *stress/cylinder* computes the stress profile along the
|
||||
The compute *stress/cylinder* computes the stress profile along the
|
||||
radial direction in cylindrical coordinates, as described in
|
||||
:ref:`(Addington)<Addington1>`. The compute *stress/spherical*
|
||||
computes the stress profile along the radial direction in spherical
|
||||
@ -79,17 +65,6 @@ coordinates, as described in :ref:`(Ikeshoji)<Ikeshoji2>`.
|
||||
Output info
|
||||
"""""""""""
|
||||
|
||||
The output columns for *stress/cartesian* are the position of the
|
||||
center of the local volume in the first and second dimensions, number
|
||||
density, :math:`P^k_{xx}`, :math:`P^k_{yy}`, :math:`P^k_{zz}`,
|
||||
:math:`P^v_{xx}`, :math:`P^v_{yy}`, and :math:`P^v_{zz}`. There are 8
|
||||
columns when one dimension is specified and 9 columns when two
|
||||
dimensions are specified. The number of bins (rows) is
|
||||
:math:`(L_1/b_1)(L_2/b_2)`, where :math:`L_1` and :math:`L_2` are the lengths
|
||||
of the simulation box in the specified dimensions and :math:`b_1` and
|
||||
:math:`b_2` are the specified bin widths. When only one dimension is
|
||||
specified, the number of bins (rows) is :math:`L_1/b_1`.
|
||||
|
||||
The default output columns for *stress/cylinder* are the radius to the
|
||||
center of the cylindrical shell, number density, :math:`P^k_{rr}`,
|
||||
:math:`P^k_{\phi\phi}`, :math:`P^k_{zz}`, :math:`P^v_{rr}`,
|
||||
@ -111,7 +86,7 @@ This array can be output with :doc:`fix ave/time <fix_ave_time>`,
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
compute p all stress/cartesian x 0.1
|
||||
compute 1 all stress/spherical 0 0 0 0.1 10
|
||||
fix 2 all ave/time 100 1 100 c_p[*] file dump_p.out mode vector
|
||||
|
||||
The values calculated by this compute are "intensive". The stress
|
||||
@ -123,16 +98,6 @@ corrections to the stress added by the
|
||||
:doc:`pair_modify tail yes <pair_modify>`
|
||||
command, since those are contributions to the global system pressure.
|
||||
|
||||
NOTE 2: The local stress profiles generated by these computes are
|
||||
similar to those obtained by the
|
||||
:doc:`method-of-planes (MOP) <compute_stress_mop>`.
|
||||
A key difference
|
||||
is that compute `stress/mop/profile <compute_stress_mop>`
|
||||
considers particles crossing a set of planes, while
|
||||
*stress/cartesian* computes averages for a set of small volumes.
|
||||
More information on the similarities and differences can be found in
|
||||
:ref:`(Ikeshoji)<Ikeshoji2>`.
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
@ -150,7 +115,8 @@ package <Build_package>` doc page for more info.
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`compute stress/atom <compute_stress_atom>`, :doc:`compute pressure <compute_pressure>`, :doc:`compute stress/mop/profile <compute_stress_mop>`
|
||||
:doc:`compute stress/atom <compute_stress_atom>`, :doc:`compute pressure <compute_pressure>`,
|
||||
:doc:`compute stress/mop/profile <compute_stress_mop>`, :doc:`compute stress/cartesian <compute_stress_cartesian>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
@ -88,11 +88,16 @@ system pressure.
|
||||
|
||||
NOTE 3: The local stress profile generated by compute
|
||||
*stress/mop/profile* is similar to that obtained by compute
|
||||
:doc:`stress/cartesian <compute_stress_profile>`. A key difference is
|
||||
that compute *stress/mop/profile* considers particles crossing a set of
|
||||
planes, while compute *stress/cartesian* computes averages for a set of
|
||||
small volumes. More information on the similarities and differences can
|
||||
be found in :ref:`(Ikeshoji)<Ikeshoji2>`.
|
||||
:doc:`stress/cartesian <compute_stress_cartesian>`.
|
||||
A key difference is that compute *stress/mop/profile*
|
||||
considers particles crossing a set of planes, while
|
||||
*stress/cartesian* computes averages for a set
|
||||
of small volumes.
|
||||
Moreover, *stress/cartesian* compute computes the diagonal components of the stress
|
||||
tensor :math:`P_{xx}`, :math:`P_{yy}`, and :math:`P_{zz}`, while
|
||||
*stress/mop/profile* computes the components
|
||||
:math:`P_{ix}`, :math:`P_{iy}`, and :math:`P_{iz}`, where :math:`i` is the
|
||||
direction normal to the plane.
|
||||
|
||||
Output info
|
||||
"""""""""""
|
||||
@ -143,9 +148,9 @@ Related commands
|
||||
|
||||
:doc:`compute stress/atom <compute_stress_atom>`,
|
||||
:doc:`compute pressure <compute_pressure>`,
|
||||
:doc:`compute stress/cartesian <compute_stress_profile>`,
|
||||
:doc:`compute stress/cylinder <compute_stress_profile>`,
|
||||
:doc:`compute stress/spherical <compute_stress_profile>`
|
||||
:doc:`compute stress/cartesian <compute_stress_cartesian>`,
|
||||
:doc:`compute stress/cylinder <compute_stress_curvilinear>`,
|
||||
:doc:`compute stress/spherical <compute_stress_curvilinear>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
|
||||
@ -14,6 +14,7 @@
|
||||
#include "compute_stress_cartesian.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "bond.h"
|
||||
#include "citeme.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
@ -61,14 +62,12 @@ ComputeStressCartesian::ComputeStressCartesian(LAMMPS *lmp, int narg, char **arg
|
||||
{
|
||||
if (lmp->citeme) lmp->citeme->add(cite_compute_stress_cartesian);
|
||||
|
||||
// narg == 5 for one-dimensional and narg == 7 for two-dimensional
|
||||
if (narg == 5)
|
||||
dims = 1;
|
||||
else if (narg == 7)
|
||||
dims = 2;
|
||||
else
|
||||
error->all(FLERR, "Illegal compute stress/cartesian command. Illegal number of arguments.");
|
||||
if (narg < 7) error->all(FLERR, "Illegal compute stress/cartesian command: illegal number of arguments.");
|
||||
|
||||
// no triclinic boxes
|
||||
if (domain->triclinic) error->all(FLERR, "Compute stress/cartesian requires an orthogonal box");
|
||||
|
||||
// Direction of first dimension
|
||||
if (strcmp(arg[3], "x") == 0)
|
||||
dir1 = 0;
|
||||
else if (strcmp(arg[3], "y") == 0)
|
||||
@ -78,14 +77,8 @@ ComputeStressCartesian::ComputeStressCartesian(LAMMPS *lmp, int narg, char **arg
|
||||
else
|
||||
error->all(FLERR, "Illegal compute stress/cartesian direction: {}", arg[3]);
|
||||
|
||||
dir2 = 0;
|
||||
bin_width1 = utils::numeric(FLERR, arg[4], false, lmp);
|
||||
bin_width2 = domain->boxhi[dir2] - domain->boxlo[dir2];
|
||||
nbins1 = (int) ((domain->boxhi[dir1] - domain->boxlo[dir1]) / bin_width1);
|
||||
nbins2 = 1;
|
||||
|
||||
// no triclinic boxes
|
||||
if (domain->triclinic) error->all(FLERR, "Compute stress/cartesian requires an orthogonal box");
|
||||
|
||||
// adjust bin width if not a perfect match
|
||||
double tmp_binwidth = (domain->boxhi[dir1] - domain->boxlo[dir1]) / nbins1;
|
||||
@ -93,14 +86,23 @@ ComputeStressCartesian::ComputeStressCartesian(LAMMPS *lmp, int narg, char **arg
|
||||
utils::logmesg(lmp, "Adjusting first bin width for compute {} from {:.6f} to {:.6f}\n", style,
|
||||
bin_width1, tmp_binwidth);
|
||||
bin_width1 = tmp_binwidth;
|
||||
invV = bin_width1;
|
||||
|
||||
if (bin_width1 <= 0.0)
|
||||
error->all(FLERR, "Illegal compute stress/cartesian command. First bin width must be > 0");
|
||||
else if (bin_width1 > domain->boxhi[dir1] - domain->boxlo[dir1])
|
||||
error->all(FLERR, "Illegal compute stress/cartesian command. First bin width > box.");
|
||||
|
||||
invV = bin_width1;
|
||||
if (dims == 2) {
|
||||
// Direction of second dimension
|
||||
if (strcmp(arg[5], "NULL") == 0) {
|
||||
dims = 1;
|
||||
dir2 = 0;
|
||||
bin_width2 = domain->boxhi[dir2] - domain->boxlo[dir2];
|
||||
nbins2 = 1;
|
||||
} else {
|
||||
dims = 2;
|
||||
|
||||
// Direction of first dimension
|
||||
if (strcmp(arg[5], "x") == 0)
|
||||
dir2 = 0;
|
||||
else if (strcmp(arg[5], "y") == 0)
|
||||
@ -149,6 +151,21 @@ ComputeStressCartesian::ComputeStressCartesian(LAMMPS *lmp, int narg, char **arg
|
||||
if (box_incompatible)
|
||||
error->all(FLERR, "Must not use compute stress/cartesian on variable box dimension");
|
||||
|
||||
// process optional args
|
||||
if (narg > 7) {
|
||||
compute_ke = false;
|
||||
compute_pair = false;
|
||||
compute_bond = false;
|
||||
int iarg = 7;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg], "ke") == 0) compute_ke = true;
|
||||
else if (strcmp(arg[iarg], "pair") == 0) compute_pair = true;
|
||||
else if (strcmp(arg[iarg], "bond") == 0) compute_bond = true;
|
||||
else error->all(FLERR, "Illegal compute stress/atom command");
|
||||
iarg++;
|
||||
}
|
||||
}
|
||||
|
||||
for (int i = 0; i < 3; i++)
|
||||
if ((dims == 1 && i != dir1) || (dims == 2 && (i != dir1 && i != dir2)))
|
||||
invV *= domain->boxhi[i] - domain->boxlo[i];
|
||||
@ -229,13 +246,6 @@ void ComputeStressCartesian::init_list(int /* id */, NeighList *ptr)
|
||||
|
||||
void ComputeStressCartesian::compute_array()
|
||||
{
|
||||
int i, j, ii, jj, inum, jnum, itype, jtype;
|
||||
int bin, bin1, bin2;
|
||||
tagint itag, jtag;
|
||||
double xtmp, ytmp, ztmp, delx, dely, delz;
|
||||
double rsq, fpair, factor_coul, factor_lj;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
|
||||
double **x = atom->x;
|
||||
double **v = atom->v;
|
||||
double *mass = atom->mass;
|
||||
@ -251,13 +261,13 @@ void ComputeStressCartesian::compute_array()
|
||||
// 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;
|
||||
int inum = list->inum;
|
||||
int *ilist = list->ilist;
|
||||
int *numneigh = list->numneigh;
|
||||
int **firstneigh = list->firstneigh;
|
||||
|
||||
// Zero arrays
|
||||
for (bin = 0; bin < nbins1 * nbins2; bin++) {
|
||||
for (int bin = 0; bin < nbins1 * nbins2; bin++) {
|
||||
tdens[bin] = 0;
|
||||
tpkxx[bin] = 0;
|
||||
tpkyy[bin] = 0;
|
||||
@ -268,105 +278,129 @@ void ComputeStressCartesian::compute_array()
|
||||
}
|
||||
|
||||
// calculate number density and kinetic contribution to pressure
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
bin1 = (int) ((x[i][dir1] - boxlo[dir1]) / bin_width1) % nbins1;
|
||||
bin2 = 0;
|
||||
if (dims == 2) bin2 = (int) ((x[i][dir2] - boxlo[dir2]) / bin_width2) % nbins2;
|
||||
if (compute_ke) {
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
int bin1 = (int) ((x[i][dir1] - boxlo[dir1]) / bin_width1) % nbins1;
|
||||
int bin2 = 0;
|
||||
if (dims == 2) bin2 = (int) ((x[i][dir2] - boxlo[dir2]) / bin_width2) % nbins2;
|
||||
|
||||
// Apply periodic boundary conditions and avoid out of range access
|
||||
if (domain->periodicity[dir1] == 1) {
|
||||
if (bin1 < 0)
|
||||
bin1 = (bin1 + nbins1) % nbins1;
|
||||
// Apply periodic boundary conditions and avoid out of range access
|
||||
if (domain->periodicity[dir1] == 1) {
|
||||
if (bin1 < 0)
|
||||
bin1 = (bin1 + nbins1) % nbins1;
|
||||
else if (bin1 >= nbins1)
|
||||
bin1 = (bin1 - nbins1) % nbins1;
|
||||
} else if (bin1 < 0)
|
||||
bin1 = 0;
|
||||
else if (bin1 >= nbins1)
|
||||
bin1 = (bin1 - nbins1) % nbins1;
|
||||
} else if (bin1 < 0)
|
||||
bin1 = 0;
|
||||
else if (bin1 >= nbins1)
|
||||
bin1 = nbins1 - 1;
|
||||
bin1 = nbins1 - 1;
|
||||
|
||||
if (domain->periodicity[dir2] == 1) {
|
||||
if (bin2 < 0)
|
||||
bin2 = (bin2 + nbins2) % nbins2;
|
||||
if (domain->periodicity[dir2] == 1) {
|
||||
if (bin2 < 0)
|
||||
bin2 = (bin2 + nbins2) % nbins2;
|
||||
else if (bin2 >= nbins2)
|
||||
bin2 = (bin2 - nbins2) % nbins2;
|
||||
} else if (bin2 < 0)
|
||||
bin2 = 0;
|
||||
else if (bin2 >= nbins2)
|
||||
bin2 = (bin2 - nbins2) % nbins2;
|
||||
} else if (bin2 < 0)
|
||||
bin2 = 0;
|
||||
else if (bin2 >= nbins2)
|
||||
bin2 = nbins2 - 1;
|
||||
bin2 = nbins2 - 1;
|
||||
|
||||
j = bin1 + bin2 * nbins1;
|
||||
tdens[j] += 1;
|
||||
tpkxx[j] += mass[type[i]] * v[i][0] * v[i][0];
|
||||
tpkyy[j] += mass[type[i]] * v[i][1] * v[i][1];
|
||||
tpkzz[j] += mass[type[i]] * v[i][2] * v[i][2];
|
||||
int j = bin1 + bin2 * nbins1;
|
||||
tdens[j] += 1;
|
||||
tpkxx[j] += mass[type[i]] * v[i][0] * v[i][0];
|
||||
tpkyy[j] += mass[type[i]] * v[i][1] * v[i][1];
|
||||
tpkzz[j] += mass[type[i]] * v[i][2] * v[i][2];
|
||||
}
|
||||
}
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
// skip if I or J are not in group
|
||||
// for newton = 0 and J = ghost atom,
|
||||
// need to ensure I,J pair is only output by one proc
|
||||
// use same itag,jtag logic as in Neighbor::neigh_half_nsq()
|
||||
// for flag = 0, just count pair interactions within force cutoff
|
||||
// for flag = 1, calculate requested output fields
|
||||
if (compute_pair && force->pair) {
|
||||
for (int ii = 0; ii < inum; ii++) {
|
||||
int i = ilist[ii];
|
||||
|
||||
Pair *pair = force->pair;
|
||||
double **cutsq = force->pair->cutsq;
|
||||
// skip if I or J are not in group
|
||||
if (!(mask[i] & groupbit)) continue;
|
||||
|
||||
double xi1, xi2;
|
||||
double xi1 = x[i][dir1] - boxlo[dir1];
|
||||
double xi2 = x[i][dir2] - boxlo[dir2];
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
if (!(mask[i] & groupbit)) continue;
|
||||
for (int jj = 0; jj < numneigh[i]; jj++) {
|
||||
int j = firstneigh[i][jj];
|
||||
double factor_lj = special_lj[sbmask(j)];
|
||||
double factor_coul = special_coul[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
if (!(mask[j] & groupbit)) continue;
|
||||
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
xi1 = x[i][dir1] - boxlo[dir1];
|
||||
xi2 = x[i][dir2] - boxlo[dir2];
|
||||
itag = tag[i];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_lj = special_lj[sbmask(j)];
|
||||
factor_coul = special_coul[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
if (!(mask[j] & groupbit)) continue;
|
||||
|
||||
// itag = jtag is possible for long cutoffs that include images of self
|
||||
// do calculation only on appropriate processor
|
||||
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;
|
||||
// for newton = 0 and J = ghost atom, need to ensure I,J pair is only output by one proc
|
||||
// use same tag[i],tag[j] logic as in Neighbor::neigh_half_nsq()
|
||||
if (newton_pair == 0 && j >= nlocal) {
|
||||
if (tag[i] > tag[j]) {
|
||||
if ((tag[i] + tag[j]) % 2 == 0) continue;
|
||||
} else if (tag[i] < tag[j]) {
|
||||
if ((tag[i] + tag[j]) % 2 == 1) continue;
|
||||
} else {
|
||||
// tag[i] = tag[j] is possible for long cutoffs that include images of self
|
||||
if (x[j][2] < x[i][2]) continue;
|
||||
if (x[j][2] == x[i][2]) {
|
||||
if (x[j][1] < x[i][1]) continue;
|
||||
if (x[j][1] == x[i][1] && x[j][0] < x[i][0]) continue;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double delx = x[j][0] - x[i][0];
|
||||
double dely = x[j][1] - x[i][1];
|
||||
double delz = x[j][2] - x[i][2];
|
||||
double rsq = delx * delx + dely * dely + delz * delz;
|
||||
|
||||
// Check if inside cut-off
|
||||
int itype = type[i];
|
||||
int jtype = type[j];
|
||||
if (rsq >= force->pair->cutsq[itype][jtype]) continue;
|
||||
|
||||
double fpair;
|
||||
force->pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair);
|
||||
compute_pressure(fpair, xi1, xi2, delx, dely, delz);
|
||||
}
|
||||
delx = x[j][0] - xtmp;
|
||||
dely = x[j][1] - ytmp;
|
||||
delz = x[j][2] - ztmp;
|
||||
}
|
||||
}
|
||||
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
jtype = type[j];
|
||||
// Loop over all bonds
|
||||
if (compute_bond && force->bond) {
|
||||
for (int i_bond = 0; i_bond < neighbor->nbondlist; i_bond++) {
|
||||
// i == atom1, j == atom2
|
||||
int i = neighbor->bondlist[i_bond][0];
|
||||
int j = neighbor->bondlist[i_bond][1];
|
||||
int btype = neighbor->bondlist[i_bond][2];
|
||||
|
||||
// Check if inside cut-off
|
||||
if (rsq >= cutsq[itype][jtype]) continue;
|
||||
pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair);
|
||||
compute_pressure(fpair, xi1, xi2, delx, dely, delz);
|
||||
// Skip if one of both atoms is not in group
|
||||
if (!(mask[i] & groupbit)) continue;
|
||||
if (!(mask[j] & groupbit)) continue;
|
||||
|
||||
// if newton_bond is off and atom2 is a ghost atom, only compute this on one processor
|
||||
if (!force->newton_bond && j >= nlocal) {
|
||||
if (tag[i] > tag[j]) {
|
||||
if ((tag[i] + tag[j]) % 2 == 0) continue;
|
||||
} else if (tag[i] < tag[j]) {
|
||||
if ((tag[i] < tag[j]) % 2 == 1) continue;
|
||||
}
|
||||
}
|
||||
|
||||
double dx = x[j][0] - x[i][0];
|
||||
double dy = x[j][1] - x[i][1];
|
||||
double dz = x[j][2] - x[i][2];
|
||||
double rsq = dx*dx + dy*dy + dz*dz;
|
||||
double xi = x[i][dir1] - boxlo[dir1];
|
||||
double yi = x[i][dir2] - boxlo[dir2];
|
||||
|
||||
double fbond;
|
||||
force->bond->single(btype, rsq, i, j, fbond);
|
||||
compute_pressure(fbond, xi, yi, dx, dy, dz);
|
||||
}
|
||||
}
|
||||
|
||||
// normalize pressure
|
||||
for (bin = 0; bin < nbins1 * nbins2; bin++) {
|
||||
for (int bin = 0; bin < nbins1 * nbins2; bin++) {
|
||||
tdens[bin] *= invV;
|
||||
tpkxx[bin] *= invV;
|
||||
tpkyy[bin] *= invV;
|
||||
@ -386,9 +420,9 @@ void ComputeStressCartesian::compute_array()
|
||||
MPI_Allreduce(tpczz, pczz, nbins1 * nbins2, MPI_DOUBLE, MPI_SUM, world);
|
||||
|
||||
// populate array to output.
|
||||
for (bin = 0; bin < nbins1 * nbins2; bin++) {
|
||||
array[bin][0] = (bin % nbins1 + 0.5) * bin_width1;
|
||||
if (dims == 2) array[bin][1] = ((int) (bin / nbins1) + 0.5) * bin_width2;
|
||||
for (int bin = 0; bin < nbins1 * nbins2; bin++) {
|
||||
array[bin][0] = (bin % nbins1 + 0.5) * bin_width1 + boxlo[dir1];
|
||||
if (dims == 2) array[bin][1] = ((int) (bin / nbins1) + 0.5) * bin_width2 + boxlo[dir2];
|
||||
array[bin][0 + dims] = dens[bin];
|
||||
array[bin][1 + dims] = pkxx[bin];
|
||||
array[bin][2 + dims] = pkyy[bin];
|
||||
@ -402,25 +436,26 @@ void ComputeStressCartesian::compute_array()
|
||||
void ComputeStressCartesian::compute_pressure(double fpair, double xi, double yi, double delx,
|
||||
double dely, double delz)
|
||||
{
|
||||
int bin1, bin2, next_bin1, next_bin2;
|
||||
double la = 0.0, lb = 0.0, l_sum = 0.0;
|
||||
double rij[3] = {delx, dely, delz};
|
||||
double l1 = 0.0, l2, rij1, rij2;
|
||||
rij1 = rij[dir1];
|
||||
rij2 = rij[dir2];
|
||||
double rij1 = rij[dir1];
|
||||
double rij2 = rij[dir2];
|
||||
|
||||
next_bin1 = (int) floor(xi / bin_width1);
|
||||
next_bin2 = (int) floor(yi / bin_width2);
|
||||
int next_bin1 = (int) floor(xi / bin_width1);
|
||||
int next_bin2 = (int) floor(yi / bin_width2);
|
||||
|
||||
// Integrating along line
|
||||
while (lb < 1.0) {
|
||||
bin1 = next_bin1;
|
||||
bin2 = next_bin2;
|
||||
int bin1 = next_bin1;
|
||||
int bin2 = next_bin2;
|
||||
|
||||
double l1;
|
||||
if (rij1 > 0)
|
||||
l1 = ((bin1 + 1) * bin_width1 - xi) / rij1;
|
||||
else
|
||||
l1 = (bin1 * bin_width1 - xi) / rij1;
|
||||
|
||||
double l2;
|
||||
if (rij2 > 0)
|
||||
l2 = ((bin2 + 1) * bin_width2 - yi) / rij2;
|
||||
else
|
||||
|
||||
@ -36,6 +36,9 @@ class ComputeStressCartesian : public Compute {
|
||||
private:
|
||||
int nbins1, nbins2, dir1, dir2, dims;
|
||||
double bin_width1, bin_width2, invV;
|
||||
bool compute_ke = true;
|
||||
bool compute_pair = true;
|
||||
bool compute_bond = true;
|
||||
|
||||
// Number density, kinetic and configurational contribution to the pressure.
|
||||
double *dens, *pkxx, *pkyy, *pkzz, *pcxx, *pcyy, *pczz;
|
||||
|
||||
Reference in New Issue
Block a user