ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13683 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2015-07-22 16:06:41 +00:00
parent 6864dfa2d1
commit a5a4d41fe4
4 changed files with 827 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

232
doc/fix_qbmsst.html Normal file
View File

@ -0,0 +1,232 @@
<HTML>
<CENTER><A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A>
</CENTER>
<HR>
<H3>fix qbmsst command
</H3>
<P><B>Syntax:</B>
</P>
<PRE>fix ID group-ID qbmsst dir shockvel keyword value ...
</PRE>
<UL><LI>ID, group-ID are documented in <A HREF = "fix.html">fix</A> command
<LI>qbmsst = style name of this fix
<LI>dir = <I>x</I> or <I>y</I> or <I>z</I>
<LI>shockvel = shock velocity (strictly positive, velocity units)
<LI>zero or more keyword/value pairs may be appended
<LI>keyword = <I>q</I> or <I>mu</I> or <I>p0</I> or <I>v0</I> or <I>e0</I> or <I>tscale</I> or <I>damp</I> or <I>seed</I>or <I>f_max</I> or <I>N_f</I> or <I>eta</I> or <I>beta</I> or <I>T_init</I>
<PRE> <I>q</I> value = cell mass-like parameter (mass^2/distance^4 units)
<I>mu</I> value = artificial viscosity (mass/distance/time units)
<I>p0</I> value = initial pressure in the shock equations (pressure units)
<I>v0</I> value = initial simulation cell volume in the shock equations (distance^3 units)
<I>e0</I> value = initial total energy (energy units)
<I>tscale</I> value = reduction in initial temperature (unitless fraction between 0.0 and 1.0)
<I>damp</I> value = damping parameter (time units) inverse of friction <i>&gamma;</i>
<I>seed</I> value = random number seed (positive integer)
<I>f_max</I> value = upper cutoff frequency of the vibration spectrum (1/time units)
<I>N_f</I> value = number of frequency bins (positive integer)
<I>eta</I> value = coupling constant between the shock system and the quantum thermal bath (positive unitless)
<I>beta</I> value = the quantum temperature is updated every beta time steps (positive integer)
<I>T_init</I> value = quantum temperature for the initial state (temperature units)
</PRE>
</UL>
<P><B>Examples:</B>
</P>
<PRE>fix 1 all qbmsst z 0.122 q 25 mu 0.9 tscale 0.01 damp 200 seed 35082 f_max 0.3 N_f 100 eta 1 beta 400 T_init 110 (liquid methane modeled with the REAX force field, real units)
fix 2 all qbmsst z 72 q 40 tscale 0.05 damp 1 seed 47508 f_max 120.0 N_f 100 eta 1.0 beta 500 T_init 300 (quartz modeled with the BKS force field, metal units)
</PRE>
<P>Two example input scripts are given, including shocked alpha quartz
and shocked liquid methane. The input script first equilibrate an
initial state with the quantum thermal bath at the target temperature
and then apply the qbmsst to simulate shock compression with quantum
nuclear correction. The following two figures plot related quantities
for shocked alpha quartz.
</P>
<CENTER><IMG SRC = "JPG/qbmsst_init.jpg">
</CENTER>
<P>Figure 1. Classical temperature <i>T</i><sup>cl</sup> = &sum;
<i>m<sub>i</sub>v<sub>i</sub><sup>2</sup>/3Nk</i><sub>B</sub> vs. time
for coupling the alpha quartz initial state with the quantum thermal
bath at target quantum temperature <i>T</i><sup>qm</sup> = 300 K. The
NpH ensemble is used for time integration while QTB provides the
colored random force. <i>T</i><sup>cl</sup> converges at the timescale
of <I>damp</I> which is set to be 1 ps.
</P>
<CENTER><IMG SRC = "JPG/qbmsst_shock.jpg">
</CENTER>
<P>Figure 2. Quantum temperature and pressure vs. time for simulating
shocked alpha quartz with the QBMSST. The shock propagates along the z
direction. Restart of the QBMSST command is demonstrated in the
example input script. Thermodynamic quantities stay continuous before
and after the restart.
</P>
<P><B>Description:</B>
</P>
<P>This command performs the Quantum-Bath coupled Multi-Scale Shock
Technique (QBMSST) integration. See <A HREF = "#Qi">(Qi)</A> for a detailed
description of this method. The QBMSST provides description of the
thermodynamics and kinetics of shock processes while incorporating
quantum nuclear effects. The <I>shockvel</I> setting determines the steady
shock velocity that will be simulated along direction <I>dir</I>.
</P>
<P>Quantum nuclear effects <A HREF = "fix_qtb.html">(fix qtb)</A> can be crucial
especially when the temperature of the initial state is below the
classical limit or there is a great change in the zero point energies
between the initial and final states. Theoretical post processing
quantum corrections of shock compressed water and methane have been
reported as much as 30% of the temperatures <A HREF = "#Goldman">(Goldman)</A>. A
self-consistent method that couples the shock to a quantum thermal
bath described by a colored noise Langevin thermostat has been
developed by Qi et al <A HREF = "#Qi">(Qi)</A> and applied to shocked methane. The
onset of chemistry is reported to be at a pressure on the shock
Hugoniot that is 40% lower than observed with classical molecular
dynamics.
</P>
<P>It is highly recommended that the system be already in an equilibrium
state with a quantum thermal bath at temperature of <I>T_init</I>. The fix
command <A HREF = "fix_qtb.html">fix qtb</A> at constant temperature <I>T_init</I> could
be used before applying this command to introduce self-consistent
quantum nuclear effects into the initial state.
</P>
<P>The parameters <I>q</I>, <I>mu</I>, <I>e0</I>, <I>p0</I>, <I>v0</I> and <I>tscale</I> are described
in the command <A HREF = "fix_msst.html">fix msst</A>. The values of <I>e0</I>, <I>p0</I>, or
<I>v0</I> will be calculated on the first step if not specified. The
parameter of <I>damp</I>, <I>f_max</I>, and <I>N_f</I> are described in the command
<A HREF = "fix_qtb.html">fix qtb</A>.
</P>
<P>The fix qbmsst command couples the shock system to a quantum thermal
bath with a rate that is proportional to the change of the total
energy of the shock system, <i>etot</i> - <i>etot</i><sub>0</sub>.
Here <i>etot</i> consists of both the system energy and a thermal
term, see <A HREF = "#Qi">(Qi)</A>, and <i>etot</i><sub>0</sub> = <I>e0</I> is the
initial total energy.
</P>
<P>The <I>eta</I> (<i>&eta;</i>) parameter is a unitless coupling constant
between the shock system and the quantum thermal bath. A small <I>eta</I>
value cannot adjust the quantum temperature fast enough during the
temperature ramping period of shock compression while large <I>eta</I>
leads to big temperature oscillation. A value of <I>eta</I> between 0.3 and
1 is usually appropriate for simulating most systems under shock
compression. We observe that different values of <I>eta</I> lead to almost
the same final thermodynamic state behind the shock, as expected.
</P>
<P>The quantum temperature is updated every <I>beta</I> (<i>&beta;</i>) steps
with an integration time interval <I>beta</I> times longer than the
simulation time step. In that case, <i>etot</i> is taken as its
average over the past <I>beta</I> steps. The temperature of the quantum
thermal bath <i>T</i><sup>qm</sup> changes dynamically according to
the following equation where &Delta;<i>t</i> is the MD time step and
<i>&gamma;</i> is the friction constant which is equal to the inverse
of the <I>damp</I> parameter.
</P>
<center><font size="4"> <i>dT</i><sup>qm</sup>/<i>dt =
&gamma;&eta;</i>&sum;<i><sup>&beta;</sup><sub>l =
1</sub></i>[<i>etot</i>(<i>t-l</i>&Delta;<i>t</i>)-<i>etot</i><sub>0</sub>]/<i>3&beta;Nk</i><sub>B</sub>
</font></center>
<P>The parameter <I>T_init</I> is the initial temperature of the quantum
thermal bath and the system before shock loading.
</P>
<P>For all pressure styles, the simulation box stays orthorhombic in
shape. Parrinello-Rahman boundary conditions (tilted box) are
supported by LAMMPS, but are not implemented for QBMSST.
</P>
<HR>
<P><B>Restart, fix_modify, output, run start/stop, minimize info:</B>
</P>
<P>Because the state of the random number generator is not written to
<A HREF = "restart.html">binary restart files</A>, this fix cannot be restarted
"exactly" in an uninterrupted fashion. However, in a statistical
sense, a restarted simulation should produce similar behaviors of the
system as if it is not interrupted. To achieve such a restart, one
should write explicitly the same value for <I>q</I>, <I>mu</I>, <I>damp</I>, <I>f_max</I>,
<I>N_f</I>, <I>eta</I>, and <I>beta</I> and set <I>tscale</I> = 0 if the system is
compressed during the first run.
</P>
<P>The progress of the QBMSST can be monitored by printing the global
scalar and global vector quantities computed by the fix. The global
vector contains five values in this order:
</P>
<P>[<I>dhugoniot</I>, <I>drayleigh</I>, <I>lagrangian_speed</I>, <I>lagrangian_position</I>,
<I>quantum_temperature</I>]
</P>
<OL><LI><I>dhugoniot</I> is the departure from the Hugoniot (temperature units).
<LI><I>drayleigh</I> is the departure from the Rayleigh line (pressure units).
<LI><I>lagrangian_speed</I> is the laboratory-frame Lagrangian speed (particle velocity) of the computational cell (velocity units).
<LI><I>lagrangian_position</I> is the computational cell position in the reference frame moving at the shock speed. This is the distance of the computational cell behind the shock front.
<LI><I>quantum_temperature</I> is the temperature of the quantum thermal bath <i>T</i><sup>qm</sup>.
</OL>
<P>To print these quantities to the log file with descriptive column
headers, the following LAMMPS commands are suggested. Here the
<A HREF = "fix_modify.html">fix_modify</A> energy command is also enabled to allow
the thermo keyword <I>etotal</I> to print the quantity <i>etot</i>. See
also the <A HREF = "thermo_style.html">thermo_style</A> command.
</P>
<PRE>fix fix_id all msst z
fix_modify fix_id energy yes
variable dhug equal f_fix_id[1]
variable dray equal f_fix_id[2]
variable lgr_vel equal f_fix_id[3]
variable lgr_pos equal f_fix_id[4]
variable T_qm equal f_fix_id[5]
thermo_style custom step temp ke pe lz pzz etotal v_dhug v_dray v_lgr_vel v_lgr_pos v_T_qm f_fix_id
</PRE>
<P>The global scalar under the entry f_fix_id is the quantity of thermo
energy as an extra part of <i>etot</i>. This global scalar and the
vector of 5 quantities can be accessed by various <A HREF = "Section_howto.html#howto_15">output
commands</A>. It is worth noting that the
temp keyword under the <A HREF = "thermo_style.html">thermo_style</A> command print
the instantaneous classical temperature <i>T</i><sup>cl</sup> as
described in the command <A HREF = "fix_qtb.html">fix qtb</A>.
</P>
<HR>
<P><B>Restrictions:</B>
</P>
<P>This fix style is part of the USER-QTB package. It is only enabled if
LAMMPS was built with that package. See the <A HREF = "Section_start.html#start_3">Making
LAMMPS</A> section for more info.
</P>
<P>All cell dimensions must be periodic. This fix can not be used with a
triclinic cell. The QBMSST fix has been tested only for the group-ID
all.
</P>
<HR>
<P><B>Related commands:</B>
</P>
<P><A HREF = "fix_qtb.html">fix qtb</A>, <A HREF = "fix_msst.html">fix msst</A>
</P>
<HR>
<P><B>Default:</B>
</P>
<P>The keyword defaults are q = 10, mu = 0, tscale = 0.01, damp = 1, seed
= 880302, f_max = 200.0, N_f = 100, eta = 1.0, beta = 100, and
T_init=300.0. e0, p0, and v0 are calculated on the first step.
</P>
<HR>
<A NAME = "Goldman"></A>
<P><B>(Goldman)</B> Goldman, Reed and Fried, J. Chem. Phys. 131, 204103 (2009)
</P>
<A NAME = "Qi"></A>
<P><B>(Qi)</B> Qi and Reed, J. Phys. Chem. A 116, 10451 (2012).
</P>
</HTML>

