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

More doc files, misc clean ups

Browse Source
This commit is contained in:
jtclemm
2024-04-10 09:47:55 -06:00
parent c63c1856ec
commit 5383bd2613
20 changed files with 187 additions and 964 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
doc/src/Howto_rheo.rst
View File

@ -90,10 +90,10 @@ criteria for creating/deleting a bond or altering force calculations).
---------- ----------
.. _howto-howto_rheo_palermo: .. _howto_rheo_palermo:
**(Palermo)** Palermo, Clemmer, Wolf, O'Connor, in preparation. **(Palermo)** Palermo, Clemmer, Wolf, O'Connor, in preparation.
.. _howto-howto_rheo_clemmer: .. _howto_rheo_clemmer:
**(Clemmer)** Clemmer, Pierce, O'Connor, Nevins, Jones, Lechman, Tencer, Appl. Math. Model., 130, 310-326 (2024). **(Clemmer)** Clemmer, Pierce, O'Connor, Nevins, Jones, Lechman, Tencer, Appl. Math. Model., 130, 310-326 (2024).

119
doc/src/bond_rheo_shell.rst
View File

@ -10,10 +10,13 @@ Syntax
bond_style rheo/shell keyword value attribute1 attribute2 ... bond_style rheo/shell keyword value attribute1 attribute2 ...
* optional keyword = *overlay/pair* or *store/local* or *smooth* or *break* * required keyword = *t/form*
* optional keyword = *store/local*
.. parsed-literal:: .. parsed-literal::
*t/form* value = formation time for a bond (time units)
*store/local* values = fix_ID N attributes ... *store/local* values = fix_ID N attributes ...
* fix_ID = ID of associated internal fix to store data * fix_ID = ID of associated internal fix to store data
* N = prepare data for output every this many timesteps * N = prepare data for output every this many timesteps
@ -24,56 +27,56 @@ Syntax
*x, y, z* = the center of mass position of the 2 atoms when the bond broke (distance units) *x, y, z* = the center of mass position of the 2 atoms when the bond broke (distance units)
*x/ref, y/ref, z/ref* = the initial center of mass position of the 2 atoms (distance units) *x/ref, y/ref, z/ref* = the initial center of mass position of the 2 atoms (distance units)
*overlay/pair* value = *yes* or *no*
bonded particles will still interact with pair forces
*smooth* value = *yes* or *no*
smooths bond forces near the breaking point
*normalize* value = *yes* or *no*
normalizes bond forces by the reference length
*break* value = *yes* or *no*
indicates whether bonds break during a run
Examples Examples
"""""""" """"""""
.. code-block:: LAMMPS .. code-block:: LAMMPS
bond_style bpm/spring bond_style rheo/shell t/form 10.0
bond_coeff 1 1.0 0.05 0.1 bond_coeff 1 1.0 0.05 0.1
bond_style bpm/spring myfix 1000 time id1 id2
dump 1 all local 1000 dump.broken f_myfix[1] f_myfix[2] f_myfix[3]
dump_modify 1 write_header no
Description Description
""""""""""" """""""""""
.. versionadded:: 4May2022 .. versionadded:: TBD
The *bpm/spring* bond style computes forces based on The *rheo/shell* bond style is designed to work with
deviations from the initial reference state of the two atoms. The :doc:`fix rheo/oxidation <fix_rheo_oxidation>` which creates candidate
reference state is stored by each bond when it is first computed in bonds between eligible surface or near-surface particles. When a bond
is first created, it computes no forces and starts a timer. Forces are
not computed until the timer reaches the specified bond formation time
and the bond is fully enabled. If the two particles move outside of the
maximum bond distance or move into the bulk before the timer reaches
the formation time, the bond automatically deletes itself. Not that
this deletion does not generate any broken bond data saved to a
*store/local* fix.
Before bonds are enabled, they are still treated as regular bonds by
all other parts of LAMMPS. This means they are written to data files
and counted in computes such as :doc:`nbond/atom <compute_nbond_atom>`.
To only count enabled bonds, use the *nbond/shell* attribute in
:doc:`compute property/atom/rheo <compute_property_atom_rheo>`.
When enabled, the bond then computes forces based on deviations from
the initial reference state of the two atoms much like a BPM style
bond (as further discussed in the :doc:`BPM howto page <Howto_bpm>`).
The reference state is stored by each bond when it is first enabled
the setup of a run. Data is then preserved across run commands and is the setup of a run. Data is then preserved across run commands and is
written to :doc:`binary restart files <restart>` such that restarting written to :doc:`binary restart files <restart>` such that restarting
the system will not reset the reference state of a bond. the system will not reset the reference state of a bond or the timer.
This bond style only applies central-body forces which conserve the This bond style is based on a model described in Ref.
translational and rotational degrees of freedom of a bonded set of :ref:`(Clemmer) <howto_rheo_clemmer>`. The force has a magnitude of
particles based on a model described by Clemmer and Robbins
:ref:`(Clemmer) <fragment-Clemmer>`. The force has a magnitude of
.. math:: .. math::
F = k (r - r_0) w F = 2 k (r - r_0) + \frac{2 * k}{r_0^2 \epsilon_c^2} (r - r_0)^3
where :math:`k` is a stiffness, :math:`r` is the current distance where :math:`k` is a stiffness, :math:`r` is the current distance
and :math:`r_0` is the initial distance between the two particles, and and :math:`r_0` is the initial distance between the two particles, and
:math:`w` is an optional smoothing factor discussed below. Bonds will :math:`\epsilon_c` is maximum strain beyond which a bond breaks. This
break at a strain of :math:`\epsilon_c`. This is done by setting is done by setting the bond type to 0 such that forces are no longer
the bond type to 0 such that forces are no longer computed. computed.
An additional damping force is applied to the bonded An additional damping force is applied to the bonded
particles. This forces is proportional to the difference in the particles. This forces is proportional to the difference in the
@ -88,15 +91,6 @@ where :math:`\gamma` is the damping strength, :math:`\hat{r}` is the
radial normal vector, and :math:`\vec{v}` is the velocity difference radial normal vector, and :math:`\vec{v}` is the velocity difference
between the two particles. between the two particles.
The smoothing factor :math:`w` can be added or removed by setting the
*smooth* keyword to *yes* or *no*, respectively. It is constructed such
that forces smoothly go to zero, avoiding discontinuities, as bonds
approach the critical strain
.. math::
w = 1.0 - \left( \frac{r - r_0}{r_0 \epsilon_c} \right)^8 .
The following coefficients must be defined for each bond type via the The following coefficients must be defined for each bond type via the
:doc:`bond_coeff <bond_coeff>` command as in the example above, or in :doc:`bond_coeff <bond_coeff>` command as in the example above, or in
the data file or restart files read by the :doc:`read_data the data file or restart files read by the :doc:`read_data
@ -106,22 +100,11 @@ the data file or restart files read by the :doc:`read_data
* :math:`\epsilon_c` (unit less) * :math:`\epsilon_c` (unit less)
* :math:`\gamma` (force/velocity units) * :math:`\gamma` (force/velocity units)
If the *normalize* keyword is set to *yes*, the elastic bond force will be Unlike other BPM-style bonds, this bond style does not update special
normalized by :math:`r_0` such that :math:`k` must be given in force units. bond settings when bonds are created or deleted. This bond style also
does not enforce specific :doc:`special_bonds <special_bonds>` settings.
By default, pair forces are not calculated between bonded particles. This behavior is purposeful such :doc:`RHEO pair forces <pair_rheo>`
Pair forces can alternatively be overlaid on top of bond forces by setting and heat flows are still calculated.
the *overlay/pair* keyword to *yes*. These settings require specific
:doc:`special_bonds <special_bonds>` settings described in the
restrictions. Further details can be found in the :doc:`how to <Howto_bpm>`
page on BPMs.
.. versionadded:: 28Mar2023
If the *break* keyword is set to *no*, LAMMPS assumes bonds should not break
during a simulation run. This will prevent some unnecessary calculation.
However, if a bond reaches a strain greater than :math:`\epsilon_c`,
it will trigger an error.
If the *store/local* keyword is used, an internal fix will track bonds that If the *store/local* keyword is used, an internal fix will track bonds that
break during the simulation. Whenever a bond breaks, data is processed break during the simulation. Whenever a bond breaks, data is processed
@ -187,39 +170,25 @@ extra quantity can be accessed by the
Restrictions Restrictions
"""""""""""" """"""""""""
This bond style is part of the BPM package. It is only enabled if This bond style is part of the RHEO package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info. <Build_package>` page for more info.
By default if pair interactions between bonded atoms are to be disabled,
this bond style requires setting
.. code-block:: LAMMPS
special_bonds lj 0 1 1 coul 1 1 1
and :doc:`newton <newton>` must be set to bond off. If the *overlay/pair*
keyword is set to *yes*, this bond style alternatively requires setting
.. code-block:: LAMMPS
special_bonds lj/coul 1 1 1
Related commands Related commands
"""""""""""""""" """"""""""""""""
:doc:`bond_coeff <bond_coeff>`, :doc:`pair bpm/spring <pair_bpm_spring>` :doc:`bond_coeff <bond_coeff>`, :doc:`fix rheo/oxidation <fix_rheo_oxidation>`
Default Default
""""""" """""""
The option defaults are *overlay/pair* = *no*, *smooth* = *yes*, *normalize* = *no*, and *break* = *yes* NA
---------- ----------
.. _fragment-Clemmer: .. _howto_rheo_clemmer:
**(Clemmer)** Clemmer and Robbins, Phys. Rev. Lett. (2022). **(Clemmer)** Clemmer, Pierce, O'Connor, Nevins, Jones, Lechman, Tencer, Appl. Math. Model., 130, 310-326 (2024).
.. _Groot4: .. _Groot4:

