diff --git a/doc/src/compute_heat_flux.rst b/doc/src/compute_heat_flux.rst
index 9642582c5e..df98c85b64 100644
--- a/doc/src/compute_heat_flux.rst
+++ b/doc/src/compute_heat_flux.rst
@@ -54,38 +54,58 @@ third calculates per-atom stress (\ *stress-ID*\ ).
    (or any group whose atoms are superset of the atoms in this compute's
    group).  LAMMPS does not check for this.
 
-The Green-Kubo formulas relate the ensemble average of the
-auto-correlation of the heat flux J to the thermal conductivity kappa:
+In case of pairwise interactions the heat flux is defined as:
 
-.. image:: Eqs/heat_flux_J.jpg
-   :align: center
+.. math::
+   \mathbf{J} &= \frac{1}{V} \left[ \sum_i e_i \mathbf{v}_i - \sum_{i} \mathbf{S}_{i} \mathbf{v}_i \right] \\
+   &= \frac{1}{V} \left[ \sum_i e_i \mathbf{v}_i + \sum_{i<j} \left( \mathbf{F}_{ij} \cdot \mathbf{v}_j \right) \mathbf{r}_{ij} \right] \\
+   &= \frac{1}{V} \left[ \sum_i e_i \mathbf{v}_i + \frac{1}{2} \sum_{i<j} \left( \mathbf{F}_{ij} \cdot \left(\mathbf{v}_i + \mathbf{v}_j \right)  \right) \mathbf{r}_{ij} \right]
 
-.. image:: Eqs/heat_flux_k.jpg
-   :align: center
-
-Ei in the first term of the equation for J is the per-atom energy
-(potential and kinetic).  This is calculated by the computes *ke-ID*
-and *pe-ID*\ .  Si in the second term of the equation for J is the
-per-atom stress tensor calculated by the compute *stress-ID*\ .  The
-tensor multiplies Vi as a 3x3 matrix-vector multiply to yield a
-vector.  Note that as discussed below, the 1/V scaling factor in the
-equation for J is NOT included in the calculation performed by this
-compute; you need to add it for a volume appropriate to the atoms
+:math:`e_i` in the first term of the equation
+is the per-atom energy (potential and kinetic).
+This is calculated by the computes *ke-ID*
+and *pe-ID*. :math:`\mathbf{S}_i` in the second term is the
+per-atom stress tensor calculated by the compute *stress-ID*.
+See :doc:`compute stress/atom <compute_stress_atom>`
+and :doc:`compute centroid/stress/atom <compute_stress_atom>`
+for possible definitions of atomic stress :math:`\mathbf{S}_i`
+in the case of thee-body and larger many-body interactions.
+The tensor multiplies :math:`\mathbf{v}_i` as a 3x3 matrix-vector multiply
+to yield a vector.
+Note that as discussed below, the :math:`\frac{1}{V}` scaling factor in the
+equation for :math:`\mathbf{J}` is NOT included in the calculation performed by
+these computes; you need to add it for a volume appropriate to the atoms
 included in the calculation.
 
 .. note::
 
-   The :doc:`compute pe/atom <compute_pe_atom>` and :doc:`compute stress/atom <compute_stress_atom>` commands have options for which
+   The :doc:`compute pe/atom <compute_pe_atom>` and
+   :doc:`compute stress/atom <compute_stress_atom>`
+   commands have options for which
    terms to include in their calculation (pair, bond, etc).  The heat
-   flux calculation will thus include exactly the same terms.  Normally
+   flux calculation will thus include exactly the same terms. Normally
    you should use :doc:`compute stress/atom virial <compute_stress_atom>`
+   or :doc:`compute centroid/stress/atom virial <compute_stress_atom>`
    so as not to include a kinetic energy term in the heat flux.
 
-This compute calculates 6 quantities and stores them in a 6-component
-vector.  The first 3 components are the x, y, z components of the full
-heat flux vector, i.e. (Jx, Jy, Jz).  The next 3 components are the x,
-y, z components of just the convective portion of the flux, i.e. the
-first term in the equation for J above.
+
+.. warning::
+
+   The compute *heat/flux* has been reported to produce unphysical
+   values for three and larger many-body interactions such
+   as angles, dihedrals and torsions,
+   when used with :doc:`compute stress/atom <compute_stress_atom>`,
+   as discussed in :ref:`(Surblys) <Surblys2>` and :ref:`(Boone) <Boone>`.
+   You are strongly advised to
+   use :doc:`compute centroid/stress/atom <compute_stress_atom>`,
+   which has been implemented specifically for such cases.
+
+The Green-Kubo formulas relate the ensemble average of the
+auto-correlation of the heat flux :math:`\mathbf{J}`
+to the thermal conductivity :math:`\kappa`:
+
+.. math::
+   \kappa  = \frac{V}{k_B T^2} \int_0^\infty \langle J_x(0)  J_x(t) \rangle \, \mathrm{d} t = \frac{V}{3 k_B T^2} \int_0^\infty \langle \mathbf{J}(0) \cdot  \mathbf{J}(t)  \rangle \, \mathrm{d}t
 
 
 ----------
@@ -109,9 +129,15 @@ result should be: average conductivity ~0.29 in W/mK.
 
 **Output info:**
 
-This compute calculates a global vector of length 6 (total heat flux
-vector, followed by convective heat flux vector), which can be
-accessed by indices 1-6.  These values can be used by any command that
+This compute calculates a global vector of length 6.
+The first 3 components are the :math:`x`, :math:`y`, :math:`z`
+components of the full heat flux vector,
+i.e. (:math:`J_x`, :math:`J_y`, :math:`J_z`).
+The next 3 components are the :math:`x`, :math:`y`, :math:`z` components
+of just the convective portion of the flux, i.e. the
+first term in the equation for :math:`\mathbf{J}`.
+Each component can be
+accessed by indices 1-6. These values can be used by any command that
 uses global vector values from a compute as input.  See the :doc:`Howto output <Howto_output>` doc page for an overview of LAMMPS output
 options.
 
@@ -212,6 +238,22 @@ Related commands
    print        "average conductivity: $k[W/mK] @ $T K, ${ndens} /A\^3"
 
 
+----------
+
+
+.. _Surblys2:
+
+
+
+**(Surblys)** Surblys, Matsubara, Kikugawa, Ohara, Phys Rev E, 99, 051301(R) (2019).
+
+.. _Boone:
+
+
+
+**(Boone)** Boone, Babaei, Wilmer, J Chem Theory Comput, 15, 5579--5587 (2019).
+
+
 .. _lws: http://lammps.sandia.gov
 .. _ld: Manual.html
 .. _lc: Commands_all.html