diff --git a/cmake/Modules/Packages/KOKKOS.cmake b/cmake/Modules/Packages/KOKKOS.cmake
index cf3d19c720..7f23a6f777 100644
--- a/cmake/Modules/Packages/KOKKOS.cmake
+++ b/cmake/Modules/Packages/KOKKOS.cmake
@@ -47,8 +47,8 @@ if(DOWNLOAD_KOKKOS)
   list(APPEND KOKKOS_LIB_BUILD_ARGS "-DCMAKE_CXX_EXTENSIONS=${CMAKE_CXX_EXTENSIONS}")
   list(APPEND KOKKOS_LIB_BUILD_ARGS "-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TOOLCHAIN_FILE}")
   include(ExternalProject)
-  set(KOKKOS_URL "https://github.com/kokkos/kokkos/archive/3.6.00.tar.gz" CACHE STRING "URL for KOKKOS tarball")
-  set(KOKKOS_MD5 "b5c44ea961031795f434002cd7b31c20" CACHE STRING "MD5 checksum of KOKKOS tarball")
+  set(KOKKOS_URL "https://github.com/kokkos/kokkos/archive/3.6.01.tar.gz" CACHE STRING "URL for KOKKOS tarball")
+  set(KOKKOS_MD5 "0ec97fc0c356dd65bd2487defe81a7bf" CACHE STRING "MD5 checksum of KOKKOS tarball")
   mark_as_advanced(KOKKOS_URL)
   mark_as_advanced(KOKKOS_MD5)
   ExternalProject_Add(kokkos_build
@@ -72,7 +72,7 @@ if(DOWNLOAD_KOKKOS)
   add_dependencies(LAMMPS::KOKKOSCORE kokkos_build)
   add_dependencies(LAMMPS::KOKKOSCONTAINERS kokkos_build)
 elseif(EXTERNAL_KOKKOS)
-  find_package(Kokkos 3.6.00 REQUIRED CONFIG)
+  find_package(Kokkos 3.6.01 REQUIRED CONFIG)
   target_link_libraries(lammps PRIVATE Kokkos::kokkos)
   target_link_libraries(lmp PRIVATE Kokkos::kokkos)
 else()
diff --git a/doc/src/compute_sna_atom.rst b/doc/src/compute_sna_atom.rst
index deb717d84d..1b6b200685 100644
--- a/doc/src/compute_sna_atom.rst
+++ b/doc/src/compute_sna_atom.rst
@@ -45,7 +45,7 @@ Syntax
 * w_1, w_2,... = list of neighbor weights, one for each type
 * nx, ny, nz = number of grid points in x, y, and z directions (positive integer)
 * zero or more keyword/value pairs may be appended
-* keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag* or *bikflag* or *switchinnerflag* or *sinner* or *dinner*
+* keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag* or *bikflag* or *switchinnerflag* or *sinner* or *dinner* or *dgradflag*
 
   .. parsed-literal::
 
@@ -68,9 +68,6 @@ Syntax
        *wselfallflag* value = *0* or *1*
           *0* = self-contribution only for element of central atom
           *1* = self-contribution for all elements
-       *bikflag* value = *0* or *1* (only implemented for compute snap)
-          *0* = per-atom bispectrum descriptors are summed over atoms
-          *1* = per-atom bispectrum descriptors are not summed over atoms
        *switchinnerflag* value = *0* or *1*
           *0* = do not use inner switching function
           *1* = use inner switching function
@@ -78,6 +75,12 @@ Syntax
           *sinnerlist* = *ntypes* values of *Sinner* (distance units)
        *dinner* values = *dinnerlist*
           *dinnerlist* = *ntypes* values of *Dinner* (distance units)
+       *bikflag* value = *0* or *1* (only implemented for compute snap)
+          *0* = descriptors are summed over atoms of each type
+          *1* = descriptors are listed separately for each atom
+       *dgradflag* value = *0* or *1* (only implemented for compute snap)
+          *0* = descriptor gradients are summed over atoms of each type
+          *1* = descriptor gradients are listed separately for each atom pair
 
 Examples
 """"""""
@@ -360,15 +363,6 @@ This option is typically used in conjunction with the *chem* keyword,
 and LAMMPS will generate a warning if both *chem* and *bnormflag*
 are not both set or not both unset.
 
-The keyword *bikflag* determines whether or not to expand the bispectrum
-rows of the global array returned by compute snap. If *bikflag* is set
-to *1* then the bispectrum row, which is typically the per-atom bispectrum
-descriptors :math:`B_{i,k}` summed over all atoms *i* to produce
-:math:`B_k`, becomes bispectrum rows equal to the number of atoms. Thus,
-the resulting bispectrum rows are :math:`B_{i,k}` instead of just
-:math:`B_k`. In this case, the entries in the final column for these rows
-are set to zero.
-
 The keyword *switchinnerflag* with value 1
 activates an additional radial switching
 function similar to :math:`f_c(r)` above, but acting to switch off
@@ -393,6 +387,36 @@ When the central atom and the neighbor atom have different types,
 the values of :math:`S_{inner}` and :math:`D_{inner}` are
 the arithmetic means of the values for both types.
 
+The keywords *bikflag* and *dgradflag* are only used by compute *snap*.
+The keyword *bikflag* determines whether or not to list the descriptors
+of each atom separately, or sum them together and list in a single row.
+If *bikflag* is set
+to *0* then a single bispectrum row is used, which contains the per-atom bispectrum
+descriptors :math:`B_{i,k}` summed over all atoms *i* to produce
+:math:`B_k`.  If *bikflag* is set
+to *1* this is replaced by a separate per-atom bispectrum row for each atom.
+In this case, the entries in the final column for these rows
+are set to zero.
+
+The keyword *dgradflag* determines whether to sum atom gradients or list
+them separately. If *dgradflag* is set to 0, the bispectrum
+descriptor gradients w.r.t. atom *j* are summed over all atoms *i'*
+of type *I* (similar to *snad/atom* above).
+If *dgradflag* is set to 1, gradients are listed separately for each pair of atoms.
+Each row corresponds
+to a single term :math:`\frac{\partial {B_{i,k}  }}{\partial {r}^a_j}`
+where :math:`{r}^a_j` is the *a-th* position coordinate of the atom with global
+index *j*. This also changes
+the number of columns to be equal to the number of bispectrum components, with 3
+additional columns representing the indices :math:`i`, :math:`j`, and :math:`a`,
+as explained more in the Output info section below. The option *dgradflag=1*
+requires that *bikflag=1*.
+
+.. note::
+
+   Using *dgradflag* = 1 produces a global array with :math:`N + 3N^2 + 1` rows
+   which becomes expensive for systems with more than 1000 atoms.
+
 .. note::
 
    If you have a bonded system, then the settings of :doc:`special_bonds
@@ -503,6 +527,42 @@ components. For the purposes of handling contributions to force, virial,
 and quadratic combinations, these :math:`N_{elem}^3` sub-blocks are
 treated as a single block of :math:`K N_{elem}^3` columns.
 
+If the *bik* keyword is set to 1, the structure of the snap array is expanded.
+The first :math:`N` rows of the snap array
+correspond to :math:`B_{i,k}` instead of a single row summed over atoms :math:`i`.
+In this case, the entries in the final column for these rows
+are set to zero. Also, each row contains only non-zero entries for the
+columns corresponding to the type of that atom. This is not true in the case
+of *dgradflag* keyword = 1 (see below).
+
+If the *dgradflag* keyword is set to 1, this changes the structure of the
+global array completely.
+Here the *snad/atom* quantities are replaced with rows corresponding to
+descriptor gradient components on single atoms:
+
+.. math::
+
+  \frac{\partial {B_{i,k}  }}{\partial {r}^a_j}
+
+where :math:`{r}^a_j` is the *a-th* position coordinate of the atom with global
+index *j*. The rows are
+organized in chunks, where each chunk corresponds to an atom with global index
+:math:`j`. The rows in an atom :math:`j` chunk correspond to
+atoms with global index :math:`i`. The total number of rows for
+these descriptor gradients is therefore :math:`3N^2`.
+The number of columns is equal to the number of bispectrum components,
+plus 3 additional left-most columns representing the global atom indices
+:math:`i`, :math:`j`,
+and Cartesian direction :math:`a`  (0, 1, 2, for x, y, z).
+The first 3 columns of the first :math:`N` rows belong to the reference
+potential force components. The remaining K columns contain the
+:math:`B_{i,k}` per-atom descriptors corresponding to the non-zero entries
+obtained when *bikflag* = 1.
+The first column of the last row, after the first
+:math:`N + 3N^2` rows, contains the reference potential
+energy. The virial components are not used with this option. The total number of
+rows is therefore :math:`N + 3N^2 + 1` and the number of columns is :math:`K + 3`.
+
 These values can be accessed by any command that uses per-atom values
 from a compute as input.  See the :doc:`Howto output <Howto_output>` doc
 page for an overview of LAMMPS output options. To see how this command
diff --git a/doc/src/dump.rst b/doc/src/dump.rst
index 4417e186ec..843d2e4883 100644
--- a/doc/src/dump.rst
+++ b/doc/src/dump.rst
@@ -783,6 +783,11 @@ To write gzipped dump files, you must either compile LAMMPS with the
 -DLAMMPS_GZIP option or use the styles from the COMPRESS package.
 See the :doc:`Build settings <Build_settings>` page for details.
 
+While a dump command is active (i.e. has not been stopped by using
+the undump command), no commands may be used that will change the
+timestep (e.g. :doc:`reset_timestep <reset_timestep>`).  LAMMPS
+will terminate with an error otherwise.
+
 The *atom/gz*, *cfg/gz*, *custom/gz*, and *xyz/gz* styles are part of
 the COMPRESS package.  They are only enabled if LAMMPS was built with
 that package.  See the :doc:`Build package <Build_package>` page for
diff --git a/doc/src/fix_reaxff_species.rst b/doc/src/fix_reaxff_species.rst
index a652777ba7..11f0e7b7e7 100644
--- a/doc/src/fix_reaxff_species.rst
+++ b/doc/src/fix_reaxff_species.rst
@@ -20,7 +20,7 @@ Syntax
 * Nfreq = calculate average bond-order every this many timesteps
 * filename = name of output file
 * zero or more keyword/value pairs may be appended
-* keyword = *cutoff* or *element* or *position*
+* keyword = *cutoff* or *element* or *position* or *delete*
 
   .. parsed-literal::
 
@@ -31,6 +31,14 @@ Syntax
        *position* value = posfreq filepos
          posfreq = write position files every this many timestep
          filepos = name of position output file
+       *delete* value = filedel keyword value
+         filedel = name of delete species output file
+         keyword = *specieslist* or *masslimit*
+           *specieslist* value = Nspecies Species1 Species2 ...
+             Nspecies = number of species in list
+           *masslimit* value = massmin massmax
+             massmin = minimum molecular weight of species to delete
+             massmax = maximum molecular weight of species to delete
 
 Examples
 """"""""
@@ -40,6 +48,7 @@ Examples
    fix 1 all reaxff/species 10 10 100 species.out
    fix 1 all reaxff/species 1 2 20 species.out cutoff 1 1 0.40 cutoff 1 2 0.55
    fix 1 all reaxff/species 1 100 100 species.out element Au O H position 1000 AuOH.pos
+   fix 1 all reaxff/species 1 100 100 species.out delete species.del masslimit 0 50
 
 Description
 """""""""""
@@ -59,13 +68,18 @@ the first line.
 .. warning::
 
    In order to compute averaged data, it is required that there are no
-   neighbor list rebuilds between the *Nfreq* steps. For that reason, fix
-   *reaxff/species* may change your neighbor list settings.  There will
-   be a warning message showing the new settings. Having an *Nfreq*
-   setting that is larger than what is required for correct computation
-   of the ReaxFF force field interactions can thus lead to incorrect
-   results.  For typical ReaxFF calculations a value of 100 is already
-   quite large.
+   neighbor list rebuilds for at least Nrepeat\*Nevery steps preceding
+   each *Nfreq* step.  For that reason, fix *reaxff/species* may
+   change your neighbor list settings.  Reneighboring will occur no
+   more frequently than every Nrepeat\*Nevery timesteps, and will
+   occur less frequently if *Nfreq* is not a multiple of
+   Nrepeat\*Nevery.  There will be a warning message showing the new
+   settings.  Having a *Nfreq* setting that is larger than what is
+   required for correct computation of the ReaxFF force field
+   interactions, in combination with certain *Nrepeat* and *Nevery*
+   settings, can thus lead to incorrect results.  For typical ReaxFF
+   calculations, reneighboring only every 100 steps is already quite a
+   low frequency.
 
 If the filename ends with ".gz", the output file is written in gzipped
 format.  A gzipped dump file will be about 3x smaller than the text version,
@@ -104,6 +118,30 @@ character appears in *filepos*, then one file per snapshot is written
 at *posfreq* and the "\*" character is replaced with the timestep
 value.  For example, AuO.pos.\* becomes AuO.pos.0, AuO.pos.1000, etc.
 
+The optional keyword *delete* enables the periodic removal of
+molecules from the system.  Criteria for deletion can be either a list
+of specific chemical formulae or a range of molecular weights.
+Molecules are deleted every *Nfreq* timesteps, and bond connectivity
+is determined using the *Nevery* and *Nrepeat* keywords.  The
+*filedel* argument is the name of the output file that records the
+species that are removed from the system.  The *specieslist* keyword
+permits specific chemical species to be deleted.  The *Nspecies*
+argument specifies how many species are eligible for deletion and is
+followed by a list of chemical formulae, whose strings are compared to
+species identified by this fix.  For example, "specieslist 2 CO CO2"
+deletes molecules that are identified as "CO" and "CO2" in the species
+output file.  When using the *specieslist* keyword, the *filedel* file
+has the following format: the first line lists the chemical formulae
+eligible for deletion, and each additional line contains the timestep
+on which a molecule deletion occurs and the number of each species
+deleted on that timestep.  The *masslimit* keyword permits deletion of
+molecules with molecular weights between *massmin* and *massmax*.
+When using the *masslimit* keyword, each line of the *filedel* file
+contains the timestep on which deletions occurs, followed by how many
+of each species are deleted (with quantities preceding chemical
+formulae).  The *specieslist* and *masslimit* keywords cannot both be
+used in the same *reaxff/species* fix.
+
 ----------
 
 The *Nevery*, *Nrepeat*, and *Nfreq* arguments specify on what
diff --git a/doc/src/pair_sw.rst b/doc/src/pair_sw.rst
index 0a6b284b96..c456b2f422 100644
--- a/doc/src/pair_sw.rst
+++ b/doc/src/pair_sw.rst
@@ -24,13 +24,16 @@ Syntax
    pair_style style keyword values
 
 * style = *sw* or *sw/mod*
-* keyword = *maxdelcs*
+* keyword = *maxdelcs* or *threebody*
 
   .. parsed-literal::
 
-       *maxdelcs* value = delta1 delta2 (optional)
+       *maxdelcs* value = delta1 delta2 (optional, sw/mod only)
          delta1 = The minimum thershold for the variation of cosine of three-body angle
          delta2 = The maximum threshold for the variation of cosine of three-body angle
+       *threebody* value = *on* or *off* (optional, sw only)
+         on (default) = Compute both the three-body and two-body terms of the potential
+         off = Compute only the two-body term of the potential
 
 Examples
 """"""""
@@ -44,6 +47,11 @@ Examples
    pair_style sw/mod maxdelcs 0.25 0.35
    pair_coeff * * tmd.sw.mod Mo S S
 
+   pair_style hybrid sw threebody on sw threebody off
+   pair_coeff * * sw 1 mW_xL.sw mW NULL
+   pair_coeff 1 2 sw 2 mW_xL.sw mW xL
+   pair_coeff 2 2 sw 2 mW_xL.sw mW xL
+
 Description
 """""""""""
 
@@ -68,22 +76,25 @@ three-body term.  The summations in the formula are over all neighbors J
 and K of atom I within a cutoff distance :math:`a `\sigma`.
 
 The *sw/mod* style is designed for simulations of materials when
-distinguishing three-body angles are necessary, such as borophene
-and transition metal dichalcogenides, which cannot be described
-by the original code for the Stillinger-Weber potential.
-For instance, there are several types of angles around each Mo atom in `MoS_2`,
-and some unnecessary angle types should be excluded in the three-body interaction.
-Such exclusion may be realized by selecting proper angle types directly.
-The exclusion of unnecessary angles is achieved here by the cut-off function (`f_C(\delta)`),
-which induces only minimum modifications for LAMMPS.
+distinguishing three-body angles are necessary, such as borophene and
+transition metal dichalcogenides, which cannot be described by the
+original code for the Stillinger-Weber potential.  For instance, there
+are several types of angles around each Mo atom in `MoS_2`, and some
+unnecessary angle types should be excluded in the three-body
+interaction.  Such exclusion may be realized by selecting proper angle
+types directly.  The exclusion of unnecessary angles is achieved here by
+the cut-off function (`f_C(\delta)`), which induces only minimum
+modifications for LAMMPS.
 
 Validation, benchmark tests, and applications of the *sw/mod* style
 can be found in :ref:`(Jiang2) <Jiang2>` and :ref:`(Jiang3) <Jiang3>`.
 
-The *sw/mod* style computes the energy E of a system of atoms, whose potential
-function is mostly the same as the Stillinger-Weber potential. The only modification
-is in the three-body term, where the value of :math:`\delta = \cos \theta_{ijk} - \cos \theta_{0ijk}`
-used in the original energy and force expression is scaled by a switching factor :math:`f_C(\delta)`:
+The *sw/mod* style computes the energy E of a system of atoms, whose
+potential function is mostly the same as the Stillinger-Weber
+potential. The only modification is in the three-body term, where the
+value of :math:`\delta = \cos \theta_{ijk} - \cos \theta_{0ijk}` used in
+the original energy and force expression is scaled by a switching factor
+:math:`f_C(\delta)`:
 
 .. math::
 
@@ -94,28 +105,46 @@ used in the original energy and force expression is scaled by a switching factor
     0 & \left| \delta \right| > \delta_2
     \end{array} \right. \\
 
-This cut-off function decreases smoothly from 1 to 0 over the range :math:`[\delta_1, \delta_2]`.
-This smoothly turns off the energy and force contributions for :math:`\left| \delta \right| > \delta_2`.
-It is suggested that :math:`\delta 1` and :math:`\delta_2` to be the value around
-:math:`0.5 \left| \cos \theta_1 - \cos \theta_2 \right|`, with
-:math:`\theta_1` and :math:`\theta_2` as the different types of angles around an atom.
-For borophene and transition metal dichalcogenides, :math:`\delta_1 = 0.25` and :math:`\delta_2 = 0.35`.
-This value enables the cut-off function to exclude unnecessary angles in the three-body SW terms.
+This cut-off function decreases smoothly from 1 to 0 over the range
+:math:`[\delta_1, \delta_2]`.  This smoothly turns off the energy and
+force contributions for :math:`\left| \delta \right| > \delta_2`.  It is
+suggested that :math:`\delta 1` and :math:`\delta_2` to be the value
+around :math:`0.5 \left| \cos \theta_1 - \cos \theta_2 \right|`, with
+:math:`\theta_1` and :math:`\theta_2` as the different types of angles
+around an atom.  For borophene and transition metal dichalcogenides,
+:math:`\delta_1 = 0.25` and :math:`\delta_2 = 0.35`.  This value enables
+the cut-off function to exclude unnecessary angles in the three-body SW
+terms.
 
 .. note::
 
-   The cut-off function is just to be used as a technique to exclude some unnecessary angles,
-   and it has no physical meaning. It should be noted that the force and potential are inconsistent
-   with each other in the decaying range of the cut-off function, as the angle dependence for the
-   cut-off function is not implemented in the force (first derivation of potential).
-   However, the angle variation is much smaller than the given threshold value for actual simulations,
-   so the inconsistency between potential and force can be neglected in actual simulations.
+   The cut-off function is just to be used as a technique to exclude
+   some unnecessary angles, and it has no physical meaning. It should be
+   noted that the force and potential are inconsistent with each other
+   in the decaying range of the cut-off function, as the angle
+   dependence for the cut-off function is not implemented in the force
+   (first derivation of potential).  However, the angle variation is
+   much smaller than the given threshold value for actual simulations,
+   so the inconsistency between potential and force can be neglected in
+   actual simulations.
 
-Only a single pair_coeff command is used with the *sw* and *sw/mod* styles
-which specifies a Stillinger-Weber potential file with parameters for all
-needed elements.  These are mapped to LAMMPS atom types by specifying
-N additional arguments after the filename in the pair_coeff command,
-where N is the number of LAMMPS atom types:
+The *threebody* keyword is optional and determines whether or not the
+three-body term of the potential is calculated.  The default value is
+"on" and it is only available for the plain *sw* pair style variants,
+but not available for the *sw/mod* and :doc:`sw/angle/table
+<pair_sw_angle_table>` pair style variants.  To turn off the threebody
+contributions all :math:`\lambda_{ijk}` parameters from the potential
+file are forcibly set to 0.  In addition the pair style implementation
+may employ code optimizations for the *threebody off* setting that can
+result in significant speedups versus the default.  These code optimizations
+are currently only available for the MANYBODY and OPENMP packages.
+
+Only a single pair_coeff command is used with the *sw* and *sw/mod*
+styles which specifies a Stillinger-Weber potential file with parameters
+for all needed elements, except for when the *threebody off* setting is
+used (see note below).  These are mapped to LAMMPS atom types by
+specifying N additional arguments after the filename in the pair_coeff
+command, where N is the number of LAMMPS atom types:
 
 * filename
 * N element names = mapping of SW elements to atom types
@@ -130,16 +159,22 @@ pair_coeff command:
 
 .. code-block:: LAMMPS
 
+   pair_style sw
    pair_coeff * * SiC.sw Si Si Si C
 
 The first 2 arguments must be \* \* so as to span all LAMMPS atom types.
-The first three Si arguments map LAMMPS atom types 1,2,3 to the Si
-element in the SW file.  The final C argument maps LAMMPS atom type 4
-to the C element in the SW file.  If a mapping value is specified as
-NULL, the mapping is not performed.  This can be used when a *sw*
+The first three Si arguments map LAMMPS atom types 1, 2, and 3 to the Si
+element in the SW file.  The final C argument maps LAMMPS atom type 4 to
+the C element in the SW file.  If an argument value is specified as
+NULL, the mapping is not performed.  This can be used when an *sw*
 potential is used as part of the *hybrid* pair style.  The NULL values
-are placeholders for atom types that will be used with other
-potentials.
+are placeholders for atom types that will be used with other potentials.
+
+.. note::
+
+   When the *threebody off* keyword is used, multiple pair_coeff commands may
+   be used to specific the pairs of atoms which don't require three-body term.
+   In these cases, the first 2 arguments are not required to be \* \*.
 
 Stillinger-Weber files in the *potentials* directory of the LAMMPS
 distribution have a ".sw" suffix.  Lines that are not blank or
@@ -243,30 +278,39 @@ described above from values in the potential file.
 This pair style does not support the :doc:`pair_modify <pair_modify>`
 shift, table, and tail options.
 
-This pair style does not write its information to :doc:`binary restart files <restart>`, since it is stored in potential files.  Thus, you
-need to re-specify the pair_style and pair_coeff commands in an input
-script that reads a restart file.
+This pair style does not write its information to :doc:`binary restart
+files <restart>`, since it is stored in potential files.  Thus, you need
+to re-specify the pair_style and pair_coeff commands in an input script
+that reads a restart file.
 
 This pair 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.
 
+The single() function of the *sw* pair style is only enabled and
+supported for the case of the *threebody off* setting.
+
 ----------
 
 Restrictions
 """"""""""""
 
-This pair style is part of the MANYBODY package.  It is only enabled
-if LAMMPS was built with that package.  See the :doc:`Build package <Build_package>` page for more info.
+This pair style is part of the MANYBODY package.  It is only enabled if
+LAMMPS was built with that package.  See the :doc:`Build package
+<Build_package>` page for more info.
 
 This pair style requires the :doc:`newton <newton>` setting to be "on"
 for pair interactions.
 
 The Stillinger-Weber potential files provided with LAMMPS (see the
 potentials directory) are parameterized for metal :doc:`units <units>`.
-You can use the SW potential with any LAMMPS units, but you would need
-to create your own SW potential file with coefficients listed in the
-appropriate units if your simulation does not use "metal" units.
+You can use the sw or sw/mod pair styles with any LAMMPS units, but you
+would need to create your own SW potential file with coefficients listed
+in the appropriate units if your simulation does not use "metal" units.
+If the potential file contains a 'UNITS:' metadata tag in the first line
+of the potential file, then LAMMPS can convert it transparently between
+"metal" and "real" units.
+
 
 Related commands
 """"""""""""""""
@@ -276,8 +320,9 @@ Related commands
 Default
 """""""
 
-The default values for the *maxdelcs* setting of the *sw/mod* pair
-style are *delta1* = 0.25 and *delta2* = 0.35`.
+The default value for the *threebody* setting of the "sw" pair style is
+"on", the default values for the "*maxdelcs* setting of the *sw/mod*
+pair style are *delta1* = 0.25 and *delta2* = 0.35`.
 
 ----------
 
diff --git a/doc/src/pair_sw_angle_table.rst b/doc/src/pair_sw_angle_table.rst
index ff88711d7a..65e73d4445 100644
--- a/doc/src/pair_sw_angle_table.rst
+++ b/doc/src/pair_sw_angle_table.rst
@@ -296,7 +296,8 @@ for pair interactions.
 Related commands
 """"""""""""""""
 
-:doc:`pair_coeff <pair_coeff>`, :doc:`pair_style sw <pair_sw>`, :doc:`pair_style threebody/table <pair_threebody_table>`
+:doc:`pair_coeff <pair_coeff>`, :doc:`pair_style sw <pair_sw>`,
+:doc:`pair_style threebody/table <pair_threebody_table>`
 
 
 ----------
diff --git a/doc/src/read_dump.rst b/doc/src/read_dump.rst
index 3a771b9c2d..c98347914b 100644
--- a/doc/src/read_dump.rst
+++ b/doc/src/read_dump.rst
@@ -24,12 +24,13 @@ Syntax
        *fx*,\ *fy*,\ *fz* = force components
 
 * zero or more keyword/value pairs may be appended
-* keyword = *nfile* or *box* or *replace* or *purge* or *trim* or *add* or *label* or *scaled* or *wrapped* or *format*
+* keyword = *nfile* or *box* or *timestep* or *replace* or *purge* or *trim* or *add* or *label* or *scaled* or *wrapped* or *format*
 
   .. parsed-literal::
 
        *nfile* value = Nfiles = how many parallel dump files exist
        *box* value = *yes* or *no* = replace simulation box with dump box
+       *timestep* value = *yes* or *no* = reset simulation timestep with dump timestep
        *replace* value = *yes* or *no* = overwrite atoms with dump atoms
        *purge* value = *yes* or *no* = delete all atoms before adding dump atoms
        *trim* value = *yes* or *no* = trim atoms not in dump snapshot
@@ -60,6 +61,7 @@ Examples
    read_dump dump.dcd 0 x y z box yes format molfile dcd
    read_dump dump.file 1000 x y z vx vy vz box yes format molfile lammpstrj /usr/local/lib/vmd/plugins/LINUXAMD64/plugins/molfile
    read_dump dump.file 5000 x y vx vy trim yes
+   read_dump dump.file 5000 x y vx vy add yes box no timestep no
    read_dump ../run7/dump.file.gz 10000 x y z box yes
    read_dump dump.xyz 10 x y z box no format molfile xyz ../plugins
    read_dump dump.dcd 0 x y z format molfile dcd
@@ -71,9 +73,9 @@ Description
 """""""""""
 
 Read atom information from a dump file to overwrite the current atom
-coordinates, and optionally the atom velocities and image flags and
-the simulation box dimensions.  This is useful for restarting a run
-from a particular snapshot in a dump file.  See the
+coordinates, and optionally the atom velocities and image flags, the
+simulation timestep, and the simulation box dimensions.  This is useful
+for restarting a run from a particular snapshot in a dump file.  See the
 :doc:`read_restart <read_restart>` and :doc:`read_data <read_data>`
 commands for alternative methods to do this.  Also see the
 :doc:`rerun <rerun>` command for a means of reading multiple snapshots
@@ -89,9 +91,9 @@ Also note that reading per-atom information from a dump snapshot is
 limited to the atom coordinates, velocities and image flags, as
 explained below.  Other atom properties, which may be necessary to run
 a valid simulation, such as atom charge, or bond topology information
-for a molecular system, are not read from (or even contained in) dump
-files.  Thus this auxiliary information should be defined in the usual
-way, e.g. in a data file read in by a :doc:`read_data <read_data>`
+for a molecular system, are not read from (or may not even be contained
+in) dump files.  Thus this auxiliary information should be defined in
+the usual way, e.g. in a data file read in by a :doc:`read_data <read_data>`
 command, before using the read_dump command, or by the :doc:`set <set>`
 command, after the dump snapshot is read.
 
@@ -165,11 +167,10 @@ variable *ntimestep*:
     uint64_t  ntimestep  5*scalar
       (0)    0 50 100 150 200
 
-Note that the *xyz*
-and *molfile* formats do not store the timestep.  For these formats,
-timesteps are numbered logically, in a sequential manner, starting
-from 0.  Thus to access the 10th snapshot in an *xyz* or *mofile*
-formatted dump file, use *Nstep* = 9.
+Note that the *xyz* and *molfile* formats do not store the timestep.
+For these formats, timesteps are numbered logically, in a sequential
+manner, starting from 0.  Thus to access the 10th snapshot in an *xyz*
+or *mofile* formatted dump file, use *Nstep* = 9.
 
 The dimensions of the simulation box for the selected snapshot are
 also read; see the *box* keyword discussion below.  For the *native*
@@ -266,8 +267,10 @@ for how this is done, determined by the specified fields and optional
 keywords.
 
 The timestep of the snapshot becomes the current timestep for the
-simulation.  See the :doc:`reset_timestep <reset_timestep>` command if
-you wish to change this after the dump snapshot is read.
+simulation unless the *timestep* keyword is specified with a *no* value
+(default setting is *yes*).  See the :doc:`reset_timestep <reset_timestep>`
+command if you wish to change this to a different value after the dump
+snapshot is read.
 
 If the *box* keyword is specified with a *yes* value, then the current
 simulation box dimensions are replaced by the dump snapshot box
@@ -391,7 +394,7 @@ Related commands
 Default
 """""""
 
-The option defaults are box = yes, replace = yes, purge = no, trim =
-no, add = no, scaled = no, wrapped = yes, and format = native.
+The option defaults are box = yes, timestep = yes, replace = yes, purge = no,
+trim = no, add = no, scaled = no, wrapped = yes, and format = native.
 
 .. _vmd: http://www.ks.uiuc.edu/Research/vmd
diff --git a/examples/snap/README.md b/examples/snap/README.md
new file mode 100644
index 0000000000..305f920ae8
--- /dev/null
+++ b/examples/snap/README.md
@@ -0,0 +1,13 @@
+This directory contains a variety of tests for the ML-SNAP package. These include:
+
+in.snap.Ta06A                     # SNAP linear Ta potential
+in.snap.W.2940                    # SNAP linear W potential
+in.snap.hybrid.WSNAP.HePair       # Hybrid overlay pair style for linear SNAP W potential and twobody tables for He-He and W-He 
+in.snap.WBe.PRB2019               # SNAP linear W/Be potential
+in.snap.InP.JCPA2020              # SNAP linear InP potential using chem keyword (explicit multi-element)
+in.snap.Mo_Chen                   # SNAP linear Mo potential
+in.snap.compute                   # SNAP compute for training a linear model
+in.snap.compute.quadratic         # SNAP compute for training a quadratic model
+in.snap.scale.Ni_Zuo_JCPA2020     # SNAP linear Ni potential with thermodynamic integration (fix adapt scale)
+
+compute_snap_dgrad.py             # SNAP compute with dgradflag (dBi/dRj) for training a non-linear model
diff --git a/examples/snap/compute_snap_dgrad.py b/examples/snap/compute_snap_dgrad.py
new file mode 100644
index 0000000000..0919af8598
--- /dev/null
+++ b/examples/snap/compute_snap_dgrad.py
@@ -0,0 +1,173 @@
+"""
+compute_snap_dgrad.py
+Purpose: Demonstrate extraction of descriptor gradient (dB/dR) array from compute snap.
+         Show that dBi/dRj components summed over neighbors i yields same output as regular compute snap with dgradflag = 0.
+         This shows that the dBi/dRj components extracted with dgradflag = 1 are correct.
+Serial syntax:
+    python compute_snap_dgrad.py
+Parallel syntax:
+    mpirun -np 4 python compute_snap_dgrad.py
+"""
+
+from __future__ import print_function
+import sys
+import ctypes
+import numpy as np
+from lammps import lammps, LMP_TYPE_ARRAY, LMP_STYLE_GLOBAL
+
+# get MPI settings from LAMMPS
+
+lmp = lammps()
+me = lmp.extract_setting("world_rank")
+nprocs = lmp.extract_setting("world_size")
+
+cmds = ["-screen", "none", "-log", "none"]
+lmp = lammps(cmdargs = cmds)
+
+def run_lammps(dgradflag):
+
+    # simulation settings
+
+    lmp.command("clear")
+    lmp.command("units metal")
+    lmp.command("boundary	p p p")
+    lmp.command("atom_modify	map hash")
+    lmp.command(f"lattice         bcc {latparam}")
+    lmp.command(f"region		box block 0 {nx} 0 {ny} 0 {nz}")
+    lmp.command(f"create_box	{ntypes} box")
+    lmp.command(f"create_atoms	{ntypes} box")
+    lmp.command("mass 		* 180.88")
+    lmp.command("displace_atoms 	all random 0.01 0.01 0.01 123456")
+
+    # potential settings
+
+    snap_options = f'{rcutfac} {rfac0} {twojmax} {radelem1} {radelem2} {wj1} {wj2} rmin0 {rmin0} quadraticflag {quadratic} bzeroflag {bzero} switchflag {switch} bikflag {bikflag} dgradflag {dgradflag}'
+    lmp.command(f"pair_style 	zero {rcutfac}")
+    lmp.command(f"pair_coeff 	* *")
+    lmp.command(f"pair_style 	zbl {zblcutinner} {zblcutouter}")
+    lmp.command(f"pair_coeff 	* * {zblz} {zblz}")
+
+    # define compute snap
+
+    lmp.command(f"compute 	snap all snap {snap_options}")
+
+    # run
+
+    lmp.command(f"thermo 		100")
+    lmp.command(f"run {nsteps}")
+
+# declare simulation/structure variables
+
+nsteps = 0
+nrep = 2
+latparam = 2.0
+ntypes = 2
+nx = nrep
+ny = nrep
+nz = nrep
+
+# declare compute snap variables
+
+twojmax = 8
+rcutfac = 1.0
+rfac0 = 0.99363
+rmin0 = 0
+radelem1 = 2.3
+radelem2 = 2.0
+wj1 = 1.0
+wj2 = 0.96
+quadratic = 0
+bzero = 0
+switch = 0
+bikflag = 1
+
+# define reference potential
+
+zblcutinner = 4.0
+zblcutouter = 4.8
+zblz = 73
+
+# number of descriptors
+
+if (twojmax % 2 == 0):
+    m  =  twojmax / 2 + 1
+    nd = int(m*(m+1)*(2*m+1)/6)
+else:
+    m  =  (twojmax + 1) / 2
+    nd = int(m*(m+1)*(m+2)/3)
+
+if me == 0:
+    print(f"Number of descriptors based on twojmax : {nd}")
+
+# run lammps with dgradflag on
+
+if me == 0:
+    print("Running with dgradflag on")
+
+dgradflag = 1
+run_lammps(dgradflag)
+
+# get global snap array
+
+lmp_snap = lmp.numpy.extract_compute("snap", LMP_STYLE_GLOBAL, LMP_TYPE_ARRAY)
+
+# print snap array to observe
+#if (me==0):
+#    np.savetxt("test_snap.dat", lmp_snap, fmt="%d %d %d %f %f %f %f %f")
+
+# take out rows with zero column
+
+# extract dBj/dRi (includes dBi/dRi)
+
+natoms = lmp.get_natoms()
+fref1 = lmp_snap[0:natoms,0:3].flatten()
+eref1 = lmp_snap[-1,0]
+dbdr_length = np.shape(lmp_snap)[0]-(natoms) - 1
+dBdR = lmp_snap[natoms:(natoms+dbdr_length),3:(nd+3)]
+force_indices = lmp_snap[natoms:(natoms+dbdr_length),0:3].astype(np.int32)
+
+# strip rows with all zero descriptor gradients to demonstrate how to save memory
+
+nonzero_rows = lmp_snap[natoms:(natoms+dbdr_length),3:(nd+3)] != 0.0
+nonzero_rows = np.any(nonzero_rows, axis=1)
+dBdR = dBdR[nonzero_rows, :]
+force_indices = force_indices[nonzero_rows,:]
+dbdr_length = np.shape(dBdR)[0]
+
+# sum over atoms i that j is a neighbor of, like dgradflag = 0 does.
+
+array1 = np.zeros((3*natoms,nd))
+for k in range(0,nd):
+    for l in range(0,dbdr_length):
+        i = force_indices[l,0]
+        j = force_indices[l,1]
+        a = force_indices[l,2]
+        #print(f"{i} {j} {a}")
+        array1[3 * j + a, k] += dBdR[l,k]
+
+# run lammps with dgradflag off
+
+if me == 0:
+    print("Running with dgradflag off")
+
+dgradflag = 0
+run_lammps(dgradflag)
+
+# get global snap array
+
+lmp_snap = lmp.numpy.extract_compute("snap", LMP_STYLE_GLOBAL, LMP_TYPE_ARRAY)
+natoms = lmp.get_natoms()
+fref2 = lmp_snap[natoms:(natoms+3*natoms),-1]
+eref2 = lmp_snap[0,-1]
+array2 = lmp_snap[natoms:natoms+(3*natoms), nd:-1]
+
+# take difference of arrays obtained from dgradflag on and off.
+
+diffm = array1 - array2
+difff = fref1 - fref2
+diffe = eref1 - eref2
+
+if me == 0:
+    print(f"Max/min difference in dSum(Bi)/dRj: {np.max(diffm)} {np.min(diffm)}")
+    print(f"Max/min difference in reference forces: {np.max(difff)} {np.min(difff)}")
+    print(f"Difference in reference energy: {diffe}")
diff --git a/lib/kokkos/CHANGELOG.md b/lib/kokkos/CHANGELOG.md
index dfbe22edde..a908507704 100644
--- a/lib/kokkos/CHANGELOG.md
+++ b/lib/kokkos/CHANGELOG.md
@@ -1,5 +1,21 @@
 # Change Log
 
+## [3.6.01](https://github.com/kokkos/kokkos/tree/3.6.01) (2022-05-23)
+[Full Changelog](https://github.com/kokkos/kokkos/compare/3.6.00...3.6.01)
+
+### Bug Fixes:
+- Fix Threads: Fix serial resizing scratch space (3.6.01 cherry-pick) [\#5109](https://github.com/kokkos/kokkos/pull/5109)
+- Fix ScatterMin/ScatterMax to use proper atomics (3.6.01 cherry-pick) [\#5046](https://github.com/kokkos/kokkos/pull/5046)
+- Fix allocating large Views [\#4907](https://github.com/kokkos/kokkos/pull/4907)
+- Fix bounds errors with Kokkos::sort [\#4980](https://github.com/kokkos/kokkos/pull/4980)
+- Fix HIP version when printing the configuration [\#4872](https://github.com/kokkos/kokkos/pull/4872)
+- Fixed `_CUDA_ARCH__` to `__CUDA_ARCH__` for CUDA LDG [\#4893](https://github.com/kokkos/kokkos/pull/4893)
+- Fixed an incorrect struct initialization [\#5028](https://github.com/kokkos/kokkos/pull/5028)
+- Fix racing condition in `HIPParallelLaunch` [\#5008](https://github.com/kokkos/kokkos/pull/5008)
+- Avoid deprecation warnings with `OpenMPExec::validate_partition` [\#4982](https://github.com/kokkos/kokkos/pull/4982)
+- Make View self-assignment not produce double-free [\#5024](https://github.com/kokkos/kokkos/pull/5024)
+
+
 ## [3.6.00](https://github.com/kokkos/kokkos/tree/3.6.00) (2022-02-18)
 [Full Changelog](https://github.com/kokkos/kokkos/compare/3.5.00...3.6.00)
 
diff --git a/lib/kokkos/CMakeLists.txt b/lib/kokkos/CMakeLists.txt
index e1c6893725..b0a54118a0 100644
--- a/lib/kokkos/CMakeLists.txt
+++ b/lib/kokkos/CMakeLists.txt
@@ -136,7 +136,7 @@ ENDIF()
 
 set(Kokkos_VERSION_MAJOR 3)
 set(Kokkos_VERSION_MINOR 6)
-set(Kokkos_VERSION_PATCH 00)
+set(Kokkos_VERSION_PATCH 01)
 set(Kokkos_VERSION "${Kokkos_VERSION_MAJOR}.${Kokkos_VERSION_MINOR}.${Kokkos_VERSION_PATCH}")
 math(EXPR KOKKOS_VERSION "${Kokkos_VERSION_MAJOR} * 10000 + ${Kokkos_VERSION_MINOR} * 100 + ${Kokkos_VERSION_PATCH}")
 
diff --git a/lib/kokkos/Makefile.kokkos b/lib/kokkos/Makefile.kokkos
index aa5f7c98f8..755831452b 100644
--- a/lib/kokkos/Makefile.kokkos
+++ b/lib/kokkos/Makefile.kokkos
@@ -12,7 +12,7 @@ endif
 
 KOKKOS_VERSION_MAJOR = 3
 KOKKOS_VERSION_MINOR = 6
-KOKKOS_VERSION_PATCH = 00
+KOKKOS_VERSION_PATCH = 01
 KOKKOS_VERSION = $(shell echo $(KOKKOS_VERSION_MAJOR)*10000+$(KOKKOS_VERSION_MINOR)*100+$(KOKKOS_VERSION_PATCH) | bc)
 
 # Options: Cuda,HIP,SYCL,OpenMPTarget,OpenMP,Threads,Serial
diff --git a/lib/kokkos/algorithms/src/Kokkos_Sort.hpp b/lib/kokkos/algorithms/src/Kokkos_Sort.hpp
index cde5e6857e..ce97de9b7d 100644
--- a/lib/kokkos/algorithms/src/Kokkos_Sort.hpp
+++ b/lib/kokkos/algorithms/src/Kokkos_Sort.hpp
@@ -422,54 +422,34 @@ class BinSort {
 
 template <class KeyViewType>
 struct BinOp1D {
-  int max_bins_;
-  double mul_;
-  typename KeyViewType::const_value_type range_;
-  typename KeyViewType::const_value_type min_;
+  int max_bins_ = {};
+  double mul_   = {};
+  double min_   = {};
 
-  BinOp1D()
-      : max_bins_(0),
-        mul_(0.0),
-        range_(typename KeyViewType::const_value_type()),
-        min_(typename KeyViewType::const_value_type()) {}
+  BinOp1D() = default;
 
   // Construct BinOp with number of bins, minimum value and maxuimum value
   BinOp1D(int max_bins__, typename KeyViewType::const_value_type min,
           typename KeyViewType::const_value_type max)
       : max_bins_(max_bins__ + 1),
-        // Cast to int64_t to avoid possible overflow when using integer
-        mul_(std::is_integral<typename KeyViewType::const_value_type>::value
-                 ? 1.0 * max_bins__ / (int64_t(max) - int64_t(min))
-                 : 1.0 * max_bins__ / (max - min)),
-        range_(max - min),
-        min_(min) {
+        // Cast to double to avoid possible overflow when using integer
+        mul_(static_cast<double>(max_bins__) /
+             (static_cast<double>(max) - static_cast<double>(min))),
+        min_(static_cast<double>(min)) {
     // For integral types the number of bins may be larger than the range
     // in which case we can exactly have one unique value per bin
     // and then don't need to sort bins.
     if (std::is_integral<typename KeyViewType::const_value_type>::value &&
-        static_cast<uint64_t>(range_) <= static_cast<uint64_t>(max_bins__)) {
+        (static_cast<double>(max) - static_cast<double>(min)) <=
+            static_cast<double>(max_bins__)) {
       mul_ = 1.;
     }
   }
 
   // Determine bin index from key value
-  template <
-      class ViewType,
-      std::enable_if_t<!std::is_integral<typename ViewType::value_type>::value,
-                       bool> = true>
+  template <class ViewType>
   KOKKOS_INLINE_FUNCTION int bin(ViewType& keys, const int& i) const {
-    return int(mul_ * (keys(i) - min_));
-  }
-
-  // Determine bin index from key value
-  template <
-      class ViewType,
-      std::enable_if_t<std::is_integral<typename ViewType::value_type>::value,
-                       bool> = true>
-  KOKKOS_INLINE_FUNCTION int bin(ViewType& keys, const int& i) const {
-    // The cast to int64_t is necessary because otherwise HIP returns the wrong
-    // result.
-    return int(mul_ * (int64_t(keys(i)) - int64_t(min_)));
+    return static_cast<int>(mul_ * (static_cast<double>(keys(i)) - min_));
   }
 
   // Return maximum bin index + 1
@@ -486,10 +466,9 @@ struct BinOp1D {
 
 template <class KeyViewType>
 struct BinOp3D {
-  int max_bins_[3];
-  double mul_[3];
-  typename KeyViewType::non_const_value_type range_[3];
-  typename KeyViewType::non_const_value_type min_[3];
+  int max_bins_[3] = {};
+  double mul_[3]   = {};
+  double min_[3]   = {};
 
   BinOp3D() = default;
 
@@ -498,15 +477,15 @@ struct BinOp3D {
     max_bins_[0] = max_bins__[0];
     max_bins_[1] = max_bins__[1];
     max_bins_[2] = max_bins__[2];
-    mul_[0]      = 1.0 * max_bins__[0] / (max[0] - min[0]);
-    mul_[1]      = 1.0 * max_bins__[1] / (max[1] - min[1]);
-    mul_[2]      = 1.0 * max_bins__[2] / (max[2] - min[2]);
-    range_[0]    = max[0] - min[0];
-    range_[1]    = max[1] - min[1];
-    range_[2]    = max[2] - min[2];
-    min_[0]      = min[0];
-    min_[1]      = min[1];
-    min_[2]      = min[2];
+    mul_[0]      = static_cast<double>(max_bins__[0]) /
+              (static_cast<double>(max[0]) - static_cast<double>(min[0]));
+    mul_[1] = static_cast<double>(max_bins__[1]) /
+              (static_cast<double>(max[1]) - static_cast<double>(min[1]));
+    mul_[2] = static_cast<double>(max_bins__[2]) /
+              (static_cast<double>(max[2]) - static_cast<double>(min[2]));
+    min_[0] = static_cast<double>(min[0]);
+    min_[1] = static_cast<double>(min[1]);
+    min_[2] = static_cast<double>(min[2]);
   }
 
   template <class ViewType>
@@ -596,9 +575,9 @@ std::enable_if_t<Kokkos::is_execution_space<ExecutionSpace>::value> sort(
   // TODO: figure out better max_bins then this ...
   int64_t max_bins = view.extent(0) / 2;
   if (std::is_integral<typename ViewType::non_const_value_type>::value) {
-    // Cast to int64_t to avoid possible overflow when using integer
-    int64_t const max_val = result.max_val;
-    int64_t const min_val = result.min_val;
+    // Cast to double to avoid possible overflow when using integer
+    auto const max_val = static_cast<double>(result.max_val);
+    auto const min_val = static_cast<double>(result.min_val);
     // using 10M as the cutoff for special behavior (roughly 40MB for the count
     // array)
     if ((max_val - min_val) < 10000000) {
@@ -606,6 +585,10 @@ std::enable_if_t<Kokkos::is_execution_space<ExecutionSpace>::value> sort(
       sort_in_bins = false;
     }
   }
+  if (std::is_floating_point<typename ViewType::non_const_value_type>::value) {
+    KOKKOS_ASSERT(std::isfinite(static_cast<double>(result.max_val) -
+                                static_cast<double>(result.min_val)));
+  }
 
   BinSort<ViewType, CompType> bin_sort(
       view, CompType(max_bins, result.min_val, result.max_val), sort_in_bins);
diff --git a/lib/kokkos/algorithms/unit_tests/TestSort.hpp b/lib/kokkos/algorithms/unit_tests/TestSort.hpp
index a03847f2b2..9108731c15 100644
--- a/lib/kokkos/algorithms/unit_tests/TestSort.hpp
+++ b/lib/kokkos/algorithms/unit_tests/TestSort.hpp
@@ -353,6 +353,55 @@ void test_issue_1160_impl() {
   }
 }
 
+template <class ExecutionSpace>
+void test_issue_4978_impl() {
+  Kokkos::View<long long*, ExecutionSpace> element_("element", 9);
+
+  auto h_element = Kokkos::create_mirror_view(element_);
+
+  h_element(0) = LLONG_MIN;
+  h_element(1) = 0;
+  h_element(2) = 3;
+  h_element(3) = 2;
+  h_element(4) = 1;
+  h_element(5) = 3;
+  h_element(6) = 6;
+  h_element(7) = 4;
+  h_element(8) = 3;
+
+  ExecutionSpace exec;
+  Kokkos::deep_copy(exec, element_, h_element);
+
+  Kokkos::sort(exec, element_);
+
+  Kokkos::deep_copy(exec, h_element, element_);
+  exec.fence();
+
+  ASSERT_EQ(h_element(0), LLONG_MIN);
+  ASSERT_EQ(h_element(1), 0);
+  ASSERT_EQ(h_element(2), 1);
+  ASSERT_EQ(h_element(3), 2);
+  ASSERT_EQ(h_element(4), 3);
+  ASSERT_EQ(h_element(5), 3);
+  ASSERT_EQ(h_element(6), 3);
+  ASSERT_EQ(h_element(7), 4);
+  ASSERT_EQ(h_element(8), 6);
+}
+
+template <class ExecutionSpace, class T>
+void test_sort_integer_overflow() {
+  // array with two extrema in reverse order to expose integer overflow bug in
+  // bin calculation
+  T a[2]  = {Kokkos::Experimental::finite_max<T>::value,
+            Kokkos::Experimental::finite_min<T>::value};
+  auto vd = Kokkos::create_mirror_view_and_copy(
+      ExecutionSpace(), Kokkos::View<T[2], Kokkos::HostSpace>(a));
+  Kokkos::sort(vd, /*force using Kokkos bin sort*/ true);
+  auto vh = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), vd);
+  EXPECT_TRUE(std::is_sorted(vh.data(), vh.data() + 2))
+      << "view (" << vh[0] << ", " << vh[1] << ") is not sorted";
+}
+
 //----------------------------------------------------------------------------
 
 template <class ExecutionSpace, typename KeyType>
@@ -376,6 +425,11 @@ void test_issue_1160_sort() {
   test_issue_1160_impl<ExecutionSpace>();
 }
 
+template <class ExecutionSpace>
+void test_issue_4978_sort() {
+  test_issue_4978_impl<ExecutionSpace>();
+}
+
 template <class ExecutionSpace, typename KeyType>
 void test_sort(unsigned int N) {
   test_1D_sort<ExecutionSpace, KeyType>(N);
@@ -385,6 +439,10 @@ void test_sort(unsigned int N) {
   test_dynamic_view_sort<ExecutionSpace, KeyType>(N);
 #endif
   test_issue_1160_sort<ExecutionSpace>();
+  test_issue_4978_sort<ExecutionSpace>();
+  test_sort_integer_overflow<ExecutionSpace, long long>();
+  test_sort_integer_overflow<ExecutionSpace, unsigned long long>();
+  test_sort_integer_overflow<ExecutionSpace, int>();
 }
 }  // namespace Impl
 }  // namespace Test
diff --git a/lib/kokkos/containers/src/Kokkos_ScatterView.hpp b/lib/kokkos/containers/src/Kokkos_ScatterView.hpp
index 024b4618a4..e4dd9531fc 100644
--- a/lib/kokkos/containers/src/Kokkos_ScatterView.hpp
+++ b/lib/kokkos/containers/src/Kokkos_ScatterView.hpp
@@ -369,18 +369,6 @@ struct ScatterValue<ValueType, Kokkos::Experimental::ScatterProd, DeviceType,
     Kokkos::atomic_div(&value, rhs);
   }
 
-  KOKKOS_FORCEINLINE_FUNCTION
-  void atomic_prod(ValueType& dest, const ValueType& src) const {
-    bool success = false;
-    while (!success) {
-      ValueType dest_old = dest;
-      ValueType dest_new = dest_old * src;
-      dest_new =
-          Kokkos::atomic_compare_exchange<ValueType>(&dest, dest_old, dest_new);
-      success = ((dest_new - dest_old) / dest_old <= 1e-15);
-    }
-  }
-
   KOKKOS_INLINE_FUNCTION
   void join(ValueType& dest, const ValueType& src) const {
     atomic_prod(&dest, src);
@@ -440,21 +428,9 @@ struct ScatterValue<ValueType, Kokkos::Experimental::ScatterMin, DeviceType,
   KOKKOS_FORCEINLINE_FUNCTION ScatterValue(ScatterValue&& other)
       : value(other.value) {}
 
-  KOKKOS_FORCEINLINE_FUNCTION
-  void atomic_min(ValueType& dest, const ValueType& src) const {
-    bool success = false;
-    while (!success) {
-      ValueType dest_old = dest;
-      ValueType dest_new = (dest_old > src) ? src : dest_old;
-      dest_new =
-          Kokkos::atomic_compare_exchange<ValueType>(&dest, dest_old, dest_new);
-      success = ((dest_new - dest_old) / dest_old <= 1e-15);
-    }
-  }
-
   KOKKOS_INLINE_FUNCTION
   void join(ValueType& dest, const ValueType& src) const {
-    atomic_min(dest, src);
+    atomic_min(&dest, src);
   }
 
   KOKKOS_INLINE_FUNCTION
@@ -511,21 +487,9 @@ struct ScatterValue<ValueType, Kokkos::Experimental::ScatterMax, DeviceType,
   KOKKOS_FORCEINLINE_FUNCTION ScatterValue(ScatterValue&& other)
       : value(other.value) {}
 
-  KOKKOS_FORCEINLINE_FUNCTION
-  void atomic_max(ValueType& dest, const ValueType& src) const {
-    bool success = false;
-    while (!success) {
-      ValueType dest_old = dest;
-      ValueType dest_new = (dest_old < src) ? src : dest_old;
-      dest_new =
-          Kokkos::atomic_compare_exchange<ValueType>(&dest, dest_old, dest_new);
-      success = ((dest_new - dest_old) / dest_old <= 1e-15);
-    }
-  }
-
   KOKKOS_INLINE_FUNCTION
   void join(ValueType& dest, const ValueType& src) const {
-    atomic_max(dest, src);
+    atomic_max(&dest, src);
   }
 
   KOKKOS_INLINE_FUNCTION
diff --git a/lib/kokkos/containers/src/Kokkos_Vector.hpp b/lib/kokkos/containers/src/Kokkos_Vector.hpp
index 88721bd89e..eddb878003 100644
--- a/lib/kokkos/containers/src/Kokkos_Vector.hpp
+++ b/lib/kokkos/containers/src/Kokkos_Vector.hpp
@@ -162,7 +162,7 @@ class vector : public DualView<Scalar*, LayoutLeft, Arg1Type> {
     }
     DV::sync_host();
     DV::modify_host();
-    if (it < begin() || it > end())
+    if (std::less<>()(it, begin()) || std::less<>()(end(), it))
       Kokkos::abort("Kokkos::vector::insert : invalid insert iterator");
     if (count == 0) return it;
     ptrdiff_t start = std::distance(begin(), it);
@@ -189,27 +189,21 @@ class vector : public DualView<Scalar*, LayoutLeft, Arg1Type> {
                           iterator>::type
   insert(iterator it, InputIterator b, InputIterator e) {
     ptrdiff_t count = std::distance(b, e);
-    if (count == 0) return it;
 
     DV::sync_host();
     DV::modify_host();
-    if (it < begin() || it > end())
+    if (std::less<>()(it, begin()) || std::less<>()(end(), it))
       Kokkos::abort("Kokkos::vector::insert : invalid insert iterator");
 
-    bool resized = false;
-    if ((size() == 0) && (it == begin())) {
-      resize(count);
-      it      = begin();
-      resized = true;
-    }
     ptrdiff_t start = std::distance(begin(), it);
     auto org_size   = size();
-    if (!resized) resize(size() + count);
-    it = begin() + start;
+
+    // Note: resize(...) invalidates it; use begin() + start instead
+    resize(size() + count);
 
     std::copy_backward(begin() + start, begin() + org_size,
                        begin() + org_size + count);
-    std::copy(b, e, it);
+    std::copy(b, e, begin() + start);
 
     return begin() + start;
   }
diff --git a/lib/kokkos/containers/unit_tests/TestVector.hpp b/lib/kokkos/containers/unit_tests/TestVector.hpp
index 57b92c38f8..c093c7b0c9 100644
--- a/lib/kokkos/containers/unit_tests/TestVector.hpp
+++ b/lib/kokkos/containers/unit_tests/TestVector.hpp
@@ -172,6 +172,23 @@ struct test_vector_insert {
       run_test(a);
       check_test(a, size);
     }
+    { test_vector_insert_into_empty(size); }
+  }
+
+  void test_vector_insert_into_empty(const size_t size) {
+    using Vector = Kokkos::vector<Scalar, Device>;
+    {
+      Vector a;
+      Vector b(size);
+      a.insert(a.begin(), b.begin(), b.end());
+      ASSERT_EQ(a.size(), size);
+    }
+
+    {
+      Vector c;
+      c.insert(c.begin(), size, Scalar{});
+      ASSERT_EQ(c.size(), size);
+    }
   }
 };
 
diff --git a/lib/kokkos/core/src/CMakeLists.txt b/lib/kokkos/core/src/CMakeLists.txt
index 88cca93f3c..793e07a841 100644
--- a/lib/kokkos/core/src/CMakeLists.txt
+++ b/lib/kokkos/core/src/CMakeLists.txt
@@ -8,6 +8,7 @@ KOKKOS_INCLUDE_DIRECTORIES(
 INSTALL (DIRECTORY
   "${CMAKE_CURRENT_SOURCE_DIR}/"
   DESTINATION ${KOKKOS_HEADER_DIR}
+  FILES_MATCHING
   PATTERN desul/src EXCLUDE
   PATTERN "*.inc"
   PATTERN "*.inc_*"
diff --git a/lib/kokkos/core/src/Cuda/Kokkos_Cuda_Instance.cpp b/lib/kokkos/core/src/Cuda/Kokkos_Cuda_Instance.cpp
index 294be2774b..aaa9ea8ad4 100644
--- a/lib/kokkos/core/src/Cuda/Kokkos_Cuda_Instance.cpp
+++ b/lib/kokkos/core/src/Cuda/Kokkos_Cuda_Instance.cpp
@@ -1007,6 +1007,15 @@ void CudaSpaceInitializer::print_configuration(std::ostream &msg,
 }
 
 }  // namespace Impl
+
+#ifdef KOKKOS_ENABLE_CXX14
+namespace Tools {
+namespace Experimental {
+constexpr DeviceType DeviceTypeTraits<Cuda>::id;
+}
+}  // namespace Tools
+#endif
+
 }  // namespace Kokkos
 
 #else
diff --git a/lib/kokkos/core/src/Cuda/Kokkos_Cuda_View.hpp b/lib/kokkos/core/src/Cuda/Kokkos_Cuda_View.hpp
index 61563a0100..dec6ef15e1 100644
--- a/lib/kokkos/core/src/Cuda/Kokkos_Cuda_View.hpp
+++ b/lib/kokkos/core/src/Cuda/Kokkos_Cuda_View.hpp
@@ -139,7 +139,7 @@ struct CudaLDGFetch {
 
   template <typename iType>
   KOKKOS_INLINE_FUNCTION ValueType operator[](const iType& i) const {
-#if defined(__CUDA_ARCH__) && (350 <= _CUDA_ARCH__)
+#if defined(__CUDA_ARCH__) && (350 <= __CUDA_ARCH__)
     AliasType v = __ldg(reinterpret_cast<const AliasType*>(&m_ptr[i]));
     return *(reinterpret_cast<ValueType*>(&v));
 #else
diff --git a/lib/kokkos/core/src/HIP/Kokkos_HIP_Instance.cpp b/lib/kokkos/core/src/HIP/Kokkos_HIP_Instance.cpp
index 4a6a3ba99e..a8a0496afe 100644
--- a/lib/kokkos/core/src/HIP/Kokkos_HIP_Instance.cpp
+++ b/lib/kokkos/core/src/HIP/Kokkos_HIP_Instance.cpp
@@ -132,7 +132,8 @@ void HIPInternal::print_configuration(std::ostream &s) const {
   s << "macro  KOKKOS_ENABLE_HIP : defined" << '\n';
 #if defined(HIP_VERSION)
   s << "macro  HIP_VERSION = " << HIP_VERSION << " = version "
-    << HIP_VERSION / 100 << "." << HIP_VERSION % 100 << '\n';
+    << HIP_VERSION_MAJOR << '.' << HIP_VERSION_MINOR << '.' << HIP_VERSION_PATCH
+    << '\n';
 #endif
 
   for (int i = 0; i < dev_info.m_hipDevCount; ++i) {
@@ -467,7 +468,6 @@ void HIPInternal::finalize() {
 }
 
 char *HIPInternal::get_next_driver(size_t driverTypeSize) const {
-  std::lock_guard<std::mutex> const lock(m_mutexWorkArray);
   if (d_driverWorkArray == nullptr) {
     KOKKOS_IMPL_HIP_SAFE_CALL(
         hipHostMalloc(&d_driverWorkArray,
diff --git a/lib/kokkos/core/src/HIP/Kokkos_HIP_KernelLaunch.hpp b/lib/kokkos/core/src/HIP/Kokkos_HIP_KernelLaunch.hpp
index 384b7ffd67..70b979e00a 100644
--- a/lib/kokkos/core/src/HIP/Kokkos_HIP_KernelLaunch.hpp
+++ b/lib/kokkos/core/src/HIP/Kokkos_HIP_KernelLaunch.hpp
@@ -490,6 +490,8 @@ struct HIPParallelLaunch<
 
       KOKKOS_ENSURE_HIP_LOCK_ARRAYS_ON_DEVICE();
 
+      std::lock_guard<std::mutex> const lock(hip_instance->m_mutexWorkArray);
+
       // Invoke the driver function on the device
       DriverType *d_driver = reinterpret_cast<DriverType *>(
           hip_instance->get_next_driver(sizeof(DriverType)));
diff --git a/lib/kokkos/core/src/HIP/Kokkos_HIP_Locks.cpp b/lib/kokkos/core/src/HIP/Kokkos_HIP_Locks.cpp
index f334d93412..e9cfbf99f7 100644
--- a/lib/kokkos/core/src/HIP/Kokkos_HIP_Locks.cpp
+++ b/lib/kokkos/core/src/HIP/Kokkos_HIP_Locks.cpp
@@ -56,8 +56,7 @@ namespace Kokkos {
 
 #ifdef KOKKOS_ENABLE_HIP_RELOCATABLE_DEVICE_CODE
 namespace Impl {
-__device__ __constant__ HIPLockArrays g_device_hip_lock_arrays = {nullptr,
-                                                                  nullptr, 0};
+__device__ __constant__ HIPLockArrays g_device_hip_lock_arrays = {nullptr, 0};
 }
 #endif
 
diff --git a/lib/kokkos/core/src/HIP/Kokkos_HIP_Space.cpp b/lib/kokkos/core/src/HIP/Kokkos_HIP_Space.cpp
index 6ade677fa8..776b7c6abe 100644
--- a/lib/kokkos/core/src/HIP/Kokkos_HIP_Space.cpp
+++ b/lib/kokkos/core/src/HIP/Kokkos_HIP_Space.cpp
@@ -464,6 +464,15 @@ void HIPSpaceInitializer::print_configuration(std::ostream& msg,
 }
 
 }  // namespace Impl
+
+#ifdef KOKKOS_ENABLE_CXX14
+namespace Tools {
+namespace Experimental {
+constexpr DeviceType DeviceTypeTraits<Kokkos::Experimental::HIP>::id;
+}
+}  // namespace Tools
+#endif
+
 }  // namespace Kokkos
 
 //==============================================================================
diff --git a/lib/kokkos/core/src/HPX/Kokkos_HPX.cpp b/lib/kokkos/core/src/HPX/Kokkos_HPX.cpp
index acf2224f02..623c7da025 100644
--- a/lib/kokkos/core/src/HPX/Kokkos_HPX.cpp
+++ b/lib/kokkos/core/src/HPX/Kokkos_HPX.cpp
@@ -199,6 +199,15 @@ void HPXSpaceInitializer::print_configuration(std::ostream &msg,
 }
 
 }  // namespace Impl
+
+#ifdef KOKKOS_ENABLE_CXX14
+namespace Tools {
+namespace Experimental {
+constexpr DeviceType DeviceTypeTraits<Kokkos::Experimental::HPX>::id;
+}
+}  // namespace Tools
+#endif
+
 }  // namespace Kokkos
 
 #else
diff --git a/lib/kokkos/core/src/Kokkos_Cuda.hpp b/lib/kokkos/core/src/Kokkos_Cuda.hpp
index 6305a1fa5d..0063b1cd1e 100644
--- a/lib/kokkos/core/src/Kokkos_Cuda.hpp
+++ b/lib/kokkos/core/src/Kokkos_Cuda.hpp
@@ -260,6 +260,7 @@ template <>
 struct DeviceTypeTraits<Cuda> {
   /// \brief An ID to differentiate (for example) Serial from OpenMP in Tooling
   static constexpr DeviceType id = DeviceType::Cuda;
+  static int device_id(const Cuda& exec) { return exec.cuda_device(); }
 };
 }  // namespace Experimental
 }  // namespace Tools
diff --git a/lib/kokkos/core/src/Kokkos_HIP_Space.hpp b/lib/kokkos/core/src/Kokkos_HIP_Space.hpp
index 1371d21d38..68869a6074 100644
--- a/lib/kokkos/core/src/Kokkos_HIP_Space.hpp
+++ b/lib/kokkos/core/src/Kokkos_HIP_Space.hpp
@@ -571,6 +571,9 @@ namespace Experimental {
 template <>
 struct DeviceTypeTraits<Kokkos::Experimental::HIP> {
   static constexpr DeviceType id = DeviceType::HIP;
+  static int device_id(const Kokkos::Experimental::HIP& exec) {
+    return exec.hip_device();
+  }
 };
 }  // namespace Experimental
 }  // namespace Tools
diff --git a/lib/kokkos/core/src/Kokkos_HPX.hpp b/lib/kokkos/core/src/Kokkos_HPX.hpp
index d2ae9c0ec2..9238ca30a7 100644
--- a/lib/kokkos/core/src/Kokkos_HPX.hpp
+++ b/lib/kokkos/core/src/Kokkos_HPX.hpp
@@ -500,6 +500,7 @@ namespace Experimental {
 template <>
 struct DeviceTypeTraits<Kokkos::Experimental::HPX> {
   static constexpr DeviceType id = DeviceType::HPX;
+  static int device_id(const Kokkos::Experimental::HPX &) { return 0; }
 };
 }  // namespace Experimental
 }  // namespace Tools
diff --git a/lib/kokkos/core/src/Kokkos_OpenMP.hpp b/lib/kokkos/core/src/Kokkos_OpenMP.hpp
index 5d76e689f2..767e5b9324 100644
--- a/lib/kokkos/core/src/Kokkos_OpenMP.hpp
+++ b/lib/kokkos/core/src/Kokkos_OpenMP.hpp
@@ -179,6 +179,7 @@ namespace Experimental {
 template <>
 struct DeviceTypeTraits<OpenMP> {
   static constexpr DeviceType id = DeviceType::OpenMP;
+  static int device_id(const OpenMP&) { return 0; }
 };
 }  // namespace Experimental
 }  // namespace Tools
diff --git a/lib/kokkos/core/src/Kokkos_OpenMPTarget.hpp b/lib/kokkos/core/src/Kokkos_OpenMPTarget.hpp
index f394f32408..373dc3d9c7 100644
--- a/lib/kokkos/core/src/Kokkos_OpenMPTarget.hpp
+++ b/lib/kokkos/core/src/Kokkos_OpenMPTarget.hpp
@@ -130,6 +130,9 @@ template <>
 struct DeviceTypeTraits<::Kokkos::Experimental::OpenMPTarget> {
   static constexpr DeviceType id =
       ::Kokkos::Profiling::Experimental::DeviceType::OpenMPTarget;
+  static int device_id(const Kokkos::Experimental::OpenMPTarget&) {
+    return omp_get_default_device();
+  }
 };
 }  // namespace Experimental
 }  // namespace Tools
diff --git a/lib/kokkos/core/src/Kokkos_SYCL.hpp b/lib/kokkos/core/src/Kokkos_SYCL.hpp
index 02095ff7b3..e29093db32 100644
--- a/lib/kokkos/core/src/Kokkos_SYCL.hpp
+++ b/lib/kokkos/core/src/Kokkos_SYCL.hpp
@@ -182,6 +182,9 @@ template <>
 struct DeviceTypeTraits<Kokkos::Experimental::SYCL> {
   /// \brief An ID to differentiate (for example) Serial from OpenMP in Tooling
   static constexpr DeviceType id = DeviceType::SYCL;
+  static int device_id(const Kokkos::Experimental::SYCL& exec) {
+    return exec.sycl_device();
+  }
 };
 }  // namespace Experimental
 }  // namespace Tools
diff --git a/lib/kokkos/core/src/Kokkos_Serial.hpp b/lib/kokkos/core/src/Kokkos_Serial.hpp
index 9aada48bf6..b2e524c374 100644
--- a/lib/kokkos/core/src/Kokkos_Serial.hpp
+++ b/lib/kokkos/core/src/Kokkos_Serial.hpp
@@ -226,6 +226,7 @@ namespace Experimental {
 template <>
 struct DeviceTypeTraits<Serial> {
   static constexpr DeviceType id = DeviceType::Serial;
+  static int device_id(const Serial&) { return 0; }
 };
 }  // namespace Experimental
 }  // namespace Tools
diff --git a/lib/kokkos/core/src/Kokkos_Threads.hpp b/lib/kokkos/core/src/Kokkos_Threads.hpp
index 45a2d0e326..5879209f12 100644
--- a/lib/kokkos/core/src/Kokkos_Threads.hpp
+++ b/lib/kokkos/core/src/Kokkos_Threads.hpp
@@ -175,6 +175,7 @@ namespace Experimental {
 template <>
 struct DeviceTypeTraits<Threads> {
   static constexpr DeviceType id = DeviceType::Threads;
+  static int device_id(const Threads&) { return 0; }
 };
 }  // namespace Experimental
 }  // namespace Tools
diff --git a/lib/kokkos/core/src/OpenMP/Kokkos_OpenMP_Exec.cpp b/lib/kokkos/core/src/OpenMP/Kokkos_OpenMP_Exec.cpp
index d2283d456f..66dbbacce9 100644
--- a/lib/kokkos/core/src/OpenMP/Kokkos_OpenMP_Exec.cpp
+++ b/lib/kokkos/core/src/OpenMP/Kokkos_OpenMP_Exec.cpp
@@ -67,8 +67,9 @@ __thread int t_openmp_hardware_id            = 0;
 __thread Impl::OpenMPExec *t_openmp_instance = nullptr;
 
 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE_3
-void OpenMPExec::validate_partition(const int nthreads, int &num_partitions,
-                                    int &partition_size) {
+void OpenMPExec::validate_partition_impl(const int nthreads,
+                                         int &num_partitions,
+                                         int &partition_size) {
   if (nthreads == 1) {
     num_partitions = 1;
     partition_size = 1;
@@ -506,6 +507,15 @@ void OpenMPSpaceInitializer::print_configuration(std::ostream &msg,
 }
 
 }  // namespace Impl
+
+#ifdef KOKKOS_ENABLE_CXX14
+namespace Tools {
+namespace Experimental {
+constexpr DeviceType DeviceTypeTraits<OpenMP>::id;
+}
+}  // namespace Tools
+#endif
+
 }  // namespace Kokkos
 
 #else
diff --git a/lib/kokkos/core/src/OpenMP/Kokkos_OpenMP_Exec.hpp b/lib/kokkos/core/src/OpenMP/Kokkos_OpenMP_Exec.hpp
index 2f647af77e..ede24d1094 100644
--- a/lib/kokkos/core/src/OpenMP/Kokkos_OpenMP_Exec.hpp
+++ b/lib/kokkos/core/src/OpenMP/Kokkos_OpenMP_Exec.hpp
@@ -93,7 +93,11 @@ class OpenMPExec {
 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE_3
   KOKKOS_DEPRECATED static void validate_partition(const int nthreads,
                                                    int& num_partitions,
-                                                   int& partition_size);
+                                                   int& partition_size) {
+    validate_partition_impl(nthreads, num_partitions, partition_size);
+  }
+  static void validate_partition_impl(const int nthreads, int& num_partitions,
+                                      int& partition_size);
 #endif
 
  private:
@@ -179,8 +183,8 @@ KOKKOS_DEPRECATED void OpenMP::partition_master(F const& f, int num_partitions,
 
     Exec* prev_instance = Impl::t_openmp_instance;
 
-    Exec::validate_partition(prev_instance->m_pool_size, num_partitions,
-                             partition_size);
+    Exec::validate_partition_impl(prev_instance->m_pool_size, num_partitions,
+                                  partition_size);
 
     OpenMP::memory_space space;
 
diff --git a/lib/kokkos/core/src/SYCL/Kokkos_SYCL_Instance.hpp b/lib/kokkos/core/src/SYCL/Kokkos_SYCL_Instance.hpp
index 907e4e9efe..45aacd7258 100644
--- a/lib/kokkos/core/src/SYCL/Kokkos_SYCL_Instance.hpp
+++ b/lib/kokkos/core/src/SYCL/Kokkos_SYCL_Instance.hpp
@@ -72,7 +72,7 @@ class SYCLInternal {
                                   bool force_shrink = false);
 
   uint32_t impl_get_instance_id() const;
-  int m_syclDev = -1;
+  int m_syclDev = 0;
 
   size_t m_maxWorkgroupSize   = 0;
   uint32_t m_maxConcurrency   = 0;
diff --git a/lib/kokkos/core/src/Threads/Kokkos_ThreadsExec.cpp b/lib/kokkos/core/src/Threads/Kokkos_ThreadsExec.cpp
index 8a7c49871b..9682564ee0 100644
--- a/lib/kokkos/core/src/Threads/Kokkos_ThreadsExec.cpp
+++ b/lib/kokkos/core/src/Threads/Kokkos_ThreadsExec.cpp
@@ -399,27 +399,68 @@ bool ThreadsExec::wake() {
 
 //----------------------------------------------------------------------------
 
+void ThreadsExec::execute_resize_scratch_in_serial() {
+  const unsigned begin = s_threads_process.m_pool_base ? 1 : 0;
+
+  auto deallocate_scratch_memory = [](ThreadsExec &exec) {
+    if (exec.m_scratch) {
+      using Record =
+          Kokkos::Impl::SharedAllocationRecord<Kokkos::HostSpace, void>;
+      Record *const r = Record::get_record(exec.m_scratch);
+      exec.m_scratch  = nullptr;
+      Record::decrement(r);
+    }
+  };
+  if (s_threads_process.m_pool_base) {
+    for (unsigned i = s_thread_pool_size[0]; begin < i;) {
+      deallocate_scratch_memory(*s_threads_exec[--i]);
+    }
+  }
+
+  s_current_function     = &first_touch_allocate_thread_private_scratch;
+  s_current_function_arg = &s_threads_process;
+
+  // Make sure function and arguments are written before activating threads.
+  memory_fence();
+
+  for (unsigned i = s_thread_pool_size[0]; begin < i;) {
+    ThreadsExec &th = *s_threads_exec[--i];
+
+    th.m_pool_state = ThreadsExec::Active;
+
+    wait_yield(th.m_pool_state, ThreadsExec::Active);
+  }
+
+  if (s_threads_process.m_pool_base) {
+    deallocate_scratch_memory(s_threads_process);
+    s_threads_process.m_pool_state = ThreadsExec::Active;
+    first_touch_allocate_thread_private_scratch(s_threads_process, nullptr);
+    s_threads_process.m_pool_state = ThreadsExec::Inactive;
+  }
+
+  s_current_function_arg = nullptr;
+  s_current_function     = nullptr;
+
+  // Make sure function and arguments are cleared before proceeding.
+  memory_fence();
+}
+
+//----------------------------------------------------------------------------
+
 void *ThreadsExec::root_reduce_scratch() {
   return s_threads_process.reduce_memory();
 }
 
-void ThreadsExec::execute_resize_scratch(ThreadsExec &exec, const void *) {
-  using Record = Kokkos::Impl::SharedAllocationRecord<Kokkos::HostSpace, void>;
-
-  if (exec.m_scratch) {
-    Record *const r = Record::get_record(exec.m_scratch);
-
-    exec.m_scratch = nullptr;
-
-    Record::decrement(r);
-  }
-
+void ThreadsExec::first_touch_allocate_thread_private_scratch(ThreadsExec &exec,
+                                                              const void *) {
   exec.m_scratch_reduce_end = s_threads_process.m_scratch_reduce_end;
   exec.m_scratch_thread_end = s_threads_process.m_scratch_thread_end;
 
   if (s_threads_process.m_scratch_thread_end) {
     // Allocate tracked memory:
     {
+      using Record =
+          Kokkos::Impl::SharedAllocationRecord<Kokkos::HostSpace, void>;
       Record *const r =
           Record::allocate(Kokkos::HostSpace(), "Kokkos::thread_scratch",
                            s_threads_process.m_scratch_thread_end);
@@ -461,7 +502,7 @@ void *ThreadsExec::resize_scratch(size_t reduce_size, size_t thread_size) {
     s_threads_process.m_scratch_reduce_end = reduce_size;
     s_threads_process.m_scratch_thread_end = reduce_size + thread_size;
 
-    execute_resize_scratch(s_threads_process, nullptr);
+    execute_resize_scratch_in_serial();
 
     s_threads_process.m_scratch = s_threads_exec[0]->m_scratch;
   }
@@ -845,6 +886,15 @@ void ThreadsSpaceInitializer::print_configuration(std::ostream &msg,
 }
 
 }  // namespace Impl
+
+#ifdef KOKKOS_ENABLE_CXX14
+namespace Tools {
+namespace Experimental {
+constexpr DeviceType DeviceTypeTraits<Threads>::id;
+}
+}  // namespace Tools
+#endif
+
 } /* namespace Kokkos */
 //----------------------------------------------------------------------------
 //----------------------------------------------------------------------------
diff --git a/lib/kokkos/core/src/Threads/Kokkos_ThreadsExec.hpp b/lib/kokkos/core/src/Threads/Kokkos_ThreadsExec.hpp
index 561b1ce292..d17f417bbc 100644
--- a/lib/kokkos/core/src/Threads/Kokkos_ThreadsExec.hpp
+++ b/lib/kokkos/core/src/Threads/Kokkos_ThreadsExec.hpp
@@ -123,12 +123,15 @@ class ThreadsExec {
   static void global_unlock();
   static void spawn();
 
-  static void execute_resize_scratch(ThreadsExec &, const void *);
+  static void first_touch_allocate_thread_private_scratch(ThreadsExec &,
+                                                          const void *);
   static void execute_sleep(ThreadsExec &, const void *);
 
   ThreadsExec(const ThreadsExec &);
   ThreadsExec &operator=(const ThreadsExec &);
 
+  static void execute_resize_scratch_in_serial();
+
  public:
   KOKKOS_INLINE_FUNCTION int pool_size() const { return m_pool_size; }
   KOKKOS_INLINE_FUNCTION int pool_rank() const { return m_pool_rank; }
diff --git a/lib/kokkos/core/src/impl/Kokkos_Profiling_Interface.hpp b/lib/kokkos/core/src/impl/Kokkos_Profiling_Interface.hpp
index 4e0e81405f..d526682056 100644
--- a/lib/kokkos/core/src/impl/Kokkos_Profiling_Interface.hpp
+++ b/lib/kokkos/core/src/impl/Kokkos_Profiling_Interface.hpp
@@ -118,11 +118,14 @@ template <typename ExecutionSpace>
 constexpr uint32_t device_id_root() {
   constexpr auto device_id =
       static_cast<uint32_t>(DeviceTypeTraits<ExecutionSpace>::id);
-  return (device_id << num_instance_bits);
+  return (device_id << (num_instance_bits + num_device_bits));
 }
 template <typename ExecutionSpace>
 inline uint32_t device_id(ExecutionSpace const& space) noexcept {
-  return device_id_root<ExecutionSpace>() + space.impl_instance_id();
+  return device_id_root<ExecutionSpace>() +
+         (DeviceTypeTraits<ExecutionSpace>::device_id(space)
+          << num_instance_bits) +
+         space.impl_instance_id();
 }
 }  // namespace Experimental
 }  // namespace Tools
diff --git a/lib/kokkos/core/src/impl/Kokkos_Serial.cpp b/lib/kokkos/core/src/impl/Kokkos_Serial.cpp
index c49e838d8f..e5917eb59d 100644
--- a/lib/kokkos/core/src/impl/Kokkos_Serial.cpp
+++ b/lib/kokkos/core/src/impl/Kokkos_Serial.cpp
@@ -233,6 +233,15 @@ void SerialSpaceInitializer::print_configuration(std::ostream& msg,
 }
 
 }  // namespace Impl
+
+#ifdef KOKKOS_ENABLE_CXX14
+namespace Tools {
+namespace Experimental {
+constexpr DeviceType DeviceTypeTraits<Serial>::id;
+}
+}  // namespace Tools
+#endif
+
 }  // namespace Kokkos
 
 #else
diff --git a/lib/kokkos/core/src/impl/Kokkos_ViewMapping.hpp b/lib/kokkos/core/src/impl/Kokkos_ViewMapping.hpp
index 09f7af0918..f606a39839 100644
--- a/lib/kokkos/core/src/impl/Kokkos_ViewMapping.hpp
+++ b/lib/kokkos/core/src/impl/Kokkos_ViewMapping.hpp
@@ -1005,15 +1005,15 @@ struct ViewOffset<
   /* Cardinality of the domain index space */
   KOKKOS_INLINE_FUNCTION
   constexpr size_type size() const {
-    return m_dim.N0 * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4 * m_dim.N5 *
-           m_dim.N6 * m_dim.N7;
+    return size_type(m_dim.N0) * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4 *
+           m_dim.N5 * m_dim.N6 * m_dim.N7;
   }
 
   /* Span of the range space */
   KOKKOS_INLINE_FUNCTION
   constexpr size_type span() const {
-    return m_dim.N0 * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4 * m_dim.N5 *
-           m_dim.N6 * m_dim.N7;
+    return size_type(m_dim.N0) * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4 *
+           m_dim.N5 * m_dim.N6 * m_dim.N7;
   }
 
   KOKKOS_INLINE_FUNCTION constexpr bool span_is_contiguous() const {
@@ -1026,23 +1026,24 @@ struct ViewOffset<
     return m_dim.N0;
   }
   KOKKOS_INLINE_FUNCTION constexpr size_type stride_2() const {
-    return m_dim.N0 * m_dim.N1;
+    return size_type(m_dim.N0) * m_dim.N1;
   }
   KOKKOS_INLINE_FUNCTION constexpr size_type stride_3() const {
-    return m_dim.N0 * m_dim.N1 * m_dim.N2;
+    return size_type(m_dim.N0) * m_dim.N1 * m_dim.N2;
   }
   KOKKOS_INLINE_FUNCTION constexpr size_type stride_4() const {
-    return m_dim.N0 * m_dim.N1 * m_dim.N2 * m_dim.N3;
+    return size_type(m_dim.N0) * m_dim.N1 * m_dim.N2 * m_dim.N3;
   }
   KOKKOS_INLINE_FUNCTION constexpr size_type stride_5() const {
-    return m_dim.N0 * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4;
+    return size_type(m_dim.N0) * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4;
   }
   KOKKOS_INLINE_FUNCTION constexpr size_type stride_6() const {
-    return m_dim.N0 * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4 * m_dim.N5;
+    return size_type(m_dim.N0) * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4 *
+           m_dim.N5;
   }
   KOKKOS_INLINE_FUNCTION constexpr size_type stride_7() const {
-    return m_dim.N0 * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4 * m_dim.N5 *
-           m_dim.N6;
+    return size_type(m_dim.N0) * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4 *
+           m_dim.N5 * m_dim.N6;
   }
 
   // Stride with [ rank ] value is the total length
@@ -1288,8 +1289,8 @@ struct ViewOffset<
   /* Cardinality of the domain index space */
   KOKKOS_INLINE_FUNCTION
   constexpr size_type size() const {
-    return m_dim.N0 * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4 * m_dim.N5 *
-           m_dim.N6 * m_dim.N7;
+    return size_type(m_dim.N0) * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4 *
+           m_dim.N5 * m_dim.N6 * m_dim.N7;
   }
 
   /* Span of the range space */
@@ -1633,15 +1634,15 @@ struct ViewOffset<
   /* Cardinality of the domain index space */
   KOKKOS_INLINE_FUNCTION
   constexpr size_type size() const {
-    return m_dim.N0 * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4 * m_dim.N5 *
-           m_dim.N6 * m_dim.N7;
+    return size_type(m_dim.N0) * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4 *
+           m_dim.N5 * m_dim.N6 * m_dim.N7;
   }
 
   /* Span of the range space */
   KOKKOS_INLINE_FUNCTION
   constexpr size_type span() const {
-    return m_dim.N0 * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4 * m_dim.N5 *
-           m_dim.N6 * m_dim.N7;
+    return size_type(m_dim.N0) * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4 *
+           m_dim.N5 * m_dim.N6 * m_dim.N7;
   }
 
   KOKKOS_INLINE_FUNCTION constexpr bool span_is_contiguous() const {
@@ -1916,14 +1917,14 @@ struct ViewOffset<
   /* Cardinality of the domain index space */
   KOKKOS_INLINE_FUNCTION
   constexpr size_type size() const {
-    return m_dim.N0 * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4 * m_dim.N5 *
-           m_dim.N6 * m_dim.N7;
+    return size_type(m_dim.N0) * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4 *
+           m_dim.N5 * m_dim.N6 * m_dim.N7;
   }
 
   /* Span of the range space */
   KOKKOS_INLINE_FUNCTION
   constexpr size_type span() const {
-    return size() > 0 ? m_dim.N0 * m_stride : 0;
+    return size() > 0 ? size_type(m_dim.N0) * m_stride : 0;
   }
 
   KOKKOS_INLINE_FUNCTION constexpr bool span_is_contiguous() const {
@@ -2066,27 +2067,29 @@ struct ViewOffset<
                 stride(/* 2 <= rank */
                        m_dim.N1 *
                        (dimension_type::rank == 2
-                            ? 1
+                            ? size_t(1)
                             : m_dim.N2 *
                                   (dimension_type::rank == 3
-                                       ? 1
+                                       ? size_t(1)
                                        : m_dim.N3 *
                                              (dimension_type::rank == 4
-                                                  ? 1
+                                                  ? size_t(1)
                                                   : m_dim.N4 *
                                                         (dimension_type::rank ==
                                                                  5
-                                                             ? 1
+                                                             ? size_t(1)
                                                              : m_dim.N5 *
                                                                    (dimension_type::
                                                                                 rank ==
                                                                             6
-                                                                        ? 1
+                                                                        ? size_t(
+                                                                              1)
                                                                         : m_dim.N6 *
                                                                               (dimension_type::
                                                                                            rank ==
                                                                                        7
-                                                                                   ? 1
+                                                                                   ? size_t(
+                                                                                         1)
                                                                                    : m_dim
                                                                                          .N7)))))))) {
   }
@@ -2447,8 +2450,8 @@ struct ViewOffset<Dimension, Kokkos::LayoutStride, void> {
   constexpr size_type size() const {
     return dimension_type::rank == 0
                ? 1
-               : m_dim.N0 * m_dim.N1 * m_dim.N2 * m_dim.N3 * m_dim.N4 *
-                     m_dim.N5 * m_dim.N6 * m_dim.N7;
+               : size_type(m_dim.N0) * m_dim.N1 * m_dim.N2 * m_dim.N3 *
+                     m_dim.N4 * m_dim.N5 * m_dim.N6 * m_dim.N7;
   }
 
  private:
diff --git a/lib/kokkos/core/src/impl/Kokkos_ViewTracker.hpp b/lib/kokkos/core/src/impl/Kokkos_ViewTracker.hpp
index fe3651886b..972b1b6d9a 100644
--- a/lib/kokkos/core/src/impl/Kokkos_ViewTracker.hpp
+++ b/lib/kokkos/core/src/impl/Kokkos_ViewTracker.hpp
@@ -91,6 +91,7 @@ struct ViewTracker {
 
   template <class RT, class... RP>
   KOKKOS_INLINE_FUNCTION void assign(const View<RT, RP...>& vt) noexcept {
+    if (this == reinterpret_cast<const ViewTracker*>(&vt.m_track)) return;
     KOKKOS_IF_ON_HOST((
         if (view_traits::is_managed && Kokkos::Impl::SharedAllocationRecord<
                                            void, void>::tracking_enabled()) {
@@ -102,6 +103,7 @@ struct ViewTracker {
 
   KOKKOS_INLINE_FUNCTION ViewTracker& operator=(
       const ViewTracker& rhs) noexcept {
+    if (this == &rhs) return *this;
     KOKKOS_IF_ON_HOST((
         if (view_traits::is_managed && Kokkos::Impl::SharedAllocationRecord<
                                            void, void>::tracking_enabled()) {
diff --git a/lib/kokkos/core/unit_test/TestViewAPI.hpp b/lib/kokkos/core/unit_test/TestViewAPI.hpp
index 21602be086..83efae6170 100644
--- a/lib/kokkos/core/unit_test/TestViewAPI.hpp
+++ b/lib/kokkos/core/unit_test/TestViewAPI.hpp
@@ -1087,6 +1087,20 @@ class TestViewAPI {
     dView4_unmanaged unmanaged_dx = dx;
     ASSERT_EQ(dx.use_count(), 1);
 
+    // Test self assignment
+#if defined(__clang__)
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wself-assign-overloaded"
+#endif
+    dx = dx;  // copy-assignment operator
+#if defined(__clang__)
+#pragma GCC diagnostic pop
+#endif
+    ASSERT_EQ(dx.use_count(), 1);
+    dx = reinterpret_cast<typename dView4::uniform_type &>(
+        dx);  // conversion assignment operator
+    ASSERT_EQ(dx.use_count(), 1);
+
     dView4_unmanaged unmanaged_from_ptr_dx = dView4_unmanaged(
         dx.data(), dx.extent(0), dx.extent(1), dx.extent(2), dx.extent(3));
 
diff --git a/lib/kokkos/core/unit_test/TestViewAPI_e.hpp b/lib/kokkos/core/unit_test/TestViewAPI_e.hpp
index d4f484a530..d1d38022a7 100644
--- a/lib/kokkos/core/unit_test/TestViewAPI_e.hpp
+++ b/lib/kokkos/core/unit_test/TestViewAPI_e.hpp
@@ -240,6 +240,35 @@ struct TestViewOverloadResolution {
 TEST(TEST_CATEGORY, view_overload_resolution) {
   TestViewOverloadResolution<TEST_EXECSPACE>::test_function_overload();
 }
+
+template <typename MemorySpace>
+struct TestViewAllocationLargeRank {
+  using ViewType = Kokkos::View<char********, MemorySpace>;
+
+  KOKKOS_FUNCTION void operator()(int) const {
+    size_t idx = v.extent(0) - 1;
+    auto& lhs  = v(idx, idx, idx, idx, idx, idx, idx, idx);
+    lhs        = 42;  // This is where it segfaulted
+  }
+
+  ViewType v;
+};
+
+TEST(TEST_CATEGORY, view_allocation_large_rank) {
+  using ExecutionSpace = typename TEST_EXECSPACE::execution_space;
+  using MemorySpace    = typename TEST_EXECSPACE::memory_space;
+  constexpr int dim    = 16;
+  using FunctorType    = TestViewAllocationLargeRank<MemorySpace>;
+  typename FunctorType::ViewType v("v", dim, dim, dim, dim, dim, dim, dim, dim);
+
+  Kokkos::parallel_for(Kokkos::RangePolicy<ExecutionSpace>(0, 1),
+                       FunctorType{v});
+  typename FunctorType::ViewType v_single(v.data() + v.size() - 1, 1, 1, 1, 1,
+                                          1, 1, 1, 1);
+  auto result =
+      Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, v_single);
+  ASSERT_EQ(result(0, 0, 0, 0, 0, 0, 0, 0), 42);
+}
 }  // namespace Test
 
 #include <TestViewIsAssignable.hpp>
diff --git a/lib/kokkos/core/unit_test/tools/TestEventCorrectness.hpp b/lib/kokkos/core/unit_test/tools/TestEventCorrectness.hpp
index 08863232ed..bb1d3156f5 100644
--- a/lib/kokkos/core/unit_test/tools/TestEventCorrectness.hpp
+++ b/lib/kokkos/core/unit_test/tools/TestEventCorrectness.hpp
@@ -238,13 +238,10 @@ TEST(kokkosp, test_id_gen) {
   using Kokkos::Tools::Experimental::DeviceTypeTraits;
   test_wrapper([&]() {
     Kokkos::DefaultExecutionSpace ex;
-    auto id      = device_id(ex);
-    auto id_ref  = identifier_from_devid(id);
-    auto success = (id_ref.instance_id == ex.impl_instance_id()) &&
-                   (id_ref.device_id ==
-                    static_cast<uint32_t>(
-                        DeviceTypeTraits<Kokkos::DefaultExecutionSpace>::id));
-    ASSERT_TRUE(success);
+    auto id     = device_id(ex);
+    auto id_ref = identifier_from_devid(id);
+    ASSERT_EQ(DeviceTypeTraits<decltype(ex)>::id, id_ref.type);
+    ASSERT_EQ(id_ref.instance_id, ex.impl_instance_id());
   });
 }
 
@@ -253,6 +250,7 @@ TEST(kokkosp, test_id_gen) {
  */
 TEST(kokkosp, test_kernel_sequence) {
   test_wrapper([&]() {
+    Kokkos::DefaultExecutionSpace ex;
     auto root = Kokkos::Tools::Experimental::device_id_root<
         Kokkos::DefaultExecutionSpace>();
     std::vector<FencePayload> expected{
@@ -260,11 +258,10 @@ TEST(kokkosp, test_kernel_sequence) {
         {"named_instance", FencePayload::distinguishable_devices::no,
          root + num_instances},
         {"test_kernel", FencePayload::distinguishable_devices::no,
-         root + num_instances}
+         Kokkos::Tools::Experimental::device_id(ex)}
 
     };
     expect_fence_events(expected, [=]() {
-      Kokkos::DefaultExecutionSpace ex;
       TestFunctor tf;
       ex.fence("named_instance");
       Kokkos::parallel_for(
diff --git a/lib/kokkos/master_history.txt b/lib/kokkos/master_history.txt
index e174b47f67..41c755a8a8 100644
--- a/lib/kokkos/master_history.txt
+++ b/lib/kokkos/master_history.txt
@@ -27,3 +27,4 @@ tag:  3.4.00     date: 04:26:2021    master: 1fb0c284    release: 5d7738d6
 tag:  3.4.01     date: 05:20:2021    master: 4b97a22f    release: 410b15c8
 tag:  3.5.00     date: 11:19:2021    master: c28a8b03    release: 21b879e4
 tag:  3.6.00     date: 04:14:2022    master: 2834f94a    release: 6ea708ff
+tag:  3.6.01     date: 06:16:2022    master: b52f8c83    release: afe9b404
diff --git a/src/EXTRA-PAIR/pair_coul_slater_long.cpp b/src/EXTRA-PAIR/pair_coul_slater_long.cpp
index 65619f0f65..6662ae3721 100644
--- a/src/EXTRA-PAIR/pair_coul_slater_long.cpp
+++ b/src/EXTRA-PAIR/pair_coul_slater_long.cpp
@@ -116,15 +116,15 @@ void PairCoulSlaterLong::compute(int eflag, int vflag)
 
       if (rsq < cut_coulsq) {
         r2inv = 1.0/rsq;
-          r = sqrt(rsq);
-          grij = g_ewald * r;
-          expm2 = exp(-grij*grij);
-          t = 1.0 / (1.0 + EWALD_P*grij);
-          erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
-          slater_term = exp(-2*r/lamda)*(1 + (2*r/lamda*(1+r/lamda)));
-          prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r;
-          forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - slater_term);
-          if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
+        r = sqrt(rsq);
+        grij = g_ewald * r;
+        expm2 = exp(-grij*grij);
+        t = 1.0 / (1.0 + EWALD_P*grij);
+        erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+        slater_term = exp(-2*r/lamda)*(1 + (2*r/lamda*(1+r/lamda)));
+        prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r;
+        forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - slater_term);
+        if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor*(1-slater_term);
 
         fpair = forcecoul * r2inv;
 
@@ -138,8 +138,8 @@ void PairCoulSlaterLong::compute(int eflag, int vflag)
         }
 
         if (eflag) {
-            ecoul = prefactor*(erfc - (1 + r/lamda)*exp(-2*r/lamda));
-          if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
+          ecoul = prefactor*(erfc - (1 + r/lamda)*exp(-2*r/lamda));
+          if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor*(1.0-(1 + r/lamda)*exp(-2*r/lamda));
         }
 
         if (evflag) ev_tally(i,j,nlocal,newton_pair,
diff --git a/src/INTEL/pair_sw_intel.cpp b/src/INTEL/pair_sw_intel.cpp
index a882f8bfba..37fe19260a 100644
--- a/src/INTEL/pair_sw_intel.cpp
+++ b/src/INTEL/pair_sw_intel.cpp
@@ -1101,7 +1101,11 @@ void PairSWIntel::allocate()
 
 void PairSWIntel::init_style()
 {
+  // there is no support for skipping threebody loops (yet)
+  bool tmp_threebody = skip_threebody_flag;
+  skip_threebody_flag = false;
   PairSW::init_style();
+  skip_threebody_flag = tmp_threebody;
 
   map[0] = map[1];
 
diff --git a/src/KOKKOS/pair_sw_kokkos.cpp b/src/KOKKOS/pair_sw_kokkos.cpp
index cae0aea0e8..a3d8d11380 100644
--- a/src/KOKKOS/pair_sw_kokkos.cpp
+++ b/src/KOKKOS/pair_sw_kokkos.cpp
@@ -392,7 +392,11 @@ void PairSWKokkos<DeviceType>::coeff(int narg, char **arg)
 template<class DeviceType>
 void PairSWKokkos<DeviceType>::init_style()
 {
+  // there is no support for skipping threebody loops (yet)
+  bool tmp_threebody = skip_threebody_flag;
+  skip_threebody_flag = false;
   PairSW::init_style();
+  skip_threebody_flag = tmp_threebody;
 
   // adjust neighbor list request for KOKKOS
 
diff --git a/src/MANYBODY/pair_sw.cpp b/src/MANYBODY/pair_sw.cpp
index 71cb1cf5bc..de2a7ac8d6 100644
--- a/src/MANYBODY/pair_sw.cpp
+++ b/src/MANYBODY/pair_sw.cpp
@@ -14,6 +14,7 @@
 
 /* ----------------------------------------------------------------------
    Contributing author: Aidan Thompson (SNL)
+   Optimizations for two-body only: Jackson Elowitt (Univ. of Utah)
 ------------------------------------------------------------------------- */
 
 #include "pair_sw.h"
@@ -44,6 +45,7 @@ PairSW::PairSW(LAMMPS *lmp) : Pair(lmp)
   manybody_flag = 1;
   centroidstressflag = CENTROID_NOTAVAIL;
   unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
+  skip_threebody_flag = false;
 
   params = nullptr;
 
@@ -137,14 +139,18 @@ void PairSW::compute(int eflag, int vflag)
       }
 
       jtag = tag[j];
-      if (itag > jtag) {
-        if ((itag+jtag) % 2 == 0) continue;
-      } else if (itag < jtag) {
-        if ((itag+jtag) % 2 == 1) continue;
-      } else {
-        if (x[j][2] < ztmp) continue;
-        if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
-        if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
+
+      // only need to skip if we have a full neighbor list
+      if (!skip_threebody_flag) {
+        if (itag > jtag) {
+          if ((itag+jtag) % 2 == 0) continue;
+        } else if (itag < jtag) {
+          if ((itag+jtag) % 2 == 1) continue;
+        } else {
+          if (x[j][2] < ztmp) continue;
+          if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
+          if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
+        }
       }
 
       twobody(&params[ijparam],rsq,fpair,eflag,evdwl);
@@ -159,9 +165,11 @@ void PairSW::compute(int eflag, int vflag)
       if (evflag) ev_tally(i,j,nlocal,newton_pair,
                            evdwl,0.0,fpair,delx,dely,delz);
     }
-
-    jnumm1 = numshort - 1;
-
+    if (skip_threebody_flag) {
+        jnumm1 = 0;
+    } else {
+        jnumm1 = numshort - 1;
+    }
     for (jj = 0; jj < jnumm1; jj++) {
       j = neighshort[jj];
       jtype = map[type[j]];
@@ -229,9 +237,21 @@ void PairSW::allocate()
    global settings
 ------------------------------------------------------------------------- */
 
-void PairSW::settings(int narg, char **/*arg*/)
+void PairSW::settings(int narg, char ** arg)
 {
-  if (narg != 0) error->all(FLERR,"Illegal pair_style command");
+  // process optional keywords
+  int iarg = 0;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"threebody") == 0) {
+      if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "pair_style sw", error);
+      skip_threebody_flag = !utils::logical(FLERR,arg[iarg+1],false,lmp);
+      // without the threebody terms we don't need to enforce
+      // pair_coeff * * and can enable the single function.
+      one_coeff = skip_threebody_flag ? 0 : 1;
+      single_enable = skip_threebody_flag ? 1 : 0;
+      iarg += 2;
+    } else error->all(FLERR, "Illegal pair_style sw keyword: {}", arg[iarg]);
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -261,9 +281,12 @@ void PairSW::init_style()
   if (force->newton_pair == 0)
     error->all(FLERR,"Pair style Stillinger-Weber requires newton pair on");
 
-  // need a full neighbor list
+  // need a full neighbor list for full threebody calculation
 
-  neighbor->add_request(this, NeighConst::REQ_FULL);
+  if (skip_threebody_flag)
+    neighbor->add_request(this);
+  else
+    neighbor->add_request(this, NeighConst::REQ_FULL);
 }
 
 /* ----------------------------------------------------------------------
@@ -279,6 +302,19 @@ double PairSW::init_one(int i, int j)
 
 /* ---------------------------------------------------------------------- */
 
+double PairSW::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
+                      double /*factor_coul*/, double /*factor_lj*/, double &fforce)
+{
+  int ijparam = elem3param[map[itype]][map[jtype]][map[jtype]];
+  double phisw = 0.0;
+  fforce = 0.0;
+
+  if (rsq < params[ijparam].cutsq) twobody(&params[ijparam],rsq,fforce,1,phisw);
+  return phisw;
+}
+
+/* ---------------------------------------------------------------------- */
+
 void PairSW::read_file(char *file)
 {
   memory->sfree(params);
@@ -291,6 +327,8 @@ void PairSW::read_file(char *file)
     PotentialFileReader reader(lmp, file, "sw", unit_convert_flag);
     char *line;
 
+    if (skip_threebody_flag) utils::logmesg(lmp, "  disabling sw potential three-body terms\n");
+
     // transparently convert units for supported conversions
 
     int unit_convert = reader.get_unit_convert();
@@ -355,6 +393,9 @@ void PairSW::read_file(char *file)
         params[nparams].epsilon *= conversion_factor;
       }
 
+      // turn off three-body term
+      if (skip_threebody_flag) params[nparams].lambda = 0;
+
       if (params[nparams].epsilon < 0.0 || params[nparams].sigma < 0.0 ||
           params[nparams].littlea < 0.0 || params[nparams].lambda < 0.0 ||
           params[nparams].gamma < 0.0 || params[nparams].biga < 0.0 ||
diff --git a/src/MANYBODY/pair_sw.h b/src/MANYBODY/pair_sw.h
index 8bae7e5282..bcc7227554 100644
--- a/src/MANYBODY/pair_sw.h
+++ b/src/MANYBODY/pair_sw.h
@@ -32,6 +32,7 @@ class PairSW : public Pair {
   void coeff(int, char **) override;
   double init_one(int, int) override;
   void init_style() override;
+  double single(int, int, int, int, double, double, double, double &) override;
 
   static constexpr int NPARAMS_PER_LINE = 14;
 
@@ -48,10 +49,11 @@ class PairSW : public Pair {
   };
 
  protected:
-  double cutmax;      // max cutoff for all elements
-  Param *params;      // parameter set for an I-J-K interaction
-  int maxshort;       // size of short neighbor list array
-  int *neighshort;    // short neighbor list array
+  double cutmax;              // max cutoff for all elements
+  Param *params;              // parameter set for an I-J-K interaction
+  int maxshort;               // size of short neighbor list array
+  int *neighshort;            // short neighbor list array
+  int skip_threebody_flag;    // whether to run threebody loop
 
   void settings(int, char **) override;
   virtual void allocate();
diff --git a/src/MANYBODY/pair_sw_angle_table.cpp b/src/MANYBODY/pair_sw_angle_table.cpp
index d66a959d52..5b9339780f 100644
--- a/src/MANYBODY/pair_sw_angle_table.cpp
+++ b/src/MANYBODY/pair_sw_angle_table.cpp
@@ -751,3 +751,13 @@ void PairSWAngleTable::uf_lookup(ParamTable *pm, double x, double &u, double &f)
       pm->angtable->deltasq6;
   }
 }
+
+
+/* ----------------------------------------------------------------------
+   global settings
+------------------------------------------------------------------------- */
+
+void PairSWAngleTable::settings(int narg, char **/*arg*/)
+{
+  if (narg != 0) error->all(FLERR,"Illegal pair_style sw/angle/table command");
+}
diff --git a/src/MANYBODY/pair_sw_angle_table.h b/src/MANYBODY/pair_sw_angle_table.h
index 6257d18283..e8ff8f82cd 100644
--- a/src/MANYBODY/pair_sw_angle_table.h
+++ b/src/MANYBODY/pair_sw_angle_table.h
@@ -69,6 +69,7 @@ class PairSWAngleTable : public PairSW {
   void spline(double *, double *, int, double, double, double *);
   double splint(double *, double *, double *, int, double);
   void uf_lookup(ParamTable *, double, double &, double &);
+  void settings(int, char **) override;
 };
 
 }    // namespace LAMMPS_NS
diff --git a/src/MANYBODY/pair_sw_mod.cpp b/src/MANYBODY/pair_sw_mod.cpp
index ce24952fc7..e6d17b0733 100644
--- a/src/MANYBODY/pair_sw_mod.cpp
+++ b/src/MANYBODY/pair_sw_mod.cpp
@@ -41,21 +41,18 @@ PairSWMOD::PairSWMOD(LAMMPS *lmp) : PairSW(lmp)
 
 void PairSWMOD::settings(int narg, char **arg)
 {
-  // process optional keywords
-
+  // process optional keywords (and not (yet) optional keywords from parent class).
   int iarg = 0;
-
   while (iarg < narg) {
     if (strcmp(arg[iarg],"maxdelcs") == 0) {
-      if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style command");
+      if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "pair_style sw/mod", error);
       delta1 = utils::numeric(FLERR,arg[iarg+1],false,lmp);
       delta2 = utils::numeric(FLERR,arg[iarg+2],false,lmp);
       iarg += 3;
       if ((delta1 < 0.0) || (delta1 > 1.0) || (delta2 < 0.0) || (delta2 > 1.0) || (delta1 > delta2))
-        error->all(FLERR,"Illegal values for maxdelcs keyword");
-    } else error->all(FLERR,"Illegal pair_style command");
+        error->all(FLERR, "Out of range value(s) for pair style sw/mod maxdelcs keyword");
+    } else error->all(FLERR, "Illegal pair_style sw/mod keyword: {}", arg[iarg]);
   }
-  PairSW::settings(narg-iarg,arg+iarg);
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/ML-SNAP/compute_snap.cpp b/src/ML-SNAP/compute_snap.cpp
index 12a5167b2d..5d2f3666bc 100644
--- a/src/ML-SNAP/compute_snap.cpp
+++ b/src/ML-SNAP/compute_snap.cpp
@@ -14,17 +14,17 @@
 
 #include "compute_snap.h"
 
-#include "sna.h"
 #include "atom.h"
-#include "update.h"
-#include "modify.h"
-#include "neighbor.h"
-#include "neigh_list.h"
-#include "force.h"
-#include "pair.h"
 #include "comm.h"
-#include "memory.h"
 #include "error.h"
+#include "force.h"
+#include "memory.h"
+#include "modify.h"
+#include "neigh_list.h"
+#include "neighbor.h"
+#include "pair.h"
+#include "sna.h"
+#include "update.h"
 
 #include <cstring>
 
@@ -58,6 +58,7 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
   bzeroflag = 1;
   quadraticflag = 0;
   bikflag = 0;
+  dgradflag = 0;
   chemflag = 0;
   bnormflag = 0;
   wselfallflag = 0;
@@ -142,8 +143,12 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
       if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
       bikflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
       iarg += 2;
-    } else if (strcmp(arg[iarg], "switchinnerflag") == 0) {
-      if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
+    } else if (strcmp(arg[iarg],"dgradflag") == 0) {
+      if (iarg + 2 > narg) error->all(FLERR,"Illegal compute snap command");
+      dgradflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"switchinnerflag") == 0) {
+      if (iarg + 2 > narg) error->all(FLERR,"Illegal compute snap command");
       switchinnerflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
       iarg += 2;
     } else if (strcmp(arg[iarg], "sinner") == 0) {
@@ -178,6 +183,12 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
         "Illegal compute {} command: switchinnerflag = 0, unexpected sinner/dinner keyword",
         style);
 
+  if (dgradflag && !bikflag)
+    error->all(FLERR,"Illegal compute snap command: dgradflag=1 requires bikflag=1");
+
+  if (dgradflag && quadraticflag)
+    error->all(FLERR,"Illegal compute snap command: dgradflag=1 not implemented for quadratic SNAP");
+
   snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag,
                    wselfallflag, nelements, switchinnerflag);
 
@@ -194,8 +205,14 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
   natoms = atom->natoms;
   bik_rows = 1;
   if (bikflag) bik_rows = natoms;
-  size_array_rows = bik_rows+ndims_force*natoms+ndims_virial;
-  size_array_cols = nvalues*atom->ntypes+1;
+  dgrad_rows = ndims_force*natoms;
+  size_array_rows = bik_rows+dgrad_rows + ndims_virial;
+  if (dgradflag) {
+    size_array_rows = bik_rows + 3*natoms*natoms + 1;
+    size_array_cols = nvalues + 3;
+    error->warning(FLERR,"dgradflag=1 creates a N^2 array, beware of large systems.");
+  }
+  else size_array_cols = nvalues*atom->ntypes + 1;
   lastcol = size_array_cols-1;
 
   ndims_peratom = ndims_force;
@@ -305,9 +322,8 @@ void ComputeSnap::compute_array()
   // clear local peratom array
 
   for (int i = 0; i < ntotal; i++)
-    for (int icoeff = 0; icoeff < size_peratom; icoeff++) {
+    for (int icoeff = 0; icoeff < size_peratom; icoeff++)
       snap_peratom[i][icoeff] = 0.0;
-    }
 
   // invoke full neighbor list (will copy or build if necessary)
 
@@ -344,7 +360,36 @@ void ComputeSnap::compute_array()
       const int typeoffset_local = ndims_peratom*nvalues*(itype-1);
       const int typeoffset_global = nvalues*(itype-1);
 
-      // insure rij, inside, and typej  are of size jnum
+      if (dgradflag) {
+
+        // dBi/dRi tags
+
+        snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 0][0] = atom->tag[i]-1;
+        snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 0][1] = atom->tag[i]-1;
+        snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 0][2] = 0;
+        snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 1][0] = atom->tag[i]-1;
+        snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 1][1] = atom->tag[i]-1;
+        snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 1][2] = 1;
+        snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 2][0] = atom->tag[i]-1;
+        snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 2][1] = atom->tag[i]-1;
+        snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 2][2] = 2;
+
+        // dBi/dRj tags
+
+        for (int j=0; j<natoms; j++) {
+          snap[bik_rows + ((j)*3*natoms) + 3*(atom->tag[i]-1) + 0][0] = atom->tag[i]-1;
+          snap[bik_rows + ((j)*3*natoms) + 3*(atom->tag[i]-1) + 0][1] = j;
+          snap[bik_rows + ((j)*3*natoms) + 3*(atom->tag[i]-1) + 0][2] = 0;
+          snap[bik_rows + ((j)*3*natoms) + 3*(atom->tag[i]-1) + 1][0] = atom->tag[i]-1;
+          snap[bik_rows + ((j)*3*natoms) + 3*(atom->tag[i]-1) + 1][1] = j;
+          snap[bik_rows + ((j)*3*natoms) + 3*(atom->tag[i]-1) + 1][2] = 1;
+          snap[bik_rows + ((j)*3*natoms) + 3*(atom->tag[i]-1) + 2][0] = atom->tag[i]-1;
+          snap[bik_rows + ((j)*3*natoms) + 3*(atom->tag[i]-1) + 2][1] = j;
+          snap[bik_rows + ((j)*3*natoms) + 3*(atom->tag[i]-1) + 2][2] = 2;
+        }
+      }
+
+      // insure rij, inside, and typej are of size jnum
 
       snaptr->grow_rij(jnum);
 
@@ -353,7 +398,9 @@ void ComputeSnap::compute_array()
       // typej = types of neighbors of I within cutoff
       // note Rij sign convention => dU/dRij = dU/dRj = -dU/dRi
 
-      int ninside = 0;
+      // assign quantities in snaptr
+
+      int ninside=0;
       for (int jj = 0; jj < jnum; jj++) {
         int j = jlist[jj];
         j &= NEIGHMASK;
@@ -382,64 +429,54 @@ void ComputeSnap::compute_array()
         }
       }
 
+      // compute bispectrum for atom i
+
       snaptr->compute_ui(ninside, ielem);
       snaptr->compute_zi();
       snaptr->compute_bi(ielem);
 
+      // loop over neighbors for descriptors derivatives
+
       for (int jj = 0; jj < ninside; jj++) {
         const int j = snaptr->inside[jj];
+
         snaptr->compute_duidrj(jj);
         snaptr->compute_dbidrj();
 
-        // Accumulate dBi/dRi, -dBi/dRj
+        // accumulate dBi/dRi, -dBi/dRj
 
-        double *snadi = snap_peratom[i]+typeoffset_local;
-        double *snadj = snap_peratom[j]+typeoffset_local;
+        if (!dgradflag) {
 
-        for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
-          snadi[icoeff] += snaptr->dblist[icoeff][0];
-          snadi[icoeff+yoffset] += snaptr->dblist[icoeff][1];
-          snadi[icoeff+zoffset] += snaptr->dblist[icoeff][2];
-          snadj[icoeff] -= snaptr->dblist[icoeff][0];
-          snadj[icoeff+yoffset] -= snaptr->dblist[icoeff][1];
-          snadj[icoeff+zoffset] -= snaptr->dblist[icoeff][2];
-        }
+          double *snadi = snap_peratom[i]+typeoffset_local;
+          double *snadj = snap_peratom[j]+typeoffset_local;
 
-        if (quadraticflag) {
-          const int quadraticoffset = ncoeff;
-          snadi += quadraticoffset;
-          snadj += quadraticoffset;
-          int ncount = 0;
           for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
-            double bi = snaptr->blist[icoeff];
-            double bix = snaptr->dblist[icoeff][0];
-            double biy = snaptr->dblist[icoeff][1];
-            double biz = snaptr->dblist[icoeff][2];
 
-            // diagonal elements of quadratic matrix
+            snadi[icoeff] += snaptr->dblist[icoeff][0];
+            snadi[icoeff+yoffset] += snaptr->dblist[icoeff][1];
+            snadi[icoeff+zoffset] += snaptr->dblist[icoeff][2];
 
-            double dbxtmp = bi*bix;
-            double dbytmp = bi*biy;
-            double dbztmp = bi*biz;
+            snadj[icoeff] -= snaptr->dblist[icoeff][0];
+            snadj[icoeff+yoffset] -= snaptr->dblist[icoeff][1];
+            snadj[icoeff+zoffset] -= snaptr->dblist[icoeff][2];
+          }
 
-            snadi[ncount] +=         dbxtmp;
-            snadi[ncount+yoffset] += dbytmp;
-            snadi[ncount+zoffset] += dbztmp;
-            snadj[ncount] -=         dbxtmp;
-            snadj[ncount+yoffset] -= dbytmp;
-            snadj[ncount+zoffset] -= dbztmp;
+          if (quadraticflag) {
+            const int quadraticoffset = ncoeff;
+            snadi += quadraticoffset;
+            snadj += quadraticoffset;
+            int ncount = 0;
+            for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+              double bi = snaptr->blist[icoeff];
+              double bix = snaptr->dblist[icoeff][0];
+              double biy = snaptr->dblist[icoeff][1];
+              double biz = snaptr->dblist[icoeff][2];
 
-            ncount++;
+              // diagonal elements of quadratic matrix
 
-            // upper-triangular elements of quadratic matrix
-
-            for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
-              double dbxtmp = bi*snaptr->dblist[jcoeff][0]
-                + bix*snaptr->blist[jcoeff];
-              double dbytmp = bi*snaptr->dblist[jcoeff][1]
-                + biy*snaptr->blist[jcoeff];
-              double dbztmp = bi*snaptr->dblist[jcoeff][2]
-                + biz*snaptr->blist[jcoeff];
+              double dbxtmp = bi*bix;
+              double dbytmp = bi*biy;
+              double dbztmp = bi*biz;
 
               snadi[ncount] +=         dbxtmp;
               snadi[ncount+yoffset] += dbytmp;
@@ -449,60 +486,119 @@ void ComputeSnap::compute_array()
               snadj[ncount+zoffset] -= dbztmp;
 
               ncount++;
+
+              // upper-triangular elements of quadratic matrix
+
+              for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
+                double dbxtmp = bi*snaptr->dblist[jcoeff][0]
+                              + bix*snaptr->blist[jcoeff];
+                double dbytmp = bi*snaptr->dblist[jcoeff][1]
+                              + biy*snaptr->blist[jcoeff];
+                double dbztmp = bi*snaptr->dblist[jcoeff][2]
+                              + biz*snaptr->blist[jcoeff];
+
+                snadi[ncount] +=         dbxtmp;
+                snadi[ncount+yoffset] += dbytmp;
+                snadi[ncount+zoffset] += dbztmp;
+                snadj[ncount] -=         dbxtmp;
+                snadj[ncount+yoffset] -= dbytmp;
+                snadj[ncount+zoffset] -= dbztmp;
+
+                ncount++;
+              }
+            }
+          }
+        } else {
+
+            for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+
+              // add to snap array for this proc
+
+              // dBi/dRj
+
+              snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 0][icoeff+3] -= snaptr->dblist[icoeff][0];
+              snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 1][icoeff+3] -= snaptr->dblist[icoeff][1];
+              snap[bik_rows + ((atom->tag[j]-1)*3*natoms) + 3*(atom->tag[i]-1) + 2][icoeff+3] -= snaptr->dblist[icoeff][2];
+
+              // dBi/dRi
+
+              snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 0][icoeff+3] += snaptr->dblist[icoeff][0];
+              snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 1][icoeff+3] += snaptr->dblist[icoeff][1];
+              snap[bik_rows + ((atom->tag[i]-1)*3*natoms) + 3*(atom->tag[i]-1) + 2][icoeff+3] += snaptr->dblist[icoeff][2];
+            }
+          }
+      } // loop over jj inside
+
+      // accumulate Bi
+
+      if (!dgradflag) {
+
+        // linear contributions
+
+        int k = typeoffset_global;
+        for (int icoeff = 0; icoeff < ncoeff; icoeff++)
+          snap[irow][k++] += snaptr->blist[icoeff];
+
+          // quadratic contributions
+
+          if (quadraticflag) {
+            for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
+              double bveci = snaptr->blist[icoeff];
+              snap[irow][k++] += 0.5*bveci*bveci;
+              for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
+                double bvecj = snaptr->blist[jcoeff];
+                snap[irow][k++] += bveci*bvecj;
+              }
             }
           }
 
+      } else {
+          int k = 3;
+          for (int icoeff = 0; icoeff < ncoeff; icoeff++)
+          snap[irow][k++] += snaptr->blist[icoeff];
         }
-      }
-
-      // Accumulate Bi
-
-      // linear contributions
-
-      int k = typeoffset_global;
-      for (int icoeff = 0; icoeff < ncoeff; icoeff++)
-        snap[irow][k++] += snaptr->blist[icoeff];
-
-      // quadratic contributions
-
-      if (quadraticflag) {
-        for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
-          double bveci = snaptr->blist[icoeff];
-          snap[irow][k++] += 0.5*bveci*bveci;
-          for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
-            double bvecj = snaptr->blist[jcoeff];
-            snap[irow][k++] += bveci*bvecj;
-          }
-        }
-      }
     }
-  }
+  } // for (int ii = 0; ii < inum; ii++) {
 
   // accumulate bispectrum force contributions to global array
 
-  for (int itype = 0; itype < atom->ntypes; itype++) {
-    const int typeoffset_local = ndims_peratom*nvalues*itype;
-    const int typeoffset_global = nvalues*itype;
-    for (int icoeff = 0; icoeff < nvalues; icoeff++) {
-      for (int i = 0; i < ntotal; i++) {
-        double *snadi = snap_peratom[i]+typeoffset_local;
-        int iglobal = atom->tag[i];
-        int irow = 3*(iglobal-1)+bik_rows;
-        snap[irow++][icoeff+typeoffset_global] += snadi[icoeff];
-        snap[irow++][icoeff+typeoffset_global] += snadi[icoeff+yoffset];
-        snap[irow][icoeff+typeoffset_global] += snadi[icoeff+zoffset];
+  if (!dgradflag) {
+    for (int itype = 0; itype < atom->ntypes; itype++) {
+      const int typeoffset_local = ndims_peratom*nvalues*itype;
+      const int typeoffset_global = nvalues*itype;
+      for (int icoeff = 0; icoeff < nvalues; icoeff++) {
+        for (int i = 0; i < ntotal; i++) {
+          double *snadi = snap_peratom[i]+typeoffset_local;
+          int iglobal = atom->tag[i];
+          int irow = 3*(iglobal-1)+bik_rows;
+          snap[irow++][icoeff+typeoffset_global] += snadi[icoeff];
+          snap[irow++][icoeff+typeoffset_global] += snadi[icoeff+yoffset];
+          snap[irow][icoeff+typeoffset_global] += snadi[icoeff+zoffset];
+        }
       }
     }
   }
 
  // accumulate forces to global array
 
-  for (int i = 0; i < atom->nlocal; i++) {
-    int iglobal = atom->tag[i];
-    int irow = 3*(iglobal-1)+bik_rows;
-    snap[irow++][lastcol] = atom->f[i][0];
-    snap[irow++][lastcol] = atom->f[i][1];
-    snap[irow][lastcol] = atom->f[i][2];
+  if (!dgradflag) {
+    for (int i = 0; i < atom->nlocal; i++) {
+      int iglobal = atom->tag[i];
+      int irow = 3*(iglobal-1)+bik_rows;
+      snap[irow++][lastcol] = atom->f[i][0];
+      snap[irow++][lastcol] = atom->f[i][1];
+      snap[irow][lastcol] = atom->f[i][2];
+    }
+  } else {
+
+    // for dgradflag=1, put forces at first 3 columns of bik rows
+
+    for (int i=0; i<atom->nlocal; i++) {
+      int iglobal = atom->tag[i];
+      snap[iglobal-1][0+0] = atom->f[i][0];
+      snap[iglobal-1][0+1] = atom->f[i][1];
+      snap[iglobal-1][0+2] = atom->f[i][2];
+    }
   }
 
   // accumulate bispectrum virial contributions to global array
@@ -514,22 +610,34 @@ void ComputeSnap::compute_array()
   MPI_Allreduce(&snap[0][0],&snapall[0][0],size_array_rows*size_array_cols,MPI_DOUBLE,MPI_SUM,world);
 
   // assign energy to last column
-  for (int i = 0; i < bik_rows; i++) snapall[i][lastcol] = 0;
-  int irow = 0;
-  double reference_energy = c_pe->compute_scalar();
-  snapall[irow][lastcol] = reference_energy;
+
+  if (!dgradflag) {
+    for (int i = 0; i < bik_rows; i++) snapall[i][lastcol] = 0;
+    int irow = 0;
+    double reference_energy = c_pe->compute_scalar();
+    snapall[irow][lastcol] = reference_energy;
+  } else {
+
+    // assign reference energy right after the dgrad rows, first column
+
+    int irow = bik_rows + 3*natoms*natoms;
+    double reference_energy = c_pe->compute_scalar();
+    snapall[irow][0] = reference_energy;
+  }
 
   // assign virial stress to last column
   // switch to Voigt notation
 
-  c_virial->compute_vector();
-  irow += 3*natoms+bik_rows;
-  snapall[irow++][lastcol] = c_virial->vector[0];
-  snapall[irow++][lastcol] = c_virial->vector[1];
-  snapall[irow++][lastcol] = c_virial->vector[2];
-  snapall[irow++][lastcol] = c_virial->vector[5];
-  snapall[irow++][lastcol] = c_virial->vector[4];
-  snapall[irow][lastcol] = c_virial->vector[3];
+  if (!dgradflag) {
+    c_virial->compute_vector();
+    int irow = 3*natoms+bik_rows;
+    snapall[irow++][lastcol] = c_virial->vector[0];
+    snapall[irow++][lastcol] = c_virial->vector[1];
+    snapall[irow++][lastcol] = c_virial->vector[2];
+    snapall[irow++][lastcol] = c_virial->vector[5];
+    snapall[irow++][lastcol] = c_virial->vector[4];
+    snapall[irow][lastcol] = c_virial->vector[3];
+  }
 
 }
 
@@ -540,7 +648,13 @@ void ComputeSnap::compute_array()
 
 void ComputeSnap::dbdotr_compute()
 {
+
+  // no virial terms for dgrad yet
+
+  if (dgradflag) return;
+
   double **x = atom->x;
+
   int irow0 = bik_rows+ndims_force*natoms;
 
   // sum over bispectrum contributions to forces
diff --git a/src/ML-SNAP/compute_snap.h b/src/ML-SNAP/compute_snap.h
index 9f4162d6c4..0a9af58519 100644
--- a/src/ML-SNAP/compute_snap.h
+++ b/src/ML-SNAP/compute_snap.h
@@ -52,8 +52,7 @@ class ComputeSnap : public Compute {
   class SNA *snaptr;
   double cutmax;
   int quadraticflag;
-  int bikflag;
-  int bik_rows;
+  int bikflag, bik_rows, dgradflag, dgrad_rows;
 
   Compute *c_pe;
   Compute *c_virial;
diff --git a/src/ML-SNAP/sna.cpp b/src/ML-SNAP/sna.cpp
index 33937b9c45..123c16e2ac 100644
--- a/src/ML-SNAP/sna.cpp
+++ b/src/ML-SNAP/sna.cpp
@@ -1595,4 +1595,3 @@ double SNA::compute_dsfac(double r, double rcut, double sinner, double dinner)
 
   return dsfac;
 }
-
diff --git a/src/OPENMP/pair_sw_omp.cpp b/src/OPENMP/pair_sw_omp.cpp
index f6d615b2a1..4f579a9e90 100644
--- a/src/OPENMP/pair_sw_omp.cpp
+++ b/src/OPENMP/pair_sw_omp.cpp
@@ -134,14 +134,16 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
       }
 
       jtag = tag[j];
-      if (itag > jtag) {
-        if ((itag+jtag) % 2 == 0) continue;
-      } else if (itag < jtag) {
-        if ((itag+jtag) % 2 == 1) continue;
-      } else {
-        if (x[j].z < ztmp) continue;
-        if (x[j].z == ztmp && x[j].y < ytmp) continue;
-        if (x[j].z == ztmp && x[j].y == ytmp && x[j].x < xtmp) continue;
+      if (!skip_threebody_flag) {
+        if (itag > jtag) {
+          if ((itag+jtag) % 2 == 0) continue;
+        } else if (itag < jtag) {
+          if ((itag+jtag) % 2 == 1) continue;
+        } else {
+          if (x[j].z < ztmp) continue;
+          if (x[j].z == ztmp && x[j].y < ytmp) continue;
+          if (x[j].z == ztmp && x[j].y == ytmp && x[j].x < xtmp) continue;
+        }
       }
 
       twobody(&params[ijparam],rsq,fpair,EFLAG,evdwl);
@@ -156,9 +158,11 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
       if (EVFLAG) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,
                                evdwl,0.0,fpair,delx,dely,delz,thr);
     }
-
-    jnumm1 = numshort - 1;
-
+    if (skip_threebody_flag) {
+      jnumm1 = 0;
+    } else {
+      jnumm1 = numshort - 1;
+    }
     for (jj = 0; jj < jnumm1; jj++) {
       j = neighshort_thr[jj];
       jtype = map[type[j]];
diff --git a/src/REAXFF/fix_reaxff_species.cpp b/src/REAXFF/fix_reaxff_species.cpp
index a89cd6f0dc..e5298bdc4d 100644
--- a/src/REAXFF/fix_reaxff_species.cpp
+++ b/src/REAXFF/fix_reaxff_species.cpp
@@ -14,11 +14,13 @@
 /* ----------------------------------------------------------------------
    Contributing authors: Ray Shan (Sandia, tnshan@sandia.gov)
                          Oleg Sergeev (VNIIA, sergeev@vniia.ru)
+                         Jacob Gissinger (NASA, jacob.r.gissinger@gmail.com), 'delete' keyword
 ------------------------------------------------------------------------- */
 
 #include "fix_reaxff_species.h"
 
 #include "atom.h"
+#include "atom_vec.h"
 #include "comm.h"
 #include "domain.h"
 #include "error.h"
@@ -45,7 +47,8 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : Fix(lmp,
 {
   if (narg < 7) error->all(FLERR, "Illegal fix reaxff/species command");
 
-  force_reneighbor = 0;
+  force_reneighbor = 1;
+  next_reneighbor = -1;
 
   vector_flag = 1;
   size_vector = 2;
@@ -125,8 +128,10 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : Fix(lmp,
   MolName = nullptr;
   MolType = nullptr;
   NMol = nullptr;
+  Mol2Spec = nullptr;
   nd = nullptr;
   molmap = nullptr;
+  mark = nullptr;
 
   nmax = 0;
   setupflag = 0;
@@ -142,10 +147,11 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : Fix(lmp,
 
   // optional args
   eletype = nullptr;
-  ele = filepos = nullptr;
+  ele = filepos = filedel = nullptr;
   eleflag = posflag = padflag = 0;
+  delflag = specieslistflag = masslimitflag = 0;
 
-  singlepos_opened = multipos_opened = 0;
+  singlepos_opened = multipos_opened = del_opened = 0;
   multipos = 0;
   posfreq = 0;
 
@@ -180,6 +186,44 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : Fix(lmp,
       eleflag = 1;
       iarg += ntypes + 1;
 
+      // delete species
+    } else if (strcmp(arg[iarg],"delete") == 0) {
+      delflag = 1;
+      filedel = new char[255];
+      strcpy(filedel,arg[iarg+1]);
+      if (me == 0) {
+        fdel = fopen(filedel, "w");
+        if (fdel == nullptr) error->one(FLERR,"Cannot open fix reaxff/species delete file");
+      }
+
+      del_opened = 1;
+
+      if (strcmp(arg[iarg+2],"masslimit") == 0) {
+        if (iarg+5 > narg) error->all(FLERR,"Illegal fix reaxff/species command");
+        masslimitflag = 1;
+        massmin = atof(arg[iarg+3]);
+        massmax = atof(arg[iarg+4]);
+        iarg += 5;
+      } else if (strcmp(arg[iarg+2],"specieslist") == 0) {
+        specieslistflag = 1;
+        ndelspec = atoi(arg[iarg+3]);
+        if (iarg+ndelspec+4 > narg) error->all(FLERR,"Illegal fix reaxff/species command");
+
+        del_species.resize(ndelspec);
+        for (int i = 0; i < ndelspec; i ++)
+          del_species[i] = arg[iarg+4+i];
+
+        if (me == 0) {
+          fprintf(fdel,"Timestep");
+          for (i = 0; i < ndelspec; i++)
+            fprintf(fdel,"\t%s",del_species[i].c_str());
+          fprintf(fdel,"\n");
+          fflush(fdel);
+        }
+
+        iarg += ndelspec + 4;
+      } else error->all(FLERR, "Illegal fix reaxff/species command");
+
       // position of molecules
     } else if (strcmp(arg[iarg], "position") == 0) {
       if (iarg + 3 > narg) error->all(FLERR, "Illegal fix reaxff/species command");
@@ -213,6 +257,9 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : Fix(lmp,
     if (ntypes > 3) ele[3] = 'N';
   }
 
+  if (delflag && specieslistflag && masslimitflag)
+    error->all(FLERR, "Illegal fix reaxff/species command");
+
   vector_nmole = 0;
   vector_nspec = 0;
 }
@@ -229,6 +276,7 @@ FixReaxFFSpecies::~FixReaxFFSpecies()
   memory->destroy(nd);
   memory->destroy(Name);
   memory->destroy(NMol);
+  memory->destroy(Mol2Spec);
   memory->destroy(MolType);
   memory->destroy(MolName);
 
@@ -358,6 +406,8 @@ void FixReaxFFSpecies::Output_ReaxFF_Bonds(bigint ntimestep, FILE * /*fp*/)
     if (me == 0) fflush(pos);
   }
 
+  if (delflag) DeleteSpecies(Nmole, Nspec);
+
   nvalid += nfreq;
 }
 
@@ -540,6 +590,11 @@ void FixReaxFFSpecies::FindSpecies(int Nmole, int &Nspec)
   memory->create(Nameall, ntypes, "reaxff/species:Nameall");
   memory->create(NMolall, Nmole, "reaxff/species:NMolall");
 
+  memory->destroy(Mol2Spec);
+  Mol2Spec = nullptr;
+  memory->create(Mol2Spec, Nmole, "reaxff/species:Mol2Spec");
+  for (m = 0; m < Nmole; m++) Mol2Spec[m] = -1;
+
   for (m = 1, Nspec = 0; m <= Nmole; m++) {
     for (n = 0; n < ntypes; n++) Name[n] = 0;
     for (n = 0, flag_mol = 0; n < nlocal; n++) {
@@ -563,11 +618,15 @@ void FixReaxFFSpecies::FindSpecies(int Nmole, int &Nspec)
         flag_spec = 0;
         for (l = 0; l < ntypes; l++)
           if (MolName[ntypes * k + l] != Name[l]) flag_spec = 1;
-        if (flag_spec == 0) NMol[k]++;
+        if (flag_spec == 0) {
+          NMol[k]++;
+          Mol2Spec[m-1] = k;
+        }
         flag_identity *= flag_spec;
       }
       if (Nspec == 0 || flag_identity == 1) {
         for (l = 0; l < ntypes; l++) MolName[ntypes * Nspec + l] = Name[l];
+        Mol2Spec[m-1] = Nspec;
         Nspec++;
       }
     }
@@ -758,6 +817,169 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec)
 
 /* ---------------------------------------------------------------------- */
 
+void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec)
+{
+  int i, j, m, n, k, itype, cid;
+  int ndel, ndelone, count, count_tmp;
+  int *Nameall;
+  int *mask = atom->mask;
+  double localmass, totalmass;
+  double **spec_atom = f_SPECBOND->array_atom;
+  std::string species_str;
+
+  AtomVec *avec = atom->avec;
+
+  mark = nullptr;
+  memory->create(mark, nlocal, "reaxff/species:mark");
+  for (i = 0; i < nlocal; i++) mark[i] = 0;
+
+  Nameall = nullptr;
+  memory->create(Nameall, ntypes, "reaxff/species:Nameall");
+
+  int ndelcomm;
+  if (masslimitflag) ndelcomm = Nspec;
+  else ndelcomm = ndelspec;
+
+  double *deletecount;
+  memory->create(deletecount, ndelcomm, "reaxff/species:deletecount");
+  for (i = 0; i < ndelcomm; i++) deletecount[i] = 0;
+
+  int nmarklist;
+  int *marklist;
+  memory->create(marklist, nlocal, "reaxff/species:marklist");
+
+  for (m = 1; m <= Nmole; m++) {
+    localmass = totalmass = count = nmarklist = 0;
+    for (n = 0; n < ntypes; n++) Name[n] = 0;
+
+    for (i = 0; i < nlocal; i++) {
+      if (!(mask[i] & groupbit)) continue;
+      cid = nint(clusterID[i]);
+      if (cid == m) {
+        itype = atom->type[i]-1;
+        Name[itype]++;
+        count++;
+        marklist[nmarklist++] = i;
+        localmass += atom->mass[atom->type[i]];
+      }
+    }
+
+    MPI_Allreduce(&count, &count_tmp, 1, MPI_INT, MPI_SUM, world);
+    count = count_tmp;
+
+    MPI_Allreduce(Name, Nameall, ntypes, MPI_INT, MPI_SUM, world);
+    for (n = 0; n < ntypes; n++) Name[n] = Nameall[n];
+
+    MPI_Allreduce(&localmass, &totalmass, 1 , MPI_DOUBLE, MPI_SUM, world);
+
+    species_str = "";
+    for (j = 0; j < ntypes; j++) {
+      if (Name[j] != 0) {
+        if (eletype) species_str += eletype[j];
+        else species_str += ele[j];
+        if (Name[j] != 1) species_str += fmt::format("{}",Name[j]);
+      }
+    }
+
+    if (masslimitflag) {
+
+      // find corresponding moltype
+
+      if (totalmass > massmin && totalmass < massmax) {
+        for (j = 0; j < nmarklist; j++) {
+          mark[marklist[j]] = 1;
+          deletecount[Mol2Spec[m-1]] += 1.0 / (double) count;
+        }
+      }
+    } else {
+      if (count > 0) {
+        for (i = 0; i < ndelspec; i++) {
+          if (del_species[i] == species_str) {
+            for (j = 0; j < nmarklist; j++) {
+              mark[marklist[j]] = 1;
+              deletecount[i] += 1.0 / (double) count;
+            }
+            break;
+          }
+        }
+      }
+    }
+  }
+
+  // delete atoms. loop in reverse order to avoid copying marked atoms
+
+  ndel = ndelone = 0;
+  for (i = atom->nlocal-1; i >= 0; i--) {
+    if (mark[i] == 1) {
+      avec->copy(atom->nlocal-1,i,1);
+      atom->nlocal--;
+      ndelone++;
+    }
+  }
+
+  MPI_Allreduce(&ndelone, &ndel, 1, MPI_INT, MPI_SUM, world);
+
+  atom->natoms -= ndel;
+
+  if (me == 0) MPI_Reduce(MPI_IN_PLACE, deletecount, ndelcomm, MPI_DOUBLE, MPI_SUM, 0, world);
+  else MPI_Reduce(deletecount, deletecount, ndelcomm, MPI_DOUBLE, MPI_SUM, 0, world);
+
+  if (me == 0) {
+    if (masslimitflag) {
+      int printflag = 0;
+      for (int m = 0; m < Nspec; m++) {
+        if (deletecount[m] > 0) {
+          if (printflag == 0) {
+            fprintf(fdel, "Timestep %lld", update->ntimestep);
+            printflag = 1;
+          }
+          fprintf(fdel, " %g ", deletecount[m]);
+          for (j = 0; j < ntypes; j ++) {
+            int itemp = MolName[ntypes * m + j];
+            if (itemp != 0) {
+              if (eletype) fprintf(fdel, "%s", eletype[j]);
+              else fprintf(fdel, "%c", ele[j]);
+              if (itemp != 1) fprintf(fdel, "%d", itemp);
+            }
+          }
+        }
+      }
+      if (printflag) {
+        fprintf(fdel, "\n");
+        fflush(fdel);
+      }
+    } else {
+      int writeflag = 0;
+      for (i = 0; i < ndelspec; i++)
+        if (deletecount[i]) writeflag = 1;
+
+      if (writeflag) {
+        fprintf(fdel, "%lld", update->ntimestep);
+        for (i = 0; i < ndelspec; i++) {
+          fprintf(fdel, "\t%g", deletecount[i]);
+        }
+        fprintf(fdel, "\n");
+        fflush(fdel);
+      }
+    }
+  }
+
+  if (ndel && (atom->map_style != Atom::MAP_NONE)) {
+    atom->nghost = 0;
+    atom->map_init();
+    atom->map_set();
+  }
+
+  next_reneighbor = update->ntimestep;
+
+  memory->destroy(Nameall);
+  memory->destroy(marklist);
+  memory->destroy(mark);
+  memory->destroy(deletecount);
+}
+
+/* ---------------------------------------------------------------------- */
+
 double FixReaxFFSpecies::compute_vector(int n)
 {
   if (n == 0) return vector_nmole;
diff --git a/src/REAXFF/fix_reaxff_species.h b/src/REAXFF/fix_reaxff_species.h
index f441c3bc92..e272ff7d6c 100644
--- a/src/REAXFF/fix_reaxff_species.h
+++ b/src/REAXFF/fix_reaxff_species.h
@@ -45,19 +45,24 @@ class FixReaxFFSpecies : public Fix {
 
  protected:
   int me, nprocs, nmax, nlocal, ntypes, ntotal;
-  int nrepeat, nfreq, posfreq, compressed;
+  int nrepeat, nfreq, posfreq, compressed, ndelspec;
   int Nmoltype, vector_nmole, vector_nspec;
-  int *Name, *MolName, *NMol, *nd, *MolType, *molmap;
+  int *Name, *MolName, *NMol, *nd, *MolType, *molmap, *mark;
+  int *Mol2Spec;
   double *clusterID;
   AtomCoord *x0;
 
   double bg_cut;
   double **BOCut;
 
-  FILE *fp, *pos;
+  std::vector<std::string> del_species;
+
+  FILE *fp, *pos, *fdel;
   int eleflag, posflag, multipos, padflag, setupflag;
-  int singlepos_opened, multipos_opened;
-  char *ele, **eletype, *filepos;
+  int delflag, specieslistflag, masslimitflag;
+  double massmin, massmax;
+  int singlepos_opened, multipos_opened, del_opened;
+  char *ele, **eletype, *filepos, *filedel;
 
   void Output_ReaxFF_Bonds(bigint, FILE *);
   AtomCoord chAnchor(AtomCoord, AtomCoord);
@@ -65,6 +70,7 @@ class FixReaxFFSpecies : public Fix {
   void SortMolecule(int &);
   void FindSpecies(int, int &);
   void WriteFormulas(int, int);
+  void DeleteSpecies(int, int);
   int CheckExistence(int, int);
 
   int nint(const double &);
diff --git a/src/input.cpp b/src/input.cpp
index c67b715bcd..2979503d30 100644
--- a/src/input.cpp
+++ b/src/input.cpp
@@ -188,17 +188,20 @@ of the file is reached.  The *infile* pointer will usually point to
 
 void Input::file()
 {
-  int m,n;
+  int m,n,mstart,ntriple,endfile;
 
   while (true) {
 
     // read a line from input script
-    // n = length of line including str terminator, 0 if end of file
+    // when done, n = length of line including str terminator, 0 if end of file
     // if line ends in continuation char '&', concatenate next line
+    // if triple quotes are used, read until closing triple quotes
 
     if (me == 0) {
-
+      ntriple = 0;
+      endfile = 0;
       m = 0;
+
       while (true) {
 
         if (infile == nullptr) {
@@ -206,38 +209,58 @@ void Input::file()
           break;
         }
 
-        if (maxline-m < 2) reallocate(line,maxline,0);
+        mstart = m;
 
-        // end of file reached, so break
-        // n == 0 if nothing read, else n = line with str terminator
+        while (1) {
+          if (maxline-m < 2) reallocate(line,maxline,0);
 
-        if (fgets(&line[m],maxline-m,infile) == nullptr) {
-          if (m) n = strlen(line) + 1;
-          else n = 0;
+          // end of file reached, so break
+          // n == 0 if nothing read, else n = line with str terminator
+
+          if (fgets(&line[m],maxline-m,infile) == nullptr) {
+            endfile = 1;
+            if (m) n = strlen(line) + 1;
+            else n = 0;
+            break;
+          }
+
+          // continue if last char read was not a newline
+          // can happen if line is very long
+
+          m += strlen(&line[m]);
+          if (line[m-1] != '\n') continue;
           break;
         }
 
-        // continue if last char read was not a newline
-        // could happen if line is very long
+        if (endfile) break;
 
-        m = strlen(line);
-        if (line[m-1] != '\n') continue;
+        // add # of triple quotes in just-read line to ntriple
 
-        // continue reading if final printable char is & char
-        // or if odd number of triple quotes
-        // else break with n = line with str terminator
+        ntriple += numtriple(&line[mstart]);
+
+        // trim whitespace from end of line
+        // line[m] = last printable char
 
         m--;
         while (m >= 0 && isspace(line[m])) m--;
-        if (m < 0 || line[m] != '&') {
-          if (numtriple(line) % 2) {
-            m += 2;
-            continue;
-          }
-          line[m+1] = '\0';
-          n = m+2;
-          break;
+
+        // continue reading if final printable char is "&"
+
+        if (m >= 0 && line[m] == '&') continue;
+
+        // continue reading if odd number of triple quotes
+
+        if (ntriple % 2) {
+          line[m+1] = '\n';
+          m += 2;
+          continue;
         }
+
+        // done, break with n = length of line with str terminator
+
+        line[m+1] = '\0';
+        n = m+2;
+        break;
       }
     }
 
@@ -371,8 +394,9 @@ char *Input::one(const std::string &single)
 }
 
 /* ----------------------------------------------------------------------
-   Send text to active echo file pointers
+   send text to active echo file pointers
 ------------------------------------------------------------------------- */
+
 void Input::write_echo(const std::string &txt)
 {
   if (me == 0) {
@@ -399,34 +423,35 @@ void Input::parse()
   if (n > maxcopy) reallocate(copy,maxcopy,n);
   strcpy(copy,line);
 
-  // strip any # comment by replacing it with 0
-  // do not strip from a # inside single/double/triple quotes
-  // quoteflag = 1,2,3 when encounter first single/double,triple quote
-  // quoteflag = 0 when encounter matching single/double,triple quote
+  // strip a # comment by replacing it with 0
+  // do not treat a # inside single/double/triple quotes as a comment
 
-  int quoteflag = 0;
+  char *ptrmatch;
   char *ptr = copy;
+
   while (*ptr) {
-    if (*ptr == '#' && !quoteflag) {
+    if (*ptr == '#') {
       *ptr = '\0';
       break;
     }
-    if (quoteflag == 0) {
+    if (*ptr == '\'') {
+      ptrmatch = strchr(ptr+1,'\'');
+      if (ptrmatch == NULL)
+        error->all(FLERR,"Unmatched single quote in command");
+      ptr = ptrmatch + 1;
+    } else if (*ptr == '"') {
       if (strstr(ptr,"\"\"\"") == ptr) {
-        quoteflag = 3;
-        ptr += 2;
+        ptrmatch = strstr(ptr+3,"\"\"\"");
+        if (ptrmatch == NULL)
+          error->all(FLERR,"Unmatched triple quote in command");
+        ptr = ptrmatch + 3;
+      } else {
+        ptrmatch = strchr(ptr+1,'"');
+        if (ptrmatch == NULL)
+          error->all(FLERR,"Unmatched double quote in command");
+        ptr = ptrmatch + 1;
       }
-      else if (*ptr == '"') quoteflag = 2;
-      else if (*ptr == '\'') quoteflag = 1;
-    } else {
-      if (quoteflag == 3 && strstr(ptr,"\"\"\"") == ptr) {
-        quoteflag = 0;
-        ptr += 2;
-      }
-      else if (quoteflag == 2 && *ptr == '"') quoteflag = 0;
-      else if (quoteflag == 1 && *ptr == '\'') quoteflag = 0;
-    }
-    ptr++;
+    } else ptr++;
   }
 
   if (utils::has_utf8(copy)) {
@@ -534,16 +559,18 @@ void Input::substitute(char *&str, char *&str2, int &max, int &max2, int flag)
 {
   // use str2 as scratch space to expand str, then copy back to str
   // reallocate str and str2 as necessary
-  // do not replace $ inside single/double/triple quotes
+  // do not replace variables inside single/double/triple quotes
   // var = pts at variable name, ended by null char
   //   if $ is followed by '{', trailing '}' becomes null char
   //   else $x becomes x followed by null char
   // beyond = points to text following variable
 
-  int i,n,paren_count;
+  int i,n,paren_count,nchars;;
   char immediate[256];
   char *var,*value,*beyond;
   int quoteflag = 0;
+  char *ptrmatch;
+
   char *ptr = str;
 
   n = strlen(str) + 1;
@@ -637,34 +664,41 @@ void Input::substitute(char *&str, char *&str2, int &max, int &max2, int flag)
         if (echo_log && logfile) fprintf(logfile,"%s%s\n",str2,beyond);
       }
 
-      continue;
-    }
+    // check for single/double/triple quotes and skip past them
 
-    // quoteflag = 1,2,3 when encounter first single/double,triple quote
-    // quoteflag = 0 when encounter matching single/double,triple quote
-    // copy 2 extra triple quote chars into str2
-
-    if (quoteflag == 0) {
+    } else if (*ptr == '\'') {
+      ptrmatch = strchr(ptr+1,'\'');
+      if (ptrmatch == NULL)
+        error->all(FLERR,"Unmatched single quote in command");
+      nchars = ptrmatch+1 - ptr;
+      strncpy(ptr2,ptr,nchars);
+      ptr += nchars;
+      ptr2 += nchars;
+    } else if (*ptr == '"') {
       if (strstr(ptr,"\"\"\"") == ptr) {
-        quoteflag = 3;
-        *ptr2++ = *ptr++;
-        *ptr2++ = *ptr++;
+        ptrmatch = strstr(ptr+3,"\"\"\"");
+        if (ptrmatch == NULL)
+          error->all(FLERR,"Unmatched triple quote in command");
+        nchars = ptrmatch+3 - ptr;
+        strncpy(ptr2,ptr,nchars);
+        ptr += nchars;
+        ptr2 += nchars;
+      } else {
+        ptrmatch = strchr(ptr+1,'"');
+        if (ptrmatch == NULL)
+          error->all(FLERR,"Unmatched double quote in command");
+        nchars = ptrmatch+1 - ptr;
+        strncpy(ptr2,ptr,nchars);
+        ptr += nchars;
+        ptr2 += nchars;
       }
-      else if (*ptr == '"') quoteflag = 2;
-      else if (*ptr == '\'') quoteflag = 1;
-    } else {
-      if (quoteflag == 3 && strstr(ptr,"\"\"\"") == ptr) {
-        quoteflag = 0;
-        *ptr2++ = *ptr++;
-        *ptr2++ = *ptr++;
-      }
-      else if (quoteflag == 2 && *ptr == '"') quoteflag = 0;
-      else if (quoteflag == 1 && *ptr == '\'') quoteflag = 0;
-    }
 
-    // copy current character into str2
+    // else copy current single character into str2
+
+    } else *ptr2++ = *ptr++;
+
+    // terminate current str2 so variable sub can perform strlen()
 
-    *ptr2++ = *ptr++;
     *ptr2 = '\0';
   }
 
diff --git a/src/lammps.cpp b/src/lammps.cpp
index 3ddfafeb5d..e78f9e2019 100644
--- a/src/lammps.cpp
+++ b/src/lammps.cpp
@@ -442,8 +442,7 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) :
   // sum of procs in all worlds must equal total # of procs
 
   if (!universe->consistent())
-    error->universe_all(FLERR,"Processor partitions do not match "
-                        "number of allocated processors");
+    error->universe_all(FLERR,"Processor partitions do not match number of allocated processors");
 
   // universe cannot use stdin for input file
 
@@ -512,10 +511,15 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) :
       else infile = fopen(arg[inflag],"r");
       if (infile == nullptr)
         error->one(FLERR,"Cannot open input script {}: {}", arg[inflag], utils::getsyserror());
+      if (!helpflag)
+        utils::logmesg(this,fmt::format("LAMMPS ({}{})\n",version,UPDATE_STRING));
+      // warn against using I/O redirection in parallel runs
+      if ((inflag == 0) && (universe->nprocs > 1))
+        error->warning(FLERR, "Using I/O redirection is unreliable with parallel runs. "
+                       "Better use -in switch to read input file.");
+      utils::flush_buffers(this);
     }
 
-    if ((universe->me == 0) && !helpflag)
-      utils::logmesg(this,fmt::format("LAMMPS ({}{})\n",version,UPDATE_STRING));
 
   // universe is one or more worlds, as setup by partition switch
   // split universe communicator into separate world communicators
diff --git a/src/read_dump.cpp b/src/read_dump.cpp
index 8af19ee601..aa9ec0bbf0 100644
--- a/src/read_dump.cpp
+++ b/src/read_dump.cpp
@@ -121,7 +121,7 @@ void ReadDump::command(int narg, char **arg)
 
   // reset timestep to nstep
 
-  update->reset_timestep(nstep, true);
+  if (timestepflag) update->reset_timestep(nstep, true);
 
   // counters
 
@@ -1202,6 +1202,7 @@ int ReadDump::fields_and_keywords(int narg, char **arg)
 
   multiproc_nfile = 0;
   boxflag = 1;
+  timestepflag = 1;
   replaceflag = 1;
   purgeflag = 0;
   trimflag = 0;
@@ -1219,6 +1220,10 @@ int ReadDump::fields_and_keywords(int narg, char **arg)
       if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
       boxflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
       iarg += 2;
+    } else if (strcmp(arg[iarg],"timestep") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
+      timestepflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
+      iarg += 2;
     } else if (strcmp(arg[iarg],"replace") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
       replaceflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
diff --git a/src/read_dump.h b/src/read_dump.h
index 1a94cd1cd3..506cd197c1 100644
--- a/src/read_dump.h
+++ b/src/read_dump.h
@@ -63,7 +63,8 @@ class ReadDump : public Command {
   int dimension;    // same as in Domain
   int triclinic;
 
-  int boxflag;                 // overwrite simulation with dump file box params
+  int boxflag;                 // overwrite simulation box with dump file box params
+  int timestepflag;            // overwrite simulation timestep with dump file timestep
   int replaceflag, addflag;    // flags for processing dump snapshot atoms
   int trimflag, purgeflag;
   int scaleflag;        // user 0/1 if dump file coords are unscaled/scaled
diff --git a/unittest/force-styles/tests/manybody-pair-sw_twobody.yaml b/unittest/force-styles/tests/manybody-pair-sw_twobody.yaml
new file mode 100644
index 0000000000..0c48d47137
--- /dev/null
+++ b/unittest/force-styles/tests/manybody-pair-sw_twobody.yaml
@@ -0,0 +1,156 @@
+---
+lammps_version: 23 Jun 2022
+date_generated: Sat Jul  2 17:42:21 2022
+epsilon: 1e-10
+skip_tests:
+prerequisites: ! |
+  pair sw
+pre_commands: ! |
+  variable newton_pair delete
+  if "$(is_active(package,gpu)) > 0.0" then "variable newton_pair index off" else "variable newton_pair index on"
+post_commands: ! ""
+input_file: in.manybody
+pair_style: sw threebody off
+pair_coeff: ! |
+  * * Si.sw Si Si Si Si Si Si Si Si
+extract: ! ""
+natoms: 64
+init_vdwl: -258.5086400674769
+init_coul: 0
+init_stress: ! |2-
+   6.4784557236140188e+00  8.9666141338671075e+00  1.6564010213468620e+01 -4.9679217055624925e+00  3.4388959220521961e+01  7.5389343797929154e-01
+init_forces: ! |2
+    1 -7.5117950925001264e-01  3.0512035938320858e+00  1.7548369060319577e+00
+    2 -2.5839943840157944e+00 -6.2407030855132184e-01 -2.0776063043681416e+00
+    3  6.4874626526429646e-01 -1.9191097130296589e-01 -2.8953907507008203e-01
+    4 -1.8780621641443120e+00  2.9563493639001268e+00  6.0045000947756955e-01
+    5 -1.3527836952861758e+00 -7.5563316408513048e-01 -4.0171785277050770e-01
+    6  3.0261027539317009e-01  3.3010495375946016e+00  1.0578432546004755e+00
+    7 -1.0220792530722718e+00 -8.9704171544868350e-01  2.2679286915120240e+00
+    8  7.3786854090375081e-02  6.8184832716858157e-01 -9.9742395408331985e-01
+    9  1.5404185888118715e-01 -2.4234656391168983e+00 -2.5547421021459220e+00
+   10  2.9152006680754144e-01 -1.3832499077045937e+00 -2.3112314100238849e+00
+   11  3.0397941651131566e+00 -3.2696818086952382e+00  2.0402858500342536e+00
+   12 -5.2427750818963892e+00 -2.9266181850189747e+00 -2.3274218711162238e+00
+   13 -3.4764983756087076e-01  5.1759636691465873e+00 -9.5777981490571462e-01
+   14 -1.3994974411693033e+00  3.6723581318644940e+00  4.9022891744156372e-01
+   15 -3.6451950102391031e+00  4.1804831124641382e+00 -2.3094319497556559e+00
+   16  1.7762576801244847e+00 -1.7769734947711013e-01  5.6118226682874788e+00
+   17  1.1626276987569719e+00  2.5434318242406255e+00 -4.1298909446437833e+00
+   18  2.6370967308167081e-01  4.8174395544262438e-01  3.2879300848711184e+00
+   19 -1.2180668170296169e+00 -2.0218982136836656e+00 -3.8444674827692227e-01
+   20 -6.1734234441618661e+00 -8.7047552939678141e-02  7.4088977580926274e-01
+   21  2.2527604735738231e+00  9.2512650808085084e-01 -2.7596546519202407e+00
+   22 -5.0794028512678571e+00  3.2609137644938517e+00 -3.9745643191135742e+00
+   23  1.8924999787954402e+00  3.3526647080652703e+00  1.2248568956854238e+00
+   24  1.5743771798508372e+00  1.3293691279058384e+00  2.6631027667815861e+00
+   25  8.3110919566309338e-01 -1.1460141214066610e-01  1.4368370341871295e+00
+   26 -4.7991049590854340e-01 -6.7362488036409351e-01  1.2327946774343894e+00
+   27  1.9900573442863836e+00 -5.0620688348406084e-01  1.5762361423313080e+00
+   28  5.7701286079414915e-01  1.3187742370100088e+00  4.1298244042734957e+00
+   29 -3.0416482921788530e+00  9.1958118398971411e-01 -1.1418981151511567e+00
+   30 -1.5986340571260498e+00  1.2172058599377520e+00 -8.9839567436920520e-01
+   31 -1.6221367664137194e+00  1.4053388522714974e+00 -4.1030835758060596e-01
+   32 -3.3999213065003993e+00  5.4646623792746152e-01 -4.9970596852221305e-01
+   33  4.3374523080961502e+00 -2.0526192390847231e+00  2.7863621822774527e+00
+   34  3.6899966300377274e-01 -9.7647298273011718e-02  3.7849094767735275e-01
+   35  1.2474956459695217e+00 -1.6786310914928160e-01  2.5255688468039841e+00
+   36  1.9423002050777016e-02 -2.5811032787101440e+00 -5.3464409483959238e-02
+   37 -1.7991755427583054e+00  1.7326288527074167e+00 -2.3172591544605829e+00
+   38 -2.6635345675769728e-01 -4.5619472453099841e-01  2.9146619578940203e-01
+   39  2.2040425723147274e+00 -1.2990665525340543e+00 -4.1031233229148558e+00
+   40  3.8636879669801210e+00 -1.9428932057790562e+00 -1.2090029697953577e+00
+   41  6.9184862237524347e-02 -5.5563025877752470e-01 -1.0104421056432380e+00
+   42 -4.1077204842477482e+00 -1.9092260092568061e+00 -2.7778884577312546e-01
+   43  6.9234276916198878e-01 -5.1959961009882676e+00 -3.8252772226000875e-01
+   44  3.1974368019332173e+00 -3.7333126721070280e+00  1.1927602690851384e+00
+   45 -1.0979421463666958e+00  1.4281871410588294e+00 -2.7844688870636198e+00
+   46 -2.3123927283410217e-01 -1.6246499267294727e+00  4.6068188624710249e+00
+   47  2.4575638270171467e+00  1.3529111936752543e+00  8.5779610605164602e-01
+   48  1.3053149069961443e+00  5.5515699432490484e-01 -3.4671208564970402e-01
+   49 -1.6274632987918105e+00 -4.7057286351454088e+00 -2.4249279782812732e+00
+   50  5.4224455439310093e-03  2.9897586979430217e+00 -1.3950914387971693e+00
+   51 -1.2459066473548792e+00  2.9456460712366876e+00  3.7288916669729959e+00
+   52  4.2616245425615205e+00 -4.4802874559504291e+00  4.4417404910061506e+00
+   53  2.7936251596841610e+00  2.8362368635929394e+00  1.5493162393308044e+00
+   54  2.4623429186012755e+00 -2.8582750315396499e+00 -1.7184511417663617e+00
+   55  8.8879734150460621e-01 -2.9956850724190431e+00 -3.5690867805221691e+00
+   56 -4.2180992335342177e-01  2.9462851508525678e-01 -2.1026878944189669e+00
+   57  1.9534125277666281e-01 -1.1222033415899244e+00 -3.0127694733977195e-01
+   58  2.7751167493305626e+00 -9.6251009716841951e-01  1.2847140239740573e+00
+   59  2.7332664305562209e+00  1.6874001147167499e+00 -1.8630365579739607e+00
+   60 -1.9920742303178285e+00 -5.1976595137487056e+00 -2.5715712256855774e+00
+   61 -3.9186643718089481e-01  1.8941052818415300e+00 -7.2852310082079819e-01
+   62 -2.6017296078018819e+00  2.4477411514482621e+00 -2.1740532348808825e+00
+   63  1.6751774499211520e+00 -2.9786131173037549e+00 -4.4746311293685642e-02
+   64  2.2350712682686380e+00  2.4856397598318205e+00  6.0442073184431777e+00
+run_vdwl: -258.499584619321
+run_coul: 0
+run_stress: ! |2-
+   6.4308150527103614e+00  8.9800084577004764e+00  1.6644352826777482e+01 -5.0072836166618568e+00  3.4219822049720790e+01  1.0492755818511172e+00
+run_forces: ! |2
+    1 -7.5937225945191122e-01  3.0453338820624811e+00  1.7588997466292056e+00
+    2 -2.5842358982706677e+00 -6.4961465564128662e-01 -2.0836136472876197e+00
+    3  6.1692053064565644e-01 -1.7923404289185974e-01 -2.5776384969836241e-01
+    4 -1.8642434201915430e+00  2.9689111966333672e+00  5.7860631771456472e-01
+    5 -1.3707214670175785e+00 -7.6882987771621347e-01 -3.8201327706453814e-01
+    6  3.5155029154762341e-01  3.3059375310476078e+00  1.0816937673465203e+00
+    7 -1.0031791908316097e+00 -8.6956258363515204e-01  2.2425245701191070e+00
+    8  5.9471905500604189e-02  6.7712638453434015e-01 -9.9455163342941122e-01
+    9  1.3241459971784919e-01 -2.4451165559600572e+00 -2.5688790543810329e+00
+   10  2.7987760793693212e-01 -1.3869445644157263e+00 -2.2737317068259308e+00
+   11  3.0093494238363268e+00 -3.2460329093478379e+00  2.0365059309648439e+00
+   12 -5.2407947220384763e+00 -2.9126835314212682e+00 -2.3181888620032565e+00
+   13 -2.9313171694052698e-01  5.1682815539279954e+00 -9.3178794848643598e-01
+   14 -1.4205609885169097e+00  3.6728355747283463e+00  4.7973399101037340e-01
+   15 -3.6438539671591643e+00  4.1738585973197502e+00 -2.2995365264114347e+00
+   16  1.7665430094094714e+00 -1.6683625143013525e-01  5.5980371284228605e+00
+   17  1.1609229005729407e+00  2.5445052813433287e+00 -4.1101936902212186e+00
+   18  2.3250564926704054e-01  5.0119516573696155e-01  3.2999642189792304e+00
+   19 -1.2455755116140121e+00 -2.0483645027853807e+00 -3.9915953456529252e-01
+   20 -6.1533476323819203e+00 -1.0336110065476412e-01  7.2207408198501155e-01
+   21  2.2580812268964414e+00  8.8665782585823316e-01 -2.7867801730661323e+00
+   22 -5.0715437260287493e+00  3.2720805913066657e+00 -3.9870515109961557e+00
+   23  1.9067384153660658e+00  3.3666024181351721e+00  1.2360966484469924e+00
+   24  1.6076366668181861e+00  1.3129049141466476e+00  2.6481286770383776e+00
+   25  8.2849608759830151e-01 -1.0946573631083545e-01  1.4137379099564085e+00
+   26 -4.9245162995242880e-01 -6.5633188499443484e-01  1.2397189435323861e+00
+   27  2.0002068209516661e+00 -5.2635863568453822e-01  1.5812202387080725e+00
+   28  5.6469856769551852e-01  1.3419655823448202e+00  4.1390158656307525e+00
+   29 -3.0414702840012948e+00  9.2717141813720239e-01 -1.1446412926158587e+00
+   30 -1.5780449202280797e+00  1.1972994610432468e+00 -9.0937832538133612e-01
+   31 -1.6135468940121711e+00  1.4097747951848163e+00 -4.1013831792153099e-01
+   32 -3.3839861094060888e+00  5.3294244523688450e-01 -5.1150418397488095e-01
+   33  4.3179202422054406e+00 -2.0530236039826524e+00  2.7909047875298811e+00
+   34  3.6044665483098232e-01 -9.6636898593039144e-02  3.8719889716901629e-01
+   35  1.2574806392912055e+00 -1.5474474915601943e-01  2.5284696614854756e+00
+   36 -9.2507293281644375e-03 -2.5688826605251358e+00 -8.8917022692741904e-02
+   37 -1.7874734690581056e+00  1.7150271178647891e+00 -2.3271057624958509e+00
+   38 -2.4661620498076103e-01 -4.4857619626472056e-01  3.0090275233366270e-01
+   39  2.1876463817083560e+00 -1.2953665685718856e+00 -4.1070251548878520e+00
+   40  3.8954619284502536e+00 -1.9483835119233155e+00 -1.2227531289486768e+00
+   41  9.7565025130135041e-02 -5.2646053136013127e-01 -9.8280124384883183e-01
+   42 -4.1063684952413535e+00 -1.9251952021816745e+00 -2.9901720065849435e-01
+   43  6.7939692714052213e-01 -5.1874569217280806e+00 -3.9562141807204243e-01
+   44  3.2038887558507829e+00 -3.7312345150786044e+00  1.1791192854596575e+00
+   45 -1.1219574440264417e+00  1.4017870123386482e+00 -2.7737869798443779e+00
+   46 -2.3209294883861276e-01 -1.6182897275443398e+00  4.6296975809710883e+00
+   47  2.4602350208971560e+00  1.3452019899041738e+00  8.6358301884597988e-01
+   48  1.3093376440414450e+00  5.7371015819685922e-01 -3.5578408564713881e-01
+   49 -1.6696474765985014e+00 -4.7474039477977401e+00 -2.4607694981777248e+00
+   50  1.2989593641085595e-02  2.9812087098985032e+00 -1.4123464138675532e+00
+   51 -1.2853374902017087e+00  2.9768731587433548e+00  3.7403754374802158e+00
+   52  4.2768627438095859e+00 -4.4665207635055904e+00  4.4526942886426610e+00
+   53  2.8263576162367916e+00  2.8711821643474851e+00  1.6198032393140254e+00
+   54  2.4557393026757914e+00 -2.8576870371250496e+00 -1.7106753514157151e+00
+   55  8.7016650440839127e-01 -2.9993943594233912e+00 -3.5406400607258477e+00
+   56 -4.2720269240408387e-01  2.6284724466018883e-01 -2.0965871396618168e+00
+   57  1.8938909101993642e-01 -1.1500277319075947e+00 -3.0413618657652775e-01
+   58  2.7772395410738167e+00 -9.4351277181421578e-01  1.2588196612829914e+00
+   59  2.7573808165218825e+00  1.7074121472099479e+00 -1.8649938272758655e+00
+   60 -2.0326029067319982e+00 -5.2224816026297454e+00 -2.6126651787431214e+00
+   61 -3.8509493120844207e-01  1.9073567822120387e+00 -7.1625738088539748e-01
+   62 -2.6004175741382847e+00  2.4418292628003133e+00 -2.1615478367945307e+00
+   63  1.6782508781720011e+00 -3.0101538006831734e+00 -7.1986308169233931e-02
+   64  2.2749536899334073e+00  2.5303495677814198e+00  6.0668040667204064e+00
+...
diff --git a/unittest/force-styles/tests/manybody-pair-sw_twothree.yaml b/unittest/force-styles/tests/manybody-pair-sw_twothree.yaml
new file mode 100644
index 0000000000..0a7c60a892
--- /dev/null
+++ b/unittest/force-styles/tests/manybody-pair-sw_twothree.yaml
@@ -0,0 +1,157 @@
+---
+lammps_version: 23 Jun 2022
+date_generated: Sat Jul  2 18:29:21 2022
+epsilon: 1e-10
+skip_tests: gpu
+prerequisites: ! |
+  pair sw
+pre_commands: ! |
+  variable newton_pair delete
+  variable newton_pair index on
+post_commands: ! ""
+input_file: in.manybody
+pair_style: hybrid sw threebody on sw threebody off
+pair_coeff: ! |
+  * *   sw 1 Si.sw Si Si Si Si NULL NULL NULL NULL
+  * 5*8 sw 2 Si.sw Si Si Si Si Si Si Si Si
+extract: ! ""
+natoms: 64
+init_vdwl: -258.5086399744112
+init_coul: 0
+init_stress: ! |2-
+   6.4784606418926112e+00  8.9666260111962721e+00  1.6564019181847929e+01 -4.9679178424853108e+00  3.4388958502311766e+01  7.5390085880741020e-01
+init_forces: ! |2
+    1 -7.5117950925001176e-01  3.0512035938320858e+00  1.7548369060319566e+00
+    2 -2.5839943840157309e+00 -6.2407030855217227e-01 -2.0776063043672743e+00
+    3  6.4874626526431012e-01 -1.9191097130277632e-01 -2.8953907507029547e-01
+    4 -1.8780621642391011e+00  2.9563493639933633e+00  6.0045000947943949e-01
+    5 -1.3527836952861756e+00 -7.5563316408513048e-01 -4.0171785277050776e-01
+    6  3.0261027539317009e-01  3.3010495375946016e+00  1.0578432546004755e+00
+    7 -1.0220792530722718e+00 -8.9704171544868350e-01  2.2679286915120240e+00
+    8  7.3786854090375081e-02  6.8184832716858157e-01 -9.9742395408331985e-01
+    9  1.5404185888294175e-01 -2.4234656392100611e+00 -2.5547421022333259e+00
+   10  2.9152004537455556e-01 -1.3832499266817839e+00 -2.3112314079357352e+00
+   11  3.0397941652061906e+00 -3.2696818086953119e+00  2.0402858501197878e+00
+   12 -5.2427766862019523e+00 -2.9266230951551324e+00 -2.3274249447095405e+00
+   13 -3.4764983756087076e-01  5.1759636691465873e+00 -9.5777981490571462e-01
+   14 -1.3994974411693033e+00  3.6723581318644940e+00  4.9022891744156372e-01
+   15 -3.6451950102391031e+00  4.1804831124641391e+00 -2.3094319497556559e+00
+   16  1.7762576801244845e+00 -1.7769734947711072e-01  5.6118226682874788e+00
+   17  1.1626295956924018e+00  2.5434335885879951e+00 -4.1298909753889967e+00
+   18  2.6370967308167081e-01  4.8174395544262394e-01  3.2879300848711188e+00
+   19 -1.2180668170296169e+00 -2.0218982136836656e+00 -3.8444674827692210e-01
+   20 -6.1734234441618661e+00 -8.7047552939678141e-02  7.4088977580926274e-01
+   21  2.2527604735738231e+00  9.2512650808085084e-01 -2.7596546519202407e+00
+   22 -5.0794028512678571e+00  3.2609137644938517e+00 -3.9745643191135742e+00
+   23  1.8924999787954402e+00  3.3526647080652703e+00  1.2248568956854238e+00
+   24  1.5743771798508370e+00  1.3293691279058382e+00  2.6631027667815861e+00
+   25  8.3110919554353102e-01 -1.1460141213930673e-01  1.4368370342964809e+00
+   26 -4.7991049590854340e-01 -6.7362488036409351e-01  1.2327946774343894e+00
+   27  1.9900570741156840e+00 -5.0620387120655386e-01  1.5762394195123659e+00
+   28  5.7701286079972869e-01  1.3187742380597833e+00  4.1298244053155067e+00
+   29 -3.0416482921788530e+00  9.1958118398971411e-01 -1.1418981151511567e+00
+   30 -1.5986340571260498e+00  1.2172058599377520e+00 -8.9839567436920520e-01
+   31 -1.6221367664137194e+00  1.4053388522714974e+00 -4.1030835758060596e-01
+   32 -3.3999213065003993e+00  5.4646623792746152e-01 -4.9970596852221305e-01
+   33  4.3374523080961529e+00 -2.0526192390847209e+00  2.7863621822774500e+00
+   34  3.6899967841270293e-01 -9.7647299317203395e-02  3.7849093017581542e-01
+   35  1.2474956459695270e+00 -1.6786310914897268e-01  2.5255688468043052e+00
+   36  1.9423002050769304e-02 -2.5811032787101440e+00 -5.3464409483951571e-02
+   37 -1.7991755427583054e+00  1.7326288527074167e+00 -2.3172591544605829e+00
+   38 -2.6635345675769728e-01 -4.5619472453099841e-01  2.9146619578940203e-01
+   39  2.2040425723147274e+00 -1.2990665525340543e+00 -4.1031233229148558e+00
+   40  3.8636879669801210e+00 -1.9428932057790562e+00 -1.2090029697953577e+00
+   41  6.9184862237524403e-02 -5.5563025877752559e-01 -1.0104421056432389e+00
+   42 -4.1077204863643084e+00 -1.9092260113246491e+00 -2.7778884595397102e-01
+   43  6.9234276922332216e-01 -5.1959961010562941e+00 -3.8252772225624965e-01
+   44  3.1974368019347907e+00 -3.7333126721157721e+00  1.1927602690968555e+00
+   45 -1.0979421463666958e+00  1.4281871410588294e+00 -2.7844688870636198e+00
+   46 -2.3123927283410217e-01 -1.6246499267294727e+00  4.6068188624710249e+00
+   47  2.4575638270171472e+00  1.3529111936752540e+00  8.5779610605164602e-01
+   48  1.3053149069961458e+00  5.5515699432490606e-01 -3.4671208564970524e-01
+   49 -1.6274632987918134e+00 -4.7057286351454106e+00 -2.4249279782812714e+00
+   50  5.4224455439262353e-03  2.9897586979430182e+00 -1.3950914387971656e+00
+   51 -1.2459066473548797e+00  2.9456460712366872e+00  3.7288916669729959e+00
+   52  4.2616245425615027e+00 -4.4802874559509265e+00  4.4417404910060432e+00
+   53  2.7936251596841610e+00  2.8362368635929394e+00  1.5493162393308042e+00
+   54  2.4623429186012760e+00 -2.8582750315396499e+00 -1.7184511417663617e+00
+   55  8.8879734150460621e-01 -2.9956850724190431e+00 -3.5690867805221691e+00
+   56 -4.2180992335342188e-01  2.9462851508525667e-01 -2.1026878944189669e+00
+   57  1.9534125277666264e-01 -1.1222033415899244e+00 -3.0127694733977195e-01
+   58  2.7751167493305635e+00 -9.6251009716841907e-01  1.2847140239740573e+00
+   59  2.7332664162886871e+00  1.6874002693437411e+00 -1.8630367163899653e+00
+   60 -1.9920742303178280e+00 -5.1976595137487056e+00 -2.5715712256855774e+00
+   61 -3.9186643718089481e-01  1.8941052818415300e+00 -7.2852310082079819e-01
+   62 -2.6017296078018801e+00  2.4477411514482652e+00 -2.1740532348808843e+00
+   63  1.6751774499211569e+00 -2.9786131173037518e+00 -4.4746311293689334e-02
+   64  2.2350712682686358e+00  2.4856397598318183e+00  6.0442073184431804e+00
+run_vdwl: -258.49958453100083
+run_coul: 0
+run_stress: ! |2-
+   6.4308197607512207e+00  8.9800199663285003e+00  1.6644361563236384e+01 -5.0072799407290933e+00  3.4219821356106827e+01  1.0492826840671752e+00
+run_forces: ! |2
+    1 -7.5937225945190234e-01  3.0453338820624740e+00  1.7588997466292104e+00
+    2 -2.5842358982706335e+00 -6.4961465564215615e-01 -2.0836136472867430e+00
+    3  6.1692053064567098e-01 -1.7923404289169326e-01 -2.5776384969855004e-01
+    4 -1.8642434203060627e+00  2.9689111967459501e+00  5.7860631771685955e-01
+    5 -1.3707214670175798e+00 -7.6882987771621170e-01 -3.8201327706454102e-01
+    6  3.5155029154762341e-01  3.3059375310476078e+00  1.0816937673465203e+00
+    7 -1.0031791908375163e+00 -8.6956258364196903e-01  2.2425245701241256e+00
+    8  5.9471905500534855e-02  6.7712638453480301e-01 -9.9455163342971664e-01
+    9  1.3241459971867986e-01 -2.4451165560724939e+00 -2.5688790544877764e+00
+   10  2.7987758741645180e-01 -1.3869445825984155e+00 -2.2737317048352144e+00
+   11  3.0093494239500198e+00 -3.2460329093479796e+00  2.0365059310692528e+00
+   12 -5.2407962483722299e+00 -2.9126882828373892e+00 -2.3181917992781393e+00
+   13 -2.9313171711099439e-01  5.1682815536825712e+00 -9.3178794843042834e-01
+   14 -1.4205609886004000e+00  3.6728355746142340e+00  4.7973399106847658e-01
+   15 -3.6438539671591443e+00  4.1738585973197599e+00 -2.2995365264114183e+00
+   16  1.7665430073015751e+00 -1.6683625231758512e-01  5.5980371264386779e+00
+   17  1.1609247141613361e+00  2.5445069687633297e+00 -4.1101937186844282e+00
+   18  2.3250564926703476e-01  5.0119516573698042e-01  3.2999642189792335e+00
+   19 -1.2455755116140015e+00 -2.0483645027853901e+00 -3.9915953456528291e-01
+   20 -6.1533476323819176e+00 -1.0336110065476101e-01  7.2207408198501200e-01
+   21  2.2580812272923181e+00  8.8665782626079470e-01 -2.7867801726837302e+00
+   22 -5.0715437260287493e+00  3.2720805913066657e+00 -3.9870515109961557e+00
+   23  1.9067384153660658e+00  3.3666024181351721e+00  1.2360966484469924e+00
+   24  1.6076366668181858e+00  1.3129049141466480e+00  2.6481286770383776e+00
+   25  8.2849608748715486e-01 -1.0946573630945597e-01  1.4137379100581711e+00
+   26 -4.9245162995242892e-01 -6.5633188499443507e-01  1.2397189435323859e+00
+   27  2.0002065581614117e+00 -5.2635572637485362e-01  1.5812234025897816e+00
+   28  5.6469856770107940e-01  1.3419655834852620e+00  4.1390158667617349e+00
+   29 -3.0414702840545953e+00  9.2717141811819714e-01 -1.1446412927789829e+00
+   30 -1.5780449202280824e+00  1.1972994610432426e+00 -9.0937832538133057e-01
+   31 -1.6135468935693842e+00  1.4097747955859368e+00 -4.1013831752390606e-01
+   32 -3.3839861094059622e+00  5.3294244523702139e-01 -5.1150418397474551e-01
+   33  4.3179202422054441e+00 -2.0530236039826488e+00  2.7909047875298785e+00
+   34  3.6044667063424402e-01 -9.6636899664461151e-02  3.8719887922396817e-01
+   35  1.2574806392912077e+00 -1.5474474915574327e-01  2.5284696614857665e+00
+   36 -9.2507293281699435e-03 -2.5688826605251358e+00 -8.8917022692736422e-02
+   37 -1.7874734690551621e+00  1.7150271178623466e+00 -2.3271057624984728e+00
+   38 -2.4661620491429920e-01 -4.4857619628292789e-01  3.0090275228821284e-01
+   39  2.1876463817081504e+00 -1.2953665685720137e+00 -4.1070251548877099e+00
+   40  3.8954619284501790e+00 -1.9483835119233022e+00 -1.2227531289487130e+00
+   41  9.7565025130134986e-02 -5.2646053136013171e-01 -9.8280124384883216e-01
+   42 -4.1063684977228627e+00 -1.9251952046358491e+00 -2.9901720087297745e-01
+   43  6.7939692719606781e-01 -5.1874569217896553e+00 -3.9562141806867668e-01
+   44  3.2038887558519624e+00 -3.7312345150848865e+00  1.1791192854681274e+00
+   45 -1.1219574440264011e+00  1.4017870123386453e+00 -2.7737869798443744e+00
+   46 -2.3209294883240861e-01 -1.6182897275509909e+00  4.6296975808382292e+00
+   47  2.4602350215523230e+00  1.3452019889417968e+00  8.6358301807951909e-01
+   48  1.3093376440442195e+00  5.7371015819454974e-01 -3.5578408565068542e-01
+   49 -1.6696474765984912e+00 -4.7474039477977534e+00 -2.4607694981777111e+00
+   50  1.2989593641085762e-02  2.9812087098985032e+00 -1.4123464138675532e+00
+   51 -1.2853374902017072e+00  2.9768731587433566e+00  3.7403754374802136e+00
+   52  4.2768627438095717e+00 -4.4665207635060362e+00  4.4526942886425624e+00
+   53  2.8263576162367916e+00  2.8711821643474851e+00  1.6198032393140251e+00
+   54  2.4557393026757937e+00 -2.8576870371250473e+00 -1.7106753514157176e+00
+   55  8.7016650440838128e-01 -2.9993943594233770e+00 -3.5406400607258579e+00
+   56 -4.2720269240408387e-01  2.6284724466018883e-01 -2.0965871396618163e+00
+   57  1.8938909101993043e-01 -1.1500277319075998e+00 -3.0413618657653302e-01
+   58  2.7772395410738175e+00 -9.4351277181421511e-01  1.2588196612829905e+00
+   59  2.7573808001442890e+00  1.7074123239008328e+00 -1.8649940082533623e+00
+   60 -2.0326029067319991e+00 -5.2224816026297418e+00 -2.6126651787431188e+00
+   61 -3.8509493124293104e-01  1.9073567822530930e+00 -7.1625738092917579e-01
+   62 -2.6004175741383238e+00  2.4418292628003044e+00 -2.1615478367946075e+00
+   63  1.6782508782162848e+00 -3.0101538006328594e+00 -7.1986308168993651e-02
+   64  2.2749536899334060e+00  2.5303495677814132e+00  6.0668040667204073e+00
+...
diff --git a/unittest/force-styles/tests/mol-pair-coul_slater_long.yaml b/unittest/force-styles/tests/mol-pair-coul_slater_long.yaml
index becaaddacb..52efdaa52d 100644
--- a/unittest/force-styles/tests/mol-pair-coul_slater_long.yaml
+++ b/unittest/force-styles/tests/mol-pair-coul_slater_long.yaml
@@ -1,6 +1,6 @@
 ---
-lammps_version: 17 Feb 2022
-date_generated: Fri Mar 18 22:17:29 2022
+lammps_version: 23 Jun 2022
+date_generated: Thu Jul  7 09:00:39 2022
 epsilon: 2e-13
 skip_tests:
 prerequisites: ! |
@@ -21,71 +21,71 @@ extract: ! |
   cut_coul 0
 natoms: 29
 init_vdwl: 0
-init_coul: 531.6959978948178
+init_coul: 254.5592842067175
 init_stress: ! |2-
-   2.5972486027624825e+02  1.6017698793801901e+02  2.4137753464923463e+02  1.1278136362375312e+02  7.5965726957592040e+01  4.3719482306226765e-01
+   1.1708139667614163e-01 -2.1626456898855512e+01 -1.2229635861086752e+01  1.1761975361937868e+01 -9.4490379133316071e-01  7.0303045871968024e+00
 init_forces: ! |2
-    1 -1.4505459775863072e+01 -1.3044138972612792e+01  2.1164399702389403e+01
-    2  1.7223357637414050e+01  1.1196601670721074e+01 -2.1963333681835362e+01
-    3  1.9690750992824077e-01  8.4027350365370912e-01  3.2676598741321244e-01
-    4 -4.1823205578686712e-01 -4.8034612329946474e-01 -6.2465767147804363e-01
-    5 -8.1067146840659532e-01 -7.0419963313519196e-01  1.0662874339675656e-01
-    6 -1.5369327833513578e+01  1.5524058497083976e+01  2.0805079862911210e+01
-    7  2.5724194203247794e+00 -6.5060840392199450e+00 -2.5818209750744270e+01
-    8  5.4989014023099685e+00 -1.3233615378795388e+01 -1.9749310147567481e+01
-    9  6.7810825735410649e+00  2.9735670644350400e+00  2.9299979326310858e+01
-   10  1.5690354534134865e+00 -3.8317369465722124e+00 -4.6200627486826261e-02
-   11 -3.3949224827476737e-01  6.3351363413160511e-01  4.2844876717729496e-01
-   12  3.5635170964422114e-01 -1.4508136789019890e+00  3.4423369459456787e+00
-   13  3.1047049296276570e+00 -1.3699917700047379e+00 -1.9294387191193760e-01
-   14 -1.8808718796901713e+00  4.5935276633752392e-01 -3.4319894559470154e+00
-   15  5.0373263520628020e-01  3.2293151135885791e+00  1.4501551118828643e-01
-   16  9.3153064058538995e+00 -7.2444035690649828e+00 -2.8031955215690214e+01
-   17 -1.2854116750038774e+01  1.2711581809604329e+01  2.2641980278104416e+01
-   18 -4.3096955761758995e+00 -1.9223912022853014e+01  7.9000623254438963e+01
-   19 -3.6557250488749723e+01 -2.4008797712603496e+01 -4.6906379653253587e+01
-   20  4.0430244682537058e+01  4.3367408075615103e+01 -3.0568173216213165e+01
-   21 -2.2233648415360879e+01 -2.5175462230434459e+01  7.7556799814543822e+01
-   22 -3.2896290497889161e+01 -7.4008238647922955e+00 -5.6779521143055455e+01
-   23  5.4738206511131686e+01  3.2962897136082404e+01 -2.0302232663948129e+01
-   24  1.6836507718223405e+01 -6.8988891903011790e+01  4.0189734824880482e+01
-   25 -4.9788824311046007e+01  1.0298368296372720e+01 -4.2822722831087503e+01
-   26  3.2480024672491567e+01  5.8484108401966807e+01  2.1410689331760921e+00
-   27  1.0631780570743651e+01 -7.5812252190257794e+01  2.6785237094554518e+01
-   28 -5.1807337133621779e+01  2.4969152478376095e+01 -3.5819948334944876e+01
-   29  4.1532654602026298e+01  5.0825271587590592e+01  9.0234792187328665e+00
+    1  1.0123638170253564e+00 -1.2509301778440569e+00  2.1945424754579066e+00
+    2  1.2954585046924219e+00 -1.6338868361150480e+00 -3.0887376983421388e+00
+    3  2.5466335168114949e-02 -1.3139527326988556e-02  3.5569205669245270e-02
+    4 -1.2001437337468715e-01 -5.6411226640158199e-02 -2.7885583471425601e-01
+    5 -5.1442751513887475e-01  3.3619784208637826e-01 -5.8612509036156044e-02
+    6  2.0722678070474260e-01 -2.0158659756241346e+00  1.2290591371946380e-01
+    7  1.2567425258842402e-01  1.6224546679223717e-01  1.2761543372504829e+00
+    8 -1.0885715086892622e+00  2.7631213623524755e+00 -1.1767360015218493e-01
+    9  1.4894653246735305e+00 -3.7617960133817596e+00  3.1053923325311867e+00
+   10 -7.1355785210118067e-02  3.1060808027765684e-02 -1.9612261171752393e-01
+   11 -6.7234571380736030e-01  7.7973238090172237e-01 -4.0949779264856478e-01
+   12  2.3138844234286426e+00 -4.3745924070374187e-01  1.4448078371866295e+00
+   13  1.4660368459556450e-01 -1.5128217536790745e-01 -1.4303100577311031e-01
+   14 -8.0170764963523022e-01  2.5412749138723761e-01 -3.4566627523769955e-01
+   15  3.3139749151276299e-01  7.0760705097098225e-02 -7.8854914191841807e-01
+   16  1.6259331841614227e-01 -9.7505266997933193e-01 -1.0324909305416741e+00
+   17 -2.8980837212603436e+00  5.6015117342873459e+00 -3.2180999995572175e+00
+   18 -1.9982604761757349e-02  2.3277591984645070e+00  5.7125759352746586e-01
+   19 -9.3976757841446423e-01 -2.1798221084742129e+00  4.9912918497113157e-01
+   20  5.2304880078763782e-01 -1.3238749831696685e-02  4.5568360647359851e-01
+   21 -4.8423714829029962e-01  1.3829719239752545e+00 -2.9605371449607998e-01
+   22 -5.7875171427304151e-01 -7.4682110680773817e-01  1.1868729457280150e-01
+   23  6.7125646044497977e-01 -2.4953977631185778e-01  6.5241242746355066e-01
+   24  2.4761790219940844e-01  1.3513641756258139e+00 -5.1709840925954764e-01
+   25 -1.1849995949271324e+00 -1.0312698215846192e+00 -5.7625135838985497e-01
+   26  4.6508977239668642e-01 -5.2650955871347083e-01  6.0143069461846999e-01
+   27 -7.4501807019330935e-01  1.4539167477546999e+00 -6.3344271356143678e-01
+   28 -2.3479546881246116e-01 -1.0024440823389091e+00  1.1214449785783520e-01
+   29  1.3369115781539251e+00 -4.6930078970690337e-01  5.1006619404609643e-01
 run_vdwl: 0
-run_coul: 531.6316938435255
+run_coul: 254.55391272325517
 run_stress: ! |2-
-   2.6009236274938775e+02  1.6019149935649892e+02  2.4091972549720990e+02  1.1287278205291214e+02  7.6421249875941328e+01  1.0477154416076908e+00
+   1.0507664945308685e-01 -2.1621055562665443e+01 -1.2221091005466242e+01  1.1755690602099170e+01 -9.5887133975896877e-01  7.0291169104254232e+00
 run_forces: ! |2
-    1 -1.4505263340624662e+01 -1.3069826884538898e+01  2.1092273240232334e+01
-    2  1.7216929699569594e+01  1.1222653354918181e+01 -2.1892026658892760e+01
-    3  1.9689659133582163e-01  8.4089709839905469e-01  3.2691229507601138e-01
-    4 -4.1534123724459276e-01 -4.7942040205303499e-01 -6.2372283685814967e-01
-    5 -8.0983629907930077e-01 -7.0364522896705939e-01  1.0613059784230094e-01
-    6 -1.5349604856791370e+01  1.5522546924799940e+01  2.0756281297832874e+01
-    7  2.5720205395173057e+00 -6.5136694191526150e+00 -2.5754699982146299e+01
-    8  5.4783011177980798e+00 -1.3193790377127209e+01 -1.9732158944509290e+01
-    9  6.7823956616450456e+00  2.9328643102381267e+00  2.9267293877977238e+01
-   10  1.5701490834675309e+00 -3.8293776378819531e+00 -4.4681502716367193e-02
-   11 -3.4080685500280089e-01  6.3394138301719993e-01  4.2521300521638838e-01
-   12  3.6335920149821210e-01 -1.4542882227067147e+00  3.4250588402471536e+00
-   13  3.1021896119350991e+00 -1.3629775427304498e+00 -1.9236601260012270e-01
-   14 -1.8791529231416744e+00  4.5416896729772194e-01 -3.4214086468019449e+00
-   15  4.9947705909602841e-01  3.2294681516348525e+00  1.5065456540899641e-01
-   16  9.3152895494380115e+00 -7.2584606007912740e+00 -2.8037105174417249e+01
-   17 -1.2853807315028448e+01  1.2730457404190695e+01  2.2646649542757135e+01
-   18 -3.8896274444740597e+00 -1.8742984005196448e+01  7.8575622292818892e+01
-   19 -3.6657681210304887e+01 -2.4147367233205287e+01 -4.6880920994817359e+01
-   20  4.0110836482195197e+01  4.3026017430590635e+01 -3.0166674909873894e+01
-   21 -2.2266717418800365e+01 -2.5007792207614962e+01  7.7439548112199361e+01
-   22 -3.2986424583514150e+01 -7.4943575521934962e+00 -5.6719145374965272e+01
-   23  5.4861504180702049e+01  3.2888028352139067e+01 -2.0245497429195499e+01
-   24  1.6994333620812245e+01 -6.9071608975734406e+01  4.0281992417369800e+01
-   25 -5.0022578934028282e+01  1.0236526699813956e+01 -4.3029340089682030e+01
-   26  3.2556772729625038e+01  5.8630191330865074e+01  2.2575959885742858e+00
-   27  1.0740675179027832e+01 -7.5810812405906574e+01  2.6655288241835144e+01
-   28 -5.1859214575097845e+01  2.4967730720910470e+01 -3.5779461664571272e+01
-   29  4.1474926685469335e+01  5.0824886566985427e+01  9.1126959066595141e+00
+    1  1.0112873572811580e+00 -1.2530276549384933e+00  2.1914171798830568e+00
+    2  1.2910756618724584e+00 -1.6320587244040410e+00 -3.0856670933130275e+00
+    3  2.5417581409305202e-02 -1.3319027044403169e-02  3.5458204071365587e-02
+    4 -1.1903572341057568e-01 -5.6295979790903755e-02 -2.7900490837896463e-01
+    5 -5.1384628561765300e-01  3.3652540945011555e-01 -5.8862099577371089e-02
+    6  2.0557510022365150e-01 -2.0123835505068053e+00  1.2614266602622692e-01
+    7  1.2735793319460478e-01  1.5853467760283904e-01  1.2701822406914272e+00
+    8 -1.0845932144371644e+00  2.7620784785572590e+00 -1.1812899137169540e-01
+    9  1.4885988571079642e+00 -3.7664983214134815e+00  3.1092838634072559e+00
+   10 -7.1576372143287825e-02  3.0962708284801813e-02 -1.9638669006856277e-01
+   11 -6.7233593581762241e-01  7.8020791757188945e-01 -4.0861239677017752e-01
+   12  2.3142942085570231e+00 -4.3672211901097824e-01  1.4445175637184764e+00
+   13  1.4585640919457624e-01 -1.5151491185814414e-01 -1.4301903387021225e-01
+   14 -8.0211823601348731e-01  2.5440367515238671e-01 -3.4418722798600709e-01
+   15  3.3243392757719303e-01  6.9345316666666712e-02 -7.8979319762464939e-01
+   16  1.6228153766652689e-01 -9.7649864654379481e-01 -1.0325498581820041e+00
+   17 -2.8985741606086464e+00  5.6081200657787251e+00 -3.2210360520327095e+00
+   18 -2.4779821561789064e-02  2.3252812699411116e+00  5.7548476355559375e-01
+   19 -9.3073282073812136e-01 -2.1761428427184564e+00  5.0466585562048227e-01
+   20  5.1962859124254412e-01 -1.3823144408831962e-02  4.4616022167444958e-01
+   21 -4.8279265992792930e-01  1.3791791897341579e+00 -2.9615403241565003e-01
+   22 -5.7432971719583037e-01 -7.4279690320955927e-01  1.2017370691438947e-01
+   23  6.6583750391632179e-01 -2.5085912738116967e-01  6.5043205083268785e-01
+   24  2.4365888962669718e-01  1.3544800679852296e+00 -5.2159886075913586e-01
+   25 -1.1794061912936245e+00 -1.0316188106119735e+00 -5.6975657126355306e-01
+   26  4.6484502577860515e-01 -5.2739147192840230e-01  6.0227258976224862e-01
+   27 -7.4688344945317542e-01  1.4562803511493001e+00 -6.3143583726328756e-01
+   28 -2.3319583763082030e-01 -1.0037908986220945e+00  1.1105104982326264e-01
+   29  1.3360518412010989e+00 -4.7065699348295009e-01  5.0895089489608403e-01
 ...