Merge remote-tracking branch 'github/develop' into small_fixes

doc/src/Build_basics.rst

@@ -203,7 +203,7 @@ LAMMPS.
       check if the detected or selected compiler is compatible with the
       C++ support requirements of LAMMPS and stop with an error, if this
       is not the case.  A C++11 compatible compiler is currently
-      required, but a transition to require C++17 is in progess and
+      required, but a transition to require C++17 is in progress and
       planned to be completed in Summer 2025. Currently, setting
       ``-DLAMMPS_CXX11=yes`` is required when configuring with CMake while
       using a C++11 compatible compiler that does not support C++17,
@@ -329,7 +329,7 @@ LAMMPS.
          either as a binary package or through compiling from source.
 
       While a C++11 compatible compiler is currently sufficient to compile
-      LAMMPS, a transition to require C++17 is in progess and planned to
+      LAMMPS, a transition to require C++17 is in progress and planned to
       be completed in Summer 2025. Currently, setting ``-DLAMMPS_CXX11``
       in the ``LMP_INC =`` line in the machine makefile is required when
       using a C++11 compatible compiler that does not support C++17.

doc/src/Commands_pair.rst

@@ -115,7 +115,9 @@ OPT.
    * :doc:`gw/zbl <pair_gw>`
    * :doc:`harmonic/cut (o) <pair_harmonic_cut>`
    * :doc:`hbond/dreiding/lj (o) <pair_hbond_dreiding>`
+   * :doc:`hbond/dreiding/lj/angleoffset (o) <pair_hbond_dreiding>`
    * :doc:`hbond/dreiding/morse (o) <pair_hbond_dreiding>`
+   * :doc:`hbond/dreiding/morse/angleoffset (o) <pair_hbond_dreiding>`
    * :doc:`hdnnp <pair_hdnnp>`
    * :doc:`hippo (g) <pair_amoeba>`
    * :doc:`ilp/graphene/hbn (t) <pair_ilp_graphene_hbn>`

doc/src/Developer_code_design.rst

@@ -203,6 +203,7 @@ processed in the expected order before types are removed from dynamic
 dispatch.
 
 .. admonition:: Important Notes
+   :class: note
 
    In order to be able to detect incompatibilities at compile time and
    to avoid unexpected behavior, it is crucial that all member functions

doc/src/Fortran.rst

@@ -323,6 +323,12 @@ of the contents of the :f:mod:`LIBLAMMPS` Fortran interface to LAMMPS.
    :ftype set_internal_variable: subroutine
    :f eval: :f:func:`eval`
    :ftype eval: function
+   :f clearstep_compute: :f:subr:`clearstep_compute`
+   :ftype clearstep_compute: subroutine
+   :f addstep_compute: :f:subr:`addstep_compute`
+   :ftype addstep_compute: subroutine
+   :f addstep_compute_all: :f:subr:`addstep_compute_all`
+   :ftype addstep_compute_all: subroutine
    :f gather_atoms: :f:subr:`gather_atoms`
    :ftype gather_atoms: subroutine
    :f gather_atoms_concat: :f:subr:`gather_atoms_concat`
@@ -956,6 +962,7 @@ Procedures Bound to the :f:type:`lammps` Derived Type
       :f:func:`extract_atom` between runs.
 
    .. admonition:: Array index order
+      :class: tip
 
       Two-dimensional arrays returned from :f:func:`extract_atom` will be
       **transposed** from equivalent arrays in C, and they will be indexed
@@ -1068,6 +1075,7 @@ Procedures Bound to the :f:type:`lammps` Derived Type
    you based on data from the :cpp:class:`Compute` class.
 
    .. admonition:: Array index order
+      :class: tip
 
       Two-dimensional arrays returned from :f:func:`extract_compute` will be
       **transposed** from equivalent arrays in C, and they will be indexed
@@ -1326,6 +1334,7 @@ Procedures Bound to the :f:type:`lammps` Derived Type
    :rtype data: polymorphic
 
    .. admonition:: Array index order
+      :class: tip
 
       Two-dimensional global, per-atom, or local array data from
       :f:func:`extract_fix` will be **transposed** from equivalent arrays in
@@ -1450,11 +1459,62 @@ Procedures Bound to the :f:type:`lammps` Derived Type
    an internal-style variable, an error is generated.
 
    :p character(len=*) name: name of the variable
-   :p read(c_double) val:  new value to assign to the variable
+   :p real(c_double) val:  new value to assign to the variable
    :to: :cpp:func:`lammps_set_internal_variable`
 
 --------
 
+.. f:function:: eval(expr)
+
+   This function is a wrapper around :cpp:func:`lammps_eval` that takes a
+   LAMMPS equal style variable string, evaluates it and returns the resulting
+   scalar value as a floating-point number.
+
+   .. versionadded:: TBD
+
+   :p character(len=\*) expr: string to be evaluated
+   :to: :cpp:func:`lammps_eval`
+   :r value [real(c_double)]: result of the evaluated string
+
+--------
+
+.. f:subroutine:: clearstep_compute()
+
+   Clear whether a compute has been invoked
+
+   .. versionadded:: TBD
+
+   :to: :cpp:func:`lammps_clearstep_compute`
+
+--------
+
+.. f:subroutine:: addstep_compute(nextstep)
+
+   Add timestep to list of future compute invocations
+   if the compute has been invoked on the current timestep
+
+   .. versionadded:: TBD
+
+   overloaded for 32-bit and 64-bit integer arguments
+
+   :p integer(kind=8 or kind=4) nextstep: next timestep
+   :to: :cpp:func:`lammps_addstep_compute`
+
+--------
+
+.. f:subroutine:: addstep_compute_all(nextstep)
+
+   Add timestep to list of future compute invocations
+
+   .. versionadded:: TBD
+
+   overloaded for 32-bit and 64-bit integer arguments
+
+   :p integer(kind=8 or kind=4) nextstep: next timestep
+   :to: :cpp:func:`lammps_addstep_compute_all`
+
+--------
+
 .. f:subroutine:: gather_atoms(name, count, data)
 
    This function calls :cpp:func:`lammps_gather_atoms` to gather the named

doc/src/Howto_barostat.rst

@@ -10,20 +10,21 @@ and/or pressure (P) is specified by the user, and the thermostat or
 barostat attempts to equilibrate the system to the requested T and/or
 P.
 
-Barostatting in LAMMPS is performed by :doc:`fixes <fix>`.  Two
+Barostatting in LAMMPS is performed by :doc:`fixes <fix>`.  Three
 barostatting methods are currently available: Nose-Hoover (npt and
-nph) and Berendsen:
+nph), Berendsen, and various linear controllers in deform/pressure:
 
 * :doc:`fix npt <fix_nh>`
 * :doc:`fix npt/sphere <fix_npt_sphere>`
 * :doc:`fix npt/asphere <fix_npt_asphere>`
 * :doc:`fix nph <fix_nh>`
 * :doc:`fix press/berendsen <fix_press_berendsen>`
+* :doc:`fix deform/pressure <fix_deform_pressure>`
 
 The :doc:`fix npt <fix_nh>` commands include a Nose-Hoover thermostat
 and barostat.  :doc:`Fix nph <fix_nh>` is just a Nose/Hoover barostat;
-it does no thermostatting.  Both :doc:`fix nph <fix_nh>` and :doc:`fix press/berendsen <fix_press_berendsen>` can be used in conjunction
+it does no thermostatting.  The fixes :doc:`nph <fix_nh>`, :doc:`press/berendsen <fix_press_berendsen>`, and :doc:`deform/pressure <fix_deform_pressure>`
-with any of the thermostatting fixes.
+can be used in conjunction with any of the thermostatting fixes.
 
 As with the :doc:`thermostats <Howto_thermostat>`, :doc:`fix npt <fix_nh>`
 and :doc:`fix nph <fix_nh>` only use translational motion of the
@@ -44,9 +45,9 @@ a temperature or pressure compute to a barostatting fix.
 .. note::
 
    As with the thermostats, the Nose/Hoover methods (:doc:`fix npt <fix_nh>` and :doc:`fix nph <fix_nh>`) perform time integration.
-   :doc:`Fix press/berendsen <fix_press_berendsen>` does NOT, so it should
+   :doc:`Fix press/berendsen <fix_press_berendsen>` and :doc:`fix deform/pressure <fix_deform_pressure>`
-   be used with one of the constant NVE fixes or with one of the NVT
+   do NOT, so they should be used with one of the constant NVE fixes or with
-   fixes.
+   one of the NVT fixes.
 
 Thermodynamic output, which can be setup via the
 :doc:`thermo_style <thermo_style>` command, often includes pressure

doc/src/Install_git.rst

@@ -52,6 +52,7 @@ your machine and "release" is one of the 3 branches listed above.
 between them at any time using "git checkout <branch name>".)
 
 .. admonition:: Saving time and disk space when using ``git clone``
+   :class: note
 
    The complete git history of the LAMMPS project is quite large because
    it contains the entire commit history of the project since fall 2006,

doc/src/Library_objects.rst

@@ -13,6 +13,9 @@ fixes, or variables in LAMMPS using the following functions:
 - :cpp:func:`lammps_set_internal_variable`
 - :cpp:func:`lammps_variable_info`
 - :cpp:func:`lammps_eval`
+- :cpp:func:`lammps_clearstep_compute`
+- :cpp:func:`lammps_addstep_compute_all`
+- :cpp:func:`lammps_addstep_compute`
 
 -----------------------
 
@@ -61,6 +64,21 @@ fixes, or variables in LAMMPS using the following functions:
 
 -----------------------
 
+.. doxygenfunction:: lammps_clearstep_compute
+   :project: progguide
+
+-----------------------
+
+.. doxygenfunction:: lammps_addstep_compute_all
+   :project: progguide
+
+-----------------------
+
+.. doxygenfunction:: lammps_addstep_compute
+   :project: progguide
+
+-----------------------
+
 .. doxygenenum:: _LMP_DATATYPE_CONST
 
 .. doxygenenum:: _LMP_STYLE_CONST

doc/src/fix_python_invoke.rst

@@ -66,6 +66,15 @@ gives access to the LAMMPS state from Python.
    from these callbacks, trying to execute input script commands will in the best
    case not work or in the worst case result in undefined behavior.
 
+Restart, fix_modify, output, run start/stop, minimize info
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+No information about this fix is written to :doc:`binary restart files <restart>`.  None of the :doc:`fix_modify <fix_modify>` options
+are relevant to this fix.  No global or per-atom quantities are stored
+by this fix for access by various :doc:`output commands <Howto_output>`.
+No parameter of this fix can be used with the *start/stop* keywords of
+the :doc:`run <run>` command.  This fix is not invoked during :doc:`energy minimization <minimize>`.
+
 Restrictions
 """"""""""""
 

doc/src/fix_reaxff_species.rst

@@ -200,8 +200,8 @@ The 2 values in the global vector are as follows:
 The per-atom vector stores the molecule ID for each atom as identified
 by the fix.  If an atom is not in a molecule, its ID will be 0.
 For atoms in the same molecule, the molecule ID for all of them
-will be the same and will be equal to the smallest atom ID of
+will be the same, and molecule IDs will range from 1 to the number
-any atom in the molecule.
+of molecules.
 
 No parameter of this fix can be used with the *start/stop* keywords of
 the :doc:`run <run>` command.

doc/src/fix_wall_gran.rst

@@ -248,11 +248,11 @@ listed in the following table.
 |     8 | Radius :math:`r` of atom                           | distance units |
 +-------+----------------------------------------------------+----------------+
 
-If a granular submodel calculates additional contact information (e.g. the
+If a granular sub-model calculates additional contact information (e.g. the
-heat submodels calculate the amount of heat exchanged), these quantities
+heat sub-models calculate the amount of heat exchanged), these quantities
 are appended to the end of this array. First, any extra values from the
-normal submodel are appended followed by the damping, tangential, rolling,
+normal sub-model are appended followed by the damping, tangential, rolling,
-twisting, then heat models. See the descriptions of granular submodels in
+twisting, then heat models. See the descriptions of granular sub-models in
 the :doc:`pair granular <pair_granular>` page for information on any extra
 quantities.
 

doc/src/fix_wall_gran_region.rst

@@ -269,11 +269,11 @@ listed in the following table.
 |     8 | Radius :math:`r` of atom                           | distance units |
 +-------+----------------------------------------------------+----------------+
 
-If a granular submodel calculates additional contact information (e.g. the
+If a granular sub-model calculates additional contact information (e.g. the
-heat submodels calculate the amount of heat exchanged), these quantities
+heat sub-models calculate the amount of heat exchanged), these quantities
 are appended to the end of this array. First, any extra values from the
-normal submodel are appended followed by the damping, tangential, rolling,
+normal sub-model are appended followed by the damping, tangential, rolling,
-twisting, then heat models. See the descriptions of granular submodels in
+twisting, then heat models. See the descriptions of granular sub-models in
 the :doc:`pair granular <pair_granular>` page for information on any extra
 quantities.
 

doc/src/pair_granular.rst

@@ -258,7 +258,7 @@ in damping model.
 
 The definition of multiple *mdr* models in the *pair_style* is currently not
 supported. Similarly, the *mdr* model cannot be combined with a different normal
-model in the *pair_style*. Physically this means that only one homogenous
+model in the *pair_style*. Physically this means that only one homogeneous
 collection of particles governed by a single *mdr* model is allowed.
 
 The *mdr* model currently only supports *fix wall/gran/region*, not
@@ -303,6 +303,7 @@ radius in the *mdr* model, the keyword/arg pair *cutoff radius* must be specifie
    simulation involving 200 particles named *in.tableting.200*.
    The second is a triaxial compaction simulation involving 12
    particles named *in.triaxial.compaction.12*.
+
 ----------
 
 In addition, the normal force is augmented by a damping term of the
@@ -941,16 +942,16 @@ particle I. The next entry (8) is the magnitude of the rolling torque.
 The next entry (9) is the magnitude of the twisting torque acting
 about the vector connecting the two particle centers.
 The next 3 (10-12) are the components of the vector connecting
