Merge branch 'develop' into mliappy_unified

.github/CONTRIBUTING.md

@@ -5,8 +5,8 @@ Thank your for considering to contribute to the LAMMPS software project.
 The following is a set of guidelines as well as explanations of policies and work flows for contributing to the LAMMPS molecular dynamics software project. These guidelines focus on submitting issues or pull requests on the LAMMPS GitHub project.
 
 Thus please also have a look at:
-* [The guide for submitting new features in the LAMMPS manual](https://lammps.sandia.gov/doc/Modify_contribute.html)
-* [The guide on programming style and requirement in the LAMMPS manual](https://lammps.sandia.gov/doc/Modify_contribute.html)
+* [The guide for submitting new features in the LAMMPS manual](https://www.lammps.org/doc/Modify_contribute.html)
+* [The guide on programming style and requirement in the LAMMPS manual](https://www.lammps.org/doc/Modify_style.html)
 * [The GitHub tutorial in the LAMMPS manual](http://lammps.sandia.gov/doc/Howto_github.html)
 
 ## Table of Contents

.github/workflows/compile-msvc.yml

@@ -26,6 +26,7 @@ jobs:
       shell: bash
       run: |
         python3 -m pip install numpy
+        python3 -m pip install pyyaml
         cmake -C cmake/presets/windows.cmake \
               -D PKG_PYTHON=on \
               -S cmake -B build \

.github/workflows/unittest-macos.yml

@@ -37,6 +37,7 @@ jobs:
       working-directory: build
       run: |
         ccache -z
+        python3 -m pip install pyyaml
         cmake -C ../cmake/presets/clang.cmake \
               -C ../cmake/presets/most.cmake \
               -D CMAKE_CXX_COMPILER_LAUNCHER=ccache \

bench/README

@@ -4,10 +4,8 @@ This directory contains 5 benchmark problems which are discussed in
 the Benchmark section of the LAMMPS documentation, and on the
 Benchmark page of the LAMMPS WWW site (https://www.lammps.org/bench.html).
 
-This directory also has several sub-directories:
+This directory also has one sub-directories:
 
-FERMI           benchmark scripts for desktop machine with Fermi GPUs (Tesla)
-KEPLER          benchmark scripts for GPU cluster with Kepler GPUs
 POTENTIALS      benchmarks scripts for various potentials in LAMMPS
 
 The results for all of these benchmarks are displayed and discussed on

cmake/CMakeLists.txt

@@ -105,8 +105,28 @@ if(CMAKE_CXX_COMPILER_ID STREQUAL "Intel")
   endif()
 endif()
 
+# silence excessive warnings for new Intel Compilers
+if(CMAKE_CXX_COMPILER_ID STREQUAL "IntelLLVM")
+  set(CMAKE_TUNE_DEFAULT "-Wno-tautological-constant-compare -Wno-unused-command-line-argument")
+endif()
+
+# silence excessive warnings for PGI/NVHPC compilers
+if((CMAKE_CXX_COMPILER_ID STREQUAL "NVHPC") OR (CMAKE_CXX_COMPILER_ID STREQUAL "PGI"))
+  set(CMAKE_TUNE_DEFAULT "-Minform=severe")
+endif()
+
+# silence nvcc warnings
+if((PKG_KOKKOS) AND (Kokkos_ENABLE_CUDA) AND NOT (CMAKE_CXX_COMPILER_ID STREQUAL "Clang"))
+  set(CMAKE_TUNE_DEFAULT "${CMAKE_TUNE_DEFAULT} -Xcudafe --diag_suppress=unrecognized_pragma")
+endif()
+
 # we require C++11 without extensions. Kokkos requires at least C++14 (currently)
-set(CMAKE_CXX_STANDARD 11)
+if(NOT CMAKE_CXX_STANDARD)
+  set(CMAKE_CXX_STANDARD 11)
+endif()
+if(CMAKE_CXX_STANDARD LESS 11)
+  message(FATAL_ERROR "C++ standard must be set to at least 11")
+endif()
 if(PKG_KOKKOS AND (CMAKE_CXX_STANDARD LESS 14))
   set(CMAKE_CXX_STANDARD 14)
 endif()
@@ -141,10 +161,12 @@ option(CMAKE_POSITION_INDEPENDENT_CODE "Create object compatible with shared lib
 option(BUILD_TOOLS "Build and install LAMMPS tools (msi2lmp, binary2txt, chain)" OFF)
 option(BUILD_LAMMPS_SHELL "Build and install the LAMMPS shell" OFF)
 
-# allow enabling clang-tidy for C++ files
+# Support using clang-tidy for C++ files with selected options
 set(ENABLE_CLANG_TIDY OFF CACHE BOOL "Include clang-tidy processing when compiling")
 if(ENABLE_CLANG_TIDY)
-  set(CMAKE_CXX_CLANG_TIDY "clang-tidy;-checks=*-header-filter=.*" CACHE STRING "")
+  set(CMAKE_CXX_CLANG_TIDY "clang-tidy;-checks=-*,performance-trivially-destructible,performance-unnecessary-copy-initialization,performance-unnecessary-value-param,readability-redundant-control-flow,readability-redundant-declaration,readability-redundant-function-ptr-dereference,readability-redundant-member-init,readability-redundant-string-cstr,readability-redundant-string-init,readability-simplify-boolean-expr,readability-static-accessed-through-instance,readability-static-definition-in-anonymous-namespace,modernize-use-override,modernize-use-bool-literals,modernize-use-emplace,modernize-return-braced-init-list,modernize-use-equals-default,modernize-use-equals-delete,modernize-replace-random-shuffle,modernize-deprecated-headers,modernize-use-nullptr,modernize-use-noexcept,modernize-redundant-void-arg;-fix;-header-filter=.*,header-filter=library.h,header-filter=fmt/*.h" CACHE STRING "clang-tidy settings")
+else()
+  unset(CMAKE_CXX_CLANG_TIDY CACHE)
 endif()
 
 include(GNUInstallDirs)
@@ -172,6 +194,7 @@ set(STANDARD_PACKAGES
   AWPMD
   BOCS
   BODY
+  BPM
   BROWNIAN
   CG-DNA
   CG-SDK
@@ -188,6 +211,7 @@ set(STANDARD_PACKAGES
   DPD-SMOOTH
   DRUDE
   EFF
+  ELECTRODE
   EXTRA-COMPUTE
   EXTRA-DUMP
   EXTRA-FIX
@@ -208,7 +232,6 @@ set(STANDARD_PACKAGES
   MDI
   MEAM
   MESONT
-  MESSAGE
   MGPT
   MISC
   ML-HDNNP
@@ -307,7 +330,9 @@ string(TOUPPER ${LAMMPS_SIZES} LAMMPS_SIZES)
 target_compile_definitions(lammps PUBLIC -DLAMMPS_${LAMMPS_SIZES})
 
 # posix_memalign is not available on Windows
-if(CMAKE_SYSTEM_NAME STREQUAL "Windows")
+# with INTEL package and Intel compilers we use TBB's aligned malloc
+if((CMAKE_SYSTEM_NAME STREQUAL "Windows")
+   AND NOT (PKG_INTEL AND ((CMAKE_CXX_COMPILER_ID STREQUAL "Intel") OR (CMAKE_CXX_COMPILER_ID STREQUAL "IntelLLVM"))))
   set(LAMMPS_MEMALIGN "0" CACHE STRING "posix_memalign() is not available on Windows" FORCE)
 else()
   set(LAMMPS_MEMALIGN "64" CACHE STRING "enables the use of the posix_memalign() call instead of malloc() when large chunks or memory are allocated by LAMMPS. Set to 0 to disable")
@@ -332,6 +357,7 @@ pkg_depends(DIELECTRIC KSPACE)
 pkg_depends(DIELECTRIC EXTRA-PAIR)
 pkg_depends(CG-DNA MOLECULE)
 pkg_depends(CG-DNA ASPHERE)
+pkg_depends(ELECTRODE KSPACE)
 
 # detect if we may enable OpenMP support by default
 set(BUILD_OMP_DEFAULT OFF)
@@ -369,7 +395,7 @@ if(BUILD_OMP)
   target_link_libraries(lmp PRIVATE OpenMP::OpenMP_CXX)
 endif()
 
-if(PKG_MSCG OR PKG_ATC OR PKG_AWPMD OR PKG_ML-QUIP OR PKG_LATTE)
+if(PKG_MSCG OR PKG_ATC OR PKG_AWPMD OR PKG_ML-QUIP OR PKG_LATTE OR PKG_ELECTRODE)
   enable_language(C)
   find_package(LAPACK)
   find_package(BLAS)
@@ -451,7 +477,7 @@ else()
 endif()
 
 foreach(PKG_WITH_INCL KSPACE PYTHON ML-IAP VORONOI COLVARS ML-HDNNP MDI MOLFILE NETCDF
-        PLUMED QMMM ML-QUIP SCAFACOS MACHDYN VTK KIM LATTE MESSAGE MSCG COMPRESS ML-PACE)
+        PLUMED QMMM ML-QUIP SCAFACOS MACHDYN VTK KIM LATTE MSCG COMPRESS ML-PACE)
   if(PKG_${PKG_WITH_INCL})
     include(Packages/${PKG_WITH_INCL})
   endif()
@@ -468,6 +494,7 @@ 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})
 endforeach()
 ########################################################################
 # Basic system tests (standard libraries, headers, functions, types)   #
@@ -572,6 +599,10 @@ foreach(PKG_LIB POEMS ATC AWPMD H5MD MESONT)
   endif()
 endforeach()
 
+if(PKG_ELECTRODE)
+  target_link_libraries(lammps PRIVATE ${LAPACK_LIBRARIES})
+endif()
+
 if(PKG_AWPMD)
   target_link_libraries(awpmd PRIVATE ${LAPACK_LIBRARIES})
 endif()
@@ -810,6 +841,7 @@ endif()
 include(FeatureSummary)
 feature_summary(DESCRIPTION "The following tools and libraries have been found and configured:" WHAT PACKAGES_FOUND)
 message(STATUS "<<< Build configuration >>>
+   LAMMPS Version:   ${PROJECT_VERSION}
    Operating System: ${CMAKE_SYSTEM_NAME} ${CMAKE_LINUX_DISTRO} ${CMAKE_DISTRO_VERSION}
    Build type:       ${LAMMPS_BUILD_TYPE}
    Install path:     ${CMAKE_INSTALL_PREFIX}
@@ -919,6 +951,12 @@ if(PKG_KSPACE)
       else()
         message(STATUS "Kokkos FFT: cuFFT")
       endif()
+    elseif(Kokkos_ENABLE_HIP)
+      if(FFT STREQUAL "KISS")
+        message(STATUS "Kokkos FFT: KISS")
+      else()
+        message(STATUS "Kokkos FFT: hipFFT")
+      endif()
     else()
       message(STATUS "Kokkos FFT: ${FFT}")
     endif()

cmake/CMakeSettings.json

@@ -151,6 +151,190 @@
           "type": "BOOL"
         }
       ]
+    },
+    {
+      "name": "x64-Debug-IntelLLVM",
+      "generator": "Ninja",
+      "configurationType": "Debug",
+      "buildRoot": "${workspaceRoot}\\build\\${name}",
+      "installRoot": "${workspaceRoot}\\install\\${name}",
+      "inheritEnvironments": [],
+      "cmakeCommandArgs": "-C ${workspaceRoot}\\cmake\\presets\\windows-intel-llvm.cmake -C ${workspaceRoot}\\cmake\\presets\\windows.cmake",
+      "variables": [
+        {
+          "name": "BUILD_SHARED_LIBS",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "BUILD_TOOLS",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "LAMMPS_EXCEPTIONS",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "PKG_PYTHON",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "ENABLE_TESTING",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "FFT",
+          "value": "MKL",
+          "type": "STRING"
+        },
+        {
+          "name": "PKG_PLUGIN",
+          "value": "True",
+          "type": "BOOL"
+        }
+      ]
+    },
+    {
+      "name": "x64-Release-IntelLLVM",
+      "generator": "Ninja",
+      "configurationType": "Release",
+      "buildRoot": "${workspaceRoot}\\build\\${name}",
+      "installRoot": "${workspaceRoot}\\install\\${name}",
+      "cmakeCommandArgs": "-C ${workspaceRoot}\\cmake\\presets\\windows-intel-llvm.cmake -C ${workspaceRoot}\\cmake\\presets\\windows.cmake",
+      "inheritEnvironments": [],
+      "variables": [
+        {
+          "name": "BUILD_SHARED_LIBS",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "BUILD_TOOLS",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "LAMMPS_EXCEPTIONS",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "PKG_PYTHON",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "ENABLE_TESTING",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "FFT",
+          "value": "MKL",
+          "type": "STRING"
+        },
+        {
+          "name": "PKG_PLUGIN",
+          "value": "True",
+          "type": "BOOL"
+        }
+      ]
+    },
+    {
+      "name": "x64-Debug-Intel-Classic",
+      "generator": "Ninja",
+      "configurationType": "Debug",
+      "buildRoot": "${workspaceRoot}\\build\\${name}",
+      "installRoot": "${workspaceRoot}\\install\\${name}",
+      "cmakeCommandArgs": "-C ${workspaceRoot}\\cmake\\presets\\windows-intel-classic.cmake -C ${workspaceRoot}\\cmake\\presets\\windows.cmake",
+      "inheritEnvironments": [],
+      "variables": [
+        {
+          "name": "BUILD_SHARED_LIBS",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "BUILD_TOOLS",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "LAMMPS_EXCEPTIONS",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "PKG_PYTHON",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "ENABLE_TESTING",
+          "value": "False",
+          "type": "BOOL"
+        },
+        {
+          "name": "FFT",
+          "value": "MKL",
+          "type": "STRING"
+        },
+        {
+          "name": "PKG_PLUGIN",
+          "value": "True",
+          "type": "BOOL"
+        }
+      ]
+    },
+    {
+      "name": "x64-Release-Intel-Classic",
+      "generator": "Ninja",
+      "configurationType": "Release",
+      "buildRoot": "${workspaceRoot}\\build\\${name}",
+      "installRoot": "${workspaceRoot}\\install\\${name}",
+      "cmakeCommandArgs": "-C ${workspaceRoot}\\cmake\\presets\\windows-intel-classic.cmake -C ${workspaceRoot}\\cmake\\presets\\windows.cmake",
+      "inheritEnvironments": [],
+      "variables": [
+        {
+          "name": "BUILD_SHARED_LIBS",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "BUILD_TOOLS",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "LAMMPS_EXCEPTIONS",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "PKG_PYTHON",
+          "value": "True",
+          "type": "BOOL"
+        },
+        {
+          "name": "ENABLE_TESTING",
+          "value": "False",
+          "type": "BOOL"
+        },
+        {
+          "name": "FFT",
+          "value": "MKL",
+          "type": "STRING"
+        },
+        {
+          "name": "PKG_PLUGIN",
+          "value": "True",
+          "type": "BOOL"
+        }
+      ]
     }
   ]
 }

cmake/Modules/Documentation.cmake

@@ -7,13 +7,13 @@ if(BUILD_DOC)
   # Sphinx 3.x requires at least Python 3.5
   if(CMAKE_VERSION VERSION_LESS 3.12)
     find_package(PythonInterp 3.5 REQUIRED)
-    set(VIRTUALENV ${PYTHON_EXECUTABLE} -m virtualenv -p ${PYTHON_EXECUTABLE})
+    set(VIRTUALENV ${PYTHON_EXECUTABLE} -m venv)
   else()
     find_package(Python3 REQUIRED COMPONENTS Interpreter)
     if(Python3_VERSION VERSION_LESS 3.5)
       message(FATAL_ERROR "Python 3.5 and up is required to build the HTML documentation")
     endif()
-    set(VIRTUALENV ${Python3_EXECUTABLE} -m virtualenv -p ${Python3_EXECUTABLE})
+    set(VIRTUALENV ${Python3_EXECUTABLE} -m venv)
   endif()
   find_package(Doxygen 1.8.10 REQUIRED)
 

cmake/Modules/Packages/COLVARS.cmake

@@ -19,6 +19,10 @@ endif()
 
 add_library(colvars STATIC ${COLVARS_SOURCES})
 target_compile_definitions(colvars PRIVATE -DCOLVARS_LAMMPS)
+separate_arguments(CMAKE_TUNE_FLAGS)
+foreach(_FLAG ${CMAKE_TUNE_FLAGS})
+  target_compile_options(colvars PRIVATE ${_FLAG})
+endforeach()
 set_target_properties(colvars PROPERTIES OUTPUT_NAME lammps_colvars${LAMMPS_MACHINE})
 target_include_directories(colvars PUBLIC ${LAMMPS_LIB_SOURCE_DIR}/colvars)
 # The line below is needed to locate math_eigen_impl.h

cmake/Modules/Packages/GPU.cmake

@@ -80,14 +80,14 @@ if(GPU_API STREQUAL "CUDA")
   # --arch translates directly instead of JIT, so this should be for the preferred or most common architecture
   set(GPU_CUDA_GENCODE "-arch=${GPU_ARCH}")
 
-  # apply the following to build "fat" CUDA binaries only for known CUDA toolkits
-  if(CUDA_VERSION VERSION_GREATER_EQUAL "12.0")
-    message(WARNING "Untested CUDA Toolkit version. Use at your own risk")
+  # apply the following to build "fat" CUDA binaries only for known CUDA toolkits since version 8.0
+  # only the Kepler achitecture and beyond is 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")
+    message(WARNING "Untested CUDA Toolkit version ${CUDA_VERSION}. Use at your own risk")
+    set(GPU_CUDA_GENCODE "-arch=all")
   else()
-    # Fermi (GPU Arch 2.x) is supported by CUDA 3.2 to CUDA 8.0
-    if((CUDA_VERSION VERSION_GREATER_EQUAL "3.2") AND (CUDA_VERSION VERSION_LESS "9.0"))
-      string(APPEND GPU_CUDA_GENCODE " -gencode arch=compute_20,code=[sm_20,compute_20] ")
-    endif()
     # 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"))
       string(APPEND GPU_CUDA_GENCODE " -gencode arch=compute_30,code=[sm_30,compute_30] ")
@@ -120,6 +120,14 @@ if(GPU_API STREQUAL "CUDA")
     if(CUDA_VERSION VERSION_GREATER_EQUAL "11.1")
       string(APPEND GPU_CUDA_GENCODE " -gencode arch=compute_86,code=[sm_86,compute_86]")
     endif()
+    # Hopper (GPU Arch 9.0) is supported by CUDA 12.0? and later
+    if(CUDA_VERSION VERSION_GREATER_EQUAL "12.0")
+      string(APPEND GPU_CUDA_GENCODE " -gencode arch=compute_90,code=[sm_90,compute_90]")
+    endif()
+    #    # Lovelace (GPU Arch 9.x) is supported by CUDA 12.0? and later
+    #if(CUDA_VERSION VERSION_GREATER_EQUAL "12.0")
+    #  string(APPEND GPU_CUDA_GENCODE " -gencode arch=compute_9x,code=[sm_9x,compute_9x]")
+    #endif()
   endif()
 
   cuda_compile_fatbin(GPU_GEN_OBJS ${GPU_LIB_CU} OPTIONS ${CUDA_REQUEST_PIC}
@@ -257,13 +265,15 @@ elseif(GPU_API STREQUAL "HIP")
     find_package(CUDA REQUIRED)
     set(HIP_ARCH "sm_50" CACHE STRING "HIP primary CUDA architecture (e.g. sm_60)")
 
+    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")
+      message(WARNING "Untested CUDA Toolkit version ${CUDA_VERSION}. Use at your own risk")
+      set(HIP_CUDA_GENCODE "-arch=all")
+    else()
       # build arch/gencode commands for nvcc based on CUDA toolkit version and use choice
       # --arch translates directly instead of JIT, so this should be for the preferred or most common architecture
       set(HIP_CUDA_GENCODE "-arch=${HIP_ARCH}")
-    # Fermi (GPU Arch 2.x) is supported by CUDA 3.2 to CUDA 8.0
-    if((CUDA_VERSION VERSION_GREATER_EQUAL "3.2") AND (CUDA_VERSION VERSION_LESS "9.0"))
-      string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_20,code=[sm_20,compute_20]")
-    endif()
       # 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"))
         string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_30,code=[sm_30,compute_30]")
@@ -292,8 +302,10 @@ elseif(GPU_API STREQUAL "HIP")
       if(CUDA_VERSION VERSION_GREATER_EQUAL "11.0")
         string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_80,code=[sm_80,compute_80]")
       endif()
+      # Hopper (GPU Arch 9.0) is supported by CUDA 12.0? and later
       if(CUDA_VERSION VERSION_GREATER_EQUAL "12.0")
-      message(WARNING "Unsupported CUDA version. Use at your own risk.")
+        string(APPEND GPU_CUDA_GENCODE " -gencode arch=compute_90,code=[sm_90,compute_90]")
+      endif()
     endif()
   endif()
 
@@ -347,6 +359,10 @@ elseif(GPU_API STREQUAL "HIP")
   target_link_libraries(gpu PRIVATE hip::host)
 
   if(HIP_USE_DEVICE_SORT)
+    if(HIP_PLATFORM STREQUAL "amd")
+      # newer version of ROCm (5.1+) require c++14 for rocprim
+      set_property(TARGET gpu PROPERTY CXX_STANDARD 14)
+    endif()
     # add hipCUB
     target_include_directories(gpu PRIVATE ${HIP_ROOT_DIR}/../include)
     target_compile_definitions(gpu PRIVATE -DUSE_HIP_DEVICE_SORT)

cmake/Modules/Packages/H5MD.cmake

@@ -1,5 +1,14 @@
 enable_language(C)
 
+# we don't use the parallel i/o interface.
+set(HDF5_PREFER_PARALLEL FALSE)
+
 find_package(HDF5 REQUIRED)
+
+# parallel HDF5 will import incompatible MPI headers with a serial build
+if((NOT BUILD_MPI) AND HDF5_IS_PARALLEL)
+  message(FATAL_ERROR "Serial LAMMPS build and parallel HDF5 library are not compatible")
+endif()
+
 target_link_libraries(h5md PRIVATE ${HDF5_LIBRARIES})
 target_include_directories(h5md PUBLIC ${HDF5_INCLUDE_DIRS})

cmake/Modules/Packages/INTEL.cmake

@@ -38,7 +38,7 @@ if(INTEL_LRT_MODE STREQUAL "C++11")
   endif()
 endif()
 
-if(CMAKE_CXX_COMPILER_ID STREQUAL "Intel")
+if((CMAKE_CXX_COMPILER_ID STREQUAL "Intel") OR (CMAKE_CXX_COMPILER_ID STREQUAL "IntelLLVM"))
   if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS 16)
     message(FATAL_ERROR "INTEL needs at least a 2016 Intel compiler, found ${CMAKE_CXX_COMPILER_VERSION}")
   endif()
@@ -46,12 +46,12 @@ else()
   message(WARNING "INTEL gives best performance with Intel compilers")
 endif()
 
-find_package(TBB_MALLOC QUIET)
+find_package(TBB_MALLOC)
 if(TBB_MALLOC_FOUND)
   target_link_libraries(lammps PRIVATE TBB::TBB_MALLOC)
 else()
   target_compile_definitions(lammps PRIVATE -DLMP_INTEL_NO_TBB)
-  if(CMAKE_CXX_COMPILER_ID STREQUAL "Intel")
+  if((CMAKE_CXX_COMPILER_ID STREQUAL "Intel") OR (CMAKE_CXX_COMPILER_ID STREQUAL "IntelLLVM"))
     message(WARNING "INTEL with Intel compilers should use TBB malloc libraries")
   endif()
 endif()
@@ -112,5 +112,9 @@ if(PKG_KSPACE)
   RegisterIntegrateStyle(${INTEL_SOURCES_DIR}/verlet_lrt_intel.h)
 endif()
 
+if(PKG_ELECTRODE)
+  list(APPEND INTEL_SOURCES ${INTEL_SOURCES_DIR}/electrode_accel_intel.cpp)
+endif()
+
 target_sources(lammps PRIVATE ${INTEL_SOURCES})
 target_include_directories(lammps PRIVATE ${INTEL_SOURCES_DIR})

cmake/Modules/Packages/KOKKOS.cmake

@@ -47,8 +47,8 @@ if(DOWNLOAD_KOKKOS)
   list(APPEND KOKKOS_LIB_BUILD_ARGS "-DCMAKE_CXX_EXTENSIONS=${CMAKE_CXX_EXTENSIONS}")
   list(APPEND KOKKOS_LIB_BUILD_ARGS "-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TOOLCHAIN_FILE}")
   include(ExternalProject)
-  set(KOKKOS_URL "https://github.com/kokkos/kokkos/archive/3.5.00.tar.gz" CACHE STRING "URL for KOKKOS tarball")
-  set(KOKKOS_MD5 "079323d973ae0e1c38c0a54a150c674e" CACHE STRING "MD5 checksum of KOKKOS tarball")
+  set(KOKKOS_URL "https://github.com/kokkos/kokkos/archive/3.6.00.tar.gz" CACHE STRING "URL for KOKKOS tarball")
+  set(KOKKOS_MD5 "b5c44ea961031795f434002cd7b31c20" CACHE STRING "MD5 checksum of KOKKOS tarball")
   mark_as_advanced(KOKKOS_URL)
   mark_as_advanced(KOKKOS_MD5)
   ExternalProject_Add(kokkos_build
@@ -72,7 +72,7 @@ if(DOWNLOAD_KOKKOS)
   add_dependencies(LAMMPS::KOKKOSCORE kokkos_build)
   add_dependencies(LAMMPS::KOKKOSCONTAINERS kokkos_build)
 elseif(EXTERNAL_KOKKOS)
-  find_package(Kokkos 3.5.00 REQUIRED CONFIG)
+  find_package(Kokkos 3.6.00 REQUIRED CONFIG)
   target_link_libraries(lammps PRIVATE Kokkos::kokkos)
   target_link_libraries(lmp PRIVATE Kokkos::kokkos)
 else()
@@ -130,6 +130,11 @@ if(PKG_KSPACE)
       target_compile_definitions(lammps PRIVATE -DFFT_CUFFT)
       target_link_libraries(lammps PRIVATE cufft)
     endif()
+  elseif(Kokkos_ENABLE_HIP)
+    if(NOT (FFT STREQUAL "KISS"))
+      target_compile_definitions(lammps PRIVATE -DFFT_HIPFFT)
+      target_link_libraries(lammps PRIVATE hipfft)
+    endif()
   endif()
 endif()
 

cmake/Modules/Packages/MACHDYN.cmake

@@ -36,3 +36,5 @@ endif()
 if((CMAKE_CXX_COMPILER_ID STREQUAL "PGI") OR (CMAKE_CXX_COMPILER_ID STREQUAL "NVHPC"))
   target_compile_definitions(lammps PRIVATE -DEIGEN_DONT_VECTORIZE)
 endif()
+
+target_compile_definitions(lammps PRIVATE -DEIGEN_NO_CUDA)

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.2.9.tar.gz" CACHE STRING "URL for MDI tarball")
-  set(MDI_MD5 "ddfa46d6ee15b4e59cfd527ec7212184" CACHE STRING "MD5 checksum for MDI tarball")
+  set(MDI_URL "https://github.com/MolSSI-MDI/MDI_Library/archive/v1.3.2.tar.gz" CACHE STRING "URL for MDI tarball")
+  set(MDI_MD5 "836f5da400d8cff0f0e4435640f9454f" CACHE STRING "MD5 checksum for MDI tarball")
   mark_as_advanced(MDI_URL)
   mark_as_advanced(MDI_MD5)
   enable_language(C)
@@ -50,6 +50,7 @@ if(DOWNLOAD_MDI)
     -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}
     -DCMAKE_MAKE_PROGRAM=${CMAKE_MAKE_PROGRAM}
     -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TOOLCHAIN_FILE}
+    -DCMAKE_POSITION_INDEPENDENT_CODE=ON
     -Dlanguage=C
     -Dlibtype=STATIC
     -Dmpi=${MDI_USE_MPI}

cmake/Modules/Packages/MESSAGE.cmake

@@ -1,31 +0,0 @@
-if(LAMMPS_SIZES STREQUAL "BIGBIG")
-  message(FATAL_ERROR "The MESSAGE Package is not compatible with -DLAMMPS_BIGBIG")
-endif()
-option(MESSAGE_ZMQ "Use ZeroMQ in MESSAGE package" OFF)
-file(GLOB_RECURSE cslib_SOURCES
-        ${LAMMPS_LIB_SOURCE_DIR}/message/cslib/[^.]*.cpp)
-
-add_library(cslib STATIC ${cslib_SOURCES})
-target_compile_definitions(cslib PRIVATE -DLAMMPS_${LAMMPS_SIZES})
-set_target_properties(cslib PROPERTIES OUTPUT_NAME lammps_cslib${LAMMPS_MACHINE})
-if(BUILD_MPI)
-  target_compile_definitions(cslib PRIVATE -DMPI_YES)
-  set_target_properties(cslib PROPERTIES OUTPUT_NAME "csmpi")
-  target_link_libraries(cslib PRIVATE MPI::MPI_CXX)
-else()
-  target_compile_definitions(cslib PRIVATE -DMPI_NO)
-  target_include_directories(cslib PRIVATE ${LAMMPS_LIB_SOURCE_DIR}/message/cslib/src/STUBS_MPI)
-  set_target_properties(cslib PROPERTIES OUTPUT_NAME "csnompi")
-endif()
-
-if(MESSAGE_ZMQ)
-  target_compile_definitions(cslib PRIVATE -DZMQ_YES)
-  find_package(ZMQ REQUIRED)
-  target_link_libraries(cslib PUBLIC ZMQ::ZMQ)
-else()
-  target_compile_definitions(cslib PRIVATE -DZMQ_NO)
-  target_include_directories(cslib PRIVATE ${LAMMPS_LIB_SOURCE_DIR}/message/cslib/src/STUBS_ZMQ)
-endif()
-
-target_link_libraries(lammps PRIVATE cslib)
-target_include_directories(lammps PRIVATE ${LAMMPS_LIB_SOURCE_DIR}/message/cslib/src)

cmake/Modules/Packages/ML-HDNNP.cmake

@@ -42,16 +42,10 @@ if(DOWNLOAD_N2P2)
   if(NOT BUILD_MPI)
     set(N2P2_PROJECT_OPTIONS "-DN2P2_NO_MPI")
   else()
-    # get path to MPI include directory when cross-compiling to windows
-    if((CMAKE_SYSTEM_NAME STREQUAL Windows) AND CMAKE_CROSSCOMPILING)
+    # get path to MPI include directory
     get_target_property(N2P2_MPI_INCLUDE MPI::MPI_CXX INTERFACE_INCLUDE_DIRECTORIES)
     set(N2P2_PROJECT_OPTIONS "-I${N2P2_MPI_INCLUDE}")
   endif()
-    if(CMAKE_CXX_COMPILER_ID STREQUAL "Intel")
-      get_target_property(N2P2_MPI_INCLUDE MPI::MPI_CXX INTERFACE_INCLUDE_DIRECTORIES)
-      set(N2P2_PROJECT_OPTIONS "-I${N2P2_MPI_INCLUDE}")
-    endif()
-  endif()
 
   # prefer GNU make, if available. N2P2 lib seems to need it.
   find_program(N2P2_MAKE NAMES gmake make)

cmake/Modules/Packages/ML-PACE.cmake

@@ -1,6 +1,6 @@
-set(PACELIB_URL "https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/v.2021.10.25.tar.gz" CACHE STRING "URL for PACE evaluator library sources")
+set(PACELIB_URL "https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/v.2021.10.25.fix.tar.gz" CACHE STRING "URL for PACE evaluator library sources")
 
-set(PACELIB_MD5 "a2ac3315c41a1a4a5c912bcb1bc9c5cc" CACHE STRING "MD5 checksum of PACE evaluator library tarball")
+set(PACELIB_MD5 "e0572de57039d4afedefb25707b6ceae" CACHE STRING "MD5 checksum of PACE evaluator library tarball")
 mark_as_advanced(PACELIB_URL)
 mark_as_advanced(PACELIB_MD5)
 
@@ -15,6 +15,10 @@ execute_process(
 )
 
 file(GLOB lib-pace ${CMAKE_BINARY_DIR}/lammps-user-pace-*)
+# enforce building libyaml-cpp as static library and turn off optional features
+set(YAML_BUILD_SHARED_LIBS OFF)
+set(YAML_CPP_BUILD_CONTRIB OFF)
+set(YAML_CPP_BUILD_TOOLS OFF)
 add_subdirectory(${lib-pace}/yaml-cpp build-yaml-cpp)
 set(YAML_CPP_INCLUDE_DIR ${lib-pace}/yaml-cpp/include)
 

cmake/Modules/Tools.cmake

@@ -3,6 +3,9 @@ if(BUILD_TOOLS)
   target_compile_definitions(binary2txt PRIVATE -DLAMMPS_${LAMMPS_SIZES})
   install(TARGETS binary2txt DESTINATION ${CMAKE_INSTALL_BINDIR})
 
+  add_executable(stl_bin2txt ${LAMMPS_TOOLS_DIR}/stl_bin2txt.cpp)
+  install(TARGETS stl_bin2txt DESTINATION ${CMAKE_INSTALL_BINDIR})
+
   include(CheckGeneratorSupport)
   if(CMAKE_GENERATOR_SUPPORT_FORTRAN)
     include(CheckLanguage)

cmake/presets/all_off.cmake

@@ -8,6 +8,7 @@ set(ALL_PACKAGES
   AWPMD
   BOCS
   BODY
+  BPM
   BROWNIAN
   CG-DNA
   CG-SDK
@@ -24,6 +25,7 @@ set(ALL_PACKAGES
   DPD-REACT
   DPD-SMOOTH
   DRUDE
+  ELECTRODE
   EFF
   EXTRA-COMPUTE
   EXTRA-DUMP
@@ -48,7 +50,6 @@ set(ALL_PACKAGES
   MDI
   MEAM
   MESONT
-  MESSAGE
   MGPT
   MISC
   ML-HDNNP

cmake/presets/all_on.cmake

@@ -10,6 +10,7 @@ set(ALL_PACKAGES
   AWPMD
   BOCS
   BODY
+  BPM
   BROWNIAN
   CG-DNA
   CG-SDK
@@ -26,6 +27,7 @@ set(ALL_PACKAGES
   DPD-REACT
   DPD-SMOOTH
   DRUDE
+  ELECTRODE
   EFF
   EXTRA-COMPUTE
   EXTRA-DUMP
@@ -50,7 +52,6 @@ set(ALL_PACKAGES
   MDI
   MEAM
   MESONT
-  MESSAGE
   MGPT
   MISC
   ML-HDNNP

cmake/presets/kokkos-cuda.cmake

@@ -3,9 +3,9 @@
 # that is compatible with all higher CC, but not the default CC 3.5
 set(PKG_KOKKOS ON CACHE BOOL "" FORCE)
 set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "" FORCE)
-set(Kokkos_ENABLE_OPENMP ON CACHE BOOL "" FORCE)
 set(Kokkos_ENABLE_CUDA   ON CACHE BOOL "" FORCE)
-set(Kokkos_ARCH_MAXWELL50 on CACHE BOOL "" FORCE)
+set(Kokkos_ARCH_PASCAL60 ON CACHE BOOL "" FORCE)
 set(BUILD_OMP ON CACHE BOOL "" FORCE)
-get_filename_component(NVCC_WRAPPER_CMD ${CMAKE_CURRENT_SOURCE_DIR}/../lib/kokkos/bin/nvcc_wrapper ABSOLUTE)
-set(CMAKE_CXX_COMPILER ${NVCC_WRAPPER_CMD} CACHE FILEPATH "" FORCE)
+
+# hide deprecation warnings temporarily for stable release
+set(Kokkos_ENABLE_DEPRECATION_WARNINGS OFF CACHE BOOL "" FORCE)

cmake/presets/kokkos-hip.cmake

@@ -0,0 +1,20 @@
+# preset that enables KOKKOS and selects HIP compilation with OpenMP
+# enabled as well. Also sets some performance related compiler flags.
+set(PKG_KOKKOS ON CACHE BOOL "" FORCE)
+set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "" FORCE)
+set(Kokkos_ENABLE_OPENMP ON CACHE BOOL "" FORCE)
+set(Kokkos_ENABLE_CUDA   OFF CACHE BOOL "" FORCE)
+set(Kokkos_ENABLE_HIP    ON CACHE BOOL "" FORCE)
+set(Kokkos_ARCH_VEGA90A on CACHE BOOL "" FORCE)
+set(Kokkos_ENABLE_HIP_MULTIPLE_KERNEL_INSTANTIATIONS ON CACHE BOOL "" FORCE)
+set(BUILD_OMP ON CACHE BOOL "" FORCE)
+
+set(CMAKE_CXX_COMPILER hipcc CACHE STRING "" FORCE)
+set(CMAKE_TUNE_FLAGS "-munsafe-fp-atomics" CACHE STRING "" FORCE)
+
+# hide deprecation warnings temporarily for stable release
+set(Kokkos_ENABLE_DEPRECATION_WARNINGS OFF CACHE BOOL "" FORCE)
+
+# these flags are needed to build with Cray MPICH on OLCF Crusher
+#-D CMAKE_CXX_FLAGS="-I/${MPICH_DIR}/include"
+#-D MPI_CXX_LIBRARIES="-L${MPICH_DIR}/lib -lmpi -L${CRAY_MPICH_ROOTDIR}/gtl/lib -lmpi_gtl_hsa"

cmake/presets/kokkos-openmp.cmake

@@ -4,3 +4,6 @@ set(Kokkos_ENABLE_SERIAL ON  CACHE BOOL "" FORCE)
 set(Kokkos_ENABLE_OPENMP ON  CACHE BOOL "" FORCE)
 set(Kokkos_ENABLE_CUDA   OFF CACHE BOOL "" FORCE)
 set(BUILD_OMP ON CACHE BOOL "" FORCE)
+
+# hide deprecation warnings temporarily for stable release
+set(Kokkos_ENABLE_DEPRECATION_WARNINGS OFF CACHE BOOL "" FORCE)

cmake/presets/kokkos-serial.cmake

@@ -3,3 +3,6 @@ set(PKG_KOKKOS ON CACHE BOOL "" FORCE)
 set(Kokkos_ENABLE_SERIAL ON  CACHE BOOL "" FORCE)
 set(Kokkos_ENABLE_OPENMP OFF CACHE BOOL "" FORCE)
 set(Kokkos_ENABLE_CUDA   OFF CACHE BOOL "" FORCE)
+
+# hide deprecation warnings temporarily for stable release
+set(Kokkos_ENABLE_DEPRECATION_WARNINGS OFF CACHE BOOL "" FORCE)

cmake/presets/kokkos-sycl.cmake

@@ -0,0 +1,18 @@
+# preset that enables KOKKOS and selects SYCL compilation with OpenMP
+# enabled as well. Also sets some performance related compiler flags.
+set(PKG_KOKKOS ON CACHE BOOL "" FORCE)
+set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "" FORCE)
+set(Kokkos_ENABLE_OPENMP ON CACHE BOOL "" FORCE)
+set(Kokkos_ENABLE_CUDA   OFF CACHE BOOL "" FORCE)
+set(Kokkos_ENABLE_SYCL   ON CACHE BOOL "" FORCE)
+set(Kokkos_ARCH_MAXWELL50 on CACHE BOOL "" FORCE)
+set(BUILD_OMP ON CACHE BOOL "" FORCE)
+
+# hide deprecation warnings temporarily for stable release
+set(Kokkos_ENABLE_DEPRECATION_WARNINGS OFF CACHE BOOL "" FORCE)
+
+set(CMAKE_CXX_COMPILER clang++ CACHE STRING "" FORCE)
+set(MPI_CXX_COMPILER "mpicxx" CACHE STRING "" FORCE)
+set(CMAKE_CXX_STANDARD 17 CACHE STRING "" FORCE)
+set(CMAKE_SHARED_LINKER_FLAGS "-Xsycl-target-frontend -O3" CACHE STRING "" FORCE)
+set(CMAKE_TUNE_FLAGS "-fgpu-inline-threshold=100000 -Xsycl-target-frontend -O3  -Xsycl-target-frontend -ffp-contract=on -Wno-unknown-cuda-version" CACHE STRING "" FORCE)

cmake/presets/mingw-cross.cmake

@@ -4,6 +4,7 @@ set(WIN_PACKAGES
   AWPMD
   BOCS
   BODY
+  BPM
   BROWNIAN
   CG-DNA
   CG-SDK
@@ -20,6 +21,7 @@ set(WIN_PACKAGES
   DPD-REACT
   DPD-SMOOTH
   DRUDE
+  ELECTRODE
   EFF
   EXTRA-COMPUTE
   EXTRA-DUMP

cmake/presets/most.cmake

@@ -6,6 +6,7 @@ set(ALL_PACKAGES
   ASPHERE
   BOCS
   BODY
+  BPM
   BROWNIAN
   CG-DNA
   CG-SDK
@@ -22,6 +23,7 @@ set(ALL_PACKAGES
   DPD-REACT
   DPD-SMOOTH
   DRUDE
+  ELECTRODE
   EFF
   EXTRA-COMPUTE
   EXTRA-DUMP

cmake/presets/nolib.cmake

@@ -6,6 +6,7 @@ set(PACKAGES_WITH_LIB
   ATC
   AWPMD
   COMPRESS
+  ELECTRODE
   GPU
   H5MD
   KIM
@@ -15,7 +16,6 @@ set(PACKAGES_WITH_LIB
   MACHDYN
   MDI
   MESONT
-  MESSAGE
   ML-HDNNP
   ML-PACE
   ML-QUIP

cmake/presets/pgi.cmake

@@ -1,4 +1,4 @@
-# preset that will enable clang/clang++ with support for MPI and OpenMP (on Linux boxes)
+# preset that will enable PGI (Nvidia) compilers with support for MPI and OpenMP (on Linux boxes)
 
 set(CMAKE_CXX_COMPILER "pgc++" CACHE STRING "" FORCE)
 set(CMAKE_C_COMPILER "pgcc" CACHE STRING "" FORCE)

cmake/presets/windows-intel-classic.cmake

@@ -0,0 +1,8 @@
+# preset that will enable Intel compilers with support for MPI and OpenMP (on Linux boxes)
+
+set(CMAKE_CXX_COMPILER "icl" CACHE STRING "" FORCE)
+set(CMAKE_C_COMPILER "icl" CACHE STRING "" FORCE)
+set(CMAKE_Fortran_COMPILER "ifort" CACHE STRING "" FORCE)
+
+unset(HAVE_OMP_H_INCLUDE CACHE)
+

cmake/presets/windows-intel-llvm.cmake

@@ -0,0 +1,8 @@
+# preset that will enable Intel compilers with support for MPI and OpenMP (on Linux boxes)
+
+set(CMAKE_CXX_COMPILER "icx" CACHE STRING "" FORCE)
+set(CMAKE_C_COMPILER "icx" CACHE STRING "" FORCE)
+set(CMAKE_Fortran_COMPILER "ifx" CACHE STRING "" FORCE)
+set(INTEL_LRT_MODE "C++11" CACHE STRING "" FORCE)
+unset(HAVE_OMP_H_INCLUDE CACHE)
+set(CMAKE_TUNE_FLAGS -Wno-unused-command-line-argument)

cmake/presets/windows.cmake

@@ -2,6 +2,7 @@ set(WIN_PACKAGES
   ASPHERE
   BOCS
   BODY
+  BPM
   BROWNIAN
   CG-DNA
   CG-SDK

doc/.gitignore

@@ -1,6 +1,6 @@
 /old
 /html
-/html-offline
+/fasthtml
 /epub
 /latex
 /mathjax

doc/Makefile

@@ -14,34 +14,22 @@ ANCHORCHECK    = $(VENV)/bin/rst_anchor_check
 SPHINXCONFIG   = $(BUILDDIR)/utils/sphinx-config
 MATHJAX        = $(SPHINXCONFIG)/_static/mathjax
 
-PYTHON         = $(shell which python3)
-DOXYGEN        = $(shell which doxygen)
-VIRTUALENV     = virtualenv
+PYTHON         = $(word 3,$(shell type python3))
+DOXYGEN        = $(word 3,$(shell type doxygen))
 HAS_PYTHON3    = NO
-HAS_VIRTUALENV = NO
 HAS_DOXYGEN    = NO
 HAS_PDFLATEX   = NO
 
-ifeq ($(shell which python3 >/dev/null 2>&1; echo $$?), 0)
+ifeq ($(shell type python3 >/dev/null 2>&1; echo $$?), 0)
 HAS_PYTHON3    = YES
 endif
 
-ifeq ($(shell which doxygen >/dev/null 2>&1; echo $$?), 0)
+ifeq ($(shell type doxygen >/dev/null 2>&1; echo $$?), 0)
 HAS_DOXYGEN    = YES
 endif
 
-ifeq ($(shell which virtualenv-3 >/dev/null 2>&1; echo $$?), 0)
-VIRTUALENV     = virtualenv-3
-HAS_VIRTUALENV = YES
-endif
-
-ifeq ($(shell which virtualenv >/dev/null 2>&1; echo $$?), 0)
-VIRTUALENV     = virtualenv
-HAS_VIRTUALENV = YES
-endif
-
-ifeq ($(shell which pdflatex >/dev/null 2>&1; echo $$?), 0)
-ifeq ($(shell which latexmk >/dev/null 2>&1; echo $$?), 0)
+ifeq ($(shell type pdflatex >/dev/null 2>&1; echo $$?), 0)
+ifeq ($(shell type latexmk >/dev/null 2>&1; echo $$?), 0)
 HAS_PDFLATEX = YES
 endif
 endif
@@ -58,7 +46,7 @@ SPHINXEXTRA = -E -j $(shell $(PYTHON) -c 'import multiprocessing;print(multiproc
 # we only want to use explicitly listed files.
 DOXYFILES      = $(shell sed -n -e 's/\#.*$$//' -e '/^ *INPUT \+=/,/^[A-Z_]\+ \+=/p' doxygen/Doxyfile.in | sed -e 's/@LAMMPS_SOURCE_DIR@/..\/src/g' -e 's/\\//g' -e 's/ \+/ /' -e 's/[A-Z_]\+ \+= *\(YES\|NO\|\)//') 
 
-.PHONY: help clean-all clean clean-spelling epub mobi rst html pdf spelling anchor_check style_check char_check xmlgen
+.PHONY: help clean-all clean clean-spelling epub mobi rst html pdf spelling anchor_check style_check char_check xmlgen fasthtml
 
 # ------------------------------------------
 
@@ -70,7 +58,8 @@ help:
 	@echo "  epub          create ePUB format manual for e-book readers"
 	@echo "  mobi          convert ePUB to MOBI format manual for e-book readers (e.g. Kindle)"
 	@echo "                      (requires ebook-convert tool from calibre)"
-	@echo "  clean         remove all intermediate RST files"
+	@echo "  fasthtml      approximate HTML page creation in fasthtml dir (for development)"
+	@echo "  clean         remove all intermediate files"
 	@echo "  clean-all     reset the entire build environment"
 	@echo "  anchor_check  scan for duplicate anchor labels"
 	@echo "  style_check   check for complete and consistent style lists"
@@ -83,7 +72,7 @@ clean-all: clean
 	rm -rf $(BUILDDIR)/docenv $(MATHJAX) $(BUILDDIR)/LAMMPS.mobi $(BUILDDIR)/LAMMPS.epub $(BUILDDIR)/Manual.pdf
 
 clean: clean-spelling
-	rm -rf $(BUILDDIR)/html $(BUILDDIR)/epub $(BUILDDIR)/latex $(BUILDDIR)/doctrees $(BUILDDIR)/doxygen/xml $(BUILDDIR)/doxygen-warn.log $(BUILDDIR)/doxygen/Doxyfile $(SPHINXCONFIG)/conf.py
+	rm -rf $(BUILDDIR)/html $(BUILDDIR)/epub $(BUILDDIR)/latex $(BUILDDIR)/doctrees $(BUILDDIR)/doxygen/xml $(BUILDDIR)/doxygen-warn.log $(BUILDDIR)/doxygen/Doxyfile $(SPHINXCONFIG)/conf.py $(BUILDDIR)/fasthtml
 
 clean-spelling:
 	rm -rf $(BUILDDIR)/spelling
@@ -118,6 +107,23 @@ html: xmlgen $(VENV) $(SPHINXCONFIG)/conf.py $(ANCHORCHECK) $(MATHJAX)
 	@rm -rf html/PDF/.[sg]*
 	@echo "Build finished. The HTML pages are in doc/html."
 
+fasthtml: xmlgen $(VENV) $(SPHINXCONFIG)/conf.py $(ANCHORCHECK) $(MATHJAX)
+	@if [ "$(HAS_BASH)" == "NO" ] ; then echo "bash was not found at $(OSHELL)! Please use: $(MAKE) SHELL=/path/to/bash" 1>&2; exit 1; fi
+	@$(MAKE) $(MFLAGS) -C graphviz all
+	@mkdir -p fasthtml
+	@(\
+		. $(VENV)/bin/activate ; env PYTHONWARNINGS= \
+		sphinx-build -j 4 -b html -c $(SPHINXCONFIG) -d $(BUILDDIR)/fasthtml/doctrees $(RSTDIR) fasthtml ;\
+		deactivate ;\
+	)
+	@rm -rf fasthtml/_sources
+	@rm -rf fasthtml/PDF
+	@rm -rf fasthtml/USER
+	@rm -rf fasthtml/JPG
+	@cp -r src/PDF fasthtml/PDF
+	@rm -rf fasthtml/PDF/.[sg]*
+	@echo "Fast HTML build finished. The HTML pages are in doc/fasthtml."
+
 spelling: xmlgen $(SPHINXCONFIG)/conf.py $(VENV) $(SPHINXCONFIG)/false_positives.txt
 	@if [ "$(HAS_BASH)" == "NO" ] ; then echo "bash was not found at $(OSHELL)! Please use: $(MAKE) SHELL=/path/to/bash" 1>&2; exit 1; fi
 	@(\
@@ -220,9 +226,8 @@ $(VENV):
 	@if [ "$(HAS_BASH)" == "NO" ] ; then echo "bash was not found at $(OSHELL)! Please use: $(MAKE) SHELL=/path/to/bash" 1>&2; exit 1; fi
 	@if [ "$(HAS_PYTHON3)" == "NO" ] ; then echo "python3 was not found! Please see README for further instructions" 1>&2; exit 1; fi
 	@if [ "$(HAS_DOXYGEN)" == "NO" ] ; then echo "doxygen was not found! Please see README for further instructions" 1>&2; exit 1; fi
-	@if [ "$(HAS_VIRTUALENV)" == "NO" ] ; then echo "virtualenv was not found! Please see README for further instructions" 1>&2; exit 1; fi
 	@( \
-		$(VIRTUALENV) -p $(PYTHON) $(VENV); \
+		$(PYTHON) -m venv $(VENV); \
 		. $(VENV)/bin/activate; \
 		pip $(PIP_OPTIONS) install --upgrade pip; \
 		pip $(PIP_OPTIONS) install -r $(BUILDDIR)/utils/requirements.txt; \

doc/doxygen/Doxyfile.in

@@ -2,7 +2,7 @@
 
 DOXYFILE_ENCODING      = UTF-8
 PROJECT_NAME           = "LAMMPS Programmer's Guide"
-PROJECT_NUMBER         = "24 August 2020"
+PROJECT_NUMBER         = "4 May 2022"
 PROJECT_BRIEF          = "Documentation of the LAMMPS library interface and Python wrapper"
 PROJECT_LOGO           = lammps-logo.png
 CREATE_SUBDIRS         = NO
@@ -437,6 +437,8 @@ INPUT                  = @LAMMPS_SOURCE_DIR@/utils.cpp                 \
                          @LAMMPS_SOURCE_DIR@/math_eigen.h              \
                          @LAMMPS_SOURCE_DIR@/platform.h                \
                          @LAMMPS_SOURCE_DIR@/platform.cpp              \
+                         @LAMMPS_SOURCE_DIR@/math_special.h            \
+                         @LAMMPS_SOURCE_DIR@/math_special.cpp          \
 
 # The EXCLUDE_SYMLINKS tag can be used to select whether or not files or
 # directories that are symbolic links (a Unix file system feature) are excluded

doc/lammps.1

@@ -1,4 +1,4 @@
-.TH LAMMPS "1" "24 March 2022" "2022-3-24"
+.TH LAMMPS "1" "4 May 2022" "2022-5-4"
 .SH NAME
 .B LAMMPS
 \- Molecular Dynamics Simulator.  Version 24 March 2022

doc/msi2lmp.1

@@ -98,7 +98,7 @@ msi2lmp decane -c 0 -f oplsaa
 
 
 .SH COPYRIGHT
-© 2003--2021 Sandia Corporation
+© 2003--2022 Sandia Corporation
 
 This package is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License version 2 as

doc/src/Build_development.rst

@@ -486,14 +486,14 @@ The following options are available.
    make fix-whitespace      # correct whitespace issues in files
    make check-homepage      # search for files with old LAMMPS homepage URLs
    make fix-homepage        # correct LAMMPS homepage URLs in files
+   make check-errordocs     # search for deprecated error docs in header files
+   make fix-errordocs       # remove error docs in header files
    make check-permissions   # search for files with permissions issues
    make fix-permissions     # correct permissions issues in files
+   make check               # run all check targets from above
 
-These should help to replace all TAB characters with blanks and remove
-any trailing whitespace.  Also all LAMMPS homepage URL references can be
-updated to the location change from Sandia to the lammps.org domain.
-And the permission check can remove executable permissions from non-executable
-files (like source code).
+These should help to make source and documentation files conforming
+to some the coding style preferences of the LAMMPS developers.
 
 Clang-format support
 --------------------

doc/src/Build_extras.rst

@@ -45,7 +45,6 @@ This is the list of packages that may require additional steps.
    * :ref:`MACHDYN <machdyn>`
    * :ref:`MDI <mdi>`
    * :ref:`MESONT <mesont>`
-   * :ref:`MESSAGE <message>`
    * :ref:`ML-HDNNP <ml-hdnnp>`
    * :ref:`ML-IAP <mliap>`
    * :ref:`ML-PACE <ml-pace>`
@@ -142,8 +141,6 @@ CMake build
 
 :code:`GPU_ARCH` settings for different GPU hardware is as follows:
 
-* sm_12 or sm_13 for GT200 (supported by CUDA 3.2 until CUDA 6.5)
-* sm_20 or sm_21 for Fermi (supported by CUDA 3.2 until CUDA 7.5)
 * sm_30 for Kepler (supported since CUDA 5 and until CUDA 10.x)
 * sm_35 or sm_37 for Kepler (supported since CUDA 5 and until CUDA 11.x)
 * sm_50 or sm_52 for Maxwell (supported since CUDA 6)
@@ -161,6 +158,12 @@ Thus the GPU_ARCH setting is merely an optimization, to have code for
 the preferred GPU architecture directly included rather than having to wait
 for the JIT compiler of the CUDA driver to translate it.
 
+When compiling for CUDA or HIP with CUDA, version 8.0 or later of the CUDA toolkit
+is required and a GPU architecture of Kepler or later, which must *also* be
+supported by the CUDA toolkit in use **and** the CUDA driver in use.
+When compiling for OpenCL, OpenCL version 1.2 or later is required and the
+GPU must be supported by the GPU driver and OpenCL runtime bundled with the driver.
+
 When building with CMake, you **must NOT** build the GPU library in ``lib/gpu``
 using the traditional build procedure. CMake will detect files generated by that
 process and will terminate with an error and a suggestion for how to remove them.
@@ -638,13 +641,27 @@ This list was last updated for version 3.5.0 of the Kokkos library.
 
          -D CMAKE_CXX_COMPILER=${HOME}/lammps/lib/kokkos/bin/nvcc_wrapper
 
-      To simplify compilation, three preset files are included in the
+      For AMD or NVIDIA GPUs using HIP, set these variables:
+
+      .. code-block:: bash
+
+         -D Kokkos_ARCH_HOSTARCH=yes   # HOSTARCH = HOST from list above
+         -D Kokkos_ARCH_GPUARCH=yes    # GPUARCH = GPU from list above
+         -D Kokkos_ENABLE_HIP=yes
+         -D Kokkos_ENABLE_OPENMP=yes
+
+      This will enable FFTs on the GPU, either by the internal KISSFFT library
+      or with the hipFFT wrapper library, which will call out to the
+      platform-appropriate vendor library: rocFFT on AMD GPUs or cuFFT on
+      NVIDIA GPUs.
+
+      To simplify compilation, five preset files are included in the
       ``cmake/presets`` folder, ``kokkos-serial.cmake``,
-      ``kokkos-openmp.cmake``, and ``kokkos-cuda.cmake``. They will
-      enable the KOKKOS package and enable some hardware choice.  So to
-      compile with OpenMP host parallelization, CUDA device
-      parallelization (for GPUs with CC 5.0 and up) with some common
-      packages enabled, you can do the following:
+      ``kokkos-openmp.cmake``, ``kokkos-cuda.cmake``,
+      ``kokkos-hip.cmake``, and ``kokkos-sycl.cmake``.  They will enable
+      the KOKKOS package and enable some hardware choice.  So to compile
+      with CUDA device parallelization (for GPUs with CC 5.0 and up)
+      with some common packages enabled, you can do the following:
 
       .. code-block:: bash
 
@@ -703,6 +720,15 @@ This list was last updated for version 3.5.0 of the Kokkos library.
          KOKKOS_ABSOLUTE_PATH = $(shell cd $(KOKKOS_PATH); pwd)
          CC = mpicxx -cxx=$(KOKKOS_ABSOLUTE_PATH)/config/nvcc_wrapper
 
+      For AMD or NVIDIA GPUs using HIP:
+
+      .. code-block:: make
+
+         KOKKOS_DEVICES = HIP
+         KOKKOS_ARCH = HOSTARCH,GPUARCH  # HOSTARCH = HOST from list above that is hosting the GPU
+                                         # GPUARCH = GPU from list above
+         FFT_INC = -DFFT_HIPFFT           # enable use of hipFFT (optional)
+         FFT_LIB = -lhipfft               # link to hipFFT library
 
 Advanced KOKKOS compilation settings
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -796,47 +822,6 @@ library.
 
 ----------
 
-.. _message:
-
-MESSAGE package
------------------------------
-
-This package can optionally include support for messaging via sockets,
-using the open-source `ZeroMQ library <http://zeromq.org>`_, which must
-be installed on your system.
-
-.. tabs::
-
-   .. tab:: CMake build
-
-      .. code-block:: bash
-
-         -D MESSAGE_ZMQ=value    # build with ZeroMQ support, value = no (default) or yes
-         -D ZMQ_LIBRARY=path     # ZMQ library file (only needed if a custom location)
-         -D ZMQ_INCLUDE_DIR=path # ZMQ include directory (only needed if a custom location)
-
-   .. tab:: Traditional make
-
-      Before building LAMMPS, you must build the CSlib library in
-      ``lib/message``\ .  You can build the CSlib library manually if
-      you prefer; follow the instructions in ``lib/message/README``\ .
-      You can also do it in one step from the ``lammps/src`` dir, using
-      a command like these, which simply invoke the
-      ``lib/message/Install.py`` script with the specified args:
-
-      .. code-block:: bash
-
-         $ make lib-message               # print help message
-         $ make lib-message args="-m -z"  # build with MPI and socket (ZMQ) support
-         $ make lib-message args="-s"     # build as serial lib with no ZMQ support
-
-      The build should produce two files: ``lib/message/cslib/src/libmessage.a``
-      and ``lib/message/Makefile.lammps``.  The latter is copied from an
-      existing ``Makefile.lammps.*`` and has settings to link with the ZeroMQ
-      library if requested in the build.
-
-----------
-
 .. _mliap:
 
 ML-IAP package
@@ -1289,6 +1274,41 @@ be built for the most part with all major versions of the C++ language.
 
 ----------
 
+.. _electrode:
+
+ELECTRODE package
+-----------------
+
+This package depends on the KSPACE package.
+
+.. tabs::
+
+   .. tab:: CMake build
+
+      No additional settings are needed besides ``-D PKG_KSPACE=yes`` and ``-D
+      PKG_ELECTRODE=yes``.
+
+   .. tab:: Traditional make
+
+      The package is activated with ``make yes-KSPACE`` and ``make
+      yes-ELECTRODE``
+
+
+      Note that the ``Makefile.lammps`` file has settings for the BLAS and
+      LAPACK linear algebra libraries.  As explained in ``lib/awpmd/README``
+      these can either exist on your system, or you can use the files provided
+      in ``lib/linalg``.  In the latter case you also need to build the library
+      in ``lib/linalg`` with a command like these:
+
+      .. code-block:: bash
+
+         $ make lib-linalg                     # print help message
+         $ make lib-linalg args="-m serial"    # build with GNU Fortran compiler (settings as with "make serial")
+         $ make lib-linalg args="-m mpi"       # build with default MPI Fortran compiler (settings as with "make mpi")
+         $ make lib-linalg args="-m gfortran"  # build with GNU Fortran compiler
+
+----------
+
 .. _ml-pace:
 
 ML-PACE package

doc/src/Build_manual.rst

@@ -78,11 +78,12 @@ folder.  The following ``make`` commands are available:
 
    make html          # generate HTML in html dir using Sphinx
    make pdf           # generate PDF  as Manual.pdf using Sphinx and PDFLaTeX
-   make fetch         # fetch HTML pages and PDF files from LAMMPS website
-                      #  and unpack into the html_www folder and Manual_www.pdf
    make epub          # generate LAMMPS.epub in ePUB format using Sphinx
    make mobi          # generate LAMMPS.mobi in MOBI format using ebook-convert
 
+   make fasthtml      # generate approximate HTML in fasthtml dir using Sphinx
+                      # some Sphinx extensions do not work correctly with this
+
    make clean         # remove intermediate RST files created by HTML build
    make clean-all     # remove entire build folder and any cached data
 
@@ -193,8 +194,13 @@ folder need to be updated or new files added. These files are written in
 `reStructuredText <rst_>`_ markup for translation with the Sphinx tool.
 
 Before contributing any documentation, please check that both the HTML
-and the PDF format documentation can translate without errors. Please also
-check the output to the console for any warnings or problems.  There will
+and the PDF format documentation can translate without errors.  During
+testing the html translation, you may use the ``make fasthtml`` command
+which does an approximate translation (i.e. not all Sphinx features and
+extensions will work), but runs very fast because it will only translate
+files that have been changed since the last ``make fasthtml`` command.
+
+Please also check the output to the console for any warnings or problems.  There will
 be multiple tests run automatically:
 
 - A test for correctness of all anchor labels and their references

doc/src/Build_package.rst

@@ -45,7 +45,6 @@ packages:
    * :ref:`KOKKOS <kokkos>`
    * :ref:`LATTE <latte>`
    * :ref:`MACHDYN <machdyn>`
-   * :ref:`MESSAGE <message>`
    * :ref:`ML-HDNNP <ml-hdnnp>`
    * :ref:`ML-PACE <ml-pace>`
    * :ref:`ML-QUIP <ml-quip>`
@@ -151,7 +150,7 @@ other files dependent on that package are also excluded.
 .. _cmake_presets:
 
 CMake presets for installing many packages
-""""""""""""""""""""""""""""""""""""""""""
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 Instead of specifying all the CMake options via the command-line,
 CMake allows initializing its settings cache using script files.
@@ -179,6 +178,11 @@ one of them as a starting point and customize it to your needs.
     cmake -C ../cmake/presets/all_off.cmake  [OPTIONS] ../cmake  # disable all packages
     mingw64-cmake -C ../cmake/presets/mingw-cross.cmake [OPTIONS] ../cmake  #  compile with MinGW cross compilers
 
+Presets that have names starting with "windows" are specifically for
+compiling LAMMPS :doc:`natively on Windows <Build_windows>` and
+presets that have names starting with "kokkos" are specifically for
+selecting configurations for compiling LAMMPS with :ref:`KOKKOS <kokkos>`.
+
 .. note::
 
    Running cmake this way manipulates the CMake settings cache in your
@@ -221,7 +225,8 @@ These commands install/un-install sets of packages:
 .. code-block:: bash
 
     make yes-all                        # install all packages
-    make no-all                         # uninstall all packages
+    make no-all                         # check for changes and uninstall all packages
+    make no-installed                   # only check and uninstall installed packages
     make yes-basic                      # install a few commonly used packages'
     make no-basic                       # remove a few commonly used packages'
     make yes-most                       # install most packages w/o libs'

doc/src/Build_windows.rst

@@ -5,6 +5,7 @@ Notes for building LAMMPS on Windows
 * :ref:`Running Linux on Windows <linux>`
 * :ref:`Using GNU GCC ported to Windows <gnu>`
 * :ref:`Using Visual Studio <msvc>`
+* :ref:`Using Intel oneAPI compilers and libraries <oneapi>`
 * :ref:`Using a cross-compiler <cross>`
 
 ----------
@@ -25,8 +26,10 @@ assistance in resolving portability issues.  This is particularly true
 for compiling LAMMPS on Windows, since this platform has significant
 differences in some low-level functionality.  As of LAMMPS version 14
 December 2021, large parts of LAMMPS can be compiled natively with the
-Microsoft Visual C++ Compilers.  This is largely facilitated by using
-the :doc:`Developer_platform` in the ``platform`` namespace.
+Microsoft Visual C++ Compilers.  As of LAMMPS version 31 May 2022, also
+the Intel oneAPI compilers can compile large parts of LAMMPS natively on
+Windows.  This is mostly facilitated by using the
+:doc:`Developer_platform` in the ``platform`` namespace and CMake.
 
 Before trying to build LAMMPS on Windows yourself, please consider the
 `pre-compiled Windows installer packages <https://packages.lammps.org/windows.html>`_
@@ -73,8 +76,9 @@ configuration should set this up automatically, but is untested.
 In case of problems, you are recommended to contact somebody with
 experience in using Cygwin.  If you do come across portability problems
 requiring changes to the LAMMPS source code, or figure out corrections
-yourself, please report them on the lammps-users mailing list, or file
-them as an issue or pull request on the LAMMPS GitHub project.
+yourself, please report them on the
+`LAMMPS forum at MatSci <https://matsci.org/c/lammps/lammps-development/>`_,
+or file them as an issue or pull request on the LAMMPS GitHub project.
 
 .. _msvc:
 
@@ -98,6 +102,10 @@ It is possible to use both the integrated CMake support of the Visual
 Studio IDE or use an external CMake installation (e.g. downloaded from
 cmake.org) to create build files and compile LAMMPS from the command line.
 
+Compilation via command line and unit tests are checked automatically
+for the LAMMPS development branch through
+`GitHub Actions <https://github.com/lammps/lammps/actions/workflows/compile-msvc.yml>`_.
+
 .. note::
 
    Versions of Visual Studio before version 17.1 may scan the entire
@@ -110,6 +118,10 @@ Please note, that for either approach CMake will create a so-called
 the command lines for building and testing LAMMPS must be adjusted
 accordingly.
 
+The LAMMPS cmake folder contains a ``CMakeSettings.json`` file with
+build configurations for MSVC compilers and the MS provided Clang
+compiler package in Debug and Release mode.
+
 To support running in parallel you can compile with OpenMP enabled using
 the OPENMP package or install Microsoft MPI (including the SDK) and compile
 LAMMPS with MPI enabled.
@@ -117,8 +129,55 @@ LAMMPS with MPI enabled.
 .. note::
 
    This is work in progress and you should contact the LAMMPS developers
-   via GitHub, the forum, or the mailing list, if you have questions or
-   LAMMPS specific problems.
+   via GitHub or the `LAMMPS forum at MatSci <https://matsci.org/c/lammps/lammps-development/>`_,
+   if you have questions or LAMMPS specific problems.
+
+.. _oneapi:
+
+Using Intel oneAPI Compilers and Libraries
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. versionadded:: 31May2022
+
+After installing the `Intel oneAPI
+<https://www.intel.com/content/www/us/en/developer/tools/oneapi/toolkits.html>`_
+base toolkit and the HPC toolkit, it is also possible to compile large
+parts of LAMMPS natively on Windows using Intel compilers.  The HPC
+toolkit provides two sets of C/C++ and Fortran compilers: the so-called
+"classic" compilers (``icl.exe`` and ``ifort.exe``) and newer, LLVM
+based compilers (``icx.exe`` and ``ifx.exe``).  In addition to the
+compilers and their dependent modules, also the thread building blocks
+(TBB) and the math kernel library (MKL) need to be installed.  Two
+presets (``cmake/presets/windows-intel-llvm.cmake`` and
+``cmake/presets/windows-intel-classic.cmake``) are provided for
+selecting the LLVM based or classic compilers, respectively. The preset
+``cmake/presets/windows.cmake`` enables compatible packages that are not
+dependent on additional features or libraries.  You **must** use the
+CMake based build procedure and use Ninja as build tool.  For compiling
+from the command prompt, thus both `CMake <https://cmake.org>`_ and
+`Ninja-build <https://ninja-build.org>`_ binaries must be installed.  It
+is also possible to use Visual Studio, if it is started (``devenv.exe``)
+from a command prompt that has the Intel oneAPI compilers enabled.  The
+Visual Studio settings file in the ``cmake`` folder contains
+configurations for both compiler variants in debug and release settings.
+Those will use the CMake and Ninja binaries bundled with Visual Studio,
+thus a separate installation is not required.
+
+.. admonition:: Known Limitations
+   :class: note
+
+   In addition to portability issues with several packages and external
+   libraries, the classic Intel compilers are currently not able to
+   compile the googletest libraries and thus enabling the ``-DENABLE_TESTING``
+   option will result in compilation failure.  The LLVM based compilers
+   are compatible.
+
+.. note::
+
+   This is work in progress and you should contact the LAMMPS developers
+   via GitHub or the `LAMMPS forum at MatSci <https://matsci.org/c/lammps/lammps-development/>`_,
+   if you have questions or LAMMPS specific problems.
+
 
 .. _cross:
 
@@ -144,14 +203,3 @@ LAMMPS developers.  We instead rely on the feedback of the users
 of these pre-compiled LAMMPS packages for Windows.  We will try to resolve
 issues to the best of our abilities if we become aware of them. However
 this is subject to time constraints and focus on HPC platforms.
-
-.. _native:
-
-Native Visual C++ support
-^^^^^^^^^^^^^^^^^^^^^^^^^
-
-Support for the Visual C++ compilers is currently not available. The
-CMake build system is capable of creating suitable a Visual Studio
-style build environment, but the LAMMPS source code itself is not
-ported to fully support Visual C++. Volunteers to take on this task
-are welcome.

doc/src/Commands_all.rst

@@ -14,7 +14,7 @@
 General commands
 ================
 
-An alphabetic list of all general LAMMPS commands.
+An alphabetic list of general LAMMPS commands.
 
 .. table_from_list::
    :columns: 5
@@ -47,36 +47,26 @@ An alphabetic list of all general LAMMPS commands.
    * :doc:`displace_atoms <displace_atoms>`
    * :doc:`dump <dump>`
    * :doc:`dump_modify <dump_modify>`
-   * :doc:`dynamical_matrix (k) <dynamical_matrix>`
    * :doc:`echo <echo>`
    * :doc:`fix <fix>`
    * :doc:`fix_modify <fix_modify>`
    * :doc:`group <group>`
-   * :doc:`group2ndx <group2ndx>`
-   * :doc:`hyper <hyper>`
    * :doc:`if <if>`
    * :doc:`improper_coeff <improper_coeff>`
    * :doc:`improper_style <improper_style>`
    * :doc:`include <include>`
    * :doc:`info <info>`
    * :doc:`jump <jump>`
-   * :doc:`kim <kim_commands>`
    * :doc:`kspace_modify <kspace_modify>`
    * :doc:`kspace_style <kspace_style>`
    * :doc:`label <label>`
    * :doc:`lattice <lattice>`
    * :doc:`log <log>`
    * :doc:`mass <mass>`
-   * :doc:`mdi/engine <mdi_engine>`
-   * :doc:`message <message>`
    * :doc:`minimize <minimize>`
    * :doc:`min_modify <min_modify>`
    * :doc:`min_style <min_style>`
-   * :doc:`min_style spin <min_spin>`
    * :doc:`molecule <molecule>`
-   * :doc:`ndx2group <group2ndx>`
-   * :doc:`neb <neb>`
-   * :doc:`neb/spin <neb_spin>`
    * :doc:`neigh_modify <neigh_modify>`
    * :doc:`neighbor <neighbor>`
    * :doc:`newton <newton>`
@@ -87,11 +77,8 @@ An alphabetic list of all general LAMMPS commands.
    * :doc:`pair_style <pair_style>`
    * :doc:`pair_write <pair_write>`
    * :doc:`partition <partition>`
-   * :doc:`plugin <plugin>`
-   * :doc:`prd <prd>`
    * :doc:`print <print>`
    * :doc:`processors <processors>`
-   * :doc:`python <python>`
    * :doc:`quit <quit>`
    * :doc:`read_data <read_data>`
    * :doc:`read_dump <read_dump>`
@@ -105,19 +92,13 @@ An alphabetic list of all general LAMMPS commands.
    * :doc:`restart <restart>`
    * :doc:`run <run>`
    * :doc:`run_style <run_style>`
-   * :doc:`server <server>`
    * :doc:`set <set>`
    * :doc:`shell <shell>`
    * :doc:`special_bonds <special_bonds>`
    * :doc:`suffix <suffix>`
-   * :doc:`tad <tad>`
-   * :doc:`temper <temper>`
-   * :doc:`temper/grem <temper_grem>`
-   * :doc:`temper/npt <temper_npt>`
    * :doc:`thermo <thermo>`
    * :doc:`thermo_modify <thermo_modify>`
    * :doc:`thermo_style <thermo_style>`
-   * :doc:`third_order (k) <third_order>`
    * :doc:`timer <timer>`
    * :doc:`timestep <timestep>`
    * :doc:`uncompute <uncompute>`
@@ -130,3 +111,27 @@ An alphabetic list of all general LAMMPS commands.
    * :doc:`write_data <write_data>`
    * :doc:`write_dump <write_dump>`
    * :doc:`write_restart <write_restart>`
+
+Additional general LAMMPS commands provided by packages.  A few
+commands have accelerated versions.  This is indicated by an
+additional letter in parenthesis: k = KOKKOS.
+
+.. table_from_list::
+   :columns: 5
+
+   * :doc:`dynamical_matrix (k) <dynamical_matrix>`
+   * :doc:`group2ndx <group2ndx>`
+   * :doc:`hyper <hyper>`
+   * :doc:`kim <kim_commands>`
+   * :doc:`mdi <mdi>`
+   * :doc:`ndx2group <group2ndx>`
+   * :doc:`neb <neb>`
+   * :doc:`neb/spin <neb_spin>`
+   * :doc:`plugin <plugin>`
+   * :doc:`prd <prd>`
+   * :doc:`python <python>`
+   * :doc:`tad <tad>`
+   * :doc:`temper <temper>`
+   * :doc:`temper/grem <temper_grem>`
+   * :doc:`temper/npt <temper_npt>`
+   * :doc:`third_order (k) <third_order>`

doc/src/Commands_bond.rst

@@ -32,6 +32,8 @@ OPT.
    *
    *
    *
+   * :doc:`bpm/rotational <bond_bpm_rotational>`
+   * :doc:`bpm/spring <bond_bpm_spring>`
    * :doc:`class2 (ko) <bond_class2>`
    * :doc:`fene (iko) <bond_fene>`
    * :doc:`fene/expand (o) <bond_fene_expand>`
@@ -86,7 +88,7 @@ OPT.
    * :doc:`dipole (o) <angle_dipole>`
    * :doc:`fourier (o) <angle_fourier>`
    * :doc:`fourier/simple (o) <angle_fourier_simple>`
-   * :doc:`gaussian <angle_gaussian>` - multicentered Gaussian-based angle potential
+   * :doc:`gaussian <angle_gaussian>`
    * :doc:`harmonic (iko) <angle_harmonic>`
    * :doc:`mm3 <angle_mm3>`
    * :doc:`quartic (o) <angle_quartic>`

doc/src/Commands_compute.rst

@@ -33,6 +33,7 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`body/local <compute_body_local>`
    * :doc:`bond <compute_bond>`
    * :doc:`bond/local <compute_bond_local>`
+   * :doc:`born/matrix <compute_born_matrix>`
    * :doc:`centro/atom <compute_centro_atom>`
    * :doc:`centroid/stress/atom <compute_stress_atom>`
    * :doc:`chunk/atom <compute_chunk_atom>`
@@ -63,6 +64,7 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`event/displace <compute_event_displace>`
    * :doc:`fabric <compute_fabric>`
    * :doc:`fep <compute_fep>`
+   * :doc:`fep/ta <compute_fep_ta>`
    * :doc:`force/tally <compute_tally>`
    * :doc:`fragment/atom <compute_cluster_atom>`
    * :doc:`global/atom <compute_global_atom>`
@@ -90,6 +92,7 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`msd <compute_msd>`
    * :doc:`msd/chunk <compute_msd_chunk>`
    * :doc:`msd/nongauss <compute_msd_nongauss>`
+   * :doc:`nbond/atom <compute_nbond_atom>`
    * :doc:`omega/chunk <compute_omega_chunk>`
    * :doc:`orientorder/atom (k) <compute_orientorder_atom>`
    * :doc:`pair <compute_pair>`

doc/src/Commands_fix.rst

@@ -51,7 +51,6 @@ OPT.
    * :doc:`bond/swap <fix_bond_swap>`
    * :doc:`box/relax <fix_box_relax>`
    * :doc:`charge/regulation <fix_charge_regulation>`
-   * :doc:`client/md <fix_client_md>`
    * :doc:`cmap <fix_cmap>`
    * :doc:`colvars <fix_colvars>`
    * :doc:`controller <fix_controller>`
@@ -67,6 +66,9 @@ OPT.
    * :doc:`edpd/source <fix_dpd_source>`
    * :doc:`efield <fix_efield>`
    * :doc:`ehex <fix_ehex>`
+   * :doc:`electrode/conp (i) <fix_electrode_conp>`
+   * :doc:`electrode/conq (i) <fix_electrode_conp>`
+   * :doc:`electrode/thermo (i) <fix_electrode_conp>`
    * :doc:`electron/stopping <fix_electron_stopping>`
    * :doc:`electron/stopping/fit <fix_electron_stopping>`
    * :doc:`enforce2d (k) <fix_enforce2d>`
@@ -101,7 +103,7 @@ OPT.
    * :doc:`lb/viscous <fix_lb_viscous>`
    * :doc:`lineforce <fix_lineforce>`
    * :doc:`manifoldforce <fix_manifoldforce>`
-   * :doc:`mdi/engine <fix_mdi_engine>`
+   * :doc:`mdi/aimd <fix_mdi_aimd>`
    * :doc:`meso/move <fix_meso_move>`
    * :doc:`mol/swap <fix_mol_swap>`
    * :doc:`momentum (k) <fix_momentum>`
@@ -142,6 +144,7 @@ OPT.
    * :doc:`nve/manifold/rattle <fix_nve_manifold_rattle>`
    * :doc:`nve/noforce <fix_nve_noforce>`
    * :doc:`nve/sphere (ko) <fix_nve_sphere>`
+   * :doc:`nve/bpm/sphere <fix_nve_bpm_sphere>`
    * :doc:`nve/spin <fix_nve_spin>`
    * :doc:`nve/tri <fix_nve_tri>`
    * :doc:`nvk <fix_nvk>`
@@ -159,7 +162,6 @@ OPT.
    * :doc:`orient/fcc <fix_orient>`
    * :doc:`orient/eco <fix_orient_eco>`
    * :doc:`pafi <fix_pafi>`
-   * :doc:`pair/tracker <fix_pair_tracker>`
    * :doc:`phonon <fix_phonon>`
    * :doc:`pimd <fix_pimd>`
    * :doc:`planeforce <fix_planeforce>`

doc/src/Commands_kspace.rst

@@ -27,6 +27,7 @@ OPT.
    * :doc:`ewald/disp/dipole <kspace_style>`
    * :doc:`ewald/dipole <kspace_style>`
    * :doc:`ewald/dipole/spin <kspace_style>`
+   * :doc:`ewald/electrode <kspace_style>`
    * :doc:`msm (o) <kspace_style>`
    * :doc:`msm/cg (o) <kspace_style>`
    * :doc:`msm/dielectric <kspace_style>`
@@ -41,4 +42,5 @@ OPT.
    * :doc:`pppm/stagger <kspace_style>`
    * :doc:`pppm/tip4p (o) <kspace_style>`
    * :doc:`pppm/dielectric <kspace_style>`
+   * :doc:`pppm/electrode (i) <kspace_style>`
    * :doc:`scafacos <kspace_style>`

doc/src/Commands_pair.rst

@@ -34,7 +34,7 @@ OPT.
    *
    *
    *
-   * :doc:`adp (o) <pair_adp>`
+   * :doc:`adp (ko) <pair_adp>`
    * :doc:`agni (o) <pair_agni>`
    * :doc:`airebo (io) <pair_airebo>`
    * :doc:`airebo/morse (io) <pair_airebo>`
@@ -53,6 +53,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:`bpm/spring <pair_bpm_spring>`
    * :doc:`brownian (o) <pair_brownian>`
    * :doc:`brownian/poly (o) <pair_brownian>`
    * :doc:`buck (giko) <pair_buck>`
@@ -88,12 +89,12 @@ OPT.
    * :doc:`coul/tt <pair_coul_tt>`
    * :doc:`coul/wolf (ko) <pair_coul>`
    * :doc:`coul/wolf/cs <pair_cs>`
-   * :doc:`dpd (gio) <pair_dpd>`
+   * :doc:`dpd (giko) <pair_dpd>`
    * :doc:`dpd/fdt <pair_dpd_fdt>`
-   * :doc:`dpd/ext <pair_dpd_ext>`
-   * :doc:`dpd/ext/tstat <pair_dpd_ext>`
+   * :doc:`dpd/ext (k) <pair_dpd_ext>`
+   * :doc:`dpd/ext/tstat (k) <pair_dpd_ext>`
    * :doc:`dpd/fdt/energy (k) <pair_dpd_fdt>`
-   * :doc:`dpd/tstat (go) <pair_dpd>`
+   * :doc:`dpd/tstat (gko) <pair_dpd>`
    * :doc:`dsmc <pair_dsmc>`
    * :doc:`e3b <pair_e3b>`
    * :doc:`drip <pair_drip>`
@@ -123,8 +124,8 @@ OPT.
    * :doc:`hbond/dreiding/lj (o) <pair_hbond_dreiding>`
    * :doc:`hbond/dreiding/morse (o) <pair_hbond_dreiding>`
    * :doc:`hdnnp <pair_hdnnp>`
-   * :doc:`ilp/graphene/hbn <pair_ilp_graphene_hbn>`
-   * :doc:`ilp/tmd <pair_ilp_tmd>`
+   * :doc:`ilp/graphene/hbn (t) <pair_ilp_graphene_hbn>`
+   * :doc:`ilp/tmd (t) <pair_ilp_tmd>`
    * :doc:`kolmogorov/crespi/full <pair_kolmogorov_crespi_full>`
    * :doc:`kolmogorov/crespi/z <pair_kolmogorov_crespi_z>`
    * :doc:`lcbop <pair_lcbop>`
@@ -151,7 +152,7 @@ OPT.
    * :doc:`lj/cut/coul/cut/dielectric (o) <pair_dielectric>`
    * :doc:`lj/cut/coul/cut/soft (o) <pair_fep_soft>`
    * :doc:`lj/cut/coul/debye (gko) <pair_lj_cut_coul>`
-   * :doc:`lj/cut/coul/debye/dielectric <pair_dielectric>`
+   * :doc:`lj/cut/coul/debye/dielectric (o) <pair_dielectric>`
    * :doc:`lj/cut/coul/dsf (gko) <pair_lj_cut_coul>`
    * :doc:`lj/cut/coul/long (gikot) <pair_lj_cut_coul>`
    * :doc:`lj/cut/coul/long/cs <pair_cs>`
@@ -230,7 +231,7 @@ OPT.
    * :doc:`oxrna2/stk <pair_oxrna2>`
    * :doc:`oxrna2/xstk <pair_oxrna2>`
    * :doc:`oxrna2/coaxstk <pair_oxrna2>`
-   * :doc:`pace <pair_pace>`
+   * :doc:`pace (k) <pair_pace>`
    * :doc:`peri/eps <pair_peri>`
    * :doc:`peri/lps (o) <pair_peri>`
    * :doc:`peri/pmb (o) <pair_peri>`
@@ -242,8 +243,10 @@ OPT.
    * :doc:`reaxff (ko) <pair_reaxff>`
    * :doc:`rebo (io) <pair_airebo>`
    * :doc:`resquared (go) <pair_resquared>`
-   * :doc:`saip/metal <pair_saip_metal>`
+   * :doc:`saip/metal (t) <pair_saip_metal>`
    * :doc:`sdpd/taitwater/isothermal <pair_sdpd_taitwater_isothermal>`
+   * :doc:`smatb <pair_smatb>`
+   * :doc:`smatb/single <pair_smatb>`
    * :doc:`smd/hertz <pair_smd_hertz>`
    * :doc:`smd/tlsph <pair_smd_tlsph>`
    * :doc:`smd/tri_surface <pair_smd_triangulated_surface>`

doc/src/Commands_parse.rst

@@ -77,18 +77,19 @@ LAMMPS:
    so that you do not have to define (or discard) a temporary variable,
    "X" in this case.
 
-   Additionally, the "immediate" variable expression may be followed by
-   a colon, followed by a C-style format string, e.g. ":%f" or ":%.10g".
-   The format string must be appropriate for a double-precision
-   floating-point value.  The format string is used to output the result
-   of the variable expression evaluation.  If a format string is not
-   specified a high-precision "%.20g" is used as the default.
+   Additionally, the entire "immediate" variable expression may be
+   followed by a colon, followed by a C-style format string,
+   e.g. ":%f" or ":%.10g".  The format string must be appropriate for
+   a double-precision floating-point value.  The format string is used
+   to output the result of the variable expression evaluation.  If a
+   format string is not specified, a high-precision "%.20g" is used as
+   the default format.
 
    This can be useful for formatting print output to a desired precision:
 
    .. code-block:: LAMMPS
 
-      print "Final energy per atom: $(pe/atoms:%10.3f) eV/atom"
+      print "Final energy per atom: $(v_ke_per_atom+v_pe_per_atom:%10.3f) eV/atom"
 
    Note that neither the curly-bracket or immediate form of variables
    can contain nested $ characters for other variables to substitute

doc/src/Developer_plugins.rst

@@ -68,7 +68,7 @@ Members of ``lammpsplugin_t``
    * - author
      - String with the name and email of the author
    * - creator.v1
-     - Pointer to factory function for pair, bond, angle, dihedral, improper, command or kspace styles
+     - Pointer to factory function for pair, bond, angle, dihedral, improper, kspace, or command styles
    * - creator.v2
      - Pointer to factory function for compute, fix, or region styles
    * - handle
@@ -262,3 +262,24 @@ A plugin may be registered under an existing style name.  In that case
 the plugin will override the existing code.  This can be used to modify
 the behavior of existing styles or to debug new versions of them without
 having to re-compile or re-install all of LAMMPS.
+
+Compiling plugins
+^^^^^^^^^^^^^^^^^
+
+Plugins need to be compiled with the same compilers and libraries
+(e.g. MPI) and compilation settings (MPI on/off, OpenMP, integer sizes)
+as the LAMMPS executable and library.  Otherwise the plugin will likely
+not load due to mismatches in the function signatures (LAMMPS is C++ so
+scope, type, and number of arguments are encoded into the symbol names
+and thus differences in them will lead to failed plugin load commands.
+Compilation of the plugin can be managed via both, CMake or traditional
+GNU makefiles.  Some examples that can be used as a template are in the
+``examples/plugins`` folder.  The CMake script code has some small
+adjustments to allow building the plugins for running unit tests with
+them. Another example that converts the KIM package into a plugin can be
+found in the ``examples/kim/plugin`` folder.  No changes to the sources
+of the KIM package themselves are needed; only the plugin interface and
+loader code needs to be added.  This example only supports building with
+CMake, but is probably a more typical example. To compile you need to
+run CMake with -DLAMMPS_SOURCE_DIR=<path/to/lammps/src/folder>.  Other
+configuration setting are identical to those for compiling LAMMPS.

doc/src/Developer_utils.rst

@@ -133,6 +133,9 @@ and parsing files or arguments.
 .. doxygenfunction:: trim_comment
    :project: progguide
 
+.. doxygenfunction:: star_subst
+   :project: progguide
+
 .. doxygenfunction:: has_utf8
    :project: progguide
 
@@ -208,6 +211,12 @@ Convenience functions
 .. doxygenfunction:: logmesg(LAMMPS *lmp, const std::string &mesg)
    :project: progguide
 
+.. doxygenfunction:: errorurl
+   :project: progguide
+
+.. doxygenfunction:: missing_cmd_args
+   :project: progguide
+
 .. doxygenfunction:: flush_buffers(LAMMPS *lmp)
    :project: progguide
 
@@ -237,6 +246,44 @@ Customized standard functions
 
 ---------------------------
 
+Special Math functions
+----------------------
+
+The ``MathSpecial`` namespace implements a selection of custom and optimized
+mathematical functions for a variety of applications.
+
+.. doxygenfunction:: factorial
+   :project: progguide
+
+.. doxygenfunction:: exp2_x86
+   :project: progguide
+
+.. doxygenfunction:: fm_exp
+   :project: progguide
+
+.. doxygenfunction:: my_erfcx
+   :project: progguide
+
+.. doxygenfunction:: expmsq
+   :project: progguide
+
+.. doxygenfunction:: square
+   :project: progguide
+
+.. doxygenfunction:: cube
+   :project: progguide
+
+.. doxygenfunction:: powsign
+   :project: progguide
+
+.. doxygenfunction:: powint
+   :project: progguide
+
+.. doxygenfunction:: powsinxx
+   :project: progguide
+
+---------------------------
+
 Tokenizer classes
 -----------------
 

doc/src/Errors.rst

@@ -11,6 +11,7 @@ them.
    :maxdepth: 1
 
    Errors_common
+   Errors_details
    Errors_bugs
    Errors_debug
    Errors_messages

doc/src/Errors_bugs.rst

@@ -17,9 +17,8 @@ the steps outlined below:
    if your issue has already been reported and if it is still open.
  * Check the `GitHub Pull Requests page <https://github.com/lammps/lammps/pulls>`_
    to see if there is already a fix for your bug pending.
- * Check the `mailing list archives <https://www.lammps.org/mail.html>`_ or
-   the `LAMMPS forum <https://www.lammps.org/forum.html>`_ to see if the
-   issue has been discussed before.
+ * Check the `LAMMPS forum at MatSci <https://matsci.org/lammps/>`_
+   to see if the issue has been discussed before.
 
 If none of these steps yields any useful information, please file a new
 bug report on the `GitHub Issue page <https://github.com/lammps/lammps/issues>`_.
@@ -38,12 +37,9 @@ generate this restart from a data file or a simple additional input.
 This input deck can be used with tools like a debugger or `valgrind
 <https://valgrind.org>`_ to further :doc:`debug the crash <Errors_debug>`.
 
-You may also send an email to the LAMMPS mailing list at
-"lammps-users at lists.sourceforge.net" describing the problem with the
-same kind of information.  The mailing list can provide a faster response,
-especially if the bug reported is actually expected behavior.  But because
-of the high volume of the mailing list, it can happen that your e-mail
-is overlooked and then forgotten.  Issues on GitHub have to be explicitly
-closed, so that will *guarantee* that at least one LAMMPS developer will
-have looked at it.
+You may also post a message in the `development category of the LAMMPS
+forum at MatSci <https://matsci.org/c/lammps/lammps-development/>`_
+describing the problem with the same kind of information.  The forum can
+provide a faster response, especially if the bug reported is actually
+expected behavior or other LAMMPS users have come across it before.
 

doc/src/Errors_details.rst

@@ -0,0 +1,27 @@
+Error and warning details
+=========================
+
+Many errors or warnings are self-explanatory and thus straightforward to
+resolve.  However, there are also cases, where there is no single cause
+and explanation, where LAMMPS can only detect symptoms of an error but
+not the exact cause, or where the explanation needs to be more detailed than
+what can be fit into a message printed by the program.  The following are
+discussions of such cases.
+
+.. _err0001:
+
+Unknown identifier in data file
+-------------------------------
+
+This error happens when LAMMPS encounters a line of text in an unexpected format
+while reading a data file. This is most commonly cause by inconsistent header and
+section data.  The header section informs LAMMPS how many entries or lines are expected in the
+various sections (like Atoms, Masses, Pair Coeffs, *etc.*\ ) of the data file.
+If there is a mismatch, LAMMPS will either keep reading beyond the end of a section
+or stop reading before the section has ended.
+
+Such a mismatch can happen unexpectedly when the first line of the data
+is *not* a comment as required by the format.  That would result in
+LAMMPS expecting, for instance, 0 atoms because the "atoms" header line
+is treated as a comment.
+

doc/src/Examples.rst

@@ -54,6 +54,8 @@ Lowercase directories
 +-------------+------------------------------------------------------------------+
 | body        | body particles, 2d system                                        |
 +-------------+------------------------------------------------------------------+
+| bpm         | BPM simulations of pouring elastic grains and plate impact       |
++-------------+------------------------------------------------------------------+
 | cmap        | CMAP 5-body contributions to CHARMM force field                  |
 +-------------+------------------------------------------------------------------+
 | colloid     | big colloid particles in a small particle solvent, 2d system     |
@@ -94,12 +96,12 @@ Lowercase directories
 +-------------+------------------------------------------------------------------+
 | latte       | examples for using fix latte for DFTB via the LATTE library      |
 +-------------+------------------------------------------------------------------+
+| mdi         | use of the MDI package and MolSSI MDI code coupling library      |
++-------------+------------------------------------------------------------------+
 | meam        | MEAM test for SiC and shear (same as shear examples)             |
 +-------------+------------------------------------------------------------------+
 | melt        | rapid melt of 3d LJ system                                       |
 +-------------+------------------------------------------------------------------+
-| message     | demos for LAMMPS client/server coupling with the MESSAGE package |
-+-------------+------------------------------------------------------------------+
 | micelle     | self-assembly of small lipid-like molecules into 2d bilayers     |
 +-------------+------------------------------------------------------------------+
 | min         | energy minimization of 2d LJ melt                                |

doc/src/Howto.rst

@@ -22,8 +22,9 @@ General howto
    Howto_replica
    Howto_library
    Howto_couple
-   Howto_client_server
    Howto_mdi
+   Howto_bpm
+   Howto_broken_bonds
 
 Settings howto
 ==============

doc/src/Howto_bpm.rst

@@ -0,0 +1,119 @@
+Bonded particle models
+======================
+
+The BPM package implements bonded particle models which can be used to
+simulate mesoscale solids.  Solids are constructed as a collection of
+particles which each represent a coarse-grained region of space much
+larger than the atomistic scale. Particles within a solid region are
+then connected by a network of bonds to provide solid elasticity.
+
+Unlike traditional bonds in molecular dynamics, the equilibrium bond
+length can vary between bonds. Bonds store the reference state.  This
+includes setting the equilibrium length equal to the initial distance
+between the two particles but can also include data on the bond
+orientation for rotational models. This produces a stress free initial
+state. Furthermore, bonds are allowed to break under large strains
+producing fracture. The examples/bpm directory has sample input scripts
+for simulations of the fragmentation of an impacted plate and the
+pouring of extended, elastic bodies.
+
+----------
+
+Bonds can be created using a :doc:`read data <read_data>` or
+:doc:`create bonds <create_bonds>` command. Alternatively, a
+:doc:`molecule <molecule>` template with bonds can be used with
+:doc:`fix deposit <fix_deposit>` or :doc:`fix pour <fix_pour>` to
+create solid grains.
+
+In this implementation, bonds store their reference state when they are
+first computed in the setup of the first simulation run. Data is then
+preserved across run commands and is written to :doc:`binary restart
+files <restart>` such that restarting the system will not reset the
+reference state of a bond. Bonds that are created midway into a run,
+such as those created by pouring grains using :doc:`fix pour
+<fix_pour>`, are initialized on that timestep.
+
+As bonds can be broken between neighbor list builds, the
+:doc:`special_bonds <special_bonds>` command works differently for BPM
+bond styles. There are two possible settings which determine how pair
+interactions work between bonded particles.  First, one can turn off
+all pair interactions between bonded particles.  Unlike :doc:`bond
+quartic <bond_quartic>`, this is not done by subtracting pair forces
+during the bond computation but rather by dynamically updating the
+special bond list. This is the default behavior of BPM bond styles and
+is done by updating the 1-2 special bond list as bonds break.  To do
+this, LAMMPS requires :doc:`newton <newton>` bond off such that all
+processors containing an atom know when a bond breaks. Additionally,
+one must do either (A) or (B).
+
+A) Use the following special bond settings
+
+   .. code-block:: LAMMPS
+
+      special_bonds lj 0 1 1 coul 1 1 1
+
+   These settings accomplish two goals. First, they turn off 1-3 and 1-4
+   special bond lists, which are not currently supported for BPMs. As
+   BPMs often have dense bond networks, generating 1-3 and 1-4 special
+   bond lists is expensive.  By setting the lj weight for 1-2 bonds to
+   zero, this turns off pairwise interactions.  Even though there are no
+   charges in BPM models, setting a nonzero coul weight for 1-2 bonds
+   ensures all bonded neighbors are still included in the neighbor list
+   in case bonds break between neighbor list builds.
+
+B) Alternatively, one can simply overlay pair interactions such that all
+   bonded particles also feel pair interactions. This can be
+   accomplished by using the *overlay/pair* keyword present in all bpm
+   bond styles and by using the following special bond settings
+
+   .. code-block:: LAMMPS
+
+      special_bonds lj/coul 1 1 1
+
+See the :doc:`Howto <Howto_broken_bonds>` page on broken bonds for
+more information.
+
+----------
+
+Currently there are two types of bonds included in the BPM
+package. The first bond style, :doc:`bond bpm/spring
+<bond_bpm_spring>`, only applies pairwise, central body forces. Point
+particles must have :doc:`bond atom style <atom_style>` and may be
+thought of as nodes in a spring network. Alternatively, the second
+bond style, :doc:`bond bpm/rotational <bond_bpm_rotational>`, resolves
+tangential forces and torques arising with the shearing, bending, and
+twisting of the bond due to rotation or displacement of particles.
+Particles are similar to those used in the :doc:`granular package
+<Howto_granular>`, :doc:`atom style sphere <atom_style>`. However,
+they must also track the current orientation of particles and store bonds
+and therefore use a :doc:`bpm/sphere atom style <atom_style>`.
+This also requires a unique integrator :doc:`fix nve/bpm/sphere
+<fix_nve_bpm_sphere>` which numerically integrates orientation similar
+to :doc:`fix nve/asphere <fix_nve_asphere>`.
+
+To monitor the fracture of bonds in the system, all BPM bond styles
+have the ability to record instances of bond breakage to output using
+the :doc:`dump local <dump>` command. Additionally, one can use
+:doc:`compute nbond/atom <compute_nbond_atom>` to tally the current
+number of bonds per atom.
+
+In addition to bond styles, a new pair style :doc:`pair bpm/spring
+<pair_bpm_spring>` was added to accompany the bpm/spring bond
+style. This pair style is simply a hookean repulsion with similar
+velocity damping as its sister bond style.
+
+----------
+
+While LAMMPS has many utilities to create and delete bonds, *only*
+the following are currently compatible with BPM bond styles:
+
+* :doc:`create_bonds <create_bonds>`
+* :doc:`delete_bonds <delete_bonds>`
+* :doc:`fix bond/create <fix_bond_create>`
+* :doc:`fix bond/break <fix_bond_break>`
+* :doc:`fix bond/swap <fix_bond_swap>`
+
+Note :doc:`create_bonds <create_bonds>` requires certain special_bonds settings.
+To subtract pair interactions, one will need to switch between different
+special_bonds settings in the input script. An example is found in
+examples/bpm/impact.

doc/src/Howto_broken_bonds.rst

@@ -0,0 +1,48 @@
+Broken Bonds
+============
+
+Typically, bond interactions persist for the duration of a simulation
+in LAMMPS. However, there are some exceptions that allow for bonds to
+break including the :doc:`quartic bond style <bond_quartic>` and the
+bond styles in the :doc:`BPM package <Howto_bpm>` which contains the
+:doc:`bpm/spring <bond_bpm_spring>` and
+:doc:`bpm/rotational <bond_bpm_rotational>` bond styles. In these cases,
+a bond can be broken if it is stretched beyond a user-defined threshold.
+LAMMPS accomplishes this by setting the bond type to zero such that the
+bond force is no longer computed.
+
+Users are normally able to weight the contribution of pair forces to atoms
+that are bonded using the :doc:`special_bonds command <special_bonds>`.
+When bonds break, this is not always the case. For the quartic bond style,
+pair forces are always turned off between bonded particles. LAMMPS does
+this via a computational sleight-of-hand. It subtracts the pairwise
+interaction as part of the bond computation. When the bond breaks, the
+subtraction stops.  For this to work, the pairwise interaction must always
+be computed by the :doc:`pair_style <pair_style>` command, whether the bond
+is broken or not.  This means that :doc:`special_bonds <special_bonds>` must
+be set to 1,1,1. After the bond breaks, the pairwise interaction between the
+two atoms is turned on, since they are no longer bonded.
+
+In the BPM package, one can either turn off all pair interactions between
+bonded particles or leave them on, overlaying pair forces on top of bond
+forces. To remove pair forces, the special bond list is dynamically
+updated. More details can be found on the :doc:`Howto BPM <Howto_bpm>`
+page.
+
+Bonds can also be broken by fixes which change bond topology, including
+:doc:`fix bond/break <fix_bond_break>` and
+:doc:`fix bond/react <fix_bond_react>`. These fixes will automatically
+trigger a rebuild of the neighbor list and update special bond data structures
+when bonds are broken.
+
+Note that when bonds are dumped to a file via the :doc:`dump local <dump>` command, bonds with type 0 are not included.  The
+:doc:`delete_bonds <delete_bonds>` command can also be used to query the
+status of broken bonds or permanently delete them, e.g.:
+
+.. code-block:: LAMMPS
+
+   delete_bonds all stats
+   delete_bonds all bond 0 remove
+
+The compute :doc:`nbond/atom <compute_nbond_atom>` can also be used
+to tally the current number of bonds per atom, excluding broken bonds.

doc/src/Howto_client_server.rst

@@ -1,163 +0,0 @@
-Using LAMMPS in client/server mode
-==================================
-
-Client/server coupling of two codes is where one code is the "client"
-and sends request messages to a "server" code.  The server responds to
-each request with a reply message.  This enables the two codes to work
-in tandem to perform a simulation.  LAMMPS can act as either a client
-or server code.
-
-Some advantages of client/server coupling are that the two codes run
-as stand-alone executables; they are not linked together.  Thus
-neither code needs to have a library interface.  This often makes it
-easier to run the two codes on different numbers of processors.  If a
-message protocol (format and content) is defined for a particular kind
-of simulation, then in principle any code that implements the
-client-side protocol can be used in tandem with any code that
-implements the server-side protocol, without the two codes needing to
-know anything more specific about each other.
-
-A simple example of client/server coupling is where LAMMPS is the
-client code performing MD timestepping.  Each timestep it sends a
-message to a server quantum code containing current coords of all the
-atoms.  The quantum code computes energy and forces based on the
-coords.  It returns them as a message to LAMMPS, which completes the
-timestep.
-
-A more complex example is where LAMMPS is the client code and
-processes a series of data files, sending each configuration to a
-quantum code to compute energy and forces.  Or LAMMPS runs dynamics
-with an atomistic force field, but pauses every N steps to ask the
-quantum code to compute energy and forces.
-
-Alternate methods for code coupling with LAMMPS are described on
-the :doc:`Howto couple <Howto_couple>` doc page.
-
-The protocol for using LAMMPS as a client is to use these 3 commands
-in this order (other commands may come in between):
-
-* :doc:`message client <message>`  # initiate client/server interaction
-* :doc:`fix client/md <fix_client_md>`   # any client fix which makes specific requests to the server
-* :doc:`message quit <message>`    # terminate client/server interaction
-
-In between the two message commands, a client fix command and
-:doc:`unfix <unfix>` command can be used multiple times.  Similarly,
-this sequence of 3 commands can be repeated multiple times, assuming
-the server program operates in a similar fashion, to initiate and
-terminate client/server communication.
-
-The protocol for using LAMMPS as a server is to use these 2 commands
-in this order (other commands may come in between):
-
-* :doc:`message server <message>`  # initiate client/server interaction
-* :doc:`server md <server_md>`    # any server command which responds to specific requests from the client
-
-This sequence of 2 commands can be repeated multiple times, assuming
-the client program operates in a similar fashion, to initiate and
-terminate client/server communication.
-
-LAMMPS support for client/server coupling is in its :ref:`MESSAGE package <PKG-MESSAGE>` which implements several
-commands that enable LAMMPS to act as a client or server, as discussed
-below.  The MESSAGE package also wraps a client/server library called
-CSlib which enables two codes to exchange messages in different ways,
-either via files, sockets, or MPI.  The CSlib is provided with LAMMPS
-in the lib/message dir.  The CSlib has its own
-`website <https://cslib.sandia.gov>`_ with documentation and test
-programs.
-
-.. note::
-
-   For client/server coupling to work between LAMMPS and another
-   code, the other code also has to use the CSlib.  This can sometimes be
-   done without any modifications to the other code by simply wrapping it
-   with a Python script that exchanges CSlib messages with LAMMPS and
-   prepares input for or processes output from the other code.  The other
-   code also has to implement a matching protocol for the format and
-   content of messages that LAMMPS exchanges with it.
-
-These are the commands currently in the MESSAGE package for two
-protocols, MD and MC (Monte Carlo).  New protocols can easily be
-defined and added to this directory, where LAMMPS acts as either the
-client or server.
-
-* :doc:`message <message>`
-* :doc:`fix client md <fix_client_md>` = LAMMPS is a client for running MD
-* :doc:`server md <server_md>` = LAMMPS is a server for computing MD forces
-* :doc:`server mc <server_mc>` = LAMMPS is a server for computing a Monte Carlo energy
-
-The server doc files give details of the message protocols
-for data that is exchanged between the client and server.
-
-These example directories illustrate how to use LAMMPS as either a
-client or server code:
-
-* examples/message
-* examples/COUPLE/README
-* examples/COUPLE/lammps_mc
-* examples/COUPLE/lammps_nwchem
-* examples/COUPLE/lammps_vasp
-
-The examples/message directory couples a client instance of LAMMPS to a
-server instance of LAMMPS.
-
-The files in the *lammps_mc* folder show how to couple LAMMPS as
-a server to a simple Monte Carlo client code as the driver.
-
-The files in the *lammps_nwchem* folder show how to couple LAMMPS
-as a client code running MD timestepping to NWChem acting as a
-server providing quantum DFT forces, through a Python wrapper script
-on NWChem.
-
-The files in the *lammps_vasp* folder show how to couple LAMMPS as
-a client code running MD timestepping to VASP acting as a server
-providing quantum DFT forces, through a Python wrapper script on VASP.
-
-Here is how to launch a client and server code together for any of the
-4 modes of message exchange that the :doc:`message <message>` command
-and the CSlib support.  Here LAMMPS is used as both the client and
-server code.  Another code could be substituted for either.
-
-The examples below show launching both codes from the same window (or
-batch script), using the "&" character to launch the first code in the
-background.  For all modes except *mpi/one*, you could also launch the
-codes in separate windows on your desktop machine.  It does not
-matter whether you launch the client or server first.
-
-In these examples either code can be run on one or more processors.
-If running in a non-MPI mode (file or zmq) you can launch a code on a
-single processor without using mpirun.
-
-IMPORTANT: If you run in mpi/two mode, you must launch both codes via
-mpirun, even if one or both of them runs on a single processor.  This
-is so that MPI can figure out how to connect both MPI processes
-together to exchange MPI messages between them.
-
-For message exchange in *file*, *zmq*, or *mpi/two* modes:
-
-.. code-block:: bash
-
-   % mpirun -np 1 lmp_mpi -log log.client < in.client &
-   % mpirun -np 2 lmp_mpi -log log.server < in.server
-
-   % mpirun -np 4 lmp_mpi -log log.client < in.client &
-   % mpirun -np 1 lmp_mpi -log log.server < in.server
-
-   % mpirun -np 2 lmp_mpi -log log.client < in.client &
-   % mpirun -np 4 lmp_mpi -log log.server < in.server
-
-For message exchange in *mpi/one* mode:
-
-Launch both codes in a single mpirun command:
-
-.. code-block:: bash
-
-   mpirun -np 2 lmp_mpi -mpicolor 0 -in in.message.client -log log.client : -np 4 lmp_mpi -mpicolor 1 -in in.message.server -log log.server
-
-The two -np values determine how many procs the client and the server
-run on.
-
-A LAMMPS executable run in this manner must use the -mpicolor color
-command-line option as their its option, where color is an integer
-label that will be used to distinguish one executable from another in
-the multiple executables that the mpirun command launches.  In this
-example the client was colored with a 0, and the server with a 1.

doc/src/Howto_couple.rst

@@ -12,16 +12,16 @@ LAMMPS can be coupled to other codes in at least 4 ways.  Each has
 advantages and disadvantages, which you will have to think about in the
 context of your application.
 
-1. Define a new :doc:`fix <fix>` command that calls the other code.  In
-   this scenario, LAMMPS is the driver code.  During timestepping,
+1. Define a new :doc:`fix <fix>` command that calls the other code.
+   In this scenario, LAMMPS is the driver code.  During timestepping,
    the fix is invoked, and can make library calls to the other code,
-   which has been linked to LAMMPS as a library.  This is the way how the
+   which has been linked to LAMMPS as a library.  This is the way the
    :ref:`LATTE <PKG-LATTE>` package, which performs density-functional
-   tight-binding calculations using the `LATTE software <https://github.com/lanl/LATTE>`_
-   to compute forces, is hooked to LAMMPS.
-   See the :doc:`fix latte <fix_latte>` command for more details.
-   Also see the :doc:`Modify <Modify>` doc pages for info on how to
-   add a new fix to LAMMPS.
+   tight-binding calculations using the `LATTE software
+   <https://github.com/lanl/LATTE>`_ to compute forces, is hooked to
+   LAMMPS.  See the :doc:`fix latte <fix_latte>` command for more
+   details.  Also see the :doc:`Modify <Modify>` doc pages for info on
+   how to add a new fix to LAMMPS.
 
 .. spacer
 
@@ -58,6 +58,12 @@ context of your application.
 
 .. spacer
 
-4. Couple LAMMPS with another code in a client/server mode.  This is
-   described on the :doc:`Howto client/server <Howto_client_server>` doc
-   page.
+4. Couple LAMMPS with another code in a client/server fashion, using
+   using the `MDI Library
+   <https://molssi-mdi.github.io/MDI_Library/html/index.html>`_
+   developed by the `Molecular Sciences Software Institute (MolSSI)
+   <https://molssi.org>`_ to run LAMMPS as either an MDI driver
+   (client) or an MDI engine (server).  The MDI driver issues commands
+   to the MDI server to exchange data between them.  See the
+   :doc:`Howto mdi <Howto_mdi>` page for more information about how
+   LAMMPS can operate in either of these modes.

doc/src/Howto_elastic.rst

@@ -18,23 +18,52 @@ At zero temperature, it is easy to estimate these derivatives by
 deforming the simulation box in one of the six directions using the
 :doc:`change_box <change_box>` command and measuring the change in the
 stress tensor. A general-purpose script that does this is given in the
-examples/elastic directory described on the :doc:`Examples <Examples>`
+examples/ELASTIC directory described on the :doc:`Examples <Examples>`
 doc page.
 
 Calculating elastic constants at finite temperature is more
 challenging, because it is necessary to run a simulation that performs
-time averages of differential properties. One way to do this is to
-measure the change in average stress tensor in an NVT simulations when
+time averages of differential properties. There are at least
+3 ways to do this in LAMMPS. The most reliable way to do this is
+by exploiting the relationship between elastic constants, stress
+fluctuations, and the Born matrix, the second derivatives of energy
+w.r.t. strain :ref:`(Ray) <Ray>`.
+The Born matrix calculation has been enabled by
+the :doc:`compute born/matrix <compute_born_matrix>` command,
+which works for any bonded or non-bonded potential in LAMMPS.
+The most expensive part of the calculation is the sampling of
+the stress fluctuations. Several examples of this method are
+provided in the examples/ELASTIC_T/BORN_MATRIX directory
+described on the :doc:`Examples <Examples>` doc page.
+
+A second way is to measure
+the change in average stress tensor in an NVT simulations when
 the cell volume undergoes a finite deformation. In order to balance
 the systematic and statistical errors in this method, the magnitude of
 the deformation must be chosen judiciously, and care must be taken to
 fully equilibrate the deformed cell before sampling the stress
-tensor. Another approach is to sample the triclinic cell fluctuations
+tensor. An example of this method is provided in the
+examples/ELASTIC_T/DEFORMATION directory
+described on the :doc:`Examples <Examples>` doc page.
+
+Another approach is to sample the triclinic cell fluctuations
 that occur in an NPT simulation. This method can also be slow to
-converge and requires careful post-processing :ref:`(Shinoda) <Shinoda1>`
+converge and requires careful post-processing :ref:`(Shinoda) <Shinoda1>`.
+We do not provide an example of this method.
+
+A nice review of the advantages and disadvantages of all of these methods
+is provided in the paper by Clavier et al. :ref:`(Clavier) <Clavier>`.
 
 ----------
 
+.. _Ray:
+
+**(Ray)** J. R. Ray and A. Rahman, J Chem Phys, 80, 4423 (1984).
+
 .. _Shinoda1:
 
 **(Shinoda)** Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004).
+
+.. _Clavier:
+
+**(Clavier)** G. Clavier, N. Desbiens, E. Bourasseau, V. Lachet, N. Brusselle-Dupend and B. Rousseau, Mol Sim, 43, 1413 (2017).

doc/src/Howto_mdi.rst

@@ -1,132 +1,144 @@
 Using LAMMPS with the MDI library for code coupling
 ===================================================
 
-.. note::
-
-  This Howto page will eventually replace the
-  :doc:`Howto client/server <Howto_client_server>` doc page.
-
-Client/server coupling of two codes is where one code is the "client"
-and sends request messages (data) to a "server" code.  The server
-responds to each request with a reply message.  This enables the two
-codes to work in tandem to perform a simulation.  LAMMPS can act as
-either a client or server code; it does this by using the `MolSSI
-Driver Interface (MDI) library
+Client/server coupling of two (or more) codes is where one code is the
+"client" and sends request messages (data) to one (or more) "server"
+code(s).  A server responds to each request with a reply message
+(data).  This enables two (or more) codes to work in tandem to perform
+a simulation.  LAMMPS can act as either a client or server code; it
+does this by using the `MolSSI Driver Interface (MDI) library
 <https://molssi-mdi.github.io/MDI_Library/html/index.html>`_,
 developed by the `Molecular Sciences Software Institute (MolSSI)
-<https://molssi.org>`_.
+<https://molssi.org>`_, which is supported by the :ref:`MDI <PKG-MDI>`
+package.
 
 Alternate methods for code coupling with LAMMPS are described on the
 :doc:`Howto couple <Howto_couple>` doc page.
 
-Some advantages of client/server coupling are that the two codes can run
+Some advantages of client/server coupling are that the codes can run
 as stand-alone executables; they need not be linked together.  Thus
-neither code needs to have a library interface.  This also makes it easy
-to run the two codes on different numbers of processors.  If a message
-protocol (format and content) is defined for a particular kind of
-simulation, then in principle any code which implements the client-side
-protocol can be used in tandem with any code which implements the
-server-side protocol.  Neither code needs to know what specific other
-code it is working with.
+neither code needs to have a library interface.  This also makes it
+easy to run the two codes on different numbers of processors.  If a
+message protocol (format and content) is defined for a particular kind
+of simulation, then in principle any code which implements the
+client-side protocol can be used in tandem with any code which
+implements the server-side protocol.  Neither code needs to know what
+specific other code it is working with.
 
 In MDI nomenclature, a client code is the "driver", and a server code is
 an "engine".  One driver code can communicate with one or more instances
 of one or more engine codes.  Driver and engine codes can be written in
 any language: C, C++, Fortran, Python, etc.
 
-In addition to allowing driver and engine(s) running to run as
-stand-alone executables, MDI also enables a server code to be a
-"plugin" to the client code.  In this scenario, server code(s) are
-compiled as shared libraries, and one (or more) instances of the
-server are instantiated by the driver code.  If the driver code runs
-in parallel, it can split its MPI communicator into multiple
-sub-communicators, and launch each plugin engine instance on a
-sub-communicator.  Driver processors in that sub-communicator exchange
-messages with that engine instance, and can also send MPI messages to
-other processors in the driver.  The driver code can also destroy
-engine instances and re-instantiate them.
+In addition to allowing driver and engine(s) to run as stand-alone
+executables, MDI also enables an engine to be a *plugin* to the client
+code.  In this scenario, server code(s) are compiled as shared
+libraries, and one (or more) instances of the server are instantiated
+by the driver code.  If the driver code runs in parallel, it can split
+its MPI communicator into multiple sub-communicators, and launch each
+plugin engine instance on a sub-communicator.  Driver processors
+within that sub-communicator exchange messages with the corresponding
+engine instance, and can also send MPI messages to other processors in
+the driver.  The driver code can also destroy engine instances and
+re-instantiate them.  LAMMPS can operate as either a stand-alone or
+plugin MDI engine.  When it operates as a driver, if can use either
+stand-alone or plugin MDI engines.
 
-The way that a driver communicates with an engine is by making
-MDI_Send() and MDI_Recv() calls, which are conceptually similar to
-MPI_Send() and MPI_Recv() calls.  Each send or receive has a string
-which identifies the command name, and optionally some data, which can
-be a single value or vector of values of any data type.  Inside the
-MDI library, data is exchanged between the driver and engine via MPI
-calls or sockets.  This a run-time choice by the user.
+The way in which an MDI driver communicates with an MDI engine is by
+making MDI_Send() and MDI_Recv() calls, which are conceptually similar
+to MPI_Send() and MPI_Recv() calls.  Each send or receive operation
+uses a string to identify the command name, and optionally some data,
+which can be a single value or vector of values of any data type.
+Inside the MDI library, data is exchanged between the driver and
+engine via MPI calls or sockets.  This a run-time choice by the user.
+
+----------
+
+The :ref:`MDI <PKG-MDI>` package provides a :doc:`mdi engine <mdi>`
+command which enables LAMMPS to operate as an MDI engine.  Its doc
+page explains the variety of standard and custom MDI commands which
+the LAMMPS engine recognizes and can respond to.
+
+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 also has a `fix mdi/aimd <fix_mdi_aimd>` command in which
+LAMMPS operates as an MDI driver to perform *ab initio* MD simulations
+in conjunction with a quantum mechanics code.  Its post_force() method
+illustrates how a driver issues MDI commands to another code.  This
+command can be used to couple to an MDI engine which is either a
+stand-alone code or a plugin library.
+
+----------
+
+The examples/mdi directory contains Python scripts and LAMMPS input
+script which use LAMMPS as either an MDI driver or engine or both.
+Three example use cases are provided:
+
+* Run ab initio MD (AIMD) using 2 instances of LAMMPS, one as driver
+  and one as an engine.  As an engine, LAMMPS is a surrogate for a
+  quantum code.
+
+* A Python script driver invokes a sequence of unrelated LAMMPS
+  calculations.  Calculations can be single-point energy/force
+  evaluations, MD runs, or energy minimizations.
+
+* Run AIMD with a Python driver code and 2 LAMMPS instances as
+  engines.  The first LAMMPS instance performs MD timestepping.  The
+  second LAMMPS instance acts as a surrogate QM code to compute
+  forces.
+
+Note that in any of these example where LAMMPS is used as an engine,
+an actual QM code (which supports MDI) could be used in its place,
+without modifying other code or scripts, except to specify the name of
+the QM code.
+
+The examples/mdi/README file explains how to launch both driver and
+engine codes so that they communicate using the MDI library via either
+MPI or sockets.
 
 -------------
 
-As an example, LAMMPS and the ``pw.x`` command from Quantum Espresso (a
-suite of quantum DFT codes), can work together via the MDI library to
-perform an ab initio MD (AIMD) simulation, where LAMMPS runs an MD
-simulation and sends a message each timestep to ``pw.x`` asking it to
-compute quantum forces on the current configuration of atoms.  Here is
-how the 2 codes are launched to communicate by MPI:
+Currently there are two quantum DFT 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 QM codes which have indirect support through QCEngine or i-PI.
+The former means they require a wrapper program (QCEngine) with MDI
+support which writes/read files to pass data to the quantum code
+itself.  The list of QCEngine-supported and i-PI-supported quantum
+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
+used in stand-alone mode:
 
 .. code-block:: bash
 
-   % mpirun -np 2 lmp_mpi -mdi "-role DRIVER -name d -method MPI" \
-     -in in.aimd : -np 16 pw.x -in qe.in -mdi "-role ENGINE -name e -method MPI"
+   % 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>`_
 
-In this case LAMMPS runs on 2 processors (MPI tasks), ``pw.x`` runs on 16
-processors.
-
-Here is how the 2 codes are launched to communicate by sockets:
+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
 
-   % mpirun -np 2 lmp_mpi -mdi "-role DRIVER -name d -method TCP -port 8021" -in in.aimd
-   % mpirun -np 16 pw.x -in qe.in -mdi "-role ENGINE -name e -method TCP -port 8021 -hostname localhost"
+   % 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
 
-These commands could be issued in different windows on a desktop
-machine.  Or in the same window, if the first command is ended with
-"&" so as to run in the background.  If "localhost" is replaced by an
-IP address, ``pw.x`` could be run on another machine on the same network, or
-even on another machine across the country.
+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:
 
-After both codes initialize themselves to model the same system, this is
-what occurs each timestep:
+.. code-block:: bash
 
-* LAMMPS send a ">COORDS" message to ``pw.x`` with a 3*N vector of current atom coords
-* ``pw.x`` receives the message/coords and computes quantum forces on all the atoms
-* LAMMPS send a "<FORCES" message to ``pw.x`` and waits for the result
-* ``pw.x`` receives the message (after its computation finishes) and sends a 3*N vector of forces
-* LAMMPS receives the forces and time integrates to complete a single timestep
-
--------------
-
-Examples scripts for using LAMMPS as an MDI engine are in the
-examples/mdi directory.  See the README file in that directory for
-instructions on how to run the examples.
-
-.. note::
-
-  Work is underway to add commands that allow LAMMPS to be used as an
-  MDI driver, e.g. for the AIMD example discussed above.  Example
-  scripts for this usage mode will be added the same directory when
-  available.
-
-If LAMMPS is used as a stand-alone engine it should set up the system
-it will be modeling in its input script, then invoke the
-:doc:`mdi/engine <mdi_engine>` command.  This will put LAMMPS into
-"engine mode" where it waits for messages and data from the driver.
-When the driver sends an "EXIT" command, LAMMPS will exit engine mode
-and the input script will continue.
-
-If LAMMPS is used as a plugin engine it operates the same way, except
-that the driver will pass LAMMPS an input script to initialize itself.
-Upon receiving the "EXIT" command, LAMMPS will exit engine mode and the
-input script will continue.  After finishing execution of the input
-script, the instance of LAMMPS will be destroyed.
-
-LAMMPS supports the full set of MD-appropriate engine commands defined
-by the MDI library.  See the :doc:`mdi/engine <mdi_engine>` page for
-a list of these.
-
-If those commands are not sufficient for a user-developed driver to use
-LAMMPS as an engine, then new commands can be easily added.  See these
-two files which implement the definition of MDI commands and the logic
-for responding to them:
-
-* src/MDI/mdi_engine.cpp
-* src/MDI/fix_mdi_engine.cpp
+   % 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

@@ -19,7 +19,7 @@ atoms and the water molecule to run a rigid SPC model.
 | LJ :math:`\sigma` of OO = 3.166
 | LJ :math:`\epsilon`, :math:`\sigma` of OH, HH = 0.0
 | :math:`r_0` of OH bond = 1.0
-| :math:`\theta` of HOH angle = 109.47\ :math:`^{\circ}`
+| :math:`\theta_0` of HOH angle = 109.47\ :math:`^{\circ}`
 |
 
 Note that as originally proposed, the SPC model was run with a 9

doc/src/Howto_structured_data.rst

@@ -55,6 +55,9 @@ JSON
 YAML format thermo_style output
 ===============================
 
+Extracting data from log file
+-----------------------------
+
 .. versionadded:: 24Mar2022
 
 LAMMPS supports the thermo style "yaml" and for "custom" style
@@ -66,7 +69,7 @@ the following style:
 .. code-block:: yaml
 
    ---
-   keywords: [Step, Temp, E_pair, E_mol, TotEng, Press, ]
+   keywords: ['Step', 'Temp', 'E_pair', 'E_mol', 'TotEng', 'Press', ]
    data:
      - [100, 0.757453103239935, -5.7585054860159, 0, -4.62236133677021, 0.207261053624721, ]
      - [110, 0.759322359337036, -5.7614668389562, 0, -4.62251889318624, 0.194314975399602, ]
@@ -79,6 +82,10 @@ This data can be extracted and parsed from a log file using python with:
 .. code-block:: python
 
    import re, yaml
+   try:
+       from yaml import CSafeLoader as Loader
+   except ImportError:
+       from yaml import SafeLoader as Loader
 
    docs = ""
    with open("log.lammps") as f:
@@ -86,7 +93,7 @@ This data can be extracted and parsed from a log file using python with:
            m = re.search(r"^(keywords:.*$|data:$|---$|\.\.\.$|  - \[.*\]$)", line)
            if m: docs += m.group(0) + '\n'
 
-   thermo = list(yaml.load_all(docs, Loader=yaml.SafeLoader))
+   thermo = list(yaml.load_all(docs, Loader=Loader))
 
    print("Number of runs: ", len(thermo))
    print(thermo[1]['keywords'][4], ' = ', thermo[1]['data'][2][4])
@@ -105,6 +112,135 @@ of that run:
    Number of runs:  2
    TotEng  =  -4.62140097780047
 
+.. versionadded:: 4May2022
+
+YAML format output has been added to multiple commands in LAMMPS,
+for example :doc:`dump yaml <dump>` or :doc:`fix ave/time <fix_ave_time>`
+Depending on the kind of data being written, organization of the data
+or the specific syntax used may change, but the principles are very
+similar and all files should be readable with a suitable YAML parser.
+
+Processing scalar data with Python
+----------------------------------
+
+.. figure:: JPG/thermo_bondeng.png
+            :figwidth: 33%
+            :align: right
+
+After reading and parsing the YAML format data, it can be easily
+imported for further processing and visualization with the `pandas
+<https://pandas.pydata.org/>`_ and `matplotlib
+<https://matplotlib.org/>`_ Python modules.  Because of the organization
+of the data in the YAML format thermo output, it needs to be told to
+process only the 'data' part of the imported data to create a pandas
+data frame, and one needs to set the column names from the 'keywords'
+entry.  The following Python script code example demonstrates this, and
+creates the image shown on the right of a simple plot of various bonded
+energy contributions versus the timestep from a run of the 'peptide'
+example input after changing the :doc:`thermo style <thermo_style>` to
+'yaml'.  The properties to be used for x and y values can be
+conveniently selected through the keywords.  Please note that those
+keywords can be changed to custom strings with the :doc:`thermo_modify
+colname <thermo_modify>` command.
+
+.. code-block:: python
+
+   import re, yaml
+   import pandas as pd
+   import matplotlib.pyplot as plt
+
+   try:
+       from yaml import CSafeLoader as Loader
+   except ImportError:
+       from yaml import SafeLoader as Loader
+
+   docs = ""
+   with open("log.lammps") as f:
+       for line in f:
+           m = re.search(r"^(keywords:.*$|data:$|---$|\.\.\.$|  - \[.*\]$)", line)
+           if m: docs += m.group(0) + '\n'
+
+   thermo = list(yaml.load_all(docs, Loader=Loader))
+
+   df = pd.DataFrame(data=thermo[0]['data'], columns=thermo[0]['keywords'])
+   fig = df.plot(x='Step', y=['E_bond', 'E_angle', 'E_dihed', 'E_impro'], ylabel='Energy in kcal/mol')
+   plt.savefig('thermo_bondeng.png')
+
+Processing vector data with Python
+----------------------------------
+
+Global *vector* data as produced by :doc:`fix ave/time <fix_ave_time>`
+uses a slightly different organization of the data. You still have the
+dictionary keys 'keywords' and 'data' for the column headers and the
+data.  But the data is a dictionary indexed by the time step and for
+each step there are multiple rows of values each with a list of the
+averaged properties.  This requires a slightly different processing,
+since the entire data cannot be directly imported into a single pandas
+DataFrame class instance.  The following Python script example
+demonstrates how to read such data.  The result will combine the data
+for the different steps into one large "multi-index" table.  The pandas
+IndexSlice class can then be used to select data from this combined data
+frame.
+
+.. code-block:: python
+
+   import re, yaml
+   import pandas as pd
+
+   try:
+       from yaml import CSafeLoader as Loader
+   except ImportError:
+       from yaml import SafeLoader as Loader
+
+   with open("ave.yaml") as f:
+       ave = yaml.load(docs, Loader=Loader)
+
+   keys = ave['keywords']
+   df = {}
+   for k in ave['data'].keys():
+       df[k] = pd.DataFrame(data=ave['data'][k], columns=keys)
+
+   # create multi-index data frame
+   df = pd.concat(df)
+
+   # output only the first 3 value for steps 200 to 300 of the column Pressure
+   idx = pd.IndexSlice
+   print(df['Pressure'].loc[idx[200:300, 0:2]])
+
+
+Processing scalar data with Perl
+--------------------------------
+
+The ease of processing YAML data is not limited to Python. Here is an
+example for extracting and processing a LAMMPS log file with Perl instead.
+
+.. code-block:: perl
+
+   use YAML::XS;
+
+   open(LOG, "log.lammps") or die("could not open log.lammps: $!");
+   my $file = "";
+   while(my $line = <LOG>) {
+       if ($line =~ /^(keywords:.*$|data:$|---$|\.\.\.$|  - \[.*\]$)/) {
+           $file .= $line;
+       }
+   }
+   close(LOG);
+
+   # convert YAML to perl as nested hash and array references
+   my $thermo = Load $file;
+
+   # convert references to real arrays
+   my @keywords = @{$thermo->{'keywords'}};
+   my @data = @{$thermo->{'data'}};
+
+   # print first two columns
+   print("$keywords[0] $keywords[1]\n");
+   foreach (@data) {
+       print("${$_}[0]  ${$_}[1]\n");
+   }
+
+
 Writing continuous data during a simulation
 ===========================================
 

doc/src/Install_mac.rst

@@ -5,7 +5,7 @@ LAMMPS can be downloaded, built, and configured for OS X on a Mac with
 `Homebrew <homebrew_>`_.  (Alternatively, see the install instructions for
 :doc:`Download an executable via Conda <Install_conda>`.)  The following LAMMPS
 packages are unavailable at this time because of additional needs not yet met:
-GPU, KOKKOS, LATTE, MSCG, MESSAGE, MPIIO POEMS VORONOI.
+GPU, KOKKOS, LATTE, MSCG, MPIIO, POEMS, VORONOI.
 
 After installing Homebrew, you can install LAMMPS on your system with
 the following commands:

doc/src/Install_tarball.rst

@@ -6,7 +6,7 @@ of the `LAMMPS website <lws_>`_.
 
 .. _download: https://www.lammps.org/download.html
 .. _bug: https://www.lammps.org/bug.html
-.. _older: https://www.lammps.org/tars
+.. _older: https://download.lammps.org/tars/
 .. _lws: https://www.lammps.org
 
 You have two choices of tarballs, either the most recent stable

doc/src/Intro_citing.rst

@@ -46,7 +46,7 @@ In addition there are DOIs for individual stable releases. Currently there are:
 
 - 3 March 2020 version: `DOI:10.5281/zenodo.3726417 <https://dx.doi.org/10.5281/zenodo.3726417>`_
 - 29 October 2020 version: `DOI:10.5281/zenodo.4157471 <https://dx.doi.org/10.5281/zenodo.4157471>`_
-
+- 29 September 2021 version: `DOI:10.5281/zenodo.6386596 <https//dx.doi.org/10.5281/zenodo.6386596>`_
 
 Home page
 ^^^^^^^^^

doc/src/Intro_website.rst

@@ -20,7 +20,6 @@ available online are listed below.
 * `Glossary of terms relevant to LAMMPS <https://www.lammps.org/glossary.html>`_
 * `LAMMPS highlights with images <https://www.lammps.org/pictures.html>`_
 * `LAMMPS highlights with movies <https://www.lammps.org/movies.html>`_
-* `Mailing list <https://www.lammps.org/mail.html>`_
 * `LAMMPS forum <https://www.lammps.org/forum.html>`_
 * `Workshops <https://www.lammps.org/workshops.html>`_
 * `Tutorials <https://www.lammps.org/tutorials.html>`_

doc/src/Manual.rst

@@ -14,11 +14,11 @@ LAMMPS is an open-source code, distributed freely under the terms of
 the GNU Public License Version 2 (GPLv2).
 
 The `LAMMPS website <lws_>`_ has a variety of information about the
-code.  It includes links to an on-line version of this manual, a
-`mailing list <https://www.lammps.org/mail.html>`_ and
-`online forum <https://www.lammps.org/forum.html>`_  where users can
-post questions, and a `GitHub site <https://github.com/lammps/lammps>`_
-where all LAMMPS development is coordinated.
+code.  It includes links to an on-line version of this manual, an
+`online forum <https://www.lammps.org/forum.html>`_ where users can post
+questions and discuss LAMMPS, and a `GitHub site
+<https://github.com/lammps/lammps>`_ where all LAMMPS development is
+coordinated.
 
 ----------
 

doc/src/Modify_atom.rst

@@ -9,34 +9,34 @@ A new atom style can be created if one of the existing atom styles
 does not define all the attributes you need to store and communicate
 with atoms.
 
-Atom_vec_atomic.cpp is the simplest example of an atom style.
+The file ``atom_vec_atomic.cpp`` is the simplest example of an atom style.
 Examining the code for others will make these instructions more clear.
 
-Note that the :doc:`atom style hybrid <atom_style>` command can be
-used to define atoms or particles which have the union of properties
-of individual styles.  Also the :doc:`fix property/atom <fix_property_atom>`
-command can be used to add a single property (e.g. charge
-or a molecule ID) to a style that does not have it.  It can also be
-used to add custom properties to an atom, with options to communicate
-them with ghost atoms or read them from a data file.  Other LAMMPS
-commands can access these custom properties, as can new pair, fix,
-compute styles that are written to work with these properties.  For
+Note that the :doc:`atom style hybrid <atom_style>` command can be used
+to define atoms or particles which have the union of properties of
+individual styles.  Also the :doc:`fix property/atom
+<fix_property_atom>` command can be used to add a single property
+(e.g. charge or a molecule ID) to a style that does not have it.  It can
+also be used to add custom properties to an atom, with options to
+communicate them with ghost atoms or read them from a data file.  Other
+LAMMPS commands can access these custom properties, as can new pair,
+fix, compute styles that are written to work with these properties.  For
 example, the :doc:`set <set>` command can be used to set the values of
 custom per-atom properties from an input script.  All of these methods
-are less work than writing code for a new atom style.
+are less work than writing and testing(!) code for a new atom style.
 
 If you follow these directions your new style will automatically work
 in tandem with others via the :doc:`atom_style hybrid <atom_style>`
 command.
 
-The first step is to define a set of strings in the constructor of the
-new derived class.  Each string will have zero or more space-separated
-variable names which are identical to those used in the atom.h header
-file for per-atom properties.  Note that some represent per-atom
+The first step is to define a set of string lists in the constructor of
+the new derived class.  Each list will have zero or more comma-separated
+strings that correspond to the variable names used in the ``atom.h``
+header file for per-atom properties.  Note that some represent per-atom
 vectors (q, molecule) while other are per-atom arrays (x,v).  For all
-but the last 2 strings you do not need to specify any of
+but the last two lists you do not need to specify any of
 (id,type,x,v,f).  Those are included automatically as needed in the
-other strings.
+other lists.
 
 .. list-table::
 
@@ -65,16 +65,16 @@ other strings.
    * - fields_data_vel
      - list of properties (in order) in the Velocities section of a data file, as read by :doc:`read_data <read_data>`
 
-In these strings you can list variable names which LAMMPS already
-defines (in some other atom style), or you can create new variable
-names.  You should not re-use a LAMMPS variable for something with
-different meaning in your atom style.  If the meaning is related, but
-interpreted differently by your atom style, then using the same
-variable name means a user should not use your style and the other
-style together in a :doc:`atom_style hybrid <atom_style>` command.
-Because there will only be one value of the variable and different
-parts of LAMMPS will then likely use it differently.  LAMMPS has
-no way of checking for this.
+In these lists you can list variable names which LAMMPS already defines
+(in some other atom style), or you can create new variable names.  You
+should not re-use a LAMMPS variable in your atom style that is used for
+something with a different meaning in another atom style.  If the
+meaning is related, but interpreted differently by your atom style, then
+using the same variable name means a user must not use your style and
+the other style together in a :doc:`atom_style hybrid <atom_style>`
+command.  Because there will only be one value of the variable and
+different parts of LAMMPS will then likely use it differently.  LAMMPS
+has no way of checking for this.
 
 If you are defining new variable names then make them descriptive and
 unique to your new atom style.  For example choosing "e" for energy is
@@ -85,32 +85,31 @@ If any of the variable names in your new atom style do not exist in
 LAMMPS, you need to add them to the src/atom.h and atom.cpp files.
 
 Search for the word "customize" or "customization" in these 2 files to
-see where to add your variable.  Adding a flag to the 2nd
-customization section in atom.h is only necessary if your code (e.g. a
-pair style) needs to check that a per-atom property is defined.  These
-flags should also be set in the constructor of the atom style child
-class.
+see where to add your variable.  Adding a flag to the 2nd customization
+section in ``atom.h`` is only necessary if your code (e.g. a pair style)
+needs to check that a per-atom property is defined.  These flags should
+also be set in the constructor of the atom style child class.
 
-In atom.cpp, aside from the constructor and destructor, there are 3
+In ``atom.cpp``, aside from the constructor and destructor, there are 3
 methods that a new variable name or flag needs to be added to.
 
-In Atom::peratom_create() when using the add_peratom() method, a
-final length argument of 0 is for per-atom vectors, a length > 1 is
-for per-atom arrays.  Note the use of an extra per-thread flag and the
-add_peratom_vary() method when last dimension of the array is
+In ``Atom::peratom_create()`` when using the ``Atom::add_peratom()``
+method, a cols argument of 0 is for per-atom vectors, a length >
+1 is for per-atom arrays.  Note the use of the extra per-thread flag and
+the add_peratom_vary() method when last dimension of the array is
 variable-length.
 
-Adding the variable name to Atom::extract() enable the per-atom data
+Adding the variable name to Atom::extract() enables the per-atom data
 to be accessed through the :doc:`LAMMPS library interface
 <Howto_library>` by a calling code, including from :doc:`Python
 <Python_head>`.
 
 The constructor of the new atom style will also typically set a few
-flags which are defined at the top of atom_vec.h.  If these are
+flags which are defined at the top of ``atom_vec.h``.  If these are
 unclear, see how other atom styles use them.
 
-The grow_pointers() method is also required to make
-a copy of peratom data pointers, as explained in the code.
+The grow_pointers() method is also required to make a copy of peratom
+data pointers, as explained in the code.
 
 There are a number of other optional methods which your atom style can
 implement.  These are only needed if you need to do something

doc/src/Modify_contribute.rst

@@ -27,11 +27,10 @@ join the `LAMMPS developers on Slack <https://lammps.slack.com>`_.  This
 slack work space is by invitation only.  Thus for access, please send an
 e-mail to ``slack@lammps.org`` explaining what part of LAMMPS you are
 working on.  Only discussions related to LAMMPS development are
-tolerated in that work space, so this is **NOT** for people that look for
-help with compiling, installing, or using LAMMPS.  Please post a message
-to the `lammps-users mailing list <https://www.lammps.org/mail.html>`_
-or the `LAMMPS forum <https://www.lammps.org/forum.html>`_ for those
-purposes.
+tolerated in that work space, so this is **NOT** for people that look
+for help with compiling, installing, or using LAMMPS.  Please post a
+message to the `LAMMPS forum <https://www.lammps.org/forum.html>`_ for
+those purposes.
 
 Packages versus individual files
 --------------------------------

doc/src/Modify_style.rst

@@ -223,6 +223,13 @@ and readable by all and no executable permissions.  Executable
 permissions (0755) should only be on shell scripts or python or similar
 scripts for interpreted script languages.
 
+You can check for these issues with the python scripts in the
+:ref:`"tools/coding_standard" <coding_standard>` folder.  When run
+normally with a source file or a source folder as argument, they will
+list all non-conforming lines.  By adding the `-f` flag to the command
+line, they will modify the flagged files to try removing the detected
+issues.
+
 Indentation and Placement of Braces (strongly preferred)
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -240,6 +247,53 @@ reformatting from clang-format yields undesirable output may be
 protected with placing a pair `// clang-format off` and `// clang-format
 on` comments around that block.
 
+Error or warning messages and explanations (preferred)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. versionchanged:: 4May2022
+
+Starting with LAMMPS version 4 May 2022 the LAMMPS developers have
+agreed on a new policy for error and warning messages.
+
+Previously, all error and warning strings were supposed to be listed in
+the class header files with an explanation.  Those would then be
+regularly "harvested" and transferred to alphabetically sorted lists in
+the manual.  To avoid excessively long lists and to reduce effort, this
+came with a requirement to have rather generic error messages (e.g.
+"Illegal ... command").  To identify the specific cause, the name of the
+source file and the line number of the error location would be printed,
+so that one could look up the cause by reading the source code.
+
+The new policy encourages more specific error messages that ideally
+indicate the cause directly and no further lookup would be needed.
+This is aided by using the `{fmt} library <https://fmt.dev>`_ to convert
+the Error class commands so that they take a variable number of arguments
+and error text will be treated like a {fmt} syntax format string.
+Error messages should still kept to a single line or two lines at the most.
+
+For more complex explanations or errors that have multiple possible
+reasons, a paragraph should be added to the `Error_details` page with an
+error code reference (e.g. ``.. _err0001:``) then the utility function
+:cpp:func:`utils::errorurl() <LAMMPS_NS::utils::errorurl>` can be used
+to generate an URL that will directly lead to that paragraph.  An error
+for missing arguments can be easily generated using the
+:cpp:func:`utils::missing_cmd_args()
+<LAMMPS_NS::utils::missing_cmd_args>` convenience function.
+
+The transformation of existing LAMMPS code to this new scheme is ongoing
+and - given the size of the LAMMPS source code - will take a significant
+amount of time until completion.  However, for new code following the
+new approach is strongly preferred.  The expectation is that the new
+scheme will make it easier for LAMMPS users, developers, and
+maintainers.
+
+An example for this approach would be the
+``src/read_data.cpp`` and ``src/atom.cpp`` files that implement the
+:doc:`read_data <read_data>` and :doc:`atom_modify <atom_modify>`
+commands and that may create :ref:`"Unknown identifier in data file" <err0001>`
+errors that seem difficult to debug for users because they may have
+one of multiple possible reasons, and thus require some additional explanations.
+
 Programming language standards (required)
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

doc/src/Packages_details.rst

@@ -32,6 +32,7 @@ page gives those details.
    * :ref:`AWPMD <PKG-AWPMD>`
    * :ref:`BOCS <PKG-BOCS>`
    * :ref:`BODY <PKG-BODY>`
+   * :ref:`BPM <PKG-BPM>`
    * :ref:`BROWNIAN <PKG-BROWNIAN>`
    * :ref:`CG-DNA <PKG-CG-DNA>`
    * :ref:`CG-SDK <PKG-CG-SDK>`
@@ -49,6 +50,7 @@ page gives those details.
    * :ref:`DPD-SMOOTH <PKG-DPD-SMOOTH>`
    * :ref:`DRUDE <PKG-DRUDE>`
    * :ref:`EFF <PKG-EFF>`
+   * :ref:`ELECTRODE <PKG-ELECTRODE>`
    * :ref:`EXTRA-COMPUTE <PKG-EXTRA-COMPUTE>`
    * :ref:`EXTRA-DUMP <PKG-EXTRA-DUMP>`
    * :ref:`EXTRA-FIX <PKG-EXTRA-FIX>`
@@ -72,7 +74,6 @@ page gives those details.
    * :ref:`MDI <PKG-MDI>`
    * :ref:`MEAM <PKG-MEAM>`
    * :ref:`MESONT <PKG-MESONT>`
-   * :ref:`MESSAGE <PKG-MESSAGE>`
    * :ref:`MGPT <PKG-MGPT>`
    * :ref:`MISC <PKG-MISC>`
    * :ref:`ML-HDNNP <PKG-ML-HDNNP>`
@@ -284,6 +285,33 @@ overview.
 
 ----------
 
+.. _PKG-BPM:
+
+BPM package
+------------
+
+**Contents:**
+
+Pair styles, bond styles, fixes, and computes for bonded particle
+models for mesoscale simulations of solids and fracture.  See the
+:doc:`Howto bpm <Howto_bpm>` page for an overview.
+
+**Authors:** Joel T. Clemmer (Sandia National Labs)
+
+**Supporting info:**
+
+* src/BPM filenames -> commands
+* :doc:`Howto_bpm <Howto_bpm>`
+* :doc:`atom_style bpm/sphere <atom_style>`
+* :doc:`bond_style bpm/rotational <bond_bpm_rotational>`
+* :doc:`bond_style bpm/spring <bond_bpm_spring>`
+* :doc:`compute nbond/atom <compute_nbond_atom>`
+* :doc:`fix nve/bpm/sphere <fix_nve_bpm_sphere>`
+* :doc:`pair_style bpm/spring <pair_bpm_spring>`
+* examples/bpm
+
+----------
+
 .. _PKG-BROWNIAN:
 
 BROWNIAN package
@@ -529,8 +557,20 @@ To use this package, also the :ref:`KSPACE <PKG-KSPACE>` and
 **Supporting info:**
 
 * src/DIELECTRIC: filenames -> commands
+* :doc:`atom_style dielectric <atom_style>`
+* :doc:`pair_style coul/cut/dielectric <pair_dielectric>`
+* :doc:`pair_style coul/long/dielectric <pair_dielectric>`
+* :doc:`pair_style lj/cut/coul/cut/dielectric <pair_dielectric>`
+* :doc:`pair_style lj/cut/coul/debye/dielectric <pair_dielectric>`
+* :doc:`pair_style lj/cut/coul/long/dielectric <pair_dielectric>`
+* :doc:`pair_style lj/cut/coul/msm/dielectric <pair_dielectric>`
+* :doc:`pair_style pppm/dielectric <kspace_style>`
+* :doc:`pair_style pppm/disp/dielectric <kspace_style>`
+* :doc:`pair_style msm/dielectric <kspace_style>`
+* :doc:`fix_style polarize/bem/icc <fix_polarize>`
+* :doc:`fix_style polarize/bem/gmres <fix_polarize>`
+* :doc:`fix_style polarize/functional <fix_polarize>`
 * :doc:`compute efield/atom  <compute_efield_atom>`
-* TODO: add all styles
 * examples/PACKAGES/dielectric
 
 ----------
@@ -770,6 +810,31 @@ tools/eff; see its README file.
 
 -------------------
 
+.. _PKG-ELECTRODE:
+
+ELECTRODE package
+-----------------
+
+**Contents:**
+
+The ELECTRODE package allows the user to enforce a constant potential method for
+groups of atoms that interact with the remaining atoms as electrolyte.
+
+**Authors:** The ELECTRODE library is written and maintained by Ludwig
+Ahrens-Iwers (TUHH, Hamburg, Germany), Shern Tee (UQ, Brisbane, Australia) and
+Robert Meissner (TUHH, Hamburg, Germany).
+
+**Install:**
+
+This package has :ref:`specific installation instructions <electrode>` on the
+:doc:`Build extras <Build_extras>` page.
+
+**Supporting info:**
+
+* :doc:`fix electrode/conp <fix_electrode_conp>`
+
+----------
+
 .. _PKG-EXTRA-COMPUTE:
 
 EXTRA-COMPUTE package
@@ -1396,17 +1461,25 @@ MDI package
 
 **Contents:**
 
-A LAMMPS command and fix to allow client-server coupling of LAMMPS to
-other atomic or molecular simulation codes via the `MolSSI Driver Interface
+A LAMMPS command and fixes to allow client-server coupling of LAMMPS
+to other atomic or molecular simulation codes or materials modeling
+workflows via the `MolSSI Driver Interface
 (MDI) library <https://molssi-mdi.github.io/MDI_Library/html/index.html>`_.
 
 **Author:** Taylor Barnes - MolSSI, taylor.a.barnes at gmail.com
 
+**Install:**
+
+This package has :ref:`specific installation instructions <mdi>` on
+the :doc:`Build extras <Build_extras>` page.
+
 **Supporting info:**
 
 * src/MDI/README
-* :doc:`mdi/engine <mdi_engine>`
-* :doc:`fix mdi/engine <fix_mdi_engine>`
+* lib/mdi/README
+* :doc:`Howto MDI <Howto_mdi>`
+* :doc:`mdi <mdi>`
+* :doc:`fix mdi/aimd <fix_mdi_aimd>`
 * examples/PACKAGES/mdi
 
 ----------
@@ -1483,32 +1556,6 @@ Philipp Kloza (U Cambridge)
 
 ----------
 
-.. _PKG-MESSAGE:
-
-MESSAGE package
----------------
-
-**Contents:**
-
-Commands to use LAMMPS as either a client or server and couple it to
-another application.
-
-**Install:**
-
-This package has :ref:`specific installation instructions <message>` on the :doc:`Build extras <Build_extras>` page.
-
-**Supporting info:**
-
-* src/MESSAGE: filenames -> commands
-* lib/message/README
-* :doc:`message <message>`
-* :doc:`fix client/md <fix_client_md>`
-* :doc:`server md <server_md>`
-* :doc:`server mc <server_mc>`
-* examples/message
-
-----------
-
 .. _PKG-MGPT:
 
 MGPT package
@@ -1564,7 +1611,6 @@ listing, "ls src/MISC", to see the list of commands.
 * :doc:`pair_style list <pair_list>`
 * :doc:`pair_style srp <pair_srp>`
 * :doc:`pair_style tracker <pair_tracker>`
-* :doc:`fix pair/tracker <fix_pair_tracker>`
 
 ----------
 
@@ -2558,18 +2604,20 @@ SMTBQ package
 
 **Contents:**
 
-A pair style which implements a Second Moment Tight Binding model with
-QEq charge equilibration (SMTBQ) potential for the description of
-ionocovalent bonds in oxides.
+Pair styles which implement Second Moment Tight Binding models.
+One with QEq charge equilibration (SMTBQ) for the description of
+ionocovalent bonds in oxides, and two more as plain SMATB models.
 
-**Authors:** Nicolas Salles, Emile Maras, Olivier Politano, and Robert
-Tetot (LAAS-CNRS, France).
+**Authors:** SMTBQ: Nicolas Salles, Emile Maras, Olivier Politano, and Robert
+Tetot (LAAS-CNRS, France);
+SMATB: Daniele Rapetti (Politecnico di Torino)
 
 **Supporting info:**
 
 * src/SMTBQ: filenames -> commands
 * src/SMTBQ/README
 * :doc:`pair_style smtbq <pair_smtbq>`
+* :doc:`pair_style smatb <pair_smatb>`, :doc:`pair_style smatb/single <pair_smatb>`
 * examples/PACKAGES/smtbq
 
 ----------

doc/src/Packages_list.rst

@@ -58,6 +58,11 @@ whether an extra library is needed to build and use the package:
      - :doc:`Howto body <Howto_body>`
      - body
      - no
+   * - :ref:`BPM <PKG-BPM>`
+     - bonded particle models
+     - :doc:`Howto bpm <Howto_bpm>`
+     - bpm
+     - no
    * - :ref:`BROWNIAN <PKG-BROWNIAN>`
      - Brownian dynamics, self-propelled particles
      - :doc:`fix brownian <fix_brownian>`, :doc:`fix propel/self <fix_propel_self>`
@@ -143,6 +148,11 @@ whether an extra library is needed to build and use the package:
      - :doc:`pair_style eff/cut <pair_eff>`
      - PACKAGES/eff
      - no
+   * - :ref:`ELECTRODE <PKG-ELECTRODE>`
+     - electrode charges to match potential
+     - :doc:`fix electrode/conp <fix_electrode_conp>`
+     - PACKAGES/electrode
+     - no
    * - :ref:`EXTRA-COMPUTE <PKG-EXTRA-COMPUTE>`
      - additional compute styles
      - :doc:`compute <compute>`
@@ -244,7 +254,7 @@ whether an extra library is needed to build and use the package:
      - n/a
      - no
    * - :ref:`MDI <PKG-MDI>`
-     - client-server coupling
+     - client-server code coupling
      - :doc:`MDI Howto <Howto_mdi>`
      - PACKAGES/mdi
      - ext
@@ -258,11 +268,6 @@ whether an extra library is needed to build and use the package:
      - pair styles :doc:`mesont/tpm <pair_mesont_tpm>`, :doc:`mesocnt <pair_mesocnt>`
      - PACKAGES/mesont
      - int
-   * - :ref:`MESSAGE <PKG-MESSAGE>`
-     - client/server messaging
-     - :doc:`message <message>`
-     - message
-     - int
    * - :ref:`MGPT <PKG-MGPT>`
      - fast MGPT multi-ion potentials
      - :doc:`pair_style mgpt <pair_mgpt>`
@@ -429,8 +434,8 @@ whether an extra library is needed to build and use the package:
      - n/a
      - no
    * - :ref:`SMTBQ <PKG-SMTBQ>`
-     - second moment tight binding potential
-     - :doc:`pair_style smtbq <pair_smtbq>`
+     - second moment tight binding potentials
+     - :doc:`pair_style smtbq <pair_smtbq>` :doc:`pair_style smatb <pair_smatb>`
      - PACKAGES/smtbq
      - no
    * - :ref:`SPH <PKG-SPH>`

doc/src/Run_options.rst

@@ -226,15 +226,6 @@ other executable(s) perform an MPI_Comm_split() with their own colors
 to shrink the MPI_COMM_WORLD communication to be the subset of
 processors they are actually running on.
 
-Currently, this is only used in LAMMPS to perform client/server
-messaging with another application.  LAMMPS can act as either a client
-or server (or both).  More details are given on the :doc:`Howto client/server <Howto_client_server>` doc page.
-
-Specifically, this refers to the "mpi/one" mode of messaging provided
-by the :doc:`message <message>` command and the CSlib library LAMMPS
-links with from the lib/message directory.  See the
-:doc:`message <message>` command for more details.
-
 ----------
 
 .. _cite:

doc/src/Tools.rst

@@ -56,6 +56,7 @@ Pre-processing tools
    * :ref:`moltemplate <moltemplate>`
    * :ref:`msi2lmp <msi>`
    * :ref:`polybond <polybond>`
+   * :ref:`stl_bin2txt <stlconvert>`
 
 
 Post-processing tools
@@ -87,7 +88,7 @@ Miscellaneous tools
 .. table_from_list::
    :columns: 6
 
-   * :ref:`CMake <cmake>`
+   * :ref:`LAMMPS coding standards <coding_standard>`
    * :ref:`emacs <emacs>`
    * :ref:`i-pi <ipi>`
    * :ref:`kate <kate>`
@@ -189,27 +190,32 @@ for the :doc:`chain benchmark <Speed_bench>`.
 
 ----------
 
-.. _cmake:
+.. _coding_standard:
 
-CMake tools
------------
+LAMMPS coding standard
+----------------------
 
-The ``cmbuild`` script is a wrapper around using ``cmake --build <dir>
---target`` and allows compiling LAMMPS in a :ref:`CMake build folder
-<cmake_build>` with a make-like syntax regardless of the actual build
-tool and the specific name of the program used (e.g. ``ninja-v1.10`` or
-``gmake``) when using ``-D CMAKE_MAKE_PROGRAM=<name>``.
+The ``coding_standard`` folder contains multiple python scripts to
+check for and apply some LAMMPS coding conventions.  The following
+scripts are available:
 
 .. parsed-literal::
 
-  Usage: cmbuild [-v] [-h] [-C <dir>] [-j <num>] [<target>]
+   permissions.py   # detects if sources have executable permissions and scripts have not
+   whitespace.py    # detects TAB characters and trailing whitespace
+   homepage.py      # detects outdated LAMMPS homepage URLs (pointing to sandia.gov instead of lammps.org)
+   errordocs.py     # detects deprecated error docs in header files
 
-  Options:
-    -h                print this message
-    -j <NUM>          allow processing of NUM concurrent tasks
-    -C DIRECTORY      execute build in folder DIRECTORY
-    -v                produce verbose output
+The tools need to be given the main folder of the LAMMPS distribution
+or individual file names as argument and will by default check them
+and report any non-compliance.  With the optional ``-f`` argument the
+corresponding script will try to change the non-compliant file(s) to
+match the conventions.
 
+For convenience this scripts can also be invoked by the make file in
+the ``src`` folder with, `make check-whitespace` or `make fix-whitespace`
+to either detect or edit the files.  Correspondingly for the other python
+scripts. `make check` will run all checks.
 
 ----------
 
@@ -1012,6 +1018,28 @@ For more details please see the README.md file in that folder.
 
 ----------
 
+.. _stlconvert:
+
+stl_bin2txt tool
+----------------
+
+The file stl_bin2txt.cpp converts binary STL files - like they are
+frequently offered for download on the web - into ASCII format STL files
+that LAMMPS can read with the :doc:`create_atoms mesh <create_atoms>` or
+the :doc:`fix smd/wall_surface <fix_smd_wall_surface>` commands.  The syntax
+for running the tool is
+
+.. code-block:: bash
+
+   stl_bin2txt infile.stl outfile.stl
+
+which creates outfile.stl from infile.stl.  This tool must be compiled
+on a platform compatible with the byte-ordering that was used to create
+the binary file.  This usually is a so-called little endian hardware
+(like x86).
+
+----------
+
 .. _swig:
 
 SWIG interface

doc/src/angle_style.rst

@@ -87,7 +87,7 @@ of (g,i,k,o,t) to indicate which accelerated styles exist.
 * :doc:`dipole <angle_dipole>` - angle that controls orientation of a point dipole
 * :doc:`fourier <angle_fourier>` - angle with multiple cosine terms
 * :doc:`fourier/simple <angle_fourier_simple>` - angle with a single cosine term
-* :doc:`gaussian <angle_gaussian>` - multicentered Gaussian-based angle potential
+* :doc:`gaussian <angle_gaussian>` - multi-centered Gaussian-based angle potential
 * :doc:`harmonic <angle_harmonic>` - harmonic angle
 * :doc:`mm3 <angle_mm3>` - anharmonic angle
 * :doc:`quartic <angle_quartic>` - angle with cubic and quartic terms

doc/src/atom_style.rst

@@ -10,7 +10,7 @@ Syntax
 
    atom_style style args
 
-* style = *angle* or *atomic* or *body* or *bond* or *charge* or *dipole* or  *dpd* or *edpd* or *electron* or *ellipsoid* or *full* or *line* or *mdpd* or *molecular* or *oxdna* or *peri* or *smd* or *sph* or *sphere* or *spin* or *tdpd* or *tri* or *template* or *hybrid*
+* style = *angle* or *atomic* or *body* or *bond* or *charge* or *dipole* or  *dpd* or *edpd* or *electron* or *ellipsoid* or *full* or *line* or *mdpd* or *molecular* or *oxdna* or *peri* or *smd* or *sph* or *sphere* or *bpm/sphere* or *spin* or *tdpd* or *tri* or *template* or *hybrid*
 
   .. parsed-literal::
 
@@ -21,6 +21,7 @@ Syntax
                          see the :doc:`Howto body <Howto_body>` doc
                          page for details
          *sphere* arg = 0/1 (optional) for static/dynamic particle radii
+         *bpm/sphere* arg = 0/1 (optional) for static/dynamic particle radii
          *tdpd* arg = Nspecies
            Nspecies = # of chemical species
          *template* arg = template-ID
@@ -120,6 +121,8 @@ quantities.
 +--------------+-----------------------------------------------------+--------------------------------------+
 | *sphere*     | diameter, mass, angular velocity                    | granular models                      |
 +--------------+-----------------------------------------------------+--------------------------------------+
+| *bpm/sphere* | diameter, mass, angular velocity, quaternion        | granular bonded particle models (BPM)|
++--------------+-----------------------------------------------------+--------------------------------------+
 | *spin*       | magnetic moment                                     | system with magnetic particles       |
 +--------------+-----------------------------------------------------+--------------------------------------+
 | *tdpd*       | chemical concentration                              | tDPD particles                       |
@@ -141,8 +144,9 @@ quantities.
    output the custom values.
 
 All of the above styles define point particles, except the *sphere*,
-*ellipsoid*, *electron*, *peri*, *wavepacket*, *line*, *tri*, and
-*body* styles, which define finite-size particles.  See the :doc:`Howto spherical <Howto_spherical>` page for an overview of using
+*bpm/sphere*, *ellipsoid*, *electron*, *peri*, *wavepacket*, *line*,
+*tri*, and *body* styles, which define finite-size particles.  See the
+:doc:`Howto spherical <Howto_spherical>` page for an overview of using
 finite-size particle models with LAMMPS.
 
 All of the point-particle styles assign mass to particles on a
@@ -150,15 +154,15 @@ per-type basis, using the :doc:`mass <mass>` command, The finite-size
 particle styles assign mass to individual particles on a per-particle
 basis.
 
-For the *sphere* style, the particles are spheres and each stores a
+For the *sphere* and *bpm/sphere* styles, the particles are spheres 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 npt/sphere
-<fix_npt_sphere>` commands, spheres can be effectively treated as 2d
+<fix_npt_sphere>` commands for the *sphere* style, spheres can be effectively treated as 2d
 discs for a 2d simulation if desired.  See also the :doc:`set
-density/disc <set>` command.  The *sphere* style takes an optional 0
+density/disc <set>` command.  The *sphere* and *bpm/sphere* styles take an optional 0
 or 1 argument.  A value of 0 means the radius of each sphere is
 constant for the duration of the simulation.  A value of 1 means the
 radii may vary dynamically during the simulation, e.g. due to use of
@@ -195,6 +199,8 @@ position, which is represented by the eradius = electron size.
 For the *peri* style, the particles are spherical and each stores a
 per-particle mass and volume.
 
+The *bpm/sphere* style is part of the BPM package.
+
 The *oxdna* style is for coarse-grained nucleotides and stores the
 3'-to-5' polarity of the nucleotide strand, which is set through
 the bond topology in the data file. The first (second) atom in a

doc/src/bond_bpm_rotational.rst

@@ -0,0 +1,258 @@
+.. index:: bond_style bpm/rotational
+
+bond_style bpm/rotational command
+=================================
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   bond_style bpm/rotational keyword value attribute1 attribute2 ...
+
+* optional keyword = *overlay/pair* or *store/local* or *smooth*
+
+  .. parsed-literal::
+
+       *store/local* values = fix_ID N attributes ...
+          * fix_ID = ID of associated internal fix to store data
+          * N = prepare data for output every this many timesteps
+          * attributes = zero or more of the below attributes may be appended
+
+            *id1, id2* = IDs of 2 atoms in the bond
+            *time* = the timestep the bond broke
+            *x, y, z* = the center of mass position of the 2 atoms when the bond broke (distance units)
+            *x/ref, y/ref, z/ref* = the initial center of mass position of the 2 atoms (distance units)
+
+       *overlay/pair* value = none
+          bonded particles will still interact with pair forces
+
+       *smooth* value = *yes* or *no*
+          smooths bond forces near the breaking point
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   bond_style bpm/rotational
+   bond_coeff 1 1.0 0.2 0.02 0.02 0.20 0.04 0.04 0.04 0.1 0.02 0.002 0.002
+
+   bond_style bpm/rotational myfix 1000 time id1 id2
+   fix myfix all store/local 1000 3
+   dump 1 all local 1000 dump.broken f_myfix[1] f_myfix[2] f_myfix[3]
+   dump_modify 1 write_header no
+
+Description
+"""""""""""
+
+The *bpm/rotational* bond style computes forces and torques 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
+written to :doc:`binary restart files <restart>` such that restarting
+the system will not reset the reference state of a bond.
+
+Forces include a normal and tangential component. The base normal force
+has a magnitude of
+
+.. math::
+
+   f_r = k_r (r - r_0)
+
+where :math:`k_r` is a stiffness and :math:`r` is the current distance and
+:math:`r_0` is the initial distance between the two particles.
+
+A tangential force is applied perpendicular to the normal direction
+which is proportional to the tangential shear displacement with a
+stiffness of :math:`k_s`. This tangential force also induces a torque.
+In addition, bending and twisting torques are also applied to
+particles which are proportional to angular bending and twisting
+displacements with stiffnesses of :math`k_b` and :math:`k_t',
+respectively.  Details on the calculations of shear displacements and
+angular displacements can be found in :ref:`(Wang) <Wang2009>` and
+:ref:`(Wang and Mora) <Wang2009b>`.
+
+Bonds will break under sufficient stress. A breaking criteria is calculated
+
+.. math::
+
+   B = \mathrm{max}\{0, \frac{f_r}{f_{r,c}} + \frac{|f_s|}{f_{s,c}} +
+       \frac{|\tau_b|}{\tau_{b,c}} + \frac{|\tau_t|}{\tau_{t,c}} \}
+
+where :math:`|f_s|` is the magnitude of the shear force and
+:math:`|\tau_b|` and :math:`|\tau_t|` are the magnitudes of the
+bending and twisting forces, respectively. The corresponding variables
+:math:`f_{r,c}` :math:`f_{s,c}`, :math:`\tau_{b,c}`, and
+:math:`\tau_{t,c}` are critical limits to each force or torque.  If
+:math:`B` is ever equal to or exceeds one, the bond will break.  This
+is done by setting by setting its type to 0 such that forces and
+torques are no longer computed.
+
+After computing the base magnitudes of the forces and torques, they
+can be optionally multiplied by an extra factor :math:`w` to smoothly
+interpolate forces and torques to zero as the bond breaks. This term
+is calculated as :math:`w = (1.0 - B^4)`. This smoothing factor can be
+added or removed using the *smooth* keyword.
+
+Finally, additional damping forces and torques are applied to the two
+particles. A force is applied proportional to the difference in the
+normal velocity of particles using a similar construction as
+dissipative particle dynamics (:ref:`(Groot) <Groot3>`):
+
+.. math::
+
+   F_D = - \gamma_n w (\hat{r} \bullet \vec{v})
+
+where :math:`\gamma_n` is the damping strength, :math:`\hat{r}` is the
+radial normal vector, and :math:`\vec{v}` is the velocity difference
+between the two particles. Similarly, tangential forces are applied to
+each atom proportional to the relative differences in sliding
+velocities with a constant prefactor :math:`\gamma_s` (:ref:`(Wang et
+al.) <Wang20152>`) along with their associated torques. The rolling and
+twisting components of the relative angular velocities of the two
+atoms are also damped by applying torques with prefactors of
+:math:`\gamma_r` and :math:`\gamma_t`, respectively.
+
+The following coefficients 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_r`           (force/distance units)
+* :math:`k_s`           (force/distance units)
+* :math:`k_t`           (force*distance/radians units)
+* :math:`k_b`           (force*distance/radians units)
+* :math:`f_{r,c}`       (force units)
+* :math:`f_{s,c}`       (force units)
+* :math:`\tau_{b,c}`    (force*distance units)
+* :math:`\tau_{t,c}`    (force*distance units)
+* :math:`\gamma_n`      (force/velocity units)
+* :math:`\gamma_s`      (force/velocity units)
+* :math:`\gamma_r`      (force*distance/velocity units)
+* :math:`\gamma_t`      (force*distance/velocity units)
+
+By default, pair forces are not calculated between bonded particles.
+Pair forces can alternatively be overlaid on top of bond forces using
+the *overlay/pair* keyword. These settings require specific
+:doc:`special_bonds <special_bonds>` settings described in the
+restrictions.  Further details can be found in the `:doc: how to
+<Howto_BPM>` page on BPMs.
+
+If the *store/local* keyword is used, this 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
+local data to be accessed by other LAMMPS commands.
+Following any optional keyword/value arguments, a list of one or more
+attributes is specified.  These include the IDs of the two atoms in
+the bond. The other attributes for the two atoms include the timestep
+during which the bond broke and the current/initial center of mass
+position of the two atoms.
+
+Data is continuously accumulated over intervals of *N*
+timesteps. At the end of each interval, all of the saved accumulated
+data is deleted to make room for new data. Individual datum may
+therefore persist anywhere between *1* to *N* timesteps depending on
+when they are saved. This data can be accessed using the *fix_ID* and a
+:doc:`dump local <dump>` command. To ensure all data is output,
+the dump frequency should correspond to the same interval of *N*
+timesteps. A dump frequency of an integer multiple of *N* can be used
+to regularly output a sample of the accumulated data.
+
+Note that when unbroken bonds are dumped to a file via the
+:doc:`dump local <dump>` command, bonds with type 0 (broken bonds)
+are not included.
+The :doc:`delete_bonds <delete_bonds>` command can also be used to
+query the status of broken bonds or permanently delete them, e.g.:
+
+.. code-block:: LAMMPS
+
+   delete_bonds all stats
+   delete_bonds all bond 0 remove
+
+
+----------
+
+Restart and other info
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+This bond style writes the reference state of each bond to
+:doc:`binary restart files <restart>`. Loading a restart file will
+properly resume bonds.
+
+The single() function of these pair styles returns 0.0 for the energy
+of a pairwise interaction, since energy is not conserved in these
+dissipative potentials.  It also returns only the normal component of
+the pairwise interaction force.
+
+The accumulated data is not written to restart files and should be
+output before a restart file is written to avoid missing data.
+
+The internal fix calculates a local vector or local array depending on the
+number of input values.  The length of the vector or number of rows in
+the array is the number of recorded, lost interactions.  If a single
+input is specified, a local vector is produced.  If two or more inputs
+are specified, a local array is produced where the number of columns =
+the number of inputs.  The vector or array can be accessed by any
+command that uses local values from a compute as input.  See the
+:doc:`Howto output <Howto_output>` page for an overview of LAMMPS
+output options.
+
+The vector or array will be floating point values that correspond to
+the specified attribute.
+
+Restrictions
+""""""""""""
+
+This bond style can only be used if LAMMPS was built with the BPM
+package. See the :doc:`Build package <Build_package>` doc page for
+more info.
+
+By default if pair interactions are to be disabled, this bond style
+requires setting
+
+.. code-block:: LAMMPS
+
+   special_bonds lj 0 1 1 coul 1 1 1
+
+and :doc:`newton <newton>` must be set to bond off.  If the
+*overlay/pair* option is used, this bond style alternatively requires
+setting
+
+.. code-block:: LAMMPS
+
+   special_bonds lj/coul 1 1 1
+
+The *bpm/rotational* style requires :doc:`atom style bpm/sphere <atom_style>`.
+
+Related commands
+""""""""""""""""
+
+:doc:`bond_coeff <bond_coeff>`, :doc:`fix nve/bpm/sphere <fix_nve_bpm_sphere>`
+
+Default
+"""""""
+
+The option defaults are *smooth* = *yes*
+
+----------
+
+.. _Wang2009:
+
+**(Wang)** Wang, Acta Geotechnica, 4,
+p 117-127 (2009).
+
+.. _Wang2009b:
+
+**(Wang and Mora)** Wang, Mora, Advances in Geocomputing,
+119, p 183-228 (2009).
+
+.. _Groot3:
+
+**(Groot)** Groot and Warren, J Chem Phys, 107, 4423-35 (1997).
+
+.. _Wang20152:
+
+**(Wang et al, 2015)** Wang, Y., Alonso-Marroquin, F., & Guo,
+W. W. (2015).  Rolling and sliding in 3-D discrete element
+models. Particuology, 23, 49-55.

doc/src/bond_bpm_spring.rst

@@ -0,0 +1,202 @@
+.. index:: bond_style bpm/spring
+
+bond_style bpm/spring command
+=============================
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   bond_style bpm/spring keyword value attribute1 attribute2 ...
+
+* optional keyword = *overlay/pair* or *store/local* or *smooth*
+
+  .. parsed-literal::
+
+       *store/local* values = fix_ID N attributes ...
+          * fix_ID = ID of associated internal fix to store data
+          * N = prepare data for output every this many timesteps
+          * attributes = zero or more of the below attributes may be appended
+
+            *id1, id2* = IDs of 2 atoms in the bond
+            *time* = the timestep the bond broke
+            *x, y, z* = the center of mass position of the 2 atoms when the bond broke (distance units)
+            *x/ref, y/ref, z/ref* = the initial center of mass position of the 2 atoms (distance units)
+
+       *overlay/pair* value = none
+          bonded particles will still interact with pair forces
+
+       *smooth* value = *yes* or *no*
+          smooths bond forces near the breaking point
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   bond_style bpm/spring
+   bond_coeff 1 1.0 0.05 0.1
+
+   bond_style bpm/spring myfix 1000 time id1 id2
+   dump 1 all local 1000 dump.broken f_myfix[1] f_myfix[2] f_myfix[3]
+   dump_modify 1 write_header no
+
+Description
+"""""""""""
+
+The *bpm/spring* bond style computes forces and torques 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
+written to :doc:`binary restart files <restart>` such that restarting
+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
+
+.. math::
+
+   F = k (r - r_0) w
+
+where :math:`k_r` is a stiffness, :math:`r` is the current distance
+and :math:`r_0` is the initial distance between the two particles, and
+:math:`w` is an optional smoothing factor discussed below. Bonds will
+break at a strain of :math:`\epsilon_c`.  This is done by setting by
+setting its type to 0 such that forces are no longer computed.
+
+An additional damping force is applied to the bonded
+particles.  This forces is proportional to the difference in the
+normal velocity of particles using a similar construction as
+dissipative particle dynamics (:ref:`(Groot) <Groot4>`):
+
+.. math::
+
+   F_D = - \gamma w (\hat{r} \bullet \vec{v})
+
+where :math:`\gamma` is the damping strength, :math:`\hat{r}` is the
+radial normal vector, and :math:`\vec{v}` is the velocity difference
+between the two particles.
+
+The smoothing factor :math:`w` can be added or removed using the
+*smooth* keyword. It is constructed such that forces smoothly go
+to zero, avoiding discontinuities, as bonds approach the critical strain
+
+.. math::
+
+   w = 1.0 - \left( \frac{r - r_0}{r_0 \epsilon_c} \right)^8 .
+
+The following coefficients 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`             (force/distance units)
+* :math:`\epsilon_c`    (unit less)
+* :math:`\gamma`        (force/velocity units)
+
+By default, pair forces are not calculated between bonded particles.
+Pair forces can alternatively be overlaid on top of bond forces using
+the *overlay/pair* keyword. These settings require specific
+:doc:`special_bonds <special_bonds>` settings described in the
+restrictions.  Further details can be found in the `:doc: how to
+<Howto_BPM>` page on BPMs.
+
+If the *store/local* keyword is used, this 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
+local data to be accessed by other LAMMPS commands.
+Following any optional keyword/value arguments, a list of one or more
+attributes is specified.  These include the IDs of the two atoms in
+the bond. The other attributes for the two atoms include the timestep
+during which the bond broke and the current/initial center of mass
+position of the two atoms.
+
+Data is continuously accumulated over intervals of *N*
+timesteps. At the end of each interval, all of the saved accumulated
+data is deleted to make room for new data. Individual datum may
+therefore persist anywhere between *1* to *N* timesteps depending on
+when they are saved. This data can be accessed using the *fix_ID* and a
+:doc:`dump local <dump>` command. To ensure all data is output,
+the dump frequency should correspond to the same interval of *N*
+timesteps. A dump frequency of an integer multiple of *N* can be used
+to regularly output a sample of the accumulated data.
+
+Note that when unbroken bonds are dumped to a file via the
+:doc:`dump local <dump>` command, bonds with type 0 (broken bonds)
+are not included.
+The :doc:`delete_bonds <delete_bonds>` command can also be used to
+query the status of broken bonds or permanently delete them, e.g.:
+
+.. code-block:: LAMMPS
+
+   delete_bonds all stats
+   delete_bonds all bond 0 remove
+
+----------
+
+Restart and other info
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+This bond style writes the reference state of each bond to
+:doc:`binary restart files <restart>`. Loading a restart
+file will properly resume bonds.
+
+The single() function of these pair styles returns 0.0 for the energy
+of a pairwise interaction, since energy is not conserved in these
+dissipative potentials.
+
+The accumulated data is not written to restart files and should be
+output before a restart file is written to avoid missing data.
+
+The internal fix calculates a local vector or local array depending on the
+number of input values.  The length of the vector or number of rows in
+the array is the number of recorded, lost interactions.  If a single
+input is specified, a local vector is produced.  If two or more inputs
+are specified, a local array is produced where the number of columns =
+the number of inputs.  The vector or array can be accessed by any
+command that uses local values from a compute as input.  See the
+:doc:`Howto output <Howto_output>` page for an overview of LAMMPS
+output options.
+
+The vector or array will be floating point values that correspond to
+the specified attribute.
+
+Restrictions
+""""""""""""
+
+This bond style can only be used if LAMMPS was built with the BPM
+package. See the :doc:`Build package <Build_package>` doc page for
+more info.
+
+By default if pair interactions are to be disabled, this bond style
+requires setting
+
+.. code-block:: LAMMPS
+
+   special_bonds lj 0 1 1 coul 1 1 1
+
+and :doc:`newton <newton>` must be set to bond off.  If the
+*overlay/pair* option is used, this bond style alternatively requires
+setting
+
+.. code-block:: LAMMPS
+
+   special_bonds lj/coul 1 1 1
+
+Related commands
+""""""""""""""""
+
+:doc:`bond_coeff <bond_coeff>`, :doc:`pair bpm/spring <pair_bpm_spring>`
+
+Default
+"""""""
+
+The option defaults are *smooth* = *yes*
+
+----------
+
+.. _Groot4:
+
+**(Groot)** Groot and Warren, J Chem Phys, 107, 4423-35 (1997).

doc/src/bond_quartic.rst

@@ -67,7 +67,8 @@ local maximum.  If a bond length ever becomes :math:`> R_c`, LAMMPS "breaks"
 the bond, which means two things.  First, the bond potential is turned
 off by setting its type to 0, and is no longer computed.  Second, a
 pairwise interaction between the two atoms is turned on, since they
-are no longer bonded.
+are no longer bonded. See the :doc:`Howto <Howto_broken_bonds>` page
+on broken bonds for more information.
 
 LAMMPS does the second task via a computational sleight-of-hand.  It
 subtracts the pairwise interaction as part of the bond computation.

doc/src/bond_style.rst

@@ -84,6 +84,8 @@ accelerated styles exist.
 * :doc:`zero <bond_zero>` - topology but no interactions
 * :doc:`hybrid <bond_hybrid>` - define multiple styles of bond interactions
 
+* :doc:`bpm/rotational <bond_bpm_rotational>` - breakable bond with forces and torques based on deviation from reference state
+* :doc:`bpm/spring <bond_bpm_spring>` - breakable bond with forces based on deviation from reference length
 * :doc:`class2 <bond_class2>` - COMPASS (class 2) bond
 * :doc:`fene <bond_fene>` - FENE (finite-extensible non-linear elastic) bond
 * :doc:`fene/expand <bond_fene_expand>` - FENE bonds with variable size particles

doc/src/commands_list.rst

@@ -59,8 +59,7 @@ Commands
    lattice
    log
    mass
-   mdi_engine
-   message
+   mdi
    min_modify
    min_spin
    min_style
@@ -96,9 +95,6 @@ Commands
    restart
    run
    run_style
-   server
-   server_mc
-   server_md
    set
    shell
    special_bonds

doc/src/compute.rst

@@ -179,6 +179,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
 * :doc:`body/local <compute_body_local>` - attributes of body sub-particles
 * :doc:`bond <compute_bond>` - energy of each bond sub-style
 * :doc:`bond/local <compute_bond_local>` - distance and energy of each bond
+* :doc:`born/matrix <compute_born_matrix>` - second derivative or potential with respect to strain
 * :doc:`centro/atom <compute_centro_atom>` - centro-symmetry parameter for each atom
 * :doc:`centroid/stress/atom <compute_stress_atom>` - centroid based stress tensor for each atom
 * :doc:`chunk/atom <compute_chunk_atom>` - assign chunk IDs to each atom
@@ -208,7 +209,8 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
 * :doc:`erotate/sphere/atom <compute_erotate_sphere_atom>` - rotational energy for each spherical particle
 * :doc:`event/displace <compute_event_displace>` - detect event on atom displacement
 * :doc:`fabric <compute_fabric>` - calculates fabric tensors from pair interactions
-* :doc:`fep <compute_fep>` -
+* :doc:`fep <compute_fep>` - compute free energies for alchemical transformation from perturbation theory
+* :doc:`fep/ta <compute_fep_ta>` - compute free energies for a test area perturbation
 * :doc:`force/tally <compute_tally>` - force between two groups of atoms via the tally callback mechanism
 * :doc:`fragment/atom <compute_cluster_atom>` - fragment ID for each atom
 * :doc:`global/atom <compute_global_atom>` -
@@ -235,6 +237,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
 * :doc:`msd <compute_msd>` - mean-squared displacement of group of atoms
 * :doc:`msd/chunk <compute_msd_chunk>` - mean-squared displacement for each chunk
 * :doc:`msd/nongauss <compute_msd_nongauss>` - MSD and non-Gaussian parameter of group of atoms
+* :doc:`nbond/atom <compute_nbond_atom>` - calculates number of bonds per atom
 * :doc:`omega/chunk <compute_omega_chunk>` - angular velocity for each chunk
 * :doc:`orientorder/atom <compute_orientorder_atom>` - Steinhardt bond orientational order parameters Ql
 * :doc:`pair <compute_pair>` - values computed by a pair style

doc/src/compute_born_matrix.rst

@@ -0,0 +1,213 @@
+.. index:: compute born/matrix
+
+compute born/matrix command
+===========================
+
+Syntax
+""""""
+
+.. parsed-literal::
+
+   compute ID group-ID born/matrix keyword value ...
+
+* ID, group-ID are documented in :doc:`compute <compute>` command
+* born/matrix = style name of this compute command
+* zero or more keyword/value pairs may be appended
+
+  .. parsed-literal::
+
+     keyword = *numdiff*
+       *numdiff* values = delta virial-ID
+         delta = magnitude of strain (dimensionless)
+         virial-ID = ID of pressure compute for virial (string)
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   compute 1 all born/matrix
+   compute 1 all born/matrix bond angle
+   compute 1 all born/matrix numdiff 1.0e-4 myvirial
+
+Description
+"""""""""""
+
+Define a compute that calculates
+:math:`\frac{\partial{}^2U}{\partial\varepsilon_{i}\partial\varepsilon_{j}}` the
+second derivatives of the potential energy :math:`U` w.r.t. strain
+tensor :math:`\varepsilon` elements. These values are related to:
+
+.. math::
+
+    C^{B}_{i,j}=\frac{1}{V}\frac{\partial{}^2U}{\partial{}\varepsilon_{i}\partial\varepsilon_{j}}
+
+also called the Born term of elastic constants in the stress-stress fluctuation
+formalism. This quantity can be used to compute the elastic constant tensor.
+Using the symmetric Voigt notation, the elastic constant tensor can be written
+as a 6x6 symmetric matrix:
+
+.. math::
+
+    C_{i,j} = \langle{}C^{B}_{i,j}\rangle
+             + \frac{V}{k_{B}T}\left(\langle\sigma_{i}\sigma_{j}\rangle\right.
+             \left.- \langle\sigma_{i}\rangle\langle\sigma_{j}\rangle\right)
+             + \frac{Nk_{B}T}{V}
+               \left(\delta_{i,j}+(\delta_{1,i}+\delta_{2,i}+\delta_{3,i})\right.
+               \left.*(\delta_{1,j}+\delta_{2,j}+\delta_{3,j})\right)
+
+In the above expression, :math:`\sigma` stands for the virial stress
+tensor, :math:`\delta` is the Kronecker delta and the usual notation apply for
+the number of particle, the temperature and volume respectively :math:`N`,
+:math:`T` and :math:`V`. :math:`k_{B}` is the Boltzmann constant.
+
+The Born term is a symmetric 6x6 matrix, as is the matrix of second derivatives
+of potential energy w.r.t strain,
+whose 21 independent elements are output in this order:
+
+.. math::
+
+    \begin{matrix}
+       C_{1}  & C_{7}   & C_{8}  & C_{9}  & C_{10} & C_{11} \\
+       C_{7}  & C_{2}   & C_{12} & C_{13} & C_{14} & C_{15} \\
+       \vdots & C_{12}  & C_{3}  & C_{16} & C_{17} & C_{18} \\
+       \vdots & C_{13}  & C_{16} & C_{4}  & C_{19} & C_{20} \\
+       \vdots & \vdots  & \vdots & C_{19} & C_{5}  & C_{21} \\
+       \vdots & \vdots  & \vdots & \vdots & C_{21} & C_{6}
+    \end{matrix}
+
+in this matrix the indices of :math:`C_{k}` value are the corresponding element
+:math:`k` in the global vector output by this compute. Each term comes from the sum
+of the derivatives of every contribution to the potential energy
+in the system as explained in :ref:`(VanWorkum)
+<VanWorkum>`.
+
+The output can be accessed using usual Lammps routines:
+
+.. code-block:: LAMMPS
+
+   compute 1 all born/matrix
+   compute 2 all pressure NULL virial
+   variable S1 equal -c_2[1]
+   variable S2 equal -c_2[2]
+   variable S3 equal -c_2[3]
+   variable S4 equal -c_2[4]
+   variable S5 equal -c_2[5]
+   variable S6 equal -c_2[6]
+   fix 1 all ave/time 1 1 1 v_S1 v_S2 v_S3 v_S4 v_S5 v_S6 c_1[*] file born.out
+
+In this example, the file *born.out* will contain the information needed to
+compute the first and second terms of the elastic constant matrix in a post
+processing procedure. The other required quantities can be accessed using any
+other *LAMMPS* usual method. Several examples of this method are
+provided in the examples/ELASTIC_T/BORN_MATRIX directory
+described on the :doc:`Examples <Examples>` doc page.
+
+NOTE: In the above :math:`C_{i,j}` computation, the fluctuation
+term involving the virial stress tensor :math:`\sigma` is the
+covariance between each elements. In a
+solid the stress fluctuations can vary rapidly, while average
+fluctuations can be slow to converge.
+A detailed analysis of the convergence rate of all the terms in
+the elastic tensor
+is provided in the paper by Clavier et al. :ref:`(Clavier) <Clavier2>`.
+
+Two different computation methods for the Born matrix are implemented in this
+compute and are mutually exclusive.
+
+The first one is a direct computation from the analytical formula from the
+different terms of the potential used for the simulations :ref:`(VanWorkum)
+<VanWorkum>`. However, the implementation of such derivations must be done
+for every potential form. This has not been done yet and can be very
+complicated for complex potentials. At the moment a warning message is
+displayed for every term that is not supporting the compute at the moment.
+This method is the default for now.
+
+The second method uses finite differences of energy to numerically approximate
+the second derivatives :ref:`(Zhen) <Zhen>`. This is useful when using
+interaction styles for which the analytical second derivatives have not been
+implemented. In this cases, the compute applies linear strain fields of
+magnitude *delta* to all the atoms relative to a point at the center of the
+box. The strain fields are in six different directions, corresponding to the
+six Cartesian components of the stress tensor defined by LAMMPS. For each
+direction it applies the strain field in both the positive and negative senses,
+and the new stress virial tensor of the entire system is calculated after each.
+The difference in these two virials divided by two times *delta*, approximates
+the corresponding components of the second derivative, after applying a
+suitable unit conversion.
+
+.. note::
+
+   It is important to choose a suitable value for delta, the magnitude of
+   strains that are used to generate finite difference
+   approximations to the exact virial stress.  For typical systems, a value in
+   the range of 1 part in 1e5 to 1e6 will be sufficient.
+   However, the best value will depend on a multitude of factors
+   including the stiffness of the interatomic potential, the thermodynamic
+   state of the material being probed, and so on. The only way to be sure
+   that you have made a good choice is to do a sensitivity study on a
+   representative atomic configuration, sweeping over a wide range of
+   values of delta. If delta is too small, the output values will vary
+   erratically due to truncation effects. If delta is increased beyond a
+   certain point, the output values will start to vary smoothly with
+   delta, due to growing contributions from higher order derivatives. In
+   between these two limits, the numerical virial values should be largely
+   independent of delta.
+
+The keyword requires the additional arguments *delta* and *virial-ID*.
+*delta* gives the size of the applied strains. *virial-ID* gives
+the ID string of the pressure compute that provides the virial stress tensor,
+requiring that it use the virial keyword e.g.
+
+.. code-block:: LAMMPS
+
+   compute myvirial all pressure NULL virial
+   compute 1 all born/matrix numdiff 1.0e-4 myvirial
+
+**Output info:**
+
+This compute calculates a global vector with 21 values that are
+the second derivatives of the potential energy w.r.t. strain.
+The values are in energy units.
+The values are ordered as explained above. These values can be used
+by any command that uses global values from a compute as input. See
+the :doc:`Howto output <Howto_output>` doc page for an overview of
+LAMMPS output options.
+
+The array values calculated by this compute are all "extensive".
+
+Restrictions
+""""""""""""
+
+This compute 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.  LAMMPS was built with that package.  See
+the :doc:`Build package <Build_package>` page for more info.
+
+The Born term can be decomposed as a product of two terms. The first one is a
+general term which depends on the configuration. The second one is specific to
+every interaction composing your force field (non-bonded, bonds, angle...).
+Currently not all LAMMPS interaction styles implement the *born_matrix* method
+giving first and second order derivatives and LAMMPS will exit with an error if
+this compute is used with such interactions unless the *numdiff* option is
+also used. The *numdiff* option cannot be used with any other keyword. In this
+situation, LAMMPS will also exit with an error.
+
+Default
+"""""""
+
+none
+
+----------
+
+.. _VanWorkum:
+
+**(Van Workum)** K. Van Workum et al., J. Chem. Phys. 125 144506 (2006)
+
+.. _Clavier2:
+
+**(Clavier)** G. Clavier, N. Desbiens, E. Bourasseau, V. Lachet, N. Brusselle-Dupend and B. Rousseau, Mol Sim, 43, 1413 (2017).
+
+.. _Zhen:
+
+**(Zhen)** Y. Zhen, C. Chu, Computer Physics Communications 183(2012)261-265

doc/src/compute_fep_ta.rst

@@ -0,0 +1,99 @@
+.. index:: compute fep/ta
+
+compute fep/ta command
+======================
+
+Syntax
+""""""
+
+.. parsed-literal::
+
+   compute ID group-ID fep/ta temp plane scale_factor keyword value ...
+
+* ID, group-ID are documented in the :doc:`compute <compute>` command
+* fep/ta = name of this compute command
+* temp = external temperature (as specified for constant-temperature run)
+* plane = *xy* or *xz* or *yz*
+* scale_factor = multiplicative factor for change in plane area
+* zero or more keyword/value pairs may be appended
+* keyword = *tail*
+
+  .. parsed-literal::
+
+       *tail* value = *no* or *yes*
+         *no* = ignore tail correction to pair energies (usually small in fep)
+         *yes* = include tail correction to pair energies
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   compute 1 all fep/ta 298 xy 1.0005
+
+Description
+"""""""""""
+
+Define a computation that calculates the change in the free energy due
+to a test-area (TA) perturbation :ref:`(Gloor) <Gloor>`. The test-area
+approach can be used to determine the interfacial tension of the system
+in a single simulation:
+
+.. math::
+
+   \gamma = \lim_{\Delta \mathcal{A} \to 0} \left( \frac{\Delta A_{0 \to 1 }}{\Delta \mathcal{A}}\right)_{N,V,T}
+   = - \frac{kT}{\Delta \mathcal{A}} \ln \left< \exp(-(U_1 - U_0)/kT) \right>_0
+
+During the perturbation, both axes of *plane* are scaled by multiplying
+:math:`\sqrt{scale\_factor}`, while the other axis divided by
+*scale_factor* such that the overall volume of the system is maintained.
+
+The *tail* keyword controls the calculation of the tail correction to
+"van der Waals" pair energies beyond the cutoff, if this has been
+activated via the :doc:`pair_modify <pair_modify>` command. If the
+perturbation is small, the tail contribution to the energy difference
+between the reference and perturbed systems should be negligible.
+
+----------
+
+Output info
+"""""""""""
+
+This compute calculates a global vector of length 3 which contains the
+energy difference ( :math:`U_1-U_0` ) as c_ID[1], the Boltzmann factor
+:math:`\exp(-(U_1-U_0)/kT)`, as c_ID[2] and the change in the *plane*
+area :math:`\Delta \mathcal{A}` as c_ID[3]. :math:`U_1` is the potential
+energy of the perturbed state and :math:`U_0` is the potential energy of
+the reference state.  The energies include kspace terms if these are
+used in the simulation.
+
+These output results can be used by any command that uses a global
+scalar or vector from a compute as input.  See the :doc:`Howto output
+<Howto_output>` page for an overview of LAMMPS output options. For
+example, the computed values can be averaged using :doc:`fix ave/time
+<fix_ave_time>`.
+
+Restrictions
+""""""""""""
+
+Constraints, like fix shake, may lead to incorrect values for energy difference.
+
+This compute is distributed as the FEP 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:`compute fep <compute_fep>`
+
+Default
+"""""""
+
+The option defaults are *tail* = *no*\ .
+
+----------
+
+.. _Gloor:
+
+**(Gloor)** Gloor, J Chem Phys, 123, 134703 (2005)

doc/src/compute_msd.rst

@@ -75,10 +75,11 @@ solids undergoing thermal motion.
 .. note::
 
    Initial coordinates are stored in "unwrapped" form, by using the
-   image flags associated with each atom.  See the :doc:`dump custom <dump>` command for a discussion of "unwrapped" coordinates.
-   See the Atoms section of the :doc:`read_data <read_data>` command for a
-   discussion of image flags and how they are set for each atom.  You can
-   reset the image flags (e.g. to 0) before invoking this compute by
+   image flags associated with each atom.  See the :doc:`dump custom
+   <dump>` command for a discussion of "unwrapped" coordinates.  See the
+   Atoms section of the :doc:`read_data <read_data>` command for a
+   discussion of image flags and how they are set for each atom.  You
+   can reset the image flags (e.g. to 0) before invoking this compute by
    using the :doc:`set image <set>` command.
 
 .. note::
@@ -108,7 +109,8 @@ distance\^2 :doc:`units <units>`.
 
 Restrictions
 """"""""""""
- none
+
+Compute *msd* cannot be used with a dynamic group.
 
 Related commands
 """"""""""""""""

doc/src/compute_msd_nongauss.rst

@@ -74,8 +74,11 @@ the third is dimensionless.
 Restrictions
 """"""""""""
 
-This compute 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 *msd/nongauss* cannot be used with a dynamic group.
+
+This compute 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.
 
 Related commands
 """"""""""""""""

doc/src/compute_nbond_atom.rst

@@ -0,0 +1,52 @@
+.. index:: compute nbond/atom
+
+compute nbond/atom command
+==========================
+
+Syntax
+""""""
+
+.. parsed-literal::
+
+   compute ID group-ID nbond/atom
+
+* ID, group-ID are documented in :doc:`compute <compute>` command
+* nbond/atom = style name of this compute command
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   compute 1 all nbond/atom
+
+Description
+"""""""""""
+
+Define a computation that computes the number of bonds each atom is
+part of.  Bonds which are broken are not counted in the tally.  See
+the :doc:`Howto broken bonds <Howto_bpm>` page for more information.
+The number of bonds will be zero for atoms not in the specified
+compute group. This compute does not depend on Newton bond settings.
+
+Output info
+"""""""""""
+
+This compute calculates a per-atom vector, which can be accessed by
+any command that uses per-atom values from a compute as input.  See
+the :doc:`Howto output <Howto_output>` doc page for an overview of
+LAMMPS output options.
+
+Restrictions
+""""""""""""
+
+This fix can only be used if LAMMPS was built with the BPM package.
+See the :doc:`Build package <Build_package>` doc page for more info.
+
+Related commands
+""""""""""""""""
+
+Default
+"""""""
+
+none

doc/src/compute_reduce.rst

@@ -23,7 +23,7 @@ Syntax
        *reduce/region* arg = region-ID
          region-ID = ID of region to use for choosing atoms
 
-* mode = *sum* or *min* or *max* or *ave* or *sumsq* or *avesq*
+* mode = *sum* or *min* or *max* or *ave* or *sumsq* or *avesq* or *sumabs* or *aveabs*
 * one or more inputs can be listed
 * input = x, y, z, vx, vy, vz, fx, fy, fz, c_ID, c_ID[N], f_ID, f_ID[N], v_name
 
@@ -77,7 +77,10 @@ option sums the square of the values in the vector into a global
 total.  The *avesq* setting does the same as *sumsq*, then divides the
 sum of squares by the number of values.  The last two options can be
 useful for calculating the variance of some quantity, e.g. variance =
-sumsq - ave\^2.
+sumsq - ave\^2.  The *sumabs* option sums the absolute values in the
+vector into a global total.  The *aveabs* setting does the same as
+*sumabs*, then divides the sum of absolute values by the number of
+values.
 
 Each listed input is operated on independently.  For per-atom inputs,
 the group specified with this command means only atoms within the
@@ -189,7 +192,7 @@ value.  If multiple inputs are specified, this compute produces a
 global vector of values, the length of which is equal to the number of
 inputs specified.
 
-As discussed below, for the *sum* and *sumsq* modes, the value(s)
+As discussed below, for the *sum*, *sumabs* and *sumsq* modes, the value(s)
 produced by this compute are all "extensive", meaning their value
 scales linearly with the number of atoms involved.  If normalized
 values are desired, this compute can be accessed by the :doc:`thermo_style custom <thermo_style>` command with :doc:`thermo_modify norm yes <thermo_modify>` set as an option.  Or it can be accessed by a
@@ -208,7 +211,7 @@ compute as input.  See the :doc:`Howto output <Howto_output>` doc page
 for an overview of LAMMPS output options.
 
 All the scalar or vector values calculated by this compute are
-"intensive", except when the *sum* or *sumsq* modes are used on
+"intensive", except when the *sum*, *sumabs* or *sumsq* modes are used on
 per-atom or local vectors, in which case the calculated values are
 "extensive".
 

doc/src/compute_rigid_local.rst

@@ -127,19 +127,16 @@ The *vx*, *vy*, *vz*, *fx*, *fy*, *fz* attributes are components of
 the COM velocity and force on the COM of the body.
 
 The *omegax*, *omegay*, and *omegaz* attributes are the angular
-velocity components of the body around its COM.
+velocity components of the body in the system frame around its COM.
 
 The *angmomx*, *angmomy*, and *angmomz* attributes are the angular
-momentum components of the body around its COM.
+momentum components of the body in the system frame around its COM.
 
 The *quatw*, *quati*, *quatj*, and *quatk* attributes are the
 components of the 4-vector quaternion representing the orientation of
 the rigid body.  See the :doc:`set <set>` command for an explanation of
 the quaternion vector.
 
-The *angmomx*, *angmomy*, and *angmomz* attributes are the angular
-momentum components of the body around its COM.
-
 The *tqx*, *tqy*, *tqz* attributes are components of the torque acting
 on the body around its COM.
 

doc/src/compute_sna_atom.rst

@@ -33,7 +33,7 @@ Syntax
 * R_1, R_2,... = list of cutoff radii, one for each type (distance units)
 * w_1, w_2,... = list of neighbor weights, one for each type
 * zero or more keyword/value pairs may be appended
-* keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag* or *bikflag* or *switchinnerflag*
+* keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag* or *bikflag* or *switchinnerflag* or *sinner* or *dinner*
 
   .. parsed-literal::
 
@@ -59,9 +59,13 @@ Syntax
        *bikflag* value = *0* or *1* (only implemented for compute snap)
           *0* = per-atom bispectrum descriptors are summed over atoms
           *1* = per-atom bispectrum descriptors are not summed over atoms
-       *switchinnerflag* values = *rinnerlist* *drinnerlist*
-          *rinnerlist* = *ntypes* values of rinner (distance units)
-          *drinnerlist* = *ntypes* values of drinner (distance units)
+       *switchinnerflag* value = *0* or *1*
+          *0* = do not use inner switching function
+          *1* = use inner switching function
+       *sinner* values = *sinnerlist*
+          *sinnerlist* = *ntypes* values of *Sinner* (distance units)
+       *dinner* values = *dinnerlist*
+          *dinnerlist* = *ntypes* values of *Dinner* (distance units)
 
 Examples
 """"""""
@@ -73,7 +77,7 @@ Examples
    compute vb all sna/atom 1.4 0.95 6 2.0 1.0
    compute snap all snap 1.4 0.95 6 2.0 1.0
    compute snap all snap 1.0 0.99363 6 3.81 3.83 1.0 0.93 chem 2 0 1
-   compute snap all snap 1.0 0.99363 6 3.81 3.83 1.0 0.93 switchinnerflag 1.1 1.3 0.5 0.6
+   compute snap all snap 1.0 0.99363 6 3.81 3.83 1.0 0.93 switchinnerflag 1 sinner 1.35 1.6 dinner 0.25 0.3
 
 Description
 """""""""""
@@ -312,25 +316,29 @@ the resulting bispectrum rows are :math:`B_{i,k}` instead of just
 :math:`B_k`. In this case, the entries in the final column for these rows
 are set to zero.
 
-The keyword *switchinnerflag* activates an additional radial switching
+The keyword *switchinnerflag* with value 1
+activates an additional radial switching
 function similar to :math:`f_c(r)` above, but acting to switch off
 smoothly contributions from neighbor atoms at short separation distances.
 This is useful when SNAP is used in combination with a simple
-repulsive potential. The keyword is followed by the *ntypes*
-values for :math:`r_{inner}` and the *ntypes*
-values for :math:`\Delta r_{inner}`. For a neighbor atom at
+repulsive potential. For a neighbor atom at
 distance :math:`r`, its contribution is scaled by a multiplicative
 factor :math:`f_{inner}(r)` defined as follows:
 
 .. math::
 
-               = & 0,  r \leq r_{inner} \\
-  f_{inner}(r) = & \frac{1}{2}(1 - \cos(\pi \frac{r-r_{inner}}{\Delta r_{inner}})), r_{inner} < r \leq r_{inner} + \Delta r_{inner} \\
-               = & 1,  r > r_{inner} + \Delta r_{inner}
+               = & 0,  r \leq S_{inner} - D_{inner} \\
+  f_{inner}(r) = & \frac{1}{2}(1 - \cos(\frac{\pi}{2} (1 + \frac{r-S_{inner}}{D_{inner}})), S_{inner} - D_{inner} < r \leq S_{inner} + D_{inner} \\
+               = & 1,  r > S_{inner} + D_{inner}
 
-The values of :math:`r_{inner}` and :math:`\Delta r_{inner}` are
-the arithmetic means of the values for the central atom of type I
-and the neighbor atom of type J.
+where the switching region is centered at :math:`S_{inner}` and it extends a distance :math:`D_{inner}`
+to the left and to the right of this.
+With this option, additional keywords *sinner* and *dinner* must be used,
+each followed by *ntypes*
+values for :math:`S_{inner}` and :math:`D_{inner}`, respectively.
+When the central atom and the neighbor atom have different types,
+the values of :math:`S_{inner}` and :math:`D_{inner}` are
+the arithmetic means of the values for both types.
 
 .. note::
 
@@ -450,7 +458,7 @@ Default
 
 The optional keyword defaults are *rmin0* = 0,
 *switchflag* = 1, *bzeroflag* = 1, *quadraticflag* = 0,
-*bnormflag* = 0, *wselfallflag* = 0
+*bnormflag* = 0, *wselfallflag* = 0, *switchinnerflag* = 0,
 
 ----------
 

doc/src/compute_tally.rst

@@ -198,7 +198,9 @@ potentials only include the pair potential portion of the EAM
 interaction when used by this compute, not the embedding term.  Also
 bonded or Kspace interactions do not contribute to this compute.
 
-The computes in this package are not compatible with dynamic groups.
+When used with dynamic groups, a :doc:`run 0 <run>` command needs to
+be inserted in order to initialize the dynamic groups before accessing
+the computes.
 
 Related commands
 """"""""""""""""

doc/src/compute_temp_profile.rst

@@ -76,21 +76,28 @@ velocity for each atom.  Note that if there is only one atom in the
 bin, its thermal velocity will thus be 0.0.
 
 After the spatially-averaged velocity field has been subtracted from
-each atom, the temperature is calculated by the formula KE = (dim\*N
-- dim\*Nx\*Ny\*Nz) k T/2, where KE = total kinetic energy of the group of
-atoms (sum of 1/2 m v\^2), dim = 2 or 3 = dimensionality of the
-simulation, N = number of atoms in the group, k = Boltzmann constant,
-and T = temperature.  The dim\*Nx\*Ny\*Nz term are degrees of freedom
-subtracted to adjust for the removal of the center-of-mass velocity in
-each of Nx\*Ny\*Nz bins, as discussed in the :ref:`(Evans) <Evans1>` paper.
+each atom, the temperature is calculated by the formula
+*KE* = (*dim\*N* - *Ns\*Nx\*Ny\*Nz* - *extra* ) *k* *T*/2, where *KE* = total
+kinetic energy of the group of atoms (sum of 1/2 *m* *v*\^2), *dim* = 2
+or 3 = dimensionality of the simulation, *Ns* = 0, 1, 2 or 3 for
+streaming velocity subtracted in 0, 1, 2 or 3 dimensions, *extra* = extra
+degrees-of-freedom, *N* = number of atoms in the group, *k* = Boltzmann
+constant, and *T* = temperature.  The *Ns\*Nx\*Ny\*Nz* term is degrees
+of freedom subtracted to adjust for the removal of the center-of-mass
+velocity in each direction of the *Nx\*Ny\*Nz* bins, as discussed in the
+:ref:`(Evans) <Evans1>` paper.  The extra term defaults to (*dim* - *Ns*)
+and accounts for overall conservation of center-of-mass velocity across
+the group in directions where streaming velocity is *not* subtracted. This
+can be altered using the *extra* option of the
+:doc:`compute_modify <compute_modify>` command.
 
 If the *out* keyword is used with a *tensor* value, which is the
 default, a kinetic energy tensor, stored as a 6-element vector, is
 also calculated by this compute for use in the computation of a
 pressure tensor.  The formula for the components of the tensor is the
-same as the above formula, except that v\^2 is replaced by vx\*vy for
-the xy component, etc.  The 6 components of the vector are ordered xx,
-yy, zz, xy, xz, yz.
+same as the above formula, except that *v*\^2 is replaced by *vx\*vy* for
+the xy component, etc.  The 6 components of the vector are ordered *xx,
+yy, zz, xy, xz, yz.*
 
 If the *out* keyword is used with a *bin* value, the count of atoms
 and computed temperature for each bin are stored for output, as an
@@ -123,10 +130,20 @@ needed, the subtracted degrees-of-freedom can be altered using the
 .. note::
 
    When using the *out* keyword with a value of *bin*, the
-   calculated temperature for each bin does not include the
-   degrees-of-freedom adjustment described in the preceding paragraph,
-   for fixes that constrain molecular motion.  It does include the
-   adjustment due to the *extra* option, which is applied to each bin.
+   calculated temperature for each bin includes the degrees-of-freedom
+   adjustment described in the preceding paragraph for fixes that
+   constrain molecular motion, as well as the adjustment due to
+   the *extra* option (which defaults to *dim* - *Ns* as described above),
+   by fractionally applying them based on the fraction of atoms in each
+   bin. As a result, the bin degrees-of-freedom summed over all bins exactly
+   equals the degrees-of-freedom used in the scalar temperature calculation,
+   :math:`\Sigma N_{DOF_i} = N_{DOF}` and the corresponding relation for temperature
+   is also satisfied :math:`\Sigma N_{DOF_i} T_i = N_{DOF} T`.
+   These relations will breakdown in cases where the adjustment
+   exceeds the actual number of degrees-of-freedom in a bin. This could happen
+   if a bin is empty or in situations where rigid molecules
+   are non-uniformly distributed, in which case the reported
+   temperature within a bin may not be accurate.
 
 See the :doc:`Howto thermostat <Howto_thermostat>` page for a
 discussion of different ways to compute temperature and perform

doc/src/compute_viscosity_cos.rst

@@ -21,12 +21,13 @@ Examples
 
 .. code-block:: LAMMPS
 
+   units    real
    compute  cos all viscosity/cos
    variable V equal c_cos[7]
-   variable A equal 0.02E-5
+   variable A equal 0.02E-5  # A/fs^2
    variable density equal density
    variable lz equal lz
-   variable reciprocalViscosity equal v_V/${A}/v_density*39.4784/v_lz/v_lz*100
+   variable reciprocalViscosity equal v_V/${A}/v_density*39.4784/v_lz/v_lz*100  # 1/(Pa*s)
 
 Description
 """""""""""

doc/src/create_atoms.rst

@@ -11,7 +11,7 @@ Syntax
    create_atoms type style args keyword values ...
 
 * type = atom type (1-Ntypes) of atoms to create (offset for molecule creation)
-* style = *box* or *region* or *single* or *random*
+* style = *box* or *region* or *single* or *mesh* or *random*
 
   .. parsed-literal::
 
@@ -20,17 +20,19 @@ Syntax
          region-ID = particles will only be created if contained in the region
        *single* args = x y z
          x,y,z = coordinates of a single particle (distance units)
+       *mesh* args = STL-file
+         STL-file = file with triangle mesh in STL format
        *random* args = N seed region-ID
          N = number of particles to create
          seed = random # seed (positive integer)
          region-ID = create atoms within this region, use NULL for entire simulation box
 
 * zero or more keyword/value pairs may be appended
-* keyword = *mol* or *basis* or *ratio* or *subset* or *remap* or *var* or *set* or *rotate* or *units*
+* keyword = *mol* or *basis* or *ratio* or *subset* or *remap* or *var* or *set* or *rotate* or *overlap* or *maxtry* or *units*
 
   .. parsed-literal::
 
-       *mol* value = template-ID seed
+       *mol* values = template-ID seed
          template-ID = ID of molecule template specified in a separate :doc:`molecule <molecule>` command
          seed = random # seed (positive integer)
        *basis* values = M itype
@@ -47,9 +49,21 @@ Syntax
        *set* values = dim name
          dim = *x* or *y* or *z*
          name = name of variable to set with x, y, or z atom position
+       *radscale* value = factor
+         factor = scale factor for setting atom radius
+       *meshmode* values = mode arg
+         mode = *bisect* or *qrand*
+         *bisect* arg = radthresh
+           radthresh = threshold value for *mesh* to determine when to split triangles (distance units)
+         *qrand* arg = density
+           density = minimum number density for atoms place on *mesh* triangles (inverse distance squared units)
        *rotate* values = theta Rx Ry Rz
          theta = rotation angle for single molecule (degrees)
          Rx,Ry,Rz = rotation vector for single molecule
+       *overlap* value = Doverlap
+         Doverlap = only insert if at least this distance from all existing atoms
+       *maxtry* value = Ntry
+         Ntry = number of attempts to insert a particle before failure
        *units* value = *lattice* or *box*
          *lattice* = the geometry is defined in lattice units
          *box* = the geometry is defined in simulation box units
@@ -64,17 +78,21 @@ Examples
    create_atoms 3 region regsphere basis 2 3 ratio 0.5 74637
    create_atoms 3 single 0 0 5
    create_atoms 1 box var v set x xpos set y ypos
+   create_atoms 2 random 50 12345 NULL overlap 2.0 maxtry 50
+   create_atoms 1 mesh open_box.stl meshmode qrand 0.1 units box
+   create_atoms 1 mesh funnel.stl meshmode bisect 4.0 units box radscale 0.9
 
 Description
 """""""""""
 
-This command creates atoms (or molecules) on a lattice, or a single
-atom (or molecule), or a random collection of atoms (or molecules), as
-an alternative to reading in their coordinates explicitly via a
-:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
-command.  A simulation box must already exist, which is typically
-created via the :doc:`create_box <create_box>` command.  Before using
-this command, a lattice must also be defined using the
+This command creates atoms (or molecules) within the simulation box,
+either on a lattice, or a single atom (or molecule), or on a surface
+defined by a triangulated mesh, or a random collection of atoms (or
+molecules).  It is an alternative to reading in atom coordinates
+explicitly via a :doc:`read_data <read_data>` or :doc:`read_restart
+<read_restart>` command.  A simulation box must already exist, which is
+typically created via the :doc:`create_box <create_box>` command.
+Before using this command, a lattice must also be defined using the
 :doc:`lattice <lattice>` command, unless you specify the *single* style
 with units = box or the *random* style.  For the remainder of this doc
 page, a created atom or molecule is referred to as a "particle".
@@ -97,58 +115,126 @@ particular dimension, LAMMPS is careful to put exactly one particle at
 the boundary (on either side of the box), not zero or two.
 
 For the *region* style, a geometric volume is filled with particles on
-the lattice.  This volume what is inside the simulation box and is
-also consistent with the region volume.  See the :doc:`region <region>`
-command for details.  Note that a region can be specified so that its
-"volume" is either inside or outside a geometric boundary.  Also note
-that if your region is the same size as a periodic simulation box (in
-some dimension), LAMMPS does not implement the same logic described
-above as for the *box* style, to insure exactly one particle at
-periodic boundaries.  if this is what you desire, you should either
-use the *box* style, or tweak the region size to get precisely the
-particles you want.
+the lattice.  This volume is what is both inside the simulation box
+and also consistent with the region volume.  See the :doc:`region
+<region>` command for details.  Note that a region can be specified so
+that its "volume" is either inside or outside its geometric boundary.
+Also note that if a region is the same size as a periodic simulation
+box (in some dimension), LAMMPS does NOT implement the same logic
+described above for the *box* style, to insure exactly one particle at
+periodic boundaries.  If this is desired, you should either use the
+*box* style, or tweak the region size to get precisely the particles
+you want.
 
 For the *single* style, a single particle is added to the system at
 the specified coordinates.  This can be useful for debugging purposes
 or to create a tiny system with a handful of particles at specified
 positions.
 
-For the *random* style, N particles are added to the system at
-randomly generated coordinates, which can be useful for generating an
-amorphous system.  The particles are created one by one using the
-specified random number *seed*, resulting in the same set of particles
-coordinates, independent of how many processors are being used in the
-simulation.  If the *region-ID* argument is specified as NULL, then
-the created particles will be anywhere in the simulation box.  If a
-*region-ID* is specified, a geometric volume is filled which is both
-inside the simulation box and is also consistent with the region
-volume.  See the :doc:`region <region>` command for details.  Note that
-a region can be specified so that its "volume" is either inside or
-outside a geometric boundary.
+.. figure:: img/marble_race.jpg
+            :figwidth: 33%
+            :align: right
+            :target: _images/marble_race.jpg
+
+For the *mesh* style, a file with a triangle mesh in `STL format
+<https://en.wikipedia.org/wiki/STL_(file_format)>`_ is read and one or
+more particles are placed into the area of each triangle.  The reader
+supports both ASCII and binary files conforming to the format on the
+Wikipedia page.  Binary STL files (e.g. as frequently offered for
+3d-printing) can also be first converted to ASCII for editing with the
+:ref:`stl_bin2txt tool <stlconvert>`.  The use of the *units box* option
+is required. There are two algorithms available for placing atoms:
+*bisect* and *qrand*. They can be selected via the *meshmode* option;
+*bisect* is the default.  If the atom style allows it, the radius will
+be set to a value depending on the algorithm and the value of the
+*radscale* parameter (see below), and the atoms created from the mesh
+are assigned a new molecule ID.
+
+In *bisect* mode a particle is created at the center of each triangle
+unless the average distance of the triangle vertices from its center is
+larger than the *radthresh* value (default is lattice spacing in
+x-direction).  In case the average distance is over the threshold, the
+triangle is recursively split into two halves along the the longest side
+until the threshold is reached. There will be at least one sphere per
+triangle. The value of *radthresh* is set as an argument to *meshmode
+bisect*.  The average distance of the vertices from the center is also
+used to set the radius.
+
+In *qrand* mode a quasi-random sequence is used to distribute particles
+on mesh triangles using an approach by :ref:`(Roberts) <Roberts2019>`.
+Particles are added to the triangle until the minimum number density is
+met or exceeded such that every triangle will have at least one
+particle.  The minimum number density is set as an argument to the
+*qrand* option.  The radius will be set so that the sum of the area of
+the radius of the particles created in place of a triangle will be equal
+to the area of that triangle.
 
 .. note::
 
-   Particles generated by the *random* style will typically be
-   highly overlapped which will cause many interatomic potentials to
-   compute large energies and forces.  Thus you should either perform an
-   :doc:`energy minimization <minimize>` or run dynamics with :doc:`fix nve/limit <fix_nve_limit>` to equilibrate such a system, before
-   running normal dynamics.
+   The atom placement algorithms in the *mesh* style benefit from meshes
+   where triangles are close to equilateral.  It is therefore
+   recommended to pre-process STL files to optimize the mesh
+   accordingly.  There are multiple open source and commercial software
+   tools available with the capability to generate optimized meshes.
 
-Note that this command adds particles to those that already exist.
-This means it can be used to add particles to a system previously read
-in from a data or restart file.  Or the create_atoms command can be
-used multiple times, to add multiple sets of particles to the
-simulation.  For example, grain boundaries can be created, by
-interleaving create_atoms with :doc:`lattice <lattice>` commands
-specifying different orientations.  By using the create_atoms command
-in conjunction with the :doc:`delete_atoms <delete_atoms>` command,
-reasonably complex geometries can be created, or a protein can be
-solvated with a surrounding box of water molecules.
+.. note::
 
-In all these cases, care should be taken to insure that new atoms do
-not overlap existing atoms inappropriately, especially if molecules
-are being added.  The :doc:`delete_atoms <delete_atoms>` command can be
-used to remove overlapping atoms or molecules.
+   In most cases the atoms created in *mesh* style will become an
+   immobile or rigid object that would not be time integrated or moved
+   by :doc:`fix move <fix_move>` or :doc:`fix rigid <fix_rigid>`.  For
+   computational efficiency *and* to avoid undesired contributions to
+   pressure and potential energy due to close contacts, it is usually
+   beneficial to exclude computing interactions between the created
+   particles using :doc:`neigh_modify exclude <neigh_modify>`.
+
+For the *random* style, *N* particles are added to the system at
+randomly generated coordinates, which can be useful for generating an
+amorphous system.  The particles are created one by one using the
+specified random number *seed*, resulting in the same set of particle
+coordinates, independent of how many processors are being used in the
+simulation.  Unless the *overlap* keyword is specified, particles
+created by the *random* style will typically be highly overlapped.
+Various additional criteria can be used to accept or reject a random
+particle insertion; see the keyword discussion below.  Multiple
+attempts per particle are made (see the *maxtry* keyword) until the
+insertion is either successful or fails.  If this command fails to add
+all requested *N* particles, a warning will be output.
+
+If the *region-ID* argument is specified as NULL, then the randomly
+created particles will be anywhere in the simulation box.  If a
+*region-ID* is specified, a geometric volume is filled which is both
+inside the simulation box and is also consistent with the region
+volume.  See the :doc:`region <region>` command for details.  Note
+that a region can be specified so that its "volume" is either inside
+or outside its geometric boundary.
+
+Note that the create_atoms command adds particles to those that
+already exist.  This means it can be used to add particles to a system
+previously read in from a data or restart file.  Or the create_atoms
+command can be used multiple times, to add multiple sets of particles
+to the simulation.  For example, grain boundaries can be created, by
+interleaving the create_atoms command with :doc:`lattice <lattice>`
+commands specifying different orientations.
+
+When this command is used, care should be taken to insure the
+resulting system does not contain particles which are highly
+overlapped.  Such overlaps will cause many interatomic potentials to
+compute huge energies and forces, leading to bad dynamics.  There are
+several strategies to avoid this problem:
+
+* Use the :doc:`delete_atoms overlap <delete_atoms>` command after
+  create_atoms.  For example, this strategy can be used to overlay and
+  surround a large protein molecule with a volume of water molecules,
+  then delete water molecules that overlap with the protein atoms.
+
+* For the *random* style, use the optional *overlap* keyword to avoid
+  overlaps when each new particle is created.
+
+* Before running dynamics on an overlapped system, perform an
+  :doc:`energy minimization <minimize>`.  Or run initial dynamics with
+  :doc:`pair_style soft <pair_soft>` or with :doc:`fix nve/limit
+  <fix_nve_limit>` to un-overlap the particles, before running normal
+  dynamics.
 
 .. note::
 
@@ -156,12 +242,13 @@ used to remove overlapping atoms or molecules.
    that are outside the simulation box; they will just be ignored by
    LAMMPS.  This is true even if you are using shrink-wrapped box
    boundaries, as specified by the :doc:`boundary <boundary>` command.
-   However, you can first use the :doc:`change_box <change_box>` command to
-   temporarily expand the box, then add atoms via create_atoms, then
-   finally use change_box command again if needed to re-shrink-wrap the
-   new atoms.  See the :doc:`change_box <change_box>` page for an
-   example of how to do this, using the create_atoms *single* style to
-   insert a new atom outside the current simulation box.
+   However, you can first use the :doc:`change_box <change_box>`
+   command to temporarily expand the box, then add atoms via
+   create_atoms, then finally use change_box command again if needed
+   to re-shrink-wrap the new atoms.  See the :doc:`change_box
+   <change_box>` doc page for an example of how to do this, using the
+   create_atoms *single* style to insert a new atom outside the
+   current simulation box.
 
 ----------
 
@@ -180,17 +267,19 @@ Using a lattice to add molecules, e.g. via the *box* or *region* or
 points, except that entire molecules are added at each point, i.e. on
 the point defined by each basis atom in the unit cell as it tiles the
 simulation box or region.  This is done by placing the geometric
-center of the molecule at the lattice point, and giving the molecule a
-random orientation about the point.  The random *seed* specified with
-the *mol* keyword is used for this operation, and the random numbers
-generated by each processor are different.  This means the coordinates
-of individual atoms (in the molecules) will be different when running
-on different numbers of processors, unlike when atoms are being
-created in parallel.
+center of the molecule at the lattice point, and (by default) giving
+the molecule a random orientation about the point.  The random *seed*
+specified with the *mol* keyword is used for this operation, and the
+random numbers generated by each processor are different.  This means
+the coordinates of individual atoms (in the molecules) will be
+different when running on different numbers of processors, unlike when
+atoms are being created in parallel.
 
-Also note that because of the random rotations, it may be important to
-use a lattice with a large enough spacing that adjacent molecules will
-not overlap, regardless of their relative orientations.
+Note that with random rotations, it may be important to use a lattice
+with a large enough spacing that adjacent molecules will not overlap,
+regardless of their relative orientations.  See the description of the
+*rotate* keyword below, which overrides the default random orientation
+and inserts all molecules at a specified orientation.
 
 .. note::
 
@@ -204,7 +293,7 @@ not overlap, regardless of their relative orientations.
 
 ----------
 
-This is the meaning of the other allowed keywords.
+This is the meaning of the other optional keywords.
 
 The *basis* keyword is only used when atoms (not molecules) are being
 created.  It specifies an atom type that will be assigned to specific
@@ -234,18 +323,24 @@ and no particle is created if its position is outside the box.
 
 The *var* and *set* keywords can be used together to provide a
 criterion for accepting or rejecting the addition of an individual
-atom, based on its coordinates.  The *name* specified for the *var*
-keyword is the name of an :doc:`equal-style variable <variable>` which
-should evaluate to a zero or non-zero value based on one or two or
-three variables which will store the x, y, or z coordinates of an atom
-(one variable per coordinate).  If used, these other variables must be
-:doc:`internal-style variables <variable>` defined in the input script;
-their initial numeric value can be anything.  They must be
+atom, based on its coordinates.  They apply to all styles except
+*single*.  The *name* specified for the *var* keyword is the name of
+an :doc:`equal-style variable <variable>` which should evaluate to a
+zero or non-zero value based on one or two or three variables which
+will store the x, y, or z coordinates of an atom (one variable per
+coordinate).  If used, these other variables must be
+:doc:`internal-style variables <variable>` defined in the input
+script; their initial numeric value can be anything.  They must be
 internal-style variables, because this command resets their values
 directly.  The *set* keyword is used to identify the names of these
 other variables, one variable for the x-coordinate of a created atom,
 one for y, and one for z.
 
+.. figure:: img/sinusoid.jpg
+            :figwidth: 50%
+            :align: right
+            :target: _images/sinusoid.jpg
+
 When an atom is created, its x,y,z coordinates become the values for
 any *set* variable that is defined.  The *var* variable is then
 evaluated.  If the returned value is 0.0, the atom is not created.  If
@@ -259,6 +354,10 @@ the sinusoid would appear to be "smoother".  Also note the use of the
 "xlat" and "ylat" :doc:`thermo_style <thermo_style>` keywords which
 converts lattice spacings to distance.
 
+.. only:: html
+
+   (Click on the image for a larger version)
+
 .. code-block:: LAMMPS
 
    dimension   2
@@ -274,13 +373,7 @@ converts lattice spacings to distance.
    create_atoms  1 box var v set x xx set y yy
    write_dump  all atom sinusoid.lammpstrj
 
-.. image:: img/sinusoid.jpg
-   :scale: 50%
-   :align: center
-
-.. raw:: html
-
-   Click on the image for a larger version.
+-----
 
 The *rotate* keyword allows specification of the orientation
 at which molecules are inserted.  The axis of rotation is
@@ -291,10 +384,79 @@ the atoms around the rotation axis is consistent with the right-hand
 rule: if your right-hand's thumb points along *R*, then your fingers
 wrap around the axis in the direction of rotation.
 
+The *radscale* keyword only applies to the *mesh* style and adjusts the
+radius of created particles (see above), provided this is supported by
+the atom style.  Its value is a prefactor (must be > 0.0, default is
+1.0) that is applied to the atom radius inferred from the size of the
+individual triangles in the triangle mesh that the particle corresponds
+to.
+
+The *overlap* keyword only applies to the *random* style.  It prevents
+newly created particles from being created closer than the specified
+*Doverlap* distance from any other particle.  When the particles being
+created are molecules, the radius of the molecule (from its geometric
+center) is added to *Doverlap*.  If particles have finite size (see
+:doc:`atom_style sphere <atom_style>` for example) *Doverlap* should
+be specified large enough to include the particle size in the
+non-overlapping criterion.
+
+.. note::
+
+   Checking for overlaps is a costly O(N(N+M)) operation for inserting
+   *N* new particles into a system with *M* existing particles.  This
+   is because distances to all *M* existing particles are computed for
+   each new particle that is added.  Thus the intended use of this
+   keyword is to add relatively small numbers of particles to systems
+   which remain at a relatively low density even after the new
+   particles are created.  Careful use of the *maxtry* keyword in
+   combination with *overlap* is recommended.  See the discussion
+   above about systems with overlapped particles for alternate
+   strategies that allow for overlapped insertions.
+
+The *maxtry* keyword only applies to the *random* style.  It limits
+the number of attempts to generate valid coordinates for a single new
+particle that satisfy all requirements imposed by the *region*, *var*,
+and *overlap* keywords.  The default is 10 attempts per particle
+before the loop over the requested *N* particles advances to the next
+particle.  Note that if insertion success is unlikely (e.g. inserting
+new particles into a dense system using the *overlap* keyword),
+setting the *maxtry* keyword to a large value may result in this
+command running for a long time.
+
+.. figure:: img/overlap.png
+            :figwidth: 30%
+            :align: right
+            :target: _images/overlap.png
+
+Here is an example for the *random* style using these commands
+
+.. code-block:: LAMMPS
+
+   units         lj
+   dimension     2
+   region        box block 0 50 0 50 -0.5 0.5
+   create_box    1 box
+   create_atoms  1 random 2000 13487 NULL overlap 1.0 maxtry 50
+   pair_style    lj/cut 2.5
+   pair_coeff    1 1 1.0 1.0 2.5
+
+to produce a system as shown in the image with 1520 particles (out of
+2000 requested) that are moderately dense and which have no overlaps
+sufficient to prevent the LJ pair_style from running properly (because
+the overlap criterion = 1.0).  The create_atoms command ran for 0.3 s
+on a single CPU core.
+
+.. only:: html
+
+   (Click on the image for a larger version)
+
+-----
+
 The *units* keyword determines the meaning of the distance units used
 to specify the coordinates of the one particle created by the *single*
-style.  A *box* value selects standard distance units as defined by
-the :doc:`units <units>` command, e.g. Angstroms for units = real or
+style, or the overlap distance *Doverlap* by the *overlap* keyword.  A
+*box* value selects standard distance units as defined by the
+:doc:`units <units>` command, e.g. Angstroms for units = real or
 metal.  A *lattice* value means the distance units are in lattice
 spacings.
 
@@ -315,9 +477,10 @@ assigned to created molecules in a similar fashion.
 
 Aside from their ID, atom type, and xyz position, other properties of
 created atoms are set to default values, depending on which quantities
-are defined by the chosen :doc:`atom style <atom_style>`.  See the :doc:`atom style <atom_style>` command for more details.  See the
-:doc:`set <set>` and :doc:`velocity <velocity>` commands for info on how
-to change these values.
+are defined by the chosen :doc:`atom style <atom_style>`.  See the
+:doc:`atom style <atom_style>` command for more details.  See the
+:doc:`set <set>` and :doc:`velocity <velocity>` commands for info on
+how to change these values.
 
 * charge = 0.0
 * dipole moment magnitude = 0.0
@@ -336,9 +499,11 @@ values specified in the file read by the :doc:`molecule <molecule>`
 command.  E.g. the file typically defines bonds (angles,etc) between
 atoms in the molecule, and can optionally define charges on each atom.
 
-Note that the *sphere* atom style sets the default particle diameter
-to 1.0 as well as the density.  This means the mass for the particle
-is not 1.0, but is PI/6 \* diameter\^3 = 0.5236.
+Note that the *sphere* atom style sets the default particle diameter to
+1.0 as well as the density.  This means the mass for the particle is not
+1.0, but is PI/6 \* diameter\^3 = 0.5236.  When using the *mesh* style,
+the particle diameter is adjusted from the size of the individual
+triangles in the triangle mesh.
 
 Note that the *ellipsoid* atom style sets the default particle shape
 to (0.0 0.0 0.0) and the density to 1.0 which means it is a point
@@ -372,5 +537,13 @@ Default
 
 The default for the *basis* keyword is that all created atoms are
 assigned the argument *type* as their atom type (when single atoms are
-being created).  The other defaults are *remap* = no, *rotate* =
-random, and *units* = lattice.
+being created).  The other defaults are *remap* = no, *rotate* = random,
+*radscale* = 1.0, *radthresh* = x-lattice spacing, *overlap* not
+checked, *maxtry* = 10, and *units* = lattice.
+
+----------
+
+.. _Roberts2019:
+
+**(Roberts)** R. Roberts (2019) "Evenly Distributing Points in a Triangle." Extreme Learning.
+`<http://extremelearning.com.au/evenly-distributing-points-in-a-triangle/>`_

doc/src/delete_atoms.rst

@@ -10,7 +10,7 @@ Syntax
 
    delete_atoms style args keyword value ...
 
-* style = *group* or *region* or *overlap* or *porosity*
+* style = *group* or *region* or *overlap* or *random* or *variable*
 
   .. parsed-literal::
 
@@ -20,12 +20,19 @@ Syntax
          cutoff = delete one atom from pairs of atoms within the cutoff (distance units)
          group1-ID = one atom in pair must be in this group
          group2-ID = other atom in pair must be in this group
-       *porosity* args = group-ID region-ID fraction seed
+       *random* args = ranstyle value eflag group-ID region-ID seed
+         ranstyle = *fraction* or *count*
+           for *fraction*:
+             value = fraction (0.0 to 1.0) of eligible atoms to delete
+             eflag =  *no* for fast approximate deletion, *yes* for exact deletion
+           for *count*:
+             value = number of atoms to delete
+             eflag = *no* for warning if count > eligible atoms, *yes* for error
          group-ID = group within which to perform deletions
          region-ID = region within which to perform deletions
                      or NULL to only impose the group criterion
-         fraction = delete this fraction of atoms
          seed = random number seed (positive integer)
+       *variable* args = variable-name
 
 * zero or more keyword/value pairs may be appended
 * keyword = *compress* or *bond* or *mol*
@@ -45,15 +52,17 @@ Examples
    delete_atoms region sphere compress no
    delete_atoms overlap 0.3 all all
    delete_atoms overlap 0.5 solvent colloid
-   delete_atoms porosity all cube 0.1 482793 bond yes
-   delete_atoms porosity polymer cube 0.1 482793 bond yes
+   delete_atoms random fraction 0.1 yes all cube 482793 bond yes
+   delete_atoms random fraction 0.3 no polymer NULL 482793 bond yes
+   delete_atoms random count 500 no ions NULL 482793
+   detele_atoms variable checkers
 
 Description
 """""""""""
 
-Delete the specified atoms.  This command can be used to carve out
-voids from a block of material or to delete created atoms that are too
-close to each other (e.g. at a grain boundary).
+Delete the specified atoms.  This command can be used, for example, to
+carve out voids from a block of material or to delete created atoms
+that are too close to each other (e.g. at a grain boundary).
 
 For style *group*, all atoms belonging to the group are deleted.
 
@@ -79,17 +88,44 @@ have occurred that no atom pairs within the cutoff will remain
 minimum number of atoms will be deleted, or that the same atoms will
 be deleted when running on different numbers of processors.
 
-For style *porosity* a specified *fraction* of atoms are deleted which
-are both in the specified group and within the specified region.  The
-region-ID can be specified as NULL to only impose the group criterion.
-Likewise, specifying the group-ID as *all* will only impose the region
-criterion.
+For style *random* a subset of eligible atoms are deleted.  Which
+atoms to delete are chosen randomly using the specified random number
+*seed*.  Which atoms are deleted may vary when running on different
+numbers of processors.
 
-For example, if fraction is 0.1, then 10% of the eligible atoms will
-be deleted.  The atoms to delete are chosen randomly.  There is no
-guarantee that the exact fraction of atoms will be deleted, or that
-the same atoms will be deleted when running on different numbers of
-processors.
+For *ranstyle* = *fraction*, the specified fractional *value* (0.0 to
+1.0) of eligible atoms are deleted.  If *eflag* is set to *no*, then
+the number of deleted atoms will be approximate, but the operation
+will be fast.  If *eflag* is set to *yes*, then the number deleted
+will match the requested fraction, but for large systems the selection
+of deleted atoms may take additional time to determine.
+
+For *ranstyle* = *count*, the specified integer *value* is the number
+of eligible atoms are deleted.  If *eflag* is set to *no*, then if the
+requested number is larger then the number of eligible atoms, a
+warning is issued and only the eligible atoms are deleted instead of
+the requested *value*.  If *eflag* is set to *yes*, an error is
+triggered instead and LAMMPS will exit.  For large systems the
+selection of atoms to delete may take additional time to determine,
+the same as for requesting an exact fraction with *pstyle* =
+*fraction*.
+
+Which atoms are eligible for deletion for style *random* is determined
+by the specified *group-ID* and *region-ID*.  To be eligible, an atom
+must be in both the specified group and region.  If *group-ID* = all,
+there is effectively no group criterion.  If *region-ID* is specified
+as NULL, no region criterion is imposed.
+
+For style *variable*, all atoms for which the atom-style variable with
+the given name evaluates to non-zero will be deleted. Additional atoms
+can be deleted if they are in a molecule for which one or more atoms
+were deleted within the region; see the *mol* keyword discussion below.
+This option allows complex selections of atoms not covered by the
+other options listed above.
+
+----------
+
+Here is the meaning of the optional keywords.
 
 If the *compress* keyword is set to *yes*, then after atoms are
 deleted, then atom IDs are re-assigned so that they run from 1 to the