Merge branch 'develop' into remove-fix-latte

.github/CODEOWNERS

@@ -150,12 +150,12 @@ tools/vim/*           @hammondkd
 unittest/*            @akohlmey
 
 # cmake
-cmake/*               @junghans @rbberger
+cmake/*               @rbberger
 cmake/Modules/LAMMPSInterfacePlugin.cmake @akohlmey
 cmake/Modules/MPI4WIN.cmake @akohlmey
 cmake/Modules/OpenCLLoader.cmake @akohlmey
-cmake/Modules/Packages/COLVARS.cmake @junghans @rbberger @giacomofiorin
+cmake/Modules/Packages/COLVARS.cmake @rbberger @giacomofiorin
-cmake/Modules/Packages/KIM.cmake @junghans @rbberger @ellio167
+cmake/Modules/Packages/KIM.cmake @rbberger @ellio167
 cmake/presets/*.cmake @akohlmey
 
 # python

cmake/CMakeLists.txt

@@ -538,7 +538,10 @@ set(CMAKE_TUNE_FLAGS "${CMAKE_TUNE_DEFAULT}" CACHE STRING "Compiler and machine
 separate_arguments(CMAKE_TUNE_FLAGS)
 foreach(_FLAG ${CMAKE_TUNE_FLAGS})
   target_compile_options(lammps PRIVATE ${_FLAG})
-  target_compile_options(lmp PRIVATE ${_FLAG})
+  # skip these flags when linking the main executable
+  if(NOT (("${_FLAG}" STREQUAL "-Xcudafe") OR (("${_FLAG}" STREQUAL "--diag_suppress=unrecognized_pragma"))))
+    target_compile_options(lmp PRIVATE ${_FLAG})
+  endif()
 endforeach()
 ########################################################################
 # Basic system tests (standard libraries, headers, functions, types)   #
@@ -838,9 +841,8 @@ if(BUILD_SHARED_LIBS)
   set(LIBLAMMPS_SHARED_BINARY ${MY_BUILD_DIR}/liblammps${LAMMPS_MACHINE}${CMAKE_SHARED_LIBRARY_SUFFIX})
   if(Python_EXECUTABLE)
     add_custom_target(
-      install-python ${CMAKE_COMMAND} -E remove_directory build
+      install-python ${Python_EXECUTABLE} ${LAMMPS_PYTHON_DIR}/install.py -p ${LAMMPS_PYTHON_DIR}/lammps
-      COMMAND ${Python_EXECUTABLE} ${LAMMPS_PYTHON_DIR}/install.py -p ${LAMMPS_PYTHON_DIR}/lammps
+      -l ${LIBLAMMPS_SHARED_BINARY} -w ${MY_BUILD_DIR} -v ${LAMMPS_SOURCE_DIR}/version.h
-              -l ${LIBLAMMPS_SHARED_BINARY} -w ${MY_BUILD_DIR}
       COMMENT "Installing LAMMPS Python module")
   else()
     add_custom_target(
@@ -853,35 +855,6 @@ else()
     ${CMAKE_COMMAND} -E echo "Must build LAMMPS as a shared library to use the Python module")
 endif()
 
-###############################################################################
-# Add LAMMPS python module to "install" target. This is taylored for building
-# LAMMPS for package managers and with different prefix settings.
-# This requires either a shared library or that the PYTHON package is included.
-###############################################################################
-if(BUILD_SHARED_LIBS OR PKG_PYTHON)
-  if(CMAKE_VERSION VERSION_LESS 3.12)
-    # adjust so we find Python 3 versions before Python 2 on old systems with old CMake
-    set(Python_ADDITIONAL_VERSIONS 3.12 3.11 3.10 3.9 3.8 3.7 3.6)
-    find_package(PythonInterp) # Deprecated since version 3.12
-    if(PYTHONINTERP_FOUND)
-      set(Python_EXECUTABLE ${PYTHON_EXECUTABLE})
-    endif()
-  else()
-    # backward compatibility
-    if(PYTHON_EXECUTABLE)
-      set(Python_EXECUTABLE ${PYTHON_EXECUTABLE})
-    endif()
-    find_package(Python COMPONENTS Interpreter)
-  endif()
-  if(Python_EXECUTABLE)
-    file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/python/lib)
-    file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/python/src)
-    file(COPY ${LAMMPS_SOURCE_DIR}/version.h  DESTINATION ${CMAKE_BINARY_DIR}/python/src)
-    file(COPY ${LAMMPS_PYTHON_DIR}/README ${LAMMPS_PYTHON_DIR}/pyproject.toml ${LAMMPS_PYTHON_DIR}/setup.py ${LAMMPS_PYTHON_DIR}/lammps  DESTINATION ${CMAKE_BINARY_DIR}/python/lib)
-    install(CODE "if(\"\$ENV{DESTDIR}\" STREQUAL \"\")\n execute_process(COMMAND ${Python_EXECUTABLE} -m pip install -v ${CMAKE_BINARY_DIR}/python/lib --prefix=${CMAKE_INSTALL_PREFIX})\n else()\n execute_process(COMMAND ${Python_EXECUTABLE} -m pip install -v ${CMAKE_BINARY_DIR}/python/lib --prefix=${CMAKE_INSTALL_PREFIX} --root=\$ENV{DESTDIR})\n endif()")
-  endif()
-endif()
-
 include(Testing)
 include(CodeCoverage)
 include(CodingStandard)

cmake/Modules/Packages/GPU.cmake

@@ -60,9 +60,9 @@ if(GPU_API STREQUAL "CUDA")
   option(CUDA_MPS_SUPPORT "Enable tweaks to support CUDA Multi-process service (MPS)" OFF)
   if(CUDA_MPS_SUPPORT)
     if(CUDPP_OPT)
-      message(FATAL_ERROR "Must use -DCUDPP_OPT=OFF with -DGPU_CUDA_MPS_SUPPORT=ON")
+      message(FATAL_ERROR "Must use -DCUDPP_OPT=OFF with -DCUDA_MPS_SUPPORT=ON")
     endif()
-    set(GPU_CUDA_MPS_FLAGS "-DCUDA_PROXY")
+    set(GPU_CUDA_MPS_FLAGS "-DCUDA_MPS_SUPPORT")
   endif()
 
   set(GPU_ARCH "sm_50" CACHE STRING "LAMMPS GPU CUDA SM primary architecture (e.g. sm_60)")
@@ -98,9 +98,11 @@ if(GPU_API STREQUAL "CUDA")
   # comparison chart according to: https://en.wikipedia.org/wiki/CUDA#GPUs_supported
   if(CUDA_VERSION VERSION_LESS 8.0)
     message(FATAL_ERROR "CUDA Toolkit version 8.0 or later is required")
-  elseif(CUDA_VERSION VERSION_GREATER_EQUAL "12.0")
+  elseif(CUDA_VERSION VERSION_GREATER_EQUAL "13.0")
     message(WARNING "Untested CUDA Toolkit version ${CUDA_VERSION}. Use at your own risk")
     set(GPU_CUDA_GENCODE "-arch=all")
+  elseif(CUDA_VERSION VERSION_GREATER_EQUAL "12.0")
+    set(GPU_CUDA_GENCODE "-arch=all")
   else()
     # Kepler (GPU Arch 3.0) is supported by CUDA 5 to CUDA 10.2
     if((CUDA_VERSION VERSION_GREATER_EQUAL "5.0") AND (CUDA_VERSION VERSION_LESS "11.0"))
@@ -175,8 +177,6 @@ if(GPU_API STREQUAL "CUDA")
     target_compile_definitions(gpu PRIVATE -DUSE_CUDPP)
   endif()
 
-  target_link_libraries(lammps PRIVATE gpu)
-
   add_executable(nvc_get_devices ${LAMMPS_LIB_SOURCE_DIR}/gpu/geryon/ucl_get_devices.cpp)
   target_compile_definitions(nvc_get_devices PRIVATE -DUCL_CUDADR)
   target_link_libraries(nvc_get_devices PRIVATE ${CUDA_LIBRARIES} ${CUDA_CUDA_LIBRARY})
@@ -249,12 +249,12 @@ elseif(GPU_API STREQUAL "OPENCL")
   else()
     target_compile_definitions(gpu PRIVATE -DMPI_GERYON -DGERYON_NUMA_FISSION -DUCL_NO_EXIT)
   endif()
-  target_link_libraries(lammps PRIVATE gpu)
 
   add_executable(ocl_get_devices ${LAMMPS_LIB_SOURCE_DIR}/gpu/geryon/ucl_get_devices.cpp)
   target_compile_definitions(ocl_get_devices PRIVATE -DUCL_OPENCL)
   target_link_libraries(ocl_get_devices PRIVATE OpenCL::OpenCL)
   add_dependencies(ocl_get_devices OpenCL::OpenCL)
+
 elseif(GPU_API STREQUAL "HIP")
   if(NOT DEFINED HIP_PATH)
       if(NOT DEFINED ENV{HIP_PATH})
@@ -285,7 +285,7 @@ elseif(GPU_API STREQUAL "HIP")
 
   set(ENV{HIP_PLATFORM} ${HIP_PLATFORM})
 
-  if(HIP_PLATFORM STREQUAL "hcc" OR HIP_PLATFORM STREQUAL "amd")
+  if(HIP_PLATFORM STREQUAL "amd")
     set(HIP_ARCH "gfx906" CACHE STRING "HIP target architecture")
   elseif(HIP_PLATFORM STREQUAL "spirv")
     set(HIP_ARCH "spirv" CACHE STRING "HIP target architecture")
@@ -358,7 +358,7 @@ elseif(GPU_API STREQUAL "HIP")
     set(CUBIN_FILE   "${LAMMPS_LIB_BINARY_DIR}/gpu/${CU_NAME}.cubin")
     set(CUBIN_H_FILE "${LAMMPS_LIB_BINARY_DIR}/gpu/${CU_NAME}_cubin.h")
 
-    if(HIP_PLATFORM STREQUAL "hcc" OR HIP_PLATFORM STREQUAL "amd")
+    if(HIP_PLATFORM STREQUAL "amd")
         configure_file(${CU_FILE} ${CU_CPP_FILE} COPYONLY)
 
         if(HIP_COMPILER STREQUAL "clang")
@@ -412,7 +412,8 @@ elseif(GPU_API STREQUAL "HIP")
       set_property(TARGET gpu PROPERTY CXX_STANDARD 14)
     endif()
     # add hipCUB
-    target_include_directories(gpu PRIVATE ${HIP_ROOT_DIR}/../include)
+    find_package(hipcub REQUIRED)
+    target_link_libraries(gpu PRIVATE hip::hipcub)
     target_compile_definitions(gpu PRIVATE -DUSE_HIP_DEVICE_SORT)
 
     if(HIP_PLATFORM STREQUAL "nvcc")
@@ -461,35 +462,25 @@ elseif(GPU_API STREQUAL "HIP")
 
   add_executable(hip_get_devices ${LAMMPS_LIB_SOURCE_DIR}/gpu/geryon/ucl_get_devices.cpp)
   target_compile_definitions(hip_get_devices PRIVATE -DUCL_HIP)
-  target_link_libraries(hip_get_devices hip::host)
+  target_link_libraries(hip_get_devices PRIVATE hip::host)
 
   if(HIP_PLATFORM STREQUAL "nvcc")
     target_compile_definitions(gpu PRIVATE -D__HIP_PLATFORM_NVCC__)
-    target_include_directories(gpu PRIVATE ${HIP_ROOT_DIR}/../include)
     target_include_directories(gpu PRIVATE ${CUDA_INCLUDE_DIRS})
     target_link_libraries(gpu PRIVATE ${CUDA_LIBRARIES} ${CUDA_CUDA_LIBRARY})
 
     target_compile_definitions(hip_get_devices PRIVATE -D__HIP_PLATFORM_NVCC__)
-    target_include_directories(hip_get_devices PRIVATE ${HIP_ROOT_DIR}/include)
     target_include_directories(hip_get_devices PRIVATE ${CUDA_INCLUDE_DIRS})
     target_link_libraries(hip_get_devices PRIVATE ${CUDA_LIBRARIES} ${CUDA_CUDA_LIBRARY})
-  elseif(HIP_PLATFORM STREQUAL "hcc")
-    target_compile_definitions(gpu PRIVATE -D__HIP_PLATFORM_HCC__)
-    target_include_directories(gpu PRIVATE ${HIP_ROOT_DIR}/../include)
-
-    target_compile_definitions(hip_get_devices PRIVATE -D__HIP_PLATFORM_HCC__)
-    target_include_directories(hip_get_devices PRIVATE ${HIP_ROOT_DIR}/../include)
-  elseif(HIP_PLATFORM STREQUAL "amd")
-    target_compile_definitions(gpu PRIVATE -D__HIP_PLATFORM_AMD__)
-    target_include_directories(gpu PRIVATE ${HIP_ROOT_DIR}/../include)
-
-    target_compile_definitions(hip_get_devices PRIVATE -D__HIP_PLATFORM_AMD__)
-    target_include_directories(hip_get_devices PRIVATE ${HIP_ROOT_DIR}/../include)
   endif()
-
-  target_link_libraries(lammps PRIVATE gpu)
 endif()
 
+if(BUILD_OMP)
+  find_package(OpenMP COMPONENTS CXX REQUIRED)
+  target_link_libraries(gpu PRIVATE OpenMP::OpenMP_CXX)
+endif()
+target_link_libraries(lammps PRIVATE gpu)
+
 set_property(GLOBAL PROPERTY "GPU_SOURCES" "${GPU_SOURCES}")
 # detect styles which have a GPU version
 RegisterStylesExt(${GPU_SOURCES_DIR} gpu GPU_SOURCES)

cmake/Modules/Packages/KIM.cmake

@@ -19,7 +19,7 @@ if(CURL_FOUND)
     target_compile_definitions(lammps PRIVATE -DLMP_NO_SSL_CHECK)
   endif()
 endif()
-set(KIM_EXTRA_UNITTESTS OFF CACHE STRING "Set extra unit tests verbose mode on/off. If on, extra tests are included.")
+option(KIM_EXTRA_UNITTESTS "Enable extra unit tests for the KIM package." OFF)
 mark_as_advanced(KIM_EXTRA_UNITTESTS)
 find_package(PkgConfig QUIET)
 set(DOWNLOAD_KIM_DEFAULT ON)

cmake/Modules/Packages/KOKKOS.cmake

@@ -72,13 +72,11 @@ if(DOWNLOAD_KOKKOS)
   set_target_properties(LAMMPS::KOKKOSCONTAINERS PROPERTIES
     IMPORTED_LOCATION "${INSTALL_DIR}/lib/libkokkoscontainers.a")
   target_link_libraries(lammps PRIVATE LAMMPS::KOKKOSCORE LAMMPS::KOKKOSCONTAINERS)
-  target_link_libraries(lmp PRIVATE LAMMPS::KOKKOSCORE LAMMPS::KOKKOSCONTAINERS)
   add_dependencies(LAMMPS::KOKKOSCORE kokkos_build)
   add_dependencies(LAMMPS::KOKKOSCONTAINERS kokkos_build)
 elseif(EXTERNAL_KOKKOS)
   find_package(Kokkos 3.7.01 REQUIRED CONFIG)
   target_link_libraries(lammps PRIVATE Kokkos::kokkos)
-  target_link_libraries(lmp PRIVATE Kokkos::kokkos)
 else()
   set(LAMMPS_LIB_KOKKOS_SRC_DIR ${LAMMPS_LIB_SOURCE_DIR}/kokkos)
   set(LAMMPS_LIB_KOKKOS_BIN_DIR ${LAMMPS_LIB_BINARY_DIR}/kokkos)
@@ -98,7 +96,6 @@ else()
                           ${LAMMPS_LIB_KOKKOS_BIN_DIR})
   target_include_directories(lammps PRIVATE ${Kokkos_INCLUDE_DIRS})
   target_link_libraries(lammps PRIVATE kokkos)
-  target_link_libraries(lmp PRIVATE kokkos)
   if(BUILD_SHARED_LIBS_WAS_ON)
     set(BUILD_SHARED_LIBS ON)
   endif()

cmake/Modules/Packages/MDI.cmake

@@ -8,8 +8,8 @@ option(DOWNLOAD_MDI "Download and compile the MDI library instead of using an al
 
 if(DOWNLOAD_MDI)
   message(STATUS "MDI download requested - we will build our own")
-  set(MDI_URL "https://github.com/MolSSI-MDI/MDI_Library/archive/v1.4.12.tar.gz" CACHE STRING "URL for MDI tarball")
+  set(MDI_URL "https://github.com/MolSSI-MDI/MDI_Library/archive/v1.4.16.tar.gz" CACHE STRING "URL for MDI tarball")
-  set(MDI_MD5 "7a222353ae8e03961d5365e6cd48baee" CACHE STRING "MD5 checksum for MDI tarball")
+  set(MDI_MD5 "407db44e2d79447ab5c1233af1965f65" CACHE STRING "MD5 checksum for MDI tarball")
   mark_as_advanced(MDI_URL)
   mark_as_advanced(MDI_MD5)
   GetFallbackURL(MDI_URL MDI_FALLBACK)

cmake/Modules/Packages/MSCG.cmake

@@ -7,8 +7,8 @@ else()
 endif()
 option(DOWNLOAD_MSCG "Download MSCG library instead of using an already installed one)" ${DOWNLOAD_MSCG_DEFAULT})
 if(DOWNLOAD_MSCG)
-  set(MSCG_URL "https://github.com/uchicago-voth/MSCG-release/archive/1.7.3.1.tar.gz" CACHE STRING "URL for MSCG tarball")
+  set(MSCG_URL "https://github.com/uchicago-voth/MSCG-release/archive/491270a73539e3f6951e76f7dbe84e258b3ebb45.tar.gz" CACHE STRING "URL for MSCG tarball")
-  set(MSCG_MD5 "8c45e269ee13f60b303edd7823866a91" CACHE STRING "MD5 checksum of MSCG tarball")
+  set(MSCG_MD5 "7ea50748fba5c3a372e0266bd31d2f11" CACHE STRING "MD5 checksum of MSCG tarball")
   mark_as_advanced(MSCG_URL)
   mark_as_advanced(MSCG_MD5)
 

cmake/Modules/Packages/PLUMED.cmake

@@ -54,8 +54,8 @@ if(DOWNLOAD_PLUMED)
     set(PLUMED_BUILD_BYPRODUCTS "<INSTALL_DIR>/lib/${CMAKE_STATIC_LIBRARY_PREFIX}plumedWrapper${CMAKE_STATIC_LIBRARY_PREFIX}")
   endif()
 
-  set(PLUMED_URL "https://github.com/plumed/plumed2/releases/download/v2.8.1/plumed-src-2.8.1.tgz" CACHE STRING "URL for PLUMED tarball")
+  set(PLUMED_URL "https://github.com/plumed/plumed2/releases/download/v2.8.2/plumed-src-2.8.2.tgz" CACHE STRING "URL for PLUMED tarball")
-  set(PLUMED_MD5 "6bfe72ebdae63dc38a9ca27d9b0e08f8" CACHE STRING "MD5 checksum of PLUMED tarball")
+  set(PLUMED_MD5 "599092b6a0aa6fff992612537ad98994" CACHE STRING "MD5 checksum of PLUMED tarball")
 
   mark_as_advanced(PLUMED_URL)
   mark_as_advanced(PLUMED_MD5)

cmake/Modules/StyleHeaderUtils.cmake

@@ -95,73 +95,76 @@ function(RegisterIntegrateStyle path)
 endfunction(RegisterIntegrateStyle)
 
 function(RegisterStyles search_path)
-    FindStyleHeaders(${search_path} ANGLE_CLASS     angle_     ANGLE     ) # angle     ) # force
+    FindStyleHeaders(${search_path} ANGLE_CLASS        angle_        ANGLE        ) # angle        ) # force
-    FindStyleHeaders(${search_path} ATOM_CLASS      atom_vec_  ATOM_VEC  ) # atom      ) # atom      atom_vec_hybrid
+    FindStyleHeaders(${search_path} ATOM_CLASS         atom_vec_     ATOM_VEC     ) # atom         ) # atom      atom_vec_hybrid
-    FindStyleHeaders(${search_path} BODY_CLASS      body_      BODY      ) # body      ) # atom_vec_body
+    FindStyleHeaders(${search_path} BODY_CLASS         body_         BODY         ) # body         ) # atom_vec_body
-    FindStyleHeaders(${search_path} BOND_CLASS      bond_      BOND      ) # bond      ) # force
+    FindStyleHeaders(${search_path} BOND_CLASS         bond_         BOND         ) # bond         ) # force
-    FindStyleHeaders(${search_path} COMMAND_CLASS   "[^.]"     COMMAND   ) # command   ) # input
+    FindStyleHeaders(${search_path} COMMAND_CLASS      "[^.]"        COMMAND      ) # command      ) # input
-    FindStyleHeaders(${search_path} COMPUTE_CLASS   compute_   COMPUTE   ) # compute   ) # modify
+    FindStyleHeaders(${search_path} COMPUTE_CLASS      compute_      COMPUTE      ) # compute      ) # modify
-    FindStyleHeaders(${search_path} DIHEDRAL_CLASS  dihedral_  DIHEDRAL  ) # dihedral  ) # force
+    FindStyleHeaders(${search_path} DIHEDRAL_CLASS     dihedral_     DIHEDRAL     ) # dihedral     ) # force
-    FindStyleHeaders(${search_path} DUMP_CLASS      dump_      DUMP      ) # dump      ) # output    write_dump
+    FindStyleHeaders(${search_path} DUMP_CLASS         dump_         DUMP         ) # dump         ) # output    write_dump
-    FindStyleHeaders(${search_path} FIX_CLASS       fix_       FIX       ) # fix       ) # modify
+    FindStyleHeaders(${search_path} FIX_CLASS          fix_          FIX          ) # fix          ) # modify
-    FindStyleHeaders(${search_path} IMPROPER_CLASS  improper_  IMPROPER  ) # improper  ) # force
+    FindStyleHeaders(${search_path} GRAN_SUB_MOD_CLASS gran_sub_mod_ GRAN_SUB_MOD ) # gran_sub_mod ) # granular_model
-    FindStyleHeaders(${search_path} INTEGRATE_CLASS "[^.]"     INTEGRATE ) # integrate ) # update
+    FindStyleHeaders(${search_path} IMPROPER_CLASS     improper_     IMPROPER     ) # improper     ) # force
-    FindStyleHeaders(${search_path} KSPACE_CLASS    "[^.]"     KSPACE    ) # kspace    ) # force
+    FindStyleHeaders(${search_path} INTEGRATE_CLASS    "[^.]"        INTEGRATE    ) # integrate    ) # update
-    FindStyleHeaders(${search_path} MINIMIZE_CLASS  min_       MINIMIZE  ) # minimize  ) # update
+    FindStyleHeaders(${search_path} KSPACE_CLASS       "[^.]"        KSPACE       ) # kspace       ) # force
-    FindStyleHeaders(${search_path} NBIN_CLASS      nbin_      NBIN      ) # nbin      ) # neighbor
+    FindStyleHeaders(${search_path} MINIMIZE_CLASS     min_          MINIMIZE     ) # minimize     ) # update
-    FindStyleHeaders(${search_path} NPAIR_CLASS     npair_     NPAIR     ) # npair     ) # neighbor
+    FindStyleHeaders(${search_path} NBIN_CLASS         nbin_         NBIN         ) # nbin         ) # neighbor
-    FindStyleHeaders(${search_path} NSTENCIL_CLASS  nstencil_  NSTENCIL  ) # nstencil  ) # neighbor
+    FindStyleHeaders(${search_path} NPAIR_CLASS        npair_        NPAIR        ) # npair        ) # neighbor
-    FindStyleHeaders(${search_path} NTOPO_CLASS     ntopo_     NTOPO     ) # ntopo     ) # neighbor
+    FindStyleHeaders(${search_path} NSTENCIL_CLASS     nstencil_     NSTENCIL     ) # nstencil     ) # neighbor
-    FindStyleHeaders(${search_path} PAIR_CLASS      pair_      PAIR      ) # pair      ) # force
+    FindStyleHeaders(${search_path} NTOPO_CLASS        ntopo_        NTOPO        ) # ntopo        ) # neighbor
-    FindStyleHeaders(${search_path} READER_CLASS    reader_    READER    ) # reader    ) # read_dump
+    FindStyleHeaders(${search_path} PAIR_CLASS         pair_         PAIR         ) # pair         ) # force
-    FindStyleHeaders(${search_path} REGION_CLASS    region_    REGION    ) # region    ) # domain
+    FindStyleHeaders(${search_path} READER_CLASS       reader_       READER       ) # reader       ) # read_dump
+    FindStyleHeaders(${search_path} REGION_CLASS       region_       REGION       ) # region       ) # domain
 endfunction(RegisterStyles)
 
 function(RegisterStylesExt search_path extension sources)
-    FindStyleHeadersExt(${search_path} ANGLE_CLASS     ${extension}  ANGLE     ${sources})
+    FindStyleHeadersExt(${search_path} ANGLE_CLASS        ${extension}  ANGLE        ${sources})
-    FindStyleHeadersExt(${search_path} ATOM_CLASS      ${extension}  ATOM_VEC  ${sources})
+    FindStyleHeadersExt(${search_path} ATOM_CLASS         ${extension}  ATOM_VEC     ${sources})
-    FindStyleHeadersExt(${search_path} BODY_CLASS      ${extension}  BODY      ${sources})
+    FindStyleHeadersExt(${search_path} BODY_CLASS         ${extension}  BODY         ${sources})
-    FindStyleHeadersExt(${search_path} BOND_CLASS      ${extension}  BOND      ${sources})
+    FindStyleHeadersExt(${search_path} BOND_CLASS         ${extension}  BOND         ${sources})
-    FindStyleHeadersExt(${search_path} COMMAND_CLASS   ${extension}  COMMAND   ${sources})
+    FindStyleHeadersExt(${search_path} COMMAND_CLASS      ${extension}  COMMAND      ${sources})
-    FindStyleHeadersExt(${search_path} COMPUTE_CLASS   ${extension}  COMPUTE   ${sources})
+    FindStyleHeadersExt(${search_path} COMPUTE_CLASS      ${extension}  COMPUTE      ${sources})
-    FindStyleHeadersExt(${search_path} DIHEDRAL_CLASS  ${extension}  DIHEDRAL  ${sources})
+    FindStyleHeadersExt(${search_path} DIHEDRAL_CLASS     ${extension}  DIHEDRAL     ${sources})
-    FindStyleHeadersExt(${search_path} DUMP_CLASS      ${extension}  DUMP      ${sources})
+    FindStyleHeadersExt(${search_path} DUMP_CLASS         ${extension}  DUMP         ${sources})
-    FindStyleHeadersExt(${search_path} FIX_CLASS       ${extension}  FIX       ${sources})
+    FindStyleHeadersExt(${search_path} FIX_CLASS          ${extension}  FIX          ${sources})
-    FindStyleHeadersExt(${search_path} IMPROPER_CLASS  ${extension}  IMPROPER  ${sources})
+    FindStyleHeadersExt(${search_path} GRAN_SUB_MOD_CLASS ${extension}  GRAN_SUB_MOD ${sources})
-    FindStyleHeadersExt(${search_path} INTEGRATE_CLASS ${extension}  INTEGRATE ${sources})
+    FindStyleHeadersExt(${search_path} IMPROPER_CLASS     ${extension}  IMPROPER     ${sources})
-    FindStyleHeadersExt(${search_path} KSPACE_CLASS    ${extension}  KSPACE    ${sources})
+    FindStyleHeadersExt(${search_path} INTEGRATE_CLASS    ${extension}  INTEGRATE    ${sources})
-    FindStyleHeadersExt(${search_path} MINIMIZE_CLASS  ${extension}  MINIMIZE  ${sources})
+    FindStyleHeadersExt(${search_path} KSPACE_CLASS       ${extension}  KSPACE       ${sources})
-    FindStyleHeadersExt(${search_path} NBIN_CLASS      ${extension}  NBIN      ${sources})
+    FindStyleHeadersExt(${search_path} MINIMIZE_CLASS     ${extension}  MINIMIZE     ${sources})
-    FindStyleHeadersExt(${search_path} NPAIR_CLASS     ${extension}  NPAIR     ${sources})
+    FindStyleHeadersExt(${search_path} NBIN_CLASS         ${extension}  NBIN         ${sources})
-    FindStyleHeadersExt(${search_path} NSTENCIL_CLASS  ${extension}  NSTENCIL  ${sources})
+    FindStyleHeadersExt(${search_path} NPAIR_CLASS        ${extension}  NPAIR        ${sources})
-    FindStyleHeadersExt(${search_path} NTOPO_CLASS     ${extension}  NTOPO     ${sources})
+    FindStyleHeadersExt(${search_path} NSTENCIL_CLASS     ${extension}  NSTENCIL     ${sources})
-    FindStyleHeadersExt(${search_path} PAIR_CLASS      ${extension}  PAIR      ${sources})
+    FindStyleHeadersExt(${search_path} NTOPO_CLASS        ${extension}  NTOPO        ${sources})
-    FindStyleHeadersExt(${search_path} READER_CLASS    ${extension}  READER    ${sources})
+    FindStyleHeadersExt(${search_path} PAIR_CLASS         ${extension}  PAIR         ${sources})
-    FindStyleHeadersExt(${search_path} REGION_CLASS    ${extension}  REGION    ${sources})
+    FindStyleHeadersExt(${search_path} READER_CLASS       ${extension}  READER       ${sources})
+    FindStyleHeadersExt(${search_path} REGION_CLASS       ${extension}  REGION       ${sources})
 endfunction(RegisterStylesExt)
 
 function(GenerateStyleHeaders output_path)
     message(STATUS "Generating style headers...")
-    GenerateStyleHeader(${output_path} ANGLE      angle     ) # force
+    GenerateStyleHeader(${output_path} ANGLE        angle        ) # force
-    GenerateStyleHeader(${output_path} ATOM_VEC   atom      ) # atom      atom_vec_hybrid
+    GenerateStyleHeader(${output_path} ATOM_VEC     atom         ) # atom      atom_vec_hybrid
-    GenerateStyleHeader(${output_path} BODY       body      ) # atom_vec_body
+    GenerateStyleHeader(${output_path} BODY         body         ) # atom_vec_body
-    GenerateStyleHeader(${output_path} BOND       bond      ) # force
+    GenerateStyleHeader(${output_path} BOND         bond         ) # force
-    GenerateStyleHeader(${output_path} COMMAND    command   ) # input
+    GenerateStyleHeader(${output_path} COMMAND      command      ) # input
-    GenerateStyleHeader(${output_path} COMPUTE    compute   ) # modify
+    GenerateStyleHeader(${output_path} COMPUTE      compute      ) # modify
-    GenerateStyleHeader(${output_path} DIHEDRAL   dihedral  ) # force
+    GenerateStyleHeader(${output_path} DIHEDRAL     dihedral     ) # force
-    GenerateStyleHeader(${output_path} DUMP       dump      ) # output    write_dump
+    GenerateStyleHeader(${output_path} DUMP         dump         ) # output    write_dump
-    GenerateStyleHeader(${output_path} FIX        fix       ) # modify
+    GenerateStyleHeader(${output_path} FIX          fix          ) # modify
-    GenerateStyleHeader(${output_path} IMPROPER   improper  ) # force
+    GenerateStyleHeader(${output_path} GRAN_SUB_MOD gran_sub_mod ) # granular_model
-    GenerateStyleHeader(${output_path} INTEGRATE  integrate ) # update
+    GenerateStyleHeader(${output_path} IMPROPER     improper     ) # force
-    GenerateStyleHeader(${output_path} KSPACE     kspace    ) # force
+    GenerateStyleHeader(${output_path} INTEGRATE    integrate    ) # update
-    GenerateStyleHeader(${output_path} MINIMIZE   minimize  ) # update
+    GenerateStyleHeader(${output_path} KSPACE       kspace       ) # force
-    GenerateStyleHeader(${output_path} NBIN       nbin      ) # neighbor
+    GenerateStyleHeader(${output_path} MINIMIZE     minimize     ) # update
-    GenerateStyleHeader(${output_path} NPAIR      npair     ) # neighbor
+    GenerateStyleHeader(${output_path} NBIN         nbin         ) # neighbor
-    GenerateStyleHeader(${output_path} NSTENCIL   nstencil  ) # neighbor
+    GenerateStyleHeader(${output_path} NPAIR        npair        ) # neighbor
-    GenerateStyleHeader(${output_path} NTOPO      ntopo     ) # neighbor
+    GenerateStyleHeader(${output_path} NSTENCIL     nstencil     ) # neighbor
-    GenerateStyleHeader(${output_path} PAIR       pair      ) # force
+    GenerateStyleHeader(${output_path} NTOPO        ntopo        ) # neighbor
-    GenerateStyleHeader(${output_path} READER     reader    ) # read_dump
+    GenerateStyleHeader(${output_path} PAIR         pair         ) # force
-    GenerateStyleHeader(${output_path} REGION     region    ) # domain
+    GenerateStyleHeader(${output_path} READER       reader       ) # read_dump
+    GenerateStyleHeader(${output_path} REGION       region       ) # domain
 endfunction(GenerateStyleHeaders)
 
 function(DetectBuildSystemConflict lammps_src_dir)

doc/Makefile

@@ -94,10 +94,11 @@ html: xmlgen $(VENV) $(SPHINXCONFIG)/conf.py $(ANCHORCHECK) $(MATHJAX)
 		rm -f $(BUILDDIR)/doxygen/xml/run.stamp;\
 		echo "############################################" ; env PYTHONWARNINGS= PYTHONDONTWRITEBYTECODE=1 \
 		rst_anchor_check src/*.rst ;\
-		python $(BUILDDIR)/utils/check-packages.py -s ../src -d src ;\
+		$(PYTHON) $(BUILDDIR)/utils/check-packages.py -s ../src -d src ;\
 		env LC_ALL=C grep -n '[^ -~]' $(RSTDIR)/*.rst ;\
 		env LC_ALL=C grep -n ' :[a-z]\+`' $(RSTDIR)/*.rst ;\
-		python $(BUILDDIR)/utils/check-styles.py -s ../src -d src ;\
+		env LC_ALL=C grep -n ' `[^`]\+<[a-z][^`]\+`[^_]' $(RSTDIR)/*.rst ;\
+		$(PYTHON) $(BUILDDIR)/utils/check-styles.py -s ../src -d src ;\
 		echo "############################################" ;\
 		deactivate ;\
 	)
@@ -174,10 +175,11 @@ pdf: xmlgen $(VENV) $(SPHINXCONFIG)/conf.py $(ANCHORCHECK)
 		rm -f $(BUILDDIR)/doxygen/xml/run.stamp;\
 		echo "############################################" ; env PYTHONWARNINGS= PYTHONDONTWRITEBYTECODE=1 \
 		rst_anchor_check src/*.rst ;\
-		python utils/check-packages.py -s ../src -d src ;\
+		$(PYTHON) utils/check-packages.py -s ../src -d src ;\
 		env LC_ALL=C grep -n '[^ -~]' $(RSTDIR)/*.rst ;\
 		env LC_ALL=C grep -n ' :[a-z]\+`' $(RSTDIR)/*.rst ;\
-		python utils/check-styles.py -s ../src -d src ;\
+		env LC_ALL=C grep -n ' `[^`]\+<[a-z][^`]\+`[^_]' $(RSTDIR)/*.rst ;\
+		$(PYTHON) utils/check-styles.py -s ../src -d src ;\
 		echo "############################################" ;\
 		deactivate ;\
 	)
@@ -208,14 +210,14 @@ anchor_check : $(ANCHORCHECK)
 style_check : $(VENV)
 	@(\
 		. $(VENV)/bin/activate ; env PYTHONWARNINGS= PYTHONDONTWRITEBYTECODE=1 \
-		python utils/check-styles.py -s ../src -d src ;\
+		$(PYTHON) utils/check-styles.py -s ../src -d src ;\
 		deactivate ;\
 	)
 
 package_check : $(VENV)
 	@(\
-		. $(VENV)/bin/activate ; env PYTHONWARNINGS= PYTHONDONTWRITEBYTECODE=1 \\
+		. $(VENV)/bin/activate ; env PYTHONWARNINGS= PYTHONDONTWRITEBYTECODE=1 \
-		python utils/check-packages.py -s ../src -d src ;\
+		$(PYTHON) utils/check-packages.py -s ../src -d src ;\
 		deactivate ;\
 	)
 
@@ -224,6 +226,7 @@ char_check :
 
 role_check :
 	@( env LC_ALL=C grep -n ' :[a-z]\+`' $(RSTDIR)/*.rst && exit 1 || : )
+	@( env LC_ALL=C grep -n ' `[^`]\+<[a-z][^`]\+`[^_]' $(RSTDIR)/*.rst && exit 1 || : )
 
 link_check : $(VENV) html
 	@(\

doc/github-development-workflow.md

@@ -123,7 +123,7 @@ request is merged.  The template for pull requests includes a header
 where connections between pull requests and issues can be listed, and
 thus where this comment should be placed.
 
-## Milestones and Release Planning
+## Release Planning
 
 LAMMPS uses a continuous release development model with incremental
 changes, i.e. significant effort is made -- including automated pre-merge
@@ -132,24 +132,28 @@ broken.  These tests are run after every update to a pull request.  More
 extensive and time-consuming tests (including regression testing) are
 performed after code is merged to the "develop" branch.  There are feature
 releases of LAMMPS made about every 4-6 weeks at a point, when the LAMMPS
-developers feel, that a sufficient number of changes has been included
+developers feel, that a sufficient number of changes have been included
 and all post-merge testing has been successful.  These feature releases are
 marked with a `patch_<version date>` tag and the "release" branch
 follows only these versions with fast-forward merges.  While "develop" may
 be temporarily broken through issues only detected by the post-merge tests,
 The "release" branch is always supposed to be of production quality.
 
-About once each year, there is a "stable" release of LAMMPS.
+About once each year, there is a "stable" release of LAMMPS.  These have
-These have seen additional, manual testing and review of
+seen additional, manual testing and review of results from testing with
-results from testing with instrumented code and static code analysis.
+instrumented code and static code analysis.  Also, the last few feature
-Also, the last few feature releases before a stable release are "release
+releases before a stable release are "release candidate" versions which
-candidate" versions which only contain bug fixes, feature additions to
+only contain bug fixes, feature additions to peripheral functionality,
-peripheral functionality, and documentation updates.  In between stable
+and documentation updates.  In between stable releases, bug fixes and
-releases, bug fixes and infrastructure updates are back-ported from the
+infrastructure updates are back-ported from the "develop" branch to the
-"develop" branch to the "maintenance" branch and occasionally merged
+"maintenance" branch and occasionally merged into "stable" and published
-into "stable" and published as update releases.
+as update releases.
+
+## Project Management
 
 For release planning and the information of code contributors, issues
-and pull requests being actively worked on are assigned a "milestone",
+and pull requests are being managed with GitHub Project Boards.  There
-which corresponds to the next stable release or the stable release after
+are currently three boards: LAMMPS Feature Requests, LAMMPS Bug Reports,
-that, with a tentative release date.
+and LAMMPS Pull Requests.  Each board is organized in columns where
+submissions are categorized.  Within each column the entries are
+(manually) sorted according their priority.

doc/lammps.1

@@ -1,7 +1,7 @@
-.TH LAMMPS "1" "8 February 2023" "2023-02-08"
+.TH LAMMPS "1" "28 March 2023" "2023-03-28"
 .SH NAME
 .B LAMMPS
-\- Molecular Dynamics Simulator.  Version 8 February 2023
+\- Molecular Dynamics Simulator.  Version 28 March 2023
 
 .SH SYNOPSIS
 .B lmp

doc/msi2lmp.1

@@ -1,11 +1,11 @@
-.TH MSI2LMP "1" "v3.9.9" "2018-11-05"
+.TH MSI2LMP "1" "v3.9.10" "2023-03-10"
 .SH NAME
 .B MSI2LMP
 \- Converter for Materials Studio files to LAMMPS
 
 .SH SYNOPSIS
 .B msi2lmp
-<ROOTNAME> [-class <I|1|II|2|O|0>] [-frc <path to frc file>] [-print #] [-ignore] [-nocenter] [-oldstyle] [-shift <x> <y> <z>]
+[-help] <ROOTNAME> [-class <I|1|II|2|O|0>] [-frc <path to frc file>] [-print #] [-ignore] [-nocenter] [-oldstyle] [-shift <x> <y> <z>]
 
 .SH DESCRIPTION
 .PP
@@ -22,6 +22,9 @@ needed between .frc and .car/.mdf files are the atom types.
 
 .SH OPTIONS
 .TP
+\fB\-h\fR, \fB\-help\fR,
+Print detailed help message to the screen and stop.
+.TP
 \fB\<ROOTNAME>\fR
 This has to be the first argument and is a
 .B mandatory

doc/src/Build_development.rst

@@ -523,6 +523,8 @@ The following options are available.
 These should help to make source and documentation files conforming
 to some the coding style preferences of the LAMMPS developers.
 
+.. _clang-format:
+
 Clang-format support
 --------------------
 

doc/src/Build_extras.rst

@@ -180,6 +180,9 @@ way no local OpenCL development headers or library needs to be present and only
 OpenCL compatible drivers need to be installed to use OpenCL.  If this is not
 desired, you can set :code:`USE_STATIC_OPENCL_LOADER` to :code:`no`.
 
+The GPU library has some multi-thread support using OpenMP.  If LAMMPS is built
+with ``-D BUILD_OMP=on`` this will also be enabled.
+
 If you are compiling with HIP, note that before running CMake you will have to
 set appropriate environment variables. Some variables such as
 :code:`HCC_AMDGPU_TARGET` (for ROCm <= 4.0) or :code:`CUDA_PATH` are necessary for :code:`hipcc`
@@ -273,10 +276,13 @@ To enable GPU binning via CUDA performance primitives set the Makefile variable
 most modern GPUs.
 
 To support the CUDA multiprocessor server you can set the define
-``-DCUDA_PROXY``.  Please note that in this case you must **not** use
+``-DCUDA_MPS_SUPPORT``.  Please note that in this case you must **not** use
 the CUDA performance primitives and thus set the variable ``CUDPP_OPT``
 to empty.
 
+The GPU library has some multi-thread support using OpenMP.  You need to add
+the compiler flag that enables OpenMP to the ``CUDR_OPTS`` Makefile variable.
+
 If the library build is successful, 3 files should be created:
 ``lib/gpu/libgpu.a``\ , ``lib/gpu/nvc_get_devices``\ , and
 ``lib/gpu/Makefile.lammps``\ .  The latter has settings that enable LAMMPS
@@ -1906,10 +1912,10 @@ OPENMP package
    Apple offers the `Xcode package and IDE
    <https://developer.apple.com/xcode/>`_ for compiling software on
    macOS, so you have likely installed it to compile LAMMPS.  Their
-   compiler is based on `Clang <https://clang.llvm.org/>`, but while it
+   compiler is based on `Clang <https://clang.llvm.org/>`_, but while it
    is capable of processing OpenMP directives, the necessary header
    files and OpenMP runtime library are missing.  The `R developers
-   <https://www.r-project.org/>` have figured out a way to build those
+   <https://www.r-project.org/>`_ have figured out a way to build those
    in a compatible fashion. One can download them from
    `https://mac.r-project.org/openmp/
    <https://mac.r-project.org/openmp/>`_.  Simply adding those files as

doc/src/Commands_bond.rst

@@ -42,6 +42,7 @@ OPT.
    * :doc:`gaussian <bond_gaussian>`
    * :doc:`gromos (o) <bond_gromos>`
    * :doc:`harmonic (iko) <bond_harmonic>`
+   * :doc:`harmonic/restrain <bond_harmonic_restrain>`
    * :doc:`harmonic/shift (o) <bond_harmonic_shift>`
    * :doc:`harmonic/shift/cut (o) <bond_harmonic_shift_cut>`
    * :doc:`lepton (o) <bond_lepton>`

doc/src/Commands_compute.rst

@@ -52,6 +52,8 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`dilatation/atom <compute_dilatation_atom>`
    * :doc:`dipole <compute_dipole>`
    * :doc:`dipole/chunk <compute_dipole_chunk>`
+   * :doc:`dipole/tip4p <compute_dipole>`
+   * :doc:`dipole/tip4p/chunk <compute_dipole_chunk>`
    * :doc:`displace/atom <compute_displace_atom>`
    * :doc:`dpd <compute_dpd>`
    * :doc:`dpd/atom <compute_dpd_atom>`
@@ -104,6 +106,7 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`pe/tally <compute_tally>`
    * :doc:`plasticity/atom <compute_plasticity_atom>`
    * :doc:`pressure <compute_pressure>`
+   * :doc:`pressure/alchemy <compute_pressure_alchemy>`
    * :doc:`pressure/uef <compute_pressure_uef>`
    * :doc:`property/atom <compute_property_atom>`
    * :doc:`property/chunk <compute_property_chunk>`

doc/src/Commands_fix.rst

@@ -29,6 +29,7 @@ OPT.
    * :doc:`adapt/fep <fix_adapt_fep>`
    * :doc:`addforce <fix_addforce>`
    * :doc:`addtorque <fix_addtorque>`
+   * :doc:`alchemy <fix_alchemy>`
    * :doc:`amoeba/bitorsion <fix_amoeba_bitorsion>`
    * :doc:`amoeba/pitorsion <fix_amoeba_pitorsion>`
    * :doc:`append/atoms <fix_append_atoms>`
@@ -69,6 +70,7 @@ OPT.
    * :doc:`dt/reset (k) <fix_dt_reset>`
    * :doc:`edpd/source <fix_dpd_source>`
    * :doc:`efield <fix_efield>`
+   * :doc:`efield/tip4p <fix_efield>`
    * :doc:`ehex <fix_ehex>`
    * :doc:`electrode/conp (i) <fix_electrode>`
    * :doc:`electrode/conq (i) <fix_electrode>`
@@ -92,6 +94,7 @@ OPT.
    * :doc:`grem <fix_grem>`
    * :doc:`halt <fix_halt>`
    * :doc:`heat <fix_heat>`
+   * :doc:`heat/flow <fix_heat_flow>`
    * :doc:`hyper/global <fix_hyper_global>`
    * :doc:`hyper/local <fix_hyper_local>`
    * :doc:`imd <fix_imd>`
@@ -107,6 +110,7 @@ OPT.
    * :doc:`lineforce <fix_lineforce>`
    * :doc:`manifoldforce <fix_manifoldforce>`
    * :doc:`mdi/qm <fix_mdi_qm>`
+   * :doc:`mdi/qmmm <fix_mdi_qmmm>`
    * :doc:`meso/move <fix_meso_move>`
    * :doc:`mol/swap <fix_mol_swap>`
    * :doc:`momentum (k) <fix_momentum>`
@@ -262,6 +266,7 @@ OPT.
    * :doc:`wall/lj1043 <fix_wall>`
    * :doc:`wall/lj126 <fix_wall>`
    * :doc:`wall/lj93 (k) <fix_wall>`
+   * :doc:`wall/lepton <fix_wall>`
    * :doc:`wall/morse <fix_wall>`
    * :doc:`wall/piston <fix_wall_piston>`
    * :doc:`wall/reflect (k) <fix_wall_reflect>`
@@ -269,4 +274,5 @@ OPT.
    * :doc:`wall/region <fix_wall_region>`
    * :doc:`wall/region/ees <fix_wall_ees>`
    * :doc:`wall/srd <fix_wall_srd>`
+   * :doc:`wall/table <fix_wall>`
    * :doc:`widom <fix_widom>`

doc/src/Commands_pair.rst

@@ -55,6 +55,7 @@ OPT.
    * :doc:`born/coul/msm (o) <pair_born>`
    * :doc:`born/coul/wolf (go) <pair_born>`
    * :doc:`born/coul/wolf/cs (g) <pair_cs>`
+   * :doc:`born/gauss <pair_born_gauss>`
    * :doc:`bpm/spring <pair_bpm_spring>`
    * :doc:`brownian (o) <pair_brownian>`
    * :doc:`brownian/poly (o) <pair_brownian>`

doc/src/Developer_notes.rst

@@ -11,6 +11,7 @@ Available topics are:
 
 - `Reading and parsing of text and text files`_
 - `Requesting and accessing neighbor lists`_
+- `Choosing between a custom atom style, fix property/atom, and fix STORE/ATOM`_
 - `Fix contributions to instantaneous energy, virial, and cumulative energy`_
 - `KSpace PPPM FFT grids`_
 
@@ -73,6 +74,8 @@ when converting "12.5", while the ValueTokenizer class will throw an
 :cpp:func:`ValueTokenizer::next_int()
 <LAMMPS_NS::ValueTokenizer::next_int>` is called on the same string.
 
+.. _request-neighbor-list:
+
 Requesting and accessing neighbor lists
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -216,6 +219,30 @@ command:
 
    neighbor->add_request(this, "delete_atoms", NeighConst::REQ_FULL);
 
+Choosing between a custom atom style, fix property/atom, and fix STORE/ATOM
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+There are multiple ways to manage per-atom data within LAMMPS.  Often
+the per-atom storage is only used locally and managed by the class that
+uses it.  If the data has to persist between multiple time steps and
+migrate with atoms when they move from sub-domain to sub-domain or
+across periodic boundaries, then using a custom atom style, or :doc:`fix
+property/atom <fix_property_atom>`, or the internal fix STORE/ATOM are
+possible options.
+
+- Using the atom style is usually the most programming effort and mostly
+  needed when the per-atom data is an integral part of the model like a
+  per-atom charge or diameter and thus should be part of the Atoms
+  section of a :doc:`data file <read_data>`.
+
+- Fix property/atom is useful if the data is optional or should be
+  entered by the user, or accessed as a (named) custom property. In this
+  case the fix should be entered as part of the input (and not
+  internally) which allows to enter and store its content with data files.
+
+- Fix STORE/ATOM should be used when the data should be accessed internally
+  only and thus the fix can be created internally.
+
 Fix contributions to instantaneous energy, virial, and cumulative energy
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

doc/src/Developer_updating.rst

@@ -24,6 +24,7 @@ Available topics in mostly chronological order are:
 - `Use of "override" instead of "virtual"`_
 - `Simplified and more compact neighbor list requests`_
 - `Split of fix STORE into fix STORE/GLOBAL and fix STORE/PERATOM`_
+- `Rename of fix STORE/PERATOM to fix STORE/ATOM and change of arguments`_
 - `Use Output::get_dump_by_id() instead of Output::find_dump()`_
 - `Refactored grid communication using Grid3d/Grid2d classes instead of GridComm`_
 
@@ -385,6 +386,34 @@ New:
 
 This change is **required** or else the code will not compile.
 
+Rename of fix STORE/PERATOM to fix STORE/ATOM and change of arguments
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. versionchanged:: 28Mar2023
+
+The available functionality of the internal fix to store per-atom
+properties was expanded to enable storing data with ghost atoms and to
+support binary restart files.  With those changes, the fix was renamed
+to fix STORE/ATOM and the number and order of (required) arguments has
+changed.
+
+Old syntax: ``ID group-ID STORE/PERATOM rflag n1 n2 [n3]``
+
+- *rflag* = 0/1, *no*/*yes* store per-atom values in restart file
+- :math:`n1 = 1, n2 = 1, \mathrm{no}\;n3 \to` per-atom vector, single value per atom
+- :math:`n1 = 1, n2 > 1, \mathrm{no}\;n3 \to` per-atom array, *n2* values per atom
+- :math:`n1 = 1, n2 > 0, n3 > 0 \to` per-atom tensor, *n2* x *n3* values per atom
+
+New syntax:  ``ID group-ID STORE/ATOM n1 n2 gflag rflag``
+
+- :math:`n1 = 1, n2 = 0 \to` per-atom vector, single value per atom
+- :math:`n1 > 1, n2 = 0 \to` per-atom array, *n1* values per atom
+- :math:`n1 > 0, n2 > 0 \to` per-atom tensor, *n1* x *n2* values per atom
+- *gflag* = 0/1, *no*/*yes* communicate per-atom values with ghost atoms
+- *rflag* = 0/1, *no*/*yes* store per-atom values in restart file
+
+Since this is an internal fix, there is no user visible change.
+
 Use Output::get_dump_by_id() instead of Output::find_dump()
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

