Merge pull request #2899 from akohlmey/next_patch_release

Step version strings for the next patch release

.github/CODEOWNERS

@@ -129,7 +129,7 @@ src/math_eigen_impl.h     @jewettaij
 tools/msi2lmp/*       @akohlmey
 tools/emacs/*         @HaoZeke
 tools/singularity/*   @akohlmey @rbberger
-tools/code_standard/* @rbberger
+tools/coding_standard/* @rbberger
 tools/valgrind/*      @akohlmey
 tools/swig/*          @akohlmey
 tools/offline/*       @rbberger

cmake/.coveragerc.in

@@ -0,0 +1,10 @@
+[run]
+source = @LAMMPS_PYTHON_DIR@
+parallel=True
+branch=True
+omit=*/install.py
+     */setup.py
+
+[paths]
+sources = python
+          @LAMMPS_PYTHON_DIR@

cmake/CMakeLists.txt

@@ -769,6 +769,13 @@ endif()
 include(Testing)
 include(CodeCoverage)
 include(CodingStandard)
+find_package(ClangFormat 8.0)
+
+if(ClangFormat_FOUND)
+  add_custom_target(format-src
+    COMMAND ${ClangFormat_EXECUTABLE} --verbose -i -style=file *.cpp *.h */*.cpp */*.h
+    WORKING_DIRECTORY ${LAMMPS_SOURCE_DIR})
+endif()
 
 get_target_property(DEFINES lammps COMPILE_DEFINITIONS)
 include(FeatureSummary)

cmake/Modules/CodeCoverage.cmake

@@ -54,6 +54,8 @@ if(ENABLE_COVERAGE)
 
     if(COVERAGE_FOUND)
         set(PYTHON_COVERAGE_HTML_DIR ${CMAKE_BINARY_DIR}/python_coverage_html)
+        configure_file(.coveragerc.in ${CMAKE_BINARY_DIR}/.coveragerc @ONLY)
+
         add_custom_command(
             OUTPUT ${CMAKE_BINARY_DIR}/unittest/python/.coverage
             COMMAND ${COVERAGE_BINARY} combine
@@ -63,16 +65,16 @@ if(ENABLE_COVERAGE)
 
         add_custom_target(
             gen_python_coverage_html
-            COMMAND ${COVERAGE_BINARY} html -d ${PYTHON_COVERAGE_HTML_DIR}
-            DEPENDS ${CMAKE_BINARY_DIR}/unittest/python/.coverage
+            COMMAND ${COVERAGE_BINARY} html --rcfile=${CMAKE_BINARY_DIR}/.coveragerc -d ${PYTHON_COVERAGE_HTML_DIR}
+            DEPENDS ${CMAKE_BINARY_DIR}/unittest/python/.coverage ${CMAKE_BINARY_DIR}/.coveragerc
             WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/unittest/python
             COMMENT "Generating HTML Python coverage report..."
         )
 
         add_custom_target(
             gen_python_coverage_xml
-            COMMAND ${COVERAGE_BINARY} xml -o ${CMAKE_BINARY_DIR}/python_coverage.xml
-            DEPENDS ${CMAKE_BINARY_DIR}/unittest/python/.coverage
+            COMMAND ${COVERAGE_BINARY} xml --rcfile=${CMAKE_BINARY_DIR}/.coveragerc -o ${CMAKE_BINARY_DIR}/python_coverage.xml
+            DEPENDS ${CMAKE_BINARY_DIR}/unittest/python/.coverage ${CMAKE_BINARY_DIR}/.coveragerc
             WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/unittest/python
             COMMENT "Generating XML Python coverage report..."
         )

cmake/Modules/FindFFTW3.cmake

@@ -38,7 +38,7 @@ if(FFTW3_FOUND)
       add_library(FFTW3::FFTW3_OMP UNKNOWN IMPORTED)
       set_target_properties(FFTW3::FFTW3_OMP PROPERTIES
         IMPORTED_LINK_INTERFACE_LANGUAGES "C"
-	IMPORTED_LOCATION "${FFTW3_OMP_LIBRARY}"
+        IMPORTED_LOCATION "${FFTW3_OMP_LIBRARY}"
         INTERFACE_INCLUDE_DIRECTORIES "${FFTW3_INCLUDE_DIRS}")
     endif()
   endif()

cmake/Modules/FindFFTW3F.cmake

@@ -37,7 +37,7 @@ if(FFTW3F_FOUND)
       add_library(FFTW3F::FFTW3F_OMP UNKNOWN IMPORTED)
       set_target_properties(FFTW3F::FFTW3F_OMP PROPERTIES
         IMPORTED_LINK_INTERFACE_LANGUAGES "C"
-	IMPORTED_LOCATION "${FFTW3F_OMP_LIBRARY}"
+        IMPORTED_LOCATION "${FFTW3F_OMP_LIBRARY}"
         INTERFACE_INCLUDE_DIRECTORIES "${FFTW3F_INCLUDE_DIRS}")
     endif()
   endif()

cmake/Modules/Packages/INTEL.cmake

@@ -3,7 +3,7 @@ if(NOT FOUND_IMMINTRIN)
   message(FATAL_ERROR "immintrin.h header not found, Intel package won't work without it")
 endif()
 
-target_compile_definitions(lammps PRIVATE -DLMP_USER_INTEL)
+target_compile_definitions(lammps PRIVATE -DLMP_INTEL)
 
 set(INTEL_ARCH "cpu" CACHE STRING "Architectures used by INTEL (cpu or knl)")
 set(INTEL_ARCH_VALUES cpu knl)

cmake/Modules/Packages/KOKKOS.cmake

@@ -74,7 +74,7 @@ else()
   target_link_libraries(lammps PRIVATE kokkos)
   target_link_libraries(lmp PRIVATE kokkos)
 endif()
-target_compile_definitions(lammps PRIVATE -DLMP_KOKKOS)
+target_compile_definitions(lammps PUBLIC $<BUILD_INTERFACE:LMP_KOKKOS>)
 
 set(KOKKOS_PKG_SOURCES_DIR ${LAMMPS_SOURCE_DIR}/KOKKOS)
 set(KOKKOS_PKG_SOURCES ${KOKKOS_PKG_SOURCES_DIR}/kokkos.cpp
@@ -127,4 +127,4 @@ endif()
 get_property(KOKKOS_PKG_SOURCES GLOBAL PROPERTY KOKKOS_PKG_SOURCES)
 
 target_sources(lammps PRIVATE ${KOKKOS_PKG_SOURCES})
-target_include_directories(lammps PRIVATE ${KOKKOS_PKG_SOURCES_DIR})
+target_include_directories(lammps PUBLIC $<BUILD_INTERFACE:${KOKKOS_PKG_SOURCES_DIR}>)

cmake/Modules/Packages/MACHDYN.cmake

@@ -7,8 +7,8 @@ endif()
 option(DOWNLOAD_EIGEN3 "Download Eigen3 instead of using an already installed one)" ${DOWNLOAD_EIGEN3_DEFAULT})
 if(DOWNLOAD_EIGEN3)
   message(STATUS "Eigen3 download requested - we will build our own")
-  set(EIGEN3_URL "https://gitlab.com/libeigen/eigen/-/archive/3.3.7/eigen-3.3.7.tar.gz" CACHE STRING "URL for Eigen3 tarball")
-  set(EIGEN3_MD5 "9e30f67e8531477de4117506fe44669b" CACHE STRING "MD5 checksum of Eigen3 tarball")
+  set(EIGEN3_URL "https://gitlab.com/libeigen/eigen/-/archive/3.3.9/eigen-3.3.9.tar.gz" CACHE STRING "URL for Eigen3 tarball")
+  set(EIGEN3_MD5 "609286804b0f79be622ccf7f9ff2b660" CACHE STRING "MD5 checksum of Eigen3 tarball")
   mark_as_advanced(EIGEN3_URL)
   mark_as_advanced(EIGEN3_MD5)
   include(ExternalProject)
@@ -30,3 +30,8 @@ else()
   endif()
   target_link_libraries(lammps PRIVATE Eigen3::Eigen)
 endif()
+
+# PGI/Nvidia compiler internals collide with vector intrinsics support in Eigen3
+if((CMAKE_CXX_COMPILER_ID STREQUAL "PGI") OR (CMAKE_CXX_COMPILER_ID STREQUAL "NVHPC"))
+  target_compile_definitions(lammps PRIVATE -DEIGEN_DONT_VECTORIZE)
+endif()

cmake/Modules/Packages/MDI.cmake

@@ -114,5 +114,5 @@ else()
   target_link_libraries(lmp PRIVATE ${mdi_LIBRARY})
 endif()
 
-target_compile_definitions(lammps PRIVATE -DLMP_USER_MDI)
-target_compile_definitions(lmp PRIVATE -DLMP_USER_MDI)
+target_compile_definitions(lammps PRIVATE -DLMP_MDI)
+target_compile_definitions(lmp PRIVATE -DLMP_MDI)

cmake/Modules/Packages/OPENMP.cmake

@@ -5,7 +5,7 @@
                        ${OPENMP_SOURCES_DIR}/fix_nh_omp.cpp
                        ${OPENMP_SOURCES_DIR}/fix_nh_sphere_omp.cpp
                        ${OPENMP_SOURCES_DIR}/domain_omp.cpp)
-  target_compile_definitions(lammps PRIVATE -DLMP_USER_OMP)
+  target_compile_definitions(lammps PRIVATE -DLMP_OPENMP)
   set_property(GLOBAL PROPERTY "OMP_SOURCES" "${OPENMP_SOURCES}")
 
   # detects styles which have OPENMP version

cmake/Modules/Tools.cmake

@@ -9,14 +9,16 @@ if(BUILD_TOOLS)
     check_language(Fortran)
     if(CMAKE_Fortran_COMPILER)
       enable_language(Fortran)
-      add_executable(chain.x ${LAMMPS_TOOLS_DIR}/chain.f)
+      add_executable(chain.x ${LAMMPS_TOOLS_DIR}/chain.f90)
       target_link_libraries(chain.x PRIVATE ${CMAKE_Fortran_IMPLICIT_LINK_LIBRARIES})
-      install(TARGETS chain.x DESTINATION ${CMAKE_INSTALL_BINDIR})
+      add_executable(micelle2d.x ${LAMMPS_TOOLS_DIR}/micelle2d.f90)
+      target_link_libraries(micelle2d.x PRIVATE ${CMAKE_Fortran_IMPLICIT_LINK_LIBRARIES})
+      install(TARGETS chain.x micelle2d.x DESTINATION ${CMAKE_INSTALL_BINDIR})
     else()
-      message(WARNING "No suitable Fortran compiler found, skipping build of 'chain.x'")
+      message(WARNING "No suitable Fortran compiler found, skipping build of 'chain.x' and 'micelle2d.x'")
     endif()
   else()
-    message(WARNING "CMake build doesn't support fortran, skipping build of 'chain.x'")
+    message(WARNING "CMake build doesn't support Fortran, skipping build of 'chain.x' and 'micelle2d.x'")
   endif()
 
   enable_language(C)

cmake/presets/clang.cmake

@@ -10,9 +10,9 @@ set(CMAKE_Fortran_COMPILER ${CLANG_FORTRAN} CACHE STRING "" FORCE)
 set(CMAKE_CXX_FLAGS_DEBUG "-Wall -Wextra -g" CACHE STRING "" FORCE)
 set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "-Wall -Wextra -g -O2 -DNDEBUG" CACHE STRING "" FORCE)
 set(CMAKE_CXX_FLAGS_RELEASE "-O3 -DNDEBUG" CACHE STRING "" FORCE)
-set(CMAKE_Fortran_FLAGS_DEBUG "-Wall -Wextra -g -std=f95" CACHE STRING "" FORCE)
-set(CMAKE_Fortran_FLAGS_RELWITHDEBINFO "-Wall -Wextra -g -O2 -DNDEBUG -std=f95" CACHE STRING "" FORCE)
-set(CMAKE_Fortran_FLAGS_RELEASE "-O3 -DNDEBUG -std=f95" CACHE STRING "" FORCE)
+set(CMAKE_Fortran_FLAGS_DEBUG "-Wall -Wextra -g -std=f2003" CACHE STRING "" FORCE)
+set(CMAKE_Fortran_FLAGS_RELWITHDEBINFO "-Wall -Wextra -g -O2 -DNDEBUG -std=f2003" CACHE STRING "" FORCE)
+set(CMAKE_Fortran_FLAGS_RELEASE "-O3 -DNDEBUG -std=f2003" CACHE STRING "" FORCE)
 set(CMAKE_C_FLAGS_DEBUG "-Wall -Wextra -g" CACHE STRING "" FORCE)
 set(CMAKE_C_FLAGS_RELWITHDEBINFO "-Wall -Wextra -g -O2 -DNDEBUG" CACHE STRING "" FORCE)
 set(CMAKE_C_FLAGS_RELEASE "-O3 -DNDEBUG" CACHE STRING "" FORCE)

cmake/presets/gcc.cmake

@@ -1,4 +1,4 @@
-# preset that will restore gcc/g++ with support for MPI and OpenMP (on Linux boxes)
+# preset that will explicitly request gcc/g++ compilers with support for MPI and OpenMP
 
 set(CMAKE_CXX_COMPILER "g++" CACHE STRING "" FORCE)
 set(CMAKE_C_COMPILER "gcc" CACHE STRING "" FORCE)
@@ -15,9 +15,9 @@ set(CMAKE_C_FLAGS_RELWITHDEBINFO "-g -O2 -DNDEBUG" CACHE STRING "" FORCE)
 set(CMAKE_C_FLAGS_RELEASE "-O3 -DNDEBUG" CACHE STRING "" FORCE)
 set(MPI_Fortran "gfortran" CACHE STRING "" FORCE)
 set(MPI_Fortran_COMPILER "mpifort" CACHE STRING "" FORCE)
-set(CMAKE_Fortran_FLAGS_DEBUG "-Wall -g" CACHE STRING "" FORCE)
-set(CMAKE_Fortran_FLAGS_RELWITHDEBINFO "-g -O2 -DNDEBUG" CACHE STRING "" FORCE)
-set(CMAKE_Fortran_FLAGS_RELEASE "-O3 -DNDEBUG" CACHE STRING "" FORCE)
+set(CMAKE_Fortran_FLAGS_DEBUG "-Wall -g -std=f2003" CACHE STRING "" FORCE)
+set(CMAKE_Fortran_FLAGS_RELWITHDEBINFO "-g -O2 -DNDEBUG -std=f2003" CACHE STRING "" FORCE)
+set(CMAKE_Fortran_FLAGS_RELEASE "-O3 -DNDEBUG -std=f2003" CACHE STRING "" FORCE)
 unset(HAVE_OMP_H_INCLUDE CACHE)
 
 set(OpenMP_C "gcc" CACHE STRING "" FORCE)

cmake/presets/hip.cmake

@@ -1,12 +1,26 @@
-# preset that will enable hipcc plus gcc with support for MPI and OpenMP (on Linux boxes)
+# preset that will enable hipcc plus gcc/gfortran with support for MPI and OpenMP (on Linux boxes)
 
 set(CMAKE_CXX_COMPILER "hipcc" CACHE STRING "" FORCE)
 set(CMAKE_C_COMPILER "gcc" CACHE STRING "" FORCE)
+set(CMAKE_Fortran_COMPILER gfortran CACHE STRING "" FORCE)
 set(CMAKE_CXX_FLAGS_DEBUG "-Wall -Wextra -g" CACHE STRING "" FORCE)
 set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "-Wall -Wextra -g -O2 -DNDEBUG" CACHE STRING "" FORCE)
-unset(HAVE_OMP_H_INCLUDE CACHE)
+set(CMAKE_CXX_FLAGS_RELEASE "-O3 -DNDEBUG" CACHE STRING "" FORCE)
+set(CMAKE_Fortran_FLAGS_DEBUG "-Wall -Wextra -g -std=f2003" CACHE STRING "" FORCE)
+set(CMAKE_Fortran_FLAGS_RELWITHDEBINFO "-Wall -Wextra -g -O2 -DNDEBUG -std=f2003" CACHE STRING "" FORCE)
+set(CMAKE_Fortran_FLAGS_RELEASE "-O3 -DNDEBUG -std=f2003" CACHE STRING "" FORCE)
+set(CMAKE_C_FLAGS_DEBUG "-Wall -Wextra -g" CACHE STRING "" FORCE)
+set(CMAKE_C_FLAGS_RELWITHDEBINFO "-Wall -Wextra -g -O2 -DNDEBUG" CACHE STRING "" FORCE)
+set(CMAKE_C_FLAGS_RELEASE "-O3 -DNDEBUG" CACHE STRING "" FORCE)
 
-set(OpenMP_CXX "hipcc" CACHE STRING "" FORCE)
+set(MPI_CXX "hipcc" CACHE STRING "" FORCE)
+set(MPI_CXX_COMPILER "mpicxx" CACHE STRING "" FORCE)
+
+unset(HAVE_OMP_H_INCLUDE CACHE)
+set(OpenMP_C "gcc" CACHE STRING "" FORCE)
+set(OpenMP_C_FLAGS "-fopenmp" CACHE STRING "" FORCE)
+set(OpenMP_C_LIB_NAMES "gomp" CACHE STRING "" FORCE)
 set(OpenMP_CXX_FLAGS "-fopenmp" CACHE STRING "" FORCE)
+set(OpenMP_CXX "hipcc" CACHE STRING "" FORCE)
 set(OpenMP_CXX_LIB_NAMES "omp" CACHE STRING "" FORCE)
 set(OpenMP_omp_LIBRARY "libomp.so" CACHE PATH "" FORCE)

doc/lammps.1

@@ -1,4 +1,4 @@
-.TH LAMMPS "30 July 2021" "2021-07-30"
+.TH LAMMPS "31 August 2021" "2021-08-31"
 .SH NAME
 .B LAMMPS
 \- Molecular Dynamics Simulator.

doc/src/Build.rst

@@ -22,4 +22,5 @@ page.
    Build_extras
    Build_manual
    Build_windows
+   Build_diskspace
    Build_development

doc/src/Build_development.rst

@@ -1,15 +1,15 @@
-Development build options (CMake only)
-======================================
+Development build options
+=========================
 
-The CMake build procedure of LAMMPS offers a few extra options which are
+The build procedures in LAMMPS offers a few extra options which are
 useful during development, testing or debugging.
 
 ----------
 
 .. _compilation:
 
-Monitor compilation flags
--------------------------
+Monitor compilation flags (CMake only)
+--------------------------------------
 
 Sometimes it is necessary to verify the complete sequence of compilation flags
 generated by the CMake build. To enable a more verbose output during
@@ -30,8 +30,8 @@ variable VERBOSE set to 1:
 
 .. _clang-tidy:
 
-Enable static code analysis with clang-tidy
--------------------------------------------
+Enable static code analysis with clang-tidy (CMake only)
+--------------------------------------------------------
 
 The `clang-tidy tool <https://clang.llvm.org/extra/clang-tidy/>`_ is a
 static code analysis tool to diagnose (and potentially fix) typical
@@ -52,8 +52,8 @@ significantly more time consuming than the compilation itself.
 
 .. _iwyu_processing:
 
-Report missing and unneeded '#include' statements
--------------------------------------------------
+Report missing and unneeded '#include' statements (CMake only)
+--------------------------------------------------------------
 
 The conventions for how and when to use and order include statements in
 LAMMPS are `documented in a separate file <https://github.com/lammps/lammps/blob/master/doc/include-file-conventions.md>`_
@@ -86,8 +86,8 @@ on recording all commands required to do the compilation.
 
 .. _sanitizer:
 
-Address, Undefined Behavior, and Thread Sanitizer Support
----------------------------------------------------------
+Address, Undefined Behavior, and Thread Sanitizer Support (CMake only)
+----------------------------------------------------------------------
 
 Compilers such as GCC and Clang support generating instrumented binaries
 which use different sanitizer libraries to detect problems in the code
@@ -116,8 +116,8 @@ compilation and linking stages.  This is done through setting the
 
 .. _testing:
 
-Code Coverage and Unit Testing
-------------------------------
+Code Coverage and Unit Testing (CMake only)
+-------------------------------------------
 
 The LAMMPS code is subject to multiple levels of automated testing
 during development: integration testing (i.e. whether the code compiles
@@ -464,7 +464,8 @@ Coding style utilities
 
 To aid with enforcing some of the coding style conventions in LAMMPS
 some additional build targets have been added. These require Python 3.5
-or later and will only work on Unix-like operating and file systems.
+or later and will only work properly on Unix-like operating and file systems.
+
 The following options are available.
 
 .. code-block:: bash
@@ -476,17 +477,43 @@ The following options are available.
    make check-permissions   # search for files with permissions issues
    make fix-permissions     # correct permissions issues in files
 
+These should help to replace all TAB characters with blanks and remove
+any trailing whitespace.  Also all LAMMPS homepage URL references can be
+updated to the location change from Sandia to the lammps.org domain.
+And the permission check can remove executable permissions from non-executable
+files (like source code).
+
+Clang-format support
+--------------------
+
 For the code in the ``unittest`` and ``src`` trees we are transitioning
 to use the `clang-format` tool to assist with having a consistent source
-code style.  The `clang-format` command bundled with Clang version 8.0
-or later is required.  The configuration is in files ``.clang-format``
-in the respective folders.  Since the modifications from `clang-format`
-can be significant and - especially for "legacy style code" - also is
-not always improving readability, a large number of files currently have
-a ``// clang-format off`` at the top, which will disable the processing.
-Over time, files will be refactored and updated to that `clang-format`
-may be applied to them (at least in part).
+code formatting style.  The `clang-format` command bundled with Clang
+version 8.0 or later is required.  The configuration is in files called
+``.clang-format`` in the respective folders.  Since the modifications
+from `clang-format` can be significant and - especially for "legacy
+style code" - they are not always improving readability, a large number
+of files currently have a ``// clang-format off`` at the top, which will
+disable the processing.  As of fall 2021 all files have been either
+"protected" this way or are enabled for full or partial `clang-format`
+processing.  Over time, the "protected" files will be refactored and
+updated so that `clang-format` may be applied to them as well.
 
-If `clang-format` is available, the source code files in the ``unittest``
-tree can be updated to conform to the formatting settings using
-``make format-tests`` and the files in ``src`` with ``make format-src``.
+It is recommended for all newly contributed files to use the clang-format
+processing while writing the code or do the coding style processing
+(including the scripts mentioned in the previous paragraph)
+
+If `clang-format` is available, files can be updated individually with
+commands like the following:
+
+.. code-block:: bash
+
+   $ clang-format -i some_file.cpp
+
+
+The following target are available for both, GNU make and CMake:
+
+.. code-block:: bash
+
+   make format-src       # apply clang-format to all files in src and the package folders
+   make format-tests     # apply clang-format to all files in the unittest tree

doc/src/Build_diskspace.rst

@@ -0,0 +1,45 @@
+Notes for saving disk space when building LAMMPS from source
+------------------------------------------------------------
+
+LAMMPS is a large software project with a large number of source files,
+extensive documentation, and a large collection of example files.
+When downloading LAMMPS by cloning the
+`git repository from GitHub <https://github.com/lammps/lammps>`_ this
+will by default also download the entire commit history since September 2006.
+Compiling LAMMPS will add the storage requirements of the compiled object
+files and libraries to the tally.
+
+In a user account on an HPC cluster with filesystem quotas or in other
+environments with restricted disk space capacity it may be needed to
+reduce the storage requirements. Here are some suggestions:
+
+- Create a so-called shallow repository by cloning only the last commit
+  instead of the full project history by using ``git clone git@github.com:lammps/lammps --depth=1 --branch=master``.
+  This reduces the downloaded size to about half.  With ``--depth=1`` it is not possible to check out different
+  versions/branches of LAMMPS, using ``--depth=1000`` will make multiple recent versions available at little
+  extra storage needs (the entire git history had nearly 30,000 commits in fall 2021).
+
+- Download a tar archive from either the `download section on the LAMMPS homepage <https://www.lammps.org/download.html>`_
+  or from the `LAMMPS releases page on GitHub <https://github.com/lammps/lammps/releases>`_ these will not
+  contain the git history at all.
+
+- Build LAMMPS without the debug flag (remove ``-g`` from the machine makefile or use ``-DCMAKE_BUILD_TYPE=Release``)
+  or use the ``strip`` command on the LAMMPS executable when no more debugging would be needed.  The strip command
+  may also be applied to the LAMMPS shared library. The static library may be deleted entirely.
+
+- Delete compiled object files and libraries after copying the LAMMPS executable to a permanent location.
+  When using the traditional build process, one may use ``make clean-<machine>`` or ``make clean-all``
+  to delete object files in the src folder.  For CMake based builds, one may use ``make clean`` or just
+  delete the entire build folder.
+
+- The folders containing the documentation tree (doc), the examples (examples) are not needed to build and
+  run LAMMPS and can be safely deleted.  Some files in the potentials folder are large and may be deleted,
+  if not needed.  The largest of those files (occupying about 120 MBytes combined) will only be downloaded on
+  demand, when the corresponding package is installed.
+
+- When using the CMake build procedure, the compilation can be done on a (local) scratch storage that will not
+  count toward the quota.  A local scratch file system may offer the additional benefit of speeding up creating
+  object files and linking with libraries compared to a networked file system.  Also with CMake (and unlike with
+  the traditional make) it is possible to compile LAMMPS executables with different settings and packages included
+  from the same source tree since all the configuration information is stored in the build folder.  So it is
+  not necessary to have multiple copies of LAMMPS.

doc/src/Commands_compute.rst

@@ -152,7 +152,7 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`temp/chunk <compute_temp_chunk>`
    * :doc:`temp/com <compute_temp_com>`
    * :doc:`temp/cs <compute_temp_cs>`
-   * :doc:`temp/deform <compute_temp_deform>`
+   * :doc:`temp/deform (k) <compute_temp_deform>`
    * :doc:`temp/deform/eff <compute_temp_deform_eff>`
    * :doc:`temp/drude <compute_temp_drude>`
    * :doc:`temp/eff <compute_temp_eff>`

doc/src/Commands_fix.rst

@@ -148,7 +148,7 @@ OPT.
    * :doc:`nvt/body <fix_nvt_body>`
    * :doc:`nvt/eff <fix_nh_eff>`
    * :doc:`nvt/manifold/rattle <fix_nvt_manifold_rattle>`
-   * :doc:`nvt/sllod (io) <fix_nvt_sllod>`
+   * :doc:`nvt/sllod (iko) <fix_nvt_sllod>`
    * :doc:`nvt/sllod/eff <fix_nvt_sllod_eff>`
    * :doc:`nvt/sphere (o) <fix_nvt_sphere>`
    * :doc:`nvt/uef <fix_nh_uef>`
@@ -236,6 +236,7 @@ OPT.
    * :doc:`ti/spring <fix_ti_spring>`
    * :doc:`tmd <fix_tmd>`
    * :doc:`ttm <fix_ttm>`
+   * :doc:`ttm/grid <fix_ttm>`
    * :doc:`ttm/mod <fix_ttm>`
    * :doc:`tune/kspace <fix_tune_kspace>`
    * :doc:`vector <fix_vector>`

doc/src/Commands_kspace.rst

@@ -24,6 +24,7 @@ OPT.
 
    * :doc:`ewald (o) <kspace_style>`
    * :doc:`ewald/disp <kspace_style>`
+   * :doc:`ewald/disp/dipole <kspace_style>`
    * :doc:`ewald/dipole <kspace_style>`
    * :doc:`ewald/dipole/spin <kspace_style>`
    * :doc:`msm (o) <kspace_style>`

doc/src/Commands_pair.rst

@@ -75,6 +75,7 @@ OPT.
    * :doc:`coul/debye (gko) <pair_coul>`
    * :doc:`coul/diel (o) <pair_coul_diel>`
    * :doc:`coul/dsf (gko) <pair_coul>`
+   * :doc:`coul/exclude <pair_coul>`
    * :doc:`coul/long (gko) <pair_coul>`
    * :doc:`coul/long/cs (g) <pair_cs>`
    * :doc:`coul/long/dielectric <pair_dielectric>`

doc/src/Developer_utils.rst

@@ -334,10 +334,11 @@ arguments of commands in LAMMPS are parsed and to make abstractions of
 repetitive tasks.
 
 The :cpp:class:`LAMMPS_NS::ArgInfo` class provides an abstraction
-for parsing references to compute or fix styles or variables. These
-would start with a "c\_", "f\_", "v\_" followed by the ID or name of
-than instance and may be postfixed with one or two array indices
-"[<number>]" with numbers > 0.
+for parsing references to compute or fix styles, variables or custom
+integer or double properties handled by :doc:`fix property/atom <fix_property_atom>`.
+These would start with a "c\_", "f\_", "v\_", "d\_", "d2\_", "i\_", or "i2\_"
+followed by the ID or name of than instance and may be postfixed with
+one or two array indices "[<number>]" with numbers > 0.
 
 A typical code segment would look like this:
 

doc/src/Howto.rst

@@ -54,6 +54,7 @@ Analysis howto
    Howto_kappa
    Howto_viscosity
    Howto_diffusion
+   Howto_structured_data
 
 Force fields howto
 ==================

doc/src/Howto_structured_data.rst

@@ -0,0 +1,154 @@
+Output structured data from LAMMPS
+##################################
+
+LAMMPS can output structured data with the :doc:`print <print>` and :doc:`fix
+print <fix_print>` command.  This gives you flexibility since you can build
+custom data formats that contain system properties, thermo data, and variables
+values. This output can be directed to the screen and/or to a file for post
+processing.
+
+Writing the current system state, thermo data, variable values
+==============================================================
+
+Use the :doc:`print <print>` command to output the current system state, which
+can include system properties, thermo data and variable values.
+
+YAML
+----
+
+.. code-block:: LAMMPS
+
+   print """---
+   timestep: $(step)
+   pe: $(pe)
+   ke: $(ke)""" file current_state.yaml screen no
+
+.. code-block:: yaml
+   :caption: current_state.yaml
+
+   ---
+   timestep: 250
+   pe: -4.7774327356321810711
+   ke: 2.4962152903997174569
+
+JSON
+----
+
+.. code-block:: LAMMPS
+
+   print """{
+     "timestep": $(step),
+     "pe": $(pe),
+     "ke": $(ke)
+   }""" file current_state.json screen no
+
+.. code-block:: JSON
+   :caption: current_state.json
+
+   {
+     "timestep": 250,
+     "pe": -4.7774327356321810711,
+     "ke": 2.4962152903997174569
+   }
+
+
+Writing continuous data during a simulation
+===========================================
+
+The :doc:`fix print <fix_print>` command allows you to output an arbitrary string at defined times during a simulation run.
+
+YAML
+----
+
+.. code-block:: LAMMPS
+
+   fix extra all print 50 """
+   - timestep: $(step)
+     pe: $(pe)
+     ke: $(ke)""" file output.yaml screen no
+
+.. code-block:: yaml
+   :caption: output.yaml
+
+   # Fix print output for fix extra
+   - timestep: 0
+     pe: -6.77336805325924729
+     ke: 4.4988750000000026219
+
+   - timestep: 50
+     pe: -4.8082494418323200591
+     ke: 2.5257981827119797558
+
+   - timestep: 100
+     pe: -4.7875608875581505686
+     ke: 2.5062598821985102582
+
+   - timestep: 150
+     pe: -4.7471033686005483787
+     ke: 2.466095925545450207
+
+   - timestep: 200
+     pe: -4.7509052858544134068
+     ke: 2.4701136792591693592
+
+   - timestep: 250
+     pe: -4.7774327356321810711
+     ke: 2.4962152903997174569
+
+Post-processing of YAML files can be easily be done with Python and other
+scripting languages. In case of Python the `yaml` package allows you to load the
+data files and obtain a list of dictionaries.
+
+.. code-block:: python
+
+   import yaml
+
+   with open("output.yaml") as f:
+      data = yaml.load(f, Loader=yaml.FullLoader)
+
+   print(data)
+
+.. code-block::
+
+   [{'timestep': 0, 'pe': -6.773368053259247, 'ke': 4.498875000000003}, {'timestep': 50, 'pe': -4.80824944183232, 'ke': 2.5257981827119798}, {'timestep': 100, 'pe': -4.787560887558151, 'ke': 2.5062598821985103}, {'timestep': 150, 'pe': -4.747103368600548, 'ke': 2.46609592554545}, {'timestep': 200, 'pe': -4.750905285854413, 'ke': 2.4701136792591694}, {'timestep': 250, 'pe': -4.777432735632181, 'ke': 2.4962152903997175}]
+
+Line Delimited JSON (LD-JSON)
+-----------------------------
+
+The JSON format itself is very strict when it comes to delimiters. For continuous
+output/streaming data it is beneficial use the *line delimited JSON* format.
+Each line represents one JSON object.
+
+.. code-block:: LAMMPS
+
+   fix extra all print 50 """{"timestep": $(step), "pe": $(pe), "ke": $(ke)}""" title "" file output.json screen no
+
+.. code-block:: json
+   :caption: output.json
+
+   {"timestep": 0, "pe": -6.77336805325924729, "ke": 4.4988750000000026219}
+   {"timestep": 50, "pe": -4.8082494418323200591, "ke": 2.5257981827119797558}
+   {"timestep": 100, "pe": -4.7875608875581505686, "ke": 2.5062598821985102582}
+   {"timestep": 150, "pe": -4.7471033686005483787, "ke": 2.466095925545450207}
+   {"timestep": 200, "pe": -4.7509052858544134068, "ke": 2.4701136792591693592}
+   {"timestep": 250, "pe": -4.7774327356321810711, "ke": 2.4962152903997174569}
+
+One simple way to load this data into a Python script is to use the `pandas`
+package. It can directly load these files into a data frame:
+
+.. code-block:: python
+
+   import pandas as pd
+
+   data = pd.read_json('output.json', lines=True)
+   print(data)
+
+.. code-block:: bash
+
+      timestep        pe        ke
+   0         0 -6.773368  4.498875
+   1        50 -4.808249  2.525798
+   2       100 -4.787561  2.506260
+   3       150 -4.747103  2.466096
+   4       200 -4.750905  2.470114
+   5       250 -4.777433  2.496215

doc/src/Intro_authors.rst

@@ -29,7 +29,7 @@ The following folks deserve special recognition.  Many of the packages
 they have written are unique for an MD code and LAMMPS would not be as
 general-purpose as it is without their expertise and efforts.
 
-* Metin Aktulga (MSU), REAXFF package for C version of ReaxFF
+* Metin Aktulga (MSU), REAXFF package for C/C++ version of ReaxFF
 * Mike Brown (Intel), GPU and INTEL packages
 * Colin Denniston (U Western Ontario), LATBOLTZ package
 * Georg Ganzenmuller (EMI), MACHDYN and SPH packages
@@ -37,9 +37,10 @@ general-purpose as it is without their expertise and efforts.
 * Reese Jones (Sandia) and colleagues, ATC package for atom/continuum coupling
 * Christoph Kloss (DCS Computing), LIGGGHTS code for granular materials, built on top of LAMMPS
 * Rudra Mukherjee (JPL), POEMS package for articulated rigid body motion
-* Trung Ngyuen (Northwestern U), GPU and RIGID and BODY packages
+* Trung Ngyuen (Northwestern U), GPU, RIGID, BODY, and DIELECTRIC packages
 * Mike Parks (Sandia), PERI package for Peridynamics
 * Roy Pollock (LLNL), Ewald and PPPM solvers
+* Julien Tranchida (Sandia), SPIN package
 * Christian Trott (Sandia), CUDA and KOKKOS packages
 * Ilya Valuev (JIHT), AWPMD package for wave packet MD
 * Greg Wagner (Northwestern U), MEAM package for MEAM potential

doc/src/Intro_features.rst

@@ -27,19 +27,19 @@ General features
 * distributed memory message-passing parallelism (MPI)
 * shared memory multi-threading parallelism (OpenMP)
 * spatial decomposition of simulation domain for MPI parallelism
-* particle decomposition inside of spatial decomposition for OpenMP parallelism
+* particle decomposition inside of spatial decomposition for OpenMP and GPU parallelism
 * GPLv2 licensed open-source distribution
 * highly portable C++-11
 * modular code with most functionality in optional packages
-* only depends on MPI library for basic parallel functionality
+* only depends on MPI library for basic parallel functionality, MPI stub for serial compilation
 * other libraries are optional and only required for specific packages
-* GPU (CUDA and OpenCL), Intel Xeon Phi, and OpenMP support for many code features
+* GPU (CUDA, OpenCL, HIP, SYCL), Intel Xeon Phi, and OpenMP support for many code features
 * easy to extend with new features and functionality
 * runs from an input script
 * syntax for defining and using variables and formulas
 * syntax for looping over runs and breaking out of loops
 * run one or multiple simulations simultaneously (in parallel) from one script
-* build as library, invoke LAMMPS through library interface or provided Python wrapper
+* build as library, invoke LAMMPS through library interface or provided Python wrapper or SWIG based wrappers
 * couple with other codes: LAMMPS calls other code, other code calls LAMMPS, umbrella code calls both
 
 .. _particle:
@@ -57,9 +57,11 @@ Particle and model types
 * granular materials
 * coarse-grained mesoscale models
 * finite-size spherical and ellipsoidal particles
-* finite-size  line segment (2d) and triangle (3d) particles
+* finite-size line segment (2d) and triangle (3d) particles
+* finite-size rounded polygons (2d) and polyhedra (3d) particles
 * point dipole particles
-* rigid collections of particles
+* particles with magnetic spin
+* rigid collections of n particles
 * hybrid combinations of these
 
 .. _ff:
@@ -74,24 +76,28 @@ commands)
 
 * pairwise potentials: Lennard-Jones, Buckingham, Morse, Born-Mayer-Huggins, Yukawa, soft, class 2 (COMPASS), hydrogen bond, tabulated
 * charged pairwise potentials: Coulombic, point-dipole
