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

.github/workflows/codeql-analysis.yml

@@ -35,7 +35,7 @@ jobs:
         python-version: '3.x'
 
     - name: Initialize CodeQL
-      uses: github/codeql-action/init@v2
+      uses: github/codeql-action/init@v3
       with:
         languages: ${{ matrix.language }}
         config-file: ./.github/codeql/${{ matrix.language }}.yml
@@ -55,4 +55,4 @@ jobs:
         cmake --build . --parallel 2
 
     - name: Perform CodeQL Analysis
-      uses: github/codeql-action/analyze@v2
+      uses: github/codeql-action/analyze@v3

.github/workflows/coverity.yml

@@ -25,7 +25,7 @@ jobs:
 
     - name: Cache Coverity
       id: cache-coverity
-      uses: actions/cache@v3
+      uses: actions/cache@v4
       with:
         path: ./download/
         key: ${{ runner.os }}-download-${{ hashFiles('**/coverity_tool.*') }}

.github/workflows/unittest-macos.yml

@@ -32,7 +32,7 @@ jobs:
       run: mkdir build
 
     - name: Set up ccache
-      uses: actions/cache@v3
+      uses: actions/cache@v4
       with:
         path: ${{ env.CCACHE_DIR }}
         key: macos-ccache-${{ github.sha }}

cmake/Modules/Packages/GPU.cmake

