diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index e89e302673a..ea50e68cdd2 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -61,6 +61,7 @@ OPT. * :doc:`controller ` * :doc:`damping/cundall ` * :doc:`deform (k) ` + * :doc:`deform/pressure ` * :doc:`deposit ` * :doc:`dpd/energy (k) ` * :doc:`drag ` diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 69a72124875..d03cab46876 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -226,6 +226,7 @@ accelerated styles exist. * :doc:`controller ` - apply control loop feedback mechanism * :doc:`damping/cundall ` - Cundall non-viscous damping for granular simulations * :doc:`deform ` - change the simulation box size/shape +* :doc:`deform/pressure ` - change the simulation box size/shape with additional loading conditions * :doc:`deposit ` - add new atoms above a surface * :doc:`dpd/energy ` - constant energy dissipative particle dynamics * :doc:`drag ` - drag atoms towards a defined coordinate diff --git a/doc/src/fix_deform.rst b/doc/src/fix_deform.rst index ee010f5645c..9146b987c86 100644 --- a/doc/src/fix_deform.rst +++ b/doc/src/fix_deform.rst @@ -4,6 +4,9 @@ fix deform command ================== +:doc:`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 ` 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 ` 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 ` 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 ` 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 ` command) + couple pressure values of various dimensions + *vol/balance/p* value = *yes* or *no* (ONLY available in :doc:`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 ` command) + rate = maximum strain rate for pressure control + *normalize/pressure* value = *yes* or *no* (ONLY available in :doc:`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 ` 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. - -This fix can be used to perform non-equilibrium MD (NEMD) simulations -of a continuously strained system. See the :doc:`fix nvt/sllod ` and :doc:`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 ` and its :doc:`fix commands `. +adjusted independently and simultaneously. + +The fix deform command allows use of all the arguments listed above, +except those flagged as available ONLY for the :doc:`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 ` 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 ` and :doc:`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 ` and its :doc:`fix +commands `. 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 ` or :doc:`fix nph ` 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 ` or :doc:`fix nph ` command. The size and shape of the simulation box at the beginning of the -simulation run were either specified by the -:doc:`create_box ` or :doc:`read_data ` or -:doc:`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 `, :doc:`read data `, and :doc:`read_restart ` commands -specify whether the simulation box is orthogonal or non-orthogonal -(triclinic) and explain the meaning of the xy,xz,yz tilt factors. 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 +` or :doc:`read_data ` or :doc:`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 `, :doc:`read data +`, and :doc:`read_restart ` commands specify +whether the simulation box is orthogonal or non-orthogonal (triclinic) +and explain the meaning of the xy,xz,yz tilt factors. 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 ` -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 ` 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 +variable with *name1* must be an :doc:`equal-style 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 +` 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 ` 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 +` 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 ` or :doc:`fix lavgevin ` - (with a bias provided by :doc:`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 ` or :doc:`fix lavgevin + ` (with a bias provided by :doc:`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 `, :doc:`compute temp/deform `, and :doc:`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 `, :doc:`compute temp/deform + `, and :doc:`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 ` - 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 ` 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 ` 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 ` 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 `. + include a background fluid and thermostat the fluid via :doc:`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 `, which allows the fix to be properly continue +This fix will restore the initial box settings from :doc:`binary +restart files `, which allows the fix to be properly continue deformation, when using the start/stop options of the :doc:`run ` command. None of the :doc:`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 ` command. -You cannot apply xy, yz, or xz deformations to a second dimension (y in -xy) that is shrink-wrapped via the :doc:`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 ` +command. Related commands """""""""""""""" +:doc:`fix deform/pressure `, :doc:`change_box ` Default diff --git a/doc/src/fix_deform_pressure.rst b/doc/src/fix_deform_pressure.rst new file mode 100644 index 00000000000..f85ad37238b --- /dev/null +++ b/doc/src/fix_deform_pressure.rst @@ -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 ` 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 ` 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 ` 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 ` 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 ` +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 ` 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 `. The target pressure +accepts either a constant numeric value or a LAMMPS :ref:`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 `. 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 ` +option. Akin to :doc:`fix nh `, pressures in different dimensions +can be coupled using the :ref:`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 ` 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 +`. 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 +`. 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 ` 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) `. + +---------- + +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 ` and :doc:`compute pressure +` commands for details. Note that the IDs of the +new computes are the fix-ID + underscore + "temp" or fix_ID ++ underscore + "press", 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 ` 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 ` command or print this +temperature or pressure during thermodynamic output via the +:doc:`thermo_style custom ` command using the +appropriate compute-ID. It also means that changing attributes of +*thermo_temp* or *thermo_press* will have no effect on this fix. + +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +This fix will restore the initial box settings from :doc:`binary +restart files `, which allows the fix to be properly continue +deformation, when using the start/stop options of the :doc:`run ` +command. No global or per-atom quantities are stored by this fix for +access by various :doc:`output commands `. + +If any pressure controls are used, the :doc:`fix_modify ` +*temp* and *press* options are supported by this fix, unlike in +:doc:`fix deform `. You can use them to assign a +:doc:`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 ` command. See the +:doc:`run ` command for details of how to do this. + +This fix is not invoked during :doc:`energy minimization `. + +Restrictions +"""""""""""" + +You cannot apply x, y, or z deformations to a dimension that is +shrink-wrapped via the :doc:`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 ` +command. + +Related commands +"""""""""""""""" + +:doc:`fix deform `, :doc:`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 diff --git a/src/.gitignore b/src/.gitignore index 1e4c5b9ddb1..9e2c36bdd00 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -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 diff --git a/src/EXTRA-FIX/fix_deform_pressure.cpp b/src/EXTRA-FIX/fix_deform_pressure.cpp new file mode 100644 index 00000000000..7f15870ef12 --- /dev/null +++ b/src/EXTRA-FIX/fix_deform_pressure.cpp @@ -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 +#include +#include + +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; +} diff --git a/src/EXTRA-FIX/fix_deform_pressure.h b/src/EXTRA-FIX/fix_deform_pressure.h new file mode 100644 index 00000000000..a52bb01c042 --- /dev/null +++ b/src/EXTRA-FIX/fix_deform_pressure.h @@ -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 diff --git a/src/KOKKOS/fix_deform_kokkos.cpp b/src/KOKKOS/fix_deform_kokkos.cpp index d49e3359866..97c78bb1a5a 100644 --- a/src/KOKKOS/fix_deform_kokkos.cpp +++ b/src/KOKKOS/fix_deform_kokkos.cpp @@ -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; diff --git a/src/fix_deform.cpp b/src/fix_deform.cpp index 02aaae59409..9ee72594826 100644 --- a/src/fix_deform.cpp +++ b/src/fix_deform.cpp @@ -34,171 +34,205 @@ #include #include +#include +#include +#include using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; -enum{NONE=0,FINAL,DELTA,SCALE,VEL,ERATE,TRATE,VOLUME,WIGGLE,VARIABLE}; -enum{ONE_FROM_ONE,ONE_FROM_TWO,TWO_FROM_ONE}; - /* ---------------------------------------------------------------------- */ FixDeform::FixDeform(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), irregular(nullptr), set(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal fix deform command"); + const std::string thiscmd = fmt::format("fix {}", style); + if (narg < 4) utils::missing_cmd_args(FLERR, thiscmd, error); no_change_box = 1; restart_global = 1; pre_exchange_migrate = 1; - nevery = utils::inumeric(FLERR,arg[3],false,lmp); - if (nevery <= 0) error->all(FLERR,"Illegal fix deform command"); + nevery = utils::inumeric(FLERR, arg[3], false, lmp); + if (nevery <= 0) error->all(FLERR, "Fix {} Nevery must be > 0", style); + + // arguments for child classes + + std::unordered_set child_parameters; + std::unordered_map child_styles; + int nskip; + if (utils::strmatch(style, "^deform/pressure")) { + child_parameters.insert("box"); + child_styles.insert({{"pressure", 4}, {"pressure/mean", 4}, {"volume", 2}}); + } // set defaults set = new Set[6]; - memset(set,0,6*sizeof(Set)); + memset(set, 0, 6 * sizeof(Set)); - // parse arguments + // parse all parameter/style arguments for this parent and also child classes + // for child classes, simply store them in leftover_iarg and skip over them triclinic = domain->triclinic; int index; int iarg = 4; + while (iarg < narg) { - if (strcmp(arg[iarg],"x") == 0 || - strcmp(arg[iarg],"y") == 0 || - strcmp(arg[iarg],"z") == 0) { + 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 (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) error->all(FLERR,"Illegal fix deform command"); - if (strcmp(arg[iarg+1],"final") == 0) { - if (iarg+4 > narg) error->all(FLERR,"Illegal fix deform command"); + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, thiscmd, error); + if (strcmp(arg[iarg + 1], "final") == 0) { + if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, thiscmd + " final", error); set[index].style = FINAL; - set[index].flo = utils::numeric(FLERR,arg[iarg+2],false,lmp); - set[index].fhi = utils::numeric(FLERR,arg[iarg+3],false,lmp); + set[index].flo = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + set[index].fhi = utils::numeric(FLERR, arg[iarg + 3], false, lmp); iarg += 4; - } else if (strcmp(arg[iarg+1],"delta") == 0) { - if (iarg+4 > narg) error->all(FLERR,"Illegal fix deform command"); + } else if (strcmp(arg[iarg + 1], "delta") == 0) { + if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, thiscmd + " delta", error); set[index].style = DELTA; - set[index].dlo = utils::numeric(FLERR,arg[iarg+2],false,lmp); - set[index].dhi = utils::numeric(FLERR,arg[iarg+3],false,lmp); + set[index].dlo = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + set[index].dhi = utils::numeric(FLERR, arg[iarg + 3], false, lmp); iarg += 4; - } else if (strcmp(arg[iarg+1],"scale") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix deform command"); + } else if (strcmp(arg[iarg + 1], "scale") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, thiscmd + " scale", error); set[index].style = SCALE; - set[index].scale = utils::numeric(FLERR,arg[iarg+2],false,lmp); + set[index].scale = utils::numeric(FLERR, arg[iarg + 2], false, lmp); iarg += 3; - } else if (strcmp(arg[iarg+1],"vel") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix deform command"); + } else if (strcmp(arg[iarg + 1], "vel") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, thiscmd + " vel", error); set[index].style = VEL; - set[index].vel = utils::numeric(FLERR,arg[iarg+2],false,lmp); + set[index].vel = utils::numeric(FLERR, arg[iarg + 2], false, lmp); iarg += 3; - } else if (strcmp(arg[iarg+1],"erate") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix deform command"); + } else if (strcmp(arg[iarg + 1], "erate") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, thiscmd + " erate", error); set[index].style = ERATE; - set[index].rate = utils::numeric(FLERR,arg[iarg+2],false,lmp); + set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp); iarg += 3; - } else if (strcmp(arg[iarg+1],"trate") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix deform command"); + } else if (strcmp(arg[iarg + 1], "trate") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, thiscmd + " trate", error); set[index].style = TRATE; - set[index].rate = utils::numeric(FLERR,arg[iarg+2],false,lmp); + set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp); iarg += 3; - } else if (strcmp(arg[iarg+1],"volume") == 0) { + } else if (strcmp(arg[iarg + 1], "volume") == 0) { set[index].style = VOLUME; iarg += 2; - } else if (strcmp(arg[iarg+1],"wiggle") == 0) { - if (iarg+4 > narg) error->all(FLERR,"Illegal fix deform command"); + } else if (strcmp(arg[iarg + 1], "wiggle") == 0) { + if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, thiscmd + " wiggle", error); set[index].style = WIGGLE; - set[index].amplitude = utils::numeric(FLERR,arg[iarg+2],false,lmp); - set[index].tperiod = utils::numeric(FLERR,arg[iarg+3],false,lmp); + set[index].amplitude = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + set[index].tperiod = utils::numeric(FLERR, arg[iarg + 3], false, lmp); if (set[index].tperiod <= 0.0) - error->all(FLERR,"Illegal fix deform command"); + error->all(FLERR, "Illegal fix {} wiggle period, must be positive", style); iarg += 4; - } else if (strcmp(arg[iarg+1],"variable") == 0) { - if (iarg+4 > narg) error->all(FLERR,"Illegal fix deform command"); + } else if (strcmp(arg[iarg + 1], "variable") == 0) { + if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, thiscmd + " variable", error); set[index].style = VARIABLE; - if (strstr(arg[iarg+2],"v_") != arg[iarg+2]) - error->all(FLERR,"Illegal fix deform command"); - if (strstr(arg[iarg+3],"v_") != arg[iarg+3]) - error->all(FLERR,"Illegal fix deform command"); + if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) + error->all(FLERR, "Illegal fix {} variable name {}", style, arg[iarg + 2]); + if (strstr(arg[iarg + 3], "v_") != arg[iarg + 3]) + error->all(FLERR, "Illegal fix {} variable name {}", style, arg[iarg + 3]); delete[] set[index].hstr; delete[] set[index].hratestr; - set[index].hstr = utils::strdup(&arg[iarg+2][2]); - set[index].hratestr = utils::strdup(&arg[iarg+3][2]); + set[index].hstr = utils::strdup(&arg[iarg + 2][2]); + set[index].hratestr = utils::strdup(&arg[iarg + 3][2]); iarg += 4; - } else error->all(FLERR,"Illegal fix deform command"); - - } else if (strcmp(arg[iarg],"xy") == 0 || - strcmp(arg[iarg],"xz") == 0 || - strcmp(arg[iarg],"yz") == 0) { - - if (triclinic == 0) - error->all(FLERR,"Fix deform tilt factors require triclinic box"); - 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) error->all(FLERR,"Illegal fix deform command"); - if (strcmp(arg[iarg+1],"final") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix deform command"); + } else if (child_styles.find(arg[iarg + 1]) != child_styles.end()) { + nskip = child_styles[arg[iarg + 1]]; + if (iarg + nskip > narg) + utils::missing_cmd_args(FLERR, fmt::format("fix {} {}", style, arg[iarg + 1]), error); + for (int i = 0; i < nskip; i++) leftover_iarg.push_back(iarg + i); + iarg += nskip; + } else error->all(FLERR, "Illegal fix {} command argument: {}", style, arg[iarg + 1]); + + } else if ((strcmp(arg[iarg], "xy") == 0) + || (strcmp(arg[iarg], "xz") == 0) + || (strcmp(arg[iarg], "yz") == 0)) { + + if (triclinic == 0) error->all(FLERR,"Fix {} tilt factors require triclinic box", style); + 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, thiscmd, error); + if (strcmp(arg[iarg + 1], "final") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, thiscmd + " final", error); set[index].style = FINAL; - set[index].ftilt = utils::numeric(FLERR,arg[iarg+2],false,lmp); + set[index].ftilt = utils::numeric(FLERR, arg[iarg + 2], false, lmp); iarg += 3; - } else if (strcmp(arg[iarg+1],"delta") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix deform command"); + } else if (strcmp(arg[iarg + 1], "delta") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, thiscmd + " delta", error); set[index].style = DELTA; - set[index].dtilt = utils::numeric(FLERR,arg[iarg+2],false,lmp); + set[index].dtilt = utils::numeric(FLERR, arg[iarg + 2], false, lmp); iarg += 3; - } else if (strcmp(arg[iarg+1],"vel") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix deform command"); + } else if (strcmp(arg[iarg + 1], "vel") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, thiscmd + " vel", error); set[index].style = VEL; - set[index].vel = utils::numeric(FLERR,arg[iarg+2],false,lmp); + set[index].vel = utils::numeric(FLERR, arg[iarg + 2], false, lmp); iarg += 3; - } else if (strcmp(arg[iarg+1],"erate") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix deform command"); + } else if (strcmp(arg[iarg + 1], "erate") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, thiscmd + " erate", error); set[index].style = ERATE; - set[index].rate = utils::numeric(FLERR,arg[iarg+2],false,lmp); + set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp); iarg += 3; - } else if (strcmp(arg[iarg+1],"trate") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix deform command"); + } else if (strcmp(arg[iarg + 1], "trate") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, thiscmd + " trate", error); set[index].style = TRATE; - set[index].rate = utils::numeric(FLERR,arg[iarg+2],false,lmp); + set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp); iarg += 3; - } else if (strcmp(arg[iarg+1],"wiggle") == 0) { - if (iarg+4 > narg) error->all(FLERR,"Illegal fix deform command"); + } else if (strcmp(arg[iarg + 1], "wiggle") == 0) { + if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, thiscmd + " wiggle", error); set[index].style = WIGGLE; - set[index].amplitude = utils::numeric(FLERR,arg[iarg+2],false,lmp); - set[index].tperiod = utils::numeric(FLERR,arg[iarg+3],false,lmp); + set[index].amplitude = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + set[index].tperiod = utils::numeric(FLERR, arg[iarg + 3], false, lmp); if (set[index].tperiod <= 0.0) - error->all(FLERR,"Illegal fix deform command"); + error->all(FLERR, "Illegal fix {} wiggle period, must be positive", style); iarg += 4; - } else if (strcmp(arg[iarg+1],"variable") == 0) { - if (iarg+4 > narg) error->all(FLERR,"Illegal fix deform command"); + } else if (strcmp(arg[iarg + 1], "variable") == 0) { + if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, thiscmd + " variable", error); set[index].style = VARIABLE; - if (strstr(arg[iarg+2],"v_") != arg[iarg+2]) - error->all(FLERR,"Illegal fix deform command"); - if (strstr(arg[iarg+3],"v_") != arg[iarg+3]) - error->all(FLERR,"Illegal fix deform command"); + if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) + error->all(FLERR, "Illegal fix {} variable name {}", style, arg[iarg + 2]); + if (strstr(arg[iarg + 3], "v_") != arg[iarg + 3]) + error->all(FLERR, "Illegal fix {} variable name {}", style, arg[iarg + 3]); delete[] set[index].hstr; delete[] set[index].hratestr; - set[index].hstr = utils::strdup(&arg[iarg+2][2]); - set[index].hratestr = utils::strdup(&arg[iarg+3][2]); + set[index].hstr = utils::strdup(&arg[iarg + 2][2]); + set[index].hratestr = utils::strdup(&arg[iarg + 3][2]); iarg += 4; - } else error->all(FLERR,"Illegal fix deform command"); - + } else if (child_styles.find(arg[iarg + 1]) != child_styles.end()) { + nskip = child_styles[arg[iarg + 1]]; + if (iarg + nskip > narg) + utils::missing_cmd_args(FLERR, fmt::format("fix {} {}", style, arg[iarg + 1]), error); + for (int i = 0; i < nskip; i++) leftover_iarg.push_back(iarg + i); + iarg += nskip; + } else error->all(FLERR, "Illegal fix {} command argument: {}", style, arg[iarg + 1]); + } else if (child_parameters.find(arg[iarg]) != child_parameters.end()) { + if (child_styles.find(arg[iarg + 1]) != child_styles.end()) { + nskip = child_styles[arg[iarg + 1]]; + if (iarg + nskip > narg) + utils::missing_cmd_args(FLERR, fmt::format("fix {} {}", style, arg[iarg + 1]), error); + for (int i = 0; i < nskip; i++) leftover_iarg.push_back(iarg + i); + iarg += nskip; + } else error->all(FLERR, "Illegal fix {} command argument: {}", style, arg[iarg + 1]); } else break; } // read options from end of input line + + iarg_options_start = iarg; + options(narg - iarg, &arg[iarg]); + // no x remap effectively moves atoms within box, so set restart_pbc - options(narg-iarg,&arg[iarg]); if (remapflag != Domain::X_REMAP) restart_pbc = 1; // setup dimflags used by other classes to check for volume-change conflicts @@ -217,28 +251,19 @@ irregular(nullptr), set(nullptr) // no tensile deformation on shrink-wrapped dims // b/c shrink wrap will change box-length - if (set[0].style && - (domain->boundary[0][0] >= 2 || domain->boundary[0][1] >= 2)) - error->all(FLERR,"Cannot use fix deform on a shrink-wrapped boundary"); - if (set[1].style && - (domain->boundary[1][0] >= 2 || domain->boundary[1][1] >= 2)) - error->all(FLERR,"Cannot use fix deform on a shrink-wrapped boundary"); - if (set[2].style && - (domain->boundary[2][0] >= 2 || domain->boundary[2][1] >= 2)) - error->all(FLERR,"Cannot use fix deform on a shrink-wrapped boundary"); + for (int i = 0; i < 3; i++) + if (set[i].style && (domain->boundary[i][0] >= 2 || domain->boundary[i][1] >= 2)) + error->all(FLERR, "Cannot use fix {} on a shrink-wrapped boundary", style); // 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 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 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 tilt on a shrink-wrapped 2nd dim"); + if (set[3].style && (domain->boundary[2][0] >= 2 || domain->boundary[2][1] >= 2)) + error->all(FLERR, "Cannot use fix {} tilt on a shrink-wrapped 2nd dim", style); + if (set[4].style && (domain->boundary[2][0] >= 2 || domain->boundary[2][1] >= 2)) + error->all(FLERR, "Cannot use fix {} tilt on a shrink-wrapped 2nd dim", style); + if (set[5].style && (domain->boundary[1][0] >= 2 || domain->boundary[1][1] >= 2)) + error->all(FLERR, "Cannot use fix {} tilt on a shrink-wrapped 2nd dim", style); // apply scaling to FINAL,DELTA,VEL,WIGGLE since they have dist/vel units @@ -247,7 +272,7 @@ irregular(nullptr), set(nullptr) if (set[i].style == FINAL || set[i].style == DELTA || set[i].style == VEL || set[i].style == WIGGLE) flag = 1; - double xscale,yscale,zscale; + double xscale, yscale, zscale; if (flag && scaleflag) { xscale = domain->lattice->xlattice; yscale = domain->lattice->ylattice; @@ -284,40 +309,40 @@ irregular(nullptr), set(nullptr) // for VOLUME, setup links to other dims // fixed, dynamic1, dynamic2 + // only check for parent, otherwise child will check - for (int i = 0; i < 3; i++) { - if (set[i].style != VOLUME) continue; - int other1 = (i+1) % 3; - int other2 = (i+2) % 3; - - if (set[other1].style == NONE) { - if (set[other2].style == NONE || set[other2].style == VOLUME) - error->all(FLERR,"Fix deform volume setting is invalid"); - set[i].substyle = ONE_FROM_ONE; - set[i].fixed = other1; - set[i].dynamic1 = other2; - } else if (set[other2].style == NONE) { - if (set[other1].style == NONE || set[other1].style == VOLUME) - error->all(FLERR,"Fix deform volume setting is invalid"); - set[i].substyle = ONE_FROM_ONE; - set[i].fixed = other2; - set[i].dynamic1 = other1; - } else if (set[other1].style == VOLUME) { - if (set[other2].style == NONE || set[other2].style == VOLUME) - error->all(FLERR,"Fix deform volume setting is invalid"); - set[i].substyle = TWO_FROM_ONE; - set[i].fixed = other1; - set[i].dynamic1 = other2; - } else if (set[other2].style == VOLUME) { + if (strcmp(style, "deform") == 0) { + 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) - error->all(FLERR,"Fix deform volume setting is invalid"); - 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; + 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; + } } } @@ -348,8 +373,6 @@ irregular(nullptr), set(nullptr) if (force_reneighbor) irregular = new Irregular(lmp); else irregular = nullptr; - - TWOPI = 2.0*MY_PI; } /* ---------------------------------------------------------------------- */ @@ -394,7 +417,7 @@ void FixDeform::init() // domain, fix nvt/sllod, compute temp/deform only work on single h_rate if (modify->get_fix_by_style("deform").size() > 1) - error->all(FLERR,"More than one fix deform"); + error->all(FLERR, "More than one fix deform"); // Kspace setting @@ -411,14 +434,14 @@ void FixDeform::init() if (set[i].style != VARIABLE) continue; set[i].hvar = input->variable->find(set[i].hstr); if (set[i].hvar < 0) - error->all(FLERR,"Variable name for fix deform does not exist"); + error->all(FLERR, "Variable name {} for fix {} does not exist", set[i].hstr, style); if (!input->variable->equalstyle(set[i].hvar)) - error->all(FLERR,"Variable for fix deform is invalid style"); + error->all(FLERR, "Variable {} for fix {} is invalid style", set[i].hstr, style); set[i].hratevar = input->variable->find(set[i].hratestr); if (set[i].hratevar < 0) - error->all(FLERR,"Variable name for fix deform does not exist"); + error->all(FLERR, "Variable name {} for fix {} does not exist", set[i].hratestr, style); if (!input->variable->equalstyle(set[i].hratevar)) - error->all(FLERR,"Variable for fix deform is invalid style"); + error->all(FLERR, "Variable {} for fix {} is invalid style", set[i].hratestr, style); } // set start/stop values for box size and shape @@ -445,30 +468,26 @@ void FixDeform::init() set[i].lo_stop = set[i].lo_start + set[i].dlo; set[i].hi_stop = set[i].hi_start + set[i].dhi; } else if (set[i].style == SCALE) { - set[i].lo_stop = 0.5*(set[i].lo_start+set[i].hi_start) - - 0.5*set[i].scale*(set[i].hi_start-set[i].lo_start); - set[i].hi_stop = 0.5*(set[i].lo_start+set[i].hi_start) + - 0.5*set[i].scale*(set[i].hi_start-set[i].lo_start); + double shift = 0.5 * set[i].scale * (set[i].hi_start - set[i].lo_start); + set[i].lo_stop = 0.5 * (set[i].lo_start + set[i].hi_start) - shift; + set[i].hi_stop = 0.5 * (set[i].lo_start + set[i].hi_start) + shift; } else if (set[i].style == VEL) { - set[i].lo_stop = set[i].lo_start - 0.5*delt*set[i].vel; - set[i].hi_stop = set[i].hi_start + 0.5*delt*set[i].vel; + set[i].lo_stop = set[i].lo_start - 0.5 * delt * set[i].vel; + set[i].hi_stop = set[i].hi_start + 0.5 * delt * set[i].vel; } else if (set[i].style == ERATE) { - set[i].lo_stop = set[i].lo_start - - 0.5*delt*set[i].rate * (set[i].hi_start-set[i].lo_start); - set[i].hi_stop = set[i].hi_start + - 0.5*delt*set[i].rate * (set[i].hi_start-set[i].lo_start); + double shift = 0.5 * delt * set[i].rate * (set[i].hi_start - set[i].lo_start); + set[i].lo_stop = set[i].lo_start - shift; + set[i].hi_stop = set[i].hi_start + shift; if (set[i].hi_stop <= set[i].lo_stop) - error->all(FLERR,"Final box dimension due to fix deform is < 0.0"); + error->all(FLERR, "Final box dimension due to fix {} is < 0.0", style); } else if (set[i].style == TRATE) { - set[i].lo_stop = 0.5*(set[i].lo_start+set[i].hi_start) - - 0.5*((set[i].hi_start-set[i].lo_start) * exp(set[i].rate*delt)); - set[i].hi_stop = 0.5*(set[i].lo_start+set[i].hi_start) + - 0.5*((set[i].hi_start-set[i].lo_start) * exp(set[i].rate*delt)); + double shift = 0.5 * ((set[i].hi_start - set[i].lo_start) * exp(set[i].rate * delt)); + set[i].lo_stop = 0.5 * (set[i].lo_start + set[i].hi_start) - shift; + set[i].hi_stop = 0.5 * (set[i].lo_start + set[i].hi_start) + shift; } else if (set[i].style == WIGGLE) { - set[i].lo_stop = set[i].lo_start - - 0.5*set[i].amplitude * sin(TWOPI*delt/set[i].tperiod); - set[i].hi_stop = set[i].hi_start + - 0.5*set[i].amplitude * sin(TWOPI*delt/set[i].tperiod); + double shift = 0.5 * set[i].amplitude * sin(MY_2PI * delt / set[i].tperiod); + set[i].lo_stop = set[i].lo_start - shift; + set[i].hi_stop = set[i].hi_start + shift; } } @@ -484,50 +503,46 @@ void FixDeform::init() } else if (set[i].style == DELTA) { set[i].tilt_stop = set[i].tilt_start + set[i].dtilt; } else if (set[i].style == VEL) { - set[i].tilt_stop = set[i].tilt_start + delt*set[i].vel; + set[i].tilt_stop = set[i].tilt_start + delt * set[i].vel; } else if (set[i].style == ERATE) { if (i == 3) set[i].tilt_stop = set[i].tilt_start + - delt*set[i].rate * (set[2].hi_start-set[2].lo_start); + delt * set[i].rate * (set[2].hi_start - set[2].lo_start); if (i == 4) set[i].tilt_stop = set[i].tilt_start + - delt*set[i].rate * (set[2].hi_start-set[2].lo_start); + delt * set[i].rate * (set[2].hi_start - set[2].lo_start); if (i == 5) set[i].tilt_stop = set[i].tilt_start + - delt*set[i].rate * (set[1].hi_start-set[1].lo_start); + delt * set[i].rate * (set[1].hi_start - set[1].lo_start); } else if (set[i].style == TRATE) { - set[i].tilt_stop = set[i].tilt_start * exp(set[i].rate*delt); + set[i].tilt_stop = set[i].tilt_start * exp(set[i].rate * delt); } else if (set[i].style == WIGGLE) { - set[i].tilt_stop = set[i].tilt_start + - set[i].amplitude * sin(TWOPI*delt/set[i].tperiod); + double shift = set[i].amplitude * sin(MY_2PI * delt / set[i].tperiod); + set[i].tilt_stop = set[i].tilt_start + shift; // compute min/max for WIGGLE = extrema tilt factor will ever reach if (set[i].amplitude >= 0.0) { - if (delt < 0.25*set[i].tperiod) { + if (delt < 0.25 * set[i].tperiod) { set[i].tilt_min = set[i].tilt_start; - set[i].tilt_max = set[i].tilt_start + - set[i].amplitude*sin(TWOPI*delt/set[i].tperiod); - } else if (delt < 0.5*set[i].tperiod) { + set[i].tilt_max = set[i].tilt_start + shift; + } else if (delt < 0.5 * set[i].tperiod) { set[i].tilt_min = set[i].tilt_start; set[i].tilt_max = set[i].tilt_start + set[i].amplitude; - } else if (delt < 0.75*set[i].tperiod) { - set[i].tilt_min = set[i].tilt_start - - set[i].amplitude*sin(TWOPI*delt/set[i].tperiod); + } else if (delt < 0.75 * set[i].tperiod) { + set[i].tilt_min = set[i].tilt_start - shift; set[i].tilt_max = set[i].tilt_start + set[i].amplitude; } else { set[i].tilt_min = set[i].tilt_start - set[i].amplitude; set[i].tilt_max = set[i].tilt_start + set[i].amplitude; } } else { - if (delt < 0.25*set[i].tperiod) { - set[i].tilt_min = set[i].tilt_start - - set[i].amplitude*sin(TWOPI*delt/set[i].tperiod); + if (delt < 0.25 * set[i].tperiod) { + set[i].tilt_min = set[i].tilt_start - shift; set[i].tilt_max = set[i].tilt_start; - } else if (delt < 0.5*set[i].tperiod) { + } else if (delt < 0.5 * set[i].tperiod) { set[i].tilt_min = set[i].tilt_start - set[i].amplitude; set[i].tilt_max = set[i].tilt_start; - } else if (delt < 0.75*set[i].tperiod) { + } else if (delt < 0.75 * set[i].tperiod) { set[i].tilt_min = set[i].tilt_start - set[i].amplitude; - set[i].tilt_max = set[i].tilt_start + - set[i].amplitude*sin(TWOPI*delt/set[i].tperiod); + set[i].tilt_max = set[i].tilt_start + shift; } else { set[i].tilt_min = set[i].tilt_start - set[i].amplitude; set[i].tilt_max = set[i].tilt_start + set[i].amplitude; @@ -540,7 +555,7 @@ void FixDeform::init() for (int i = 3; i < 6; i++) if (set[i].style == TRATE && set[i].tilt_start == 0.0) - error->all(FLERR,"Cannot use fix deform trate on a box with zero tilt"); + error->all(FLERR, "Cannot use fix {} trate on a box with zero tilt", style); // if yz changes and will cause box flip, then xy cannot be changing // yz = [3], xy = [5] @@ -555,20 +570,20 @@ void FixDeform::init() int flag = 0; double lo,hi; if (flipflag && set[3].style == VARIABLE) - error->all(FLERR,"Fix deform cannot use yz variable with xy"); + error->all(FLERR, "Fix {} cannot use yz variable with xy", style); if (set[3].style == WIGGLE) { lo = set[3].tilt_min; hi = set[3].tilt_max; } else lo = hi = set[3].tilt_stop; if (flipflag) { - if (lo/(set[1].hi_start-set[1].lo_start) < -0.5 || - hi/(set[1].hi_start-set[1].lo_start) > 0.5) flag = 1; + if (lo / (set[1].hi_start - set[1].lo_start) < -0.5 || + hi / (set[1].hi_start - set[1].lo_start) > 0.5) flag = 1; if (set[1].style) { - if (lo/(set[1].hi_stop-set[1].lo_stop) < -0.5 || - hi/(set[1].hi_stop-set[1].lo_stop) > 0.5) flag = 1; + if (lo / (set[1].hi_stop - set[1].lo_stop) < -0.5 || + hi / (set[1].hi_stop - set[1].lo_stop) > 0.5) flag = 1; } if (flag) - error->all(FLERR,"Fix deform is changing yz too much with xy"); + error->all(FLERR, "Fix {} is changing yz too much with xy", style); } } @@ -584,7 +599,7 @@ void FixDeform::init() if (set[i].style == FINAL || set[i].style == DELTA || set[i].style == SCALE || set[i].style == VEL || set[i].style == ERATE) { - double dlo_dt,dhi_dt; + double dlo_dt, dhi_dt; if (delt != 0.0) { dlo_dt = (set[i].lo_stop - set[i].lo_start) / delt; dhi_dt = (set[i].hi_stop - set[i].hi_start) / delt; @@ -633,7 +648,7 @@ void FixDeform::pre_exchange() domain->set_global_box(); domain->set_local_box(); - domain->image_flip(flipxy,flipxz,flipyz); + domain->image_flip(flipxy, flipxz, flipyz); double **x = atom->x; imageint *image = atom->image; @@ -651,104 +666,72 @@ void FixDeform::pre_exchange() void FixDeform::end_of_step() { - int i; - - double delta = update->ntimestep - update->beginstep; - if (delta != 0.0) delta /= update->endstep - update->beginstep; - // wrap variable evaluations with clear/add if (varflag) modify->clearstep_compute(); - // set new box size + // set new box size for strain-based dims + + apply_strain(); + + // set new box size for VOLUME dims that are linked to other dims + // NOTE: still need to set h_rate for these dims + + apply_volume(); + + if (varflag) modify->addstep_compute(update->ntimestep + nevery); + + update_domain(); + + // redo KSpace coeffs since box has changed + + if (kspace_flag) force->kspace->setup(); +} + +/* ---------------------------------------------------------------------- + apply strain controls +------------------------------------------------------------------------- */ + +void FixDeform::apply_strain() +{ // for NONE, target is current box size // for TRATE, set target directly based on current time, also set h_rate // for WIGGLE, set target directly based on current time, also set h_rate // for VARIABLE, set target directly via variable eval, also set h_rate // for others except VOLUME, target is linear value between start and stop - for (i = 0; i < 3; i++) { + double delta = update->ntimestep - update->beginstep; + if (delta != 0.0) delta /= update->endstep - update->beginstep; + + for (int i = 0; i < 3; i++) { if (set[i].style == NONE) { set[i].lo_target = domain->boxlo[i]; set[i].hi_target = domain->boxhi[i]; } else if (set[i].style == TRATE) { double delt = (update->ntimestep - update->beginstep) * update->dt; - set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) - - 0.5*((set[i].hi_start-set[i].lo_start) * exp(set[i].rate*delt)); - set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) + - 0.5*((set[i].hi_start-set[i].lo_start) * exp(set[i].rate*delt)); + double shift = 0.5 * ((set[i].hi_start - set[i].lo_start) * exp(set[i].rate * delt)); + 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; h_rate[i] = set[i].rate * domain->h[i]; - h_ratelo[i] = -0.5*h_rate[i]; + h_ratelo[i] = -0.5 * h_rate[i]; } 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); - 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); - h_ratelo[i] = -0.5*h_rate[i]; + double shift = 0.5 * set[i].amplitude * sin(MY_2PI * delt / set[i].tperiod); + set[i].lo_target = set[i].lo_start - shift; + set[i].hi_target = set[i].hi_start + shift; + 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); - set[i].lo_target = set[i].lo_start - 0.5*del; - set[i].hi_target = set[i].hi_start + 0.5*del; + set[i].lo_target = set[i].lo_start - 0.5 * del; + set[i].hi_target = set[i].hi_start + 0.5 * del; h_rate[i] = input->variable->compute_equal(set[i].hratevar); - h_ratelo[i] = -0.5*h_rate[i]; - } else if (set[i].style != VOLUME) { - set[i].lo_target = set[i].lo_start + - delta*(set[i].lo_stop - set[i].lo_start); - set[i].hi_target = set[i].hi_start + - delta*(set[i].hi_stop - set[i].hi_start); - } - } - - // set new box size for VOLUME dims that are linked to other dims - // NOTE: still need to set h_rate for these dims - - for (i = 0; i < 3; i++) { - if (set[i].style != VOLUME) continue; - - if (set[i].substyle == ONE_FROM_ONE) { - set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) - - 0.5*(set[i].vol_start / - (set[set[i].dynamic1].hi_target - - set[set[i].dynamic1].lo_target) / - (set[set[i].fixed].hi_start-set[set[i].fixed].lo_start)); - set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) + - 0.5*(set[i].vol_start / - (set[set[i].dynamic1].hi_target - - set[set[i].dynamic1].lo_target) / - (set[set[i].fixed].hi_start-set[set[i].fixed].lo_start)); - - } else if (set[i].substyle == ONE_FROM_TWO) { - set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) - - 0.5*(set[i].vol_start / - (set[set[i].dynamic1].hi_target - - set[set[i].dynamic1].lo_target) / - (set[set[i].dynamic2].hi_target - - set[set[i].dynamic2].lo_target)); - set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) + - 0.5*(set[i].vol_start / - (set[set[i].dynamic1].hi_target - - set[set[i].dynamic1].lo_target) / - (set[set[i].dynamic2].hi_target - - set[set[i].dynamic2].lo_target)); - - } else if (set[i].substyle == TWO_FROM_ONE) { - set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) - - 0.5*sqrt(set[i].vol_start / - (set[set[i].dynamic1].hi_target - - set[set[i].dynamic1].lo_target) / - (set[set[i].fixed].hi_start - - set[set[i].fixed].lo_start) * - (set[i].hi_start - set[i].lo_start)); - set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) + - 0.5*sqrt(set[i].vol_start / - (set[set[i].dynamic1].hi_target - - set[set[i].dynamic1].lo_target) / - (set[set[i].fixed].hi_start - - set[set[i].fixed].lo_start) * - (set[i].hi_start - set[i].lo_start)); + h_ratelo[i] = -0.5 * h_rate[i]; + } else if (set[i].style == FINAL || set[i].style == DELTA || set[i].style == SCALE || + set[i].style == VEL || set[i].style == ERATE) { + set[i].lo_target = set[i].lo_start + delta * (set[i].lo_stop - set[i].lo_start); + set[i].hi_target = set[i].hi_start + delta * (set[i].hi_stop - set[i].hi_start); } } @@ -760,55 +743,97 @@ void FixDeform::end_of_step() // for other styles, target is linear value between start and stop values if (triclinic) { - double *h = domain->h; - - for (i = 3; i < 6; i++) { + for (int i = 3; i < 6; i++) { if (set[i].style == NONE) { if (i == 5) set[i].tilt_target = domain->xy; else if (i == 4) set[i].tilt_target = domain->xz; else if (i == 3) set[i].tilt_target = domain->yz; } else if (set[i].style == TRATE) { double delt = (update->ntimestep - update->beginstep) * update->dt; - set[i].tilt_target = set[i].tilt_start * exp(set[i].rate*delt); + set[i].tilt_target = set[i].tilt_start * exp(set[i].rate * delt); h_rate[i] = set[i].rate * domain->h[i]; } 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; h_rate[i] = input->variable->compute_equal(set[i].hratevar); } else { - set[i].tilt_target = set[i].tilt_start + - delta*(set[i].tilt_stop - set[i].tilt_start); + set[i].tilt_target = set[i].tilt_start + delta * (set[i].tilt_stop - set[i].tilt_start); } + } + } +} + +/* ---------------------------------------------------------------------- + apply volume controls +------------------------------------------------------------------------- */ + +void FixDeform::apply_volume() +{ + 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) { + 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)); + } + + 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; + } +} + +/* ---------------------------------------------------------------------- + Update box domain +------------------------------------------------------------------------- */ - // tilt_target can be large positive or large negative value - // add/subtract box lengths until tilt_target is closest to current value +void FixDeform::update_domain() +{ + // tilt_target can be large positive or large negative value + // add/subtract box lengths until tilt_target is closest to current value + if (triclinic) { + double *h = domain->h; + for (int i = 3; i < 6; i++) { int idenom = 0; if (i == 5) idenom = 0; else if (i == 4) idenom = 0; else if (i == 3) idenom = 1; double denom = set[idenom].hi_target - set[idenom].lo_target; - double current = h[i]/h[idenom]; + double current = h[i] / h[idenom]; - while (set[i].tilt_target/denom - current > 0.0) + while (set[i].tilt_target / denom - current > 0.0) set[i].tilt_target -= denom; - while (set[i].tilt_target/denom - current < 0.0) + while (set[i].tilt_target / denom - current < 0.0) set[i].tilt_target += denom; - if (fabs(set[i].tilt_target/denom - 1.0 - current) < - fabs(set[i].tilt_target/denom - current)) + if (fabs(set[i].tilt_target / denom - 1.0 - current) < + fabs(set[i].tilt_target / denom - current)) set[i].tilt_target -= denom; } } - if (varflag) modify->addstep_compute(update->ntimestep + nevery); - // if any tilt ratios exceed 0.5, set flip = 1 and compute new tilt values // do not flip in x or y if non-periodic (can tilt but not flip) // this is b/c the box length would be changed (dramatically) by flip @@ -823,12 +848,12 @@ void FixDeform::end_of_step() double yprd = set[1].hi_target - set[1].lo_target; double xprdinv = 1.0 / xprd; double yprdinv = 1.0 / yprd; - if (set[3].tilt_target*yprdinv < -0.5 || - set[3].tilt_target*yprdinv > 0.5 || - set[4].tilt_target*xprdinv < -0.5 || - set[4].tilt_target*xprdinv > 0.5 || - set[5].tilt_target*xprdinv < -0.5 || - set[5].tilt_target*xprdinv > 0.5) { + if (set[3].tilt_target * yprdinv < -0.5 || + set[3].tilt_target * yprdinv > 0.5 || + set[4].tilt_target * xprdinv < -0.5 || + set[4].tilt_target * xprdinv > 0.5 || + set[5].tilt_target * xprdinv < -0.5 || + set[5].tilt_target * xprdinv > 0.5) { set[3].tilt_flip = set[3].tilt_target; set[4].tilt_flip = set[4].tilt_target; set[5].tilt_flip = set[5].tilt_target; @@ -836,30 +861,30 @@ void FixDeform::end_of_step() flipxy = flipxz = flipyz = 0; if (domain->yperiodic) { - if (set[3].tilt_flip*yprdinv < -0.5) { + if (set[3].tilt_flip * yprdinv < -0.5) { set[3].tilt_flip += yprd; set[4].tilt_flip += set[5].tilt_flip; flipyz = 1; - } else if (set[3].tilt_flip*yprdinv > 0.5) { + } else if (set[3].tilt_flip * yprdinv > 0.5) { set[3].tilt_flip -= yprd; set[4].tilt_flip -= set[5].tilt_flip; flipyz = -1; } } if (domain->xperiodic) { - if (set[4].tilt_flip*xprdinv < -0.5) { + if (set[4].tilt_flip * xprdinv < -0.5) { set[4].tilt_flip += xprd; flipxz = 1; } - if (set[4].tilt_flip*xprdinv > 0.5) { + if (set[4].tilt_flip * xprdinv > 0.5) { set[4].tilt_flip -= xprd; flipxz = -1; } - if (set[5].tilt_flip*xprdinv < -0.5) { + if (set[5].tilt_flip * xprdinv < -0.5) { set[5].tilt_flip += xprd; flipxy = 1; } - if (set[5].tilt_flip*xprdinv > 0.5) { + if (set[5].tilt_flip * xprdinv > 0.5) { set[5].tilt_flip -= xprd; flipxy = -1; } @@ -878,9 +903,9 @@ void FixDeform::end_of_step() int *mask = atom->mask; int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) + for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) - domain->x2lamda(x[i],x[i]); + domain->x2lamda(x[i], x[i]); for (auto &ifix : rfix) ifix->deform(0); @@ -889,22 +914,22 @@ void FixDeform::end_of_step() // reset global and local box to new size/shape // only if deform fix is controlling the dimension - if (set[0].style) { + if (dimflag[0]) { domain->boxlo[0] = set[0].lo_target; domain->boxhi[0] = set[0].hi_target; } - if (set[1].style) { + if (dimflag[1]) { domain->boxlo[1] = set[1].lo_target; domain->boxhi[1] = set[1].hi_target; } - if (set[2].style) { + if (dimflag[2]) { domain->boxlo[2] = set[2].lo_target; domain->boxhi[2] = set[2].hi_target; } if (triclinic) { - if (set[3].style) domain->yz = set[3].tilt_target; - if (set[4].style) domain->xz = set[4].tilt_target; - if (set[5].style) domain->xy = set[5].tilt_target; + if (dimflag[3]) domain->yz = set[3].tilt_target; + if (dimflag[4]) domain->xz = set[4].tilt_target; + if (dimflag[5]) domain->xy = set[5].tilt_target; } domain->set_global_box(); @@ -917,17 +942,13 @@ void FixDeform::end_of_step() int *mask = atom->mask; int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) + for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) - domain->lamda2x(x[i],x[i]); + domain->lamda2x(x[i], x[i]); for (auto &ifix : rfix) ifix->deform(1); } - - // redo KSpace coeffs since box has changed - - if (kspace_flag) force->kspace->setup(); } /* ---------------------------------------------------------------------- @@ -937,9 +958,9 @@ void FixDeform::end_of_step() void FixDeform::write_restart(FILE *fp) { if (comm->me == 0) { - int size = 6*sizeof(Set); - fwrite(&size,sizeof(int),1,fp); - fwrite(set,sizeof(Set),6,fp); + int size = 6 * sizeof(Set); + fwrite(&size, sizeof(int), 1, fp); + fwrite(set, sizeof(Set), 6, fp); } } @@ -951,7 +972,7 @@ void FixDeform::restart(char *buf) { int samestyle = 1; Set *set_restart = (Set *) buf; - for (int i=0; i<6; ++i) { + 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; @@ -964,39 +985,60 @@ void FixDeform::restart(char *buf) samestyle = 0; } if (!samestyle) - error->all(FLERR,"Fix deform settings not consistent with restart"); + error->all(FLERR, "Fix {} settings not consistent with restart", style); } /* ---------------------------------------------------------------------- */ void FixDeform::options(int narg, char **arg) { - if (narg < 0) error->all(FLERR,"Illegal fix deform command"); + const std::string thiscmd = fmt::format("fix {}", style); + if (narg < 0) utils::missing_cmd_args(FLERR, thiscmd, error); remapflag = Domain::X_REMAP; scaleflag = 1; flipflag = 1; + // arguments for child classes + + std::unordered_map child_options; + int nskip; + if (utils::strmatch(style, "^deform/pressure")) { + child_options.insert({{"couple", 2}, + {"max/rate", 2}, + {"normalize/pressure", 2}, + {"vol/balance/p", 2}}); + } + + // parse all optional arguments for this parent and also child classes + // for child classes, simply store them in leftover_iarg and skip over them + int iarg = 0; while (iarg < narg) { - if (strcmp(arg[iarg],"remap") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix deform command"); - if (strcmp(arg[iarg+1],"x") == 0) remapflag = Domain::X_REMAP; - else if (strcmp(arg[iarg+1],"v") == 0) remapflag = Domain::V_REMAP; - else if (strcmp(arg[iarg+1],"none") == 0) remapflag = Domain::NO_REMAP; - else error->all(FLERR,"Illegal fix deform command"); + if (strcmp(arg[iarg], "remap") == 0) { + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, thiscmd + " remap", error); + if (strcmp(arg[iarg + 1], "x") == 0) remapflag = Domain::X_REMAP; + else if (strcmp(arg[iarg + 1], "v") == 0) remapflag = Domain::V_REMAP; + else if (strcmp(arg[iarg + 1], "none") == 0) remapflag = Domain::NO_REMAP; + else error->all(FLERR, "Illegal fix {} remap command: {}", style, arg[iarg + 1]); iarg += 2; - } else if (strcmp(arg[iarg],"units") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix deform command"); - if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; - else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; - else error->all(FLERR,"Illegal fix deform command"); + } else if (strcmp(arg[iarg], "units") == 0) { + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, thiscmd + " units", error); + if (strcmp(arg[iarg + 1], "box") == 0) scaleflag = 0; + else if (strcmp(arg[iarg + 1], "lattice") == 0) scaleflag = 1; + else error->all(FLERR, "Illegal fix {} units command: {}", style, arg[iarg + 1]); iarg += 2; - } else if (strcmp(arg[iarg],"flip") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix deform command"); - flipflag = utils::logical(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "flip") == 0) { + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, thiscmd + " flip", error); + flipflag = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else error->all(FLERR,"Illegal fix deform command"); + } else if (child_options.find(arg[iarg]) != child_options.end()) { + nskip = child_options[arg[iarg]]; + if (iarg + nskip > narg) + utils::missing_cmd_args(FLERR, fmt::format("fix {} {}", style, arg[iarg]), error); + for (int i = 0; i < nskip; i++) leftover_iarg.push_back(iarg + i); + iarg += nskip; + } else error->all(FLERR, "Unknown fix {} keyword: {}", style, arg[iarg]); } } diff --git a/src/fix_deform.h b/src/fix_deform.h index 20f6ac59014..b1337294448 100644 --- a/src/fix_deform.h +++ b/src/fix_deform.h @@ -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 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 leftover_iarg; + int iarg_options_start; + void options(int, char **); + void virtual apply_volume(); + void apply_strain(); + void update_domain(); }; } // namespace LAMMPS_NS