39
doc/src/compute_rheo_property_atom.rst
View File

@ -27,8 +27,8 @@ Syntax
.. parsed-literal:: .. parsed-literal::
*phase* = atom phase status *phase* = atom phase state
*chi* = atom phase neighborhood metric *chi* = atom local phase metric
*surface* = atom surface status *surface* = atom surface status
*surface/r* = atom distance from the surface *surface/r* = atom distance from the surface
*surface/divr* = divergence of position at atom position *surface/divr* = divergence of position at atom position
@ -45,6 +45,7 @@ Syntax
*status* = atom full status *status* = atom full status
*rho* = atom density *rho* = atom density
*grad/v/\** = atom velocity gradient *grad/v/\** = atom velocity gradient
*nbond/shell* = number of oxide bonds
Examples Examples
"""""""" """"""""
@ -52,21 +53,34 @@ Examples
.. code-block:: LAMMPS .. code-block:: LAMMPS
compute 1 all rheo/property/atom phase surface/r pressure compute 1 all rheo/property/atom phase surface/r pressure
compute 2 all rheo/property/atom shift/v/x grad/v/xx
Description Description
""""""""""" """""""""""
Define a computation that simply stores atom attributes specific to the .. versionadded:: TBD
RHEO package for each atom in the group. This is useful so that the
values can be used by other :doc:`output commands <Howto_output>` that Define a computation that stores atom attributes specific to the RHEO
take computes as inputs. See for example, the :doc:`compute reduce package for each atom in the group. This is useful so that the values
<compute_reduce>`, :doc:`fix ave/atom <fix_ave_atom>`, :doc:`fix can be used by other :doc:`output commands <Howto_output>` that take
ave/histo <fix_ave_histo>`, :doc:`fix ave/chunk <fix_ave_chunk>`, computes as inputs. See for example, the
and :doc:`atom-style variable <variable>` commands. :doc:`compute reduce <compute_reduce>`,
:doc:`fix ave/atom <fix_ave_atom>`,
:doc:`fix ave/histo <fix_ave_histo>`,
:doc:`fix ave/chunk <fix_ave_chunk>`, and
:doc:`atom-style variable <variable>` commands.
The possible attributes are described in more detail in other RHEO doc The possible attributes are described in more detail in other RHEO doc
pages include :doc:`fix rheo <fix_rheo>`, :doc:`pair rheo <pair_rheo>`, pages including :doc:`the RHEO howto page <Howto_rheo>`. Many
and :doc:`the RHEO howto page <Howto_rheo>`. properties require their respective fixes, listed below in related
commands, be defined.
The *surface/n/\** and *shift/v/\** attributes are vectors that require
specification of the *x*, *y*, or *z* component, e.g. *surface/n/x*.
The *grad/v/\** attribute is a tensor and requires specification of
the *xx*, *yy*, *zz*, *xy*, *xz*, *yx*, *yz*, *zx*, or *zy* component,
e.g. *grad/v/xy*.
The values are stored in a per-atom vector or array as discussed The values are stored in a per-atom vector or array as discussed
below. Zeroes are stored for atoms not in the specified group or for below. Zeroes are stored for atoms not in the specified group or for
@ -98,7 +112,8 @@ Related commands
:doc:`fix rheo/viscosity <fix_rheo_viscosity>`, :doc:`fix rheo/viscosity <fix_rheo_viscosity>`,
:doc:`fix rheo/pressure <fix_rheo_pressure>`, :doc:`fix rheo/pressure <fix_rheo_pressure>`,
:doc:`fix rheo/thermal <fix_rheo_thermal>`, :doc:`fix rheo/thermal <fix_rheo_thermal>`,
:doc:`pair rheo <pair_rheo>` :doc:`fix rheo/oxdiation <fix_rheo_oxidation>`,
:doc:`fix rheo <fix_rheo>`
Default Default
""""""" """""""

