Commit d4769d99 authored by sjplimp's avatar sjplimp
Browse files

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5439 f3b2605a-c512-4ea7-a41b-209d697bcdaa
parent c49be34d
Loading
Loading
Loading
Loading
+8 −0
Original line number Diff line number Diff line
@@ -159,6 +159,10 @@ listed above.
  4.18 <A HREF = "Section_howto.html#4_18">Elastic constants</A> 
<BR>
  4.19 <A HREF = "Section_howto.html#4_19">Library interface to LAMMPS</A> 
<BR>
  4.20 <A HREF = "Section_howto.html#4_20">Calculating thermal conductivity</A> 
<BR>
  4.21 <A HREF = "Section_howto.html#4_21">Calculating viscosity</A> 
<BR></UL>
<LI><A HREF = "Section_example.html">Example problems</A> 

@@ -295,6 +299,10 @@ listed above.










+5 −1
Original line number Diff line number Diff line
@@ -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)
+186 −1
Original line number Diff line number Diff line
@@ -33,6 +33,8 @@ LAMMPS.
4.17 <A HREF = "#4_17">Walls</A><BR>
4.18 <A HREF = "#4_18">Elastic constants</A><BR>
4.19 <A HREF = "#4_19">Library interface to LAMMPS</A><BR>
4.20 <A HREF = "#4_20">Calculating thermal conductivity</A><BR>
4.21 <A HREF = "#4_21">Calculating viscosity</A> <BR>

<P>The example input scripts included in the LAMMPS distribution and
highlighted in <A HREF = "Section_example.html">this section</A> also show how to
@@ -1707,6 +1709,184 @@ grab data from LAMMPS, change it, and put it back into LAMMPS.
</P>
<HR>

