update compute heat/flux documentation

   (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:

.. image:: Eqs/heat_flux_J.jpg

   :align: center

.. 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

In case of pairwise interactions the heat flux is defined as:

.. 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]

: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

   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

----------

**Output info:**

This compute calculates a global vector of length 6 (total heat flux

vector, followed by convective heat flux vector), which can be

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.

   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