2
doc/src/fix_rheo.rst
View File

@ -43,6 +43,8 @@ Examples
Description Description
""""""""""" """""""""""
.. versionadded:: TBD
Perform time integration for RHEO particles, updating positions, velocities, Perform time integration for RHEO particles, updating positions, velocities,
and densities. For a detailed breakdown of the integration timestep and and densities. For a detailed breakdown of the integration timestep and
numerical details, see :ref:`(Palermo) <howto_rheo_palermo>`. For an numerical details, see :ref:`(Palermo) <howto_rheo_palermo>`. For an

53
doc/src/fix_rheo_oxidation.rst
View File

@ -8,43 +8,42 @@ Syntax
.. parsed-literal:: .. parsed-literal::
fix ID group-ID rheo/oxidation cut btype fix ID group-ID rheo/oxidation cut btype rsurf
* ID, group-ID are documented in :doc:`fix <fix>` command * ID, group-ID are documented in :doc:`fix <fix>` command
* rheo/oxidation = style name of this fix command * rheo/oxidation = style name of this fix command
* cut = maximum bond length (distance units) * cut = maximum bond length (distance units)
* btype = type of bonds created * btype = type of bonds created
* rsurf = distance from surface to create bonds (distance units)
Examples Examples
"""""""" """"""""
.. code-block:: LAMMPS .. code-block:: LAMMPS
fix 1 all rheo/oxidation 1.5 2 fix 1 all rheo/oxidation 1.5 2 0.0
fix 1 all rheo/oxidation 1.0 1 2.0
Description Description
""""""""""" """""""""""
This fix... .. versionadded:: TBD
Each list consists of a series of type This fix dynamically creates bonds on the surface of fluids to
ranges separated by commas. The range can be specified as a represent physical processes such as oxidation. It is intended
single numeric value, or a wildcard asterisk can be used to specify a range for use with bond style :doc:`bond rheo/shell <bond_rheo_shell>`.
of values. This takes the form "\*" or "\*n" or "n\*" or "m\*n". For
example, if M = the number of atom types, then an asterisk with no numeric
values means all types from 1 to M. A leading asterisk means all types
from 1 to n (inclusive). A trailing asterisk means all types from n to M
(inclusive). A middle asterisk means all types from m to n (inclusive).
Note that all atom types must be included in exactly one of the N collections.
While the *Tfreeze* keyword is optional, the *conductivity* and Every timestep, particles check neighbors within a distance of *cut*.
*specific/heat* keywords are mandatory. This distance must be smaller than the kernel length defined in
:doc:`fix rheo <fix_rheo>`. If both particles are on the fluid surface,
or within a distance of *rsurf* from the surface, a bond of type
*btype* is created between the two particles. This process is
further described in Ref. :ref:`(Clemmer) <howto_rheo_clemmer>`.
Multiple instances of this fix may be defined to apply different If used in conjunction with solid bodies, such as those generated
properties to different groups. However, the union of fix groups by the *react* option of :doc:`fix rheo/thermal <fix_rheo_thermal>`,
across all instances of fix rheo/thermal must cover all atoms. it is recommended that one uses a :doc:`hybrid bond style <bond_hybrid>`
If there are multiple instances of this fix, any intersections in with different bond types for solid and oxide bonds.
the fix groups will lead to incorrect thermal integration.
Restart, fix_modify, output, run start/stop, minimize info Restart, fix_modify, output, run start/stop, minimize info
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
@ -58,10 +57,8 @@ the :doc:`run <run>` command. This fix is not invoked during :doc:`energy minim
Restrictions Restrictions
"""""""""""" """"""""""""
This fix must be used with an atom style that includes temperature, This fix must be used with an bond style :doc:`rheo/shell <bond_rheo_shell>`
heatflow, and conductivity such as atom_tyle rheo/thermal This fix and :doc:`fix rheo <fix_rheo>` with surface detection enabled.
must be used in conjuction with :doc:`fix rheo <fix_rheo>` with the
*thermal* setting.
This fix is part of the RHEO package. It is only enabled if This fix is part of the RHEO package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package <Build_package>` page for more info. LAMMPS was built with that package. See the :doc:`Build package <Build_package>` page for more info.
@ -70,12 +67,16 @@ Related commands
"""""""""""""""" """"""""""""""""
:doc:`fix rheo <fix_rheo>`, :doc:`fix rheo <fix_rheo>`,
:doc:`fix rheo/viscosity <fix_rheo_viscosity>`, :doc:`bond rheo/shell <bond_rheo_shell>`,
:doc:`fix rheo/pressure <fix_rheo_pressure>`,
:doc:`pair rheo <pair_rheo>`,
:doc:`compute rheo/property/atom <compute_rheo_property_atom>` :doc:`compute rheo/property/atom <compute_rheo_property_atom>`
Default Default
""""""" """""""
none none
----------
.. _howto_rheo_clemmer:
**(Clemmer)** Clemmer, Pierce, O'Connor, Nevins, Jones, Lechman, Tencer, Appl. Math. Model., 130, 310-326 (2024).

2
doc/src/fix_rheo_pressure.rst
View File

@ -33,6 +33,8 @@ Examples
Description Description
""""""""""" """""""""""
.. versionadded:: TBD
This fix defines a pressure equation of state for RHEO particles. One can This fix defines a pressure equation of state for RHEO particles. One can
define different equations of state for different atom types, but an define different equations of state for different atom types, but an
equation must be specified for every atom type. equation must be specified for every atom type.

2
doc/src/fix_rheo_thermal.rst
View File

@ -48,6 +48,8 @@ Examples
Description Description
""""""""""" """""""""""
.. versionadded:: TBD
This fix performs time integration of temperature evolution for atom style This fix performs time integration of temperature evolution for atom style
rheo/thermal. In addition, it defines multiple thermal properties of rheo/thermal. In addition, it defines multiple thermal properties of
particles and handles melting/solidification, if applicable. For more details particles and handles melting/solidification, if applicable. For more details