-the centers of the two particles (x_I - x_J). If a granular submodel
+the centers of the two particles (x_I - x_J). If a granular sub-model
-calculates additional contact information (e.g. the heat submodels
+calculates additional contact information (e.g. the heat sub-models
 calculate the amount of heat exchanged), these quantities are appended
-to the end of this list. First, any extra values from the normal submodel
+to the end of this list. First, any extra values from the normal sub-model
 are appended followed by the damping, tangential, rolling, twisting, then
-heat models. See the descriptions of specific granular submodels above
+heat models. See the descriptions of specific granular sub-models above
 for information on any extra quantities. If two or more models are
 defined by pair coefficients, the size of the array is set by the
 maximum number of extra quantities in a model but the order of quantities
-is determined by each model's specific set of submodels. Any unused
+is determined by each model's specific set of sub-models. Any unused
 quantities are zeroed.
 
 These extra quantities can be accessed by the :doc:`compute pair/local <compute_pair_local>` command, as *p1*, *p2*, ...,
@@ -1046,7 +1047,7 @@ a bulk elastic response. Journal of the Mechanics and Physics of Solids,
 **(Zunker et al, 2025)** Zunker, W., Dunatunga, S., Thakur, S.,
 Tang, P., & Kamrin, K. (2025). Experimentally validated DEM for large
 deformation powder compaction: mechanically-derived contact model and
-screening of non-physical contacts. engrXiv.
+screening of non-physical contacts.
 
 .. _Luding2008:
 

doc/src/pair_hbond_dreiding.rst

@@ -1,30 +1,46 @@
 .. index:: pair_style hbond/dreiding/lj
 .. index:: pair_style hbond/dreiding/lj/omp
+.. index:: pair_style hbond/dreiding/lj/angleoffset
+.. index:: pair_style hbond/dreiding/lj/angleoffset/omp
 .. index:: pair_style hbond/dreiding/morse
 .. index:: pair_style hbond/dreiding/morse/omp
+.. index:: pair_style hbond/dreiding/morse/angleoffset
+.. index:: pair_style hbond/dreiding/morse/angleoffset/omp
 
 pair_style hbond/dreiding/lj command
 ====================================
 
 Accelerator Variants: *hbond/dreiding/lj/omp*
 
+pair_style hbond/dreiding/lj/angleoffset command
+================================================
+
+Accelerator Variants: *hbond/dreiding/lj/angleoffset/omp*
+
 pair_style hbond/dreiding/morse command
 =======================================
 
 Accelerator Variants: *hbond/dreiding/morse/omp*
 
+pair_style hbond/dreiding/morse/angleoffset command
+===================================================
+
+Accelerator Variants: *hbond/dreiding/morse/angleoffset/omp*
+
+
 Syntax
 """"""
 
 .. code-block:: LAMMPS
 
-   pair_style style N inner_distance_cutoff outer_distance_cutoff angle_cutoff
+   pair_style style N inner_distance_cutoff outer_distance_cutoff angle_cutoff equilibrium_angle
 
-* style = *hbond/dreiding/lj* or *hbond/dreiding/morse*
+* style = *hbond/dreiding/lj* or *hbond/dreiding/morse* or *hbond/dreiding/lj/angleoffset* or *hbond/dreiding/morse/angleoffset*
 * N = power of cosine of angle theta (integer)
 * inner_distance_cutoff = global inner cutoff for Donor-Acceptor interactions (distance units)
 * outer_distance_cutoff = global cutoff for Donor-Acceptor interactions (distance units)
 * angle_cutoff = global angle cutoff for Acceptor-Hydrogen-Donor interactions (degrees)
+* (with style angleoffset) equilibrium_angle = global equilibrium angle for Acceptor-Hydrogen-Donor interactions (degrees)
 
 Examples
 """"""""
@@ -40,6 +56,9 @@ Examples
    labelmap atom 1 C 2 O 3 H
    pair_coeff C O hbond/dreiding/morse H i 3.88 1.7241379 2.9 2 9.0 11.0 90.0
 
+   pair_style hybrid/overlay lj/cut 10.0 hbond/dreiding/lj 4 9.0 11.0 90 170.0
+   pair_coeff 1 2 hbond/dreiding/lj 3 i 9.5 2.75 4 9.0 11.0 90.0
+
 Description
 """""""""""
 
@@ -74,42 +93,53 @@ hydrogen (H) and the donor atoms:
 .. image:: JPG/dreiding_hbond.jpg
    :align: center
 
-These 3-body interactions can be defined for pairs of acceptor and
+These 3-body interactions can be defined for pairs of acceptor and donor
-donor atoms, based on atom types.  For each donor/acceptor atom pair,
+atoms, based on atom types.  For each donor/acceptor atom pair, the
-the third atom in the interaction is a hydrogen permanently bonded to
+third atom in the interaction is a hydrogen permanently bonded to the
-the donor atom, e.g. in a bond list read in from a data file via the
+donor atom, e.g. in a bond list read in from a data file via the
 :doc:`read_data <read_data>` command.  The atom types of possible
 hydrogen atoms for each donor/acceptor type pair are specified by the
 :doc:`pair_coeff <pair_coeff>` command (see below).
 
 Style *hbond/dreiding/lj* is the original DREIDING potential of
-:ref:`(Mayo) <pair-Mayo>`.  It uses a LJ 12/10 functional for the Donor-Acceptor
+:ref:`(Mayo) <pair-Mayo>`.  It uses a LJ 12/10 functional for the
-interactions. To match the results in the original paper, use n = 4.
+Donor-Acceptor interactions. To match the results in the original paper,
+use n = 4.
 
 Style *hbond/dreiding/morse* is an improved version using a Morse
 potential for the Donor-Acceptor interactions. :ref:`(Liu) <Liu>` showed
 that the Morse form gives improved results for Dendrimer simulations,
 when n = 2.
 
+.. versionadded:: TBD
+
+The style variants *hbond/dreiding/lj/angleoffset* and
+*hbond/dreiding/lj/angleoffset* take the equilibrium angle of the AHD as
+input, allowing it to reach 180 degrees. This variant option was added
+to account for cases (especially in some coarse-grained models) in which
+the equilibrium state of the bonds may equal the minimum energy state.
+
 See the :doc:`Howto bioFF <Howto_bioFF>` page for more information
 on the DREIDING force field.
 
 .. note::
 
-   Because the Dreiding hydrogen bond potential is only one portion
+   Because the Dreiding hydrogen bond potential is only one portion of
-   of an overall force field which typically includes other pairwise
+   an overall force field which typically includes other pairwise
-   interactions, it is common to use it as a sub-style in a :doc:`pair_style hybrid/overlay <pair_hybrid>` command, where another pair style
+   interactions, it is common to use it as a sub-style in a
-   provides the repulsive core interaction between pairs of atoms, e.g. a
+   :doc:`pair_style hybrid/overlay <pair_hybrid>` command, where another
-   1/r\^12 Lennard-Jones repulsion.
+   pair style provides the repulsive core interaction between pairs of
+   atoms, e.g. a 1/r\^12 Lennard-Jones repulsion.
 
 .. note::
 
-   When using the hbond/dreiding pair styles with :doc:`pair_style hybrid/overlay <pair_hybrid>`, you should explicitly define pair
+   When using the hbond/dreiding pair styles with :doc:`pair_style
+   hybrid/overlay <pair_hybrid>`, you should explicitly define pair
    interactions between the donor atom and acceptor atoms, (as well as
    between these atoms and ALL other atoms in your system).  Whenever
-   :doc:`pair_style hybrid/overlay <pair_hybrid>` is used, ordinary mixing
+   :doc:`pair_style hybrid/overlay <pair_hybrid>` is used, ordinary
-   rules are not applied to atoms like the donor and acceptor atoms
+   mixing rules are not applied to atoms like the donor and acceptor
-   because they are typically referenced in multiple pair styles.
+   atoms because they are typically referenced in multiple pair styles.
    Neglecting to do this can cause difficult-to-detect physics problems.
 
 .. note::
@@ -119,6 +149,13 @@ on the DREIDING force field.
    special_bonds command (e.g. "special_bonds lj 0.0 0.0 1.0") to turn
    these interactions on.
 
+.. note::
+
+   For the *angleoffset* variants, the referenced angle offset is the
+   supplementary angle of the equilibrium angle parameter. It means if
+   the equilibrium angle is 166.6 degrees, the calculated angle offset
+   is 13.4 degrees.
+
 ----------
 
 The following coefficients must be defined for pairs of eligible
@@ -169,7 +206,10 @@ follows:
 * distance cutoff :math:`r_{out}` (distance units)
 * angle cutoff (degrees)
 
-A single hydrogen atom type K can be specified, or a wild-card asterisk
+For both the *hbond/dreiding/lj/angleoffset* and *hbond/dreiding/morse/angleoffset* styles an additional parameter is added:
+* equilibrium angle (degrees)
+
+For all styles, a single hydrogen atom type K can be specified, or a wild-card asterisk
 can be used in place of or in conjunction with the K arguments to
 select multiple types as hydrogen atoms.  This takes the form
 "\*" or "\*n" or "n\*" or "m\*n".  See the :doc:`pair_coeff <pair_coeff>`
@@ -245,8 +285,7 @@ heading) the following commands could be included in an input script:
 Restrictions
 """"""""""""
 
-This pair style can only be used if LAMMPS was built with the
+The base pair styles can only be used if LAMMPS was built with the MOLECULE package.  The *angleoffset* variant also requires the EXTRA-MOLECULE package.  See the :doc:`Build package <Build_package>` doc page
-MOLECULE package.  See the :doc:`Build package <Build_package>` doc page
 for more info.
 
 Related commands

doc/src/pair_mliap.rst

@@ -145,6 +145,7 @@ per line.
 The detail of *nn* module implementation can be found at :ref:`(Yanxon) <Yanxon2020>`.
 
 .. admonition:: Notes on mliappy models
+   :class: note
 
    When the *model* keyword is *mliappy*, if the filename ends in '.pt',
    or '.pth', it will be loaded using pytorch; otherwise, it will be

doc/src/pair_reaxff.rst

@@ -158,11 +158,36 @@ drops to zero.
 Optional keywords *safezone*, *mincap*, and *minhbonds* are used
 for allocating reaxff arrays.  Increasing these values can avoid memory
 problems, such as segmentation faults and bondchk failed errors, that
-could occur under certain conditions. These keywords are not used by
+could occur under certain conditions. These keywords are **not** used by
 the Kokkos version, which instead uses a more robust memory allocation
 scheme that checks if the sizes of the arrays have been exceeded and
 automatically allocates more memory.
 
+.. admonition:: Memory management problems with ReaxFF
+   :class: tip
+
+   The LAMMPS implementation of ReaxFF is adapted from a standalone MD
+   program written in C called `PuReMD
+   <https://github.com/msu-sparta/PuReMD>`_.  It inherits from this code
+   a heuristic memory management that is different from what the rest of
+   LAMMPS uses.  It assumes that a system is dense and already well
+   equilibrated, so that there are no large changes in how many and what
+   types of neighbors atoms have.  However, not all systems are like
+   that, and thus there can be errors or segmentation faults if the
+   system changes too much.  If you run into problems, here are three
+   options to avoid them:
+
+   - Use the KOKKOS version of ReaxFF (KOKKOS is not only for GPUs,
+     but can also be compiled for serial or OpenMP execution) which
+     uses a different memory management approach.
+   - Break down a run command during which memory related errors happen
+     into multiple smaller segments so that the memory management
+     heuristics are re-initialized for each segment before they become
+     invalid.
+   - Increase the values for *safezone*, *mincap*, and *minhbonds* as
+     needed.  This can lead to significant increase of memory consumption
+     through.
+
 The keyword *tabulate* controls the size of interpolation table for
 Lennard-Jones and Coulomb interactions. Tabulation may also be set in the
 control file (see below). If tabulation is set in both the input script and the

doc/src/pair_style.rst

@@ -207,7 +207,9 @@ accelerated styles exist.
 * :doc:`gw/zbl <pair_gw>` - Gao-Weber potential with a repulsive ZBL core
 * :doc:`harmonic/cut <pair_harmonic_cut>` - repulsive-only harmonic potential
 * :doc:`hbond/dreiding/lj <pair_hbond_dreiding>` - DREIDING hydrogen bonding LJ potential
+* :doc:`hbond/dreiding/lj/angleoffset <pair_hbond_dreiding>` - DREIDING hydrogen bonding LJ potential with offset for hbond angle
 * :doc:`hbond/dreiding/morse <pair_hbond_dreiding>` - DREIDING hydrogen bonding Morse potential
+* :doc:`hbond/dreiding/morse/angleoffset <pair_hbond_dreiding>` - DREIDING hydrogen bonding Morse potential with offset for hbond angle
 * :doc:`hdnnp <pair_hdnnp>` - High-dimensional neural network potential
 * :doc:`hippo <pair_amoeba>` -
 * :doc:`ilp/graphene/hbn <pair_ilp_graphene_hbn>` - registry-dependent interlayer potential (ILP)

doc/utils/sphinx-config/false_positives.txt

@@ -108,6 +108,7 @@ Andrienko
 Andzelm
 Ang
 anglegrad
+angleoffset
 angletangrad
 angmom
 angmomx
@@ -1583,6 +1584,7 @@ Impropers
 imulator
 includelink
 incompressible
+incompressibility
 incrementing
 indenter
 indenters
@@ -1762,6 +1764,7 @@ Kadiri
 Kai
 Kalia
 Kamberaj
+Kamrin
 Kantorovich
 Kapfer
 Kapil
@@ -2532,6 +2535,7 @@ Nevery
 Nevins
 newfile
 Newns
+newstep
 newtype
 nextsort
 Neyts
@@ -3731,6 +3735,7 @@ tgnpt
 tgnvt
 th
 Thakkar
+Thakur
 Thaokar
 thb
 thei
@@ -3769,6 +3774,7 @@ Tigran
 Tij
 Tildesley
 Timan
+timeflag
 timeI
 timespan
 timestamp
@@ -3831,6 +3837,7 @@ Tref
 Tretyakov
 tri
 triangleflag
+triaxial
 Tribello
 triclinic
 Triclinic

examples/COUPLE/plugin/liblammpsplugin.c

@@ -118,6 +118,9 @@ liblammpsplugin_t *liblammpsplugin_load(const char *lib)
   ADDSYM(set_internal_variable);
   ADDSYM(variable_info);
   ADDSYM(eval);
+  ADDSYM(clearstep_compute);
+  ADDSYM(addstep_compute);
+  ADDSYM(addstep_compute_all);
 
   ADDSYM(gather_atoms);
   ADDSYM(gather_atoms_concat);

examples/COUPLE/plugin/liblammpsplugin.h

@@ -164,6 +164,9 @@ struct _liblammpsplugin {
   int (*set_internal_variable)(void *, const char *, double);
   int (*variable_info)(void *, int, char *, int);
   double (*eval)(void *, const char *);
+  void (*clearstep_compute)(void *);
+  void (*addstep_compute)(void *, void *);
+  void (*addstep_compute_all)(void *, void *);
 
   void (*gather_atoms)(void *, const char *, int, int, void *);
   void (*gather_atoms_concat)(void *, const char *, int, int, void *);

examples/PACKAGES/stressprofile/in.flat

@@ -32,7 +32,7 @@ fix             1 all nvt temp 0.7 0.7 0.2
 #dump_modify    3 pad 3
 
 fix 2 all recenter NULL NULL 15 units lattice
-compute p1 all stress/cartesian z 0.5
+compute p1 all stress/cartesian z 0.5 NULL 0
 fix 3 all ave/time 100 1 100 c_p1[*] file flat.out mode vector
 
 thermo          50

examples/multi/in.granular

@@ -1,4 +1,4 @@
-# Big colloid particles and small LJ particles
+# Binary granular system
 
 units           lj
 atom_style      sphere

fortran/lammps.f90

@@ -127,6 +127,16 @@ MODULE LIBLAMMPS
     PROCEDURE :: set_string_variable    => lmp_set_string_variable
     PROCEDURE :: set_internal_variable  => lmp_set_internal_variable
     PROCEDURE :: eval                   => lmp_eval
+
+    PROCEDURE :: clearstep_compute      => lmp_clearstep_compute
+    PROCEDURE, PRIVATE :: lmp_addstep_compute_smallint
+    PROCEDURE, PRIVATE :: lmp_addstep_compute_bigint
+    GENERIC :: addstep_compute          => lmp_addstep_compute_smallint, lmp_addstep_compute_bigint
+    PROCEDURE, PRIVATE :: lmp_addstep_compute_all_smallint
+    PROCEDURE, PRIVATE :: lmp_addstep_compute_all_bigint
+    GENERIC :: addstep_compute_all      => lmp_addstep_compute_all_smallint, &
+                                           lmp_addstep_compute_all_bigint
+
     PROCEDURE, PRIVATE :: lmp_gather_atoms_int
     PROCEDURE, PRIVATE :: lmp_gather_atoms_double
     GENERIC   :: gather_atoms           => lmp_gather_atoms_int, &
@@ -626,6 +636,24 @@ MODULE LIBLAMMPS
       REAL(c_double) :: lammps_eval
     END FUNCTION lammps_eval
 
+    SUBROUTINE lammps_clearstep_compute(handle) BIND(C)
+      IMPORT :: c_ptr
+      IMPLICIT NONE
+      TYPE(c_ptr), VALUE :: handle
+    END SUBROUTINE lammps_clearstep_compute
+
+    SUBROUTINE lammps_addstep_compute(handle, step) BIND(C)
+      IMPORT :: c_ptr
+      IMPLICIT NONE
+      TYPE(c_ptr), VALUE :: handle, step
+    END SUBROUTINE lammps_addstep_compute
+
+    SUBROUTINE lammps_addstep_compute_all(handle, step) BIND(C)
+      IMPORT :: c_ptr
+      IMPLICIT NONE
+      TYPE(c_ptr), VALUE :: handle, step
+    END SUBROUTINE lammps_addstep_compute_all
+
     SUBROUTINE lammps_gather_atoms(handle, name, TYPE, count, DATA) BIND(C)
       IMPORT :: c_int, c_ptr
       IMPLICIT NONE
@@ -1846,6 +1874,80 @@ CONTAINS
     CALL lammps_free(Cexpr)
   END FUNCTION lmp_eval
 
+  ! equivalent subroutine to lammps_clearstep_compute
+  SUBROUTINE lmp_clearstep_compute(self)
+    CLASS(lammps), INTENT(IN) :: self
+    CALL lammps_clearstep_compute(self%handle)
+  END SUBROUTINE lmp_clearstep_compute
+
+  ! equivalent subroutine to lammps_addstep_compute
+  SUBROUTINE lmp_addstep_compute_bigint(self, nextstep)
+    CLASS(lammps), INTENT(IN) :: self
+    INTEGER(kind=8), INTENT(IN) :: nextstep
+    INTEGER(c_int), TARGET :: smallstep
+    INTEGER(c_int64_t), TARGET :: bigstep
+    TYPE(c_ptr) :: ptrstep
+    IF (SIZE_BIGINT == 4_c_int) THEN
+        smallstep = INT(nextstep,kind=c_int)
+        ptrstep = C_LOC(smallstep)
+    ELSE
+        bigstep = nextstep
+        ptrstep = C_LOC(bigstep)
+    END IF
+    CALL lammps_addstep_compute(self%handle, ptrstep)
+  END SUBROUTINE lmp_addstep_compute_bigint
+
+  ! equivalent subroutine to lammps_addstep_compute
+  SUBROUTINE lmp_addstep_compute_smallint(self, nextstep)
+    CLASS(lammps), INTENT(IN) :: self
+    INTEGER(kind=4), INTENT(IN) :: nextstep
+    INTEGER(c_int), TARGET :: smallstep
+    INTEGER(c_int64_t), TARGET :: bigstep
+    TYPE(c_ptr) :: ptrstep
+    IF (SIZE_BIGINT == 4_c_int) THEN
+        smallstep = nextstep
+        ptrstep = C_LOC(smallstep)
+    ELSE
+        bigstep = nextstep
+        ptrstep = C_LOC(bigstep)
+    END IF
+    CALL lammps_addstep_compute(self%handle, ptrstep)
+  END SUBROUTINE lmp_addstep_compute_smallint
+
+  ! equivalent subroutine to lammps_addstep_compute_all
+  SUBROUTINE lmp_addstep_compute_all_bigint(self, nextstep)
+    CLASS(lammps), INTENT(IN) :: self
+    INTEGER(kind=8), INTENT(IN) :: nextstep
+    INTEGER(c_int), TARGET :: smallstep
+    INTEGER(c_int64_t), TARGET :: bigstep
+    TYPE(c_ptr) :: ptrstep
+    IF (SIZE_BIGINT == 4_c_int) THEN
+        smallstep = INT(nextstep,kind=c_int)
+        ptrstep = C_LOC(smallstep)
+    ELSE
+        bigstep = nextstep
+        ptrstep = C_LOC(bigstep)
+    END IF
+    CALL lammps_addstep_compute_all(self%handle, ptrstep)
+  END SUBROUTINE lmp_addstep_compute_all_bigint
+
+  ! equivalent subroutine to lammps_addstep_compute_all
+  SUBROUTINE lmp_addstep_compute_all_smallint(self, nextstep)
+    CLASS(lammps), INTENT(IN) :: self
+    INTEGER(kind=4), INTENT(IN) :: nextstep
+    INTEGER(c_int), TARGET :: smallstep
+    INTEGER(c_int64_t), TARGET :: bigstep
+    TYPE(c_ptr) :: ptrstep
+    IF (SIZE_BIGINT == 4_c_int) THEN
+        smallstep = nextstep
+        ptrstep = C_LOC(smallstep)
+    ELSE
+        bigstep = nextstep
+        ptrstep = C_LOC(bigstep)
+    END IF
+    CALL lammps_addstep_compute_all(self%handle, ptrstep)
+  END SUBROUTINE lmp_addstep_compute_all_smallint
+
   ! equivalent function to lammps_gather_atoms (for integers)
   SUBROUTINE lmp_gather_atoms_int(self, name, count, data)
     CLASS(lammps), INTENT(IN) :: self

python/lammps/core.py

@@ -422,6 +422,10 @@ class lammps(object):
     self.lib.lammps_extract_variable_datatype.argtypes = [c_void_p, c_char_p]
     self.lib.lammps_extract_variable_datatype.restype = c_int
 
+    self.lib.lammps_clearstep_compute.argtype = [c_void_p]
+    self.lib.lammps_addstep_compute.argtype = [c_void_p, c_void_p]
+    self.lib.lammps_addstep_compute_all.argtype = [c_void_p, c_void_p]
+
     self.lib.lammps_eval.argtypes = [c_void_p, c_char_p]
     self.lib.lammps_eval.restype = c_double
 
@@ -1594,6 +1598,26 @@ class lammps(object):
 
   # -------------------------------------------------------------------------
 
+  def clearstep_compute(self, nextstep):
+    with ExceptionCheck(self):
+      return self.lib.lammps_clearstep_compute(self.lmp)
+
+  # -------------------------------------------------------------------------
+
+  def addstep_compute(self, nextstep):
+    with ExceptionCheck(self):
+      nextstep = self.c_bigint(nextstep)
+      return self.lib.lammps_addstep_compute(self.lmp, POINTER(nextstep))
+
+  # -------------------------------------------------------------------------
+
+  def addstep_compute_all(self, nextstep):
+    with ExceptionCheck(self):
+      nextstep = self.c_bigint(nextstep)
+      return self.lib.lammps_addstep_compute_all(self.lmp, POINTER(nextstep))
+
+  # -------------------------------------------------------------------------
+
   def flush_buffers(self):
     """Flush output buffers
 

src/.gitignore

@@ -1281,6 +1281,10 @@
 /pair_hbond_dreiding_lj.h
 /pair_hbond_dreiding_morse.cpp
 /pair_hbond_dreiding_morse.h
+/pair_hbond_dreiding_lj_angleoffset.cpp
+/pair_hbond_dreiding_lj_angleoffset.h
+/pair_hbond_dreiding_morse_angleoffset.cpp
+/pair_hbond_dreiding_morse_angleoffset.h
 /pair_hdnnp.cpp
 /pair_hdnnp.h
 /pair_ilp_graphene_hbn.cpp

src/EXTRA-FIX/fix_electron_stopping_fit.cpp

@@ -68,7 +68,8 @@ FixElectronStoppingFit::FixElectronStoppingFit(LAMMPS *lmp, int narg, char **arg
      error->all(FLERR,"Incorrect number of fix electron/stopping/fit arguments");
   }
 
-  scalar_flag = 1;
+  scalar_flag = 1; // intensive total energy loss since start of run
+  extscalar = 0;
   global_freq = 1;
 
   energy_coh_in = new double[atom->ntypes+1];

src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp

@@ -0,0 +1,131 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Tod A Pascal (Caltech), Don Xu/EiPi Fun
+------------------------------------------------------------------------- */
+
+#include "pair_hbond_dreiding_lj_angleoffset.h"
+
+#include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
+#include "error.h"
+#include "force.h"
+#include "math_const.h"
+#include "math_special.h"
+#include "memory.h"
+#include "molecule.h"
+#include "neigh_list.h"
+#include "neighbor.h"
+
+#include <cmath>
+#include <cstring>
+
+static constexpr double SMALL = 0.001;
+static constexpr int CHUNK = 8;
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+// using namespace MathSpecial;
+
+/* ---------------------------------------------------------------------- */
+
+PairHbondDreidingLJAngleoffset::PairHbondDreidingLJAngleoffset(LAMMPS *lmp)
+    : PairHbondDreidingLJ(lmp) {
+
+  angle_offset_flag = 1;
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairHbondDreidingLJAngleoffset::coeff(int narg, char **arg)
+{
+  if (narg < 6 || narg > 11)
+    error->all(FLERR,"Incorrect args for pair coefficients");
+  if (!allocated) allocate();
+
+  int ilo,ihi,jlo,jhi,klo,khi;
+  utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
+  utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
+  utils::bounds_typelabel(FLERR, arg[2], 1, atom->ntypes, klo, khi, lmp, Atom::ATOM);
+
+  int donor_flag;
+  if (strcmp(arg[3],"i") == 0) donor_flag = 0;
+  else if (strcmp(arg[3],"j") == 0) donor_flag = 1;
+  else error->all(FLERR,"Incorrect args for pair coefficients");
+
+  double epsilon_one = utils::numeric(FLERR, arg[4], false, lmp);
+  double sigma_one = utils::numeric(FLERR, arg[5], false, lmp);
+
+  int ap_one = ap_global;
+  if (narg > 6) ap_one = utils::inumeric(FLERR, arg[6], false, lmp);
+  double cut_inner_one = cut_inner_global;
+  double cut_outer_one = cut_outer_global;
+  if (narg > 8) {
+    cut_inner_one = utils::numeric(FLERR, arg[7], false, lmp);
+    cut_outer_one = utils::numeric(FLERR, arg[8], false, lmp);
+  }
+  if (cut_inner_one>cut_outer_one)
+    error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff");
+  double cut_angle_one = cut_angle_global;
+  if (narg > 9) cut_angle_one = utils::numeric(FLERR, arg[9], false, lmp) * MY_PI/180.0;
+  double angle_offset_one = angle_offset_global;
+  if (narg == 11) angle_offset_one = (180.0 - utils::numeric(FLERR, arg[10], false, lmp)) * MY_PI/180.0;
+  if (angle_offset_one < 0.0 || angle_offset_one > 90.0 * MY_PI/180.0)
+    error->all(FLERR,"Illegal angle offset");
+
+  // grow params array if necessary
+
+  if (nparams == maxparam) {
+    maxparam += CHUNK;
+    params = (Param *) memory->srealloc(params, maxparam*sizeof(Param),
+                                        "pair:params");
+
+    // make certain all addional allocated storage is initialized
+    // to avoid false positives when checking with valgrind
+
+    memset(params + nparams, 0, CHUNK*sizeof(Param));
+  }
+
+  params[nparams].epsilon = epsilon_one;
+  params[nparams].sigma = sigma_one;
+  params[nparams].ap = ap_one;
+  params[nparams].cut_inner = cut_inner_one;
+  params[nparams].cut_outer = cut_outer_one;
+  params[nparams].cut_innersq = cut_inner_one*cut_inner_one;
+  params[nparams].cut_outersq = cut_outer_one*cut_outer_one;
+  params[nparams].cut_angle = cut_angle_one;
+  params[nparams].angle_offset = angle_offset_one;
+  params[nparams].denom_vdw =
+    (params[nparams].cut_outersq-params[nparams].cut_innersq) *
+    (params[nparams].cut_outersq-params[nparams].cut_innersq) *
+    (params[nparams].cut_outersq-params[nparams].cut_innersq);
+
+  // flag type2param with either i,j = D,A or j,i = D,A
+
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++)
+    for (int j = MAX(jlo,i); j <= jhi; j++)
+      for (int k = klo; k <= khi; k++) {
+        if (donor_flag == 0) type2param[i][j][k] = nparams;
+        else type2param[j][i][k] = nparams;
+        count++;
+      }
+  nparams++;
+
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}

src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h

@@ -0,0 +1,38 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+// clang-format off
+PairStyle(hbond/dreiding/lj/angleoffset,PairHbondDreidingLJAngleoffset);
+// clang-format on
+#else
+
+#ifndef LMP_PAIR_HBOND_DREIDING_LJ_ANGLEOFFSET_H
+#define LMP_PAIR_HBOND_DREIDING_LJ_ANGLEOFFSET_H
+
+#include "pair_hbond_dreiding_lj.h"
+
+namespace LAMMPS_NS {
+
+class PairHbondDreidingLJAngleoffset : public PairHbondDreidingLJ {
+
+ public:
+  PairHbondDreidingLJAngleoffset(class LAMMPS *);
+  void coeff(int, char **) override;
+
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif

src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp

@@ -0,0 +1,129 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+   ------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Tod A Pascal (Caltech), Don Xu/EiPi Fun
+------------------------------------------------------------------------- */
+
+#include "pair_hbond_dreiding_morse_angleoffset.h"
+
+#include "atom.h"
+#include "atom_vec.h"
+#include "domain.h"
+#include "error.h"
+#include "force.h"
+#include "math_const.h"
+#include "math_special.h"
+#include "memory.h"
+#include "molecule.h"
+#include "neigh_list.h"
+#include "neighbor.h"
+
+#include <cmath>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+static constexpr int CHUNK = 8;
+
+/* ---------------------------------------------------------------------- */
+
+PairHbondDreidingMorseAngleoffset::PairHbondDreidingMorseAngleoffset(LAMMPS *lmp) :
+  PairHbondDreidingMorse(lmp) {
+
+  angle_offset_flag = 1;
+}
+
+/* ----------------------------------------------------------------------
+ *    set coeffs for one or more type pairs
+ * ---------------------------------------------------------------------- */
+
+void PairHbondDreidingMorseAngleoffset::coeff(int narg, char **arg)
+{
+  if (narg < 7 || narg > 12)
+    error->all(FLERR,"Incorrect args for pair coefficients");
+  if (!allocated) allocate();
+
+  int ilo,ihi,jlo,jhi,klo,khi;
+  utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
+  utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
+  utils::bounds_typelabel(FLERR, arg[2], 1, atom->ntypes, klo, khi, lmp, Atom::ATOM);
+
+  int donor_flag;
+  if (strcmp(arg[3],"i") == 0) donor_flag = 0;
+  else if (strcmp(arg[3],"j") == 0) donor_flag = 1;
+  else error->all(FLERR,"Incorrect args for pair coefficients");
+
+  double d0_one = utils::numeric(FLERR, arg[4], false, lmp);
+  double alpha_one = utils::numeric(FLERR, arg[5], false, lmp);
+  double r0_one = utils::numeric(FLERR, arg[6], false, lmp);
+
+  int ap_one = ap_global;
+  if (narg > 7) ap_one = utils::inumeric(FLERR, arg[7], false, lmp);
+  double cut_inner_one = cut_inner_global;
+  double cut_outer_one = cut_outer_global;
+  if (narg > 9) {
+    cut_inner_one = utils::numeric(FLERR, arg[8], false, lmp);
+    cut_outer_one = utils::numeric(FLERR, arg[9], false, lmp);
+  }
+  if (cut_inner_one>cut_outer_one)
+    error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff");
+  double cut_angle_one = cut_angle_global;
+  if (narg > 10) cut_angle_one = utils::numeric(FLERR, arg[10], false, lmp) * MY_PI/180.0;
+  double angle_offset_one = angle_offset_global;
+  if (narg == 12) angle_offset_one = (180.0 - utils::numeric(FLERR, arg[11], false, lmp)) * MY_PI/180.0;
+  if (angle_offset_one < 0.0 || angle_offset_one > 90.0 * MY_PI/180.0)
+    error->all(FLERR,"Illegal angle offset {}", angle_offset_one);
+
+  // grow params array if necessary
+
+  if (nparams == maxparam) {
+    maxparam += CHUNK;
+    params = (Param *) memory->srealloc(params, maxparam*sizeof(Param),"pair:params");
+
+    // make certain all addional allocated storage is initialized
+    // to avoid false positives when checking with valgrind
+
+    memset(params + nparams, 0, CHUNK*sizeof(Param));
+  }
+
+  params[nparams].d0 = d0_one;
+  params[nparams].alpha = alpha_one;
+  params[nparams].r0 = r0_one;
+  params[nparams].ap = ap_one;
+  params[nparams].cut_inner = cut_inner_one;
+  params[nparams].cut_outer = cut_outer_one;
+  params[nparams].cut_innersq = cut_inner_one*cut_inner_one;
+  params[nparams].cut_outersq = cut_outer_one*cut_outer_one;
+  params[nparams].cut_angle = cut_angle_one;
+  params[nparams].angle_offset = angle_offset_one;
+  params[nparams].denom_vdw = (params[nparams].cut_outersq-params[nparams].cut_innersq) *
+    (params[nparams].cut_outersq-params[nparams].cut_innersq) *
+    (params[nparams].cut_outersq-params[nparams].cut_innersq);
+
+  // flag type2param with either i,j = D,A or j,i = D,A
+
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++)
+    for (int j = MAX(jlo,i); j <= jhi; j++)
+      for (int k = klo; k <= khi; k++) {
+        if (donor_flag == 0) type2param[i][j][k] = nparams;
+        else type2param[j][i][k] = nparams;
+        count++;
+      }
+  nparams++;
+
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}

src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h

@@ -0,0 +1,38 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+// clang-format off
+PairStyle(hbond/dreiding/morse/angleoffset,PairHbondDreidingMorseAngleoffset);
+// clang-format on
+#else
+
+#ifndef LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_H
+#define LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_H
+
+#include "pair_hbond_dreiding_morse.h"
+
+namespace LAMMPS_NS {
+
+class PairHbondDreidingMorseAngleoffset : public PairHbondDreidingMorse {
+
+ public:
+  PairHbondDreidingMorseAngleoffset(class LAMMPS *);
+  void coeff(int, char **) override;
+
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif

src/MISC/pair_srp.cpp

@@ -397,10 +397,11 @@ void PairSRP::coeff(int narg, char **arg)
     error->all(FLERR,"PairSRP: Incorrect args for pair coeff");
   if (!allocated) allocate();
 
-  if (btype_str.size() > 0)
+  if (btype_str.size() > 0) {
     btype = utils::expand_type_int(FLERR, btype_str, Atom::BOND, lmp);
     if ((btype > atom->nbondtypes) || (btype <= 0))
       error->all(FLERR,"Invalid bond type {} for pair style srp", btype);
+  }
 
   if (bptype_str.size() > 0)
     bptype = utils::expand_type_int(FLERR, bptype_str, Atom::ATOM, lmp);

src/MOLECULE/pair_hbond_dreiding_lj.cpp

@@ -13,7 +13,7 @@
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
-   Contributing author: Tod A Pascal (Caltech)
+   Contributing authors: Tod A Pascal (Caltech), Don Xu/EiPi Fun
 ------------------------------------------------------------------------- */
 
 #include "pair_hbond_dreiding_lj.h"
@@ -55,6 +55,9 @@ PairHbondDreidingLJ::PairHbondDreidingLJ(LAMMPS *lmp) : Pair(lmp)
 
   nextra = 2;
   pvector = new double[2];
+
+  angle_offset_flag = 0;
+  angle_offset_global = 0.0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -82,7 +85,7 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag)
   tagint tagprev;
   double delx,dely,delz,rsq,rsq1,rsq2,r1,r2;
   double factor_hb,force_angle,force_kernel,evdwl,eng_lj,ehbond,force_switch;
-  double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2,d;
+  double c,s,a,b,d,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
   double fi[3],fj[3],delr1[3],delr2[3];
   double r2inv,r10inv;
   double switch1,switch2;
@@ -178,6 +181,13 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag)
           if (c < -1.0) c = -1.0;
           ac = acos(c);
 
+          if (angle_offset_flag){
+            ac = ac + pm.angle_offset;
+            c = cos(ac);
+            if (c > 1.0) c = 1.0;
+            if (c < -1.0) c = -1.0;
+          }
+
           if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) {
             s = sqrt(1.0 - c*c);
             if (s < SMALL) s = SMALL;
@@ -186,15 +196,12 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag)
 
             r2inv = 1.0/rsq;
             r10inv = r2inv*r2inv*r2inv*r2inv*r2inv;
-            force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv *
+            force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * powint(c,pm.ap);
-              powint(c,pm.ap);
+            force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * powint(c,pm.ap-1)*s;
-            force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) *
+            force_switch = 0.0;
-              powint(c,pm.ap-1)*s;
 
             eng_lj = r10inv*(pm.lj3*r2inv - pm.lj4);
 
