From ecca46acf9fd283e2c40e57a3135618c9108a9eb Mon Sep 17 00:00:00 2001 From: Lars Veldscholte Date: Tue, 9 May 2023 17:11:06 +0200 Subject: [PATCH 1/8] Include bond interactions in force --- .../compute_stress_cartesian.cpp | 33 +++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp b/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp index 53ff70c561..d2b9fe7bd7 100644 --- a/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp @@ -14,6 +14,7 @@ #include "compute_stress_cartesian.h" #include "atom.h" +#include "bond.h" #include "citeme.h" #include "comm.h" #include "domain.h" @@ -365,6 +366,38 @@ void ComputeStressCartesian::compute_array() } } + // Loop over all bonds + 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 type = neighbor->bondlist[i_bond][2]; + + // Skip if one of both atoms is not in group + if (!(atom->mask[i] & groupbit)) continue; + if (!(atom->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 >= atom->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(type, rsq, i, j, fbond); + compute_pressure(fbond, xi, yi, dx, dy, dz); + } + // normalize pressure for (bin = 0; bin < nbins1 * nbins2; bin++) { tdens[bin] *= invV; From e246864682f00655749a3a251965f82a81dae3d3 Mon Sep 17 00:00:00 2001 From: Lars Veldscholte Date: Wed, 24 May 2023 16:17:20 +0200 Subject: [PATCH 2/8] Refactor compute_array() and compute_pressure(): Remove unnecessary copies of variables, declare variables locally so they are properly scoped --- .../compute_stress_cartesian.cpp | 136 +++++++----------- 1 file changed, 52 insertions(+), 84 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp b/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp index d2b9fe7bd7..2f6e75b528 100644 --- a/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp @@ -230,35 +230,16 @@ 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; tagint *tag = atom->tag; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - double *special_coul = force->special_coul; - double *special_lj = force->special_lj; - int newton_pair = force->newton_pair; double *boxlo = domain->boxlo; // 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; - // 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; @@ -269,9 +250,9 @@ 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; + for (int i = 0; i < atom->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 @@ -295,73 +276,59 @@ void ComputeStressCartesian::compute_array() else if (bin2 >= nbins2) bin2 = nbins2 - 1; - j = bin1 + bin2 * nbins1; + 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]; + tpkxx[j] += atom->mass[atom->type[i]] * v[i][0] * v[i][0]; + tpkyy[j] += atom->mass[atom->type[i]] * v[i][1] * v[i][1]; + tpkzz[j] += atom->mass[atom->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 + for (int ii = 0; ii < list->inum; ii++) { + int i = list->ilist[ii]; - Pair *pair = force->pair; - double **cutsq = force->pair->cutsq; + // skip if I or J are not in group + if (!(atom->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; - - 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)]; + for (int jj = 0; jj < list->numneigh[i]; jj++) { + int j = list->firstneigh[i][jj]; + double factor_lj = force->special_lj[sbmask(j)]; + double factor_coul = force->special_coul[sbmask(j)]; j &= NEIGHMASK; - if (!(mask[j] & groupbit)) continue; + if (!(atom->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; + // 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 (force->newton_pair == 0 && j >= atom->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 { - 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; + // 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; } } } - 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]; + 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 - if (rsq >= cutsq[itype][jtype]) continue; - pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair); + int itype = atom->type[i]; + int jtype = atom->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); } } @@ -399,7 +366,7 @@ void ComputeStressCartesian::compute_array() } // 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; @@ -419,7 +386,7 @@ 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++) { + for (int 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; array[bin][0 + dims] = dens[bin]; @@ -435,25 +402,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 From ddc34e03d66452402993ec20203b12a089f2fe21 Mon Sep 17 00:00:00 2001 From: Lars Veldscholte Date: Wed, 31 May 2023 09:50:50 +0200 Subject: [PATCH 3/8] Revert removal of copies of pointers --- .../compute_stress_cartesian.cpp | 52 ++++++++++++------- 1 file changed, 32 insertions(+), 20 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp b/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp index 2f6e75b528..106ffef8fb 100644 --- a/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp @@ -232,12 +232,24 @@ void ComputeStressCartesian::compute_array() { double **x = atom->x; double **v = atom->v; + double *mass = atom->mass; tagint *tag = atom->tag; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; double *boxlo = domain->boxlo; // invoke half neighbor list (will copy or build if necessary) neighbor->build_one(list); + int inum = list->inum; + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + // Zero arrays for (int bin = 0; bin < nbins1 * nbins2; bin++) { tdens[bin] = 0; @@ -250,7 +262,7 @@ void ComputeStressCartesian::compute_array() } // calculate number density and kinetic contribution to pressure - for (int i = 0; i < atom->nlocal; i++) { + 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; @@ -278,31 +290,31 @@ void ComputeStressCartesian::compute_array() int j = bin1 + bin2 * nbins1; tdens[j] += 1; - tpkxx[j] += atom->mass[atom->type[i]] * v[i][0] * v[i][0]; - tpkyy[j] += atom->mass[atom->type[i]] * v[i][1] * v[i][1]; - tpkzz[j] += atom->mass[atom->type[i]] * v[i][2] * v[i][2]; + 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 - for (int ii = 0; ii < list->inum; ii++) { - int i = list->ilist[ii]; + for (int ii = 0; ii < inum; ii++) { + int i = ilist[ii]; // skip if I or J are not in group - if (!(atom->mask[i] & groupbit)) continue; + if (!(mask[i] & groupbit)) continue; double xi1 = x[i][dir1] - boxlo[dir1]; double xi2 = x[i][dir2] - boxlo[dir2]; - for (int jj = 0; jj < list->numneigh[i]; jj++) { - int j = list->firstneigh[i][jj]; - double factor_lj = force->special_lj[sbmask(j)]; - double factor_coul = force->special_coul[sbmask(j)]; + 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 (!(atom->mask[j] & groupbit)) continue; + if (!(mask[j] & groupbit)) 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 (force->newton_pair == 0 && j >= atom->nlocal) { + 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]) { @@ -323,8 +335,8 @@ void ComputeStressCartesian::compute_array() double rsq = delx * delx + dely * dely + delz * delz; // Check if inside cut-off - int itype = atom->type[i]; - int jtype = atom->type[j]; + int itype = type[i]; + int jtype = type[j]; if (rsq >= force->pair->cutsq[itype][jtype]) continue; double fpair; @@ -338,14 +350,14 @@ void ComputeStressCartesian::compute_array() // i == atom1, j == atom2 int i = neighbor->bondlist[i_bond][0]; int j = neighbor->bondlist[i_bond][1]; - int type = neighbor->bondlist[i_bond][2]; + int btype = neighbor->bondlist[i_bond][2]; // Skip if one of both atoms is not in group - if (!(atom->mask[i] & groupbit)) continue; - if (!(atom->mask[j] & groupbit)) continue; + 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 >= atom->nlocal) { + 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]) { @@ -361,7 +373,7 @@ void ComputeStressCartesian::compute_array() double yi = x[i][dir2] - boxlo[dir2]; double fbond; - force->bond->single(type, rsq, i, j, fbond); + force->bond->single(btype, rsq, i, j, fbond); compute_pressure(fbond, xi, yi, dx, dy, dz); } From 045b230587c3fcf1eb4ea7a0fcf00f0f08b1d19c Mon Sep 17 00:00:00 2001 From: Lars Veldscholte Date: Wed, 7 Jun 2023 10:43:57 +0200 Subject: [PATCH 4/8] Fix shifted coordinates: Add `boxlo` to the bin centers --- src/EXTRA-COMPUTE/compute_stress_cartesian.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp b/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp index 106ffef8fb..0aae797111 100644 --- a/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp @@ -399,8 +399,8 @@ void ComputeStressCartesian::compute_array() // populate array to output. for (int 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; + 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]; From f2f8e139d8c23e750ce603ef4cd2a95b245bdb5a Mon Sep 17 00:00:00 2001 From: Lars Veldscholte Date: Wed, 7 Jun 2023 16:03:32 +0200 Subject: [PATCH 5/8] Add optional keywords to arguments for ke/pair/bond forces --- .../compute_stress_cartesian.cpp | 232 ++++++++++-------- src/EXTRA-COMPUTE/compute_stress_cartesian.h | 3 + 2 files changed, 130 insertions(+), 105 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp b/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp index 0aae797111..3835b8bea9 100644 --- a/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp @@ -62,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) @@ -79,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; @@ -94,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) @@ -150,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]; @@ -262,119 +278,125 @@ void ComputeStressCartesian::compute_array() } // calculate number density and kinetic contribution to pressure - 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; + 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; - 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]; + 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 - for (int ii = 0; ii < inum; ii++) { - int i = ilist[ii]; + if (compute_pair && force->pair) { + for (int ii = 0; ii < inum; ii++) { + int i = ilist[ii]; - // skip if I or J are not in group - if (!(mask[i] & groupbit)) continue; + // skip if I or J are not in group + if (!(mask[i] & groupbit)) continue; - double xi1 = x[i][dir1] - boxlo[dir1]; - double xi2 = x[i][dir2] - boxlo[dir2]; + double xi1 = x[i][dir1] - boxlo[dir1]; + double xi2 = x[i][dir2] - boxlo[dir2]; - 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; + 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; - // 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; + // 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); } - - 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); } } // Loop over all bonds - 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]; + 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]; - // Skip if one of both atoms is not in group - if (!(mask[i] & groupbit)) continue; - if (!(mask[j] & groupbit)) continue; + // 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; + // 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); } - - 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 diff --git a/src/EXTRA-COMPUTE/compute_stress_cartesian.h b/src/EXTRA-COMPUTE/compute_stress_cartesian.h index 0954d9ca71..40f2e9d4af 100644 --- a/src/EXTRA-COMPUTE/compute_stress_cartesian.h +++ b/src/EXTRA-COMPUTE/compute_stress_cartesian.h @@ -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; From c851c7304c02d8b38424c240e3ee84671a402f7a Mon Sep 17 00:00:00 2001 From: Lars Veldscholte Date: Thu, 8 Jun 2023 11:03:21 +0200 Subject: [PATCH 6/8] Update documentation for compute stress/cartesian, and split the doc page compute_stress_profile into compute_stress_cartesian and compute_stress_curvilinear --- doc/src/compute_stress_cartesian.rst | 118 ++++++++++++++++++ ...ile.rst => compute_stress_curvilinear.rst} | 48 ++----- 2 files changed, 125 insertions(+), 41 deletions(-) create mode 100644 doc/src/compute_stress_cartesian.rst rename doc/src/{compute_stress_profile.rst => compute_stress_curvilinear.rst} (69%) diff --git a/doc/src/compute_stress_cartesian.rst b/doc/src/compute_stress_cartesian.rst new file mode 100644 index 0000000000..f52b360aa6 --- /dev/null +++ b/doc/src/compute_stress_cartesian.rst @@ -0,0 +1,118 @@ +.. 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 ` 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)`. 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. + +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 `, + +.. 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 `. The number density +values are in inverse volume :doc:`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 ` +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) `. +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, 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)`. + +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 ` doc page for more info. + +Related commands +"""""""""""""""" + +:doc:`compute stress/atom `, :doc:`compute pressure `, +:doc:`compute stress/mop/profile `, :doc:`compute stress/spherical `, +:doc:`compute stress/cylinder ` + +---------- + +.. _Ikeshoji2: + +**(Ikeshoji)** Ikeshoji, Hafskjold, Furuholt, Mol Sim, 29, 101-109, (2003). \ No newline at end of file diff --git a/doc/src/compute_stress_profile.rst b/doc/src/compute_stress_curvilinear.rst similarity index 69% rename from doc/src/compute_stress_profile.rst rename to doc/src/compute_stress_curvilinear.rst index 7a4b6a38e6..f7cbf5efbe 100644 --- a/doc/src/compute_stress_profile.rst +++ b/doc/src/compute_stress_curvilinear.rst @@ -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 ` 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)`. 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)`. The compute *stress/spherical* computes the stress profile along the radial direction in spherical @@ -79,17 +65,6 @@ coordinates, as described in :ref:`(Ikeshoji)`. 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 `, .. 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 ` 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) `. -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. -More information on the similarities and differences can be found in -:ref:`(Ikeshoji)`. - Restrictions """""""""""" @@ -150,7 +115,8 @@ package ` doc page for more info. Related commands """""""""""""""" -:doc:`compute stress/atom `, :doc:`compute pressure `, :doc:`compute stress/mop/profile ` +:doc:`compute stress/atom `, :doc:`compute pressure `, :doc:`compute +stress/mop/profile `, :doc:`compute stress/cartesian ` Default """"""" From 8aeb059ce867b279f1b30fe359475b6fc1b3b243 Mon Sep 17 00:00:00 2001 From: Lars Veldscholte Date: Wed, 14 Jun 2023 15:03:35 +0200 Subject: [PATCH 7/8] Update doc/src/compute_stress_cartesian.rst Co-authored-by: Axel Kohlmeyer --- doc/src/compute_stress_cartesian.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/doc/src/compute_stress_cartesian.rst b/doc/src/compute_stress_cartesian.rst index f52b360aa6..b9f090ef58 100644 --- a/doc/src/compute_stress_cartesian.rst +++ b/doc/src/compute_stress_cartesian.rst @@ -41,6 +41,10 @@ split into a kinetic contribution :math:`P^k` and a virial contribution 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 From e6cd79e0e928071fbff2617ac9fd172657acc5d3 Mon Sep 17 00:00:00 2001 From: Lars Veldscholte Date: Wed, 14 Jun 2023 15:04:35 +0200 Subject: [PATCH 8/8] Fix doc links --- doc/src/Commands_compute.rst | 6 +++--- doc/src/compute.rst | 6 +++--- doc/src/compute_stress_cartesian.rst | 6 +++--- doc/src/compute_stress_curvilinear.rst | 4 ++-- doc/src/compute_stress_mop.rst | 22 +++++++++++++++------- 5 files changed, 26 insertions(+), 18 deletions(-) diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 755000c976..9066531459 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -153,11 +153,11 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`sph/t/atom ` * :doc:`spin ` * :doc:`stress/atom ` - * :doc:`stress/cartesian ` - * :doc:`stress/cylinder ` + * :doc:`stress/cartesian ` + * :doc:`stress/cylinder ` * :doc:`stress/mop ` * :doc:`stress/mop/profile ` - * :doc:`stress/spherical ` + * :doc:`stress/spherical ` * :doc:`stress/tally ` * :doc:`tdpd/cc/atom ` * :doc:`temp (k) ` diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 950487cb72..4e5d8e6cb4 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -307,11 +307,11 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`sph/t/atom ` - per-atom internal temperature of Smooth-Particle Hydrodynamics atoms * :doc:`spin ` - magnetic quantities for a system of atoms having spins * :doc:`stress/atom ` - stress tensor for each atom -* :doc:`stress/cartesian ` - stress tensor in cartesian coordinates -* :doc:`stress/cylinder ` - stress tensor in cylindrical coordinates +* :doc:`stress/cartesian ` - stress tensor in cartesian coordinates +* :doc:`stress/cylinder ` - stress tensor in cylindrical coordinates * :doc:`stress/mop ` - normal components of the local stress tensor using the method of planes * :doc:`stress/mop/profile ` - profile of the normal components of the local stress tensor using the method of planes -* :doc:`stress/spherical ` - stress tensor in spherical coordinates +* :doc:`stress/spherical ` - stress tensor in spherical coordinates * :doc:`stress/tally ` - stress between two groups of atoms via the tally callback mechanism * :doc:`tdpd/cc/atom ` - per-atom chemical concentration of a specified species for each tDPD particle * :doc:`temp ` - temperature of group of atoms diff --git a/doc/src/compute_stress_cartesian.rst b/doc/src/compute_stress_cartesian.rst index b9f090ef58..3e1e9d9a8a 100644 --- a/doc/src/compute_stress_cartesian.rst +++ b/doc/src/compute_stress_cartesian.rst @@ -83,13 +83,13 @@ 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) `. -A key difference -is that compute `stress/mop/profile ` +A key difference is that compute +:doc:`stress/mop/profile ` 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 +*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. diff --git a/doc/src/compute_stress_curvilinear.rst b/doc/src/compute_stress_curvilinear.rst index f7cbf5efbe..edd1df9217 100644 --- a/doc/src/compute_stress_curvilinear.rst +++ b/doc/src/compute_stress_curvilinear.rst @@ -115,8 +115,8 @@ package ` doc page for more info. Related commands """""""""""""""" -:doc:`compute stress/atom `, :doc:`compute pressure `, :doc:`compute -stress/mop/profile `, :doc:`compute stress/cartesian ` +:doc:`compute stress/atom `, :doc:`compute pressure `, +:doc:`compute stress/mop/profile `, :doc:`compute stress/cartesian ` Default """"""" diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index 4ad2261bb0..0b31805fbf 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -76,12 +76,16 @@ command, since those are contributions to the global system pressure. NOTE 3: The local stress profile generated by compute *stress/mop/profile* is similar to that obtained by compute -:doc:`stress/cartesian `. -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)`. +:doc:`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 """"""""""" @@ -123,7 +127,11 @@ intra-molecular interactions, and long range (kspace) interactions. Related commands """""""""""""""" -:doc:`compute stress/atom `, :doc:`compute pressure `, :doc:`compute stress/cartesian `, :doc:`compute stress/cylinder `, :doc:`compute stress/spherical ` +:doc:`compute stress/atom `, +:doc:`compute pressure `, +:doc:`compute stress/cartesian `, +:doc:`compute stress/cylinder `, +:doc:`compute stress/spherical ` Default """""""