doc/src/Developer_write.rst

@@ -6,250 +6,8 @@ be extended by writing new classes that derive from existing
 parent classes in LAMMPS.  Here, some specific coding
 details are provided for writing code for LAMMPS.
 
-Writing a new fix style
+.. toctree::
-^^^^^^^^^^^^^^^^^^^^^^^
+   :maxdepth: 1
-
-Writing fixes is a flexible way of extending LAMMPS.  Users can
-implement many things using fixes:
-
-- changing particles attributes (positions, velocities, forces, etc.). Examples: FixNVE, FixFreeze.
-- reading/writing data. Example: FixRestart.
-- adding or modifying properties due to geometry. Example: FixWall.
-- interacting with other subsystems or external code: Examples: FixTTM, FixExternal
-- saving information for analysis or future use (previous positions,
-  for instance). Examples: Fix AveTime, FixStoreState.
-
-
-All fixes are derived from the Fix base class and must have a
-constructor with the signature: ``FixPrintVel(class LAMMPS *, int, char **)``.
-
-Every fix must be registered in LAMMPS by writing the following lines
-of code in the header before include guards:
-
-.. code-block:: c++
-
-   #ifdef FIX_CLASS
-   // clang-format off
-   FixStyle(print/vel,FixPrintVel);
-   // clang-format on
-   #else
-   /* the definition of the FixPrintVel class comes here */
-   ...
-   #endif
-
-Where ``print/vel`` is the style name of your fix in the input script and
-``FixPrintVel`` is the name of the class. The header file would be called
-``fix_print_vel.h`` and the implementation file ``fix_print_vel.cpp``.
-These conventions allow LAMMPS to automatically integrate it into the
-executable when compiling and associate your new fix class with the designated
-keyword when it parses the input script.
-
-Let's write a simple fix which will print the average velocity at the end
-of each timestep. First of all, implement a constructor:
-
-.. code-block:: c++
-
-   FixPrintVel::FixPrintVel(LAMMPS *lmp, int narg, char **arg)
-   : Fix(lmp, narg, arg)
-   {
-     if (narg < 4)
-       error->all(FLERR,"Illegal fix print/vel command");
-
-     nevery = utils::inumeric(FLERR,arg[3],false,lmp);
-     if (nevery <= 0)
-       error->all(FLERR,"Illegal fix print/vel command");
-   }
-
-In the constructor you should parse your fix arguments which are
-specified in the script. All fixes have pretty much the same syntax:
-``fix <fix-ID> <fix group> <fix name> <fix arguments ...>``. The
-first 3 parameters are parsed by Fix base class constructor, while
-``<fix arguments>`` should be parsed by you. In our case, we need to
-specify how often we want to print an average velocity. For instance,
-once in 50 timesteps: ``fix 1 print/vel 50``. There is a special variable
-in the Fix class called ``nevery`` which specifies how often the method
-``end_of_step()`` is called. Thus all we need to do is just set it up.
-
-The next method we need to implement is ``setmask()``:
-
-.. code-block:: c++
-
-   int FixPrintVel::setmask()
-   {
-     int mask = 0;
-     mask |= FixConst::END_OF_STEP;
-     return mask;
-   }
-
-Here the user specifies which methods of your fix should be called
-during execution. The constant ``END_OF_STEP`` corresponds to the
-``end_of_step()`` method. The most important available methods that
-are called during a timestep and the order in which they are called
-are shown in the previous section.
-
-.. code-block:: c++
-
-   void FixPrintVel::end_of_step()
-   {
-     // for add3, scale3
-     using namespace MathExtra;
-
-     double** v = atom->v;
-     int nlocal = atom->nlocal;
-     double localAvgVel[4]; // 4th element for particles count
-     memset(localAvgVel, 0, 4 * sizeof(double));
-     for (int particleInd = 0; particleInd < nlocal; ++particleInd) {
-       add3(localAvgVel, v[particleInd], localAvgVel);
-     }
-     localAvgVel[3] = nlocal;
-     double globalAvgVel[4];
-     memset(globalAvgVel, 0, 4 * sizeof(double));
-     MPI_Allreduce(localAvgVel, globalAvgVel, 4, MPI_DOUBLE, MPI_SUM, world);
-     scale3(1.0 / globalAvgVel[3], globalAvgVel);
-     if ((comm->me == 0) && screen) {
-       fmt::print(screen,"{}, {}, {}\n",
-                  globalAvgVel[0], globalAvgVel[1], globalAvgVel[2]);
-     }
-   }
-
-In the code above, we use MathExtra routines defined in
-``math_extra.h``.  There are bunch of math functions to work with
-arrays of doubles as with math vectors.  It is also important to note
-that LAMMPS code should always assume to be run in parallel and that
-atom data is thus distributed across the MPI ranks.  Thus you can
-only process data from local atoms directly and need to use MPI library
-calls to combine or exchange data.  For serial execution, LAMMPS
-comes bundled with the MPI STUBS library that contains the MPI library
-function calls in dummy versions that only work for a single MPI rank.
-
-In this code we use an instance of Atom class. This object is stored
-in the Pointers class (see ``pointers.h``) which is the base class of
-the Fix base class. This object contains references to various class
-instances (the original instances are created and held by the LAMMPS
-class) with all global information about the simulation system.
-Data from the Pointers class is available to all classes inherited from
-it using protected inheritance. Hence when you write you own class,
-which is going to use LAMMPS data, don't forget to inherit from Pointers
-or pass an Pointer to it to all functions that need access. When writing
-fixes we inherit from class Fix which is inherited from Pointers so
-there is no need to inherit from it directly.
-
-The code above computes average velocity for all particles in the
-simulation.  Yet you have one unused parameter in fix call from the
-script: ``group_name``.  This parameter specifies the group of atoms
-used in the fix. So we should compute average for all particles in the
-simulation only if ``group_name == "all"``, but it can be any group.
-The group membership information of an atom is contained in the *mask*
-property of and atom and the bit corresponding to a given group is
-stored in the groupbit variable which is defined in Fix base class:
-
-.. code-block:: c++
-
-   for (int i = 0; i < nlocal; ++i) {
-     if (atom->mask[i] & groupbit) {
-     // Do all job here
-     }
-   }
-
-Class Atom encapsulates atoms positions, velocities, forces, etc. User
-can access them using particle index. Note, that particle indexes are
-usually changed every few timesteps because of neighbor list rebuilds
-and spatial sorting (to improve cache efficiency).
-
-Let us consider another Fix example: We want to have a fix which stores
-atoms position from previous time step in your fix. The local atoms
-indexes may not be valid on the next iteration. In order to handle
-this situation there are several methods which should be implemented:
-
-- ``double memory_usage()``: return how much memory the fix uses (optional)
-- ``void grow_arrays(int)``: do reallocation of the per particle arrays in your fix
-- ``void copy_arrays(int i, int j, int delflag)``: copy i-th per-particle
-  information to j-th. Used when atom sorting is performed. if delflag is set
-  and atom j owns a body, move the body information to atom i.
-- ``void set_arrays(int i)``: sets i-th particle related information to zero
-
-Note, that if your class implements these methods, it must call add calls of
-add_callback and delete_callback to constructor and destructor. Since we want
-to store positions of atoms from previous timestep, we need to add
-``double** xold`` to the header file. Than add allocation code
-to the constructor:
-
-.. code-block:: c++
-
-   FixSavePos::FixSavePos(LAMMPS *lmp, int narg, char **arg), xold(nullptr)
-   {
-   //...
-     memory->create(xold, atom->nmax, 3, "FixSavePos:x");
-     atom->add_callback(0);
-   }
-
-   FixSavePos::~FixSavePos() {
-     atom->delete_callback(id, 0);
-     memory->destroy(xold);
-   }
-
-Implement the aforementioned methods:
-
-.. code-block:: c++
-
-   double FixSavePos::memory_usage()
-   {
-     int nmax = atom->nmax;
-     double bytes = 0.0;
-     bytes += nmax * 3 * sizeof(double);
-     return bytes;
-   }
-
-   void FixSavePos::grow_arrays(int nmax)
-   {
-     memory->grow(xold, nmax, 3, "FixSavePos:xold");
-   }
-
-   void FixSavePos::copy_arrays(int i, int j, int delflag)
-   {
-     memcpy(xold[j], xold[i], sizeof(double) * 3);
-   }
-
-   void FixSavePos::set_arrays(int i)
-   {
-     memset(xold[i], 0, sizeof(double) * 3);
-   }
-
-   int FixSavePos::pack_exchange(int i, double *buf)
-   {
-     int m = 0;
-     buf[m++] = xold[i][0];
-     buf[m++] = xold[i][1];
-     buf[m++] = xold[i][2];
-
-     return m;
-   }
-
-   int FixSavePos::unpack_exchange(int nlocal, double *buf)
-   {
-     int m = 0;
-     xold[nlocal][0] = buf[m++];
-     xold[nlocal][1] = buf[m++];
-     xold[nlocal][2] = buf[m++];
-
-     return m;
-   }
-
-Now, a little bit about memory allocation. We use the Memory class which
-is just a bunch of template functions for allocating 1D and 2D
-arrays. So you need to add include ``memory.h`` to have access to them.
-
-Finally, if you need to write/read some global information used in
-your fix to the restart file, you might do it by setting flag
-``restart_global = 1`` in the constructor and implementing methods void
-``write_restart(FILE *fp)`` and ``void restart(char *buf)``.
-If, in addition, you want to write the per-atom property to restart
-files additional settings and functions are needed:
-
-- a fix flag indicating this needs to be set ``restart_peratom = 1;``
-- ``atom->add_callback()`` and ``atom->delete_callback()`` must be called
-  a second time with the final argument set to 1 instead of 0 (indicating
-  restart processing instead of per-atom data memory management).
-- the functions ``void pack_restart(int i, double *buf)`` and
-  ``void unpack_restart(int nlocal, int nth)`` need to be implemented
 
+   Developer_write_pair
+   Developer_write_fix

doc/src/Developer_write_fix.rst

@@ -0,0 +1,245 @@
+Writing a new fix style
+^^^^^^^^^^^^^^^^^^^^^^^
+
+Writing fix styles is a flexible way of extending LAMMPS.  Users can
+implement many things using fixes.  Some fix styles are only used
+internally to support compute styles or pair styles:
+
+- change particles attributes (positions, velocities, forces, etc.). Examples: ``FixNVE``, ``FixFreeze``.
+- read or write data.  Example: ``FixRestart``.
+- adding or modifying properties due to geometry. Example: ``FixWall``.
+- interacting with other subsystems or external code: Examples: ``FixTTM``, ``FixExternal``, ``FixMDI``
+- saving information for analysis or future use (previous positions,
+  for instance). Examples: ``FixAveTime``, ``FixStoreState``.
+
+All fixes are derived from the ``Fix`` base class and must have a
+constructor with the signature: ``FixPrintVel(class LAMMPS *, int, char **)``.
+
+Every fix must be registered in LAMMPS by writing the following lines
+of code in the header before include guards:
+
+.. code-block:: c++
+
+   #ifdef FIX_CLASS
+   // clang-format off
+   FixStyle(print/vel,FixPrintVel);
+   // clang-format on
+   #else
+   /* the definition of the FixPrintVel class comes here */
+   ...
+   #endif
+
+Where ``print/vel`` is the style name of your fix in the input script and
+``FixPrintVel`` is the name of the class. The header file would be called
+``fix_print_vel.h`` and the implementation file ``fix_print_vel.cpp``.
+These conventions allow LAMMPS to automatically integrate it into the
+executable when compiling and associate your new fix class with the designated
+keyword when it parses the input script.
+
+Let's write a simple fix which will print the average velocity at the end
+of each timestep. First of all, implement a constructor:
+
+.. code-block:: c++
+
+   FixPrintVel::FixPrintVel(LAMMPS *lmp, int narg, char **arg)
+   : Fix(lmp, narg, arg)
+   {
+     if (narg < 4) utils::missing_cmd_args(FLERR, "fix print/vel", error);
+
+     nevery = utils::inumeric(FLERR,arg[3],false,lmp);
+     if (nevery <= 0)
+       error->all(FLERR,"Illegal fix print/vel nevery value: {}", nevery);
+   }
+
+In the constructor you should parse the fix arguments which are
+specified in the script.  All fixes have pretty much the same syntax:
+``fix <fix-ID> <fix group> <fix name> <fix arguments ...>``. The first 3
+parameters are parsed by Fix base class constructor, while ``<fix
+arguments>`` should be parsed by you.  In our case, we need to specify
+how often we want to print an average velocity. For instance, once in 50
+timesteps: ``fix 1 print/vel 50``. There is a special variable in the
+Fix class called ``nevery`` which specifies how often the method
+``end_of_step()`` is called.  Thus all we need to do is just set it up.
+
+The next method we need to implement is ``setmask()``:
+
+.. code-block:: c++
+
+   int FixPrintVel::setmask()
+   {
+     int mask = 0;
+     mask |= FixConst::END_OF_STEP;
+     return mask;
+   }
+
+Here the we specify which methods of the fix should be called during
+:doc:`execution of a timestep <Developer_flow>`.  The constant
+``END_OF_STEP`` corresponds to the ``end_of_step()`` method. The most
+important available methods that are called during a timestep.
+
+.. code-block:: c++
+
+   void FixPrintVel::end_of_step()
+   {
+     // for add3, scale3
+     using namespace MathExtra;
+
+     double** v = atom->v;
+     int nlocal = atom->nlocal;
+     double localAvgVel[4]; // 4th element for particles count
+     memset(localAvgVel, 0, 4 * sizeof(double));
+     for (int particleInd = 0; particleInd < nlocal; ++particleInd) {
+       add3(localAvgVel, v[particleInd], localAvgVel);
+     }
+     localAvgVel[3] = nlocal;
+     double globalAvgVel[4];
+     memset(globalAvgVel, 0, 4 * sizeof(double));
+     MPI_Allreduce(localAvgVel, globalAvgVel, 4, MPI_DOUBLE, MPI_SUM, world);
+     scale3(1.0 / globalAvgVel[3], globalAvgVel);
+     if ((comm->me == 0) && screen) {
+       fmt::print(screen,"{}, {}, {}\n",
+                  globalAvgVel[0], globalAvgVel[1], globalAvgVel[2]);
+     }
+   }
+
+In the code above, we use MathExtra routines defined in
+``math_extra.h``.  There are bunch of math functions to work with
+arrays of doubles as with math vectors.  It is also important to note
+that LAMMPS code should always assume to be run in parallel and that
+atom data is thus distributed across the MPI ranks.  Thus you can
+only process data from local atoms directly and need to use MPI library
+calls to combine or exchange data.  For serial execution, LAMMPS
+comes bundled with the MPI STUBS library that contains the MPI library
+function calls in dummy versions that only work for a single MPI rank.
+
+In this code we use an instance of Atom class. This object is stored
+in the Pointers class (see ``pointers.h``) which is the base class of
+the Fix base class. This object contains references to various class
+instances (the original instances are created and held by the LAMMPS
+class) with all global information about the simulation system.
+Data from the Pointers class is available to all classes inherited from
+it using protected inheritance. Hence when you write you own class,
+which is going to use LAMMPS data, don't forget to inherit from Pointers
+or pass a Pointer to it to all functions that need access. When writing
+fixes we inherit from class Fix which is inherited from Pointers so
+there is no need to inherit from it directly.
+
+The code above computes average velocity for all particles in the
+simulation.  Yet you have one unused parameter in fix call from the
+script: ``group_name``.  This parameter specifies the group of atoms
+used in the fix. So we should compute average for all particles in the
+simulation only if ``group_name == "all"``, but it can be any group.
+The group membership information of an atom is contained in the *mask*
+property of an atom and the bit corresponding to a given group is
+stored in the groupbit variable which is defined in Fix base class:
+
+.. code-block:: c++
+
+   for (int i = 0; i < nlocal; ++i) {
+     if (atom->mask[i] & groupbit) {
+     // Do all job here
+     }
+   }
+
+Class Atom encapsulates atoms positions, velocities, forces, etc. Users
+can access them using particle index. Note, that particle indexes are
+usually changed every few timesteps because of neighbor list rebuilds
+and spatial sorting (to improve cache efficiency).
+
+Let us consider another Fix example: We want to have a fix which stores
+atoms position from the previous time step in your fix. The local atoms
+indexes may not be valid on the next iteration. In order to handle
+this situation there are several methods which should be implemented:
+
+- ``double memory_usage()``: return how much memory the fix uses (optional)
+- ``void grow_arrays(int)``: do reallocation of the per-particle arrays in your fix
+- ``void copy_arrays(int i, int j, int delflag)``: copy i-th per-particle
+  information to j-th. Used when atom sorting is performed. if delflag is set
+  and atom j owns a body, move the body information to atom i.
+- ``void set_arrays(int i)``: sets i-th particle related information to zero
+
+Note, that if your class implements these methods, it must add calls of
+add_callback and delete_callback to the constructor and destructor. Since we want
+to store positions of atoms from the previous timestep, we need to add
+``double** xold`` to the header file. Than add allocation code
+to the constructor:
+
+.. code-block:: c++
+
+   FixSavePos::FixSavePos(LAMMPS *lmp, int narg, char **arg), xold(nullptr)
+   {
+   //...
+     memory->create(xold, atom->nmax, 3, "FixSavePos:x");
+     atom->add_callback(0);
+   }
+
+   FixSavePos::~FixSavePos() {
+     atom->delete_callback(id, 0);
+     memory->destroy(xold);
+   }
+
+Implement the aforementioned methods:
+
+.. code-block:: c++
+
+   double FixSavePos::memory_usage()
+   {
+     int nmax = atom->nmax;
+     double bytes = 0.0;
+     bytes += nmax * 3 * sizeof(double);
+     return bytes;
+   }
+
+   void FixSavePos::grow_arrays(int nmax)
+   {
+     memory->grow(xold, nmax, 3, "FixSavePos:xold");
+   }
+
+   void FixSavePos::copy_arrays(int i, int j, int delflag)
+   {
+     memcpy(xold[j], xold[i], sizeof(double) * 3);
+   }
+
+   void FixSavePos::set_arrays(int i)
+   {
+     memset(xold[i], 0, sizeof(double) * 3);
+   }
+
+   int FixSavePos::pack_exchange(int i, double *buf)
+   {
+     int m = 0;
+     buf[m++] = xold[i][0];
+     buf[m++] = xold[i][1];
+     buf[m++] = xold[i][2];
+
+     return m;
+   }
+
+   int FixSavePos::unpack_exchange(int nlocal, double *buf)
+   {
+     int m = 0;
+     xold[nlocal][0] = buf[m++];
+     xold[nlocal][1] = buf[m++];
+     xold[nlocal][2] = buf[m++];
+
+     return m;
+   }
+
+Now, a little bit about memory allocation. We use the Memory class which
+is just a bunch of template functions for allocating 1D and 2D
+arrays. So you need to add include ``memory.h`` to have access to them.
+
+Finally, if you need to write/read some global information used in
+your fix to the restart file, you might do it by setting flag
+``restart_global = 1`` in the constructor and implementing methods void
+``write_restart(FILE *fp)`` and ``void restart(char *buf)``.
+If, in addition, you want to write the per-atom property to restart
+files additional settings and functions are needed:
+
+- a fix flag indicating this needs to be set ``restart_peratom = 1;``
+- ``atom->add_callback()`` and ``atom->delete_callback()`` must be called
+  a second time with the final argument set to 1 instead of 0 (indicating
+  restart processing instead of per-atom data memory management).
+- the functions ``void pack_restart(int i, double *buf)`` and
+  ``void unpack_restart(int nlocal, int nth)`` need to be implemented
+

doc/src/Howto.rst

@@ -70,6 +70,7 @@ Force fields howto
    Howto_amoeba
    Howto_tip3p
    Howto_tip4p
+   Howto_tip5p
    Howto_spc
 
 Packages howto

doc/src/Howto_granular.rst

@@ -43,6 +43,15 @@ The fix style *freeze* zeroes both the force and torque of frozen
 atoms, and should be used for granular system instead of the fix style
 *setforce*\ .
 
+To model heat conduction, one must add the temperature and heatflow
+atom variables with:
+* :doc:`fix property/atom <fix_property_atom>`
+a temperature integration fix
+* :doc:`fix heat/flow <fix_heat_flow>`
+and a heat conduction option defined in both
+* :doc:`pair_style granular <pair_granular>`
+* :doc:`fix wall/gran <fix_wall_gran>`
+
 For computational efficiency, you can eliminate needless pairwise
 computations between frozen atoms by using this command:
 
@@ -55,3 +64,6 @@ computations between frozen atoms by using this command:
    will be the same as in 3d.  If you wish to model granular particles in
    2d as 2d discs, see the note on this topic on the :doc:`Howto 2d <Howto_2d>`
    doc page, where 2d simulations are discussed.
+
+To add custom granular contact models, see the
+:doc:`modifying granular sub-models page <Modify_gran_sub_mod>`.