@@ -1,3 +1,10 @@
+
+# Silence CMake warnings about FindCUDA being obsolete.
+# We may need to eventually rewrite this section to use enable_language(CUDA)
+if(POLICY CMP0146)
+  cmake_policy(SET CMP0146 OLD)
+endif()
+
 set(GPU_SOURCES_DIR ${LAMMPS_SOURCE_DIR}/GPU)
 set(GPU_SOURCES ${GPU_SOURCES_DIR}/gpu_extra.h
                 ${GPU_SOURCES_DIR}/fix_gpu.h

cmake/Modules/Packages/KOKKOS.cmake

@@ -16,11 +16,6 @@ endif()
 if(Kokkos_ENABLE_OPENMP)
   if(NOT BUILD_OMP)
     message(FATAL_ERROR "Must enable BUILD_OMP with Kokkos_ENABLE_OPENMP")
-  else()
-    # NVHPC/(AMD)Clang does not seem to provide a detectable OpenMP version, but is far beyond version 3.1
-    if((OpenMP_CXX_VERSION VERSION_LESS 3.1) AND NOT ((CMAKE_CXX_COMPILER_ID STREQUAL "NVHPC") OR (CMAKE_CXX_COMPILER_ID STREQUAL "Clang")))
-      message(FATAL_ERROR "Compiler must support OpenMP 3.1 or later with Kokkos_ENABLE_OPENMP")
-    endif()
   endif()
 endif()
 ########################################################################

cmake/Modules/Packages/MDI.cmake

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

cmake/presets/oneapi.cmake

@@ -18,11 +18,11 @@ set(MPI_CXX_COMPILER "mpicxx" CACHE STRING "" FORCE)
 
 unset(HAVE_OMP_H_INCLUDE CACHE)
 set(OpenMP_C "icx" CACHE STRING "" FORCE)
-set(OpenMP_C_FLAGS "-qopenmp -qopenmp-simd" CACHE STRING "" FORCE)
+set(OpenMP_C_FLAGS "-qopenmp;-qopenmp-simd" CACHE STRING "" FORCE)
 set(OpenMP_C_LIB_NAMES "omp" CACHE STRING "" FORCE)
 set(OpenMP_CXX "icpx" CACHE STRING "" FORCE)
-set(OpenMP_CXX_FLAGS "-qopenmp -qopenmp-simd" CACHE STRING "" FORCE)
+set(OpenMP_CXX_FLAGS "-qopenmp;-qopenmp-simd" CACHE STRING "" FORCE)
 set(OpenMP_CXX_LIB_NAMES "omp" CACHE STRING "" FORCE)
-set(OpenMP_Fortran_FLAGS "-qopenmp -qopenmp-simd" CACHE STRING "" FORCE)
+set(OpenMP_Fortran_FLAGS "-qopenmp;-qopenmp-simd" CACHE STRING "" FORCE)
 set(OpenMP_omp_LIBRARY "libiomp5.so" CACHE PATH "" FORCE)
 

doc/Makefile

@@ -100,6 +100,7 @@ html: xmlgen $(VENV) $(SPHINXCONFIG)/conf.py $(ANCHORCHECK) $(MATHJAX)
 		env LC_ALL=C grep -n ' :[a-z]\+`' $(RSTDIR)/*.rst ;\
 		env LC_ALL=C grep -n ' `[^`]\+<[a-z][^`]\+`[^_]' $(RSTDIR)/*.rst ;\
 		env LC_ALL=C grep -n ':\(ref\|doc\):[^`]' $(RSTDIR)/*.rst ;\
+		env LC_ALL=C grep -n '\(ref\|doc\)`[^`]' $(RSTDIR)/*.rst ;\
 		$(PYTHON) $(BUILDDIR)/utils/check-styles.py -s ../src -d src ;\
 		echo "############################################" ;\
 		deactivate ;\
@@ -182,6 +183,7 @@ pdf: xmlgen $(VENV) $(SPHINXCONFIG)/conf.py $(ANCHORCHECK)
 		env LC_ALL=C grep -n ' :[a-z]\+`' $(RSTDIR)/*.rst ;\
 		env LC_ALL=C grep -n ' `[^`]\+<[a-z][^`]\+`[^_]' $(RSTDIR)/*.rst ;\
 		env LC_ALL=C grep -n ':\(ref\|doc\):[^`]' $(RSTDIR)/*.rst ;\
+		env LC_ALL=C grep -n '\(ref\|doc\)`[^`]' $(RSTDIR)/*.rst ;\
 		$(PYTHON) utils/check-styles.py -s ../src -d src ;\
 		echo "############################################" ;\
 		deactivate ;\
@@ -231,6 +233,7 @@ role_check :
 	@( env LC_ALL=C grep -n ' :[a-z]\+`' $(RSTDIR)/*.rst && exit 1 || : )
 	@( env LC_ALL=C grep -n ' `[^`]\+<[a-z][^`]\+`[^_]' $(RSTDIR)/*.rst && exit 1 || : )
 	@( env LC_ALL=C grep -n ':\(ref\|doc\):[^`]' $(RSTDIR)/*.rst && exit 1 || : )
+	@( env LC_ALL=C grep -n '\(ref\|doc\)`[^`]' $(RSTDIR)/*.rst && exit 1 || : )
 
 link_check : $(VENV) html
 	@(\

doc/github-development-workflow.md

@@ -36,10 +36,10 @@ requests.
 MUST be submitted as a pull request to GitHub.  All changes to the
 "develop" branch must be made exclusively through merging pull requests.
 The "release" and "stable" branches, respectively, are only to be
-updated upon feature or stable releases based on the associated
+updated upon "feature releases" or "stable releases" based on the
-tags.  Updates to the stable release in between stable releases
+associated tags.  Updates to the stable release in between stable releases
 (for example, back-ported bug fixes) are first merged into the "maintenance"
-branch and then into the "stable" branch as update releases.
+branch and then into the "stable" branch as "stable update releases".
 
 Pull requests may also be submitted to (long-running) feature branches
 created by LAMMPS developers inside the LAMMPS project, if needed. Those
@@ -131,7 +131,7 @@ testing -- that the code in the branch "develop" does not get easily
 broken.  These tests are run after every update to a pull request.  More
 extensive and time-consuming tests (including regression testing) are
 performed after code is merged to the "develop" branch.  There are feature
-releases of LAMMPS made about every 4-6 weeks at a point, when the LAMMPS
+releases of LAMMPS made about every 4-8 weeks at a point, when the LAMMPS
 developers feel, that a sufficient number of changes have been included
 and all post-merge testing has been successful.  These feature releases are
 marked with a `patch_<version date>` tag and the "release" branch

doc/graphviz/Makefile

@@ -16,8 +16,11 @@ clean:
 	rm -f $(IMGSVG) $(IMGPDF) $(IMGPNG) *~
 
 ifeq ($(HAS_DOT),YES)
-$(IMGDIR)/%.png: %.dot
+$(IMGDIR)/lammps-classes.png : lammps-classes.dot
 	dot -Tpng -Kneato -o $@ $<
+
+$(IMGDIR)/%.png: %.dot
+	dot -Tpng -Kdot -o $@ $<
 endif
 
 ifeq ($(HAS_DOT),NO)

doc/graphviz/lammps-releases.dot

@@ -0,0 +1,34 @@
+// LAMMPS branches and releases
+digraph releases {
+    rankdir="LR";
+    github [shape="box" label="Pull Requests\non GitHub" height=0.75];
+    github -> develop [label="Merge commits"];
+    {
+        rank = "same";
+        work [shape="none" label="Development branches:"]
+        develop [label="'develop' branch" height=0.75];
+        maintenance [label="'maintenance' branch" height=0.75];
+    };
+    {
+        rank = "same";
+        upload [shape="none" label="Release branches:"]
+        release [label="'release' branch" height=0.75];
+        stable [label="'stable' branch" height=0.75];
+    };
+    develop -> release [label="Feature release\n(every 4-8 weeks)"];
+    release -> stable [label="Stable release\n(once per year)"];
+    stable -> maintenance [label="Reset on stable release" style="setlinewidth(2)"];
+    develop -> maintenance [label="Backports of bugfixes" style="dashed"];
+    maintenance -> stable [label="Updates to stable release"];
+    {
+        rank = "same";
+        tag [shape="none" label="Applied tags:"];
+        patchtag [shape="box" label="patch_<date>"];
+        stabletag [shape="box" label="stable_<date>"];
+        updatetag [shape="box" label="stable_<date>_update<num>"];
+    };
+    release -> patchtag [label="feature release" style="dotted"];
+    stable -> stabletag [label="stable release" style="dotted"];
+    stable -> updatetag [label="update release" style="dotted"];
+}
+

doc/src/Build_development.rst

@@ -255,16 +255,18 @@ A test run is then a a collection multiple individual test runs each
 with many comparisons to reference results based on template input
 files, individual command settings, relative error margins, and
 reference data stored in a YAML format file with ``.yaml``
-suffix. Currently the programs ``test_pair_style``, ``test_bond_style``, and
+suffix. Currently the programs ``test_pair_style``, ``test_bond_style``,
-``test_angle_style`` are implemented.  They will compare forces, energies and
+``test_angle_style``, ``test_dihedral_style``, and
-(global) stress for all atoms after a ``run 0`` calculation and after a
+``test_improper_style`` are implemented.  They will compare forces,
-few steps of MD with :doc:`fix nve <fix_nve>`, each in multiple variants
+energies and (global) stress for all atoms after a ``run 0`` calculation
-with different settings and also for multiple accelerated styles. If a
+and after a few steps of MD with :doc:`fix nve <fix_nve>`, each in
-prerequisite style or package is missing, the individual tests are
+multiple variants with different settings and also for multiple
-skipped.  All tests will be executed on a single MPI process, so using
+accelerated styles. If a prerequisite style or package is missing, the
-the CMake option ``-D BUILD_MPI=off`` can significantly speed up testing,
+individual tests are skipped.  All force style tests will be executed on
-since this will skip the MPI initialization for each test run.
+a single MPI process, so using the CMake option ``-D BUILD_MPI=off`` can
-Below is an example command and output:
+significantly speed up testing, since this will skip the MPI
+initialization for each test run.  Below is an example command and
+output:
 
 .. code-block:: console
 
@@ -416,15 +418,16 @@ When compiling LAMMPS with enabled tests, most test executables will
 need to be linked against the LAMMPS library.  Since this can be a very
 large library with many C++ objects when many packages are enabled, link
 times can become very long on machines that use the GNU BFD linker (e.g.
-Linux systems).  Alternatives like the ``lld`` linker of the LLVM project
+Linux systems).  Alternatives like the ``mold`` linker, the ``lld``
-or the ``gold`` linker available with GNU binutils can speed up this step
+linker of the LLVM project, or the ``gold`` linker available with GNU
-substantially. CMake will by default test if any of the two can be
+binutils can speed up this step substantially (in this order).  CMake
-enabled and use it when ``ENABLE_TESTING`` is active.  It can also be
+will by default test if any of the three can be enabled and use it when
-selected manually through the ``CMAKE_CUSTOM_LINKER`` CMake variable.
+``ENABLE_TESTING`` is active.  It can also be selected manually through
-Allowed values are ``lld``, ``gold``, ``bfd``, or ``default``.  The
+the ``CMAKE_CUSTOM_LINKER`` CMake variable.  Allowed values are
-``default`` option will use the system default linker otherwise, the
+``mold``, ``lld``, ``gold``, ``bfd``, or ``default``.  The ``default``
-linker is chosen explicitly.  This option is only available for the
+option will use the system default linker otherwise, the linker is
-GNU or Clang C++ compiler.
+chosen explicitly.  This option is only available for the GNU or Clang
+C++ compilers.
 
 Tests for other components and utility functions
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

doc/src/Commands_bond.rst

@@ -124,7 +124,7 @@ OPT.
    *
    *
    * :doc:`charmm (iko) <dihedral_charmm>`
-   * :doc:`charmmfsw <dihedral_charmm>`
+   * :doc:`charmmfsw (k) <dihedral_charmm>`
    * :doc:`class2 (ko) <dihedral_class2>`
    * :doc:`cosine/shift/exp (o) <dihedral_cosine_shift_exp>`
    * :doc:`fourier (io) <dihedral_fourier>`

doc/src/Commands_pair.rst

@@ -146,7 +146,7 @@ OPT.
    * :doc:`lj/charmm/coul/long/soft (o) <pair_fep_soft>`
    * :doc:`lj/charmm/coul/msm (o) <pair_charmm>`
    * :doc:`lj/charmmfsw/coul/charmmfsh <pair_charmm>`
-   * :doc:`lj/charmmfsw/coul/long <pair_charmm>`
+   * :doc:`lj/charmmfsw/coul/long (k) <pair_charmm>`
    * :doc:`lj/class2 (gko) <pair_class2>`
    * :doc:`lj/class2/coul/cut (ko) <pair_class2>`
    * :doc:`lj/class2/coul/cut/soft <pair_fep_soft>`

doc/src/Developer_platform.rst

@@ -70,6 +70,9 @@ File and path functions and global constants
 .. doxygenfunction:: is_console
    :project: progguide
 
+.. doxygenfunction:: disk_free
+   :project: progguide
+
 .. doxygenfunction:: path_is_directory
    :project: progguide
 

doc/src/Developer_unittest.rst

@@ -121,7 +121,7 @@ will be suppressed and only a summary printed, but adding
 the '-V' option will then produce output from the tests
 above like the following:
 
-.. code-block::
+.. code-block:: console
 
    [...]
    1: [ RUN      ] Tokenizer.empty_string
@@ -526,3 +526,102 @@ The ``unittest/tools`` folder contains tests for programs in the
 shell, which are implemented as a python scripts using the ``unittest``
 Python module and launching the tool commands through the ``subprocess``
 Python module.
+
+
+Troubleshooting failed unit tests
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The are by default no unit tests for newly added features (e.g. pair, fix,
+or compute styles) unless your pull request also includes tests for the
+added features.  If you are modifying some features, you may see failures
+for existing tests, if your modifications have some unexpected side effects
+or your changes render the existing text invalid.  If you are adding an
+accelerated version of an existing style, then only tests for INTEL,
+KOKKOS (with OpenMP only), OPENMP, and OPT will be run automatically.
+Tests for the GPU package are time consuming and thus are only run
+*after* a merge, or when a special label, ``gpu_unit_tests`` is added
+to the pull request.  After the test has started, it is often best to
+remove the label since every PR activity will re-trigger the test (that
+is a limitation of triggering a test with a label).  Support for unit
+tests with using KOKKOS with GPU acceleration is currently not supported.
+
+When you see a failed build on GitHub, click on ``Details`` to be taken
+to the corresponding LAMMPS Jenkins CI web page.  Click on the "Exit"
+symbol near the ``Logout`` button on the top right of that page to go to
+the "classic view".  In the classic view, there is a list of the
+individual runs that make up this test run (they are shown but cannot be
+inspected in the default view).  You can click on any of those.
+Clicking on ``Test Result`` will display the list of failed tests. Click
+on the "Status" column to sort the tests based on their Failed or Passed
+status.  Then click on the failed test to expand its output.
+
+For example, the following output snippet shows the failed unit test
+
+.. code-block:: console
+
+   [ RUN      ] PairStyle.gpu
+   /home/builder/workspace/dev/pull_requests/ubuntu_gpu/unit_tests/cmake_gpu_opencl_mixed_smallbig_clang_static/unittest/force-styles/test_main.cpp:63: Failure
+   Expected: (err) <= (epsilon)
+   Actual: 0.00018957912910606503 vs 0.0001
+   Google Test trace:
+   /home/builder/workspace/dev/pull_requests/ubuntu_gpu/unit_tests/cmake_gpu_opencl_mixed_smallbig_clang_static/unittest/force-styles/test_main.cpp:56: EXPECT_FORCES: init_forces (newton off)
+   /home/builder/workspace/dev/pull_requests/ubuntu_gpu/unit_tests/cmake_gpu_opencl_mixed_smallbig_clang_static/unittest/force-styles/test_main.cpp:64: Failure
+   Expected: (err) <= (epsilon)
+   Actual: 0.00022892713393549854 vs 0.0001
+
+The failed assertions provide line numbers in the test source
+(e.g. ``test_main.cpp:56``), from which one can understand what
+specific assertion failed.
+
+Note that the force style engine runs one of a small number of systems
+in a rather off-equilibrium configuration with a few atoms for a few
+steps, writes data and restart files, uses :doc:`the clear command
+<clear>` to reset LAMMPS, and then runs from those files with different
+settings (e.g. newton on/off) and integrators (e.g. verlet vs. respa).
+Beyond potential issues/bugs in the source code, the mismatch between
+the expected and actual values could be that force arrays are not
+properly cleared between multiple run commands or that class members are
+not correctly initialized or written to or read from a data or restart
+file.
+
+While the epsilon (relative precision) for a single, `IEEE 754 compliant
+<https://en.wikipedia.org/wiki/IEEE_754>`_, double precision floating
+point operation is at about 2.2e-16, the achievable precision for the
+tests is lower due to most numbers being sums over intermediate results
+and the non-associativity of floating point math leading to larger
+errors.  In some cases specific properties of the tested style.  As a
+rule of thumb, the test epsilon can often be in the range 5.0e-14 to
+1.0e-13.  But for "noisy" force kernels, e.g. those a larger amount of
+arithmetic operations involving `exp()`, `log()` or `sin()` functions,
+and also due to the effect of compiler optimization or differences
+between compilers or platforms, epsilon may need to be further relaxed,
+sometimes epsilon can be relaxed to 1.0e-12. If interpolation or lookup
+tables are used, epsilon may need to be set to 1.0e-10 or even higher.
+For tests of accelerated styles, the per-test epsilon is multiplied
+by empirical factors that take into account the differences in the order
+of floating point operations or that some or most intermediate operations
+may be done using approximations or with single precision floating point
+math.
+
+To rerun the failed unit test individually, change to the ``build`` directory
+and run the test with verbose output. For example,
+
+.. code-block:: bash
+
+    env TEST_ARGS=-v ctest -R ^MolPairStyle:lj_cut_coul_long -V
+
+``ctest`` with the ``-V`` flag also shows the exact command line
+of the test. One can then use ``gdb --args`` to further debug and
+catch exceptions with the test command, for example,
+
+.. code-block:: bash
+
+    gdb --args /path/to/lammps/build/test_pair_style /path/to/lammps/unittest/force-styles/tests/mol-pair-lj_cut_coul_long.yaml
+
+
+It is recommended to configure the build with ``-D
+BUILD_SHARED_LIBS=on`` and use a custom linker to shorten the build time
+during recompilation.  Installing `ccache` in your development
+environment helps speed up recompilation by caching previous
+compilations and detecting when the same compilation is being done
+again.  Please see :doc:`Build_development` for further details.

doc/src/Developer_updating.rst

@@ -20,6 +20,7 @@ Available topics in mostly chronological order are:
 - `Use ev_init() to initialize variables derived from eflag and vflag`_
 - `Use utils::numeric() functions instead of force->numeric()`_
 - `Use utils::open_potential() function to open potential files`_
+- `Use symbolic Atom and AtomVec constants instead of numerical values`_
 - `Simplify customized error messages`_
 - `Use of "override" instead of "virtual"`_
 - `Simplified and more compact neighbor list requests`_
@@ -196,6 +197,71 @@ New:
 
    fp = utils::open_potential(filename, lmp);
 
+Use symbolic Atom and AtomVec constants instead of numerical values
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. versionchanged:: 18Sep2020
+
+Properties in LAMMPS that were represented by integer values (0, 1,
+2, 3) to indicate settings in the ``Atom`` and ``AtomVec`` classes (or
+classes derived from it) (and its derived classes) have been converted
+to use scoped enumerators instead.
+
+.. list-table::
+   :header-rows: 1
+   :widths: auto
+
+   * - Symbolic Constant
+     - Value
+     - Symbolic Constant
+     - Value
+   * - Atom::GROW
+     - 0
+     - Atom::MAP_NONE
+     - 0
+   * - Atom::RESTART
+     - 1
+     - Atom::MAP_ARRAY
+     - 1
+   * - Atom::BORDER
+     - 2
+     - Atom::MAP_HASH
+     - 2
+   * - Atom::ATOMIC
+     - 0
+     - Atom::MAP_YES
+     - 3
+   * - Atom::MOLECULAR
+     - 1
+     - AtomVec::PER_ATOM
+     - 0
+   * - Atom::TEMPLATE
+     - 2
+     - AtomVec::PER_TYPE
+     - 1
+
+Old:
+
+.. code-block:: c++
+
+   molecular = 0;
+   mass_type = 1;
+   if (atom->molecular == 2)
+   if (atom->map_style == 2)
+   atom->add_callback(0);
+   atom->delete_callback(id,1);
+
+New:
+
+.. code-block:: c++
+
+   molecular = Atom::ATOMIC;
+   mass_type = AtomVec::PER_TYPE;
+   if (atom->molecular == Atom::TEMPLATE)
+   if (atom->map_style == Atom::MAP_HASH)
+   atom->add_callback(Atom::GROW);
+   atom->delete_callback(id,Atom::RESTART);
+
 Simplify customized error messages
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

doc/src/Fortran.rst

@@ -315,6 +315,10 @@ of the contents of the :f:mod:`LIBLAMMPS` Fortran interface to LAMMPS.
    :ftype extract_variable: function
    :f set_variable: :f:subr:`set_variable`
    :ftype set_variable: subroutine
+   :f set_string_variable: :f:subr:`set_set_string_variable`
+   :ftype set_string_variable: subroutine
+   :f set_internal_variable: :f:subr:`set_internal_variable`
+   :ftype set_internal_variable: subroutine
    :f gather_atoms: :f:subr:`gather_atoms`
    :ftype gather_atoms: subroutine
    :f gather_atoms_concat: :f:subr:`gather_atoms_concat`
@@ -1398,7 +1402,28 @@ Procedures Bound to the :f:type:`lammps` Derived Type
 
    Set the value of a string-style variable.
 
-   .. versionadded:: 3Nov2022
+   .. deprecated:: TBD
+
+   This function assigns a new value from the string *str* to the string-style
+   variable *name*\ . If *name* does not exist or is not a string-style
+   variable, an error is generated.
+
+   .. warning::
+
+      This subroutine is deprecated and :f:subr:`set_string_variable`
+      should be used instead.
+
+   :p character(len=*) name: name of the variable
+   :p character(len=*) str:  new value to assign to the variable
+   :to: :cpp:func:`lammps_set_variable`
+
+--------
+
+.. f:subroutine:: set_string_variable(name, str)
+
+   Set the value of a string-style variable.
+
+   .. versionadded:: TBD
 
    This function assigns a new value from the string *str* to the string-style
    variable *name*\ . If *name* does not exist or is not a string-style
@@ -1406,7 +1431,23 @@ Procedures Bound to the :f:type:`lammps` Derived Type
 
    :p character(len=*) name: name of the variable
    :p character(len=*) str:  new value to assign to the variable
-   :to: :cpp:func:`lammps_set_variable`
+   :to: :cpp:func:`lammps_set_string_variable`
+
+--------
+
+.. f:subroutine:: set_internal_variable(name, val)
+
+   Set the value of a internal-style variable.
+
+   .. versionadded:: TBD
+
+   This function assigns a new value from the floating-point number *val* to
+   the internal-style variable *name*\ . If *name* does not exist or is not
+   an internal-style variable, an error is generated.
+
+   :p character(len=*) name: name of the variable
+   :p read(c_double) val:  new value to assign to the variable
+   :to: :cpp:func:`lammps_set_internal_variable`
 
 --------
 

doc/src/Howto_github.rst

@@ -480,11 +480,11 @@ Some recent changes to the workflow are not captured in this tutorial.
 For example, in addition to the *develop* branch, to which all new
 features should be submitted, there is also a *release*, a *stable*, and
 a *maintenance* branch; the *release* branch is updated from the
-*develop* as part of a feature release, and *stable* (together with
+*develop* branch as part of a "feature release", and *stable* (together
-*release*) are updated from *develop* when a stable release is made. In
+with *release*) are updated from *develop* when a "stable release" is
-between stable releases, selected bug fixes and infrastructure updates
+made. In between stable releases, selected bug fixes and infrastructure
-are back-ported from the *develop* branch to the *maintenance* branch
+updates are back-ported from the *develop* branch to the *maintenance*
-and occasionally merged to *stable* as an update release.
+branch and occasionally merged to *stable* as an update release.
 
 Furthermore, the naming of the release tags now follow the pattern
 "patch_<Day><Month><Year>" to simplify comparisons between releases.

doc/src/Install_git.rst

@@ -28,16 +28,16 @@ provides `limited support for subversion clients <svn_>`_.
 
 You can follow the LAMMPS development on 4 different git branches:
 
-* **release**  :  this branch is updated with every patch or feature release;
+* **develop** : this branch follows the ongoing development and is
+  updated with every merge commit of a pull request
+* **release** : this branch is updated with every "feature release";
    updates are always "fast-forward" merges from *develop*
-* **develop**  :  this branch follows the ongoing development and
-  is updated with every merge commit of a pull request
-* **stable**   :  this branch is updated from the *release* branch with
-  every stable release version and also has selected bug fixes with every
-  update release when the *maintenance* branch is merged into it
 * **maintenance** : this branch collects back-ported bug fixes from the
-  *develop* branch to the *stable* branch. It is used to update *stable*
+  *develop* branch to the *stable* branch.  It is used to update the
-  for update releases and it synchronized with *stable* at each stable release.
+  *stable* branch for "stable update releases".
+* **stable** : this branch is updated from the *release* branch with
+  every "stable release" version and also has selected bug fixes with
+  every "update release" when the *maintenance* branch is merged into it
 
 To access the git repositories on your box, use the clone command to
 create a local copy of the LAMMPS repository with a command like:

doc/src/Library_objects.rst

@@ -9,6 +9,8 @@ fixes, or variables in LAMMPS using the following functions:
 - :cpp:func:`lammps_extract_variable_datatype`
 - :cpp:func:`lammps_extract_variable`
 - :cpp:func:`lammps_set_variable`
+- :cpp:func:`lammps_set_string_variable`
+- :cpp:func:`lammps_set_internal_variable`
 - :cpp:func:`lammps_variable_info`
 
 -----------------------
@@ -38,6 +40,16 @@ fixes, or variables in LAMMPS using the following functions:
 
 -----------------------
 
+.. doxygenfunction:: lammps_set_string_variable
+   :project: progguide
+
+-----------------------
+
+.. doxygenfunction:: lammps_set_internal_variable
+   :project: progguide
+
+-----------------------
+
 .. doxygenfunction:: lammps_variable_info
    :project: progguide
 

doc/src/Manual_version.rst

@@ -3,45 +3,25 @@ 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, and we aim to keep it working
-correctly and reliably at all times.  You can follow its development
+correctly and reliably at all times.  Also, several variants of static
-in a public `git repository on GitHub <https://github.com/lammps/lammps>`_.
+code analysis are run regularly to maintain or improve the overall code
-
+quality, consistency, and compliance with programming standards, best
-Modifications of the LAMMPS source code (like bug fixes, code refactors,
+practices and style conventions.  You can follow its development in a
-updates to existing features, or addition of new features) are organized
+public `git repository on GitHub <https://github.com/lammps/lammps>`_.
-into pull requests.  Pull requests will be merged into the *develop*
-branch of the git repository after they pass automated testing and code
-review by the LAMMPS developers.  When a sufficient number of changes
-have accumulated *and* the *develop* branch version passes an extended
-set of automated tests, we release it as a *feature release*, which are
-currently made every 4 to 8 weeks.  The *release* branch of the git
-repository is updated with every such release.  A summary of the most
-important changes of the patch releases are on `this website page
-<https://www.lammps.org/bug.html>`_.  More detailed release notes are
-`available on GitHub <https://github.com/lammps/lammps/releases/>`_.
-
-Once or twice a year, we have a "stabilization period" where we apply
-only bug fixes and small, non-intrusive changes to the *develop*
-branch.  At the same time, the code is subjected to more detailed and
-thorough manual testing than the default automated testing.  Also,
-several variants of static code analysis are run to improve the overall
-code quality, consistency, and compliance with programming standards,
-best practices and style conventions.
-
-The release after such a stabilization period is called a *stable*
-version and both, the *release* and the *stable* branches are updated
-with it.  Between stable releases, we collect back-ported bug fixes and
-updates from the *develop* branch in the *maintenance* branch.  From the
-*maintenance* branch we make occasional update releases and update the
-*stable* branch accordingly.
 
 Each version of LAMMPS contains all the documented *features* up to and
 including its version date.  For recently added features, we add markers
 to the documentation at which specific LAMMPS version a feature or
 keyword was added or significantly changed.
 
-The version date is printed to the screen and log file every time you run
+Identifying the Version
-LAMMPS.  It is also in the file src/version.h and in the LAMMPS
+^^^^^^^^^^^^^^^^^^^^^^^
-directory name created when you unpack a tarball.  And it is on the
+
+The version date is printed to the screen and log file every time you
+run LAMMPS.  There also is an indication, if a LAMMPS binary was
+compiled from version with modifications **after** a release.
+It is also visible in the file src/version.h and in the LAMMPS directory
+name created when you unpack a downloaded tarball.  And it is on the
 first page of the :doc:`manual <Manual>`.
 
 * If you browse the HTML pages of the online version of the LAMMPS
@@ -53,3 +33,56 @@ first page of the :doc:`manual <Manual>`.
 * If you browse the HTML pages included in your downloaded tarball, they
   describe the version you have, which may be older than the online
   version.
+
+LAMMPS releases, branches, and tags
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. figure:: JPG/lammps-releases.png
+   :figclass: align-center
+
+   Relations between releases, main branches, and tags in the LAMMPS git repository
+
+Development
+"""""""""""
+
+Modifications of the LAMMPS source code (like bug fixes, code
+refactoring, updates to existing features, or addition of new features)
+are organized into pull requests.  Pull requests will be merged into the
+*develop* branch of the git repository after they pass automated testing
+and code review by the LAMMPS developers.
+
+Feature Releases
+""""""""""""""""
+
+When a sufficient number of new features and updates have accumulated
+*and* the LAMMPS version on the *develop* branch passes an extended set
+of automated tests, we release it as a *feature release*, which are
+currently made every 4 to 8 weeks.  The *release* branch of the git
+repository is updated with every such *feature release* and a tag in the
+format ``patch_1May2014`` is added.  A summary of the most important
+changes of these releases for the current year are posted on `this
+website page <https://www.lammps.org/bug.html>`_.  More detailed release
+notes are `available on GitHub
+<https://github.com/lammps/lammps/releases/>`_.
+
+Stable Releases
+"""""""""""""""
+
+About once a year, we release a *stable release* version of LAMMPS.
+This is done after a "stabilization period" where we apply only bug
+fixes and small, non-intrusive changes to the *develop* branch but no
+new features.  At the same time, the code is subjected to more detailed
+and thorough manual testing than the default automated testing.
+After such a *stable release*, both the *release* and the *stable*
+branches are updated and two tags are applied, a ``patch_1May2014`` format
+and a ``stable_1May2014`` format tag.
+
+Stable Release Updates
+""""""""""""""""""""""
+
+Between *stable releases*, we collect bug fixes and updates back-ported
+from the *develop* branch in a branch called *maintenance*.  From the
+*maintenance* branch we make occasional *stable update releases* and
+update the *stable* branch accordingly.  The first update to the
+``stable_1May2014`` release would be tagged as
+``stable_1May2014_update1``.  These updates contain no new features.

doc/src/angle_charmm.rst

@@ -70,7 +70,9 @@ for more info.
 Related commands
 """"""""""""""""
 
-:doc:`angle_coeff <angle_coeff>`
+:doc:`angle_coeff <angle_coeff>`, :doc:`pair_style lj/charmm variants <pair_charmm>`,
+:doc:`dihedral_style charmm <dihedral_charmm>`,
+:doc:`dihedral_style charmmfsw <dihedral_charmm>`, :doc:`fix cmap <fix_cmap>`
 
 Default
 """""""

doc/src/angle_lepton.rst

@@ -11,7 +11,16 @@ Syntax
 
 .. code-block:: LAMMPS
 
-   angle_style lepton
+   angle_style style args
+
+* style = *lepton*
+* args = optional arguments
+
+.. parsed-literal::
+
+   args = *auto_offset* or *no_offset*
+     *auto_offset* = offset the potential energy so that the value at theta0 is 0.0 (default)
+     *no_offset* = do not offset the potential energy
 
 Examples
 """"""""
@@ -19,6 +28,7 @@ Examples
 .. code-block:: LAMMPS
 
    angle_style lepton
+   angle_style lepton no_offset
 
    angle_coeff  1  120.0  "k*theta^2; k=250.0"
    angle_coeff  2   90.0  "k2*theta^2 + k3*theta^3 + k4*theta^4; k2=300.0; k3=-100.0; k4=50.0"
@@ -41,6 +51,13 @@ angle coefficient.  For example `"200.0*theta^2"` represents a
 
    U_{angle,i} = K (\theta_i - \theta_0)^2 = K \theta^2 \qquad \theta = \theta_i - \theta_0
 
+.. versionchanged:: TBD
+
+By default the potential energy U is shifted so that the value U is 0.0
+for $theta = theta_0$.  This is equivalent to using the optional keyword
+*auto_offset*.  When using the keyword *no_offset* instead, the
+potential energy is not shifted.
+
 The `Lepton library <https://simtk.org/projects/lepton>`_, that the
 *lepton* angle style interfaces with, evaluates this expression string
 at run time to compute the pairwise energy.  It also creates an

doc/src/atom_style.rst

@@ -49,248 +49,221 @@ Examples
 Description
 """""""""""
 
-Define what style of atoms to use in a simulation.  This determines
+The *atom_style* command selects which per-atom attributes are
-what attributes are associated with the atoms.  This command must be
+associated with atoms in a LAMMPS simulation and thus stored and
-used before a simulation is setup via a :doc:`read_data <read_data>`,
+communicated with those atoms as well as read from and stored in data
-:doc:`read_restart <read_restart>`, or :doc:`create_box <create_box>`
+and restart files.  Different models (e.g. :doc:`pair styles
-command.
+<pair_style>`) require access to specific per-atom attributes and thus
+require a specific atom style.  For example, to compute Coulomb
+interactions, the atom must have a "charge" (aka "q") attribute.
+
+A number of distinct atom styles exist that combine attributes.  Some
+atom styles are a superset of other atom styles.  Further attributes
+may be added to atoms either via using a hybrid style which provides a
+union of the attributes of the sub-styles, or via the :doc:`fix
+property/atom <fix_property_atom>` command.  The *atom_style* command
+must be used before a simulation is setup via a :doc:`read_data
+<read_data>`, :doc:`read_restart <read_restart>`, or :doc:`create_box
+<create_box>` command.
 
 .. note::
 
-   Many of the atom styles discussed here are only enabled if
+   Many of the atom styles discussed here are only enabled if LAMMPS was
-   LAMMPS was built with a specific package, as listed below in the
+   built with a specific package, as listed below in the Restrictions
-   Restrictions section.
+   section.
 
-Once a style is assigned, it cannot be changed, so use a style general
+Once a style is selected and the simulation box defined, it cannot be
-enough to encompass all attributes.  E.g. with style *bond*, angular
+changed but only augmented with the :doc:`fix property/atom
-terms cannot be used or added later to the model.  It is OK to use a
+<fix_property_atom>` command.  So one should select an atom style
-style more general than needed, though it may be slightly inefficient.
+general enough to encompass all attributes required.  E.g. with atom
+style *bond*, it is not possible to define angles and use angle styles.
 
-The choice of style affects what quantities are stored by each atom,
+It is OK to use a style more general than needed, though it may be
-what quantities are communicated between processors to enable forces
+slightly inefficient because it will allocate and communicate
-to be computed, and what quantities are listed in the data file read
+additional unused data.
-by the :doc:`read_data <read_data>` command.
 
-These are the additional attributes of each style and the typical
+Atom style attributes
-kinds of physical systems they are used to model.  All styles store
+"""""""""""""""""""""
-coordinates, velocities, atom IDs and types.  See the
+
+The atom style *atomic* has the minimum subset of per-atom attributes
+and is also the default setting.  It encompasses the following per-atom
+attributes (name of the vector or array in the :doc:`Atom class
+<Classes_atom>` is given in parenthesis): atom-ID (tag), type (type),
+position (x), velocities (v), forces (f), image flags (image), group
+membership (mask).  Since all atom styles are a superset of atom style
+*atomic*\ , they all include these attributes.
+
+This table lists all the available atom styles, which attributes they
+provide, which :doc:`package <Packages>` is required to use them, and
+what the typical applications are that use them.  See the
 :doc:`read_data <read_data>`, :doc:`create_atoms <create_atoms>`, and
-:doc:`set <set>` commands for info on how to set these various
+:doc:`set <set>` commands for details on how to set these various
-quantities.
+quantities.  More information about many of the styles is provided in
+the Additional Information section below.
 
-+--------------+-----------------------------------------------------+--------------------------------------+
+.. list-table::
-| *amoeba*     | molecular + charge + 1/5 neighbors                  | AMOEBA/HIPPO polarized force fields  |
+   :header-rows: 1
-+--------------+-----------------------------------------------------+--------------------------------------+
+   :widths: auto
-| *angle*      | bonds and angles                                    | bead-spring polymers with stiffness  |
+
-+--------------+-----------------------------------------------------+--------------------------------------+
+   * - Atom style
-| *atomic*     | only the default values                             | coarse-grain liquids, solids, metals |
+     - Attributes
-+--------------+-----------------------------------------------------+--------------------------------------+
+     - Required package
-| *body*       | mass, inertia moments, quaternion, angular momentum | arbitrary bodies                     |
+     - Applications
-+--------------+-----------------------------------------------------+--------------------------------------+
+   * - *amoeba*
-| *bond*       | bonds                                               | bead-spring polymers                 |
+     - *full* + "1-5 special neighbor data"
-+--------------+-----------------------------------------------------+--------------------------------------+
+     - :ref:`AMOEBA <PKG-AMOEBA>`
-| *charge*     | charge                                              | atomic system with charges           |
+     - AMOEBA/HIPPO force fields
-+--------------+-----------------------------------------------------+--------------------------------------+
+   * - *angle*
-| *dielectric* | normx normy normz area/patch ed em epsilon curv     | system with surface polarization     |
+     - *bond* + "angle data"
-+--------------+-----------------------------------------------------+--------------------------------------+
+     - :ref:`MOLECULE <PKG-MOLECULE>`
-| *dipole*     | charge and dipole moment                            | system with dipolar particles        |
+     - bead-spring polymers with stiffness
-+--------------+-----------------------------------------------------+--------------------------------------+
+   * - *atomic*
-| *dpd*        | internal temperature and internal energies          | DPD particles                        |
+     - tag, type, x, v, f, image, mask
-+--------------+-----------------------------------------------------+--------------------------------------+
+     -
-| *edpd*       | temperature and heat capacity                       | eDPD particles                       |
+     - atomic liquids, solids, metals
-+--------------+-----------------------------------------------------+--------------------------------------+
+   * - *body*
-| *electron*   | charge and spin and eradius                         | electronic force field               |
+     - *atomic* + radius, rmass, angmom, torque, body
-+--------------+-----------------------------------------------------+--------------------------------------+
+     - :ref:`BODY <PKG-BODY>`
-| *ellipsoid*  | shape, quaternion, angular momentum                 | aspherical particles                 |
+     - arbitrary bodies, see :doc:`body howto <Howto_body>`
-+--------------+-----------------------------------------------------+--------------------------------------+
+   * - *bond*
-| *full*       | molecular + charge                                  | bio-molecules                        |
+     - *atomic* + molecule, nspecial, special + "bond data"
-+--------------+-----------------------------------------------------+--------------------------------------+
+     - :ref:`MOLECULE <PKG-MOLECULE>`
-| *line*       | end points, angular velocity                        | rigid bodies                         |
+     - bead-spring polymers
-+--------------+-----------------------------------------------------+--------------------------------------+
+   * - *bpm/sphere*
-| *mdpd*       | density                                             | mDPD particles                       |
+     - *bond* + radius, rmass, omega, torque, quat
-+--------------+-----------------------------------------------------+--------------------------------------+
+     - :ref:`BPM <PKG-BPM>`
-| *molecular*  | bonds, angles, dihedrals, impropers                 | uncharged molecules                  |
+     - granular bonded particle models, see :doc:`BPM howto <Howto_bpm>`
-+--------------+-----------------------------------------------------+--------------------------------------+
+   * - *charge*
-| *oxdna*      | nucleotide polarity                                 | coarse-grained DNA and RNA models    |
+     - *atomic* + q
-+--------------+-----------------------------------------------------+--------------------------------------+
+     -
-| *peri*       | mass, volume                                        | mesoscopic Peridynamic models        |
+     - atomic systems with charges
-+--------------+-----------------------------------------------------+--------------------------------------+
+   * - *dielectric*
-| *smd*        | volume, kernel diameter, contact radius, mass       | solid and fluid SPH particles        |
+     - *full* + mu, area, ed, em, epsilon, curvature, q_scaled
-+--------------+-----------------------------------------------------+--------------------------------------+
+     - :ref:`DIELECTRIC <PKG-DIELECTRIC>`
-| *sph*        | rho, esph, cv                                       | SPH particles                        |
+     - systems with surface polarization
-+--------------+-----------------------------------------------------+--------------------------------------+
+   * - *dipole*
-| *sphere*     | diameter, mass, angular velocity                    | granular models                      |
+     - *charge* + mu
-+--------------+-----------------------------------------------------+--------------------------------------+
+     - :ref:`DIPOLE <PKG-DIPOLE>`
-| *bpm/sphere* | diameter, mass, angular velocity, quaternion        | granular bonded particle models (BPM)|
+     - atomic systems with charges and point dipoles
-+--------------+-----------------------------------------------------+--------------------------------------+
+   * - *dpd*
-| *spin*       | magnetic moment                                     | system with magnetic particles       |
+     - *atomic* + rho + "reactive DPD data"
-+--------------+-----------------------------------------------------+--------------------------------------+
+     - :ref:`DPD-REACT <PKG-DPD-REACT>`
-| *tdpd*       | chemical concentration                              | tDPD particles                       |
+     - reactive DPD
-+--------------+-----------------------------------------------------+--------------------------------------+
+   * - *edpd*
-| *template*   | template index, template atom                       | small molecules with fixed topology  |
+     - *atomic* + "eDPD data"
-+--------------+-----------------------------------------------------+--------------------------------------+
+     - :ref:`DPD-MESO <PKG-DPD-MESO>`
-| *tri*        | corner points, angular momentum                     | rigid bodies                         |
+     - Energy conservative DPD (eDPD)
-+--------------+-----------------------------------------------------+--------------------------------------+
+   * - *electron*
-| *wavepacket* | charge, spin, eradius, etag, cs_re, cs_im           | AWPMD                                |
+     - *charge* + espin, eradius, ervel, erforce
-+--------------+-----------------------------------------------------+--------------------------------------+
+     - :ref:`EFF <PKG-EFF>`
+     - Electron force field systems
+   * - *ellipsoid*
+     - *atomic* + rmass, angmom, torque, ellipsoid
+     -
+     - aspherical particles
+   * - *full*
+     - *molecular* + q
+     - :ref:`MOLECULE <PKG-MOLECULE>`
+     - molecular force fields
+   * - *line*
+     - *atomic* + molecule, radius, rmass, omega, torque, line
+     -
+     - 2-d rigid body particles
+   * - *mdpd*
+     - *atomic* + rho, drho, vest
+     - :ref:`DPD-MESO <PKG-DPD-MESO>`
+     - Many-body DPD (mDPD)
+   * - *molecular*
+     - *angle* + "dihedral and improper data"
+     - :ref:`MOLECULE <PKG-MOLECULE>`
+     - apolar and uncharged molecules
+   * - *oxdna*
+     - *atomic* + id5p
+     - :ref:`CG-DNA <PKG-CG-DNA>`
+     - coarse-grained DNA and RNA models
+   * - *peri*
+     - *atomic* + rmass, vfrac, s0, x0
+     - :ref:`PERI <PKG-PERI>`
+     - mesoscopic Peridynamics models
+   * - *smd*
+     - *atomic* + molecule, radius, rmass + "smd data"
+     - :ref:`MACHDYN <PKG-MACHDYN>`
+     - Smooth Mach Dynamics models
+   * - *sph*
+     - *atomic* + "sph data"
+     - :ref:`SPH <PKG-SPH>`
+     - Smoothed particle hydrodynamics models
+   * - *sphere*
+     - *atomic* + radius, rmass, omega, torque
+     -
+     - finite size spherical particles, e.g. granular models
+   * - *spin*
+     - *atomic* + "magnetic moment data"
+     - :ref:`SPIN <PKG-SPIN>`
+     - magnetic particles
+   * - *tdpd*
+     - *atomic* + cc, cc_flux, vest
+     - :ref:`DPD-MESO <PKG-DPD-MESO>`
+     - Transport DPD (tDPD)
+   * - *template*
+     - *atomic* + molecule, molindex, molatom
+     - :ref:`MOLECULE <PKG-MOLECULE>`
+     - molecular systems where attributes are taken from :doc:`molecule files <molecule>`
+   * - *tri*
+     - *sphere* + molecule, angmom, tri
+     -
+     - 3-d triangulated rigid body LJ particles
+   * - *wavepacket*
+     - *charge* + "wavepacket data"
+     - :ref:`AWPMD <PKG-AWPMD>`
+     - Antisymmetrized wave packet MD
 
 .. note::
 
-   It is possible to add some attributes, such as a molecule ID, to
+   It is possible to add some attributes, such as a molecule ID and
-   atom styles that do not have them via the :doc:`fix property/atom
+   charge, to atom styles that do not have them built in using the
-   <fix_property_atom>` command.  This command also allows new custom
+   :doc:`fix property/atom <fix_property_atom>` command.  This command
-   attributes consisting of extra integer or floating-point values to
+   also allows new custom-named attributes consisting of extra integer
-   be added to atoms.  See the :doc:`fix property/atom
+   or floating-point values or vectors to be added to atoms.  See the
-   <fix_property_atom>` page for examples of cases where this is
+   :doc:`fix property/atom <fix_property_atom>` page for examples of
-   useful and details on how to initialize, access, and output the
+   cases where this is useful and details on how to initialize,
-   custom values.
+   access, and output these custom values.
 
-All of the above styles define point particles, except the *sphere*,
+----------
-*bpm/sphere*, *ellipsoid*, *electron*, *peri*, *wavepacket*, *line*,
-*tri*, and *body* styles, which define finite-size particles.  See the
-:doc:`Howto spherical <Howto_spherical>` page for an overview of using
-finite-size particle models with LAMMPS.
 
-All of the point-particle styles assign mass to particles on a
+Particle size and mass
-per-type basis, using the :doc:`mass <mass>` command, The finite-size
+""""""""""""""""""""""
-particle styles assign mass to individual particles on a per-particle
-basis.
 
-For the *sphere* and *bpm/sphere* styles, the particles are spheres
+All of the atom styles define point particles unless they (1) define
-and each stores a per-particle diameter and mass.  If the diameter >
+finite-size spherical particles via the *radius* attribute, or (2)
-0.0, the particle is a finite-size sphere.  If the diameter = 0.0, it
+define finite-size aspherical particles (e.g. the *body*, *ellipsoid*,
-is a point particle.  Note that by use of the *disc* keyword with the
+*line*, and *tri* styles).  Most of these styles can also be used with
-:doc:`fix nve/sphere <fix_nve_sphere>`, :doc:`fix nvt/sphere
+mixtures of point and finite-size particles.
-<fix_nvt_sphere>`, :doc:`fix nph/sphere <fix_nph_sphere>`,
-:doc:`fix npt/sphere <fix_npt_sphere>` commands for the *sphere* style,
-spheres can be effectively treated as 2d discs for a 2d simulation if
-desired.  See also the :doc:`set density/disc <set>` command.  These
-styles take an optional 0 or 1 argument.  A value of 0 means the
-radius of each sphere is constant for the duration of the simulation.
-A value of 1 means the radii may vary dynamically during the simulation,
-e.g. due to use of the :doc:`fix adapt <fix_adapt>` command.
 
-For the *ellipsoid* style, the particles are ellipsoids and each
+Note that the *radius* property may need to be provided as a
-stores a flag which indicates whether it is a finite-size ellipsoid or
+*diameter* (e.g. in :doc:`molecule files <molecule>` or :doc:`data
-a point particle.  If it is an ellipsoid, it also stores a shape
+files <read_data>`).  See the :doc:`Howto spherical <Howto_spherical>`
-vector with the 3 diameters of the ellipsoid and a quaternion 4-vector
+page for an overview of using finite-size spherical and aspherical
-with its orientation.
+particle models with LAMMPS.
 
-For the *dielectric* style, each particle can be either a physical
+Unless an atom style defines the per-atom *rmass* attribute, particle
-particle (e.g. an ion), or an interface particle representing a boundary
+masses are defined on a per-type basis, using the :doc:`mass <mass>`
-element between two regions of different dielectric constant. For
+command.  This means each particle's mass is indexed by its atom
-interface particles, in addition to the properties associated with
+*type*.
-atom_style full, each particle also should be assigned a normal unit
-vector (defined by normx, normy, normz), an area (area/patch), the
-difference and mean of the dielectric constants of two sides of the
-interface along the direction of the normal vector (ed and em), the
-local dielectric constant at the boundary element (epsilon), and a mean
-local curvature (curv).  Physical particles must be assigned these
-values, as well, but only their local dielectric constants will be used;
-see documentation for associated :doc:`pair styles <pair_dielectric>`
-and :doc:`fixes <fix_polarize>`.  The distinction between the physical
-and interface particles is only meaningful when :doc:`fix polarize
-<fix_polarize>` commands are applied to the interface particles. This
-style is part of the DIELECTRIC package.
 
-For the *dipole* style, a point dipole is defined for each point
+A few styles define the per-atom *rmass* attribute which can also be
-particle.  Note that if you wish the particles to be finite-size
+added using the :doc:`fix property/atom <fix_property_atom>` command.
-spheres as in a Stockmayer potential for a dipolar fluid, so that the
+In this case each particle stores its own mass.  Atom styles that have
-particles can rotate due to dipole-dipole interactions, then you need
+a per-atom rmass may define it indirectly through setting particle
-to use atom_style hybrid sphere dipole, which will assign both a
+diameter and density on a per-particle basis.  If both per-type mass
-diameter and dipole moment to each particle.
+and per-atom *rmass* are defined (e.g. in a hybrid style), the
+per-atom mass will take precedence in any operation which which works
+with both flavors of mass.
 
-For the *electron* style, the particles representing electrons are 3d
+----------
-Gaussians with a specified position and bandwidth or uncertainty in
-position, which is represented by the eradius = electron size.
 
-For the *peri* style, the particles are spherical and each stores a
+Additional information about specific atom styles
-per-particle mass and volume.
+"""""""""""""""""""""""""""""""""""""""""""""""""
-
-The *bpm/sphere* style is part of the BPM package.
-
-The *oxdna* style is for coarse-grained nucleotides and stores the
-3'-to-5' polarity of the nucleotide strand, which is set through
-the bond topology in the data file. The first (second) atom in a
-bond definition is understood to point towards the 3'-end (5'-end)
-of the strand. Note that this style is part of the CG-DNA package.
-
-The *dpd* style is for dissipative particle dynamics (DPD) particles.
-Note that it is part of the DPD-REACT package, and is not for use with
-the :doc:`pair_style dpd or dpd/stat <pair_dpd>` commands, which can
-simply use atom_style atomic.  Atom_style dpd extends DPD particle
-properties with internal temperature (dpdTheta), internal conductive
-energy (uCond), internal mechanical energy (uMech), and internal
-chemical energy (uChem).
-
-The *edpd* style is for energy-conserving dissipative particle
-dynamics (eDPD) particles which store a temperature (edpd_temp), and
-heat capacity(edpd_cv).
-
-The *mdpd* style is for many-body dissipative particle dynamics (mDPD)
-particles which store a density (rho) for considering
-density-dependent many-body interactions.
-
-The *tdpd* style is for transport dissipative particle dynamics (tDPD)
-particles which store a set of chemical concentration. An integer
-"cc_species" is required to specify the number of chemical species
-involved in a tDPD system.
-
-The *sph* style is for smoothed particle hydrodynamics (SPH)
-particles which store a density (rho), energy (esph), and heat capacity
-(cv).
-
-The *smd* style is for a general formulation of Smooth Particle
-Hydrodynamics.  Both fluids and solids can be modeled.  Particles
-store the mass and volume of an integration point, a kernel diameter
-used for calculating the field variables (e.g. stress and deformation)
-and a contact radius for calculating repulsive forces which prevent
-individual physical bodies from penetrating each other.
-
-For the *spin* style, a magnetic spin is associated to each atom.
-Those spins have a norm (their magnetic moment) and a direction.
-
-The *wavepacket* style is similar to *electron*, but the electrons may
-consist of several Gaussian wave packets, summed up with coefficients
-cs= (cs_re,cs_im).  Each of the wave packets is treated as a separate
-particle in LAMMPS, wave packets belonging to the same electron must
-have identical *etag* values.
-
-For the *line* style, the particles are idealized line segments and
-each stores a per-particle mass and length and orientation (i.e. the
-end points of the line segment).
-
-For the *tri* style, the particles are planar triangles and each
-stores a per-particle mass and size and orientation (i.e. the corner
-points of the triangle).
-
-The *template* style allows molecular topology (bonds,angles,etc) to be
-defined via a molecule template using the :doc:`molecule <molecule>`
-command.  The template stores one or more molecules with a single copy
-of the topology info (bonds,angles,etc) of each.  Individual atoms
-only store a template index and template atom to identify which
-molecule and which atom-within-the-molecule they represent.  Using the
-*template* style instead of the *bond*, *angle*, *molecular* styles
-can save memory for systems comprised of a large number of small
-molecules, all of a single type (or small number of types).  See the
-paper by Grime and Voth, in :ref:`(Grime) <Grime>`, for examples of how this
-can be advantageous for large-scale coarse-grained systems.
-The ``examples/template`` directory has a few demo inputs and examples
-showing the use of the *template* atom style versus *molecular*.
-
-.. note::
-
-   When using the *template* style with a :doc:`molecule template
-   <molecule>` that contains multiple molecules, you should ensure the
-   atom types, bond types, angle_types, etc in all the molecules are
-   consistent.  E.g. if one molecule represents H2O and another CO2,
-   then you probably do not want each molecule file to define 2 atom
-   types and a single bond type, because they will conflict with each
-   other when a mixture system of H2O and CO2 molecules is defined,
-   e.g. by the :doc:`read_data <read_data>` command.  Rather the H2O
-   molecule should define atom types 1 and 2, and bond type 1.  And
-   the CO2 molecule should define atom types 3 and 4 (or atom types 3
-   and 2 if a single oxygen type is desired), and bond type 2.
 
 For the *body* style, the particles are arbitrary bodies with internal
 attributes defined by the "style" of the bodies, which is specified by
@@ -309,6 +282,148 @@ Note that there may be additional arguments required along with the
 *bstyle* specification, in the atom_style body command.  These
 arguments are described on the :doc:`Howto body <Howto_body>` doc page.
 
+For the *dielectric* style, each particle can be either a physical
+particle (e.g. an ion), or an interface particle representing a boundary
+element between two regions of different dielectric constant. For
+interface particles, in addition to the properties associated with
+atom_style full, each particle also should be assigned a normal unit
+vector (defined by normx, normy, normz), an area (area/patch), the
+difference and mean of the dielectric constants of two sides of the
+interface along the direction of the normal vector (ed and em), the
+local dielectric constant at the boundary element (epsilon), and a mean
+local curvature (curv).  Physical particles must be assigned these
+values, as well, but only their local dielectric constants will be used;
+see documentation for associated :doc:`pair styles <pair_dielectric>`
+and :doc:`fixes <fix_polarize>`.  The distinction between the physical
+and interface particles is only meaningful when :doc:`fix polarize
+<fix_polarize>` commands are applied to the interface particles. This
+style is part of the DIELECTRIC package.
+
+For the *dipole* style, a point dipole vector mu is defined for each
+point particle.  Note that if you wish the particles to be finite-size
+spheres as in a Stockmayer potential for a dipolar fluid, so that the
+particles can rotate due to dipole-dipole interactions, then you need
+to use the command `atom_style hybrid sphere dipole`, which will
+assign both a diameter and dipole moment to each particle.  This also
+requires using an integrator with a "/sphere" suffix like :doc:`fix
+nve/sphere <fix_nve_sphere>` or :doc:`fix nvt/sphere <fix_nvt_sphere>`
+and the "update dipole" or "update dlm" parameters to the fix
+commands.
+
+The *dpd* style is for reactive dissipative particle dynamics (DPD)
+particles.  Note that it is part of the DPD-REACT package, and is not
+required for use with the :doc:`pair_style dpd or dpd/stat <pair_dpd>`
+commands, which only require the attributes from atom_style *atomic*.
+Atom_style *dpd* extends DPD particle properties with internal
+temperature (dpdTheta), internal conductive energy (uCond), internal
+mechanical energy (uMech), and internal chemical energy (uChem).
+
+The *edpd* style is for energy-conserving dissipative particle
+dynamics (eDPD) particles which store a temperature (edpd_temp), and
+heat capacity (edpd_cv).
+
+For the *electron* style, the particles representing electrons are 3d
+Gaussians with a specified position and bandwidth or uncertainty in
+position, which is represented by the eradius = electron size.
+
+For the *ellipsoid* style, particles can be ellipsoids which each
+stores a shape vector with the 3 diameters of the ellipsoid and a
+quaternion 4-vector with its orientation.  Each particle stores a flag
+in the ellipsoid vector which indicates whether it is an ellipsoid (1)
+or a point particle (0).
+
+For the *line* style, particles can be are idealized line segments
+which store a per-particle mass and length and orientation (i.e. the
+end points of the line segment).  Each particle stores a flag in the
+line vector which indicates whether it is a line segment (1) or a
+point particle (0).
+
+The *mdpd* style is for many-body dissipative particle dynamics (mDPD)
+particles which store a density (rho) for considering density-dependent
+many-body interactions.
+
+The *oxdna* style is for coarse-grained nucleotides and stores the
+3'-to-5' polarity of the nucleotide strand, which is set through
+the bond topology in the data file. The first (second) atom in a
+bond definition is understood to point towards the 3'-end (5'-end)
+of the strand.
+
+For the *peri* style, the particles are spherical and each stores a
+per-particle mass and volume.
+
+The *smd* style is for Smooth Particle Mach dynamics.  Both fluids and
+solids can be modeled.  Particles store the mass and volume of an
+integration point, a kernel diameter used for calculating the field
+variables (e.g. stress and deformation) and a contact radius for
+calculating repulsive forces which prevent individual physical bodies
+from penetrating each other.
+
+The *sph* style is for smoothed particle hydrodynamics (SPH) particles
+which store a density (rho), energy (esph), and heat capacity (cv).
+
+For the *spin* style, a magnetic spin is associated with each atom.
+Those spins have a norm (their magnetic moment) and a direction.
+
+The *tdpd* style is for transport dissipative particle dynamics (tDPD)
+particles which store a set of chemical concentration. An integer
+"cc_species" is required to specify the number of chemical species
+involved in a tDPD system.
+
+The *wavepacket* style is similar to the *electron* style, but the
+electrons may consist of several Gaussian wave packets, summed up with
+coefficients cs= (cs_re,cs_im).  Each of the wave packets is treated
+as a separate particle in LAMMPS, wave packets belonging to the same
+electron must have identical *etag* values.
+
+The *sphere* and *bpm/sphere* styles allow particles to be either point
+particles or finite-size particles.  If the *radius* attribute is >
+0.0, the particle is a finite-size sphere.  If the diameter = 0.0, it
+is a point particle.  Note that by using the *disc* keyword with the
+:doc:`fix nve/sphere <fix_nve_sphere>`, :doc:`fix nvt/sphere
+<fix_nvt_sphere>`, :doc:`fix nph/sphere <fix_nph_sphere>`, :doc:`fix
+npt/sphere <fix_npt_sphere>` commands for the *sphere* style, spheres
+can be effectively treated as 2d discs for a 2d simulation if desired.
+See also the :doc:`set density/disc <set>` command.  These styles also
+take an optional 0 or 1 argument.  A value of 0 means the radius of
+each sphere is constant for the duration of the simulation (this is
+the default).  A value of 1 means the radii may vary dynamically
+during the simulation, e.g. due to use of the :doc:`fix adapt
+<fix_adapt>` command.
+
+The *template* style allows molecular topology (bonds,angles,etc) to be
+defined via a molecule template using the :doc:`molecule <molecule>`
+command.  The template stores one or more molecules with a single copy
+of the topology info (bonds,angles,etc) of each.  Individual atoms only
+store a template index and template atom to identify which molecule and
+which atom-within-the-molecule they represent.  Using the *template*
+style instead of the *bond*, *angle*, *molecular* styles can save memory
+for systems comprised of a large number of small molecules, all of a
+single type (or small number of types).  See the paper by Grime and
+Voth, in :ref:`(Grime) <Grime>`, for examples of how this can be
+advantageous for large-scale coarse-grained systems.  The
+``examples/template`` directory has a few demo inputs and examples
+showing the use of the *template* atom style versus *molecular*.
+
+.. note::
+
+   When using the *template* style with a :doc:`molecule template
+   <molecule>` that contains multiple molecules, you should ensure the
+   atom types, bond types, angle_types, etc in all the molecules are
+   consistent.  E.g. if one molecule represents H2O and another CO2,
+   then you probably do not want each molecule file to define 2 atom
+   types and a single bond type, because they will conflict with each
+   other when a mixture system of H2O and CO2 molecules is defined,
+   e.g. by the :doc:`read_data <read_data>` command.  Rather the H2O
+   molecule should define atom types 1 and 2, and bond type 1.  And
+   the CO2 molecule should define atom types 3 and 4 (or atom types 3
+   and 2 if a single oxygen type is desired), and bond type 2.
+
+For the *tri* style, particles can be planar triangles which each
+stores a per-particle mass and size and orientation (i.e. the corner
+points of the triangle).  Each particle stores a flag in the tri
+vector which indicates whether it is a triangle (1) or a point
+particle (0).
+
 ----------
 
 Typically, simulations require only a single (non-hybrid) atom style.
@@ -326,11 +441,12 @@ dipole".  When a hybrid style is used, atoms store and communicate the
 union of all quantities implied by the individual styles.
 
 When using the *hybrid* style, you cannot combine the *template* style
-with another molecular style that stores bond,angle,etc info on a
+with another molecular style that stores bond, angle, etc info on a
 per-atom basis.
 
-LAMMPS can be extended with new atom styles as well as new body
+LAMMPS can be extended with new atom styles as well as new body styles;
-styles; see the :doc:`Modify <Modify>` doc page.
+see the corresponding manual page on :doc:`modifying & extending LAMMPS
+<Modify_atom>`.
 
 ----------
 
@@ -346,54 +462,20 @@ This command cannot be used after the simulation box is defined by a
 
 Many of the styles listed above are only enabled if LAMMPS was built
 with a specific package, as listed below.  See the :doc:`Build package
-<Build_package>` page for more info.
+<Build_package>` page for more info.  The table above lists which package
-
+is required for individual atom styles.
-The *amoeba* style is part of the AMOEBA package.
-
-The *angle*, *bond*, *full*, *molecular*, and *template* styles are
-part of the MOLECULE package.
-
-The *line* and *tri* styles are part of the ASPHERE package.
-
-The *body* style is part of the BODY package.
-
-The *dipole* style is part of the DIPOLE package.
-
-The *peri* style is part of the PERI package for Peridynamics.
-
-The *oxdna* style is part of the CG-DNA package for coarse-grained
-simulation of DNA and RNA.
-
-The *electron* style is part of the EFF package for :doc:`electronic
-force fields <pair_eff>`.
-
-The *dpd* style is part of the DPD-REACT package for dissipative
-particle dynamics (DPD).
-
-The *edpd*, *mdpd*, and *tdpd* styles are part of the DPD-MESO package
-for energy-conserving dissipative particle dynamics (eDPD), many-body
-dissipative particle dynamics (mDPD), and transport dissipative particle
-dynamics (tDPD), respectively.
-
-The *sph* style is part of the SPH package for smoothed particle
-hydrodynamics (SPH).  See `this PDF guide
-<PDF/SPH_LAMMPS_userguide.pdf>`_ to using SPH in LAMMPS.
-
-The *spin* style is part of the SPIN package.
-
-The *wavepacket* style is part of the AWPMD package for the
-:doc:`antisymmetrized wave packet MD method <pair_awpmd>`.
 
 Related commands
 """"""""""""""""
 
-:doc:`read_data <read_data>`, :doc:`pair_style <pair_style>`
+:doc:`read_data <read_data>`, :doc:`pair_style <pair_style>`,
+:doc:`fix property/atom <fix_property_atom>`, :doc:`set <set>`
 
 Default
 """""""
 
-The default atom style is atomic.  If atom_style sphere is used its
+The default atom style is *atomic*.  If atom_style *sphere* or
-default argument is 0.
+*bpm/sphere* is used, its default argument is 0.
 
 ----------
 

doc/src/bond_bpm_rotational.rst

@@ -147,8 +147,8 @@ By default, pair forces are not calculated between bonded particles.
 Pair forces can alternatively be overlaid on top of bond forces by setting
 the *overlay/pair* keyword to *yes*. These settings require specific
 :doc:`special_bonds <special_bonds>` settings described in the
-restrictions.  Further details can be found in the :doc:`how to
+restrictions.  Further details can be found in the :doc:`how to <Howto_bpm>`
-<Howto_bpm>` page on BPMs.
+page on BPMs.
 
 .. versionadded:: 28Mar2023
 

doc/src/bond_bpm_spring.rst

@@ -113,8 +113,8 @@ By default, pair forces are not calculated between bonded particles.
 Pair forces can alternatively be overlaid on top of bond forces by setting
 the *overlay/pair* keyword to *yes*. These settings require specific
 :doc:`special_bonds <special_bonds>` settings described in the
-restrictions.  Further details can be found in the :doc:`how to
+restrictions.  Further details can be found in the :doc:`how to <Howto_bpm>`
-<Howto_bpm>` page on BPMs.
+page on BPMs.
 
 .. versionadded:: 28Mar2023
 

doc/src/bond_lepton.rst

@@ -11,7 +11,16 @@ Syntax
 
 .. code-block:: LAMMPS
 
-   bond_style lepton
+   bond_style style args
+
+* style = *lepton*
+* args = optional arguments
+
+.. parsed-literal::
+
+   args = *auto_offset* or *no_offset*
+     *auto_offset* = offset the potential energy so that the value at r0 is 0.0 (default)
+     *no_offset* = do not offset the potential energy
 
 Examples
 """"""""
@@ -19,6 +28,7 @@ Examples
 .. code-block:: LAMMPS
 
    bond_style lepton
+   bond_style lepton no_offset
 
    bond_coeff  1  1.5 "k*r^2; k=250.0"
    bond_coeff  2  1.1 "k2*r^2 + k3*r^3 + k4*r^4; k2=300.0; k3=-100.0; k4=50.0"
@@ -40,6 +50,13 @@ constant *K* of 200.0 energy units:
 
    U_{bond,i} = K (r_i - r_0)^2 = K r^2 \qquad r = r_i - r_0
 
+.. versionchanged:: TBD
+
+By default the potential energy U is shifted so that he value U is 0.0
+for $r = r_0$.  This is equivalent to using the optional keyword
+*auto_offset*.  When using the keyword *no_offset* instead, the
+potential energy is not shifted.
+
 The `Lepton library <https://simtk.org/projects/lepton>`_, that the
 *lepton* bond style interfaces with, evaluates this expression string at
 run time to compute the pairwise energy.  It also creates an analytical

doc/src/dihedral_charmm.rst

@@ -3,6 +3,7 @@
 .. index:: dihedral_style charmm/kk
 .. index:: dihedral_style charmm/omp
 .. index:: dihedral_style charmmfsw
+.. index:: dihedral_style charmmfsw/kk
 
 dihedral_style charmm command
 =============================
@@ -12,6 +13,8 @@ Accelerator Variants: *charmm/intel*, *charmm/kk*, *charmm/omp*
 dihedral_style charmmfsw command
 ================================
 
+Accelerator Variants: *charmmfsw/kk*
+
 Syntax
 """"""
 
@@ -144,7 +147,9 @@ for more info.
 Related commands
 """"""""""""""""
 
-:doc:`dihedral_coeff <dihedral_coeff>`
+:doc:`dihedral_coeff <dihedral_coeff>`,
+:doc:`pair_style lj/charmm variants <pair_charmm>`,
+:doc:`angle_style charmm <angle_charmm>`, :doc:`fix cmap <fix_cmap>`
 
 Default
 """""""

doc/src/fix_halt.rst

@@ -183,4 +183,4 @@ Related commands
 Default
 """""""
 
-The option defaults are error = hard, message = yes, and path = ".".
+The option defaults are error = soft, message = yes, and path = ".".

doc/src/fix_neb.rst

@@ -109,7 +109,7 @@ Note that in this case the specified *Kspring* is in
 force/distance units.
 
 With a value of *ideal*, the spring force is computed as suggested in
-ref`(WeinanE) <WeinanE>`
+:ref:`(WeinanE) <WeinanE>`
 
 .. math::
 
@@ -120,18 +120,18 @@ and :math:`RD_{ideal}` is the ideal *RD* for which all the images are
 equally spaced.  I.e. :math:`RD_{ideal} = (i-1) \cdot meanDist` when the
 climbing replica is off, where *i* is the replica number).  The
 *meanDist* is the average distance between replicas.  Note that in this
-case the specified *Kspring* is in force units. When the climbing replica
+case the specified *Kspring* is in force units. When the climbing
-is on, :math:`RD_{ideal}` and :math:`meanDist` are calculated separately
+replica is on, :math:`RD_{ideal}` and :math:`meanDist` are calculated
-each side of the climbing image. Note that the *ideal* form of nudging
+separately each side of the climbing image. Note that the *ideal* form
-can often be more effective at keeping the replicas equally spaced before
+of nudging can often be more effective at keeping the replicas equally
-climbing, then equally spaced either side of the climbing image whilst
+spaced before climbing, then equally spaced either side of the climbing
-climbing.
+image whilst climbing.
 
-With a value of *equal* the spring force is computed as for *ideal*
+With a value of *equal* the spring force is computed as for *ideal* when
-when the climbing replica is off, promoting equidistance. When the climbing
+the climbing replica is off, promoting equidistance. When the climbing
 replica is on, the spring force is computed to promote equidistant
-absolute differences in energy, rather than distance, each side of
+absolute differences in energy, rather than distance, each side of the
-the climbing image:
+climbing image:
 
 .. math::
 
@@ -143,23 +143,22 @@ where *ED* is the cumulative sum of absolute energy differences:
 
    ED = \sum_{i<N} \left|E(R_{i+1}) - E(R_i)\right|,
 
-*meanEdist* is the average absolute energy difference between
+*meanEdist* is the average absolute energy difference between replicas
-replicas up to the climbing image or from the climbing image
+up to the climbing image or from the climbing image to the final image,
-to the final image, for images before or after the climbing
+for images before or after the climbing image
-image respectively. :math:`ED_{ideal}` is the corresponding
+respectively. :math:`ED_{ideal}` is the corresponding cumulative sum of
-cumulative sum of average absolute energy differences in
+average absolute energy differences in each case, in close analogy to
-each case, in close analogy to *ideal*. This form of nudging
+*ideal*. This form of nudging is to aid schemes which integrate forces
-is to aid schemes which integrate forces along, or near to,
+along, or near to, NEB pathways such as :doc:`fix_pafi <fix_pafi>`.
-NEB pathways such as :doc:`fix_pafi <fix_pafi>`.
 
 ----------
 
-The keyword *perp* specifies if and how a perpendicular nudging force
+The keyword *perp* specifies if and how a perpendicular nudging force is
-is computed.  It adds a spring force perpendicular to the path in
+computed.  It adds a spring force perpendicular to the path in order to
-order to prevent the path from becoming too strongly kinked.  It can
+prevent the path from becoming too strongly kinked.  It can
 significantly improve the convergence of the NEB calculation when the
-resolution is poor.  I.e. when few replicas are used; see
+resolution is poor.  I.e. when few replicas are used; see :ref:`(Maras)
-:ref:`(Maras) <Maras1>` for details.
+<Maras1>` for details.
 
 The perpendicular spring force is given by
 
@@ -181,10 +180,10 @@ force is added.
 
 By default, no additional forces act on the first and last replicas
 during the NEB relaxation, so these replicas simply relax toward their
-respective local minima.  By using the key word *end*, additional
+respective local minima.  By using the key word *end*, additional forces
-forces can be applied to the first and/or last replicas, to enable
+can be applied to the first and/or last replicas, to enable them to
-them to relax toward a MEP while constraining their energy E to the
+relax toward a MEP while constraining their energy E to the target
-target energy ETarget.
+energy ETarget.
 
 If :math:`E_{Target} > E`, the interatomic force :math:`F_i` for the
 specified replica becomes:
@@ -197,33 +196,33 @@ specified replica becomes:
 The "spring" constant on the difference in energies is the specified
 *Kspring3* value.
 
-When *estyle* is specified as *first*, the force is applied to the
+When *estyle* is specified as *first*, the force is applied to the first
-first replica.  When *estyle* is specified as *last*, the force is
+replica.  When *estyle* is specified as *last*, the force is applied to
-applied to the last replica.  Note that the *end* keyword can be used
+the last replica.  Note that the *end* keyword can be used twice to add
-twice to add forces to both the first and last replicas.
+forces to both the first and last replicas.
 
 For both these *estyle* settings, the target energy *ETarget* is set
 to the initial energy of the replica (at the start of the NEB
 calculation).
 
 If the *estyle* is specified as *last/efirst* or *last/efirst/middle*,
-force is applied to the last replica, but the target energy *ETarget*
+force is applied to the last replica, but the target energy *ETarget* is
-is continuously set to the energy of the first replica, as it evolves
+continuously set to the energy of the first replica, as it evolves
 during the NEB relaxation.
 
 The difference between these two *estyle* options is as follows.  When
 *estyle* is specified as *last/efirst*, no change is made to the
-inter-replica force applied to the intermediate replicas (neither
+inter-replica force applied to the intermediate replicas (neither first
-first or last).  If the initial path is too far from the MEP, an
+or last).  If the initial path is too far from the MEP, an intermediate
-intermediate replica may relax "faster" and reach a lower energy than
+replica may relax "faster" and reach a lower energy than the last
-the last replica.  In this case the intermediate replica will be
+replica.  In this case the intermediate replica will be relaxing toward
-relaxing toward its own local minima.  This behavior can be prevented
+its own local minima.  This behavior can be prevented by specifying
-by specifying *estyle* as *last/efirst/middle* which will alter the
+*estyle* as *last/efirst/middle* which will alter the inter-replica
-inter-replica force applied to intermediate replicas by removing the
+force applied to intermediate replicas by removing the contribution of
-contribution of the gradient to the inter-replica force.  This will
+the gradient to the inter-replica force.  This will only be done if a
-only be done if a particular intermediate replica has a lower energy
+particular intermediate replica has a lower energy than the first
-than the first replica.  This should effectively prevent the
+replica.  This should effectively prevent the intermediate replicas from
-intermediate replicas from over-relaxing.
+over-relaxing.
 
 After converging a NEB calculation using an *estyle* of
 *last/efirst/middle*, you should check that all intermediate replicas
@@ -237,9 +236,10 @@ target energy.
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
-No information about this fix is written to :doc:`binary restart files <restart>`.  None of the :doc:`fix_modify <fix_modify>` options
+No information about this fix is written to :doc:`binary restart files
-are relevant to this fix.  No global or per-atom quantities are stored
+<restart>`.  None of the :doc:`fix_modify <fix_modify>` options are
-by this fix for access by various :doc:`output commands <Howto_output>`.
+relevant to this fix.  No global or per-atom quantities are stored by
+this fix for access by various :doc:`output commands <Howto_output>`.
 No parameter of this fix can be used with the *start/stop* keywords of
 the :doc:`run <run>` command.
 

doc/src/fix_qeq.rst

@@ -232,8 +232,6 @@ These fixes are part of the QEQ package.  They are only enabled if
 LAMMPS was built with that package.  See the :doc:`Build package
 <Build_package>` page for more info.
 
-These qeq fixes are not compatible with the GPU and USER-INTEL packages.
-
 These qeq fixes will ignore electric field contributions from
 :doc:`fix efield <fix_efield>`.
 

doc/src/molecule.rst

@@ -126,14 +126,50 @@ molecule (header keyword = inertia).
 Format of a molecule file
 """""""""""""""""""""""""
 
-The format of an individual molecule file is similar but
+The format of an individual molecule file looks similar but is
-(not identical) to the data file read by the :doc:`read_data <read_data>`
+different than that of a data file read by the :doc:`read_data <read_data>`
-commands, and is as follows.
+commands.  Here is a simple example for a TIP3P water molecule:
+
+.. code-block::
+
+   # Water molecule. TIP3P geometry
+   # header section:
+   3 atoms
+   2 bonds
+   1 angles
+
+   # body section:
+   Coords
+
+   1    0.00000  -0.06556   0.00000
+   2    0.75695   0.52032   0.00000
+   3   -0.75695   0.52032   0.00000
+
+   Types
+
+   1        1   # O
+   2        2   # H
+   3        2   # H
+
+   Charges
+
+   1       -0.834
+   2        0.417
+   3        0.417
+
+   Bonds
+
+   1   1      1      2
+   2   1      1      3
+
+   Angles
+
+   1   1      2      1      3
 
 A molecule file has a header and a body.  The header appears first.  The
-first line of the header and thus of the molecule file is *always* skipped;
+first line of the header and thus of the molecule file is *always*
-it typically contains a description of the file or a comment from the software
+skipped; it typically contains a description of the file or a comment
-that created the file.
+from the software that created the file.
 
 Then lines are read one line at a time.  Lines can have a trailing
 comment starting with '#' that is ignored.  There *must* be at least one
@@ -158,25 +194,62 @@ appear if the value(s) are different than the default, except when
 defining a *body* particle, which requires setting the number of
 *atoms* to 1, and setting the *inertia* in a specific section (see below).
 
-* N *atoms* = # of atoms N in molecule, default = 0
+   .. list-table::
-* Nb *bonds* = # of bonds Nb in molecule, default = 0
+      :header-rows: 1
-* Na *angles* = # of angles Na in molecule, default = 0
+      :widths: auto
-* Nd *dihedrals* = # of dihedrals Nd in molecule, default = 0
-* Ni *impropers* = # of impropers Ni in molecule, default = 0
-* Nf *fragments* = # of fragments Nf in molecule, default = 0
-* Ninteger Ndouble *body* = # of integer and floating-point values
-  in body particle, default = 0
-* Mtotal *mass* = total mass of molecule
-* Xc Yc Zc *com* = coordinates of center-of-mass of molecule
-* Ixx Iyy Izz Ixy Ixz Iyz *inertia* = 6 components of inertia tensor of molecule
 
-For *mass*, *com*, and *inertia*, the default is for LAMMPS to
+      * - Number(s)
-calculate this quantity itself if needed, assuming the molecules
+        - Keyword
-consist of a set of point particles or finite-size particles (with a
+        - Meaning
-non-zero diameter) that do not overlap.  If finite-size particles in
+        - Default Value
-the molecule do overlap, LAMMPS will not account for the overlap
+      * - N
-effects when calculating any of these 3 quantities, so you should
+        - atoms
-pre-compute them yourself and list the values in the file.
+        - # of atoms N in molecule
+        - 0
+      * - Nb
+        - bonds
+        - # of bonds Nb in molecule
+        - 0
+      * - Na
+        - angles
+        - # of angles Na in molecule
+        - 0
+      * - Nd
+        - dihedrals
+        - # of dihedrals Nd in molecule
+        - 0
+      * - Ni
+        - impropers
+        - # of impropers Ni in molecule
+        - 0
+      * - Nf
+        - fragments
+        - # of fragments Nf in molecule
+        - 0
+      * - Ninteger Ndouble
+        - body
+        - # of integer and floating-point values in body particle
+        - 0
+      * - Mtotal
+        - mass
+        - total mass of molecule
+        - computed
+      * - Xc Yc Zc
+        - com
+        - coordinates of center-of-mass of molecule
+        - computed
+      * - Ixx Iyy Izz Ixy Ixz Iyz
+        - inertia
+        - 6 components of inertia tensor of molecule
+        - computed
+
+For *mass*, *com*, and *inertia*, the default is for LAMMPS to calculate
+this quantity itself if needed, assuming the molecules consist of a set
+of point particles or finite-size particles (with a non-zero diameter)
+that do **not** overlap.  If finite-size particles in the molecule
+**do** overlap, LAMMPS will not account for the overlap effects when
+calculating any of these 3 quantities, so you should pre-compute them
+yourself and list the values in the file.
 
 The mass and center-of-mass coordinates (Xc,Yc,Zc) are
 self-explanatory.  The 6 moments of inertia (ixx,iyy,izz,ixy,ixz,iyz)
@@ -188,7 +261,7 @@ internally.
 
 These are the allowed section keywords for the body of the file.
 
-* *Coords, Types, Molecules, Fragments, Charges, Diameters, Masses* = atom-property sections
+* *Coords, Types, Molecules, Fragments, Charges, Diameters, Dipoles, Masses* = atom-property sections
 * *Bonds, Angles, Dihedrals, Impropers* = molecular topology sections
 * *Special Bond Counts, Special Bonds* = special neighbor info
 * *Shake Flags, Shake Atoms, Shake Bond Types* = SHAKE info
@@ -303,6 +376,21 @@ not listed, the default diameter of each atom in the molecule is 1.0.
 
 ----------
 
+.. versionadded:: TBD
+
+*Dipoles* section:
+
+* one line per atom
+* line syntax: ID mux muy muz
+* mux,muy,muz = x-, y-, and z-component of point dipole vector of atom
+
+This section is only allowed for :doc:`atom styles <atom_style>` that
+support particles with point dipoles, e.g. atom_style dipole.  If not
+listed, the default dipole component of each atom in the molecule is set
+to 0.0.
+
+----------
+
 *Masses* section:
 
 * one line per atom

doc/src/neb.rst

@@ -10,7 +10,7 @@ Syntax
 
    neb etol ftol N1 N2 Nevery file-style arg keyword values
 
-* etol = stopping tolerance for energy (energy units)
+* etol = stopping tolerance for energy (dimensionless)
 * ftol = stopping tolerance for force (force units)
 * N1 = max # of iterations (timesteps) to run initial NEB
 * N2 = max # of iterations (timesteps) to run barrier-climbing NEB
@@ -89,10 +89,11 @@ potentials, and the starting configuration when the neb command is
 issued should be the same for every replica.
 
 In a NEB calculation each replica is connected to other replicas by
-inter-replica nudging forces.  These forces are imposed by the :doc:`fix neb <fix_neb>` command, which must be used in conjunction with the
+inter-replica nudging forces.  These forces are imposed by the
-neb command.  The group used to define the fix neb command defines the
+:doc:`fix neb <fix_neb>` command, which must be used in conjunction
-NEB atoms which are the only ones that inter-replica springs are
+with the neb command.  The group used to define the fix neb command
-applied to.  If the group does not include all atoms, then non-NEB
+defines the NEB atoms which are the only ones that inter-replica springs
+are applied to.  If the group does not include all atoms, then non-NEB
 atoms have no inter-replica springs and the forces they feel and their
 motion is computed in the usual way due only to other atoms within
 their replica.  Conceptually, the non-NEB atoms provide a background
@@ -445,7 +446,7 @@ Related commands
 """"""""""""""""
 
 :doc:`prd <prd>`, :doc:`temper <temper>`, :doc:`fix langevin <fix_langevin>`,
-:doc:`fix viscous <fix_viscous>`
+:doc:`fix viscous <fix_viscous>`, :doc:`fix neb <fix_neb>`
 
 Default
 """""""

doc/src/package.rst

@@ -344,12 +344,10 @@ specify additional flags for the runtime build.
 
 ----------
 
-The *intel* style invokes settings associated with the use of the
+The *intel* style invokes settings associated with the use of the INTEL
-INTEL package.  All of its settings, except the *omp* and *mode*
+package.  The keywords *balance*, *ghost*, *tpc*, and *tptask* are
-keywords, are ignored if LAMMPS was not built with Xeon Phi
+**only** applicable if LAMMPS was built with Xeon Phi co-processor
-co-processor support.  All of its settings, including the *omp* and
+support and are otherwise ignored.
-*mode* keyword are applicable if LAMMPS was built with co-processor
-support.
 
 The *Nphi* argument sets the number of co-processors per node.
 This can be set to any value, including 0, if LAMMPS was not
@@ -474,13 +472,13 @@ If the *neigh/thread* keyword is set to *off*, then the KOKKOS package
 threads only over atoms. However, for small systems, this may not expose
 enough parallelism to keep a GPU busy. When this keyword is set to *on*,
 the KOKKOS package threads over both atoms and neighbors of atoms. When
-using *neigh/thread* *on*, a full neighbor list must also be used. Using
+using *neigh/thread* *on*, the :doc:`newton pair <newton>` setting must
-*neigh/thread* *on* may be slower for large systems, so this this option
+be "off". Using *neigh/thread* *on* may be slower for large systems, so
-is turned on by default only when there are 16K atoms or less owned by
+this this option is turned on by default only when running on one or
-an MPI rank and when using a full neighbor list. Not all KOKKOS-enabled
+more GPUs and there are 16k atoms or less owned by an MPI rank. Not all
-potentials support this keyword yet, and only thread over atoms. Many
+KOKKOS-enabled potentials support this keyword yet, and only thread over
-simple pairwise potentials such as Lennard-Jones do support threading
+atoms. Many simple pairwise potentials such as Lennard-Jones do support
-over both atoms and neighbors.
+threading over both atoms and neighbors.
 
 If the *neigh/transpose* keyword is set to *off*, then the KOKKOS
 package will use the same memory layout for building the neighbor list on
@@ -732,7 +730,7 @@ comm = device, sort = device, neigh/transpose = off, gpu/aware = on. When
 LAMMPS can safely detect that GPU-aware MPI is not available, the default value
 of gpu/aware becomes "off". For CPUs or Xeon Phis, the option defaults are
 neigh = half, neigh/qeq = half, newton = on, binsize = 0.0, comm = no, and sort
-= no.  The option neigh/thread = on when there are 16K atoms or less on an MPI
+= no. For GPUs, option neigh/thread = on when there are 16k atoms or less on an MPI
 rank, otherwise it is "off". These settings are made automatically by the
 required "-k on" :doc:`command-line switch <Run_options>`. You can change them
 by using the package kokkos command in your input script or via the :doc:`-pk

doc/src/pair_aip_water_2dm.rst

@@ -22,13 +22,24 @@ Examples
 .. code-block:: LAMMPS
 
    pair_style  hybrid/overlay aip/water/2dm 16.0 1
-   pair_coeff  * * aip/water/2dm  COH.aip.water.2dm C Ow Hw
+   pair_coeff  * * aip/water/2dm  CBNOH.aip.water.2dm C Ow Hw
 
    pair_style  hybrid/overlay aip/water/2dm 16.0 lj/cut/tip4p/long 2 3 1 1 0.1546 10 8.5
    pair_coeff  2 2   lj/cut/tip4p/long    8.0313e-3  3.1589      # O-O
    pair_coeff  2 3   lj/cut/tip4p/long    0.0        0.0         # O-H
    pair_coeff  3 3   lj/cut/tip4p/long    0.0        0.0         # H-H
-   pair_coeff  * *   aip/water/2dm        COH.aip.water.2dm    C Ow Hw
+   pair_coeff  * *   aip/water/2dm        CBNOH.aip.water.2dm    C Ow Hw
+
+   pair_style  hybrid/overlay aip/water/2dm 16.0 lj/cut/tip4p/long 3 4 1 1 0.1546 10 8.5 coul/shield 16.0 1
+   pair_coeff  1*2 1*2   none
+   pair_coeff  3 3   lj/cut/tip4p/long    8.0313e-3  3.1589      # O-O
+   pair_coeff  3 4   lj/cut/tip4p/long    0.0        0.0         # O-H
+   pair_coeff  4 4   lj/cut/tip4p/long    0.0        0.0         # H-H
+   pair_coeff  * *   aip/water/2dm        CBNOH.aip.water.2dm  B N Ow Hw
+   pair_coeff  1 3   coul/shield          1.333
+   pair_coeff  1 4   coul/shield          1.333
+   pair_coeff  2 3   coul/shield          1.333
+   pair_coeff  2 4   coul/shield          1.333
 
 Description
 """""""""""
@@ -37,7 +48,7 @@ Description
 
 The *aip/water/2dm* style computes the anisotropic interfacial potential
 (AIP) potential for interfaces of water with two-dimensional (2D)
-materials as described in :ref:`(Feng) <Feng>`.
+materials as described in :ref:`(Feng1) <Feng1>` and :ref:`(Feng2) <Feng2>`.
 
 .. math::
 
@@ -62,12 +73,12 @@ larger than :math:`r_c` :doc:`pair_style ilp_graphene_hbn
 .. note::
 
    This pair style uses the atomic normal vector definition from
-   :ref:`(Feng) <Feng>`), where the atomic normal vectors of the
+   :ref:`(Feng1) <Feng1>`), where the atomic normal vectors of the
    hydrogen atoms are assumed to lie along the corresponding
    oxygen-hydrogen bonds and the normal vector of the central oxygen
    atom is defined as their average.
 
-The provided parameter file, ``COH.aip.water.2dm``, is intended for use
+The provided parameter file, ``CBNOH.aip.water.2dm``, is intended for use
 with *metal* :doc:`units <units>`, with energies in meV.  Two additional
 parameters, *S*, and *rcut* are included in the parameter file. *S* is
 designed to facilitate scaling of energies; *rcut* is the cutoff for an
@@ -77,7 +88,7 @@ the calculation of the normals for all atom pairs.
 .. note::
 
    The parameters presented in the provided parameter file,
-   ``COH.aip.water.2dm``, are fitted with the taper function enabled by
+   ``CBNOH.aip.water.2dm``, are fitted with the taper function enabled by
    setting the cutoff equal to 16.0 Angstrom.  Using a different cutoff
    or taper function setting should be carefully checked as they can
    lead to significant errors.  These parameters provide a good
@@ -134,7 +145,7 @@ if LAMMPS was built with that package.  See the :doc:`Build package
 This pair style requires the newton setting to be *on* for pair
 interactions.
 
-The ``COH.aip.water.2dm`` potential file provided with LAMMPS is
+The ``CBNOH.aip.water.2dm`` potential file provided with LAMMPS is
 parameterized for *metal* units.  You can use this pair style with any
 LAMMPS units, but you would need to create your own potential file with
 parameters in the appropriate units, if your simulation does not use
@@ -162,6 +173,10 @@ tap_flag = 1
 
 ----------
 
-.. _Feng:
+.. _Feng1:
 
-**(Feng)** Z. Feng and W. Ouyang et al., J. Phys. Chem. C. 127, 8704-8713 (2023).
+**(Feng1)** Z. Feng, ..., and W. Ouyang, J. Phys. Chem. C. 127(18), 8704-8713 (2023).
+
+.. _Feng2:
+
+**(Feng2)** Z. Feng, ..., and W. Ouyang, Langmuir 39(50), 18198-18207 (2023).

doc/src/pair_charmm.rst

@@ -16,6 +16,7 @@
 .. index:: pair_style lj/charmm/coul/msm/omp
 .. index:: pair_style lj/charmmfsw/coul/charmmfsh
 .. index:: pair_style lj/charmmfsw/coul/long
+.. index:: pair_style lj/charmmfsw/coul/long/kk
 
 pair_style lj/charmm/coul/charmm command
 ========================================
@@ -43,6 +44,8 @@ pair_style lj/charmmfsw/coul/charmmfsh command
 pair_style lj/charmmfsw/coul/long command
 =========================================
 
+Accelerator Variants: *lj/charmmfsw/coul/long/kk*
+
 Syntax
 """"""
 
@@ -281,7 +284,9 @@ page for more info.
 Related commands
 """"""""""""""""
 
-:doc:`pair_coeff <pair_coeff>`
+:doc:`pair_coeff <pair_coeff>`, :doc:`angle_style charmm <angle_charmm>`,
+:doc:`dihedral_style charmm <dihedral_charmm>`,
+:doc:`dihedral_style charmmfsw <dihedral_charmm>`, :doc:`fix cmap <fix_cmap>`
 
 Default
 """""""

doc/src/pair_ilp_tmd.rst

@@ -22,12 +22,12 @@ Examples
 .. code-block:: LAMMPS
 
    pair_style  hybrid/overlay ilp/tmd 16.0 1
-   pair_coeff  * * ilp/tmd  MoS2.ILP Mo S S
+   pair_coeff  * * ilp/tmd  TMD.ILP Mo S S
 
    pair_style  hybrid/overlay sw/mod sw/mod ilp/tmd 16.0
    pair_coeff  * * sw/mod 1  tmd.sw.mod Mo S S NULL NULL NULL
-   pair_coeff  * * sw/mod 2  tmd.sw.mod NULL NULL NULL Mo S S
+   pair_coeff  * * sw/mod 2  tmd.sw.mod NULL NULL NULL W Se Se
-   pair_coeff  * * ilp/tmd   MoS2.ILP   Mo S S Mo S S
+   pair_coeff  * * ilp/tmd   TMD.ILP   Mo S S W Se Se
 
 Description
 """""""""""
@@ -36,7 +36,7 @@ Description
 
 The *ilp/tmd* style computes the registry-dependent interlayer
 potential (ILP) potential for transition metal dichalcogenides (TMD)
-as described in :ref:`(Ouyang7) <Ouyang7>`.
+as described in :ref:`(Ouyang7) <Ouyang7>` and :ref:`(Jiang) <Jiang>`.
 
 .. math::
 
@@ -69,7 +69,7 @@ calculating the normals.
    each atom `i`, its six nearest neighboring atoms belonging to the same
    sub-layer are chosen to define the normal vector `{\bf n}_i`.
 
-The parameter file (e.g. MoS2.ILP), is intended for use with *metal*
+The parameter file (e.g. TMD.ILP), is intended for use with *metal*
 :doc:`units <units>`, with energies in meV. Two additional parameters,
 *S*, and *rcut* are included in the parameter file. *S* is designed to
 facilitate scaling of energies. *rcut* is designed to build the neighbor
@@ -77,7 +77,7 @@ list for calculating the normals for each atom pair.
 
 .. note::
 
-   The parameters presented in the parameter file (e.g. MoS2.ILP),
+   The parameters presented in the parameter file (e.g. TMD.ILP),
    are fitted with taper function by setting the cutoff equal to 16.0
    Angstrom.  Using different cutoff or taper function should be careful.
    These parameters provide a good description in both short- and long-range
@@ -133,10 +133,10 @@ if LAMMPS was built with that package.  See the :doc:`Build package
 This pair style requires the newton setting to be *on* for pair
 interactions.
 
-The MoS2.ILP potential file provided with LAMMPS (see the potentials
+The TMD.ILP potential file provided with LAMMPS (see the potentials
 directory) are parameterized for *metal* units.  You can use this
 potential with any LAMMPS units, but you would need to create your own
-custom MoS2.ILP potential file with coefficients listed in the appropriate
+custom TMD.ILP potential file with coefficients listed in the appropriate
 units, if your simulation does not use *metal* units.
 
 Related commands
@@ -164,3 +164,7 @@ tap_flag = 1
 .. _Ouyang7:
 
 **(Ouyang7)** W. Ouyang, et al., J. Chem. Theory Comput. 17, 7237 (2021).
+
+.. _Jiang:
+
+**(Jiang)** W. Jiang, et al., J. Phys. Chem. A, 127, 46, 9820-9830 (2023).

doc/src/pair_lepton.rst

@@ -72,7 +72,7 @@ interactions between particles which depend on the distance and have a
 cutoff.  The potential function must be provided as an expression string
 using "r" as the distance variable.  With pair style *lepton/coul* one
 may additionally reference the charges of the two atoms of the pair with
-"qi" and "qj", respectively.  With pair style *lepton/coul* one may
+"qi" and "qj", respectively.  With pair style *lepton/sphere* one may
 instead reference the radii of the two atoms of the pair with "radi" and
 "radj", respectively; this is half of the diameter that can be set in
 :doc:`data files <read_data>` or the :doc:`set command <set>`.
@@ -166,8 +166,8 @@ mixing.  Thus, expressions for *all* I,J pairs must be specified
 explicitly.
 
 Only pair style *lepton* supports the :doc:`pair_modify shift <pair_modify>`
-option for shifting the energy of the pair interaction so that it is
+option for shifting the potential energy of the pair interaction so that
-0 at the cutoff, pair styles *lepton/coul* and *lepton/sphere* do *not*.
+it is 0 at the cutoff, pair styles *lepton/coul* and *lepton/sphere* do *not*.
 
 The :doc:`pair_modify table <pair_modify>` options are not relevant for
 the these pair styles.

doc/src/run_style.rst

@@ -329,7 +329,8 @@ Restrictions
 The *verlet/split* style can only be used if LAMMPS was built with the
 REPLICA package. Correspondingly the *respa/omp* style is available
 only if the OPENMP package was included. See the :doc:`Build package
-<Build_package>` page for more info.
+<Build_package>` page for more info.  It is not compatible with
+kspace styles from the INTEL package.
 
 Whenever using rRESPA, the user should experiment with trade-offs in
 speed and accuracy for their system, and verify that they are

doc/utils/sphinx-config/false_positives.txt

@@ -125,6 +125,7 @@ antisymmetry
 anton
 Antonelli
 api
+apolar
 Apoorva
 Appl
 Appshaw
@@ -151,6 +152,7 @@ asphericity
 Asq
 assignee
 assively
+associativity
 Asta
 Astart
 Astop
@@ -792,11 +794,13 @@ dispersionflag
 dissipative
 Dissipative
 distharm
+distutils
 dl
 dlabel
 dlambda
 DLAMMPS
 dll
+dlm
 dlopen
 dm
 dmax
@@ -1014,6 +1018,7 @@ Ercolessi
 Erdmann
 erf
 erfc
+erforce
 Erhart
 erorate
 erose
@@ -2234,8 +2239,10 @@ Mohd
 Mohles
 mol
 Mol
+molatom
 molfile
 Molfile
+molindex
 MolPairStyle
 moltemplate
 momb
@@ -2568,6 +2575,7 @@ ns
 Ns
 Nsample
 Nskip
+nspecial
 Nspecies
 nsq
 Nstart
@@ -3876,6 +3884,7 @@ versa
 Verstraelen
 ves
 vflag
+vfrac
 vhi
 vibrational
 Vij

examples/COUPLE/plugin/liblammpsplugin.c

@@ -110,6 +110,8 @@ liblammpsplugin_t *liblammpsplugin_load(const char *lib)
   ADDSYM(extract_variable);
   ADDSYM(extract_variable_datatype);
   ADDSYM(set_variable);
+  ADDSYM(set_string_variable);
+  ADDSYM(set_internal_variable);
   ADDSYM(variable_info);
 
   ADDSYM(gather_atoms);

examples/COUPLE/plugin/liblammpsplugin.h

@@ -152,9 +152,11 @@ struct _liblammpsplugin {
 
   void *(*extract_compute)(void *, const char *, int, int);
   void *(*extract_fix)(void *, const char *, int, int, int, int);
-  void *(*extract_variable)(void *, const char *, char *);
+  void *(*extract_variable)(void *, const char *, const char *);
   int (*extract_variable_datatype)(void *, const char *);
-  int (*set_variable)(void *, char *, char *);
+  int (*set_variable)(void *, const char *, const char *);
+  int (*set_string_variable)(void *, const char *, const char *);
+  int (*set_internal_variable)(void *, const char *, double);
   int (*variable_info)(void *, int, char *, int);
 
   void (*gather_atoms)(void *, const char *, int, int, void *);

examples/COUPLE/simple/simple.cpp

@@ -67,7 +67,7 @@ int main(int narg, char **arg)
   FILE *fp;
   if (me == 0) {
     fp = fopen(arg[2],"r");
-    if (fp == NULL) {
+    if (fp == nullptr) {
       printf("ERROR: Could not open LAMMPS input script\n");
       MPI_Abort(MPI_COMM_WORLD,1);
     }
@@ -78,14 +78,14 @@ int main(int narg, char **arg)
   // (could just send it to proc 0 of comm_lammps and let it Bcast)
   // all LAMMPS procs call input->one() on the line
   
-  LAMMPS *lmp = NULL;
+  LAMMPS *lmp = nullptr;
-  if (lammps == 1) lmp = new LAMMPS(0,NULL,comm_lammps);
+  if (lammps == 1) lmp = new LAMMPS(0,nullptr,comm_lammps);
 
   int n;
   char line[1024];
-  while (1) {
+  while (true) {
     if (me == 0) {
-      if (fgets(line,1024,fp) == NULL) n = 0;
+      if (fgets(line,1024,fp) == nullptr) n = 0;
       else n = strlen(line) + 1;
       if (n == 0) fclose(fp);
     }
@@ -101,8 +101,8 @@ int main(int narg, char **arg)
   // put coords back into LAMMPS
   // run a single step with changed coords
 
-  double *x = NULL;
+  double *x = nullptr;
-  double *v = NULL;
+  double *v = nullptr;
 
   if (lammps == 1) {
     lmp->input->one("run 10");
@@ -147,7 +147,7 @@ int main(int narg, char **arg)
   // create_atoms() to create new ones with old coords, vels
   // initial thermo should be same as step 20
 
-  int *type = NULL;
+  int *type = nullptr;
 
   if (lammps == 1) {
     int natoms = static_cast<int> (lmp->atom->natoms);
@@ -155,7 +155,7 @@ int main(int narg, char **arg)
     for (int i = 0; i < natoms; i++) type[i] = 1;
 
     lmp->input->one("delete_atoms group all");
-    lammps_create_atoms(lmp,natoms,NULL,type,x,v,NULL,0);
+    lammps_create_atoms(lmp,natoms,nullptr,type,x,v,nullptr,0);
     lmp->input->one("run 10");
   }
 

examples/PACKAGES/cgdna/util/generate.py

@@ -22,6 +22,10 @@
 """
 Import basic modules
 """
+
+# for python2/3 compatibility
+from __future__ import print_function
+
 import sys, os, timeit
 
 from timeit import default_timer as timer
@@ -33,7 +37,7 @@ which needs to be provided
 try:
     import numpy as np
 except:
-    print >> sys.stderr, "numpy not found. Exiting."
+    print("numpy not found. Exiting.", file=sys.stderr)
     sys.exit(1)
 
 """
@@ -45,8 +49,8 @@ try:
     box_length = float(sys.argv[2])
     infile = sys.argv[3]
 except:
-    print >> sys.stderr, "Usage: %s <%s> <%s> <%s>" % (sys.argv[0], \
+    print( "Usage: %s <%s> <%s> <%s>" % (sys.argv[0], \
-	"box offset", "box length", "file with sequences")
+	"box offset", "box length", "file with sequences"), file=sys.stderr)
     sys.exit(1)
 box = np.array ([box_length, box_length, box_length])
 
@@ -57,8 +61,7 @@ try:
     inp = open (infile, 'r')
     inp.close()
 except:
-    print >> sys.stderr, "Could not open file '%s' for reading. \
+    print( "Could not open file '%s' for reading. Aborting." % infile, file=sys.stderr)
-					      Aborting." % infile
     sys.exit(2)
 
 # return parts of a string
@@ -174,21 +177,21 @@ def add_strands (mynewpositions, mynewa1s, mynewa3s):
     # placed particles i we check whether it overlaps with any of the
     # newly created particles j
 
-    print >> sys.stdout, "## Checking for overlaps"
+    print( "## Checking for overlaps", file=sys.stdout)
 
-    for i in xrange(len(positions)):
+    for i in range(len(positions)):
 
         p = positions[i]
         pa1 = a1s[i]
 
-	for j in xrange (len(mynewpositions)):
+        for j in range (len(mynewpositions)):
 
             q = mynewpositions[j]
             qa1 = mynewa1s[j]
 
             # skip particles that are anyway too far away
             dr = p - q
-	    dr -= box * np.rint (dr / box)
+            dr -= box * np.rint(dr / box)
             if np.dot(dr, dr) > RC2:
                 continue
 
@@ -200,13 +203,13 @@ def add_strands (mynewpositions, mynewa1s, mynewa3s):
 
             # check for no overlap between the two backbone sites
             dr = p_pos_back - q_pos_back
-            dr -= box * np.rint (dr / box)
+            dr -= box * np.rint(dr / box)
             if np.dot(dr, dr) < RC2_BACK:
                 overlap = True
 
             # check for no overlap between the two base sites
             dr = p_pos_base -  q_pos_base
-            dr -= box * np.rint (dr / box)
+            dr -= box * np.rint(dr / box)
             if np.dot(dr, dr) < RC2_BASE:
                 overlap = True
 
@@ -238,7 +241,7 @@ def add_strands (mynewpositions, mynewa1s, mynewa3s):
         for p in mynewa3s:
             a3s.append (p)
         # calculate quaternion from local body frame and append
-	for ia in xrange(len(mynewpositions)):
+        for ia in range(len(mynewpositions)):
             mynewquaternions = exyz_to_quat(mynewa1s[ia],mynewa3s[ia])
             quaternions.append(mynewquaternions)
 
@@ -301,13 +304,12 @@ def generate_strand(bp, sequence=None, start_pos=np.array([0, 0, 0]), \
     elif len(sequence) != bp:
         n = bp - len(sequence)
         sequence += np.random.randint(1, 5, n)
-	print >> sys.stderr, "sequence is too short, adding %d random bases" % n
+        print( "sequence is too short, adding %d random bases" % n, file=sys.stderr)
 
     # normalize direction
     dir_norm = np.sqrt(np.dot(dir,dir))
     if dir_norm < 1e-10:
-	print >> sys.stderr, "direction must be a valid vector, \
+        print( "direction must be a valid vector, defaulting to (0, 0, 1)", file=sys.stderr)
-			      defaulting to (0, 0, 1)"
         dir = np.array([0, 0, 1])
     else: dir /= dir_norm
 
@@ -391,7 +393,7 @@ def read_strands(filename):
     try:
         infile = open (filename)
     except:
-        print >> sys.stderr, "Could not open file '%s'. Aborting." % filename
+        print( "Could not open file '%s'. Aborting." % filename, file=sys.stderr )
         sys.exit(2)
 
     # This block works out the number of nucleotides and strands by reading
@@ -406,28 +408,27 @@ def read_strands(filename):
         if line[:6] == 'DOUBLE':
             line = line.split()[1]
             length = len(line)
-            print >> sys.stdout, "## Found duplex of %i base pairs" % length
+            print( "## Found duplex of %i base pairs" % length, file=sys.stdout)
             nnucl += 2*length
             nstrands += 2
             nbonds += (2*length-2)
         else:
             line = line.split()[0]
             length = len(line)
-            print >> sys.stdout, \
+            print( "## Found single strand of %i bases" % length, file=sys.stdout)
-		    "## Found single strand of %i bases" % length
             nnucl += length
             nstrands += 1
             nbonds += length-1
     # rewind the sequence input file
     infile.seek(0)
 
-    print >> sys.stdout, "## nstrands, nnucl = ", nstrands, nnucl
+    print( "## nstrands, nnucl = ", nstrands, nnucl, file=sys.stdout)
 
     # generate the data file in LAMMPS format
     try:
         out = open ("data.oxdna", "w")
     except:
-        print >> sys.stderr, "Could not open data file for writing. Aborting."
+        print( "Could not open data file for writing. Aborting.", file=sys.stderr)
         sys.exit(2)
 
     lines = infile.readlines()
@@ -452,11 +453,11 @@ def read_strands(filename):
             seq = [(base_to_number[x]) for x in line]
 
             myns += 1
-	    for b in xrange(length):
+            for b in range(length):
                 basetype.append(seq[b])
                 strandnum.append(myns)
 
-	    for b in xrange(length-1):
+            for b in range(length-1):
                 bondpair = [noffset + b, noffset + b + 1]
                 bonds.append(bondpair)
             noffset += length
@@ -467,16 +468,16 @@ def read_strands(filename):
             seq2.reverse()
 
             myns += 1
-	    for b in xrange(length):
+            for b in range(length):
                 basetype.append(seq2[b])
                 strandnum.append(myns)
 
-	    for b in xrange(length-1):
+            for b in range(length-1):
                 bondpair = [noffset + b, noffset + b + 1]
                 bonds.append(bondpair)
             noffset += length
 
-            print >> sys.stdout, "## Created duplex of %i bases" % (2*length)
+            print( "## Created duplex of %i bases" % (2*length), file=sys.stdout)
 
             # generate random position of the first nucleotide
             cdm = box_offset + np.random.random_sample(3) * box
@@ -499,10 +500,10 @@ def read_strands(filename):
                 axis /= np.sqrt(np.dot(axis, axis))
                 newpositions, newa1s, newa3s = generate_strand(len(line), \
                     sequence=seq, dir=axis, start_pos=cdm, double=True)
-                print >> sys.stdout, "## Trying %i" % i
+                print( "## Trying %i" % i, file=sys.stdout)
             end = timer()
-            print >> sys.stdout, "## Added duplex of %i bases (line %i/%i) in %.2fs, now at %i/%i" % \
+            print( "## Added duplex of %i bases (line %i/%i) in %.2fs, now at %i/%i" % \
-				      (2*length, i, nlines, end-start, len(positions), nnucl)
+				      (2*length, i, nlines, end-start, len(positions), nnucl), file=sys.stdout)
 
         # block for single strands: last argument of the generate function
         # is set to 'False'
@@ -511,11 +512,11 @@ def read_strands(filename):
             seq = [(base_to_number[x]) for x in line]
 
             myns += 1
-	    for b in xrange(length):
+            for b in range(length):
                 basetype.append(seq[b])
                 strandnum.append(myns)
 
-	    for b in xrange(length-1):
+            for b in range(length-1):
                 bondpair = [noffset + b, noffset + b + 1]
                 bonds.append(bondpair)
             noffset += length
@@ -527,8 +528,7 @@ def read_strands(filename):
             axis = np.random.random_sample(3)
             axis /= np.sqrt(np.dot(axis, axis))
 
-            print >> sys.stdout, \
+            print("## Created single strand of %i bases" % length, file=sys.stdout)
-		      "## Created single strand of %i bases" % length
 
             newpositions, newa1s, newa3s = generate_strand(length, \
 		      sequence=seq, dir=axis, start_pos=cdm, double=False)
@@ -541,14 +541,14 @@ def read_strands(filename):
                     sequence=seq, dir=axis, start_pos=cdm, double=False)
                 print >> sys.stdout, "## Trying  %i" % (i)
             end = timer()
-            print >> sys.stdout, "## Added single strand of %i bases (line %i/%i) in %.2fs, now at %i/%i" % \
+            print( "## Added single strand of %i bases (line %i/%i) in %.2fs, now at %i/%i" % \
-				      (length, i, nlines, end-start,len(positions), nnucl)
+				      (length, i, nlines, end-start,len(positions), nnucl), file=sys.stdout)
 
         i += 1
 
     # sanity check
     if not len(positions) == nnucl:
-        print len(positions), nnucl
+        print( len(positions), nnucl )
         raise AssertionError
 
     out.write('# LAMMPS data file\n')
@@ -580,18 +580,16 @@ def read_strands(filename):
     out.write('Atoms\n')
     out.write('\n')
 
-    for i in xrange(nnucl):
+    for i in range(nnucl):
         out.write('%d %d %22.15le %22.15le %22.15le %d 1 1\n' \
-		  % (i+1, basetype[i], \
+            % (i+1, basetype[i], positions[i][0], positions[i][1], positions[i][2], strandnum[i]))
-		     positions[i][0], positions[i][1], positions[i][2], \
-		     strandnum[i]))
 
     out.write('\n')
     out.write('# Atom-ID, translational, rotational velocity\n')
     out.write('Velocities\n')
     out.write('\n')
 
-    for i in xrange(nnucl):
+    for i in range(nnucl):
         out.write("%d %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le\n" \
             % (i+1,0.0,0.0,0.0,0.0,0.0,0.0))
 
@@ -600,9 +598,8 @@ def read_strands(filename):
     out.write('Ellipsoids\n')
     out.write('\n')
 
-    for i in xrange(nnucl):
+    for i in range(nnucl):
-	out.write(\
+        out.write("%d %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le\n"  \
-    "%d %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le\n"  \
             % (i+1,1.1739845031423408,1.1739845031423408,1.1739845031423408, \
             quaternions[i][0],quaternions[i][1], quaternions[i][2],quaternions[i][3]))
 
@@ -611,13 +608,13 @@ def read_strands(filename):
     out.write('Bonds\n')
     out.write('\n')
 
-    for i in xrange(nbonds):
+    for i in range(nbonds):
         out.write("%d  %d  %d  %d\n" % (i+1,1,bonds[i][0],bonds[i][1]))
 
     out.close()
 
-    print >> sys.stdout, "## Wrote data to 'data.oxdna'"
+    print("## Wrote data to 'data.oxdna'", file=sys.stdout)
-    print >> sys.stdout, "## DONE"
+    print("## DONE", file=sys.stdout)
 
 # call the above main() function, which executes the program
 read_strands (infile)
@@ -627,4 +624,6 @@ runtime = end_time-start_time
 hours = runtime/3600
 minutes = (runtime-np.rint(hours)*3600)/60
 seconds = (runtime-np.rint(hours)*3600-np.rint(minutes)*60)%60
-print >> sys.stdout, "## Total runtime %ih:%im:%.2fs" % (hours,minutes,seconds)
+print( "## Total runtime %ih:%im:%.2fs" % (hours,minutes,seconds), file=sys.stdout)
+
+

examples/PACKAGES/cgdna/util/generate_simple.py

@@ -1,5 +1,8 @@
 # Setup tool for oxDNA input in LAMMPS format.
 
+# for python2/3 compatibility
+from __future__ import print_function
+
 import math,numpy as np,sys,os
 
 # system size

examples/PACKAGES/interlayer/aip_water_2dm/CBNOH.aip.water.2dm

@@ -0,0 +1 @@
+../../../../potentials/CBNOH.aip.water.2dm

examples/PACKAGES/interlayer/aip_water_2dm/COH.aip.water.2dm

@@ -1 +0,0 @@
-../../../../potentials/COH.aip.water.2dm

examples/PACKAGES/interlayer/ilp_tmds/MoS2.ILP

@@ -1 +0,0 @@
-../../../../potentials/MoS2.ILP

examples/PACKAGES/interlayer/ilp_tmds/TMD.ILP

@@ -0,0 +1 @@
+../../../../potentials/TMD.ILP

examples/PACKAGES/interlayer/ilp_tmds/in.mos2

@@ -12,7 +12,7 @@ mass            4 95.94
 pair_style     	hybrid/overlay sw/mod sw/mod ilp/tmd 16.0
 pair_coeff     	* * sw/mod 1 tmd.sw.mod Mo S S NULL NULL NULL
 pair_coeff     	* * sw/mod 2 tmd.sw.mod NULL NULL NULL Mo S S
-pair_coeff     	* * ilp/tmd  MoS2.ILP Mo S S Mo S S
+pair_coeff     	* * ilp/tmd  TMD.ILP Mo S S Mo S S
 
 # Calculate the pair potential
 compute   	0 all pair ilp/tmd

examples/PACKAGES/reaction/create_atoms_polystyrene/in.grow_styrene

@@ -40,7 +40,7 @@ fix 1 statted_grp_REACT nvt temp $T $T 100
 
 fix 4 bond_react_MASTER_group temp/rescale 1 $T $T 1 1
 
-thermo_style custom step temp press density f_myrxns[1]
+thermo_style custom step temp press density f_myrxns[*]
 
 thermo 100
 

examples/PACKAGES/reaction/nylon,6-6_melt/in.large_nylon_melt

@@ -50,7 +50,7 @@ fix 1 statted_grp_REACT nvt temp 800 800 100
 # you can use the internally created 'bond_react_MASTER_group', like so:
 # fix 2 bond_react_MASTER_group temp/rescale 1 800 800 10 1
 
-thermo_style custom step temp press density f_myrxns[1] f_myrxns[2] # cumulative reaction counts
+thermo_style custom step temp press density f_myrxns[*] # cumulative reaction counts
 
 # restart 100 restart1 restart2
 

examples/PACKAGES/reaction/tiny_epoxy/in.tiny_epoxy.stabilized

@@ -20,7 +20,8 @@ improper_style class2
 special_bonds lj/coul 0 0 1
 pair_modify tail yes mix sixthpower
 
-read_data tiny_epoxy.data
+read_data tiny_epoxy.data &
+  extra/special/per/atom 25
 
 velocity all create 300.0 4928459 dist gaussian
 
@@ -44,7 +45,7 @@ fix rxns all bond/react stabilization yes statted_grp .03 &
 
 fix 1 statted_grp_REACT nvt temp 300 300 100
 
-thermo_style custom step temp f_rxns[1] f_rxns[2] f_rxns[3] f_rxns[4]
+thermo_style custom step temp f_rxns[*]
 
 run 2000
 

examples/PACKAGES/reaction/tiny_nylon/in.tiny_nylon.stabilized

@@ -50,7 +50,7 @@ fix 1 statted_grp_REACT nvt temp 300 300 100
 # by using the internally-created 'bond_react_MASTER_group', like so:
 fix 4 bond_react_MASTER_group temp/rescale 1 300 300 10 1
 
-thermo_style custom step temp press density f_myrxns[1] f_myrxns[2]
+thermo_style custom step temp press density f_myrxns[*]
 
 # restart 100 restart1 restart2
 

examples/PACKAGES/reaction/tiny_nylon/in.tiny_nylon.stabilized_variable_probability

@@ -54,7 +54,7 @@ fix 1 statted_grp_REACT nvt temp 300 300 100
 # by using the internally-created 'bond_react_MASTER_group', like so:
 fix 4 bond_react_MASTER_group temp/rescale 1 300 300 10 1
 
-thermo_style custom step temp press density v_prob1 v_prob2 f_myrxns[1] f_myrxns[2]
+thermo_style custom step temp press density v_prob1 v_prob2 f_myrxns[*]
 
 # restart 100 restart1 restart2
 

examples/PACKAGES/reaction/tiny_nylon/in.tiny_nylon.unstabilized

@@ -47,7 +47,7 @@ fix myrxns all bond/react stabilization no &
 
 fix 1 all nve/limit .03
 
-thermo_style custom step temp press density f_myrxns[1] f_myrxns[2]
+thermo_style custom step temp press density f_myrxns[*]
 
 # restart 100 restart1 restart2
 

examples/PACKAGES/reaction/tiny_polystyrene/in.tiny_polystyrene.stabilized

@@ -51,7 +51,7 @@ fix 1 statted_grp_REACT nvt temp $T $T 100
 
 fix 4 bond_react_MASTER_group temp/rescale 1 $T $T 1 1
 
-thermo_style custom step temp press density f_rxn1[1] f_rxn1[2] f_rxn1[3]
+thermo_style custom step temp press density f_rxn1[*]
 
 run 10000
 

fortran/lammps.f90

@@ -118,6 +118,8 @@ MODULE LIBLAMMPS
     PROCEDURE :: extract_fix            => lmp_extract_fix
     PROCEDURE :: extract_variable       => lmp_extract_variable
     PROCEDURE :: set_variable           => lmp_set_variable
+    PROCEDURE :: set_string_variable    => lmp_set_string_variable
+    PROCEDURE :: set_internal_variable  => lmp_set_internal_variable
     PROCEDURE, PRIVATE :: lmp_gather_atoms_int
     PROCEDURE, PRIVATE :: lmp_gather_atoms_double
     GENERIC   :: gather_atoms           => lmp_gather_atoms_int, &
@@ -557,6 +559,21 @@ MODULE LIBLAMMPS
       INTEGER(c_int) :: lammps_set_variable
     END FUNCTION lammps_set_variable
 
+    FUNCTION lammps_set_string_variable(handle, name, str) BIND(C)
+      IMPORT :: c_int, c_ptr
+      IMPLICIT NONE
+      TYPE(c_ptr), VALUE :: handle, name, str
+      INTEGER(c_int) :: lammps_set_string_variable
+    END FUNCTION lammps_set_string_variable
+
+    FUNCTION lammps_set_internal_variable(handle, name, val) BIND(C)
+      IMPORT :: c_int, c_ptr, c_double
+      IMPLICIT NONE
+      TYPE(c_ptr), VALUE :: handle, name
+      REAL(c_double), VALUE :: val
+      INTEGER(c_int) :: lammps_set_internal_variable
+    END FUNCTION lammps_set_internal_variable
+
     SUBROUTINE lammps_gather_atoms(handle, name, type, count, data) BIND(C)
       IMPORT :: c_int, c_ptr
       IMPLICIT NONE
@@ -1631,6 +1648,43 @@ CONTAINS
     END IF
   END SUBROUTINE lmp_set_variable
 
+  ! equivalent function to lammps_set_string_variable
+  SUBROUTINE lmp_set_string_variable(self, name, str)
+    CLASS(lammps), INTENT(IN) :: self
+    CHARACTER(LEN=*), INTENT(IN) :: name, str
+    INTEGER :: err
+    TYPE(c_ptr) :: Cstr, Cname
+
+    Cstr = f2c_string(str)
+    Cname = f2c_string(name)
+    err = lammps_set_string_variable(self%handle, Cname, Cstr)
+    CALL lammps_free(Cname)
+    CALL lammps_free(Cstr)
+    IF (err /= 0) THEN
+      CALL lmp_error(self, LMP_ERROR_WARNING + LMP_ERROR_WORLD, &
+        'WARNING: unable to set string variable "' // name &
+        // '" [Fortran/set_variable]')
+    END IF
+  END SUBROUTINE lmp_set_string_variable
+
+  ! equivalent function to lammps_set_internal_variable
+  SUBROUTINE lmp_set_internal_variable(self, name, val)
+    CLASS(lammps), INTENT(IN) :: self
+    CHARACTER(LEN=*), INTENT(IN) :: name
+    REAL(KIND=c_double), INTENT(IN) :: val
+    INTEGER :: err
+    TYPE(c_ptr) :: Cstr, Cname
+
+    Cname = f2c_string(name)
+    err = lammps_set_internal_variable(self%handle, Cname, val)
+    CALL lammps_free(Cname)
+    IF (err /= 0) THEN
+      CALL lmp_error(self, LMP_ERROR_WARNING + LMP_ERROR_WORLD, &
+        'WARNING: unable to set internal variable "' // name &
+        // '" [Fortran/set_variable]')
+    END IF
+  END SUBROUTINE lmp_set_internal_variable
+
   ! equivalent function to lammps_gather_atoms (for integers)
   SUBROUTINE lmp_gather_atoms_int(self, name, count, data)
     CLASS(lammps), INTENT(IN) :: self

lib/gpu/geryon/nvd_device.h

@@ -138,7 +138,7 @@ class UCL_Device {
   /** \note You cannot delete the default stream **/
   inline void pop_command_queue() {
     if (_cq.size()<2) return;
-    CU_SAFE_CALL_NS(cuStreamDestroy(_cq.back()));
+    cuStreamDestroy(_cq.back());
     _cq.pop_back();
   }
 
@@ -426,8 +426,8 @@ void UCL_Device::clear() {
   if (_device>-1) {
     for (int i=1; i<num_queues(); i++) pop_command_queue();
 #if GERYON_NVD_PRIMARY_CONTEXT
-    CU_SAFE_CALL_NS(cuCtxSetCurrent(_old_context));
+    cuCtxSetCurrent(_old_context);
-    CU_SAFE_CALL_NS(cuDevicePrimaryCtxRelease(_cu_device));
+    cuDevicePrimaryCtxRelease(_cu_device);
 #else
     cuCtxDestroy(_context);
 #endif

lib/gpu/lal_amoeba.cu

@@ -2033,13 +2033,13 @@ __kernel void k_amoeba_special15(__global int * dev_nbor,
                           const __global tagint *restrict special15,
                           const int inum, const int nall, const int nbor_pitch,
                           const int t_per_atom) {
-  int tid, ii, offset, n_stride, i;
+  int tid, ii, offset, n_stride, j;
   atom_info(t_per_atom,ii,tid,offset);
 
   if (ii<inum) {
 
     int numj, nbor, nbor_end;
-    nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
+    nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,j,numj,
               n_stride,nbor_end,nbor);
 
     int n15 = nspecial15[ii];
@@ -2048,7 +2048,7 @@ __kernel void k_amoeba_special15(__global int * dev_nbor,
 
       int sj=dev_packed[nbor];
       int which = sj >> SBBITS & 3;
-      int j = sj & NEIGHMASK;
+      j = sj & NEIGHMASK;
       tagint jtag = tag[j];
 
       if (!which) {

lib/gpu/lal_base_sph.cpp

@@ -195,7 +195,7 @@ void BaseSPHT::compute(const int f_ago, const int inum_full, const int nall,
                        int **firstneigh, const bool eflag_in, const bool vflag_in,
                        const bool eatom, const bool vatom, int &host_start,
                        const double cpu_time, bool &success, tagint *tag,
-                       double **host_v, const int nlocal) {
+                       double **host_v) {
   acc_timers();
   int eflag, vflag;
   if (eatom) eflag=2;

lib/gpu/lal_base_sph.h

@@ -15,7 +15,7 @@
  ***************************************************************************/
 
 #ifndef LAL_BASE_SPH_H
-#define LAL_BASE_DPD_H
+#define LAL_BASE_SPH_H
 
 #include "lal_device.h"
 #include "lal_balance.h"
@@ -132,7 +132,7 @@ class BaseSPH {
                int **firstneigh, const bool eflag, const bool vflag,
                const bool eatom, const bool vatom, int &host_start,
                const double cpu_time, bool &success, tagint *tag,
-               double **v, const int nlocal);
+               double **v);
 
   /// Pair loop with device neighboring
   int** compute(const int ago, const int inum_full, const int nall,

lib/gpu/lal_coul_slater_long.cu

@@ -102,6 +102,7 @@ __kernel void k_coul_slater_long(const __global numtyp4 *restrict x_,
         numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*grij);
         _erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
         fetch(prefactor,j,q_tex);
+        prefactor *= qqrd2e * scale[mtype] * qtmp/r;
         numtyp rlamdainv = r * lamdainv;
         numtyp exprlmdainv = ucl_exp((numtyp)-2.0*rlamdainv);
         numtyp slater_term = exprlmdainv*((numtyp)1.0 + ((numtyp)2.0*rlamdainv*((numtyp)1.0+rlamdainv)));

lib/gpu/lal_coul_slater_long.h

@@ -13,8 +13,8 @@
     email                : ndactrung@gmail.com
  ***************************************************************************/
 
-#ifndef LAL_Coul_Slater_Long_H
+#ifndef LAL_COUL_SLATER_LONG_H
-#define LAL_Coul_Slater_Long_H
+#define LAL_COUL_SLATER_LONG_H
 
 #include "lal_base_charge.h"
 

lib/gpu/lal_eam.cpp

@@ -303,7 +303,7 @@ double EAMT::host_memory_usage() const {
 }
 
 // ---------------------------------------------------------------------------
-// Copy nbor list from host if necessary and then compute atom energies/forces
+// Copy nbor list from host if necessary and then compute per-atom fp
 // ---------------------------------------------------------------------------
 template <class numtyp, class acctyp>
 void EAMT::compute(const int f_ago, const int inum_full, const int nlocal,
@@ -379,7 +379,7 @@ void EAMT::compute(const int f_ago, const int inum_full, const int nlocal,
 }
 
 // ---------------------------------------------------------------------------
-// Reneighbor on GPU and then compute per-atom densities
+// Reneighbor on GPU and then compute per-atom fp
 // ---------------------------------------------------------------------------
 template <class numtyp, class acctyp>
 int** EAMT::compute(const int ago, const int inum_full, const int nall,
@@ -461,7 +461,7 @@ int** EAMT::compute(const int ago, const int inum_full, const int nall,
 }
 
 // ---------------------------------------------------------------------------
-// Copy nbor list from host if necessary and then calculate forces, virials,..
+// Update per-atom fp, and then calculate forces, virials,..
 // ---------------------------------------------------------------------------
 template <class numtyp, class acctyp>
 void EAMT::compute2(int *ilist, const bool eflag, const bool vflag,

lib/gpu/lal_edpd.cu

@@ -401,10 +401,8 @@ __kernel void k_edpd_fast(const __global numtyp4 *restrict x_,
   __local numtyp4 sc[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
   __local numtyp4 kc[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
   __local numtyp sp_lj[4];
-  __local numtyp sp_sqrt[4];
   if (tid<4) {
     sp_lj[tid]=sp_lj_in[tid];
-    sp_sqrt[tid]=sp_sqrt_in[tid];
   }
   if (tid<MAX_SHARED_TYPES*MAX_SHARED_TYPES) {
     coeff[tid]=coeff_in[tid];

lib/gpu/lal_hippo.cu

@@ -1225,7 +1225,9 @@ __kernel void k_hippo_udirect2b(const __global numtyp4 *restrict x_,
   atom_info(t_per_atom,ii,tid,offset);
 
   int n_stride;
+#if (SHUFFLE_AVAIL == 0)
   local_allocate_store_charge();
+#endif
 
   acctyp _fieldp[6];
   for (int l=0; l<6; l++) _fieldp[l]=(acctyp)0;
@@ -1410,7 +1412,9 @@ __kernel void k_hippo_umutual2b(const __global numtyp4 *restrict x_,
   atom_info(t_per_atom,ii,tid,offset);
 
   int n_stride;
+#if (SHUFFLE_AVAIL == 0)
   local_allocate_store_charge();
+#endif
 
   acctyp _fieldp[6];
   for (int l=0; l<6; l++) _fieldp[l]=(acctyp)0;
@@ -2452,13 +2456,13 @@ __kernel void k_hippo_special15(__global int * dev_nbor,
                           const __global tagint *restrict special15,
                           const int inum, const int nall, const int nbor_pitch,
                           const int t_per_atom) {
-  int tid, ii, offset, n_stride, i;
+  int tid, ii, offset, n_stride, j;
   atom_info(t_per_atom,ii,tid,offset);
 
   if (ii<inum) {
 
     int numj, nbor, nbor_end;
-    nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
+    nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,j,numj,
               n_stride,nbor_end,nbor);
 
     int n15 = nspecial15[ii];
@@ -2467,7 +2471,7 @@ __kernel void k_hippo_special15(__global int * dev_nbor,
 
       int sj=dev_packed[nbor];
       int which = sj >> SBBITS & 3;
-      int j = sj & NEIGHMASK;
+      j = sj & NEIGHMASK;
       tagint jtag = tag[j];
 
       if (!which) {

lib/gpu/lal_mdpd.cu

@@ -210,13 +210,12 @@ __kernel void k_mdpd(const __global numtyp4 *restrict x_,
 
     const numtyp rhoi = extra[i].x;
 
-    numtyp factor_dpd, factor_sqrt;
+    numtyp factor_dpd;
     for ( ; nbor<nbor_end; nbor+=n_stride) {
       ucl_prefetch(dev_packed+nbor+n_stride);
 
       int j=dev_packed[nbor];
       factor_dpd = sp_lj[sbmask(j)];
-      factor_sqrt = sp_sqrt[sbmask(j)];
       j &= NEIGHMASK;
 
       numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
@@ -247,7 +246,6 @@ __kernel void k_mdpd(const __global numtyp4 *restrict x_,
         numtyp sigmaij = coeff[mtype].w;
         numtyp cutij =   coeff2[mtype].x;
         numtyp cut_rij = coeff2[mtype].y;
-        numtyp cutsqij = coeff2[mtype].z;
 
         numtyp wc = (numtyp)1.0 - r/cutij;
         numtyp wc_r = (numtyp)1.0 - r/cut_rij;
@@ -322,10 +320,8 @@ __kernel void k_mdpd_fast(const __global numtyp4 *restrict x_,
   __local numtyp4 coeff[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
   __local numtyp4 coeff2[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
   __local numtyp sp_lj[4];
-  __local numtyp sp_sqrt[4];
   if (tid<4) {
     sp_lj[tid]=sp_lj_in[tid];
-    sp_sqrt[tid]=sp_sqrt_in[tid];
   }
   if (tid<MAX_SHARED_TYPES*MAX_SHARED_TYPES) {
     coeff[tid]=coeff_in[tid];
@@ -339,7 +335,6 @@ __kernel void k_mdpd_fast(const __global numtyp4 *restrict x_,
   const numtyp sigmaij=coeff_in[ONETYPE].w;
   const numtyp cutij=coeff2_in[ONETYPE].x;
   const numtyp cut_rij=coeff2_in[ONETYPE].y;
-  const numtyp cutsqij=coeff2_in[ONETYPE].z;
   const numtyp cutsq_p=cutsq[ONETYPE];
   #endif
 
@@ -412,7 +407,6 @@ __kernel void k_mdpd_fast(const __global numtyp4 *restrict x_,
         numtyp sigmaij = coeff[mtype].w;
         numtyp cutij =   coeff2[mtype].x;
         numtyp cut_rij = coeff2[mtype].y;
-        numtyp cutsqij = coeff2[mtype].z;
         #endif
 
         numtyp wc = (numtyp)1.0 - r/cutij;

lib/gpu/lal_sph_heatconduction.cu

@@ -70,7 +70,9 @@ __kernel void k_sph_heatconduction(const __global numtyp4 *restrict x_,
   atom_info(t_per_atom,ii,tid,offset);
 
   int n_stride;
+#if (SHUFFLE_AVAIL == 0)
   local_allocate_store_pair();
+#endif
 
   acctyp dEacc = (acctyp)0;
 
@@ -171,7 +173,9 @@ __kernel void k_sph_heatconduction_fast(const __global numtyp4 *restrict x_,
   #endif
 
   int n_stride;
+#if (SHUFFLE_AVAIL == 0)
   local_allocate_store_pair();
+#endif
 
   acctyp dEacc = (acctyp)0;
 

lib/gpu/lal_sph_heatconduction.h

@@ -13,8 +13,8 @@
     email                : ndactrung@gmail.com
  ***************************************************************************/
 
-#ifndef LAL_SPH_LJ_H
+#ifndef LAL_SPH_HEATCONDUCTION_H
-#define LAL_SPH_LJ_H
+#define LAL_SPH_HEATCONDUCTION_H
 
 #include "lal_base_sph.h"
 

lib/gpu/lal_sph_heatconduction_ext.cpp

@@ -110,10 +110,10 @@ void sph_heatconduction_gpu_compute(const int ago, const int inum_full, const in
                         int **firstneigh, const bool eflag, const bool vflag,
                         const bool eatom, const bool vatom, int &host_start,
                         const double cpu_time, bool &success, tagint *host_tag,
-                        double **host_v, const int nlocal) {
+                        double **host_v) {
   SPHHeatConductionMF.compute(ago, inum_full, nall, host_x, host_type, ilist, numj,
                   firstneigh, eflag, vflag, eatom, vatom, host_start, cpu_time, success,
-                  host_tag, host_v, nlocal);
+                  host_tag, host_v);
 }
 
 void sph_heatconduction_gpu_get_extra_data(double *host_rho, double *host_esph) {

lib/gpu/lal_sph_lj_ext.cpp

@@ -110,10 +110,10 @@ void sph_lj_gpu_compute(const int ago, const int inum_full, const int nall,
                         int **firstneigh, const bool eflag, const bool vflag,
                         const bool eatom, const bool vatom, int &host_start,
                         const double cpu_time, bool &success, tagint *host_tag,
-                        double **host_v, const int nlocal) {
+                        double **host_v) {
   SPHLJMF.compute(ago, inum_full, nall, host_x, host_type, ilist, numj,
                   firstneigh, eflag, vflag, eatom, vatom, host_start, cpu_time, success,
-                  host_tag, host_v, nlocal);
+                  host_tag, host_v);
 }
 
 void sph_lj_gpu_get_extra_data(double *host_rho, double *host_esph, double *host_cv) {

lib/gpu/lal_sph_taitwater_ext.cpp

@@ -114,10 +114,10 @@ void sph_taitwater_gpu_compute(const int ago, const int inum_full, const int nal
                         int **firstneigh, const bool eflag, const bool vflag,
                         const bool eatom, const bool vatom, int &host_start,
                         const double cpu_time, bool &success, tagint *host_tag,
-                        double **host_v, const int nlocal) {
+                        double **host_v) {
   SPHTaitwaterMF.compute(ago, inum_full, nall, host_x, host_type, ilist, numj,
                   firstneigh, eflag, vflag, eatom, vatom, host_start, cpu_time, success,
-                  host_tag, host_v, nlocal);
+                  host_tag, host_v);
 }
 
 void sph_taitwater_gpu_get_extra_data(double *host_rho) {

lib/kokkos/core/src/Kokkos_Printf.hpp

@@ -31,7 +31,7 @@ namespace Kokkos {
 // backends. The GPU backends always return 1 and NVHPC only compiles if we
 // don't ask for the return value.
 template <typename... Args>
-KOKKOS_FUNCTION void printf(const char* format, Args... args) {
+KOKKOS_FORCEINLINE_FUNCTION void printf(const char* format, Args... args) {
 #ifdef KOKKOS_ENABLE_SYCL
   // Some compilers warn if "args" is empty and format is not a string literal
   if constexpr (sizeof...(Args) == 0)

lib/kokkos/core/src/OpenMP/Kokkos_OpenMP_Instance.hpp

@@ -219,8 +219,6 @@ KOKKOS_DEPRECATED void OpenMP::partition_master(F const& f, int num_partitions,
     Exec::validate_partition_impl(prev_instance->m_pool_size, num_partitions,
                                   partition_size);
 
-    OpenMP::memory_space space;
-
 #pragma omp parallel num_threads(num_partitions)
     {
       Exec thread_local_instance(partition_size);

lib/mdi/Install.py

@@ -32,7 +32,7 @@ make lib-mdi args="-m mpi" # build MDI lib with same settings as in the mpi Make
 
 # settings
 
-version = "1.4.16"
+version = "1.4.26"
 url = "https://github.com/MolSSI-MDI/MDI_Library/archive/v%s.tar.gz" % version
 
 # known checksums for different MDI versions. used to validate the download.
@@ -41,6 +41,7 @@ checksums = { \
               '1.4.12' : '7a222353ae8e03961d5365e6cd48baee', \
               '1.4.14' : '7a059bb12535360fdcb7de8402f9a0fc', \
               '1.4.16' : '407db44e2d79447ab5c1233af1965f65', \
+              '1.4.26' : '3124bb85259471e2a53a891f04bf697a', \
               }
 
 # print error message or help

potentials/CBNOH.aip.water.2dm

@@ -0,0 +1,58 @@
+# DATE: 2023-12-20 UNITS: metal CONTRIBUTOR: Wengen Ouyang w.g.ouyang@gmail.com
+# CITATION: Z. Feng, ..., and W. Ouyang, J. Phys. Chem. C 127(18), 8704 (2023).
+# CITATION: Z. Feng, ..., and W. Ouyang, Langmuir 39(50), 18198-18207 (2023).
+# Anisotropic Potential (AIP) for water/graphene and water/hBN heterojunctions
+# The parameters below are fitted against the PBE + MBD-NL (graphene/water) and SCAN (hBN/water) DFT reference data.
+
+# The parameters for bilayer graphene/graphene, graphene/hBN and hBN/hBN junctions are taken from
+# CITATION: Ouyang, Mandelli, Urbakh, Hod, Nano Letters 18, 6009-6016 (2018).
+#
+# -------------------Repulsion Potential ------------------++++++++++++++ Vdw Potential ++++++++++++++++************
+#    beta(A)    alpha    delta(A)  epsilon(meV)  C(meV)        d         sR      reff(A)   C6(meV*A^6)    S    rcut
+#
+# For water-graphene
+C Ow  5.453696  6.181724  1.250255   3.349092     0.687806  9.057065   1.232495   2.775772  100226.555031  1.0  2.0
+C Hw  2.553809  9.686644  1.964892  41.776171   -16.300128  9.015685   0.744155   2.415456    7409.128564  1.0  2.0
+Ow C  5.453696  6.181724  1.250255   3.349092     0.687806  9.057065   1.232495   2.775772  100226.555031  1.0  1.2
+Hw C  2.553809  9.686644  1.964892  41.776171   -16.300128  9.015685   0.744155   2.415456    7409.128564  1.0  1.2
+
+# For water-hBN
+N Ow  3.530598  16.377816  1.285374  1.717537    1.339337  24.797794   0.771411   3.928357   33589.850651  1.0  2.0
+N Hw  4.029390   5.360546  0.950352 15.945549   -1.486701  10.797276   1.352684   2.293775   41247.181447  1.0  2.0
+B Ow  3.907514   7.842519  2.380078 32.122737    1.190485  17.482482   0.788174   2.368217  139539.370785  1.0  2.0
+B Hw  3.804966   2.356248  1.114761  9.193309   -5.922514   9.000572   1.334703   1.746122   43796.489158  1.0  2.0
+Ow N  3.530598  16.377816  1.285374  1.717537    1.339337  24.797794   0.771411   3.928357   33589.850651  1.0  1.2
+Hw N  4.029390   5.360546  0.950352 15.945549   -1.486701  10.797276   1.352684   2.293775   41247.181447  1.0  1.2
+Ow B  3.907514   7.842519  2.380078 32.122737    1.190485  17.482482   0.788174   2.368217  139539.370785  1.0  1.2
+Hw B  3.804966   2.356248  1.114761  9.193309   -5.922514   9.000572   1.334703   1.746122   43796.489158  1.0  1.2
+
+# For graphene and hydrocarbons
+C C   3.205843  7.511126  1.235334  1.528338E-5  37.530428  15.499947  0.7954443  3.681440  25.714535E3   1.0   2.0
+H H   3.974540   6.53799  1.080633  0.6700556    0.8333833  15.022371  0.7490632  2.767223  1.6159581E3   1.0   1.2
+C H   2.642950  12.91410  1.020257  0.9750012    25.340996  15.222927  0.8115998  3.887324  5.6874617E3   1.0   1.5
+H C   2.642950  12.91410  1.020257  0.9750012    25.340996  15.222927  0.8115998  3.887324  5.6874617E3   1.0   1.5
+
+# For hBN
+B B   3.143737  9.825139  1.936405   2.7848400   14.495957   15.199263  0.7834022  3.682950  49.498013E3   1.0  2.0
+N N   3.443196  7.084490  1.747349   2.9139991   46.508553   15.020370  0.8008370  3.551843  14.810151E3   1.0  2.0
+B N   3.295257  7.224311  2.872667   1.3715032    0.4347152  14.594578  0.8044028  3.765728  24.669996E3   1.0  2.0
+B H   2.718657  9.214551  3.273063  14.015714    14.760509   15.084752  0.7768383  3.640866  7.9642467E3   1.0  1.5
+N B   3.295257  7.224311  2.872667   1.3715032    0.4347152  14.594578  0.8044028  3.765728  24.669996E3   1.0  2.0
+H B   2.718657  9.214551  3.273063  14.015714    14.760509   15.084752  0.7768383  3.640866  7.9642467E3   1.0  1.5
+
+# For graphene-hBN
+C B  3.303662  10.54415   2.926741  16.719972     0.3571734  15.305254  0.7001581  3.097327  30.162869E3   1.0  2.0
+C N  3.253564   8.825921  1.059550  18.344740    21.913573   15.000000  0.7234983  3.013117  19.063095E3   1.0  2.0
+B C  3.303662  10.54415   2.926741  16.719972     0.3571734  15.305254  0.7001581  3.097327  30.162869E3   1.0  2.0
+N C  3.253564   8.825921  1.059550  18.344740    21.913573   15.000000  0.7234983  3.013117  19.063095E3   1.0  2.0
+
+# The AIPs for other elements are turned off
+H Ow  5.453696   6.181724  1.250255  0.000000   0.000000    9.057065   1.232495   2.775772      0.000000   1.0  1.2
+H Hw  5.453696   6.181724  1.250255  0.000000   0.000000    9.057065   1.232495   2.775772      0.000000   1.0  1.2
+Ow H  5.453696   6.181724  1.250255  0.000000   0.000000    9.057065   1.232495   2.775772      0.000000   1.0  1.2
+Hw H  5.453696   6.181724  1.250255  0.000000   0.000000    9.057065   1.232495   2.775772      0.000000   1.0  1.2
+
+Ow Ow 5.453696   6.181724  1.250255  0.000000   0.000000    9.057065   1.232495   2.775772      0.000000   1.0  1.2
+Hw Hw 5.453696   6.181724  1.250255  0.000000   0.000000    9.057065   1.232495   2.775772      0.000000   1.0  1.2
+Ow Hw 5.453696   6.181724  1.250255  0.000000   0.000000    9.057065   1.232495   2.775772      0.000000   1.0  1.2
+Hw Ow 5.453696   6.181724  1.250255  0.000000   0.000000    9.057065   1.232495   2.775772      0.000000   1.0  1.2

potentials/COH.aip.water.2dm

@@ -1,28 +0,0 @@
-# DATE: 2022-12-02 UNITS: metal CONTRIBUTOR: Wengen Ouyang w.g.ouyang@gmail.com CITATION: Z. Feng, ..., and W. Ouyang, J. Phys. Chem. C 127, 8704 (2023).
-# Anisotropic Interfacial Potential (AIP) parameters for water/graphene heterojunctions
-# The parameters below are fitted against the PBE + MBD-NL DFT reference data from 2.5 A to 15 A.
-#
-# ----------------- Repulsion Potential ------------------++++++++++++++ Vdw Potential ++++++++++++++++************
-#    beta(A)    alpha    delta(A)  epsilon(meV)  C(meV)        d         sR      reff(A)   C6(meV*A^6)    S    rcut
-# For graphene and hydrocarbons
-C C  3.205843  7.511126  1.235334  1.528338E-5  37.530428  15.499947  0.7954443  3.681440  25.714535E3   1.0   2.0
-H H  3.974540   6.53799  1.080633  0.6700556    0.8333833  15.022371  0.7490632  2.767223  1.6159581E3   1.0   1.2
-C H  2.642950  12.91410  1.020257  0.9750012    25.340996  15.222927  0.8115998  3.887324  5.6874617E3   1.0   1.5
-H C  2.642950  12.91410  1.020257  0.9750012    25.340996  15.222927  0.8115998  3.887324  5.6874617E3   1.0   1.5
-
-# For water-graphene
-C Ow  5.45369612 6.18172364 1.25025450  3.34909245   0.68780636 9.05706482 1.23249498 2.77577173 100226.55503127   1.0   2.0
-C Hw  2.55380862 9.68664390 1.96489198 41.77617053 -16.30012807 9.01568534 0.74415463 2.41545571   7409.12856378   1.0   2.0
-Ow C  5.45369612 6.18172364 1.25025450  3.34909245   0.68780636 9.05706482 1.23249498 2.77577173 100226.55503127   1.0   1.2
-Hw C  2.55380862 9.68664390 1.96489198 41.77617053 -16.30012807 9.01568534 0.74415463 2.41545571   7409.12856378   1.0   1.2
-
-# # The ILPs for other systems are set to zero
-H Ow  5.45369612 6.18172364 1.25025450 0.00000000 0.00000000 9.05706482 1.23249498 2.77577173 0.00000000    1.0     1.2
-H Hw  5.45369612 6.18172364 1.25025450 0.00000000 0.00000000 9.05706482 1.23249498 2.77577173 0.00000000    1.0     1.2
-Ow H  5.45369612 6.18172364 1.25025450 0.00000000 0.00000000 9.05706482 1.23249498 2.77577173 0.00000000    1.0     1.2
-Hw H  5.45369612 6.18172364 1.25025450 0.00000000 0.00000000 9.05706482 1.23249498 2.77577173 0.00000000    1.0     1.2
-
-Ow Ow 5.45369612 6.18172364 1.25025450 0.00000000 0.00000000 9.05706482 1.23249498 2.77577173 0.00000000    1.0     1.2
-Hw Hw 5.45369612 6.18172364 1.25025450 0.00000000 0.00000000 9.05706482 1.23249498 2.77577173 0.00000000    1.0     1.2
-Ow Hw 5.45369612 6.18172364 1.25025450 0.00000000 0.00000000 9.05706482 1.23249498 2.77577173 0.00000000    1.0     1.2
-Hw Ow 5.45369612 6.18172364 1.25025450 0.00000000 0.00000000 9.05706482 1.23249498 2.77577173 0.00000000    1.0     1.2

potentials/TMD.ILP

@@ -0,0 +1,25 @@
+# DATE: 2021-12-02 UNITS: metal CONTRIBUTOR: Wengen Ouyang w.g.ouyang@gmail.com
+# CITATION: W. Ouyang, et al., J. Chem. Theory Comput. 17, 7237 (2021).
+# CITATION: W.  Jiang, et al., J. Phys. Chem. A, 127, 46, 9820–9830 (2023).
+# Interlayer Potential (ILP) for bilayer and bulk Group-VI Transition Metal Dichalcogenides.
+# The parameters below are fitted against the HSE + MBD-NL DFT reference data.
+#
+# -------------------- Repulsion Potential -------------------++++++++++++++++ Vdw Potential ++++++++++++++++*********
+#       beta(A)     alpha    delta(A)  epsilon(meV)  C(meV)         d         sR     reff(A)    C6(meV*A^6)    S   rcut
+Mo Mo  5.579450   9.377662   2.027222  144.151775   97.978570   89.437597  2.059031  5.122055  491850.316195  1.0  4.0
+W   W  5.530854   6.624992   1.983208    0.271792  140.174059  107.392585  1.356333  4.437591  691850.243962  1.0  4.0
+S   S  3.161402   8.093263   1.953140    4.586764  118.065466   58.809416  0.215367  4.299600  148811.243409  1.0  4.0
+Se Se  3.938627  10.515924   2.415783    3.012583   22.400612  116.864517  0.151121  5.884241  112506.195626  1.0  4.0
+Mo  W  5.412298   8.647128   2.108665   51.177950  184.342860  201.281256  2.547743  2.492287   99996.913401  1.0  4.0
+Mo  S  3.627152  19.971375   7.585031   76.101931    3.317496   45.720328  0.947470  4.410425  150597.857716  1.0  4.0
+Mo Se  6.196447   4.844134  14.362005    7.407221    0.058823   27.156223  0.976771  3.979186  786029.840651  1.0  4.0
+W   S  3.680136  11.163004  32.254117  110.019679   79.381335  138.340438  0.900750  8.875776  250600.809034  1.0  4.0
+W  Se  3.559392  20.638856   1.202717   20.478669  197.422484   10.005271  1.052738  3.815817  288321.561114  1.0  4.0
+S  Se  2.820092   7.491151   1.933323  141.532559  293.127817   90.470904  0.390492  4.170885  117688.987069  1.0  4.0
+# Symmetric Atom Pair
+W  Mo  5.412298   8.647128   2.108665   51.177950  184.342860  201.281256  2.547743  2.492287   99996.913401  1.0  4.0
+S  Mo  3.627152  19.971375   7.585031   76.101931    3.317496   45.720328  0.947470  4.410425  150597.857716  1.0  4.0
+Se Mo  6.196447   4.844134  14.362005    7.407221    0.058823   27.156223  0.976771  3.979186  786029.840651  1.0  4.0
+S   W  3.680136  11.163004  32.254117  110.019679   79.381335  138.340438  0.900750  8.875776  250600.809034  1.0  4.0
+Se  W  3.559392  20.638856   1.202717   20.478669  197.422484   10.005271  1.052738  3.815817  288321.561114  1.0  4.0
+Se  S  2.820092   7.491151   1.933323  141.532559  293.127817   90.470904  0.390492  4.170885  117688.987069  1.0  4.0

python/lammps/core.py

@@ -282,6 +282,8 @@ class lammps(object):
     self.lib.lammps_config_accelerator.argtypes = [c_char_p, c_char_p, c_char_p]
 
     self.lib.lammps_set_variable.argtypes = [c_void_p, c_char_p, c_char_p]
+    self.lib.lammps_set_string_variable.argtypes = [c_void_p, c_char_p, c_char_p]
+    self.lib.lammps_set_internal_variable.argtypes = [c_void_p, c_char_p, c_double]
 
     self.lib.lammps_has_style.argtypes = [c_void_p, c_char_p, c_char_p]
 
@@ -1252,6 +1254,8 @@ class lammps(object):
   def set_variable(self,name,value):
     """Set a new value for a LAMMPS string style variable
 
+    .. deprecated:: TBD
+
     This is a wrapper around the :cpp:func:`lammps_set_variable`
     function of the C-library interface.
 
@@ -1271,6 +1275,52 @@ class lammps(object):
 
   # -------------------------------------------------------------------------
 
+  def set_string_variable(self,name,value):
+    """Set a new value for a LAMMPS string style variable
+
+    .. versionadded:: TBD
+
+    This is a wrapper around the :cpp:func:`lammps_set_string_variable`
+    function of the C-library interface.
+
+    :param name: name of the variable
+    :type name: string
+    :param value: new variable value
+    :type value: any. will be converted to a string
+    :return: either 0 on success or -1 on failure
+    :rtype: int
+    """
+    if name: name = name.encode()
+    else: return -1
+    if value: value = str(value).encode()
+    else: return -1
+    with ExceptionCheck(self):
+      return self.lib.lammps_set_string_variable(self.lmp,name,value)
+
+  # -------------------------------------------------------------------------
+
+  def set_internal_variable(self,name,value):
+    """Set a new value for a LAMMPS internal style variable
+
+    .. versionadded:: TBD
+
+    This is a wrapper around the :cpp:func:`lammps_set_internal_variable`
+    function of the C-library interface.
+
+    :param name: name of the variable
+    :type name: string
+    :param value: new variable value
+    :type value: float or compatible. will be converted to float
+    :return: either 0 on success or -1 on failure
+    :rtype: int
+    """
+    if name: name = name.encode()
+    else: return -1
+    with ExceptionCheck(self):
+      return self.lib.lammps_set_internal_variable(self.lmp,name,value)
+
+  # -------------------------------------------------------------------------
+
   # return vector of atom properties gathered across procs
   # 3 variants to match src/library.cpp
   # name = atom property recognized by LAMMPS in atom->extract()

src/.gitignore

@@ -348,6 +348,8 @@
 /compute_nbond_atom.h
 /fix_nve_bpm_sphere.cpp
 /fix_nve_bpm_sphere.h
+/fix_update_special_bonds.cpp
+/fix_update_special_bonds.h
 /pair_bpm_spring.cpp
 /pair_bpm_spring.h
 

src/ADIOS/reader_adios.cpp

@@ -30,7 +30,7 @@
 
 using namespace LAMMPS_NS;
 
-#define SMALL 1.0e-6
+static constexpr double SMALL = 1.0e-6;
 
 namespace LAMMPS_NS {
 class ReadADIOSInternal {

src/AMOEBA/amoeba_convolution.cpp

@@ -47,13 +47,8 @@ enum{MPOLE_GRID,POLAR_GRID,POLAR_GRIDC,DISP_GRID,INDUCE_GRID,INDUCE_GRIDC};
 //#define SCALE 1
 #define SCALE 0
 
-#ifdef FFT_SINGLE
+static constexpr FFT_SCALAR ZEROF = 0.0;
-#define ZEROF 0.0f
+static constexpr FFT_SCALAR ONEF =  1.0;
-#define ONEF  1.0f
-#else
-#define ZEROF 0.0
-#define ONEF  1.0
-#endif
 
 /* ----------------------------------------------------------------------
    partition an FFT grid across processors

src/AMOEBA/amoeba_file.cpp

@@ -25,23 +25,23 @@
 
 using namespace LAMMPS_NS;
 
-enum{UNKNOWN,FFIELD,LITERATURE,ATOMTYPE,VDWL,VDWLPAIR,BSTRETCH,SBEND,ABEND,
+enum { UNKNOWN, FFIELD, LITERATURE, ATOMTYPE, VDWL, VDWLPAIR, BSTRETCH, SBEND, ABEND,
-     PAULI,DISPERSION,UB,OUTPLANE,TORSION,PITORSION,ATOMMULT,
+  PAULI, DISPERSION, UB, OUTPLANE, TORSION, PITORSION, ATOMMULT, QPENETRATION, DIPPOLAR,
-     QPENETRATION,DIPPOLAR,QTRANSFER,END_OF_FILE};
+  QTRANSFER, END_OF_FILE };
-enum{ALLINGER,BUFFERED_14_7};
+enum { ALLINGER, BUFFERED_14_7 };
-enum{ARITHMETIC,GEOMETRIC,CUBIC_MEAN,R_MIN,SIGMA,DIAMETER,HARMONIC,HHG,W_H};
+enum { ARITHMETIC, GEOMETRIC, CUBIC_MEAN, R_MIN, SIGMA, DIAMETER, HARMONIC, HHG, W_H };
-enum{MUTUAL,OPT,TCG,DIRECT};
+enum { MUTUAL, OPT, TCG, DIRECT };
-enum{NOFRAME,ZONLY,ZTHENX,BISECTOR,ZBISECT,THREEFOLD};
+enum { NOFRAME, ZONLY, ZTHENX, BISECTOR, ZBISECT, THREEFOLD };
-enum{GEAR,ASPC,LSQR};
+enum { GEAR, ASPC, LSQR };
 
-#define MAXLINE 65536              // crazy big for TORSION-TORSION section
+static constexpr int MAXLINE = 65536;             // crazy big for TORSION-TORSION section
-#define MAX_TYPE_PER_GROUP 6       // max types per AMOEBA group
+static constexpr int MAX_TYPE_PER_GROUP = 6;     // max types per AMOEBA group
-#define MAX_FRAME_PER_TYPE 32      // max multipole frames for any AMOEBA type
+static constexpr int MAX_FRAME_PER_TYPE = 32;    // max multipole frames for any AMOEBA type
 
-#define DELTA_TYPE_CLASS 32
+static constexpr int DELTA_TYPE_CLASS = 32;
-#define DELTA_VDWL_PAIR 16
+static constexpr int DELTA_VDWL_PAIR = 16;
 
-#define BOHR 0.52917721067         // Bohr in Angstroms
+static constexpr double BOHR = 0.52917721067;    // Bohr in Angstroms
 
 // methods to read, parse, and store info from force field file
 
@@ -79,7 +79,7 @@ void PairAmoeba::read_prmfile(char *filename)
 
   int me = comm->me;
   FILE *fptr;
-  char line[MAXLINE];
+  char line[MAXLINE] = {'\0'};
 
   if (me == 0) {
     fptr = utils::open_potential(filename, lmp, nullptr);
@@ -179,8 +179,7 @@ void PairAmoeba::read_prmfile(char *filename)
       for (int i = 1; i <= n_amtype; i++) nmultiframe[i] = 0;
     }
 
-    char next[MAXLINE];
+    char next[MAXLINE] = {'\0'};
-    next[0] = '\0';
     bool has_next = false;
     int n;
     while (true) {
@@ -381,7 +380,7 @@ void PairAmoeba::read_keyfile(char *filename)
 
   int me = comm->me;
   FILE *fptr;
-  char line[MAXLINE];
+  char line[MAXLINE] = {'\0'};
   if (me == 0) {
     fptr = utils::open_potential(filename, lmp, nullptr);
     if (fptr == nullptr)

src/AMOEBA/amoeba_induce.cpp

@@ -41,7 +41,7 @@ enum{GEAR,ASPC,LSQR};
 enum{BUILD,APPLY};
 enum{GORDON1,GORDON2};
 
-#define DEBYE 4.80321    // conversion factor from q-Angs (real units) to Debye
+static constexpr double DEBYE = 4.80321;    // conversion factor from q-Angs (real units) to Debye
 
 /* ----------------------------------------------------------------------
    induce = induced dipole moments via pre-conditioned CG solver

src/AMOEBA/amoeba_multipole.cpp

@@ -35,11 +35,11 @@ enum{FIELD,ZRSD,TORQUE,UFLD};                          // reverse comm
 enum{VDWL,REPULSE,QFER,DISP,MPOLE,POLAR,USOLV,DISP_LONG,MPOLE_LONG,POLAR_LONG};
 
 #ifdef FFT_SINGLE
-#define ZEROF 0.0f
+static constexpr FFT_SCALAR ZEROF = 0.0f;
-#define ONEF  1.0f
+static constexpr FFT_SCALAR ONEF =  1.0f;
 #else
-#define ZEROF 0.0
+static constexpr FFT_SCALAR ZEROF = 0.0;
-#define ONEF  1.0
+static constexpr FFT_SCALAR ONEF =  1.0;
 #endif
 
 /* ----------------------------------------------------------------------

src/AMOEBA/angle_amoeba.cpp

@@ -30,7 +30,7 @@
 using namespace LAMMPS_NS;
 using namespace MathConst;
 
-#define SMALL 0.001
+static constexpr double SMALL = 0.001;
 
 /* ---------------------------------------------------------------------- */
 

src/AMOEBA/fix_amoeba_bitorsion.cpp

@@ -32,10 +32,10 @@ using namespace LAMMPS_NS;
 using namespace FixConst;
 using namespace MathConst;
 
-#define BITORSIONMAX 6   // max # of BiTorsion terms stored by one atom
+static constexpr int BITORSIONMAX = 6;   // max # of BiTorsion terms stored by one atom
-#define LISTDELTA 10000
+static constexpr int LISTDELTA = 10000;
-#define LB_FACTOR 1.5
+static constexpr double LB_FACTOR = 1.5;
-#define MAXLINE 1024
+static constexpr int MAXLINE = 1024;
 
 // spline weighting factors
 
@@ -724,7 +724,7 @@ double FixAmoebaBiTorsion::compute_scalar()
 
 void FixAmoebaBiTorsion::read_grid_data(char *bitorsion_file)
 {
-  char line[MAXLINE];
+  char line[MAXLINE] = {'\0'};
   char *eof;
 
   FILE *fp = nullptr;

src/AMOEBA/fix_amoeba_pitorsion.cpp

@@ -32,9 +32,9 @@ using namespace LAMMPS_NS;
 using namespace FixConst;
 using namespace MathConst;
 
-#define PITORSIONMAX 6   // max # of PiTorsion terms stored by one atom
+static constexpr int PITORSIONMAX = 6;   // max # of PiTorsion terms stored by one atom
-#define LISTDELTA 8196
+static constexpr int LISTDELTA = 8196;
-#define LB_FACTOR 1.5
+static constexpr double LB_FACTOR = 1.5;
 
 /* ---------------------------------------------------------------------- */
 

src/AMOEBA/improper_amoeba.cpp

@@ -28,8 +28,8 @@
 using namespace LAMMPS_NS;
 using namespace MathConst;
 
-#define TOLERANCE 0.05
+static constexpr double TOLERANCE = 0.05;
-#define SMALL     0.001
+static constexpr double SMALL =     0.001;
 
 /* ---------------------------------------------------------------------- */
 

src/AMOEBA/pair_amoeba.cpp

@@ -1,4 +1,3 @@
-// clang-format off
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    https://www.lammps.org/ Sandia National Laboratories
@@ -39,17 +38,19 @@ using namespace LAMMPS_NS;
 
 using MathSpecial::powint;
 
-enum{INDUCE,RSD,SETUP_AMOEBA,SETUP_HIPPO,KMPOLE,AMGROUP,PVAL};  // forward comm
+enum { INDUCE, RSD, SETUP_AMOEBA, SETUP_HIPPO, KMPOLE, AMGROUP, PVAL };    // forward comm
-enum{FIELD,ZRSD,TORQUE,UFLD};                                   // reverse comm
+enum { FIELD, ZRSD, TORQUE, UFLD };                                        // reverse comm
-enum{ARITHMETIC,GEOMETRIC,CUBIC_MEAN,R_MIN,SIGMA,DIAMETER,HARMONIC,HHG,W_H};
+enum { ARITHMETIC, GEOMETRIC, CUBIC_MEAN, R_MIN, SIGMA, DIAMETER, HARMONIC, HHG, W_H };
-enum{HAL,REPULSE,QFER,DISP,MPOLE,POLAR,USOLV,DISP_LONG,MPOLE_LONG,POLAR_LONG};
+enum { HAL, REPULSE, QFER, DISP, MPOLE, POLAR, USOLV, DISP_LONG, MPOLE_LONG, POLAR_LONG };
-enum{MPOLE_GRID,POLAR_GRID,POLAR_GRIDC,DISP_GRID,INDUCE_GRID,INDUCE_GRIDC};
+enum { MPOLE_GRID, POLAR_GRID, POLAR_GRIDC, DISP_GRID, INDUCE_GRID, INDUCE_GRIDC };
-enum{MUTUAL,OPT,TCG,DIRECT};
+enum { MUTUAL, OPT, TCG, DIRECT };
-enum{GEAR,ASPC,LSQR};
+enum { GEAR, ASPC, LSQR };
 
-#define DELTASTACK 16
+static constexpr int DELTASTACK = 16;
 #define DEBUG_AMOEBA 0
 
+// clang-format off
+
 /* ---------------------------------------------------------------------- */
 
 PairAmoeba::PairAmoeba(LAMMPS *lmp) : Pair(lmp)
@@ -429,7 +430,7 @@ void PairAmoeba::compute(int eflag, int vflag)
   // output FF settings to screen and logfile
   //   delay until here because RMS force accuracy is computed based on rpole
 
-  if (first_flag_compute && (comm->me == 0)) print_settings();
+  if (first_flag_compute) print_settings();
   first_flag_compute = 0;
 
   if (amoeba) pbc_xred();
@@ -827,28 +828,36 @@ void PairAmoeba::init_style()
 
   // check if all custom atom arrays were set via fix property/atom
 
-  int flag,cols;
+  // clang-format on
+  const char *names[6] = {"amtype", "amgroup", "redID", "xyzaxis", "polaxe", "pval"};
+  const int flag_check[6] = {0, 0, 1, 1, 0, 1};     // correct type (0 int, 1 dbl)
+  const int cols_check[6] = {0, 0, 0, 3, 0, 0};     // xyzaxis 3 cols, all others 0
+  const int ghost_check[6] = {1, 1, 1, 0, 0, 1};    // which types need ghost; TO-DO: check
+  int flag, cols, ghost, index[6];
 
-  index_amtype = atom->find_custom("amtype",flag,cols);
+  // clang-format off
-  if (index_amtype < 0 || flag || cols)
-    error->all(FLERR,"Pair {} amtype is not defined", mystyle);
-  index_amgroup = atom->find_custom("amgroup",flag,cols);
-  if (index_amgroup < 0 || flag || cols)
-    error->all(FLERR,"Pair {} amgroup is not defined", mystyle);
 
-  index_redID = atom->find_custom("redID",flag,cols);
+  for (int i = 0; i < 6; i++) {
-  if (index_redID < 0 || !flag || cols)
+    if (ghost_check[i]) {
-    error->all(FLERR,"Pair {} redID is not defined", mystyle);
+      index[i] = atom->find_custom_ghost(names[i], flag, cols, ghost);
-  index_xyzaxis = atom->find_custom("xyzaxis",flag,cols);
+    } else {
-  if (index_xyzaxis < 0 || !flag || cols == 0)
+      index[i] = atom->find_custom(names[i], flag, cols);
-    error->all(FLERR,"Pair {} xyzaxis is not defined", mystyle);
+    }
+    std::string err = "";
+    if (index[i] < 0) err = "was not defined";
+    else if (flag_check[i] != flag) err = "has the wrong type";
+    else if (cols_check[i] != cols) err = "has the wrong number of columns";
+    else if (ghost_check[i] && !ghost) err = "must be set by fix property/atom with ghost yes";
+    if (err != "")
+      error->all(FLERR,"Pair {} per-atom variable {} {}", mystyle, names[i], err);
+  }
 
-  index_polaxe = atom->find_custom("polaxe",flag,cols);
+  index_amtype  = index[0];
-  if (index_polaxe < 0 || flag || cols)
+  index_amgroup = index[1];
-    error->all(FLERR,"Pair {} polaxe is not defined", mystyle);
+  index_redID   = index[2];
-  index_pval = atom->find_custom("pval",flag,cols);
+  index_xyzaxis = index[3];
-  if (index_pval < 0 || !flag || cols)
+  index_polaxe  = index[4];
-    error->all(FLERR,"Pair {} pval is not defined", mystyle);
+  index_pval    = index[5];
 
   // -------------------------------------------------------------------
   // one-time initializations
@@ -1069,7 +1078,14 @@ void PairAmoeba::init_style()
 void PairAmoeba::print_settings()
 {
   std::string mesg = utils::uppercase(mystyle) + " force field settings\n";
+  double estimated_mpole_accuracy = 0.0;
 
+  if (use_ewald) {
+    choose(MPOLE_LONG);
+    estimated_mpole_accuracy = final_accuracy_mpole();
+  }
+
+  if (comm->me == 0) {
     if (amoeba) {
       choose(HAL);
       mesg += fmt::format("  hal: cut {} taper {} vscale {} {} {} {}\n", sqrt(off2),sqrt(cut2),
@@ -1098,13 +1114,12 @@ void PairAmoeba::print_settings()
 
     if (use_ewald) {
       choose(MPOLE_LONG);
-    double estimated_accuracy = final_accuracy_mpole();
+      mesg += fmt::format("  multipole: cut {} aewald {} bsorder {} FFT {} {} {}\n"
-    mesg += fmt::format("  multipole: cut {} aewald {} bsorder {} FFT {} {} {}; "
+                          "             estimated absolute RMS force accuracy = {:.8g}\n"
-                        "estimated absolute RMS force accuracy = {:.8g}; "
+                          "             estimated relative RMS force accuracy = {:.8g}\n"
-                        "estimated relative RMS force accuracy = {:.8g}; "
+                          "             mscale {} {} {} {}\n",
-                        "mscale {} {} {} {}\n",
                           sqrt(off2),aewald,bseorder,nefft1,nefft2,nefft3,
-                        estimated_accuracy,estimated_accuracy/two_charge_force,
+                          estimated_mpole_accuracy,estimated_mpole_accuracy/two_charge_force,
                           special_mpole[1],special_mpole[2],special_mpole[3],special_mpole[4]);
     } else {
       choose(MPOLE);
@@ -1142,6 +1157,7 @@ void PairAmoeba::print_settings()
     choose(USOLV);
     mesg += fmt::format("  precondition: cut {}\n",sqrt(off2));
     utils::logmesg(lmp, mesg);
+  }
 }
 
 /* ----------------------------------------------------------------------

src/ASPHERE/compute_temp_asphere.cpp

@@ -33,7 +33,7 @@ using namespace LAMMPS_NS;
 
 enum{ROTATE,ALL};
 
-#define INERTIA 0.2          // moment of inertia prefactor for ellipsoid
+static constexpr double INERTIA = 0.2;          // moment of inertia prefactor for ellipsoid
 
 /* ---------------------------------------------------------------------- */
 

src/ASPHERE/fix_nve_asphere.cpp

@@ -26,7 +26,7 @@
 using namespace LAMMPS_NS;
 using namespace FixConst;
 
-#define INERTIA 0.2          // moment of inertia prefactor for ellipsoid
+static constexpr double INERTIA = 0.2;          // moment of inertia prefactor for ellipsoid
 
 /* ---------------------------------------------------------------------- */
 

src/ASPHERE/fix_nve_line.cpp

@@ -23,7 +23,7 @@ using namespace LAMMPS_NS;
 using namespace FixConst;
 using namespace MathConst;
 
-#define INERTIA (1.0/12.0)     // moment of inertia prefactor for line segment
+static constexpr double INERTIA = (1.0/12.0);     // moment of inertia prefactor for line segment
 
 /* ---------------------------------------------------------------------- */