ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'develop' into collected-small-fixes

Browse Source
This commit is contained in:
Axel Kohlmeyer
2024-02-22 04:56:43 -05:00
parent a466d74f31 e51a65696d
commit 00615d94ad
13 changed files with 1880 additions and 428 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
doc/src/Commands_fix.rst
View File

@ -61,6 +61,7 @@ OPT.
* :doc:`controller <fix_controller>`
* :doc:`damping/cundall <fix_damping_cundall>`
* :doc:`deform (k) <fix_deform>`
* :doc:`deform/pressure <fix_deform_pressure>`
* :doc:`deposit <fix_deposit>`
* :doc:`dpd/energy (k) <fix_dpd_energy>`
* :doc:`drag <fix_drag>`

1
doc/src/fix.rst
View File

@ -226,6 +226,7 @@ accelerated styles exist.
* :doc:`controller <fix_controller>` - apply control loop feedback mechanism
* :doc:`damping/cundall <fix_damping_cundall>` - Cundall non-viscous damping for granular simulations
* :doc:`deform <fix_deform>` - change the simulation box size/shape
* :doc:`deform/pressure <fix_deform_pressure>` - change the simulation box size/shape with additional loading conditions
* :doc:`deposit <fix_deposit>` - add new atoms above a surface
* :doc:`dpd/energy <fix_dpd_energy>` - constant energy dissipative particle dynamics
* :doc:`drag <fix_drag>` - drag atoms towards a defined coordinate

181
doc/src/fix_deform.rst
View File