-            force_switch=0.0;
-
             if (rsq > pm.cut_innersq) {
               switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) *
                         (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) /
@@ -300,12 +307,19 @@ void PairHbondDreidingLJ::allocate()
 
 void PairHbondDreidingLJ::settings(int narg, char **arg)
 {
-  if (narg != 4) error->all(FLERR,"Illegal pair_style command");
+
+  // narg = 4 for standard form, narg = 5 or 6 if angleoffset LJ or Morse variants respectively (from EXTRA-MOLECULE)
+  if (narg != 4 && narg != 5) error->all(FLERR,"Illegal pair_style command");
 
   ap_global = utils::inumeric(FLERR,arg[0],false,lmp);
   cut_inner_global = utils::numeric(FLERR,arg[1],false,lmp);
   cut_outer_global = utils::numeric(FLERR,arg[2],false,lmp);
   cut_angle_global = utils::numeric(FLERR,arg[3],false,lmp) * MY_PI/180.0;
+
+  // update when using angleoffset variant
+  if (angle_offset_flag) {
+    angle_offset_global = (180.0 - utils::numeric(FLERR, arg[4], false, lmp)) * MY_PI/180.0;
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -314,8 +328,14 @@ void PairHbondDreidingLJ::settings(int narg, char **arg)
 
 void PairHbondDreidingLJ::coeff(int narg, char **arg)
 {
-  if (narg < 6 || narg > 10)
+  // account for angleoffset variant in EXTRA-MOLECULE
+  int maxarg = 10;
+  if (angle_offset_flag == 1) maxarg = 11;
+
+  // check settings
+  if (narg < 6 || narg > maxarg)
     error->all(FLERR,"Incorrect args for pair coefficients");
+
   if (!allocated) allocate();
 
   int ilo,ihi,jlo,jhi,klo,khi;
@@ -343,16 +363,15 @@ void PairHbondDreidingLJ::coeff(int narg, char **arg)
     error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff");
   double cut_angle_one = cut_angle_global;
   if (narg == 10) cut_angle_one = utils::numeric(FLERR, arg[9], false, lmp) * MY_PI/180.0;
+
   // grow params array if necessary
 
   if (nparams == maxparam) {
     maxparam += CHUNK;
-    params = (Param *) memory->srealloc(params, maxparam*sizeof(Param),
+    params = (Param *) memory->srealloc(params, maxparam*sizeof(Param), "pair:params");
-                                        "pair:params");
 
     // make certain all addional allocated storage is initialized
     // to avoid false positives when checking with valgrind
-
     memset(params + nparams, 0, CHUNK*sizeof(Param));
   }
 
@@ -540,6 +559,13 @@ double PairHbondDreidingLJ::single(int i, int j, int itype, int jtype,
     if (c < -1.0) c = -1.0;
     ac = acos(c);
 
+    if (angle_offset_flag){
+      ac = ac + pm.angle_offset;
+      c = cos(ac);
+      if (c > 1.0) c = 1.0;
+      if (c < -1.0) c = -1.0;
+    }
+
     if (ac < pm.cut_angle || ac > (2.0*MY_PI - pm.cut_angle)) return 0.0;
     s = sqrt(1.0 - c*c);
     if (s < SMALL) s = SMALL;

src/MOLECULE/pair_hbond_dreiding_lj.h

@@ -25,6 +25,7 @@ PairStyle(hbond/dreiding/lj,PairHbondDreidingLJ);
 namespace LAMMPS_NS {
 
 class PairHbondDreidingLJ : public Pair {
+
  public:
   PairHbondDreidingLJ(class LAMMPS *);
   ~PairHbondDreidingLJ() override;
@@ -38,6 +39,8 @@ class PairHbondDreidingLJ : public Pair {
  protected:
   double cut_inner_global, cut_outer_global, cut_angle_global;
   int ap_global;
+  int angle_offset_flag;          // 1 if angle offset variant used
+  double angle_offset_global;          // updated if angle offset variant used
 
   struct Param {
     double epsilon, sigma;
@@ -45,7 +48,7 @@ class PairHbondDreidingLJ : public Pair {
     double d0, alpha, r0;
     double morse1;
     double denom_vdw;
-    double cut_inner, cut_outer, cut_innersq, cut_outersq, cut_angle, offset;
+    double cut_inner, cut_outer, cut_innersq, cut_outersq, cut_angle, offset, angle_offset;
     int ap;
   };
 

src/MOLECULE/pair_hbond_dreiding_morse.cpp

@@ -148,6 +148,13 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
           if (c < -1.0) c = -1.0;
           ac = acos(c);
 
+          if (angle_offset_flag){
+            ac = ac + pm.angle_offset;
+            c = cos(ac);
+            if (c > 1.0) c = 1.0;
+            if (c < -1.0) c = -1.0;
+          }
+
           if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) {
             s = sqrt(1.0 - c*c);
             if (s < SMALL) s = SMALL;
@@ -158,6 +165,7 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
             dr = r - pm.r0;
             dexp = exp(-pm.alpha * dr);
             eng_morse = pm.d0 * (dexp*dexp - 2.0*dexp);
+
             force_kernel = pm.morse1*(dexp*dexp - dexp)/r * powint(c,pm.ap);
             force_angle = pm.ap * eng_morse * powint(c,pm.ap-1)*s;
             force_switch = 0.0;
@@ -196,12 +204,12 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
             vz1 = a11*delr1[2] + a12*delr2[2];
             vz2 = a22*delr2[2] + a12*delr1[2];
 
-            fi[0] = vx1 + (b+d)*delx;
+            fi[0] = vx1 + b*delx + d*delx;
-            fi[1] = vy1 + (b+d)*dely;
+            fi[1] = vy1 + b*dely + d*dely;
-            fi[2] = vz1 + (b+d)*delz;
+            fi[2] = vz1 + b*delz + d*delz;
-            fj[0] = vx2 - (b+d)*delx;
+            fj[0] = vx2 - b*delx - d*delx;
-            fj[1] = vy2 - (b+d)*dely;
+            fj[1] = vy2 - b*dely - d*dely;
-            fj[2] = vz2 - (b+d)*delz;
+            fj[2] = vz2 - b*delz - d*delz;
 
             f[i][0] += fi[0];
             f[i][1] += fi[1];
@@ -238,7 +246,9 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
 
 void PairHbondDreidingMorse::coeff(int narg, char **arg)
 {
-  if (narg < 7 || narg > 11)
+  int maxarg = 12;
+  if (angle_offset_flag == 1) maxarg = 12;
+  if (narg < 7 || narg > maxarg)
     error->all(FLERR,"Incorrect args for pair coefficients");
   if (!allocated) allocate();
 
@@ -443,6 +453,13 @@ double PairHbondDreidingMorse::single(int i, int j, int itype, int jtype,
     if (c < -1.0) c = -1.0;
     ac = acos(c);
 
+    if (angle_offset_flag){
+      ac = ac + pm.angle_offset;
+      c = cos(ac);
+      if (c > 1.0) c = 1.0;
+      if (c < -1.0) c = -1.0;
+    }
+
     if (ac < pm.cut_angle || ac > (2.0*MY_PI - pm.cut_angle)) return 0.0;
     s = sqrt(1.0 - c*c);
     if (s < SMALL) s = SMALL;

src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp

@@ -0,0 +1,127 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   This software is distributed under the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Axel Kohlmeyer (Temple U), Don Xu/EiPi Fun
+------------------------------------------------------------------------- */
+
+#include "pair_hbond_dreiding_lj_angleoffset_omp.h"
+
+#include "atom.h"
+#include "atom_vec.h"
+#include "comm.h"
+#include "domain.h"
+#include "force.h"
+#include "math_const.h"
+#include "math_special.h"
+#include "memory.h"
+#include "molecule.h"
+#include "neigh_list.h"
+#include "suffix.h"
+
+#include <cmath>
+
+#include "omp_compat.h"
+using namespace LAMMPS_NS;
+using namespace MathConst;
+using namespace MathSpecial;
+
+static constexpr int CHUNK = 8;
+
+/* ---------------------------------------------------------------------- */
+
+PairHbondDreidingLJAngleoffsetOMP::PairHbondDreidingLJAngleoffsetOMP(LAMMPS *lmp) :
+  PairHbondDreidingLJOMP(lmp) {
+  angle_offset_flag = 1;
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairHbondDreidingLJAngleoffsetOMP::coeff(int narg, char **arg)
+{
+  auto mylmp = PairHbondDreidingLJ::lmp;
+  if (narg < 6 || narg > 11)
+    error->all(FLERR,"Incorrect args for pair coefficients");
+  if (!allocated) allocate();
+
+  int ilo,ihi,jlo,jhi,klo,khi;
+  utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
+  utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
+  utils::bounds_typelabel(FLERR, arg[2], 1, atom->ntypes, klo, khi, mylmp, Atom::ATOM);
+
+  int donor_flag;
+  if (strcmp(arg[3],"i") == 0) donor_flag = 0;
+  else if (strcmp(arg[3],"j") == 0) donor_flag = 1;
+  else error->all(FLERR,"Incorrect args for pair coefficients");
+
+  double epsilon_one = utils::numeric(FLERR, arg[4], false, mylmp);
+  double sigma_one = utils::numeric(FLERR, arg[5], false, mylmp);
+
+  int ap_one = ap_global;
+  if (narg > 6) ap_one = utils::inumeric(FLERR, arg[6], false, mylmp);
+  double cut_inner_one = cut_inner_global;
+  double cut_outer_one = cut_outer_global;
+  if (narg > 8) {
+    cut_inner_one = utils::numeric(FLERR, arg[7], false, mylmp);
+    cut_outer_one = utils::numeric(FLERR, arg[8], false, mylmp);
+  }
+  if (cut_inner_one>cut_outer_one)
+    error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff");
+  double cut_angle_one = cut_angle_global;
+  if (narg > 9) cut_angle_one = utils::numeric(FLERR, arg[9], false, mylmp) * MY_PI/180.0;
+  double angle_offset_one = angle_offset_global;
+  if (narg == 11) angle_offset_one = (180.0 - utils::numeric(FLERR, arg[10], false, mylmp)) * MY_PI/180.0;
+  if (angle_offset_one < 0.0 || angle_offset_one > 90.0 * MY_PI/180.0)
+    error->all(FLERR,"Illegal angle offset");
+
+  // grow params array if necessary
+
+  if (nparams == maxparam) {
+    maxparam += CHUNK;
+    params = (Param *) memory->srealloc(params, maxparam*sizeof(Param),
+                                        "pair:params");
+
+    // make certain all addional allocated storage is initialized
+    // to avoid false positives when checking with valgrind
+
+    memset(params + nparams, 0, CHUNK*sizeof(Param));
+  }
+
+  params[nparams].epsilon = epsilon_one;
+  params[nparams].sigma = sigma_one;
+  params[nparams].ap = ap_one;
+  params[nparams].cut_inner = cut_inner_one;
+  params[nparams].cut_outer = cut_outer_one;
+  params[nparams].cut_innersq = cut_inner_one*cut_inner_one;
+  params[nparams].cut_outersq = cut_outer_one*cut_outer_one;
+  params[nparams].cut_angle = cut_angle_one;
+  params[nparams].angle_offset = angle_offset_one;
+  params[nparams].denom_vdw =
+    (params[nparams].cut_outersq-params[nparams].cut_innersq) *
+    (params[nparams].cut_outersq-params[nparams].cut_innersq) *
+    (params[nparams].cut_outersq-params[nparams].cut_innersq);
+
+  // flag type2param with either i,j = D,A or j,i = D,A
+
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++)
+    for (int j = MAX(jlo,i); j <= jhi; j++)
+      for (int k = klo; k <= khi; k++) {
+        if (donor_flag == 0) type2param[i][j][k] = nparams;
+        else type2param[j][i][k] = nparams;
+        count++;
+      }
+  nparams++;
+
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}

src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h

@@ -0,0 +1,41 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Axel Kohlmeyer (Temple U), Don Xu/EiPi Fun
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+// clang-format off
+PairStyle(hbond/dreiding/lj/angleoffset/omp,PairHbondDreidingLJAngleoffsetOMP);
+// clang-format on
+#else
+
+#ifndef LMP_PAIR_HBOND_DREIDING_LJ_ANGLEOFFSET_OMP_H
+#define LMP_PAIR_HBOND_DREIDING_LJ_ANGLEOFFSET_OMP_H
+
+#include "pair_hbond_dreiding_lj_omp.h"
+
+namespace LAMMPS_NS {
+
+class PairHbondDreidingLJAngleoffsetOMP : public PairHbondDreidingLJOMP {
+
+ public:
+  PairHbondDreidingLJAngleoffsetOMP(class LAMMPS *);
+  void coeff(int, char **) override;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif

src/OPENMP/pair_hbond_dreiding_lj_omp.cpp

@@ -120,15 +120,16 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr)
   int i,j,k,m,ii,jj,kk,jnum,knum,itype,jtype,ktype,iatom,imol;
   tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq,rsq1,rsq2,r1,r2;
-  double factor_hb,force_angle,force_kernel,evdwl,eng_lj;
+  double factor_hb,force_angle,force_kernel,evdwl,eng_lj,ehbond,force_switch;
-  double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
+  double c,s,a,b,d,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
   double fi[3],fj[3],delr1[3],delr2[3];
   double r2inv,r10inv;
   double switch1,switch2;
   int *ilist,*jlist,*numneigh,**firstneigh;
   const tagint *klist;
 
-  evdwl = 0.0;
+  evdwl = ehbond = 0.0;
+  int hbcount = 0;
 
   const auto * _noalias const x = (dbl3_t *) atom->x[0];
   auto * _noalias const f = (dbl3_t *) thr->get_f()[0];
@@ -151,9 +152,6 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr)
   // jj = loop over acceptors
   // kk = loop over hydrogens bonded to donor
 
-  int hbcount = 0;
-  double hbeng = 0.0;
-
   for (ii = iifrom; ii < iito; ++ii) {
     i = ilist[ii];
     itype = type[i];
@@ -222,6 +220,13 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr)
           if (c < -1.0) c = -1.0;
           ac = acos(c);
 
+          if (angle_offset_flag){
+            ac = ac + pm.angle_offset;
+            c = cos(ac);
+            if (c > 1.0) c = 1.0;
+            if (c < -1.0) c = -1.0;
+          }
+
           if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) {
             s = sqrt(1.0 - c*c);
             if (s < SMALL) s = SMALL;
@@ -230,30 +235,34 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr)
 
             r2inv = 1.0/rsq;
             r10inv = r2inv*r2inv*r2inv*r2inv*r2inv;
-            force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv *
+            force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * powint(c,pm.ap);
-              powint(c,pm.ap);
+            force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * powint(c,pm.ap-1)*s;
-            force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) *
+            force_switch = 0.0;
-              powint(c,pm.ap-1)*s;
 
             eng_lj = r10inv*(pm.lj3*r2inv - pm.lj4);
+
             if (rsq > pm.cut_innersq) {
               switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) *
                         (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) /
                         pm.denom_vdw;
               switch2 = 12.0*rsq * (pm.cut_outersq-rsq) *
                         (rsq-pm.cut_innersq) / pm.denom_vdw;
-              force_kernel = force_kernel*switch1 + eng_lj*switch2/rsq;
+
+              force_kernel *= switch1;
               force_angle  *= switch1;
+              force_switch  = eng_lj*switch2/rsq;
               eng_lj       *= switch1;
             }
 
             if (EFLAG) {
               evdwl = eng_lj * powint(c,pm.ap);
               evdwl *= factor_hb;
+              ehbond += evdwl;
             }
 
             a = factor_hb*force_angle/s;
             b = factor_hb*force_kernel;
+            d = factor_hb*force_switch;
 
             a11 = a*c / rsq1;
             a12 = -a / (r1*r2);
@@ -266,12 +275,12 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr)
             vz1 = a11*delr1[2] + a12*delr2[2];
             vz2 = a22*delr2[2] + a12*delr1[2];
 
-            fi[0] = vx1 + b*delx;
+            fi[0] = vx1 + b*delx + d*delx;
-            fi[1] = vy1 + b*dely;
+            fi[1] = vy1 + b*dely + d*dely;
-            fi[2] = vz1 + b*delz;
+            fi[2] = vz1 + b*delz + d*delz;
-            fj[0] = vx2 - b*delx;
+            fj[0] = vx2 - b*delx - d*delx;
-            fj[1] = vy2 - b*dely;
+            fj[1] = vy2 - b*dely - d*dely;
-            fj[2] = vz2 - b*delz;
+            fj[2] = vz2 - b*delz - d*delz;
 
             fxtmp += fi[0];
             fytmp += fi[1];
@@ -290,7 +299,7 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr)
             if (EVFLAG) ev_tally3_thr(this,k,i,j,evdwl,0.0,fi,fj,delr1,delr2,thr);
             if (EFLAG) {
               hbcount++;
-              hbeng += evdwl;
+              ehbond += evdwl;
             }
           }
         }
@@ -302,7 +311,7 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr)
   }
   const int tid = thr->get_tid();
   hbcount_thr[tid] = static_cast<double>(hbcount);
-  hbeng_thr[tid] = hbeng;
+  hbeng_thr[tid] = ehbond;
 }
 
 /* ---------------------------------------------------------------------- */

