tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

more thorough tests and PBC handling for compute stress/cartesian

Browse Source
This commit is contained in:
Axel Kohlmeyer
2023-05-04 03:46:18 -04:00
parent de45437cc9
commit a0d9854e11
2 changed files with 60 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
doc/src/compute_stress_profile.rst
View File

@ -136,12 +136,12 @@ More information on the similarities and differences can be found in
Restrictions
""""""""""""
These computes calculate the stress tensor contributions for pair
styles only (i.e., no bond, angle, dihedral, etc. contributions, and in
the presence of bonded interactions, the result will be incorrect due to
exclusions for special bonds) and requires pairwise force calculations
not available for most many-body pair styles.
Note that :math:`k`-space calculations are also excluded.
These computes calculate the stress tensor contributions for pair styles
only (i.e., no bond, angle, dihedral, etc. contributions, and in the
presence of bonded interactions, the result may be incorrect due to
exclusions for :doc:`special bonds <special_bonds>` excluding pairs of atoms
completely). It requires pairwise force calculations not available for most
many-body pair styles. Note that :math:`k`-space calculations are also excluded.
These computes are part of the EXTRA-COMPUTE package. They are only
enabled if LAMMPS was built with that package. See the :doc:`Build

62
src/EXTRA-COMPUTE/compute_stress_cartesian.cpp
View File

@ -18,8 +18,10 @@
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "fix.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "neigh_list.h"
#include "neighbor.h"
#include "pair.h"
@ -57,7 +59,6 @@ ComputeStressCartesian::ComputeStressCartesian(LAMMPS *lmp, int narg, char **arg
pcxx(nullptr), pcyy(nullptr), pczz(nullptr), tdens(nullptr), tpkxx(nullptr), tpkyy(nullptr),
tpkzz(nullptr), tpcxx(nullptr), tpcyy(nullptr), tpczz(nullptr), list(nullptr)
{
if (lmp->citeme) lmp->citeme->add(cite_compute_stress_cartesian);
// narg == 5 for one-dimensional and narg == 7 for two-dimensional
@ -75,7 +76,7 @@ ComputeStressCartesian::ComputeStressCartesian(LAMMPS *lmp, int narg, char **arg
else if (strcmp(arg[3], "z") == 0)
dir1 = 2;
else
error->all(FLERR, "Illegal compute stress/cartesian command.");
error->all(FLERR, "Illegal compute stress/cartesian direction: {}", arg[3]);
dir2 = 0;
bin_width1 = utils::numeric(FLERR, arg[4], false, lmp);
@ -83,17 +84,20 @@ ComputeStressCartesian::ComputeStressCartesian(LAMMPS *lmp, int narg, char **arg
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;
if ((fabs(tmp_binwidth - bin_width1) > SMALL) && (comm->me == 0))
utils::logmesg(lmp, "Adjusting second bin width for compute {} from {:.6f} to {:.6f}\n", style,
utils::logmesg(lmp, "Adjusting first bin width for compute {} from {:.6f} to {:.6f}\n", style,
bin_width1, tmp_binwidth);
bin_width1 = tmp_binwidth;
if (bin_width1 <= 0.0)
error->all(FLERR, "Illegal compute stress/cartesian command. Bin width must be > 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. Bin width larger than box.");
error->all(FLERR, "Illegal compute stress/cartesian command. First bin width > box.");
invV = bin_width1;
if (dims == 2) {
@ -104,7 +108,7 @@ ComputeStressCartesian::ComputeStressCartesian(LAMMPS *lmp, int narg, char **arg
else if (strcmp(arg[5], "z") == 0)
dir2 = 2;
else
error->all(FLERR, "Illegal compute stress/cartesian command.");
error->all(FLERR, "Illegal compute stress/cartesian direction {}", arg[5]);
bin_width2 = utils::numeric(FLERR, arg[6], false, lmp);
nbins2 = (int) ((domain->boxhi[dir2] - domain->boxlo[dir2]) / bin_width2);
@ -119,11 +123,32 @@ ComputeStressCartesian::ComputeStressCartesian(LAMMPS *lmp, int narg, char **arg
invV *= bin_width2;
if (bin_width2 <= 0.0)
error->all(FLERR, "Illegal compute stress/cartesian command. Bin width must be > 0");
error->all(FLERR, "Illegal compute stress/cartesian command. Second bin width must be > 0");
else if (bin_width2 > domain->boxhi[dir2] - domain->boxlo[dir2])
error->all(FLERR, "Illegal compute stress/cartesian command. Bin width larger than box");
error->all(FLERR, "Illegal compute stress/cartesian command. Second bin width > box");
}
// check for overflow
if (((bigint) nbins1 * nbins2) > MAXSMALLINT)
error->all(FLERR, "Too many bins in compute stress/cartesian");
// check for variable box dimension
int box_incompatible = 0;
for (auto ifix : modify->get_fix_list()) {
if (((dir1 == 0) && (ifix->box_change & Fix::BOX_CHANGE_X)) ||
((dir1 == 1) && (ifix->box_change & Fix::BOX_CHANGE_Y)) ||
((dir1 == 2) && (ifix->box_change & Fix::BOX_CHANGE_Z)))
box_incompatible = 1;
if (dims == 2) {
if (((dir2 == 0) && (ifix->box_change & Fix::BOX_CHANGE_X)) ||
((dir2 == 1) && (ifix->box_change & Fix::BOX_CHANGE_Y)) ||
((dir2 == 2) && (ifix->box_change & Fix::BOX_CHANGE_Z)))
box_incompatible = 1;
}
}
if (box_incompatible)
error->all(FLERR, "Must not use compute stress/cartesian on variable box dimension");
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];
@ -248,6 +273,27 @@ void ComputeStressCartesian::compute_array()
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;
else if (bin1 >= nbins1)
bin1 = (bin1 - nbins1) % nbins1;
} else if (bin1 < 0)
bin1 = 0;
else if (bin1 >= nbins1)
bin1 = nbins1 - 1;
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 = nbins2 - 1;
j = bin1 + bin2 * nbins1;
tdens[j] += 1;
tpkxx[j] += mass[type[i]] * v[i][0] * v[i][0];