+Syntax: +
+fix ID group-ID qbmsst dir shockvel keyword value ... ++
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) ++ +
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) ++
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. +
+
+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. +
+
+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) 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) 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). 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) 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 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. 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. +
+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), 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. +
+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. +
+Restart, fix_modify, output, run start/stop, minimize info: +
+Because the state of the random number generator is not written to +binary restart files, 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] +
+To print these quantities to the log file with descriptive column +headers, the following LAMMPS commands are suggested. Here the +fix_modify energy command is also enabled to allow +the thermo keyword etotal to print the quantity etot. See +also the thermo_style 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 ++
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. It is worth noting that the +temp keyword under the thermo_style command print +the instantaneous classical temperature Tcl as +described in the command fix qtb. +
+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 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. +
+Related commands: +
+ +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. +
+(Goldman) Goldman, Reed and Fried, J. Chem. Phys. 131, 204103 (2009) +
+ + +(Qi) Qi and Reed, J. Phys. Chem. A 116, 10451 (2012). +
+ diff --git a/doc/fix_qbmsst.txt b/doc/fix_qbmsst.txt new file mode 100644 index 0000000000..da18130c23 --- /dev/null +++ b/doc/fix_qbmsst.txt @@ -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 γ + {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 +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. + +Syntax: +
+fix ID group-ID qtb keyword value ... ++
temp value = target quantum temperature (temperature units) + 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) ++ +
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) ++
Description: +
+This command performs the quantum thermal bath scheme proposed by +(Dammak) to include self-consistent quantum nuclear effects, +when used in conjunction with the fix nve or fix +nph 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: +
+Here mi, ai, fi +, Ri, γ and vi +represent mass, acceleration, force exerted by all other atoms, random +force, frictional coefficient (the inverse of damping parameter damp), +and velocity. The random force Ri is "colored" so +that any vibrational mode with frequency ω will have a +temperature-sensitive energy θ(ω,T) which +resembles the energy expectation for a quantum harmonic oscillator +with the same natural frequency: +
+To efficiently generate the random forces, we employ the method +of (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 Rtot +does not necessarily vanish for a finite number of atoms, +Ri is replaced by Ri - Rtot/Ntot +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 Tcl derived from +the atom velocities at thermal equilibrium. A non-zero +Tcl will be present even when the quantum +temperature approaches zero. This is associated with zero-point energy +at low temperatures. +
+The damp parameter is specified in time units, and it equals the +inverse of the frictional coefficient γ. γ +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 γ is too large, it gives an energy spectrum that +differs from the desired Bose-Einstein spectrum. When γ +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 γ 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 Δt as the time interval for the +MD integration, f_max is always reset by the code to make +α = (int)(2f_maxΔt)-1 a +positive integer and print out relative information. An appropriate +value for the cutoff frequency f_max would be around 2~3 +fD, where fD 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. 3×2 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 +γ, N_f should be big enough to converge the classical +temperature Tcl 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 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 or fix nph 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 and fix temp/rescale. +
+Restart, fix_modify, output, run start/stop, minimizie info: +
+No information about this fix is written to binary restart +files. 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. +
+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 for more info. +
+Related commands: +
+fix nve, fix nph, fix +langevin, fix qbmsst +
+Default: +
+The keyword defaults are temp = 300, damp = 1, seed = 880302, +f_max=200.0 and N_f = 100. +
+(Dammak) Dammak, Chalopin, Laroche, Hayoun, and Greffet, Phys Rev +Lett, 103, 190601 (2009). +
+ + +(Barrat) Barrat and Rodney, J. Stat. Phys, 144, 679 (2011). +
+ diff --git a/doc/fix_qtb.txt b/doc/fix_qtb.txt new file mode 100644 index 0000000000..837dd514eb --- /dev/null +++ b/doc/fix_qtb.txt @@ -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 γ + {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: + +