Merge pull request #2841 from donatas-surblys/heat-flux-virial-tally

A new tally compute to obtain virial heat flow from group 2 to group 1

doc/src/Commands_compute.rst

@@ -72,6 +72,7 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`gyration/shape/chunk <compute_gyration_shape_chunk>`
    * :doc:`heat/flux <compute_heat_flux>`
    * :doc:`heat/flux/tally <compute_tally>`
+   * :doc:`heat/flux/virial/tally <compute_tally>`
    * :doc:`hexorder/atom <compute_hexorder_atom>`
    * :doc:`hma <compute_hma>`
    * :doc:`improper <compute_improper>`

doc/src/compute.rst

@@ -208,7 +208,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` doc
 * :doc:`event/displace <compute_event_displace>` - detect event on atom displacement
 * :doc:`fabric <compute_fabric>` - calculates fabric tensors from pair interactions
 * :doc:`fep <compute_fep>` -
-* :doc:`force/tally <compute_tally>` -
+* :doc:`force/tally <compute_tally>` - force between two groups of atoms via the tally callback mechanism
 * :doc:`fragment/atom <compute_cluster_atom>` - fragment ID for each atom
 * :doc:`global/atom <compute_global_atom>` -
 * :doc:`group/group <compute_group_group>` - energy/force between two groups of atoms
@@ -217,7 +217,8 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` doc
 * :doc:`gyration/shape <compute_gyration_shape>` - shape parameters from gyration tensor
 * :doc:`gyration/shape/chunk <compute_gyration_shape_chunk>` - shape parameters from gyration tensor for each chunk
 * :doc:`heat/flux <compute_heat_flux>` - heat flux through a group of atoms
-* :doc:`heat/flux/tally <compute_tally>` -
+* :doc:`heat/flux/tally <compute_tally>` - heat flux through a group of atoms via the tally callback mechanism
+* :doc:`heat/flux/virial/tally <compute_tally>` - virial heat flux between two groups via the tally callback mechanism
 * :doc:`hexorder/atom <compute_hexorder_atom>` - bond orientational order parameter q6
 * :doc:`hma <compute_hma>` - harmonically mapped averaging for atomic crystals
 * :doc:`improper <compute_improper>` - energy of each improper sub-style
