diff --git a/doc/src/compute_tally.rst b/doc/src/compute_tally.rst
index 32e3e31030..d2e4569820 100644
--- a/doc/src/compute_tally.rst
+++ b/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
-two groups can be the same. This is similar to :doc:`compute group/group <compute_group_group>` only that the data is
+atoms by accumulating them from pairwise non-bonded computations.
+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 flux vector,
+and the next 3 components are the corresponding components
+of just the convective portion of the flux, 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".