Merge pull request #2587 from akohlmey/next-patch-version

Step version strings for next patch release
Merge pull request #2596 from lammps/create-bonds-bugs
2021-02-10 21:30:20 -05:00 · 2021-02-10 20:35:55 -05:00 · 2021-02-10 20:14:19 -05:00 · 2021-02-10 18:40:25 -05:00 · 2021-02-10 16:21:25 -07:00 · 2021-02-10 18:19:21 -05:00
2164 changed files with 76581 additions and 69394 deletions
--- a/.github/workflows/codeql-analysis.yml
+++ b/.github/workflows/codeql-analysis.yml
@ -0,0 +1,47 @@
+# GitHub action to run static code analysis on C++ and Python code
+name: "CodeQL Code Analysis"
+
+on:
+  push:
+    branches: [master]
+
+jobs:
+  analyze:
+    name: Analyze
+    if: ${{ github.repository == 'lammps/lammps' }}
+    runs-on: ubuntu-latest
+
+    strategy:
+      fail-fast: false
+      matrix:
+        language: ['cpp', 'python']
+
+    steps:
+    - name: Checkout repository
+      uses: actions/checkout@v2
+      with:
+        fetch-depth: 2
+
+    - name: Setup Python
+      uses: actions/setup-python@v2
+      with:
+        python-version: '3.x'
+
+    - name: Initialize CodeQL
+      uses: github/codeql-action/init@v1
+      with:
+        languages: ${{ matrix.language }}
+
+    - name: Create Build Environment
+      run: cmake -E make_directory ${{github.workspace}}/build
+
+    - name: Building LAMMPS via CMake
+      if: ${{ matrix.language == 'cpp' }}
+      shell: bash
+      working-directory: ${{github.workspace}}/build
+      run: |
+        cmake -C $GITHUB_WORKSPACE/cmake/presets/most.cmake $GITHUB_WORKSPACE/cmake
+        cmake --build . --parallel 2
+
+    - name: Perform CodeQL Analysis
+      uses: github/codeql-action/analyze@v1
--- a/.github/workflows/unittest-macos.yml
+++ b/.github/workflows/unittest-macos.yml
@ -0,0 +1,34 @@
+# GitHub action to build LAMMPS on MacOS and run unit tests
+name: "Unittest for MacOS"
+
+on:
+  push:
+    branches: [master]
+
+jobs:
+  build:
+    name: MacOS Unit Test
+    if: ${{ github.repository == 'lammps/lammps' }}
+    runs-on: macos-latest
+
+    steps:
+    - name: Checkout repository
+      uses: actions/checkout@v2
+      with:
+        fetch-depth: 2
+
+    - name: Create Build Environment
+      run: cmake -E make_directory ${{github.workspace}}/build
+
+    - name: Building LAMMPS via CMake
+      shell: bash
+      working-directory: ${{github.workspace}}/build
+      run: |
+        cmake -C $GITHUB_WORKSPACE/cmake/presets/most.cmake $GITHUB_WORKSPACE/cmake \
+              -DENABLE_TESTING=ON -DBUILD_SHARED_LIBS=ON -DLAMMPS_EXCEPTIONS=ON
+        cmake --build . --parallel 2
+
+    - name: Run Tests
+      working-directory: ${{github.workspace}}/build
+      shell: bash
+      run: ctest -V
--- a/cmake/CMakeLists.txt
+++ b/cmake/CMakeLists.txt
@ -25,7 +25,7 @@ set(LAMMPS_POTENTIALS_DIR ${LAMMPS_DIR}/potentials)
 find_package(Git)

 # by default, install into $HOME/.local (not /usr/local), so that no root access (and sudo!!) is needed
-if (CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
+if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
  set(CMAKE_INSTALL_PREFIX "$ENV{HOME}/.local" CACHE PATH "default install path" FORCE )
 endif()

@ -33,7 +33,7 @@ endif()
 set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/Modules)

 # make sure LIBRARY_PATH is set if environment variable is set
-if (DEFINED ENV{LIBRARY_PATH})
+if(DEFINED ENV{LIBRARY_PATH})
  list(APPEND CMAKE_LIBRARY_PATH "$ENV{LIBRARY_PATH}")
  message(STATUS "Appending $ENV{LIBRARY_PATH} to CMAKE_LIBRARY_PATH: ${CMAKE_LIBRARY_PATH}")
 endif()
@ -107,13 +107,15 @@ option(CMAKE_VERBOSE_MAKEFILE "Generate verbose Makefiles" OFF)
 set(STANDARD_PACKAGES ASPHERE BODY CLASS2 COLLOID COMPRESS DIPOLE
  GRANULAR KSPACE LATTE MANYBODY MC MESSAGE MISC MLIAP MOLECULE PERI POEMS
  QEQ REPLICA RIGID SHOCK SPIN SNAP SRD KIM PYTHON MSCG MPIIO VORONOI
-  USER-ATC USER-AWPMD USER-BOCS USER-CGDNA USER-MESODPD USER-CGSDK USER-COLVARS
-  USER-DIFFRACTION USER-DPD USER-DRUDE USER-EFF USER-FEP USER-H5MD USER-LB
-  USER-MANIFOLD USER-MEAMC USER-MESONT USER-MGPT USER-MISC USER-MOFFF USER-MOLFILE
-  USER-NETCDF USER-PHONON USER-PLUMED USER-PTM USER-QTB USER-REACTION
-  USER-REAXC USER-SCAFACOS USER-SDPD USER-SMD USER-SMTBQ USER-SPH USER-TALLY
-  USER-UEF USER-VTK USER-QUIP USER-QMMM USER-YAFF USER-ADIOS)
-set(SUFFIX_PACKAGES CORESHELL USER-OMP KOKKOS OPT USER-INTEL GPU)
+  USER-ADIOS USER-ATC USER-AWPMD USER-BOCS USER-CGDNA USER-MESODPD USER-CGSDK
+  USER-COLVARS USER-DIFFRACTION USER-DPD USER-DRUDE USER-EFF USER-FEP USER-H5MD
+  USER-LB USER-MANIFOLD USER-MEAMC USER-MESONT USER-MGPT USER-MISC USER-MOFFF
+  USER-MOLFILE USER-NETCDF USER-PHONON USER-PLUMED USER-PTM USER-QTB
+  USER-REACTION USER-REAXC USER-SCAFACOS USER-SDPD USER-SMD USER-SMTBQ USER-SPH
+  USER-TALLY USER-UEF USER-VTK USER-QUIP USER-QMMM USER-YAFF)
+
+set(SUFFIX_PACKAGES CORESHELL GPU KOKKOS OPT USER-INTEL USER-OMP)
+
 foreach(PKG ${STANDARD_PACKAGES} ${SUFFIX_PACKAGES})
  option(PKG_${PKG} "Build ${PKG} Package" OFF)
 endforeach()
@ -220,6 +222,7 @@ if(BUILD_OMP)
  endif()

  if (((CMAKE_CXX_COMPILER_ID STREQUAL "GNU") AND (CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 9.0)) OR
+      (CMAKE_CXX_COMPILER_ID STREQUAL "PGI") OR
      ((CMAKE_CXX_COMPILER_ID STREQUAL "Clang") AND (CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 10.0)) OR
      ((CMAKE_CXX_COMPILER_ID STREQUAL "Intel") AND (CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 19.0)))
    # GCC 9.x and later plus Clang 10.x and later implement strict OpenMP 4.0 semantics for consts.
@ -372,7 +375,7 @@ else()
  set(CUDA_REQUEST_PIC)
 endif()

