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

Merge branch 'collected-small-changes' into doxygen-integration

Browse Source
This commit is contained in:
Axel Kohlmeyer
2020-08-26 11:51:45 -04:00
parent 1c10aa6a4d bdc952ac88
commit 8137ada848
13 changed files with 479 additions and 227 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
doc/src/Errors_messages.rst
View File

@ -502,7 +502,8 @@ Doc page with :doc:`WARNING messages <Errors_warnings>`
*Bond/react: Unknown section in map file* *Bond/react: Unknown section in map file*
Please ensure reaction map files are properly formatted. Please ensure reaction map files are properly formatted.
*Bond/react: Atom affected by reaction too close to template edge* *Bond/react: Atom type affected by reaction too close to template edge*
*Bond/react: Bond type affected by reaction too close to template edge*
This means an atom which changes type or connectivity during the This means an atom which changes type or connectivity during the
reaction is too close to an 'edge' atom defined in the map reaction is too close to an 'edge' atom defined in the map
file. This could cause incorrect assignment of bonds, angle, etc. file. This could cause incorrect assignment of bonds, angle, etc.

31
doc/src/fix_bond_react.rst
View File

@ -14,19 +14,22 @@ Syntax
react react-ID react-group-ID Nevery Rmin Rmax template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ... react react-ID react-group-ID Nevery Rmin Rmax template-ID(pre-reacted) template-ID(post-reacted) map_file individual_keyword values ...
... ...
* ID, group-ID are documented in :doc:`fix <fix>` command. Group-ID is ignored. * ID, group-ID are documented in :doc:`fix <fix>` command.
* bond/react = style name of this fix command * bond/react = style name of this fix command
* the common keyword/values may be appended directly after 'bond/react' * the common keyword/values may be appended directly after 'bond/react'
* this applies to all reaction specifications (below) * this applies to all reaction specifications (below)
* common_keyword = *stabilization* * common_keyword = *stabilization* or *reset_mol_ids*
.. parsed-literal:: .. parsed-literal::
*stabilization* values = *no* or *yes* *group-ID* *xmax* *stabilization* values = *no* or *yes* *group-ID* *xmax*
*no* = no reaction site stabilization *no* = no reaction site stabilization (default)
*yes* = perform reaction site stabilization *yes* = perform reaction site stabilization
*group-ID* = user-assigned prefix for the dynamic group of atoms not currently involved in a reaction *group-ID* = user-assigned prefix for the dynamic group of atoms not currently involved in a reaction
*xmax* = xmax value that is used by an internally-created :doc:`nve/limit <fix_nve_limit>` integrator *xmax* = xmax value that is used by an internally-created :doc:`nve/limit <fix_nve_limit>` integrator
*reset_mol_ids* values = *yes* or *no*
*yes* = update molecule IDs based on new global topology (default)
*no* = do not update molecule IDs
* react = mandatory argument indicating new reaction specification * react = mandatory argument indicating new reaction specification
* react-ID = user-assigned name for the reaction * react-ID = user-assigned name for the reaction
@ -50,9 +53,9 @@ Syntax
*stabilize_steps* value = timesteps *stabilize_steps* value = timesteps
timesteps = number of timesteps to apply the internally-created :doc:`nve/limit <fix_nve_limit>` fix to reacting atoms timesteps = number of timesteps to apply the internally-created :doc:`nve/limit <fix_nve_limit>` fix to reacting atoms
*update_edges* value = *none* or *charges* or *custom* *update_edges* value = *none* or *charges* or *custom*
none = do not update topology near the edges of reaction templates *none* = do not update topology near the edges of reaction templates
charges = update atomic charges of all atoms in reaction templates *charges* = update atomic charges of all atoms in reaction templates
custom = force the update of user-specified atomic charges *custom* = force the update of user-specified atomic charges
Examples Examples
"""""""" """"""""
@ -154,6 +157,13 @@ due to the internal dynamic grouping performed by fix bond/react.
If the group-ID is an existing static group, react-group-IDs If the group-ID is an existing static group, react-group-IDs
should also be specified as this static group, or a subset. should also be specified as this static group, or a subset.
The *reset_mol_ids* keyword invokes the :doc:`reset_mol_ids <reset_mol_ids>`
command after a reaction occurs, to ensure that molecule IDs are
consistent with the new bond topology. The group-ID used for
:doc:`reset_mol_ids <reset_mol_ids>` is the group-ID for this fix.
Resetting molecule IDs is necessarily a global operation, and so can
be slow for very large systems.
The following comments pertain to each *react* argument (in other The following comments pertain to each *react* argument (in other
words, can be customized for each reaction, or reaction step): words, can be customized for each reaction, or reaction step):
@ -203,9 +213,10 @@ surrounding topology. As described below, the bonding atom pairs of
the pre-reacted template are specified by atom ID in the map file. The the pre-reacted template are specified by atom ID in the map file. The
pre-reacted molecule template should contain as few atoms as possible pre-reacted molecule template should contain as few atoms as possible
while still completely describing the topology of all atoms affected while still completely describing the topology of all atoms affected
by the reaction. For example, if the force field contains dihedrals, by the reaction (which includes all atoms that change atom type or
the pre-reacted template should contain any atom within three bonds of connectivity, and all bonds that change bond type). For example, if
reacting atoms. the force field contains dihedrals, the pre-reacted template should
contain any atom within three bonds of reacting atoms.
Some atoms in the pre-reacted template that are not reacting may have Some atoms in the pre-reacted template that are not reacting may have
missing topology with respect to the simulation. For example, the missing topology with respect to the simulation. For example, the
@ -554,7 +565,7 @@ Default
""""""" """""""
The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60, The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60,
update_edges = none reset_mol_ids = yes, update_edges = none
---------- ----------

143
doc/src/pair_granular.rst
View File