doc/src/Howto_mdi.rst

@@ -63,22 +63,29 @@ The package also provides a :doc:`mdi plugin <mdi>` command, which
 enables LAMMPS to operate as an MDI driver and load an MDI engine as a
 plugin library.
 
-The package furthermore includes a `fix mdi/qm <fix_mdi_qm>` command, in
+The package furthermore includes a :doc:`fix mdi/qm <fix_mdi_qm>`
-which LAMMPS operates as an MDI driver in conjunction with a quantum
+command, in which LAMMPS operates as an MDI driver in conjunction with a
-mechanics code as an MDI engine.  The post_force() method of the
+quantum mechanics code as an MDI engine.  The post_force() method of the
-``fix_mdi_qm.cpp`` file shows how a driver issues MDI commands to another
+``fix_mdi_qm.cpp`` file shows how a driver issues MDI commands to
-code.  This command can be used to couple to an MDI engine, which is
+another code.  This command can be used to couple to an MDI engine,
-either a stand-alone code or a plugin library.
+which is either a stand-alone code or a plugin library.
 
-As explained in the `fix mdi/qm <fix_mdi_qm>` command documentation, it
+As explained in the :doc:`fix mdi/qm <fix_mdi_qm>` command
-can be used to perform *ab initio* MD simulations or energy
+documentation, it can be used to perform *ab initio* MD simulations or
-minimizations, or to evaluate the quantum energy and forces for a series
+energy minimizations, or to evaluate the quantum energy and forces for a
-of independent systems.  The ``examples/mdi`` directory has example
+series of independent systems.  The ``examples/mdi`` directory has
-input scripts for all of these use cases.
+example input scripts for all of these use cases.
+
+The package also has a :doc:`fix mdi/qmmm <fix_mdi_qmmm>` command in
+which LAMMPS operates as an MDI driver in conjunction with a quantum
+mechanics code as an MDI engine to perform QM/MM simulations.  The
+LAMMPS input script partitions the system into QM and MM (molecular
+mechanics) atoms.  As described below the ``examples/QUANTUM`` directory
+has examples for coupling to 3 different quantum codes in this manner.
 
 ----------
 
-The examples/mdi directory contains Python scripts and LAMMPS input
+The ``examples/mdi`` directory contains Python scripts and LAMMPS input
 script which use LAMMPS as either an MDI driver or engine, or both.
 Currently, 5 example use cases are provided:
 
@@ -119,45 +126,26 @@ as a plugin library.
 
 -------------
 
-Currently, there are at least two quantum DFT codes which have direct MDI
+As of March 2023, these are quantum codes with MDI support provided via
+Python wrapper scripts included in the LAMMPS distribution.  These can
+be used with the fix mdi/qm and fix mdi/qmmm commands to perform QM
+calculations of an entire system (e.g. AIMD) or QM/MM simulations.  See
+the ``examples/QUANTUM`` sub-directories for more details:
+
+* LATTE - AIMD only
+* PySCF - QM/MM only
+* NWChem - AIMD or QM/MM
+
+There are also at least two quantum codes which have direct MDI
 support, `Quantum ESPRESSO (QE) <https://www.quantum-espresso.org/>`_
-and `INQ <https://qsg.llnl.gov/node/101.html>`_.  There are also several
+and `INQ <https://qsg.llnl.gov/node/101.html>`_.  There are also
-QM codes which have indirect support through QCEngine or i-PI.  The
+several QM codes which have indirect support through QCEngine or i-PI.
-former means they require a wrapper program (QCEngine) with MDI support
+The former means they require a wrapper program (QCEngine) with MDI
-which writes/read files to pass data to the quantum code itself.  The
+support which writes/read files to pass data to the quantum code
-list of QCEngine-supported and i-PI-supported quantum codes is on the
+itself.  The list of QCEngine-supported and i-PI-supported quantum
-`MDI webpage
+codes is on the `MDI webpage
 <https://molssi-mdi.github.io/MDI_Library/html/index.html>`_.
 
-Here is how to build QE as a stand-alone ``pw.x`` file which can be
+These direct- and indirect-support codes should be usable for full
-used in stand-alone mode:
+system calculations (e.g. AIMD).  Whether they support QM/MM models
-
+depends on the individual QM code.
-.. code-block:: bash
-
-   git clone --branch mdi_plugin https://github.com/MolSSI-MDI/q-e.git <base_path>/q-e
-   build the executable pw.x, following the `QE build guide <https://gitlab.com/QEF/q-e/-/wikis/Developers/CMake-build-system>`_
-
-Here is how to build QE as a shared library which can be used in plugin mode,
-which results in a ``libqemdi.so`` file in ``<base_path>/q-e/MDI/src``:
-
-.. code-block:: bash
-
-   git clone --branch mdi_plugin https://github.com/MolSSI-MDI/q-e.git <base_path>/q-e
-   cd <base_path>/q-e
-   ./configure --enable-parallel --enable-openmp --enable-shared FFLAGS="-fPIC" FCFLAGS="-fPIC" CFLAGS="-fPIC" foxflags="-fPIC" try_foxflags="-fPIC"
-   make -j 4 mdi
-
-INQ cannot be built as a stand-alone code; it is by design a library.
-Here is how to build INQ as a shared library which can be used in
-plugin mode, which results in a ``libinqmdi.so`` file in
-``<base_path>/inq/build/examples``:
-
-.. code-block:: bash
-
-   git clone --branch mdi --recurse-submodules https://gitlab.com/taylor-a-barnes/inq.git <base_path>/inq
-   cd <base_path>/inq
-   mkdir -p build
-   cd build
-   ../configure --prefix=<install_path>/install
-   make -j 4
-   make install

doc/src/Howto_spc.rst

@@ -20,7 +20,6 @@ atoms and the water molecule to run a rigid SPC model.
 | LJ :math:`\epsilon`, :math:`\sigma` of OH, HH = 0.0
 | :math:`r_0` of OH bond = 1.0
 | :math:`\theta_0` of HOH angle = 109.47\ :math:`^{\circ}`
-|
 
 Note that as originally proposed, the SPC model was run with a 9
 Angstrom cutoff for both LJ and Coulomb terms.  It can also be used
@@ -33,16 +32,123 @@ the partial charge assignments change:
 
 | O charge = -0.8476
 | H charge = 0.4238
-|
 
 See the :ref:`(Berendsen) <howto-Berendsen>` reference for more details on both
 the SPC and SPC/E models.
 
+Below is the code for a LAMMPS input file and a molecule file
+(``spce.mol``) of SPC/E water for use with the :doc:`molecule command
+<molecule>` demonstrating how to set up a small bulk water system for
+SPC/E with rigid bonds.
+
+.. code-block:: LAMMPS
+
+   units real
+   atom_style full
+   region box block -5 5 -5 5 -5 5
+   create_box 2 box  bond/types 1 angle/types 1 &
+                   extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2
+
+   mass 1 15.9994
+   mass 2 1.008
+
+   pair_style lj/cut/coul/cut 10.0
+   pair_coeff 1 1 0.1553 3.166
+   pair_coeff 1 2 0.0    1.0
+   pair_coeff 2 2 0.0    1.0
+
+   bond_style zero
+   bond_coeff 1 1.0
+
+   angle_style zero
+   angle_coeff 1 109.47
+
+   molecule water spce.mol
+   create_atoms 0 random 33 34564 NULL mol water 25367 overlap 1.33
+
+   timestep 1.0
+   fix rigid     all shake 0.0001 10 10000 b 1 a 1
+   minimize 0.0 0.0 1000 10000
+   run 0 post no
+   reset_timestep 0
+   velocity all create 300.0 5463576
+   fix integrate all nvt temp 300.0 300.0 1.0
+
+   thermo_style custom step temp press etotal density pe ke
+   thermo 1000
+   run 20000 upto
+   write_data tip4p.data nocoeff
+
+.. _spce_molecule:
+.. code-block::
+
+   # Water molecule. SPC/E geometry
+
+   3 atoms
+   2 bonds
+   1 angles
+
+   Coords
+
+   1    0.00000  -0.06461   0.00000
+   2    0.81649   0.51275   0.00000
+   3   -0.81649   0.51275   0.00000
+
+   Types
+
+   1        1   # O
+   2        2   # H
+   3        2   # H
+
+   Charges
+
+   1       -0.8476
+   2        0.4238
+   3        0.4238
+
+   Bonds
+
+   1   1      1      2
+   2   1      1      3
+
+   Angles
+
+   1   1      2      1      3
+
+   Shake Flags
+
+   1 1
+   2 1
+   3 1
+
+   Shake Atoms
+
+   1 1 2 3
+   2 1 2 3
+   3 1 2 3
+
+   Shake Bond Types
+
+   1 1 1 1
+   2 1 1 1
+   3 1 1 1
+
+   Special Bond Counts
+
+   1 2 0 0
+   2 1 1 0
+   3 1 1 0
+
+   Special Bonds
+
+   1 2 3
+   2 1 3
+   3 1 2
+
 Wikipedia also has a nice article on `water models <https://en.wikipedia.org/wiki/Water_model>`_.
 
 ----------
 
 .. _howto-Berendsen:
 
-**(Berendsen)** Berendsen, Grigera, Straatsma, J Phys Chem, 91,
+**(Berendsen)** Berendsen, Grigera, Straatsma, J Phys Chem, 91, 6269-6271 (1987).
-6269-6271 (1987).

doc/src/Howto_tip3p.rst

@@ -1,53 +1,211 @@
 TIP3P water model
 =================
 
-The TIP3P water model as implemented in CHARMM
+The TIP3P water model as implemented in CHARMM :ref:`(MacKerell)
-:ref:`(MacKerell) <howto-tip3p>` specifies a 3-site rigid water molecule with
+<howto-tip3p>` specifies a 3-site rigid water molecule with charges and
-charges and Lennard-Jones parameters assigned to each of the 3 atoms.
+Lennard-Jones parameters assigned to each of the 3 atoms.
-In LAMMPS the :doc:`fix shake <fix_shake>` command can be used to hold
-the two O-H bonds and the H-O-H angle rigid.  A bond style of
-*harmonic* and an angle style of *harmonic* or *charmm* should also be
-used.
 
-These are the additional parameters (in real units) to set for O and H
+A suitable pair style with cutoff Coulomb would be:
-atoms and the water molecule to run a rigid TIP3P-CHARMM model with a
-cutoff.  The K values can be used if a flexible TIP3P model (without
-fix shake) is desired.  If the LJ epsilon and sigma for HH and OH are
-set to 0.0, it corresponds to the original 1983 TIP3P model
-:ref:`(Jorgensen) <Jorgensen1>`.
 
-| O mass = 15.9994
+* :doc:`pair_style lj/cut/coul/cut <pair_lj_cut_coul>`
-| H mass = 1.008
-| O charge = -0.834
-| H charge = 0.417
-| LJ :math:`\epsilon` of OO = 0.1521
-| LJ :math:`\sigma` of OO = 3.1507
-| LJ :math:`\epsilon` of HH = 0.0460
-| LJ :math:`\sigma` of HH = 0.4000
-| LJ :math:`\epsilon` of OH = 0.0836
-| LJ :math:`\sigma` of OH = 1.7753
-| K of OH bond = 450
-| :math:`r_0` of OH bond = 0.9572
-| K of HOH angle = 55
-| :math:`\theta` of HOH angle = 104.52\ :math:`^{\circ}`
-|
 
-These are the parameters to use for TIP3P with a long-range Coulomb
+or these commands for a long-range Coulomb model:
-solver (e.g. Ewald or PPPM in LAMMPS), see :ref:`(Price) <Price1>` for
+
-details:
+* :doc:`pair_style lj/cut/coul/long <pair_lj_cut_coul>`
+* :doc:`pair_style lj/cut/coul/long/soft <pair_fep_soft>`
+* :doc:`kspace_style pppm <kspace_style>`
+* :doc:`kspace_style pppm/disp <kspace_style>`
+
+In LAMMPS the :doc:`fix shake or fix rattle <fix_shake>` command can be
+used to hold the two O-H bonds and the H-O-H angle rigid.  A bond style
+of :doc:`harmonic <bond_harmonic>` and an angle style of :doc:`harmonic
+<angle_harmonic>` or :doc:`charmm <angle_charmm>` should also be used.
+In case of rigid bonds also bond style :doc:`zero <bond_zero>` and angle
+style :doc:`zero <angle_zero>` can be used.
+
+The table below lists the force field parameters (in real :doc:`units
+<units>`) to for the water molecule atoms to run a rigid or flexible
+TIP3P-CHARMM model with a cutoff, the original 1983 TIP3P model
+:ref:`(Jorgensen) <Jorgensen1>`, or a TIP3P model with parameters
+optimized for a long-range Coulomb solver (e.g. Ewald or PPPM in LAMMPS)
+:ref:`(Price) <Price1>`.   The K values can be used if a flexible TIP3P
+model (without fix shake) is desired, for rigid bonds/angles they are
+ignored.
+
+   .. list-table::
+      :header-rows: 1
+      :widths: auto
+
+      * - Parameter
+        - TIP3P-CHARMM
+        - TIP3P (original)
+        - TIP3P (Ewald)
+      * - O mass (amu)
+        - 15.9994
+        - 15.9994
+        - 15.9994
+      * - H mass (amu)
+        - 1.008
+        - 1.008
+        - 1.008
+      * - O charge (:math:`e`)
+        - -0.834
+        - -0.834
+        - -0.834
+      * - H charge (:math:`e`)
+        - 0.417
+        - 0.417
+        - 0.417
+      * - LJ :math:`\epsilon` of OO (kcal/mole)
+        - 0.1521
+        - 0.1521
+        - 0.1020
+      * - LJ :math:`\sigma` of OO (:math:`\AA`)
+        - 3.1507
+        - 3.1507
+        - 3.188
+      * - LJ :math:`\epsilon` of HH (kcal/mole)
+        - 0.0460
+        - 0.0
+        - 0.0
+      * - LJ :math:`\sigma` of HH (:math:`\AA`)
+        - 0.4
+        - 1.0
+        - 1.0
+      * - LJ :math:`\epsilon` of OH (kcal/mole)
+        - 0.0836
+        - 0.0
+        - 0.0
+      * - LJ :math:`\sigma` of OH (:math:`\AA`)
+        - 1.7753
+        - 1.0
+        - 1.0
+      * - K of OH bond (kcal/mole/:math:`\AA^2`)
+        - 450
+        - 450
+        - 450
+      * - :math:`r_0` of OH bond (:math:`\AA`)
+        - 0.9572
+        - 0.9572
+        - 0.9572
+      * - K of HOH angle (kcal/mole)
+        - 55.0
+        - 55.0
+        - 55.0
+      * - :math:`\theta_0` of HOH angle
+        - 104.52\ :math:`^{\circ}`
+        - 104.52\ :math:`^{\circ}`
+        - 104.52\ :math:`^{\circ}`
+
+Below is the code for a LAMMPS input file and a molecule file
+(``tip3p.mol``) of TIP3P water for use with the :doc:`molecule command
+<molecule>` demonstrating how to set up a small bulk water system for
+TIP3P with rigid bonds.
+
+.. code-block:: LAMMPS
+
+    units real
+    atom_style full
+    region box block -5 5 -5 5 -5 5
+    create_box 2 box bond/types 1 angle/types 1 &
+                extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2
+
+    mass 1 15.9994
+    mass 2 1.008
+
+    pair_style lj/cut/coul/cut 8.0
+    pair_coeff 1 1 0.1521 3.1507
+    pair_coeff 2 2 0.0    1.0
+
+    bond_style zero
+    bond_coeff 1 0.9574
+
+    angle_style zero
+    angle_coeff 1 104.52
+
+    molecule water tip3p.mol
+    create_atoms 0 random 33 34564 NULL mol water 25367 overlap 1.33
+
+    fix rigid all shake 0.001 10 10000 b 1 a 1
+    minimize 0.0 0.0 1000 10000
+    run 0 post no
+
+    reset_timestep 0
+    velocity all create 300.0 5463576
+    fix integrate all nvt temp 300 300 1.0
+
+    thermo_style custom step temp press etotal pe
+
+    thermo 1000
+    run 20000
+    write_data tip3p.data nocoeff
+
+.. _tip3p_molecule:
+.. code-block::
+
+   # Water molecule. TIP3P geometry
+
+   3 atoms
+   2 bonds
+   1 angles
+
+   Coords
+
+   1    0.00000  -0.06556   0.00000
+   2    0.75695   0.52032   0.00000
+   3   -0.75695   0.52032   0.00000
+
+   Types
+
+   1        1   # O
+   2        2   # H
+   3        2   # H
+
+   Charges
+
+   1       -0.834
+   2        0.417
+   3        0.417
+
+   Bonds
+
+   1   1      1      2
+   2   1      1      3
+
+   Angles
+
+   1   1      2      1      3
+
+   Shake Flags
+
+   1 1
+   2 1
+   3 1
+
+   Shake Atoms
+
+   1 1 2 3
+   2 1 2 3
+   3 1 2 3
+
+   Shake Bond Types
+
+   1 1 1 1
+   2 1 1 1
+   3 1 1 1
+
+   Special Bond Counts
+
+   1 2 0 0
+   2 1 1 0
+   3 1 1 0
+
+   Special Bonds
+
+   1 2 3
+   2 1 3
+   3 1 2
 
-| O mass = 15.9994
-| H mass = 1.008
-| O charge = -0.830
-| H charge = 0.415
-| LJ :math:`\epsilon` of OO = 0.102
-| LJ :math:`\sigma` of OO = 3.188
-| LJ :math:`\epsilon`, :math:`\sigma` of OH, HH = 0.0
-| K of OH bond = 450
-| :math:`r_0` of OH bond = 0.9572
-| K of HOH angle = 55
-| :math:`\theta` of HOH angle = 104.52\ :math:`^{\circ}`
-|
 
 Wikipedia also has a nice article on `water models <https://en.wikipedia.org/wiki/Water_model>`_.
 

doc/src/Howto_tip4p.rst

@@ -2,100 +2,138 @@ TIP4P water model
 =================
 
 The four-point TIP4P rigid water model extends the traditional
-three-point TIP3P model by adding an additional site, usually
+:doc:`three-point TIP3P <Howto_tip3p>` model by adding an additional
-massless, where the charge associated with the oxygen atom is placed.
+site M, usually massless, where the charge associated with the oxygen
-This site M is located at a fixed distance away from the oxygen along
+atom is placed.  This site M is located at a fixed distance away from
-the bisector of the HOH bond angle.  A bond style of *harmonic* and an
+the oxygen along the bisector of the HOH bond angle.  A bond style of
-angle style of *harmonic* or *charmm* should also be used.
+:doc:`harmonic <bond_harmonic>` and an angle style of :doc:`harmonic
+<angle_harmonic>` or :doc:`charmm <angle_charmm>` should also be used.
+In case of rigid bonds also bond style :doc:`zero <bond_zero>` and angle
+style :doc:`zero <angle_zero>` can be used.
 
-A TIP4P model is run with LAMMPS using either these commands
+There are two ways to implement TIP4P water in LAMMPS:
-for a cutoff model:
 
-* :doc:`pair_style tip4p/cut <pair_lj_cut_tip4p>`
+#. Use a specially written pair style that uses the :ref:`TIP3P geometry
-* :doc:`pair_style lj/cut/tip4p/cut <pair_lj_cut_tip4p>`
+   <tip3p_molecule>` without the point M. The point M location is then
+   implicitly derived from the other atoms or each water molecule and
+   used during the force computation.  The forces on M are then
+   projected on the oxygen and the two hydrogen atoms.  This is
+   computationally very efficient, but the charge distribution in space
+   is only correct within the tip4p labeled styles.  So all other
+   computations using charges will "see" the negative charge incorrectly
+   on the oxygen atom.
 
-or these commands for a long-range model:
+   This can be done with the following pair styles for Coulomb with a cutoff:
 
-* :doc:`pair_style tip4p/long <pair_coul>`
+   * :doc:`pair_style tip4p/cut <pair_lj_cut_tip4p>`
-* :doc:`pair_style lj/cut/tip4p/long <pair_lj_cut_tip4p>`
+   * :doc:`pair_style lj/cut/tip4p/cut <pair_lj_cut_tip4p>`
-* :doc:`pair_style lj/long/tip4p/long <pair_lj_long>`
-* :doc:`pair_style tip4p/long/soft <pair_fep_soft>`
-* :doc:`pair_style lj/cut/tip4p/long/soft <pair_fep_soft>`
-* :doc:`kspace_style pppm/tip4p <kspace_style>`
-* :doc:`kspace_style pppm/disp/tip4p <kspace_style>`
 
-The bond lengths and bond angles should be held fixed using the
+   or these commands for a long-range Coulomb treatment:
-:doc:`fix shake <fix_shake>` or :doc:`fix rattle <fix_shake>` command,
-unless a parameterization for a flexible TIP4P model is used.  The
-parameter sets listed below are all for rigid TIP4P model variants and
-thus the bond and angle force constants are not used and can be set to
-any legal value; only equilibrium length and angle are used.
 
-These are the additional parameters (in real units) to set for O and H
+   * :doc:`pair_style tip4p/long <pair_coul>`
-atoms and the water molecule to run a rigid TIP4P model with a cutoff
+   * :doc:`pair_style lj/cut/tip4p/long <pair_lj_cut_tip4p>`
-:ref:`(Jorgensen) <Jorgensen5>`.  Note that the OM distance is specified in
+   * :doc:`pair_style lj/long/tip4p/long <pair_lj_long>`
-the :doc:`pair_style <pair_style>` command, not as part of the pair
+   * :doc:`pair_style tip4p/long/soft <pair_fep_soft>`
-coefficients.
+   * :doc:`pair_style lj/cut/tip4p/long/soft <pair_fep_soft>`
+   * :doc:`kspace_style pppm/tip4p <kspace_style>`
+   * :doc:`kspace_style pppm/disp/tip4p <kspace_style>`
 
-| O mass = 15.9994
+   The bond lengths and bond angles should be held fixed using the
-| H mass = 1.008
+   :doc:`fix shake <fix_shake>` or :doc:`fix rattle <fix_shake>` command,
-| O charge = -1.040
+   unless a parameterization for a flexible TIP4P model is used.  The
-| H charge = 0.520
+   parameter sets listed below are all for rigid TIP4P model variants and
-| :math:`r_0` of OH bond = 0.9572
+   thus the bond and angle force constants are not used and can be set to
-| :math:`\theta` of HOH angle = 104.52\ :math:`^{\circ}`
+   any legal value; only equilibrium length and angle are used.
-| OM distance = 0.15
-| LJ :math:`\epsilon` of O-O = 0.1550
-| LJ :math:`\sigma` of O-O = 3.1536
-| LJ :math:`\epsilon`, :math:`\sigma` of OH, HH = 0.0
-| Coulomb cutoff = 8.5
-|
 
-For the TIP4/Ice model (J Chem Phys, 122, 234511 (2005);
+#. Use an :ref:`explicit 4 point TIP4P geometry <tip4p_molecule>` where
-https://doi.org/10.1063/1.1931662) these values can be used:
+   the oxygen atom carries no charge and the M point no Lennard-Jones
+   interactions.  Since :doc:`fix shake <fix_shake>` or :doc:`fix rattle
+   <fix_shake>` may not be applied to this kind of geometry, :doc:`fix
+   rigid or fix rigid/small <fix_rigid>` or its thermostatted variants
+   are required to maintain a rigid geometry.  This avoids some of the
+   issues with respect to analysis and non-tip4p styles, but it is a
+   more costly force computation (more atoms in the same volume and thus
+   more neighbors in the neighbor lists) and requires a much shorter
+   timestep for stable integration of the rigid body motion.  Since no
+   bonds or angles are required, they do not need to be defined and atom
+   style charge would be sufficient for a bulk TIP4P water system.  In
+   order to avoid that LAMMPS produces an error due to the massless M
+   site a tiny non-zero mass needs to be assigned.
 
-| O mass = 15.9994
+The table below lists the force field parameters (in real :doc:`units
-| H mass =  1.008
+<units>`) to for a selection of popular variants of the TIP4P model.
-| O charge = -1.1794
+There is the rigid TIP4P model with a cutoff :ref:`(Jorgensen)
-| H charge =  0.5897
+<Jorgensen5>`, the TIP4/Ice model :ref:`(Abascal1) <Abascal1>`, the
-| :math:`r_0` of OH bond = 0.9572
+TIP4P/2005 model :ref:`(Abascal2) <Abascal2>` and a version of TIP4P
-| :math:`\theta` of HOH angle = 104.52\ :math:`^{\circ}`
+parameters adjusted for use with a long-range Coulombic solver
-| OM distance = 0.1577
+(e.g. Ewald or PPPM in LAMMPS).  Note that for implicit TIP4P models the
-| LJ :math:`\epsilon` of O-O = 0.21084
+OM distance is specified in the :doc:`pair_style <pair_style>` command,
-| LJ :math:`\sigma` of O-O = 3.1668
+not as part of the pair coefficients.
-| LJ :math:`\epsilon`, :math:`\sigma` of OH, HH = 0.0
-| Coulomb cutoff = 8.5
-|
 
-For the TIP4P/2005 model (J Chem Phys, 123, 234505 (2005);
+   .. list-table::
-https://doi.org/10.1063/1.2121687), these values can be used:
+      :header-rows: 1
+      :widths: auto
 
-| O mass = 15.9994
+      * - Parameter
-| H mass =  1.008
+        - TIP4P (original)
-| O charge = -1.1128
+        - TIP4P/Ice
-| H charge = 0.5564
+        - TIP4P/2005
-| :math:`r_0` of OH bond = 0.9572
+        - TIP4P (Ewald)
-| :math:`\theta` of HOH angle = 104.52\ :math:`^{\circ}`
+      * - O mass (amu)
-| OM distance = 0.1546
+        - 15.9994
-| LJ :math:`\epsilon` of O-O = 0.1852
+        - 15.9994
-| LJ :math:`\sigma` of O-O = 3.1589
+        - 15.9994
-| LJ :math:`\epsilon`, :math:`\sigma` of OH, HH = 0.0
+        - 15.9994
-| Coulomb cutoff = 8.5
+      * - H mass (amu)
-|
+        - 1.008
-
+        - 1.008
-These are the parameters to use for TIP4P with a long-range Coulombic
+        - 1.008
-solver (e.g. Ewald or PPPM in LAMMPS):
+        - 1.008
-
+      * - O or M charge (:math:`e`)
-| O mass = 15.9994
+        - -1.040
-| H mass = 1.008
+        - -1.1794
-| O charge = -1.0484
+        - -1.1128
-| H charge = 0.5242
+        - -1.04844
-| :math:`r_0` of OH bond = 0.9572
+      * - H charge (:math:`e`)
-| :math:`\theta` of HOH angle = 104.52\ :math:`^{\circ}`
+        - 0.520
-| OM distance = 0.1250
+        - 0.5897
-| LJ :math:`\epsilon` of O-O = 0.16275
+        - 0.5564
-| LJ :math:`\sigma` of O-O = 3.16435
+        - 0.52422
-| LJ :math:`\epsilon`, :math:`\sigma` of OH, HH = 0.0
+      * - LJ :math:`\epsilon` of OO (kcal/mole)
-|
+        - 0.1550
+        - 0.1577
+        - 0.1852
+        - 0.16275
+      * - LJ :math:`\sigma` of OO (:math:`\AA`)
+        - 3.1536
+        - 3.1668
+        - 3.1589
+        - 3.16435
+      * - LJ :math:`\epsilon` of HH, MM, OH, OM, HM (kcal/mole)
+        - 0.0
+        - 0.0
+        - 0.0
+        - 0.0
+      * - LJ :math:`\sigma` of HH, MM, OH, OM, HM (:math:`\AA`)
+        - 1.0
+        - 1.0
+        - 1.0
+        - 1.0
+      * - :math:`r_0` of OH bond (:math:`\AA`)
+        - 0.9572
+        - 0.9572
+        - 0.9572
+        - 0.9572
+      * - :math:`\theta_0` of HOH angle
+        - 104.52\ :math:`^{\circ}`
+        - 104.52\ :math:`^{\circ}`
+        - 104.52\ :math:`^{\circ}`
+        - 104.52\ :math:`^{\circ}`
+      * - OM distance (:math:`\AA`)
+        - 0.15
+        - 0.1577
+        - 0.1546
+        - 0.1250
 
 Note that the when using the TIP4P pair style, the neighbor list cutoff
 for Coulomb interactions is effectively extended by a distance 2 \* (OM
@@ -108,6 +146,117 @@ trade-off for your model.  The OM distance and the LJ and Coulombic
 cutoffs are set in the :doc:`pair_style lj/cut/tip4p/long
 <pair_lj_cut_tip4p>` command.
 
+Below is the code for a LAMMPS input file using the implicit method and
+the :ref:`TIP3P molecule file <tip3p_molecule>`.  Because the TIP4P
+charges are different from TIP3P they need to be reset (or the molecule
+file changed):
+
+.. code-block:: LAMMPS
+
+    units real
+    atom_style full
+    region box block -5 5 -5 5 -5 5
+    create_box 2 box bond/types 1 angle/types 1 &
+                extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2
+
+    mass 1 15.9994
+    mass 2 1.008
+
+    pair_style lj/cut/tip4p/cut 1 2 1 1 0.15 8.0
+    pair_coeff 1 1 0.1550 3.1536
+    pair_coeff 2 2 0.0    1.0
+
+    bond_style zero
+    bond_coeff 1 0.9574
+
+    angle_style zero
+    angle_coeff 1 104.52
+
+    molecule water tip3p.mol  # this uses the TIP3P geometry
+    create_atoms 0 random 33 34564 NULL mol water 25367 overlap 1.33
+    # must change charges for TIP4P
+    set type 1 charge -1.040
+    set type 2 charge  0.520
+
+    fix rigid all shake 0.001 10 10000 b 1 a 1
+    minimize 0.0 0.0 1000 10000
+    run 0 post no
+
+    reset_timestep 0
+    velocity all create 300.0 5463576
+    fix integrate all nvt temp 300 300 1.0
+
+    thermo_style custom step temp press etotal pe
+
+    thermo 1000
+    run 20000
+    write_data tip3p.data nocoeff
+
+Below is the code for a LAMMPS input file using the explicit method and
+a TIP4P molecule file.  Because of using :doc:`fix rigid/nvt/small
+<fix_rigid>` no bonds need to be defined and thus no extra storage needs
+to be reserved for them, but we need to switch to atom style full or use
+:doc:`fix property/atom mol <fix_property_atom>` so that fix
+rigid/nvt/small can identify rigid bodies by their molecule ID:
+
+.. code-block:: LAMMPS
+
+    units real
+    atom_style charge
+    region box block -5 5 -5 5 -5 5
+    create_box 3 box
+
+    mass 1 15.9994
+    mass 2 1.008
+    mass 3 1.0e-100
+
+    pair_style lj/cut/coul/cut 8.0
+    pair_coeff 1 1 0.1550 3.1536
+    pair_coeff 2 2 0.0    1.0
+    pair_coeff 3 3 0.0    1.0
+
+    fix mol all property/atom mol
+    molecule water tip4p.mol
+    create_atoms 0 random 33 34564 NULL mol water 25367 overlap 1.33
+
+    timestep 0.1
+    fix integrate all rigid/nvt/small molecule temp 300.0 300.0 1.0
+    velocity all create 300.0 5463576
+
+    thermo_style custom step temp press etotal density pe ke
+    thermo 1000
+    run 20000
+    write_data tip4p.data nocoeff
+
+.. _tip4p_molecule:
+.. code-block::
+
+   # Water molecule. Explicit TIP4P geometry for use with fix rigid
+
+   4 atoms
+
+   Coords
+
+   1    0.00000  -0.06556   0.00000
+   2    0.75695   0.52032   0.00000
+   3   -0.75695   0.52032   0.00000
+   4    0.00000   0.08444   0.00000
+
+   Types
+
+   1        1   # O
+   2        2   # H
+   3        2   # H
+   4        3   # M
+
+   Charges
+
+   1        0.000
+   2        0.520
+   3        0.520
+   4       -1.040
+
+
 Wikipedia also has a nice article on `water models <https://en.wikipedia.org/wiki/Water_model>`_.
 
 ----------
@@ -116,3 +265,13 @@ Wikipedia also has a nice article on `water models <https://en.wikipedia.org/wik
 
 **(Jorgensen)** Jorgensen, Chandrasekhar, Madura, Impey, Klein, J Chem
 Phys, 79, 926 (1983).
+
+.. _Abascal1:
+
+**(Abascal1)** Abascal, Sanz, Fernandez, Vega, J Chem Phys, 122, 234511 (2005)
+   https://doi.org/10.1063/1.1931662
+
+.. _Abascal2:
+
+**(Abascal2)** Abascal, J Chem Phys, 123, 234505 (2005)
+   https://doi.org/10.1063/1.2121687

doc/src/Howto_tip5p.rst

@@ -0,0 +1,161 @@
+TIP5P water model
+=================
+
+The five-point TIP5P rigid water model extends the :doc:`three-point
+TIP3P model <Howto_tip3p>` by adding two additional sites L, usually
+massless, where the charge associated with the oxygen atom is placed.
+These sites L are located at a fixed distance away from the oxygen atom,
+forming a tetrahedral angle that is rotated by 90 degrees from the HOH
+plane.  Those sites thus somewhat approximate lone pairs of the oxygen
+and consequently improve the water structure to become even more
+"tetrahedral" in comparison to the :doc:`four-point TIP4P model
+<Howto_tip4p>`.
+
+A suitable pair style with cutoff Coulomb would be:
+
+* :doc:`pair_style lj/cut/coul/cut <pair_lj_cut_coul>`
+
+or these commands for a long-range model:
+
+* :doc:`pair_style lj/cut/coul/long <pair_lj_cut_coul>`
+* :doc:`pair_style lj/cut/coul/long/soft <pair_fep_soft>`
+* :doc:`kspace_style pppm <kspace_style>`
+* :doc:`kspace_style pppm/disp <kspace_style>`
+
+A TIP5P model *must* be run using a :doc:`rigid fix <fix_rigid>` since
+there is no other option to keep this kind of structure rigid in LAMMPS.
+In order to avoid that LAMMPS produces an error due to the massless L
+sites, those need to be assigned a tiny non-zero mass.
+
+The table below lists the force field parameters (in real :doc:`units
+<units>`) to for a the TIP5P model with a cutoff :ref:`(Mahoney)
+<Mahoney>` and the TIP5P-E model :ref:`(Rick) <Rick>` for use with a
+long-range Coulombic solver (e.g. Ewald or PPPM in LAMMPS).
+
+   .. list-table::
+      :header-rows: 1
+      :widths: auto
+
+      * - Parameter
+        - TIP5P
+        - TIP5P-E
+      * - O mass (amu)
+        - 15.9994
+        - 15.9994
+      * - H mass (amu)
+        - 1.008
+        - 1.008
+      * - O charge (:math:`e`)
+        - 0.0
+        - 0.0
+      * - L charge (:math:`e`)
+        - -0.241
+        - -0.241
+      * - H charge (:math:`e`)
+        - 0.241
+        - 0.241
+      * - LJ :math:`\epsilon` of OO (kcal/mole)
+        - 0.1600
+        - 0.1780
+      * - LJ :math:`\sigma` of OO (:math:`\AA`)
+        - 3.1200
+        - 3.0970
+      * - LJ :math:`\epsilon` of HH, LL, OH, OL, HL (kcal/mole)
+        - 0.0
+        - 0.0
+      * - LJ :math:`\sigma` of HH, LL, OH, OL, HL (:math:`\AA`)
+        - 1.0
+        - 1.0
+      * - :math:`r_0` of OH bond (:math:`\AA`)
+        - 0.9572
+        - 0.9572
+      * - :math:`\theta_0` of HOH angle
+        - 104.52\ :math:`^{\circ}`
+        - 104.52\ :math:`^{\circ}`
+      * - OL distance (:math:`\AA`)
+        - 0.70
+        - 0.70
+      * - :math:`\theta_0` of LOL angle
+        - 109.47\ :math:`^{\circ}`
+        - 109.47\ :math:`^{\circ}`
+
+Below is the code for a LAMMPS input file for setting up a simulation of
+TIP5P water with a molecule file.  Because of using :doc:`fix
+rigid/nvt/small <fix_rigid>` no bonds need to be defined and thus no
+extra storage needs to be reserved for them, but we need to switch to
+atom style full or use :doc:`fix property/atom mol <fix_property_atom>`
+so that fix rigid/nvt/small can identify rigid bodies by their molecule
+ID:
+
+.. code-block:: LAMMPS
+
+    units real
+    atom_style charge
+    region box block -5 5 -5 5 -5 5
+    create_box 3 box
+
+    mass 1 15.9994
+    mass 2 1.008
+    mass 3 1.0e-100
+
+    pair_style lj/cut/coul/cut 8.0
+    pair_coeff 1 1 0.160  3.12
+    pair_coeff 2 2 0.0    1.0
+    pair_coeff 3 3 0.0    1.0
+
+    fix mol all property/atom mol
+    molecule water tip5p.mol
+    create_atoms 0 random 33 34564 NULL mol water 25367 overlap 1.33
+
+    timestep 0.20
+    fix integrate all rigid/nvt/small molecule temp 300.0 300.0 1.0
+    reset_timestep 0
+    velocity all create 300.0 5463576
+
+    thermo_style custom step temp press etotal density pe ke
+    thermo 1000
+    run 20000
+    write_data tip5p.data nocoeff
+
+.. _tip5p_molecule:
+.. code-block::
+
+   # Water molecule. Explicit TIP5P geometry for use with fix rigid
+
+   5 atoms
+
+   Coords
+
+   1    0.00000  -0.06556   0.00000
+   2    0.75695   0.52032   0.00000
+   3   -0.75695   0.52032   0.00000
+   4    0.00000  -0.46971   0.57154
+   5    0.00000  -0.46971  -0.57154
+
+   Types
+
+   1        1   # O
+   2        2   # H
+   3        2   # H
+   4        3   # L
+   5        3   # L
+
+   Charges
+
+   1        0.000
+   2        0.241
+   3        0.241
+   4       -0.241
+   5       -0.241
+
+Wikipedia also has a nice article on `water models <https://en.wikipedia.org/wiki/Water_model>`_.
+
+----------
+
+.. _Mahoney:
+
+**(Mahoney)** Mahoney, Jorgensen, J Chem Phys 112, 8910 (2000)
+
+.. _Rick:
+
+**(Rick)** Rick, J Chem Phys 120, 6085 (2004)