-foreach(PKG_WITH_INCL KSPACE PYTHON VORONOI USER-COLVARS USER-MOLFILE USER-NETCDF USER-PLUMED USER-QMMM
+foreach(PKG_WITH_INCL KSPACE PYTHON MLIAP VORONOI USER-COLVARS USER-MOLFILE USER-NETCDF USER-PLUMED USER-QMMM
        USER-QUIP USER-SCAFACOS USER-SMD USER-VTK KIM LATTE MESSAGE MSCG COMPRESS)
  if(PKG_${PKG_WITH_INCL})
    include(Packages/${PKG_WITH_INCL})
@ -579,7 +582,7 @@ add_dependencies(lammps gitversion)
 ############################################
 get_property(LANGUAGES GLOBAL PROPERTY ENABLED_LANGUAGES)
 list (FIND LANGUAGES "Fortran" _index)
-if (${_index} GREATER -1)
+if(${_index} GREATER -1)
  target_link_libraries(lammps PRIVATE ${CMAKE_Fortran_IMPLICIT_LINK_LIBRARIES})
 endif()
 set(LAMMPS_CXX_HEADERS angle.h atom.h bond.h citeme.h comm.h compute.h dihedral.h domain.h error.h fix.h force.h group.h improper.h
@ -649,7 +652,7 @@ install(
 if(BUILD_SHARED_LIBS)
  if(CMAKE_VERSION VERSION_LESS 3.12)
    # adjust so we find Python 3 versions before Python 2 on old systems with old CMake
-    set(Python_ADDITIONAL_VERSIONS 3.8 3.7 3.6 3.5)
+    set(Python_ADDITIONAL_VERSIONS 3.9 3.8 3.7 3.6 3.5)
    find_package(PythonInterp) # Deprecated since version 3.12
    if(PYTHONINTERP_FOUND)
        set(Python_EXECUTABLE ${PYTHON_EXECUTABLE})
@ -659,10 +662,10 @@ if(BUILD_SHARED_LIBS)
  endif()
  if (Python_EXECUTABLE)
    add_custom_target(
-      install-python
-      ${Python_EXECUTABLE} install.py -v ${LAMMPS_SOURCE_DIR}/version.h
-      -m ${LAMMPS_PYTHON_DIR}/lammps.py
-      -l ${CMAKE_BINARY_DIR}/liblammps${CMAKE_SHARED_LIBRARY_SUFFIX}
+      install-python ${CMAKE_COMMAND} -E remove_directory build
+      COMMAND ${Python_EXECUTABLE} install.py -v ${LAMMPS_SOURCE_DIR}/version.h
+      -p ${LAMMPS_PYTHON_DIR}/lammps
+      -l ${CMAKE_BINARY_DIR}/liblammps${LAMMPS_MACHINE}${CMAKE_SHARED_LIBRARY_SUFFIX}
      WORKING_DIRECTORY  ${LAMMPS_PYTHON_DIR}
      COMMENT "Installing LAMMPS Python module")
  else()
@ -691,11 +694,8 @@ if(BUILD_SHARED_LIBS OR PKG_PYTHON)
    find_package(Python COMPONENTS Interpreter)
  endif()
  if (Python_EXECUTABLE)
-    execute_process(COMMAND ${Python_EXECUTABLE}
-      -c "import distutils.sysconfig as cg; print(cg.get_python_lib(1,0,prefix='${CMAKE_INSTALL_PREFIX}'))"
-      OUTPUT_VARIABLE PYTHON_DEFAULT_INSTDIR OUTPUT_STRIP_TRAILING_WHITESPACE)
-    set(PYTHON_INSTDIR ${PYTHON_DEFAULT_INSTDIR} CACHE PATH "Installation folder for LAMMPS Python module")
-    install(FILES ${LAMMPS_PYTHON_DIR}/lammps.py DESTINATION ${PYTHON_INSTDIR})
+    file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/python)
+    install(CODE "execute_process(COMMAND ${Python_EXECUTABLE} setup.py build -b ${CMAKE_BINARY_DIR}/python install --prefix=${CMAKE_INSTALL_PREFIX} --root=\$ENV{DESTDIR}/ WORKING_DIRECTORY ${LAMMPS_PYTHON_DIR})")
  endif()
 endif()

@ -739,14 +739,14 @@ if(OPTIONS)
 endif()
 get_property(LANGUAGES GLOBAL PROPERTY ENABLED_LANGUAGES)
 list (FIND LANGUAGES "Fortran" _index)
-if (${_index} GREATER -1)
+if(${_index} GREATER -1)
  message(STATUS "Fortran Compiler: ${CMAKE_Fortran_COMPILER}
      Type:          ${CMAKE_Fortran_COMPILER_ID}
      Version:       ${CMAKE_Fortran_COMPILER_VERSION}
      Fortran Flags:${CMAKE_Fortran_FLAGS} ${CMAKE_Fortran_FLAGS_${BTYPE}}")
 endif()
 list (FIND LANGUAGES "C" _index)
-if (${_index} GREATER -1)
+if(${_index} GREATER -1)
  message(STATUS "C compiler:       ${CMAKE_C_COMPILER}
      Type:          ${CMAKE_C_COMPILER_ID}
      Version:       ${CMAKE_C_COMPILER_VERSION}
--- a/cmake/Modules/CodingStandard.cmake
+++ b/cmake/Modules/CodingStandard.cmake
@ -8,7 +8,7 @@ else()
    find_package(Python3 COMPONENTS Interpreter QUIET)
 endif()

-if (Python3_EXECUTABLE)
+if(Python3_EXECUTABLE)
    if(Python3_VERSION VERSION_GREATER_EQUAL 3.5)
        add_custom_target(
          check-whitespace
--- a/cmake/Modules/FindCythonize.cmake
+++ b/cmake/Modules/FindCythonize.cmake
@ -0,0 +1,30 @@
+# Find the Cythonize tool.
+#
+# This code sets the following variables:
+#
+#  Cythonize_EXECUTABLE
+#
+# adapted from https://github.com/cmarshall108/cython-cmake-example/blob/master/cmake/FindCython.cmake
+#=============================================================================
+
+if(CMAKE_VERSION VERSION_LESS 3.12)
+    find_package(PythonInterp 3.6 QUIET) # Deprecated since version 3.12
+    if(PYTHONINTERP_FOUND)
+        set(Python3_EXECUTABLE ${PYTHON_EXECUTABLE})
+    endif()
+else()
+    find_package(Python3 3.6 COMPONENTS Interpreter QUIET)
+endif()
+
+# Use the Cython executable that lives next to the Python executable
+# if it is a local installation.
+if(Python3_EXECUTABLE)
+  get_filename_component(_python_path ${Python3_EXECUTABLE} PATH)
+  find_program(Cythonize_EXECUTABLE
+    NAMES cythonize3 cythonize cythonize.bat
+    HINTS ${_python_path})
+endif()
+
+include(FindPackageHandleStandardArgs)
+FIND_PACKAGE_HANDLE_STANDARD_ARGS(Cythonize REQUIRED_VARS Cythonize_EXECUTABLE)
+mark_as_advanced(Cythonize_EXECUTABLE)
--- a/cmake/Modules/LAMMPSUtils.cmake
+++ b/cmake/Modules/LAMMPSUtils.cmake
@ -50,6 +50,7 @@ function(check_for_autogen_files source_dir)
    file(GLOB SRC_AUTOGEN_FILES ${source_dir}/style_*.h)
    file(GLOB SRC_AUTOGEN_PACKAGES ${source_dir}/packages_*.h)
    list(APPEND SRC_AUTOGEN_FILES ${SRC_AUTOGEN_PACKAGES} ${source_dir}/lmpinstalledpkgs.h ${source_dir}/lmpgitversion.h)
+    list(APPEND SRC_AUTOGEN_FILES ${SRC_AUTOGEN_PACKAGES} ${source_dir}/mliap_model_python_couple.h ${source_dir}/mliap_model_python_couple.cpp)
    foreach(_SRC ${SRC_AUTOGEN_FILES})
      get_filename_component(FILENAME "${_SRC}" NAME)
      if(EXISTS ${source_dir}/${FILENAME})
--- a/cmake/Modules/MPI4WIN.cmake
+++ b/cmake/Modules/MPI4WIN.cmake
@ -1,7 +1,7 @@
 # Download and configure custom MPICH files for Windows
 message(STATUS "Downloading and configuring MPICH-1.4.1 for Windows")
 include(ExternalProject)
-if (CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64")
+if(CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64")
  ExternalProject_Add(mpi4win_build
    URL https://download.lammps.org/thirdparty/mpich2-win64-devel.tar.gz
    URL_MD5 4939fdb59d13182fd5dd65211e469f14
--- a/cmake/Modules/Packages/GPU.cmake
+++ b/cmake/Modules/Packages/GPU.cmake
@ -2,6 +2,7 @@ set(GPU_SOURCES_DIR ${LAMMPS_SOURCE_DIR}/GPU)
 set(GPU_SOURCES ${GPU_SOURCES_DIR}/gpu_extra.h
                ${GPU_SOURCES_DIR}/fix_gpu.h
                ${GPU_SOURCES_DIR}/fix_gpu.cpp)
+target_compile_definitions(lammps PRIVATE -DLMP_GPU)

 set(GPU_API "opencl" CACHE STRING "API used by GPU package")
 set(GPU_API_VALUES opencl cuda hip)
@ -35,6 +36,9 @@ if(GPU_API STREQUAL "CUDA")
  option(CUDPP_OPT "Enable CUDPP_OPT" ON)
  option(CUDA_MPS_SUPPORT "Enable tweaks to support CUDA Multi-process service (MPS)" OFF)
  if(CUDA_MPS_SUPPORT)
+    if(CUDPP_OPT)
+      message(FATAL_ERROR "Must use -DCUDPP_OPT=OFF with -DGPU_CUDA_MPS_SUPPORT=ON")
+    endif()
    set(GPU_CUDA_MPS_FLAGS "-DCUDA_PROXY")
  endif()

@ -93,9 +97,9 @@ if(GPU_API STREQUAL "CUDA")
  if(CUDA_VERSION VERSION_GREATER_EQUAL "10.0")
    string(APPEND GPU_CUDA_GENCODE " -gencode arch=compute_75,code=[sm_75,compute_75]")
  endif()
-  # Ampere (GPU Arch 8.0) is supported by CUDA 11 and later
+  # Ampere (GPU Arch 8.0 and 8.6) is supported by CUDA 11 and later
  if(CUDA_VERSION VERSION_GREATER_EQUAL "11.0")
-    string(APPEND GPU_CUDA_GENCODE " -gencode arch=compute_80,code=[sm_80,compute_80]")
+    string(APPEND GPU_CUDA_GENCODE " -gencode arch=compute_80,code=[sm_80,compute_80] -gencode arch=compute_86,code=[sm_86,compute_86]")
  endif()
  if(CUDA_VERSION VERSION_GREATER_EQUAL "12.0")
    message(WARNING "Unsupported CUDA version. Use at your own risk.")
--- a/cmake/Modules/Packages/KIM.cmake
+++ b/cmake/Modules/Packages/KIM.cmake
@ -19,6 +19,8 @@ if(CURL_FOUND)
    target_compile_definitions(lammps PRIVATE -DLMP_NO_SSL_CHECK)
  endif()
 endif()
+set(KIM_EXTRA_UNITTESTS OFF CACHE STRING "Set extra unit tests verbose mode on/off. If on, extra tests are included.")
+mark_as_advanced(KIM_EXTRA_UNITTESTS)
 find_package(PkgConfig QUIET)
 set(DOWNLOAD_KIM_DEFAULT ON)
 if(PKG_CONFIG_FOUND)
@ -34,8 +36,8 @@ if(DOWNLOAD_KIM)
  enable_language(C)
  enable_language(Fortran)
  ExternalProject_Add(kim_build
-    URL https://s3.openkim.org/kim-api/kim-api-2.2.0.txz
-    URL_MD5 e7f944e1593cffd7444679a660607f6c
+    URL https://s3.openkim.org/kim-api/kim-api-2.2.1.txz
+    URL_MD5 ae1ddda2ef7017ea07934e519d023dca
    BINARY_DIR build
    CMAKE_ARGS ${CMAKE_REQUEST_PIC}
               -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}
--- a/cmake/Modules/Packages/KOKKOS.cmake
+++ b/cmake/Modules/Packages/KOKKOS.cmake
@ -1,4 +1,7 @@
 ########################################################################
+# As of version 3.3.0 Kokkos requires C++14
+set(CMAKE_CXX_STANDARD 14)
+########################################################################
 # consistency checks and Kokkos options/settings required by LAMMPS
 if(Kokkos_ENABLE_CUDA)
  message(STATUS "KOKKOS: Enabling CUDA LAMBDA function support")
@ -35,8 +38,8 @@ if(DOWNLOAD_KOKKOS)
  list(APPEND KOKKOS_LIB_BUILD_ARGS "-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TOOLCHAIN_FILE}")
  include(ExternalProject)
  ExternalProject_Add(kokkos_build
-    URL https://github.com/kokkos/kokkos/archive/3.2.01.tar.gz
-    URL_MD5 ba72440e285ccde05b403694ea0c92e5
+    URL https://github.com/kokkos/kokkos/archive/3.3.01.tar.gz
+    URL_MD5 08201d1c7cf5bc458ce0f5b44a629d5a
    CMAKE_ARGS ${KOKKOS_LIB_BUILD_ARGS}
    BUILD_BYPRODUCTS <INSTALL_DIR>/lib/libkokkoscore.a
  )
@ -50,7 +53,7 @@ if(DOWNLOAD_KOKKOS)
  target_link_libraries(lammps PRIVATE LAMMPS::KOKKOS)
  add_dependencies(LAMMPS::KOKKOS kokkos_build)
 elseif(EXTERNAL_KOKKOS)
-  find_package(Kokkos 3.2.01 REQUIRED CONFIG)
+  find_package(Kokkos 3.3.01 REQUIRED CONFIG)
  target_link_libraries(lammps PRIVATE Kokkos::kokkos)
 else()
  set(LAMMPS_LIB_KOKKOS_SRC_DIR ${LAMMPS_LIB_SOURCE_DIR}/kokkos)
--- a/cmake/Modules/Packages/MLIAP.cmake
+++ b/cmake/Modules/Packages/MLIAP.cmake
@ -0,0 +1,40 @@
+# if PYTHON package is included we may also include Python support in MLIAP
+set(MLIAP_ENABLE_PYTHON_DEFAULT OFF)
+if(PKG_PYTHON)
+  find_package(Cythonize QUIET)
+  if(Cythonize_FOUND)
+    set(MLIAP_ENABLE_PYTHON_DEFAULT ON)
+  endif()
+endif()
+
+option(MLIAP_ENABLE_PYTHON "Build MLIAP package with Python support" ${MLIAP_ENABLE_PYTHON_DEFAULT})
+
+if(MLIAP_ENABLE_PYTHON)
+  find_package(Cythonize REQUIRED)
+  if(NOT PKG_PYTHON)
+    message(FATAL_ERROR "Must enable PYTHON package for including Python support in MLIAP")
+  endif()
+  if(CMAKE_VERSION VERSION_LESS 3.12)
+    if(PYTHONLIBS_VERSION_STRING VERSION_LESS 3.6)
+      message(FATAL_ERROR "Python support in MLIAP requires Python 3.6 or later")
+    endif()
+  else()
+    if(Python_VERSION VERSION_LESS 3.6)
+      message(FATAL_ERROR "Python support in MLIAP requires Python 3.6 or later")
+    endif()
+  endif()
+
+  set(MLIAP_BINARY_DIR ${CMAKE_BINARY_DIR}/cython)
+  set(MLIAP_CYTHON_SRC ${LAMMPS_SOURCE_DIR}/MLIAP/mliap_model_python_couple.pyx)
+  get_filename_component(MLIAP_CYTHON_BASE ${MLIAP_CYTHON_SRC} NAME_WE)
+  file(MAKE_DIRECTORY ${MLIAP_BINARY_DIR})
+  add_custom_command(OUTPUT  ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.cpp ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.h
+          COMMAND            ${CMAKE_COMMAND} -E copy_if_different ${MLIAP_CYTHON_SRC} ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.pyx
+          COMMAND            ${Cythonize_EXECUTABLE} -3 ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.pyx
+          WORKING_DIRECTORY  ${MLIAP_BINARY_DIR}
+          MAIN_DEPENDENCY    ${MLIAP_CYTHON_SRC}
+          COMMENT "Generating C++ sources with cythonize...")
+  target_compile_definitions(lammps PRIVATE -DMLIAP_PYTHON)
+  target_sources(lammps PRIVATE ${MLIAP_BINARY_DIR}/${MLIAP_CYTHON_BASE}.cpp)
+  target_include_directories(lammps PRIVATE ${MLIAP_BINARY_DIR})
+endif()
--- a/cmake/Modules/Packages/USER-PLUMED.cmake
+++ b/cmake/Modules/Packages/USER-PLUMED.cmake
@ -55,8 +55,8 @@ if(DOWNLOAD_PLUMED)
  endif()
  include(ExternalProject)
  ExternalProject_Add(plumed_build
-    URL https://github.com/plumed/plumed2/releases/download/v2.6.1/plumed-src-2.6.1.tgz
-    URL_MD5 89a9a450fc6025299fe16af235957163
+    URL https://github.com/plumed/plumed2/releases/download/v2.7.0/plumed-src-2.7.0.tgz
+    URL_MD5 95f29dd0c067577f11972ff90dfc7d12
    BUILD_IN_SOURCE 1
    CONFIGURE_COMMAND <SOURCE_DIR>/configure --prefix=<INSTALL_DIR>
                                             ${CONFIGURE_REQUEST_PIC}
--- a/cmake/presets/pgi.cmake
+++ b/cmake/presets/pgi.cmake
@ -0,0 +1,16 @@
+# preset that will enable clang/clang++ 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)
+set(CMAKE_Fortran_COMPILER "pgfortran" CACHE STRING "" FORCE)
+set(MPI_CXX "pgc++" CACHE STRING "" FORCE)
+set(MPI_CXX_COMPILER "mpicxx" CACHE STRING "" FORCE)
+unset(HAVE_OMP_H_INCLUDE CACHE)
+
+set(OpenMP_C "pgcc" CACHE STRING "" FORCE)
+set(OpenMP_C_FLAGS "-mp" CACHE STRING "" FORCE)
+set(OpenMP_C_LIB_NAMES "omp" CACHE STRING "" FORCE)
+set(OpenMP_CXX "pgc++" CACHE STRING "" FORCE)
+set(OpenMP_CXX_FLAGS "-mp" CACHE STRING "" FORCE)
+set(OpenMP_CXX_LIB_NAMES "omp" CACHE STRING "" FORCE)
+set(OpenMP_omp_LIBRARY "libomp.so" CACHE PATH "" FORCE)
--- a/doc/Makefile
+++ b/doc/Makefile
@ -94,7 +94,7 @@ $(SPHINXCONFIG)/conf.py: $(SPHINXCONFIG)/conf.py.in
 	    -e 's,@LAMMPS_PYTHON_DIR@,$(BUILDDIR)/../python,g' \
 	    -e 's,@LAMMPS_DOC_DIR@,$(BUILDDIR),g' $< > $@

-html: xmlgen $(SPHINXCONFIG)/conf.py $(ANCHORCHECK) $(MATHJAX)
+html: 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
 	@(\
@ -118,7 +118,7 @@ html: xmlgen $(SPHINXCONFIG)/conf.py $(ANCHORCHECK) $(MATHJAX)
 	@rm -rf html/PDF/.[sg]*
 	@echo "Build finished. The HTML pages are in doc/html."

-spelling: xmlgen $(VENV) $(SPHINXCONFIG)/false_positives.txt
+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
 	@(\
 		. $(VENV)/bin/activate ; env PYTHONWARNINGS= \
@ -229,7 +229,7 @@ $(VENV):
 		$(VIRTUALENV) -p $(PYTHON) $(VENV); \
 		. $(VENV)/bin/activate; \
 		pip install --upgrade pip; \
-		pip install --use-feature=2020-resolver -r $(BUILDDIR)/utils/requirements.txt; \
+		pip install -r $(BUILDDIR)/utils/requirements.txt; \
 		deactivate;\
 	)

--- a/doc/doxygen/Doxyfile.in
+++ b/doc/doxygen/Doxyfile.in
@ -424,6 +424,8 @@ INPUT                  = @LAMMPS_SOURCE_DIR@/utils.cpp                 \
                         @LAMMPS_SOURCE_DIR@/input.h                   \
                         @LAMMPS_SOURCE_DIR@/tokenizer.cpp             \
                         @LAMMPS_SOURCE_DIR@/tokenizer.h               \
+                         @LAMMPS_SOURCE_DIR@/arg_info.cpp              \
+                         @LAMMPS_SOURCE_DIR@/arg_info.h                \
                         @LAMMPS_SOURCE_DIR@/text_file_reader.cpp      \
                         @LAMMPS_SOURCE_DIR@/text_file_reader.h        \
                         @LAMMPS_SOURCE_DIR@/potential_file_reader.cpp \
--- a/doc/lammps.1
+++ b/doc/lammps.1
@ -1,4 +1,4 @@
-.TH LAMMPS "30 November 2020" "2020-10-29"
+.TH LAMMPS "10 February 2021" "2021-02-10"
 .SH NAME
 .B LAMMPS
 \- Molecular Dynamics Simulator.
--- a/doc/src/Build_basics.rst
+++ b/doc/src/Build_basics.rst
@ -236,12 +236,15 @@ LAMMPS.
         cmake ../cmake -DCMAKE_C_COMPILER=icc -DCMAKE_CXX_COMPILER=icpc -DCMAKE_Fortran_COMPILER=ifort
         # Building with LLVM/Clang Compilers:
         cmake ../cmake -DCMAKE_C_COMPILER=clang -DCMAKE_CXX_COMPILER=clang++ -DCMAKE_Fortran_COMPILER=flang
+         # Building with PGI/Nvidia Compilers:
+         cmake ../cmake -DCMAKE_C_COMPILER=pgcc -DCMAKE_CXX_COMPILER=pgc++ -DCMAKE_Fortran_COMPILER=pgfortran

      For compiling with the Clang/LLVM compilers a CMake preset is
      provided that can be loaded with
      `-C ../cmake/presets/clang.cmake`.  Similarly,
      `-C ../cmake/presets/intel.cmake` should switch the compiler
-      toolchain to the Intel compilers.
+      toolchain to the Intel compilers and `-C ../cmake/presets/pgi.cmake`
+      should switch the compiler to the PGI compilers.

      In addition you can set ``CMAKE_TUNE_FLAGS`` to specifically add
      compiler flags to tune for optimal performance on given hosts. By
--- a/doc/src/Build_extras.rst
+++ b/doc/src/Build_extras.rst
@ -37,6 +37,7 @@ This is the list of packages that may require additional steps.
   * :ref:`KOKKOS <kokkos>`
   * :ref:`LATTE <latte>`
   * :ref:`MESSAGE <message>`
+   * :ref:`MLIAP <mliap>`
   * :ref:`MSCG <mscg>`
   * :ref:`OPT <opt>`
   * :ref:`POEMS <poems>`
@ -130,7 +131,7 @@ CMake build
   -D HIP_USE_DEVICE_SORT=value # enables GPU sorting
                                # value = yes (default) or no
   -D CUDPP_OPT=value           # optimization setting for GPU_API=cuda
-                                # enables CUDA Performance Primitives Optimizations
+                                # enables CUDA Performance Primitives Optimizations, must be "no" for CUDA_MPS_SUPPORT=yes
                                # value = yes (default) or no
   -D CUDA_MPS_SUPPORT=value    # enables some tweaks required to run with active nvidia-cuda-mps daemon
                                # value = yes or no (default)
@ -218,11 +219,19 @@ Makefile if desired:
 * ``CUDA_PRECISION`` = precision (double, mixed, single)
 * ``EXTRAMAKE`` = which Makefile.lammps.\* file to copy to Makefile.lammps

-The file Makefile.linux_multi is set up to include support for multiple
+The file Makefile.cuda is set up to include support for multiple
 GPU architectures as supported by the CUDA toolkit in use. This is done
 through using the "--gencode " flag, which can be used multiple times and
 thus support all GPU architectures supported by your CUDA compiler.

+To include CUDA performance primitives set the Makefile variable
+``CUDPP_OPT = -DUSE_CUDPP -Icudpp_mini``.
+
+To support the CUDA multiprocessor server you can set the define
+``-DCUDA_PROXY``.  Please note that in this case you should **not** use
+the CUDA performance primitives and thus set the variable ``CUDPP_OPT``
+to empty.
+
 If the library build is successful, 3 files should be created:
 ``lib/gpu/libgpu.a``\ , ``lib/gpu/nvc_get_devices``\ , and
 ``lib/gpu/Makefile.lammps``\ .  The latter has settings that enable LAMMPS
@ -282,6 +291,7 @@ minutes to hours) to build.  Of course you only need to do that once.)
         -D DOWNLOAD_KIM=value           # download OpenKIM API v2 for build, value = no (default) or yes
         -D LMP_DEBUG_CURL=value         # set libcurl verbose mode on/off, value = off (default) or on
         -D LMP_NO_SSL_CHECK=value       # tell libcurl to not verify the peer, value = no (default) or yes
+         -D KIM_EXTRA_UNITTESTS=value    # enables extra unit tests, value = no (default) or yes

      If ``DOWNLOAD_KIM`` is set to *yes* (or *on*), the KIM API library
      will be downloaded and built inside the CMake build directory.  If
@ -290,6 +300,11 @@ minutes to hours) to build.  Of course you only need to do that once.)
      ``PKG_CONFIG_PATH`` environment variable so that libkim-api can be
      found, or run the command ``source kim-api-activate``.

+      Extra unit tests can only be available if they are explicitly requested
+      (``KIM_EXTRA_UNITTESTS`` is set to *yes* (or *on*)) and the prerequisites
+      are met. See :ref:`KIM Extra unit tests <kim_extra_unittests>` for
+      more details on this.
+
   .. tab:: Traditional make

      You can download and build the KIM library manually if you prefer;
@ -338,6 +353,38 @@ specify your own CA cert path by setting the environment variable
 ``CURL_CA_BUNDLE`` to the path of your choice.  A call to the KIM web
 query would get this value from the environment variable.

+.. _kim_extra_unittests:
+
+KIM Extra unit tests (CMake only)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+During development, testing, or debugging, if
+:doc:`unit testing <Build_development>` is enabled in LAMMPS, one can also
+enable extra tests on :doc:`KIM commands <kim_commands>` by setting the
+``KIM_EXTRA_UNITTESTS`` to *yes* (or *on*).
+
+Enabling the extra unit tests have some requirements,
+
+* It requires to have internet access.
+* It requires to have libcurl installed with the matching development headers
+  and the curl-config tool.
+* It requires to build LAMMPS with the PYTHON package installed and linked to
+  Python 3.6 or later. See the :ref:`PYTHON package build info <python>` for
+  more details on this.
+* It requires to have ``kim-property`` Python package installed, which can be
+  easily done using *pip* as ``pip install kim-property``, or from the
+  *conda-forge* channel as ``conda install kim-property`` if LAMMPS is built in
+  Conda. More detailed information is available at:
+  `kim-property installation <https://github.com/openkim/kim-property#installing-kim-property>`_.
+* It is also necessary to install
+  ``EAM_Dynamo_Mendelev_2007_Zr__MO_848899341753_000``, and
+  ``EAM_Dynamo_ErcolessiAdams_1994_Al__MO_123629422045_005`` KIM models.
+  See `Obtaining KIM Models <http://openkim.org/doc/usage/obtaining-models>`_
+  to learn how to install a pre-build binary of the OpenKIM Repository of
+  Models or see
+  `Installing KIM Models <https://openkim.org/doc/usage/obtaining-models/#installing_models>`_
+  to learn how to install the specific KIM models.
+
 ----------

 .. _kokkos:
@ -482,11 +529,14 @@ They must be specified in uppercase.
   *  - VEGA906
      - GPU
      - AMD GPU MI50/MI60 GFX906
+   *  - VEGA908
+      - GPU
+      - AMD GPU GFX908
   *  - INTEL_GEN
      - GPU
      - Intel GPUs Gen9+

-This list was last updated for version 3.2 of the Kokkos library.
+This list was last updated for version 3.3 of the Kokkos library.

 .. tabs::

@ -732,6 +782,54 @@ be installed on your system.

 ----------

+.. _mliap:
+
+MLIAP package
+---------------------------
+
+Building the MLIAP package requires including the :ref:`SNAP <PKG-SNAP>`
+package.  There will be an error message if this requirement is not satisfied.
+Using the *mliappy* model also requires enabling Python support, which
+in turn requires the :ref:`PYTHON <PKG-PYTHON>`
+package **and** requires you have the `cython <https://cython.org>`_ software
+installed and with it a working ``cythonize`` command.  This feature requires
+compiling LAMMPS with Python version 3.6 or later.
+
+.. tabs::
+
+   .. tab:: CMake build
+
+      .. code-block:: bash
+
+         -D MLIAP_ENABLE_PYTHON=value   # enable mliappy model (default is autodetect)
+
+      Without this setting, CMake will check whether it can find a
+      suitable Python version and the ``cythonize`` command and choose
+      the default accordingly.  During the build procedure the provided
+      .pyx file(s) will be automatically translated to C++ code and compiled.
+      Please do **not** run ``cythonize`` manually in the ``src/MLIAP`` folder,
+      as that can lead to compilation errors if Python support is not enabled.
+      If you did by accident, please remove the generated .cpp and .h files.
+
+   .. tab:: Traditional make
+
+      The build uses the ``lib/python/Makefile.mliap_python`` file in the
+      compile/link process to add a rule to update the files generated by
+      the ``cythonize`` command in case the corresponding .pyx file(s) were
+      modified.  You may need to modify ``lib/python/Makefile.lammps``
+      if the LAMMPS build fails.
+      To manually enforce building MLIAP with Python support enabled,
+      you can add
+      ``-DMLIAP_PYTHON`` to the ``LMP_INC`` variable in your machine makefile.
+      You may have to manually run the ``cythonize`` command on .pyx file(s)
+      in the ``src`` folder, if this is not automatically done during
+      installing the MLIAP package.  Please do **not** run ``cythonize``
+      in the ``src/MLIAP`` folder, as that can lead to compilation errors
+      if Python support is not enabled.
+      If you did by accident, please remove the generated .cpp and .h files.
+
+----------
+
 .. _mscg:

 MSCG package
--- a/doc/src/Build_package.rst
+++ b/doc/src/Build_package.rst
@ -160,6 +160,7 @@ one of them as a starting point and customize it to your needs.
    cmake -C ../cmake/presets/clang.cmake    [OPTIONS] ../cmake  # change settings to use the Clang compilers by default
    cmake -C ../cmake/presets/gcc.cmake      [OPTIONS] ../cmake  # change settings to use the GNU compilers by default
    cmake -C ../cmake/presets/intel.cmake    [OPTIONS] ../cmake  # change settings to use the Intel compilers by default
+    cmake -C ../cmake/presets/pgi.cmake      [OPTIONS] ../cmake  # change settings to use the PGI compilers by default
    cmake -C ../cmake/presets/all_on.cmake   [OPTIONS] ../cmake  # enable all packages
    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
--- a/doc/src/Commands_fix.rst
+++ b/doc/src/Commands_fix.rst
@ -221,6 +221,8 @@ OPT.
   * :doc:`temp/rescale <fix_temp_rescale>`
   * :doc:`temp/rescale/eff <fix_temp_rescale_eff>`
   * :doc:`tfmc <fix_tfmc>`
+   * :doc:`tgnpt/drude <fix_tgnh_drude>`
+   * :doc:`tgnvt/drude <fix_tgnh_drude>`
   * :doc:`thermal/conductivity <fix_thermal_conductivity>`
   * :doc:`ti/spring <fix_ti_spring>`
   * :doc:`tmd <fix_tmd>`
--- a/doc/src/Commands_pair.rst
+++ b/doc/src/Commands_pair.rst
@ -96,6 +96,7 @@ OPT.
   * :doc:`eam/cd <pair_eam>`
   * :doc:`eam/cd/old <pair_eam>`
   * :doc:`eam/fs (gikot) <pair_eam>`
+   * :doc:`eam/he <pair_eam>`
   * :doc:`edip (o) <pair_edip>`
   * :doc:`edip/multi <pair_edip>`
   * :doc:`edpd <pair_mesodpd>`
@ -262,6 +263,7 @@ OPT.
   * :doc:`ufm (got) <pair_ufm>`
   * :doc:`vashishta (gko) <pair_vashishta>`
   * :doc:`vashishta/table (o) <pair_vashishta>`
+   * :doc:`wf/cut <pair_wf_cut>`
   * :doc:`yukawa (gko) <pair_yukawa>`
   * :doc:`yukawa/colloid (go) <pair_yukawa_colloid>`
   * :doc:`zbl (gko) <pair_zbl>`
--- a/doc/src/Commands_parse.rst
+++ b/doc/src/Commands_parse.rst
@ -162,3 +162,26 @@ LAMMPS:
   triple quotes can be nested in the usual manner.  See the doc pages
   for those commands for examples.  Only one of level of nesting is
   allowed, but that should be sufficient for most use cases.
+
+.. admonition:: ASCII versus UTF-8
+        :class: note
+
+   LAMMPS expects and processes 7-bit ASCII format text internally.
+   Many modern environments use UTF-8 encoding, which is a superset
+   of the 7-bit ASCII character table and thus mostly compatible.
+   However, there are several non-ASCII characters that can look
+   very similar to their ASCII equivalents or are invisible (so they
+   look like a blank), but are encoded differently.  Web browsers,
+   PDF viewers, document editors are known to sometimes replace one
+   with the other for a better looking output.  However, that can
+   lead to problems, for instance, when using cut-n-paste of input
+   file examples from web pages, or when using a document editor
+   (not a dedicated plain text editor) for writing LAMMPS inputs.
+   LAMMPS will try to detect this and substitute the non-ASCII
+   characters with their ASCII equivalents where known.  There also
+   is going to be a warning printed, if this occurs.  It is
+   recommended to avoid such characters altogether in LAMMPS input,
+   data and potential files.  The replacement tables are likely
+   incomplete and dependent on users reporting problems processing
+   correctly looking input containing UTF-8 encoded non-ASCII
+   characters.
--- a/doc/src/Developer.rst
+++ b/doc/src/Developer.rst
@ -13,6 +13,7 @@ of time and requests from the LAMMPS user community.
   Developer_org
   Developer_flow
   Developer_write
+   Developer_notes
   Developer_unittest
   Classes
   Developer_utils
--- a/doc/src/Developer_notes.rst
+++ b/doc/src/Developer_notes.rst
@ -0,0 +1,227 @@
+Notes for developers and code maintainers
+-----------------------------------------
+
+This section documents how some of the code functionality within
+LAMMPS works at a conceptual level.  Comments on code in source files
+typically document what a variable stores, what a small section of
+code does, or what a function does and its input/outputs.  The topics
+on this page are intended to document code functionality at a higher level.
+
+Fix contributions to instantaneous energy, virial, and cumulative energy
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Fixes can calculate contributions to the instantaneous energy and/or
+virial of the system, both in a global and peratom sense.  Fixes that
+perform thermostatting or barostatting can calculate the cumulative
+energy they add to or subtract from the system, which is accessed by
+the *ecouple* and *econserve* thermodynamic keywords.  This subsection
+explains how both work and what flags to set in a new fix to enable
+this functionality.
+
+Let's start with thermostatting and barostatting fixes.  Examples are
+the :doc:`fix langevin <fix_langevin>` and :doc:`fix npt <fix_nh>`
+commands.  Here is what the fix needs to do:
+
+* Set the variable *ecouple_flag* = 1 in the constructor.  Also set
+  *scalar_flag* = 1, *extscalar* = 1, and *global_freq* to a timestep
+  increment which matches how often the fix is invoked.
+* Implement a compute_scalar() method that returns the cumulative
+  energy added or subtracted by the fix, e.g. by rescaling the
+  velocity of atoms.  The sign convention is that subtracted energy is
+  positive, added energy is negative.  This must be the total energy
+  added to the entire system, i.e. an "extensive" quantity, not a
+  per-atom energy.  Cumulative means the summed energy since the fix
+  was instantiated, even across multiple runs.  This is because the
+  energy is used by the *econserve* thermodynamic keyword to check
+  that the fix is conserving the total energy of the system,
+  i.e. potential energy + kinetic energy + coupling energy = a
+  constant.
+
+And here is how the code operates:
+
+* The Modify class makes a list of all fixes that set *ecouple_flag* = 1.
+* The :doc:`thermo_style custom <thermo_style>` command defines
+  *ecouple* and *econserve* keywords.
+* These keywords sum the energy contributions from all the
+  *ecouple_flag* = 1 fixes by invoking the energy_couple() method in
+  the Modify class, which calls the compute_scalar() method of each
+  fix in the list.
+
+------------------
+
+Next, here is how a fix contributes to the instantaneous energy and
+virial of the system.  First, it sets any or all of these flags to a
+value of 1 in their constructor:
+
+* *energy_global_flag* to contribute to global energy, example: :doc:`fix indent <fix_indent>`
+* *energy_peratom_flag* to contribute to peratom energy, :doc:`fix cmap <fix_cmap>`
+* *virial_global_flag* to contribute to global virial, example: :doc:`fix wall <fix_wall>`
+* *virial_peratom_flag* to contribute to peratom virial, example: :doc:`fix wall <fix_wall>`
+
+The fix must also do the following:
+
+* For global energy, implement a compute_scalar() method that returns
+  the energy added or subtracted on this timestep.  Here the sign
+  convention is that added energy is positive, subtracted energy is
+  negative.
+* For peratom energy, invoke the ev_init(eflag,vflag) function each
+  time the fix is invoked, which initializes per-atom energy storage.
+  The value of eflag may need to be stored from an earlier call to the
+  fix during the same timestep.  See how the :doc:`fix cmap
+  <fix_cmap>` command does this in src/MOLECULE/fix_cmap.cpp.  When an
+  energy for one or more atoms is calculated, invoke the ev_tally()
+  function to tally the contribution to each atom.  Both the ev_init()
+  and ev_tally() methods are in the parent Fix class.
+* For global and/or peratom virial, invoke the v_init(vflag) function
+  each time the fix is invoked, which initializes virial storage.
+  When forces on one or more atoms are calculated, invoke the
+  v_tally() function to tally the contribution.  Both the v_init() and
+  v_tally() methods are in the parent Fix class.  Note that there are
+  several variants of v_tally(); choose the one appropriate to your
+  fix.
+
+.. note::
+
+   The ev_init() and ev_tally() methods also account for global and
+   peratom virial contributions.  Thus you do not need to invoke the
+   v_init() and v_tally() methods, if the fix also calculates peratom
+   energies.
+
+The fix must also specify whether (by default) to include or exclude
+these contributions to the global/peratom energy/virial of the system.
+For the fix to include the contributions, set either of both of these
+variables in the constructor:
+
+* *thermo_energy* = 1, for global and peratom energy
+* *thermo_virial* = 1, for global and peratom virial
+
+Note that these variables are zeroed in fix.cpp.  Thus if you don't
+set the variables, the contributions will be excluded (by default)
+
+However, the user has ultimate control over whether to include or
+exclude the contributions of the fix via the :doc:`fix modify
+<fix_modify>` command:
+
+* fix modify *energy yes* to include global and peratom energy contributions
+* fix modify *virial yes* to include global and peratom virial contributions
+
+If the fix contributes to any of the global/peratom energy/virial
+values for the system, it should be explained on the fix doc page,
+along with the default values for the *energy yes/no* and *virial
+yes/no* settings of the :doc:`fix modify <fix_modify>` command.
+
+Finally, these 4 contributions are included in the output of 4
+computes:
+
+* global energy in :doc:`compute pe <compute_pe>`
+* peratom energy in :doc:`compute pe/atom <compute_pe_atom>`
+* global virial in :doc:`compute pressure <compute_pressure>`
+* peratom virial in :doc:`compute stress/atom <compute_stress_atom>`
+
+These computes invoke a method of the Modify class to include
+contributions from fixes that have the corresponding flags set,
+e.g. *energy_peratom_flag* and *thermo_energy* for :doc:`compute
+pe/atom <compute_pe_atom>`.
+
+Note that each compute has an optional keyword to either include or
+exclude all contributions from fixes.  Also note that :doc:`compute pe
+<compute_pe>` and :doc:`compute pressure <compute_pressure>` are what
+is used (by default) by :doc:`thermodynamic output <thermo_style>` to
+calculate values for its *pe* and *press* keywords.
+
+KSpace PPPM FFT grids
+^^^^^^^^^^^^^^^^^^^^^
+
+The various :doc:`KSpace PPPM <kspace_style>` styles in LAMMPS use
+FFTs to solve Poisson's equation.  This subsection describes:
+
+* how FFT grids are defined
+* how they are decomposed across processors
+* how they are indexed by each processor
+* how particle charge and electric field values are mapped to/from
+  the grid
+
+An FFT grid cell is a 3d volume; grid points are corners of a grid
+cell and the code stores values assigned to grid points in vectors or
+3d arrays.  A global 3d FFT grid has points indexed 0 to N-1 inclusive
+in each dimension.
+
+Each processor owns two subsets of the grid, each subset is
+brick-shaped.  Depending on how it is used, these subsets are
+allocated as a 1d vector or 3d array.  Either way, the ordering of
+values within contiguous memory x fastest, then y, z slowest.
+
+For the ``3d decomposition`` of the grid, the global grid is
+partitioned into bricks that correspond to the sub-domains of the
+simulation box that each processor owns.  Often, this is a regular 3d
+array (Px by Py by Pz) of bricks, where P = number of processors =
+Px * Py * Pz.  More generally it can be a tiled decomposition, where
+each processor owns a brick and the union of all the bricks is the
+global grid.  Tiled decompositions are produced by load balancing with
+the RCB algorithm; see the :doc:`balance rcb <balance>` command.
+
+For the ``FFT decompostion`` of the grid, each processor owns a brick
+that spans the entire x dimension of the grid while the y and z
+dimensions are partitioned as a regular 2d array (P1 by P2), where P =
+P1 * P2.
+
+The following indices store the inclusive bounds of the brick a
+processor owns, within the global grid:
+
+.. parsed-literal::
+
+   nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in = 3d decomposition brick
+   nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft = FFT decomposition brick
+   nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out = 3d decomposition brick + ghost cells
+
+The ``in`` and ``fft`` indices are from 0 to N-1 inclusive in each
+dimension, where N is the grid size.
+
+The ``out`` indices index an array which stores the ``in`` subset of
+the grid plus ghost cells that surround it.  These indices can thus be
+< 0 or >= N.
+
+The number of ghost cells a processor owns in each of the 6 directions
+is a function of:
+
+.. parsed-literal::
+
+   neighbor skin distance (since atoms can move outside a proc subdomain)
+   qdist = offset or charge from atom due to TIP4P fictitious charge
+   order = mapping stencil size
+   shift = factor used when order is an even number (see below)
+
+Here is an explanation of how the PPPM variables ``order``,
+``nlower`` / ``nupper``, ``shift``, and ``OFFSET`` work. They are the
+relevant variables that determine how atom charge is mapped to grid
+points and how field values are mapped from grid points to atoms:
+
+.. parsed-literal::
+
+   order = # of nearby grid points in each dim that atom charge/field are mapped to/from
+   nlower,nupper = extent of stencil around the grid point an atom is assigned to
+   OFFSET = large integer added/subtracted when mapping to avoid int(-0.75) = 0 when -1 is the desired result
+
+The particle_map() method assigns each atom to a grid point.
+
+If order is even, say 4:
+
+.. parsed-literal::
+
+   atom is assigned to grid point to its left (in each dim)
+   shift = OFFSET
+   nlower = -1, nupper = 2, which are offsets from assigned grid point
+   window of mapping grid pts is thus 2 grid points to left of atom, 2 to right
+
+If order is odd, say 5:
+
+.. parsed-literal::
+
+   atom is assigned to left/right grid pt it is closest to (in each dim)
+   shift = OFFSET + 0.5
+   nlower = 2, nupper = 2
+   if point is in left half of cell, then window of affected grid pts is 3 grid points to left of atom, 2 to right
+   if point is in right half of cell, then window of affected grid pts is 2 grid points to left of atom, 3 to right
+
+These settings apply to each dimension, so that if order = 5, an
+atom's charge is mapped to 125 grid points that surround the atom.
--- a/doc/src/Developer_org.rst
+++ b/doc/src/Developer_org.rst
@ -1,68 +1,75 @@
 Source files
 ------------

-The source files of the LAMMPS code are found in two
-directories of the distribution: ``src`` and ``lib``.
-Most of the code is C++ but there are small numbers of files
-in several other languages.
+The source files of the LAMMPS code are found in two directories of the
+distribution: ``src`` and ``lib``.  Most of the code is written in C++
+but there are small a number of files in several other languages like C,
+Fortran, Shell script, or Python.

-The core of the code is located in the
-``src`` folder and its sub-directories.
-A sizable number of these files are in the ``src`` directory
-itself, but there are plenty of :doc:`packages <Packages>`, which can be
-included or excluded when LAMMPS is built.  See the :doc:`Include
-packages in build <Build_package>` section of the manual for more
-information about that part of the build process.  LAMMPS currently
-supports building with :doc:`conventional makefiles <Build_make>` and
-through :doc:`CMake <Build_cmake>` which differ in how packages are
-enabled or disabled for a LAMMPS binary.  The source files for each
+The core of the code is located in the ``src`` folder and its
+sub-directories.  A sizable number of these files are in the ``src``
+directory itself, but there are plenty of :doc:`packages <Packages>`,
+which can be included or excluded when LAMMPS is built.  See the
+:doc:`Include packages in build <Build_package>` section of the manual
+for more information about that part of the build process.  LAMMPS
+currently supports building with :doc:`conventional makefiles
+<Build_make>` and through :doc:`CMake <Build_cmake>`.  Those procedures
+differ in how packages are enabled or disabled for inclusion into a
+LAMMPS binary so they cannot be mixed.  The source files for each
 package are in all-uppercase sub-directories of the ``src`` folder, for
 example ``src/MOLECULE`` or ``src/USER-MISC``.  The ``src/STUBS``
 sub-directory is not a package but contains a dummy MPI library, that is
 used when building a serial version of the code. The ``src/MAKE``
-directory contains makefiles with settings and flags for a variety of
-configuration and machines for the build process with traditional
-makefiles.
+directory and its sub-directories contain makefiles with settings and
+flags for a variety of configuration and machines for the build process
+with traditional makefiles.

 The ``lib`` directory contains the source code for several supporting
 libraries or files with configuration settings to use globally installed
-libraries, that are required by some of the optional packages.
-Each sub-directory, like ``lib/poems`` or ``lib/gpu``, contains the
-source files, some of which are in different languages such as Fortran
-or CUDA. These libraries are linked to during a LAMMPS build, if the
-corresponding package is installed.
+libraries, that are required by some of the optional packages.  They may
+include python scripts that can transparently download additional source
+code on request.  Each sub-directory, like ``lib/poems`` or ``lib/gpu``,
+contains the source files, some of which are in different languages such
+as Fortran or CUDA. These libraries included in the LAMMPS build,
+if the corresponding package is installed.

 LAMMPS C++ source files almost always come in pairs, such as
 ``src/run.cpp`` (implementation file) and ``src/run.h`` (header file).
-Each pair of files defines a C++
-class, for example the :cpp:class:`LAMMPS_NS::Run` class which contains
-the code invoked by the :doc:`run <run>` command in a LAMMPS input script.
-As this example illustrates, source file and class names often have a
-one-to-one correspondence with a command used in a LAMMPS input script.
-Some source files and classes do not have a corresponding input script
+Each pair of files defines a C++ class, for example the
+:cpp:class:`LAMMPS_NS::Run` class which contains the code invoked by the
+:doc:`run <run>` command in a LAMMPS input script.  As this example
+illustrates, source file and class names often have a one-to-one
+correspondence with a command used in a LAMMPS input script.  Some
+source files and classes do not have a corresponding input script
 command, e.g. ``src/force.cpp`` and the :cpp:class:`LAMMPS_NS::Force`
 class.  They are discussed in the next section.

-A small number of C++ classes and utility functions are implemented with
-only a ``.h`` file. Examples are the Pointer class or the MathVec functions.
+The names of all source files are in lower case and may use the
+underscore character '_' to separate words. Outside of bundled libraries
+which may have different conventions, all C and C++ header files have a
+``.h`` extension, all C++ files have a ``.cpp`` extension, and C files a
+``.c`` extension. A small number of C++ classes and utility functions
+are implemented with only a ``.h`` file. Examples are the Pointer class
+or the MathVec functions.

 Class topology
 --------------

 Though LAMMPS has a lot of source files and classes, its class topology
-is relative flat, as outlined in the :ref:`class-topology` figure.  Each
-name refers to a class and has a pair of associated source files in the
-``src`` folder, for example the class :cpp:class:`LAMMPS_NS::Memory`
-corresponds to the files ``memory.cpp`` and ``memory.h``, or the class
-:cpp:class:`LAMMPS_NS::AtomVec` corresponds to the files
-``atom_vec.cpp`` and ``atom_vec.h``.  Full lines in the figure represent
-compositing: that is the class to the left holds a pointer to an
-instance of the class to the right.  Dashed lines instead represent
-inheritance: the class to the right is derived from the class on the
-left. Classes with a red boundary are not instantiated directly, but
-they represent the base classes for "styles".  Those "styles" make up
-the bulk of the LAMMPS code and only a few typical examples are included
-in the figure for demonstration purposes.
+is not very deep, which can be seen from the :ref:`class-topology`
+figure.  In that figure, each name refers to a class and has a pair of
+associated source files in the ``src`` folder, for example the class
+:cpp:class:`LAMMPS_NS::Memory` corresponds to the files ``memory.cpp``
+and ``memory.h``, or the class :cpp:class:`LAMMPS_NS::AtomVec`
+corresponds to the files ``atom_vec.cpp`` and ``atom_vec.h``.  Full
+lines in the figure represent compositing: that is the class at the base
+of the arrow holds a pointer to an instance of the class at the tip.
+Dashed lines instead represent inheritance: the class to the tip of the
+arrow is derived from the class at the base.  Classes with a red boundary
+are not instantiated directly, but they represent the base classes for
+"styles".  Those "styles" make up the bulk of the LAMMPS code and only
+a few representative examples are included in the figure so it remains
+readable.

 .. _class-topology:
 .. figure:: JPG/lammps-classes.png
@ -82,69 +89,76 @@ in the figure for demonstration purposes.
   derived classes, which may also hold instances of other classes.

 The :cpp:class:`LAMMPS_NS::LAMMPS` class is the topmost class and
-represents what is referred to an "instance" of LAMMPS.  It is a
-composite holding references to instances of other core classes
+represents what is generally referred to an "instance" of LAMMPS.  It is
+a composite holding pointers to instances of other core classes
 providing the core functionality of the MD engine in LAMMPS and through
 them abstractions of the required operations.  The constructor of the
 LAMMPS class will instantiate those instances, process the command line
 flags, initialize MPI (if not already done) and set up file pointers for
-input and output. The destructor will shut everything down and free all
+input and output.  The destructor will shut everything down and free all
 associated memory.  Thus code for the standalone LAMMPS executable in
 ``main.cpp`` simply initializes MPI, instantiates a single instance of
-LAMMPS, and passes it the command line flags and input script. It
+LAMMPS while passing it the command line flags and input script. It
 deletes the LAMMPS instance after the method reading the input returns
 and shuts down the MPI environment before it exits the executable.

 The :cpp:class:`LAMMPS_NS::Pointers` is not shown in the
-:ref:`class-topology` figure, it holds references to members of the
-`LAMMPS_NS::LAMMPS`, so that all classes derived from
-:cpp:class:`LAMMPS_NS::Pointers` have direct access to those reference.
-From the class topology all classes with blue boundary are referenced in
-this class and all classes in the second and third columns, that are not
-listed as derived classes are instead derived from
-:cpp:class:`LAMMPS_NS::Pointers`.
+:ref:`class-topology` figure for clarity.  It holds references to many
+of the members of the `LAMMPS_NS::LAMMPS`, so that all classes derived
+from :cpp:class:`LAMMPS_NS::Pointers` have direct access to those
+reference.  From the class topology all classes with blue boundary are
+referenced in the Pointers class and all classes in the second and third
+columns, that are not listed as derived classes are instead derived from
+:cpp:class:`LAMMPS_NS::Pointers`.  To initialize the pointer references
+in Pointers, a pointer to the LAMMPS class instance needs to be passed
+to the constructor and thus all constructors for classes derived from it
+must do so and pass this pointer to the constructor for Pointers.

-Since all storage is encapsulated, the LAMMPS class can also be
-instantiated multiple times by a calling code, and that can be either
-simultaneously or consecutively.  When running in parallel with MPI,
-care has to be taken, that suitable communicators are used to not
-create conflicts between different instances.
+Since all storage is supposed to be encapsulated (there are a few
+exceptions), the LAMMPS class can also be instantiated multiple times by
+a calling code.  Outside of the aforementioned exceptions, those LAMMPS
+instances can be used alternately.  As of the time of this writing
+(early 2021) LAMMPS is not yet sufficiently thread-safe for concurrent
+execution.  When running in parallel with MPI, care has to be taken,
+that suitable copies of communicators are used to not create conflicts
+between different instances.

-The LAMMPS class currently holds instances of 19 classes representing
-different core functionalities There are a handful of virtual parent
-classes in LAMMPS that define what LAMMPS calls ``styles``.  They are
-shaded red in the :ref:`class-topology` figure.  Each of these are
+The LAMMPS class currently (early 2021) holds instances of 19 classes
+representing the core functionality.  There are a handful of virtual
+parent classes in LAMMPS that define what LAMMPS calls ``styles``.  They
+are shaded red in the :ref:`class-topology` figure.  Each of these are
 parents of a number of child classes that implement the interface
 defined by the parent class.  There are two main categories of these
 ``styles``: some may only have one instance active at a time (e.g. atom,
 pair, bond, angle, dihedral, improper, kspace, comm) and there is a
-dedicated pointer variable in the composite class that manages them.
+dedicated pointer variable for each of them in the composite class.
 Setups that require a mix of different such styles have to use a
-*hybrid* class that manages and forwards calls to the corresponding
-sub-styles for the designated subset of atoms or data. or the composite
-class may have lists of class instances, e.g. Modify handles lists of
-compute and fix styles, while Output handles dumps class instances.
+*hybrid* class that takes the place of the one allowed instance and then
+manages and forwards calls to the corresponding sub-styles for the
+designated subset of atoms or data.  The composite class may also have
+lists of class instances, e.g. Modify handles lists of compute and fix
+styles, while Output handles a list of dump class instances.

-The exception to this scheme are the ``command`` style classes. These
-implement specific commands that can be invoked before, after, or between
-runs or are commands which launch a simulation.  For these an instance
-of the class is created, its command() method called and then, after
-completion, the class instance deleted.  Examples for this are the
-create_box, create_atoms, minimize, run, or velocity command styles.
+The exception to this scheme are the ``command`` style classes.  These
+implement specific commands that can be invoked before, after, or in
+between runs.  For these an instance of the class is created, its
+command() method called and then, after completion, the class instance
+deleted.  Examples for this are the create_box, create_atoms, minimize,
+run, or velocity command styles.

 For all those ``styles`` certain naming conventions are employed: for
-the fix nve command the class is called FixNVE and the files are
+the fix nve command the class is called FixNVE and the source files are
 ``fix_nve.h`` and ``fix_nve.cpp``. Similarly for fix ave/time we have
-FixAveTime and ``fix_ave_time.h`` and ``fix_ave_time.cpp``. Style names
+FixAveTime and ``fix_ave_time.h`` and ``fix_ave_time.cpp``.  Style names
 are lower case and without spaces or special characters. A suffix or
-multiple appended with a forward slash '/' denotes a variant of the
-corresponding class without the suffix. To connect the style name and
-the class name, LAMMPS uses macros like the following ATOM\_CLASS,
-PAIR\_CLASS, BOND\_CLASS, REGION\_CLASS, FIX\_CLASS, COMPUTE\_CLASS,
-or DUMP\_CLASS in the corresponding header file.  During compilation
-files with the pattern ``style_name.h`` are created that contain include
-statements including all headers of all styles of a given type that
-are currently active (or "installed).
+words are appended with a forward slash '/' which denotes a variant of
+the corresponding class without the suffix.  To connect the style name
+and the class name, LAMMPS uses macros like: ``AtomStyle()``,
+``PairStyle()``, ``BondStyle()``, ``RegionStyle()``, and so on in the
+corresponding header file.  During configuration or compilation files
+with the pattern ``style_<name>.h`` are created that consist of a list
+of include statements including all headers of all styles of a given
+type that are currently active (or "installed).


 More details on individual classes in the :ref:`class-topology` are as
@ -152,11 +166,11 @@ follows:

 - The Memory class handles allocation of all large vectors and arrays.

- The Error class prints all error and warning messages.
+- The Error class prints all (terminal) error and warning messages.

- The Universe class sets up partitions of processors so that multiple
-  simulations can be run, each on a subset of the processors allocated
-  for a run, e.g. by the mpirun command.
+- The Universe class sets up one or more partitions of processors so
+  that one or multiple simulations can be run, on the processors
+  allocated for a run, e.g. by the mpirun command.

 - The Input class reads and processes input input strings and files,
  stores variables, and invokes :doc:`commands <Commands_all>`.
@ -241,7 +255,8 @@ follows:
 .. TODO section on "Spatial decomposition and parallel operations"
 ..       diagram of 3d processor grid, brick vs. tiled. local vs. ghost
 ..       atoms, 6-way communication with pack/unpack functions,
-..       PBC as part of the communication
+..       PBC as part of the communication, forward and reverse communication
+..       rendezvous communication, ring communication.

 .. TODO section on "Fixes, Computes, and Variables"
 ..      how and when data is computed and provided and how it is
--- a/doc/src/Developer_utils.rst
+++ b/doc/src/Developer_utils.rst
@ -71,12 +71,21 @@ and parsing files or arguments.

 ----------

+.. doxygenfunction:: strdup
+   :project: progguide
+
 .. doxygenfunction:: trim
   :project: progguide

 .. doxygenfunction:: trim_comment
   :project: progguide

+.. doxygenfunction:: has_utf8
+   :project: progguide
+
+.. doxygenfunction:: utf8_subst
+   :project: progguide
+
 .. doxygenfunction:: count_words(const char *text)
   :project: progguide

@ -286,6 +295,50 @@ This code example should produce the following output:

 ----------

+
+Argument parsing classes
+---------------------------
+
+The purpose of argument parsing classes it to simplify and unify how
+arguments of commands in LAMMPS are parsed and to make abstractions of
+repetitive tasks.
+
+The :cpp:class:`LAMMPS_NS::ArgInfo` class provides an abstraction
+for parsing references to compute or fix styles or variables. These
+would start with a "c\_", "f\_", "v\_" followed by the ID or name of
+than instance and may be postfixed with one or two array indices
+"[<number>]" with numbers > 0.
+
+A typical code segment would look like this:
+
+.. code-block:: C++
+   :caption: Usage example for ArgInfo class
+
+   int nvalues = 0;
+   for (iarg = 0; iarg < nargnew; iarg++) {
+     ArgInfo argi(arg[iarg]);
+
+     which[nvalues] = argi.get_type();
+     argindex[nvalues] = argi.get_index1();
+     ids[nvalues] = argi.copy_name();
+
+     if ((which[nvalues] == ArgInfo::UNKNOWN)
+          || (which[nvalues] == ArgInfo::NONE)
+          || (argi.get_dim() > 1))
+       error->all(FLERR,"Illegal compute XXX command");
+
+     nvalues++;
+   }
+
+----------
+
+.. doxygenclass:: LAMMPS_NS::ArgInfo
+   :project: progguide
+   :members:
+
+
+----------
+
 File reader classes
 -------------------

--- a/doc/src/Howto_drude.rst
+++ b/doc/src/Howto_drude.rst
@ -42,10 +42,11 @@ screening. It may be necessary to use the *extra/special/per/atom*
 keyword of the :doc:`read_data <read_data>` command. If using :doc:`fix shake <fix_shake>`, make sure no Drude particle is in this fix
 group.

-There are two ways to thermostat the Drude particles at a low
+There are three ways to thermostat the Drude particles at a low
 temperature: use either :doc:`fix langevin/drude <fix_langevin_drude>`
 for a Langevin thermostat, or :doc:`fix drude/transform/\* <fix_drude_transform>` for a Nose-Hoover
-thermostat. The former requires use of the command :doc:`comm_modify vel yes <comm_modify>`. The latter requires two separate integration
+thermostat, or :doc:`fix tgnvt/drude <fix_tgnh_drude>` for a temperature-grouped Nose-Hoover thermostat.
+The first and third require use of the command :doc:`comm_modify vel yes <comm_modify>`. The second requires two separate integration
 fixes like *nvt* or *npt*\ . The correct temperatures of the reduced
 degrees of freedom can be calculated using the :doc:`compute temp/drude <compute_temp_drude>`. This requires also to use the
 command *comm_modify vel yes*.
--- a/doc/src/Howto_drude2.rst
+++ b/doc/src/Howto_drude2.rst
@ -221,6 +221,14 @@ modification of forces but no position/velocity updates), the fix

   fix NVE all nve

+To avoid the flying ice cube artifact, where the atoms progressively freeze and the
+center of mass of the whole system drifts faster and faster, the *fix momentum*
+can be used. For instance:
+
+.. code-block:: LAMMPS
+
+   fix MOMENTUM all momentum 100 linear 1 1 1
+
 Finally, do not forget to update the atom type elements if you use
 them in a *dump_modify ... element ...* command, by adding the element
 type of the DPs. Here for instance
@ -376,14 +384,7 @@ For our phenol example, the groups would be defined as

 Note that with the fixes *drude/transform*\ , it is not required to
 specify *comm_modify vel yes* because the fixes do it anyway (several
-times and for the forces also).  To avoid the flying ice cube artifact
-:ref:`(Lamoureux and Roux) <Lamoureux2>`, where the atoms progressively freeze and the
-center of mass of the whole system drifts faster and faster, the *fix
-momentum* can be used. For instance:
-
-.. code-block:: LAMMPS
-
-   fix MOMENTUM all momentum 100 linear 1 1 1
+times and for the forces also).

 It is a bit more tricky to run a NPT simulation with Nose-Hoover
 barostat and thermostat.  First, the volume should be integrated only
@ -404,6 +405,31 @@ instructions for thermostatting and barostatting will look like
   fix NVT DRUDES nvt temp 1. 1. 20
   fix INVERSE all drude/transform/inverse

+Another option for thermalizing the Drude model is to use the
+temperature-grouped Nose-Hoover (TGNH) thermostat proposed by :ref:`(Son) <TGNH-SON>`.
+This is implemented as :doc:`fix tgnvt/drude <fix_tgnh_drude>` and :doc:`fix tgnpt/drude <fix_tgnh_drude>`.
+It separates the kinetic energy into three contributions:
+the molecular center of mass (COM) motion, the motion of atoms or atom-Drude pairs relative to molecular COMs,
+and the relative motion of atom-Drude pairs.
+An independent Nose-Hoover chain is applied to each type of motion.
+When TGNH is used, the temperatures of molecular, atomic and Drude motion can be printed out with :doc:`thermo_style` command.
+
+NVT simulation with TGNH thermostat
+
+.. code-block:: LAMMPS
+
+   comm_modify vel yes
+   fix TGNVT all tgnvt/drude temp 300. 300. 100 1. 20
+   thermo_style custom f_TGNVT[1] f_TGNVT[2] f_TGNVT[3]
+
+NPT simulation with TGNH thermostat
+
+.. code-block:: LAMMPS
+
+   comm_modify vel yes
+   fix TGNPT all tgnpt/drude temp 300. 300. 100 1. 20 iso 1. 1. 500
+   thermo_style custom f_TGNPT[1] f_TGNPT[2] f_TGNPT[3]
+
 ----------

 **Rigid bodies**
@ -480,3 +506,7 @@ NPT ensemble using Nose-Hoover thermostat:

 **(SWM4-NDP)** Lamoureux, Harder, Vorobyov, Roux, MacKerell, Chem Phys
 Let, 418, 245-249 (2006)
+
+.. _TGNH-Son:
+
+**(Son)** Son, McDaniel, Cui and Yethiraj, J Phys Chem Lett, 10, 7523 (2019).
--- a/doc/src/Howto_pylammps.rst
+++ b/doc/src/Howto_pylammps.rst
@ -9,8 +9,8 @@ Overview
 ``PyLammps`` is a Python wrapper class for LAMMPS which can be created
 on its own or use an existing lammps Python object.  It creates a simpler,
 more "pythonic" interface to common LAMMPS functionality, in contrast to
-the ``lammps.py`` wrapper for the C-style LAMMPS library interface which
-is written using `Python ctypes <ctypes_>`_.  The ``lammps.py`` wrapper
+the ``lammps`` wrapper for the C-style LAMMPS library interface which
+is written using `Python ctypes <ctypes_>`_.  The ``lammps`` wrapper
 is discussed on the :doc:`Python_head` doc page.

 Unlike the flat ``ctypes`` interface, PyLammps exposes a discoverable
--- a/doc/src/Intro_citing.rst
+++ b/doc/src/Intro_citing.rst
@ -38,17 +38,18 @@ In addition there are DOIs for individual stable releases. Currently there are:
 Home page
 ^^^^^^^^^

-The LAMMPS website at `https://lammps.sandia.gov/ <https://lammps.sandia.gov>`_ is the canonical
-location for information about LAMMPS and more detailed lists of publications
-using LAMMPS and contributing features.
+The LAMMPS website at `https://lammps.sandia.gov/
+<https://lammps.sandia.gov>`_ is the canonical location for information
+about LAMMPS and its features.

 Citing contributions
 ^^^^^^^^^^^^^^^^^^^^

-LAMMPS has many features and uses previously published methods and
-algorithms or novel features. It also includes potential parameter
-filed for specific models.  You can look up relevant publications either
-in the LAMMPS output to the screen, the ``log.cite`` file (which is
-populated with references to relevant papers through embedding them into
-the source code) and in the documentation of the :doc:`corresponding commands
+LAMMPS has many features and that use either previously published
+methods and algorithms or novel features.  It also includes potential
+parameter filed for specific models.  Where available, a reminder about
+references for optional features used in a specific run is printed to
+the screen and log file.  Style and output location can be selected with
+the :ref:`-cite command-line switch <cite>`.  Additional references are
+given in the documentation of the :doc:`corresponding commands
 <Commands_all>` or in the :doc:`Howto tutorials <Howto>`.
--- a/doc/src/JPG/WF_LJ.jpg
+++ b/doc/src/JPG/WF_LJ.jpg
--- a/doc/src/Library_config.rst
+++ b/doc/src/Library_config.rst
@ -14,6 +14,7 @@ This section documents the following functions:
 - :cpp:func:`lammps_config_has_package`
 - :cpp:func:`lammps_config_package_count`
 - :cpp:func:`lammps_config_package_name`
+- :cpp:func:`lammps_config_accelerator`
 - :cpp:func:`lammps_has_style`
 - :cpp:func:`lammps_style_count`
 - :cpp:func:`lammps_style_name`
@ -126,6 +127,11 @@ approach.

 -----------------------

+.. doxygenfunction:: lammps_config_accelerator
+   :project: progguide
+
+-----------------------
+
 .. doxygenfunction:: lammps_has_style
   :project: progguide

--- a/doc/src/Modify_contribute.rst
+++ b/doc/src/Modify_contribute.rst
@ -206,16 +206,22 @@ packages in the src directory for examples. If you are uncertain, please ask.
  algorithm/science behind the feature itself, or its initial usage, or
  its implementation in LAMMPS), you can add the citation to the \*.cpp
  source file.  See src/USER-EFF/atom_vec_electron.cpp for an example.
-  A LaTeX citation is stored in a variable at the top of the file and a
-  single line of code that references the variable is added to the
-  constructor of the class.  Whenever a user invokes your feature from
-  their input script, this will cause LAMMPS to output the citation to a
-  log.cite file and prompt the user to examine the file.  Note that you
-  should only use this for a paper you or your group authored.
-  E.g. adding a cite in the code for a paper by Nose and Hoover if you
-  write a fix that implements their integrator is not the intended
-  usage.  That kind of citation should just be in the doc page you
-  provide.
+  A LaTeX citation is stored in a variable at the top of the file and
+  a single line of code registering this variable is added to the
+  constructor of the class.  If there is additional functionality (which
+  may have been added later) described in a different publication,
+  additional citation descriptions may be added for as long as they
+  are only registered when the corresponding keyword activating this
+  functionality is used.  With these options it is possible to have
+  LAMMPS output a specific citation reminder whenever a user invokes
+  your feature from their input script.  Note that you should only use
+  this for the most relevant paper for a feature and a publication that
+  you or your group authored.  E.g. adding a citation in the code for
+  a paper by Nose and Hoover if you write a fix that implements their
+  integrator is not the intended usage.  That kind of citation should
+  just be included in the documentation page you provide describing
+  your contribution.  If you are not sure what the best option would
+  be, please contact the LAMMPS developers for advice.

 Finally, as a general rule-of-thumb, the more clear and
 self-explanatory you make your documentation and README files, and the
--- a/doc/src/Packages_details.rst
+++ b/doc/src/Packages_details.rst
@ -662,19 +662,31 @@ MLIAP package

 **Contents:**

-A general interface for machine-learning interatomic potentials.
+A general interface for machine-learning interatomic potentials, including PyTorch.

 **Install:**

-To use this package, also the :ref:`SNAP package <PKG-SNAP>` needs to be installed.
+To use this package, also the :ref:`SNAP package <PKG-SNAP>` package needs
+to be installed.  To make the *mliappy* model available, also the
+:ref:`PYTHON package <PKG-PYTHON>` package needs to be installed, the version
+of Python must be 3.6 or later, and the `cython <https://cython.org/>`_ software
+must be installed.

-**Author:** Aidan Thompson (Sandia).
+**Author:** Aidan Thompson (Sandia), Nicholas Lubbers (LANL).

 **Supporting info:**

 * src/MLIAP: filenames -> commands
+* src/MLIAP/README
 * :doc:`pair_style mliap <pair_mliap>`
-* examples/mliap
+* :doc:`compute_style mliap <compute_mliap>`
+* examples/mliap (see README)
+
+When built with the *mliappy* model this package includes an extension for
+coupling with Python models, including PyTorch. In this case, the Python
+interpreter linked to LAMMPS will need the ``cython`` and ``numpy`` modules
+installed.  The provided examples build models with PyTorch, which would
+therefore also needs to be installed to run those examples.

 ----------

--- a/doc/src/Python_config.rst
+++ b/doc/src/Python_config.rst
@ -39,7 +39,9 @@ about compile time settings and included packages and styles.
 * :py:attr:`lammps.has_jpeg_support <lammps.lammps.has_jpeg_support>`
 * :py:attr:`lammps.has_ffmpeg_support <lammps.lammps.has_ffmpeg_support>`

-* :py:attr:`lammps.installed_packages <lammps.lammps.installed_pages>`
+* :py:attr:`lammps.installed_packages <lammps.lammps.installed_packages>`
+
+* :py:meth:`lammps.get_accelerator_config <lammps.lammps.accelerator_config>`

 * :py:meth:`lammps.has_style() <lammps.lammps.has_style()>`
 * :py:meth:`lammps.available_styles() <lammps.lammps.available_styles()>`
--- a/doc/src/Python_ext.rst
+++ b/doc/src/Python_ext.rst
@ -9,7 +9,7 @@ This means you can extend the Python wrapper by following these steps:
 * Add a new interface function to ``src/library.cpp`` and
  ``src/library.h``.
 * Rebuild LAMMPS as a shared library.
-* Add a wrapper method to ``python/lammps.py`` for this interface
+* Add a wrapper method to ``python/lammps/core.py`` for this interface
  function.
 * Define the corresponding ``argtypes`` list and ``restype``
  in the ``lammps.__init__()`` function.
--- a/doc/src/Python_install.rst
+++ b/doc/src/Python_install.rst
@ -8,9 +8,9 @@ module.  Because of the dynamic loading, it is required that LAMMPS is
 compiled in :ref:`"shared" mode <exe>`.  It is also recommended to
 compile LAMMPS with :ref:`C++ exceptions <exceptions>` enabled.

-Two files are necessary for Python to be able to invoke LAMMPS code:
+Two components are necessary for Python to be able to invoke LAMMPS code:

-* The LAMMPS Python Module (``lammps.py``) from the ``python`` folder
+* The LAMMPS Python Package (``lammps``) from the ``python`` folder
 * The LAMMPS Shared Library (``liblammps.so``, ``liblammps.dylib`` or
  ``liblammps.dll``) from the folder where you compiled LAMMPS.

@ -25,10 +25,10 @@ Installing the LAMMPS Python Module and Shared Library
 ======================================================

 Making LAMMPS usable within Python and vice versa requires putting the
-LAMMPS Python module file (``lammps.py``) into a location where the
+LAMMPS Python package (``lammps``) into a location where the
 Python interpreter can find it and installing the LAMMPS shared library
-into a folder that the dynamic loader searches or into the same folder
-where the ``lammps.py`` file is.  There are multiple ways to achieve
+into a folder that the dynamic loader searches or inside of the installed
+``lammps`` package folder.  There are multiple ways to achieve
 this.

 #. Do a full LAMMPS installation of libraries, executables, selected
@ -36,13 +36,13 @@ this.
   available via CMake), which can also be either system-wide or into
   user specific folders.

-#. Install both files into a Python ``site-packages`` folder, either
+#. Install both components into a Python ``site-packages`` folder, either
   system-wide or in the corresponding user-specific folder. This way no
   additional environment variables need to be set, but the shared
   library is otherwise not accessible.

-#. Do an installation into a virtual environment. This can either be
-   an installation of the python module only or a full installation.
+#. Do an installation into a virtual environment. This can either be an
+   installation of the Python package only or a full installation of LAMMPS.

 #. Leave the files where they are in the source/development tree and
   adjust some environment variables.
@ -69,7 +69,7 @@ this.
         cd build

         # configure LAMMPS compilation
-         cmake -C cmake/presets/minimal.cmake -D BUILD_SHARED_LIBS=on \
+         cmake -C ../cmake/presets/minimal.cmake -D BUILD_SHARED_LIBS=on \
               -D LAMMPS_EXCEPTIONS=on -D PKG_PYTHON=on ../cmake

         # compile LAMMPS
@ -81,40 +81,42 @@ this.

      This leads to an installation to the following locations:

-      +------------------------+-----------------------------------------------------------+-------------------------------------------------------------+
-      | File                   | Location                                                  | Notes                                                       |
-      +========================+===========================================================+=============================================================+
-      | LAMMPS Python Module   | * ``$HOME/.local/lib/pythonX.Y/site-packages/`` (32bit)   | ``X.Y`` depends on the installed Python version             |
-      |                        | * ``$HOME/.local/lib64/pythonX.Y/site-packages/`` (64bit) |                                                             |
-      +------------------------+-----------------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS shared library  | * ``$HOME/.local/lib/`` (32bit)                           |                                                             |
-      |                        | * ``$HOME/.local/lib64/`` (64bit)                         |                                                             |
-      +------------------------+-----------------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS executable      | * ``$HOME/.local/bin/``                                   |                                                             |
-      +------------------------+-----------------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS potential files | * ``$HOME/.local/share/lammps/potentials/``               | Set ``LAMMPS_POTENTIALS`` environment variable to this path |
-      +------------------------+-----------------------------------------------------------+-------------------------------------------------------------+
+      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+
+      | File                   | Location                                                        | Notes                                                       |
+      +========================+=================================================================+=============================================================+
+      | LAMMPS Python package  | * ``$HOME/.local/lib/pythonX.Y/site-packages/lammps`` (32bit)   | ``X.Y`` depends on the installed Python version             |
+      |                        | * ``$HOME/.local/lib64/pythonX.Y/site-packages/lammps`` (64bit) |                                                             |
+      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+
+      | LAMMPS shared library  | * ``$HOME/.local/lib/`` (32bit)                                 | Set shared loader environment variable to this path         |
+      |                        | * ``$HOME/.local/lib64/`` (64bit)                               | (see below for more info on this)                           |
+      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+
+      | LAMMPS executable      | * ``$HOME/.local/bin/``                                         |                                                             |
+      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+
+      | LAMMPS potential files | * ``$HOME/.local/share/lammps/potentials/``                     | Set ``LAMMPS_POTENTIALS`` environment variable to this path |
+      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+

      For a system-wide installation you need to set
      ``CMAKE_INSTALL_PREFIX`` to a system folder like ``/usr`` (or
-      ``/usr/local``).  The installation step (**not** the
+      ``/usr/local``); the default is ``${HOME}/.local``.  The
+      installation step for a system folder installation (**not** the
      configuration/compilation) needs to be done with superuser
      privilege, e.g. by using ``sudo cmake --install .``.  The
-      installation folders will then by changed to:
+      installation folders will then be changed to (assuming ``/usr`` as
+      prefix):

-      +------------------------+---------------------------------------------------+-------------------------------------------------------------+
-      | File                   | Location                                          | Notes                                                       |
-      +========================+===================================================+=============================================================+
-      | LAMMPS Python Module   | * ``/usr/lib/pythonX.Y/site-packages/`` (32bit)   | ``X.Y`` depends on the installed Python version             |
-      |                        | * ``/usr/lib64/pythonX.Y/site-packages/`` (64bit) |                                                             |
-      +------------------------+---------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS shared library  | * ``/usr/lib/`` (32bit)                           |                                                             |
-      |                        | * ``/usr/lib64/`` (64bit)                         |                                                             |
-      +------------------------+---------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS executable      | * ``/usr/bin/``                                   |                                                             |
-      +------------------------+---------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS potential files | * ``/usr/share/lammps/potentials/``               |                                                             |
-      +------------------------+---------------------------------------------------+-------------------------------------------------------------+
+      +------------------------+---------------------------------------------------------+-------------------------------------------------------------+
+      | File                   | Location                                                | Notes                                                       |
+      +========================+=========================================================+=============================================================+
+      | LAMMPS Python package  | * ``/usr/lib/pythonX.Y/site-packages/lammps`` (32bit)   | ``X.Y`` depends on the installed Python version             |
+      |                        | * ``/usr/lib64/pythonX.Y/site-packages/lammps`` (64bit) |                                                             |
+      +------------------------+---------------------------------------------------------+-------------------------------------------------------------+
+      | LAMMPS shared library  | * ``/usr/lib/`` (32bit)                                 |                                                             |
+      |                        | * ``/usr/lib64/`` (64bit)                               |                                                             |
+      +------------------------+---------------------------------------------------------+-------------------------------------------------------------+
+      | LAMMPS executable      | * ``/usr/bin/``                                         |                                                             |
+      +------------------------+---------------------------------------------------------+-------------------------------------------------------------+
+      | LAMMPS potential files | * ``/usr/share/lammps/potentials/``                     |                                                             |
+      +------------------------+---------------------------------------------------------+-------------------------------------------------------------+

      To be able to use the "user" installation you have to ensure that
      the folder containing the LAMMPS shared library is either included
@ -146,7 +148,7 @@ this.
      necessary due to files installed in system folders that are loaded
      automatically when a login shell is started.

-   .. tab:: Python module only
+   .. tab:: Python package only

      Compile LAMMPS with either :doc:`CMake <Build_cmake>` or the
      :doc:`traditional make <Build_make>` procedure in :ref:`shared
@ -157,37 +159,37 @@ this.

         make install-python

-      This will try to install (only) the shared library and the python
-      module into a system folder and if that fails (due to missing
+      This will try to install (only) the shared library and the Python
+      package into a system folder and if that fails (due to missing
      write permissions) will instead do the installation to a user
      folder under ``$HOME/.local``.  For a system-wide installation you
      would have to gain superuser privilege, e.g. though ``sudo``

-      +------------------------+-----------------------------------------------------------+-------------------------------------------------------------+
-      | File                   | Location                                                  | Notes                                                       |
-      +========================+===========================================================+=============================================================+
-      | LAMMPS Python Module   | * ``$HOME/.local/lib/pythonX.Y/site-packages/`` (32bit)   | ``X.Y`` depends on the installed Python version             |
-      |                        | * ``$HOME/.local/lib64/pythonX.Y/site-packages/`` (64bit) |                                                             |
-      +------------------------+-----------------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS shared library  | * ``$HOME/.local/lib/pythonX.Y/site-packages/`` (32bit)   | ``X.Y`` depends on the installed Python version             |
-      |                        | * ``$HOME/.local/lib64/pythonX.Y/site-packages/`` (64bit) |                                                             |
-      +------------------------+-----------------------------------------------------------+-------------------------------------------------------------+
+      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+
+      | File                   | Location                                                        | Notes                                                       |
+      +========================+=================================================================+=============================================================+
+      | LAMMPS Python package  | * ``$HOME/.local/lib/pythonX.Y/site-packages/lammps`` (32bit)   | ``X.Y`` depends on the installed Python version             |
+      |                        | * ``$HOME/.local/lib64/pythonX.Y/site-packages/lammps`` (64bit) |                                                             |
+      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+
+      | LAMMPS shared library  | * ``$HOME/.local/lib/pythonX.Y/site-packages/lammps`` (32bit)   | ``X.Y`` depends on the installed Python version             |
+      |                        | * ``$HOME/.local/lib64/pythonX.Y/site-packages/lammps`` (64bit) |                                                             |
+      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+

      For a system-wide installation those folders would then become.

-      +------------------------+---------------------------------------------------+-------------------------------------------------------------+
-      | File                   | Location                                          | Notes                                                       |
-      +========================+===================================================+=============================================================+
-      | LAMMPS Python Module   | * ``/usr/lib/pythonX.Y/site-packages/`` (32bit)   | ``X.Y`` depends on the installed Python version             |
-      |                        | * ``/usr/lib64/pythonX.Y/site-packages/`` (64bit) |                                                             |
-      +------------------------+---------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS shared library  | * ``/usr/lib/pythonX.Y/site-packages/`` (32bit)   | ``X.Y`` depends on the installed Python version             |
-      |                        | * ``/usr/lib64/pythonX.Y/site-packages/`` (64bit) |                                                             |
-      +------------------------+---------------------------------------------------+-------------------------------------------------------------+
+      +------------------------+---------------------------------------------------------+-------------------------------------------------------------+
+      | File                   | Location                                                | Notes                                                       |
+      +========================+=========================================================+=============================================================+
+      | LAMMPS Python package  | * ``/usr/lib/pythonX.Y/site-packages/lammps`` (32bit)   | ``X.Y`` depends on the installed Python version             |
+      |                        | * ``/usr/lib64/pythonX.Y/site-packages/lammps`` (64bit) |                                                             |
+      +------------------------+---------------------------------------------------------+-------------------------------------------------------------+
+      | LAMMPS shared library  | * ``/usr/lib/pythonX.Y/site-packages/lammps`` (32bit)   | ``X.Y`` depends on the installed Python version             |
+      |                        | * ``/usr/lib64/pythonX.Y/site-packages/lammps`` (64bit) |                                                             |
+      +------------------------+---------------------------------------------------------+-------------------------------------------------------------+

      No environment variables need to be set for those, as those
      folders are searched by default by Python or the LAMMPS Python
-      module.
+      package.

      For the traditional make process you can override the python
      version to version x.y when calling ``make`` with
@ -199,9 +201,9 @@ this.

      .. code-block:: bash

-         $ python install.py -m <python module> -l <shared library> -v <version.h file> [-d <pydir>]
+         $ python install.py -p <python package> -l <shared library> -v <version.h file> [-d <pydir>]

-      * The ``-m`` flag points to the ``lammps.py`` python module file to be installed,
+      * The ``-p`` flag points to the ``lammps`` Python package folder to be installed,
      * the ``-l`` flag points to the LAMMPS shared library file to be installed,
      * the ``-v`` flag points to the ``version.h`` file in the LAMMPS source
      * and the optional ``-d`` flag to a custom (legacy) installation folder
@ -249,38 +251,38 @@ this.
      When using CMake to build LAMMPS, you need to set
      ``CMAKE_INSTALL_PREFIX`` to the value of the ``$VIRTUAL_ENV``
      environment variable during the configuration step. For the
-      traditional make procedure, not additional steps are needed.
-      After compiling LAMMPS you can do a "Python module only"
+      traditional make procedure, no additional steps are needed.
+      After compiling LAMMPS you can do a "Python package only"
      installation with ``make install-python`` and the LAMMPS Python
-      module and the shared library file are installed into the
+      package and the shared library file are installed into the
      following locations:

-      +------------------------+-----------------------------------------------------------+-------------------------------------------------------------+
-      | File                   | Location                                                  | Notes                                                       |
-      +========================+===========================================================+=============================================================+
-      | LAMMPS Python Module   | * ``$VIRTUAL_ENV/lib/pythonX.Y/site-packages/`` (32bit)   | ``X.Y`` depends on the installed Python version             |
-      |                        | * ``$VIRTUAL_ENV/lib64/pythonX.Y/site-packages/`` (64bit) |                                                             |
-      +------------------------+-----------------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS shared library  | * ``$VIRTUAL_ENV/lib/pythonX.Y/site-packages/`` (32bit)   | ``X.Y`` depends on the installed Python version             |
-      |                        | * ``$VIRTUAL_ENV/lib64/pythonX.Y/site-packages/`` (64bit) |                                                             |
-      +------------------------+-----------------------------------------------------------+-------------------------------------------------------------+
+      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+
+      | File                   | Location                                                        | Notes                                                       |
+      +========================+=================================================================+=============================================================+
+      | LAMMPS Python Module   | * ``$VIRTUAL_ENV/lib/pythonX.Y/site-packages/lammps`` (32bit)   | ``X.Y`` depends on the installed Python version             |
+      |                        | * ``$VIRTUAL_ENV/lib64/pythonX.Y/site-packages/lammps`` (64bit) |                                                             |
+      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+
+      | LAMMPS shared library  | * ``$VIRTUAL_ENV/lib/pythonX.Y/site-packages/lammps`` (32bit)   | ``X.Y`` depends on the installed Python version             |
+      |                        | * ``$VIRTUAL_ENV/lib64/pythonX.Y/site-packages/lammps`` (64bit) |                                                             |
+      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+

      If you do a full installation (CMake only) with "install", this
      leads to the following installation locations:

-      +------------------------+-----------------------------------------------------------+-------------------------------------------------------------+
-      | File                   | Location                                                  | Notes                                                       |
-      +========================+===========================================================+=============================================================+
-      | LAMMPS Python Module   | * ``$VIRTUAL_ENV/lib/pythonX.Y/site-packages/`` (32bit)   | ``X.Y`` depends on the installed Python version             |
-      |                        | * ``$VIRTUAL_ENV/lib64/pythonX.Y/site-packages/`` (64bit) |                                                             |
-      +------------------------+-----------------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS shared library  | * ``$VIRTUAL_ENV/lib/`` (32bit)                           |                                                             |
-      |                        | * ``$VIRTUAL_ENV/lib64/`` (64bit)                         |                                                             |
-      +------------------------+-----------------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS executable      | * ``$VIRTUAL_ENV/bin/``                                   |                                                             |
-      +------------------------+-----------------------------------------------------------+-------------------------------------------------------------+
-      | LAMMPS potential files | * ``$VIRTUAL_ENV/share/lammps/potentials/``               |                                                             |
-      +------------------------+-----------------------------------------------------------+-------------------------------------------------------------+
+      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+
+      | File                   | Location                                                        | Notes                                                       |
+      +========================+=================================================================+=============================================================+
+      | LAMMPS Python Module   | * ``$VIRTUAL_ENV/lib/pythonX.Y/site-packages/lammps`` (32bit)   | ``X.Y`` depends on the installed Python version             |
+      |                        | * ``$VIRTUAL_ENV/lib64/pythonX.Y/site-packages/lammps`` (64bit) |                                                             |
+      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+
+      | LAMMPS shared library  | * ``$VIRTUAL_ENV/lib/`` (32bit)                                 | Set shared loader environment variable to this path         |
+      |                        | * ``$VIRTUAL_ENV/lib64/`` (64bit)                               | (see below for more info on this)                           |
+      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+
+      | LAMMPS executable      | * ``$VIRTUAL_ENV/bin/``                                         |                                                             |
+      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+
+      | LAMMPS potential files | * ``$VIRTUAL_ENV/share/lammps/potentials/``                     | Set ``LAMMPS_POTENTIALS`` environment variable to this path |
+      +------------------------+-----------------------------------------------------------------+-------------------------------------------------------------+

      In that case you need to modify the ``$HOME/myenv/bin/activate``
      script in a similar fashion you need to update your
@ -296,15 +298,15 @@ this.
         echo 'export LD_LIBRARY_PATH=$VIRTUAL_ENV/lib:$LD_LIBRARY_PATH' >> $HOME/myenv/bin/activate

         # MacOS
-         echo 'export DYLD_LIBRARY_PATH=$VIRTUAL_ENV/lib:$LD_LIBRARY_PATH' >> $HOME/myenv/bin/activate
+         echo 'export DYLD_LIBRARY_PATH=$VIRTUAL_ENV/lib:$DYLD_LIBRARY_PATH' >> $HOME/myenv/bin/activate

   .. tab:: In place usage

      You can also :doc:`compile LAMMPS <Build>` as usual in
      :ref:`"shared" mode <exe>` leave the shared library and Python
-      module files inside the source/compilation folders. Instead of
+      package inside the source/compilation folders. Instead of
      copying the files where they can be found, you need to set the environment
-      variables ``PYTHONPATH`` (for the Python module) and
+      variables ``PYTHONPATH`` (for the Python package) and
      ``LD_LIBRARY_PATH`` (or ``DYLD_LIBRARY_PATH`` on MacOS

      For Bourne shells (bash, ksh and similar) the commands are:
@ -325,6 +327,10 @@ this.
      You can make those changes permanent by editing your ``$HOME/.bashrc``
      or ``$HOME/.login`` files, respectively.

+      .. note::
+
+         The ``PYTHONPATH`` needs to point to the parent folder that contains the ``lammps`` package!
+

 To verify if LAMMPS can be successfully started from Python, start the
 Python interpreter, load the ``lammps`` Python module and create a
@ -346,7 +352,7 @@ output similar to the following:
 .. note::

   Unless you opted for "In place use", you will have to rerun the installation
-   any time you recompile LAMMPS to ensure the latest Python module and shared
+   any time you recompile LAMMPS to ensure the latest Python package and shared
   library are installed and used.

 .. note::
--- a/doc/src/Python_module.rst
+++ b/doc/src/Python_module.rst
@ -3,12 +3,12 @@ The ``lammps`` Python module

 .. py:module:: lammps

-The LAMMPS Python interface is implemented as a module called
-:py:mod:`lammps` in the ``lammps.py`` file in the ``python`` folder of
-the LAMMPS source code distribution.  After compilation of LAMMPS, the
-module can be installed into a Python system folder or a user folder
-with ``make install-python``.  Components of the module can then loaded
-into a Python session with the ``import`` command.
+The LAMMPS Python interface is implemented as a module called :py:mod:`lammps`
+which is defined in the ``lammps`` package in the ``python`` folder of the
+LAMMPS source code distribution.  After compilation of LAMMPS, the module can
+be installed into a Python system folder or a user folder with ``make
+install-python``.  Components of the module can then loaded into a Python
+session with the ``import`` command.

 There are multiple Python interface classes in the :py:mod:`lammps` module:

@ -26,6 +26,23 @@ There are multiple Python interface classes in the :py:mod:`lammps` module:

 .. _mpi4py_url: https://mpi4py.readthedocs.io

+.. admonition:: Version check
+   :class: note
+
+   The :py:mod:`lammps` module stores the version number of the LAMMPS
+   version it is installed from.  When initializing the
+   :py:class:`lammps <lammps.lammps>` class, this version is checked to
+   be the same as the result from :py:func:`lammps.version`, the version
+   of the LAMMPS shared library that the module interfaces to.  If the
+   they are not the same an AttributeError exception is raised since a
+   mismatch of versions (e.g.  due to incorrect use of the
+   ``LD_LIBRARY_PATH`` or ``PYTHONPATH`` environment variables can lead
+   to crashes or data corruption and otherwise incorrect behavior.
+
+.. automodule:: lammps
+   :members:
+   :noindex:
+
 ----------

 The ``lammps`` class API
@ -44,7 +61,7 @@ functions. Below is a detailed documentation of the API.
 .. autoclass:: lammps.lammps
   :members:

-.. autoclass:: lammps.numpy_wrapper
+.. autoclass:: lammps.numpy::numpy_wrapper
   :members:

 ----------
@ -117,8 +134,8 @@ Style Constants
   to request from computes or fixes. See :cpp:enum:`_LMP_STYLE_CONST`
   for the equivalent constants in the C library interface. Used in
   :py:func:`lammps.extract_compute`, :py:func:`lammps.extract_fix`, and their NumPy variants
-   :py:func:`lammps.numpy.extract_compute() <numpy_wrapper.extract_compute>` and
-   :py:func:`lammps.numpy.extract_fix() <numpy_wrapper.extract_fix>`.
+   :py:func:`lammps.numpy.extract_compute() <lammps.numpy.numpy_wrapper.extract_compute>` and
+   :py:func:`lammps.numpy.extract_fix() <lammps.numpy.numpy_wrapper.extract_fix>`.

 .. _py_type_constants:

@ -132,8 +149,8 @@ Type Constants
   to request  from computes  or fixes.  See :cpp:enum:`_LMP_TYPE_CONST`
   for the equivalent constants in the C library interface. Used in
   :py:func:`lammps.extract_compute`, :py:func:`lammps.extract_fix`, and their NumPy variants
-   :py:func:`lammps.numpy.extract_compute() <numpy_wrapper.extract_compute>` and
-   :py:func:`lammps.numpy.extract_fix() <numpy_wrapper.extract_fix>`.
+   :py:func:`lammps.numpy.extract_compute() <lammps.numpy.numpy_wrapper.extract_compute>` and
+   :py:func:`lammps.numpy.extract_fix() <lammps.numpy.numpy_wrapper.extract_fix>`.

 .. _py_vartype_constants:

@ -153,6 +170,6 @@ Classes representing internal objects
   :members:
   :no-undoc-members:

-.. autoclass:: lammps.NumPyNeighList
+.. autoclass:: lammps.numpy::NumPyNeighList
   :members:
   :no-undoc-members:
--- a/doc/src/Python_overview.rst
+++ b/doc/src/Python_overview.rst
@ -2,9 +2,9 @@ Overview
 ========

 The LAMMPS distribution includes a ``python`` directory with the Python
-code needed to run LAMMPS from Python.  The ``python/lammps.py``
-contains :doc:`the "lammps" Python <Python_module>` that wraps the
-LAMMPS C-library interface.  This file makes it is possible to do the
+code needed to run LAMMPS from Python.  The ``python/lammps`` package
+contains :doc:`the "lammps" Python module <Python_module>` that wraps the
+LAMMPS C-library interface.  This module makes it is possible to do the
 following either from a Python script, or interactively from a Python
 prompt:

@ -20,8 +20,8 @@ have a version of Python that extends Python to enable multiple
 instances of Python to read what you type.

 To do all of this, you must build LAMMPS in :ref:`"shared" mode <exe>`
-and make certain that your Python interpreter can find the ``lammps.py``
-file and the LAMMPS shared library file.
+and make certain that your Python interpreter can find the ``lammps``
+Python package and the LAMMPS shared library file.

 .. _ctypes: https://docs.python.org/3/library/ctypes.html

--- a/doc/src/Python_trouble.rst
+++ b/doc/src/Python_trouble.rst
@ -33,7 +33,7 @@ the constructor call as follows (see :ref:`python_create_lammps` for more detail
   >>> lmp = lammps(name='mpi')

 You can also test the load directly in Python as follows, without
-first importing from the lammps.py file:
+first importing from the ``lammps`` module:

 .. code-block:: python

--- a/doc/src/Run_options.rst
+++ b/doc/src/Run_options.rst
@ -11,6 +11,7 @@ letter abbreviation can be used:
 * :ref:`-k or -kokkos <run-kokkos>`
 * :ref:`-l or -log <log>`
 * :ref:`-m or -mpicolor <mpicolor>`
+* :ref:`-c or -cite <cite>`
 * :ref:`-nc or -nocite <nocite>`
 * :ref:`-pk or -package <package>`
 * :ref:`-p or -partition <partition>`
@ -220,14 +221,31 @@ links with from the lib/message directory.  See the

 ----------

+.. _cite:
+
+**-cite style or file name**
+
+Select how and where to output a reminder about citing contributions
+to the LAMMPS code that were used during the run. Available styles are
+"both", "none", "screen", or "log".  Any flag will be considered a file
+name to write the detailed citation info to.  Default is the "log" style
+where there is a short summary in the screen output and detailed citations
+in BibTeX format in the logfile.  The option "both" selects the detailed
+output for both, "none", the short output for both, and "screen" will
+write the detailed info to the screen and the short version to the log
+file.  If a dedicated citation info file is requested, the screen and
+log file output will be in the short format (same as with "none").
+
+See the :doc:`citation page <Intro_citing>` for more details on
+how to correctly reference and cite LAMMPS.
+
+----------
+
 .. _nocite:

 **-nocite**

-Disable writing the log.cite file which is normally written to list
-references for specific cite-able features used during a LAMMPS run.
-See the `citation page <https://lammps.sandia.gov/cite.html>`_ for more
-details.
+Disable generating a citation reminder (see above) at all.

 ----------

--- a/doc/src/Speed_kokkos.rst
+++ b/doc/src/Speed_kokkos.rst
@ -26,6 +26,15 @@ task). These are Serial (MPI-only for CPUs and Intel Phi), OpenMP
 GPUs) and HIP (for AMD GPUs). You choose the mode at build time to
 produce an executable compatible with a specific hardware.

+.. admonition:: C++14 support
+      :class: note
+
+   Kokkos requires using a compiler that supports the c++14 standard. For
+   some compilers, it may be necessary to add a flag to enable c++14 support.
+   For example, the GNU compiler uses the -std=c++14 flag. For a list of
+   compilers that have been tested with the Kokkos library, see the Kokkos
+   `README <https://github.com/kokkos/kokkos/blob/master/README.md>`_.
+
 .. admonition:: NVIDIA CUDA support
   :class: note

@ -38,14 +47,14 @@ produce an executable compatible with a specific hardware.
   :class: note

   Kokkos with CUDA currently implicitly assumes that the MPI library is
-   CUDA-aware. This is not always the case, especially when using
+   GPU-aware. This is not always the case, especially when using
   pre-compiled MPI libraries provided by a Linux distribution. This is
   not a problem when using only a single GPU with a single MPI
   rank. When running with multiple MPI ranks, you may see segmentation
-   faults without CUDA-aware MPI support. These can be avoided by adding
-   the flags :doc:`-pk kokkos cuda/aware off <Run_options>` to the
+   faults without GPU-aware MPI support. These can be avoided by adding
+   the flags :doc:`-pk kokkos gpu/aware off <Run_options>` to the
   LAMMPS command line or by using the command :doc:`package kokkos
-   cuda/aware off <package>` in the input file.
+   gpu/aware off <package>` in the input file.

 .. admonition:: AMD GPU support
   :class: note
@ -242,8 +251,8 @@ case, also packing/unpacking communication buffers on the host may give
 speedup (see the KOKKOS :doc:`package <package>` command). Using CUDA MPS
 is recommended in this scenario.

-Using a CUDA-aware MPI library is highly recommended. CUDA-aware MPI use can be
-avoided by using :doc:`-pk kokkos cuda/aware no <package>`. As above for
+Using a GPU-aware MPI library is highly recommended. GPU-aware MPI use can be
+avoided by using :doc:`-pk kokkos gpu/aware off <package>`. As above for
 multi-core CPUs (and no GPU), if N is the number of physical cores/node,
 then the number of MPI tasks/node should not exceed N.

--- a/doc/src/balance.rst
+++ b/doc/src/balance.rst
@ -288,7 +288,7 @@ adjacent planes are closer together than the neighbor skin distance
 (as specified by the :doc`neigh_modify <neigh_modify>` command), then
 the plane positions are shifted to separate them by at least this
 amount.  This is to prevent particles being lost when dynamics are run
-with processor subdomains that are too narrow in one or more
+with processor sub-domains that are too narrow in one or more
 dimensions.

 Once the re-balancing is complete and final processor sub-domains
--- a/doc/src/comm_modify.rst
+++ b/doc/src/comm_modify.rst
@ -84,13 +84,15 @@ information is available, then also a heuristic based on that bond length
 is computed. It is used as communication cutoff, if there is no pair
 style present and no *comm_modify cutoff* command used. Otherwise a
 warning is printed, if this bond based estimate is larger than the
-communication cutoff used. A
+communication cutoff used.

 The *cutoff/multi* option is equivalent to *cutoff*\ , but applies to
 communication mode *multi* instead. Since in this case the communication
 cutoffs are determined per atom type, a type specifier is needed and
 cutoff for one or multiple types can be extended. Also ranges of types
-using the usual asterisk notation can be given.
+using the usual asterisk notation can be given. For granular pair styles,
+the default cutoff is set to the sum of the current maximum atomic radii
+for each type.

 These are simulation scenarios in which it may be useful or even
 necessary to set a ghost cutoff > neighbor cutoff:
--- a/doc/src/compute_displace_atom.rst
+++ b/doc/src/compute_displace_atom.rst
@ -46,11 +46,12 @@ the compute command was issued.  The value of the displacement will be
 .. 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
-   using the :doc:`set image <set>` command.
+   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::

--- a/doc/src/compute_mliap.rst
+++ b/doc/src/compute_mliap.rst
@ -18,7 +18,7 @@ Syntax
  .. parsed-literal::

       *model* values = style
-         style = *linear* or *quadratic*
+         style = *linear* or *quadratic* or *mliappy*
       *descriptor* values = style filename
         style = *sna*
         filename = name of file containing descriptor definitions
@ -56,13 +56,15 @@ and it is also straightforward to add new descriptor styles.
 The compute *mliap* command must be followed by two keywords
 *model* and *descriptor* in either order.

-The *model* keyword is followed by a model style, currently limited to
-either *linear* or *quadratic*.
+The *model* keyword is followed by the model style (*linear*, *quadratic* or *mliappy*).
+The *mliappy* model is only available
+if lammps is built with MLIAPPY package.

 The *descriptor* keyword is followed by a descriptor style, and additional arguments.
-Currently the only descriptor style is *sna*, indicating the bispectrum component
-descriptors used by the Spectral Neighbor Analysis Potential (SNAP) potentials of
-:doc:`pair_style snap <pair_snap>`.
+The compute currently supports just one descriptor style, but it is
+is straightforward to add new descriptor styles.
+The SNAP descriptor style *sna* is the same as that used by :doc:`pair_style snap <pair_snap>`,
+including the linear, quadratic, and chem variants.
 A single additional argument specifies the descriptor filename
 containing the parameters and setting used by the SNAP descriptor.
 The descriptor filename usually ends in the *.mliap.descriptor* extension.
@ -162,9 +164,10 @@ potentials, see the examples in `FitSNAP <https://github.com/FitSNAP/FitSNAP>`_.
 Restrictions
 """"""""""""

-This compute is part of the MLIAP package.  It is only enabled if
-LAMMPS was built with that package.  In addition, building LAMMPS with the MLIAP package
+This compute is part of the MLIAP package.  It is only enabled if LAMMPS
+was built with that package. In addition, building LAMMPS with the MLIAP package
 requires building LAMMPS with the SNAP package.
+The *mliappy* model requires building LAMMPS with the PYTHON package.
 See the :doc:`Build package <Build_package>` doc page for more info.

 Related commands
--- a/doc/src/compute_orientorder_atom.rst
+++ b/doc/src/compute_orientorder_atom.rst
@ -115,8 +115,8 @@ The optional keyword *chunksize* is only applicable when using the
 the KOKKOS package and is ignored otherwise. This keyword controls
 the number of atoms in each pass used to compute the bond-orientational
 order parameters and is used to avoid running out of memory. For example
-if there are 4000 atoms in the simulation and the *chunksize*
-is set to 2000, the parameter calculation will be broken up
+if there are 32768 atoms in the simulation and the *chunksize*
+is set to 16384, the parameter calculation will be broken up
 into two passes.

 The value of :math:`Q_l` is set to zero for atoms not in the
@ -193,7 +193,7 @@ Default

 The option defaults are *cutoff* = pair style cutoff, *nnn* = 12,
 *degrees* = 5 4 6 8 10 12 i.e. :math:`Q_4`, :math:`Q_6`, :math:`Q_8`, :math:`Q_{10}`, and :math:`Q_{12}`,
-*wl* = no, *wl/hat* = no, *components* off, and *chunksize* = 2000
+*wl* = no, *wl/hat* = no, *components* off, and *chunksize* = 16384

 ----------

--- a/doc/src/compute_pressure.rst
+++ b/doc/src/compute_pressure.rst
@ -122,8 +122,11 @@ Output info
 This compute calculates a global scalar (the pressure) and a global
 vector of length 6 (pressure tensor), which can be accessed by indices
 1-6.  These values can be used by any command that uses global scalar
-or vector values from a compute as input.  See the :doc:`Howto output <Howto_output>` doc page for an overview of LAMMPS output
-options.
+or vector values from a compute as input.  See the :doc:`Howto output
+<Howto_output>` doc page for an overview of LAMMPS output options.
+
+The ordering of values in the symmetric pressure tensor is as follows:
+pxx, pyy, pzz, pxy, pxz, pyz.

 The scalar and vector values calculated by this compute are
 "intensive".  The scalar and vector values will be in pressure
--- a/doc/src/compute_reduce_chunk.rst
+++ b/doc/src/compute_reduce_chunk.rst
@ -30,7 +30,7 @@ Examples

 .. code-block:: LAMMPS

-   compute 1 all reduce/chunk/atom mychunk min c_cluster
+   compute 1 all reduce/chunk mychunk min c_cluster

 Description
 """""""""""
--- a/doc/src/compute_stress_atom.rst
+++ b/doc/src/compute_stress_atom.rst
@ -216,6 +216,11 @@ an identical manner to compute *stress/atom*.  See the :doc:`Howto
 output <Howto_output>` doc page for an overview of LAMMPS output
 options.

+The ordering of the 6 columns for *stress/atom* is as follows: xx, yy,
+zz, xy, xz, yz.  The ordering of the 9 columns for
+*centroid/stress/atom* is as follows: xx, yy, zz, xy, xz, yz, yx, zx,
+zy.
+
 The per-atom array values will be in pressure\*volume :doc:`units
 <units>` as discussed above.

--- a/doc/src/create_bonds.rst
+++ b/doc/src/create_bonds.rst
@ -125,6 +125,16 @@ cannot appear in the neighbor list, to avoid creation of duplicate
 bonds.  The neighbor list for all atom type pairs must also extend to
 a distance that encompasses the *rmax* for new bonds to create.

+.. note::
+
+   If you want to create bonds between pairs of 1-3 or 1-4 atoms in
+   the current bond topology, then you need to use :doc:`special_bonds
+   lj 0 1 1 <special_bonds>` to insure those pairs appear in the
+   neighbor list.  They will not appear with the default special_bonds
+   settings which are zero for 1-2, 1-3, and 1-4 atoms.  1-3 or 1-4
+   atoms are those which are 2 hops or 3 hops apart in the bond
+   topology.
+
 An additional requirement for this style is that your system must be
 ready to perform a simulation.  This means, for example, that all
 :doc:`pair_style <pair_style>` coefficients be set via the
--- a/doc/src/fix.rst
+++ b/doc/src/fix.rst
@ -364,6 +364,8 @@ accelerated styles exist.
 * :doc:`temp/rescale <fix_temp_rescale>` - temperature control by velocity rescaling
 * :doc:`temp/rescale/eff <fix_temp_rescale_eff>` - temperature control by velocity rescaling in the electron force field model
 * :doc:`tfmc <fix_tfmc>` - perform force-bias Monte Carlo with time-stamped method
+* :doc:`tgnvt/drude <fix_tgnh_drude>` - NVT time integration for Drude polarizable model via temperature-grouped Nose-Hoover
+* :doc:`tgnpt/drude <fix_tgnh_drude>` - NPT time integration for Drude polarizable model via temperature-grouped Nose-Hoover
 * :doc:`thermal/conductivity <fix_thermal_conductivity>` - Muller-Plathe kinetic energy exchange for thermal conductivity calculation
 * :doc:`ti/spring <fix_ti_spring>` -
 * :doc:`tmd <fix_tmd>` - guide a group of atoms to a new configuration
--- a/doc/src/fix_addforce.rst
+++ b/doc/src/fix_addforce.rst
@ -118,32 +118,38 @@ converge properly.
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-No information about this fix is written to :doc:`binary restart files <restart>`.
+No information about this fix is written to :doc:`binary restart files
+<restart>`.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the potential "energy" inferred by the added force to the
-system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.  This is a fictitious quantity but is
-needed so that the :doc:`minimize <minimize>` command can include the
-forces added by this fix in a consistent manner.  I.e. there is a
-decrease in potential energy when atoms move in the direction of the
-added force.
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the potential energy inferred by the added force to
+the global potential energy of the system as part of
+:doc:`thermodynamic output <thermo_style>`.  The default setting for
+this fix is :doc:`fix_modify energy no <fix_modify>`.  Note that this
+energy is a fictitious quantity but is needed so that the
+:doc:`minimize <minimize>` command can include the forces added by
+this fix in a consistent manner.  I.e. there is a decrease in
+potential energy when atoms move in the direction of the added force.

-The :doc:`fix_modify <fix_modify>` *virial* option is supported by this
-fix to add the contribution due to the added forces on atoms to the
-system's virial as part of :doc:`thermodynamic output <thermo_style>`.
-The default is *virial no*
+The :doc:`fix_modify <fix_modify>` *virial* option is supported by
+this fix to add the contribution due to the added forces on atoms to
+both the global pressure and per-atom stress of the system via the
+:doc:`compute pressure <compute_pressure>` and :doc:`compute
+stress/atom <compute_stress_atom>` commands.  The former can be
+accessed by :doc:`thermodynamic output <thermo_style>`.  The default
+setting for this fix is :doc:`fix_modify virial no <fix_modify>`.

 The :doc:`fix_modify <fix_modify>` *respa* option is supported by this
-fix. This allows to set at which level of the :doc:`r-RESPA <run_style>`
-integrator the fix is adding its forces. Default is the outermost
-level.
+fix. This allows to set at which level of the :doc:`r-RESPA
+<run_style>` integrator the fix is adding its forces. Default is the
+outermost level.

 This fix computes a global scalar and a global 3-vector of forces,
-which can be accessed by various :doc:`output commands <Howto_output>`.
-The scalar is the potential energy discussed above.  The vector is the
-total force on the group of atoms before the forces on individual
-atoms are changed by the fix.  The scalar and vector values calculated
-by this fix are "extensive".
+which can be accessed by various :doc:`output commands
+<Howto_output>`.  The scalar is the potential energy discussed above.
+The vector is the total force on the group of atoms before the forces
+on individual atoms are changed by the fix.  The scalar and vector
+values calculated by this fix are "extensive".

 No parameter of this fix can be used with the *start/stop* keywords of
 the :doc:`run <run>` command.
--- a/doc/src/fix_addtorque.rst
+++ b/doc/src/fix_addtorque.rst
@ -55,13 +55,15 @@ Restart, fix_modify, output, run start/stop, minimize info

 No information about this fix is written to :doc:`binary restart files <restart>`.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the potential "energy" inferred by the added forces to the
-system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.  This is a fictitious quantity but is
-needed so that the :doc:`minimize <minimize>` command can include the
-forces added by this fix in a consistent manner.  I.e. there is a
-decrease in potential energy when atoms move in the direction of the
-added forces.
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the potential "energy" inferred by the added torques
+to the global potential energy of the system as part of
+:doc:`thermodynamic output <thermo_style>`.  The default setting for
+this fix is :doc:`fix_modify energy no <fix_modify>`.  Note that this
+is a fictitious quantity but is needed so that the :doc:`minimize
+<minimize>` command can include the forces added by this fix in a
+consistent manner.  I.e. there is a decrease in potential energy when
+atoms move in the direction of the added forces.

 The :doc:`fix_modify <fix_modify>` *respa* option is supported by
 this fix. This allows to set at which level of the :doc:`r-RESPA <run_style>`
@ -78,16 +80,28 @@ No parameter of this fix can be used with the *start/stop* keywords of
 the :doc:`run <run>` command.

 The forces due to this fix are imposed during an energy minimization,
-invoked by the :doc:`minimize <minimize>` command.  You should not
-specify force components with a variable that has time-dependence for
-use with a minimizer, since the minimizer increments the timestep as
-the iteration count during the minimization.
+invoked by the :doc:`minimize <minimize>` command.
+
+.. note::
+
+   If you want the fictitious potential energy associated with the
+   added forces to be included in the total potential energy of the
+   system (the quantity being minimized), you MUST enable the
+   :doc:`fix_modify <fix_modify>` *energy* option for this fix.
+
+.. note::
+
+   You should not specify force components with a variable that has
+   time-dependence for use with a minimizer, since the minimizer
+   increments the timestep as the iteration count during the
+   minimization.

 Restrictions
 """"""""""""

 This fix is part of the USER-MISC package.  It is only enabled if
-LAMMPS was built with that package.  See the :doc:`Build package <Build_package>` doc page for more info.
+LAMMPS was built with that package.  See the :doc:`Build package
+<Build_package>` doc page for more info.

 Related commands
 """"""""""""""""
--- a/doc/src/fix_atc.rst
+++ b/doc/src/fix_atc.rst
@ -122,10 +122,22 @@ Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

 No information about this fix is written to :doc:`binary restart files
-<restart>`.  The :doc:`fix_modify <fix_modify>` options relevant to this
-fix are listed below.  No global scalar or vector or per-atom quantities
-are stored by this fix for access by various :doc:`output commands
-<Howto_output>`.  No parameter of this fix can be used with the
+<restart>`.
+
+The :doc:`fix_modify <fix_modify>` *energy* option is not supported by
+this fix, but this fix does add the kinetic energy imparted to atoms
+by the momentum coupling mode of the AtC package to the global
+potential energy of the system as part of :doc:`thermodynamic output
+<thermo_style>`.
+
+Additional :doc:`fix_modify <fix_modify>` options relevant to this
+fix are listed below.
+
+This fix computes a global scalar which can be accessed by various
+:doc:`output commands <Howto_output>`.  The scalar is the energy
+discussed in the previous paragraph.  The scalar value is "extensive".
+
+No parameter of this fix can be used with the
 *start/stop* keywords of the :doc:`run <run>` command.  This fix is not
 invoked during :doc:`energy minimization <minimize>`.

@ -145,10 +157,10 @@ one, e.g. nve, nvt, etc. In addition, currently:
 Related commands
 """"""""""""""""

-After specifying this fix in your input script, several other
-:doc:`fix_modify <fix_modify>` commands are used to setup the problem,
-e.g. define the finite element mesh and prescribe initial and boundary
-conditions.
+After specifying this fix in your input script, several
+:doc:`fix_modify AtC <fix_modify>` commands are used to setup the
+problem, e.g. define the finite element mesh and prescribe initial and
+boundary conditions.  Each of these options has its own doc page.

 *fix_modify* commands for setup:

@ -240,7 +252,8 @@ miscellaneous *fix_modify* commands:
 * :doc:`fix_modify AtC remove_species <atc_remove_species>`
 * :doc:`fix_modify AtC remove_molecule <atc_remove_molecule>`

-Note: a set of example input files with the attendant material files are included in the ``examples/USER/atc`` folders.
+Note: a set of example input files with the attendant material files
+are included in the ``examples/USER/atc`` folders.

 Default
 """""""
@ -252,30 +265,52 @@ For detailed exposition of the theory and algorithms please see:

 .. _Wagner:

-**(Wagner)** Wagner, GJ; Jones, RE; Templeton, JA; Parks, MA, "An atomistic-to-continuum coupling method for heat transfer in solids." Special Issue of Computer Methods and Applied Mechanics (2008) 197:3351.
+**(Wagner)** Wagner, GJ; Jones, RE; Templeton, JA; Parks, MA, "An
+ atomistic-to-continuum coupling method for heat transfer in solids."
+ Special Issue of Computer Methods and Applied Mechanics (2008)
+ 197:3351.

 .. _Zimmeman2004:

-**(Zimmerman2004)** Zimmerman, JA; Webb, EB; Hoyt, JJ;. Jones, RE; Klein, PA; Bammann, DJ, "Calculation of stress in atomistic simulation." Special Issue of Modelling and Simulation in Materials Science and Engineering (2004), 12:S319.
+**(Zimmerman2004)** Zimmerman, JA; Webb, EB; Hoyt, JJ;. Jones, RE;
+ Klein, PA; Bammann, DJ, "Calculation of stress in atomistic
+ simulation." Special Issue of Modelling and Simulation in Materials
+ Science and Engineering (2004), 12:S319.

 .. _Zimmerman2010:

-**(Zimmerman2010)** Zimmerman, JA; Jones, RE; Templeton, JA, "A material frame approach for evaluating continuum variables in atomistic simulations." Journal of Computational Physics (2010), 229:2364.
+**(Zimmerman2010)** Zimmerman, JA; Jones, RE; Templeton, JA, "A
+ material frame approach for evaluating continuum variables in
+ atomistic simulations." Journal of Computational Physics (2010),
+ 229:2364.

 .. _Templeton2010:

-**(Templeton2010)** Templeton, JA; Jones, RE; Wagner, GJ, "Application of a field-based method to spatially varying thermal transport problems in molecular dynamics." Modelling and Simulation in Materials Science and Engineering (2010), 18:085007.
+**(Templeton2010)** Templeton, JA; Jones, RE; Wagner, GJ, "Application
+ of a field-based method to spatially varying thermal transport
+ problems in molecular dynamics." Modelling and Simulation in
+ Materials Science and Engineering (2010), 18:085007.

 .. _Jones:

-**(Jones)** Jones, RE; Templeton, JA; Wagner, GJ; Olmsted, D; Modine, JA, "Electron transport enhanced molecular dynamics for metals and semi-metals." International Journal for Numerical Methods in Engineering (2010), 83:940.
+**(Jones)** Jones, RE; Templeton, JA; Wagner, GJ; Olmsted, D; Modine,
+ JA, "Electron transport enhanced molecular dynamics for metals and
+ semi-metals." International Journal for Numerical Methods in
+ Engineering (2010), 83:940.

 .. _Templeton2011:

-**(Templeton2011)** Templeton, JA; Jones, RE; Lee, JW; Zimmerman, JA; Wong, BM, "A long-range electric field solver for molecular dynamics based on atomistic-to-continuum modeling." Journal of Chemical Theory and Computation (2011), 7:1736.
+**(Templeton2011)** Templeton, JA; Jones, RE; Lee, JW; Zimmerman, JA;
+ Wong, BM, "A long-range electric field solver for molecular dynamics
+ based on atomistic-to-continuum modeling." Journal of Chemical Theory
+ and Computation (2011), 7:1736.

 .. _Mandadapu:

-**(Mandadapu)** Mandadapu, KK; Templeton, JA; Lee, JW, "Polarization as a field variable from molecular dynamics simulations." Journal of Chemical Physics (2013), 139:054115.
+**(Mandadapu)** Mandadapu, KK; Templeton, JA; Lee, JW, "Polarization
+ as a field variable from molecular dynamics simulations." Journal of
+ Chemical Physics (2013), 139:054115.

-Please refer to the standard finite element (FE) texts, e.g. T.J.R Hughes " The finite element method ", Dover 2003, for the basics of FE simulation.
+Please refer to the standard finite element (FE) texts, e.g. T.J.R
+Hughes " The finite element method ", Dover 2003, for the basics of FE
+simulation.
--- a/doc/src/fix_ave_correlate.rst
+++ b/doc/src/fix_ave_correlate.rst
@ -93,7 +93,7 @@ from a compute, fix, or variable, then see the :doc:`fix ave/chunk <fix_ave_chun
 :doc:`fix ave/histo <fix_ave_histo>` commands.  If you wish to convert a
 per-atom quantity into a single global value, see the :doc:`compute reduce <compute_reduce>` command.

-The input values must either be all scalars.  What kinds of
+The input values must be all scalars.  What kinds of
 correlations between input values are calculated is determined by the
 *type* keyword as discussed below.

--- a/doc/src/fix_bocs.rst
+++ b/doc/src/fix_bocs.rst
@ -75,6 +75,39 @@ Note that *V_avg* and *Coeff_i* should all be in the proper units, e.g. if you
 are using *units real*\ , *V_avg* should be in cubic angstroms, and the
 coefficients should all be in atmospheres \* cubic angstroms.

+----------
+
+Restart, fix_modify, output, run start/stop, minimize info
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+This fix writes the cumulative global energy change to :doc:`binary
+restart files <restart>`.  See the :doc:`read_restart <read_restart>`
+command for info on how to re-specify a fix in an input script that
+reads a restart file, so that the fix continues in an uninterrupted
+fashion.
+
+The :doc:`fix_modify <fix_modify>` *temp* option is supported by this
+fix.  You can use it to assign a temperature :doc:`compute <compute>`
+you have defined to this fix which will be used in its thermostatting
+procedure, as described above.  For consistency, the group used by
+this fix and by the compute should be the same.
+
+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*.  See the :doc:`thermo_style <thermo_style>`
+doc page for details.
+
+This fix computes a global scalar which can be accessed by various
+:doc:`output commands <Howto_output>`.  The scalar is the same
+cumulative energy change due to this fix described in the previous
+paragraph.  The scalar value calculated by this fix is "extensive".
+
+This fix can ramp its target temperature over multiple runs, using the
+*start* and *stop* keywords of the :doc:`run <run>` command.  See the
+:doc:`run <run>` command for details of how to do this.
+
+This fix is not invoked during :doc:`energy minimization <minimize>`.
+
 Restrictions
 """"""""""""

--- a/doc/src/fix_bond_react.rst
+++ b/doc/src/fix_bond_react.rst
@ -41,7 +41,7 @@ Syntax
 * template-ID(post-reacted) = ID of a molecule template containing post-reaction topology
 * map_file = name of file specifying corresponding atom-IDs in the pre- and post-reacted templates
 * zero or more individual keyword/value pairs may be appended to each react argument
-* individual_keyword = *prob* or *max_rxn* or *stabilize_steps* or *custom_charges*
+* individual_keyword = *prob* or *max_rxn* or *stabilize_steps* or *custom_charges* or *molecule* or *modify_create*

  .. parsed-literal::

@ -59,6 +59,12 @@ Syntax
           off = allow both inter- and intramolecular reactions (default)
           inter = search for reactions between molecules with different IDs
           intra = search for reactions within the same molecule
+         *modify_create* keyword values
+           *fit* value = *all* or *fragmentID*
+             all = use all eligible atoms for create-atoms fit (default)
+             fragmentID = ID of molecule fragment used for create-atoms fit
+           *overlap* value = R
+             R = only insert atom/molecule if further than R from existing particles (distance units)

 Examples
 """"""""
@ -89,7 +95,9 @@ documentation. Topology changes are defined in pre- and post-reaction
 molecule templates and can include creation and deletion of bonds,
 angles, dihedrals, impropers, bond types, angle types, dihedral types,
 atom types, or atomic charges. In addition, reaction by-products or
-other molecules can be identified and deleted.
+other molecules can be identified and deleted. Finally, atoms can be
+created and inserted at specific positions relative to the reaction
+site.

 Fix bond/react does not use quantum mechanical (eg. fix qmmm) or
 pairwise bond-order potential (eg. Tersoff or AIREBO) methods to
@ -262,14 +270,14 @@ command page.

 The post-reacted molecule template contains a sample of the reaction
 site and its surrounding topology after the reaction has occurred. It
-must contain the same number of atoms as the pre-reacted template. A
-one-to-one correspondence between the atom IDs in the pre- and
-post-reacted templates is specified in the map file as described
-below. Note that during a reaction, an atom, bond, etc. type may
-change to one that was previously not present in the simulation. These
-new types must also be defined during the setup of a given simulation.
-A discussion of correctly handling this is also provided on the
-:doc:`molecule <molecule>` command page.
+must contain the same number of atoms as the pre-reacted template
+(unless there are created atoms). A one-to-one correspondence between
+the atom IDs in the pre- and post-reacted templates is specified in
+the map file as described below. Note that during a reaction, an atom,
+bond, etc. type may change to one that was previously not present in
+the simulation. These new types must also be defined during the setup
+of a given simulation. A discussion of correctly handling this is also
+provided on the :doc:`molecule <molecule>` command page.

 .. note::

@ -283,7 +291,7 @@ A discussion of correctly handling this is also provided on the
 The map file is a text document with the following format:

 A map file has a header and a body. The header of map file the
-contains one mandatory keyword and four optional keywords. The
+contains one mandatory keyword and five optional keywords. The
 mandatory keyword is 'equivalences':

 .. parsed-literal::
@ -296,11 +304,12 @@ The optional keywords are 'edgeIDs', 'deleteIDs', 'chiralIDs' and
 .. parsed-literal::

   N *edgeIDs* = # of edge atoms N in the pre-reacted molecule template
-   N *deleteIDs* = # of atoms N that are specified for deletion
-   N *chiralIDs* = # of specified chiral centers N
-   N *constraints* = # of specified reaction constraints N
+   N *deleteIDs* = # of atoms N that are deleted
+   N *createIDs* = # of atoms N that are created
+   N *chiralIDs* = # of chiral centers N
+   N *constraints* = # of reaction constraints N

-The body of the map file contains two mandatory sections and four
+The body of the map file contains two mandatory sections and five
 optional sections. The first mandatory section begins with the keyword
 'InitiatorIDs' and lists the two atom IDs of the initiator atom pair
 in the pre-reacted molecule template. The second mandatory section
@ -313,8 +322,10 @@ the keyword 'EdgeIDs' and lists the atom IDs of edge atoms in the
 pre-reacted molecule template. The second optional section begins with
 the keyword 'DeleteIDs' and lists the atom IDs of pre-reaction
 template atoms to delete. The third optional section begins with the
+keyword 'CreateIDs' and lists the atom IDs of the post-reaction
+template atoms to create. The fourth optional section begins with the
 keyword 'ChiralIDs' lists the atom IDs of chiral atoms whose
-handedness should be enforced. The fourth optional section begins with
+handedness should be enforced. The fifth optional section begins with
 the keyword 'Constraints' and lists additional criteria that must be
 satisfied in order for the reaction to occur. Currently, there are
 five types of constraints available, as discussed below: 'distance',
@ -353,6 +364,38 @@ A sample map file is given below:

 ----------

+A user-specified set of atoms can be deleted by listing their
+pre-reaction template IDs in the DeleteIDs section. A deleted atom
+must still be included in the post-reaction molecule template, in
+which it cannot be bonded to an atom that is not deleted. In addition
+to deleting unwanted reaction by-products, this feature can be used to
+remove specific topologies, such as small rings, that may be otherwise
+indistinguishable.
+
+Atoms can be created by listing their post-reaction template IDs in
+the CreateIDs section. A created atom should not be included in the
+pre-reaction template. The inserted positions of created atoms are
+determined by the coordinates of the post-reaction template, after
+optimal translation and rotation of the post-reaction template to the
+reaction site (using a fit with atoms that are neither created nor
+deleted). The *modify_create* keyword can be used to modify the
+default behavior when creating atoms. The *modify_create* keyword has
+two sub-keywords, *fit* and *overlap*. One or more of the sub-keywords
+may be used after the *modify_create* keyword. The *fit* sub-keyword
+can be used to specify which post-reaction atoms are used for the
+optimal translation and rotation of the post-reaction template. The
+*fragmentID* value of the *fit* sub-keyword must be the name of a
+molecule fragment defined in the post-reaction :doc:`molecule
+<molecule>` template, and only atoms in this fragment are used for the
+fit. Atoms are created only if no current atom in the simulation is
+within a distance R of any created atom, including the effect of
+periodic boundary conditions if applicable. R is defined by the
+*overlap* sub-keyword. Note that the default value for R is 0.0, which
+will allow atoms to strongly overlap if you are inserting where other
+atoms are present. The velocity of each created atom is initialized in
+a random direction with a magnitude calculated from the instantaneous
+temperature of the reaction site.
+
 The handedness of atoms that are chiral centers can be enforced by
 listing their IDs in the ChiralIDs section. A chiral atom must be
 bonded to four atoms with mutually different atom types. This feature
@ -457,6 +500,23 @@ example, the molecule fragment could consist of only the backbone
 atoms of a polymer chain. This constraint can be used to enforce a
 specific relative position and orientation between reacting molecules.

+By default, all constraints must be satisfied for the reaction to
+occur. In other words, constraints are evaluated as a series of
+logical values using the logical AND operator "&&". More complex logic
+can be achieved by explicitly adding the logical AND operator "&&" or
+the logical OR operator "||" after a given constraint command. If a
+logical operator is specified after a constraint, it must be placed
+after all constraint parameters, on the same line as the constraint
+(one per line). Similarly, parentheses can be used to group
+constraints. The expression that results from concatenating all
+constraints should be a valid logical expression that can be read by
+the :doc:`variable <variable>` command after converting each
+constraint to a logical value. Because exactly one constraint is
+allowed per line, having a valid logical expression implies that left
+parentheses "(" should only appear before a constraint, and right
+parentheses ")" should only appear after a constraint and before any
+logical operator.
+
 Once a reaction site has been successfully identified, data structures
 within LAMMPS that store bond topology are updated to reflect the
 post-reacted molecule template. All force fields with fixed bonds,
@ -511,15 +571,6 @@ the same molecule ID are considered for the reaction.

 A few other considerations:

-Many reactions result in one or more atoms that are considered
-unwanted by-products. Therefore, bond/react provides the option to
-delete a user-specified set of atoms. These pre-reaction atoms are
-identified in the map file. A deleted atom must still be included in
-the post-reaction molecule template, in which it cannot be bonded to
-an atom that is not deleted. In addition to deleting unwanted reaction
-by-products, this feature can be used to remove specific topologies,
-such as small rings, that may be otherwise indistinguishable.
-
 Optionally, you can enforce additional behaviors on reacting atoms.
 For example, it may be beneficial to force reacting atoms to remain at
 a certain temperature. For this, you can use the internally-created
@ -593,14 +644,14 @@ Default
 """""""

 The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60,
-reset_mol_ids = yes, custom_charges = no, molecule = off
+reset_mol_ids = yes, custom_charges = no, molecule = off, modify_create = no

 ----------

 .. _Gissinger:

-**(Gissinger)** Gissinger, Jensen and Wise, Polymer, 128, 211 (2017).
+**(Gissinger2017)** Gissinger, Jensen and Wise, Polymer, 128, 211-217 (2017).

 .. _Gissinger2020:

-**(Gissinger)** Gissinger, Jensen and Wise, Macromolecules (2020, in press).
+**(Gissinger2020)** Gissinger, Jensen and Wise, Macromolecules, 53, 22, 9953-9961 (2020).
--- a/doc/src/fix_client_md.rst
+++ b/doc/src/fix_client_md.rst
@ -47,14 +47,15 @@ for running *ab initio* MD with quantum forces.
 The group associated with this fix is ignored.

 The protocol and :doc:`units <units>` for message format and content
-that LAMMPS exchanges with the server code is defined on the :doc:`server md <server_md>` doc page.
+that LAMMPS exchanges with the server code is defined on the
+:doc:`server md <server_md>` doc page.

 Note that when using LAMMPS as an MD client, your LAMMPS input script
 should not normally contain force field commands, like a
 :doc:`pair_style <pair_style>`, :doc:`bond_style <bond_style>`, or
-:doc:`kspace_style <kspace_style>` command.  However it is possible for
-a server code to only compute a portion of the full force-field, while
-LAMMPS computes the remaining part.  Your LAMMPS script can also
+:doc:`kspace_style <kspace_style>` command.  However it is possible
+for a server code to only compute a portion of the full force-field,
+while LAMMPS computes the remaining part.  Your LAMMPS script can also
 specify boundary conditions or force constraints in the usual way,
 which will be added to the per-atom forces returned by the server
 code.
@ -69,16 +70,21 @@ LAMMPS and another code in tandem to perform a coupled simulation.
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-No information about this fix is written to :doc:`binary restart files <restart>`.
+No information about this fix is written to :doc:`binary restart files
+<restart>`.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the potential energy computed by the server application to
-the system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the potential energy set by the server application to
+the global potential energy of the system as part of
+:doc:`thermodynamic output <thermo_style>`.  The default setting for
+this fix is :doc:`fix_modify energy yes <fix_modify>`.

-The :doc:`fix_modify <fix_modify>` *virial* option is supported by this
-fix to add the server application's contribution to the system's
-virial as part of :doc:`thermodynamic output <thermo_style>`.  The
-default is *virial yes*
+The :doc:`fix_modify <fix_modify>` *virial* option is supported by
+this fix to add the contribution computed by the server application to
+the global pressure of the system via the :doc:`compute pressure
+<compute_pressure>` command.  This can be accessed by
+:doc:`thermodynamic output <thermo_style>`.  The default setting for
+this fix is :doc:`fix_modify virial yes <fix_modify>`.

 This fix computes a global scalar which can be accessed by various
 :doc:`output commands <Howto_output>`.  The scalar is the potential
@ -86,13 +92,16 @@ energy discussed above.  The scalar value calculated by this fix is
 "extensive".

 No parameter of this fix can be used with the *start/stop* keywords of
-the :doc:`run <run>` command.  This fix is not invoked during :doc:`energy minimization <minimize>`.
+the :doc:`run <run>` command.
+
+This fix is not invoked during :doc:`energy minimization <minimize>`.

 Restrictions
 """"""""""""

 This fix is part of the MESSAGE package.  It is only enabled if LAMMPS
-was built with that package.  See the :doc:`Build package <Build_package>` doc page for more info.
+was built with that package.  See the :doc:`Build package
+<Build_package>` doc page for more info.

 A script that uses this command must also use the
 :doc:`message <message>` command to setup and shut down the messaging
--- a/doc/src/fix_cmap.rst
+++ b/doc/src/fix_cmap.rst
@ -29,11 +29,12 @@ Description
 This command enables CMAP cross-terms to be added to simulations which
 use the CHARMM force field.  These are relevant for any CHARMM model
 of a peptide or protein sequences that is 3 or more amino-acid
-residues long; see :ref:`(Buck) <Buck>` and :ref:`(Brooks) <Brooks2>` for details,
-including the analytic energy expressions for CMAP interactions.  The
-CMAP cross-terms add additional potential energy contributions to pairs
-of overlapping phi-psi dihedrals of amino-acids, which are important
-to properly represent their conformational behavior.
+residues long; see :ref:`(Buck) <Buck>` and :ref:`(Brooks) <Brooks2>`
+for details, including the analytic energy expressions for CMAP
+interactions.  The CMAP cross-terms add additional potential energy
+contributions to pairs of overlapping phi-psi dihedrals of
+amino-acids, which are important to properly represent their
+conformational behavior.

 The examples/cmap directory has a sample input script and data file
 for a small peptide, that illustrates use of the fix cmap command.
@ -93,19 +94,27 @@ the note below about how to include the CMAP energy when performing an
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-This fix writes the list of CMAP cross-terms to :doc:`binary restart files <restart>`.  See the :doc:`read_restart <read_restart>` command
+This fix writes the list of CMAP cross-terms to :doc:`binary restart
+files <restart>`.  See the :doc:`read_restart <read_restart>` command
 for info on how to re-specify a fix in an input script that reads a
 restart file, so that the operation of the fix continues in an
 uninterrupted fashion.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the potential "energy" of the CMAP interactions system's
-potential energy as part of :doc:`thermodynamic output <thermo_style>`.
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the potential energy of the CMAP interactions to both
+the global potential energy and peratom potential energies of the
+system as part of :doc:`thermodynamic output <thermo_style>` or
+output by the :doc:`compute pe/atom <compute_pe_atom>` command.  The
+default setting for this fix is :doc:`fix_modify energy yes
+<fix_modify>`.

-The :doc:`fix_modify <fix_modify>` *virial* option is supported by this
-fix to add the contribution due to the interaction between atoms to
-the system's virial as part of :doc:`thermodynamic output <thermo_style>`.
-The default is *virial yes*
+The :doc:`fix_modify <fix_modify>` *virial* option is supported by
+this fix to add the contribution due to the CMAP interactions to both
+the global pressure and per-atom stress of the system via the
+:doc:`compute pressure <compute_pressure>` and :doc:`compute
+stress/atom <compute_stress_atom>` commands.  The former can be
+accessed by :doc:`thermodynamic output <thermo_style>`.  The default
+setting for this fix is :doc:`fix_modify virial yes <fix_modify>`.

 This fix computes a global scalar which can be accessed by various
 :doc:`output commands <Howto_output>`.  The scalar is the potential
@ -121,8 +130,8 @@ invoked by the :doc:`minimize <minimize>` command.
 .. note::

   If you want the potential energy associated with the CMAP terms
-   forces to be included in the total potential energy of the system (the
-   quantity being minimized), you MUST enable the
+   forces to be included in the total potential energy of the system
+   (the quantity being minimized), you MUST not disable the
   :doc:`fix_modify <fix_modify>` *energy* option for this fix.

 Restrictions
--- a/doc/src/fix_colvars.rst
+++ b/doc/src/fix_colvars.rst
@ -35,12 +35,12 @@ Examples
 Description
 """""""""""

-This fix interfaces LAMMPS to the collective variables "Colvars"
-library, which allows to calculate potentials of mean force
-(PMFs) for any set of colvars, using different sampling methods:
-currently implemented are the Adaptive Biasing Force (ABF) method,
-metadynamics, Steered Molecular Dynamics (SMD) and Umbrella Sampling
-(US) via a flexible harmonic restraint bias.
+This fix interfaces LAMMPS to the collective variables (Colvars)
+library, which allows to calculate potentials of mean force (PMFs) for
+any set of colvars, using different sampling methods: currently
+implemented are the Adaptive Biasing Force (ABF) method, metadynamics,
+Steered Molecular Dynamics (SMD) and Umbrella Sampling (US) via a
+flexible harmonic restraint bias.

 This documentation describes only the fix colvars command itself and
 LAMMPS specific parts of the code.  The full documentation of the
@ -98,9 +98,11 @@ This fix writes the current status of the colvars module into
 mode status file that is written by the colvars module itself and the
 kind of information in both files is identical.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy change from the biasing force added by the fix
-to the system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the energy change from the biasing force added by
+Colvars to the global potential energy of the system as part of
+:doc:`thermodynamic output <thermo_style>`.  The default setting for
+this fix is :doc:`fix_modify energy no <fix_modify>`.

 The *fix_modify configfile <config file>* option allows to add settings
 from an additional config file to the colvars module. This option can
@ -113,15 +115,16 @@ in a pair of double quotes ("), or can span multiple lines when bracketed
 by a pair of triple double quotes (""", like python embedded documentation).

 This fix computes a global scalar which can be accessed by various
-:doc:`output commands <Howto_output>`.  The scalar is the cumulative
-energy change due to this fix.  The scalar value calculated by this
-fix is "extensive".
+:doc:`output commands <Howto_output>`.  The scalar is the Colvars
+energy mentioned above.  The scalar value calculated by this fix is
+"extensive".

 Restrictions
 """"""""""""

 This fix is part of the USER-COLVARS package.  It is only enabled if
-LAMMPS was built with that package.  See the :doc:`Build package <Build_package>` doc page for more info.
+LAMMPS was built with that package.  See the :doc:`Build package
+<Build_package>` doc page for more info.

 There can only be one colvars fix active at a time. Since the interface
 communicates only the minimum amount of information and colvars module
--- a/doc/src/fix_drude.rst
+++ b/doc/src/fix_drude.rst
@ -41,12 +41,12 @@ Restrictions
 """"""""""""

 This fix should be invoked before any other commands that implement
-the Drude oscillator model, such as :doc:`fix langevin/drude <fix_langevin_drude>`, :doc:`fix drude/transform <fix_drude_transform>`, :doc:`compute temp/drude <compute_temp_drude>`, :doc:`pair_style thole <pair_thole>`.
+the Drude oscillator model, such as :doc:`fix langevin/drude <fix_langevin_drude>`, :doc:`fix tgnvt/drude <fix_tgnh_drude>`, :doc:`fix drude/transform <fix_drude_transform>`, :doc:`compute temp/drude <compute_temp_drude>`, :doc:`pair_style thole <pair_thole>`.

 Related commands
 """"""""""""""""

-:doc:`fix langevin/drude <fix_langevin_drude>`, :doc:`fix drude/transform <fix_drude_transform>`, :doc:`compute temp/drude <compute_temp_drude>`, :doc:`pair_style thole <pair_thole>`
+:doc:`fix langevin/drude <fix_langevin_drude>`, :doc:`fix tgnvt/drude <fix_tgnh_drude>`, :doc:`fix drude/transform <fix_drude_transform>`, :doc:`compute temp/drude <compute_temp_drude>`, :doc:`pair_style thole <pair_thole>`

 Default
 """""""
--- a/doc/src/fix_efield.rst
+++ b/doc/src/fix_efield.rst
@ -100,9 +100,9 @@ minimize the orientation of dipoles in an applied electric field.

 The *energy* keyword specifies the name of an atom-style
 :doc:`variable <variable>` which is used to compute the energy of each
-atom as function of its position.  Like variables used for *ex*\ , *ey*\ ,
-*ez*\ , the energy variable is specified as v_name, where name is the
-variable name.
+atom as function of its position.  Like variables used for *ex*\ ,
+*ey*\ , *ez*\ , the energy variable is specified as v_name, where name
+is the variable name.

 Note that when the *energy* keyword is used during an energy
 minimization, you must insure that the formula defined for the
@ -117,31 +117,38 @@ minimization will not converge properly.
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-No information about this fix is written to :doc:`binary restart files <restart>`.
+No information about this fix is written to :doc:`binary restart files
+<restart>`.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the potential "energy" inferred by the added force due to
-the electric field to the system's potential energy as part of
-:doc:`thermodynamic output <thermo_style>`.  This is a fictitious
-quantity but is needed so that the :doc:`minimize <minimize>` command
-can include the forces added by this fix in a consistent manner.
-I.e. there is a decrease in potential energy when atoms move in the
-direction of the added force due to the electric field.
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the potential energy inferred by the added force due
+to the electric field to the global potential energy of the system as
+part of :doc:`thermodynamic output <thermo_style>`.  The default
+setting for this fix is :doc:`fix_modify energy no <fix_modify>`.
+Note that this energy is a fictitious quantity but is needed so that
+the :doc:`minimize <minimize>` command can include the forces added by
+this fix in a consistent manner.  I.e. there is a decrease in
+potential energy when atoms move in the direction of the added force
+due to the electric field.

-The :doc:`fix_modify <fix_modify>` *virial* option is supported by this
-fix to add the contribution due to the added forces on atoms to the
-system's virial as part of :doc:`thermodynamic output <thermo_style>`.
-The default is *virial no*
+The :doc:`fix_modify <fix_modify>` *virial* option is supported by
+this fix to add the contribution due to the added forces on atoms to
+both the global pressure and per-atom stress of the system via the
+:doc:`compute pressure <compute_pressure>` and :doc:`compute
+stress/atom <compute_stress_atom>` commands.  The former can be
+accessed by :doc:`thermodynamic output <thermo_style>`.  The default
+setting for this fix is :doc:`fix_modify virial no <fix_modify>`.

 The :doc:`fix_modify <fix_modify>` *respa* option is supported by this
-fix. This allows to set at which level of the :doc:`r-RESPA <run_style>`
-integrator the fix adding its forces. Default is the outermost level.
+fix. This allows to set at which level of the :doc:`r-RESPA
+<run_style>` integrator the fix adding its forces. Default is the
+outermost level.

 This fix computes a global scalar and a global 3-vector of forces,
-which can be accessed by various :doc:`output commands <Howto_output>`.
-The scalar is the potential energy discussed above.  The vector is the
-total force added to the group of atoms.  The scalar and vector values
-calculated by this fix are "extensive".
+which can be accessed by various :doc:`output commands
+<Howto_output>`.  The scalar is the potential energy discussed above.
+The vector is the total force added to the group of atoms.  The scalar
+and vector values calculated by this fix are "extensive".

 No parameter of this fix can be used with the *start/stop* keywords of
 the :doc:`run <run>` command.
--- a/doc/src/fix_external.rst
+++ b/doc/src/fix_external.rst
@ -84,7 +84,8 @@ code `Quest <quest_>`_.

 If mode is *pf/array* then the fix simply stores force values in an
 array.  The fix adds these forces to each atom in the group, once
-every *Napply* steps, similar to the way the :doc:`fix addforce <fix_addforce>` command works.
+every *Napply* steps, similar to the way the :doc:`fix addforce
+<fix_addforce>` command works.

 The name of the public force array provided by the FixExternal
 class is
@ -150,19 +151,27 @@ of properties that the caller code may want to communicate  to LAMMPS
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-No information about this fix is written to :doc:`binary restart files <restart>`.
+No information about this fix is written to :doc:`binary restart files
+<restart>`.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the potential "energy" set by the external driver to the
-system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.  This is a fictitious quantity but is
-needed so that the :doc:`minimize <minimize>` command can include the
-forces added by this fix in a consistent manner.  I.e. there is a
-decrease in potential energy when atoms move in the direction of the
-added force.
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the potential energy set by the external driver to
+both the global potential energy and peratom potential energies of the
+system as part of :doc:`thermodynamic output <thermo_style>` or output
+by the :doc:`compute pe/atom <compute_pe_atom>` command.  The default
+setting for this fix is :doc:`fix_modify energy yes <fix_modify>`.
+Note that this energy may be a fictitious quantity but it is needed so
+that the :doc:`minimize <minimize>` command can include the forces
+added by this fix in a consistent manner.  I.e. there is a decrease in
+potential energy when atoms move in the direction of the added force.

-The :doc:`fix_modify <fix_modify>` *virial* option is supported by this
-fix to add the contribution due to the interactions computed by the
-external program to the system's virial as part of :doc:`thermodynamic output <thermo_style>`. The default is *virial yes*
+The :doc:`fix_modify <fix_modify>` *virial* option is supported by
+this fix to add the contribution computed by the external program to
+both the global pressure and per-atom stress of the system via the
+:doc:`compute pressure <compute_pressure>` and :doc:`compute
+stress/atom <compute_stress_atom>` commands.  The former can be
+accessed by :doc:`thermodynamic output <thermo_style>`.  The default
+setting for this fix is :doc:`fix_modify virial yes <fix_modify>`.

 This fix computes a global scalar which can be accessed by various
 :doc:`output commands <Howto_output>`.  The scalar is the potential
@ -178,7 +187,7 @@ invoked by the :doc:`minimize <minimize>` command.

   If you want the fictitious potential energy associated with the
   added forces to be included in the total potential energy of the
-   system (the quantity being minimized), you MUST enable the
+   system (the quantity being minimized), you MUST not disable the
   :doc:`fix_modify <fix_modify>` *energy* option for this fix.

 Restrictions
--- a/doc/src/fix_ffl.rst
+++ b/doc/src/fix_ffl.rst
@ -76,28 +76,31 @@ Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

 The instantaneous values of the extended variables are written to
-:doc:`binary restart files <restart>`.  Because the state of the random
-number generator is not saved in restart files, this means you cannot
-do "exact" restarts with this fix, where the simulation continues on
-the same as if no restart had taken place. However, in a statistical
-sense, a restarted simulation should produce the same behavior.
-Note however that you should use a different seed each time you
-restart, otherwise the same sequence of random numbers will be used
-each time, which might lead to stochastic synchronization and
+:doc:`binary restart files <restart>`.  Because the state of the
+random number generator is not saved in restart files, this means you
+cannot do "exact" restarts with this fix, where the simulation
+continues on the same as if no restart had taken place. However, in a
+statistical sense, a restarted simulation should produce the same
+behavior.  Note however that you should use a different seed each time
+you restart, otherwise the same sequence of random numbers will be
+used each time, which might lead to stochastic synchronization and
 subtle artifacts in the sampling.

+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*.  See the :doc:`thermo_style <thermo_style>`
+doc page for details.
+
+This fix computes a global scalar which can be accessed by various
+:doc:`output commands <Howto_output>`.  The scalar is the same
+cumulative energy change due to this fix described in the previous
+paragraph.  The scalar value calculated by this fix is "extensive".
+
 This fix can ramp its target temperature over multiple runs, using the
 *start* and *stop* keywords of the :doc:`run <run>` command.  See the
 :doc:`run <run>` command for details of how to do this.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy change induced by Langevin thermostatting to the
-system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.
-
-This fix computes a global scalar which can be accessed by various
-:doc:`output commands <Howto_output>`.  The scalar is the cumulative
-energy change due to this fix.  The scalar value calculated by this
-fix is "extensive".
+This fix is not invoked during :doc:`energy minimization <minimize>`.

 Restrictions
 """"""""""""
--- a/doc/src/fix_flow_gauss.rst
+++ b/doc/src/fix_flow_gauss.rst
@ -55,17 +55,22 @@ momentum.
 GD applies an external fluctuating gravitational field that acts as a
 driving force to keep the system away from equilibrium. To maintain
 steady state, a profile-unbiased thermostat must be implemented to
-dissipate the heat that is added by the driving force. :doc:`Compute temp/profile <compute_temp_profile>` can be used to implement a
+dissipate the heat that is added by the driving force. :doc:`Compute
+temp/profile <compute_temp_profile>` can be used to implement a
 profile-unbiased thermostat.

 A common use of this fix is to compute a pressure drop across a pipe,
 pore, or membrane. The pressure profile can be computed in LAMMPS with
-:doc:`compute stress/atom <compute_stress_atom>` and :doc:`fix ave/chunk <fix_ave_chunk>`, or with the hardy method in :doc:`fix atc <fix_atc>`. Note that the simple :doc:`compute stress/atom <compute_stress_atom>` method is only accurate away
-from inhomogeneities in the fluid, such as fixed wall atoms. Further,
-the computed pressure profile must be corrected for the acceleration
+:doc:`compute stress/atom <compute_stress_atom>` and :doc:`fix
+ave/chunk <fix_ave_chunk>`, or with the hardy method in :doc:`fix atc
+<fix_atc>`. Note that the simple :doc:`compute stress/atom
+<compute_stress_atom>` method is only accurate away from
+inhomogeneities in the fluid, such as fixed wall atoms. Further, the
+computed pressure profile must be corrected for the acceleration
 applied by GD before computing a pressure drop or comparing it to
-other methods, such as the pump method :ref:`(Zhu) <Zhu>`. The pressure
-correction is discussed and described in :ref:`(Strong) <Strong>`.
+other methods, such as the pump method :ref:`(Zhu) <Zhu>`. The
+pressure correction is discussed and described in :ref:`(Strong)
+<Strong>`.

 For a complete example including the considerations discussed
 above, see the examples/USER/flow_gauss directory.
@ -102,14 +107,15 @@ computed by the fix will return zero.
   of wall atoms fixed, such as :doc:`fix spring/self <fix_spring_self>`.

 If this fix is used in a simulation with the :doc:`rRESPA <run_style>`
-integrator, the applied acceleration must be computed and applied at the same
-rRESPA level as the interactions between the flowing fluid and the obstacle.
-The rRESPA level at which the acceleration is applied can be changed using
-the :doc:`fix_modify <fix_modify>` *respa* option discussed below. If the
-flowing fluid and the obstacle interact through multiple interactions that are
-computed at different rRESPA levels, then there must be a separate flow/gauss
-fix for each level. For example, if the flowing fluid and obstacle interact
-through pairwise and long-range Coulomb interactions, which are computed at
+integrator, the applied acceleration must be computed and applied at
+the same rRESPA level as the interactions between the flowing fluid
+and the obstacle.  The rRESPA level at which the acceleration is
+applied can be changed using the :doc:`fix_modify <fix_modify>`
+*respa* option discussed below. If the flowing fluid and the obstacle
+interact through multiple interactions that are computed at different
+rRESPA levels, then there must be a separate flow/gauss fix for each
+level. For example, if the flowing fluid and obstacle interact through
+pairwise and long-range Coulomb interactions, which are computed at
 rRESPA levels 3 and 4, respectively, then there must be two separate
 flow/gauss fixes, one that specifies *fix_modify respa 3* and one with
 *fix_modify respa 4*.
@ -119,38 +125,49 @@ flow/gauss fixes, one that specifies *fix_modify respa 3* and one with
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-This fix is part of the USER-MISC package.  It is only enabled if
-LAMMPS was built with that package.  See the :doc:`Build package <Build_package>` doc page for more info.
+No information about this fix is written to :doc:`binary restart files
+<restart>`.

-No information about this fix is written to :doc:`binary restart files <restart>`.
-
-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to subtract the work done from the
-system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the potential energy added by the fix to the global
+potential energy of the system as part of :doc:`thermodynamic output
+<thermo_style>`.  The default setting for this fix is :doc:`fix_modify
+energy no <fix_modify>`.

 The :doc:`fix_modify <fix_modify>` *respa* option is supported by this
-fix. This allows the user to set at which level of the :doc:`rRESPA <run_style>`
-integrator the fix computes and adds the external acceleration. Default is the
-outermost level.
+fix. This allows the user to set at which level of the :doc:`rRESPA
+<run_style>` integrator the fix computes and adds the external
+acceleration. Default is the outermost level.

 This fix computes a global scalar and a global 3-vector of forces,
-which can be accessed by various :doc:`output commands <Howto_output>`.
-The scalar is the negative of the work done on the system, see above
-discussion.  The vector is the total force that this fix applied to
-the group of atoms on the current timestep.  The scalar and vector
-values calculated by this fix are "extensive".
+which can be accessed by various :doc:`output commands
+<Howto_output>`.  The scalar is the negative of the work done on the
+system, see the discussion above.  It is only calculated if the
+*energy* keyword is enabled or :doc:`fix_modify energy yes
+<fix_modify>` is set.
+
+The vector is the total force that this fix applied to the group of
+atoms on the current timestep.  The scalar and vector values
+calculated by this fix are "extensive".

 No parameter of this fix can be used with the *start/stop* keywords of
 the :doc:`run <run>` command.

+This fix is not invoked during :doc:`energy minimization <minimize>`.
+
 Restrictions
 """"""""""""
- none
+
+This fix is part of the USER-MISC package.  It is only enabled if
+LAMMPS was built with that package.  See the :doc:`Build package
+<Build_package>` doc page for more info.

 Related commands
 """"""""""""""""

-:doc:`fix addforce <fix_addforce>`, :doc:`compute temp/profile <compute_temp_profile>`, :doc:`velocity <velocity>`
+:doc:`fix addforce <fix_addforce>`,
+:doc:`compute temp/profile <compute_temp_profile>`,
+:doc:`velocity <velocity>`

 Default
 """""""
--- a/doc/src/fix_gcmc.rst
+++ b/doc/src/fix_gcmc.rst
@ -398,12 +398,13 @@ adds all inserted atoms of the specified type to the
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-This fix writes the state of the fix to :doc:`binary restart files <restart>`.  This includes information about the random
-number generator seed, the next timestep for MC exchanges,  the number
-of MC step attempts and successes etc.  See
-the :doc:`read_restart <read_restart>` command for info on how to
-re-specify a fix in an input script that reads a restart file, so that
-the operation of the fix continues in an uninterrupted fashion.
+This fix writes the state of the fix to :doc:`binary restart files
+<restart>`.  This includes information about the random number
+generator seed, the next timestep for MC exchanges, the number of MC
+step attempts and successes etc.  See the :doc:`read_restart
+<read_restart>` command for info on how to re-specify a fix in an
+input script that reads a restart file, so that the operation of the
+fix continues in an uninterrupted fashion.

 .. note::

@ -411,8 +412,8 @@ the operation of the fix continues in an uninterrupted fashion.
   after reading the restart with :doc:`reset_timestep <reset_timestep>`.
   The fix will try to detect it and stop with an error.

-None of the :doc:`fix_modify <fix_modify>` options are relevant to this
-fix.
+None of the :doc:`fix_modify <fix_modify>` options are relevant to
+this fix.

 This fix computes a global vector of length 8, which can be accessed
 by various :doc:`output commands <Howto_output>`.  The vector values are
@ -430,7 +431,8 @@ the following global cumulative quantities:
 The vector values calculated by this fix are "extensive".

 No parameter of this fix can be used with the *start/stop* keywords of
-the :doc:`run <run>` command.  This fix is not invoked during :doc:`energy minimization <minimize>`.
+the :doc:`run <run>` command.  This fix is not invoked during
+:doc:`energy minimization <minimize>`.

 Restrictions
 """"""""""""
--- a/doc/src/fix_gle.rst
+++ b/doc/src/fix_gle.rst
@ -103,28 +103,31 @@ Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

 The instantaneous values of the extended variables are written to
-:doc:`binary restart files <restart>`.  Because the state of the random
-number generator is not saved in restart files, this means you cannot
-do "exact" restarts with this fix, where the simulation continues on
-the same as if no restart had taken place. However, in a statistical
-sense, a restarted simulation should produce the same behavior.
-Note however that you should use a different seed each time you
-restart, otherwise the same sequence of random numbers will be used
-each time, which might lead to stochastic synchronization and
+:doc:`binary restart files <restart>`.  Because the state of the
+random number generator is not saved in restart files, this means you
+cannot do "exact" restarts with this fix, where the simulation
+continues on the same as if no restart had taken place. However, in a
+statistical sense, a restarted simulation should produce the same
+behavior.  Note however that you should use a different seed each time
+you restart, otherwise the same sequence of random numbers will be
+used each time, which might lead to stochastic synchronization and
 subtle artifacts in the sampling.

+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*.  See the :doc:`thermo_style <thermo_style>`
+doc page for details.
+
+This fix computes a global scalar which can be accessed by various
+:doc:`output commands <Howto_output>`.  The scalar is the same
+cumulative energy change due to this fix described in the previous
+paragraph.  The scalar value calculated by this fix is "extensive".
+
 This fix can ramp its target temperature over multiple runs, using the
 *start* and *stop* keywords of the :doc:`run <run>` command.  See the
 :doc:`run <run>` command for details of how to do this.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy change induced by Langevin thermostatting to the
-system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.
-
-This fix computes a global scalar which can be accessed by various
-:doc:`output commands <Howto_output>`.  The scalar is the cumulative
-energy change due to this fix.  The scalar value calculated by this
-fix is "extensive".
+This fix is not invoked during :doc:`energy minimization <minimize>`.

 Restrictions
 """"""""""""
--- a/doc/src/fix_gravity.rst
+++ b/doc/src/fix_gravity.rst
@ -103,23 +103,27 @@ Restart, fix_modify, output, run start/stop, minimize info

 No information about this fix is written to :doc:`binary restart files <restart>`.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the gravitational potential energy of the system to the
-system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the gravitational potential energy of the system to
+the global potential energy of the system as part of
+:doc:`thermodynamic output <thermo_style>`.  The default setting for
+this fix is :doc:`fix_modify energy no <fix_modify>`.

 The :doc:`fix_modify <fix_modify>` *respa* option is supported by this
-fix. This allows to set at which level of the :doc:`r-RESPA <run_style>`
-integrator the fix is adding its forces. Default is the outermost level.
+fix. This allows to set at which level of the :doc:`r-RESPA
+<run_style>` integrator the fix is adding its forces. Default is the
+outermost level.

 This fix computes a global scalar which can be accessed by various
-:doc:`output commands <Howto_output>`.  This scalar is the gravitational
-potential energy of the particles in the defined field, namely mass \*
-(g dot x) for each particles, where x and mass are the particles
-position and mass, and g is the gravitational field.  The scalar value
-calculated by this fix is "extensive".
+:doc:`output commands <Howto_output>`.  This scalar is the
+gravitational potential energy of the particles in the defined field,
+namely mass \* (g dot x) for each particles, where x and mass are the
+particles position and mass, and g is the gravitational field.  The
+scalar value calculated by this fix is "extensive".

 No parameter of this fix can be used with the *start/stop* keywords of
-the :doc:`run <run>` command.  This fix is not invoked during :doc:`energy minimization <minimize>`.
+the :doc:`run <run>` command.  This fix is not invoked during
+:doc:`energy minimization <minimize>`.

 Restrictions
 """"""""""""
--- a/doc/src/fix_hyper_global.rst
+++ b/doc/src/fix_hyper_global.rst
@ -200,16 +200,20 @@ algorithm.
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-No information about this fix is written to :doc:`binary restart files <restart>`.
+No information about this fix is written to :doc:`binary restart files
+<restart>`.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy of the bias potential to the system's
-potential energy as part of :doc:`thermodynamic output <thermo_style>`.
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the energy of the bias potential to the global
+potential energy of the system as part of :doc:`thermodynamic output
+<thermo_style>`.  The default setting for this fix is :doc:`fix_modify
+energy no <fix_modify>`.

-This fix computes a global scalar and global vector of length 12, which
-can be accessed by various :doc:`output commands <Howto_output>`.  The
-scalar is the magnitude of the bias potential (energy units) applied on
-the current timestep.  The vector stores the following quantities:
+This fix computes a global scalar and global vector of length 12,
+which can be accessed by various :doc:`output commands
+<Howto_output>`.  The scalar is the magnitude of the bias potential
+(energy units) applied on the current timestep.  The vector stores the
+following quantities:

 * 1 = boost factor on this step (unitless)
 * 2 = max strain :math:`E_{ij}` of any bond on this step (absolute value, unitless)
@ -253,7 +257,8 @@ The scalar and vector values calculated by this fix are all
 "intensive".

 No parameter of this fix can be used with the *start/stop* keywords of
-the :doc:`run <run>` command.  This fix is not invoked during :doc:`energy minimization <minimize>`.
+the :doc:`run <run>` command.  This fix is not invoked during
+:doc:`energy minimization <minimize>`.

 Restrictions
 """"""""""""
--- a/doc/src/fix_hyper_local.rst
+++ b/doc/src/fix_hyper_local.rst
@ -370,15 +370,17 @@ Restart, fix_modify, output, run start/stop, minimize info

 No information about this fix is written to :doc:`binary restart files <restart>`.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy of the bias potential to the system's potential
-energy as part of :doc:`thermodynamic output <thermo_style>`.
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the energy of the bias potential to the global
+potential energy of the system as part of :doc:`thermodynamic output
+<thermo_style>`.  The default setting for this fix is :doc:`fix_modify
+energy no <fix_modify>`.

 This fix computes a global scalar and global vector of length 28,
-which can be accessed by various :doc:`output commands <Howto_output>`.
-The scalar is the magnitude of the bias potential (energy units)
-applied on the current timestep, summed over all biased bonds.  The
-vector stores the following quantities:
+which can be accessed by various :doc:`output commands
+<Howto_output>`.  The scalar is the magnitude of the bias potential
+(energy units) applied on the current timestep, summed over all biased
+bonds.  The vector stores the following quantities:

 * 1 = average boost for all bonds on this step (unitless)
 * 2 = # of biased bonds on this step
@ -510,8 +512,8 @@ Value 27 computes the average boost for biased bonds only on this step.
 Value 28 is the count of bonds with an absolute value of strain >= q
 on this step.

-The scalar and vector values calculated by this fix are all
-"intensive".
+The scalar value is an "extensive" quantity since it grows with the
+system size; the vector values are all "intensive".

 This fix also computes a local vector of length the number of bonds
 currently in the system.  The value for each bond is its :math:`C_{ij}`
@ -524,7 +526,8 @@ close to 1.0, which indicates a good choice of :math:`V^{max}`.
 The local values calculated by this fix are unitless.

 No parameter of this fix can be used with the *start/stop* keywords of
-the :doc:`run <run>` command.  This fix is not invoked during :doc:`energy minimization <minimize>`.
+the :doc:`run <run>` command.  This fix is not invoked during
+:doc:`energy minimization <minimize>`.

 Restrictions
 """"""""""""
--- a/doc/src/fix_indent.rst
+++ b/doc/src/fix_indent.rst
@ -179,20 +179,25 @@ contains *xlat*\ , *ylat*\ , *zlat* keywords of the
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-No information about this fix is written to :doc:`binary restart files <restart>`.
+No information about this fix is written to :doc:`binary restart files
+<restart>`.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy of interaction between atoms and the indenter to
-the system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.  The energy of each particle interacting
-with the indenter is K/3 (r - R)\^3.
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the energy of interaction between atoms and the
+indenter to the global potential energy of the system as part of
+:doc:`thermodynamic output <thermo_style>`.  The default setting for
+this fix is :doc:`fix_modify energy no <fix_modify>`. The energy of
+each particle interacting with the indenter is K/3 (r - R)\^3.

 The :doc:`fix_modify <fix_modify>` *respa* option is supported by this
-fix. This allows to set at which level of the :doc:`r-RESPA <run_style>`
-integrator the fix is adding its forces. Default is the outermost level.
+fix. This allows to set at which level of the :doc:`r-RESPA
+<run_style>` integrator the fix is adding its forces. Default is the
+outermost level.

 This fix computes a global scalar energy and a global 3-vector of
-forces (on the indenter), which can be accessed by various :doc:`output commands <Howto_output>`.  The scalar and vector values calculated
-by this fix are "extensive".
+forces (on the indenter), which can be accessed by various
+:doc:`output commands <Howto_output>`.  The scalar and vector values
+calculated by this fix are "extensive".

 The forces due to this fix are imposed during an energy minimization,
 invoked by the :doc:`minimize <minimize>` command.  Note that if you
@ -204,10 +209,10 @@ check if you have done this.

 .. note::

-   If you want the atom/indenter interaction energy to be included
-   in the total potential energy of the system (the quantity being
-   minimized), you must enable the :doc:`fix_modify <fix_modify>` *energy*
-   option for this fix.
+   If you want the atom/indenter interaction energy to be included in
+   the total potential energy of the system (the quantity being
+   minimized), you must enable the :doc:`fix_modify <fix_modify>`
+   *energy* option for this fix.

 Restrictions
 """"""""""""
--- a/doc/src/fix_langevin.rst
+++ b/doc/src/fix_langevin.rst
@ -230,7 +230,7 @@ conservation.

 .. note::

-   this accumulated energy does NOT include kinetic energy removed
+   This accumulated energy does NOT include kinetic energy removed
   by the *zero* flag. LAMMPS will print a warning when both options are
   active.

@ -244,7 +244,8 @@ to zero by subtracting off an equal part of it from each atom in the
 group.  As a result, the center-of-mass of a system with zero initial
 momentum will not drift over time.

-The keyword *gjf* can be used to run the :ref:`Gronbech-Jensen/Farago <Gronbech-Jensen>` time-discretization of the Langevin model.  As
+The keyword *gjf* can be used to run the :ref:`Gronbech-Jensen/Farago
+<Gronbech-Jensen>` time-discretization of the Langevin model.  As
 described in the papers cited below, the purpose of this method is to
 enable longer timesteps to be used (up to the numerical stability
 limit of the integrator), while still producing the correct Boltzmann
@ -252,19 +253,20 @@ distribution of atom positions.

 The current implementation provides the user with the option to output
 the velocity in one of two forms: *vfull* or *vhalf*\ , which replaces
-the outdated option *yes*\ . The *gjf* option *vfull* outputs the on-site
-velocity given in :ref:`Gronbech-Jensen/Farago <Gronbech-Jensen>`; this velocity
-is shown to be systematically lower than the target temperature by a small
-amount, which grows quadratically with the timestep.
-The *gjf* option *vhalf* outputs the 2GJ half-step velocity given in
-:ref:`Gronbech Jensen/Gronbech-Jensen <2Gronbech-Jensen>`; for linear systems,
-this velocity is shown to not have any statistical errors for any stable time step.
-An overview of statistically correct Boltzmann and Maxwell-Boltzmann
-sampling of true on-site and true half-step velocities is given in
-:ref:`Gronbech-Jensen <1Gronbech-Jensen>`.
-Regardless of the choice of output velocity, the sampling of the configurational
-distribution of atom positions is the same, and linearly consistent with the
-target temperature.
+the outdated option *yes*\ . The *gjf* option *vfull* outputs the
+on-site velocity given in :ref:`Gronbech-Jensen/Farago
+<Gronbech-Jensen>`; this velocity is shown to be systematically lower
+than the target temperature by a small amount, which grows
+quadratically with the timestep.  The *gjf* option *vhalf* outputs the
+2GJ half-step velocity given in :ref:`Gronbech Jensen/Gronbech-Jensen
+<2Gronbech-Jensen>`; for linear systems, this velocity is shown to not
+have any statistical errors for any stable time step.  An overview of
+statistically correct Boltzmann and Maxwell-Boltzmann sampling of true
+on-site and true half-step velocities is given in
+:ref:`Gronbech-Jensen <1Gronbech-Jensen>`.  Regardless of the choice
+of output velocity, the sampling of the configurational distribution
+of atom positions is the same, and linearly consistent with the target
+temperature.

 ----------

@ -287,16 +289,18 @@ you have defined to this fix which will be used in its thermostatting
 procedure, as described above.  For consistency, the group used by
 this fix and by the compute should be the same.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy change induced by Langevin thermostatting to the
-system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.  Note that use of this option requires
-setting the *tally* keyword to *yes*\ .
+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*\ , but only if the *tally* keyword to set to
+*yes*\ .  See the :doc:`thermo_style <thermo_style>` doc page for
+details.

 This fix computes a global scalar which can be accessed by various
-:doc:`output commands <Howto_output>`.  The scalar is the cumulative
-energy change due to this fix.  The scalar value calculated by this
-fix is "extensive".  Note that calculation of this quantity requires
-setting the *tally* keyword to *yes*\ .
+:doc:`output commands <Howto_output>`.  The scalar is the same
+cumulative energy change due to this fix described in the previous
+paragraph.  The scalar value calculated by this fix is "extensive".
+Note that calculation of this quantity also requires setting the
+*tally* keyword to *yes*\ .

 This fix can ramp its target temperature over multiple runs, using the
 *start* and *stop* keywords of the :doc:`run <run>` command.  See the
--- a/doc/src/fix_langevin_eff.rst
+++ b/doc/src/fix_langevin_eff.rst
@ -81,16 +81,18 @@ you have defined to this fix which will be used in its thermostatting
 procedure, as described above.  For consistency, the group used by
 this fix and by the compute should be the same.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy change induced by Langevin thermostatting to the
-system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.  Note that use of this option requires
-setting the *tally* keyword to *yes*\ .
+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*\ , but only if the *tally* keyword to set to
+*yes*\ .  See the :doc:`thermo_style <thermo_style>` doc page for
+details.

 This fix computes a global scalar which can be accessed by various
-:doc:`output commands <Howto_output>`.  The scalar is the cumulative
-energy change due to this fix.  The scalar value calculated by this
-fix is "extensive".  Note that calculation of this quantity requires
-setting the *tally* keyword to *yes*\ .
+:doc:`output commands <Howto_output>`.  The scalar is the same
+cumulative energy change due to this fix described in the previous
+paragraph.  The scalar value calculated by this fix is "extensive".
+Note that calculation of this quantity also requires setting the
+*tally* keyword to *yes*\ .

 This fix can ramp its target temperature over multiple runs, using the
 *start* and *stop* keywords of the :doc:`run <run>` command.  See the
--- a/doc/src/fix_latte.rst
+++ b/doc/src/fix_latte.rst
@ -107,16 +107,27 @@ larger system sizes and longer time scales
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-No information about this fix is written to :doc:`binary restart files <restart>`.
+No information about this fix is written to :doc:`binary restart files
+<restart>`.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the potential energy computed by LATTE to the system's
-potential energy as part of :doc:`thermodynamic output <thermo_style>`.
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the potential energy computed by LATTE to the global
+potential energy of the system as part of :doc:`thermodynamic output
+<thermo_style>`.  The default setting for this fix is :doc:`fix_modify
+energy yes <fix_modify>`.

-The :doc:`fix_modify <fix_modify>` *virial* option is supported by this
-fix to add the LATTE DFTB contribution to the system's virial as part
-of :doc:`thermodynamic output <thermo_style>`.  The default is *virial
-yes*
+The :doc:`fix_modify <fix_modify>` *virial* option is supported by
+this fix to add the contribution compute by LATTE to the global
+pressure of the system via the :doc:`compute pressure
+<compute_pressure>` command.  This can be accessed by
+:doc:`thermodynamic output <thermo_style>`.  The default setting for
+this fix is :doc:`fix_modify virial yes <fix_modify>`.
+
+The :doc:`fix_modify <fix_modify>` *virial* option is supported by
+this fix to add the contribution computed by LATTE to the global
+pressure of the system as part of :doc:`thermodynamic output
+<thermo_style>`.  The default setting for this fix is :doc:`fix_modify
+virial yes <fix_modify>`.

 This fix computes a global scalar which can be accessed by various
 :doc:`output commands <Howto_output>`.  The scalar is the potential
@ -127,20 +138,22 @@ No parameter of this fix can be used with the *start/stop* keywords of
 the :doc:`run <run>` command.

 The DFTB forces computed by LATTE via this fix are imposed during an
-energy minimization, invoked by the :doc:`minimize <minimize>` command.
+energy minimization, invoked by the :doc:`minimize <minimize>`
+command.

 .. note::

   If you want the potential energy associated with the DFTB
   forces to be included in the total potential energy of the system (the
-   quantity being minimized), you MUST enable the
+   quantity being minimized), you MUST not disable the
   :doc:`fix_modify <fix_modify>` *energy* option for this fix.

 Restrictions
 """"""""""""

 This fix is part of the LATTE package.  It is only enabled if LAMMPS
-was built with that package.  See the :doc:`Build package <Build_package>` doc page for more info.
+was built with that package.  See the :doc:`Build package
+<Build_package>` doc page for more info.

 You must use metal units, as set by the :doc:`units <units>` command to
 use this fix.
--- a/doc/src/fix_lb_rigid_pc_sphere.rst
+++ b/doc/src/fix_lb_rigid_pc_sphere.rst
@ -50,7 +50,8 @@ This fix is based on the :doc:`fix rigid <fix_rigid>` command, and was
 created to be used in place of that fix, to integrate the equations of
 motion of spherical rigid bodies when a lattice-Boltzmann fluid is
 present with a user-specified value of the force-coupling constant.
-The fix uses the integration algorithm described in :ref:`Mackay et al. <Mackay>` to update the positions, velocities, and orientations of
+The fix uses the integration algorithm described in :ref:`Mackay et
+al. <Mackay>` to update the positions, velocities, and orientations of
 a set of spherical rigid bodies experiencing velocity dependent
 hydrodynamic forces.  The spherical bodies are assumed to rotate as
 solid, uniform density spheres, with moments of inertia calculated
@ -88,9 +89,18 @@ Restart, fix_modify, output, run start/stop, minimize info
 No information about the *rigid* and *rigid/nve* fixes are written to
 :doc:`binary restart files <restart>`.

+The :doc:`fix_modify <fix_modify>` *virial* option is supported by
+this fix to add the contribution due to the added forces on atoms to
+both the global pressure and per-atom stress of the system via the
+:doc:`compute pressure <compute_pressure>` and :doc:`compute
+stress/atom <compute_stress_atom>` commands.  The former can be
+accessed by :doc:`thermodynamic output <thermo_style>`.  The default
+setting for this fix is :doc:`fix_modify virial yes <fix_modify>`.
+
 Similar to the :doc:`fix rigid <fix_rigid>` command: The rigid fix
-computes a global scalar which can be accessed by various :doc:`output commands <Howto_output>`.  The scalar value calculated by these
-fixes is "intensive".  The scalar is the current temperature of the
+computes a global scalar which can be accessed by various :doc:`output
+commands <Howto_output>`.  The scalar value calculated by these fixes
+is "intensive".  The scalar is the current temperature of the
 collection of rigid bodies.  This is averaged over all rigid bodies
 and their translational and rotational degrees of freedom.  The
 translational energy of a rigid body is 1/2 m v\^2, where m = total
@ -130,7 +140,8 @@ Restrictions
 """"""""""""

 This fix is part of the USER-LB package.  It is only enabled if LAMMPS
-was built with that package.  See the :doc:`Build package <Build_package>` doc page for more info.
+was built with that package.  See the :doc:`Build package
+<Build_package>` doc page for more info.

 Can only be used if a lattice-Boltzmann fluid has been created via the
 :doc:`fix lb/fluid <fix_lb_fluid>` command, and must come after this
--- a/doc/src/fix_modify.rst
+++ b/doc/src/fix_modify.rst
@ -59,43 +59,58 @@ define their own compute by default, as described in their
 documentation.  Thus this option allows the user to override the
 default method for computing P.

-The *energy* keyword can be used with fixes that support it.
-*energy yes* adds a contribution to the potential energy of the
-system. The fix's global and per-atom
-energy is included in the calculation performed by the :doc:`compute pe <compute_pe>` or :doc:`compute pe/atom <compute_pe_atom>`
-commands.  See the :doc:`thermo_style <thermo_style>` command for info
-on how potential energy is output.  For fixes that tally a global
-energy, it can be printed by using the keyword f_ID in the
-thermo_style custom command, where ID is the fix-ID of the appropriate
-fix.
+The *energy* keyword can be used with fixes that support it, which is
+explained at the bottom of their doc page.  *Energy yes* will add a
+contribution to the potential energy of the system.  More
+specifically, the fix's global or per-atom energy is included in the
+calculation performed by the :doc:`compute pe <compute_pe>` or
+:doc:`compute pe/atom <compute_pe_atom>` commands.  The former is what
+is used the :doc:`thermo_style <thermo_style>` command for output of
+any quantity that includes the global potential energy of the system.
+Note that the :doc:`compute pe <compute_pe>` and :doc:`compute pe/atom
+<compute_pe_atom>` commands also have an option to include or exclude
+the contribution from fixes.  For fixes that tally a global energy, it
+can also be printed with thermodynamic output by using the keyword
+f_ID in the thermo_style custom command, where ID is the fix-ID of the
+appropriate fix.

 .. note::

-   You must also specify the *energy yes* setting for a fix if you
-   are using it when performing an :doc:`energy minimization <minimize>`
-   and if you want the energy and forces it produces to be part of the
-   optimization criteria.
-
-The *virial* keyword can be used with fixes that support it.
-*virial yes* adds a contribution to the virial of the
-system. The fix's global and per-atom
-virial is included in the calculation performed by the :doc:`compute pressure <compute_pressure>` or
-:doc:`compute stress/atom <compute_stress_atom>`
-commands.  See the :doc:`thermo_style <thermo_style>` command for info
-on how pressure is output.
+   If you are performing an :doc:`energy minimization <minimize>` with
+   one of these fixes and want the energy and forces it produces to be
+   part of the optimization criteria, you must specify the *energy
+   yes* setting.

 .. note::

-   You must specify the *virial yes* setting for a fix if you
-   are doing :doc:`box relaxation <fix_box_relax>` and
-   if you want virial contribution of the fix to be part of the
-   relaxation criteria, although this seems unlikely.
+   For most fixes that support the *energy* keyword, the default
+   setting is *no*.  For a few it is *yes*, when a user would expect
+   that to be the case.  The doc page of each fix gives the default.
+
+The *virial* keyword can be used with fixes that support it, which is
+explained at the bottom of their doc page.  *Virial yes* will add a
+contribution to the virial of the system.  More specifically, the
+fix's global or per-atom virial is included in the calculation
+performed by the :doc:`compute pressure <compute_pressure>` or
+:doc:`compute stress/atom <compute_stress_atom>` commands.  The former
+is what is used the :doc:`thermo_style <thermo_style>` command for
+output of any quantity that includes the global pressure of the
+system.  Note that the :doc:`compute pressure <compute_pressure>` and
+:doc:`compute stress/atom <compute_stress_atom>` commands also have an
+option to include or exclude the contribution from fixes.

 .. note::

-   This option is only supported by fixes that explicitly say
-   so. For some of these (e.g. the :doc:`fix shake <fix_shake>` command)
-   the default setting is *virial yes*\ , for others it is *virial no*\ .
+   If you are performing an :doc:`energy minimization <minimize>` with
+   :doc:`box relaxation <fix_box_relax>` and one of these fixes and
+   want the virial contribution of the fix to be part of the
+   optimization criteria, you must specify the *virial yes* setting.
+
+.. note::
+
+   For most fixes that support the *virial* keyword, the default
+   setting is *no*.  For a few it is *yes*, when a user would expect
+   that to be the case.  The doc page of each fix gives the default.

 For fixes that set or modify forces, it may be possible to select at
 which :doc:`r-RESPA <run_style>` level the fix operates via the *respa*
@ -112,13 +127,15 @@ The *dynamic/dof* keyword determines whether the number of atoms N in
 the fix group and their associated degrees of freedom are re-computed
 each time a temperature is computed.  Only fix styles that calculate
 their own internal temperature use this option.  Currently this is
-only the :doc:`fix rigid/nvt/small <fix_rigid>` and :doc:`fix rigid/npt/small <fix_rigid>` commands for the purpose of
+only the :doc:`fix rigid/nvt/small <fix_rigid>` and :doc:`fix
+rigid/npt/small <fix_rigid>` commands for the purpose of
 thermostatting rigid body translation and rotation.  By default, N and
 their DOF are assumed to be constant.  If you are adding atoms or
-molecules to the system (see the :doc:`fix pour <fix_pour>`, :doc:`fix deposit <fix_deposit>`, and :doc:`fix gcmc <fix_gcmc>` commands) or
+molecules to the system (see the :doc:`fix pour <fix_pour>`, :doc:`fix
+deposit <fix_deposit>`, and :doc:`fix gcmc <fix_gcmc>` commands) or
 expect atoms or molecules to be lost (e.g. due to exiting the
-simulation box or via :doc:`fix evaporate <fix_evaporate>`), then
-this option should be used to insure the temperature is correctly
+simulation box or via :doc:`fix evaporate <fix_evaporate>`), then this
+option should be used to insure the temperature is correctly
 normalized.

 .. note::
