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 new file mode 100644 index 0000000000..3e1e9d9a8a --- /dev/null +++ b/doc/src/compute_stress_cartesian.rst @@ -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 ` 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. + +.. 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 `, + +.. 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 +: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 +: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..edd1df9217 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 """"""" diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index 5662fc76d4..ec654dc61c 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -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 `. 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 """"""""""" @@ -143,9 +148,9 @@ Related commands :doc:`compute stress/atom `, :doc:`compute pressure `, -:doc:`compute stress/cartesian `, -:doc:`compute stress/cylinder `, -:doc:`compute stress/spherical ` +:doc:`compute stress/cartesian `, +:doc:`compute stress/cylinder `, +:doc:`compute stress/spherical ` Default """"""" diff --git a/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp b/src/EXTRA-COMPUTE/compute_stress_cartesian.cpp index 53ff70c561..3835b8bea9 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" @@ -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 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;