219
doc/fix_qbmsst.txt Normal file
View File

@ -0,0 +1,219 @@
"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
fix qbmsst command :h3
[Syntax:]
fix ID group-ID qbmsst dir shockvel keyword value ... :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
qbmsst = style name of this fix :l
dir = {x} or {y} or {z} :l
shockvel = shock velocity (strictly positive, velocity units) :l
zero or more keyword/value pairs may be appended :l
keyword = {q} or {mu} or {p0} or {v0} or {e0} or {tscale} or {damp} or {seed}or {f_max} or {N_f} or {eta} or {beta} or {T_init} :l
{q} value = cell mass-like parameter (mass^2/distance^4 units)
{mu} value = artificial viscosity (mass/distance/time units)
{p0} value = initial pressure in the shock equations (pressure units)
{v0} value = initial simulation cell volume in the shock equations (distance^3 units)
{e0} value = initial total energy (energy units)
{tscale} value = reduction in initial temperature (unitless fraction between 0.0 and 1.0)
{damp} value = damping parameter (time units) inverse of friction <i>&gamma;</i>
{seed} value = random number seed (positive integer)
{f_max} value = upper cutoff frequency of the vibration spectrum (1/time units)
{N_f} value = number of frequency bins (positive integer)
{eta} value = coupling constant between the shock system and the quantum thermal bath (positive unitless)
{beta} value = the quantum temperature is updated every beta time steps (positive integer)
{T_init} value = quantum temperature for the initial state (temperature units) :pre
:ule
[Examples:]
fix 1 all qbmsst z 0.122 q 25 mu 0.9 tscale 0.01 damp 200 seed 35082 f_max 0.3 N_f 100 eta 1 beta 400 T_init 110 (liquid methane modeled with the REAX force field, real units)
fix 2 all qbmsst z 72 q 40 tscale 0.05 damp 1 seed 47508 f_max 120.0 N_f 100 eta 1.0 beta 500 T_init 300 (quartz modeled with the BKS force field, metal units) :pre
Two example input scripts are given, including shocked alpha quartz
and shocked liquid methane. The input script first equilibrate an
initial state with the quantum thermal bath at the target temperature
and then apply the qbmsst to simulate shock compression with quantum
nuclear correction. The following two figures plot related quantities
for shocked alpha quartz.
:c,image(JPG/qbmsst_init.jpg)
Figure 1. Classical temperature <i>T</i><sup>cl</sup> = &sum;
<i>m<sub>i</sub>v<sub>i</sub><sup>2</sup>/3Nk</i><sub>B</sub> vs. time
for coupling the alpha quartz initial state with the quantum thermal
bath at target quantum temperature <i>T</i><sup>qm</sup> = 300 K. The
NpH ensemble is used for time integration while QTB provides the
colored random force. <i>T</i><sup>cl</sup> converges at the timescale
of {damp} which is set to be 1 ps.
:c,image(JPG/qbmsst_shock.jpg)
Figure 2. Quantum temperature and pressure vs. time for simulating
shocked alpha quartz with the QBMSST. The shock propagates along the z
direction. Restart of the QBMSST command is demonstrated in the
example input script. Thermodynamic quantities stay continuous before
and after the restart.
[Description:]
This command performs the Quantum-Bath coupled Multi-Scale Shock
Technique (QBMSST) integration. See "(Qi)"_#Qi for a detailed
description of this method. The QBMSST provides description of the
thermodynamics and kinetics of shock processes while incorporating
quantum nuclear effects. The {shockvel} setting determines the steady
shock velocity that will be simulated along direction {dir}.
Quantum nuclear effects "(fix qtb)"_fix_qtb.html can be crucial
especially when the temperature of the initial state is below the
classical limit or there is a great change in the zero point energies
between the initial and final states. Theoretical post processing
quantum corrections of shock compressed water and methane have been
reported as much as 30% of the temperatures "(Goldman)"_#Goldman. A
self-consistent method that couples the shock to a quantum thermal
bath described by a colored noise Langevin thermostat has been
developed by Qi et al "(Qi)"_#Qi and applied to shocked methane. The
onset of chemistry is reported to be at a pressure on the shock
Hugoniot that is 40% lower than observed with classical molecular
dynamics.
It is highly recommended that the system be already in an equilibrium
state with a quantum thermal bath at temperature of {T_init}. The fix
command "fix qtb"_fix_qtb.html at constant temperature {T_init} could
be used before applying this command to introduce self-consistent
quantum nuclear effects into the initial state.
The parameters {q}, {mu}, {e0}, {p0}, {v0} and {tscale} are described
in the command "fix msst"_fix_msst.html. The values of {e0}, {p0}, or
{v0} will be calculated on the first step if not specified. The
parameter of {damp}, {f_max}, and {N_f} are described in the command
"fix qtb"_fix_qtb.html.
The fix qbmsst command couples the shock system to a quantum thermal
bath with a rate that is proportional to the change of the total
energy of the shock system, <i>etot</i> - <i>etot</i><sub>0</sub>.
Here <i>etot</i> consists of both the system energy and a thermal
term, see "(Qi)"_#Qi, and <i>etot</i><sub>0</sub> = {e0} is the
initial total energy.
The {eta} (<i>&eta;</i>) parameter is a unitless coupling constant
between the shock system and the quantum thermal bath. A small {eta}
value cannot adjust the quantum temperature fast enough during the
temperature ramping period of shock compression while large {eta}
leads to big temperature oscillation. A value of {eta} between 0.3 and
1 is usually appropriate for simulating most systems under shock
compression. We observe that different values of {eta} lead to almost
the same final thermodynamic state behind the shock, as expected.
The quantum temperature is updated every {beta} (<i>&beta;</i>) steps
with an integration time interval {beta} times longer than the
simulation time step. In that case, <i>etot</i> is taken as its
average over the past {beta} steps. The temperature of the quantum
thermal bath <i>T</i><sup>qm</sup> changes dynamically according to
the following equation where &Delta;<i>t</i> is the MD time step and
<i>&gamma;</i> is the friction constant which is equal to the inverse
of the {damp} parameter.
<center><font size="4"> <i>dT</i><sup>qm</sup>/<i>dt =
&gamma;&eta;</i>&sum;<i><sup>&beta;</sup><sub>l =
1</sub></i>[<i>etot</i>(<i>t-l</i>&Delta;<i>t</i>)-<i>etot</i><sub>0</sub>]/<i>3&beta;Nk</i><sub>B</sub>
</font></center>
The parameter {T_init} is the initial temperature of the quantum
thermal bath and the system before shock loading.
For all pressure styles, the simulation box stays orthorhombic in
shape. Parrinello-Rahman boundary conditions (tilted box) are
supported by LAMMPS, but are not implemented for QBMSST.
:line
[Restart, fix_modify, output, run start/stop, minimize info:]
Because the state of the random number generator is not written to
"binary restart files"_restart.html, this fix cannot be restarted
"exactly" in an uninterrupted fashion. However, in a statistical
sense, a restarted simulation should produce similar behaviors of the
system as if it is not interrupted. To achieve such a restart, one
should write explicitly the same value for {q}, {mu}, {damp}, {f_max},
{N_f}, {eta}, and {beta} and set {tscale} = 0 if the system is
compressed during the first run.
The progress of the QBMSST can be monitored by printing the global
scalar and global vector quantities computed by the fix. The global
vector contains five values in this order:
\[{dhugoniot}, {drayleigh}, {lagrangian_speed}, {lagrangian_position},
{quantum_temperature}\]
{dhugoniot} is the departure from the Hugoniot (temperature units).
{drayleigh} is the departure from the Rayleigh line (pressure units).
{lagrangian_speed} is the laboratory-frame Lagrangian speed (particle velocity) of the computational cell (velocity units).
{lagrangian_position} is the computational cell position in the reference frame moving at the shock speed. This is the distance of the computational cell behind the shock front.
{quantum_temperature} is the temperature of the quantum thermal bath <i>T</i><sup>qm</sup>. :ol
To print these quantities to the log file with descriptive column
headers, the following LAMMPS commands are suggested. Here the
"fix_modify"_fix_modify.html energy command is also enabled to allow
the thermo keyword {etotal} to print the quantity <i>etot</i>. See
also the "thermo_style"_thermo_style.html command.
fix fix_id all msst z
fix_modify fix_id energy yes
variable dhug equal f_fix_id\[1\]
variable dray equal f_fix_id\[2\]
variable lgr_vel equal f_fix_id\[3\]
variable lgr_pos equal f_fix_id\[4\]
variable T_qm equal f_fix_id\[5\]
thermo_style custom step temp ke pe lz pzz etotal v_dhug v_dray v_lgr_vel v_lgr_pos v_T_qm f_fix_id :pre
The global scalar under the entry f_fix_id is the quantity of thermo
energy as an extra part of <i>etot</i>. This global scalar and the
vector of 5 quantities can be accessed by various "output
commands"_Section_howto.html#howto_15. It is worth noting that the
temp keyword under the "thermo_style"_thermo_style.html command print
the instantaneous classical temperature <i>T</i><sup>cl</sup> as
described in the command "fix qtb"_fix_qtb.html.
:line
[Restrictions:]
This fix style is part of the USER-QTB package. It is only enabled if
LAMMPS was built with that package. See the "Making
LAMMPS"_Section_start.html#start_3 section for more info.
All cell dimensions must be periodic. This fix can not be used with a
triclinic cell. The QBMSST fix has been tested only for the group-ID
all.
:line
[Related commands:]
"fix qtb"_fix_qtb.html, "fix msst"_fix_msst.html
:line
[Default:]
The keyword defaults are q = 10, mu = 0, tscale = 0.01, damp = 1, seed
= 880302, f_max = 200.0, N_f = 100, eta = 1.0, beta = 100, and
T_init=300.0. e0, p0, and v0 are calculated on the first step.
:line
:link(Goldman)
[(Goldman)] Goldman, Reed and Fried, J. Chem. Phys. 131, 204103 (2009)
:link(Qi)
[(Qi)] Qi and Reed, J. Phys. Chem. A 116, 10451 (2012).