--- a/doc/src/fix_msst.rst
+++ b/doc/src/fix_msst.rst
@ -131,20 +131,29 @@ command for info on how to re-specify a fix in an input script that
 reads a restart file, so that the operation of the fix continues in an
 uninterrupted fashion.

+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*.  See the :doc:`thermo_style <thermo_style>`
+doc page for details.
+
+This fix computes a global scalar which can be accessed by various
+:doc:`output commands <Howto_output>`.  The scalar is the same
+cumulative energy change due to this fix described in the previous
+paragraph.  The scalar value calculated by this fix is "extensive".
+
 The progress of the MSST can be monitored by printing the global
 scalar and global vector quantities computed by the fix.

-The scalar is the cumulative energy change due to the fix. This is
-also the energy added to the potential energy by the
-:doc:`fix_modify <fix_modify>` *energy* command.  With this command, the
-thermo keyword *etotal* prints the conserved quantity of the MSST
-dynamic equations. This can be used to test if the MD timestep is
-sufficiently small for accurate integration of the dynamic
-equations. See also :doc:`thermo_style <thermo_style>` command.
+As mentioned above, the scalar is the cumulative energy change due to
+the fix.  By monitoring the thermodynamic *econserve* output, this can
+be used to test if the MD timestep is sufficiently small for accurate
+integration of the dynamic equations.