src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp

@@ -0,0 +1,127 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   This software is distributed under the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Axel Kohlmeyer (Temple U), Don Xu/EiPi Fun
+------------------------------------------------------------------------- */
+
+#include "pair_hbond_dreiding_morse_angleoffset_omp.h"
+
+#include "atom.h"
+#include "atom_vec.h"
+#include "comm.h"
+#include "domain.h"
+#include "force.h"
+#include "math_const.h"
+#include "math_special.h"
+#include "memory.h"
+#include "molecule.h"
+#include "neigh_list.h"
+#include "suffix.h"
+
+#include <cmath>
+
+#include "omp_compat.h"
+using namespace LAMMPS_NS;
+using namespace MathConst;
+using namespace MathSpecial;
+
+static constexpr int CHUNK = 8;
+
+/* ---------------------------------------------------------------------- */
+
+PairHbondDreidingMorseAngleoffsetOMP::PairHbondDreidingMorseAngleoffsetOMP(LAMMPS *lmp) :
+  PairHbondDreidingMorseOMP(lmp) {
+  angle_offset_flag = 1;
+}
+
+/* ----------------------------------------------------------------------
+ *    set coeffs for one or more type pairs
+ * ---------------------------------------------------------------------- */
+
+void PairHbondDreidingMorseAngleoffsetOMP::coeff(int narg, char **arg)
+{
+  auto mylmp = PairHbondDreidingMorse::lmp;
+  if (narg < 7 || narg > 12)
+    error->all(FLERR,"Incorrect args for pair coefficients");
+  if (!allocated) allocate();
+
+  int ilo,ihi,jlo,jhi,klo,khi;
+  utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
+  utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
+  utils::bounds_typelabel(FLERR, arg[2], 1, atom->ntypes, klo, khi, mylmp, Atom::ATOM);
+
+  int donor_flag;
+  if (strcmp(arg[3],"i") == 0) donor_flag = 0;
+  else if (strcmp(arg[3],"j") == 0) donor_flag = 1;
+  else error->all(FLERR,"Incorrect args for pair coefficients");
+
+  double d0_one = utils::numeric(FLERR, arg[4], false, mylmp);
+  double alpha_one = utils::numeric(FLERR, arg[5], false, mylmp);
+  double r0_one = utils::numeric(FLERR, arg[6], false, mylmp);
+
+  int ap_one = ap_global;
+  if (narg > 7) ap_one = utils::inumeric(FLERR, arg[7], false, mylmp);
+  double cut_inner_one = cut_inner_global;
+  double cut_outer_one = cut_outer_global;
+  if (narg > 9) {
+    cut_inner_one = utils::numeric(FLERR, arg[8], false, mylmp);
+    cut_outer_one = utils::numeric(FLERR, arg[9], false, mylmp);
+  }
+  if (cut_inner_one>cut_outer_one)
+    error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff");
+  double cut_angle_one = cut_angle_global;
+  if (narg > 10) cut_angle_one = utils::numeric(FLERR, arg[10], false, mylmp) * MY_PI/180.0;
+  double angle_offset_one = angle_offset_global;
+  if (narg == 12) angle_offset_one = (180.0 - utils::numeric(FLERR, arg[11], false, mylmp)) * MY_PI/180.0;
+  if (angle_offset_one < 0.0 || angle_offset_one > 90.0 * MY_PI/180.0)
+    error->all(FLERR,"Illegal angle offset {}", angle_offset_one);
+
+  // grow params array if necessary
+
+  if (nparams == maxparam) {
+    maxparam += CHUNK;
+    params = (Param *) memory->srealloc(params, maxparam*sizeof(Param),"pair:params");
+
+    // make certain all addional allocated storage is initialized
+    // to avoid false positives when checking with valgrind
+
+    memset(params + nparams, 0, CHUNK*sizeof(Param));
+  }
+
+  params[nparams].d0 = d0_one;
+  params[nparams].alpha = alpha_one;
+  params[nparams].r0 = r0_one;
+  params[nparams].ap = ap_one;
+  params[nparams].cut_inner = cut_inner_one;
+  params[nparams].cut_outer = cut_outer_one;
+  params[nparams].cut_innersq = cut_inner_one*cut_inner_one;
+  params[nparams].cut_outersq = cut_outer_one*cut_outer_one;
+  params[nparams].cut_angle = cut_angle_one;
+  params[nparams].angle_offset = angle_offset_one;
+  params[nparams].denom_vdw = (params[nparams].cut_outersq-params[nparams].cut_innersq) *
+    (params[nparams].cut_outersq-params[nparams].cut_innersq) *
+    (params[nparams].cut_outersq-params[nparams].cut_innersq);
+
+  // flag type2param with either i,j = D,A or j,i = D,A
+
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++)
+    for (int j = MAX(jlo,i); j <= jhi; j++)
+      for (int k = klo; k <= khi; k++) {
+        if (donor_flag == 0) type2param[i][j][k] = nparams;
+        else type2param[j][i][k] = nparams;
+        count++;
+      }
+  nparams++;
+
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}

