Adding more oxidation features + doc pages

doc/src/Howto_rheo.rst

@@ -1,4 +1,99 @@
 Reproducing hydrodynamics and elastic objects (RHEO)
-======================
+====================================================
 
-Text
+The RHEO package is built around an implementation of smoothed particle
+hydrodynamics (SPH) coupled to the :doc:`BPM package <Howto_bpm>` to model
+solid elements of a system. The SPH solver supports many advanced options
+including reproducing kernels, particle shifting, free surface identification,
+and solid surface reconstruction. To model fluid-solid systems, the status of
+particles can dynamically change between a fluid and solid state, e.g. during
+melting/solidification, which determines how they interact and their physical
+behavior. The package is designed with modularity in mind, so one can easily
+turn various features on/off, adjust physical details of the system, or
+develop new capabilities. Additional numerical details can be found in
+:ref:`(Palermo) <howto_rheo_palermo>` and :ref:`(Clemmer) <howto_rheo_clemmer>`.
+
+----------
+
+At the core of the package is :doc:`fix rheo <fix_rheo>` which integrates
+particle trajectories and controls many optional features (e.g. the use
+of reproducing kernels). In conjunction to fix rheo, one must specify an
+instance of :doc:`fix rheo/pressure <fix_rheo_pressure>` and
+:doc:`fix rheo/viscosity <fix_rheo_viscosity>` to define a pressure equation
+of state and viscosity model, respectively. Optionally, one can model
+a heat equation with :doc:`fix rheo/thermal`, which also allows the user
+to specify equations for a particle's thermal conductivity,  specific heat,
+latent heat, and melting temperature. Fix rheo must be defined prior to all
+other RHEO fixes.
+
+Typically, RHEO requires atom style rheo. In addition to typical atom
+properties like positions and forces, particles store a local density,
+viscosity, pressure, and status. If thermal evolution is modeled, one must
+use atom style rheo/thermal which also include a local temperature and
+conductivity. The status variable uses bitmasking to track various
+properties of a particle such as its current phase (fluid or solid) and its
+location relative to a surface. Many of these properties (and others) can
+be easily accessed using
+:doc:`compute rheo/property/atom <fix_rheo_property_atom>`.
+
+Fluid interactions, including pressure forces, viscous forces, and heat exchange,
+are calculated using :doc:`pair rheo <pair_rheo>`. Unlike typical pair styles,
+pair rheo ignores the :doc:`special bond <special_bonds>` settings. Instead,
+it determines whether to calculate forces based on the status of particles:
+hydrodynamic forces are only calculated if a fluid particle is involved.
+
+----------
+
+To model elastic objects, there are current two mechanisms in RHEO, one designed
+for bulk solid bodies and the other for thin shells. Both mechanisms rely on
+overlaying bonds and therefore require a hybrid of atom style bond and rheo
+(or rheo/thermal).
+
+To create an elastic solid body, one has to (a) change the status of constituent
+particles to solid (e.g. with the :doc:`set <set>` command), (b) create bpm
+bonds between the particles (see the :doc:`bpm howto <Howto_bpm>` page for
+more details), and (c) use :doc:`pair rheo/solid <pair_rheo_solid>` to
+apply repulsive contact forces between distinct solid bodies. Akin to pair rheo,
+looks at a particles fluid/solid status to determine whether to apply forces.
+However, unlike pair rheo, pair rheo/solid does obey special bond settings such
+that contact forces do not have to be calculated between two bonded solid particles
+in the same elastic body.
+
+In systems with thermal evolution, fix rheo/thermal can optionally set a
+melting/solidification temperature allowing particles to dynamically swap their
+state between fluid and solid. Using the *react* option, one can specify a maximum
+bond length and a bond type. Then, when solidifying, particles will search their
+local neighbors and automatically create bonds with any neighboring solid particles
+in range. For BPM bond styles, bonds will then use the immediate position of the two
+particles to calculate a reference state. When melting, particles will then delete
+any bonds of the specified type when reverting to a fluid state. Special bonds are
+updated as bonds are created/broken.
+
+The other option for elastic objects is an elastic shell that is nominally much
+thinner than a particle diameter, e.g. a oxide skin which gradually forms over time
+on the surface of a fluid. Currently, this is implemented using
+:doc:`fix rheo/oxidaton <fix_rheo_oxidation>` and bond style
+:doc:`rheo/shell <bond_rheo_shell>`. Essentially, fix rheo/oxidaton creates candidate
+bonds of a specified type between surface fluid particles within a specified distance.
+a newly created rheo/shell bond will then start a timer. While the timer is counting
+down, the bond will delete itself if particles move too far apart or move away from the
+surface. However, if the timer reaches a user-defined threshold, then the bond will
+activate and apply additional forces to the fluid particles. Bond style rheo/shell
+then operates very similarly to a BPM bond style, storing a reference length and
+breaking if stretched too far. Unlike the above method, this option does not remove
+the underlying fluid interactions (although particle shifting is turned off) and does
+not modify special bond settings of particles.
+
+While these two options are not expected to be appropriate for every multiphase system,
+either framework can be modified to create more suitable models (e.g. by changing the
+criteria for creating/deleting a bond or altering force calculations).
+
+----------
+
+.. _howto-howto_rheo_palermo:
+
+**(Palermo)** Palermo, Clemmer, Wolf, O'Connor, in preparation.
+
+.. _howto-howto_rheo_clemmer:
+
+**(Clemmer)** Clemmer, Pierce, O'Connor, Nevins, Jones, Lechman, Tencer, Appl. Math. Model., 130, 310-326 (2024).