<A NAME = "4_20"></A><H4>4.20 Calculating thermal conductivity 
</H4>
<P>The thermal conductivity kappa of a material can be measured in at
least 3 ways using various options in LAMMPS.  (See <A HREF = "Section_howto.html#4_21">this
section</A> 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
</P>
<P>J = -kappa grad(T)
</P>
<P>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.
</P>
<P>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 <A HREF = "Section_howto.html#4_13">thermostatting fix</A>, 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 <A HREF = "#Ikeshoji">Ikeshoji and Hafskjold</A> for
details of this idea.  Note that thermostatting fixes such as <A HREF = "fix_nh.html">fix
nvt</A>, <A HREF = "fix_langevin.html">fix langevin</A>, and <A HREF = "fix_temp_rescale.html">fix
temp/rescale</A> store the cumulative energy they
add/subtract.  Alternatively, the <A HREF = "fix_heat.html">fix heat</A> 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 <A HREF = "fix_ave_spatial.html">fix ave/spatial</A> and
<A HREF = "compute_ke_atom.html">compute ke/atom</A> commands.
</P>
<P>The second method is to perform a reverse non-equilibrium MD
simulation using the <A HREF = "fix_thermal_conductivity.html">fix
thermal/conductivity</A> 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 <A HREF = "fix_ave_spatial.html">fix ave/spatial</A> and
<A HREF = "compute_ke_atom.html">compute ke/atom</A> commands.  The fix tallies the
cumulative energy transfer that it performs.  See the <A HREF = "fix_thermal_conductivity.html">fix
thermal/conductivity</A> command for
details.
</P>
<P>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.
</P>
<P>The <A HREF = "compute_heat_flux.html">compute heat/flux</A> command can calculate
the needed heat flux and describes how to implement the Green_Kubo
formalism using additional LAMMPS commands, such as the <A HREF = "fix_ave_correlate.html">fix
ave/correlate</A> command to calculate the needed
auto-correlation.  See the doc page for the <A HREF = "compute_heat_flux.html">compute
heat/flux</A> command for an example input script
that calculates the thermal conductivity of solid Ar via the GK
formalism.
</P>
<HR>

<A NAME = "4_21"></A><H4>4.21 Calculating viscosity 
</H4>
<P>The shear viscosity eta of a fluid can be measured in at least 3 ways
using various options in LAMMPS.  (See <A HREF = "Section_howto.html#4_">this
section</A>??? 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
</P>
<P>J = -eta grad(Vstream)
</P>
<P>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.
</P>
<P>The first method is to perform a non-equlibrium MD (NEMD) simulation
by shearing the simulation box via the <A HREF = "fix_deform.html">fix deform</A>
command, and using the <A HREF = "fix_nvt_sllod.html">fix nvt/sllod</A> 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
<A HREF = "fix_ave_spatial.html">fix ave/spatial</A> 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 <A HREF = "compute_pressure.html">compute
pressure</A> command, can also be monitored, which
is the J term in the equation above.  See <A HREF = "Section_howto.html#4_13">this
section</A> of the manual for details on NEMD
simulations.
</P>
<P>The second method is to perform a reverse non-equilibrium MD
simulation using the <A HREF = "fix_viscosity.html">fix viscosity</A> 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 <A HREF = "fix_ave_spatial.html">fix
ave/spatial</A> command.  The fix tallies the
cummulative momentum transfer that it performs.  See the <A HREF = "fix_viscosity.html">fix
viscosity</A> command for details.
</P>
<P>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.
</P>
<P>Here is an example input script that calculates the viscosity of
liquid Ar via the GK formalism:
</P>
<PRE># Sample LAMMPS input script for viscosity of liquid Ar 
</PRE>
<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>
<PRE># convert from LAMMPS real units to SI 
</PRE>
<PRE>variable    kB equal 1.3806504e-23    # [J/K/</B> 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>
<PRE># setup problem 
</PRE>
<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>
<PRE># equilibration and thermalization 
</PRE>
<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>
<PRE># viscosity calculation, switch to NVE if desired 
</PRE>
<PRE>#unfix       NVT
#fix         NVE all nve 
</PRE>
<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/</B>)*${scale}
variable     v22 equal trap(f_SS[4/</B>)*${scale}
variable     v33 equal trap(f_SS[5/</B>)*${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/</B> @ $T K, ${ndens} /A^3" 
</PRE>
<HR>

<HR>

<A NAME = "Berendsen"></A>
@@ -1724,6 +1904,11 @@ Spellmeyer, Fox, Caldwell, Kollman, JACS 117, 5179-5197 (1995).
<P><B>(Horn)</B> Horn, Swope, Pitera, Madura, Dick, Hura, and Head-Gordon,
J Chem Phys, 120, 9665 (2004).
</P>
<A NAME = "Ikeshoji"></A>

<P><B>(Ikeshoji)</B> Ikeshoji and Hafskjold, Molecular Physics, 81, 251-261
(1994).
</P>
<A NAME = "MacKerell"></A>

<P><B>(MacKerell)</B> MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field,
+185 −1
Original line number Diff line number Diff line
@@ -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).
+46 −40
Original line number Diff line number Diff line
@@ -89,8 +89,8 @@ integral estimated, and the Green-Kubo formula above evaluated.
the autocorrelation.  The trap() function in the
<A HREF = "variable.html">variable</A> command can calculate the integral.
</P>
<P>An an example LAMMPS input script for solid Ar is appended below.  
The result should be: average conductivity ~0.29 in W/mK.
<P>An example LAMMPS input script for solid Ar is appended below.  The
result should be: average conductivity ~0.29 in W/mK.
</P>
<HR>

@@ -128,12 +128,9 @@ energy/area/time <A HREF = "units.html">units</A>
</P>
<HR>

<H4>Sample LAMMPS input script 
</H4>
<PRE>atom_style  atomic
units       real
variable    kB equal 1.3806504e-23 # [J/K] Boltzmann
variable    kCal2J equal 4186.0/6.02214e23
<PRE># Sample LAMMPS input script for thermal conductivity of solid Ar 
</PRE>
<PRE>units       real
variable    T equal 70
variable    V equal vol
variable    dt equal 4.0
@@ -141,8 +138,16 @@ variable p equal 200 # correlation length
variable    s equal 10      # sample interval
variable    d equal $p*$s   # dump interval 
</PRE>
<P># ---------------------------------------------------------
</P>
<PRE># convert from LAMMPS real units to SI 
</PRE>
<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>
<PRE># setup problem 
</PRE>
<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
@@ -155,14 +160,14 @@ pair_coeff * * 0.2381 3.405
timestep     ${dt}
thermo	     $d 
</PRE>
<P># ------------- equilibration and thermalization ----------------
</P>
<PRE># equilibration and thermalization 
</PRE>
<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>
<P># -------------- flux calculation, switch to NVE if desired --------
</P>
<PRE># thermal conductivity calculation, switch to NVE if desired 
</PRE>
<PRE>#unfix       NVT
#fix         NVE all nve 
</PRE>
@@ -176,13 +181,14 @@ 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     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]" 
variable     ndens equal count(all)/vol
print        "average conductivity: $k[W/mK] @ $T K, ${ndens} /A^3" 
</PRE>
</HTML>
Loading