diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index bf730eda7c..2a9350640a 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -115,6 +115,7 @@ OPT. * :doc:`npt (iko) ` * :doc:`npt/asphere (o) ` * :doc:`npt/body ` + * :doc:`npt/cauchy ` * :doc:`npt/eff ` * :doc:`npt/sphere (o) ` * :doc:`npt/uef ` diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 916c78b2ee..ae0c3d625a 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -267,6 +267,7 @@ accelerated styles exist. * :doc:`npt ` - constant NPT time integration via Nose/Hoover * :doc:`npt/asphere ` - NPT for aspherical particles * :doc:`npt/body ` - NPT for body particles +* :doc:`npt/cauchy ` - NPT with Cauchy stress * :doc:`npt/eff ` - NPT for nuclei and electrons in the electron force field model * :doc:`npt/sphere ` - NPT for spherical particles * :doc:`npt/uef ` - NPT style time integration with diagonal flow diff --git a/doc/src/fix_nh.rst b/doc/src/fix_nh.rst index 6f3d31b9cf..0f75125058 100644 --- a/doc/src/fix_nh.rst +++ b/doc/src/fix_nh.rst @@ -132,7 +132,7 @@ integrators derived by Tuckerman et al in :ref:`(Tuckerman) `. ---------- -The thermostat parameters for fix styles *nvt* and *npt* is specified +The thermostat parameters for fix styles *nvt* and *npt* are specified using the *temp* keyword. Other thermostat-related keywords are *tchain*\ , *tloop* and *drag*\ , which are discussed below. diff --git a/doc/src/fix_npt_cauchy.rst b/doc/src/fix_npt_cauchy.rst new file mode 100644 index 0000000000..9d0e04b0c0 --- /dev/null +++ b/doc/src/fix_npt_cauchy.rst @@ -0,0 +1,679 @@ +.. index:: fix npt/cauchy + +fix npt/cauchy command +====================== + +Syntax +"""""" + + +.. parsed-literal:: + + fix ID group-ID style_name keyword value ... + +* ID, group-ID are documented in :doc:`fix ` command +* style\_name = *npt/cauchy* +* one or more keyword/value pairs may be appended +* 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* or *update* + + .. parsed-literal:: + + *temp* values = Tstart Tstop Tdamp + Tstart,Tstop = external temperature at start/end of run + Tdamp = temperature damping parameter (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* values = 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 + *drag* value = Df + Df = drag factor added to barostat/thermostat (0.0 = no drag) + *dilate* value = dilate-group-ID + dilate-group-ID = only dilate atoms in this group due to barostat volume changes + *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 + *cauchystat* cauchystat values = alpha continue + alpha = strength of Cauchy stress control parameter + continue = *yes* or *no* = whether of not to continue from a previous run + *fixedpoint* values = x y z + x,y,z = perform barostat dilation/contraction around this point (distance units) + + + +Examples +"""""""" + +fix 1 water npt/cauchy temp 300.0 300.0 100.0 iso 0.0 0.0 1000.0 + +Description +""""""""""" + +This command performs time integration on Nose-Hoover style +non-Hamiltonian equations of motion which are designed to generate +positions and velocities sampled from the isothermal-isobaric (npt) +ensembles. This updates the position and velocity for atoms in the +group each timestep and the box dimensions. + +The thermostatting and barostatting is achieved by adding some dynamic +variables which are coupled to the particle velocities +(thermostatting) and simulation domain dimensions (barostatting). In +addition to basic thermostatting and barostatting, this fix can +also create a chain of thermostats coupled to the particle thermostat, +and another chain of thermostats coupled to the barostat +variables. The barostat can be coupled to the overall box volume, or +to individual dimensions, including the *xy*\ , *xz* and *yz* tilt +dimensions. The external pressure of the barostat can be specified as +either a scalar pressure (isobaric ensemble) or as components of a +symmetric stress tensor (constant stress ensemble). When used +correctly, the time-averaged temperature and stress tensor of the +particles will match the target values specified by Tstart/Tstop and +Pstart/Pstop. + +The equations of motion used are those of Shinoda et al in +:ref:`(Shinoda) `, which combine the hydrostatic equations of +Martyna, Tobias and Klein in :ref:`(Martyna) ` with the strain +energy proposed by Parrinello and Rahman in +:ref:`(Parrinello) `. The time integration schemes closely +follow the time-reversible measure-preserving Verlet and rRESPA +integrators derived by Tuckerman et al in :ref:`(Tuckerman) `. + + +---------- + + +The thermostat parameters are specified using the *temp* keyword. +Other thermostat-related keywords are *tchain*\ , *tloop* and *drag*\ , +which are discussed below. + +The thermostat is applied to only the translational degrees of freedom +for the particles. The translational degrees of freedom can also have +a bias velocity removed before thermostatting takes place; see the +description below. The desired temperature at each timestep 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 +fmsec or psec - see the :doc:`units ` command). The atoms in the +fix group are the only ones whose velocities and positions are updated +by the velocity/position update portion of the integration. + +.. 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. Note that this is NOT the same as 100 time units for + most :doc:`units ` settings. + + +---------- + + +The barostat parameters are 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*\ , *drag*\ , and *dilate*\ , which are discussed below. + +Orthogonal simulation boxes have 3 adjustable dimensions (x,y,z). +Triclinic (non-orthogonal) simulation boxes have 6 adjustable +dimensions (x,y,z,xy,xz,yz). The :doc:`create\_box `, :doc:`read data `, and :doc:`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. + +The target pressures for each of the 6 components of the stress tensor +can be specified independently via the *x*\ , *y*\ , *z*\ , *xy*\ , *xz*\ , *yz* +keywords, which correspond to the 6 simulation box dimensions. For +each component, the external pressure or tensor component at each +timestep is a ramped value during the run from *Pstart* to *Pstop*\ . +If a target pressure is specified for a component, then the +corresponding box dimension will change during a simulation. For +example, if the *y* keyword is used, the y-box length will change. If +the *xy* keyword is used, the xy tilt factor will change. A box +dimension will not change if that component is not specified, although +you have the option to change that dimension via the :doc:`fix deform ` command. + +Note that in order to use the *xy*\ , *xz*\ , or *yz* keywords, the +simulation box must be triclinic, even if its initial tilt factors are +0.0. + +For all barostat keywords, the *Pdamp* parameter operates like the +*Tdamp* parameter, determining the time scale on which pressure is +relaxed. For example, a value of 10.0 means to relax the pressure in +a timespan of (roughly) 10 time units (e.g. tau or fmsec or psec - see +the :doc:`units ` command). + +.. note:: + + A Nose-Hoover barostat will not work well for arbitrary values + of *Pdamp*\ . If *Pdamp* is too small, the pressure and volume can + fluctuate wildly; if it is too large, the pressure will take a very + long time to equilibrate. A good choice for many models is a *Pdamp* + of around 1000 timesteps. However, note that *Pdamp* is specified in + time units, and that timesteps are NOT the same as time units for most + :doc:`units ` settings. + +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, unless the keyword +*dilate* is specified with a *dilate-group-ID* for a group that +represents a subset of the atoms. This can be useful, for example, to +leave the coordinates of atoms in a solid substrate unchanged and +controlling the pressure of a surrounding fluid. This option should +be used with care, since it can be unphysical to dilate some atoms and +not others, because it can introduce large, instantaneous +displacements between a pair of atoms (one dilated, one not) that are +far from the dilation origin. Also note that for atoms not in the fix +group, a separate time integration fix like :doc:`fix nve ` or +:doc:`fix nvt ` can be used on them, independent of whether they +are dilated or not. + + +---------- + + +The *couple* keyword allows two or three of the diagonal components of +the pressure tensor to be "coupled" together. The value specified +with the keyword determines which are coupled. For example, *xz* +means the *Pxx* and *Pzz* components of the stress tensor are coupled. +*Xyz* means all 3 diagonal components are coupled. Coupling means two +things: the instantaneous stress will be computed as an average of the +corresponding diagonal components, and the coupled box dimensions will +be changed together in lockstep, meaning coupled dimensions will be +dilated or contracted by the same percentage every timestep. The +*Pstart*\ , *Pstop*\ , *Pdamp* parameters for any coupled dimensions must +be identical. *Couple xyz* can be used for a 2d simulation; the *z* +dimension is simply ignored. + + +---------- + + +The *iso*\ , *aniso*\ , and *tri* keywords are simply shortcuts that are +equivalent to specifying several other keywords together. + +The keyword *iso* means couple all 3 diagonal components together when +pressure is computed (hydrostatic pressure), and dilate/contract the +dimensions together. Using "iso Pstart Pstop Pdamp" is the same as +specifying these 4 keywords: + + +.. parsed-literal:: + + x Pstart Pstop Pdamp + y Pstart Pstop Pdamp + z Pstart Pstop Pdamp + couple xyz + +The keyword *aniso* means *x*\ , *y*\ , and *z* dimensions are controlled +independently using the *Pxx*\ , *Pyy*\ , and *Pzz* components of the +stress tensor as the driving forces, and the specified scalar external +pressure. Using "aniso Pstart Pstop Pdamp" is the same as specifying +these 4 keywords: + + +.. parsed-literal:: + + x Pstart Pstop Pdamp + y Pstart Pstop Pdamp + z Pstart Pstop Pdamp + couple none + +The keyword *tri* means *x*\ , *y*\ , *z*\ , *xy*\ , *xz*\ , and *yz* dimensions +are controlled independently using their individual stress components +as the driving forces, and the specified scalar pressure as the +external normal stress. Using "tri Pstart Pstop Pdamp" is the same as +specifying these 7 keywords: + + +.. parsed-literal:: + + x Pstart Pstop Pdamp + y Pstart Pstop Pdamp + z Pstart Pstop Pdamp + xy 0.0 0.0 Pdamp + yz 0.0 0.0 Pdamp + xz 0.0 0.0 Pdamp + couple none + + +---------- + + +In some cases (e.g. for solids) the pressure (volume) and/or +temperature of the system can oscillate undesirably when a Nose/Hoover +barostat and thermostat is applied. The optional *drag* keyword will +damp these oscillations, although it alters the Nose/Hoover equations. +A value of 0.0 (no drag) leaves the Nose/Hoover formalism unchanged. +A non-zero value adds a drag term; the larger the value specified, the +greater the damping effect. Performing a short run and monitoring the +pressure and temperature is the best way to determine if the drag term +is working. Typically a value between 0.2 to 2.0 is sufficient to +damp oscillations after a few periods. Note that use of the drag +keyword will interfere with energy conservation and will also change +the distribution of positions and velocities so that they do not +correspond to the nominal NVT, NPT, or NPH ensembles. + +An alternative way to control initial oscillations is to use chain +thermostats. The keyword *tchain* determines the number of thermostats +in the particle thermostat. A value of 1 corresponds to the original +Nose-Hoover thermostat. The keyword *pchain* specifies the number of +thermostats in the chain thermostatting the barostat degrees of +freedom. A value of 0 corresponds to no thermostatting of the +barostat variables. + +The *mtk* keyword controls whether or not the correction terms due to +Martyna, Tuckerman, and Klein are included in the equations of motion +:ref:`(Martyna) `. Specifying *no* reproduces the original +Hoover barostat, whose volume probability distribution function +differs from the true NPT and NPH ensembles by a factor of 1/V. Hence +using *yes* is more correct, but in many cases the difference is +negligible. + +The keyword *tloop* can be used to improve the accuracy of integration +scheme at little extra cost. The initial and final updates of the +thermostat variables are broken up into *tloop* sub-steps, each of +length *dt*\ /\ *tloop*\ . This corresponds to using a first-order +Suzuki-Yoshida scheme :ref:`(Tuckerman) `. The keyword *ploop* +does the same thing for the barostat thermostat. + +The keyword *nreset* controls how often the reference dimensions used +to define the strain energy are reset. If this keyword is not used, +or is given a value of zero, then the reference dimensions are set to +those of the initial simulation domain and are never changed. If the +simulation domain changes significantly during the simulation, then +the final average pressure tensor will differ significantly from the +specified values of the external stress tensor. A value of *nstep* +means that every *nstep* timesteps, the reference dimensions are set +to those of the current simulation domain. + +The *scaleyz*\ , *scalexz*\ , and *scalexy* keywords control whether or +not the corresponding tilt factors are scaled with the associated box +dimensions when barostatting triclinic periodic cells. The default +values *yes* will turn on scaling, which corresponds to adjusting the +linear dimensions of the cell while preserving its shape. Choosing +*no* ensures that the tilt factors are not scaled with the box +dimensions. See below for restrictions and default values in different +situations. In older versions of LAMMPS, scaling of tilt factors was +not performed. The old behavior can be recovered by setting all three +scale keywords to *no*\ . + +The *flip* keyword allows the tilt factors for a triclinic box to +exceed half the distance of the parallel box length, as discussed +below. If the *flip* value is set to *yes*\ , the bound is enforced by +flipping the box when it is exceeded. If the *flip* value is set to +*no*\ , the tilt will continue to change without flipping. Note that if +applied stress induces large deformations (e.g. in a liquid), this +means the box shape can tilt dramatically and LAMMPS will run less +efficiently, due to the large volume of communication needed to +acquire ghost atoms around a processor's irregular-shaped sub-domain. +For extreme values of tilt, LAMMPS may also lose atoms and generate an +error. + +The *fixedpoint* keyword specifies the fixed point for barostat volume +changes. By default, it is the center of the box. Whatever point is +chosen will not move during the simulation. For example, if the lower +periodic boundaries pass through (0,0,0), and this point is provided +to *fixedpoint*\ , then the lower periodic boundaries will remain at +(0,0,0), while the upper periodic boundaries will move twice as +far. In all cases, the particle trajectories are unaffected by the +chosen value, except for a time-dependent constant translation of +positions. + + +---------- + + +.. note:: + + Using a barostat coupled to tilt dimensions *xy*\ , *xz*\ , *yz* can + sometimes result in arbitrarily large values of the tilt dimensions, + i.e. a dramatically deformed simulation box. LAMMPS allows the tilt + factors to grow a small amount beyond the normal limit of half the box + length (0.6 times the box length), and then performs a box "flip" to + an equivalent periodic cell. See the discussion of the *flip* keyword + above, to allow this bound to be exceeded, if desired. + +The flip operation is described in more detail in the doc page for +:doc:`fix deform `. Both the barostat dynamics and the atom +trajectories are unaffected by this operation. However, if a tilt +factor is incremented by a large amount (1.5 times the box length) on +a single timestep, LAMMPS can not accommodate this event and will +terminate the simulation with an error. This error typically indicates +that there is something badly wrong with how the simulation was +constructed, such as specifying values of *Pstart* that are too far +from the current stress value, or specifying a timestep that is too +large. Triclinic barostatting should be used with care. This also is +true for other barostat styles, although they tend to be more +forgiving of insults. In particular, it is important to recognize that +equilibrium liquids can not support a shear stress and that +equilibrium solids can not support shear stresses that exceed the +yield stress. + +One exception to this rule is if the 1st dimension in the tilt factor +(x for xy) is non-periodic. In that case, the limits on the tilt +factor are not enforced, since flipping the box in that dimension does +not change the atom positions due to non-periodicity. In this mode, +if you tilt the system to extreme angles, the simulation will simply +become inefficient due to the highly skewed simulation box. + +.. note:: + + Unlike the :doc:`fix temp/berendsen ` command + which performs thermostatting but NO time integration, this fix + performs thermostatting/barostatting AND time integration. Thus you + should not use any other time integration fix, such as :doc:`fix nve ` on atoms to which this fix is applied. Likewise, + fix npt/cauchy should not normally be used on atoms that also + have their temperature controlled by another fix - e.g. by :doc:`fix langevin ` or :doc:`fix temp/rescale ` + commands. + +See the :doc:`Howto thermostat ` and :doc:`Howto barostat ` doc pages for a discussion of different +ways to compute temperature and perform thermostatting and +barostatting. + + +---------- + + +This fix compute a temperature and pressure each timestep. To do +this, the fix creates its own computes of style "temp" and "pressure", +as if one of these sets of commands had been issued: + + +.. parsed-literal:: + + compute fix-ID_temp all temp + compute fix-ID_press all pressure fix-ID_temp + +The group for both the new temperature and pressure compute is "all" +since pressure is computed for the entire system. See the :doc:`compute temp ` and :doc:`compute pressure ` +commands for details. Note that the IDs of the new computes are the +fix-ID + underscore + "temp" or fix\_ID + underscore + "press". + +Note that these are NOT the computes used by thermodynamic output (see +the :doc:`thermo\_style ` command) with ID = *thermo\_temp* +and *thermo\_press*. This means you can change the attributes of these +fix's temperature or pressure via the +:doc:`compute\_modify ` command. Or you can print this +temperature or pressure during thermodynamic output via the +:doc:`thermo\_style custom ` command using the appropriate +compute-ID. It also means that changing attributes of *thermo\_temp* +or *thermo\_press* will have no effect on this fix. + +Like other fixes that perform thermostatting, fix npt/cauchy can +be used with :doc:`compute commands ` that calculate a +temperature after removing a "bias" from the atom velocities. +E.g. removing the center-of-mass velocity from a group of atoms or +only calculating temperature on the x-component of velocity or only +calculating temperature for atoms in a geometric region. This is not +done by default, but only if the :doc:`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 :doc:`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 degrees of freedom, and the bias is added back in. + + +---------- + + +This fix can be used with either the *verlet* or *respa* +:doc:`integrators `. When using this fix +with *respa*\ , LAMMPS uses an integrator constructed +according to the following factorization of the Liouville propagator +(for two rRESPA levels): + +.. image:: Eqs/fix_nh1.jpg + :align: center + +This factorization differs somewhat from that of Tuckerman et al, in +that the barostat is only updated at the outermost rRESPA level, +whereas Tuckerman's factorization requires splitting the pressure into +pieces corresponding to the forces computed at each rRESPA level. In +theory, the latter method will exhibit better numerical stability. In +practice, because Pdamp is normally chosen to be a large multiple of +the outermost rRESPA timestep, the barostat dynamics are not the +limiting factor for numerical stability. Both factorizations are +time-reversible and can be shown to preserve the phase space measure +of the underlying non-Hamiltonian equations of motion. + +.. note:: + + Under NPT dynamics, for a system with zero initial total linear + momentum, the total momentum fluctuates close to zero. It may occasionally + undergo brief excursions to non-negligible values, before returning close + to zero. Over long simulations, this has the effect of causing the + center-of-mass to undergo a slow random walk. This can be mitigated by + resetting the momentum at infrequent intervals using the + :doc:`fix momentum ` command. + + +---------- + + +**Restart, fix\_modify, output, run start/stop, minimize info:** + +This fix writes the state of all the thermostat and barostat +variables to :doc:`binary restart files `. See the +:doc:`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 :doc:`fix\_modify ` *temp* and *press* options are +supported by this fix. You can use them to assign a +:doc:`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 :doc:`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 :doc:`fix\_modify ` *energy* option is supported by this +fix to add the energy change induced by Nose/Hoover thermostatting +and barostatting to the system's potential energy as part of +:doc:`thermodynamic output `. + +This fix computes a global scalar and a global vector of quantities, +which can be accessed by various :doc:`output commands `. +The scalar value calculated by this fix is "extensive"; the vector +values are "intensive". + +The scalar is the cumulative energy change due to the fix. + +The vector stores internal Nose/Hoover thermostat and barostat +variables. The number and meaning of the vector values depends on +which fix is used and the settings for keywords *tchain* and *pchain*\ , +which specify the number of Nose/Hoover chains for the thermostat and +barostat. If no thermostatting is done, then *tchain* is 0. If no +barostatting is done, then *pchain* is 0. In the following list, +"ndof" is 0, 1, 3, or 6, and is the number of degrees of freedom in +the barostat. Its value is 0 if no barostat is used, else its value +is 6 if any off-diagonal stress tensor component is barostatted, else +its value is 1 if *couple xyz* is used or *couple xy* for a 2d +simulation, otherwise its value is 3. + +The order of values in the global vector and their meaning is as +follows. The notation means there are tchain values for eta, followed +by tchain for eta\_dot, followed by ndof for omega, etc: + +* eta[tchain] = particle thermostat displacements (unitless) +* eta\_dot[tchain] = particle thermostat velocities (1/time units) +* omega[ndof] = barostat displacements (unitless) +* omega\_dot[ndof] = barostat velocities (1/time units) +* etap[pchain] = barostat thermostat displacements (unitless) +* etap\_dot[pchain] = barostat thermostat velocities (1/time units) +* PE\_eta[tchain] = potential energy of each particle thermostat displacement (energy units) +* KE\_eta\_dot[tchain] = kinetic energy of each particle thermostat velocity (energy units) +* PE\_omega[ndof] = potential energy of each barostat displacement (energy units) +* KE\_omega\_dot[ndof] = kinetic energy of each barostat velocity (energy units) +* PE\_etap[pchain] = potential energy of each barostat thermostat displacement (energy units) +* KE\_etap\_dot[pchain] = kinetic energy of each barostat thermostat velocity (energy units) +* PE\_strain[1] = scalar strain energy (energy units) + +This fix can ramp its external temperature and pressure over +multiple runs, using the *start* and *stop* keywords of the +:doc:`run ` command. See the :doc:`run ` command for details of +how to do this. + +This fix is not invoked during :doc:`energy minimization `. + + +---------- + + +Restrictions +"""""""""""" + + +This fix is part of the USER-MISC package. It is only enabled if +LAMMPS was built with that package. See the :doc:`Build package ` doc page for more info. + +*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 2nd dimension in the keyword +(\ *y* dimension in *xy*\ ) is periodic. *Z*\ , *xz*\ , and *yz*\ , cannot be +barostatted for 2D simulations. The :doc:`create\_box `, +:doc:`read data `, and :doc:`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* and *scalexz yes* keyword/value pairs can not be used +for 2D simulations. *scaleyz yes*\ , *scalexz yes*\ , and *scalexy yes* options +can only be used if the 2nd dimension in the keyword is periodic, +and if the tilt factor is not coupled to the barostat via keywords +*tri*\ , *yz*\ , *xz*\ , and *xy*\ . + +Without the *cauchystat* keyword, the barostat algorithm +controls the Second-Piola Kirchhoff stress, which is a stress measure +referred to the unmodified (initial) simulation box. If the box +deforms substantially during the equilibration, the difference between +the set values and the final true (Cauchy) stresses can be +considerable. + +The *cauchystat* keyword modifies the barostat as per Miller et +al. (Miller)\_"#nc-Miller" so that the Cauchy stress is controlled. +*alpha* is the non-dimensional parameter, typically set to 0.001 or +0.01 that determines how aggressively the algorithm drives the system +towards the set Cauchy stresses. Larger values of *alpha* will modify +the system more quickly, but can lead to instabilities. Smaller +values will lead to longer convergence time. Since *alpha* also +influences how much the stress fluctuations deviate from the +equilibrium fluctuations, it should be set as small as possible. + +A *continue* value of *yes* indicates that the fix is subsequent to a +previous run with the npt/cauchy fix, and the intention is to continue +from the converged stress state at the end of the previous run. This +may be required, for example, when implementing a multi-step loading/unloading +sequence over several fixes. + +Setting *alpha* to zero is not permitted. To "turn off" the +cauchystat control and thus restore the equilibrium stress +fluctuations, two subsequent fixes should be used. In the first, the +cauchystat flag is used and the simulation box equilibrates to the +correct shape for the desired stresses. In the second, the *fix* +statement is identical except that the *cauchystat* keyword is removed +(along with related *alpha* and *continue* values). This restores the +original Parrinello-Rahman algorithm, but now with the correct simulation +box shape from the first fix. + +This fix can be used with dynamic groups as defined by the +:doc:`group ` command. Likewise it can be used with groups to +which atoms are added or deleted over time, e.g. a deposition +simulation. However, the conservation properties of the thermostat +and barostat are defined for systems with a static set of atoms. You +may observe odd behavior if the atoms in a group vary dramatically +over time or the atom count becomes very small. + +Related commands +"""""""""""""""" + +:doc:`fix nve `, :doc:`fix\_modify `, +:doc:`run\_style ` + +Default +""""""" + +The keyword defaults are tchain = 3, pchain = 3, mtk = yes, tloop = +ploop = 1, nreset = 0, drag = 0.0, dilate = all, couple = none, +cauchystat = no, +scaleyz = scalexz = scalexy = yes if periodic in 2nd dimension and +not coupled to barostat, otherwise no. + + +---------- + + +.. _nc-Martyna: + + + +**(Martyna)** Martyna, Tobias and Klein, J Chem Phys, 101, 4177 (1994). + +.. _nc-Parrinello: + + + +**(Parrinello)** Parrinello and Rahman, J Appl Phys, 52, 7182 (1981). + +.. _nc-Tuckerman: + + + +**(Tuckerman)** Tuckerman, Alejandre, Lopez-Rendon, Jochim, and +Martyna, J Phys A: Math Gen, 39, 5629 (2006). + +.. _nc-Shinoda: + + + +**(Shinoda)** Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004). + +.. _nc-Miller: + + + +**(Miller)** Miller, Tadmor, Gibson, Bernstein and Pavia, J Chem Phys, +144, 184107 (2016). + + +.. _lws: http://lammps.sandia.gov +.. _ld: Manual.html +.. _lc: Commands_all.html diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 7832abcdc8..c9f1c0cd8e 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -311,6 +311,8 @@ cartesian CasP Caswell Cates +cauchy +cauchystat Cavium Cawkwell cbecker @@ -781,6 +783,7 @@ equi equil equilibrate equilibrated +equilibrates equilibrating equilibration Equilibria @@ -2129,6 +2132,7 @@ pathangle Patomtrans Pattnaik Pavese +Pavia Paxton pbc pc @@ -2188,6 +2192,7 @@ piecewise Pieniazek Pieter pimd +Piola Pisarev Pishevar Pitera diff --git a/examples/USER/misc/cauchy/NiAlH_jea.eam.alloy b/examples/USER/misc/cauchy/NiAlH_jea.eam.alloy new file mode 120000 index 0000000000..9d763bffe4 --- /dev/null +++ b/examples/USER/misc/cauchy/NiAlH_jea.eam.alloy @@ -0,0 +1 @@ +../../../../potentials/NiAlH_jea.eam.alloy \ No newline at end of file diff --git a/examples/USER/misc/cauchy/README b/examples/USER/misc/cauchy/README new file mode 100644 index 0000000000..73b94ec153 --- /dev/null +++ b/examples/USER/misc/cauchy/README @@ -0,0 +1,10 @@ +Run this example by executing: + +% lmp -in in.cauchystat + +Note that this example use an EAM potential, and therefore must be +run with a LAMMPS executable built with the MANYBODY package. + +The first cauchystat fix equilibrates the temperature at zero stress, +the second fix applies a shear stress. Output in avg.txt shows +convergence to correct values. diff --git a/examples/USER/misc/cauchy/in.cauchystat b/examples/USER/misc/cauchy/in.cauchystat new file mode 100644 index 0000000000..169151a2af --- /dev/null +++ b/examples/USER/misc/cauchy/in.cauchystat @@ -0,0 +1,68 @@ +units metal +atom_style atomic +atom_modify map array + +# Box and atom positions: +boundary p p p + +# Defining lattice and creating simulation +# box with atoms inside +lattice fcc 4.05 +region simbox prism 0 5 0 5 0 5 0 0 0 units lattice +create_box 2 simbox +create_atoms 2 box + +# Atomic mass: +mass 1 58.69 +mass 2 26.98154 + +# Potential, Al fcc crystal +pair_style eam/alloy +pair_coeff * * NiAlH_jea.eam.alloy Ni Al +neigh_modify delay 5 + +thermo 100 +thermo_style custom step temp pxx pyy pzz pxy pxz pyz +compute cna all cna/atom 2.8 + +fix 1 all npt/cauchy temp 600.0 600.0 1.0 & + x 0.0 0.0 0.1 & + y 0.0 0.0 0.1 & + z 0.0 0.0 0.1 & + couple none alpha 0.001 continue no + +# dump 1 all cfg 1000 test*.cfg mass type xs ys zs type c_cna + +timestep 0.002 + +variable px equal pxx +variable py equal pyy +variable pz equal pzz +variable sxy equal pxy +variable sxz equal pxz +variable syz equal pyz +variable t equal temp + +fix avg all ave/time 1 100 100 v_t v_px v_py v_pz v_sxy v_sxz v_syz file avg.txt + +variable lx equal lx +variable ly equal ly +variable lz equal ly +variable xy equal xy +variable xz equal xz +variable yz equal yz + +fix box all ave/time 1 100 100 v_lx v_ly v_lz v_xy v_xz v_yz file box.txt + +velocity all create 1200 4928459 rot yes dist gaussian + +run 1000 + +fix 1 all npt/cauchy temp 600.0 600.0 1.0 & + x 0.0 0.0 0.1 & + y 0.0 0.0 0.1 & + z 0.0 0.0 0.1 & + xy -10000.0 -10000.0 0.1 & + couple none alpha 0.001 continue yes + +run 1000 diff --git a/examples/USER/misc/cauchy/log.9Jan2020.cauchystat.g++.1 b/examples/USER/misc/cauchy/log.9Jan2020.cauchystat.g++.1 new file mode 100644 index 0000000000..0ca25b1a95 --- /dev/null +++ b/examples/USER/misc/cauchy/log.9Jan2020.cauchystat.g++.1 @@ -0,0 +1,172 @@ +LAMMPS (09 Jan 2020) +OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:93) + using 1 OpenMP thread(s) per MPI task +units metal +atom_style atomic +atom_modify map array + +# Box and atom positions: +boundary p p p + +# Defining lattice and creating simulation +# box with atoms inside +lattice fcc 4.05 +Lattice spacing in x,y,z = 4.05 4.05 4.05 +region simbox prism 0 5 0 5 0 5 0 0 0 units lattice +create_box 2 simbox +Created triclinic box = (0 0 0) to (20.25 20.25 20.25) with tilt (0 0 0) + 1 by 1 by 1 MPI processor grid +create_atoms 2 box +Created 500 atoms + create_atoms CPU = 0.000987053 secs + +# Atomic mass: +mass 1 58.69 +mass 2 26.98154 + +# Potential, Al fcc crystal +pair_style eam/alloy +pair_coeff * * NiAlH_jea.eam.alloy Ni Al +Reading potential file NiAlH_jea.eam.alloy with DATE: 2007-11-30 +neigh_modify delay 5 + +thermo 100 +thermo_style custom step temp pxx pyy pzz pxy pxz pyz +compute cna all cna/atom 2.8 + +fix 1 all npt/cauchy temp 600.0 600.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1 couple none alpha 0.001 continue no +Using fix npt/cauchy with alpha=0.001000 + this is NOT a continuation run + +# dump 1 all cfg 1000 test*.cfg mass type xs ys zs type c_cna + +timestep 0.002 + +variable px equal pxx +variable py equal pyy +variable pz equal pzz +variable sxy equal pxy +variable sxz equal pxz +variable syz equal pyz +variable t equal temp + +fix avg all ave/time 1 100 100 v_t v_px v_py v_pz v_sxy v_sxz v_syz file avg.txt + +variable lx equal lx +variable ly equal ly +variable lz equal ly +variable xy equal xy +variable xz equal xz +variable yz equal yz + +fix box all ave/time 1 100 100 v_lx v_ly v_lz v_xy v_xz v_yz file box.txt + +velocity all create 1200 4928459 rot yes dist gaussian + +run 1000 +Neighbor list info ... + update every 1 steps, delay 5 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 7.65 + ghost atom cutoff = 7.65 + binsize = 3.825, bins = 6 6 6 + 2 neighbor lists, perpetual/occasional/extra = 1 1 0 + (1) pair eam/alloy, perpetual + attributes: half, newton on + pair build: half/bin/newton/tri + stencil: half/bin/3d/newton/tri + bin: standard + (2) compute cna/atom, occasional + attributes: full, newton on + pair build: full/bin/atomonly + stencil: full/bin/3d + bin: standard +Per MPI rank memory allocation (min/avg/max) = 4.04 | 4.04 | 4.04 Mbytes +Step Temp Pxx Pyy Pzz Pxy Pxz Pyz + 0 1200 9859.2374 9729.7389 10279.526 -110.10907 -391.60768 295.10918 + 100 461.95579 11262.405 9918.4702 7373.1896 1389.9833 -165.54737 -128.04989 + 200 452.7497 4758.0631 6285.2022 9593.9725 389.15901 835.71435 -1853.9679 + 300 451.50974 7980.6036 7524.3514 9584.5276 297.33672 -154.88768 -1927.573 + 400 461.52812 5074.9544 4877.0864 2689.9029 389.66084 224.44814 563.12739 + 500 458.17416 7672.6668 5358.5073 4670.0236 -1251.047 1175.8268 -373.96822 + 600 461.28593 3629.8562 7265.1611 6970.1746 523.3139 1295.8252 -121.17116 + 700 466.86592 5224.2421 4121.434 4368.4226 230.85768 -65.765274 -1271.8354 + 800 491.38828 -233.79818 2799.6028 5023.998 919.08469 -411.66796 422.33219 + 900 473.16465 6486.5426 4028.6955 2503.9771 451.96928 1309.8322 -557.83472 + 1000 472.85932 4303.6923 4674.969 5268.2263 94.551286 1425.2222 -1352.0883 +Loop time of 1.24705 on 1 procs for 1000 steps with 500 atoms + +Performance: 138.567 ns/day, 0.173 hours/ns, 801.892 timesteps/s +94.0% CPU use with 1 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 1.112 | 1.112 | 1.112 | 0.0 | 89.17 +Neigh | 0.063329 | 0.063329 | 0.063329 | 0.0 | 5.08 +Comm | 0.01994 | 0.01994 | 0.01994 | 0.0 | 1.60 +Output | 0.0014246 | 0.0014246 | 0.0014246 | 0.0 | 0.11 +Modify | 0.045429 | 0.045429 | 0.045429 | 0.0 | 3.64 +Other | | 0.004881 | | | 0.39 + +Nlocal: 500 ave 500 max 500 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 2017 ave 2017 max 2017 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 24689 ave 24689 max 24689 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +FullNghs: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 24689 +Ave neighs/atom = 49.378 +Neighbor list builds = 34 +Dangerous builds = 0 + +fix 1 all npt/cauchy temp 600.0 600.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1 xy -10000.0 -10000.0 0.1 couple none alpha 0.001 continue yes +Using fix npt/cauchy with alpha=0.001000 + this is a continuation run + +run 1000 +Per MPI rank memory allocation (min/avg/max) = 4.056 | 4.056 | 4.056 Mbytes +Step Temp Pxx Pyy Pzz Pxy Pxz Pyz + 1000 472.85932 4303.6923 4674.969 5268.2263 94.551286 1425.2222 -1352.0883 + 1100 471.04772 5593.1614 5874.9867 3608.9922 -1861.938 459.86813 -813.36882 + 1200 473.34727 2337.4765 2050.4694 4330.2198 -3590.2198 -1285.2197 748.05137 + 1300 465.46145 4909.5722 2880.9183 4995.0091 -2860.6934 -895.40937 -382.07531 + 1400 508.53262 92.57534 3722.1136 557.50974 -3121.7615 349.6147 194.5089 + 1500 498.34579 -5755.2352 -3798.1466 -1445.2047 -3218.0887 1733.9103 -555.96371 + 1600 546.45882 -257.80132 407.73403 -39.803803 -3578.1137 1438.3526 -1710.3139 + 1700 570.72785 -2951.9658 -622.89945 1138.4113 -4573.7982 -984.65235 2906.3144 + 1800 650.75622 6086.1524 1111.2919 1726.5115 -3504.716 1140.9767 414.81284 + 1900 690.32264 2763.2044 -609.41535 289.85307 -3788.8761 -1306.8569 760.00116 + 2000 724.01451 -675.72484 522.04263 -468.58167 -6603.3906 -1712.7317 47.61212 +Loop time of 1.27211 on 1 procs for 1000 steps with 500 atoms + +Performance: 135.837 ns/day, 0.177 hours/ns, 786.093 timesteps/s +93.9% CPU use with 1 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 1.0529 | 1.0529 | 1.0529 | 0.0 | 82.77 +Neigh | 0.13671 | 0.13671 | 0.13671 | 0.0 | 10.75 +Comm | 0.018747 | 0.018747 | 0.018747 | 0.0 | 1.47 +Output | 0.00075722 | 0.00075722 | 0.00075722 | 0.0 | 0.06 +Modify | 0.057984 | 0.057984 | 0.057984 | 0.0 | 4.56 +Other | | 0.00499 | | | 0.39 + +Nlocal: 500 ave 500 max 500 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 2040 ave 2040 max 2040 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 23757 ave 23757 max 23757 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +FullNghs: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 23757 +Ave neighs/atom = 47.514 +Neighbor list builds = 78 +Dangerous builds = 0 +Total wall time: 0:00:02 diff --git a/src/.gitignore b/src/.gitignore index 74400a44fb..fabd3b4e4b 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -56,6 +56,9 @@ /fix_qeq*.cpp /fix_qeq*.h +/fix_*cauchy.cpp +/fix_*cauchy.h + /compute_test_nbl.cpp /compute_test_nbl.h /pair_multi_lucy.cpp diff --git a/src/USER-MISC/README b/src/USER-MISC/README index e75e110aa2..ab646aab5a 100644 --- a/src/USER-MISC/README +++ b/src/USER-MISC/README @@ -57,6 +57,7 @@ fix gle, Michele Ceriotti (EPFL Lausanne), michele.ceriotti at gmail.com, 24 Nov fix grem, David Stelter, dstelter@bu.edu, 22 Nov 16 fix imd, Axel Kohlmeyer, akohlmey at gmail.com, 9 Nov 2009 fix ipi, Michele Ceriotti (EPFL Lausanne), michele.ceriotti at gmail.com, 24 Nov 2014 +fix npt/cauchy, R. E. Miller (Carleton University), F. Pavia and S. Pattamatta, 12 Jan 2020 fix nvk, Efrem Braun (UC Berkeley), efrem.braun at gmail.com, https://github.com/lammps/lammps/pull/310 fix pimd, Yuxing Peng (U Chicago), yuxing at uchicago.edu, 24 Nov 2014 fix rhok, Ulf Pedersen (Roskilde U), ulf at urp.dk, 25 Sep 2017 diff --git a/src/USER-MISC/fix_npt_cauchy.cpp b/src/USER-MISC/fix_npt_cauchy.cpp new file mode 100644 index 0000000000..1b0deb5ae3 --- /dev/null +++ b/src/USER-MISC/fix_npt_cauchy.cpp @@ -0,0 +1,2822 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: +------------------------------------------------------------------------- */ + +#include +#include +#include +#include "fix_npt_cauchy.h" +#include "math_extra.h" +#include "atom.h" +#include "force.h" +#include "group.h" +#include "comm.h" +#include "neighbor.h" +#include "irregular.h" +#include "modify.h" +#include "fix_deform.h" +#include "fix_store.h" +#include "compute.h" +#include "kspace.h" +#include "update.h" +#include "respa.h" +#include "domain.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +#define DELTAFLIP 0.1 +#define TILTMAX 1.5 + +enum{NOBIAS,BIAS}; +enum{NONE,XYZ,XY,YZ,XZ}; +enum{ISO,ANISO,TRICLINIC}; + +/* ---------------------------------------------------------------------- + NVT,NPH,NPT integrators for improved Nose-Hoover equations of motion + ---------------------------------------------------------------------- */ + +FixNPTCauchy::FixNPTCauchy(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), + rfix(NULL), id_dilate(NULL), irregular(NULL), id_temp(NULL), id_press(NULL), + eta(NULL), eta_dot(NULL), eta_dotdot(NULL), + eta_mass(NULL), etap(NULL), etap_dot(NULL), etap_dotdot(NULL), + etap_mass(NULL), id_store(NULL),init_store(NULL) +{ + if (narg < 4) error->all(FLERR,"Illegal fix npt/cauchy command"); + + dynamic_group_allow = 1; + time_integrate = 1; + scalar_flag = 1; + vector_flag = 1; + global_freq = 1; + extscalar = 1; + extvector = 0; + + // for CauchyStat + + alpha=0.001; + initRUN = 0; + restartPK = 0; + restart_global = 1; + restart_stored = 0; + + // default values + + pcouple = NONE; + drag = 0.0; + allremap = 1; + mtchain = mpchain = 3; + nc_tchain = nc_pchain = 1; + mtk_flag = 1; + deviatoric_flag = 0; + nreset_h0 = 0; + eta_mass_flag = 1; + omega_mass_flag = 0; + etap_mass_flag = 0; + flipflag = 1; + dipole_flag = 0; + dlm_flag = 0; + + tcomputeflag = 0; + pcomputeflag = 0; + + // turn on tilt factor scaling, whenever applicable + + dimension = domain->dimension; + + scaleyz = scalexz = scalexy = 0; + if (domain->yperiodic && domain->xy != 0.0) scalexy = 1; + if (domain->zperiodic && dimension == 3) { + if (domain->yz != 0.0) scaleyz = 1; + if (domain->xz != 0.0) scalexz = 1; + } + + // set fixed-point to default = center of cell + + fixedpoint[0] = 0.5*(domain->boxlo[0]+domain->boxhi[0]); + fixedpoint[1] = 0.5*(domain->boxlo[1]+domain->boxhi[1]); + fixedpoint[2] = 0.5*(domain->boxlo[2]+domain->boxhi[2]); + + // used by FixNVTSllod to preserve non-default value + + mtchain_default_flag = 1; + + tstat_flag = 0; + double t_period = 0.0; + + double p_period[6]; + for (int i = 0; i < 6; i++) { + p_start[i] = p_stop[i] = p_period[i] = p_target[i] = 0.0; + p_flag[i] = 0; + } + + // process keywords + + int iarg = 3; + + while (iarg < narg) { + if (strcmp(arg[iarg],"temp") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + tstat_flag = 1; + t_start = force->numeric(FLERR,arg[iarg+1]); + t_target = t_start; + t_stop = force->numeric(FLERR,arg[iarg+2]); + t_period = force->numeric(FLERR,arg[iarg+3]); + if (t_start <= 0.0 || t_stop <= 0.0) + error->all(FLERR, + "Target temperature for fix npt/cauchy cannot be 0.0"); + iarg += 4; + + } else if (strcmp(arg[iarg],"iso") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + pcouple = XYZ; + p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]); + p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); + p_period[0] = p_period[1] = p_period[2] = + force->numeric(FLERR,arg[iarg+3]); + p_flag[0] = p_flag[1] = p_flag[2] = 1; + if (dimension == 2) { + p_start[2] = p_stop[2] = p_period[2] = 0.0; + p_flag[2] = 0; + } + iarg += 4; + } else if (strcmp(arg[iarg],"aniso") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + pcouple = NONE; + p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]); + p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); + p_period[0] = p_period[1] = p_period[2] = + force->numeric(FLERR,arg[iarg+3]); + p_flag[0] = p_flag[1] = p_flag[2] = 1; + if (dimension == 2) { + p_start[2] = p_stop[2] = p_period[2] = 0.0; + p_flag[2] = 0; + } + iarg += 4; + } else if (strcmp(arg[iarg],"tri") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + pcouple = NONE; + scalexy = scalexz = scaleyz = 0; + p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]); + p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); + p_period[0] = p_period[1] = p_period[2] = + force->numeric(FLERR,arg[iarg+3]); + p_flag[0] = p_flag[1] = p_flag[2] = 1; + p_start[3] = p_start[4] = p_start[5] = 0.0; + p_stop[3] = p_stop[4] = p_stop[5] = 0.0; + p_period[3] = p_period[4] = p_period[5] = + force->numeric(FLERR,arg[iarg+3]); + p_flag[3] = p_flag[4] = p_flag[5] = 1; + if (dimension == 2) { + p_start[2] = p_stop[2] = p_period[2] = 0.0; + p_flag[2] = 0; + p_start[3] = p_stop[3] = p_period[3] = 0.0; + p_flag[3] = 0; + p_start[4] = p_stop[4] = p_period[4] = 0.0; + p_flag[4] = 0; + } + iarg += 4; + } else if (strcmp(arg[iarg],"x") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + p_start[0] = force->numeric(FLERR,arg[iarg+1]); + p_stop[0] = force->numeric(FLERR,arg[iarg+2]); + p_period[0] = force->numeric(FLERR,arg[iarg+3]); + p_flag[0] = 1; + deviatoric_flag = 1; + iarg += 4; + } else if (strcmp(arg[iarg],"y") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + p_start[1] = force->numeric(FLERR,arg[iarg+1]); + p_stop[1] = force->numeric(FLERR,arg[iarg+2]); + p_period[1] = force->numeric(FLERR,arg[iarg+3]); + p_flag[1] = 1; + deviatoric_flag = 1; + iarg += 4; + } else if (strcmp(arg[iarg],"z") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + p_start[2] = force->numeric(FLERR,arg[iarg+1]); + p_stop[2] = force->numeric(FLERR,arg[iarg+2]); + p_period[2] = force->numeric(FLERR,arg[iarg+3]); + p_flag[2] = 1; + deviatoric_flag = 1; + iarg += 4; + if (dimension == 2) + error->all(FLERR,"Invalid fix npt/cauchy command for a 2d simulation"); + + } else if (strcmp(arg[iarg],"yz") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + p_start[3] = force->numeric(FLERR,arg[iarg+1]); + p_stop[3] = force->numeric(FLERR,arg[iarg+2]); + p_period[3] = force->numeric(FLERR,arg[iarg+3]); + p_flag[3] = 1; + deviatoric_flag = 1; + scaleyz = 0; + iarg += 4; + if (dimension == 2) + error->all(FLERR,"Invalid fix npt/cauchy command for a 2d simulation"); + } else if (strcmp(arg[iarg],"xz") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + p_start[4] = force->numeric(FLERR,arg[iarg+1]); + p_stop[4] = force->numeric(FLERR,arg[iarg+2]); + p_period[4] = force->numeric(FLERR,arg[iarg+3]); + p_flag[4] = 1; + deviatoric_flag = 1; + scalexz = 0; + iarg += 4; + if (dimension == 2) + error->all(FLERR,"Invalid fix npt/cauchy command for a 2d simulation"); + } else if (strcmp(arg[iarg],"xy") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + p_start[5] = force->numeric(FLERR,arg[iarg+1]); + p_stop[5] = force->numeric(FLERR,arg[iarg+2]); + p_period[5] = force->numeric(FLERR,arg[iarg+3]); + p_flag[5] = 1; + deviatoric_flag = 1; + scalexy = 0; + iarg += 4; + + } else if (strcmp(arg[iarg],"couple") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + if (strcmp(arg[iarg+1],"xyz") == 0) pcouple = XYZ; + else if (strcmp(arg[iarg+1],"xy") == 0) pcouple = XY; + else if (strcmp(arg[iarg+1],"yz") == 0) pcouple = YZ; + else if (strcmp(arg[iarg+1],"xz") == 0) pcouple = XZ; + else if (strcmp(arg[iarg+1],"none") == 0) pcouple = NONE; + else error->all(FLERR,"Illegal fix npt/cauchy command"); + iarg += 2; + + } else if (strcmp(arg[iarg],"drag") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + drag = force->numeric(FLERR,arg[iarg+1]); + if (drag < 0.0) error->all(FLERR,"Illegal fix npt/cauchy command"); + iarg += 2; + } else if (strcmp(arg[iarg],"dilate") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + if (strcmp(arg[iarg+1],"all") == 0) allremap = 1; + else { + allremap = 0; + delete [] id_dilate; + int n = strlen(arg[iarg+1]) + 1; + id_dilate = new char[n]; + strcpy(id_dilate,arg[iarg+1]); + int idilate = group->find(id_dilate); + if (idilate == -1) + error->all(FLERR,"Fix npt/cauchy dilate group ID does not exist"); + } + iarg += 2; + + } else if (strcmp(arg[iarg],"tchain") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + mtchain = force->inumeric(FLERR,arg[iarg+1]); + // used by FixNVTSllod to preserve non-default value + mtchain_default_flag = 0; + if (mtchain < 1) error->all(FLERR,"Illegal fix npt/cauchy command"); + iarg += 2; + } else if (strcmp(arg[iarg],"pchain") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + mpchain = force->inumeric(FLERR,arg[iarg+1]); + if (mpchain < 0) error->all(FLERR,"Illegal fix npt/cauchy command"); + iarg += 2; + } else if (strcmp(arg[iarg],"mtk") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + if (strcmp(arg[iarg+1],"yes") == 0) mtk_flag = 1; + else if (strcmp(arg[iarg+1],"no") == 0) mtk_flag = 0; + else error->all(FLERR,"Illegal fix npt/cauchy command"); + iarg += 2; + } else if (strcmp(arg[iarg],"tloop") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + nc_tchain = force->inumeric(FLERR,arg[iarg+1]); + if (nc_tchain < 0) error->all(FLERR,"Illegal fix npt/cauchy command"); + iarg += 2; + } else if (strcmp(arg[iarg],"ploop") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + nc_pchain = force->inumeric(FLERR,arg[iarg+1]); + if (nc_pchain < 0) error->all(FLERR,"Illegal fix npt/cauchy command"); + iarg += 2; + } else if (strcmp(arg[iarg],"nreset") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + nreset_h0 = force->inumeric(FLERR,arg[iarg+1]); + if (nreset_h0 < 0) error->all(FLERR,"Illegal fix npt/cauchy command"); + iarg += 2; + } else if (strcmp(arg[iarg],"scalexy") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + if (strcmp(arg[iarg+1],"yes") == 0) scalexy = 1; + else if (strcmp(arg[iarg+1],"no") == 0) scalexy = 0; + else error->all(FLERR,"Illegal fix npt/cauchy command"); + iarg += 2; + } else if (strcmp(arg[iarg],"scalexz") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + if (strcmp(arg[iarg+1],"yes") == 0) scalexz = 1; + else if (strcmp(arg[iarg+1],"no") == 0) scalexz = 0; + else error->all(FLERR,"Illegal fix npt/cauchy command"); + iarg += 2; + } else if (strcmp(arg[iarg],"scaleyz") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + if (strcmp(arg[iarg+1],"yes") == 0) scaleyz = 1; + else if (strcmp(arg[iarg+1],"no") == 0) scaleyz = 0; + else error->all(FLERR,"Illegal fix npt/cauchy command"); + iarg += 2; + } else if (strcmp(arg[iarg],"flip") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + if (strcmp(arg[iarg+1],"yes") == 0) flipflag = 1; + else if (strcmp(arg[iarg+1],"no") == 0) flipflag = 0; + else error->all(FLERR,"Illegal fix npt/cauchy command"); + iarg += 2; + } else if (strcmp(arg[iarg],"update") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + if (strcmp(arg[iarg+1],"dipole") == 0) dipole_flag = 1; + else if (strcmp(arg[iarg+1],"dipole/dlm") == 0) { + dipole_flag = 1; + dlm_flag = 1; + } else error->all(FLERR,"Illegal fix npt/cauchy command"); + iarg += 2; + } else if (strcmp(arg[iarg],"alpha") == 0) { + alpha = force->numeric(FLERR,arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"continue") == 0) { + if (strcmp(arg[iarg+1],"yes") != 0 && strcmp(arg[iarg+1],"no") != 0) + error->all(FLERR,"Illegal cauchystat continue value. " + "Must be 'yes' or 'no'"); + restartPK = !strcmp(arg[iarg+1],"yes"); + iarg += 2; + } else if (strcmp(arg[iarg],"fixedpoint") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command"); + fixedpoint[0] = force->numeric(FLERR,arg[iarg+1]); + fixedpoint[1] = force->numeric(FLERR,arg[iarg+2]); + fixedpoint[2] = force->numeric(FLERR,arg[iarg+3]); + iarg += 4; + + // disc keyword is also parsed in fix/nh/sphere + + } else if (strcmp(arg[iarg],"disc") == 0) { + iarg++; + + // keywords erate, strain, and ext are also parsed in fix/nh/uef + + } else if (strcmp(arg[iarg],"erate") == 0) { + iarg += 3; + } else if (strcmp(arg[iarg],"strain") == 0) { + iarg += 3; + } else if (strcmp(arg[iarg],"ext") == 0) { + iarg += 2; + + } else error->all(FLERR,"Illegal fix npt/cauchy command"); + } + + // error checks + + if (dimension == 2 && (p_flag[2] || p_flag[3] || p_flag[4])) + error->all(FLERR,"Invalid fix npt/cauchy command for a 2d simulation"); + if (dimension == 2 && (pcouple == YZ || pcouple == XZ)) + error->all(FLERR,"Invalid fix npt/cauchy command for a 2d simulation"); + if (dimension == 2 && (scalexz == 1 || scaleyz == 1 )) + error->all(FLERR,"Invalid fix npt/cauchy command for a 2d simulation"); + + if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0)) + error->all(FLERR,"Invalid fix npt/cauchy command pressure settings"); + if (pcouple == XYZ && dimension == 3 && p_flag[2] == 0) + error->all(FLERR,"Invalid fix npt/cauchy command pressure settings"); + if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0)) + error->all(FLERR,"Invalid fix npt/cauchy command pressure settings"); + if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0)) + error->all(FLERR,"Invalid fix npt/cauchy command pressure settings"); + if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0)) + error->all(FLERR,"Invalid fix npt/cauchy command pressure settings"); + + // require periodicity in tensile dimension + + if (p_flag[0] && domain->xperiodic == 0) + error->all(FLERR,"Cannot use fix npt/cauchy on a non-periodic dimension"); + if (p_flag[1] && domain->yperiodic == 0) + error->all(FLERR,"Cannot use fix npt/cauchy on a non-periodic dimension"); + if (p_flag[2] && domain->zperiodic == 0) + error->all(FLERR,"Cannot use fix npt/cauchy on a non-periodic dimension"); + + // require periodicity in 2nd dim of off-diagonal tilt component + + if (p_flag[3] && domain->zperiodic == 0) + error->all(FLERR, + "Cannot use fix npt/cauchy on a 2nd non-periodic dimension"); + if (p_flag[4] && domain->zperiodic == 0) + error->all(FLERR, + "Cannot use fix npt/cauchy on a 2nd non-periodic dimension"); + if (p_flag[5] && domain->yperiodic == 0) + error->all(FLERR, + "Cannot use fix npt/cauchy on a 2nd non-periodic dimension"); + + if (scaleyz == 1 && domain->zperiodic == 0) + error->all(FLERR,"Cannot use fix npt/cauchy " + "with yz scaling when z is non-periodic dimension"); + if (scalexz == 1 && domain->zperiodic == 0) + error->all(FLERR,"Cannot use fix npt/cauchy " + "with xz scaling when z is non-periodic dimension"); + if (scalexy == 1 && domain->yperiodic == 0) + error->all(FLERR,"Cannot use fix npt/cauchy " + "with xy scaling when y is non-periodic dimension"); + + if (p_flag[3] && scaleyz == 1) + error->all(FLERR,"Cannot use fix npt/cauchy with " + "both yz dynamics and yz scaling"); + if (p_flag[4] && scalexz == 1) + error->all(FLERR,"Cannot use fix npt/cauchy with " + "both xz dynamics and xz scaling"); + if (p_flag[5] && scalexy == 1) + error->all(FLERR,"Cannot use fix npt/cauchy with " + "both xy dynamics and xy scaling"); + + if (!domain->triclinic && (p_flag[3] || p_flag[4] || p_flag[5])) + error->all(FLERR,"Can not specify Pxy/Pxz/Pyz in " + "fix npt/cauchy with non-triclinic box"); + + if (pcouple == XYZ && dimension == 3 && + (p_start[0] != p_start[1] || p_start[0] != p_start[2] || + p_stop[0] != p_stop[1] || p_stop[0] != p_stop[2] || + p_period[0] != p_period[1] || p_period[0] != p_period[2])) + error->all(FLERR,"Invalid fix npt/cauchy pressure settings"); + if (pcouple == XYZ && dimension == 2 && + (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] || + p_period[0] != p_period[1])) + error->all(FLERR,"Invalid fix npt/cauchy pressure settings"); + if (pcouple == XY && + (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] || + p_period[0] != p_period[1])) + error->all(FLERR,"Invalid fix npt/cauchy pressure settings"); + if (pcouple == YZ && + (p_start[1] != p_start[2] || p_stop[1] != p_stop[2] || + p_period[1] != p_period[2])) + error->all(FLERR,"Invalid fix npt/cauchy pressure settings"); + if (pcouple == XZ && + (p_start[0] != p_start[2] || p_stop[0] != p_stop[2] || + p_period[0] != p_period[2])) + error->all(FLERR,"Invalid fix npt/cauchy pressure settings"); + + if (dipole_flag) { + if (!atom->sphere_flag) + error->all(FLERR,"Using update dipole flag requires atom style sphere"); + if (!atom->mu_flag) + error->all(FLERR,"Using update dipole flag requires atom attribute mu"); + } + + if ((tstat_flag && t_period <= 0.0) || + (p_flag[0] && p_period[0] <= 0.0) || + (p_flag[1] && p_period[1] <= 0.0) || + (p_flag[2] && p_period[2] <= 0.0) || + (p_flag[3] && p_period[3] <= 0.0) || + (p_flag[4] && p_period[4] <= 0.0) || + (p_flag[5] && p_period[5] <= 0.0)) + error->all(FLERR,"Fix npt/cauchy damping parameters must be > 0.0"); + + // set pstat_flag and box change and restart_pbc variables + + pre_exchange_flag = 0; + pstat_flag = 0; + pstyle = ISO; + + for (int i = 0; i < 6; i++) + if (p_flag[i]) pstat_flag = 1; + + if (pstat_flag) { + if (p_flag[0] || p_flag[1] || p_flag[2]) box_change_size = 1; + if (p_flag[3] || p_flag[4] || p_flag[5]) box_change_shape = 1; + no_change_box = 1; + if (allremap == 0) restart_pbc = 1; + + // pstyle = TRICLINIC if any off-diagonal term is controlled -> 6 dof + // else pstyle = ISO if XYZ coupling or XY coupling in 2d -> 1 dof + // else pstyle = ANISO -> 3 dof + + if (p_flag[3] || p_flag[4] || p_flag[5]) pstyle = TRICLINIC; + else if (pcouple == XYZ || (dimension == 2 && pcouple == XY)) pstyle = ISO; + else pstyle = ANISO; + + // pre_exchange only required if flips can occur due to shape changes + + if (flipflag && (p_flag[3] || p_flag[4] || p_flag[5])) + pre_exchange_flag = 1; + if (flipflag && (domain->yz != 0.0 || domain->xz != 0.0 || + domain->xy != 0.0)) + pre_exchange_flag = 1; + } + + // convert input periods to frequencies + + t_freq = 0.0; + p_freq[0] = p_freq[1] = p_freq[2] = p_freq[3] = p_freq[4] = p_freq[5] = 0.0; + + if (tstat_flag) t_freq = 1.0 / t_period; + if (p_flag[0]) p_freq[0] = 1.0 / p_period[0]; + if (p_flag[1]) p_freq[1] = 1.0 / p_period[1]; + if (p_flag[2]) p_freq[2] = 1.0 / p_period[2]; + if (p_flag[3]) p_freq[3] = 1.0 / p_period[3]; + if (p_flag[4]) p_freq[4] = 1.0 / p_period[4]; + if (p_flag[5]) p_freq[5] = 1.0 / p_period[5]; + + // Nose/Hoover temp and pressure init + + size_vector = 0; + + if (tstat_flag) { + int ich; + eta = new double[mtchain]; + + // add one extra dummy thermostat, set to zero + + eta_dot = new double[mtchain+1]; + eta_dot[mtchain] = 0.0; + eta_dotdot = new double[mtchain]; + for (ich = 0; ich < mtchain; ich++) { + eta[ich] = eta_dot[ich] = eta_dotdot[ich] = 0.0; + } + eta_mass = new double[mtchain]; + size_vector += 2*2*mtchain; + } + + if (pstat_flag) { + omega[0] = omega[1] = omega[2] = 0.0; + omega_dot[0] = omega_dot[1] = omega_dot[2] = 0.0; + omega_mass[0] = omega_mass[1] = omega_mass[2] = 0.0; + omega[3] = omega[4] = omega[5] = 0.0; + omega_dot[3] = omega_dot[4] = omega_dot[5] = 0.0; + omega_mass[3] = omega_mass[4] = omega_mass[5] = 0.0; + if (pstyle == ISO) size_vector += 2*2*1; + else if (pstyle == ANISO) size_vector += 2*2*3; + else if (pstyle == TRICLINIC) size_vector += 2*2*6; + + if (mpchain) { + int ich; + etap = new double[mpchain]; + + // add one extra dummy thermostat, set to zero + + etap_dot = new double[mpchain+1]; + etap_dot[mpchain] = 0.0; + etap_dotdot = new double[mpchain]; + for (ich = 0; ich < mpchain; ich++) { + etap[ich] = etap_dot[ich] = + etap_dotdot[ich] = 0.0; + } + etap_mass = new double[mpchain]; + size_vector += 2*2*mpchain; + } + + if (deviatoric_flag) size_vector += 1; + } + + nrigid = 0; + rfix = NULL; + + if (pre_exchange_flag) irregular = new Irregular(lmp); + else irregular = NULL; + + // initialize vol0,t0 to zero to signal uninitialized + // values then assigned in init(), if necessary + + vol0 = t0 = 0.0; + + if (!tstat_flag) + error->all(FLERR,"Temperature control must be used with fix npt/cauchy"); + if (!pstat_flag) + error->all(FLERR,"Pressure control must be used with fix npt/cauchy"); + + // create a new compute temp style + // id = fix-ID + temp + // compute group = all since pressure is always global (group all) + // and thus its KE/temperature contribution should use group all + + int n = strlen(id) + 6; + id_temp = new char[n]; + strcpy(id_temp,id); + strcat(id_temp,"_temp"); + + char **newarg = new char*[3]; + newarg[0] = id_temp; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "temp"; + + modify->add_compute(3,newarg); + delete [] newarg; + tcomputeflag = 1; + + // create a new compute pressure style + // id = fix-ID + press, compute group = all + // pass id_temp as 4th arg to pressure constructor + + n = strlen(id) + 7; + id_press = new char[n]; + strcpy(id_press,id); + strcat(id_press,"_press"); + + newarg = new char*[4]; + newarg[0] = id_press; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "pressure"; + newarg[3] = id_temp; + modify->add_compute(4,newarg); + delete [] newarg; + pcomputeflag = 1; +} + +/* ---------------------------------------------------------------------- */ + +void FixNPTCauchy::post_constructor() +{ + CauchyStat_init(); +} + +/* ---------------------------------------------------------------------- */ + +FixNPTCauchy::~FixNPTCauchy() +{ + if (copymode) return; + + delete [] id_dilate; + delete [] rfix; + + delete [] id_store; + delete irregular; + + // delete temperature and pressure if fix created them + + if (tcomputeflag) modify->delete_compute(id_temp); + delete [] id_temp; + + if (tstat_flag) { + delete [] eta; + delete [] eta_dot; + delete [] eta_dotdot; + delete [] eta_mass; + } + + if (pstat_flag) { + if (pcomputeflag) modify->delete_compute(id_press); + delete [] id_press; + if (mpchain) { + delete [] etap; + delete [] etap_dot; + delete [] etap_dotdot; + delete [] etap_mass; + } + } +} + +/* ---------------------------------------------------------------------- */ + +int FixNPTCauchy::setmask() +{ + int mask = 0; + mask |= INITIAL_INTEGRATE; + mask |= FINAL_INTEGRATE; + mask |= THERMO_ENERGY; + mask |= INITIAL_INTEGRATE_RESPA; + mask |= FINAL_INTEGRATE_RESPA; + if (pre_exchange_flag) mask |= PRE_EXCHANGE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixNPTCauchy::init() +{ + // recheck that dilate group has not been deleted + + if (allremap == 0) { + int idilate = group->find(id_dilate); + if (idilate == -1) + error->all(FLERR,"Fix npt/cauchy dilate group ID does not exist"); + dilate_group_bit = group->bitmask[idilate]; + } + + // ensure no conflict with fix deform + + if (pstat_flag) + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"deform") == 0) { + int *dimflag = ((FixDeform *) modify->fix[i])->dimflag; + if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) || + (p_flag[2] && dimflag[2]) || (p_flag[3] && dimflag[3]) || + (p_flag[4] && dimflag[4]) || (p_flag[5] && dimflag[5])) + error->all(FLERR,"Cannot use fix npt/cauchy and fix deform " + "on same component of stress tensor"); + } + + // set temperature and pressure ptrs + + int icompute = modify->find_compute(id_temp); + if (icompute < 0) + error->all(FLERR,"Temperature ID for fix npt/cauchy does not exist"); + temperature = modify->compute[icompute]; + + if (temperature->tempbias) which = BIAS; + else which = NOBIAS; + + if (pstat_flag) { + icompute = modify->find_compute(id_press); + if (icompute < 0) + error->all(FLERR,"Pressure ID for fix npt/cauchy does not exist"); + pressure = modify->compute[icompute]; + } + + // set timesteps and frequencies + + dtv = update->dt; + dtf = 0.5 * update->dt * force->ftm2v; + dthalf = 0.5 * update->dt; + dt4 = 0.25 * update->dt; + dt8 = 0.125 * update->dt; + dto = dthalf; + + p_freq_max = 0.0; + if (pstat_flag) { + p_freq_max = MAX(p_freq[0],p_freq[1]); + p_freq_max = MAX(p_freq_max,p_freq[2]); + if (pstyle == TRICLINIC) { + p_freq_max = MAX(p_freq_max,p_freq[3]); + p_freq_max = MAX(p_freq_max,p_freq[4]); + p_freq_max = MAX(p_freq_max,p_freq[5]); + } + pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain); + } + + if (tstat_flag) + tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain); + + // tally the number of dimensions that are barostatted + // set initial volume and reference cell, if not already done + + if (pstat_flag) { + pdim = p_flag[0] + p_flag[1] + p_flag[2]; + if (vol0 == 0.0) { + if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd; + else vol0 = domain->xprd * domain->yprd; + h0_inv[0] = domain->h_inv[0]; + h0_inv[1] = domain->h_inv[1]; + h0_inv[2] = domain->h_inv[2]; + h0_inv[3] = domain->h_inv[3]; + h0_inv[4] = domain->h_inv[4]; + h0_inv[5] = domain->h_inv[5]; + } + } + + boltz = force->boltz; + nktv2p = force->nktv2p; + + if (force->kspace) kspace_flag = 1; + else kspace_flag = 0; + + if (strstr(update->integrate_style,"respa")) { + nlevels_respa = ((Respa *) update->integrate)->nlevels; + step_respa = ((Respa *) update->integrate)->step; + dto = 0.5*step_respa[0]; + } + + // detect if any rigid fixes exist so rigid bodies move when box is remapped + // rfix[] = indices to each fix rigid + + delete [] rfix; + nrigid = 0; + rfix = NULL; + + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->rigid_flag) nrigid++; + if (nrigid) { + rfix = new int[nrigid]; + nrigid = 0; + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i; + } +} + +/* ---------------------------------------------------------------------- + compute T,P before integrator starts +------------------------------------------------------------------------- */ + +void FixNPTCauchy::setup(int /*vflag*/) +{ + // tdof needed by compute_temp_target() + + t_current = temperature->compute_scalar(); + tdof = temperature->dof; + + // t_target is needed by NVT and NPT in compute_scalar() + // If no thermostat or using fix nphug, + // t_target must be defined by other means. + + if (tstat_flag && strstr(style,"nphug") == NULL) { + compute_temp_target(); + } else if (pstat_flag) { + + // t0 = reference temperature for masses + // cannot be done in init() b/c temperature cannot be called there + // is b/c Modify::init() inits computes after fixes due to dof dependence + // guesstimate a unit-dependent t0 if actual T = 0.0 + // if it was read in from a restart file, leave it be + + if (t0 == 0.0) { + t0 = temperature->compute_scalar(); + if (t0 == 0.0) { + if (strcmp(update->unit_style,"lj") == 0) t0 = 1.0; + else t0 = 300.0; + } + } + t_target = t0; + } + + if (pstat_flag) compute_press_target(); + + if (pstat_flag) { + if (pstyle == ISO) pressure->compute_scalar(); + else pressure->compute_vector(); + couple(); + pressure->addstep(update->ntimestep+1); + } + + // masses and initial forces on thermostat variables + + if (tstat_flag) { + eta_mass[0] = tdof * boltz * t_target / (t_freq*t_freq); + for (int ich = 1; ich < mtchain; ich++) + eta_mass[ich] = boltz * t_target / (t_freq*t_freq); + for (int ich = 1; ich < mtchain; ich++) { + eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1] - + boltz * t_target) / eta_mass[ich]; + } + } + + // masses and initial forces on barostat variables + + if (pstat_flag) { + double kt = boltz * t_target; + double nkt = (atom->natoms + 1) * kt; + + for (int i = 0; i < 3; i++) + if (p_flag[i]) + omega_mass[i] = nkt/(p_freq[i]*p_freq[i]); + + if (pstyle == TRICLINIC) { + for (int i = 3; i < 6; i++) + if (p_flag[i]) omega_mass[i] = nkt/(p_freq[i]*p_freq[i]); + } + + // masses and initial forces on barostat thermostat variables + + if (mpchain) { + etap_mass[0] = boltz * t_target / (p_freq_max*p_freq_max); + for (int ich = 1; ich < mpchain; ich++) + etap_mass[ich] = boltz * t_target / (p_freq_max*p_freq_max); + for (int ich = 1; ich < mpchain; ich++) + etap_dotdot[ich] = + (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] - + boltz * t_target) / etap_mass[ich]; + } + } +} + +/* ---------------------------------------------------------------------- + 1st half of Verlet update +------------------------------------------------------------------------- */ + +void FixNPTCauchy::initial_integrate(int /*vflag*/) +{ + // update eta_press_dot + + if (pstat_flag && mpchain) nhc_press_integrate(); + + // update eta_dot + + if (tstat_flag) { + compute_temp_target(); + nhc_temp_integrate(); + } + + // need to recompute pressure to account for change in KE + // t_current is up-to-date, but compute_temperature is not + // compute appropriately coupled elements of mvv_current + + if (pstat_flag) { + if (pstyle == ISO) { + temperature->compute_scalar(); + pressure->compute_scalar(); + } else { + temperature->compute_vector(); + pressure->compute_vector(); + } + couple(); + pressure->addstep(update->ntimestep+1); + } + + if (pstat_flag) { + compute_press_target(); + nh_omega_dot(); + nh_v_press(); + } + + nve_v(); + + // remap simulation box by 1/2 step + + if (pstat_flag) remap(); + + nve_x(); + + // remap simulation box by 1/2 step + // redo KSpace coeffs since volume has changed + + if (pstat_flag) { + remap(); + if (kspace_flag) force->kspace->setup(); + } +} + +/* ---------------------------------------------------------------------- + 2nd half of Verlet update +------------------------------------------------------------------------- */ + +void FixNPTCauchy::final_integrate() +{ + nve_v(); + + // re-compute temp before nh_v_press() + // only needed for temperature computes with BIAS on reneighboring steps: + // b/c some biases store per-atom values (e.g. temp/profile) + // per-atom values are invalid if reneigh/comm occurred + // since temp->compute() in initial_integrate() + + if (which == BIAS && neighbor->ago == 0) + t_current = temperature->compute_scalar(); + + if (pstat_flag) nh_v_press(); + + // compute new T,P after velocities rescaled by nh_v_press() + // compute appropriately coupled elements of mvv_current + + t_current = temperature->compute_scalar(); + tdof = temperature->dof; + + if (pstat_flag) { + if (pstyle == ISO) pressure->compute_scalar(); + else pressure->compute_vector(); + couple(); + pressure->addstep(update->ntimestep+1); + } + + if (pstat_flag) nh_omega_dot(); + + // update eta_dot + // update eta_press_dot + + if (tstat_flag) nhc_temp_integrate(); + if (pstat_flag && mpchain) nhc_press_integrate(); +} + +/* ---------------------------------------------------------------------- */ + +void FixNPTCauchy::initial_integrate_respa(int /*vflag*/, int ilevel, int /*iloop*/) +{ + // set timesteps by level + + dtv = step_respa[ilevel]; + dtf = 0.5 * step_respa[ilevel] * force->ftm2v; + dthalf = 0.5 * step_respa[ilevel]; + + // outermost level - update eta_dot and omega_dot, apply to v + // all other levels - NVE update of v + // x,v updates only performed for atoms in group + + if (ilevel == nlevels_respa-1) { + + // update eta_press_dot + + if (pstat_flag && mpchain) nhc_press_integrate(); + + // update eta_dot + + if (tstat_flag) { + compute_temp_target(); + nhc_temp_integrate(); + } + + // recompute pressure to account for change in KE + // t_current is up-to-date, but compute_temperature is not + // compute appropriately coupled elements of mvv_current + + if (pstat_flag) { + if (pstyle == ISO) { + temperature->compute_scalar(); + pressure->compute_scalar(); + } else { + temperature->compute_vector(); + pressure->compute_vector(); + } + couple(); + pressure->addstep(update->ntimestep+1); + } + + if (pstat_flag) { + compute_press_target(); + nh_omega_dot(); + nh_v_press(); + } + + nve_v(); + + } else nve_v(); + + // innermost level - also update x only for atoms in group + // if barostat, perform 1/2 step remap before and after + + if (ilevel == 0) { + if (pstat_flag) remap(); + nve_x(); + if (pstat_flag) remap(); + } + + // if barostat, redo KSpace coeffs at outermost level, + // since volume has changed + + if (ilevel == nlevels_respa-1 && kspace_flag && pstat_flag) + force->kspace->setup(); +} + +/* ---------------------------------------------------------------------- */ + +void FixNPTCauchy::final_integrate_respa(int ilevel, int /*iloop*/) +{ + // set timesteps by level + + dtf = 0.5 * step_respa[ilevel] * force->ftm2v; + dthalf = 0.5 * step_respa[ilevel]; + + // outermost level - update eta_dot and omega_dot, apply via final_integrate + // all other levels - NVE update of v + + if (ilevel == nlevels_respa-1) final_integrate(); + else nve_v(); +} + +/* ---------------------------------------------------------------------- */ + +void FixNPTCauchy::couple() +{ + double *tensor = pressure->vector; + + if (pstyle == ISO) + p_current[0] = p_current[1] = p_current[2] = pressure->scalar; + else if (pcouple == XYZ) { + double ave = 1.0/3.0 * (tensor[0] + tensor[1] + tensor[2]); + p_current[0] = p_current[1] = p_current[2] = ave; + } else if (pcouple == XY) { + double ave = 0.5 * (tensor[0] + tensor[1]); + p_current[0] = p_current[1] = ave; + p_current[2] = tensor[2]; + } else if (pcouple == YZ) { + double ave = 0.5 * (tensor[1] + tensor[2]); + p_current[1] = p_current[2] = ave; + p_current[0] = tensor[0]; + } else if (pcouple == XZ) { + double ave = 0.5 * (tensor[0] + tensor[2]); + p_current[0] = p_current[2] = ave; + p_current[1] = tensor[1]; + } else { + p_current[0] = tensor[0]; + p_current[1] = tensor[1]; + p_current[2] = tensor[2]; + } + + if (!std::isfinite(p_current[0]) || !std::isfinite(p_current[1]) || !std::isfinite(p_current[2])) + error->all(FLERR,"Non-numeric pressure - simulation unstable"); + + // switch order from xy-xz-yz to Voigt + + if (pstyle == TRICLINIC) { + p_current[3] = tensor[5]; + p_current[4] = tensor[4]; + p_current[5] = tensor[3]; + + if (!std::isfinite(p_current[3]) || !std::isfinite(p_current[4]) || !std::isfinite(p_current[5])) + error->all(FLERR,"Non-numeric pressure - simulation unstable"); + } +} + +/* ---------------------------------------------------------------------- + change box size + remap all atoms or dilate group atoms depending on allremap flag + if rigid bodies exist, scale rigid body centers-of-mass +------------------------------------------------------------------------- */ + +void FixNPTCauchy::remap() +{ + int i; + double oldlo,oldhi; + double expfac; + + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + double *h = domain->h; + + // omega is not used, except for book-keeping + + for (int i = 0; i < 6; i++) omega[i] += dto*omega_dot[i]; + + // convert pertinent atoms and rigid bodies to lamda coords + + if (allremap) domain->x2lamda(nlocal); + else { + for (i = 0; i < nlocal; i++) + if (mask[i] & dilate_group_bit) + domain->x2lamda(x[i],x[i]); + } + + if (nrigid) + for (i = 0; i < nrigid; i++) + modify->fix[rfix[i]]->deform(0); + + // reset global and local box to new size/shape + + // this operation corresponds to applying the + // translate and scale operations + // corresponding to the solution of the following ODE: + // + // h_dot = omega_dot * h + // + // where h_dot, omega_dot and h are all upper-triangular + // 3x3 tensors. In Voigt notation, the elements of the + // RHS product tensor are: + // h_dot = [0*0, 1*1, 2*2, 1*3+3*2, 0*4+5*3+4*2, 0*5+5*1] + // + // Ordering of operations preserves time symmetry. + + double dto2 = dto/2.0; + double dto4 = dto/4.0; + double dto8 = dto/8.0; + + // off-diagonal components, first half + + if (pstyle == TRICLINIC) { + + if (p_flag[4]) { + expfac = exp(dto8*omega_dot[0]); + h[4] *= expfac; + h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]); + h[4] *= expfac; + } + + if (p_flag[3]) { + expfac = exp(dto4*omega_dot[1]); + h[3] *= expfac; + h[3] += dto2*(omega_dot[3]*h[2]); + h[3] *= expfac; + } + + if (p_flag[5]) { + expfac = exp(dto4*omega_dot[0]); + h[5] *= expfac; + h[5] += dto2*(omega_dot[5]*h[1]); + h[5] *= expfac; + } + + if (p_flag[4]) { + expfac = exp(dto8*omega_dot[0]); + h[4] *= expfac; + h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]); + h[4] *= expfac; + } + } + + // scale diagonal components + // scale tilt factors with cell, if set + + if (p_flag[0]) { + oldlo = domain->boxlo[0]; + oldhi = domain->boxhi[0]; + expfac = exp(dto*omega_dot[0]); + domain->boxlo[0] = (oldlo-fixedpoint[0])*expfac + fixedpoint[0]; + domain->boxhi[0] = (oldhi-fixedpoint[0])*expfac + fixedpoint[0]; + } + + if (p_flag[1]) { + oldlo = domain->boxlo[1]; + oldhi = domain->boxhi[1]; + expfac = exp(dto*omega_dot[1]); + domain->boxlo[1] = (oldlo-fixedpoint[1])*expfac + fixedpoint[1]; + domain->boxhi[1] = (oldhi-fixedpoint[1])*expfac + fixedpoint[1]; + if (scalexy) h[5] *= expfac; + } + + if (p_flag[2]) { + oldlo = domain->boxlo[2]; + oldhi = domain->boxhi[2]; + expfac = exp(dto*omega_dot[2]); + domain->boxlo[2] = (oldlo-fixedpoint[2])*expfac + fixedpoint[2]; + domain->boxhi[2] = (oldhi-fixedpoint[2])*expfac + fixedpoint[2]; + if (scalexz) h[4] *= expfac; + if (scaleyz) h[3] *= expfac; + } + + // off-diagonal components, second half + + if (pstyle == TRICLINIC) { + + if (p_flag[4]) { + expfac = exp(dto8*omega_dot[0]); + h[4] *= expfac; + h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]); + h[4] *= expfac; + } + + if (p_flag[3]) { + expfac = exp(dto4*omega_dot[1]); + h[3] *= expfac; + h[3] += dto2*(omega_dot[3]*h[2]); + h[3] *= expfac; + } + + if (p_flag[5]) { + expfac = exp(dto4*omega_dot[0]); + h[5] *= expfac; + h[5] += dto2*(omega_dot[5]*h[1]); + h[5] *= expfac; + } + + if (p_flag[4]) { + expfac = exp(dto8*omega_dot[0]); + h[4] *= expfac; + h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]); + h[4] *= expfac; + } + + } + + domain->yz = h[3]; + domain->xz = h[4]; + domain->xy = h[5]; + + // tilt factor to cell length ratio can not exceed TILTMAX in one step + + if (domain->yz < -TILTMAX*domain->yprd || + domain->yz > TILTMAX*domain->yprd || + domain->xz < -TILTMAX*domain->xprd || + domain->xz > TILTMAX*domain->xprd || + domain->xy < -TILTMAX*domain->xprd || + domain->xy > TILTMAX*domain->xprd) + error->all(FLERR,"Fix npt/cauchy has tilted box too far in one step - " + "periodic cell is too far from equilibrium state"); + + domain->set_global_box(); + domain->set_local_box(); + + // convert pertinent atoms and rigid bodies back to box coords + + if (allremap) domain->lamda2x(nlocal); + else { + for (i = 0; i < nlocal; i++) + if (mask[i] & dilate_group_bit) + domain->lamda2x(x[i],x[i]); + } + + if (nrigid) + for (i = 0; i < nrigid; i++) + modify->fix[rfix[i]]->deform(1); +} + +/* ---------------------------------------------------------------------- + pack entire state of Fix into one write +------------------------------------------------------------------------- */ + +void FixNPTCauchy::write_restart(FILE *fp) +{ + int nsize = size_restart_global(); + + double *list; + memory->create(list,nsize,"nh:list"); + + pack_restart_data(list); + + if (comm->me == 0) { + int size = nsize * sizeof(double); + fwrite(&size,sizeof(int),1,fp); + fwrite(list,sizeof(double),nsize,fp); + } + + memory->destroy(list); +} + +/* ---------------------------------------------------------------------- + calculate the number of data to be packed +------------------------------------------------------------------------- */ + +int FixNPTCauchy::size_restart_global() +{ + int nsize = 2; + if (tstat_flag) nsize += 1 + 2*mtchain; + if (pstat_flag) { + nsize += 16 + 2*mpchain; + if (deviatoric_flag) nsize += 6; + } + + return nsize; +} + +/* ---------------------------------------------------------------------- + pack restart data +------------------------------------------------------------------------- */ + +int FixNPTCauchy::pack_restart_data(double *list) +{ + int n = 0; + + list[n++] = tstat_flag; + if (tstat_flag) { + list[n++] = mtchain; + for (int ich = 0; ich < mtchain; ich++) + list[n++] = eta[ich]; + for (int ich = 0; ich < mtchain; ich++) + list[n++] = eta_dot[ich]; + } + + list[n++] = pstat_flag; + if (pstat_flag) { + list[n++] = omega[0]; + list[n++] = omega[1]; + list[n++] = omega[2]; + list[n++] = omega[3]; + list[n++] = omega[4]; + list[n++] = omega[5]; + list[n++] = omega_dot[0]; + list[n++] = omega_dot[1]; + list[n++] = omega_dot[2]; + list[n++] = omega_dot[3]; + list[n++] = omega_dot[4]; + list[n++] = omega_dot[5]; + list[n++] = vol0; + list[n++] = t0; + list[n++] = mpchain; + if (mpchain) { + for (int ich = 0; ich < mpchain; ich++) + list[n++] = etap[ich]; + for (int ich = 0; ich < mpchain; ich++) + list[n++] = etap_dot[ich]; + } + + list[n++] = deviatoric_flag; + if (deviatoric_flag) { + list[n++] = h0_inv[0]; + list[n++] = h0_inv[1]; + list[n++] = h0_inv[2]; + list[n++] = h0_inv[3]; + list[n++] = h0_inv[4]; + list[n++] = h0_inv[5]; + } + } + + return n; +} + +/* ---------------------------------------------------------------------- + use state info from restart file to restart the Fix +------------------------------------------------------------------------- */ + +void FixNPTCauchy::restart(char *buf) +{ + int n = 0; + double *list = (double *) buf; + int flag = static_cast (list[n++]); + if (flag) { + int m = static_cast (list[n++]); + if (tstat_flag && m == mtchain) { + for (int ich = 0; ich < mtchain; ich++) + eta[ich] = list[n++]; + for (int ich = 0; ich < mtchain; ich++) + eta_dot[ich] = list[n++]; + } else n += 2*m; + } + flag = static_cast (list[n++]); + if (flag) { + omega[0] = list[n++]; + omega[1] = list[n++]; + omega[2] = list[n++]; + omega[3] = list[n++]; + omega[4] = list[n++]; + omega[5] = list[n++]; + omega_dot[0] = list[n++]; + omega_dot[1] = list[n++]; + omega_dot[2] = list[n++]; + omega_dot[3] = list[n++]; + omega_dot[4] = list[n++]; + omega_dot[5] = list[n++]; + vol0 = list[n++]; + t0 = list[n++]; + int m = static_cast (list[n++]); + if (pstat_flag && m == mpchain) { + for (int ich = 0; ich < mpchain; ich++) + etap[ich] = list[n++]; + for (int ich = 0; ich < mpchain; ich++) + etap_dot[ich] = list[n++]; + } else n+=2*m; + flag = static_cast (list[n++]); + if (flag) { + h0_inv[0] = list[n++]; + h0_inv[1] = list[n++]; + h0_inv[2] = list[n++]; + h0_inv[3] = list[n++]; + h0_inv[4] = list[n++]; + h0_inv[5] = list[n++]; + } + } +} + +/* ---------------------------------------------------------------------- */ + +int FixNPTCauchy::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"temp") == 0) { + if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); + if (tcomputeflag) { + modify->delete_compute(id_temp); + tcomputeflag = 0; + } + delete [] id_temp; + int n = strlen(arg[1]) + 1; + id_temp = new char[n]; + strcpy(id_temp,arg[1]); + + int icompute = modify->find_compute(arg[1]); + if (icompute < 0) + error->all(FLERR,"Could not find fix_modify temperature ID"); + temperature = modify->compute[icompute]; + + if (temperature->tempflag == 0) + error->all(FLERR, + "Fix_modify temperature ID does not compute temperature"); + if (temperature->igroup != 0 && comm->me == 0) + error->warning(FLERR,"Temperature for fix modify is not for group all"); + + // reset id_temp of pressure to new temperature ID + + if (pstat_flag) { + icompute = modify->find_compute(id_press); + if (icompute < 0) + error->all(FLERR,"Pressure ID for fix modify does not exist"); + modify->compute[icompute]->reset_extra_compute_fix(id_temp); + } + + return 2; + + } else if (strcmp(arg[0],"press") == 0) { + if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); + if (!pstat_flag) error->all(FLERR,"Illegal fix_modify command"); + if (pcomputeflag) { + modify->delete_compute(id_press); + pcomputeflag = 0; + } + delete [] id_press; + int n = strlen(arg[1]) + 1; + id_press = new char[n]; + strcpy(id_press,arg[1]); + + int icompute = modify->find_compute(arg[1]); + if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID"); + pressure = modify->compute[icompute]; + + if (pressure->pressflag == 0) + error->all(FLERR,"Fix_modify pressure ID does not compute pressure"); + return 2; + } + + return 0; +} + +/* ---------------------------------------------------------------------- */ + +double FixNPTCauchy::compute_scalar() +{ + int i; + double volume; + double energy; + double kt = boltz * t_target; + double lkt_press = 0.0; + int ich; + if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd; + else volume = domain->xprd * domain->yprd; + + energy = 0.0; + + // thermostat chain energy is equivalent to Eq. (2) in + // Martyna, Tuckerman, Tobias, Klein, Mol Phys, 87, 1117 + // Sum(0.5*p_eta_k^2/Q_k,k=1,M) + L*k*T*eta_1 + Sum(k*T*eta_k,k=2,M), + // where L = tdof + // M = mtchain + // p_eta_k = Q_k*eta_dot[k-1] + // Q_1 = L*k*T/t_freq^2 + // Q_k = k*T/t_freq^2, k > 1 + + if (tstat_flag) { + energy += ke_target * eta[0] + 0.5*eta_mass[0]*eta_dot[0]*eta_dot[0]; + for (ich = 1; ich < mtchain; ich++) + energy += kt * eta[ich] + 0.5*eta_mass[ich]*eta_dot[ich]*eta_dot[ich]; + } + + // barostat energy is equivalent to Eq. (8) in + // Martyna, Tuckerman, Tobias, Klein, Mol Phys, 87, 1117 + // Sum(0.5*p_omega^2/W + P*V), + // where N = natoms + // p_omega = W*omega_dot + // W = N*k*T/p_freq^2 + // sum is over barostatted dimensions + + if (pstat_flag) { + for (i = 0; i < 3; i++) { + if (p_flag[i]) { + energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i] + + p_hydro*(volume-vol0) / (pdim*nktv2p); + lkt_press += kt; + } + } + + if (pstyle == TRICLINIC) { + for (i = 3; i < 6; i++) { + if (p_flag[i]) { + energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i]; + lkt_press += kt; + } + } + } + + // extra contributions from thermostat chain for barostat + + if (mpchain) { + energy += lkt_press * etap[0] + 0.5*etap_mass[0]*etap_dot[0]*etap_dot[0]; + for (ich = 1; ich < mpchain; ich++) + energy += kt * etap[ich] + + 0.5*etap_mass[ich]*etap_dot[ich]*etap_dot[ich]; + } + + // extra contribution from strain energy + + if (deviatoric_flag) energy += compute_strain_energy(); + } + + return energy; +} + +/* ---------------------------------------------------------------------- + return a single element of the following vectors, in this order: + eta[tchain], eta_dot[tchain], omega[ndof], omega_dot[ndof] + etap[pchain], etap_dot[pchain], PE_eta[tchain], KE_eta_dot[tchain] + PE_omega[ndof], KE_omega_dot[ndof], PE_etap[pchain], KE_etap_dot[pchain] + PE_strain[1] + if no thermostat exists, related quantities are omitted from the list + if no barostat exists, related quantities are omitted from the list + ndof = 1,3,6 degrees of freedom for pstyle = ISO,ANISO,TRI +------------------------------------------------------------------------- */ + +double FixNPTCauchy::compute_vector(int n) +{ + int ilen; + + if (tstat_flag) { + ilen = mtchain; + if (n < ilen) return eta[n]; + n -= ilen; + ilen = mtchain; + if (n < ilen) return eta_dot[n]; + n -= ilen; + } + + if (pstat_flag) { + if (pstyle == ISO) { + ilen = 1; + if (n < ilen) return omega[n]; + n -= ilen; + } else if (pstyle == ANISO) { + ilen = 3; + if (n < ilen) return omega[n]; + n -= ilen; + } else { + ilen = 6; + if (n < ilen) return omega[n]; + n -= ilen; + } + + if (pstyle == ISO) { + ilen = 1; + if (n < ilen) return omega_dot[n]; + n -= ilen; + } else if (pstyle == ANISO) { + ilen = 3; + if (n < ilen) return omega_dot[n]; + n -= ilen; + } else { + ilen = 6; + if (n < ilen) return omega_dot[n]; + n -= ilen; + } + + if (mpchain) { + ilen = mpchain; + if (n < ilen) return etap[n]; + n -= ilen; + ilen = mpchain; + if (n < ilen) return etap_dot[n]; + n -= ilen; + } + } + + double volume; + double kt = boltz * t_target; + double lkt_press = kt; + int ich; + if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd; + else volume = domain->xprd * domain->yprd; + + if (tstat_flag) { + ilen = mtchain; + if (n < ilen) { + ich = n; + if (ich == 0) + return ke_target * eta[0]; + else + return kt * eta[ich]; + } + n -= ilen; + ilen = mtchain; + if (n < ilen) { + ich = n; + if (ich == 0) + return 0.5*eta_mass[0]*eta_dot[0]*eta_dot[0]; + else + return 0.5*eta_mass[ich]*eta_dot[ich]*eta_dot[ich]; + } + n -= ilen; + } + + if (pstat_flag) { + if (pstyle == ISO) { + ilen = 1; + if (n < ilen) + return p_hydro*(volume-vol0) / nktv2p; + n -= ilen; + } else if (pstyle == ANISO) { + ilen = 3; + if (n < ilen) { + if (p_flag[n]) + return p_hydro*(volume-vol0) / (pdim*nktv2p); + else + return 0.0; + } + n -= ilen; + } else { + ilen = 6; + if (n < ilen) { + if (n > 2) return 0.0; + else if (p_flag[n]) + return p_hydro*(volume-vol0) / (pdim*nktv2p); + else + return 0.0; + } + n -= ilen; + } + + if (pstyle == ISO) { + ilen = 1; + if (n < ilen) + return pdim*0.5*omega_dot[n]*omega_dot[n]*omega_mass[n]; + n -= ilen; + } else if (pstyle == ANISO) { + ilen = 3; + if (n < ilen) { + if (p_flag[n]) + return 0.5*omega_dot[n]*omega_dot[n]*omega_mass[n]; + else return 0.0; + } + n -= ilen; + } else { + ilen = 6; + if (n < ilen) { + if (p_flag[n]) + return 0.5*omega_dot[n]*omega_dot[n]*omega_mass[n]; + else return 0.0; + } + n -= ilen; + } + + if (mpchain) { + ilen = mpchain; + if (n < ilen) { + ich = n; + if (ich == 0) return lkt_press * etap[0]; + else return kt * etap[ich]; + } + n -= ilen; + ilen = mpchain; + if (n < ilen) { + ich = n; + if (ich == 0) + return 0.5*etap_mass[0]*etap_dot[0]*etap_dot[0]; + else + return 0.5*etap_mass[ich]*etap_dot[ich]*etap_dot[ich]; + } + n -= ilen; + } + + if (deviatoric_flag) { + ilen = 1; + if (n < ilen) + return compute_strain_energy(); + n -= ilen; + } + } + + return 0.0; +} + +/* ---------------------------------------------------------------------- */ + +void FixNPTCauchy::reset_target(double t_new) +{ + t_target = t_start = t_stop = t_new; +} + +/* ---------------------------------------------------------------------- */ + +void FixNPTCauchy::reset_dt() +{ + dtv = update->dt; + dtf = 0.5 * update->dt * force->ftm2v; + dthalf = 0.5 * update->dt; + dt4 = 0.25 * update->dt; + dt8 = 0.125 * update->dt; + dto = dthalf; + + // If using respa, then remap is performed in innermost level + + if (strstr(update->integrate_style,"respa")) + dto = 0.5*step_respa[0]; + + if (pstat_flag) + pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain); + + if (tstat_flag) + tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain); +} + +/* ---------------------------------------------------------------------- + extract thermostat properties +------------------------------------------------------------------------- */ + +void *FixNPTCauchy::extract(const char *str, int &dim) +{ + dim=0; + if (tstat_flag && strcmp(str,"t_target") == 0) { + return &t_target; + } else if (tstat_flag && strcmp(str,"t_start") == 0) { + return &t_start; + } else if (tstat_flag && strcmp(str,"t_stop") == 0) { + return &t_stop; + } else if (tstat_flag && strcmp(str,"mtchain") == 0) { + return &mtchain; + } else if (pstat_flag && strcmp(str,"mpchain") == 0) { + return &mtchain; + } + dim=1; + if (tstat_flag && strcmp(str,"eta") == 0) { + return η + } else if (pstat_flag && strcmp(str,"etap") == 0) { + return η + } else if (pstat_flag && strcmp(str,"p_flag") == 0) { + return &p_flag; + } else if (pstat_flag && strcmp(str,"p_start") == 0) { + return &p_start; + } else if (pstat_flag && strcmp(str,"p_stop") == 0) { + return &p_stop; + } else if (pstat_flag && strcmp(str,"p_target") == 0) { + return &p_target; + } + return NULL; +} + +/* ---------------------------------------------------------------------- + perform half-step update of chain thermostat variables +------------------------------------------------------------------------- */ + +void FixNPTCauchy::nhc_temp_integrate() +{ + int ich; + double expfac; + double kecurrent = tdof * boltz * t_current; + + // Update masses, to preserve initial freq, if flag set + + if (eta_mass_flag) { + eta_mass[0] = tdof * boltz * t_target / (t_freq*t_freq); + for (int ich = 1; ich < mtchain; ich++) + eta_mass[ich] = boltz * t_target / (t_freq*t_freq); + } + + if (eta_mass[0] > 0.0) + eta_dotdot[0] = (kecurrent - ke_target)/eta_mass[0]; + else eta_dotdot[0] = 0.0; + + double ncfac = 1.0/nc_tchain; + for (int iloop = 0; iloop < nc_tchain; iloop++) { + + for (ich = mtchain-1; ich > 0; ich--) { + expfac = exp(-ncfac*dt8*eta_dot[ich+1]); + eta_dot[ich] *= expfac; + eta_dot[ich] += eta_dotdot[ich] * ncfac*dt4; + eta_dot[ich] *= tdrag_factor; + eta_dot[ich] *= expfac; + } + + expfac = exp(-ncfac*dt8*eta_dot[1]); + eta_dot[0] *= expfac; + eta_dot[0] += eta_dotdot[0] * ncfac*dt4; + eta_dot[0] *= tdrag_factor; + eta_dot[0] *= expfac; + + factor_eta = exp(-ncfac*dthalf*eta_dot[0]); + nh_v_temp(); + + // rescale temperature due to velocity scaling + // should not be necessary to explicitly recompute the temperature + + t_current *= factor_eta*factor_eta; + kecurrent = tdof * boltz * t_current; + + if (eta_mass[0] > 0.0) + eta_dotdot[0] = (kecurrent - ke_target)/eta_mass[0]; + else eta_dotdot[0] = 0.0; + + for (ich = 0; ich < mtchain; ich++) + eta[ich] += ncfac*dthalf*eta_dot[ich]; + + eta_dot[0] *= expfac; + eta_dot[0] += eta_dotdot[0] * ncfac*dt4; + eta_dot[0] *= expfac; + + for (ich = 1; ich < mtchain; ich++) { + expfac = exp(-ncfac*dt8*eta_dot[ich+1]); + eta_dot[ich] *= expfac; + eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1] + - boltz * t_target)/eta_mass[ich]; + eta_dot[ich] += eta_dotdot[ich] * ncfac*dt4; + eta_dot[ich] *= expfac; + } + } +} + +/* ---------------------------------------------------------------------- + perform half-step update of chain thermostat variables for barostat + scale barostat velocities +------------------------------------------------------------------------- */ + +void FixNPTCauchy::nhc_press_integrate() +{ + int ich,i,pdof; + double expfac,factor_etap,kecurrent; + double kt = boltz * t_target; + double lkt_press; + + // Update masses, to preserve initial freq, if flag set + + if (omega_mass_flag) { + double nkt = (atom->natoms + 1) * kt; + for (int i = 0; i < 3; i++) + if (p_flag[i]) + omega_mass[i] = nkt/(p_freq[i]*p_freq[i]); + + if (pstyle == TRICLINIC) { + for (int i = 3; i < 6; i++) + if (p_flag[i]) omega_mass[i] = nkt/(p_freq[i]*p_freq[i]); + } + } + + if (etap_mass_flag) { + if (mpchain) { + etap_mass[0] = boltz * t_target / (p_freq_max*p_freq_max); + for (int ich = 1; ich < mpchain; ich++) + etap_mass[ich] = boltz * t_target / (p_freq_max*p_freq_max); + for (int ich = 1; ich < mpchain; ich++) + etap_dotdot[ich] = + (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] - + boltz * t_target) / etap_mass[ich]; + } + } + + kecurrent = 0.0; + pdof = 0; + for (i = 0; i < 3; i++) + if (p_flag[i]) { + kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i]; + pdof++; + } + + if (pstyle == TRICLINIC) { + for (i = 3; i < 6; i++) + if (p_flag[i]) { + kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i]; + pdof++; + } + } + + lkt_press = pdof * kt; + etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0]; + + double ncfac = 1.0/nc_pchain; + for (int iloop = 0; iloop < nc_pchain; iloop++) { + + for (ich = mpchain-1; ich > 0; ich--) { + expfac = exp(-ncfac*dt8*etap_dot[ich+1]); + etap_dot[ich] *= expfac; + etap_dot[ich] += etap_dotdot[ich] * ncfac*dt4; + etap_dot[ich] *= pdrag_factor; + etap_dot[ich] *= expfac; + } + + expfac = exp(-ncfac*dt8*etap_dot[1]); + etap_dot[0] *= expfac; + etap_dot[0] += etap_dotdot[0] * ncfac*dt4; + etap_dot[0] *= pdrag_factor; + etap_dot[0] *= expfac; + + for (ich = 0; ich < mpchain; ich++) + etap[ich] += ncfac*dthalf*etap_dot[ich]; + + factor_etap = exp(-ncfac*dthalf*etap_dot[0]); + for (i = 0; i < 3; i++) + if (p_flag[i]) omega_dot[i] *= factor_etap; + + if (pstyle == TRICLINIC) { + for (i = 3; i < 6; i++) + if (p_flag[i]) omega_dot[i] *= factor_etap; + } + + kecurrent = 0.0; + for (i = 0; i < 3; i++) + if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i]; + + if (pstyle == TRICLINIC) { + for (i = 3; i < 6; i++) + if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i]; + } + + etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0]; + + etap_dot[0] *= expfac; + etap_dot[0] += etap_dotdot[0] * ncfac*dt4; + etap_dot[0] *= expfac; + + for (ich = 1; ich < mpchain; ich++) { + expfac = exp(-ncfac*dt8*etap_dot[ich+1]); + etap_dot[ich] *= expfac; + etap_dotdot[ich] = + (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] - boltz*t_target) / + etap_mass[ich]; + etap_dot[ich] += etap_dotdot[ich] * ncfac*dt4; + etap_dot[ich] *= expfac; + } + } +} + +/* ---------------------------------------------------------------------- + perform half-step barostat scaling of velocities +-----------------------------------------------------------------------*/ + +void FixNPTCauchy::nh_v_press() +{ + double factor[3]; + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + factor[0] = exp(-dt4*(omega_dot[0]+mtk_term2)); + factor[1] = exp(-dt4*(omega_dot[1]+mtk_term2)); + factor[2] = exp(-dt4*(omega_dot[2]+mtk_term2)); + + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + v[i][0] *= factor[0]; + v[i][1] *= factor[1]; + v[i][2] *= factor[2]; + if (pstyle == TRICLINIC) { + v[i][0] += -dthalf*(v[i][1]*omega_dot[5] + v[i][2]*omega_dot[4]); + v[i][1] += -dthalf*v[i][2]*omega_dot[3]; + } + v[i][0] *= factor[0]; + v[i][1] *= factor[1]; + v[i][2] *= factor[2]; + } + } + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + v[i][0] *= factor[0]; + v[i][1] *= factor[1]; + v[i][2] *= factor[2]; + if (pstyle == TRICLINIC) { + v[i][0] += -dthalf*(v[i][1]*omega_dot[5] + v[i][2]*omega_dot[4]); + v[i][1] += -dthalf*v[i][2]*omega_dot[3]; + } + v[i][0] *= factor[0]; + v[i][1] *= factor[1]; + v[i][2] *= factor[2]; + temperature->restore_bias(i,v[i]); + } + } + } +} + +/* ---------------------------------------------------------------------- + perform half-step update of velocities +-----------------------------------------------------------------------*/ + +void FixNPTCauchy::nve_v() +{ + double dtfm; + double **v = atom->v; + double **f = atom->f; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; + } + } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; + } + } + } +} + +/* ---------------------------------------------------------------------- + perform full-step update of positions +-----------------------------------------------------------------------*/ + +void FixNPTCauchy::nve_x() +{ + double **x = atom->x; + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + // x update by full step only for atoms in group + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + } + } +} + +/* ---------------------------------------------------------------------- + perform half-step thermostat scaling of velocities +-----------------------------------------------------------------------*/ + +void FixNPTCauchy::nh_v_temp() +{ + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + v[i][0] *= factor_eta; + v[i][1] *= factor_eta; + v[i][2] *= factor_eta; + } + } + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + v[i][0] *= factor_eta; + v[i][1] *= factor_eta; + v[i][2] *= factor_eta; + temperature->restore_bias(i,v[i]); + } + } + } +} + +/* ---------------------------------------------------------------------- + compute sigma tensor + needed whenever p_target or h0_inv changes +-----------------------------------------------------------------------*/ + +void FixNPTCauchy::compute_sigma() +{ + // if nreset_h0 > 0, reset vol0 and h0_inv + // every nreset_h0 timesteps + + if (nreset_h0 > 0) { + int delta = update->ntimestep - update->beginstep; + if (delta % nreset_h0 == 0) { + if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd; + else vol0 = domain->xprd * domain->yprd; + h0_inv[0] = domain->h_inv[0]; + h0_inv[1] = domain->h_inv[1]; + h0_inv[2] = domain->h_inv[2]; + h0_inv[3] = domain->h_inv[3]; + h0_inv[4] = domain->h_inv[4]; + h0_inv[5] = domain->h_inv[5]; + } + } + + // generate upper-triangular half of + // sigma = vol0*h0inv*(p_target-p_hydro)*h0inv^t + // units of sigma are are PV/L^2 e.g. atm.A + // + // [ 0 5 4 ] [ 0 5 4 ] [ 0 5 4 ] [ 0 - - ] + // [ 5 1 3 ] = [ - 1 3 ] [ 5 1 3 ] [ 5 1 - ] + // [ 4 3 2 ] [ - - 2 ] [ 4 3 2 ] [ 4 3 2 ] + + sigma[0] = + vol0*(h0_inv[0]*((p_target[0]-p_hydro)*h0_inv[0] + + p_target[5]*h0_inv[5]+p_target[4]*h0_inv[4]) + + h0_inv[5]*(p_target[5]*h0_inv[0] + + (p_target[1]-p_hydro)*h0_inv[5]+p_target[3]*h0_inv[4]) + + h0_inv[4]*(p_target[4]*h0_inv[0]+p_target[3]*h0_inv[5] + + (p_target[2]-p_hydro)*h0_inv[4])); + sigma[1] = + vol0*(h0_inv[1]*((p_target[1]-p_hydro)*h0_inv[1] + + p_target[3]*h0_inv[3]) + + h0_inv[3]*(p_target[3]*h0_inv[1] + + (p_target[2]-p_hydro)*h0_inv[3])); + sigma[2] = + vol0*(h0_inv[2]*((p_target[2]-p_hydro)*h0_inv[2])); + sigma[3] = + vol0*(h0_inv[1]*(p_target[3]*h0_inv[2]) + + h0_inv[3]*((p_target[2]-p_hydro)*h0_inv[2])); + sigma[4] = + vol0*(h0_inv[0]*(p_target[4]*h0_inv[2]) + + h0_inv[5]*(p_target[3]*h0_inv[2]) + + h0_inv[4]*((p_target[2]-p_hydro)*h0_inv[2])); + sigma[5] = + vol0*(h0_inv[0]*(p_target[5]*h0_inv[1]+p_target[4]*h0_inv[3]) + + h0_inv[5]*((p_target[1]-p_hydro)*h0_inv[1]+p_target[3]*h0_inv[3]) + + h0_inv[4]*(p_target[3]*h0_inv[1]+(p_target[2]-p_hydro)*h0_inv[3])); +} + +/* ---------------------------------------------------------------------- + compute strain energy +-----------------------------------------------------------------------*/ + +double FixNPTCauchy::compute_strain_energy() +{ + // compute strain energy = 0.5*Tr(sigma*h*h^t) in energy units + + double* h = domain->h; + double d0,d1,d2; + + d0 = + sigma[0]*(h[0]*h[0]+h[5]*h[5]+h[4]*h[4]) + + sigma[5]*( h[1]*h[5]+h[3]*h[4]) + + sigma[4]*( h[2]*h[4]); + d1 = + sigma[5]*( h[5]*h[1]+h[4]*h[3]) + + sigma[1]*( h[1]*h[1]+h[3]*h[3]) + + sigma[3]*( h[2]*h[3]); + d2 = + sigma[4]*( h[4]*h[2]) + + sigma[3]*( h[3]*h[2]) + + sigma[2]*( h[2]*h[2]); + + double energy = 0.5*(d0+d1+d2)/nktv2p; + return energy; +} + +/* ---------------------------------------------------------------------- + compute deviatoric barostat force = h*sigma*h^t +-----------------------------------------------------------------------*/ + +void FixNPTCauchy::compute_deviatoric() +{ + // generate upper-triangular part of h*sigma*h^t + // units of fdev are are PV, e.g. atm*A^3 + // [ 0 5 4 ] [ 0 5 4 ] [ 0 5 4 ] [ 0 - - ] + // [ 5 1 3 ] = [ - 1 3 ] [ 5 1 3 ] [ 5 1 - ] + // [ 4 3 2 ] [ - - 2 ] [ 4 3 2 ] [ 4 3 2 ] + + double* h = domain->h; + + fdev[0] = + h[0]*(sigma[0]*h[0]+sigma[5]*h[5]+sigma[4]*h[4]) + + h[5]*(sigma[5]*h[0]+sigma[1]*h[5]+sigma[3]*h[4]) + + h[4]*(sigma[4]*h[0]+sigma[3]*h[5]+sigma[2]*h[4]); + fdev[1] = + h[1]*( sigma[1]*h[1]+sigma[3]*h[3]) + + h[3]*( sigma[3]*h[1]+sigma[2]*h[3]); + fdev[2] = + h[2]*( sigma[2]*h[2]); + fdev[3] = + h[1]*( sigma[3]*h[2]) + + h[3]*( sigma[2]*h[2]); + fdev[4] = + h[0]*( sigma[4]*h[2]) + + h[5]*( sigma[3]*h[2]) + + h[4]*( sigma[2]*h[2]); + fdev[5] = + h[0]*( sigma[5]*h[1]+sigma[4]*h[3]) + + h[5]*( sigma[1]*h[1]+sigma[3]*h[3]) + + h[4]*( sigma[3]*h[1]+sigma[2]*h[3]); +} + +/* ---------------------------------------------------------------------- + compute target temperature and kinetic energy +-----------------------------------------------------------------------*/ + +void FixNPTCauchy::compute_temp_target() +{ + double delta = update->ntimestep - update->beginstep; + if (delta != 0.0) delta /= update->endstep - update->beginstep; + + t_target = t_start + delta * (t_stop-t_start); + ke_target = tdof * boltz * t_target; +} + +/* ---------------------------------------------------------------------- + compute hydrostatic target pressure +-----------------------------------------------------------------------*/ + +void FixNPTCauchy::compute_press_target() +{ + double delta = update->ntimestep - update->beginstep; + if (delta != 0.0) delta /= update->endstep - update->beginstep; + + p_hydro = 0.0; + for (int i = 0; i < 3; i++) + if (p_flag[i]) { + p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]); + p_hydro += p_target[i]; + } + if (pdim > 0) p_hydro /= pdim; + + if (pstyle == TRICLINIC) + for (int i = 3; i < 6; i++) + p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]); + + // CauchyStat: call CauchyStat to modify p_target[i] and p_hydro, + // if CauchyStat enabled and pressure->vector computation has been initiated + + if (initRUN == 1) CauchyStat(); + if (initRUN == 0) { + for (int i=0; i < 6; i++) { + h_old[i]=domain->h[i]; + } + } + + // when run is initialized tensor[] not available (pressure on cell wall) + + initRUN=1; + + // if deviatoric, recompute sigma each time p_target changes + + if (deviatoric_flag) compute_sigma(); +} + +/* ---------------------------------------------------------------------- + update omega_dot, omega +-----------------------------------------------------------------------*/ + +void FixNPTCauchy::nh_omega_dot() +{ + double f_omega,volume; + + if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd; + else volume = domain->xprd*domain->yprd; + + if (deviatoric_flag) compute_deviatoric(); + + mtk_term1 = 0.0; + if (mtk_flag) { + if (pstyle == ISO) { + mtk_term1 = tdof * boltz * t_current; + mtk_term1 /= pdim * atom->natoms; + } else { + double *mvv_current = temperature->vector; + for (int i = 0; i < 3; i++) + if (p_flag[i]) + mtk_term1 += mvv_current[i]; + mtk_term1 /= pdim * atom->natoms; + } + } + + for (int i = 0; i < 3; i++) + if (p_flag[i]) { + f_omega = (p_current[i]-p_hydro)*volume / + (omega_mass[i] * nktv2p) + mtk_term1 / omega_mass[i]; + if (deviatoric_flag) f_omega -= fdev[i]/(omega_mass[i] * nktv2p); + omega_dot[i] += f_omega*dthalf; + omega_dot[i] *= pdrag_factor; + } + + mtk_term2 = 0.0; + if (mtk_flag) { + for (int i = 0; i < 3; i++) + if (p_flag[i]) + mtk_term2 += omega_dot[i]; + if (pdim > 0) mtk_term2 /= pdim * atom->natoms; + } + + if (pstyle == TRICLINIC) { + for (int i = 3; i < 6; i++) { + if (p_flag[i]) { + f_omega = p_current[i]*volume/(omega_mass[i] * nktv2p); + if (deviatoric_flag) + f_omega -= fdev[i]/(omega_mass[i] * nktv2p); + omega_dot[i] += f_omega*dthalf; + omega_dot[i] *= pdrag_factor; + } + } + } +} + +/* ---------------------------------------------------------------------- + if any tilt ratios exceed limits, set flip = 1 and compute new tilt values + do not flip in x or y if non-periodic (can tilt but not flip) + this is b/c the box length would be changed (dramatically) by flip + if yz tilt exceeded, adjust C vector by one B vector + if xz tilt exceeded, adjust C vector by one A vector + if xy tilt exceeded, adjust B vector by one A vector + check yz first since it may change xz, then xz check comes after + if any flip occurs, create new box in domain + image_flip() adjusts image flags due to box shape change induced by flip + remap() puts atoms outside the new box back into the new box + perform irregular on atoms in lamda coords to migrate atoms to new procs + important that image_flip comes before remap, since remap may change + image flags to new values, making eqs in doc of Domain:image_flip incorrect +------------------------------------------------------------------------- */ + +void FixNPTCauchy::pre_exchange() +{ + double xprd = domain->xprd; + double yprd = domain->yprd; + + // flip is only triggered when tilt exceeds 0.5 by DELTAFLIP + // this avoids immediate re-flipping due to tilt oscillations + + double xtiltmax = (0.5+DELTAFLIP)*xprd; + double ytiltmax = (0.5+DELTAFLIP)*yprd; + + int flipxy,flipxz,flipyz; + flipxy = flipxz = flipyz = 0; + + if (domain->yperiodic) { + if (domain->yz < -ytiltmax) { + domain->yz += yprd; + domain->xz += domain->xy; + flipyz = 1; + } else if (domain->yz >= ytiltmax) { + domain->yz -= yprd; + domain->xz -= domain->xy; + flipyz = -1; + } + } + + if (domain->xperiodic) { + if (domain->xz < -xtiltmax) { + domain->xz += xprd; + flipxz = 1; + } else if (domain->xz >= xtiltmax) { + domain->xz -= xprd; + flipxz = -1; + } + if (domain->xy < -xtiltmax) { + domain->xy += xprd; + flipxy = 1; + } else if (domain->xy >= xtiltmax) { + domain->xy -= xprd; + flipxy = -1; + } + } + + int flip = 0; + if (flipxy || flipxz || flipyz) flip = 1; + + if (flip) { + domain->set_global_box(); + domain->set_local_box(); + + domain->image_flip(flipxy,flipxz,flipyz); + + double **x = atom->x; + imageint *image = atom->image; + int nlocal = atom->nlocal; + for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]); + + domain->x2lamda(atom->nlocal); + irregular->migrate_atoms(); + domain->lamda2x(atom->nlocal); + } +} + +/* ---------------------------------------------------------------------- + memory usage of Irregular +------------------------------------------------------------------------- */ + +double FixNPTCauchy::memory_usage() +{ + double bytes = 0.0; + if (irregular) bytes += irregular->memory_usage(); + return bytes; +} + +/* ---------------------------------------------------------------------- + initialize Cauchystat +------------------------------------------------------------------------- */ + +void FixNPTCauchy::CauchyStat_init() +{ + if (comm->me == 0) { + if (screen) { + fprintf(screen,"Using fix npt/cauchy with alpha=%f\n",alpha); + if (restartPK==1) { + fprintf(screen," (this is a continuation run)\n"); + } else { + fprintf(screen," (this is NOT a continuation run)\n"); + } + } + if (logfile) { + fprintf(logfile,"Using fix npt/cauchy with alpha=%f\n",alpha); + if (restartPK==1) { + fprintf(logfile," this is a continuation run\n"); + } else { + fprintf(logfile," this is NOT a continuation run\n"); + } + } + } + + if (!id_store) { + int n = strlen(id) + 14; + id_store = new char[n]; + strcpy(id_store,id); + strcat(id_store,"_FIX_NH_STORE"); + } + restart_stored = modify->find_fix(id_store); + + if (restartPK==1 && restart_stored < 0) + error->all(FLERR,"Illegal npt/cauchy command. Continuation run" + " must follow a previously equilibrated npt/cauchy run"); + + if (alpha<=0.0) + error->all(FLERR,"Illegal fix npt/cauchy command: " + " Alpha cannot be zero or negative."); + + if (restart_stored < 0) { + char **newarg = new char *[6]; + newarg[0] = id_store; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "STORE"; + newarg[3] = (char *) "global"; + newarg[4] = (char *) "1"; + newarg[5] = (char *) "6"; + modify->add_fix(6,newarg); + delete[] newarg; + restart_stored = modify->find_fix(id_store); + } + init_store = (FixStore *)modify->fix[restart_stored]; + + initRUN = 0; + initPK = 1; + +#define H0(row,col) (H0[(row-1)][(col-1)]) +#define invH0(row,col) (invH0[(row-1)][(col-1)]) + + // initialize H0 to current box shape + + double *h = domain->h; + double *h_inv = domain->h_inv; + + H0(1,1)=h[0]; H0(1,2)=h[5]; H0(1,3)=h[4]; + H0(2,1)=0.0; H0(2,2)=h[1]; H0(2,3)=h[3]; + H0(3,1)=0.0; H0(3,2)=0.0; H0(3,3)=h[2]; + + invH0(1,1)=h_inv[0]; invH0(1,2)=h_inv[5]; invH0(1,3)=h_inv[4]; + invH0(2,1)=0.0; invH0(2,2)=h_inv[1]; invH0(2,3)=h_inv[3]; + invH0(3,1)=0.0; invH0(3,2)=0.0; invH0(3,3)=h_inv[2]; + + CSvol0 = abs(MathExtra::det3(H0)); //find reference volume + +#undef H0 +#undef invH0 +} + +/* ---------------------------------------------------------------------- + Cauchystat +------------------------------------------------------------------------- */ + +void FixNPTCauchy::CauchyStat() +{ + double* h = domain->h; // shape matrix in Voigt notation + double* h_rate = domain->h_rate; // rate of box size/shape change in Voigt notation + double H[3][3]; + double Hdot[3][3]; + +#define H(row,col) (H[(row-1)][(col-1)]) +#define Hdot(row,col) (Hdot[(row-1)][(col-1)]) +#define F(row,col) (F[(row-1)][(col-1)]) +#define Fi(row,col) (Fi[(row-1)][(col-1)]) +#define Fdot(row,col) (Fdot[(row-1)][(col-1)]) +#define invH0(row,col) (invH0[(row-1)][(col-1)]) +#define cauchy(row,col) (cauchy[(row-1)][(col-1)]) +#define setcauchy(row,col) (setcauchy[(row-1)][(col-1)]) +#define setPK(row,col) (setPK[(row-1)][(col-1)]) +#define sigmahat(row,col) (sigmahat[(row-1)][(col-1)]) + + H(1,1)=h[0]; H(1,2)=0.0; H(1,3)=0.0; + H(2,1)=0.0; H(2,2)=h[1]; H(2,3)=0.0; + H(3,1)=0.0; H(3,2)=0.0; H(3,3)=h[2]; + + for (int i=0;i<6;i++) { + h_rate[i]=(h[i]-h_old[i])/update->dt; + h_old[i]=h[i]; + } + + Hdot(1,1)=h_rate[0]; Hdot(1,2)=0.0; Hdot(1,3)=0.0; + Hdot(2,1)=0.0; Hdot(2,2)=h_rate[1]; Hdot(2,3)=0.0; + Hdot(3,1)=0.0; Hdot(3,2)=0.0; Hdot(3,3)=h_rate[2]; + + if (domain->triclinic) { + H(1,2)=h[5]; H(1,3)=h[4]; H(2,3)=h[3]; + Hdot(1,2)=h_rate[5]; Hdot(1,3)=h_rate[4]; Hdot(2,3)=h_rate[3]; + } + + double F[3][3]={{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}}; + double Fi[3][3]={{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}}; + double Fdot[3][3]={{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}}; + double J,vol; + + MathExtra::times3(H,invH0,F); + MathExtra::times3(Hdot,invH0,Fdot); + MathExtra::invert3(F,Fi); + J = MathExtra::det3(F); + vol=CSvol0*J; + + double deltat; + deltat = update->dt; //increment of time step + + // Current pressure on the cell boundaries: + //tensor: 0 1 2 3 4 5 + // x y z xy xz yz + + double *tensor = pressure->vector; + + double cauchy[3][3]; + cauchy(1,1)=-tensor[0]; cauchy(1,2)=0.0; cauchy(1,3)=0.0; + cauchy(2,1)=0.0; cauchy(2,2)=-tensor[1]; cauchy(2,3)=0.0; + cauchy(3,1)=0.0; cauchy(3,2)=0.0; cauchy(3,3)=-tensor[2]; + + if (domain->triclinic) { + cauchy(1,2)=-tensor[3]; cauchy(1,3)=-tensor[4]; + cauchy(2,1)=-tensor[3]; cauchy(2,3)=-tensor[5]; + cauchy(3,1)=-tensor[4]; cauchy(3,2)=-tensor[5]; + } + + // target pressure on the cell boundaries: + // p_target: 0 1 2 3 4 5 + // x y z yz xz xy + + double setcauchy[3][3]; + setcauchy(1,1)=-p_target[0]; setcauchy(1,2)=0.0; setcauchy(1,3)=0.0; + setcauchy(2,1)=0.0; setcauchy(2,2)=-p_target[1]; setcauchy(2,3)=0.0; + setcauchy(3,1)=0.0; setcauchy(3,2)=0.0; setcauchy(3,3)=-p_target[2]; + + if (domain->triclinic) { + setcauchy(1,2)=-p_target[5]; setcauchy(1,3)=-p_target[4]; + setcauchy(2,1)=-p_target[5]; setcauchy(2,3)=-p_target[3]; + setcauchy(3,1)=-p_target[4]; setcauchy(3,2)=-p_target[3]; + } + + //initialize: + if (initPK==1) { + if (restartPK==1) { + double *setPKinit = init_store->astore[0]; + setPK(1,1)=setPKinit[0]; setPK(1,2)=setPKinit[1]; setPK(1,3)=setPKinit[2]; + setPK(2,1)=setPKinit[1]; setPK(2,2)=setPKinit[3]; setPK(2,3)=setPKinit[4]; + setPK(3,1)=setPKinit[2]; setPK(3,2)=setPKinit[4]; setPK(3,3)=setPKinit[5]; + }else { + setPK(1,1)=cauchy(1,1); setPK(1,2)=cauchy(1,2); setPK(1,3)=cauchy(1,3); + setPK(2,1)=cauchy(2,1); setPK(2,2)=cauchy(2,2); setPK(2,3)=cauchy(2,3); + setPK(3,1)=cauchy(3,1); setPK(3,2)=cauchy(3,2); setPK(3,3)=cauchy(3,3); + } + initPK=0; + } + + CauchyStat_Step(Fi,Fdot,cauchy,setcauchy,setPK,vol,CSvol0,deltat,alpha); + + // use currentPK as new target: + // p_target: 0 1 2 3 4 5 + // x y z yz xz xy + + p_target[0]=-setPK(1,1); + p_target[1]=-setPK(2,2); + p_target[2]=-setPK(3,3); + if (pstyle == TRICLINIC) { + p_target[3]=-setPK(2,3); + p_target[4]=-setPK(1,3); + p_target[5]=-setPK(1,2); + } + + p_hydro = 0.0; + for (int i = 0; i < 3; i++) + if (p_flag[i]) { + p_hydro += p_target[i]; + } + p_hydro /= pdim; + + // save information for Cauchystat restart + double *setPKinit = init_store->astore[0]; + setPKinit[0] = setcauchy(1,1); + setPKinit[1] = setcauchy(1,2); + setPKinit[2] = setcauchy(1,3); + setPKinit[3] = setcauchy(2,2); + setPKinit[4] = setcauchy(2,3); + setPKinit[5] = setcauchy(3,3); + +#undef H +#undef Hdot +#undef F +#undef Fi +#undef Fdot +#undef invH0 +#undef cauchy +#undef setcauchy +#undef setPK +#undef sigmahat + +} + +/* ---------------------------------------------------------------------- + CauchyStat_Step + + Inputs: + Fi(3,3) : inverse of the deformation gradient + Fdot(3,3) : time derivative of the deformation gradient (velocity) + cauchy(3,3) : current cauchy stress tensor + setcauchy(3,3) : target cauchy stress tensor + alpha : parameter =0.01 + setPK(3,3) : current PK stress tensor, at initialization (time=0) setPK=cauchy + volume : current volume of the periodic box + volume0 : initial volume of the periodic box + deltat : simulation step increment + alpha : parameter + + Outputs: + setPK(3,3) : PK stress tensor at the next time step + ------------------------------------------------------------------------- */ + +void FixNPTCauchy::CauchyStat_Step(double (&Fi)[3][3], double (&Fdot)[3][3], + double (&cauchy)[3][3], double (&setcauchy)[3][3], + double (&setPK)[3][3], double volume, + double volume0, double deltat, double alpha) +{ + + //macros to go from c to fortran style for arrays: +#define F(row,col) (F[(row-1)][(col-1)]) +#define Fi(row,col) (Fi[(row-1)][(col-1)]) +#define Fdot(row,col) (Fdot[(row-1)][(col-1)]) +#define cauchy(row,col) (cauchy[(row-1)][(col-1)]) +#define setcauchy(row,col) (setcauchy[(row-1)][(col-1)]) +#define setPK(row,col) (setPK[(row-1)][(col-1)]) +#define uv(row,col) (uv[(row-1)][(col-1)]) +#define deltastress(row) (deltastress[(row-1)]) +#define fdotvec(row) (fdotvec[(row-1)]) +#define dsdf(row,col) (dsdf[(row-1)][(col-1)]) +#define dsds(row,col) (dsds[(row-1)][(col-1)]) +#define deltaF(row) (deltaF[(row-1)]) +#define deltaPK(row) (deltaPK[(row-1)]) + + //initialize local variables: + int i,j,m,n; + + double deltastress[6],fdotvec[6]; + double dsdf[6][6]; + double dsds[6][6]; + double jac; + double deltaF[6]; + double deltaPK[6]; + + // zero arrays + for (i = 0; i < 6; ++i) { + deltaF[i] = deltaPK[i] = deltastress[i] = fdotvec[i] = 0.0; + for (j = 0; j < 6; ++j) + dsdf[i][j] = dsds[i][j] = 0.0; + } + + int uv[6][2]; + uv(1,1)=1; uv(1,2)=1; + uv(2,1)=2; uv(2,2)=2; + uv(3,1)=3; uv(3,2)=3; + uv(4,1)=2; uv(4,2)=3; + uv(5,1)=1; uv(5,2)=3; + uv(6,1)=1; uv(6,2)=2; + + for(int ii = 1;ii <= 6;ii++) { + i=uv(ii,1); + j=uv(ii,2); + deltastress(ii)=setcauchy(i,j)-cauchy(i,j); + if(ii>3) deltastress(ii)=deltastress(ii)*2.0; + fdotvec(ii)=Fdot(i,j)*deltat; + } + + for(int ii = 1;ii <= 6;ii++) { + i=uv(ii,1); + j=uv(ii,2); + for(int jj = 1;jj <= 6;jj++) { + m=uv(jj,1); + n=uv(jj,2); + dsds(ii,jj) = Fi(i,m)*Fi(j,n) + Fi(i,n)*Fi(j,m) + Fi(j,m)*Fi(i,n) + Fi(j,n)*Fi(i,m); + for(int l = 1;l <= 3;l++) { + for(int k = 1;k <= 3;k++) { + dsdf(ii,jj) = dsdf(ii,jj) + cauchy(k,l)* + ( Fi(i,k)*Fi(j,l)*Fi(n,m) - Fi(i,m)*Fi(j,l)*Fi(n,k) - Fi(i,k)*Fi(j,m)*Fi(n,l) ); + } + } + } + } + + jac=volume/volume0; + for(int ii = 1;ii <= 6;ii++) { + for(int jj = 1;jj <= 6;jj++) { + dsds(ii,jj)=dsds(ii,jj)*jac/4.0; + dsdf(ii,jj)=dsdf(ii,jj)*jac; + } + } + + for(int ii = 1;ii <= 6;ii++) { + for(int jj = 1;jj <= 6;jj++) { + deltaF(ii)=deltaF(ii)+dsdf(ii,jj)*fdotvec(jj); + } + } + + for(int ii = 1;ii <= 6;ii++) { + for(int jj = 1;jj <= 6;jj++) { + deltaPK(ii)=deltaPK(ii)+alpha*dsds(ii,jj)*deltastress(jj); + } + deltaPK(ii)=deltaPK(ii)+alpha*deltaF(ii); + } + + setPK(1,1)=setPK(1,1)+deltaPK(1); + setPK(2,2)=setPK(2,2)+deltaPK(2); + setPK(3,3)=setPK(3,3)+deltaPK(3); + setPK(2,3)=setPK(2,3)+deltaPK(4); + setPK(3,2)=setPK(3,2)+deltaPK(4); + setPK(1,3)=setPK(1,3)+deltaPK(5); + setPK(3,1)=setPK(3,1)+deltaPK(5); + setPK(1,2)=setPK(1,2)+deltaPK(6); + setPK(2,1)=setPK(2,1)+deltaPK(6); + +#undef F +#undef Fi +#undef Fdot +#undef cauchy +#undef setcauchy +#undef setPK +#undef uv +#undef deltastress +#undef fdotvec +#undef dsdf +#undef dsds +#undef deltaF +#undef deltaPK + +} diff --git a/src/USER-MISC/fix_npt_cauchy.h b/src/USER-MISC/fix_npt_cauchy.h new file mode 100644 index 0000000000..052a5686ad --- /dev/null +++ b/src/USER-MISC/fix_npt_cauchy.h @@ -0,0 +1,316 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(npt/cauchy,FixNPTCauchy) + +#else + +#ifndef LMP_FIX_NPT_CAUCHY_H +#define LMP_FIX_NPT_CAUCHY_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixNPTCauchy : public Fix { + public: + FixNPTCauchy(class LAMMPS *, int, char **); + virtual void post_constructor(); + virtual ~FixNPTCauchy(); + int setmask(); + virtual void init(); + virtual void setup(int); + virtual void initial_integrate(int); + virtual void final_integrate(); + void initial_integrate_respa(int, int, int); + void final_integrate_respa(int, int); + virtual void pre_exchange(); + double compute_scalar(); + virtual double compute_vector(int); + void write_restart(FILE *); + virtual int pack_restart_data(double *); // pack restart data + virtual void restart(char *); + int modify_param(int, char **); + void reset_target(double); + void reset_dt(); + virtual void *extract(const char*,int &); + double memory_usage(); + + protected: + int dimension,which; + double dtv,dtf,dthalf,dt4,dt8,dto; + double boltz,nktv2p,tdof; + double vol0; // reference volume + double t0; // reference temperature + // used for barostat mass + double t_start,t_stop; + double t_current,t_target,ke_target; + double t_freq; + + int tstat_flag; // 1 if control T + int pstat_flag; // 1 if control P + + int pstyle,pcouple,allremap; + int p_flag[6]; // 1 if control P on this dim, 0 if not + double p_start[6],p_stop[6]; + double p_freq[6],p_target[6]; + double omega[6],omega_dot[6]; + double omega_mass[6]; + double p_current[6]; + double drag,tdrag_factor; // drag factor on particle thermostat + double pdrag_factor; // drag factor on barostat + int kspace_flag; // 1 if KSpace invoked, 0 if not + int nrigid; // number of rigid fixes + int dilate_group_bit; // mask for dilation group + int *rfix; // indices of rigid fixes + char *id_dilate; // group name to dilate + class Irregular *irregular; // for migrating atoms after box flips + + int nlevels_respa; + double *step_respa; + + char *id_temp,*id_press; + class Compute *temperature,*pressure; + int tcomputeflag,pcomputeflag; // 1 = compute was created by fix + // 0 = created externally + + double *eta,*eta_dot; // chain thermostat for particles + double *eta_dotdot; + double *eta_mass; + int mtchain; // length of chain + int mtchain_default_flag; // 1 = mtchain is default + + double *etap; // chain thermostat for barostat + double *etap_dot; + double *etap_dotdot; + double *etap_mass; + int mpchain; // length of chain + + int mtk_flag; // 0 if using Hoover barostat + int pdim; // number of barostatted dims + double p_freq_max; // maximum barostat frequency + + double p_hydro; // hydrostatic target pressure + + int nc_tchain,nc_pchain; + double factor_eta; + double sigma[6]; // scaled target stress + double fdev[6]; // deviatoric force on barostat + int deviatoric_flag; // 0 if target stress tensor is hydrostatic + double h0_inv[6]; // h_inv of reference (zero strain) box + int nreset_h0; // interval for resetting h0 + + double mtk_term1,mtk_term2; // Martyna-Tobias-Klein corrections + + int eta_mass_flag; // 1 if eta_mass updated, 0 if not. + int omega_mass_flag; // 1 if omega_mass updated, 0 if not. + int etap_mass_flag; // 1 if etap_mass updated, 0 if not. + int dipole_flag; // 1 if dipole is updated, 0 if not. + int dlm_flag; // 1 if using the DLM rotational integrator, 0 if not + + int scaleyz; // 1 if yz scaled with lz + int scalexz; // 1 if xz scaled with lz + int scalexy; // 1 if xy scaled with ly + int flipflag; // 1 if box flips are invoked as needed + + int pre_exchange_flag; // set if pre_exchange needed for box flips + + double fixedpoint[3]; // location of dilation fixed-point + + void couple(); + virtual void remap(); + void nhc_temp_integrate(); + void nhc_press_integrate(); + + virtual void nve_x(); // may be overwritten by child classes + virtual void nve_v(); + virtual void nh_v_press(); + virtual void nh_v_temp(); + virtual void compute_temp_target(); + virtual int size_restart_global(); + + void compute_sigma(); + void compute_deviatoric(); + double compute_strain_energy(); + void compute_press_target(); + void nh_omega_dot(); + + // Implementation of CauchyStat + char *id_store; // fix id of the STORE fix for retaining data + class FixStore *init_store; // fix pointer to STORE fix + double H0[3][3]; // shape matrix for the undeformed cell + double h_old[6]; // previous time step shape matrix for + // the undeformed cell + double invH0[3][3]; // inverse of H0; + double CSvol0; // volume of undeformed cell + double setPK[3][3]; // current set values of the PK stress + // (this is modified until the cauchy + // stress converges) + double alpha; // integration parameter for the cauchystat + int initPK; // 1 if setPK needs to be initialized either + // from cauchy or restart, else 0 + int restartPK; // Read PK stress from the previous run + int restart_stored; // values of PK stress from the previous step stored + int initRUN; // 0 if run not initialized + // (pressure->vector not computed yet), + // else 1 (pressure->vector available) + + void CauchyStat_init(); + void CauchyStat_cleanup(); + void CauchyStat(); + void CauchyStat_Step(double (&Fi)[3][3], double (&Fdot)[3][3], + double (&cauchy)[3][3], double (&setcauchy)[3][3], + double (&setPK)[3][3], double volume, double volume0, + double deltat, double alpha); + +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Target temperature for fix npt/cauchy cannot be 0.0 + +Self-explanatory. + +E: Invalid fix npt/cauchy command for a 2d simulation + +Cannot control z dimension in a 2d model. + +E: Fix npt/cauchy dilate group ID does not exist + +Self-explanatory. + +E: Invalid fix npt/cauchy command pressure settings + +If multiple dimensions are coupled, those dimensions must be +specified. + +E: Cannot use fix npt/cauchy on a non-periodic dimension + +When specifying a diagonal pressure component, the dimension must be +periodic. + +E: Cannot use fix npt/cauchy on a 2nd non-periodic dimension + +When specifying an off-diagonal pressure component, the 2nd of the two +dimensions must be periodic. E.g. if the xy component is specified, +then the y dimension must be periodic. + +E: Cannot use fix npt/cauchy with yz scaling when z is non-periodic dimension + +The 2nd dimension in the barostatted tilt factor must be periodic. + +E: Cannot use fix npt/cauchy with xz scaling when z is non-periodic dimension + +The 2nd dimension in the barostatted tilt factor must be periodic. + +E: Cannot use fix npt/cauchy with xy scaling when y is non-periodic dimension + +The 2nd dimension in the barostatted tilt factor must be periodic. + +E: Cannot use fix npt/cauchy with both yz dynamics and yz scaling + +Self-explanatory. + +E: Cannot use fix npt/cauchy with both xz dynamics and xz scaling + +Self-explanatory. + +E: Cannot use fix npt/cauchy with both xy dynamics and xy scaling + +Self-explanatory. + +E: Can not specify Pxy/Pxz/Pyz in fix npt/cauchy with non-triclinic box + +Only triclinic boxes can be used with off-diagonal pressure components. +See the region prism command for details. + +E: Invalid fix npt/cauchy pressure settings + +Settings for coupled dimensions must be the same. + +E: Using update dipole flag requires atom style sphere + +Self-explanatory. + +E: Using update dipole flag requires atom attribute mu + +Self-explanatory. + +E: Fix npt/cauchy damping parameters must be > 0.0 + +Self-explanatory. + +E: Cannot use fix npt/cauchy and fix deform on same component of stress tensor + +This would be changing the same box dimension twice. + +E: Temperature ID for fix npt/cauchy does not exist + +Self-explanatory. + +E: Pressure ID for fix npt/cauchy does not exist + +Self-explanatory. + +E: Non-numeric pressure - simulation unstable + +UNDOCUMENTED + +E: Fix npt/cauchy has tilted box too far in one step - periodic cell is too far from equilibrium state + +Self-explanatory. The change in the box tilt is too extreme +on a short timescale. + +E: Could not find fix_modify temperature ID + +The compute ID for computing temperature does not exist. + +E: Fix_modify temperature ID does not compute temperature + +The compute ID assigned to the fix must compute temperature. + +W: Temperature for fix modify is not for group all + +The temperature compute is being used with a pressure calculation +which does operate on group all, so this may be inconsistent. + +E: Pressure ID for fix modify does not exist + +Self-explanatory. + +E: Could not find fix_modify pressure ID + +The compute ID for computing pressure does not exist. + +E: Fix_modify pressure ID does not compute pressure + +The compute ID assigned to the fix must compute pressure. + +U: The dlm flag must be used with update dipole + +Self-explanatory. + +*/