doc/src/Install_linux.rst

@@ -42,7 +42,7 @@ static linkage, there is no ``liblammps.so`` library file and thus also the
 LAMMPS python module, which depends on it, is not included.
 
 The compressed tar archives available for download have names following
-the pattern `lammps-linux-x86_64-<version>.tar.gz` and will all unpack
+the pattern ``lammps-linux-x86_64-<version>.tar.gz`` and will all unpack
 into a ``lammps-static`` folder.  The executables are then in the
 ``lammps-static/bin/`` folder.  Since they do not depend on any other
 software, they may be freely moved or copied around.
@@ -172,7 +172,7 @@ Pre-built EPEL Linux executable
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 Pre-built LAMMPS (and KIM) packages for stable releases are available
-in the `Extra Packages for Enterprise Linux (EPEL) repository <https://fedoraproject.org/wiki/EPEL>`_
+in the `Extra Packages for Enterprise Linux (EPEL) repository <https://docs.fedoraproject.org/en-US/epel/>`_
 for use with Red Hat Enterprise Linux (RHEL) or CentOS version 7.x
 and compatible Linux distributions. Names of packages, executable,
 and content are the same as described above for Fedora Linux.

doc/src/Install_windows.rst

@@ -17,11 +17,12 @@ install the Windows MPI package (MPICH2 from Argonne National Labs),
 needed to run in parallel with MPI.
 
 The LAMMPS binaries contain *all* :doc:`optional packages <Packages>`
-included in the source distribution except: KIM, KOKKOS, MSCG, PYTHON,
+included in the source distribution except: ADIOS, H5MD, KIM, ML-PACE,
-ADIOS, H5MD, NETCDF, QMMM, ML-QUIP, and VTK.
+ML-QUIP, MSCG, NETCDF, PLUMED, QMMM, SCAFACOS, and VTK.  The serial
-The serial version also does not include the MPIIO and
+version also does not include the MPIIO and LATBOLTZ packages.  The
-LATBOLTZ packages.  The GPU package is compiled for OpenCL with
+PYTHON package is only available in the Python installers that bundle a
-mixed precision kernels.
+Python runtime.  The GPU package is compiled for OpenCL with mixed
+precision kernels.
 
 The LAMMPS library is compiled as a shared library and the
 :doc:`LAMMPS Python module <Python_module>` is installed, so that

doc/src/Intro_features.rst

@@ -195,7 +195,7 @@ Multi-replica models
 * :doc:`parallel replica dynamics <prd>`
 * :doc:`temperature accelerated dynamics <tad>`
 * :doc:`parallel tempering <temper>`
-* path-integral MD: `first variant <fix_pimd>`, `second variant <fix_ipi>`
+* path-integral MD: :doc:`first variant <fix_pimd>`, :doc:`second variant <fix_ipi>`
 * multi-walker collective variables with :doc:`Colvars <fix_colvars>` and :doc:`Plumed <fix_plumed>`
 
 .. _prepost:

doc/src/Modify.rst

@@ -34,5 +34,6 @@ style requirements and recommendations <Modify_style>`.
    Modify_min
    Modify_region
    Modify_body
+   Modify_gran_sub_mod
    Modify_thermo
    Modify_variable

doc/src/Modify_gran_sub_mod.rst

@@ -0,0 +1,177 @@
+Granular Sub-Model styles
+===============================
+
+In granular models, particles are spheres with a finite radius and rotational
+degrees of freedom as further described in the
+:doc:`Howto granular page <Howto_granular>`. Interactions between pair of
+particles or particles and walls may therefore depend on many different modes
+of motion as described in :doc:`pair granular <pair_granular>` and
+:doc:`fix wall/gran <fix_wall_gran>`. In both cases, the exchange of forces,
+torques, and heat flow between two types of bodies is defined using a
+GranularModel class. The GranularModel class organizes the details of an
+interaction using a series of granular sub-models each of which describe a
+particular interaction mode (e.g. normal forces or rolling friction). From a
+parent GranSubMod class, several types of sub-model classes are derived:
+
+* GranSubModNormal: normal force sub-model
+* GranSubModDamping: normal damping sub-model
+* GranSubModTangential: tangential forces and sliding friction sub-model
+* GranSubModRolling: rolling friction sub-model
+* GranSubModTwisting: twisting friction sub-model
+* GranSubModHeat: heat conduction sub-model
+
+For each type of sub-model, more classes are further derived, each describing a
+specific implementation. For instance, from the GranSubModNormal class the
+GranSubModNormalHooke, GranSubModNormalHertz, and GranSubModNormalJKR classes
+are derived which calculate Hookean, Hertzian, or JKR normal forces,
+respectively.  This modular structure simplifies the addition of new granular
+contact models as one only needs to create a new GranSubMod class without
+having to modify the more complex PairGranular, FixGranWall, and GranularModel
+classes. Most GranSubMod methods are also already defined by the parent classes
+so new contact models typically only require edits to a few relevant methods
+(e.g. methods that define coefficients and calculate forces).
+
+Each GranSubMod class has a pointer to both the LAMMPS class and the GranularModel
+class which owns it, ``lmp`` and ``gm``, respectively. The GranularModel class
+includes several public variables that describe the geometry/dynamics of the
+contact such as
+
+.. list-table::
+
+   * - ``xi`` and ``xj``
+     - Positions of the two contacting bodies
+   * - ``vi`` and ``vj``
+     - Velocities of the two contacting bodies
+   * - ``omegai`` and ``omegaj``
+     - Angular velocities of the two contacting bodies
+   * - ``dx`` and ``nx``
+     - The displacement and normalized displacement vectors
+   * - ``r``, ``rsq``, and ``rinv``
+     - The distance, distance squared, and inverse distance
+   * - ``radsum``
+     - The sum of particle radii
+   * - ``vr``, ``vn``, and ``vt``
+     - The relative velocity vector and its normal and tangential components
+   * - ``wr``
+     - The relative rotational velocity
+
+These quantities, among others, are calculated in the ``GranularModel->check_contact()``
+and ``GranularModel->calculate_forces()`` methods which can be referred to for more
+details.
+
+To create a new GranSubMod class, it is recommended that one first looks at similar
+GranSubMod classes. All GranSubMod classes share several general methods which may
+need to be defined
+
+.. list-table::
+
+   * - ``GranSubMod->mix_coeff()``
+     - Optional method to define how coefficients are mixed for different atom types. By default, coefficients are mixed using a geometric mean.
+   * - ``GranSubMod->coeffs_to_local()``
+     - Parses coefficients to define local variables. Run once at model construction.
+   * - ``GranSubMod->init()``
+     - Optional method to define local variables after other GranSubMod types were created. For instance, this method may be used by a tangential model that derives parameters from the normal model.
+
+There are also several type-specific methods
+
+.. list-table::
+
+   * - ``GranSubModNormal->touch()``
+     - Optional method to test when particles are in contact. By default, this is when particles overlap.
+   * - ``GranSubModNormal->pulloff_distance()``
+     - Optional method to return the distance at which particles stop interacting. By default, this is when particles no longer overlap.
+   * - ``GranSubModNormal->calculate_area()``
+     - Optional method to return the surface area of the contact. By default, this returns the geometric cross section.
+   * - ``GranSubModNormal->set_fncrit()``
+     - Optional method that defines the critical force to break the contact used by some tangential, rolling, and twisting sub-models. By default, this is the current total normal force including damping.
+   * - ``GranSubModNormal->calculate_forces()``
+     - Required method that returns the normal contact force
+   * - ``GranSubModDamping->calculate_forces()``
+     - Required method that returns the normal damping force
+   * - ``GranSubModTangential->calculate_forces()``
+     - Required method that calculates tangential forces/torques
+   * - ``GranSubModTwisting->calculate_forces()``
+     - Required method that calculates twisting friction forces/torques
+   * - ``GranSubModRolling->calculate_forces()``
+     - Required method that calculates rolling friction forces/torques
+   * - ``GranSubModHeat->calculate_heat()``
+     - Required method that returns the rate of heat flow
+
+As an example, say one wanted to create a new normal force option that consisted
+of a Hookean force with a piecewise stiffness. This could be done by adding a new
+set of files ``gran_sub_mod_custom.h``:
+
+.. code-block:: c++
+
+   #ifdef GranSubMod_CLASS
+   // clang-format off
+   GranSubModStyle(hooke/piecewise,
+            GranSubModNormalHookePiecewise,
+            NORMAL);
+   // clang-format on
+   #else
+
+   #ifndef GRAN_SUB_MOD_CUSTOM_H_
+   #define GRAN_SUB_MOD_CUSTOM_H_
+
+   #include "gran_sub_mod.h"
+   #include "gran_sub_mod_normal.h"
+
+   namespace LAMMPS_NS {
+   namespace Granular_NS {
+
+   class GranSubModNormalHookePiecewise : public GranSubModNormal {
+    public:
+     GranSubModNormalHookePiecewise(class GranularModel *, class LAMMPS *);
+     void coeffs_to_local() override;
+     double calculate_forces();
+    protected:
+     double k1, k2, delta_switch;
+   };
+
+   }    // namespace Granular_NS
+   }    // namespace LAMMPS_NS
+
+   #endif /*GRAN_SUB_MOD_CUSTOM_H_ */
+   #endif /*GRAN_SUB_MOD_CLASS_H_ */
+
+
+and ``gran_sub_mod_custom.cpp``
+
+.. code-block:: c++
+
+   #include "gran_sub_mod_custom.h"
+   #include "gran_sub_mod_normal.h"
+   #include "granular_model.h"
+
+   using namespace LAMMPS_NS;
+   using namespace Granular_NS;
+
+   GranSubModNormalHookePiecewise::GranSubModNormalHookePiecewise(GranularModel *gm, LAMMPS *lmp) :  GranSubModNormal(gm, lmp)
+   {
+     num_coeffs = 4;
+   }
+
+   /* ---------------------------------------------------------------------- */
+
+   void GranSubModNormalHookePiecewise::coeffs_to_local()
+   {
+     k1 = coeffs[0];
+     k2 = coeffs[1];
+     damp = coeffs[2];
+     delta_switch = coeffs[3];
+   }
+
+   /* ---------------------------------------------------------------------- */
+
+   double GranSubModNormalHookePiecewise::calculate_forces()
+   {
+     double Fne;
+     if (gm->delta >= delta_switch) {
+       Fne = k1 * delta_switch + k2 * (gm->delta - delta_switch);
+     } else {
+       Fne = k1 * gm->delta;
+     }
+     return Fne;
+   }
+

doc/src/Modify_pair.rst

@@ -2,11 +2,11 @@ Pair styles
 ===========
 
 Classes that compute pairwise non-bonded interactions are derived from
-the Pair class.  In LAMMPS, pairwise calculation include many-body
+the ``Pair`` class.  In LAMMPS, pairwise force calculations include
-potentials such as EAM, Tersoff, or ReaxFF where particles interact
+many-body potentials such as EAM, Tersoff, or ReaxFF where particles
-without an explicit bond topology but include interactions beyond
+interact without an explicit bond topology but include interactions
-pairwise non-bonded contributions.  New styles can be created to add
+beyond pairwise non-bonded contributions.  New styles can be created to
-support for additional pair potentials to LAMMPS.  When the
+add support for additional pair potentials to LAMMPS.  When the
 modifications are small, sometimes it is more effective to derive from
 an existing pair style class.  This latter approach is also used by
 :doc:`Accelerator packages <Speed_packages>` where the accelerated style
@@ -15,10 +15,13 @@ names differ from their base classes by an appended suffix.
 The file ``src/pair_lj_cut.cpp`` is an example of a Pair class with a
 very simple potential function.  It includes several optional methods to
 enable its use with :doc:`run_style respa <run_style>` and :doc:`compute
-group/group <compute_group_group>`.
+group/group <compute_group_group>`.  :doc:`Developer_write_pair` contains
+a detailed discussion of writing new pair styles from scratch, and how
+simple and more complex pair styles can be implemented with examples
+from existing pair styles.
 
 Here is a brief list of some the class methods in the Pair class that
-*must* be or *may* be overridden in a derived class.
+*must* be or *may* be overridden in a derived class for a new pair style.
 
 +---------------------------------+---------------------------------------------------------------------+
 | Required                        | "pure" methods that *must* be overridden in a derived class         |

doc/src/Python_install.rst

@@ -27,140 +27,36 @@ interpreter can find it and installing the LAMMPS shared library into a
 folder that the dynamic loader searches or inside of the installed
 ``lammps`` package folder.  There are multiple ways to achieve this.
 
-#. Do a full LAMMPS installation of libraries, executables, selected
-   headers, documentation (if enabled), and supporting files (only
-   available via CMake), which can also be either system-wide or into
-   user specific folders.
-
 #. Install both components into a Python ``site-packages`` folder, either
    system-wide or in the corresponding user-specific folder. This way no
    additional environment variables need to be set, but the shared
    library is otherwise not accessible.
 
-#. Do an installation into a virtual environment. This can either be an
+#. Do an installation into a virtual environment.
-   installation of the Python package only or a full installation of LAMMPS.
 
 #. Leave the files where they are in the source/development tree and
    adjust some environment variables.
 
 .. tabs::
 