doc/src/Packages_details.rst

@@ -101,6 +101,7 @@ page gives those details.
    * :ref:`QEQ <PKG-QEQ>`
    * :ref:`QMMM <PKG-QMMM>`
    * :ref:`QTB <PKG-QTB>`
+   * :ref:`RHEO <PKG-RHEO>`
    * :ref:`REACTION <PKG-REACTION>`
    * :ref:`REAXFF <PKG-REAXFF>`
    * :ref:`REPLICA <PKG-REPLICA>`
@@ -2571,6 +2572,41 @@ another set.
 
 ----------
 
+.. _PKG-RHEO:
+
+RHEO package
+------------
+
+**Contents:**
+
+Pair styles, bond styles, fixes, and computes for reproducing
+hydrodynamics and elastic objects. See the
+:doc:`Howto rheo <Howto_rheo>` page for an overview.
+
+**Authors:** Joel T. Clemmer (Sandia National Labs),
+Thomas C. O'Connor (Carnegie Mellon University)
+
+.. versionadded:: TBD
+
+**Supporting info:**
+
+* src/RHEO filenames -> commands
+* :doc:`Howto_rheo <Howto_rheo>`
+* :doc:`atom_style rheo <atom_style>`
+* :doc:`atom_style rheo/thermal <atom_style>`
+* :doc:`bond_style rheo/shell <bond_rheo_shell>`
+* :doc:`compute rheo/property/atom <compute_rheo_property_atom>`
+* :doc:`fix rheo <fix_rheo>`
+* :doc:`fix rheo/oxidation <fix_rheo_oxidation>`
+* :doc:`fix rheo/pressure <fix_rheo_pressure>`
+* :doc:`fix rheo/thermal <fix_rheo_thermal>`
+* :doc:`fix rheo/viscosity <fix_rheo_viscosity>`
+* :doc:`pair_style rheo <pair_rheo>`
+* :doc:`pair_style rheo/solid <pair_rheo_solid>`
+* examples/rheo
+
+----------
+
 .. _PKG-RIGID:
 
 RIGID package

doc/src/Packages_list.rst

@@ -403,6 +403,11 @@ whether an extra library is needed to build and use the package:
      - :doc:`fix qtb <fix_qtb>` :doc:`fix qbmsst <fix_qbmsst>`
      - qtb
      - no
+   * - :ref:`RHEO <PKG-RHEO>`
+     - reproducing hydrodynamics and elastic objects
+     - :doc:`Howto rheo <Howto_rheo>`
+     - rheo
+     - no
    * - :ref:`REACTION <PKG-REACTION>`
      - chemical reactions in classical MD
      - :doc:`fix bond/react <fix_bond_react>`

