fix tgnvt/drude command¶
fix tgnpt/drude command¶
Syntax¶
fix ID group-ID style_name keyword values ...
ID, group-ID are documented in fix command
style_name = tgnvt/drude or tgnpt/drude
one or more keyword/values pairs may be appended
keyword = temp iso or aniso or tri or x or y or z or xy or yz or xz or couple or tchain or pchain or mtk or tloop or ploop or nreset or scalexy or scaleyz or scalexz or flip or fixedpoint temp values = Tstart Tstop Tdamp Tdrude Tdamp_drude Tstart, Tstop = external temperature at start/end of run (temperature units) Tdamp = temperature damping parameter (time units) Tdrude = desired temperature of Drude oscillators (temperature units) Tdamp_drude = temperature damping parameter for Drude oscillators (time units) iso or aniso or tri values = Pstart Pstop Pdamp Pstart,Pstop = scalar external pressure at start/end of run (pressure units) Pdamp = pressure damping parameter (time units) x or y or z or xy or yz or xz values = Pstart Pstop Pdamp Pstart,Pstop = external stress tensor component at start/end of run (pressure units) Pdamp = stress damping parameter (time units) couple = none or xyz or xy or yz or xz tchain value = N N = length of thermostat chain (1 = single thermostat) pchain value = N N length of thermostat chain on barostat (0 = no thermostat) mtk value = yes or no = add in MTK adjustment term or not tloop value = M M = number of sub-cycles to perform on thermostat ploop value = M M = number of sub-cycles to perform on barostat thermostat nreset value = reset reference cell every this many timesteps scalexy value = yes or no = scale xy with ly scaleyz value = yes or no = scale yz with lz scalexz value = yes or no = scale xz with lz flip value = yes or no = allow or disallow box flips when it becomes highly skewed fixedpoint values = x y z x,y,z = perform barostat dilation/contraction around this point (distance units)
Examples¶
comm_modify vel yes
fix 1 all tgnvt/drude temp 300.0 300.0 100.0 1.0 20.0
fix 1 water tgnpt/drude temp 300.0 300.0 100.0 1.0 20.0 iso 0.0 0.0 1000.0
fix 2 jello tgnpt/drude temp 300.0 300.0 100.0 1.0 20.0 tri 5.0 5.0 1000.0
fix 2 ice tgnpt/drude temp 250.0 250.0 100.0 1.0 20.0 x 1.0 1.0 0.5 y 2.0 2.0 0.5 z 3.0 3.0 0.5 yz 0.1 0.1 0.5 xz 0.2 0.2 0.5 xy 0.3 0.3 0.5 nreset 1000
Example input scripts available: examples/PACKAGES/drude
Description¶
These commands are variants of the Nose-Hoover fix styles fix nvt and fix npt for thermalized Drude polarizable models. They apply temperature-grouped Nose-Hoover thermostat (TGNH) proposed by (Son). When there are fast vibrational modes with frequencies close to Drude oscillators (e.g. double bonds or out-of-plane torsions), this thermostat can provide better kinetic energy equipartitioning.
The difference between TGNH and the original Nose-Hoover thermostat is that, TGNH separates the kinetic energy of the group into three contributions: molecular center of mass (COM) motion, motion of COM of atom-Drude pairs or non-polarizable atoms relative to molecular COM, and relative motion of atom-Drude pairs. An independent Nose-Hoover chain is applied to each type of motion. The temperatures for these three types of motion are denoted as molecular translational temperature (\(T_\mathrm{M}\)), real atomic temperature (\(T_\mathrm{R}\)) and Drude temperature (\(T_\mathrm{D}\)), which are defined in terms of their associated degrees of freedom (DOF):
Here \(N_\mathrm{mol}\) and \(N_\mathrm{mol,sys}\) are the numbers of molecules in the group and in the whole system, respectively. \(N_\mathrm{real}\) is the number of atom-Drude pairs and non-polarizable atoms in the group. \(N_\mathrm{drude}\) is the number of Drude particles in the group. \(N_\mathrm{DOF}\) is the DOF of the group. \(M_i\) and \(V_i\) are the mass and the COM velocity of the i-th molecule. \(m_i\) is the mass of the i-th atom-Drude pair or non-polarizable atom. \(v_i\) is the velocity of COM of i-th atom-Drude pair or non-polarizable atom. \(v_{M,i}\) is the COM velocity of the molecule the i-th atom-Drude pair or non-polarizable atom belongs to. \(m_i^\prime\) and \(v_i^\prime\) are the reduced mass and the relative velocity of the i-th atom-Drude pair.
Note
These fixes require that each atom knows whether it is a Drude particle or not. You must therefore use the fix drude command to specify the Drude status of each atom type.
Because the TGNH thermostat thermostats the molecular COM motion, all atoms belonging to the same molecule must be in the same group. That is, these fixes can not be applied to a subset of a molecule.
For this fix to act correctly, ghost atoms need to know their velocity. You must use the comm_modify command to enable this.
These fixes assume that the translational DOF of the whole system is removed. It is therefore recommended to invoke fix momentum command so that the \(T_\mathrm{M}\) is calculated correctly.
The thermostat parameters are specified using the temp keyword. The thermostat is applied to only the translational DOF for the particles. The translational DOF can also have a bias velocity removed before thermostatting takes place; see the description below. The desired temperature for molecular and real atomic motion is a ramped value during the run from Tstart to Tstop. The Tdamp parameter is specified in time units and determines how rapidly the temperature is relaxed. For example, a value of 10.0 means to relax the temperature in a timespan of (roughly) 10 time units (e.g. \(\tau\) or fs or ps - see the units command). The parameter Tdrude is the desired temperature for Drude motion at each timestep. Similar to Tdamp, the Tdamp_drude parameter determines the relaxation speed for Drude motion. Fix group are the only ones whose velocities and positions are updated by the velocity/position update portion of the integration. Other thermostat-related keywords are tchain and tloop, which are detailed in fix nvt.
Note
A Nose-Hoover thermostat will not work well for arbitrary values of Tdamp. If Tdamp is too small, the temperature can fluctuate wildly; if it is too large, the temperature will take a very long time to equilibrate. A good choice for many models is a Tdamp of around 100 timesteps. A smaller Tdamp_drude value would be required to maintain Drude motion at low temperature.
fix 1 all nvt temp 300.0 300.0 $(100.0*dt) 1.0 $(20.0*dt)
The barostat parameters for fix style tgnpt/drude is specified using one or more of the iso, aniso, tri, x, y, z, xy, xz, yz, and couple keywords. These keywords give you the ability to specify all 6 components of an external stress tensor, and to couple various of these components together so that the dimensions they represent are varied together during a constant-pressure simulation. Other barostat-related keywords are pchain, mtk, ploop, nreset, scalexy, scaleyz, scalexz, flipand fixedpoint. The meaning of barostat parameters are detailed in fix npt.
Regardless of what atoms are in the fix group (the only atoms which are time integrated), a global pressure or stress tensor is computed for all atoms. Similarly, when the size of the simulation box is changed, all atoms are re-scaled to new positions.
Note
Unlike the fix temp/berendsen command which performs thermostatting but NO time integration, these fixes perform thermostatting/barostatting AND time integration. Thus you should not use any other time integration fix, such as fix nve on atoms to which this fix is applied. Likewise, these fixes should not be used on atoms that also have their temperature controlled by another fix - e.g. by fix langevin/drude command.
See the Howto thermostat and Howto barostat doc pages for a discussion of different ways to compute temperature and perform thermostatting and barostatting.
Like other fixes that perform thermostatting, these fixes can be used with compute commands that calculate a temperature after removing a “bias” from the atom velocities. This is not done by default, but only if the fix_modify command is used to assign a temperature compute to this fix that includes such a bias term. See the doc pages for individual compute commands to determine which ones include a bias. In this case, the thermostat works in the following manner: the current temperature is calculated taking the bias into account, bias is removed from each atom, thermostatting is performed on the remaining thermal DOF, and the bias is added back in.
Note
However, not all temperature compute commands are valid to be used with these fixes. Precisely, only temperature compute that does not modify the DOF of the group can be used. E.g. compute temp/ramp and compute viscosity/cos compute the kinetic energy after remove a velocity gradient without affecting the DOF of the group, then they can be invoked in this way. In contrast, compute temp/partial may remove the DOF at one or more dimensions, therefore it cannot be used with these fixes.
Restart, fix_modify, output, run start/stop, minimize info¶
These fixes writes the state of all the thermostat and barostat variables to binary restart files. See the read_restart 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 fix_modify temp and press options are supported by these fixes. You can use them to assign a compute you have defined to this fix which will be used in its thermostatting or barostatting procedure, as described above. If you do this, note that the kinetic energy derived from the compute temperature should be consistent with the virial term computed using all atoms for the pressure. LAMMPS will warn you if you choose to compute temperature on a subset of atoms.
Note
If both the temp and press keywords are used in a single thermo_modify command (or in two separate commands), then the order in which the keywords are specified is important. Note that a pressure compute defines its own temperature compute as an argument when it is specified. The temp keyword will override this (for the pressure compute being used by fix npt), but only if the temp keyword comes after the press keyword. If the temp keyword comes before the press keyword, then the new pressure compute specified by the press keyword will be unaffected by the temp setting.
The cumulative energy change in the system imposed by these fixes, due to thermostatting and/or barostatting, are included in the thermodynamic output keywords ecouple and econserve. See the thermo_style page for details.
These fixes compute a global scalar which can be accessed by various output commands. The scalar is the same cumulative energy change due to this fix described in the previous paragraph. The scalar value calculated by this fix is “extensive”.
These fixes also compute a global vector of quantities, which can be accessed by various output commands. The vector values are “intensive”. The vector stores the three temperatures \(T_\mathrm{M}\), \(T_\mathrm{R}\) and \(T_\mathrm{D}\).
These fixes can ramp their external temperature and pressure over multiple runs, using the start and stop keywords of the run command. See the run command for details of how to do this.
These fixes are not invoked during energy minimization.
Restrictions¶
These fixes are only available when LAMMPS was built with the DRUDE package. These fixes cannot be used with dynamic groups as defined by the group command. These fixes cannot be used in 2D simulations.
X, y, z cannot be barostatted if the associated dimension is not periodic. Xy, xz, and yz can only be barostatted if the simulation domain is triclinic and the second dimension in the keyword (y dimension in xy) is periodic. The create_box, read data, and read_restart commands specify whether the simulation box is orthogonal or non-orthogonal (triclinic) and explain the meaning of the xy,xz,yz tilt factors.
For the temp keyword, the final Tstop cannot be 0.0 since it would make the external T = 0.0 at some timestep during the simulation which is not allowed in the Nose/Hoover formulation.
The scaleyz yes, scalexz yes, and scalexy yes options can only be used if the second dimension in the keyword is periodic, and if the tilt factor is not coupled to the barostat via keywords tri, yz, xz, and xy.
Default¶
The keyword defaults are tchain = 3, pchain = 3, mtk = yes, tloop = 1, ploop = 1, nreset = 0, couple = none, flip = yes, scaleyz = scalexz = scalexy = yes if periodic in second dimension and not coupled to barostat, otherwise no.
(Son) Son, McDaniel, Cui and Yethiraj, J Phys Chem Lett, 10, 7523 (2019).