2
doc/src/fix_rheo_viscosity.rst
View File

@ -32,6 +32,8 @@ Examples
Description Description
""""""""""" """""""""""
.. versionadded:: TBD
This fix defines a viscosity for RHEO particles. One can define different This fix defines a viscosity for RHEO particles. One can define different
viscosities for different atom types, but a viscosity must be specified for viscosities for different atom types, but a viscosity must be specified for
every atom type. every atom type.

2
doc/src/pair_rheo.rst
View File

@ -31,6 +31,8 @@ Examples
Description Description
""""""""""" """""""""""
.. versionadded:: TBD
Pair style *rheo* computes pressure and viscous forces between particles Pair style *rheo* computes pressure and viscous forces between particles
in the :doc:`rheo package <Howto_rheo>`. If thermal evolution is turned in the :doc:`rheo package <Howto_rheo>`. If thermal evolution is turned
on in :doc:`fix rheo <fix_rheo>`, then the pair style also calculates on in :doc:`fix rheo <fix_rheo>`, then the pair style also calculates

67
doc/src/pair_rheo_react.rst
View File

@ -1,67 +0,0 @@
.. index:: pair_style rheo/react
pair_style rheo/react command
=========================
Syntax
""""""
.. code-block:: LAMMPS
pair_style rheo/react
Examples
""""""""
.. code-block:: LAMMPS
pair_style rheo/react
pair_coeff * * 1.0 1.5 1.0 0.05 1.0 100 2.0
Description
"""""""""""
pair style...
The following coefficients must be defined for each pair of atom types
via the :doc:`pair_coeff <pair_coeff>` command as in the example above,
or in the data file or restart files read by the
:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
commands, or by mixing as described below:
* :math:`k` (force/distance units)
* :math:`r_max` (distance units)
* :math:`\epsilon` (unitless)
* :math:`\gamma` (force/velocity units)
* :math:`t_form` (time units)
* :math:`r_from_surface` (distance units)
----------
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This style does not support the :doc:`pair_modify <pair_modify>`
shift, table, and tail options.
This style does not write information to :doc:`binary restart files <restart>`. Thus, you need to re-specify the pair_style and
pair_coeff commands in an input script that reads a restart file.
This style can only be used via the *pair* keyword of the :doc:`run_style respa <run_style>` command. It does not support the *inner*, *middle*, *outer* keywords.
Restrictions
""""""""""""
This fix is part of the RHEO package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package <Build_package>` page for more info.
Related commands
""""""""""""""""
:doc:`fix rheo <fix_rheo>`,
:doc:`compute rheo/property/atom <compute_rheo_property_atom>`
Default
"""""""
none

2
doc/src/pair_rheo_solid.rst
View File

@ -21,6 +21,8 @@ Examples
Description Description
""""""""""" """""""""""
.. versionadded:: TBD
Style *rheo/solid* is effectively a copy of pair style Style *rheo/solid* is effectively a copy of pair style
:doc:`bpm/spring <pair_bpm_spring>` except it only applies forces :doc:`bpm/spring <pair_bpm_spring>` except it only applies forces
between solid RHEO particles, determined by checking the status of between solid RHEO particles, determined by checking the status of

54
src/RHEO/compute_rheo_property_atom.cpp
View File