-   .. tab:: Full install (CMake-only)
+   .. tab:: Python package
-
-      :ref:`Build the LAMMPS executable and library <library>` with
-      ``-DBUILD_SHARED_LIBS=on``, ``-DLAMMPS_EXCEPTIONS=on`` and
-      ``-DPKG_PYTHON=on`` (The first option is required, the other two
-      are optional by recommended).  The exact file name of the shared
-      library depends on the platform (Unix/Linux, macOS, Windows) and
-      the build configuration being used.  The installation base folder
-      is already set by default to the ``$HOME/.local`` directory, but
-      it can be changed to a custom location defined by the
-      ``CMAKE_INSTALL_PREFIX`` CMake variable.  This uses a folder
-      called ``build`` to store files generated during compilation.
-
-      .. code-block:: bash
-
-         # create build folder
-         mkdir build
-         cd build
-
-         # configure LAMMPS compilation
-         cmake -C ../cmake/presets/basic.cmake -D BUILD_SHARED_LIBS=on \
-               -D LAMMPS_EXCEPTIONS=on -D PKG_PYTHON=on ../cmake
-
-         # compile LAMMPS
-         cmake --build .
-
-         # install LAMMPS into $HOME/.local
-         cmake --install .
-
-
-      This leads to an installation to the following locations:
-
-      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+
-      | File                   | Location                                                        | Notes                                                       |
-      +========================+=================================================================+=============================================================+
-      | LAMMPS Python package  | * ``$HOME/.local/lib/pythonX.Y/site-packages/lammps`` (32bit)   | ``X.Y`` depends on the installed Python version             |
-      |                        | * ``$HOME/.local/lib64/pythonX.Y/site-packages/lammps`` (64bit) |                                                             |
-      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS shared library  | * ``$HOME/.local/lib/`` (32bit)                                 | Set shared loader environment variable to this path         |
-      |                        | * ``$HOME/.local/lib64/`` (64bit)                               | (see below for more info on this)                           |
-      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS executable      | * ``$HOME/.local/bin/``                                         |                                                             |
-      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS potential files | * ``$HOME/.local/share/lammps/potentials/``                     | Set ``LAMMPS_POTENTIALS`` environment variable to this path |
-      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+
-
-      For a system-wide installation you need to set
-      ``CMAKE_INSTALL_PREFIX`` to a system folder like ``/usr`` (or
-      ``/usr/local``); the default is ``${HOME}/.local``.  The
-      installation step for a system folder installation (**not** the
-      configuration/compilation) needs to be done with superuser
-      privilege, e.g. by using ``sudo cmake --install .``.  The
-      installation folders will then be changed to (assuming ``/usr`` as
-      prefix):
-
-      +------------------------+---------------------------------------------------------+-------------------------------------------------------------+
-      | File                   | Location                                                | Notes                                                       |
-      +========================+=========================================================+=============================================================+
-      | LAMMPS Python package  | * ``/usr/lib/pythonX.Y/site-packages/lammps`` (32bit)   | ``X.Y`` depends on the installed Python version             |
-      |                        | * ``/usr/lib64/pythonX.Y/site-packages/lammps`` (64bit) |                                                             |
-      +------------------------+---------------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS shared library  | * ``/usr/lib/`` (32bit)                                 |                                                             |
-      |                        | * ``/usr/lib64/`` (64bit)                               |                                                             |
-      +------------------------+---------------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS executable      | * ``/usr/bin/``                                         |                                                             |
-      +------------------------+---------------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS potential files | * ``/usr/share/lammps/potentials/``                     |                                                             |
-      +------------------------+---------------------------------------------------------+-------------------------------------------------------------+
-
-      To be able to use the "user" installation you have to ensure that
-      the folder containing the LAMMPS shared library is either included
-      in a path searched by the shared linker (e.g. like
-      ``/usr/lib64/``) or part of the ``LD_LIBRARY_PATH`` environment
-      variable (or ``DYLD_LIBRARY_PATH`` on macOS).  Otherwise you will
-      get an error when trying to create a LAMMPS object through the
-      Python module.
-
-      .. code-block:: bash
-
-         # Unix/Linux
-         export LD_LIBRARY_PATH=$HOME/.local/lib:$LD_LIBRARY_PATH
-
-         # macOS
-         export DYLD_LIBRARY_PATH=$HOME/.local/lib:$DYLD_LIBRARY_PATH
-
-      If you plan to use the LAMMPS executable (e.g., ``lmp``), you may
-      also need to adjust the ``PATH`` environment variable (but many
-      newer Linux distributions already have ``$HOME/.local/bin``
-      included). Example:
-
-      .. code-block:: bash
-
-         export PATH=$HOME/.local/bin:$PATH
-
-      To make those changes permanent, you can add the commands to your
-      ``$HOME/.bashrc`` file.  For a system-wide installation is is not
-      necessary due to files installed in system folders that are loaded
-      automatically when a login shell is started.
-
-   .. tab:: Python package only
 
       Compile LAMMPS with either :doc:`CMake <Build_cmake>` or the
       :doc:`traditional make <Build_make>` procedure in :ref:`shared
-      mode <exe>`.  After compilation has finished type (in the
+      mode <exe>`.  After compilation has finished, type (in the
       compilation folder):
 
       .. code-block:: bash
 
          make install-python
 
-      This will try to build a so-called (binary) 'wheel', a compressed
+      This will try to build a so-called (binary) wheel file, a
-      binary python package and then install it with the python package
+      compressed binary python package and then install it with the
-      manager 'pip'.  Installation will be attempted into a system-wide
+      python package manager 'pip'.  Installation will be attempted into
-      ``site-packages`` folder and if that fails into the corresponding
+      a system-wide ``site-packages`` folder and if that fails into the
-      folder in the user's home directory.  For a system-wide installation you
+      corresponding folder in the user's home directory.  For a
-      would have to gain superuser privilege, e.g. though ``sudo``
+      system-wide installation you usually would have to gain superuser
+      privilege first, e.g. though ``sudo``
 
       +------------------------+----------------------------------------------------------+-------------------------------------------------------------+
       | File                   | Location                                                 | Notes                                                       |
@@ -211,10 +107,11 @@ folder that the dynamic loader searches or inside of the installed
 
       .. code-block:: bash
 
-         python install.py -p <python package> -l <shared library> [-n]
+         python install.py -p <python package> -l <shared library> -v <version.h file> [-n]
 
       * The ``-p`` flag points to the ``lammps`` Python package folder to be installed,
       * the ``-l`` flag points to the LAMMPS shared library file to be installed,
+      * the ``-v`` flag points to the LAMMPS version header file to extract the version date,
       * and the optional ``-n`` instructs the script to only build a wheel file
         but not attempt to install it.
 
@@ -272,38 +169,6 @@ folder that the dynamic loader searches or inside of the installed
       | LAMMPS shared library  | * ``$VIRTUAL_ENV/lib/pythonX.Y/site-packages/lammps``  | ``X.Y`` depends on the installed Python version             |
       +------------------------+--------------------------------------------------------+-------------------------------------------------------------+
 
-      If you do a full installation (CMake only) with "install", this
-      leads to the following installation locations:
-
-      +------------------------+--------------------------------------------------------+-------------------------------------------------------------+
-      | File                   | Location                                               | Notes                                                       |
-      +========================+========================================================+=============================================================+
-      | LAMMPS Python Module   | * ``$VIRTUAL_ENV/lib/pythonX.Y/site-packages/lammps``  | ``X.Y`` depends on the installed Python version             |
-      +------------------------+--------------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS shared library  | * ``$VIRTUAL_ENV/lib/`` (32bit)                        | Set shared loader environment variable to this path         |
-      |                        | * ``$VIRTUAL_ENV/lib64/`` (64bit)                      | (see below for more info on this)                           |
-      +------------------------+--------------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS executable      | * ``$VIRTUAL_ENV/bin/``                                |                                                             |
-      +------------------------+--------------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS potential files | * ``$VIRTUAL_ENV/share/lammps/potentials/``            | Set ``LAMMPS_POTENTIALS`` environment variable to this path |
-      +------------------------+--------------------------------------------------------+-------------------------------------------------------------+
-
-      In that case you need to modify the ``$HOME/myenv/bin/activate``
-      script in a similar fashion you need to update your
-      ``$HOME/.bashrc`` file to include the shared library and
-      executable locations in ``LD_LIBRARY_PATH`` (or
-      ``DYLD_LIBRARY_PATH`` on macOS) and ``PATH``, respectively.
-
-      For example with:
-
-      .. code-block:: bash
-
-         # Unix/Linux
-         echo 'export LD_LIBRARY_PATH=$VIRTUAL_ENV/lib:$LD_LIBRARY_PATH' >> $HOME/myenv/bin/activate
-
-         # macOS
-         echo 'export DYLD_LIBRARY_PATH=$VIRTUAL_ENV/lib:$DYLD_LIBRARY_PATH' >> $HOME/myenv/bin/activate
-
    .. tab:: In place usage
 
       You can also :doc:`compile LAMMPS <Build>` as usual in

doc/src/Tools.rst

@@ -320,7 +320,8 @@ eam generate tool
 -----------------------------
 
 The tools/eam_generate directory contains several one-file C programs
-that convert an analytic formula into a tabulated :doc:`embedded atom method (EAM) <pair_eam>` setfl potential file.  The potentials they
+that convert an analytic formula into a tabulated :doc:`embedded atom
+method (EAM) <pair_eam>` setfl potential file.  The potentials they
 produce are in the potentials directory, and can be used with the
 :doc:`pair_style eam/alloy <pair_eam>` command.
 

doc/src/atom_style.rst

@@ -161,15 +161,14 @@ and each stores a per-particle diameter and mass.  If the diameter >
 0.0, the particle is a finite-size sphere.  If the diameter = 0.0, it
 is a point particle.  Note that by use of the *disc* keyword with the
 :doc:`fix nve/sphere <fix_nve_sphere>`, :doc:`fix nvt/sphere
-<fix_nvt_sphere>`, :doc:`fix nph/sphere <fix_nph_sphere>`, :doc:`fix
+<fix_nvt_sphere>`, :doc:`fix nph/sphere <fix_nph_sphere>`,
-npt/sphere <fix_npt_sphere>` commands for the *sphere* style, spheres
+:doc:`fix npt/sphere <fix_npt_sphere>` commands for the *sphere* style,
-can be effectively treated as 2d discs for a 2d simulation if desired.
+spheres can be effectively treated as 2d discs for a 2d simulation if
-See also the :doc:`set density/disc <set>` command.  The *sphere* and
+desired.  See also the :doc:`set density/disc <set>` command.  These
-*bpm/sphere* styles take an optional 0 or 1 argument.  A value of 0
+styles take an optional 0 or 1 argument.  A value of 0 means the
-means the radius of each sphere is constant for the duration of the
+radius of each sphere is constant for the duration of the simulation.
-simulation.  A value of 1 means the radii may vary dynamically during
+A value of 1 means the radii may vary dynamically during the simulation,
-the simulation, e.g. due to use of the :doc:`fix adapt <fix_adapt>`
+e.g. due to use of the :doc:`fix adapt <fix_adapt>` command.
-command.
 
 For the *ellipsoid* style, the particles are ellipsoids and each
 stores a flag which indicates whether it is a finite-size ellipsoid or

doc/src/balance.rst

@@ -53,6 +53,7 @@ Syntax
              name = name of the atom-style variable
            *store* name = store weight in custom atom property defined by :doc:`fix property/atom <fix_property_atom>` command
              name = atom property name (without d\_ prefix)
+       *sort* arg = *no* or *yes*
        *out* arg = filename
          filename = write each processor's subdomain to a file
 
@@ -492,6 +493,14 @@ different kinds of custom atom vectors or arrays as arguments.
 
 ----------
 
+The *sort* keyword determines whether the communication of per-atom
+data to other processors during load-balancing will be random or
+deterministic.  Random is generally faster; deterministic will ensure
+the new ordering of atoms on each processor is the same each time the
+same simulation is run.  This can be useful for debugging purposes.
+Since the balance command is a one-time operation, the default is
+*yes* to perform sorting.
+
 The *out* keyword writes a text file to the specified *filename* with
 the results of the balancing operation.  The file contains the bounds
 of the subdomain for each processor after the balancing operation
@@ -569,4 +578,5 @@ Related commands
 Default
 """""""
 
-none
+The default setting is sort = yes.
+

doc/src/bond_bpm_rotational.rst

@@ -10,7 +10,7 @@ Syntax
 
    bond_style bpm/rotational keyword value attribute1 attribute2 ...
 
-* optional keyword = *overlay/pair* or *store/local* or *smooth*
+* optional keyword = *overlay/pair* or *store/local* or *smooth* or *break/no*
 
   .. parsed-literal::
 
@@ -30,6 +30,9 @@ Syntax
        *smooth* value = *yes* or *no*
           smooths bond forces near the breaking point
 
+       *break/no*
+          indicates that bonds should not break during a run
+
 Examples
 """"""""
 
@@ -140,6 +143,12 @@ the *overlay/pair* keyword. These settings require specific
 restrictions.  Further details can be found in the `:doc: how to
 <Howto_BPM>` page on BPMs.
 
+.. versionadded:: 28Mar2023
+
+If the *break/no* keyword is used, then LAMMPS assumes bonds should not break
+during a simulation run. This will prevent some unnecessary calculation.
+However, if a bond does break, it will trigger an error.
+
 If the *store/local* keyword is used, an internal fix will track bonds that
 break during the simulation. Whenever a bond breaks, data is processed
 and transferred to an internal fix labeled *fix_ID*. This allows the

doc/src/bond_bpm_spring.rst

@@ -10,7 +10,7 @@ Syntax
 
    bond_style bpm/spring keyword value attribute1 attribute2 ...
 
-* optional keyword = *overlay/pair* or *store/local* or *smooth*
+* optional keyword = *overlay/pair* or *store/local* or *smooth* or *break/no*
 
   .. parsed-literal::
 
@@ -30,6 +30,9 @@ Syntax
        *smooth* value = *yes* or *no*
           smooths bond forces near the breaking point
 
+       *break/no*
+          indicates that bonds should not break during a run
+
 Examples
 """"""""
 
@@ -47,7 +50,7 @@ Description
 
 .. versionadded:: 4May2022
 
-The *bpm/spring* bond style computes forces and torques based on
+The *bpm/spring* bond style computes forces based on
 deviations from the initial reference state of the two atoms.  The
 reference state is stored by each bond when it is first computed in
 the setup of a run. Data is then preserved across run commands and is
@@ -56,7 +59,8 @@ the system will not reset the reference state of a bond.
 
 This bond style only applies central-body forces which conserve the
 translational and rotational degrees of freedom of a bonded set of
-particles. The force has a magnitude of
+particles based on a model described by Clemmer and Robbins
+:ref:`(Clemmer) <fragment-Clemmer>`. The force has a magnitude of
 
 .. math::
 
@@ -105,6 +109,12 @@ the *overlay/pair* keyword. These settings require specific
 restrictions.  Further details can be found in the `:doc: how to
 <Howto_BPM>` page on BPMs.
 
+.. versionadded:: 28Mar2023
+
+If the *break/no* keyword is used, then LAMMPS assumes bonds should not break
+during a simulation run. This will prevent some unnecessary calculation.
+However, if a bond does break, it will trigger an error.
+
 If the *store/local* keyword is used, an internal fix will track bonds that
 break during the simulation. Whenever a bond breaks, data is processed
 and transferred to an internal fix labeled *fix_ID*. This allows the
@@ -200,6 +210,10 @@ The option defaults are *smooth* = *yes*
 
 ----------
 
+.. _fragment-Clemmer:
+
+**(Clemmer)** Clemmer and Robbins, Phys. Rev. Lett. (2022).
+
 .. _Groot4:
 
 **(Groot)** Groot and Warren, J Chem Phys, 107, 4423-35 (1997).

doc/src/bond_fene.rst

@@ -51,7 +51,7 @@ in the same form as in pair style :doc:`nm/cut <pair_nm>`. The bond energy is th
 
 .. math::
 
-  E = -0.5 K r_0^2  \ln \left[ 1 - \left(\frac{r}{R_0}\right)^2\right] + \frac{E_0}{(n-m)} \left[ m \left(\frac{r_0}{r}\right)^n - n \left(\frac{r_0}{r}\right)^m \right]
+  E = -0.5 K R_0^2  \ln \left[ 1 - \left(\frac{r}{R_0}\right)^2\right] + \frac{E_0}{(n-m)} \left[ m \left(\frac{r_0}{r}\right)^n - n \left(\frac{r_0}{r}\right)^m \right]
 
 Similar to the *fene* style, the generalized Lennard-Jones is cut off at
 the potential minimum, :math:`r_0`, to be repulsive only.  The following

doc/src/bond_harmonic_restrain.rst

@@ -0,0 +1,90 @@
+.. index:: bond_style harmonic/restrain
+
+bond_style harmonic/restrain command
+====================================
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   bond_style harmonic/restrain
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   bond_style harmonic
+   bond_coeff 5 80.0
+
+Description
+"""""""""""
+
+.. versionadded:: 28Mar2023
+
+The *harmonic/restrain* bond style uses the potential
+
+.. math::
+
+   E = K (r - r_{t=0})^2
+
+where :math:`r_{t=0}` is the distance between the bonded atoms at the
+beginning of the first :doc:`run <run>` or :doc:`minimize <minimize>`
+command after the bond style has been defined (*t=0*).  Note that the
+usual 1/2 factor is included in :math:`K`.  This will effectively
+restrain bonds to their initial length, whatever that is.  This is where
+this bond style differs from :doc:`bond style harmonic <bond_harmonic>`
+where the bond length is set through the per bond type coefficients.
+
+The following coefficient must be defined for each bond type via the
+:doc:`bond_coeff <bond_coeff>` command as in the example above, or in
+the data file or restart files read by the :doc:`read_data <read_data>`
+or :doc:`read_restart <read_restart>` commands
+
+* :math:`K` (energy/distance\^2)
+
+This bond style differs from other options to add harmonic restraints
+like :doc:`fix restrain <fix_restrain>` or :doc:`pair style list
+<pair_list>` or :doc:`fix colvars <fix_colvars>` in that it requires a
+bond topology, and thus the defined bonds will trigger exclusion of
+special neighbors from the neighbor list according to the
+:doc:`special_bonds <special_bonds>` settings.
+
+Restart info
+""""""""""""
+
+This bond style supports the :doc:`write_restart <write_restart>` and
+:doc:`read_restart <read_restart>` commands. The state of the initial
+bond lengths is stored with restart files and read back.
+
+Restrictions
+""""""""""""
+
+This bond style can only be used if LAMMPS was built with the
+EXTRA-MOLECULE package.  See the :doc:`Build package <Build_package>`
+page for more info.
+
+This bond style maintains internal data to determine the original bond
+lengths :math:`r_{t=0}`.  This information will be written to
+:doc:`binary restart files <write_restart>` but **not** to :doc:`data
+files <write_data>`.  Thus, continuing a simulation is *only* possible
+with :doc:`read_restart <read_restart>`.  When using the :doc:`read_data
+command <read_data>`, the reference bond lengths :math:`r_{t=0}` will be
+re-initialized from the current geometry.
+
+This bond style cannot be used with :doc:`fix shake or fix rattle
+<fix_shake>`, with :doc:`fix filter/corotate <fix_filter_corotate>`, or
+any :doc:`tip4p pair style <pair_lj_cut_tip4p>` since there is no specific
+equilibrium distance for a given bond type.
+
+Related commands
+""""""""""""""""
+
+:doc:`bond_coeff <bond_coeff>`,  :doc:`bond_harmonic <bond_harmonic>`,
+:doc:`fix restrain <fix_restrain>`, :doc:`pair style list <pair_list>`
+
+Default
+"""""""
+
+none

doc/src/bond_style.rst

@@ -10,7 +10,7 @@ Syntax
 
    bond_style style args
 
-* style = *none* or *zero* or *hybrid* or *bpm/rotational* or *bpm/spring* or *class2* or *fene* or *fene/expand* or *fene/nm* or *gaussian* or *gromos* or *harmonic* or *harmonic/shift* or *harmonic/shift/cut* or *lepton* or *morse* or *nonlinear* or *oxdna/fene* or *oxdena2/fene* or *oxrna2/fene* or *quartic* or *special* or *table*
+* style = *none* or *zero* or *hybrid* or *bpm/rotational* or *bpm/spring* or *class2* or *fene* or *fene/expand* or *fene/nm* or *gaussian* or *gromos* or *harmonic* or *harmonic/restrain* *harmonic/shift* or *harmonic/shift/cut* or *lepton* or *morse* or *nonlinear* or *oxdna/fene* or *oxdena2/fene* or *oxrna2/fene* or *quartic* or *special* or *table*
 
 * args = none for any style except *hybrid*
 
@@ -93,6 +93,7 @@ accelerated styles exist.
 * :doc:`gaussian <bond_gaussian>` - multicentered Gaussian-based bond potential
 * :doc:`gromos <bond_gromos>` - GROMOS force field bond
 * :doc:`harmonic <bond_harmonic>` - harmonic bond
+* :doc:`harmonic/restrain <bond_harmonic_restrain>` - harmonic bond to restrain to original bond distance
 * :doc:`harmonic/shift <bond_harmonic_shift>` - shifted harmonic bond
 * :doc:`harmonic/shift/cut <bond_harmonic_shift_cut>` - shifted harmonic bond with a cutoff
 * :doc:`lepton <bond_lepton>` - bond potential from evaluating a string

doc/src/bond_table.rst

@@ -112,8 +112,9 @@ are estimated (less accurately) by the first two and last two force
 values in the table.
 
 The "EQ" parameter is also optional.  If used, it is followed by a the
-equilibrium bond length, which is used, for example, by the :doc:`fix shake <fix_shake>` command.  If not used, the equilibrium bond
+equilibrium bond length, which is used, for example, by the :doc:`fix
-length is to the distance in the table with the lowest potential energy.
+shake <fix_shake>` command.  If not used, the equilibrium bond length is
+to the distance in the table with the lowest potential energy.
 
 Following a blank line, the next N lines list the tabulated values.
 On each line, the first value is the index from 1 to N, the second value is
@@ -135,16 +136,15 @@ one that matches the specified keyword.
 
 ----------
 
-Restart, fix_modify, output, run start/stop, minimize info
+Restart info
-"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+""""""""""""
 
-This bond style writes the settings for the "bond_style table"
+This bond style writes the settings for the "bond_style table" command
-command to :doc:`binary restart files <restart>`, so a bond_style
+to :doc:`binary restart files <restart>`, so a bond_style command does
-command does not need to specified in an input script that reads a
+not need to specified in an input script that reads a restart file.
-restart file.  However, the coefficient information is not stored in
+However, the coefficient information is not stored in the restart file,
-the restart file, since it is tabulated in the potential files.  Thus,
+since it is tabulated in the potential files.  Thus, bond_coeff commands
-bond_coeff commands do need to be specified in the restart input
+do need to be specified in the restart input script.
-script.
 
 Restrictions
 """"""""""""

doc/src/bond_write.rst

@@ -70,7 +70,7 @@ be specified even if the potential has a finite value at r = 0.0.
 Related commands
 """"""""""""""""
 
-:doc:`bond_style table <bond_table>`, `angle_write <angle_write>`,
+:doc:`bond_style table <bond_table>`, :doc:`angle_write <angle_write>`,
 :doc:`bond_style <bond_style>`, :doc:`bond_coeff <bond_coeff>`
 
 Default

doc/src/comm_modify.rst

@@ -69,6 +69,7 @@ For many systems this is an efficient algorithm, but for systems with
 widely varying cutoffs for different type pairs, the *multi* or *multi/old* mode can
 be faster.  In *multi*, each atom is assigned to a collection which should
 correspond to a set of atoms with similar interaction cutoffs.
+See the :doc:`neighbor <neighbor>` command for a detailed description of collections.
 In this case, each atom collection is assigned its own distance
 cutoff for communication purposes, and fewer atoms will be
 communicated. in *multi/old*, a similar technique is used but atoms

doc/src/compute.rst

@@ -206,6 +206,8 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
 * :doc:`dilatation/atom <compute_dilatation_atom>` - Peridynamic dilatation for each atom
 * :doc:`dipole <compute_dipole>` - dipole vector and total dipole
 * :doc:`dipole/chunk <compute_dipole_chunk>` - dipole vector and total dipole for each chunk
+* :doc:`dipole/tip4p <compute_dipole>` - dipole vector and total dipole with TIP4P pair style
+* :doc:`dipole/tip4p/chunk <compute_dipole_chunk>` - dipole vector and total dipole for each chunk with TIP4P pair style
 * :doc:`displace/atom <compute_displace_atom>` - displacement of each atom
 * :doc:`dpd <compute_dpd>` - total values of internal conductive energy, internal mechanical energy, chemical energy, and harmonic average of internal temperature
 * :doc:`dpd/atom <compute_dpd_atom>` - per-particle values of internal conductive energy, internal mechanical energy, chemical energy, and internal temperature
@@ -258,6 +260,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
 * :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:`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:`property/atom <compute_property_atom>` - convert atom attributes to per-atom vectors/arrays
 * :doc:`property/chunk <compute_property_chunk>` - extract various per-chunk attributes

doc/src/compute_dipole.rst

@@ -1,6 +1,10 @@
 .. index:: compute dipole
+.. index:: compute dipole/tip4p
 
 compute dipole command
+======================
+
+compute dipole/tip4p command
 ============================
 
 Syntax
@@ -8,10 +12,10 @@ Syntax
 
 .. code-block:: LAMMPS
 
-   compute ID group-ID dipole arg
+   compute ID group-ID style arg
 
 * ID, group-ID are documented in :doc:`compute <compute>` command
-* dipole = style name of this compute command
+* style = *dipole* or *dipole/tip4p*
 * arg = *mass* or *geometry* = use COM or geometric center for charged chunk correction (optional)
 
 Examples
@@ -21,6 +25,7 @@ Examples
 
    compute 1 fluid dipole
    compute dw water dipole geometry
+   compute dw water dipole/tip4p
 
 Description
 """""""""""
@@ -28,13 +33,20 @@ Description
 Define a computation that calculates the dipole vector and total dipole
 for a group of atoms.
 
-This compute calculates the x,y,z coordinates of the dipole vector
+These computes calculate the x,y,z coordinates of the dipole vector and
-and the total dipole moment for the atoms in the compute group.
+the total dipole moment for the atoms in the compute group.  This
-This includes all effects due to atoms passing through periodic boundaries.
+includes all effects due to atoms passing through periodic boundaries.
-For a group with a net charge the resulting dipole is made position independent
+For a group with a net charge the resulting dipole is made position
-by subtracting the position vector of the center of mass or geometric center
+independent by subtracting the position vector of the center of mass or
-times the net charge from the computed dipole vector. Both per-atom charges
+geometric center times the net charge from the computed dipole
-and per-atom dipole moments, if present, contribute to the computed dipole.
+vector.  Both per-atom charges and per-atom dipole moments, if present,
+contribute to the computed dipole.
+
+.. versionadded:: 28Mar2023
+
+Compute *dipole/tip4p* includes adjustments for the charge carrying
+point M in molecules with TIP4P water geometry.  The corresponding
+parameters are extracted from the pair style.
 
 .. note::
 
@@ -49,10 +61,10 @@ and per-atom dipole moments, if present, contribute to the computed dipole.
 Output info
 """""""""""
 
-This compute calculations a global scalar containing the magnitude of
+These computes calculate a global scalar containing the magnitude of the
-the computed dipole moment and a global vector of length 3 with the
+computed dipole moment and a global vector of length 3 with the dipole
-dipole vector.  See the :doc:`Howto output <Howto_output>` page for
+vector.  See the :doc:`Howto output <Howto_output>` page for an overview
-an overview of LAMMPS output options.
+of LAMMPS output options.
 
 The computed values are "intensive".  The array values will be in
 dipole units (i.e., charge :doc:`units <units>` times distance
@@ -60,7 +72,12 @@ dipole units (i.e., charge :doc:`units <units>` times distance
 
 Restrictions
 """"""""""""
- none
+
+Compute style *dipole/tip4p* is part of the EXTRA-COMPUTE package. It is
+only enabled if LAMMPS was built with that package.  See the :doc:`Build
+package <Build_package>` page for more info.
+
+Compute style *dipole/tip4p* can only be used with tip4p pair styles.
 
 Related commands
 """"""""""""""""

doc/src/compute_dipole_chunk.rst

@@ -1,17 +1,21 @@
 .. index:: compute dipole/chunk
+.. index:: compute dipole/tip4p/chunk
 
 compute dipole/chunk command
 ============================
 
+compute dipole/tip4p/chunk command
+==================================
+
 Syntax
 """"""
 
 .. code-block:: LAMMPS
 
-   compute ID group-ID dipole/chunk chunkID arg
+   compute ID group-ID style chunkID arg
 
 * ID, group-ID are documented in :doc:`compute <compute>` command
-* dipole/chunk = style name of this compute command
+* style = *dipole/chunk* or *dipole/tip4p/chunk*
 * chunkID = ID of :doc:`compute chunk/atom <compute_chunk_atom>` command
 * arg = *mass* or *geometry* = use COM or geometric center for charged chunk correction (optional)
 
@@ -38,13 +42,20 @@ or atoms in a spatial bin.  See the :doc:`compute chunk/atom
 details of how chunks can be defined and examples of how they can be
 used to measure properties of a system.
 
-This compute calculates the :math:`(x,y,z)` coordinates of the dipole vector
+These computes calculate the :math:`(x,y,z)` coordinates of the dipole
-and the total dipole moment for each chunk, which includes all effects due
+vector and the total dipole moment for each chunk, which includes all
-to atoms passing through periodic boundaries. For chunks with a net
+effects due to atoms passing through periodic boundaries.  For chunks
-charge the resulting dipole is made position independent by subtracting
+with a net charge the resulting dipole is made position independent by
-the position vector of the center of mass or geometric center times the
+subtracting the position vector of the center of mass or geometric
-net charge from the computed dipole vector. Both per-atom charges and
+center times the net charge from the computed dipole vector.  Both
-per-atom dipole moments, if present, contribute to the computed dipole.
+per-atom charges and per-atom dipole moments, if present, contribute to
+the computed dipole.
+
+.. versionadded:: 28Mar2023
+
+Compute *dipole/tip4p/chunk* includes adjustments for the charge
+carrying point M in molecules with TIP4P water geometry.  The
+corresponding parameters are extracted from the pair style.
 
 Note that only atoms in the specified group contribute to the
 calculation.  The :doc:`compute chunk/atom <compute_chunk_atom>` command
@@ -78,12 +89,12 @@ command, for example:
 Output info
 """""""""""
 
-This compute calculates a global array where the number of rows = the
+These computes calculate a global array where the number of rows = the
 number of chunks *Nchunk* as calculated by the specified :doc:`compute
-chunk/atom <compute_chunk_atom>` command.  The number of columns is 4 for
+chunk/atom <compute_chunk_atom>` command.  The number of columns is 4
-the :math:`(x,y,z)` dipole vector components and the total dipole of each
+for the :math:`(x,y,z)` dipole vector components and the total dipole of
-chunk. These values can be accessed by any command that uses global
+each chunk.  These values can be accessed by any command that uses
-array values from a compute as input.  See the :doc:`Howto output
+global array values from a compute as input.  See the :doc:`Howto output
 <Howto_output>` page for an overview of LAMMPS output options.
 
 The array values are "intensive".  The array values will be in
@@ -92,7 +103,13 @@ dipole units (i.e., charge :doc:`units <units>` times distance
 
 Restrictions
 """"""""""""
- none
+
+Compute style *dipole/tip4p/chunk* is part of the EXTRA-COMPUTE
+package. It is only enabled if LAMMPS was built with that package.  See
+the :doc:`Build package <Build_package>` page for more info.
+
+Compute style *dipole/tip4p/chunk* can only be used with tip4p pair
+styles.
 
 Related commands
 """"""""""""""""

doc/src/compute_pressure_alchemy.rst

@@ -0,0 +1,80 @@
+.. index:: compute pressure/alchemy
+
+compute pressure/alchemy command
+================================
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   compute ID group-ID pressure/alchemy fix-ID
+
+* ID, group-ID are documented in :doc:`compute <compute>` command
+* pressure/alchemy = style name of this compute command
+* fix-ID = ID of :doc:`fix alchemy <fix_alchemy>` command
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   fix trans all alchemy
+   compute mixed all pressure/alchemy trans
+   thermo_modify press mixed
+
+Description
+"""""""""""
+
+.. versionadded:: 28Mar2023
+
+Define a compute style that makes the "mixed" system pressure available
+for a system that uses the :doc:`fix alchemy <fix_alchemy>` command to
+transform one topology to another.  This can be used in combination with
+either :doc:`thermo_modify press <thermo_modify>` or :doc:`fix_modify
+press <fix_modify>` to output and access a pressure consistent with the
+simulated combined two topology system.
+
+The actual pressure is determined with :doc:`compute pressure
+<compute_pressure>` commands that are internally used by :doc:`fix
+alchemy <fix_alchemy>` for each topology individually and then combined.
+This command just extracts the information from the fix.
+
+The ``examples/PACKAGES/alchemy`` folder contains an example input for this command.
+
+----------
+
+Output info
+"""""""""""
+
+This compute calculates a global scalar (the pressure) and a global
+vector of length 6 (the pressure tensor), which can be accessed by
+indices 1--6.  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 ordering of values in the symmetric pressure tensor is as follows:
+:math:`p_{xx},` :math:`p_{yy},` :math:`p_{zz},` :math:`p_{xy},`
+:math:`p_{xz},` :math:`p_{yz}.`
+
+The scalar and vector values calculated by this compute are "intensive".
+The scalar and vector values will be in pressure :doc:`units <units>`.
+
+Restrictions
+""""""""""""
+
+This compute is part of the REPLICA package.  It is only enabled if
+LAMMPS was built with that package.  See the :doc:`Build package
+<Build_package>` page for more info.
+
+
+Related commands
+""""""""""""""""
+
+:doc:`fix alchemy <fix_alchemy>`, :doc:`compute pressure <compute_pressure>`,
+:doc:`thermo_modify <thermo_modify>`, :doc:`fix_modify <fix_modify>`
+
+Default
+"""""""
+
+none

doc/src/dump.rst

@@ -111,6 +111,7 @@ Syntax
                                q, mux, muy, muz, mu,
                                radius, diameter, omegax, omegay, omegaz,
                                angmomx, angmomy, angmomz, tqx, tqy, tqz,
+                               heatflow, temperature,
                                c_ID, c_ID[I], f_ID, f_ID[I], v_name,
                                i_name, d_name, i2_name[I], d2_name[I]
 
@@ -133,10 +134,12 @@ Syntax
            q = atom charge
            mux,muy,muz = orientation of dipole moment of atom
            mu = magnitude of dipole moment of atom
-           radius,diameter = radius,diameter of spherical particle
+           radius,diameter = radius, diameter of spherical particle
            omegax,omegay,omegaz = angular velocity of spherical particle
            angmomx,angmomy,angmomz = angular momentum of aspherical particle
            tqx,tqy,tqz = torque on finite-size particles
+           heatflow = rate of heat flow into particle
+           temperature = temperature of particle
            c_ID = per-atom vector calculated by a compute with ID
            c_ID[I] = Ith column of per-atom array calculated by a compute with ID, I can include wildcard (see below)
            f_ID = per-atom vector calculated by a fix with ID
@@ -375,7 +378,7 @@ output with each snapshot:
    nx ny nz
 
 The value dim will be 2 or 3 for 2d or 3d simulations.  It is included
-so that post-processing tools like `OVITO <https://www.ovito.org>`,
+so that post-processing tools like `OVITO <https://www.ovito.org>`_,
 which can visualize grid-based quantities know how to draw each grid
 cell.  The grid size will match the input script parameters for
 grid(s) created by the computes or fixes which are referenced by the

doc/src/fix.rst

@@ -181,6 +181,7 @@ accelerated styles exist.
 * :doc:`adapt/fep <fix_adapt_fep>` - enhanced version of fix adapt
 * :doc:`addforce <fix_addforce>` - add a force to each atom
 * :doc:`addtorque <fix_addtorque>` - add a torque to a group of atoms
+* :doc:`alchemy <fix_alchemy>` - perform an "alchemical transformation" between two partitions
 * :doc:`amoeba/bitorsion <fix_amoeba_bitorsion>` - torsion/torsion terms in AMOEBA force field
 * :doc:`amoeba/pitorsion <fix_amoeba_pitorsion>` - 6-body terms in AMOEBA force field
 * :doc:`append/atoms <fix_append_atoms>` - append atoms to a running simulation
@@ -221,6 +222,7 @@ accelerated styles exist.
 * :doc:`dt/reset <fix_dt_reset>` - reset the timestep based on velocity, forces
 * :doc:`edpd/source <fix_dpd_source>` - add heat source to eDPD simulations
 * :doc:`efield <fix_efield>` - impose electric field on system
+* :doc:`efield/tip4p <fix_efield>` - impose electric field on system with TIP4P molecules
 * :doc:`ehex <fix_ehex>` - enhanced heat exchange algorithm
 * :doc:`electrode/conp <fix_electrode>` - impose electric potential
 * :doc:`electrode/conq <fix_electrode>` - impose total electric charge
@@ -244,6 +246,7 @@ accelerated styles exist.
 * :doc:`grem <fix_grem>` - implements the generalized replica exchange method
 * :doc:`halt <fix_halt>` - terminate a dynamics run or minimization
 * :doc:`heat <fix_heat>` - add/subtract momentum-conserving heat
+* :doc:`heat/flow <fix_heat_flow>` - plain time integration of heat flow with per-atom temperature updates
 * :doc:`hyper/global <fix_hyper_global>` - global hyperdynamics
 * :doc:`hyper/local <fix_hyper_local>` - local hyperdynamics
 * :doc:`imd <fix_imd>` - implements the "Interactive MD" (IMD) protocol
@@ -258,7 +261,8 @@ accelerated styles exist.
 * :doc:`lb/viscous <fix_lb_viscous>` - :doc:`fix viscous <fix_viscous>` replacement for use with a lattice-Boltzmann fluid
 * :doc:`lineforce <fix_lineforce>` - constrain atoms to move in a line
 * :doc:`manifoldforce <fix_manifoldforce>` - restrain atoms to a manifold during minimization
-* :doc:`mdi/qm <fix_mdi_qm>` - LAMMPS operates as driver for a quantum code via the MolSSI Driver Interface (MDI)
+* :doc:`mdi/qm <fix_mdi_qm>` - LAMMPS operates as a client for a quantum code via the MolSSI Driver Interface (MDI)
+* :doc:`mdi/qmmm <fix_mdi_qmmm>` - LAMMPS operates as client for QM/MM simulation with a quantum code via the MolSSI Driver Interface (MDI)
 * :doc:`meso/move <fix_meso_move>` - move mesoscopic SPH/SDPD particles in a prescribed fashion
 * :doc:`mol/swap <fix_mol_swap>` - Monte Carlo atom type swapping with a molecule
 * :doc:`momentum <fix_momentum>` - zero the linear and/or angular momentum of a group of atoms
@@ -414,6 +418,7 @@ accelerated styles exist.
 * :doc:`wall/lj1043 <fix_wall>` - Lennard-Jones 10--4--3 wall
 * :doc:`wall/lj126 <fix_wall>` - Lennard-Jones 12--6 wall
 * :doc:`wall/lj93 <fix_wall>` - Lennard-Jones 9--3 wall
+* :doc:`wall/lepton <fix_wall>` - Custom Lepton expression wall
 * :doc:`wall/morse <fix_wall>` - Morse potential wall
 * :doc:`wall/piston <fix_wall_piston>` - moving reflective piston wall
 * :doc:`wall/reflect <fix_wall_reflect>` - reflecting wall(s)
@@ -421,6 +426,7 @@ accelerated styles exist.
 * :doc:`wall/region <fix_wall_region>` - use region surface as wall
 * :doc:`wall/region/ees <fix_wall_ees>` - use region surface as wall for ellipsoidal particles
 * :doc:`wall/srd <fix_wall_srd>` - slip/no-slip wall for SRD particles
+* :doc:`wall/table <fix_wall>` - Tabulated potential wall wall
 * :doc:`widom <fix_widom>` - Widom insertions of atoms or molecules
 
 Restrictions

doc/src/fix_adapt.rst

@@ -133,13 +133,15 @@ formulas for the meaning of these parameters:
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`born/coul/long, born/coul/msm <pair_born>`                             | coulombic_cutoff                                 | type global |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
+| :doc:`born/gauss <pair_born_gauss>`                                          | biga0,biga1,r0                                   | type pairs  |
++------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`buck, buck/coul/cut  <pair_buck>`                                      | a,c                                              | type pairs  |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`buck/coul/long, buck/coul/msm <pair_buck>`                             | a,c,coulombic_cutoff                             | type pairs  |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`buck/mdf <pair_mdf>`                                                   | a,c                                              | type pairs  |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
-| :doc:`coul/cut <pair_coul>`                                                  | scale                                            | type pairs  |
+| :doc:`coul/cut, coul/cut/global <pair_coul>`                                 | scale                                            | type pairs  |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`coul/cut/soft <pair_fep_soft>`                                         | lambda                                           | type pairs  |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
@@ -151,10 +153,16 @@ formulas for the meaning of these parameters:
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`coul/long/soft <pair_fep_soft>`                                        | scale, lambda, coulombic_cutoff                  | type pairs  |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
+| :doc:`coul/slater/long <pair_coul_slater>`                                   | scale                                            | type pairs  |
++------------------------------------------------------------------------------+--------------------------------------------------+-------------+
+| :doc:`coul/streitz <pair_coul>`                                              | scale                                            | type pairs  |
++------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`eam, eam/alloy, eam/fs <pair_eam>`                                     | scale                                            | type pairs  |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`gauss <pair_gauss>`                                                    | a                                                | type pairs  |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
+| :doc:`harmonic/cut <pair_harmonic_cut>`                                      | k, cutoff                                        | type pairs  |
++------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`lennard/mdf <pair_mdf>`                                                | A,B                                              | type pairs  |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`lj/class2 <pair_class2>`                                               | epsilon,sigma                                    | type pairs  |
@@ -181,6 +189,8 @@ formulas for the meaning of these parameters:
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`lubricate <pair_lubricate>`                                            | mu                                               | global      |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
+| :doc:`meam <pair_meam>`                                                      | scale                                            | type pairs  |
++------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`mie/cut <pair_mie>`                                                    | epsilon,sigma,gamma_repulsive,gamma_attractive   | type pairs  |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`morse, morse/smooth/linear <pair_morse>`                               | D0,R0,alpha                                      | type pairs  |
@@ -191,7 +201,7 @@ formulas for the meaning of these parameters:
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`nm/cut/coul/cut, nm/cut/coul/long <pair_nm>`                           | E0,R0,m,n,coulombic_cutoff                       | type pairs  |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
-| :doc:`reaxff <pair_reaxff>`                                                  | chi, eta, gamma                                  | type global |
+| :doc:`pace, pace/extrapolation <pair_pace>`                                  | scale                                            | type pairs  |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`snap <pair_snap>`                                                      | scale                                            | type pairs  |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
@@ -203,11 +213,13 @@ formulas for the meaning of these parameters:
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`spin/neel <pair_spin_neel>`                                            | coulombic_cutoff                                 | type global |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
+| :doc:`soft <pair_soft>`                                                      | a                                                | type pairs  |
++------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`table <pair_table>`                                                    | table_cutoff                                     | type pairs  |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 | :doc:`ufm <pair_ufm>`                                                        | epsilon,sigma                                    | type pairs  |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
-| :doc:`soft <pair_soft>`                                                      | a                                                | type pairs  |
+| :doc:`wf/cut <pair_wf_cut>`                                                  | epsilon,sigma,nu,mu                              | type pairs  |
 +------------------------------------------------------------------------------+--------------------------------------------------+-------------+
 
 .. note::

doc/src/fix_alchemy.rst

@@ -0,0 +1,173 @@
+.. index:: fix alchemy
+
+fix alchemy command
+===================
+
+Syntax
+""""""
+
+.. parsed-literal::
+
+   fix ID group-ID alchemy v_name
+
+* ID, group-ID are documented in :doc:`fix <fix>` command
+* alchemy = style name of this fix command
+* v_name = variable with name that determines the :math:`\lambda_R` value
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   fix trans all alchemy v_ramp
+
+Description
+"""""""""""
+
+.. versionadded:: 28Mar2023
+
+This fix command enables an "alchemical transformation" to be performed
+between two systems, whereby one system slowly transforms into the other
+over the course of a molecular dynamics run.  This is useful for
+measuring thermodynamic differences between two different systems.  It
+also allows transformations that are not easily possible with the
+:doc:`pair style hybrid/scaled <pair_hybrid>`, :doc:`fix adapt
+<fix_adapt>` or :doc:`fix adapt/fep <fix_adapt_fep>` commands.
+
+Example inputs are included in the ``examples/PACKAGES/alchemy``
+directory for (a) transforming a pure copper system into a
+copper/aluminum bronze alloy and (b) transforming two water molecules
+in a box of water into a hydronium and a hydroxyl ion.
+
+The two systems must be defined as :doc:`separate replica
+<Howto_replica>` and run in separate partitions of processors using the
+:doc:`-partition <Run_options>` command-line switch.  Exactly two
+partitions must be specified, and each partition must use the same number
+of processors and the same domain decomposition.
+
+Because the forces applied to the atoms are the same mix of the forces
+from each partition and the simulation starts with the same atom
+positions across both partitions, they will generate the same trajectory
+of coordinates for each atom, and the same simulation box size and
+shape.  The latter two conditions are *enforced* by this fix; it
+exchanges coordinates and box information between the replicas.  This is
+not strictly required, but since MD simulations are an example of a
+chaotic system, even the tiniest random difference will eventually grow
+exponentially into an unwanted divergence.
+
+Otherwise, the properties of each atom (type, charge, bond and angle
+partners, etc.), as well as energy and forces between interacting atoms
+(pair, bond, angle styles, etc.) can be different in the two systems.
+
+This can be initialized in the same input script by using commands which
+only apply to one or the other replica.  The example scripts use a
+world-style :doc:`variable <variable>` command along with
+:doc:`if/then/else <if>` commands for this purpose.  The
+:doc:`partition <partition>` command can also be used.
+
+.. code-block:: LAMMPS
+
+   create_box 2 box
+   create_atoms 1 box
+   pair_style eam/alloy
+   pair_coeff * * AlCu.eam.alloy Cu Al
+
+   # replace 5% of copper with aluminum on the second partition only
+
+   variable name world pure alloy
+   if "${name} == alloy" then &
+     "set type 1 type/fraction 2 0.05 6745234"
+
+Both replicas must define an instance of this fix, but with a different
+*v_name* variable.  The named variable must be an equal-style or
+equivalent :doc:`variable <variable>`.  The two variables should be
+defined so that one ramps *down* from 1.0 to 0.0 for the *first* replica
+(*R=0*) and the other ramps *up* from 0.0 to 1.0 for the *second*
+replica (*R=1*).  A simple way is to do this is linearly, which can be
+done using the ramp() function of the :doc:`variable <variable>`
+command.  You could also define a variable which returns a value between
+0.0 and 1.0 as a non-linear function of the timestep.  Here is a linear
+example:
+
+.. code-block:: LAMMPS
+
+   partition yes 1 variable ramp equal ramp(1.0,0.0)
+   partition yes 2 variable ramp equal ramp(0.0,1.0)
+   fix 2 all alchemy v_ramp
+
+.. note::
+
+   For an alchemical transformation, the two variables should sum to
+   exactly 1.0 at any timestep.  LAMMPS does *NOT* check that this is
+   the case.
+
+If you use the ``ramp()`` function to define the two variables, this fix
+can easily be used across successive runs in the same input script by
+ensuring each instance of the :doc:`run <run>` command specifies the
+appropriate *start* or *stop* options.
+
+At each timestep of an MD run, the two instances of this fix evaluate
+their respective variables as a :math:`\lambda_R` factor, where *R* = 0
+or 1 for each replica.  The forces used by each system for the
+propagation of their atoms is set to the sum of the forces for the two
+systems, each scaled by their respective :math:`\lambda_R` factor.  Thus,
+during the MD run, the system will transform incrementally from the
+first system to the second system.
+
+.. note::
+
+   As mentioned above, the coordinates of the atoms and box size/shape
+   must be exactly the same in the two replicas.  Therefore, it is
+   generally not a good idea to initialize the two replicas by reading
+   different data files or creating them individually from scratch.
+   Rather, a single system should be initialized and then desired
+   modifications applied to the system to either replica.  If your
+   input script somehow induces the two systems to become different
+   (e.g. by performing :doc:`atom_modify sort <atom_modify>`
+   differently, or by adding or depositing a different number of atoms),
+   then LAMMPS will detect the mismatch and generate an error.  This is
+   done by ensuring that each step the number and ordering of atoms is
+   identical within each pair of processors in the two replicas.
+
+----------
+
+Restart, fix_modify, output, run start/stop, minimize info
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+No information about this fix is written to :doc:`binary restart files
+<restart>`.  None of the :doc:`fix_modify <fix_modify>` options are
+relevant to this fix.
+
+This fix stores a global scalar (the current value of :math:`\lambda_R`)
+and a global vector of length 3 which contains the potential energy of
+the first partition, the second partition and the combined value,
+respectively. The global scalar is unitless and "intensive", the vector
+is in :doc:`energy units <units>` and "extensive".  These values can be
+used by any command that uses a global value from a fix as input.  See
+the :doc:`output howto <Howto_output>` page for an overview of LAMMPS
+output options.
+
+This fix is not invoked during :doc:`energy minimization <minimize>`.
+
+Restrictions
+""""""""""""
+
+This fix is part of the REPLICA package.  It is only enabled if LAMMPS
+was built with that package.  See the :doc:`Build package
+<Build_package>` page for more info.
+
+There may be only one instance of this fix in use at a time within
+each replica.
+
+
+Related commands
+""""""""""""""""
+
+:doc:`compute pressure/alchemy <compute_pressure_alchemy>` command,
+:doc:`fix adapt <fix_adapt>` command, :doc:`fix adapt/fep <fix_adapt_fep>`
+command, :doc:`pair_style hybrid/scaled <pair_hybrid>` command.
+
+Default
+"""""""
+
+none

doc/src/fix_ave_correlate_long.rst

@@ -15,7 +15,7 @@ Syntax
 * Nevery = use input values every this many time steps
 * Nfreq = save state of the time correlation functions every this many time steps
 * one or more input values can be listed
-* value = c_ID, c_ID[N], f_ID, f_ID[N], v_name
+* value = c_ID, c_ID[N], f_ID, f_ID[N], v_name, v_name[I]
 
   .. parsed-literal::
 
@@ -24,6 +24,7 @@ Syntax
        f_ID = global scalar calculated by a fix with ID
        f_ID[I] = Ith component of global vector calculated by a fix with ID
        v_name = global value calculated by an equal-style variable with name
+       v_name[I] = Ith component of global vector calculated by a vector-style variable with name
 
 * zero or more keyword/arg pairs may be appended
 * keyword = *type* or *start* or *file* or *overwrite* or *title1* or *title2* or *ncorr* or *nlen* or *ncount*

doc/src/fix_balance.rst

@@ -43,6 +43,7 @@ Syntax
              name = name of the atom-style variable
            *store* name = store weight in custom atom property defined by :doc:`fix property/atom <fix_property_atom>` command
              name = atom property name (without d\_ prefix)
+       *sort* arg = *no* or *yes*
        *out* arg = filename
          filename = write each processor's subdomain to a file, at each re-balancing
 
@@ -308,6 +309,14 @@ in that sub-box.
 
 ----------
 
+The *sort* keyword determines whether the communication of per-atom
+data to other processors during load-balancing will be random or
+deterministic.  Random is generally faster; deterministic will ensure
+the new ordering of atoms on each processor is the same each time the
+same simulation is run.  This can be useful for debugging purposes.
+Since the fix balance command is performed during timestepping, the
+default is *no* so that sorting is not performed.
+
 The *out* keyword writes text to the specified *filename* with the
 results of each re-balancing operation.  The file contains the bounds
 of the subdomain for each processor after the balancing operation
@@ -415,4 +424,4 @@ Related commands
 Default
 """""""
 
-none
+The default setting is sort = no.

doc/src/fix_deposit.rst

@@ -220,11 +220,11 @@ ignored if the *global* or *local* keywords are used, since those
 options choose a z-coordinate for insertion independently.
 
 The vx, vy, and vz components of velocity for the inserted particle
-are set using the values specified for the *vx*, *vy*, and *vz*
+are set by sampling a uniform distribution between the bounds set by
-keywords.  Note that normally, new particles should be a assigned a
+the values specified for the *vx*, *vy*, and *vz* keywords. Note that
-negative vertical velocity so that they move towards the surface.  For
+normally, new particles should be a assigned a negative vertical
-molecules, the same velocity is given to every particle (no rotation
+velocity so that they move towards the surface.  For molecules, the
-or bond vibration).
+same velocity is given to every particle (no rotation or bond vibration).
 
 If the *target* option is used, the velocity vector of the inserted
 particle is changed so that it points from the insertion position

doc/src/fix_efield.rst

@@ -1,17 +1,21 @@
 .. index:: fix efield
+.. index:: fix efield/tip4p
 
 fix efield command
 ==================
 
+fix efield/tip4p command
+========================
+
 Syntax
 """"""
 
 .. parsed-literal::
 
-   fix ID group-ID efield ex ey ez keyword value ...
+   fix ID group-ID style ex ey ez keyword value ...
 
 * ID, group-ID are documented in :doc:`fix <fix>` command
-* efield = style name of this fix command
+* style = *efield* or *efield/tip4p*
 * ex,ey,ez = E-field component values (electric field units)
 * any of ex,ey,ez can be a variable (see below)
 * zero or more keyword/value pairs may be appended to args
@@ -31,27 +35,36 @@ Examples
 
    fix kick external-field efield 1.0 0.0 0.0
    fix kick external-field efield 0.0 0.0 v_oscillate
+   fix kick external-field efield/tip4p 1.0 0.0 0.0
 
 Description
 """""""""""
 
-Add a force F = qE to each charged atom in the group due to an
+Add a force :math:`\vec{F} = q\vec{E}` to each charged atom in the group due to an
 external electric field being applied to the system.  If the system
 contains point-dipoles, also add a torque on the dipoles due to the
 external electric field.
 
-For charges, any of the 3 quantities defining the E-field components
+.. versionadded:: 28Mar2023
-can be specified as an equal-style or atom-style
+
-:doc:`variable <variable>`, namely *ex*, *ey*, *ez*\ .  If the value is a
+When the *efield/tip4p* style is used, the E-field will be applied to
-variable, it should be specified as v_name, where name is the variable
+the position of the virtual charge site M of a TIP4P molecule instead of
-name.  In this case, the variable will be evaluated each timestep, and
+the oxygen position as it is defined by a corresponding :doc:`TIP4P pair
-its value used to determine the E-field component.
+style <pair_lj_cut_tip4p>`.  The forces on the M site due to the
+external field are projected on the oxygen and hydrogen atoms of the
+TIP4P molecules.
+
+For charges, any of the 3 quantities defining the E-field components can
+be specified as an equal-style or atom-style :doc:`variable <variable>`,
+namely *ex*, *ey*, *ez*\ .  If the value is a variable, it should be
+specified as v_name, where name is the variable name.  In this case, the
+variable will be evaluated each timestep, and its value used to
+determine the E-field component.
 
 For point-dipoles, equal-style variables can be used, but atom-style
 variables are not currently supported, since they imply a spatial
 gradient in the electric field which means additional terms with
-gradients of the field are required for the force and torque on
+gradients of the field are required for the force and torque on dipoles.
-dipoles.
 
 Equal-style variables can specify formulas with various mathematical
 functions, and include :doc:`thermo_style <thermo_style>` command
@@ -81,10 +94,18 @@ self-consistent minimization problem (see below).
 The *energy* keyword is not allowed if the added field is a constant
 vector (ex,ey,ez), with all components defined as numeric constants
 and not as variables.  This is because LAMMPS can compute the energy
-for each charged particle directly as E = -x dot qE = -q (x\*ex + y\*ey
+for each charged particle directly as
-+ z\*ez), so that -Grad(E) = F.  Similarly for point-dipole particles
+
-the energy can be computed as E = -mu dot E = -(mux\*ex + muy\*ey +
+.. math::
-muz\*ez).
+
+   U_{efield} = -\vec{x} \cdot q\vec{E} = -q (x\cdot E_x + y\cdot E_y + z\cdot Ez),
+
+so that :math:`-\nabla U_{efield} = \vec{F}`.  Similarly for point-dipole particles
+the energy can be computed as
+
+.. math::
+
+   U_{efield} = -\vec{\mu} \cdot \vec{E} = -\mu_x\cdot E_x + \mu_y\cdot E_y + \mu_z\cdot E_z
 
 The *energy* keyword is optional if the added force is defined with
 one or more variables, and if you are performing dynamics via the
@@ -120,29 +141,28 @@ Restart, fix_modify, output, run start/stop, minimize info
 No information about this fix is written to :doc:`binary restart files
 <restart>`.
 
-The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-this fix to add the potential energy inferred by the added force due
+fix to add the potential energy inferred by the added force due to the
-to the electric field to the global potential energy of the system as
+electric field to the global potential energy of the system as part of
-part of :doc:`thermodynamic output <thermo_style>`.  The default
+:doc:`thermodynamic output <thermo_style>`.  The default setting for
-setting for this fix is :doc:`fix_modify energy no <fix_modify>`.
+this fix is :doc:`fix_modify energy no <fix_modify>`.  Note that this
-Note that this energy is a fictitious quantity but is needed so that
+energy is a fictitious quantity but is needed so that the :doc:`minimize
-the :doc:`minimize <minimize>` command can include the forces added by
+<minimize>` command can include the forces added by this fix in a
-this fix in a consistent manner.  I.e. there is a decrease in
+consistent manner.  I.e. there is a decrease in potential energy when
-potential energy when atoms move in the direction of the added force
+atoms move in the direction of the added force due to the electric
-due to the electric field.
+field.
 
-The :doc:`fix_modify <fix_modify>` *virial* option is supported by
+The :doc:`fix_modify <fix_modify>` *virial* option is supported by this
-this fix to add the contribution due to the added forces on atoms to
+fix to add the contribution due to the added forces on atoms to both the
-both the global pressure and per-atom stress of the system via the
+global pressure and per-atom stress of the system via the :doc:`compute
-:doc:`compute pressure <compute_pressure>` and :doc:`compute
+pressure <compute_pressure>` and :doc:`compute stress/atom
-stress/atom <compute_stress_atom>` commands.  The former can be
+<compute_stress_atom>` commands.  The former can be accessed by
-accessed by :doc:`thermodynamic output <thermo_style>`.  The default
+:doc:`thermodynamic output <thermo_style>`.  The default setting for
-setting for this fix is :doc:`fix_modify virial no <fix_modify>`.
+this fix is :doc:`fix_modify virial no <fix_modify>`.
 
 The :doc:`fix_modify <fix_modify>` *respa* option is supported by this
-fix. This allows to set at which level of the :doc:`r-RESPA
+fix. This allows to set at which level of the :doc:`r-RESPA <run_style>`
-<run_style>` integrator the fix adding its forces. Default is the
+integrator the fix adding its forces. Default is the outermost level.
-outermost level.
 
 This fix computes a global scalar and a global 3-vector of forces,
 which can be accessed by various :doc:`output commands
@@ -169,7 +189,11 @@ the iteration count during the minimization.
 Restrictions
 """"""""""""
 
-None
+Fix style *efield/tip4p* 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.
+
+Fix style *efield/tip4p* can only be used with tip4p pair styles.
 
 Related commands
 """"""""""""""""

doc/src/fix_heat_flow.rst

@@ -0,0 +1,70 @@
+.. index:: fix heat/flow
+
+fix heat/flow command
+==========================
+
+Syntax
+""""""
+
+.. parsed-literal::
+
+   fix ID group-ID heat/flow style values ...
+
+* ID, group-ID are documented in :doc:`fix <fix>` command
+* heat/flow = style name of this fix command
+* one style with corresponding value(s) needs to be listed
+
+  .. parsed-literal::
+
+     style = *constant* or *type*
+       *constant* = cp
+         cp = value of specifc heat (energy/(mass * temperature) units)
+       *type* = cp1 ... cpN
+         cpN = value of specifc heat for type N (energy/(mass * temperature) units)
+
+*
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   fix 1 all heat/flow constant 1.0
+   fix 1 all heat/flow type 1.0 0.5
+
+Description
+"""""""""""
+
+Perform plain time integration to update temperature for atoms in the
+group each timestep. The specific heat of atoms can be defined using either
+the *constant* or *type* keywords. For style *constant*, the specific heat
+is a constant value *cp* for all atoms. For style *type*, *N* different values
+of the specific heat are defined, one for each of the *N* types of atoms.
+
+----------
+
+Restart, fix_modify, output, run start/stop, minimize info
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+No information about this fix is written to :doc:`binary restart files <restart>`.
+None of the :doc:`fix_modify <fix_modify>` options are relevant to this fix.
+No global or per-atom quantities are stored by this fix for access by various
+:doc:`output commands <Howto_output>`. No parameter of this fix can be used
+with the *start/stop* keywords of the :doc:`run <run>` command.  This fix is
+not invoked during :doc:`energy minimization <minimize>`.
+
+Restrictions
+""""""""""""
+
+This fix requires that atoms store temperature and heat flow
+as defined by the :doc:`fix property/atom <fix_property_atom>` command.
+
+Related commands
+""""""""""""""""
+
+:doc:`pair granular <pair_granular>`, :doc:`fix property/atom <fix_property_atom>`
+
+Default
+"""""""
+
+none

doc/src/fix_mdi_qm.rst

@@ -8,12 +8,12 @@ Syntax
 
 .. parsed-literal::
 
-   fix ID group-ID mdi/qm keyword
+   fix ID group-ID mdi/qm keyword value(s) keyword value(s) ...
 
 * ID, group-ID are documented in :doc:`fix <fix>` command
 * mdi/qm = style name of this fix command
 * zero or more keyword/value pairs may be appended
-* keyword = *virial* or *add* or *every* or *connect* or *elements*
+* keyword = *virial* or *add* or *every* or *connect* or *elements* or *mc*
 
   .. parsed-literal::
 
@@ -29,7 +29,9 @@ Syntax
          yes = perform a one-time connection to the MDI engine code
          no = do not perform the connection operation
        *elements* args = N_1 N_2 ... N_ntypes
-         N_1,N_2,...N_ntypes = atomic number for each of ntypes LAMMPS atom types
+         N_1,N_2,...N_ntypes = chemical symbol for each of ntypes LAMMPS atom types
+       *mc* args = mcfixID
+         mcfixID = ID of a Monte Carlo fix designed to work with this fix
 
 Examples
 """"""""
@@ -38,7 +40,7 @@ Examples
 
    fix 1 all mdi/qm
    fix 1 all mdi/qm virial yes
-   fix 1 all mdi/qm add no every 100 elements 13 29
+   fix 1 all mdi/qm add no every 100 elements C C H O
 
 Description
 """""""""""
@@ -57,12 +59,27 @@ The server code must support use of the `MDI Library
 <https://molssi-mdi.github.io/MDI_Library/html/index.html>`_ as
 explained below.
 
+Typically, to use this fix, the input script should not define any
+other classical force field components, e.g. a pair style, bond style,
+etc.
+
 These are example use cases for this fix, discussed further below:
 
 * perform an ab initio MD (AIMD) simulation with quantum forces
 * perform an energy minimization with quantum forces
 * perform a nudged elastic band (NEB) calculation with quantum forces
-* perform a QM calculation for a series of independent systems which LAMMPS reads or generates
+* perform a QM calculation for a series of independent systems which
+  LAMMPS reads or generates once
+* run a classical MD simulation and calculate QM energy/forces once
+  every N steps on the current configuration
+
+More generally any command which calculates per-atom forces can instead
+use quantum forces by defining this fix.  Examples are the Monte Carlo
+commands :doc:`fix gcmc <fix_gcmc>` and :doc:`fix atom/swap
+<fix_atom_swap>`, as well as the :doc:`compute born/matrix
+<compute_born_matrix>` command.  The only requirement is that internally
+the command invokes the post_force() method of fixes such as this one,
+which will trigger the quantum calculation.
 
 The code coupling performed by this command is done via the `MDI
 Library <https://molssi-mdi.github.io/MDI_Library/html/index.html>`_.
@@ -72,27 +89,32 @@ for MDI.  See the :doc:`Howto mdi <Howto_mdi>` page for more
 information about how LAMMPS can operate as either an MDI driver or
 engine.
 
-The examples/mdi directory contains input scripts using this fix in
+The ``examples/mdi`` directory contains input scripts using this fix in
 the various use cases discussed below.  In each case, two instances of
-LAMMPS are used, once as an MDI driver, once as an MDI engine
+LAMMPS are used, once as an MDI driver, once as an MDI engine (surrogate
-(surrogate for a QM code).  The examples/mdi/README file explains how
+for a QM code).  The ``examples/mdi/README`` file explains how to launch
-to launch two codes so that they communicate via the MDI library using
+two codes so that they communicate via the MDI library using either MPI
-either MPI or sockets.  Any QM code that supports MDI could be used in
+or sockets.  Any QM code that supports MDI could be used in place of
-place of LAMMPS acting as a QM surrogate.  See the :doc:`Howto mdi
+LAMMPS acting as a QM surrogate.  See the :doc:`Howto mdi <Howto_mdi>`
-<Howto_mdi>` page for a current list (March 2022) of such QM codes.
+page for a current list (March 2022) of such QM codes.  The
+``examples/QUANTUM`` directory has examples for coupling LAMMPS to 3 QM
+codes either via this fix or the :doc:`fix mdi/qmmm <fix_mdi_qmmm>`
+command.
 
-Note that an engine code can support MDI in either or both of two
+Note that an engine code can support MDI in either or both of two modes.
-modes.  It can be used as a stand-alone code, launched at the same
+It can be used as a stand-alone code, launched at the same time as
-time as LAMMPS.  Or it can be used as a plugin library, which LAMMPS
+LAMMPS.  Or it can be used as a plugin library, which LAMMPS loads.  See
-loads.  See the :doc:`mdi plugin <mdi>` command for how to trigger
+the :doc:`mdi plugin <mdi>` command for how to trigger LAMMPS to load a
-LAMMPS to load a plugin library.  The examples/mdi/README file
+plugin library.  The ``examples/mdi/README`` file and
-explains how to launch the two codes in either mode.
+``examples/QUANTUM/QM-code/README`` files explain how to launch the two
+codes in either mode.
 
 ----------
 
-The *virial* keyword setting of yes or no determines whether
+The *virial* keyword setting of yes or no determines whether LAMMPS
-LAMMPS will request the QM code to also compute and return
+will request the QM code to also compute and return the QM
-a 6-element symmetric virial tensor for the system.
+contribution to a stress tensor for the system which LAMMPS will
+convert to a 6-element symmetric virial tensor.
 
 The *add* keyword setting of *yes* or *no* determines whether the
 energy and forces and virial returned by the QM code will be added to
@@ -109,25 +131,27 @@ commands.  See details below.
 The *every* keyword determines how often the QM code will be invoked
 during a dynamics run with the current LAMMPS simulation box and
 configuration of atoms.  The QM code will be called once every
-*Nevery* timesteps.
+*Nevery* timesteps.  By default *Nevery* = 1.
 
 The *connect* keyword determines whether this fix performs a one-time
-connection to the QM code.  The default is *yes*.  The only time a
+connection to the QM code.  The default is *yes*.  The only time a *no*
-*no* is needed is if this command is used multiple times in an input
+is needed is if this command is used multiple times in an input script
-script.  E.g. if it used inside a loop which also uses the :doc:`clear
+and the MDI coupling is between two stand-alone codes (not plugin mode).
-<clear>` command to destroy the system (including any defined fixes).
+E.g. if it used inside a loop which also uses the :doc:`clear <clear>`
-See the examples/mdi/in.series.driver script as an example of this,
+command to destroy the system (including this fix).  See the
-where LAMMPS is using the QM code to compute energy and forces for a
+``examples/mdi/in.series.driver`` script as an example of this, where
-series of system configurations.  In this use case *connect no*
+LAMMPS is using the QM code to compute energy and forces for a series of
-is used along with the :doc:`mdi connect and exit <mdi>` command
+system configurations.  In this use case *connect no* is used along with
-to one-time initiate/terminate the connection outside the loop.
+the :doc:`mdi connect and exit <mdi>` command to one-time
+initiate/terminate the connection outside the loop.
 
 The *elements* keyword allows specification of what element each
-LAMMPS atom type corresponds to.  This is specified by the atomic
+LAMMPS atom type corresponds to.  This is specified by the chemical
-number of the element, e.g. 13 for Al.  An atomic number must be
+symbol of the element, e.g. C or Al or Si.  A symbol must be specified
-specified for each of the ntypes LAMMPS atom types.  Ntypes is
+for each of the ntypes LAMMPS atom types.  Multiple LAMMPS types can
-typically specified via the create_box command or in the data file
+represent the same element.  Ntypes is typically specified via the
-read by the read_data command.
+:doc:`create_box <create_box>` command or in the data file read by the
+:doc:`read_data <read_data>` command.
 
 If this keyword is specified, then this fix will send the MDI
 ">ELEMENTS" command to the engine, to ensure the two codes are
@@ -136,10 +160,18 @@ not specified, then this fix will send the MDI >TYPES command to the
 engine.  This is fine if both the LAMMPS driver and the MDI engine are
 initialized so that the atom type values are consistent in both codes.
 
+The *mc* keyword enables this fix to be used with a Monte Carlo (MC)
+fix to calculate before/after quantum energies as part of the MC
+accept/reject criterion.  The :doc:`fix gcmc <fix_gcmc>` and :doc:`fix
+atom/swap <fix_atom_swap>` commands can be used in this manner.
+Specify the ID of the MC fix following the *mc* keyword.  This allows
+the two fixes to coordinate when MC events are being calculated versus
+MD timesteps between the MC events.
+
 ----------
 
-The following 3 example use cases are illustrated in the examples/mdi
+The following 3 example use cases are illustrated in the
-directory.  See its README file for more details.
+``examples/mdi`` directory.  See its README file for more details.
 
 (1) To run an ab initio MD (AIMD) dynamics simulation, or an energy
 minimization with QM forces, or a multi-replica NEB calculation, use
@@ -252,12 +284,6 @@ This command is part of the MDI package.  It is only enabled if
 LAMMPS was built with that package.  See the :doc:`Build package
 <Build_package>` page for more info.
 
-The QM code does not currently compute and return per-atom energy or
-per-atom virial contributions.  So they will not show up as part of
-the calculations performed by the :doc:`compute pe/atom
-<compute_pe_atom>` or :doc:`compute stress/atom <compute_stress_atom>`
-commands.
-
 To use LAMMPS as an MDI driver in conjunction with other MDI-enabled
 codes (MD or QM codes), the :doc:`units <units>` command should be
 used to specify *real* or *metal* units.  This will ensure the correct
@@ -265,12 +291,15 @@ unit conversions between LAMMPS and MDI units.  The other code will
 also perform similar unit conversions into its preferred units.
 
 LAMMPS can also be used as an MDI driver in other unit choices it
-supports, e.g. *lj*, but then no unit conversion is performed.
+supports, e.g. *lj*, but then no unit conversion to MDI units is
+performed.
 
 Related commands
 """"""""""""""""
 
-:doc:`mdi plugin <mdi>`, :doc:`mdi engine <mdi>`
+:doc:`mdi plugin <mdi>`,
+:doc:`mdi engine <mdi>`,
+:doc:`fix mdi/qmmm <fix_mdi_qmmm>`
 
 Default
 """""""

doc/src/fix_mdi_qmmm.rst

@@ -0,0 +1,274 @@
+.. index:: fix mdi/qmmm
+
+fix mdi/qmmm command
+====================
+
+Syntax
+""""""
+
+.. parsed-literal::
+
+   fix ID group-ID mdi/qmmm mode keyword value(s) keyword value(s) ...
+
+* ID, group-ID are documented in :doc:`fix <fix>` command
+* mdi/qmmm = style name of this fix command
+* mode = *direct* or *potential*
+* zero or more keyword/value pairs may be appended
+* keyword = *virial* or *add* or *every* or *connect* or *elements*
+
+  .. parsed-literal::
+
+       *virial* args = *yes* or *no*
+         yes = request virial tensor from server code
+         no = do not request virial tensor from server code
+       *connect* args = *yes* or *no*
+         yes = perform a one-time connection to the MDI engine code
+         no = do not perform the connection operation
+       *elements* args = N_1 N_2 ... N_ntypes
+         N_1,N_2,...N_ntypes = chemical symbol for each of ntypes LAMMPS atom types
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   fix 1 all mdi/qmmm direct
+   fix 1 all mdi/qmmm potential virial yes
+   fix 1 all mdi/qmmm potential virial yes elements 13 29
+
+Description
+"""""""""""
+
+.. versionadded:: 28Mar2023
+
+This command enables LAMMPS to act as a client with another server code
+to perform a coupled QM/MM (quantum-mechanics/molecular-mechanics)
+simulation.  LAMMPS will perform classical MD (molecular mechanics
+or MM) for the (typically larger) MM portion of the system.  A quantum
+mechanics code will calculate quantum energy and forces for the QM
+portion of the system.  The two codes work together to calculate the
+energy and forces due to the cross interactions between QM and MM atoms.
+The QM server code must support use of the `MDI Library
+<https://molssi-mdi.github.io/MDI_Library/html/index.html>`_ as
+explained below.
+
+The partitioning of the system between QM and MM atoms is as follows.
+Atoms in the specified group are QM atoms; the remaining atoms are MM
+atoms.  The input script should thus define this partitioning.
+See additional information below about other requirements for an input
+script to use this fix and perform a QM/MM simulation.
+
+The code coupling performed by this command is done via the `MDI
+Library <https://molssi-mdi.github.io/MDI_Library/html/index.html>`_.
+LAMMPS runs as an MDI driver (client), and sends MDI commands to an
+external MDI engine code (server), in this case a QM code which has
+support for MDI.  See the :doc:`Howto mdi <Howto_mdi>` page for more
+information about how LAMMPS can operate as either an MDI driver or
+engine.
+
+The ``examples/QUANTUM`` directory has sub-directories with example
+input scripts using this fix in tandem with different QM codes.  The
+README files in the sub-directories explain how to download and build
+the various QM codes.  They also explain how to launch LAMMPS and the QM
+code so that they communicate via the MDI library using either MPI or
+sockets.  Any QM code that supports MDI could be used in addition to
+those discussed in the sub-directories.  See the :doc:`Howto mdi
+<Howto_mdi>` page for a current list (March 2022) of such QM codes.
+
+Note that an engine code can support MDI in either or both of two modes.
+It can be used as a stand-alone code, launched at the same time as
+LAMMPS.  Or it can be used as a plugin library, which LAMMPS loads.  See
+the :doc:`mdi plugin <mdi>` command for how to trigger LAMMPS to load a
+plugin library.  The ``examples/QUANTUM`` sub-directory README files
+explains how to launch the two codes in either mode.
+
+----------
+
+The *mode* setting determines which QM/MM coupling algorithm is used.
+LAMMPS currently supports *direct* and *potential* algorithms, based
+on the *mode* setting.  Both algorithms should give reasonably
+accurate results, but some QM codes support only one of the two modes.
+E.g. in the ``examples/QUANTUM`` directory, PySCF supports only *direct*,
+NWChem supports only *potential*, and LATTE currently supports
+neither, so it cannot be used for QM/MM simulations using this fix.
+
+The *direct* option passes the coordinates and charges of each MM atom
+to the quantum code, in addition to the coordinates of each QM atom.
+The quantum code returns forces on each QM atom as well as forces on
+each MM atom.  The latter is effectively the force on MM atoms due to
+the QM atoms.
+
+The input script for performing a *direct* mode QM/MM simulation should
+do the following:
+
+* delete all bonds (angles, dihedrals, etc) between QM atoms
+* set the charge on each QM atom to zero
+* define no bonds (angles, dihedrals, etc) which involve both QM and MM atoms
+* define a force field (pair, bonds, angles, optional kspace) for the entire system
+
+The first two bullet can be performed using the :doc:`delete_bonds
+<delete_bonds>` and :doc:`set <set>` commands.
+
+The third bullet is required to have a consistent model, but is not
+checked by LAMMPS.
+
+The fourth bullet implies that non-bonded non-Coulombic interactions
+(e.g. van der Waals) between QM/QM and QM/MM pairs of atoms are
+computed by LAMMPS.
+
+See the ``examples/QUANTUM/PySCF/in.*`` files for examples of input
+scripts for QM/MM simulations using the *direct* mode.
+
+The *potential* option passes the coordinates of each QM atom and a
+Coulomb potential for each QM atom to the quantum code.  The latter is
+calculated by performing a Coulombics-only calculation for the entire
+system, subtracting all QM/QM pairwise Coulombic terms, and dividing
+the Coulomb energy on each QM atom by the charge of the QM atom.  The
+potential value represents the Coulombic influence of all the MM atoms
+on each QM atom.
+
+The quantum code returns forces and charge on each QM atom.  The new
+charges on the QM atom are used to re-calculate the MM force field,
+resulting in altered forces on the MM atoms.
+
+The input script for performing a *potential* mode QM/MM simulation
+should do the following:
+
+* delete all bonds (angles, dihedrals, etc) between QM atoms
+* define a hybrid pair style which includes a Coulomb-only pair sub-style
+* define no bonds (angles, dihedrals, etc) which involve both QM and MM atoms
+* define a force field (pair, bonds, angles, optional kspace) for the entire system
+
+The first operation can be performed using the :doc:`delete_bonds
+<delete_bonds>` command.  See the ``examples/QUANTUM/NWChem/in.*`` files
+for examples of how to do this.
+
+The second operation is necessary so that this fix can calculate the
+Coulomb potential for the QM atoms.
+
+The third bullet is required to have a consistent model, but is not
+checked by LAMMPS.
+
+The fourth bullet implies that non-bonded non-Coulombic interactions
+(e.g. van der Waals) between QM/QM and QM/MM pairs of atoms are computed
+by LAMMPS.  However, some QM codes do not want the MM code (LAMMPS) to
+compute QM/QM van der Waals interactions.  NWChem is an example.  In
+this case, the coefficients for those interactions need to be turned
+off, which typically requires the atom types for the QM atoms be
+different than those for the MM atoms.
+
+See the ``examples/QUANTUM/NWChem/in.*`` files for examples of input
+scripts for QM/MM simulations using the *potential* mode.  Those scripts
+also illustrate how to turn off QM/QM van der Waals interactions.
+
+----------
+
+The *virial* keyword setting of yes or no determines whether LAMMPS
+will request the QM code to also compute and return the QM
+contribution to a stress tensor for the system which LAMMPS will
+convert to a 6-element symmetric virial tensor.
+
+The *connect* keyword determines whether this fix performs a one-time
+connection to the QM code.  The default is *yes*.  The only time a
+*no* is needed is if this command is used multiple times in an input
+script.  E.g. if it used inside a loop which also uses the :doc:`clear
+<clear>` command to destroy the system (including this fix).  As
+example would be a script which loop over a series of independent QM/MM
+simulations, e.g. each with their own data file.  In this use case
+*connect no* could be used along with the :doc:`mdi connect and exit
+<mdi>` command to one-time initiate/terminate the connection outside
+the loop.
+
+The *elements* keyword allows specification of what element each
+LAMMPS atom type corresponds to.  This is specified by the chemical
+symbol of the element, e.g. C or Al or Si.  A symbol must be specified
+for each of the ntypes LAMMPS atom types.  Multiple LAMMPS types can
+represent the same element.  Ntypes is typically specified via the
+:doc:`create_box <create_box>` command or in the data file read by the
+:doc:`read_data <read_data>` command.
+
+If this keyword is specified, then this fix will send the MDI
+">ELEMENTS" command to the engine, to insure the two codes are
+consistent in their definition of atomic species.  If this keyword is
+not specified, then this fix will send the MDI >TYPES command to the
+engine.  This is fine if both the LAMMPS driver and the MDI engine are
+initialized so that the atom type values are consistent in both codes.
+
+----------
+
+Restart, fix_modify, output, run start/stop, minimize info
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+No information about this fix is written to :doc:`binary restart files
+<restart>`.
+
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the potential energy computed by the QM code to the
+global potential energy of the system as part of :doc:`thermodynamic
+output <thermo_style>`.  The default setting for this fix is
+:doc:`fix_modify energy yes <fix_modify>`.
+
+The :doc:`fix_modify <fix_modify>` *virial* option is supported by
+this fix to add the contribution computed by the QM code to the global
+pressure of the system as part of :doc:`thermodynamic output
+<thermo_style>`.  The default setting for this fix is :doc:`fix_modify
+virial yes <fix_modify>`.
+
+This fix computes a global scalar which can be accessed by various
+:doc:`output commands <Howto_output>`.  The scalar is the energy
+returned by the QM code.  The scalar value calculated by this fix is
+"extensive".
+
+This fix also computes a global vector with of length 6 which contains
+the symmetric virial tensor values returned by the QM code.  It can
+likewise be accessed by various :doc:`output commands <Howto_output>`.
+
+The ordering of values in the symmetric virial tensor is as follows:
+vxx, vyy, vzz, vxy, vxz, vyz.  The values will be in pressure
+:doc:`units <units>`.
+
+This fix also computes a peratom array with 3 columns which contains
+the peratom forces returned by the QM code.  It can likewise be
+accessed by various :doc:`output commands <Howto_output>`.  Note that
+for *direct* mode this will be quantum forces on both QM and MM atoms.
+For *potential* mode it will only be quantum forces on QM atoms; the
+forces for MM atoms will be zero.
+
+No parameter of this fix can be used with the *start/stop* keywords of
+the :doc:`run <run>` command.
+
+The forces computed by the QM code are used during an energy
+minimization, invoked by the :doc:`minimize <minimize>` command.
+
+.. note::
+
+   If you want the potential energy associated with the QM forces to
+   be included in the total potential energy of the system (the
+   quantity being minimized), you MUST not disable the
+   :doc:`fix_modify <fix_modify>` *energy* option for this fix.
+
+
+Restrictions
+""""""""""""
+
+This command is part of the MDI package.  It is only enabled if
+LAMMPS was built with that package.  See the :doc:`Build package
+<Build_package>` page for more info.
+
+To use LAMMPS as an MDI driver in conjunction with other MDI-enabled
+codes (MD or QM codes), the :doc:`units <units>` command should be
+used to specify *real* or *metal* units.  This will ensure the correct
+unit conversions between LAMMPS and MDI units.  The other code will
+also perform similar unit conversions into its preferred units.
+
+Related commands
+""""""""""""""""
+
+:doc:`mdi plugin <mdi>`,
+:doc:`mdi engine <mdi>`,
+:doc:`fix mdi/qm <fix_mdi_qm>`
+
+Default
+"""""""
+
+The default for the optional keywords are virial = no and connect = yes.

doc/src/fix_mscg.rst

@@ -40,15 +40,15 @@ Description
 
 This fix applies the Multi-Scale Coarse-Graining (MSCG) method to
 snapshots from a dump file to generate potentials for coarse-grained
-simulations from all-atom simulations, using a force-matching
+simulations from all-atom simulations, using a force-matching technique
-technique (:ref:`Izvekov <Izvekov>`, :ref:`Noid <Noid>`).
+(:ref:`Izvekov <Izvekov>`, :ref:`Noid <Noid>`).
 
-It makes use of the MS-CG library, written and maintained by Greg
+It makes use of the MS-CG library, written and maintained by Greg Voth's
-Voth's group at the University of Chicago, which is freely available
+group at the University of Chicago, which is freely available on their
-on their `MS-CG GitHub site <https://github.com/uchicago-voth/MSCG-release>`_.  See instructions
+`MS-CG GitHub site <https://github.com/uchicago-voth/MSCG-release>`_.
-on obtaining and installing the MS-CG library in the src/MSCG/README
+See instructions on obtaining and installing the MS-CG library in the
-file, which must be done before you build LAMMPS with this fix command
+src/MSCG/README file, which must be done before you build LAMMPS with
-and use the command in a LAMMPS input script.
+this fix command and use the command in a LAMMPS input script.
 
 An example script using this fix is provided the examples/mscg
 directory.
@@ -65,15 +65,18 @@ simulations is as follows:
 6. Check the results of the force matching.
 7. Run coarse-grained simulations using the new coarse-grained potentials.
 
-This fix can perform the range finding and force matching steps 4 and
+This fix can perform the range finding and force matching steps 4 and 5
-5 of the above workflow when used in conjunction with the
+of the above workflow when used in conjunction with the :doc:`rerun
-:doc:`rerun <rerun>` command.  It does not perform steps 1-3 and 6-7.
+<rerun>` command.  It does not perform steps 1-3 and 6-7.
 
-Step 2 can be performed using a Python script (what is the name?)
+Step 2 can be performed using a Python script (cgmap), which defines the
-provided with the MS-CG library which defines the coarse-grained model
+coarse-grained model and converts a standard LAMMPS dump file for an
-and converts a standard LAMMPS dump file for an all-atom simulation
+all-atom simulation (step 1) into a LAMMPS dump file which has the
-(step 1) into a LAMMPS dump file which has the positions of and forces
+positions of and forces on the coarse-grained beads.  To use cgmap the
-on the coarse-grained beads.
+following repositories need to be downloaded and installed.
+
+#. The custom lammpsdata branch of mdtraj from https://github.com/hockyg/mdtraj/tree/lammpsdata
+#. The master branch of cgmap from https://github.com/uchicago-voth/cgmap
 
 In step 3, an input file named "control.in" is needed by the MS-CG
 library which sets parameters for the range finding and force matching
@@ -83,12 +86,12 @@ info on this file.
 
 When this fix is used to perform steps 4 and 5, the MS-CG library also
 produces additional output files.  The range finder functionality
-(step 4) outputs files defining pair and bonded interaction ranges.
+(step 4) outputs files defining pair and bonded interaction ranges.  The
-The force matching functionality (step 5) outputs tabulated force
+force matching functionality (step 5) outputs tabulated force files for
-files for every interaction in the system. Other diagnostic files can
+every interaction in the system. Other diagnostic files can also be
-also be output depending on the parameters in the MS-CG library input
+output depending on the parameters in the MS-CG library input script.
-script.  Again, see the documentation provided with the MS-CG library
+Again, see the documentation provided with the MS-CG library for more
-for more info.
+info.
 
 ----------
 
@@ -97,8 +100,8 @@ be invoked. If *on*, the step 4 range finder functionality is invoked.
 *off*, the step 5 force matching functionality is invoked.
 
 If the *name* keyword is used, string names are defined to associate
-with the integer atom types in LAMMPS.  *Ntype* names must be
+with the integer atom types in LAMMPS.  *Ntype* names must be provided,
-provided, one for each atom type (1-Ntype).
+one for each atom type (1-Ntype).
 
 The *max* keyword specifies the maximum number of bonds, angles, and
 dihedrals a bead can have in the coarse-grained model.
@@ -107,16 +110,13 @@ Restrictions
 """"""""""""
 
 This fix is part of the MSCG package. It is only enabled if LAMMPS was
-built with that package.  See the :doc:`Build package <Build_package>`
+built with that package.  Building the MSCG package also requires
-doc page for more info.
+external libraries. See the :doc:`Build_package` and :doc:`Build_extras`
+pages for more info.
 
-The MS-CG library uses C++11, which may not be supported by older
+Currently, the MS-CG library is not set up to run in parallel with MPI,
-compilers. The MS-CG library also has some additional numeric library
+so this fix can only be used in a serial LAMMPS build and run on a
-dependencies, which are described in its documentation.
+single processor.
-
-Currently, the MS-CG library is not setup to run in parallel with MPI,
-so this fix can only be used in a serial LAMMPS build and run
-on a single processor.
 
 Related commands
 """"""""""""""""

doc/src/fix_mvv_dpd.rst

@@ -65,33 +65,37 @@ a default value of 0.5 is used, which effectively reproduces the
 standard velocity-Verlet (VV) scheme.  For more details, see
 :ref:`Groot <Groot2>`.
 
-Fix *mvv/dpd* updates the position and velocity of each atom.  It can
+Fix *mvv/dpd* updates the position and velocity of each atom.  It can be
-be used with the :doc:`pair_style mdpd <pair_mesodpd>` command or other
+used with the :doc:`pair_style mdpd <pair_mesodpd>` command or other
 pair styles such as :doc:`pair dpd <pair_dpd>`.
 
-Fix *mvv/edpd* updates the per-atom temperature, in addition to
+Fix *mvv/edpd* updates the per-atom temperature, in addition to position
-position and velocity, and must be used with the :doc:`pair_style edpd <pair_mesodpd>` command.
+and velocity, and must be used with the :doc:`pair_style edpd
+<pair_mesodpd>` command.
 
-Fix *mvv/tdpd* updates the per-atom chemical concentration, in
+Fix *mvv/tdpd* updates the per-atom chemical concentration, in addition
-addition to position and velocity, and must be used with the
+to position and velocity, and must be used with the :doc:`pair_style
-:doc:`pair_style tdpd <pair_mesodpd>` command.
+tdpd <pair_mesodpd>` command.
 
 ----------
 
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
-No information about this fix is written to :doc:`binary restart files <restart>`.  None of the :doc:`fix_modify <fix_modify>` options
+No information about this fix is written to :doc:`binary restart files
-are relevant to this fix.  No global or per-atom quantities are stored
+<restart>`.  None of the :doc:`fix_modify <fix_modify>` options are
-by this fix for access by various :doc:`output commands <Howto_output>`.
+relevant to this fix.  No global or per-atom quantities are stored by
+this fix for access by various :doc:`output commands <Howto_output>`.
 No parameter of this fix can be used with the *start/stop* keywords of
-the :doc:`run <run>` command.  This fix is not invoked during :doc:`energy minimization <minimize>`.
+the :doc:`run <run>` command.  This fix is not invoked during
+:doc:`energy minimization <minimize>`.
 
 Restrictions
 """"""""""""
 
-This fix is part of the DPD-MESO package. It is only enabled if
+These fixes are part of the DPD-MESO package. They are only enabled if
-LAMMPS was built with that package. See the :doc:`Build package <Build_package>` page for more info.
+LAMMPS was built with that package. See the :doc:`Build package
+<Build_package>` page for more info.
 
 Related commands
 """"""""""""""""

doc/src/fix_nh_uef.rst

@@ -218,10 +218,11 @@ use :doc:`change_box <change_box>` before invoking the fix.
 Related commands
 """"""""""""""""
 
-:doc:`fix nvt <fix_nh>`, :doc:`fix npt <fix_nh>`, `fix nvt/sllod
+:doc:`fix nvt <fix_nh>`, :doc:`fix npt <fix_nh>`,
-:doc:<fix_nvt_sllod>`, `compute temp/uef <compute_temp_uef>`,
+:doc:`fix nvt/sllod <fix_nvt_sllod>`,
-:doc::doc:`compute pressure/uef <compute_pressure_uef>`, `dump cfg/uef
+:doc:`compute temp/uef <compute_temp_uef>`,
-:doc:<dump_cfg_uef>`
+:doc:`compute pressure/uef <compute_pressure_uef>`,
+:doc:`dump cfg/uef <dump_cfg_uef>`
 
 Default
 """""""

doc/src/fix_pimd.rst

@@ -1,7 +1,7 @@
 .. index:: fix pimd/nvt
 
 fix pimd/nvt command
-================
+====================
 
 Syntax
 """"""
@@ -33,7 +33,7 @@ Examples
 Description
 """""""""""
 
-.. versionchanged:: TBD
+.. versionchanged:: 28Mar2023
 
 Fix pimd was renamed to fix pimd/nvt.
 

doc/src/fix_polarize.rst

@@ -31,7 +31,7 @@ Examples
 
    fix 2 interface polarize/bem/gmres 5 0.0001
    fix 1 interface polarize/bem/icc 1 0.0001
-   fix 3 interface polarize/functional 1 0.001
+   fix 3 interface polarize/functional 1 0.0001
 
 
 Used in input scripts:
@@ -69,8 +69,9 @@ along the normal vector is then 78 - 4 = 74, the mean dielectric value
 is (78 + 4) / 2 = 41. Each boundary element also has its area and the
 local mean curvature, which is used by these fixes for computing a
 correction term in the local electric field.  To model charged
-interfaces, the interface particle will have a non-zero charge value,
+interfaces, an interface particle will have a non-zero charge value,
-coming from its area and surface charge density.
+coming from its area and surface charge density, and its local dielectric
+constant set to the mean dielectric value.
 
 For non-interface particles such as atoms and charged particles, the
 interface normal vectors, element area, and dielectric mismatch are
@@ -211,6 +212,8 @@ Note that the *polarize/bem/gmres* and *polarize/bem/icc* fixes only
 support :doc:`units <units>` *lj*, *real*, *metal*, *si* and *nano* at
 the moment.
 
+Note that *polarize/functional* does not yet support charged interfaces.
+
 
 Related commands
 """"""""""""""""
@@ -223,7 +226,7 @@ Related commands
 Default
 """""""
 
-*iter_max* = 20
+*iter_max* = 50
 
 *kspace* = yes
 

doc/src/fix_press_berendsen.rst

@@ -61,7 +61,7 @@ unchanged and controlling the pressure of a surrounding fluid.
    atoms.  This fix can be used in conjunction with thermostatting fixes
    to control the temperature, such as :doc:`fix nvt <fix_nh>` or :doc:`fix langevin <fix_langevin>` or :doc:`fix temp/berendsen <fix_temp_berendsen>`.
 
-See the :doc:`Howto baroostat <Howto_barostat>` page for a
+See the :doc:`Howto barostat <Howto_barostat>` page for a
 discussion of different ways to perform barostatting.
 
 ----------

doc/src/fix_shake.rst

@@ -63,7 +63,7 @@ however, can *only* be applied during molecular dynamics runs.
 
 .. versionchanged:: 15Sep2022
 
-These fixes may still be used during minimization.  In that case the
+These fixes may now also be used during minimization.  In that case the
 constraints are *approximated* by strong harmonic restraints.
 
 **SHAKE vs RATTLE:**
@@ -133,9 +133,9 @@ constraint lists atom types.  All bonds connected to an atom of the
 specified type will be constrained.  The *m* constraint lists atom
 masses.  All bonds connected to atoms of the specified masses will be
 constrained (within a fudge factor of MASSDELTA specified in
-fix_shake.cpp).  The *a* constraint lists angle types.  If both bonds
+``src/RIGID/fix_shake.cpp``).  The *a* constraint lists angle types.  If
-in the angle are constrained then the angle will also be constrained
+both bonds in the angle are constrained then the angle will also be
-if its type is in the list.
+constrained if its type is in the list.
 
 For all constraints, a particular bond is only constrained if both
 atoms in the bond are in the group specified with the SHAKE fix.
