diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS
index c3b3ddbce19..331c9569ff5 100644
--- a/.github/CODEOWNERS
+++ b/.github/CODEOWNERS
@@ -174,6 +174,10 @@
 
 # docs
 /doc/*                  @akohlmey
+/doc/src/Howto_viz.rst  @akohlmey
+/doc/src/Howto_spc.rst  @akohlmey
+/doc/src/Howto_tip*.rst @akohlmey
+/examples/GRAPHICS/     @akohlmey
 /examples/plugin/       @akohlmey
 /examples/PACKAGES/pace/plugin/  @akohlmey
 
diff --git a/.github/release_steps.md b/.github/release_steps.md
index 6a3bda821cf..76863f74e87 100644
--- a/.github/release_steps.md
+++ b/.github/release_steps.md
@@ -212,6 +212,26 @@ gh release upload patch_4Feb2025 LAMMPS-Win10-64bit-GUI-4Feb2025.exe
 The symbolic link is used to have a consistent naming scheme for the packages
 attached to the GitHub release page.
 
+#### LAMMPS Online Manual
+
+Creating the online docs requires setting some environment variables to have
+the extra box at the bottom of the navigation bar included that allows to
+switch between release, stable, and develop branch versions of the manual.
+Also the search box should use Google search instead of the javascript search.
+
+``` sh
+cd lammps-release
+make -C doc clean
+make -C doc upgrade
+env LAMMPS_WEBSITE_BUILD=1 LAMMPS_WEBSITE_BUILD_VERSION=release LAMMPS_WEBSITE_BASEURL="https://docs.lammps.org/" \
+    make -C doc html WEB_SEARCH=YES
+make -C doc pdf
+mv doc/Manual.pdf doc/html
+rsync -arpv doc/html/ www.lammps.org:/var/www/lammps/docs/release-new
+
+Then log into www.lammps.org and move the current folder away and
+the new folder in its place and update permissions.
+
 #### Clean up:
 
 ``` sh
@@ -243,6 +263,7 @@ and installed (e.g. to `$HOME/.local`) as static libraries only:
 - jpeg
 - zlib
 - png
+- voro++ (for VORONOI package)
 - Qt (for LAMMPS-GUI)
 
 When configuring LAMMPS the `cmake/presets/clang.cmake` should be used
@@ -298,8 +319,8 @@ Check out the LAMMPS website repo
 https://github.com/lammps/lammps-website.git and edit the file
 `src/download.txt` for the new release.  Test translation with `make
 html` and review `html/download.html` Then add and commit to git and
-push the changes to GitHub.  The Temple Jenkis cluster will
-automatically update https://www.lammps.org/download.html accordingly.
+push the changes to GitHub.  A cron job will automatically update
+https://www.lammps.org/ accordingly if there are changes.
 
 Also notify Steve of the release so he can update `src/bug.txt` on the
 website from the available release notes.
@@ -373,11 +394,10 @@ git tag -s -m 'Update 2 for Stable LAMMPS version 29 August 2024' stable_29Aug20
 git push git@github.com:lammps/lammps.git --tags maintenance stable
 ```
 
-Associate draft release notes with new tag and publish as "latest release".
-
-On https://ci.lammps.org/ go to "dev", "stable" and manually execute
-the "update\_release" task. This will update https://docs.lammps.org/stable
-and prepare a stable tarball.
+Associate draft release notes with new tag and prepare and upload
+various source and binary packages. Publish as "latest release"
+only when all uploads are complete as they cannot be changed after
+they are published since we are creating "immutable" releases now.
 
 ### Build and upload binary packages and source tarball to GitHub
 
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
index d28aae24256..a9a4bd861f8 100644
--- a/cmake/CMakeLists.txt
+++ b/cmake/CMakeLists.txt
@@ -193,7 +193,10 @@ endif()
 
 # do not include the (obsolete) MPI C++ bindings which makes for leaner object files
 # and avoids namespace conflicts. Put this early to increase its visbility.
-set(MPI_CXX_SKIP_MPICXX TRUE CACHE BOOL "Skip MPI C++ Bindings" FORCE)
+# do not apply when compiling SCAFACOS since it uses the obsolete bindings
+if (NOT PKG_SCAFACOS)
+  set(MPI_CXX_SKIP_MPICXX TRUE CACHE BOOL "Skip MPI C++ Bindings" FORCE)
+endif()
 
 ########################################################################
 # User input options                                                   #
diff --git a/cmake/Modules/Packages/SCAFACOS.cmake b/cmake/Modules/Packages/SCAFACOS.cmake
index 2905a207b0b..b0455c5f006 100644
--- a/cmake/Modules/Packages/SCAFACOS.cmake
+++ b/cmake/Modules/Packages/SCAFACOS.cmake
@@ -3,7 +3,7 @@ enable_language(C)
 
 find_package(GSL REQUIRED)
 find_package(PkgConfig QUIET)
-find_package(MPI REQUIRED)
+find_package(MPI REQUIRED COMPONENTS C MPICXX CXX Fortran)
 set(DOWNLOAD_SCAFACOS_DEFAULT ON)
 if(PKG_CONFIG_FOUND)
   pkg_check_modules(SCAFACOS QUIET scafacos)
@@ -19,6 +19,8 @@ if(DOWNLOAD_SCAFACOS)
   mark_as_advanced(SCAFACOS_URL)
   mark_as_advanced(SCAFACOS_MD5)
   GetFallbackURL(SCAFACOS_URL SCAFACOS_FALLBACK)
+  set(SCAFACOS_CXX_FLAGS "${CMAKE_CXX_FLAGS_${CMAKE_BUILD_TYPE}} ${CMAKE_CXX_FLAGS}")
+  set(SCAFACOS_C_FLAGS "${CMAKE_C_FLAGS_${CMAKE_BUILD_TYPE}} ${CMAKE_C_FLAGS}")
 
   include(ExternalProject)
   ExternalProject_Add(scafacos_build
@@ -32,6 +34,8 @@ if(DOWNLOAD_SCAFACOS)
                                              CXX=${CMAKE_MPI_CXX_COMPILER}
                                              CC=${CMAKE_MPI_C_COMPILER}
                                              F77=
+                                             CFLAGS=${SCAFACOS_C_FLAGS}
+                                             CXXFLAGS=${SCAFACOS_CXX_FLAGS}
     BUILD_BYPRODUCTS
       <INSTALL_DIR>/lib/libfcs.a
       <INSTALL_DIR>/lib/libfcs_direct.a
@@ -55,7 +59,7 @@ if(DOWNLOAD_SCAFACOS)
   set_target_properties(LAMMPS::SCAFACOS PROPERTIES
     IMPORTED_LOCATION "${INSTALL_DIR}/lib/libfcs.a"
     INTERFACE_INCLUDE_DIRECTORIES "${INSTALL_DIR}/include"
-    INTERFACE_LINK_LIBRARIES "${INSTALL_DIR}/lib/libfcs.a;${INSTALL_DIR}/lib/libfcs_direct.a;${INSTALL_DIR}/lib/libfcs_ewald.a;${INSTALL_DIR}/lib/libfcs_fmm.a;${INSTALL_DIR}/lib/libfcs_p2nfft.a;${INSTALL_DIR}/lib/libfcs_p3m.a;GSL::gsl;${INSTALL_DIR}/lib/libfcs_near.a;${INSTALL_DIR}/lib/libfcs_gridsort.a;${INSTALL_DIR}/lib/libfcs_resort.a;${INSTALL_DIR}/lib/libfcs_redist.a;${INSTALL_DIR}/lib/libfcs_common.a;${INSTALL_DIR}/lib/libfcs_pnfft.a;${INSTALL_DIR}/lib/libfcs_pfft.a;${INSTALL_DIR}/lib/libfcs_fftw3_mpi.a;${INSTALL_DIR}/lib/libfcs_fftw3.a;MPI::MPI_Fortran;MPI::MPI_C")
+    INTERFACE_LINK_LIBRARIES "${INSTALL_DIR}/lib/libfcs.a;${INSTALL_DIR}/lib/libfcs_direct.a;${INSTALL_DIR}/lib/libfcs_ewald.a;${INSTALL_DIR}/lib/libfcs_fmm.a;${INSTALL_DIR}/lib/libfcs_p2nfft.a;${INSTALL_DIR}/lib/libfcs_p3m.a;GSL::gsl;${INSTALL_DIR}/lib/libfcs_near.a;${INSTALL_DIR}/lib/libfcs_gridsort.a;${INSTALL_DIR}/lib/libfcs_resort.a;${INSTALL_DIR}/lib/libfcs_redist.a;${INSTALL_DIR}/lib/libfcs_common.a;${INSTALL_DIR}/lib/libfcs_pnfft.a;${INSTALL_DIR}/lib/libfcs_pfft.a;${INSTALL_DIR}/lib/libfcs_fftw3_mpi.a;${INSTALL_DIR}/lib/libfcs_fftw3.a;MPI::MPI_CXX;MPI::MPI_Fortran;MPI::MPI_C;${CMAKE_Fortran_IMPLICIT_LINK_LIBRARIES}")
   target_link_libraries(lammps PRIVATE LAMMPS::SCAFACOS)
   add_dependencies(LAMMPS::SCAFACOS scafacos_build)
 else()
diff --git a/doc/doxygen/Doxyfile.in b/doc/doxygen/Doxyfile.in
index 1439f5b54f6..48f456f6d6a 100644
--- a/doc/doxygen/Doxyfile.in
+++ b/doc/doxygen/Doxyfile.in
@@ -431,6 +431,8 @@ INPUT                  = @LAMMPS_SOURCE_DIR@/utils.cpp                 \
                          @LAMMPS_SOURCE_DIR@/text_file_reader.h        \
                          @LAMMPS_SOURCE_DIR@/potential_file_reader.cpp \
                          @LAMMPS_SOURCE_DIR@/potential_file_reader.h   \
+                         @LAMMPS_SOURCE_DIR@/label_map.cpp             \
+                         @LAMMPS_SOURCE_DIR@/label_map.h               \
                          @LAMMPS_SOURCE_DIR@/my_page.cpp               \
                          @LAMMPS_SOURCE_DIR@/my_page.h                 \
                          @LAMMPS_SOURCE_DIR@/my_pool_chunk.cpp         \
diff --git a/doc/src/Developer_utils.rst b/doc/src/Developer_utils.rst
index 866945fc881..cd9073f43e8 100644
--- a/doc/src/Developer_utils.rst
+++ b/doc/src/Developer_utils.rst
@@ -558,6 +558,142 @@ A file that would be parsed by the reader code fragment looks like this:
 
 ----------
 
+Type label support
+------------------
+
+Overview
+^^^^^^^^
+
+The :cpp:class:`LabelMap <LAMMPS_NS::LabelMap>` class provides a two way
+mapping between symbolic type labels in input and output files and
+numeric types as they are used by LAMMPS internally.  Instead of
+changing the numeric types in files to satisfy the requirements from
+LAMMPS for a given application, the symbolic types can remain and only
+the label map needs to be adjusted.  When following the convention that
+the labels for bonded interactions are created by joining the
+constituent atom types with hyphens, this can significantly improve
+readability, maintainability, and re-usability of inputs and reduces the
+chance of errors.
+
+The LabelMap class also provides automatic type inference for bonded
+interactions based on their constituent atom types.  For instance, based
+on the atom type labels, the corresponding bond, angle, dihedral, or
+improper types can be inferred provided the corresponding type labels
+follow the convention that they are composed of the symbolic atom types
+connected by hyphens.
+
+Integration with utils namespace
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Several utility functions in the ``utils`` namespace work with type labels and
+interact with the LabelMap class:
+
+* :cpp:func:`utils::is_type() <LAMMPS_NS::utils::is_type>` - Validates whether
+  a string is a valid type label.
+
+* :cpp:func:`utils::expand_type() <LAMMPS_NS::utils::expand_type>` - Converts
+  a type label string to its numeric equivalent using the LabelMap.
+
+* :cpp:func:`utils::bounds_typelabel() <LAMMPS_NS::utils::bounds_typelabel>` -
+  Extended version of ``utils::bounds()`` that accepts type labels in addition
+  to numeric ranges. Uses ``expand_type()`` internally to convert labels to
+  numeric bounds before processing.
+
+These functions enable seamless integration of type labels throughout LAMMPS,
+allowing commands that accept type specifications to work with both numeric
+indices and symbolic labels. Below are some code examples.
+
+Finding types from labels and vice versa
+""""""""""""""""""""""""""""""""""""""""
+
+.. code-block:: c++
+
+   #include "label_map.h"
+   #include "atom.h"
+
+   // assuming this code is used inside a class that is derived from LAMMPS_NS::Pointers
+   LabelMap *lmap = atom->lmap;
+
+   // Forward lookup: Get numeric type from label
+   int ctype = lmap->find_type("C", Atom::ATOM);    // Returns atom type for "C"
+   int htype = lmap->find_type("H", Atom::ATOM);    // Returns atom type for "H"
+   int missing = lmap->find_type("X", Atom::ATOM);  // Returns -1 (not found)
+
+   // Reverse lookup: Get label from numeric type
+   const std::string &label1 = lmap->find_label(1, Atom::ATOM);  // Returns label for type 1
+   const std::string &label2 = lmap->find_label(2, Atom::BOND);  // Returns bond label for type 2
+
+   // Check if all types have labels
+   bool complete = lmap->is_complete(Atom::ATOM);  // Returns true if all atom types labeled
+
+Inferring bonded types from atom types
+""""""""""""""""""""""""""""""""""""""
+
+.. code-block:: c++
+
+   #include "label_map.h"
+   #include "atom.h"
+
+   LabelMap *lmap = atom->lmap;
+
+   // Assume we have: labelmap atom 1 C 2 H 3 N
+   // And: labelmap bond 1 C-H 2 C-C 3 C-N
+
+   // Infer bond type from numeric atom types
+   int bt1 = lmap->infer_bondtype(1, 2);  // Returns 1 (C-H bond)
+   int bt2 = lmap->infer_bondtype(1, 1);  // Returns 2 (C-C bond)
+   int bt3 = lmap->infer_bondtype(3, 1);  // Returns 3 (C-N bond, symmetric match)
+
+   // Infer bond type from atom type labels (handles symmetry automatically)
+   int bt4 = lmap->infer_bondtype({"C", "H"});  // Returns 1 (C-H)
+   int bt5 = lmap->infer_bondtype({"H", "C"});  // Returns 1 (symmetric match)
+   int bt6 = lmap->infer_bondtype({"C", "N"});  // Returns 3 (C-N)
+
+Validating and expanding type labels
+""""""""""""""""""""""""""""""""""""
+
+.. code-block:: c++
+
+   #include "utils.h"
+   #include "lammps.h"
+
+   using namespace LAMMPS_NS;
+
+   LAMMPS *lmp = /* ... */;
+
+   // Validate type label strings
+   int result1 = utils::is_type("C");     // Returns 1 (valid label)
+   int result2 = utils::is_type("123");   // Returns 0 (numeric type)
+   int result3 = utils::is_type("*");     // Returns -1 (invalid - starts with *)
+   int result4 = utils::is_type("C H");   // Returns -1 (invalid - contains whitespace)
+
+   // Convert type label to numeric string
+   char *numstr = utils::expand_type(FLERR, "C", Atom::ATOM, lmp);
+   if (numstr) {
+       // Use the numeric type string
+       delete[] numstr;  // Must delete after use
+   }
+
+   // Convert type label to integer
+   int type = utils::expand_type_int(FLERR, "C", Atom::ATOM, lmp, true);
+   // The 'true' argument enables range verification
+
+   // Use bounds_typelabel for ranges with label support
+   int lo, hi;
+   utils::bounds_typelabel(FLERR, "C:H", 1, 10, lo, hi, lmp, Atom::ATOM);
+   // Expands "C:H" to numeric range, e.g., "1:2" -> lo=1, hi=2
+
+----------
+
+LabelMap class reference
+^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. doxygenclass:: LAMMPS_NS::LabelMap
+   :project: progguide
+   :members:
+
+----------
+
 Memory pool classes
 -------------------
 
diff --git a/doc/src/Howto_viz.rst b/doc/src/Howto_viz.rst
index ea80da342c9..6f4239bcbf9 100644
--- a/doc/src/Howto_viz.rst
+++ b/doc/src/Howto_viz.rst
@@ -34,6 +34,83 @@ rarely needed these days.
 
 ------------------------
 
+Basic workflow for loading LAMMPS trajectories in VMD
+=====================================================
+
+VMD can read native LAMMPS dump files (in text format not binary) and
+several other dump styles.  The native LAMMPS format is preferred since
+it contains simulation box information.  VMD does not support reading
+data files directly, but `the TopoTools plugin
+<https://www.ks.uiuc.edu/Research/vmd/plugins/topotools/>`_ does.  This
+step is usually necessary, since the LAMMPS dump files do not contain
+information about the molecular topology and elements and that
+information can be obtained from the data file.  There also is the
+`pbctools plugin
+<https://www.ks.uiuc.edu/Research/vmd/plugins/pbctools/>`_ that can be
+used for "repairing" bonds that are broken because one of its atom
+passed through a periodic boundary.  Because of using the plugins, the
+following commands need to be typed into the VMD console to read a
+trajectory.  The following commands are based on the ``peptide`` example
+to which the command ``dump 1 all atom 10 dump.peptide`` was added.
+Note the use of the group "all" in contrast to the commented out dump
+commands that use the group "peptide".  This is to match the number of
+atoms in the data file.  VMD does not support cases where the number of
+atoms change or you are trying to read in a subset of a system unless
+*all* files contain the same subset.
+
+.. code-block:: Tcl
+
+   # load data file including its coordinates
+   topo readlammpsdata data.peptide full
+   # atom style full is the default ^^^^ and can be omitted
+
+   # try to infer missing atom properties from available information
+   topo guessatom lammps data
+
+   # now add one of more dump files
+   mol addfile dump.peptide type lammpstrj waitfor all
+   # we tell VMD which file type ^^^^^^^^^ if the filename does not end in .lammpstrj
+
+   # jump to the first frame with the data file coordinates
+   animate goto 0
+   # if needed  use the following command to try and repair broken bonds
+   # pbc join fragment -now
+
+   # unwrap the entire trajectory to repair all broken bonds
+   # this is faster and more accurate than running `pbc join` for all frames
+   pbc unwrap -all
+   pbc wrap -all -compound fragment -orthorhombic
+   #       remove for triclinic box ^^^^^^^^^^^^^
+
+   # delete the coordinate frame from the data file
+   animate delete beg 0 end 0 top
+   # write out the topology information without coordinates to a PSF file
+   animate write psf peptide.psf top
+   # and the processed trajectory data to a DCD file for future use
+   animate write dcd peptide.dcd top waitfor all
+
+   ###############################################################
+   # the steps up to this line only need to be run once
+   ###############################################################
+   # set visualization defaults
+   display perspective orthographic
+   mol default style Licorice
+   mol default material Diffuse
+   display backgroundgradient on
+
+   # delete molecule and load psf and dcd file
+   mol delete top
+   mol new peptide.psf
+   mol addfile peptide.dcd waitfor all
+   pbc box
+
+To automate this process, you can also save these commands to a file,
+e.g. ``peptide.vmd`` and then load this file from the VMD "File" menu
+with "Load Visualization State..."  or type in the command console
+``play peptide.vmd``.
+
+------------------------
+
 Advanced graphics features in the *dump image* command
 ======================================================
 
diff --git a/doc/src/Manual.rst b/doc/src/Manual.rst
index c8ada64aeb3..1bf7ae11b78 100644
--- a/doc/src/Manual.rst
+++ b/doc/src/Manual.rst
@@ -81,9 +81,9 @@ The manual is organized into three parts:
 
 .. _user_documentation:
 
-************
+**********
 User Guide
-************
+**********
 
 .. toctree::
    :maxdepth: 2
@@ -106,9 +106,9 @@ User Guide
 
 .. _programmer_documentation:
 
-******************
+****************
 Programmer Guide
-******************
+****************
 
 .. toctree::
    :maxdepth: 2
diff --git a/doc/src/create_atoms.rst b/doc/src/create_atoms.rst
index 07e19a30732..2ccf82ba5b4 100644
--- a/doc/src/create_atoms.rst
+++ b/doc/src/create_atoms.rst
@@ -80,10 +80,10 @@ Examples
    create_atoms Pt box
 
    labelmap atom 1 C 2 Si
-   create_atoms C region regsphere basis Si C
+   create_atoms C region regsphere basis 1 Si
 
-   create_atoms 3 region regsphere basis 2 3
-   create_atoms 3 region regsphere basis 2 3 ratio 0.5 74637
+   create_atoms 2 region regsphere basis 2 3
+   create_atoms 1 region regsphere basis 2 2 ratio 0.5 74637
    create_atoms 3 single 0 0 5 group newatom
    create_atoms 1 box var v set x xpos set y ypos
    create_atoms 2 random 50 12345 NULL overlap 2.0 maxtry 50
@@ -120,7 +120,8 @@ added to the specified atom *type* (e.g., if *type* = 2 and the file
 specifies atom types 1, 2, and 3, then each created molecule will have
 atom types 3, 4, and 5).
 
-.. note::
+.. admonition:: Atoms outside the box
+   :class: note
 
    You cannot use this command to create atoms that are outside the
    simulation box; they will just be ignored by LAMMPS.  This is true
@@ -378,10 +379,20 @@ This is the meaning of the other optional keywords.
 
 The *basis* keyword is only used when atoms (not molecules) are being
 created.  It specifies an atom type that will be assigned to specific
-basis atoms as they are created.  See the :doc:`lattice <lattice>`
-command for specifics on how basis atoms are defined for the unit cell
-of the lattice.  By default, all created atoms are assigned the
-argument *type* as their atom type.
+basis atoms as they are created that *overrides* the default atom type.
+See the :doc:`lattice <lattice>` command for specifics on how basis
+atoms are defined for the unit cell of the lattice.  By default, all
+created atoms are assigned the argument *type* as their atom type.  To
+illustrate this please see the following 4 example *create_atom* command
+lines that all create the same system:
+
+.. code-block:: LAMMPS
+
+   lattice diamond 4.36
+   create_atoms 1 box basis 1 1 basis 2 2 basis 3 1 basis 4 2 basis 5 1 basis 6 2 basis 7 1 basis 8 2
+   create_atoms 2 box basis 1 1 basis 2 2 basis 3 1 basis 4 2 basis 5 1 basis 6 2 basis 7 1 basis 8 2
+   create_atoms 1 box basis 2 2 basis 4 2 basis 6 2 basis 8 2
+   create_atoms 2 box basis 1 1 basis 3 1 basis 5 1 basis 7 1
 
 The *ratio* and *subset* keywords can be used in conjunction with the
 *box* or *region* styles to limit the total number of particles
diff --git a/doc/src/dump_image.rst b/doc/src/dump_image.rst
index 9e22a005f15..90b5e0e391a 100644
--- a/doc/src/dump_image.rst
+++ b/doc/src/dump_image.rst
@@ -635,7 +635,7 @@ commands are in the :doc:`Howto_viz` howto.
 
 .. versionchanged:: TBD
 
-   draw style *transparency* was added
+   draw style *transparent* was added
 
 The *region* keyword can be used to create a graphical representation of
 a :doc:`region <region>`.  This can be helpful in debugging the location
@@ -645,7 +645,7 @@ region-ID, the color for drawing the region, the draw style, and
 possible additional arguments as required by the draw style.
 
 Four draw styles of representing a region are available: *filled*\,
-*transparency*\, *frame*\, and *points*.  With draw style *filled* the
+*transparent*\, *frame*\, and *points*.  With draw style *filled* the
 surface of the region is triangulated and drawn.  For region styles that
 support open faces, surfaces for such open faces are skipped.  The style
 *transparent* is like *filled* but takes an additional parameter in the
@@ -844,15 +844,16 @@ color map.  The color map is used to assign a specific RGB
 based on the atom's attribute, which is a numeric value, e.g. its
 x-component of velocity if the atom-attribute "vx" was specified.
 
-The basic idea of a color map is that the atom-attribute will be
-within a range of values, and that range is associated with a series
-of colors (e.g. red, blue, green).  An atom's specific value (vx =
--3.2) can then mapped to the series of colors (e.g. halfway between
-red and blue), and a specific color is determined via an interpolation
-procedure.
+The basic idea of a color map is that the atom-attribute will be within
+a range of values, and that range is associated with a series of colors
+(e.g. red, blue, green).  An atom's specific value (vx = -3.2) can then
+mapped to the series of colors (e.g. halfway between red and blue), and
+a specific color is determined via an interpolation procedure.  There
+are some example command lines and resulting images at the end of this
+paragraph.
 
-There are many possible options for the color map, enabled by the
-*amap* keyword.  Here are the details.
+There are many possible options for the color map, enabled by the *amap*
+keyword.  Here are the details.
 
 The *lo* and *hi* settings determine the range of values allowed for
 the atom attribute.  If numeric values are used for *lo* and/or *hi*,
@@ -946,23 +947,71 @@ green.  The color of the atom is the color of its bin.  Note that the
 sequential color map is really a shorthand way of defining a discrete
 color map without having to specify where all the bin boundaries are.
 