@ -57,7 +57,8 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
if (nvalues == 1) size_peratom_cols = 0; if (nvalues == 1) size_peratom_cols = 0;
else size_peratom_cols = nvalues; else size_peratom_cols = nvalues;
pressure_flag = thermal_flag = interface_flag = surface_flag = shift_flag = shell_flag = 0; pressure_flag = thermal_flag = interface_flag = 0;
surface_flag = shift_flag = shell_flag = 0;
// parse input values // parse input values
// customize a new keyword by adding to if statement // customize a new keyword by adding to if statement
@ -70,32 +71,34 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
for (int iarg = 3; iarg < narg; iarg++) { for (int iarg = 3; iarg < narg; iarg++) {
i = iarg-3; i = iarg-3;
if (strcmp(arg[iarg],"phase") == 0) { if (strcmp(arg[iarg], "phase") == 0) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_phase; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_phase;
} else if (strcmp(arg[iarg],"rho") == 0) { } else if (strcmp(arg[iarg], "rho") == 0) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_rho; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_rho;
} else if (strcmp(arg[iarg],"chi") == 0) { } else if (strcmp(arg[iarg], "chi") == 0) {
interface_flag = 1; interface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_chi; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_chi;
} else if (strcmp(arg[iarg],"surface") == 0) { } else if (strcmp(arg[iarg], "surface") == 0) {
surface_flag = 1; surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface;
} else if (strcmp(arg[iarg],"surface/r") == 0) { } else if (strcmp(arg[iarg], "surface/r") == 0) {
surface_flag = 1; surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_r; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_r;
} else if (strcmp(arg[iarg],"surface/divr") == 0) { } else if (strcmp(arg[iarg], "surface/divr") == 0) {
surface_flag = 1; surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_divr; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_divr;
} else if (utils::strmatch(arg[iarg], "^surface/n")) { } else if (utils::strmatch(arg[iarg], "^surface/n")) {
surface_flag = 1; surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_n; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_n;
col_index[i] = get_vector_index(arg[iarg]); col_index[i] = get_vector_index(arg[iarg]);
} else if (strcmp(arg[iarg],"coordination") == 0) { } else if (strcmp(arg[iarg], "coordination") == 0) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_coordination; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_coordination;
} else if (strcmp(arg[iarg],"pressure") == 0) { } else if (strcmp(arg[iarg], "pressure") == 0) {
pressure_flag = 1; pressure_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_pressure; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_pressure;
} else if (strcmp(arg[iarg],"cv") == 0) { } else if (strcmp(arg[iarg], "viscosity") == 0) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_viscosity;
} else if (strcmp(arg[iarg], "cv") == 0) {
thermal_flag = 1; thermal_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_cv; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_cv;
} else if (utils::strmatch(arg[iarg], "^shift/v")) { } else if (utils::strmatch(arg[iarg], "^shift/v")) {
@ -113,7 +116,7 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
arg[iarg], atom->get_style()); arg[iarg], atom->get_style());
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_atom_style; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_atom_style;
thermal_flag = 1; thermal_flag = 1;
} else if (strcmp(arg[iarg],"nbond/shell") == 0) { } else if (strcmp(arg[iarg], "nbond/shell") == 0) {
shell_flag = 1; shell_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_nbond_shell; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_nbond_shell;
} else { } else {
@ -124,9 +127,9 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
arg[iarg], atom->get_style()); arg[iarg], atom->get_style());
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_atom_style; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_atom_style;
if (strcmp(arg[iarg],"temperature") == 0) thermal_flag = 1; if (strcmp(arg[iarg], "temperature") == 0) thermal_flag = 1;
if (strcmp(arg[iarg],"heatflow") == 0) thermal_flag = 1; if (strcmp(arg[iarg], "heatflow") == 0) thermal_flag = 1;
if (strcmp(arg[iarg],"conductivity") == 0) thermal_flag = 1; if (strcmp(arg[iarg], "conductivity") == 0) thermal_flag = 1;
} }
} }
@ -149,7 +152,7 @@ ComputeRHEOPropertyAtom::~ComputeRHEOPropertyAtom()
void ComputeRHEOPropertyAtom::init() void ComputeRHEOPropertyAtom::init()
{ {
auto fixes = modify->get_fix_by_style("^rheo$"); auto fixes = modify->get_fix_by_style("^rheo$");
if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/viscosity"); if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use compute rheo/property/atom");
fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]); fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
if (interface_flag && !(fix_rheo->interface_flag)) if (interface_flag && !(fix_rheo->interface_flag))
@ -390,6 +393,21 @@ void ComputeRHEOPropertyAtom::pack_pressure(int n)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_viscosity(int n)
{
int *mask = atom->mask;
double *viscosity = atom->viscosity;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = viscosity[i];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_nbond_shell(int n) void ComputeRHEOPropertyAtom::pack_nbond_shell(int n)
{ {
int *nbond = fix_oxidation->nbond; int *nbond = fix_oxidation->nbond;
@ -490,11 +508,11 @@ int ComputeRHEOPropertyAtom::get_tensor_index(char* option)
int ComputeRHEOPropertyAtom::get_vector_index(char* option) int ComputeRHEOPropertyAtom::get_vector_index(char* option)
{ {
int index; int index;
if (utils::strmatch(option,"x$")) { if (utils::strmatch(option, "x$")) {
index = 0; index = 0;
} else if (utils::strmatch(option,"y$")) { } else if (utils::strmatch(option, "y$")) {
index = 1; index = 1;
} else if (utils::strmatch(option,"z$")) { } else if (utils::strmatch(option, "z$")) {
if (domain->dimension == 2) if (domain->dimension == 2)
error->all(FLERR, "Invalid compute rheo/property/atom property {} in 2D", option); error->all(FLERR, "Invalid compute rheo/property/atom property {} in 2D", option);
index = 2; index = 2;

4
src/RHEO/compute_rheo_property_atom.h
View File

@ -35,7 +35,8 @@ class ComputeRHEOPropertyAtom : public Compute {
private: private:
int nvalues, nmax; int nvalues, nmax;
int pressure_flag, thermal_flag, interface_flag, surface_flag, shift_flag, shell_flag; int pressure_flag, thermal_flag, interface_flag;
int surface_flag, shift_flag, shell_flag;
int *avec_index; int *avec_index;
int *col_index; int *col_index;
double *buf; double *buf;
@ -55,6 +56,7 @@ class ComputeRHEOPropertyAtom : public Compute {
void pack_shift_v(int); void pack_shift_v(int);
void pack_gradv(int); void pack_gradv(int);
void pack_pressure(int); void pack_pressure(int);
void pack_viscosity(int);
void pack_nbond_shell(int); void pack_nbond_shell(int);
void pack_atom_style(int); void pack_atom_style(int);

6
src/RHEO/fix_rheo.cpp
View File

@ -19,6 +19,7 @@
#include "fix_rheo.h" #include "fix_rheo.h"
#include "atom.h" #include "atom.h"
#include "citeme.h"
#include "compute_rheo_grad.h" #include "compute_rheo_grad.h"
#include "compute_rheo_interface.h" #include "compute_rheo_interface.h"
#include "compute_rheo_surface.h" #include "compute_rheo_surface.h"
@ -37,6 +38,9 @@ using namespace LAMMPS_NS;
using namespace RHEO_NS; using namespace RHEO_NS;
using namespace FixConst; using namespace FixConst;
static const char cite_rheo[] =
"TBD\n\n";
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
@ -137,6 +141,8 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
} }
iarg += 1; iarg += 1;
} }
if (lmp->citeme) lmp->citeme->add(cite_rheo);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */

15
src/RHEO/fix_rheo_oxidation.cpp
View File

@ -20,6 +20,7 @@
#include "atom.h" #include "atom.h"
#include "atom_vec.h" #include "atom_vec.h"
#include "citeme.h"
#include "compute_rheo_surface.h" #include "compute_rheo_surface.h"
#include "error.h" #include "error.h"
#include "fix_rheo.h" #include "fix_rheo.h"
@ -34,6 +35,18 @@ using namespace RHEO_NS;
using namespace FixConst; using namespace FixConst;
enum {NONE, CONSTANT}; enum {NONE, CONSTANT};
static const char cite_rheo_oxide[] =
"@article{ApplMathModel.130.310,\n"
" title = {A hybrid smoothed-particle hydrodynamics model of oxide skins on molten aluminum},\n"
" journal = {Applied Mathematical Modelling},\n"
" volume = {130},\n"
" pages = {310-326},\n"
" year = {2024},\n"
" issn = {0307-904X},\n"
" doi = {https://doi.org/10.1016/j.apm.2024.02.027},\n"
" author = {Joel T. Clemmer and Flint Pierce and Thomas C. O'Connor and Thomas D. Nevins and Elizabeth M.C. Jones and Jeremy B. Lechman and John Tencer},\n"
"}\n\n";
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixRHEOOxidation::FixRHEOOxidation(LAMMPS *lmp, int narg, char **arg) : FixRHEOOxidation::FixRHEOOxidation(LAMMPS *lmp, int narg, char **arg) :
@ -51,6 +64,8 @@ FixRHEOOxidation::FixRHEOOxidation(LAMMPS *lmp, int narg, char **arg) :
if (rsurf <= 0.0) error->all(FLERR, "Illegal surface distance {} in fix rheo/oxidation", cut); if (rsurf <= 0.0) error->all(FLERR, "Illegal surface distance {} in fix rheo/oxidation", cut);
cutsq = cut * cut; cutsq = cut * cut;
if (lmp->citeme) lmp->citeme->add(cite_rheo_oxide);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */

137
src/RHEO/fix_rheo_stress.cpp
View File

@ -1,137 +0,0 @@
/* ----------------------------------------------------------------------
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 authors: Joel Clemmer (SNL)
----------------------------------------------------------------------- */
#include "fix_rheo_stress.h"
#include "atom.h"
#include "comm.h"
#include "compute.h"
#include "domain.h"
#include "fix_store_atom.h"
#include "group.h"
#include "error.h"
#include "modify.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixRHEOStress::FixRHEOStress(LAMMPS *lmp, int narg, char **arg) :
id_compute(nullptr), id_fix(nullptr), stress_compute(nullptr), store_fix(nullptr), Fix(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR,"Illegal fix rheo/stress command");
comm_forward = 6;
}
/* ---------------------------------------------------------------------- */
FixRHEOStress::~FixRHEOStress()
{
modify->delete_compute(id_compute);
modify->delete_fix(id_fix);
}
/* ---------------------------------------------------------------------- */
void FixRHEOStress::post_constructor()
{
id_fix = utils::strdup(std::string(id) + "_store");
store_fix = dynamic_cast<FixStoreAtom *>(modify->add_fix(fmt::format("{} {} STORE/ATOM d_pxx d_pyy d_pzz d_pxy d_pxz d_pyz", id_fix, group->names[igroup])));
array_atom = store_fix->astore;
id_compute = utils::strdup(std::string(id) + "_compute");
stress_compute = modify->add_compute(fmt::format("{} {} stress/atom NULL ke pair bond", id_compute, group->names[igroup]));
}
/* ---------------------------------------------------------------------- */
int FixRHEOStress::setmask()
{
int mask = 0;
mask |= PRE_FORCE;
mask |= END_OF_STEP;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixRHEOStress::init()
{
stress_compute->addstep(update->ntimestep+1);
}
/* ---------------------------------------------------------------------- */
void FixRHEOStress::pre_force(int vflag)
{
// add pre-force and forward to ghosts (not done in store/atom)
comm->forward_comm(this);
}
/* ---------------------------------------------------------------------- */
void FixRHEOStress::end_of_step()
{
stress_compute->compute_peratom();
// copy compute to fix property atom
double **saved_stress = store_fix->astore;
double **stress = stress_compute->array_atom;
int ntotal = atom->nlocal+atom->nghost;
for (int i = 0; i < ntotal; i++)
for (int a = 0; a < 6; a++)
saved_stress[i][a] = stress[i][a];
stress_compute->addstep(update->ntimestep + 1);
}
/* ---------------------------------------------------------------------- */
int FixRHEOStress::pack_forward_comm(int n, int *list, double *buf,
int /*pbc_flag*/, int * /*pbc*/)
{
int i, j, a, m;
double **saved_stress = store_fix->astore;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
for (a = 0; a < 6; a++)
buf[m++] = saved_stress[j][a];
}
return m;
}
/* ---------------------------------------------------------------------- */
void FixRHEOStress::unpack_forward_comm(int n, int first, double *buf)
{
int i, a, m, last;
double **saved_stress = store_fix->astore;
m = 0;
last = first + n;
for (i = first; i < last; i++)
for (a = 0; a < 6; a++)
saved_stress[i][a] = buf[m++];
}

48
src/RHEO/fix_rheo_stress.h
View File

@ -1,48 +0,0 @@
/* -*- 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(rheo/stress,FixRHEOStress);
// clang-format on
#else
#ifndef LMP_FIX_RHEO_STRESS_H
#define LMP_FIX_RHEO_STRESS_H
#include "fix.h"
namespace LAMMPS_NS {
class FixRHEOStress : public Fix {
public:
FixRHEOStress(class LAMMPS *, int, char **);
~FixRHEOStress() override;
void post_constructor() override;
int setmask() override;
void init() override;
void pre_force(int) override;
void end_of_step() override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
private:
char *id_compute, *id_fix;
class Compute *stress_compute;
class FixStoreAtom *store_fix;
};
} // namespace LAMMPS_NS
#endif
#endif

15
src/RHEO/fix_rheo_thermal.cpp
View File

@ -20,6 +20,7 @@
#include "atom.h" #include "atom.h"
#include "atom_vec.h" #include "atom_vec.h"
#include "citeme.h"
#include "comm.h" #include "comm.h"
#include "compute_rheo_grad.h" #include "compute_rheo_grad.h"
#include "compute_rheo_vshift.h" #include "compute_rheo_vshift.h"
@ -43,6 +44,18 @@ using namespace RHEO_NS;
using namespace FixConst; using namespace FixConst;
enum {NONE, CONSTANT}; enum {NONE, CONSTANT};
static const char cite_rheo_oxide[] =
"@article{ApplMathModel.130.310,\n"
" title = {A hybrid smoothed-particle hydrodynamics model of oxide skins on molten aluminum},\n"
" journal = {Applied Mathematical Modelling},\n"
" volume = {130},\n"
" pages = {310-326},\n"
" year = {2024},\n"
" issn = {0307-904X},\n"
" doi = {https://doi.org/10.1016/j.apm.2024.02.027},\n"
" author = {Joel T. Clemmer and Flint Pierce and Thomas C. O'Connor and Thomas D. Nevins and Elizabeth M.C. Jones and Jeremy B. Lechman and John Tencer},\n"
"}\n\n";
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
@ -193,6 +206,8 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
if (Tc_style[i] == NONE && L_style[i] != NONE) if (Tc_style[i] == NONE && L_style[i] != NONE)
error->all(FLERR, "Must specify critical temperature for atom type {} to use latent heat in fix rheo/thermal", i); error->all(FLERR, "Must specify critical temperature for atom type {} to use latent heat in fix rheo/thermal", i);
} }
if (lmp->citeme) lmp->citeme->add(cite_rheo_oxide);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */

515
src/RHEO/pair_rheo_react.cpp
View File

@ -1,515 +0,0 @@
/* ----------------------------------------------------------------------
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 authors:
Joel Clemmer (SNL)
----------------------------------------------------------------------- */
#include "pair_rheo_react.h"
#include "atom.h"
#include "comm.h"
#include "compute_rheo_surface.h"
#include "error.h"
#include "fix.h"
#include "fix_dummy.h"
#include "fix_neigh_history.h"
#include "fix_rheo.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "utils.h"
using namespace LAMMPS_NS;
using namespace RHEO_NS;
/* ---------------------------------------------------------------------- */
PairRHEOReact::PairRHEOReact(LAMMPS *lmp) : Pair(lmp),
dbond(NULL)
{
single_enable = 0;
size_history = 2;
beyond_contact = 1;
comm_reverse = 1;
nondefault_history_transfer = 1;
// create dummy fix as placeholder for FixNeighHistory
// this is so final order of Modify:fix will conform to input script
fix_history = nullptr;
fix_dummy = dynamic_cast<FixDummy *>(
modify->add_fix("NEIGH_HISTORY_RHEO_REACT_DUMMY" + std::to_string(instance_me) + " all DUMMY"));
// For nbond, create an instance of fix property atom
// Need restarts + exchanging with neighbors since it needs to persist
// between timesteps (fix property atom will handle callbacks)
int tmp1, tmp2;
index_nb = atom->find_custom("react_nbond", tmp1, tmp2);
if (index_nb == -1) {
id_fix = utils::strdup("pair_rheo_react_fix_property_atom");
modify->add_fix(fmt::format("{} all property/atom i_react_nbond", id_fix));
index_nb = atom->find_custom("react_nbond", tmp1, tmp2);
}
nbond = atom->ivector[index_nb];
//Store non-persistent per atom quantities, intermediate
nmax_store = atom->nmax;
memory->create(dbond, nmax_store, "rheo/react:dbond");
}
/* ---------------------------------------------------------------------- */
PairRHEOReact::~PairRHEOReact()
{
if (modify->nfix && fix_history) modify->delete_fix("NEIGH_HISTORY_RHEO_REACT" + std::to_string(instance_me));
if (modify->nfix && fix_dummy) modify->delete_fix("NEIGH_HISTORY_RHEO_REACT_DUMMY" + std::to_string(instance_me));
if (modify->nfix) modify->delete_fix(id_fix);
delete[] id_fix;
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cutbsq);
memory->destroy(cutbond);
memory->destroy(k);
memory->destroy(eps);
memory->destroy(gamma);
memory->destroy(t_form);
memory->destroy(rlimit);
}
memory->destroy(dbond);
}
/* ---------------------------------------------------------------------- */
void PairRHEOReact::compute(int eflag, int vflag)
{
int i, j, ii, jj, inum, jnum, fluidi, fluidj;
double xtmp, ytmp, ztmp, delx, dely, delz;
double vxtmp, vytmp, vztmp, delvx, delvy, delvz;
double rsq, r, rinv, r0, fpair, dot, smooth;
int itype, jtype;
int *ilist, *jlist, *numneigh, **firstneigh;
int *saved, **firstsaved;
double *data, *alldata, **firstdata;
ev_init(eflag,vflag);
int bondupdate = 1;
if (update->setupflag) bondupdate = 0;
double dt = update->dt;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
int *type = atom->type;
int *status = atom->status;
int *mask = atom->mask;
int *nbond = atom->ivector[index_nb];
double *rsurf = compute_surface->rsurface;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
firstsaved = fix_history->firstflag;
firstdata = fix_history->firstvalue;
if (atom->nmax > nmax_store){
nmax_store = atom->nmax;
memory->destroy(dbond);
memory->create(dbond, nmax_store, "rheo/react:dbond");
}
size_t nbytes = nmax_store * sizeof(int);
memset(&dbond[0], 0, nbytes);
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
vxtmp = v[i][0];
vytmp = v[i][1];
vztmp = v[i][2];
fluidi = !(status[i] & PHASECHECK);
saved = firstsaved[i];
alldata = firstdata[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
fluidj = !(status[j] & PHASECHECK);
data = &alldata[2*jj];
// If not bonded and there's an internal fluid particle, unsave any data and skip
if (!(saved[jj] == 1 && data[0] > 0)) {
if ((fluidi && (rsurf[i] > rlimit[itype][jtype])) || (fluidj && (rsurf[j] > rlimit[itype][jtype]))) {
saved[jj] = 0;
continue;
}
}
// If both are solid, unbond and skip
if (!fluidi && !fluidj) {
//If bonded, deincrement
if (saved[jj] == 1 && data[0] > 0) {
dbond[i] --;
dbond[j] --;
}
saved[jj] = 0;
continue;
}
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
// If unbonded and beyond bond distance, unsave and skip
if (data[0] == -1 && rsq > cutbsq[itype][jtype]) {
saved[jj] = 0;
continue;
}
r = sqrt(rsq);
// Initialize data if not currently saved since all could bond if they are on the surface
if (saved[jj] == 0) {
data[0] = -1;
data[1] = 0;
saved[jj] = 1;
}
// Check for bond formation (unbonded) or breakage (bonded)
if (data[0] == -1) {
// If unbonded, check if we count down to bonding if both on surface (not given for r or s)
if (bondupdate && (rsurf[i] <= rlimit[itype][jtype]) && (rsurf[j] <= rlimit[itype][jtype])) {
data[1] += dt;
if (data[1] >= t_form[itype][jtype]) {
data[0] = r;
dbond[i] ++;
dbond[j] ++;
data[1] = 0;
}
}
} else {
// If bonded, check if breaks in tension
r0 = data[0];
if (r > ((1.0 + eps[itype][jtype]) * r0)) {
saved[jj] = 0;
dbond[i] --;
dbond[j] --;
data[0] = -1;
}
}
// Skip if unbonded
if (data[0] <= 0) continue;
delvx = vxtmp - v[j][0];
delvy = vytmp - v[j][1];
delvz = vztmp - v[j][2];
rinv = 1.0 / r;
r0 = data[0];
fpair = k[itype][jtype] * (r0 - r);
dot = delx * delvx + dely * delvy + delz * delvz;
fpair -= gamma[itype][jtype] * dot * rinv;
smooth = 1.0;
if (r > r0) {
smooth = (r - r0) / (r0 * eps[itype][jtype]);
smooth *= smooth;
smooth *= smooth;
smooth = 1 - smooth;
}
fpair *= rinv * smooth;
f[i][0] += delx * fpair;
f[i][1] += dely * fpair;
f[i][2] += delz * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx * fpair;
f[j][1] -= dely * fpair;
f[j][2] -= delz * fpair;
}
if (evflag) ev_tally(i, j, nlocal, newton_pair, 0.0, 0.0, fpair, delx, dely, delz);
}
}
// Communicate changes in nbond
if (newton_pair) comm->reverse_comm(this);
for(i = 0; i < nlocal; i++) {
fluidi = !(status[i] & PHASECHECK);
nbond[i] += dbond[i];
// If it has bonds it is reactive (no shifting)
// If a reactive particle breaks all bonds, return to fluid
// Keep it non-shifting for this timestep to be safe
if (nbond[i] != 0 && fluidi) status[i] |= STATUS_NO_SHIFT;
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairRHEOReact::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag, n + 1, n + 1, "pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cutbond, n + 1, n + 1, "pair:cutbond");
memory->create(cutbsq, n + 1, n + 1, "pair:cutbsq");
memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
memory->create(k, n + 1, n + 1, "pair:k");
memory->create(eps, n + 1, n + 1, "pair:eps");
memory->create(gamma, n + 1, n + 1, "pair:gamma");
memory->create(t_form, n + 1, n + 1, "pair:t_form");
memory->create(rlimit, n + 1, n + 1, "pair:rlimit");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairRHEOReact::settings(int narg, char **arg)
{
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairRHEOReact::coeff(int narg, char **arg)
{
if (narg != 9) error->all(FLERR, "Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo, ihi, jlo, jhi;
utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi,error);
utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi,error);
double k_one = utils::numeric(FLERR, arg[2], false, lmp);
double cutb_one = utils::numeric(FLERR, arg[3], false, lmp);
double eps_one = utils::numeric(FLERR, arg[4], false, lmp);
double gamma_one = utils::numeric(FLERR, arg[5], false, lmp);
double t_form_one = utils::numeric(FLERR, arg[6], false, lmp);
double rlimit_one = utils::numeric(FLERR, arg[7], false, lmp);
if (k_one < 0.0 || eps_one < 0.0 || t_form_one < 0.0)
error->all(FLERR, "Illegal pair_style command");
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
k[i][j] = k_one;
cutbond[i][j] = cutb_one;
eps[i][j] = eps_one;
gamma[i][j] = gamma_one;
t_form[i][j] = t_form_one;
rlimit[i][j] = rlimit_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairRHEOReact::init_style()
{
int irequest = neighbor->request(this, instance_me);
//neighbor->requests[irequest]->history = 1;
if (fix_history == nullptr) {
auto cmd = fmt::format("NEIGH_HISTORY_RHEO_REACT{} all NEIGH_HISTORY {}", instance_me, size_history);
fix_history = dynamic_cast<FixNeighHistory *>(
modify->replace_fix("NEIGH_HISTORY_RHEO_REACT_DUMMY" + std::to_string(instance_me), cmd, 1));
fix_history->pair = this;
fix_dummy = nullptr;
}
}
/* ----------------------------------------------------------------------
setup specific to this pair style
------------------------------------------------------------------------- */
void PairRHEOReact::setup()
{
auto fixes = modify->get_fix_by_style("^rheo$");
if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/tension");
fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
if (!fix_rheo->surface_flag) error->all(FLERR,
"Pair rheo/react requires surface calculation in fix rheo");
compute_surface = fix_rheo->compute_surface;
if (force->newton_pair == 0) error->all(FLERR,
"Pair rheo/react needs newton pair on for bond changes to be consistent");
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairRHEOReact::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR, "All pair coeffs are not set");
cutbsq[i][j] = cutbond[i][j] * cutbond[i][j];
cutbsq[j][i] = cutbsq[i][j];
cutbond[j][i] = cutbond[i][j];
k[j][i] = k[i][j];
eps[j][i] = eps[i][j];
gamma[j][i] = gamma[i][j];
t_form[j][i] = t_form[i][j];
rlimit[j][i] = rlimit[i][j];
double cut = cutbond[i][j] * (1.0 + eps[i][j]);
return cut;
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairRHEOReact::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++)
fwrite(&setflag[i][j], sizeof(int), 1, fp);
if (setflag[i][j]) {
fwrite(&k[i][j], sizeof(double), 1, fp);
fwrite(&cutbond[i][j], sizeof(double), 1, fp);
fwrite(&eps[i][j], sizeof(double), 1, fp);
fwrite(&gamma[i][j], sizeof(double), 1, fp);
fwrite(&t_form[i][j], sizeof(double), 1, fp);
fwrite(&rlimit[i][j], sizeof(double), 1, fp);
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairRHEOReact::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
if (setflag[i][j]) {
if (me == 0) {
utils::sfread(FLERR, &k[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &cutbond[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &eps[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &gamma[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &t_form[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &rlimit[i][j], sizeof(double), 1, fp, nullptr, error);
}
MPI_Bcast(&k[i][j], 1,MPI_DOUBLE, 0, world);
MPI_Bcast(&cutbond[i][j], 1,MPI_DOUBLE, 0, world);
MPI_Bcast(&eps[i][j], 1,MPI_DOUBLE, 0, world);
MPI_Bcast(&gamma[i][j], 1,MPI_DOUBLE, 0, world);
MPI_Bcast(&t_form[i][j], 1,MPI_DOUBLE, 0, world);
MPI_Bcast(&rlimit[i][j], 1,MPI_DOUBLE, 0, world);
}
}
}
/* ----------------------------------------------------------------------
transfer history during fix/neigh/history exchange - transfer same sign
------------------------------------------------------------------------- */
void PairRHEOReact::transfer_history(double* source, double* target)
{
for (int i = 0; i < size_history; i++)
target[i] = source[i];
}
/* ---------------------------------------------------------------------- */
int PairRHEOReact::pack_reverse_comm(int n, int first, double *buf)
{
int i, m, last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
buf[m++] = dbond[i];
}
return m;
}
/* ---------------------------------------------------------------------- */
void PairRHEOReact::unpack_reverse_comm(int n, int *list, double *buf)
{
int i, j, m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
dbond[j] += buf[m++];
}
}

63
src/RHEO/pair_rheo_react.h
View File

@ -1,63 +0,0 @@
/* -*- 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 PAIR_CLASS
// clang-format off
PairStyle(rheo/react,PairRHEOReact)
// clang-format on
#else
#ifndef LMP_PAIR_RHEO_REACT_H
#define LMP_PAIR_RHEO_REACT_H
#include "pair.h"
namespace LAMMPS_NS {
class PairRHEOReact : public Pair {
public:
PairRHEOReact(class LAMMPS *);
~PairRHEOReact() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override;
void setup() override;
double init_one(int, int) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
int pack_reverse_comm(int, int, double *) override;
void unpack_reverse_comm(int, int *, double *) override;
protected:
double **cutbond, **cutbsq, **k, **eps, **gamma, **t_form, **rlimit;
void allocate();
void transfer_history(double*, double*);
int size_history;
int *dbond, *nbond;
int index_nb, nmax_store;
char *id_fix;
class FixDummy *fix_dummy;
class FixNeighHistory *fix_history;
class FixRHEO *fix_rheo;
class ComputeRHEOSurface *compute_surface;
};
} // namespace LAMMPS_NS
#endif
#endif