From d4769d99a324231e1c60262d1138c2884bd270c4 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Mon, 3 Jan 2011 16:42:55 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5439 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- doc/Manual.html | 8 ++ doc/Manual.txt | 6 +- doc/Section_howto.html | 187 ++++++++++++++++++++++++++++++++++++- doc/Section_howto.txt | 186 +++++++++++++++++++++++++++++++++++- doc/compute_heat_flux.html | 86 +++++++++-------- doc/compute_heat_flux.txt | 76 ++++++++------- 6 files changed, 471 insertions(+), 78 deletions(-) diff --git a/doc/Manual.html b/doc/Manual.html index 1d118439ed..a795349fd3 100644 --- a/doc/Manual.html +++ b/doc/Manual.html @@ -159,6 +159,10 @@ listed above. 4.18 Elastic constants
4.19 Library interface to LAMMPS +
+ 4.20 Calculating thermal conductivity +
+ 4.21 Calculating viscosity
  • Example problems @@ -295,6 +299,10 @@ listed above. + + + + diff --git a/doc/Manual.txt b/doc/Manual.txt index f086010e9a..3ff50e844b 100644 --- a/doc/Manual.txt +++ b/doc/Manual.txt @@ -113,7 +113,9 @@ listed above. 4.16 "Thermostatting, barostatting, and compute temperature"_4_16 :b 4.17 "Walls"_4_17 :b 4.18 "Elastic constants"_4_18 :b - 4.19 "Library interface to LAMMPS"_4_19 :ule,b + 4.19 "Library interface to LAMMPS"_4_19 :b + 4.20 "Calculating thermal conductivity"_4_20 :b + 4.21 "Calculating viscosity"_4_21 :ule,b "Example problems"_Section_example.html :l "Performance & scalability"_Section_perf.html :l "Additional tools"_Section_tools.html :l @@ -176,6 +178,8 @@ listed above. :link(4_17,Section_howto.html#4_17) :link(4_18,Section_howto.html#4_18) :link(4_19,Section_howto.html#4_19) +:link(4_20,Section_howto.html#4_20) +:link(4_21,Section_howto.html#4_21) :link(9_1,Section_python.html#9_1) :link(9_2,Section_python.html#9_2) diff --git a/doc/Section_howto.html b/doc/Section_howto.html index 4d729e8f0a..91ce3d9730 100644 --- a/doc/Section_howto.html +++ b/doc/Section_howto.html @@ -32,7 +32,9 @@ LAMMPS. 4.16 Thermostatting, barostatting and computing temperature
    4.17 Walls
    4.18 Elastic constants
    -4.19 Library interface to LAMMPS
    +4.19 Library interface to LAMMPS
    +4.20 Calculating thermal conductivity
    +4.21 Calculating viscosity

    The example input scripts included in the LAMMPS distribution and highlighted in this section also show how to @@ -1707,6 +1709,184 @@ grab data from LAMMPS, change it, and put it back into LAMMPS.

    +

    4.20 Calculating thermal conductivity +

    +

    The thermal conductivity kappa of a material can be measured in at +least 3 ways using various options in LAMMPS. (See this +section of the manual for an analogous +discussion for viscosity). The thermal conducitivity tensor kappa is +a measure of the propensity of a material to transmit heat energy in a +diffusive manner as given by Fourier's law +

    +

    J = -kappa grad(T) +

    +

    where J is the heat flux in units of energy per area per time and +grad(T) is the spatial gradient of temperature. The thermal +conductivity thus has units of energy per distance per time per degree +K and is often approximated as an isotropic quantity, i.e. as a +scalar. +

    +

    The first method is to setup two thermostatted regions at opposite +ends of a simulation box, or one in the middle and one at the end of a +periodic box. By holding the two regions at different temperatures +with a thermostatting fix, the energy added +to the hot region should equal the energy subtracted from the cold +region and be proportional to the heat flux moving between the +regions. See the paper by Ikeshoji and Hafskjold for +details of this idea. Note that thermostatting fixes such as fix +nvt, fix langevin, and fix +temp/rescale store the cumulative energy they +add/subtract. Alternatively, the fix heat command can +used in place of thermostats on each of two regions, and the resulting +temperatures of the two regions monitored with the "compute +temp/region" command or the temperature profile of the intermediate +region monitored with the fix ave/spatial and +compute ke/atom commands. +

    +

    The second method is to perform a reverse non-equilibrium MD +simulation using the fix +thermal/conductivity command which +implements the rNEMD algorithm of Muller-Plathe. Kinetic energy is +swapped between atoms in two different layers of the simulation box. +This induces a temperature gradient between the two layers which can +be monitored with the fix ave/spatial and +compute ke/atom commands. The fix tallies the +cumulative energy transfer that it performs. See the fix +thermal/conductivity command for +details. +

    +

    The third method is based on the Green-Kubo (GK) formula which relates +the ensemble average of the auto-correlation of the heat flux to +kappa. The heat flux can be calculated from the fluctuations of +per-atom potential and kinetic energies and per-atom stress tensor in +a steady-state equilibrated simulation. This is in contrast to the +two preceding non-equilibrium methods, where energy flows continuously +between hot and cold regions of the simulation box. +

    +

    The compute heat/flux command can calculate +the needed heat flux and describes how to implement the Green_Kubo +formalism using additional LAMMPS commands, such as the fix +ave/correlate command to calculate the needed +auto-correlation. See the doc page for the compute +heat/flux command for an example input script +that calculates the thermal conductivity of solid Ar via the GK +formalism. +

    +
    + +

    4.21 Calculating viscosity +

    +

    The shear viscosity eta of a fluid can be measured in at least 3 ways +using various options in LAMMPS. (See this +section??? of the manual for an analogous +discussion for thermal conductivity). Eta is a measure of the +propensity of a fluid to transmit momentum in a direction +perpendicular to the direction of velocity or momentum flow. +Alternatively it is the resistance the fluid has to being sheared. It +is given by +

    +

    J = -eta grad(Vstream) +

    +

    where J is the momentum flux in units of momentum per area per time. +and grad(Vstream) is the spatial gradient of the velocity of the fluid +moving in another direction, normal to the area through which the +momentum flows. Viscosity thus has units of pressure-time. +

    +

    The first method is to perform a non-equlibrium MD (NEMD) simulation +by shearing the simulation box via the fix deform +command, and using the fix nvt/sllod command to +thermostat the fluid via the SLLOD equations of motion. The velocity +profile setup in the fluid by this procedure can be monitored by the +fix ave/spatial command, which determines +grad(Vstream) in the equation above. E.g. the derivative in the +y-direction of the Vx component of fluid motion or grad(Vstream) = +dVx/dy. In this case, the Pxy off-diagonal component of the pressure +or stress tensor, as calculated by the compute +pressure command, can also be monitored, which +is the J term in the equation above. See this +section of the manual for details on NEMD +simulations. +

    +

    The second method is to perform a reverse non-equilibrium MD +simulation using the fix viscosity command which +implements the rNEMD algorithm of Muller-Plathe. Momentum in one +dimension is swapped between atoms in two different layers of the +simulation box in a different dimension. This induces a velocity +gradient which can be monitored with the fix +ave/spatial command. The fix tallies the +cummulative momentum transfer that it performs. See the fix +viscosity command for details. +

    +

    The third method is based on the Green-Kubo (GK) formula which relates +the ensemble average of the auto-correlation of the stress/pressure +tensor to eta. This can be done in a steady-state equilibrated +simulation which is in contrast to the two preceding non-equilibrium +methods, where momentum flows continuously through the simulation box. +

    +

    Here is an example input script that calculates the viscosity of +liquid Ar via the GK formalism: +

    +
    # Sample LAMMPS input script for viscosity of liquid Ar 
+
    +
    units       real
+variable    T equal 86.4956
+variable    V equal vol
+variable    dt equal 4.0
+variable    p equal 400     # correlation length
+variable    s equal 5       # sample interval
+variable    d equal $p*$s   # dump interval 
+
    +
    # convert from LAMMPS real units to SI 
+
    +
    variable    kB equal 1.3806504e-23    # [J/K/ Boltzmann
+variable    atm2Pa equal 101325.0
+variable    A2m equal 1.0e-10
+variable    fs2s equal 1.0e-15
+variable    convert equal ${atm2Pa}*${atm2Pa}*${fs2s}*${A2m}*${A2m}*${A2m} 
+
    +
    # setup problem 
+
    +
    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 
+
    +
    # equilibration and thermalization 
+
    +
    velocity     all create $T 102486 mom yes rot yes dist gaussian
+fix          NVT all nvt temp $T $T 10 drag 0.2
+run          8000 
+
    +
    # viscosity calculation, switch to NVE if desired 
+
    +
    #unfix       NVT
+#fix         NVE all nve 
+
    +
    reset_timestep 0
+variable     pxy equal pxy
+variable     pxz equal pxz
+variable     pyz equal pyz
+fix          SS all ave/correlate $s $p $d &
+             v_pxy v_pxz v_pyz type auto file S0St.dat ave running
+variable     scale equal ${convert}/(${kB}*$T)*$V*$s*${dt}
+variable     v11 equal trap(f_SS[3/)*${scale}
+variable     v22 equal trap(f_SS[4/)*${scale}
+variable     v33 equal trap(f_SS[5/)*${scale}
+thermo_style custom step temp press v_pxy v_pxz v_pyz v_v11 v_v22 v_v33
+run          100000
+variable     v equal (v_v11+v_v22+v_v33)/3.0
+variable     ndens equal count(all)/vol
+print        "average viscosity: $v [Pa.s/ @ $T K, ${ndens} /A^3" 
+
    +
    +
    @@ -1724,6 +1904,11 @@ Spellmeyer, Fox, Caldwell, Kollman, JACS 117, 5179-5197 (1995).

    (Horn) Horn, Swope, Pitera, Madura, Dick, Hura, and Head-Gordon, J Chem Phys, 120, 9665 (2004).

    + + +

    (Ikeshoji) Ikeshoji and Hafskjold, Molecular Physics, 81, 251-261 +(1994). +

    (MacKerell) MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field, diff --git a/doc/Section_howto.txt b/doc/Section_howto.txt index baf09fc10d..832ef25ae5 100644 --- a/doc/Section_howto.txt +++ b/doc/Section_howto.txt @@ -29,7 +29,9 @@ LAMMPS. 4.16 "Thermostatting, barostatting and computing temperature"_#4_16 4.17 "Walls"_#4_17 4.18 "Elastic constants"_#4_18 -4.19 "Library interface to LAMMPS"_#4_19 :all(b) +4.19 "Library interface to LAMMPS"_#4_19 +4.20 "Calculating thermal conductivity"_#4_20 +4.21 "Calculating viscosity"_#4_21 :all(b) The example input scripts included in the LAMMPS distribution and highlighted in "this section"_Section_example.html also show how to @@ -1693,6 +1695,184 @@ have example C++ and C and Python codes which show how a driver code can link to LAMMPS as a library, run LAMMPS on a subset of processors, grab data from LAMMPS, change it, and put it back into LAMMPS. +:line + +4.20 Calculating thermal conductivity :link(4_20),h4 + +The thermal conductivity kappa of a material can be measured in at +least 3 ways using various options in LAMMPS. (See "this +section"_Section_howto.html#4_21 of the manual for an analogous +discussion for viscosity). The thermal conducitivity tensor kappa is +a measure of the propensity of a material to transmit heat energy in a +diffusive manner as given by Fourier's law + +J = -kappa grad(T) + +where J is the heat flux in units of energy per area per time and +grad(T) is the spatial gradient of temperature. The thermal +conductivity thus has units of energy per distance per time per degree +K and is often approximated as an isotropic quantity, i.e. as a +scalar. + +The first method is to setup two thermostatted regions at opposite +ends of a simulation box, or one in the middle and one at the end of a +periodic box. By holding the two regions at different temperatures +with a "thermostatting fix"_Section_howto.html#4_13, the energy added +to the hot region should equal the energy subtracted from the cold +region and be proportional to the heat flux moving between the +regions. See the paper by "Ikeshoji and Hafskjold"_#Ikeshoji for +details of this idea. Note that thermostatting fixes such as "fix +nvt"_fix_nh.html, "fix langevin"_fix_langevin.html, and "fix +temp/rescale"_fix_temp_rescale.html store the cumulative energy they +add/subtract. Alternatively, the "fix heat"_fix_heat.html command can +used in place of thermostats on each of two regions, and the resulting +temperatures of the two regions monitored with the "compute +temp/region" command or the temperature profile of the intermediate +region monitored with the "fix ave/spatial"_fix_ave_spatial.html and +"compute ke/atom"_compute_ke_atom.html commands. + +The second method is to perform a reverse non-equilibrium MD +simulation using the "fix +thermal/conductivity"_fix_thermal_conductivity.html command which +implements the rNEMD algorithm of Muller-Plathe. Kinetic energy is +swapped between atoms in two different layers of the simulation box. +This induces a temperature gradient between the two layers which can +be monitored with the "fix ave/spatial"_fix_ave_spatial.html and +"compute ke/atom"_compute_ke_atom.html commands. The fix tallies the +cumulative energy transfer that it performs. See the "fix +thermal/conductivity"_fix_thermal_conductivity.html command for +details. + +The third method is based on the Green-Kubo (GK) formula which relates +the ensemble average of the auto-correlation of the heat flux to +kappa. The heat flux can be calculated from the fluctuations of +per-atom potential and kinetic energies and per-atom stress tensor in +a steady-state equilibrated simulation. This is in contrast to the +two preceding non-equilibrium methods, where energy flows continuously +between hot and cold regions of the simulation box. + +The "compute heat/flux"_compute_heat_flux.html command can calculate +the needed heat flux and describes how to implement the Green_Kubo +formalism using additional LAMMPS commands, such as the "fix +ave/correlate"_fix_ave_correlate.html command to calculate the needed +auto-correlation. See the doc page for the "compute +heat/flux"_compute_heat_flux.html command for an example input script +that calculates the thermal conductivity of solid Ar via the GK +formalism. + +:line + +4.21 Calculating viscosity :link(4_21),h4 + +The shear viscosity eta of a fluid can be measured in at least 3 ways +using various options in LAMMPS. (See "this +section"_Section_howto.html#4_??? of the manual for an analogous +discussion for thermal conductivity). Eta is a measure of the +propensity of a fluid to transmit momentum in a direction +perpendicular to the direction of velocity or momentum flow. +Alternatively it is the resistance the fluid has to being sheared. It +is given by + +J = -eta grad(Vstream) + +where J is the momentum flux in units of momentum per area per time. +and grad(Vstream) is the spatial gradient of the velocity of the fluid +moving in another direction, normal to the area through which the +momentum flows. Viscosity thus has units of pressure-time. + +The first method is to perform a non-equlibrium MD (NEMD) simulation +by shearing the simulation box via the "fix deform"_fix_deform.html +command, and using the "fix nvt/sllod"_fix_nvt_sllod.html command to +thermostat the fluid via the SLLOD equations of motion. The velocity +profile setup in the fluid by this procedure can be monitored by the +"fix ave/spatial"_fix_ave_spatial.html command, which determines +grad(Vstream) in the equation above. E.g. the derivative in the +y-direction of the Vx component of fluid motion or grad(Vstream) = +dVx/dy. In this case, the Pxy off-diagonal component of the pressure +or stress tensor, as calculated by the "compute +pressure"_compute_pressure.html command, can also be monitored, which +is the J term in the equation above. See "this +section"_Section_howto.html#4_13 of the manual for details on NEMD +simulations. + +The second method is to perform a reverse non-equilibrium MD +simulation using the "fix viscosity"_fix_viscosity.html command which +implements the rNEMD algorithm of Muller-Plathe. Momentum in one +dimension is swapped between atoms in two different layers of the +simulation box in a different dimension. This induces a velocity +gradient which can be monitored with the "fix +ave/spatial"_fix_ave_spatial.html command. The fix tallies the +cummulative momentum transfer that it performs. See the "fix +viscosity"_fix_viscosity.html command for details. + +The third method is based on the Green-Kubo (GK) formula which relates +the ensemble average of the auto-correlation of the stress/pressure +tensor to eta. This can be done in a steady-state equilibrated +simulation which is in contrast to the two preceding non-equilibrium +methods, where momentum flows continuously through the simulation box. + +Here is an example input script that calculates the viscosity of +liquid Ar via the GK formalism: + +# Sample LAMMPS input script for viscosity of liquid Ar :pre + +units real +variable T equal 86.4956 +variable V equal vol +variable dt equal 4.0 +variable p equal 400 # correlation length +variable s equal 5 # 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 atm2Pa equal 101325.0 +variable A2m equal 1.0e-10 +variable fs2s equal 1.0e-15 +variable convert equal $\{atm2Pa\}*$\{atm2Pa\}*$\{fs2s\}*$\{A2m\}*$\{A2m\}*$\{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 + +# viscosity calculation, switch to NVE if desired :pre + +#unfix NVT +#fix NVE all nve :pre + +reset_timestep 0 +variable pxy equal pxy +variable pxz equal pxz +variable pyz equal pyz +fix SS all ave/correlate $s $p $d & + v_pxy v_pxz v_pyz type auto file S0St.dat ave running +variable scale equal $\{convert\}/($\{kB\}*$T)*$V*$s*$\{dt\} +variable v11 equal trap(f_SS\[3/])*$\{scale\} +variable v22 equal trap(f_SS\[4/])*$\{scale\} +variable v33 equal trap(f_SS\[5/])*$\{scale\} +thermo_style custom step temp press v_pxy v_pxz v_pyz v_v11 v_v22 v_v33 +run 100000 +variable v equal (v_v11+v_v22+v_v33)/3.0 +variable ndens equal count(all)/vol +print "average viscosity: $v \[Pa.s/] @ $T K, $\{ndens\} /A^3" :pre + :line :line @@ -1708,6 +1888,10 @@ Spellmeyer, Fox, Caldwell, Kollman, JACS 117, 5179-5197 (1995). [(Horn)] Horn, Swope, Pitera, Madura, Dick, Hura, and Head-Gordon, J Chem Phys, 120, 9665 (2004). +:link(Ikeshoji) +[(Ikeshoji)] Ikeshoji and Hafskjold, Molecular Physics, 81, 251-261 +(1994). + :link(MacKerell) [(MacKerell)] MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field, Fischer, Gao, Guo, Ha, et al, J Phys Chem, 102, 3586 (1998). diff --git a/doc/compute_heat_flux.html b/doc/compute_heat_flux.html index d930637f3c..84f8b28351 100644 --- a/doc/compute_heat_flux.html +++ b/doc/compute_heat_flux.html @@ -89,8 +89,8 @@ integral estimated, and the Green-Kubo formula above evaluated. the autocorrelation. The trap() function in the variable command can calculate the integral.

    -

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

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

    @@ -128,24 +128,29 @@ energy/area/time units

    -

    Sample LAMMPS input script -

    -
    atom_style  atomic
-units       real
-variable    kB equal 1.3806504e-23 # [J/K] Boltzmann
-variable    kCal2J equal 4186.0/6.02214e23
+
    # Sample LAMMPS input script for thermal conductivity of solid Ar 
+
    
+
    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 
+variable    p equal 200     # correlation length
+variable    s equal 10      # sample interval
+variable    d equal $p*$s   # dump interval 
+
    
+
    # convert from LAMMPS real units to SI 
+
    
+
    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} 
+
    
+
    # setup problem 
 
    
-
    # ---------------------------------------------------------
-
    
 
    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
+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
@@ -155,34 +160,35 @@ pair_coeff   * * 0.2381 3.405
 timestep     ${dt}
 thermo	     $d 
 
    
-
    # ------------- equilibration and thermalization ----------------
-
    
-
    velocity  all create $T 102486 mom yes rot yes dist gaussian
-fix       NVT all nvt temp $T $T 10 drag 0.2
-run       8000 
+
    # equilibration and thermalization 
 
    
-
    # -------------- flux calculation, switch to NVE if desired --------
-
    
-
    #unfix     NVT
-#fix       NVE all nve 
+
    velocity     all create $T 102486 mom yes rot yes dist gaussian
+fix          NVT all nvt temp $T $T 10 drag 0.2
+run          8000 
 
    
-
    reset_timestep  0
-compute   myKE all ke/atom
-compute   myPE all pe/atom
-compute   myStress all stress/atom 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 ${kCal2J}*${kCal2J}/${kB}/$T/$T/$V*$s*${dt}*1.0e25
-variable  k11 equal trap(f_JJ[3])*${scale}
-variable  k22 equal trap(f_JJ[4])*${scale}
-variable  k33 equal trap(f_JJ[5])*${scale}
+
    # thermal conductivity calculation, switch to NVE if desired 
+
    
+
    #unfix       NVT
+#fix         NVE all nve 
+
    
+
    reset_timestep 0
+compute      myKE all ke/atom
+compute      myPE all pe/atom
+compute      myStress all stress/atom 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
-print     "average conductivity: $k [W/mK]" 
+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" 
 
    
 
diff --git a/doc/compute_heat_flux.txt b/doc/compute_heat_flux.txt
index 3826b24a42..ef4907cdd3 100644
--- a/doc/compute_heat_flux.txt
+++ b/doc/compute_heat_flux.txt
@@ -86,8 +86,8 @@ The "fix ave/correlate"_fix_ave_correlate.html command can calclate
 the autocorrelation.  The trap() function in the
 "variable"_variable.html command can calculate the integral.
 
-An an example LAMMPS input script for solid Ar is appended below.  
-The result should be: average conductivity ~0.29 in W/mK.
+An example LAMMPS input script for solid Ar is appended below.  The
+result should be: average conductivity ~0.29 in W/mK.
 
 :line
 
@@ -125,24 +125,29 @@ energy/area/time "units"_units.html
 
 :line
 
-Sample LAMMPS input script :h4
+# Sample LAMMPS input script for thermal conductivity of solid Ar :pre
 
-atom_style  atomic
 units       real
-variable    kB equal 1.3806504e-23 # \[J/K\] Boltzmann
-variable    kCal2J equal 4186.0/6.02214e23
 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
+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
+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
@@ -152,32 +157,33 @@ pair_coeff   * * 0.2381 3.405
 timestep     $\{dt\}
 thermo	     $d :pre
 
-# ------------- equilibration and thermalization ----------------
+# 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
+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
 
-# -------------- flux calculation, switch to NVE if desired --------
+# thermal conductivity calculation, switch to NVE if desired :pre
 
-#unfix     NVT
-#fix       NVE all nve :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 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 $\{kCal2J\}*$\{kCal2J\}/$\{kB\}/$T/$T/$V*$s*$\{dt\}*1.0e25
-variable  k11 equal trap(f_JJ\[3\])*$\{scale\}
-variable  k22 equal trap(f_JJ\[4\])*$\{scale\}
-variable  k33 equal trap(f_JJ\[5\])*$\{scale\}
+reset_timestep 0
+compute      myKE all ke/atom
+compute      myPE all pe/atom
+compute      myStress all stress/atom 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
-print     "average conductivity: $k \[W/mK\]" :pre
+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