-Here is an example of using a sequential color map to color all the
-atoms in individual molecules with a different color.  See the
-examples/pour/in.pour.2d.molecule input script for an example of how
-this is used.
+Here is an example for using a sequential color map to color all the
+atoms in individual molecules with a different color.  See below for how
+this can be used in the ``examples/pour/in.pour.2d.molecule`` input
+script.
 
 .. code-block:: LAMMPS
 
-   variable        colors string &
-                   "red green blue yellow white &
-                   purple pink orange lime gray"
-   variable        mol atom mol%10
-   dump            1 all image 250 image.*.jpg v_mol type &
-                   zoom 1.6 adiam 1.5
-   dump_modify     1 pad 5 amap 0 10 sa 1 10 ${colors}
+   variable    colors string "red green blue yellow white purple pink orange lime gray"
+   variable    mol2 atom mol%10
+   dump        2 all image 250 image.*.png v_mol2 type region slab black frame 0.25 &
+                               zoom 3.5 adiam 1.4 size 1200 600 fsaa yes shiny 0.2
+   dump_modify 2 pad 5 amap 0 10 sa 1 10 ${colors} backcolor darkgray boxcolor silver
 
-In this case, 10 colors are defined, and molecule IDs are
-mapped to one of the colors, even if there are 1000s of molecules.
+In this case, 10 colors are defined, and molecule IDs are mapped to one
+of the colors, even if there are 1000s of molecules.
+
+Here is an example for coloring the atoms in the "melt" example by their
+velocity with a custom continuous color map and using :doc:`fix
+graphics/labels <fix_graphics_labels>` to generate a colormap legend:
+
+.. code-block:: LAMMPS
+
+   # compute atom velocity
+   variable vel atom sqrt(vx*vx+vy*vy+vz*vz)
+
+   # overlay the top of the image with a horizontal color scale legend
+   fix obj all graphics/labels 100 colorscale "viz" "Atom Velocity (sigma/tau)" 300.0 560.0 0.0 size 24 &
+               transcolor none framecolor darkgray backcolor darkgray length 800
+
+   # output images and set atom color by the value of the variable "vel"
+   dump viz all image 100 melt-*.png v_vel type size 600 600 zoom 1.4 shiny 0.2 view 85 -5 &
+                          fsaa yes box yes 0.025 center s 0.5 0.5 0.6 fix obj const 1 0
+   dump_modify viz pad 6 boxcolor lightskyblue backcolor darkgray backcolor2 silver adiam * 1.2
+
+   # customize the color map using a continuous map with fractions
+   dump_modify viz amap 0.0 8 cf 0.0 6 min red 0.2 organge 0.4 green 0.6 darkcyan 0.8 blue max purple
+
+This is an altered *dump_modify* command line to generate a sequential color map:
+
+.. code-block:: LAMMPS
+
+   dump_modify viz amap 0.5 5.5 sf 0.167 6 red orange green darkcyan blue purple
+
+And another altered *dump_modify* command line to generate a discrete color map using absolute values:
+
+.. code-block:: LAMMPS
+
+   dump_modify viz amap 0.5 5.5 da 0.0 6 min 1.0 red 1.0 2.0 orange 2.0 3.0 green 3.0 4.0 darkcyan
+
+.. |amap1| image:: img/amap1.png
+   :width: 38%
+.. |amap2| image:: img/amap2.png
+   :width: 19%
+.. |amap3| image:: img/amap3.png
+   :width: 19%
+.. |amap4| image:: img/amap4.png
+   :width: 19%
+
+Here are images of the examples from above.
+
+|amap1|  |amap2|  |amap3|  |amap4|
+
+.. raw:: html
+
+   <center>(Click to see the full-size images)</center>
 
 ----------
 
diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst
index 5103722f0d2..d30786fcfac 100644
--- a/doc/src/fix_bond_react.rst
+++ b/doc/src/fix_bond_react.rst
@@ -396,20 +396,21 @@ Wildcard atoms match to any atom type in the simulation.  Wildcard atoms
 can be used to reduce the number of reaction templates needed to model a
 set of similar reactions.  Wildcard atoms are specified in the Wildcards
 section of the map file.  The atom types of wildcard atoms in the
-simulation are not updated.  Any bond, angle, dihedral, or improper, that
-is defined in the reaction templates and contains a wildcard atom, will be
-updated by inferring its type from its constituent atom types.  To use
-wildcard atoms, a specific :doc:`type label <Howto_type_labels>` format is
-necessary to infer the types of higher-order interactions.  Bond, angle,
-dihedral, and improper type labels must contain their constituent atom
-types delimited by hyphens, e.g., ‘c2-c2-c2-n’ for a dihedral that contains
-three atoms of type 'c2' and one atom of 'n'.  Certain symmetries are
-considered to account for equivalent ways of writing higher-order
-interactions.  Type labels for bonds, angles, and dihedrals are assumed to
-be equivalent to those written in reverse order.  For example, an angle
-with type label 'c1-c2-n' is equivalent to 'n-c2-c1'.  Symmetries for
-impropers are more complex and are described on the doc page for each
-improper style in the 'Symmetry convention' section.
+simulation are not updated.  Any bond, angle, dihedral, or improper,
+that is defined in the reaction templates and contains a wildcard atom,
+will be updated by inferring its type from its constituent atom types.
+To use wildcard atoms, a specific :doc:`type label <Howto_type_labels>`
+format is necessary to infer the types of higher-order interactions.
+Bond, angle, dihedral, and improper type labels must contain their
+constituent atom types delimited by hyphens, e.g., 'c2-c2-c2-n' for a
+dihedral that contains three atoms of type 'c2' and one atom of 'n'.
+Certain symmetries are considered to account for equivalent ways of
+writing higher-order interactions.  Type labels for bonds, angles, and
+dihedrals are assumed to be equivalent to those written in reverse
+order.  For example, an angle with type label 'c1-c2-n' is equivalent to
+'n-c2-c1'.  Symmetries for impropers are more complex and are described
+on the doc page for each improper style in the 'Symmetry convention'
+section.
 
 The post-reacted molecule template contains a sample of the reaction
 site and its surrounding topology after the reaction has occurred. It
diff --git a/doc/src/fix_graphics_labels.rst b/doc/src/fix_graphics_labels.rst
index 87cd57ae0c2..986e13cb383 100644
--- a/doc/src/fix_graphics_labels.rst
+++ b/doc/src/fix_graphics_labels.rst
@@ -14,7 +14,7 @@ Syntax
 * graphics/labels = style name of this fix command
 * Nevery = update graphics information every this many time steps
 * zero or more keyword/args groups may be appended
-* keyword = *image* or *text*
+* keyword = *image* or *text* or *colorscale*
 
   .. parsed-literal::
 
@@ -35,7 +35,7 @@ Syntax
         x, y, z  = position where the center of the text is located in the visualization
         any of x, y, or z can be a variable (see below)
 
-        keyword = *fontcolor* or *framecolor* or *backcolor* or *transcolor* or *size*
+        keyword = *fontcolor* or *framecolor* or *backcolor* or *transcolor* or *size* or *horizontal* or *vertical*
           *fontcolor* arg = select color for text: *white* (default) or *black* or *r/g/b*
              *white* = uses white
              *black* = uses black
@@ -60,6 +60,42 @@ Syntax
              *none* = disables transparency
              *r/g/b* = provide three integers in the range 0 to 255
           *size* value = set the size of the characters (default 24), can be a variable (see below)
+          *horizontal* = create horizontal text label
+          *vertical* = create vertical text label
+
+     *colorscale* dump-ID titletext x y z keyword args = display a colormap label in visualization
+        labeltext = text for the legend of the colormap label, must be quoted if it contains whitespace
+        x, y, z  = position where the center of the colormap label is located in the visualization
+        any of x, y, or z can be a variable (see below)
+
+        keyword = *fontcolor* or *framecolor* or *backcolor* or *transcolor* or *size* or *horizontal* or *vertical*
+          *fontcolor* arg = select color for text: *white* (default) or *black* or *r/g/b*
+             *white* = uses white
+             *black* = uses black
+             *r/g/b* = provide three integers in the range 0 to 255
+          *framecolor* arg = select color for frame around text: *silver* (default) or *darkgray* or *white* or *black* or *r/g/b*
+             *silver* = uses a very light gray
+             *darkgray* = uses a very dark gray
+             *white* = uses white
+             *black* = uses black
+             *r/g/b* = provide three integers in the range 0 to 255
+          *backcolor* arg = select color for background of the text: *silver* (default) or *darkgray* or *white* or *black* *r/g/b*
+             *silver* = uses a very light gray
+             *darkgray* = uses a very dark gray
+             *white* = uses white
+             *black* = uses black
+             *r/g/b* = provide three integers in the range 0 to 255
+          *transcolor* arg = select color for transparency: *silver* (default) or *darkgray* or *white* or *black* or *none* or *r/g/b*
+             *silver* = uses a very light gray
+             *darkgray* = uses a very dark gray
+             *white* = uses white
+             *black* = uses black
+             *none* = disables transparency
+             *r/g/b* = provide three integers in the range 0 to 255
+          *size* value = set the size of the characters (default 24), can be a variable (see below)
+          *length* value = approximate minimal length of the colorscale label
+          *horizontal* = create horizontal text label
+          *vertical* = create vertical text label
 
 Examples
 """"""""
@@ -70,6 +106,8 @@ Examples
    fix pot all graphics/labels 100 image teapot.ppm 1.0 v_ypos v_zpos scale v_prog transcolor 19/92/192
    fix lbl all graphics/labels 1000 text "LAMMPS graphics demo" 5.0 -1.0 -2.0 backcolor darkgray framecolor black
    fix info all graphics/labels 1000 text "Step: $(step)  Angle: ${rot}" 5.0 -1.0 -2.0 size 32