-* many-body potentials: EAM, Finnis/Sinclair EAM, modified EAM (MEAM), embedded ion method (EIM), EDIP, ADP, Stillinger-Weber, Tersoff, REBO, AIREBO, ReaxFF, COMB, SNAP, Streitz-Mintmire, 3-body polymorphic
-* long-range interactions for charge, point-dipoles, and LJ dispersion:     Ewald, Wolf, PPPM (similar to particle-mesh Ewald)
+* many-body potentials: EAM, Finnis/Sinclair EAM, modified EAM (MEAM), embedded ion method (EIM), EDIP, ADP, Stillinger-Weber, Tersoff, REBO, AIREBO, ReaxFF, COMB, Streitz-Mintmire, 3-body polymorphic, BOP, Vashishta
+* machine learning potentials: SNAP, GAP, ACE, N2P2, RANN, AGNI
+* long-range interactions for charge, point-dipoles, and LJ dispersion:  Ewald, Wolf, PPPM (similar to particle-mesh Ewald), MSM
 * polarization models: :doc:`QEq <fix_qeq>`,     :doc:`core/shell model <Howto_coreshell>`,     :doc:`Drude dipole model <Howto_drude>`
 * charge equilibration (QEq via dynamic, point, shielded, Slater methods)
 * coarse-grained potentials: DPD, GayBerne, REsquared, colloidal, DLVO
-* mesoscopic potentials: granular, Peridynamics, SPH
+* mesoscopic potentials: granular, Peridynamics, SPH, mesoscopic tubular potential (MESONT)
+* semi-empirical potentials: multi-ion generalized pseudopotential theory (MGPT), second moment tight binding + QEq (SMTB-Q), density functional tight-binding (LATTE)
 * electron force field (eFF, AWPMD)
-* bond potentials: harmonic, FENE, Morse, nonlinear, class 2,     quartic (breakable)
-* angle potentials: harmonic, CHARMM, cosine, cosine/squared, cosine/periodic,     class 2 (COMPASS)
-* dihedral potentials: harmonic, CHARMM, multi-harmonic, helix,     class 2 (COMPASS), OPLS
-* improper potentials: harmonic, cvff, umbrella, class 2 (COMPASS)
+* bond potentials: harmonic, FENE, Morse, nonlinear, class 2, quartic (breakable), tabulated
+* angle potentials: harmonic, CHARMM, cosine, cosine/squared, cosine/periodic, class 2 (COMPASS), tabulated
+* dihedral potentials: harmonic, CHARMM, multi-harmonic, helix, class 2 (COMPASS), OPLS, tabulated
+* improper potentials: harmonic, cvff, umbrella, class 2 (COMPASS), tabulated
 * polymer potentials: all-atom, united-atom, bead-spring, breakable
-* water potentials: TIP3P, TIP4P, SPC
+* water potentials: TIP3P, TIP4P, SPC, SPC/E and variants
+* interlayer potentials for graphene and analogues
+* metal-organic framework potentials (QuickFF, MO-FF)
 * implicit solvent potentials: hydrodynamic lubrication, Debye
-* force-field compatibility with common CHARMM, AMBER, DREIDING,     OPLS, GROMACS, COMPASS options
+* force-field compatibility with common CHARMM, AMBER, DREIDING, OPLS, GROMACS, COMPASS options
 * access to the `OpenKIM Repository <http://openkim.org>`_ of potentials via     :doc:`kim command <kim_commands>`
-* hybrid potentials: multiple pair, bond, angle, dihedral, improper     potentials can be used in one simulation
-* overlaid potentials: superposition of multiple pair potentials
+* hybrid potentials: multiple pair, bond, angle, dihedral, improper potentials can be used in one simulation
+* overlaid potentials: superposition of multiple pair potentials (including many-body) with optional scale factor
 
 .. _create:
 
@@ -124,9 +130,10 @@ Ensembles, constraints, and boundary conditions
 * harmonic (umbrella) constraint forces
 * rigid body constraints
 * SHAKE bond and angle constraints
-* Monte Carlo bond breaking, formation, swapping
+* motion constraints to manifold surfaces
+* Monte Carlo bond breaking, formation, swapping, template based reaction modeling
 * atom/molecule insertion and deletion
-* walls of various kinds
+* walls of various kinds, static and moving
 * non-equilibrium molecular dynamics (NEMD)
 * variety of additional boundary conditions and constraints
 
@@ -150,6 +157,7 @@ Diagnostics
 ^^^^^^^^^^^
 
 * see various flavors of the :doc:`fix <fix>` and :doc:`compute <compute>` commands
+* introspection command for system, simulation, and compile time settings and configurations
 
 .. _output:
 
@@ -164,8 +172,9 @@ Output
 * parallel I/O of dump and restart files
 * per-atom quantities (energy, stress, centro-symmetry parameter, CNA, etc)
 * user-defined system-wide (log file) or per-atom (dump file) calculations
-* spatial and time averaging of per-atom quantities
-* time averaging of system-wide quantities
+* custom partitioning (chunks) for binning, and static or dynamic grouping of atoms for analysis
+* spatial, time, and per-chunk averaging of per-atom quantities
+* time averaging and histogramming of system-wide quantities
 * atom snapshots in native, XYZ, XTC, DCD, CFG formats
 
 .. _replica1:
@@ -178,7 +187,7 @@ Multi-replica models
 * :doc:`parallel replica dynamics <prd>`
 * :doc:`temperature accelerated dynamics <tad>`
 * :doc:`parallel tempering <temper>`
-* :doc:`path-integral MD <fix_pimd>`
+* path-integral MD: `first variant <fix_pimd>`, `second variant <fix_ipi>`
 * multi-walker collective variables with :doc:`Colvars <fix_colvars>` and :doc:`Plumed <fix_plumed>`
 
 .. _prepost:
@@ -210,11 +219,12 @@ page for details.
 These are LAMMPS capabilities which you may not think of as typical
 classical MD options:
 
-* :doc:`static <balance>` and :doc:`dynamic load-balancing <fix_balance>`
+* :doc:`static <balance>` and :doc:`dynamic load-balancing <fix_balance>`, optional with recursive bisectioning decomposition
 * :doc:`generalized aspherical particles <Howto_body>`
 * :doc:`stochastic rotation dynamics (SRD) <fix_srd>`
-* :doc:`real-time visualization and interactive MD <fix_imd>`
+* :doc:`real-time visualization and interactive MD <fix_imd>`, :doc:`built-in renderer for images and movies <dump_image>`
 * calculate :doc:`virtual diffraction patterns <compute_xrd>`
+* calculate :doc:`finite temperature phonon dispersion <fix_phonon>` and the :doc:`dynamical matrix of minimized structures <dynamical_matrix>`
 * :doc:`atom-to-continuum coupling <fix_atc>` with finite elements
 * coupled rigid body integration via the :doc:`POEMS <fix_poems>` library
 * :doc:`QM/MM coupling <fix_qmmm>`

doc/src/Intro_opensource.rst

@@ -1,40 +1,61 @@
 LAMMPS open-source license
 --------------------------
 
-LAMMPS is a freely-available open-source code, distributed under the
-terms of the `GNU Public License Version 2 <gpl_>`_, which means you can
-use or modify the code however you wish for your own purposes, but have
-to adhere to certain rules when redistributing it or software derived
+GPL version of LAMMPS
+^^^^^^^^^^^^^^^^^^^^^
+
+LAMMPS is an open-source code, available free-of-charge, and distributed
+under the terms of the `GNU Public License Version 2 <gpl_>`_ (GPLv2),
+which means you can use or modify the code however you wish for your own
+purposes, but have to adhere to certain rules when redistributing it -
+specifically in binary form - or are distributing software derived
 from it or that includes parts of it.
 
-LAMMPS comes with no warranty of any kind.  As each source file states
-in its header, it is a copyrighted code that is distributed free-of-
-charge, under the terms of the `GNU Public License Version 2 <gpl_>`_
-(GPLv2).  This is often referred to as open-source distribution - see
-`www.gnu.org <gnuorg_>`_ or `www.opensource.org <opensource_>`_.  The
-legal text of the GPL is in the LICENSE file included in the LAMMPS
-distribution.
+LAMMPS comes with no warranty of any kind.
+
+As each source file states in its header, it is a copyrighted code, and
+thus not in the public domain. For more information about open-source
+software and open-source distribution, see `www.gnu.org <gnuorg_>`_
+or `www.opensource.org <opensource_>`_.  The legal text of the GPL as it
+applies to LAMMPS is in the LICENSE file included in the LAMMPS distribution.
 
 .. _gpl: https://github.com/lammps/lammps/blob/master/LICENSE
 
+.. _lgpl: https://www.gnu.org/licenses/old-licenses/lgpl-2.1.html
+
 .. _gnuorg: http://www.gnu.org
 
 .. _opensource: http://www.opensource.org
 
-Here is a summary of what the GPL means for LAMMPS users:
+Here is a more specific summary of what the GPL means for LAMMPS users:
 
-(1) Anyone is free to use, modify, or extend LAMMPS in any way they
+(1) Anyone is free to use, copy, modify, or extend LAMMPS in any way they
 choose, including for commercial purposes.
 
 (2) If you **distribute** a modified version of LAMMPS, it must remain
-open-source, meaning you distribute **all** of it under the terms of
-the GPL.  You should clearly annotate such a code as a derivative version
-of LAMMPS.
+open-source, meaning you are required to distribute **all** of it under
+the terms of the GPL.  You should clearly annotate such a modified code
+as a derivative version of LAMMPS.
 
 (3) If you release any code that includes or uses LAMMPS source code,
 then it must also be open-sourced, meaning you distribute it under
-the terms of the GPL.
+the terms of the GPL.  You may write code that interfaces LAMMPS to
+a differently licensed library.  In that case the code that provides
+the interface must be licensed GPL, but not necessarily that library
+unless you are distributing binaries that require the library to run.
 
 (4) If you give LAMMPS files to someone else, the GPL LICENSE file and
 source file headers (including the copyright and GPL notices) should
 remain part of the code.
+
+
+LGPL version of LAMMPS
+^^^^^^^^^^^^^^^^^^^^^^
+
+We occasionally make stable LAMMPS releases available under the `GNU
+Lesser Public License v2.1 <lgpl_>`_.  This is on request only and with
+non-LGPL compliant files removed.  This allows uses linking non-GPL
+compatible software with the (otherwise unmodified) LAMMPS library
+or loading it dynamically at runtime.  Any **modifications** to
+the LAMMPS code however, even with the LGPL licensed version, must still
+be made available under the same open source terms as LAMMPS itself.

doc/src/Intro_overview.rst

@@ -10,24 +10,26 @@ conditions.  It can model 2d or 3d systems with only a few particles
 up to millions or billions.
 
 LAMMPS can be built and run on a laptop or desktop machine, but is
-designed for parallel computers.  It will run on any parallel machine
-that supports the `MPI <mpi_>`_ message-passing library.  This includes
-shared-memory boxes and distributed-memory clusters and
-supercomputers.
+designed for parallel computers.  It will run in serial and on any
+parallel machine that supports the `MPI <mpi_>`_ message-passing
+library.  This includes shared-memory boxes and distributed-memory
+clusters and supercomputers. Parts of LAMMPS also support
+`OpenMP multi-threading <omp_>`_, vectorization and GPU acceleration.
 
 .. _mpi: https://en.wikipedia.org/wiki/Message_Passing_Interface
 .. _lws: https://www.lammps.org
+.. _omp: https://www.openmp.org
 
 LAMMPS is written in C++ and requires a compiler that is at least
-compatible with the C++-11 standard.
-Earlier versions were written in F77 and F90.  See the `History page
+compatible with the C++-11 standard.  Earlier versions were written in
+F77, F90, and C++-98.  See the `History page
 <https://www.lammps.org/history.html>`_ of the website for details.  All
-versions can be downloaded from the `LAMMPS website <lws_>`_.
+versions can be downloaded as source code from the `LAMMPS website
+<lws_>`_.
 
-LAMMPS is designed to be easy to modify or extend with new
-capabilities, such as new force fields, atom types, boundary
-conditions, or diagnostics.  See the :doc:`Modify <Modify>` page for
-more details.
+LAMMPS is designed to be easy to modify or extend with new capabilities,
+such as new force fields, atom types, boundary conditions, or
+diagnostics.  See the :doc:`Modify <Modify>` page for more details.
 
 In the most general sense, LAMMPS integrates Newton's equations of
 motion for a collection of interacting particles.  A single particle
@@ -47,4 +49,5 @@ MPI parallelization to partition the simulation domain into small
 sub-domains of equal computational cost, one of which is assigned to
 each processor.  Processors communicate and store "ghost" atom
 information for atoms that border their sub-domain.  Multi-threading
-parallelization with with particle-decomposition can be used in addition.
+parallelization and GPU acceleration with with particle-decomposition
+can be used in addition.

doc/src/Manual_version.rst

@@ -2,12 +2,21 @@ What does a LAMMPS version mean
 -------------------------------
 
 The LAMMPS "version" is the date when it was released, such as 1 May
-2014. LAMMPS is updated continuously.  Whenever we fix a bug or add a
-feature, we release it in the next *patch* release, which are
-typically made every couple of weeks.  Info on patch releases are on
-`this website page <https://www.lammps.org/bug.html>`_. Every few
-months, the latest patch release is subjected to more thorough testing
-and labeled as a *stable* version.
+2014.  LAMMPS is updated continuously and we aim to keep it working
+correctly and reliably at all times.  You can follow its development
+in a public `git repository on GitHub <https://github.com/lammps/lammps>`_.
+
+Whenever we fix a bug or update or add a feature, it will be merged into
+the `master` branch of the git repository.  When a sufficient number of
+changes have accumulated *and* the software passes a set of automated
+tests, we release it in the next *patch* release, which are made every
+few weeks.  Info on patch releases are on `this website page
+<https://www.lammps.org/bug.html>`_.
+
+Once or twice a year, only bug fixes and small, non-intrusive changes are
+included for a period of time, and the code is subjected to more detailed
+and thorough testing than the default automated testing.  The latest
+patch release after such a period is then labeled as a *stable* version.
 
 Each version of LAMMPS contains all the features and bug-fixes up to
 and including its version date.

doc/src/Tools.rst

@@ -172,7 +172,7 @@ Chris Lorenz (chris.lorenz at kcl.ac.uk), King's College London.
 chain tool
 ----------------------
 
-The file chain.f creates a LAMMPS data file containing bead-spring
+The file chain.f90 creates a LAMMPS data file containing bead-spring
 polymer chains and/or monomer solvent atoms.  It uses a text file
 containing chain definition parameters as an input.  The created
 chains and solvent atoms can strongly overlap, so LAMMPS needs to run

doc/src/angle_cosine_shift.rst

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

doc/src/angle_cosine_shift_exp.rst

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

doc/src/angle_fourier.rst

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

doc/src/angle_fourier_simple.rst

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

doc/src/angle_gaussian.rst

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

doc/src/angle_quartic.rst

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

doc/src/balance.rst

@@ -477,12 +477,18 @@ atom-style variables can reference the position of a particle, its
 velocity, the volume of its Voronoi cell, etc.
 
 The *store* weight style does not compute a weight factor.  Instead it
-stores the current accumulated weights in a custom per-atom property
-specified by *name*\ .  This must be a property defined as *d_name* via
-the :doc:`fix property/atom <fix_property_atom>` command.  Note that
-these custom per-atom properties can be output in a :doc:`dump <dump>`
-file, so this is a way to examine, debug, or visualize the
-per-particle weights computed during the load-balancing operation.
+stores the current accumulated weights in a custom per-atom vector
+specified by *name*\ .  This must be a vector defined as *d_name* via
+the :doc:`fix property/atom <fix_property_atom>` command.  This means
+the values in the vector can be read as part of a data file with the
+:doc:`read_data <read_data>` command or specified with the :doc:`set
+<set>` command.  These weights can also be output in a :doc:`dump
+<dump>` file, so this is a way to examine, debug, or visualize the
+per-particle weights used during the load-balancing operation.
+
+Note that the name of the custom per-atom vector is specified just
+as *name*, not as *d_name* as it is for other commands that use
+different kinds of custom atom vectors or arrays as arguments.
 
 ----------
 

doc/src/bond_gaussian.rst

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

doc/src/bond_harmonic_shift.rst

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

doc/src/bond_harmonic_shift_cut.rst

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

doc/src/compute_property_atom.rst

@@ -20,7 +20,8 @@ Syntax
                              x, y, z, xs, ys, zs, xu, yu, zu, ix, iy, iz,
                              vx, vy, vz, fx, fy, fz,
                              q, mux, muy, muz, mu,
-                             sp, spx, spy, spz, fmx, fmy, fmz,
+                             spx, spy, spz, sp, fmx, fmy, fmz,
+                             nbonds,
                              radius, diameter, omegax, omegay, omegaz,
                              angmomx, angmomy, angmomz,
                              shapex,shapey, shapez,
@@ -29,42 +30,42 @@ Syntax
                              corner1x, corner1y, corner1z,
                              corner2x, corner2y, corner2z,
                              corner3x, corner3y, corner3z,
-                             nbonds,
-                             buckling,
-                             vfrac, s0,
-                             spin, eradius, ervel, erforce,
-                             rho, drho, e, de, cv,
-                             i_name, d_name
+                             i_name, d_name, i2_name[I], d2_name[I],
+                             vfrac, s0, spin, eradius, ervel, erforce,
+                             rho, drho, e, de, cv, buckling,
 
   .. parsed-literal::
 
-           id = atom ID
-           mol = molecule ID
-           proc = ID of processor that owns atom
-           type = atom type
-           mass = atom mass
-           x,y,z = unscaled atom coordinates
-           xs,ys,zs = scaled atom coordinates
-           xu,yu,zu = unwrapped atom coordinates
-           ix,iy,iz = box image that the atom is in
-           vx,vy,vz = atom velocities
-           fx,fy,fz = forces on atoms
-           q = atom charge
-           mux,muy,muz = orientation of dipole moment of atom
-           mu = magnitude of dipole moment of atom
-           sp = atomic magnetic spin moment
-           spx, spy, spz = direction of the atomic magnetic spin
-           fmx, fmy, fmz = magnetic force
-           radius,diameter = radius,diameter of spherical particle
-           omegax,omegay,omegaz = angular velocity of spherical particle
-           angmomx,angmomy,angmomz = angular momentum of aspherical particle
-           shapex,shapey,shapez = 3 diameters of aspherical particle
-           quatw,quati,quatj,quatk = quaternion components for aspherical or body particles
-           tqx,tqy,tqz = torque on finite-size particles
-           end12x, end12y, end12z = end points of line segment
-           corner123x, corner123y, corner123z = corner points of triangle
-           nbonds = number of bonds assigned to an atom
-           buckling = buckling flag used in mesoscopic simulation of nanotubes
+           *id* = atom ID
+           *mol* = molecule ID
+           *proc* = ID of processor that owns atom
+           *type* = atom type
+           *mass* = atom mass
+           *x,y,z* = unscaled atom coordinates
+           *xs,ys,zs* = scaled atom coordinates
+           *xu,yu,zu* = unwrapped atom coordinates
+           *ix,iy,iz* = box image that the atom is in
+           *vx,vy,vz* = atom velocities
+           *fx,fy,fz* = forces on atoms
+           *q* = atom charge
+           *mux,muy,muz* = orientation of dipole moment of atom
+           *mu* = magnitude of dipole moment of atom
+           *spx, spy, spz* = direction of the atomic magnetic spin
+           *sp* = magintude of atomic magnetic spin moment
+           *fmx, fmy, fmz* = magnetic force
+           *nbonds* = number of bonds assigned to an atom
+           *radius,diameter* = radius,diameter of spherical particle
+           *omegax,omegay,omegaz* = angular velocity of spherical particle
+           *angmomx,angmomy,angmomz* = angular momentum of aspherical particle
+           *shapex,shapey,shapez* = 3 diameters of aspherical particle
+           *quatw,quati,quatj,quatk* = quaternion components for aspherical or body particles
+           *tqx,tqy,tqz* = torque on finite-size particles
+           *end12x, end12y, end12z* = end points of line segment
+           *corner123x, corner123y, corner123z* = corner points of triangle
+           *i_name* = custom integer vector with name
+           *d_name* = custom floating point vector with name
+           *i2_name[I]* = Ith column of custom integer array with name
+           *d2_name[I]* = Ith column of custom floating-point array with name
 
   .. parsed-literal::
 
@@ -91,9 +92,8 @@ Syntax
 
   .. parsed-literal::
 
-           :doc:`fix property/atom <fix_property_atom>` per-atom properties:
-           i_name = custom integer vector with name
-           d_name = custom integer vector with name
+           MESONT package per-atom properties:
+           buckling = buckling flag used in mesoscopic simulation of nanotubes
 
 Examples
 """"""""
@@ -104,6 +104,7 @@ Examples
    compute 2 all property/atom type
    compute 1 all property/atom ix iy iz
    compute 3 all property/atom sp spx spy spz
+   compute 1 all property/atom i_myFlag d_Sxyz[1] d_Sxyz[3]
 
 Description
 """""""""""
@@ -116,20 +117,37 @@ ave/atom <fix_ave_atom>`, :doc:`fix ave/histo <fix_ave_histo>`,
 :doc:`fix ave/chunk <fix_ave_chunk>`, and :doc:`atom-style variable
 <variable>` commands.
 
-The list of possible attributes is the same as that used by the
-:doc:`dump custom <dump>` command, which describes their meaning, with
-some additional quantities that are only defined for certain
-:doc:`atom styles <atom_style>`.  Basically, this augmented list gives
-an input script access to any per-atom quantity stored by LAMMPS.
+The list of possible attributes is essentially the same as that used
+by the :doc:`dump custom <dump>` command, which describes their
+meaning, with some additional quantities that are only defined for
+certain :doc:`atom styles <atom_style>`.  The goal of this augmented
+list gives an input script access to any per-atom quantity stored by
+LAMMPS.
 
 The values are stored in a per-atom vector or array as discussed
 below.  Zeroes are stored for atoms not in the specified group or for
 quantities that are not defined for a particular particle in the group
 (e.g. *shapex* if the particle is not an ellipsoid).
 
+Attributes *i_name*, *d_name*, *i2_name*, *d2_name* refer to custom
+per-atom integer and floating-point vectors or arrays that have been
+added via the :doc:`fix property/atom <fix_property_atom>` command.
+When that command is used specific names are given to each attribute
+which are the "name" portion of these attributes.  For arrays *i2_name*
+and *d2_name*, the column of the array must also be included following
+the name in brackets: e.g. d2_xyz[2], i2_mySpin[3].
+
 The additional quantities only accessible via this command, and not
 directly via the :doc:`dump custom <dump>` command, are as follows.
 
+*Nbonds* is available for all molecular atom styles and refers to the
+number of explicit bonds assigned to an atom.  Note that if the
+:doc:`newton bond <newton>` command is set to *on*\ , which is the
+default, then every bond in the system is assigned to only one of the
+two atoms in the bond.  Thus a bond between atoms I,J may be tallied
+for either atom I or atom J.  If :doc:`newton bond off <newton>` is
+set, it will be tallied with both atom I and atom J.
+
 *Shapex*, *shapey*, and *shapez* are defined for ellipsoidal particles
 and define the 3d shape of each particle.
 
@@ -146,19 +164,8 @@ line segment.
 *corner2z*, *corner3x*, *corner3y*, *corner3z*, are defined for
 triangular particles and define the corner points of each triangle.
 
-*Nbonds* is available for all molecular atom styles and refers to the
-number of explicit bonds assigned to an atom.  Note that if the
-:doc:`newton bond <newton>` command is set to *on*, which is the
-default, then every bond in the system is assigned to only one of the
-two atoms in the bond.  Thus a bond between atoms I,J may be tallied
-for either atom I or atom J.  If :doc:`newton bond off <newton>` is
-set, it will be tallied with both atom I and atom J.
-
-The *i_name* and *d_name* attributes refer to custom integer and
-floating-point properties that have been added to each atom via the
-:doc:`fix property/atom <fix_property_atom>` command.  When that
-command is used specific names are given to each attribute which are
-what is specified as the "name" portion of *i_name* or *d_name*.
+In addition, the various per-atom quantities listed above for specific
+packages are only accessible by this command.
 
 Output info
 """""""""""

doc/src/compute_temp_deform.rst

@@ -1,8 +1,11 @@
 .. index:: compute temp/deform
+.. index:: compute temp/deform/kk
 
 compute temp/deform command
 ===========================
 
+Accelerator Variants: *temp/deform/kk*
+
 Syntax
 """"""
 

doc/src/dihedral_fourier.rst

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

doc/src/dihedral_hybrid.rst

@@ -57,7 +57,7 @@ If *class2* is one of the dihedral hybrid styles, the same rule holds
 for specifying additional AngleTorsion (and EndBondTorsion, etc)
 coefficients either via the input script or in the data file.
 I.e. *class2* must be added to each line after the dihedral type.  For
-lines in the AngleTorsion (or EndBondTorsion, etc) section of the data
+lines in the AngleTorsion (or EndBondTorsion, etc) Coeffs section of the data
 file for dihedral types that are not *class2*, you must use an
 dihedral style of *skip* as a placeholder, e.g.
 

doc/src/dihedral_nharmonic.rst

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

doc/src/dihedral_quadratic.rst

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

doc/src/dihedral_spherical.rst

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

doc/src/dihedral_table.rst

@@ -244,9 +244,10 @@ script after reading the restart file.
 Restrictions
 """"""""""""
 
-These dihedral styles can only be used if LAMMPS was built with the
-MOLECULE package.  See the :doc:`Build package <Build_package>` doc
-page for more info.
+The *table* dihedral style can only be used if LAMMPS was built with the
+MOLECULE package.  The *table/cut* dihedral style can only be used if
+LAMMPS was built with the EXTRA-MOLECULE package. See the
+:doc:`Build package <Build_package>` doc page for more info.
 
 Related commands
 """"""""""""""""

doc/src/dump.rst

@@ -80,7 +80,8 @@ Syntax
                                q, mux, muy, muz, mu,
                                radius, diameter, omegax, omegay, omegaz,
                                angmomx, angmomy, angmomz, tqx, tqy, tqz,
-                               c_ID, c_ID[N], f_ID, f_ID[N], v_name
+                               c_ID, c_ID[I], f_ID, f_ID[I], v_name,
+                               i_name, d_name, i2_name[I], d2_name[I]
 
   .. parsed-literal::
 
@@ -110,8 +111,10 @@ Syntax
            f_ID = per-atom vector calculated by a fix with ID
            f_ID[I] = Ith column of per-atom array calculated by a fix with ID, I can include wildcard (see below)
            v_name = per-atom vector calculated by an atom-style variable with name
-           d_name = per-atom floating point vector with name, managed by fix property/atom
-           i_name = per-atom integer vector with name, managed by fix property/atom
+           i_name = custom integer vector with name
+           d_name = custom floating point vector with name
+           i2_name[I] = Ith column of custom integer array with name, I can include wildcard (see below)
+           d2_name[I] = Ith column of custom floating point vector with name, I can include wildcard (see below)
 
 * *local* or *local/gz* or *local/zstd* args = list of local attributes
 
@@ -474,16 +477,15 @@ styles.
 ----------
 
 Note that in the discussion which follows, for styles which can
-reference values from a compute or fix, like the *custom*, *cfg*, or
-*local* styles, the bracketed index I can be specified using a
-wildcard asterisk with the index to effectively specify multiple
-values.  This takes the form "\*" or "\*n" or "n\*" or "m\*n".  If N = the
-size of the vector (for *mode* = scalar) or the number of columns in
-the array (for *mode* = vector), then an asterisk with no numeric
-values means all indices from 1 to N.  A leading asterisk means all
-indices from 1 to n (inclusive).  A trailing asterisk means all
-indices from n to N (inclusive).  A middle asterisk means all indices
-from m to n (inclusive).
+reference values from a compute or fix or custom atom property, like
+the *custom*\ , *cfg*\ , or *local* styles, the bracketed index I can
+be specified using a wildcard asterisk with the index to effectively
+specify multiple values.  This takes the form "\*" or "\*n" or "n\*"
+or "m\*n".  If N = the number of columns in the array, then an
+asterisk with no numeric values means all column indices from 1 to N.
+A leading asterisk means all indices from 1 to n (inclusive).  A
+trailing asterisk means all indices from n to N (inclusive).  A middle
+asterisk means all indices from m to n (inclusive).
 
 Using a wildcard is the same as if the individual columns of the array
 had been listed one by one.  E.g. these 2 dump commands are
@@ -521,8 +523,9 @@ bonds and angles.
 
 Note that computes which calculate global or per-atom quantities, as
 opposed to local quantities, cannot be output in a dump local command.
-Instead, global quantities can be output by the :doc:`thermo_style custom <thermo_style>` command, and per-atom quantities can be
-output by the dump custom command.
+Instead, global quantities can be output by the :doc:`thermo_style
+custom <thermo_style>` command, and per-atom quantities can be output
+by the dump custom command.
 
 If *c_ID* is used as a attribute, then the local vector calculated by
 the compute is printed.  If *c_ID[I]* is used, then I must be in the
@@ -566,10 +569,11 @@ Nprocs-1) that currently owns the atom.  *Procp1* is the proc ID+1,
 which can be convenient in place of a *type* attribute (1 to Ntypes)
 for coloring atoms in a visualization program.  *Type* is the atom
 type (1 to Ntypes).  *Element* is typically the chemical name of an
-element, which you must assign to each type via the :doc:`dump_modify element <dump_modify>` command.  More generally, it can be any
-string you wish to associated with an atom type.  *Mass* is the atom
-mass.  *Vx*, *vy*, *vz*, *fx*, *fy*, *fz*, and *q* are components of
-atom velocity and force and atomic charge.
+element, which you must assign to each type via the :doc:`dump_modify
+element <dump_modify>` command.  More generally, it can be any string
+you wish to associated with an atom type.  *Mass* is the atom mass.
+*Vx*, *vy*, *vz*, *fx*, *fy*, *fz*, and *q* are components of atom
+velocity and force and atomic charge.
 
 There are several options for outputting atom coordinates.  The *x*,
 *y*, *z* attributes write atom coordinates "unscaled", in the
@@ -643,11 +647,12 @@ above for how I can be specified with a wildcard asterisk to
 effectively specify multiple values.
 
 The *f_ID* and *f_ID[I]* attributes allow vector or array per-atom
-quantities calculated by a :doc:`fix <fix>` to be output.  The ID in the
-attribute should be replaced by the actual ID of the fix that has been
-defined previously in the input script.  The :doc:`fix ave/atom <fix_ave_atom>` command is one that calculates per-atom
-quantities.  Since it can time-average per-atom quantities produced by
-any :doc:`compute <compute>`, :doc:`fix <fix>`, or atom-style
+quantities calculated by a :doc:`fix <fix>` to be output.  The ID in
+the attribute should be replaced by the actual ID of the fix that has
+been defined previously in the input script.  The :doc:`fix ave/atom
+<fix_ave_atom>` command is one that calculates per-atom quantities.
+Since it can time-average per-atom quantities produced by any
+:doc:`compute <compute>`, :doc:`fix <fix>`, or atom-style
 :doc:`variable <variable>`, this allows those time-averaged results to
 be written to a dump file.
 
@@ -664,14 +669,21 @@ should be replaced by the actual name of the variable that has been
 defined previously in the input script.  Only an atom-style variable
 can be referenced, since it is the only style that generates per-atom
 values.  Variables of style *atom* can reference individual atom
-attributes, per-atom attributes, thermodynamic keywords, or
-invoke other computes, fixes, or variables when they are evaluated, so
-this is a very general means of creating quantities to output to a
-dump file.
+attributes, per-atom attributes, thermodynamic keywords, or invoke
+other computes, fixes, or variables when they are evaluated, so this
+is a very general means of creating quantities to output to a dump
+file.
 
-The *d_name* and *i_name* attributes allow to output custom per atom
-floating point or integer properties that are managed by
-:doc:`fix property/atom <fix_property_atom>`.
+The *i_name*, *d_name*, *i2_name*, *d2_name* attributes refer to
+per-atom integer and floating-point vectors or arrays that have been
+added via the :doc:`fix property/atom <fix_property_atom>` command.
+When that command is used specific names are given to each attribute
+which are the "name" portion of these keywords.  For arrays *i2_name*
+and *d2_name*, the column of the array must also be included following
+the name in brackets: e.g. d2_xyz[I], i2_mySpin[I], where I is in the
+range from 1-M, where M is the number of columns in the custom array.
+See the discussion above for how I can be specified with a wildcard
+asterisk to effectively specify multiple values.
 
 See the :doc:`Modify <Modify>` page for information on how to add
 new compute and fix styles to LAMMPS to calculate per-atom quantities

doc/src/fix.rst

@@ -378,7 +378,8 @@ accelerated styles exist.
 * :doc:`thermal/conductivity <fix_thermal_conductivity>` - Muller-Plathe kinetic energy exchange for thermal conductivity calculation
 * :doc:`ti/spring <fix_ti_spring>` -
 * :doc:`tmd <fix_tmd>` - guide a group of atoms to a new configuration
-* :doc:`ttm <fix_ttm>` - two-temperature model for electronic/atomic coupling
+* :doc:`ttm <fix_ttm>` - two-temperature model for electronic/atomic coupling (replicated grid)
+* :doc:`ttm/grid <fix_ttm>` - two-temperature model for electronic/atomic coupling (distributed grid)
 * :doc:`ttm/mod <fix_ttm>` - enhanced two-temperature model with additional options
 * :doc:`tune/kspace <fix_tune_kspace>` - auto-tune KSpace parameters
 * :doc:`vector <fix_vector>` - accumulate a global vector every N timesteps

doc/src/fix_bond_swap.rst

@@ -27,21 +27,26 @@ Examples
 Description
 """""""""""
 
-In a simulation of polymer chains, this command attempts to swap bonds
-between two different chains, effectively grafting the end of one
-chain onto another chain and vice versa.  This is done via Monte Carlo
-rules using the Boltzmann acceptance criterion.  The purpose is to
-equilibrate the polymer chain conformations more rapidly than dynamics
-alone would do it, by enabling instantaneous large conformational
-changes in a dense polymer melt.  The polymer chains should thus more
-rapidly converge to the proper end-to-end distances and radii of
-gyration.  It is designed for use with systems of
-:doc:`FENE <bond_fene>` or :doc:`harmonic <bond_harmonic>` bead-spring
-polymer chains where each polymer is a linear chain of monomers, but
-LAMMPS does not enforce this requirement, i.e. any
-:doc:`bond_style <bond_style>` can be used.
+In a simulation of polymer chains this command attempts to swap a pair
+of bonds, as illustrated below.  This is done via Monte Carlo rules
+using the Boltzmann acceptance criterion, typically with the goal of
+equilibrating the polymer system more quickly.  This fix is designed
+for use with idealized bead-spring polymer chains where each polymer
+is a linear chain of monomers, but LAMMPS does not check that is the
+case for your system.
 