@ -93,7 +93,7 @@ on particle *i* due to contact with particle *j* is given by:
.. math:: .. math::
\mathbf{F}_{ne, Hooke} = k_N \delta_{ij} \mathbf{n} \mathbf{F}_{ne, Hooke} = k_n \delta_{ij} \mathbf{n}
Where :math:`\delta_{ij} = R_i + R_j - \|\mathbf{r}_{ij}\|` is the particle Where :math:`\delta_{ij} = R_i + R_j - \|\mathbf{r}_{ij}\|` is the particle
overlap, :math:`R_i, R_j` are the particle radii, :math:`\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j` is the vector separating the two overlap, :math:`R_i, R_j` are the particle radii, :math:`\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j` is the vector separating the two
@ -106,7 +106,7 @@ For the *hertz* model, the normal component of force is given by:
.. math:: .. math::
\mathbf{F}_{ne, Hertz} = k_N R_{eff}^{1/2}\delta_{ij}^{3/2} \mathbf{n} \mathbf{F}_{ne, Hertz} = k_n R_{eff}^{1/2}\delta_{ij}^{3/2} \mathbf{n}
Here, :math:`R_{eff} = \frac{R_i R_j}{R_i + R_j}` is the effective Here, :math:`R_{eff} = \frac{R_i R_j}{R_i + R_j}` is the effective
radius, denoted for simplicity as *R* from here on. For *hertz*\ , the radius, denoted for simplicity as *R* from here on. For *hertz*\ , the
@ -123,7 +123,7 @@ Here, :math:`E_{eff} = E = \left(\frac{1-\nu_i^2}{E_i} + \frac{1-\nu_j^2}{E_j}\r
modulus, with :math:`\nu_i, \nu_j` the Poisson ratios of the particles of modulus, with :math:`\nu_i, \nu_j` the Poisson ratios of the particles of
types *i* and *j*\ . Note that if the elastic modulus and the shear types *i* and *j*\ . Note that if the elastic modulus and the shear
modulus of the two particles are the same, the *hertz/material* model modulus of the two particles are the same, the *hertz/material* model
is equivalent to the *hertz* model with :math:`k_N = 4/3 E_{eff}` is equivalent to the *hertz* model with :math:`k_n = 4/3 E_{eff}`
The *dmt* model corresponds to the The *dmt* model corresponds to the
:ref:`(Derjaguin-Muller-Toporov) <DMT1975>` cohesive model, where the force :ref:`(Derjaguin-Muller-Toporov) <DMT1975>` cohesive model, where the force
@ -140,7 +140,7 @@ where the force is computed as:
\mathbf{F}_{ne, jkr} = \left(\frac{4Ea^3}{3R} - 2\pi a^2\sqrt{\frac{4\gamma E}{\pi a}}\right)\mathbf{n} \mathbf{F}_{ne, jkr} = \left(\frac{4Ea^3}{3R} - 2\pi a^2\sqrt{\frac{4\gamma E}{\pi a}}\right)\mathbf{n}
Here, *a* is the radius of the contact zone, related to the overlap Here, :math:`a` is the radius of the contact zone, related to the overlap
:math:`\delta` according to: :math:`\delta` according to:
.. math:: .. math::
@ -167,7 +167,7 @@ following general form:
\mathbf{F}_{n,damp} = -\eta_n \mathbf{v}_{n,rel} \mathbf{F}_{n,damp} = -\eta_n \mathbf{v}_{n,rel}
Here, :math:`\mathbf{v}_{n,rel} = (\mathbf{v}_j - \mathbf{v}_i) \cdot \mathbf{n} \mathbf{n}` is the component of relative velocity along Here, :math:`\mathbf{v}_{n,rel} = (\mathbf{v}_j - \mathbf{v}_i) \cdot \mathbf{n}\ \mathbf{n}` is the component of relative velocity along
:math:`\mathbf{n}`. :math:`\mathbf{n}`.
The optional *damping* keyword to the *pair_coeff* command followed by The optional *damping* keyword to the *pair_coeff* command followed by
@ -259,7 +259,9 @@ tangential model choices and their expected parameters are as follows:
1. *linear_nohistory* : :math:`x_{\gamma,t}`, :math:`\mu_s` 1. *linear_nohistory* : :math:`x_{\gamma,t}`, :math:`\mu_s`
2. *linear_history* : :math:`k_t`, :math:`x_{\gamma,t}`, :math:`\mu_s` 2. *linear_history* : :math:`k_t`, :math:`x_{\gamma,t}`, :math:`\mu_s`
3. *mindlin* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s` 3. *mindlin* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s`
4. *mindlin_rescale* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s` 4. *mindlin/force* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s`
5. *mindlin_rescale* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s`
6. *mindlin_rescale/force* : :math:`k_t` or NULL, :math:`x_{\gamma,t}`, :math:`\mu_s`
Here, :math:`x_{\gamma,t}` is a dimensionless multiplier for the normal Here, :math:`x_{\gamma,t}` is a dimensionless multiplier for the normal
damping :math:`\eta_n` that determines the magnitude of the tangential damping :math:`\eta_n` that determines the magnitude of the tangential
@ -268,11 +270,11 @@ coefficient, and :math:`k_t` is the tangential stiffness coefficient.
For *tangential linear_nohistory*, a simple velocity-dependent Coulomb For *tangential linear_nohistory*, a simple velocity-dependent Coulomb
friction criterion is used, which mimics the behavior of the *pair friction criterion is used, which mimics the behavior of the *pair
gran/hooke* style. The tangential force (\mathbf{F}_t\) is given by: gran/hooke* style. The tangential force :math:`\mathbf{F}_t` is given by:
.. math:: .. math::
\mathbf{F}_t = -min(\mu_t F_{n0}, \|\mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} \mathbf{F}_t = -\min(\mu_t F_{n0}, \|\mathbf{F}_\mathrm{t,damp}\|) \mathbf{t}
The tangential damping force :math:`\mathbf{F}_\mathrm{t,damp}` is given by: The tangential damping force :math:`\mathbf{F}_\mathrm{t,damp}` is given by:
@ -294,8 +296,8 @@ keyword also affects the tangential damping. The parameter
literature use :math:`x_{\gamma,t} = 1` (:ref:`Marshall <Marshall2009>`, literature use :math:`x_{\gamma,t} = 1` (:ref:`Marshall <Marshall2009>`,
:ref:`Tsuji et al <Tsuji1992>`, :ref:`Silbert et al <Silbert2001>`). The relative :ref:`Tsuji et al <Tsuji1992>`, :ref:`Silbert et al <Silbert2001>`). The relative
tangential velocity at the point of contact is given by tangential velocity at the point of contact is given by
:math:`\mathbf{v}_{t, rel} = \mathbf{v}_{t} - (R_i\Omega_i + R_j\Omega_j) \times \mathbf{n}`, where :math:`\mathbf{v}_{t} = \mathbf{v}_r - \mathbf{v}_r\cdot\mathbf{n}{n}`, :math:`\mathbf{v}_{t, rel} = \mathbf{v}_{t} - (R_i\mathbf{\Omega}_i + R_j\mathbf{\Omega}_j) \times \mathbf{n}`, where :math:`\mathbf{v}_{t} = \mathbf{v}_r - \mathbf{v}_r\cdot\mathbf{n}\ \mathbf{n}`,
:math:`\mathbf{v}_r = \mathbf{v}_j - \mathbf{v}_i`. :math:`\mathbf{v}_r = \mathbf{v}_j - \mathbf{v}_i` .
The direction of the applied force is :math:`\mathbf{t} = \mathbf{v_{t,rel}}/\|\mathbf{v_{t,rel}}\|` . The direction of the applied force is :math:`\mathbf{t} = \mathbf{v_{t,rel}}/\|\mathbf{v_{t,rel}}\|` .
The normal force value :math:`F_{n0}` used to compute the critical force The normal force value :math:`F_{n0}` used to compute the critical force
@ -314,21 +316,24 @@ form:
.. math:: .. math::
F_{n0} = \|\mathbf{F}_ne + 2 F_{pulloff}\| F_{n0} = \|\mathbf{F}_{ne} + 2 F_{pulloff}\|
Where :math:`F_{pulloff} = 3\pi \gamma R` for *jkr*\ , and Where :math:`F_{pulloff} = 3\pi \gamma R` for *jkr*\ , and
:math:`F_{pulloff} = 4\pi \gamma R` for *dmt*\ . :math:`F_{pulloff} = 4\pi \gamma R` for *dmt*\ .
The remaining tangential options all use accumulated tangential The remaining tangential options all use accumulated tangential
displacement (i.e. contact history). This is discussed below in the displacement (i.e. contact history), except for the options
context of the *linear_history* option, but the same treatment of the *mindlin/force* and *mindlin_rescale/force*, that use accumulated
accumulated displacement applies to the other options as well. tangential force instead, and are discussed further below.
The accumulated tangential displacement is discussed in details below
in the context of the *linear_history* option. The same treatment of
the accumulated displacement applies to the other options as well.
For *tangential linear_history*, the tangential force is given by: For *tangential linear_history*, the tangential force is given by:
.. math:: .. math::
\mathbf{F}_t = -min(\mu_t F_{n0}, \|-k_t\mathbf{\xi} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} \mathbf{F}_t = -\min(\mu_t F_{n0}, \|-k_t\mathbf{\xi} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t}
Here, :math:`\mathbf{\xi}` is the tangential displacement accumulated Here, :math:`\mathbf{\xi}` is the tangential displacement accumulated
during the entire duration of the contact: during the entire duration of the contact:
@ -356,7 +361,7 @@ work:
.. math:: .. math::
\mathbf{\xi} = \left(\mathbf{\xi'} - (\mathbf{n} \cdot \mathbf{\xi'})\mathbf{n}\right) \frac{\|\mathbf{\xi'}\|}{\|\mathbf{\xi'}\| - \mathbf{n}\cdot\mathbf{\xi'}} \mathbf{\xi} = \left(\mathbf{\xi'} - (\mathbf{n} \cdot \mathbf{\xi'})\mathbf{n}\right) \frac{\|\mathbf{\xi'}\|}{\|\mathbf{\xi'} - (\mathbf{n}\cdot\mathbf{\xi'})\mathbf{n}\|}
Here, :math:`\mathbf{\xi'}` is the accumulated displacement prior to the Here, :math:`\mathbf{\xi'}` is the accumulated displacement prior to the
current time step and :math:`\mathbf{\xi}` is the corrected current time step and :math:`\mathbf{\xi}` is the corrected
@ -372,7 +377,7 @@ discussion):
.. math:: .. math::
\mathbf{\xi} = -\frac{1}{k_t}\left(\mu_t F_{n0}\mathbf{t} + \mathbf{F}_{t,damp}\right) \mathbf{\xi} = -\frac{1}{k_t}\left(\mu_t F_{n0}\mathbf{t} - \mathbf{F}_{t,damp}\right)
The tangential force is added to the total normal force (elastic plus The tangential force is added to the total normal force (elastic plus
damping) to produce the total force on the particle. The tangential damping) to produce the total force on the particle. The tangential
@ -387,27 +392,68 @@ overlap region) to induce a torque on each particle according to:
\mathbf{\tau}_j = -(R_j - 0.5 \delta) \mathbf{n} \times \mathbf{F}_t \mathbf{\tau}_j = -(R_j - 0.5 \delta) \mathbf{n} \times \mathbf{F}_t
For *tangential mindlin*\ , the :ref:`Mindlin <Mindlin1949>` no-slip solution is used, which differs from the *linear_history* For *tangential mindlin*\ , the :ref:`Mindlin <Mindlin1949>` no-slip solution
option by an additional factor of *a*\ , the radius of the contact region. The tangential force is given by: is used which differs from the *linear_history* option by an additional factor
of :math:`a`, the radius of the contact region. The tangential force is given by:
.. math:: .. math::
\mathbf{F}_t = -min(\mu_t F_{n0}, \|-k_t a \mathbf{\xi} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t} \mathbf{F}_t = -\min(\mu_t F_{n0}, \|-k_t a \mathbf{\xi} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t}
Here, *a* is the radius of the contact region, given by :math:`a =\sqrt{R\delta}`
Here, :math:`a` is the radius of the contact region, given by :math:`a =\sqrt{R\delta}`
for all normal contact models, except for *jkr*\ , where it is given for all normal contact models, except for *jkr*\ , where it is given
implicitly by :math:`\delta = a^2/R - 2\sqrt{\pi \gamma a/E}`, see implicitly by :math:`\delta = a^2/R - 2\sqrt{\pi \gamma a/E}`, see
discussion above. To match the Mindlin solution, one should set :math:`k_t = 4G/(2-\nu)`, where :math:`G` is the shear modulus, related to Young's modulus discussion above. To match the Mindlin solution, one should set
:math:`E` by :math:`G = E/(2(1+\nu))`, where :math:`\nu` is Poisson's ratio. This :math:`k_t = 8G_{eff}`, where :math:`G_{eff}` is the effective shear modulus given by:
can also be achieved by specifying *NULL* for :math:`k_t`, in which case a
.. math::
G_{eff} = \left(\frac{2-\nu_i}{G_i} + \frac{2-\nu_j}{G_j}\right)^{-1}
where :math:`G` is the shear modulus, related to Young's modulus :math:`E`
and Poisson's ratio :math:`\nu` by :math:`G = E/(2(1+\nu))`. This can also be
achieved by specifying *NULL* for :math:`k_t`, in which case a
normal contact model that specifies material parameters :math:`E` and normal contact model that specifies material parameters :math:`E` and
:math:`\nu` is required (e.g. *hertz/material*\ , *dmt* or *jkr*\ ). In this :math:`\nu` is required (e.g. *hertz/material*\ , *dmt* or *jkr*\ ). In this
case, mixing of the shear modulus for different particle types *i* and case, mixing of the shear modulus for different particle types *i* and
*j* is done according to: *j* is done according to the formula above.
.. note::
The radius of the contact region :math:`a` depends on the normal overlap.
As a result, the tangential force for *mindlin* can change due to
a variation in normal overlap, even with no change in tangential displacement.
For *tangential mindlin/force*, the accumulated elastic tangential force
characterizes the contact history, instead of the accumulated tangential
displacement. This prevents the dependence of the tangential force on the
normal overlap as noted above. The tangential force is given by:
.. math:: .. math::
1/G = 2(2-\nu_i)(1+\nu_i)/E_i + 2(2-\nu_j)(1+\nu_j)/E_j \mathbf{F}_t = -\min(\mu_t F_{n0}, \|\mathbf{F}_{te} + \mathbf{F}_\mathrm{t,damp}\|) \mathbf{t}
The increment of the elastic component of the tangential force
:math:`\mathbf{F}_{te}` is given by:
.. math::
\mathrm{d}\mathbf{F}_{te} = -k_t a \mathbf{v}_{t,rel} \mathrm{d}\tau
The changes in frame of reference of the contacting pair of particles during
contact are accounted for by the same formula as above, replacing the
accumulated tangential displacement :math:`\xi`, by the accumulated tangential
elastic force :math:`F_{te}`. When the tangential force exceeds the critical
force, the tangential force is directly re-scaled to match the value for
the critical force:
.. math::
\mathbf{F}_{te} = - \mu_t F_{n0}\mathbf{t} + \mathbf{F}_{t,damp}
The same rules as those described for *mindlin* apply regarding the tangential
stiffness and mixing of the shear modulus for different particle types.
The *mindlin_rescale* option uses the same form as *mindlin*\ , but the The *mindlin_rescale* option uses the same form as *mindlin*\ , but the
magnitude of the tangential displacement is re-scaled as the contact magnitude of the tangential displacement is re-scaled as the contact
@ -421,9 +467,32 @@ Here, :math:`t_{n-1}` indicates the value at the previous time
step. This rescaling accounts for the fact that a decrease in the step. This rescaling accounts for the fact that a decrease in the
contact area upon unloading leads to the contact being unable to contact area upon unloading leads to the contact being unable to
support the previous tangential loading, and spurious energy is support the previous tangential loading, and spurious energy is
created without the rescaling above (:ref:`Walton <WaltonPC>` ). See also created without the rescaling above (:ref:`Walton <WaltonPC>` ).
discussion in :ref:`Thornton et al, 2013 <Thornton2013>` , particularly
equation 18(b) of that work and associated discussion. .. note::
For *mindlin*, a decrease in the tangential force already occurs as the
contact unloads, due to the dependence of the tangential force on the normal
force described above. By re-scaling :math:`\xi`, *mindlin_rescale*
effectively re-scales the tangential force twice, i.e., proportionally to
:math:`a^2`. This peculiar behavior results from use of the accumulated
tangential displacement to characterize the contact history. Although
*mindlin_rescale* remains available for historic reasons and backward
compatibility purposes, it should be avoided in favor of *mindlin_rescale/force*.
The *mindlin_rescale/force* option uses the same form as *mindlin/force*,
but the magnitude of the tangential elastic force is re-scaled as the contact
unloads, i.e. if :math:`a < a_{t_{n-1}}`:
.. math::
\mathbf{F}_{te} = \mathbf{F}_{te, t_{n-1}} \frac{a}{a_{t_{n-1}}}
This approach provides a better approximation of the :ref:`Mindlin-Deresiewicz <Mindlin1953>`
laws and is more consistent than *mindlin_rescale*. See discussions in
:ref:`Thornton et al, 2013 <Thornton2013>`, particularly equation 18(b) of that
work and associated discussion, and :ref:`Agnolin and Roux, 2007 <AgnolinRoux2007>`,
particularly Appendix A.
---------- ----------
@ -460,7 +529,7 @@ exceeds a critical value:
.. math:: .. math::
\mathbf{F}_{roll} = min(\mu_{roll} F_{n,0}, \|\mathbf{F}_{roll,0}\|)\mathbf{k} \mathbf{F}_{roll} = \min(\mu_{roll} F_{n,0}, \|\mathbf{F}_{roll,0}\|)\mathbf{k}
Here, :math:`\mathbf{k} = \mathbf{v}_{roll}/\|\mathbf{v}_{roll}\|` is the direction of Here, :math:`\mathbf{k} = \mathbf{v}_{roll}/\|\mathbf{v}_{roll}\|` is the direction of
the pseudo-force. As with tangential displacement, the rolling the pseudo-force. As with tangential displacement, the rolling
@ -512,7 +581,7 @@ is then truncated according to:
.. math:: .. math::
\tau_{twist} = min(\mu_{twist} F_{n,0}, \tau_{twist,0}) \tau_{twist} = \min(\mu_{twist} F_{n,0}, \tau_{twist,0})
Similar to the sliding and rolling displacement, the angular Similar to the sliding and rolling displacement, the angular
displacement is rescaled so that it corresponds to the critical value displacement is rescaled so that it corresponds to the critical value
@ -763,3 +832,15 @@ Technology, 233, 30-46.
.. _WaltonPC: .. _WaltonPC:
**(Otis R. Walton)** Walton, O.R., Personal Communication **(Otis R. Walton)** Walton, O.R., Personal Communication
.. _Mindlin1953:
**(Mindlin and Deresiewicz, 1953)** Mindlin, R.D., & Deresiewicz, H (1953).
Elastic Spheres in Contact under Varying Oblique Force.
J. Appl. Mech., ASME 20, 327-344.
.. _AgnolinRoux2007:
**(Agnolin and Roux 2007)** Agnolin, I. & Roux, J-N. (2007).
Internal states of model isotropic granular packings.
I. Assembling process, geometry, and contact networks. Phys. Rev. E, 76, 061302.