src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h

@@ -0,0 +1,42 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Axel Kohlmeyer (Temple U), Don Xu/EiPi Fun
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+// clang-format off
+PairStyle(hbond/dreiding/morse/angleoffset/omp,PairHbondDreidingMorseAngleoffsetOMP);
+// clang-format on
+#else
+
+#ifndef LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_OMP_H
+#define LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_OMP_H
+
+#include "pair_hbond_dreiding_morse_omp.h"
+
+namespace LAMMPS_NS {
+
+class PairHbondDreidingMorseAngleoffsetOMP :
+    public PairHbondDreidingMorseOMP {
+
+ public:
+  PairHbondDreidingMorseAngleoffsetOMP(class LAMMPS *);
+  void coeff(int, char **) override;
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif

src/OPENMP/pair_hbond_dreiding_morse_omp.cpp

@@ -120,14 +120,14 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr)
   int i,j,k,m,ii,jj,kk,jnum,knum,itype,jtype,ktype,imol,iatom;
   tagint tagprev;
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq,rsq1,rsq2,r1,r2;
-  double factor_hb,force_angle,force_kernel,evdwl;
+  double factor_hb,force_angle,force_kernel,force_switch,evdwl,ehbond;
-  double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
+  double c,s,a,b,d,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
   double fi[3],fj[3],delr1[3],delr2[3];
   double r,dr,dexp,eng_morse,switch1,switch2;
   int *ilist,*jlist,*numneigh,**firstneigh;
   const tagint *klist;
 
-  evdwl = 0.0;
+  evdwl = ehbond = 0.0;
 
   const auto * _noalias const x = (dbl3_t *) atom->x[0];
   auto * _noalias const f = (dbl3_t *) thr->get_f()[0];
@@ -151,7 +151,6 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr)
   // kk = loop over hydrogens bonded to donor
 
   int hbcount = 0;
-  double hbeng = 0.0;
 
   for (ii = iifrom; ii < iito; ++ii) {
 
@@ -222,6 +221,13 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr)
           if (c < -1.0) c = -1.0;
           ac = acos(c);
 
+          if (angle_offset_flag){
+            ac = ac + pm.angle_offset;
+            c = cos(ac);
+            if (c > 1.0) c = 1.0;
+            if (c < -1.0) c = -1.0;
+          }
+
           if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) {
             s = sqrt(1.0 - c*c);
             if (s < SMALL) s = SMALL;
@@ -232,8 +238,10 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr)
             dr = r - pm.r0;
             dexp = exp(-pm.alpha * dr);
             eng_morse = pm.d0 * (dexp*dexp - 2.0*dexp);
+
             force_kernel = pm.morse1*(dexp*dexp - dexp)/r * powint(c,pm.ap);
             force_angle = pm.ap * eng_morse * powint(c,pm.ap-1)*s;
+            force_switch = 0.0;
 
             if (rsq > pm.cut_innersq) {
               switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) *
@@ -241,18 +249,22 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr)
                         pm.denom_vdw;
               switch2 = 12.0*rsq * (pm.cut_outersq-rsq) *
                         (rsq-pm.cut_innersq) / pm.denom_vdw;
-              force_kernel = force_kernel*switch1 + eng_morse*switch2/rsq;
+
+              force_kernel *= switch1;
               force_angle  *= switch1;
+              force_switch  = eng_morse*switch2/rsq;
               eng_morse    *= switch1;
             }
 
             if (EFLAG) {
               evdwl = eng_morse * powint(c,pm.ap);
               evdwl *= factor_hb;
+              ehbond += evdwl;
             }
 
             a = factor_hb*force_angle/s;
             b = factor_hb*force_kernel;
+            d = factor_hb*force_switch;
 
             a11 = a*c / rsq1;
             a12 = -a / (r1*r2);
@@ -265,12 +277,12 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr)
             vz1 = a11*delr1[2] + a12*delr2[2];
             vz2 = a22*delr2[2] + a12*delr1[2];
 
-            fi[0] = vx1 + b*delx;
+            fi[0] = vx1 + b*delx + d*delx;
-            fi[1] = vy1 + b*dely;
+            fi[1] = vy1 + b*dely + d*dely;
-            fi[2] = vz1 + b*delz;
+            fi[2] = vz1 + b*delz + d*delz;
-            fj[0] = vx2 - b*delx;
+            fj[0] = vx2 - b*delx - d*delx;
-            fj[1] = vy2 - b*dely;
+            fj[1] = vy2 - b*dely - d*dely;
-            fj[2] = vz2 - b*delz;
+            fj[2] = vz2 - b*delz - d*delz;
 
             fxtmp += fi[0];
             fytmp += fi[1];
@@ -289,7 +301,7 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr)
             if (EVFLAG) ev_tally3_thr(this,k,i,j,evdwl,0.0,fi,fj,delr1,delr2,thr);
             if (EFLAG) {
               hbcount++;
-              hbeng += evdwl;
+              ehbond += evdwl;
             }
           }
         }
@@ -301,7 +313,7 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr)
   }
   const int tid = thr->get_tid();
   hbcount_thr[tid] = static_cast<double>(hbcount);
-  hbeng_thr[tid] = hbeng;
+  hbeng_thr[tid] = ehbond;
 }
 
 /* ---------------------------------------------------------------------- */

src/PYTHON/fix_python_invoke.cpp

@@ -22,6 +22,7 @@
 #include "lmppython.h"
 #include "python_compat.h"
 #include "python_utils.h"
+#include "modify.h"
 #include "update.h"
 
 #include <cstring>
@@ -70,6 +71,12 @@ FixPythonInvoke::FixPythonInvoke(LAMMPS *lmp, int narg, char **arg) :
   }
 
   lmpPtr = PY_VOID_POINTER(lmp);
+
+  // nvalid = next step on which end_of_step or post_force does something
+  // add nextvalid() to all computes that store invocation times
+  // since we don't know a priori which are invoked by python code
+  nvalid = nextvalid();
+  modify->addstep_compute_all(nvalid);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -89,8 +96,23 @@ int FixPythonInvoke::setmask()
 
 /* ---------------------------------------------------------------------- */
 
+void FixPythonInvoke::init()
+{
+  // need to reset nvalid if nvalid < ntimestep b/c minimize was performed
+
+  if (nvalid < update->ntimestep) {
+    nvalid = nextvalid();
+    modify->addstep_compute_all(nvalid);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixPythonInvoke::end_of_step()
 {
+  // python code may invoke computes so wrap with clear/add
+  modify->clearstep_compute();
+
   PyUtils::GIL lock;
 
   PyObject * result = PyObject_CallFunction((PyObject*)pFunc, (char *)"O", (PyObject*)lmpPtr);
@@ -101,6 +123,9 @@ void FixPythonInvoke::end_of_step()
   }
 
   Py_CLEAR(result);
+
+  nvalid = nextvalid();
+  modify->addstep_compute(nvalid);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -116,6 +141,9 @@ void FixPythonInvoke::post_force(int vflag)
 {
   if (update->ntimestep % nevery != 0) return;
 
+  // python code may invoke computes so wrap with clear/add
+  modify->clearstep_compute();
+
   PyUtils::GIL lock;
   char fmt[] = "Oi";
 
@@ -127,4 +155,14 @@ void FixPythonInvoke::post_force(int vflag)
   }
 
   Py_CLEAR(result);
+
+  nvalid = nextvalid();
+  modify->addstep_compute(nvalid);
+}
+
+/* ---------------------------------------------------------------------- */
+
+bigint FixPythonInvoke::nextvalid()
+{
+  return (update->ntimestep/nevery + 1)*nevery;
 }

src/PYTHON/fix_python_invoke.h

@@ -31,6 +31,7 @@ class FixPythonInvoke : public Fix {
   ~FixPythonInvoke() override;
   int setmask() override;
   void setup(int) override;
+  void init() override;
   void end_of_step() override;
   void post_force(int) override;
 
@@ -38,6 +39,8 @@ class FixPythonInvoke : public Fix {
   void *lmpPtr;
   void *pFunc;
   int selected_callback;
+  bigint nextvalid();
+  bigint nvalid;
 };
 
 }    // namespace LAMMPS_NS

src/REAXFF/fix_reaxff_species.cpp

@@ -26,6 +26,7 @@
 #include "domain.h"
 #include "error.h"
 #include "fix_ave_atom.h"
+#include "fix_property_atom.h"
 #include "force.h"
 #include "group.h"
 #include "input.h"
@@ -141,13 +142,6 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) :
   }
 
   x0 = nullptr;
-  clusterID = nullptr;
-
-  int ntmp = atom->nmax;
-  memory->create(x0, ntmp, "reaxff/species:x0");
-  memory->create(clusterID, ntmp, "reaxff/species:clusterID");
-  memset(clusterID, 0, sizeof(double) * ntmp);
-  vector_atom = clusterID;
 
   nmax = 0;
   setupflag = 0;
@@ -304,7 +298,6 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) :
 FixReaxFFSpecies::~FixReaxFFSpecies()
 {
   memory->destroy(BOCut);
-  memory->destroy(clusterID);
   memory->destroy(x0);
 
   memory->destroy(nd);
@@ -330,6 +323,7 @@ FixReaxFFSpecies::~FixReaxFFSpecies()
   try {
     modify->delete_compute(fmt::format("SPECATOM_{}", id));
     modify->delete_fix(fmt::format("SPECBOND_{}", id));
+    modify->delete_fix(fmt::format("clusterID_{}", id));
   } catch (std::exception &) {
   }
 }
@@ -372,9 +366,6 @@ void FixReaxFFSpecies::init()
 
   reaxff->fixspecies_flag = 1;
 