@@ -205,11 +205,11 @@ LAMMPS closely follows (:ref:`Andersen (1983) <Andersen3>`).
    The *fix rattle* command modifies forces and velocities and thus
    should be defined after all other integration fixes in your input
    script.  If you define other fixes that modify velocities or forces
-   after *fix rattle* operates, then *fix rattle* will not take them into
+   after *fix rattle* operates, then *fix rattle* will not take them
-   account and the overall time integration will typically not satisfy
+   into account and the overall time integration will typically not
-   the RATTLE constraints.  You can check whether the constraints work
+   satisfy the RATTLE constraints.  You can check whether the
-   correctly by setting the value of RATTLE_DEBUG in src/fix_rattle.cpp
+   constraints work correctly by setting the value of RATTLE_DEBUG in
-   to 1 and recompiling LAMMPS.
+   ``src/RIGID/fix_rattle.cpp`` to 1 and recompiling LAMMPS.
 
 ----------
 
@@ -275,8 +275,8 @@ reducing the :doc:`timestep <timestep>`.
 Related commands
 """"""""""""""""
 
-`fix rigid <fix_rigid>`, `fix ehex <fix_ehex>`,
+:doc:`fix rigid <fix_rigid>`, :doc:`fix ehex <fix_ehex>`,
-`fix nve/manifold/rattle <fix_nve_manifold_rattle>`
+:doc:`fix nve/manifold/rattle <fix_nve_manifold_rattle>`
 
 
 Default

doc/src/fix_wall.rst

@@ -4,7 +4,9 @@
 .. index:: fix wall/lj1043
 .. index:: fix wall/colloid
 .. index:: fix wall/harmonic
+.. index:: fix wall/lepton
 .. index:: fix wall/morse
+.. index:: fix wall/table
 
 fix wall/lj93 command
 =====================
@@ -23,20 +25,31 @@ fix wall/colloid command
 fix wall/harmonic command
 =========================
 
+fix wall/lepton command
+=========================
+
 fix wall/morse command
 ======================
 
+fix wall/table command
+=========================
+
 Syntax
 """"""
 
 .. parsed-literal::
 
-   fix ID group-ID style face args ... keyword value ...
+   fix ID group-ID style [tabstyle] [N] face args ... keyword value ...
 
 * ID, group-ID are documented in :doc:`fix <fix>` command
-* style = *wall/lj93* or *wall/lj126* or *wall/lj1043* or *wall/colloid* or *wall/harmonic* or *wall/morse*
+* style = *wall/lj93* or *wall/lj126* or *wall/lj1043* or *wall/colloid* or *wall/harmonic* or *wall/lepton* or *wall/morse* or *wall/table*
+* tabstyle = *linear* or *spline* = method of table interpolation (only applies to *wall/table*)
+* N = use N values in *linear* or *spline* interpolation (only applies to *wall/table*)
 * one or more face/arg pairs may be appended
 * face = *xlo* or *xhi* or *ylo* or *yhi* or *zlo* or *zhi*
+
+.. spacer
+
 * args for styles *lj93* or *lj126* or *lj1043* or *colloid* or *harmonic*
 
   .. parsed-literal::
@@ -50,7 +63,19 @@ Syntax
            epsilon can be a variable (see below)
          sigma = size factor for wall-particle interaction (distance units)
            sigma can be a variable (see below)
-         cutoff = distance from wall at which wall-particle interaction is cut off (distance units)
+         cutoff = distance from wall at which wall-particle interactions are cut off (distance units)
+
+* args for style *lepton*
+
+  .. parsed-literal::
+
+         args = coord expression cutoff
+         coord = position of wall = EDGE or constant or variable
+           EDGE = current lo or hi edge of simulation box
+           constant = number like 0.0 or -30.0 (distance units)
+           variable = :doc:`equal-style variable <variable>` like v_x or v_wiggle
+         expression = Lepton expression for the potential  (energy units)
+         cutoff = distance from wall at which wall-particle interactions are cut off (distance units)
 
 * args for style *morse*
 
@@ -67,7 +92,20 @@ Syntax
            alpha can be a variable (see below)
          r_0 = distance of the potential minimum from the face of region (distance units)
            r_0 can be a variable (see below)
-         cutoff = distance from wall at which wall-particle interaction is cut off (distance units)
+         cutoff = distance from wall at which wall-particle interactions are cut off (distance units)
+
+* args for style *table*
+
+  .. parsed-literal::
+
+         args = coord filename keyword cutoff
+         coord = position of wall = EDGE or constant or variable
+           EDGE = current lo or hi edge of simulation box
+           constant = number like 0.0 or -30.0 (distance units)
+           variable = :doc:`equal-style variable <variable>` like v_x or v_wiggle
+         filename = file containing tabulated energy and force values
+         keyword = section identifier to select a specific table in table file
+         cutoff = distance from wall at which wall-particle interactions are cut off (distance units)
 
 * zero or more keyword/value pairs may be appended
 * keyword = *units* or *fld* or *pbc*
@@ -91,9 +129,13 @@ Examples
 
    fix wallhi all wall/lj93 xlo -1.0 1.0 1.0 2.5 units box
    fix wallhi all wall/lj93 xhi EDGE 1.0 1.0 2.5
+   fix wallhi all wall/harmonic xhi EDGE 100.0 0.0 4.0 units box
    fix wallhi all wall/morse xhi EDGE 1.0 1.0 1.0 2.5 units box
    fix wallhi all wall/lj126 v_wiggle 23.2 1.0 1.0 2.5
    fix zwalls all wall/colloid zlo 0.0 1.0 1.0 0.858 zhi 40.0 1.0 1.0 0.858
+   fix xwall mobile wall/table spline 200 EDGE -5.0 walltab.dat HARMONIC 4.0
+   fix xwalls mobile wall/lepton xlo -5.0 "k*(r-rc)^2;k=100.0" 4.0 xhi 5.0 "k*(r-rc)^2;k=100.0" 4.0
+
 
 Description
 """""""""""