3
doc/utils/sphinx-config/false_positives.txt
View File

@ -43,6 +43,7 @@ Afshar
agilio agilio
Agilio Agilio
agni agni
Agnolin
Ai Ai
Aidan Aidan
aij aij
@ -599,6 +600,7 @@ Dequidt
der der
dereference dereference
derekt derekt
Deresiewicz
Derjagin Derjagin
Derjaguin Derjaguin
Derlet Derlet
@ -2219,6 +2221,7 @@ oxdna
oxrna oxrna
oxDNA oxDNA
oxRNA oxRNA
packings
padua padua
Padua Padua
palegoldenrod palegoldenrod

161
src/GRANULAR/pair_granular.cpp
View File

@ -55,7 +55,8 @@ using namespace MathSpecial;
enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR};
enum {VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI}; enum {VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI};
enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY,
TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE}; TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE,
TANGENTIAL_MINDLIN_FORCE, TANGENTIAL_MINDLIN_RESCALE_FORCE};
enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL}; enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL};
enum {ROLL_NONE, ROLL_SDS}; enum {ROLL_NONE, ROLL_SDS};
@ -175,7 +176,8 @@ void PairGranular::compute(int eflag, int vflag)
double signtwist, magtwist, magtortwist, Mtcrit; double signtwist, magtwist, magtortwist, Mtcrit;
double tortwist1, tortwist2, tortwist3; double tortwist1, tortwist2, tortwist3;
double shrmag,rsht; double shrmag,rsht,prjmag;
bool frameupdate;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist,*jlist,*numneigh,**firstneigh;
int *touch,**firsttouch; int *touch,**firsttouch;
double *history,*allhistory,**firsthistory; double *history,*allhistory,**firsthistory;
@ -372,6 +374,12 @@ void PairGranular::compute(int eflag, int vflag)
// tangential force, including history effects // tangential force, including history effects
//**************************************** //****************************************
// For linear, mindlin, mindlin_rescale:
// history = cumulative tangential displacement
//
// For mindlin/force, mindlin_rescale/force:
// history = cumulative tangential elastic force
// tangential component // tangential component
vt1 = vr1 - vn1; vt1 = vr1 - vn1;
vt2 = vr2 - vn2; vt2 = vr2 - vn2;
@ -414,12 +422,15 @@ void PairGranular::compute(int eflag, int vflag)
damp_normal_prefactor; damp_normal_prefactor;
if (tangential_history) { if (tangential_history) {
if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN) { if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN ||
tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_FORCE) {
k_tangential *= a; k_tangential *= a;
} else if (tangential_model[itype][jtype] == } else if (tangential_model[itype][jtype] ==
TANGENTIAL_MINDLIN_RESCALE) { TANGENTIAL_MINDLIN_RESCALE ||
tangential_model[itype][jtype] ==
TANGENTIAL_MINDLIN_RESCALE_FORCE) {
k_tangential *= a; k_tangential *= a;
// on unloading, rescale the shear displacements // on unloading, rescale the shear displacements/force
if (a < history[3]) { if (a < history[3]) {
double factor = a/history[3]; double factor = a/history[3];
history[0] *= factor; history[0] *= factor;
@ -427,37 +438,66 @@ void PairGranular::compute(int eflag, int vflag)
history[2] *= factor; history[2] *= factor;
} }
} }
// rotate and update displacements. // rotate and update displacements / force.
// see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235
if (historyupdate) { if (historyupdate) {
rsht = history[0]*nx + history[1]*ny + history[2]*nz; rsht = history[0]*nx + history[1]*ny + history[2]*nz;
if (fabs(rsht) < EPSILON) rsht = 0; if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_FORCE ||
if (rsht > 0) { tangential_model[itype][jtype] ==
TANGENTIAL_MINDLIN_RESCALE_FORCE)
frameupdate = fabs(rsht) < EPSILON*Fscrit;
else
frameupdate = fabs(rsht)*k_tangential < EPSILON*Fscrit;
if (frameupdate) {
shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + shrmag = sqrt(history[0]*history[0] + history[1]*history[1] +
history[2]*history[2]); history[2]*history[2]);
// if rsht == shrmag, contacting pair has rotated 90 deg // projection
// in one step, in which case you deserve a crash!
scalefac = shrmag/(shrmag - rsht);
history[0] -= rsht*nx; history[0] -= rsht*nx;
history[1] -= rsht*ny; history[1] -= rsht*ny;
history[2] -= rsht*nz; history[2] -= rsht*nz;
// also rescale to preserve magnitude // also rescale to preserve magnitude
prjmag = sqrt(history[0]*history[0] + history[1]*history[1] +
history[2]*history[2]);
if (prjmag > 0) scalefac = shrmag/prjmag;
else scalefac = 0;
history[0] *= scalefac; history[0] *= scalefac;
history[1] *= scalefac; history[1] *= scalefac;
history[2] *= scalefac; history[2] *= scalefac;
} }
// update history // update history
history[0] += vtr1*dt; if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY ||
history[1] += vtr2*dt; tangential_model[itype][jtype] == TANGENTIAL_MINDLIN ||
history[2] += vtr3*dt; tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) // tangential displacement
history[0] += vtr1*dt;
history[1] += vtr2*dt;
history[2] += vtr3*dt;
} else {
// tangential force
// see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46
history[0] -= k_tangential*vtr1*dt;
history[1] -= k_tangential*vtr2*dt;
history[2] -= k_tangential*vtr3*dt;
}
if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE ||
tangential_model[itype][jtype] ==
TANGENTIAL_MINDLIN_RESCALE_FORCE)
history[3] = a; history[3] = a;
} }
// tangential forces = history + tangential velocity damping // tangential forces = history + tangential velocity damping
fs1 = -k_tangential*history[0] - damp_tangential*vtr1; if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY ||
fs2 = -k_tangential*history[1] - damp_tangential*vtr2; tangential_model[itype][jtype] == TANGENTIAL_MINDLIN ||
fs3 = -k_tangential*history[2] - damp_tangential*vtr3; tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
fs1 = -k_tangential*history[0] - damp_tangential*vtr1;
fs2 = -k_tangential*history[1] - damp_tangential*vtr2;
fs3 = -k_tangential*history[2] - damp_tangential*vtr3;
} else {
fs1 = history[0] - damp_tangential*vtr1;
fs2 = history[1] - damp_tangential*vtr2;
fs3 = history[2] - damp_tangential*vtr3;
}
// rescale frictional displacements and forces if needed // rescale frictional displacements and forces if needed
fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
@ -465,12 +505,21 @@ void PairGranular::compute(int eflag, int vflag)
shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + shrmag = sqrt(history[0]*history[0] + history[1]*history[1] +
history[2]*history[2]); history[2]*history[2]);
if (shrmag != 0.0) { if (shrmag != 0.0) {
history[0] = -1.0/k_tangential*(Fscrit*fs1/fs + if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY ||
damp_tangential*vtr1); tangential_model[itype][jtype] == TANGENTIAL_MINDLIN ||
history[1] = -1.0/k_tangential*(Fscrit*fs2/fs + tangential_model[itype][jtype] ==
damp_tangential*vtr2); TANGENTIAL_MINDLIN_RESCALE) {
history[2] = -1.0/k_tangential*(Fscrit*fs3/fs + history[0] = -1.0/k_tangential*(Fscrit*fs1/fs +
damp_tangential*vtr3); damp_tangential*vtr1);
history[1] = -1.0/k_tangential*(Fscrit*fs2/fs +
damp_tangential*vtr2);
history[2] = -1.0/k_tangential*(Fscrit*fs3/fs +
damp_tangential*vtr3);
} else {
history[0] = Fscrit*fs1/fs + damp_tangential*vtr1;
history[1] = Fscrit*fs2/fs + damp_tangential*vtr2;
history[2] = Fscrit*fs3/fs + damp_tangential*vtr3;
}
fs1 *= Fscrit/fs; fs1 *= Fscrit/fs;
fs2 *= Fscrit/fs; fs2 *= Fscrit/fs;
fs3 *= Fscrit/fs; fs3 *= Fscrit/fs;
@ -512,18 +561,27 @@ void PairGranular::compute(int eflag, int vflag)
int rhist1 = rhist0 + 1; int rhist1 = rhist0 + 1;
int rhist2 = rhist1 + 1; int rhist2 = rhist1 + 1;
rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; k_roll = roll_coeffs[itype][jtype][0];
damp_roll = roll_coeffs[itype][jtype][1];
Frcrit = roll_coeffs[itype][jtype][2] * Fncrit;
if (historyupdate) { if (historyupdate) {
if (fabs(rolldotn) < EPSILON) rolldotn = 0; rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz;
if (rolldotn > 0) { // rotate into tangential plane frameupdate = fabs(rolldotn)*k_roll < EPSILON*Frcrit;
if (frameupdate) { // rotate into tangential plane
rollmag = sqrt(history[rhist0]*history[rhist0] + rollmag = sqrt(history[rhist0]*history[rhist0] +
history[rhist1]*history[rhist1] + history[rhist1]*history[rhist1] +
history[rhist2]*history[rhist2]); history[rhist2]*history[rhist2]);
scalefac = rollmag/(rollmag - rolldotn); // projection
history[rhist0] -= rolldotn*nx; history[rhist0] -= rolldotn*nx;
history[rhist1] -= rolldotn*ny; history[rhist1] -= rolldotn*ny;
history[rhist2] -= rolldotn*nz; history[rhist2] -= rolldotn*nz;
// also rescale to preserve magnitude // also rescale to preserve magnitude
prjmag = sqrt(history[rhist0]*history[rhist0] +
history[rhist1]*history[rhist1] +
history[rhist2]*history[rhist2]);
if (prjmag > 0) scalefac = rollmag/prjmag;
else scalefac = 0;
history[rhist0] *= scalefac; history[rhist0] *= scalefac;
history[rhist1] *= scalefac; history[rhist1] *= scalefac;
history[rhist2] *= scalefac; history[rhist2] *= scalefac;
@ -533,14 +591,11 @@ void PairGranular::compute(int eflag, int vflag)
history[rhist2] += vrl3*dt; history[rhist2] += vrl3*dt;
} }
k_roll = roll_coeffs[itype][jtype][0];
damp_roll = roll_coeffs[itype][jtype][1];
fr1 = -k_roll*history[rhist0] - damp_roll*vrl1; fr1 = -k_roll*history[rhist0] - damp_roll*vrl1;
fr2 = -k_roll*history[rhist1] - damp_roll*vrl2; fr2 = -k_roll*history[rhist1] - damp_roll*vrl2;
fr3 = -k_roll*history[rhist2] - damp_roll*vrl3; fr3 = -k_roll*history[rhist2] - damp_roll*vrl3;
// rescale frictional displacements and forces if needed // rescale frictional displacements and forces if needed
Frcrit = roll_coeffs[itype][jtype][2] * Fncrit;
fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3);
if (fr > Frcrit) { if (fr > Frcrit) {
@ -734,7 +789,8 @@ void PairGranular::coeff(int narg, char **arg)
//Defaults //Defaults
normal_model_one = tangential_model_one = -1; normal_model_one = tangential_model_one = -1;
roll_model_one = twist_model_one = 0; roll_model_one = ROLL_NONE;
twist_model_one = TWIST_NONE;
damping_model_one = VISCOELASTIC; damping_model_one = VISCOELASTIC;
int iarg = 2; int iarg = 2;
@ -820,7 +876,9 @@ void PairGranular::coeff(int narg, char **arg)
iarg += 4; iarg += 4;
} else if ((strcmp(arg[iarg+1], "linear_history") == 0) || } else if ((strcmp(arg[iarg+1], "linear_history") == 0) ||
(strcmp(arg[iarg+1], "mindlin") == 0) || (strcmp(arg[iarg+1], "mindlin") == 0) ||
(strcmp(arg[iarg+1], "mindlin_rescale") == 0)) { (strcmp(arg[iarg+1], "mindlin_rescale") == 0) ||
(strcmp(arg[iarg+1], "mindlin/force") == 0) ||
(strcmp(arg[iarg+1], "mindlin_rescale/force") == 0)) {
if (iarg + 4 >= narg) if (iarg + 4 >= narg)
error->all(FLERR,"Illegal pair_coeff command, " error->all(FLERR,"Illegal pair_coeff command, "
"not enough parameters provided for tangential model"); "not enough parameters provided for tangential model");
@ -830,9 +888,15 @@ void PairGranular::coeff(int narg, char **arg)
tangential_model_one = TANGENTIAL_MINDLIN; tangential_model_one = TANGENTIAL_MINDLIN;
else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0) else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0)
tangential_model_one = TANGENTIAL_MINDLIN_RESCALE; tangential_model_one = TANGENTIAL_MINDLIN_RESCALE;
else if (strcmp(arg[iarg+1], "mindlin/force") == 0)
tangential_model_one = TANGENTIAL_MINDLIN_FORCE;
else if (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0)
tangential_model_one = TANGENTIAL_MINDLIN_RESCALE_FORCE;
tangential_history = 1; tangential_history = 1;
if ((tangential_model_one == TANGENTIAL_MINDLIN || if ((tangential_model_one == TANGENTIAL_MINDLIN ||
tangential_model_one == TANGENTIAL_MINDLIN_RESCALE) && tangential_model_one == TANGENTIAL_MINDLIN_RESCALE ||
tangential_model_one == TANGENTIAL_MINDLIN_FORCE ||
tangential_model_one == TANGENTIAL_MINDLIN_RESCALE_FORCE) &&
(strcmp(arg[iarg+2], "NULL") == 0)) { (strcmp(arg[iarg+2], "NULL") == 0)) {
if (normal_model_one == HERTZ || normal_model_one == HOOKE) { if (normal_model_one == HERTZ || normal_model_one == HOOKE) {
error->all(FLERR, "NULL setting for Mindlin tangential " error->all(FLERR, "NULL setting for Mindlin tangential "
@ -1014,7 +1078,8 @@ void PairGranular::init_style()
} }
for (int i = 1; i <= atom->ntypes; i++) for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++) for (int j = i; j <= atom->ntypes; j++)
if (tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE) { if (tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE ||
tangential_model[i][j] == TANGENTIAL_MINDLIN_RESCALE_FORCE) {
size_history += 1; size_history += 1;
roll_history_index += 1; roll_history_index += 1;
twist_history_index += 1; twist_history_index += 1;
@ -1484,6 +1549,12 @@ double PairGranular::single(int i, int j, int itype, int jtype,
// tangential force, including history effects // tangential force, including history effects
//**************************************** //****************************************
// For linear, mindlin, mindlin_rescale:
// history = cumulative tangential displacement
//
// For mindlin/force, mindlin_rescale/force:
// history = cumulative tangential elastic force
// tangential component // tangential component
vt1 = vr1 - vn1; vt1 = vr1 - vn1;
vt2 = vr2 - vn2; vt2 = vr2 - vn2;
@ -1519,9 +1590,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor; damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal_prefactor;
if (tangential_history) { if (tangential_history) {
if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN) { if (tangential_model[itype][jtype] != TANGENTIAL_HISTORY) {
k_tangential *= a;
} else if (tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
k_tangential *= a; k_tangential *= a;
} }
@ -1529,9 +1598,17 @@ double PairGranular::single(int i, int j, int itype, int jtype,
history[2]*history[2]); history[2]*history[2]);
// tangential forces = history + tangential velocity damping // tangential forces = history + tangential velocity damping
fs1 = -k_tangential*history[0] - damp_tangential*vtr1; if (tangential_model[itype][jtype] == TANGENTIAL_HISTORY ||
fs2 = -k_tangential*history[1] - damp_tangential*vtr2; tangential_model[itype][jtype] == TANGENTIAL_MINDLIN ||
fs3 = -k_tangential*history[2] - damp_tangential*vtr3; tangential_model[itype][jtype] == TANGENTIAL_MINDLIN_RESCALE) {
fs1 = -k_tangential*history[0] - damp_tangential*vtr1;
fs2 = -k_tangential*history[1] - damp_tangential*vtr2;
fs3 = -k_tangential*history[2] - damp_tangential*vtr3;
} else {
fs1 = history[0] - damp_tangential*vtr1;
fs2 = history[1] - damp_tangential*vtr2;
fs3 = history[2] - damp_tangential*vtr3;
}
// rescale frictional forces if needed // rescale frictional forces if needed
fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
@ -1540,7 +1617,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
fs1 *= Fscrit/fs; fs1 *= Fscrit/fs;
fs2 *= Fscrit/fs; fs2 *= Fscrit/fs;
fs3 *= Fscrit/fs; fs3 *= Fscrit/fs;
fs *= Fscrit/fs; fs *= Fscrit/fs;
} else fs1 = fs2 = fs3 = fs = 0.0; } else fs1 = fs2 = fs3 = fs = 0.0;
} }

44
src/USER-REACTION/fix_bond_react.cpp
View File

@ -33,6 +33,7 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu)
#include "neigh_list.h" #include "neigh_list.h"
#include "neigh_request.h" #include "neigh_request.h"
#include "random_mars.h" #include "random_mars.h"
#include "reset_mol_ids.h"
#include "molecule.h" #include "molecule.h"
#include "group.h" #include "group.h"
#include "citeme.h" #include "citeme.h"
@ -92,6 +93,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
fix1 = NULL; fix1 = NULL;
fix2 = NULL; fix2 = NULL;
fix3 = NULL; fix3 = NULL;
reset_mol_ids = NULL;
if (narg < 8) error->all(FLERR,"Illegal fix bond/react command: " if (narg < 8) error->all(FLERR,"Illegal fix bond/react command: "
"too few arguments"); "too few arguments");
@ -144,7 +146,8 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
int iarg = 3; int iarg = 3;
stabilization_flag = 0; stabilization_flag = 0;
int num_common_keywords = 1; reset_mol_ids_flag = 1;
int num_common_keywords = 2;
for (int m = 0; m < num_common_keywords; m++) { for (int m = 0; m < num_common_keywords; m++) {
if (strcmp(arg[iarg],"stabilization") == 0) { if (strcmp(arg[iarg],"stabilization") == 0) {
if (strcmp(arg[iarg+1],"no") == 0) { if (strcmp(arg[iarg+1],"no") == 0) {
@ -162,11 +165,23 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
nve_limit_xmax = arg[iarg+3]; nve_limit_xmax = arg[iarg+3];
iarg += 4; iarg += 4;
} }
} else if (strcmp(arg[iarg],"reset_mol_ids") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: "
"'reset_mol_ids' keyword has too few arguments");
if (strcmp(arg[iarg+1],"yes") == 0) ; // default
if (strcmp(arg[iarg+1],"no") == 0) reset_mol_ids_flag = 0;
iarg += 2;
} else if (strcmp(arg[iarg],"react") == 0) { } else if (strcmp(arg[iarg],"react") == 0) {
break; break;
} else error->all(FLERR,"Illegal fix bond/react command: unknown keyword"); } else error->all(FLERR,"Illegal fix bond/react command: unknown keyword");
} }
if (reset_mol_ids_flag) {
delete reset_mol_ids;
reset_mol_ids = new ResetMolIDs(lmp);
reset_mol_ids->create_computes(id,group->names[igroup]);
}
// set up common variables as vectors of length 'nreacts' // set up common variables as vectors of length 'nreacts'
// nevery, cutoff, onemol, twomol, superimpose file // nevery, cutoff, onemol, twomol, superimpose file
@ -229,11 +244,11 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
int n = strlen(arg[iarg]) + 1; int n = strlen(arg[iarg]) + 1;
if (n > MAXLINE) error->all(FLERR,"Reaction name (react-ID) is too long (limit: 256 characters)"); if (n > MAXLINE) error->all(FLERR,"Reaction name (react-ID) is too long (limit: 256 characters)");
strncpy(rxn_name[rxn],arg[iarg++],n); strcpy(rxn_name[rxn],arg[iarg++]);
int igroup = group->find(arg[iarg++]); int groupid = group->find(arg[iarg++]);
if (igroup == -1) error->all(FLERR,"Could not find fix group ID"); if (groupid == -1) error->all(FLERR,"Could not find fix group ID");
groupbits[rxn] = group->bitmask[igroup]; groupbits[rxn] = group->bitmask[groupid];
if (strncmp(arg[iarg],"v_",2) == 0) { if (strncmp(arg[iarg],"v_",2) == 0) {
n = strlen(&arg[iarg][2]) + 1; n = strlen(&arg[iarg][2]) + 1;
@ -499,6 +514,8 @@ FixBondReact::~FixBondReact()
} }
delete [] random; delete [] random;
delete reset_mol_ids;
memory->destroy(partner); memory->destroy(partner);
memory->destroy(finalpartner); memory->destroy(finalpartner);
memory->destroy(ncreate); memory->destroy(ncreate);
@ -623,9 +640,9 @@ void FixBondReact::post_constructor()
group->assign(cmd); group->assign(cmd);
if (stabilization_flag == 1) { if (stabilization_flag == 1) {
int igroup = group->find(exclude_group); int groupid = group->find(exclude_group);
// create exclude_group if not already existing, or use as parent group if static // create exclude_group if not already existing, or use as parent group if static
if (igroup == -1 || group->dynamic[igroup] == 0) { if (groupid == -1 || group->dynamic[groupid] == 0) {
// create stabilization per-atom property // create stabilization per-atom property
cmd = std::string("bond_react_stabilization_internal"); cmd = std::string("bond_react_stabilization_internal");
id_fix3 = new char[cmd.size()+1]; id_fix3 = new char[cmd.size()+1];
@ -655,7 +672,7 @@ void FixBondReact::post_constructor()
strcat(exclude_group,"_REACT"); strcat(exclude_group,"_REACT");
group->find_or_create(exclude_group); group->find_or_create(exclude_group);
if (igroup == -1) if (groupid == -1)
cmd = fmt::format("{} dynamic all property statted_tags",exclude_group); cmd = fmt::format("{} dynamic all property statted_tags",exclude_group);
else else
cmd = fmt::format("{} dynamic {} property statted_tags",exclude_group,exclude_PARENT_group); cmd = fmt::format("{} dynamic {} property statted_tags",exclude_group,exclude_PARENT_group);
@ -2061,7 +2078,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
for (int i = 0; i < twomol->natoms; i++) { for (int i = 0; i < twomol->natoms; i++) {
if (twomol->type[i] != onemol->type[equivalences[i][1][myrxn]-1] && landlocked_atoms[i][myrxn] == 0) { if (twomol->type[i] != onemol->type[equivalences[i][1][myrxn]-1] && landlocked_atoms[i][myrxn] == 0) {
char str[128]; char str[128];
snprintf(str,128,"Bond/react: Atom affected by reaction %s too close to template edge",rxn_name[myrxn]); snprintf(str,128,"Bond/react: Atom type affected by reaction %s too close to template edge",rxn_name[myrxn]);
error->all(FLERR,str); error->all(FLERR,str);
} }
} }
@ -2080,7 +2097,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
if (onemol_batom == equivalences[twomol_atomj-1][1][myrxn]) { if (onemol_batom == equivalences[twomol_atomj-1][1][myrxn]) {
if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomi-1][m]) { if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomi-1][m]) {
char str[128]; char str[128];
snprintf(str,128,"Bond/react: Atom affected by reaction %s too close to template edge",rxn_name[myrxn]); snprintf(str,128,"Bond/react: Bond type affected by reaction %s too close to template edge",rxn_name[myrxn]);
error->all(FLERR,str); error->all(FLERR,str);
} }
} }
@ -2092,7 +2109,7 @@ void FixBondReact::find_landlocked_atoms(int myrxn)
if (onemol_batom == equivalences[i][1][myrxn]) { if (onemol_batom == equivalences[i][1][myrxn]) {
if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomj-1][m]) { if (twomol->bond_type[i][j] != onemol->bond_type[onemol_atomj-1][m]) {
char str[128]; char str[128];
snprintf(str,128,"Bond/react: Atom affected by reaction %s too close to template edge",rxn_name[myrxn]); snprintf(str,128,"Bond/react: Bond type affected by reaction %s too close to template edge",rxn_name[myrxn]);
error->all(FLERR,str); error->all(FLERR,str);
} }
} }
@ -2503,7 +2520,7 @@ void FixBondReact::ghost_glovecast()
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
update charges, types, special lists and all topology update molecule IDs, charges, types, special lists and all topology
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void FixBondReact::update_everything() void FixBondReact::update_everything()
@ -3042,6 +3059,9 @@ void FixBondReact::update_everything()
atom->natoms -= ndel; atom->natoms -= ndel;
// done deleting atoms // done deleting atoms
// reset mol ids
if (reset_mol_ids_flag) reset_mol_ids->reset();
// something to think about: this could done much more concisely if // something to think about: this could done much more concisely if
// all atom-level info (bond,angles, etc...) were kinda inherited from a common data struct --JG // all atom-level info (bond,angles, etc...) were kinda inherited from a common data struct --JG