-  // reset next output timestep if not yet set or timestep has been reset
-  if (nvalid != update->ntimestep) nvalid = update->ntimestep + nfreq;
-
   if (!setupflag) {
     // create a compute to store properties
     modify->add_compute(fmt::format("SPECATOM_{} all SPEC/ATOM q x y z vx vy vz abo01 abo02 "
@@ -387,11 +378,25 @@ void FixReaxFFSpecies::init()
     auto fixcmd = fmt::format("SPECBOND_{} all ave/atom {} {} {}", id, nevery, nrepeat, nfreq);
     for (int i = 1; i < 32; ++i) fixcmd += fmt::format(" c_SPECATOM_{}[{}]", id, i);
     f_SPECBOND = dynamic_cast<FixAveAtom *>(modify->add_fix(fixcmd));
+
+    // create a fix to point to fix_property_atom for storing clusterID
+    fixcmd = fmt::format("clusterID_{} all property/atom d_clusterID ghost yes", id);
+    f_clusterID = dynamic_cast<FixPropertyAtom *>(modify->add_fix(fixcmd));
+
+    // per-atom property for clusterID
+    int flag,cols;
+    int index1 = atom->find_custom("clusterID",flag,cols);
+    clusterID = atom->dvector[index1];
+    vector_atom = clusterID;
+
+    int ntmp = atom->nmax;
+    memory->create(x0, ntmp, "reaxff/species:x0");
+
     setupflag = 1;
   }
 
   // check for valid variable name for delete Nlimit keyword
-  if (delete_Nsteps > 0) {
+  if (delete_Nsteps > 0 && delete_Nlimit_varid > -1) {
     delete_Nlimit_varid = input->variable->find(delete_Nlimit_varname.c_str());
     if (delete_Nlimit_varid < 0)
       error->all(FLERR, "Fix reaxff/species: Variable name {} does not exist",
@@ -424,10 +429,16 @@ void FixReaxFFSpecies::Output_ReaxFF_Bonds(bigint ntimestep, FILE * /*fp*/)
 {
   int Nmole, Nspec;
 
+  // per-atom property for clusterID
+  int flag,cols;
+  int index1 = atom->find_custom("clusterID",flag,cols);
+  clusterID = atom->dvector[index1];
+  vector_atom = clusterID;
+
   // point to fix_ave_atom
   f_SPECBOND->end_of_step();
 
-  if (ntimestep != nvalid) {
+  if (ntimestep != nvalid && nvalid != -1) {
     // push back delete_Tcount on every step
     if (delete_Nsteps > 0)
       for (int i = delete_Nsteps - 1; i > 0; i--) delete_Tcount[i] = delete_Tcount[i - 1];
@@ -439,11 +450,7 @@ void FixReaxFFSpecies::Output_ReaxFF_Bonds(bigint ntimestep, FILE * /*fp*/)
   if (atom->nmax > nmax) {
     nmax = atom->nmax;
     memory->destroy(x0);
-    memory->destroy(clusterID);
     memory->create(x0, nmax, "reaxff/species:x0");
-    memory->create(clusterID, nmax, "reaxff/species:clusterID");
-    memset(clusterID, 0, sizeof(double) * nmax);
-    vector_atom = clusterID;
   }
 
   for (int i = 0; i < nmax; i++) { x0[i].x = x0[i].y = x0[i].z = 0.0; }
@@ -464,9 +471,14 @@ void FixReaxFFSpecies::Output_ReaxFF_Bonds(bigint ntimestep, FILE * /*fp*/)
     if (comm->me == 0) fflush(pos);
   }
 
-  if (delflag) DeleteSpecies(Nmole, Nspec);
+  if (delflag && nvalid != -1) {
+    DeleteSpecies(Nmole, Nspec);
 
-  nvalid += nfreq;
+    // reset molecule ID to index from 1
+    SortMolecule(Nmole);
+  }
+
+  nvalid = ntimestep + nfreq;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -826,7 +838,8 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec)
   int count, count_tmp, m, n, k;
   int *Nameall;
   int *mask = atom->mask;
-  double avq, avq_tmp, avx[3], avx_tmp, box[3], halfbox[3];
+  double *rmass = atom->rmass;
+  double totq, totq_tmp, com[3], com_tmp, thism, totm, box[3], halfbox[3];
   double **spec_atom = f_SPECBOND->array_atom;
 
   if (multipos) OpenPos();
@@ -844,7 +857,7 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec)
                update->ntimestep, Nmole, Nspec, domain->boxlo[0], domain->boxhi[0],
                domain->boxlo[1], domain->boxhi[1], domain->boxlo[2], domain->boxhi[2]);
 
-    fprintf(pos, "ID\tAtom_Count\tType\tAve_q\t\tCoM_x\t\tCoM_y\t\tCoM_z\n");
+    fprintf(pos, "ID\tAtom_Count\tType\tTot_q\t\tCoM_x\t\tCoM_y\t\tCoM_z\n");
   }
 
   Nameall = nullptr;
@@ -853,8 +866,9 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec)
   for (m = 1; m <= Nmole; m++) {
 
     count = 0;
-    avq = 0.0;
+    totq = 0.0;
-    for (n = 0; n < 3; n++) avx[n] = 0.0;
+    totm = 0.0;
+    for (n = 0; n < 3; n++) com[n] = 0.0;
     for (n = 0; n < nutypes; n++) Name[n] = 0;
 
     for (i = 0; i < nlocal; i++) {
@@ -864,30 +878,37 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec)
         itype = ele2uele[atom->type[i] - 1];
         Name[itype]++;
         count++;
-        avq += spec_atom[i][0];
+        totq += spec_atom[i][0];
         if ((x0[i].x - spec_atom[i][1]) > halfbox[0]) spec_atom[i][1] += box[0];
         if ((spec_atom[i][1] - x0[i].x) > halfbox[0]) spec_atom[i][1] -= box[0];
         if ((x0[i].y - spec_atom[i][2]) > halfbox[1]) spec_atom[i][2] += box[1];
         if ((spec_atom[i][2] - x0[i].y) > halfbox[1]) spec_atom[i][2] -= box[1];
         if ((x0[i].z - spec_atom[i][3]) > halfbox[2]) spec_atom[i][3] += box[2];
         if ((spec_atom[i][3] - x0[i].z) > halfbox[2]) spec_atom[i][3] -= box[2];
-        for (n = 0; n < 3; n++) avx[n] += spec_atom[i][n + 1];
+        if (rmass) thism = rmass[i];
+        else thism = atom->mass[atom->type[i]];
+        for (n = 0; n < 3; n++) com[n] += spec_atom[i][n+1]*thism;
+        totm += thism;
       }
     }
 
-    avq_tmp = 0.0;
+    totq_tmp = 0.0;
-    MPI_Allreduce(&avq, &avq_tmp, 1, MPI_DOUBLE, MPI_SUM, world);
+    MPI_Allreduce(&totq, &totq_tmp, 1, MPI_DOUBLE, MPI_SUM, world);
-    avq = avq_tmp;
+    totq = totq_tmp;
 
     for (n = 0; n < 3; n++) {
-      avx_tmp = 0.0;
+      com_tmp = 0.0;
-      MPI_Reduce(&avx[n], &avx_tmp, 1, MPI_DOUBLE, MPI_SUM, 0, world);
+      MPI_Reduce(&com[n], &com_tmp, 1, MPI_DOUBLE, MPI_SUM, 0, world);
-      avx[n] = avx_tmp;
+      com[n] = com_tmp;
     }
 
     MPI_Reduce(&count, &count_tmp, 1, MPI_INT, MPI_SUM, 0, world);
     count = count_tmp;
 
+    com_tmp = 0.0;
+    MPI_Reduce(&totm, &com_tmp, 1, MPI_DOUBLE, MPI_SUM, 0, world);
+    totm = com_tmp;
+
     MPI_Reduce(Name, Nameall, nutypes, MPI_INT, MPI_SUM, 0, world);
     for (n = 0; n < nutypes; n++) Name[n] = Nameall[n];
 
@@ -900,16 +921,15 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec)
         }
       }
       if (count > 0) {
-        avq /= count;
         for (k = 0; k < 3; k++) {
-          avx[k] /= count;
+          com[k] /= totm;
-          if (avx[k] >= domain->boxhi[k]) avx[k] -= box[k];
+          if (com[k] >= domain->boxhi[k]) com[k] -= box[k];
-          if (avx[k] < domain->boxlo[k]) avx[k] += box[k];
+          if (com[k] < domain->boxlo[k]) com[k] += box[k];
 
-          avx[k] -= domain->boxlo[k];
+          com[k] -= domain->boxlo[k];
-          avx[k] /= box[k];
+          com[k] /= box[k];
         }
-        fprintf(pos, "\t%.8f \t%.8f \t%.8f \t%.8f", avq, avx[0], avx[1], avx[2]);
+        fprintf(pos, "\t%.8f \t%.8f \t%.8f \t%.8f", totq, com[0], com[1], com[2]);
       }
       fprintf(pos, "\n");
     }
@@ -922,21 +942,29 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec)
 
 void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec)
 {
-  int ndeletions;
+  int i, ndeletions;
   int headroom = -1;
   if (delete_Nsteps > 0) {
-    if (delete_Tcount[delete_Nsteps - 1] == -1) return;
+    if (delete_Tcount[delete_Nsteps - 1] == -1) {
+      for (i = delete_Nsteps - 1; i > 0; i--) delete_Tcount[i] = delete_Tcount[i - 1];
+      return;
+    }
     ndeletions = delete_Tcount[0] - delete_Tcount[delete_Nsteps - 1];
     if (delete_Nlimit_varid > -1)
       delete_Nlimit = input->variable->compute_equal(delete_Nlimit_varid);
     headroom = MAX(0, delete_Nlimit - ndeletions);
-    if (headroom == 0) return;
+    if (headroom == 0) {
+      for (i = delete_Nsteps - 1; i > 0; i--) delete_Tcount[i] = delete_Tcount[i - 1];
+      return;
+    }
   }
 
-  int i, j, m, n, itype, cid;
+  int j, m, n, itype, cid;
   int ndel, ndelone, count, count_tmp;
   int *Nameall;
   int *mask = atom->mask;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
   double localmass, totalmass;
   std::string species_str;
 
@@ -989,7 +1017,8 @@ void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec)
         Name[itype]++;
         count++;
         marklist[nmarklist++] = i;
-        localmass += atom->mass[atom->type[i]];
+        if (rmass) localmass += rmass[i];
+        else localmass += atom->mass[atom->type[i]];
       }
     }
 
@@ -1177,7 +1206,7 @@ double FixReaxFFSpecies::memory_usage()
 {
   double bytes;
 
-  bytes = 4 * nmax * sizeof(double);    // clusterID + x0
+  bytes = 3 * nmax * sizeof(double);    // x0
 
   return bytes;
 }

src/REAXFF/fix_reaxff_species.h

@@ -88,6 +88,7 @@ class FixReaxFFSpecies : public Fix {
 
   class NeighList *list;
   class FixAveAtom *f_SPECBOND;
+  class FixPropertyAtom *f_clusterID;
   class PairReaxFF *reaxff;
 };
 }    // namespace LAMMPS_NS

src/RHEO/fix_rheo_pressure.cpp

@@ -273,11 +273,13 @@ double FixRHEOPressure::calc_rho(double p, int i)
     error->one(FLERR,
                "Rho calculation from pressure not yet supported for cubic pressure equation");
   } else if (pressure_style[type] == TAITWATER) {
-    rho = pow(7.0 * p + csq[type] * rho0[type], SEVENTH);
+    double tmp = 7.0 * p + csq[type] * rho0[type];
+    rho = pow(MAX(0.0, tmp), SEVENTH);
     rho *= pow(rho0[type], 6.0 * SEVENTH);
     rho *= pow(csq[type], -SEVENTH);
   } else if (pressure_style[type] == TAITGENERAL) {
-    rho = pow(tpower[type] * p + csq[type] * rho0[type], 1.0 / tpower[type]);
+    double tmp = tpower[type] * p + csq[type] * rho0[type];
+    rho = pow(MAX(0.0, tmp), 1.0 / tpower[type]);
     rho *= pow(rho0[type], 1.0 - 1.0 / tpower[type]);
     rho *= pow(csq[type], -1.0 / tpower[type]);
   } else if (pressure_style[type] == IDEAL) {

src/library.cpp

@@ -2925,6 +2925,8 @@ int lammps_variable_info(void *handle, int idx, char *buffer, int buf_size) {
   return 0;
 }
 
+/* ---------------------------------------------------------------------- */
+
 /** Evaluate an immediate variable expression
  *
 \verbatim embed:rst
@@ -2958,6 +2960,86 @@ double lammps_eval(void *handle, const char *expr)
   return result;
 }
 
+/* ---------------------------------------------------------------------- */
+
+/** Clear whether a compute has been invoked.
+ *
+\verbatim embed:rst
+
+.. versionadded:: TBD
+
+   This function clears the invoked flag of all computes.
+   Called everywhere that computes are used, before computes are invoked.
+   The invoked flag is used to avoid re-invoking same compute multiple times
+   and to flag computes that store invocation times as having been invoked
+
+*See also*
+    :cpp:func:`lammps_addstep_compute_all`
+    :cpp:func:`lammps_addstep_compute`
+
+\endverbatim
+
+ * \param handle   pointer to a previously created LAMMPS instance cast to ``void *``.
+ */
+void lammps_clearstep_compute(void *handle) {
+  auto lmp = (LAMMPS *) handle;
+  lmp->modify->clearstep_compute();
+}
+
+/* ---------------------------------------------------------------------- */
+
+/** Add next timestep to all computes
+ *
+\verbatim embed:rst
+
+.. versionadded:: TBD
+
+   loop over all computes
+   schedule next invocation for those that store invocation times
+   called when not sure what computes will be needed on newstep
+   do not loop only over n_timeflag, since may not be set yet
+
+*See also*
+    :cpp:func:`lammps_clearstep_compute`
+    :cpp:func:`lammps_addstep_compute`
+
+\endverbatim
+
+ * \param handle   pointer to a previously created LAMMPS instance cast to ``void *``.
+ * \param newstep  pointer to bigint of next timestep the compute will be invoked
+ */
+void lammps_addstep_compute_all(void *handle, void *newstep) {
+  auto lmp = (LAMMPS *) handle;
+  auto ns = (bigint *) newstep;
+  if (lmp && lmp->modify && ns) lmp->modify->addstep_compute_all(*ns);
+}
+/* ---------------------------------------------------------------------- */
+
+/** Add next timestep to compute if it has been invoked in the current timestep
+ *
+\verbatim embed:rst
+
+.. versionadded:: TBD
+
+   loop over computes that store invocation times
+   if its invoked flag set on this timestep, schedule next invocation
+   called everywhere that computes are used, after computes are invoked
+
+*See also*
+    :cpp:func:`lammps_addstep_compute_all`
+    :cpp:func:`lammps_clearstep_compute`
+
+\endverbatim
+
+ * \param handle   pointer to a previously created LAMMPS instance cast to ``void *``.
+ * \param newstep  next timestep the compute will be invoked
+ */
+void lammps_addstep_compute(void *handle, void *newstep) {
+  auto lmp = (LAMMPS *) handle;
+  auto ns = (bigint *) newstep;
+  if (lmp && lmp->modify && ns) lmp->modify->addstep_compute(*ns);
+}
+
 // ----------------------------------------------------------------------
 // Library functions for scatter/gather operations of data
 // ----------------------------------------------------------------------

src/library.h

@@ -191,6 +191,10 @@ int lammps_set_internal_variable(void *handle, const char *name, double value);
 int lammps_variable_info(void *handle, int idx, char *buf, int bufsize);
 double lammps_eval(void *handle, const char *expr);
 
+void lammps_clearstep_compute(void *handle);
+void lammps_addstep_compute_all(void *handle, void * nextstep);
+void lammps_addstep_compute(void *handle, void * nextstep);
+
 /* ----------------------------------------------------------------------
  * Library functions for scatter/gather operations of data
  * ---------------------------------------------------------------------- */

tools/swig/lammps.i

@@ -143,6 +143,9 @@ extern int    lammps_set_string_variable(void *, const char *, const char *);
 extern int    lammps_set_internal_variable(void *, const char *, double);
 extern int    lammps_variable_info(void *handle, int idx, char *buf, int bufsize);
 extern double lammps_eval(void *handle, const char *expr);
+extern void   lammps_clearstep_compute(void *handle);
+extern void   lammps_addstep_compute(void *handle, void *nstep);
+extern void   lammps_addstep_compute_all(void *handle, void *nstep);
 
 extern void   lammps_gather_atoms(void *, char *, int, int, void *);
 extern void   lammps_gather_atoms_concat(void *, char *, int, int, void *);
@@ -336,6 +339,9 @@ extern int    lammps_set_string_variable(void *, const char *, const char *);
 extern int    lammps_set_internal_variable(void *, const char *, double);
 extern int    lammps_variable_info(void *handle, int idx, char *buf, int bufsize);
 extern double lammps_eval(void *handle, const char *expr);
+extern void   lammps_clearstep_compute(void *handle);
+extern void   lammps_addstep_compute(void *handle, void *nstep);
+extern void   lammps_addstep_compute_all(void *handle, void *nstep);
 
 extern void   lammps_gather_atoms(void *, char *, int, int, void *);
 extern void   lammps_gather_atoms_concat(void *, char *, int, int, void *);

unittest/c-library/test_library_properties.cpp

@@ -3,9 +3,12 @@
 #include "library.h"
 
 #include "atom.h"
+#include "compute.h"
 #include "lammps.h"
 #include "lmptype.h"
+#include "modify.h"
 #include "platform.h"
+
 #include <string>
 #include <vector>
 
@@ -668,6 +671,77 @@ TEST_F(LibraryProperties, neighlist)
     }
 };
 