-A schematic of the kinds of bond swaps that can occur is shown here:
+Here are two use cases for this fix.
+
+The first use case is for swapping bonds on two different chains,
+effectively grafting the end of one chain onto the other chain and
+vice versa.  The purpose is to equilibrate the polymer chain
+conformations more rapidly than dynamics alone would do it, by
+enabling instantaneous large conformational changes in a dense polymer
+melt.  The polymer chains should thus more rapidly converge to the
+proper end-to-end distances and radii of gyration.
+
+A schematic of the kinds of bond swaps that can occur in this use case
+is shown here:
 
 .. image:: JPG/bondswap.jpg
    :align: center
@@ -53,38 +58,76 @@ attempt to delete the A1-A2 and B1-B2 bonds and replace them with
 A1-B2 and B1-A2 bonds.  If the swap is energetically favorable, the
 two chains on the right are the result and each polymer chain has
 undergone a dramatic conformational change.  This reference,
-:ref:`(Sides) <Sides>` provides more details on how the algorithm works and
-its application:
+:ref:`(Sides) <Sides>` provides more details on the algorithm's
+effectiveness for this use case.
 
-The bond swapping operation is invoked every *Nevery* timesteps.  If
-any bond is swapped, a re-build of the neighbor lists is triggered,
-since a swap alters the list of which neighbors are considered for
-pairwise interaction.  At each invocation, each processor considers a
-random specified *fraction* of its atoms as potential swapping
-monomers for this timestep.  Choosing a small *fraction* value can
-reduce the likelihood of a reverse swap occurring soon after an
-initial swap.
+The second use case is a collection of polymer chains with some
+fraction of their sites identified as "sticker" sites.  Initially each
+polymer chain is isolated from the others in a topological sense, and
+there is an intra-chain bond between every pair of sticker sites on
+the same chain.  Over time, bonds swap so that inter-molecular sticker
+bonds are created.  This models a vitrification-style process whereby
+the polymer chains all become interconnected.  For this use case, if
+angles are defined they should not include bonds between sticker
+sites.
 
-For each monomer A1, its neighbors are examined to find a possible B1
-monomer.  Both A1 and B1 must be in the fix group, their separation
-must be less than the specified *cutoff*, and the molecule IDs of A1
-and B1 must be the same (see below).  If a suitable partner is found,
-the energy change due to swapping the 2 bonds is computed.  This
-includes changes in pairwise, bond, and angle energies due to the
-altered connectivity of the 2 chains.  Dihedral and improper
-interactions are not allowed to be defined when this fix is used.
+----------
+
+The bond swapping operation is invoked once every *Nevery* timesteps.
+If any bond in the entire system is swapped, a re-build of the
+neighbor lists is triggered, since a swap alters the list of which
+neighbors are considered for pairwise interaction.  At each
+invocation, each processor considers a random specified *fraction* of
+its atoms as potential swapping monomers for this timestep.  Choosing
+a small *fraction* value can reduce the likelihood of a reverse swap
+occurring soon after an initial swap.
+
+For each monomer A1, its neighbors are looped over as B1 monomers.
+For each A1,B1 an additional double loop of bond partners A2 of A1,
+and bond partners B2 of B1 a is performed.  For each pair of A1-A2 and
+B1-B2 bonds to be eligible for swapping, the following 4 criteria must
+be met:
+
+(1) All 4 monomers must be in the fix group.
+
+(2) All 4 monomers must be owned by the processor (not ghost atoms).
+This insures that another processor does not attempt to swap bonds
+involving the same atoms on the same timestep.  Note that this also
+means that bond pairs which straddle processor boundaries are not
+eligible for swapping on this step.
+
+(3) The distances between 4 pairs of atoms -- (A1,A2), (B1,B2),
+(A1,B2), (B1,A2) -- must all be less than the specified *cutoff*\ .
+
+(4) The molecule IDs of A1 and B1 must be the same (see below).
+
+If an eligible B1 partner is found, the energy change due to swapping
+the 2 bonds is computed.  This includes changes in pairwise, bond, and
+angle energies due to the altered connectivity of the 2 chains.
+Dihedral and improper interactions are not allowed to be defined when
+this fix is used.
 
 If the energy decreases due to the swap operation, the bond swap is
 accepted.  If the energy increases it is accepted with probability
 exp(-delta/kT) where delta is the increase in energy, k is the
 Boltzmann constant, and T is the current temperature of the system.
-Whether the swap is accepted or rejected, no other swaps are attempted
-by this processor on this timestep.
 
-The criterion for matching molecule IDs is how bond swaps performed by
-this fix conserve chain length.  To use this features you must setup
-the molecule IDs for your polymer chains in a certain way, typically
-in the data file, read by the :doc:`read_data <read_data>` command.
+.. note::
+
+   IMPORTANT: Whether the swap is accepted or rejected, no other swaps
+   are attempted by this processor on this timestep.  No other
+   eligible 4-tuples of atoms are considered.  This means that each
+   processor will perform either a single swap or none on timesteps
+   this fix is invoked.
+
+----------
+
+The criterion for matching molecule IDs is how the first use case
+described above can be simulated while conserving chain lengths.  This
+is done by setting up the molecule IDs for the polymer chains in a
+specific way, typically in the data file, read by the :doc:`read_data
+<read_data>` command.
+
 Consider a system of 6-mer chains.  You have 2 choices.  If the
 molecule IDs for monomers on each chain are set to 1,2,3,4,5,6 then
 swaps will conserve chain length.  For a particular monomer there will
@@ -113,6 +156,9 @@ ends of a chain swap with each other.
    running dynamics, but can affect calculation of some diagnostic
    quantities or the printing of unwrapped coordinates to a dump file.
 
+For the second use case described above, the molecule IDs for all
+sticker sites should be the same.
+
 ----------
 
 This fix computes a temperature each time it is invoked for use by the
@@ -129,27 +175,28 @@ appended and the group for the new compute is "all", so that the
 temperature of the entire system is used.
 
 Note that this is NOT the compute used by thermodynamic output (see
-the :doc:`thermo_style <thermo_style>` command) with ID = *thermo_temp*.
-This means you can change the attributes of this fix's temperature
-(e.g. its degrees-of-freedom) via the
-:doc:`compute_modify <compute_modify>` command or print this temperature
-during thermodynamic output via the :doc:`thermo_style custom <thermo_style>` command using the appropriate compute-ID.
-It also means that changing attributes of *thermo_temp* will have no
-effect on this fix.
+the :doc:`thermo_style <thermo_style>` command) with ID =
+*thermo_temp*.  This means you can change the attributes of this fix's
+temperature (e.g. its degrees-of-freedom) via the :doc:`compute_modify
+<compute_modify>` command or print this temperature during
+thermodynamic output via the :doc:`thermo_style custom <thermo_style>`
+command using the appropriate compute-ID.  It also means that changing
+attributes of *thermo_temp* will have no effect on this fix.
 
 ----------
 
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
-No information about this fix is written to :doc:`binary restart files <restart>`.  Because the state of the random number generator
-is not saved in restart files, this means you cannot do "exact"
-restarts with this fix, where the simulation continues on the same as
-if no restart had taken place.  However, in a statistical sense, a
-restarted simulation should produce the same behavior.  Also note that
-each processor generates possible swaps independently of other
-processors.  Thus if you repeat the same simulation on a different number
-of processors, the specific swaps performed will be different.
+No information about this fix is written to :doc:`binary restart files
+<restart>`.  Because the state of the random number generator is not
+saved in restart files, this means you cannot do "exact" restarts with
+this fix, where the simulation continues on the same as if no restart
+had taken place.  However, in a statistical sense, a restarted
+simulation should produce the same behavior.  Also note that each
+processor generates possible swaps independently of other processors.
+Thus if you repeat the same simulation on a different number of
+processors, the specific swaps performed will be different.
 
 The :doc:`fix_modify <fix_modify>` *temp* option is supported by this
 fix.  You can use it to assign a :doc:`compute <compute>` you have
@@ -157,16 +204,18 @@ defined to this fix which will be used to compute the temperature for
 the Boltzmann criterion.
 
 This fix computes two statistical quantities as a global 2-vector of
-output, which can be accessed by various :doc:`output commands <Howto_output>`.  The first component of the vector is the
-cumulative number of swaps performed by all processors.  The second
-component of the vector is the cumulative number of swaps attempted
-(whether accepted or rejected).  Note that a swap "attempt" only
-occurs when swap partners meeting the criteria described above are
-found on a particular timestep.  The vector values calculated by this
-fix are "intensive".
+output, which can be accessed by various :doc:`output commands
+<Howto_output>`.  The first component of the vector is the cumulative
+number of swaps performed by all processors.  The second component of
+the vector is the cumulative number of swaps attempted (whether
+accepted or rejected).  Note that a swap "attempt" only occurs when
+swap partners meeting the criteria described above are found on a
+particular timestep.  The vector values calculated by this fix are
+"intensive".
 
 No parameter of this fix can be used with the *start/stop* keywords of
-the :doc:`run <run>` command.  This fix is not invoked during :doc:`energy minimization <minimize>`.
+the :doc:`run <run>` command.  This fix is not invoked during
+:doc:`energy minimization <minimize>`.
 
 Restrictions
 """"""""""""

doc/src/fix_halt.rst

@@ -19,7 +19,7 @@ Syntax
 
        bondmax = length of longest bond in the system (in length units)
        tlimit = elapsed CPU time (in seconds)
-       diskfree = free disk space (in megabytes)
+       diskfree = free disk space (in MBytes)
        v_name = name of :doc:`equal-style variable <variable>`
 
 * operator = "<" or "<=" or ">" or ">=" or "==" or "!=" or "\|\^"
@@ -81,7 +81,7 @@ the timer frequently across a large number of processors may be
 non-negligible.
 
 The *diskfree* attribute will check for available disk space (in
-megabytes) on supported operating systems. By default it will
+MBytes) on supported operating systems. By default it will
 check the file system of the current working directory.  This
 can be changed with the optional *path* keyword, which will take
 the path to a file or folder on the file system to be checked

doc/src/fix_nvt_sllod.rst

@@ -1,11 +1,12 @@
 .. index:: fix nvt/sllod
 .. index:: fix nvt/sllod/intel
 .. index:: fix nvt/sllod/omp
+.. index:: fix nvt/sllod/kk
 
 fix nvt/sllod command
 =====================
 
-Accelerator Variants: *nvt/sllod/intel*, *nvt/sllod/omp*
+Accelerator Variants: *nvt/sllod/intel*, *nvt/sllod/omp*, *nvt/sllod/kk*
 
 Syntax
 """"""

doc/src/fix_property_atom.rst

@@ -11,11 +11,11 @@ Syntax
 
 .. parsed-literal::
 
-   fix ID group-ID property/atom vec1 vec2 ... keyword value ...
+   fix ID group-ID property/atom name1 name2 ... keyword value ...
 
 * ID, group-ID are documented in :doc:`fix <fix>` command
 * property/atom = style name of this fix command
-* vec1,vec2,... = *mol* or *q* or *rmass* or *i_name* or *d_name*
+* name1,name2,... = *mol* or *q* or *rmass* or *i_name* or *d_name* or *i2_name* or *d2_name*
 
   .. parsed-literal::
 
@@ -24,6 +24,10 @@ Syntax
        *rmass* = per-atom mass
        *i_name* = new integer vector referenced by name
        *d_name* = new floating-point vector referenced by name
+       *i2_name* = new integer array referenced by name
+          i2_name arg = N = number of columns in the array
+       *d2_name* = new floating-point array referenced by name
+          d2_name arg = N = number of columns in the array
 
 * zero of more keyword/value pairs may be appended
 * keyword = *ghost*
@@ -39,58 +43,64 @@ Examples
 
    fix 1 all property/atom mol
    fix 1 all property/atom i_myflag1 i_myflag2
-   fix 1 all property/atom d_sx d_sy d_sz
+   fix 1 all property/atom d2_sxyz 3 ghost yes
 
 Description
 """""""""""
 
-Create one or more additional per-atom vectors to store information
-about atoms and to use during a simulation.  The specified *group-ID*
-is ignored by this fix.
+Create one or more additional per-atom vectors or arrays to store
+information about atoms and to use during a simulation.  The specified
+*group-ID* is ignored by this fix.
 
 The atom style used for a simulation defines a set of per-atom
 properties, as explained on the :doc:`atom_style <atom_style>` and
-:doc:`read_data <read_data>` doc pages.  The latter command allows these
-properties to be defined for each atom in the system when a data file
-is read.  This fix will augment the set of properties with new custom
+:doc:`read_data <read_data>` doc pages.  The latter command defines
+these properties for each atom in the system when a data file is read.
+This fix augments the set of per-atom properties with new custom
 ones. This can be useful in several scenarios.
 
-If the atom style does not define molecule IDs, per-atom charge,
-or per-atom mass, they can be added using the *mol*, *q* or *rmass*
-keywords.  This can be useful, e.g, to define "molecules" to use as
-rigid bodies with the :doc:`fix rigid <fix_rigid>` command, or just to
-carry around an extra flag with the atoms (stored as a molecule ID)
-that can be used to group atoms without having to use the group
+If the atom style does not define molecule IDs, per-atom charge, or
+per-atom mass, they can be added using the *mol*\ , *q* or *rmass*
+keywords.  This could be useful to define "molecules" to use as rigid
+bodies with the :doc:`fix rigid <fix_rigid>` command, or to carry
+around an extra flag with atoms (stored as a molecule ID) that can be
+used by various commands like :doc:`compute chunk/atom
+<compute_chunk_atom>` to group atoms without having to use the group
 command (which is limited to a total of 32 groups including *all*\ ).
 
-Another application would be to use the *rmass* flag in order to have
-per-atom masses instead of per-type masses, for example this can be
-useful to study isotope effects with partial isotope substitution.
-Please :ref:`see below <isotopes>` for an example of simulating a mixture
-of light and heavy water with the TIP4P water potential.
+Another application is to use the *rmass* flag in order to have
+per-atom masses instead of per-type masses.  This could be used to
+study isotope effects with partial isotope substitution. :ref:`See
+below <isotopes>` for an example of simulating a mixture of light and
+heavy water with the TIP4P water potential.
 
-An alternative to using fix *property/atom* in these ways is to
+An alternative to using fix *property/atom* for these examples is to
 use an atom style that does define molecule IDs or charge or per-atom
 mass (indirectly via diameter and density) or to use a hybrid atom
-style that combines two or more atom styles
-to provide the union of all atom properties. However, this has two
-practical drawbacks:  first, it typically necessitates changing the
-format of the data file, which can be tedious for large systems;
-and second, it may define additional properties that are not needed
-such as bond lists, which has some overhead when there are no bonds.
+style that combines two or more atom styles to provide the union of
+all their atom properties. However, this has two practical drawbacks:
+first, it typically necessitates changing the format of the Atoms
+section in the data file and second, it may define additional
+properties that are not needed such as bond lists, which incurs some
+overhead when there are no bonds.
 
-In the future, we may add additional per-atom properties similar to
-*mol*, *q* or *rmass*, which "turn-on" specific properties defined
-by some atom styles, so they can be used by atom styles that do not
-define them.
+In the future, we may add additional existing per-atom properties to
+fix property/atom, similar to *mol*\ , *q* or *rmass*\ , which
+"turn-on" specific properties defined by some atom styles, so they can
+be easily used by atom styles that do not define them.
 
-More generally, the *i_name* and *d_name* vectors allow one or more
-new custom per-atom properties to be defined.  Each name must be
-unique and can use alphanumeric or underscore characters.  These
-vectors can store whatever values you decide are useful in your
-simulation.  As explained below there are several ways to initialize
-and access and output these values, both via input script commands and
-in new code that you add to LAMMPS.
+More generally, the *i_name* and *d_name* options allow one or more
+new custom per-atom vectors to be defined.  Likewise the *i2_name* and
+*d2_name* options allow one or more custom per-atom arrays to be
+defined.  The *i2_name* and *d2_name* options take an argument *N*
+which specifies the number of columns in the per-atom array, i.e. the
+number of attributes associated with each atom.  *N* >= 1 is required.
+
+Each name must be unique and can use alphanumeric or underscore
+characters.  These vectors and arrays can store whatever values you
+decide are useful in your simulation.  As explained below there are
+several ways to initialize, access, and output these values, via input
+script commands, data files, and in new code you add to LAMMPS.
 
 This is effectively a simple way to add per-atom properties to a model
 without needing to write code for a new :doc:`atom style <atom_style>`
@@ -108,43 +118,39 @@ new properties are also defined for the ghost atoms.
 .. admonition:: Properties on ghost atoms
    :class: note
 
-   If you use this command with the *mol*, *q* or *rmass* vectors,
-   then you most likely want to set *ghost* yes, since these properties
-   are stored with ghost atoms if you use an :doc:`atom_style <atom_style>`
-   that defines them, and many LAMMPS operations that use molecule IDs or
-   charge, such as neighbor lists and pair styles, will expect ghost
-   atoms to have these values.  LAMMPS will issue a warning it you define
-   those vectors but do not set *ghost* yes.
+   If you use the *mol*\ , *q* or *rmass* names, you most likely want
+   to set *ghost* yes, since these properties are stored with ghost
+   atoms if you use an :doc:`atom_style <atom_style>` that defines
+   them.  Many LAMMPS operations that use molecule IDs or charge, such
+   as neighbor lists and pair styles, will expect ghost atoms to have
+   these values.  LAMMPS will issue a warning it you define those
+   vectors but do not set *ghost* yes.
 
 .. admonition:: Limitations on ghost atom properties
    :class: note
 
-   The properties for ghost atoms are not updated every timestep,
-   but only once every few steps when neighbor lists are re-built.  Thus
-   the *ghost* keyword is suitable for static properties, like molecule
-   IDs, but not for dynamic properties that change every step.  For the
-   latter, the code you add to LAMMPS to change the properties will also
-   need to communicate their new values to/from ghost atoms, an operation
-   that can be invoked from within a :doc:`pair style <pair_style>` or
-   :doc:`fix <fix>` or :doc:`compute <compute>` that you write.
+   The specified properties for ghost atoms are not updated every
+   timestep, but only once every few steps when neighbor lists are
+   re-built.  Thus the *ghost* keyword is suitable for static
+   properties, like molecule IDs, but not for dynamic properties that
+   change every step.  For the latter, the code you add to LAMMPS to
+   change the properties will also need to communicate their new
+   values to/from ghost atoms, an operation that can be invoked from
+   within a :doc:`pair style <pair_style>` or :doc:`fix <fix>` or
+   :doc:`compute <compute>` that you write.
 
 ----------
 
 This fix is one of a small number that can be defined in an input
 script before the simulation box is created or atoms are defined.
-This is so it can be used with the :doc:`read_data <read_data>` command
-as described below.
+This is so it can be used with the :doc:`read_data <read_data>`
+command as described next.
 
-.. note::
-
-   If this fix is defined **after** the simulation box is created,
-   a 'run 0' command may be needed to properly initialize the storage
-   created by this fix.
-
-Per-atom properties that are defined by the :doc:`atom style <atom_style>` are initialized when atoms are created, e.g. by
-the :doc:`read_data <read_data>` or :doc:`create_atoms <create_atoms>`
+Per-atom properties that are defined by the :doc:`atom style
+<atom_style>` are initialized when atoms are created, e.g. by the
+:doc:`read_data <read_data>` or :doc:`create_atoms <create_atoms>`
 commands.  The per-atom properties defined by this fix are not.  So
-you need to initialize them explicitly.  This can be done by the
+you need to initialize them explicitly.  One way to do this is
 :doc:`read_data <read_data>` command, using its *fix* keyword and
 passing it the fix-ID of this fix.
 
@@ -169,15 +175,24 @@ would allow a data file to have a section like this:
    ...
    N 763 4.5
 
-where N is the number of atoms, and the first field on each line is
-the atom-ID, followed by a molecule-ID and a floating point value that
-will be stored in a new property called "flag".  Note that the list of
-per-atom properties can be in any order.
+where N is the number of atoms, the first field on each line is the
+atom-ID, the next two are a molecule-ID and a floating point value
+that will be stored in a new property called "flag".  If a per-atom
+array was specified in the fix property/atom command then the *N*
+values for that array must be specified consecutively for that
+property on each line.  Note that the order of values on each line
+corresponds to the order of custom names in the fix property/atom
+command.
 
-Another way of initializing the new properties is via the
-:doc:`set <set>` command.  For example, if you wanted molecules
-defined for every set of 10 atoms, based on their atom-IDs,
-these commands could be used:
+Note that the the lines of per-atom properties can be listed in any
+order.  Also note that all the per-atom properties specified by the
+fix ID (prop in this case) must be included on each line in the
+specified data file section (Molecules in this case).
+
+Another way of initializing the new properties is via the :doc:`set
+<set>` command.  For example, if you wanted molecules defined for
+every set of 10 atoms, based on their atom-IDs, these commands could
+be used:
 
 .. code-block:: LAMMPS
 
@@ -187,53 +202,59 @@ these commands could be used:
 
 The :doc:`atom-style variable <variable>` will create values for atoms
 with IDs 31,32,33,...40 that are 4.0,4.1,4.2,...,4.9.  When the
-:doc:`set <set>` commands assigns them to the molecule ID for each atom,
-they will be truncated to an integer value, so atoms 31-40 will all be
-assigned a molecule ID of 4.
+:doc:`set <set>` commands assigns them to the molecule ID for each
+atom, they will be truncated to an integer value, so atoms 31-40 will
+all be assigned a molecule ID of 4.
 
-Note that :doc:`atomfile-style variables <variable>` can also be used in
-place of atom-style variables, which means in this case that the
+Note that :doc:`atomfile-style variables <variable>` can also be used
+in place of atom-style variables, which means in this case that the
 molecule IDs could be read-in from a separate file and assigned by the
 :doc:`set <set>` command.  This allows you to initialize new per-atom
 properties in a completely general fashion.
 
 ----------
 
-For new atom properties specified as *i_name* or *d_name*, the
-:doc:`compute property/atom <compute_property_atom>` command can access
-their values.  This means that the values can be output via the :doc:`dump custom <dump>` command, accessed by fixes like :doc:`fix ave/atom <fix_ave_atom>`, accessed by other computes like :doc:`compute reduce <compute_reduce>`, or used in :doc:`atom-style variables <variable>`.
+For new atom properties specified as *i_name*, *d_name*, *i2_name*, or
+*d2_name*, the :doc:`dump custom <dump>` and :doc:`compute
+property/atom <compute_property_atom>` commands can access their
+values.  This means that the values can be used accessed by fixes like
+:doc:`fix ave/atom <fix_ave_atom>`, accessed by other computes like
+:doc:`compute reduce <compute_reduce>`, or used in :doc:`atom-style
+variables <variable>`.
 
-For example, these commands will output two new properties to a custom
-dump file:
+For example, these commands will output both the instantaneous and
+time-averaged values of two new properties to a custom dump file:
 
 .. code-block:: LAMMPS
 
-   fix prop all property/atom i_flag1 d_flag2
+   fix myprops all property/atom i_flag1 d_flag2
    compute 1 all property/atom i_flag1 d_flag2
-   dump 1 all custom 100 tmp.dump id x y z c_1[1] c_1[2]
+   fix 1 all ave/atom 10 10 100 c_1[1] c_1[2]
+   dump 1 all custom 100 tmp.dump id x y z i_flag1 d_flag2 f_1[1] f_1[2]
 
 ----------
 
-If you wish to add new :doc:`pair styles <pair_style>`,
-:doc:`fixes <fix>`, or :doc:`computes <compute>` that use the per-atom
-properties defined by this fix, see the :doc:`Modify atom <Modify_atom>`
-doc page which has details on how the properties can be accessed from
-added classes.
+If you wish to add new :doc:`pair styles <pair_style>`, :doc:`fixes
+<fix>`, or :doc:`computes <compute>` that use the per-atom properties
+defined by this fix, see the :doc:`Modify atom <Modify_atom>` doc page
+which has details on how the custom properties of this fix can be
+accessed from added classes.
 
 ----------
 
 .. _isotopes:
 
-Example for using per-atom masses with TIP4P water to
-study isotope effects. When setting up simulations with the :doc:`TIP4P pair styles <Howto_tip4p>` for water, you have to provide exactly
-one atom type each to identify the water oxygen and hydrogen
-atoms. Since the atom mass is normally tied to the atom type, this
-makes it impossible to study multiple isotopes in the same simulation.
-With *fix property/atom rmass* however, the per-type masses are
-replaced by per-atom masses. Asumming you have a working input deck
-for regular TIP4P water, where water oxygen is atom type 1 and water
-hydrogen is atom type 2, the following lines of input script convert
-this to using per-atom masses:
+Here is an example of using per-atom masses with TIP4P water to study
+isotope effects. When setting up simulations with the :doc:`TIP4P pair
+styles <Howto_tip4p>` for water, you have to provide exactly one atom
+type each to identify the water oxygen and hydrogen atoms. Since the
+atom mass is normally tied to the atom type, this makes it impossible
+to study multiple isotopes in the same simulation.  With *fix
+property/atom rmass* however, the per-type masses are replaced by
+per-atom masses. Asumming you have a working input deck for regular
+TIP4P water, where water oxygen is atom type 1 and water hydrogen is
+atom type 2, the following lines of input script convert this to using
+per-atom masses:
 
 .. code-block:: LAMMPS
 
@@ -241,22 +262,22 @@ this to using per-atom masses:
    set type 1 mass 15.9994
    set type 2 mass 1.008
 
-When writing out the system data with the :doc:`write_data <write_data>`
-command, there will be a new section named with the fix-ID
-(i.e. *Isotopes* in this case). Alternatively, you can take an
-existing data file and just add this *Isotopes* section with
-one line per atom containing atom-ID and mass. Either way, the
-extended data file can be read back with:
+When writing out the system data with the :doc:`write_data
+<write_data>` command, there will be a new section named with the
+fix-ID (i.e. *Isotopes* in this case). Alternatively, you can take an
+existing data file and just add this *Isotopes* section with one line
+per atom containing atom-ID and mass. Either way, the extended data
+file can be read back with:
 
 .. code-block:: LAMMPS
 
    fix Isotopes all property/atom rmass ghost yes
    read_data tip4p-isotopes.data fix Isotopes NULL Isotopes
 
-Please note that the first *Isotopes* refers to the fix-ID
-and the second to the name of the section. The following input
-script code will now change the first 100 water molecules in this
-example to heavy water:
+Please note that the first *Isotopes* refers to the fix-ID and the
+second to the name of the section. The following input script code
+will now change the first 100 water molecules in this example to heavy
+water:
 
 .. code-block:: LAMMPS
 
@@ -276,24 +297,27 @@ Restart, fix_modify, output, run start/stop, minimize info
 This fix writes the per-atom values it stores to :doc:`binary restart
 files <restart>`, so that the values can be restored when a simulation
 is restarted.  See the :doc:`read_restart <read_restart>` command for
-info on how to re-specify a fix in an input script that reads a restart
-file, so that the operation of the fix continues in an uninterrupted
-fashion.
+info on how to re-specify a fix in an input script that reads a
+restart file, so that the operation of the fix continues in an
+uninterrupted fashion.
 
 .. warning::
 
-   When reading data from a restart, the fix command has to be specified
-   **exactly** the same way as before. LAMMPS will only check whether a
-   fix is of the same style and has the same fix ID and in case of a match
-   will then try to initialize the fix with the data stored in the binary
-   restart file.  If the fix property/atom command does not match exactly,
-   data can be corrupted or LAMMPS may crash.
+   When reading data from a restart file, this fix command has to be
+   specified **exactly** the same was in the input script that created
+   the restart file. LAMMPS will only check whether a fix is of the
+   same style and has the same fix ID and in case of a match will then
+   try to initialize the fix with the data stored in the binary
+   restart file.  If the names and associated date types in the new
+   fix property/atom command do not match the old one exactly, data
+   can be corrupted or LAMMPS may crash.
 
-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
+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>`.
+:doc:`run <run>` command.  This fix is not invoked during :doc:`energy
+minimization <minimize>`.
 
 Restrictions
 """"""""""""
@@ -302,9 +326,10 @@ Restrictions
 Related commands
 """"""""""""""""
 
-:doc:`read_data <read_data>`, :doc:`set <set>`, :doc:`compute property/atom <compute_property_atom>`
+:doc:`read_data <read_data>`, :doc:`set <set>`,
+:doc:`compute property/atom <compute_property_atom>`
 
 Default
 """""""
 
-The default keyword values are ghost = no.
+The default keyword value is ghost = no.

doc/src/fix_qtb.rst

@@ -128,7 +128,7 @@ spectrum while consumes more memory.  With fixed *f_max* and
 :math:`\gamma`, *N_f* should be big enough to converge the classical
 temperature :math:`T^{cl}` as a function of target quantum bath
 temperature. Memory usage per processor could be from 10 to 100
-Mbytes.
+MBytes.
 
 .. note::
 

doc/src/fix_rigid.rst

@@ -73,7 +73,7 @@ Syntax
        *single* args = none
        *molecule* args = none
        *custom* args = *i_propname* or *v_varname*
-         i_propname = an integer property defined via fix property/atom
+         i_propname = a custom integer vector defined via fix property/atom
          v_varname  = an atom-style or atomfile-style variable
        *group* args = N groupID1 groupID2 ...
          N = # of groups
@@ -296,15 +296,16 @@ includes atoms you want to be part of rigid bodies.
 
 Bodystyle *custom* is similar to bodystyle *molecule* except that it
 is more flexible in using other per-atom properties to define the sets
-of atoms that form rigid bodies.  An integer vector defined by the
-:doc:`fix property/atom <fix_property_atom>` command can be used.  Or an
-:doc:`atom-style or atomfile-style variable <variable>` can be used; the
-floating-point value produced by the variable is rounded to an
-integer.  As with bodystyle *molecule*, each set of atoms in the fix
-groups with the same integer value is treated as a different rigid
-body.  Since fix property/atom vectors and atom-style variables
-produce values for all atoms, you should be careful to use a fix group
-that only includes atoms you want to be part of rigid bodies.
+of atoms that form rigid bodies.  A custom per-atom integer vector
+defined by the :doc:`fix property/atom <fix_property_atom>` command
+can be used.  Or an :doc:`atom-style or atomfile-style variable
+<variable>` can be used; the floating-point value produced by the
+variable is rounded to an integer.  As with bodystyle *molecule*\ ,
+each set of atoms in the fix groups with the same integer value is
+treated as a different rigid body.  Since fix property/atom custom
+vectors and atom-style variables produce values for all atoms, you
+should be careful to use a fix group that only includes atoms you want
+to be part of rigid bodies.
 
 .. note::
 

doc/src/fix_store_state.rst

@@ -23,8 +23,8 @@ Syntax
                              q, mux, muy, muz, mu,
                              radius, diameter, omegax, omegay, omegaz,
                              angmomx, angmomy, angmomz, tqx, tqy, tqz,
-                             c_ID, c_ID[N], f_ID, f_ID[N], v_name,
-                             d_name, i_name
+                             c_ID, c_ID[I], f_ID, f_ID[I], v_name,
+                             d_name, i_name, i2_name[I], d2_name[I],
 
   .. parsed-literal::
 
@@ -46,13 +46,15 @@ Syntax
            omegax,omegay,omegaz = angular velocity of spherical particle
            angmomx,angmomy,angmomz = angular momentum of aspherical particle
            tqx,tqy,tqz = torque on finite-size particles
-           c_ID = per-atom vector calculated by a compute with ID
-           c_ID[I] = Ith column of per-atom array calculated by a compute with ID
-           f_ID = per-atom vector calculated by a fix with ID
-           f_ID[I] = Ith column of per-atom array calculated by a fix with ID
-           v_name = per-atom vector calculated by an atom-style variable with name
-           d_name = per-atom floating point vector name, managed by fix property/atom
-           i_name = per-atom integer vector name, managed by fix property/atom
+           *c_ID* = per-atom vector calculated by a compute with ID
+           *c_ID[I]* = Ith column of per-atom array calculated by a compute with ID
+           *f_ID* = per-atom vector calculated by a fix with ID
+           *f_ID[I]* = Ith column of per-atom array calculated by a fix with ID
+           *v_name* = per-atom vector calculated by an atom-style variable with name
+           *i_name* = custom integer vector with name
+           *d_name* = custom floating point vector with name
+           *i2_name[I]* = Ith column of custom integer array with name
+           *d2_name[I]* = Ith column of custom floating-point array with name
 
 * zero or more keyword/value pairs may be appended
 * keyword = *com*
@@ -92,7 +94,8 @@ steps.
    those attributes may require quantities that are not defined in
    between runs.
 
-The list of possible attributes is the same as that used by the :doc:`dump custom <dump>` command, which describes their meaning.
+The list of possible attributes is the same as that used by the
+:doc:`dump custom <dump>` command, which describes their meaning.
 
 If the *com* keyword is set to *yes* then the *xu*, *yu*, and *zu*
 inputs store the position of each atom relative to the center-of-mass
@@ -105,15 +108,16 @@ group.
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
-This fix writes the per-atom values it stores to :doc:`binary restart files <restart>`, so that the values can be restored when a
-simulation is restarted.  See the :doc:`read_restart <read_restart>`
-command for info on how to re-specify a fix in an input script that
-reads a restart file, so that the operation of the fix continues in an
+This fix writes the per-atom values it stores to :doc:`binary restart
+files <restart>`, so that the values can be restored when a simulation
+is restarted.  See the :doc:`read_restart <read_restart>` command for
+info on how to re-specify a fix in an input script that reads a
+restart file, so that the operation of the fix continues in an
 uninterrupted fashion.
 
 .. warning::
 
-   When reading data from a restart, the fix command has to be specified
+   When reading data from a restart file, this fix command has to be specified
    **exactly** the same way as before. LAMMPS will only check whether a
    fix is of the same style and has the same fix ID and in case of a match
    will then try to initialize the fix with the data stored in the binary
@@ -130,7 +134,8 @@ can be accessed by various :doc:`output commands <Howto_output>`.  The
 per-atom values be accessed on any timestep.
 
 No parameter of this fix can be used with the *start/stop* keywords of
-the :doc:`run <run>` command.  This fix is not invoked during :doc:`energy minimization <minimize>`.
+the :doc:`run <run>` command.  This fix is not invoked during
+:doc:`energy minimization <minimize>`.
 
 Restrictions
 """"""""""""

doc/src/fix_ttm.rst

@@ -1,9 +1,13 @@
 .. index:: fix ttm
+.. index:: fix ttm/grid
 .. index:: fix ttm/mod
 
 fix ttm command
 ===============
 
+fix ttm/grid command
+====================
+
 fix ttm/mod command
 ===================
 
@@ -12,13 +16,13 @@ Syntax
 
 .. parsed-literal::
 
-   fix ID group-ID ttm seed C_e rho_e kappa_e gamma_p gamma_s v_0 Nx Ny Nz T_infile N T_outfile
-   fix ID group-ID ttm/mod seed init_file Nx Ny Nz T_infile N T_outfile
+   fix ID group-ID ttm seed C_e rho_e kappa_e gamma_p gamma_s v_0 Nx Ny Nz keyword value ...
+   fix ID group-ID ttm/mod seed init_file Nx Ny Nz keyword value ...
 
 * ID, group-ID are documented in :doc:`fix <fix>` command
-* style = *ttm* or *ttm_mod*
+* style = *ttm* or *ttm/grid* or *ttm/mod*
 * seed = random number seed to use for white noise (positive integer)
-* remaining arguments for fix ttm:
+* remaining arguments for fix ttm or fix ttm/grid
 
   .. parsed-literal::
 
@@ -31,9 +35,6 @@ Syntax
        Nx = number of thermal solve grid points in the x-direction (positive integer)
        Ny = number of thermal solve grid points in the y-direction (positive integer)
        Nz = number of thermal solve grid points in the z-direction (positive integer)
-       T_infile = filename to read initial electronic temperature from
-       N = dump TTM temperatures every this many timesteps, 0 = no dump
-       T_outfile = filename to write TTM temperatures to (only needed if N > 0)
 
 * remaining arguments for fix ttm/mod:
 
@@ -43,18 +44,29 @@ Syntax
        Nx = number of thermal solve grid points in the x-direction (positive integer)
        Ny = number of thermal solve grid points in the y-direction (positive integer)
        Nz = number of thermal solve grid points in the z-direction (positive integer)
