diff --git a/doc/src/fix_msst.txt b/doc/src/fix_msst.txt index 025c73389..310692669 100644 --- a/doc/src/fix_msst.txt +++ b/doc/src/fix_msst.txt @@ -1,193 +1,193 @@ "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 msst command :h3 [Syntax:] fix ID group-ID msst dir shockvel keyword value ... :pre ID, group-ID are documented in "fix"_fix.html command :ulb,l msst = style name of this fix :l dir = {x} or {y} or {z} :l shockvel = shock velocity (strictly positive, distance/time 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 {beta} or {dftb} :l {q} value = cell mass-like parameter (mass^2/distance^4 units) {mu} value = artificial viscosity (mass/length/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) {dftb} value = {yes} or {no} for whether using MSST in conjunction with DFTB+ {beta} value = scale factor for improved energy conservation :pre :ule [Examples:] fix 1 all msst y 100.0 q 1.0e5 mu 1.0e5 fix 2 all msst z 50.0 q 1.0e4 mu 1.0e4 v0 4.3419e+03 p0 3.7797e+03 e0 -9.72360e+02 tscale 0.01 fix 1 all msst y 100.0 q 1.0e5 mu 1.0e5 dftb yes beta 0.5 :pre [Description:] This command performs the Multi-Scale Shock Technique (MSST) integration to update positions and velocities each timestep to mimic a compressive shock wave passing over the system. See "(Reed)"_#Reed for a detailed description of this method. The MSST varies the cell volume and temperature in such a way as to restrain the system to the shock Hugoniot and the Rayleigh line. These restraints correspond to the macroscopic conservation laws dictated by a shock front. {shockvel} determines the steady shock velocity that will be simulated. To perform a simulation, choose a value of {q} that provides volume compression on the timescale of 100 fs to 1 ps. If the volume is not compressing, either the shock speed is chosen to be below the material sound speed or {p0} has been chosen inaccurately. Volume compression at the start can be sped up by using a non-zero value of {tscale}. Use the smallest value of {tscale} that results in compression. Under some special high-symmetry conditions, the pressure (volume) and/or temperature of the system may oscillate for many cycles even with an appropriate choice of mass-like parameter {q}. Such oscillations have physical significance in some cases. The optional {mu} keyword adds an artificial viscosity that helps break the system symmetry to equilibrate to the shock Hugoniot and Rayleigh line more rapidly in such cases. The keyword {tscale} is a factor between 0 and 1 that determines what fraction of thermal kinetic energy is converted to compressive strain kinetic energy at the start of the simulation. Setting this parameter to a non-zero value may assist in compression at the start of simulations where it is slow to occur. If keywords {e0}, {p0},or {v0} are not supplied, these quantities will be calculated on the first step, after the energy specified by {tscale} is removed. The value of {e0} is not used in the dynamical equations, but is used in calculating the deviation from the Hugoniot. The keyword {beta} is a scaling term that can be added to the MSST ionic equations of motion to account for drift in the conserved quantity during long timescale simulations, similar to a Berendson -thermostat. See "(Reed)"_#Reed and "(Goldman)"_#Goldman for more +thermostat. See "(Reed)"_#Reed and "(Goldman)"_#Goldman2 for more details. The value of {beta} must be between 0.0 and 1.0 inclusive. A value of 0.0 means no contribution, a value of 1.0 means a full contribution. Values of shockvel less than a critical value determined by the material response will not have compressive solutions. This will be reflected in lack of significant change of the volume in the MSST. For all pressure styles, the simulation box stays orthogonal in shape. Parrinello-Rahman boundary conditions (tilted box) are supported by LAMMPS, but are not implemented for MSST. This fix computes a temperature and pressure and potential energy each timestep. To do this, the fix creates its own computes of style "temp" "pressure", and "pe", as if these commands had been issued: compute fix-ID_MSST_temp all temp compute fix-ID_MSST_press all pressure fix-ID_MSST_temp :pre compute fix-ID_MSST_pe all pe :pre See the "compute temp"_compute_temp.html and "compute pressure"_compute_pressure.html commands for details. Note that the IDs of the new computes are the fix-ID + "_MSST_temp" or "_MSST_press" or "_MSST_pe". The group for the new computes is "all". :line The {dftb} keyword is to allow this fix to be used when LAMMPS is being driven by DFTB+, a density-functional tight-binding code. If the keyword {dftb} is used with a value of {yes}, then the MSST equations are altered to account for the electron entropy contribution to the Hugonio relations and total energy. See "(Reed2)"_#Reed2 and -"(Goldman)"_#Goldman for details on this contribution. In this case, +"(Goldman)"_#Goldman2 for details on this contribution. In this case, you must define a "fix external"_fix_external.html command in your input script, which is used to callback to DFTB+ during the LAMMPS timestepping. DFTB+ will communicate its info to LAMMPS via that fix. :line [Restart, fix_modify, output, run start/stop, minimize info:] This fix writes the state of all internal variables to "binary restart files"_restart.html. See the "read_restart"_read_restart.html command for info on how to re-specify a fix in an input script that reads a restart file, so that the operation of the fix continues in an uninterrupted fashion. The progress of the MSST can be monitored by printing the global scalar and global vector quantities computed by the fix. The scalar is the cumulative energy change due to the fix. This is also the energy added to the potential energy by the "fix_modify"_fix_modify.html {energy} command. With this command, the thermo keyword {etotal} prints the conserved quantity of the MSST dynamic equations. This can be used to test if the MD timestep is sufficiently small for accurate integration of the dynamic equations. See also "thermo_style"_thermo_style.html command. The global vector contains four values in this order: \[{dhugoniot}, {drayleigh}, {lagrangian_speed}, {lagrangian_position}\] {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 usually a good estimate of distance of the computational cell behind the shock front. :ol To print these quantities to the log file with descriptive column headers, the following LAMMPS commands are suggested: fix msst all msst z fix_modify msst energy yes variable dhug equal f_msst\[1\] variable dray equal f_msst\[2\] variable lgr_vel equal f_msst\[3\] variable lgr_pos equal f_msst\[4\] thermo_style custom step temp ke pe lz pzz etotal v_dhug v_dray v_lgr_vel v_lgr_pos f_msst :pre These fixes compute a global scalar and a global vector of 4 quantities, which can be accessed by various "output commands"_Section_howto.html#howto_15. The scalar values calculated by this fix are "extensive"; the vector values are "intensive". [Restrictions:] This fix style is part of the SHOCK 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 MSST fix has been tested only for the group-ID all. [Related commands:] "fix nphug"_fix_nphug.html, "fix deform"_fix_deform.html [Default:] The keyword defaults are q = 10, mu = 0, tscale = 0.01, dftb = no, beta = 0.0. Note that p0, v0, and e0 are calculated on the first timestep. :line :link(Reed) [(Reed)] Reed, Fried, and Joannopoulos, Phys. Rev. Lett., 90, 235503 (2003). :link(Reed2) [(Reed2)] Reed, J. Phys. Chem. C, 116, 2205 (2012). -:link(Goldman) +:link(Goldman2) [(Goldman)] Goldman, Srinivasan, Hamel, Fried, Gaus, and Elstner, J. Phys. Chem. C, 117, 7885 (2013). diff --git a/doc/src/fix_qbmsst.txt b/doc/src/fix_qbmsst.txt index 468206a57..2c116fb0f 100644 --- a/doc/src/fix_qbmsst.txt +++ b/doc/src/fix_qbmsst.txt @@ -1,219 +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 γ {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 Tcl = ∑ mivi2/3NkB vs. time for coupling the alpha quartz initial state with the quantum thermal bath at target quantum temperature Tqm = 300 K. The NpH ensemble is used for time integration while QTB provides the colored random force. Tcl 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 +reported as much as 30% of the temperatures "(Goldman)"_#Goldman1. 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, etot - etot0. Here etot consists of both the system energy and a thermal term, see "(Qi)"_#Qi, and etot0 = {e0} is the initial total energy. The {eta} (η) 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} (β) steps with an integration time interval {beta} times longer than the simulation time step. In that case, etot is taken as its average over the past {beta} steps. The temperature of the quantum thermal bath Tqm changes dynamically according to the following equation where Δt is the MD time step and γ is the friction constant which is equal to the inverse of the {damp} parameter.
dTqm/dt = γη∑βl = 1[etot(t-lΔt)-etot0]/3βNkB
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 Tqm. :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 etot. 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 etot. 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 Tcl 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) +:link(Goldman1) [(Goldman)] Goldman, Reed and Fried, J. Chem. Phys. 131, 204103 (2009) :link(Qi) [(Qi)] Qi and Reed, J. Phys. Chem. A 116, 10451 (2012).