+static constexpr char lj_setup[] = "lattice         fcc 0.8442\n"
+                                   "region          box block 0 10 0 10 0 10\n"
+                                   "create_box      1 box\n"
+                                   "create_atoms    1 box\n"
+                                   "mass            1 1.0\n"
+                                   "pair_style      lj/cut 2.5\n"
+                                   "pair_coeff      1 1 1.0 1.0\n"
+                                   "fix             1 all nve\n";
+
+TEST_F(LibraryProperties, step_compute)
+{
+    ::testing::internal::CaptureStdout();
+    lammps_commands_string(lmp, lj_setup);
+    lammps_command(lmp, "compute pr all pressure thermo_temp");
+    lammps_command(lmp, "fix av all ave/time 2 1 2 c_pr mode scalar");
+    lammps_command(lmp, "run 2 post no");
+    std::string output = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+    if (lammps_has_error(lmp)) {
+        char buf[2048];
+        lammps_get_last_error_message(lmp, buf, 2048);
+        FAIL() << buf << "\n";
+    }
+    auto lammps = (LAMMPS_NS::LAMMPS *)lmp;
+    auto icomp  = lammps->modify->get_compute_by_id("pr");
+    EXPECT_EQ(icomp->ntime, 2);
+    EXPECT_EQ(icomp->tlist[0], 4);
+    EXPECT_EQ(icomp->tlist[1], 2);
+    EXPECT_EQ(icomp->invoked_flag, 0);
+    EXPECT_EQ(icomp->invoked_scalar, 2);
+    EXPECT_EQ(icomp->invoked_vector, -1);
+    lammps_clearstep_compute(lmp);
+    EXPECT_EQ(icomp->invoked_flag, 0);
+    EXPECT_EQ(icomp->invoked_scalar, 2);
+    EXPECT_EQ(icomp->invoked_vector, -1);
+    bigint nextstep = 6;
+    lammps_addstep_compute(lmp, (void *)&nextstep);
+    EXPECT_EQ(icomp->ntime, 3);
+    EXPECT_EQ(icomp->tlist[0], 6);
+    EXPECT_EQ(icomp->tlist[1], 4);
+    EXPECT_EQ(icomp->tlist[2], 2);
+    EXPECT_EQ(icomp->invoked_flag, 0);
+    EXPECT_EQ(icomp->invoked_scalar, 2);
+    EXPECT_EQ(icomp->invoked_vector, -1);
+    lammps_command(lmp, "run 4 post no");
+    EXPECT_EQ(icomp->ntime, 2);
+    EXPECT_EQ(icomp->tlist[0], 8);
+    EXPECT_EQ(icomp->tlist[1], 6);
+    EXPECT_EQ(icomp->invoked_flag, 0);
+    EXPECT_EQ(icomp->invoked_scalar, 6);
+    EXPECT_EQ(icomp->invoked_vector, -1);
+    lammps_command(lmp, "run 2 post no");
+    EXPECT_EQ(icomp->ntime, 2);
+    EXPECT_EQ(icomp->tlist[0], 10);
+    EXPECT_EQ(icomp->tlist[1], 8);
+    EXPECT_EQ(icomp->invoked_flag, 0);
+    EXPECT_EQ(icomp->invoked_scalar, 8);
+    EXPECT_EQ(icomp->invoked_vector, -1);
+    nextstep = 9;
+    lammps_addstep_compute(lmp, (void *)&nextstep);
+    lammps_command(lmp, "run 1 post no");
+    EXPECT_EQ(icomp->ntime, 2);
+    EXPECT_EQ(icomp->tlist[0], 10);
+    EXPECT_EQ(icomp->tlist[1], 9);
+    EXPECT_EQ(icomp->invoked_flag, 0);
+    EXPECT_EQ(icomp->invoked_scalar, -1);
+    EXPECT_EQ(icomp->invoked_vector, -1);
+    icomp->compute_scalar();
+    EXPECT_EQ(icomp->invoked_scalar, 9);
+}
+
 TEST_F(LibraryProperties, has_error)
 {
     EXPECT_EQ(lammps_has_error(lmp), 0);

unittest/force-styles/CMakeLists.txt

@@ -246,7 +246,7 @@ if(MLIAP_ENABLE_PYTHON AND (NOT WIN32))
   add_executable(test_mliappy_unified test_mliappy_unified.cpp)
   target_link_libraries(test_mliappy_unified PRIVATE lammps GTest::GMockMain)
   add_test(NAME TestMliapPyUnified COMMAND test_mliappy_unified)
-  set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "${FORCE_TEST_ENVIRONMENT}")
+  set_tests_properties(TestMliapPyUnified PROPERTIES ENVIRONMENT "${FORCE_TEST_ENVIRONMENT}")
 endif()
 
 add_executable(test_pair_list test_pair_list.cpp)

unittest/force-styles/test_angle_style.cpp

@@ -53,15 +53,16 @@ using ::testing::StartsWith;
 
 using namespace LAMMPS_NS;
 
-void cleanup_lammps(LAMMPS *lmp, const TestConfig &cfg)
+void cleanup_lammps(LAMMPS *&lmp, const TestConfig &cfg)
 {
     platform::unlink(cfg.basename + ".restart");
     platform::unlink(cfg.basename + ".data");
     platform::unlink(cfg.basename + "-coeffs.in");
     delete lmp;
+    lmp = nullptr;
 }
 
-LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton = true)
+LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton)
 {
     LAMMPS *lmp;
 
@@ -92,21 +93,7 @@ LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton
 
     // utility lambdas to improve readability
     auto command = [&](const std::string &line) {
-        try {
         lmp->input->one(line);
-        } catch (LAMMPSAbortException &ae) {
-            fprintf(stderr, "LAMMPS Error: %s\n", ae.what());
-            exit(2);
-        } catch (LAMMPSException &e) {
-            fprintf(stderr, "LAMMPS Error: %s\n", e.what());
-            exit(3);
-        } catch (fmt::format_error &fe) {
-            fprintf(stderr, "fmt::format_error: %s\n", fe.what());
-            exit(4);
-        } catch (std::exception &e) {
-            fprintf(stderr, "General exception: %s\n", e.what());
-            exit(5);
-        }
     };
     auto parse_input_script = [&](const std::string &filename) {
         lmp->input->file(filename.c_str());
@@ -230,7 +217,12 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
     // initialize system geometry
     LAMMPS::argv args = {"AngleStyle", "-log", "none", "-echo", "screen", "-nocite"};
 
-    LAMMPS *lmp = init_lammps(args, config);
+    LAMMPS *lmp = nullptr;
+    try {
+        lmp = init_lammps(args, config, true);
+    } catch (std::exception &e) {
+        FAIL() << e.what();
+    }
     if (!lmp) {
         std::cerr << "One or more prerequisite styles are not available "
                      "in this LAMMPS configuration:\n";
@@ -321,8 +313,14 @@ TEST(AngleStyle, plain)
     LAMMPS::argv args = {"AngleStyle", "-log", "none", "-echo", "screen", "-nocite"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -371,7 +369,12 @@ TEST(AngleStyle, plain)
 
     if (!verbose) ::testing::internal::CaptureStdout();
     cleanup_lammps(lmp, test_config);
+    try {
         lmp = init_lammps(args, test_config, false);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     // skip over these tests if newton bond is forced to be on
@@ -439,8 +442,14 @@ TEST(AngleStyle, omp)
                          "-pk",        "omp",  "4",    "-sf",   "omp"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -493,7 +502,12 @@ TEST(AngleStyle, omp)
 
     if (!verbose) ::testing::internal::CaptureStdout();
     cleanup_lammps(lmp, test_config);
+    try {
         lmp = init_lammps(args, test_config, false);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     // skip over these tests if newton bond is forced to be on
@@ -540,14 +554,21 @@ TEST(AngleStyle, kokkos_omp)
     // if KOKKOS has GPU support enabled, it *must* be used. We cannot test OpenMP only.
     if (Info::has_accelerator_feature("KOKKOS", "api", "cuda") ||
         Info::has_accelerator_feature("KOKKOS", "api", "hip") ||
-        Info::has_accelerator_feature("KOKKOS", "api", "sycl")) GTEST_SKIP();
+        Info::has_accelerator_feature("KOKKOS", "api", "sycl"))
+        GTEST_SKIP() << "Cannot test KOKKOS/OpenMP with GPU support enabled";
 
     LAMMPS::argv args = {"AngleStyle", "-log", "none", "-echo", "screen", "-nocite",
                          "-k",         "on",   "t",    "4",     "-sf",    "kk"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -597,7 +618,12 @@ TEST(AngleStyle, kokkos_omp)
 
     if (!verbose) ::testing::internal::CaptureStdout();
     cleanup_lammps(lmp, test_config);
+    try {
         lmp = init_lammps(args, test_config, false);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     // skip over these tests if newton bond is forced to be on
@@ -664,8 +690,14 @@ TEST(AngleStyle, numdiff)
     LAMMPS::argv args = {"AngleStyle", "-log", "none", "-echo", "screen", "-nocite"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -717,7 +749,13 @@ TEST(AngleStyle, single)
 
     // create a LAMMPS instance with standard settings to detect the number of atom types
     if (!verbose) ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config);
+    LAMMPS *lmp = nullptr;
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     if (!lmp) {
@@ -865,7 +903,13 @@ TEST(AngleStyle, extract)
     LAMMPS::argv args = {"AngleStyle", "-log", "none", "-echo", "screen", "-nocite"};
 
     if (!verbose) ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     if (!lmp) {

unittest/force-styles/test_bond_style.cpp

@@ -53,15 +53,16 @@ using ::testing::StartsWith;
 
 using namespace LAMMPS_NS;
 
-void cleanup_lammps(LAMMPS *lmp, const TestConfig &cfg)
+void cleanup_lammps(LAMMPS *&lmp, const TestConfig &cfg)
 {
     platform::unlink(cfg.basename + ".restart");
     platform::unlink(cfg.basename + ".data");
     platform::unlink(cfg.basename + "-coeffs.in");
     delete lmp;
+    lmp = nullptr;
 }
 
-LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton = true)
+LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton)
 {
     LAMMPS *lmp;
 
@@ -92,21 +93,7 @@ LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton
 
     // utility lambdas to improve readability
     auto command = [&](const std::string &line) {
-        try {
         lmp->input->one(line);
-        } catch (LAMMPSAbortException &ae) {
-            fprintf(stderr, "LAMMPS Error: %s\n", ae.what());
-            exit(2);
-        } catch (LAMMPSException &e) {
-            fprintf(stderr, "LAMMPS Error: %s\n", e.what());
-            exit(3);
-        } catch (fmt::format_error &fe) {
-            fprintf(stderr, "fmt::format_error: %s\n", fe.what());
-            exit(4);
-        } catch (std::exception &e) {
-            fprintf(stderr, "General exception: %s\n", e.what());
-            exit(5);
-        }
     };
     auto parse_input_script = [&](const std::string &filename) {
         lmp->input->file(filename.c_str());
@@ -230,7 +217,13 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
     // initialize system geometry
     LAMMPS::argv args = {"BondStyle", "-log", "none", "-echo", "screen", "-nocite"};
 
-    LAMMPS *lmp = init_lammps(args, config);
+    LAMMPS *lmp = nullptr;
+    try {
+        lmp = init_lammps(args, config, true);
+    } catch (std::exception &e) {
+        FAIL() << e.what();
+    }
+
     if (!lmp) {
         std::cerr << "One or more prerequisite styles are not available "
                      "in this LAMMPS configuration:\n";
@@ -321,8 +314,14 @@ TEST(BondStyle, plain)
     LAMMPS::argv args = {"BondStyle", "-log", "none", "-echo", "screen", "-nocite"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -371,7 +370,12 @@ TEST(BondStyle, plain)
 
     if (!verbose) ::testing::internal::CaptureStdout();
     cleanup_lammps(lmp, test_config);
+    try {
         lmp = init_lammps(args, test_config, false);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     // skip over these tests if newton bond is forced to be on
@@ -441,8 +445,14 @@ TEST(BondStyle, omp)
                          "-pk",       "omp",  "4",    "-sf",   "omp"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -495,7 +505,12 @@ TEST(BondStyle, omp)
 
     if (!verbose) ::testing::internal::CaptureStdout();
     cleanup_lammps(lmp, test_config);
+    try {
         lmp = init_lammps(args, test_config, false);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     // skip over these tests if newton bond is forced to be on
@@ -542,14 +557,21 @@ TEST(BondStyle, kokkos_omp)
     // if KOKKOS has GPU support enabled, it *must* be used. We cannot test OpenMP only.
     if (Info::has_accelerator_feature("KOKKOS", "api", "cuda") ||
         Info::has_accelerator_feature("KOKKOS", "api", "hip") ||
-        Info::has_accelerator_feature("KOKKOS", "api", "sycl")) GTEST_SKIP();
+        Info::has_accelerator_feature("KOKKOS", "api", "sycl"))
+        GTEST_SKIP() << "Cannot test KOKKOS/OpenMP with GPU support enabled";
 
     LAMMPS::argv args = {"BondStyle", "-log", "none", "-echo", "screen", "-nocite",
                          "-k",        "on",   "t",    "4",     "-sf",    "kk"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -603,7 +625,12 @@ TEST(BondStyle, kokkos_omp)
 
     if (!verbose) ::testing::internal::CaptureStdout();
     cleanup_lammps(lmp, test_config);
+    try {
         lmp = init_lammps(args, test_config, false);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     // skip over these tests if newton bond is forced to be on
@@ -652,8 +679,14 @@ TEST(BondStyle, numdiff)
     LAMMPS::argv args = {"BondStyle", "-log", "none", "-echo", "screen", "-nocite"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -705,7 +738,13 @@ TEST(BondStyle, single)
 
     // create a LAMMPS instance with standard settings to detect the number of atom types
     if (!verbose) ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config);
+    LAMMPS *lmp = nullptr;
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     if (!lmp) {
@@ -959,7 +998,13 @@ TEST(BondStyle, extract)
     LAMMPS::argv args = {"BondStyle", "-log", "none", "-echo", "screen", "-nocite"};
 
     if (!verbose) ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     if (!lmp) {

unittest/force-styles/test_dihedral_style.cpp

@@ -53,15 +53,16 @@ using ::testing::StartsWith;
 
 using namespace LAMMPS_NS;
 
-void cleanup_lammps(LAMMPS *lmp, const TestConfig &cfg)
+void cleanup_lammps(LAMMPS *&lmp, const TestConfig &cfg)
 {
     platform::unlink(cfg.basename + ".restart");
     platform::unlink(cfg.basename + ".data");
     platform::unlink(cfg.basename + "-coeffs.in");
     delete lmp;
+    lmp = nullptr;
 }
 
-LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton = true)
+LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton)
 {
     auto *lmp = new LAMMPS(args, MPI_COMM_WORLD);
 
@@ -237,7 +238,13 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
     // initialize system geometry
     LAMMPS::argv args = {"DihedralStyle", "-log", "none", "-echo", "screen", "-nocite"};
 
-    LAMMPS *lmp = init_lammps(args, config);
+    LAMMPS *lmp = nullptr;
+    try {
+        lmp = init_lammps(args, config, true);
+    } catch (std::exception &e) {
+        FAIL() << e.what();
+    }
+
     if (!lmp) {
         std::cerr << "One or more prerequisite styles are not available "
                      "in this LAMMPS configuration:\n";
@@ -322,8 +329,14 @@ TEST(DihedralStyle, plain)
     LAMMPS::argv args = {"DihedralStyle", "-log", "none", "-echo", "screen", "-nocite"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -372,7 +385,12 @@ TEST(DihedralStyle, plain)
 
     if (!verbose) ::testing::internal::CaptureStdout();
     cleanup_lammps(lmp, test_config);
+    try {
         lmp = init_lammps(args, test_config, false);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     // skip over these tests if newton bond is forced to be on
@@ -442,8 +460,14 @@ TEST(DihedralStyle, omp)
                          "-pk",           "omp",  "4",    "-sf",   "omp"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -497,7 +521,12 @@ TEST(DihedralStyle, omp)
 
     if (!verbose) ::testing::internal::CaptureStdout();
     cleanup_lammps(lmp, test_config);
+    try {
         lmp = init_lammps(args, test_config, false);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     // skip over these tests if newton bond is forced to be on
@@ -544,15 +573,22 @@ TEST(DihedralStyle, kokkos_omp)
     // if KOKKOS has GPU support enabled, it *must* be used. We cannot test OpenMP only.
     if (Info::has_accelerator_feature("KOKKOS", "api", "cuda") ||
         Info::has_accelerator_feature("KOKKOS", "api", "hip") ||
-        Info::has_accelerator_feature("KOKKOS", "api", "sycl")) GTEST_SKIP();
+        Info::has_accelerator_feature("KOKKOS", "api", "sycl"))
+        GTEST_SKIP() << "Cannot test KOKKOS/OpenMP with GPU support enabled";
 
     LAMMPS::argv args = {"DihedralStyle", "-log", "none", "-echo", "screen",
                          "-nocite",       "-k",   "on",   "t",     "4",
                          "-sf",           "kk"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -606,7 +642,12 @@ TEST(DihedralStyle, kokkos_omp)
 
     if (!verbose) ::testing::internal::CaptureStdout();
     cleanup_lammps(lmp, test_config);
+    try {
         lmp = init_lammps(args, test_config, false);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     // skip over these tests if newton bond is forced to be on
@@ -655,8 +696,14 @@ TEST(DihedralStyle, numdiff)
     LAMMPS::argv args = {"DihedralStyle", "-log", "none", "-echo", "screen", "-nocite"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 

unittest/force-styles/test_fix_timestep.cpp

@@ -56,13 +56,14 @@ using ::testing::StartsWith;
 
 using namespace LAMMPS_NS;
 
-void cleanup_lammps(LAMMPS *lmp, const TestConfig &cfg)
+void cleanup_lammps(LAMMPS *&lmp, const TestConfig &cfg)
 {
     platform::unlink(cfg.basename + ".restart");
     delete lmp;
+    lmp = nullptr;
 }
 
-LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool use_respa = false)
+LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool use_respa)
 {
     LAMMPS *lmp;
 
@@ -178,7 +179,12 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
 {
     // initialize system geometry
     LAMMPS::argv args = {"FixIntegrate", "-log", "none", "-echo", "screen", "-nocite"};
-    LAMMPS *lmp       = init_lammps(args, config);
+    LAMMPS *lmp       = nullptr;
+    try {
+        lmp = init_lammps(args, config, false);
+    } catch (std::exception &e) {
+        FAIL() << e.what();
+    }
     if (!lmp) {
         std::cerr << "One or more prerequisite styles are not available "
                      "in this LAMMPS configuration:\n";
@@ -272,7 +278,14 @@ TEST(FixTimestep, plain)
     LAMMPS::argv args = {"FixTimestep", "-log", "none", "-echo", "screen", "-nocite"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp        = init_lammps(args, test_config);
+    LAMMPS *lmp = nullptr;
+    try {
+        lmp = init_lammps(args, test_config, false);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -430,14 +443,20 @@ TEST(FixTimestep, plain)
     // fix nve/limit cannot work with r-RESPA
     ifix = lmp->modify->get_fix_by_id("test");
     if (ifix && !utils::strmatch(ifix->style, "^rigid") &&
-        !utils::strmatch(ifix->style, "^nve/limit") &&
+        !utils::strmatch(ifix->style, "^nve/limit") && !utils::strmatch(ifix->style, "^recenter")) {
-        !utils::strmatch(ifix->style, "^recenter")) {
         if (!verbose) ::testing::internal::CaptureStdout();
         cleanup_lammps(lmp, test_config);
         if (!verbose) ::testing::internal::GetCapturedStdout();
 
         ::testing::internal::CaptureStdout();
+        LAMMPS *lmp = nullptr;
+        try {
             lmp = init_lammps(args, test_config, true);
+        } catch (std::exception &e) {
+            output = ::testing::internal::GetCapturedStdout();
+            if (verbose) std::cout << output;
+            FAIL() << e.what();
+        }
         output = ::testing::internal::GetCapturedStdout();
         if (verbose) std::cout << output;
 
@@ -572,15 +591,26 @@ TEST(FixTimestep, omp)
                          "-pk",         "omp",  "4",    "-sf",   "omp"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp        = init_lammps(args, test_config);
+    LAMMPS *lmp = nullptr;
+    try {
+        lmp = init_lammps(args, test_config, false);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
     if (!lmp) {
-        std::cerr << "One or more prerequisite styles are not available "
+        std::cerr << "One or more prerequisite styles with /omp suffix are not available "
                      "in this LAMMPS configuration:\n";
         for (auto &prerequisite : test_config.prerequisites) {
+            if (prerequisite.first == "atom") {
                 std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
+            } else {
+                std::cerr << prerequisite.first << "_style " << prerequisite.second << "/omp\n";
+            }
         }
         GTEST_SKIP();
     }
@@ -731,7 +761,13 @@ TEST(FixTimestep, omp)
         if (!verbose) ::testing::internal::GetCapturedStdout();
 
         ::testing::internal::CaptureStdout();
+        try {
             lmp = init_lammps(args, test_config, true);
+        } catch (std::exception &e) {
+            output = ::testing::internal::GetCapturedStdout();
+            if (verbose) std::cout << output;
+            FAIL() << e.what();
+        }
         output = ::testing::internal::GetCapturedStdout();
         if (verbose) std::cout << output;
 
@@ -861,13 +897,21 @@ TEST(FixTimestep, kokkos_omp)
     // if KOKKOS has GPU support enabled, it *must* be used. We cannot test OpenMP only.
     if (Info::has_accelerator_feature("KOKKOS", "api", "cuda") ||
         Info::has_accelerator_feature("KOKKOS", "api", "hip") ||
-        Info::has_accelerator_feature("KOKKOS", "api", "sycl")) GTEST_SKIP();
+        Info::has_accelerator_feature("KOKKOS", "api", "sycl")) {
-
+        GTEST_SKIP() << "Cannot test KOKKOS/OpenMP with GPU support enabled";
+    }
     LAMMPS::argv args = {"FixTimestep", "-log", "none", "-echo", "screen", "-nocite",
                          "-k",          "on",   "t",    "4",     "-sf",    "kk"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp        = init_lammps(args, test_config);
+    LAMMPS *lmp = nullptr;
+    try {
+        lmp = init_lammps(args, test_config, false);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -875,7 +919,7 @@ TEST(FixTimestep, kokkos_omp)
         std::cerr << "One or more prerequisite styles with /kk suffix\n"
                      "are not available in this LAMMPS configuration:\n";
         for (auto &prerequisite : test_config.prerequisites) {
-            std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
+            std::cerr << prerequisite.first << "_style " << prerequisite.second << "/kk\n";
         }
         GTEST_SKIP();
     }

unittest/force-styles/test_improper_style.cpp

@@ -53,15 +53,16 @@ using ::testing::StartsWith;
 
 using namespace LAMMPS_NS;
 
-void cleanup_lammps(LAMMPS *lmp, const TestConfig &cfg)
+void cleanup_lammps(LAMMPS *&lmp, const TestConfig &cfg)
 {
     platform::unlink(cfg.basename + ".restart");
     platform::unlink(cfg.basename + ".data");
     platform::unlink(cfg.basename + "-coeffs.in");
     delete lmp;
+    lmp = nullptr;
 }
 
-LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton = true)
+LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton)
 {
     LAMMPS *lmp;
 
@@ -92,21 +93,7 @@ LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton
 
     // utility lambdas to improve readability
     auto command = [&](const std::string &line) {
-        try {
         lmp->input->one(line);
-        } catch (LAMMPSAbortException &ae) {
-            fprintf(stderr, "LAMMPS Error: %s\n", ae.what());
-            exit(2);
-        } catch (LAMMPSException &e) {
-            fprintf(stderr, "LAMMPS Error: %s\n", e.what());
-            exit(3);
-        } catch (fmt::format_error &fe) {
-            fprintf(stderr, "fmt::format_error: %s\n", fe.what());
-            exit(4);
-        } catch (std::exception &e) {
-            fprintf(stderr, "General exception: %s\n", e.what());
-            exit(5);
-        }
     };
     auto parse_input_script = [&](const std::string &filename) {
         lmp->input->file(filename.c_str());
@@ -230,7 +217,12 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
     // initialize system geometry
     LAMMPS::argv args = {"ImproperStyle", "-log", "none", "-echo", "screen", "-nocite"};
 
-    LAMMPS *lmp = init_lammps(args, config);
+    LAMMPS *lmp = nullptr;
+    try {
+        lmp = init_lammps(args, config, true);
+    } catch (std::exception &e) {
+        FAIL() << e.what();
+    }
     if (!lmp) {
         std::cerr << "One or more prerequisite styles are not available "
                      "in this LAMMPS configuration:\n";
@@ -315,8 +307,14 @@ TEST(ImproperStyle, plain)
     LAMMPS::argv args = {"ImproperStyle", "-log", "none", "-echo", "screen", "-nocite"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -365,7 +363,12 @@ TEST(ImproperStyle, plain)
 
     if (!verbose) ::testing::internal::CaptureStdout();
     cleanup_lammps(lmp, test_config);
+    try {
         lmp = init_lammps(args, test_config, false);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     // skip over these tests if newton bond is forced to be on
@@ -435,8 +438,14 @@ TEST(ImproperStyle, omp)
                          "-pk",           "omp",  "4",    "-sf",   "omp"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -490,7 +499,12 @@ TEST(ImproperStyle, omp)
 
     if (!verbose) ::testing::internal::CaptureStdout();
     cleanup_lammps(lmp, test_config);
+    try {
         lmp = init_lammps(args, test_config, false);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     // skip over these tests if newton bond is forced to be on
@@ -536,15 +550,22 @@ TEST(ImproperStyle, kokkos_omp)
     // if KOKKOS has GPU support enabled, it *must* be used. We cannot test OpenMP only.
     if (Info::has_accelerator_feature("KOKKOS", "api", "cuda") ||
         Info::has_accelerator_feature("KOKKOS", "api", "hip") ||
-        Info::has_accelerator_feature("KOKKOS", "api", "sycl")) GTEST_SKIP();
+        Info::has_accelerator_feature("KOKKOS", "api", "sycl"))
+        GTEST_SKIP() << "Cannot test KOKKOS/OpenMP with GPU support enabled";
 
     LAMMPS::argv args = {"ImproperStyle", "-log", "none", "-echo", "screen",
                          "-nocite",       "-k",   "on",   "t",     "4",
                          "-sf",           "kk"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -597,7 +618,12 @@ TEST(ImproperStyle, kokkos_omp)
 
     if (!verbose) ::testing::internal::CaptureStdout();
     cleanup_lammps(lmp, test_config);
+    try {
         lmp = init_lammps(args, test_config, false);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     // skip over these tests if newton bond is forced to be on
@@ -644,8 +670,14 @@ TEST(ImproperStyle, numdiff)
     LAMMPS::argv args = {"ImproperStyle", "-log", "none", "-echo", "screen", "-nocite"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 

unittest/force-styles/test_pair_style.cpp

@@ -54,15 +54,16 @@ using ::testing::StartsWith;
 
 using namespace LAMMPS_NS;
 
-void cleanup_lammps(LAMMPS *lmp, const TestConfig &cfg)
+void cleanup_lammps(LAMMPS *&lmp, const TestConfig &cfg)
 {
     platform::unlink(cfg.basename + ".restart");
     platform::unlink(cfg.basename + ".data");
     platform::unlink(cfg.basename + "-coeffs.in");
     delete lmp;
+    lmp = nullptr;
 }
 
-LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton = true)
+LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton)
 {
     LAMMPS *lmp;
 
@@ -93,21 +94,7 @@ LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton
 
     // utility lambdas to improve readability
     auto command = [&](const std::string &line) {
-        try {
         lmp->input->one(line);
-        } catch (LAMMPSAbortException &ae) {
-            fprintf(stderr, "LAMMPS Error: %s\n", ae.what());
-            exit(2);
-        } catch (LAMMPSException &e) {
-            fprintf(stderr, "LAMMPS Error: %s\n", e.what());
-            exit(3);
-        } catch (fmt::format_error &fe) {
-            fprintf(stderr, "fmt::format_error: %s\n", fe.what());
-            exit(4);
-        } catch (std::exception &e) {
-            fprintf(stderr, "General exception: %s\n", e.what());
-            exit(5);
-        }
     };
 
     auto parse_input_script = [&](const std::string &filename) {
@@ -242,8 +229,12 @@ void generate_yaml_file(const char *outfile, const TestConfig &config)
 {
     // initialize system geometry
     LAMMPS::argv args = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite"};
-
+    LAMMPS *lmp       = nullptr;
-    LAMMPS *lmp = init_lammps(args, config);
+    try {
+        lmp = init_lammps(args, config, true);
+    } catch (std::exception &e) {
+        FAIL() << e.what();
+    }
     if (!lmp) {
         std::cerr << "One or more prerequisite styles are not available "
                      "in this LAMMPS configuration:\n";
@@ -340,8 +331,14 @@ TEST(PairStyle, plain)
     LAMMPS::argv args = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -399,7 +396,12 @@ TEST(PairStyle, plain)
 
     if (!verbose) ::testing::internal::CaptureStdout();
     cleanup_lammps(lmp, test_config);
+    try {
         lmp = init_lammps(args, test_config, false);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     // skip over these tests if newton pair is forced to be on
@@ -480,9 +482,14 @@ TEST(PairStyle, plain)
     if (pair->respa_enable) {
         if (!verbose) ::testing::internal::CaptureStdout();
         cleanup_lammps(lmp, test_config);
+        try {
             lmp = init_lammps(args, test_config, false);
             lmp->input->one("run_style respa 2 1 inner 1 4.8 5.5 outer 2");
             run_lammps(lmp);
+        } catch (std::exception &e) {
+            if (!verbose) ::testing::internal::GetCapturedStdout();
+            FAIL() << e.what();
+        }
         if (!verbose) ::testing::internal::GetCapturedStdout();
 
         // need to relax error by a large amount with tabulation, since
@@ -519,8 +526,14 @@ TEST(PairStyle, omp)
     if (utils::strmatch(test_config.pair_style, "^dpd")) args[8] = "1";
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -578,7 +591,12 @@ TEST(PairStyle, omp)
 
         if (!verbose) ::testing::internal::CaptureStdout();
         cleanup_lammps(lmp, test_config);
+        try {
             lmp = init_lammps(args, test_config, false);
+        } catch (std::exception &e) {
+            if (!verbose) ::testing::internal::GetCapturedStdout();
+            FAIL() << e.what();
+        }
         if (!verbose) ::testing::internal::GetCapturedStdout();
 
         pair = lmp->force->pair;
@@ -636,7 +654,9 @@ TEST(PairStyle, kokkos_omp)
     // if KOKKOS has GPU support enabled, it *must* be used. We cannot test OpenMP only.
     if (Info::has_accelerator_feature("KOKKOS", "api", "cuda") ||
         Info::has_accelerator_feature("KOKKOS", "api", "hip") ||
-        Info::has_accelerator_feature("KOKKOS", "api", "sycl")) GTEST_SKIP();
+        Info::has_accelerator_feature("KOKKOS", "api", "sycl")) {
+        GTEST_SKIP() << "Cannot test KOKKOS/OpenMP with GPU support enabled";
+    }
 
     LAMMPS::argv args = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite",
                          "-k",        "on",   "t",    "4",     "-sf",    "kk"};
@@ -655,8 +675,14 @@ TEST(PairStyle, kokkos_omp)
         args[9] = "1";
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -713,7 +739,12 @@ TEST(PairStyle, kokkos_omp)
     if (lmp->force->newton_pair == 0) {
         if (!verbose) ::testing::internal::CaptureStdout();
         cleanup_lammps(lmp, test_config);
+        try {
             lmp = init_lammps(args, test_config, false);
+        } catch (std::exception &e) {
+            if (!verbose) ::testing::internal::GetCapturedStdout();
+            FAIL() << e.what();
+        }
         if (!verbose) ::testing::internal::GetCapturedStdout();
 
         pair = lmp->force->pair;
@@ -788,8 +819,14 @@ TEST(PairStyle, gpu)
     }
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, false);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, false);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -868,8 +905,14 @@ TEST(PairStyle, intel)
     if (utils::strmatch(test_config.pair_style, "^dpd")) args[12] = "1";
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -948,8 +991,14 @@ TEST(PairStyle, opt)
     LAMMPS::argv args = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite", "-sf", "opt"};
 
     ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config);
+    LAMMPS *lmp = nullptr;
-
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        std::string output = ::testing::internal::GetCapturedStdout();
+        if (verbose) std::cout << output;
+        FAIL() << e.what();
+    }
     std::string output = ::testing::internal::GetCapturedStdout();
     if (verbose) std::cout << output;
 
@@ -1032,7 +1081,13 @@ TEST(PairStyle, single)
 
     // create a LAMMPS instance with standard settings to detect the number of atom types
     if (!verbose) ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config);
+    LAMMPS *lmp = nullptr;
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     if (!lmp) {
@@ -1276,7 +1331,13 @@ TEST(PairStyle, extract)
     LAMMPS::argv args = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite"};
 
     if (!verbose) ::testing::internal::CaptureStdout();
-    LAMMPS *lmp = init_lammps(args, test_config, true);
+    LAMMPS *lmp = nullptr;
+    try {
+        lmp = init_lammps(args, test_config, true);
+    } catch (std::exception &e) {
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        FAIL() << e.what();
+    }
     if (!verbose) ::testing::internal::GetCapturedStdout();
 
     if (!lmp) {

unittest/force-styles/tests/fix-timestep-efield_lepton.yaml

@@ -6,7 +6,7 @@ epsilon: 2e-13
 skip_tests:
 prerequisites: ! |
   atom full
-  fix efield
+  fix efield/lepton
 pre_commands: ! ""
 post_commands: ! |
   region half block 0 EDGE EDGE EDGE EDGE EDGE

unittest/force-styles/tests/fix-timestep-efield_lepton_region.yaml

@@ -6,7 +6,7 @@ epsilon: 2e-13
 skip_tests:
 prerequisites: ! |
   atom full
-  fix efield
+  fix efield/lepton
 pre_commands: ! ""
 post_commands: ! |
   region half block 0 EDGE EDGE EDGE EDGE EDGE