1.. index:: fix npt/cauchy 2 3fix npt/cauchy command 4====================== 5 6Syntax 7"""""" 8 9.. parsed-literal:: 10 11 fix ID group-ID style_name keyword value ... 12 13* ID, group-ID are documented in :doc:`fix <fix>` command 14* style_name = *npt/cauchy* 15* one or more keyword/value pairs may be appended 16* keyword = *temp* or *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 *drag* or *dilate* or *scalexy* or *scaleyz* or *scalexz* or *flip* or *fixedpoint* 17 18 .. parsed-literal:: 19 20 *temp* values = Tstart Tstop Tdamp 21 Tstart,Tstop = external temperature at start/end of run 22 Tdamp = temperature damping parameter (time units) 23 *iso* or *aniso* or *tri* values = Pstart Pstop Pdamp 24 Pstart,Pstop = scalar external pressure at start/end of run (pressure units) 25 Pdamp = pressure damping parameter (time units) 26 *x* or *y* or *z* or *xy* or *yz* or *xz* values = Pstart Pstop Pdamp 27 Pstart,Pstop = external stress tensor component at start/end of run (pressure units) 28 Pdamp = stress damping parameter (time units) 29 *couple* = *none* or *xyz* or *xy* or *yz* or *xz* 30 *tchain* value = N 31 N = length of thermostat chain (1 = single thermostat) 32 *pchain* values = N 33 N length of thermostat chain on barostat (0 = no thermostat) 34 *mtk* value = *yes* or *no* = add in MTK adjustment term or not 35 *tloop* value = M 36 M = number of sub-cycles to perform on thermostat 37 *ploop* value = M 38 M = number of sub-cycles to perform on barostat thermostat 39 *nreset* value = reset reference cell every this many timesteps 40 *drag* value = Df 41 Df = drag factor added to barostat/thermostat (0.0 = no drag) 42 *dilate* value = dilate-group-ID 43 dilate-group-ID = only dilate atoms in this group due to barostat volume changes 44 *scalexy* value = *yes* or *no* = scale xy with ly 45 *scaleyz* value = *yes* or *no* = scale yz with lz 46 *scalexz* value = *yes* or *no* = scale xz with lz 47 *flip* value = *yes* or *no* = allow or disallow box flips when it becomes highly skewed 48 *alpha* value = strength of Cauchy stress control parameter 49 *continue* value = *yes* or *no* = whether of not to continue from a previous run 50 *fixedpoint* values = x y z 51 x,y,z = perform barostat dilation/contraction around this point (distance units) 52 53Examples 54"""""""" 55 56.. code-block:: LAMMPS 57 58 fix 1 water npt/cauchy temp 300.0 300.0 100.0 iso 0.0 0.0 1000.0 alpha 0.001 59 60Description 61""""""""""" 62 63This command performs time integration on Nose-Hoover style 64non-Hamiltonian equations of motion which are designed to generate 65positions and velocities sampled from the isothermal-isobaric (npt) 66ensembles. This updates the position and velocity for atoms in the 67group each timestep and the box dimensions. 68 69The thermostatting and barostatting is achieved by adding some dynamic 70variables which are coupled to the particle velocities 71(thermostatting) and simulation domain dimensions (barostatting). In 72addition to basic thermostatting and barostatting, this fix can 73also create a chain of thermostats coupled to the particle thermostat, 74and another chain of thermostats coupled to the barostat 75variables. The barostat can be coupled to the overall box volume, or 76to individual dimensions, including the *xy*, *xz* and *yz* tilt 77dimensions. The external pressure of the barostat can be specified as 78either a scalar pressure (isobaric ensemble) or as components of a 79symmetric stress tensor (constant stress ensemble). When used 80correctly, the time-averaged temperature and stress tensor of the 81particles will match the target values specified by Tstart/Tstop and 82Pstart/Pstop. 83 84The equations of motion used are those of Shinoda et al in 85:ref:`(Shinoda) <nc-Shinoda>`, which combine the hydrostatic equations of 86Martyna, Tobias and Klein in :ref:`(Martyna) <nc-Martyna>` with the strain 87energy proposed by Parrinello and Rahman in 88:ref:`(Parrinello) <nc-Parrinello>`. The time integration schemes closely 89follow the time-reversible measure-preserving Verlet and rRESPA 90integrators derived by Tuckerman et al in :ref:`(Tuckerman) <nc-Tuckerman>`. 91 92---------- 93 94The thermostat parameters are specified using the *temp* keyword. 95Other thermostat-related keywords are *tchain*, *tloop* and *drag*, 96which are discussed below. 97 98The thermostat is applied to only the translational degrees of freedom 99for the particles. The translational degrees of freedom can also have 100a bias velocity removed before thermostatting takes place; see the 101description below. The desired temperature at each timestep is a 102ramped value during the run from *Tstart* to *Tstop*\ . The *Tdamp* 103parameter is specified in time units and determines how rapidly the 104temperature is relaxed. For example, a value of 10.0 means to relax 105the temperature in a timespan of (roughly) 10 time units (e.g. :math:`\tau` 106or fs or ps - see the :doc:`units <units>` command). The atoms in the 107fix group are the only ones whose velocities and positions are updated 108by the velocity/position update portion of the integration. 109 110.. note:: 111 112 A Nose-Hoover thermostat will not work well for arbitrary values 113 of *Tdamp*\ . If *Tdamp* is too small, the temperature can fluctuate 114 wildly; if it is too large, the temperature will take a very long time 115 to equilibrate. A good choice for many models is a *Tdamp* of around 116 100 timesteps. Note that this is NOT the same as 100 time units for 117 most :doc:`units <units>` settings. 118 119---------- 120 121The barostat parameters are specified using one or more of the *iso*, 122*aniso*, *tri*, *x*, *y*, *z*, *xy*, *xz*, *yz*, and *couple* keywords. 123These keywords give you the ability to specify all 6 components of an 124external stress tensor, and to couple various of these components 125together so that the dimensions they represent are varied together 126during a constant-pressure simulation. 127 128Other barostat-related keywords are *pchain*, *mtk*, *ploop*, 129*nreset*, *drag*, and *dilate*, which are discussed below. 130 131Orthogonal simulation boxes have 3 adjustable dimensions (x,y,z). 132Triclinic (non-orthogonal) simulation boxes have 6 adjustable 133dimensions (x,y,z,xy,xz,yz). The :doc:`create_box <create_box>`, :doc:`read data <read_data>`, and :doc:`read_restart <read_restart>` commands 134specify whether the simulation box is orthogonal or non-orthogonal 135(triclinic) and explain the meaning of the xy,xz,yz tilt factors. 136 137The target pressures for each of the 6 components of the stress tensor 138can be specified independently via the *x*, *y*, *z*, *xy*, *xz*, *yz* 139keywords, which correspond to the 6 simulation box dimensions. For 140each component, the external pressure or tensor component at each 141timestep is a ramped value during the run from *Pstart* to *Pstop*\ . 142If a target pressure is specified for a component, then the 143corresponding box dimension will change during a simulation. For 144example, if the *y* keyword is used, the y-box length will change. If 145the *xy* keyword is used, the xy tilt factor will change. A box 146dimension will not change if that component is not specified, although 147you have the option to change that dimension via the :doc:`fix deform <fix_deform>` command. 148 149Note that in order to use the *xy*, *xz*, or *yz* keywords, the 150simulation box must be triclinic, even if its initial tilt factors are 1510.0. 152 153For all barostat keywords, the *Pdamp* parameter operates like the 154*Tdamp* parameter, determining the time scale on which pressure is 155relaxed. For example, a value of 10.0 means to relax the pressure in 156a timespan of (roughly) 10 time units (e.g. :math:`\tau` or fs or ps 157- see the :doc:`units <units>` command). 158 159.. note:: 160 161 A Nose-Hoover barostat will not work well for arbitrary values 162 of *Pdamp*\ . If *Pdamp* is too small, the pressure and volume can 163 fluctuate wildly; if it is too large, the pressure will take a very 164 long time to equilibrate. A good choice for many models is a *Pdamp* 165 of around 1000 timesteps. However, note that *Pdamp* is specified in 166 time units, and that timesteps are NOT the same as time units for most 167 :doc:`units <units>` settings. 168 169Regardless of what atoms are in the fix group (the only atoms which 170are time integrated), a global pressure or stress tensor is computed 171for all atoms. Similarly, when the size of the simulation box is 172changed, all atoms are re-scaled to new positions, unless the keyword 173*dilate* is specified with a *dilate-group-ID* for a group that 174represents a subset of the atoms. This can be useful, for example, to 175leave the coordinates of atoms in a solid substrate unchanged and 176controlling the pressure of a surrounding fluid. This option should 177be used with care, since it can be unphysical to dilate some atoms and 178not others, because it can introduce large, instantaneous 179displacements between a pair of atoms (one dilated, one not) that are 180far from the dilation origin. Also note that for atoms not in the fix 181group, a separate time integration fix like :doc:`fix nve <fix_nve>` or 182:doc:`fix nvt <fix_nh>` can be used on them, independent of whether they 183are dilated or not. 184 185---------- 186 187The *couple* keyword allows two or three of the diagonal components of 188the pressure tensor to be "coupled" together. The value specified 189with the keyword determines which are coupled. For example, *xz* 190means the *Pxx* and *Pzz* components of the stress tensor are coupled. 191*Xyz* means all 3 diagonal components are coupled. Coupling means two 192things: the instantaneous stress will be computed as an average of the 193corresponding diagonal components, and the coupled box dimensions will 194be changed together in lockstep, meaning coupled dimensions will be 195dilated or contracted by the same percentage every timestep. The 196*Pstart*, *Pstop*, *Pdamp* parameters for any coupled dimensions must 197be identical. *Couple xyz* can be used for a 2d simulation; the *z* 198dimension is simply ignored. 199 200---------- 201 202The *iso*, *aniso*, and *tri* keywords are simply shortcuts that are 203equivalent to specifying several other keywords together. 204 205The keyword *iso* means couple all 3 diagonal components together when 206pressure is computed (hydrostatic pressure), and dilate/contract the 207dimensions together. Using "iso Pstart Pstop Pdamp" is the same as 208specifying these 4 keywords: 209 210.. parsed-literal:: 211 212 x Pstart Pstop Pdamp 213 y Pstart Pstop Pdamp 214 z Pstart Pstop Pdamp 215 couple xyz 216 217The keyword *aniso* means *x*, *y*, and *z* dimensions are controlled 218independently using the *Pxx*, *Pyy*, and *Pzz* components of the 219stress tensor as the driving forces, and the specified scalar external 220pressure. Using "aniso Pstart Pstop Pdamp" is the same as specifying 221these 4 keywords: 222 223.. parsed-literal:: 224 225 x Pstart Pstop Pdamp 226 y Pstart Pstop Pdamp 227 z Pstart Pstop Pdamp 228 couple none 229 230The keyword *tri* means *x*, *y*, *z*, *xy*, *xz*, and *yz* dimensions 231are controlled independently using their individual stress components 232as the driving forces, and the specified scalar pressure as the 233external normal stress. Using "tri Pstart Pstop Pdamp" is the same as 234specifying these 7 keywords: 235 236.. parsed-literal:: 237 238 x Pstart Pstop Pdamp 239 y Pstart Pstop Pdamp 240 z Pstart Pstop Pdamp 241 xy 0.0 0.0 Pdamp 242 yz 0.0 0.0 Pdamp 243 xz 0.0 0.0 Pdamp 244 couple none 245 246---------- 247 248In some cases (e.g. for solids) the pressure (volume) and/or 249temperature of the system can oscillate undesirably when a Nose/Hoover 250barostat and thermostat is applied. The optional *drag* keyword will 251damp these oscillations, although it alters the Nose/Hoover equations. 252A value of 0.0 (no drag) leaves the Nose/Hoover formalism unchanged. 253A non-zero value adds a drag term; the larger the value specified, the 254greater the damping effect. Performing a short run and monitoring the 255pressure and temperature is the best way to determine if the drag term 256is working. Typically a value between 0.2 to 2.0 is sufficient to 257damp oscillations after a few periods. Note that use of the drag 258keyword will interfere with energy conservation and will also change 259the distribution of positions and velocities so that they do not 260correspond to the nominal NVT, NPT, or NPH ensembles. 261 262An alternative way to control initial oscillations is to use chain 263thermostats. The keyword *tchain* determines the number of thermostats 264in the particle thermostat. A value of 1 corresponds to the original 265Nose-Hoover thermostat. The keyword *pchain* specifies the number of 266thermostats in the chain thermostatting the barostat degrees of 267freedom. A value of 0 corresponds to no thermostatting of the 268barostat variables. 269 270The *mtk* keyword controls whether or not the correction terms due to 271Martyna, Tuckerman, and Klein are included in the equations of motion 272:ref:`(Martyna) <nc-Martyna>`. Specifying *no* reproduces the original 273Hoover barostat, whose volume probability distribution function 274differs from the true NPT and NPH ensembles by a factor of 1/V. Hence 275using *yes* is more correct, but in many cases the difference is 276negligible. 277 278The keyword *tloop* can be used to improve the accuracy of integration 279scheme at little extra cost. The initial and final updates of the 280thermostat variables are broken up into *tloop* sub-steps, each of 281length *dt*\ /\ *tloop*\ . This corresponds to using a first-order 282Suzuki-Yoshida scheme :ref:`(Tuckerman) <nc-Tuckerman>`. The keyword *ploop* 283does the same thing for the barostat thermostat. 284 285The keyword *nreset* controls how often the reference dimensions used 286to define the strain energy are reset. If this keyword is not used, 287or is given a value of zero, then the reference dimensions are set to 288those of the initial simulation domain and are never changed. If the 289simulation domain changes significantly during the simulation, then 290the final average pressure tensor will differ significantly from the 291specified values of the external stress tensor. A value of *nstep* 292means that every *nstep* timesteps, the reference dimensions are set 293to those of the current simulation domain. 294 295The *scaleyz*, *scalexz*, and *scalexy* keywords control whether or 296not the corresponding tilt factors are scaled with the associated box 297dimensions when barostatting triclinic periodic cells. The default 298values *yes* will turn on scaling, which corresponds to adjusting the 299linear dimensions of the cell while preserving its shape. Choosing 300*no* ensures that the tilt factors are not scaled with the box 301dimensions. See below for restrictions and default values in different 302situations. In older versions of LAMMPS, scaling of tilt factors was 303not performed. The old behavior can be recovered by setting all three 304scale keywords to *no*\ . 305 306The *flip* keyword allows the tilt factors for a triclinic box to 307exceed half the distance of the parallel box length, as discussed 308below. If the *flip* value is set to *yes*, the bound is enforced by 309flipping the box when it is exceeded. If the *flip* value is set to 310*no*, the tilt will continue to change without flipping. Note that if 311applied stress induces large deformations (e.g. in a liquid), this 312means the box shape can tilt dramatically and LAMMPS will run less 313efficiently, due to the large volume of communication needed to 314acquire ghost atoms around a processor's irregular-shaped sub-domain. 315For extreme values of tilt, LAMMPS may also lose atoms and generate an 316error. 317 318The *fixedpoint* keyword specifies the fixed point for barostat volume 319changes. By default, it is the center of the box. Whatever point is 320chosen will not move during the simulation. For example, if the lower 321periodic boundaries pass through (0,0,0), and this point is provided 322to *fixedpoint*, then the lower periodic boundaries will remain at 323(0,0,0), while the upper periodic boundaries will move twice as 324far. In all cases, the particle trajectories are unaffected by the 325chosen value, except for a time-dependent constant translation of 326positions. 327 328---------- 329 330.. note:: 331 332 Using a barostat coupled to tilt dimensions *xy*, *xz*, *yz* can 333 sometimes result in arbitrarily large values of the tilt dimensions, 334 i.e. a dramatically deformed simulation box. LAMMPS allows the tilt 335 factors to grow a small amount beyond the normal limit of half the box 336 length (0.6 times the box length), and then performs a box "flip" to 337 an equivalent periodic cell. See the discussion of the *flip* keyword 338 above, to allow this bound to be exceeded, if desired. 339 340The flip operation is described in more detail in the page for 341:doc:`fix deform <fix_deform>`. Both the barostat dynamics and the atom 342trajectories are unaffected by this operation. However, if a tilt 343factor is incremented by a large amount (1.5 times the box length) on 344a single timestep, LAMMPS can not accommodate this event and will 345terminate the simulation with an error. This error typically indicates 346that there is something badly wrong with how the simulation was 347constructed, such as specifying values of *Pstart* that are too far 348from the current stress value, or specifying a timestep that is too 349large. Triclinic barostatting should be used with care. This also is 350true for other barostat styles, although they tend to be more 351forgiving of insults. In particular, it is important to recognize that 352equilibrium liquids can not support a shear stress and that 353equilibrium solids can not support shear stresses that exceed the 354yield stress. 355 356One exception to this rule is if the first dimension in the tilt factor 357(x for xy) is non-periodic. In that case, the limits on the tilt 358factor are not enforced, since flipping the box in that dimension does 359not change the atom positions due to non-periodicity. In this mode, 360if you tilt the system to extreme angles, the simulation will simply 361become inefficient due to the highly skewed simulation box. 362 363.. note:: 364 365 Unlike the :doc:`fix temp/berendsen <fix_temp_berendsen>` command 366 which performs thermostatting but NO time integration, this fix 367 performs thermostatting/barostatting AND time integration. Thus you 368 should not use any other time integration fix, such as :doc:`fix nve <fix_nve>` on atoms to which this fix is applied. Likewise, 369 fix npt/cauchy should not normally be used on atoms that also 370 have their temperature controlled by another fix - e.g. by :doc:`fix langevin <fix_nh>` or :doc:`fix temp/rescale <fix_temp_rescale>` 371 commands. 372 373See the :doc:`Howto thermostat <Howto_thermostat>` and :doc:`Howto barostat <Howto_barostat>` doc pages for a discussion of different 374ways to compute temperature and perform thermostatting and 375barostatting. 376 377---------- 378 379This fix compute a temperature and pressure each timestep. To do 380this, the fix creates its own computes of style "temp" and "pressure", 381as if one of these sets of commands had been issued: 382 383.. code-block:: LAMMPS 384 385 compute fix-ID_temp all temp 386 compute fix-ID_press all pressure fix-ID_temp 387 388The group for both the new temperature and pressure compute is "all" 389since pressure is computed for the entire system. See the :doc:`compute temp <compute_temp>` and :doc:`compute pressure <compute_pressure>` 390commands for details. Note that the IDs of the new computes are the 391fix-ID + underscore + "temp" or fix_ID + underscore + "press". 392 393Note that these are NOT the computes used by thermodynamic output (see 394the :doc:`thermo_style <thermo_style>` command) with ID = *thermo_temp* 395and *thermo_press*. This means you can change the attributes of these 396fix's temperature or pressure via the 397:doc:`compute_modify <compute_modify>` command. Or you can print this 398temperature or pressure during thermodynamic output via the 399:doc:`thermo_style custom <thermo_style>` command using the appropriate 400compute-ID. It also means that changing attributes of *thermo_temp* 401or *thermo_press* will have no effect on this fix. 402 403Like other fixes that perform thermostatting, fix npt/cauchy can 404be used with :doc:`compute commands <compute>` that calculate a 405temperature after removing a "bias" from the atom velocities. 406E.g. removing the center-of-mass velocity from a group of atoms or 407only calculating temperature on the x-component of velocity or only 408calculating temperature for atoms in a geometric region. This is not 409done by default, but only if the :doc:`fix_modify <fix_modify>` command 410is used to assign a temperature compute to this fix that includes such 411a bias term. See the doc pages for individual :doc:`compute commands <compute>` to determine which ones include a bias. In 412this case, the thermostat works in the following manner: the current 413temperature is calculated taking the bias into account, bias is 414removed from each atom, thermostatting is performed on the remaining 415thermal degrees of freedom, and the bias is added back in. 416 417---------- 418 419This fix can be used with either the *verlet* or *respa* 420:doc:`integrators <run_style>`. When using this fix 421with *respa*, LAMMPS uses an integrator constructed 422according to the following factorization of the Liouville propagator 423(for two rRESPA levels): 424 425.. math:: 426 427 \exp \left(\mathrm{i} L \Delta t \right) = & \hat{E} 428 \exp \left(\mathrm{i} L_{\rm T\textrm{-}baro} \frac{\Delta t}{2} \right) 429 \exp \left(\mathrm{i} L_{\rm T\textrm{-}part} \frac{\Delta t}{2} \right) 430 \exp \left(\mathrm{i} L_{\epsilon , 2} \frac{\Delta t}{2} \right) 431 \exp \left(\mathrm{i} L_{2}^{(2)} \frac{\Delta t}{2} \right) \\ 432 &\times \left[ 433 \exp \left(\mathrm{i} L_{2}^{(1)} \frac{\Delta t}{2n} \right) 434 \exp \left(\mathrm{i} L_{\epsilon , 1} \frac{\Delta t}{2n} \right) 435 \exp \left(\mathrm{i} L_1 \frac{\Delta t}{n} \right) 436 \exp \left(\mathrm{i} L_{\epsilon , 1} \frac{\Delta t}{2n} \right) 437 \exp \left(\mathrm{i} L_{2}^{(1)} \frac{\Delta t}{2n} \right) 438 \right]^n \\ 439 &\times 440 \exp \left(\mathrm{i} L_{2}^{(2)} \frac{\Delta t}{2} \right) 441 \exp \left(\mathrm{i} L_{\epsilon , 2} \frac{\Delta t}{2} \right) 442 \exp \left(\mathrm{i} L_{\rm T\textrm{-}part} \frac{\Delta t}{2} \right) 443 \exp \left(\mathrm{i} L_{\rm T\textrm{-}baro} \frac{\Delta t}{2} \right) \\ 444 &+ \mathcal{O} \left(\Delta t^3 \right) 445 446This factorization differs somewhat from that of Tuckerman et al, in 447that the barostat is only updated at the outermost rRESPA level, 448whereas Tuckerman's factorization requires splitting the pressure into 449pieces corresponding to the forces computed at each rRESPA level. In 450theory, the latter method will exhibit better numerical stability. In 451practice, because Pdamp is normally chosen to be a large multiple of 452the outermost rRESPA timestep, the barostat dynamics are not the 453limiting factor for numerical stability. Both factorizations are 454time-reversible and can be shown to preserve the phase space measure 455of the underlying non-Hamiltonian equations of motion. 456 457.. note:: 458 459 Under NPT dynamics, for a system with zero initial total linear 460 momentum, the total momentum fluctuates close to zero. It may occasionally 461 undergo brief excursions to non-negligible values, before returning close 462 to zero. Over long simulations, this has the effect of causing the 463 center-of-mass to undergo a slow random walk. This can be mitigated by 464 resetting the momentum at infrequent intervals using the 465 :doc:`fix momentum <fix_momentum>` command. 466 467---------- 468 469Restart, fix_modify, output, run start/stop, minimize info 470""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" 471 472This fix writes the state of all the thermostat and barostat 473variables to :doc:`binary restart files <restart>`. See the 474:doc:`read_restart <read_restart>` command for info on how to re-specify 475a fix in an input script that reads a restart file, so that the 476operation of the fix continues in an uninterrupted fashion. 477 478The :doc:`fix_modify <fix_modify>` *temp* and *press* options are 479supported by this fix. You can use them to assign a 480:doc:`compute <compute>` you have defined to this fix which will be used 481in its thermostatting or barostatting procedure, as described above. 482If you do this, note that the kinetic energy derived from the compute 483temperature should be consistent with the virial term computed using 484all atoms for the pressure. LAMMPS will warn you if you choose to 485compute temperature on a subset of atoms. 486 487.. note:: 488 489 If both the *temp* and *press* keywords are used in a single 490 thermo_modify command (or in two separate commands), then the order 491 in which the keywords are specified is important. Note that a 492 :doc:`pressure compute <compute_pressure>` defines its own 493 temperature compute as an argument when it is specified. The 494 *temp* keyword will override this (for the pressure compute being 495 used by fix npt), but only if the *temp* keyword comes after the 496 *press* keyword. If the *temp* keyword comes before the *press* 497 keyword, then the new pressure compute specified by the *press* 498 keyword will be unaffected by the *temp* setting. 499 500The cumulative energy change in the system imposed by this fix, due to 501thermostatting and/or barostatting, is included in the 502:doc:`thermodynamic output <thermo_style>` keywords *ecouple* and 503*econserve*. See the :doc:`thermo_style <thermo_style>` page for 504details. 505 506This fix computes a global scalar which can be accessed by various 507:doc:`output commands <Howto_output>`. The scalar is the same 508cumulative energy change due to this fix described in the previous 509paragraph. The scalar value calculated by this fix is "extensive". 510 511This fix also computes a global vector of quantities, which can be 512accessed by various :doc:`output commands <Howto_output>`. Rhe vector 513values are "intensive". 514 515The vector stores internal Nose/Hoover thermostat and barostat 516variables. The number and meaning of the vector values depends on 517which fix is used and the settings for keywords *tchain* and *pchain*, 518which specify the number of Nose/Hoover chains for the thermostat and 519barostat. If no thermostatting is done, then *tchain* is 0. If no 520barostatting is done, then *pchain* is 0. In the following list, 521"ndof" is 0, 1, 3, or 6, and is the number of degrees of freedom in 522the barostat. Its value is 0 if no barostat is used, else its value 523is 6 if any off-diagonal stress tensor component is barostatted, else 524its value is 1 if *couple xyz* is used or *couple xy* for a 2d 525simulation, otherwise its value is 3. 526 527The order of values in the global vector and their meaning is as 528follows. The notation means there are tchain values for eta, followed 529by tchain for eta_dot, followed by ndof for omega, etc: 530 531* eta[tchain] = particle thermostat displacements (unitless) 532* eta_dot[tchain] = particle thermostat velocities (1/time units) 533* omega[ndof] = barostat displacements (unitless) 534* omega_dot[ndof] = barostat velocities (1/time units) 535* etap[pchain] = barostat thermostat displacements (unitless) 536* etap_dot[pchain] = barostat thermostat velocities (1/time units) 537* PE_eta[tchain] = potential energy of each particle thermostat displacement (energy units) 538* KE_eta_dot[tchain] = kinetic energy of each particle thermostat velocity (energy units) 539* PE_omega[ndof] = potential energy of each barostat displacement (energy units) 540* KE_omega_dot[ndof] = kinetic energy of each barostat velocity (energy units) 541* PE_etap[pchain] = potential energy of each barostat thermostat displacement (energy units) 542* KE_etap_dot[pchain] = kinetic energy of each barostat thermostat velocity (energy units) 543* PE_strain[1] = scalar strain energy (energy units) 544 545This fix can ramp its external temperature and pressure over 546multiple runs, using the *start* and *stop* keywords of the 547:doc:`run <run>` command. See the :doc:`run <run>` command for details of 548how to do this. 549 550This fix is not invoked during :doc:`energy minimization <minimize>`. 551 552---------- 553 554Restrictions 555"""""""""""" 556 557This fix is part of the EXTRA-FIX package. It is only enabled if 558LAMMPS was built with that package. See the :doc:`Build package <Build_package>` page for more info. 559 560*X*, *y*, *z* cannot be barostatted if the associated dimension is not 561periodic. *Xy*, *xz*, and *yz* can only be barostatted if the 562simulation domain is triclinic and the second dimension in the keyword 563(\ *y* dimension in *xy*\ ) is periodic. *Z*, *xz*, and *yz*, cannot be 564barostatted for 2D simulations. The :doc:`create_box <create_box>`, 565:doc:`read data <read_data>`, and :doc:`read_restart <read_restart>` 566commands specify whether the simulation box is orthogonal or 567non-orthogonal (triclinic) and explain the meaning of the xy,xz,yz 568tilt factors. 569 570For the *temp* keyword, the final Tstop cannot be 0.0 since it would 571make the external T = 0.0 at some timestep during the simulation which 572is not allowed in the Nose/Hoover formulation. 573 574The *scaleyz yes* and *scalexz yes* keyword/value pairs can not be used 575for 2D simulations. *scaleyz yes*, *scalexz yes*, and *scalexy yes* options 576can only be used if the second dimension in the keyword is periodic, 577and if the tilt factor is not coupled to the barostat via keywords 578*tri*, *yz*, *xz*, and *xy*\ . 579 580The *alpha* keyword modifies the barostat as per Miller et 581al. (Miller)_"#nc-Miller" so that the Cauchy stress is controlled. 582*alpha* is the non-dimensional parameter, typically set to 0.001 or 5830.01 that determines how aggressively the algorithm drives the system 584towards the set Cauchy stresses. Larger values of *alpha* will modify 585the system more quickly, but can lead to instabilities. Smaller 586values will lead to longer convergence time. Since *alpha* also 587influences how much the stress fluctuations deviate from the 588equilibrium fluctuations, it should be set as small as possible. 589 590A *continue* value of *yes* indicates that the fix is subsequent to a 591previous run with the npt/cauchy fix, and the intention is to continue 592from the converged stress state at the end of the previous run. This 593may be required, for example, when implementing a multi-step loading/unloading 594sequence over several fixes. 595 596Setting *alpha* to zero is not permitted. To "turn off" the 597cauchystat control and thus restore the equilibrium stress 598fluctuations, two subsequent fixes should be used. In the first, fix 599npt/cauchy is used and the simulation box equilibrates to the 600correct shape for the desired stresses. In the second, 601:doc:`fix npt <fix_nh>` is used instead which uses the 602original Parrinello-Rahman algorithm, but now with the corrected 603simulation box shape from using fix npt/cauchy. 604 605This fix can be used with dynamic groups as defined by the 606:doc:`group <group>` command. Likewise it can be used with groups to 607which atoms are added or deleted over time, e.g. a deposition 608simulation. However, the conservation properties of the thermostat 609and barostat are defined for systems with a static set of atoms. You 610may observe odd behavior if the atoms in a group vary dramatically 611over time or the atom count becomes very small. 612 613Related commands 614"""""""""""""""" 615 616:doc:`fix nve <fix_nve>`, :doc:`fix_modify <fix_modify>`, 617:doc:`run_style <run_style>` 618 619Default 620""""""" 621 622The keyword defaults are tchain = 3, pchain = 3, mtk = yes, tloop = 623ploop = 1, nreset = 0, drag = 0.0, dilate = all, couple = none, 624cauchystat = no, 625scaleyz = scalexz = scalexy = yes if periodic in second dimension and 626not coupled to barostat, otherwise no. 627 628---------- 629 630.. _nc-Martyna: 631 632**(Martyna)** Martyna, Tobias and Klein, J Chem Phys, 101, 4177 (1994). 633 634.. _nc-Parrinello: 635 636**(Parrinello)** Parrinello and Rahman, J Appl Phys, 52, 7182 (1981). 637 638.. _nc-Tuckerman: 639 640**(Tuckerman)** Tuckerman, Alejandre, Lopez-Rendon, Jochim, and 641Martyna, J Phys A: Math Gen, 39, 5629 (2006). 642 643.. _nc-Shinoda: 644 645**(Shinoda)** Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004). 646 647.. _nc-Miller: 648 649**(Miller)** Miller, Tadmor, Gibson, Bernstein and Pavia, J Chem Phys, 650144, 184107 (2016). 651