-       T_infile = filename to read initial electronic temperature from
-       N = dump TTM temperatures every this many timesteps, 0 = no dump
-       T_outfile = filename to write TTM temperatures to (only needed if N > 0)
+
+* zero or more keyword/value(s) pairs may be appended
+* keyword = *set* or *infile* or *outfile*
+
+  .. parsed-literal::
+
+       *set* value = Tinit
+         Tinit = initial electronic temperature at all grid points (temperature units)
+       *infile* value = file.in with grid values for electronic temperatures
+       *outfile* values = Nout file.out
+         Nout = dump grid temperatures every this many timesteps
+         file.out = filename to write grid temperatures to
 
 Examples
 """"""""
 
 .. code-block:: LAMMPS
 
-   fix 2 all ttm 699489 1.0 1.0 10 0.1 0.0 2.0 1 12 1 initialTs 1000 T.out
-   fix 2 all ttm 123456 1.0 1.0 1.0 1.0 1.0 5.0 5 5 5 Te.in 1 Te.out
-   fix 2 all ttm/mod 34277 parameters.txt 5 5 5 T_init 10 T_out
+   fix 2 all ttm 699489 1.0 1.0 10 0.1 0.0 2.0 1 12 1 infile initial outfile 1000 T.out
+   fix 3 all ttm/grid 123456 1.0 1.0 1.0 1.0 1.0 5.0 5 5 5 infile Te.in
+   fix 4 all ttm/mod 34277 parameters.txt 5 5 5 infile T_init outfile 10 T_out
+
+Example input scripts using these commands can be found in examples/ttm.
 
 Description
 """""""""""
@@ -62,36 +74,48 @@ Description
 Use a two-temperature model (TTM) to represent heat transfer through
 and between electronic and atomic subsystems.  LAMMPS models the
 atomic subsystem as usual with a molecular dynamics model and the
-classical force field specified by the user, but the electronic
-subsystem is modeled as a continuum, or a background "gas", on a
-regular grid.  Energy can be transferred spatially within the grid
-representing the electrons.  Energy can also be transferred between
-the electronic and the atomic subsystems.  The algorithm underlying
-this fix was derived by D. M.  Duffy and A. M. Rutherford and is
-discussed in two J Physics: Condensed Matter papers: :ref:`(Duffy) <Duffy>`
-and :ref:`(Rutherford) <Rutherford>`.  They used this algorithm in cascade
-simulations where a primary knock-on atom (PKA) was initialized with a
-high velocity to simulate a radiation event.
+classical force field specified by the user.  The electronic subsystem
+is modeled as a continuum, or a background "gas", on a regular grid
+which overlays the simulation domain.  Energy can be transferred
+spatially within the grid representing the electrons.  Energy can also
+be transferred between the electronic and atomic subsystems.  The
+algorithm underlying this fix was derived by D. M.  Duffy
+and A. M. Rutherford and is discussed in two J Physics: Condensed
+Matter papers: :ref:`(Duffy) <Duffy>` and :ref:`(Rutherford)
+<Rutherford>`.  They used this algorithm in cascade simulations where
+a primary knock-on atom (PKA) was initialized with a high velocity to
+simulate a radiation event.
 
-The description in this sub-section applies to both fix ttm and fix
-ttm/mod.  Fix ttm/mod adds options to account for external heat
-sources (e.g. at a surface) and for specifying parameters that allow
-the electronic heat capacity to depend strongly on electronic
-temperature.  It is more expensive computationally than fix ttm
-because it treats the thermal diffusion equation as non-linear.  More
-details on fix ttm/mod are given below.
+The description in this sub-section applies to all 3 fix styles:
+*ttm*, *ttm/grid*, and *ttm/mod*.
+
+Fix *ttm/grid* distributes the regular grid across processors consistent
+with the sub-domains of atoms owned by each processor, but is otherwise
+identical to fix ttm.  Note that fix *ttm* stores a copy of the grid on
+each processor, which is acceptable when the overall grid is reasonably
+small.  For larger grids you should use fix *ttm/grid* instead.
+
+Fix *ttm/mod* adds options to account for external heat sources (e.g. at
+a surface) and for specifying parameters that allow the electronic
+heat capacity to depend strongly on electronic temperature.  It is
+more expensive computationally than fix *ttm* because it treats the
+thermal diffusion equation as non-linear.  More details on fix *ttm/mod*
+are given below.
 
 Heat transfer between the electronic and atomic subsystems is carried
-out via an inhomogeneous Langevin thermostat.  This thermostat differs
-from the regular Langevin thermostat (:doc:`fix langevin <fix_langevin>`) in three important ways.  First, the
-Langevin thermostat is applied uniformly to all atoms in the
+out via an inhomogeneous Langevin thermostat.  Only atoms in the fix
+group contribute to and are affected by this heat transfer.
+
+This thermostatting differs from the regular Langevin thermostat
+(:doc:`fix langevin <fix_langevin>`) in three important ways.  First,
+the Langevin thermostat is applied uniformly to all atoms in the
 user-specified group for a single target temperature, whereas the TTM
-fix applies Langevin thermostatting locally to atoms within the
+fixes apply Langevin thermostatting locally to atoms within the
 volumes represented by the user-specified grid points with a target
 temperature specific to that grid point.  Second, the Langevin
 thermostat couples the temperature of the atoms to an infinite heat
-reservoir, whereas the heat reservoir for fix TTM is finite and
-represents the local electrons.  Third, the TTM fix allows users to
+reservoir, whereas the heat reservoir for the TTM fixes is finite and
+represents the local electrons.  Third, the TTM fixes allow users to
 specify not just one friction coefficient, but rather two independent
 friction coefficients: one for the electron-ion interactions
 (*gamma_p*), and one for electron stopping (*gamma_s*).
@@ -123,29 +147,59 @@ as that in equation 6 of :ref:`(Duffy) <Duffy>`, with the exception that the
 electronic density is explicitly represented, rather than being part
 of the specific heat parameter.
 
-Currently, fix ttm assumes that none of the user-supplied parameters
-will vary with temperature. Note that :ref:`(Duffy) <Duffy>` used a tanh()
-functional form for the temperature dependence of the electronic
-specific heat, but ignored temperature dependencies of any of the
-other parameters.  See more discussion below for fix ttm/mod.
+Currently, the TTM fixes assume that none of the user-supplied
+parameters will vary with temperature. Note that :ref:`(Duffy)
+<Duffy>` used a tanh() functional form for the temperature dependence
+of the electronic specific heat, but ignored temperature dependencies
+of any of the other parameters.  See more discussion below for fix
+ttm/mod.
 
-These fixes require use of periodic boundary conditions and a 3D
-simulation.  Periodic boundary conditions are also used in the heat
-equation solve for the electronic subsystem.  This varies from the
-approach of :ref:`(Rutherford) <Rutherford>` where the atomic subsystem was
+..note::
+
+  These fixes do not perform time integration of the atoms in the fix
+  group, they only rescale their velocities.  Thus a time integration
+  fix such as :doc:`fix nve <fix_nve>` should be used in conjunction
+  with these fixes.  These fixes should not normally be used on atoms
+  that have their temperature controlled by another thermostatting
+  fix, e.g. :doc:`fix nvt <fix_nh>` or :doc:`fix langevin
+  <fix_langevin>`.
+
+..note::
+
+  These fixes require use of an orthogonal 3d simulation box with
+  periodic boundary conditions in all dimensions.  They also require
+  that the size and shape of the simulation box do not vary
+  dynamically, e.g. due to use of the :doc:`fix npt <fix_nh>` command.
+  Likewise, the size/shape of processor sub-domains cannot vary due to
+  dynamic load-balancing via use of the :doc:`fix balance
+  <fix_balance>` command.  It is possible however to load balance
+  before the simulation starts using the :doc:`balance <balance>`
+  command, so that each processor has a different size sub-domain.
+
+Periodic boundary conditions are also used in the heat equation solve
+for the electronic subsystem.  This varies from the approach of
+:ref:`(Rutherford) <Rutherford>` where the atomic subsystem was
 embedded within a larger continuum representation of the electronic
 subsystem.
 
-The initial electronic temperature input file, *T_infile*, is a text
-file LAMMPS reads in with no header and with four numeric columns
-(ix,iy,iz,Temp) and with a number of rows equal to the number of
-user-specified grid points (Nx by Ny by Nz).  The ix,iy,iz are node
-indices from 0 to nxnodes-1, etc.  For example, the initial electronic
-temperatures on a 1 by 2 by 3 grid could be specified in a *T_infile*
-as follows:
+The *set* keyword specifies a *Tinit* temperature value to initialize
+the value stored on all grid points.
+
+The *infile* keyword specifies an input file of electronic temperatures
+for each grid point to be read in to initialize the grid.  By default
+the temperatures are all zero when the grid is created.  The input file
+is a text file which may have comments starting with the '#' character.
+Each line contains four numeric columns: ix,iy,iz,Temperature.  Empty
+or comment-only lines will be ignored. The
+number of lines must be equal to the number of user-specified grid
+points (Nx by Ny by Nz).  The ix,iy,iz are grid point indices ranging
+from 0 to nxnodes-1 inclusive in each dimension.  The lines can appear
+in any order.  For example, the initial electronic temperatures on a 1
+by 2 by 3 grid could be specified in the file as follows:
 
 .. parsed-literal::
 
+   # UNITS: metal COMMENT: initial electron temperature
    0 0 0 1.0
    0 0 1 1.0
    0 0 2 1.0
@@ -155,40 +209,37 @@ as follows:
 
 where the electronic temperatures along the y=0 plane have been set to
 1.0, and the electronic temperatures along the y=1 plane have been set
-to 2.0.  The order of lines in this file is no important.  If all the
-nodal values are not specified, LAMMPS will generate an error.
+to 2.0.  If all the grid point values are not specified, LAMMPS will
+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.
 
-The temperature output file, *T_oufile*, is created and written by
-this fix.  Temperatures for both the electronic and atomic subsystems
-at every node and every N timesteps are output.  If N is specified as
-zero, no output is generated, and no output filename is needed.  The
-format of the output is as follows.  One long line is written every
-output timestep.  The timestep itself is given in the first column.
-The next Nx\*Ny\*Nz columns contain the temperatures for the atomic
-subsystem, and the final Nx\*Ny\*Nz columns contain the temperatures for
-the electronic subsystem.  The ordering of the Nx\*Ny\*Nz columns is
-with the z index varying fastest, y the next fastest, and x the
-slowest.
+..note::
 
-These fixes do not change the coordinates of their atoms; they only
-scales their velocities.  Thus a time integration fix (e.g. :doc:`fix nve <fix_nve>`) should still be used to time integrate the affected
-atoms.  The fixes should not normally be used on atoms that have their
-temperature controlled by another fix - e.g. :doc:`fix nvt <fix_nh>` or
-:doc:`fix langevin <fix_langevin>`.
+  The electronic temperature at each grid point must be a non-zero
+  positive value, both initially, and as the temperature evovles over
+  time.  Thus you must use either the *set* or *infile* keyword or be
+  restarting a simulation that used this fix previously.
 
-.. note::
+The *outfile* keyword has 2 values.  The first value *Nout* triggers
+output of the electronic temperatures for each grid point every Nout
+timesteps.  The second value is the filename for output which will
+be suffixed by the timestep.  The format of each output file is exactly
+the same as the input temperature file. It will contain a comment in
+the first line reporting the date the file was created, the LAMMPS
+units setting in use, grid size and the current timestep.
 
-   The current implementations of these fixes create a copy of the
-   electron grid that overlays the entire simulation domain, for each
-   processor.  Values on the grid are summed across all processors.  Thus
-   you should insure that this grid is not too large, else your
-   simulation could incur high memory and communication costs.
+Note that the atomic temperature for atoms in each grid cell can also
+be computed and output by the :doc:`fix ave/chunk <fix_ave_chunk>`
+command using the :doc:`compute chunk/atom <compute_chunk_atom>`
+command to create a 3d array of chunks consistent with the grid used
+by this fix.
 
 ----------
 
 **Additional details for fix ttm/mod**
 
-Fix ttm/mod uses the heat diffusion equation with possible external
+Fix *ttm/mod* uses the heat diffusion equation with possible external
 heat sources (e.g. laser heating in ablation simulations):
 
 .. math::
@@ -222,7 +273,8 @@ acting on an ion is:
 
 .. math::
 
-  {\vec F}_i = - \partial U / \partial {\vec r}_i + {\vec F}_{langevin} - \nabla P_e/n_{ion}
+  {\vec F}_i = - \partial U / \partial {\vec r}_i + {\vec
+  F}_{langevin} - \nabla P_e/n_{ion}
 
 where F_langevin is a force from Langevin thermostat simulating
 electron-phonon coupling, and nabla P_e/n_ion is the electron blast
@@ -246,7 +298,9 @@ is calculated as
 
 .. math::
 
-  \nabla_x P_e = \left[\frac{C_e{}T_e(x)\lambda}{(x+\lambda)^2} + \frac{x}{x+\lambda}\frac{(C_e{}T_e)_{x+\Delta x}-(C_e{}T_e)_{x}}{\Delta x} \right]
+  \nabla_x P_e = \left[\frac{C_e{}T_e(x)\lambda}{(x+\lambda)^2} +
+  \frac{x}{x+\lambda}\frac{(C_e{}T_e)_{x+\Delta
+  x}-(C_e{}T_e)_{x}}{\Delta x} \right]
 
 where lambda is the electron mean free path (see :ref:`(Norman) <Norman>`,
 :ref:`(Pisarev) <Pisarev>`)
@@ -286,10 +340,12 @@ Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
 These fixes write the state of the electronic subsystem and the energy
-exchange between the subsystems to :doc:`binary restart files <restart>`.  See the :doc:`read_restart <read_restart>` command
-for info on how to re-specify a fix in an input script that reads a
+exchange between the subsystems to :doc:`binary restart files
+<restart>`.  See the :doc:`read_restart <read_restart>` command for
+info on how to re-specify a fix in an input script that reads a
 restart file, so that the operation of the fix continues in an
-uninterrupted fashion.
+uninterrupted fashion.  Note that the restart script must define the
+same size grid as the original script.
 
 Because the state of the random number generator is not saved in the
 restart files, this means you cannot do "exact" restarts with this
@@ -297,16 +353,16 @@ fix, where the simulation continues on the same as if no restart had
 taken place.  However, in a statistical sense, a restarted simulation
 should produce the same behavior.
 
-None of the :doc:`fix_modify <fix_modify>` options are relevant to these
-fixes.
+None of the :doc:`fix_modify <fix_modify>` options are relevant to
+these fixes.
 
-Both fixes compute 2 output quantities stored in a vector of length 2,
-which can be accessed by various :doc:`output commands <Howto_output>`.
-The first quantity is the total energy of the electronic
-subsystem. The second quantity is the energy transferred from the
-electronic to the atomic subsystem on that timestep. Note that the
-velocity verlet integrator applies the fix ttm forces to the atomic
-subsystem as two half-step velocity updates: one on the current
+These fixes compute 2 output quantities stored in a vector of length
+2, which can be accessed by various :doc:`output commands
+<Howto_output>`.  The first quantity is the total energy of the
+electronic subsystem.  The second quantity is the energy transferred
+from the electronic to the atomic subsystem on that timestep. Note
+that the velocity verlet integrator applies the fix ttm forces to the
+atomic subsystem as two half-step velocity updates: one on the current
 timestep and one on the subsequent timestep.  Consequently, the change
 in the atomic subsystem energy is lagged by half a timestep relative
 to the change in the electronic subsystem energy. As a result of this,
@@ -322,13 +378,12 @@ of the :doc:`run <run>` command.  The fixes are not invoked during
 Restrictions
 """"""""""""
 
-Fix *ttm* and *ttm/mod* are part of the EXTRA-FIX package. They are
-only enabled if LAMMPS was built with that package.
-See the :doc:`Build package <Build_package>` page for more info.
+All these fixes are part of the EXTRA-FIX package. They are only
+enabled if LAMMPS was built with that package.  See the :doc:`Build
+package <Build_package>` page for more info.
 
-These fixes can only be used for 3d simulations and orthogonal
-simulation boxes.  You must also use periodic
-:doc:`boundary <boundary>` conditions.
+As mentioned above, these fixes require 3d simulations and orthogonal
+simulation boxes periodic in all 3 dimensions.
 
 Related commands
 """"""""""""""""

doc/src/group.rst

@@ -41,7 +41,7 @@ Syntax
          keyword = *region* or *var* or *every*
            *region* value = region-ID
            *var* value = name of variable
-           *property* value = name of per-atom property
+           *property* value = name of custom integer or floating point vector
            *every* value = N = update group every this many timesteps
        *static* = no args
 
@@ -226,18 +226,33 @@ simulation runs.  This is in contrast to static groups where atoms are
 permanently assigned to the group.  The way the assignment occurs is
 as follows.  Only atoms in the group specified as the parent group via
 the parent-ID are assigned to the dynamic group before the following
-conditions are applied.  If the *region* keyword is used, atoms not in
-the specified region are removed from the dynamic group.  If the *var*
-keyword is used, the variable name must be an atom-style or
-atomfile-style variable.  The variable is evaluated and atoms whose
-per-atom values are 0.0, are removed from the dynamic group. If the *property*
-keyword is used, the per-atom property name must be a previously defined
-per-atom property.  The per-atom property is evaluated and atoms whose
-values are 0.0 are removed from the dynamic group.
+conditions are applied.
+
+If the *region* keyword is used, atoms not in the specified region are
+removed from the dynamic group.
+
+If the *var* keyword is used, the variable name must be an atom-style
+or atomfile-style variable.  The variable is evaluated and atoms whose
+per-atom values are 0.0, are removed from the dynamic group.
+
+If the *property* keyword is used, the name refers to a custom integer
+or floating point per-atom vector defined via the :doc:`fix
+property/atom <fix_property_atom>` command.  This means the values in
+the vector can be read as part of a data file with the :doc:`read_data
+<read_data>` command or specified with the :doc:`set <set>` command.
+Or accessed and changed via the :doc:`library interface to LAMMPS
+<Howto_library>`, or by styles you add to LAMMPS (pair, fix, compute,
+etc) which access the custom vector and modify its values.  Which
+means the values can be modified between or during simulations.  Atoms
+whose values in the custom vector are zero are removed from the
+dynamic group.  Note that the name of the custom per-atom vector is
+specified just as *name*, not as *i_name* or *d_name* as it is for
+other commands that use different kinds of custom atom vectors or
+arrays as arguments.
 
 The assignment of atoms to a dynamic group is done at the beginning of
-each run and on every timestep that is a multiple of *N*, which is the
-argument for the *every* keyword (N = 1 is the default).  For an
+each run and on every timestep that is a multiple of *N*\ , which is
+the argument for the *every* keyword (N = 1 is the default).  For an
 energy minimization, via the :doc:`minimize <minimize>` command, an
 assignment is made at the beginning of the minimization, but not
 during the iterations of the minimizer.

doc/src/improper_hybrid.rst

@@ -40,13 +40,38 @@ potential with coefficients 120.0, 30 for :math:`K`, :math:`\chi_0`.
 Improper type 2 would be computed with a *cvff* potential with coefficients
 20.0, -1, 2 for K, d, and n, respectively.
 
-If the improper *class2* potential is one of the hybrid styles, it
-requires additional AngleAngle coefficients be specified in the data
-file.  These lines must also have an additional "class2" argument
-added after the improper type.  For improper types which are assigned
-to other hybrid styles, use the style name (e.g. "harmonic")
-appropriate to that style.  The AngleAngle coeffs for that improper
-type will then be ignored.
+If improper coefficients are specified in the data file read via the
+:doc:`read_data <read_data>` command, then the same rule applies.
+E.g. "harmonic" or "cvff", must be added after the improper type, for
+each line in the "Improper Coeffs" section, e.g.
+
+.. parsed-literal::
+
+   Improper Coeffs
+
+   1 harmonic 120.0 30
+   2 cvff 20.0 -1 2
+   ...
+
+If *class2* is one of the improper hybrid styles, the same rule holds
+for specifying additional AngleAngle coefficients either via the input
+script or in the data file. I.e. *class2* must be added to each line
+after the improper type.  For
+lines in the AngleAngle Coeffs section of the data
+file for dihedral types that are not *class2*, you must use an
+improper style of *skip* as a placeholder, e.g.
+
+.. parsed-literal::
+
+   AngleAngle Coeffs
+
+   1 skip
+   2 class2 0.0 0.0 0.0 115.06 130.01 115.06
+   ...
+
+Note that it is not necessary to use the improper style *skip* in the
+input script, since AngleAngle coefficients
+need not be specified at all for improper types that are not *class2*.
 
 An improper style of *none* can be specified as the second argument to
 the improper_coeff command, if you desire to turn off certain improper

doc/src/kspace_style.rst

@@ -2,6 +2,7 @@
 .. index:: kspace_style ewald/dipole
 .. index:: kspace_style ewald/dipole/spin
 .. index:: kspace_style ewald/disp
+.. index:: kspace_style ewald/disp/dipole
 .. index:: kspace_style ewald/omp
 .. index:: kspace_style pppm
 .. index:: kspace_style pppm/kk
@@ -39,7 +40,7 @@ Syntax
 
    kspace_style style value
 
-* style = *none* or *ewald* or *ewald/dipole* or *ewald/dipole/spin* or *ewald/disp* or *ewald/omp* or *pppm* or *pppm/cg* or *pppm/disp* or *pppm/tip4p* or *pppm/stagger* or *pppm/disp/tip4p* or *pppm/gpu* or *pppm/intel* or *pppm/disp/intel* or *pppm/kk* or *pppm/omp* or *pppm/cg/omp* or *pppm/disp/tip4p/omp* or *pppm/tip4p/omp* or *pppm/dielectic* or *pppm/disp/dielectric* or *msm* or *msm/cg* or *msm/omp* or *msm/cg/omp* or *msm/dielectric* or *scafacos*
+* style = *none* or *ewald* or *ewald/dipole* or *ewald/dipole/spin* or *ewald/disp* or *ewald/disp/dipole* or *ewald/omp* or *pppm* or *pppm/cg* or *pppm/disp* or *pppm/tip4p* or *pppm/stagger* or *pppm/disp/tip4p* or *pppm/gpu* or *pppm/intel* or *pppm/disp/intel* or *pppm/kk* or *pppm/omp* or *pppm/cg/omp* or *pppm/disp/tip4p/omp* or *pppm/tip4p/omp* or *pppm/dielectic* or *pppm/disp/dielectric* or *msm* or *msm/cg* or *msm/omp* or *msm/cg/omp* or *msm/dielectric* or *scafacos*
 
   .. parsed-literal::
 
@@ -52,6 +53,8 @@ Syntax
          accuracy = desired relative error in forces
        *ewald/disp* value = accuracy
          accuracy = desired relative error in forces
+       *ewald/disp/dipole* value = accuracy
+         accuracy = desired relative error in forces
        *ewald/omp* value = accuracy
          accuracy = desired relative error in forces
        *pppm* value = accuracy
@@ -156,6 +159,8 @@ matching keyword to the name of the KSpace style, as in this table:
 +----------------------+-----------------------+
 | tip4p/long           | tip4p                 |
 +----------------------+-----------------------+
+| dipole/long          | dipole                |
++----------------------+-----------------------+
 
 ----------
 
@@ -168,7 +173,8 @@ The *ewald/disp* style adds a long-range dispersion sum option for
 but in a more efficient manner than the *ewald* style.  The :math:`1/r^6`
 capability means that Lennard-Jones or Buckingham potentials can be
 used without a cutoff, i.e. they become full long-range potentials.
-The *ewald/disp* style can also be used with point-dipoles, see
+
+The *ewald/disp/dipole* style can also be used with point-dipoles, see
 :ref:`(Toukmaji) <Toukmaji>`.
 
 The *ewald/dipole* style adds long-range standard Ewald summations

doc/src/pair_airebo.rst

@@ -229,6 +229,11 @@ but you would need to create your own AIREBO or AIREBO-M potential file
 with coefficients listed in the appropriate units, if your simulation
 does not use "metal" units.
 
