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.
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. +
+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"
+
+(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.
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"