-The global vector contains four values in this order:
+The global vector contains four values in the following order.  The
+vector values output by this fix are "intensive".

-[\ *dhugoniot*\ , *drayleigh*\ , *lagrangian_speed*, *lagrangian_position*]
+[\ *dhugoniot*\ , *drayleigh*\ , *lagrangian_speed*,
+*lagrangian_position*]

 1. *dhugoniot* is the departure from the Hugoniot (temperature units).
 2. *drayleigh* is the departure from the Rayleigh line (pressure units).
@ -157,17 +166,11 @@ headers, the following LAMMPS commands are suggested:
 .. code-block:: LAMMPS

   fix              msst all msst z
-   fix_modify       msst energy yes
   variable dhug    equal f_msst[1]
   variable dray    equal f_msst[2]
   variable lgr_vel equal f_msst[3]
   variable lgr_pos equal f_msst[4]
-   thermo_style     custom step temp ke pe lz pzz etotal v_dhug v_dray v_lgr_vel v_lgr_pos f_msst
-
-These fixes compute a global scalar and a global vector of 4
-quantities, which can be accessed by various :doc:`output commands
-<Howto_output>`.  The scalar values calculated by this fix are
-"extensive"; the vector values are "intensive".
+   thermo_style     custom step temp ke pe lz pzz econserve v_dhug v_dray v_lgr_vel v_lgr_pos f_msst

 Restrictions
 """"""""""""