@ -4,6 +4,9 @@
fix deform command
==================
:doc:`fix deform/pressure <fix_deform_pressure>` command
========================================================
Accelerator Variants: *deform/kk*
Syntax
@ -11,18 +14,18 @@ Syntax
.. code-block:: LAMMPS
fix ID group-ID deform N parameter args ... keyword value ...
fix ID group-ID fix_style N parameter style args ... keyword value ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* deform = style name of this fix command
* fix_style = *deform* or *deform/pressure*
* N = perform box deformation every this many timesteps
* one or more parameter/arg pairs may be appended
* one or more parameter/style/args sequences of arguments may be appended
.. parsed-literal::
parameter = *x* or *y* or *z* or *xy* or *xz* or *yz*
*x*, *y*, *z* args = style value(s)
style = *final* or *delta* or *scale* or *vel* or *erate* or *trate* or *volume* or *wiggle* or *variable*
style = *final* or *delta* or *scale* or *vel* or *erate* or *trate* or *volume* or *wiggle* or *variable* or *pressure* or *pressure/mean*
*final* values = lo hi
lo hi = box boundaries at end of run (distance units)
*delta* values = dlo dhi
@ -43,8 +46,15 @@ Syntax
*variable* values = v_name1 v_name2
v_name1 = variable with name1 for box length change as function of time
v_name2 = variable with name2 for change rate as function of time
*pressure* values = target gain (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
target = target pressure (pressure units)
gain = proportional gain constant (1/(time * pressure) or 1/time units)
*pressure/mean* values = target gain (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
target = target pressure (pressure units)
gain = proportional gain constant (1/(time * pressure) or 1/time units)
*xy*, *xz*, *yz* args = style value
style = *final* or *delta* or *vel* or *erate* or *trate* or *wiggle*
style = *final* or *delta* or *vel* or *erate* or *trate* or *wiggle* or *variable*
*final* value = tilt
tilt = tilt factor at end of run (distance units)
*delta* value = dtilt
@ -62,9 +72,12 @@ Syntax
*variable* values = v_name1 v_name2
v_name1 = variable with name1 for tilt change as function of time
v_name2 = variable with name2 for change rate as function of time
*pressure* values = target gain (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
target = target pressure (pressure units)
gain = proportional gain constant (1/(time * pressure) or 1/time units)
* zero or more keyword/value pairs may be appended
* keyword = *remap* or *flip* or *units*
* keyword = *remap* or *flip* or *units* or *couple* or *vol/balance/p* or *max/rate* or *normalize/pressure*
.. parsed-literal::
@ -77,6 +90,15 @@ Syntax
*units* value = *lattice* or *box*
lattice = distances are defined in lattice units
box = distances are defined in simulation box units
*couple* value = *none* or *xyz* or *xy* or *yz* or *xz* (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
couple pressure values of various dimensions
*vol/balance/p* value = *yes* or *no* (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
Modifies the behavior of the *volume* option to try and balance pressures
*max/rate* value = *rate* (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
rate = maximum strain rate for pressure control
*normalize/pressure* value = *yes* or *no* (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
Modifies pressure controls such that the deviation in pressure is normalized by the target pressure
Examples
""""""""
@ -88,6 +110,8 @@ Examples
fix 1 all deform 1 xy erate 0.001 remap v
fix 1 all deform 10 y delta -0.5 0.5 xz vel 1.0
See examples for :doc:`fix deform/pressure <fix_deform_pressure>` on its doc page
Description
"""""""""""
@ -95,29 +119,46 @@ Change the volume and/or shape of the simulation box during a dynamics
run. Orthogonal simulation boxes have 3 adjustable parameters
(x,y,z). Triclinic (non-orthogonal) simulation boxes have 6
adjustable parameters (x,y,z,xy,xz,yz). Any or all of them can be
adjusted independently and simultaneously by this command.
adjusted independently and simultaneously.
This fix can be used to perform non-equilibrium MD (NEMD) simulations
of a continuously strained system. See the :doc:`fix nvt/sllod <fix_nvt_sllod>` and :doc:`compute temp/deform <compute_temp_deform>` commands for more details. Note
that simulation of a continuously extended system (extensional flow)
can be modeled using the :ref:`UEF package <PKG-UEF>` and its :doc:`fix commands <fix_nh_uef>`.
The fix deform command allows use of all the arguments listed above,
except those flagged as available ONLY for the :doc:`fix
deform/pressure <fix_deform_pressure>` command, which are
pressure-based controls. The fix deform/pressure command allows use
of all the arguments listed above.
The rest of this doc page explains the options common to both
commands. The :doc:`fix deform/pressure <fix_deform_pressure>` doc
page explains the options available ONLY with the fix deform/pressure
command. Note that a simulation can define only a single deformation
command: fix deform or fix deform/pressure.
Both these fixes can be used to perform non-equilibrium MD (NEMD)
simulations of a continuously strained system. See the :doc:`fix
nvt/sllod <fix_nvt_sllod>` and :doc:`compute temp/deform
<compute_temp_deform>` commands for more details. Note that
simulation of a continuously extended system (extensional flow) can be
modeled using the :ref:`UEF package <PKG-UEF>` and its :doc:`fix
commands <fix_nh_uef>`.
For the *x*, *y*, *z* parameters, the associated dimension cannot be
shrink-wrapped. For the *xy*, *yz*, *xz* parameters, the associated
second dimension cannot be shrink-wrapped. Dimensions not varied by this
command can be periodic or non-periodic. Dimensions corresponding to
unspecified parameters can also be controlled by a :doc:`fix npt <fix_nh>` or :doc:`fix nph <fix_nh>` command.
second dimension cannot be shrink-wrapped. Dimensions not varied by
this command can be periodic or non-periodic. Dimensions
corresponding to unspecified parameters can also be controlled by a
:doc:`fix npt <fix_nh>` or :doc:`fix nph <fix_nh>` command.
The size and shape of the simulation box at the beginning of the
simulation run were either specified by the
:doc:`create_box <create_box>` or :doc:`read_data <read_data>` or
:doc:`read_restart <read_restart>` command used to setup the simulation
initially if it is the first run, or they are the values from the end
of the previous run. 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. If
fix deform changes the xy,xz,yz tilt factors, then the simulation box
must be triclinic, even if its initial tilt factors are 0.0.
simulation run were either specified by the :doc:`create_box
<create_box>` or :doc:`read_data <read_data>` or :doc:`read_restart
<read_restart>` command used to setup the simulation initially if it
is the first run, or they are the values from the end of the previous
run. 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. If fix deform
changes the xy,xz,yz tilt factors, then the simulation box must be
triclinic, even if its initial tilt factors are 0.0.
As described below, the desired simulation box size and shape at the
end of the run are determined by the parameters of the fix deform
@ -258,21 +299,22 @@ of the units keyword below.
The *variable* style changes the specified box length dimension by
evaluating a variable, which presumably is a function of time. The
variable with *name1* must be an :doc:`equal-style variable <variable>`
and should calculate a change in box length in units of distance.
Note that this distance is in box units, not lattice units; see the
discussion of the *units* keyword below. The formula associated with
variable *name1* can reference the current timestep. Note that it
should return the "change" in box length, not the absolute box length.
This means it should evaluate to 0.0 when invoked on the initial
timestep of the run following the definition of fix deform. It should
evaluate to a value > 0.0 to dilate the box at future times, or a
value < 0.0 to compress the box.
variable with *name1* must be an :doc:`equal-style variable
<variable>` and should calculate a change in box length in units of
distance. Note that this distance is in box units, not lattice units;
see the discussion of the *units* keyword below. The formula
associated with variable *name1* can reference the current timestep.
Note that it should return the "change" in box length, not the
absolute box length. This means it should evaluate to 0.0 when
invoked on the initial timestep of the run following the definition of
fix deform. It should evaluate to a value > 0.0 to dilate the box at
future times, or a value < 0.0 to compress the box.
The variable *name2* must also be an :doc:`equal-style variable <variable>` and should calculate the rate of box length
change, in units of distance/time, i.e. the time-derivative of the
*name1* variable. This quantity is used internally by LAMMPS to reset
atom velocities when they cross periodic boundaries. It is computed
The variable *name2* must also be an :doc:`equal-style variable
<variable>` and should calculate the rate of box length change, in
units of distance/time, i.e. the time-derivative of the *name1*
variable. This quantity is used internally by LAMMPS to reset atom
velocities when they cross periodic boundaries. It is computed
internally for the other styles, but you must provide it when using an
arbitrary variable.
@ -414,12 +456,13 @@ can reference the current timestep. Note that it should return the
should evaluate to 0.0 when invoked on the initial timestep of the run
following the definition of fix deform.
The variable *name2* must also be an :doc:`equal-style variable <variable>` and should calculate the rate of tilt change,
in units of distance/time, i.e. the time-derivative of the *name1*
variable. This quantity is used internally by LAMMPS to reset atom
velocities when they cross periodic boundaries. It is computed
internally for the other styles, but you must provide it when using an
arbitrary variable.
The variable *name2* must also be an :doc:`equal-style variable
<variable>` and should calculate the rate of tilt change, in units of
distance/time, i.e. the time-derivative of the *name1* variable. This
quantity is used internally by LAMMPS to reset atom velocities when
they cross periodic boundaries. It is computed internally for the
other styles, but you must provide it when using an arbitrary
variable.
Here is an example of using the *variable* style to perform the same
box deformation as the *wiggle* style formula listed above, where we
@ -510,33 +553,40 @@ box without explicit remapping of their coordinates.
.. note::
For non-equilibrium MD (NEMD) simulations using "remap v" it is
usually desirable that the fluid (or flowing material, e.g. granular
particles) stream with a velocity profile consistent with the
deforming box. As mentioned above, using a thermostat such as :doc:`fix nvt/sllod <fix_nvt_sllod>` or :doc:`fix lavgevin <fix_langevin>`
(with a bias provided by :doc:`compute temp/deform <compute_temp_deform>`), will typically accomplish
that. If you do not use a thermostat, then there is no driving force
pushing the atoms to flow in a manner consistent with the deforming
box. E.g. for a shearing system the box deformation velocity may vary
usually desirable that the fluid (or flowing material,
e.g. granular particles) stream with a velocity profile consistent
with the deforming box. As mentioned above, using a thermostat
such as :doc:`fix nvt/sllod <fix_nvt_sllod>` or :doc:`fix lavgevin
<fix_langevin>` (with a bias provided by :doc:`compute temp/deform
<compute_temp_deform>`), will typically accomplish that. If you do
not use a thermostat, then there is no driving force pushing the
atoms to flow in a manner consistent with the deforming box.
E.g. for a shearing system the box deformation velocity may vary
from 0 at the bottom to 10 at the top of the box. But the stream
velocity profile of the atoms may vary from -5 at the bottom to +5 at
the top. You can monitor these effects using the :doc:`fix ave/chunk <fix_ave_chunk>`, :doc:`compute temp/deform <compute_temp_deform>`, and :doc:`compute temp/profile <compute_temp_profile>` commands. One way to induce
atoms to stream consistent with the box deformation is to give them an
velocity profile of the atoms may vary from -5 at the bottom to +5
at the top. You can monitor these effects using the :doc:`fix
ave/chunk <fix_ave_chunk>`, :doc:`compute temp/deform
<compute_temp_deform>`, and :doc:`compute temp/profile
<compute_temp_profile>` commands. One way to induce atoms to
stream consistent with the box deformation is to give them an
initial velocity profile, via the :doc:`velocity ramp <velocity>`
command, that matches the box deformation rate. This also typically
helps the system come to equilibrium more quickly, even if a
thermostat is used.
command, that matches the box deformation rate. This also
typically helps the system come to equilibrium more quickly, even
if a thermostat is used.
.. note::
If a :doc:`fix rigid <fix_rigid>` is defined for rigid bodies, and
*remap* is set to *x*, then the center-of-mass coordinates of rigid
bodies will be remapped to the changing simulation box. This will be
done regardless of whether atoms in the rigid bodies are in the fix
deform group or not. The velocity of the centers of mass are not
remapped even if *remap* is set to *v*, since :doc:`fix nvt/sllod <fix_nvt_sllod>` does not currently do anything special
bodies will be remapped to the changing simulation box. This will
be done regardless of whether atoms in the rigid bodies are in the
fix deform group or not. The velocity of the centers of mass are
not remapped even if *remap* is set to *v*, since :doc:`fix
nvt/sllod <fix_nvt_sllod>` does not currently do anything special
for rigid particles. If you wish to perform a NEMD simulation of
rigid particles, you can either thermostat them independently or
include a background fluid and thermostat the fluid via :doc:`fix nvt/sllod <fix_nvt_sllod>`.
include a background fluid and thermostat the fluid via :doc:`fix
nvt/sllod <fix_nvt_sllod>`.
The *flip* keyword allows the tilt factors for a triclinic box to
exceed half the distance of the parallel box length, as discussed
@ -568,7 +618,8 @@ command if you want to include lattice spacings in a variable formula.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This fix will restore the initial box settings from :doc:`binary restart files <restart>`, which allows the fix to be properly continue
This fix will restore the initial box settings from :doc:`binary
restart files <restart>`, which allows the fix to be properly continue
deformation, when using the start/stop options of the :doc:`run <run>`
command. None of the :doc:`fix_modify <fix_modify>` options are
relevant to this fix. No global or per-atom quantities are stored by
@ -586,12 +637,14 @@ Restrictions
You cannot apply x, y, or z deformations to a dimension that is
shrink-wrapped via the :doc:`boundary <boundary>` command.
You cannot apply xy, yz, or xz deformations to a second dimension (y in
xy) that is shrink-wrapped via the :doc:`boundary <boundary>` command.
You cannot apply xy, yz, or xz deformations to a second dimension (y
in xy) that is shrink-wrapped via the :doc:`boundary <boundary>`
command.
Related commands
""""""""""""""""
:doc:`fix deform/pressure <fix_deform_pressure>`,
:doc:`change_box <change_box>`
Default

319
doc/src/fix_deform_pressure.rst Normal file
View File

@ -0,0 +1,319 @@
.. index:: fix deform/pressure
fix deform/pressure command
===========================
Syntax
""""""
.. parsed-literal::
fix ID group-ID deform/pressure N parameter style args ... keyword value ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* deform/pressure = style name of this fix command
* N = perform box deformation every this many timesteps
* one or more parameter/arg sequences may be appended
.. parsed-literal::
parameter = *x* or *y* or *z* or *xy* or *xz* or *yz* or *box*
*x*, *y*, *z* args = style value(s)
style = *final* or *delta* or *scale* or *vel* or *erate* or *trate* or *volume* or *wiggle* or *variable* or *pressure* or *pressure/mean*
*pressure* values = target gain
target = target pressure (pressure units)
gain = proportional gain constant (1/(time * pressure) or 1/time units)
*pressure/mean* values = target gain
target = target pressure (pressure units)
gain = proportional gain constant (1/(time * pressure) or 1/time units)
NOTE: All other styles are documented by the :doc:`fix deform <fix_deform>` command
*xy*, *xz*, *yz* args = style value
style = *final* or *delta* or *vel* or *erate* or *trate* or *wiggle* or *variable* or *pressure*
*pressure* values = target gain
target = target pressure (pressure units)
gain = proportional gain constant (1/(time * pressure) or 1/time units)
NOTE: All other styles are documented by the :doc:`fix deform <fix_deform>` command
*box* = style value
style = *volume* or *pressure*
*volume* value = none = isotropically adjust system to preserve volume of system
*pressure* values = target gain
target = target mean pressure (pressure units)
gain = proportional gain constant (1/(time * pressure) or 1/time units)
* zero or more keyword/value pairs may be appended
* keyword = *remap* or *flip* or *units* or *couple* or *vol/balance/p* or *max/rate* or *normalize/pressure*
.. parsed-literal::
*couple* value = *none* or *xyz* or *xy* or *yz* or *xz*
couple pressure values of various dimensions
*vol/balance/p* value = *yes* or *no*
Modifies the behavior of the *volume* option to try and balance pressures
*max/rate* value = *rate*
rate = maximum strain rate for pressure control
*normalize/pressure* value = *yes* or *no*
Modifies pressure controls such that the deviation in pressure is normalized by the target pressure
NOTE: All other keywords are documented by the :doc:`fix deform <fix_deform>` command
Examples
""""""""
.. code-block:: LAMMPS
fix 1 all deform/pressure 1 x pressure 2.0 0.1 normalize/pressure yes max/rate 0.001
fix 1 all deform/pressure 1 x trate 0.1 y volume z volume vol/balance/p yes
fix 1 all deform/pressure 1 x trate 0.1 y pressure/mean 0.0 1.0 z pressure/mean 0.0 1.0
Description
"""""""""""
.. versionadded:: TBD
This fix is an extension of the :doc:`fix deform <fix_deform>`
command, which allows all of its options to be used as well as new
pressure-based controls implemented by this command.
All arguments described on the :doc:`fix deform <fix_deform>` doc page
also apply to this fix unless otherwise noted below. The rest of this
doc page explains the arguments specific to this fix. Note that a
simulation can define only a single deformation command: fix deform or
fix deform/pressure.
----------
For the *x*, *y*, and *z* parameters, this is the meaning of the
styles and values provided by this fix.
The *pressure* style adjusts a dimension's box length to control the
corresponding component of the pressure tensor. This option attempts to
maintain a specified target pressure using a linear controller where the
box length :math:`L` evolves according to the equation
.. parsed-literal::
\frac{d L(t)}{dt} = L(t) k (P_t - P)
where :math:`k` is a proportional gain constant, :math:`P_t` is the target
pressure, and :math:`P` is the current pressure along that dimension. This
approach is similar to the method used to control the pressure by
:doc:`fix press/berendsen <fix_press_berendsen>`. The target pressure
accepts either a constant numeric value or a LAMMPS :ref:`variable <variable>`.
Notably, this variable can be a function of time or other components of
the pressure tensor. By default, :math:`k` has units of 1/(time * pressure)
although this will change if the *normalize/pressure* option is set as
:ref:`discussed below <deform_normalize>`. There is no proven method
to choosing an appropriate value of :math:`k` as it will depend on the
specific details of a simulation. Testing different values is recommended.
By default, there is no limit on the resulting strain rate in any dimension.
A maximum limit can be applied using the :ref:`max/rate <deform_max_rate>`
option. Akin to :doc:`fix nh <fix_nh>`, pressures in different dimensions
can be coupled using the :ref:`couple <deform_couple>` option. This means
the instantaneous pressure along coupled dimensions are averaged and the box
strains identically along the coupled dimensions.
The *pressure/mean* style changes a dimension's box length to maintain
a constant mean pressure defined as the trace of the pressure tensor.
This option has identical arguments to the *pressure* style and a similar
functional equation, except the current and target pressures refer to the
mean trace of the pressure tensor. All options for the *pressure* style
also apply to the *pressure/mean* style except for the
:ref:`couple <deform_couple>` option.
Note that while this style can be identical to coupled *pressure* styles,
it is generally not the same. For instance in 2D, a coupled *pressure*
style in the *x* and *y* dimensions would be equivalent to using the
*pressure/mean* style with identical settings in each dimension. However,
it would not be the same if settings (e.g. gain constants) were used in
the *x* and *y* dimensions or if the *pressure/mean* command was only applied
along one dimension.
----------
For the *xy*, *xz*, and *yz* parameters, this is the meaning of the
styles and values provided by this fix. Note that changing the
tilt factors of a triclinic box does not change its volume.
The *pressure* style adjusts a tilt factor to control the corresponding
off-diagonal component of the pressure tensor. This option attempts to
maintain a specified target value using a linear controller where the
tilt factor T evolves according to the equation
.. parsed-literal::
\frac{d T(t)}{dt} = L(t) k (P - P_t)
where :math:`k` is a proportional gain constant, :math:`P_t` is the
target pressure, :math:`P` is the current pressure, and :math:`L` is
the perpendicular box length. The target pressure accepts either a
constant numeric value or a LAMMPS :ref:`variable
<variable>`. Notably, this variable can be a function of time or other
components of the pressure tensor. By default, :math:`k` has units of
1/(time * pressure) although this will change if the
*normalize/pessure* option is set as :ref:`discussed below
<deform_normalize>`. There is no proven method to choosing an
appropriate value of :math:`k` as it will depend on the specific
details of a simulation and testing different values is
recommended. One can also apply a maximum limit to the magnitude of
the applied strain using the :ref:`max/rate <deform_max_rate>` option.
----------
The *box* parameter provides an additional control over the *x*, *y*,
and *z* box lengths by isotropically dilating or contracting the box
to either maintain a fixed mean pressure or volume. This isotropic
scaling is applied after the box is deformed by the above *x*, *y*,
*z*, *xy*, *xz*, and *yz* styles, acting as a second deformation
step. This parameter will change the overall strain rate in the *x*,
*y*, or *z* dimensions. This parameter can only be used in
combination with the *x*, *y*, or *z* commands: *vel*, *erate*,
*trate*, *pressure*, or *wiggle*. This is the meaning of its styles
and values.
The *volume* style isotropically scales box lengths to maintain a constant
box volume in response to deformation from other parameters. This style
may be useful in scenarios where one wants to apply a constant deviatoric
pressure using *pressure* styles in the *x*, *y*, and *z* dimensions (
deforming the shape of the box), while maintaining a constant volume.
The *pressure* style isotropically scales box lengths in an attempt to
maintain a target mean pressure (the trace of the pressure tensor) of the
system. This is accomplished by isotropically scaling all box lengths
:math:`L` by an additional factor of :math:`k (P_t - P_m)` where :math:`k`
is the proportional gain constant, :math:`P_t` is the target pressure, and
:math:`P_m` is the current mean pressure. This style may be useful in
scenarios where one wants to apply a constant deviatoric strain rate
using various strain-based styles (e.g. *trate*) along the *x*, *y*, and *z*
dimensions (deforming the shape of the box), while maintaining a mean pressure.
----------
The optional keywords provided by this fix are described below.
.. _deform_normalize:
The *normalize/pressure* keyword changes how box dimensions evolve when
using the *pressure* or *pressure/mean* deformation styles. If the
*deform/normalize* value is set to *yes*, then the deviation from the
target pressure is normalized by the absolute value of the target
pressure such that the proportional gain constant scales a percentage
error and has units of 1/time. If the target pressure is ever zero, this
will produce an error unless the *max/rate* keyword is defined,
described below, which will cap the divergence.
.. _deform_max_rate:
The *max/rate* keyword sets an upper threshold, *rate*, that limits the
maximum magnitude of the instantaneous strain rate applied in any dimension.
This keyword only applies to the *pressure* and *pressure/mean* options. If
a pressure-controlled rate is used for both *box* and either *x*, *y*, or
*z*, then this threshold will apply separately to each individual controller
such that the cumulative strain rate on a box dimension may be up to twice
the value of *rate*.
.. _deform_couple:
The *couple* keyword allows two or three of the diagonal components of
the pressure tensor to be "coupled" together for the *pressure* option.
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. If a *pressure*
style is defined for more than one coupled dimension, the target pressures
and gain constants must be identical. Alternatively, if a *pressure*
style is only defined for one of the coupled dimensions, its settings are
copied to other dimensions with undefined styles. *Couple xyz* can be used
for a 2d simulation; the *z* dimension is simply ignored.
.. _deform_balance:
The *vol/balance/p* keyword modifies the behavior of the *volume* style when
applied to two of the *x*, *y*, and *z* dimensions. Instead of straining
the two dimensions in lockstep, the two dimensions are allowed to
separately dilate or contract in a manner to maintain a constant
volume while simultaneously trying to keep the pressure along each
dimension equal using a method described in :ref:`(Huang2014) <Huang2014>`.
----------
If any pressure controls are used, this fix computes a temperature and
pressure each timestep. To do this, the fix creates its own computes
of style "temp" and "pressure", as if these commands had been issued:
.. code-block:: LAMMPS
compute fix-ID_temp group-ID temp
compute fix-ID_press group-ID pressure fix-ID_temp
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", and the group for the new computes is the same
as the fix group.
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 this fix's temperature or pressure via the
:doc:`compute_modify <compute_modify>` command or 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.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This fix will restore the initial box settings from :doc:`binary
restart files <restart>`, which allows the fix to be properly continue
deformation, when using the start/stop options of the :doc:`run <run>`
command. No global or per-atom quantities are stored by this fix for
access by various :doc:`output commands <Howto_output>`.
If any pressure controls are used, the :doc:`fix_modify <fix_modify>`
*temp* and *press* options are supported by this fix, unlike in
:doc:`fix deform <fix_deform>`. You can use them to assign a
:doc:`compute <compute>` you have defined to this fix which will be
used in its temperature and pressure calculations. 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.
This fix can perform deformation 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
""""""""""""
You cannot apply x, y, or z deformations to a dimension that is
shrink-wrapped via the :doc:`boundary <boundary>` command.
You cannot apply xy, yz, or xz deformations to a second dimension (y
in xy) that is shrink-wrapped via the :doc:`boundary <boundary>`
command.
Related commands
""""""""""""""""
:doc:`fix deform <fix_deform>`, :doc:`change_box <change_box>`
Default
"""""""
The option defaults are normalize/pressure = no.
----------
.. _Huang2014:
**(Huang2014)** X. Huang, "Exploring critical-state behavior using DEM",
Doctoral dissertation, Imperial College. (2014). https://doi.org/10.25560/25316

2
src/.gitignore vendored
View File

@ -791,6 +791,8 @@
/fix_cmap.h
/fix_damping_cundall.cpp
/fix_damping_cundall.h
/fix_deform_pressure.cpp
/fix_deform_pressure.h
/fix_dpd_energy.cpp
/fix_dpd_energy.h
/fix_efield_tip4p.cpp

944
src/EXTRA-FIX/fix_deform_pressure.cpp Normal file
View File

@ -0,0 +1,944 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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 author: Joel Clemmer (SNL)
------------------------------------------------------------------------- */
#include "fix_deform_pressure.h"
#include "atom.h"
#include "comm.h"
#include "compute.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "input.h"
#include "irregular.h"
#include "kspace.h"
#include "lattice.h"
#include "math_const.h"
#include "modify.h"
#include "update.h"
#include "variable.h"
#include <cmath>
#include <cstring>
#include <unordered_map>
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
enum { NOCOUPLE = 0, XYZ, XY, YZ, XZ };
/* ---------------------------------------------------------------------- */
FixDeformPressure::FixDeformPressure(LAMMPS *lmp, int narg, char **arg) :
FixDeform(lmp, narg, arg), id_temp(nullptr), id_press(nullptr)
{
// set defaults
set_extra = new SetExtra[7];
memset(set_extra, 0, 7 * sizeof(SetExtra));
memset(&set_box, 0, sizeof(Set));
// parse only parameter/style arguments specific to this child class
int index, iarg;
int i = 0;
while (i < leftover_iarg.size()) {
iarg = leftover_iarg[i];
if (strcmp(arg[iarg], "x") == 0 ||
strcmp(arg[iarg], "y") == 0 ||
strcmp(arg[iarg], "z") == 0) {
if (strcmp(arg[iarg], "x") == 0) index = 0;
else if (strcmp(arg[iarg], "y") == 0) index = 1;
else if (strcmp(arg[iarg], "z") == 0) index = 2;
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure", error);
if (strcmp(arg[iarg + 1], "pressure") == 0) {
if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure pressure", error);
set[index].style = PRESSURE;
if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) {
set_extra[index].ptarget = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
} else {
set_extra[index].pstr = utils::strdup(&arg[iarg + 2][2]);
set_extra[index].pvar_flag = 1;
}
set_extra[index].pgain = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
i += 4;
} else if (strcmp(arg[iarg + 1], "pressure/mean") == 0) {
if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure pressure/mean", error);
set[index].style = PMEAN;
if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) {
set_extra[index].ptarget = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
} else {
set_extra[index].pstr = utils::strdup(&arg[iarg + 2][2]);
set_extra[index].pvar_flag = 1;
}
set_extra[index].pgain = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
i += 4;
} else error->all(FLERR, "Illegal fix deform/pressure command argument: {}", arg[iarg + 1]);
} else if (strcmp(arg[iarg], "xy") == 0 ||
strcmp(arg[iarg], "xz") == 0 ||
strcmp(arg[iarg], "yz") == 0) {
if (strcmp(arg[iarg], "xy") == 0) index = 5;
else if (strcmp(arg[iarg], "xz") == 0) index = 4;
else if (strcmp(arg[iarg], "yz") == 0) index = 3;
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure", error);
if (strcmp(arg[iarg + 1], "pressure") == 0) {
if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure pressure", error);
set[index].style = PRESSURE;
if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) {
set_extra[index].ptarget = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
} else {
set_extra[index].pstr = utils::strdup(&arg[iarg + 2][2]);
set_extra[index].pvar_flag = 1;
}
set_extra[index].pgain = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
i += 4;
} else error->all(FLERR, "Illegal fix deform/pressure command: {}", arg[iarg + 1]);
} else if (strcmp(arg[iarg], "box") == 0) {
if (strcmp(arg[iarg + 1], "volume") == 0) {
set_box.style = VOLUME;
i += 2;
} else if (strcmp(arg[iarg + 1], "pressure") == 0) {
if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure pressure", error);
set_box.style = PRESSURE;
if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) {
set_extra[6].ptarget = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
} else {
set_extra[6].pstr = utils::strdup(&arg[iarg + 2][2]);
set_extra[6].pvar_flag = 1;
}
set_extra[6].pgain = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
i += 4;
} else error->all(FLERR, "Illegal fix deform/pressure command argument: {}", arg[iarg + 1]);
} else break;
}
// read options from end of input line
// shift arguments before reading
iarg = iarg_options_start;
options(i, narg - iarg, &arg[iarg]);
// repeat: setup dimflags used by other classes to check for volume-change conflicts
for (int i = 0; i < 6; i++)
if (set[i].style == NONE) dimflag[i] = 0;
else dimflag[i] = 1;
if (set_box.style != NONE) {
dimflag[0] = 1;
dimflag[1] = 1;
dimflag[2] = 1;
}
if (dimflag[0]) box_change |= BOX_CHANGE_X;
if (dimflag[1]) box_change |= BOX_CHANGE_Y;
if (dimflag[2]) box_change |= BOX_CHANGE_Z;
if (dimflag[3]) box_change |= BOX_CHANGE_YZ;
if (dimflag[4]) box_change |= BOX_CHANGE_XZ;
if (dimflag[5]) box_change |= BOX_CHANGE_XY;
// repeat: no tensile deformation on shrink-wrapped dims
// b/c shrink wrap will change box-length
for (int i = 0; i < 3; i++)
if (set_box.style && (domain->boundary[i][0] >= 2 || domain->boundary[i][1] >= 2))
error->all(FLERR, "Cannot use fix deform/pressure on a shrink-wrapped boundary");
// repeat: no tilt deformation on shrink-wrapped 2nd dim
// b/c shrink wrap will change tilt factor in domain::reset_box()
if (set[3].style && (domain->boundary[2][0] >= 2 || domain->boundary[2][1] >= 2))
error->all(FLERR, "Cannot use fix deform/pressure tilt on a shrink-wrapped 2nd dim");
if (set[4].style && (domain->boundary[2][0] >= 2 || domain->boundary[2][1] >= 2))
error->all(FLERR, "Cannot use fix deform/pressure tilt on a shrink-wrapped 2nd dim");
if (set[5].style && (domain->boundary[1][0] >= 2 || domain->boundary[1][1] >= 2))
error->all(FLERR, "Cannot use fix deform/pressure tilt on a shrink-wrapped 2nd dim");
// for VOLUME, setup links to other dims
// fixed, dynamic1, dynamic2
for (int i = 0; i < 3; i++) {
if (set[i].style != VOLUME) continue;
int other1 = (i + 1) % 3;
int other2 = (i + 2) % 3;
// Cannot use VOLUME option without at least one deformed dimension
if (set[other1].style == NONE || set[other1].style == VOLUME)
if (set[other2].style == NONE || set[other2].style == VOLUME)
error->all(FLERR, "Fix {} volume setting is invalid", style);
if (set[other1].style == NONE) {
set[i].substyle = ONE_FROM_ONE;
set[i].fixed = other1;
set[i].dynamic1 = other2;
} else if (set[other2].style == NONE) {
set[i].substyle = ONE_FROM_ONE;
set[i].fixed = other2;
set[i].dynamic1 = other1;
} else if (set[other1].style == VOLUME) {
set[i].substyle = TWO_FROM_ONE;
set[i].fixed = other1;
set[i].dynamic1 = other2;
} else if (set[other2].style == VOLUME) {
set[i].substyle = TWO_FROM_ONE;
set[i].fixed = other2;
set[i].dynamic1 = other1;
} else {
set[i].substyle = ONE_FROM_TWO;
set[i].dynamic1 = other1;
set[i].dynamic2 = other2;
}
}
// repeat: set varflag
for (int i = 0; i < 7; i++)
if (set_extra[i].pvar_flag) varflag = 1;
// repeat: reneighboring only forced if flips can occur due to shape changes
if ((!force_reneighbor) && flipflag && (set[3].style || set[4].style || set[5].style)) {
force_reneighbor = 1;
irregular = new Irregular(lmp);
}
// set initial values at time fix deform/pressure is issued
set_box.vol_initial = domain->xprd * domain->yprd * domain->zprd;
// populate coupled pressure controls
if (pcouple != NOCOUPLE) {
if (dimension == 2)
if (pcouple == XYZ || pcouple == XZ || pcouple == YZ)
error->all(FLERR, "Cannot couple Z dimension in fix deform/pressure in 2D");
int coupled_indices[3] = {0};
int j = -1;
if (pcouple == XYZ || pcouple == XY || pcouple == XZ)
coupled_indices[0] = 1;
if (pcouple == XYZ || pcouple == XY || pcouple == YZ)
coupled_indices[1] = 1;
if (pcouple == XYZ || pcouple == XZ || pcouple == YZ)
coupled_indices[2] = 1;
// Check coupled styles and find reference
for (int i = 0; i < 3; i++) {
if (coupled_indices[i]) {
set_extra[i].coupled_flag = 1;
if (set[i].style != PRESSURE && set[i].style != NONE)
error->all(FLERR, "Cannot couple non-pressure-controlled dimensions");
if (set[i].style == PRESSURE)
j = i;
}
}
if (j == -1)
error->all(FLERR, "Must specify deformation style for at least one coupled dimension");
// Copy or compare data for each coupled dimension
for (int i = 0; i < 3; i++) {
if (coupled_indices[i]) {
// Copy coupling information if dimension style is undefined
if (set[i].style == NONE) {
set[i].style = PRESSURE;
dimflag[i] = 1;
set_extra[i].pgain = set_extra[j].pgain;
if (set_extra[j].pvar_flag) {
set_extra[i].pstr = set_extra[j].pstr;
set_extra[i].pvar_flag = 1;
} else {
set_extra[i].ptarget = set_extra[j].ptarget;
}
} else {
// Check for incompatibilities in style
if (set[j].style != set[i].style && set[i].style != NONE)
error->all(FLERR, "Cannot couple dimensions with different control options");
if (set[j].style != PRESSURE) continue;
// If pressure controlled, check for incompatibilities in parameters
if (set_extra[i].pgain != set_extra[j].pgain || set_extra[i].pvar_flag != set_extra[j].pvar_flag ||
set_extra[i].ptarget != set_extra[j].ptarget)
error->all(FLERR, "Coupled dimensions must have identical gain parameters");
if (set_extra[j].pvar_flag)
if (strcmp(set_extra[i].pstr, set_extra[j].pstr) != 0)
error->all(FLERR, "Coupled dimensions must have the same target pressure");
}
}
}
}
// if vol/balance/p used, must have 2 free dimensions
if (vol_balance_flag) {
for (int i = 0; i < 3; i++) {
if (set[i].style != VOLUME) continue;
if (set[i].substyle != TWO_FROM_ONE)
error->all(FLERR, "Two dimensions must maintain constant volume to use the vol/balance/p option");
}
}
// set strain_flag
strain_flag = 0;
for (int i = 0; i < 6; i++)
if (set[i].style != NONE && set[i].style != VOLUME &&
set[i].style != PRESSURE && set[i].style != PMEAN)
strain_flag = 1;
// set pressure_flag
pressure_flag = 0;
for (int i = 0; i < 6; i++) {
if (set[i].style == PRESSURE || set[i].style == PMEAN) {
pressure_flag = 1;
if (set_extra[i].pgain <= 0.0)
error->all(FLERR, "Illegal fix deform/pressure gain constant, must be positive");
}
if (set_extra[i].coupled_flag) pressure_flag = 1;
}
if (set_box.style == PRESSURE) pressure_flag = 1;
if (vol_balance_flag) pressure_flag = 1;
// check conflict between constant volume/pressure
volume_flag = 0;
for (int i = 0; i < 3; i++)
if (set[i].style == VOLUME)
volume_flag = 1;
if (volume_flag)
for (int i = 0; i < 6; i++)
if (set[i].style == PMEAN)
error->all(FLERR, "Cannot use fix deform/pressure to assign constant volume and pressure");
// check conflicts between x,y,z styles and box
if (set_box.style)
for (int i = 0; i < 3; i++)
if (set[i].style == FINAL || set[i].style == DELTA || set[i].style == SCALE || set[i].style == PMEAN || set[i].style == VARIABLE)
error->all(FLERR, "Cannot use fix deform/pressure box parameter with x, y, or z styles other than vel, erate, trate, pressure, and wiggle");
// check pressure used for max rate and normalize error flag
if (!pressure_flag && max_h_rate != 0)
error->all(FLERR, "Can only assign a maximum strain rate using pressure-controlled dimensions");
if (!pressure_flag && normalize_pressure_flag)
error->all(FLERR, "Can only normalize error using pressure-controlled dimensions");
// Create pressure compute, if needed
pflag = 0;
tflag = 0;
if (pressure_flag) {
// 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
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} all temp",id_temp));
tflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
id_press = utils::strdup(std::string(id) + "_press");
modify->add_compute(fmt::format("{} all pressure {}",id_press, id_temp));
pflag = 1;
}
}
/* ---------------------------------------------------------------------- */
FixDeformPressure::~FixDeformPressure()
{
if (set_extra)
for (int i = 0; i < 7; i++)
delete[] set_extra[i].pstr;
delete[] set_extra;
delete[] set_box.hstr;
delete[] set_box.hratestr;
// delete temperature and pressure if fix created them
if (tflag) modify->delete_compute(id_temp);
if (pflag) modify->delete_compute(id_press);
delete [] id_temp;
delete [] id_press;
}
/* ---------------------------------------------------------------------- */
void FixDeformPressure::init()
{
FixDeform::init();
set_box.vol_start = domain->xprd * domain->yprd * domain->zprd;
// check optional variables for PRESSURE or PMEAN style
for (int i = 0; i < 7; i++) {
if (!set_extra[i].pvar_flag) continue;
set_extra[i].pvar = input->variable->find(set_extra[i].pstr);
if (set_extra[i].pvar < 0)
error->all(FLERR, "Variable name {} for fix deform/pressure does not exist", set_extra[i].pstr);
if (!input->variable->equalstyle(set_extra[i].pvar))
error->all(FLERR, "Variable {} for fix deform/pressure is invalid style", set_extra[i].pstr);
}
// Find pressure/temp computes if needed
if (pressure_flag) {
int icompute = modify->find_compute(id_temp);
if (icompute < 0) error->all(FLERR, "Temperature ID for fix deform/pressure does not exist");
temperature = modify->compute[icompute];
icompute = modify->find_compute(id_press);
if (icompute < 0) error->all(FLERR, "Pressure ID for fix deform/pressure does not exist");
pressure = modify->compute[icompute];
}
}
/* ----------------------------------------------------------------------
compute T,P if needed before integrator starts
------------------------------------------------------------------------- */
void FixDeformPressure::setup(int /*vflag*/)
{
// trigger virial computation on next timestep
if (pressure_flag) pressure->addstep(update->ntimestep+1);
}
/* ---------------------------------------------------------------------- */
void FixDeformPressure::end_of_step()
{
// wrap variable evaluations with clear/add
if (varflag) modify->clearstep_compute();
// set new box size for strain-based dims
if (strain_flag) FixDeform::apply_strain();
// set new box size for pressure-based dims
if (pressure_flag) {
temperature->compute_vector();
pressure->compute_vector();
pressure->compute_scalar();
for (int i = 0; i < 3; i++) {
if (!set_extra[i].saved) {
set_extra[i].saved = 1;
set_extra[i].prior_rate = 0.0;
set_extra[i].prior_pressure = pressure->vector[i];
}
}
apply_pressure();
}
// set new box size for VOLUME dims that are linked to other dims
// NOTE: still need to set h_rate for these dims
if (volume_flag) apply_volume();
// apply any final box scalings
if (set_box.style) apply_box();
// Save pressure/strain rate if required
if (pressure_flag) {
for (int i = 0; i < 3; i++) {
set_extra[i].prior_pressure = pressure->vector[i];
set_extra[i].prior_rate = ((set[i].hi_target - set[i].lo_target) /
(domain->boxhi[i] - domain->boxlo[i]) - 1.0) / update->dt;
}
}
if (varflag) modify->addstep_compute(update->ntimestep + nevery);
FixDeform::update_domain();
// trigger virial computation, if needed, on next timestep
if (pressure_flag)
pressure->addstep(update->ntimestep+1);
}
/* ----------------------------------------------------------------------
apply pressure controls
------------------------------------------------------------------------- */
void FixDeformPressure::apply_pressure()
{
// If variable pressure, calculate current target
for (int i = 0; i < 6; i++)
if (set[i].style == PRESSURE)
if (set_extra[i].pvar_flag)
set_extra[i].ptarget = input->variable->compute_equal(set_extra[i].pvar);
// Find current (possibly coupled/hydrostatic) pressure for X, Y, Z
double *tensor = pressure->vector;
double scalar = pressure->scalar;
double p_current[3];
if (pcouple == XYZ) {
double ave = THIRD * (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 {
if (set[0].style == PRESSURE) p_current[0] = tensor[0];
else if (set[0].style == PMEAN) p_current[0] = scalar;
if (set[1].style == PRESSURE) p_current[1] = tensor[1];
else if (set[1].style == PMEAN) p_current[1] = scalar;
if (set[2].style == PRESSURE) p_current[2] = tensor[2];
else if (set[2].style == PMEAN) p_current[2] = scalar;
}
for (int i = 0; i < 3; i++) {
if (set[i].style != PRESSURE && set[i].style != PMEAN) continue;
h_rate[i] = set_extra[i].pgain * (p_current[i] - set_extra[i].ptarget);
if (normalize_pressure_flag) {
if (set_extra[i].ptarget == 0) {
if (max_h_rate == 0) {
error->all(FLERR, "Cannot normalize error for zero pressure without defining a max rate");
} else h_rate[i] = max_h_rate * h_rate[i] / fabs(h_rate[i]);
} else h_rate[i] /= fabs(set_extra[i].ptarget);
}
if (max_h_rate != 0)
if (fabs(h_rate[i]) > max_h_rate)
h_rate[i] = max_h_rate * h_rate[i] / fabs(h_rate[i]);
h_ratelo[i] = -0.5 * h_rate[i];
double offset = 0.5 * (domain->boxhi[i] - domain->boxlo[i]) * (1.0 + update->dt * h_rate[i]);
set[i].lo_target = 0.5 * (set[i].lo_start + set[i].hi_start) - offset;
set[i].hi_target = 0.5 * (set[i].lo_start + set[i].hi_start) + offset;
}
for (int i = 3; i < 6; i++) {
if (set[i].style != PRESSURE) continue;
double L, tilt, pcurrent;
if (i == 3) {
L = domain->zprd;
tilt = domain->yz;
pcurrent = tensor[5];
} else if (i == 4) {
L = domain->zprd;
tilt = domain->xz + update->dt;
pcurrent = tensor[4];
} else {
L = domain->yprd;
tilt = domain->xy;
pcurrent = tensor[3];
}
h_rate[i] = L * set_extra[i].pgain * (pcurrent - set_extra[i].ptarget);
if (normalize_pressure_flag) {
if (set_extra[i].ptarget == 0) {
if (max_h_rate == 0) {
error->all(FLERR, "Cannot normalize error for zero pressure without defining a max rate");
} else h_rate[i] = max_h_rate * h_rate[i] / fabs(h_rate[i]);
} else h_rate[i] /= fabs(set_extra[i].ptarget);
}
if (max_h_rate != 0)
if (fabs(h_rate[i]) > max_h_rate)
h_rate[i] = max_h_rate * h_rate[i] / fabs(h_rate[i]);
set[i].tilt_target = tilt + update->dt * h_rate[i];
}
}
/* ----------------------------------------------------------------------
apply volume controls
------------------------------------------------------------------------- */
void FixDeformPressure::apply_volume()
{
double e1, e2;
int linked_pressure = 0;
for (int i = 0; i < 3; i++) {
if (set[i].style != VOLUME) continue;
int dynamic1 = set[i].dynamic1;
int dynamic2 = set[i].dynamic2;
int fixed = set[i].fixed;
double v0 = set[i].vol_start;
double shift;
if (set[i].substyle == ONE_FROM_ONE) {
shift = 0.5 * (v0 / (set[dynamic1].hi_target - set[dynamic1].lo_target) /
(set[fixed].hi_start-set[fixed].lo_start));
} else if (set[i].substyle == ONE_FROM_TWO) {
shift = 0.5 * (v0 / (set[dynamic1].hi_target - set[dynamic1].lo_target) /
(set[dynamic2].hi_target - set[dynamic2].lo_target));
} else if (set[i].substyle == TWO_FROM_ONE) {
if (!vol_balance_flag) {
shift = 0.5 * sqrt(v0 * (set[i].hi_start - set[i].lo_start) /
(set[dynamic1].hi_target - set[dynamic1].lo_target) /
(set[fixed].hi_start - set[fixed].lo_start));
} else {
double dt = update->dt;
double e1i = set_extra[i].prior_rate;
double e2i = set_extra[fixed].prior_rate;
double L1i = domain->boxhi[i] - domain->boxlo[i];
double L2i = domain->boxhi[fixed] - domain->boxlo[fixed];
double L3i = domain->boxhi[dynamic1] - domain->boxlo[dynamic1];
double L3 = (set[dynamic1].hi_target - set[dynamic1].lo_target);
double Vi = L1i * L2i * L3i;
double V = L3 * L1i * L2i;
double e3 = (L3 / L3i - 1.0) / dt;
double p1 = pressure->vector[i];
double p2 = pressure->vector[fixed];
double p1i = set_extra[i].prior_pressure;
double p2i = set_extra[fixed].prior_pressure;
double denominator;
if (e3 == 0) {
e1 = 0.0;
e2 = 0.0;
shift = 0.5 * L1i;
} else if (e1i == 0 || e2i == 0 || (p2 == p2i && p1 == p1i)) {
// If no prior strain or no change in pressure (initial step) just scale shift by relative box lengths
shift = 0.5 * sqrt(v0 * L1i / L3 / L2i);
} else {
if (!linked_pressure) {
// Calculate first strain rate by expanding stress to linear order, p1(t+dt) = p2(t+dt)
// Calculate second strain rate to preserve volume
denominator = p2 - p2i + e2i * ((p1 - p1i) / e1i);
if (denominator != 0.0 && e1i != 0.0) {
e1 = (((p2 - p2i) * (Vi - V) / (V * dt)) - e2i * (p1 - p2)) / denominator;
} else {
e1 = e2i;
}
e2 = (Vi - V * (1 + e1 * dt)) / (V * (1 + e1 * dt) * dt);
// If strain rate exceeds limit in either dimension, cap it at the maximum compatible rate
if (max_h_rate != 0) {
if ((fabs(e1) > max_h_rate) || (fabs(e2) > max_h_rate)) {
if (fabs(e1) > fabs(e2))
adjust_linked_rates(e1, e2, e3, Vi, V);
else
adjust_linked_rates(e2, e1, e3, Vi, V);
}
}
shift = 0.5 * L1i * (1.0 + e1 * dt);
linked_pressure = 1;
} else {
// Already calculated value of e2
shift = 0.5 * L1i * (1.0 + e2 * dt);
}
}
}
}
h_rate[i] = (2.0 * shift / (domain->boxhi[i] - domain->boxlo[i]) - 1.0) / update->dt;
h_ratelo[i] = -0.5 * h_rate[i];
set[i].lo_target = 0.5 * (set[i].lo_start + set[i].hi_start) - shift;
set[i].hi_target = 0.5 * (set[i].lo_start + set[i].hi_start) + shift;
}
}
/* ----------------------------------------------------------------------
Rescale volume preserving strain rates to enforce max rate
------------------------------------------------------------------------- */
void FixDeformPressure::adjust_linked_rates(double &e_larger, double &e_smaller, double e3, double Vi, double V)
{
double dt = update->dt;
double e_lim_positive = (Vi - V * (1 + max_h_rate * dt)) / (V * (1 + max_h_rate * dt) * dt);
double e_lim_negative = (Vi - V * (1 - max_h_rate * dt)) / (V * (1 - max_h_rate * dt) * dt);
if ((e_larger * e3) >= 0) {
if (e_larger > 0.0) {
// Same sign as primary strain rate, cap third dimension
e_smaller = -max_h_rate;
e_larger = e_lim_negative;
} else {
e_smaller = max_h_rate;
e_larger = e_lim_positive;
}
} else {
// Opposite sign, set to maxrate.
if (e_larger > 0.0) {
e_larger = max_h_rate;
e_smaller = e_lim_positive;
} else {
e_larger = -max_h_rate;
e_smaller = e_lim_negative;
}
}
}
/* ----------------------------------------------------------------------
apply box controls
------------------------------------------------------------------------- */
void FixDeformPressure::apply_box()
{
int i;
double scale, shift;
double v_rate;
if (set_box.style == VOLUME) {
double v0 = set_box.vol_start;
double v = 1.0;
for (i = 0; i < 3; i++)
v *= (set[i].hi_target - set[i].lo_target);
scale = std::pow(v0 / v, THIRD);
for (i = 0; i < 3; i++) {
shift = 0.5 * (set[i].hi_target - set[i].lo_target) * scale;
set[i].lo_target = 0.5 * (set[i].lo_start + set[i].hi_start) - shift;
set[i].hi_target = 0.5 * (set[i].lo_start + set[i].hi_start) + shift;
// Recalculate h_rate
h_rate[i] = (set[i].hi_target - set[i].lo_target) / (domain->boxhi[i] - domain->boxlo[i]) - 1.0;
h_rate[i] /= update->dt;
h_ratelo[i] = -0.5 * h_rate[i];
}
} else if (set_box.style == PRESSURE) {
// If variable pressure, calculate current target
if (set_extra[6].pvar_flag)
set_extra[6].ptarget = input->variable->compute_equal(set_extra[6].pvar);
v_rate = set_extra[6].pgain * (pressure->scalar - set_extra[6].ptarget);
if (normalize_pressure_flag) {
if (set_extra[6].ptarget == 0) {
if (max_h_rate == 0) {
error->all(FLERR, "Cannot normalize error for zero pressure without defining a max rate");
} else v_rate = max_h_rate * v_rate / fabs(v_rate);
} else v_rate /= fabs(set_extra[6].ptarget);
}
if (max_h_rate != 0)
if (fabs(v_rate) > max_h_rate)
v_rate = max_h_rate * v_rate / fabs(v_rate);
set_extra[6].cumulative_strain += update->dt * v_rate;
scale = (1.0 + set_extra[6].cumulative_strain);
for (i = 0; i < 3; i++) {
shift = 0.5 * (set[i].hi_target - set[i].lo_target) * scale;
set[i].lo_target = 0.5 * (set[i].lo_start + set[i].hi_start) - shift;
set[i].hi_target = 0.5 * (set[i].lo_start + set[i].hi_start) + shift;
// Recalculate h_rate
h_rate[i] = (set[i].hi_target - set[i].lo_target) / (domain->boxhi[i] - domain->boxlo[i]) - 1.0;
h_rate[i] /= update->dt;
h_ratelo[i] = -0.5 * h_rate[i];
}
}
}
/* ----------------------------------------------------------------------
write Set data to restart file
------------------------------------------------------------------------- */
void FixDeformPressure::write_restart(FILE *fp)
{
if (comm->me == 0) {
int size = 9 * sizeof(double) + 7 * sizeof(Set) + 7 * sizeof(SetExtra);
fwrite(&size, sizeof(int), 1, fp);
fwrite(set, sizeof(Set), 6, fp);
fwrite(&set_box, sizeof(Set), 1, fp);
fwrite(set_extra, sizeof(SetExtra), 7, fp);
}
}
/* ----------------------------------------------------------------------
use selected state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixDeformPressure::restart(char *buf)
{
int n = 0;
auto list = (double *) buf;
for (int i = 0; i < 6; i++)
h_rate[i] = list[n++];
for (int i = 0; i < 3; i++)
h_ratelo[i] = list[n++];
n = n * sizeof(double);
int samestyle = 1;
Set *set_restart = (Set *) &buf[n];
for (int i = 0; i < 6; ++i) {
// restore data from initial state
set[i].lo_initial = set_restart[i].lo_initial;
set[i].hi_initial = set_restart[i].hi_initial;
set[i].vol_initial = set_restart[i].vol_initial;
set[i].tilt_initial = set_restart[i].tilt_initial;
// check if style settings are consistent (should do the whole set?)
if (set[i].style != set_restart[i].style)
samestyle = 0;
if (set[i].substyle != set_restart[i].substyle)
samestyle = 0;
}
n += 6 * sizeof(Set);
// Only restore relevant box variables & check consistency
Set set_box_restart;
memcpy(&set_box_restart, (Set *) &buf[n], sizeof(Set));
set_box.vol_initial = set_box_restart.vol_initial;
if (set_box.style != set_box_restart.style)
samestyle = 0;
if (!samestyle)
error->all(FLERR, "Fix deform/pressure settings not consistent with restart");
n += sizeof(Set);
SetExtra *set_extra_restart = (SetExtra *) &buf[n];
for (int i = 0; i < 7; ++i) {
set_extra[i].saved = set_extra_restart[i].saved;
set_extra[i].prior_rate = set_extra_restart[i].prior_rate;
set_extra[i].prior_pressure = set_extra_restart[i].prior_pressure;
set_extra[i].cumulative_strain = set_extra_restart[i].cumulative_strain;
}
}
/* ---------------------------------------------------------------------- */
void FixDeformPressure::options(int i, int narg, char **arg)
{
pcouple = NOCOUPLE;
max_h_rate = 0.0;
vol_balance_flag = 0;
normalize_pressure_flag = 0;
// parse only options not handled by parent class
int iarg, nskip;
while (i < leftover_iarg.size()) {
iarg = leftover_iarg[i];
if (strcmp(arg[iarg], "couple") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure couple", error);
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 = NOCOUPLE;
else error->all(FLERR, "Illegal fix deform/pressure couple command: {}", arg[iarg + 1]);
i += 2;
} else if (strcmp(arg[iarg], "max/rate") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure max/rate", error);
max_h_rate = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
if (max_h_rate <= 0.0)
error->all(FLERR, "Maximum strain rate must be a positive, non-zero value");
i += 2;
} else if (strcmp(arg[iarg], "normalize/pressure") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure normalize/pressure", error);
normalize_pressure_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
i += 2;
} else if (strcmp(arg[iarg], "vol/balance/p") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure vol/balance/p", error);
vol_balance_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
i += 2;
} else error->all(FLERR, "Illegal fix deform/pressure command: {}", arg[iarg]);
}
}
/* ---------------------------------------------------------------------- */
int FixDeformPressure::modify_param(int narg, char **arg)
{
if (strcmp(arg[0], "temp") == 0) {
if (narg < 2) error->all(FLERR, "Illegal fix_modify command");
if (tflag) {
modify->delete_compute(id_temp);
tflag = 0;
}
delete[] id_temp;
id_temp = utils::strdup(arg[1]);
temperature = modify->get_compute_by_id(arg[1]);
if (!temperature)
error->all(FLERR, "Could not find fix_modify temperature compute ID: ", arg[1]);
if (temperature->tempflag == 0)
error->all(FLERR, "Fix_modify temperature compute {} does not compute temperature", arg[1]);
if (temperature->igroup != 0 && comm->me == 0)
error->warning(FLERR, "Temperature compute {} for fix {} is not for group all: {}",
arg[1], style, group->names[temperature->igroup]);
// reset id_temp of pressure to new temperature ID
auto icompute = modify->get_compute_by_id(id_press);
if (!icompute)
error->all(FLERR, "Pressure compute ID {} for fix {} does not exist", id_press, style);
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 (pflag) {
modify->delete_compute(id_press);
pflag = 0;
}
delete[] id_press;
id_press = utils::strdup(arg[1]);
pressure = modify->get_compute_by_id(arg[1]);
if (!pressure) error->all(FLERR, "Could not find fix_modify pressure compute ID: {}", arg[1]);
if (pressure->pressflag == 0)
error->all(FLERR, "Fix_modify pressure compute {} does not compute pressure", arg[1]);
return 2;
}
return 0;
}

74
src/EXTRA-FIX/fix_deform_pressure.h Normal file
View File

@ -0,0 +1,74 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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
// clang-format off
FixStyle(deform/pressure,FixDeformPressure);
// clang-format on
#else
#ifndef LMP_FIX_DEFORM_PRESSURE_H
#define LMP_FIX_DEFORM_PRESSURE_H
#include "fix_deform.h"
namespace LAMMPS_NS {
class FixDeformPressure : public FixDeform {
public:
FixDeformPressure(class LAMMPS *, int, char **);
~FixDeformPressure() override;
void init() override;
void setup(int) override;
void end_of_step() override;
void write_restart(FILE *) override;
void restart(char *buf) override;
int modify_param(int, char **) override;
protected:
int pcouple, dimension;
double max_h_rate;
int strain_flag; // 1 if strain-based option is used, 0 if not
int pressure_flag; // 1 if pressure tensor used, 0 if not
int volume_flag; // 1 if VOLUME option is used, 0 if not
int normalize_pressure_flag; // 1 if normalize pressure deviation by target
int vol_balance_flag; // 1 if pressures balanced when maintaining const vol
char *id_temp, *id_press;
class Compute *temperature, *pressure;
int tflag, pflag;
struct SetExtra {
double ptarget, pgain;
double prior_pressure, prior_rate;
double cumulative_strain;
int saved;
char *pstr;
int pvar, pvar_flag;
int coupled_flag;
};
SetExtra *set_extra;
Set set_box;
void options(int, int, char **);
void apply_volume() override;
void apply_pressure();
void apply_box();
void couple();
void adjust_linked_rates(double&, double&, double, double, double);
};
} // namespace LAMMPS_NS
#endif
#endif

11
src/KOKKOS/compute_reaxff_atom_kokkos.cpp
View File

@ -67,10 +67,10 @@ void ComputeReaxFFAtomKokkos<DeviceType>::init()
template<class DeviceType>
void ComputeReaxFFAtomKokkos<DeviceType>::compute_bonds()
{
if (atom->nlocal > nlocal) {
if (atom->nmax > nmax) {
memory->destroy(array_atom);
nlocal = atom->nlocal;
memory->create(array_atom, nlocal, 3, "reaxff/atom:array_atom");
nmax = atom->nmax;
memory->create(array_atom, nmax, 3, "reaxff/atom:array_atom");
}
// retrieve bond information from kokkos pair style. the data potentially
@ -85,6 +85,7 @@ void ComputeReaxFFAtomKokkos<DeviceType>::compute_bonds()
else
host_pair()->FindBond(maxnumbonds, groupbit);
const int nlocal = atom->nlocal;
nbuf = ((store_bonds ? maxnumbonds*2 : 0) + 3)*nlocal;
if (!buf || ((int)k_buf.extent(0) < nbuf)) {
@ -135,6 +136,7 @@ void ComputeReaxFFAtomKokkos<DeviceType>::compute_local()
int b = 0;
int j = 0;
auto tag = atom->tag;
const int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; ++i) {
const int numbonds = static_cast<int>(buf[j+2]);
@ -161,6 +163,7 @@ void ComputeReaxFFAtomKokkos<DeviceType>::compute_peratom()
compute_bonds();
// extract peratom bond information from buffer
const int nlocal = atom->nlocal;
int j = 0;
for (int i = 0; i < nlocal; ++i) {
@ -180,7 +183,7 @@ void ComputeReaxFFAtomKokkos<DeviceType>::compute_peratom()
template<class DeviceType>
double ComputeReaxFFAtomKokkos<DeviceType>::memory_usage()
{
double bytes = (double)(nlocal*3) * sizeof(double);
double bytes = (double)(nmax*3) * sizeof(double);
if (store_bonds)
bytes += (double)(nbonds*3) * sizeof(double);
bytes += (double)(nbuf > 0 ? nbuf * sizeof(double) : 0);

14
src/KOKKOS/fix_deform_kokkos.cpp
View File

@ -120,11 +120,11 @@ void FixDeformKokkos::end_of_step()
} else if (set[i].style == WIGGLE) {
double delt = (update->ntimestep - update->beginstep) * update->dt;
set[i].lo_target = set[i].lo_start -
0.5*set[i].amplitude * sin(TWOPI*delt/set[i].tperiod);
0.5*set[i].amplitude * sin(MY_2PI*delt/set[i].tperiod);
set[i].hi_target = set[i].hi_start +
0.5*set[i].amplitude * sin(TWOPI*delt/set[i].tperiod);
h_rate[i] = TWOPI/set[i].tperiod * set[i].amplitude *
cos(TWOPI*delt/set[i].tperiod);
0.5*set[i].amplitude * sin(MY_2PI*delt/set[i].tperiod);
h_rate[i] = MY_2PI/set[i].tperiod * set[i].amplitude *
cos(MY_2PI*delt/set[i].tperiod);
h_ratelo[i] = -0.5*h_rate[i];
} else if (set[i].style == VARIABLE) {
double del = input->variable->compute_equal(set[i].hvar);
@ -212,9 +212,9 @@ void FixDeformKokkos::end_of_step()
} else if (set[i].style == WIGGLE) {
double delt = (update->ntimestep - update->beginstep) * update->dt;
set[i].tilt_target = set[i].tilt_start +
set[i].amplitude * sin(TWOPI*delt/set[i].tperiod);
h_rate[i] = TWOPI/set[i].tperiod * set[i].amplitude *
cos(TWOPI*delt/set[i].tperiod);
set[i].amplitude * sin(MY_2PI*delt/set[i].tperiod);
h_rate[i] = MY_2PI/set[i].tperiod * set[i].amplitude *
cos(MY_2PI*delt/set[i].tperiod);
} else if (set[i].style == VARIABLE) {
double delta_tilt = input->variable->compute_equal(set[i].hvar);
set[i].tilt_target = set[i].tilt_start + delta_tilt;

26
src/REAXFF/compute_reaxff_atom.cpp
View File

@ -43,7 +43,7 @@ ComputeReaxFFAtom::ComputeReaxFFAtom(LAMMPS *lmp, int narg, char **arg) :
// initialize output
nlocal = -1;
nmax = -1;
nbonds = 0;
prev_nbonds = -1;
@ -162,20 +162,22 @@ void ComputeReaxFFAtom::compute_bonds()
{
invoked_bonds = update->ntimestep;
if (atom->nlocal > nlocal) {
if (atom->nmax > nmax) {
memory->destroy(abo);
memory->destroy(neighid);
memory->destroy(bondcount);
memory->destroy(array_atom);
nlocal = atom->nlocal;
nmax = atom->nmax;
if (store_bonds) {
memory->create(abo, nlocal, MAXREAXBOND, "reaxff/atom:abo");
memory->create(neighid, nlocal, MAXREAXBOND, "reaxff/atom:neighid");
memory->create(abo, nmax, MAXREAXBOND, "reaxff/atom:abo");
memory->create(neighid, nmax, MAXREAXBOND, "reaxff/atom:neighid");
}
memory->create(bondcount, nlocal, "reaxff/atom:bondcount");
memory->create(array_atom, nlocal, 3, "reaxff/atom:array_atom");
memory->create(bondcount, nmax, "reaxff/atom:bondcount");
memory->create(array_atom, nmax, 3, "reaxff/atom:array_atom");
}
const int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
bondcount[i] = 0;
for (int j = 0; store_bonds && j < MAXREAXBOND; j++) {
@ -208,6 +210,8 @@ void ComputeReaxFFAtom::compute_local()
int b = 0;
const int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; ++i) {
const int numbonds = bondcount[i];
@ -230,6 +234,8 @@ void ComputeReaxFFAtom::compute_peratom()
compute_bonds();
}
const int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; ++i) {
auto ptr = array_atom[i];
ptr[0] = reaxff->api->workspace->total_bond_order[i];
@ -244,10 +250,10 @@ void ComputeReaxFFAtom::compute_peratom()
double ComputeReaxFFAtom::memory_usage()
{
double bytes = (double)(nlocal*3) * sizeof(double);
bytes += (double)(nlocal) * sizeof(int);
double bytes = (double)(nmax*3) * sizeof(double);
bytes += (double)(nmax) * sizeof(int);
if (store_bonds) {
bytes += (double)(2*nlocal*MAXREAXBOND) * sizeof(double);
bytes += (double)(2*nmax*MAXREAXBOND) * sizeof(double);
bytes += (double)(nbonds*3) * sizeof(double);
}
return bytes;

2
src/REAXFF/compute_reaxff_atom.h
View File

@ -40,7 +40,7 @@ class ComputeReaxFFAtom : public Compute {
protected:
bigint invoked_bonds; // last timestep on which compute_bonds() was invoked
int nlocal;
int nmax;
int nbonds;
int prev_nbonds;
int nsub;

718
src/fix_deform.cpp
View File

File diff suppressed because it is too large Load Diff

15
src/fix_deform.h
View File

@ -29,14 +29,17 @@ class FixDeform : public Fix {
int remapflag; // whether x,v are remapped across PBC
int dimflag[6]; // which dims are deformed
enum { NONE, FINAL, DELTA, SCALE, VEL, ERATE, TRATE, VOLUME, WIGGLE, VARIABLE, PRESSURE, PMEAN };
enum { ONE_FROM_ONE, ONE_FROM_TWO, TWO_FROM_ONE };
FixDeform(class LAMMPS *, int, char **);
~FixDeform() override;
int setmask() override;
void init() override;
void pre_exchange() override;
void end_of_step() override;
void write_restart(FILE *) override;
void restart(char *buf) override;
void virtual write_restart(FILE *) override;
void virtual restart(char *buf) override;
double memory_usage() override;
protected:
@ -48,8 +51,6 @@ class FixDeform : public Fix {
std::vector<Fix *> rfix; // pointers to rigid fixes
class Irregular *irregular; // for migrating atoms after box flips
double TWOPI;
struct Set {
int style, substyle;
double flo, fhi, ftilt;
@ -67,7 +68,13 @@ class FixDeform : public Fix {
};
Set *set;
std::vector<int> leftover_iarg;
int iarg_options_start;
void options(int, char **);
void virtual apply_volume();
void apply_strain();
void update_domain();
};
} // namespace LAMMPS_NS