doc/src/fix_rheo.rst

@@ -8,40 +8,107 @@ Syntax
 
 .. parsed-literal::
 
-   fix ID group-ID rheo cut kstyle keyword values...
+   fix ID group-ID rheo cut kstyle zmin keyword values...
 
 * ID, group-ID are documented in :doc:`fix <fix>` command
 * rheo = style name of this fix command
-* cut = *quintic* or *RK0* or *RK1* or *RK2*
+* cut = cutoff for the kernel (distance)
+* kstyle = *quintic* or *RK0* or *RK1* or *RK2*
+* zmin = minimal number of neighbors for reproducing kernels
 * zero or more keyword/value pairs may be appended to args
-* keyword = *shift* or *thermal* or *surface/detection* or *interface/reconstruction* or
+* keyword = *thermal* or *interface/reconstruct* or *surface/detection* or
-            *rho/sum* or *density* or *sound/squared*
+            *shift* or *rho/sum* or *density* or *speed/sound*
 
   .. parsed-literal::
 
-       *shift* values = none, turns on velocity shifting
        *thermal* values = none, turns on thermal evolution
+       *interface/reconstruct* values = none, reconstructs interfaces with solid particles
        *surface/detection* values = *sdstyle* *limit* *limit/splash*
          *sdstyle* = *coordination* or *divergence*
-         *limit* = threshold for surface particles (unitless)
+         *limit* = threshold for surface particles
-         *limit/splash* = threshold for splash particles (unitless)
+         *limit/splash* = threshold for splash particles
-       *interface/reconstruct* values = none, reconstructs interfaces with solid particles
+       *shift* values = none, turns on velocity shifting
        *rho/sum* values = none, uses the kernel to compute the density of particles
-       *density* values = *rho0* (density)
+       *density* values = *rho01*, ... *rho0N* (density)
-       *sound/squared* values = *csq* (velocity\^2)
+       *speed/sound* values = *cs0*, ... *csN* (velocity)
 
 Examples
 """"""""
 
 .. code-block:: LAMMPS
 
