"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c

:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)

:line

compute heat/flux command :h3

[Syntax:]

compute ID group-ID heat/flux ke-ID pe-ID stress-ID :pre

ID, group-ID are documented in "compute"_compute.html command
heat/flux = style name of this compute command
ke-ID = ID of a compute that calculates per-atom kinetic energy
pe-ID = ID of a compute that calculates per-atom potential energy
stress-ID = ID of a compute that calculates per-atom stress :ul

[Examples:]

compute myFlux all heat/flux myKE myPE myStress :pre

[Description:]

Define a computation that calculates the heat flux vector based on
contributions from atoms in the specified group.  This can be used by
itself to measure the heat flux into or out of a reservoir of atoms,
or to calculate a thermal conductivity using the Green-Kubo formalism.

See the "fix thermal/conductivity"_fix_thermal_conductivity.html
command for details on how to compute thermal conductivity in an
alternate way, via the Muller-Plathe method.  See the "fix
heat"_fix_heat.html command for a way to control the heat added or
subtracted to a group of atoms.

The compute takes three arguments which are IDs of other
"computes"_compute.html.  One calculates per-atom kinetic energy
({ke-ID}), one calculates per-atom potential energy ({pe-ID)}, and the
third calculates per-atom stress ({stress-ID}).

NOTE: These other computes should provide values for all the atoms in
the group this compute specifies.  That means the other computes could
use the same group as this compute, or they can just use group "all"
(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:

:c,image(Eqs/heat_flux_J.jpg)

:c,image(Eqs/heat_flux_k.jpg)

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
included in the calculation.

NOTE: The "compute pe/atom"_compute_pe_atom.html and "compute
stress/atom"_compute_stress_atom.html 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 "compute stress/atom virial"_compute_stress_atom.html
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.

:line

The heat flux can be output every so many timesteps (e.g. via the
"thermo_style custom"_thermo_style.html command).  Then as a
post-processing operation, an autocorrelation can be performed, its
integral estimated, and the Green-Kubo formula above evaluated.

The "fix ave/correlate"_fix_ave_correlate.html command can calculate
the autocorrelation.  The trap() function in the
"variable"_variable.html command can calculate the integral.

An example LAMMPS input script for solid Ar is appended below.  The
result should be: average conductivity ~0.29 in W/mK.

:line

[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
uses global vector values from a compute as input.  See "this
section"_Section_howto.html#howto_15 for an overview of LAMMPS output
options.

The vector values calculated by this compute are "extensive", meaning
they scale with the number of atoms in the simulation.  They can be
divided by the appropriate volume to get a flux, which would then be
an "intensive" value, meaning independent of the number of atoms in
the simulation.  Note that if the compute is "all", then the
appropriate volume to divide by is the simulation box volume.
However, if a sub-group is used, it should be the volume containing
those atoms.

The vector values will be in energy*velocity "units"_units.html.  Once
divided by a volume the units will be that of flux, namely
energy/area/time "units"_units.html

[Restrictions:] none

[Related commands:]

"fix thermal/conductivity"_fix_thermal_conductivity.html,
"fix ave/correlate"_fix_ave_correlate.html,
"variable"_variable.html

[Default:] none

:line

# Sample LAMMPS input script for thermal conductivity of solid Ar :pre

units       real
variable    T equal 70
variable    V equal vol
variable    dt equal 4.0
variable    p equal 200     # correlation length
variable    s equal 10      # sample interval
variable    d equal $p*$s   # dump interval :pre

# convert from LAMMPS real units to SI :pre

variable    kB equal 1.3806504e-23    # \[J/K\] Boltzmann
variable    kCal2J equal 4186.0/6.02214e23
variable    A2m equal 1.0e-10
variable    fs2s equal 1.0e-15
variable    convert equal $\{kCal2J\}*$\{kCal2J\}/$\{fs2s\}/$\{A2m\} :pre

# setup problem :pre

dimension    3
boundary     p p p
lattice      fcc 5.376 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
region       box block 0 4 0 4 0 4
create_box   1 box
create_atoms 1 box
mass         1 39.948
pair_style   lj/cut 13.0
pair_coeff   * * 0.2381 3.405
timestep     $\{dt\}
thermo       $d :pre

# equilibration and thermalization :pre

velocity     all create $T 102486 mom yes rot yes dist gaussian
fix          NVT all nvt temp $T $T 10 drag 0.2
run          8000 :pre

# thermal conductivity calculation, switch to NVE if desired :pre

#unfix       NVT
#fix         NVE all nve :pre

reset_timestep 0
compute      myKE all ke/atom
compute      myPE all pe/atom
compute      myStress all stress/atom NULL virial
compute      flux all heat/flux myKE myPE myStress
variable     Jx equal c_flux\[1\]/vol
variable     Jy equal c_flux\[2\]/vol
variable     Jz equal c_flux\[3\]/vol
fix          JJ all ave/correlate $s $p $d &
             c_flux\[1\] c_flux\[2\] c_flux\[3\] type auto file J0Jt.dat ave running
variable     scale equal $\{convert\}/$\{kB\}/$T/$T/$V*$s*$\{dt\}
variable     k11 equal trap(f_JJ\[3\])*$\{scale\}
variable     k22 equal trap(f_JJ\[4\])*$\{scale\}
variable     k33 equal trap(f_JJ\[5\])*$\{scale\}
thermo_style custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33
run          100000
variable     k equal (v_k11+v_k22+v_k33)/3.0
variable     ndens equal count(all)/vol
print        "average conductivity: $k\[W/mK\] @ $T K, $\{ndens\} /A^3" :pre