--- a/doc/src/fix_nh.rst
+++ b/doc/src/fix_nh.rst
@ -594,17 +594,20 @@ compute temperature on a subset of atoms.
   specified by the *press* keyword will be unaffected by the *temp*
   setting.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by these
-fixes to add the energy change induced by Nose/Hoover thermostatting
-and barostatting to the system's potential energy as part of
-:doc:`thermodynamic output <thermo_style>`.
+The cumulative energy change in the system imposed by these fixes, via
+either thermostatting and/or barostatting, is included in the
+:doc:`thermodynamic output <thermo_style>` keywords *ecouple* and
+*econserve*.  See the :doc:`thermo_style <thermo_style>` doc page for
+details.

-These fixes compute a global scalar and a global vector of quantities,
-which can be accessed by various :doc:`output commands <Howto_output>`.
-The scalar value calculated by these fixes is "extensive"; the vector
-values are "intensive".
+These fixes compute a global scalar which can be accessed by various
+:doc:`output commands <Howto_output>`.  The scalar is the same
+cumulative energy change due to this fix described in the previous
+paragraph.  The scalar value calculated by this fix is "extensive".

-The scalar is the cumulative energy change due to the fix.
+These fixes compute also compute a global vector of quantities, which
+can be accessed by various :doc:`output commands <Howto_output>`.  The
+vector values are "intensive".

 The vector stores internal Nose/Hoover thermostat and barostat
 variables.  The number and meaning of the vector values depends on