-   fix 1 all rheo 1.0 quintic thermal density 0.1 sound/squared 10.0
+   fix 1 all rheo 1.0 quintic thermal density 0.1 speed/sound 10.0
    fix 1 all rheo 1.0 RK1 shift surface/detection coordination 40
 
 Description
 """""""""""
 
-Fix description...
+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) <howto_rheo_palermo>`. For an
+overview of other features available in the RHEO package, see
+:doc:`the RHEO howto <Howto_rheo>`.
+
+The type of kernel is specified using *kstyle* and the cutoff is *cut*. Four
+kernels are currently available. The *quintic* kernel is a standard quintic
+spline function commonly used in SPH. The other options, *RK0*, *RK1*, and
+*RK2*, are zeroth, first, and second order reproducing. To generate a reproducing kernel, a particle must have sufficient neighbors inside the
+kernel cutoff distance (a coordination number) to accurately calculate
+moments. This threshold is set by *zmin*. If reproducing kernels are
+requested but a particle has fewer neighbors, then it will revert to a
+non-reproducing quintic kernel until it gains more neighbors.
+
+To model temperature evolution, one must specify the *thermal* keyword,
+define a separate instance of :doc:`fix rheo/thermal <fix_rheo_thermal>`,
+and use atom style rheo/thermal.
+
+By default, the density of solid RHEO particles does not evolve and forces
+with fluid particles are calculated using the current velocity of the solid
+particle. If the *interface/reconstruct* keyword is used, then the density
+and velocity of solid particles are alternatively reconstructed for every
+fluid-solid interaction to ensure no-slip and pressure-balanced boundaries.
+This is done by estimating the location of the fluid-solid interface and
+extrapolating fluid particle properties across the interface to calculate a
+temporary apparent density and velocity for a solid particle. The numerical
+details are the same as those described in
+:ref:`(Palermo) <howto_rheo_palermo>` except there is an additional
+restriction that the reconstructed solid density cannot be less than the
+equilibrium density. This prevents fluid particles from sticking to solid
+surfaces.
+
+A modified form of Fickian particle shifting can be enabled with the
+*shift* keyword. This effectively shifts particle positions to generate a
+more uniform spatial distribution. In systems with free surfaces, the
+*surface/detection* keyword can be used to classify the location of
+particles as being within the bulk fluid, on a free surface, or isolated
+from other particles in a splash or droplet. Shifting is then disabled in
+the direction away from the free surface to prevent it from diffusing
+particles away from the bulk fluid. Surface detection can also be used
+to control surface-nucleated effects like oxidation when used in combination
+with :doc:`fix rheo/oxidation <fix_rheo_oxidation>`.
+
+The *surface/detection* keyword takes three arguments: *sdstyle*, *limit*,
+and *limi/splash*. The first, *sdstyle*, specifies whether surface particles
+are identified using a coordination number (*coordination*) or the divergence
+of the local particle positions (*divergence*). The threshold value for a
+surface particle for either of these criteria is set by the numerical value
+of *limit*. Additionally, if a particle's coordination number is too low,
+i.e. if it has separated off from the bulk in a droplet, it is not possible
+to define surfaces and a particle is classified as a splash. The coordination
+threshold for this classification is set by the numerical value of
+*limit/splash*.
+
+By default, RHEO integrates particles' densities using a mass diffusion
+equation. Alternatively, one can update densities every timestep by performing
+a kernel summation of the masses of neighboring particles by specifying the *rho/sum* keyword.
+
+The *density* is used to specify the equilbrium density of each of the N
+particle types. It must be followed by N numerical values specifying each
+type's equilibrium density *rho0*.
+
+The *density* is used to specify the speed of sound of each of the N particle
+types. It must be followed by N numerical values specifying each type's speed
+of sound *cs*.
 
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
@@ -55,13 +122,14 @@ the :doc:`run <run>` command.  This fix is not invoked during :doc:`energy minim
 Restrictions
 """"""""""""
 
-This fix must be used with an atom style that includes density
+This fix must be used with atom style rheo or rheo/thermal.
-such as atom_style rheo or rheo/thermal. This fix must be used in
+This fix must be used in conjuction with
-conjuction with :doc:`fix rheo/pressure <fix_rheo_pressure>`. and
+:doc:`fix rheo/pressure <fix_rheo_pressure>`. and
 :doc:`fix rheo/viscosity <fix_rheo_viscosity>`, If the *thermal*
 setting is used, there must also be an instance of
 :doc:`fix rheo/thermal <fix_rheo_thermal>`. The fix group must be
-set to all.
+set to all. Only one instance of fix rheo may be defined and it
+must be defined prior to all other RHEO fixes.
 
 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.
@@ -78,4 +146,10 @@ Related commands
 Default
 """""""
 
-*rho0* and *csq* are set to 1.0.
+*rho0* and *cs* are set to 1.0 for all atom types.
+
+----------
+
+.. _howto-howto_rheo_palermo:
+
+**(Palermo)** Palermo, Clemmer, Wolf, O'Connor, in preparation.

doc/src/pair_rheo.rst

@@ -16,8 +16,8 @@ Syntax
 
 .. parsed-literal::
 
-     *rho/damp* args = density damping prefactor :math:`\xi` (units?)
+     *rho/damp* args = density damping prefactor :math:`\xi`
-     *artificial/visc* args = artificial viscosity prefactor :math:`\zeta` (units?)
+     *artificial/visc* args = artificial viscosity prefactor :math:`\zeta`
      *harmonic/means* args = none
 
 Examples
@@ -31,7 +31,27 @@ Examples
 Description
 """""""""""
 
