Merge remote-tracking branch 'lammps/develop' into electrode

.github/CONTRIBUTING.md

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

cmake/CMakeLists.txt

@@ -105,8 +105,28 @@ if(CMAKE_CXX_COMPILER_ID STREQUAL "Intel")
   endif()
 endif()
 
+# silence excessive warnings for new Intel Compilers
+if(CMAKE_CXX_COMPILER_ID STREQUAL "IntelLLVM")
+  set(CMAKE_TUNE_DEFAULT "-Wno-tautological-constant-compare")
+endif()
+
+# silence excessive warnings for PGI/NVHPC compilers
+if((CMAKE_CXX_COMPILER_ID STREQUAL "NVHPC") OR (CMAKE_CXX_COMPILER_ID STREQUAL "PGI"))
+  set(CMAKE_TUNE_DEFAULT "-Minform=severe")
+endif()
+
+# silence nvcc warnings
+if((PKG_KOKKOS) AND (Kokkos_ENABLE_CUDA))
+  set(CMAKE_TUNE_DEFAULT "${CMAKE_TUNE_DEFAULT} -Xcudafe --diag_suppress=unrecognized_pragma")
+endif()
+
 # we require C++11 without extensions. Kokkos requires at least C++14 (currently)
-set(CMAKE_CXX_STANDARD 11)
+if(NOT CMAKE_CXX_STANDARD)
+  set(CMAKE_CXX_STANDARD 11)
+endif()
+if(CMAKE_CXX_STANDARD LESS 11)
+  message(FATAL_ERROR "C++ standard must be set to at least 11")
+endif()
 if(PKG_KOKKOS AND (CMAKE_CXX_STANDARD LESS 14))
   set(CMAKE_CXX_STANDARD 14)
 endif()