194
doc/fix_qtb.html Normal file
View File

@ -0,0 +1,194 @@
<HTML>
<CENTER><A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A>
</CENTER>
<HR>
<H3>fix qtb command
</H3>
<P><B>Syntax:</B>
</P>
<PRE>fix ID group-ID qtb keyword value ...
</PRE>
<UL><LI>ID, group-ID are documented in <A HREF = "fix.html">fix</A> command
<LI>qtb = style name of this fix
<LI>zero or more keyword/value pairs may be appended
<LI>keyword = <I>temp</I> or <I>damp</I> or <I>seed</I> or <I>f_max</I> or <I>N_f</I>
<PRE> <I>temp</I> value = target quantum temperature (temperature units)
<I>damp</I> value = damping parameter (time units) inverse of friction <i>&gamma</i>;
<I>seed</I> value = random number seed (positive integer)
<I>f_max</I> value = upper cutoff frequency of the vibration spectrum (1/time units)
<I>N_f</I> value = number of frequency bins (positive integer)
</PRE>
</UL>
<P><B>Examples:</B>
</P>
<PRE>fix 1 all nve
fix 1 all qtb temp 110 damp 200 seed 35082 f_max 0.3 N_f 100 (liquid methane modeled with the REAX force field, real units)
fix 2 all nph iso 1.01325 1.01325 1
fix 2 all qtb temp 300 damp 1 seed 47508 f_max 120.0 N_f 100 (quartz modeled with the BKS force field, metal units)
</PRE>
<P><B>Description:</B>
</P>
<P>This command performs the quantum thermal bath scheme proposed by
<A HREF = "#Dammak">(Dammak)</A> to include self-consistent quantum nuclear effects,
when used in conjunction with the <A HREF = "fix_nve.html">fix nve</A> or <A HREF = "fix_nh.html">fix
nph</A> commands.
</P>
<P>Classical molecular dynamics simulation does not include any quantum
nuclear effect. Quantum treatment of the vibrational modes will
introduce zero point energy into the system, alter the energy power
spectrum and bias the heat capacity from the classical limit. Missing
all the quantum nuclear effects, classical MD cannot model systems at
temperatures lower than their classical limits. This effect is
especially important for materials with a large population of hydrogen
atoms and thus higher classical limits.
</P>
<P>The equation of motion implemented by this command follows a Langevin
form:
</P>
<center><font size="4"><i> m<sub>i</sub>a<sub>i</sub> = f<sub>i</sub>
+ R<sub>i</sub> -
m<sub>i</sub>&gamma;v<sub>i</sub>. </i></font></center>
<P>Here <i>m<sub>i</sub></i>, <i>a<sub>i</sub></i>, <i>f<sub>i</sub>
</i>, <i>R<sub>i</sub></i>, <i>&gamma;</i> and <i>v<sub>i</sub> </i>
represent mass, acceleration, force exerted by all other atoms, random
force, frictional coefficient (the inverse of damping parameter damp),
and velocity. The random force <i>R<sub>i</sub></i> is "colored" so
that any vibrational mode with frequency <i>&omega;</i> will have a
temperature-sensitive energy <i>&theta;</i>(<i>&omega;,T</i>) which
resembles the energy expectation for a quantum harmonic oscillator
with the same natural frequency:
</P>
<center><font size="4"> <i>&theta;</i>(<i>&omega;,T</i>) =
&#8463;&omega;/2 +
&#8463;&omega;[</i>exp(<i>&#8463;&omega;/k</i><sub>B</sub><i>T</i>)<i>-1</i>]<i><sup>-1</sup></i>
</font></center>
<P>To efficiently generate the random forces, we employ the method
of <A HREF = "#Barrat">(Barrat)</A>, that circumvents the need to generate all
random forces for all times before the simulation. The memory
requirement of this approach is less demanding and independent
of the simulation duration. Since the total random force <i>R</i><sub>tot</sub>
does not necessarily vanish for a finite number of atoms,
<i>R<sub>i</sub></i> is replaced by <i>R<sub>i</sub></i> - <i>R</i><sub>tot</sub>/<i>N</i><sub>tot</sub>
to avoid collective motion of the system.
</P>
<P>The <I>temp</I> parameter sets the target quantum temperature. LAMMPS will
still have an output temperature in its thermo style. That is the
instantaneous classical temperature <i>T</i><sup>cl</sup> derived from
the atom velocities at thermal equilibrium. A non-zero
<i>T</i><sup>cl</sup> will be present even when the quantum
temperature approaches zero. This is associated with zero-point energy
at low temperatures.
</P>
<center><font size="4"> <i>T</i><sup>cl</sup> = &sum;
<i>m<sub>i</sub>v<sub>i</sub><sup>2</sup>/3Nk</i><sub>B</sub>
</font></center>
<P>The <I>damp</I> parameter is specified in time units, and it equals the
inverse of the frictional coefficient <i>&gamma;</i>. <i>&gamma;</i>
should be as small as possible but slightly larger than the timescale
of anharmonic coupling in the system which is about 10 ps to 100
ps. When <i>&gamma;</i> is too large, it gives an energy spectrum that
differs from the desired Bose-Einstein spectrum. When <i>&gamma;</i>
is too small, the quantum thermal bath coupling to the system will be
less significant than anharmonic effects, reducing to a classical
limit. We find that setting <i>&gamma;</i> between 5 THz and 1 THz
could be appropriate depending on the system.
</P>
<P>The random number <I>seed</I> is a positive integer used to initiate a
Marsaglia random number generator. Each processor uses the input seed
to generate its own unique seed and its own stream of random
numbers. Thus the dynamics of the system will not be identical on two
runs on different numbers of processors.
</P>
<P>The <I>f_max</I> parameter truncate the noise frequency domain so that
vibrational modes with frequencies higher than <I>f_max</I> will not be
modulated. If we denote &Delta;<i>t</i> as the time interval for the
MD integration, <I>f_max</I> is always reset by the code to make
<i>&alpha;</i> = (int)(2<I>f_max</I>&Delta;<i>t</i>)<sup><i>-1</i></sup> a
positive integer and print out relative information. An appropriate
value for the cutoff frequency <I>f_max</I> would be around 2~3
<i>f</i><sub>D</sub>, where <i>f</i><sub>D</sub> is the Debye
frequency.
</P>
<P>The <I>N_f</I> parameter is the frequency grid size, the number of points
from 0 to <I>f_max</I> in the frequency domain that will be
sampled. <i>3&times;2</i> <I>N_f</I> per-atom random numbers are required
in the random force generation and there could be as many atoms as in
the whole simulation that can migrate into every individual
processor. A larger <I>N_f</I> provides a more accurate sampling of the
spectrum while consumes more memory. With fixed <I>f_max</I> and
<i>&gamma;</i>, <I>N_f</I> should be big enough to converge the classical
temperature <i>T</i><sup>cl</sup> as a function of target quantum bath
temperature. Memory usage per processor could be from 10 to 100
Mbytes.
</P>
<P>IMPORTANT NOTE: Unlike the <A HREF = "fix_nh.html">fix nvt</A> command which
performs Nose/Hoover thermostatting AND time integration, this fix
does NOT perform time integration. It only modifies forces to a
colored thermostat. Thus you must use a separate time integration fix,
like <A HREF = "fix_nve.html">fix nve</A> or <A HREF = "fix_nh.html">fix nph</A> to actually
update the velocities and positions of atoms (as shown in the
examples). Likewise, this fix should not normally be used with other
fixes or commands that also specify system temperatures , e.g. <A HREF = "fix_nh.html">fix
nvt</A> and <A HREF = "fix_temp_rescale.html">fix temp/rescale</A>.
</P>
<HR>
<P><B>Restart, fix_modify, output, run start/stop, minimizie info:</B>
</P>
<P>No information about this fix is written to <A HREF = "restart.html">binary restart
files</A>. Because the state of the random number generator
is not saved in restart files, this means you cannot do "exact"
restarts with this fix. However, in a statistical sense, a restarted
simulation should produce similar behaviors of the system.
</P>
<P>This fix is not invoked during <A HREF = "minimize.html">energy minimization</A>.
</P>
<HR>
<P><B>Restrictions:</B>
</P>
<P>This fix style is part of the USER-QTB package. It is only enabled if
LAMMPS was built with that package. See the <A HREF = "Section_start.html#start_3">Making
LAMMPS</A> section for more info.
</P>
<HR>
<P><B>Related commands:</B>
</P>
<P><A HREF = "fix_nve.html">fix nve</A>, <A HREF = "fix_nh.html">fix nph</A>, <A HREF = "fix_langevin.html">fix
langevin</A>, <A HREF = "fix_qbmsst.html">fix qbmsst</A>
</P>
<HR>
<P><B>Default:</B>
</P>
<P>The keyword defaults are temp = 300, damp = 1, seed = 880302,
f_max=200.0 and N_f = 100.
</P>
<HR>
<A NAME = "Dammak"></A>
<P><B>(Dammak)</B> Dammak, Chalopin, Laroche, Hayoun, and Greffet, Phys Rev
Lett, 103, 190601 (2009).
</P>
<A NAME = "Barrat"></A>
<P><B>(Barrat)</B> Barrat and Rodney, J. Stat. Phys, 144, 679 (2011).
</P>
</HTML>