-pair style...
+Pair style *rheo* computes pressure and viscous forces between particles
+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
+heat exchanged between particles.
+
+The *artificial/viscosity* keyword is used to specify the magnitude
+:math:`\zeta` of an optional artificial viscosity contribution to forces.
+This factor can help stabilize simulations by smoothing out small length
+scale variations in velocity fields.
+
+The *rho/damp* keyword is used to specify the magnitude :math:`\xi` of
+an optional pairwise damping term between the density of particles. This
+factor can help stabilize simulations by smoothing out small length
+scale variations in density fields.
+
+If particles have different viscosities or conductivities, the
+*harmonic/means* keyword changes how they are averaged before calculating
+pairwise forces or heat exchanges. By default, an arithmetic averaged is
+used, however, a harmonic mean may improve stability in multiphase systems
+with large disparities in viscosities. This keyword has no effect on
+results if viscosities and conductivities are constant.
 
 No coefficients are defined for each pair of atoms types via the
 :doc:`pair_coeff <pair_coeff>` command as in the examples

src/RHEO/atom_vec_rheo.cpp

@@ -13,7 +13,7 @@
 
 /* ----------------------------------------------------------------------
    Contributing authors:
-   Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU)
+   Joel Clemmer (SNL)
 ----------------------------------------------------------------------- */
 
 #include "atom_vec_rheo.h"

src/RHEO/atom_vec_rheo_thermal.cpp

@@ -13,7 +13,7 @@
 
 /* ----------------------------------------------------------------------
    Contributing authors:
-   Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU)
+   Joel Clemmer (SNL)
 ----------------------------------------------------------------------- */
 
 #include "atom_vec_rheo_thermal.h"

src/RHEO/bond_rheo_shell.cpp

@@ -11,6 +11,11 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
+/* ----------------------------------------------------------------------
+   Contributing authors:
+   Joel Clemmer (SNL)
+----------------------------------------------------------------------- */
+
 #include "bond_rheo_shell.h"
 
 #include "atom.h"
@@ -40,17 +45,38 @@ BondRHEOShell::BondRHEOShell(LAMMPS *_lmp) :
     BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr)
 {
   partial_flag = 1;
+  comm_reverse = 1;
 
   tform = rmax = -1;
 
   single_extra = 1;
   svector = new double[1];
+
+  // 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("shell_nbond", tmp1, tmp2);
+  if (index_nb == -1) {
+    id_fix = utils::strdup("bond_rheo_shell_fix_property_atom");
+    modify->add_fix(fmt::format("{} all property/atom i_shell_nbond", id_fix));
+    index_nb = atom->find_custom("shell_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");
 }
 
 /* ---------------------------------------------------------------------- */
 
 BondRHEOShell::~BondRHEOShell()
 {
+  if (modify->nfix) modify->delete_fix(id_fix);
+  delete[] id_fix;
   delete[] svector;
 
   if (allocated) {
@@ -59,6 +85,8 @@ BondRHEOShell::~BondRHEOShell()
     memory->destroy(ecrit);
     memory->destroy(gamma);
   }
+
+  memory->destroy(dbond);
 }
 
 /* ----------------------------------------------------------------------
@@ -151,6 +179,15 @@ void BondRHEOShell::compute(int eflag, int vflag)
 
   double **bondstore = fix_bond_history->bondstore;
 
+  if (atom->nmax > nmax_store){
+    nmax_store = atom->nmax;
+    memory->destroy(dbond);
+    memory->create(dbond, nmax_store, "rheo/shell:dbond");
+  }
+
+  size_t nbytes = nmax_store * sizeof(int);
+  memset(&dbond[0], 0, nbytes);
+
   for (n = 0; n < nbondlist; n++) {
 
     // skip bond if already broken
@@ -196,6 +233,8 @@ void BondRHEOShell::compute(int eflag, int vflag)
       if (bondstore[n][1] >= tform) {
         bondstore[n][0] = r;
         r0 = r;
+        if (newton_bond || i1 < nlocal) dbond[i1] ++;
+        if (newton_bond || i2 < nlocal) dbond[i2] ++;
       } else {
         continue;
       }
@@ -205,6 +244,8 @@ void BondRHEOShell::compute(int eflag, int vflag)
     if (fabs(e) > ecrit[type]) {
       bondlist[n][2] = 0;
       process_broken(i1, i2);
+      if (newton_bond || i1 < nlocal) dbond[i1] --;
+      if (newton_bond || i2 < nlocal) dbond[i2] --;
       continue;
     }
 
@@ -233,6 +274,17 @@ void BondRHEOShell::compute(int eflag, int vflag)
 
     if (evflag) ev_tally(i1, i2, nlocal, newton_bond, 0.0, fbond, delx, dely, delz);
   }
+
+
+  // Communicate changes in nbond
+  if (newton_bond) comm->reverse_comm(this);
+
+  for(i = 0; i < nlocal; i++) {
+    nbond[i] += dbond[i];
+
+    // If it has bonds, no shifting
+    if (nbond[i] != 0) status[i] |= STATUS_NO_SHIFT;
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -419,6 +471,34 @@ void BondRHEOShell::read_restart_settings(FILE *fp)
   MPI_Bcast(&rmax, 1, MPI_DOUBLE, 0, world);
 }
 
+
+/* ---------------------------------------------------------------------- */
+
+int BondRHEOShell::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 BondRHEOShell::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++];
+  }
+}
+
 /* ---------------------------------------------------------------------- */
 
 double BondRHEOShell::single(int type, double rsq, int i, int j, double &fforce)