--- a/doc/src/fix_nph_asphere.rst
+++ b/doc/src/fix_nph_asphere.rst
@ -87,7 +87,8 @@ It also means that changing attributes of *thermo_temp* or
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-This fix writes the state of the Nose/Hoover barostat to :doc:`binary restart files <restart>`.  See the :doc:`read_restart <read_restart>`
+This fix writes the state of the Nose/Hoover barostat to :doc:`binary
+restart files <restart>`.  See the :doc:`read_restart <read_restart>`
 command for info on how to re-specify a fix in an input script that
 reads a restart file, so that the operation of the fix continues in an
 uninterrupted fashion.
@ -101,9 +102,10 @@ consistent with the virial term computed using all atoms for the
 pressure.  LAMMPS will warn you if you choose to compute temperature
 on a subset of atoms.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy change induced by Nose/Hoover barostatting to
-the system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.
+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*.  See the :doc:`thermo_style <thermo_style>`
+doc page for details.

 This fix computes the same global scalar and global vector of
 quantities as does the :doc:`fix nph <fix_nh>` command.
--- a/doc/src/fix_nph_body.rst
+++ b/doc/src/fix_nph_body.rst
@ -84,7 +84,8 @@ It also means that changing attributes of *thermo_temp* or
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-This fix writes the state of the Nose/Hoover barostat to :doc:`binary restart files <restart>`.  See the :doc:`read_restart <read_restart>`
+This fix writes the state of the Nose/Hoover barostat to :doc:`binary
+restart files <restart>`.  See the :doc:`read_restart <read_restart>`
 command for info on how to re-specify a fix in an input script that
 reads a restart file, so that the operation of the fix continues in an
 uninterrupted fashion.