182
doc/fix_qtb.txt Normal file
View File

@ -0,0 +1,182 @@
"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
fix qtb command :h3
[Syntax:]
fix ID group-ID qtb keyword value ... :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
qtb = style name of this fix :l
zero or more keyword/value pairs may be appended :l
keyword = {temp} or {damp} or {seed} or {f_max} or {N_f} :l
{temp} value = target quantum temperature (temperature units)
{damp} value = damping parameter (time units) inverse of friction <i>&gamma</i>;
{seed} value = random number seed (positive integer)
{f_max} value = upper cutoff frequency of the vibration spectrum (1/time units)
{N_f} value = number of frequency bins (positive integer) :pre
:ule
[Examples:]
fix 1 all nve
fix 1 all qtb temp 110 damp 200 seed 35082 f_max 0.3 N_f 100 (liquid methane modeled with the REAX force field, real units)
fix 2 all nph iso 1.01325 1.01325 1
fix 2 all qtb temp 300 damp 1 seed 47508 f_max 120.0 N_f 100 (quartz modeled with the BKS force field, metal units) :pre
[Description:]
This command performs the quantum thermal bath scheme proposed by
"(Dammak)"_#Dammak to include self-consistent quantum nuclear effects,
when used in conjunction with the "fix nve"_fix_nve.html or "fix
nph"_fix_nh.html commands.
Classical molecular dynamics simulation does not include any quantum
nuclear effect. Quantum treatment of the vibrational modes will
introduce zero point energy into the system, alter the energy power
spectrum and bias the heat capacity from the classical limit. Missing
all the quantum nuclear effects, classical MD cannot model systems at
temperatures lower than their classical limits. This effect is
especially important for materials with a large population of hydrogen
atoms and thus higher classical limits.
The equation of motion implemented by this command follows a Langevin
form:
<center><font size="4"><i> m<sub>i</sub>a<sub>i</sub> = f<sub>i</sub>
+ R<sub>i</sub> -
m<sub>i</sub>&gamma;v<sub>i</sub>. </i></font></center>
Here <i>m<sub>i</sub></i>, <i>a<sub>i</sub></i>, <i>f<sub>i</sub>
</i>, <i>R<sub>i</sub></i>, <i>&gamma;</i> and <i>v<sub>i</sub> </i>
represent mass, acceleration, force exerted by all other atoms, random
force, frictional coefficient (the inverse of damping parameter damp),
and velocity. The random force <i>R<sub>i</sub></i> is "colored" so
that any vibrational mode with frequency <i>&omega;</i> will have a
temperature-sensitive energy <i>&theta;</i>(<i>&omega;,T</i>) which
resembles the energy expectation for a quantum harmonic oscillator
with the same natural frequency:
<center><font size="4"> <i>&theta;</i>(<i>&omega;,T</i>) =
&#8463;&omega;/2 +
&#8463;&omega;[</i>exp(<i>&#8463;&omega;/k</i><sub>B</sub><i>T</i>)<i>-1</i>]<i><sup>-1</sup></i>
</font></center>
To efficiently generate the random forces, we employ the method
of "(Barrat)"_#Barrat, that circumvents the need to generate all
random forces for all times before the simulation. The memory
requirement of this approach is less demanding and independent
of the simulation duration. Since the total random force <i>R</i><sub>tot</sub>
does not necessarily vanish for a finite number of atoms,
<i>R<sub>i</sub></i> is replaced by <i>R<sub>i</sub></i> - <i>R</i><sub>tot</sub>/<i>N</i><sub>tot</sub>
to avoid collective motion of the system.
The {temp} parameter sets the target quantum temperature. LAMMPS will
still have an output temperature in its thermo style. That is the
instantaneous classical temperature <i>T</i><sup>cl</sup> derived from
the atom velocities at thermal equilibrium. A non-zero
<i>T</i><sup>cl</sup> will be present even when the quantum
temperature approaches zero. This is associated with zero-point energy
at low temperatures.
<center><font size="4"> <i>T</i><sup>cl</sup> = &sum;
<i>m<sub>i</sub>v<sub>i</sub><sup>2</sup>/3Nk</i><sub>B</sub>
</font></center>
The {damp} parameter is specified in time units, and it equals the
inverse of the frictional coefficient <i>&gamma;</i>. <i>&gamma;</i>
should be as small as possible but slightly larger than the timescale
of anharmonic coupling in the system which is about 10 ps to 100
ps. When <i>&gamma;</i> is too large, it gives an energy spectrum that
differs from the desired Bose-Einstein spectrum. When <i>&gamma;</i>
is too small, the quantum thermal bath coupling to the system will be
less significant than anharmonic effects, reducing to a classical
limit. We find that setting <i>&gamma;</i> between 5 THz and 1 THz
could be appropriate depending on the system.
The random number {seed} is a positive integer used to initiate a
Marsaglia random number generator. Each processor uses the input seed
to generate its own unique seed and its own stream of random
numbers. Thus the dynamics of the system will not be identical on two
runs on different numbers of processors.
The {f_max} parameter truncate the noise frequency domain so that
vibrational modes with frequencies higher than {f_max} will not be
modulated. If we denote &Delta;<i>t</i> as the time interval for the
MD integration, {f_max} is always reset by the code to make
<i>&alpha;</i> = (int)(2{f_max}&Delta;<i>t</i>)<sup><i>-1</i></sup> a
positive integer and print out relative information. An appropriate
value for the cutoff frequency {f_max} would be around 2~3
<i>f</i><sub>D</sub>, where <i>f</i><sub>D</sub> is the Debye
frequency.
The {N_f} parameter is the frequency grid size, the number of points
from 0 to {f_max} in the frequency domain that will be
sampled. <i>3&times;2</i> {N_f} per-atom random numbers are required
in the random force generation and there could be as many atoms as in
the whole simulation that can migrate into every individual
processor. A larger {N_f} provides a more accurate sampling of the
spectrum while consumes more memory. With fixed {f_max} and
<i>&gamma;</i>, {N_f} should be big enough to converge the classical
temperature <i>T</i><sup>cl</sup> as a function of target quantum bath
temperature. Memory usage per processor could be from 10 to 100
Mbytes.
IMPORTANT NOTE: Unlike the "fix nvt"_fix_nh.html command which
performs Nose/Hoover thermostatting AND time integration, this fix
does NOT perform time integration. It only modifies forces to a
colored thermostat. Thus you must use a separate time integration fix,
like "fix nve"_fix_nve.html or "fix nph"_fix_nh.html to actually
update the velocities and positions of atoms (as shown in the
examples). Likewise, this fix should not normally be used with other
fixes or commands that also specify system temperatures , e.g. "fix
nvt"_fix_nh.html and "fix temp/rescale"_fix_temp_rescale.html.
:line
[Restart, fix_modify, output, run start/stop, minimizie info:]
No information about this fix is written to "binary restart
files"_restart.html. Because the state of the random number generator
is not saved in restart files, this means you cannot do "exact"
restarts with this fix. However, in a statistical sense, a restarted
simulation should produce similar behaviors of the system.
This fix is not invoked during "energy minimization"_minimize.html.
:line
[Restrictions:]
This fix style is part of the USER-QTB package. It is only enabled if
LAMMPS was built with that package. See the "Making
LAMMPS"_Section_start.html#start_3 section for more info.
:line
[Related commands:]
"fix nve"_fix_nve.html, "fix nph"_fix_nh.html, "fix
langevin"_fix_langevin.html, "fix qbmsst"_fix_qbmsst.html
:line
[Default:]
The keyword defaults are temp = 300, damp = 1, seed = 880302,
f_max=200.0 and N_f = 100.
:line
:link(Dammak)
[(Dammak)] Dammak, Chalopin, Laroche, Hayoun, and Greffet, Phys Rev
Lett, 103, 190601 (2009).
:link(Barrat)
[(Barrat)] Barrat and Rodney, J. Stat. Phys, 144, 679 (2011).