@@ -103,7 +145,7 @@ wall that interacts with the atoms in the group by generating a force
 on the atom in a direction perpendicular to the wall.  The energy of
 wall-particle interactions depends on the style.
 
-For style *wall/lj93*, the energy E is given by the 9/3 potential:
+For style *wall/lj93*, the energy E is given by the 9-3 Lennard-Jones potential:
 
 .. math::
 
@@ -111,7 +153,7 @@ For style *wall/lj93*, the energy E is given by the 9/3 potential:
                        \left(\frac{\sigma}{r}\right)^3 \right]
                        \qquad r < r_c
 
-For style *wall/lj126*, the energy E is given by the 12/6 potential:
+For style *wall/lj126*, the energy E is given by the 12-6 Lennard-Jones potential:
 
 .. math::
 
@@ -119,7 +161,7 @@ For style *wall/lj126*, the energy E is given by the 12/6 potential:
                        \left(\frac{\sigma}{r}\right)^6 \right]
                        \qquad r < r_c
 
-For style *wall/lj1043*, the energy E is given by the 10/4/3 potential:
+For style *wall/lj1043*, the energy E is given by the 10-4-3 Lennard-Jones potential:
 
 .. math::
 
@@ -138,8 +180,8 @@ of the :doc:`pair_style colloid <pair_colloid>` potential:
     & \left. - \frac{1}{6} \left(\frac{2R(D+R) + D(D+2R)
     \left[ \ln D - \ln (D+2R) \right]}{D(D+2R)} \right) \right] \qquad r < r_c
 
-For style *wall/harmonic*, the energy E is given by a harmonic spring
+For style *wall/harmonic*, the energy E is given by a repulsive-only harmonic
-potential:
+spring potential:
 
 .. math::
 
@@ -152,6 +194,60 @@ For style *wall/morse*, the energy E is given by a Morse potential:
    E = D_0 \left[ e^{- 2 \alpha (r - r_0)} - 2 e^{- \alpha (r - r_0)} \right]
        \qquad r < r_c
 
+.. versionadded:: 28Mar2023
+
+For style *wall/lepton*, the energy E is provided as an Lepton
+expression string using "r" as the distance variable.  The `Lepton
+library <https://simtk.org/projects/lepton>`_, that the *wall/lepton*
+style interfaces with, evaluates this expression string at run time to
+compute the wall-particle energy.  It also creates an analytical
+representation of the first derivative of this expression with respect
+to "r" and then uses that to compute the force between the wall and
+atoms in the fix group.  The Lepton expression must be either enclosed
+in quotes or must not contain any whitespace so that LAMMPS recognizes
+it as a single keyword.
+
+Optionally, the expression may use "rc" to refer to the cutoff distance
+for the given wall.  Further constants in the expression can be defined
+in the same string as additional expressions separated by semi-colons.
+The expression "k*(r-rc)^2;k=100.0" represents a repulsive-only harmonic
+spring as in fix *wall/harmonic* with a force constant *K* (same as
+:math:`\epsilon` above) of 100 energy units.  More details on the Lepton
+expression strings are given below.
+
+.. versionadded:: 28Mar2023
+
+For style *wall/table*, the energy E and forces are determined from
+interpolation tables listed in one or more files as a function of
+distance.  The interpolation tables are used to evaluate energy and
+forces between particles and the wall similar to how analytic formulas
+are used for the other wall styles.
+
+The interpolation tables are created as a pre-computation by fitting
+cubic splines to the file values and interpolating energy and force
+values at each of *N* distances.  During a simulation, the tables are
+used to interpolate energy and force values as needed for each wall and
+particle separated by a distance *R*\ .  The interpolation is done in
+one of two styles: *linear* or *spline*\ .
+
+For the *linear* style, the distance *R* is used to find the 2
+surrounding table values from which an energy or force is computed by
+linear interpolation.
+
+For the *spline* style, cubic spline coefficients are computed and
+stored for each of the *N* values in the table, one set of splines for
+energy, another for force.  Note that these splines are different than
+the ones used to pre-compute the *N* values.  Those splines were fit
+to the *Nfile* values in the tabulated file, where often *Nfile* <
+*N*\ .  The distance *R* is used to find the appropriate set of spline
+coefficients which are used to evaluate a cubic polynomial which
+computes the energy or force.
+
+For each wall a filename and a keyword must be provided as in the
+examples above.  The filename specifies a file containing tabulated
+energy and force values.  The keyword specifies a section of the file.
+The format of this file is described below.
+
 In all cases, *r* is the distance from the particle to the wall at
 position *coord*, and :math:`r_c` is the *cutoff* distance at which the
 particle and wall no longer interact.  The energy of the wall
@@ -180,11 +276,12 @@ box parameters and timestep and elapsed time.  Thus it is easy to
 specify a time-dependent wall position.  See examples below.
 
 For the *wall/lj93* and *wall/lj126* and *wall/lj1043* styles,
-:math:`\epsilon` and :math:`\sigma` are the usual Lennard-Jones parameters, which
+:math:`\epsilon` and :math:`\sigma` are the usual Lennard-Jones
-determine the strength and size of the particle as it interacts with
+parameters, which determine the strength and size of the particle as it
-the wall.  Epsilon has energy units.  Note that this :math:`\epsilon` and
+interacts with the wall.  Epsilon has energy units.  Note that this
-:math:`\sigma` may be different than any :math:`\epsilon` or :math:`\sigma` values defined
+:math:`\epsilon` and :math:`\sigma` may be different than any
-for a pair style that computes particle-particle interactions.
+:math:`\epsilon` or :math:`\sigma` values defined for a pair style that
+computes particle-particle interactions.
 
 The *wall/lj93* interaction is derived by integrating over a 3d
 half-lattice of Lennard-Jones 12/6 particles.  The *wall/lj126*
@@ -207,11 +304,11 @@ are the number density of the constituent particles, in the wall and
 colloid respectively, in units of 1/volume.
 
 The *wall/colloid* interaction is derived by integrating over
-constituent LJ particles of size :math:`\sigma` within the colloid particle
+constituent LJ particles of size :math:`\sigma` within the colloid
-and a 3d half-lattice of Lennard-Jones 12/6 particles of size :math:`\sigma`
+particle and a 3d half-lattice of Lennard-Jones 12/6 particles of size
-in the wall.  As mentioned in the preceding paragraph, the density of
+:math:`\sigma` in the wall.  As mentioned in the preceding paragraph,
-particles in the wall and colloid can be different, as specified by
+the density of particles in the wall and colloid can be different, as
-the :math:`\epsilon` prefactor.
+specified by the :math:`\epsilon` prefactor.
 
 For the *wall/harmonic* style, :math:`\epsilon` is effectively the spring
 constant K, and has units (energy/distance\^2).  The input parameter
@@ -220,20 +317,21 @@ spring is at the *cutoff*\ .  This is a repulsive-only spring since the
 interaction is truncated at the *cutoff*
 
 For the *wall/morse* style, the three parameters are in this order:
-:math:`D_0` the depth of the potential, :math:`\alpha` the width parameter, and
+:math:`D_0` the depth of the potential, :math:`\alpha` the width
-:math:`r_0` the location of the minimum.  :math:`D_0` has energy units, :math:`\alpha`
+parameter, and :math:`r_0` the location of the minimum.  :math:`D_0` has
-inverse distance units, and :math:`r_0` distance units.
+energy units, :math:`\alpha` inverse distance units, and :math:`r_0`
+distance units.
 
-For any wall, the :math:`\epsilon` and/or :math:`\sigma` and/or :math:`\alpha` parameter can
+For any wall that supports them, the :math:`\epsilon` and/or
-be specified
+:math:`\sigma` and/or :math:`\alpha` parameter can be specified as an
-as an :doc:`equal-style variable <variable>`, in which case it should be
+:doc:`equal-style variable <variable>`, in which case it should be
 specified as v_name, where name is the variable name.  As with a
-variable wall position, the variable is evaluated each timestep and
+variable wall position, the variable is evaluated each timestep and the
-the result becomes the current epsilon or sigma of the wall.
+result becomes the current epsilon or sigma of the wall.  Equal-style
-Equal-style variables can specify formulas with various mathematical
+variables can specify formulas with various mathematical functions, and
-functions, and include :doc:`thermo_style <thermo_style>` command
+include :doc:`thermo_style <thermo_style>` command keywords for the
-keywords for the simulation box parameters and timestep and elapsed
+simulation box parameters and timestep and elapsed time.  Thus it is
-time.  Thus it is easy to specify a time-dependent wall interaction.
+easy to specify a time-dependent wall interaction.
 
 .. note::
 
@@ -266,20 +364,19 @@ define the lattice spacings.
 The *fld* keyword can be used with a *yes* setting to invoke the wall
 constraint before pairwise interactions are computed.  This allows an
 implicit FLD model using :doc:`pair_style lubricateU <pair_lubricateU>`
-to include the wall force in its calculations.  If the setting is
+to include the wall force in its calculations.  If the setting is *no*,
-*no*, wall forces are imposed after pairwise interactions, in the
+wall forces are imposed after pairwise interactions, in the usual
-usual manner.
+manner.
 
-The *pbc* keyword can be used with a *yes* setting to allow walls to
+The *pbc* keyword can be used with a *yes* setting to allow walls to be
-be specified in a periodic dimension.  See the
+specified in a periodic dimension.  See the :doc:`boundary <boundary>`
-:doc:`boundary <boundary>` command for options on simulation box
+command for options on simulation box boundaries.  The default for *pbc*
-boundaries.  The default for *pbc* is *no*, which means the system
+is *no*, which means the system must be non-periodic when using a wall.
-must be non-periodic when using a wall.  But you may wish to use a
+But you may wish to use a periodic box.  E.g. to allow some particles to
-periodic box.  E.g. to allow some particles to interact with the wall
+interact with the wall via the fix group-ID, and others to pass through
-via the fix group-ID, and others to pass through it and wrap around a
+it and wrap around a periodic box.  In this case you should ensure that
-periodic box.  In this case you should ensure that the wall if
+the wall is sufficiently far enough away from the box boundary.  If you
-sufficiently far enough away from the box boundary.  If you do not,
+do not, then particles may interact with both the wall and with periodic
-then particles may interact with both the wall and with periodic
 images on the other side of the box, which is probably not what you
 want.
 
@@ -328,6 +425,57 @@ perturbation on the particles:
 
 ----------
 
+.. include:: lepton_expression.rst
+
+----------
+
+Table file format
+"""""""""""""""""
+
+Suitable tables for use with fix *wall/table* can be created by the
+Python code in the ``tools/tabulate`` folder of the LAMMPS source code
+distribution.
+
+The format of a tabulated file is as follows (without the parenthesized
+comments):
+
+.. parsed-literal::
+
+   # Tabulated wall potential UNITS: real
+
+   HARMONIC                       (keyword is the first text on a line)
+   N 100 FP 200 200
+                                  (blank line)
+   1  0.04    1568.16    792.00   (index, distance to wall, energy, force)
+   2  0.08    1536.64    784.00
+   3  0.12    1505.44    776.00
+   ...
+   99  3.96      0.16      8.00
+   100  4.00     0         0
+
+A section begins with a non-blank line whose first character is not a
+"#"; blank lines or lines starting with "#" can be used as comments
+between sections.  The first line begins with a keyword which identifies
+the section.  The line can contain additional text, but the initial text
+must match the argument specified in the fix *wall/table* command.  The
+next line lists (in any order) one or more parameters for the table.
+Each parameter is a keyword followed by one or more numeric values.
+
+The parameter "N" is required and its value is the number of table
+entries that follow.  Note that this may be different than the *N*
+specified in the fix *wall/table* command.  Let Ntable = *N* in the fix
+command, and Nfile = "N" in the tabulated file.  What LAMMPS does is a
+preliminary interpolation by creating splines using the Nfile tabulated
+values as nodal points.  It uses these to interpolate as needed to
+generate energy and force values at Ntable different points.  The
+resulting tables of length Ntable are then used as described above, when
+computing energy and force for wall-particle interactions.  This means that
+if you want the interpolation tables of length Ntable to match exactly
+what is in the tabulated file (with effectively no preliminary
+interpolation), you should set Ntable = Nfile.
+
+----------
+
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
@@ -354,16 +502,15 @@ fix. This allows to set at which level of the :doc:`r-RESPA
 <run_style>` integrator the fix is adding its forces. Default is the
 outermost level.
 
-This fix computes a global scalar energy and a global vector of
+This fix computes a global scalar energy and a global vector of forces,
-forces, which can be accessed by various :doc:`output commands
+which can be accessed by various :doc:`output commands <Howto_output>`.
-<Howto_output>`.  Note that the scalar energy is the sum of
+Note that the scalar energy is the sum of interactions with all defined
-interactions with all defined walls.  If you want the energy on a
+walls.  If you want the energy on a per-wall basis, you need to use
-per-wall basis, you need to use multiple fix wall commands.  The
+multiple fix wall commands.  The length of the vector is equal to the
-length of the vector is equal to the number of walls defined by the
+number of walls defined by the fix.  Each vector value is the normal
-fix.  Each vector value is the normal force on a specific wall.  Note
+force on a specific wall.  Note that an outward force on a wall will be
-that an outward force on a wall will be a negative value for *lo*
+a negative value for *lo* walls and a positive value for *hi* walls.
-walls and a positive value for *hi* walls.  The scalar and vector
+The scalar and vector values calculated by this fix are "extensive".
-values calculated by this fix are "extensive".
 
 No parameter of this fix can be used with the *start/stop* keywords of
 the :doc:`run <run>` command.
@@ -386,7 +533,11 @@ invoked by the :doc:`minimize <minimize>` command.
 
 Restrictions
 """"""""""""
- none
+
+Fix *wall/lepton* is part of the LEPTON package and only enabled if
+LAMMPS was built with this package.  See the :doc:`Build package
+<Build_package>` page for more info.
+
 
 Related commands
 """"""""""""""""

doc/src/fix_wall_gran.rst

@@ -46,7 +46,7 @@ Syntax
          radius = cylinder radius (distance units)
 
 * zero or more keyword/value pairs may be appended to args
-* keyword = *wiggle* or *shear* or *contacts*
+* keyword = *wiggle* or *shear* or *contacts* or *temperature*
 
   .. parsed-literal::
 
@@ -59,6 +59,8 @@ Syntax
          vshear = magnitude of shear velocity (velocity units)
       *contacts* value = none
          generate contact information for each particle
+      *temperature* value = temperature
+         specify temperature of wall
 
 
 Examples
@@ -71,7 +73,7 @@ Examples
    fix 2 all wall/gran hooke 100000.0 20000.0 50.0 30.0 0.5 1 zcylinder 15.0 wiggle z 3.0 2.0
    fix 3 all wall/gran/region granular hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 damping velocity region myBox
    fix 4 all wall/gran/region granular jkr 1e5 1500.0 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall region myCone
-   fix 5 all wall/gran/region granular dmt 1e5 0.2 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall damping tsuji region myCone
+   fix 5 all wall/gran/region granular dmt 1e5 0.2 0.3 10.0 tangential mindlin NULL 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall damping tsuji heat 10 region myCone temperature 1.0
    fix 6 all wall/gran hooke  200000.0 NULL 50.0 NULL 0.5 0 xplane -10.0 10.0 contacts
 
 Description
@@ -177,6 +179,16 @@ the clockwise direction for *vshear* > 0 or counter-clockwise for
 *vshear* < 0.  In this case, *vshear* is the tangential velocity of
 the wall at whatever *radius* has been defined.
 
+The *temperature* keyword is used to assign a temperature to the wall.
+The following value can either be a numeric value or an equal-style
+:doc:`variable <variable>`.  If the value is a variable, it should be
+specified as v_name, where name is the variable name.  In this case, the
+variable will be evaluated each timestep, and its value used to determine
+the temperature. This option must be used in conjunction with a heat
+conduction model defined in :doc:`pair_style granular <pair_granular>`,
+:doc:`fix property/atom <fix_property_atom>` to store temperature and a
+heat flow, and :doc:`fix heat/flow <fix_heat_flow>` to integrate heat
+flow.
 
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

doc/src/fix_wall_gran_region.rst

@@ -36,12 +36,14 @@ Syntax
 
 * wallstyle = region (see :doc:`fix wall/gran <fix_wall_gran>` for options for other kinds of walls)
 * region-ID = region whose boundary will act as wall
-* keyword = *contacts*
+* keyword = *contacts* or *temperature*
 
   .. parsed-literal::
 
       *contacts* value = none
          generate contact information for each particle
+      *temperature* value = temperature
+         specify temperature of wall
 
 Examples
 """"""""
@@ -200,6 +202,17 @@ values for the 6 wall/particle coefficients than for particle/particle
 interactions.  E.g. if you wish to model the wall as a different
 material.
 
+The *temperature* keyword is used to assign a temperature to the wall.
+The following value can either be a numeric value or an equal-style
+:doc:`variable <variable>`.  If the value is a variable, it should be
+specified as v_name, where name is the variable name.  In this case, the
+variable will be evaluated each timestep, and its value used to determine
+the temperature. This option must be used in conjunction with a heat
+conduction model defined in :doc:`pair_style granular <pair_granular>`,
+:doc:`fix property/atom <fix_property_atom>` to store temperature and a
+heat flow, and :doc:`fix heat/flow <fix_heat_flow>` to integrate heat
+flow.
+
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 

doc/src/improper_amoeba.rst

@@ -1,7 +1,7 @@
 .. index:: improper_style amoeba
 
-improper_style harmonic command
+improper_style amoeba command
-===============================
+=============================
 
 Syntax
 """"""

doc/src/mdi.rst

@@ -17,7 +17,7 @@ Syntax
      *engine* args = zero or more keyword/args pairs
        keywords = *elements*
          *elements* args = N_1 N_2 ... N_ntypes
-           N_1,N_2,...N_ntypes = atomic number for each of ntypes LAMMPS atom types
+           N_1,N_2,...N_ntypes = chemical symbol for each of ntypes LAMMPS atom types
      *plugin* args = name keyword value keyword value ...
        name = name of plugin library (e.g., *lammps* means a liblammps.so library will be loaded)
        keyword/value pairs in any order, some are required, some are optional
@@ -35,7 +35,7 @@ Examples
 .. code-block:: LAMMPS
 
    mdi engine
-   mdi engine elements 13 29
+   mdi engine elements Al Cu
    mdi plugin lammps mdi "-role ENGINE -name lammps -method LINK" &
               infile in.aimd.engine extra "-log log.aimd.engine.plugin" &
               command "run 5"
@@ -173,13 +173,16 @@ commands, which are described further below.
    atom type values are consistent in both codes, then the >TYPES
    command can be used.  If not, the optional *elements* keyword can
    be used to specify what element each LAMMPS atom type corresponds
-   to.  This is specified by the atomic number of the element (e.g., 13
+   to.  This is specified by the chemical symbol of the element,
-   for Al).  An atomic number must be specified for each of the ntypes
+   e.g. C or Al or Si.  A symbol must be specified for each of the
-   LAMMPS atom types.  Ntypes is typically specified via the
+   ntypes LAMMPS atom types.  Each LAMMPS type must map to a unique
-   create_box command or in the data file read by the read_data
+   element; two or more types cannot map to the same element.  Ntypes
-   command.  In this has been done, the MDI driver can send an
+   is typically specified via the :doc:`create_box <create_box>`
-   >ELEMENTS command to the LAMMPS driver with the atomic number of
+   command or in the data file read by the :doc:`read_data
-   each atom.
+   <read_data>` command.  Once this has been done, the MDI driver can
+   send an >ELEMENTS command to the LAMMPS driver with the atomic
+   number of each atom and the LAMMPS engine will be able to map it to
+   a LAMMPS atom type.
 
 The MD and OPTG commands perform an entire MD simulation or energy
 minimization (to convergence) with no communication from the driver

doc/src/minimize.rst

@@ -49,19 +49,19 @@ and forces) by pushing the atoms off of each other.
 
    The distance that atoms can move during individual minimization steps
    can be quite large, especially at the beginning of a minimization.
-   Thus `neighbor list settings <neigh_modify>` of *every = 1* and
+   Thus :doc:`neighbor list settings <neigh_modify>` of *every = 1* and
    *delay = 0* are **required**.  This may be combined with either
    *check = no* (always update the neighbor list) or *check = yes* (only
    update the neighbor list if at least one atom has moved more than
-   half the `neighbor list skin <neighbor>` distance since the last
+   half the :doc:`neighbor list skin <neighbor>` distance since the last
    reneighboring).  Using *check = yes* is recommended since it avoids
    unneeded reneighboring steps when the system is closer to the minimum
-   and thus atoms move only small distances.  Using *check = no* may
+   and thus atoms move only small distances.  Using *check = no* may be
-   be required for debugging or when coupling LAMMPS with external
+   required for debugging or when coupling LAMMPS with external codes
-   codes that require a predictable sequence of neighbor list updates.
+   that require a predictable sequence of neighbor list updates.
 
-   If the settings are **not** *every = 1* and *delay = 0*, LAMMPS
+   If the settings are **not** *every = 1* and *delay = 0*, LAMMPS will
-   will temporarily apply a `neigh_modify every 1 delay 0 check yes
+   temporarily apply a :doc:`neigh_modify every 1 delay 0 check yes
    <neigh_modify>` setting during the minimization and restore the
    original setting at the end of the minimization.  A corresponding
    message will be printed to the screen and log file, if this happens.

doc/src/molecule.rst

@@ -121,6 +121,11 @@ molecule (header keyword = inertia).
    ensure space is allocated for storing topology info for molecules that
    are added later.
 
+----------
+
+Format of a molecule file
+"""""""""""""""""""""""""
+
 The format of an individual molecule file is similar but
 (not identical) to the data file read by the :doc:`read_data <read_data>`
 commands, and is as follows.

doc/src/neighbor.rst

@@ -59,9 +59,21 @@ long cutoff, but other type pairs have a much shorter cutoff. The
 sized particles, where "size" may mean the physical size of the particle
 or its cutoff distance for interacting with other particles. Different
 sets of bins are then used to construct the neighbor lists as as further
-described by Shire, Hanley, and Stratford :ref:`(Shire) <bytype-Shire>`.
+described by Shire, Hanley, and Stratford :ref:`(Shire) <multi-Shire>`
-This imposes some extra setup overhead, but the searches themselves may
+and Monti et al. :ref:`(Monti) <multi-Monti>`. This imposes some extra
-be much faster. By default, each atom type defines a separate collection
+setup overhead, but the searches themselves may be much faster.
+
+For instance in a dense binary system in d-dimensions with a ratio of the size
+of the largest to smallest collection bin :math:`\lambda`, the computational
+costs of building a default neighbor list grows as :math:`\lambda^{2d}` while
+the costs for *multi* grows as :math:`\lambda^d`, equivalent to the cost
+of force evaluations, as argued in Monti et al. :ref:`(Monti) <multi-Monti>`.
+In other words, the neighboring costs of *multi* are expected to scale the
+same as force calculations, such that its relative cost is independent of
+the particle size ratio. This is not the case for the default style which
+becomes substantially more expensive with increasing size ratios.
+
+By default in *multi*, each atom type defines a separate collection
 of particles. For systems where two or more atom types have the same
 size (either physical size or cutoff distance), the definition of
 collections can be customized, which can result in less overhead and
@@ -75,8 +87,11 @@ An alternate style, *multi/old*, sets the bin size to 1/2 of the shortest
 cutoff distance and multiple sets of bins are defined to search over for
 different atom types. This algorithm used to be the default *multi*
 algorithm in LAMMPS but was found to be significantly slower than the new
-approach. For now we are keeping the old option in case there are use cases
+approach. For the dense binary system, computational costs still grew as
-where multi/old outperforms the new multi style.
+:math:`\lambda^{2d}` at large enough :math:`\lambda`. This is equivalent
+to the default style, albeit with a smaller prefactor. For now we are
+keeping the old option in case there are use cases where multi/old
+outperforms the new multi style.
 
 .. note::
 
@@ -118,6 +133,10 @@ Default
 
 ----------
 
-.. _bytype-Shire:
+.. _multi-Shire:
 
-**(Shire)** Shire, Hanley and Stratford, Comp Part Mech, (2020).
+**(Shire)** Shire, Hanley and Stratford, Comp. Part. Mech., (2020).
+
+.. _multi-Monti:
+
+**(Monti)** Monti, Clemmer, Srivastava, Silbert, Grest, and Lechman, Phys. Rev. E, (2022).

doc/src/package.rst

@@ -319,7 +319,7 @@ CONFIG_ID, SIMD_SIZE, MEM_THREADS, SHUFFLE_AVAIL, FAST_MATH,
 THREADS_PER_ATOM, THREADS_PER_CHARGE, THREADS_PER_THREE, BLOCK_PAIR,
 BLOCK_BIO_PAIR, BLOCK_ELLIPSE, PPPM_BLOCK_1D, BLOCK_NBOR_BUILD,
 BLOCK_CELL_2D, BLOCK_CELL_ID, MAX_SHARED_TYPES, MAX_BIO_SHARED_TYPES,
-PPPM_MAX_SPLINE.
+PPPM_MAX_SPLINE, NBOR_PREFETCH.
 
 CONFIG_ID can be 0. SHUFFLE_AVAIL in {0,1} indicates that inline-PTX
 (NVIDIA) or OpenCL extensions (Intel) should be used for horizontal

doc/src/pair_born_gauss.rst

@@ -0,0 +1,103 @@
+.. index:: pair_style born/gauss
+
+pair_style born/gauss command
+=============================
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   pair_style born/gauss cutoff
+
+* born/gauss = name of the pair style
+* cutoff = global cutoff (distance units)
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   pair_style born/gauss 10.0
+   pair_coeff 1 1 1 1 8.2464e13 12.48 0.042644277 0.44 3.56
+
+Description
+"""""""""""
+
+.. versionadded:: 28Mar2023
+
+Pair style *born/gauss* computes pairwise interactions from a combination of a Born-Mayer
+repulsive term and a Gaussian attractive term according to :ref:`(Bomont) <Bomont>`:
+
+.. math::
+
+   E = A_0 \exp \left( -\alpha r \right) - A_1 \exp\left[ -\beta \left(r - r_0 \right)^2 \right]
+       \qquad r < r_c
+
+:math:`r_c` is the cutoff.
+
+The following coefficients 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:
+
+* :math:`A_0` (energy units)
+* :math:`\alpha` (1/distance units)
+* :math:`A_1` (energy units)
+* :math:`\beta` (1/(distance units)^2)
+* :math:`r_0` (distance units)
+* cutoff (distance units)
+
+The last coefficient is optional.  If not specified, the global cutoff is used.
+
+----------
+
+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 supports the :doc:`pair_modify <pair_modify>` shift
+option for the energy of the pair interaction.
+
+The :doc:`pair_modify <pair_modify>` table options are 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.
+
+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 pair style is only enabled if LAMMPS was built with the EXTRA-PAIR
+package.  See the :doc:`Build package <Build_package>` page for more
+info.
+
+Related commands
+""""""""""""""""
+
+:doc:`pair_coeff <pair_coeff>`, :doc:`pair_style born <pair_born>`
+
+Default
+"""""""
+
+none
+
+--------------
+
+.. _Bomont:
+
+**(Bomont)** Bomont, Bretonnet, J. Chem. Phys. 124, 054504 (2006)

doc/src/pair_eam.rst

@@ -116,11 +116,12 @@ are parameterized in terms of LAMMPS :doc:`metal units <units>`.
 
 .. note::
 
-   Note that unlike for other potentials, cutoffs for EAM
+   Note that unlike for other potentials, cutoffs for EAM potentials are not
-   potentials are not set in the pair_style or pair_coeff command; they
+   set in the pair_style or pair_coeff command; they are specified in the EAM
-   are specified in the EAM potential files themselves.  Likewise, the
+   potential files themselves.  Likewise, valid EAM potential files usually
-   EAM potential files list atomic masses; thus you do not need to use
+   contain atomic masses; thus you may not need to use the :doc:`mass <mass>`
-   the :doc:`mass <mass>` command to specify them.
+   command to specify them, unless the potential file uses a dummy value (e.g.
+   0.0). LAMMPS will print a warning, if this is the case.
 
 There are web sites that distribute and document EAM potentials stored
 in DYNAMO or other formats:

doc/src/pair_gauss.rst

@@ -103,12 +103,13 @@ Mixing, shift, table, tail correction, restart, rRESPA info
 
 For atom type pairs I,J and I != J, the A, B, H, sigma_h, r_mh
 parameters, and the cutoff distance for these pair styles can be mixed:
-A (energy units)
+
-sqrt(1/B) (distance units, see below)
+* A (energy units)
-H (energy units)
+* :math:`\sqrt{\frac{1}{B}}` (distance units, see below)
-sigma_h (distance units)
+* H (energy units)
-r_mh (distance units)
+* :math:`r_{mh}` (distance units)
-cutoff (distance units):ul
+* :math:`\sigma_h` (distance units)
+* cutoff (distance units)
 
 The default mix value is *geometric*\ .
 Only *arithmetic* and *geometric* mix values are supported.
@@ -169,7 +170,7 @@ The *gauss* and *gauss/cut* styles are part of the EXTRA-PAIR package.
 They are only enabled if LAMMPS is build with that package.  See the
 :doc:`Build package <Build_package>` page for more info.
 
-.. versionchanged:: TBD
+.. versionchanged:: 28Mar2023
 
 Prior to this version, the *gauss* pair style did not apply
 :doc:`special_bonds <special_bonds>` factors.

doc/src/pair_gran.rst

@@ -279,6 +279,11 @@ 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.
 
+Accumulated values for individual contacts are saved to to restart
+files but are not saved to data files. Therefore, forces may
+differ significantly when a system is reloaded using A
+:doc:`read_data <read_data>` command.
+
 Related commands
 """"""""""""""""
 

doc/src/pair_granular.rst

@@ -37,6 +37,9 @@ Examples
    pair_coeff 1 1 dmt 1000.0 50.0 0.3 0.0 tangential mindlin NULL 0.5 0.5 rolling sds 500.0 200.0 0.5 twisting marshall
    pair_coeff 2 2 dmt 1000.0 50.0 0.3 10.0 tangential mindlin NULL 0.5 0.1 rolling sds 500.0 200.0 0.1 twisting marshall
 
+   pair_style granular
+   pair_coeff * * hertz 1000.0 50.0 tangential mindlin 1000.0 1.0 0.4 heat area 0.1
+
 Description
 """""""""""
 
@@ -223,16 +226,19 @@ for all models except *jkr*, for which it is given implicitly according
 to :math:`\delta = a^2/R - 2\sqrt{\pi \gamma a/E}`.  For *damping viscoelastic*,
 :math:`\eta_{n0}` is in units of 1/(\ *time*\ \*\ *distance*\ ).
 
-The *tsuji* model is based on the work of :ref:`(Tsuji et al) <Tsuji1992>`. Here, the damping coefficient specified as part of
+The *tsuji* model is based on the work of :ref:`(Tsuji et al) <Tsuji1992>`.
-the normal model is interpreted as a restitution coefficient
+Here, the damping coefficient specified as part of the normal model is interpreted
-:math:`e`. The damping constant :math:`\eta_n` is given by:
+as a restitution coefficient :math:`e`. The damping constant :math:`\eta_n` is
+given by:
 
 .. math::
 
    \eta_n = \alpha (m_{eff}k_n)^{1/2}
 
-For normal contact models based on material parameters, :math:`k_n = 4/3Ea`.  The parameter :math:`\alpha` is related to the restitution
+For normal contact models based on material parameters, :math:`k_n = 4/3Ea`. This
-coefficient *e* according to:
+damping model is not compatible with cohesive normal models such as *JKR* or *DMT*.
+The parameter :math:`\alpha` is related to the restitution coefficient *e*
+according to:
 
 .. math::
 
@@ -631,6 +637,34 @@ attractive force. This keyword cannot be used with the JKR or DMT models.
 
 ----------
 
+The optional *heat* keyword enables heat conduction. The options currently
+supported are:
+
+1. *none*
+2. *area* : :math:`k_{s}`
+
+If the *heat* keyword is not specified, the model defaults to *none*.
+
+For *heat* *area*, the heat
+:math:`Q` conducted between two particles is given by
+
+.. math::
+
+   Q = k_{s} a \Delta T
+
+
+
+where :math:`\Delta T` is the difference in the two particles' temperature,
+:math:`k_{s}` is a non-negative numeric value for the conductivity, and
+:math:`a` is the area of the contact and depends on the normal force model.
+
+Note that the option *none* must either be used in all or none of of the
+*pair_coeff* calls. See :doc:`fix heat/flow <fix_heat_flow>` and
+:doc:`fix property/atom <fix_property_atom>` for more information on this
+option.
+
+----------
+
 The *granular* pair style can reproduce the behavior of the
 *pair gran/\** styles with the appropriate settings (some very
 minor differences can be expected due to corrections in
@@ -641,7 +675,7 @@ is equivalent to *pair gran/hooke 1000.0 NULL 50.0 50.0 0.4 1*\ .
 The second example is equivalent to
 *pair gran/hooke/history 1000.0 500.0 50.0 50.0 0.4 1*\ .
 The third example is equivalent to
-*pair gran/hertz/history 1000.0 500.0 50.0 50.0 0.4 1*\ .
+*pair gran/hertz/history 1000.0 500.0 50.0 50.0 0.4 1 limit_damping*\ .
 
 ----------
 
@@ -733,22 +767,30 @@ These extra quantities can be accessed by the :doc:`compute pair/local <compute_
 Restrictions
 """"""""""""
 