@ -98,9 +99,10 @@ consistent with the virial term computed using all atoms for the
 pressure.  LAMMPS will warn you if you choose to compute temperature
 on a subset of atoms.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy change induced by Nose/Hoover barostatting to
-the system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.
+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*.  See the :doc:`thermo_style <thermo_style>`
+doc page for details.

 This fix computes the same global scalar and global vector of
 quantities as does the :doc:`fix nph <fix_nh>` command.
--- a/doc/src/fix_nph_sphere.rst
+++ b/doc/src/fix_nph_sphere.rst
@ -100,7 +100,8 @@ It also means that changing attributes of *thermo_temp* or
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-This fix writes the state of the Nose/Hoover barostat to :doc:`binary restart files <restart>`.  See the :doc:`read_restart <read_restart>`
+This fix writes the state of the Nose/Hoover barostat to :doc:`binary
+restart files <restart>`.  See the :doc:`read_restart <read_restart>`
 command for info on how to re-specify a fix in an input script that
 reads a restart file, so that the operation of the fix continues in an
 uninterrupted fashion.
@ -114,9 +115,10 @@ consistent with the virial term computed using all atoms for the
 pressure.  LAMMPS will warn you if you choose to compute temperature
 on a subset of atoms.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy change induced by Nose/Hoover barostatting to
-the system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.
+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*.  See the :doc:`thermo_style <thermo_style>`
+doc page for details.

 This fix computes the same global scalar and global vector of
 quantities as does the :doc:`fix nph <fix_nh>` command.
--- a/doc/src/fix_nphug.rst
+++ b/doc/src/fix_nphug.rst
@ -168,43 +168,47 @@ specified, then the instantaneous value in the system at the start of
 the simulation is used.

 The :doc:`fix_modify <fix_modify>` *temp* and *press* options are
-supported by these fixes.  You can use them to assign a
-:doc:`compute <compute>` you have defined to this fix which will be used
-in its thermostatting or barostatting procedure, as described above.
-If you do this, note that the kinetic energy derived from the compute
+supported by this fix.  You can use them to assign a :doc:`compute
+<compute>` you have defined to this fix which will be used in its
+thermostatting or barostatting procedure, as described above.  If you
+do this, note that the kinetic energy derived from the compute
 temperature should be consistent with the virial term computed using
 all atoms for the pressure.  LAMMPS will warn you if you choose to
 compute temperature on a subset of atoms.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by these
-fixes to add the energy change induced by Nose/Hoover thermostatting
-and barostatting to the system's potential energy as part of
-:doc:`thermodynamic output <thermo_style>`. Either way, this energy is \*not\*
-included in the definition of internal energy E when calculating the value
-of Delta in the above equation.
+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*.  See the :doc:`thermo_style <thermo_style>`
+doc page for details.  Note that this energy is \*not\* included in
+the definition of internal energy E when calculating the value of
+Delta in the above equation.

-These fixes compute a global scalar and a global vector of quantities,
-which can be accessed by various :doc:`output commands <Howto_output>`.
-The scalar value calculated by these fixes is "extensive"; the vector
-values are "intensive".
+This fix computes a global scalar which can be accessed by various
+:doc:`output commands <Howto_output>`.  The scalar is the same
+cumulative energy change due to this fix described in the previous
+paragraph.  The scalar value calculated by this fix is "extensive".

-The scalar is the cumulative energy change due to the fix.
+This fix also computes a global vector of quantities, which can be
+accessed by various :doc:`output commands <Howto_output>`.  The scalar
+The vector values are "intensive".

-The vector stores three quantities unique to this fix (:math:`\Delta`, Us, and up),
-followed by all the internal Nose/Hoover thermostat and barostat
-variables defined for :doc:`fix npt <fix_nh>`. Delta is the deviation
-of the temperature from the target temperature, given by the above equation.
-Us and up are the shock and particle velocity corresponding to a steady
-shock calculated from the RH conditions. They have units of distance/time.
+The vector stores three quantities unique to this fix (:math:`\Delta`,
+Us, and up), followed by all the internal Nose/Hoover thermostat and
+barostat variables defined for :doc:`fix npt <fix_nh>`. Delta is the
+deviation of the temperature from the target temperature, given by the
+above equation.  Us and up are the shock and particle velocity
+corresponding to a steady shock calculated from the RH
+conditions. They have units of distance/time.

 Restrictions
 """"""""""""

 This fix style is part of the SHOCK package.  It is only enabled if
-LAMMPS was built with that package. See the :doc:`Build package <Build_package>` doc page for more info.
+LAMMPS was built with that package. See the :doc:`Build package
+<Build_package>` doc page for more info.

-All the usual restrictions for :doc:`fix npt <fix_nh>` apply,
-plus the additional ones mentioned above.
+All the usual restrictions for :doc:`fix npt <fix_nh>` apply, plus the
+additional ones mentioned above.

 Related commands
 """"""""""""""""
--- a/doc/src/fix_npt_asphere.rst
+++ b/doc/src/fix_npt_asphere.rst
@ -124,10 +124,10 @@ consistent with the virial term computed using all atoms for the
 pressure.  LAMMPS will warn you if you choose to compute temperature
 on a subset of atoms.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy change induced by Nose/Hoover thermostatting and
-barostatting to the system's potential energy as part of
-:doc:`thermodynamic output <thermo_style>`.
+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*.  See the :doc:`thermo_style <thermo_style>`
+doc page for details.

 This fix computes the same global scalar and global vector of
 quantities as does the :doc:`fix npt <fix_nh>` command.
--- a/doc/src/fix_npt_body.rst
+++ b/doc/src/fix_npt_body.rst
@ -45,7 +45,8 @@ can also have a bias velocity removed from them before thermostatting
 takes place; see the description below.

 Additional parameters affecting the thermostat and barostat are
-specified by keywords and values documented with the :doc:`fix npt <fix_nh>` command.  See, for example, discussion of the *temp*\ ,
+specified by keywords and values documented with the :doc:`fix npt
+<fix_nh>` command.  See, for example, discussion of the *temp*\ ,
 *iso*\ , *aniso*\ , and *dilate* keywords.

 The particles in the fix group are the only ones whose velocities and
@ -121,10 +122,10 @@ consistent with the virial term computed using all atoms for the
 pressure.  LAMMPS will warn you if you choose to compute temperature
 on a subset of atoms.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy change induced by Nose/Hoover thermostatting and
-barostatting to the system's potential energy as part of
-:doc:`thermodynamic output <thermo_style>`.
+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*.  See the :doc:`thermo_style <thermo_style>`
+doc page for details.

 This fix computes the same global scalar and global vector of
 quantities as does the :doc:`fix npt <fix_nh>` command.
--- a/doc/src/fix_npt_cauchy.rst
+++ b/doc/src/fix_npt_cauchy.rst
@ -487,27 +487,31 @@ compute temperature on a subset of atoms.
 .. note::

   If both the *temp* and *press* keywords are used in a single
-   thermo_modify command (or in two separate commands), then the order in
-   which the keywords are specified is important.  Note that a :doc:`pressure compute <compute_pressure>` defines its own temperature compute as
-   an argument when it is specified.  The *temp* keyword will override
-   this (for the pressure compute being used by fix npt), but only if the
-   *temp* keyword comes after the *press* keyword.  If the *temp* keyword
-   comes before the *press* keyword, then the new pressure compute
-   specified by the *press* keyword will be unaffected by the *temp*
-   setting.
+   thermo_modify command (or in two separate commands), then the order
+   in which the keywords are specified is important.  Note that a
+   :doc:`pressure compute <compute_pressure>` defines its own
+   temperature compute as an argument when it is specified.  The
+   *temp* keyword will override this (for the pressure compute being
+   used by fix npt), but only if the *temp* keyword comes after the
+   *press* keyword.  If the *temp* keyword comes before the *press*
+   keyword, then the new pressure compute specified by the *press*
+   keyword will be unaffected by the *temp* setting.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy change induced by Nose/Hoover thermostatting
-and barostatting to the system's potential energy as part of
-:doc:`thermodynamic output <thermo_style>`.
+The cumulative energy change in the system imposed by this fix, due to
+thermostatting and/or barostatting, is included in the
+:doc:`thermodynamic output <thermo_style>` keywords *ecouple* and
+*econserve*.  See the :doc:`thermo_style <thermo_style>` doc page for
+details.

-This fix computes a global scalar and a global vector of quantities,
-which can be accessed by various :doc:`output commands <Howto_output>`.
-The scalar value calculated by this fix is "extensive"; the vector
+This fix computes a global scalar which can be accessed by various
+:doc:`output commands <Howto_output>`.  The scalar is the same
+cumulative energy change due to this fix described in the previous
+paragraph.  The scalar value calculated by this fix is "extensive".
+
+This fix also computes a global vector of quantities, which can be
+accessed by various :doc:`output commands <Howto_output>`.  Rhe vector
 values are "intensive".

-The scalar is the cumulative energy change due to the fix.
-
 The vector stores internal Nose/Hoover thermostat and barostat
 variables.  The number and meaning of the vector values depends on
 which fix is used and the settings for keywords *tchain* and *pchain*\ ,
--- a/doc/src/fix_npt_sphere.rst
+++ b/doc/src/fix_npt_sphere.rst
@ -93,13 +93,14 @@ IDs of the new computes are the fix-ID + underscore + "temp" or fix_ID
 since pressure is computed for the entire system.

 Note that these are NOT the computes used by thermodynamic output (see
-the :doc:`thermo_style <thermo_style>` command) with ID = *thermo_temp*
-and *thermo_press*.  This means you can change the attributes of this
-fix's temperature or pressure via the
-:doc:`compute_modify <compute_modify>` command or print this temperature
-or pressure during thermodynamic output via the :doc:`thermo_style custom <thermo_style>` command using the appropriate compute-ID.
-It also means that changing attributes of *thermo_temp* or
-*thermo_press* will have no effect on this fix.
+the :doc:`thermo_style <thermo_style>` command) with ID =
+*thermo_temp* and *thermo_press*.  This means you can change the
+attributes of this fix's temperature or pressure via the
+:doc:`compute_modify <compute_modify>` command or print this
+temperature or pressure during thermodynamic output via the
+:doc:`thermo_style custom <thermo_style>` command using the
+appropriate compute-ID.  It also means that changing attributes of
+*thermo_temp* or *thermo_press* will have no effect on this fix.

 Like other fixes that perform thermostatting, this fix can be used
 with :doc:`compute commands <compute>` that calculate a temperature
@ -129,18 +130,18 @@ a fix in an input script that reads a restart file, so that the
 operation of the fix continues in an uninterrupted fashion.

 The :doc:`fix_modify <fix_modify>` *temp* and *press* options are
-supported by this fix.  You can use them to assign a
-:doc:`compute <compute>` you have defined to this fix which will be used
-in its thermostatting or barostatting procedure.  If you do this, note
-that the kinetic energy derived from the compute temperature should be
+supported by this fix.  You can use them to assign a :doc:`compute
+<compute>` you have defined to this fix which will be used in its
+thermostatting or barostatting procedure.  If you do this, note that
+the kinetic energy derived from the compute temperature should be
 consistent with the virial term computed using all atoms for the
 pressure.  LAMMPS will warn you if you choose to compute temperature
 on a subset of atoms.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy change induced by Nose/Hoover thermostatting and
-barostatting to the system's potential energy as part of
-:doc:`thermodynamic output <thermo_style>`.
+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*.  See the :doc:`thermo_style <thermo_style>`
+doc page for details.

 This fix computes the same global scalar and global vector of
 quantities as does the :doc:`fix npt <fix_nh>` command.
--- a/doc/src/fix_nvt_asphere.rst
+++ b/doc/src/fix_nvt_asphere.rst
@ -92,19 +92,21 @@ thermal degrees of freedom, and the bias is added back in.
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-This fix writes the state of the Nose/Hoover thermostat to :doc:`binary restart files <restart>`.  See the :doc:`read_restart <read_restart>`
-command for info on how to re-specify a fix in an input script that
-reads a restart file, so that the operation of the fix continues in an
-uninterrupted fashion.
+This fix writes the state of the Nose/Hoover thermostat to
+:doc:`binary restart files <restart>`.  See the :doc:`read_restart
+<read_restart>` command for info on how to re-specify a fix in an
+input script that reads a restart file, so that the operation of the
+fix continues in an uninterrupted fashion.

 The :doc:`fix_modify <fix_modify>` *temp* option is supported by this
 fix.  You can use it to assign a :doc:`compute <compute>` you have
 defined to this fix which will be used in its thermostatting
 procedure.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy change induced by Nose/Hoover thermostatting to
-the system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.
+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*.  See the :doc:`thermo_style <thermo_style>`
+doc page for details.

 This fix computes the same global scalar and global vector of
 quantities as does the :doc:`fix nvt <fix_nh>` command.
--- a/doc/src/fix_nvt_body.rst
+++ b/doc/src/fix_nvt_body.rst
@ -89,19 +89,21 @@ thermal degrees of freedom, and the bias is added back in.
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-This fix writes the state of the Nose/Hoover thermostat to :doc:`binary restart files <restart>`.  See the :doc:`read_restart <read_restart>`
-command for info on how to re-specify a fix in an input script that
-reads a restart file, so that the operation of the fix continues in an
-uninterrupted fashion.
+This fix writes the state of the Nose/Hoover thermostat to
+:doc:`binary restart files <restart>`.  See the :doc:`read_restart
+<read_restart>` command for info on how to re-specify a fix in an
+input script that reads a restart file, so that the operation of the
+fix continues in an uninterrupted fashion.

 The :doc:`fix_modify <fix_modify>` *temp* option is supported by this
 fix.  You can use it to assign a :doc:`compute <compute>` you have
 defined to this fix which will be used in its thermostatting
 procedure.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy change induced by Nose/Hoover thermostatting to
-the system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.
+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*.  See the :doc:`thermo_style <thermo_style>`
+doc page for details.

 This fix computes the same global scalar and global vector of
 quantities as does the :doc:`fix nvt <fix_nh>` command.
--- a/doc/src/fix_nvt_sllod.rst
+++ b/doc/src/fix_nvt_sllod.rst
@ -122,19 +122,21 @@ thermal degrees of freedom, and the bias is added back in.
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-This fix writes the state of the Nose/Hoover thermostat to :doc:`binary restart files <restart>`.  See the :doc:`read_restart <read_restart>`
-command for info on how to re-specify a fix in an input script that
-reads a restart file, so that the operation of the fix continues in an
-uninterrupted fashion.
+This fix writes the state of the Nose/Hoover thermostat to
+:doc:`binary restart files <restart>`.  See the :doc:`read_restart
+<read_restart>` command for info on how to re-specify a fix in an
+input script that reads a restart file, so that the operation of the
+fix continues in an uninterrupted fashion.

 The :doc:`fix_modify <fix_modify>` *temp* option is supported by this
 fix.  You can use it to assign a :doc:`compute <compute>` you have
 defined to this fix which will be used in its thermostatting
 procedure.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy change induced by Nose/Hoover thermostatting to
-the system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.
+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*.  See the :doc:`thermo_style <thermo_style>`
+doc page for details.

 This fix computes the same global scalar and global vector of
 quantities as does the :doc:`fix nvt <fix_nh>` command.
--- a/doc/src/fix_nvt_sllod_eff.rst
+++ b/doc/src/fix_nvt_sllod_eff.rst
@ -41,19 +41,21 @@ velocity.
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-This fix writes the state of the Nose/Hoover thermostat to :doc:`binary restart files <restart>`.  See the :doc:`read_restart <read_restart>`
-command for info on how to re-specify a fix in an input script that
-reads a restart file, so that the operation of the fix continues in an
-uninterrupted fashion.
+This fix writes the state of the Nose/Hoover thermostat to
+:doc:`binary restart files <restart>`.  See the :doc:`read_restart
+<read_restart>` command for info on how to re-specify a fix in an
+input script that reads a restart file, so that the operation of the
+fix continues in an uninterrupted fashion.

 The :doc:`fix_modify <fix_modify>` *temp* option is supported by this
 fix.  You can use it to assign a :doc:`compute <compute>` you have
 defined to this fix which will be used in its thermostatting
 procedure.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy change induced by Nose/Hoover thermostatting to
-the system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.
+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*.  See the :doc:`thermo_style <thermo_style>`
+doc page for details.

 This fix computes the same global scalar and global vector of
 quantities as does the :doc:`fix nvt/eff <fix_nh_eff>` command.
--- a/doc/src/fix_nvt_sphere.rst
+++ b/doc/src/fix_nvt_sphere.rst
@ -106,19 +106,21 @@ thermal degrees of freedom, and the bias is added back in.
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-This fix writes the state of the Nose/Hoover thermostat to :doc:`binary restart files <restart>`.  See the :doc:`read_restart <read_restart>`
-command for info on how to re-specify a fix in an input script that
-reads a restart file, so that the operation of the fix continues in an
-uninterrupted fashion.
+This fix writes the state of the Nose/Hoover thermostat to
+:doc:`binary restart files <restart>`.  See the :doc:`read_restart
+<read_restart>` command for info on how to re-specify a fix in an
+input script that reads a restart file, so that the operation of the
+fix continues in an uninterrupted fashion.

 The :doc:`fix_modify <fix_modify>` *temp* option is supported by this
 fix.  You can use it to assign a :doc:`compute <compute>` you have
 defined to this fix which will be used in its thermostatting
 procedure.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the energy change induced by Nose/Hoover thermostatting to
-the system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.
+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*.  See the :doc:`thermo_style <thermo_style>`
+doc page for details.

 This fix computes the same global scalar and global vector of
 quantities as does the :doc:`fix nvt <fix_nh>` command.
--- a/doc/src/fix_orient.rst
+++ b/doc/src/fix_orient.rst
@ -144,16 +144,20 @@ writing the orientation files is given in :ref:`(Wicaksono2) <Wicaksono2>`
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-No information about this fix is written to :doc:`binary restart files <restart>`.
+No information about this fix is written to :doc:`binary restart files
+<restart>`.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the potential energy of atom interactions with the grain
-boundary driving force to the system's potential energy as part of
-:doc:`thermodynamic output <thermo_style>`.
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the potential energy of atom interactions with the
+grain boundary driving force to the global potential energy of the
+system as part of :doc:`thermodynamic output <thermo_style>`.  The
+default setting for this fix is :doc:`fix_modify energy no
+<fix_modify>`.

-The :doc:`fix_modify <fix_modify>` *respa* option is supported by these
-fixes. This allows to set at which level of the :doc:`r-RESPA <run_style>`
-integrator a fix is adding its forces. Default is the outermost level.
+The :doc:`fix_modify <fix_modify>` *respa* option is supported by
+these fixes. This allows to set at which level of the :doc:`r-RESPA
+<run_style>` integrator a fix is adding its forces. Default is the
+outermost level.

 This fix calculates a global scalar which can be accessed by various
 :doc:`output commands <Howto_output>`.  The scalar is the potential
@ -166,13 +170,16 @@ order parameter Xi and normalized order parameter (0 to 1) for each
 atom.  The per-atom values can be accessed on any timestep.

 No parameter of this fix can be used with the *start/stop* keywords of
-the :doc:`run <run>` command.  This fix is not invoked during :doc:`energy minimization <minimize>`.
+the :doc:`run <run>` command.
+
+This fix is not invoked during :doc:`energy minimization <minimize>`.

 Restrictions
 """"""""""""

 This fix is part of the MISC package.  It is only enabled if LAMMPS
-was built with that package.  See the :doc:`Build package <Build_package>` doc page for more info.
+was built with that package.  See the :doc:`Build package
+<Build_package>` doc page for more info.

 This fix should only be used with fcc or bcc lattices.

--- a/doc/src/fix_orient_eco.rst
+++ b/doc/src/fix_orient_eco.rst
@ -26,22 +26,24 @@ Examples
 Description
 """""""""""