+   fix obj all graphics/labels 200 colorscale viz "Atom Velocity" 20.0 6.5 13.0 size 32 length 1000 &
+                                     transcolor none framecolor white backcolor darkgray tics 12
 
 Description
 """""""""""
@@ -79,11 +117,23 @@ Description
 This fix allows to add either images or text as "labels" to :doc:`dump
 image <dump_image>` created images by using the *fix* keyword.  This can
 be useful to augment images with additional graphics or text directly
-and without having to post-process the images.  Since the positioning
-uses the coordinate system of the simulation and because the graphics
-objects use the depth buffer of the image rasterizer, atoms and other
-graphics in the "scene" can be located before or behind any text or
-image label.
+and without having to post-process the images.  The positions can be
+either interpreted as coordinates in the simulation box or as
+coordinates in the coordinate system of the image.  The selection is
+made by setting the *fflag1* keyword in the :doc:`dump image fix
+<dump_image>` command (see the "Dump image info" section below).  When
+the positioning uses the coordinate system of the simulation the
+distance of the graphics objects from the camera is determined from the
+given z-coordinate and atoms or other graphics objects in the "scene"
+can be located in front of or behind any *image*, *text* or *colorscale*
+label.  The label is *always* parallel to the image plane.
+
+When the image coordinate system is used, the labels are *always* on
+top, and if two labels are overlapping, the label that is added to the
+image *first* will be on top of the other.  That order cannot be changed
+within the same fix, but you can use multiple fix commands and then the
+order of the *fix* keywords in the *dump image"* command line determines
+the order and thus which label is drawn on top of the other.
 
 The *group-id* is ignored by this fix.
 
@@ -96,14 +146,16 @@ compatible.
 The *image* keyword reads an image file and adds it to the visualization
 centered around the provided position and optionally scaled by the
 provided scale factor.  The filename suffix determines whether LAMMPS
-will try to read a file in JPEG, PNG, or PPM format.  If the suffix is
-".jpg" or ".jpeg", then LAMMPS attempts to read the image in `JPEG
-format <jpeg_format_>`_, if the suffix is ".png", then lammps attempts
-to read the image in `PNG format <png_format_>`_.  Otherwise LAMMPS will
-try to read the image in `ppm (aka netpbm) format <ppm_format_>`_.  Not
-all variants of those file formats are compatible with image reader code
-in LAMMPS.  If LAMMPS encounters an incompatible or unrecognizable file
-format or a corrupted file, it will stop with an error.
+will try to read a file in JPEG, PNG, TGA, or PPM format.  If the suffix
+is ".jpg" or ".jpeg", then LAMMPS attempts to read the image in `JPEG
+format <jpeg_format_>`_, if the suffix is ".png", then LAMMPS attempts
+to read the image in `PNG format <png_format_>`_, and if the suffix is
+".tga" then LAMMPS will read the file in `TGA format <tga_format_>`_.
+Otherwise LAMMPS will try to read the image in `ppm (aka netpbm) format
+<ppm_format_>`_.  Not all variants of those file formats are compatible
+with image reader code in LAMMPS.  If LAMMPS encounters an incompatible
+or unrecognizable file format or a corrupted file, it will stop with an
+error.
 
 If LAMMPS detects during a run that the file has been changed, it will
 re-read it.  This allows for instance to create a plot using internal
@@ -166,11 +218,94 @@ are required arguments.  Optional keyword / value pairs may be added:
     processing.
 
   When rendering text with transparent background it is recommended to
-  select a similar color but slightly darker or brighter color as background.
-  This will reduce unwanted color effects at the edges due to anti-aliasing.
+  select a similar color but slightly darker or brighter color as
+  background.  This will reduce unwanted color effects at the edges due
+  to anti-aliasing.
+
+  The *horizontal* keyword selects creating a horizontal text label
+  (this is the default setting).  The *vertical* keyword selects
+  creating a vertical text label instead.
+
+The *colorscale* keyword will create a colormap legend indicating the
+mapping of values to the color of atoms in the :doc:`dump image
+<dump_image>` instance with the given dump-ID and adds it to the
+visualization centered around the provided position in a similar fashion
+as with the *image* or *text* keywords.  The requirements for the text
+argument are the same as in the :doc:`fix print <fix_print>` command: it
+must be a single argument, so text with whitespace must be quoted; and
+the text may contain equal style or immediate variables using the
+``${name}`` or ``$(expression)`` format.  The variables are evaluated
+and expanded at every *Nevery* time step.  The text is shown in the
+center of and above the colormap.  To the left from the text is the
+lower boundary value and to the right the upper boundary value.  The
+colors are created by a linear interpolation between the lower and upper
+boundary value and writing out pixels in the corresponding color.  The
+fix will receive the actual values from the dump with the given
+*dump-ID*.
+
+.. admonition:: Dynamic color maps
+   :class: note
+
+  When using a dynamic color map with "min" or "max" as the upper or
+  lower range values of the map, the dump will execute only *after* the
+  fix, and thus the upper and lower boundary values will be those from
+  the *previous* step where the dump created an image. will be
+  determined every time the fix is executed and the numbers updated
+  accordingly.  Thus when adding a *colorscale* label with this fix it
+  is generally recommended to use a map with a fixed range. This is
+  especially true when creating movies as a fixed range prevents the
+  color scale label to shrink or grow due to the different width of
+  characters.
+
+When using the *colorscale* keyword, the dump-ID, text and its position
+in the "scene" are required arguments.  Optional keyword / value pairs
+may be added:
+
+  The *size* value determines the size of the letters in the text in
+  pixels (approximately) and values between 4 and 512 are accepted.  The
+  default value is 24. The size (height and width) of the colorbar
+  follows the size of the text.
+
+  The *length* value allows to set a minimal length of the colorscale
+  label.  For technical reasons, this is not exactly enforced, but
+  rather a rough approximation that is used to determine the amount of
+  padding in the text.
 
-There may be multiple *image* or *text* keywords with their arguments
-in a single fix *graphics/labels* command.
+  The *tics* value determines how many "tics" or lines separating the
+  colors are drawn.  This can simplify determining which value a
+  specific color corresponds to.
+
+  There are four color settings: *fontcolor* or *framecolor* or
+  *backcolor* or *transcolor*\ .  The color can be specified for all of
+  those either as an R/G/B triple with values ranging from 0 to 255 for
+  each channel (e.g. yellow would be "255/255/0").  There are also a few
+  shortcuts for common choices: *silver*, *darkgray*, *white*, *black*.
+  The default *backcolor* value is *silver*.
+
+  - *fontcolor* selects the color for the text, the border of the
+    colorbar and the tics. default is *white*
+  - *backcolor* selects the color for the background, default is
+    *silver*
+  - *framecolor* selects the color for the frame around the background,
+    default is *silver*.
+  - *transcolor* value selects a color for transparency, default is
+    *silver*.  If this color is the same as any of the other color
+    settings, those pixels are not drawn.  Thus with the default
+    settings, the text will be rendered in white without background or
+    frame.  The *none* setting for *transcolor* disables transparency
+    processing.
+
+  When rendering text with transparent background it is recommended to
+  select a similar color but slightly darker or brighter color as
+  background.  This will reduce unwanted color effects at the edges due
+  to anti-aliasing.
+
+  The *horizontal* keyword selects creating a horizontal colorscale label
+  (this is the default setting).  The *vertical* keyword selects
+  creating a vertical text label instead.
+
+There may be multiple *image* or *text* or *colorscale* keywords with
+their arguments in a single fix *graphics/labels* command.
 
 The arguments for the positions of an *image* or *text* and the *scale*
 factor of an *image* or the *size* of a *text* can be specified as
@@ -185,6 +320,7 @@ more detailed discussion on using variables with graphics objects.
 .. _jpeg_format: https://jpeg.org/jpeg/
 .. _png_format: https://en.wikipedia.org/wiki/portable_network_graphics
 .. _ppm_format: https://en.wikipedia.org/wiki/netpbm
+.. _tga_format: https://en.wikipedia.org/wiki/Truevision_TGA
 
 -----------
 
@@ -200,7 +336,16 @@ The color style setting for the fix in the :doc:`dump image
 transparency is by default fully opaque and can be changed with
 *dump\_modify ftrans*\ .
 
-The *fflag1* and *fflag2* settings of *dump image fix* are ignored.
+The *fflag1* setting of *dump image fix* determines how the coordinates
+for the location of the center of the image or the center of the text
+label are interpreted.  Setting *fflag1* to 0 uses the simulation box
+coordinate system (x, y, and z) while setting *fflag1* to 1 uses the
+image coordinate system where (0,0) is the location of the lower left
+corner and (<image width>, <image height>) the upper right corner.  In
+the latter case, the z-coordinate is ignored and the image or label is
+placed on top of everything.
+
+The *fflag2* settings of *dump image fix* is ignored.
 
 --------------
 
@@ -311,4 +456,4 @@ Related commands
 Default
 """""""
 
-transcolor = "none" for *image* and "silver" for *text*, scale = 1.0, fontcolor = white, backcolor = silver, framecolor = silver, size = 24
+transcolor = "none" for *image* and "silver" for *text*, scale = 1.0, fontcolor = white, backcolor = silver, framecolor = silver, size = 24, horizontal, tics = 0
diff --git a/doc/src/img/amap1.png b/doc/src/img/amap1.png
new file mode 100644
index 00000000000..bdd1238e928
Binary files /dev/null and b/doc/src/img/amap1.png differ
diff --git a/doc/src/img/amap2.png b/doc/src/img/amap2.png
new file mode 100644
index 00000000000..dd5faf17bff
Binary files /dev/null and b/doc/src/img/amap2.png differ
diff --git a/doc/src/img/amap3.png b/doc/src/img/amap3.png
new file mode 100644
index 00000000000..8d513a0ce5f
Binary files /dev/null and b/doc/src/img/amap3.png differ
diff --git a/doc/src/img/amap4.png b/doc/src/img/amap4.png
new file mode 100644
index 00000000000..1a6feb9a320
Binary files /dev/null and b/doc/src/img/amap4.png differ
diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt
index 63336b4101e..9274ebc864a 100644
--- a/doc/utils/sphinx-config/false_positives.txt
+++ b/doc/utils/sphinx-config/false_positives.txt
@@ -118,6 +118,7 @@ anglegrad
 anglelist
 angleoffset
 angletangrad
+angletype
 angmom
 angmomx
 angmomy
@@ -585,6 +586,9 @@ coleman
 Colliex
 collinear
 collisional
+colorbar
+colormap
+colorscale
 Columic
 colvars
 Colvars
@@ -856,6 +860,7 @@ diffusivity
 dihedral
 dihedrallist
 dihedrals
+dihedraltype
 Dihedrals
 dihydride
 Dij
@@ -2114,6 +2119,7 @@ llammps
 lld
 LLVM
 lm
+lmap
 lmp
 lmpptr
 lmpqst
diff --git a/examples/GRAPHICS/README b/examples/GRAPHICS/README
index 138701e5e95..2bd67d1a270 100644
--- a/examples/GRAPHICS/README
+++ b/examples/GRAPHICS/README
@@ -11,5 +11,10 @@ automatically as they are created.
 Here is a list of the inputs and which graphics features they demonstrate
 
 - in.cubes-and-pyramids:    visualization of body particles, uses progress bar
+                            https://www.youtube.com/shorts/OYn_VVodnIg
+
 - in.breakable:             using autobond, color-by-property, adding objects and text
+                            https://www.youtube.com/watch?v=f4hfPs7aCmI
+
 - in.water-arrows:          box of water molecules, uses many kinds of arrows
+                            https://www.youtube.com/shorts/4Cm5p0SfgNU
diff --git a/src/BROWNIAN/fix_brownian.cpp b/src/BROWNIAN/fix_brownian.cpp
index 0b73858f29e..df3cc5f4477 100644
--- a/src/BROWNIAN/fix_brownian.cpp
+++ b/src/BROWNIAN/fix_brownian.cpp
@@ -34,10 +34,9 @@ using namespace FixConst;
 FixBrownian::FixBrownian(LAMMPS *lmp, int narg, char **arg) : FixBrownianBase(lmp, narg, arg)
 {
   if (dipole_flag || gamma_t_eigen_flag || gamma_r_eigen_flag || gamma_r_flag || rot_temp_flag ||
-      planar_rot_flag) {
-    error->all(FLERR, "Illegal fix brownian command.");
-  }
-  if (!gamma_t_flag) { error->all(FLERR, "Illegal fix brownian command."); }
+      planar_rot_flag)
+    error->all(FLERR, "Using keyword(s) not compatible with fix {}", style);
+  if (!gamma_t_flag) error->all(FLERR, "Keyword gamma_t is required for fix {}", style);
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/BROWNIAN/fix_brownian_asphere.cpp b/src/BROWNIAN/fix_brownian_asphere.cpp
index 5ce57df7ce8..0dd4da95e7a 100644
--- a/src/BROWNIAN/fix_brownian_asphere.cpp
+++ b/src/BROWNIAN/fix_brownian_asphere.cpp
@@ -37,17 +37,14 @@ using namespace FixConst;
 FixBrownianAsphere::FixBrownianAsphere(LAMMPS *lmp, int narg, char **arg) :
     FixBrownianBase(lmp, narg, arg), avec(nullptr)
 {
-  if (!gamma_t_eigen_flag || !gamma_r_eigen_flag) {
-    error->all(FLERR, "Illegal fix brownian command.");
-  }
-
-  if (gamma_t_flag || gamma_r_flag) error->all(FLERR, "Illegal fix brownian command.");
+  if (gamma_t_flag || gamma_r_flag)
+    error->all(FLERR, "Keywords gamma_t or gamma_r are not compatible with fix {}", style);
+  if (!gamma_t_eigen_flag || !gamma_r_eigen_flag)
+    error->all(FLERR, "Keywords gamma_t_eigen and gamma_r_eigen are required for fix {}", style);
 
   if (dipole_flag && !atom->mu_flag)
-    error->all(FLERR, "Fix brownian/asphere dipole requires atom attribute mu");
-
-  if (!atom->ellipsoid_flag)
-    error->all(FLERR, "Fix brownian/asphere requires atom style ellipsoid");
+    error->all(FLERR, "Fix {} dipole requires atom attribute mu", style);
+  if (!atom->ellipsoid_flag) error->all(FLERR, "Fix {} requires atom style ellipsoid", style);
 
   if (planar_rot_flag && (comm->me == 0)) {
     error->warning(FLERR, "Ignoring first two entries of gamma_r_eigen since rotation is planar.");
diff --git a/src/BROWNIAN/fix_brownian_base.cpp b/src/BROWNIAN/fix_brownian_base.cpp
index 6b00ec29854..35a8424e747 100644
--- a/src/BROWNIAN/fix_brownian_base.cpp
+++ b/src/BROWNIAN/fix_brownian_base.cpp
@@ -48,18 +48,19 @@ FixBrownianBase::FixBrownianBase(LAMMPS *lmp, int narg, char **arg) :
   planar_rot_flag = 0;
   g2 = 0.0;
 
-  if (narg < 5) utils::missing_cmd_args(FLERR, "fix brownian", error);
+  std::string mystyle = fmt::format("fix {}", style);
+  if (narg < 5) utils::missing_cmd_args(FLERR, mystyle, error);
 
   temp = utils::numeric(FLERR, arg[3], false, lmp);
-  if (temp <= 0) error->all(FLERR, "Fix brownian temp must be > 0.0");
+  if (temp <= 0) error->all(FLERR, 3, "Fix {} temp must be > 0.0", style);
 
   seed = utils::inumeric(FLERR, arg[4], false, lmp);
-  if (seed <= 0) error->all(FLERR, "Fix brownian seed must be > 0");
+  if (seed <= 0) error->all(FLERR, 4, "Fix {} seed must be > 0", style);
 
   int iarg = 5;
   while (iarg < narg) {
     if (strcmp(arg[iarg], "rng") == 0) {
-      if (narg < iarg + 1) utils::missing_cmd_args(FLERR, "fix brownian rng", error);
+      if (narg < iarg + 2) utils::missing_cmd_args(FLERR, mystyle + " rng", error);
       if (strcmp(arg[iarg + 1], "uniform") == 0) {
         noise_flag = 1;
       } else if (strcmp(arg[iarg + 1], "gaussian") == 0) {
@@ -68,11 +69,11 @@ FixBrownianBase::FixBrownianBase(LAMMPS *lmp, int narg, char **arg) :
       } else if (strcmp(arg[iarg + 1], "none") == 0) {
         noise_flag = 0;
       } else {
-        error->all(FLERR, "Unknown fix brownian rng keyword {}", arg[iarg + 1]);
+        error->all(FLERR, iarg + 1, "Unknown fix {} rng keyword {}", style, arg[iarg + 1]);
       }
       iarg = iarg + 2;
     } else if (strcmp(arg[iarg], "dipole") == 0) {
-      if (narg < iarg + 3) utils::missing_cmd_args(FLERR, "fix brownian dipole", error);
+      if (narg < iarg + 4) utils::missing_cmd_args(FLERR, mystyle + " dipole", error);
 
       dipole_flag = 1;
       delete[] dipole_body;
@@ -84,27 +85,31 @@ FixBrownianBase::FixBrownianBase(LAMMPS *lmp, int narg, char **arg) :
       iarg = iarg + 4;
 
     } else if (strcmp(arg[iarg], "gamma_t_eigen") == 0) {
-      if (narg < iarg + 3) utils::missing_cmd_args(FLERR, "fix brownian gamma_t_eigen", error);
+      if (narg < iarg + 4) utils::missing_cmd_args(FLERR, mystyle + " gamma_t_eigen", error);
+
+      double gamma_t_tmp[3];
+      gamma_t_tmp[0] = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
+      gamma_t_tmp[1] = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
+      if (domain->dimension == 2) {
+        if (strcmp(arg[iarg + 3], "inf") != 0)
+          error->all(FLERR, iarg + 3, "Fix {} gamma_t_eigen third value must be inf for 2D system",
+                     style);
+        gamma_t_tmp[2] = 1.0;
+      } else {
+        gamma_t_tmp[2] = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
+      }
+      if ((gamma_t_tmp[0] <= 0.0) || (gamma_t_tmp[1] <= 0.0) || (gamma_t_tmp[2] <= 0.0))
+        error->all(FLERR, iarg, "Fix {} gamma_t_eigen values must be > 0", style);
 
       gamma_t_eigen_flag = 1;
       delete[] gamma_t_inv;
       delete[] gamma_t_invsqrt;
       gamma_t_inv = new double[3];
       gamma_t_invsqrt = new double[3];
-      gamma_t_inv[0] = 1. / utils::numeric(FLERR, arg[iarg + 1], false, lmp);
-      gamma_t_inv[1] = 1. / utils::numeric(FLERR, arg[iarg + 2], false, lmp);
-
-      if (domain->dimension == 2) {
-        if (strcmp(arg[iarg + 3], "inf") != 0) {
-          error->all(FLERR, "Fix brownian gamma_t_eigen third value must be inf for 2D system.");
-        }
-        gamma_t_inv[2] = 0;
-      } else {
-        gamma_t_inv[2] = 1.0 / utils::numeric(FLERR, arg[iarg + 3], false, lmp);
-      }
-
-      if (gamma_t_inv[0] < 0 || gamma_t_inv[1] < 0 || gamma_t_inv[2] < 0)
-        error->all(FLERR, "Fix brownian gamma_t_eigen values must be > 0.");
+      gamma_t_inv[0] = 1.0 / gamma_t_tmp[0];
+      gamma_t_inv[1] = 1.0 / gamma_t_tmp[1];
+      gamma_t_inv[2] = 1.0 / gamma_t_tmp[2];
+      if (domain->dimension == 2) gamma_t_inv[2] = 0.0;
 
       gamma_t_invsqrt[0] = sqrt(gamma_t_inv[0]);
       gamma_t_invsqrt[1] = sqrt(gamma_t_inv[1]);
@@ -112,7 +117,27 @@ FixBrownianBase::FixBrownianBase(LAMMPS *lmp, int narg, char **arg) :
       iarg = iarg + 4;
 
     } else if (strcmp(arg[iarg], "gamma_r_eigen") == 0) {
-      if (narg == iarg + 3) error->all(FLERR, "Illegal fix brownian command.");
+      if (narg < iarg + 4) utils::missing_cmd_args(FLERR, mystyle + " gamma_r_eigen", error);
+
+      double gamma_r_tmp[3];
+      if (domain->dimension == 2) {
+        if (strcmp(arg[iarg + 1], "inf") != 0)
+          error->all(FLERR, iarg + 1, "Fix {} gamma_r_eigen first value must be inf for 2D system",
+                     style);
+        gamma_r_tmp[0] = 1.0;
+
+        if (strcmp(arg[iarg + 2], "inf") != 0)
+          error->all(FLERR, iarg + 2, "Fix {} gamma_r_eigen second value must be inf for 2D system",
+                     style);
+        gamma_r_tmp[1] = 1.0;
+      } else {
+        gamma_r_tmp[0] = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
+        gamma_r_tmp[1] = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
+      }
+      gamma_r_tmp[2] = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
+
+      if ((gamma_r_tmp[0] <= 0.0) || (gamma_r_tmp[1] <= 0) || (gamma_r_tmp[2] <= 0))
+        error->all(FLERR, iarg, "Fix {} gamma_r_eigen values must be > 0", style);
 
       gamma_r_eigen_flag = 1;
       delete[] gamma_r_inv;
@@ -120,65 +145,52 @@ FixBrownianBase::FixBrownianBase(LAMMPS *lmp, int narg, char **arg) :
       gamma_r_inv = new double[3];
       gamma_r_invsqrt = new double[3];
 
+      gamma_r_inv[0] = 1.0 / gamma_r_tmp[0];
+      gamma_r_inv[1] = 1.0 / gamma_r_tmp[1];
+      gamma_r_inv[2] = 1.0 / gamma_r_tmp[2];
       if (domain->dimension == 2) {
-        if (strcmp(arg[iarg + 1], "inf") != 0) {
-          error->all(FLERR, "Fix brownian gamma_r_eigen first value must be inf for 2D system.");
-        }
-        gamma_r_inv[0] = 0;
-
-        if (strcmp(arg[iarg + 2], "inf") != 0) {
-          error->all(FLERR, "Fix brownian gamma_r_eigen second value must be inf for 2D system.");
-        }
-        gamma_r_inv[1] = 0;
-      } else {
-
-        gamma_r_inv[0] = 1. / utils::numeric(FLERR, arg[iarg + 1], false, lmp);
-        gamma_r_inv[1] = 1. / utils::numeric(FLERR, arg[iarg + 2], false, lmp);
+        gamma_r_inv[0] = 0.0;
+        gamma_r_inv[1] = 0.0;
       }
-
-      gamma_r_inv[2] = 1. / utils::numeric(FLERR, arg[iarg + 3], false, lmp);
-
-      if (gamma_r_inv[0] < 0 || gamma_r_inv[1] < 0 || gamma_r_inv[2] < 0)
-        error->all(FLERR, "Fix brownian gamma_r_eigen values must be > 0.");
-
       gamma_r_invsqrt[0] = sqrt(gamma_r_inv[0]);
       gamma_r_invsqrt[1] = sqrt(gamma_r_inv[1]);
       gamma_r_invsqrt[2] = sqrt(gamma_r_inv[2]);
       iarg = iarg + 4;
 
     } else if (strcmp(arg[iarg], "gamma_t") == 0) {
-      if (narg == iarg + 1) { error->all(FLERR, "Illegal fix brownian command."); }
+      if (narg < iarg + 2) utils::missing_cmd_args(FLERR, mystyle + " gamma_t", error);
 
       gamma_t_flag = 1;
       gamma_t = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
-      if (gamma_t <= 0) error->all(FLERR, "Fix brownian gamma_t must be > 0.");
+      if (gamma_t <= 0.0) error->all(FLERR, iarg + 1, "Fix {} gamma_t value must be > 0", style);
       iarg = iarg + 2;
 
     } else if (strcmp(arg[iarg], "gamma_r") == 0) {
-      if (narg == iarg + 1) { error->all(FLERR, "Illegal fix brownian command."); }
+      if (narg < iarg + 2) utils::missing_cmd_args(FLERR, mystyle + " gamma_r", error);
 
       gamma_r_flag = 1;
       gamma_r = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
-      if (gamma_r <= 0) error->all(FLERR, "Fix brownian gamma_r must be > 0.");
+      if (gamma_r <= 0) error->all(FLERR, iarg + 1, "Fix {} gamma_r value must be > 0", style);
       iarg = iarg + 2;
 
     } else if (strcmp(arg[iarg], "rotation_temp") == 0) {
-      if (narg == iarg + 1) { error->all(FLERR, "Illegal fix brownian command."); }
+      if (narg < iarg + 2) utils::missing_cmd_args(FLERR, mystyle + " rotation_temp", error);
 
       rot_temp_flag = 1;
       rot_temp = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
-      if (rot_temp <= 0) error->all(FLERR, "Fix brownian rotation_temp must be > 0.");
+      if (rot_temp <= 0)
+        error->all(FLERR, iarg + 1, "Fix {} rotation_temp value must be > 0", style);
       iarg = iarg + 2;
 
     } else if (strcmp(arg[iarg], "planar_rotation") == 0) {
 
       planar_rot_flag = 1;
       if (domain->dimension == 2)
-        error->all(FLERR, "The planar_rotation keyword is not allowed for 2D simulations");
+        error->all(FLERR, iarg, "The planar_rotation keyword is not allowed for 2D simulations");
       iarg = iarg + 1;
 
     } else {
-      error->all(FLERR, "Illegal fix brownian command.");
+      error->all(FLERR, iarg, "Unknown fix {} keyword {}", style, arg[iarg]);
     }
   }
   if (!rot_temp_flag) rot_temp = temp;
diff --git a/src/BROWNIAN/fix_brownian_sphere.cpp b/src/BROWNIAN/fix_brownian_sphere.cpp
index 79e3858f00b..52d89091687 100644
--- a/src/BROWNIAN/fix_brownian_sphere.cpp
+++ b/src/BROWNIAN/fix_brownian_sphere.cpp
@@ -35,12 +35,12 @@ using namespace FixConst;
 FixBrownianSphere::FixBrownianSphere(LAMMPS *lmp, int narg, char **arg) :
     FixBrownianBase(lmp, narg, arg)
 {
-  if (gamma_t_eigen_flag || gamma_r_eigen_flag) {
-    error->all(FLERR, "Illegal fix brownian/sphere command.");
-  }
-
-  if (!gamma_t_flag || !gamma_r_flag) error->all(FLERR, "Illegal fix brownian/sphere command.");
-  if (!atom->mu_flag) error->all(FLERR, "Fix brownian/sphere requires atom attribute mu");
+  if (gamma_t_eigen_flag || gamma_r_eigen_flag)
+    error->all(FLERR, "Keywords gamma_t_eigen or gamma_r_eigen are not compatible with fix {}",
+               style);
+  if (!gamma_t_flag || !gamma_r_flag)
+    error->all(FLERR, "Keywords gamma_t and gamma_r are required for fix {}", style);
+  if (!atom->mu_flag) error->all(FLERR, "Fix {} requires atom attribute mu", style);
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/GRAPHICS/dump_image.cpp b/src/GRAPHICS/dump_image.cpp
index 7f1870eb820..4d4c74f69b8 100644
--- a/src/GRAPHICS/dump_image.cpp
+++ b/src/GRAPHICS/dump_image.cpp
@@ -101,23 +101,23 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
 
   // set filetype based on filename suffix
 
-  if (utils::strmatch(filename, "\\.jpg$") || utils::strmatch(filename, "\\.JPG$")
-      || utils::strmatch(filename, "\\.jpeg$") || utils::strmatch(filename, "\\.JPEG$"))
+  if (utils::strmatch(filename, R"(\.jpg$)") || utils::strmatch(filename, R"(\.JPG$)") ||
+      utils::strmatch(filename, R"(\.jpeg$)") || utils::strmatch(filename, R"(\.JPEG$)"))
     filetype = JPG;
-  else if (utils::strmatch(filename, "\\.png$") || utils::strmatch(filename, "\\.PNG$"))
+  else if (utils::strmatch(filename, R"(\.png$)") || utils::strmatch(filename, R"(\.PNG$)"))
     filetype = PNG;
-  else if (utils::strmatch(filename, "\\.tga$") || utils::strmatch(filename, "\\.TGA$"))
+  else if (utils::strmatch(filename, R"(\.tga$)") || utils::strmatch(filename, R"(\.TGA$)"))
     filetype = TGA;
-  else if (compressed && (utils::strmatch(filename, "\\.jpg\\.\\w+$") ||
-                          utils::strmatch(filename, "\\.JPG\\.\\w+$") ||
-                          utils::strmatch(filename, "\\.jpeg\\.\\w+$") ||
-                          utils::strmatch(filename, "\\.JPEG\\.\\w+$")))
+  else if (compressed && (utils::strmatch(filename, R"(\.jpg\.\w+$)") ||
+                          utils::strmatch(filename, R"(\.JPG\.\w+$)") ||
+                          utils::strmatch(filename, R"(\.jpeg\.\w+$)") ||
+                          utils::strmatch(filename, R"(\.JPEG\.\w+$)")))
     error->all(FLERR, Error::NOLASTLINE, "Cannot use compression with JPEG images");
-  else if (compressed && (utils::strmatch(filename, "\\.png\\.\\w+$") ||
-                           utils::strmatch(filename, "\\.PNG\\.\\w+$")))
+  else if (compressed && (utils::strmatch(filename, R"(\.png\.\w+$)") ||
+                           utils::strmatch(filename, R"(\.PNG\.\w+$)")))
     error->all(FLERR, Error::NOLASTLINE, "Cannot use compression with PNG images");
-  else if (compressed && (utils::strmatch(filename, "\\.tga\\.\\w+$") ||
-                          utils::strmatch(filename, "\\.TGA\\.\\w+$")))
+  else if (compressed && (utils::strmatch(filename, R"(\.tga\.\w+$)") ||
+                          utils::strmatch(filename, R"(\.TGA\.\w+$)")))
     error->all(FLERR, Error::NOLASTLINE, "Cannot use compression with TGA images");
   else filetype = PPM;
 
@@ -1639,8 +1639,12 @@ void DumpImage::create_image()
         // get pointer to pixmap buffer and get background transparency color
         const auto *pixmap = (const unsigned char *) ubuf(fixarray[i][6]).i;
         double transcolor[3] = {fixarray[i][7], fixarray[i][8], fixarray[i][9]};
-        image->draw_pixmap(&fixarray[i][1], (int) fixarray[i][4], (int) fixarray[i][5], pixmap,
-                           transcolor, fixarray[i][10], opacity);
+        if (ifix.flag1 == 0.0)    // coordinates are in box coordinates
+          image->draw_pixmap(&fixarray[i][1], (int) fixarray[i][4], (int) fixarray[i][5], pixmap,
+                             transcolor, fixarray[i][10], opacity);
+        else    // coordinates are in image coordinates, ignore z
+          image->draw_pixmap((int) fixarray[i][1], (int) fixarray[i][2], (int) fixarray[i][4],
+                             (int) fixarray[i][5], pixmap, transcolor, fixarray[i][10], opacity);
       } else if (fixvec[i] == Graphics::BOND) {
         int type1 = static_cast<int>(fixarray[i][0] - 1.0) % ntypes + 1;
         int type2 = static_cast<int>(fixarray[i][1] - 1.0) % ntypes + 1;
@@ -2320,6 +2324,8 @@ int DumpImage::modify_param(int narg, char **arg)
     if (narg < 3) utils::missing_cmd_args(FLERR, "dump_modify bcolor", error);
     if (atom->nbondtypes == 0)
       error->all(FLERR,"Dump modify bcolor not allowed with no bond types");
+    // ignore if bonds are not displayed
+    if (bondflag == NO) return 3;
     int nlo,nhi;
     utils::bounds_typelabel(FLERR,arg[1],1,atom->nbondtypes,nlo,nhi,lmp,Atom::BOND);
 
@@ -2343,6 +2349,8 @@ int DumpImage::modify_param(int narg, char **arg)
     if (narg < 3) utils::missing_cmd_args(FLERR, "dump_modify bdiam", error);
     if (atom->nbondtypes == 0)
       error->all(FLERR, argoff, "Dump modify bdiam not allowed with no bond types");
+    // ignore if bonds are not displayed
+    if (bondflag == NO) return 3;
     int nlo,nhi;
     utils::bounds_typelabel(FLERR,arg[1],1,atom->nbondtypes,nlo,nhi,lmp,Atom::BOND);
     double diam = utils::numeric(FLERR,arg[2],false,lmp);
@@ -2356,6 +2364,8 @@ int DumpImage::modify_param(int narg, char **arg)
     if (narg < 3) utils::missing_cmd_args(FLERR, "dump_modify btrans", error);
     if (atom->nbondtypes == 0)
       error->all(FLERR,"Dump modify btrans not allowed with no bond types");
+    // ignore if bonds are not displayed
+    if (bondflag == NO) return 3;
     int nlo,nhi;
     utils::bounds_typelabel(FLERR,arg[1],1,atom->nbondtypes,nlo,nhi,lmp,Atom::BOND);
     double opacity = utils::numeric(FLERR,arg[2],false,lmp);
diff --git a/src/GRAPHICS/fix_graphics_labels.cpp b/src/GRAPHICS/fix_graphics_labels.cpp
index 61344730258..bd1b2a42d34 100644
--- a/src/GRAPHICS/fix_graphics_labels.cpp
+++ b/src/GRAPHICS/fix_graphics_labels.cpp
@@ -15,11 +15,14 @@
 
 #include "atom.h"
 #include "comm.h"
+#include "dump_image.h"
 #include "error.h"
 #include "graphics.h"
+#include "image.h"
 #include "input.h"
 #include "memory.h"
 #include "modify.h"
+#include "output.h"
 #include "respa.h"
 #include "text_file_reader.h"
 #include "tokenizer.h"
@@ -80,6 +83,21 @@ void get_color(const std::string &color, unsigned char *rgb)
   }
 }
 
+struct TGAHeader {
+  unsigned char idlength;
+  unsigned char colormaptype;
+  unsigned char datatypecode;
+  unsigned char colormaporigin[2];
+  unsigned char colormaplength[2];
+  unsigned char colormapdepth;
+  unsigned char x_origin[2];
+  unsigned char y_origin[2];
+  unsigned char width[2];
+  unsigned char height[2];
+  unsigned char bitsperpixel;
+  unsigned char imagedescriptor;
+};
+
 // read image into buffer that is locally allocated with new
 // return null pointer if incompatible format or not supported
 
@@ -89,8 +107,8 @@ unsigned char *read_image(FILE *fp, int &width, int &height, const std::string &
   if (!fp) return nullptr;
   unsigned char *pixmap = nullptr;
 
-  if (utils::strmatch(filename, "\\.jpg$") || utils::strmatch(filename, "\\.JPG$") ||
-      utils::strmatch(filename, "\\.jpeg$") || utils::strmatch(filename, "\\.JPEG$")) {
+  if (utils::strmatch(filename, R"(\.jpg$)") || utils::strmatch(filename, R"(\.JPG$)") ||
+      utils::strmatch(filename, R"(\.jpeg$)") || utils::strmatch(filename, R"(\.JPEG$)")) {
 
 #if defined(LAMMPS_JPEG)
     struct jpeg_decompress_struct cinfo;
@@ -129,7 +147,7 @@ unsigned char *read_image(FILE *fp, int &width, int &height, const std::string &
     return nullptr;
 #endif
 
-  } else if (utils::strmatch(filename, "\\.png$") || utils::strmatch(filename, "\\.PNG$")) {
+  } else if (utils::strmatch(filename, R"(\.png$)") || utils::strmatch(filename, R"(\.PNG$)")) {
 
 #if defined(LAMMPS_PNG)
     png_structp png_ptr = nullptr;
@@ -205,7 +223,107 @@ unsigned char *read_image(FILE *fp, int &width, int &height, const std::string &
     info = "PNG image format not supported in this LAMMPS binary";
     return nullptr;
 #endif
+  } else if (utils::strmatch(filename, R"(\.tga$)") || utils::strmatch(filename, R"(\.TGA$)")) {
 
+    TGAHeader header;
+    auto rv = fread(&header, sizeof(TGAHeader), 1, fp);
+    if (rv != 1) {
+      info = "Short TGA file";
+      return nullptr;
+    }
+
+    bool compressed = false;
+    if (header.datatypecode == 10) compressed = true;
+    if ((header.datatypecode != 2) && (header.datatypecode != 10)) {
+      info = "Unsupported TGA file type";
+      return nullptr;
+    }
+    width = (int) header.width[0] + ((int) header.width[1]) * 256;
+    height = (int) header.height[0] + ((int) header.height[1]) * 256;
+
+    bool right2left = (header.imagedescriptor & 0x10) ? true : false;
+    bool fromtop = (header.imagedescriptor & 0x20) ? true : false;
+
+    if (((header.imagedescriptor & 0xC0) != 0) || ((header.imagedescriptor & 0x0F) != 0) ||
+        (header.bitsperpixel != 3 * 8 * sizeof(unsigned char))) {
+      info = "Unsupported TGA file type";
+      return nullptr;
+    }
+    char *id = nullptr;
+    if (header.idlength > 0) {
+      id = new char[header.idlength + 1];
+      if (fread(id, header.idlength, 1, fp) != 1) {
+        delete[] id;
+        return nullptr;
+      }
+      id[header.idlength] = '\0';
+    }
+
+    info = fmt::format("{}x{} TGA file, {}-bit RGB", width, height, (int) header.bitsperpixel / 3);
+    if (right2left) info += ", right-to-left";
+    if (fromtop) info += ", top-to-bottom";
+    if (compressed) info += ", RLE-encoded";
+    if (header.idlength) info += id;
+    delete[] id;
+
+    pixmap = new unsigned char[3 * width * height];
+    if (compressed) {
+      unsigned char len;
+      unsigned char pix[3];
+      int i = 0;
+      while (i < 3 * width * height) {
+        if (0 == fread(&len, 1, 1, fp)) break;
+        if (len < 128) {
+          ++len;
+          for (int j = 0; j < len; ++j) {
+            int y = (fromtop) ? (height - 1 - i / (3 * width)) : i / (3 * width);
+            int x = (right2left) ? (width - 1 - (i - 3 * y * width) / 3) : (i - 3 * y * width) / 3;
+            if (fread(pix, sizeof(unsigned char), 3, fp) != 3) {
+              delete[] pixmap;
+              info = "Short TGA file";
+              return nullptr;
+            }
+            pixmap[y * 3 * width + 3 * x] = pix[2];
+            pixmap[y * 3 * width + 3 * x + 1] = pix[1];
+            pixmap[y * 3 * width + 3 * x + 2] = pix[0];
+            i += 3;
+          }
+        } else {
+          len -= 127;
+          if (fread(pix, sizeof(unsigned char), 3, fp) != 3) {
+            delete[] pixmap;
+            info = "Short TGA file";
+            return nullptr;
+          }
+          for (int j = 0; j < len; ++j) {
+            int y = (fromtop) ? (height - 1 - i / (3 * width)) : i / (3 * width);
+            int x = (right2left) ? (width - 1 - (i - 3 * y * width) / 3) : (i - 3 * y * width) / 3;
+            pixmap[y * 3 * width + 3 * x] = pix[2];
+            pixmap[y * 3 * width + 3 * x + 1] = pix[1];
+            pixmap[y * 3 * width + 3 * x + 2] = pix[0];
+            i += 3;
+          }
+        }
+      }
+    } else {
+      unsigned char pix[3];
+      for (int i = 0; i < height; ++i) {
+        for (int j = 0; j < width; ++j) {
+          int y = fromtop ? (height - 1 - i) : i;
+          int x = right2left ? (width - 1 - j) : j;
+          if (fread(pix, sizeof(unsigned char), 3, fp) != 3) {
+            delete[] pixmap;
+            info = "Short TGA file";
+            return nullptr;
+          }
+          // swap BGR to RGB
+          pixmap[y * 3 * width + 3 * x] = pix[2];
+          pixmap[y * 3 * width + 3 * x + 1] = pix[1];
+          pixmap[y * 3 * width + 3 * x + 2] = pix[0];
+        }
+      }
+    }
+    return pixmap;
   } else {
 
     // read file in NetPBM binary or ASCII format
@@ -341,8 +459,10 @@ FixGraphicsLabels::FixGraphicsLabels(LAMMPS *lmp, int narg, char **arg) :
 
       // check remaining arguments for optional image arguments
       while (iarg < narg) {
-        // next argument is next keyword; exit loop
-        if ((strcmp(arg[iarg], "image") == 0) || (strcmp(arg[iarg], "text") == 0)) break;
+        // if next argument is next keyword; exit loop
+        if ((strcmp(arg[iarg], "image") == 0) || (strcmp(arg[iarg], "text") == 0) ||
+            (strcmp(arg[iarg], "colorscale") == 0))
+          break;
 
         if (strcmp(arg[iarg], "scale") == 0) {
           if (iarg + 2 > narg)
@@ -386,7 +506,7 @@ FixGraphicsLabels::FixGraphicsLabels(LAMMPS *lmp, int narg, char **arg) :
 
       // clang-format off
       TextInfo txt{"", {0.0, 0.0, 0.0}, 0, 0, nullptr, {255, 255, 255}, {192, 192, 192},
-                   {192, 192, 192}, {192, 192, 192}, false, 48.0, 0.5,
+                   {192, 192, 192}, {192, 192, 192}, false, true, 48.0, 0.5,
                    -1, -1, -1, -1, nullptr, nullptr, nullptr, nullptr};
       // clang-format on
 
@@ -400,8 +520,10 @@ FixGraphicsLabels::FixGraphicsLabels(LAMMPS *lmp, int narg, char **arg) :
 
       // check remaining arguments for optional image arguments
       while (iarg < narg) {
-        // next argument is next keyword; exit loop
-        if ((strcmp(arg[iarg], "image") == 0) || (strcmp(arg[iarg], "text") == 0)) break;
+        // if next argument is next keyword; exit loop
+        if ((strcmp(arg[iarg], "image") == 0) || (strcmp(arg[iarg], "text") == 0) ||
+            (strcmp(arg[iarg], "colorscale") == 0))
+          break;
 
         if (strcmp(arg[iarg], "size") == 0) {
           if (iarg + 2 > narg)
@@ -464,17 +586,134 @@ FixGraphicsLabels::FixGraphicsLabels(LAMMPS *lmp, int narg, char **arg) :
             }
           }
           iarg += 2;
+        } else if (strcmp(arg[iarg], "horizontal") == 0) {
+          txt.horizontal = true;
+          ++iarg;
+        } else if (strcmp(arg[iarg], "vertical") == 0) {
+          txt.horizontal = false;
+          ++iarg;
         } else {
           error->all(FLERR, iarg, "Unknown fix graphics/labels text keyword: {}", arg[iarg]);
         }
       }
       texts.emplace_back(txt);
+    } else if (strcmp(arg[iarg], "colorscale") == 0) {
+      if (iarg + 6 > narg) utils::missing_cmd_args(FLERR, "fix graphics/labels colorscale", error);
+
+      // clang-format off
+      ScaleInfo scale{"", "", {0.0, 0.0, 0.0}, 0, 0, nullptr, {255, 255, 255}, {192, 192, 192},
+                      {192, 192, 192}, {192, 192, 192}, false, true, 48.0, 0.5, 0, 0,
+                      -1, -1, -1, -1, nullptr, nullptr, nullptr, nullptr};
+      // clang-format on
+      scale.dumpid = arg[iarg + 1];
+      scale.text = arg[iarg + 2];
+
+      // we always need to trigger computes in case of dynamic color scales
+      varflag = 1;
+
+      PARSE_VARIABLE(scale.pos[0], scale.xstr, iarg + 3);
+      PARSE_VARIABLE(scale.pos[1], scale.ystr, iarg + 4);
+      PARSE_VARIABLE(scale.pos[2], scale.zstr, iarg + 5);
+      iarg += 6;
+
+      // check remaining arguments for optional image arguments
+      while (iarg < narg) {
+        // if next argument is next keyword; exit loop
+        if ((strcmp(arg[iarg], "image") == 0) || (strcmp(arg[iarg], "text") == 0) ||
+            (strcmp(arg[iarg], "colorscale") == 0))
+          break;
+
+        if (strcmp(arg[iarg], "size") == 0) {
+          if (iarg + 2 > narg)
+            utils::missing_cmd_args(FLERR, "fix graphics/labels colorscale size", error);
+          PARSE_VARIABLE(scale.size, scale.sstr, iarg + 1);
+          // for sizes 4 to 64, text is rendered at 2x2 size and scaled down for anti-aliasing.
+          // for larger sizes, the image is rendered at max supported size and scaled as needed.
+          scale.size *= 2.0;
+          if ((scale.size < 8.0) || (scale.size > 1024.0))
+            error->all(FLERR, iarg + 1, "Invalid fix graphics/labels colorscale size value: {}",
+                       scale.size * 0.5);
+          if (scale.size > 128.0) {
+            scale.scale = scale.size / 256.0;
+            scale.size = 128.0;
+          } else {
+            scale.scale = 0.5;
+          }
+          iarg += 2;
+        } else if (strcmp(arg[iarg], "length") == 0) {
+          if (iarg + 2 > narg)
+            utils::missing_cmd_args(FLERR, "fix graphics/labels colorscale length", error);
+          scale.length = 2.0 * utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
+          iarg += 2;
+        } else if (strcmp(arg[iarg], "tics") == 0) {
+          if (iarg + 2 > narg)
+            utils::missing_cmd_args(FLERR, "fix graphics/labels colorscale tics", error);
+          scale.tics = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
+          if (scale.tics < 0) error->all(FLERR, iarg + 1, "Invalid tics value");
+          iarg += 2;
+        } else if (strcmp(arg[iarg], "fontcolor") == 0) {
+          if (iarg + 2 > narg)
+            utils::missing_cmd_args(FLERR, "fix graphics/labels colorscale fontcolor", error);
+          try {
+            get_color(arg[iarg + 1], scale.fontcolor);
+          } catch (TokenizerException &e) {
+            error->all(FLERR, iarg + 1, "Error parsing RGB font color value {}: {}", arg[iarg + 1],
+                       e.what());
+          }
+          iarg += 2;
+        } else if (strcmp(arg[iarg], "backcolor") == 0) {
+          if (iarg + 2 > narg)
+            utils::missing_cmd_args(FLERR, "fix graphics/labels colorscale backcolor", error);
+          try {
+            get_color(arg[iarg + 1], scale.backcolor);
+          } catch (TokenizerException &e) {
+            error->all(FLERR, iarg + 1, "Error parsing RGB font color value {}: {}", arg[iarg + 1],
+                       e.what());
+          }
+          iarg += 2;
+        } else if (strcmp(arg[iarg], "framecolor") == 0) {
+          if (iarg + 2 > narg)
+            utils::missing_cmd_args(FLERR, "fix graphics/labels colorscale framecolor", error);
+          try {
+            get_color(arg[iarg + 1], scale.framecolor);
+          } catch (TokenizerException &e) {
+            error->all(FLERR, iarg + 1, "Error parsing RGB font color value {}: {}", arg[iarg + 1],
+                       e.what());
+          }
+          iarg += 2;
+        } else if (strcmp(arg[iarg], "transcolor") == 0) {
+          if (iarg + 2 > narg)
+            utils::missing_cmd_args(FLERR, "fix graphics/labels colorscale transcolor", error);
+          if (strcmp(arg[iarg + 1], "none") == 0) {
+            scale.notrans = true;
+          } else {
+            try {
+              get_color(arg[iarg + 1], scale.transcolor);
+            } catch (TokenizerException &e) {
+              error->all(FLERR, iarg + 1, "Error parsing RGB font color value {}: {}",
+                         arg[iarg + 1], e.what());
+            }
+          }
+          iarg += 2;
+        } else if (strcmp(arg[iarg], "horizontal") == 0) {
+          scale.horizontal = true;
+          ++iarg;
+        } else if (strcmp(arg[iarg], "vertical") == 0) {
+          scale.horizontal = false;
+          ++iarg;
+        } else {
+          error->all(FLERR, iarg, "Unknown fix graphics/labels colorscale keyword: {}", arg[iarg]);
+        }
+      }
+      scales.emplace_back(scale);
     } else {
       error->all(FLERR, iarg, "Unknown fix graphics/labels keyword: {}", arg[iarg]);
     }
   }
-}
 
+  if (varflag) modify->addstep_compute_all(update->ntimestep);
+}
+#undef PARSE_VARIABLE
 /* ---------------------------------------------------------------------- */
 
 FixGraphicsLabels::~FixGraphicsLabels()
@@ -495,6 +734,14 @@ FixGraphicsLabels::~FixGraphicsLabels()
     delete[] txt.sstr;
   }
 
+  for (auto &scale : scales) {
+    delete[] scale.pixmap;
+    delete[] scale.xstr;
+    delete[] scale.ystr;
+    delete[] scale.zstr;
+    delete[] scale.sstr;
+  }
+
   memory->destroy(imgobjs);
   memory->destroy(imgparms);
 }
@@ -507,96 +754,59 @@ int FixGraphicsLabels::setmask()
 }
 
 /* ---------------------------------------------------------------------- */
+#define CHECK_VARIABLE(index, name)                                                    \
+  if (name) {                                                                          \
+    int ivar = input->variable->find(name);                                            \
+    if (ivar < 0)                                                                      \
+      error->all(FLERR, Error::NOLASTLINE,                                             \
+                 "Variable name {} for fix graphics/labels does not exist", name);     \
+    if (input->variable->equalstyle(ivar) == 0)                                        \
+      error->all(FLERR, Error::NOLASTLINE,                                             \
+                 "Fix graphics/labels variable {} is not equal-style variable", name); \
+    index = ivar;                                                                      \
+  }
 
 void FixGraphicsLabels::init()
 {
   for (auto &pix : pixmaps) {
-    if (pix.xstr) {
-      int ivar = input->variable->find(pix.xstr);
-      if (ivar < 0)
-        error->all(FLERR, Error::NOLASTLINE,
-                   "Variable name {} for fix graphics/labels does not exist", pix.xstr);
-      if (input->variable->equalstyle(ivar) == 0)
-        error->all(FLERR, Error::NOLASTLINE,
-                   "fix graphics/labels variable {} is not equal-style variable", pix.xstr);
-      pix.xvar = ivar;
-    }
-    if (pix.ystr) {
-      int ivar = input->variable->find(pix.ystr);
-      if (ivar < 0)
-        error->all(FLERR, Error::NOLASTLINE,
-                   "Variable name {} for fix graphics/labels does not exist", pix.ystr);
-      if (input->variable->equalstyle(ivar) == 0)
-        error->all(FLERR, Error::NOLASTLINE,
-                   "fix graphics/labels variable {} is not equal-style variable", pix.ystr);
-      pix.yvar = ivar;
-    }
-    if (pix.zstr) {
-      int ivar = input->variable->find(pix.zstr);
-      if (ivar < 0)
-        error->all(FLERR, Error::NOLASTLINE,
-                   "Variable name {} for fix graphics/labels does not exist", pix.zstr);
-      if (input->variable->equalstyle(ivar) == 0)
-        error->all(FLERR, Error::NOLASTLINE,
-                   "fix graphics/labels variable {} is not equal-style variable", pix.zstr);
-      pix.zvar = ivar;
-    }
-    if (pix.sstr) {
-      int ivar = input->variable->find(pix.sstr);
-      if (ivar < 0)
-        error->all(FLERR, Error::NOLASTLINE,
-                   "Variable name {} for fix graphics/labels does not exist", pix.sstr);
-      if (input->variable->equalstyle(ivar) == 0)
-        error->all(FLERR, Error::NOLASTLINE,
-                   "fix graphics/labels variable {} is not equal-style variable", pix.sstr);
-      pix.svar = ivar;
-    }
+    CHECK_VARIABLE(pix.xvar, pix.xstr);
+    CHECK_VARIABLE(pix.yvar, pix.ystr);
+    CHECK_VARIABLE(pix.zvar, pix.zstr);
+    CHECK_VARIABLE(pix.svar, pix.sstr);
   }
 
   for (auto &txt : texts) {
-    if (txt.xstr) {
-      int ivar = input->variable->find(txt.xstr);
-      if (ivar < 0)
-        error->all(FLERR, Error::NOLASTLINE,
-                   "Variable name {} for fix graphics/labels does not exist", txt.xstr);
-      if (input->variable->equalstyle(ivar) == 0)
-        error->all(FLERR, Error::NOLASTLINE,
-                   "fix graphics/labels variable {} is not equal-style variable", txt.xstr);
-      txt.xvar = ivar;
-    }
-    if (txt.ystr) {
-      int ivar = input->variable->find(txt.ystr);
-      if (ivar < 0)
-        error->all(FLERR, Error::NOLASTLINE,
-                   "Variable name {} for fix graphics/labels does not exist", txt.ystr);
-      if (input->variable->equalstyle(ivar) == 0)
-        error->all(FLERR, Error::NOLASTLINE,
-                   "fix graphics/labels variable {} is not equal-style variable", txt.ystr);
-      txt.yvar = ivar;
-    }
-    if (txt.zstr) {
-      int ivar = input->variable->find(txt.zstr);
-      if (ivar < 0)
-        error->all(FLERR, Error::NOLASTLINE,
-                   "Variable name {} for fix graphics/labels does not exist", txt.zstr);
-      if (input->variable->equalstyle(ivar) == 0)
-        error->all(FLERR, Error::NOLASTLINE,
-                   "fix graphics/labels variable {} is not equal-style variable", txt.zstr);
-      txt.zvar = ivar;
-    }
-    if (txt.sstr) {
-      int ivar = input->variable->find(txt.sstr);
-      if (ivar < 0)
-        error->all(FLERR, Error::NOLASTLINE,
-                   "Variable name {} for fix graphics/labels does not exist", txt.sstr);
-      if (input->variable->equalstyle(ivar) == 0)
-        error->all(FLERR, Error::NOLASTLINE,
-                   "fix graphics/labels variable {} is not equal-style variable", txt.sstr);
-      txt.svar = ivar;
-    }
+    CHECK_VARIABLE(txt.xvar, txt.xstr);
+    CHECK_VARIABLE(txt.yvar, txt.ystr);
+    CHECK_VARIABLE(txt.zvar, txt.zstr);
+    CHECK_VARIABLE(txt.svar, txt.sstr);
   }
-}
 
+  for (auto &scale : scales) {
+    CHECK_VARIABLE(scale.xvar, scale.xstr);
+    CHECK_VARIABLE(scale.yvar, scale.ystr);
+    CHECK_VARIABLE(scale.zvar, scale.zstr);
+    CHECK_VARIABLE(scale.svar, scale.sstr);
+
+    // check if dump exists and if the color map is dynamic
+    auto *dump = dynamic_cast<DumpImage *>(output->get_dump_by_id(scale.dumpid));
+    if (!dump)
+      error->all(FLERR, Error::NOLASTLINE,
+                 "Dump ID {} for colorscale not found or not dump style image", scale.dumpid);
+    int dim = 0;
+    auto *image = static_cast<Image *>(dump->extract("image", dim));
+    if (!image || (dim != 0))
+      error->all(FLERR, Error::NOLASTLINE, "Could not extract color scale info from dump {}",
+                 scale.dumpid);
+    double lo, hi;
+    if (image->map_info(0, lo, hi) && (comm->me == 0))
+      error->warning(FLERR,
+                     "Dump {} uses a dynamic color map. "
+                     "Color scale can only use data from previous dump output\n",
+                     scale.dumpid);
+  }
+}
+#undef CHECK_VARIABLE
 /* ---------------------------------------------------------------------- */
 
 void FixGraphicsLabels::setup(int)
@@ -608,7 +818,7 @@ void FixGraphicsLabels::setup(int)
 
 void FixGraphicsLabels::end_of_step()
 {
-  numobjs = pixmaps.size() + texts.size();
+  numobjs = pixmaps.size() + texts.size() + scales.size();
   if (numobjs == 0) return;
 
   if (varflag) modify->clearstep_compute();
@@ -620,6 +830,9 @@ void FixGraphicsLabels::end_of_step()
 
   int n = 0;
   for (auto &pix : pixmaps) {
+
+    // update values from variables
+
     if (pix.xstr) pix.pos[0] = input->variable->compute_equal(pix.xvar);
     if (pix.ystr) pix.pos[1] = input->variable->compute_equal(pix.yvar);
     if (pix.zstr) pix.pos[2] = input->variable->compute_equal(pix.zvar);
@@ -664,10 +877,16 @@ void FixGraphicsLabels::end_of_step()
 
   // initialize font renderer and load in-memory font
 
+  SSFN::ScalableFont renderfont;
+
   try {
-    SSFN::ScalableFont renderfont;
+
+    // process text labels
 
     for (auto &txt : texts) {
+
+      // update values from variables
+
       if (txt.xstr) txt.pos[0] = input->variable->compute_equal(txt.xvar);
       if (txt.ystr) txt.pos[1] = input->variable->compute_equal(txt.yvar);
       if (txt.zstr) txt.pos[2] = input->variable->compute_equal(txt.zvar);
@@ -704,8 +923,8 @@ void FixGraphicsLabels::end_of_step()
         }
 
         delete[] txt.pixmap;
-        txt.pixmap = renderfont.create_pixmap(expanded, txt.width, txt.height, txt.fontcolor,
-                                              txt.framecolor, txt.backcolor);
+        txt.pixmap = renderfont.create_label(expanded, txt.width, txt.height, txt.fontcolor,
+                                             txt.framecolor, txt.backcolor, txt.horizontal);
       }
       imgobjs[n] = Graphics::PIXMAP;
       imgparms[n][0] = 1;
@@ -730,6 +949,82 @@ void FixGraphicsLabels::end_of_step()
     error->all(FLERR, Error::NOLASTLINE, "Error during font rendering: {}", e.what());
   }
 
+  // process color scales
+  try {
+    for (auto &scale : scales) {
+
+      auto *dump = dynamic_cast<DumpImage *>(output->get_dump_by_id(scale.dumpid));
+      if (!dump)
+        error->all(FLERR, "Dump ID {} for colorscale not found or not dump style image",
+                   scale.dumpid);
+      int dim = 0;
+      auto *image = static_cast<Image *>(dump->extract("image", dim));
+      if (!image || (dim != 0))
+        error->all(FLERR, "Could not extract color scale info from dump {}", scale.dumpid);
+
+      // update values from variables
+
+      if (scale.xstr) scale.pos[0] = input->variable->compute_equal(scale.xvar);
+      if (scale.ystr) scale.pos[1] = input->variable->compute_equal(scale.yvar);
+      if (scale.zstr) scale.pos[2] = input->variable->compute_equal(scale.zvar);
+
+      // text is rasterized at twice the size for some anti-aliasing. clamp to avoid crashes.
+      if (scale.sstr) {
+        scale.size = 2.0 * input->variable->compute_equal(scale.svar);
+        if (scale.size > 128.0) {
+          scale.scale = scale.size / 256.0;
+          scale.size = 128.0;
+        } else {
+          scale.size = MAX(scale.size, 8.0);
+          scale.scale = 0.5;
+        }
+      }
+
+      renderfont.select_font(SSFN::FAMILY_SANS, SSFN::STYLE_REGULAR, (int) scale.size);
+
+      auto expanded = scale.text;
+
+      // substitute variables in text
+      if (expanded.find('$') != std::string::npos) {
+        int ncopy = expanded.length() + 1;
+        int nwork = ncopy;
+        char *copy = (char *) memory->smalloc(ncopy * sizeof(char), "fix/graphics/labels:copy");
+        char *work = (char *) memory->smalloc(nwork * sizeof(char), "fix/graphics/labels:work");
+        strncpy(copy, expanded.c_str(), ncopy);
+        input->substitute(copy, work, ncopy, nwork, 0);
+        expanded = copy;
+        memory->sfree(copy);
+        memory->sfree(work);
+      }
+
+      delete[] scale.pixmap;
+      scale.pixmap = renderfont.create_colorscale(
+          expanded, scale.width, scale.height, scale.fontcolor, scale.framecolor, scale.backcolor,
+          scale.horizontal, scale.length, image, 0, scale.tics);
+
+      imgobjs[n] = Graphics::PIXMAP;
+      imgparms[n][0] = 1;
+      imgparms[n][1] = scale.pos[0];
+      imgparms[n][2] = scale.pos[1];
+      imgparms[n][3] = scale.pos[2];
+      imgparms[n][4] = scale.width;
+      imgparms[n][5] = scale.height;
+      imgparms[n][6] = ubuf((int64_t) scale.pixmap).d;
+      if (scale.notrans) {
+        imgparms[n][7] = imgparms[n][8] = imgparms[n][9] = -1.0;
+      } else {
+        imgparms[n][7] = (double) scale.transcolor[0] / 255.0;
+        imgparms[n][8] = (double) scale.transcolor[1] / 255.0;
+        imgparms[n][9] = (double) scale.transcolor[2] / 255.0;
+      }
+
+      imgparms[n][10] = scale.scale;
+      ++n;
+    }
+  } catch (const SSFN::SSFNException &e) {
+    error->all(FLERR, Error::NOLASTLINE, "Error during font rendering: {}", e.what());
+  }
+
   if (varflag) modify->addstep_compute((update->ntimestep / nevery) * nevery + nevery);
 }
 
diff --git a/src/GRAPHICS/fix_graphics_labels.h b/src/GRAPHICS/fix_graphics_labels.h
index 6143367fbbe..064548e149a 100644
--- a/src/GRAPHICS/fix_graphics_labels.h
+++ b/src/GRAPHICS/fix_graphics_labels.h
@@ -63,6 +63,7 @@ class FixGraphicsLabels : public Fix {
     unsigned char framecolor[3];
     unsigned char transcolor[3];
     bool notrans;
+    bool horizontal;
     double size;
     double scale;
     int xvar, yvar, zvar, svar;
@@ -70,6 +71,28 @@ class FixGraphicsLabels : public Fix {
   };
   std::vector<TextInfo> texts;
 
+  struct ScaleInfo {
+    std::string dumpid;
+    std::string text;
+    double pos[3];
+    int width;
+    int height;
+    unsigned char *pixmap;
+    unsigned char fontcolor[3];
+    unsigned char backcolor[3];
+    unsigned char framecolor[3];
+    unsigned char transcolor[3];
+    bool notrans;
+    bool horizontal;
+    double size;
+    double scale;
+    int length;
+    int tics;
+    int xvar, yvar, zvar, svar;
+    char *xstr, *ystr, *zstr, *sstr;
+  };
+  std::vector<ScaleInfo> scales;
+
   int varflag;
   int numobjs;
   int *imgobjs;
diff --git a/src/GRAPHICS/image.cpp b/src/GRAPHICS/image.cpp
index 4455cbcd941..9ab3e584900 100644
--- a/src/GRAPHICS/image.cpp
+++ b/src/GRAPHICS/image.cpp
@@ -13,7 +13,8 @@
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
-   Contributing author: Nathan Fabian (Sandia)
+   Contributing authors: Nathan Fabian (Sandia)
+                     and Axel Kohlmeyer (Temple)
 ------------------------------------------------------------------------- */
 
 #include "image.h"
@@ -825,11 +826,34 @@ void Image::draw_pixmap(const double *x, int pixwidth, int pixheight, const unsi
   xc += width / 2;
   yc += height / 2;
 
+  // convert back to non-FSAA image coordinates, so we can re-use the pixmap drawing code
+  if (fsaa) {
+    xc /= 2;
+    yc /= 2;
+  }
+
+  draw_pixmap(xc, yc, pixwidth, pixheight, pixmap, transcolor, scale, opacity, dist);
+}
+
+/* ----------------------------------------------------------------------
+   scale and add pixmap centered at location xc, yc in image coordinates to image
+   background color indicates transparency and pixels in that color are skipped
+------------------------------------------------------------------------- */
+
+void Image::draw_pixmap(int xc, int yc, int pixwidth, int pixheight, const unsigned char *pixmap,
+                        double *transcolor, double scale, double opacity, double dist)
+{
   const unsigned char *mypixmap = pixmap;
   unsigned char *npixmap = nullptr;
 
-  // adjust scale factor for FSAA and only scale as much as needed.
-  if (fsaa) scale *= 2.0;
+  // adjust scale factor and image location for FSAA
+  if (fsaa) {
+    scale *= 2.0;
+    xc *= 2;
+    yc *= 2;
+  }
+
+  // only scale as much as needed.
   if (scale != 1.0) {
     int nwidth = std::lround(scale * pixwidth + 0.5);
     int nheight = std::lround(scale * pixheight + 0.5);
@@ -1729,6 +1753,15 @@ int Image::map_minmax(int index, double mindynamic, double maxdynamic)
   return maps[index]->minmax(mindynamic,maxdynamic);
 }
 
+/* ----------------------------------------------------------------------
+   get min/max bounds of dynamic color map index and return 1 if dynamic
+------------------------------------------------------------------------- */
+
+int Image::map_info(int index, double &min, double &max)
+{
+  return maps[index]->info(min, max);
+}
+
 /* ----------------------------------------------------------------------
    return 3-vector color corresponding to value from color map index
 ------------------------------------------------------------------------- */
@@ -2290,8 +2323,8 @@ ColorMap::ColorMap(LAMMPS *lmp, Image *caller) : Pointers(lmp)
 
   dynamic = 1;
 
-  mlo = MINVALUE;
-  mhi = MAXVALUE;
+  locurrent = mlo = MINVALUE;
+  hicurrent = mhi = MAXVALUE;
   mstyle = CONTINUOUS;
   mrange = FRACTIONAL;
 
@@ -2355,6 +2388,7 @@ int ColorMap::reset(int narg, char **arg)
   if (nentry < 1) return 5;
   delete [] mentry;
   mentry = new MapEntry[nentry];
+  mentry[0].svalue = 0.0;
 
   int n = 5;
   for (int i = 0; i < nentry; i++) {
@@ -2467,6 +2501,13 @@ int ColorMap::minmax(double mindynamic, double maxdynamic)
   return 0;
 }
 
+int ColorMap::info(double &min, double &max)
+{
+  min = locurrent;
+  max = hicurrent;
+  return dynamic;
+}
+
 /* ----------------------------------------------------------------------
    convert value into an RGB color via color map
    return pointer to 3-vector
diff --git a/src/GRAPHICS/image.h b/src/GRAPHICS/image.h
index 02b3382d7f5..ac87d2b5cb1 100644
--- a/src/GRAPHICS/image.h
+++ b/src/GRAPHICS/image.h
@@ -61,10 +61,13 @@ class Image : protected Pointers {
   void draw_axes(double (*)[3], double, double opacity = 1.0);
   void draw_pixmap(const double *, int, int, const unsigned char *, double *, double scale = 1.0,
                    double opacity = 1.0);
+  void draw_pixmap(int, int, int, int, const unsigned char *, double *, double scale = 1.0,
+                   double opacity = 1.0, double depth = 0.0);
 
   int map_dynamic(int);
   int map_reset(int, int, char **);
   int map_minmax(int, double, double);
+  int map_info(int, double &, double &);
   double *map_value2color(int, double);
 
   int addcolor(char *, double, double, double);
@@ -164,6 +167,7 @@ class ColorMap : protected Pointers {
   ~ColorMap() override;
   int reset(int, char **);
   int minmax(double, double);
+  int info(double &, double &);
   double *value2color(double);
 
  private:
diff --git a/src/GRAPHICS/scalable_font.cpp b/src/GRAPHICS/scalable_font.cpp
index 9fc52f083e6..d00f6c22293 100644
--- a/src/GRAPHICS/scalable_font.cpp
+++ b/src/GRAPHICS/scalable_font.cpp
@@ -45,6 +45,7 @@
 #include <cstring>
 
 #include "fmt/format.h"
+#include "image.h"
 #include "scalable_font.h"
 
 static constexpr int SSFN_DATA_MAX = 65536;
@@ -156,16 +157,20 @@ using ssfn_t = struct _ssfn_t {
 /*** normal renderer (ca. 22k, fully featured with error checking) ***/
 
 namespace {
-/* error codes */
-#define SSFN_OK 0           /* success */
-#define SSFN_ERR_ALLOC 1    /* allocation error */
-#define SSFN_ERR_NOFACE 2   /* no font face selected */
-#define SSFN_ERR_INVINP 3   /* invalid input */
-#define SSFN_ERR_BADFILE 4  /* bad SSFN file format */
-#define SSFN_ERR_BADSTYLE 5 /* bad style */
-#define SSFN_ERR_BADSIZE 6  /* bad size */
-#define SSFN_ERR_BADMODE 7  /* bad mode */
-#define SSFN_ERR_NOGLYPH 8  /* glyph (or kerning info) not found */
+  /* error codes */
+  enum {
+    SSFN_OK = 0,       /* success */
+    SSFN_ERR_ALLOC,    /* allocation error */
+    SSFN_ERR_NOFACE,   /* no font face selected */
+    SSFN_ERR_INVINP,   /* invalid input */
+    SSFN_ERR_BADFILE,  /* bad SSFN file format */
+    SSFN_ERR_BADSTYLE, /* bad style */
+    SSFN_ERR_BADSIZE,  /* bad size */
+    SSFN_ERR_BADMODE,  /* bad mode */
+    SSFN_ERR_NOGLYPH,  /* glyph (or kerning info) not found */
+    SSFN_ERR_NOMAP,    /* could not access color map information */
+    SSFN_ERR_LAST
+  };
 
   /**
    * Error code strings
@@ -178,7 +183,8 @@ namespace {
                                "Invalid style",
                                "Invalid size",
                                "Invalid mode",
-                               "Glyph not found"};
+                               "Glyph not found",
+                               "Color map not available"};
 
 // include font data as constant in memory byte sequence
 #include "scalable_sans_font.h"
@@ -1055,9 +1061,9 @@ void ScalableFont::select_font(int family, int style, int size)
   _ssfn_select((ssfn_t *) ctx, family, style, size);
 }
 
-unsigned char *ScalableFont::create_pixmap(const std::string &text, int &width, int &height,
-                                           const unsigned char *font, const unsigned char *frame,
-                                           const unsigned char *back)
+unsigned char *ScalableFont::create_label(const std::string &text, int &width, int &height,
+                                          const unsigned char *font, const unsigned char *frame,
+                                          const unsigned char *back, bool horizontal)
 {
   auto *ctxptr = (ssfn_t *) ctx;
   ssfn_glyph_t *g;
@@ -1140,13 +1146,222 @@ unsigned char *ScalableFont::create_pixmap(const std::string &text, int &width,
     penx += g->adv_x;
     free(g);
   }
+
+  // for vertical text copy the rotated pixmap accordingly and swap width and height
+  if (!horizontal) {
+    auto *flipped = new unsigned char[width * height * 3];
+    for (int y = 0; y < height; ++y) {
+      for (int x = 0; x < width; ++x) {
+        for (int i = 0; i < 3; ++i) {
+          flipped[x * 3 * height + 3 * (height - 1 - y) + i] = pixmap[y * 3 * width + 3 * x + i];
+        }
+      }
+    }
+    delete[] pixmap;
+    pixmap = flipped;
+    std::swap(width, height);
+  }
+  return pixmap;
+}
+
+unsigned char *ScalableFont::create_colorscale(const std::string &text, int &width, int &height,
+                                               const unsigned char *font,
+                                               const unsigned char *frame,
+                                               const unsigned char *back, bool horizontal,
+                                               int minwidth, LAMMPS_NS::Image *image, int mapidx,
+                                               int tics)
+{
+  auto *ctxptr = (ssfn_t *) ctx;
+  ssfn_glyph_t *g;
+
+  // get a font size specific spacing for a border
+  g = _ssfn_render(ctxptr, ' ');
+  int xspace = g->adv_x;
+  free(g);
+  // get a font size specific spacing for a angular bracket
+  g = _ssfn_render(ctxptr, '<');
+  int xfill = g->adv_x;
+  free(g);
+
+  double lo, hi;
+  image->map_info(mapidx, lo, hi);
+  auto newtext = fmt::format("{:3.3}  {}  {:3.3}", lo, text, hi);
+
+  // dry run to determine size of pixmap
+  width = 0;
+  int miny = 1073741824;
+  int maxy = 0;
+  for (auto c : newtext + "gll") {    // append these characters for consistent spacing
+    if (c == '_') c = ' ';            // ugly hack to work around font issue
+
+    // render character and apply its width
+    g = _ssfn_render(ctxptr, c);
+    width += g->adv_x;
+
+    // loop over bitmap to find minimum and maximum y position
+    for (int y = 0; y < g->h; ++y) {
+      const int ypos = g->h - 1 - y + g->baseline;
+      for (int x = 0, i = 0, m = 1; x < g->w; ++x, m <<= 1) {
+        if (m > 0x80) {
+          m = 1;
+          ++i;
+        }
+        if (g->data[y * g->pitch + i] & m) {
+          miny = std::min(miny, ypos);
+          maxy = std::max(maxy, ypos);
+        }
+      }
+    }
+    free(g);
+  }
+
+  int xhalf = xspace / 2;
+  height = maxy - miny + 1 + 6 * xspace;
+
+  if (minwidth > width) {
+    int wextra = (minwidth - width) / xfill / 4;
+    newtext =
+        fmt::format("{:3.3} {:<<{}} {} {:>>{}} {:3.3}", lo, '<', wextra, text, '>', wextra, hi);
+    width = minwidth;
+  }
+
+  // re-do dry run to determine size of the pixmap with the padded range string
+  width = 0;
+  miny = 1073741824;
+  maxy = 0;
+  for (auto c : newtext + "gll") {    // append these characters for consistent spacing
+    if (c == '_') c = ' ';            // ugly hack to work around font issue
+
+    // render character and apply its width
+    g = _ssfn_render(ctxptr, c);
+    width += g->adv_x;
+
+    // loop over bitmap to find minimum and maximum y position
+    for (int y = 0; y < g->h; ++y) {
+      const int ypos = g->h - 1 - y + g->baseline;
+      for (int x = 0, i = 0, m = 1; x < g->w; ++x, m <<= 1) {
+        if (m > 0x80) {
+          m = 1;
+          ++i;
+        }
+        if (g->data[y * g->pitch + i] & m) {
+          miny = std::min(miny, ypos);
+          maxy = std::max(maxy, ypos);
+        }
+      }
+    }
+    free(g);
+  }
+
+  // allocate and fill pixmap with background and frame color
+
+  auto *pixmap = new unsigned char[width * height * 3];
+  for (int y = 0; y < height; ++y) {
+    int yoffs = 3 * y * width;
+    for (int x = 0; x < width; ++x) {
+      if ((y < xhalf) || (y >= height - xhalf) || (x < xhalf) || (x >= width - xhalf)) {
+        pixmap[yoffs + 3 * x] = frame[0];
+        pixmap[yoffs + 3 * x + 1] = frame[1];
+        pixmap[yoffs + 3 * x + 2] = frame[2];
+      } else {
+        pixmap[yoffs + 3 * x] = back[0];
+        pixmap[yoffs + 3 * x + 1] = back[1];
+        pixmap[yoffs + 3 * x + 2] = back[2];
+      }
+    }
+  }
+
+  // draw colormap
+  double delta = (hi - lo) / (width - 2 * xspace);
+  int ticinc = tics ? (width - 2 * xspace) / (tics + 1) : 1 << 31;
+  int ticmin = ticinc + xspace + 1;
+  int ticmax = ticmin + xhalf;
+  for (int x = xspace; x < width - xspace; ++x) {
+    for (int y = xhalf + xhalf / 2; y < height - 4 * xspace - xhalf; ++y) {
+      int offs = 3 * y * width + 3 * x;
+      double val = lo + delta * static_cast<double>(x - xspace);
+      auto *color = image->map_value2color(mapidx, val);
+      if (color) {
+        if (tics > 0) {
+          if ((((x >= ticmin) && (x <= ticmax)) || (x <= xspace + xhalf) ||
+               (x >= width - xspace - xhalf)) ||
+              (y <= xspace + xhalf / 2) || (y >= height - 5 * xspace)) {
+
+            pixmap[offs] = font[0];
+            pixmap[offs + 1] = font[1];
+            pixmap[offs + 2] = font[2];
+          } else {
+            pixmap[offs] = color[0] * 255;
+            pixmap[offs + 1] = color[1] * 255;
+            pixmap[offs + 2] = color[2] * 255;
+          }
+          if (x == ticmax) {
+            ticmin += ticinc;
+            ticmax += ticinc;
+          }
+        } else {
+          if ((x <= xspace + xhalf) || (x >= width - xspace - xhalf) || (y <= xspace + xhalf / 2) ||
+              (y >= height - 5 * xspace)) {
+            pixmap[offs] = font[0];
+            pixmap[offs + 1] = font[1];
+            pixmap[offs + 2] = font[2];
+          } else {
+            pixmap[offs] = color[0] * 255;
+            pixmap[offs + 1] = color[1] * 255;
+            pixmap[offs + 2] = color[2] * 255;
+          }
+        }
+      }
+    }
+  }
+
+  // now render each character again and change the pixels in the pixmap accordingly
+  int penx = 2 * xspace;
+  int peny = 3 * xspace + xhalf;
+  for (auto c : newtext) {
+    if (c == '_') c = ' ';    // ugly hack to work around font issue
+
+    g = _ssfn_render(ctxptr, c);
+    for (int y = 0; y < g->h; ++y) {
+      const int yoffs = (g->h - 1 - y + peny + g->baseline - miny + xspace + xhalf / 2) * width * 3;
+      for (int x = 0, i = 0, m = 1; x < g->w; ++x, m <<= 1) {
+        if (m > 0x80) {
+          m = 1;
+          ++i;
+        }
+        const int xoffs = (penx + x) * 3;
+        if (g->data[y * g->pitch + i] & m) {
+          pixmap[yoffs + xoffs] = font[0];
+          pixmap[yoffs + xoffs + 1] = font[1];
+          pixmap[yoffs + xoffs + 2] = font[2];
+        }
+      }
+    }
+    penx += g->adv_x;
+    free(g);
+  }
+
+  // for vertical text copy the rotated pixmap accordingly and swap width and height
+  if (!horizontal) {
+    auto *flipped = new unsigned char[width * height * 3];
+    for (int y = 0; y < height; ++y) {
+      for (int x = 0; x < width; ++x) {
+        for (int i = 0; i < 3; ++i) {
+          flipped[x * 3 * height + 3 * (height - 1 - y) + i] = pixmap[y * 3 * width + 3 * x + i];
+        }
+      }
+    }
+    delete[] pixmap;
+    pixmap = flipped;
+    std::swap(width, height);
+  }
   return pixmap;
 }
 
 SSFNException::SSFNException(const std::string &file, int line, int flag)
 {
   message = fmt::format("In file {}:{} ", truncpath(file), line);
-  if ((flag < SSFN_OK) || (flag > SSFN_ERR_NOGLYPH)) {
+  if ((flag < SSFN_OK) || (flag >= SSFN_ERR_LAST)) {
     message.append("Unknown Error");
   } else {
     message.append(ssfn_errstr[flag]);
diff --git a/src/GRAPHICS/scalable_font.h b/src/GRAPHICS/scalable_font.h
index ef09e0d8965..9efc398daa1 100644
--- a/src/GRAPHICS/scalable_font.h
+++ b/src/GRAPHICS/scalable_font.h
@@ -35,9 +35,13 @@ class ScalableFont {
   ~ScalableFont();
 
   void select_font(int family, int style, int size);
-  unsigned char *create_pixmap(const std::string &text, int &width, int &height,
-                               const unsigned char *font, const unsigned char *frame,
-                               const unsigned char *back);
+  unsigned char *create_label(const std::string &text, int &width, int &height,
+                              const unsigned char *font, const unsigned char *frame,
+                              const unsigned char *back, bool horizontal = true);
+  unsigned char *create_colorscale(const std::string &text, int &width, int &height,
+                                   const unsigned char *font, const unsigned char *frame,
+                                   const unsigned char *back, bool horizontal, int minwidth,
+                                   class LAMMPS_NS::Image *image, int mapidx, int tics);
 
  private:
   void *ctx;
diff --git a/src/KOKKOS/atom_kokkos.cpp b/src/KOKKOS/atom_kokkos.cpp
index d9d85e2a3e7..d77f8cc5a55 100644
--- a/src/KOKKOS/atom_kokkos.cpp
+++ b/src/KOKKOS/atom_kokkos.cpp
@@ -214,13 +214,21 @@ void AtomKokkos::sort()
       if (!fix_iextra->sort_device) {
         flag = 0;
         if (comm->me == 0)
-          error->warning(FLERR,"Fix {} not compatible with Kokkos sorting on device", fix_iextra->style);
+          error->warning(FLERR,"Fix {} not (yet) compatible with Kokkos sorting on device", fix_iextra->style);
         break;
       }
     }
     if (!flag) {
       if (comm->me == 0) {
-        error->warning(FLERR,"Fix with atom-based arrays not compatible with Kokkos sorting on device, "
+        error->warning(FLERR,"Fix with atom-based arrays not (yet) compatible with Kokkos sorting on device, "
+                           "switching to legacy host sorting");
+      }
+      sort_legacy = true;
+    }
+
+    if (hybrid_flag) {
+      if (comm->me == 0) {
+        error->warning(FLERR,"Atom_style hybrid not (yet) compatible with Kokkos sorting on device, "
                            "switching to legacy host sorting");
       }
       sort_legacy = true;
@@ -415,12 +423,12 @@ AtomVec *AtomKokkos::new_avec(const std::string &style, int trysuffix, int &sfla
 {
   // check if avec already exists, if so this is a hybrid substyle
 
-  int hybrid_substyle_flag = (avec != nullptr);
+  hybrid_flag = (avec != nullptr);
 
   AtomVec *avec = Atom::new_avec(style, trysuffix, sflag);
   if (!avec->kokkosable) error->all(FLERR, "KOKKOS package requires a Kokkos-enabled atom_style");
 
-  if (!hybrid_substyle_flag)
+  if (!hybrid_flag)
     avecKK = dynamic_cast<AtomVecKokkos*>(avec);
 
   return avec;
diff --git a/src/KOKKOS/atom_kokkos.h b/src/KOKKOS/atom_kokkos.h
index 07aa1df57e4..80ec4f7e78b 100644
--- a/src/KOKKOS/atom_kokkos.h
+++ b/src/KOKKOS/atom_kokkos.h
@@ -25,7 +25,7 @@ namespace LAMMPS_NS {
 class AtomKokkos : public Atom {
  public:
   bool sort_legacy;
-  int nprop_atom;
+  int nprop_atom, hybrid_flag;
   class FixPropertyAtomKokkos **fix_prop_atom;
 
   DAT::tdual_tagint_1d k_tag;
diff --git a/src/KOKKOS/comm_kokkos.cpp b/src/KOKKOS/comm_kokkos.cpp
index 5163283fa88..6d14e6099da 100644
--- a/src/KOKKOS/comm_kokkos.cpp
+++ b/src/KOKKOS/comm_kokkos.cpp
@@ -174,6 +174,7 @@ void CommKokkos::forward_comm_device()
           if (size_forward_recv[iswap]) {
             buf = (double*)atomKK->k_x.view<DeviceType>().data() +
               firstrecv[iswap]*atomKK->k_x.view<DeviceType>().extent(1);
+            DeviceType().fence();
             MPI_Irecv(buf,size_forward_recv[iswap],MPI_DOUBLE,
                       recvproc[iswap],0,world,&request);
           }
@@ -193,6 +194,7 @@ void CommKokkos::forward_comm_device()
 
         } else if (ghost_velocity) {
           if (size_forward_recv[iswap]) {
+            DeviceType().fence();
             MPI_Irecv(k_buf_recv.view<DeviceType>().data(),
                       size_forward_recv[iswap],MPI_DOUBLE,
                       recvproc[iswap],0,world,&request);
@@ -211,10 +213,12 @@ void CommKokkos::forward_comm_device()
           }
           atomKK->avecKK->unpack_comm_vel_kokkos(recvnum[iswap],firstrecv[iswap],k_buf_recv);
         } else {
-          if (size_forward_recv[iswap])
+          if (size_forward_recv[iswap]) {
+            DeviceType().fence();
             MPI_Irecv(k_buf_recv.view<DeviceType>().data(),
                       size_forward_recv[iswap],MPI_DOUBLE,
                       recvproc[iswap],0,world,&request);
+          }
           auto k_sendlist_iswap = Kokkos::subview(k_sendlist,iswap,Kokkos::ALL);
           n = atomKK->avecKK->pack_comm_kokkos(sendnum[iswap],k_sendlist_iswap,
                                      k_buf_send,pbc_flag[iswap],pbc[iswap]);
@@ -285,9 +289,11 @@ void CommKokkos::reverse_comm_device()
   for (int iswap = nswap-1; iswap >= 0; iswap--) {
     if (sendproc[iswap] != me) {
       if (comm_f_only && !atomKK->k_f.NEED_TRANSFORM) {
-        if (size_reverse_recv[iswap])
-            MPI_Irecv(k_buf_recv.view<DeviceType>().data(),size_reverse_recv[iswap],MPI_DOUBLE,
+        if (size_reverse_recv[iswap]) {
+          DeviceType().fence();
+          MPI_Irecv(k_buf_recv.view<DeviceType>().data(),size_reverse_recv[iswap],MPI_DOUBLE,
                     sendproc[iswap],0,world,&request);
+        }
         if (size_reverse_send[iswap]) {
           buf = (double*)atomKK->k_f.view<DeviceType>().data() +
             firstrecv[iswap]*atomKK->k_f.view<DeviceType>().extent(1);
@@ -302,10 +308,12 @@ void CommKokkos::reverse_comm_device()
         }
 
       } else {
-        if (size_reverse_recv[iswap])
+        if (size_reverse_recv[iswap]) {
+          DeviceType().fence();
           MPI_Irecv(k_buf_recv.view<DeviceType>().data(),
                     size_reverse_recv[iswap],MPI_DOUBLE,
                     sendproc[iswap],0,world,&request);
+        }
         n = atomKK->avecKK->pack_reverse_kokkos(recvnum[iswap],firstrecv[iswap],k_buf_send);
         if (n) {
           DeviceType().fence();
@@ -395,6 +403,7 @@ void CommKokkos::forward_comm_device(Fix *fix, int size)
       }
 
       if (recvnum[iswap]) {
+        DeviceType().fence();
         MPI_Irecv(buf_recv_fix,nsize*recvnum[iswap],MPI_DOUBLE,
                   recvproc[iswap],0,world,&request);
       }
@@ -576,6 +585,7 @@ void CommKokkos::forward_comm_device(Pair *pair, int size)
       }
 
       if (recvnum[iswap]) {
+        DeviceType().fence();
         MPI_Irecv(buf_recv_pair,nsize*recvnum[iswap],MPI_DOUBLE,
                   recvproc[iswap],0,world,&request);
       }
@@ -675,8 +685,10 @@ void CommKokkos::reverse_comm_device(Pair *pair, int size)
     }
 
     if (sendproc[iswap] != me) {
-      if (sendnum[iswap])
+      if (sendnum[iswap]) {
+        DeviceType().fence();
         MPI_Irecv(buf_recv_pair,nsize*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0,world,&request);
+      }
       if (recvnum[iswap]) {
         DeviceType().fence();
         MPI_Send(buf_send_pair,n,MPI_DOUBLE,recvproc[iswap],0,world);
@@ -934,6 +946,7 @@ void CommKokkos::exchange_device()
         }
         if (nrecv > maxrecv) grow_recv_kokkos(nrecv);
 
+        DeviceType().fence();
         MPI_Irecv(k_buf_recv.view<DeviceType>().data(),nrecv1,
                   MPI_DOUBLE,procneigh[dim][1],0,
                   world,&request);
@@ -1003,6 +1016,7 @@ void CommKokkos::exchange_device()
 
             if (nextrarecv > maxrecv) grow_recv_kokkos(nextrarecv);
 
+            DeviceType().fence();
             MPI_Irecv(k_buf_recv.view<DeviceType>().data(),nextrarecv1,
                       MPI_DOUBLE,procneigh[dim][1],0,
                       world,&request);
@@ -1313,9 +1327,12 @@ void CommKokkos::borders_device() {
         MPI_Sendrecv(&nsend,1,MPI_INT,sendproc[iswap],0,
                      &nrecv,1,MPI_INT,recvproc[iswap],0,world,MPI_STATUS_IGNORE);
         if (nrecv*size_border > maxrecv) grow_recv_kokkos(nrecv*size_border);
-        if (nrecv) MPI_Irecv(k_buf_recv.view<DeviceType>().data(),
-                             nrecv*size_border,MPI_DOUBLE,
-                             recvproc[iswap],0,world,&request);
+        if (nrecv) {
+          DeviceType().fence();
+          MPI_Irecv(k_buf_recv.view<DeviceType>().data(),
+                    nrecv*size_border,MPI_DOUBLE,
+                    recvproc[iswap],0,world,&request);
+        }
         if (n) {
           DeviceType().fence();
           MPI_Send(k_buf_send.view<DeviceType>().data(),n,
diff --git a/src/KOKKOS/comm_tiled_kokkos.cpp b/src/KOKKOS/comm_tiled_kokkos.cpp
index f05fadfae4c..ad4fb343cda 100644
--- a/src/KOKKOS/comm_tiled_kokkos.cpp
+++ b/src/KOKKOS/comm_tiled_kokkos.cpp
@@ -138,6 +138,7 @@ void CommTiledKokkos::forward_comm_device()
         for (i = 0; i < nrecv; i++) {
           buf = (double*)atomKK->k_x.view<DeviceType>().data() +
             firstrecv[iswap][i]*atomKK->k_x.view<DeviceType>().extent(1);
+          DeviceType().fence();
           MPI_Irecv(buf,size_forward_recv[iswap][i],
                     MPI_DOUBLE,recvproc[iswap][i],0,world,&requests[i]);
         }
@@ -166,6 +167,7 @@ void CommTiledKokkos::forward_comm_device()
         for (i = 0; i < nrecv; i++) {
           buf = k_buf_recv.view<DeviceType>().data() +
             forward_recv_offset[iswap][i]*k_buf_recv.view<DeviceType>().extent(1);
+          DeviceType().fence();
           MPI_Irecv(buf,
                     size_forward_recv[iswap][i],MPI_DOUBLE,recvproc[iswap][i],0,world,&requests[i]);
         }
@@ -201,6 +203,7 @@ void CommTiledKokkos::forward_comm_device()
         for (i = 0; i < nrecv; i++) {
           buf = k_buf_recv.view<DeviceType>().data() +
             forward_recv_offset[iswap][i]*k_buf_recv.view<DeviceType>().extent(1);
+          DeviceType().fence();
           MPI_Irecv(buf,
                     size_forward_recv[iswap][i],MPI_DOUBLE,recvproc[iswap][i],0,world,&requests[i]);
         }
@@ -279,6 +282,7 @@ void CommTiledKokkos::reverse_comm_device()
         for (i = 0; i < nsend; i++) {
           buf = k_buf_recv.view<DeviceType>().data() +
             reverse_recv_offset[iswap][i]*k_buf_recv.view<DeviceType>().extent(1);
+          DeviceType().fence();
           MPI_Irecv(buf,
                     size_reverse_recv[iswap][i],MPI_DOUBLE,sendproc[iswap][i],0,world,&requests[i]);
         }
@@ -313,6 +317,7 @@ void CommTiledKokkos::reverse_comm_device()
         for (i = 0; i < nsend; i++) {
           buf = k_buf_recv.view<DeviceType>().data() +
             reverse_recv_offset[iswap][i]*k_buf_recv.view<DeviceType>().extent(1);
+          DeviceType().fence();
           MPI_Irecv(buf,
                     size_reverse_recv[iswap][i],MPI_DOUBLE,sendproc[iswap][i],0,world,&requests[i]);
         }
diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp
index 181100b11bf..79fbeb7cc2f 100644
--- a/src/REPLICA/fix_neb.cpp
+++ b/src/REPLICA/fix_neb.cpp
@@ -764,7 +764,7 @@ void FixNEB::inter_replica_comm()
   if (ireplica > 0) {
     if (me == 0) MPI_Waitall(2, requests, statuses);
 
-    MPI_Bcast(tagrecvall, nebatoms, MPI_INT, 0, world);
+    MPI_Bcast(tagrecvall, nebatoms, MPI_LMP_TAGINT, 0, world);
     MPI_Bcast(xrecvall[0], 3 * nebatoms, MPI_DOUBLE, 0, world);
 
     for (i = 0; i < nebatoms; i++) {
@@ -790,7 +790,7 @@ void FixNEB::inter_replica_comm()
   if (ireplica < nreplica - 1) {
     if (me == 0) MPI_Waitall(2, requests, statuses);
 
-    MPI_Bcast(tagrecvall, nebatoms, MPI_INT, 0, world);
+    MPI_Bcast(tagrecvall, nebatoms, MPI_LMP_TAGINT, 0, world);
     MPI_Bcast(xrecvall[0], 3 * nebatoms, MPI_DOUBLE, 0, world);
     MPI_Bcast(frecvall[0], 3 * nebatoms, MPI_DOUBLE, 0, world);
 
diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp
index 090de40ee31..cf917279b21 100644
--- a/src/RIGID/fix_shake.cpp
+++ b/src/RIGID/fix_shake.cpp
@@ -109,10 +109,10 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
   bool allow_typelabels = (atom->labelmapflag != 0);
   if (allow_typelabels) {
     for (int i = Atom::ATOM; i < Atom::DIHEDRAL; ++i) {
-      if ((atom->lmap->find("b", i) >= 0) ||
-          (atom->lmap->find("a", i) >= 0) ||
-          (atom->lmap->find("t", i) >= 0) ||
-          (atom->lmap->find("m", i) >= 0)) allow_typelabels = false;
+      if ((atom->lmap->find_type("b", i) >= 0) ||
+          (atom->lmap->find_type("a", i) >= 0) ||
+          (atom->lmap->find_type("t", i) >= 0) ||
+          (atom->lmap->find_type("m", i) >= 0)) allow_typelabels = false;
     }
     if (!allow_typelabels && (comm->me == 0))
       error->warning(FLERR, "At least one typelabel conflicts with a fix shake option: "
diff --git a/src/SPIN/fix_neb_spin.cpp b/src/SPIN/fix_neb_spin.cpp
index 2517938b9f6..74afd9cb875 100644
--- a/src/SPIN/fix_neb_spin.cpp
+++ b/src/SPIN/fix_neb_spin.cpp
@@ -890,7 +890,7 @@ void FixNEBSpin::inter_replica_comm()
   if (ireplica > 0) {
     if (me == 0) MPI_Waitall(2,requests,statuses);
 
-    MPI_Bcast(tagrecvall,nebatoms,MPI_INT,0,world);
+    MPI_Bcast(tagrecvall,nebatoms,MPI_LMP_TAGINT,0,world);
     MPI_Bcast(xrecvall[0],3*nebatoms,MPI_DOUBLE,0,world);
     MPI_Bcast(sprecvall[0],3*nebatoms,MPI_DOUBLE,0,world);
 
@@ -925,7 +925,7 @@ void FixNEBSpin::inter_replica_comm()
   if (ireplica < nreplica-1) {
     if (me == 0) MPI_Waitall(2,requests,statuses);
 
-    MPI_Bcast(tagrecvall,nebatoms,MPI_INT,0,world);
+    MPI_Bcast(tagrecvall,nebatoms,MPI_LMP_TAGINT,0,world);
     MPI_Bcast(xrecvall[0],3*nebatoms,MPI_DOUBLE,0,world);
     MPI_Bcast(frecvall[0],3*nebatoms,MPI_DOUBLE,0,world);
     MPI_Bcast(sprecvall[0],3*nebatoms,MPI_DOUBLE,0,world);
diff --git a/src/atom.cpp b/src/atom.cpp
index d1de419315b..5d7c36743ea 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -1267,7 +1267,7 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
           case 1: {    // type label
             if (!labelmapflag)
               error->one(FLERR, "Invalid line in {}: {}", location, utils::trim(buf));
-            type[nlocal - 1] = lmap->find(typestr, Atom::ATOM);
+            type[nlocal - 1] = lmap->find_type(typestr, Atom::ATOM);
             if (type[nlocal - 1] == -1)
               error->one(FLERR, "Invalid line in {}: {}", location, utils::trim(buf));
             break;
@@ -1374,7 +1374,7 @@ void Atom::data_bonds(int n, char *buf, int *count, tagint id_offset,
         }
         case 1: {    // type label
           if (!atom->labelmapflag) error->all(FLERR, "Invalid {}: {}", location, utils::trim(buf));
-          itype = lmap->find(typestr, Atom::BOND);
+          itype = lmap->find_type(typestr, Atom::BOND);
           if (itype == -1) error->all(FLERR, "Invalid {}: {}", location, utils::trim(buf));
           break;
         }
@@ -1469,7 +1469,7 @@ void Atom::data_angles(int n, char *buf, int *count, tagint id_offset,
         }
         case 1: {    // type label
           if (!atom->labelmapflag) error->all(FLERR, "Invalid {}: {}", location, utils::trim(buf));
-          itype = lmap->find(typestr, Atom::ANGLE);
+          itype = lmap->find_type(typestr, Atom::ANGLE);
           if (itype == -1) error->all(FLERR, "Invalid {}: {}", location, utils::trim(buf));
           break;
         }
@@ -1581,7 +1581,7 @@ void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset,
         }
         case 1: {    // type label
           if (!atom->labelmapflag) error->all(FLERR, "Invalid {}: {}", location, utils::trim(buf));
-          itype = lmap->find(typestr, Atom::DIHEDRAL);
+          itype = lmap->find_type(typestr, Atom::DIHEDRAL);
           if (itype == -1) error->all(FLERR, "Invalid {}: {}", location, utils::trim(buf));
           break;
         }
@@ -1709,7 +1709,7 @@ void Atom::data_impropers(int n, char *buf, int *count, tagint id_offset,
         }
         case 1: {    // type label
           if (!atom->labelmapflag) error->all(FLERR, "Invalid {}: {}", location, utils::trim(buf));
-          itype = lmap->find(typestr, Atom::IMPROPER);
+          itype = lmap->find_type(typestr, Atom::IMPROPER);
           if (itype == -1) error->all(FLERR, "Invalid {}: {}", location, utils::trim(buf));
           break;
         }
@@ -1968,7 +1968,7 @@ void Atom::set_mass(const char *file, int line, const char *str, int type_offset
     case 1: {    // type label
       if (!atom->labelmapflag)
         error->all(file, line, "Invalid atom type in {}: {}", location, utils::trim(str));
-      itype = lmap->find(typestr, Atom::ATOM);
+      itype = lmap->find_type(typestr, Atom::ATOM);
       if (itype == -1)
         error->all(file, line, "Unknown atom type {} in {}: {}", typestr, location,
                    utils::trim(str));
diff --git a/src/graphics.h b/src/graphics.h
index e5a14e011a4..81e9ebf932b 100644
--- a/src/graphics.h
+++ b/src/graphics.h
@@ -16,28 +16,26 @@
 
 // common definitions and declarations for graphics support in LAMMPS
 
-namespace LAMMPS_NS {
-namespace Graphics {
-  enum {
-    NONE,
-    SPHERE,    // a single sphere with radius provided
-    LINE,      // a cylinder with diameter given through fflag2
-    TRI,       // a surface mesh as triangles or cylinder mesh based on fflag1, fflag2 sets diameter
-    CYLINDER,    // a cylinder with diameter given by fix, fflag1 choose caps, fflag2 adjusts diameter
-    TRIANGLE,    // a regular triangle, no settings apply
-    BOND,     // two connected cylinders with bond diameter, colored by atom types, fflag1 sets cap
-    ARROW,    // a cylinder with a conical tip and a flat cap at the bottom
-    CONE,     // a truncated cone with flat caps, fflag1 sets caps
-    PIXMAP    // a pointer to a pixmap buffer at x,y,z location
-  };    // used by some Body and Fix child classes
+namespace LAMMPS_NS::Graphics {
+enum {
+  NONE,
+  SPHERE,      // a single sphere with radius provided
+  LINE,        // a cylinder with diameter given through fflag2
+  TRI,         // a surface mesh as triangles or cylinder mesh based on fflag1, fflag2 sets diameter
+  CYLINDER,    // a cylinder with diameter given by fix, fflag1 choose caps, fflag2 adjusts diameter
+  TRIANGLE,    // a regular triangle, no settings apply
+  BOND,        // two connected cylinders with bond diameter, colored by atom types, fflag1 sets cap
+  ARROW,       // a cylinder with a conical tip and a flat cap at the bottom
+  CONE,        // a truncated cone with flat caps, fflag1 sets caps
+  PIXMAP       // a pointer to a pixmap buffer at x,y,z location
+};    // used by some Body and Fix child classes
 
-  // definitions for rendering caps and sides of a truncated cone
-  enum {
-    CONE_TOP = 1 << 0,                            // draw top cap of cone/cylinder
-    CONE_BOT = 1 << 1,                            // draw bottom cap of cone/cylinder
-    CONE_SIDE = 1 << 2,                           // draw side of cone/cylinder
-    CONE_ALL = CONE_TOP | CONE_BOT | CONE_SIDE    // all of the above
-  };
-}    // namespace Graphics
-}    // namespace LAMMPS_NS
+// definitions for rendering caps and sides of a truncated cone
+enum {
+  CONE_TOP = 1 << 0,                            // draw top cap of cone/cylinder
+  CONE_BOT = 1 << 1,                            // draw bottom cap of cone/cylinder
+  CONE_SIDE = 1 << 2,                           // draw side of cone/cylinder
+  CONE_ALL = CONE_TOP | CONE_BOT | CONE_SIDE    // all of the above
+};
+}    // namespace LAMMPS_NS::Graphics
 #endif
diff --git a/src/label_map.cpp b/src/label_map.cpp
index a4bc67d9a04..1e176037ed9 100644
--- a/src/label_map.cpp
+++ b/src/label_map.cpp
@@ -19,6 +19,7 @@
 #include "error.h"
 #include "force.h"
 #include "improper.h"
+#include "tokenizer.h"
 
 #include <algorithm>
 #include <cstring>
@@ -271,7 +272,7 @@ int LabelMap::find_or_create(const std::string &mylabel, std::vector<std::string
    return -1 if type not yet defined
 ------------------------------------------------------------------------- */
 
-int LabelMap::find(const std::string &mylabel, int mode) const
+int LabelMap::find_type(const std::string &mylabel, int mode) const
 {
   switch (mode) {
     case Atom::ATOM:
@@ -299,41 +300,36 @@ int LabelMap::find(const std::string &mylabel, int mode) const
    return "" if type label does not exist
 ------------------------------------------------------------------------- */
 
-const std::string &LabelMap::find(int i, int mode) const
+const std::string &LabelMap::find_label(int i, int mode) const
 {
   switch (mode) {
-  case Atom::ATOM:
-    if ((i > 0) && (i <= atom->ntypes)) {
-      if (is_complete(mode))
-        return typelabel[i-1];
-    }
-    break;
-  case Atom::BOND:
-    if ((i > 0) && (i <= atom->nbondtypes)) {
-      if (is_complete(mode))
-        return btypelabel[i-1];
-    }
-    break;
-  case Atom::ANGLE:
-    if ((i > 0) && (i <= atom->nangletypes)) {
-      if (is_complete(mode))
-        return atypelabel[i-1];
-    }
-    break;
-  case Atom::DIHEDRAL:
-    if ((i > 0) && (i <= atom->ndihedraltypes)) {
-      if (is_complete(mode))
-        return dtypelabel[i-1];
-    }
-    break;
-  case Atom::IMPROPER:
-    if ((i > 0) && (i <= atom->nimpropertypes)) {
-      if (is_complete(mode))
-        return itypelabel[i-1];
-    }
-    break;
-  default:
-    return empty;
+    case Atom::ATOM:
+      if ((i > 0) && (i <= atom->ntypes)) {
+        if (is_complete(mode)) return typelabel[i - 1];
+      }
+      break;
+    case Atom::BOND:
+      if ((i > 0) && (i <= atom->nbondtypes)) {
+        if (is_complete(mode)) return btypelabel[i - 1];
+      }
+      break;
+    case Atom::ANGLE:
+      if ((i > 0) && (i <= atom->nangletypes)) {
+        if (is_complete(mode)) return atypelabel[i - 1];
+      }
+      break;
+    case Atom::DIHEDRAL:
+      if ((i > 0) && (i <= atom->ndihedraltypes)) {
+        if (is_complete(mode)) return dtypelabel[i - 1];
+      }
+      break;
+    case Atom::IMPROPER:
+      if ((i > 0) && (i <= atom->nimpropertypes)) {
+        if (is_complete(mode)) return itypelabel[i - 1];
+      }
+      break;
+    default:
+      return empty;
   }
   return empty;
 }
@@ -385,9 +381,13 @@ bool LabelMap::is_complete(int mode) const
 
 int LabelMap::infer_bondtype(int type1, int type2)
 {
+  // check for out of range input
+  if ((type1 < 1) || (type1 > natomtypes) || (type2 < 1) || (type2 > natomtypes)) return -1;
+
+  // convert numeric atom types to type label
   std::vector<std::string> mytypes(2);
-  mytypes[0] = typelabel[type1-1];
-  mytypes[1] = typelabel[type2-1];
+  mytypes[0] = typelabel[type1 - 1];
+  mytypes[1] = typelabel[type2 - 1];
   if (mytypes[0].empty() || mytypes[1].empty()) return -1;
 
   return infer_bondtype(mytypes);
@@ -398,21 +398,20 @@ int LabelMap::infer_bondtype(int type1, int type2)
    assumes bond types are of the form "a-b" for atom types 'a' and 'b'
 ------------------------------------------------------------------------- */
 
-int LabelMap::infer_bondtype(std::vector<std::string> mytypes)
+int LabelMap::infer_bondtype(const std::vector<std::string> &mytypes)
 {
   // search for matching bond type label with symmetry considerations
-
   std::vector<std::string> btypes(2);
   for (int i = 0; i < nbondtypes; i++) {
     int status = parse_typelabel(2, btypelabel[i], btypes);
-    if (status != -1)
+    if ((status != -1) && (btypes.size() == 2))
       if ((mytypes[0] == btypes[0] && mytypes[1] == btypes[1]) ||
-          (mytypes[0] == btypes[1] && mytypes[1] == btypes[0])) return i+1;
+          (mytypes[0] == btypes[1] && mytypes[1] == btypes[0]))
+        return i + 1;
   }
   return -1;
 }
 
-
 /* ----------------------------------------------------------------------
    infer angle type from three atom types
    input/output is numeric types, uses type labels internally
@@ -421,12 +420,16 @@ int LabelMap::infer_bondtype(std::vector<std::string> mytypes)
 
 int LabelMap::infer_angletype(int type1, int type2, int type3)
 {
-  // convert numeric atom types to type label
+  // check for out of range input
+  if ((type1 < 1) || (type1 > natomtypes) || (type2 < 1) || (type2 > natomtypes) || (type3 < 1) ||
+      (type3 > natomtypes))
+    return -1;
 
+  // convert numeric atom types to type label
   std::vector<std::string> mytypes(3);
-  mytypes[0] = typelabel[type1-1];
-  mytypes[1] = typelabel[type2-1];
-  mytypes[2] = typelabel[type3-1];
+  mytypes[0] = typelabel[type1 - 1];
+  mytypes[1] = typelabel[type2 - 1];
+  mytypes[2] = typelabel[type3 - 1];
   for (size_t i = 0; i < 3; i++)
     if (mytypes[i].empty()) return -1;
 
@@ -439,7 +442,7 @@ int LabelMap::infer_angletype(int type1, int type2, int type3)
    assumes angle types of the form "a-b-c" for atom types 'a', 'b', 'c'
 ------------------------------------------------------------------------- */
 
-int LabelMap::infer_angletype(std::vector<std::string> mytypes)
+int LabelMap::infer_angletype(const std::vector<std::string> &mytypes)
 {
   // search for matching angle type label, with symmetry considerations
 
@@ -449,12 +452,12 @@ int LabelMap::infer_angletype(std::vector<std::string> mytypes)
     status = parse_typelabel(3, atypelabel[i], atypes);
     if (status != -1 && mytypes[1] == atypes[1])
       if ((mytypes[0] == atypes[0] && mytypes[2] == atypes[2]) ||
-          (mytypes[0] == atypes[2] && mytypes[2] == atypes[0])) return i+1;
+          (mytypes[0] == atypes[2] && mytypes[2] == atypes[0]))
+        return i + 1;
   }
   return -1;
 }
 
-
 /* ----------------------------------------------------------------------
    infer dihedral type from four atom types
    input/output is numeric types, uses type labels internally
@@ -463,13 +466,17 @@ int LabelMap::infer_angletype(std::vector<std::string> mytypes)
 
 int LabelMap::infer_dihedraltype(int type1, int type2, int type3, int type4)
 {
-  // convert numeric atom types to type label
+  // check for out of range input
+  if ((type1 < 1) || (type1 > natomtypes) || (type2 < 1) || (type2 > natomtypes) || (type3 < 1) ||
+      (type3 > natomtypes) || (type4 < 1) || (type4 > natomtypes))
+    return -1;
 
+  // convert numeric atom types to type label
   std::vector<std::string> mytypes(4);
-  mytypes[0] = typelabel[type1-1];
-  mytypes[1] = typelabel[type2-1];
-  mytypes[2] = typelabel[type3-1];
-  mytypes[3] = typelabel[type4-1];
+  mytypes[0] = typelabel[type1 - 1];
+  mytypes[1] = typelabel[type2 - 1];
+  mytypes[2] = typelabel[type3 - 1];
+  mytypes[3] = typelabel[type4 - 1];
   for (size_t i = 0; i < 4; i++)
     if (mytypes[i].empty()) return -1;
 
@@ -482,7 +489,7 @@ int LabelMap::infer_dihedraltype(int type1, int type2, int type3, int type4)
    assumes dihedral types of the form "a-b-c-d"
 ------------------------------------------------------------------------- */
 
-int LabelMap::infer_dihedraltype(std::vector<std::string> mytypes)
+int LabelMap::infer_dihedraltype(const std::vector<std::string> &mytypes)
 {
   // search for matching dihedral type label
 
@@ -491,10 +498,11 @@ int LabelMap::infer_dihedraltype(std::vector<std::string> mytypes)
   for (int i = 0; i < ndihedraltypes; i++) {
     status = parse_typelabel(4, dtypelabel[i], dtypes);
     if (status != -1)
-      if ((mytypes[0] == dtypes[0] && mytypes[1] == dtypes[1] &&
-          mytypes[2] == dtypes[2] && mytypes[3] == dtypes[3]) ||
-          (mytypes[3] == dtypes[0] && mytypes[2] == dtypes[1] &&
-           mytypes[1] == dtypes[2] && mytypes[0] == dtypes[3])) return i+1;
+      if ((mytypes[0] == dtypes[0] && mytypes[1] == dtypes[1] && mytypes[2] == dtypes[2] &&
+           mytypes[3] == dtypes[3]) ||
+          (mytypes[3] == dtypes[0] && mytypes[2] == dtypes[1] && mytypes[1] == dtypes[2] &&
+           mytypes[0] == dtypes[3]))
+        return i + 1;
   }
   return -1;
 }
@@ -508,13 +516,17 @@ int LabelMap::infer_dihedraltype(std::vector<std::string> mytypes)
 
 int LabelMap::infer_impropertype(int type1, int type2, int type3, int type4)
 {
-  // convert numeric atom types to type label
+  // check for out of range input
+  if ((type1 < 1) || (type1 > natomtypes) || (type2 < 1) || (type2 > natomtypes) || (type3 < 1) ||
+      (type3 > natomtypes) || (type4 < 1) || (type4 > natomtypes))
+    return -1;
 
+  // convert numeric atom types to type label
   std::vector<std::string> mytypes(4);
-  mytypes[0] = typelabel[type1-1];
-  mytypes[1] = typelabel[type2-1];
-  mytypes[2] = typelabel[type3-1];
-  mytypes[3] = typelabel[type4-1];
+  mytypes[0] = typelabel[type1 - 1];
+  mytypes[1] = typelabel[type2 - 1];
+  mytypes[2] = typelabel[type3 - 1];
+  mytypes[3] = typelabel[type4 - 1];
   for (int i = 0; i < 4; i++)
     if (mytypes[i].empty()) return -1;
 
@@ -528,7 +540,7 @@ int LabelMap::infer_impropertype(int type1, int type2, int type3, int type4)
    the symmetry of the improper is encoded in improper.symmatoms
 ------------------------------------------------------------------------- */
 
-int LabelMap::infer_impropertype(std::vector<std::string> mytypes)
+int LabelMap::infer_impropertype(const std::vector<std::string> &mytypes)
 {
   // search for matching improper type label
 
@@ -541,7 +553,7 @@ int LabelMap::infer_impropertype(std::vector<std::string> mytypes)
     status = parse_typelabel(4, itypelabel[i], itypes);
     if (status != -1) {
       for (int j = 0; j < 4; j++) {
-        if (force->improper->symmatoms[j] == 1) {
+        if (force->improper && (force->improper->symmatoms[j] == 1)) {
           if (mytypes[j] != itypes[j]) {
             status = -1;
             break;
@@ -552,14 +564,14 @@ int LabelMap::infer_impropertype(std::vector<std::string> mytypes)
         }
       }
       if (status == -1) continue;
-      std::sort(list1.begin(),list1.end());
-      std::sort(list2.begin(),list2.end());
+      std::sort(list1.begin(), list1.end());
+      std::sort(list2.begin(), list2.end());
       for (int j = 0; j < nlist; j++)
         if (list1[j] != list2[j]) {
           status = -1;
           break;
         }
-      if (status != -1) return i+1;
+      if (status != -1) return i + 1;
     }
   }
   return -1;
@@ -569,18 +581,11 @@ int LabelMap::infer_impropertype(std::vector<std::string> mytypes)
    return -1 if number of parsed strings is not equal to ntypes input
 ------------------------------------------------------------------------- */
 
-int LabelMap::parse_typelabel(int ntypes, std::string label, std::vector<std::string> &types)
+int LabelMap::parse_typelabel(int ntypes, const std::string &label, std::vector<std::string> &types)
 {
-  std::vector<std::string> out;
-  size_t start = label.find_first_not_of('-');
-
-  while (start != std::string::npos) {
-      size_t end = label.find('-', start);
-      out.emplace_back(label.substr(start, end - start));
-      start = label.find_first_not_of('-', end);
-  }
-  if (out.size() != ntypes) return -1;
-  types = out;
+  auto out = Tokenizer(label,"-").as_vector();
+  if ((int)out.size() != ntypes) return -1;
+  types = std::move(out);
   return 1;
 }
 
diff --git a/src/label_map.h b/src/label_map.h
index 4b3632fb2d0..c736fdcdae3 100644
--- a/src/label_map.h
+++ b/src/label_map.h
@@ -11,6 +11,8 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
+/*! \file label_map.h */
+
 #ifndef LMP_LABEL_MAP_H
 #define LMP_LABEL_MAP_H
 
@@ -20,6 +22,25 @@
 
 namespace LAMMPS_NS {
 
+/*! \class LabelMap
+ *  \brief Manage type labels for atoms, bonds, angles, dihedrals, and impropers
+ *
+ * The LabelMap class provides functionality to map between string labels and
+ * numeric type indices for atoms, bonds, angles, dihedrals, and impropers in LAMMPS.
+ * This enables users to reference types by symbolic names (e.g., "C", "H", "C-H")
+ * instead of numeric indices, improving readability and maintainability of input scripts.
+ *
+ * Type labels for bonded interactions *may* (but are not required to) use a
+ * hyphen-delimited format indicating the types of the constituent atoms.  Examples:
+ * - Bond types: "atom1-atom2" (e.g., "C-H", "N-O")
+ * - Angle types: "atom1-atom2-atom3" (e.g., "H-C-H", "C-N-C")
+ * - Dihedral types: "atom1-atom2-atom3-atom4" (e.g., "C-C-N-H")
+ * - Improper types: "atom1-atom2-atom3-atom4" (e.g., "C-N-C-C")
+ *
+ * The class supports bidirectional lookup (label <-> type) and can infer bonded
+ * interaction types from constituent atom types when using a hyphen-delimited
+ * format convention. */
+
 class LabelMap : protected Pointers {
   friend class AtomVec;
   friend class DumpCustom;
@@ -28,71 +49,255 @@ class LabelMap : protected Pointers {
   friend class ReadData;
 
  public:
+  /*! Construct a LabelMap instance
+   *
+   * \param  lmp              Pointer to LAMMPS instance
+   * \param  natomtypes       Number of atom types in map
+   * \param  nbondtypes       Number of bond types in map
+   * \param  nangletypes      Number of angle types in map
+   * \param  ndihedraltypes   Number of dihedral types in map
+   * \param  nimpropertypes   Number of improper types in map */
   LabelMap(LAMMPS *lmp, int, int, int, int, int);
   ~LabelMap() override;
 
-  void modify_lmap(int, char **);              // labelmap command in the input script
-  void merge_lmap(LabelMap *, int);            // copy another lmap into this one
-  void create_lmap2lmap(LabelMap *, int);      // index mapping between two lmaps
-  int find(const std::string &, int) const;    // find numeric type of type label
-  const std::string &find(int, int) const;     // find type label for numeric type
-  bool is_complete(int) const;                 // check if all types are assigned
+  /*! Process labelmap command from input script
+   *
+   \verbatim embed:rst
+
+Add or modify type label mappings from the LAMMPS
+:doc:`labelmap <labelmap>` input command.
+
+   \endverbatim
+   *
+   * \param  narg  Number of arguments
+   * \param  arg   Array of argument strings */
+  void modify_lmap(int, char **);
+
+  /*! Copy another LabelMap into this one
+   *
+   \verbatim embed:rst
+
+Merge type labels from another LabelMap instance into the current one.
+Currently used when combining data from multiple sources with
+:doc:`read_data add <read_data>` or when replicating the system with
+:doc:`replicate <replicate>`.
+
+   \endverbatim
+   *
+   * \param  lmap  Pointer to source LabelMap
+   * \param  mode  Merge mode flag */
+  void merge_lmap(LabelMap *, int);
+
+  /*! Create index mapping between two LabelMaps
+   *
+   * Build a mapping structure (lmap2lmap) that translates type indices
+   * from another LabelMap to the current one based on matching labels.
+   *
+   * \param  lmap  Pointer to source LabelMap
+   * \param  mode  Mapping mode flag */
+  void create_lmap2lmap(LabelMap *, int);
+
+  /*! Find numeric type from type label
+   *
+   * Look up the numeric type index corresponding to a type label string.
+   *
+   * \param  mylabel  Type label string to search for
+   * \param  mode     Type category: Atom::ATOM, Atom::BOND, Atom::ANGLE,
+   *                  Atom::DIHEDRAL, or Atom::IMPROPER
+   * \return          Numeric type index (1-based), or -1 if not found */
+  int find_type(const std::string &, int) const;
+
+  /*! Find type label from numeric type
+   *
+   * Reverse lookup: retrieve the type label string for a given numeric type.
+   *
+   * \param  i     Numeric type index (1-based)
+   * \param  mode  Type category: Atom::ATOM, Atom::BOND, Atom::ANGLE,
+   *               Atom::DIHEDRAL, or Atom::IMPROPER
+   * \return       Reference to type label string, or empty string if not found */
+  const std::string &find_label(int, int) const;
 
-  // infer interaction types from standard hyphen-delimited format
+  /*! Check if all types have assigned labels
+   *
+   * Verify that every type in the specified category has a corresponding label.
+   *
+   * \param  mode  Type category: Atom::ATOM, Atom::BOND, Atom::ANGLE,
+   *               Atom::DIHEDRAL, or Atom::IMPROPER
+   * \return       True if all types have labels, false otherwise */
+  bool is_complete(int) const;
 
-  int infer_bondtype(int, int);                      // infer bond type from two atom types
-  int infer_bondtype(std::vector<std::string>);      // infer bond type from two atom type labels
+  /*! \name Interaction type inference from hyphen-delimited labels
+   *
+   * These methods infer bonded interaction types (bonds, angles, dihedrals, impropers)
+   * from constituent atom types using a hyphen-delimited format.  This requires that
+   * type labels for bonded interactions were entered following this convention.
+   * The inference functions consider the symmetry of the interaction and thus atom
+   * types may be swapped accordingly and the bonded type will still be matched.
+   * @{ */
 
-  int infer_angletype(int, int, int);                // infer angle type from three atom types
-  int infer_angletype(std::vector<std::string>);     // infer angle type from three atom type labels
+  /*! Infer bond type from two numeric atom types
+   *
+   * Look up or create a bond type from two atom type indices by constructing
+   * a hyphen-delimited label (e.g., "C-H") and searching the bond type labels.
+   *
+   * \param  atype1  First atom type index
+   * \param  atype2  Second atom type index
+   * \return         Bond type index, or -1 if not found */
+  int infer_bondtype(int, int);
 
-  int infer_dihedraltype(int, int, int, int);        // infer dihedral type from four atom types
-  int infer_dihedraltype(std::vector<std::string>);  // infer dihedral type from four atom type labels
+  /*! Infer bond type from atom type labels
+   *
+   * \overload
+   *
+   * Look up a bond type from two atom type labels.
+   *
+   * \param  labels  Vector of two atom type label strings
+   * \return         Bond type index, or -1 if not found */
+  int infer_bondtype(const std::vector<std::string> &);
 
-  int infer_impropertype(int, int, int, int);        // infer improper type from four atom types
-  int infer_impropertype(std::vector<std::string>);  // infer improper type from four atom type labels
+  /*! Infer angle type from three numeric atom types
+   *
+   * Look up or create an angle type from three atom type indices by
+   * constructing a hyphen-delimited label (e.g., "H1-C1-H2").
+   *
+   * \param  atype1  First atom type index
+   * \param  atype2  Second atom type index (center atom)
+   * \param  atype3  Third atom type index
+   * \return         Angle type index, or -1 if not found */
+  int infer_angletype(int, int, int);
 
-  int parse_typelabel(int, std::string, std::vector<std::string> &); // get strings within hyphen delimiters
+  /*! Infer angle type from three atom type labels
+   *
+   * \overload
+   *
+   * Look up an angle type from three atom type labels.
+   *
+   * \param  labels  Vector of three atom type label strings
+   * \return         Angle type index, or -1 if not found */
+  int infer_angletype(const std::vector<std::string> &);
 
-  // input/output for atom class label map
+  /*! Infer dihedral type from four numeric atom types
+   *
+   * Look up a dihedral type from four atom type indices by
+   * constructing a hyphen-delimited label (e.g., "C-C-N-H").
+   *
+   * \param  atype1  First atom type index
+   * \param  atype2  Second atom type index
+   * \param  atype3  Third atom type index
+   * \param  atype4  Fourth atom type index
+   * \return         Dihedral type index, or -1 if not found */
+  int infer_dihedraltype(int, int, int, int);
 
+  /*! Infer dihedral type from atom type labels
+   *
+   * \overload
+   *
+   * Look up a dihedral type from four atom type labels.
+   *
+   * \param  labels  Vector of four atom type label strings
+   * \return         Dihedral type index, or -1 if not found */
+  int infer_dihedraltype(const std::vector<std::string> &);
+
+  /*! Infer improper type from four numeric atom types
+   *
+   * Look up an improper type from four atom type indices by
+   * constructing a hyphen-delimited label (e.g., "C-N-C-C").
+   *
+   * \param  atype1  First atom type index (center atom)
+   * \param  atype2  Second atom type index
+   * \param  atype3  Third atom type index
+   * \param  atype4  Fourth atom type index
+   * \return         Improper type index, or -1 if not found */
+  int infer_impropertype(int, int, int, int);
+
+  /*! Infer improper type from atom type labels
+   *
+   * \overload
+   *
+   * Look up an improper type from four atom type labels.
+   *
+   * \param  labels  Vector of four atom type label strings
+   * \return         Improper type index, or -1 if not found */
+  int infer_impropertype(const std::vector<std::string> &);
+
+  /*! @} */
+
+  /*! Parse hyphen-delimited type label into components
+   *
+   * Split a hyphen-delimited label (e.g., "C-N-H") into individual type strings.
+   * Validates that the number of components matches the expected count.
+   *
+   * \param  ntypes  Expected number of components
+   * \param  label   Hyphen-delimited label string
+   * \param  types   Output vector to store component strings
+   * \return         0 on success, -1 if component count doesn't match ntypes */
+  int parse_typelabel(int, const std::string &, std::vector<std::string> &);
+
+  /*! \name I/O methods for label map persistence
+   * @{ */
+
+  /*! Write label map to data file
+   *
+   * Output all type labels as sections to a LAMMPS data file.
+   *
+   * \param  fp  File pointer for writing */
   void write_data(FILE *);
+
+  /*! Read label map from restart file
+   *
+   * Restore label map data from a LAMMPS restart file.
+   *
+   * \param  fp  File pointer for reading */
   void read_restart(FILE *fp);
+
+  /*! Write label map to restart file
+   *
+   * Save label map data to a LAMMPS restart file for later restoration.
+   *
+   * \param  fp  File pointer for writing */
   void write_restart(FILE *);
 
-protected:
-  int natomtypes, nbondtypes, nangletypes, ndihedraltypes, nimpropertypes;
-  std::vector<std::string> typelabel, btypelabel, atypelabel;
-  std::vector<std::string> dtypelabel, itypelabel;
+  /*! @} */
+
+ protected:
+  int natomtypes, nbondtypes, nangletypes, ndihedraltypes, nimpropertypes;    //!< Type counts
+  std::vector<std::string> typelabel, btypelabel,
+      atypelabel;                                     //!< Label storage (atoms, bonds, angles)
+  std::vector<std::string> dtypelabel, itypelabel;    //!< Label storage (dihedrals, impropers)
 
-  std::unordered_map<std::string, int> typelabel_map;
-  std::unordered_map<std::string, int> btypelabel_map;
-  std::unordered_map<std::string, int> atypelabel_map;
-  std::unordered_map<std::string, int> dtypelabel_map;
-  std::unordered_map<std::string, int> itypelabel_map;
+  std::unordered_map<std::string, int> typelabel_map;     //!< Atom label → type mapping
+  std::unordered_map<std::string, int> btypelabel_map;    //!< Bond label → type mapping
+  std::unordered_map<std::string, int> atypelabel_map;    //!< Angle label → type mapping
+  std::unordered_map<std::string, int> dtypelabel_map;    //!< Dihedral label → type mapping
+  std::unordered_map<std::string, int> itypelabel_map;    //!< Improper label → type mapping
 
-  // per-type data struct mapping this label map to another
+  /*! \struct Lmap2Lmap
+   *  \brief Mapping structure between two LabelMaps
+   *
+   * Stores per-type index mappings from another LabelMap to this one,
+   * enabling type translation when merging or comparing different label maps. */
 
   struct Lmap2Lmap {
-    int *atom;
-    int *bond;
-    int *angle;
-    int *dihedral;
-    int *improper;
+    int *atom;        //!< Atom type mapping array
+    int *bond;        //!< Bond type mapping array
+    int *angle;       //!< Angle type mapping array
+    int *dihedral;    //!< Dihedral type mapping array
+    int *improper;    //!< Improper type mapping array
   };
 
-  Lmap2Lmap lmap2lmap;
+  Lmap2Lmap lmap2lmap;    //!< Instance of inter-map translation data
 
-  void reset_type_labels();
-  int find_or_create(const std::string &, std::vector<std::string> &,
-                     std::unordered_map<std::string, int> &);    // look up type or create new type
+  void reset_type_labels();    //!< Clear all type labels
+  int find_or_create(
+      const std::string &, std::vector<std::string> &,
+      std::unordered_map<std::string, int> &);    //!< Look up type or create new type
   int search(const std::string &,
-             const std::unordered_map<std::string, int> &) const;    // look up type index
-  char *read_string(FILE *);
-  void write_string(const std::string &, FILE *);
-  int read_int(FILE *);
+             const std::unordered_map<std::string, int> &) const;    //!< Look up type index
+  char *read_string(FILE *);                         //!< Read string from binary file
+  void write_string(const std::string &, FILE *);    //!< Write string to binary file
+  int read_int(FILE *);                              //!< Read integer from binary file
 
-  void write_map(const std::string &);
+  void write_map(const std::string &);    //!< Write label map to file for debugging
 };
 
 }    // namespace LAMMPS_NS
diff --git a/src/molecule.cpp b/src/molecule.cpp
index 5a06144795a..e3a6f41a3df 100644
--- a/src/molecule.cpp
+++ b/src/molecule.cpp
@@ -671,7 +671,7 @@ void Molecule::from_json(const std::string &molid, const json &moldata)
           error->all(FLERR, Error::NOLASTLINE,
                      "Molecule template {}: invalid atom type in \"types\" JSON section", id,
                      typestr);
-        type[iatom] = atom->lmap->find(typestr, Atom::ATOM);
+        type[iatom] = atom->lmap->find_type(typestr, Atom::ATOM);
         if (type[iatom] == -1)
           error->all(FLERR, Error::NOLASTLINE,
                      "Molecule template {}: Unknown atom type {} in \"types\" JSON section", id,
@@ -1056,7 +1056,7 @@ void Molecule::from_json(const std::string &molid, const json &moldata)
               error->all(FLERR, Error::NOLASTLINE,
                          "Molecule template {}: invalid bond type in \"bonds\" JSON section", id,
                          typestr);
-            itype = atom->lmap->find(typestr, Atom::BOND);
+            itype = atom->lmap->find_type(typestr, Atom::BOND);
             if (itype == -1)
               error->all(FLERR, Error::NOLASTLINE,
                          "Molecule template {}: Unknown bond type {} in \"bonds\" JSON section", id,
@@ -1148,7 +1148,7 @@ void Molecule::from_json(const std::string &molid, const json &moldata)
               error->all(FLERR, Error::NOLASTLINE,
                          "Molecule template {}: invalid angle type in \"angles\" JSON section", id,
                          typestr);
-            itype = atom->lmap->find(typestr, Atom::ANGLE);
+            itype = atom->lmap->find_type(typestr, Atom::ANGLE);
             if (itype == -1)
               error->all(FLERR, Error::NOLASTLINE,
                          "Molecule template {}: Unknown angle type {} in \"angles\" JSON section",
@@ -1259,7 +1259,7 @@ void Molecule::from_json(const std::string &molid, const json &moldata)
                   FLERR, Error::NOLASTLINE,
                   "Molecule template {}: invalid dihedral type in \"dihedrals\" JSON section", id,
                   typestr);
-            itype = atom->lmap->find(typestr, Atom::DIHEDRAL);
+            itype = atom->lmap->find_type(typestr, Atom::DIHEDRAL);
             if (itype == -1)
               error->all(
                   FLERR, Error::NOLASTLINE,
@@ -1384,7 +1384,7 @@ void Molecule::from_json(const std::string &molid, const json &moldata)
                   FLERR, Error::NOLASTLINE,
                   "Molecule template {}: invalid improper type in \"impropers\" JSON section", id,
                   typestr);
-            itype = atom->lmap->find(typestr, Atom::IMPROPER);
+            itype = atom->lmap->find_type(typestr, Atom::IMPROPER);
             if (itype == -1)
               error->all(
                   FLERR, Error::NOLASTLINE,
@@ -1882,7 +1882,7 @@ json Molecule::to_json() const
     moldata["types"]["format"] = {"atom-id", "type"};
     if (atom->labelmapflag && atom->lmap->is_complete(Atom::ATOM)) {
       for (int i = 0; i < natoms; ++i) {
-        moldata["types"]["data"][i] = {i + 1, atom->lmap->find(type[i], Atom::ATOM)};
+        moldata["types"]["data"][i] = {i + 1, atom->lmap->find_label(type[i], Atom::ATOM)};
       }
     } else {
       for (int i = 0; i < natoms; ++i) moldata["types"]["data"][i] = {i + 1, type[i]};
@@ -1939,7 +1939,7 @@ json Molecule::to_json() const
       for (int j = 0; j < num_bond[i]; j++) {
         if (has_newton_bond || (i + 1 < bond_atom[i][j])) {
           if (has_typelabels) {
-            moldata["bonds"]["data"][idx] = {atom->lmap->find(bond_type[i][j], Atom::BOND), i + 1,
+            moldata["bonds"]["data"][idx] = {atom->lmap->find_label(bond_type[i][j], Atom::BOND), i + 1,
                                              bond_atom[i][j]};
           } else {
             moldata["bonds"]["data"][idx] = {bond_type[i][j], i + 1, bond_atom[i][j]};
@@ -1958,7 +1958,7 @@ json Molecule::to_json() const
       for (int j = 0; j < num_angle[i]; j++) {
         if (has_newton_bond || (i + 1 == angle_atom2[i][j])) {
           if (has_typelabels) {
-            moldata["angles"]["data"][idx] = {atom->lmap->find(angle_type[i][j], Atom::ANGLE),
+            moldata["angles"]["data"][idx] = {atom->lmap->find_label(angle_type[i][j], Atom::ANGLE),
                                               angle_atom1[i][j], angle_atom2[i][j],
                                               angle_atom3[i][j]};
           } else {
@@ -1980,7 +1980,7 @@ json Molecule::to_json() const
         if (has_newton_bond || (i + 1 == dihedral_atom2[i][j])) {
           if (has_typelabels) {
             moldata["dihedrals"]["data"][idx] = {
-                atom->lmap->find(dihedral_type[i][j], Atom::DIHEDRAL), dihedral_atom1[i][j],
+                atom->lmap->find_label(dihedral_type[i][j], Atom::DIHEDRAL), dihedral_atom1[i][j],
                 dihedral_atom2[i][j], dihedral_atom3[i][j], dihedral_atom4[i][j]};
           } else {
             moldata["dihedrals"]["data"][idx] = {dihedral_type[i][j], dihedral_atom1[i][j],
@@ -2002,7 +2002,7 @@ json Molecule::to_json() const
         if (has_newton_bond || (i + 1 == improper_atom2[i][j])) {
           if (has_typelabels) {
             moldata["impropers"]["data"][idx] = {
-                atom->lmap->find(improper_type[i][j], Atom::IMPROPER), improper_atom1[i][j],
+                atom->lmap->find_label(improper_type[i][j], Atom::IMPROPER), improper_atom1[i][j],
                 improper_atom2[i][j], improper_atom3[i][j], improper_atom4[i][j]};
           } else {
             moldata["impropers"]["data"][idx] = {improper_type[i][j], improper_atom1[i][j],
@@ -2043,9 +2043,9 @@ json Molecule::to_json() const
                                                      shake_atom[i][2]};
           if (has_typelabels) {
             moldata["shake"]["types"]["data"][i][1] = {
-                atom->lmap->find(shake_type[i][0], Atom::BOND),
-                atom->lmap->find(shake_type[i][1], Atom::BOND),
-                atom->lmap->find(shake_type[i][2], Atom::ANGLE)};
+                atom->lmap->find_label(shake_type[i][0], Atom::BOND),
+                atom->lmap->find_label(shake_type[i][1], Atom::BOND),
+                atom->lmap->find_label(shake_type[i][2], Atom::ANGLE)};
           } else {
             moldata["shake"]["types"]["data"][i][1] = {shake_type[i][0], shake_type[i][1],
                                                        shake_type[i][2]};
@@ -2055,8 +2055,8 @@ json Molecule::to_json() const
           moldata["shake"]["atoms"]["data"][i][1] = {shake_atom[i][0], shake_atom[i][1]};
           if (has_typelabels) {
             moldata["shake"]["types"]["data"][i][1] = {
-                atom->lmap->find(shake_type[i][0], Atom::BOND),
-                atom->lmap->find(shake_type[i][1], Atom::BOND)};
+                atom->lmap->find_label(shake_type[i][0], Atom::BOND),
+                atom->lmap->find_label(shake_type[i][1], Atom::BOND)};
           } else {
             moldata["shake"]["types"]["data"][i][1] = {shake_type[i][0], shake_type[i][1]};
           }
@@ -2066,9 +2066,9 @@ json Molecule::to_json() const
                                                      shake_atom[i][2]};
           if (has_typelabels) {
             moldata["shake"]["types"]["data"][i][1] = {
-                atom->lmap->find(shake_type[i][0], Atom::BOND),
-                atom->lmap->find(shake_type[i][1], Atom::BOND),
-                atom->lmap->find(shake_type[i][2], Atom::BOND)};
+                atom->lmap->find_label(shake_type[i][0], Atom::BOND),
+                atom->lmap->find_label(shake_type[i][1], Atom::BOND),
+                atom->lmap->find_label(shake_type[i][2], Atom::BOND)};
           } else {
             moldata["shake"]["types"]["data"][i][1] = {shake_type[i][0], shake_type[i][1],
                                                        shake_type[i][2]};
@@ -2079,10 +2079,10 @@ json Molecule::to_json() const
                                                      shake_atom[i][2], shake_atom[i][3]};
           if (has_typelabels) {
             moldata["shake"]["types"]["data"][i][1] = {
-                atom->lmap->find(shake_type[i][0], Atom::BOND),
-                atom->lmap->find(shake_type[i][1], Atom::BOND),
-                atom->lmap->find(shake_type[i][2], Atom::BOND),
-                atom->lmap->find(shake_type[i][3], Atom::BOND)};
+                atom->lmap->find_label(shake_type[i][0], Atom::BOND),
+                atom->lmap->find_label(shake_type[i][1], Atom::BOND),
+                atom->lmap->find_label(shake_type[i][2], Atom::BOND),
+                atom->lmap->find_label(shake_type[i][3], Atom::BOND)};
           } else {
             moldata["shake"]["types"]["data"][i][1] = {shake_type[i][0], shake_type[i][1],
                                                        shake_type[i][2], shake_type[i][3]};
@@ -2772,7 +2772,7 @@ void Molecule::types(char *line)
         if (!atom->labelmapflag)
           error->all(FLERR, fileiarg, "Invalid atom type {} in {}: {}", typestr, location,
                      utils::trim(line));
-        type[iatom] = atom->lmap->find(typestr, Atom::ATOM);
+        type[iatom] = atom->lmap->find_type(typestr, Atom::ATOM);
         if (type[iatom] == -1)
           error->all(FLERR, fileiarg, "Unknown atom type {} in {}: {}", typestr, location,
                      utils::trim(line));
@@ -3050,7 +3050,7 @@ void Molecule::bonds(int flag, char *line)
       case 1: {    // type label
         if (!atom->labelmapflag)
           error->all(FLERR, fileiarg, "Invalid bond type {} in {}: {}", typestr, location, utils::trim(line));
-        itype = atom->lmap->find(typestr, Atom::BOND);
+        itype = atom->lmap->find_type(typestr, Atom::BOND);
         if (itype == -1)
           error->all(FLERR, fileiarg, "Unknown bond type {} in {}: {}", typestr, location, utils::trim(line));
         break;
@@ -3136,7 +3136,7 @@ void Molecule::angles(int flag, char *line)
       case 1: {    // type label
         if (!atom->labelmapflag)
           error->all(FLERR, fileiarg, "Invalid angle type {} in {}: {}", typestr, location, utils::trim(line));
-        itype = atom->lmap->find(typestr, Atom::ANGLE);
+        itype = atom->lmap->find_type(typestr, Atom::ANGLE);
         if (itype == -1)
           error->all(FLERR, fileiarg, "Unknown angle type {} in {}: {}", typestr, location, utils::trim(line));
         break;
@@ -3237,7 +3237,7 @@ void Molecule::dihedrals(int flag, char *line)
       case 1: {    // type label
         if (!atom->labelmapflag)
           error->all(FLERR, fileiarg, "Invalid dihedral type {} in {}: {}", typestr, location, utils::trim(line));
-        itype = atom->lmap->find(typestr, Atom::DIHEDRAL);
+        itype = atom->lmap->find_type(typestr, Atom::DIHEDRAL);
         if (itype == -1)
           error->all(FLERR, fileiarg, "Unknown dihedral type {} in {}: {}", typestr, location, utils::trim(line));
         break;
@@ -3352,7 +3352,7 @@ void Molecule::impropers(int flag, char *line)
       case 1: {    // type label
         if (!atom->labelmapflag)
           error->all(FLERR, fileiarg, "Invalid improper type {} in {}: {}", typestr, location, utils::trim(line));
-        itype = atom->lmap->find(typestr, Atom::IMPROPER);
+        itype = atom->lmap->find_type(typestr, Atom::IMPROPER);
         if (itype == -1)
           error->all(FLERR, fileiarg, "Unknown improper type {} in {}: {}", typestr, location, utils::trim(line));
         break;
@@ -4274,7 +4274,7 @@ void Molecule::print(FILE *fp)
     fputs("\nTypes\n\n", fp);
     if (atom->labelmapflag && atom->lmap->is_complete(Atom::ATOM)) {
       for (int i = 0; i < natoms; i++)
-        utils::print(fp, " {} {}\n", i + 1, atom->lmap->find(type[i], Atom::ATOM));
+        utils::print(fp, " {} {}\n", i + 1, atom->lmap->find_label(type[i], Atom::ATOM));
     } else {
       for (int i = 0; i < natoms; i++) utils::print(fp, " {}  {}\n", i + 1, type[i]);
     }
@@ -4329,7 +4329,7 @@ void Molecule::print(FILE *fp)
         if (has_newton_bond || (i + 1 < bond_atom[i][j])) {
           ++idx;
           if (has_typelabels) {
-            utils::print(fp, " {}  {}", idx, atom->lmap->find(bond_type[i][j], Atom::BOND));
+            utils::print(fp, " {}  {}", idx, atom->lmap->find_label(bond_type[i][j], Atom::BOND));
           } else {
             utils::print(fp, " {}  {}", idx, bond_type[i][j]);
           }
@@ -4348,7 +4348,7 @@ void Molecule::print(FILE *fp)
         if (has_newton_bond || (i + 1 == angle_atom2[i][j])) {
           ++idx;
           if (has_typelabels) {
-            utils::print(fp, " {}  {}", idx, atom->lmap->find(angle_type[i][j], Atom::ANGLE));
+            utils::print(fp, " {}  {}", idx, atom->lmap->find_label(angle_type[i][j], Atom::ANGLE));
           } else {
             utils::print(fp, " {}  {}", idx, angle_type[i][j]);
           }
@@ -4367,7 +4367,7 @@ void Molecule::print(FILE *fp)
         if (has_newton_bond || (i + 1 == dihedral_atom2[i][j])) {
           ++idx;
           if (has_typelabels) {
-            utils::print(fp, " {}  {}", idx, atom->lmap->find(dihedral_type[i][j], Atom::DIHEDRAL));
+            utils::print(fp, " {}  {}", idx, atom->lmap->find_label(dihedral_type[i][j], Atom::DIHEDRAL));
           } else {
             utils::print(fp, " {}  {}", idx, dihedral_type[i][j]);
           }
@@ -4387,7 +4387,7 @@ void Molecule::print(FILE *fp)
         if (has_newton_bond || (i + 1 == improper_atom2[i][j])) {
           ++idx;
           if (has_typelabels) {
-            utils::print(fp, " {}  {}", idx, atom->lmap->find(improper_type[i][j], Atom::IMPROPER));
+            utils::print(fp, " {}  {}", idx, atom->lmap->find_label(improper_type[i][j], Atom::IMPROPER));
           } else {
             utils::print(fp, " {}  {}", idx, improper_type[i][j]);
           }
@@ -4452,9 +4452,9 @@ void Molecule::print(FILE *fp)
         case 1:
           has_typelabels = has_typelabels && atom->lmap->is_complete(Atom::ANGLE);
           if (has_typelabels) {
-            utils::print(fp, " {} {} {}\n", atom->lmap->find(shake_type[i][0], Atom::BOND),
-                         atom->lmap->find(shake_type[i][1], Atom::BOND),
-                         atom->lmap->find(shake_type[i][2], Atom::ANGLE));
+            utils::print(fp, " {} {} {}\n", atom->lmap->find_label(shake_type[i][0], Atom::BOND),
+                         atom->lmap->find_label(shake_type[i][1], Atom::BOND),
+                         atom->lmap->find_label(shake_type[i][2], Atom::ANGLE));
           } else {
             utils::print(fp, " {} {} {}\n", shake_type[i][0], shake_type[i][1], shake_type[i][2]);
           }
@@ -4462,8 +4462,8 @@ void Molecule::print(FILE *fp)
 
         case 2:
           if (has_typelabels) {
-            utils::print(fp, " {} {}\n", atom->lmap->find(shake_type[i][0], Atom::BOND),
-                         atom->lmap->find(shake_type[i][1], Atom::BOND));
+            utils::print(fp, " {} {}\n", atom->lmap->find_label(shake_type[i][0], Atom::BOND),
+                         atom->lmap->find_label(shake_type[i][1], Atom::BOND));
           } else {
             utils::print(fp, " {} {}\n", shake_type[i][0], shake_type[i][1]);
           }
@@ -4471,9 +4471,9 @@ void Molecule::print(FILE *fp)
 
         case 3:
           if (has_typelabels) {
-            utils::print(fp, " {} {} {}\n", atom->lmap->find(shake_type[i][0], Atom::BOND),
-                         atom->lmap->find(shake_type[i][1], Atom::BOND),
-                         atom->lmap->find(shake_type[i][2], Atom::BOND));
+            utils::print(fp, " {} {} {}\n", atom->lmap->find_label(shake_type[i][0], Atom::BOND),
+                         atom->lmap->find_label(shake_type[i][1], Atom::BOND),
+                         atom->lmap->find_label(shake_type[i][2], Atom::BOND));
           } else {
             utils::print(fp, " {} {} {}\n", shake_type[i][0], shake_type[i][1], shake_type[i][2]);
           }
@@ -4481,10 +4481,10 @@ void Molecule::print(FILE *fp)
 
         case 4:
           if (has_typelabels) {
-            utils::print(fp, " {} {} {} {}\n", atom->lmap->find(shake_type[i][0], Atom::BOND),
-                         atom->lmap->find(shake_type[i][1], Atom::BOND),
-                         atom->lmap->find(shake_type[i][2], Atom::BOND),
-                         atom->lmap->find(shake_type[i][3], Atom::BOND));
+            utils::print(fp, " {} {} {} {}\n", atom->lmap->find_label(shake_type[i][0], Atom::BOND),
+                         atom->lmap->find_label(shake_type[i][1], Atom::BOND),
+                         atom->lmap->find_label(shake_type[i][2], Atom::BOND),
+                         atom->lmap->find_label(shake_type[i][3], Atom::BOND));
           } else {
             utils::print(fp, " {} {} {} {}\n", shake_type[i][0], shake_type[i][1], shake_type[i][2],
                          shake_type[i][3]);
diff --git a/src/output.cpp b/src/output.cpp
index f7ac581b4dd..9fa7ecff494 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -745,7 +745,7 @@ void Output::write_molecule_json(FILE *fp, int json_level, int printflag, int *i
         for (auto myatom : atoms_root) {
           int mytype = myatom.type;
           std::string typestr = std::to_string(mytype);
-          if (atom->labelmapflag) typestr = atom->lmap->find(mytype, Atom::ATOM);
+          if (atom->labelmapflag) typestr = atom->lmap->find_label(mytype, Atom::ATOM);
           utils::print(fp, "{}[{}, \"{}\"]", indent, myatom.tag, typestr);
           if (std::next(it) == atoms_root.end()) fprintf(fp, "\n");
           else fprintf(fp, ",\n");
diff --git a/src/utils.cpp b/src/utils.cpp
index 45e65c4b328..038a34b9cf5 100644
--- a/src/utils.cpp
+++ b/src/utils.cpp
@@ -1150,7 +1150,7 @@ char *utils::expand_type(const char *file, int line, const std::string &str, int
       lmp->error->all(file, line, "{} type string {} cannot be used without a labelmap",
                       labeltypes[mode], typestr);
 
-    int type = lmp->atom->lmap->find(typestr, mode);
+    int type = lmp->atom->lmap->find_type(typestr, mode);
     if (type == -1)
       lmp->error->all(file, line, "{} type string {} not found in labelmap", labeltypes[mode],
                       typestr);
diff --git a/src/variable.cpp b/src/variable.cpp
index afe3b840204..5b0492dd1c8 100644
--- a/src/variable.cpp
+++ b/src/variable.cpp
@@ -4543,15 +4543,15 @@ int Variable::special_function(const std::string &word, char *contents, Tree **t
 
     int value = -1;
     if (kind == "atom") {
-      value = atom->lmap->find(typestr,Atom::ATOM);
+      value = atom->lmap->find_type(typestr,Atom::ATOM);
     } else if (kind == "bond") {
-      value = atom->lmap->find(typestr,Atom::BOND);
+      value = atom->lmap->find_type(typestr,Atom::BOND);
     } else if (kind == "angle") {
-      value = atom->lmap->find(typestr,Atom::ANGLE);
+      value = atom->lmap->find_type(typestr,Atom::ANGLE);
     } else if (kind == "dihedral") {
-      value = atom->lmap->find(typestr,Atom::DIHEDRAL);
+      value = atom->lmap->find_type(typestr,Atom::DIHEDRAL);
     } else if (kind == "improper") {
-      value = atom->lmap->find(typestr,Atom::IMPROPER);
+      value = atom->lmap->find_type(typestr,Atom::IMPROPER);
     } else {
       print_var_error(FLERR, fmt::format("Invalid kind {} in {}() in variable", kind, word),ivar);
     }
diff --git a/unittest/commands/test_labelmap.cpp b/unittest/commands/test_labelmap.cpp
index 1025b46e84f..cab34ea61fd 100644
--- a/unittest/commands/test_labelmap.cpp
+++ b/unittest/commands/test_labelmap.cpp
@@ -92,11 +92,18 @@ TEST_F(LabelMapTest, Atoms)
     command("mass \"O#\" 10.0");
     END_HIDE_OUTPUT();
     EXPECT_TRUE(atom->lmap->is_complete(Atom::ATOM));
-    EXPECT_EQ(atom->lmap->find("C1", Atom::ATOM), 1);
-    EXPECT_EQ(atom->lmap->find("N2", Atom::ATOM), 2);
-    EXPECT_EQ(atom->lmap->find("O#", Atom::ATOM), 3);
-    EXPECT_EQ(atom->lmap->find("H", Atom::ATOM), 4);
-    EXPECT_EQ(atom->lmap->find("X", Atom::ATOM), -1);
+    EXPECT_EQ(atom->lmap->find_type("C1", Atom::ATOM), 1);
+    EXPECT_EQ(atom->lmap->find_type("N2", Atom::ATOM), 2);
+    EXPECT_EQ(atom->lmap->find_type("O#", Atom::ATOM), 3);
+    EXPECT_EQ(atom->lmap->find_type("H", Atom::ATOM), 4);
+    EXPECT_EQ(atom->lmap->find_type("X", Atom::ATOM), -1);
+    EXPECT_EQ(atom->lmap->find_type("", Atom::ATOM), -1);
+    EXPECT_THAT(atom->lmap->find_label(1, Atom::ATOM), StrEq("C1"));
+    EXPECT_THAT(atom->lmap->find_label(2, Atom::ATOM), StrEq("N2"));
+    EXPECT_THAT(atom->lmap->find_label(3, Atom::ATOM), StrEq("O#"));
+    EXPECT_THAT(atom->lmap->find_label(4, Atom::ATOM), StrEq("H"));
+    EXPECT_THAT(atom->lmap->find_label(-1, Atom::ATOM), StrEq(""));
+    EXPECT_THAT(atom->lmap->find_label(5, Atom::ATOM), StrEq(""));
     EXPECT_DOUBLE_EQ(atom->mass[3], 10.0);
 
     EXPECT_EQ(utils::expand_type(FLERR, "1", Atom::ATOM, lmp), nullptr);
@@ -161,18 +168,22 @@ TEST_F(LabelMapTest, Atoms)
     command("labelmap atom 3 C1 2 N2");
     END_HIDE_OUTPUT();
     EXPECT_FALSE(atom->lmap->is_complete(Atom::ATOM));
-    EXPECT_EQ(atom->lmap->find("C1", Atom::ATOM), 3);
-    EXPECT_EQ(atom->lmap->find("N2", Atom::ATOM), 2);
+    EXPECT_EQ(atom->lmap->find_type("C1", Atom::ATOM), 3);
+    EXPECT_EQ(atom->lmap->find_type("N2", Atom::ATOM), 2);
 
     BEGIN_HIDE_OUTPUT();
     command("labelmap clear");
     command(R"(labelmap atom 1 "C1'" 2 'C2"' 3 """C1'-C2" """ 4 """ C2"-C1'""")");
     END_HIDE_OUTPUT();
     EXPECT_TRUE(atom->lmap->is_complete(Atom::ATOM));
-    EXPECT_EQ(atom->lmap->find("C1'", Atom::ATOM), 1);
-    EXPECT_EQ(atom->lmap->find("C2\"", Atom::ATOM), 2);
-    EXPECT_EQ(atom->lmap->find("C1'-C2\"", Atom::ATOM), 3);
-    EXPECT_EQ(atom->lmap->find("C2\"-C1'", Atom::ATOM), 4);
+    EXPECT_EQ(atom->lmap->find_type("C1'", Atom::ATOM), 1);
+    EXPECT_EQ(atom->lmap->find_type(R"(C2")", Atom::ATOM), 2);
+    EXPECT_EQ(atom->lmap->find_type(R"(C1'-C2")", Atom::ATOM), 3);
+    EXPECT_EQ(atom->lmap->find_type(R"(C2"-C1')", Atom::ATOM), 4);
+    EXPECT_THAT(atom->lmap->find_label(1, Atom::ATOM), StrEq("C1'"));
+    EXPECT_THAT(atom->lmap->find_label(2, Atom::ATOM), StrEq(R"(C2")"));
+    EXPECT_THAT(atom->lmap->find_label(3, Atom::ATOM), StrEq(R"(C1'-C2")"));
+    EXPECT_THAT(atom->lmap->find_label(4, Atom::ATOM), StrEq(R"(C2"-C1')"));
 }
 
 TEST_F(LabelMapTest, Topology)
@@ -220,17 +231,33 @@ TEST_F(LabelMapTest, Topology)
     EXPECT_TRUE(atom->lmap->is_complete(Atom::ANGLE));
     EXPECT_TRUE(atom->lmap->is_complete(Atom::DIHEDRAL));
     EXPECT_TRUE(atom->lmap->is_complete(Atom::IMPROPER));
-    EXPECT_EQ(atom->lmap->find("C1", Atom::ATOM), 1);
-    EXPECT_EQ(atom->lmap->find("N2'", Atom::ATOM), 2);
-    EXPECT_EQ(atom->lmap->find("C1-N2", Atom::BOND), 1);
-    EXPECT_EQ(atom->lmap->find("[C1][C1]", Atom::BOND), 2);
-    EXPECT_EQ(atom->lmap->find("N2=N2", Atom::BOND), 3);
-    EXPECT_EQ(atom->lmap->find("C1-N2-C1", Atom::ANGLE), 1);
-    EXPECT_EQ(atom->lmap->find("N2'-C1\"-N2'", Atom::ANGLE), 2);
-    EXPECT_EQ(atom->lmap->find("C1-N2-C1-N2", Atom::DIHEDRAL), 1);
-    EXPECT_EQ(atom->lmap->find("C1-N2-C1-N2", Atom::IMPROPER), 1);
-    EXPECT_EQ(atom->lmap->find("X", Atom::ATOM), -1);
-    EXPECT_EQ(atom->lmap->find("N2'-C1\"-N2'", Atom::BOND), -1);
+    EXPECT_EQ(atom->lmap->find_type("C1", Atom::ATOM), 1);
+    EXPECT_EQ(atom->lmap->find_type("N2'", Atom::ATOM), 2);
+    EXPECT_EQ(atom->lmap->find_type("C1-N2", Atom::BOND), 1);
+    EXPECT_EQ(atom->lmap->find_type("[C1][C1]", Atom::BOND), 2);
+    EXPECT_EQ(atom->lmap->find_type("N2=N2", Atom::BOND), 3);
+    EXPECT_EQ(atom->lmap->find_type("C1-N2-C1", Atom::ANGLE), 1);
+    EXPECT_EQ(atom->lmap->find_type("N2'-C1\"-N2'", Atom::ANGLE), 2);
+    EXPECT_EQ(atom->lmap->find_type("C1-N2-C1-N2", Atom::DIHEDRAL), 1);
+    EXPECT_EQ(atom->lmap->find_type("C1-N2-C1-N2", Atom::IMPROPER), 1);
+    EXPECT_EQ(atom->lmap->find_type("X", Atom::ATOM), -1);
+    EXPECT_EQ(atom->lmap->find_type("N2'-C1\"-N2'", Atom::BOND), -1);
+
+    EXPECT_THAT(atom->lmap->find_label(1, Atom::BOND), StrEq("C1-N2"));
+    EXPECT_THAT(atom->lmap->find_label(2, Atom::BOND), StrEq("[C1][C1]"));
+    EXPECT_THAT(atom->lmap->find_label(3, Atom::BOND), StrEq("N2=N2"));
+    EXPECT_THAT(atom->lmap->find_label(1, Atom::ANGLE), StrEq("C1-N2-C1"));
+    EXPECT_THAT(atom->lmap->find_label(2, Atom::ANGLE), StrEq(R"(N2'-C1"-N2')"));
+    EXPECT_THAT(atom->lmap->find_label(1, Atom::DIHEDRAL), StrEq("C1-N2-C1-N2"));
+    EXPECT_THAT(atom->lmap->find_label(1, Atom::IMPROPER), StrEq("C1-N2-C1-N2"));
+    EXPECT_THAT(atom->lmap->find_label(0, Atom::BOND), StrEq(""));
+    EXPECT_THAT(atom->lmap->find_label(4, Atom::BOND), StrEq(""));
+    EXPECT_THAT(atom->lmap->find_label(-1, Atom::ANGLE), StrEq(""));
+    EXPECT_THAT(atom->lmap->find_label(3, Atom::ANGLE), StrEq(""));
+    EXPECT_THAT(atom->lmap->find_label(0, Atom::DIHEDRAL), StrEq(""));
+    EXPECT_THAT(atom->lmap->find_label(2, Atom::DIHEDRAL), StrEq(""));
+    EXPECT_THAT(atom->lmap->find_label(-1, Atom::IMPROPER), StrEq(""));
+    EXPECT_THAT(atom->lmap->find_label(0, Atom::IMPROPER), StrEq(""));
     EXPECT_DOUBLE_EQ(atom->mass[1], 12.0);
     EXPECT_DOUBLE_EQ(atom->mass[2], 14.0);
 
@@ -255,17 +282,17 @@ TEST_F(LabelMapTest, Topology)
     EXPECT_TRUE(atom->lmap->is_complete(Atom::ANGLE));
     EXPECT_TRUE(atom->lmap->is_complete(Atom::DIHEDRAL));
     EXPECT_TRUE(atom->lmap->is_complete(Atom::IMPROPER));
-    EXPECT_EQ(atom->lmap->find("C1", Atom::ATOM), 1);
-    EXPECT_EQ(atom->lmap->find("N2'", Atom::ATOM), 2);
-    EXPECT_EQ(atom->lmap->find("C1-N2", Atom::BOND), 1);
-    EXPECT_EQ(atom->lmap->find("[C1][C1]", Atom::BOND), 2);
-    EXPECT_EQ(atom->lmap->find("N2=N2", Atom::BOND), 3);
-    EXPECT_EQ(atom->lmap->find("C1-N2-C1", Atom::ANGLE), 1);
-    EXPECT_EQ(atom->lmap->find("N2'-C1\"-N2'", Atom::ANGLE), 2);
-    EXPECT_EQ(atom->lmap->find("C1-N2-C1-N2", Atom::DIHEDRAL), 1);
-    EXPECT_EQ(atom->lmap->find("C1-N2-C1-N2", Atom::IMPROPER), 1);
-    EXPECT_EQ(atom->lmap->find("X", Atom::ATOM), -1);
-    EXPECT_EQ(atom->lmap->find("N2'-C1\"-N2'", Atom::BOND), -1);
+    EXPECT_EQ(atom->lmap->find_type("C1", Atom::ATOM), 1);
+    EXPECT_EQ(atom->lmap->find_type("N2'", Atom::ATOM), 2);
+    EXPECT_EQ(atom->lmap->find_type("C1-N2", Atom::BOND), 1);
+    EXPECT_EQ(atom->lmap->find_type("[C1][C1]", Atom::BOND), 2);
+    EXPECT_EQ(atom->lmap->find_type("N2=N2", Atom::BOND), 3);
+    EXPECT_EQ(atom->lmap->find_type("C1-N2-C1", Atom::ANGLE), 1);
+    EXPECT_EQ(atom->lmap->find_type("N2'-C1\"-N2'", Atom::ANGLE), 2);
+    EXPECT_EQ(atom->lmap->find_type("C1-N2-C1-N2", Atom::DIHEDRAL), 1);
+    EXPECT_EQ(atom->lmap->find_type("C1-N2-C1-N2", Atom::IMPROPER), 1);
+    EXPECT_EQ(atom->lmap->find_type("X", Atom::ATOM), -1);
+    EXPECT_EQ(atom->lmap->find_type("N2'-C1\"-N2'", Atom::BOND), -1);
     platform::unlink("labelmap_topology.inc");
 
     auto *expanded = utils::expand_type(FLERR, "N2'", Atom::ATOM, lmp);
@@ -291,6 +318,57 @@ TEST_F(LabelMapTest, Topology)
                  utils::expand_type(FLERR, "XX", Atom::DIHEDRAL, lmp););
     TEST_FAILURE(".*ERROR: Improper type string XX not found in labelmap.*",
                  utils::expand_type(FLERR, "XX", Atom::IMPROPER, lmp););
+
+    // check for inference
+
+    BEGIN_HIDE_OUTPUT();
+    command(R"(labelmap atom 2 "N2")");
+    command(R"(labelmap bond 2 "C1-C1")");
+    command(R"(labelmap bond 3 "N2-N2")");
+    command(R"(labelmap angle 2 "N2-C1-N2")");
+    END_HIDE_OUTPUT();
+
+    EXPECT_EQ(atom->lmap->infer_bondtype(1, 1), 2);
+    EXPECT_EQ(atom->lmap->infer_bondtype(1, 2), 1);
+    EXPECT_EQ(atom->lmap->infer_bondtype(2, 2), 3);
+    EXPECT_EQ(atom->lmap->infer_bondtype(2, 3), -1);
+    EXPECT_EQ(atom->lmap->infer_bondtype({"N2", "N2"}), 3);
+    EXPECT_EQ(atom->lmap->infer_bondtype({"C1", "N2"}), 1);
+    EXPECT_EQ(atom->lmap->infer_bondtype({"C1", "C1"}), 2);
+    EXPECT_EQ(atom->lmap->infer_bondtype({"N2", "C1"}), 1);
+    EXPECT_EQ(atom->lmap->infer_bondtype({"N3", "C1"}), -1);
+    EXPECT_EQ(atom->lmap->infer_angletype(1, 2, 1), 1);
+    EXPECT_EQ(atom->lmap->infer_angletype(2, 1, 2), 2);
+    EXPECT_EQ(atom->lmap->infer_angletype(1, 1, 1), -1);
+    EXPECT_EQ(atom->lmap->infer_angletype(3, 1, 1), -1);
+    EXPECT_EQ(atom->lmap->infer_angletype(1, 3, 1), -1);
+    EXPECT_EQ(atom->lmap->infer_angletype(1, 1, 3), -1);
+    EXPECT_EQ(atom->lmap->infer_angletype({"C1", "N2", "C1"}), 1);
+    EXPECT_EQ(atom->lmap->infer_angletype({"N2", "C1", "N2"}), 2);
+    EXPECT_EQ(atom->lmap->infer_angletype({"C1", "N2", "N2"}), -1);
+    EXPECT_EQ(atom->lmap->infer_angletype({"C1", "N3", "C1"}), -1);
+    EXPECT_EQ(atom->lmap->infer_dihedraltype(1, 2, 1, 2), 1);
+    EXPECT_EQ(atom->lmap->infer_dihedraltype(2, 1, 2, 2), -1);
+    EXPECT_EQ(atom->lmap->infer_dihedraltype(2, 1, 2, 1), 1);
+    EXPECT_EQ(atom->lmap->infer_dihedraltype(3, 1, 2, 1), -1);
+    EXPECT_EQ(atom->lmap->infer_dihedraltype(2, 3, 2, 1), -1);
+    EXPECT_EQ(atom->lmap->infer_dihedraltype(2, 1, 3, 1), -1);
+    EXPECT_EQ(atom->lmap->infer_dihedraltype(2, 1, 2, 3), -1);
+    EXPECT_EQ(atom->lmap->infer_dihedraltype({"C1", "N2", "C1", "N2"}), 1);
+    EXPECT_EQ(atom->lmap->infer_dihedraltype({"N2", "C1", "N2", "C1"}), 1);
+    EXPECT_EQ(atom->lmap->infer_dihedraltype({"C1", "N2", "N2", "C1"}), -1);
+    EXPECT_EQ(atom->lmap->infer_dihedraltype({"N3", "C1", "N2", "C1"}), -1);
+    EXPECT_EQ(atom->lmap->infer_impropertype(1, 2, 1, 2), 1);
+    EXPECT_EQ(atom->lmap->infer_impropertype(2, 1, 2, 2), -1);
+    EXPECT_EQ(atom->lmap->infer_impropertype(2, 1, 2, 1), 1);
+    EXPECT_EQ(atom->lmap->infer_impropertype(3, 1, 2, 1), -1);
+    EXPECT_EQ(atom->lmap->infer_impropertype(2, 3, 2, 1), -1);
+    EXPECT_EQ(atom->lmap->infer_impropertype(2, 1, 3, 1), -1);
+    EXPECT_EQ(atom->lmap->infer_impropertype(2, 1, 2, 3), -1);
+    EXPECT_EQ(atom->lmap->infer_impropertype({"C1", "N2", "C1", "N2"}), 1);
+    EXPECT_EQ(atom->lmap->infer_impropertype({"N2", "C1", "N2", "C1"}), 1);
+    EXPECT_EQ(atom->lmap->infer_impropertype({"C1", "N2", "N2", "C1"}), 1);
+    EXPECT_EQ(atom->lmap->infer_impropertype({"N3", "C1", "N2", "C1"}), -1);
 }
 } // namespace LAMMPS_NS