@@ -240,8 +241,8 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` doc
 * :doc:`pe <compute_pe>` - potential energy
 * :doc:`pe/atom <compute_pe_atom>` - potential energy for each atom
 * :doc:`mesont <compute_mesont>` - Nanotube bending,stretching, and intertube energies
-* :doc:`pe/mol/tally <compute_tally>` -
+* :doc:`pe/mol/tally <compute_tally>` - potential energy between two groups of atoms separated into intermolecular and intramolecular components via the tally callback mechanism
-* :doc:`pe/tally <compute_tally>` -
+* :doc:`pe/tally <compute_tally>` - potential energy between two groups of atoms via the tally callback mechanism
 * :doc:`plasticity/atom <compute_plasticity_atom>` - Peridynamic plasticity for each atom
 * :doc:`pressure <compute_pressure>` - total pressure and pressure tensor
 * :doc:`pressure/cylinder <compute_pressure_cylinder>` - pressure tensor in cylindrical coordinates
@@ -289,7 +290,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` doc
 * :doc:`stress/atom <compute_stress_atom>` - stress tensor for each atom
 * :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/tally <compute_tally>` -
+* :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
 * :doc:`temp/asphere <compute_temp_asphere>` - temperature of aspherical particles

doc/src/compute_tally.rst

@@ -1,5 +1,6 @@
 .. index:: compute force/tally
 .. index:: compute heat/flux/tally
+.. index:: compute heat/flux/virial/tally
 .. index:: compute pe/tally
 .. index:: compute pe/mol/tally
 .. index:: compute stress/tally
@@ -10,6 +11,9 @@ compute force/tally command
 compute heat/flux/tally command
 ===============================
 
+compute heat/flux/virial/tally command
+======================================
+
 compute pe/tally command
 ========================
 
@@ -27,7 +31,7 @@ Syntax
    compute ID group-ID style group2-ID
 
 * ID, group-ID are documented in :doc:`compute <compute>` command
-* style = *force/tally* or *pe/tally* or *pe/mol/tally* or *stress/tally*
+* style = *force/tally* or *heat/flux/tally* or *heat/flux/virial/tally* or * or *pe/tally* or *pe/mol/tally* or *stress/tally*
 * group2-ID = group ID of second (or same) group
 
 Examples
@@ -38,13 +42,17 @@ Examples
    compute 1 lower force/tally upper
    compute 1 left pe/tally right
    compute 1 lower stress/tally lower
+   compute 1 subregion heat/flux/tally all
+   compute 1 liquid heat/flux/virial/tally solid
 
 Description
 """""""""""
 
 Define a computation that calculates properties between two groups of
-atoms by accumulating them from pairwise non-bonded computations.  The
+atoms by accumulating them from pairwise non-bonded computations.
-two groups can be the same. This is similar to :doc:`compute group/group <compute_group_group>` only that the data is
+Except for *heat/flux/virial/tally*, the two groups can be the same.
+This is similar to :doc:`compute group/group <compute_group_group>`
+only that the data is
 accumulated directly during the non-bonded force computation. The
 computes *force/tally*\ , *pe/tally*\ , *stress/tally*\ , and
 *heat/flux/tally* are primarily provided as example how to program
@@ -57,6 +65,76 @@ the based classes of LAMMPS.
 
 ----------
 
+Compute *heat/flux/tally* obtains the heat flux
+(strictly speaking, heat flow) inside the first group,
+which is the sum of the convective contribution
+due to atoms in the first group and the virial contribution
+due to interaction between the first and second groups:
+
+.. math::
+
+   \mathbf{Q}=  \sum_{i \in \text{group 1}} e_i \mathbf{v}_i + \frac{1}{2} \sum_{i \in \text{group 1}} \sum_{\substack{j \in \text{group 2} \\ j \neq i } } \left( \mathbf{F}_{ij} \cdot \mathbf{v}_j \right) \mathbf{r}_{ij}
+
+When the second group in *heat/flux/tally* is set to "all",
+the resulting values will be identical
+to that obtained by :doc:`compute heat/flux <compute_heat_flux>`,
+provided only pairwise interactions exist.
+
+Compute *heat/flux/virial/tally* obtains the total virial heat flux
+(strictly speaking, heat flow) into the first group due to interaction
+with the second group, and is defined as:
+
+.. math::
+
+   Q = \frac{1}{2} \sum_{i \in \text{group 1}} \sum_{j \in \text{group 2}} \mathbf{F}_{ij} \cdot \left(\mathbf{v}_i + \mathbf{v}_j \right)
+
+Although, the *heat/flux/virial/tally* compute
+does not include the convective term,
+it can be used to obtain the total heat flux over control surfaces,
+when there are no particles crossing over,
+such as is often in solid-solid and solid-liquid interfaces.
+This would be identical to the method of planes method.
+Note that the *heat/flux/virial/tally* compute is distinctly different
+from the *heat/flux* and *heat/flux/tally* computes,
+that are essentially volume averaging methods.
+The following example demonstrates the difference:
+
+.. code-block:: LAMMPS
+
+   # System with only pairwise interactions.
+   # Non-periodic boundaries in the x direction.
+   # Has LeftLiquid and RightWall groups along x direction.
+
+   # Heat flux over the solid-liquid interface
+   compute hflow_hfvt LeftLiquid heat/flux/virial/tally RightWall
+   variable hflux_hfvt equal c_hflow_hfvt/(ly*lz)
+
+   # x component of approximate heat flux vector inside the liquid region,
+   # two approaches.
+   #
+   compute myKE all ke/atom
+   compute myPE all pe/atom
+   compute myStress all stress/atom NULL virial
+   compute hflow_hf LeftLiquid heat/flux myKE myPE myStress
+   variable hflux_hf equal c_hflow_hf[1]/${volLiq}
+   #
+   compute hflow_hft LeftLiquid heat/flux/tally all
+   variable hflux_hft equal c_hflow_hft[1]/${volLiq}
+
+   # Pressure over the solid-liquid interface, three approaches.
+   #
+   compute force_gg RightWall group/group LeftLiquid
+   variable press_gg equal c_force_gg[1]/(ly*lz)
+   #
+   compute force_ft RightWall force/tally LeftLiquid
+   compute rforce_ft RightWall reduce sum c_force_ft[1]
+   variable press_ft equal c_rforce_ft/(ly*lz)
+   #
+   compute rforce_hfvt all reduce sum c_hflow_hfvt[1]
+   variable press_hfvt equal -c_rforce_hfvt/(ly*lz)
+
+----------
+
 The pairwise contributions are computing via a callback that the
 compute registers with the non-bonded pairwise force computation.
 This limits the use to systems that have no bonds, no Kspace, and no
@@ -83,7 +161,17 @@ magnitude) and a per atom 3-element vector (force contribution from
 each atom).  Compute *stress/tally* calculates a global scalar
 (average of the diagonal elements of the stress tensor) and a per atom
 vector (the 6 elements of stress tensor contributions from the
-individual atom).
+individual atom). As in :doc:`compute heat/flux <compute_heat_flux>`,
+compute *heat/flux/tally* calculates a global vector of length 6,
+where the first 3 components are the :math:`x`, :math:`y`, :math:`z`
+components of the full heat flow vector,
+and the next 3 components are the corresponding components
+of just the convective portion of the flow, i.e. the
+first term in the equation for :math:`\mathbf{Q}`.
+Compute *heat/flux/virial/tally* calculates a global scalar (heat flow)
+and a per atom 3-element vector
+(contribution to the force acting over atoms in the first group
+from individual atoms in both groups).
 
 Both the scalar and vector values calculated by this compute are
 "extensive".

src/TALLY/compute_heat_flux_virial_tally.cpp

@@ -0,0 +1,238 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "compute_heat_flux_virial_tally.h"
+
+#include "atom.h"
+#include "comm.h"
+#include "error.h"
+#include "force.h"
+#include "group.h"
+#include "memory.h"
+#include "pair.h"
+#include "update.h"
+
+#include <cmath>
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeHeatFluxVirialTally::ComputeHeatFluxVirialTally(LAMMPS *lmp, int narg, char **arg) :
+    Compute(lmp, narg, arg)
+{
+  if (narg < 4) error->all(FLERR, "Illegal compute heat/flux/virial/tally command");
+
+  igroup2 = group->find(arg[3]);
+  if (igroup2 == -1)
+    error->all(FLERR, "Could not find compute heat/flux/virial/tally second group ID");
+  groupbit2 = group->bitmask[igroup2];
+
+  scalar_flag = 1;
+  vector_flag = 0;
+  peratom_flag = 1;
+  timeflag = 1;
+  dynamic_group_allow = 0;
+
+  comm_reverse = size_peratom_cols = 3;
+  extscalar = 1;
+  peflag = 1;    // we need Pair::ev_tally() to be run
+
+  did_setup = invoked_peratom = invoked_scalar = -1;
+  nmax = -1;
+  fatom = nullptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeHeatFluxVirialTally::~ComputeHeatFluxVirialTally()
+{
+  if (force && force->pair) force->pair->del_tally_callback(this);
+  memory->destroy(fatom);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeHeatFluxVirialTally::init()
+{
+  if (force->pair == nullptr)
+    error->all(FLERR, "Trying to use compute heat/flux/virial/tally without pair style");
+  else
+    force->pair->add_tally_callback(this);
+
+  if (comm->me == 0) {
+    if (force->pair->single_enable == 0 || force->pair->manybody_flag)
+      error->warning(FLERR, "Compute heat/flux/virial/tally used with incompatible pair style");
+
+    if (force->bond || force->angle || force->dihedral || force->improper || force->kspace)
+      error->warning(FLERR, "Compute heat/flux/virial/tally only called from pair style");
+  }
+
+  // error out if any atoms are in both groups
+  for (int i = 0; i < atom->nlocal; i++) {
+    if ((atom->mask[i] & groupbit) && (atom->mask[i] & groupbit2)) {
+      if (atom->tag_enable) {
+        error->all(FLERR,
+                   "Atom {} belonging to both groups is not allowed"
+                   " with compute heat/flux/virial/tally",
+                   atom->tag[i]);
+      } else {
+        error->all(FLERR,
+                   "Atom belonging to both groups is not allowed"
+                   " with compute heat/flux/virial/tally");
+      }
+    }
+  }
+
+  did_setup = -1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeHeatFluxVirialTally::pair_setup_callback(int, int)
+{
+  // run setup only once per time step.
+  // we may be called from multiple pair styles
+
+  if (did_setup == update->ntimestep) return;
+
+  const int ntotal = atom->nlocal + atom->nghost;
+
+  // grow per-atom storage, if needed
+
+  if (atom->nmax > nmax) {
+    memory->destroy(fatom);
+    nmax = atom->nmax;
+    memory->create(fatom, nmax, size_peratom_cols, "heat/flux/virial/tally:fatom");
+    array_atom = fatom;
+  }
+
+  // clear storage
+
+  for (int i = 0; i < ntotal; ++i)
+    for (int j = 0; j < size_peratom_cols; ++j) fatom[i][j] = 0.0;
+
+  did_setup = update->ntimestep;
+}
+
+/* ---------------------------------------------------------------------- */
+void ComputeHeatFluxVirialTally::pair_tally_callback(int i, int j, int nlocal, int newton, double,
+                                                     double, double fpair, double dx, double dy,
+                                                     double dz)
+{
+  const int *const mask = atom->mask;
+
+  const bool ingroup1 = (mask[i] & groupbit);
+  if ((ingroup1 && (mask[j] & groupbit2)) || ((mask[i] & groupbit2) && (mask[j] & groupbit))) {
+
+    // signs set to calculate heat flux from group2 into group1
+    const double fx = (ingroup1 ? 0.5 : -0.5) * dx * fpair;
+    const double fy = (ingroup1 ? 0.5 : -0.5) * dy * fpair;
+    const double fz = (ingroup1 ? 0.5 : -0.5) * dz * fpair;
+
+    if (newton || i < nlocal) {
+      fatom[i][0] += fx;
+      fatom[i][1] += fy;
+      fatom[i][2] += fz;
+    }
+    if (newton || j < nlocal) {
+      fatom[j][0] += fx;
+      fatom[j][1] += fy;
+      fatom[j][2] += fz;
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int ComputeHeatFluxVirialTally::pack_reverse_comm(int n, int first, double *buf)
+{
+  int i, m, last;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) {
+    buf[m++] = fatom[i][0];
+    buf[m++] = fatom[i][1];
+    buf[m++] = fatom[i][2];
+  }
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeHeatFluxVirialTally::unpack_reverse_comm(int n, int *list, double *buf)
+{
+  int i, j, m;
+
+  m = 0;
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    fatom[j][0] += buf[m++];
+    fatom[j][1] += buf[m++];
+    fatom[j][2] += buf[m++];
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double ComputeHeatFluxVirialTally::compute_scalar()
+{
+  // need to collect contributions from ghost atoms
+  // because we need to use local velocities to compute heat flux
+  if (invoked_peratom != update->ntimestep) compute_peratom();
+
+  invoked_scalar = update->ntimestep;
+  if ((did_setup != invoked_scalar) || (update->eflag_global != invoked_scalar))
+    error->all(FLERR, "Energy was not tallied on needed timestep");
+
+  // sum heat flux across procs
+  double hflux = 0.0;
+  for (int i = 0; i < atom->nlocal; i++)
+    hflux +=
+        fatom[i][0] * atom->v[i][0] + fatom[i][1] * atom->v[i][1] + fatom[i][2] * atom->v[i][2];
+
+  MPI_Allreduce(&hflux, &scalar, 1, MPI_DOUBLE, MPI_SUM, world);
+
+  return scalar;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeHeatFluxVirialTally::compute_peratom()
+{
+  invoked_peratom = update->ntimestep;
+  if ((did_setup != invoked_peratom) || (update->eflag_global != invoked_peratom))
+    error->all(FLERR, "Energy was not tallied on needed timestep");
+
+  // collect contributions from ghost atoms
+
+  if (force->newton_pair) {
+    comm->reverse_comm_compute(this);
+
+    // clear out ghost atom data after it has been collected to local atoms
+    const int nall = atom->nlocal + atom->nghost;
+    for (int i = atom->nlocal; i < nall; ++i)
+      for (int j = 0; j < size_peratom_cols; ++j) fatom[i][j] = 0.0;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double ComputeHeatFluxVirialTally::memory_usage()
+{
+  double bytes = (nmax < 0) ? 0 : nmax * size_peratom_cols * sizeof(double);
+  return bytes;
+}

src/TALLY/compute_heat_flux_virial_tally.h

@@ -0,0 +1,64 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+// clang-format off
+ComputeStyle(heat/flux/virial/tally,ComputeHeatFluxVirialTally);
+// clang-format on
+#else
+
+#ifndef LMP_COMPUTE_HEAT_FLUX_VIRIAL_TALLY_H
+#define LMP_COMPUTE_HEAT_FLUX_VIRIAL_TALLY_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeHeatFluxVirialTally : public Compute {
+
+ public:
+  ComputeHeatFluxVirialTally(class LAMMPS *, int, char **);
+  virtual ~ComputeHeatFluxVirialTally();
+
+  void init();
+
+  double compute_scalar();
+  void compute_peratom();
+
+  int pack_reverse_comm(int, int, double *);
+  void unpack_reverse_comm(int, int *, double *);
+  double memory_usage();
+
+  void pair_setup_callback(int, int);
+  void pair_tally_callback(int, int, int, int, double, double, double, double, double, double);
+
+ private:
+  bigint did_setup;
+  int nmax, igroup2, groupbit2;
+  double **fatom;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+*/