6
src/USER-REACTION/fix_bond_react.h
View File

@ -61,6 +61,7 @@ class FixBondReact : public Fix {
int *max_rxn,*nlocalskips,*nghostlyskips; int *max_rxn,*nlocalskips,*nghostlyskips;
tagint lastcheck; tagint lastcheck;
int stabilization_flag; int stabilization_flag;
int reset_mol_ids_flag;
int custom_exclude_flag; int custom_exclude_flag;
int *stabilize_steps_flag; int *stabilize_steps_flag;
int *update_edges_flag; int *update_edges_flag;
@ -93,6 +94,7 @@ class FixBondReact : public Fix {
class RanMars **random; // random number for 'prob' keyword class RanMars **random; // random number for 'prob' keyword
class RanMars **rrhandom; // random number for Arrhenius constraint class RanMars **rrhandom; // random number for Arrhenius constraint
class NeighList *list; class NeighList *list;
class ResetMolIDs *reset_mol_ids; // class for resetting mol IDs
int *reacted_mol,*unreacted_mol; int *reacted_mol,*unreacted_mol;
int *limit_duration; // indicates how long to relax int *limit_duration; // indicates how long to relax
@ -240,8 +242,10 @@ E: Bond/react: Invalid template atom ID in map file
Atom IDs in molecule templates range from 1 to the number of atoms in the template. Atom IDs in molecule templates range from 1 to the number of atoms in the template.
E or W: Bond/react: Atom affected by reaction %s too close to template edge E or W: Bond/react: Atom affected by reaction %s too close to template edge
Bond/react: Atom type affected by reaction %s too close to template edge
Bond/react: Bond type affected by reaction %s too close to template edge
This means an atom which changes type or connectivity during the This means an atom (or bond) that changes type or connectivity during the
reaction is too close to an 'edge' atom defined in the map file. This reaction is too close to an 'edge' atom defined in the map file. This
could cause incorrect assignment of bonds, angle, etc. Generally, this could cause incorrect assignment of bonds, angle, etc. Generally, this
means you must include more atoms in your templates, such that there means you must include more atoms in your templates, such that there

134
src/reset_mol_ids.cpp
View File

@ -17,7 +17,6 @@
#include "reset_mol_ids.h" #include "reset_mol_ids.h"
#include <mpi.h> #include <mpi.h>
#include <string>
#include "atom.h" #include "atom.h"
#include "domain.h" #include "domain.h"
#include "comm.h" #include "comm.h"
@ -33,10 +32,32 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
ResetMolIDs::ResetMolIDs(LAMMPS *lmp) : Pointers(lmp) {} ResetMolIDs::ResetMolIDs(LAMMPS *lmp) : Pointers(lmp) {
cfa = NULL;
cca = NULL;
// default settings
compressflag = 1;
singleflag = 0;
offset = -1;
idfrag.clear();
idchunk.clear();
}
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
ResetMolIDs::~ResetMolIDs()
{
if (!idfrag.empty()) modify->delete_compute(idfrag);
if (compressflag && !idchunk.empty()) modify->delete_compute(idchunk);
}
/* ----------------------------------------------------------------------
called as reset_mol_ids command in input script
------------------------------------------------------------------------- */
void ResetMolIDs::command(int narg, char **arg) void ResetMolIDs::command(int narg, char **arg)
{ {
if (domain->box_exist == 0) if (domain->box_exist == 0)
@ -49,13 +70,7 @@ void ResetMolIDs::command(int narg, char **arg)
// process args // process args
if (narg < 1) error->all(FLERR,"Illegal reset_mol_ids command"); if (narg < 1) error->all(FLERR,"Illegal reset_mol_ids command");
int igroup = group->find(arg[0]); char *groupid = arg[0];
if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID");
int groupbit = group->bitmask[igroup];
int compressflag = 1;
int singleflag = 0;
tagint offset = -1;
int iarg = 1; int iarg = 1;
while (iarg < narg) { while (iarg < narg) {
@ -86,20 +101,6 @@ void ResetMolIDs::command(int narg, char **arg)
MPI_Barrier(world); MPI_Barrier(world);
double time1 = MPI_Wtime(); double time1 = MPI_Wtime();
// create instances of compute fragment/atom, compute reduce (if needed),
// and compute chunk/atom. all use the group-ID for this command
const std::string idfrag = "reset_mol_ids_FRAGMENT_ATOM";
if (singleflag)
modify->add_compute(fmt::format("{} {} fragment/atom single yes",idfrag,arg[0]));
else
modify->add_compute(fmt::format("{} {} fragment/atom single no",idfrag,arg[0]));
const std::string idchunk = "reset_mol_ids_CHUNK_ATOM";
if (compressflag)
modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes",
idchunk,arg[0]));
// initialize system since comm->borders() will be invoked // initialize system since comm->borders() will be invoked
lmp->init(); lmp->init();
@ -116,12 +117,73 @@ void ResetMolIDs::command(int narg, char **arg)
comm->borders(); comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
// create computes
create_computes((char *) "COMMAND",groupid);
// reset molecule IDs
reset();
// total time
MPI_Barrier(world);
if (comm->me == 0) {
if (nchunk < 0)
utils::logmesg(lmp,fmt::format(" number of new molecule IDs = unknown\n"));
else
utils::logmesg(lmp,fmt::format(" number of new molecule IDs = {}\n",nchunk));
utils::logmesg(lmp,fmt::format(" reset_mol_ids CPU = {:.3f} seconds\n",
MPI_Wtime()-time1));
}
}
/* ----------------------------------------------------------------------
create computes used by reset_mol_ids
------------------------------------------------------------------------- */
void ResetMolIDs::create_computes(char *fixid, char *groupid)
{
int igroup = group->find(groupid);
if (igroup == -1) error->all(FLERR,"Could not find reset_mol_ids group ID");
groupbit = group->bitmask[igroup];
// create instances of compute fragment/atom, compute reduce (if needed),
// and compute chunk/atom. all use the group-ID for this command.
// 'fixid' allows for creating independent instances of the computes
idfrag = fmt::format("{}_reset_mol_ids_FRAGMENT_ATOM",fixid);
if (singleflag)
modify->add_compute(fmt::format("{} {} fragment/atom single yes",idfrag,groupid));
else
modify->add_compute(fmt::format("{} {} fragment/atom single no",idfrag,groupid));
idchunk = fmt::format("{}_reset_mol_ids_CHUNK_ATOM",fixid);
if (compressflag)
modify->add_compute(fmt::format("{} {} chunk/atom molecule compress yes",
idchunk,groupid));
int icompute = modify->find_compute(idfrag);
cfa = (ComputeFragmentAtom *) modify->compute[icompute];
if (compressflag) {
icompute = modify->find_compute(idchunk);
cca = (ComputeChunkAtom *) modify->compute[icompute];
}
}
/* ----------------------------------------------------------------------
called from command() and directly from fixes that update molecules
------------------------------------------------------------------------- */
void ResetMolIDs::reset()
{
// invoke peratom method of compute fragment/atom // invoke peratom method of compute fragment/atom
// walks bond connectivity and assigns each atom a fragment ID // walks bond connectivity and assigns each atom a fragment ID
// if singleflag = 0, atoms w/out bonds will be assigned fragID = 0 // if singleflag = 0, atoms w/out bonds will be assigned fragID = 0
int icompute = modify->find_compute(idfrag);
ComputeFragmentAtom *cfa = (ComputeFragmentAtom *) modify->compute[icompute];
cfa->compute_peratom(); cfa->compute_peratom();
double *fragIDs = cfa->vector_atom; double *fragIDs = cfa->vector_atom;
@ -138,15 +200,13 @@ void ResetMolIDs::command(int narg, char **arg)
// if compressflag = 0, done // if compressflag = 0, done
// set nchunk = -1 since cannot easily determine # of new molecule IDs // set nchunk = -1 since cannot easily determine # of new molecule IDs
int nchunk = -1; nchunk = -1;
// if compressflag = 1, invoke peratom method of compute chunk/atom // if compressflag = 1, invoke peratom method of compute chunk/atom
// will compress new molecule IDs to be contiguous 1 to Nmol // will compress new molecule IDs to be contiguous 1 to Nmol
// NOTE: use of compute chunk/atom limits Nmol to a 32-bit int // NOTE: use of compute chunk/atom limits Nmol to a 32-bit int
if (compressflag) { if (compressflag) {
icompute = modify->find_compute(idchunk);
ComputeChunkAtom *cca = (ComputeChunkAtom *) modify->compute[icompute];
cca->compute_peratom(); cca->compute_peratom();
double *chunkIDs = cca->vector_atom; double *chunkIDs = cca->vector_atom;
nchunk = cca->nchunk; nchunk = cca->nchunk;
@ -193,22 +253,4 @@ void ResetMolIDs::command(int narg, char **arg)
} }
} }
} }
// clean up
modify->delete_compute(idfrag);
if (compressflag) modify->delete_compute(idchunk);
// total time
MPI_Barrier(world);
if (comm->me == 0) {
if (nchunk < 0)
utils::logmesg(lmp,fmt::format(" number of new molecule IDs = unknown\n"));
else
utils::logmesg(lmp,fmt::format(" number of new molecule IDs = {}\n",nchunk));
utils::logmesg(lmp,fmt::format(" reset_mol_ids CPU = {:.3f} seconds\n",
MPI_Wtime()-time1));
}
} }

