forked from lijiext/lammps
Merge pull request #1096 from ProfessorMiller/master
Changes to the NH fix enabling Cauchy stress control (Cauhchystat) du…
This commit is contained in:
Axel Kohlmeyer
GitHub
commit
f1c79fb914
|
|
@ -115,6 +115,7 @@ OPT.
|
|||
* :doc:`npt (iko) <fix_nh>`
|
||||
* :doc:`npt/asphere (o) <fix_npt_asphere>`
|
||||
* :doc:`npt/body <fix_npt_body>`
|
||||
* :doc:`npt/cauchy <fix_npt_cauchy>`
|
||||
* :doc:`npt/eff <fix_nh_eff>`
|
||||
* :doc:`npt/sphere (o) <fix_npt_sphere>`
|
||||
* :doc:`npt/uef <fix_nh_uef>`
|
||||
|
|
|
|||
|
|
@ -267,6 +267,7 @@ accelerated styles exist.
|
|||
* :doc:`npt <fix_nh>` - constant NPT time integration via Nose/Hoover
|
||||
* :doc:`npt/asphere <fix_npt_asphere>` - NPT for aspherical particles
|
||||
* :doc:`npt/body <fix_npt_body>` - NPT for body particles
|
||||
* :doc:`npt/cauchy <fix_npt_cauchy>` - NPT with Cauchy stress
|
||||
* :doc:`npt/eff <fix_nh_eff>` - NPT for nuclei and electrons in the electron force field model
|
||||
* :doc:`npt/sphere <fix_npt_sphere>` - NPT for spherical particles
|
||||
* :doc:`npt/uef <fix_nh_uef>` - NPT style time integration with diagonal flow
|
||||
|
|
|
|||
|
|
@ -132,7 +132,7 @@ integrators derived by Tuckerman et al in :ref:`(Tuckerman) <nh-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.
|
||||
|
||||
|
|
|
|||
|
|
@ -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 <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) <nc-Shinoda>`, which combine the hydrostatic equations of
|
||||
Martyna, Tobias and Klein in :ref:`(Martyna) <nc-Martyna>` with the strain
|
||||
energy proposed by Parrinello and Rahman in
|
||||
:ref:`(Parrinello) <nc-Parrinello>`. The time integration schemes closely
|
||||
follow the time-reversible measure-preserving Verlet and rRESPA
|
||||
integrators derived by Tuckerman et al in :ref:`(Tuckerman) <nc-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 <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 <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 <create_box>`, :doc:`read data <read_data>`, and :doc:`read\_restart <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 <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 <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 <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 <fix_nve>` or
|
||||
:doc:`fix nvt <fix_nh>` 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) <nc-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) <nc-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 <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 <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 <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 <fix_nh>` or :doc:`fix temp/rescale <fix_temp_rescale>`
|
||||
commands.
|
||||
|
||||
See the :doc:`Howto thermostat <Howto_thermostat>` and :doc:`Howto barostat <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 <compute_temp>` and :doc:`compute pressure <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 <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 <compute_modify>` command. Or you can print this
|
||||
temperature or pressure during thermodynamic output via the
|
||||
:doc:`thermo\_style custom <thermo_style>` 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 <compute>` 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 <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 <compute>` 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 <run_style>`. 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 <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 <restart>`. See the
|
||||
:doc:`read\_restart <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 <fix_modify>` *temp* and *press* options are
|
||||
supported by this fix. You can use them to assign a
|
||||
:doc:`compute <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 <compute_pressure>` 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 <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 <thermo_style>`.
|
||||
|
||||
This fix computes a global scalar and a global vector of quantities,
|
||||
which can be accessed by various :doc:`output commands <Howto_output>`.
|
||||
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 <run>` command. See the :doc:`run <run>` command for details of
|
||||
how to do this.
|
||||
|
||||
This fix is not invoked during :doc:`energy minimization <minimize>`.
|
||||
|
||||
|
||||
----------
|
||||
|
||||
|
||||
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 <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 <create_box>`,
|
||||
:doc:`read data <read_data>`, and :doc:`read\_restart <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 <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 <fix_nve>`, :doc:`fix\_modify <fix_modify>`,
|
||||
:doc:`run\_style <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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -0,0 +1 @@
|
|||
../../../../potentials/NiAlH_jea.eam.alloy
|
||||
|
|
@ -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.
|
||||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
|
|
@ -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.
|
||||
|
||||
*/
|
||||
Loading…
Reference in New Issue