From 5383bd26139bb0a2584ee96a9f7ff1e9ba1735e5 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Wed, 10 Apr 2024 09:47:55 -0600 Subject: [PATCH] More doc files, misc clean ups --- doc/src/Howto_rheo.rst | 4 +- doc/src/bond_rheo_shell.rst | 119 ++---- doc/src/compute_rheo_property_atom.rst | 39 +- doc/src/fix_rheo.rst | 2 + doc/src/fix_rheo_oxidation.rst | 53 +-- doc/src/fix_rheo_pressure.rst | 2 + doc/src/fix_rheo_thermal.rst | 2 + doc/src/fix_rheo_viscosity.rst | 2 + doc/src/pair_rheo.rst | 2 + doc/src/pair_rheo_react.rst | 67 --- doc/src/pair_rheo_solid.rst | 2 + src/RHEO/compute_rheo_property_atom.cpp | 54 ++- src/RHEO/compute_rheo_property_atom.h | 4 +- src/RHEO/fix_rheo.cpp | 6 + src/RHEO/fix_rheo_oxidation.cpp | 15 + src/RHEO/fix_rheo_stress.cpp | 137 ------- src/RHEO/fix_rheo_stress.h | 48 --- src/RHEO/fix_rheo_thermal.cpp | 15 + src/RHEO/pair_rheo_react.cpp | 515 ------------------------ src/RHEO/pair_rheo_react.h | 63 --- 20 files changed, 187 insertions(+), 964 deletions(-) delete mode 100644 doc/src/pair_rheo_react.rst delete mode 100644 src/RHEO/fix_rheo_stress.cpp delete mode 100644 src/RHEO/fix_rheo_stress.h delete mode 100644 src/RHEO/pair_rheo_react.cpp delete mode 100644 src/RHEO/pair_rheo_react.h diff --git a/doc/src/Howto_rheo.rst b/doc/src/Howto_rheo.rst index 146716ba18..db710d2366 100644 --- a/doc/src/Howto_rheo.rst +++ b/doc/src/Howto_rheo.rst @@ -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. -.. _howto-howto_rheo_clemmer: +.. _howto_rheo_clemmer: **(Clemmer)** Clemmer, Pierce, O'Connor, Nevins, Jones, Lechman, Tencer, Appl. Math. Model., 130, 310-326 (2024). diff --git a/doc/src/bond_rheo_shell.rst b/doc/src/bond_rheo_shell.rst index 7f6bab1eab..713bd84136 100644 --- a/doc/src/bond_rheo_shell.rst +++ b/doc/src/bond_rheo_shell.rst @@ -10,10 +10,13 @@ Syntax 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:: + *t/form* value = formation time for a bond (time units) + *store/local* values = fix_ID N attributes ... * fix_ID = ID of associated internal fix to store data * 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/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 """""""" .. 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_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 """"""""""" -.. versionadded:: 4May2022 +.. versionadded:: TBD -The *bpm/spring* bond style computes forces based on -deviations from the initial reference state of the two atoms. The -reference state is stored by each bond when it is first computed in +The *rheo/shell* bond style is designed to work with +:doc:`fix rheo/oxidation ` which creates candidate +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 `. +To only count enabled bonds, use the *nbond/shell* attribute in +:doc:`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 `). +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 written to :doc:`binary restart files ` 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 -translational and rotational degrees of freedom of a bonded set of -particles based on a model described by Clemmer and Robbins -:ref:`(Clemmer) `. The force has a magnitude of +This bond style is based on a model described in Ref. +:ref:`(Clemmer) `. The force has a magnitude of .. 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 and :math:`r_0` is the initial distance between the two particles, and -:math:`w` is an optional smoothing factor discussed below. Bonds will -break at a strain of :math:`\epsilon_c`. This is done by setting -the bond type to 0 such that forces are no longer computed. +:math:`\epsilon_c` is maximum strain beyond which a bond breaks. This +is done by setting the bond type to 0 such that forces are no longer +computed. An additional damping force is applied to the bonded 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 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 :doc:`bond_coeff ` command as in the example above, or in 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:`\gamma` (force/velocity units) -If the *normalize* keyword is set to *yes*, the elastic bond force will be -normalized by :math:`r_0` such that :math:`k` must be given in force units. - -By default, pair forces are not calculated between bonded particles. -Pair forces can alternatively be overlaid on top of bond forces by setting -the *overlay/pair* keyword to *yes*. These settings require specific -:doc:`special_bonds ` settings described in the -restrictions. Further details can be found in the :doc:`how to ` -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. +Unlike other BPM-style bonds, this bond style does not update special +bond settings when bonds are created or deleted. This bond style also +does not enforce specific :doc:`special_bonds ` settings. +This behavior is purposeful such :doc:`RHEO pair forces ` +and heat flows are still calculated. 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 @@ -187,39 +170,25 @@ extra quantity can be accessed by the 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 ` 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 ` 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 """""""""""""""" -:doc:`bond_coeff `, :doc:`pair bpm/spring ` +:doc:`bond_coeff `, :doc:`fix rheo/oxidation ` 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: diff --git a/doc/src/compute_rheo_property_atom.rst b/doc/src/compute_rheo_property_atom.rst index 7f5de17c3b..8d1cf47f58 100644 --- a/doc/src/compute_rheo_property_atom.rst +++ b/doc/src/compute_rheo_property_atom.rst @@ -27,8 +27,8 @@ Syntax .. parsed-literal:: - *phase* = atom phase status - *chi* = atom phase neighborhood metric + *phase* = atom phase state + *chi* = atom local phase metric *surface* = atom surface status *surface/r* = atom distance from the surface *surface/divr* = divergence of position at atom position @@ -45,6 +45,7 @@ Syntax *status* = atom full status *rho* = atom density *grad/v/\** = atom velocity gradient + *nbond/shell* = number of oxide bonds Examples """""""" @@ -52,21 +53,34 @@ Examples .. code-block:: LAMMPS compute 1 all rheo/property/atom phase surface/r pressure + compute 2 all rheo/property/atom shift/v/x grad/v/xx Description """"""""""" -Define a computation that simply stores atom attributes specific to the -RHEO package for each atom in the group. This is useful so that the -values can be used by other :doc:`output commands ` that -take computes as inputs. See for example, the :doc:`compute reduce -`, :doc:`fix ave/atom `, :doc:`fix -ave/histo `, :doc:`fix ave/chunk `, -and :doc:`atom-style variable ` commands. +.. versionadded:: TBD + +Define a computation that stores atom attributes specific to the RHEO +package for each atom in the group. This is useful so that the values +can be used by other :doc:`output commands ` that take +computes as inputs. See for example, the +:doc:`compute reduce `, +:doc:`fix ave/atom `, +:doc:`fix ave/histo `, +:doc:`fix ave/chunk `, and +:doc:`atom-style variable ` commands. The possible attributes are described in more detail in other RHEO doc -pages include :doc:`fix rheo `, :doc:`pair rheo `, -and :doc:`the RHEO howto page `. +pages including :doc:`the RHEO howto page `. Many +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 below. Zeroes are stored for atoms not in the specified group or for @@ -98,7 +112,8 @@ Related commands :doc:`fix rheo/viscosity `, :doc:`fix rheo/pressure `, :doc:`fix rheo/thermal `, -:doc:`pair rheo ` +:doc:`fix rheo/oxdiation `, +:doc:`fix rheo ` Default """"""" diff --git a/doc/src/fix_rheo.rst b/doc/src/fix_rheo.rst index 76d4ae3972..f48dd3fcde 100644 --- a/doc/src/fix_rheo.rst +++ b/doc/src/fix_rheo.rst @@ -43,6 +43,8 @@ Examples Description """"""""""" +.. versionadded:: TBD + Perform time integration for RHEO particles, updating positions, velocities, and densities. For a detailed breakdown of the integration timestep and numerical details, see :ref:`(Palermo) `. For an diff --git a/doc/src/fix_rheo_oxidation.rst b/doc/src/fix_rheo_oxidation.rst index a747ec582f..42e2b6d59e 100644 --- a/doc/src/fix_rheo_oxidation.rst +++ b/doc/src/fix_rheo_oxidation.rst @@ -8,43 +8,42 @@ Syntax .. 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 ` command * rheo/oxidation = style name of this fix command * cut = maximum bond length (distance units) * btype = type of bonds created +* rsurf = distance from surface to create bonds (distance units) Examples """""""" .. 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 """"""""""" -This fix... +.. versionadded:: TBD -Each list consists of a series of type -ranges separated by commas. The range can be specified as a -single numeric value, or a wildcard asterisk can be used to specify a range -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. +This fix dynamically creates bonds on the surface of fluids to +represent physical processes such as oxidation. It is intended +for use with bond style :doc:`bond rheo/shell `. -While the *Tfreeze* keyword is optional, the *conductivity* and -*specific/heat* keywords are mandatory. +Every timestep, particles check neighbors within a distance of *cut*. +This distance must be smaller than the kernel length defined in +:doc:`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) `. -Multiple instances of this fix may be defined to apply different -properties to different groups. However, the union of fix groups -across all instances of fix rheo/thermal must cover all atoms. -If there are multiple instances of this fix, any intersections in -the fix groups will lead to incorrect thermal integration. +If used in conjunction with solid bodies, such as those generated +by the *react* option of :doc:`fix rheo/thermal `, +it is recommended that one uses a :doc:`hybrid bond style ` +with different bond types for solid and oxide bonds. Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" @@ -58,10 +57,8 @@ the :doc:`run ` command. This fix is not invoked during :doc:`energy minim Restrictions """""""""""" -This fix must be used with an atom style that includes temperature, -heatflow, and conductivity such as atom_tyle rheo/thermal This fix -must be used in conjuction with :doc:`fix rheo ` with the -*thermal* setting. +This fix must be used with an bond style :doc:`rheo/shell ` +and :doc:`fix rheo ` with surface detection enabled. This fix is part of the RHEO package. It is only enabled if LAMMPS was built with that package. See the :doc:`Build package ` page for more info. @@ -70,12 +67,16 @@ Related commands """""""""""""""" :doc:`fix rheo `, -:doc:`fix rheo/viscosity `, -:doc:`fix rheo/pressure `, -:doc:`pair rheo `, +:doc:`bond rheo/shell `, :doc:`compute rheo/property/atom ` Default """"""" none + +---------- + +.. _howto_rheo_clemmer: + +**(Clemmer)** Clemmer, Pierce, O'Connor, Nevins, Jones, Lechman, Tencer, Appl. Math. Model., 130, 310-326 (2024). diff --git a/doc/src/fix_rheo_pressure.rst b/doc/src/fix_rheo_pressure.rst index 2edd703299..f2a994da88 100644 --- a/doc/src/fix_rheo_pressure.rst +++ b/doc/src/fix_rheo_pressure.rst @@ -33,6 +33,8 @@ Examples Description """"""""""" +.. versionadded:: TBD + This fix defines a pressure equation of state for RHEO particles. One can define different equations of state for different atom types, but an equation must be specified for every atom type. diff --git a/doc/src/fix_rheo_thermal.rst b/doc/src/fix_rheo_thermal.rst index 2ece5521bc..63d9f817ad 100644 --- a/doc/src/fix_rheo_thermal.rst +++ b/doc/src/fix_rheo_thermal.rst @@ -48,6 +48,8 @@ Examples Description """"""""""" +.. versionadded:: TBD + This fix performs time integration of temperature evolution for atom style rheo/thermal. In addition, it defines multiple thermal properties of particles and handles melting/solidification, if applicable. For more details diff --git a/doc/src/fix_rheo_viscosity.rst b/doc/src/fix_rheo_viscosity.rst index 64c6539acc..8c403f8d0b 100644 --- a/doc/src/fix_rheo_viscosity.rst +++ b/doc/src/fix_rheo_viscosity.rst @@ -32,6 +32,8 @@ Examples Description """"""""""" +.. versionadded:: TBD + This fix defines a viscosity for RHEO particles. One can define different viscosities for different atom types, but a viscosity must be specified for every atom type. diff --git a/doc/src/pair_rheo.rst b/doc/src/pair_rheo.rst index f6c3d9e3ba..9ad5586b97 100644 --- a/doc/src/pair_rheo.rst +++ b/doc/src/pair_rheo.rst @@ -31,6 +31,8 @@ Examples Description """"""""""" +.. versionadded:: TBD + Pair style *rheo* computes pressure and viscous forces between particles in the :doc:`rheo package `. If thermal evolution is turned on in :doc:`fix rheo `, then the pair style also calculates diff --git a/doc/src/pair_rheo_react.rst b/doc/src/pair_rheo_react.rst deleted file mode 100644 index 6e7eb49c9d..0000000000 --- a/doc/src/pair_rheo_react.rst +++ /dev/null @@ -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 ` command as in the example above, -or in the data file or restart files read by the -:doc:`read_data ` or :doc:`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 ` -shift, table, and tail options. - -This style does not write information to :doc:`binary restart files `. 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 ` 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 ` page for more info. - -Related commands -"""""""""""""""" - -:doc:`fix rheo `, -:doc:`compute rheo/property/atom ` - -Default -""""""" - -none diff --git a/doc/src/pair_rheo_solid.rst b/doc/src/pair_rheo_solid.rst index ee86992776..4bac9b726f 100644 --- a/doc/src/pair_rheo_solid.rst +++ b/doc/src/pair_rheo_solid.rst @@ -21,6 +21,8 @@ Examples Description """"""""""" +.. versionadded:: TBD + Style *rheo/solid* is effectively a copy of pair style :doc:`bpm/spring ` except it only applies forces between solid RHEO particles, determined by checking the status of diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp index 61b83a4542..11b03623b1 100644 --- a/src/RHEO/compute_rheo_property_atom.cpp +++ b/src/RHEO/compute_rheo_property_atom.cpp @@ -57,7 +57,8 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a if (nvalues == 1) size_peratom_cols = 0; 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 // 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++) { i = iarg-3; - if (strcmp(arg[iarg],"phase") == 0) { + if (strcmp(arg[iarg], "phase") == 0) { 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; - } else if (strcmp(arg[iarg],"chi") == 0) { + } else if (strcmp(arg[iarg], "chi") == 0) { interface_flag = 1; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_chi; - } else if (strcmp(arg[iarg],"surface") == 0) { + } else if (strcmp(arg[iarg], "surface") == 0) { surface_flag = 1; 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; 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; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_divr; } else if (utils::strmatch(arg[iarg], "^surface/n")) { surface_flag = 1; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_n; 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; - } else if (strcmp(arg[iarg],"pressure") == 0) { + } else if (strcmp(arg[iarg], "pressure") == 0) { pressure_flag = 1; 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; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_cv; } 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()); pack_choice[i] = &ComputeRHEOPropertyAtom::pack_atom_style; thermal_flag = 1; - } else if (strcmp(arg[iarg],"nbond/shell") == 0) { + } else if (strcmp(arg[iarg], "nbond/shell") == 0) { shell_flag = 1; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_nbond_shell; } else { @@ -124,9 +127,9 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a arg[iarg], atom->get_style()); pack_choice[i] = &ComputeRHEOPropertyAtom::pack_atom_style; - if (strcmp(arg[iarg],"temperature") == 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], "temperature") == 0) thermal_flag = 1; + if (strcmp(arg[iarg], "heatflow") == 0) thermal_flag = 1; + if (strcmp(arg[iarg], "conductivity") == 0) thermal_flag = 1; } } @@ -149,7 +152,7 @@ ComputeRHEOPropertyAtom::~ComputeRHEOPropertyAtom() void ComputeRHEOPropertyAtom::init() { 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(fixes[0]); 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) { int *nbond = fix_oxidation->nbond; @@ -490,11 +508,11 @@ int ComputeRHEOPropertyAtom::get_tensor_index(char* option) int ComputeRHEOPropertyAtom::get_vector_index(char* option) { int index; - if (utils::strmatch(option,"x$")) { + if (utils::strmatch(option, "x$")) { index = 0; - } else if (utils::strmatch(option,"y$")) { + } else if (utils::strmatch(option, "y$")) { index = 1; - } else if (utils::strmatch(option,"z$")) { + } else if (utils::strmatch(option, "z$")) { if (domain->dimension == 2) error->all(FLERR, "Invalid compute rheo/property/atom property {} in 2D", option); index = 2; diff --git a/src/RHEO/compute_rheo_property_atom.h b/src/RHEO/compute_rheo_property_atom.h index fc5c5d1bd0..a61455f0c5 100644 --- a/src/RHEO/compute_rheo_property_atom.h +++ b/src/RHEO/compute_rheo_property_atom.h @@ -35,7 +35,8 @@ class ComputeRHEOPropertyAtom : public Compute { private: 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 *col_index; double *buf; @@ -55,6 +56,7 @@ class ComputeRHEOPropertyAtom : public Compute { void pack_shift_v(int); void pack_gradv(int); void pack_pressure(int); + void pack_viscosity(int); void pack_nbond_shell(int); void pack_atom_style(int); diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 04b7d33541..98382b07b5 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -19,6 +19,7 @@ #include "fix_rheo.h" #include "atom.h" +#include "citeme.h" #include "compute_rheo_grad.h" #include "compute_rheo_interface.h" #include "compute_rheo_surface.h" @@ -37,6 +38,9 @@ using namespace LAMMPS_NS; using namespace RHEO_NS; using namespace FixConst; +static const char cite_rheo[] = + "TBD\n\n"; + /* ---------------------------------------------------------------------- */ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : @@ -137,6 +141,8 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : } iarg += 1; } + + if (lmp->citeme) lmp->citeme->add(cite_rheo); } /* ---------------------------------------------------------------------- */ diff --git a/src/RHEO/fix_rheo_oxidation.cpp b/src/RHEO/fix_rheo_oxidation.cpp index f536df7842..bc31593653 100644 --- a/src/RHEO/fix_rheo_oxidation.cpp +++ b/src/RHEO/fix_rheo_oxidation.cpp @@ -20,6 +20,7 @@ #include "atom.h" #include "atom_vec.h" +#include "citeme.h" #include "compute_rheo_surface.h" #include "error.h" #include "fix_rheo.h" @@ -34,6 +35,18 @@ using namespace RHEO_NS; using namespace FixConst; 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) : @@ -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); cutsq = cut * cut; + + if (lmp->citeme) lmp->citeme->add(cite_rheo_oxide); } /* ---------------------------------------------------------------------- */ diff --git a/src/RHEO/fix_rheo_stress.cpp b/src/RHEO/fix_rheo_stress.cpp deleted file mode 100644 index b391527f1c..0000000000 --- a/src/RHEO/fix_rheo_stress.cpp +++ /dev/null @@ -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 -#include - -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(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++]; -} diff --git a/src/RHEO/fix_rheo_stress.h b/src/RHEO/fix_rheo_stress.h deleted file mode 100644 index 4bf522793f..0000000000 --- a/src/RHEO/fix_rheo_stress.h +++ /dev/null @@ -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 diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index a133428d4a..9fbdb8c8f6 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -20,6 +20,7 @@ #include "atom.h" #include "atom_vec.h" +#include "citeme.h" #include "comm.h" #include "compute_rheo_grad.h" #include "compute_rheo_vshift.h" @@ -43,6 +44,18 @@ using namespace RHEO_NS; using namespace FixConst; 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) : @@ -193,6 +206,8 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : 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); } + + if (lmp->citeme) lmp->citeme->add(cite_rheo_oxide); } /* ---------------------------------------------------------------------- */ diff --git a/src/RHEO/pair_rheo_react.cpp b/src/RHEO/pair_rheo_react.cpp deleted file mode 100644 index 4709ea169e..0000000000 --- a/src/RHEO/pair_rheo_react.cpp +++ /dev/null @@ -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( - 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( - 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(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++]; - } -} diff --git a/src/RHEO/pair_rheo_react.h b/src/RHEO/pair_rheo_react.h deleted file mode 100644 index 88d5dbeb0e..0000000000 --- a/src/RHEO/pair_rheo_react.h +++ /dev/null @@ -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