src/RHEO/bond_rheo_shell.h

@@ -37,12 +37,18 @@ class BondRHEOShell : public BondBPM {
   void read_restart(FILE *) override;
   void write_restart_settings(FILE *) override;
   void read_restart_settings(FILE *) override;
+  int pack_reverse_comm(int, int, double *) override;
+  void unpack_reverse_comm(int, int *, double *) override;
   double single(int, double, int, int, double &) override;
 
  protected:
   double *k, *ecrit, *gamma;
   double tform, rmax;
 
+  int *dbond, *nbond;
+  int index_nb, nmax_store;
+  char *id_fix;
+
   void process_ineligibility(int, int);
   void allocate();
   void store_data();

src/RHEO/compute_rheo_grad.cpp

@@ -13,7 +13,7 @@
 
 /* ----------------------------------------------------------------------
    Contributing authors:
-   Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU)
+   Joel Clemmer (SNL), Thomas O'Connor (CMU)
 ----------------------------------------------------------------------- */
 
 #include "compute_rheo_grad.h"

src/RHEO/compute_rheo_kernel.cpp

@@ -13,7 +13,7 @@
 
 /* ----------------------------------------------------------------------
    Contributing authors:
-   Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU)
+   Joel Clemmer (SNL), Thomas O'Connor (CMU)
 ----------------------------------------------------------------------- */
 
 #include "compute_rheo_kernel.h"

src/RHEO/compute_rheo_property_atom.cpp

@@ -56,7 +56,7 @@ 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 = 0;
+  pressure_flag = thermal_flag = interface_flag = surface_flag = shift_flag = shell_flag = 0;
 
   // parse input values
   // customize a new keyword by adding to if statement
@@ -112,6 +112,9 @@ 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) {
+      shell_flag = 1;
+      pack_choice[i] = &ComputeRHEOPropertyAtom::pack_nbond_shell;
     } else {
       avec_index[i] = atom->avec->property_atom(arg[iarg]);
       if (avec_index[i] < 0)
@@ -156,6 +159,8 @@ void ComputeRHEOPropertyAtom::init()
     error->all(FLERR, "Cannot request velocity shifting property without corresponding option in fix rheo");
   if (thermal_flag && !(fix_rheo->thermal_flag))
     error->all(FLERR, "Cannot request thermal property without fix rheo/thermal");
+  if (shell_flag && !(fix_rheo->oxidation_flag))
+    error->all(FLERR, "Cannot request number of shell bonds without fix rheo/oxidation");
 
   compute_interface = fix_rheo->compute_interface;
   compute_kernel = fix_rheo->compute_kernel;

src/RHEO/compute_rheo_property_atom.h

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

src/RHEO/compute_rheo_rho_sum.cpp

@@ -11,6 +11,11 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
+/* ----------------------------------------------------------------------
+   Contributing authors:
+   Joel Clemmer (SNL)
+----------------------------------------------------------------------- */
+
 #include "compute_rheo_rho_sum.h"
 
 #include "atom.h"

src/RHEO/compute_rheo_vshift.cpp

@@ -147,6 +147,8 @@ void ComputeRHEOVShift::compute_peratom()
 
       fluidj = !(status[j] & PHASECHECK);
       if ((!fluidi) && (!fluidj)) continue;
+
+      // Will skip shifting in FixRHEO initial integrate, but also skip here to save time
       if ((status[i] & STATUS_NO_SHIFT) && (status[j] & STATUS_NO_SHIFT)) continue;
 
       dx[0] = xtmp - x[j][0];

src/RHEO/fix_rheo.cpp

@@ -221,7 +221,7 @@ void FixRHEO::setup_pre_force(int /*vflag*/)
 {
   // Check to confirm accessory fixes do not preceed FixRHEO
   // Note: fixes set this flag in setup_pre_force()
-  if (viscosity_fix_defined || pressure_fix_defined || thermal_fix_defined)
+  if (viscosity_fix_defined || pressure_fix_defined || thermal_fix_defined || oxidation_fix_defined)
     error->all(FLERR, "Fix RHEO must be defined before all other RHEO fixes");
 
   // Calculate surfaces

src/RHEO/fix_rheo.h

@@ -55,6 +55,7 @@ class FixRHEO : public Fix {
   int viscosity_fix_defined;
   int pressure_fix_defined;
   int thermal_fix_defined;
+  int oxidation_fix_defined;
 
   class ComputeRHEOGrad *compute_grad;
   class ComputeRHEOKernel *compute_kernel;

src/RHEO/fix_rheo_oxidation.cpp

@@ -60,6 +60,7 @@ FixRHEOOxidation::~FixRHEOOxidation()
 int FixRHEOOxidation::setmask()
 {
   int mask = 0;
+  mask |= PRE_FORCE;
   mask |= POST_INTEGRATE;
   return mask;
 }
@@ -70,7 +71,7 @@ void FixRHEOOxidation::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/oxidation");
-  class FixRHEO *fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
+  fix_rheo = dynamic_cast<FixRHEO *>(fixes[0]);
   double cut_kernel = fix_rheo->h;
 
   if (cut > cut_kernel)
@@ -102,6 +103,26 @@ void FixRHEOOxidation::init_list(int /*id*/, NeighList *ptr)
 
 /* ---------------------------------------------------------------------- */
 
+void FixRHEOOxidation::setup_pre_force(int /*vflag*/)
+{
+  // Not strictly required that this fix be after FixRHEO,
+  // but enforce to be consistent with other RHEO fixes
+  fix_rheo->oxidation_fix_defined = 1;
+
+  if (!fix_rheo->surface_flag) error->all(FLERR,
+      "fix rheo/oxidation requires surface calculation in fix rheo");
+
+  pre_force(0);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixRHEOThermal::pre_force(int /*vflag*/)
+{
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixRHEOOxidation::post_integrate()
 {
   int i, j, n, ii, jj, inum, jnum, bflag;
@@ -187,4 +208,10 @@ void FixRHEOOxidation::post_integrate()
       }
     }
   }
+
+  //todo:
+  // allow option to create near-surface bonds
+  // extract # of bonds in property/atom
+  // check bond style shell used, get index to bonds, share with compute property atom
+  // add type option to compute nbond/atom
 }

src/RHEO/fix_rheo_oxidation.h

@@ -33,6 +33,8 @@ class FixRHEOOxidation : public Fix {
   int setmask() override;
   void init() override;
   void init_list(int, class NeighList *) override;
+  void setup_pre_force(int) override;
+  void pre_force(int) override;
   void post_integrate() override;
 
  private:
@@ -40,6 +42,7 @@ class FixRHEOOxidation : public Fix {
   double cut, cutsq;
   class NeighList *list;
 
+  class FixRHEO *fix_rheo;
   //class FixBondHistory *fix_bond_history;
 };
 

src/RHEO/fix_rheo_pressure.cpp

@@ -13,7 +13,7 @@
 
 /* ----------------------------------------------------------------------
    Contributing authors:
-   Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU)
+   Joel Clemmer (SNL), Thomas O'Connor (CMU)
 ----------------------------------------------------------------------- */
 
 #include "fix_rheo_pressure.h"

src/RHEO/fix_rheo_thermal.cpp

@@ -498,21 +498,36 @@ void FixRHEOThermal::break_bonds()
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
 
-  // Rapidly delete all bonds for local atoms that melt (no shifting)
+  // Rapidly delete all bonds for local atoms that melt of a given type
   for (int i = 0; i < nlocal; i++) {
     if (!(status[i] & STATUS_MELTING)) continue;
-    for (m = 0; m < num_bond[i]; m++) {
+    for (m = (num_bond[i] - 1); m >= 0; m--) {
-      j = atom->map(bond_atom[i][m]);
+      if (bond_type[i][m] != btype) continue;
-      bond_type[i][m] = 0;
 
-      if (n_histories > 0)
+      j = atom->map(bond_atom[i][m]);
-        for (auto &ihistory: histories)
+
-          dynamic_cast<FixBondHistory *>(ihistory)->delete_history(i, m);
+      nmax = num_bond[i] - 1;
+      if (m == nmax) {
+        if (n_histories > 0)
+          for (auto &ihistory: histories)
+            dynamic_cast<FixBondHistory *>(ihistory)->delete_history(i, m);
+      } else {
+        bond_type[i][m] = bond_type[i][nmax];
+        bond_atom[i][m] = bond_atom[i][nmax];
+        if (n_histories > 0) {
+          for (auto &ihistory: histories) {
+            auto fix_bond_history = dynamic_cast<FixBondHistory *>  (ihistory);
+            fix_bond_history->shift_history(i, m, nmax);
+            fix_bond_history->delete_history(i, nmax);
+          }
+        }
+      }
+      bond_type[i][nmax] = 0;
+      num_bond[i]--;
 
       if (fix_update_special_bonds)
         fix_update_special_bonds->add_broken_bond(i, j);
     }
-    num_bond[i] = 0;
   }
 
   // Update bond list and break solid-melted bonds
@@ -530,7 +545,7 @@ void FixRHEOThermal::break_bonds()
     // Delete bonds for non-melted local atoms (shifting)
     if (i < nlocal) {
       for (m = 0; m < num_bond[i]; m++) {
-        if (bond_atom[i][m] == tag[j]) {
+        if (bond_atom[i][m] == tag[j] && bond_type[i][m] == btype) {
           nmax = num_bond[i] - 1;
           bond_type[i][m] = bond_type[i][nmax];
           bond_atom[i][m] = bond_atom[i][nmax];
@@ -549,7 +564,7 @@ void FixRHEOThermal::break_bonds()
 
     if (j < nlocal) {
       for (m = 0; m < num_bond[j]; m++) {
-        if (bond_atom[j][m] == tag[i]) {
+        if (bond_atom[j][m] == tag[i] && bond_type[j][m] == btype) {
           nmax = num_bond[j] - 1;
           bond_type[j][m] = bond_type[j][nmax];
           bond_atom[j][m] = bond_atom[j][nmax];

src/RHEO/fix_rheo_viscosity.cpp

@@ -13,7 +13,7 @@
 
 /* ----------------------------------------------------------------------
    Contributing authors:
-   Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU)
+   Joel Clemmer (SNL), Thomas O'Connor (CMU)
 ----------------------------------------------------------------------- */
 
 #include "fix_rheo_viscosity.h"

src/RHEO/pair_rheo.cpp

@@ -13,7 +13,7 @@
 
 /* ----------------------------------------------------------------------
    Contributing authors:
-   Joel Clemmer (SNL), Thomas O'Connor (CMU), Eric Palermo (CMU)
+   Joel Clemmer (SNL), Thomas O'Connor (CMU)
 ----------------------------------------------------------------------- */
 
 #include "pair_rheo.h"

src/RHEO/pair_rheo_solid.cpp

@@ -11,6 +11,11 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
+/* ----------------------------------------------------------------------
+   Contributing authors:
+   Joel Clemmer (SNL)
+----------------------------------------------------------------------- */
+
 #include "pair_rheo_solid.h"
 
 #include "atom.h"