@@ -172,6 +192,7 @@ set(STANDARD_PACKAGES
   AWPMD
   BOCS
   BODY
+  BPM
   BROWNIAN
   CG-DNA
   CG-SDK
@@ -469,6 +490,7 @@ set(CMAKE_TUNE_FLAGS "${CMAKE_TUNE_DEFAULT}" CACHE STRING "Compiler and machine
 separate_arguments(CMAKE_TUNE_FLAGS)
 foreach(_FLAG ${CMAKE_TUNE_FLAGS})
   target_compile_options(lammps PRIVATE ${_FLAG})
+  target_compile_options(lmp PRIVATE ${_FLAG})
 endforeach()
 ########################################################################
 # Basic system tests (standard libraries, headers, functions, types)   #

cmake/Modules/Documentation.cmake

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

cmake/Modules/Packages/COLVARS.cmake

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

cmake/Modules/Packages/GPU.cmake

@@ -80,14 +80,14 @@ if(GPU_API STREQUAL "CUDA")
   # --arch translates directly instead of JIT, so this should be for the preferred or most common architecture
   set(GPU_CUDA_GENCODE "-arch=${GPU_ARCH}")
 
-  # apply the following to build "fat" CUDA binaries only for known CUDA toolkits
+  # apply the following to build "fat" CUDA binaries only for known CUDA toolkits since version 8.0
-  if(CUDA_VERSION VERSION_GREATER_EQUAL "12.0")
+  # only the Kepler achitecture and beyond is supported
-    message(WARNING "Untested CUDA Toolkit version. Use at your own risk")
+  if(CUDA_VERSION VERSION_LESS 8.0)
+    message(FATAL_ERROR "CUDA Toolkit version 8.0 or later is required")
+  elseif(CUDA_VERSION VERSION_GREATER_EQUAL "12.0")
+    message(WARNING "Untested CUDA Toolkit version ${CUDA_VERSION}. Use at your own risk")
+    set(GPU_CUDA_GENCODE "-arch=all")
   else()
-    # Fermi (GPU Arch 2.x) is supported by CUDA 3.2 to CUDA 8.0
-    if((CUDA_VERSION VERSION_GREATER_EQUAL "3.2") AND (CUDA_VERSION VERSION_LESS "9.0"))
-      string(APPEND GPU_CUDA_GENCODE " -gencode arch=compute_20,code=[sm_20,compute_20] ")
-    endif()
     # Kepler (GPU Arch 3.0) is supported by CUDA 5 to CUDA 10.2
     if((CUDA_VERSION VERSION_GREATER_EQUAL "5.0") AND (CUDA_VERSION VERSION_LESS "11.0"))
       string(APPEND GPU_CUDA_GENCODE " -gencode arch=compute_30,code=[sm_30,compute_30] ")
@@ -120,6 +120,14 @@ if(GPU_API STREQUAL "CUDA")
     if(CUDA_VERSION VERSION_GREATER_EQUAL "11.1")
       string(APPEND GPU_CUDA_GENCODE " -gencode arch=compute_86,code=[sm_86,compute_86]")
     endif()
+    # Hopper (GPU Arch 9.0) is supported by CUDA 12.0? and later
+    if(CUDA_VERSION VERSION_GREATER_EQUAL "12.0")
+      string(APPEND GPU_CUDA_GENCODE " -gencode arch=compute_90,code=[sm_90,compute_90]")
+    endif()
+    #    # Lovelace (GPU Arch 9.x) is supported by CUDA 12.0? and later
+    #if(CUDA_VERSION VERSION_GREATER_EQUAL "12.0")
+    #  string(APPEND GPU_CUDA_GENCODE " -gencode arch=compute_9x,code=[sm_9x,compute_9x]")
+    #endif()
   endif()
 
   cuda_compile_fatbin(GPU_GEN_OBJS ${GPU_LIB_CU} OPTIONS ${CUDA_REQUEST_PIC}
@@ -257,43 +265,47 @@ elseif(GPU_API STREQUAL "HIP")
     find_package(CUDA REQUIRED)
     set(HIP_ARCH "sm_50" CACHE STRING "HIP primary CUDA architecture (e.g. sm_60)")
 
-    # build arch/gencode commands for nvcc based on CUDA toolkit version and use choice
+    if(CUDA_VERSION VERSION_LESS 8.0)
-    # --arch translates directly instead of JIT, so this should be for the preferred or most common architecture
+      message(FATAL_ERROR "CUDA Toolkit version 8.0 or later is required")
-    set(HIP_CUDA_GENCODE "-arch=${HIP_ARCH}")
+    elseif(CUDA_VERSION VERSION_GREATER_EQUAL "12.0")
-    # Fermi (GPU Arch 2.x) is supported by CUDA 3.2 to CUDA 8.0
+      message(WARNING "Untested CUDA Toolkit version ${CUDA_VERSION}. Use at your own risk")
-    if((CUDA_VERSION VERSION_GREATER_EQUAL "3.2") AND (CUDA_VERSION VERSION_LESS "9.0"))
+      set(HIP_CUDA_GENCODE "-arch=all")
-      string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_20,code=[sm_20,compute_20]")
+    else()
-    endif()
+      # build arch/gencode commands for nvcc based on CUDA toolkit version and use choice
-    # Kepler (GPU Arch 3.0) is supported by CUDA 5 to CUDA 10.2
+      # --arch translates directly instead of JIT, so this should be for the preferred or most common architecture
-    if((CUDA_VERSION VERSION_GREATER_EQUAL "5.0") AND (CUDA_VERSION VERSION_LESS "11.0"))
+      set(HIP_CUDA_GENCODE "-arch=${HIP_ARCH}")
-      string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_30,code=[sm_30,compute_30]")
+      # Kepler (GPU Arch 3.0) is supported by CUDA 5 to CUDA 10.2
-    endif()
+      if((CUDA_VERSION VERSION_GREATER_EQUAL "5.0") AND (CUDA_VERSION VERSION_LESS "11.0"))
-    # Kepler (GPU Arch 3.5) is supported by CUDA 5 to CUDA 11.0
+        string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_30,code=[sm_30,compute_30]")
-    if((CUDA_VERSION VERSION_GREATER_EQUAL "5.0") AND (CUDA_VERSION VERSION_LESS "12.0"))
+      endif()
-      string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_35,code=[sm_35,compute_35]")
+      # Kepler (GPU Arch 3.5) is supported by CUDA 5 to CUDA 11.0
-    endif()
+      if((CUDA_VERSION VERSION_GREATER_EQUAL "5.0") AND (CUDA_VERSION VERSION_LESS "12.0"))
-    # Maxwell (GPU Arch 5.x) is supported by CUDA 6 and later
+        string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_35,code=[sm_35,compute_35]")
-    if(CUDA_VERSION VERSION_GREATER_EQUAL "6.0")
+      endif()
-      string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_50,code=[sm_50,compute_50] -gencode arch=compute_52,code=[sm_52,compute_52]")
+      # Maxwell (GPU Arch 5.x) is supported by CUDA 6 and later
-    endif()
+      if(CUDA_VERSION VERSION_GREATER_EQUAL "6.0")
-    # Pascal (GPU Arch 6.x) is supported by CUDA 8 and later
+        string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_50,code=[sm_50,compute_50] -gencode arch=compute_52,code=[sm_52,compute_52]")
-    if(CUDA_VERSION VERSION_GREATER_EQUAL "8.0")
+      endif()
-      string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_60,code=[sm_60,compute_60] -gencode arch=compute_61,code=[sm_61,compute_61]")
+      # Pascal (GPU Arch 6.x) is supported by CUDA 8 and later
-    endif()
+      if(CUDA_VERSION VERSION_GREATER_EQUAL "8.0")
-    # Volta (GPU Arch 7.0) is supported by CUDA 9 and later
+        string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_60,code=[sm_60,compute_60] -gencode arch=compute_61,code=[sm_61,compute_61]")
-    if(CUDA_VERSION VERSION_GREATER_EQUAL "9.0")
+      endif()
-      string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_70,code=[sm_70,compute_70]")
+      # Volta (GPU Arch 7.0) is supported by CUDA 9 and later
-    endif()
+      if(CUDA_VERSION VERSION_GREATER_EQUAL "9.0")
-    # Turing (GPU Arch 7.5) is supported by CUDA 10 and later
+        string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_70,code=[sm_70,compute_70]")
-    if(CUDA_VERSION VERSION_GREATER_EQUAL "10.0")
+      endif()
-      string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_75,code=[sm_75,compute_75]")
+      # Turing (GPU Arch 7.5) is supported by CUDA 10 and later
-    endif()
+      if(CUDA_VERSION VERSION_GREATER_EQUAL "10.0")
-    # Ampere (GPU Arch 8.0) is supported by CUDA 11 and later
+        string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_75,code=[sm_75,compute_75]")
-    if(CUDA_VERSION VERSION_GREATER_EQUAL "11.0")
+      endif()
-      string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_80,code=[sm_80,compute_80]")
+      # Ampere (GPU Arch 8.0) is supported by CUDA 11 and later
-    endif()
+      if(CUDA_VERSION VERSION_GREATER_EQUAL "11.0")
-    if(CUDA_VERSION VERSION_GREATER_EQUAL "12.0")
+        string(APPEND HIP_CUDA_GENCODE " -gencode arch=compute_80,code=[sm_80,compute_80]")
-      message(WARNING "Unsupported CUDA version. Use at your own risk.")
+      endif()
+      # Hopper (GPU Arch 9.0) is supported by CUDA 12.0? and later
+      if(CUDA_VERSION VERSION_GREATER_EQUAL "12.0")
+        string(APPEND GPU_CUDA_GENCODE " -gencode arch=compute_90,code=[sm_90,compute_90]")
+      endif()
     endif()
   endif()
 

cmake/Modules/Packages/H5MD.cmake

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

cmake/Modules/Packages/MACHDYN.cmake

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

cmake/Modules/Packages/MDI.cmake

@@ -8,8 +8,8 @@ option(DOWNLOAD_MDI "Download and compile the MDI library instead of using an al
 
 if(DOWNLOAD_MDI)
   message(STATUS "MDI download requested - we will build our own")
-  set(MDI_URL "https://github.com/MolSSI-MDI/MDI_Library/archive/v1.3.0.tar.gz" CACHE STRING "URL for MDI tarball")
+  set(MDI_URL "https://github.com/MolSSI-MDI/MDI_Library/archive/v1.3.2.tar.gz" CACHE STRING "URL for MDI tarball")
-  set(MDI_MD5 "8a8da217148bd9b700083b67d795af5e" CACHE STRING "MD5 checksum for MDI tarball")
+  set(MDI_MD5 "836f5da400d8cff0f0e4435640f9454f" CACHE STRING "MD5 checksum for MDI tarball")
   mark_as_advanced(MDI_URL)
   mark_as_advanced(MDI_MD5)
   enable_language(C)

cmake/presets/all_off.cmake

@@ -8,6 +8,7 @@ set(ALL_PACKAGES
   AWPMD
   BOCS
   BODY
+  BPM
   BROWNIAN
   CG-DNA
   CG-SDK

cmake/presets/all_on.cmake

@@ -10,6 +10,7 @@ set(ALL_PACKAGES
   AWPMD
   BOCS
   BODY
+  BPM
   BROWNIAN
   CG-DNA
   CG-SDK

cmake/presets/kokkos-cuda.cmake

@@ -5,7 +5,5 @@ set(PKG_KOKKOS ON CACHE BOOL "" FORCE)
 set(Kokkos_ENABLE_SERIAL ON CACHE BOOL "" FORCE)
 set(Kokkos_ENABLE_OPENMP ON CACHE BOOL "" FORCE)
 set(Kokkos_ENABLE_CUDA   ON CACHE BOOL "" FORCE)
-set(Kokkos_ARCH_MAXWELL50 on CACHE BOOL "" FORCE)
+set(Kokkos_ARCH_PASCAL60 ON CACHE BOOL "" FORCE)
 set(BUILD_OMP ON CACHE BOOL "" FORCE)
-get_filename_component(NVCC_WRAPPER_CMD ${CMAKE_CURRENT_SOURCE_DIR}/../lib/kokkos/bin/nvcc_wrapper ABSOLUTE)
-set(CMAKE_CXX_COMPILER ${NVCC_WRAPPER_CMD} CACHE FILEPATH "" FORCE)

cmake/presets/mingw-cross.cmake

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

cmake/presets/most.cmake

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

cmake/presets/pgi.cmake

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

cmake/presets/windows.cmake

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

doc/.gitignore

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

doc/Makefile

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

doc/msi2lmp.1

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

doc/src/Build_development.rst

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

doc/src/Build_extras.rst

@@ -141,8 +141,6 @@ CMake build
 
 :code:`GPU_ARCH` settings for different GPU hardware is as follows:
 
-* sm_12 or sm_13 for GT200 (supported by CUDA 3.2 until CUDA 6.5)
-* sm_20 or sm_21 for Fermi (supported by CUDA 3.2 until CUDA 7.5)
 * sm_30 for Kepler (supported since CUDA 5 and until CUDA 10.x)
 * sm_35 or sm_37 for Kepler (supported since CUDA 5 and until CUDA 11.x)
 * sm_50 or sm_52 for Maxwell (supported since CUDA 6)
@@ -160,6 +158,12 @@ Thus the GPU_ARCH setting is merely an optimization, to have code for
 the preferred GPU architecture directly included rather than having to wait
 for the JIT compiler of the CUDA driver to translate it.
 
+When compiling for CUDA or HIP with CUDA, version 8.0 or later of the CUDA toolkit
+is required and a GPU architecture of Kepler or later, which must *also* be
+supported by the CUDA toolkit in use **and** the CUDA driver in use.
+When compiling for OpenCL, OpenCL version 1.2 or later is required and the
+GPU must be supported by the GPU driver and OpenCL runtime bundled with the driver.
+
 When building with CMake, you **must NOT** build the GPU library in ``lib/gpu``
 using the traditional build procedure. CMake will detect files generated by that
 process and will terminate with an error and a suggestion for how to remove them.

doc/src/Build_manual.rst

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

doc/src/Commands_bond.rst

@@ -32,6 +32,8 @@ OPT.
    *
    *
    *
+   * :doc:`bpm/rotational <bond_bpm_rotational>`
+   * :doc:`bpm/spring <bond_bpm_spring>`
    * :doc:`class2 (ko) <bond_class2>`
    * :doc:`fene (iko) <bond_fene>`
    * :doc:`fene/expand (o) <bond_fene_expand>`

doc/src/Commands_compute.rst

@@ -33,6 +33,7 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`body/local <compute_body_local>`
    * :doc:`bond <compute_bond>`
    * :doc:`bond/local <compute_bond_local>`
+   * :doc:`born/matrix <compute_born_matrix>`
    * :doc:`centro/atom <compute_centro_atom>`
    * :doc:`centroid/stress/atom <compute_stress_atom>`
    * :doc:`chunk/atom <compute_chunk_atom>`
@@ -91,6 +92,7 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`msd <compute_msd>`
    * :doc:`msd/chunk <compute_msd_chunk>`
    * :doc:`msd/nongauss <compute_msd_nongauss>`
+   * :doc:`nbond/atom <compute_nbond_atom>`
    * :doc:`omega/chunk <compute_omega_chunk>`
    * :doc:`orientorder/atom (k) <compute_orientorder_atom>`
    * :doc:`pair <compute_pair>`

doc/src/Commands_fix.rst

@@ -144,6 +144,7 @@ OPT.
    * :doc:`nve/manifold/rattle <fix_nve_manifold_rattle>`
    * :doc:`nve/noforce <fix_nve_noforce>`
    * :doc:`nve/sphere (ko) <fix_nve_sphere>`
+   * :doc:`nve/bpm/sphere <fix_nve_bpm_sphere>`
    * :doc:`nve/spin <fix_nve_spin>`
    * :doc:`nve/tri <fix_nve_tri>`
    * :doc:`nvk <fix_nvk>`
@@ -161,7 +162,6 @@ OPT.
    * :doc:`orient/fcc <fix_orient>`
    * :doc:`orient/eco <fix_orient_eco>`
    * :doc:`pafi <fix_pafi>`
-   * :doc:`pair/tracker <fix_pair_tracker>`
    * :doc:`phonon <fix_phonon>`
    * :doc:`pimd <fix_pimd>`
    * :doc:`planeforce <fix_planeforce>`

doc/src/Commands_pair.rst

@@ -53,6 +53,7 @@ OPT.
    * :doc:`born/coul/msm (o) <pair_born>`
    * :doc:`born/coul/wolf (go) <pair_born>`
    * :doc:`born/coul/wolf/cs (g) <pair_cs>`
+   * :doc:`bpm/spring <pair_bpm_spring>`
    * :doc:`brownian (o) <pair_brownian>`
    * :doc:`brownian/poly (o) <pair_brownian>`
    * :doc:`buck (giko) <pair_buck>`
@@ -123,7 +124,7 @@ OPT.
    * :doc:`hbond/dreiding/lj (o) <pair_hbond_dreiding>`
    * :doc:`hbond/dreiding/morse (o) <pair_hbond_dreiding>`
    * :doc:`hdnnp <pair_hdnnp>`
-   * :doc:`ilp/graphene/hbn <pair_ilp_graphene_hbn>`
+   * :doc:`ilp/graphene/hbn (t) <pair_ilp_graphene_hbn>`
    * :doc:`ilp/tmd <pair_ilp_tmd>`
    * :doc:`kolmogorov/crespi/full <pair_kolmogorov_crespi_full>`
    * :doc:`kolmogorov/crespi/z <pair_kolmogorov_crespi_z>`
@@ -151,7 +152,7 @@ OPT.
    * :doc:`lj/cut/coul/cut/dielectric (o) <pair_dielectric>`
    * :doc:`lj/cut/coul/cut/soft (o) <pair_fep_soft>`
    * :doc:`lj/cut/coul/debye (gko) <pair_lj_cut_coul>`
-   * :doc:`lj/cut/coul/debye/dielectric <pair_dielectric>`
+   * :doc:`lj/cut/coul/debye/dielectric (o) <pair_dielectric>`
    * :doc:`lj/cut/coul/dsf (gko) <pair_lj_cut_coul>`
    * :doc:`lj/cut/coul/long (gikot) <pair_lj_cut_coul>`
    * :doc:`lj/cut/coul/long/cs <pair_cs>`

doc/src/Developer_utils.rst

@@ -211,6 +211,12 @@ Convenience functions
 .. doxygenfunction:: logmesg(LAMMPS *lmp, const std::string &mesg)
    :project: progguide
 
+.. doxygenfunction:: errorurl
+   :project: progguide
+
+.. doxygenfunction:: missing_cmd_args
+   :project: progguide
+
 .. doxygenfunction:: flush_buffers(LAMMPS *lmp)
    :project: progguide
 

doc/src/Errors.rst

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

doc/src/Errors_details.rst

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

doc/src/Examples.rst

@@ -54,6 +54,8 @@ Lowercase directories
 +-------------+------------------------------------------------------------------+
 | body        | body particles, 2d system                                        |
 +-------------+------------------------------------------------------------------+
+| bpm         | BPM simulations of pouring elastic grains and plate impact       |
++-------------+------------------------------------------------------------------+
 | cmap        | CMAP 5-body contributions to CHARMM force field                  |
 +-------------+------------------------------------------------------------------+
 | colloid     | big colloid particles in a small particle solvent, 2d system     |

doc/src/Howto.rst

@@ -23,6 +23,8 @@ General howto
    Howto_library
    Howto_couple
    Howto_mdi
+   Howto_bpm
+   Howto_broken_bonds
 
 Settings howto
 ==============

doc/src/Howto_bpm.rst

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

doc/src/Howto_broken_bonds.rst

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

doc/src/Howto_elastic.rst

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

doc/src/Modify_atom.rst

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

doc/src/Modify_style.rst

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

doc/src/Packages_details.rst

@@ -32,6 +32,7 @@ page gives those details.
    * :ref:`AWPMD <PKG-AWPMD>`
    * :ref:`BOCS <PKG-BOCS>`
    * :ref:`BODY <PKG-BODY>`
+   * :ref:`BPM <PKG-BPM>`
    * :ref:`BROWNIAN <PKG-BROWNIAN>`
    * :ref:`CG-DNA <PKG-CG-DNA>`
    * :ref:`CG-SDK <PKG-CG-SDK>`
@@ -284,6 +285,33 @@ overview.
 
 ----------
 
+.. _PKG-BPM:
+
+BPM package
+------------
+
+**Contents:**
+
+Pair styles, bond styles, fixes, and computes for bonded particle
+models for mesoscale simulations of solids and fracture.  See the
+:doc:`Howto bpm <Howto_bpm>` page for an overview.
+
+**Authors:** Joel T. Clemmer (Sandia National Labs)
+
+**Supporting info:**
+
+* src/BPM filenames -> commands
+* :doc:`Howto_bpm <Howto_bpm>`
+* :doc:`atom_style bpm/sphere <atom_style>`
+* :doc:`bond_style bpm/rotational <bond_bpm_rotational>`
+* :doc:`bond_style bpm/spring <bond_bpm_spring>`
+* :doc:`compute nbond/atom <compute_nbond_atom>`
+* :doc:`fix nve/bpm/sphere <fix_nve_bpm_sphere>`
+* :doc:`pair_style bpm/spring <pair_bpm_spring>`
+* examples/bpm
+
+----------
+
 .. _PKG-BROWNIAN:
 
 BROWNIAN package
@@ -529,8 +557,20 @@ To use this package, also the :ref:`KSPACE <PKG-KSPACE>` and
 **Supporting info:**
 
 * src/DIELECTRIC: filenames -> commands
+* :doc:`atom_style dielectric <atom_style>`
+* :doc:`pair_style coul/cut/dielectric <pair_dielectric>`
+* :doc:`pair_style coul/long/dielectric <pair_dielectric>`
+* :doc:`pair_style lj/cut/coul/cut/dielectric <pair_dielectric>`
+* :doc:`pair_style lj/cut/coul/debye/dielectric <pair_dielectric>`
+* :doc:`pair_style lj/cut/coul/long/dielectric <pair_dielectric>`
+* :doc:`pair_style lj/cut/coul/msm/dielectric <pair_dielectric>`
+* :doc:`pair_style pppm/dielectric <kspace_style>`
+* :doc:`pair_style pppm/disp/dielectric <kspace_style>`
+* :doc:`pair_style msm/dielectric <kspace_style>`
+* :doc:`fix_style polarize/bem/icc <fix_polarize>`
+* :doc:`fix_style polarize/bem/gmres <fix_polarize>`
+* :doc:`fix_style polarize/functional <fix_polarize>`
 * :doc:`compute efield/atom  <compute_efield_atom>`
-* TODO: add all styles
 * examples/PACKAGES/dielectric
 
 ----------
@@ -1571,7 +1611,6 @@ listing, "ls src/MISC", to see the list of commands.
 * :doc:`pair_style list <pair_list>`
 * :doc:`pair_style srp <pair_srp>`
 * :doc:`pair_style tracker <pair_tracker>`
-* :doc:`fix pair/tracker <fix_pair_tracker>`
 
 ----------
 

doc/src/Packages_list.rst

@@ -58,6 +58,11 @@ whether an extra library is needed to build and use the package:
      - :doc:`Howto body <Howto_body>`
      - body
      - no
+   * - :ref:`BPM <PKG-BPM>`
+     - bonded particle models
+     - :doc:`Howto bpm <Howto_bpm>`
+     - bpm
+     - no
    * - :ref:`BROWNIAN <PKG-BROWNIAN>`
      - Brownian dynamics, self-propelled particles
      - :doc:`fix brownian <fix_brownian>`, :doc:`fix propel/self <fix_propel_self>`

doc/src/Tools.rst

@@ -87,7 +87,7 @@ Miscellaneous tools
 .. table_from_list::
    :columns: 6
 
-   * :ref:`CMake <cmake>`
+   * :ref:`LAMMPS coding standards <coding_standard>`
    * :ref:`emacs <emacs>`
    * :ref:`i-pi <ipi>`
    * :ref:`kate <kate>`
@@ -189,27 +189,32 @@ for the :doc:`chain benchmark <Speed_bench>`.
 
 ----------
 
-.. _cmake:
+.. _coding_standard:
 
-CMake tools
+LAMMPS coding standard
------------
+----------------------
 
-The ``cmbuild`` script is a wrapper around using ``cmake --build <dir>
+The ``coding_standard`` folder contains multiple python scripts to
---target`` and allows compiling LAMMPS in a :ref:`CMake build folder
+check for and apply some LAMMPS coding conventions.  The following
-<cmake_build>` with a make-like syntax regardless of the actual build
+scripts are available:
-tool and the specific name of the program used (e.g. ``ninja-v1.10`` or
-``gmake``) when using ``-D CMAKE_MAKE_PROGRAM=<name>``.
 
 .. parsed-literal::
 
-  Usage: cmbuild [-v] [-h] [-C <dir>] [-j <num>] [<target>]
+   permissions.py   # detects if sources have executable permissions and scripts have not
+   whitespace.py    # detects TAB characters and trailing whitespace
+   homepage.py      # detects outdated LAMMPS homepage URLs (pointing to sandia.gov instead of lammps.org)
+   errordocs.py     # detects deprecated error docs in header files
 
-  Options:
+The tools need to be given the main folder of the LAMMPS distribution
-    -h                print this message
+or individual file names as argument and will by default check them
-    -j <NUM>          allow processing of NUM concurrent tasks
+and report any non-compliance.  With the optional ``-f`` argument the
-    -C DIRECTORY      execute build in folder DIRECTORY
+corresponding script will try to change the non-compliant file(s) to
-    -v                produce verbose output
+match the conventions.
 
+For convenience this scripts can also be invoked by the make file in
+the ``src`` folder with, `make check-whitespace` or `make fix-whitespace`
+to either detect or edit the files.  Correspondingly for the other python
+scripts. `make check` will run all checks.
 
 ----------
 

doc/src/atom_style.rst

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

doc/src/bond_bpm_rotational.rst

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

doc/src/bond_bpm_spring.rst

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

doc/src/bond_quartic.rst

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

doc/src/bond_style.rst

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

doc/src/compute.rst

@@ -179,6 +179,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
 * :doc:`body/local <compute_body_local>` - attributes of body sub-particles
 * :doc:`bond <compute_bond>` - energy of each bond sub-style
 * :doc:`bond/local <compute_bond_local>` - distance and energy of each bond
+* :doc:`born/matrix <compute_born_matrix>` - second derivative or potential with respect to strain
 * :doc:`centro/atom <compute_centro_atom>` - centro-symmetry parameter for each atom
 * :doc:`centroid/stress/atom <compute_stress_atom>` - centroid based stress tensor for each atom
 * :doc:`chunk/atom <compute_chunk_atom>` - assign chunk IDs to each atom
@@ -236,6 +237,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
 * :doc:`msd <compute_msd>` - mean-squared displacement of group of atoms
 * :doc:`msd/chunk <compute_msd_chunk>` - mean-squared displacement for each chunk
 * :doc:`msd/nongauss <compute_msd_nongauss>` - MSD and non-Gaussian parameter of group of atoms
+* :doc:`nbond/atom <compute_nbond_atom>` - calculates number of bonds per atom
 * :doc:`omega/chunk <compute_omega_chunk>` - angular velocity for each chunk
 * :doc:`orientorder/atom <compute_orientorder_atom>` - Steinhardt bond orientational order parameters Ql
 * :doc:`pair <compute_pair>` - values computed by a pair style

doc/src/compute_born_matrix.rst

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

doc/src/compute_nbond_atom.rst

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

doc/src/compute_rigid_local.rst

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

doc/src/compute_temp_profile.rst

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

doc/src/dump_modify.rst

@@ -474,8 +474,9 @@ The *fileper* keyword is documented below with the *nfile* keyword.
 
 The *header* keyword toggles whether the dump file will include a
 header.  Excluding a header will reduce the size of the dump file for
-fixes such as :doc:`fix pair/tracker <fix_pair_tracker>` which do not
+data produced by :doc:`pair tracker <pair_tracker>` or
-require the information typically written to the header.
+:doc:`bpm bond styles <Howto_bpm>` which may not require the
+information typically written to the header.
 
 ----------
 

doc/src/fix.rst

@@ -287,6 +287,7 @@ accelerated styles exist.
 * :doc:`nve/manifold/rattle <fix_nve_manifold_rattle>` -
 * :doc:`nve/noforce <fix_nve_noforce>` - NVE without forces (v only)
 * :doc:`nve/sphere <fix_nve_sphere>` - NVE for spherical particles
+* :doc:`nve/bpm/sphere <fix_nve_bpm_sphere>` - NVE for spherical particles used in the BPM package
 * :doc:`nve/spin <fix_nve_spin>` - NVE for a spin or spin-lattice system
 * :doc:`nve/tri <fix_nve_tri>` - NVE for triangles
 * :doc:`nvk <fix_nvk>` - constant kinetic energy time integration
@@ -304,7 +305,6 @@ accelerated styles exist.
 * :doc:`orient/fcc <fix_orient>` - add grain boundary migration force for FCC
 * :doc:`orient/eco <fix_orient_eco>` - add generalized grain boundary migration force
 * :doc:`pafi <fix_pafi>` - constrained force averages on hyper-planes to compute free energies (PAFI)
-* :doc:`pair/tracker <fix_pair_tracker>` - track properties of pairwise interactions
 * :doc:`phonon <fix_phonon>` - calculate dynamical matrix from MD simulations
 * :doc:`pimd <fix_pimd>` - Feynman path integral molecular dynamics
 * :doc:`planeforce <fix_planeforce>` - constrain atoms to move in a plane

doc/src/fix_adapt.rst

@@ -14,7 +14,7 @@ Syntax
 * adapt = style name of this fix command
 * N = adapt simulation settings every this many timesteps
 * one or more attribute/arg pairs may be appended
-* attribute = *pair* or *bond* or *kspace* or *atom*
+* attribute = *pair* or *bond* or *angle* or *kspace* or *atom*
 
   .. parsed-literal::
 
@@ -28,11 +28,16 @@ Syntax
          bparam = parameter to adapt over time
          I = type bond to set parameter for
          v_name = variable with name that calculates value of bparam
+       *angle* args = astyle aparam I v_name
+         astyle = angle style name, e.g. harmonic
+         aparam = parameter to adapt over time
+         I = type angle to set parameter for
+         v_name = variable with name that calculates value of aparam
        *kspace* arg = v_name
          v_name = variable with name that calculates scale factor on K-space terms
-       *atom* args = aparam v_name
+       *atom* args = atomparam v_name
-         aparam = parameter to adapt over time
+         atomparam = parameter to adapt over time
-         v_name = variable with name that calculates value of aparam
+         v_name = variable with name that calculates value of atomparam
 
 * zero or more keyword/value pairs may be appended
 * keyword = *scale* or *reset* or *mass*
@@ -283,30 +288,62 @@ operates. The only difference is that now a bond coefficient for a
 given bond type is adapted.
 
 A wild-card asterisk can be used in place of or in conjunction with
-the bond type argument to set the coefficients for multiple bond types.
+the bond type argument to set the coefficients for multiple bond
-This takes the form "\*" or "\*n" or "n\*" or "m\*n".  If N = the number of
+types.  This takes the form "\*" or "\*n" or "n\*" or "m\*n".  If N =
-atom types, then an asterisk with no numeric values means all types
+the number of bond types, then an asterisk with no numeric values
-from 1 to N.  A leading asterisk means all types from 1 to n (inclusive).
+means all types from 1 to N.  A leading asterisk means all types from
-A trailing asterisk means all types from n to N (inclusive).  A middle
+1 to n (inclusive).  A trailing asterisk means all types from n to N
-asterisk means all types from m to n (inclusive).
+(inclusive).  A middle asterisk means all types from m to n
+(inclusive).
 
 Currently *bond* does not support bond_style hybrid nor bond_style
-hybrid/overlay as bond styles. The only bonds that currently are
+hybrid/overlay as bond styles. The bond styles that currently work
-working with fix_adapt are
+with fix_adapt are
 
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`class2 <bond_class2>`        | r0    | type bonds |
+| :doc:`class2 <bond_class2>`        | r0         | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`fene <bond_fene>`            | k, r0 | type bonds |
+| :doc:`fene <bond_fene>`            | k,r0       | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`gromos <bond_gromos>`        | k, r0 | type bonds |
+| :doc:`fene/nm <bond_fene>`         | k,r0       | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`harmonic <bond_harmonic>`    | k,r0  | type bonds |
+| :doc:`gromos <bond_gromos>`        | k,r0       | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`morse <bond_morse>`          | r0    | type bonds |
+| :doc:`harmonic <bond_harmonic>`    | k,r0       | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
-| :doc:`nonlinear <bond_nonlinear>`  | r0    | type bonds |
+| :doc:`morse <bond_morse>`          | r0         | type bonds |
-+------------------------------------+-------+------------+
++------------------------------------+-------+-----------------+
+| :doc:`nonlinear <bond_nonlinear>`  | epsilon,r0 | type bonds |
++------------------------------------+-------+-----------------+
+
+----------
+
+The *angle* keyword uses the specified variable to change the value of
+an angle coefficient over time, very similar to how the *pair* keyword
+operates. The only difference is that now an angle coefficient for a
+given angle type is adapted.
+
+A wild-card asterisk can be used in place of or in conjunction with
+the angle type argument to set the coefficients for multiple angle
+types.  This takes the form "\*" or "\*n" or "n\*" or "m\*n".  If N =
+the number of angle types, then an asterisk with no numeric values
+means all types from 1 to N.  A leading asterisk means all types from
+1 to n (inclusive).  A trailing asterisk means all types from n to N
+(inclusive).  A middle asterisk means all types from m to n
+(inclusive).
+
+Currently *angle* does not support angle_style hybrid nor angle_style
+hybrid/overlay as angle styles. The angle styles that currently work
+with fix_adapt are
+
++------------------------------------+-------+-----------------+
+| :doc:`harmonic <angle_harmonic>`    | k,theta0 | type angles |
++------------------------------------+-------+-----------------+
+| :doc:`cosine <angle_cosine>`        | k        | type angles |
++------------------------------------+-------+-----------------+
+
+Note that internally, theta0 is stored in radians, so the variable
+this fix uses to reset theta0 needs to generate values in radians.
 
 ----------
 

doc/src/fix_bond_break.rst

@@ -113,6 +113,9 @@ You can dump out snapshots of the current bond topology via the :doc:`dump local
    may need to thermostat your system to compensate for energy changes
    resulting from broken bonds (and angles, dihedrals, impropers).
 
+See the :doc:`Howto <Howto_broken_bonds>` page on broken bonds for more
+information on related features in LAMMPS.
+
 ----------
 
 Restart, fix_modify, output, run start/stop, minimize info

doc/src/fix_nve.rst

@@ -35,6 +35,10 @@ consistent with the microcanonical ensemble (NVE) provided there
 are (full) periodic boundary conditions and no other "manipulations"
 of the system (e.g. fixes that modify forces or velocities).
 
+This fix invokes the velocity form of the
+Stoermer-Verlet time integration algorithm (velocity-Verlet). Other
+time integration options can be invoked using the :doc:`run_style <run_style>` command.
+
 ----------
 
 .. include:: accel_styles.rst
@@ -57,7 +61,7 @@ Restrictions
 Related commands
 """"""""""""""""
 
-:doc:`fix nvt <fix_nh>`, :doc:`fix npt <fix_nh>`
+:doc:`fix nvt <fix_nh>`, :doc:`fix npt <fix_nh>`, :doc:`run_style <run_style>`
 
 Default
 """""""

doc/src/fix_nve_bpm_sphere.rst

@@ -0,0 +1,87 @@
+.. index:: fix nve/bpm/sphere
+
+fix nve/bpm/sphere command
+==========================
+
+Syntax
+""""""
+
+.. parsed-literal::
+
+   fix ID group-ID nve/bpm/sphere
+
+* ID, group-ID are documented in :doc:`fix <fix>` command
+* nve/bpm/sphere = style name of this fix command
+* zero or more keyword/value pairs may be appended
+* keyword = *disc*
+
+  .. parsed-literal::
+
+       *disc* value = none = treat particles as 2d discs, not spheres
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   fix 1 all nve/bpm/sphere
+   fix 1 all nve/bpm/sphere disc
+
+Description
+"""""""""""
+
+Perform constant NVE integration to update position, velocity, angular
+velocity, and quaternion orientation for finite-size spherical
+particles in the group each timestep.  V is volume; E is energy.  This
+creates a system trajectory consistent with the microcanonical
+ensemble.
+
+This fix differs from the :doc:`fix nve <fix_nve>` command, which
+assumes point particles and only updates their position and velocity.
+It also differs from the :doc:`fix nve/sphere <fix_nve_sphere>`
+command which assumes finite-size spheroid particles which do not
+store a quaternion.  It thus does not update a particle's orientation
+or quaternion.
+
+If the *disc* keyword is used, then each particle is treated as a 2d
+disc (circle) instead of as a sphere.  This is only possible for 2d
+simulations, as defined by the :doc:`dimension <dimension>` keyword.
+The only difference between discs and spheres in this context is their
+moment of inertia, as used in the time integration.
+
+----------
+
+Restart, fix_modify, output, run start/stop, minimize info
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+No information about this fix is written to :doc:`binary restart files
+<restart>`.  None of the :doc:`fix_modify <fix_modify>` options are
+relevant to this fix.  No global or per-atom quantities are stored by
+this fix for access by various :doc:`output commands <Howto_output>`.
+No parameter of this fix can be used with the *start/stop* keywords of
+the :doc:`run <run>` command.  This fix is not invoked during
+:doc:`energy minimization <minimize>`.
+
+Restrictions
+""""""""""""
+
+This fix requires that atoms store torque, angular velocity (omega), a
+radius, and a quaternion as defined by the :doc:`atom_style bpm/sphere
+<atom_style>` command.
+
+All particles in the group must be finite-size spheres with
+quaternions.  They cannot be point particles.
+
+Use of the *disc* keyword is only allowed for 2d simulations, as
+defined by the :doc:`dimension <dimension>` keyword.
+
+Related commands
+""""""""""""""""
+
+:doc:`fix nve <fix_nve>`, :doc:`fix nve/sphere <fix_nve_sphere>`
+
+Default
+"""""""
+
+none
+

doc/src/fix_pair_tracker.rst

@@ -1,124 +0,0 @@
-.. index:: fix pair/tracker
-
-fix pair/tracker command
-========================
-
-Syntax
-""""""
-
-.. parsed-literal::
-
-   fix ID group-ID pair/tracker N attribute1 attribute2 ... keyword values ...
-
-* ID, group-ID are documented in :doc:`fix <fix>` command
-* pair/tracker = style name of this fix command
-* N = prepare data for output every this many timesteps
-* one or more attributes may be appended
-
-  .. parsed-literal::
-
-       possible attributes = id1 id2 time/created time/broken time/total
-                             rmin rave x y z
-
-  .. parsed-literal::
-
-          id1, id2 = IDs of the 2 atoms in each pair interaction
-          time/created = the time that the 2 atoms began interacting
-          time/broken = the time that the 2 atoms stopped interacting
-          time/total = the total time the 2 atoms interacted
-          r/min = the minimum radial distance between the 2 atoms during the interaction
-          r/ave = the average radial distance between the 2 atoms during the interaction
-          x, y, z = the center of mass position of the 2 atoms when they stopped interacting
-
-* zero or more keyword/value pairs may be appended
-* keyword = *time/min* or *type/include*
-
-  .. parsed-literal::
-
-       *time/min* value = T
-         T = minimum interaction time
-       *type/include* value = arg1 arg2
-         arg = separate lists of types (see below)
-
-Examples
-""""""""
-
-.. code-block:: LAMMPS
-
-   fix 1 all pair/tracker 1000 id1 id2 time/min 100
-   fix 1 all pair/tracker 1000 time/created time/broken type/include 1 * type/include 2 3,4
-
-Description
-"""""""""""
-
-Tracks properties of pairwise interactions between two atoms and records data
-whenever the atoms move beyond the interaction cutoff.
-Must be used in conjunction with :doc:`pair tracker <pair_tracker>`.
-Data is accumulated over a span of *N* timesteps before being deleted.
-The number of datums generated, aggregated across all processors, equals
-the number of broken interactions. Interactions are only included if both
-atoms are included in the specified fix group. Additional filters can be
-applied using the *time/min* or *type/include* keywords described below.
-
-.. note::
-
-   For extremely long-lived interactions, the calculation of *r/ave* may not be
-   correct due to double overflow.
-
-The *time/min* keyword sets a minimum amount of time that an interaction must
-persist to be included. This setting can be used to censor short-lived interactions.
-The *type/include* keyword filters interactions based on the types of the two atoms.
-Data is only saved for interactions between atoms with types in the two lists.
-Each list consists of a series of type
-ranges separated by commas. The range can be specified as a
-single numeric value, or a wildcard asterisk can be used to specify a range
-of values.  This takes the form "\*" or "\*n" or "n\*" or "m\*n".  For
-example, if M = the number of atom types, then an asterisk with no numeric
-values means all types from 1 to M.  A leading asterisk means all types
-from 1 to n (inclusive).  A trailing asterisk means all types from n to M
-(inclusive).  A middle asterisk means all types from m to n (inclusive).
-Multiple *type/include* keywords may be added.
-
-----------
-
-Restart, fix_modify, run start/stop, minimize info
-"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
-
-No information about this fix is written to :doc:`binary restart files <restart>`.
-None of the :doc:`fix_modify <fix_modify>` options are
-relevant to this fix.
-No parameter of this fix can be used with the *start/stop* keywords of
-the :doc:`run <run>` command.
-
-Output info
-"""""""""""
-
-This compute calculates a local vector or local array depending on the
-number of input values.  The length of the vector or number of rows in
-the array is the number of recorded, lost interactions.  If a single input is
-specified, a local vector is produced.  If two or more inputs are
-specified, a local array is produced where the number of columns = the
-number of inputs.  The vector or array can be accessed by any command
-that uses local values from a compute as input.  See the :doc:`Howto output <Howto_output>` page for an overview of LAMMPS output
-options.
-
-The vector or array values will be doubles that correspond to the
-specified attribute.
-
-Restrictions
-""""""""""""
-
-Must be used in conjunction with :doc:`pair style tracker <pair_tracker>`.
-
-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>` page for more info.
-
-Related commands
-""""""""""""""""
-
-:doc:`pair tracker <pair_tracker>`
-
-Default
-"""""""
-
-none

doc/src/fix_polarize.rst

@@ -80,7 +80,7 @@ It is assumed that the particles cannot pass through the interface during the si
 so that its local dielectric constant value does not change.
 
 There are some example scripts for using these fixes
-with LAMMPS in the ``examples/PACKAGES/dielectric directory``. The README file
+with LAMMPS in the ``examples/PACKAGES/dielectric`` directory. The README file
 therein contains specific details on the system setup. Note that the example data files
 show the additional fields (columns) needed for :doc:`atom_style dielectric <atom_style>`
 beyond the conventional fields *id*, *mol*, *type*, *q*, *x*, *y*, and *z*.
@@ -91,7 +91,7 @@ For fix *polarize/bem/gmres* and fix *polarize/bem/icc* the induced
 charges of the atoms in the specified group, which are the vertices on
 the interface, are computed using the equation:
 
-..math::
+.. math::
 
   \sigma_b(\mathbf{s}) = \dfrac{1 - \bar{\epsilon}}{\bar{\epsilon}}
      \sigma_f(\mathbf{s}) - \epsilon_0 \dfrac{\Delta \epsilon}{\bar{\epsilon}}
@@ -154,6 +154,9 @@ if LAMMPS was built with that package, which requires that also the
 KSPACE package is installed.  See the :doc:`Build package
 <Build_package>` page for more info.
 
+Note that the *polarize/bem/gmres* and *polarize/bem/icc* fixes only support
+:doc:`units <units>` *lj*, *real*, *metal*, *si* and *nano* at the moment.
+
 
 Related commands
 """"""""""""""""

doc/src/fix_property_atom.rst

@@ -304,13 +304,15 @@ uninterrupted fashion.
 .. warning::
 
    When reading data from a restart file, this fix command has to be
-   specified **exactly** the same was in the input script that created
+   specified **after** the *read_restart* command and **exactly** the
-   the restart file. LAMMPS will only check whether a fix is of the
+   same was in the input script that created the restart file.  LAMMPS
-   same style and has the same fix ID and in case of a match will then
+   will only check whether a fix is of the same style and has the same
-   try to initialize the fix with the data stored in the binary
+   fix ID and in case of a match will then try to initialize the fix
-   restart file.  If the names and associated date types in the new
+   with the data stored in the binary restart file.  If the names and
-   fix property/atom command do not match the old one exactly, data
+   associated date types in the new fix property/atom command do not
-   can be corrupted or LAMMPS may crash.
+   match the old one exactly, data can be corrupted or LAMMPS may crash.
+   If the fix is specified **before** the *read_restart* command its
+   data will not be restored.
 
 None of the :doc:`fix_modify <fix_modify>` options are relevant to
 this fix.  No global or per-atom quantities are stored by this fix for

doc/src/fix_ttm.rst

@@ -214,10 +214,10 @@ generate an error. LAMMPS will check if a "UNITS:" tag is in the first
 line and stop with an error, if there is a mismatch with the current
 units used.
 
-..note::
+.. note::
 
   The electronic temperature at each grid point must be a non-zero
-  positive value, both initially, and as the temperature evovles over
+  positive value, both initially, and as the temperature evolves over
   time.  Thus you must use either the *set* or *infile* keyword or be
   restarting a simulation that used this fix previously.
 

doc/src/group.rst

@@ -258,11 +258,17 @@ assignment is made at the beginning of the minimization, but not
 during the iterations of the minimizer.
 
 The point in the timestep at which atoms are assigned to a dynamic
-group is after the initial stage of velocity Verlet time integration
+group is after interatomic forces have been computed, but before any
-has been performed, and before neighbor lists or forces are computed.
+fixes which alter forces or otherwise update the system have been
-This is the point in the timestep where atom positions have just
+invoked.  This means that atom positions have been updated, neighbor
-changed due to the time integration, so the region criterion should be
+lists and ghost atoms are current, and both intermolecular and
-accurate, if applied.
+intramolecular forces have been calculated based on the new
+coordinates.  Thus the region criterion, if applied, should be
+accurate.  Also, any computes invoked by an atom-style variable should
+use updated information for that timestep, e.g. potential energy/atom
+or coordination number/atom.  Similarly, fixes or computes which are
+invoked after that point in the timestep, should operate on the new
+group of atoms.
 
 .. note::
 

doc/src/pair_bpm_spring.rst

@@ -0,0 +1,114 @@
+.. index:: pair_style bpm/spring
+
+pair_style bpm/spring command
+=============================
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   pair_style bpm/spring
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   pair_style bpm/spring
+   pair_coeff * * 1.0 1.0 1.0
+   pair_coeff 1 1 1.0 1.0 1.0
+
+Description
+"""""""""""
+
+Style *bpm/spring* computes pairwise forces with the formula
+
+.. math::
+
+   F = k (r - r_c)
+
+where :math:`k` is a stiffness and :math:`r_c` is the cutoff length.
+An additional damping force is also applied to interacting
+particles. The force is proportional to the difference in the
+normal velocity of particles
+
+.. math::
+
+   F_D = - \gamma w (\hat{r} \bullet \vec{v})
+
+where :math:`\gamma` is the damping strength, :math:`\hat{r}` is the
+radial normal vector, :math:`\vec{v}` is the velocity difference
+between the two particles, and :math:`w` is a smoothing factor.
+This smoothing factor is constructed such that damping forces go to zero
+as particles come out of contact to avoid discontinuities. It is
+given by
+
+.. math::
+
+   w = 1.0 - \left( \frac{r}{r_c} \right)^8 .
+
+This pair style is designed for use in a spring-based bonded particle
+model.  It mirrors the construction of the :doc:`bpm/spring
+<bond_bpm_spring>` bond style.
+
+This pair interaction is always applied to pairs of non-bonded particles
+that are within the interaction distance. For pairs of bonded particles
+that are within the interaction distance, there is the option to either
+include this pair interaction and overlay the pair force over the bond
+force or to exclude this pair interaction such that the two particles
+only interact via the bond force. See discussion of the *overlay/pair*
+option for BPM bond styles and the :doc:`special_bonds <special_bonds>`
+command in the `:doc: how to <Howto_BPM>` page on BPMs for more details.
+
+The following coefficients must be defined for each pair of atom types
+via the :doc:`pair_coeff <pair_coeff>` command as in the examples
+above, or in the data file or restart files read by the
+:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
+commands, or by mixing as described below:
+
+* :math:`k`             (force/distance units)
+* :math:`r_c`           (distance units)
+* :math:`\gamma`        (force/velocity units)
+
+
+----------
+
+Mixing, shift, table, tail correction, restart, rRESPA info
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+For atom type pairs I,J and I != J, the A coefficient and cutoff
+distance for this pair style can be mixed.  A is always mixed via a
+*geometric* rule.  The cutoff is mixed according to the pair_modify
+mix value.  The default mix value is *geometric*\ .  See the
+"pair_modify" command for details.
+
+This pair style does not support the :doc:`pair_modify <pair_modify>`
+shift option, since the pair interaction goes to 0.0 at the cutoff.
+
+The :doc:`pair_modify <pair_modify>` table and tail options are not
+relevant for this pair style.
+
+This pair style writes its information to :doc:`binary restart files
+<restart>`, so pair_style and pair_coeff commands do not need to be
+specified in an input script that reads a restart file.
+
+This pair style can only be used via the *pair* keyword of the
+:doc:`run_style respa <run_style>` command.  It does not support the
+*inner*, *middle*, *outer* keywords.
+
+----------
+
+Restrictions
+""""""""""""
+ none
+
+Related commands
+""""""""""""""""
+
+:doc:`pair_coeff <pair_coeff>`, :doc:`bond bpm/spring <bond_bpm_spring>`
+
+Default
+"""""""
+
+none

doc/src/pair_dielectric.rst

@@ -3,6 +3,7 @@
 .. index:: pair_style lj/cut/coul/cut/dielectric
 .. index:: pair_style lj/cut/coul/cut/dielectric/omp
 .. index:: pair_style lj/cut/coul/debye/dielectric
+.. index:: pair_style lj/cut/coul/debye/dielectric/omp
 .. index:: pair_style lj/cut/coul/long/dielectric
 .. index:: pair_style lj/cut/coul/long/dielectric/omp
 .. index:: pair_style lj/cut/coul/msm/dielectric
@@ -22,6 +23,8 @@ Accelerator Variants: *lj/cut/coul/cut/dielectric/omp*
 pair_style lj/cut/coul/debye/dielectric command
 ===============================================
 
+Accelerator Variants: *lj/cut/coul/debye/dielectric/omp*
+
 pair_style lj/cut/coul/long/dielectric command
 ==============================================
 

doc/src/pair_ilp_graphene_hbn.rst

@@ -1,8 +1,11 @@
 .. index:: pair_style ilp/graphene/hbn
+.. index:: pair_style ilp/graphene/hbn/opt
 
 pair_style ilp/graphene/hbn command
 ===================================
 
+Accelerator Variant: *ilp/graphene/hbn/opt*
+
 Syntax
 """"""
 
@@ -125,6 +128,10 @@ headings) the following commands could be included in an input script:
 
 ----------
 
+.. include:: accel_styles.rst
+
+----------
+
 Mixing, shift, table, tail correction, restart, rRESPA info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 

doc/src/pair_lj_long.rst

@@ -238,7 +238,7 @@ none
 
 .. _Veld2:
 
-**(In 't Veld)** In 't Veld, Ismail, Grest, J Chem Phys (accepted) (2007).
+**(In 't Veld)** In 't Veld, Ismail, Grest, J Chem Phys, 127, 144711 (2007).
 
 .. _Jorgensen4:
 

doc/src/pair_style.rst

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

doc/src/pair_tracker.rst

@@ -8,89 +8,182 @@ Syntax
 
 .. code-block:: LAMMPS
 
-   pair_style tracker keyword
+   pair_style tracker fix_ID N keyword values attribute1 attribute2 ...
 
-* zero or more keyword/arg pairs may be appended
+* fix_ID = ID of associated internal fix to store data
-* keyword = *finite*
+* N = prepare data for output every this many timesteps
+* zero or more keywords may be appended
+* keyword = *finite* or *time/min* or *type/include*
 
   .. parsed-literal::
 
       *finite* value = none
          pair style uses atomic diameters to identify contacts
+      *time/min* value = T
+         T = minimum number of timesteps of interaction
+      *type/include* value = list1 list2
+         list1,list2 = separate lists of types (see below)
+
+* one or more attributes may be appended
+
+  .. parsed-literal::
+
+       possible attributes = id1 id2 time/created time/broken time/total
+                             r/min r/ave x y z
+
+  .. parsed-literal::
+
+          id1, id2 = IDs of the 2 atoms in each pair interaction
+          time/created = the timestep that the 2 atoms began interacting
+          time/broken = the timestep that the 2 atoms stopped interacting
+          time/total = the total number of timesteps the 2 atoms interacted
+          r/min = the minimum radial distance between the 2 atoms during the interaction (distance units)
+          r/ave = the average radial distance between the 2 atoms during the interaction (distance units)
+          x, y, z = the center of mass position of the 2 atoms when they stopped interacting (distance units)
 
 Examples
 """"""""
 
 .. code-block:: LAMMPS
 
-   pair_style hybrid/overlay tracker ...
+   pair_style hybrid/overlay tracker myfix 1000 id1 id2 type/include 1 * type/include 2 3,4  lj/cut 2.5
    pair_coeff 1 1 tracker 2.0
 
-   pair_style hybrid/overlay tracker finite ...
+   pair_style hybrid/overlay tracker myfix 1000 finite x y z time/min 100 granular
    pair_coeff * * tracker
 
-   fix 1 all pair/tracker 1000 time/created time/broken
+   dump 1 all local 1000 dump.local f_myfix[1] f_myfix[2] f_myfix[3]
-   dump 1 all local 1000 dump.local f_1[1] f_1[2]
    dump_modify 1 write_header no
 
 Description
 """""""""""
 
-Style *tracker* monitors information about pairwise interactions.
+Style *tracker* monitors information about pairwise interactions.  It
-It does not calculate any forces on atoms.
+does not calculate any forces on atoms. :doc:`Pair hybrid/overlay
-:doc:`Pair hybrid/overlay <pair_hybrid>` can be used to combine this pair
+<pair_hybrid>` can be used to combine this pair style with any other
-style with another pair style. Style *tracker*  must be used in conjunction
+pair style, as shown in the examples above.
-with about :doc:`fix pair_tracker <fix_pair_tracker>` which contains
-information on what data can be output.
 
-If the *finite* keyword is not defined, the following coefficients must be
+At each timestep, if two neighboring atoms move beyond the interaction
-defined for each pair of atom types via the :doc:`pair_coeff <pair_coeff>`
+cutoff, pairwise data is processed and transferred to an internal fix
-command as in the examples above, or in the data file or restart files
+labeled *fix_ID*. This allows the local data to be accessed by other
-read by the :doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
+LAMMPS commands. Additional
-commands, or by mixing as described below:
+filters can be applied using the *time/min* or *type/include* keywords
+described below.  Note that this is the interaction cutoff defined by
+this pair style, not the short-range cutoff defined by the pair style
+that is calculating forces on atoms.
+
+Following any optional keyword/value arguments, a list of one or more
+attributes is specified.  These include the IDs of the two atoms in
+the pair.  The other attributes for the pair of atoms are the
+duration of time they were "interacting" or at the point in time they
+started or stopped interacting.  In this context, "interacting" means
+the time window during which the two atoms were closer than the
+interaction cutoff distance.  The attributes for time/* refer to
+timesteps.
+
+Data is continuously accumulated by the internal fix over intervals of *N*
+timesteps. At the end of each interval, all of the saved accumulated
+data is deleted to make room for new data. Individual datum may
+therefore persist anywhere between *1* to *N* timesteps depending on
+when they are saved. This data can be accessed using the *fix_ID* and a
+:doc:`dump local <dump>` command. To ensure all data is output,
+the dump frequency should correspond to the same interval of *N*
+timesteps. A dump frequency of an integer multiple of *N* can be used
+to regularly output a sample of the accumulated data.
+
+----------
+
+The following optional keywords may be used.
+
+If the *finite* keyword is not used, the following coefficients must
+be defined for each pair of atom types via the :doc:`pair_coeff
+<pair_coeff>` command as in the examples above, or in the data file or
+restart files read by the :doc:`read_data <read_data>` or
+:doc:`read_restart <read_restart>` commands, or by mixing as described
+below:
 
 * cutoff (distance units)
 
-If the *finite* keyword is defined, no coefficients may be defined.
+If the *finite* keyword is used, there are no additional coefficients
-Interaction cutoffs are alternatively calculated based on the
+to set for each pair of atom types via the
-diameter of finite particles.
+:doc:`pair_coeff <pair_coeff>` command. Interaction cutoffs are
+instead calculated based on the diameter of finite particles. However
+you must still use the :doc:`pair_coeff <pair_coeff>` for all atom
+types. For example the command
 
+.. code-block:: LAMMPS
+
+   pair_coeff * *
+
+should be used.
+
+The *time/min* keyword sets a minimum amount of time that an
+interaction must persist to be included.  This setting can be used to
+censor short-lived interactions.
+
+The *type/include* keyword filters interactions based on the types of
+the two atoms.  Data is only saved for interactions between atoms
+whose two atom types appear in *list1* and *list2*.  Atom type 1 must
+be in list1 and atom type 2 in list2.  Or vice versa.
+
+Each type list consists of a series of type ranges separated by
+commas.  Each range can be specified as a single numeric value, or a
+wildcard asterisk can be used to specify a range of values.  This
+takes the form "\*" or "\*n" or "n\*" or "m\*n".  For example, if M =
+the number of atom types, then an asterisk with no numeric values
+means all types from 1 to M.  A leading asterisk means all types from
+1 to n (inclusive).  A trailing asterisk means all types from n to M
+(inclusive).  A middle asterisk means all types from m to n
+(inclusive).  Note that the *type/include* keyword can be specified
+multiple times.
 
 Mixing, shift, table, tail correction, restart, rRESPA info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
 For atom type pairs I,J and I != J, the cutoff coefficient and cutoff
-distance for this pair style can be mixed.  The cutoff is always mixed via a
+distance for this pair style can be mixed.  The cutoff is always mixed
-*geometric* rule.  The cutoff is mixed according to the pair_modify
+via a *geometric* rule.  The cutoff is mixed according to the
-mix value.  The default mix value is *geometric*\ .  See the
+pair_modify mix value.  The default mix value is *geometric*\ .  See
-"pair_modify" command for details.
+the "pair_modify" command for details.
 
-This pair style writes its information to :doc:`binary restart files <restart>`, so
+This pair style writes its information to :doc:`binary restart files
-pair_style and pair_coeff commands do not need
+<restart>`, so pair_style and pair_coeff commands do not need to be
-to be specified in an input script that reads a restart file.
+specified in an input script that reads a restart file.
 
 The :doc:`pair_modify <pair_modify>` shift, table, and tail options
 are not relevant for this pair style.
 
+The accumulated data is not written to restart files and should be
+output before a restart file is written to avoid missing data.
+
+The internal fix calculates a local vector or local array depending on the
+number of input values.  The length of the vector or number of rows in
+the array is the number of recorded, lost interactions.  If a single
+input is specified, a local vector is produced.  If two or more inputs
+are specified, a local array is produced where the number of columns =
+the number of inputs.  The vector or array can be accessed by any
+command that uses local values from a compute as input.  See the
+:doc:`Howto output <Howto_output>` page for an overview of LAMMPS
+output options.
+
+The vector or array will be floating point values that correspond to
+the specified attribute.
+
 ----------
 
 Restrictions
 """"""""""""
 
-A corresponding :doc:`fix pair_tracker <fix_pair_tracker>` must be defined
-to use this pair style.
-
-This pair style is currently incompatible with granular pair styles that extend
-beyond the contact (e.g. JKR and DMT).
-
 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>` page for more info.
+was built with that package.  See the :doc:`Build package
+<Build_package>` page for more info.
+
+This pair style is currently incompatible with granular pair styles
+that extend beyond the contact (e.g. JKR and DMT).
 
 Related commands
 """"""""""""""""
 
-:doc:`fix pair_tracker <fix_pair_tracker>`
-
 Default
 """""""
 

doc/src/read_data.rst

@@ -650,6 +650,8 @@ of analysis.
      - atom-ID atom-type rho esph cv x y z
    * - sphere
      - atom-ID atom-type diameter density x y z
+   * - bpm/sphere
+     - atom-ID molecule-ID atom-type diameter density x y z
    * - spin
      - atom-ID atom-type x y z spx spy spz sp
    * - tdpd

doc/src/run_style.rst

@@ -67,7 +67,8 @@ Description
 Choose the style of time integrator used for molecular dynamics
 simulations performed by LAMMPS.
 
-The *verlet* style is a standard velocity-Verlet integrator.
+The *verlet* style is the velocity form of the
+Stoermer-Verlet time integration algorithm (velocity-Verlet)
 
 ----------
 

doc/src/set.rst

@@ -22,7 +22,7 @@ Syntax
   *volume* or *image* or *bond* or *angle* or *dihedral* or
   *improper* or *sph/e* or *sph/cv* or *sph/rho* or
   *smd/contact/radius* or *smd/mass/density* or *dpd/theta* or
-  *edpd/temp* or *edpd/cv* or *cc* or
+  *edpd/temp* or *edpd/cv* or *cc* or *epsilon* or
   *i_name* or *d_name* or *i2_name* or *d2_name*
 
   .. parsed-literal::
@@ -124,6 +124,7 @@ Syntax
        *cc* values = index cc
          index = index of a chemical species (1 to Nspecies)
          cc = chemical concentration of tDPD particles for a species (mole/volume units)
+       *epsilon* value = dielectric constant of the medium where the atoms reside
        *i_name* value = custom integer vector with name
        *d_name* value = custom floating-point vector with name
        *i2_name* value = column of a custom integer array with name
@@ -372,13 +373,13 @@ vector of the particles is set to the 3 specified components.
 
 Keyword *omega* sets the angular velocity of selected atoms.  The
 particles must be spheres as defined by the :doc:`atom_style sphere
-<atom_style>` command.  The angular velocity vector of the particles
+<atom_style>` command.  The angular velocity vector of the particles is
-is set to the 3 specified components.
+set to the 3 specified components.
 
 Keyword *mass* sets the mass of all selected particles.  The particles
-must have a per-atom mass attribute, as defined by the
+must have a per-atom mass attribute, as defined by the :doc:`atom_style
-:doc:`atom_style <atom_style>` command.  See the "mass" command for
+<atom_style>` command.  See the "mass" command for how to set mass
-how to set mass values on a per-type basis.
+values on a per-type basis.
 
 Keyword *density* or *density/disc* also sets the mass of all selected
 particles, but in a different way.  The particles must have a per-atom
@@ -387,22 +388,21 @@ command.  If the atom has a radius attribute (see :doc:`atom_style
 sphere <atom_style>`) and its radius is non-zero, its mass is set from
 the density and particle volume for 3d systems (the input density is
 assumed to be in mass/distance\^3 units).  For 2d, the default is for
-LAMMPS to model particles with a radius attribute as spheres.
+LAMMPS to model particles with a radius attribute as spheres.  However,
-However, if the *density/disc* keyword is used, then they can be
+if the *density/disc* keyword is used, then they can be modeled as 2d
-modeled as 2d discs (circles).  Their mass is set from the density and
+discs (circles).  Their mass is set from the density and particle area
-particle area (the input density is assumed to be in mass/distance\^2
+(the input density is assumed to be in mass/distance\^2 units).
-units).
 
 If the atom has a shape attribute (see :doc:`atom_style ellipsoid
-<atom_style>`) and its 3 shape parameters are non-zero, then its mass
+<atom_style>`) and its 3 shape parameters are non-zero, then its mass is
-is set from the density and particle volume (the input density is
+set from the density and particle volume (the input density is assumed
-assumed to be in mass/distance\^3 units).  The *density/disc* keyword
+to be in mass/distance\^3 units).  The *density/disc* keyword has no
-has no effect; it does not (yet) treat 3d ellipsoids as 2d ellipses.
+effect; it does not (yet) treat 3d ellipsoids as 2d ellipses.
 
 If the atom has a length attribute (see :doc:`atom_style line
-<atom_style>`) and its length is non-zero, then its mass is set from
+<atom_style>`) and its length is non-zero, then its mass is set from the
-the density and line segment length (the input density is assumed to
+density and line segment length (the input density is assumed to be in
-be in mass/distance units).  If the atom has an area attribute (see
+mass/distance units).  If the atom has an area attribute (see
 :doc:`atom_style tri <atom_style>`) and its area is non-zero, then its
 mass is set from the density and triangle area (the input density is
 assumed to be in mass/distance\^2 units).
@@ -410,84 +410,91 @@ assumed to be in mass/distance\^2 units).
 If none of these cases are valid, then the mass is set to the density
 value directly (the input density is assumed to be in mass units).
 
-Keyword *volume* sets the volume of all selected particles.
+Keyword *volume* sets the volume of all selected particles.  Currently,
-Currently, only the :doc:`atom_style peri <atom_style>` command defines
+only the :doc:`atom_style peri <atom_style>` command defines particles
-particles with a volume attribute.  Note that this command does not
+with a volume attribute.  Note that this command does not adjust the
-adjust the particle mass.
+particle mass.
 
 Keyword *image* sets which image of the simulation box the atom is
 considered to be in.  An image of 0 means it is inside the box as
-defined.  A value of 2 means add 2 box lengths to get the true value.
+defined.  A value of 2 means add 2 box lengths to get the true value.  A
-A value of -1 means subtract 1 box length to get the true value.
+value of -1 means subtract 1 box length to get the true value.  LAMMPS
-LAMMPS updates these flags as atoms cross periodic boundaries during
+updates these flags as atoms cross periodic boundaries during the
-the simulation.  The flags can be output with atom snapshots via the
+simulation.  The flags can be output with atom snapshots via the
-:doc:`dump <dump>` command.  If a value of NULL is specified for any
+:doc:`dump <dump>` command.  If a value of NULL is specified for any of
-of nx,ny,nz, then the current image value for that dimension is
+nx,ny,nz, then the current image value for that dimension is unchanged.
-unchanged.  For non-periodic dimensions only a value of 0 can be
+For non-periodic dimensions only a value of 0 can be specified.  This
-specified.  This command can be useful after a system has been
+command can be useful after a system has been equilibrated and atoms
-equilibrated and atoms have diffused one or more box lengths in
+have diffused one or more box lengths in various directions.  This
-various directions.  This command can then reset the image values for
+command can then reset the image values for atoms so that they are
-atoms so that they are effectively inside the simulation box, e.g if a
+effectively inside the simulation box, e.g if a diffusion coefficient is
-diffusion coefficient is about to be measured via the :doc:`compute
+about to be measured via the :doc:`compute msd <compute_msd>` command.
-msd <compute_msd>` command.  Care should be taken not to reset the
+Care should be taken not to reset the image flags of two atoms in a bond
-image flags of two atoms in a bond to the same value if the bond
+to the same value if the bond straddles a periodic boundary (rather they
-straddles a periodic boundary (rather they should be different by +/-
+should be different by +/- 1).  This will not affect the dynamics of a
-1).  This will not affect the dynamics of a simulation, but may mess
+simulation, but may mess up analysis of the trajectories if a LAMMPS
-up analysis of the trajectories if a LAMMPS diagnostic or your own
+diagnostic or your own analysis relies on the image flags to unwrap a
-analysis relies on the image flags to unwrap a molecule which
+molecule which straddles the periodic box.
-straddles the periodic box.
 
-Keywords *bond*, *angle*, *dihedral*, and *improper*, set the bond
+Keywords *bond*, *angle*, *dihedral*, and *improper*, set the bond type
-type (angle type, etc) of all bonds (angles, etc) of selected atoms to
+(angle type, etc) of all bonds (angles, etc) of selected atoms to the
-the specified value from 1 to nbondtypes (nangletypes, etc).  All
+specified value from 1 to nbondtypes (nangletypes, etc).  All atoms in a
-atoms in a particular bond (angle, etc) must be selected atoms in
+particular bond (angle, etc) must be selected atoms in order for the
-order for the change to be made.  The value of nbondtype (nangletypes,
+change to be made.  The value of nbondtype (nangletypes, etc) was set by
-etc) was set by the *bond types* (\ *angle types*, etc) field in the
+the *bond types* (\ *angle types*, etc) field in the header of the data
-header of the data file read by the :doc:`read_data <read_data>`
+file read by the :doc:`read_data <read_data>` command.  These keywords
-command.  These keywords do not allow use of an atom-style variable.
+do not allow use of an atom-style variable.
 
-Keywords *sph/e*, *sph/cv*, and *sph/rho* set the energy, heat
+Keywords *sph/e*, *sph/cv*, and *sph/rho* set the energy, heat capacity,
-capacity, and density of smoothed particle hydrodynamics (SPH)
+and density of smoothed particle hydrodynamics (SPH) particles.  See
-particles.  See `this PDF guide <PDF/SPH_LAMMPS_userguide.pdf>`_
+`this PDF guide <PDF/SPH_LAMMPS_userguide.pdf>`_ to using SPH in LAMMPS.
-to using SPH in LAMMPS.
 
-Keyword *smd/mass/density* sets the mass of all selected particles,
+Keyword *smd/mass/density* sets the mass of all selected particles, but
-but it is only applicable to the Smooth Mach Dynamics package
+it is only applicable to the Smooth Mach Dynamics package MACHDYN.  It
-MACHDYN.  It assumes that the particle volume has already been
+assumes that the particle volume has already been correctly set and
-correctly set and calculates particle mass from the provided mass
+calculates particle mass from the provided mass density value.
-density value.
 
-Keyword *smd/contact/radius* only applies to simulations with the
+Keyword *smd/contact/radius* only applies to simulations with the Smooth
-Smooth Mach Dynamics package MACHDYN.  Itsets an interaction radius
+Mach Dynamics package MACHDYN.  Itsets an interaction radius for
-for computing short-range interactions, e.g. repulsive forces to
+computing short-range interactions, e.g. repulsive forces to prevent
-prevent different individual physical bodies from penetrating each
+different individual physical bodies from penetrating each other. Note
-other. Note that the SPH smoothing kernel diameter used for computing
+that the SPH smoothing kernel diameter used for computing long range,
-long range, nonlocal interactions, is set using the *diameter*
+nonlocal interactions, is set using the *diameter* keyword.
-keyword.
 
 Keyword *dpd/theta* sets the internal temperature of a DPD particle as
-defined by the DPD-REACT package.  If the specified value is a number
+defined by the DPD-REACT package.  If the specified value is a number it
-it must be >= 0.0.  If the specified value is NULL, then the kinetic
+must be >= 0.0.  If the specified value is NULL, then the kinetic
-temperature Tkin of each particle is computed as 3/2 k Tkin = KE = 1/2
+temperature Tkin of each particle is computed as 3/2 k Tkin = KE = 1/2 m
-m v\^2 = 1/2 m (vx\*vx+vy\*vy+vz\*vz).  Each particle's internal
+v\^2 = 1/2 m (vx\*vx+vy\*vy+vz\*vz).  Each particle's internal
 temperature is set to Tkin.  If the specified value is an atom-style
-variable, then the variable is evaluated for each particle.  If a
+variable, then the variable is evaluated for each particle.  If a value
-value >= 0.0, the internal temperature is set to that value.  If it is
+>= 0.0, the internal temperature is set to that value.  If it is < 0.0,
-< 0.0, the computation of Tkin is performed and the internal
+the computation of Tkin is performed and the internal temperature is set
-temperature is set to that value.
+to that value.
 
 Keywords *edpd/temp* and *edpd/cv* set the temperature and volumetric
 heat capacity of an eDPD particle as defined by the DPD-MESO package.
 Currently, only :doc:`atom_style edpd <atom_style>` defines particles
-with these attributes. The values for the temperature and heat
+with these attributes. The values for the temperature and heat capacity
-capacity must be positive.
+must be positive.
 
 Keyword *cc* sets the chemical concentration of a tDPD particle for a
 specified species as defined by the DPD-MESO package. Currently, only
 :doc:`atom_style tdpd <atom_style>` defines particles with this
 attribute. An integer for "index" selects a chemical species (1 to
-Nspecies) where Nspecies is set by the atom_style command. The value
+Nspecies) where Nspecies is set by the atom_style command. The value for
-for the chemical concentration must be >= 0.0.
+the chemical concentration must be >= 0.0.
+
+Keyword *epsilon* sets the dielectric constant of a particle, precisely
+of the medium where the particle resides as defined by the DIELECTRIC
+package. Currently, only :doc:`atom_style dielectric <atom_style>`
+defines particles with this attribute. The value for the dielectric
+constant must be >= 0.0.  Note that the set command with this keyword
+will rescale the particle charge accordingly so that the real charge
+(e.g., as read from a data file) stays intact. To change the real
+charges, one needs to use the set command with the *charge*
+keyword. Care must be taken to ensure that the real and scaled charges,
+and dielectric constants are consistent.
 
 Keywords *i_name*, *d_name*, *i2_name*, *d2_name* refer to custom
 per-atom integer and floating-point vectors or arrays that have been

doc/src/thermo_modify.rst

@@ -252,6 +252,6 @@ flush = no, and temp/press = compute IDs defined by thermo_style.
 
 The defaults for the line and format options depend on the thermo style.
 For styles "one" and "custom", the line and format defaults are "one",
-"%10d", and "%12.8g".  For style "multi", the line and format defaults
+"%10d", and "%14.8g".  For style "multi", the line and format defaults
 are "multi", "%14d", and "%14.4f". For style "yaml", the line and format
 defaults are "%d" and "%.15g".

doc/src/thermo_style.rst

@@ -10,7 +10,7 @@ Syntax
 
    thermo_style style args
 
-* style = *one* or *multi* *yaml* or *custom*
+* style = *one* or *multi* or *yaml* or *custom*
 * args = list of arguments for a particular style
 
   .. parsed-literal::

doc/utils/requirements.txt

@@ -6,3 +6,4 @@ breathe
 Pygments
 six
 pyyaml
+wheel

doc/utils/sphinx-config/false_positives.txt

@@ -316,6 +316,7 @@ borophene
 Botero
 Botu
 Bouguet
+Bourasseau
 Bourne
 boxcolor
 boxhi
@@ -329,6 +330,8 @@ boxzlo
 bp
 bpclermont
 bpls
+bpm
+BPM
 br
 Branduardi
 Branicio
@@ -341,6 +344,7 @@ Broglie
 brownian
 brownw
 Broyden
+Brusselle
 Bryantsev
 Btarget
 btype
@@ -676,6 +680,7 @@ Deresiewicz
 Derjagin
 Derjaguin
 Derlet
+Desbiens
 Deserno
 Destree
 destructor
@@ -789,6 +794,7 @@ Dullweber
 dumpfile
 Dunbrack
 Dunweg
+Dupend
 Dupont
 dUs
 dV
@@ -1155,8 +1161,10 @@ gdot
 GeC
 Geier
 gencode
+Geocomputing
 georg
 Georg
+Geotechnica
 germain
 Germann
 Germano
@@ -1326,6 +1334,7 @@ Holm
 holonomic
 Homebrew
 hooke
+hookean
 Hookean
 hostname
 hotpink
@@ -1671,6 +1680,7 @@ Kusters
 Kutta
 Kuznetsov
 kx
+Lachet
 Lackmann
 Ladd
 lagrangian
@@ -2082,6 +2092,7 @@ monopole
 monovalent
 Montalenti
 Montero
+Mora
 Morefoo
 Morfill
 Mori
@@ -2211,6 +2222,7 @@ Nbin
 Nbins
 nbody
 Nbody
+nbond
 nbonds
 nbondtype
 Nbondtype
@@ -3185,12 +3197,14 @@ Steinhauser
 Stepaniants
 stepwise
 Stesmans
+stiffnesses
 Stillinger
 stk
 stochastically
 stochasticity
 Stockmayer
 Stoddard
+Stoermer
 stoichiometric
 stoichiometry
 Stokesian
@@ -3602,6 +3616,7 @@ Voronoi
 VORONOI
 Vorselaars
 Voth
+Voyiatzis
 vpz
 vratio
 Vries
@@ -3668,6 +3683,7 @@ wn
 Wolde
 workflow
 workflows
+Workum
 Worley
 Wriggers
 Wuppertal

examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/born_matrix.out

@@ -0,0 +1,7 @@
+# Cij Matrix from post process computation
+ 3.36316  1.87373  1.87607 -0.00346 -0.00172 -0.00104
+ 1.87373  3.36170  1.87425  0.00443  0.00033  0.00014
+ 1.87607  1.87425  3.36573  0.00143  0.00155  0.00127
+-0.00346  0.00443  0.00143  1.87425  0.00127  0.00033
+-0.00172  0.00033  0.00155  0.00127  1.87607 -0.00346
+-0.00104  0.00014  0.00127  0.00033 -0.00346  1.87373

examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/compute_born.py

@@ -0,0 +1,118 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import sys
+import numpy as np
+
+def reduce_Born(Cf):
+    C = np.zeros((6,6), dtype=np.float64)
+    C[0,0] = Cf[0]
+    C[1,1] = Cf[1]
+    C[2,2] = Cf[2]
+    C[3,3] = Cf[3]
+    C[4,4] = Cf[4]
+    C[5,5] = Cf[5]
+    C[0,1] = Cf[6]
+    C[0,2] = Cf[7]
+    C[0,3] = Cf[8]
+    C[0,4] = Cf[9]
+    C[0,5] = Cf[10]
+    C[1,2] = Cf[11]
+    C[1,3] = Cf[12]
+    C[1,4] = Cf[13]
+    C[1,5] = Cf[14]
+    C[2,3] = Cf[15]
+    C[2,4] = Cf[16]
+    C[2,5] = Cf[17]
+    C[3,4] = Cf[18]
+    C[3,5] = Cf[19]
+    C[4,5] = Cf[20]
+    C = np.where(C,C,C.T)
+    return C
+
+def compute_delta():
+    D = np.zeros((3,3,3,3))
+    for a in range(3):
+        for b in range(3):
+            for m in range(3):
+                for n in range(3):
+                    D[a,b,m,n] = (a==m)*(b==n) + (a==n)*(b==m)
+    d = np.zeros((6,6))
+    d[0,0] = D[0,0,0,0]
+    d[1,1] = D[1,1,1,1]
+    d[2,2] = D[2,2,2,2]
+    d[3,3] = D[1,2,1,2]
+    d[4,4] = D[0,2,0,2]
+    d[5,5] = D[0,1,0,1]
+    d[0,1] = D[0,0,1,1]
+    d[0,2] = D[0,0,2,2]
+    d[0,3] = D[0,0,1,2]
+    d[0,4] = D[0,0,0,2]
+    d[0,5] = D[0,0,0,1]
+    d[1,2] = D[1,1,2,2]
+    d[1,3] = D[1,1,1,2]
+    d[1,4] = D[1,1,0,2]
+    d[1,5] = D[1,1,0,1]
+    d[2,3] = D[2,2,1,2]
+    d[2,4] = D[2,2,0,2]
+    d[2,5] = D[2,2,0,1]
+    d[3,4] = D[1,2,0,2]
+    d[3,5] = D[1,2,0,1]
+    d[4,5] = D[0,2,0,1]
+    d = np.where(d,d,d.T)
+    return d
+
+
+def write_matrix(C, filename):
+    with open(filename, 'w') as f:
+        f.write("# Cij Matrix from post process computation\n")
+        for i in C:
+            f.write("{:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f}\n".format(
+                i[0]*10**-9, i[1]*10**-9, i[2]*10**-9, i[3]*10**-9, i[4]*10**-9, i[5]*10**-9,
+                )
+            )
+    return
+
+def main():
+
+    N = 500
+    vol = 27.047271**3 * 10**-30  # m^3
+    T = 60  # K
+    kb = 1.380649 * 10**-23  # J/K
+    kbT = T*kb  # J
+    kcalmol2J = 4183.9954/(6.022*10**23)
+
+    born = np.loadtxt('born.out')
+    stre = np.loadtxt('vir.out')
+    stre[:, 1:] = -stre[:, 1:]*101325  # -> Pa
+    try:
+        mean_born = np.mean(born[:, 1:], axis=0)
+    except IndexError:
+        mean_born = born[1:]
+
+    CB = kcalmol2J/vol*reduce_Born(mean_born)  # -> J/m^3=Pa
+    Cs = vol/kbT*np.cov(stre[:,1:].T)
+    Ct = N*kbT/vol * compute_delta()
+    C = CB - Cs + Ct
+    write_matrix(CB, 'born_matrix.out')
+    write_matrix(Cs, 'stre_matrix.out')
+    write_matrix(Ct, 'temp_matrix.out')
+    write_matrix(C, 'full_matrix.out')
+    C11 = np.mean([C[0,0], C[1,1], C[2,2]]) * 10**-9
+    C12 = np.mean([C[0,1], C[0,2], C[1,2]]) * 10**-9
+    C44 = np.mean([C[3,3], C[4,4], C[5,5]]) * 10**-9
+    eC11 = np.std([C[0,0], C[1,1], C[2,2]]) * 10**-9
+    eC12 = np.std([C[0,1], C[0,2], C[1,2]]) * 10**-9
+    eC44 = np.std([C[3,3], C[4,4], C[5,5]]) * 10**-9
+    print(C*10**-9)
+    print("C11 = {:f} ± {:f}; C12 = {:f} ± {:f}; C44 = {:f} ± {:f}".format(C11, eC11, C12, eC12, C44, eC44))
+
+    return
+
+
+if __name__ == "__main__":
+    try:
+        main()
+    except KeyboardInterrupt:
+        raise SystemExit("User interruption.")
+

examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/full_matrix.out

@@ -0,0 +1,7 @@
+# Cij Matrix from post process computation
+ 2.18161  1.13726  1.16596 -0.01607 -0.02637  0.00291
+ 1.13726  2.20242  1.16714  0.00386 -0.05820  0.02644
+ 1.16596  1.16714  2.24704 -0.00354 -0.00368  0.02714
+-0.01607  0.00386 -0.00354  1.43706  0.00210  0.01003
+-0.02637 -0.05820 -0.00368  0.00210  1.37530  0.01401
+ 0.00291  0.02644  0.02714  0.01003  0.01401  1.42403

examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/in.ljcov

@@ -0,0 +1,154 @@
+# Analytical calculation
+# of Born matrix
+
+# Note that because of cubic symmetry and central forces, we have:
+# C11, pure axial == positive mean value: 1,2,3
+# C44==C23, pure shear == positive mean value, exactly match in pairs: (4,12),(5,8),(6,7)
+# C14==C56, shear/axial(normal) == zero mean, exactly match in pairs: (9,21),(14,20),(18,19)
+# C15, shear/axial(in-plane) == zero mean: 10,11,13,15,16,17
+
+# adjustable parameters
+
+units           real
+variable        nsteps index 100000    # length of run
+variable        nthermo index  1000    # thermo output interval 
+variable        nlat equal 5           # size of box
+variable        T    equal 60.         # Temperature in K
+variable        rho  equal 5.405       # Lattice spacing in A
+
+atom_style      atomic
+
+lattice         fcc ${rho}
+region          box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
+create_box      1 box
+create_atoms    1 box
+
+mass            * 39.948
+
+velocity        all create ${T} 87287 loop geom
+velocity        all zero linear
+
+pair_style      lj/cut 12.0
+pair_coeff      1 1 0.238067 3.405
+
+neighbor        0.0 bin
+neigh_modify    every 1 delay 0 check no
+
+variable vol equal vol
+thermo 100
+fix aL all ave/time 1 1 1 v_vol ave running
+fix NPT all npt temp $T $T 100 aniso 1. 1. 1000 fixedpoint 0. 0. 0.
+
+run 20000
+
+unfix NPT
+
+variable newL equal "f_aL^(1./3.)"
+change_box all x final 0 ${newL} y final 0. ${newL} z final 0. ${newL} remap units box
+
+unfix aL
+
+reset_timestep 0
+
+# Conversion variables
+variable kb        equal 1.38065e-23    # J/K
+variable Myvol     equal "vol*10^-30" # Volume in m^3
+variable kbt       equal "v_kb*v_T"
+variable Nat       equal atoms
+variable Rhokbt    equal "v_kbt*v_Nat/v_Myvol"
+variable at2Pa     equal 101325
+variable kcalmol2J equal "4183.9954/(6.022e23)"
+variable C1        equal "v_kcalmol2J/v_Myvol" # Convert Cb from energy to pressure units
+variable C2        equal "v_Myvol/v_kbt"       # Factor for Cfl terms
+variable Pa2GPa    equal 1e-9
+
+# Born compute giving <C^b> terms
+compute     born all born/matrix
+# The six virial stress component to compute <C^fl>
+compute     VIR all pressure NULL virial
+variable s1 equal "-c_VIR[1]*v_at2Pa"
+variable s2 equal "-c_VIR[2]*v_at2Pa"
+variable s3 equal "-c_VIR[3]*v_at2Pa"
+variable s6 equal "-c_VIR[4]*v_at2Pa"
+variable s5 equal "-c_VIR[5]*v_at2Pa"
+variable s4 equal "-c_VIR[6]*v_at2Pa"
+variable press equal press
+
+
+# Average of Born term and vector to store stress
+# for post processing
+fix CB all ave/time 1 ${nthermo} ${nthermo} c_born[*] ave running file born.out overwrite
+fix CPR all ave/time 1 1 1 c_VIR[*] file vir.out
+fix APR all ave/time 1 1 1 v_press ave running
+fix VEC all vector 1 v_s1 v_s2 v_s3 v_s4 v_s5 v_s6
+
+thermo      ${nthermo}
+thermo_style    custom step temp press f_APR c_born[1] f_CB[1] c_born[12] f_CB[12] c_born[4] f_CB[4]
+thermo_modify line multi
+
+fix     1 all nvt temp $T $T 100
+
+run         ${nsteps}
+
+# Compute vector averages
+# Note the indice switch.
+# LAMMPS convention is NOT the Voigt notation.
+variable aves1 equal "ave(f_VEC[1])"
+variable aves2 equal "ave(f_VEC[2])"
+variable aves3 equal "ave(f_VEC[3])"
+variable aves4 equal "ave(f_VEC[6])"
+variable aves5 equal "ave(f_VEC[5])"
+variable aves6 equal "ave(f_VEC[4])"
+
+# Computing the covariance through the <s_{i}s_{j}>-<s_i><s_j>
+# is numerically instable. Here we go through the <(s-<s>)^2>
+# definition.
+
+# Computing difference relative to average values
+variable ds1 vector "f_VEC[1]-v_aves1"
+variable ds2 vector "f_VEC[2]-v_aves2"
+variable ds3 vector "f_VEC[3]-v_aves3"
+variable ds4 vector "f_VEC[4]-v_aves4"
+variable ds5 vector "f_VEC[5]-v_aves5"
+variable ds6 vector "f_VEC[6]-v_aves6"
+
+# Squaring and averaging
+variable dds1 vector "v_ds1*v_ds1"
+variable dds2 vector "v_ds2*v_ds2"
+variable dds3 vector "v_ds3*v_ds3"
+variable vars1 equal "ave(v_dds1)"
+variable vars2 equal "ave(v_dds2)"
+variable vars3 equal "ave(v_dds3)"
+variable C11 equal "v_Pa2GPa*(v_C1*f_CB[1] - v_C2*v_vars1 + 2*v_Rhokbt)"
+variable C22 equal "v_Pa2GPa*(v_C1*f_CB[2] - v_C2*v_vars2 + 2*v_Rhokbt)"
+variable C33 equal "v_Pa2GPa*(v_C1*f_CB[3] - v_C2*v_vars3 + 2*v_Rhokbt)"
+
+variable dds12 vector "v_ds1*v_ds2"
+variable dds13 vector "v_ds1*v_ds3"
+variable dds23 vector "v_ds2*v_ds3"
+variable vars12 equal "ave(v_dds12)"
+variable vars13 equal "ave(v_dds13)"
+variable vars23 equal "ave(v_dds23)"
+variable C12 equal "v_Pa2GPa*(v_C1*f_CB[7] - v_C2*v_vars12)"
+variable C13 equal "v_Pa2GPa*(v_C1*f_CB[8] - v_C2*v_vars13)"
+variable C23 equal "v_Pa2GPa*(v_C1*f_CB[12] - v_C2*v_vars23)"
+
+variable dds4 vector "v_ds4*v_ds4"
+variable dds5 vector "v_ds5*v_ds5"
+variable dds6 vector "v_ds6*v_ds6"
+variable vars4 equal "ave(v_dds4)"
+variable vars5 equal "ave(v_dds5)"
+variable vars6 equal "ave(v_dds6)"
+variable C44 equal "v_Pa2GPa*(v_C1*f_CB[4] - v_C2*v_vars4 + v_Rhokbt)"
+variable C55 equal "v_Pa2GPa*(v_C1*f_CB[5] - v_C2*v_vars5 + v_Rhokbt)"
+variable C66 equal "v_Pa2GPa*(v_C1*f_CB[6] - v_C2*v_vars6 + v_Rhokbt)"
+
+variable aC11 equal "(v_C11 + v_C22 + v_C33)/3."
+variable aC12 equal "(v_C12 + v_C13 + v_C23)/3."
+variable aC44 equal "(v_C44 + v_C55 + v_C66)/3."
+
+print """
+C11 = ${aC11}
+C12 = ${aC12}
+C44 = ${aC44}
+"""

examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/stre_matrix.out

@@ -0,0 +1,7 @@
+# Cij Matrix from post process computation
+ 1.22342  0.73647  0.71011  0.01261  0.02465 -0.00395
+ 0.73647  1.20115  0.70711  0.00057  0.05854 -0.02630
+ 0.71011  0.70711  1.16055  0.00497  0.00524 -0.02587
+ 0.01261  0.00057  0.00497  0.45813 -0.00083 -0.00969
+ 0.02465  0.05854  0.00524 -0.00083  0.52170 -0.01747
+-0.00395 -0.02630 -0.02587 -0.00969 -0.01747  0.47064

examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/temp_matrix.out

@@ -0,0 +1,7 @@
+# Cij Matrix from post process computation
+ 0.04187  0.00000  0.00000  0.00000  0.00000  0.00000
+ 0.00000  0.04187  0.00000  0.00000  0.00000  0.00000
+ 0.00000  0.00000  0.04187  0.00000  0.00000  0.00000
+ 0.00000  0.00000  0.00000  0.02093  0.00000  0.00000
+ 0.00000  0.00000  0.00000  0.00000  0.02093  0.00000
+ 0.00000  0.00000  0.00000  0.00000  0.00000  0.02093

examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/born_matrix.out

@@ -0,0 +1,7 @@
+# Cij Matrix from post process computation
+ 3.35855  1.86892  1.87139  0.00233  0.00218 -0.00179
+ 1.86892  3.37104  1.87285  0.00112  0.00085 -0.00007
+ 1.87139  1.87285  3.37707 -0.00058  0.00038 -0.00057
+ 0.00233  0.00112 -0.00058  1.88326 -0.00039  0.00065
+ 0.00218  0.00085  0.00038 -0.00039  1.88229  0.00242
+-0.00179 -0.00007 -0.00057  0.00065  0.00242  1.87968

examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/compute_born.py

@@ -0,0 +1,118 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import sys
+import numpy as np
+
+def reduce_Born(Cf):
+    C = np.zeros((6,6), dtype=np.float64)
+    C[0,0] = Cf[0]
+    C[1,1] = Cf[1]
+    C[2,2] = Cf[2]
+    C[3,3] = Cf[3]
+    C[4,4] = Cf[4]
+    C[5,5] = Cf[5]
+    C[0,1] = Cf[6]
+    C[0,2] = Cf[7]
+    C[0,3] = Cf[8]
+    C[0,4] = Cf[9]
+    C[0,5] = Cf[10]
+    C[1,2] = Cf[11]
+    C[1,3] = Cf[12]
+    C[1,4] = Cf[13]
+    C[1,5] = Cf[14]
+    C[2,3] = Cf[15]
+    C[2,4] = Cf[16]
+    C[2,5] = Cf[17]
+    C[3,4] = Cf[18]
+    C[3,5] = Cf[19]
+    C[4,5] = Cf[20]
+    C = np.where(C,C,C.T)
+    return C
+
+def compute_delta():
+    D = np.zeros((3,3,3,3))
+    for a in range(3):
+        for b in range(3):
+            for m in range(3):
+                for n in range(3):
+                    D[a,b,m,n] = (a==m)*(b==n) + (a==n)*(b==m)
+    d = np.zeros((6,6))
+    d[0,0] = D[0,0,0,0]
+    d[1,1] = D[1,1,1,1]
+    d[2,2] = D[2,2,2,2]
+    d[3,3] = D[1,2,1,2]
+    d[4,4] = D[0,2,0,2]
+    d[5,5] = D[0,1,0,1]
+    d[0,1] = D[0,0,1,1]
+    d[0,2] = D[0,0,2,2]
+    d[0,3] = D[0,0,1,2]
+    d[0,4] = D[0,0,0,2]
+    d[0,5] = D[0,0,0,1]
+    d[1,2] = D[1,1,2,2]
+    d[1,3] = D[1,1,1,2]
+    d[1,4] = D[1,1,0,2]
+    d[1,5] = D[1,1,0,1]
+    d[2,3] = D[2,2,1,2]
+    d[2,4] = D[2,2,0,2]
+    d[2,5] = D[2,2,0,1]
+    d[3,4] = D[1,2,0,2]
+    d[3,5] = D[1,2,0,1]
+    d[4,5] = D[0,2,0,1]
+    d = np.where(d,d,d.T)
+    return d
+
+
+def write_matrix(C, filename):
+    with open(filename, 'w') as f:
+        f.write("# Cij Matrix from post process computation\n")
+        for i in C:
+            f.write("{:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f}\n".format(
+                i[0]*10**-9, i[1]*10**-9, i[2]*10**-9, i[3]*10**-9, i[4]*10**-9, i[5]*10**-9,
+                )
+            )
+    return
+
+def main():
+
+    N = 500
+    vol = 27.047271**3 * 10**-30  # m^3
+    T = 60  # K
+    kb = 1.380649 * 10**-23  # J/K
+    kbT = T*kb  # J
+    kcalmol2J = 4183.9954/(6.022*10**23)
+
+    born = np.loadtxt('born.out')
+    stre = np.loadtxt('vir.out')
+    stre[:, 1:] = -stre[:, 1:]*101325  # -> Pa
+    try:
+        mean_born = np.mean(born[:, 1:], axis=0)
+    except IndexError:
+        mean_born = born[1:]
+
+    CB = kcalmol2J/vol*reduce_Born(mean_born)  # -> J/m^3=Pa
+    Cs = vol/kbT*np.cov(stre[:,1:].T)
+    Ct = N*kbT/vol * compute_delta()
+    C = CB - Cs + Ct
+    write_matrix(CB, 'born_matrix.out')
+    write_matrix(Cs, 'stre_matrix.out')
+    write_matrix(Ct, 'temp_matrix.out')
+    write_matrix(C, 'full_matrix.out')
+    C11 = np.mean([C[0,0], C[1,1], C[2,2]]) * 10**-9
+    C12 = np.mean([C[0,1], C[0,2], C[1,2]]) * 10**-9
+    C44 = np.mean([C[3,3], C[4,4], C[5,5]]) * 10**-9
+    eC11 = np.std([C[0,0], C[1,1], C[2,2]]) * 10**-9
+    eC12 = np.std([C[0,1], C[0,2], C[1,2]]) * 10**-9
+    eC44 = np.std([C[3,3], C[4,4], C[5,5]]) * 10**-9
+    print(C*10**-9)
+    print("C11 = {:f} ± {:f}; C12 = {:f} ± {:f}; C44 = {:f} ± {:f}".format(C11, eC11, C12, eC12, C44, eC44))
+
+    return
+
+
+if __name__ == "__main__":
+    try:
+        main()
+    except KeyboardInterrupt:
+        raise SystemExit("User interruption.")
+

examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/full_matrix.out

@@ -0,0 +1,7 @@
+# Cij Matrix from post process computation
+ 2.29922  1.17439  1.20854 -0.01653 -0.01684  0.01188
+ 1.17439  2.20673  1.21718 -0.00781 -0.00753  0.00867
+ 1.20854  1.21718  2.30804  0.01535 -0.01596  0.00426
+-0.01653 -0.00781  0.01535  1.47647 -0.01355 -0.01601
+-0.01684 -0.00753 -0.01596 -0.01355  1.37905  0.01975
+ 0.01188  0.00867  0.00426 -0.01601  0.01975  1.40170

examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/in.ljcov

@@ -0,0 +1,154 @@
+# Numerical difference calculation
+# of Born matrix
+
+# Note that because of cubic symmetry and central forces, we have:
+# C11, pure axial == positive mean value: 1,2,3
+# C44==C23, pure shear == positive mean value, exactly match in pairs: (4,12),(5,8),(6,7)
+# C14==C56, shear/axial(normal) == zero mean, exactly match in pairs: (9,21),(14,20),(18,19)
+# C15, shear/axial(in-plane) == zero mean: 10,11,13,15,16,17
+
+# adjustable parameters
+
+units           real
+variable        nsteps index 100000    # length of run
+variable        nthermo index  1000    # thermo output interval 
+variable        nlat equal 5           # size of box
+variable        T    equal 60.         # Temperature in K
+variable        rho  equal 5.405       # Lattice spacing in A
+
+atom_style      atomic
+
+lattice         fcc ${rho}
+region          box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
+create_box      1 box
+create_atoms    1 box
+
+mass            * 39.948
+
+velocity        all create ${T} 87287 loop geom
+velocity        all zero linear
+
+pair_style      lj/cut 12.0
+pair_coeff      1 1 0.238067 3.405
+
+neighbor        0.0 bin
+neigh_modify    every 1 delay 0 check no
+
+variable vol equal vol
+thermo 100
+fix aL all ave/time 1 1 1 v_vol ave running
+fix NPT all npt temp $T $T 100 aniso 1. 1. 1000 fixedpoint 0. 0. 0.
+
+run 20000
+
+unfix NPT
+
+variable newL equal "f_aL^(1./3.)"
+change_box all x final 0 ${newL} y final 0. ${newL} z final 0. ${newL} remap units box
+
+unfix aL
+
+reset_timestep 0
+
+# Conversion variables
+variable kb        equal 1.38065e-23    # J/K
+variable Myvol     equal "vol*10^-30" # Volume in m^3
+variable kbt       equal "v_kb*v_T"
+variable Nat       equal atoms
+variable Rhokbt    equal "v_kbt*v_Nat/v_Myvol"
+variable at2Pa     equal 101325
+variable kcalmol2J equal "4183.9954/(6.022e23)"
+variable C1        equal "v_kcalmol2J/v_Myvol" # Convert Cb from energy to pressure units
+variable C2        equal "v_Myvol/v_kbt"       # Factor for Cfl terms
+variable Pa2GPa    equal 1e-9
+
+# Born compute giving <C^b> terms
+# The six virial stress component to compute <C^fl>
+compute     VIR all pressure NULL virial
+compute     born all born/matrix numdiff 1e-6 VIR
+variable s1 equal "-c_VIR[1]*v_at2Pa"
+variable s2 equal "-c_VIR[2]*v_at2Pa"
+variable s3 equal "-c_VIR[3]*v_at2Pa"
+variable s6 equal "-c_VIR[4]*v_at2Pa"
+variable s5 equal "-c_VIR[5]*v_at2Pa"
+variable s4 equal "-c_VIR[6]*v_at2Pa"
+variable press equal press
+
+
+# Average of Born term and vector to store stress
+# for post processing
+fix CB all ave/time 1 ${nthermo} ${nthermo} c_born[*] ave running file born.out overwrite
+fix CPR all ave/time 1 1 1 c_VIR[*] file vir.out
+fix APR all ave/time 1 1 1 v_press ave running
+fix VEC all vector 1 v_s1 v_s2 v_s3 v_s4 v_s5 v_s6
+
+thermo      ${nthermo}
+thermo_style    custom step temp press f_APR c_born[1] f_CB[1] c_born[12] f_CB[12] c_born[4] f_CB[4]
+thermo_modify line multi
+
+fix     1 all nvt temp $T $T 100
+
+run         ${nsteps}
+
+# Compute vector averages
+# Note the indice switch.
+# LAMMPS convention is NOT the Voigt notation.
+variable aves1 equal "ave(f_VEC[1])"
+variable aves2 equal "ave(f_VEC[2])"
+variable aves3 equal "ave(f_VEC[3])"
+variable aves4 equal "ave(f_VEC[6])"
+variable aves5 equal "ave(f_VEC[5])"
+variable aves6 equal "ave(f_VEC[4])"
+
+# Computing the covariance through the <s_{i}s_{j}>-<s_i><s_j>
+# is numerically instable. Here we go through the <(s-<s>)^2>
+# definition.
+
+# Computing difference relative to average values
+variable ds1 vector "f_VEC[1]-v_aves1"
+variable ds2 vector "f_VEC[2]-v_aves2"
+variable ds3 vector "f_VEC[3]-v_aves3"
+variable ds4 vector "f_VEC[4]-v_aves4"
+variable ds5 vector "f_VEC[5]-v_aves5"
+variable ds6 vector "f_VEC[6]-v_aves6"
+
+# Squaring and averaging
+variable dds1 vector "v_ds1*v_ds1"
+variable dds2 vector "v_ds2*v_ds2"
+variable dds3 vector "v_ds3*v_ds3"
+variable vars1 equal "ave(v_dds1)"
+variable vars2 equal "ave(v_dds2)"
+variable vars3 equal "ave(v_dds3)"
+variable C11 equal "v_Pa2GPa*(v_C1*f_CB[1] - v_C2*v_vars1 + 2*v_Rhokbt)"
+variable C22 equal "v_Pa2GPa*(v_C1*f_CB[2] - v_C2*v_vars2 + 2*v_Rhokbt)"
+variable C33 equal "v_Pa2GPa*(v_C1*f_CB[3] - v_C2*v_vars3 + 2*v_Rhokbt)"
+
+variable dds12 vector "v_ds1*v_ds2"
+variable dds13 vector "v_ds1*v_ds3"
+variable dds23 vector "v_ds2*v_ds3"
+variable vars12 equal "ave(v_dds12)"
+variable vars13 equal "ave(v_dds13)"
+variable vars23 equal "ave(v_dds23)"
+variable C12 equal "v_Pa2GPa*(v_C1*f_CB[7] - v_C2*v_vars12)"
+variable C13 equal "v_Pa2GPa*(v_C1*f_CB[8] - v_C2*v_vars13)"
+variable C23 equal "v_Pa2GPa*(v_C1*f_CB[12] - v_C2*v_vars23)"
+
+variable dds4 vector "v_ds4*v_ds4"
+variable dds5 vector "v_ds5*v_ds5"
+variable dds6 vector "v_ds6*v_ds6"
+variable vars4 equal "ave(v_dds4)"
+variable vars5 equal "ave(v_dds5)"
+variable vars6 equal "ave(v_dds6)"
+variable C44 equal "v_Pa2GPa*(v_C1*f_CB[4] - v_C2*v_vars4 + v_Rhokbt)"
+variable C55 equal "v_Pa2GPa*(v_C1*f_CB[5] - v_C2*v_vars5 + v_Rhokbt)"
+variable C66 equal "v_Pa2GPa*(v_C1*f_CB[6] - v_C2*v_vars6 + v_Rhokbt)"
+
+variable aC11 equal "(v_C11 + v_C22 + v_C33)/3."
+variable aC12 equal "(v_C12 + v_C13 + v_C23)/3."
+variable aC44 equal "(v_C44 + v_C55 + v_C66)/3."
+
+print """
+C11 = ${aC11}
+C12 = ${aC12}
+C44 = ${aC44}
+"""

examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/stre_matrix.out

@@ -0,0 +1,7 @@
+# Cij Matrix from post process computation
+ 1.10120  0.69454  0.66285  0.01885  0.01902 -0.01367
+ 0.69454  1.20617  0.65567  0.00893  0.00839 -0.00873
+ 0.66285  0.65567  1.11090 -0.01593  0.01634 -0.00483
+ 0.01885  0.00893 -0.01593  0.42772  0.01316  0.01666
+ 0.01902  0.00839  0.01634  0.01316  0.52416 -0.01733
+-0.01367 -0.00873 -0.00483  0.01666 -0.01733  0.49891

examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/temp_matrix.out

@@ -0,0 +1,7 @@
+# Cij Matrix from post process computation
+ 0.04187  0.00000  0.00000  0.00000  0.00000  0.00000
+ 0.00000  0.04187  0.00000  0.00000  0.00000  0.00000
+ 0.00000  0.00000  0.04187  0.00000  0.00000  0.00000
+ 0.00000  0.00000  0.00000  0.02093  0.00000  0.00000
+ 0.00000  0.00000  0.00000  0.00000  0.02093  0.00000
+ 0.00000  0.00000  0.00000  0.00000  0.00000  0.02093

examples/ELASTIC_T/BORN_MATRIX/Argon/README.md

@@ -0,0 +1,13 @@
+This repository is a test case for the compute born/matrix. It provides short
+scripts creating argon fcc crystal and computing the Born term using the two
+methods described in the documentation.
+
+In the __Analytical__ directory the terms are computed using the analytical
+derivation of the Born term for the lj/cut pair style.
+
+In the __Numdiff__ directory, the Born term is evaluated through small
+numerical differences of the stress tensor. This method can be used with any
+interaction potential.
+
+Both script show examples on how to compute the full Cij elastic stiffness
+tensor in LAMMPS.

examples/ELASTIC_T/BORN_MATRIX/Silicon/Si.sw

@@ -0,0 +1 @@
+../../../../potentials/Si.sw

examples/ELASTIC_T/BORN_MATRIX/Silicon/final_output.in

@@ -0,0 +1,68 @@
+# Average moduli for cubic crystals
+
+variable C11cubic equal (${C11}+${C22}+${C33})/3.0
+variable C12cubic equal (${C12}+${C13}+${C23})/3.0
+variable C44cubic equal (${C44}+${C55}+${C66})/3.0
+
+variable bulkmodulus equal (${C11cubic}+2*${C12cubic})/3.0
+variable shearmodulus1 equal ${C44cubic}
+variable shearmodulus2 equal (${C11cubic}-${C12cubic})/2.0
+variable poissonratio equal 1.0/(1.0+${C11cubic}/${C12cubic})
+  
+# For Stillinger-Weber silicon, the analytical results
+# are known to be (E. R. Cowley, 1988):
+#               C11 = 151.4 GPa
+#               C12 = 76.4 GPa
+#               C44 = 56.4 GPa
+
+#print "========================================="
+#print "Components of the Elastic Constant Tensor"
+#print "========================================="
+
+print "Elastic Constant C11 = ${C11} ${cunits}"
+print "Elastic Constant C22 = ${C22} ${cunits}"
+print "Elastic Constant C33 = ${C33} ${cunits}"
+
+print "Elastic Constant C12 = ${C12} ${cunits}"
+print "Elastic Constant C13 = ${C13} ${cunits}"
+print "Elastic Constant C23 = ${C23} ${cunits}"
+
+print "Elastic Constant C44 = ${C44} ${cunits}"
+print "Elastic Constant C55 = ${C55} ${cunits}"
+print "Elastic Constant C66 = ${C66} ${cunits}"
+
+print "Elastic Constant C14 = ${C14} ${cunits}"
+print "Elastic Constant C15 = ${C15} ${cunits}"
+print "Elastic Constant C16 = ${C16} ${cunits}"
+
+print "Elastic Constant C24 = ${C24} ${cunits}"
+print "Elastic Constant C25 = ${C25} ${cunits}"
+print "Elastic Constant C26 = ${C26} ${cunits}"
+
+print "Elastic Constant C34 = ${C34} ${cunits}"
+print "Elastic Constant C35 = ${C35} ${cunits}"
+print "Elastic Constant C36 = ${C36} ${cunits}"
+
+print "Elastic Constant C45 = ${C45} ${cunits}"
+print "Elastic Constant C46 = ${C46} ${cunits}"
+print "Elastic Constant C56 = ${C56} ${cunits}"
+
+print "========================================="
+print "Average properties for a cubic crystal"
+print "========================================="
+
+print "Bulk Modulus = ${bulkmodulus} ${cunits}"
+print "Shear Modulus 1 = ${shearmodulus1} ${cunits}"
+print "Shear Modulus 2 = ${shearmodulus2} ${cunits}"
+print "Poisson Ratio = ${poissonratio}"
+
+# summarize sampling protocol
+
+variable tmp equal atoms
+print "Number of atoms = ${tmp}"
+print "Stress sampling interval = ${nevery}"
+variable tmp equal ${nrun}/${nevery}
+print "Stress sample count = ${tmp}"
+print "Born sampling interval = ${neveryborn}"
+variable tmp equal ${nrun}/${neveryborn}
+print "Born sample count = ${tmp}"

examples/ELASTIC_T/BORN_MATRIX/Silicon/in.elastic

@@ -0,0 +1,25 @@
+# Compute elastic constant tensor for a crystal at finite temperature
+# These settings replicate the 1477~K benchmark of
+# Kluge, Ray, and Rahman (1986) that is Ref.[15] in:
+# Y. Zhen, C. Chu, Computer Physics Communications 183(2012) 261-265
+
+# here: Y. Zhen, C. Chu, Computer Physics Communications 183(2012) 261-265
+
+include init.in
+
+# Compute initial state
+
+include potential.in
+thermo_style custom step temp pe press density
+run ${nequil}
+
+# Run dynamics
+
+include potential.in
+include output.in
+
+run ${nrun}
+
+# Output final values
+
+include final_output.in

examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in

@@ -0,0 +1,59 @@
+# NOTE: This script can be modified for different atomic structures, 
+# units, etc. See in.elastic for more info.
+#
+
+# Define MD parameters
+# These can be modified by the user
+# These settings replicate the 1477~K benchmark of
+# Kluge, Ray, and Rahman (1986) that is Ref.[15] in:
+# Y. Zhen, C. Chu, Computer Physics Communications 183(2012) 261-265
+
+# select temperature and pressure (lattice constant)
+
+variable temp index 1477.0                 # temperature of initial sample
+variable a index 5.457                     # lattice constant
+
+# select sampling parameters, important for speed/convergence
+
+variable nthermo index 1500               # interval for thermo output
+variable nevery index 10                  # stress sampling interval
+variable neveryborn index 100             # Born sampling interval
+variable timestep index 0.000766          # timestep
+variable nlat index 3                     # number of lattice unit cells
+
+# other settings
+
+variable mass1 index 28.06                # mass
+variable tdamp index 0.01                 # time constant for thermostat
+variable seed index 123457                # seed for thermostat
+variable thermostat index 1               # 0 if NVE, 1 if NVT
+variable delta index 1.0e-6               # Born numdiff strain magnitude 
+
+# hard-coded rules-of-thumb for run length, etc.
+
+variable nfreq equal ${nthermo}           # interval for averaging output
+variable nrepeat equal floor(${nfreq}/${nevery}) # number of samples
+variable nrepeatborn equal floor(${nfreq}/${neveryborn}) # number of samples
+variable nequil equal 10*${nthermo}       # length of equilibration run
+variable nrun equal 100*${nthermo}        # length of equilibrated run
+
+# generate the box and atom positions using a diamond lattice
+
+units		metal
+
+boundary	p p p
+
+# this generates a standard 8-atom cubic cell
+
+lattice         diamond $a
+region		box prism 0 1 0 1 0 1 0 0 0
+
+# this generates a 2-atom triclinic cell
+#include tri.in
+
+create_box	1 box
+create_atoms	1 box
+mass 1 ${mass1}
+replicate ${nlat} ${nlat} ${nlat}
+velocity	all create ${temp} 87287
+

examples/ELASTIC_T/BORN_MATRIX/Silicon/output.in

@@ -0,0 +1,140 @@
+# Setup output
+
+# Stress fluctuation term F
+
+compute stress all pressure thermo_temp
+variable s1 equal c_stress[1]
+variable s2 equal c_stress[2]
+variable s3 equal c_stress[3]
+variable s4 equal c_stress[6]
+variable s5 equal c_stress[5]
+variable s6 equal c_stress[4]
+
+variable s11 equal v_s1*v_s1
+variable s22 equal v_s2*v_s2
+variable s33 equal v_s3*v_s3
+variable s44 equal v_s4*v_s4
+variable s55 equal v_s5*v_s5
+variable s66 equal v_s6*v_s6
+variable s33 equal v_s3*v_s3
+variable s12 equal v_s1*v_s2
+variable s13 equal v_s1*v_s3
+variable s14 equal v_s1*v_s4
+variable s15 equal v_s1*v_s5
+variable s16 equal v_s1*v_s6
+variable s23 equal v_s2*v_s3
+variable s24 equal v_s2*v_s4
+variable s25 equal v_s2*v_s5
+variable s26 equal v_s2*v_s6
+variable s34 equal v_s3*v_s4
+variable s35 equal v_s3*v_s5
+variable s36 equal v_s3*v_s6
+variable s45 equal v_s4*v_s5
+variable s46 equal v_s4*v_s6
+variable s56 equal v_s5*v_s6
+
+variable mytemp equal temp
+variable mypress equal press
+variable mype equal pe/atoms
+fix avt all ave/time ${nevery} ${nrepeat} ${nfreq} v_mytemp ave running
+fix avp all ave/time ${nevery} ${nrepeat} ${nfreq} v_mypress ave running
+fix avpe all ave/time ${nevery} ${nrepeat} ${nfreq} v_mype ave running
+fix avs all ave/time ${nevery} ${nrepeat} ${nfreq} v_s1 v_s2 v_s3 v_s4 v_s5 v_s6 ave running
+fix avssq all ave/time ${nevery} ${nrepeat} ${nfreq} &
+v_s11 v_s22 v_s33 v_s44 v_s55 v_s66 &
+v_s12 v_s13 v_s14 v_s15 v_s16 &
+v_s23 v_s24 v_s25 v_s26 &
+v_s34 v_s35 v_s36 &
+v_s45 v_s46 &
+v_s56 &
+ave running
+
+# bar to GPa
+variable pconv equal 1.0e5/1.0e9 
+variable cunits index GPa
+# metal unit constants from LAMMPS
+# force->nktv2p = 1.6021765e6;
+# force->boltz = 8.617343e-5;
+variable boltz equal 8.617343e-5
+variable nktv2p equal 1.6021765e6
+variable vkt equal vol/(${boltz}*${temp})/${nktv2p}
+variable ffac equal ${pconv}*${vkt}
+
+variable F11 equal -(f_avssq[1]-f_avs[1]*f_avs[1])*${ffac}
+variable F22 equal -(f_avssq[2]-f_avs[2]*f_avs[2])*${ffac}
+variable F33 equal -(f_avssq[3]-f_avs[3]*f_avs[3])*${ffac}
+variable F44 equal -(f_avssq[4]-f_avs[4]*f_avs[4])*${ffac}
+variable F55 equal -(f_avssq[5]-f_avs[5]*f_avs[5])*${ffac}
+variable F66 equal -(f_avssq[6]-f_avs[6]*f_avs[6])*${ffac}
+
+variable F12 equal -(f_avssq[7]-f_avs[1]*f_avs[2])*${ffac}
+variable F13 equal -(f_avssq[8]-f_avs[1]*f_avs[3])*${ffac}
+variable F14 equal -(f_avssq[9]-f_avs[1]*f_avs[4])*${ffac}
+variable F15 equal -(f_avssq[10]-f_avs[1]*f_avs[5])*${ffac}
+variable F16 equal -(f_avssq[11]-f_avs[1]*f_avs[6])*${ffac}
+
+variable F23 equal -(f_avssq[12]-f_avs[2]*f_avs[3])*${ffac}
+variable F24 equal -(f_avssq[13]-f_avs[2]*f_avs[4])*${ffac}
+variable F25 equal -(f_avssq[14]-f_avs[2]*f_avs[5])*${ffac}
+variable F26 equal -(f_avssq[15]-f_avs[2]*f_avs[6])*${ffac}
+
+variable F34 equal -(f_avssq[16]-f_avs[3]*f_avs[4])*${ffac}
+variable F35 equal -(f_avssq[17]-f_avs[3]*f_avs[5])*${ffac}
+variable F36 equal -(f_avssq[18]-f_avs[3]*f_avs[6])*${ffac}
+
+variable F45 equal -(f_avssq[19]-f_avs[4]*f_avs[5])*${ffac}
+variable F46 equal -(f_avssq[20]-f_avs[4]*f_avs[6])*${ffac}
+
+variable F56 equal -(f_avssq[21]-f_avs[5]*f_avs[6])*${ffac}
+
+# Born term
+
+compute     virial all pressure NULL virial
+compute     born all born/matrix numdiff ${delta} virial
+fix 	    avborn all ave/time ${neveryborn} ${nrepeatborn} ${nfreq} c_born[*] ave running
+
+variable bfac equal ${pconv}*${nktv2p}/vol
+variable B vector f_avborn*${bfac}
+
+# Kinetic term
+
+variable kfac equal ${pconv}*${nktv2p}*atoms*${boltz}*${temp}/vol
+variable K11 equal 4.0*${kfac}
+variable K22 equal 4.0*${kfac}
+variable K33 equal 4.0*${kfac}
+variable K44 equal 2.0*${kfac}
+variable K55 equal 2.0*${kfac}
+variable K66 equal 2.0*${kfac}
+
+# Add F, K, and B together
+
+variable C11 equal v_F11+v_B[1]+v_K11
+variable C22 equal v_F22+v_B[2]+v_K22
+variable C33 equal v_F33+v_B[3]+v_K33
+variable C44 equal v_F44+v_B[4]+v_K44
+variable C55 equal v_F55+v_B[5]+v_K55
+variable C66 equal v_F66+v_B[6]+v_K66
+
+variable C12 equal v_F12+v_B[7] 
+variable C13 equal v_F13+v_B[8] 
+variable C14 equal v_F14+v_B[9] 
+variable C15 equal v_F15+v_B[10]
+variable C16 equal v_F16+v_B[11]
+
+variable C23 equal v_F23+v_B[12] 
+variable C24 equal v_F24+v_B[13] 
+variable C25 equal v_F25+v_B[14]
+variable C26 equal v_F26+v_B[15]
+
+variable C34 equal v_F34+v_B[16]
+variable C35 equal v_F35+v_B[17]
+variable C36 equal v_F36+v_B[18]
+
+variable C45 equal v_F45+v_B[19]
+variable C46 equal v_F46+v_B[20]
+
+variable C56 equal v_F56+v_B[21]
+
+thermo		${nthermo}
+thermo_style custom step temp pe press density f_avt f_avp f_avpe v_F11 v_F22 v_F33 v_F44 v_F55 v_F66 v_F12 v_F13 v_F23 v_B[1] v_B[2] v_B[3] v_B[4] v_B[5] v_B[6] v_B[7] v_B[8] v_B[12]
+thermo_modify norm no

examples/ELASTIC_T/BORN_MATRIX/Silicon/potential.in

@@ -0,0 +1,21 @@
+# NOTE: This script can be modified for different pair styles 
+# See in.elastic for more info.
+
+reset_timestep 0
+
+# Choose potential
+pair_style	sw
+pair_coeff * * Si.sw Si
+
+# Setup neighbor style
+neighbor 1.0 nsq
+neigh_modify once no every 1 delay 0 check yes
+
+# Setup MD
+
+timestep ${timestep}
+fix 4 all nve
+if "${thermostat} == 1" then &
+   "fix 5 all langevin ${temp} ${temp} ${tdamp} ${seed}"
+
+

examples/ELASTIC_T/BORN_MATRIX/Silicon/tri.in

@@ -0,0 +1,26 @@
+# this generates a 2-atom triclinic cell
+# due to rotation on to x-axis,
+# elastic constant analysis is not working yet
+
+# unit lattice vectors are
+# a1 = (1 0 0)
+# a2 = (1/2 sqrt3/2 0)
+# a3 = (1/2 1/(2sqrt3) sqrt2/sqrt3)
+
+variable        a1x equal 1
+variable 	a2x equal 1/2
+variable 	a2y equal sqrt(3)/2
+variable 	a3x equal 1/2
+variable 	a3y equal 1/(2*sqrt(3))
+variable	a3z equal sqrt(2/3)
+variable 	l equal $a/sqrt(2)
+
+lattice         custom ${l}             &
+                a1 ${a1x} 0 0           &
+                a2 ${a2x} ${a2y} 0.0    &
+                a3 ${a3x} ${a3y} ${a3z} &
+                basis 0 0 0             &
+                basis 0.25 0.25 0.25    &
+		spacing 1 1 1
+
+region		box prism 0 ${a1x} 0 ${a2y} 0 ${a3z} ${a2x} ${a3x} ${a3y}

examples/ELASTIC_T/DEFORMATION/Silicon/Si.sw

@@ -0,0 +1 @@
+../../../../potentials/Si.sw

examples/ELASTIC_T/DEFORMATION/Silicon/potential.mod

@@ -2,7 +2,7 @@
 # See in.elastic for more info.
 
 # we must undefine any fix ave/* fix before using reset_timestep
-if "$(is_defined(fix,avp)" then "unfix avp"
+if "$(is_defined(fix,avp))" then "unfix avp"
 reset_timestep 0
 
 # Choose potential

examples/ELASTIC_T/Si.sw

@@ -1,18 +0,0 @@
-# DATE: 2007-06-11 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov CITATION: Stillinger and Weber,  Phys Rev B, 31, 5262, (1985)
-# Stillinger-Weber parameters for various elements and mixtures
-# multiple entries can be added to this file, LAMMPS reads the ones it needs
-# these entries are in LAMMPS "metal" units:
-#   epsilon = eV; sigma = Angstroms
-#   other quantities are unitless
-
-# format of a single entry (one or more lines):
-#   element 1, element 2, element 3, 
-#   epsilon, sigma, a, lambda, gamma, costheta0, A, B, p, q, tol
-
-# Here are the original parameters in metal units, for Silicon from:
-#
-# Stillinger and Weber,  Phys. Rev. B, v. 31, p. 5262, (1985)
-#
-
-Si Si Si 2.1683  2.0951  1.80  21.0  1.20  -0.333333333333
-         7.049556277  0.6022245584  4.0  0.0 0.0

examples/PACKAGES/dielectric/data.confined

@@ -1,4 +1,4 @@
-LAMMPS data file: two oppositely charged ions confined between two walls
+LAMMPS data file: two oppositely charged ions confined between two walls epsilon1=2 | epsilon2=10 | epsilon1=2
 
 4002 atoms
 3 atom types
@@ -4015,5 +4015,5 @@ Atoms # dielectric: id mol type q x y z normx normy normz area_per_patch ed em e
 3998 0 1 0 38.5078 42.4438 30.002  0 0 -1 0.866 8 6 6 0
 3999 0 1 0 39.0079 41.5776 30.002  0 0 -1 0.866 8 6 6 0
 4000 0 1 0 39.508 42.4438 30.002  0 0 -1 0.866 8 6 6 0
-4001 0 2 1 15 20 15 0 0 1 0.866 8 6 10 0
+4001 0 2 1  15 20 15 0 0 1 1.0 8 6 10 0
-4002 0 3 -1 25 20 25 0 0 1 0.866 8 6 10 0
+4002 0 3 -1 25 20 25 0 0 1 1.0 8 6 10 0

examples/PACKAGES/dielectric/in.confined

@@ -4,6 +4,8 @@
 # top interface:    n = (0, 0, -1)
 # so that ed's are the same for both interfaces
 
+# Dielectric constants can be set to be different from the input data file
+
 variable epsilon1 index 20
 variable epsilon2 index 8
 
@@ -21,6 +23,8 @@ variable    method index gmres  # gmres = BEM/GMRES
                                 # dof   = Direct optimization of the functional
                                 # none
 
+# compute the relevant values for the interface particles
+
 variable ed      equal "v_epsilon2 - v_epsilon1"
 variable em      equal "(v_epsilon2 + v_epsilon1)/2"
 variable epsilon equal 1.0                             # epsilon at the patch, not used for now
@@ -34,12 +38,10 @@ group ions type 2 3
 group cations type 2
 group anions  type 3
 
-# 1.0 = q * epsilon2 = qreal for cations
+# set the dielectric constant of the medium where the ions reside
-# -1.0 = q * epsilon2 = qreal for anions
+
-variable qscale equal "1.0 / v_epsilon2"
+set group cations epsilon ${epsilon2}
-set group cations charge ${qscale}
+set group anions  epsilon ${epsilon2}
-variable qscale equal "-1.0 / v_epsilon2"
-set group anions charge ${qscale}
 
 pair_style      lj/cut/coul/long/dielectric 1.122 10.0
 pair_coeff      * * 1.0 1.0
@@ -59,6 +61,8 @@ dump            3 ions custom 100 ions.dump id mol type q x y z fx fy fz #c_ef[1
 
 fix             1 ions nve
 
+# fix modify is used to set the properties of the interface particle group
+
 if "${method} == gmres" then &
   "fix             3 interface polarize/bem/gmres 1 1.0e-4" &
   "fix_modify      3 itr_max 50 dielectrics ${ed} ${em} ${epsilon} ${area} NULL" &