-All the granular pair styles are part of the GRANULAR package.  It is
+This pair style is part of the GRANULAR package.  It is
-only enabled if LAMMPS was built with that package.  See the :doc:`Build package <Build_package>` page for more info.
+only enabled if LAMMPS was built with that package.
+See the :doc:`Build package <Build_package>` page for more info.
 
-These pair styles require that atoms store torque and angular velocity
+This pair style requires that atoms store per-particle radius,
-(omega) as defined by the :doc:`atom_style <atom_style>`.  They also
+torque, and angular velocity (omega) as defined by the
-require a per-particle radius is stored.  The *sphere* atom style does
+:doc:`atom_style sphere <atom_style>`.
-all of this.
 
-This pair style requires you to use the :doc:`comm_modify vel yes <comm_modify>` command so that velocities are stored by ghost
+This pair style requires you to use the :doc:`comm_modify vel yes <comm_modify>`
-atoms.
+command so that velocities are stored by ghost atoms.
 
-These pair styles will not restart exactly when using the
+This pair style will not restart exactly when using the
-:doc:`read_restart <read_restart>` command, though they should provide
+:doc:`read_restart <read_restart>` command, though it should provide
-statistically similar results.  This is because the forces they
+statistically similar results.  This is because the forces it
-compute depend on atom velocities.  See the
+computes depend on atom velocities and the atom velocities have
-:doc:`read_restart <read_restart>` command for more details.
+been propagated half a timestep between the force computation and
+when the restart is written, due to using Velocity Verlet time
+integration. See the :doc:`read_restart <read_restart>` command
+for more details.
+
+Accumulated values for individual contacts are saved to restart
+files but are not saved to data files. Therefore, forces may
+differ significantly when a system is reloaded using the
+:doc:`read_data <read_data>` command.
 
 Related commands
 """"""""""""""""

doc/src/pair_hybrid.rst

@@ -91,7 +91,8 @@ contributions from sub-styles are weighted by their scale factors, which
 may be fractional or even negative.  Furthermore the scale factors may
 be variables that may change during a simulation.  This enables
 switching smoothly between two different pair styles or two different
-parameter sets during a run.
+parameter sets during a run in a similar fashion as could be done
+with :doc:`fix adapt <fix_adapt>` or :doc:`fix alchemy <fix_alchemy>`.
 
 All pair styles that will be used are listed as "sub-styles" following
 the *hybrid* or *hybrid/overlay* keyword, in any order.  In case of the

doc/src/pair_lj.rst

@@ -68,7 +68,7 @@ Note that :math:`\sigma` is defined in the LJ formula as the zero-crossing
 distance for the potential, not as the energy minimum at :math:`2^{\frac{1}{6}} \sigma`.
 
 The last coefficient is optional.  If not specified, the global
-LJ cutoff specified in the pair_style command are used.
+LJ cutoff specified in the pair_style command is used.
 
 ----------
 
@@ -103,11 +103,12 @@ portion of the pair interaction.
 All of the *lj/cut* pair styles write their 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.
 
-The *lj/cut* pair styles support the use of the
+The *lj/cut* pair styles support the use of the *inner*, *middle*, and
-*inner*, *middle*, and *outer* keywords of the :doc:`run_style respa <run_style>` command, meaning the pairwise forces can be
+*outer* keywords of the :doc:`run_style respa <run_style>` command,
-partitioned by distance at different levels of the rRESPA hierarchy.
+meaning the pairwise forces can be partitioned by distance at different
-The other styles only support the *pair* keyword of run_style respa.
+levels of the rRESPA hierarchy.  The other styles only support the
-See the :doc:`run_style <run_style>` command for details.
+*pair* keyword of run_style respa.  See the :doc:`run_style <run_style>`
+command for details.
 
 ----------
 

doc/src/pair_mesodpd.rst

@@ -299,7 +299,7 @@ Restrictions
 """"""""""""
 
 The pair styles *edpd*, *mdpd*, *mdpd/rhosum* and *tdpd* are part of
-the DPD-MESO package. It is only enabled if LAMMPS was built with
+the DPD-MESO package. They are only enabled if LAMMPS was built with
 that package.  See the :doc:`Build package <Build_package>` page for
 more info.
 

doc/src/pair_style.rst

@@ -132,6 +132,7 @@ accelerated styles exist.
 * :doc:`born/coul/msm <pair_born>` - Born with long-range MSM Coulomb
 * :doc:`born/coul/wolf <pair_born>` - Born with Wolf potential for Coulomb
 * :doc:`born/coul/wolf/cs <pair_cs>` - Born with Wolf potential for Coulomb and core/shell model
+* :doc:`born/gauss <pair_born_gauss>` - Born-Mayer / Gaussian potential
 * :doc:`bpm/spring <pair_bpm_spring>` - repulsive harmonic force with damping
 * :doc:`brownian <pair_brownian>` - Brownian potential for Fast Lubrication Dynamics
 * :doc:`brownian/poly <pair_brownian>` - Brownian potential for Fast Lubrication Dynamics with polydispersity

doc/src/pair_tersoff_mod.rst

@@ -145,7 +145,7 @@ coefficients in the formulae above:
 * c3
 * c4
 * c5
-* c0 (energy units, tersoff/mod/c only):ul
+* c0 (energy units, tersoff/mod/c only)
 
 The n, :math:`\eta`, :math:`\lambda_2`, B, :math:`\lambda_1`, and A parameters are only used for
 two-body interactions.  The :math:`\beta`, :math:`\alpha`, c1, c2, c3, c4, c5, h

doc/src/pair_ylz.rst

@@ -142,10 +142,11 @@ the :doc:`atom_style ellipsoid <atom_style>` command.
 Related commands
 """"""""""""""""
 
-:doc:`pair_coeff <pair_coeff>`, :doc:`fix nve/asphere
+:doc:`pair_coeff <pair_coeff>`,
-:doc:<fix_nve_asphere>`, `compute temp/asphere <compute_temp_asphere>`,
+:doc:`fix nve/asphere <fix_nve_asphere>`,
-:doc::doc:`pair_style resquared <pair_resquared>`, :doc:`pair_style
+:doc:`compute temp/asphere <compute_temp_asphere>`,
-:doc:gayberne <pair_gayberne>`
+:doc:`pair_style resquared <pair_resquared>`,
+:doc:`pair_style gayberne <pair_gayberne>`
 
 Default
 """""""

doc/src/read_data.rst

@@ -1477,10 +1477,12 @@ The *Triangles* section must appear after the *Atoms* section.
 
 where the keywords have these meanings:
 
-vx,vy,vz = translational velocity of atom
+   .. parsed-literal::
-lx,ly,lz = angular momentum of aspherical atom
+
-wx,wy,wz = angular velocity of spherical atom
+      vx,vy,vz = translational velocity of atom
-ervel = electron radial velocity (0 for fixed-core):ul
+      lx,ly,lz = angular momentum of aspherical atom
+      wx,wy,wz = angular velocity of spherical atom
+      ervel = electron radial velocity (0 for fixed-core)
 
 The velocity lines can appear in any order.  This section can only be
 used after an *Atoms* section.  This is because the *Atoms* section

doc/src/run_style.rst

@@ -78,27 +78,41 @@ processors.  See the :doc:`-partition command-line switch <Run_options>`
 for info on how to run LAMMPS with multiple partitions.
 
 Specifically, this style performs all computation except the
-:doc:`kspace_style <kspace_style>` portion of the force field on the first
+:doc:`kspace_style <kspace_style>` portion of the force field on the
-partition.  This include the :doc:`pair style <pair_style>`, :doc:`bond style <bond_style>`, :doc:`neighbor list building <neighbor>`,
+first partition.  This include the :doc:`pair style <pair_style>`,
-:doc:`fixes <fix>` including time integration, and output.  The
+:doc:`bond style <bond_style>`, :doc:`neighbor list building
-:doc:`kspace_style <kspace_style>` portion of the calculation is
+<neighbor>`, :doc:`fixes <fix>` including time integration, and output.
+The :doc:`kspace_style <kspace_style>` portion of the calculation is
 performed on the second partition.
 
-This is most useful for the PPPM kspace_style when its performance on
+This can lead to a significant speedup, if the number of processors can
-a large number of processors degrades due to the cost of communication
+be easily increased and the fraction of time is spent in computing
-in its 3d FFTs.  In this scenario, splitting your P total processors
+Kspace interactions is significant, too.  The two partitions may have a
-into 2 subsets of processors, P1 in the first partition and P2 in the
+different number of processors.  This is most useful for the PPPM
-second partition, can enable your simulation to run faster.  This is
+kspace_style when its performance on a large number of processors
-because the long-range forces in PPPM can be calculated at the same
+degrades due to the cost of communication in its 3d FFTs.  In this
-time as pairwise and bonded forces are being calculated, and the FFTs
+scenario, splitting your P total processors into 2 subsets of
-can actually speed up when running on fewer processors.
+processors, P1 in the first partition and P2 in the second partition,
+can enable your simulation to run faster.  This is because the
+long-range forces in PPPM can be calculated at the same time as pairwise
+and bonded forces are being calculated *and* the parallel 3d FFTs can be
+faster to compute when running on fewer processors.  Please note that
+the scenario of using fewer MPI processes to reduce communication
+overhead can also be implemented through using MPI with OpenMP threads
+via the INTEL, KOKKOS, or OPENMP package.  This alternative option is
+typically more effective in case of a fixed number of available
+processors and less complex to execute.
 
-To use this style, you must define 2 partitions where P1 is a multiple
+To use the *verlet/split* style, you must define 2 partitions with the
-of P2.  Typically having P1 be 3x larger than P2 is a good choice.
+:doc:`-partition command-line switch <Run_options>`, where partition P1
-The 3d processor layouts in each partition must overlay in the
+is either the same size or an integer multiple of the size of the
-following sense.  If P1 is a Px1 by Py1 by Pz1 grid, and P2 = Px2 by
+partition P2.  Typically having P1 be 3x larger than P2 is a good
-Py2 by Pz2, then Px1 must be an integer multiple of Px2, and similarly
+choice, since the (serial) performance of LAMMPS is often best if the
-for Py1 a multiple of Py2, and Pz1 a multiple of Pz2.
+time spent in the ``Pair`` computation versus ``Kspace`` is a 3:1 split.
+The 3d processor layouts in each partition must overlay in the following
+sense.  If P1 is a Px1 by Py1 by Pz1 grid, and P2 = Px2 by Py2 by Pz2,
+then Px1 must be an integer multiple of Px2, and similarly for Py1 a
+multiple of Py2, and Pz1 a multiple of Pz2.
 
 Typically the best way to do this is to let the first partition choose
 its own optimal layout, then require the second partition's layout to
@@ -122,9 +136,10 @@ of 60 and 15 processors each:
 When you run in 2-partition mode with the *verlet/split* style, the
 thermodynamic data for the entire simulation will be output to the log
 and screen file of the first partition, which are log.lammps.0 and
-screen.0 by default; see the :doc:`-plog and -pscreen command-line switches <Run_options>` to change this.  The log and screen file
+screen.0 by default; see the :doc:`-plog and -pscreen command-line
-for the second partition will not contain thermodynamic output beyond the
+switches <Run_options>` to change this.  The log and screen file for the
-first timestep of the run.
+second partition will not contain thermodynamic output beyond the first
+timestep of the run.
 
 See the :doc:`Accelerator packages <Speed_packages>` page for
 performance details of the speed-up offered by the *verlet/split*
@@ -137,70 +152,73 @@ options to support this, and strategies are discussed in :doc:`Section
 ----------
 
 The *respa* style implements the rRESPA multi-timescale integrator
-:ref:`(Tuckerman) <Tuckerman3>` with N hierarchical levels, where level 1 is
+:ref:`(Tuckerman) <Tuckerman3>` with N hierarchical levels, where level
-the innermost loop (shortest timestep) and level N is the outermost
+1 is the innermost loop (shortest timestep) and level N is the outermost
 loop (largest timestep).  The loop factor arguments specify what the
-looping factor is between levels.  N1 specifies the number of
+looping factor is between levels.  N1 specifies the number of iterations
-iterations of level 1 for a single iteration of level 2, N2 is the
+of level 1 for a single iteration of level 2, N2 is the iterations of
-iterations of level 2 per iteration of level 3, etc.  N-1 looping
+level 2 per iteration of level 3, etc.  N-1 looping parameters must be
-parameters must be specified.
+specified.
 
-Thus with a 4-level respa setting of "2 2 2" for the 3 loop factors,
+Thus with a 4-level respa setting of "2 2 2" for the 3 loop factors, you
-you could choose to have bond interactions computed 8x per large
+could choose to have bond interactions computed 8x per large timestep,
-timestep, angle interactions computed 4x, pair interactions computed
+angle interactions computed 4x, pair interactions computed 2x, and
-2x, and long-range interactions once per large timestep.
+long-range interactions once per large timestep.
 
 The :doc:`timestep <timestep>` command sets the large timestep for the
 outermost rRESPA level.  Thus if the 3 loop factors are "2 2 2" for
-4-level rRESPA, and the outer timestep is set to 4.0 fs, then the
+4-level rRESPA, and the outer timestep is set to 4.0 fs, then the inner
-inner timestep would be 8x smaller or 0.5 fs.  All other LAMMPS
+timestep would be 8x smaller or 0.5 fs.  All other LAMMPS commands that
-commands that specify number of timesteps (e.g. :doc:`thermo <thermo>`
+specify number of timesteps (e.g. :doc:`thermo <thermo>` for thermo
-for thermo output every N steps, :doc:`neigh_modify delay/every <neigh_modify>` parameters, :doc:`dump <dump>` every N
+output every N steps, :doc:`neigh_modify delay/every <neigh_modify>`
-steps, etc) refer to the outermost timesteps.
+parameters, :doc:`dump <dump>` every N steps, etc) refer to the
+outermost timesteps.
 
-The rRESPA keywords enable you to specify at what level of the
+The rRESPA keywords enable you to specify at what level of the hierarchy
-hierarchy various forces will be computed.  If not specified, the
+various forces will be computed.  If not specified, the defaults are
-defaults are that bond forces are computed at level 1 (innermost
+that bond forces are computed at level 1 (innermost loop), angle forces
-loop), angle forces are computed where bond forces are, dihedral
+are computed where bond forces are, dihedral forces are computed where
-forces are computed where angle forces are, improper forces are
+angle forces are, improper forces are computed where dihedral forces
-computed where dihedral forces are, pair forces are computed at the
+are, pair forces are computed at the outermost level, and kspace forces
-outermost level, and kspace forces are computed where pair forces are.
+are computed where pair forces are.  The inner, middle, outer forces
-The inner, middle, outer forces have no defaults.
+have no defaults.
 
 For fixes that support it, the rRESPA level at which a given fix is
-active, can be selected through the :doc:`fix_modify <fix_modify>` command.
+active, can be selected through the :doc:`fix_modify <fix_modify>`
+command.
 
-The *inner* and *middle* keywords take additional arguments for
+The *inner* and *middle* keywords take additional arguments for cutoffs
-cutoffs that are used by the pairwise force computations.  If the 2
+that are used by the pairwise force computations.  If the 2 cutoffs for
-cutoffs for *inner* are 5.0 and 6.0, this means that all pairs up to
+*inner* are 5.0 and 6.0, this means that all pairs up to 6.0 apart are
-6.0 apart are computed by the inner force.  Those between 5.0 and 6.0
+computed by the inner force.  Those between 5.0 and 6.0 have their force
-have their force go ramped to 0.0 so the overlap with the next regime
+go ramped to 0.0 so the overlap with the next regime (middle or outer)
-(middle or outer) is smooth.  The next regime (middle or outer) will
+is smooth.  The next regime (middle or outer) will compute forces for
-compute forces for all pairs from 5.0 outward, with those from 5.0 to
+all pairs from 5.0 outward, with those from 5.0 to 6.0 having their
-6.0 having their value ramped in an inverse manner.
+value ramped in an inverse manner.
 
 Note that you can use *inner* and *outer* without using *middle* to
 split the pairwise computations into two portions instead of three.
-Unless you are using a very long pairwise cutoff, a 2-way split is
+Unless you are using a very long pairwise cutoff, a 2-way split is often
-often faster than a 3-way split, since it avoids too much duplicate
+faster than a 3-way split, since it avoids too much duplicate
 computation of pairwise interactions near the intermediate cutoffs.
 
-Also note that only a few pair potentials support the use of the
+Also note that only a few pair potentials support the use of the *inner*
-*inner* and *middle* and *outer* keywords.  If not, only the *pair*
+and *middle* and *outer* keywords.  If not, only the *pair* keyword can
-keyword can be used with that pair style, meaning all pairwise forces
+be used with that pair style, meaning all pairwise forces are computed
-are computed at the same rRESPA level.  See the doc pages for
+at the same rRESPA level.  See the doc pages for individual pair styles
-individual pair styles for details.
+for details.
 
 Another option for using pair potentials with rRESPA is with the
-*hybrid* keyword, which requires the use of the :doc:`pair_style hybrid or hybrid/overlay <pair_hybrid>` command.  In this scenario, different
+*hybrid* keyword, which requires the use of the :doc:`pair_style hybrid
+or hybrid/overlay <pair_hybrid>` command.  In this scenario, different
 sub-styles of the hybrid pair style are evaluated at different rRESPA
-levels.  This can be useful, for example, to set different timesteps
+levels.  This can be useful, for example, to set different timesteps for
-for hybrid coarse-grained/all-atom models.  The *hybrid* keyword
+hybrid coarse-grained/all-atom models.  The *hybrid* keyword requires as
-requires as many level assignments as there are hybrid sub-styles,
+many level assignments as there are hybrid sub-styles, which assigns
-which assigns each sub-style to a rRESPA level, following their order
+each sub-style to a rRESPA level, following their order of definition in
-of definition in the pair_style command. Since the *hybrid* keyword
+the pair_style command. Since the *hybrid* keyword operates on pair
-operates on pair style computations, it is mutually exclusive with
+style computations, it is mutually exclusive with either the *pair* or
-either the *pair* or the *inner*\ /\ *middle*\ /\ *outer* keywords.
+the *inner*\ /\ *middle*\ /\ *outer* keywords.
 
 When using rRESPA (or for any MD simulation) care must be taken to
 choose a timestep size(s) that ensures the Hamiltonian for the chosen

doc/src/set.rst

@@ -11,7 +11,16 @@ Syntax
    set style ID keyword values ...
 
 * style = *atom* or *type* or *mol* or *group* or *region*
-* ID = atom ID range or type range or mol ID range or group ID or region ID
+* ID = depends on style
+
+.. parsed-literal::
+
+       for style = *atom*, ID = a range of atom IDs
+       for style = *type*, ID = a range of numeric types or a single type label
+       for style = *mol*, ID = a range of molecule IDs
+       for style = *group*, ID = a group ID
+       for style = *region*, ID = a region ID
+
 * one or more keyword/value pairs may be appended
 * keyword = *type* or *type/fraction* or *type/ratio* or *type/subset*
   or *mol* or *x* or *y* or *z* or *vx* or *vy* or *vz* or *charge* or
@@ -19,7 +28,7 @@ Syntax
   *spin/electron* or *radius/electron* or
   *quat* or *quat/random* or *diameter* or *shape* or *length* or *tri* or
   *theta* or *theta/random* or *angmom* or *omega* or
-  *mass* or *density* or *density/disc* or
+  *mass* or *density* or *density/disc* or *temperature* or
   *volume* or *image* or *bond* or *angle* or *dihedral* or
   *improper* or *sph/e* or *sph/cv* or *sph/rho* or
   *smd/contact/radius* or *smd/mass/density* or *dpd/theta* or
@@ -28,22 +37,22 @@ Syntax
 
   .. parsed-literal::
 
-       *type* value = atom type
+       *type* value = numeric atom type or type label
          value can be an atom-style variable (see below)
        *type/fraction* values = type fraction seed
-         type = new atom type
+         type = numeric atom type or type label
          fraction = approximate fraction of selected atoms to set to new atom type
          seed = random # seed (positive integer)
        *type/ratio* values = type fraction seed
-         type = new atom type
+         type = numeric atom type or type label
          fraction = exact fraction of selected atoms to set to new atom type
          seed = random # seed (positive integer)
        *type/subset* values = type Nsubset seed
-         type = new atom type
+         type = numeric atom type or type label
          Nsubset = exact number of selected atoms to set to new atom type
          seed = random # seed (positive integer)
        *mol* value = molecule ID
-         value can be an atom-style variable (see below)
+       value can be an atom-style variable (see below)
        *x*,\ *y*,\ *z* value = atom coordinate (distance units)
          value can be an atom-style variable (see below)
        *vx*,\ *vy*,\ *vz* value = atom velocity (velocity units)
@@ -100,15 +109,17 @@ Syntax
          value can be an atom-style variable (see below)
        *density/disc* value = particle density for a 2d disc or ellipse (mass/distance\^2 units)
          value can be an atom-style variable (see below)
+       *temperature* value = temperature for finite-size particles (temperature units)
+         value can be an atom-style variable (see below)
        *volume* value = particle volume for Peridynamic particle (distance\^3 units)
          value can be an atom-style variable (see below)
        *image* nx ny nz
          nx,ny,nz = which periodic image of the simulation box the atom is in
          any of nx,ny,nz can be an atom-style variable (see below)
-       *bond* value = bond type for all bonds between selected atoms
+       *bond* value = numeric bond type or bond type label, for all bonds between selected atoms
-       *angle* value = angle type for all angles between selected atoms
+       *angle* value = numeric angle type or angle type label, for all angles between selected atoms
-       *dihedral* value = dihedral type for all dihedrals between selected atoms
+       *dihedral* value = numeric dihedral type or dihedral type label, for all dihedrals between selected atoms
-       *improper* value = improper type for all impropers between selected atoms
+       *improper* value = numeric improper type or improper type label, for all impropers between selected atoms
        *sph/e* value = energy of SPH particles (need units)
          value can be an atom-style variable (see below)
        *sph/cv* value = heat capacity of SPH particles (need units)
@@ -143,15 +154,19 @@ Examples
 .. code-block:: LAMMPS
 
    set group solvent type 2
+   set group solvent type C
    set group solvent type/fraction 2 0.5 12393
+   set group solvent type/fraction C 0.5 12393
    set group edge bond 4
    set region half charge 0.5
    set type 3 charge 0.5
+   set type H charge 0.5
    set type 1*3 charge 0.5
    set atom * charge v_atomfile
    set atom 100*200 x 0.5 y 1.0
    set atom 100 vx 0.0 vy 0.0 vz -1.0
    set atom 1492 type 3
+   set atom 1492 type H
    set atom * i_myVal 5
    set atom * d2_Sxyz[1] 6.4
 
@@ -181,9 +196,17 @@ properties to reset and what the new values are.  Some strings like
 This section describes how to select which atoms to change
 the properties of, via the *style* and *ID* arguments.
 
-The style *atom* selects all the atoms in a range of atom IDs.  The
+.. versionchanged:: 28Mar2023
-style *type* selects all the atoms in a range of types.  The style
+
-*mol* selects all the atoms in a range of molecule IDs.
+   Support for type labels was added for selecting atoms by type
+
+The style *atom* selects all the atoms in a range of atom IDs.
+
+The style *type* selects all the atoms in a range of types or type
+labels.  The style *type* selects atoms in one of two ways.  A range
+of numeric atom types can be specified.  Or a single atom type label
+can be specified, e.g. "C".  The style *mol* selects all the atoms in
+a range of molecule IDs.
 
 In each of the range cases, the range can be specified as a single
 numeric value, or a wildcard asterisk can be used to specify a range
@@ -235,14 +258,23 @@ from a file.
    such as the molecule ID, then the floating point value is truncated to
    its integer portion, e.g. a value of 2.6 would become 2.
 
-Keyword *type* sets the atom type for all selected atoms.  The
+.. versionchanged:: 28Mar2023
-specified value must be from 1 to ntypes, where ntypes was set by the
-:doc:`create_box <create_box>` command or the *atom types* field in the
-header of the data file read by the :doc:`read_data <read_data>`
-command.
 
-Keyword *type/fraction* sets the atom type for a fraction of the
+   Support for type labels was added for setting atom, bond, angle,
-selected atoms.  The actual number of atoms changed is not guaranteed
+   dihedral, and improper types
+
+Keyword *type* sets the atom type for all selected atoms.  A specified
+value can be either a numeric atom type or an atom type label. When
+using a numeric type, the specified value must be from 1 to ntypes,
+where ntypes was set by the :doc:`create_box <create_box>` command or
+the *atom types* field in the header of the data file read by the
+:doc:`read_data <read_data>` command.  When using a type label it must
+have been defined previously.  See the :doc:`Howto type labels
+<Howto_type_labels>` doc page for the allowed syntax of type labels
+and a general discussion of how type labels can be used.
+
+Keyword *type/fraction* sets the atom type for a fraction of the selected
+atoms.  The actual number of atoms changed is not guaranteed
 to be exactly the specified fraction (0 <= *fraction* <= 1), but
 should be statistically close.  Random numbers are used in such a way
 that a particular atom is changed or not changed, regardless of how
@@ -429,6 +461,12 @@ assumed to be in mass/distance\^2 units).
 If none of these cases are valid, then the mass is set to the density
 value directly (the input density is assumed to be in mass units).
 
+Keyword *temperature* sets the temperature of a finite-size particle.
+Currently, only the GRANULAR package supports this attribute. The
+temperature must be added using an instance of
+:doc:`fix property/atom <fix_property_atom>` The values for the
+temperature must be positive.
+
 Keyword *volume* sets the volume of all selected particles.  Currently,
 only the :doc:`atom_style peri <atom_style>` command defines particles
 with a volume attribute.  Note that this command does not adjust the
@@ -455,14 +493,18 @@ simulation, but may mess up analysis of the trajectories if a LAMMPS
 diagnostic or your own analysis relies on the image flags to unwrap a
 molecule which straddles the periodic box.
 
-Keywords *bond*, *angle*, *dihedral*, and *improper*, set the bond type
+Keywords *bond*, *angle*, *dihedral*, and *improper*, set the bond
-(angle type, etc) of all bonds (angles, etc) of selected atoms to the
+type (angle type, etc) of all bonds (angles, etc) of selected atoms to
-specified value from 1 to nbondtypes (nangletypes, etc).  All atoms in a
+the specified value.  The value can be a numeric type from 1 to
-particular bond (angle, etc) must be selected atoms in order for the
+nbondtypes (nangletypes, etc).  Or it can be a type label (bond type
-change to be made.  The value of nbondtype (nangletypes, etc) was set by
+label, angle type label, etc).  See the :doc:`Howto type labels
-the *bond types* (\ *angle types*, etc) field in the header of the data
+<Howto_type_labels>` doc page for the allowed syntax of type labels
-file read by the :doc:`read_data <read_data>` command.  These keywords
+and a general discussion of how type labels can be used.  All atoms in
-do not allow use of an atom-style variable.
+a particular bond (angle, etc) must be selected atoms in order for the
+change to be made.  The value of nbondtypes (nangletypes, etc) was set
+by the *bond types* (\ *angle types*, etc) field in the header of the
+data file read by the :doc:`read_data <read_data>` command.  These
+keywords do not allow use of an atom-style variable.
 
 Keywords *sph/e*, *sph/cv*, and *sph/rho* set the energy, heat capacity,
 and density of smoothed particle hydrodynamics (SPH) particles.  See

doc/utils/check-packages.py

@@ -49,7 +49,9 @@ pkgs = []
 # folder, and is not called 'MAKE' is a package 
 
 for d in pkgdirs:
-    pkg = dirs.match(d).group(1)
+    match = dirs.match(d)
+    if not match: continue
+    pkg = match.group(1)
     if not os.path.isdir(os.path.join(src_dir, pkg)): continue
     if pkg in ['DEPEND','MAKE','STUBS']: continue
     pkgs.append(pkg)

doc/utils/requirements.txt

@@ -1,5 +1,6 @@
-Sphinx<6.0.0
+Sphinx < 6.0.0
 sphinxcontrib-spelling
+sphinxcontrib-jquery >=3.0.0
 git+https://github.com/akohlmey/sphinx-fortran@parallel-read
 sphinx_tabs
 breathe

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

@@ -47,6 +47,7 @@ extensions = [
     'sphinx.ext.mathjax',
     'sphinx.ext.imgmath',
     'sphinx.ext.autodoc',
+    'sphinxcontrib.jquery',
     'sphinxfortran.fortran_domain',
     'sphinx_tabs.tabs',
     'table_from_list',
@@ -326,6 +327,7 @@ latex_elements = {
 % Make ToC number fields wider to accommodate sections with >= 100 subsections
 %   or >= 10 subsections with >= 10 subsubsections
 \makeatletter
+\@addtoreset{chapter}{part}
 \renewcommand*{\sphinxtableofcontentshook}{%
   \renewcommand*\l@section{\@dottedtocline{1}{1.5em}{3.1em}}
   \renewcommand*\l@subsection{\@dottedtocline{2}{4.6em}{4.5em}}

doc/utils/sphinx-config/false_positives.txt

@@ -1,5 +1,6 @@
 aa
 aat
+Abascal
 abc
 abf
 ABI
@@ -270,6 +271,7 @@ Bialke
 biaxial
 bicrystal
 Biersack
+biga
 bigbig
 bigint
 Bij
@@ -317,6 +319,7 @@ Bogaerts
 Bogusz
 Bohrs
 boltz
+Bomont
 BondAngle
 BondBond
 bondchk
@@ -355,6 +358,7 @@ br
 Branduardi
 Branicio
 brennan
+Bretonnet
 Briels
 Brien
 Brilliantov
@@ -437,6 +441,7 @@ cfile
 CFL
 cgdna
 CGDNA
+cgmap
 cgs
 cgsdk
 CGSDK
@@ -629,6 +634,7 @@ cutlo
 cutmax
 cutoffA
 cutoffC
+cutsq
 cuu
 cv
 Cv
@@ -1370,6 +1376,7 @@ hdnnp
 HDNNP
 Hearn
 heatconduction
+heatflow
 Hebbeker
 Hebenstreit
 Hecht
@@ -1442,6 +1449,7 @@ hux
 hwloc
 hx
 hy
+hydronium
 hydrophobicity
 hydrostatic
 hydrostatically
@@ -1791,6 +1799,7 @@ lammps
 Lammps
 LAMMPS
 lammpsbin
+lammpsdata
 lammpsplot
 lammpsplugin
 Lamoureux
@@ -1974,6 +1983,7 @@ machdyn
 MACHDYN
 Mackay
 Mackrodt
+MacLaurin
 macOS
 Macromolecules
 macroparticle
@@ -1985,6 +1995,7 @@ magelec
 Maginn
 magneton
 magnetons
+Mahoney
 mainboard
 mainboards
 makefile
@@ -2070,6 +2081,7 @@ mdi
 MDI
 mdpd
 mDPD
+mdtraj
 meam
 MEAM
 meamf
@@ -2206,6 +2218,7 @@ monopole
 monovalent
 Montalenti
 Montero
+Monti
 Mora
 Morefoo
 Morfill
@@ -2463,6 +2476,7 @@ nodeless
 nodeset
 nodesets
 Noehring
+nofix
 Noffset
 noforce
 noguess
@@ -2566,7 +2580,7 @@ nvidia
 nvk
 nvt
 Nwait
-nwchem
+NWChem
 nx
 Nx
 nxlo
@@ -2900,6 +2914,7 @@ PyLammps
 pymbar
 pymodule
 pymol
+PySCF
 pythonic
 pytorch
 pyy
@@ -3186,6 +3201,7 @@ Sandia
 sandybrown
 sanitizer
 Sanyal
+Sanz
 Sarath
 sc
 scafacos
@@ -3243,6 +3259,7 @@ Sep
 seqdep
 Serpico
 setfl
+setflag
 setforce
 Sethna
 setmask
@@ -3478,6 +3495,7 @@ Sx
 sy
 Sy
 symplectic
+symmetrize
 Synechococcus
 sys
 sysdim
@@ -3487,6 +3505,7 @@ sz
 Sz
 Tabbernor
 tabinner
+tabstyle
 Tadmor
 Tafipolsky
 tagID

examples/PACKAGES/alchemy/AlCu.eam.alloy

@@ -0,0 +1 @@
+../../../potentials/AlCu.eam.alloy

examples/PACKAGES/alchemy/h2o.mol

@@ -0,0 +1,62 @@
+# Water molecule. SPC/E model.
+
+3 atoms
+2 bonds
+1 angles
+
+Coords
+
+1    1.12456   0.09298   1.27452
+2    1.53683   0.75606   1.89928
+3    0.49482   0.56390   0.65678
+
+Types
+
+1        1
+2        2
+3        2
+
+Charges
+
+1       -0.8472
+2        0.4236
+3        0.4236
+
+Bonds
+
+1   1      1      2
+2   1      1      3
+
+Angles
+
+1   1      2      1      3
+
+Shake Flags
+
+1 1
+2 1
+3 1
+
+Shake Atoms
+
+1 1 2 3
+2 1 2 3
+3 1 2 3
+
+Shake Bond Types
+
+1 1 1 1
+2 1 1 1
+3 1 1 1
+
+Special Bond Counts
+
+1 2 0 0
+2 1 1 0
+3 1 1 0
+
+Special Bonds
+
+1 2 3
+2 1 3
+3 1 2