15
src/reset_mol_ids.h
View File

@ -21,13 +21,28 @@ CommandStyle(reset_mol_ids,ResetMolIDs)
#define LMP_RESET_MOL_IDS_H #define LMP_RESET_MOL_IDS_H
#include "pointers.h" #include "pointers.h"
#include <string>
namespace LAMMPS_NS { namespace LAMMPS_NS {
class ResetMolIDs : protected Pointers { class ResetMolIDs : protected Pointers {
public: public:
ResetMolIDs(class LAMMPS *); ResetMolIDs(class LAMMPS *);
~ResetMolIDs();
void command(int, char **); void command(int, char **);
void create_computes(char *, char *);
void reset();
private:
std::string idfrag, idchunk;
int nchunk;
int groupbit;
int compressflag; // 1 = contiguous values for new IDs
int singleflag; // 0 = mol IDs of single atoms set to 0
tagint offset; // offset for contiguous mol ID values
class ComputeFragmentAtom *cfa;
class ComputeChunkAtom *cca;
}; };
} }

32
src/utils.cpp
View File

@ -352,7 +352,7 @@ tagint utils::tnumeric(const char *file, int line, const char *str,
Return string without leading or trailing whitespace Return string without leading or trailing whitespace
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
std::string utils::trim(const std::string & line) { std::string utils::trim(const std::string &line) {
int beg = re_match(line.c_str(),"\\S+"); int beg = re_match(line.c_str(),"\\S+");
int end = re_match(line.c_str(),"\\s+$"); int end = re_match(line.c_str(),"\\s+$");
if (beg < 0) beg = 0; if (beg < 0) beg = 0;
@ -365,7 +365,7 @@ std::string utils::trim(const std::string & line) {
Return string without trailing # comment Return string without trailing # comment
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
std::string utils::trim_comment(const std::string & line) { std::string utils::trim_comment(const std::string &line) {
auto end = line.find_first_of("#"); auto end = line.find_first_of("#");
if (end != std::string::npos) { if (end != std::string::npos) {
return line.substr(0, end); return line.substr(0, end);
@ -377,7 +377,7 @@ std::string utils::trim_comment(const std::string & line) {
return number of words return number of words
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
size_t utils::count_words(const char * text) { size_t utils::count_words(const char *text) {
size_t count = 0; size_t count = 0;
const char * buf = text; const char * buf = text;
char c = *buf; char c = *buf;
@ -406,7 +406,7 @@ size_t utils::count_words(const char * text) {
return number of words return number of words
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
size_t utils::count_words(const std::string & text) { size_t utils::count_words(const std::string &text) {
return utils::count_words(text.c_str()); return utils::count_words(text.c_str());
} }
@ -414,7 +414,7 @@ size_t utils::count_words(const std::string & text) {
Return number of words Return number of words
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
size_t utils::count_words(const std::string & text, const std::string & separators) { size_t utils::count_words(const std::string &text, const std::string &separators) {
size_t count = 0; size_t count = 0;
size_t start = text.find_first_not_of(separators); size_t start = text.find_first_not_of(separators);
@ -435,7 +435,7 @@ size_t utils::count_words(const std::string & text, const std::string & separato
Trim comment from string and return number of words Trim comment from string and return number of words
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
size_t utils::trim_and_count_words(const std::string & text, const std::string & separators) { size_t utils::trim_and_count_words(const std::string &text, const std::string &separators) {
return utils::count_words(utils::trim_comment(text), separators); return utils::count_words(utils::trim_comment(text), separators);
} }
@ -526,7 +526,7 @@ std::vector<std::string> utils::split_words(const std::string &text)
Return whether string is a valid integer number Return whether string is a valid integer number
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
bool utils::is_integer(const std::string & str) { bool utils::is_integer(const std::string &str) {
if (str.size() == 0) { if (str.size() == 0) {
return false; return false;
} }
@ -542,7 +542,7 @@ bool utils::is_integer(const std::string & str) {
Return whether string is a valid floating-point number Return whether string is a valid floating-point number
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
bool utils::is_double(const std::string & str) { bool utils::is_double(const std::string &str) {
if (str.size() == 0) { if (str.size() == 0) {
return false; return false;
} }
@ -560,7 +560,7 @@ bool utils::is_double(const std::string & str) {
strip off leading part of path, return just the filename strip off leading part of path, return just the filename
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
std::string utils::path_basename(const std::string & path) { std::string utils::path_basename(const std::string &path) {
#if defined(_WIN32) #if defined(_WIN32)
size_t start = path.find_last_of("/\\"); size_t start = path.find_last_of("/\\");
#else #else
@ -580,7 +580,7 @@ std::string utils::path_basename(const std::string & path) {
join two paths join two paths
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
std::string utils::path_join(const std::string & a, const std::string & b) { std::string utils::path_join(const std::string &a, const std::string &b) {
#if defined(_WIN32) #if defined(_WIN32)
return fmt::format("{}\\{}", a, b); return fmt::format("{}\\{}", a, b);
#else #else
@ -592,7 +592,7 @@ std::string utils::path_join(const std::string & a, const std::string & b) {
try to open file for reading try to open file for reading
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
bool utils::file_is_readable(const std::string & path) { bool utils::file_is_readable(const std::string &path) {
FILE * fp = fopen(path.c_str(), "r"); FILE * fp = fopen(path.c_str(), "r");
if(fp) { if(fp) {
fclose(fp); fclose(fp);
@ -607,7 +607,7 @@ bool utils::file_is_readable(const std::string & path) {
specified specified
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
std::string utils::get_potential_file_path(const std::string& path) { std::string utils::get_potential_file_path(const std::string &path) {
std::string filepath = path; std::string filepath = path;
std::string filename = utils::path_basename(path); std::string filename = utils::path_basename(path);
@ -634,7 +634,7 @@ std::string utils::get_potential_file_path(const std::string& path) {
if it has a DATE field, return the following word if it has a DATE field, return the following word
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
std::string utils::get_potential_date(const std::string & path, const std::string & potential_name) { std::string utils::get_potential_date(const std::string &path, const std::string &potential_name) {
TextFileReader reader(path, potential_name); TextFileReader reader(path, potential_name);
reader.ignore_comments = false; reader.ignore_comments = false;
@ -657,7 +657,7 @@ std::string utils::get_potential_date(const std::string & path, const std::strin
if it has UNITS field, return following word if it has UNITS field, return following word
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
std::string utils::get_potential_units(const std::string & path, const std::string & potential_name) { std::string utils::get_potential_units(const std::string &path, const std::string &potential_name) {
TextFileReader reader(path, potential_name); TextFileReader reader(path, potential_name);
reader.ignore_comments = false; reader.ignore_comments = false;
@ -710,7 +710,7 @@ double utils::get_conversion_factor(const int property, const int conversion)
the strings "off" and "unlimited" result in -1.0; the strings "off" and "unlimited" result in -1.0;
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
double utils::timespec2seconds(const std::string & timespec) double utils::timespec2seconds(const std::string &timespec)
{ {
double vals[3]; double vals[3];
int i = 0; int i = 0;
@ -728,7 +728,7 @@ double utils::timespec2seconds(const std::string & timespec)
if (!values.has_next()) break; if (!values.has_next()) break;
vals[i] = values.next_int(); vals[i] = values.next_int();
} }
} catch (TokenizerException & e) { } catch (TokenizerException &e) {
return -1.0; return -1.0;
} }

130
src/utils.h
View File

@ -29,7 +29,7 @@ namespace LAMMPS_NS {
namespace utils { namespace utils {
/** \brief Match text against a simplified regex pattern /** Match text against a simplified regex pattern
* *
* \param text the text to be matched against the pattern * \param text the text to be matched against the pattern
* \param pattern the search pattern, which may contain regexp markers * \param pattern the search pattern, which may contain regexp markers
@ -37,14 +37,14 @@ namespace LAMMPS_NS {
*/ */
bool strmatch(const std::string &text, const std::string &pattern); bool strmatch(const std::string &text, const std::string &pattern);
/** \brief Send message to screen and logfile, if available /** Send message to screen and logfile, if available
* *
* \param lmp pointer to LAMMPS class instance * \param lmp pointer to LAMMPS class instance
* \param mesg message to be printed * \param mesg message to be printed
*/ */
void logmesg(LAMMPS *lmp, const std::string &mesg); void logmesg(LAMMPS *lmp, const std::string &mesg);
/** \brief return a string representing the current system error status /** return a string representing the current system error status
* *
* This is a wrapper around calling strerror(errno). * This is a wrapper around calling strerror(errno).
* *
@ -52,7 +52,7 @@ namespace LAMMPS_NS {
*/ */
std::string getsyserror(); std::string getsyserror();
/** \brief safe wrapper around fgets() which aborts on errors /** safe wrapper around fgets() which aborts on errors
* or EOF and prints a suitable error message to help debugging * or EOF and prints a suitable error message to help debugging
* *
* \param srcname name of the calling source file (from FLERR macro) * \param srcname name of the calling source file (from FLERR macro)
@ -66,7 +66,7 @@ namespace LAMMPS_NS {
void sfgets(const char *srcname, int srcline, char *s, int size, void sfgets(const char *srcname, int srcline, char *s, int size,
FILE *fp, const char *filename, Error *error); FILE *fp, const char *filename, Error *error);
/** \brief safe wrapper around fread() which aborts on errors /** safe wrapper around fread() which aborts on errors
* or EOF and prints a suitable error message to help debugging * or EOF and prints a suitable error message to help debugging
* *
* \param srcname name of the calling source file (from FLERR macro) * \param srcname name of the calling source file (from FLERR macro)
@ -81,7 +81,7 @@ namespace LAMMPS_NS {
void sfread(const char *srcname, int srcline, void *s, size_t size, void sfread(const char *srcname, int srcline, void *s, size_t size,
size_t num, FILE *fp, const char *filename, Error *error); size_t num, FILE *fp, const char *filename, Error *error);
/** \brief Report if a requested style is in a package or may have a typo /** Report if a requested style is in a package or may have a typo
* *
* \param style type of style that is to be checked for * \param style type of style that is to be checked for
* \param name name of style that was not found * \param name name of style that was not found
@ -91,8 +91,8 @@ namespace LAMMPS_NS {
std::string check_packages_for_style(const std::string &style, std::string check_packages_for_style(const std::string &style,
const std::string &name, LAMMPS *lmp); const std::string &name, LAMMPS *lmp);
/** \brief Convert a string to a floating point number while checking /** Convert a string to a floating point number while checking
if it is a valid floating point or integer number * if it is a valid floating point or integer number
* *
* \param file name of source file for error message * \param file name of source file for error message
* \param line in source file for error message * \param line in source file for error message
@ -104,8 +104,8 @@ namespace LAMMPS_NS {
double numeric(const char *file, int line, const char *str, double numeric(const char *file, int line, const char *str,
bool do_abort, LAMMPS *lmp); bool do_abort, LAMMPS *lmp);
/** \brief Convert a string to an integer number while checking /** Convert a string to an integer number while checking
if it is a valid integer number (regular int) * if it is a valid integer number (regular int)
* *
* \param file name of source file for error message * \param file name of source file for error message
* \param line in source file for error message * \param line in source file for error message
@ -117,8 +117,8 @@ namespace LAMMPS_NS {
int inumeric(const char *file, int line, const char *str, int inumeric(const char *file, int line, const char *str,
bool do_abort, LAMMPS *lmp); bool do_abort, LAMMPS *lmp);
/** \brief Convert a string to an integer number while checking /** Convert a string to an integer number while checking
if it is a valid integer number (bigint) * if it is a valid integer number (bigint)
* *
* \param file name of source file for error message * \param file name of source file for error message
* \param line in source file for error message * \param line in source file for error message
@ -130,8 +130,8 @@ namespace LAMMPS_NS {
bigint bnumeric(const char *file, int line, const char *str, bigint bnumeric(const char *file, int line, const char *str,
bool do_abort, LAMMPS *lmp); bool do_abort, LAMMPS *lmp);
/** \brief Convert a string to an integer number while checking /** Convert a string to an integer number while checking
if it is a valid integer number (tagint) * if it is a valid integer number (tagint)
* *
* \param file name of source file for error message * \param file name of source file for error message
* \param line in source file for error message * \param line in source file for error message
@ -143,54 +143,51 @@ namespace LAMMPS_NS {
tagint tnumeric(const char *file, int line, const char *str, tagint tnumeric(const char *file, int line, const char *str,
bool do_abort, LAMMPS *lmp); bool do_abort, LAMMPS *lmp);
/** /** Trim leading and trailing whitespace. Like TRIM() in Fortran.
* \brief Trim leading and trailing whitespace. Like TRIM() in Fortran. *
* \param line string that should be trimmed * \param line string that should be trimmed
* \return new string without whitespace (string) * \return new string without whitespace (string)
*/ */
std::string trim(const std::string &line); std::string trim(const std::string &line);
/** /** Trim anything from '#' onward
* \brief Trim anything from '#' onward *
* \param line string that should be trimmed * \param line string that should be trimmed
* \return new string without comment (string) * \return new string without comment (string)
*/ */
std::string trim_comment(const std::string &line); std::string trim_comment(const std::string &line);
/** /** Count words in string
* \brief Count words in string *
* \param text string that should be searched * \param text string that should be searched
* \param separators string containing characters that will be treated as whitespace * \param separators string containing characters that will be treated as whitespace
* \return number of words found * \return number of words found
*/ */
size_t count_words(const std::string & text, const std::string & separators); size_t count_words(const std::string &text, const std::string &separators);
/** /** Count words in string, ignore any whitespace matching " \t\r\n\f"
* \brief Count words in string, ignore any whitespace matching " \t\r\n\f" *
* \param text string that should be searched * \param text string that should be searched
* \param separators string containing characters that will be treated as whitespace
* \return number of words found * \return number of words found
*/ */
size_t count_words(const std::string & text); size_t count_words(const std::string &text);
/** /** Count words in C-string, ignore any whitespace matching " \t\r\n\f"
* \brief Count words in C-string, ignore any whitespace matching " \t\r\n\f" *
* \param text string that should be searched * \param text string that should be searched
* \param separators string containing characters that will be treated as whitespace
* \return number of words found * \return number of words found
*/ */
size_t count_words(const char * text); size_t count_words(const char *text);
/** /** Count words in a single line, trim anything from '#' onward
* \brief Count words in a single line, trim anything from '#' onward *
* \param text string that should be trimmed and searched * \param text string that should be trimmed and searched
* \param separators string containing characters that will be treated as whitespace * \param separators string containing characters that will be treated as whitespace
* \return number of words found * \return number of words found
*/ */
size_t trim_and_count_words(const std::string & text, const std::string & separators = " \t\r\n\f"); size_t trim_and_count_words(const std::string &text, const std::string &separators = " \t\r\n\f");
/** /** Take text and split into non-whitespace words.
* \brief Take text and split into non-whitespace words.
* *
* This can handle single and double quotes, escaped quotes, * This can handle single and double quotes, escaped quotes,
* and escaped codes within quotes, but due to using an STL * and escaped codes within quotes, but due to using an STL
@ -203,21 +200,21 @@ namespace LAMMPS_NS {
*/ */
std::vector<std::string> split_words(const std::string &text); std::vector<std::string> split_words(const std::string &text);
/** /** Check if string can be converted to valid integer
* \brief Check if string can be converted to valid integer *
* \param text string that should be checked * \param str string that should be checked
* \return true, if string contains valid integer, false otherwise * \return true, if string contains valid integer, false otherwise
*/ */
bool is_integer(const std::string & str); bool is_integer(const std::string &str);
/** /** Check if string can be converted to valid floating-point number
* \brief Check if string can be converted to valid floating-point number *
* \param text string that should be checked * \param str string that should be checked
* \return true, if string contains valid floating-point number, false otherwise * \return true, if string contains valid floating-point number, false otherwise
*/ */
bool is_double(const std::string & str); bool is_double(const std::string &str);
/** \brief try to detect pathname from FILE pointer. /** Try to detect pathname from FILE pointer.
* *
* Currently only supported on Linux, otherwise will report "(unknown)". * Currently only supported on Linux, otherwise will report "(unknown)".
* *
@ -228,12 +225,12 @@ namespace LAMMPS_NS {
*/ */
const char *guesspath(char *buf, int len, FILE *fp); const char *guesspath(char *buf, int len, FILE *fp);
/** /** Strip off leading part of path, return just the filename
* \brief Strip off leading part of path, return just the filename *
* \param path file path * \param path file path
* \return file name * \return file name
*/ */
std::string path_basename(const std::string & path); std::string path_basename(const std::string &path);
/** /**
* \brief Join two paths * \brief Join two paths
@ -241,64 +238,65 @@ namespace LAMMPS_NS {
* \param b second path * \param b second path
* \return combined path * \return combined path
*/ */
std::string path_join(const std::string & a, const std::string & b); std::string path_join(const std::string &a, const std::string &b);
/** /**
* \brief Check if file exists and is readable * \brief Check if file exists and is readable
* \param path file path * \param path file path
* \return true if file exists and is readable * \return true if file exists and is readable
*/ */
bool file_is_readable(const std::string & path); bool file_is_readable(const std::string &path);
/** /** Determine full path of potential file. If file is not found in current directory,
* \brief Determine full path of potential file * search LAMMPS_POTENTIALS folder
* If file is not found in current directory, search LAMMPS_POTENTIALS folder *
* \param path file path * \param path file path
* \return full path to potential file * \return full path to potential file
*/ */
std::string get_potential_file_path(const std::string& path); std::string get_potential_file_path(const std::string &path);
/** /** Read potential file and return DATE field if it is present
* \brief Read potential file and return DATE field if it is present *
* \param path file path * \param path file path
* \param potential_name name of potential that is being read * \param potential_name name of potential that is being read
* \return DATE field if present * \return DATE field if present
*/ */
std::string get_potential_date(const std::string & path, const std::string & potential_name); std::string get_potential_date(const std::string &path, const std::string &potential_name);
/** /** Read potential file and return UNITS field if it is present
* \brief Read potential file and return UNITS field if it is present *
* \param path file path * \param path file path
* \param potential_name name of potential that is being read * \param potential_name name of potential that is being read
* \return UNITS field if present * \return UNITS field if present
*/ */
std::string get_potential_units(const std::string & path, const std::string & potential_name); std::string get_potential_units(const std::string &path, const std::string &potential_name);
enum { NOCONVERT = 0, METAL2REAL = 1, REAL2METAL = 1<<1 }; enum { NOCONVERT = 0, METAL2REAL = 1, REAL2METAL = 1<<1 };
enum { UNKNOWN = 0, ENERGY }; enum { UNKNOWN = 0, ENERGY };
/** /** Return bitmask of available conversion factors for a given property
* \brief Return bitmask of available conversion factors for a given propert *
* \param property property to be converted * \param property property to be converted
* \return bitmask indicating available conversions * \return bitmask indicating available conversions
*/ */
int get_supported_conversions(const int property); int get_supported_conversions(const int property);
/** /** Return unit conversion factor for given property and selected from/to units
* \brief Return unit conversion factor for given property and selected from/to units *
* \param property property to be converted * \param property property to be converted
* \param conversion constant indicating the conversion * \param conversion constant indicating the conversion
* \return conversion factor * \return conversion factor
*/ */
double get_conversion_factor(const int property, const int conversion); double get_conversion_factor(const int property, const int conversion);
/** /** Convert a time string to seconds
* \brief Convert a time string to seconds *
* The strings "off" and "unlimited" result in -1 * The strings "off" and "unlimited" result in -1
*
* \param timespec a string in the following format: ([[HH:]MM:]SS) * \param timespec a string in the following format: ([[HH:]MM:]SS)
* \return total in seconds * \return total in seconds
*/ */
double timespec2seconds(const std::string & timespec); double timespec2seconds(const std::string &timespec);
} }
} }

2
unittest/force-styles/tests/atomic-pair-reax_c.yaml
View File

@ -1,7 +1,7 @@
--- ---
lammps_version: 21 Aug 2020 lammps_version: 21 Aug 2020
date_generated: Sun Aug 23 06:35:54 202 date_generated: Sun Aug 23 06:35:54 202
epsilon: 5e-12 epsilon: 7.5e-12
prerequisites: ! | prerequisites: ! |
pair reax/c pair reax/c
fix qeq/reax fix qeq/reax

2
unittest/force-styles/tests/atomic-pair-reax_c_lgvdw.yaml
View File

@ -1,7 +1,7 @@
--- ---
lammps_version: 21 Aug 2020 lammps_version: 21 Aug 2020
date_generated: Sun Aug 23 06:41:04 202 date_generated: Sun Aug 23 06:41:04 202
epsilon: 5e-11 epsilon: 1e-10
prerequisites: ! | prerequisites: ! |
pair reax/c pair reax/c
fix qeq/reax fix qeq/reax