-The fix applies a synthetic driving force to a grain boundary which can
-be used for the investigation of grain boundary motion. The affiliation
-of atoms to either of the two grains forming the grain boundary is
-determined from an orientation-dependent order parameter as described
-in :ref:`(Ulomek) <Ulomek>`. The potential energy of atoms is either increased by an amount
-of 0.5*\ *u0* or -0.5*\ *u0* according to the orientation of the surrounding
-crystal. This creates a potential energy gradient which pushes atoms near
-the grain boundary to orient according to the energetically favorable
-grain orientation. This fix is designed for applications in bicrystal system
-with one grain boundary and open ends, or two opposite grain boundaries in
-a periodic system. In either case, the entire system can experience a
-displacement during the simulation which needs to be accounted for in the
-evaluation of the grain boundary velocity. While the basic method is
-described in :ref:`(Ulomek) <Ulomek>`, the implementation follows the efficient
-implementation from :ref:`(Schratt & Mohles) <Schratt>`. The synthetic potential energy added to an
-atom j is given by the following formulas
+The fix applies a synthetic driving force to a grain boundary which
+can be used for the investigation of grain boundary motion. The
+affiliation of atoms to either of the two grains forming the grain
+boundary is determined from an orientation-dependent order parameter
+as described in :ref:`(Ulomek) <Ulomek>`. The potential energy of
+atoms is either increased by an amount of 0.5*\ *u0* or -0.5*\ *u0*
+according to the orientation of the surrounding crystal. This creates
+a potential energy gradient which pushes atoms near the grain boundary
+to orient according to the energetically favorable grain
+orientation. This fix is designed for applications in bicrystal system
+with one grain boundary and open ends, or two opposite grain
+boundaries in a periodic system. In either case, the entire system can
+experience a displacement during the simulation which needs to be
+accounted for in the evaluation of the grain boundary velocity. While
+the basic method is described in :ref:`(Ulomek) <Ulomek>`, the
+implementation follows the efficient implementation from
+:ref:`(Schratt & Mohles) <Schratt>`. The synthetic potential energy
+added to an atom j is given by the following formulas

 .. math::

@ -60,60 +62,69 @@ atom j is given by the following formulas
 which are fully explained in :ref:`(Ulomek) <Ulomek>`
 and :ref:`(Schratt & Mohles) <Schratt>`.

-The force on each atom is the negative gradient of the synthetic potential energy. It
-depends on the surrounding of this atom. An atom far from the grain boundary does not
-experience a synthetic force as its surrounding is that of an oriented single crystal
-and thermal fluctuations are masked by the parameter *eta*\ . Near the grain boundary
-however, the gradient is nonzero and synthetic force terms are computed.
-The orientationsFile specifies the perfect oriented crystal basis vectors for the
-two adjoining crystals. The first three lines (line=row vector) for the energetically penalized and the
-last three lines for the energetically favored grain assuming *u0* is positive. For
-negative *u0*, this is reversed. With the *cutoff* parameter, the size of the region around
-each atom which is used in the order parameter computation is defined. The cutoff must be
-smaller than the interaction range of the MD potential. It should at
-least include the nearest neighbor shell. For high temperatures or low angle
-grain boundaries, it might be beneficial to increase the cutoff in order to get a more
-precise identification of the atoms surrounding. However, computation time will
-increase as more atoms are considered in the order parameter and force computation.
-It is also worth noting that the cutoff radius must not exceed the communication
-distance for ghost atoms in LAMMPS. With orientationsFile, the
-6 oriented crystal basis vectors is specified. Each line of the input file
-contains the three components of a primitive lattice vector oriented according to
-the grain orientation in the simulation box. The first (last) three lines correspond
-to the primitive lattice vectors of the first (second) grain. An example for
-a :math:`\Sigma\langle001\rangle` mis-orientation is given at the end.
-
-If no synthetic energy difference between the grains is created, :math:`u0=0`, the
-force computation is omitted. In this case, still, the order parameter of the
-driving force is computed and can be used to track the grain boundary motion throughout the
-simulation.
-
+The force on each atom is the negative gradient of the synthetic
+potential energy. It depends on the surrounding of this atom. An atom
+far from the grain boundary does not experience a synthetic force as
+its surrounding is that of an oriented single crystal and thermal
+fluctuations are masked by the parameter *eta*\ . Near the grain
+boundary however, the gradient is nonzero and synthetic force terms
+are computed.  The orientationsFile specifies the perfect oriented
+crystal basis vectors for the two adjoining crystals. The first three
+lines (line=row vector) for the energetically penalized and the last
+three lines for the energetically favored grain assuming *u0* is
+positive. For negative *u0*, this is reversed. With the *cutoff*
+parameter, the size of the region around each atom which is used in
+the order parameter computation is defined. The cutoff must be smaller
+than the interaction range of the MD potential. It should at least
+include the nearest neighbor shell. For high temperatures or low angle
+grain boundaries, it might be beneficial to increase the cutoff in
+order to get a more precise identification of the atoms
+surrounding. However, computation time will increase as more atoms are
+considered in the order parameter and force computation.  It is also
+worth noting that the cutoff radius must not exceed the communication
+distance for ghost atoms in LAMMPS. With orientationsFile, the 6
+oriented crystal basis vectors is specified. Each line of the input
+file contains the three components of a primitive lattice vector
+oriented according to the grain orientation in the simulation box. The
+first (last) three lines correspond to the primitive lattice vectors
+of the first (second) grain. An example for a
+:math:`\Sigma\langle001\rangle` mis-orientation is given at the end.

+If no synthetic energy difference between the grains is created,
+:math:`u0=0`, the force computation is omitted. In this case, still,
+the order parameter of the driving force is computed and can be used
+to track the grain boundary motion throughout the simulation.

 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-No information about this fix is written to :doc: `binary restart files <restart>`.
+No information about this fix is written to :doc: `binary restart
+files <restart>`.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this fix to
-add the potential energy of atom interactions with the grain boundary
-driving force to the system's potential energy as part of thermodynamic output.
-The total sum of added synthetic potential energy is computed and can be accessed
-by various output options. The order parameter as well as the thermally masked
-output parameter are stored in per-atom arrays and can also be accessed by various
-:doc:`output commands <Howto_output>`.
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the potential energy of atom interactions with the
+grain boundary driving force to the global potential energy of the
+system as part of :doc:`thermodynamic output <thermo_style>`.  The
+default setting for this fix is :doc:`fix_modify energy no
+<fix_modify>`.

-No parameter of this fix can be used with the start/stop keywords of the run command. This fix is
-not invoked during energy minimization.
+This fix calculates a per-atom array with 2 columns, which can be
+accessed by indices 1-1 by any command that uses per-atom values from
+a fix as input.  See the :doc:`Howto output <Howto_output>` doc page
+for an overview of LAMMPS output options.

+The first column is the order parameter for each atom; the second is
+the thermal masking value for each atom.  Both are described above.

+No parameter of this fix can be used with the start/stop keywords of
+the run command. This fix is not invoked during energy minimization.

 Restrictions
 """"""""""""

-This fix is part of the USER-MISC package. It is only enabled if LAMMPS was
-built with that package. See the :doc:`Build package <Build_package>` doc page for more info.
-
+This fix is part of the USER-MISC package. It is only enabled if
+LAMMPS was built with that package. See the :doc:`Build package
+<Build_package>` doc page for more info.


 Related commands
--- a/doc/src/fix_plumed.rst
+++ b/doc/src/fix_plumed.rst
@ -66,7 +66,8 @@ plumed fix in the LAMMPS input.

 The *plumedfile* keyword allows the user to specify the name of the
 PLUMED input file.  Instructions as to what should be included in a
-plumed input file can be found in the `documentation for PLUMED <plumeddocs_>`_
+plumed input file can be found in the `documentation for PLUMED
+<plumeddocs_>`_

 The *outfile* keyword allows the user to specify the name of a file in
 which to output the PLUMED log.  This log file normally just repeats the
@ -78,31 +79,45 @@ be specified by the user in the PLUMED input file.
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-When performing a restart of a calculation that involves PLUMED you must
-include a RESTART command in the PLUMED input file as detailed in the
-`PLUMED documentation <plumeddocs_>`_.  When the restart command is found in
-the PLUMED input PLUMED will append to the files that were generated in
-the run that was performed previously.  No part of the PLUMED restart
-data is included in the LAMMPS restart files.  Furthermore, any history
-dependent bias potentials that were accumulated in previous calculations
-will be read in when the RESTART command is included in the PLUMED
-input.
+When performing a restart of a calculation that involves PLUMED you
+must include a RESTART command in the PLUMED input file as detailed in
+the `PLUMED documentation <plumeddocs_>`_.  When the restart command
+is found in the PLUMED input PLUMED will append to the files that were
+generated in the run that was performed previously.  No part of the
+PLUMED restart data is included in the LAMMPS restart files.
+Furthermore, any history dependent bias potentials that were
+accumulated in previous calculations will be read in when the RESTART
+command is included in the PLUMED input.

-The :doc:`fix_modify <fix_modify>` *energy* option is not supported by
-this fix.
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the energy change from the biasing force added by
+PLUMED to the global potential energy of the system as part of
+:doc:`thermodynamic output <thermo_style>`.  The default setting for
+this fix is :doc:`fix_modify energy yes <fix_modify>`.

-Nothing is computed by this fix that can be accessed by any of the
-:doc:`output commands <Howto_output>` within LAMMPS.  All the quantities
-of interest can be output by commands that are native to PLUMED,
-however.
+The :doc:`fix_modify <fix_modify>` *virial* option is supported by
+this fix to add the contribution from the biasing force to the global
+pressure of the system via the :doc:`compute pressure
+<compute_pressure>` command.  This can be accessed by
+:doc:`thermodynamic output <thermo_style>`.  The default setting for
+this fix is :doc:`fix_modify virial yes <fix_modify>`.
+
+This fix computes a global scalar which can be accessed by various
+:doc:`output commands <Howto_output>`.  The scalar is the PLUMED
+energy mentioned above.  The scalar value calculated by this fix is
+"extensive".
+
+Note that other quantities of interest can be output by commands that
+are native to PLUMED.

 Restrictions
 """"""""""""

 This fix is part of the USER-PLUMED package.  It is only enabled if
-LAMMPS was built with that package.  See the :doc:`Build package <Build_package>` doc page for more info.
+LAMMPS was built with that package.  See the :doc:`Build package
+<Build_package>` doc page for more info.

-There can only be one plumed fix active at a time.
+There can only be one fix plumed command active at a time.

 Related commands
 """"""""""""""""
--- a/doc/src/fix_poems.rst
+++ b/doc/src/fix_poems.rst
@ -39,14 +39,15 @@ useful for treating a large biomolecule as a collection of connected,
 coarse-grained particles.

 The coupling, associated motion constraints, and time integration is
-performed by the software package `Parallelizable Open source Efficient Multibody Software (POEMS)` which computes the
-constrained rigid-body motion of articulated (jointed) multibody
-systems :ref:`(Anderson) <Anderson>`.  POEMS was written and is distributed
-by Prof Kurt Anderson, his graduate student Rudranarayan Mukherjee,
-and other members of his group at Rensselaer Polytechnic Institute
-(RPI).  Rudranarayan developed the LAMMPS/POEMS interface.  For
-copyright information on POEMS and other details, please refer to the
-documents in the poems directory distributed with LAMMPS.
+performed by the software package `Parallelizable Open source
+Efficient Multibody Software (POEMS)` which computes the constrained
+rigid-body motion of articulated (jointed) multibody systems
+:ref:`(Anderson) <Anderson>`.  POEMS was written and is distributed by
+Prof Kurt Anderson, his graduate student Rudranarayan Mukherjee, and
+other members of his group at Rensselaer Polytechnic Institute (RPI).
+Rudranarayan developed the LAMMPS/POEMS interface.  For copyright
+information on POEMS and other details, please refer to the documents
+in the poems directory distributed with LAMMPS.

 This fix updates the positions and velocities of the rigid atoms with
 a constant-energy time integration, so you should not update the same
@ -107,7 +108,16 @@ off, and there is only a single fix poems defined.
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-No information about this fix is written to :doc:`binary restart files <restart>`.
+No information about this fix is written to :doc:`binary restart files
+<restart>`.
+
+The :doc:`fix_modify <fix_modify>` *virial* option is supported by
+this fix to add the contribution due to the added forces and torques
+on atoms to both the global pressure and per-atom stress of the system
+via the :doc:`compute pressure <compute_pressure>` and :doc:`compute
+stress/atom <compute_stress_atom>` commands.  The former can be
+accessed by :doc:`thermodynamic output <thermo_style>`.  The default
+setting for this fix is :doc:`fix_modify virial yes <fix_modify>`.

 The :doc:`fix_modify <fix_modify>` *bodyforces* option is supported by
 this fix style to set whether per-body forces and torques are computed
@ -115,16 +125,18 @@ early or late in a timestep, i.e. at the post-force stage or at the
 final-integrate stage, respectively.

 No global or per-atom quantities are stored by this fix for access by
-various :doc:`output commands <Howto_output>`.  No parameter of this fix
-can be used with the *start/stop* keywords of the :doc:`run <run>`
-command.  This fix is not invoked during :doc:`energy minimization <minimize>`.
+various :doc:`output commands <Howto_output>`.  No parameter of this
+fix can be used with the *start/stop* keywords of the :doc:`run <run>`
+command.  This fix is not invoked during :doc:`energy minimization
+<minimize>`.

 Restrictions
 """"""""""""

-This fix is part of the :ref:`POEMS <PKG-POEMS>` package.  It is only enabled if LAMMPS
-was built with that package, which also requires the POEMS library be
-built and linked with LAMMPS.  See the :doc:`Build package <Build_package>` doc page for more info.
+This fix is part of the :ref:`POEMS <PKG-POEMS>` package.  It is only
+enabled if LAMMPS was built with that package, which also requires the
+POEMS library be built and linked with LAMMPS.  See the :doc:`Build
+package <Build_package>` doc page for more info.

 Related commands
 """"""""""""""""
--- a/doc/src/fix_precession_spin.rst
+++ b/doc/src/fix_precession_spin.rst
@ -42,11 +42,12 @@ Examples
 Description
 """""""""""

-This fix applies a precession torque to each magnetic spin in the group.
+This fix applies a precession torque to each magnetic spin in the
+group.

-Style *zeeman* is used for the simulation of the interaction
-between the magnetic spins in the defined group and an external
-magnetic field:
+Style *zeeman* is used for the simulation of the interaction between
+the magnetic spins in the defined group and an external magnetic
+field:

 .. math::

@ -62,17 +63,17 @@ with:

 The field value in Tesla is multiplied by the gyromagnetic
 ratio, :math:`g \cdot \mu_B/\hbar`, converting it into a precession frequency in
-rad.THz (in metal units and with :math:`\mu_B = 5.788\cdot 10^{-5}` eV/T).
+rad.THz (in metal units and with :math:`\mu_B = 5.788\cdot 10^{-5}`
+eV/T).

-As a comparison, the figure below displays the simulation of a
-single spin (of norm :math:`\mu_i = 1.0`) submitted to an external
-magnetic field of :math:`\vert B_{ext}\vert = 10.0\; \mathrm{Tesla}` (and oriented along the z
-axis).
-The upper plot shows the average magnetization along the
-external magnetic field axis and the lower plot the Zeeman
-energy, both as a function of temperature.
-The reference result is provided by the plot of the Langevin
-function for the same parameters.
+As a comparison, the figure below displays the simulation of a single
+spin (of norm :math:`\mu_i = 1.0`) submitted to an external magnetic
+field of :math:`\vert B_{ext}\vert = 10.0\; \mathrm{Tesla}` (and
+oriented along the z axis).  The upper plot shows the average
+magnetization along the external magnetic field axis and the lower
+plot the Zeeman energy, both as a function of temperature.  The
+reference result is provided by the plot of the Langevin function for
+the same parameters.

 .. image:: JPG/zeeman_langevin.jpg
   :align: center
@ -88,10 +89,12 @@ for the magnetic spins in the defined group:

 .. math::

-   H_{aniso}  = -\sum_{{ i}=1}^{N} K_{an}(\mathbf{r}_{i})\, \left( \vec{s}_{i} \cdot \vec{n}_{i} \right)^2
+   H_{aniso} = -\sum_{{ i}=1}^{N} K_{an}(\mathbf{r}_{i})\, \left(
+   \vec{s}_{i} \cdot \vec{n}_{i} \right)^2

-with :math:`n` defining the direction of the anisotropy, and :math:`K` (in eV) its intensity.
-If :math:`K > 0`, an easy axis is defined, and if :math:`K < 0`, an easy plane is defined.
+with :math:`n` defining the direction of the anisotropy, and :math:`K`
+(in eV) its intensity.  If :math:`K > 0`, an easy axis is defined, and
+if :math:`K < 0`, an easy plane is defined.

 Style *cubic* is used to simulate a cubic anisotropy, with three
 possible easy axis for the magnetic spins in the defined group:
@ -110,17 +113,17 @@ possible easy axis for the magnetic spins in the defined group:
   \left(\vec{s}_{i} \cdot \vec{n_2} \right)^2
   \left(\vec{s}_{i} \cdot \vec{n_3} \right)^2

-with :math:`K_1` and :math:`K_{2c}` (in eV) the intensity coefficients and
-:math:`\vec{n}_1`, :math:`\vec{n}_2` and :math:`\vec{n}_3` defining the three anisotropic directions
-defined by the command (from *n1x* to *n3z*).
-For :math:`\vec{n}_1 = (1 0 0)`, :math:`\vec{n}_2 = (0 1 0)`, and :math:`\vec{n}_3 = (0 0 1)`, :math:`K_1 < 0` defines an
+with :math:`K_1` and :math:`K_{2c}` (in eV) the intensity coefficients
+and :math:`\vec{n}_1`, :math:`\vec{n}_2` and :math:`\vec{n}_3`
+defining the three anisotropic directions defined by the command (from
+*n1x* to *n3z*).  For :math:`\vec{n}_1 = (1 0 0)`, :math:`\vec{n}_2 =
+(0 1 0)`, and :math:`\vec{n}_3 = (0 0 1)`, :math:`K_1 < 0` defines an
 iron type anisotropy (easy axis along the :math:`(0 0 1)`-type cube
-edges), and :math:`K_1 > 0` defines a nickel type anisotropy (easy axis
-along the :math:`(1 1 1)`-type cube diagonals).
-:math:`K_2^c > 0` also defines easy axis along the :math:`(1 1 1)`-type cube
-diagonals.
-See chapter 2 of :ref:`(Skomski) <Skomski1>` for more details on cubic
-anisotropies.
+edges), and :math:`K_1 > 0` defines a nickel type anisotropy (easy
+axis along the :math:`(1 1 1)`-type cube diagonals).  :math:`K_2^c >
+0` also defines easy axis along the :math:`(1 1 1)`-type cube
+diagonals.  See chapter 2 of :ref:`(Skomski) <Skomski1>` for more
+details on cubic anisotropies.

 In all cases, the choice of :math:`(x y z)` only imposes the vector
 directions for the forces. Only the direction of the vector is
@ -134,32 +137,35 @@ Those styles can be combined within one single command line.
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-By default, the energy associated to this fix is not added to the potential
-energy of the system.
-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this fix
-to add this magnetic potential energy to the potential energy of the system,
+No information about this fix is written to :doc:`binary restart files
+<restart>`.

-.. code-block:: LAMMPS
-
-   fix             1 all precession/spin zeeman 1.0 0.0 0.0 1.0
-   fix_modify      1 energy yes
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the energy associated with the spin precession
+torque to the global potential energy of the system as part of
+:doc:`thermodynamic output <thermo_style>`.  The default setting for
+this fix is :doc:`fix_modify energy no <fix_modify>`.

 This fix computes a global scalar which can be accessed by various
-:doc:`output commands <Howto_output>`.
+:doc:`output commands <Howto_output>`.  The scalar is the potential
+energy (in energy units) discussed in the previous paragraph.  The
+scalar value is an "extensive" quantity.

-No information about this fix is written to :doc:`binary restart files <restart>`.
+No information about this fix is written to :doc:`binary restart files
+<restart>`.

 Restrictions
 """"""""""""

 The *precession/spin* style is part of the SPIN package.  This style
 is only enabled if LAMMPS was built with this package, and if the
-atom_style "spin" was declared.  See the :doc:`Build package <Build_package>` doc page for more info.
+atom_style "spin" was declared.  See the :doc:`Build package
+<Build_package>` doc page for more info.

-The *precession/spin* style can only be declared once. If more
-than one precession type (for example combining an anisotropy and a Zeeman interactions)
-has to be declared, they have to be chained in the same command
-line (as shown in the examples above).
+The *precession/spin* style can only be declared once. If more than
+one precession type (for example combining an anisotropy and a Zeeman
+interactions) has to be declared, they have to be chained in the same
+command line (as shown in the examples above).

 Related commands
 """"""""""""""""
--- a/doc/src/fix_qbmsst.rst
+++ b/doc/src/fix_qbmsst.rst
@ -152,13 +152,30 @@ Because the state of the random number generator is not written to
 "exactly" in an uninterrupted fashion. However, in a statistical
 sense, a restarted simulation should produce similar behaviors of the
 system as if it is not interrupted.  To achieve such a restart, one
-should write explicitly the same value for *q*\ , *mu*\ , *damp*\ , *f_max*,
-*N_f*, *eta*\ , and *beta* and set *tscale* = 0 if the system is
-compressed during the first run.
+should write explicitly the same value for *q*\ , *mu*\ , *damp*\ ,
+*f_max*, *N_f*, *eta*\ , and *beta* and set *tscale* = 0 if the system
+is compressed during the first run.
+
+The cumulative energy change in the system imposed by this fix is
+included in the :doc:`thermodynamic output <thermo_style>` keywords
+*ecouple* and *econserve*.  See the :doc:`thermo_style <thermo_style>`
+doc page for details.
+
+This fix computes a global scalar which can be accessed by various
+:doc:`output commands <Howto_output>`.  The scalar is the same
+cumulative energy change due to this fix described in the previous
+paragraph.  The scalar value calculated by this fix is "extensive".

 The progress of the QBMSST can be monitored by printing the global
-scalar and global vector quantities computed by the fix.  The global
-vector contains five values in this order:
+scalar and global vector quantities computed by the fix.
+
+As mentioned above, the scalar is the cumulative energy change due to
+the fix.  By monitoring the thermodynamic *econserve* output, this can
+be used to test if the MD timestep is sufficiently small for accurate
+integration of the dynamic equations.
+
+The global vector contains five values in the following order.  The
+vector values output by this fix are "intensive".

 [\ *dhugoniot*\ , *drayleigh*\ , *lagrangian_speed*, *lagrangian_position*,
 *quantum_temperature*]
@ -170,29 +187,21 @@ vector contains five values in this order:
 5. *quantum_temperature* is the temperature of the quantum thermal bath :math:`T^{qm}`.

 To print these quantities to the log file with descriptive column
-headers, the following LAMMPS commands are suggested. Here the
-:doc:`fix_modify <fix_modify>` energy command is also enabled to allow
-the thermo keyword *etotal* to print the quantity :math:`E^{tot}`.  See
-also the :doc:`thermo_style <thermo_style>` command.
+headers, the following LAMMPS commands are suggested.

 .. parsed-literal::

   fix             fix_id all msst z
-   fix_modify      fix_id energy yes
   variable        dhug    equal f_fix_id[1]
   variable        dray    equal f_fix_id[2]
   variable        lgr_vel equal f_fix_id[3]
   variable        lgr_pos equal f_fix_id[4]
   variable        T_qm    equal f_fix_id[5]
-   thermo_style    custom  step temp ke pe lz pzz etotal v_dhug v_dray v_lgr_vel v_lgr_pos v_T_qm f_fix_id
+   thermo_style    custom  step temp ke pe lz pzz econserve v_dhug v_dray v_lgr_vel v_lgr_pos v_T_qm f_fix_id

-The global scalar under the entry f_fix_id is the quantity of thermo
-energy as an extra part of :math:`E^{tot}`. This global scalar and the
-vector of 5 quantities can be accessed by various :doc:`output commands <Howto_output>`.
-It is worth noting that the temp keyword
-under the :doc:`thermo_style <thermo_style>` command print the
-instantaneous classical temperature :math:`T^{cl}` as described
-in the command :doc:`fix qtb <fix_qtb>`.
+It is worth noting that the temp keyword for the :doc:`thermo_style
+<thermo_style>` command prints the instantaneous classical temperature
+:math:`T^{cl}` as described by the :doc:`fix qtb <fix_qtb>` command.

 ----------

@ -200,7 +209,8 @@ Restrictions
 """"""""""""

 This fix style is part of the USER-QTB package.  It is only enabled if
-LAMMPS was built with that package. See the :doc:`Build package <Build_package>` doc page for more info.
+LAMMPS was built with that package. See the :doc:`Build package
+<Build_package>` doc page for more info.

 All cell dimensions must be periodic. This fix can not be used with a
 triclinic cell.  The QBMSST fix has been tested only for the group-ID
--- a/doc/src/fix_restrain.rst
+++ b/doc/src/fix_restrain.rst
@ -219,15 +219,19 @@ current dihedral angle :math:`\phi` is equal to :math:`\phi_0`.
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

-No information about this fix is written to :doc:`binary restart files <restart>`.
+No information about this fix is written to :doc:`binary restart files
+<restart>`.

-The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
-fix to add the potential energy associated with this fix to the
-system's potential energy as part of :doc:`thermodynamic output <thermo_style>`.
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the potential energy associated with this fix to the
+global potential energy of the system as part of :doc:`thermodynamic
+output <thermo_style>` The default setting for this fix is
+:doc:`fix_modify energy no <fix_modify>`.

 The :doc:`fix_modify <fix_modify>` *respa* option is supported by this
-fix. This allows to set at which level of the :doc:`r-RESPA <run_style>`
-integrator the fix is adding its forces. Default is the outermost level.
+fix. This allows to set at which level of the :doc:`r-RESPA
+<run_style>` integrator the fix is adding its forces. Default is the
+outermost level.

 .. note::

--- a/doc/src/fix_rhok.rst
+++ b/doc/src/fix_rhok.rst
@ -45,11 +45,34 @@ temperatures :ref:`(Pedersen) <Pedersen>`.
 An example of using the interface pinning method is located in the
 *examples/USER/misc/rhok* directory.

+Restart, fix_modify, output, run start/stop, minimize info
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+No information about this fix is written to :doc:`binary restart files
+<restart>`.
+
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by
+this fix to add the potential energy calculated by the fix to the
+global potential energy of the system as part of :doc:`thermodynamic
+output <thermo_style>`.  The default setting for this fix is
+:doc:`fix_modify energy no <fix_modify>`.
+
+This fix computes a global scalar which can be accessed by various
+:doc:`output commands <Howto_output>`.  The scalar is the potential
+energy discussed in the preceding paragraph.  The scalar stored by
+this fix is "extensive".
+
+No parameter of this fix can be used with the *start/stop* keywords of
+the :doc:`run <run>` command.
+
+This fix is not invoked during :doc:`energy minimization <minimize>`.
+
 Restrictions
 """"""""""""

 This fix is part of the USER-MISC package.  It is only enabled if
-LAMMPS was built with that package.  See the :doc:`Build package <Build_package>` doc page for more info.
+LAMMPS was built with that package.  See the :doc:`Build package
+<Build_package>` doc page for more info.

 Related commands
 """"""""""""""""
--- a/Show More
+++ b/Show More