Merge pull request #4205 from akohlmey/next-release

Update version tags for next feature release

.github/CODEOWNERS

@@ -38,6 +38,7 @@ src/ML-HDNNP/*        @singraber
 src/ML-IAP/*          @athomps
 src/ML-PACE/*         @yury-lysogorskiy
 src/ML-POD/*          @exapde
+src/ML-UF3/*          @monk-04
 src/MOFFF/*           @hheenen
 src/MOLFILE/*         @akohlmey
 src/NETCDF/*          @pastewka
@@ -58,7 +59,8 @@ src/VTK/*             @rbberger
 
 # individual files in packages
 src/GPU/pair_vashishta_gpu.*        @andeplane
-src/KOKKOS/pair_vashishta_kokkos.*  @andeplane
+src/KOKKOS/pair_vashishta_kokkos.*  @andeplane @stanmoore1
+src/KOSSOS/pair_pod_kokkos.*        @exapde @stanmoore1
 src/MANYBODY/pair_vashishta_table.* @andeplane
 src/MANYBODY/pair_atm.*             @sergeylishchuk
 src/MANYBODY/pair_nb3b_screened.*   @flodesani
@@ -72,6 +74,8 @@ src/MC/fix_sgcmc.*                  @athomps
 src/REAXFF/compute_reaxff_atom.*    @rbberger
 src/KOKKOS/compute_reaxff_atom_kokkos.*    @rbberger
 src/REPLICA/fix_pimd_langevin.*     @Yi-FanLi
+src/DPD-BASIC/pair_dpd_coul_slater_long.* @Eddy-Barraud
+src/GPU/pair_dpd_coul_slater_long.* @Eddy-Barraud
 
 # core LAMMPS classes
 src/lammps.*              @sjplimp

.github/workflows/unittest-macos.yml

@@ -15,7 +15,7 @@ jobs:
   build:
     name: MacOS Unit Test
     if: ${{ github.repository == 'lammps/lammps' }}
-    runs-on: macos-latest
+    runs-on: macos-13
     env:
       CCACHE_DIR: ${{ github.workspace }}/.ccache
 
@@ -43,6 +43,8 @@ jobs:
       working-directory: build
       run: |
         ccache -z
+        python3 -m venv macosenv
+        source macosenv/bin/activate
         python3 -m pip install numpy
         python3 -m pip install pyyaml
         cmake -C ../cmake/presets/clang.cmake \

.gitignore

@@ -43,12 +43,12 @@ Thumbs.db
 
 #cmake
 /build*
-/CMakeCache.txt
-/CMakeFiles/
-/Testing
+CMakeCache.txt
+CMakeFiles
 /Makefile
-/Testing
-/cmake_install.cmake
+Testing
+Temporary
+cmake_install.cmake
 /lmp
 out/Debug
 out/RelWithDebInfo
@@ -60,3 +60,4 @@ src/Makefile.package.settings-e
 /cmake/build/x64-Debug-Clang
 /install/x64-GUI-MSVC
 /install
+.Rhistory

cmake/CMakeLists.txt

@@ -23,6 +23,7 @@ project(lammps CXX)
 set(SOVERSION 0)
 get_property(BUILD_IS_MULTI_CONFIG GLOBAL PROPERTY GENERATOR_IS_MULTI_CONFIG)
 
+include(GNUInstallDirs)
 get_filename_component(LAMMPS_DIR ${CMAKE_CURRENT_SOURCE_DIR}/.. ABSOLUTE)
 get_filename_component(LAMMPS_LIB_BINARY_DIR ${CMAKE_BINARY_DIR}/lib ABSOLUTE)
 # collect all executables and shared libs in the top level build folder
@@ -208,7 +209,7 @@ else()
   unset(CMAKE_CXX_CLANG_TIDY CACHE)
 endif()
 
-include(GNUInstallDirs)
+
 file(GLOB ALL_SOURCES CONFIGURE_DEPENDS ${LAMMPS_SOURCE_DIR}/[^.]*.cpp)
 file(GLOB MAIN_SOURCES CONFIGURE_DEPENDS ${LAMMPS_SOURCE_DIR}/main.cpp)
 list(REMOVE_ITEM ALL_SOURCES ${MAIN_SOURCES})
@@ -256,6 +257,7 @@ set(STANDARD_PACKAGES
   DRUDE
   EFF
   ELECTRODE
+  EXTRA-COMMAND
   EXTRA-COMPUTE
   EXTRA-DUMP
   EXTRA-FIX
@@ -281,10 +283,11 @@ set(STANDARD_PACKAGES
   ML-HDNNP
   ML-IAP
   ML-PACE
+  ML-POD
   ML-QUIP
   ML-RANN
   ML-SNAP
-  ML-POD
+  ML-UF3
   MOFFF
   MOLECULE
   MOLFILE
@@ -689,7 +692,7 @@ endif()
 # packages which selectively include variants based on enabled styles
 # e.g. accelerator packages
 ######################################################################
-foreach(PKG_WITH_INCL CORESHELL DPD-SMOOTH MC MISC PHONON QEQ OPENMP KOKKOS OPT INTEL GPU)
+foreach(PKG_WITH_INCL CORESHELL DPD-BASIC DPD-SMOOTH MC MISC PHONON QEQ OPENMP KOKKOS OPT INTEL GPU)
   if(PKG_${PKG_WITH_INCL})
     include(Packages/${PKG_WITH_INCL})
   endif()

cmake/Modules/OpenCLLoader.cmake

@@ -1,6 +1,6 @@
 message(STATUS "Downloading and building OpenCL loader library")
-set(OPENCL_LOADER_URL "${LAMMPS_THIRDPARTY_URL}/opencl-loader-2024.02.09.tar.gz" CACHE STRING "URL for OpenCL loader tarball")
-set(OPENCL_LOADER_MD5 "f3573cf9daa3558ba46fd5866517f38f" CACHE STRING "MD5 checksum of OpenCL loader tarball")
+set(OPENCL_LOADER_URL "${LAMMPS_THIRDPARTY_URL}/opencl-loader-2024.05.09.tar.gz" CACHE STRING "URL for OpenCL loader tarball")
+set(OPENCL_LOADER_MD5 "e7796826b21c059224fabe997e0f2075" CACHE STRING "MD5 checksum of OpenCL loader tarball")
 mark_as_advanced(OPENCL_LOADER_URL)
 mark_as_advanced(OPENCL_LOADER_MD5)
 

cmake/Modules/Packages/DPD-BASIC.cmake

@@ -0,0 +1,9 @@
+# pair style dpd/coul/slater/long may only be installed if also KSPACE is installed
+if(NOT PKG_KSPACE)
+  get_property(LAMMPS_PAIR_HEADERS GLOBAL PROPERTY PAIR)
+  list(REMOVE_ITEM LAMMPS_PAIR_HEADERS ${LAMMPS_SOURCE_DIR}/DPD-BASIC/pair_dpd_coul_slater_long.h)
+  set_property(GLOBAL PROPERTY PAIR "${LAMMPS_PAIR_HEADERS}")
+  get_target_property(LAMMPS_SOURCES lammps SOURCES)
+  list(REMOVE_ITEM LAMMPS_SOURCES ${LAMMPS_SOURCE_DIR}/DPD-BASIC/pair_dpd_coul_slater_long.cpp)
+  set_property(TARGET lammps PROPERTY SOURCES "${LAMMPS_SOURCES}")
+endif()

cmake/Modules/Packages/KOKKOS.cmake

@@ -45,8 +45,8 @@ if(DOWNLOAD_KOKKOS)
   list(APPEND KOKKOS_LIB_BUILD_ARGS "-DCMAKE_CXX_EXTENSIONS=${CMAKE_CXX_EXTENSIONS}")
   list(APPEND KOKKOS_LIB_BUILD_ARGS "-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TOOLCHAIN_FILE}")
   include(ExternalProject)
-  set(KOKKOS_URL "https://github.com/kokkos/kokkos/archive/4.3.00.tar.gz" CACHE STRING "URL for KOKKOS tarball")
-  set(KOKKOS_MD5 "889dcea2b5ced3debdc5b0820044bdc4" CACHE STRING "MD5 checksum of KOKKOS tarball")
+  set(KOKKOS_URL "https://github.com/kokkos/kokkos/archive/4.3.01.tar.gz" CACHE STRING "URL for KOKKOS tarball")
+  set(KOKKOS_MD5 "243de871b3dc2cf3990c1c404032df83" CACHE STRING "MD5 checksum of KOKKOS tarball")
   mark_as_advanced(KOKKOS_URL)
   mark_as_advanced(KOKKOS_MD5)
   GetFallbackURL(KOKKOS_URL KOKKOS_FALLBACK)
@@ -71,7 +71,7 @@ if(DOWNLOAD_KOKKOS)
   add_dependencies(LAMMPS::KOKKOSCORE kokkos_build)
   add_dependencies(LAMMPS::KOKKOSCONTAINERS kokkos_build)
 elseif(EXTERNAL_KOKKOS)
-  find_package(Kokkos 4.3.00 REQUIRED CONFIG)
+  find_package(Kokkos 4.3.01 REQUIRED CONFIG)
   target_link_libraries(lammps PRIVATE Kokkos::kokkos)
 else()
   set(LAMMPS_LIB_KOKKOS_SRC_DIR ${LAMMPS_LIB_SOURCE_DIR}/kokkos)

cmake/Modules/Packages/PLUMED.cmake

@@ -1,5 +1,9 @@
 # Plumed2 support for PLUMED package
 
+# for supporting multiple concurrent plumed2 installations for debugging and testing
+set(PLUMED_SUFFIX "" CACHE STRING "Suffix for Plumed2 library")
+mark_as_advanced(PLUMED_SUFFIX)
+
 if(BUILD_MPI)
   set(PLUMED_CONFIG_MPI "--enable-mpi")
   set(PLUMED_CONFIG_CC  ${CMAKE_MPI_C_COMPILER})
@@ -21,9 +25,11 @@ else()
   set(PLUMED_CONFIG_OMP "--disable-openmp")
 endif()
 
-set(PLUMED_URL "https://github.com/plumed/plumed2/releases/download/v2.8.3/plumed-src-2.8.3.tgz"
+# Note: must also adjust check for supported API versions in
+# fix_plumed.cpp when version changes from v2.n.x to v2.n+1.y
+set(PLUMED_URL "https://github.com/plumed/plumed2/releases/download/v2.9.1/plumed-src-2.9.1.tgz"
   CACHE STRING "URL for PLUMED tarball")
-set(PLUMED_MD5 "76d23cd394eba9e6530316ed1184e219" CACHE STRING "MD5 checksum of PLUMED tarball")
+set(PLUMED_MD5 "c3b2d31479c1e9ce211719d40e9efbd7" CACHE STRING "MD5 checksum of PLUMED tarball")
 
 mark_as_advanced(PLUMED_URL)
 mark_as_advanced(PLUMED_MD5)
@@ -151,15 +157,15 @@ else()
     file(MAKE_DIRECTORY ${INSTALL_DIR}/include)
   else()
     find_package(PkgConfig REQUIRED)
-    pkg_check_modules(PLUMED REQUIRED plumed)
+    pkg_check_modules(PLUMED REQUIRED plumed${PLUMED_SUFFIX})
     add_library(LAMMPS::PLUMED INTERFACE IMPORTED)
     if(PLUMED_MODE STREQUAL "STATIC")
-      include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.static)
+      include(${PLUMED_LIBDIR}/plumed${PLUMED_SUFFIX}/src/lib/Plumed.cmake.static)
     elseif(PLUMED_MODE STREQUAL "SHARED")
-      include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.shared)
+      include(${PLUMED_LIBDIR}/plumed${PLUMED_SUFFIX}/src/lib/Plumed.cmake.shared)
     elseif(PLUMED_MODE STREQUAL "RUNTIME")
-      set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_COMPILE_DEFINITIONS "__PLUMED_DEFAULT_KERNEL=${PLUMED_LIBDIR}/${CMAKE_SHARED_LIBRARY_PREFIX}plumedKernel${CMAKE_SHARED_LIBRARY_SUFFIX}")
-      include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.runtime)
+      set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_COMPILE_DEFINITIONS "__PLUMED_DEFAULT_KERNEL=${PLUMED_LIBDIR}/${CMAKE_SHARED_LIBRARY_PREFIX}plumed${PLUMED_SUFFIX}Kernel${CMAKE_SHARED_LIBRARY_SUFFIX}")
+      include(${PLUMED_LIBDIR}/plumed${PLUMED_SUFFIX}/src/lib/Plumed.cmake.runtime)
     endif()
     set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_LINK_LIBRARIES "${PLUMED_LOAD}")
     set_target_properties(LAMMPS::PLUMED PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${PLUMED_INCLUDE_DIRS}")

cmake/packaging/build_linux_tgz.sh

@@ -59,12 +59,14 @@ done
 
 echo "Set up wrapper script"
 MYDIR=$(dirname "$0")
+cp ${MYDIR}/xdg-open ${DESTDIR}/bin
 cp ${MYDIR}/linux_wrapper.sh ${DESTDIR}/bin
 for s in ${DESTDIR}/bin/*
 do \
         EXE=$(basename $s)
         test ${EXE} = linux_wrapper.sh && continue
         test ${EXE} = qt.conf && continue
+        test ${EXE} = xdg-open && continue
         ln -s bin/linux_wrapper.sh ${DESTDIR}/${EXE}
 done
 

cmake/packaging/linux_wrapper.sh

@@ -4,15 +4,17 @@
 # reset locale to avoid problems with decimal numbers
 export LC_ALL=C
 
-BASEDIR=$(dirname "$0")
-EXENAME=$(basename "$0")
+BASEDIR="$(dirname "$0")"
+EXENAME="$(basename "$0")"
+
+PATH="${BASEDIR}/bin:${PATH}"
 
 # append to LD_LIBRARY_PATH to prefer local (newer) libs
-LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${BASEDIR}/lib
+LD_LIBRARY_PATH="${LD_LIBRARY_PATH}:${BASEDIR}/lib"
 
 # set some environment variables for LAMMPS etc.
-LAMMPS_POTENTIALS=${BASEDIR}/share/lammps/potentials
-MSI2LMP_LIBRARY=${BASEDIR}/share/lammps/frc_files
-export LD_LIBRARY_PATH LAMMPS_POTENTIALS MSI2LMP_LIBRARY
+LAMMPS_POTENTIALS="${BASEDIR}/share/lammps/potentials"
+MSI2LMP_LIBRARY="${BASEDIR}/share/lammps/frc_files"
+export LD_LIBRARY_PATH LAMMPS_POTENTIALS MSI2LMP_LIBRARY PATH
 
 exec "${BASEDIR}/bin/${EXENAME}" "$@"

cmake/presets/all_off.cmake

@@ -26,8 +26,9 @@ set(ALL_PACKAGES
   DPD-REACT
   DPD-SMOOTH
   DRUDE
-  ELECTRODE
   EFF
+  ELECTRODE
+  EXTRA-COMMAND
   EXTRA-COMPUTE
   EXTRA-DUMP
   EXTRA-FIX
@@ -60,6 +61,7 @@ set(ALL_PACKAGES
   ML-QUIP
   ML-RANN
   ML-SNAP
+  ML-UF3
   MOFFF
   MOLECULE
   MOLFILE

cmake/presets/all_on.cmake

@@ -28,8 +28,9 @@ set(ALL_PACKAGES
   DPD-REACT
   DPD-SMOOTH
   DRUDE
-  ELECTRODE
   EFF
+  ELECTRODE
+  EXTRA-COMMAND
   EXTRA-COMPUTE
   EXTRA-DUMP
   EXTRA-FIX
@@ -62,6 +63,7 @@ set(ALL_PACKAGES
   ML-QUIP
   ML-RANN
   ML-SNAP
+  ML-UF3
   MOFFF
   MOLECULE
   MOLFILE

cmake/presets/mingw-cross.cmake

@@ -22,8 +22,9 @@ set(WIN_PACKAGES
   DPD-REACT
   DPD-SMOOTH
   DRUDE
-  ELECTRODE
   EFF
+  ELECTRODE
+  EXTRA-COMMAND
   EXTRA-COMPUTE
   EXTRA-DUMP
   EXTRA-FIX
@@ -50,6 +51,7 @@ set(WIN_PACKAGES
   ML-POD
   ML-RANN
   ML-SNAP
+  ML-UF3
   MOFFF
   MOLECULE
   MOLFILE

cmake/presets/most.cmake

@@ -26,6 +26,7 @@ set(ALL_PACKAGES
   DRUDE
   EFF
   ELECTRODE
+  EXTRA-COMMAND
   EXTRA-COMPUTE
   EXTRA-DUMP
   EXTRA-FIX
@@ -45,6 +46,7 @@ set(ALL_PACKAGES
   ML-IAP
   ML-POD
   ML-SNAP
+  ML-UF3
   MOFFF
   MOLECULE
   OPENMP

cmake/presets/windows.cmake

@@ -22,6 +22,7 @@ set(WIN_PACKAGES
   DRUDE
   EFF
   ELECTRODE
+  EXTRA-COMMAND
   EXTRA-COMPUTE
   EXTRA-DUMP
   EXTRA-FIX
@@ -42,6 +43,7 @@ set(WIN_PACKAGES
   ML-IAP
   ML-POD
   ML-SNAP
+  ML-UF3
   MOFFF
   MOLECULE
   MOLFILE
@@ -50,8 +52,8 @@ set(WIN_PACKAGES
   ORIENT
   PERI
   PHONON
-  POEMS
   PLUGIN
+  POEMS
   PTM
   QEQ
   QTB

doc/lammps.1

@@ -1,7 +1,7 @@
-.TH LAMMPS "1" "17 April 2024" "2024-04-17"
+.TH LAMMPS "1" "27 June 2024" "2024-06-27"
 .SH NAME
 .B LAMMPS
-\- Molecular Dynamics Simulator.  Version 17 April 2024
+\- Molecular Dynamics Simulator.  Version 27 June 2024
 
 .SH SYNOPSIS
 .B lmp

doc/src/Commands_bond.rst

@@ -27,7 +27,7 @@ OPT.
 
    * :doc:`none <bond_none>`
    * :doc:`zero <bond_zero>`
-   * :doc:`hybrid <bond_hybrid>`
+   * :doc:`hybrid (k) <bond_hybrid>`
    *
    *
    *

doc/src/Commands_compute.rst

@@ -108,6 +108,10 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`pe/mol/tally <compute_tally>`
    * :doc:`pe/tally <compute_tally>`
    * :doc:`plasticity/atom <compute_plasticity_atom>`
+   * :doc:`pod/atom <compute_pod_atom>`
+   * :doc:`podd/atom <compute_pod_atom>`
+   * :doc:`pod/local <compute_pod_atom>`
+   * :doc:`pod/global <compute_pod_atom>`
    * :doc:`pressure <compute_pressure>`
    * :doc:`pressure/alchemy <compute_pressure_alchemy>`
    * :doc:`pressure/uef <compute_pressure_uef>`

doc/src/Commands_pair.rst

@@ -25,16 +25,16 @@ OPT.
 
    * :doc:`none <pair_none>`
    * :doc:`zero <pair_zero>`
-   * :doc:`hybrid (k) <pair_hybrid>`
-   * :doc:`hybrid/overlay (k) <pair_hybrid>`
-   * :doc:`hybrid/scaled <pair_hybrid>`
+   * :doc:`hybrid (ko) <pair_hybrid>`
+   * :doc:`hybrid/molecular (o) <pair_hybrid>`
+   * :doc:`hybrid/overlay (ko) <pair_hybrid>`
+   * :doc:`hybrid/scaled (o) <pair_hybrid>`
    * :doc:`kim <pair_kim>`
    * :doc:`list <pair_list>`
    * :doc:`tracker <pair_tracker>`
    *
    *
    *
-   *
    * :doc:`adp (ko) <pair_adp>`
    * :doc:`agni (o) <pair_agni>`
    * :doc:`aip/water/2dm (t) <pair_aip_water_2dm>`
@@ -94,9 +94,10 @@ OPT.
    * :doc:`coul/wolf (ko) <pair_coul>`
    * :doc:`coul/wolf/cs <pair_cs>`
    * :doc:`dpd (giko) <pair_dpd>`
-   * :doc:`dpd/fdt <pair_dpd_fdt>`
+   * :doc:`dpd/coul/slater/long (g) <pair_dpd_coul_slater_long>`
    * :doc:`dpd/ext (ko) <pair_dpd_ext>`
    * :doc:`dpd/ext/tstat (ko) <pair_dpd_ext>`
+   * :doc:`dpd/fdt <pair_dpd_fdt>`
    * :doc:`dpd/fdt/energy (k) <pair_dpd_fdt>`
    * :doc:`dpd/tstat (gko) <pair_dpd>`
    * :doc:`dsmc <pair_dsmc>`
@@ -246,7 +247,7 @@ OPT.
    * :doc:`pace (k) <pair_pace>`
    * :doc:`pace/extrapolation (k) <pair_pace>`
    * :doc:`pedone (o) <pair_pedone>`
-   * :doc:`pod <pair_pod>`
+   * :doc:`pod (k) <pair_pod>`
    * :doc:`peri/eps <pair_peri>`
    * :doc:`peri/lps (o) <pair_peri>`
    * :doc:`peri/pmb (o) <pair_peri>`
@@ -269,7 +270,7 @@ OPT.
    * :doc:`smd/ulsph <pair_smd_ulsph>`
    * :doc:`smtbq <pair_smtbq>`
    * :doc:`snap (ik) <pair_snap>`
-   * :doc:`soft (go) <pair_soft>`
+   * :doc:`soft (gko) <pair_soft>`
    * :doc:`sph/heatconduction (g) <pair_sph_heatconduction>`
    * :doc:`sph/idealgas <pair_sph_idealgas>`
    * :doc:`sph/lj (g) <pair_sph_lj>`
@@ -303,6 +304,7 @@ OPT.
    * :doc:`tip4p/long/soft (o) <pair_fep_soft>`
    * :doc:`tri/lj <pair_tri_lj>`
    * :doc:`ufm (got) <pair_ufm>`
+   * :doc:`uf3 (k) <pair_uf3>`
    * :doc:`vashishta (gko) <pair_vashishta>`
    * :doc:`vashishta/table (o) <pair_vashishta>`
    * :doc:`wf/cut <pair_wf_cut>`

doc/src/Commands_removed.rst

@@ -148,6 +148,14 @@ performance characteristics on NVIDIA GPUs. Both, the KOKKOS
 and the :ref:`GPU package <PKG-GPU>` are maintained
 and allow running LAMMPS with GPU acceleration.
 
+i-PI tool
+---------
+
+.. versionchanged:: 27June2024
+
+The i-PI tool has been removed from the LAMMPS distribution.  Instead,
+instructions to install i-PI from PyPI via pip are provided.
+
 restart2data tool
 -----------------
 

doc/src/Developer_utils.rst

@@ -211,6 +211,9 @@ Argument processing
 .. doxygenfunction:: bounds
    :project: progguide
 
+.. doxygenfunction:: bounds_typelabel
+   :project: progguide
+
 .. doxygenfunction:: expand_args
    :project: progguide
 

doc/src/Fortran.rst

@@ -305,6 +305,8 @@ of the contents of the :f:mod:`LIBLAMMPS` Fortran interface to LAMMPS.
    :ftype extract_setting: function
    :f extract_global: :f:func:`extract_global`
    :ftype extract_global: function
+   :f map_atom: :f:func:`map_atom`
+   :ftype map_atom: function
    :f extract_atom: :f:func:`extract_atom`
    :ftype extract_atom: function
    :f extract_compute: :f:func:`extract_compute`

doc/src/Howto_bioFF.rst

@@ -1,6 +1,10 @@
 CHARMM, AMBER, COMPASS, and DREIDING force fields
 =================================================
 
+A compact summary of the concepts, definitions, and properties of
+force fields with explicit bonded interactions (like the ones discussed
+in this HowTo) is given in :ref:`(Gissinger) <Typelabel2>`.
+
 A force field has 2 parts: the formulas that define it and the
 coefficients used for a particular system.  Here we only discuss
 formulas implemented in LAMMPS that correspond to formulas commonly used
@@ -11,12 +15,42 @@ commands like :doc:`pair_coeff <pair_coeff>` or :doc:`bond_coeff
 <bond_coeff>` and so on.  See the :doc:`Tools <Tools>` doc page for
 additional tools that can use CHARMM, AMBER, or Materials Studio
 generated files to assign force field coefficients and convert their
-output into LAMMPS input.
+output into LAMMPS input. LAMMPS input scripts can also be generated by
+`charmm-gui.org <https://charmm-gui.org/>`_.
 
-See :ref:`(MacKerell) <howto-MacKerell>` for a description of the CHARMM
-force field.  See :ref:`(Cornell) <howto-Cornell>` for a description of
-the AMBER force field.  See :ref:`(Sun) <howto-Sun>` for a description
-of the COMPASS force field.
+CHARMM and AMBER
+----------------
+
+The `CHARMM force field
+<https://mackerell.umaryland.edu/charmm_ff.shtml>`_ :ref:`(MacKerell)
+<howto-MacKerell>` and `AMBER force field
+<https://ambermd.org/AmberModels.php>`_ :ref:`(Cornell) <howto-Cornell>`
+have potential energy function of the form
+
+.. math::
+
+  V & = \sum_{bonds} E_b + \sum_{angles} \!E_a + \!\overbrace{\sum_{dihedral} \!\!E_d}^{\substack{
+         \text{charmm} \\
+        \text{charmmfsw}
+      }} +\!\!\! \sum_{impropers} \!\!\!E_i \\[.6em]
+      & \quad + \!\!\!\!\!\!\!\!\!\!\underbrace{~\sum_{pairs} \left(E_{LJ}+E_{coul}\right)}_{\substack{
+         \text{lj/charmm/coul/charmm} \\
+        \text{lj/charmm/coul/charmm/implicit} \\
+        \text{lj/charmm/coul/long} \\
+        \text{lj/charmm/coul/msm} \\
+         \text{lj/charmmfsw/coul/charmmfsh} \\
+        \text{lj/charmmfsw/coul/long}
+      }} \!\!\!\!\!\!\!\!+ \!\!\sum_{special}\! E_s + \!\!\!\!\sum_{residues} \!\!\!{\scriptstyle\mathrm{CMAP}(\phi,\psi)}
+
+
+The terms are computed by bond styles (relationship between 2 atoms),
+angle styles (between 3 atoms) , dihedral/improper styles (between 4
+atoms), pair styles (non-covalently bonded pair interactions) and
+special bonds. The CMAP term (see :doc:`fix cmap <fix_cmap>` command for
+details) corrects for pairs of dihedral angles ("Correction MAP") to
+significantly improve the structural and dynamic properties of proteins
+in crystalline and solution environments :ref:`(Brooks)
+<howto-Brooks>`. The AMBER force field does not include the CMAP term.
 
 The interaction styles listed below compute force field formulas that
 are consistent with common options in CHARMM or AMBER.  See each
@@ -31,10 +65,81 @@ command's documentation for the formula it computes.
 * :doc:`pair_style <pair_charmm>` lj/charmm/coul/charmm
 * :doc:`pair_style <pair_charmm>` lj/charmm/coul/charmm/implicit
 * :doc:`pair_style <pair_charmm>` lj/charmm/coul/long
-
 * :doc:`special_bonds <special_bonds>` charmm
 * :doc:`special_bonds <special_bonds>` amber
 
+The pair styles compute Lennard Jones (LJ) and Coulombic interactions
+with additional switching or shifting functions that ramp the energy
+and/or force smoothly to zero between an inner :math:`(a)` and outer
+:math:`(b)` cutoff. The older styles with *charmm* (not *charmmfsw* or
+*charmmfsh*\ ) in their name compute the LJ and Coulombic interactions
+with an energy switching function (esw) S(r) which ramps the energy
+smoothly to zero between the inner and outer cutoff. This can cause
+irregularities in pairwise forces (due to the discontinuous second
+derivative of energy at the boundaries of the switching region), which
+in some cases can result in complications in energy minimization and
+detectable artifacts in MD simulations.
+
+.. grid:: 1 1 2 2
+
+   .. grid-item::
+
+      .. math::
+
+         LJ(r) &= 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} -
+                  \left(\frac{\sigma}{r}\right)^6 \right]\\[.6em]
+         C(r) &= \frac{C q_i q_j}{ \epsilon r}\\[.6em]
+         S(r) &=  \frac{ \left(b^2 - r^2\right)^2 \left(b^2 + 2r^2 - 3{a^2}\right)}
+                 { \left(b^2 - a^2\right)^3 }\\[.6em]
+         E_{LJ}(r) &=  \begin{cases}
+           LJ(r), & r \leq a \\
+           LJ(r) S(r), & a < r \leq b \\
+           0, &r > b
+         \end{cases} \\[.6em]
+         E_{coul}(r) &=  \begin{cases}
+           C(r), & r \leq a \\
+           C(r) S(r), & a < r \leq b \\
+           0, & r > b
+         \end{cases}
+
+   .. grid-item::
+
+      .. image:: img/howto_charmm_ELJ.png
+         :align: center
+
+The newer styles with *charmmfsw* or *charmmfsh* in their name replace
+energy switching with force switching (fsw) for LJ interactions and
+force shifting (fsh) functions for Coulombic interactions
+:ref:`(Steinbach) <howto-Steinbach>`
+
+.. grid:: 1 1 2 2
+
+   .. grid-item::
+
+      .. math::
+
+           E_{LJ}(r) = & \begin{cases}
+       4  \epsilon \sigma^6  \left(\frac{\displaystyle\sigma
+         ^6-r^6}{\displaystyle r^{12}}-\frac{\displaystyle\sigma ^6}{\displaystyle a^6
+         b^6}+\frac{\displaystyle 1}{\displaystyle a^3 b^3}\right) & r\leq a \\
+       \frac{\displaystyle 4 \epsilon \sigma^6   \left(\sigma ^6
+         \left(b^6-r^6\right)^2-b^3 r^6 \left(a^3+b^3\right)
+         \left(b^3-r^3\right)^2\right)}{\displaystyle b^6 r^{12}
+         \left(b^6-a^6\right)} & a<r \leq b\\
+         0, & r>b
+        \end{cases}\\[.6em]
+         E_{coul}(r) & =  \begin{cases}
+              C(r) \frac{\displaystyle (b-r)^2}{\displaystyle r b^2}, &  r \leq b \\
+              0, & r > b
+            \end{cases}
+
+   .. grid-item::
+      .. image:: img/howto_charmmfsw_ELJ.png
+         :align: center
+
+These styles are used by LAMMPS input scripts generated by
+https://charmm-gui.org/ :ref:`(Brooks) <howto-Brooks>`.
+
 .. note::
 
    For CHARMM, newer *charmmfsw* or *charmmfsh* styles were released in
@@ -43,17 +148,33 @@ command's documentation for the formula it computes.
    <pair_charmm>` and :doc:`dihedral charmm <dihedral_charmm>` doc
    pages.
 
+.. note::
+
+  The TIP3P water model is strongly recommended for use with the CHARMM
+  force field. In fact, `"using the SPC model with CHARMM parameters is
+  a bad idea"
+  <https://matsci.org/t/using-spc-water-with-charmm-ff/24715>`_ and `"to
+  enable TIP4P style water in CHARMM, you would have to write a new pair
+  style"
+  <https://matsci.org/t/hybrid-pair-styles-for-charmm-and-tip4p-ew/32609>`_
+  . LAMMPS input scripts generated by Solution Builder on https://charmm-gui.org
+  use TIP3P molecules for solvation.  Any other water model can and
+  probably will lead to false conclusions.
+
+COMPASS
+-------
+
 COMPASS is a general force field for atomistic simulation of common
 organic molecules, inorganic small molecules, and polymers which was
-developed using ab initio and empirical parameterization techniques.
-See the :doc:`Tools <Tools>` page for the msi2lmp tool for creating
-LAMMPS template input and data files from BIOVIA's Materials Studio
-files.  Please note that the msi2lmp tool is very old and largely
-unmaintained, so it does not support all features of Materials Studio
-provided force field files, especially additions during the last decade.
-You should watch the output carefully and compare results, where
-possible.  See :ref:`(Sun) <howto-Sun>` for a description of the COMPASS force
-field.
+developed using ab initio and empirical parameterization techniques
+:ref:`(Sun) <howto-Sun>`.  See the :doc:`Tools <Tools>` page for the
+msi2lmp tool for creating LAMMPS template input and data files from
+BIOVIA's Materials Studio files.  Please note that the msi2lmp tool is
+very old and largely unmaintained, so it does not support all features
+of Materials Studio provided force field files, especially additions
+during the last decade.  You should watch the output carefully and
+compare results, where possible.  See :ref:`(Sun) <howto-Sun>` for a
+description of the COMPASS force field.
 
 These interaction styles listed below compute force field formulas that
 are consistent with the COMPASS force field.  See each command's
@@ -70,14 +191,21 @@ documentation for the formula it computes.
 
 * :doc:`special_bonds <special_bonds>` lj/coul 0 0 1
 
-DREIDING is a generic force field developed by the `Goddard group <http://www.wag.caltech.edu>`_ at Caltech and is useful for
-predicting structures and dynamics of organic, biological and main-group
-inorganic molecules. The philosophy in DREIDING is to use general force
-constants and geometry parameters based on simple hybridization
-considerations, rather than individual force constants and geometric
-parameters that depend on the particular combinations of atoms involved
-in the bond, angle, or torsion terms. DREIDING has an :doc:`explicit hydrogen bond term <pair_hbond_dreiding>` to describe interactions involving a
-hydrogen atom on very electronegative atoms (N, O, F).
+DREIDING
+--------
+
+DREIDING is a generic force field developed by the `Goddard group
+<http://www.wag.caltech.edu>`_ at Caltech and is useful for predicting
+structures and dynamics of organic, biological and main-group inorganic
+molecules.  The philosophy in DREIDING is to use general force constants
+and geometry parameters based on simple hybridization considerations,
+rather than individual force constants and geometric parameters that
+depend on the particular combinations of atoms involved in the bond,
+angle, or torsion terms.  DREIDING has an :doc:`explicit hydrogen bond
+term <pair_hbond_dreiding>` to describe interactions involving a
+hydrogen atom on very electronegative atoms (N, O, F).  Unlike CHARMM
+or AMBER, the DREIDING force field has not been parameterized for
+considering solvents (like water).
 
 See :ref:`(Mayo) <howto-Mayo>` for a description of the DREIDING force field
 
@@ -110,21 +238,31 @@ documentation for the formula it computes.
 
 ----------
 
+.. _Typelabel2:
+
+**(Gissinger)** J. R. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M. Choi, D. S. Karls, A. Stukowski, W. Im, H. Heinz, A. Kohlmeyer, and E. B. Tadmor, J Phys Chem B, 128, 3282-3297 (2024).
+
 .. _howto-MacKerell:
 
-**(MacKerell)** MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field,
-Fischer, Gao, Guo, Ha, et al, J Phys Chem, 102, 3586 (1998).
+**(MacKerell)** MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field, Fischer, Gao, Guo, Ha, et al (1998).  J Phys Chem, 102, 3586 . https://doi.org/10.1021/jp973084f
 
 .. _howto-Cornell:
 
-**(Cornell)** Cornell, Cieplak, Bayly, Gould, Merz, Ferguson,
-Spellmeyer, Fox, Caldwell, Kollman, JACS 117, 5179-5197 (1995).
+**(Cornell)** Cornell, Cieplak, Bayly, Gould, Merz, Ferguson, Spellmeyer, Fox, Caldwell, Kollman (1995).  JACS 117, 5179-5197. https://doi.org/10.1021/ja00124a002
+
+.. _howto-Steinbach:
+
+**(Steinbach)** Steinbach, Brooks (1994). J Comput Chem, 15, 667. https://doi.org/10.1002/jcc.540150702
+
+.. _howto-Brooks:
+
+**(Brooks)** Brooks, et al (2009). J Comput Chem, 30, 1545. https://onlinelibrary.wiley.com/doi/10.1002/jcc.21287
 
 .. _howto-Sun:
 
-**(Sun)** Sun, J. Phys. Chem. B, 102, 7338-7364 (1998).
+**(Sun)** Sun (1998). J. Phys. Chem. B, 102, 7338-7364. https://doi.org/10.1021/jp980939v
 
 .. _howto-Mayo:
 
-**(Mayo)** Mayo, Olfason, Goddard III, J Phys Chem, 94, 8897-8909
-(1990).
+**(Mayo)** Mayo, Olfason, Goddard III (1990). J Phys Chem, 94, 8897-8909. https://doi.org/10.1021/j100389a010
+

doc/src/Howto_bpm.rst

@@ -15,7 +15,8 @@ orientation for rotational models. This produces a stress-free initial
 state. Furthermore, bonds are allowed to break under large strains,
 producing fracture. The examples/bpm directory has sample input scripts
 for simulations of the fragmentation of an impacted plate and the
-pouring of extended, elastic bodies.
+pouring of extended, elastic bodies. See :ref:`(Clemmer) <howto-Clemmer>`
+for more general information on the approach and the LAMMPS implementation.
 
 ----------
 
@@ -150,3 +151,9 @@ the following are currently compatible with BPM bond styles:
    interactions, one will need to switch between different *special_bonds*
    settings in the input script. An example is found in
    ``examples/bpm/impact``.
+
+----------
+
+.. _howto-Clemmer:
+
+**(Clemmer)** Clemmer, Monti, Lechman, Soft Matter, 20, 1702 (2024).

doc/src/Library_properties.rst

@@ -13,6 +13,7 @@ This section documents the following functions:
 - :cpp:func:`lammps_extract_setting`
 - :cpp:func:`lammps_extract_global_datatype`
 - :cpp:func:`lammps_extract_global`
+- :cpp:func:`lammps_map_atom`
 
 --------------------
 
@@ -120,3 +121,8 @@ subdomains and processors.
 .. doxygenfunction:: lammps_extract_global
    :project: progguide
 
+-----------------------
+
+.. doxygenfunction:: lammps_map_atom
+   :project: progguide
+

doc/src/Packages_details.rst

@@ -52,6 +52,7 @@ page gives those details.
    * :ref:`DRUDE <PKG-DRUDE>`
    * :ref:`EFF <PKG-EFF>`
    * :ref:`ELECTRODE <PKG-ELECTRODE>`
+   * :ref:`EXTRA-COMMAND <PKG-EXTRA-COMMAND>`
    * :ref:`EXTRA-COMPUTE <PKG-EXTRA-COMPUTE>`
    * :ref:`EXTRA-DUMP <PKG-EXTRA-DUMP>`
    * :ref:`EXTRA-FIX <PKG-EXTRA-FIX>`
@@ -84,6 +85,7 @@ page gives those details.
    * :ref:`ML-QUIP <PKG-ML-QUIP>`
    * :ref:`ML-RANN <PKG-ML-RANN>`
    * :ref:`ML-SNAP <PKG-ML-SNAP>`
+   * :ref:`ML-UF3 <PKG-ML-UF3>`
    * :ref:`MOFFF <PKG-MOFFF>`
    * :ref:`MOLECULE <PKG-MOLECULE>`
    * :ref:`MOLFILE <PKG-MOLFILE>`
@@ -403,6 +405,7 @@ and :ref:`ASPHERE <PKG-ASPHERE>` packages are installed.
 * :doc:`bond_style oxdna2/\* <bond_oxdna>`
 * :doc:`bond_style oxrna2/\* <bond_oxdna>`
 * :doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`
+* examples/PACKAGES/cgdna
 
 ----------
 
@@ -676,7 +679,12 @@ DPD-BASIC package
 Pair styles for the basic dissipative particle dynamics (DPD) method
 and DPD thermostatting.
 
-**Author:** Kurt Smith (U Pittsburgh), Martin Svoboda, Martin Lisal (ICPF and UJEP)
+Pair style :doc:`dpd/coul/slater/long <pair_dpd_coul_slater_long>` also
+includes smeared charges for coulomb interactions and thus requires the
+:ref:`KSPACE <PKG-KSPACE>` package to be installed to handle the long-range
+Coulomb part of the interactions.
+
+**Authors:** Kurt Smith (U Pittsburgh), Martin Svoboda, Martin Lisal (ICPF and UJEP), Eddy Barraud (IFPEN)
 
 **Supporting info:**
 
@@ -685,6 +693,7 @@ and DPD thermostatting.
 * :doc:`pair_style dpd/tstat <pair_dpd>`
 * :doc:`pair_style dpd/ext <pair_dpd_ext>`
 * :doc:`pair_style dpd/ext/tstat <pair_dpd_ext>`
+* :doc:`pair_style dpd/coul/slater/long <pair_dpd_coul_slater_long>`
 * examples/PACKAGES/dpd-basic
 
 ----------
@@ -886,6 +895,22 @@ This package has :ref:`specific installation instructions <electrode>` on the
 
 ----------
 
+.. _PKG-EXTRA-COMMAND:
+
+EXTRA-COMMAND package
+---------------------
+
+**Contents:**
+
+Additional command styles that are less commonly used.
+
+**Supporting info:**
+
+* src/EXTRA-COMMAND: filenames -> commands
+* :doc:`general commands <Commands_all>`
+
+----------
+
 .. _PKG-EXTRA-COMPUTE:
 
 EXTRA-COMPUTE package
@@ -1925,6 +1950,31 @@ computes which analyze attributes of the potential.
 
 ----------
 
+.. _PKG-ML-UF3:
+
+ML-UF3 package
+--------------
+
+**Contents:**
+
+A pair style for the ultra-fast force field potentials (UF3). UF3 is a
+methodology for deriving a highly accurate classical potential which is
+fast to evaluate and is fitted to a large archives of quantum mechanical
+(DFT) data.  The use of b-spline basis set in UF3 enables the rapid
+evaluation of 2-body and 3-body interactions.
+
+**Authors:** Ajinkya C Hire (University of Florida),
+Hendrik Krass (University of Constance),
+Matthias Rupp (Luxembourg Institute of Science and Technology),
+Richard Hennig (University of Florida)
+
+**Supporting info:**
+
+* src/ML-UF3: filenames -> commands
+* :doc:`pair_style uf3 <pair_uf3>`
+* examples/uf3
+* https://github.com/uf3/uf3
+
 .. _PKG-MOFFF:
 
 MOFFF package

doc/src/Packages_list.rst

@@ -158,6 +158,11 @@ whether an extra library is needed to build and use the package:
      - :doc:`fix electrode/conp <fix_electrode>`
      - PACKAGES/electrode
      - no
+   * - :ref:`EXTRA-COMMAND <PKG-EXTRA-COMMAND>`
+     - additional command styles
+     - :doc:`general commands <Commands_all>`
+     - n/a
+     - no
    * - :ref:`EXTRA-COMPUTE <PKG-EXTRA-COMPUTE>`
      - additional compute styles
      - :doc:`compute <compute>`
@@ -318,6 +323,11 @@ whether an extra library is needed to build and use the package:
      - :doc:`pair_style snap <pair_snap>`
      - snap
      - no
+   * - :ref:`ML-UF3 <PKG-ML-UF3>`
+     - quantum-fitted ultra fast potentials
+     - :doc:`pair_style uf3 <pair_uf3>`
+     - PACKAGES/uf3
+     - no
    * - :ref:`MOFFF <PKG-MOFFF>`
      - styles for `MOF-FF <MOFplus_>`_ force field
      - :doc:`pair_style buck6d/coul/gauss <pair_buck6d_coul_gauss>`

doc/src/Tools.rst

@@ -90,7 +90,7 @@ Miscellaneous tools
 
    * :ref:`LAMMPS coding standards <coding_standard>`
    * :ref:`emacs <emacs>`
-   * :ref:`i-pi <ipi>`
+   * :ref:`i-PI <ipi>`
    * :ref:`kate <kate>`
    * :ref:`LAMMPS shell <lammps_shell>`
    * :ref:`LAMMPS GUI <lammps_gui>`
@@ -376,21 +376,40 @@ See README file in the tools/fep directory.
 
 .. _ipi:
 
-i-pi tool
+i-PI tool
 -------------------
 
-The tools/i-pi directory contains a version of the i-PI package, with
-all the LAMMPS-unrelated files removed.  It is provided so that it can
-be used with the :doc:`fix ipi <fix_ipi>` command to perform
-path-integral molecular dynamics (PIMD).
+.. versionchanged:: 27June2024
+
+The tools/i-pi directory used to contain a bundled version of the i-PI
+software package for use with LAMMPS.  This version, however, was
+removed in 06/2024.
 
 The i-PI package was created and is maintained by Michele Ceriotti,
 michele.ceriotti at gmail.com, to interface to a variety of molecular
 dynamics codes.
 
-See the tools/i-pi/manual.pdf file for an overview of i-PI, and the
-:doc:`fix ipi <fix_ipi>` page for further details on running PIMD
-calculations with LAMMPS.
+i-PI is now available via PyPI using the pip package manager at:
+https://pypi.org/project/ipi/
+
+Here are the commands to set up a virtual environment and install
+i-PI into it with all its dependencies.
+
+.. code-block:: sh
+
+   python -m venv ipienv
+   source ipienv/bin/activate
+   pip install --upgrade pip
+   pip install ipi
+
+To install the development version from GitHub, please use:
+
+.. code-block:: sh
+
+   pip install git+https://github.com/i-pi/i-pi.git
+
+For further information, please consult the [i-PI home
+page](https://ipi-code.org).
 
 ----------
 
@@ -709,8 +728,8 @@ CMake is required.
 The LAMMPS GUI has been successfully compiled and tested on:
 
 - Ubuntu Linux 20.04LTS x86_64 using GCC 9, Qt version 5.12
-- Fedora Linux 38 x86\_64 using GCC 13 and Clang 16, Qt version 5.15LTS
-- Fedora Linux 38 x86\_64 using GCC 13, Qt version 6.5LTS
+- Fedora Linux 40 x86\_64 using GCC 14 and Clang 17, Qt version 5.15LTS
+- Fedora Linux 40 x86\_64 using GCC 14, Qt version 6.5LTS
 - Apple macOS 12 (Monterey) and macOS 13 (Ventura) with Xcode on arm64 and x86\_64, Qt version 5.15LTS
 - Windows 10 and 11 x86_64 with Visual Studio 2022 and Visual C++ 14.36, Qt version 5.15LTS
 - Windows 10 and 11 x86_64 with MinGW / GCC 10.0 cross-compiler on Fedora 38, Qt version 5.15LTS
@@ -752,22 +771,23 @@ if necessary.  When both, Qt5 and Qt6 are available, Qt6 will be preferred
 unless ``-D LAMMPS_GUI_USE_QT5=yes`` is set.
 
 It should be possible to build the LAMMPS GUI as a standalone
-compilation (e.g. when LAMMPS has been compiled with traditional make),
-then the CMake configuration needs to be told where to find the LAMMPS
+compilation (e.g. when LAMMPS has been compiled with traditional make).
+Then the CMake configuration needs to be told where to find the LAMMPS
 headers and the LAMMPS library, via ``-D
 LAMMPS_SOURCE_DIR=/path/to/lammps/src``.  CMake will try to guess a
 build folder with the LAMMPS library from that path, but it can also be
 set with ``-D LAMMPS_LIB_DIR=/path/to/lammps/lib``.
 
 Rather than linking to the LAMMPS library during compilation, it is also
-possible to compile the GUI with a plugin loader library that will load
+possible to compile the GUI with a plugin loader that will load
 the LAMMPS library dynamically at runtime during the start of the GUI
 from a shared library; e.g. ``liblammps.so`` or ``liblammps.dylib`` or
 ``liblammps.dll`` (depending on the operating system).  This has the
-advantage that the LAMMPS library can be updated LAMMPS without having
-to recompile the GUI.  The ABI of the LAMMPS C-library interface is very
-stable and generally backward compatible.  This feature is enabled by
-setting ``-D LAMMPS_GUI_USE_PLUGIN=on`` and then ``-D
+advantage that the LAMMPS library can be built from updated or modified
+LAMMPS source without having to recompile the GUI.  The ABI of the
+LAMMPS C-library interface is very stable and generally backward
+compatible.  This feature is enabled by setting
+``-D LAMMPS_GUI_USE_PLUGIN=on`` and then ``-D
 LAMMPS_PLUGINLIB_DIR=/path/to/lammps/plugin/loader``. Typically, this
 would be the ``examples/COUPLE/plugin`` folder of the LAMMPS
 distribution.
@@ -779,8 +799,8 @@ macOS
 """""
 
 When building on macOS, the build procedure will try to manufacture a
-drag-n-drop installer, LAMMPS-macOS-multiarch.dmg, when using the 'dmg'
-target (i.e. ``cmake --build <build dir> --target dmg`` or ``make dmg``.
+drag-n-drop installer, ``LAMMPS-macOS-multiarch.dmg``, when using the
+'dmg' target (i.e. ``cmake --build <build dir> --target dmg`` or ``make dmg``.
 
 To build multi-arch executables that will run on both, arm64 and x86_64
 architectures natively, it is necessary to set the CMake variable ``-D
@@ -819,12 +839,12 @@ and LAMMPS GUI can be launched from anywhere from the command line.
 
 The standard CMake build procedure can be applied and the
 ``mingw-cross.cmake`` preset used. By using ``mingw64-cmake`` the CMake
-command will automatically include a suitable CMake toolset file (the
-regular cmake command can be used after that).  After building the
-libraries and executables, you can build the target 'zip'
-(i.e. ``cmake --build <build dir> --target zip`` or ``make zip``
-to stage all installed files into a LAMMPS_GUI folder and then
-run a script to copy all required dependencies, some other files,
+command will automatically include a suitable CMake toolchain file (the
+regular cmake command can be used after that to modify the configuration
+settings, if needed).  After building the libraries and executables,
+you can build the target 'zip' (i.e. ``cmake --build <build dir> --target zip``
+or ``make zip`` to stage all installed files into a LAMMPS_GUI folder
+and then run a script to copy all required dependencies, some other files,
 and create a zip file from it.
 
 Linux

doc/src/bond_hybrid.rst

@@ -1,8 +1,11 @@
 .. index:: bond_style hybrid
+.. index:: bond_style hybrid/kk
 
 bond_style hybrid command
 =========================
 
+Accelerator Variants: *hybrid/kk*
+
 Syntax
 """"""
 
@@ -15,7 +18,7 @@ Syntax
 Examples
 """"""""
 
-.. code-block: LAMMPS
+.. code-block:: LAMMPS
 
    bond_style hybrid harmonic fene
    bond_coeff 1 harmonic 80.0 1.2
@@ -60,6 +63,10 @@ bond types.
 
 ----------
 
+.. include:: accel_styles.rst
+
+----------
+
 Restrictions
 """"""""""""
 

doc/src/bond_oxdna.rst

@@ -27,6 +27,7 @@ Examples
 
 .. code-block:: LAMMPS
 
+   # LJ units
    bond_style oxdna/fene
    bond_coeff * 2.0 0.25 0.7525
 
@@ -36,6 +37,32 @@ Examples
    bond_style oxrna2/fene
    bond_coeff * 2.0 0.25 0.76107
 
+   bond_style oxdna/fene
+   bond_coeff * oxdna_lj.cgdna
+
+   # Real units
+   bond_style oxdna/fene
+   bond_coeff * 11.92337812042065 2.1295 6.409795
+
+   bond_style oxdna2/fene
+   bond_coeff * 11.92337812042065 2.1295 6.4430152
+
+   bond_style oxrna2/fene
+   bond_coeff * 11.92337812042065 2.1295 6.482800913
+
+   bond_style oxrna2/fene
+   bond_coeff * oxrna2_real.cgdna
+
+.. note::
+
+   The coefficients in the above examples have to be kept fixed and
+   cannot be changed without reparameterizing the entire model. They are
+   provided in forms compatible with both *units lj* and *units real*
+   (see documentation of :doc:`units <units>`).  These can also be read
+   from a potential file with correct unit style by specifying the name
+   of the file. Several potential files for each unit style are included
+   in the ``potentials`` directory of the LAMMPS distribution.
+
 Description
 """""""""""
 
@@ -46,15 +73,14 @@ The *oxdna/fene*, *oxdna2/fene*, and *oxrna2/fene* bond styles use the potential
    E = - \frac{\epsilon}{2} \ln \left[ 1 - \left(\frac{r-r_0}{\Delta}\right)^2\right]
 
 to define a modified finite extensible nonlinear elastic (FENE)
-potential :ref:`(Ouldridge) <Ouldridge0>` to model the connectivity of the
-phosphate backbone in the oxDNA/oxRNA force field for coarse-grained
+potential :ref:`(Ouldridge) <Ouldridge0>` to model the connectivity of
+the phosphate backbone in the oxDNA/oxRNA force field for coarse-grained
 modelling of DNA/RNA.
 
 The following coefficients must be defined for the bond type via the
 :doc:`bond_coeff <bond_coeff>` command as given in the above example, or
-in the data file or restart files read by the
-:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
-commands:
+in the data file or restart files read by the :doc:`read_data
+<read_data>` or :doc:`read_restart <read_restart>` commands:
 
 * :math:`\epsilon` (energy)
 * :math:`\Delta` (distance)
@@ -62,41 +88,40 @@ commands:
 
 .. note::
 
-   The oxDNA bond style has to be used together with the
-   corresponding oxDNA pair styles for excluded volume interaction
-   *oxdna/excv* , stacking *oxdna/stk* , cross-stacking *oxdna/xstk* and
-   coaxial stacking interaction *oxdna/coaxstk* as well as
-   hydrogen-bonding interaction *oxdna/hbond* (see also documentation of
-   :doc:`pair_style oxdna/excv <pair_oxdna>`). For the oxDNA2
-   :ref:`(Snodin) <Snodin0>` bond style the analogous pair styles
-   *oxdna2/excv* , *oxdna2/stk* , *oxdna2/xstk* , *oxdna2/coaxstk* ,
-   *oxdna2/hbond* and an additional Debye-Hueckel pair style
-   *oxdna2/dh* have to be defined. The same applies to the oxRNA2
-   :ref:`(Sulc1) <Sulc01>` styles.
-   The coefficients in the above example have to be kept fixed and cannot
-   be changed without reparameterizing the entire model.
+   The oxDNA bond style has to be used together with the corresponding
+   oxDNA pair styles for excluded volume interaction *oxdna/excv* ,
+   stacking *oxdna/stk* , cross-stacking *oxdna/xstk* and coaxial
+   stacking interaction *oxdna/coaxstk* as well as hydrogen-bonding
+   interaction *oxdna/hbond* (see also documentation of :doc:`pair_style
+   oxdna/excv <pair_oxdna>`). For the oxDNA2 :ref:`(Snodin) <Snodin0>`
+   bond style the analogous pair styles *oxdna2/excv* , *oxdna2/stk* ,
+   *oxdna2/xstk* , *oxdna2/coaxstk* , *oxdna2/hbond* and an additional
+   Debye-Hueckel pair style *oxdna2/dh* have to be defined. The same
+   applies to the oxRNA2 :ref:`(Sulc1) <Sulc01>` styles.
 
 .. note::
 
-   This bond style has to be used with the *atom_style hybrid bond ellipsoid oxdna*
-   (see documentation of :doc:`atom_style <atom_style>`). The *atom_style oxdna*
-   stores the 3'-to-5' polarity of the nucleotide strand, which is set through
-   the bond topology in the data file. The first (second) atom in a bond definition
-   is understood to point towards the 3'-end (5'-end) of the strand.
+   This bond style has to be used with the *atom_style hybrid bond
+   ellipsoid oxdna* (see documentation of :doc:`atom_style
+   <atom_style>`). The *atom_style oxdna* stores the 3'-to-5' polarity
+   of the nucleotide strand, which is set through the bond topology in
+   the data file. The first (second) atom in a bond definition is
+   understood to point towards the 3'-end (5'-end) of the strand.
 
 .. warning::
 
-   If data files are produced with :doc:`write_data <write_data>`, then the
-   :doc:`newton <newton>` command should be set to *newton on* or *newton off on*.
-   Otherwise the data files will not have the same 3'-to-5' polarity as the
-   initial data file. This limitation does not apply to binary restart files
-   produced with :doc:`write_restart <write_restart>`.
+   If data files are produced with :doc:`write_data <write_data>`, then
+   the :doc:`newton <newton>` command should be set to *newton on* or
+   *newton off on*.  Otherwise the data files will not have the same
+   3'-to-5' polarity as the initial data file. This limitation does not
+   apply to binary restart files produced with :doc:`write_restart
+   <write_restart>`.
 
 Example input and data files for DNA and RNA duplexes can be found in
-examples/PACKAGES/cgdna/examples/oxDNA/ , /oxDNA2/ and /oxRNA2/.  A simple python
-setup tool which creates single straight or helical DNA strands, DNA/RNA
-duplexes or arrays of DNA/RNA duplexes can be found in
-examples/PACKAGES/cgdna/util/.
+``examples/PACKAGES/cgdna/examples/oxDNA/`, `.../oxDNA2/`` and
+``.../oxRNA2/``.  A simple python setup tool which creates single
+straight or helical DNA strands, DNA/RNA duplexes or arrays of DNA/RNA
+duplexes can be found in ``examples/PACKAGES/cgdna/util/``.
 
 Please cite :ref:`(Henrich) <Henrich0>` in any publication that uses
 this implementation. An updated documentation that contains general information
@@ -113,6 +138,39 @@ and for sequence-specific hydrogen-bonding and stacking interactions
 
 ----------
 
+Potential file reading
+""""""""""""""""""""""
+
+For each style oxdna, oxdna2 and oxrna2, the first parameter argument
+can be a filename, and if it is, no further arguments should be
+supplied. Therefore the following command:
+
+.. code-block:: LAMMPS
+
+   bond_style oxdna/fene
+   bond_coeff * oxdna_lj.cgdna
+
+will be interpreted as a request to read the (FENE) potential
+:ref:`(Ouldridge) <Ouldridge0>` parameters from the file with the given
+name.  The file can define multiple potential parameters for both bonded
+and pair interactions, but for the above bonded interactions there must
+exist in the file a line of the form:
+
+.. code-block:: LAMMPS
+
+   *   fene    epsilon delta r0
+
+There are sample potential files for each unit style in the
+``potentials`` directory of the LAMMPS distribution. The potential file
+unit system must align with the units defined via the :doc:`units
+<units>` command. For conversion between different *LJ* and *real* unit
+systems for oxDNA, the python tool *lj2real.py* located in the
+``examples/PACKAGES/cgdna/util/`` directory can be used. This tool
+assumes similar file structure to the examples found in
+``examples/PACKAGES/cgdna/examples/``.
+
+----------
+
 Restrictions
 """"""""""""
 

doc/src/compute.rst

@@ -272,6 +272,10 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
 * :doc:`pe/mol/tally <compute_tally>` - potential energy between two groups of atoms separated into intermolecular and intramolecular components via the tally callback mechanism
 * :doc:`pe/tally <compute_tally>` - potential energy between two groups of atoms via the tally callback mechanism
 * :doc:`plasticity/atom <compute_plasticity_atom>` - Peridynamic plasticity for each atom
+* :doc:`pod/atom <compute_pod_atom>` - POD descriptors for each atom
+* :doc:`podd/atom <compute_pod_atom>` - derivative of POD descriptors for each atom
+* :doc:`pod/local <compute_pod_atom>` - local POD descriptors and their derivatives
+* :doc:`pod/global <compute_pod_atom>` - global POD descriptors and their derivatives
 * :doc:`pressure <compute_pressure>` - total pressure and pressure tensor
 * :doc:`pressure/alchemy <compute_pressure_alchemy>` - mixed system total pressure and pressure tensor for :doc:`fix alchemy <fix_alchemy>` runs
 * :doc:`pressure/uef <compute_pressure_uef>` - pressure tensor in the reference frame of an applied flow field

doc/src/compute_count_type.rst

@@ -12,7 +12,7 @@ Syntax
 
 * ID, group-ID are documented in :doc:`compute <compute>` command
 * count/type = style name of this compute command
-* mode = {atom} or {bond} or {angle} or {dihedral} or {improper}
+* mode = *atom* or *bond* or *angle* or *dihedral* or *improper*
 
 Examples
 """"""""
@@ -69,29 +69,42 @@ for each type:
 
 ----------
 
-If the {mode} setting is {atom} then the count of atoms for each atom
+If the *mode* setting is *atom* then the count of atoms for each atom
 type is tallied.  Only atoms in the specified group are counted.
 
-If the {mode} setting is {bond} then the count of bonds for each bond
+The atom count for each type can be normalized by the total number of
+atoms like so:
+
+.. code-block:: LAMMPS
+
+   compute typevec all count/type atom # number of atoms of each type
+   variable normtypes vector c_typevec/atoms # divide by total number of atoms
+   variable ntypes equal extract_setting(ntypes) # number of atom types
+   thermo_style custom step v_normtypes[*${ntypes}] # vector variable needs upper limit
+
+Similarly, bond counts can be normalized by the total number of bonds.
+The same goes for angles, dihedrals, and impropers (see below).
+
+If the *mode* setting is *bond* then the count of bonds for each bond
 type is tallied.  Only bonds with both atoms in the specified group
 are counted.
 
-For {mode} = {bond}, broken bonds with a bond type of zero are also
+For *mode* = *bond*, broken bonds with a bond type of zero are also
 counted.  The :doc:`bond_style quartic <bond_quartic>` and :doc:`BPM
-bond styles <Howto_bpm>` break bonds by doing this.  See the :doc:`
-Howto broken bonds <Howto_broken_bonds>` doc page for more details.
+bond styles <Howto_bpm>` break bonds by doing this.  See the
+:doc:`Howto broken bonds <Howto_broken_bonds>` doc page for more details.
 Note that the group setting is ignored for broken bonds; all broken
 bonds in the system are counted.
 
-If the {mode} setting is {angle} then the count of angles for each
+If the *mode* setting is *angle* then the count of angles for each
 angle type is tallied.  Only angles with all 3 atoms in the specified
 group are counted.
 
-If the {mode} setting is {dihedral} then the count of dihedrals for
+If the *mode* setting is *dihedral* then the count of dihedrals for
 each dihedral type is tallied.  Only dihedrals with all 4 atoms in the
 specified group are counted.
 
-If the {mode} setting is {improper} then the count of impropers for
+If the *mode* setting is *improper* then the count of impropers for
 each improper type is tallied.  Only impropers with all 4 atoms in the
 specified group are counted.
 
@@ -101,18 +114,19 @@ Output info
 """""""""""
 
 This compute calculates a global vector of counts.  If the mode is
-{atom} or {bond} or {angle} or {dihedral} or {improper}, then the
+*atom* or *bond* or *angle* or *dihedral* or *improper*, then the
 vector length is the number of atom types or bond types or angle types
 or dihedral types or improper types, respectively.
 
-If the mode is {bond} this compute also calculates a global scalar
+If the mode is *bond* this compute also calculates a global scalar
 which is the number of broken bonds with type = 0, as explained above.
 
 These values can be used by any command that uses global scalar or
 vector values from a compute as input.  See the :doc:`Howto output
 <Howto_output>` page for an overview of LAMMPS output options.
 
-The scalar and vector values calculated by this compute are "extensive".
+The scalar and vector values calculated by this compute are both
+"intensive".
 
 Restrictions
 """"""""""""

doc/src/compute_pod_atom.rst

@@ -0,0 +1,145 @@
+.. index:: compute pod/atom
+.. index:: compute podd/atom
+.. index:: compute pod/local
+.. index:: compute pod/global
+
+compute pod/atom command
+========================
+
+compute podd/atom command
+=========================
+
+compute pod/local command
+=========================
+
+compute pod/global command
+==========================
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   compute ID group-ID pod/atom param.pod coefficients.pod
+   compute ID group-ID podd/atom param.pod coefficients.pod
+   compute ID group-ID pod/local param.pod coefficients.pod
+   compute ID group-ID pod/global param.pod coefficients.pod
+
+* ID, group-ID are documented in :doc:`compute <compute>` command
+* pod/atom = style name of this compute command
+* param.pod = the parameter file specifies parameters of the POD descriptors
+* coefficients.pod = the coefficient file specifies coefficients of the POD potential
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   compute d all pod/atom Ta_param.pod
+   compute dd all podd/atom Ta_param.pod
+   compute ldd all pod/local Ta_param.pod
+   compute gdd all podd/global Ta_param.pod
+   compute d all pod/atom Ta_param.pod Ta_coefficients.pod
+   compute dd all podd/atom Ta_param.pod Ta_coefficients.pod
+   compute ldd all pod/local Ta_param.pod Ta_coefficients.pod
+   compute gdd all podd/global Ta_param.pod Ta_coefficients.pod
+
+Description
+"""""""""""
+
+.. versionadded:: 27June2024
+
+Define a computation that calculates a set of quantities related to the
+POD descriptors of the atoms in a group. These computes are used
+primarily for calculating the dependence of energy and force components
+on the linear coefficients in the :doc:`pod pair_style <pair_pod>`,
+which is useful when training a POD potential to match target data. POD
+descriptors of an atom are characterized by the radial and angular
+distribution of neighbor atoms. The detailed mathematical definition is
+given in the papers by :ref:`(Nguyen and Rohskopf) <Nguyen20222c>`,
+:ref:`(Nguyen2023) <Nguyen20232c>`, :ref:`(Nguyen2024) <Nguyen20242c>`,
+and :ref:`(Nguyen and Sema) <Nguyen20243c>`.
+
+Compute *pod/atom* calculates the per-atom POD descriptors.
+
+Compute *podd/atom* calculates derivatives of the per-atom POD
+descriptors with respect to atom positions.
+
+Compute *pod/local* calculates the per-atom POD descriptors and their
+derivatives with respect to atom positions.
+
+Compute *pod/global* calculates the global POD descriptors and their
+derivatives with respect to atom positions.
+
+Examples how to use Compute POD commands are found in the directory
+``examples/PACKAGES/pod``.
+
+
+.. warning::
+
+   All of these compute styles produce *very* large per-atom output
+   arrays that scale with the total number of atoms in the system.
+   This will result in *very* large memory consumption for systems
+   with a large number of atoms.
+
+----------
+
+Output info
+"""""""""""
+
+Compute *pod/atom* produces an 2D array of size :math:`N \times M`,
+where :math:`N` is the number of atoms and :math:`M` is the number of
+descriptors. Each column corresponds to a particular POD descriptor.
+
+Compute *podd/atom* produces an 2D array of size :math:`N \times (M * 3
+N)`. Each column corresponds to a particular derivative of a POD
+descriptor.
+
+Compute *pod/local* produces an 2D array of size :math:`(1 + 3N) \times
+(M * N)`.  The first row contains the per-atom descriptors, and the last
+3N rows contain the derivatives of the per-atom descriptors with respect
+to atom positions.
+
+Compute *pod/global* produces an 2D array of size :math:`(1 + 3N) \times
+(M)`.  The first row contains the global descriptors, and the last 3N
+rows contain the derivatives of the global descriptors with respect to
+atom positions.
+
+Restrictions
+""""""""""""
+
+These computes are part of the ML-POD package.  They are only enabled
+if LAMMPS was built with that package.  See the :doc:`Build package
+<Build_package>` page for more info.
+
+Related commands
+""""""""""""""""
+
+:doc:`fitpod <fitpod_command>`,
+:doc:`pair_style pod <pair_pod>`
+
+
+Default
+"""""""
+
+none
+
+----------
+
+.. _Nguyen20222c:
+
+**(Nguyen and Rohskopf)** Nguyen and Rohskopf,  Journal of Computational Physics, 480, 112030, (2023).
+
+.. _Nguyen20232c:
+
+**(Nguyen2023)** Nguyen, Physical Review B, 107(14), 144103, (2023).
+
+.. _Nguyen20242c:
+
+**(Nguyen2024)** Nguyen, Journal of Computational Physics, 113102, (2024).
+
+.. _Nguyen20243c:
+
+**(Nguyen and Sema)** Nguyen and Sema, https://arxiv.org/abs/2405.00306, (2024).
+
+

doc/src/compute_rdf.rst

@@ -13,8 +13,8 @@ Syntax
 * ID, group-ID are documented in :doc:`compute <compute>` command
 * rdf = style name of this compute command
 * Nbin = number of RDF bins
-* itypeN = central atom type for Nth RDF histogram (see asterisk form below)
-* jtypeN = distribution atom type for Nth RDF histogram (see asterisk form below)
+* itypeN = central atom type for Nth RDF histogram (integer, type label, or asterisk form)
+* jtypeN = distribution atom type for Nth RDF histogram (integer, type label, or asterisk form)
 * zero or more keyword/value pairs may be appended
 * keyword = *cutoff*
 
@@ -96,14 +96,16 @@ is computed for :math:`g(r)` between all atom types.  If one or more pairs are
 listed, then a separate histogram is generated for each
 *itype*,\ *jtype* pair.
 
-The *itypeN* and *jtypeN* settings can be specified in one of two
-ways.  An explicit numeric value can be used, as in the fourth example
-above.  Or a wild-card asterisk can be used to specify a range of atom
-types.  This takes the form "\*" or "\*n" or "m\*" or "m\*n".  If
-:math:`N` is the number of atom types, then an asterisk with no numeric values
-means all types from 1 to :math:`N`. A leading asterisk means all types from 1
-to n (inclusive).  A trailing asterisk means all types from m to :math:`N`
-(inclusive).  A middle asterisk means all types from m to n (inclusive).
+The *itypeN* and *jtypeN* settings can be specified in one of three
+ways.  One or both of the types in the I,J pair can be a
+:doc:`type label <Howto_type_labels>`.  Or an explicit numeric value can be
+used, as in the fourth example above. Or a wild-card asterisk can be used
+to specify a range of atom types.  This takes the form "\*" or "\*n" or
+"m\*" or "m\*n". If :math:`N` is the number of atom types, then an asterisk
+with no numeric values means all types from 1 to :math:`N`. A leading
+asterisk means all types from 1 to n (inclusive).  A trailing asterisk
+means all types from m to :math:`N` (inclusive).  A middle asterisk means
+all types from m to n (inclusive).
 
 If both *itypeN* and *jtypeN* are single values, as in the fourth example
 above, this means that a :math:`g(r)` is computed where atoms of type *itypeN*

doc/src/compute_stress_mop.rst

@@ -126,7 +126,7 @@ These styles are part of the EXTRA-COMPUTE package. They are only
 enabled if LAMMPS is built with that package. See the :doc:`Build
 package <Build_package>` doc page on for more info.
 
-The method is only implemented for 3d orthogonal simulation boxes whose
+The method is implemented for orthogonal simulation boxes whose
 size does not change in time, and axis-aligned planes.
 
 The method only works with two-body pair interactions, because it

doc/src/create_atoms.rst

@@ -10,7 +10,7 @@ Syntax
 
    create_atoms type style args keyword values ...
 
-* type = atom type (1-Ntypes) of atoms to create (offset for molecule creation)
+* type = atom type (1-Ntypes or type label) of atoms to create (offset for molecule creation)
 * style = *box* or *region* or *single* or *mesh* or *random*
 
   .. parsed-literal::
@@ -37,7 +37,7 @@ Syntax
          seed = random # seed (positive integer)
        *basis* values = M itype
          M = which basis atom
-         itype = atom type (1-N) to assign to this basis atom
+         itype = atom type (1-Ntypes or type label) to assign to this basis atom
        *ratio* values = frac seed
          frac = fraction of lattice sites (0 to 1) to populate randomly
          seed = random # seed (positive integer)
@@ -74,6 +74,13 @@ Examples
 .. code-block:: LAMMPS
 
    create_atoms 1 box
+
+   labelmap atom 1 Pt
+   create_atoms Pt box
+
+   labelmap atom 1 C 2 Si
+   create_atoms C region regsphere basis Si C
+
    create_atoms 3 region regsphere basis 2 3
    create_atoms 3 region regsphere basis 2 3 ratio 0.5 74637
    create_atoms 3 single 0 0 5
@@ -468,12 +475,13 @@ to.
 
 The *overlap* keyword only applies to the *random* style.  It prevents
 newly created particles from being created closer than the specified
-*Doverlap* distance from any other particle.  When the particles being
-created are molecules, the radius of the molecule (from its geometric
-center) is added to *Doverlap*.  If particles have finite size (see
-:doc:`atom_style sphere <atom_style>` for example) *Doverlap* should
-be specified large enough to include the particle size in the
-non-overlapping criterion.
+*Doverlap* distance from any other particle.  If particles have finite
+size (see :doc:`atom_style sphere <atom_style>` for example) *Doverlap*
+should be specified large enough to include the particle size in the
+non-overlapping criterion.  If molecules are being randomly inserted, then
+an insertion is only accepted if each particle in the molecule meets the
+overlap criterion with respect to other particles (not including particles
+in the molecule itself).
 
 .. note::
 

doc/src/delete_bonds.rst

@@ -43,6 +43,9 @@ Examples
    delete_bonds all bond 0*3 special
    delete_bonds all stats
 
+   labelmap atom 4 hc
+   delete_bonds all atom hc special
+
 Description
 """""""""""
 
@@ -59,19 +62,20 @@ For all styles, by default, an interaction is only turned off (or on)
 if all the atoms involved are in the specified group.  See the *any*
 keyword to change the behavior.
 
-Several of the styles (\ *atom*, *bond*, *angle*, *dihedral*,
-*improper*\ ) take a *type* as an argument.  The specified *type* should
-be an integer from 0 to :math:`N`, where :math:`N` is the number of relevant
+Several of the styles (\ *atom*, *bond*, *angle*, *dihedral*, *improper*\ )
+take a *type* as an argument.  The specified *type* can be a
+:doc:`type label <Howto_type_labels>`.  Otherwise, the type should be an
+integer from 0 to :math:`N`, where :math:`N` is the number of relevant
 types (atom types, bond types, etc.).  A value of 0 is only relevant for
-style *bond*\ ; see details below.  In all cases, a wildcard asterisk
+style *bond*\ ; see details below.  For numeric types, a wildcard asterisk
 can be used in place of or in conjunction with the *type* argument to
 specify a range of types.  This takes the form "\*" or "\*n" or "m\*" or
-"m\*n".  If :math:`N` is the number of types, then an asterisk with no numeric
-values means all types from 0 to :math:`N`.  A leading asterisk means all
-types from 0 to n (inclusive).  A trailing asterisk means all types
-from m to N (inclusive).  A middle asterisk means all types from m to
-n (inclusive).  Note that it is fine to include a type of 0 for
-non-bond styles; it will simply be ignored.
+"m\*n".  If :math:`N` is the number of types, then an asterisk with no
+numeric values means all types from 0 to :math:`N`.  A leading asterisk
+means all types from 0 to n (inclusive).  A trailing asterisk means all
+types from m to N (inclusive).  A middle asterisk means all types from m to
+n (inclusive).  Note that it is fine to include a type of 0 for non-bond
+styles; it will simply be ignored.
 
 For style *multi* all bond, angle, dihedral, and improper interactions
 of any type, involving atoms in the group, are turned off.

doc/src/dump.rst

@@ -114,6 +114,7 @@ Syntax
            proc = ID of processor that owns atom
            procp1 = ID+1 of processor that owns atom
            type = atom type
+           typelabel = atom :doc:`type label <Howto_type_labels>`
            element = name of atom element, as defined by :doc:`dump_modify <dump_modify>` command
            mass = atom mass
            x,y,z = unscaled atom coordinates
@@ -470,8 +471,9 @@ followed by one line per atom with the atom type and the :math:`x`-,
 :math:`y`-, and :math:`z`-coordinate of that atom.  You can use the
 :doc:`dump_modify element <dump_modify>` option to change the output
 from using the (numerical) atom type to an element name (or some other
-label). This will help many visualization programs to guess bonds and
-colors.
+label). This option will help many visualization programs to guess bonds
+and colors. You can use the :doc:`dump_modify types labels <dump_modify>`
+option to replace numeric atom types with :doc:`type labels <Howto_type_labels>`.
 
 .. versionadded:: 22Dec2022
 
@@ -774,21 +776,21 @@ command creates a per-atom array with six columns:
 
 Per-atom attributes used as arguments to the *custom* and *cfg* styles:
 
-The *id*, *mol*, *proc*, *procp1*, *type*, *element*, *mass*, *vx*,
-*vy*, *vz*, *fx*, *fy*, *fz*, *q* attributes are self-explanatory.
+The *id*, *mol*, *proc*, *procp1*, *type*, *typelabel*, *element*, *mass*,
+*vx*, *vy*, *vz*, *fx*, *fy*, *fz*, *q* attributes are self-explanatory.
 
-*Id* is the atom ID.  *Mol* is the molecule ID, included in the data
-file for molecular systems.  *Proc* is the ID of the processor (0 to
+*Id* is the atom ID.  *Mol* is the molecule ID, included in the data file
+for molecular systems.  *Proc* is the ID of the processor (0 to
 :math:`N_\text{procs}-1`) that currently owns the atom.  *Procp1* is the
 proc ID+1, which can be convenient in place of a *type* attribute (1 to
 :math:`N_\text{types}`) for coloring atoms in a visualization program.
-*Type* is the atom type (1 to :math:`N_\text{types}`).  *Element* is
-typically the chemical name of an element, which you must assign to each
-type via the :doc:`dump_modify element <dump_modify>` command.  More
-generally, it can be any string you wish to associated with an atom
-type.  *Mass* is the atom mass. The quantities *vx*, *vy*, *vz*, *fx*,
-*fy*, *fz*, and *q* are components of atom velocity and force and atomic
-charge.
+*Type* is the atom type (1 to :math:`N_\text{types}`).  *Typelabel* is the
+atom :doc:`type label <Howto_type_labels>`.  *Element* is typically the
+chemical name of an element, which you must assign to each type via the
+:doc:`dump_modify element <dump_modify>` command.  More generally, it can
+be any string you wish to associated with an atom type.  *Mass* is the atom
+mass.  The quantities *vx*, *vy*, *vz*, *fx*, *fy*, *fz*, and *q* are
+components of atom velocity and force and atomic charge.
 
 There are several options for outputting atom coordinates.  The *x*,
 *y*, and *z* attributes write atom coordinates "unscaled", in the

doc/src/dump_modify.rst

@@ -17,7 +17,7 @@ Syntax
 * one or more keyword/value pairs may be appended
 
 * these keywords apply to various dump styles
-* keyword = *append* or *at* or *balance* or *buffer* or *colname* or *delay* or *element* or *every* or *every/time* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *skip* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *triclinic/general* or *units* or *unwrap*
+* keyword = *append* or *at* or *balance* or *buffer* or *colname* or *delay* or *element* or *every* or *every/time* or *fileper* or *first* or *flush* or *format* or *header* or *image* or *label* or *maxfiles* or *nfile* or *pad* or *pbc* or *precision* or *region* or *refresh* or *scale* or *sfactor* or *skip* or *sort* or *tfactor* or *thermo* or *thresh* or *time* or *triclinic/general* or *types* or *units* or *unwrap*
 
   .. parsed-literal::
 
@@ -81,6 +81,7 @@ Syntax
          these 3 args can be replaced by the word "none" to turn off thresholding
        *time* arg = *yes* or *no*
        *triclinic/general* arg = *yes* or *no*
+       *types* value = *numeric* or *labels*
        *units* arg = *yes* or *no*
        *unwrap* arg = *yes* or *no*
 
@@ -849,6 +850,13 @@ The default setting is *no*\ .
 
 ----------
 
+The *types* keyword applies only to the dump xyz style. If this keyword is
+used with a value of *numeric*, then numeric atom types are printed in the
+xyz file (default). If the value *labels* is specified, then
+:doc:`type labels <Howto_type_labels>` are printed for atom types.
+
+----------
+
 The *triclinic/general* keyword only applies to the dump *atom* and
 *custom* styles.  It can only be used with a value of *yes* if the
 simulation box was created as a general triclinic box.  See the
@@ -960,11 +968,11 @@ The option defaults are
 * sort = id for dump styles *dcd*, *xtc*, and *xyz*
 * thresh = none
 * time = no
-* triclinic/general no
+* triclinic/general = no
+* types = numeric
 * units = no
 * unwrap = no
 
 * compression_level = 9 (gz variants)
 * compression_level = 0 (zstd variants)
 * checksum = yes (zstd variants)
-

doc/src/fitpod_command.rst

@@ -1,18 +1,19 @@
 .. index:: fitpod
 
 fitpod command
-======================
+==============
 
 Syntax
 """"""
 
 .. code-block:: LAMMPS
 
-   fitpod Ta_param.pod Ta_data.pod
+   fitpod Ta_param.pod Ta_data.pod Ta_coefficients.pod
 
 * fitpod = style name of this command
 * Ta_param.pod = an input file that describes proper orthogonal descriptors (PODs)
 * Ta_data.pod = an input file that specifies DFT data used to fit a POD potential
+* Ta_coefficients.pod (optional) = an input file that specifies trainable coefficients of a POD potential
 
 Examples
 """"""""
@@ -20,20 +21,26 @@ Examples
 .. code-block:: LAMMPS
 
    fitpod Ta_param.pod Ta_data.pod
+   fitpod Ta_param.pod Ta_data.pod Ta_coefficients.pod
 
 Description
 """""""""""
 .. versionadded:: 22Dec2022
 
 Fit a machine-learning interatomic potential (ML-IAP) based on proper
-orthogonal descriptors (POD).  Two input files are required for this
-command. The first input file describes a POD potential parameter
-settings, while the second input file specifies the DFT data used for
-the fitting procedure.
+orthogonal descriptors (POD); please see :ref:`(Nguyen and Rohskopf)
+<Nguyen20222a>`, :ref:`(Nguyen2023) <Nguyen20232a>`, :ref:`(Nguyen2024)
+<Nguyen20242a>`, and :ref:`(Nguyen and Sema) <Nguyen20243a>` for details.
+The fitted POD potential can be used to run MD simulations via
+:doc:`pair_style pod <pair_pod>`.
 
-The table below has one-line descriptions of all the keywords that can
-be used in the first input file (i.e. ``Ta_param.pod`` in the example
-above):
+Two input files are required for this command. The first input file
+describes a POD potential parameter settings, while the second input
+file specifies the DFT data used for the fitting procedure. All keywords
+except *species* have default values. If a keyword is not set in the
+input file, its default value is used. The table below has one-line
+descriptions of all the keywords that can be used in the first input
+file (i.e. ``Ta_param.pod``)
 
 .. list-table::
    :header-rows: 1
@@ -52,7 +59,7 @@ above):
      - INT
      - three integer constants specify boundary conditions
    * - rin
-     - 1.0
+     - 0.5
      - REAL
      - a real number specifies the inner cut-off radius
    * - rcut
@@ -60,46 +67,75 @@ above):
      - REAL
      - a real number specifies the outer cut-off radius
    * - bessel_polynomial_degree
-     - 3
+     - 4
      - INT
      - the maximum degree of Bessel polynomials
    * - inverse_polynomial_degree
-     - 6
+     - 8
      - INT
      - the maximum degree of inverse radial basis functions
+   * - number_of_environment_clusters
+     - 1
+     - INT
+     - the number of clusters for environment-adaptive potentials
+   * - number_of_principal_components
+     - 2
+     - INT
+     - the number of principal components for dimensionality reduction
    * - onebody
      - 1
      - BOOL
      - turns on/off one-body potential
    * - twobody_number_radial_basis_functions
-     - 6
+     - 8
      - INT
      - number of radial basis functions for two-body potential
    * - threebody_number_radial_basis_functions
-     - 5
+     - 6
      - INT
      - number of radial basis functions for three-body potential
-   * - threebody_number_angular_basis_functions
+   * - threebody_angular_degree
      - 5
      - INT
-     - number of angular basis functions for three-body potential
-   * - fourbody_snap_twojmax
+     - angular degree for three-body potential
+   * - fourbody_number_radial_basis_functions
+     - 4
+     - INT
+     - number of radial basis functions for four-body potential
+   * - fourbody_angular_degree
+     - 3
+     - INT
+     - angular degree for four-body potential
+   * - fivebody_number_radial_basis_functions
      - 0
      - INT
-     - band limit for SNAP bispectrum components (0,2,4,6,8... allowed)
-   * - fourbody_snap_chemflag
+     - number of radial basis functions for five-body potential
+   * - fivebody_angular_degree
      - 0
-     - BOOL
-     - turns on/off the explicit multi-element variant of the SNAP bispectrum components
-   * - quadratic_pod_potential
+     - INT
+     - angular degree for five-body potential
+   * - sixbody_number_radial_basis_functions
      - 0
-     - BOOL
-     - turns on/off quadratic POD potential
+     - INT
+     - number of radial basis functions for six-body potential
+   * - sixbody_angular_degree
+     - 0
+     - INT
+     - angular degree for six-body potential
+   * - sevenbody_number_radial_basis_functions
+     - 0
+     - INT
+     - number of radial basis functions for seven-body potential
+   * - sevenbody_angular_degree
+     - 0
+     - INT
+     - angular degree for seven-body potential
+
+Note that both the number of radial basis functions and angular degree
+must decrease as the body order increases. The next table describes all
+keywords that can be used in the second input file (i.e. ``Ta_data.pod``
+in the example above):
 
-All keywords except *species* have default values. If a keyword is not
-set in the input file, its default value is used.  The next table
-describes all keywords that can be used in the second input file
-(i.e. ``Ta_data.pod`` in the example above):
 
 .. list-table::
    :header-rows: 1
@@ -125,6 +161,10 @@ describes all keywords that can be used in the second input file
      - ""
      - STRING
      - specifies the path to test data files in double quotes
+   * - path_to_environment_configuration_set
+     - ""
+     - STRING
+     - specifies the path to environment configuration files in double quotes
    * - fraction_training_data_set
      - 1.0
      - REAL
@@ -133,6 +173,14 @@ describes all keywords that can be used in the second input file
      - 0
      - BOOL
      - turns on/off randomization of the training set
+   * - fraction_test_data_set
+     - 1.0
+     - REAL
+     - a real number (<= 1.0) specifies the fraction of the test set used to validate POD
+   * - randomize_test_data_set
+     - 0
+     - BOOL
+     - turns on/off randomization of the test set
    * - fitting_weight_energy
      - 100.0
      - REAL
@@ -161,6 +209,10 @@ describes all keywords that can be used in the second input file
      - 8
      - INT
      - number of digits after the decimal points for numbers in the coefficient file
+   * - group_weights
+     - global
+     - STRING
+     - ``table`` uses group weights defined for each group named by filename
 
 All keywords except *path_to_training_data_set* have default values. If
 a keyword is not set in the input file, its default value is used.  After
@@ -172,14 +224,44 @@ successful training, a number of output files are produced, if enabled:
 * ``<basename>_test_analysis.pod`` reports detailed errors for all test configurations
 * ``<basename>_coefficients.pod`` contains the coefficients of the POD potential
 
-After training the POD potential, ``Ta_param.pod`` and ``<basename>_coefficients.pod``
-are the two files needed to use the POD potential in LAMMPS. See
-:doc:`pair_style pod <pair_pod>` for using the POD potential. Examples
-about training and using POD potentials are found in the directory
-lammps/examples/PACKAGES/pod.
+After training the POD potential, ``Ta_param.pod`` and
+``<basename>_coefficients.pod`` are the two files needed to use the POD
+potential in LAMMPS.  See :doc:`pair_style pod <pair_pod>` for using the
+POD potential. Examples about training and using POD potentials are
+found in the directory lammps/examples/PACKAGES/pod and the Github repo
+https://github.com/cesmix-mit/pod-examples.
 
-Parameterized Potential Energy Surface
-""""""""""""""""""""""""""""""""""""""
+Loss Function Group Weights
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The *group_weights* keyword in the ``data.pod`` file is responsible for
+weighting certain groups of configurations in the loss function. For
+example:
+
+.. code-block:: LAMMPS
+
+    group_weights table
+    Displaced_A15 100.0 1.0
+    Displaced_BCC 100.0 1.0
+    Displaced_FCC 100.0 1.0
+    Elastic_BCC   100.0 1.0
+    Elastic_FCC   100.0 1.0
+    GSF_110       100.0 1.0
+    GSF_112       100.0 1.0
+    Liquid        100.0 1.0
+    Surface       100.0 1.0
+    Volume_A15    100.0 1.0
+    Volume_BCC    100.0 1.0
+    Volume_FCC    100.0 1.0
+
+This will apply an energy weight of ``100.0`` and a force weight of
+``1.0`` for all groups in the ``Ta`` example. The groups are named by
+their respective filename. If certain groups are left out of this table,
+then the globally defined weights from the ``fitting_weight_energy`` and
+``fitting_weight_force`` keywords will be used.
+
+POD Potential
+"""""""""""""
 
 We consider a multi-element system of *N* atoms with :math:`N_{\rm e}`
 unique elements.  We denote by :math:`\boldsymbol r_n` and :math:`Z_n`
@@ -187,535 +269,82 @@ position vector and type of an atom *n* in the system,
 respectively. Note that we have :math:`Z_n \in \{1, \ldots, N_{\rm e}
 \}`, :math:`\boldsymbol R = (\boldsymbol r_1, \boldsymbol r_2, \ldots,
 \boldsymbol r_N) \in \mathbb{R}^{3N}`, and :math:`\boldsymbol Z = (Z_1,
-Z_2, \ldots, Z_N) \in \mathbb{N}^{N}`. The potential energy surface
-(PES) of the system can be expressed as a many-body expansion of the
-form
+Z_2, \ldots, Z_N) \in \mathbb{N}^{N}`. The total energy of the
+POD potential is expressed as :math:`E(\boldsymbol R, \boldsymbol Z) =
+\sum_{i=1}^N E_i(\boldsymbol R_i, \boldsymbol Z_i)`, where
 
 .. math::
 
-    E(\boldsymbol R, \boldsymbol Z, \boldsymbol{\eta}, \boldsymbol{\mu}) \ = \ & \sum_{i} V^{(1)}(\boldsymbol r_i, Z_i, \boldsymbol \mu^{(1)} ) + \frac12 \sum_{i,j} V^{(2)}(\boldsymbol r_i, \boldsymbol r_j, Z_i, Z_j, \boldsymbol \eta, \boldsymbol \mu^{(2)})  \\
-    & + \frac16 \sum_{i,j,k} V^{(3)}(\boldsymbol r_i, \boldsymbol r_j, \boldsymbol r_k, Z_i, Z_j, Z_k, \boldsymbol \eta, \boldsymbol \mu^{(3)}) + \ldots
-
-where :math:`V^{(1)}` is the one-body potential often used for
-representing external field or energy of isolated elements, and the
-higher-body potentials :math:`V^{(2)}, V^{(3)}, \ldots` are symmetric,
-uniquely defined, and zero if two or more indices take identical values.
-The superscript on each potential denotes its body order. Each *q*-body
-potential :math:`V^{(q)}` depends on :math:`\boldsymbol \mu^{(q)}` which
-are sets of parameters to fit the PES. Note that :math:`\boldsymbol \mu`
-is a collection of all potential parameters :math:`\boldsymbol
-\mu^{(1)}`, :math:`\boldsymbol \mu^{(2)}`, :math:`\boldsymbol
-\mu^{(3)}`, etc, and that :math:`\boldsymbol \eta` is a set of
-hyper-parameters such as inner cut-off radius :math:`r_{\rm in}` and
-outer cut-off radius :math:`r_{\rm cut}`.
-
-Interatomic potentials rely on parameters to learn relationship between
-atomic environments and interactions.  Since interatomic potentials are
-approximations by nature, their parameters need to be set to some
-reference values or fitted against data by necessity.  Typically,
-potential fitting finds optimal parameters, :math:`\boldsymbol \mu^*`,
-to minimize a certain loss function of the predicted quantities and
-data. Since the fitted potential depends on the data set used to fit it,
-different data sets will yield different optimal parameters and thus
-different fitted potentials. When fitting the same functional form on
-*Q* different data sets, we would obtain *Q* different optimized
-potentials, :math:`E(\boldsymbol R,\boldsymbol Z, \boldsymbol \eta,
-\boldsymbol \mu_q^*), 1 \le q \le Q`.  Consequently, there exist many
-different sets of optimized parameters for empirical interatomic
-potentials.
-
-Instead of optimizing the potential parameters, inspired by the reduced
-basis method :ref:`(Grepl) <Grepl20072>` for parameterized partial
-differential equations, we view the parameterized PES as a parametric
-manifold of potential energies
-
-.. math::
-
-    \mathcal{M} = \{E(\boldsymbol R, \boldsymbol Z, \boldsymbol \eta, \boldsymbol \mu) \ | \  \boldsymbol \mu \in \Omega^{\boldsymbol \mu} \}
-
-where :math:`\Omega^{\boldsymbol \mu}` is a parameter domain in which
-:math:`\boldsymbol \mu` resides.  The parametric manifold
-:math:`\mathcal{M}` contains potential energy surfaces for all values of
-:math:`\boldsymbol \mu \in \Omega^{\boldsymbol \mu}`.  Therefore, the
-parametric manifold yields a much richer and more transferable atomic
-representation than any particular individual PES :math:`E(\boldsymbol
-R, \boldsymbol Z, \boldsymbol \eta, \boldsymbol \mu^*)`.
-
-We propose specific forms of the parameterized potentials for one-body,
-two-body, and three-body interactions. We apply the Karhunen-Loeve
-expansion to snapshots of the parameterized potentials to obtain sets of
-orthogonal basis functions. These basis functions are aggregated
-according to the chemical elements of atoms, thus leading to
-multi-element proper orthogonal descriptors.
-
-Proper Orthogonal Descriptors
-"""""""""""""""""""""""""""""
-
-Proper orthogonal descriptors are finger prints characterizing the
-radial and angular distribution of a system of atoms. The detailed
-mathematical definition is given in the paper by Nguyen and Rohskopf
-:ref:`(Nguyen) <Nguyen20222>`.
-
-The descriptors for the one-body interaction are used to capture energy
-of isolated elements and defined as follows
-
-.. math::
-
-    D_{ip}^{(1)} =  \left\{
-        \begin{array}{ll}
-        1, & \mbox{if } Z_i = p \\
-        0, & \mbox{if } Z_i \neq p
-        \end{array}
-    \right.
-
-for :math:`1 \le i \le N, 1 \le p \le N_{\rm e}`. The number of one-body
-descriptors per atom is equal to the number of elements. The one-body
-descriptors are independent of atom positions, but dependent on atom
-types. The one-body descriptors are active only when the keyword
-*onebody* is set to 1.
-
-We adopt the usual assumption that the direct interaction between two
-atoms vanishes smoothly when their distance is greater than the outer
-cutoff distance :math:`r_{\rm cut}`. Furthermore, we assume that two
-atoms can not get closer than the inner cutoff distance :math:`r_{\rm
-in}` due to the Pauli repulsion principle. Let :math:`r \in (r_{\rm in},
-r_{\rm cut})`, we introduce the following parameterized radial functions
-
-.. math::
-
-    \phi(r, r_{\rm in}, r_{\rm cut}, \alpha, \beta)  = \frac{\sin (\alpha \pi x) }{r - r_{\rm in}}, \qquad  \varphi(r, \gamma)  = \frac{1}{r^\gamma} ,
-
-where the scaled distance function :math:`x` is defined below to enrich the two-body manifold
-
-.. math::
-
-    x(r, r_{\rm in}, r_{\rm cut}, \beta) = \frac{e^{-\beta(r - r_{\rm in})/(r_{\rm cut} - r_{\rm in})} - 1}{e^{-\beta} - 1} .
-
-We introduce the following function as a convex combination of the two functions
-
-.. math::
-
-    \psi(r, r_{\rm in}, r_{\rm cut}, \alpha, \beta, \gamma, \kappa)  = \kappa \phi(r, r_{\rm in}, r_{\rm cut}, \alpha, \beta) + (1- \kappa)  \varphi(r, \gamma) .
-
-We see that :math:`\psi` is a function of distance :math:`r`, cut-off
-distances :math:`r_{\rm in}` and :math:`r_{\rm cut}`, and parameters
-:math:`\alpha, \beta, \gamma, \kappa`. Together these parameters allow
-the function :math:`\psi` to characterize a diverse spectrum of two-body
-interactions within the cut-off interval :math:`(r_{\rm in}, r_{\rm
-cut})`.
-
-Next, we introduce the following parameterized potential
-
-.. math::
-
-    W^{(2)}(r_{ij}, \boldsymbol \eta, \boldsymbol \mu^{(2)})  = f_{\rm c}(r_{ij}, \boldsymbol \eta) \psi(r_{ij}, \boldsymbol \eta, \boldsymbol \mu^{(2)})
-
-where :math:`\eta_1 = r_{\rm in}, \eta_2 = r_{\rm cut}, \mu_1^{(2)} =
-\alpha, \mu_2^{(2)} = \beta, \mu_3^{(2)} = \gamma`, and
-:math:`\mu_4^{(2)} = \kappa`. Here the cut-off function :math:`f_{\rm
-c}(r_{ij}, \boldsymbol \eta)` proposed in [refs] is used to ensure the
-smooth vanishing of the potential and its derivative for :math:`r_{ij}
-\ge r_{\rm cut}`:
-
-.. math::
-
-    f_{\rm c}(r_{ij},  r_{\rm in}, r_{\rm cut})  =  \exp \left(1 -\frac{1}{\sqrt{\left(1 - \frac{(r-r_{\rm in})^3}{(r_{\rm cut} - r_{\rm in})^3} \right)^2 + 10^{-6}}} \right)
-
-Based on the parameterized potential, we form a set of snapshots as
-follows.  We assume that we are given :math:`N_{\rm s}` parameter tuples
-:math:`\boldsymbol \mu^{(2)}_\ell, 1 \le \ell \le N_{\rm s}`. We
-introduce the following set of snapshots on :math:`(r_{\rm in}, r_{\rm
-cut})`:
-
-.. math::
-
-    \xi_\ell(r_{ij}, \boldsymbol \eta) =  W^{(2)}(r_{ij}, \boldsymbol \eta, \boldsymbol \mu^{(2)}_\ell),  \quad \ell = 1, \ldots, N_{\rm s} .
-
-To ensure adequate sampling of the PES for different parameters, we
-choose :math:`N_{\rm s}` parameter points :math:`\boldsymbol
-\mu^{(2)}_\ell = (\alpha_\ell, \beta_\ell, \gamma_\ell, \kappa_\ell), 1
-\le \ell \le N_{\rm s}` as follows. The parameters :math:`\alpha \in [1,
-N_\alpha]` and :math:`\gamma \in [1, N_\gamma]` are integers, where
-:math:`N_\alpha` and :math:`N_\gamma` are the highest degrees for
-:math:`\alpha` and :math:`\gamma`, respectively. We next choose
-:math:`N_\beta` different values of :math:`\beta` in the interval
-:math:`[\beta_{\min}, \beta_{\max}]`, where :math:`\beta_{\min} = 0` and
-:math:`\beta_{\max} = 4`. The parameter :math:`\kappa` can be set either
-0 or 1.  Hence, the total number of parameter points is :math:`N_{\rm s}
-= N_\alpha N_\beta + N_\gamma`.  Although :math:`N_\alpha, N_\beta,
-N_\gamma` can be chosen conservatively large, we find that
-:math:`N_\alpha = 6, N_\beta = 3, N_\gamma = 8` are adequate for most
-problems.  Note that :math:`N_\alpha` and :math:`N_\gamma` correspond to
-*bessel_polynomial_degree* and *inverse_polynomial_degree*,
-respectively.
-
-We employ the Karhunen-Loeve (KL) expansion to generate an orthogonal
-basis set which is known to be optimal for representation of the
-snapshot family :math:`\{\xi_\ell\}_{\ell=1}^{N_{\rm s}}`. The two-body
-orthogonal basis functions are computed as follows
-
-.. math::
-
-    U^{(2)}_m(r_{ij}, \boldsymbol \eta) = \sum_{\ell = 1}^{N_{\rm s}} A_{\ell m}(\boldsymbol \eta) \,  \xi_\ell(r_{ij}, \boldsymbol \eta), \qquad m = 1, \ldots, N_{\rm 2b} ,
-
-where the matrix :math:`\boldsymbol A \in \mathbb{R}^{N_{\rm s} \times
-N_{\rm s}}` consists of eigenvectors of the eigenvalue problem
-
-.. math::
-
-    \boldsymbol C \boldsymbol a = \lambda \boldsymbol a
-
-with the entries of :math:`\boldsymbol C \in \mathbb{R}^{N_{\rm s} \times N_{\rm s}}` being given by
-
-.. math::
-
-    C_{ij}  = \frac{1}{N_{\rm s}} \int_{r_{\rm in}}^{r_{\rm cut}} \xi_i(x, \boldsymbol \eta) \xi_j(x, \boldsymbol \eta) dx, \quad 1 \le i, j \le N_{\rm s}
-
-Note that the eigenvalues :math:`\lambda_\ell, 1 \le \ell \le N_{\rm
-s}`, are ordered such that :math:`\lambda_1 \ge \lambda_2 \ge \ldots \ge
-\lambda_{N_{\rm s}}`, and that the matrix :math:`\boldsymbol A` is
-pe-computed and stored for any given :math:`\boldsymbol \eta`.  Owing to
-the rapid convergence of the KL expansion, only a small number of
-orthogonal basis functions is needed to obtain accurate
-approximation. The value of :math:`N_{\rm 2b}` corresponds to
-*twobody_number_radial_basis_functions*.
-
-The two-body proper orthogonal descriptors at each atom *i* are computed
-by summing the orthogonal basis functions over the neighbors of atom *i*
-and numerating on the atom types as follows
-
-.. math::
-
-    D^{(2)}_{im l(p, q) }(\boldsymbol \eta)  = \left\{
-    \begin{array}{ll}
-    \displaystyle \sum_{\{j | Z_j = q\}} U^{(2)}_m(r_{ij},  \boldsymbol \eta), & \mbox{if } Z_i = p \\
-    0, & \mbox{if } Z_i \neq p
-    \end{array}
-    \right.
-
-for :math:`1 \le i \le N, 1 \le m \le N_{\rm 2b}, 1 \le q, p \le N_{\rm
-e}`. Here :math:`l(p,q)` is a symmetric index mapping such that
-
-.. math::
-
-    l(p,q)  = \left\{
-    \begin{array}{ll}
-    q + (p-1) N_{\rm e} - p(p-1)/2, & \mbox{if } q \ge p \\
-    p + (q-1) N_{\rm e} - q(q-1)/2, & \mbox{if } q < p .
-    \end{array}
-    \right.
-
-The number of two-body descriptors per atom is thus :math:`N_{\rm 2b}
-N_{\rm e}(N_{\rm e}+1)/2`.
-
-It is important to note that the orthogonal basis functions do not
-depend on the atomic numbers :math:`Z_i` and :math:`Z_j`. Therefore, the
-cost of evaluating the basis functions and their derivatives with
-respect to :math:`r_{ij}` is independent of the number of elements
-:math:`N_{\rm e}`. Consequently, even though the two-body proper
-orthogonal descriptors depend on :math:`\boldsymbol Z`, their
-computational complexity is independent of :math:`N_{\rm e}`.
-
-In order to provide proper orthogonal descriptors for three-body
-interactions, we need to introduce a three-body parameterized
-potential. In particular, the three-body potential is defined as a
-product of radial and angular functions as follows
-
-.. math::
-
-    W^{(3)}(r_{ij}, r_{ik}, \theta_{ijk}, \boldsymbol \eta, \boldsymbol \mu^{(3)})  =  \psi(r_{ij}, r_{\rm min}, r_{\rm max}, \alpha, \beta, \gamma, \kappa) f_{\rm c}(r_{ij}, r_{\rm min}, r_{\rm max}) \\
-    \psi(r_{ik}, r_{\rm min}, r_{\rm max}, \alpha, \beta, \gamma, \kappa) f_{\rm c}(r_{ik}, r_{\rm min}, r_{\rm max}) \\
-    \cos (\sigma \theta_{ijk} + \zeta)
-
-where :math:`\sigma` is the periodic multiplicity, :math:`\zeta` is the
-equilibrium angle, :math:`\boldsymbol \mu^{(3)} = (\alpha, \beta,
-\gamma, \kappa, \sigma, \zeta)`. The three-body potential provides an
-angular fingerprint of the atomic environment through the bond angles
-:math:`\theta_{ijk}` formed with each pair of neighbors :math:`j` and
-:math:`k`.  Compared to the two-body potential, the three-body potential
-has two extra parameters :math:`(\sigma, \zeta)` associated with the
-angular component.
-
-Let :math:`\boldsymbol \varrho = (\alpha, \beta, \gamma, \kappa)`. We
-assume that we are given :math:`L_{\rm r}` parameter tuples
-:math:`\boldsymbol \varrho_\ell, 1 \le \ell \le L_{\rm r}`.  We
-introduce the following set of snapshots on :math:`(r_{\min},
-r_{\max})`:
-
-.. math::
-
-    \zeta_\ell(r_{ij}, r_{\rm min}, r_{\rm max} ) =  \psi(r_{ij}, r_{\rm min}, r_{\rm max}, \boldsymbol \varrho_\ell) f_{\rm c}(r_{ij}, r_{\rm min},  r_{\rm max}), \quad 1 \le \ell \le L_{\rm r} .
-
-We apply the Karhunen-Loeve (KL) expansion to this set of snapshots to
-obtain orthogonal basis functions as follows
-
-.. math::
-
-    U^{r}_m(r_{ij}, r_{\rm min}, r_{\rm max} ) = \sum_{\ell = 1}^{L_{\rm r}} A_{\ell m} \,  \zeta_\ell(r_{ij}, r_{\rm min}, r_{\rm max} ), \qquad m = 1, \ldots, N_{\rm r} ,
-
-where the matrix :math:`\boldsymbol A \in \mathbb{R}^{L_{\rm r} \times L_{\rm r}}` consists
-of eigenvectors of the eigenvalue problem. For the parameterized angular function,
-we consider angular basis functions
-
-.. math::
-
-    U^{a}_n(\theta_{ijk}) = \cos ((n-1) \theta_{ijk}), \qquad  n = 1,\ldots, N_{\rm a},
-
-where :math:`N_{\rm a}` is the number of angular basis functions. The orthogonal
-basis functions for the parameterized potential are computed as follows
-
-.. math::
-
-    U^{(3)}_{mn}(r_{ij}, r_{ik}, \theta_{ijk}, \boldsymbol \eta) = U^{r}_m(r_{ij}, \boldsymbol \eta) U^{r}_m(r_{ik}, \boldsymbol \eta) U^{a}_n(\theta_{ijk}),
-
-for :math:`1 \le m \le N_{\rm r}, 1 \le n \le N_{\rm a}`. The number of three-body
-orthogonal basis functions is equal to :math:`N_{\rm 3b} = N_{\rm r} N_{\rm a}` and
-independent of the number of elements. The value of :math:`N_{\rm r}` corresponds to
-*threebody_number_radial_basis_functions*, while that of :math:`N_{\rm a}` to
-*threebody_number_angular_basis_functions*.
-
-The three-body proper orthogonal descriptors at each atom *i*
-are obtained by summing over the neighbors *j* and *k* of atom *i* as
-
-.. math::
-
-    D^{(3)}_{imn \ell(p, q, s)}(\boldsymbol \eta)  = \left\{
-    \begin{array}{ll}
-    \displaystyle \sum_{\{j | Z_j = q\}} \sum_{\{k | Z_k = s\}} U^{(3)}_{mn}(r_{ij}, r_{ik}, \theta_{ijk}, \boldsymbol \eta), & \mbox{if } Z_i = p \\
-    0, & \mbox{if } Z_i \neq p
-    \end{array}
-    \right.
-
-for :math:`1 \le i \le N, 1 \le m \le N_{\rm r}, 1 \le n \le N_{\rm a}, 1 \le q, p, s \le N_{\rm e}`,
-where
-
-.. math::
-
-    \ell(p,q,s)  = \left\{
-    \begin{array}{ll}
-    s + (q-1) N_{\rm e} - q(q-1)/2 + (p-1)N_{\rm e}(1+N_{\rm e})/2 , & \mbox{if } s \ge q \\
-    q + (s-1) N_{\rm e} - s(s-1)/2 + (p-1)N_{\rm e}(1+N_{\rm e})/2, & \mbox{if } s < q .
-    \end{array}
-    \right.
-
-The number of three-body descriptors per atom is thus :math:`N_{\rm 3b} N_{\rm e}^2(N_{\rm e}+1)/2`.
-While the number of three-body PODs is cubic function of the number of elements,
-the computational complexity of the three-body PODs is independent of the number of elements.
-
-Four-Body SNAP Descriptors
-""""""""""""""""""""""""""
-
-In addition to the proper orthogonal descriptors described above, we also employ
-the spectral neighbor analysis potential (SNAP) descriptors. SNAP uses bispectrum components
-to characterize the local neighborhood of each atom in a very general way. The mathematical definition
-of the bispectrum calculation and its derivatives w.r.t. atom positions is described in
-:doc:`compute snap <compute_sna_atom>`. In SNAP, the
-total energy is decomposed into a sum over atom energies. The energy of
-atom *i* is expressed as a weighted sum over bispectrum components.
-
-.. math::
-
-   E_i^{\rm SNAP} = \sum_{k=1}^{N_{\rm 4b}} \sum_{p=1}^{N_{\rm e}} c_{kp}^{(4)} D_{ikp}^{(4)}
-
-
-where the SNAP descriptors are related to the bispectrum components by
-
-.. math::
-
-    D^{(4)}_{ikp}  = \left\{
-    \begin{array}{ll}
-    \displaystyle B_{ik}, & \mbox{if } Z_i = p \\
-    0, & \mbox{if } Z_i \neq p
-    \end{array}
-    \right.
-
-Here :math:`B_{ik}` is the *k*\ -th bispectrum component of atom *i*. The number of
-bispectrum components :math:`N_{\rm 4b}` depends on the value of *fourbody_snap_twojmax* :math:`= 2 J_{\rm max}`
-and *fourbody_snap_chemflag*. If *fourbody_snap_chemflag* = 0
-then :math:`N_{\rm 4b} = (J_{\rm max}+1)(J_{\rm max}+2)(J_{\rm max}+1.5)/3`.
-If *fourbody_snap_chemflag* = 1 then :math:`N_{\rm 4b} = N_{\rm e}^3 (J_{\rm max}+1)(J_{\rm max}+2)(J_{\rm max}+1.5)/3`.
-The bispectrum calculation is described in more detail in :doc:`compute sna/atom <compute_sna_atom>`.
-
-Linear Proper Orthogonal Descriptor Potentials
-""""""""""""""""""""""""""""""""""""""""""""""
-
-The proper orthogonal descriptors and SNAP descriptors are used to define the atomic energies
-in the following expansion
-
-.. math::
-
-    E_{i}(\boldsymbol \eta) = \sum_{p=1}^{N_{\rm e}} c^{(1)}_p D^{(1)}_{ip} + \sum_{m=1}^{N_{\rm 2b}}  \sum_{l=1}^{N_{\rm e}(N_{\rm e}+1)/2} c^{(2)}_{ml} D^{(2)}_{iml}(\boldsymbol \eta) + \sum_{m=1}^{N_{\rm r}} \sum_{n=1}^{N_{\rm a}}  \sum_{\ell=1}^{N_{\rm e}^2(N_{\rm e}+1)/2} c^{(3)}_{mn\ell} D^{(3)}_{imn\ell}(\boldsymbol \eta) + \sum_{k=1}^{N_{\rm 4b}} \sum_{p=1}^{N_{\rm e}} c_{kp}^{(4)} D_{ikp}^{(4)}(\boldsymbol \eta),
-
-where :math:`D^{(1)}_{ip}, D^{(2)}_{iml}, D^{(3)}_{imn\ell}, D^{(4)}_{ikp}` are the  one-body, two-body, three-body, four-body descriptors,
-respectively, and :math:`c^{(1)}_p, c^{(2)}_{ml}, c^{(3)}_{mn\ell}, c^{(4)}_{kp}` are their respective expansion
-coefficients. In a more compact notation that implies summation over descriptor indices
-the atomic energies can be written as
-
-.. math::
-
-    E_i(\boldsymbol \eta) =  \sum_{m=1}^{N_{\rm e}} c^{(1)}_m D^{(1)}_{im} +  \sum_{m=1}^{N_{\rm d}^{(2)}} c^{(2)}_k D^{(2)}_{im} + \sum_{m=1}^{N_{\rm d}^{(3)}} c^{(3)}_m D^{(3)}_{im} + \sum_{m=1}^{N_{\rm d}^{(4)}} c^{(4)}_m D^{(4)}_{im}
-
-where :math:`N_{\rm d}^{(2)} = N_{\rm 2b} N_{\rm e} (N_{\rm e}+1)/2`,
-:math:`N_{\rm d}^{(3)} = N_{\rm 3b} N_{\rm e}^2 (N_{\rm e}+1)/2`, and
-:math:`N_{\rm d}^{(4)} = N_{\rm 4b} N_{\rm e}` are
-the number of two-body, three-body, and four-body descriptors, respectively.
-
-The potential energy is then obtained by summing local atomic energies :math:`E_i`
-for all atoms :math:`i` in the system
-
-.. math::
-
-    E(\boldsymbol \eta) = \sum_{i}^N E_{i}(\boldsymbol \eta)
-
-Because the descriptors are one-body, two-body, and three-body terms,
-the resulting POD potential is a three-body PES. We can express the potential
-energy as a linear combination of the global descriptors as follows
-
-.. math::
-
-    E(\boldsymbol \eta) = \sum_{m=1}^{N_{\rm e}} c^{(1)}_m d^{(1)}_{m} +  \sum_{m=1}^{N_{\rm d}^{(2)}} c^{(2)}_m d^{(2)}_{m} + \sum_{m=1}^{N_{\rm d}^{(3)}} c^{(3)}_m d^{(3)}_{m} + \sum_{m=1}^{N_{\rm d}^{(4)}} c^{(4)}_m d^{(4)}_{m}
-
-where  the global descriptors are given by
-
-.. math::
-
-    d_{m}^{(1)}(\boldsymbol \eta) = \sum_{i=1}^N D_{im}^{(1)}(\boldsymbol \eta), \quad d_{m}^{(2)}(\boldsymbol \eta) = \sum_{i=1}^N D_{im}^{(2)}(\boldsymbol \eta), \quad d_{m}^{(3)}(\boldsymbol \eta) = \sum_{i=1}^N D_{im}^{(3)}(\boldsymbol \eta), \quad d_{m}^{(4)}(\boldsymbol \eta) = \sum_{i=1}^N D_{im}^{(4)}(\boldsymbol \eta)
-
-Hence, we obtain the atomic forces as
-
-.. math::
-
-    \boldsymbol F = -\nabla E(\boldsymbol \eta) = - \sum_{m=1}^{N_{\rm d}^{(2)}}  c^{(2)}_m  \nabla d_m^{(2)} - \sum_{m=1}^{N_{\rm d}^{(3)}}  c^{(3)}_m \nabla d_m^{(3)} - \sum_{m=1}^{N_{\rm d}^{(4)}}  c^{(4)}_m \nabla d_m^{(4)}
-
-where :math:`\nabla d_m^{(2)}`, :math:`\nabla d_m^{(3)}` and :math:`\nabla d_m^{(4)}` are derivatives of the two-body
-three-body, and four-body global descriptors with respect to atom positions, respectively.
-Note that since the first-body global descriptors are constant, their derivatives are zero.
-
-Quadratic Proper Orthogonal Descriptor Potentials
-"""""""""""""""""""""""""""""""""""""""""""""""""
-
-We recall two-body PODs :math:`D^{(2)}_{ik}, 1 \le k \le N_{\rm d}^{(2)}`,
-and three-body PODs :math:`D^{(3)}_{im}, 1 \le m \le N_{\rm d}^{(3)}`,
-with :math:`N_{\rm d}^{(2)} = N_{\rm 2b} N_{\rm e} (N_{\rm e}+1)/2` and
-:math:`N_{\rm d}^{(3)} = N_{\rm 3b} N_{\rm e}^2 (N_{\rm e}+1)/2` being
-the number of descriptors per atom for the two-body PODs and three-body PODs,
-respectively. We employ them to define a new set of atomic descriptors as follows
-
-.. math::
-
-    D^{(2*3)}_{ikm} = \frac{1}{2N}\left( D^{(2)}_{ik} \sum_{j=1}^N D^{(3)}_{jm} + D^{(3)}_{im} \sum_{j=1}^N D^{(2)}_{jk}  \right)
-
-for :math:`1 \le i \le N, 1 \le k \le N_{\rm d}^{(2)}, 1 \le m \le N_{\rm d}^{(3)}`.
-The new descriptors are four-body because they involve central atom :math:`i` together
-with three neighbors :math:`j, k` and :math:`l`. The total number of new  descriptors per atom is equal to
-
-.. math::
-
-    N_{\rm d}^{(2*3)} = N_{\rm d}^{(2)} * N_{\rm d}^{(3)} = N_{\rm 2b} N_{\rm 3b} N_{\rm e}^3 (N_{\rm e}+1)^2/4 .
-
-The new global descriptors are calculated as
-
-.. math::
-
-    d^{(2*3)}_{km} = \sum_{i=1}^N D^{(2*3)}_{ikm} = \left( \sum_{i=1}^N D^{(2)}_{ik} \right) \left( \sum_{i=1}^N D^{(3)}_{im} \right) = d^{(2)}_{k} d^{(3)}_m,
-
-for :math:`1 \le k \le N_{\rm d}^{(2)}, 1 \le m \le N_{\rm d}^{(3)}`. Hence, the gradient
-of the new global descriptors with respect to atom positions is calculated as
-
-.. math::
-
-    \nabla d^{(2*3)}_{km} = d^{(3)}_m \nabla d^{(2)}_{k}  +  d^{(2)}_{k} \nabla d^{(3)}_m, \quad 1 \le k \le N_{\rm d}^{(2)}, 1 \le m \le N_{\rm d}^{(3)} .
-
-The quadratic  POD potential is defined as a linear combination of the
-original and new global descriptors as follows
-
-.. math::
-
-    E^{(2*3)} = \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} \sum_{m=1}^{N_{\rm 3d}^{(2*3)}} c^{(2*3)}_{km} d^{(2*3)}_{km} .
-
-It thus follows that
-
-.. math::
-
-    E^{(2*3)} = 0.5 \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} \left( \sum_{m=1}^{N_{\rm 3d}^{(2*3)}} c^{(2*3)}_{km} d_m^{(3)} \right) d_k^{(2)} + 0.5 \sum_{m=1}^{N_{\rm 3d}^{(2*3)}} \left( \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} c^{(2*3)}_{km} d_k^{(2)} \right) d_m^{(3)}  ,
-
-which is simplified to
-
-.. math::
-
-    E^{(2*3)} =  0.5 \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} b_k^{(2)} d_k^{(2)} +  0.5 \sum_{m=1}^{N_{\rm 3d}^{(2*3)}}   b_m^{(3)} d_m^{(3)}
-
-where
-
-.. math::
-
-    b_k^{(2)} & = \sum_{m=1}^{N_{\rm 3d}^{(2*3)}} c^{(2*3)}_{km} d_m^{(3)}, \quad k = 1,\ldots, N_{\rm 2d}^{(2*3)}, \\
-    b_m^{(3)} & = \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} c^{(2*3)}_{km} d_k^{(2)}, \quad m = 1,\ldots, N_{\rm 3d}^{(2*3)} .
-
-The quadratic POD potential results in the following atomic forces
-
-.. math::
-
-    \boldsymbol F^{(2*3)} = - \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} \sum_{m=1}^{N_{\rm 3d}^{(2*3)}} c^{(2*3)}_{km}  \nabla d^{(2*3)}_{km} .
-
-It can be shown that
-
-.. math::
-
-    \boldsymbol F^{(2*3)} = - \sum_{k=1}^{N_{\rm 2d}^{(2*3)}}   b^{(2)}_k \nabla d_k^{(2)} - \sum_{m=1}^{N_{\rm 3d}^{(2*3)}}  b^{(3)}_m  \nabla d_m^{(3)} .
-
-The calculation of the atomic forces for the quadratic POD  potential
-only requires the extra calculation of :math:`b_k^{(2)}` and :math:`b_m^{(3)}` which can be negligible.
-As a result, the quadratic  POD potential does not increase the computational complexity.
-
+    E_i(\boldsymbol R_i, \boldsymbol Z_i) \ = \ \sum_{m=1}^M c_m \mathcal{D}_{im}(\boldsymbol R_i, \boldsymbol Z_i)
+
+
+Here :math:`c_m` are trainable coefficients and
+:math:`\mathcal{D}_{im}(\boldsymbol R_i, \boldsymbol Z_i)` are per-atom
+POD descriptors. Summing the per-atom descriptors over :math:`i` yields
+the global descriptors :math:`d_m(\boldsymbol R, \boldsymbol Z) =
+\sum_{i=1}^N \mathcal{D}_{im}(\boldsymbol R_i, \boldsymbol Z_i)`.  It
+thus follows that :math:`E(\boldsymbol R, \boldsymbol Z) = \sum_{m=1}^M
+c_m d_m(\boldsymbol R, \boldsymbol Z)`.
+
+The per-atom POD descriptors include one, two, three, four, five, six,
+and seven-body descriptors, which can be specified in the first input
+file. Furthermore, the per-atom POD descriptors also depend on the
+number of environment clusters specified in the first input file.
+Please see :ref:`(Nguyen2024) <Nguyen20242a>` and :ref:`(Nguyen and Sema)
+<Nguyen20243a>` for the detailed description of the per-atom POD
+descriptors.
 
 Training
 """"""""
 
-POD potentials are trained using the least-squares regression against
+A POD potential is trained using the least-squares regression against
 density functional theory (DFT) data.  Let :math:`J` be the number of
 training configurations, with :math:`N_j` being the number of atoms in
-the j-th configuration. Let :math:`\{E^{\star}_j\}_{j=1}^{J}` and
-:math:`\{\boldsymbol F^{\star}_j\}_{j=1}^{J}` be the DFT energies and
-forces for :math:`J` configurations. Next, we calculate the global
-descriptors and their derivatives for all training configurations. Let
-:math:`d_{jm}, 1 \le m \le M`, be the global descriptors associated with
-the j-th configuration, where :math:`M` is the number of global
-descriptors. We then form a matrix :math:`\boldsymbol A \in
-\mathbb{R}^{J \times M}` with entries :math:`A_{jm} = d_{jm}/ N_j` for
-:math:`j=1,\ldots,J` and :math:`m=1,\ldots,M`.  Moreover, we form a
-matrix :math:`\boldsymbol B \in \mathbb{R}^{\mathcal{N} \times M}` by
-stacking the derivatives of the global descriptors for all training
-configurations from top to bottom, where :math:`\mathcal{N} =
-3\sum_{j=1}^{J} N_j`.
+the j-th configuration. The training configurations are extracted from
+the extended XYZ files located in a directory (i.e.,
+path_to_training_data_set in the second input file).  Let
+:math:`\{E^{\star}_j\}_{j=1}^{J}` and :math:`\{\boldsymbol
+F^{\star}_j\}_{j=1}^{J}` be the DFT energies and forces for :math:`J`
+configurations. Next, we calculate the global descriptors and their
+derivatives for all training configurations. Let :math:`d_{jm}, 1 \le m
+\le M`, be the global descriptors associated with the j-th
+configuration, where :math:`M` is the number of global descriptors. We
+then form a matrix :math:`\boldsymbol A \in \mathbb{R}^{J \times M}`
+with entries :math:`A_{jm} = d_{jm}/ N_j` for :math:`j=1,\ldots,J` and
+:math:`m=1,\ldots,M`.  Moreover, we form a matrix :math:`\boldsymbol B
+\in \mathbb{R}^{\mathcal{N} \times M}` by stacking the derivatives of
+the global descriptors for all training configurations from top to
+bottom, where :math:`\mathcal{N} = 3\sum_{j=1}^{J} N_j`.
 
 The coefficient vector :math:`\boldsymbol c` of the POD potential is
 found by solving the following least-squares problem
 
 .. math::
 
-    {\min}_{\boldsymbol c \in \mathbb{R}^{M}} \ w_E \|\boldsymbol A(\boldsymbol \eta) \boldsymbol c - \bar{\boldsymbol E}^{\star} \|^2 + w_F \|\boldsymbol B(\boldsymbol \eta) \boldsymbol c + \boldsymbol F^{\star} \|^2 + w_R \|\boldsymbol c \|^2,
+    {\min}_{\boldsymbol c \in \mathbb{R}^{M}} \ w_E \|\boldsymbol A \boldsymbol c - \bar{\boldsymbol E}^{\star} \|^2 + w_F \|\boldsymbol B \boldsymbol c + \boldsymbol F^{\star} \|^2 + w_R \|\boldsymbol c \|^2,
 
 where :math:`w_E` and :math:`w_F` are weights for the energy
 (*fitting_weight_energy*) and force (*fitting_weight_force*),
-respectively; and :math:`w_R` is the regularization parameter (*fitting_regularization_parameter*).  Here :math:`\bar{\boldsymbol E}^{\star} \in
-\mathbb{R}^{J}` is a vector of with entries :math:`\bar{E}^{\star}_j =
-E^{\star}_j/N_j` and :math:`\boldsymbol F^{\star}` is a vector of
-:math:`\mathcal{N}` entries obtained by stacking :math:`\{\boldsymbol
-F^{\star}_j\}_{j=1}^{J}` from top to bottom.
+respectively; and :math:`w_R` is the regularization parameter
+(*fitting_regularization_parameter*).  Here :math:`\bar{\boldsymbol
+E}^{\star} \in \mathbb{R}^{J}` is a vector of with entries
+:math:`\bar{E}^{\star}_j = E^{\star}_j/N_j` and :math:`\boldsymbol
+F^{\star}` is a vector of :math:`\mathcal{N}` entries obtained by
+stacking :math:`\{\boldsymbol F^{\star}_j\}_{j=1}^{J}` from top to
+bottom.
 
-The training procedure is the same for both the linear and quadratic POD
-potentials.  However, since the quadratic POD potential has a
-significantly larger number of the global descriptors, it is more
-expensive to train the linear POD potential. This is because the
-training of the quadratic POD potential still requires us to calculate
-and store the quadratic global descriptors and their
-gradient. Furthermore, the quadratic POD potential may require more
-training data in order to prevent over-fitting. In order to reduce the
-computational cost of fitting the quadratic POD potential and avoid
-over-fitting, we can use subsets of two-body and three-body PODs for
-constructing the new descriptors.
+Validation
+""""""""""
 
+POD potential can be validated on a test dataset in a directory
+specified by setting path_to_test_data_set in the second input file.  It
+is possible to validate the POD potential after the training is
+complete.  This is done by providing the coefficient file as an input to
+:doc:`fitpod <fitpod_command>`, for example,
+
+.. code-block:: LAMMPS
+
+   fitpod Ta_param.pod Ta_data.pod Ta_coefficients.pod
 
 Restrictions
 """"""""""""
@@ -727,7 +356,11 @@ LAMMPS was built with that package. See the :doc:`Build package
 Related commands
 """"""""""""""""
 
-:doc:`pair_style pod <pair_pod>`
+:doc:`pair_style pod <pair_pod>`,
+:doc:`compute pod/atom <compute_pod_atom>`,
+:doc:`compute podd/atom <compute_pod_atom>`,
+:doc:`compute pod/local <compute_pod_atom>`,
+:doc:`compute pod/global <compute_pod_atom>`
 
 Default
 """""""
@@ -736,10 +369,20 @@ The keyword defaults are also given in the description of the input files.
 
 ----------
 
-.. _Grepl20072:
+.. _Nguyen20222a:
 
-**(Grepl)** Grepl, Maday, Nguyen, and Patera, ESAIM: Mathematical Modelling and Numerical Analysis 41(3), 575-605, (2007).
+**(Nguyen and Rohskopf)** Nguyen and Rohskopf,  Journal of Computational Physics, 480, 112030, (2023).
+
+.. _Nguyen20232a:
+
+**(Nguyen2023)** Nguyen, Physical Review B, 107(14), 144103, (2023).
+
+.. _Nguyen20242a:
+
+**(Nguyen2024)** Nguyen, Journal of Computational Physics, 113102, (2024).
+
+.. _Nguyen20243a:
+
+**(Nguyen and Sema)** Nguyen and Sema, https://arxiv.org/abs/2405.00306, (2024).
 
-.. _Nguyen20222:
 
-**(Nguyen)** Nguyen and Rohskopf, arXiv preprint arXiv:2209.02362 (2022).

doc/src/fix_adapt.rst

@@ -21,17 +21,17 @@ Syntax
        *pair* args = pstyle pparam I J v_name
          pstyle = pair style name (e.g., lj/cut)
          pparam = parameter to adapt over time
-         I,J = type pair(s) to set parameter for
+         I,J = type pair(s) to set parameter for (integer or type label)
          v_name = variable with name that calculates value of pparam
        *bond* args = bstyle bparam I v_name
          bstyle = bond style name (e.g., harmonic)
          bparam = parameter to adapt over time
-         I = type bond to set parameter for
+         I = type bond to set parameter for (integer or type label)
          v_name = variable with name that calculates value of bparam
        *angle* args = astyle aparam I v_name
          astyle = angle style name (e.g., harmonic)
          aparam = parameter to adapt over time
-         I = type angle to set parameter for
+         I = type angle to set parameter for (integer or type label)
          v_name = variable with name that calculates value of aparam
        *kspace* arg = v_name
          v_name = variable with name that calculates scale factor on :math:`k`-space terms
@@ -67,6 +67,9 @@ Examples
    variable ramp_up equal "ramp(0.01,0.5)"
    fix stretch all adapt 1 bond harmonic r0 1 v_ramp_up
 
+   labelmap atom 1 c1
+   fix 1 all adapt 1 pair soft a c1 c1 v_prefactor
+
 Description
 """""""""""
 
@@ -254,10 +257,12 @@ should be specified to indicate which type pairs to apply it to.  If a global
 parameter is specified, the :math:`I` and :math:`J` settings still need to be
 specified, but are ignored.
 
-Similar to the :doc:`pair_coeff command <pair_coeff>`, :math:`I` and :math:`J`
-can be specified in one of two ways.  Explicit numeric values can be used for
-each, as in the first example above.  :math:`I \le J` is required.  LAMMPS sets
-the coefficients for the symmetric :math:`J,I` interaction to the same values.
+Similar to the :doc:`pair_coeff command <pair_coeff>`, :math:`I` and
+:math:`J` can be specified in one of several ways.  Explicit numeric values
+can be used for each, as in the first example above.  Or, one or both of
+the types in the I,J pair can be a :doc:`type label <Howto_type_labels>`.
+LAMMPS sets the coefficients for the symmetric :math:`J,I` interaction to
+the same values.
 
 A wild-card asterisk can be used in place of or in conjunction with
 the :math:`I,J` arguments to set the coefficients for multiple pairs of atom
@@ -266,8 +271,9 @@ is the number of atom types, then an asterisk with no numeric values
 means all types from 1 to :math:`N`.  A leading asterisk means all types from
 1 to n (inclusive).  A trailing asterisk means all types from m to :math:`N`
 (inclusive).  A middle asterisk means all types from m to n
-(inclusive).  Note that only type pairs with :math:`I \le J` are considered; if
-asterisks imply type pairs where :math:`J < I`, they are ignored.
+(inclusive).  For the asterisk syntax, note that only type pairs with
+:math:`I \le J` are considered; if asterisks imply type pairs where
+:math:`J < I`, they are ignored.
 
 IMPORTANT NOTE: If :doc:`pair_style hybrid or hybrid/overlay
 <pair_hybrid>` is being used, then the *pstyle* will be a sub-style

doc/src/fix_adapt_fep.rst

@@ -21,13 +21,13 @@ Syntax
        *pair* args = pstyle pparam I J v_name
          pstyle = pair style name (e.g., lj/cut)
          pparam = parameter to adapt over time
-         I,J = type pair(s) to set parameter for
+         I,J = type pair(s) to set parameter for (integer or type label)
          v_name = variable with name that calculates value of pparam
        *kspace* arg = v_name
          v_name = variable with name that calculates scale factor on K-space terms
        *atom* args = aparam v_name
          aparam = parameter to adapt over time
-         I = type(s) to set parameter for
+         I = type(s) to set parameter for (integer or type label)
          v_name = variable with name that calculates value of aparam
 
 * zero or more keyword/value pairs may be appended
@@ -56,6 +56,9 @@ Examples
    fix 1 all adapt/fep 1 pair lj/cut epsilon * * v_scale1 coul/cut scale 3 3 v_scale2 scale yes reset yes
    fix 1 all adapt/fep 10 atom diameter 1 v_size
 
+   labelmap atom 1 c1
+   fix 1 all adapt/fep 1 pair soft a c1 c1 v_prefactor
+
 
 Example input scripts available: examples/PACKAGES/fep
 
@@ -218,10 +221,12 @@ be specified to indicate which type pairs to apply it to.  If a global
 parameter is specified, the *I* and *J* settings still need to be
 specified, but are ignored.
 
-Similar to the :doc:`pair_coeff command <pair_coeff>`, I and J can be
-specified in one of two ways.  Explicit numeric values can be used for
-each, as in the first example above.  :math:`I \le J` is required.  LAMMPS sets
-the coefficients for the symmetric J,I interaction to the same values.
+Similar to the :doc:`pair_coeff command <pair_coeff>`, :math:`I` and
+:math:`J` can be specified in one of several ways.  Explicit numeric values
+can be used for each, as in the first example above.  Or, one or both of
+the types in the I,J pair can be a :doc:`type label <Howto_type_labels>`.
+LAMMPS sets the coefficients for the symmetric :math:`J,I` interaction to
+the same values.
 
 A wild-card asterisk can be used in place of or in conjunction with
 the :math:`I,J` arguments to set the coefficients for multiple pairs of atom
@@ -230,8 +235,9 @@ the number of atom types, then an asterisk with no numeric values means
 all types from 1 to :math:`N`.  A leading asterisk means all types from 1 to n
 (inclusive).  A trailing asterisk means all types from m to :math:`N`
 (inclusive).  A middle asterisk means all types from m to n
-(inclusive).  Note that only type pairs with :math:`I \le J` are considered; if
-asterisks imply type pairs where :math:`J < I`, they are ignored.
+(inclusive).  For the asterisk syntax, note that only type pairs with
+:math:`I \le J` are considered; if asterisks imply type pairs where
+:math:`J < I`, they are ignored.
 
 IMPROTANT NOTE: If :doc:`pair_style hybrid or hybrid/overlay <pair_hybrid>` is
 being used, then the *pstyle* will be a sub-style name.  You must specify

doc/src/fix_amoeba_bitorsion.rst

@@ -35,7 +35,11 @@ the implementation of AMOEBA and HIPPO in LAMMPS.
 
 Bitorsion interactions add additional potential energy contributions
 to pairs of overlapping phi-psi dihedrals of amino-acids, which are
-important to properly represent their conformational behavior.
+important to properly represent their conformational behavior.  Each
+bitorsion interaction is thus defined for a 5-tuple of atoms
+:math:`IJKLM` with bonds between successive atoms in the list,
+i.e. two overlapping dihedral interactions for atoms :math:`IJKL` and
+:math:`JKLM`.
 
 The examples/amoeba directory has a sample input script and data file
 for ubiquitin, which illustrates use of the fix amoeba/bitorsion
@@ -68,14 +72,15 @@ lines:
           [...]
           N       3     314     315     317      318    330
 
-The first column is an index from 1 to :math:`N` to enumerate the bitorsion
-5-atom tuples; it is ignored by LAMMPS.  The second column is the
-*type* of the interaction; it is an index into the bitorsion force
-field file.  The remaining 5 columns are the atom IDs of the atoms in
-the two 4-atom dihedrals that overlap to create the bitorsion 5-body
-interaction.  Note that the *bitorsions* and *BiTorsions* keywords for
-the header and body sections match those specified in the
-:doc:`read_data <read_data>` command following the data file name.
+The first column is an index from 1 to :math:`N` to enumerate the
+bitorsion 5-atom tuples; it is ignored by LAMMPS.  The second column
+is the *type* of the interaction; it is an index into the bitorsion
+force field file.  The remaining 5 columns are the atom IDs of the
+atoms (in order) for the 5-tuple :math:`IJKLM`, as described above.
+
+Note that the *bitorsions* and *BiTorsions* keywords for the header
+and body sections match those specified in the :doc:`read_data
+<read_data>` command following the data file name.
 
 The data file should be generated by using the
 tools/tinker/tinker2lmp.py conversion script which creates a LAMMPS

doc/src/fix_amoeba_pitorsion.rst

@@ -57,7 +57,7 @@ should have two lines like these in its header section:
    M pitorsion types
    N pitorsions
 
-where :math:`N` is the number of pitorsion 5-body interactions and :math:`M` is
+where :math:`N` is the number of pitorsion 6-body interactions and :math:`M` is
 the number of pitorsion types.  It should also have two sections in the body
 of the data file like these with :math:`M` and :math:`N` lines each:
 
@@ -74,21 +74,20 @@ of the data file like these with :math:`M` and :math:`N` lines each:
 
    PiTorsions
 
-          1       1       8      10      12      18      20
-          2       5      18      20      22      25      27
+          1       1       2      4      3      20      21      24
+          2       5      21     23     22      37      38      41
           [...]
-          N       3     314     315     317      318    330
+          N       7      27     29     28      30      35      36
 
-For PiTorsion Coeffs, the first column is an index from 1 to :math:`M` to
-enumerate the pitorsion types.  The second column is the single
+For PiTorsion Coeffs, the first column is an index from 1 to :math:`M`
+to enumerate the pitorsion types.  The second column is the single
 prefactor coefficient needed for each type.
 
-For PiTorsions, the first column is an index from 1 to :math:`N` to enumerate
-the pitorsion 5-atom tuples; it is ignored by LAMMPS.  The second
-column is the "type" of the interaction; it is an index into the
-PiTorsion Coeffs.  The remaining 5 columns are the atom IDs of the
-atoms in the two 4-atom dihedrals that overlap to create the pitorsion
-5-body interaction.
+For PiTorsions, the first column is an index from 1 to :math:`N` to
+enumerate the pitorsion 6-atom tuples; it is ignored by LAMMPS.  The
+second column is the "type" of the interaction; it is an index into
+the PiTorsion Coeffs.  The remaining 6 columns are the atom IDs of the
+atoms (in order) for the 6-tuple :math:`IJKLMN`, as described above.
 
 Note that the *pitorsion types* and *pitorsions* and *PiTorsion
 Coeffs* and *PiTorsions* keywords for the header and body sections of

doc/src/fix_atom_swap.rst

@@ -21,7 +21,7 @@ Syntax
 
   .. parsed-literal::
 
-       *types* values = two or more atom types
+       *types* values = two or more atom types (1-Ntypes or type label)
        *mu* values = chemical potential of swap types (energy units)
        *ke* value = *no* or *yes*
          *no* = no conservation of kinetic energy after atom swaps
@@ -168,7 +168,7 @@ the following global cumulative quantities:
 * 1 = swap attempts
 * 2 = swap accepts
 
-The vector values calculated by this fix are "extensive".
+The vector values calculated by this fix are "intensive".
 
 No parameter of this fix can be used with the *start/stop* keywords of
 the :doc:`run <run>` command.  This fix is not invoked during

doc/src/fix_bond_break.rst

@@ -13,7 +13,7 @@ Syntax
 * ID, group-ID are documented in :doc:`fix <fix>` command
 * bond/break = style name of this fix command
 * Nevery = attempt bond breaking every this many steps
-* bondtype = type of bonds to break
+* bondtype = type of bonds to break (integer or type label)
 * Rmax = bond longer than Rmax can break (distance units)
 * zero or more keyword/value pairs may be appended
 * keyword = *prob*

doc/src/fix_bond_create.rst

@@ -17,9 +17,9 @@ Syntax
 * ID, group-ID are documented in :doc:`fix <fix>` command
 * bond/create = style name of this fix command
 * Nevery = attempt bond creation every this many steps
-* itype,jtype = atoms of itype can bond to atoms of jtype
+* itype,jtype = atoms of itype can bond to atoms of jtype (1-Ntypes or type label)
 * Rmin = 2 atoms separated by less than Rmin can bond (distance units)
-* bondtype = type of created bonds
+* bondtype = type of created bonds (integer or type label)
 * zero or more keyword/value pairs may be appended to args
 * keyword = *iparam* or *jparam* or *prob* or *atype* or *dtype* or *itype* or *aconstrain*
 
@@ -27,19 +27,19 @@ Syntax
 
        *iparam* values = maxbond, newtype
          maxbond = max # of bonds of bondtype the itype atom can have
-         newtype = change the itype atom to this type when maxbonds exist
+         newtype = change the itype atom to this type when maxbonds exist (1-Ntypes or type label)
        *jparam* values = maxbond, newtype
          maxbond = max # of bonds of bondtype the jtype atom can have
-         newtype = change the jtype atom to this type when maxbonds exist
+         newtype = change the jtype atom to this type when maxbonds exist (1-Ntypes or type label)
        *prob* values = fraction seed
          fraction = create a bond with this probability if otherwise eligible
          seed = random number seed (positive integer)
        *atype* value = angletype
-         angletype = type of created angles
+         angletype = type of created angles (integer or type label)
        *dtype* value = dihedraltype
-         dihedraltype = type of created dihedrals
+         dihedraltype = type of created dihedrals (integer or type label)
        *itype* value = impropertype
-         impropertype = type of created impropers
+         impropertype = type of created impropers (integer or type label)
        *aconstrain* value = amin amax
          amin = minimal angle at which new bonds can be created
          amax = maximal angle at which new bonds can be created
@@ -54,6 +54,10 @@ Examples
    fix 5 all bond/create 1 3 3 0.8 1 prob 0.5 85784 iparam 2 3 atype 1 dtype 2
    fix 5 all bond/create/angle 10 1 2 1.122 1 aconstrain 120 180 prob 1 4928459 iparam 2 1 jparam 2 2
 
+   labelmap atom 1 c1 2 n2
+   labelmap bond 1 c1-n2
+   fix 5 all bond/create 10 c1 n2 0.8 c1-n2
+
 Description
 """""""""""
 

doc/src/fix_charge_regulation.rst

@@ -13,8 +13,8 @@ Syntax
 
 * ID, group-ID are documented in fix command
 * charge/regulation = style name of this fix command
-* cation_type = atom type of free cations
-* anion_type = atom type of free anions
+* cation_type = atom type of free cations (integer or type label)
+* anion_type = atom type of free anions (integer or type label)
 
 * zero or more keyword/value pairs may be appended
 
@@ -27,8 +27,8 @@ Syntax
      *pIp* value = activity (effective concentration) of free cations (in the -log10 representation)
      *pIm* value = activity (effective concentration) of free anions (in the -log10 representation)
      *pKs* value = solvent self-dissociation constant (in the -log10 representation)
-     *acid_type* = atom type of acid groups
-     *base_type*  = atom type of base groups
+     *acid_type* = atom type of acid groups (integer or type label)
+     *base_type*  = atom type of base groups (integer or type label)
      *lunit_nm* value = unit length used by LAMMPS (# in the units of nanometers)
      *temp* value = temperature
      *tempfixid* value = fix ID of temperature thermostat
@@ -51,6 +51,9 @@ Examples
     fix chareg all charge/regulation 1 2 acid_type 3 base_type 4 pKa 5.0 pKb 6.0 pH 7.0 pIp 3.0 pIm 3.0 nevery 200 nmc 200 seed 123 tempfixid fT
     fix chareg all charge/regulation 1 2 pIp 3 pIm 3 onlysalt yes 2 -1 seed 123 tag yes temp 1.0
 
+    labelmap atom 1 H+ 2 OH-
+    fix chareg all charge/regulation H+ OH- pIp 3 pIm 3 onlysalt yes 2 -1 seed 123 tag yes temp 1.0
+
 Description
 """""""""""
 

doc/src/fix_deform.rst

@@ -64,6 +64,8 @@ Syntax
                  effectively an engineering shear strain rate
            *erate* value = R
              R = engineering shear strain rate (1/time units)
+           *erate/rescale* value = R (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
+             R = engineering shear strain rate (1/time units)
            *trate* value = R
              R = true shear strain rate (1/time units)
            *wiggle* values = A Tp

doc/src/fix_deform_pressure.rst

@@ -29,10 +29,12 @@ Syntax
            NOTE: All other styles are documented by the :doc:`fix deform <fix_deform>` command
 
        *xy*, *xz*, *yz* args = style value
-         style = *final* or *delta* or *vel* or *erate* or *trate* or *wiggle* or *variable* or *pressure*
+         style = *final* or *delta* or *vel* or *erate* or *trate* or *wiggle* or *variable* or *pressure* or *erate/rescale*
            *pressure* values = target gain
              target = target pressure (pressure units)
              gain = proportional gain constant (1/(time * pressure) or 1/time units)
+           *erate/rescale* value = R
+             R = engineering shear strain rate (1/time units)
            NOTE: All other styles are documented by the :doc:`fix deform <fix_deform>` command
 
        *box* = style value
@@ -159,6 +161,21 @@ details of a simulation and testing different values is
 recommended. One can also apply a maximum limit to the magnitude of
 the applied strain using the :ref:`max/rate <deform_max_rate>` option.
 
+The *erate/rescale* style operates similarly to the *erate* style with
+a specified strain rate in units of 1/time. The difference is that
+the change in the tilt factor will depend on the current length of
+the box perpendicular to the shear direction, L, instead of the
+original length, L0. The tilt factor T as a function of time will
+change as
+
+.. parsed-literal::
+
+   T(t) = T(t-1) + L\*erate\* \Delta t
+
+where T(t-1) is the tilt factor on the previous timestep and :math:`\Delta t`
+is the timestep size. This option may be useful in scenarios where
+L changes in time.
+
 ----------
 
 The *box* parameter provides an additional control over the *x*, *y*,
@@ -294,6 +311,10 @@ This fix is not invoked during :doc:`energy minimization <minimize>`.
 Restrictions
 """"""""""""
 
+This fix is part of the EXTRA-FIX package.  It is only enabled if LAMMPS
+was built with that package.  See the :doc:`Build package <Build_package>`
+page for more info.
+
 You cannot apply x, y, or z deformations to a dimension that is
 shrink-wrapped via the :doc:`boundary <boundary>` command.
 

doc/src/fix_deposit.rst

@@ -13,7 +13,7 @@ Syntax
 * ID, group-ID are documented in :doc:`fix <fix>` command
 * deposit = style name of this fix command
 * N = # of atoms or molecules to insert
-* type = atom type to assign to inserted atoms (offset for molecule insertion)
+* type = atom type (1-Ntypes or type label) to assign to inserted atoms (offset for molecule insertion)
 * M = insert a single atom or molecule every M steps
 * seed = random # seed (positive integer)
 * one or more keyword/value pairs may be appended to args
@@ -76,6 +76,9 @@ Examples
    fix 4 sputter deposit 1000 2 500 12235 region sphere vz -1.0 -1.0 target 5.0 5.0 0.0 units lattice
    fix 5 insert deposit 200 2 100 777 region disk gaussian 5.0 5.0 9.0 1.0 units box
 
+   labelmap atom 1 Au
+   fix 4 sputter deposit 1000 Au 500 12235 region sphere vz -1.0 -1.0 target 5.0 5.0 0.0 units lattice
+
 Description
 """""""""""
 

doc/src/fix_gcmc.rst

@@ -15,7 +15,7 @@ Syntax
 * N = invoke this fix every N steps
 * X = average number of GCMC exchanges to attempt every N steps
 * M = average number of MC moves to attempt every N steps
-* type = atom type for inserted atoms (must be 0 if mol keyword used)
+* type = atom type (1-Ntypes or type label) for inserted atoms (must be 0 if mol keyword used)
 * seed = random # seed (positive integer)
 * T = temperature of the ideal gas reservoir (temperature units)
 * mu = chemical potential of the ideal gas reservoir (energy units)
@@ -45,7 +45,7 @@ Syntax
        *group* value = group-ID
          group-ID = group-ID for inserted atoms (string)
        *grouptype* values = type group-ID
-         type = atom type (int)
+         type = atom type (1-Ntypes or type label)
          group-ID = group-ID for inserted atoms (string)
        *intra_energy* value = intramolecular energy (energy units)
        *tfac_insert* value = scale up/down temperature of inserted atoms (unitless)
@@ -62,52 +62,47 @@ Examples
    fix 3 water gcmc 10 100 100 0 3456543 3.0 -2.5 0.1 mol my_one_water maxangle 180 full_energy
    fix 4 my_gas gcmc 1 10 10 1 123456543 300.0 -12.5 1.0 region disk
 
+   labelmap atom 1 Li
+   fix 2 ion gcmc 10 1000 1000 Li 29494 298.0 -0.5 0.01
+
 Description
 """""""""""
 
-This fix performs grand canonical Monte Carlo (GCMC) exchanges of
-atoms or molecules with an imaginary ideal gas
-reservoir at the specified T and chemical potential (mu) as discussed
-in :ref:`(Frenkel) <Frenkel2>`. It also
-attempts  Monte Carlo (MC) moves (translations and molecule
-rotations) within the simulation cell or
-region. If used with the :doc:`fix nvt <fix_nh>`
+This fix performs grand canonical Monte Carlo (GCMC) exchanges of atoms or
+molecules with an imaginary ideal gas reservoir at the specified T and
+chemical potential (mu) as discussed in :ref:`(Frenkel) <Frenkel2>`.  It
+also attempts Monte Carlo (MC) moves (translations and molecule rotations)
+within the simulation cell or region.  If used with the :doc:`fix nvt <fix_nh>`
 command, simulations in the grand canonical ensemble (muVT, constant
 chemical potential, constant volume, and constant temperature) can be
 performed.  Specific uses include computing isotherms in microporous
 materials, or computing vapor-liquid coexistence curves.
 
-Every N timesteps the fix attempts both GCMC exchanges
-(insertions or deletions) and MC moves of gas atoms or molecules.
-On those timesteps, the average number of attempted GCMC exchanges is X,
-while the average number of attempted MC moves is M.
-For GCMC exchanges of either molecular or atomic gasses,
-these exchanges can be either deletions or insertions,
-with equal probability.
+Every N timesteps the fix attempts both GCMC exchanges (insertions or
+deletions) and MC moves of gas atoms or molecules.  On those timesteps, the
+average number of attempted GCMC exchanges is X, while the average number
+of attempted MC moves is M.  For GCMC exchanges of either molecular or
+atomic gasses, these exchanges can be either deletions or insertions, with
+equal probability.
 
-The possible choices for MC moves are translation of an atom,
-translation of a molecule, and rotation of a molecule.
-The relative amounts of each are determined by the optional
-*mcmoves* keyword (see below).
-The default behavior is as follows.
-If the *mol* keyword is used, only molecule translations
-and molecule rotations are performed with equal probability.
-Conversely, if the *mol* keyword is not used, only atom
-translations are performed.
-M should typically be
-chosen to be approximately equal to the expected number of gas atoms
-or molecules of the given type within the simulation cell or region,
-which will result in roughly one MC move per atom or molecule
-per MC cycle.
+The possible choices for MC moves are translation of an atom, translation
+of a molecule, and rotation of a molecule.  The relative amounts of each are
+determined by the optional *mcmoves* keyword (see below).  The default
+behavior is as follows. If the *mol* keyword is used, only molecule
+translations and molecule rotations are performed with equal probability.
+Conversely, if the *mol* keyword is not used, only atom translations are
+performed.  M should typically be chosen to be approximately equal to the
+expected number of gas atoms or molecules of the given type within the
+simulation cell or region, which will result in roughly one MC move per
+atom or molecule per MC cycle.
 
-All inserted particles are always added to two groups: the default
-group "all" and the fix group specified in the fix command.
-In addition, particles are also added to any groups
-specified by the *group* and *grouptype* keywords.  If inserted
-particles are individual atoms, they are assigned the atom type given
-by the type argument.  If they are molecules, the type argument has no
-effect and must be set to zero. Instead, the type of each atom in the
-inserted molecule is specified in the file read by the
+All inserted particles are always added to two groups: the default group
+"all" and the fix group specified in the fix command.  In addition,
+particles are also added to any groups specified by the *group* and
+*grouptype* keywords.  If inserted particles are individual atoms, they are
+assigned the atom type given by the type argument.  If they are molecules,
+the type argument has no effect and must be set to zero. Instead, the type
+of each atom in the inserted molecule is specified in the file read by the
 :doc:`molecule <molecule>` command.
 
 .. note::
@@ -427,7 +422,7 @@ the following global cumulative quantities:
 * 7 = rotation attempts
 * 8 = rotation successes
 
-The vector values calculated by this fix are "extensive".
+The vector values calculated by this fix are "intensive".
 
 No parameter of this fix can be used with the *start/stop* keywords of
 the :doc:`run <run>` command.  This fix is not invoked during

doc/src/fix_hyper_local.rst

@@ -512,8 +512,7 @@ Value 27 computes the average boost for biased bonds only on this step.
 Value 28 is the count of bonds with an absolute value of strain >= q
 on this step.
 
-The scalar value is an "extensive" quantity since it grows with the
-system size; the vector values are all "intensive".
+The scalar value and vector values are all "intensive".
 
 This fix also computes a local vector of length the number of bonds
 currently in the system.  The value for each bond is its :math:`C_{ij}`

doc/src/fix_ipi.rst

@@ -35,23 +35,24 @@ Description
 """""""""""
 
 This fix enables LAMMPS to be run as a client for the i-PI Python
-wrapper :ref:`(IPI) <ipihome>` for performing a path integral molecular dynamics
-(PIMD) simulation.  The philosophy behind i-PI is described in the
-following publication :ref:`(IPI-CPC) <IPICPC>`.
+wrapper :ref:`(IPI) <ipihome>`. i-PI is a universal force engine,
+designed to perform advanced molecular simulations, with a special
+focus on path integral molecular dynamics (PIMD) simulation.
+The philosophy behind i-PI is to separate the evaluation of the
+energy and forces, which is delegated to the client, and the evolution
+of the dynamics, that is the responsibility of i-PI. This approach also
+simplifies combining energies computed from different codes, which
+can for instance be useful to mix first-principles calculations,
+empirical force fields or machine-learning potentials.
+The following publication :ref:`(IPI-CPC-2014) <IPICPC>` discusses the
+overall implementation of i-PI, and focuses on path-integral techniques,
+while a later release :ref:`(IPI-CPC-2019) <IPICPC2>` introduces several
+additional features and simulation schemes.
 
-A version of the i-PI package, containing only files needed for use
-with LAMMPS, is provided in the tools/i-pi directory.  See the
-tools/i-pi/manual.pdf for an introduction to i-PI.  The
-examples/PACKAGES/i-pi directory contains example scripts for using i-PI
-with LAMMPS.
-
-In brief, the path integral molecular dynamics is performed by the
-Python wrapper, while the client (LAMMPS in this case) simply computes
-forces and energy for each configuration. The communication between
-the two components takes place using sockets, and is reduced to the
-bare minimum. All the parameters of the dynamics are specified in the
-input of i-PI, and all the parameters of the force field must be
-specified as LAMMPS inputs, preceding the *fix ipi* command.
+The communication between i-PI and LAMMPS takes place using sockets,
+and is reduced to the bare minimum. All the parameters of the dynamics
+are specified in the input of i-PI, and all the parameters of the force
+field must be specified as LAMMPS inputs, preceding the *fix ipi* command.
 
 The server address must be specified by the *address* argument, and
 can be either the IP address, the fully-qualified name of the server,
@@ -75,6 +76,20 @@ If the cell varies too wildly, it may be advisable to re-initialize
 these interactions at each call. This behavior can be requested by
 setting the *reset* switch.
 
+Obtaining i-PI
+""""""""""""""
+
+Here are the commands to set up a virtual environment and install
+i-PI into it with all its dependencies via the PyPI repository and
+the pip package manager.
+
+.. code-block:: sh
+
+   python -m venv ipienv
+   source ipienv/bin/activate
+   pip install --upgrade pip
+   pip install ipi
+
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
@@ -111,9 +126,14 @@ Related commands
 
 .. _IPICPC:
 
-**(IPI-CPC)** Ceriotti, More and Manolopoulos, Comp Phys Comm, 185,
+**(IPI-CPC-2014)** Ceriotti, More and Manolopoulos, Comp Phys Comm 185,
 1019-1026 (2014).
 
+.. _IPICPC2:
+
+**(IPI-CPC-2019)** Kapil et al., Comp Phys Comm 236, 214-223 (2019).
+
+
 .. _ipihome:
 
 **(IPI)**

doc/src/fix_mol_swap.rst

@@ -14,7 +14,7 @@ Syntax
 * atom/swap = style name of this fix command
 * N = invoke this fix every N steps
 * X = number of swaps to attempt every N steps
-* itype,jtype = two atom types to swap with each other
+* itype,jtype = two atom types (1-Ntypes or type label) to swap with each other
 * seed = random # seed (positive integer)
 * T = scaling temperature of the MC swaps (temperature units)
 * zero or more keyword/value pairs may be appended to args
@@ -34,6 +34,9 @@ Examples
    fix 2 all mol/swap 100 1 2 3 29494 300.0 ke no
    fix mySwap fluid mol/swap 500 10 1 2 482798 1.0
 
+   labelmap atom 1 A 2 B
+   fix mySwap fluid mol/swap 500 10 A B 482798 1.0
+
 Description
 """""""""""
 
@@ -146,7 +149,7 @@ the following global cumulative quantities:
 * 1 = swap attempts
 * 2 = swap accepts
 
-The vector values calculated by this fix are "extensive".
+The vector values calculated by this fix are "intensive".
 
 No parameter of this fix can be used with the *start/stop* keywords of
 the :doc:`run <run>` command.  This fix is not invoked during

doc/src/fix_nonaffine_displacement.rst

@@ -8,7 +8,7 @@ Syntax
 
 .. parsed-literal::
 
-   fix ID group nonaffine/displacement style args reference/style nstep
+   fix ID group nonaffine/displacement style args reference/style nstep keyword values
 
 * ID, group are documented in :doc:`fix <fix>` command
 * nonaffine/displacement = style name of this fix command
@@ -32,6 +32,13 @@ Syntax
        *update* = update the reference frame every *nstep* timesteps
        *offset* = update the reference frame *nstep* timesteps before calculating the nonaffine displacement
 
+* zero or more keyword/value pairs may be appended
+
+  .. parsed-literal::
+
+       *z/min* values = zmin
+         zmin = minimum coordination number to calculate D2min
+
 Examples
 """"""""
 
@@ -76,6 +83,12 @@ is the identity tensor. This calculation is only performed on timesteps that
 are a multiple of *nevery* (including timestep zero). Data accessed before
 this occurs will simply be zeroed.
 
+For particles with low coordination numbers, calculations of :math:`D^2_\mathrm{min}`
+may not be accurate. An optional minimum coordination number can be defined using
+the *z/min* keyword. If any particle has fewer than the specified number of particles
+in the cutoff distance or in contact, the above calculations will be skipped and the
+corresponding peratom array entries will be zero.
+
 The *integrated* style simply integrates the velocity of particles
 every timestep to calculate a displacement. This style only works if
 used in conjunction with another fix that deforms the box and displaces

doc/src/fix_reaxff_species.rst

@@ -20,7 +20,7 @@ Syntax
 * Nfreq = calculate average bond-order every this many timesteps
 * filename = name of output file
 * zero or more keyword/value pairs may be appended
-* keyword = *cutoff* or *element* or *position* or *delete*
+* keyword = *cutoff* or *element* or *position* or *delete* or *delete_rate_limit*
 
   .. parsed-literal::
 
@@ -110,10 +110,10 @@ all types from 1 to :math:`N`.  A leading asterisk means all types from
 The optional keyword *element* can be used to specify the chemical
 symbol printed for each LAMMPS atom type. The number of symbols must
 match the number of LAMMPS atom types and each symbol must consist of
-1 or 2 alphanumeric characters. Normally, these symbols should be
-chosen to match the chemical identity of each LAMMPS atom type, as
-specified using the :doc:`reaxff pair_coeff <pair_reaxff>` command and
-the ReaxFF force field file.
+1 or 2 alphanumeric characters. By default, these symbols are the same
+as the chemical identity of each LAMMPS atom type, as specified by the
+:doc:`ReaxFF pair_coeff <pair_reaxff>` command and the ReaxFF force
+field file.
 
 The optional keyword *position* writes center-of-mass positions of
 each identified molecules to file *filepos* every *posfreq* timesteps.
@@ -134,36 +134,34 @@ value.  For example, AuO.pos.\* becomes AuO.pos.0, AuO.pos.1000, etc.
 
 .. versionadded:: 3Aug2022
 
-The optional keyword *delete* enables the periodic removal of
-molecules from the system.  Criteria for deletion can be either a list
-of specific chemical formulae or a range of molecular weights.
-Molecules are deleted every *Nfreq* timesteps, and bond connectivity
-is determined using the *Nevery* and *Nrepeat* keywords.  The
-*filedel* argument is the name of the output file that records the
-species that are removed from the system.  The *specieslist* keyword
-permits specific chemical species to be deleted.  The *Nspecies*
-argument specifies how many species are eligible for deletion and is
-followed by a list of chemical formulae, whose strings are compared to
-species identified by this fix.  For example, "specieslist 2 CO CO2"
-deletes molecules that are identified as "CO" and "CO2" in the species
-output file.  When using the *specieslist* keyword, the *filedel* file
-has the following format: the first line lists the chemical formulae
-eligible for deletion, and each additional line contains the timestep
-on which a molecule deletion occurs and the number of each species
-deleted on that timestep.  The *masslimit* keyword permits deletion of
-molecules with molecular weights between *massmin* and *massmax*.
-When using the *masslimit* keyword, each line of the *filedel* file
-contains the timestep on which deletions occurs, followed by how many
-of each species are deleted (with quantities preceding chemical
-formulae).  The *specieslist* and *masslimit* keywords cannot both be
-used in the same *reaxff/species* fix.  The *delete_rate_limit*
-keyword can enforce an upper limit on the overall rate of molecule
-deletion.  The number of deletion occurrences is limited to Nlimit
-within an interval of Nsteps timesteps.   Nlimit can be specified with
-an equal-style :doc:`variable <variable>`.  When using the
-*delete_rate_limit* keyword, no deletions are permitted to occur
-within the first Nsteps timesteps of the first run (after reading a
-either a data or restart file).
+The optional keyword *delete* enables the periodic removal of molecules
+from the system :ref:`(Gissinger) <Delete>`.  Criteria for deletion can
+be either a list of specific chemical formulae or a range of molecular
+weights.  Molecules are deleted every *Nfreq* timesteps, and bond
+connectivity is determined using the *Nevery* and *Nrepeat* keywords.  The
+*filedel* argument is the name of the output file that records the species
+that are removed from the system.  The *specieslist* keyword permits
+specific chemical species to be deleted.  The *Nspecies* argument specifies
+how many species are eligible for deletion and is followed by a list of
+chemical formulae, whose strings are compared to species identified by this
+fix.  For example, "specieslist 2 CO CO2" deletes molecules that are
+identified as "CO" and "CO2" in the species output file.  When using the
+*specieslist* keyword, the *filedel* file has the following format: the
+first line lists the chemical formulae eligible for deletion, and each
+additional line contains the timestep on which a molecule deletion occurs
+and the number of each species deleted on that timestep.  The *masslimit*
+keyword permits deletion of molecules with molecular weights between
+*massmin* and *massmax*.  When using the *masslimit* keyword, each line of
+the *filedel* file contains the timestep on which deletions occurs,
+followed by how many of each species are deleted (with quantities preceding
+chemical formulae).  The *specieslist* and *masslimit* keywords cannot both
+be used in the same *reaxff/species* fix.  The *delete_rate_limit* keyword
+can enforce an upper limit on the overall rate of molecule deletion.  The
+number of deletion occurrences is limited to Nlimit within an interval of
+Nsteps timesteps.  Nlimit can be specified with an equal-style
+:doc:`variable <variable>`.  When using the *delete_rate_limit* keyword, no
+deletions are permitted to occur within the first Nsteps timesteps of the
+first run (after reading a either a data or restart file).
 
 ----------
 
@@ -233,5 +231,9 @@ Default
 """""""
 
 The default values for bond-order cutoffs are 0.3 for all I-J pairs.
-The default element symbols are C, H, O, N.
+The default element symbols are taken from the ReaxFF pair_coeff command.
 Position files are not written by default.
+
+.. _Delete:
+
+**(Gissinger)** Jacob R. Gissinger, Scott R. Zavada, Joseph G. Smith, Josh Kemppainen, Ivan Gallegos, Gregory M. Odegard, Emilie J. Siochi, and Kristopher E. Wise, Carbon, 202, 336-347 (2023).

doc/src/fix_sgcmc.rst

@@ -15,7 +15,7 @@ Syntax
 * every_nsteps = number of MD steps between MC cycles
 * swap_fraction = fraction of a full MC cycle carried out at each call (a value of 1.0 will perform as many trial moves as there are atoms)
 * temperature = temperature that enters Boltzmann factor in Metropolis criterion (usually the same as MD temperature)
-* deltamu = chemical potential difference(s) (`N-1` values must be provided, with `N` being the number of elements)
+* deltamu = `N-1` chemical potential differences :math:`\mu_1-\mu_2, \ldots, \mu_1-\mu_N` (`N` is the number of atom types)
 * Zero or more keyword/value pairs may be appended to fix definition line:
 
   .. parsed-literal::
@@ -23,7 +23,7 @@ Syntax
      keyword = *variance* or *randseed* or *window_moves* or *window_size*
        *variance* kappa conc1 [conc2] ... [concN]
          kappa = variance constraint parameter
-         conc1,conc2,...  = target concentration(s) in the range 0.0-1.0 (*N-1* values must be provided, with *N* being the number of elements)
+         `c_2`, `c_3`,..., `c_N` = `N-1` target concentration fractions
        *randseed* N
          N = seed for pseudo random number generator
        *window_moves* N
@@ -90,11 +90,10 @@ the simulation, e.g., to speed up equilibration at low temperatures.
 
 ------------
 
-The parameter *deltamu* is used to set the chemical potential difference
-in the SGC MC algorithm (see Eq. 16 in :ref:`Sadigh1 <Sadigh1>`). By
-convention it is the difference of the chemical potentials of elements
-`B`, `C` ..., with respect to element A. When the simulation includes
-`N` elements, `N-1` values must be specified.
+The parameter *deltamu* is used to set the chemical potential differences
+in the SGC MC algorithm (see Eq. 16 in :ref:`Sadigh1 <Sadigh1>`).
+The `N-1` differences are defined as :math:`\mu_1-\mu_2, \ldots, \mu_1-\mu_N`,
+where `N` is the number of atom types.
 
 ------------
 
@@ -104,12 +103,12 @@ the effective average constraint in the parallel VC-SGC MC algorithm
 (parameter :math:`\delta\mu_0` in Eq. (20) of :ref:`Sadigh1
 <Sadigh1>`). The parameter *kappa* specifies the variance constraint
 (see Eqs. (20-21) in :ref:`Sadigh1 <Sadigh1>`).
-
-The parameter *conc* sets the target concentration (parameter
-:math:`c_0` in Eqs.  (20-21) of :ref:`Sadigh1 <Sadigh1>`). The atomic
-concentrations refer to components `B`, `C` ..., with `A` being set
-automatically. When the simulation includes `N` elements, `N-1`
-concentration values must be specified.
+The parameter *conc* sets the `N-1` target atomic concentration
+fractions (parameter :math:`c_0` in Eqs.  (20-21) of :ref:`Sadigh1 <Sadigh1>`)
+:math:`0 \le c_2, \ldots, c_N \le 1`, with
+:math:`c_1 = 1 - \Sigma_{i=2}^N c_i`.
+When the simulation includes `N` atom types (elements),
+`N-1` concentration values must be specified.
 
 ------------
 
@@ -143,10 +142,12 @@ components of the vector represent the following quantities:
 
 * 1 = The absolute number of accepted trial swaps during the last MC step
 * 2 = The absolute number of rejected trial swaps during the last MC step
-* 3 = The current global concentration of species *A* (= number of atoms of type 1 / total number of atoms)
-* 4 = The current global concentration of species *B* (= number of atoms of type 2 / total number of atoms)
+* 3 = Current global concentration `c_1` (= number of atoms of type 1 / total number of atoms)
+* 4 = Current global concentration `c_2` (= number of atoms of type 2 / total number of atoms)
 * ...
-* N+2: The current global concentration of species *X* (= number of atoms of type *N* / total number of atoms)
+* N+2: Current global concentration `c_N` (= number of atoms of type *N* / total number of atoms)
+
+The vector values calculated by this fix are "intensive".
 
 Restrictions
 """"""""""""

doc/src/fix_wall_gran.rst

@@ -115,6 +115,18 @@ friction and twisting friction supported by the :doc:`pair_style granular <pair_
 supported for walls. These are discussed in greater detail on the doc
 page for :doc:`pair_style granular <pair_granular>`.
 
+.. note::
+   When *fstyle* *granular* is specified, the associated *fstyle_params* are taken as
+   those for a wall/particle interaction. For example, for the *hertz/material* normal
+   contact model with :math:`E = 960` and :math:`\nu = 0.2`, the effective Young's
+   modulus for a wall/particle interaction is computed as
+   :math:`E_{eff} = \frac{960}{2(1-0.2^2)} = 500`. Any pair coefficients defined by
+   :doc:`pair_style granular <pair_granular>` are not taken into consideration. To
+   model different wall/particle interactions for particles of different material
+   types, the user may define multiple fix wall/gran commands operating on separate
+   groups (e.g. based on particle type) each with a different wall/particle effective
+   Young's modulus.
+
 Note that you can choose a different force styles and/or different
 values for the wall/particle coefficients than for particle/particle
 interactions.  E.g. if you wish to model the wall as a different

doc/src/fix_widom.rst

@@ -14,7 +14,7 @@ Syntax
 * widom = style name of this fix command
 * N = invoke this fix every N steps
 * M = number of Widom insertions to attempt every N steps
-* type = atom type for inserted atoms (must be 0 if mol keyword used)
+* type = atom type (1-Ntypes or type label) for inserted atoms (must be 0 if mol keyword used)
 * seed = random # seed (positive integer)
 * T = temperature of the system (temperature units)
 * zero or more keyword/value pairs may be appended to args
@@ -38,6 +38,9 @@ Examples
    fix 2 gas widom 1 50000 1 19494 2.0
    fix 3 water widom 1000 100 0 29494 300.0 mol h2omol full_energy
 
+   labelmap atom 1 Li
+   fix 2 ion widom 1 50000 Li 19494 2.0
+
 Description
 """""""""""
 
@@ -179,7 +182,7 @@ the following global cumulative quantities:
 * 2 = average difference in potential energy on each timestep
 * 3 = volume of the insertion region
 
-The vector values calculated by this fix are "extensive".
+The vector values calculated by this fix are "intensive".
 
 No parameter of this fix can be used with the *start/stop* keywords of
 the :doc:`run <run>` command.  This fix is not invoked during

doc/src/group.rst

@@ -20,13 +20,13 @@ Syntax
        *empty* = no args
        *region* args = region-ID
        *type* or *id* or *molecule*
-         args = list of one or more atom types, atom IDs, or molecule IDs
-           any entry in list can be a sequence formatted as A:B or A:B:C where
+         args = list of one or more atom types (1-Ntypes or type label), atom IDs, or molecule IDs
+           any numeric entry in list can be a sequence formatted as A:B or A:B:C where
            A = starting index, B = ending index,
            C = increment between indices, 1 if not specified
          args = logical value
            logical = "<" or "<=" or ">" or ">=" or "==" or "!="
-           value = an atom type or atom ID or molecule ID (depending on *style*\ )
+           value = an atom type (1-Ntypes or type label) or atom ID or molecule ID (depending on *style*\ )
          args = logical value1 value2
            logical = "<>"
            value1,value2 = atom types or atom IDs or molecule IDs (depending on *style*\ )
@@ -52,6 +52,7 @@ Examples
 
    group edge region regstrip
    group water type 3 4
+   group water type OW HT
    group sub id 10 25 50
    group sub id 10 25 50 500:1000
    group sub id 100:10000:10
@@ -119,7 +120,7 @@ three styles can use arguments specified in one of two formats.
 
 The first format is a list of values (types or IDs).  For example, the
 second command in the examples above puts all atoms of type 3 or 4 into
-the group named *water*\ .  Each entry in the list can be a
+the group named *water*\ .  Each numeric entry in the list can be a
 colon-separated sequence ``A:B`` or ``A:B:C``, as in two of the examples
 above.  A "sequence" generates a sequence of values (types or IDs),
 with an optional increment.  The first example with ``500:1000`` has the
@@ -135,7 +136,8 @@ except ``<>`` take a single argument.  The third example above adds all
 atoms with IDs from 1 to 150 to the group named *sub*\ .  The logical ``<>``
 means "between" and takes 2 arguments.  The fourth example above adds all
 atoms belonging to molecules with IDs from 50 to 250 (inclusive) to
-the group named polyA.
+the group named polyA.  For the *type* style, type labels are converted into
+numeric types before being evaluated.
 
 The *variable* style evaluates a variable to determine which atoms to
 add to the group.  It must be an :doc:`atom-style variable <variable>`

doc/src/group2ndx.rst

@@ -34,32 +34,66 @@ Description
 Write or read a Gromacs style index file in text format that associates
 atom IDs with the corresponding group definitions. This index file can be
 used with in combination with Gromacs analysis tools or to import group
-definitions into the :doc:`fix colvars <fix_colvars>` input file. It can
-also be used to save and restore group definitions for static groups.
+definitions into the :doc:`fix colvars <fix_colvars>` input file.
+
+It can also be used to save and restore group definitions for static groups
+using the individual atom IDs. This may be important if the original
+group definition depends on a region or otherwise on the geometry and thus
+cannot be easily recreated.
+
+Another application would be to import atom groups defined for Gromacs
+simulation into LAMMPS.  When translating Gromacs topology and geometry
+data to LAMMPS.
 
 The *group2ndx* command will write group definitions to an index file.
-Without specifying any group IDs, all groups will be written to the index
-file. When specifying group IDs, only those groups will be written to the
-index file. In order to follow the Gromacs conventions, the group *all*
-will be renamed to *System* in the index file.
+Without specifying any group IDs, all groups will be written to the
+index file.  When specifying group IDs, only those groups will be
+written to the index file.  In order to follow the Gromacs conventions,
+the group *all* will be renamed to *System* in the index file.
 
-The *ndx2group* command will create of update group definitions from those
-stored in an index file. Without specifying any group IDs, all groups except
-*System* will be read from the index file and the corresponding groups
-recreated. If a group of the same name already exists, it will be completely
-reset. When specifying group IDs, those groups, if present, will be read
-from the index file and restored.
+The *ndx2group* command will create of update group definitions from
+those stored in an index file.  Without specifying any group IDs, all
+groups except *System* will be read from the index file and the
+corresponding groups recreated.  If a group of the same name already
+exists, it will be completely reset.  When specifying group IDs, those
+groups, if present, will be read from the index file and restored.
+
+File Format
+"""""""""""
+
+The file format is equivalent and compatible with what is produced by
+the `Gromacs make_ndx command <https://manual.gromacs.org/current/onlinehelp/gmx-make_ndx.html>`_.
+and follows the `Gromacs definition of an ndx file <https://manual.gromacs.org/current/reference-manual/file-formats.html#ndx>`_
+
+Each group definition begins with the group name in square brackets with
+blanks, e.g. \[ water \] and is then followed by the list of atom
+indices, which may be spread over multiple lines.  Here is a small
+example file:
+
+.. code-block:: ini
+
+   [ Oxygen ]
+   1  4  7
+   [ Hydrogen ]
+   2  3  5  6
+   8  9
+   [ Water ]
+   1 2 3 4 5 6 7 8 9
+
+The index file defines 3 groups: Oxygen, Hydrogen, and Water and the
+latter happens to be the union of the first two.
 
 ----------
 
 Restrictions
 """"""""""""
 
-This command requires that atoms have atom IDs, since this is the
+These commands require that atoms have atom IDs, since this is the
 information that is written to the index file.
 
-These commands are part of the COLVARS package.  They are only
-enabled if LAMMPS was built with that package.  See the :doc:`Build package <Build_package>` page for more info.
+These commands are part of the EXTRA-COMMAND package.  They are only
+enabled if LAMMPS was built with that package.  See the
+:doc:`Build package <Build_package>` page for more info.
 
 Related commands
 """"""""""""""""

doc/src/improper_cossq.rst

@@ -64,8 +64,8 @@ Restrictions
 """"""""""""
 
 This improper style can only be used if LAMMPS was built with the
-MOLECULE package.  See the :doc:`Build package <Build_package>` doc
-page for more info.
+EXTRA-MOLECULE package.  See the :doc:`Build package <Build_package>`
+doc page for more info.
 
 Related commands
 """"""""""""""""

doc/src/improper_distance.rst

@@ -54,8 +54,8 @@ Restrictions
 """"""""""""
 
 This improper style can only be used if LAMMPS was built with the
-MOLECULE package.  See the :doc:`Build package <Build_package>` doc
-page for more info.
+EXTRA-MOLECULE package.  See the :doc:`Build package <Build_package>`
+doc page for more info.
 
 Related commands
 """"""""""""""""

doc/src/improper_fourier.rst

@@ -60,8 +60,8 @@ Restrictions
 """"""""""""
 
 This angle style can only be used if LAMMPS was built with the
-MOLECULE package.  See the :doc:`Build package <Build_package>` doc
-page for more info.
+EXTRA-MOLECULE package.  See the :doc:`Build package <Build_package>`
+doc page for more info.
 
 Related commands
 """"""""""""""""

doc/src/improper_ring.rst

@@ -72,8 +72,8 @@ Restrictions
 """"""""""""
 
 This improper style can only be used if LAMMPS was built with the
-MOLECULE package.  See the :doc:`Build package <Build_package>` doc
-page for more info.
+EXTRA-MOLECULE package.  See the :doc:`Build package <Build_package>`
+doc page for more info.
 
 Related commands
 """"""""""""""""

doc/src/labelmap.rst

@@ -24,6 +24,7 @@ Examples
 
 .. code-block:: LAMMPS
 
+   labelmap atom 1 c1 2 hc 3 cp 4 nt
    labelmap atom 3 carbon 4 'c3"' 5 "c1'" 6 "c#"
    labelmap atom $(label2type(atom,carbon)) C  # change type label from 'carbon' to 'C'
    labelmap clear
@@ -43,8 +44,8 @@ The label map can also be defined by the :doc:`read_data <read_data>`
 command when it reads these sections in a data file: Atom Type Labels,
 Bond Type Labels, etc.  See the :doc:`Howto type labels
 <Howto_type_labels>` doc page for a general discussion of how type
-labels can be used.  See :ref:`(Gissinger) <Typelabel>` for a discussion
-of the type label implementation in LAMMPS and its uses.
+labels can be used.  See :ref:`(Gissinger) <Typelabel1>` for a
+discussion of the type label implementation in LAMMPS and its uses.
 
 Valid type labels can contain any alphanumeric character, but must not
 start with a number, a '#', or a '*' character.  They can contain other
@@ -102,6 +103,6 @@ none
 
 -----------
 
-.. _Typelabel:
+.. _Typelabel1:
 
 **(Gissinger)** J. R. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M. Choi, D. S. Karls, A. Stukowski, W. Im, H. Heinz, A. Kohlmeyer, and E. B. Tadmor, J Phys Chem B, 128, 3282-3297 (2024).

doc/src/pair_charmm.rst

@@ -112,26 +112,22 @@ Description
 These pair styles compute Lennard Jones (LJ) and Coulombic
 interactions with additional switching or shifting functions that ramp
 the energy and/or force smoothly to zero between an inner and outer
-cutoff.  They are implementations of the widely used CHARMM force
-field used in the `CHARMM <https://www.charmm.org>`_ MD code (and
-others).  See :ref:`(MacKerell) <pair-MacKerell>` for a description of the
-CHARMM force field.
+cutoff.  They implement the widely used CHARMM force field, see
+:doc:`Howto discussion on biomolecular force fields <Howto_bioFF>` for
+details.
 
 The styles with *charmm* (not *charmmfsw* or *charmmfsh*\ ) in their
 name are the older, original LAMMPS implementations.  They compute the
-LJ and Coulombic interactions with an energy switching function (esw,
-shown in the formula below as S(r)), which ramps the energy smoothly
-to zero between the inner and outer cutoff.  This can cause
-irregularities in pairwise forces (due to the discontinuous second
-derivative of energy at the boundaries of the switching region), which
-in some cases can result in detectable artifacts in an MD simulation.
+LJ and Coulombic interactions with an energy switching function which
+ramps the energy smoothly to zero between the inner and outer cutoff.
+This can cause irregularities in pairwise forces (due to the discontinuous
+second derivative of energy at the boundaries of the switching region),
+which in some cases can result in detectable artifacts in an MD simulation.
 
 The newer styles with *charmmfsw* or *charmmfsh* in their name replace
 the energy switching with force switching (fsw) and force shifting
 (fsh) functions, for LJ and Coulombic interactions respectively.
-These follow the formulas and description given in
-:ref:`(Steinbach) <Steinbach>` and :ref:`(Brooks) <Brooks1>` to minimize these
-artifacts.
+
 
 .. note::
 
@@ -152,26 +148,6 @@ artifacts.
    the CHARMM force field energies and forces, when using one of these
    two CHARMM pair styles.
 
-.. math::
-
-   E = & LJ(r) \qquad \qquad \qquad r < r_{\rm in} \\
-     = & S(r) * LJ(r) \qquad \qquad r_{\rm in} < r < r_{\rm out} \\
-     = & 0 \qquad \qquad \qquad \qquad r > r_{\rm out} \\
-   E = & C(r) \qquad \qquad \qquad r < r_{\rm in} \\
-     = & S(r) * C(r) \qquad \qquad r_{\rm in} < r < r_{\rm out} \\
-     = & 0 \qquad \qquad \qquad \qquad r > r_{\rm out} \\
-   LJ(r) = & 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} -
-           \left(\frac{\sigma}{r}\right)^6 \right] \\
-   C(r) = & \frac{C q_i q_j}{ \epsilon r} \\
-   S(r) = & \frac{ \left[r_{\rm out}^2 - r^2\right]^2
-     \left[r_{\rm out}^2 + 2r^2 - 3{r_{\rm in}^2}\right]}
-   { \left[r_{\rm out}^2 - {r_{\rm in}}^2\right]^3 }
-
-where S(r) is the energy switching function mentioned above for the
-*charmm* styles.  See the :ref:`(Steinbach) <Steinbach>` paper for the
-functional forms of the force switching and force shifting functions
-used in the *charmmfsw* and *charmmfsh* styles.
-
 When using the *lj/charmm/coul/charmm styles*, both the LJ and
 Coulombic terms require an inner and outer cutoff. They can be the
 same for both formulas or different depending on whether 2 or 4

doc/src/pair_coul.rst

@@ -379,10 +379,11 @@ These pair styles can only be used via the *pair* keyword of the
 Restrictions
 """"""""""""
 
-The *coul/cut/global*, *coul/long*, *coul/msm*, *coul/streitz*, and *tip4p/long* styles
-are part of the KSPACE package.  They are only enabled if LAMMPS was built
-with that package.  See the :doc:`Build package <Build_package>` doc page
-for more info.
+The *coul/long*, *coul/msm*, *coul/streitz*, and *tip4p/long* styles are
+part of the KSPACE package.  The *coul/cut/global* and *coul/exclude* are
+part of the EXTRA-PAIR package.  A pair style is only enabled if LAMMPS was
+built with its corresponding package.  See the :doc:`Build package <Build_package>`
+doc page for more info.
 
 Related commands
 """"""""""""""""

doc/src/pair_dpd_coul_slater_long.rst

@@ -0,0 +1,196 @@
+.. index:: pair_style dpd/coul/slater/long
+.. index:: pair_style dpd/coul/slater/long/gpu
+
+pair_style dpd/coul/slater/long command
+=======================================
+
+Accelerator Variants: *dpd/coul/slater/long/gpu*
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   pair_style dpd/coul/slater/long T cutoff_DPD seed lambda cutoff_coul
+
+* T = temperature (temperature units)
+* cutoff_DPD = global cutoff for DPD interactions (distance units)
+* seed = random # seed (positive integer)
+* lambda = decay length of the charge (distance units)
+* cutoff_coul = global cutoff for Coulombic interactions (distance units)
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   pair_style dpd/coul/slater/long 1.0 2.5 34387 0.25 3.0
+
+   pair_coeff 1 1 78.0 4.5          # not charged by default
+   pair_coeff 2 2 78.0 4.5 yes
+
+
+Description
+"""""""""""
+
+.. versionadded:: 27June2024
+
+Style *dpd/coul/slater/long* computes a force field for dissipative
+particle dynamics (DPD) following the exposition in :ref:`(Groot)
+<Groot5>`.  It also allows for the use of charged particles in the
+model by adding a long-range Coulombic term to the DPD interactions.
+The short-range portion of the Coulombics is calculated by this pair
+style.  The long-range Coulombics are computed by use of the
+:doc:`kspace_style <kspace_style>` command, e.g. using the Ewald or
+PPPM styles.
+
+Coulombic forces in mesoscopic models such as DPD employ potentials
+without explicit excluded-volume interactions.  The goal is to prevent
+artificial ionic pair formation by including a charge distribution in
+the Coulomb potential, following the formulation in :ref:`(Melchor1)
+<Melchor1>`.
+
+.. note::
+
+   This pair style is effectively the combination of the
+   :doc:`pair_style dpd <pair_dpd>` and :doc:`pair_style
+   coul/slater/long <pair_coul_slater>` commands, but should be more
+   efficient (especially on GPUs) than using :doc:`pair_style
+   hybrid/overlay dpd coul/slater/long <pair_hybrid>`. That is
+   particularly true for the GPU package version of the pair style since
+   this version is compatible with computing neighbor lists on the GPU
+   instead of the CPU as is required for hybrid styles.
+
+In the charged DPD model, the force on bead I due to bead J is given
+as a sum of 4 terms:
+
+.. math::
+
+   \vec{f}  = & (F^C + F^D + F^R + F^E) \hat{r_{ij}} \\
+   F^C      = & A w(r) \qquad \qquad \qquad \qquad \qquad r < r_{DPD} \\
+   F^D      = & - \gamma w^2(r) (\hat{r_{ij}} \bullet \vec{v}_{ij}) \qquad \qquad r < r_{DPD} \\
+   F^R      = & \sigma w(r) \alpha (\Delta t)^{-1/2} \qquad \qquad \qquad r < r_{DPD} \\
+   w(r)     = & 1 - \frac{r}{r_{DPD}} \\
+   F^E      = & \frac{C q_iq_j}{\epsilon r^2} \left( 1- exp\left( \frac{2r_{ij}}{\lambda} \right) \left( 1 + \frac{2r_{ij}}{\lambda} \left( 1 + \frac{r_{ij}}{\lambda} \right)\right) \right)
+
+where :math:`F^C` is a conservative force, :math:`F^D` is a
+dissipative force, :math:`F^R` is a random force, and :math:`F^E` is
+an electrostatic force.  :math:`\hat{r_{ij}}` is a unit vector in the
+direction :math:`r_i - r_j`, :math:`\vec{v}_{ij}` is the vector
+difference in velocities of the two atoms :math:`\vec{v}_i -
+\vec{v}_j`, :math:`\alpha` is a Gaussian random number with zero mean
+and unit variance, *dt* is the timestep size, and :math:`w(r)` is a
+weighting factor that varies between 0 and 1.
+
+:math:`\sigma` is set equal to :math:`\sqrt{2 k_B T \gamma}`, where
+:math:`k_B` is the Boltzmann constant and *T* is the temperature
+parameter in the pair_style command.
+
+:math:`r_{DPD}` is the pairwise cutoff for the first 3 DPD terms in
+the formula as specified by *cutoff_DPD*.  For the :math:`F^E` term,
+pairwise interactions within the specified *cutoff_coul* distance are
+computed directly; interactions beyond that distance are computed in
+reciprocal space.  *C* is the same Coulomb conversion factor used in
+the Coulombic formulas described on the :doc:`pair_coul <pair_coul>`
+doc page.
+
+The following parameters must be defined for each pair of atoms types
+via the :doc:`pair_coeff <pair_coeff>` command as in the examples
+above, or in the data file or restart files read by the
+:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
+commands:
+
+* A (force units)
+* :math:`\gamma` (force/velocity units)
+* is_charged (optional boolean, default = no)
+
+The *is_charged* parameter is optional and can be specified as *yes* or
+*no*.  *Yes* should be used for interactions between two types of
+charged particles.  *No* is the default and should be used for
+interactions between two types of particles when one or both are
+uncharged.
+
+----------
+
+.. include:: accel_styles.rst
+
+----------
+
+Mixing, shift, table, tail correction, restart, rRESPA info
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+This pair style does not support mixing.  Thus, coefficients for all
+I,J pairs must be specified explicitly.
+
+This pair style does not support the :doc:`pair_modify <pair_modify>`
+shift option for the energy of the pair interaction.
+
+The :doc:`pair_modify <pair_modify>` table option is not relevant
+for this pair style.
+
+This pair style does not support the :doc:`pair_modify <pair_modify>`
+tail option for adding long-range tail corrections to energy and
+pressure.
+
+This pair style writes its information to :doc:`binary restart files
+<restart>`, so pair_style and pair_coeff commands do not need to be
+specified in an input script that reads a restart file.  Note that the
+user-specified random number seed is stored in the restart file, so
+when a simulation is restarted, each processor will re-initialize its
+random number generator the same way it did initially.  This means the
+random forces will be random, but will not be the same as they would
+have been if the original simulation had continued past the restart
+time.
+
+This pair style can only be used via the *pair* keyword of the
+:doc:`run_style respa <run_style>` command.  It does not support the
+*inner*, *middle*, *outer* keywords.
+
+----------
+
+Restrictions
+""""""""""""
+
+This style is part of the DPD-BASIC package.  It is only enabled if
+LAMMPS was built with that package.  See the :doc:`Build package
+<Build_package>` page for more info.
+
+The default frequency for rebuilding neighbor lists is every 10 steps
+(see the :doc:`neigh_modify <neigh_modify>` command). This may be too
+infrequent since particles move rapidly and can overlap by large
+amounts.  If this setting yields a non-zero number of "dangerous"
+reneighborings (printed at the end of a simulation), you should
+experiment with forcing reneighboring more often and see if system
+energies/trajectories change.
+
+This pair style requires use of the :doc:`comm_modify vel yes
+<comm_modify>` command so that velocities are stored by ghost atoms.
+
+This pair style also requires use of a long-range solvers from the
+KSPACE package.
+
+This pair style will not restart exactly when using the
+:doc:`read_restart <read_restart>` command, though they should provide
+statistically similar results.  This is because the forces they compute
+depend on atom velocities.  See the :doc:`read_restart <read_restart>`
+command for more details.
+
+Related commands
+""""""""""""""""
+
+:doc:`pair_style dpd <pair_dpd>`, :doc:`pair_style coul/slater/long <pair_coul_slater>`,
+
+Default
+"""""""
+
+For the pair_coeff command, the default is is_charged = no.
+
+----------
+
+.. _Groot5:
+
+**(Groot)** Groot and Warren, J Chem Phys, 107, 4423-35 (1997).
+
+.. _Melchor1:
+
+**(Melchor)** Gonzalez-Melchor, Mayoral, Velazquez, and Alejandre, J Chem Phys, 125, 224107 (2006).

doc/src/pair_granular.rst

@@ -111,7 +111,7 @@ For the *hertz* model, the normal component of force is given by:
 
    \mathbf{F}_{ne, Hertz} = k_n R_{eff}^{1/2}\delta_{ij}^{3/2} \mathbf{n}
 
-Here, :math:`R_{eff} = \frac{R_i R_j}{R_i + R_j}` is the effective
+Here, :math:`R_{eff} = R = \frac{R_i R_j}{R_i + R_j}` is the effective
 radius, denoted for simplicity as *R* from here on.  For *hertz*, the
 units of the spring constant :math:`k_n` are *force*\ /\ *length*\ \^2, or
 equivalently *pressure*\ .
@@ -120,13 +120,14 @@ For the *hertz/material* model, the force is given by:
 
 .. math::
 
-   \mathbf{F}_{ne, Hertz/material} = \frac{4}{3} E_{eff} R_{eff}^{1/2}\delta_{ij}^{3/2} \mathbf{n}
+   \mathbf{F}_{ne, Hertz/material} = \frac{4}{3} E_{eff} R^{1/2}\delta_{ij}^{3/2} \mathbf{n}
 
-Here, :math:`E_{eff} = E = \left(\frac{1-\nu_i^2}{E_i} + \frac{1-\nu_j^2}{E_j}\right)^{-1}` is the effective Young's
-modulus, with :math:`\nu_i, \nu_j` the Poisson ratios of the particles of
-types *i* and *j*\ . Note that if the elastic modulus and the shear
-modulus of the two particles are the same, the *hertz/material* model
-is equivalent to the *hertz* model with :math:`k_n = 4/3 E_{eff}`
+Here, :math:`E_{eff} = E = \left(\frac{1-\nu_i^2}{E_i} + \frac{1-\nu_j^2}{E_j}\right)^{-1}`
+is the effective Young's modulus, with :math:`\nu_i, \nu_j` the Poisson ratios
+of the particles of types *i* and *j*. :math:`E_{eff}` is denoted as *E* from here on.
+Note that if the elastic modulus and the shear modulus of the two particles are the
+same, the *hertz/material* model is equivalent to the *hertz* model with
+:math:`k_n = 4/3 E`
 
 The *dmt* model corresponds to the
 :ref:`(Derjaguin-Muller-Toporov) <DMT1975>` cohesive model, where the force
@@ -187,6 +188,7 @@ for the damping model currently supported are:
 2. *mass_velocity*
 3. *viscoelastic*
 4. *tsuji*
+5. *coeff_restitution*
 
 If the *damping* keyword is not specified, the *viscoelastic* model is
 used by default.
@@ -248,6 +250,30 @@ The dimensionless coefficient of restitution :math:`e` specified as part
 of the normal contact model parameters should be between 0 and 1, but
 no error check is performed on this.
 
+The *coeff_restitution* model is useful when a specific normal coefficient of
+restitution :math:`e` is required. In these models, the normal coefficient of
+restitution :math:`e` is specified as an input. Following the approach of
+:ref:`(Brilliantov et al) <Brill1996>`, when using the *hooke* normal model,
+*coeff_restitution* calculates the damping coefficient as:
+
+.. math::
+
+   \eta_n = \sqrt{\frac{4m_{eff}k_n}{1+\left( \frac{\pi}{\log(e)}\right)^2}} ,
+
+For any other normal model, e.g. the *hertz* and *hertz/material* models, the damping
+coefficient is:
+
+.. math::
+
+   \eta_n = -2\sqrt{\frac{5}{6}}\frac{\log(e)}{\sqrt{\pi^2+(\log(e))^2}}(R_{eff} \delta_{ij})^{\frac{1}{4}}\sqrt{\frac{3}{2}k_n m_{eff}} ,
+
+where :math:`k_n = \frac{4}{3} E_{eff}` for the *hertz/material* model. Since
+*coeff_restitution* accounts for the effective mass, effective radius, and
+pairwise overlaps (except when used with the *hooke* normal model) when calculating
+the damping coefficient, it accurately reproduces the specified coefficient of
+restitution for both monodisperse and polydisperse particle pairs.  This damping
+model is not compatible with cohesive normal models such as *JKR* or *DMT*.
+
 The total normal force is computed as the sum of the elastic and
 damping components:
 
@@ -417,11 +443,11 @@ discussion above. To match the Mindlin solution, one should set
 
    G_{eff} = \left(\frac{2-\nu_i}{G_i} + \frac{2-\nu_j}{G_j}\right)^{-1}
 
-where :math:`G` is the shear modulus, related to Young's modulus :math:`E`
-and Poisson's ratio :math:`\nu` by :math:`G = E/(2(1+\nu))`. This can also be
-achieved by specifying *NULL* for :math:`k_t`, in which case a
-normal contact model that specifies material parameters :math:`E` and
-:math:`\nu` is required (e.g. *hertz/material*, *dmt* or *jkr*\ ). In this
+where :math:`G_i` is the shear modulus of a particle of type :math:`i`, related to Young's
+modulus :math:`E_i` and Poisson's ratio :math:`\nu_i` by :math:`G_i = E_i/(2(1+\nu_i))`.
+This can also be achieved by specifying *NULL* for :math:`k_t`, in which case a
+normal contact model that specifies material parameters :math:`E_i` and
+:math:`\nu_i` is required (e.g. *hertz/material*, *dmt* or *jkr*\ ). In this
 case, mixing of the shear modulus for different particle types *i* and
 *j* is done according to the formula above.
 
@@ -551,7 +577,7 @@ opposite torque on each particle, according to:
 
 .. math::
 
-   \tau_{roll,i} =  R_{eff} \mathbf{n} \times \mathbf{F}_{roll}
+   \tau_{roll,i} =  R \mathbf{n} \times \mathbf{F}_{roll}
 
 .. math::
 

doc/src/pair_hybrid.rst

@@ -1,28 +1,41 @@
 .. index:: pair_style hybrid
 .. index:: pair_style hybrid/kk
+.. index:: pair_style hybrid/omp
+.. index:: pair_style hybrid/molecular
+.. index:: pair_style hybrid/molecular/omp
 .. index:: pair_style hybrid/overlay
+.. index:: pair_style hybrid/overlay/omp
 .. index:: pair_style hybrid/overlay/kk
 .. index:: pair_style hybrid/scaled
+.. index:: pair_style hybrid/scaled/omp
 
 pair_style hybrid command
 =========================
 
-Accelerator Variants: *hybrid/kk*
+Accelerator Variants: *hybrid/kk*, *hybrid/omp*
+
+pair_style hybrid/molecular command
+===================================
+
+Accelerator Variant: *hybrid/molecular/omp*
 
 pair_style hybrid/overlay command
 =================================
 
-Accelerator Variants: *hybrid/overlay/kk*
+Accelerator Variants: *hybrid/overlay/kk*, *hybrid/overlay/omp*
 
 pair_style hybrid/scaled command
 ==================================
 
+Accelerator Variant: *hybrid/scaled/omp*
+
 Syntax
 """"""
 
 .. code-block:: LAMMPS
 
    pair_style hybrid style1 args style2 args ...
+   pair_style hybrid/molecular factor1 style1 args factor2 style 2 args
    pair_style hybrid/overlay style1 args style2 args ...
    pair_style hybrid/scaled factor1 style1 args factor2 style 2 args ...
 
@@ -47,6 +60,10 @@ Examples
    pair_coeff * * tersoff Si.tersoff Si
    pair_coeff * * sw Si.sw Si
 
+   pair_style hybrid/molecular lj/cut 2.5 lj/cut 2.5
+   pair_coeff * * lj/cut 1 1.0 1.0
+   pair_coeff * * lj/cut 2 1.5 1.0
+
    variable one equal ramp(1.0,0.0)
    variable two equal 1.0-v_one
    pair_style hybrid/scaled v_one lj/cut 2.5 v_two morse 2.5
@@ -56,17 +73,26 @@ Examples
 Description
 """""""""""
 
-The *hybrid*, *hybrid/overlay*, and *hybrid/scaled* styles enable the
-use of multiple pair styles in one simulation.  With the *hybrid* style,
-exactly one pair style is assigned to each pair of atom types.  With the
-*hybrid/overlay* and *hybrid/scaled* styles, one or more pair styles can
-be assigned to each pair of atom types.  The assignment of pair styles
-to type pairs is made via the :doc:`pair_coeff <pair_coeff>` command.
-The major difference between the *hybrid/overlay* and *hybrid/scaled*
-styles is that the *hybrid/scaled* adds a scale factor for each
-sub-style contribution to forces, energies and stresses.  Because of the
-added complexity, the *hybrid/scaled* style has more overhead and thus
-may be slower than *hybrid/overlay*.
+The *hybrid*, *hybrid/overlay*, *hybrid/molecular*, and *hybrid/scaled*
+styles enable the use of multiple pair styles in one simulation.  With
+the *hybrid* style, exactly one pair style is assigned to each pair of
+atom types.  With the *hybrid/overlay* and *hybrid/scaled* styles, one
+or more pair styles can be assigned to each pair of atom types.  With
+the hybrid/molecular style, pair styles are assigned to either intra-
+or inter-molecular interactions.
+
+The assignment of pair styles to type pairs is made via the
+:doc:`pair_coeff <pair_coeff>` command.  The major difference between
+the *hybrid/overlay* and *hybrid/scaled* styles is that the
+*hybrid/scaled* adds a scale factor for each sub-style contribution to
+forces, energies and stresses.  Because of the added complexity, the
+*hybrid/scaled* style has more overhead and thus may be slower than
+*hybrid/overlay*.
+
+The *hybrid/molecular* pair style accepts *only* two sub-styles: the
+first is assigned to intra-molecular interactions (i.e. both atoms
+have the same molecule ID), the second to inter-molecular interactions
+(i.e. interacting atoms have different molecule IDs).
 
 Here are two examples of hybrid simulations.  The *hybrid* style could
 be used for a simulation of a metal droplet on a LJ surface.  The metal
@@ -476,6 +502,8 @@ the same or else LAMMPS will generate an error.
 Pair style *hybrid/scaled* currently only works for non-accelerated
 pair styles and pair styles from the OPT package.
 
+Pair style *hybrid/molecular* is not compatible with manybody potentials.
+
 When using pair styles from the GPU package they must not be listed
 multiple times.  LAMMPS will detect this and abort.
 

doc/src/pair_oxdna.rst

@@ -37,18 +37,19 @@ Syntax
 
      *oxdna/stk* args = seq T xi kappa 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
        seq = seqav (for average sequence stacking strength) or seqdep (for sequence-dependent stacking strength)
-       T = temperature (oxDNA units, 0.1 = 300 K)
-       xi = 1.3448 (temperature-independent coefficient in stacking strength)
-       kappa = 2.6568 (coefficient of linear temperature dependence in stacking strength)
+       T = temperature (LJ units: 0.1 = 300 K, real units: 300 = 300 K)
+       xi = 1.3448 (LJ units) or 8.01727944817084 (real units), temperature-independent coefficient in stacking strength
+       kappa = 2.6568 (LJ units) or 0.005279604 (real units), coefficient of linear temperature dependence in stacking strength
      *oxdna/hbond* args = seq eps 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
        seq = seqav (for average sequence base-pairing strength) or seqdep (for sequence-dependent base-pairing strength)
-       eps = 1.077 (between base pairs A-T and C-G) or 0 (all other pairs)
+       eps = 1.077 (LJ units) or 6.42073911784652 (real units), average hydrogen bonding strength between A-T and C-G Watson-Crick base pairs, 0 between all other pairs
 
 Examples
 """"""""
 
 .. code-block:: LAMMPS
 
+   # LJ units
    pair_style hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk
    pair_coeff * * oxdna/excv    2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
    pair_coeff * * oxdna/stk     seqdep 0.1 1.3448 2.6568 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
@@ -58,55 +59,105 @@ Examples
    pair_coeff * * oxdna/xstk    47.5 0.575 0.675 0.495 0.655 2.25 0.791592653589793 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0 0.65 1.7 0.875 0.68 1.7 0.875 0.68
    pair_coeff * * oxdna/coaxstk 46.0 0.4 0.6 0.22 0.58 2.0 2.541592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 -0.65 2.0 -0.65
 
+   pair_style hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk
+   pair_coeff * * oxdna/excv    oxdna_lj.cgdna
+   pair_coeff * * oxdna/stk     seqav 0.1 1.3448 2.6568 oxdna_lj.cgdna
+   pair_coeff * * oxdna/hbond   seqav oxdna_lj.cgdna
+   pair_coeff 1 4 oxdna/hbond   seqav oxdna_lj.cgdna
+   pair_coeff 2 3 oxdna/hbond   seqav oxdna_lj.cgdna
+   pair_coeff * * oxdna/xstk    oxdna_lj.cgdna
+   pair_coeff * * oxdna/coaxstk oxdna_lj.cgdna
+
+   # Real units
+   pair_style hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk
+   pair_coeff * * oxdna/excv    11.92337812042065 5.9626 5.74965 11.92337812042065 4.38677 4.259 11.92337812042065 2.81094 2.72576
+   pair_coeff * * oxdna/stk     seqdep 300.0 8.01727944817084 0.005279604 0.70439070204273 3.4072 7.6662 2.72576 6.3885 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 2.0 0.65 2.0 0.65
+   pair_coeff * * oxdna/hbond   seqdep 0.0 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45
+   pair_coeff 1 4 oxdna/hbond   seqdep 6.42073911784652 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45
+   pair_coeff 2 3 oxdna/hbond   seqdep 6.42073911784652 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45
+   pair_coeff * * oxdna/xstk    3.9029021145006 4.89785 5.74965 4.21641 5.57929 2.25 0.791592654 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0.0 0.65 1.7 0.875 0.68 1.7 0.875 0.68
+   pair_coeff * * oxdna/coaxstk 3.77965257404268 3.4072 5.1108 1.87396 4.94044 2.0 2.541592654 0.65 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 2.0 -0.65 2.0 -0.65
+
+   pair_style hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk
+   pair_coeff * * oxdna/excv    oxdna_real.cgdna
+   pair_coeff * * oxdna/stk     seqav 300.0 8.01727944817084 0.005279604 oxdna_real.cgdna
+   pair_coeff * * oxdna/hbond   seqav oxdna_real.cgdna
+   pair_coeff 1 4 oxdna/hbond   seqav oxdna_real.cgdna
+   pair_coeff 2 3 oxdna/hbond   seqav oxdna_real.cgdna
+   pair_coeff * * oxdna/xstk    oxdna_real.cgdna
+   pair_coeff * * oxdna/coaxstk oxdna_real.cgdna
+
+.. note::
+
+   The coefficients in the above examples are provided in forms
+   compatible with both *units lj* and *units real* (see documentation
+   of :doc:`units <units>`).  These can also be read from a potential
+   file with correct unit style by specifying the name of the
+   file. Several potential files for each unit style are included in the
+   ``potentials`` directory of the LAMMPS distribution.
+
 Description
 """""""""""
 
-The *oxdna* pair styles compute the pairwise-additive parts of the oxDNA force field
-for coarse-grained modelling of DNA. The effective interaction between the nucleotides consists of potentials for the
-excluded volume interaction *oxdna/excv*, the stacking *oxdna/stk*, cross-stacking *oxdna/xstk*
-and coaxial stacking interaction *oxdna/coaxstk* as well
-as the hydrogen-bonding interaction *oxdna/hbond* between complementary pairs of nucleotides on
-opposite strands. Average sequence or sequence-dependent stacking and base-pairing strengths
-are supported :ref:`(Sulc) <Sulc1>`. Quasi-unique base-pairing between nucleotides can be achieved by using
-more complementary pairs of atom types like 5-8 and 6-7, 9-12 and 10-11, 13-16 and 14-15, etc.
-This prevents the hybridization of in principle complementary bases within Ntypes/4 bases
-up and down along the backbone.
+The *oxdna* pair styles compute the pairwise-additive parts of the oxDNA
+force field for coarse-grained modelling of DNA. The effective
+interaction between the nucleotides consists of potentials for the
+excluded volume interaction *oxdna/excv*, the stacking *oxdna/stk*,
+cross-stacking *oxdna/xstk* and coaxial stacking interaction
+*oxdna/coaxstk* as well as the hydrogen-bonding interaction
+*oxdna/hbond* between complementary pairs of nucleotides on opposite
+strands. Average sequence or sequence-dependent stacking and
+base-pairing strengths are supported :ref:`(Sulc) <Sulc1>`. Quasi-unique
+base-pairing between nucleotides can be achieved by using more
+complementary pairs of atom types like 5-8 and 6-7, 9-12 and 10-11,
+13-16 and 14-15, etc.  This prevents the hybridization of in principle
+complementary bases within Ntypes/4 bases up and down along the
+backbone.
 
-The exact functional form of the pair styles is rather complex.
-The individual potentials consist of products of modulation factors,
-which themselves are constructed from a number of more basic potentials
-(Morse, Lennard-Jones, harmonic angle and distance) as well as quadratic smoothing and modulation terms.
-We refer to :ref:`(Ouldridge-DPhil) <Ouldridge-DPhil1>` and :ref:`(Ouldridge) <Ouldridge1>`
-for a detailed description of the oxDNA force field.
+The exact functional form of the pair styles is rather complex.  The
+individual potentials consist of products of modulation factors, which
+themselves are constructed from a number of more basic potentials
+(Morse, Lennard-Jones, harmonic angle and distance) as well as quadratic
+smoothing and modulation terms.  We refer to :ref:`(Ouldridge-DPhil)
+<Ouldridge-DPhil1>` and :ref:`(Ouldridge) <Ouldridge1>` for a detailed
+description of the oxDNA force field.
 
 .. note::
 
-   These pair styles have to be used together with the related oxDNA bond style
-   *oxdna/fene* for the connectivity of the phosphate backbone (see also documentation of
-   :doc:`bond_style oxdna/fene <bond_oxdna>`). Most of the coefficients
-   in the above example have to be kept fixed and cannot be changed without reparameterizing the entire model.
-   Exceptions are the first four coefficients after *oxdna/stk* (seq=seqdep, T=0.1, xi=1.3448 and kappa=2.6568 in the above example)
-   and the first coefficient after *oxdna/hbond* (seq=seqdep in the above example).
-   When using a Langevin thermostat, e.g. through :doc:`fix langevin <fix_langevin>`
-   or :doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`
-   the temperature coefficients have to be matched to the one used in the fix.
+   These pair styles have to be used together with the related oxDNA
+   bond style *oxdna/fene* for the connectivity of the phosphate
+   backbone (see also documentation of :doc:`bond_style oxdna/fene
+   <bond_oxdna>`). Most of the coefficients in the above example have to
+   be kept fixed and cannot be changed without reparameterizing the
+   entire model.  Exceptions are the first four coefficients after
+   *oxdna/stk* (seq=seqdep, T=0.1, xi=1.3448 and kappa=2.6568 and
+   corresponding *real unit* equivalents in the above examples) and the
+   first coefficient after *oxdna/hbond* (seq=seqdep in the above
+   example).  When using a Langevin thermostat, e.g. through :doc:`fix
+   langevin <fix_langevin>` or :doc:`fix nve/dotc/langevin
+   <fix_nve_dotc_langevin>` the temperature coefficients have to be
+   matched to the one used in the fix.
 
 .. note::
 
-   These pair styles have to be used with the *atom_style hybrid bond ellipsoid oxdna*
-   (see documentation of :doc:`atom_style <atom_style>`). The *atom_style oxdna*
-   stores the 3'-to-5' polarity of the nucleotide strand, which is set through
-   the bond topology in the data file. The first (second) atom in a bond definition
-   is understood to point towards the 3'-end (5'-end) of the strand.
+   These pair styles have to be used with the *atom_style hybrid bond
+   ellipsoid oxdna* (see documentation of :doc:`atom_style
+   <atom_style>`). The *atom_style oxdna* stores the 3'-to-5' polarity
+   of the nucleotide strand, which is set through the bond topology in
+   the data file. The first (second) atom in a bond definition is
+   understood to point towards the 3'-end (5'-end) of the strand.
 
-Example input and data files for DNA duplexes can be found in examples/PACKAGES/cgdna/examples/oxDNA/ and /oxDNA2/.
-A simple python setup tool which creates single straight or helical DNA strands,
-DNA duplexes or arrays of DNA duplexes can be found in examples/PACKAGES/cgdna/util/.
+Example input and data files for DNA duplexes can be found in
+``examples/PACKAGES/cgdna/examples/oxDNA/`` and ``.../oxDNA2/``.  A
+simple python setup tool which creates single straight or helical DNA
+strands, DNA duplexes or arrays of DNA duplexes can be found in
+``examples/PACKAGES/cgdna/util/``.
 
 Please cite :ref:`(Henrich) <Henrich1>` in any publication that uses
-this implementation. An updated documentation that contains general information
-on the model, its implementation and performance as well as the structure of
-the data and input file can be found `here <PDF/CG-DNA.pdf>`_.
+this implementation. An updated documentation that contains general
+information on the model, its implementation and performance as well as
+the structure of the data and input file can be found `here
+<PDF/CG-DNA.pdf>`_.
 
 Please cite also the relevant oxDNA publications
 :ref:`(Ouldridge) <Ouldridge1>`,
@@ -115,6 +166,57 @@ and :ref:`(Sulc) <Sulc1>`.
 
 ----------
 
+Potential file reading
+""""""""""""""""""""""
+
+For each pair style above the first non-modifiable argument can be a
+filename, and if it is, no further arguments should be
+supplied. Therefore the following command:
+
+.. code-block:: LAMMPS
+
+   pair_coeff 1 4 oxdna/hbond   seqav oxdna_lj.cgdna
+
+will be interpreted as a request to read the corresponding hydrogen
+bonding potential parameters from the file with the given name. The file
+can define multiple potential parameters for both bonded and pair
+interactions, but for the example pair interaction above there must
+exist in the file a line of the form:
+
+.. code-block:: LAMMPS
+
+   1 4 hbond     <coefficients>
+
+If potential customization is required, the potential file reading can
+be mixed with the manual specification of the potential parameters. For
+example, the following command:
+
+.. code-block:: LAMMPS
+
+   pair_style hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk
+   pair_coeff * * oxdna/excv    oxdna_lj.cgdna
+   pair_coeff * * oxdna/stk     seqav 0.1 1.3448 2.6568 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
+   pair_coeff * * oxdna/hbond   seqav oxdna_lj.cgdna
+   pair_coeff 1 4 oxdna/hbond   seqav oxdna_lj.cgdna
+   pair_coeff 2 3 oxdna/hbond   seqav oxdna_lj.cgdna
+   pair_coeff * * oxdna/xstk    oxdna_lj.cgdna
+   pair_coeff * * oxdna/coaxstk 46.0 0.4 0.6 0.22 0.58 2.0 2.541592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 -0.65 2.0 -0.65
+
+will read the stacking and coaxial stacking potential parameters from
+the manual specification and all others from the potential file
+*oxdna_lj.cgdna*.
+
+There are sample potential files for each unit style in the
+``potentials`` directory of the LAMMPS distribution. The potential file
+unit system must align with the units defined via the :doc:`units
+<units>` command. For conversion between different *LJ* and *real* unit
+systems for oxDNA, the python tool *lj2real.py* located in the
+``examples/PACKAGES/cgdna/util/`` directory can be used. This tool
+assumes similar file structure to the examples found in
+``examples/PACKAGES/cgdna/examples/``.
+
+----------
+
 Restrictions
 """"""""""""
 

doc/src/pair_oxdna2.rst

@@ -41,14 +41,14 @@ Syntax
 
      *oxdna2/stk* args = seq T xi kappa 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
        seq = seqav (for average sequence stacking strength) or seqdep (for sequence-dependent stacking strength)
-       T = temperature (oxDNA units, 0.1 = 300 K)
-       xi = 1.3523 (temperature-independent coefficient in stacking strength)
-       kappa = 2.6717 (coefficient of linear temperature dependence in stacking strength)
+       T = temperature (LJ units: 0.1 = 300 K, real units: 300 = 300 K)
+       xi = 1.3523 (LJ units) or 8.06199211612242 (real units), temperature-independent coefficient in stacking strength
+       kappa = 2.6717 (LJ units) or 0.005309213 (real units), coefficient of linear temperature dependence in stacking strength
      *oxdna2/hbond* args = seq eps 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
        seq = seqav (for average sequence base-pairing strength) or seqdep (for sequence-dependent base-pairing strength)
-       eps = 1.0678 (between base pairs A-T and C-G) or 0 (all other pairs)
+       eps = 1.0678 (LJ units) or 6.36589157849259 (real units), average hydrogen bonding strength between A-T and C-G Watson-Crick base pairs, 0 between all other pairs
      *oxdna2/dh* args = T rhos qeff
-       T = temperature (oxDNA units, 0.1 = 300 K)
+       T = temperature (LJ units: 0.1 = 300 K, real units: 300 = 300 K)
        rhos = salt concentration (mole per litre)
        qeff = 0.815 (effective charge in elementary charges)
 
@@ -57,6 +57,7 @@ Examples
 
 .. code-block:: LAMMPS
 
+   # LJ units
    pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
    pair_coeff * * oxdna2/excv    2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
    pair_coeff * * oxdna2/stk     seqdep 0.1 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
@@ -67,61 +68,169 @@ Examples
    pair_coeff * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793
    pair_coeff * * oxdna2/dh      0.1 0.5 0.815
 
+   pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
+   pair_coeff * * oxdna2/excv    oxdna2_lj.cgdna
+   pair_coeff * * oxdna2/stk     seqdep 0.1 1.3523 2.6717 oxdna2_lj.cgdna
+   pair_coeff * * oxdna2/hbond   seqdep oxdna2_lj.cgdna
+   pair_coeff 1 4 oxdna2/hbond   seqdep oxdna2_lj.cgdna
+   pair_coeff 2 3 oxdna2/hbond   seqdep oxdna2_lj.cgdna
+   pair_coeff * * oxdna2/xstk    oxdna2_lj.cgdna
+   pair_coeff * * oxdna2/coaxstk oxdna2_lj.cgdna
+   pair_coeff * * oxdna2/dh      0.1 0.5 oxdna2_lj.cgdna
+
+   # Real units
+   pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
+   pair_coeff * * oxdna2/excv    11.92337812042065 5.9626 5.74965 11.92337812042065 4.38677 4.259 11.92337812042065 2.81094 2.72576
+   pair_coeff * * oxdna2/stk     seqdep 300.0 8.06199211612242 0.005309213 0.70439070204273 3.4072 7.6662 2.72576 6.3885 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 2.0 0.65 2.0 0.65
+   pair_coeff * * oxdna2/hbond   seqdep 0.0 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45
+   pair_coeff 1 4 oxdna2/hbond   seqdep 6.36589157849259 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45
+   pair_coeff 2 3 oxdna2/hbond   seqdep 6.36589157849259 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45
+   pair_coeff * * oxdna2/xstk    3.9029021145006 4.89785 5.74965 4.21641 5.57929 2.25 0.791592654 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0.0 0.65 1.7 0.875 0.68 1.7 0.875 0.68
+   pair_coeff * * oxdna2/coaxstk 4.80673207785863 3.4072 5.1108 1.87396 4.94044 2.0 2.891592653589793 0.65 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 40.0 3.116592653589793
+   pair_coeff * * oxdna2/dh      300.0 0.5 0.815
+
+   pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
+   pair_coeff * * oxdna2/excv    oxdna2_real.cgdna
+   pair_coeff * * oxdna2/stk     seqdep 300.0 8.06199211612242 0.005309213 oxdna2_real.cgdna
+   pair_coeff * * oxdna2/hbond   seqdep oxdna2_real.cgdna
+   pair_coeff 1 4 oxdna2/hbond   seqdep oxdna2_real.cgdna
+   pair_coeff 2 3 oxdna2/hbond   seqdep oxdna2_real.cgdna
+   pair_coeff * * oxdna2/xstk    oxdna2_real.cgdna
+   pair_coeff * * oxdna2/coaxstk oxdna2_real.cgdna
+   pair_coeff * * oxdna2/dh      300.0 0.5 oxdna2_real.cgdna
+
+.. note::
+
+   The coefficients in the above examples are provided in forms
+   compatible with both *units lj* and *units real* (see documentation
+   of :doc:`units <units>`).  These can also be read from a potential
+   file with correct unit style by specifying the name of the
+   file. Several potential files for each unit style are included in the
+   ``potentials`` directory of the LAMMPS distribution.
+
 Description
 """""""""""
 
-The *oxdna2* pair styles compute the pairwise-additive parts of the oxDNA force field
-for coarse-grained modelling of DNA. The effective interaction between the nucleotides consists of potentials for the
-excluded volume interaction *oxdna2/excv*, the stacking *oxdna2/stk*, cross-stacking *oxdna2/xstk*
-and coaxial stacking interaction *oxdna2/coaxstk*, electrostatic Debye-Hueckel interaction *oxdna2/dh*
-as well as the hydrogen-bonding interaction *oxdna2/hbond* between complementary pairs of nucleotides on
-opposite strands. Average sequence or sequence-dependent stacking and base-pairing strengths
-are supported :ref:`(Sulc) <Sulc2>`. Quasi-unique base-pairing between nucleotides can be achieved by using
-more complementary pairs of atom types like 5-8 and 6-7, 9-12 and 10-11, 13-16 and 14-15, etc.
-This prevents the hybridization of in principle complementary bases within Ntypes/4 bases
+The *oxdna2* pair styles compute the pairwise-additive parts of the
+oxDNA force field for coarse-grained modelling of DNA. The effective
+interaction between the nucleotides consists of potentials for the
+excluded volume interaction *oxdna2/excv*, the stacking *oxdna2/stk*,
+cross-stacking *oxdna2/xstk* and coaxial stacking interaction
+*oxdna2/coaxstk*, electrostatic Debye-Hueckel interaction *oxdna2/dh* as
+well as the hydrogen-bonding interaction *oxdna2/hbond* between
+complementary pairs of nucleotides on opposite strands. Average sequence
+or sequence-dependent stacking and base-pairing strengths are supported
+:ref:`(Sulc) <Sulc2>`. Quasi-unique base-pairing between nucleotides can
+be achieved by using more complementary pairs of atom types like 5-8 and
+6-7, 9-12 and 10-11, 13-16 and 14-15, etc.  This prevents the
+hybridization of in principle complementary bases within Ntypes/4 bases
 up and down along the backbone.
 
-The exact functional form of the pair styles is rather complex.
-The individual potentials consist of products of modulation factors,
-which themselves are constructed from a number of more basic potentials
-(Morse, Lennard-Jones, harmonic angle and distance) as well as quadratic smoothing and modulation terms.
-We refer to :ref:`(Snodin) <Snodin2>` and the original oxDNA publications :ref:`(Ouldridge-DPhil) <Ouldridge-DPhil2>`
-and  :ref:`(Ouldridge) <Ouldridge2>` for a detailed description of the oxDNA2 force field.
+The exact functional form of the pair styles is rather complex.  The
+individual potentials consist of products of modulation factors, which
+themselves are constructed from a number of more basic potentials
+(Morse, Lennard-Jones, harmonic angle and distance) as well as quadratic
+smoothing and modulation terms.  We refer to :ref:`(Snodin) <Snodin2>`
+and the original oxDNA publications :ref:`(Ouldridge-DPhil)
+<Ouldridge-DPhil2>` and :ref:`(Ouldridge) <Ouldridge2>` for a detailed
+description of the oxDNA2 force field.
 
 .. note::
 
-   These pair styles have to be used together with the related oxDNA2 bond style
-   *oxdna2/fene* for the connectivity of the phosphate backbone (see also documentation of
-   :doc:`bond_style oxdna2/fene <bond_oxdna>`). Most of the coefficients
-   in the above example have to be kept fixed and cannot be changed without reparameterizing the entire model.
-   Exceptions are the first four coefficients after *oxdna2/stk* (seq=seqdep, T=0.1, xi=1.3523 and kappa=2.6717 in the above example),
-   the first coefficient after *oxdna2/hbond* (seq=seqdep in the above example) and the three coefficients
-   after *oxdna2/dh* (T=0.1, rhos=0.5, qeff=0.815 in the above example). When using a Langevin thermostat
-   e.g. through :doc:`fix langevin <fix_langevin>` or :doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`
-   the temperature coefficients have to be matched to the one used in the fix.
+   These pair styles have to be used together with the related oxDNA2
+   bond style *oxdna2/fene* for the connectivity of the phosphate
+   backbone (see also documentation of :doc:`bond_style oxdna2/fene
+   <bond_oxdna>`). Most of the coefficients in the above example have to
+   be kept fixed and cannot be changed without reparameterizing the
+   entire model.  Exceptions are the first four coefficients after
+   *oxdna2/stk* (seq=seqdep, T=0.1, xi=1.3523 and kappa=2.6717 and
+   corresponding *real unit* equivalents in the above examples).  the
+   first coefficient after *oxdna2/hbond* (seq=seqdep in the above
+   example) and the three coefficients after *oxdna2/dh* (T=0.1,
+   rhos=0.5, qeff=0.815 in the above example). When using a Langevin
+   thermostat e.g. through :doc:`fix langevin <fix_langevin>` or
+   :doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>` the temperature
+   coefficients have to be matched to the one used in the fix.
 
 .. note::
 
-   These pair styles have to be used with the *atom_style hybrid bond ellipsoid oxdna*
-   (see documentation of :doc:`atom_style <atom_style>`). The *atom_style oxdna*
-   stores the 3'-to-5' polarity of the nucleotide strand, which is set through
-   the bond topology in the data file. The first (second) atom in a bond definition
-   is understood to point towards the 3'-end (5'-end) of the strand.
+   These pair styles have to be used with the *atom_style hybrid bond
+   ellipsoid oxdna* (see documentation of :doc:`atom_style
+   <atom_style>`). The *atom_style oxdna* stores the 3'-to-5' polarity
+   of the nucleotide strand, which is set through the bond topology in
+   the data file. The first (second) atom in a bond definition is
+   understood to point towards the 3'-end (5'-end) of the strand.
 
-Example input and data files for DNA duplexes can be found in examples/PACKAGES/cgdna/examples/oxDNA/ and /oxDNA2/.
-A simple python setup tool which creates single straight or helical DNA strands,
-DNA duplexes or arrays of DNA duplexes can be found in examples/PACKAGES/cgdna/util/.
+Example input and data files for DNA duplexes can be found in
+``examples/PACKAGES/cgdna/examples/oxDNA/`` and ``.../oxDNA2/``.  A
+simple python setup tool which creates single straight or helical DNA
+strands, DNA duplexes or arrays of DNA duplexes can be found in
+``examples/PACKAGES/cgdna/util/``.
 
 Please cite :ref:`(Henrich) <Henrich2>` in any publication that uses
-this implementation. An updated documentation that contains general information
-on the model, its implementation and performance as well as the structure of
-the data and input file can be found `here <PDF/CG-DNA.pdf>`_.
+this implementation. An updated documentation that contains general
+information on the model, its implementation and performance as well as
+the structure of the data and input file can be found `here
+<PDF/CG-DNA.pdf>`_.
 
 Please cite also the relevant oxDNA2 publications
 :ref:`(Snodin) <Snodin2>` and :ref:`(Sulc) <Sulc2>`.
 
 ----------
 
+Potential file reading
+""""""""""""""""""""""
+
+For each pair style above the first non-modifiable argument can be a
+filename (with exception of Debye-Hueckel, for which the effective
+charge argument can be a filename), and if it is, no further arguments
+should be supplied.  Therefore the following command:
+
+.. code-block:: LAMMPS
+
+   pair_coeff 1 4 oxdna2/hbond   seqdep oxdna_real.cgdna
+
+will be interpreted as a request to read the corresponding hydrogen
+bonding potential parameters from the file with the given name.  The
+file can define multiple potential parameters for both bonded and pair
+interactions, but for the example pair interaction above there must
+exist in the file a line of the form:
+
+.. code-block:: LAMMPS
+
+  1 4 hbond     <coefficients>
+
+If potential customization is required, the potential file reading can
+be mixed with the manual specification of the potential parameters. For
+example, the following command:
+
+.. code-block:: LAMMPS
+
+   pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
+   pair_coeff * * oxdna2/excv    2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
+   pair_coeff * * oxdna2/stk     seqdep 0.1 1.3523 2.6717 oxdna2_lj.cgdna
+   pair_coeff * * oxdna2/hbond   seqdep oxdna2_lj.cgdna
+   pair_coeff 1 4 oxdna2/hbond   seqdep oxdna2_lj.cgdna
+   pair_coeff 2 3 oxdna2/hbond   seqdep oxdna2_lj.cgdna
+   pair_coeff * * oxdna2/xstk    oxdna2_lj.cgdna
+   pair_coeff * * oxdna2/coaxstk oxdna2_lj.cgdna
+   pair_coeff * * oxdna2/dh      0.1 0.5 0.815
+
+will read the excluded volume and Debye-Hueckel effective charge *qeff*
+parameters from the manual specification and all others from the
+potential file *oxdna2_lj.cgdna*.
+
+There are sample potential files for each unit style in the ``potentials``
+directory of the LAMMPS distribution. The potential file unit system
+must align with the units defined via the :doc:`units <units>`
+command. For conversion between different *LJ* and *real* unit systems
+for oxDNA, the python tool *lj2real.py* located in the
+``examples/PACKAGES/cgdna/util/`` directory can be used. This tool assumes
+similar file structure to the examples found in
+``examples/PACKAGES/cgdna/examples/``.
+
+----------
+
 Restrictions
 """"""""""""
 

doc/src/pair_oxrna2.rst

@@ -41,14 +41,14 @@ Syntax
 
      *oxrna2/stk* args = seq T xi kappa 6.0 0.43 0.93 0.35 0.78 0.9 0 0.95 0.9 0 0.95 1.3 0 0.8 1.3 0 0.8 2.0 0.65 2.0 0.65
        seq = seqav (for average sequence stacking strength) or seqdep (for sequence-dependent stacking strength)
-       T = temperature (oxDNA units, 0.1 = 300 K)
-       xi = 1.40206 (temperature-independent coefficient in stacking strength)
-       kappa = 2.77 (coefficient of linear temperature dependence in stacking strength)
+       T = temperature (LJ units: 0.1 = 300 K, real units: 300 = 300 K)
+       xi = 1.40206 (LJ units) or 8.35864576375849 (real units), temperature-independent coefficient in stacking strength
+       kappa = 2.77 (LJ units) or 0.005504556 (real units), coefficient of linear temperature dependence in stacking strength
      *oxrna2/hbond* args = seq eps 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
        seq = seqav (for average sequence base-pairing strength) or seqdep (for sequence-dependent base-pairing strength)
-       eps = 0.870439 (between base pairs A-T, C-G and G-T) or 0 (all other pairs)
+       eps = 0.870439 (LJ units) or 5.18928666388042 (real units), average hydrogen bonding strength between A-U and C-G Watson-Crick and G-U wobble base pairs, 0 between all other pairs
      *oxrna2/dh* args = T rhos qeff
-       T = temperature (oxDNA units, 0.1 = 300 K)
+       T = temperature (LJ units: 0.1 = 300 K, real units: 300 = 300 K)
        rhos = salt concentration (mole per litre)
        qeff = 1.02455 (effective charge in elementary charges)
 
@@ -57,6 +57,7 @@ Examples
 
 .. code-block:: LAMMPS
 
+   # LJ units
    pair_style hybrid/overlay oxrna2/excv oxrna2/stk oxrna2/hbond oxrna2/xstk oxrna2/coaxstk oxrna2/dh
    pair_coeff * * oxrna2/excv    2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
    pair_coeff * * oxrna2/stk     seqdep 0.1 1.40206 2.77 6.0 0.43 0.93 0.35 0.78 0.9 0 0.95 0.9 0 0.95 1.3 0 0.8 1.3 0 0.8 2.0 0.65 2.0 0.65
@@ -68,58 +69,168 @@ Examples
    pair_coeff * * oxrna2/coaxstk 80 0.5 0.6 0.42 0.58 2.0 2.592 0.65 1.3 0.151 0.8 0.9 0.685 0.95 0.9 0.685 0.95 2.0 -0.65 2.0 -0.65
    pair_coeff * * oxrna2/dh      0.1 0.5 1.02455
 
+   pair_style hybrid/overlay oxrna2/excv oxrna2/stk oxrna2/hbond oxrna2/xstk oxrna2/coaxstk oxrna2/dh
+   pair_coeff * * oxrna2/excv    oxrna2_lj.cgdna
+   pair_coeff * * oxrna2/stk     seqdep 0.1 1.40206 2.77 oxrna2_lj.cgdna
+   pair_coeff * * oxrna2/hbond   seqdep oxrna2_lj.cgdna
+   pair_coeff 1 4 oxrna2/hbond   seqdep oxrna2_lj.cgdna
+   pair_coeff 2 3 oxrna2/hbond   seqdep oxrna2_lj.cgdna
+   pair_coeff 3 4 oxrna2/hbond   seqdep oxrna2_lj.cgdna
+   pair_coeff * * oxrna2/xstk    oxrna2_lj.cgdna
+   pair_coeff * * oxrna2/coaxstk oxrna2_lj.cgdna
+   pair_coeff * * oxrna2/dh      0.1 0.5 oxrna2_lj.cgdna
+
+   # Real units
+   pair_style hybrid/overlay oxrna2/excv oxrna2/stk oxrna2/hbond oxrna2/xstk oxrna2/coaxstk oxrna2/dh
+   pair_coeff * * oxrna2/excv    11.92337812042065 5.9626 5.74965 11.92337812042065 4.38677 4.259 11.92337812042065 2.81094 2.72576
+   pair_coeff * * oxrna2/stk     seqdep 300.0 8.35864576375849 0.005504556 0.70439070204273 3.66274 7.92174 2.9813 6.64404 0.9 0.0 0.95 0.9 0.0 0.95 1.3 0.0 0.8 1.3 0.0 0.8 2.0 0.65 2.0 0.65
+   pair_coeff * * oxrna2/hbond   seqdep 0.0 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
+   pair_coeff 1 4 oxrna2/hbond   seqdep 5.18928666388042 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
+   pair_coeff 2 3 oxrna2/hbond   seqdep 5.18928666388042 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
+   pair_coeff 3 4 oxrna2/hbond   seqdep 5.18928666388042 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
+   pair_coeff * * oxrna2/xstk    4.92690859644113 4.259 5.1108 3.57756 4.94044 2.25 0.505 0.58 1.7 1.266 0.68 1.7 1.266 0.68 1.7 0.309 0.68 1.7 0.309 0.68
+   pair_coeff * * oxrna2/coaxstk 6.57330882442206 4.259 5.1108 3.57756 4.94044 2.0 2.592 0.65 1.3 0.151 0.8 0.9 0.685 0.95 0.9 0.685 0.95 2.0 -0.65 2.0 -0.65
+   pair_coeff * * oxrna2/dh      300.0 0.5 1.02455
+
+   pair_style hybrid/overlay oxrna2/excv oxrna2/stk oxrna2/hbond oxrna2/xstk oxrna2/coaxstk oxrna2/dh
+   pair_coeff * * oxrna2/excv    oxrna2_real.cgdna
+   pair_coeff * * oxrna2/stk     seqdep 300.0 8.35864576375849 0.005504556 oxrna2_real.cgdna
+   pair_coeff * * oxrna2/hbond   seqdep oxrna2_real.cgdna
+   pair_coeff 1 4 oxrna2/hbond   seqdep oxrna2_real.cgdna
+   pair_coeff 2 3 oxrna2/hbond   seqdep oxrna2_real.cgdna
+   pair_coeff 3 4 oxrna2/hbond   seqdep oxrna2_real.cgdna
+   pair_coeff * * oxrna2/xstk    oxrna2_real.cgdna
+   pair_coeff * * oxrna2/coaxstk oxrna2_real.cgdna
+   pair_coeff * * oxrna2/dh      300.0 0.5 oxrna2_real.cgdna
+
+.. note::
+
+   The coefficients in the above examples are provided in forms
+   compatible with both *units lj* and *units real* (see documentation
+   of :doc:`units <units>`).  These can also be read from a potential
+   file with correct unit style by specifying the name of the
+   file. Several potential files for each unit style are included in the
+   ``potentials`` directory of the LAMMPS distribution.
+
 Description
 """""""""""
 
-The *oxrna2* pair styles compute the pairwise-additive parts of the oxDNA force field
-for coarse-grained modelling of DNA. The effective interaction between the nucleotides consists of potentials for the
-excluded volume interaction *oxrna2/excv*, the stacking *oxrna2/stk*, cross-stacking *oxrna2/xstk*
-and coaxial stacking interaction *oxrna2/coaxstk*, electrostatic Debye-Hueckel interaction *oxrna2/dh*
-as well as the hydrogen-bonding interaction *oxrna2/hbond* between complementary pairs of nucleotides on
-opposite strands. Average sequence or sequence-dependent stacking and base-pairing strengths
-are supported :ref:`(Sulc2) <Sulc32>`. Quasi-unique base-pairing between nucleotides can be achieved by using
-more complementary pairs of atom types like 5-8 and 6-7, 9-12 and 10-11, 13-16 and 14-15, etc.
-This prevents the hybridization of in principle complementary bases within Ntypes/4 bases
+The *oxrna2* pair styles compute the pairwise-additive parts of the
+oxDNA force field for coarse-grained modelling of RNA. The effective
+interaction between the nucleotides consists of potentials for the
+excluded volume interaction *oxrna2/excv*, the stacking *oxrna2/stk*,
+cross-stacking *oxrna2/xstk* and coaxial stacking interaction
+*oxrna2/coaxstk*, electrostatic Debye-Hueckel interaction *oxrna2/dh* as
+well as the hydrogen-bonding interaction *oxrna2/hbond* between
+complementary pairs of nucleotides on opposite strands. Average sequence
+or sequence-dependent stacking and base-pairing strengths are supported
+:ref:`(Sulc2) <Sulc32>`. Quasi-unique base-pairing between nucleotides
+can be achieved by using more complementary pairs of atom types like 5-8
+and 6-7, 9-12 and 10-11, 13-16 and 14-15, etc.  This prevents the
+hybridization of in principle complementary bases within Ntypes/4 bases
 up and down along the backbone.
 
-The exact functional form of the pair styles is rather complex.
-The individual potentials consist of products of modulation factors,
-which themselves are constructed from a number of more basic potentials
-(Morse, Lennard-Jones, harmonic angle and distance) as well as quadratic smoothing and modulation terms.
-We refer to :ref:`(Sulc1) <Sulc31>` and the original oxDNA publications :ref:`(Ouldridge-DPhil) <Ouldridge-DPhil3>`
-and  :ref:`(Ouldridge) <Ouldridge3>` for a detailed description of the oxRNA2 force field.
+The exact functional form of the pair styles is rather complex.  The
+individual potentials consist of products of modulation factors, which
+themselves are constructed from a number of more basic potentials
+(Morse, Lennard-Jones, harmonic angle and distance) as well as quadratic
+smoothing and modulation terms.  We refer to :ref:`(Sulc1) <Sulc31>` and
+the original oxDNA publications :ref:`(Ouldridge-DPhil)
+<Ouldridge-DPhil3>` and :ref:`(Ouldridge) <Ouldridge3>` for a detailed
+description of the oxRNA2 force field.
 
 .. note::
 
-   These pair styles have to be used together with the related oxDNA2 bond style
-   *oxrna2/fene* for the connectivity of the phosphate backbone (see also documentation of
-   :doc:`bond_style oxrna2/fene <bond_oxdna>`). Most of the coefficients
-   in the above example have to be kept fixed and cannot be changed without reparameterizing the entire model.
-   Exceptions are the first four coefficients after *oxrna2/stk* (seq=seqdep, T=0.1, xi=1.40206 and kappa=2.77 in the above example),
-   the first coefficient after *oxrna2/hbond* (seq=seqdep in the above example) and the three coefficients
-   after *oxrna2/dh* (T=0.1, rhos=0.5, qeff=1.02455 in the above example). When using a Langevin thermostat
-   e.g. through :doc:`fix langevin <fix_langevin>` or :doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`
-   the temperature coefficients have to be matched to the one used in the fix.
+   These pair styles have to be used together with the related oxDNA2
+   bond style *oxrna2/fene* for the connectivity of the phosphate
+   backbone (see also documentation of :doc:`bond_style oxrna2/fene
+   <bond_oxdna>`). Most of the coefficients in the above example have to
+   be kept fixed and cannot be changed without reparameterizing the
+   entire model.  Exceptions are the first four coefficients after
+   *oxrna2/stk* (seq=seqdep, T=0.1, xi=1.40206 and kappa=2.77 and
+   corresponding *real unit* equivalents in the above examples), the
+   first coefficient after *oxrna2/hbond* (seq=seqdep in the above
+   example) and the three coefficients after *oxrna2/dh* (T=0.1,
+   rhos=0.5, qeff=1.02455 in the above example). When using a Langevin
+   thermostat e.g. through :doc:`fix langevin <fix_langevin>` or
+   :doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>` the temperature
+   coefficients have to be matched to the one used in the fix.
 
 .. note::
 
-   These pair styles have to be used with the *atom_style hybrid bond ellipsoid oxdna*
-   (see documentation of :doc:`atom_style <atom_style>`). The *atom_style oxdna*
-   stores the 3'-to-5' polarity of the nucleotide strand, which is set through
-   the bond topology in the data file. The first (second) atom in a bond definition
-   is understood to point towards the 3'-end (5'-end) of the strand.
+   These pair styles have to be used with the *atom_style hybrid bond
+   ellipsoid oxdna* (see documentation of :doc:`atom_style
+   <atom_style>`). The *atom_style oxdna* stores the 3'-to-5' polarity
+   of the nucleotide strand, which is set through the bond topology in
+   the data file. The first (second) atom in a bond definition is
+   understood to point towards the 3'-end (5'-end) of the strand.
 
-Example input and data files for DNA duplexes can be found in examples/PACKAGES/cgdna/examples/oxDNA/ and /oxDNA2/.
-A simple python setup tool which creates single straight or helical DNA strands,
-DNA duplexes or arrays of DNA duplexes can be found in examples/PACKAGES/cgdna/util/.
+Example input and data files for DNA duplexes can be found in
+``examples/PACKAGES/cgdna/examples/oxDNA/`` and ``.../oxDNA2/``.  A simple python
+setup tool which creates single straight or helical DNA strands, DNA
+duplexes or arrays of DNA duplexes can be found in
+``examples/PACKAGES/cgdna/util/``.
 
 Please cite :ref:`(Henrich) <Henrich3>` in any publication that uses
-this implementation.  The article contains general information
-on the model, its implementation and performance as well as the structure of
-the data and input file. The preprint version of the article can be found
-`here <PDF/CG-DNA.pdf>`_.
-Please cite also the relevant oxRNA2 publications
-:ref:`(Sulc1) <Sulc31>` and :ref:`(Sulc2) <Sulc32>`.
+this implementation.  The article contains general information on the
+model, its implementation and performance as well as the structure of
+the data and input file. The preprint version of the article can be
+found `here <PDF/CG-DNA.pdf>`_.  Please cite also the relevant oxRNA2
+publications :ref:`(Sulc1) <Sulc31>` and :ref:`(Sulc2) <Sulc32>`.
+
+----------
+
+Potential file reading
+""""""""""""""""""""""
+
+For each pair style above the first non-modifiable argument can be a
+filename (with exception of Debye-Hueckel, for which the effective
+charge argument can be a filename), and if it is, no further arguments
+should be supplied.  Therefore the following command:
+
+.. code-block:: LAMMPS
+
+   pair_coeff 3 4 oxrna2/hbond   seqdep oxrna2_lj.cgdna
+
+will be interpreted as a request to read the corresponding hydrogen
+bonding potential parameters from the file with the given name.  The
+file can define multiple potential parameters for both bonded and pair
+interactions, but for the example pair interaction above there must
+exist in the file a line of the form:
+
+.. code-block:: LAMMPS
+
+  3 4 hbond     <coefficients>
+
+If potential customization is required, the potential file reading can
+be mixed with the manual specification of the potential parameters. For
+example, the following command:
+
+.. code-block:: LAMMPS
+
+   pair_style hybrid/overlay oxrna2/excv oxrna2/stk oxrna2/hbond oxrna2/xstk oxrna2/coaxstk oxrna2/dh
+   pair_coeff * * oxrna2/excv    2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
+   pair_coeff * * oxrna2/stk     seqdep 0.1 1.40206 2.77 oxrna2_lj.cgdna
+   pair_coeff * * oxrna2/hbond   seqdep oxrna2_lj.cgdna
+   pair_coeff 1 4 oxrna2/hbond   seqdep oxrna2_lj.cgdna
+   pair_coeff 2 3 oxrna2/hbond   seqdep oxrna2_lj.cgdna
+   pair_coeff 3 4 oxrna2/hbond   seqdep oxrna2_lj.cgdna
+   pair_coeff * * oxrna2/xstk    oxrna2_lj.cgdna
+   pair_coeff * * oxrna2/coaxstk oxrna2_lj.cgdna
+   pair_coeff * * oxrna2/dh      0.1 0.5 1.02455
+
+will read the excluded volume and Debye-Hueckel effective charge *qeff*
+parameters from the manual specification and all others from the
+potential file *oxrna2_lj.cgdna*.
+
+There are sample potential files for each unit style in the
+``potentials`` directory of the LAMMPS distribution. The potential file
+unit system must align with the units defined via the :doc:`units
+<units>` command. For conversion between different *LJ* and *real* unit
+systems for oxDNA, the python tool *lj2real.py* located in the
+``examples/PACKAGES/cgdna/util/`` directory can be used. This tool
+assumes similar file structure to the examples found in
+``examples/PACKAGES/cgdna/examples/``.
 
 ----------
 

doc/src/pair_pod.rst

@@ -1,8 +1,11 @@
 .. index:: pair_style pod
+.. index:: pair_style pod/kk
 
 pair_style pod command
 ========================
 
+Accelerator Variants: *pod/kk*
+
 Syntax
 """"""
 
@@ -24,29 +27,33 @@ Description
 .. versionadded:: 22Dec2022
 
 Pair style *pod* defines the proper orthogonal descriptor (POD)
-potential :ref:`(Nguyen) <Nguyen20221>`.  The mathematical definition of
-the POD potential is described from :doc:`fitpod <fitpod_command>`, which is
-used to fit the POD potential to *ab initio* energy and force data.
+potential :ref:`(Nguyen and Rohskopf) <Nguyen20222b>`,
+:ref:`(Nguyen2023) <Nguyen20232b>`, :ref:`(Nguyen2024) <Nguyen20242b>`,
+and :ref:`(Nguyen and Sema) <Nguyen20243b>`.  The :doc:`fitpod
+<fitpod_command>` is used to fit the POD potential.
 
 Only a single pair_coeff command is used with the *pod* style which
-specifies a POD parameter file followed by a coefficient file.
+specifies a POD parameter file followed by a coefficient file, a
+projection matrix file, and a centroid file.
 
-The coefficient file (``Ta_coefficients.pod``) contains coefficients for the
-POD potential. The top of the coefficient file can contain any number of
-blank and comment lines (start with #), but follows a strict format
-after that. The first non-blank non-comment line must contain:
+The POD parameter file (``Ta_param.pod``) can contain blank and comment
+lines (start with #) anywhere. Each non-blank non-comment line must
+contain one keyword/value pair. See :doc:`fitpod <fitpod_command>` for
+the description of all the keywords that can be assigned in the
+parameter file.
 
-* POD_coefficients: *ncoeff*
+The coefficient file (``Ta_coefficients.pod``) contains coefficients for
+the POD potential. The top of the coefficient file can contain any
+number of blank and comment lines (start with #), but follows a strict
+format after that. The first non-blank non-comment line must contain:
 
-This is followed by *ncoeff* coefficients, one per line. The coefficient
+* model_coefficients: *ncoeff* *nproj* *ncentroid*
+
+This is followed by *ncoeff* coefficients, *nproj* projection matrix entries,
+and *ncentroid* centroid coordinates, one per line. The coefficient
 file is generated after training the POD potential using :doc:`fitpod
 <fitpod_command>`.
 
-The POD parameter file (``Ta_param.pod``) can contain blank and comment lines
-(start with #) anywhere. Each non-blank non-comment line must contain
-one keyword/value pair. See :doc:`fitpod <fitpod_command>` for the description
-of all the keywords that can be assigned in the parameter file.
-
 As an example, if a LAMMPS indium phosphide simulation has 4 atoms
 types, with the first two being indium and the third and fourth being
 phophorous, the pair_coeff command would look like this:
@@ -67,7 +74,33 @@ the *hybrid* pair style.  The NULL values are placeholders for atom
 types that will be used with other potentials.
 
 Examples about training and using POD potentials are found in the
-directory lammps/examples/PACKAGES/pod.
+directory lammps/examples/PACKAGES/pod and the Github repo https://github.com/cesmix-mit/pod-examples.
+
+----------
+
+Mixing, shift, table, tail correction, restart, rRESPA info
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+For atom type pairs I,J and I != J, where types I and J correspond to
+two different element types, mixing is performed by LAMMPS with
+user-specifiable parameters as described above.  You never need to
+specify a pair_coeff command with I != J arguments for this style.
+
+This pair style does not support the :doc:`pair_modify <pair_modify>`
+shift, table, and tail options.
+
+This pair style does not write its information to :doc:`binary restart
+files <restart>`, since it is stored in potential files.  Thus, you need
+to re-specify the pair_style and pair_coeff commands in an input script
+that reads a restart file.
+
+This pair style can only be used via the *pair* keyword of the
+:doc:`run_style respa <run_style>` command.  It does not support the
+*inner*, *middle*, *outer* keywords.
+
+----------
+
+.. include:: accel_styles.rst
 
 ----------
 
@@ -78,12 +111,14 @@ This style is part of the ML-POD package.  It is only enabled if LAMMPS
 was built with that package. See the :doc:`Build package
 <Build_package>` page for more info.
 
-This pair style does not compute per-atom energies and per-atom stresses.
-
 Related commands
 """"""""""""""""
 
 :doc:`fitpod <fitpod_command>`,
+:doc:`compute pod/atom <compute_pod_atom>`,
+:doc:`compute podd/atom <compute_pod_atom>`,
+:doc:`compute pod/local <compute_pod_atom>`,
+:doc:`compute pod/global <compute_pod_atom>`
 
 Default
 """""""
@@ -92,6 +127,20 @@ none
 
 ----------
 
-.. _Nguyen20221:
+.. _Nguyen20222b:
+
+**(Nguyen and Rohskopf)** Nguyen and Rohskopf,  Journal of Computational Physics, 480, 112030, (2023).
+
+.. _Nguyen20232b:
+
+**(Nguyen2023)** Nguyen, Physical Review B, 107(14), 144103, (2023).
+
+.. _Nguyen20242b:
+
+**(Nguyen2024)** Nguyen, Journal of Computational Physics, 113102, (2024).
+
+.. _Nguyen20243b:
+
+**(Nguyen and Sema)** Nguyen and Sema, https://arxiv.org/abs/2405.00306, (2024).
+
 
-**(Nguyen)** Nguyen and Rohskopf, arXiv preprint arXiv:2209.02362 (2022).

doc/src/pair_soft.rst

@@ -1,11 +1,12 @@
 .. index:: pair_style soft
 .. index:: pair_style soft/gpu
+.. index:: pair_style soft/kk
 .. index:: pair_style soft/omp
 
 pair_style soft command
 =======================
 
-Accelerator Variants: *soft/gpu*, *soft/omp*
+Accelerator Variants: *soft/gpu*, *soft/kk*, *soft/omp*
 
 Syntax
 """"""

doc/src/pair_style.rst

@@ -108,6 +108,7 @@ accelerated styles exist.
 
 * :doc:`none <pair_none>` - turn off pairwise interactions
 * :doc:`hybrid <pair_hybrid>` - multiple styles of pairwise interactions
+* :doc:`hybrid/molecular <pair_hybrid>` - different pair styles for intra- and inter-molecular interactions
 * :doc:`hybrid/overlay <pair_hybrid>` - multiple styles of superposed pairwise interactions
 * :doc:`hybrid/scaled <pair_hybrid>` - multiple styles of scaled superposed pairwise interactions
 * :doc:`zero <pair_zero>` - neighbor list but no interactions
@@ -171,6 +172,7 @@ accelerated styles exist.
 * :doc:`coul/wolf <pair_coul>` - Coulomb via Wolf potential
 * :doc:`coul/wolf/cs <pair_cs>` - Coulomb via Wolf potential with core/shell adjustments
 * :doc:`dpd <pair_dpd>` - dissipative particle dynamics (DPD)
+* :doc:`dpd/coul/slater/long <pair_dpd_coul_slater_long>` - dissipative particle dynamics (DPD) with electrostatic interactions
 * :doc:`dpd/ext <pair_dpd_ext>` - generalized force field for DPD
 * :doc:`dpd/ext/tstat <pair_dpd_ext>` - pairwise DPD thermostatting  with generalized force field
 * :doc:`dpd/fdt <pair_dpd_fdt>` - DPD for constant temperature and pressure
@@ -382,6 +384,7 @@ accelerated styles exist.
 * :doc:`tracker <pair_tracker>` - monitor information about pairwise interactions
 * :doc:`tri/lj <pair_tri_lj>` - LJ potential between triangles
 * :doc:`ufm <pair_ufm>` -
+* :doc:`uf3 <pair_uf3>` - UF3 machine-learning potential
 * :doc:`vashishta <pair_vashishta>` - Vashishta 2-body and 3-body potential
 * :doc:`vashishta/table <pair_vashishta>` -
 * :doc:`wf/cut <pair_wf_cut>` - Wang-Frenkel Potential for short-ranged interactions

doc/src/pair_uf3.rst

@@ -0,0 +1,213 @@
+.. index:: pair_style uf3
+.. index:: pair_style uf3/kk
+
+pair_style uf3 command
+======================
+
+Accelerator Variants: *uf3/kk*
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+    pair_style style BodyFlag
+
+* style = *uf3* or *uf3/kk*
+
+  .. parsed-literal::
+
+       BodyFlag = Indicates whether to calculate only 2-body or 2 and 3-body interactions. Possible values: 2 or 3
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+    pair_style uf3 3
+    pair_coeff * * Nb.uf3 Nb
+
+    pair_style uf3 2
+    pair_coeff * * NbSn.uf3 Nb Sn
+
+    pair_style uf3 3
+    pair_coeff * * NbSn.uf3 Nb Sn
+
+Description
+"""""""""""
+
+.. versionadded:: 27June2024
+
+The *uf3* style computes the :ref:`Ultra-Fast Force Fields (UF3)
+<Xie23>` potential, a machine-learning interatomic potential. In UF3,
+the total energy of the system is defined via two- and three-body
+interactions:
+
+.. math::
+
+   E & = \sum_{i,j} V_2(r_{ij}) + \sum_{i,j,k} V_3 (r_{ij},r_{ik},r_{jk}) \\
+   V_2(r_{ij}) & = \sum_{n=0}^N c_n B_n(r_{ij}) \\
+   V_3 (r_{ij},r_{ik},r_{jk}) & = \sum_{l=0}^{N_l} \sum_{m=0}^{N_m} \sum_{n=0}^{N_n} c_{l,m,n} B_l(r_{ij}) B_m(r_{ik}) B_n(r_{jk})
+
+where :math:`V_2(r_{ij})` and :math:`V_3 (r_{ij},r_{ik},r_{jk})` are the
+two- and three-body interactions, respectively. For the two-body the
+summation is over all neighbors J and for the three-body the summation
+is over all neighbors J and K of atom I within a cutoff distance
+determined from the potential files. :math:`B_n(r_{ij})` are the cubic
+b-spline basis, :math:`c_n` and :math:`c_{l,m,n}` are the machine-learned
+interaction parameters and :math:`N`, :math:`N_l`, :math:`N_m`, and
+:math:`N_n` denote the number of basis functions per spline or tensor
+spline dimension.
+
+With *uf3* style only a single pair_coeff command is used to indicate the
+UF3 LAMMPS potential file containing all the two- and three-body interactions
+followed by N additional arguments specifying the mapping of UF3 elements to
+LAMMPS atom types, where N is the number of LAMMPS atom types:
+
+* UF3 LAMMPS potential file
+* N elements names = mapping of UF3 elements to atom types
+
+As an example, if a LAMMPS simulation contains 2 atom types (elements
+'A' and 'B'), the pair_coeff command will be:
+
+.. code-block:: LAMMPS
+
+   pair_style uf3 3
+   pair_coeff * * AB.uf3 A B
+
+The AB.uf3 file should contain all two-body (A-A, A-B, B-B) and three-body
+(A-A-A, A-A-B, A-B-B, B-A-A, B-A-B, B-B-B).
+
+If a value of "2" is specified in the :code:`pair_style uf3` command,
+only the two-body potentials are needed. For 3-body interaction the
+first atom type is the central atom. We recommend using the
+:code:`generate_uf3_lammps_pots.py` script (found `here
+<https://github.com/uf3/uf3/tree/develop/lammps_plugin/scripts>`_) for
+generating the UF3 LAMMPS potential file from the UF3 JSON potentials.
+
+----------
+
+UF3 LAMMPS potential file in the *potentials* directory of the LAMMPS
+distribution have a ".uf3" suffix. The interaction block in UF3 LAMMPS potential
+file should start with :code:`#UF3 POT` and end with :code:`#` characters.
+Following shows the format of a generic 2-body and 3-body potential block in
+UF3 LAMMPS potential file-
+
+.. code-block:: LAMMPS
+
+   #UF3 POT UNITS: units DATE: POT_GEN_DATE AUTHOR: AUTHOR_NAME CITATION: CITE
+   2B ELEMENT1 ELEMENT2 LEADING_TRIM TRAILING_TRIM
+   Rij_CUTOFF NUM_OF_KNOTS
+   BSPLINE_KNOTS
+   NUM_OF_COEFF
+   COEFF
+   #
+   #UF3 POT UNITS: units DATE: POT_GEN_DATE AUTHOR: AUTHOR_NAME CITATION: CITE
+   3B ELEMENT1 ELEMENT2 ELEMENT3 LEADING_TRIM TRAILING_TRIM
+   Rjk_CUTOFF Rik_CUTOFF Rij_CUTOFF NUM_OF_KNOTS_JK NUM_OF_KNOTS_IK NUM_OF_KNOTS_IJ
+   BSPLINE_KNOTS_FOR_JK
+   BSPLINE_KNOTS_FOR_IK
+   BSPLINE_KNOTS_FOR_IJ
+   SHAPE_OF_COEFF_MATRIX[I][J][K]
+   COEFF_MATRIX[0][0][K]
+   COEFF_MATRIX[0][1][K]
+   COEFF_MATRIX[0][2][K]
+   .
+   .
+   .
+   COEFF_MATRIX[1][0][K]
+   COEFF_MATRIX[1][1][K]
+   COEFF_MATRIX[1][2][K]
+   .
+   .
+   .
+   #
+
+The second line indicates whether the block contains data for 2-body
+(:code:`2B`) or 3-body (:code:`3B`) interaction. This is followed by element
+combination interaction, :code:`LEADING_TRIM` and :code:`TRAILING_TRIM`
+number on the same line. The current implementation is only tested for
+:code:`LEADING_TRIM=0` and :code:`TRAILING_TRIM=3`.
+If other values are used LAMMPS is terminated after issuing an error message.
+The :code:`Rij_CUTOFF` sets the 2-body cutoff for the interaction described
+by the potential block. :code:`NUM_OF_KNOTS` is the number of knots
+(or the length of the knot vector) present on the very next line. The
+:code:`BSPLINE_KNOTS` line should contain all the knots in ascending order.
+:code:`NUM_OF_COEFF` is the number of coefficients in the :code:`COEFF` line.
+All the numbers in the BSPLINE_KNOTS and COEFF line should be space-separated.
+Similar to the 2-body potential block, the third line sets the cutoffs and
+length of the knots. The cutoff distance between atom-type I and J is
+:code:`Rij_CUTOFF`, atom-type I and K is :code:`Rik_CUTOFF` and between
+J and K is :code:`Rjk_CUTOFF`.
+
+.. note::
+
+   The current implementation only works for UF3 potentials with cutoff
+   distances for 3-body interactions that follows
+   :code:`2Rij_CUTOFF=2Rik_CUTOFF=Rjk_CUTOFF` relation.
+
+The :code:`BSPLINE_KNOTS_FOR_JK`, :code:`BSPLINE_KNOTS_FOR_IK`, and
+:code:`BSPLINE_KNOTS_FOR_IJ` lines (note the order) contain the knots in
+increasing order for atoms J and K, I and K, and atoms I and J
+respectively. The number of knots is defined by the
+:code:`NUM_OF_KNOTS_*` characters in the previous line.  The shape of
+the coefficient matrix is defined on the
+:code:`SHAPE_OF_COEFF_MATRIX[I][J][K]` line followed by the columns of
+the coefficient matrix, one per line, as shown above. For example, if
+the coefficient matrix has the shape of 8x8x13, then
+:code:`SHAPE_OF_COEFF_MATRIX[I][J][K]` will be :code:`8 8 13` followed
+by 64 (8x8) lines each containing 13 coefficients separated by space.
+
+----------
+
+.. include:: accel_styles.rst
+
+----------
+
+Mixing, shift, table, tail correction, restart, rRESPA info
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+For atom type pairs I,J and I != J, where types I and J correspond to
+two different element types, mixing is performed by LAMMPS as described
+above from values in the potential file.
+
+This pair style does not support the :doc:`pair_modify <pair_modify>`
+shift, table, and tail options.
+
+This pair style does not write its information to :doc:`binary restart
+files <restart>`, since it is stored in potential file.
+
+This pair style can only be used via the *pair* keyword of the
+:doc:`run_style respa <run_style>` command.  It does not support the
+*inner*, *middle*, *outer* keywords.
+
+Restrictions
+""""""""""""
+
+The 'uf3' pair style is part of the ML-UF3 package. It is only enabled
+if LAMMPS was built with that package. See the :doc:`Build package
+<Build_package>` page for more info.
+
+This pair style requires the :doc:`newton <newton>` setting to be "on".
+
+The UF3 LAMMPS potential file provided with LAMMPS (see the potentials
+directory) are parameterized for metal :doc:`units <units>`.
+
+The single() function of 'uf3' pair style only return the 2-body
+interaction energy.
+
+Related commands
+""""""""""""""""
+
+:doc:`pair_coeff <pair_coeff>`
+
+Default
+"""""""
+
+none
+
+----------
+
+.. _Xie23:
+
+**(Xie23)** Xie, S.R., Rupp, M. & Hennig, R.G. Ultra-fast interpretable machine-learning potentials. npj Comput Mater 9, 162 (2023). https://doi.org/10.1038/s41524-023-01092-7

doc/src/replicate.rst

@@ -8,56 +8,118 @@ Syntax
 
 .. code-block:: LAMMPS
 
-   replicate nx ny nz *keyword*
+   replicate nx ny nz keyword ...
 
 nx,ny,nz = replication factors in each dimension
 
-* optional *keyword* = *bbox*
+* zero or more keywords may be appended
+* keyword = *bbox* or *bond/periodic*
 
   .. parsed-literal::
 
-       *bbox* = only check atoms in replicas that overlap with a processor's subdomain
+       *bbox* = use a bounding-box algorithm which is faster for large proc counts
+       *bond/periodic* = use an algorithm that correctly replicates periodic bond loops
 
 Examples
 """"""""
 
+For examples of replicating simple linear polymer chains (periodic or
+non-periodic) or periodic carbon nanotubes, see examples/replicate.
+
 .. code-block:: LAMMPS
 
    replicate 2 3 2
+   replicate 2 3 2 bbox
+   replicate 2 3 2 bond/periodic
 
 Description
 """""""""""
 
-Replicate the current simulation one or more times in each dimension.
-For example, replication factors of 2,2,2 will create a simulation
-with 8x as many atoms by doubling the simulation domain in each
-dimension.  A replication factor of 1 in a dimension leaves the
-simulation domain unchanged.  When the new simulation box is created
-it is also partitioned into a regular 3d grid of rectangular bricks,
-one per processor, based on the number of processors being used and
-the settings of the :doc:`processors <processors>` command.  The
-partitioning can later be changed by the :doc:`balance <balance>` or
-:doc:`fix balance <fix_balance>` commands.
+Replicate the current system one or more times in each dimension.  For
+example, replication factors of 2,2,2 will create a simulation with 8x
+as many atoms by doubling the size of the simulation box in each
+dimension.  A replication factor of 1 leaves the simulation domain
+unchanged in that dimension.
 
-All properties of the atoms are replicated, including their
-velocities, which may or may not be desirable.  New atom IDs are
-assigned to new atoms, as are molecule IDs.  Bonds and other topology
-interactions are created between pairs of new atoms as well as between
-old and new atoms.  This is done by using the image flag for each atom
-to "unwrap" it out of the periodic box before replicating it.
+When the new simulation box is created it is partitioned into a
+regular 3d grid of rectangular bricks, one per processor, based on the
+number of processors being used and the settings of the
+:doc:`processors <processors>` command.  The partitioning can be
+changed by subsequent :doc:`balance <balance>` or :doc:`fix balance
+<fix_balance>` commands.
 
-This means that any molecular bond you specify in the original data
-file that crosses a periodic boundary should be between two atoms with
-image flags that differ by 1.  This will allow the bond to be
-unwrapped appropriately.
+All properties of each atom are replicated (except per-atom fix data,
+see the Restrictions section below).  This includes their velocities,
+which may or may not be desirable.  New atom IDs are assigned to new
+atoms, as are new molecule IDs.  Bonds and other topology interactions
+are created between pairs of new atoms as well as between old and new
+atoms.
 
-The optional keyword *bbox* uses a bounding box to only check atoms in
-replicas that overlap with a processor's subdomain when assigning
-atoms to processors.  It typically results in a substantial speedup
-when using the replicate command on a large number of processors.  It
-does require temporary use of more memory, specifically that each
-processor can store all atoms in the entire system before it is
-replicated.
+.. note::
+
+   The bond discussion which follows only refers to models with
+   permanent covalent bonds typically defined in LAMMPS via a data
+   file.  It is not relevant to systems modeled with many-body
+   potentials which can define bonds on-the-fly, based on the current
+   positions of nearby atoms, e.g. models using the :doc:`AIREBO
+   <pair_airebo>` or :doc:`ReaxFF <pair_reaxff>` potentials.
+
+If the *bond/periodic* keyword is not specified, bond replication is
+done by using the image flag for each atom to "unwrap" it out of the
+periodic box before replicating it.  After replication is performed,
+atoms outside the new periodic box are wrapped back into it.  This
+assigns correct images flags to all atoms in the system.  For this to
+work, all original atoms in the original simulation box must have
+consistent image flags.  This means that if two atoms have a bond
+between them which crosses a periodic boundary, their respective image
+flags will differ by 1 in that dimension.
+
+Image flag consistency is not possible if a system has a periodic bond
+loop, meaning there is a chain of bonds which crosses an entire
+dimension and re-connects to itself across a periodic boundary.  In
+this case you MUST use the *bond/periodic* keyword to correctly
+replicate the system.  This option zeroes the image flags for all
+atoms and uses a different algorithm to find new (nearby) bond
+neighbors in the replicated system.  In the final replicated system
+all image flags are zero (in each dimension).
+
+.. note::
+
+   LAMMPS does not check for image flag consistency before performing
+   the replication (it does issue a warning about this before a
+   simulation is run).  If the original image flags are inconsistent,
+   the replicated system will also have inconsistent image flags, but
+   will otherwise be correctly replicated.  This is NOT the case if
+   there is a periodic bond loop.  See the next note.
+
+.. note::
+
+   LAMMPS does not check for periodic bond loops.  If you use the
+   *bond/periodic* keyword for a system without periodic bond loops,
+   the system will be correctly replicated, but image flag information
+   will be lost (which may or may not be important to your model).  If
+   you do not use the *bond/periodic* keyword for a system with
+   periodic bond loops, the replicated system will have invalid bonds
+   (typically very long), resulting in bad dynamics.
+
+If possible, the *bbox* keyword should be used when running on a large
+number of processors, as it can result in a substantial speed-up for
+the replication operation.  It uses a bounding box to only check atoms
+in replicas that overlap with each processor's new subdomain when
+assigning atoms to processors.  It also preserves image flag
+information.  The only drawback to the *bbox* option is that it
+requires a temporary use of more memory.  Each processor must be able
+to store all atoms (and their per-atom data) in the original system,
+before it is replicated.
+
+.. note::
+
+  The algorithm used by the *bond/periodic* keyword builds on the
+  algorithm used by the *bbox* keyword and thus has the same memory
+  requirements.  If you specify only the *bond/peridoic* keyword it
+  will internally set the *bbox* keyword as well.
+
+----------
 
 Restrictions
 """"""""""""
@@ -65,49 +127,30 @@ Restrictions
 A 2d simulation cannot be replicated in the z dimension.
 
 If a simulation is non-periodic in a dimension, care should be used
-when replicating it in that dimension, as it may put atoms nearly on
-top of each other.
-
-.. note::
-
-   You cannot use the replicate command on a system which has a
-   molecule that spans the box and is bonded to itself across a periodic
-   boundary, so that the molecule is effectively a loop.  A simple
-   example would be a linear polymer chain that spans the simulation box
-   and bonds back to itself across the periodic boundary.  More realistic
-   examples would be a CNT (meant to be an infinitely long CNT) or a
-   graphene sheet or a bulk periodic crystal where there are explicit
-   bonds specified between near neighbors.  (Note that this only applies
-   to systems that have permanent bonds as specified in the data file.  A
-   CNT that is just atoms modeled with the :doc:`AIREBO potential <pair_airebo>` has no such permanent bonds, so it can be
-   replicated.)  The reason replication does not work with those systems
-   is that the image flag settings described above cannot be made
-   consistent.  I.e. it is not possible to define images flags so that
-   when every pair of bonded atoms is unwrapped (using the image flags),
-   they will be close to each other.  The only way the replicate command
-   could work in this scenario is for it to break a bond, insert more
-   atoms, and re-connect the loop for the larger simulation box.  But it
-   is not clever enough to do this.  So you will have to construct a
-   larger version of your molecule as a pre-processing step and input a
-   new data file to LAMMPS.
+when replicating it in that dimension, as it may generate atoms nearly
+on top of each other.
 
 If the current simulation was read in from a restart file (before a
-run is performed), there must not be any fix information stored in
-the file for individual atoms.  Similarly, no fixes can be defined at
-the time the replicate command is used that require vectors of atom
+run is performed), there must not be any fix information stored in the
+file for individual atoms.  Similarly, no fixes can be defined at the
+time the replicate command is used that require vectors of atom
 information to be stored.  This is because the replicate command does
 not know how to replicate that information for new atoms it creates.
-To work around this restriction, restart files may be converted into
-data files and fixes may be undefined via the :doc:`unfix <unfix>`
-command before and redefined after the replicate command.
+
+To work around this restriction two options are possible.  (1) Fixes
+which use the stored data in the restart file can be defined before
+replication and then deleted via the :doc:`unfix <unfix>` command and
+re-defined after it.  Or (2) the restart file can be converted to a
+data file (which deletes the stored fix information) and fixes defined
+after the replicate command.  In both these scenarios, the per-atom
+fix information in the restart file is lost.
 
 Related commands
 """"""""""""""""
 
 none
 
-
 Default
 """""""
 
-none
+No settings for using the *bbox* or *bond/periodic* algorithms.

doc/src/run_style.rst

@@ -327,10 +327,12 @@ Restrictions
 """"""""""""
 
 The *verlet/split* style can only be used if LAMMPS was built with the
-REPLICA package. Correspondingly the *respa/omp* style is available
+REPLICA package.  Correspondingly the *respa/omp* style is available
 only if the OPENMP package was included. See the :doc:`Build package
-<Build_package>` page for more info.  It is not compatible with
-kspace styles from the INTEL package.
+<Build_package>` page for more info.
+
+Run style *verlet/split* is not compatible with kspace styles from
+the INTEL package and it is not compatible with any TIP4P styles.
 
 Whenever using rRESPA, the user should experiment with trade-offs in
 speed and accuracy for their system, and verify that they are

doc/src/variable.rst

@@ -67,7 +67,7 @@ Syntax
                            bound(group,dir,region), gyration(group,region), ke(group,reigon),
                            angmom(group,dim,region), torque(group,dim,region),
                            inertia(group,dimdim,region), omega(group,dim,region)
-         special functions = sum(x), min(x), max(x), ave(x), trap(x), slope(x), gmask(x), rmask(x), grmask(x,y), next(x), is_file(name), is_os(name), extract_setting(name), label2type(kind,label), is_typelabel(kind,label)
+         special functions = sum(x), min(x), max(x), ave(x), trap(x), slope(x), sort(x), rsort(x), gmask(x), rmask(x), grmask(x,y), next(x), is_file(name), is_os(name), extract_setting(name), label2type(kind,label), is_typelabel(kind,label)
          feature functions = is_available(category,feature), is_active(category,feature), is_defined(category,id)
          atom value = id[i], mass[i], type[i], mol[i], x[i], y[i], z[i], vx[i], vy[i], vz[i], fx[i], fy[i], fz[i], q[i]
          atom vector = id, mass, type, mol, radius, q, x, y, z, vx, vy, vz, fx, fy, fz
@@ -547,7 +547,7 @@ variables.
 +------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
 | Region functions       | count(ID,IDR), mass(ID,IDR), charge(ID,IDR), xcm(ID,dim,IDR), vcm(ID,dim,IDR), fcm(ID,dim,IDR), bound(ID,dir,IDR), gyration(ID,IDR), ke(ID,IDR), angmom(ID,dim,IDR), torque(ID,dim,IDR), inertia(ID,dimdim,IDR), omega(ID,dim,IDR)                                                                                                                |
 +------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
-| Special functions      | sum(x), min(x), max(x), ave(x), trap(x), slope(x), gmask(x), rmask(x), grmask(x,y), next(x), is_file(name), is_os(name), extract_setting(name), label2type(kind,label), is_typelabel(kind,label)                                                                                                                                                  |
+| Special functions      | sum(x), min(x), max(x), ave(x), trap(x), slope(x), sort(x), rsort(x), gmask(x), rmask(x), grmask(x,y), next(x), is_file(name), is_os(name), extract_setting(name), label2type(kind,label), is_typelabel(kind,label)                                                                                                                               |
 +------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
 | Feature functions      | is_available(category,feature), is_active(category,feature), is_defined(category,id)                                                                                                                                                                                                                                                              |
 +------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
@@ -913,23 +913,27 @@ Special Functions
 Special functions take specific kinds of arguments, meaning their
 arguments cannot be formulas themselves.
 
-The sum(x), min(x), max(x), ave(x), trap(x), and slope(x) functions
-each take 1 argument which is of the form "c_ID" or "c_ID[N]" or
-"f_ID" or "f_ID[N]" or "v_name".  The first two are computes and the
-second two are fixes; the ID in the reference should be replaced by
-the ID of a compute or fix defined elsewhere in the input script.  The
-compute or fix must produce either a global vector or array.  If it
-produces a global vector, then the notation without "[N]" should be
-used.  If it produces a global array, then the notation with "[N]"
-should be used, when N is an integer, to specify which column of the
-global array is being referenced.  The last form of argument "v_name"
-is for a vector-style variable where "name" is replaced by the name of
-the variable.
+The sum(x), min(x), max(x), ave(x), trap(x), slope(x), sort(x), and
+rsort(x) functions each take 1 argument which is of the form "c_ID" or
+"c_ID[N]" or "f_ID" or "f_ID[N]" or "v_name".  The first two are
+computes and the second two are fixes; the ID in the reference should be
+replaced by the ID of a compute or fix defined elsewhere in the input
+script.  The compute or fix must produce either a global vector or
+array.  If it produces a global vector, then the notation without "[N]"
+should be used.  If it produces a global array, then the notation with
+"[N]" should be used, where N is an integer, to specify which column of
+the global array is being referenced.  The last form of argument
+"v_name" is for a vector-style variable where "name" is replaced by the
+name of the variable.
 
-These functions operate on a global vector of inputs and reduce it to
-a single scalar value.  This is analogous to the operation of the
-:doc:`compute reduce <compute_reduce>` command, which performs similar
-operations on per-atom and local vectors.
+The sum(x), min(x), max(x), ave(x), trap(x), and slope(x) functions
+operate on a global vector of inputs and reduce it to a single scalar
+value.  This is analogous to the operation of the :doc:`compute reduce
+<compute_reduce>` command, which performs similar operations on per-atom
+and local vectors.
+
+The sort(x) and rsort(x) functions operate on a global vector of inputs
+and return a global vector of the same length.
 
 The sum() function calculates the sum of all the vector elements.  The
 min() and max() functions find the minimum and maximum element
@@ -953,6 +957,12 @@ of points, equally spaced by 1 in their x coordinate: (1,V1), (2,V2),
 length N.  The returned value is the slope of the line.  If the line
 has a single point or is vertical, it returns 1.0e20.
 
+.. versionadded:: 27June2024
+
+The sort(x) and rsort(x) functions sort the data of the input vector by
+their numeric value: sort(x) sorts in ascending order, rsort(x) sorts
+in descending order.
+
 The gmask(x) function takes 1 argument which is a group ID.  It
 can only be used in atom-style variables.  It returns a 1 for
 atoms that are in the group, and a 0 for atoms that are not.

doc/src/write_data.rst

@@ -12,14 +12,14 @@ Syntax
 
 * file = name of data file to write out
 * zero or more keyword/value pairs may be appended
-* keyword = *pair* or *nocoeff* or *nofix* or *nolabelmap*
+* keyword = *nocoeff* or *nofix* or *nolabelmap* or *triclinic/general* or *types* or *pair*
 
   .. parsed-literal::
 
        *nocoeff* = do not write out force field info
        *nofix* = do not write out extra sections read by fixes
        *nolabelmap* = do not write out type labels
-       *triclinic/general = write data file in general triclinic format
+       *triclinic/general* = write data file in general triclinic format
        *types* value = *numeric* or *labels*
        *pair* value = *ii* or *ij*
          *ii* = write one line of pair coefficient info per atom type
@@ -189,4 +189,4 @@ Related commands
 Default
 """""""
 
-The option defaults are pair = ii and types_style = numeric.
+The option defaults are pair = ii and types = numeric.

doc/utils/requirements.txt

@@ -1,6 +1,7 @@
-Sphinx >= 5.3.0, <7.3
+Sphinx >= 5.3.0, <8.0
 sphinxcontrib-spelling
 sphinxcontrib-jquery
+sphinx-design
 git+https://github.com/akohlmey/sphinx-fortran@parallel-read
 sphinx-tabs>=3.4.1
 breathe

doc/utils/sphinx-config/conf.py.in

@@ -57,6 +57,7 @@ extensions = [
     'table_from_list',
     'tab_or_note',
     'breathe',
+    'sphinx_design'
 ]
 
 images_config = {

doc/utils/sphinx-config/false_positives.txt

@@ -992,6 +992,7 @@ emax
 Emax
 Embt
 emi
+Emilie
 Emmrich
 emol
 eN
@@ -1172,6 +1173,7 @@ finitecutflag
 Finnis
 Fiorin
 fitpod
+fivebody
 fixID
 fj
 Fji
@@ -1433,6 +1435,7 @@ Hendrik
 Henin
 Henkelman
 Henkes
+Hennig
 henrich
 Henrich
 Hermitian
@@ -1594,6 +1597,7 @@ interlayer
 intermolecular
 interoperable
 Interparticle
+interpretable
 interstitials
 intertube
 Intr
@@ -1732,6 +1736,7 @@ Kalia
 Kamberaj
 Kantorovich
 Kapfer
+Kapil
 Karhunen
 Karls
 Karlsruhe
@@ -1763,8 +1768,10 @@ keflag
 Keir
 Kelchner
 Kelkar
+Kemppainen
 Kemper
 kepler
+Kemppainen
 keV
 Keyes
 keyfile
@@ -1813,6 +1820,7 @@ Koziol
 Kp
 kradius
 Kraker
+Krass
 Kraus
 Kremer
 Kress
@@ -2483,6 +2491,7 @@ Nevery
 newfile
 Newns
 newtype
+nextsort
 Neyts
 Nf
 nfft
@@ -2672,6 +2681,7 @@ nzlo
 ocl
 octahedral
 octants
+Odegard
 Ohara
 O'Hearn
 ohenrich
@@ -2888,6 +2898,7 @@ Pmolrotate
 Pmoltrans
 pN
 png
+podd
 Podhorszki
 Poiseuille
 poisson
@@ -2953,6 +2964,7 @@ Priya
 proc
 Proc
 procs
+procgrid
 progguide
 Prony
 ps
@@ -3252,6 +3264,7 @@ rRESPA
 Rsi
 Rso
 Rspace
+rsort
 rsq
 rst
 rstyle
@@ -3264,6 +3277,7 @@ Rudranarayan
 Rudzinski
 Runge
 runtime
+Rupp
 Rutuparna
 rx
 rxd
@@ -3287,6 +3301,7 @@ Saidi
 saip
 Salanne
 Salles
+sametag
 sandia
 Sandia
 sandybrown
@@ -3357,6 +3372,7 @@ setmask
 Setmask
 setpoint
 setvel
+sevenbody
 sfftw
 sfree
 Sg
@@ -3404,8 +3420,10 @@ sinh
 sinusoid
 sinusoidally
 SiO
+Siochi
 Sirk
 Sival
+sixbody
 sizeI
 sizeJ
 sizex
@@ -3451,6 +3469,7 @@ solvated
 solvation
 someuser
 Sorensen
+sortfreq
 soundspeed
 sourceforge
 Souza
@@ -3720,7 +3739,6 @@ tokyo
 tol
 tomic
 toolchain
-toolset
 topologies
 Toporov
 Torder
@@ -3798,6 +3816,7 @@ typeJ
 typelabel
 typeN
 typesafe
+typestr
 Tz
 Tzou
 ub
@@ -3810,6 +3829,8 @@ uChem
 uCond
 uef
 UEF
+uf
+uf3
 ufm
 Uhlenbeck
 Ui
@@ -4149,6 +4170,7 @@ yy
 yz
 Zagaceta
 Zannoni
+Zavada
 Zavattieri
 zbl
 ZBL
@@ -4162,6 +4184,7 @@ zenodo
 Zepeda
 zflag
 Zhang
+Zhao
 Zhen
 zhi
 Zhigilei

examples/COUPLE/plugin/liblammpsplugin.c

@@ -101,6 +101,7 @@ liblammpsplugin_t *liblammpsplugin_load(const char *lib)
   ADDSYM(extract_setting);
   ADDSYM(extract_global_datatype);
   ADDSYM(extract_global);
+  ADDSYM(map_atom);
 
   ADDSYM(extract_atom_datatype);
   ADDSYM(extract_atom);

examples/COUPLE/plugin/liblammpsplugin.h

@@ -146,6 +146,7 @@ struct _liblammpsplugin {
   int (*extract_setting)(void *, const char *);
   int *(*extract_global_datatype)(void *, const char *);
   void *(*extract_global)(void *, const char *);
+  void *(*map_atom)(void *, const void *);
 
   int *(*extract_atom_datatype)(void *, const char *);
   void *(*extract_atom)(void *, const char *);