+The pair styles provided here **only** support potential files parameterized
+for the elements carbon and hydrogen (designated with "C" and "H" in the
+*pair_coeff* command.  Using potential files for other elements will trigger
+an error.
+
 Related commands
 """"""""""""""""
 

doc/src/pair_coul.rst

@@ -10,6 +10,7 @@
 .. index:: pair_style coul/dsf/gpu
 .. index:: pair_style coul/dsf/kk
 .. index:: pair_style coul/dsf/omp
+.. index:: pair_style coul/exclude
 .. index:: pair_style coul/cut/global
 .. index:: pair_style coul/cut/global/omp
 .. index:: pair_style coul/long
@@ -42,6 +43,9 @@ pair_style coul/dsf command
 
 Accelerator Variants: *coul/dsf/gpu*, *coul/dsf/kk*, *coul/dsf/omp*
 
+pair_style coul/exclude command
+===============================
+
 pair_style coul/cut/global command
 ==================================
 
@@ -83,6 +87,7 @@ Syntax
    pair_style coul/cut cutoff
    pair_style coul/debye kappa cutoff
    pair_style coul/dsf alpha cutoff
+   pair_style coul/exclude cutoff
    pair_style coul/cut/global cutoff
    pair_style coul/long cutoff
    pair_style coul/wolf alpha cutoff
@@ -110,6 +115,9 @@ Examples
    pair_style coul/dsf 0.05 10.0
    pair_coeff * *
 
+   pair_style hybrid/overlay coul/exclude 10.0 ...
+   pair_coeff * * coul/exclude
+
    pair_style coul/long 10.0
    pair_coeff * *
 
@@ -257,6 +265,19 @@ as style *coul/cut* except that it allows only a single global cutoff
 and thus makes it compatible for use in combination with long-range
 coulomb styles in :doc:`hybrid pair styles <pair_hybrid>`.
 
+Pair style *coul/exclude* computes Coulombic interactions like *coul/cut*
+but **only** applies them to excluded pairs using a scaling factor
+of :math:`\gamma - 1.0` with :math:`\gamma` being the factor assigned
+to that excluded pair via the :doc:`special_bonds coul <special_bonds>`
+setting. With this it is possible to treat Coulomb interactions for
+molecular systems with :doc:`kspace style scafacos <kspace_style>`,
+which always computes the *full* Coulomb interactions without exclusions.
+Pair style *coul/exclude* will then *subtract* the excluded interactions
+accordingly. So to achieve the same forces as with ``pair_style lj/cut/coul/long 12.0``
+with ``kspace_style pppm 1.0e-6``, one would use
+``pair_style hybrid/overlay lj/cut 12.0 coul/exclude 12.0`` with
+``kspace_style scafacos p3m 1.0e-6``.
+
 Styles *coul/long* and *coul/msm* compute the same Coulombic
 interactions as style *coul/cut* except that an additional damping
 factor is applied so it can be used in conjunction with the

doc/src/pair_snap.rst

@@ -135,7 +135,7 @@ with #) anywhere. Each non-blank non-comment line must contain one
 keyword/value pair. The required keywords are *rcutfac* and
 *twojmax*\ . Optional keywords are *rfac0*, *rmin0*,
 *switchflag*, *bzeroflag*, *quadraticflag*, *chemflag*,
-*bnormflag*, *wselfallflag*, and *chunksize*\ .
+*bnormflag*, *wselfallflag*, *chunksize*, and *parallelthresh*\ .
 
 The default values for these keywords are
 
@@ -147,7 +147,8 @@ The default values for these keywords are
 * *chemflag* = 0
 * *bnormflag* = 0
 * *wselfallflag* = 0
-* *chunksize* = 4096
+* *chunksize* = 32768
+* *parallelthresh* = 8192
 
 If *quadraticflag* is set to 1, then the SNAP energy expression includes
 additional quadratic terms that have been shown to increase the overall
@@ -188,14 +189,21 @@ corresponding *K*-vector of linear coefficients for element
 which must equal the number of unique elements appearing in the LAMMPS
 pair_coeff command, to avoid ambiguity in the number of coefficients.
 
-The keyword *chunksize* is only applicable when using the
-pair style *snap* with the KOKKOS package and is ignored otherwise.
-This keyword controls
-the number of atoms in each pass used to compute the bispectrum
-components and is used to avoid running out of memory. For example
-if there are 8192 atoms in the simulation and the *chunksize*
-is set to 4096, the bispectrum calculation will be broken up
-into two passes.
+The keywords *chunksize* and *parallelthresh* are only applicable when
+using the pair style *snap* with the KOKKOS package on GPUs and are
+ignored otherwise.  The *chunksize* keyword controls the number of atoms
+in each pass used to compute the bispectrum components and is used to
+avoid running out of memory.  For example if there are 8192 atoms in the
+simulation and the *chunksize* is set to 4096, the bispectrum
+calculation will be broken up into two passes (running on a single GPU).
+The *parallelthresh* keyword controls a crossover threshold for
+performing extra parallelism.  For small systems, exposing additional
+parallelism can be beneficial when there is not enough work to fully
+saturate the GPU threads otherwise.  However, the extra parallelism also
+leads to more divergence and can hurt performance when the system is
+already large enough to saturate the GPU threads.  Extra parallelism
+will be performed if the *chunksize* (or total number of atoms per GPU)
+is smaller than *parallelthresh*.
 
 Detailed definitions for all the other keywords
 are given on the :doc:`compute sna/atom <compute_sna_atom>` doc page.

doc/src/pair_style.rst

@@ -139,6 +139,7 @@ accelerated styles exist.
 * :doc:`coul/debye <pair_coul>` - cutoff Coulomb potential with Debye screening
 * :doc:`coul/diel <pair_coul_diel>` - Coulomb potential with dielectric permittivity
 * :doc:`coul/dsf <pair_coul>` - Coulomb with damped-shifted-force model
+* :doc:`coul/exclude <pair_coul>` - subtract Coulomb potential for excluded pairs
 * :doc:`coul/long <pair_coul>` - long-range Coulomb potential
 * :doc:`coul/long/cs <pair_cs>` - long-range Coulomb potential and core/shell
 * :doc:`coul/long/dielectric <pair_dielectric>` -

doc/src/read_data.rst

@@ -66,8 +66,8 @@ simulation.  The file can be ASCII text or a gzipped text file
 (detected by a .gz suffix).  This is one of 3 ways to specify initial
 atom coordinates; see the :doc:`read_restart <read_restart>` and
 :doc:`create_atoms <create_atoms>` commands for alternative methods.
-Also see the explanation of the :doc:`-restart command-line switch <Run_options>` which can convert a restart file to a data
-file.
+Also see the explanation of the :doc:`-restart command-line switch
+<Run_options>` which can convert a restart file to a data file.
 
 This command can be used multiple times to add new atoms and their
 properties to an existing system by using the *add*, *offset*, and
@@ -246,22 +246,22 @@ appear in any order, with a few exceptions as noted below.
 
 The keyword *fix* can be used one or more times.  Each usage specifies
 a fix that will be used to process a specific portion of the data
-file.  Any header line containing *header-string* and any section with
-a name containing *section-string* will be passed to the specified
+file.  Any header line containing *header-string* and any section that
+is an exact match to *section-string* will be passed to the specified
 fix.  See the :doc:`fix property/atom <fix_property_atom>` command for
-an example of a fix that operates in this manner.  The page for
-the fix defines the syntax of the header line(s) and section(s) that
-it reads from the data file.  Note that the *header-string* can be
+an example of a fix that operates in this manner.  The doc page for
+the fix defines the syntax of the header line(s) and section that it
+reads from the data file.  Note that the *header-string* can be
 specified as NULL, in which case no header lines are passed to the
-fix.  This means that it can infer the length of its Section from
+fix.  This means the fix can infer the length of its Section from
 standard header settings, such as the number of atoms.
 
 The formatting of individual lines in the data file (indentation,
 spacing between words and numbers) is not important except that header
 and section keywords (e.g. atoms, xlo xhi, Masses, Bond Coeffs) must
-be capitalized as shown and can't have extra white-space between their
-words - e.g. two spaces or a tab between the 2 words in "xlo xhi" or
-the 2 words in "Bond Coeffs", is not valid.
+be capitalized as shown and cannot have extra white-space between
+their words - e.g. two spaces or a tab between the 2 words in "xlo
+xhi" or the 2 words in "Bond Coeffs", is not valid.
 
 ----------
 
@@ -636,12 +636,14 @@ of analysis.
      - atom-ID molecule-ID atom-type lineflag density x y z
    * - mdpd
      - atom-ID atom-type rho x y z
+   * - mesont
+     - atom-ID molecule-ID atom-type bond_nt mass mradius mlength buckling x y z
    * - molecular
      - atom-ID molecule-ID atom-type x y z
    * - peri
      - atom-ID atom-type volume density x y z
    * - smd
-     - atom-ID atom-type molecule volume mass kernel-radius contact-radius x0 y0 z0 x y z
+     - atom-ID atom-type molecule volume mass kradius cradius x0 y0 z0 x y z
    * - sph
      - atom-ID atom-type rho esph cv x y z
    * - sphere
@@ -664,8 +666,10 @@ The per-atom values have these meanings and units, listed alphabetically:
 * atom-ID = integer ID of atom
 * atom-type = type of atom (1-Ntype)
 * bodyflag = 1 for body particles, 0 for point particles
+* bond_nt = bond NT factor for MESONT particles (?? units)
+* buckling = buckling factor for MESONT particles (?? units)
 * ccN = chemical concentration for tDPD particles for each species (mole/volume units)
-* contact-radius = ??? (distance units)
+* cradius = contact radius for SMD particles (distance units)
 * cs_re,cs_im = real/imaginary parts of wave packet coefficients
 * cv = heat capacity (need units) for SPH particles
 * density = density of particle (mass/distance\^3 or mass/distance\^2 or mass/distance units, depending on dimensionality of particle)
@@ -676,10 +680,12 @@ The per-atom values have these meanings and units, listed alphabetically:
 * ellipsoidflag = 1 for ellipsoidal particles, 0 for point particles
 * eradius = electron radius (or fixed-core radius)
 * etag = integer ID of electron that each wave packet belongs to
-* kernel-radius = ??? (distance units)
+* kradius = kernel radius for SMD particles (distance units)
 * lineflag = 1 for line segment particles, 0 for point or spherical particles
 * mass = mass of particle (mass units)
+* mlength = ?? length for MESONT particles (distance units)
 * molecule-ID = integer ID of molecule the atom belongs to
+* mradius = ?? radius for MESONT particles (distance units)
 * mux,muy,muz = components of dipole moment of atom (dipole units)
 * q = charge on atom (charge units)
 * rho = density (need units) for SPH particles

doc/src/rerun.rst

@@ -175,24 +175,28 @@ include and perform all the usual operations of an input script that
 uses the :doc:`run <run>` command.  There are a few exceptions and
 points to consider, as discussed here.
 
-Fixes that perform time integration, such as :doc:`fix nve <fix_nve>` or
-:doc:`fix npt <fix_nh>` are not invoked, since no time integration is
-performed.  Fixes that perturb or constrain the forces on atoms will
-be invoked, just as they would during a normal run.  Examples are :doc:`fix indent <fix_indent>` and :doc:`fix langevin <fix_langevin>`.  So you
-should think carefully as to whether that makes sense for the manner
-in which you are reprocessing the dump snapshots.
+Fixes that perform time integration, such as :doc:`fix nve <fix_nve>`
+or :doc:`fix npt <fix_nh>` are not invoked, since no time integration
+is performed.  Fixes that perturb or constrain the forces on atoms
+will be invoked, just as they would during a normal run.  Examples are
+:doc:`fix indent <fix_indent>` and :doc:`fix langevin <fix_langevin>`.
+So you should think carefully as to whether that makes sense for the
+manner in which you are reprocessing the dump snapshots.
 
-If you only want the rerun script to perform an analysis that does
-not involve pair interactions, such as use compute msd to calculated
-displacements over time, you do not need to define a :doc:`pair style <pair_style>`, which may also mean neighbor lists will not
-need to be calculated which saves time.  The :doc:`comm_modify cutoff <comm_modify>` command can also be used to insure ghost
-atoms are acquired from far enough away for operations like bond and
-angle evaluations, if no pair style is being used.
+If you only want the rerun script to perform an analysis that does not
+involve pair interactions, such as use compute msd to calculated
+displacements over time, you do not need to define a :doc:`pair style
+<pair_style>`, which may also mean neighbor lists will not need to be
+calculated which saves time.  The :doc:`comm_modify cutoff
+<comm_modify>` command can also be used to insure ghost atoms are
+acquired from far enough away for operations like bond and angle
+evaluations, if no pair style is being used.
 
 Every time a snapshot is read, the timestep for the simulation is
 reset, as if the :doc:`reset_timestep <reset_timestep>` command were
 used.  This command has some restrictions as to what fixes can be
-defined.  See its page for details.  For example, the :doc:`fix deposit <fix_deposit>` and :doc:`fix dt/reset <fix_dt_reset>` fixes
+defined.  See its page for details.  For example, the :doc:`fix
+deposit <fix_deposit>` and :doc:`fix dt/reset <fix_dt_reset>` fixes
 are in this category.  They also make no sense to use with a rerun
 command.
 

doc/src/set.rst

@@ -13,16 +13,17 @@ Syntax
 * style = *atom* or *type* or *mol* or *group* or *region*
 * ID = atom ID range or type range or mol ID range or group ID or region ID
 * one or more keyword/value pairs may be appended
-* keyword = *type* or *type/fraction* or *type/ratio* or *type/subset* or *mol*
-  or *x* or *y* or *z* or *vx* or *vy* or *vz*
-  or *charge* or *dipole* or *dipole/random* or *spin* or *spin/random*
-  or *quat* *quat/random* or *diameter* or *shape* or *length* or *tri*
-  or *theta* or *theta/random* or *angmom* or *omega* or *mass*
-  or *density* or *density/disc* or *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 *i_name* or *d_name*
+* keyword = *type* or *type/fraction* or *type/ratio* or *type/subset*
+  or *mol* or *x* or *y* or *z* or *charge* or *dipole* or
+  *dipole/random* or *quat* or *spin* or *spin/random* or
+  *quat* or *quat/random* or *diameter* or *shape* or
+  *length* or *tri* or *theta* or *theta/random* or *angmom* or
+  *omega* or *mass* or *density* or *density/disc* or
+  *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
+  *i_name* or *d_name* or *i2_name* or *d2_name*
 
   .. parsed-literal::
 
@@ -123,8 +124,12 @@ 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)
-       *i_name* value = value for custom integer vector with name
-       *d_name* value = value for custom floating-point vector with name
+       *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
+                         column specified as i2_name[N] where N is 1 to Ncol
+       *d2_name* value = column of a custom floating-point array with name
+                         column specified as d2_name[N] where N is 1 to Ncol
 
 Examples
 """"""""
@@ -141,6 +146,8 @@ Examples
    set atom 100*200 x 0.5 y 1.0
    set atom 100 vx 0.0 vy 0.0 vz -1.0
    set atom 1492 type 3
+   set atom * i_myVal 5
+   set atom * d2_Sxyz[1] 6.4
 
 Description
 """""""""""
@@ -482,11 +489,13 @@ attribute. An integer for "index" selects a chemical species (1 to
 Nspecies) where Nspecies is set by the atom_style command. The value
 for the chemical concentration must be >= 0.0.
 
-Keywords *i_name* and *d_name* refer to custom integer and
-floating-point properties that have been added to each atom via the
-:doc:`fix property/atom <fix_property_atom>` command.  When that command
-is used specific names are given to each attribute which are what is
-specified as the "name" portion of *i_name* or *d_name*.
+Keywords *i_name*, *d_name*, *i2_name*, *d2_name* refer to custom
+per-atom integer and floating-point vectors or arrays that have been
+added via the :doc:`fix property/atom <fix_property_atom>` command.
+When that command is used specific names are given to each attribute
+which are the "name" portion of these keywords.  For arrays *i2_name*
+and *d2_name*, the column of the array must also be included following
+the name in brackets: e.g. d2_xyz[2], i2_mySpin[3].
 
 Restrictions
 """"""""""""

doc/src/velocity.rst

@@ -92,25 +92,25 @@ velocity field.
 The *scale* style computes the current temperature of the group of
 atoms and then rescales the velocities to the specified temperature.
 
-The *ramp* style is similar to that used by the :doc:`compute temp/ramp <compute_temp_ramp>` command.  Velocities ramped
-uniformly from vlo to vhi are applied to dimension vx, or vy, or vz.
-The value assigned to a particular atom depends on its relative
-coordinate value (in dim) from clo to chi.  For the example above, an
-atom with y-coordinate of 10 (1/4 of the way from 5 to 25), would be
-assigned a x-velocity of 1.25 (1/4 of the way from 0.0 to 5.0).  Atoms
-outside the coordinate bounds (less than 5 or greater than 25 in this
-case), are assigned velocities equal to vlo or vhi (0.0 or 5.0 in this
-case).
+The *ramp* style is similar to that used by the :doc:`compute
+temp/ramp <compute_temp_ramp>` command.  Velocities ramped uniformly
+from vlo to vhi are applied to dimension vx, or vy, or vz.  The value
+assigned to a particular atom depends on its relative coordinate value
+(in dim) from clo to chi.  For the example above, an atom with
+y-coordinate of 10 (1/4 of the way from 5 to 25), would be assigned a
+x-velocity of 1.25 (1/4 of the way from 0.0 to 5.0).  Atoms outside
+the coordinate bounds (less than 5 or greater than 25 in this case),
+are assigned velocities equal to vlo or vhi (0.0 or 5.0 in this case).
 
 The *zero* style adjusts the velocities of the group of atoms so that
 the aggregate linear or angular momentum is zero.  No other changes
 are made to the velocities of the atoms.  If the *rigid* option is
 specified (see below), then the zeroing is performed on individual
-rigid bodies, as defined by the :doc:`fix rigid or fix rigid/small <fix_rigid>` commands.  In other words, zero linear
-will set the linear momentum of each rigid body to zero, and zero
-angular will set the angular momentum of each rigid body to zero.
-This is done by adjusting the velocities of the atoms in each rigid
-body.
+rigid bodies, as defined by the :doc:`fix rigid or fix rigid/small
+<fix_rigid>` commands.  In other words, zero linear will set the
+linear momentum of each rigid body to zero, and zero angular will set
+the angular momentum of each rigid body to zero.  This is done by
+adjusting the velocities of the atoms in each rigid body.
 
 All temperatures specified in the velocity command are in temperature
 units; see the :doc:`units <units>` command.  The units of velocities and

doc/utils/sphinx-config/false_positives.txt

@@ -531,6 +531,7 @@ covalently
 covariance
 cpp
 cpu
+cradius
 createatoms
 createAtoms
 CreateIDs
@@ -1173,6 +1174,7 @@ googletest
 Gordan
 Goudeau
 GPa
+GPL
 gpu
 gpuID
 gpus
@@ -1598,6 +1600,7 @@ Kosovan
 Koster
 Kosztin
 Kp
+kradius
 Kraker
 Kraus
 Kremer
@@ -1687,6 +1690,7 @@ Lett
 Leuven
 Leven
 Lewy
+LGPL
 lgvdw
 Liang
 libatc
@@ -1887,7 +1891,6 @@ maxX
 Mayergoyz
 Mayoral
 mbt
-Mbytes
 MBytes
 mc
 McLachlan
@@ -1993,6 +1996,7 @@ Mj
 mK
 mkdir
 mkv
+mlength
 mliap
 mliappy
 mlparks
@@ -2049,6 +2053,7 @@ mpiio
 mpirun
 mplayer
 mps
+mradius
 Mrovec
 Mryglod
 mscg
@@ -2107,6 +2112,7 @@ myIndex
 mylammps
 MyPool
 mysocket
+mySpin
 myTemp
 myVec
 na
@@ -2275,6 +2281,7 @@ normals
 Noskov
 noslip
 noticable
+Nout
 noutcol
 Noutput
 noutrow
@@ -2317,6 +2324,7 @@ nsub
 Nswap
 Nt
 Ntable
+nt
 ntheta
 nthreads
 ntimestep

examples/PACKAGES/interlayer/kolmogorov_crespi_full/CC.KC-full

@@ -1 +0,0 @@
-../../../../potentials/CC.KC-full

examples/PACKAGES/interlayer/kolmogorov_crespi_full/CH.KC

@@ -1 +0,0 @@
-../../../../potentials/CH.KC

examples/PACKAGES/interlayer/lebedeva/lebedeva00.plot

@@ -28,7 +28,7 @@ set yrange [-0.002:0.001]
 #set yrange [*:*]
 
 plot \
-	"LamppsResult.dat" u 2:5 t "Leb LAMMPS",\
+	"LammpsResult.dat" u 2:5 t "Leb LAMMPS",\
 	"PerlResult.dat" u 1:($4*0.001/2.) w l t "Leb Perl"
 
 exit

examples/README

@@ -114,6 +114,7 @@ tad:      temperature-accelerated dynamics of vacancy diffusion in bulk Si
 template: examples for using atom_style template and comparing to atom style molecular
 tersoff:  regression test input for Tersoff variants
 threebody: regression test input for a variety of threebody potentials
+ttm:      two-temeperature model examples
 vashishta: models using the Vashishta potential
 voronoi:  Voronoi tesselation via compute voronoi/atom command
 wall:     use of reflective walls with different stochastic models

examples/ttm/in.ttm

@@ -0,0 +1,43 @@
+units           metal
+atom_style      atomic
+boundary        p p p
+
+variable        latc  equal 2.87
+lattice         bcc ${latc}
+variable 	xmax equal   10.0
+variable 	xmin equal  -10.0
+variable 	ymax equal   10.0
+variable 	ymin equal  -10.0
+variable 	zmax equal   10.0
+variable 	zmin equal  -10.0
+
+region 		sim_box block ${xmin} ${xmax} ${ymin} ${ymax} ${zmin} ${zmax} &
+                units lattice
+create_box 	1 sim_box
+region 		atom_box block ${xmin} ${xmax} ${ymin} ${ymax} ${zmin} ${zmax} &
+                units lattice
+create_atoms 	1 region atom_box
+
+mass            1 55.845
+
+pair_style      eam/fs
+pair_coeff      * * FeVoter-ChenRecheck.fs Fe
+
+neighbor        2.0 bin
+neigh_modify    every 5 delay 0 check yes
+
+fix             1 all nve
+
+fix             twotemp all ttm 342785 1.2470e-5 0.087614 &
+                0.005365 29.5917 47.5679 58.4613 10 10 10 set 1800.0
+
+compute         pe all pe/atom
+compute         ke all ke/atom
+
+timestep        0.0001
+thermo          100
+
+thermo_style    custom step temp etotal f_twotemp[1] f_twotemp[2]
+                thermo_modify format float "%20.16g"
+
+run             1000

examples/ttm/in.ttm.grid

@@ -0,0 +1,43 @@
+units           metal
+atom_style      atomic
+boundary        p p p
+
+variable        latc  equal 2.87
+lattice         bcc ${latc}
+variable 	xmax equal   10.0
+variable 	xmin equal  -10.0
+variable 	ymax equal   10.0
+variable 	ymin equal  -10.0
+variable 	zmax equal   10.0
+variable 	zmin equal  -10.0
+
+region 		sim_box block ${xmin} ${xmax} ${ymin} ${ymax} ${zmin} ${zmax} &
+                units lattice
+create_box 	1 sim_box
+region 		atom_box block ${xmin} ${xmax} ${ymin} ${ymax} ${zmin} ${zmax} &
+                units lattice
+create_atoms 	1 region atom_box
+
+mass            1 55.845
+
+pair_style      eam/fs
+pair_coeff      * * FeVoter-ChenRecheck.fs Fe
+
+neighbor        2.0 bin
+neigh_modify    every 5 delay 0 check yes
+
+fix             1 all nve
+
+fix             twotemp all ttm/grid 342785 1.2470e-5 0.087614 &
+                0.005365 29.5917 47.5679 58.4613 10 10 10 set 1800.0
+
+compute         pe all pe/atom
+compute         ke all ke/atom
+
+timestep        0.0001
+thermo          100
+
+thermo_style    custom step temp etotal f_twotemp[1] f_twotemp[2]
+                thermo_modify format float "%20.16g"
+
+run             1000

examples/ttm/log.26Aug21.ttm.g++.1

@@ -0,0 +1,113 @@
+LAMMPS (30 Jul 2021)
+units           metal
+atom_style      atomic
+boundary        p p p
+
+variable        latc  equal 2.87
+lattice         bcc ${latc}
+lattice         bcc 2.87
+Lattice spacing in x,y,z = 2.8700000 2.8700000 2.8700000
+variable 	xmax equal   10.0
+variable 	xmin equal  -10.0
+variable 	ymax equal   10.0
+variable 	ymin equal  -10.0
+variable 	zmax equal   10.0
+variable 	zmin equal  -10.0
+
+region 		sim_box block ${xmin} ${xmax} ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 ${xmax} ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 10 ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 10 -10 ${ymax} ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 10 -10 10 ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 10 -10 10 -10 ${zmax}                 units lattice
+region 		sim_box block -10 10 -10 10 -10 10                 units lattice
+create_box 	1 sim_box
+Created orthogonal box = (-28.700000 -28.700000 -28.700000) to (28.700000 28.700000 28.700000)
+  1 by 1 by 1 MPI processor grid
+region 		atom_box block ${xmin} ${xmax} ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 ${xmax} ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 10 ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 10 -10 ${ymax} ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 10 -10 10 ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 10 -10 10 -10 ${zmax}                 units lattice
+region 		atom_box block -10 10 -10 10 -10 10                 units lattice
+create_atoms 	1 region atom_box
+Created 16000 atoms
+  using lattice units in orthogonal box = (-28.700000 -28.700000 -28.700000) to (28.700000 28.700000 28.700000)
+  create_atoms CPU = 0.005 seconds
+
+mass            1 55.845
+
+include		pot_iron.mod
+pair_style      eam/fs
+pair_coeff      * * FeVoter-ChenRecheck.fs Fe
+
+neighbor        2.0 bin
+neigh_modify    every 5 delay 0 check yes
+
+fix             1 all nve
+
+fix             twotemp all ttm 342785 1.2470e-5 0.087614                 0.005365 29.5917 47.5679 58.4613 10 10 10 set 1800.0
+
+compute         pe all pe/atom
+compute         ke all ke/atom
+
+timestep        0.0001
+thermo          100
+
+thermo_style    custom step temp etotal f_twotemp[1] f_twotemp[2]
+                thermo_modify format float "%20.16g"
+
+run             1000
+Neighbor list info ...
+  update every 5 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 6.524
+  ghost atom cutoff = 6.524
+  binsize = 3.262, bins = 18 18 18
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair eam/fs, perpetual
+      attributes: half, newton on
+      pair build: half/bin/atomonly/newton
+      stencil: half/bin/3d
+      bin: standard
+Per MPI rank memory allocation (min/avg/max) = 10.97 | 10.97 | 10.97 Mbytes
+Step Temp TotEng f_twotemp[1] f_twotemp[2] 
+       0                    0   -68483.52254543516    371.9188105082105                    0 
+     100    17.01353086098387   -68446.50228930202    334.6217068813629   0.3763710887774045 
+     200    27.91331236535322   -68413.16008042906    301.3181773007303   0.3165912892484031 
+     300    32.20115656493127   -68383.19634217303    271.3756381838044   0.2901111802983097 
+     400    33.46056398887347   -68355.73057141017    243.9344715501159   0.2548133388123376 
+     500     35.5346204243821   -68331.63060947017      219.84946888619   0.2388591367999412 
+     600    40.61692458441595   -68309.36124792948    197.5527667607885   0.3056696014124333 
+     700    46.20303146200327   -68290.12727395596    178.3775768561404   0.1982123493608405 
+     800    50.43750181899325   -68272.72651051797     160.995046695269   0.1708386295858842 
+     900     52.1701171463511   -68257.85059865141    146.1567281868866   0.1032829304640773 
+    1000    53.49296457217385   -68244.38715993935    132.7166474251701  0.06428993394665879 
+Loop time of 17.1447 on 1 procs for 1000 steps with 16000 atoms
+
+Performance: 0.504 ns/day, 47.624 hours/ns, 58.327 timesteps/s
+99.9% CPU use with 1 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 15.811     | 15.811     | 15.811     |   0.0 | 92.22
+Neigh   | 0          | 0          | 0          |   0.0 |  0.00
+Comm    | 0.094539   | 0.094539   | 0.094539   |   0.0 |  0.55
+Output  | 0.00093974 | 0.00093974 | 0.00093974 |   0.0 |  0.01
+Modify  | 1.1898     | 1.1898     | 1.1898     |   0.0 |  6.94
+Other   |            | 0.04797    |            |       |  0.28
+
+Nlocal:        16000.0 ave       16000 max       16000 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:        13449.0 ave       13449 max       13449 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:        896000.0 ave      896000 max      896000 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 896000
+Ave neighs/atom = 56.000000
+Neighbor list builds = 0
+Dangerous builds = 0
+Total wall time: 0:00:17

examples/ttm/log.26Aug21.ttm.g++.4

@@ -0,0 +1,113 @@
+LAMMPS (30 Jul 2021)
+units           metal
+atom_style      atomic
+boundary        p p p
+
+variable        latc  equal 2.87
+lattice         bcc ${latc}
+lattice         bcc 2.87
+Lattice spacing in x,y,z = 2.8700000 2.8700000 2.8700000
+variable 	xmax equal   10.0
+variable 	xmin equal  -10.0
+variable 	ymax equal   10.0
+variable 	ymin equal  -10.0
+variable 	zmax equal   10.0
+variable 	zmin equal  -10.0
+
+region 		sim_box block ${xmin} ${xmax} ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 ${xmax} ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 10 ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 10 -10 ${ymax} ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 10 -10 10 ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 10 -10 10 -10 ${zmax}                 units lattice
+region 		sim_box block -10 10 -10 10 -10 10                 units lattice
+create_box 	1 sim_box
+Created orthogonal box = (-28.700000 -28.700000 -28.700000) to (28.700000 28.700000 28.700000)
+  1 by 2 by 2 MPI processor grid
+region 		atom_box block ${xmin} ${xmax} ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 ${xmax} ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 10 ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 10 -10 ${ymax} ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 10 -10 10 ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 10 -10 10 -10 ${zmax}                 units lattice
+region 		atom_box block -10 10 -10 10 -10 10                 units lattice
+create_atoms 	1 region atom_box
+Created 16000 atoms
+  using lattice units in orthogonal box = (-28.700000 -28.700000 -28.700000) to (28.700000 28.700000 28.700000)
+  create_atoms CPU = 0.002 seconds
+
+mass            1 55.845
+
+include		pot_iron.mod
+pair_style      eam/fs
+pair_coeff      * * FeVoter-ChenRecheck.fs Fe
+
+neighbor        2.0 bin
+neigh_modify    every 5 delay 0 check yes
+
+fix             1 all nve
+
+fix             twotemp all ttm 342785 1.2470e-5 0.087614                 0.005365 29.5917 47.5679 58.4613 10 10 10 set 1800.0
+
+compute         pe all pe/atom
+compute         ke all ke/atom
+
+timestep        0.0001
+thermo          100
+
+thermo_style    custom step temp etotal f_twotemp[1] f_twotemp[2]
+                thermo_modify format float "%20.16g"
+
+run             1000
+Neighbor list info ...
+  update every 5 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 6.524
+  ghost atom cutoff = 6.524
+  binsize = 3.262, bins = 18 18 18
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair eam/fs, perpetual
+      attributes: half, newton on
+      pair build: half/bin/atomonly/newton
+      stencil: half/bin/3d
+      bin: standard
+Per MPI rank memory allocation (min/avg/max) = 4.845 | 4.845 | 4.845 Mbytes
+Step Temp TotEng f_twotemp[1] f_twotemp[2] 
+       0                    0   -68483.52254530673    371.9188105082105                    0 
+     100    16.95536995775683   -68446.64765713879    334.7745598327934   0.3602932995006091 
+     200    27.82619298359641   -68413.48663012494    301.6568409464845   0.2922875754523596 
+     300      32.286609763559   -68383.41369945828    271.6030085280584   0.2698738824780399 
+     400    33.33119316198579   -68356.74598240001    244.9747750036312   0.2061586600914007 
+     500    35.14534756499593   -68332.73504057307    220.9328922343961   0.2800368538794571 
+     600    39.58922469808519   -68311.03191758461    199.2602622784512    0.231030319616688 
+     700    45.34652315787152   -68291.65247941406    179.9297699858465   0.1438135463248857 
+     800    49.66707856481077   -68274.98092841901    163.2540575286428   0.1600890300738259 
+     900    52.17692450487316    -68259.8031091165    148.1017576370546   0.1177316234407944 
+    1000    54.24228199265477   -68245.58589458199    133.8816957314364   0.1314999893461338 
+Loop time of 5.03135 on 4 procs for 1000 steps with 16000 atoms
+
+Performance: 1.717 ns/day, 13.976 hours/ns, 198.754 timesteps/s
+98.8% CPU use with 4 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 4.1173     | 4.2634     | 4.3858     |   5.4 | 84.74
+Neigh   | 0          | 0          | 0          |   0.0 |  0.00
+Comm    | 0.2218     | 0.34547    | 0.49422    |  19.4 |  6.87
+Output  | 0.00031185 | 0.00036952 | 0.00044986 |   0.0 |  0.01
+Modify  | 0.39294    | 0.39605    | 0.39877    |   0.4 |  7.87
+Other   |            | 0.02604    |            |       |  0.52
+
+Nlocal:        4000.00 ave        4000 max        4000 min
+Histogram: 4 0 0 0 0 0 0 0 0 0
+Nghost:        6329.00 ave        6329 max        6329 min
+Histogram: 4 0 0 0 0 0 0 0 0 0
+Neighs:        224000.0 ave      227436 max      220450 min
+Histogram: 1 0 1 0 0 0 0 1 0 1
+
+Total # of neighbors = 896000
+Ave neighs/atom = 56.000000
+Neighbor list builds = 0
+Dangerous builds = 0
+Total wall time: 0:00:05

examples/ttm/log.26Aug21.ttm.grid.g++.1

@@ -0,0 +1,113 @@
+LAMMPS (30 Jul 2021)
+units           metal
+atom_style      atomic
+boundary        p p p
+
+variable        latc  equal 2.87
+lattice         bcc ${latc}
+lattice         bcc 2.87
+Lattice spacing in x,y,z = 2.8700000 2.8700000 2.8700000
+variable 	xmax equal   10.0
+variable 	xmin equal  -10.0
+variable 	ymax equal   10.0
+variable 	ymin equal  -10.0
+variable 	zmax equal   10.0
+variable 	zmin equal  -10.0
+
+region 		sim_box block ${xmin} ${xmax} ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 ${xmax} ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 10 ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 10 -10 ${ymax} ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 10 -10 10 ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 10 -10 10 -10 ${zmax}                 units lattice
+region 		sim_box block -10 10 -10 10 -10 10                 units lattice
+create_box 	1 sim_box
+Created orthogonal box = (-28.700000 -28.700000 -28.700000) to (28.700000 28.700000 28.700000)
+  1 by 1 by 1 MPI processor grid
+region 		atom_box block ${xmin} ${xmax} ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 ${xmax} ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 10 ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 10 -10 ${ymax} ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 10 -10 10 ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 10 -10 10 -10 ${zmax}                 units lattice
+region 		atom_box block -10 10 -10 10 -10 10                 units lattice
+create_atoms 	1 region atom_box
+Created 16000 atoms
+  using lattice units in orthogonal box = (-28.700000 -28.700000 -28.700000) to (28.700000 28.700000 28.700000)
+  create_atoms CPU = 0.005 seconds
+
+mass            1 55.845
+
+include		pot_iron.mod
+pair_style      eam/fs
+pair_coeff      * * FeVoter-ChenRecheck.fs Fe
+
+neighbor        2.0 bin
+neigh_modify    every 5 delay 0 check yes
+
+fix             1 all nve
+
+fix             twotemp all ttm/grid 342785 1.2470e-5 0.087614                 0.005365 29.5917 47.5679 58.4613 10 10 10 set 1800.0
+
+compute         pe all pe/atom
+compute         ke all ke/atom
+
+timestep        0.0001
+thermo          100
+
+thermo_style    custom step temp etotal f_twotemp[1] f_twotemp[2]
+                thermo_modify format float "%20.16g"
+
+run             1000
+Neighbor list info ...
+  update every 5 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 6.524
+  ghost atom cutoff = 6.524
+  binsize = 3.262, bins = 18 18 18
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair eam/fs, perpetual
+      attributes: half, newton on
+      pair build: half/bin/atomonly/newton
+      stencil: half/bin/3d
+      bin: standard
+Per MPI rank memory allocation (min/avg/max) = 10.99 | 10.99 | 10.99 Mbytes
+Step Temp TotEng f_twotemp[1] f_twotemp[2] 
+       0                    0   -68483.52254543516    371.9188105082105                    0 
+     100    17.01353086098387   -68446.50228930202    334.6217068813629   0.3763710887774046 
+     200    27.91331236535322   -68413.16008042906    301.3181773007303   0.3165912892484031 
+     300    32.20115656493125   -68383.19634217303    271.3756381838045   0.2901111802983097 
+     400    33.46056398887347   -68355.73057141017    243.9344715501159   0.2548133388123378 
+     500     35.5346204243821   -68331.63060947017      219.84946888619   0.2388591367999414 
+     600    40.61692458441596   -68309.36124792948    197.5527667607886   0.3056696014124338 
+     700    46.20303146200326   -68290.12727395598    178.3775768561405   0.1982123493608406 
+     800     50.4375018189932   -68272.72651051797     160.995046695269   0.1708386295858845 
+     900    52.17011714635106   -68257.85059865142    146.1567281868867   0.1032829304640776 
+    1000    53.49296457217382   -68244.38715993936    132.7166474251702  0.06428993394665769 
+Loop time of 14.8767 on 1 procs for 1000 steps with 16000 atoms
+
+Performance: 0.581 ns/day, 41.324 hours/ns, 67.219 timesteps/s
+100.0% CPU use with 1 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 13.759     | 13.759     | 13.759     |   0.0 | 92.49
+Neigh   | 0          | 0          | 0          |   0.0 |  0.00
+Comm    | 0.087429   | 0.087429   | 0.087429   |   0.0 |  0.59
+Output  | 0.00063941 | 0.00063941 | 0.00063941 |   0.0 |  0.00
+Modify  | 0.98561    | 0.98561    | 0.98561    |   0.0 |  6.63
+Other   |            | 0.04396    |            |       |  0.30
+
+Nlocal:        16000.0 ave       16000 max       16000 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:        13449.0 ave       13449 max       13449 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:        896000.0 ave      896000 max      896000 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 896000
+Ave neighs/atom = 56.000000
+Neighbor list builds = 0
+Dangerous builds = 0
+Total wall time: 0:00:14

examples/ttm/log.26Aug21.ttm.grid.g++.4

@@ -0,0 +1,113 @@
+LAMMPS (30 Jul 2021)
+units           metal
+atom_style      atomic
+boundary        p p p
+
+variable        latc  equal 2.87
+lattice         bcc ${latc}
+lattice         bcc 2.87
+Lattice spacing in x,y,z = 2.8700000 2.8700000 2.8700000
+variable 	xmax equal   10.0
+variable 	xmin equal  -10.0
+variable 	ymax equal   10.0
+variable 	ymin equal  -10.0
+variable 	zmax equal   10.0
+variable 	zmin equal  -10.0
+
+region 		sim_box block ${xmin} ${xmax} ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 ${xmax} ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 10 ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 10 -10 ${ymax} ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 10 -10 10 ${zmin} ${zmax}                 units lattice
+region 		sim_box block -10 10 -10 10 -10 ${zmax}                 units lattice
+region 		sim_box block -10 10 -10 10 -10 10                 units lattice
+create_box 	1 sim_box
+Created orthogonal box = (-28.700000 -28.700000 -28.700000) to (28.700000 28.700000 28.700000)
+  1 by 2 by 2 MPI processor grid
+region 		atom_box block ${xmin} ${xmax} ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 ${xmax} ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 10 ${ymin} ${ymax} ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 10 -10 ${ymax} ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 10 -10 10 ${zmin} ${zmax}                 units lattice
+region 		atom_box block -10 10 -10 10 -10 ${zmax}                 units lattice
+region 		atom_box block -10 10 -10 10 -10 10                 units lattice
+create_atoms 	1 region atom_box
+Created 16000 atoms
+  using lattice units in orthogonal box = (-28.700000 -28.700000 -28.700000) to (28.700000 28.700000 28.700000)
+  create_atoms CPU = 0.002 seconds
+
+mass            1 55.845
+
+include		pot_iron.mod
+pair_style      eam/fs
+pair_coeff      * * FeVoter-ChenRecheck.fs Fe
+
+neighbor        2.0 bin
+neigh_modify    every 5 delay 0 check yes
+
+fix             1 all nve
+
+fix             twotemp all ttm/grid 342785 1.2470e-5 0.087614                 0.005365 29.5917 47.5679 58.4613 10 10 10 set 1800.0
+
+compute         pe all pe/atom
+compute         ke all ke/atom
+
+timestep        0.0001
+thermo          100
+
+thermo_style    custom step temp etotal f_twotemp[1] f_twotemp[2]
+                thermo_modify format float "%20.16g"
+
+run             1000
+Neighbor list info ...
+  update every 5 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 6.524
+  ghost atom cutoff = 6.524
+  binsize = 3.262, bins = 18 18 18
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair eam/fs, perpetual
+      attributes: half, newton on
+      pair build: half/bin/atomonly/newton
+      stencil: half/bin/3d
+      bin: standard
+Per MPI rank memory allocation (min/avg/max) = 4.843 | 4.843 | 4.843 Mbytes
+Step Temp TotEng f_twotemp[1] f_twotemp[2] 
+       0                    0   -68483.52254530673    371.9188105082186                    0 
+     100    16.95536995775684   -68446.64765713879    334.7745598327931   0.3602932995006087 
+     200    27.82619298359641   -68413.48663012494    301.6568409464847   0.2922875754523593 
+     300    32.28660976355901   -68383.41369945828    271.6030085280586     0.26987388247804 
+     400    33.33119316198579   -68356.74598240001    244.9747750036311   0.2061586600914003 
+     500    35.14534756499593   -68332.73504057307    220.9328922343961   0.2800368538794578 
+     600    39.58922469808521   -68311.03191758461    199.2602622784512   0.2310303196166884 
+     700    45.34652315787151   -68291.65247941404    179.9297699858464   0.1438135463248855 
+     800    49.66707856481075   -68274.98092841901    163.2540575286425   0.1600890300738265 
+     900    52.17692450487317    -68259.8031091165    148.1017576370548   0.1177316234407941 
+    1000    54.24228199265479   -68245.58589458198    133.8816957314364   0.1314999893461343 
+Loop time of 4.95173 on 4 procs for 1000 steps with 16000 atoms
+
+Performance: 1.745 ns/day, 13.755 hours/ns, 201.949 timesteps/s
+98.2% CPU use with 4 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 4.1159     | 4.2665     | 4.4446     |   7.1 | 86.16
+Neigh   | 0          | 0          | 0          |   0.0 |  0.00
+Comm    | 0.14618    | 0.32518    | 0.47663    |  25.8 |  6.57
+Output  | 0.00030633 | 0.00034695 | 0.00044693 |   0.0 |  0.01
+Modify  | 0.33258    | 0.3333     | 0.33402    |   0.1 |  6.73
+Other   |            | 0.0264     |            |       |  0.53
+
+Nlocal:        4000.00 ave        4000 max        4000 min
+Histogram: 4 0 0 0 0 0 0 0 0 0
+Nghost:        6329.00 ave        6329 max        6329 min
+Histogram: 4 0 0 0 0 0 0 0 0 0
+Neighs:        224000.0 ave      227436 max      220450 min
+Histogram: 1 0 1 0 0 0 0 1 0 1
+
+Total # of neighbors = 896000
+Ave neighs/atom = 56.000000
+Neighbor list builds = 0
+Dangerous builds = 0
+Total wall time: 0:00:05

lib/README

@@ -3,7 +3,7 @@ LAMMPS, if particular packages are included in the LAMMPS build.
 
 Most of these directories contain code for the library; some contain
 a Makefile.lammps file that points to where the library is installed
-elsewhere on your system.  
+elsewhere on your system.
 
 In either case, the library itself must be installed and/or built
 first, so that the appropriate library files exist for LAMMPS to link

lib/atc/ATC_Coupling.h

@@ -20,16 +20,16 @@ namespace ATC {
 
   /**
    *  @class  ATC_Coupling
-   *  @brief  Base class for atom-continuum coupling 
+   *  @brief  Base class for atom-continuum coupling
    */
 
   class ATC_Coupling : public ATC_Method {
 
   public: /** methods */
-    
+
     friend class ExtrinsicModel; // friend is not inherited
     friend class ExtrinsicModelTwoTemperature;
-    friend class ExtrinsicModelDriftDiffusion; 
+    friend class ExtrinsicModelDriftDiffusion;
     friend class ExtrinsicModelDriftDiffusionConvection;
     friend class ExtrinsicModelElectrostatic;
     friend class ExtrinsicModelElectrostaticMomentum;
@@ -72,13 +72,13 @@ namespace ATC {
     virtual void pre_init_integrate();
 
     /** Predictor phase, Verlet first step for velocity and position */
-    virtual void init_integrate(); 
+    virtual void init_integrate();
 
     /** Predictor phase, executed after Verlet */
     virtual void post_init_integrate();
 
     /** Corrector phase, executed after Verlet*/
-    
+
     virtual void post_final_integrate();
 
     /** pre/post atomic force calculation in minimize */
@@ -95,9 +95,9 @@ namespace ATC {
     const std::set<std::string> & boundary_face_names() {return boundaryFaceNames_;};
     /** access to boundary integration method */
     int boundary_integration_type() {return bndyIntType_;};
-    void set_boundary_integration_type(int boundaryIntegrationType) 
+    void set_boundary_integration_type(int boundaryIntegrationType)
     {bndyIntType_ = boundaryIntegrationType;};
-    void set_boundary_face_set(const std::set< std::pair<int,int> > * boundaryFaceSet) 
+    void set_boundary_face_set(const std::set< std::pair<int,int> > * boundaryFaceSet)
     {bndyFaceSet_ = boundaryFaceSet;};
     BoundaryIntegrationType parse_boundary_integration
       (int narg, char **arg, const std::set< std::pair<int,int> > * boundaryFaceSet);
@@ -107,10 +107,10 @@ namespace ATC {
 //--------------------------------------------------------
 
     /** access to all boundary fluxes */
-    FIELDS &boundary_fluxes() {return boundaryFlux_;}; 
+    FIELDS &boundary_fluxes() {return boundaryFlux_;};
     /** wrapper for FE_Engine's compute_boundary_flux functions */
     void compute_boundary_flux(const Array2D<bool> & rhs_mask,
-                               const FIELDS &fields, 
+                               const FIELDS &fields,
                                FIELDS &rhs,
                                const Array< std::set <int> > atomMaterialGroups,
                                const VectorDependencyManager<SPAR_MAT * > * shpFcnDerivs,
@@ -128,12 +128,12 @@ namespace ATC {
 
     DENS_MAN &field_rhs(FieldName thisField) { return rhs_[thisField]; };
     /** allow FE_Engine to construct ATC structures after mesh is constructed */
-    virtual void initialize_mesh_data(void);  
+    virtual void initialize_mesh_data(void);
 
 // public for FieldIntegrator
-    bool source_atomic_quadrature(FieldName /* field */)  
+    bool source_atomic_quadrature(FieldName /* field */)
       { return (sourceIntegration_ == FULL_DOMAIN_ATOMIC_QUADRATURE_SOURCE); }
-    ATC::IntegrationDomainType source_integration() 
+    ATC::IntegrationDomainType source_integration()
       { return sourceIntegration_; }
 
     /** wrapper for FE_Engine's compute_sources */
@@ -143,25 +143,25 @@ namespace ATC {
                                   FIELD_MATS & atomicSources);
     /** computes tangent matrix using atomic quadrature near FE region */
    void masked_atom_domain_rhs_tangent(const std::pair<FieldName,FieldName> row_col,
-                                       const RHS_MASK & rhsMask,      
-                                       const FIELDS & fields,                  
+                                       const RHS_MASK & rhsMask,
+                                       const FIELDS & fields,
                                        SPAR_MAT & stiffness,
                                        const PhysicsModel * physicsModel);
     /** wrapper for FE_Engine's compute_rhs_vector functions */
     void compute_rhs_vector(const RHS_MASK & rhs_mask,
-                            const FIELDS &fields, 
+                            const FIELDS &fields,
                             FIELDS &rhs,
                             const IntegrationDomainType domain, // = FULL_DOMAIN
                             const PhysicsModel * physicsModel=nullptr);
    /** wrapper for FE_Engine's compute_tangent_matrix */
    void compute_rhs_tangent(const std::pair<FieldName,FieldName> row_col,
-                            const RHS_MASK & rhsMask,      
-                            const FIELDS & fields,                  
+                            const RHS_MASK & rhsMask,
+                            const FIELDS & fields,
                             SPAR_MAT & stiffness,
                             const IntegrationDomainType integrationType,
                             const PhysicsModel * physicsModel=nullptr);
    void tangent_matrix(const std::pair<FieldName,FieldName> row_col,
-                            const RHS_MASK & rhsMask,      
+                            const RHS_MASK & rhsMask,
                             const PhysicsModel * physicsModel,
                             SPAR_MAT & stiffness);
 
@@ -172,7 +172,7 @@ namespace ATC {
 
 // public for ImplicitSolveOperator
     /** return pointer to PrescribedDataManager */
-    PrescribedDataManager * prescribed_data_manager() 
+    PrescribedDataManager * prescribed_data_manager()
       { return prescribedDataMgr_; }
 // public for Kinetostat
     // TODO rename to "mass_matrix"
@@ -228,7 +228,7 @@ namespace ATC {
     void set_mass_mat_time_filter(FieldName thisField,TimeFilterManager::FilterIntegrationType filterIntegrationType);
 
     /** return reference to ExtrinsicModelManager */
-    ExtrinsicModelManager & extrinsic_model_manager() 
+    ExtrinsicModelManager & extrinsic_model_manager()
       { return extrinsicModelManager_; }
     /** access to time integrator */
     const TimeIntegrator * time_integrator(const FieldName & field) const {
@@ -242,7 +242,7 @@ namespace ATC {
     //---------------------------------------------------------------
     /*@{*/
     /** allow FE_Engine to construct data manager after mesh is constructed */
-    void construct_prescribed_data_manager (void); 
+    void construct_prescribed_data_manager (void);
     /** method to create physics model */
     void create_physics_model(const PhysicsType & physicsType,
                               std::string matFileName);
@@ -289,7 +289,7 @@ namespace ATC {
     void set_sources();
     /** assemble various contributions to the heat flux in the atomic region */
     void compute_atomic_sources(const Array2D<bool> & rhs_mask,
-                                const FIELDS &fields, 
+                                const FIELDS &fields,
                                 FIELDS &atomicSources);
 
     DENS_MAT &get_source(FieldName thisField){return sources_[thisField].set_quantity();};
@@ -313,25 +313,25 @@ namespace ATC {
     //---------------------------------------------------------------
     /*@{*/
     /** access for field mask */
-    Array2D<bool> &field_mask() {return fieldMask_;}; 
+    Array2D<bool> &field_mask() {return fieldMask_;};
     /** create field mask */
     void reset_flux_mask();
     /** field mask for intrinsic integration */
     Array2D<bool> intrinsicMask_;
     /** wrapper for FE_Engine's compute_flux functions */
     void compute_flux(const Array2D<bool> & rhs_mask,
-                      const FIELDS &fields, 
+                      const FIELDS &fields,
                       GRAD_FIELD_MATS &flux,
                       const PhysicsModel * physicsModel=nullptr,
                       const bool normalize = false);
     /** evaluate rhs on the atomic domain which is near the FE region */
     void masked_atom_domain_rhs_integral(const Array2D<bool> & rhs_mask,
-                                         const FIELDS &fields, 
+                                         const FIELDS &fields,
                                          FIELDS &rhs,
                                          const PhysicsModel * physicsModel);
     /** evaluate rhs on a specified domain defined by mask and physics model */
     void evaluate_rhs_integral(const Array2D<bool> & rhs_mask,
-                           const FIELDS &fields, 
+                           const FIELDS &fields,
                            FIELDS &rhs,
                            const IntegrationDomainType domain,
                            const PhysicsModel * physicsModel=nullptr);
@@ -422,8 +422,8 @@ namespace ATC {
     /*@{*/
     Array<int> elementToMaterialMap_;  // ATOMIC_TAG * elementToMaterialMap_;
     /** atomic ATC material tag */
-    
-    
+
+
     Array< std::set <int> > atomMaterialGroups_;  // ATOMIC_TAG*atomMaterialGroups_;
     Array< std::set <int> > atomMaterialGroupsMask_;  // ATOMIC_TAG*atomMaterialGroupsMask_;
     /*@}*/
@@ -432,7 +432,7 @@ namespace ATC {
     /** computational geometry */
     //---------------------------------------------------------------
     /*@{*/
-    bool atomQuadForInternal_;  
+    bool atomQuadForInternal_;
     MatrixDependencyManager<DenseMatrix, bool> * elementMask_;
     MatrixDependencyManager<DenseMatrix, bool> * elementMaskMass_;
     MatrixDependencyManager<DenseMatrix, bool> * elementMaskMassMd_;
@@ -465,7 +465,7 @@ namespace ATC {
     DIAG_MAN * atomicWeightsMask_;
     SPAR_MAN * shpFcnMask_;
     VectorDependencyManager<SPAR_MAT * > * shpFcnDerivsMask_;
-    Array<bool> atomMask_;  
+    Array<bool> atomMask_;
     SPAR_MAN atomToOverlapMat_;
     DIAG_MAN nodalMaskMat_;
     /*@}*/
@@ -475,7 +475,7 @@ namespace ATC {
     //---------------------------------------------------------------
     /*@{*/
     /** mask for computation of fluxes */
-    Array2D<bool> fieldMask_; 
+    Array2D<bool> fieldMask_;
     DIAG_MAT fluxMask_;
     DIAG_MAT fluxMaskComplement_;
     /** sources */

lib/atc/ATC_CouplingEnergy.cpp

@@ -20,7 +20,7 @@ namespace ATC {
   //  Class ATC_CouplingEnergy
   //--------------------------------------------------------
   //--------------------------------------------------------
-  
+
   //--------------------------------------------------------
   //  Constructor
   //--------------------------------------------------------
@@ -33,18 +33,18 @@ namespace ATC {
       nodalAtomicKineticTemperature_(nullptr),
       nodalAtomicConfigurationalTemperature_(nullptr)
   {
-    // Allocate PhysicsModel 
+    // Allocate PhysicsModel
     create_physics_model(THERMAL, matParamFile);
 
     // create extrinsic physics model
     if (extrinsicModel != NO_MODEL) {
-      extrinsicModelManager_.create_model(extrinsicModel,matParamFile);  
+      extrinsicModelManager_.create_model(extrinsicModel,matParamFile);
     }
 
     // Defaults
     set_time();
     bndyIntType_ = FE_INTERPOLATION;
-  
+
     // set up field data based on physicsModel
     physicsModel_->num_fields(fieldSizes_,fieldMask_);
 
@@ -182,16 +182,16 @@ namespace ATC {
     // register the per-atom quantity for the temperature definition
     interscaleManager_.add_per_atom_quantity(atomEnergyForTemperature,
                                              "AtomicEnergyForTemperature");
-    
+
     // nodal restriction of the atomic energy quantity for the temperature definition
     AtfShapeFunctionRestriction * nodalAtomicEnergy = new AtfShapeFunctionRestriction(this,
                                                                                       atomEnergyForTemperature,
                                                                                       shpFcn_);
     interscaleManager_.add_dense_matrix(nodalAtomicEnergy,
                                         "NodalAtomicEnergy");
-    
+
     // nodal atomic temperature field
-    
+
     AtfShapeFunctionMdProjection * nodalAtomicTemperature = new AtfShapeFunctionMdProjection(this,
                                                                                              nodalAtomicEnergy,
                                                                                              TEMPERATURE);
@@ -210,18 +210,18 @@ namespace ATC {
   //---------------------------------------------------------
   void ATC_CouplingEnergy::init_filter()
   {
-    
+
     TimeIntegrator::TimeIntegrationType timeIntegrationType = timeIntegrators_[TEMPERATURE]->time_integration_type();
-    
-    
-    
-    
-    if (timeFilterManager_.end_equilibrate()) { 
+
+
+
+
+    if (timeFilterManager_.end_equilibrate()) {
       if (timeIntegrationType==TimeIntegrator::GEAR) {
         if (equilibriumStart_) {
-          
-          
-          
+
+
+
           if (atomicRegulator_->regulator_target()==AtomicRegulator::DYNAMICS) { // based on FE equation
             DENS_MAT vdotflamMat(-2.*(nodalAtomicFields_[TEMPERATURE].quantity())); // note 2 is for 1/2 vdotflam addition
             atomicRegulator_->reset_lambda_contribution(vdotflamMat);
@@ -249,7 +249,7 @@ namespace ATC {
   {
     bool foundMatch = false;
     int argIndx = 0;
-    
+
     // check to see if input is a transfer class command
     // check derived class before base class
 
@@ -311,19 +311,19 @@ namespace ATC {
     // output[2] = average temperature
 
     double mvv2e = lammpsInterface_->mvv2e(); // convert to lammps energy units
-  
+
     if (n == 0) {
       Array<FieldName> mask(1);
       FIELD_MATS energy;
       mask(0) = TEMPERATURE;
-      
-      feEngine_->compute_energy(mask, 
+
+      feEngine_->compute_energy(mask,
                                 fields_,
                                 physicsModel_,
                                 elementToMaterialMap_,
                                 energy,
                                 &(elementMask_->quantity()));
-      
+
       double phononEnergy = mvv2e * energy[TEMPERATURE].col_sum();
       return phononEnergy;
     }
@@ -353,9 +353,9 @@ namespace ATC {
         _keTemp_ = nodalAtomicKineticTemperature_->quantity();
       if (nodalAtomicConfigurationalTemperature_)
       _peTemp_ = nodalAtomicConfigurationalTemperature_->quantity();
-      
+
       OUTPUT_LIST outputData;
-        
+
       // base class output
       ATC_Method::output();
 
@@ -370,7 +370,7 @@ namespace ATC {
       }
       atomicRegulator_->output(outputData);
       extrinsicModelManager_.output(outputData);
-      
+
       DENS_MAT & temperature(nodalAtomicFields_[TEMPERATURE].set_quantity());
       DENS_MAT & dotTemperature(dot_fields_[TEMPERATURE].set_quantity());
       DENS_MAT & ddotTemperature(ddot_fields_[TEMPERATURE].set_quantity());
@@ -384,18 +384,18 @@ namespace ATC {
         feEngine_->add_global("temperature_std_dev",  T_stddev);
         double Ta_mean =  (nodalAtomicFields_[TEMPERATURE].quantity()).col_sum(0)/nNodes_;
         feEngine_->add_global("atomic_temperature_mean",  Ta_mean);
-        double Ta_stddev =  (nodalAtomicFields_[TEMPERATURE].quantity()).col_stdev(0); 
+        double Ta_stddev =  (nodalAtomicFields_[TEMPERATURE].quantity()).col_stdev(0);
         feEngine_->add_global("atomic_temperature_std_dev",  Ta_stddev);
 
         // different temperature measures, if appropriate
         if (nodalAtomicKineticTemperature_)
           outputData["kinetic_temperature"] = & _keTemp_;
-        
+
         if (nodalAtomicConfigurationalTemperature_) {
           _peTemp_ *= 2; // account for full temperature
           outputData["configurational_temperature"] = & _peTemp_;
         }
-        
+
         // mesh data
         outputData["NodalAtomicTemperature"] = &temperature;
         outputData["dot_temperature"] = &dotTemperature;

lib/atc/ATC_CouplingEnergy.h

@@ -25,9 +25,9 @@ namespace ATC {
   class ATC_CouplingEnergy : public ATC_Coupling {
 
   public:
-  
+
     // constructor
-    ATC_CouplingEnergy(std::string groupName, 
+    ATC_CouplingEnergy(std::string groupName,
                        double ** & perAtomArray,
                        LAMMPS_NS::Fix * thisFix,
                        std::string matParamFile,
@@ -64,7 +64,7 @@ namespace ATC {
     /** physics specific filter initialization */
     void init_filter();
 
-    
+
     double compute_lambda_power(int gid);
 
     /** kinetic temperature for post-processing */
@@ -77,6 +77,6 @@ namespace ATC {
     DENS_MAT _keTemp_, _peTemp_;
 
   };
-    
+
 };
 #endif

lib/atc/ATC_CouplingMass.cpp

@@ -31,11 +31,11 @@ namespace ATC {
   //  Class ATC_CouplingMass
   //--------------------------------------------------------
   //--------------------------------------------------------
-  
+
   //--------------------------------------------------------
   //  Constructor
   //--------------------------------------------------------
-  ATC_CouplingMass::ATC_CouplingMass(string groupName, 
+  ATC_CouplingMass::ATC_CouplingMass(string groupName,
                                      double **& perAtomArray,
                                      LAMMPS_NS::Fix * thisFix,
                                      string matParamFile,
@@ -43,19 +43,19 @@ namespace ATC {
     : ATC_Coupling(groupName,perAtomArray,thisFix),
       resetNlocal_(false)
   {
-    // Allocate PhysicsModel 
-    create_physics_model(SPECIES, matParamFile); 
+    // Allocate PhysicsModel
+    create_physics_model(SPECIES, matParamFile);
 
     // create extrinsic physics model
     if (extrinsicModel != NO_MODEL) {
-      extrinsicModelManager_.create_model(extrinsicModel,matParamFile);  
+      extrinsicModelManager_.create_model(extrinsicModel,matParamFile);
     }
 
     // Defaults
     set_time();
     bndyIntType_ = NO_QUADRATURE;
-    
-  
+
+
     // set up field data based on physicsModel
     physicsModel_->num_fields(fieldSizes_,fieldMask_);
 
@@ -87,13 +87,13 @@ namespace ATC {
 
   //--------------------------------------------------------
   //  modify
-  //    parses inputs and modifies state 
+  //    parses inputs and modifies state
   //--------------------------------------------------------
   bool ATC_CouplingMass::modify(int narg, char **arg)
   {
     bool match = false;
     // check to see if it is a transfer class command
-    
+
     // check derived class before base class
     int argIndex = 0;
     // pass-through to concentration regulator
@@ -117,7 +117,7 @@ namespace ATC {
   //--------------------------------------------------------
   void ATC_CouplingMass::initialize()
   {
-    
+
     fieldSizes_[SPECIES_CONCENTRATION] = ntracked();
 
     // Base class initalizations
@@ -135,9 +135,9 @@ namespace ATC {
     FieldManager fmgr(this);
     atomicFields_[MASS_DENSITY]  = fmgr.nodal_atomic_field(MASS_DENSITY, field_to_intrinsic_name(MASS_DENSITY));
 
-    if (has_tracked_species()) { 
+    if (has_tracked_species()) {
       atomicFields_[SPECIES_CONCENTRATION]  = fmgr.nodal_atomic_field(SPECIES_CONCENTRATION, field_to_intrinsic_name(SPECIES_CONCENTRATION));
-      
+
       //if (atomicRegulator_->needs_temperature()) {
 
         atomicFields_[TEMPERATURE]  = fmgr.nodal_atomic_field(KINETIC_TEMPERATURE, field_to_intrinsic_name(TEMPERATURE));
@@ -146,7 +146,7 @@ namespace ATC {
         //}
     }
     else {
-      
+
       throw ATC_Error("ATC_CouplingMass: no tracked species");
     }
 
@@ -161,14 +161,14 @@ namespace ATC {
       const string moleculeName = molecule->first;
       SmallMoleculeSet * smallMoleculeSet = interscaleManager_.small_molecule_set(moleculeName);
       SPAR_MAN * shpFcnMol = interscaleManager_.sparse_matrix("ShapeFunction"+moleculeName);
-      AtomToSmallMoleculeTransfer<double> * moleculeMass = 
+      AtomToSmallMoleculeTransfer<double> * moleculeMass =
         new AtomToSmallMoleculeTransfer<double>(this,mass,smallMoleculeSet);
       interscaleManager_.add_dense_matrix(moleculeMass,"MoleculeMass"+moleculeName);
-      MotfShapeFunctionRestriction * nodalAtomicMoleculeMass = 
+      MotfShapeFunctionRestriction * nodalAtomicMoleculeMass =
         new MotfShapeFunctionRestriction(moleculeMass,shpFcnMol);
       interscaleManager_.add_dense_matrix(nodalAtomicMoleculeMass,"NodalMoleculeMass"+moleculeName);
 
-      
+
       AtfShapeFunctionMdProjection * nodalAtomicMoleculeMassDensity =
         new AtfShapeFunctionMdProjection(this,nodalAtomicMoleculeMass,MASS_DENSITY);
       interscaleManager_.add_dense_matrix(nodalAtomicMoleculeMassDensity,"NodalMoleculeMassDensity"+moleculeName);
@@ -180,7 +180,7 @@ namespace ATC {
 
   void ATC_CouplingMass::init_filter()
   {
-    
+
     ATC_Coupling::init_filter();
   }
 
@@ -192,11 +192,11 @@ namespace ATC {
   void ATC_CouplingMass::pre_exchange()
   {
     ATC_Coupling::pre_exchange();
-    
+
     //if (atomicRegulator_->needs_temperature()) {
-      field(TEMPERATURE) = atomicFields_[TEMPERATURE]->quantity(); 
+      field(TEMPERATURE) = atomicFields_[TEMPERATURE]->quantity();
 ///}
-    atomicRegulator_->pre_exchange(); 
+    atomicRegulator_->pre_exchange();
     if (resetNlocal_) {
       this->reset_nlocal();
       resetNlocal_ = false;
@@ -208,7 +208,7 @@ namespace ATC {
   //    does post-processing steps and outputs data
   //--------------------------------------------------------
   void ATC_CouplingMass::output()
-  { 
+  {
     if (output_now()) {
       feEngine_->departition_mesh();
       OUTPUT_LIST outputData;
@@ -226,10 +226,10 @@ namespace ATC {
         (_tiIt_->second)->output(outputData);
       }
       extrinsicModelManager_.output(outputData);
-      atomicRegulator_->output(outputData); 
+      atomicRegulator_->output(outputData);
 
       FIELD_POINTERS::iterator itr;
-      for (itr=atomicFields_.begin(); itr!=atomicFields_.end();itr++) { 
+      for (itr=atomicFields_.begin(); itr!=atomicFields_.end();itr++) {
         FieldName name = itr->first;
         const DENS_MAT & data = (itr->second)->quantity();
         outputData[field_to_intrinsic_name(name)] = & data;

lib/atc/ATC_CouplingMass.h

@@ -32,20 +32,20 @@ namespace ATC {
   class ATC_CouplingMass : public ATC_Coupling {
 
   public:
-  
+
     // constructor
-    ATC_CouplingMass(std::string groupName, 
+    ATC_CouplingMass(std::string groupName,
                      double **& perAtomArray,
                      LAMMPS_NS::Fix * thisFix,
                      std::string matParamFile,
                      ExtrinsicModelType extrinsic = NO_MODEL);
-      
+
     // destructor
     virtual ~ATC_CouplingMass();
 
     /** parser/modifier */
     virtual bool modify(int narg, char **arg);
-  
+
     /** pre time integration */
     virtual void initialize();
 
@@ -78,8 +78,8 @@ namespace ATC {
 
     bool resetNlocal_;
 
-    
-    
+
+
     //      i.e. we only need the correct shape function matrix for restriction
   };
 

lib/atc/ATC_CouplingMomentum.cpp

@@ -22,11 +22,11 @@ namespace ATC {
   //  Class ATC_CouplingMomentum
   //--------------------------------------------------------
   //--------------------------------------------------------
-  
+
   //--------------------------------------------------------
   //  Constructor
   //--------------------------------------------------------
-  ATC_CouplingMomentum::ATC_CouplingMomentum(string groupName, 
+  ATC_CouplingMomentum::ATC_CouplingMomentum(string groupName,
                                              double **& perAtomArray,
                                              LAMMPS_NS::Fix * thisFix,
                                              string matParamFile,
@@ -35,14 +35,14 @@ namespace ATC {
     : ATC_Coupling(groupName,perAtomArray,thisFix),
       refPE_(0)
   {
-    // Allocate PhysicsModel 
+    // Allocate PhysicsModel
     create_physics_model(intrinsicModel, matParamFile);
 
     // create extrinsic physics model
     if (extrinsicModel != NO_MODEL) {
-      extrinsicModelManager_.create_model(extrinsicModel,matParamFile);  
+      extrinsicModelManager_.create_model(extrinsicModel,matParamFile);
     }
-  
+
     // set up field data based on physicsModel
     physicsModel_->num_fields(fieldSizes_,fieldMask_);
 
@@ -166,14 +166,14 @@ namespace ATC {
         atomicMassWeightedDisplacement = new AtomicMassWeightedDisplacement(this);
       interscaleManager_.add_per_atom_quantity(atomicMassWeightedDisplacement,
                                                "AtomicMassWeightedDisplacement");
-      
+
       // nodal (RHS) mass-weighted displacement
       AtfShapeFunctionRestriction * nodalAtomicMassWeightedDisplacement = new AtfShapeFunctionRestriction(this,
                                                                                                  atomicMassWeightedDisplacement,
                                                                                                  shpFcn_);
       interscaleManager_.add_dense_matrix(nodalAtomicMassWeightedDisplacement,
                                           "NodalAtomicMassWeightedDisplacement");
-      
+
       // nodal displacement derived only from atoms
       AtfShapeFunctionMdProjection * nodalAtomicDisplacement = new AtfShapeFunctionMdProjection(this,
                                                                                                 nodalAtomicMassWeightedDisplacement,
@@ -194,7 +194,7 @@ namespace ATC {
   //---------------------------------------------------------
   void ATC_CouplingMomentum::init_filter()
   {
-    
+
     ATC_Coupling::init_filter();
 
     if (timeFilterManager_.end_equilibrate() && equilibriumStart_) // set up correct initial lambda forces to enforce initial accerlation
@@ -288,35 +288,35 @@ namespace ATC {
     ATC_Method::min_post_force();
 
     // Set sources
-    
+
     prescribedDataMgr_->set_sources(time(),sources_);
     extrinsicModelManager_.set_sources(fields_,extrinsicSources_);
     extrinsicModelManager_.pre_final_integrate();
 
 
 
-    
+
     if (outputNow_) {
       update_time(1.0);
       update_step();
       output();
       outputNow_ = false;
     }
-    
-    
+
+
     localStep_ += 1;
   }
-  
+
   //--------------------------------------------------------
   //  output
   //    does post-processing steps and outputs data
   //--------------------------------------------------------
-  void ATC_CouplingMomentum::output() 
+  void ATC_CouplingMomentum::output()
   {
     if (output_now()) {
       feEngine_->departition_mesh();
       OUTPUT_LIST outputData;
-      
+
       // base class output
       ATC_Method::output();
 
@@ -331,7 +331,7 @@ namespace ATC {
       }
       atomicRegulator_->output(outputData);
       extrinsicModelManager_.output(outputData);
-      
+
       DENS_MAT & velocity(nodalAtomicFields_[VELOCITY].set_quantity());
       DENS_MAT & rhs(rhs_[VELOCITY].set_quantity());
       if (lammpsInterface_->rank_zero()) {
@@ -341,7 +341,7 @@ namespace ATC {
         if (trackDisplacement_) {
           outputData["NodalAtomicDisplacement"] = & nodalAtomicFields_[DISPLACEMENT].set_quantity();
         }
-        
+
         feEngine_->write_data(output_index(), fields_, & outputData);
       }
       // force optional variables to reset to keep in sync
@@ -376,7 +376,7 @@ namespace ATC {
       kineticEnergy += v.dot(M*v);
     }
     if (domain == FE_DOMAIN) {
-      
+
       Array<FieldName> massMask(1);
       massMask(0) = VELOCITY;
       feEngine_->compute_lumped_mass_matrix(massMask,fields_,physicsModel_,atomMaterialGroups_,
@@ -388,20 +388,20 @@ namespace ATC {
         kineticEnergy -= v.dot(Ma*v);
       }
     }
-    double mvv2e = lammpsInterface_->mvv2e(); 
+    double mvv2e = lammpsInterface_->mvv2e();
     kineticEnergy *= 0.5*mvv2e; // convert to LAMMPS units
 
     return kineticEnergy;
   }
   //--------------------------------------------------------------------
-  //     potential/strain energy 
+  //     potential/strain energy
   //--------------------------------------------------------------------
   double ATC_CouplingMomentum::potential_energy(const IntegrationDomainType domain) const
   {
     Array<FieldName> mask(1);
     mask(0) = VELOCITY;
     FIELD_MATS energy;
-    feEngine_->compute_energy(mask, 
+    feEngine_->compute_energy(mask,
                                 fields_,
                                 physicsModel_,
                                 elementToMaterialMap_,
@@ -424,7 +424,7 @@ namespace ATC {
     // output[3] = total coarse scale energy
     // output[4] = fe-only coarse scale kinetic energy
     // output[5] = fe-only coarse scale potential energy
-  
+
 
 
     if (n == 0) {

lib/atc/ATC_CouplingMomentum.h

@@ -26,21 +26,21 @@ namespace ATC {
   class ATC_CouplingMomentum : public ATC_Coupling {
 
   public:
-  
+
     // constructor
-    ATC_CouplingMomentum(std::string groupName, 
+    ATC_CouplingMomentum(std::string groupName,
                          double **& perAtomArray,
                          LAMMPS_NS::Fix * thisFix,
                          std::string matParamFile,
                          PhysicsType intrinsicModel,
                          ExtrinsicModelType extrinsicModel = NO_MODEL);
-      
+
     // destructor
     virtual ~ATC_CouplingMomentum();
 
     /** parser/modifier */
     virtual bool modify(int narg, char **arg);
-  
+
     /** pre time integration */
     virtual void initialize();
 

lib/atc/ATC_CouplingMomentumEnergy.cpp

@@ -22,7 +22,7 @@ namespace ATC {
   //  Class ATC_CouplingMomentumEnergy
   //--------------------------------------------------------
   //--------------------------------------------------------
-  
+
   //--------------------------------------------------------
   //  Constructor
   //--------------------------------------------------------
@@ -36,12 +36,12 @@ namespace ATC {
       nodalAtomicConfigurationalTemperature_(nullptr),
       refPE_(0)
   {
-    // Allocate PhysicsModel 
+    // Allocate PhysicsModel
     create_physics_model(THERMO_ELASTIC, matParamFile);
 
     // create extrinsic physics model
     if (extrinsicModel != NO_MODEL) {
-      extrinsicModelManager_.create_model(extrinsicModel,matParamFile);  
+      extrinsicModelManager_.create_model(extrinsicModel,matParamFile);
     }
 
     // set up field data based on physicsModel
@@ -129,14 +129,14 @@ namespace ATC {
     AtomicMomentum * atomicMomentum = new AtomicMomentum(this);
     interscaleManager_.add_per_atom_quantity(atomicMomentum,
                                              "AtomicMomentum");
-    
+
     // nodal momentum for RHS
     AtfShapeFunctionRestriction * nodalAtomicMomentum = new AtfShapeFunctionRestriction(this,
                                                                                         atomicMomentum,
                                                                                         shpFcn_);
     interscaleManager_.add_dense_matrix(nodalAtomicMomentum,
                                         "NodalAtomicMomentum");
-    
+
     // nodal forces
     FundamentalAtomQuantity * atomicForce = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_FORCE);
     AtfShapeFunctionRestriction * nodalAtomicForce = new AtfShapeFunctionRestriction(this,
@@ -144,14 +144,14 @@ namespace ATC {
                                                                                      shpFcn_);
     interscaleManager_.add_dense_matrix(nodalAtomicForce,
                                         "NodalAtomicForce");
-    
+
     // nodal velocity derived only from atoms
     AtfShapeFunctionMdProjection * nodalAtomicVelocity = new AtfShapeFunctionMdProjection(this,
                                                                                           nodalAtomicMomentum,
                                                                                           VELOCITY);
     interscaleManager_.add_dense_matrix(nodalAtomicVelocity,
                                         "NodalAtomicVelocity");
-    
+
     if (trackDisplacement_) {
       // mass-weighted (center-of-mass) displacement of each atom
       AtomicMassWeightedDisplacement * atomicMassWeightedDisplacement;
@@ -166,14 +166,14 @@ namespace ATC {
         atomicMassWeightedDisplacement = new AtomicMassWeightedDisplacement(this);
       interscaleManager_.add_per_atom_quantity(atomicMassWeightedDisplacement,
                                                "AtomicMassWeightedDisplacement");
-      
+
       // nodal (RHS) mass-weighted displacement
       AtfShapeFunctionRestriction * nodalAtomicMassWeightedDisplacement = new AtfShapeFunctionRestriction(this,
                                                                                                           atomicMassWeightedDisplacement,
                                                                                                           shpFcn_);
       interscaleManager_.add_dense_matrix(nodalAtomicMassWeightedDisplacement,
                                           "NodalAtomicMassWeightedDisplacement");
-      
+
       // nodal displacement derived only from atoms
       AtfShapeFunctionMdProjection * nodalAtomicDisplacement = new AtfShapeFunctionMdProjection(this,
                                                                                                 nodalAtomicMassWeightedDisplacement,
@@ -183,7 +183,7 @@ namespace ATC {
     }
 
     // always need fluctuating velocity and kinetic energy
-    
+
     FtaShapeFunctionProlongation * atomicMeanVelocity = new FtaShapeFunctionProlongation(this,&fields_[VELOCITY],shpFcn_);
     interscaleManager_.add_per_atom_quantity(atomicMeanVelocity,
                                              field_to_prolongation_name(VELOCITY));
@@ -267,16 +267,16 @@ namespace ATC {
     // register the per-atom quantity for the temperature definition
     interscaleManager_.add_per_atom_quantity(atomEnergyForTemperature,
                                              "AtomicEnergyForTemperature");
-    
+
     // nodal restriction of the atomic energy quantity for the temperature definition
     AtfShapeFunctionRestriction * nodalAtomicEnergy = new AtfShapeFunctionRestriction(this,
                                                                                       atomEnergyForTemperature,
                                                                                       shpFcn_);
     interscaleManager_.add_dense_matrix(nodalAtomicEnergy,
                                         "NodalAtomicEnergy");
-    
+
     // nodal atomic temperature field
-    
+
     AtfShapeFunctionMdProjection * nodalAtomicTemperature = new AtfShapeFunctionMdProjection(this,
                                                                                              nodalAtomicEnergy,
                                                                                              TEMPERATURE);
@@ -299,17 +299,17 @@ namespace ATC {
       throw ATC_Error("ATC_CouplingMomentumEnergy::initialize - method only valid with fractional step time integration");
     }
 
-    
+
     ATC_Coupling::init_filter();
-    
-    
-    
-    
+
+
+
+
     if (timeFilterManager_.end_equilibrate() && equilibriumStart_) {
       if (atomicRegulator_->coupling_mode(VELOCITY)==AtomicRegulator::FLUX || atomicRegulator_->coupling_mode(VELOCITY)==AtomicRegulator::GHOST_FLUX)
         // nothing needed in other cases since kinetostat force is balanced by boundary flux in FE equations
        atomicRegulator_->reset_lambda_contribution(nodalAtomicFieldsRoc_[VELOCITY].quantity(),VELOCITY);
- 
+
       DENS_MAT powerMat(-1.*(nodalAtomicFields_[TEMPERATURE].quantity()));
       atomicRegulator_->reset_lambda_contribution(powerMat,TEMPERATURE);
     }
@@ -360,7 +360,7 @@ namespace ATC {
     Array<FieldName> mask(1);
     mask(0) = VELOCITY;
     FIELD_MATS energy;
-    feEngine_->compute_energy(mask, 
+    feEngine_->compute_energy(mask,
                                 fields_,
                                 physicsModel_,
                                 elementToMaterialMap_,
@@ -371,7 +371,7 @@ namespace ATC {
     potentialEnergy *= mvv2e; // convert to LAMMPS units
     return potentialEnergy-refPE_;
   }
-  
+
   //--------------------------------------------------------------------
   //     compute_vector
   //--------------------------------------------------------------------
@@ -399,14 +399,14 @@ namespace ATC {
       Array<FieldName> mask(1);
       FIELD_MATS energy;
       mask(0) = TEMPERATURE;
-      
-      feEngine_->compute_energy(mask, 
+
+      feEngine_->compute_energy(mask,
                                 fields_,
                                 physicsModel_,
                                 elementToMaterialMap_,
                                 energy,
                                 &(elementMask_->quantity()));
-      
+
       double phononEnergy = mvv2e * energy[TEMPERATURE].col_sum();
       return phononEnergy;
     }
@@ -435,9 +435,9 @@ namespace ATC {
         _keTemp_ = nodalAtomicKineticTemperature_->quantity();
       if (nodalAtomicConfigurationalTemperature_)
         _peTemp_ = nodalAtomicConfigurationalTemperature_->quantity();
-      
+
       OUTPUT_LIST outputData;
-      
+
       // base class output
       ATC_Method::output();
 
@@ -447,13 +447,13 @@ namespace ATC {
       }
 
       // auxiliary data
-      
+
       for (_tiIt_ = timeIntegrators_.begin(); _tiIt_ != timeIntegrators_.end(); ++_tiIt_) {
         (_tiIt_->second)->output(outputData);
       }
       atomicRegulator_->output(outputData);
       extrinsicModelManager_.output(outputData);
-        
+
       DENS_MAT & velocity(nodalAtomicFields_[VELOCITY].set_quantity());
       DENS_MAT & rhs(rhs_[VELOCITY].set_quantity());
       DENS_MAT & temperature(nodalAtomicFields_[TEMPERATURE].set_quantity());
@@ -469,18 +469,18 @@ namespace ATC {
         feEngine_->add_global("temperature_std_dev",  T_stddev);
         double Ta_mean =  (nodalAtomicFields_[TEMPERATURE].quantity()).col_sum(0)/nNodes_;
         feEngine_->add_global("atomic_temperature_mean",  Ta_mean);
-        double Ta_stddev =  (nodalAtomicFields_[TEMPERATURE].quantity()).col_stdev(0); 
+        double Ta_stddev =  (nodalAtomicFields_[TEMPERATURE].quantity()).col_stdev(0);
         feEngine_->add_global("atomic_temperature_std_dev",  Ta_stddev);
 
         // different temperature measures, if appropriate
         if (nodalAtomicKineticTemperature_)
           outputData["kinetic_temperature"] = & _keTemp_;
-        
+
         if (nodalAtomicConfigurationalTemperature_) {
           _peTemp_ *= 2; // account for full temperature
           outputData["configurational_temperature"] = & _peTemp_;
         }
-        
+
         // mesh data
         outputData["NodalAtomicVelocity"] = &velocity;
         outputData["FE_Force"] = &rhs;
@@ -491,10 +491,10 @@ namespace ATC {
         outputData["ddot_temperature"] = &ddotTemperature;
         outputData["NodalAtomicPower"] = &rocTemperature;
         outputData["fePower"] = &fePower;
-      
+
         feEngine_->write_data(output_index(), fields_, & outputData);
       }
-      
+
       //      hence propagation is performed on proc 0 but not others.
       //      The real fix is to have const data in the output list
       // force optional variables to reset to keep in sync

lib/atc/ATC_CouplingMomentumEnergy.h

@@ -25,9 +25,9 @@ namespace ATC {
   class ATC_CouplingMomentumEnergy : public ATC_Coupling {
 
   public:
-  
+
     // constructor
-    ATC_CouplingMomentumEnergy(std::string groupName, 
+    ATC_CouplingMomentumEnergy(std::string groupName,
                                double ** & perAtomArray,
                                LAMMPS_NS::Fix * thisFix,
                                std::string matParamFile,
@@ -78,6 +78,6 @@ namespace ATC {
     // data
     double refPE_;
   };
-    
+
 };
 #endif

lib/atc/ATC_Method.h

@@ -66,12 +66,12 @@ namespace ATC {
 
     /** Predictor phase, executed before Verlet */
     virtual void pre_init_integrate() {
-      feEngine_->partition_mesh(); 
+      feEngine_->partition_mesh();
       update_step();
     };
 
     /** Predictor phase, Verlet first step for velocity and position */
-    virtual void init_integrate(); 
+    virtual void init_integrate();
 
     /** Predictor phase, executed after Verlet */
     virtual void post_init_integrate();
@@ -127,7 +127,7 @@ namespace ATC {
     /** access to interscale manager */
     InterscaleManager & interscale_manager() {return interscaleManager_;};
     /** access to lammps interface */
-    
+
     LammpsInterface const * lammps_interface() const {return lammpsInterface_;};
     /** access to time filter */
     TimeFilterManager * time_filter_manager() {return &timeFilterManager_;};
@@ -142,7 +142,7 @@ namespace ATC {
     /** compute vector for output */
     virtual double compute_vector(int /* n */) {return 0.;}
     /** compute vector for output */
-    virtual double compute_array(int /* irow */, int /* icol */) {return 0.;}; 
+    virtual double compute_array(int /* irow */, int /* icol */) {return 0.;};
     int scalar_flag() const {return scalarFlag_;}
     int vector_flag() const {return vectorFlag_;}
     int size_vector() const {return sizeVector_;}
@@ -164,16 +164,16 @@ namespace ATC {
 
     /** time/step functions */
     bool sample_now(void) const
-    { 
+    {
       int s = step();
       bool now =  ( (sampleFrequency_ > 0) && (s % sampleFrequency_ == 0));
       return now;
     }
     bool output_now(void) const
-    { 
+    {
       int s = step();
 
-      bool now = ( (outputFrequency_ > 0) && (s == 1 || s % outputFrequency_ == 0) ); 
+      bool now = ( (outputFrequency_ > 0) && (s == 1 || s % outputFrequency_ == 0) );
       now = now || outputNow_;
       return now;
     }
@@ -328,7 +328,7 @@ namespace ATC {
     const std::map<std::string,std::pair<MolSize,int> > & molecule_ids() const {return moleculeIds_;};
     /** access to internal element set */
     const std::string & internal_element_set() {return internalElementSet_;};
-    
+
 
     //----------------------------------------------------------------
     /** \name mass matrix operations */
@@ -339,7 +339,7 @@ namespace ATC {
     //      inverted using GMRES
     void apply_inverse_mass_matrix(MATRIX & data, FieldName thisField)
     {
-      
+
       if (useConsistentMassMatrix_(thisField)) {
         //data = consistentMassInverse_*data;
         data = (consistentMassMatsInv_[thisField].quantity())*data;
@@ -351,7 +351,7 @@ namespace ATC {
     void apply_inverse_mass_matrix(const MATRIX & data_in, MATRIX & data_out,
                                    FieldName thisField)
     {
-      if (useConsistentMassMatrix_(thisField)) { 
+      if (useConsistentMassMatrix_(thisField)) {
         //data_out = consistentMassInverse_*data_in;
         data_out = (consistentMassMatsInv_[thisField].quantity())*data_in;
         return;
@@ -365,7 +365,7 @@ namespace ATC {
 
     DIAG_MAN &mass_mat(FieldName thisField)
       { return massMats_[thisField];};
-    
+
 
     //---------------------------------------------------------------
     /** \name mass matrices  */
@@ -383,13 +383,13 @@ namespace ATC {
     void register_mass_matrix_dependency(DependencyManager * dependent,
                                          FieldName thisField)
     {
-      if (useConsistentMassMatrix_(thisField)) { 
+      if (useConsistentMassMatrix_(thisField)) {
         consistentMassMatsInv_[thisField].register_dependence(dependent);
         return;
       }
       massMatsInv_[thisField].register_dependence(dependent);
     };
-    
+
     void apply_inverse_md_mass_matrix(MATRIX & data, FieldName thisField)
     { data = (massMatsMdInv_[thisField].quantity())*data; };
     void register_md_mass_matrix_dependency(DependencyManager * dependent,
@@ -416,11 +416,11 @@ namespace ATC {
        return man->second;
      };
     /*@}*/
-                
+
     //----------------------------------------------------------------
     /** \name Interscale operators */
     //----------------------------------------------------------------
-    bool use_md_mass_normalization() const { return mdMassNormalization_;}  
+    bool use_md_mass_normalization() const { return mdMassNormalization_;}
     bool kernel_based() { return kernelBased_; }
     bool kernel_on_the_fly() const { return kernelOnTheFly_;}
     bool has_kernel_function() { return kernelFunction_ != nullptr; }
@@ -439,20 +439,20 @@ namespace ATC {
 
     double ke_scale() { return keScale_; }
     double pe_scale() { return peScale_; }
-    
+
     /** from a atom group, find the nodes that have non-zero shape function contributions */
 
     bool nodal_influence(const int groupbit, std::set<int>& nset, std::set<int>& aset, double tol =1.e-8);
-    int nodal_influence(const int groupbit, std::set<int>& nset, std::set<int>& aset, 
+    int nodal_influence(const int groupbit, std::set<int>& nset, std::set<int>& aset,
       bool ghost,
       double tol =1.e-8);
     /*@{*/
 
-    
+
 
  /** Restrict based on atomic volume integration for volumetric quantities : given w_\alpha, w_I = \sum_\alpha N_{I\alpha} w_\alpha */
     void restrict_volumetric_quantity(const MATRIX &atomData,
-                                      MATRIX &nodeData); 
+                                      MATRIX &nodeData);
     void restrict_volumetric_quantity(const MATRIX &atomData,
                                       MATRIX &nodeData,
                                       const SPAR_MAT & shpFcn);
@@ -474,7 +474,7 @@ namespace ATC {
   protected: /** methods */
     /** time functions */
     void set_time(double t=0) {simTime_=t;};
-    void update_time(double alpha = 1.0)  
+    void update_time(double alpha = 1.0)
     {
       double dt = lammpsInterface_->dt();
       simTime_ += alpha*dt;
@@ -506,7 +506,7 @@ namespace ATC {
 
     virtual void  read_restart_data(std::string fileName_, RESTART_LIST & data);
     virtual void write_restart_data(std::string fileName_, RESTART_LIST & data);
-    void pack_fields(RESTART_LIST & data); 
+    void pack_fields(RESTART_LIST & data);
 
    /** mass matrices */
     DIAG_MAT  massMatInv_;
@@ -564,11 +564,11 @@ namespace ATC {
 
     void reset_fields();
 
-  private: /** methods */ 
+  private: /** methods */
     ATC_Method(); // do not define
 
   protected: /** data */
-    
+
     /* parsed input requires changes */
     bool needReset_;
 
@@ -609,7 +609,7 @@ namespace ATC {
     PerAtomQuantity<double> * atomProcGhostCoarseGrainingPositions_;
     PerAtomQuantity<double> * atomReferencePositions_;
 
-    
+
     /** number of unique FE nodes */
     int nNodes_;
 
@@ -631,17 +631,17 @@ namespace ATC {
     bool trackDisplacement_;
 
     /** map from reference positions to element id, pointer is to internal only */
-    
+
     bool needsAtomToElementMap_;
     PerAtomQuantity<int> * atomElement_;
     PerAtomQuantity<int> * atomGhostElement_;
 
     /* use element sets to define internal and/or ghost regions */
     std::string internalElementSet_;
-    
+
     /** atomic ATC material tag */
-    
-    
+
+
     double Xprd_,Yprd_,Zprd_; // lengths of periodic box in reference frame
     double XY_,YZ_,XZ_;
     double boxXlo_,boxXhi_; // lo/hi bounds of periodic box in reference frame
@@ -671,22 +671,22 @@ namespace ATC {
     /** base name for output files */
     std::string outputPrefix_;
 
-    /** output flag */ 
-    
+    /** output flag */
+
     bool outputNow_;
     /** output time or step (for lammps compatibility) */
     bool outputTime_;
-  
+
     /** output frequency */
     int outputFrequency_;
-  
+
     /** sample frequency */
     int sampleFrequency_;
-  
+
     /** sample counter */
     int sampleCounter_;
 
-    TAG_FIELDS filteredData_; 
+    TAG_FIELDS filteredData_;
 
     double peScale_,keScale_;
 
@@ -702,7 +702,7 @@ namespace ATC {
     int sizeVector_;              // N = size of global vector
     int scalarVectorFreq_;        // frequency compute s/v data is available at
     int sizePerAtomCols_;         // N = size of per atom vector to dump
-    
+
     double **perAtomOutput_;      // per atom data
     double **&perAtomArray_;      // per atom data
     int extScalar_;               // 0/1 if scalar is intensive/extensive
@@ -724,15 +724,15 @@ namespace ATC {
     /** \name time integration and filtering fields */
     //---------------------------------------------------------------
     /*@{*/
-    
-    FIELDS dot_fields_; 
+
+    FIELDS dot_fields_;
     FIELDS ddot_fields_;
     FIELDS dddot_fields_;
 
     /** Restricted Fields */
 
     FIELDS nodalAtomicFields_;  // replaces fieldNdFiltered_
-    FIELDS nodalAtomicFieldsRoc_; 
+    FIELDS nodalAtomicFieldsRoc_;
 
     /*@}*/
 
@@ -740,9 +740,9 @@ namespace ATC {
     /** \name quadrature weights */
     //---------------------------------------------------------------
     /*@{*/
-    
-    DIAG_MAT NodeVolumes_; 
-    DIAG_MAN invNodeVolumes_; 
+
+    DIAG_MAT NodeVolumes_;
+    DIAG_MAN invNodeVolumes_;
     /** atomic quadrature integration weights (V_\alpha) */
     ProtectedAtomDiagonalMatrix<double> * atomVolume_;
     std::string atomicWeightsFile_;
@@ -770,7 +770,7 @@ namespace ATC {
     bool needProcGhost_;
     std::string groupTag_;
     std::string groupTagGhost_;
-  
+
     /** number of atoms of correct type,
         ghosts are atoms outside our domain of interest
         boundary are atoms contributing to boundary flux terms */
@@ -824,7 +824,7 @@ namespace ATC {
     SPAR_MAN  kernelAccumulantMol_; // KKM add
     SPAR_MAN  kernelAccumulantMolGrad_; // KKM add
     DIAG_MAN* accumulantWeights_;
-    DIAG_MAN* accumulantInverseVolumes_;  
+    DIAG_MAN* accumulantInverseVolumes_;
     int accumulantBandwidth_;
     /*@}*/
 
@@ -845,7 +845,7 @@ namespace ATC {
 
 
     //---------------------------------------------------------------
-    /** \name reference data */ 
+    /** \name reference data */
     //---------------------------------------------------------------
     bool hasRefPE_;
     bool setRefPE_;

lib/atc/ATC_Transfer.cpp

@@ -68,7 +68,7 @@ namespace ATC {
                              LAMMPS_NS::Fix * thisFix,
                              string matParamFile)
     : ATC_Method(groupName,perAtomArray,thisFix),
-      xPointer_(nullptr), 
+      xPointer_(nullptr),
       outputStepZero_(true),
       neighborReset_(false),
       pairMap_(nullptr),
@@ -95,7 +95,7 @@ namespace ATC {
       lmp->lattice(cb.cell_vectors, cb.basis_vectors);
       cb.inv_atom_volume = 1.0 / lmp->volume_per_atom();
       cb.e2mvv           = 1.0 / lmp->mvv2e();
-      cb.atom_mass       = lmp->atom_mass(1); 
+      cb.atom_mass       = lmp->atom_mass(1);
       cb.boltzmann       = lmp->boltz();
       cb.hbar            = lmp->hbar();
       cauchyBornStress_ = new StressCauchyBorn(fileId, cb);
@@ -103,7 +103,7 @@ namespace ATC {
 
     // Defaults
     set_time();
- 
+
     outputFlags_.reset(NUM_TOTAL_FIELDS);
     outputFlags_ = false;
     fieldFlags_.reset(NUM_TOTAL_FIELDS);
@@ -117,20 +117,20 @@ namespace ATC {
     for (int i = 0; i < NUM_TOTAL_FIELDS; i++) { outputFields_[i] = nullptr; }
 
     // Hardy requires ref positions for processor ghosts for bond list
-    
+
     //needXrefProcessorGhosts_ = true;
   }
 
   //-------------------------------------------------------------------
   ATC_Transfer::~ATC_Transfer()
   {
-    interscaleManager_.clear(); 
-    if (cauchyBornStress_) delete cauchyBornStress_; 
+    interscaleManager_.clear();
+    if (cauchyBornStress_) delete cauchyBornStress_;
   }
 
   //-------------------------------------------------------------------
   // called before the beginning of a "run"
-  void ATC_Transfer::initialize() 
+  void ATC_Transfer::initialize()
   {
     if (kernelOnTheFly_ && !readRefPE_ && !setRefPEvalue_) {
       if (setRefPE_) {
@@ -143,7 +143,7 @@ namespace ATC {
 
     ATC_Method::initialize();
 
-    if (!initialized_) { 
+    if (!initialized_) {
       if (cauchyBornStress_) cauchyBornStress_->initialize();
     }
 
@@ -164,16 +164,16 @@ namespace ATC {
     ghostManager_.initialize();
 
     // initialize bond matrix B_Iab
-    if ((! bondOnTheFly_) 
-       && ( ( fieldFlags_(STRESS) 
-           || fieldFlags_(ESHELBY_STRESS) 
+    if ((! bondOnTheFly_)
+       && ( ( fieldFlags_(STRESS)
+           || fieldFlags_(ESHELBY_STRESS)
            || fieldFlags_(HEAT_FLUX) ) ) ) {
       try {
-        compute_bond_matrix(); 
-      } 
-      catch(bad_alloc&) { 
+        compute_bond_matrix();
+      }
+      catch(bad_alloc&) {
         ATC::LammpsInterface::instance()->print_msg("stress/heat_flux will be computed on-the-fly");
-        
+
         bondOnTheFly_ = true;
       }
     }
@@ -181,7 +181,7 @@ namespace ATC {
     // set sample frequency to output if sample has not be specified
     if (sampleFrequency_ == 0) sampleFrequency_ = outputFrequency_;
 
-    // output for step 0 
+    // output for step 0
     if (!initialized_) {
       if (outputFrequency_ > 0) {
         // initialize filtered data
@@ -224,9 +224,9 @@ namespace ATC {
 
     lammpsInterface_->computes_addstep(lammpsInterface_->ntimestep()+sampleFrequency_);
 
-    
+
     //remap_ghost_ref_positions();
-    update_peratom_output(); 
+    update_peratom_output();
   }
 
   //-------------------------------------------------------------------
@@ -242,7 +242,7 @@ namespace ATC {
   void ATC_Transfer::construct_time_integration_data()
   {
     if (!initialized_) {
-      
+
       // size arrays for requested/required  fields
       for(int index=0; index < NUM_TOTAL_FIELDS; ++index) {
         if (fieldFlags_(index)) {
@@ -298,13 +298,13 @@ namespace ATC {
   {
     // interpolant
     if (!(kernelOnTheFly_)) {
-      // finite element shape functions for interpolants 
+      // finite element shape functions for interpolants
       PerAtomShapeFunction * atomShapeFunctions = new PerAtomShapeFunction(this);
       interscaleManager_.add_per_atom_sparse_matrix(atomShapeFunctions,"Interpolant");
       shpFcn_ = atomShapeFunctions;
     }
     // accummulant and weights
-    
+
     this->create_atom_volume();
     // accumulants
     if (kernelFunction_) {
@@ -312,7 +312,7 @@ namespace ATC {
       if (kernelOnTheFly_) {
         ConstantQuantity<double> * atomCount = new ConstantQuantity<double>(this,1.);
         interscaleManager_.add_per_atom_quantity(atomCount,"AtomCount");
-        OnTheFlyKernelAccumulation * myWeights 
+        OnTheFlyKernelAccumulation * myWeights
            = new OnTheFlyKernelAccumulation(this,
              atomCount, kernelFunction_, atomCoarseGrainingPositions_);
         interscaleManager_.add_dense_matrix(myWeights,
@@ -337,7 +337,7 @@ namespace ATC {
       if (kernelOnTheFly_) {
         ConstantQuantity<double> * atomCount = new ConstantQuantity<double>(this,1.);
         interscaleManager_.add_per_atom_quantity(atomCount,"AtomCount");
-        OnTheFlyMeshAccumulation * myWeights 
+        OnTheFlyMeshAccumulation * myWeights
            = new OnTheFlyMeshAccumulation(this,
              atomCount, atomCoarseGrainingPositions_);
         interscaleManager_.add_dense_matrix(myWeights,
@@ -363,13 +363,13 @@ namespace ATC {
     // molecule centroid, molecule charge, dipole moment and quadrupole moment calculations KKM add
     if (!moleculeIds_.empty()) {
       map<string,pair<MolSize,int> >::const_iterator molecule;
-      InterscaleManager & interscaleManager = this->interscale_manager(); // KKM add, may be we do not need this as interscaleManager_ already exists.  
+      InterscaleManager & interscaleManager = this->interscale_manager(); // KKM add, may be we do not need this as interscaleManager_ already exists.
       PerAtomQuantity<double> * atomProcGhostCoarseGrainingPositions_ = interscaleManager.per_atom_quantity("AtomicProcGhostCoarseGrainingPositions");
       FundamentalAtomQuantity * mass = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_MASS,PROC_GHOST);
       molecule = moleculeIds_.begin();
       int groupbit = (molecule->second).second;
       smallMoleculeSet_ = new SmallMoleculeSet(this,groupbit);
-      smallMoleculeSet_->initialize(); // KKM add, why should we? 
+      smallMoleculeSet_->initialize(); // KKM add, why should we?
       interscaleManager_.add_small_molecule_set(smallMoleculeSet_,"MoleculeSet");
       moleculeCentroid_ = new SmallMoleculeCentroid(this,mass,smallMoleculeSet_,atomProcGhostCoarseGrainingPositions_);
       interscaleManager_.add_dense_matrix(moleculeCentroid_,"MoleculeCentroid");
@@ -393,7 +393,7 @@ namespace ATC {
 
     // set pointer to positions
     // REFACTOR use method's handling of xref/xpointer
-    set_xPointer(); 
+    set_xPointer();
 
     ATC_Method::construct_transfers();
 
@@ -412,7 +412,7 @@ namespace ATC {
     }
 
     // for hardy-based fluxes
-    
+
     bool needsBondMatrix =  (! bondOnTheFly_ ) &&
              (fieldFlags_(STRESS)
            || fieldFlags_(ESHELBY_STRESS)
@@ -434,7 +434,7 @@ namespace ATC {
 
       const FE_Mesh * fe_mesh = feEngine_->fe_mesh();
       if (!kernelBased_) {
-        bondMatrix_ = new BondMatrixPartitionOfUnity(lammpsInterface_,*pairMap_,xPointer_,fe_mesh,accumulantInverseVolumes_); 
+        bondMatrix_ = new BondMatrixPartitionOfUnity(lammpsInterface_,*pairMap_,xPointer_,fe_mesh,accumulantInverseVolumes_);
       }
       else {
         bondMatrix_ = new BondMatrixKernel(lammpsInterface_,*pairMap_,xPointer_,fe_mesh,kernelFunction_);
@@ -470,7 +470,7 @@ namespace ATC {
 
     FieldManager fmgr(this);
 
-//  for(int index=0; index < NUM_TOTAL_FIELDS; ++index) 
+//  for(int index=0; index < NUM_TOTAL_FIELDS; ++index)
     for(int i=0; i < numFields_; ++i) {
       FieldName index = indices_[i];
       if (fieldFlags_(index)) {
@@ -492,9 +492,9 @@ namespace ATC {
       interscaleManager_.add_per_atom_quantity(c,tag);
       int projection = iter->second;
       DIAG_MAN * w = nullptr;
-      if      (projection == VOLUME_NORMALIZATION ) 
+      if      (projection == VOLUME_NORMALIZATION )
          { w = accumulantInverseVolumes_; }
-      else if (projection == NUMBER_NORMALIZATION ) 
+      else if (projection == NUMBER_NORMALIZATION )
          { w = accumulantWeights_; }
       if (kernelFunction_ && kernelOnTheFly_) {
         OnTheFlyKernelAccumulationNormalized * C = new OnTheFlyKernelAccumulationNormalized(this, c, kernelFunction_, atomCoarseGrainingPositions_, w);
@@ -507,7 +507,7 @@ namespace ATC {
         outputFieldsTagged_[tag] = C;
       }
     }
-  
+
   }
 
   //-------------------------------------------------------------------
@@ -519,18 +519,18 @@ namespace ATC {
     if ((!initialized_) || timeFilterManager_.need_reset()) {
       timeFilters_.reset(NUM_TOTAL_FIELDS+nComputes_);
       sampleCounter_ = 0;
-      
+
       // for filtered  fields
       for(int index=0; index < NUM_TOTAL_FIELDS; ++index) {
         if (fieldFlags_(index)) {
           string name = field_to_string((FieldName) index);
           filteredData_[name]  = 0.0;
-          timeFilters_(index) = timeFilterManager_.construct(); 
+          timeFilters_(index) = timeFilterManager_.construct();
         }
       }
-      
+
       // for filtered projected computes
-      
+
       //       lists/accessing of fields ( & computes)
       map <string,int>::const_iterator iter;
       int index = NUM_TOTAL_FIELDS;
@@ -546,7 +546,7 @@ namespace ATC {
 
   //-------------------------------------------------------------------
   // called after the end of a "run"
-  void ATC_Transfer::finish() 
+  void ATC_Transfer::finish()
   {
     // base class
     ATC_Method::finish();
@@ -560,7 +560,7 @@ namespace ATC {
 
     int argIdx = 0;
     // check to see if it is a transfer class command
-      /*! \page man_hardy_fields fix_modify AtC fields 
+      /*! \page man_hardy_fields fix_modify AtC fields
         \section syntax
         fix_modify AtC fields <all | none> \n
         fix_modify AtC fields <add | delete> <list_of_fields> \n
@@ -575,13 +575,13 @@ namespace ATC {
         temperature :  temperature derived from the relative atomic kinetic energy (as done by ) \n
         kinetic_temperature : temperature derived from the full kinetic energy  \n
         number_density : simple kernel estimation of number of atoms per unit volume \n
-        stress : 
+        stress :
         Cauchy stress tensor for eulerian analysis (atom_element_map), or
         1st Piola-Kirchhoff stress tensor for lagrangian analysis    \n
-        transformed_stress : 
-        1st Piola-Kirchhoff stress tensor for eulerian analysis (atom_element_map), or 
+        transformed_stress :
+        1st Piola-Kirchhoff stress tensor for eulerian analysis (atom_element_map), or
                              Cauchy stress tensor for lagrangian analysis    \n
-        heat_flux : spatial heat flux vector for eulerian, 
+        heat_flux : spatial heat flux vector for eulerian,
         or referential heat flux vector for lagrangian \n
         potential_energy : potential energy per unit volume \n
         kinetic_energy : kinetic energy per unit volume \n
@@ -590,23 +590,23 @@ namespace ATC {
         energy : total energy (potential + kinetic) per unit volume \n
         number_density : number of atoms per unit volume \n
         eshelby_stress: configurational stress (energy-momentum) tensor defined by Eshelby
-          [References: Philos. Trans. Royal Soc. London A, Math. Phys. Sci., Vol. 244, 
+          [References: Philos. Trans. Royal Soc. London A, Math. Phys. Sci., Vol. 244,
           No. 877 (1951) pp. 87-112; J. Elasticity, Vol. 5, Nos. 3-4 (1975) pp. 321-335] \n
         vacancy_concentration: volume fraction of vacancy content \n
         type_concentration: volume fraction of a specific atom type \n
         \section examples
         <TT> fix_modify AtC fields add velocity temperature </TT>
         \section description
-        Allows modification of the fields calculated and output by the 
+        Allows modification of the fields calculated and output by the
         transfer class. The commands are cumulative, e.g.\n
-        <TT> fix_modify AtC fields none </TT> \n 
-        followed by \n 
+        <TT> fix_modify AtC fields none </TT> \n
+        followed by \n
         <TT> fix_modify AtC fields add velocity temperature </TT> \n
         will only output the velocity and temperature fields.
         \section restrictions
         Must be used with the hardy/field type of AtC fix, see \ref man_fix_atc.
-        Currently, the stress and heat flux formulas are only correct for 
-        central force potentials, e.g. Lennard-Jones and EAM 
+        Currently, the stress and heat flux formulas are only correct for
+        central force potentials, e.g. Lennard-Jones and EAM
         but not Stillinger-Weber.
         \section related
         See \ref man_hardy_gradients , \ref man_hardy_rates  and \ref man_hardy_computes
@@ -615,30 +615,30 @@ namespace ATC {
       */
       if (strcmp(arg[argIdx],"fields")==0) {
         argIdx++;
-        if (strcmp(arg[argIdx],"all")==0) { 
+        if (strcmp(arg[argIdx],"all")==0) {
           outputFlags_ = true;
           match = true;
-        } 
-        else if (strcmp(arg[argIdx],"none")==0) { 
+        }
+        else if (strcmp(arg[argIdx],"none")==0) {
           outputFlags_ = false;
           match = true;
-        } 
-        else if (strcmp(arg[argIdx],"add")==0) { 
+        }
+        else if (strcmp(arg[argIdx],"add")==0) {
           argIdx++;
           for (int i = argIdx; i < narg; ++i) {
             FieldName field_name = string_to_field(arg[i]);
-            outputFlags_(field_name) = true; 
+            outputFlags_(field_name) = true;
           }
           match = true;
-        } 
-        else if (strcmp(arg[argIdx],"delete")==0) { 
+        }
+        else if (strcmp(arg[argIdx],"delete")==0) {
           argIdx++;
           for (int i = argIdx; i < narg; ++i) {
             FieldName field_name = string_to_field(arg[i]);
-            outputFlags_(field_name) = false; 
+            outputFlags_(field_name) = false;
           }
           match = true;
-        } 
+        }
         check_field_dependencies();
         if (fieldFlags_(DISPLACEMENT)) { trackDisplacement_ = true; }
       }
@@ -649,17 +649,17 @@ namespace ATC {
         - add | delete (keyword) = add or delete the calculation of gradients for the listed output fields \n
         - fields (keyword) =  \n
         gradients can be calculated for all fields listed in \ref man_hardy_fields
- 
+
         \section examples
         <TT> fix_modify AtC gradients add temperature velocity stress </TT> \n
         <TT> fix_modify AtC gradients delete velocity </TT> \n
         \section description
         Requests calculation and output of gradients of the fields from the
         transfer class. These gradients will be with regard to spatial or material
-        coordinate for eulerian or lagrangian analysis, respectively, as specified by 
+        coordinate for eulerian or lagrangian analysis, respectively, as specified by
         atom_element_map (see \ref man_atom_element_map )
         \section restrictions
-        Must be used with the hardy/field type of AtC fix 
+        Must be used with the hardy/field type of AtC fix
         ( see \ref man_fix_atc )
         \section related
         \section default
@@ -667,33 +667,33 @@ namespace ATC {
       */
       else if (strcmp(arg[argIdx],"gradients")==0) {
         argIdx++;
-        if (strcmp(arg[argIdx],"add")==0) { 
+        if (strcmp(arg[argIdx],"add")==0) {
           argIdx++;
           FieldName field_name;
           for (int i = argIdx; i < narg; ++i) {
             field_name = string_to_field(arg[i]);
-            gradFlags_(field_name) = true; 
+            gradFlags_(field_name) = true;
           }
           match = true;
-        } 
-        else if (strcmp(arg[argIdx],"delete")==0) { 
+        }
+        else if (strcmp(arg[argIdx],"delete")==0) {
           argIdx++;
           FieldName field_name;
           for (int i = argIdx; i < narg; ++i) {
             field_name = string_to_field(arg[i]);
-            gradFlags_(field_name) = false; 
+            gradFlags_(field_name) = false;
           }
           match = true;
-        } 
+        }
       }
 
-      /*! \page man_hardy_rates fix_modify AtC rates 
+      /*! \page man_hardy_rates fix_modify AtC rates
         \section syntax
         fix_modify AtC rates <add | delete> <list_of_fields> \n
         - add | delete (keyword) = add or delete the calculation of rates (time derivatives) for the listed output fields \n
         - fields (keyword) =  \n
         rates can be calculated for all fields listed in \ref man_hardy_fields
- 
+
         \section examples
         <TT> fix_modify AtC rates add temperature velocity stress </TT> \n
         <TT> fix_modify AtC rates delete stress </TT> \n
@@ -703,7 +703,7 @@ namespace ATC {
         are the partial time derivatives of the nodal fields, not the full (material) time
         derivatives. \n
         \section restrictions
-        Must be used with the hardy/field type of AtC fix 
+        Must be used with the hardy/field type of AtC fix
         ( see \ref man_fix_atc )
         \section related
         \section default
@@ -711,16 +711,16 @@ namespace ATC {
       */
       else if (strcmp(arg[argIdx],"rates")==0) {
         argIdx++;
-        if (strcmp(arg[argIdx],"add")==0) { 
+        if (strcmp(arg[argIdx],"add")==0) {
           argIdx++;
           FieldName field_name;
           for (int i = argIdx; i < narg; ++i) {
             field_name = string_to_field(arg[i]);
-            rateFlags_(field_name) = true; 
+            rateFlags_(field_name) = true;
           }
           match = true;
-        } 
-        else if (strcmp(arg[argIdx],"delete")==0) { 
+        }
+        else if (strcmp(arg[argIdx],"delete")==0) {
           argIdx++;
           FieldName field_name;
           for (int i = argIdx; i < narg; ++i) {
@@ -728,7 +728,7 @@ namespace ATC {
             rateFlags_(field_name) = false;
           }
           match = true;
-        } 
+        }
       }
 
 
@@ -736,7 +736,7 @@ namespace ATC {
         \section syntax
         fix_modify AtC pair_interactions <on|off> \n
         fix_modify AtC bond_interactions <on|off> \n
-      
+
         \section examples
         <TT> fix_modify AtC bond_interactions on </TT> \n
         \section description
@@ -748,27 +748,27 @@ namespace ATC {
       */
       if (strcmp(arg[argIdx],"pair_interactions")==0) {    // default true
         argIdx++;
-        if (strcmp(arg[argIdx],"on")==0) { hasPairs_ = true; }      
+        if (strcmp(arg[argIdx],"on")==0) { hasPairs_ = true; }
         else                             { hasPairs_ = false;}
         match = true;
-      }  
+      }
       if (strcmp(arg[argIdx],"bond_interactions")==0) {    // default false
         argIdx++;
-        if (strcmp(arg[argIdx],"on")==0) { hasBonds_ = true; }      
+        if (strcmp(arg[argIdx],"on")==0) { hasBonds_ = true; }
         else                             { hasBonds_ = false;}
         match = true;
-      }  
- 
-      /*! \page man_hardy_computes fix_modify AtC computes 
+      }
+
+      /*! \page man_hardy_computes fix_modify AtC computes
         \section syntax
         fix_modify AtC computes <add | delete> [per-atom compute id] <volume | number> \n
         - add | delete (keyword) = add or delete the calculation of an equivalent continuum field
         for the specified per-atom compute as volume or number density quantity \n
-        - per-atom compute id =  name/id for per-atom compute, 
+        - per-atom compute id =  name/id for per-atom compute,
         fields can be calculated for all per-atom computes available from LAMMPS \n
         - volume | number (keyword) = field created is a per-unit-volume quantity
-        or a per-atom quantity as weighted by kernel functions \n 
- 
+        or a per-atom quantity as weighted by kernel functions \n
+
         \section examples
         <TT> compute virial all stress/atom </TT> \n
         <TT> fix_modify AtC computes add virial volume </TT> \n
@@ -782,24 +782,24 @@ namespace ATC {
         Must be used with the hardy/field type of AtC fix ( see \ref man_fix_atc ) \n
         Per-atom compute must be specified before corresponding continuum field can be requested \n
         \section related
-        See manual page for compute 
+        See manual page for compute
         \section default
         No defaults exist for this command
       */
       else if (strcmp(arg[argIdx],"computes")==0) {
         argIdx++;
-        if (strcmp(arg[argIdx],"add")==0) { 
+        if (strcmp(arg[argIdx],"add")==0) {
           argIdx++;
           string tag(arg[argIdx++]);
           int normalization = NO_NORMALIZATION;
           if (narg > argIdx) {
-            if      (strcmp(arg[argIdx],"volume")==0) { 
+            if      (strcmp(arg[argIdx],"volume")==0) {
               normalization = VOLUME_NORMALIZATION;
             }
-            else if (strcmp(arg[argIdx],"number")==0) { 
+            else if (strcmp(arg[argIdx],"number")==0) {
               normalization = NUMBER_NORMALIZATION;
             }
-            else if (strcmp(arg[argIdx],"mass")==0) { 
+            else if (strcmp(arg[argIdx],"mass")==0) {
               normalization = MASS_NORMALIZATION;
               throw ATC_Error("mass normalized not implemented");
             }
@@ -807,8 +807,8 @@ namespace ATC {
           computes_[tag] = normalization;
           nComputes_++;
           match = true;
-        } 
-        else if (strcmp(arg[argIdx],"delete")==0) { 
+        }
+        else if (strcmp(arg[argIdx],"delete")==0) {
           argIdx++;
           string tag(arg[argIdx]);
           if (computes_.find(tag) != computes_.end()) {
@@ -816,10 +816,10 @@ namespace ATC {
             nComputes_--;
           }
           else {
-            throw ATC_Error(tag+" compute is not in list"); 
+            throw ATC_Error(tag+" compute is not in list");
           }
           match = true;
-        } 
+        }
       }
 
 
@@ -833,7 +833,7 @@ namespace ATC {
         Specifies a frequency at which fields are computed for the case
         where time filters are being applied.
         \section restrictions
-        Must be used with the hardy/field AtC fix ( see \ref man_fix_atc ) 
+        Must be used with the hardy/field AtC fix ( see \ref man_fix_atc )
         and is only relevant when time filters are being used.
         \section related
         \section default
@@ -869,11 +869,11 @@ namespace ATC {
   // REFACTOR move this to post_neighbor
   void ATC_Transfer::pre_final_integrate()
   {
-    // update time 
+    // update time
     update_time(); // time uses step if dt = 0
 
 
-    
+
     if ( neighborReset_ && sample_now() ) {
       if (! kernelOnTheFly_ ) {
         if (!moleculeIds_.empty()) compute_kernel_matrix_molecule(); //KKM add
@@ -889,7 +889,7 @@ namespace ATC {
     // compute spatially smoothed quantities
     double dt = lammpsInterface_->dt();
     if ( sample_now() ) {
-      
+
       bool needsBond =  (! bondOnTheFly_ ) &&
              (fieldFlags_(STRESS)
            || fieldFlags_(ESHELBY_STRESS)
@@ -898,7 +898,7 @@ namespace ATC {
       if ( needsBond ) {
         if (pairMap_->need_reset()) {
 //        ATC::LammpsInterface::instance()->print_msg("Recomputing bond matrix due to atomReset_ value");
-          compute_bond_matrix(); 
+          compute_bond_matrix();
         }
       }
       time_filter_pre (dt);
@@ -923,8 +923,8 @@ namespace ATC {
   //-------------------------------------------------------------------
   void ATC_Transfer::compute_fields(void)
   {
-    
-    // keep per-atom computes fresh. JAZ and REJ not sure why; 
+
+    // keep per-atom computes fresh. JAZ and REJ not sure why;
     // need to confer with JAT. (JAZ, 4/5/12)
     interscaleManager_.lammps_force_reset();
 
@@ -937,16 +937,16 @@ namespace ATC {
       }
     }
 
-    if (fieldFlags_(STRESS)) 
+    if (fieldFlags_(STRESS))
       compute_stress(hardyData_["stress"].set_quantity());
-    if (fieldFlags_(HEAT_FLUX)) 
+    if (fieldFlags_(HEAT_FLUX))
       compute_heatflux(hardyData_["heat_flux"].set_quantity());
 // molecule data
     if (fieldFlags_(DIPOLE_MOMENT))
-      compute_dipole_moment(hardyData_["dipole_moment"].set_quantity()); 
+      compute_dipole_moment(hardyData_["dipole_moment"].set_quantity());
     if (fieldFlags_(QUADRUPOLE_MOMENT))
       compute_quadrupole_moment(hardyData_["quadrupole_moment"].set_quantity());
-    if (fieldFlags_(DISLOCATION_DENSITY)) 
+    if (fieldFlags_(DISLOCATION_DENSITY))
       compute_dislocation_density(hardyData_["dislocation_density"].set_quantity());
 
     // (2) derived quantities
@@ -963,7 +963,7 @@ namespace ATC {
         }
       }
     }
-    // compute: eshelby stress 
+    // compute: eshelby stress
     if (fieldFlags_(ESHELBY_STRESS)) {
       {
       compute_eshelby_stress(hardyData_["eshelby_stress"].set_quantity(),
@@ -985,18 +985,18 @@ namespace ATC {
                              E,hardyData_["stress"].quantity(),
                              hardyData_["displacement_gradient"].quantity());
     }
-    // compute: cauchy born stress 
+    // compute: cauchy born stress
     if (fieldFlags_(CAUCHY_BORN_STRESS)) {
       ATOMIC_DATA::const_iterator tfield = hardyData_.find("temperature");
       const DENS_MAT *temp = tfield==hardyData_.end() ? nullptr : &((tfield->second).quantity());
       cauchy_born_stress(hardyData_["displacement_gradient"].quantity(),
                          hardyData_["cauchy_born_stress"].set_quantity(), temp);
     }
-    // compute: cauchy born energy 
+    // compute: cauchy born energy
     if (fieldFlags_(CAUCHY_BORN_ENERGY)) {
       ATOMIC_DATA::const_iterator tfield = hardyData_.find("temperature");
       const DENS_MAT *temp = tfield==hardyData_.end() ? nullptr : &((tfield->second).quantity());
-      cauchy_born_energy(hardyData_["displacement_gradient"].quantity(), 
+      cauchy_born_energy(hardyData_["displacement_gradient"].quantity(),
                          hardyData_["cauchy_born_energy"].set_quantity(), temp);
     }
     // 1st PK transformed to cauchy (lag) or cauchy transformed to 1st PK (eul)
@@ -1014,13 +1014,13 @@ namespace ATC {
       compute_electric_potential(
         hardyData_[field_to_string(ELECTRIC_POTENTIAL)].set_quantity());
     }
-    // compute: rotation and/or stretch from deformation gradient 
+    // compute: rotation and/or stretch from deformation gradient
     if (fieldFlags_(ROTATION) || fieldFlags_(STRETCH)) {
       compute_polar_decomposition(hardyData_["rotation"].set_quantity(),
                                   hardyData_["stretch"].set_quantity(),
                                   hardyData_["displacement_gradient"].quantity());
     }
-    // compute: rotation and/or stretch from deformation gradient 
+    // compute: rotation and/or stretch from deformation gradient
     if (fieldFlags_(CAUCHY_BORN_ELASTIC_DEFORMATION_GRADIENT)) {
       compute_elastic_deformation_gradient2(hardyData_["elastic_deformation_gradient"].set_quantity(),
                                            hardyData_["stress"].quantity(),
@@ -1082,9 +1082,9 @@ namespace ATC {
             F(0,0) += 1.0; F(1,1) += 1.0; F(2,2) += 1.0;
             FT = F.transpose();
             FTINV = inv(FT);
-            
+
             //       volumes are already reference volumes.
-            PK1 = CAUCHY*FTINV; 
+            PK1 = CAUCHY*FTINV;
             matrix_to_vector(k,PK1,myData);
           }
         }
@@ -1094,7 +1094,7 @@ namespace ATC {
       }
 #endif
     }
-    
+
   }// end of compute_fields routine
 
   //-------------------------------------------------------------------
@@ -1156,7 +1156,7 @@ namespace ATC {
         filteredData_[grad_field] = hardyData_[grad_field];
       }
     }
-    
+
     //       lists/accessing of fields ( & computes)
     map <string,int>::const_iterator iter;
     int index = NUM_TOTAL_FIELDS;
@@ -1177,7 +1177,7 @@ namespace ATC {
   void ATC_Transfer::output()
   {
     feEngine_->departition_mesh();
-    
+
     for(int index=0; index < NUM_TOTAL_FIELDS; ++index) {
       if (outputFlags_(index)) {
         FieldName fName = (FieldName) index;
@@ -1185,7 +1185,7 @@ namespace ATC {
         fields_[fName] = filteredData_[name];
       }
     }
-    
+
     ATC_Method::output();
     if (lammpsInterface_->comm_rank() == 0) {
       // data
@@ -1209,7 +1209,7 @@ namespace ATC {
           output_data[grad_name] = & ( filteredData_[grad_name].set_quantity());
         }
       }
-      
+
       //       lists/accessing of fields ( & computes)
       map <string,int>::const_iterator iter;
       for (iter = computes_.begin(); iter != computes_.end(); iter++) {
@@ -1226,7 +1226,7 @@ namespace ATC {
       output_data["NodalInverseVolumes"] = &nodalInverseVolumes;
 
       // output
-      feEngine_->write_data(output_index(), & output_data); 
+      feEngine_->write_data(output_index(), & output_data);
     }
     feEngine_->partition_mesh();
   }
@@ -1235,7 +1235,7 @@ namespace ATC {
 
 /////// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   //-------------------------------------------------------------------
-  // computes nodeData = N*atomData 
+  // computes nodeData = N*atomData
   void ATC_Transfer::project(const DENS_MAT & atomData,
                                   DENS_MAT & nodeData)
   {
@@ -1290,8 +1290,8 @@ namespace ATC {
   void ATC_Transfer::project_count_normalized(const DENS_MAT & atomData,
                                   DENS_MAT & nodeData)
   {
-    DENS_MAT tmp; 
-    project(atomData,tmp); 
+    DENS_MAT tmp;
+    project(atomData,tmp);
     nodeData = (accumulantWeights_->quantity())*tmp;
   }
 
@@ -1301,7 +1301,7 @@ namespace ATC {
                                         DENS_MAT & nodeData)
   {
     DENS_MAT tmp;
-    project(atomData,tmp); 
+    project(atomData,tmp);
     nodeData = (accumulantInverseVolumes_->quantity())*tmp;
   }
 
@@ -1310,8 +1310,8 @@ namespace ATC {
   void ATC_Transfer::project_volume_normalized_molecule(const DENS_MAT & molData,
                                         DENS_MAT & nodeData)
   {
-    DENS_MAT tmp; 
-    project_molecule(molData,tmp); 
+    DENS_MAT tmp;
+    project_molecule(molData,tmp);
     nodeData = (accumulantInverseVolumes_->quantity())*tmp;
   }
 
@@ -1320,8 +1320,8 @@ namespace ATC {
   void ATC_Transfer::project_volume_normalized_molecule_gradient(const DENS_MAT & molData,
                                         DENS_MAT & nodeData)
   {
-    DENS_MAT tmp; 
-    project_molecule_gradient(molData,tmp); 
+    DENS_MAT tmp;
+    project_molecule_gradient(molData,tmp);
     nodeData = (accumulantInverseVolumes_->quantity())*tmp;
   }
 
@@ -1354,14 +1354,14 @@ namespace ATC {
 
   //-------------------------------------------------------------------
   // computes "virial" part of heat flux
-  // This is correct ONLY for pair potentials. 
+  // This is correct ONLY for pair potentials.
   void ATC_Transfer::compute_heat_matrix()
   {
     atomicHeatMatrix_ = pairHeatFlux_->quantity();
   }
 
   //-------------------------------------------------------------------
-  // set xPointer_ to xref or xatom depending on Lagrangian/Eulerian analysis 
+  // set xPointer_ to xref or xatom depending on Lagrangian/Eulerian analysis
   void ATC_Transfer::set_xPointer()
   {
     xPointer_ = xref_;
@@ -1386,7 +1386,7 @@ namespace ATC {
       fieldFlags_(DISPLACEMENT) = true;
     }
     if (fieldFlags_(CAUCHY_BORN_STRESS)
-        || fieldFlags_(CAUCHY_BORN_ENERGY) 
+        || fieldFlags_(CAUCHY_BORN_ENERGY)
         || fieldFlags_(CAUCHY_BORN_ESHELBY_STRESS)
         || fieldFlags_(CAUCHY_BORN_ELASTIC_DEFORMATION_GRADIENT))  {
       if (! (cauchyBornStress_) ) {
@@ -1419,7 +1419,7 @@ namespace ATC {
       fieldFlags_(KINETIC_ENERGY) = true;
     }
     if (fieldFlags_(TEMPERATURE) || fieldFlags_(HEAT_FLUX) ||
-        fieldFlags_(KINETIC_ENERGY) || fieldFlags_(THERMAL_ENERGY) || 
+        fieldFlags_(KINETIC_ENERGY) || fieldFlags_(THERMAL_ENERGY) ||
         fieldFlags_(ENERGY) || fieldFlags_(INTERNAL_ENERGY) ||
         fieldFlags_(KINETIC_ENERGY) || (fieldFlags_(STRESS) &&
         atomToElementMapType_ == EULERIAN) ) {
@@ -1438,15 +1438,15 @@ namespace ATC {
       fieldFlags_(NUMBER_DENSITY) = true;
     }
 
-    if (fieldFlags_(ROTATION) || 
+    if (fieldFlags_(ROTATION) ||
         fieldFlags_(STRETCH)) {
       fieldFlags_(DISPLACEMENT) = true;
     }
     if (fieldFlags_(ESHELBY_STRESS)
-       || fieldFlags_(CAUCHY_BORN_STRESS) 
-       || fieldFlags_(CAUCHY_BORN_ENERGY) 
-       || fieldFlags_(CAUCHY_BORN_ESHELBY_STRESS) 
-       || fieldFlags_(CAUCHY_BORN_ELASTIC_DEFORMATION_GRADIENT) 
+       || fieldFlags_(CAUCHY_BORN_STRESS)
+       || fieldFlags_(CAUCHY_BORN_ENERGY)
+       || fieldFlags_(CAUCHY_BORN_ESHELBY_STRESS)
+       || fieldFlags_(CAUCHY_BORN_ELASTIC_DEFORMATION_GRADIENT)
        || fieldFlags_(VACANCY_CONCENTRATION)
        || fieldFlags_(ROTATION)
        || fieldFlags_(STRETCH) ) {
@@ -1459,8 +1459,8 @@ namespace ATC {
         throw ATC_Error("Calculation of  stress field not possible with selected pair type.");
       }
     }
-    
-  } 
+
+  }
 
 //============== THIN WRAPPERS ====================================
 // OBSOLETE
@@ -1492,7 +1492,7 @@ namespace ATC {
 
     // calculate kinetic energy tensor part of stress for Eulerian analysis
     if (atomToElementMapType_ == EULERIAN && nLocal_>0) {
-      compute_kinetic_stress(stress); 
+      compute_kinetic_stress(stress);
     }
     else {
       // zero stress table for Lagrangian analysis or if nLocal_ = 0
@@ -1511,7 +1511,7 @@ namespace ATC {
         compute_force_matrix();
         // calculate force part of stress tensor
         local_potential_hardy_stress = atomicBondMatrix_*atomicForceMatrix_;
-        local_potential_hardy_stress *= 0.5; 
+        local_potential_hardy_stress *= 0.5;
       }
     }
     // global summation of potential part of stress tensor
@@ -1570,7 +1570,7 @@ namespace ATC {
       compute_kinetic_heatflux(flux);
     }
     else {
-      flux.zero(); // zero stress table for Lagrangian analysis 
+      flux.zero(); // zero stress table for Lagrangian analysis
     }
     // add potential part of heat flux vector
     int nrows = flux.nRows();
@@ -1628,11 +1628,11 @@ namespace ATC {
     // - e^0_I v_I + \sigma^T_I v_I
     for (int i = 0; i < nNodes_; i++) {
       double e_i = energy(i,0);
-      flux(i,0) += (e_i + stress(i,0))*velocity(i,0) 
+      flux(i,0) += (e_i + stress(i,0))*velocity(i,0)
                  + stress(i,3)*velocity(i,1)+ stress(i,4)*velocity(i,2);
-      flux(i,1) += (e_i + stress(i,1))*velocity(i,1) 
+      flux(i,1) += (e_i + stress(i,1))*velocity(i,1)
                  + stress(i,3)*velocity(i,0)+ stress(i,5)*velocity(i,2);
-      flux(i,2) += (e_i + stress(i,2))*velocity(i,2) 
+      flux(i,2) += (e_i + stress(i,2))*velocity(i,2)
                  + stress(i,4)*velocity(i,0)+ stress(i,5)*velocity(i,1);
     }
   }
@@ -1643,7 +1643,7 @@ namespace ATC {
     const DENS_MAT & rho = (restrictedCharge_->quantity());
     SPAR_MAT K;
     feEngine_->stiffness_matrix(K);
-    double permittivity = lammpsInterface_->dielectric(); 
+    double permittivity = lammpsInterface_->dielectric();
     permittivity *= LammpsInterface::instance()->epsilon0();
     K *= permittivity;
     BC_SET bcs;
@@ -1670,7 +1670,7 @@ namespace ATC {
 
     for (int i = 0; i < nLocal_; i++) {
       int atomIdx = internalToAtom_(i);
-      if (type[atomIdx] != 13) { 
+      if (type[atomIdx] != 13) {
         atomCnt(i,0) = myAtomicWeights(i,i);
         atomic_weight_sum += myAtomicWeights(i,i);
         number_atoms++;
@@ -1725,7 +1725,7 @@ namespace ATC {
 #ifndef H_BASED
         F(0,0) += 1.0; F(1,1) += 1.0; F(2,2) += 1.0;
 #endif
-        FT = F.transpose(); 
+        FT = F.transpose();
       }
       else if (atomToElementMapType_ == EULERIAN) {
         vector_to_symm_matrix(i,S,P);
@@ -1741,9 +1741,9 @@ namespace ATC {
         // Q stores (1-H)
         Q -= FT.transpose();
         DENS_MAT F(3,3);
-        F = inv(Q); 
+        F = inv(Q);
         FT = F.transpose();
-        ESH = FT*ESH; 
+        ESH = FT*ESH;
       }
       // copy to global
       matrix_to_vector(i,ESH,M);
@@ -1761,7 +1761,7 @@ namespace ATC {
     DENS_MAT_VEC &h = hField[DISPLACEMENT];
     h.assign(nsd_, DENS_MAT(nNodes_,nsd_));
     tField.assign(nsd_, DENS_MAT(nNodes_,nsd_));
-    // each row is the CB stress at a node stored in voigt form 
+    // each row is the CB stress at a node stored in voigt form
     T.reset(nNodes_,FieldSizes[CAUCHY_BORN_STRESS]);
     const double nktv2p = lammpsInterface_->nktv2p();
     const double fact = -lammpsInterface_->mvv2e()*nktv2p;
@@ -1779,7 +1779,7 @@ namespace ATC {
     DENS_MAT S(nNodes_,6);
     symm_dens_mat_vec_to_vector(tField,S);
     S *= fact;
-    
+
     // tField/S holds the 2nd P-K stress tensor. Transform to
     // Cauchy for EULERIAN analysis, transform to 1st P-K
     // for LAGRANGIAN analysis.
@@ -1799,7 +1799,7 @@ namespace ATC {
         FT  = transpose(F);
         double J = det(F);
         STRESS = F*PK2*FT;
-        STRESS *= 1/J; 
+        STRESS *= 1/J;
         symm_matrix_to_vector(i,STRESS,T);
       }
       else{
@@ -1810,7 +1810,7 @@ namespace ATC {
         STRESS = F*PK2;
         matrix_to_vector(i,STRESS,T);
       }
- 
+
     }
   }
   //---------------------------------------------------------------------------
@@ -1861,7 +1861,7 @@ namespace ATC {
   void ATC_Transfer::cauchy_born_entropic_energy(const DENS_MAT &H, DENS_MAT &E, const DENS_MAT &T)
   {
     FIELD_MATS uField;      // uField should contain temperature.
-    uField[TEMPERATURE] = T; 
+    uField[TEMPERATURE] = T;
     GRAD_FIELD_MATS hField;
     DENS_MAT_VEC &h = hField[DISPLACEMENT];
     h.assign(nsd_, DENS_MAT(nNodes_,nsd_));
@@ -1916,13 +1916,13 @@ namespace ATC {
           vector_to_matrix(i,H,F);
 
           F(0,0) += 1.0; F(1,1) += 1.0; F(2,2) += 1.0;
-          FT = F.transpose(); 
+          FT = F.transpose();
         }
         //
         double J = det(FT);
         FT *= 1/J;
         if (atomToElementMapType_ == EULERIAN) {
-          FT = inv(FT); 
+          FT = inv(FT);
         }
         S = P*FT;
         matrix_to_vector(i,S,stress);
@@ -1933,10 +1933,10 @@ namespace ATC {
     DENS_MAT & stretch, const DENS_MAT & H)
   {
     DENS_MAT F(3,3),R(3,3),U(3,3);
-    for (int i = 0; i < nNodes_; i++) { 
+    for (int i = 0; i < nNodes_; i++) {
       vector_to_matrix(i,H,F);
       F(0,0) += 1.0; F(1,1) += 1.0; F(2,2) += 1.0;
-      if (atomToElementMapType_ == EULERIAN) { 
+      if (atomToElementMapType_ == EULERIAN) {
         polar_decomposition(F,R,U,false); } // F = V R
       else  {
         polar_decomposition(F,R,U); } // F = R U
@@ -1953,12 +1953,12 @@ namespace ATC {
   //--------------------------------------------------------------------
   void ATC_Transfer::compute_elastic_deformation_gradient(DENS_MAT & Fe,
     const DENS_MAT & P, const DENS_MAT & H)
-  
+
   {
     // calculate Fe for every node
     const double nktv2p = lammpsInterface_->nktv2p();
     const double fact = 1.0/ ( lammpsInterface_->mvv2e()*nktv2p );
-    for (int i = 0; i < nNodes_; i++) { 
+    for (int i = 0; i < nNodes_; i++) {
       DENS_VEC Pv = global_vector_to_vector(i,P);
       Pv *= fact;
       CBElasticTangentOperator tangent(cauchyBornStress_, Pv);
@@ -1977,11 +1977,11 @@ namespace ATC {
     const double nktv2p = lammpsInterface_->nktv2p();
     const double fact = 1.0/ ( lammpsInterface_->mvv2e()*nktv2p );
     DENS_MAT F(3,3),R(3,3),U(3,3),PP(3,3),S(3,3);
-    for (int i = 0; i < nNodes_; i++) { 
+    for (int i = 0; i < nNodes_; i++) {
       // get F = RU
       vector_to_matrix(i,H,F);
       F(0,0) += 1.0; F(1,1) += 1.0; F(2,2) += 1.0;
-      if (atomToElementMapType_ == EULERIAN) { 
+      if (atomToElementMapType_ == EULERIAN) {
         polar_decomposition(F,R,U,false); } // F = V R
       else  {
         polar_decomposition(F,R,U); } // F = R U
@@ -1989,7 +1989,7 @@ namespace ATC {
       vector_to_matrix(i,P,PP);
       //S = PP*transpose(F);
       S = inv(F)*PP;
-      
+
       S += S.transpose(); S *= 0.5; // symmetrize
       DENS_VEC Sv = to_voigt(S);
       Sv *= fact;

lib/atc/ATC_Transfer.h

@@ -22,9 +22,9 @@ class TimeFilter;
 class ATC_Transfer : public ATC_Method {
 
  public:
-  
+
   // constructor
-  ATC_Transfer(std::string groupName, 
+  ATC_Transfer(std::string groupName,
                double **& perAtomArray,
                LAMMPS_NS::Fix * thisFix,
                std::string matParamFile = "none");
@@ -52,7 +52,7 @@ class ATC_Transfer : public ATC_Method {
   virtual void pre_neighbor() {ATC_Method::pre_neighbor(); neighborReset_ = true;};
 
   /** output function */
-  virtual void output(); 
+  virtual void output();
 
   /** external access to hardy data and other information*/
   const DENS_MAT * hardy_data(std::string field) { return &hardyData_[field].quantity(); }
@@ -63,25 +63,25 @@ class ATC_Transfer : public ATC_Method {
   double ** xPointer_;
 
   /** data */
-  TAG_FIELDS hardyData_; 
+  TAG_FIELDS hardyData_;
   SmallMoleculeSet * smallMoleculeSet_; // KKM add
   SmallMoleculeCentroid * moleculeCentroid_; // KKM add
   SmallMoleculeDipoleMoment * dipoleMoment_; // KKM add
   SmallMoleculeQuadrupoleMoment * quadrupoleMoment_; // KKM add
   /** container for dependency managed data */
   std::vector < DENS_MAN * > outputFields_;
-  
+
   std::map < std::string, DENS_MAN * > outputFieldsTagged_;
 
   DENS_MAN * restrictedCharge_; // WIP/TEMP
 
-  /** work space */ 
+  /** work space */
   DENS_MAT atomicScalar_;
   DENS_MAT atomicVector_;
   DENS_MAT atomicTensor_;
 
   /** calculation flags */
-  Array<bool> fieldFlags_; 
+  Array<bool> fieldFlags_;
   Array<bool> outputFlags_;
   Array<bool> gradFlags_;
   Array<bool> rateFlags_;
@@ -117,7 +117,7 @@ class ATC_Transfer : public ATC_Method {
   void compute_heatflux(DENS_MAT & flux);
   /** derived quantities: compute nodal to nodal quantities */
   void compute_eshelby_stress(DENS_MAT & eshebly_stress,
-    const DENS_MAT &  energy, const DENS_MAT & stress, 
+    const DENS_MAT &  energy, const DENS_MAT & stress,
     const DENS_MAT & displacement_gradient);
   void cauchy_born_stress(const DENS_MAT &dudx, DENS_MAT &T, const DENS_MAT *temp=0);
   void cauchy_born_energy(const DENS_MAT &dudx, DENS_MAT &T, const DENS_MAT *temp=0);
@@ -151,15 +151,15 @@ class ATC_Transfer : public ATC_Method {
   virtual void compute_dislocation_density(DENS_MAT & dislocation_density) = 0;
 
   /** compute smooth fields */
-  void compute_fields(void); 
-  void time_filter_pre (double dt); 
-  void time_filter_post(double dt); 
+  void compute_fields(void);
+  void time_filter_pre (double dt);
+  void time_filter_post(double dt);
 
   /** mapping of atomic pairs to pair index value */
-  class PairMap * pairMap_; 
-  class BondMatrix * bondMatrix_; 
-  class PairVirial * pairVirial_; 
-  class PairPotentialHeatFlux * pairHeatFlux_; 
+  class PairMap * pairMap_;
+  class BondMatrix * bondMatrix_;
+  class PairVirial * pairVirial_;
+  class PairPotentialHeatFlux * pairHeatFlux_;
 
   /** routine to calculate matrix of force & position dyads */
   void compute_force_matrix();
@@ -176,7 +176,7 @@ class ATC_Transfer : public ATC_Method {
                                   DENS_MAT & nodeData) = 0;
 
   /** routine to calculate matrix of bond functions */
-  virtual void compute_bond_matrix(); 
+  virtual void compute_bond_matrix();
 
   /** routine to set xPointer to xref or xatom */
   void set_xPointer();
@@ -200,21 +200,21 @@ class ATC_Transfer : public ATC_Method {
   void project_count_normalized(const DENS_MAT & atomData,
                                       DENS_MAT & nodeData);
 
-  /** hardy_project (volume density): given w_\alpha, 
-               w_I = 1/\Omega_I \sum_\alpha N_{I\alpha} w_\alpha 
+  /** hardy_project (volume density): given w_\alpha,
+               w_I = 1/\Omega_I \sum_\alpha N_{I\alpha} w_\alpha
                where   \Omega_I = \int_{support region of node I} N_{I} dV  */
   // REFACTOR AtfNodeWeightedShapeFunctionRestriction
   void project_volume_normalized(const DENS_MAT & atomData,
                                        DENS_MAT & nodeData);
   void project_volume_normalized_molecule(const DENS_MAT & molData,
-                                                DENS_MAT & nodeData); // KKM add    
+                                                DENS_MAT & nodeData); // KKM add
   void project_volume_normalized_molecule_gradient(const DENS_MAT & molData,
-                                                DENS_MAT & nodeData); // KKM add    
-  
-  
-  /** gradient_compute: given w_I, 
-               w_J = \sum_I N_I'{xJ} \dyad w_I 
-               where N_I'{xJ} is the gradient of the normalized 
+                                                DENS_MAT & nodeData); // KKM add
+
+
+  /** gradient_compute: given w_I,
+               w_J = \sum_I N_I'{xJ} \dyad w_I
+               where N_I'{xJ} is the gradient of the normalized
                               shape function of node I evaluated at node J */
   // REFACTOR MatToGradBySparse
   void gradient_compute(const DENS_MAT & inNodeData,
@@ -226,7 +226,7 @@ class ATC_Transfer : public ATC_Method {
   /** workset data */
   VectorDependencyManager<SPAR_MAT * > * gradientMatrix_;
 
-  
+
   SPAR_MAT atomicBondMatrix_;
   DENS_MAT atomicForceMatrix_;
   DENS_MAT atomicHeatMatrix_;
@@ -247,7 +247,7 @@ class ATC_Transfer : public ATC_Method {
   Array<TimeFilter *> timeFilters_;
 
   /** check consistency of fieldFlags_ */
-  void check_field_dependencies(); 
+  void check_field_dependencies();
 
 };
 

lib/atc/ATC_TransferKernel.cpp

@@ -1,4 +1,4 @@
-// ATC headers 
+// ATC headers
 #include "ATC_TransferKernel.h"
 #include "ATC_Error.h"
 #include "FE_Engine.h"
@@ -44,10 +44,10 @@ using ATC_Utility::to_string;
   {
     bool match = false;
 
-    /*! \page man_hardy_kernel fix_modify AtC kernel 
+    /*! \page man_hardy_kernel fix_modify AtC kernel
       \section syntax
       fix_modify AtC kernel <type> <parameters>
-      - type (keyword) = step, cell, cubic_bar, cubic_cylinder, cubic_sphere, 
+      - type (keyword) = step, cell, cubic_bar, cubic_cylinder, cubic_sphere,
                          quartic_bar, quartic_cylinder, quartic_sphere \n
       - parameters :\n
       step = radius (double) \n
@@ -62,7 +62,7 @@ using ATC_Utility::to_string;
       <TT> fix_modify AtC kernel cell 1.0 1.0 1.0 </TT> \n
       <TT> fix_modify AtC kernel quartic_sphere 10.0 </TT>
       \section description
-      
+
       \section restrictions
       Must be used with the hardy AtC fix \n
       For bar kernel types, half-width oriented along x-direction \n
@@ -92,28 +92,28 @@ using ATC_Utility::to_string;
       SPAR_MAT & dN(kernelAccumulantMolGrad_.set_quantity());
       dN.reset(nLocalMol,nNodes_);
       DENS_VEC derivKer(nsd_);
-      DENS_VEC xI(nsd_),xm(nsd_),xmI(nsd_); 
+      DENS_VEC xI(nsd_),xm(nsd_),xmI(nsd_);
       const DENS_MAT & centroidMolMatrix(moleculeCentroid_->quantity());
-      ATC::LammpsInterface::instance()->stream_msg_once("computing kernel matrix molecule ",true,false); 
-      int heartbeatFreq = (nNodes_ <= 10 ? 1 : (int) nNodes_ / 10); 
+      ATC::LammpsInterface::instance()->stream_msg_once("computing kernel matrix molecule ",true,false);
+      int heartbeatFreq = (nNodes_ <= 10 ? 1 : (int) nNodes_ / 10);
       for (int i = 0; i < nNodes_; i++) {
-        if (i % heartbeatFreq == 0 ) { 
+        if (i % heartbeatFreq == 0 ) {
           ATC::LammpsInterface::instance()->stream_msg_once(".",false,false);
-        }  
+        }
         xI = (feEngine_->fe_mesh())->nodal_coordinates(i);
         for (int j = 0; j < nLocalMol; j++) {
           for (int k = 0; k < nsd_; k++) {
             xm(k) = centroidMolMatrix(j,k);
           }
-          xmI = xm - xI; 
+          xmI = xm - xI;
           lammpsInterface_->periodicity_correction(xmI.ptr());
           double val = kernelFunction_->value(xmI);
           if (val > 0) N.add(j,i,val);
           kernelFunction_->derivative(xmI,derivKer);
           double val_grad = derivKer(2);
           if (val_grad!= 0) dN.add(j,i,val_grad);
-        }   
-      } 
+        }
+      }
       // reset kernelShpFunctions with the weights of molecules on processors
       DENS_VEC fractions(N.nRows());
       DENS_VEC fractions_deriv(dN.nRows());
@@ -126,10 +126,10 @@ using ATC_Utility::to_string;
         dN.compress();
       if (lammpsInterface_->rank_zero()) {
         ATC::LammpsInterface::instance()->stream_msg_once("done",false,true);
-      }  
-    }   
+      }
+    }
   }
-  
+
   //-------------------------------------------------------------------
   void ATC_TransferKernel::compute_projection(const DENS_MAT & atomData,
                                               DENS_MAT & nodeData)
@@ -183,7 +183,7 @@ using ATC_Utility::to_string;
     int **firstneigh = lammpsInterface_->neighbor_list_firstneigh();
     double ** xatom    = lammpsInterface_->xatom();
     double lam1,lam2;
-    double bond_value; 
+    double bond_value;
     // process differently for mesh vs translation-invariant kernels
     ATC::LammpsInterface::instance()->stream_msg_once("computing potential stress: ",true,false);
     int heartbeatFreq = (nNodes_ <= 10 ? 1 : (int) nNodes_ / 10);
@@ -202,7 +202,7 @@ using ATC_Utility::to_string;
       int inode = i;
       for (int j = 0; j < nLocal_; j++) {
         // second (neighbor) atom location
-        int lammps_j = internalToAtom_(j); 
+        int lammps_j = internalToAtom_(j);
         xa.copy(xPointer_[lammps_j],3);
         // difference vector
         xaI = xa - xI;
@@ -217,8 +217,8 @@ using ATC_Utility::to_string;
           kernelFunction_->bond_intercepts(xaI,xbI,lam1,lam2);
           // compute virial
           if (lam1 < lam2) {
-            bond_value 
-              = kernel_inv_vol*(kernelFunction_->bond(xaI,xbI,lam1,lam2)); 
+            bond_value
+              = kernel_inv_vol*(kernelFunction_->bond(xaI,xbI,lam1,lam2));
             double delx = xatom[lammps_j][0] - xatom[lammps_k][0];
             double dely = xatom[lammps_j][1] - xatom[lammps_k][1];
             double delz = xatom[lammps_j][2] - xatom[lammps_k][2];
@@ -227,9 +227,9 @@ using ATC_Utility::to_string;
             lammpsInterface_->pair_force(lammps_j,lammps_k,rsq,fforce);
             fforce *= 0.5; // dbl count
             if (atomToElementMapType_ == LAGRANGIAN) {
-              double delX = xref_[lammps_j][0] - xref_[lammps_k][0]; 
-              double delY = xref_[lammps_j][1] - xref_[lammps_k][1]; 
-              double delZ = xref_[lammps_j][2] - xref_[lammps_k][2]; 
+              double delX = xref_[lammps_j][0] - xref_[lammps_k][0];
+              double delY = xref_[lammps_j][1] - xref_[lammps_k][1];
+              double delZ = xref_[lammps_j][2] - xref_[lammps_k][2];
               stress(inode,0) +=-delx*fforce*delX*bond_value;
               stress(inode,1) +=-delx*fforce*delY*bond_value;
               stress(inode,2) +=-delx*fforce*delZ*bond_value;
@@ -266,9 +266,9 @@ using ATC_Utility::to_string;
     int **firstneigh = lammpsInterface_->neighbor_list_firstneigh();
     double ** xatom    = lammpsInterface_->xatom();
     double ** vatom    = lammpsInterface_->vatom();
-    
+
     double lam1,lam2;
-    double bond_value; 
+    double bond_value;
     // process differently for mesh vs translation-invariant kernels
     // "normal" kernel functions
     DENS_VEC xa(nsd_),xI(nsd_),xaI(nsd_),xb(nsd_),xbI(nsd_),xba(nsd_);
@@ -281,7 +281,7 @@ using ATC_Utility::to_string;
         continue;
       }
       for (int j = 0; j < nLocal_; j++) {
-        int lammps_j = internalToAtom_(j); 
+        int lammps_j = internalToAtom_(j);
         xa.copy(xPointer_[lammps_j],3);
         // difference vector
         xaI = xa - xI;
@@ -296,8 +296,8 @@ using ATC_Utility::to_string;
           kernelFunction_->bond_intercepts(xaI,xbI,lam1,lam2);
           // compute virial
           if (lam1 < lam2) {
-            bond_value 
-              = kernel_inv_vol*(kernelFunction_->bond(xaI,xbI,lam1,lam2)); 
+            bond_value
+              = kernel_inv_vol*(kernelFunction_->bond(xaI,xbI,lam1,lam2));
             double delx = xatom[lammps_j][0] - xatom[lammps_k][0];
             double dely = xatom[lammps_j][1] - xatom[lammps_k][1];
             double delz = xatom[lammps_j][2] - xatom[lammps_k][2];
@@ -308,9 +308,9 @@ using ATC_Utility::to_string;
             double * v = vatom[lammps_j];
             fforce *= (delx*v[0] + dely*v[1] + delz*v[2]);
             if (atomToElementMapType_ == LAGRANGIAN) {
-              double delX = xref_[lammps_j][0] - xref_[lammps_k][0]; 
-              double delY = xref_[lammps_j][1] - xref_[lammps_k][1]; 
-              double delZ = xref_[lammps_j][2] - xref_[lammps_k][2]; 
+              double delX = xref_[lammps_j][0] - xref_[lammps_k][0];
+              double delY = xref_[lammps_j][1] - xref_[lammps_k][1];
+              double delZ = xref_[lammps_j][2] - xref_[lammps_k][2];
               flux(inode,0) +=fforce*delX*bond_value;
               flux(inode,1) +=fforce*delY*bond_value;
               flux(inode,2) +=fforce*delZ*bond_value;
@@ -327,7 +327,7 @@ using ATC_Utility::to_string;
   }
 
   //-------------------------------------------------------------------
-  // calculation of the dislocation density tensor 
+  // calculation of the dislocation density tensor
   void ATC_TransferKernel::compute_dislocation_density(DENS_MAT & A)
   {
    A.reset(nNodes_,9);
@@ -348,7 +348,7 @@ using ATC_Utility::to_string;
    lammpsInterface_->int_allsum(&localNumberLines,&totalNumberLines,1);
    if (totalNumberLines == 0) {
      ATC::LammpsInterface::instance()->print_msg_once("no dislocation lines found");
-     return; 
+     return;
    }
 
    // for output
@@ -366,7 +366,7 @@ using ATC_Utility::to_string;
 
 
    DENS_MAT local_A(nNodes_,9);
-   local_A.zero(); 
+   local_A.zero();
    DENS_VEC xa(nsd_),xI(nsd_),xaI(nsd_),xb(nsd_),xbI(nsd_),xba(nsd_);
    double kernel_inv_vol = kernelFunction_->inv_vol();
    int iPt = 0, iSeg= 0;
@@ -393,7 +393,7 @@ using ATC_Utility::to_string;
          xa(k) = x1[k];
          xb(k) = x2[k];
          xba(k) = delta[k];
-       } 
+       }
        for (int I = 0; I < nNodes_; I++) {
          xI = (feEngine_->fe_mesh())->nodal_coordinates(I);
          if (!kernelFunction_->node_contributes(xI)) {
@@ -405,7 +405,7 @@ using ATC_Utility::to_string;
          double lam1=0,lam2=0;
          kernelFunction_->bond_intercepts(xaI,xbI,lam1,lam2);
          if (lam1 < lam2) {
-           double bond_value 
+           double bond_value
              = kernel_inv_vol*(kernelFunction_->bond(xaI,xbI,lam1,lam2));
            local_A(I,0) += xba(0)*burgers[0]*bond_value;
            local_A(I,1) += xba(0)*burgers[1]*bond_value;
@@ -449,7 +449,7 @@ using ATC_Utility::to_string;
    lammpsInterface_->int_allsum(&nSeg,&totalNumberSegments,1);
 
    // output
-   double volume = lammpsInterface_->domain_volume(); 
+   double volume = lammpsInterface_->domain_volume();
    stringstream ss;
    ss << "total dislocation line length = " << totalDislocationDensity;
    ss << " lines = " << totalNumberLines << " segments = " << totalNumberSegments;
@@ -474,10 +474,10 @@ using ATC_Utility::to_string;
      segOutput.write_geometry(&segCoor,&segConn);
      OUTPUT_LIST segOut;
      segOut["burgers_vector"] = &segBurg;
-     segOutput.write_data(0,&segOut); 
+     segOutput.write_data(0,&segOut);
    }
 #else
-  throw ATC_Error("unimplemented function compute_dislocation_density (DXA support not included"); 
+  throw ATC_Error("unimplemented function compute_dislocation_density (DXA support not included");
 #endif
   }