Merge branch 'develop' of github.com:lammps/lammps into fix_viscous_kokkos

.github/CODEOWNERS

@@ -37,6 +37,7 @@ src/MESONT/*          @iafoss
 src/ML-HDNNP/*        @singraber
 src/ML-IAP/*          @athomps
 src/ML-PACE/*         @yury-lysogorskiy
+src/ML-POD/*          @exapde @rohskopf
 src/MOFFF/*           @hheenen
 src/MOLFILE/*         @akohlmey
 src/NETCDF/*          @pastewka
@@ -63,6 +64,8 @@ src/MANYBODY/pair_atm.*             @sergeylishchuk
 src/REPLICA/*_grem.*                @dstelter92
 src/EXTRA-COMPUTE/compute_stress_mop*.* @RomainVermorel
 src/MISC/*_tracker.*                @jtclemm
+src/MC/fix_gcmc.*                   @athomps
+src/MC/fix_sgcmc.*                  @athomps
 
 # core LAMMPS classes
 src/lammps.*              @sjplimp

.github/workflows/compile-msvc.yml

@@ -26,7 +26,7 @@ jobs:
     - name: Select Python version
       uses: actions/setup-python@v4
       with:
-        python-version: '3.10'
+        python-version: '3.11'
 
     - name: Building LAMMPS via CMake
       shell: bash
@@ -37,6 +37,8 @@ jobs:
         nuget install MSMPIDIST
         cmake -C cmake/presets/windows.cmake \
               -D PKG_PYTHON=on \
+              -D WITH_PNG=off \
+              -D WITH_JPEG=off \
               -S cmake -B build \
               -D BUILD_SHARED_LIBS=on \
               -D LAMMPS_EXCEPTIONS=on \

cmake/CMakeLists.txt

@@ -266,6 +266,7 @@ set(STANDARD_PACKAGES
   ML-QUIP
   ML-RANN
   ML-SNAP
+  ML-POD  
   MOFFF
   MOLECULE
   MOLFILE
@@ -432,7 +433,7 @@ if(BUILD_OMP)
   target_link_libraries(lmp PRIVATE OpenMP::OpenMP_CXX)
 endif()
 
-if(PKG_MSCG OR PKG_ATC OR PKG_AWPMD OR PKG_ML-QUIP OR PKG_LATTE OR PKG_ELECTRODE)
+if(PKG_MSCG OR PKG_ATC OR PKG_AWPMD OR PKG_ML-QUIP OR PKG_ML-POD OR PKG_LATTE OR PKG_ELECTRODE)
   enable_language(C)
   if (NOT USE_INTERNAL_LINALG)
     find_package(LAPACK)
@@ -638,7 +639,7 @@ foreach(PKG_LIB POEMS ATC AWPMD H5MD MESONT)
   endif()
 endforeach()
 
-if(PKG_ELECTRODE)
+if(PKG_ELECTRODE OR PKG_ML-POD)
   target_link_libraries(lammps PRIVATE ${LAPACK_LIBRARIES})
 endif()
 
@@ -667,7 +668,7 @@ endif()
 # packages which selectively include variants based on enabled styles
 # e.g. accelerator packages
 ######################################################################
-foreach(PKG_WITH_INCL CORESHELL DPD-SMOOTH MISC PHONON QEQ OPENMP KOKKOS OPT INTEL GPU)
+foreach(PKG_WITH_INCL CORESHELL DPD-SMOOTH MC MISC PHONON QEQ OPENMP KOKKOS OPT INTEL GPU)
   if(PKG_${PKG_WITH_INCL})
     include(Packages/${PKG_WITH_INCL})
   endif()

cmake/Modules/Packages/MC.cmake

@@ -0,0 +1,9 @@
+# fix sgcmc may only be installed if also the EAM pair style from MANYBODY is installed
+if(NOT PKG_MANYBODY)
+  get_property(LAMMPS_FIX_HEADERS GLOBAL PROPERTY FIX)
+  list(REMOVE_ITEM LAMMPS_FIX_HEADERS ${LAMMPS_SOURCE_DIR}/MC/fix_sgcmc.h)
+  set_property(GLOBAL PROPERTY FIX "${LAMMPS_FIX_HEADERS}")
+  get_target_property(LAMMPS_SOURCES lammps SOURCES)
+  list(REMOVE_ITEM LAMMPS_SOURCES ${LAMMPS_SOURCE_DIR}/MC/fix_sgcmc.cpp)
+  set_property(TARGET lammps PROPERTY SOURCES "${LAMMPS_SOURCES}")
+endif()

cmake/presets/all_off.cmake

@@ -56,6 +56,7 @@ set(ALL_PACKAGES
   ML-HDNNP
   ML-IAP
   ML-PACE
+  ML-POD
   ML-QUIP
   ML-RANN
   ML-SNAP

cmake/presets/all_on.cmake

@@ -58,6 +58,7 @@ set(ALL_PACKAGES
   ML-HDNNP
   ML-IAP
   ML-PACE
+  ML-POD
   ML-QUIP
   ML-RANN
   ML-SNAP

cmake/presets/mingw-cross.cmake

@@ -47,6 +47,7 @@ set(WIN_PACKAGES
   MISC
   ML-HDNNP
   ML-IAP
+  ML-POD
   ML-RANN
   ML-SNAP
   MOFFF

cmake/presets/most.cmake

@@ -41,6 +41,7 @@ set(ALL_PACKAGES
   MEAM
   MISC
   ML-IAP
+  ML-POD
   ML-SNAP
   MOFFF
   MOLECULE

doc/src/Build_extras.rst

@@ -36,6 +36,7 @@ This is the list of packages that may require additional steps.
    * :ref:`AWPMD <awpmd>`
    * :ref:`COLVARS <colvars>`
    * :ref:`COMPRESS <compress>`
+   * :ref:`ELECTRODE <electrode>`
    * :ref:`GPU <gpu>`
    * :ref:`H5MD <h5md>`
    * :ref:`INTEL <intel>`
@@ -48,6 +49,7 @@ This is the list of packages that may require additional steps.
    * :ref:`ML-HDNNP <ml-hdnnp>`
    * :ref:`ML-IAP <mliap>`
    * :ref:`ML-PACE <ml-pace>`
+   * :ref:`ML-POD <ml-pod>`
    * :ref:`ML-QUIP <ml-quip>`
    * :ref:`MOLFILE <molfile>`
    * :ref:`MSCG <mscg>`
@@ -1411,6 +1413,49 @@ at: `https://github.com/ICAMS/lammps-user-pace/ <https://github.com/ICAMS/lammps
 
 ----------
 
+.. _ml-pod:
+
+ML-POD package
+-----------------------------
+
+.. tabs::
+
+   .. tab:: CMake build
+
+      No additional settings are needed besides ``-D PKG_ML-POD=yes``.
+
+   .. tab:: Traditional make
+
+      Before building LAMMPS, you must configure the ML-POD support
+      settings in ``lib/mlpod``.  You can do this manually, if you
+      prefer, or do it in one step from the ``lammps/src`` dir, using a
+      command like the following, which simply invoke the
+      ``lib/mlpod/Install.py`` script with the specified args:
+
+      .. code-block:: bash
+
+         $ make lib-mlpod                   # print help message
+         $ make lib-mlpod args="-m serial"  # build with GNU g++ compiler and MPI STUBS (settings as with "make serial")
+         $ make lib-mlpod args="-m mpi"     # build with default MPI compiler (settings as with "make mpi")
+         $ make lib-mlpod args="-m mpi -e linalg"   # same as above but use the bundled linalg lib
+
+      Note that the ``Makefile.lammps`` file has settings to use the BLAS
+      and LAPACK linear algebra libraries.  These can either exist on
+      your system, or you can use the files provided in ``lib/linalg``.
+      In the latter case you also need to build the library in
+      ``lib/linalg`` with a command like these:
+
+      .. code-block:: bash
+
+         $ make lib-linalg                     # print help message
+         $ make lib-linalg args="-m serial"    # build with GNU Fortran compiler (settings as with "make serial")
+         $ make lib-linalg args="-m mpi"       # build with default MPI Fortran compiler (settings as with "make mpi")
+         $ make lib-linalg args="-m gfortran"  # build with GNU Fortran compiler
+
+      The package itself is activated with ``make yes-ML-POD``.
+
+----------
+
 .. _plumed:
 
 PLUMED package

doc/src/Commands_all.rst

@@ -31,7 +31,6 @@ table above.
    * :doc:`bond_style <bond_style>`
    * :doc:`bond_write <bond_write>`
    * :doc:`boundary <boundary>`
-   * :doc:`box <box>`
    * :doc:`change_box <change_box>`
    * :doc:`clear <clear>`
    * :doc:`comm_modify <comm_modify>`
@@ -90,8 +89,7 @@ table above.
    * :doc:`region <region>`
    * :doc:`replicate <replicate>`
    * :doc:`rerun <rerun>`
-   * :doc:`reset_atom_ids <reset_atom_ids>`
+   * :doc:`reset_atoms <reset_atoms>`
-   * :doc:`reset_mol_ids <reset_mol_ids>`
    * :doc:`reset_timestep <reset_timestep>`
    * :doc:`restart <restart>`
    * :doc:`run <run>`
@@ -127,6 +125,7 @@ additional letter in parenthesis: k = KOKKOS.
    * :doc:`group2ndx <group2ndx>`
    * :doc:`hyper <hyper>`
    * :doc:`kim <kim_commands>`
+   * :doc:`fitpod <fitpod_command>`
    * :doc:`mdi <mdi>`
    * :doc:`ndx2group <group2ndx>`
    * :doc:`neb <neb>`

doc/src/Commands_category.rst

@@ -25,7 +25,6 @@ Setup simulation box:
    :columns: 4
 
    * :doc:`boundary <boundary>`
-   * :doc:`box <box>`
    * :doc:`change_box <change_box>`
    * :doc:`create_box <create_box>`
    * :doc:`dimension <dimension>`

doc/src/Commands_fix.rst

@@ -69,9 +69,9 @@ OPT.
    * :doc:`edpd/source <fix_dpd_source>`
    * :doc:`efield <fix_efield>`
    * :doc:`ehex <fix_ehex>`
-   * :doc:`electrode/conp (i) <fix_electrode_conp>`
+   * :doc:`electrode/conp (i) <fix_electrode>`
-   * :doc:`electrode/conq (i) <fix_electrode_conp>`
+   * :doc:`electrode/conq (i) <fix_electrode>`
-   * :doc:`electrode/thermo (i) <fix_electrode_conp>`
+   * :doc:`electrode/thermo (i) <fix_electrode>`
    * :doc:`electron/stopping <fix_electron_stopping>`
    * :doc:`electron/stopping/fit <fix_electron_stopping>`
    * :doc:`enforce2d (k) <fix_enforce2d>`
@@ -213,6 +213,7 @@ OPT.
    * :doc:`saed/vtk <fix_saed_vtk>`
    * :doc:`setforce (k) <fix_setforce>`
    * :doc:`setforce/spin <fix_setforce>`
+   * :doc:`sgcmc <fix_sgcmc>`
    * :doc:`shake (k) <fix_shake>`
    * :doc:`shardlow (k) <fix_shardlow>`
    * :doc:`smd <fix_smd>`

doc/src/Commands_pair.rst

@@ -237,6 +237,7 @@ OPT.
    * :doc:`oxrna2/coaxstk <pair_oxrna2>`
    * :doc:`pace (k) <pair_pace>`
    * :doc:`pace/extrapolation <pair_pace>`
+   * :doc:`pod <pair_pod>`
    * :doc:`peri/eps <pair_peri>`
    * :doc:`peri/lps (o) <pair_peri>`
    * :doc:`peri/pmb (o) <pair_peri>`

doc/src/Commands_removed.rst

@@ -2,14 +2,17 @@ Removed commands and packages
 =============================
 
 This page lists LAMMPS commands and packages that have been removed from
-the distribution and provides suggestions for alternatives or replacements.
+the distribution and provides suggestions for alternatives or
-LAMMPS has special dummy styles implemented, that will stop LAMMPS and
+replacements.  LAMMPS has special dummy styles implemented, that will
-print a suitable error message in most cases, when a style/command is used
+stop LAMMPS and print a suitable error message in most cases, when a
-that has been removed.
+style/command is used that has been removed or will replace the command
+with the direct alternative (if available) and print a warning.
 
 Fix ave/spatial and fix ave/spatial/sphere
 ------------------------------------------
 
+.. deprecated:: 11Dec2015
+
 The fixes ave/spatial and ave/spatial/sphere have been removed from LAMMPS
 since they were superseded by the more general and extensible "chunk
 infrastructure".  Here the system is partitioned in one of many possible
@@ -17,10 +20,23 @@ ways through the :doc:`compute chunk/atom <compute_chunk_atom>` command
 and then averaging is done using :doc:`fix ave/chunk <fix_ave_chunk>`.
 Please refer to the :doc:`chunk HOWTO <Howto_chunk>` section for an overview.
 
-Reset_ids command
+Box command
------------------
+-----------
 
-The reset_ids command has been renamed to :doc:`reset_atom_ids <reset_atom_ids>`.
+.. deprecated:: TBD
+
+The *box* command has been removed and the LAMMPS code changed so it won't
+be needed.  If present, LAMMPS will ignore the command and print a warning.
+
+Reset_ids, reset_atom_ids, reset_mol_ids commands
+-------------------------------------------------
+
+.. deprecated:: TBD
+
+The *reset_ids*, *reset_atom_ids*, and *reset_mol_ids* commands have
+been folded into the :doc:`reset_atoms <reset_atoms>` command.  If
+present, LAMMPS will replace the commands accordingly and print a
+warning.
 
 MEAM package
 ------------
@@ -30,18 +46,21 @@ The code in the :ref:`MEAM package <PKG-MEAM>` is a translation of the
 Fortran code of MEAM into C++, which removes several restrictions
 (e.g. there can be multiple instances in hybrid pair styles) and allows
 for some optimizations leading to better performance.  The pair style
-:doc:`meam <pair_meam>` has the exact same syntax.
+:doc:`meam <pair_meam>` has the exact same syntax.  For a transition
+period the C++ version of MEAM was called USER-MEAMC so it could
+coexist with the Fortran version.
 
 REAX package
 ------------
 
 The REAX package has been removed since it was superseded by the
-:ref:`REAXFF package <PKG-REAXFF>`.  The REAXFF
+:ref:`REAXFF package <PKG-REAXFF>`.  The REAXFF package has been tested
-package has been tested to yield equivalent results to the REAX package,
+to yield equivalent results to the REAX package, offers better
-offers better performance, supports OpenMP multi-threading via OPENMP,
+performance, supports OpenMP multi-threading via OPENMP, and GPU and
-and GPU and threading parallelization through KOKKOS.  The new pair styles
+threading parallelization through KOKKOS.  The new pair styles are not
-are not syntax compatible with the removed reax pair style, so input
+syntax compatible with the removed reax pair style, so input files will
-files will have to be adapted.
+have to be adapted.  The REAXFF package was originally called
+USER-REAXC.
 
 USER-CUDA package
 -----------------
@@ -60,5 +79,6 @@ restart2data tool
 The functionality of the restart2data tool has been folded into the
 LAMMPS executable directly instead of having a separate tool.  A
 combination of the commands :doc:`read_restart <read_restart>` and
-:doc:`write_data <write_data>` can be used to the same effect.  For added
+:doc:`write_data <write_data>` can be used to the same effect.  For
-convenience this conversion can also be triggered by :doc:`command line flags <Run_options>`
+added convenience this conversion can also be triggered by
+:doc:`command line flags <Run_options>`

doc/src/Developer_code_design.rst

@@ -50,7 +50,7 @@ parallel each MPI process creates such an instance.  This can be seen
 in the ``main.cpp`` file where the core steps of running a LAMMPS
 simulation are the following 3 lines of code:
 
-.. code-block:: C++
+.. code-block:: c++
 
     LAMMPS *lammps = new LAMMPS(argc, argv, lammps_comm);
     lammps->input->file();
@@ -78,7 +78,7 @@ LAMMPS makes extensive use of the object oriented programming (OOP)
 principles of *compositing* and *inheritance*. Classes like the
 ``LAMMPS`` class are a **composite** containing pointers to instances
 of other classes like ``Atom``, ``Comm``, ``Force``, ``Neighbor``,
-``Modify``, and so on.  Each of these classes implement certain
+``Modify``, and so on.  Each of these classes implements certain
 functionality by storing and manipulating data related to the
 simulation and providing member functions that trigger certain
 actions.  Some of those classes like ``Force`` are themselves
@@ -87,9 +87,9 @@ interactions.  Similarly the ``Modify`` class contains a list of
 ``Fix`` and ``Compute`` classes.  If the input commands that
 correspond to these classes include the word *style*, then LAMMPS
 stores only a single instance of that class.  E.g. *atom_style*,
-*comm_style*, *pair_style*, *bond_style*.  It the input command does
+*comm_style*, *pair_style*, *bond_style*.  If the input command does
-not include the word *style*, there can be many instances of that
+**not** include the word *style*, then there may be many instances of
-class defined.  E.g. *region*, *fix*, *compute*, *dump*.
+that class defined, for example *region*, *fix*, *compute*, *dump*.
 
 **Inheritance** enables creation of *derived* classes that can share
 common functionality in their base class while providing a consistent
@@ -232,7 +232,7 @@ macro ``PairStyle()`` will associate the style name "lj/cut"
 with a factory function creating an instance of the ``PairLJCut``
 class.
 
-.. code-block:: C++
+.. code-block:: c++
 
    // from force.h
    typedef Pair *(*PairCreator)(LAMMPS *);
@@ -360,7 +360,7 @@ characters; "{:<8}" would do this as left aligned, "{:^8}" as centered,
 argument type must be compatible or else the {fmt} formatting code will
 throw an exception. Some format string examples are given below:
 
-.. code-block:: C
+.. code-block:: c++
 
    auto mesg = fmt::format("  CPU time: {:4d}:{:02d}:{:02d}\n", cpuh, cpum, cpus);
    mesg = fmt::format("{:<8s}| {:<10.5g} | {:<10.5g} | {:<10.5g} |{:6.1f} |{:6.2f}\n",

doc/src/Developer_notes.rst

@@ -105,7 +105,7 @@ list, where each pair of atoms is listed only once (except when the
 pairs straddling sub-domains or periodic boundaries will be listed twice).
 Thus these are the default settings when a neighbor list request is created in:
 
-.. code-block:: C++
+.. code-block:: c++
 
    void Pair::init_style()
    {
@@ -129,7 +129,7 @@ neighbor list request to the specific needs of a style an additional
 request flag is needed.  The :doc:`tersoff <pair_tersoff>` pair style,
 for example, needs a "full" neighbor list:
 
-.. code-block:: C++
+.. code-block:: c++
 
    void PairTersoff::init_style()
    {
@@ -141,7 +141,7 @@ When a pair style supports r-RESPA time integration with different cutoff region
 the request flag may depend on the corresponding r-RESPA settings. Here an example
 from pair style lj/cut:
 
-.. code-block:: C++
+.. code-block:: c++
 
    void PairLJCut::init_style()
    {
@@ -160,7 +160,7 @@ Granular pair styles need neighbor lists based on particle sizes and not cutoff
 and also may require to have the list of previous neighbors available ("history").
 For example with:
 
-.. code-block:: C++
+.. code-block:: c++
 
    if (use_history) neighbor->add_request(this, NeighConst::REQ_SIZE | NeighConst::REQ_HISTORY);
    else neighbor->add_request(this, NeighConst::REQ_SIZE);
@@ -170,7 +170,7 @@ settings each request can set an id which is then used in the corresponding
 ``init_list()`` function to assign it to the suitable pointer variable. This is
 done for example by the :doc:`pair style meam <pair_meam>`:
 
-.. code-block:: C++
+.. code-block:: c++
 
    void PairMEAM::init_style()
    {
@@ -189,7 +189,7 @@ just once) and this can also be indicated by a flag.  As an example here
 is the request from the ``FixPeriNeigh`` class which is created
 internally by :doc:`Peridynamics pair styles <pair_peri>`:
 
-.. code-block:: C++
+.. code-block:: c++
 
    neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL);
 
@@ -198,7 +198,7 @@ than what is usually inferred from the pair style settings (largest cutoff of
 all pair styles plus neighbor list skin).  The following is used in the
 :doc:`compute rdf <compute_rdf>` command implementation:
 
-.. code-block:: C++
+.. code-block:: c++
 
   if (cutflag)
     neighbor->add_request(this, NeighConst::REQ_OCCASIONAL)->set_cutoff(mycutneigh);
@@ -212,7 +212,7 @@ for printing the neighbor list summary the name of the requesting command
 should be set.  Below is the request from the :doc:`delete atoms <delete_atoms>`
 command:
 
-.. code-block:: C++
+.. code-block:: c++
 
    neighbor->add_request(this, "delete_atoms", NeighConst::REQ_FULL);
 

doc/src/Developer_plugins.rst

@@ -95,7 +95,7 @@ a class ``PairMorse2`` in the files ``pair_morse2.h`` and
 ``pair_morse2.cpp`` with the factory function and initialization
 function would look like this:
 
-.. code-block:: C++
+.. code-block:: c++
 
   #include "lammpsplugin.h"
   #include "version.h"
@@ -141,7 +141,7 @@ list of argument strings), then the pointer type is ``lammpsplugin_factory2``
 and it must be assigned to the *creator.v2* member of the plugin struct.
 Below is an example for that:
 
-.. code-block:: C++
+.. code-block:: c++
 
   #include "lammpsplugin.h"
   #include "version.h"
@@ -176,7 +176,7 @@ demonstrated in the following example, which also shows that the
 implementation of the plugin class may be within the same source
 file as the plugin interface code:
 
-.. code-block:: C++
+.. code-block:: c++
 
    #include "lammpsplugin.h"
 

doc/src/Developer_unittest.rst

@@ -194,7 +194,7 @@ macro. These tests operate by capturing the screen output when executing
 the failing command and then comparing that with a provided regular
 expression string pattern.  Example:
 
-.. code-block:: C++
+.. code-block:: c++
 
    TEST_F(SimpleCommandsTest, UnknownCommand)
    {
@@ -226,9 +226,9 @@ The following test programs are currently available:
    * - ``test_kim_commands.cpp``
      - KimCommands
      - Tests for several commands from the :ref:`KIM package <PKG-KIM>`
-   * - ``test_reset_ids.cpp``
+   * - ``test_reset_atoms.cpp``
-     - ResetIDs
+     - ResetAtoms
-     - Tests to validate the :doc:`reset_atom_ids <reset_atom_ids>` and :doc:`reset_mol_ids <reset_mol_ids>` commands
+     - Tests to validate the :doc:`reset_atoms <reset_atoms>` sub-commands
 
 
 Tests for the C-style library interface
@@ -249,7 +249,7 @@ MPI support.  These include tests where LAMMPS is run in multi-partition
 mode or only on a subset of the MPI world communicator.  The CMake
 script code for adding this kind of test looks like this:
 
-.. code-block:: CMake
+.. code-block:: cmake
 
    if (BUILD_MPI)
      add_executable(test_library_mpi test_library_mpi.cpp)

doc/src/Developer_updating.rst

@@ -61,7 +61,7 @@ header file needs to be updated accordingly.
 
 Old:
 
-.. code-block:: C++
+.. code-block:: c++
 
    int PairEAM::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
    {
@@ -75,7 +75,7 @@ Old:
 
 New:
 
-.. code-block:: C++
+.. code-block:: c++
 
    int PairEAM::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
    {
@@ -112,14 +112,14 @@ Example from a pair style:
 
 Old:
 
-.. code-block:: C++
+.. code-block:: c++
 
    if (eflag || vflag) ev_setup(eflag, vflag);
    else evflag = vflag_fdotr = eflag_global = eflag_atom = 0;
 
 New:
 
-.. code-block:: C++
+.. code-block:: c++
 
    ev_init(eflag, vflag);
 
@@ -142,14 +142,14 @@ when they are called from only one or a subset of the MPI processes.
 
 Old:
 
-.. code-block:: C++
+.. code-block:: c++
 
     val = force->numeric(FLERR, arg[1]);
     num = force->inumeric(FLERR, arg[2]);
 
 New:
 
-.. code-block:: C++
+.. code-block:: c++
 
     val = utils::numeric(FLERR, true, arg[1], lmp);
     num = utils::inumeric(FLERR, false, arg[2], lmp);
@@ -183,14 +183,14 @@ copy them around for simulations.
 
 Old:
 
-.. code-block:: C++
+.. code-block:: c++
 
    fp = force->open_potential(filename);
    fp = fopen(filename, "r");
 
 New:
 
-.. code-block:: C++
+.. code-block:: c++
 
    fp = utils::open_potential(filename, lmp);
 
@@ -207,7 +207,7 @@ Example:
 
 Old:
 
-.. code-block:: C++
+.. code-block:: c++
 
    if (fptr == NULL) {
      char str[128];
@@ -217,7 +217,7 @@ Old:
 
 New:
 
-.. code-block:: C++
+.. code-block:: c++
 
    if (fptr == nullptr)
      error->one(FLERR, "Cannot open AEAM potential file {}: {}", filename, utils::getsyserror());
@@ -237,7 +237,7 @@ an example from the ``FixWallReflect`` class:
 
 Old:
 
-.. code-block:: C++
+.. code-block:: c++
 
    FixWallReflect(class LAMMPS *, int, char **);
    virtual ~FixWallReflect();
@@ -247,7 +247,7 @@ Old:
 
 New:
 
-.. code-block:: C++
+.. code-block:: c++
 
    FixWallReflect(class LAMMPS *, int, char **);
    ~FixWallReflect() override;
@@ -271,7 +271,7 @@ the type of the "this" pointer argument.
 
 Old:
 
-.. code-block:: C++
+.. code-block:: c++
 
    comm->forward_comm_pair(this);
    comm->forward_comm_fix(this);
@@ -284,7 +284,7 @@ Old:
 
 New:
 
-.. code-block:: C++
+.. code-block:: c++
 
    comm->forward_comm(this);
    comm->reverse_comm(this);
@@ -304,7 +304,7 @@ requests can be :doc:`found here <Developer_notes>`. Example from the
 
 Old:
 
-.. code-block:: C++
+.. code-block:: c++
 
    int irequest = neighbor->request(this,instance_me);
    neighbor->requests[irequest]->pair = 0;
@@ -317,7 +317,7 @@ Old:
 
 New:
 
-.. code-block:: C++
+.. code-block:: c++
 
    auto req = neighbor->add_request(this, NeighConst::REQ_OCCASIONAL);
    if (cutflag) req->set_cutoff(mycutneigh);
@@ -340,7 +340,7 @@ these are internal fixes, there is no user visible change.
 
 Old:
 
-.. code-block:: C++
+.. code-block:: c++
 
    #include "fix_store.h"
 
@@ -351,7 +351,7 @@ Old:
 
 New:
 
-.. code-block:: C++
+.. code-block:: c++
 
    #include "fix_store_peratom.h"
 
@@ -362,7 +362,7 @@ New:
 
 Old:
 
-.. code-block:: C++
+.. code-block:: c++
 
    #include "fix_store.h"
 
@@ -373,7 +373,7 @@ Old:
 
 New:
 
-.. code-block:: C++
+.. code-block:: c++
 
    #include "fix_store_global.h"
 
@@ -396,7 +396,7 @@ the dump directly.  Example:
 
 Old:
 
-.. code-block:: C++
+.. code-block:: c++
 
    int idump = output->find_dump(arg[iarg+1]);
    if (idump < 0)
@@ -412,7 +412,7 @@ Old:
 
 New:
 
-.. code-block:: C++
+.. code-block:: c++
 
    auto idump = output->get_dump_by_id(arg[iarg+1]);
    if (!idump) error->all(FLERR,"Dump ID {} in hyper command does not exist", arg[iarg+1]);

doc/src/Developer_utils.rst

@@ -133,6 +133,9 @@ and parsing files or arguments.
 .. doxygenfunction:: trim_comment
    :project: progguide
 
+.. doxygenfunction:: strip_style_suffix
+   :project: progguide
+
 .. doxygenfunction:: star_subst
    :project: progguide
 
@@ -317,7 +320,7 @@ are all "whitespace" characters, i.e. the space character, the tabulator
 character, the carriage return character, the linefeed character, and
 the form feed character.
 
-.. code-block:: C++
+.. code-block:: c++
    :caption: Tokenizer class example listing entries of the PATH environment variable
 
    #include "tokenizer.h"
@@ -349,7 +352,7 @@ tokenizer into a ``try`` / ``catch`` block to handle errors.  The
 when a (type of) number is requested as next token that is not
 compatible with the string representing the next word.
 
-.. code-block:: C++
+.. code-block:: c++
    :caption: ValueTokenizer class example with exception handling
 
    #include "tokenizer.h"
@@ -427,7 +430,7 @@ one or two array indices "[<number>]" with numbers > 0.
 
 A typical code segment would look like this:
 
-.. code-block:: C++
+.. code-block:: c++
    :caption: Usage example for ArgInfo class
 
    int nvalues = 0;
@@ -476,7 +479,7 @@ open the file, and will call the :cpp:class:`LAMMPS_NS::Error` class in
 case of failures to read or to convert numbers, so that LAMMPS will be
 aborted.
 
-.. code-block:: C++
+.. code-block:: c++
    :caption: Use of PotentialFileReader class in pair style coul/streitz
 
     PotentialFileReader reader(lmp, file, "coul/streitz");
@@ -555,7 +558,7 @@ chunk size needs to be known in advance, 2) with :cpp:func:`MyPage::vget()
 its size is registered later with :cpp:func:`MyPage::vgot()
 <LAMMPS_NS::MyPage::vgot>`.
 
-.. code-block:: C++
+.. code-block:: c++
    :caption: Example of using :cpp:class:`MyPage <LAMMPS_NS::MyPage>`
 
       #include "my_page.h"

doc/src/Developer_write.rst

@@ -26,7 +26,7 @@ constructor with the signature: ``FixPrintVel(class LAMMPS *, int, char **)``.
 Every fix must be registered in LAMMPS by writing the following lines
 of code in the header before include guards:
 
-.. code-block:: c
+.. code-block:: c++
 
    #ifdef FIX_CLASS
    // clang-format off
@@ -47,7 +47,7 @@ keyword when it parses the input script.
 Let's write a simple fix which will print the average velocity at the end
 of each timestep. First of all, implement a constructor:
 
-.. code-block:: C++
+.. code-block:: c++
 
    FixPrintVel::FixPrintVel(LAMMPS *lmp, int narg, char **arg)
    : Fix(lmp, narg, arg)
@@ -72,7 +72,7 @@ in the Fix class called ``nevery`` which specifies how often the method
 
 The next method we need to implement is ``setmask()``:
 
-.. code-block:: C++
+.. code-block:: c++
 
    int FixPrintVel::setmask()
    {
@@ -87,7 +87,7 @@ during execution. The constant ``END_OF_STEP`` corresponds to the
 are called during a timestep and the order in which they are called
 are shown in the previous section.
 
-.. code-block:: C++
+.. code-block:: c++
 
    void FixPrintVel::end_of_step()
    {
@@ -143,7 +143,7 @@ The group membership information of an atom is contained in the *mask*
 property of and atom and the bit corresponding to a given group is
 stored in the groupbit variable which is defined in Fix base class:
 
-.. code-block:: C++
+.. code-block:: c++
 
    for (int i = 0; i < nlocal; ++i) {
      if (atom->mask[i] & groupbit) {
@@ -174,7 +174,7 @@ to store positions of atoms from previous timestep, we need to add
 ``double** xold`` to the header file. Than add allocation code
 to the constructor:
 
-.. code-block:: C++
+.. code-block:: c++
 
    FixSavePos::FixSavePos(LAMMPS *lmp, int narg, char **arg), xold(nullptr)
    {
@@ -190,7 +190,7 @@ to the constructor:
 
 Implement the aforementioned methods:
 
-.. code-block:: C++
+.. code-block:: c++
 
    double FixSavePos::memory_usage()
    {

doc/src/Howto_pylammps.rst

@@ -152,14 +152,14 @@ Creating a new instance of PyLammps
 To create a PyLammps object you need to first import the class from the lammps
 module. By using the default constructor, a new *lammps* instance is created.
 
-.. code-block:: Python
+.. code-block:: python
 
    from lammps import PyLammps
    L = PyLammps()
 
 You can also initialize PyLammps on top of this existing *lammps* object:
 
-.. code-block:: Python
+.. code-block:: python
 
    from lammps import lammps, PyLammps
    lmp = lammps()
@@ -180,14 +180,14 @@ For instance, let's take the following LAMMPS command:
 In the original interface this command can be executed with the following
 Python code if *L* was a lammps instance:
 
-.. code-block:: Python
+.. code-block:: python
 
    L.command("region box block 0 10 0 5 -0.5 0.5")
 
 With the PyLammps interface, any command can be split up into arbitrary parts
 separated by white-space, passed as individual arguments to a region method.
 
-.. code-block:: Python
+.. code-block:: python
 
    L.region("box block", 0, 10, 0, 5, -0.5, 0.5)
 
@@ -199,14 +199,14 @@ The benefit of this approach is avoiding redundant command calls and easier
 parameterization. In the original interface parameterization needed to be done
 manually by creating formatted strings.
 
-.. code-block:: Python
+.. code-block:: python
 
    L.command("region box block %f %f %f %f %f %f" % (xlo, xhi, ylo, yhi, zlo, zhi))
 
 In contrast, methods of PyLammps accept parameters directly and will convert
 them automatically to a final command string.
 
-.. code-block:: Python
+.. code-block:: python
 
    L.region("box block", xlo, xhi, ylo, yhi, zlo, zhi)
 
@@ -256,7 +256,7 @@ LAMMPS variables can be both defined and accessed via the PyLammps interface.
 
 To define a variable you can use the :doc:`variable <variable>` command:
 
-.. code-block:: Python
+.. code-block:: python
 
    L.variable("a index 2")
 
@@ -265,14 +265,14 @@ A dictionary of all variables is returned by L.variables
 you can access an individual variable by retrieving a variable object from the
 L.variables dictionary by name
 
-.. code-block:: Python
+.. code-block:: python
 
    a = L.variables['a']
 
 The variable value can then be easily read and written by accessing the value
 property of this object.
 
-.. code-block:: Python
+.. code-block:: python
 
    print(a.value)
    a.value = 4
@@ -284,7 +284,7 @@ LAMMPS expressions can be immediately evaluated by using the eval method. The
 passed string parameter can be any expression containing global thermo values,
 variables, compute or fix data.
 
-.. code-block:: Python
+.. code-block:: python
 
    result = L.eval("ke") # kinetic energy
    result = L.eval("pe") # potential energy
@@ -298,7 +298,7 @@ All atoms in the current simulation can be accessed by using the L.atoms list.
 Each element of this list is an object which exposes its properties (id, type,
 position, velocity, force, etc.).
 
-.. code-block:: Python
+.. code-block:: python
 
    # access first atom
    L.atoms[0].id
@@ -311,7 +311,7 @@ position, velocity, force, etc.).
 
 Some properties can also be used to set:
 
-.. code-block:: Python
+.. code-block:: python
 
    # set position in 2D simulation
    L.atoms[0].position = (1.0, 0.0)
@@ -328,7 +328,7 @@ after a run via the L.runs list. This list contains a growing list of run data.
 The first element is the output of the first run, the second element that of
 the second run.
 
-.. code-block:: Python
+.. code-block:: python
 
    L.run(1000)
    L.runs[0] # data of first 1000 time steps
@@ -339,14 +339,14 @@ the second run.
 Each run contains a dictionary of all trajectories. Each trajectory is
 accessible through its thermo name:
 
-.. code-block:: Python
+.. code-block:: python
 
    L.runs[0].thermo.Step # list of time steps in first run
    L.runs[0].thermo.Ke   # list of kinetic energy values in first run
 
 Together with matplotlib plotting data out of LAMMPS becomes simple:
 
-.. code-block:: Python
+.. code-block:: python
 
    import matplotlib.plot as plt
    steps = L.runs[0].thermo.Step
@@ -406,7 +406,7 @@ Four atoms are placed in the simulation and the dihedral potential is applied on
 them using a datafile. Then one of the atoms is rotated along the central axis by
 setting its position from Python, which changes the dihedral angle.
 
-.. code-block:: Python
+.. code-block:: python
 
    phi = [d \* math.pi / 180 for d in range(360)]
 
@@ -439,7 +439,7 @@ Initially, a 2D system is created in a state with minimal energy.
 
 It is then disordered by moving each atom by a random delta.
 
-.. code-block:: Python
+.. code-block:: python
 
    random.seed(27848)
    deltaperturb = 0.2
@@ -458,7 +458,7 @@ It is then disordered by moving each atom by a random delta.
 Finally, the Monte Carlo algorithm is implemented in Python. It continuously
 moves random atoms by a random delta and only accepts certain moves.
 
-.. code-block:: Python
+.. code-block:: python
 
    estart = L.eval("pe")
    elast = estart
@@ -517,7 +517,7 @@ PyLammps can be run in parallel using mpi4py. This python package can be install
 The following is a short example which reads in an existing LAMMPS input file and
 executes it in parallel.  You can find in.melt in the examples/melt folder.
 
-.. code-block:: Python
+.. code-block:: python
 
    from mpi4py import MPI
    from lammps import PyLammps

doc/src/Howto_structured_data.rst

@@ -43,7 +43,7 @@ JSON
      "ke": $(ke)
    }""" file current_state.json screen no
 
-.. code-block:: JSON
+.. code-block:: json
    :caption: current_state.json
 
    {

doc/src/Howto_triclinic.rst

@@ -144,11 +144,6 @@ does not change the atom positions due to non-periodicity.  In this
 mode, if you tilt the system to extreme angles, the simulation will
 simply become inefficient, due to the highly skewed simulation box.
 
-The limitation on not creating a simulation box with a tilt factor
-skewing the box more than half the distance of the parallel box length
-can be overridden via the :doc:`box <box>` command.  Setting the *tilt*
-keyword to *large* allows any tilt factors to be specified.
-
 Box flips that may occur using the :doc:`fix deform <fix_deform>` or
 :doc:`fix npt <fix_nh>` commands can be turned off using the *flip no*
 option with either of the commands.

doc/src/Library_config.rst

@@ -39,7 +39,7 @@ crashes within LAMMPS may be recovered from by enabling
 :ref:`exceptions <exceptions>`, avoiding them proactively is a safer
 approach.
 
-.. code-block:: C
+.. code-block:: c
    :caption: Example for using configuration settings functions
 
    #include "library.h"

doc/src/Library_create.rst

@@ -22,7 +22,7 @@ as the "handle" argument in subsequent function calls until that
 instance is destroyed by calling :cpp:func:`lammps_close`.  Here is a
 simple example demonstrating its use:
 
-.. code-block:: C
+.. code-block:: c
 
    #include "library.h"
    #include <stdio.h>

doc/src/Library_execute.rst

@@ -30,7 +30,7 @@ be included in the file or strings, and expansion of variables with
 ``${name}`` or ``$(expression)`` syntax is performed.
 Below is a short example using some of these functions.
 
-.. code-block:: C
+.. code-block:: c
 
    /* define to make the otherwise hidden prototype for "lammps_open()" visible */
    #define LAMMPS_LIB_MPI

doc/src/Library_properties.rst

@@ -32,7 +32,7 @@ indexed accordingly. Per-atom data can change sizes and ordering at
 every neighbor list rebuild or atom sort event as atoms migrate between
 sub-domains and processors.
 
-.. code-block:: C
+.. code-block:: c
 
    #include "library.h"
    #include <stdio.h>

doc/src/Modify_bond.rst

@@ -13,24 +13,65 @@ Here is a brief description of common methods you define in your
 new derived class.  See bond.h, angle.h, dihedral.h, and improper.h
 for details and specific additional methods.
 
-+-----------------------+---------------------------------------------------------------------------------------+
++-----------------------+---------------------------------------------------------------------+
-| init                  | check if all coefficients are set, calls *init_style* (optional)                      |
+| Required              | "pure" methods that *must* be overridden in a derived class         |
-+-----------------------+---------------------------------------------------------------------------------------+
++=======================+=====================================================================+
-| init_style            | check if style specific conditions are met (optional)                                 |
+| compute               | compute the molecular interactions for all listed items             |
-+-----------------------+---------------------------------------------------------------------------------------+
++-----------------------+---------------------------------------------------------------------+
-| compute               | compute the molecular interactions (required)                                         |
+| coeff                 | set coefficients for one type                                       |
-+-----------------------+---------------------------------------------------------------------------------------+
++-----------------------+---------------------------------------------------------------------+
-| settings              | apply global settings for all types (optional)                                        |
+| equilibrium_distance  | length of bond, used by SHAKE (bond styles only)                    |
-+-----------------------+---------------------------------------------------------------------------------------+
++-----------------------+---------------------------------------------------------------------+
-| coeff                 | set coefficients for one type (required)                                              |
+| equilibrium_angle     | opening of angle, used by SHAKE (angle styles only)                 |
-+-----------------------+---------------------------------------------------------------------------------------+
++-----------------------+---------------------------------------------------------------------+
-| equilibrium_distance  | length of bond, used by SHAKE (required, bond only)                                   |
+| write & read_restart  | writes/reads coeffs to restart files                                |
-+-----------------------+---------------------------------------------------------------------------------------+
++-----------------------+---------------------------------------------------------------------+
-| equilibrium_angle     | opening of angle, used by SHAKE (required, angle only)                                |
+| single                | force/r (bond styles only) and energy of a single bond or angle     |
-+-----------------------+---------------------------------------------------------------------------------------+
++-----------------------+---------------------------------------------------------------------+
-| write & read_restart  | writes/reads coeffs to restart files (required)                                       |
+
-+-----------------------+---------------------------------------------------------------------------------------+
+
-| single                | force (bond only) and energy of a single bond or angle (required, bond or angle only) |
++--------------------------------+----------------------------------------------------------------------+
-+-----------------------+---------------------------------------------------------------------------------------+
+| Optional                       | methods that have a default or dummy implementation                  |
-| memory_usage          | tally memory allocated by the style (optional)                                        |
++================================+======================================================================+
-+-----------------------+---------------------------------------------------------------------------------------+
+| init                           | check if all coefficients are set, calls init_style()                |
++--------------------------------+----------------------------------------------------------------------+
+| init_style                     | check if style specific conditions are met                           |
++--------------------------------+----------------------------------------------------------------------+
+| settings                       | apply global settings for all types                                  |
++--------------------------------+----------------------------------------------------------------------+
+| write & read_restart_settings  | writes/reads global style settings to restart files                  |
++--------------------------------+----------------------------------------------------------------------+
+| write_data                     | write corresponding Coeffs section(s) in data file                   |
++--------------------------------+----------------------------------------------------------------------+
+| memory_usage                   | tally memory allocated by the style                                  |
++--------------------------------+----------------------------------------------------------------------+
+| extract                        | provide access to internal data  (bond or angle styles only)         |
++--------------------------------+----------------------------------------------------------------------+
+| reinit                         | reset all type-based parameters, called by fix adapt (bonds only)    |
++--------------------------------+----------------------------------------------------------------------+
+| pack & unpack_forward_comm     | copy data to and from buffer in forward communication (bonds only)   |
++--------------------------------+----------------------------------------------------------------------+
+| pack & unpack_reverse_comm     | copy data to and from buffer in reverse communication (bonds only)   |
++--------------------------------+----------------------------------------------------------------------+
+
+Here is a list of flags or settings that should be set in the
+constructor of the derived class when they differ from the default
+setting.
+
++---------------------------------+------------------------------------------------------------------------------+---------+
+| Name of flag                    | Description                                                                  | default |
++=================================+==============================================================================+=========+
+| writedata                       | 1 if write_data() is implemented                                             | 1       |
++---------------------------------+------------------------------------------------------------------------------+---------+
+| single_extra                    | number of extra single values calculated (bond styles only)                  | 0       |
++---------------------------------+------------------------------------------------------------------------------+---------+
+| partial_flag                    | 1 if bond type can be set to 0 and deleted (bond styles only)                | 0       |
++---------------------------------+------------------------------------------------------------------------------+---------+
+| reinitflag                      | 1 if style has reinit() and is compatible with fix adapt                     | 1       |
++---------------------------------+------------------------------------------------------------------------------+---------+
+| comm_forward                    | size of buffer (in doubles) for forward communication  (bond styles only)    | 0       |
++---------------------------------+------------------------------------------------------------------------------+---------+
+| comm_reverse                    | size of buffer (in doubles) for reverse communication  (bond styles only)    | 0       |
++---------------------------------+------------------------------------------------------------------------------+---------+
+| comm_reverse_off                | size of buffer for reverse communication with newton off (bond styles only)  | 0       |
++---------------------------------+------------------------------------------------------------------------------+---------+

doc/src/Modify_pair.rst

@@ -1,35 +1,121 @@
 Pair styles
 ===========
 
-Classes that compute pairwise interactions are derived from the Pair
+Classes that compute pairwise non-bonded interactions are derived from
-class.  In LAMMPS, pairwise calculation include many-body potentials
+the Pair class.  In LAMMPS, pairwise calculation include many-body
-such as EAM or Tersoff where particles interact without a static bond
+potentials such as EAM, Tersoff, or ReaxFF where particles interact
-topology.  New styles can be created to add new pair potentials to
+without an explicit bond topology but include interactions beyond
-LAMMPS.
+pairwise non-bonded contributions.  New styles can be created to add
+support for additional pair potentials to LAMMPS.  When the
+modifications are small, sometimes it is more effective to derive from
+an existing pair style class.  This latter approach is also used by
+:doc:`Accelerator packages <Speed_packages>` where the accelerated style
+names differ from their base classes by an appended suffix.
 
-Pair_lj_cut.cpp is a simple example of a Pair class, though it
+The file ``src/pair_lj_cut.cpp`` is an example of a Pair class with a
-includes some optional methods to enable its use with rRESPA.
+very simple potential function.  It includes several optional methods to
+enable its use with :doc:`run_style respa <run_style>` and :doc:`compute
+group/group <compute_group_group>`.
 
-Here is a brief description of the class methods in pair.h:
+Here is a brief list of some the class methods in the Pair class that
+*must* be or *may* be overridden in a derived class.
 
 +---------------------------------+---------------------------------------------------------------------+
+| Required                        | "pure" methods that *must* be overridden in a derived class         |
++=================================+=====================================================================+
 | compute                         | workhorse routine that computes pairwise interactions               |
 +---------------------------------+---------------------------------------------------------------------+
-| settings                        | reads the input script line with arguments you define               |
+| settings                        | processes the arguments to the pair_style command                   |
 +---------------------------------+---------------------------------------------------------------------+
-| coeff                           | set coefficients for one i,j type pair                              |
+| coeff                           | set coefficients for one i,j type pair, called from pair_coeff      |
-+---------------------------------+---------------------------------------------------------------------+
-| init_one                        | perform initialization for one i,j type pair                        |
-+---------------------------------+---------------------------------------------------------------------+
-| init_style                      | initialization specific to this pair style                          |
-+---------------------------------+---------------------------------------------------------------------+
-| write & read_restart            | write/read i,j pair coeffs to restart files                         |
-+---------------------------------+---------------------------------------------------------------------+
-| write & read_restart_settings   | write/read global settings to restart files                         |
-+---------------------------------+---------------------------------------------------------------------+
-| single                          | force/r and energy of a single pairwise interaction between 2 atoms |
-+---------------------------------+---------------------------------------------------------------------+
-| compute_inner/middle/outer      | versions of compute used by rRESPA                                  |
 +---------------------------------+---------------------------------------------------------------------+
 
-The inner/middle/outer routines are optional.
++---------------------------------+----------------------------------------------------------------------+
+| Optional                        | methods that have a default or dummy implementation                  |
++=================================+======================================================================+
+| init_one                        | perform initialization for one i,j type pair                         |
++---------------------------------+----------------------------------------------------------------------+
+| init_style                      | style initialization: request neighbor list(s), error checks         |
++---------------------------------+----------------------------------------------------------------------+
+| init_list                       | Neighbor class callback function to pass neighbor list to pair style |
++---------------------------------+----------------------------------------------------------------------+
+| single                          | force/r and energy of a single pairwise interaction between 2 atoms  |
++---------------------------------+----------------------------------------------------------------------+
+| compute_inner/middle/outer      | versions of compute used by rRESPA                                   |
++---------------------------------+----------------------------------------------------------------------+
+| memory_usage                    | return estimated amount of memory used by the pair style             |
++---------------------------------+----------------------------------------------------------------------+
+| modify_params                   | process arguments to pair_modify command                             |
++---------------------------------+----------------------------------------------------------------------+
+| extract                         | provide access to internal scalar or per-type data like cutoffs      |
++---------------------------------+----------------------------------------------------------------------+
+| extract_peratom                 | provide access to internal per-atom data                             |
++---------------------------------+----------------------------------------------------------------------+
+| setup                           | initialization at the beginning of a run                             |
++---------------------------------+----------------------------------------------------------------------+
+| finish                          | called at the end of a run, e.g. to print                            |
++---------------------------------+----------------------------------------------------------------------+
+| write & read_restart            | write/read i,j pair coeffs to restart files                          |
++---------------------------------+----------------------------------------------------------------------+
+| write & read_restart_settings   | write/read global settings to restart files                          |
++---------------------------------+----------------------------------------------------------------------+
+| write_data                      | write Pair Coeffs section to data file                               |
++---------------------------------+----------------------------------------------------------------------+
+| write_data_all                  | write PairIJ Coeffs section to data file                             |
++---------------------------------+----------------------------------------------------------------------+
+| pack & unpack_forward_comm      | copy data to and from buffer if style uses forward communication     |
++---------------------------------+----------------------------------------------------------------------+
+| pack & unpack_reverse_comm      | copy data to and from buffer if style uses reverse communication     |
++---------------------------------+----------------------------------------------------------------------+
+| reinit                          | reset all type-based parameters, called by fix adapt for example     |
++---------------------------------+----------------------------------------------------------------------+
+| reset_dt                        | called when the time step is changed by timestep or fix reset/dt     |
++---------------------------------+----------------------------------------------------------------------+
+
+Here is a list of flags or settings that should be set in the
+constructor of the derived pair class when they differ from the default
+setting.
+
++---------------------------------+-------------------------------------------------------------+---------+
+| Name of flag                    | Description                                                 | default |
++=================================+=============================================================+=========+
+| single_enable                   | 1 if single() method is implemented, 0 if missing           | 1       |
++---------------------------------+-------------------------------------------------------------+---------+
+| respa_enable                    | 1 if pair style has compute_inner/middle/outer()            | 0       |
++---------------------------------+-------------------------------------------------------------+---------+
+| restartinfo                     | 1 if pair style writes its settings to a restart            | 1       |
++---------------------------------+-------------------------------------------------------------+---------+
+| one_coeff                       | 1 if only a pair_coeff * * command is allowed               | 0       |
++---------------------------------+-------------------------------------------------------------+---------+
+| manybody_flag                   | 1 if pair style is a manybody potential                     | 0       |
++---------------------------------+-------------------------------------------------------------+---------+
+| unit_convert_flag               | value != 0 indicates support for unit conversion            | 0       |
++---------------------------------+-------------------------------------------------------------+---------+
+| no_virial_fdotr_compute         | 1 if pair style does not call virial_fdotr_compute()        | 0       |
++---------------------------------+-------------------------------------------------------------+---------+
+| writedata                       | 1 if write_data() and write_data_all() are implemented      | 0       |
++---------------------------------+-------------------------------------------------------------+---------+
+| comm_forward                    | size of buffer (in doubles) for forward communication       | 0       |
++---------------------------------+-------------------------------------------------------------+---------+
+| comm_reverse                    | size of buffer (in doubles) for reverse communication       | 0       |
++---------------------------------+-------------------------------------------------------------+---------+
+| ghostneigh                      | 1 if cutghost is set and style uses neighbors of ghosts     | 0       |
++---------------------------------+-------------------------------------------------------------+---------+
+| finitecutflag                   | 1 if cutoff depends on diameter of atoms                    | 0       |
++---------------------------------+-------------------------------------------------------------+---------+
+| reinitflag                      | 1 if style has reinit() and is compatible with fix adapt    | 0       |
++---------------------------------+-------------------------------------------------------------+---------+
+| ewaldflag                       | 1 if compatible with kspace_style ewald                     | 0       |
++---------------------------------+-------------------------------------------------------------+---------+
+| pppmflag                        | 1 if compatible with kspace_style pppm                      | 0       |
++---------------------------------+-------------------------------------------------------------+---------+
+| msmflag                         | 1 if compatible with kspace_style msm                       | 0       |
++---------------------------------+-------------------------------------------------------------+---------+
+| dispersionflag                  | 1 if compatible with ewald/disp or pppm/disp                | 0       |
++---------------------------------+-------------------------------------------------------------+---------+
+| tip4pflag                       | 1 if compatible with kspace_style pppm/tip4p                | 0       |
++---------------------------------+-------------------------------------------------------------+---------+
+| dipoleflag                      | 1 if compatible with dipole kspace_style                    | 0       |
++---------------------------------+-------------------------------------------------------------+---------+
+| spinflag                        | 1 if compatible with spin kspace_style                      | 0       |
++---------------------------------+-------------------------------------------------------------+---------+

doc/src/Packages_details.rst

@@ -80,6 +80,7 @@ page gives those details.
    * :ref:`ML-HDNNP <PKG-ML-HDNNP>`
    * :ref:`ML-IAP <PKG-ML-IAP>`
    * :ref:`ML-PACE <PKG-ML-PACE>`
+   * :ref:`ML-POD <PKG-ML-POD>`
    * :ref:`ML-QUIP <PKG-ML-QUIP>`
    * :ref:`ML-RANN <PKG-ML-RANN>`
    * :ref:`ML-SNAP <PKG-ML-SNAP>`
@@ -865,7 +866,7 @@ ELECTRODE package
 The ELECTRODE package allows the user to enforce a constant potential method for
 groups of atoms that interact with the remaining atoms as electrolyte.
 
-**Authors:** The ELECTRODE library is written and maintained by Ludwig
+**Authors:** The ELECTRODE package is written and maintained by Ludwig
 Ahrens-Iwers (TUHH, Hamburg, Germany), Shern Tee (UQ, Brisbane, Australia) and
 Robert Meissner (TUHH, Hamburg, Germany).
 
@@ -878,7 +879,7 @@ This package has :ref:`specific installation instructions <electrode>` on the
 
 **Supporting info:**
 
-* :doc:`fix electrode/conp <fix_electrode_conp>`
+* :doc:`fix electrode/conp <fix_electrode>`
 
 ----------
 
@@ -1485,8 +1486,9 @@ MC package
 
 Several fixes and a pair style that have Monte Carlo (MC) or MC-like
 attributes.  These include fixes for creating, breaking, and swapping
-bonds, for performing atomic swaps, and performing grand-canonical MC
+bonds, for performing atomic swaps, and performing grand canonical
-(GCMC) or similar processes in conjunction with dynamics.
+MC (GCMC), semi-grand canonical MC (SGCMC), or similar processes in
+conjunction with molecular dynamics (MD).
 
 **Supporting info:**
 
@@ -1498,6 +1500,7 @@ bonds, for performing atomic swaps, and performing grand-canonical MC
 * :doc:`fix bond/swap <fix_bond_swap>`
 * :doc:`fix charge/regulation <fix_charge_regulation>`
 * :doc:`fix gcmc <fix_gcmc>`
+* :doc:`fix sgcmc <fix_sgcmc>`
 * :doc:`fix tfmc <fix_tfmc>`
 * :doc:`fix widom <fix_widom>`
 * :doc:`pair_style dsmc <pair_dsmc>`
@@ -1796,6 +1799,39 @@ This package has :ref:`specific installation instructions <ml-pace>` on the
 
 ----------
 
+.. _PKG-ML-POD:
+
+ML-POD package
+-------------------
+
+**Contents:**
+
+A pair style and fitpod style for Proper Orthogonal Descriptors
+(POD). POD is a methodology for deriving descriptors based on the proper
+orthogonal decomposition. The ML-POD package provides an efficient
+implementation for running simulations with POD potentials, along with
+fitting the potentials natively in LAMMPS.
+
+**Authors:**
+
+Ngoc Cuong Nguyen (MIT), Andrew Rohskopf (Sandia)
+
+.. versionadded:: TBD
+
+**Install:**
+
+This package has :ref:`specific installation instructions <ml-pod>` on the
+:doc:`Build extras <Build_extras>` page.
+
+**Supporting info:**
+
+* src/ML-POD: filenames -> commands
+* :doc:`pair_style pod <pair_pod>`
+* :doc:`command_style fitpod <fitpod_command>`
+* examples/PACKAGES/pod
+
+----------
+
 .. _PKG-ML-QUIP:
 
 ML-QUIP package

doc/src/Packages_list.rst

@@ -155,7 +155,7 @@ whether an extra library is needed to build and use the package:
      - no
    * - :ref:`ELECTRODE <PKG-ELECTRODE>`
      - electrode charges to match potential
-     - :doc:`fix electrode/conp <fix_electrode_conp>`
+     - :doc:`fix electrode/conp <fix_electrode>`
      - PACKAGES/electrode
      - no
    * - :ref:`EXTRA-COMPUTE <PKG-EXTRA-COMPUTE>`
@@ -298,6 +298,11 @@ whether an extra library is needed to build and use the package:
      - :doc:`pair pace <pair_pace>`
      - PACKAGES/pace
      - ext
+   * - :ref:`ML-POD <PKG-ML-POD>`
+     - Proper orthogonal decomposition potentials
+     - :doc:`pair pod <pair_pod>`
+     - pod
+     - ext
    * - :ref:`ML-QUIP <PKG-ML-QUIP>`
      - QUIP/libatoms interface
      - :doc:`pair_style quip <pair_quip>`

doc/src/Python_atoms.rst

@@ -58,7 +58,7 @@ against invalid accesses.
       Each element of this list is a :py:class:`Atom <lammps.Atom>` or :py:class:`Atom2D <lammps.Atom2D>` object. The attributes of
       these objects provide access to their data (id, type, position, velocity, force, etc.):
 
-      .. code-block:: Python
+      .. code-block:: python
 
          # access first atom
          L.atoms[0].id
@@ -71,7 +71,7 @@ against invalid accesses.
 
       Some attributes can be changed:
 
-      .. code-block:: Python
+      .. code-block:: python
 
          # set position in 2D simulation
          L.atoms[0].position = (1.0, 0.0)

doc/src/Python_config.rst

@@ -4,7 +4,7 @@ Configuration information
 The following methods can be used to query the LAMMPS library
 about compile time settings and included packages and styles.
 
-.. code-block:: Python
+.. code-block:: python
    :caption: Example for using configuration settings functions
 
    from lammps import lammps

doc/src/Python_create.rst

@@ -74,7 +74,7 @@ Here are simple examples using all three Python interfaces:
       :py:class:`PyLammps <lammps.PyLammps>` objects can also be created on top of an existing
       :py:class:`lammps <lammps.lammps>` object:
 
-      .. code-block:: Python
+      .. code-block:: python
 
          from lammps import lammps, PyLammps
          ...
@@ -113,7 +113,7 @@ Here are simple examples using all three Python interfaces:
 
       You can also initialize IPyLammps on top of an existing :py:class:`lammps` or :py:class:`PyLammps` object:
 
-      .. code-block:: Python
+      .. code-block:: python
 
          from lammps import lammps, IPyLammps
          ...
@@ -142,7 +142,7 @@ the MPI and/or Kokkos environment if enabled and active.
 Note that you can create multiple LAMMPS objects in your Python
 script, and coordinate and run multiple simulations, e.g.
 
-.. code-block:: Python
+.. code-block:: python
 
    from lammps import lammps
    lmp1 = lammps()

doc/src/Python_error.rst

@@ -7,7 +7,7 @@ current Python process with an error message.  C++ exceptions allow capturing
 them on the C++ side and rethrowing them on the Python side. This way
 LAMMPS errors can be handled through the Python exception handling mechanism.
 
-.. code-block:: Python
+.. code-block:: python
 
    from lammps import lammps, MPIAbortException
 

doc/src/Python_execute.rst

@@ -60,7 +60,7 @@ it is possible to "compute" what the next LAMMPS command should be.
       can be executed using with the lammps API with the following Python code if ``lmp`` is an
       instance of :py:class:`lammps <lammps.lammps>`:
 
-      .. code-block:: Python
+      .. code-block:: python
 
          from lammps import lammps
 
@@ -73,7 +73,7 @@ it is possible to "compute" what the next LAMMPS command should be.
       The arguments of the command can be passed as one string, or
       individually.
 
-      .. code-block:: Python
+      .. code-block:: python
 
          from lammps import PyLammps
 
@@ -93,14 +93,14 @@ it is possible to "compute" what the next LAMMPS command should be.
       parameterization. In the lammps API parameterization needed to be done
       manually by creating formatted command strings.
 
-      .. code-block:: Python
+      .. code-block:: python
 
          lmp.command("region box block %f %f %f %f %f %f" % (xlo, xhi, ylo, yhi, zlo, zhi))
 
       In contrast, methods of PyLammps accept parameters directly and will convert
       them automatically to a final command string.
 
-      .. code-block:: Python
+      .. code-block:: python
 
          L.region("box block", xlo, xhi, ylo, yhi, zlo, zhi)
 

doc/src/Python_launch.rst

@@ -56,7 +56,7 @@ and you should see the same output as if you had typed
 
 Note that without the mpi4py specific lines from ``test.py``
 
-.. code-block:: Python
+.. code-block:: python
 
    from lammps import lammps
    lmp = lammps()

doc/src/Python_objects.rst

@@ -76,7 +76,7 @@ computes, fixes, or variables in LAMMPS using the :py:mod:`lammps` module.
 
       To define a variable you can use the :doc:`variable <variable>` command:
 
-      .. code-block:: Python
+      .. code-block:: python
 
          L.variable("a index 2")
 
@@ -85,14 +85,14 @@ computes, fixes, or variables in LAMMPS using the :py:mod:`lammps` module.
       you can access an individual variable by retrieving a variable object from the
       ``L.variables`` dictionary by name
 
-      .. code-block:: Python
+      .. code-block:: python
 
          a = L.variables['a']
 
       The variable value can then be easily read and written by accessing the value
       property of this object.
 
-      .. code-block:: Python
+      .. code-block:: python
 
          print(a.value)
          a.value = 4

doc/src/Python_properties.rst

@@ -105,7 +105,7 @@ against invalid accesses.
       variables, compute or fix data (see :doc:`Howto_output`):
 
 
-      .. code-block:: Python
+      .. code-block:: python
 
          result = L.eval("ke") # kinetic energy
          result = L.eval("pe") # potential energy

doc/src/Python_scatter.rst

@@ -1,7 +1,7 @@
 Scatter/gather operations
 =========================
 
-.. code-block:: Python
+.. code-block:: python
 
    data = lmp.gather_atoms(name,type,count)  # return per-atom property of all atoms gathered into data, ordered by atom ID
                                              # name = "x", "charge", "type", etc
@@ -42,7 +42,7 @@ For the scatter methods, the array of coordinates passed to must be a
 ctypes vector of ints or doubles, allocated and initialized something
 like this:
 
-.. code-block:: Python
+.. code-block:: python
 
    from ctypes import c_double
    natoms = lmp.get_natoms()

doc/src/Run_options.rst

@@ -262,6 +262,8 @@ Disable generating a citation reminder (see above) at all.
 
 **-nonbuf**
 
+.. versionadded:: 15Sep2022
+
 Turn off buffering for screen and logfile output.  For performance
 reasons, output to the screen and logfile is usually buffered, i.e.
 output is only written to a file if its buffer - typically 4096 bytes -

doc/src/Tools.rst

@@ -57,7 +57,7 @@ Pre-processing tools
    * :ref:`msi2lmp <msi>`
    * :ref:`polybond <polybond>`
    * :ref:`stl_bin2txt <stlconvert>`
-
+   * :ref:`tabulate <tabulate>`
 
 Post-processing tools
 =====================
@@ -1159,13 +1159,27 @@ For illustration purposes below is a part of the Tcl example script.
 
 ----------
 
+.. _tabulate:
+
+tabulate tool
+--------------
+
+.. versionadded:: TBD
+
+The ``tabulate`` folder contains Python scripts scripts to generate tabulated
+potential files for LAMMPS.  The bulk of the code is in the ``tabulate`` module
+in the ``tabulate.py`` file.  Some example files demonstrating its use are
+included.  See the README file for more information.
+
+----------
+
 .. _vim:
 
 vim tool
 ------------------
 
-The files in the tools/vim directory are add-ons to the VIM editor
+The files in the ``tools/vim`` directory are add-ons to the VIM editor
-that allow easier editing of LAMMPS input scripts.  See the README.txt
+that allow easier editing of LAMMPS input scripts.  See the ``README.txt``
 file for details.
 
 These files were provided by Gerolf Ziegenhain (gerolf at

doc/src/angle_gaussian.rst

@@ -25,23 +25,25 @@ The *gaussian* angle style uses the potential:
 
 .. math::
 
-   E = -k_B T ln\left(\sum_{i=1}^{n} \frac{A_i}{w_i \sqrt{\pi/2}} exp\left( \frac{-(\theta-\theta_{i})^2}{w_i^2})\right) \right)
+   E = -k_B T ln\left(\sum_{i=1}^{n} \frac{A_i}{w_i \sqrt{\pi/2}} exp\left( \frac{-2(\theta-\theta_{i})^2}{w_i^2}\right) \right)
+
+This analytical form is a suitable potential for obtaining mesoscale
+effective force fields which can reproduce target atomistic
+distributions :ref:`(Milano) <Milano1>`.
 
-This analytical form is a suitable potential for obtaining
-mesoscale effective force fields which can reproduce target atomistic distributions :ref:`(Milano) <Milano1>`
 The following coefficients must be defined for each angle type via the
 :doc:`angle_coeff <angle_coeff>` command as in the example above, or in
 the data file or restart files read by the :doc:`read_data <read_data>`
 or :doc:`read_restart <read_restart>` commands:
 
-* T temperature at which the potential was derived
+* :math:`T` temperature at which the potential was derived
 * :math:`n` (integer >=1)
-* :math:`A_1` (-)
+* :math:`A_1` (> 0, radians)
-* :math:`w_1` (-)
+* :math:`w_1` (> 0, radians)
 * :math:`\theta_1` (degrees)
 * ...
-* :math:`A_n` (-)
+* :math:`A_n` (> 0, radians)
-* :math:`w_n` (-)
+* :math:`w_n` (> 0, radians)
 * :math:`\theta_n` (degrees)
 
 

doc/src/angle_table.rst

@@ -59,6 +59,10 @@ format of this file is described below.
 
 ----------
 
+Suitable tables for use with this angle style can be created using the
+Python code in the ``tools/tabulate`` folder of the LAMMPS source code
+distribution.
+
 The format of a tabulated file is as follows (without the
 parenthesized comments):
 

doc/src/bond_bpm_rotational.rst

@@ -69,7 +69,7 @@ which is proportional to the tangential shear displacement with a
 stiffness of :math:`k_s`. This tangential force also induces a torque.
 In addition, bending and twisting torques are also applied to
 particles which are proportional to angular bending and twisting
-displacements with stiffnesses of :math`k_b` and :math:`k_t',
+displacements with stiffnesses of :math:`k_b` and :math:`k_t`,
 respectively.  Details on the calculations of shear displacements and
 angular displacements can be found in :ref:`(Wang) <Wang2009>` and
 :ref:`(Wang and Mora) <Wang2009b>`.

doc/src/bond_gaussian.rst

@@ -25,33 +25,34 @@ The *gaussian* bond style uses the potential:
 
 .. math::
 
-   E = -k_B T ln\left(\sum_{i=1}^{n} \frac{A_i}{w_i \sqrt{\pi/2}} exp\left( \frac{-(r-r_{i})^2}{w_i^2})\right) \right)
+   E = -k_B T ln\left(\sum_{i=1}^{n} \frac{A_i}{w_i \sqrt{\pi/2}} exp\left( \frac{-2(r-r_{i})^2}{w_i^2}\right)\right)
 
-This analytical form is a suitable potential for obtaining
+This analytical form is a suitable potential for obtaining mesoscale
-mesoscale effective force fields which can reproduce target atomistic distributions :ref:`(Milano) <Milano0>`
+effective force fields which can reproduce target atomistic
+distributions :ref:`(Milano) <Milano0>`
 
 The following coefficients must be defined for each bond type via the
 :doc:`bond_coeff <bond_coeff>` command as in the example above, or in
 the data file or restart files read by the :doc:`read_data <read_data>`
 or :doc:`read_restart <read_restart>` commands:
 
-* T temperature at which the potential was derived
+* :math:`T` temperature at which the potential was derived
 * :math:`n` (integer >=1)
-* :math:`A_1` (-)
+* :math:`A_1` (> 0, distance)
-* :math:`w_1` (-)
+* :math:`w_1` (> 0, distance)
-* :math:`r_1` (length)
+* :math:`r_1` (>= 0, distance)
 * ...
-* :math:`A_n` (-)
+* :math:`A_n` (> 0, distance)
-* :math:`w_n` (-)
+* :math:`w_n` (> 0, distance)
-* :math:`r_n` (length)
+* :math:`r_n` (>= 0, distance)
 
 
 Restrictions
 """"""""""""
 
 This bond style can only be used if LAMMPS was built with the
-EXTRA-MOLECULE package.  See the :doc:`Build package <Build_package>` doc
+EXTRA-MOLECULE package.  See the :doc:`Build package <Build_package>`
-page for more info.
+doc page for more info.
 
 Related commands
 """"""""""""""""

doc/src/bond_table.rst

@@ -59,6 +59,13 @@ this file is described below.
 
 ----------
 
+Suitable tables for use with this bond style can be created by LAMMPS
+itself from existing bond styles using the :doc:`bond_write
+<bond_write>` command.  This can be useful to have a template file for
+testing the bond style settings and to build a compatible custom file.
+Another option to generate tables is the Python code in the
+``tools/tabulate`` folder of the LAMMPS source code distribution.
+
 The format of a tabulated file is as follows (without the
 parenthesized comments):
 
@@ -149,7 +156,8 @@ info.
 Related commands
 """"""""""""""""
 
-:doc:`bond_coeff <bond_coeff>`, :doc:`delete_bonds <delete_bonds>`
+:doc:`bond_coeff <bond_coeff>`, :doc:`delete_bonds <delete_bonds>`,
+:doc:`bond_write <bond_write>`
 
 Default
 """""""

doc/src/box.rst

@@ -1,70 +0,0 @@
-.. index:: box
-
-box command
-===========
-
-Syntax
-""""""
-
-.. code-block:: LAMMPS
-
-   box keyword value ...
-
-* one or more keyword/value pairs may be appended
-* keyword = *tilt*
-
-  .. parsed-literal::
-
-       *tilt* value = *small* or *large*
-
-Examples
-""""""""
-
-.. code-block:: LAMMPS
-
-   box tilt large
-   box tilt small
-
-Description
-"""""""""""
-
-Set attributes of the simulation box.
-
-For triclinic (non-orthogonal) simulation boxes, the *tilt* keyword
-allows simulation domains to be created with arbitrary tilt factors,
-e.g. via the :doc:`create_box <create_box>` or
-:doc:`read_data <read_data>` commands.  Tilt factors determine how
-skewed the triclinic box is; see the :doc:`Howto triclinic <Howto_triclinic>` page for a discussion of triclinic
-boxes in LAMMPS.
-
-LAMMPS normally requires that no tilt factor can skew the box more
-than half the distance of the parallel box length, which is the first
-dimension in the tilt factor (x for xz).  If *tilt* is set to
-*small*, which is the default, then an error will be
-generated if a box is created which exceeds this limit.  If *tilt*
-is set to *large*, then no limit is enforced.  You can create
-a box with any tilt factors you wish.
-
-Note that if a simulation box has a large tilt factor, LAMMPS will run
-less efficiently, due to the large volume of communication needed to
-acquire ghost atoms around a processor's irregular-shaped sub-domain.
-For extreme values of tilt, LAMMPS may also lose atoms and generate an
-error.
-
-Restrictions
-""""""""""""
-
-This command cannot be used after the simulation box is defined by a
-:doc:`read_data <read_data>` or :doc:`create_box <create_box>` command or
-:doc:`read_restart <read_restart>` command.
-
-Related commands
-""""""""""""""""
-
-none
-
-
-Default
-"""""""
-
-The default value is tilt = small.

doc/src/commands_list.rst

@@ -13,7 +13,6 @@ Commands
    bond_style
    bond_write
    boundary
-   box
    change_box
    clear
    comm_modify
@@ -43,6 +42,7 @@ Commands
    echo
    fix
    fix_modify
+   fitpod_command
    group
    group2ndx
    hyper
@@ -90,8 +90,7 @@ Commands
    region
    replicate
    rerun
-   reset_atom_ids
+   reset_atoms
-   reset_mol_ids
    reset_timestep
    restart
    run

doc/src/compute_pair_local.rst

@@ -59,7 +59,7 @@ commands.
 The value *dist* is the distance between the pair of atoms.
 The values *dx*, *dy*, and *dz* are the :math:`(x,y,z)` components of the
 *distance* between the pair of atoms. This value is always the
-distance from the atom of lower to the one with the higher id.
+distance from the atom of higher to the one with the lower atom ID.
 
 The value *eng* is the interaction energy for the pair of atoms.
 

doc/src/create_box.rst

@@ -66,20 +66,21 @@ positive or negative values and are called "tilt factors" because they
 are the amount of displacement applied to faces of an originally
 orthogonal box to transform it into the parallelepiped.
 
-By default, a *prism* region used with the create_box command must
+By default, a *prism* region used with the create_box command must have
-have tilt factors :math:`(xy,xz,yz)` that do not skew the box more than half
+tilt factors :math:`(xy,xz,yz)` that do not skew the box more than half
 the distance of the parallel box length.  For example, if
 :math:`x_\text{lo} = 2` and :math:`x_\text{hi} = 12`, then the :math:`x`
-box length is 10 and the :math:`xy` tilt factor must be between :math:`-5` and
+box length is 10 and the :math:`xy` tilt factor must be between
-:math:`5`.  Similarly, both :math:`xz` and :math:`yz` must be between
+:math:`-5` and :math:`5`.  Similarly, both :math:`xz` and :math:`yz`
-:math:`-(x_\text{hi}-x_\text{lo})/2` and :math:`+(y_\text{hi}-y_\text{lo})/2`.
+must be between :math:`-(x_\text{hi}-x_\text{lo})/2` and
-Note that this is not a limitation, since if the maximum tilt factor is 5 (as
+:math:`+(y_\text{hi}-y_\text{lo})/2`.  Note that this is not a
-in this example), then configurations with tilt :math:`= \dots, -15`,
+limitation, since if the maximum tilt factor is 5 (as in this example),
-:math:`-5`, :math:`5`, :math:`15`, :math:`25, \dots`
+then configurations with tilt :math:`= \dots, -15`, :math:`-5`,
-are all geometrically equivalent.  If you wish to define a box with tilt
+:math:`5`, :math:`15`, :math:`25, \dots` are all geometrically
-factors that exceed these limits, you can use the :doc:`box tilt <box>`
+equivalent.  Simulations with large tilt factors will run inefficiently,
-command, with a setting of *large*\ ; a setting of *small* is the
+since they require more ghost atoms and thus more communication.  With
-default.
+very large tilt factors, LAMMPS will eventually produce incorrect
+trajectories and stop with errors due to lost atoms or similar.
 
 See the :doc:`Howto triclinic <Howto_triclinic>` page for a
 geometric description of triclinic boxes, as defined by LAMMPS, and

doc/src/delete_atoms.rst

@@ -135,7 +135,7 @@ number of atoms in the system.  Note that this is not done for
 molecular systems (see the :doc:`atom_style <atom_style>` command),
 regardless of the *compress* setting, since it would foul up the bond
 connectivity that has already been assigned.  However, the
-:doc:`reset_atom_ids <reset_atom_ids>` command can be used after this
+:doc:`reset_atoms id <reset_atoms>` command can be used after this
 command to accomplish the same thing.
 
 Note that the re-assignment of IDs is not really a compression, where
@@ -203,7 +203,7 @@ using molecule template files via the :doc:`molecule <molecule>` and
 Related commands
 """"""""""""""""
 
-:doc:`create_atoms <create_atoms>`, :doc:`reset_atom_ids <reset_atom_ids>`
+:doc:`create_atoms <create_atoms>`, :doc:`reset_atoms id <reset_atoms>`
 
 Default
 """""""

doc/src/dihedral_table.rst

@@ -114,6 +114,10 @@ below.
 
 ----------
 
+Suitable tables for use with this dihedral style can be created using
+the Python code in the ``tools/tabulate`` folder of the LAMMPS source
+code distribution.
+
 The format of a tabulated file is as follows (without the
 parenthesized comments).  It can begin with one or more comment
 or blank lines.

doc/src/dump_image.rst

@@ -69,7 +69,7 @@ Syntax
          yes/no = do or do not draw simulation box lines
          diam = diameter of box lines as fraction of shortest box length
        *axes* values = axes length diam = draw xyz axes
-         axes = *yes* or *no = do or do not draw xyz axes lines next to simulation box
+         axes = *yes* or *no* = do or do not draw xyz axes lines next to simulation box
          length = length of axes lines as fraction of respective box lengths
          diam = diameter of axes lines as fraction of shortest box length
        *subbox* values = lines diam = draw outline of processor sub-domains

doc/src/fitpod_command.rst

@@ -0,0 +1,745 @@
+.. index:: fitpod
+
+fitpod command
+======================
+
+Syntax
+""""""
+
+.. parsed-literal::
+
+   fitpod Ta_param.pod Ta_data.pod
+
+* fitpod = style name of this command
+* Ta_param.pod = an input file that describes proper orthogonal descriptors (PODs)
+* Ta_data.pod = an input file that specifies DFT data used to fit a POD potential
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   fitpod Ta_param.pod Ta_data.pod
+
+Description
+"""""""""""
+.. versionadded:: TBD
+
+Fit a machine-learning interatomic potential (ML-IAP) based on proper
+orthogonal descriptors (POD).  Two input files are required for this
+command. The first input file describes a POD potential parameter
+settings, while the second input file specifies the DFT data used for
+the fitting procedure.
+
+The table below has one-line descriptions of all the keywords that can
+be used in the first input file (i.e. ``Ta_param.pod`` in the example
+above):
+
+.. list-table::
+   :header-rows: 1
+   :widths: auto
+
+   * - Keyword
+     - Default
+     - Type
+     - Description
+   * - species
+     - (none)
+     - STRING
+     - Chemical symbols for all elements in the system and have to match XYZ training files.
+   * - pbc
+     - 1 1 1
+     - INT
+     - three integer constants specify boundary conditions
+   * - rin
+     - 1.0
+     - REAL
+     - a real number specifies the inner cut-off radius
+   * - rcut
+     - 5.0
+     - REAL
+     - a real number specifies the outer cut-off radius
+   * - bessel_polynomial_degree
+     - 3
+     - INT
+     - the maximum degree of Bessel polynomials
+   * - inverse_polynomial_degree
+     - 6
+     - INT
+     - the maximum degree of inverse radial basis functions
+   * - onebody
+     - 1
+     - BOOL
+     - turns on/off one-body potential
+   * - twobody_number_radial_basis_functions
+     - 6
+     - INT
+     - number of radial basis functions for two-body potential
+   * - threebody_number_radial_basis_functions
+     - 5
+     - INT
+     - number of radial basis functions for three-body potential
+   * - threebody_number_angular_basis_functions
+     - 5
+     - INT
+     - number of angular basis functions for three-body potential
+   * - fourbody_snap_twojmax
+     - 0
+     - INT
+     - band limit for SNAP bispectrum components (0,2,4,6,8... allowed)
+   * - fourbody_snap_chemflag
+     - 0
+     - BOOL
+     - turns on/off the explicit multi-element variant of the SNAP bispectrum components
+   * - quadratic_pod_potential
+     - 0
+     - BOOL
+     - turns on/off quadratic POD potential
+
+All keywords except *species* have default values. If a keyword is not
+set in the input file, its default value is used.  The next table
+describes all keywords that can be used in the second input file
+(i.e. ``Ta_data.pod`` in the example above):
+
+.. list-table::
+   :header-rows: 1
+   :widths: auto
+
+   * - Keyword
+     - Default
+     - Type
+     - Description
+   * - file_format
+     - extxyz
+     - STRING
+     - only the extended xyz format (extxyz) is currently supported
+   * - file_extension
+     - xyz
+     - STRING
+     - extension of the data files
+   * - path_to_training_data_set
+     - (none)
+     - STRING
+     - specifies the path to training data files in double quotes
+   * - path_to_test_data_set
+     - ""
+     - STRING
+     - specifies the path to test data files in double quotes
+   * - fraction_training_data_set
+     - 1.0
+     - REAL
+     - a real number (<= 1.0) specifies the fraction of the training set used to fit POD
+   * - randomize_training_data_set
+     - 0
+     - BOOL
+     - turns on/off randomization of the training set
+   * - fitting_weight_energy
+     - 100.0
+     - REAL
+     - a real constant specifies the weight for energy in the least-squares fit
+   * - fitting_weight_force
+     - 1.0
+     - REAL
+     - a real constant specifies the weight for force in the least-squares fit
+   * - fitting_regularization_parameter
+     - 1.0e-10
+     - REAL
+     - a real constant specifies the regularization parameter in the least-squares fit
+   * - error_analysis_for_training_data_set
+     - 0
+     - BOOL
+     - turns on/off error analysis for the training data set
+   * - error_analysis_for_test_data_set
+     - 0
+     - BOOL
+     - turns on/off error analysis for the test data set
+   * - basename_for_output_files
+     - pod
+     - STRING
+     - a basename string added to the output files
+   * - precision_for_pod_coefficients
+     - 8
+     - INT
+     - number of digits after the decimal points for numbers in the coefficient file
+
+All keywords except *path_to_training_data_set* have default values. If
+a keyword is not set in the input file, its default value is used.  After
+successful training, a number of output files are produced, if enabled:
+
+* ``<basename>_training_errors.pod``  reports the errors in energy and forces for the training data set
+* ``<basename>_training_analysis.pod`` reports detailed errors for all training configurations
+* ``<basename>_test_errors.pod`` reports errors for the test data set
+* ``<basename>_test_analysis.pod`` reports detailed errors for all test configurations
+* ``<basename>_coefficients.pod`` contains the coefficients of the POD potential
+
+After training the POD potential, ``Ta_param.pod`` and ``<basename>_coefficients.pod``
+are the two files needed to use the POD potential in LAMMPS. See
+:doc:`pair_style pod <pair_pod>` for using the POD potential. Examples
+about training and using POD potentials are found in the directory
+lammps/examples/PACKAGES/pod.
+
+Parameterized Potential Energy Surface
+""""""""""""""""""""""""""""""""""""""
+
+We consider a multi-element system of *N* atoms with :math:`N_{\rm e}`
+unique elements.  We denote by :math:`\boldsymbol r_n` and :math:`Z_n`
+position vector and type of an atom *n* in the system,
+respectively. Note that we have :math:`Z_n \in \{1, \ldots, N_{\rm e}
+\}`, :math:`\boldsymbol R = (\boldsymbol r_1, \boldsymbol r_2, \ldots,
+\boldsymbol r_N) \in \mathbb{R}^{3N}`, and :math:`\boldsymbol Z = (Z_1,
+Z_2, \ldots, Z_N) \in \mathbb{N}^{N}`. The potential energy surface
+(PES) of the system can be expressed as a many-body expansion of the
+form
+
+.. math::
+
+    E(\boldsymbol R, \boldsymbol Z, \boldsymbol{\eta}, \boldsymbol{\mu}) \ = \ & \sum_{i} V^{(1)}(\boldsymbol r_i, Z_i, \boldsymbol \mu^{(1)} ) + \frac12 \sum_{i,j} V^{(2)}(\boldsymbol r_i, \boldsymbol r_j, Z_i, Z_j, \boldsymbol \eta, \boldsymbol \mu^{(2)})  \\
+    & + \frac16 \sum_{i,j,k} V^{(3)}(\boldsymbol r_i, \boldsymbol r_j, \boldsymbol r_k, Z_i, Z_j, Z_k, \boldsymbol \eta, \boldsymbol \mu^{(3)}) + \ldots
+
+where :math:`V^{(1)}` is the one-body potential often used for
+representing external field or energy of isolated elements, and the
+higher-body potentials :math:`V^{(2)}, V^{(3)}, \ldots` are symmetric,
+uniquely defined, and zero if two or more indices take identical values.
+The superscript on each potential denotes its body order. Each *q*-body
+potential :math:`V^{(q)}` depends on :math:`\boldsymbol \mu^{(q)}` which
+are sets of parameters to fit the PES. Note that :math:`\boldsymbol \mu`
+is a collection of all potential parameters :math:`\boldsymbol
+\mu^{(1)}`, :math:`\boldsymbol \mu^{(2)}`, :math:`\boldsymbol
+\mu^{(3)}`, etc, and that :math:`\boldsymbol \eta` is a set of
+hyper-parameters such as inner cut-off radius :math:`r_{\rm in}` and
+outer cut-off radius :math:`r_{\rm cut}`.
+
+Interatomic potentials rely on parameters to learn relationship between
+atomic environments and interactions.  Since interatomic potentials are
+approximations by nature, their parameters need to be set to some
+reference values or fitted against data by necessity.  Typically,
+potential fitting finds optimal parameters, :math:`\boldsymbol \mu^*`,
+to minimize a certain loss function of the predicted quantities and
+data. Since the fitted potential depends on the data set used to fit it,
+different data sets will yield different optimal parameters and thus
+different fitted potentials. When fitting the same functional form on
+*Q* different data sets, we would obtain *Q* different optimized
+potentials, :math:`E(\boldsymbol R,\boldsymbol Z, \boldsymbol \eta,
+\boldsymbol \mu_q^*), 1 \le q \le Q`.  Consequently, there exist many
+different sets of optimized parameters for empirical interatomic
+potentials.
+
+Instead of optimizing the potential parameters, inspired by the reduced
+basis method :ref:`(Grepl) <Grepl20072>` for parameterized partial
+differential equations, we view the parameterized PES as a parametric
+manifold of potential energies
+
+.. math::
+
+    \mathcal{M} = \{E(\boldsymbol R, \boldsymbol Z, \boldsymbol \eta, \boldsymbol \mu) \ | \  \boldsymbol \mu \in \Omega^{\boldsymbol \mu} \}
+
+where :math:`\Omega^{\boldsymbol \mu}` is a parameter domain in which
+:math:`\boldsymbol \mu` resides.  The parametric manifold
+:math:`\mathcal{M}` contains potential energy surfaces for all values of
+:math:`\boldsymbol \mu \in \Omega^{\boldsymbol \mu}`.  Therefore, the
+parametric manifold yields a much richer and more transferable atomic
+representation than any particular individual PES :math:`E(\boldsymbol
+R, \boldsymbol Z, \boldsymbol \eta, \boldsymbol \mu^*)`.
+
+We propose specific forms of the parameterized potentials for one-body,
+two-body, and three-body interactions. We apply the Karhunen-Loeve
+expansion to snapshots of the parameterized potentials to obtain sets of
+orthogonal basis functions. These basis functions are aggregated
+according to the chemical elements of atoms, thus leading to
+multi-element proper orthogonal descriptors.
+
+Proper Orthogonal Descriptors
+"""""""""""""""""""""""""""""
+
+Proper orthogonal descriptors are finger prints characterizing the
+radial and angular distribution of a system of atoms. The detailed
+mathematical definition is given in the paper by Nguyen and Rohskopf
+:ref:`(Nguyen) <Nguyen20222>`.
+
+The descriptors for the one-body interaction are used to capture energy
+of isolated elements and defined as follows
+
+.. math::
+
+    D_{ip}^{(1)} =  \left\{
+        \begin{array}{ll}
+        1, & \mbox{if } Z_i = p \\
+        0, & \mbox{if } Z_i \neq p
+        \end{array}
+    \right.
+
+for :math:`1 \le i \le N, 1 \le p \le N_{\rm e}`. The number of one-body
+descriptors per atom is equal to the number of elements. The one-body
+descriptors are independent of atom positions, but dependent on atom
+types. The one-body descriptors are active only when the keyword
+*onebody* is set to 1.
+
+We adopt the usual assumption that the direct interaction between two
+atoms vanishes smoothly when their distance is greater than the outer
+cutoff distance :math:`r_{\rm cut}`. Furthermore, we assume that two
+atoms can not get closer than the inner cutoff distance :math:`r_{\rm
+in}` due to the Pauli repulsion principle. Let :math:`r \in (r_{\rm in},
+r_{\rm cut})`, we introduce the following parameterized radial functions
+
+.. math::
+
+    \phi(r, r_{\rm in}, r_{\rm cut}, \alpha, \beta)  = \frac{\sin (\alpha \pi x) }{r - r_{\rm in}}, \qquad  \varphi(r, \gamma)  = \frac{1}{r^\gamma} ,
+
+where the scaled distance function :math:`x` is defined below to enrich the two-body manifold
+
+.. math::
+
+    x(r, r_{\rm in}, r_{\rm cut}, \beta) = \frac{e^{-\beta(r - r_{\rm in})/(r_{\rm cut} - r_{\rm in})} - 1}{e^{-\beta} - 1} .
+
+We introduce the following function as a convex combination of the two functions
+
+.. math::
+
+    \psi(r, r_{\rm in}, r_{\rm cut}, \alpha, \beta, \gamma, \kappa)  = \kappa \phi(r, r_{\rm in}, r_{\rm cut}, \alpha, \beta) + (1- \kappa)  \varphi(r, \gamma) .
+
+We see that :math:`\psi` is a function of distance :math:`r`, cut-off
+distances :math:`r_{\rm in}` and :math:`r_{\rm cut}`, and parameters
+:math:`\alpha, \beta, \gamma, \kappa`. Together these parameters allow
+the function :math:`\psi` to characterize a diverse spectrum of two-body
+interactions within the cut-off interval :math:`(r_{\rm in}, r_{\rm
+cut})`.
+
+Next, we introduce the following parameterized potential
+
+.. math::
+
+    W^{(2)}(r_{ij}, \boldsymbol \eta, \boldsymbol \mu^{(2)})  = f_{\rm c}(r_{ij}, \boldsymbol \eta) \psi(r_{ij}, \boldsymbol \eta, \boldsymbol \mu^{(2)})
+
+where :math:`\eta_1 = r_{\rm in}, \eta_2 = r_{\rm cut}, \mu_1^{(2)} =
+\alpha, \mu_2^{(2)} = \beta, \mu_3^{(2)} = \gamma`, and
+:math:`\mu_4^{(2)} = \kappa`. Here the cut-off function :math:`f_{\rm
+c}(r_{ij}, \boldsymbol \eta)` proposed in [refs] is used to ensure the
+smooth vanishing of the potential and its derivative for :math:`r_{ij}
+\ge r_{\rm cut}`:
+
+.. math::
+
+    f_{\rm c}(r_{ij},  r_{\rm in}, r_{\rm cut})  =  \exp \left(1 -\frac{1}{\sqrt{\left(1 - \frac{(r-r_{\rm in})^3}{(r_{\rm cut} - r_{\rm in})^3} \right)^2 + 10^{-6}}} \right)
+
+Based on the parameterized potential, we form a set of snapshots as
+follows.  We assume that we are given :math:`N_{\rm s}` parameter tuples
+:math:`\boldsymbol \mu^{(2)}_\ell, 1 \le \ell \le N_{\rm s}`. We
+introduce the following set of snapshots on :math:`(r_{\rm in}, r_{\rm
+cut})`:
+
+.. math::
+
+    \xi_\ell(r_{ij}, \boldsymbol \eta) =  W^{(2)}(r_{ij}, \boldsymbol \eta, \boldsymbol \mu^{(2)}_\ell),  \quad \ell = 1, \ldots, N_{\rm s} .
+
+To ensure adequate sampling of the PES for different parameters, we
+choose :math:`N_{\rm s}` parameter points :math:`\boldsymbol
+\mu^{(2)}_\ell = (\alpha_\ell, \beta_\ell, \gamma_\ell, \kappa_\ell), 1
+\le \ell \le N_{\rm s}` as follows. The parameters :math:`\alpha \in [1,
+N_\alpha]` and :math:`\gamma \in [1, N_\gamma]` are integers, where
+:math:`N_\alpha` and :math:`N_\gamma` are the highest degrees for
+:math:`\alpha` and :math:`\gamma`, respectively. We next choose
+:math:`N_\beta` different values of :math:`\beta` in the interval
+:math:`[\beta_{\min}, \beta_{\max}]`, where :math:`\beta_{\min} = 0` and
+:math:`\beta_{\max} = 4`. The parameter :math:`\kappa` can be set either
+0 or 1.  Hence, the total number of parameter points is :math:`N_{\rm s}
+= N_\alpha N_\beta + N_\gamma`.  Although :math:`N_\alpha, N_\beta,
+N_\gamma` can be chosen conservatively large, we find that
+:math:`N_\alpha = 6, N_\beta = 3, N_\gamma = 8` are adequate for most
+problems.  Note that :math:`N_\alpha` and :math:`N_\gamma` correspond to
+*bessel_polynomial_degree* and *inverse_polynomial_degree*,
+respectively.
+
+We employ the Karhunen-Loeve (KL) expansion to generate an orthogonal
+basis set which is known to be optimal for representation of the
+snapshot family :math:`\{\xi_\ell\}_{\ell=1}^{N_{\rm s}}`. The two-body
+orthogonal basis functions are computed as follows
+
+.. math::
+
+    U^{(2)}_m(r_{ij}, \boldsymbol \eta) = \sum_{\ell = 1}^{N_{\rm s}} A_{\ell m}(\boldsymbol \eta) \,  \xi_\ell(r_{ij}, \boldsymbol \eta), \qquad m = 1, \ldots, N_{\rm 2b} ,
+
+where the matrix :math:`\boldsymbol A \in \mathbb{R}^{N_{\rm s} \times
+N_{\rm s}}` consists of eigenvectors of the eigenvalue problem
+
+.. math::
+
+    \boldsymbol C \boldsymbol a = \lambda \boldsymbol a
+
+with the entries of :math:`\boldsymbol C \in \mathbb{R}^{N_{\rm s} \times N_{\rm s}}` being given by
+
+.. math::
+
+    C_{ij}  = \frac{1}{N_{\rm s}} \int_{r_{\rm in}}^{r_{\rm cut}} \xi_i(x, \boldsymbol \eta) \xi_j(x, \boldsymbol \eta) dx, \quad 1 \le i, j \le N_{\rm s}
+
+Note that the eigenvalues :math:`\lambda_\ell, 1 \le \ell \le N_{\rm
+s}`, are ordered such that :math:`\lambda_1 \ge \lambda_2 \ge \ldots \ge
+\lambda_{N_{\rm s}}`, and that the matrix :math:`\boldsymbol A` is
+pe-computed and stored for any given :math:`\boldsymbol \eta`.  Owing to
+the rapid convergence of the KL expansion, only a small number of
+orthogonal basis functions is needed to obtain accurate
+approximation. The value of :math:`N_{\rm 2b}` corresponds to
+*twobody_number_radial_basis_functions*.
+
+The two-body proper orthogonal descriptors at each atom *i* are computed
+by summing the orthogonal basis functions over the neighbors of atom *i*
+and numerating on the atom types as follows
+
+.. math::
+
+    D^{(2)}_{im l(p, q) }(\boldsymbol \eta)  = \left\{
+    \begin{array}{ll}
+    \displaystyle \sum_{\{j | Z_j = q\}} U^{(2)}_m(r_{ij},  \boldsymbol \eta), & \mbox{if } Z_i = p \\
+    0, & \mbox{if } Z_i \neq p
+    \end{array}
+    \right.
+
+for :math:`1 \le i \le N, 1 \le m \le N_{\rm 2b}, 1 \le q, p \le N_{\rm
+e}`. Here :math:`l(p,q)` is a symmetric index mapping such that
+
+.. math::
+
+    l(p,q)  = \left\{
+    \begin{array}{ll}
+    q + (p-1) N_{\rm e} - p(p-1)/2, & \mbox{if } q \ge p \\
+    p + (q-1) N_{\rm e} - q(q-1)/2, & \mbox{if } q < p .
+    \end{array}
+    \right.
+
+The number of two-body descriptors per atom is thus :math:`N_{\rm 2b}
+N_{\rm e}(N_{\rm e}+1)/2`.
+
+It is important to note that the orthogonal basis functions do not
+depend on the atomic numbers :math:`Z_i` and :math:`Z_j`. Therefore, the
+cost of evaluating the basis functions and their derivatives with
+respect to :math:`r_{ij}` is independent of the number of elements
+:math:`N_{\rm e}`. Consequently, even though the two-body proper
+orthogonal descriptors depend on :math:`\boldsymbol Z`, their
+computational complexity is independent of :math:`N_{\rm e}`.
+
+In order to provide proper orthogonal descriptors for three-body
+interactions, we need to introduce a three-body parameterized
+potential. In particular, the three-body potential is defined as a
+product of radial and angular functions as follows
+
+.. math::
+
+    W^{(3)}(r_{ij}, r_{ik}, \theta_{ijk}, \boldsymbol \eta, \boldsymbol \mu^{(3)})  =  \psi(r_{ij}, r_{\rm min}, r_{\rm max}, \alpha, \beta, \gamma, \kappa) f_{\rm c}(r_{ij}, r_{\rm min}, r_{\rm max}) \\
+    \psi(r_{ik}, r_{\rm min}, r_{\rm max}, \alpha, \beta, \gamma, \kappa) f_{\rm c}(r_{ik}, r_{\rm min}, r_{\rm max}) \\
+    \cos (\sigma \theta_{ijk} + \zeta)
+
+where :math:`\sigma` is the periodic multiplicity, :math:`\zeta` is the
+equilibrium angle, :math:`\boldsymbol \mu^{(3)} = (\alpha, \beta,
+\gamma, \kappa, \sigma, \zeta)`. The three-body potential provides an
+angular fingerprint of the atomic environment through the bond angles
+:math:`\theta_{ijk}` formed with each pair of neighbors :math:`j` and
+:math:`k`.  Compared to the two-body potential, the three-body potential
+has two extra parameters :math:`(\sigma, \zeta)` associated with the
+angular component.
+
+Let :math:`\boldsymbol \varrho = (\alpha, \beta, \gamma, \kappa)`. We
+assume that we are given :math:`L_{\rm r}` parameter tuples
+:math:`\boldsymbol \varrho_\ell, 1 \le \ell \le L_{\rm r}`.  We
+introduce the following set of snapshots on :math:`(r_{\min},
+r_{\max})`:
+
+.. math::
+
+    \zeta_\ell(r_{ij}, r_{\rm min}, r_{\rm max} ) =  \psi(r_{ij}, r_{\rm min}, r_{\rm max}, \boldsymbol \varrho_\ell) f_{\rm c}(r_{ij}, r_{\rm min},  r_{\rm max}), \quad 1 \le \ell \le L_{\rm r} .
+
+We apply the Karhunen-Loeve (KL) expansion to this set of snapshots to
+obtain orthogonal basis functions as follows
+
+.. math::
+
+    U^{r}_m(r_{ij}, r_{\rm min}, r_{\rm max} ) = \sum_{\ell = 1}^{L_{\rm r}} A_{\ell m} \,  \zeta_\ell(r_{ij}, r_{\rm min}, r_{\rm max} ), \qquad m = 1, \ldots, N_{\rm r} ,
+
+where the matrix :math:`\boldsymbol A \in \mathbb{R}^{L_{\rm r} \times L_{\rm r}}` consists
+of eigenvectors of the eigenvalue problem. For the parameterized angular function,
+we consider angular basis functions
+
+.. math::
+
+    U^{a}_n(\theta_{ijk}) = \cos ((n-1) \theta_{ijk}), \qquad  n = 1,\ldots, N_{\rm a},
+
+where :math:`N_{\rm a}` is the number of angular basis functions. The orthogonal
+basis functions for the parameterized potential are computed as follows
+
+.. math::
+
+    U^{(3)}_{mn}(r_{ij}, r_{ik}, \theta_{ijk}, \boldsymbol \eta) = U^{r}_m(r_{ij}, \boldsymbol \eta) U^{r}_m(r_{ik}, \boldsymbol \eta) U^{a}_n(\theta_{ijk}),
+
+for :math:`1 \le m \le N_{\rm r}, 1 \le n \le N_{\rm a}`. The number of three-body
+orthogonal basis functions is equal to :math:`N_{\rm 3b} = N_{\rm r} N_{\rm a}` and
+independent of the number of elements. The value of :math:`N_{\rm r}` corresponds to
+*threebody_number_radial_basis_functions*, while that of :math:`N_{\rm a}` to
+*threebody_number_angular_basis_functions*.
+
+The three-body proper orthogonal descriptors at each atom *i*
+are obtained by summing over the neighbors *j* and *k* of atom *i* as
+
+.. math::
+
+    D^{(3)}_{imn \ell(p, q, s)}(\boldsymbol \eta)  = \left\{
+    \begin{array}{ll}
+    \displaystyle \sum_{\{j | Z_j = q\}} \sum_{\{k | Z_k = s\}} U^{(3)}_{mn}(r_{ij}, r_{ik}, \theta_{ijk}, \boldsymbol \eta), & \mbox{if } Z_i = p \\
+    0, & \mbox{if } Z_i \neq p
+    \end{array}
+    \right.
+
+for :math:`1 \le i \le N, 1 \le m \le N_{\rm r}, 1 \le n \le N_{\rm a}, 1 \le q, p, s \le N_{\rm e}`,
+where
+
+.. math::
+
+    \ell(p,q,s)  = \left\{
+    \begin{array}{ll}
+    s + (q-1) N_{\rm e} - q(q-1)/2 + (p-1)N_{\rm e}(1+N_{\rm e})/2 , & \mbox{if } s \ge q \\
+    q + (s-1) N_{\rm e} - s(s-1)/2 + (p-1)N_{\rm e}(1+N_{\rm e})/2, & \mbox{if } s < q .
+    \end{array}
+    \right.
+
+The number of three-body descriptors per atom is thus :math:`N_{\rm 3b} N_{\rm e}^2(N_{\rm e}+1)/2`.
+While the number of three-body PODs is cubic function of the number of elements,
+the computational complexity of the three-body PODs is independent of the number of elements.
+
+Four-Body SNAP Descriptors
+""""""""""""""""""""""""""
+
+In addition to the proper orthogonal descriptors described above, we also employ
+the spectral neighbor analysis potential (SNAP) descriptors. SNAP uses bispectrum components
+to characterize the local neighborhood of each atom in a very general way. The mathematical definition
+of the bispectrum calculation and its derivatives w.r.t. atom positions is described in
+:doc:`compute snap <compute_sna_atom>`. In SNAP, the
+total energy is decomposed into a sum over atom energies. The energy of
+atom *i* is expressed as a weighted sum over bispectrum components.
+
+.. math::
+
+   E_i^{\rm SNAP} = \sum_{k=1}^{N_{\rm 4b}} \sum_{p=1}^{N_{\rm e}} c_{kp}^{(4)} D_{ikp}^{(4)}
+
+
+where the SNAP descriptors are related to the bispectrum components by
+
+.. math::
+
+    D^{(4)}_{ikp}  = \left\{
+    \begin{array}{ll}
+    \displaystyle B_{ik}, & \mbox{if } Z_i = p \\
+    0, & \mbox{if } Z_i \neq p
+    \end{array}
+    \right.
+
+Here :math:`B_{ik}` is the *k*\ -th bispectrum component of atom *i*. The number of
+bispectrum components :math:`N_{\rm 4b}` depends on the value of *fourbody_snap_twojmax* :math:`= 2 J_{\rm max}`
+and *fourbody_snap_chemflag*. If *fourbody_snap_chemflag* = 0
+then :math:`N_{\rm 4b} = (J_{\rm max}+1)(J_{\rm max}+2)(J_{\rm max}+1.5)/3`.
+If *fourbody_snap_chemflag* = 1 then :math:`N_{\rm 4b} = N_{\rm e}^3 (J_{\rm max}+1)(J_{\rm max}+2)(J_{\rm max}+1.5)/3`.
+The bispectrum calculation is described in more detail in :doc:`compute sna/atom <compute_sna_atom>`.
+
+Linear Proper Orthogonal Descriptor Potentials
+""""""""""""""""""""""""""""""""""""""""""""""
+
+The proper orthogonal descriptors and SNAP descriptors are used to define the atomic energies
+in the following expansion
+
+.. math::
+
+    E_{i}(\boldsymbol \eta) = \sum_{p=1}^{N_{\rm e}} c^{(1)}_p D^{(1)}_{ip} + \sum_{m=1}^{N_{\rm 2b}}  \sum_{l=1}^{N_{\rm e}(N_{\rm e}+1)/2} c^{(2)}_{ml} D^{(2)}_{iml}(\boldsymbol \eta) + \sum_{m=1}^{N_{\rm r}} \sum_{n=1}^{N_{\rm a}}  \sum_{\ell=1}^{N_{\rm e}^2(N_{\rm e}+1)/2} c^{(3)}_{mn\ell} D^{(3)}_{imn\ell}(\boldsymbol \eta) + \sum_{k=1}^{N_{\rm 4b}} \sum_{p=1}^{N_{\rm e}} c_{kp}^{(4)} D_{ikp}^{(4)}(\boldsymbol \eta),
+
+where :math:`D^{(1)}_{ip}, D^{(2)}_{iml}, D^{(3)}_{imn\ell}, D^{(4)}_{ikp}` are the  one-body, two-body, three-body, four-body descriptors,
+respectively, and :math:`c^{(1)}_p, c^{(2)}_{ml}, c^{(3)}_{mn\ell}, c^{(4)}_{kp}` are their respective expansion
+coefficients. In a more compact notation that implies summation over descriptor indices
+the atomic energies can be written as
+
+.. math::
+
+    E_i(\boldsymbol \eta) =  \sum_{m=1}^{N_{\rm e}} c^{(1)}_m D^{(1)}_{im} +  \sum_{m=1}^{N_{\rm d}^{(2)}} c^{(2)}_k D^{(2)}_{im} + \sum_{m=1}^{N_{\rm d}^{(3)}} c^{(3)}_m D^{(3)}_{im} + \sum_{m=1}^{N_{\rm d}^{(4)}} c^{(4)}_m D^{(4)}_{im}
+
+where :math:`N_{\rm d}^{(2)} = N_{\rm 2b} N_{\rm e} (N_{\rm e}+1)/2`,
+:math:`N_{\rm d}^{(3)} = N_{\rm 3b} N_{\rm e}^2 (N_{\rm e}+1)/2`, and
+:math:`N_{\rm d}^{(4)} = N_{\rm 4b} N_{\rm e}` are
+the number of two-body, three-body, and four-body descriptors, respectively.
+
+The potential energy is then obtained by summing local atomic energies :math:`E_i`
+for all atoms :math:`i` in the system
+
+.. math::
+
+    E(\boldsymbol \eta) = \sum_{i}^N E_{i}(\boldsymbol \eta)
+
+Because the descriptors are one-body, two-body, and three-body terms,
+the resulting POD potential is a three-body PES. We can express the potential
+energy as a linear combination of the global descriptors as follows
+
+.. math::
+
+    E(\boldsymbol \eta) = \sum_{m=1}^{N_{\rm e}} c^{(1)}_m d^{(1)}_{m} +  \sum_{m=1}^{N_{\rm d}^{(2)}} c^{(2)}_m d^{(2)}_{m} + \sum_{m=1}^{N_{\rm d}^{(3)}} c^{(3)}_m d^{(3)}_{m} + \sum_{m=1}^{N_{\rm d}^{(4)}} c^{(4)}_m d^{(4)}_{m}
+
+where  the global descriptors are given by
+
+.. math::
+
+    d_{m}^{(1)}(\boldsymbol \eta) = \sum_{i=1}^N D_{im}^{(1)}(\boldsymbol \eta), \quad d_{m}^{(2)}(\boldsymbol \eta) = \sum_{i=1}^N D_{im}^{(2)}(\boldsymbol \eta), \quad d_{m}^{(3)}(\boldsymbol \eta) = \sum_{i=1}^N D_{im}^{(3)}(\boldsymbol \eta), \quad d_{m}^{(4)}(\boldsymbol \eta) = \sum_{i=1}^N D_{im}^{(4)}(\boldsymbol \eta)
+
+Hence, we obtain the atomic forces as
+
+.. math::
+
+    \boldsymbol F = -\nabla E(\boldsymbol \eta) = - \sum_{m=1}^{N_{\rm d}^{(2)}}  c^{(2)}_m  \nabla d_m^{(2)} - \sum_{m=1}^{N_{\rm d}^{(3)}}  c^{(3)}_m \nabla d_m^{(3)} - \sum_{m=1}^{N_{\rm d}^{(4)}}  c^{(4)}_m \nabla d_m^{(4)}
+
+where :math:`\nabla d_m^{(2)}`, :math:`\nabla d_m^{(3)}` and :math:`\nabla d_m^{(4)}` are derivatives of the two-body
+three-body, and four-body global descriptors with respect to atom positions, respectively.
+Note that since the first-body global descriptors are constant, their derivatives are zero.
+
+Quadratic Proper Orthogonal Descriptor Potentials
+"""""""""""""""""""""""""""""""""""""""""""""""""
+
+We recall two-body PODs :math:`D^{(2)}_{ik}, 1 \le k \le N_{\rm d}^{(2)}`,
+and three-body PODs :math:`D^{(3)}_{im}, 1 \le m \le N_{\rm d}^{(3)}`,
+with :math:`N_{\rm d}^{(2)} = N_{\rm 2b} N_{\rm e} (N_{\rm e}+1)/2` and
+:math:`N_{\rm d}^{(3)} = N_{\rm 3b} N_{\rm e}^2 (N_{\rm e}+1)/2` being
+the number of descriptors per atom for the two-body PODs and three-body PODs,
+respectively. We employ them to define a new set of atomic descriptors as follows
+
+.. math::
+
+    D^{(2*3)}_{ikm} = \frac{1}{2N}\left( D^{(2)}_{ik} \sum_{j=1}^N D^{(3)}_{jm} + D^{(3)}_{im} \sum_{j=1}^N D^{(2)}_{jk}  \right)
+
+for :math:`1 \le i \le N, 1 \le k \le N_{\rm d}^{(2)}, 1 \le m \le N_{\rm d}^{(3)}`.
+The new descriptors are four-body because they involve central atom :math:`i` together
+with three neighbors :math:`j, k` and :math:`l`. The total number of new  descriptors per atom is equal to
+
+.. math::
+
+    N_{\rm d}^{(2*3)} = N_{\rm d}^{(2)} * N_{\rm d}^{(3)} = N_{\rm 2b} N_{\rm 3b} N_{\rm e}^3 (N_{\rm e}+1)^2/4 .
+
+The new global descriptors are calculated as
+
+.. math::
+
+    d^{(2*3)}_{km} = \sum_{i=1}^N D^{(2*3)}_{ikm} = \left( \sum_{i=1}^N D^{(2)}_{ik} \right) \left( \sum_{i=1}^N D^{(3)}_{im} \right) = d^{(2)}_{k} d^{(3)}_m,
+
+for :math:`1 \le k \le N_{\rm d}^{(2)}, 1 \le m \le N_{\rm d}^{(3)}`. Hence, the gradient
+of the new global descriptors with respect to atom positions is calculated as
+
+.. math::
+
+    \nabla d^{(2*3)}_{km} = d^{(3)}_m \nabla d^{(2)}_{k}  +  d^{(2)}_{k} \nabla d^{(3)}_m, \quad 1 \le k \le N_{\rm d}^{(2)}, 1 \le m \le N_{\rm d}^{(3)} .
+
+The quadratic  POD potential is defined as a linear combination of the
+original and new global descriptors as follows
+
+.. math::
+
+    E^{(2*3)} = \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} \sum_{m=1}^{N_{\rm 3d}^{(2*3)}} c^{(2*3)}_{km} d^{(2*3)}_{km} .
+
+It thus follows that
+
+.. math::
+
+    E^{(2*3)} = 0.5 \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} \left( \sum_{m=1}^{N_{\rm 3d}^{(2*3)}} c^{(2*3)}_{km} d_m^{(3)} \right) d_k^{(2)} + 0.5 \sum_{m=1}^{N_{\rm 3d}^{(2*3)}} \left( \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} c^{(2*3)}_{km} d_k^{(2)} \right) d_m^{(3)}  ,
+
+which is simplified to
+
+.. math::
+
+    E^{(2*3)} =  0.5 \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} b_k^{(2)} d_k^{(2)} +  0.5 \sum_{m=1}^{N_{\rm 3d}^{(2*3)}}   b_m^{(3)} d_m^{(3)}
+
+where
+
+.. math::
+
+    b_k^{(2)} & = \sum_{m=1}^{N_{\rm 3d}^{(2*3)}} c^{(2*3)}_{km} d_m^{(3)}, \quad k = 1,\ldots, N_{\rm 2d}^{(2*3)}, \\
+    b_m^{(3)} & = \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} c^{(2*3)}_{km} d_k^{(2)}, \quad m = 1,\ldots, N_{\rm 3d}^{(2*3)} .
+
+The quadratic POD potential results in the following atomic forces
+
+.. math::
+
+    \boldsymbol F^{(2*3)} = - \sum_{k=1}^{N_{\rm 2d}^{(2*3)}} \sum_{m=1}^{N_{\rm 3d}^{(2*3)}} c^{(2*3)}_{km}  \nabla d^{(2*3)}_{km} .
+
+It can be shown that
+
+.. math::
+
+    \boldsymbol F^{(2*3)} = - \sum_{k=1}^{N_{\rm 2d}^{(2*3)}}   b^{(2)}_k \nabla d_k^{(2)} - \sum_{m=1}^{N_{\rm 3d}^{(2*3)}}  b^{(3)}_m  \nabla d_m^{(3)} .
+
+The calculation of the atomic forces for the quadratic POD  potential
+only requires the extra calculation of :math:`b_k^{(2)}` and :math:`b_m^{(3)}` which can be negligible.
+As a result, the quadratic  POD potential does not increase the computational complexity.
+
+
+Training
+""""""""
+
+POD potentials are trained using the least-squares regression against
+density functional theory (DFT) data.  Let :math:`J` be the number of
+training configurations, with :math:`N_j` being the number of atoms in
+the j-th configuration. Let :math:`\{E^{\star}_j\}_{j=1}^{J}` and
+:math:`\{\boldsymbol F^{\star}_j\}_{j=1}^{J}` be the DFT energies and
+forces for :math:`J` configurations. Next, we calculate the global
+descriptors and their derivatives for all training configurations. Let
+:math:`d_{jm}, 1 \le m \le M`, be the global descriptors associated with
+the j-th configuration, where :math:`M` is the number of global
+descriptors. We then form a matrix :math:`\boldsymbol A \in
+\mathbb{R}^{J \times M}` with entries :math:`A_{jm} = d_{jm}/ N_j` for
+:math:`j=1,\ldots,J` and :math:`m=1,\ldots,M`.  Moreover, we form a
+matrix :math:`\boldsymbol B \in \mathbb{R}^{\mathcal{N} \times M}` by
+stacking the derivatives of the global descriptors for all training
+configurations from top to bottom, where :math:`\mathcal{N} =
+3\sum_{j=1}^{J} N_j`.
+
+The coefficient vector :math:`\boldsymbol c` of the POD potential is
+found by solving the following least-squares problem
+
+.. math::
+
+    {\min}_{\boldsymbol c \in \mathbb{R}^{M}} \ w_E \|\boldsymbol A(\boldsymbol \eta) \boldsymbol c - \bar{\boldsymbol E}^{\star} \|^2 + w_F \|\boldsymbol B(\boldsymbol \eta) \boldsymbol c + \boldsymbol F^{\star} \|^2 + w_R \|\boldsymbol c \|^2,
+
+where :math:`w_E` and :math:`w_F` are weights for the energy
+(*fitting_weight_energy*) and force (*fitting_weight_force*),
+respectively; and :math:`w_R` is the regularization parameter (*fitting_regularization_parameter*).  Here :math:`\bar{\boldsymbol E}^{\star} \in
+\mathbb{R}^{J}` is a vector of with entries :math:`\bar{E}^{\star}_j =
+E^{\star}_j/N_j` and :math:`\boldsymbol F^{\star}` is a vector of
+:math:`\mathcal{N}` entries obtained by stacking :math:`\{\boldsymbol
+F^{\star}_j\}_{j=1}^{J}` from top to bottom.
+
+The training procedure is the same for both the linear and quadratic POD
+potentials.  However, since the quadratic POD potential has a
+significantly larger number of the global descriptors, it is more
+expensive to train the linear POD potential. This is because the
+training of the quadratic POD potential still requires us to calculate
+and store the quadratic global descriptors and their
+gradient. Furthermore, the quadratic POD potential may require more
+training data in order to prevent over-fitting. In order to reduce the
+computational cost of fitting the quadratic POD potential and avoid
+over-fitting, we can use subsets of two-body and three-body PODs for
+constructing the new descriptors.
+
+
+Restrictions
+""""""""""""
+
+This command is part of the ML-POD package.  It is only enabled if
+LAMMPS was built with that package. See the :doc:`Build package
+<Build_package>` page for more info.
+
+Related commands
+""""""""""""""""
+
+:doc:`pair_style pod <pair_pod>`
+
+Default
+"""""""
+
+The keyword defaults are also given in the description of the input files.
+
+----------
+
+.. _Grepl20072:
+
+**(Grepl)** Grepl, Maday, Nguyen, and Patera, ESAIM: Mathematical Modelling and Numerical Analysis 41(3), 575-605, (2007).
+
+.. _Nguyen20222:
+
+**(Nguyen)** Nguyen and Rohskopf, arXiv preprint arXiv:2209.02362 (2022).

doc/src/fix.rst

@@ -216,9 +216,9 @@ accelerated styles exist.
 * :doc:`edpd/source <fix_dpd_source>` - add heat source to eDPD simulations
 * :doc:`efield <fix_efield>` - impose electric field on system
 * :doc:`ehex <fix_ehex>` - enhanced heat exchange algorithm
-* :doc:`electrode/conp <fix_electrode_conp>` - impose electric potential
+* :doc:`electrode/conp <fix_electrode>` - impose electric potential
-* :doc:`electrode/conq <fix_electrode_conp>` - impose total electric charge
+* :doc:`electrode/conq <fix_electrode>` - impose total electric charge
-* :doc:`electrode/thermo <fix_electrode_conp>` - apply thermo-potentiostat
+* :doc:`electrode/thermo <fix_electrode>` - apply thermo-potentiostat
 * :doc:`electron/stopping <fix_electron_stopping>` - electronic stopping power as a friction force
 * :doc:`electron/stopping/fit <fix_electron_stopping>` - electronic stopping power as a friction force
 * :doc:`enforce2d <fix_enforce2d>` - zero out *z*-dimension velocity and force
@@ -360,6 +360,7 @@ accelerated styles exist.
 * :doc:`saed/vtk <fix_saed_vtk>` - time-average the intensities from :doc:`compute saed <compute_saed>`
 * :doc:`setforce <fix_setforce>` - set the force on each atom
 * :doc:`setforce/spin <fix_setforce>` - set magnetic precession vectors on each atom
+* :doc:`sgcmc <fix_sgcmc>` - fix for hybrid semi-grand canonical MD/MC simulations
 * :doc:`shake <fix_shake>` - SHAKE constraints on bonds and/or angles
 * :doc:`shardlow <fix_shardlow>` - integration of DPD equations of motion using the Shardlow splitting
 * :doc:`smd <fix_smd>` - applied a steered MD force to a group

doc/src/fix_ave_spatial.rst

@@ -0,0 +1,9 @@
+Fix ave/spatial command
+=======================
+
+.. meta::
+   :http-equiv=Refresh: 5; url='https://docs.lammps.org/Commands_removed.html#fix-ave-spatial-and-fix-ave-spatial-sphere'
+
+.. deprecated:: 11Dec2015
+
+The `fix ave/spatial` command has been superseded by :doc:`fix ave/chunk <fix_ave_chunk>`.

doc/src/fix_ave_spatial_sphere.rst

@@ -0,0 +1,9 @@
+Fix ave/spatial command
+=======================
+
+.. meta::
+   :http-equiv=Refresh: 5; url='https://docs.lammps.org/Commands_removed.html#fix-ave-spatial-and-fix-ave-spatial-sphere'
+
+.. deprecated:: 11Dec2015
+
+The `fix ave/spatial/sphere` command has been superseded by :doc:`fix ave/chunk <fix_ave_chunk>`.

doc/src/fix_bond_react.rst

@@ -42,13 +42,16 @@ Syntax
 * template-ID(post-reacted) = ID of a molecule template containing post-reaction topology
 * map_file = name of file specifying corresponding atom-IDs in the pre- and post-reacted templates
 * zero or more individual keyword/value pairs may be appended to each react argument
-* individual_keyword = *prob* or *max_rxn* or *stabilize_steps* or *custom_charges* or *molecule* or *modify_create*
+* individual_keyword = *prob* or *rate_limit* or *max_rxn* or *stabilize_steps* or *custom_charges* or *rescale_charges* or *molecule* or *modify_create*
 
   .. parsed-literal::
 
          *prob* values = fraction seed
            fraction = initiate reaction with this probability if otherwise eligible
            seed = random number seed (positive integer)
+         *rate_limit* = Nlimit Nsteps
+           Nlimit = maximum number of reactions allowed to occur within interval
+           Nsteps = the interval (number of timesteps) over which to count reactions
          *max_rxn* value = N
            N = maximum number of reactions allowed to occur
          *stabilize_steps* value = timesteps
@@ -56,6 +59,9 @@ Syntax
          *custom_charges* value = *no* or fragment-ID
            *no* = update all atomic charges (default)
            fragment-ID = ID of molecule fragment whose charges are updated
+         *rescale_charges* value = *no* or *yes*
+           *no* = do not rescale atomic charges (default)
+           *yes* = rescale charges such that total charge does not change during reaction
          *molecule* value = *off* or *inter* or *intra*
            *off* = allow both inter- and intramolecular reactions (default)
            *inter* = search for reactions between molecules with different IDs
@@ -171,12 +177,12 @@ due to the internal dynamic grouping performed by fix bond/react.
    If the group-ID is an existing static group, react-group-IDs
    should also be specified as this static group or a subset.
 
-The *reset_mol_ids* keyword invokes the :doc:`reset_mol_ids <reset_mol_ids>`
+The *reset_mol_ids* keyword invokes the :doc:`reset_atoms mol
-command after a reaction occurs, to ensure that molecule IDs are
+<reset_atoms>` command after a reaction occurs, to ensure that
-consistent with the new bond topology. The group-ID used for
+molecule IDs are consistent with the new bond topology. The group-ID
-:doc:`reset_mol_ids <reset_mol_ids>` is the group-ID for this fix.
+used for :doc:`reset_atoms mol <reset_atoms>` is the group-ID for this
-Resetting molecule IDs is necessarily a global operation, so it can
+fix.  Resetting molecule IDs is necessarily a global operation, so it
-be slow for very large systems.
+can be slow for very large systems.
 
 The following comments pertain to each *react* argument (in other
 words, they can be customized for each reaction, or reaction step):
@@ -514,28 +520,40 @@ example, the molecule fragment could consist of only the backbone
 atoms of a polymer chain. This constraint can be used to enforce a
 specific relative position and orientation between reacting molecules.
 
+.. versionchanged:: TBD
+
 The constraint of type "custom" has the following syntax:
 
 .. parsed-literal::
 
    custom *varstring*
 
-where "custom" is the required keyword, and *varstring* is a
+where 'custom' is the required keyword, and *varstring* is a variable
-variable expression. The expression must be a valid equal-style
+expression. The expression must be a valid equal-style variable
-variable formula that can be read by the :doc:`variable <variable>` command,
+formula that can be read by the :doc:`variable <variable>` command,
 after any special reaction functions are evaluated. If the resulting
 expression is zero, the reaction is prevented from occurring;
-otherwise, it is permitted to occur. There are two special reaction
+otherwise, it is permitted to occur. There are three special reaction
-functions available, "rxnsum" and "rxnave". These functions operate
+functions available, 'rxnbond', 'rxnsum', and 'rxnave'. The 'rxnbond'
-over the atoms in a given reaction site, and have one mandatory
+function allows per-bond values to be included in the variable strings
-argument and one optional argument. The mandatory argument is the
+of the custom constraint. The 'rxnbond' function has two mandatory
-identifier for an atom-style variable. The second, optional argument
+arguments. The first argument is the ID of a previously defined
-is the name of a molecule fragment in the pre-reaction template, and
+'compute bond/local' command. This 'compute bond/local' must compute
-can be used to operate over a subset of atoms in the reaction site.
+only one value, e.g. bond force. This value is returned by the
-The "rxnsum" function sums the atom-style variable over the reaction
+'rxnbond' function. The second argument is the name of a molecule
-site, while the "rxnave" returns the average value. For example, a
+fragment in the pre-reaction template. The fragment must contain
-constraint on the total potential energy of atoms involved in the
+exactly two atoms, corresponding to the atoms involved in the bond
-reaction can be imposed as follows:
+whose value should be calculated. An example of a constraint that uses
+the force experienced by a bond is provided below. The 'rxnsum' and
+'rxnave' functions operate over the atoms in a given reaction site,
+and have one mandatory argument and one optional argument. The
+mandatory argument is the identifier for an atom-style variable. The
+second, optional argument is the name of a molecule fragment in the
+pre-reaction template, and can be used to operate over a subset of
+atoms in the reaction site. The 'rxnsum' function sums the atom-style
+variable over the reaction site, while the 'rxnave' returns the
+average value. For example, a constraint on the total potential energy
+of atoms involved in the reaction can be imposed as follows:
 
 .. code-block:: LAMMPS
 
@@ -547,11 +565,32 @@ reaction can be imposed as follows:
    custom "rxnsum(v_my_pe) > 100" # in Constraints section of map file
 
 The above example prevents the reaction from occurring unless the
-total potential energy of the reaction site is above 100. The variable
+total potential energy of the reaction site is above 100. As a second
-expression can be interpreted as the probability of the reaction
+example, this time using the 'rxnbond' function, consider a modified
-occurring by using an inequality and the :doc:`random(x,y,z) <variable>`
+Arrhenius constraint that depends on the bond force of a specific bond:
-function available for equal-style variables, similar to the 'arrhenius'
+
-constraint above.
+.. code-block:: LAMMPS
+
+   # in LAMMPS input script
+
+   compute bondforce all bond/local force
+
+   compute ke_atom all ke/atom
+   variable ke atom c_ke_atom
+
+   variable E_a equal 100.0 # activation energy
+   variable l0 equal 1.0 # characteristic length
+
+
+.. code-block:: LAMMPS
+
+   # in Constraints section of map file
+
+   custom "exp(-(v_E_a-rxnbond(c_bondforce,bond1frag)*v_l0)/(2/3*rxnave(v_ke))) < random(0,1,12345)"
+
+By using an inequality and the 'random(x,y,z)' function, the left-hand
+side can be interpreted as the probability of the reaction occurring,
+similar to the 'arrhenius' constraint above.
 
 By default, all constraints must be satisfied for the reaction to
 occur. In other words, constraints are evaluated as a series of
@@ -598,6 +637,15 @@ eligible reaction only occurs if the random number is less than the
 fraction. Up to :math:`N` reactions are permitted to occur, as optionally
 specified by the *max_rxn* keyword.
 
+.. versionadded:: TBD
+
+The *rate_limit* keyword can enforce an upper limit on the overall
+rate of the reaction. The number of reaction occurrences is limited to
+Nlimit within an interval of Nsteps timesteps. No reactions are
+permitted to occur within the first Nsteps timesteps of the first run
+after reading a data file. Nlimit can be specified with an equal-style
+:doc:`variable <variable>`.
+
 The *stabilize_steps* keyword allows for the specification of how many
 time steps a reaction site is stabilized before being returned to the
 overall system thermostat. In order to produce the most physical
@@ -616,6 +664,19 @@ charges are updated to those specified by the post-reaction template
 fragment defined in the pre-reaction molecule template. In this case,
 only the atomic charges of atoms in the molecule fragment are updated.
 
+.. versionadded:: TBD
+
+The *rescale_charges* keyword can be used to ensure the total charge
+of the system does not change as reactions occur. When the argument is
+set to *yes*\ , a fixed value is added to the charges of post-reaction
+atoms such that their total charge equals that of the pre-reaction
+site. If only a subset of atomic charges are updated via the
+*custom_charges* keyword, this rescaling is applied to the subset.
+This keyword could be useful for systems that contain different
+molecules with the same reactive site, if the partial charges on the
+reaction site vary from molecule to molecule, or when removing
+reaction by-products.
+
 The *molecule* keyword can be used to force the reaction to be
 intermolecular, intramolecular or either. When the value is set to
 *off*\ , molecule IDs are not considered when searching for reactions

doc/src/fix_electrode.rst

@@ -0,0 +1,421 @@
+.. index:: fix electrode/conp
+.. index:: fix electrode/conq
+.. index:: fix electrode/thermo
+.. index:: fix electrode/conp/intel
+.. index:: fix electrode/conq/intel
+.. index:: fix electrode/thermo/intel
+
+fix electrode/conp command
+==========================
+
+Accelerator Variant: *electrode/conp/intel*
+
+fix electrode/conq command
+==========================
+
+Accelerator Variant: *electrode/conq/intel*
+
+fix electrode/thermo command
+============================
+
+Accelerator Variant: *electrode/thermo/intel*
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   fix ID group-ID style args keyword value ...
+
+* ID, group-ID are documented in :doc:`fix <fix>` command
+* style = *electrode/conp* or *electrode/conq* or *electrode/thermo*
+* args = arguments used by a particular style
+
+  .. parsed-literal::
+
+       *electrode/conp* args = potential eta
+       *electrode/conq* args = charge eta
+       *electrode/thermo* args = potential eta *temp* values
+            potential = electrode potential
+            charge = electrode charge
+            eta = reciprocal width of electrode charge smearing
+            *temp* values = T_v tau_v rng_v
+                T_v = temperature of thermo-potentiostat
+                tau_v = time constant of thermo-potentiostat
+                rng_v = integer used to initialize random number generator
+
+* zero or more keyword/value pairs may be appended
+* keyword = *algo* or *symm* or *couple* or *etypes* or *ffield* or *write_mat* or *write_inv* or *read_mat* or *read_inv*
+
+.. parsed-literal::
+
+    *algo* values = *mat_inv* or *mat_cg* tol or *cg* tol
+        specify the algorithm used to compute the electrode charges
+    *symm* value = *on* or *off*
+        turn on/off charge neutrality constraint for the electrodes
+    *couple* values = group-ID val
+        group-ID = group of atoms treated as additional electrode
+        val = electric potential or charge on this electrode
+    *etypes* value = *on* or *off*
+        turn on/off type-based optimized neighbor lists (electrode and electrolyte types may not overlap)
+    *ffield* value = *on* or *off*
+        turn on/off finite-field implementation
+    *write_mat* value = filename
+        filename = file to which to write elastance matrix
+    *write_inv* value = filename
+        filename = file to which to write inverted matrix
+    *read_mat* value = filename
+        filename = file from which to read elastance matrix
+    *read_inv* value = filename
+        filename = file from which to read inverted matrix
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   fix fxconp bot electrode/conp -1.0 1.805 couple top 1.0 couple ref 0.0 write_inv inv.csv symm on
+   fix fxconp electrodes electrode/conq 0.0 1.805 algo cg 1e-5
+   fix fxconp bot electrode/thermo -1.0 1.805 temp 298 100 couple top 1.0
+
+Description
+"""""""""""
+
+The *electrode* fixes implement the constant potential method (CPM)
+(:ref:`Siepmann <Siepmann>`, :ref:`Reed <Reed3>`), and modern variants,
+to accurately model electrified, conductive electrodes. This is
+primarily useful for studying electrode-electrolyte interfaces,
+especially at high potential differences or ionicities, with non-planar
+electrodes such as nanostructures or nanopores, and to study dynamic
+phenomena such as charging or discharging time scales or conductivity or
+ionic diffusivities.
+
+Each *electrode* fix allows users to set additional electrostatic
+relationships between the specified groups which model useful
+electrostatic configurations:
+
+* *electrode/conp* sets potentials or potential differences between electrodes
+
+  *  (resulting in changing electrode total charges)
+
+* *electrode/conq* sets the total charge on each electrode
+
+  *  (resulting in changing electrode potentials)
+
+* *electrode/thermo* sets a thermopotentiostat
+  :ref:`(Deissenbeck)<Deissenbeck>` between two electrodes
+
+  * (resulting in changing charges and potentials with appropriate
+     average potential difference and thermal variance)
+
+The first group-ID provided to each fix specifies the first electrode
+group, and more group(s) are added using the *couple* keyword for each
+additional group.  While *electrode/thermo* only accepts two groups,
+*electrode/conp* and *electrode/conq* accept any number of groups, up to
+LAMMPS's internal restrictions (see Restrictions below). Electrode
+groups must not overlap, i.e.  the fix will issue an error if any
+particle is detected to belong to at least two electrode groups.
+
+CPM involves updating charges on groups of electrode particles, per time
+step, so that the system's total energy is minimized with respect to
+those charges.  From basic electrostatics, this is equivalent to making
+each group conductive, or imposing an equal electrostatic potential on
+every particle in the same group (hence the name CPM).  The charges are
+usually modelled as a Gaussian distribution to make the charge-charge
+interaction matrix invertible (:ref:`Gingrich <Gingrich>`).  The keyword
+*eta* specifies the distribution's width in units of inverse length.
+
+.. versionadded:: TBD
+
+Three algorithms are available to minimize the energy, varying in how
+matrices are pre-calculated before a run to provide computational
+speedup. These algorithms can be selected using the keyword *algo*:
+
+* *algo mat_inv* pre-calculates the capacitance matrix and obtains the
+  charge configuration in one matrix-vector calculation per time step
+
+* *algo mat_cg* pre-calculates the elastance matrix (inverse of
+  capacitance matrix) and obtains the charge configuration using a
+  conjugate gradient solver in multiple matrix-vector calculations per
+  time step
+
+* *algo cg* does not perform any pre-calculation and obtains the charge
+  configuration using a conjugate gradient solver and multiple
+  calculations of the electric potential per time step.
+
+For both *cg* methods, the command must specify the conjugate gradient
+tolerance. *fix electrode/thermo* currently only supports the *mat_inv*
+algorithm.
+
+The keyword *symm* can be set *on* (or *off*) to turn on (or turn off)
+the capacitance matrix constraint that sets total electrode charge to be
+zero.  This has slightly different effects for each *fix electrode*
+variant.  For *fix electrode/conp*, with *symm off*, the potentials
+specified are absolute potentials, but the charge configurations
+satisfying them may add up to an overall non-zero, varying charge for
+the electrodes (and thus the simulation box). With *symm on*, the total
+charge over all electrode groups is constrained to zero, and potential
+differences rather than absolute potentials are the physically relevant
+quantities.
+
+For *fix electrode/conq*, with *symm off*, overall neutrality is
+explicitly obeyed or violated by the user input (which is not
+checked!). With *symm on*, overall neutrality is ensured by ignoring the
+user-input charge for the last listed electrode (instead, its charge
+will always be minus the total sum of all other electrode charges). For
+*fix electrode/thermo*, overall neutrality is always automatically
+imposed for any setting of *symm*, but *symm on* allows finite-field
+mode (*ffield on*, described below) for faster simulations.
+
+For all three fixes, any potential (or charge for *conq*) can be
+specified as an equal-style variable prefixed with "v\_". For example,
+the following code will ramp the potential difference between electrodes
+from 0.0V to 2.0V over the course of the simulation:
+
+.. code-block:: LAMMPS
+
+   fix fxconp bot electrode/conp 0.0 1.805 couple top v_v symm on
+   variable v equal ramp(0.0, 2.0)
+
+Note that these fixes only parse their supplied variable name when
+starting a run, and so these fixes will accept equal-style variables
+defined *after* the fix definition, including variables dependent on the
+fix's own output. This is useful, for example, in the fix's internal
+finite-field commands (see below).  For an advanced example of this see
+the in.conq2 input file in the directory
+``examples/PACKAGES/electrode/graph-il``.
+
+This fix necessitates the use of a long range solver that calculates and
+provides the matrix of electrode-electrode interactions and a vector of
+electrode-electrolyte interactions.  The Kspace styles
+*ewald/electrode*, *pppm/electrode* and *pppm/electrode/intel* are
+created specifically for this task :ref:`(Ahrens-Iwers) <Ahrens-Iwers>`.
+
+For systems with non-periodic boundaries in one or two directions dipole
+corrections are available with the :doc:`kspace_modify <kspace_modify>`.
+For ewald/electrode a two-dimensional Ewald summation :ref:`(Hu) <Hu>`
+can be used by setting "slab ew2d":
+
+.. code-block:: LAMMPS
+
+   kspace_modify slab <slab_factor>
+   kspace_modify wire <wire_factor>
+   kspace_modify slab ew2d
+
+Two implementations for the calculation of the elastance matrix are
+available with pppm and can be selected using the *amat onestep/twostep*
+keyword.  *onestep* is the default; *twostep* can be faster for large
+electrodes and a moderate mesh size but requires more memory.
+
+.. code-block:: LAMMPS
+
+   kspace_modify amat onestep/twostep
+
+For all versions of the fix, the keyword-value *ffield on* enables the
+finite-field mode (:ref:`Dufils <Dufils>`, :ref:`Tee <Tee>`), which uses
+an electric field across a periodic cell instead of non-periodic
+boundary conditions to impose a potential difference between the two
+electrodes bounding the cell. The fix (with name *fix-ID*) detects which
+of the two electrodes is "on top" (has the larger maximum *z*-coordinate
+among all particles).  Assuming the first electrode group is on top, it
+then issues the following commands internally:
+
+.. code-block:: LAMMPS
+
+   variable fix-ID_ffield_zfield equal (f_fix-ID[2]-f_fix-ID[1])/lz
+   efield fix-ID_efield all efield 0.0 0.0 v_fix-ID_ffield_zfield
+
+which implements the required electric field as the potential difference
+divided by cell length.  The internal commands use variable so that the
+electric field will correctly vary with changing potentials in the
+correct way (for example with equal-style potential difference or with
+*fix electrode/conq*).  This keyword requires two electrodes and will
+issue an error with any other number of electrodes. This keyword
+requires electroneutrality to be imposed (*symm on*) and will issue an
+error otherwise.
+
+.. versionchanged:: TBD
+
+For all versions of the fix, the keyword-value *etypes on* enables
+type-based optimized neighbor lists. With this feature enabled, LAMMPS
+provides the fix with an occasional neighbor list restricted to
+electrode-electrode interactions for calculating the electrode matrix,
+and a perpetual neighbor list restricted to electrode-electrolyte
+interactions for calculating the electrode potentials, using particle
+types to list only desired interactions, and typically resulting in
+5--10\% less computational time.  Without this feature the fix will
+simply use the active pair style's neighbor list.  This feature cannot
+be enabled if any electrode particle has the same type as any
+electrolyte particle (which would be unusual in a typical simulation)
+and the fix will issue an error in that case.
+
+Restart, fix_modify, output, run start/stop, minimize info
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+This fix currently does not write any information to restart files.
+
+The *fix_modify tf* option enables the Thomas-Fermi metallicity model
+(:ref:`Scalfi <Scalfi>`) and allows parameters to be set for each atom type.
+
+.. code-block:: LAMMPS
+
+   fix_modify ID tf type length voronoi
+
+
+If this option is used parameters must be set for all atom types of the
+electrode.
+
+The *fix_modify timer* option turns on (off) additional timer outputs in the log
+file, for code developers to track optimization.
+
+.. code-block:: LAMMPS
+
+   fix_modify ID timer on/off
+
+----------
+
+These fixes compute a global (extensive) scalar, a global (intensive)
+vector, and a global array, which can be accessed by various
+:doc:`output commands <Howto_output>`.
+
+The global scalar outputs the energy added to the system by this fix,
+which is the negative of the total charge on each electrode multiplied
+by that electrode's potential.
+
+The global vector outputs the potential on each electrode (and thus has
+*N* entries if the fix manages *N* electrode groups), in :doc:`units
+<units>` of electric field multiplied by distance (thus volts for *real*
+and *metal* units).  The electrode groups' ordering follows the order in
+which they were input in the fix command using *couple*. The global
+vector output is useful for *fix electrode/conq* and *fix
+electrode/thermo*, where potential is dynamically updated based on
+electrolyte configuration instead of being directly set.
+
+The global array has *N* rows and *2N+1* columns, where the fix manages
+*N* electrode groups managed by the fix. For the *I*-th row of the
+array, the elements are:
+
+* array[I][1] = total charge that group *I* would have had *if it were
+  at 0 V applied potential* * array[I][2 to *N* + 1] = the *N* entries
+  of the *I*-th row of the electrode capacitance matrix (definition
+  follows) * array[I][*N* + 2 to *2N* + 1] = the *N* entries of the
+  *I*-th row of the electrode elastance matrix (the inverse of the
+  electrode capacitance matrix)
+
+The :math:`N \times N` electrode capacitance matrix, denoted :math:`\mathbf{C}`
+in the following equation, summarizes how the total charge induced on each
+electrode (:math:`\mathbf{Q}` as an *N*-vector) is related to the potential on
+each electrode, :math:`\mathbf{V}`, and the charge-at-0V :math:`\mathbf{Q}_{0V}`
+(which is influenced by the local electrolyte structure):
+
+.. math::
+
+   \mathbf{Q} = \mathbf{Q}_{0V} + \mathbf{C} \cdot \mathbf{V}
+
+The charge-at-0V, electrode capacitance and elastance matrices are internally
+used to calculate the potentials required to induce the specified total
+electrode charges in *fix electrode/conq* and *fix electrode/thermo*. With the
+*symm on* option, the electrode capacitance matrix would be singular, and thus
+its last row is replaced with *N* copies of its top-left entry
+(:math:`\mathbf{C}_{11}`) for invertibility.
+
+The global array output is mainly useful for quickly determining the 'vacuum
+capacitance' of the system (capacitance with only electrodes, no electrolyte),
+and can also be used for advanced simulations setting the potential as some
+function of the charge-at-0V (such as the ``in.conq2`` example mentioned above).
+
+Please cite :ref:`(Ahrens-Iwers2022) <Ahrens-Iwers2>` in any publication that
+uses this implementation.  Please cite also the publication on the combination
+of the CPM with PPPM if you use *pppm/electrode* :ref:`(Ahrens-Iwers)
+<Ahrens-Iwers>`.
+
+----------
+
+Restrictions
+""""""""""""
+
+For algorithms that use a matrix for the electrode-electrode
+interactions, positions of electrode particles have to be immobilized at
+all times.
+
+With *ffield off* (i.e. the default), the box geometry is expected to be
+*z*-non-periodic (i.e. *boundary p p f*), and this fix will issue an
+error if the box is *z*-periodic. With *ffield on*, the box geometry is
+expected to be *z*-periodic, and this fix will issue an error if the box
+is *z*-non-periodic.
+
+The parallelization for the fix works best if electrode atoms are evenly
+distributed across processors. For a system with two electrodes at the bottom
+and top of the cell this can be achieved with *processors * * 2*, or with the
+line
+
+.. code-block:: LAMMPS
+
+   if "$(extract_setting(world_size) % 2) == 0" then "processors * * 2"
+
+which avoids an error if the script is run on an odd number of
+processors (such as on just one processor for testing).
+
+The fix creates an additional group named *[fix-ID]_group* which is the
+union of all electrode groups supplied to LAMMPS. This additional group
+counts towards LAMMPS's limitation on the total number of groups
+(currently 32), which may not allow scripts that use that many groups to
+run with this fix.
+
+The matrix-based algorithms (*algo mat_inv* and *algo mat_cg*) currently
+store an interaction matrix (either elastance or capacitance) of *N* by
+*N* doubles for each MPI process. This memory requirement may be
+prohibitive for large electrode groups.  The fix will issue a warning if
+it expects to use more than 0.5 GiB of memory.
+
+Default
+"""""""
+
+The default keyword-option settings are *algo mat_inv*, *symm off*,
+*etypes off* and *ffield off*.
+
+----------
+
+.. include:: accel_styles.rst
+
+----------
+
+.. _Siepmann:
+
+**(Siepmann)** Siepmann and Sprik, J. Chem. Phys. 102, 511 (1995).
+
+.. _Reed3:
+
+**(Reed)** Reed *et al.*, J. Chem. Phys. 126, 084704 (2007).
+
+.. _Deissenbeck:
+
+**(Deissenbeck)** Deissenbeck *et al.*, Phys. Rev. Letters 126, 136803 (2021).
+
+.. _Gingrich:
+
+**(Gingrich)** Gingrich, `MSc thesis` <https://gingrich.chem.northwestern.edu/papers/ThesiswCorrections.pdf>` (2010).
+
+.. _Ahrens-Iwers:
+
+**(Ahrens-Iwers)** Ahrens-Iwers and Meissner, J. Chem. Phys. 155, 104104 (2021).
+
+.. _Hu:
+
+**(Hu)** Hu, J. Chem. Theory Comput. 10, 5254 (2014).
+
+.. _Dufils:
+
+**(Dufils)** Dufils *et al.*, Phys. Rev. Letters 123, 195501 (2019).
+
+.. _Tee:
+
+**(Tee)** Tee and Searles, J. Chem. Phys. 156, 184101 (2022).
+
+.. _Scalfi:
+
+**(Scalfi)** Scalfi *et al.*, J. Chem. Phys., 153, 174704 (2020).
+
+.. _Ahrens-Iwers2:
+
+**(Ahrens-Iwers2022)** Ahrens-Iwers *et al.*, J. Chem. Phys. 157, 084801 (2022).

doc/src/fix_electrode_conp.rst

@@ -1,230 +0,0 @@
-.. index:: fix electrode/conp
-.. index:: fix electrode/conq
-.. index:: fix electrode/thermo
-.. index:: fix electrode/conp/intel
-.. index:: fix electrode/conq/intel
-.. index:: fix electrode/thermo/intel
-
-fix electrode/conp command
-==========================
-
-Accelerator Variant: *electrode/conp/intel*
-
-fix electrode/conq command
-==========================
-
-Accelerator Variant: *electrode/conq/intel*
-
-fix electrode/thermo command
-============================
-
-Accelerator Variant: *electrode/thermo/intel*
-
-Syntax
-""""""
-
-.. parsed-literal::
-
-   fix ID group-ID style args keyword value ...
-
-* ID, group-ID are documented in :doc:`fix <fix>` command
-* style = *electrode/conp* or *electrode/conq* or *electrode/thermo*
-* args = arguments used by a particular style
-
-  .. parsed-literal::
-
-       *electrode/conp* args = potential eta
-       *electrode/conq* args = charge eta
-       *electrode/thermo* args = potential eta *temp* values
-            potential = electrode potential
-            charge = electrode charge
-            eta = reciprocal width of electrode charge smearing
-            *temp* values = T_v tau_v rng_v
-                T_v = temperature of thermo-potentiostat
-                tau_v = time constant of thermo-potentiostat
-                rng_v = integer used to initialize random number generator
-
-* zero or more keyword/value pairs may be appended
-* keyword = *symm* or *couple* or *etypes* or *ffield* or *write_mat* or *write_inv* or *read_mat* or *read_inv*
-
-.. parsed-literal::
-
-    *symm* value = *on* or *off*
-        turn on/off charge neutrality constraint for the electrodes
-    *couple* values = group-ID val
-        group-ID = group of atoms treated as additional electrode
-        val = electric potential or charge on this electrode
-    *etypes* values = type
-        type = atom type (can be a range) exclusive to the electrode for optimized neighbor lists
-    *ffield* value = *on* or *off*
-        turn on/off finite-field implementation
-    *write_mat* value = filename
-        filename = file to which to write elastance matrix
-    *write_inv* value = filename
-        filename = file to which to write inverted matrix
-    *read_mat* value = filename
-        filename = file from which to read elastance matrix
-    *read_inv* value = filename
-        filename = file from which to read inverted matrix
-
-Examples
-""""""""
-
-.. code-block:: LAMMPS
-
-   fix fxconp bot electrode/conp -1.0 1.805 couple top 1.0 couple ref 0.0 write_inv inv.csv symm on
-   fix fxconp electrodes electrode/conq 0.0 1.805
-   fix fxconp bot electrode/thermo -1.0 1.805 temp 298 100 couple top 1.0
-
-Description
-"""""""""""
-
-fix electrode/conp mode implements a constant potential method (CPM)
-(:ref:`Siepmann <Siepmann>`, :ref:`Reed <Reed3>`). Charges of groups specified
-via group-ID and optionally with the `couple` keyword are adapted to meet their respective
-potential at every time step.  An arbitrary number of electrodes can be set but
-the respective groups may not overlap. Electrode charges have a Gaussian charge
-distribution with reciprocal width eta. The energy minimization is achieved via
-matrix inversion :ref:`(Wang) <Wang5>`.
-
-fix electrode/conq enforces a total charge specified in the input on each electrode. The energy is
-minimized w.r.t. the charge distribution within the electrode.
-
-fix electrode/thermo implements a thermo-potentiostat :ref:`(Deissenbeck)
-<Deissenbeck>`. Temperature and time constant of the thermo-potentiostat need
-to be specified using the temp keyword. Currently, only two electrodes are possible with
-this style.
-
-This fix necessitates the use of a long range solver that calculates and provides the matrix
-of electrode-electrode interactions and a vector of electrode-electrolyte
-interactions.  The Kspace styles *ewald/electrode*, *pppm/electrode* and
-*pppm/electrode/intel* are created specifically for this task
-:ref:`(Ahrens-Iwers) <Ahrens-Iwers>`.
-
-For systems with non-periodic boundaries in one or two directions dipole
-corrections are available with the :doc:`kspace_modify <kspace_modify>`.  For
-ewald/electrode a two-dimensional Ewald summation :ref:`(Hu) <Hu>` can be used
-by setting "slab ew2d":
-
-.. code-block:: LAMMPS
-
-   kspace_modify slab <slab_factor>
-   kspace_modify wire <wire_factor>
-   kspace_modify slab ew2d
-
-Two implementations for the calculation of the elastance matrix are available
-with pppm and can be selected using the *amat onestep/twostep* keyword.
-*onestep* is the default; *twostep* can be faster for large electrodes and a
-moderate mesh size but requires more memory.
-
-.. code-block:: LAMMPS
-
-   kspace_modify amat onestep/twostep
-
-
-The *fix_modify tf* option enables the Thomas-Fermi metallicity model
-(:ref:`Scalfi <Scalfi>`) and allows parameters to be set for each atom type.
-
-.. code-block:: LAMMPS
-
-   fix_modify ID tf type length voronoi
-
-
-If this option is used parameters must be set for all atom types of the electrode.
-
-The *fix_modify timer* option turns on (off) additional timer outputs in the log
-file, for code developers to track optimization.
-
-.. code-block:: LAMMPS
-
-   fix_modify ID timer on/off
-
-The *fix_modify set* options allow calculated quantities to be accessed via
-internal variables. Currently four types of quantities can be accessed:
-
-.. code-block:: LAMMPS
-
-   fix-modify ID set v group-ID variablename
-   fix-modify ID set qsb group-ID variablename
-   fix-modify ID set mc group-ID1 group-ID2 variablename
-   fix-modify ID set me group-ID1 group-ID2 variablename
-
-One use case is to output the potential that is internally calculated and
-applied to each electrode group by *fix electrode/conq* or *fix electrode/thermo*.
-For that case the *v* option makes *fix electrode* update the variable
-*variablename* with the potential applied to group *group-ID*, where *group-ID*
-must be a group whose charges are updated by *fix electrode* and *variablename*
-must be an internal-style variable:
-
-.. code-block:: LAMMPS
-
-   fix conq bot electrode/conq -1.0 1.979 couple top 1.0
-   variable vbot internal 0.0
-   fix_modify conq set v bot vbot
-
-The *qsb* option similarly outputs the total updated charge of the group if its
-potential were 0.0V. The *mc* option requires two *group-IDs*, and outputs the
-entry \{*group-ID1*, *group-ID2*\} of the (symmetric) *macro-capacitance* matrix
-(MC) which relates the electrodes' applied potentials (V), total charges (Q), and
-total charges at 0.0 V (Qsb):
-
-.. math::
-
-   \mathbf{Q} = \mathbf{Q}_{SB} + \mathbf{MC} \cdot \mathbf{V}
-
-Lastly, the *me* option also requires two *group-IDs* and outputs the entry
-\{*group-ID1*, *group-ID2*\} of the *macro-elastance* matrix, which is the
-inverse of the macro-capacitance matrix. (As the names denote, the
-macro-capacitance matrix gives electrode charges from potentials, and the
-macro-elastance matrix gives electrode potentials from charges).
-
-.. warning::
-
-   Positions of electrode particles have to be immobilized at all times.
-
-The parallelization for the fix works best if electrode atoms are evenly
-distributed across processors. For a system with two electrodes at the bottom
-and top of the cell this can be achieved with *processors * * 2*, or with the
-line
-
-.. code-block:: LAMMPS
-
-   if "$(extract_setting(world_size) % 2) == 0" then "processors * * 2"
-
-which avoids an error if the script is run on an odd number of processors (such
-as on just one processor for testing).
-
-----------
-
-.. include:: accel_styles.rst
-
-----------
-
-.. _Siepmann:
-
-**(Siepmann)** Siepmann and Sprik, J. Chem. Phys. 102, 511 (1995).
-
-.. _Reed3:
-
-**(Reed)** Reed *et al.*, J. Chem. Phys. 126, 084704 (2007).
-
-.. _Wang5:
-
-**(Wang)** Wang *et al.*, J. Chem. Phys. 141, 184102 (2014).
-
-.. _Deissenbeck:
-
-**(Deissenbeck)** Deissenbeck *et al.*, Phys. Rev. Letters 126, 136803 (2021).
-
-.. _Ahrens-Iwers:
-
-**(Ahrens-Iwers)** Ahrens-Iwers and Meissner, J. Chem. Phys. 155, 104104 (2021).
-
-.. _Hu:
-
-**(Hu)** Hu, J. Chem. Theory Comput. 10, 5254 (2014).
-
-.. _Scalfi:
-
-**(Scalfi)** Scalfi *et al.*, J. Chem. Phys., 153, 174704 (2020).
-

doc/src/fix_ipi.rst

@@ -90,6 +90,12 @@ coordinates are transferred. However, one could use this strategy to
 define an external potential acting on the atoms that are moved by
 i-PI.
 
+Since the i-PI code uses atomic units internally, this fix needs to
+convert LAMMPS data to and from its :doc:`specified units <units>`
+accordingly when communicating with i-PI.  This is not possible for
+reduced units ("units lj") and thus *fix ipi* will stop with an error in
+this case.
+
 This fix is part of the MISC package.  It is only enabled if
 LAMMPS was built with that package.  See the
 :doc:`Build package <Build_package>` page for more info.

doc/src/fix_nh_uef.rst

@@ -44,19 +44,23 @@ Examples
 Description
 """""""""""
 
-This fix can be used to simulate non-equilibrium molecular dynamics
+These fixes can be used to simulate non-equilibrium molecular dynamics
-(NEMD) under diagonal flow fields, including uniaxial and bi-axial
+(NEMD) under diagonal flow fields, including uniaxial and bi-axial flow.
-flow.  Simulations under continuous extensional flow may be carried
+Simulations under continuous extensional flow may be carried out for an
-out for an indefinite amount of time.  It is an implementation of the
+indefinite amount of time.  It is an implementation of the boundary
-boundary conditions from :ref:`(Dobson) <Dobson>`, and also uses numerical
+conditions from :ref:`(Dobson) <Dobson>`, and also uses numerical
 lattice reduction as was proposed by :ref:`(Hunt) <Hunt>`. The lattice
-reduction algorithm is from :ref:`(Semaev) <Semaev>`. The fix is intended for
+reduction algorithm is from :ref:`(Semaev) <Semaev>`. The fix is
-simulations of homogeneous flows, and integrates the SLLOD equations
+intended for simulations of homogeneous flows, and integrates the SLLOD
-of motion, originally proposed by Hoover and Ladd (see :ref:`(Evans and Morriss) <Sllod>`).  Additional detail about this implementation can be
+equations of motion, originally proposed by Hoover and Ladd (see
-found in :ref:`(Nicholson and Rutledge) <Nicholson>`.
+:ref:`(Evans and Morriss) <Sllod>`).  Additional detail about this
+implementation can be found in :ref:`(Nicholson and Rutledge)
+<Nicholson>`.
 
 Note that NEMD simulations of a continuously strained system can be
-performed using the :doc:`fix deform <fix_deform>`, :doc:`fix nvt/sllod <fix_nvt_sllod>`, and :doc:`compute temp/deform <compute_temp_deform>` commands.
+performed using the :doc:`fix deform <fix_deform>`, :doc:`fix nvt/sllod
+<fix_nvt_sllod>`, and :doc:`compute temp/deform <compute_temp_deform>`
+commands.
 
 The applied flow field is set by the *eps* keyword. The values
 *edot_x* and *edot_y* correspond to the strain rates in the xx and yy
@@ -73,11 +77,11 @@ to -(*edot_x* + *edot_y*).
 The boundary conditions require a simulation box that does not have a
 consistent alignment relative to the applied flow field. Since LAMMPS
 utilizes an upper-triangular simulation box, it is not possible to
-express the evolving simulation box in the same coordinate system as
+express the evolving simulation box in the same coordinate system as the
-the flow field.  This fix keeps track of two coordinate systems: the
+flow field.  These fixes keep track of two coordinate systems: the flow
-flow frame, and the upper triangular LAMMPS frame. The coordinate
+frame, and the upper triangular LAMMPS frame. The coordinate systems are
-systems are related to each other through the QR decomposition, as is
+related to each other through the QR decomposition, as is illustrated in
-illustrated in the image below.
+the image below.
 
 .. image:: JPG/uef_frames.jpg
    :align: center
@@ -99,12 +103,12 @@ using the dump command will be in the LAMMPS frame unless the
 ----------
 
 Temperature control is achieved with the default Nose-Hoover style
-thermostat documented in :doc:`fix npt <fix_nh>`. When this fix is
+thermostat documented in :doc:`fix nvt <fix_nh>`.  When this fix is
 active, only the peculiar velocity of each atom is stored, defined as
-the velocity relative to the streaming velocity. This is in contrast
+the velocity relative to the streaming velocity. This is in contrast to
-to :doc:`fix nvt/sllod <fix_nvt_sllod>`, which uses a lab-frame
+:doc:`fix nvt/sllod <fix_nvt_sllod>`, which uses a lab-frame velocity,
-velocity, and removes the contribution from the streaming velocity in
+and removes the contribution from the streaming velocity in order to
-order to compute the temperature.
+compute the temperature.
 
 Pressure control is achieved using the default Nose-Hoover barostat
 documented in :doc:`fix npt <fix_nh>`. There are two ways to control the
@@ -156,8 +160,8 @@ The following commands will not work:
 
 ----------
 
-These fix computes a temperature and pressure each timestep.  To do
+These fixes compute a temperature and pressure each timestep.  To do
-this, it creates its own computes of style "temp/uef" and
+this, they create their own computes of style "temp/uef" and
 "pressure/uef", as if one of these two sets of commands had been
 issued:
 
@@ -169,18 +173,19 @@ issued:
    compute fix-ID_temp all temp/uef
    compute fix-ID_press all pressure/uef fix-ID_temp
 
-See the :doc:`compute temp/uef <compute_temp_uef>` and :doc:`compute pressure/uef <compute_pressure_uef>` commands for details.  Note
+See the :doc:`compute temp/uef <compute_temp_uef>` and :doc:`compute
-that the IDs of the new computes are the fix-ID + underscore + "temp"
+pressure/uef <compute_pressure_uef>` commands for details.  Note that
-or fix_ID + underscore + "press".
+the IDs of the new computes are the fix-ID + underscore + "temp" or
+fix_ID + underscore + "press".
 
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
 The fix writes the state of all the thermostat and barostat variables,
-as well as the cumulative strain applied, to :doc:`binary restart files <restart>`.  See the :doc:`read_restart <read_restart>` command
+as well as the cumulative strain applied, to :doc:`binary restart files
-for info on how to re-specify a fix in an input script that reads a
+<restart>`.  See the :doc:`read_restart <read_restart>` command for info
-restart file, so that the operation of the fix continues in an
+on how to re-specify a fix in an input script that reads a restart file,
-uninterrupted fashion.
+so that the operation of the fix continues in an uninterrupted fashion.
 
 .. note::
 
@@ -189,43 +194,41 @@ uninterrupted fashion.
    not contain the cumulative applied strain, will this keyword be
    necessary.
 
-This fix can be used with the :doc:`fix_modify <fix_modify>` *temp* and
+These fixes can be used with the :doc:`fix_modify <fix_modify>` *temp*
-*press* options. The temperature and pressure computes used must be of
+and *press* options. The temperature and pressure computes used must be
-type *temp/uef* and *pressure/uef*\ .
+of type *temp/uef* and *pressure/uef*\ .
 
-This fix computes the same global scalar and vector quantities as :doc:`fix npt <fix_nh>`.
+These fixes compute the same global scalar and vector quantities as
+:doc:`fix nvt andnpt <fix_nh>`.
 
-The fix is not invoked during :doc:`energy minimization <minimize>`.
+These fixes are not invoked during :doc:`energy minimization <minimize>`.
 
 Restrictions
 """"""""""""
 
-This fix is part of the UEF package. It is only enabled if LAMMPS
+These fixes are part of the UEF package. They are only enabled if LAMMPS
-was built with that package. See the :doc:`Build package <Build_package>` page for more info.
+was built with that package. See the :doc:`Build package
+<Build_package>` page for more info.
 
 Due to requirements of the boundary conditions, when the *strain*
 keyword is set to zero (or unset), the initial simulation box must be
 cubic and have style triclinic. If the box is initially of type ortho,
 use :doc:`change_box <change_box>` before invoking the fix.
 
-.. note::
-
-   When resuming from restart files, you may need to use :doc:`box tilt
-   large <box>` since LAMMPS has internal criteria from lattice
-   reduction that are not the same as the criteria in the numerical
-   lattice reduction algorithm.
-
 Related commands
 """"""""""""""""
 
-:doc:`fix nvt <fix_nh>`, :doc:`fix nvt/sllod <fix_nvt_sllod>`, :doc:`compute temp/uef <compute_temp_uef>`, :doc:`compute pressure/uef <compute_pressure_uef>`, :doc:`dump cfg/uef <dump_cfg_uef>`
+:doc:`fix nvt <fix_nh>`, :doc:`fix npt <fix_nh>`, `fix nvt/sllod
+:doc:<fix_nvt_sllod>`, `compute temp/uef <compute_temp_uef>`,
+:doc::doc:`compute pressure/uef <compute_pressure_uef>`, `dump cfg/uef
+:doc:<dump_cfg_uef>`
 
 Default
 """""""
 
-The default keyword values specific to this fix are exy = xyz, strain
+The default keyword values specific to these fixes are exy = xyz, strain
-= 0 0.  The remaining defaults are the same as for :doc:`fix npt <fix_nh>`
+= 0 0.  The remaining defaults are the same as for :doc:`fix nvt or npt
-except tchain = 1.  The reason for this change is given in
+<fix_nh>` except tchain = 1.  The reason for this change is given in
 :doc:`fix nvt/sllod <fix_nvt_sllod>`.
 
 ----------

doc/src/fix_pimd.rst

@@ -156,6 +156,8 @@ This fix is part of the REPLICA package.  It is only enabled if
 LAMMPS was built with that package.  See the
 :doc:`Build package <Build_package>` page for more info.
 
+Fix pid cannot be used with :doc:`lj units <units>`.
+
 A PIMD simulation can be initialized with a single data file read via
 the :doc:`read_data <read_data>` command.  However, this means all
 quasi-beads in a ring polymer will have identical positions and

doc/src/fix_sgcmc.rst

@@ -0,0 +1,188 @@
+.. index:: fix sgcmc
+
+fix sgcmc command
+=================
+
+Syntax
+""""""
+
+.. parsed-literal::
+
+   fix ID group-ID sgcmc every_nsteps swap_fraction temperature deltamu ...
+
+* ID, group-ID are documented in :doc:`fix <fix>` command
+* sgcmc = style name of this fix command
+* every_nsteps = number of MD steps between MC cycles
+* swap_fraction = fraction of a full MC cycle carried out at each call (a value of 1.0 will perform as many trial moves as there are atoms)
+* temperature = temperature that enters Boltzmann factor in Metropolis criterion (usually the same as MD temperature)
+* deltamu = chemical potential difference(s) (`N-1` values must be provided, with `N` being the number of elements)
+* Zero or more keyword/value pairs may be appended to fix definition line:
+
+  .. parsed-literal::
+
+     keyword = *variance* or *randseed* or *window_moves* or *window_size*
+       *variance* kappa conc1 [conc2] ... [concN]
+         kappa = variance constraint parameter
+         conc1,conc2,...  = target concentration(s) in the range 0.0-1.0 (*N-1* values must be provided, with *N* being the number of elements)
+       *randseed* N
+         N = seed for pseudo random number generator
+       *window_moves* N
+         N = number of times sampling window is moved during one MC cycle
+       *window_size* frac
+         frac = size of sampling window (must be between 0.5 and 1.0)
+
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   fix mc all sgcmc 50 0.1 400.0 -0.55
+   fix vc all sgcmc 20 0.2 700.0 -0.7 randseed 324234 variance 2000.0 0.05
+   fix 2  all sgcmc 20 0.1 700.0 -0.7 window_moves 20
+
+Description
+"""""""""""
+
+.. versionadded:: TBD
+
+This command allows to carry out parallel hybrid molecular
+dynamics/Monte Carlo (MD/MC) simulations using the algorithms described
+in :ref:`(Sadigh1) <Sadigh1>`.  Simulations can be carried out in either
+the semi-grand canonical (SGC) or variance constrained semi-grand
+canonical (VC-SGC) ensemble :ref:`(Sadigh2) <Sadigh2>`. Only atom type
+swaps are performed by the SGCMC fix. Relaxations are accounted for by
+the molecular dynamics integration steps.
+
+This fix can be used with standard multi-element EAM potentials
+(:doc:`pair styles eam/alloy or eam/fs <pair_eam>`)
+
+The SGCMC fix can handle Finnis/Sinclair type EAM potentials where
+:math:`\rho(r)` is atom-type specific, such that different elements can
+contribute differently to the total electron density at an atomic site
+depending on the identity of the element at that atomic site.
+
+------------
+
+If this fix is applied, the regular MD simulation will be interrupted in
+defined intervals to carry out a fraction of a Monte Carlo (MC)
+cycle. The interval is set using the parameter *every_nsteps* which
+determines how many MD integrator steps are taken between subsequent
+calls to the MC routine.
+
+It is possible to carry out pure lattice MC simulations by setting
+*every_nsteps* to 1 and not defining an integration fix such as NVE,
+NPT etc.  In that case, the particles will not move and only the MC
+routine will be called to perform atom type swaps.
+
+The parameter *swap_fraction* determines how many MC trial steps are carried
+out every time the MC routine is entered. It is measured in units of full MC
+cycles where one full cycle, *swap_fraction=1*, corresponds to as many MC
+trial steps as there are atoms.
+
+------------
+
+The parameter *temperature* specifies the temperature that is used
+to evaluate the Metropolis acceptance criterion. While it usually
+should be set to the same value as the MD temperature there are cases
+when it can be useful to use two different values for at least part of
+the simulation, e.g., to speed up equilibration at low temperatures.
+
+------------
+
+The parameter *deltamu* is used to set the chemical potential difference
+in the SGC MC algorithm (see Eq. 16 in :ref:`Sadigh1 <Sadigh1>`). By
+convention it is the difference of the chemical potentials of elements
+`B`, `C` ..., with respect to element A. When the simulation includes
+`N` elements, `N-1` values must be specified.
+
+------------
+
+The variance-constrained SGC MC algorithm is activated if the keyword
+*variance* is used. In that case the fix parameter *deltamu* determines
+the effective average constraint in the parallel VC-SGC MC algorithm
+(parameter :math:`\delta\mu_0` in Eq. (20) of :ref:`Sadigh1
+<Sadigh1>`). The parameter *kappa* specifies the variance constraint
+(see Eqs. (20-21) in :ref:`Sadigh1 <Sadigh1>`).
+
+The parameter *conc* sets the target concentration (parameter
+:math:`c_0` in Eqs.  (20-21) of :ref:`Sadigh1 <Sadigh1>`). The atomic
+concentrations refer to components `B`, `C` ..., with `A` being set
+automatically. When the simulation includes `N` elements, `N-1`
+concentration values must be specified.
+
+------------
+
+There are several technical parameters that can be set via optional flags.
+
+*randseed* is expected to be a positive integer number and is used
+to initialize the random number generator on each processor.
+
+*window_size* controls the size of the sampling window in a parallel MC
+simulation. The size has to lie between 0.5 and 1.0. Normally, this
+parameter should be left unspecified which instructs the code to choose
+the optimal window size automatically (see Sect. III.B and Figure 6 in
+:ref:`Sadigh1 <Sadigh1>` for details).
+
+The number of times the window is moved during a MC cycle is set using
+the parameter *window_moves* (see Sect. III.B in :ref:`Sadigh1
+<Sadigh1>` for details).
+
+------------
+
+Restart, fix_modify, output, run start/stop, minimize info
+==========================================================
+
+No information about this fix is written to restart files.
+
+The MC routine keeps track of the global concentration(s) as well as the
+number of accepted and rejected trial swaps during each MC step. These
+values are provided by the sgcmc fix in the form of a global vector that
+can be accessed by various :doc:`output commands <Howto_output>`
+components of the vector represent the following quantities:
+
+* 1 = The absolute number of accepted trial swaps during the last MC step
+* 2 = The absolute number of rejected trial swaps during the last MC step
+* 3 = The current global concentration of species *A* (= number of atoms of type 1 / total number of atoms)
+* 4 = The current global concentration of species *B* (= number of atoms of type 2 / total number of atoms)
+* ...
+* N+2: The current global concentration of species *X* (= number of atoms of type *N* / total number of atoms)
+
+Restrictions
+============
+
+This fix is part of the MC package. It is only enabled if LAMMPS was
+built with that package.  See the :doc:`Build package <Build_package>`
+page for more info.
+
+At present the fix provides optimized subroutines for EAM type
+potentials (see above) that calculate potential energy changes due to
+*local* atom type swaps very efficiently.  Other potentials are
+supported by using the generic potential functions. This, however, will
+lead to exceedingly slow simulations since it implies that the
+energy of the *entire* system is recomputed at each MC trial step.  If
+other potentials are to be used it is strongly recommended to modify and
+optimize the existing generic potential functions for this purpose.
+Also, the generic energy calculation can not be used for parallel
+execution i.e. it only works with a single MPI process.
+
+------------
+
+Default
+=======
+
+The optional parameters default to the following values:
+
+* *randseed* = 324234
+* *window_moves* = 8
+* *window_size* = automatic
+
+------------
+
+.. _Sadigh1:
+
+**(Sadigh1)** B. Sadigh, P. Erhart, A. Stukowski, A. Caro, E. Martinez, and L. Zepeda-Ruiz, Phys. Rev. B **85**, 184203 (2012)
+
+.. _Sadigh2:
+
+**(Sadigh2)** B. Sadigh and P. Erhart, Phys. Rev. B **86**, 134204 (2012)

doc/src/kspace_style.rst

@@ -283,7 +283,7 @@ parameters and how to choose them is described in
 ----------
 
 The *electrode* styles add methods that are required for the constant potential
-method implemented in :doc:`fix electrode/* <fix_electrode_conp>`.  The styles
+method implemented in :doc:`fix electrode/* <fix_electrode>`.  The styles
 *ewald/electrode*, *pppm/electrode* and *pppm/electrode/intel* are available.
 These styles do not support the `kspace_modify slab nozforce` command.
 

doc/src/pair_coul.rst

@@ -174,11 +174,11 @@ shifted force model described in :ref:`Fennell <Fennell1>`, given by:
    E = q_iq_j \left[ \frac{\mbox{erfc} (\alpha r)}{r} -  \frac{\mbox{erfc} (\alpha r_c)}{r_c} +
    \left( \frac{\mbox{erfc} (\alpha r_c)}{r_c^2} +  \frac{2\alpha}{\sqrt{\pi}}\frac{\exp (-\alpha^2    r^2_c)}{r_c} \right)(r-r_c) \right] \qquad r < r_c
 
-where :math:`\alpha` is the damping parameter and erfc() is the
+where :math:`\alpha` is the damping parameter and *erfc()* is the
-complementary error-function. The potential corrects issues in the
+complementary error-function. The potential corrects issues in the Wolf
-Wolf model (described below) to provide consistent forces and energies
+model (described below) to provide consistent forces and energies (the
-(the Wolf potential is not differentiable at the cutoff) and smooth
+Wolf potential is not differentiable at the cutoff) and smooth decay to
-decay to zero.
+zero.
 
 ----------
 
@@ -192,30 +192,32 @@ summation method, described in :ref:`Wolf <Wolf1>`, given by:
    \frac{1}{2} \sum_{j \neq i}
    \frac{q_i q_j {\rm erf}(\alpha r_{ij})}{r_{ij}} \qquad r < r_c
 
-where :math:`\alpha` is the damping parameter, and erc() and erfc() are
+where :math:`\alpha` is the damping parameter, and *erf()* and *erfc()*
-error-function and complementary error-function terms.  This potential
+are error-function and complementary error-function terms.  This
-is essentially a short-range, spherically-truncated,
+potential is essentially a short-range, spherically-truncated,
 charge-neutralized, shifted, pairwise *1/r* summation.  With a
 manipulation of adding and subtracting a self term (for i = j) to the
-first and second term on the right-hand-side, respectively, and a
+first and second term on the right-hand-side, respectively, and a small
-small enough :math:`\alpha` damping parameter, the second term shrinks and
+enough :math:`\alpha` damping parameter, the second term shrinks and the
-the potential becomes a rapidly-converging real-space summation.  With
+potential becomes a rapidly-converging real-space summation.  With a
-a long enough cutoff and small enough :math:`\alpha` parameter, the energy and
+long enough cutoff and small enough :math:`\alpha` parameter, the energy
-forces calculated by the Wolf summation method approach those of the
+and forces calculated by the Wolf summation method approach those of the
 Ewald sum.  So it is a means of getting effective long-range
 interactions with a short-range potential.
 
 ----------
 
-Style *coul/streitz* is the Coulomb pair interaction defined as part
+Style *coul/streitz* is the Coulomb pair interaction defined as part of
-of the Streitz-Mintmire potential, as described in :ref:`this paper <Streitz2>`, in which charge distribution about an atom is modeled
+the Streitz-Mintmire potential, as described in :ref:`this paper
-as a Slater 1\ *s* orbital.  More details can be found in the referenced
+<Streitz2>`, in which charge distribution about an atom is modeled as a
+Slater 1\ *s* orbital.  More details can be found in the referenced
 paper.  To fully reproduce the published Streitz-Mintmire potential,
-which is a variable charge potential, style *coul/streitz* must be
+which is a variable charge potential, style *coul/streitz* must be used
-used with :doc:`pair_style eam/alloy <pair_eam>` (or some other
+with :doc:`pair_style eam/alloy <pair_eam>` (or some other short-range
-short-range potential that has been parameterized appropriately) via
+potential that has been parameterized appropriately) via the
-the :doc:`pair_style hybrid/overlay <pair_hybrid>` command.  Likewise,
+:doc:`pair_style hybrid/overlay <pair_hybrid>` command.  Likewise,
-charge equilibration must be performed via the :doc:`fix qeq/slater <fix_qeq>` command. For example:
+charge equilibration must be performed via the :doc:`fix qeq/slater
+<fix_qeq>` command. For example:
 
 .. code-block:: LAMMPS
 

doc/src/pair_pod.rst

@@ -0,0 +1,97 @@
+.. index:: pair_style pod
+
+pair_style pod command
+========================
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   pair_style pod
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   pair_style pod
+   pair_coeff * * Ta_param.pod Ta_coefficients.pod Ta
+
+Description
+"""""""""""
+
+.. versionadded:: TBD
+
+Pair style *pod* defines the proper orthogonal descriptor (POD)
+potential :ref:`(Nguyen) <Nguyen20221>`.  The mathematical definition of
+the POD potential is described from :doc:`fitpod <fitpod_command>`, which is
+used to fit the POD potential to *ab initio* energy and force data.
+
+Only a single pair_coeff command is used with the *pod* style which
+specifies a POD parameter file followed by a coefficient file.
+
+The coefficient file (``Ta_coefficients.pod``) contains coefficients for the
+POD potential. The top of the coefficient file can contain any number of
+blank and comment lines (start with #), but follows a strict format
+after that. The first non-blank non-comment line must contain:
+
+* POD_coefficients: *ncoeff*
+
+This is followed by *ncoeff* coefficients, one per line. The coefficient
+file is generated after training the POD potential using :doc:`fitpod
+<fitpod_command>`.
+
+The POD parameter file (``Ta_param.pod``) can contain blank and comment lines
+(start with #) anywhere. Each non-blank non-comment line must contain
+one keyword/value pair. See :doc:`fitpod <fitpod_command>` for the description
+of all the keywords that can be assigned in the parameter file.
+
+As an example, if a LAMMPS indium phosphide simulation has 4 atoms
+types, with the first two being indium and the third and fourth being
+phophorous, the pair_coeff command would look like this:
+
+.. code-block:: LAMMPS
+
+   pair_coeff * * pod InP_param.pod InP_coefficients.pod In In P P
+
+The first 2 arguments must be \* \* so as to span all LAMMPS atom types.
+The two filenames are for the parameter and coefficient files, respectively.
+The two trailing 'In' arguments map LAMMPS atom types 1 and 2 to the
+POD 'In' element. The two trailing 'P' arguments map LAMMPS atom types
+3 and 4 to the POD 'P' element.
+
+If a POD mapping value is specified as NULL, the mapping is not
+performed.  This can be used when a *pod* potential is used as part of
+the *hybrid* pair style.  The NULL values are placeholders for atom
+types that will be used with other potentials.
+
+Examples about training and using POD potentials are found in the
+directory lammps/examples/PACKAGES/pod.
+
+----------
+
+Restrictions
+""""""""""""
+
+This style is part of the ML-POD package.  It is only enabled if LAMMPS
+was built with that package. See the :doc:`Build package
+<Build_package>` page for more info.
+
+This pair style does not compute per-atom energies and per-atom stresses.
+
+Related commands
+""""""""""""""""
+
+:doc:`fitpod <fitpod_command>`,
+
+Default
+"""""""
+
+none
+
+----------
+
+.. _Nguyen20221:
+
+**(Nguyen)** Nguyen and Rohskopf, arXiv preprint arXiv:2209.02362 (2022).

doc/src/pair_style.rst

@@ -314,6 +314,7 @@ accelerated styles exist.
 * :doc:`oxrna2/xstk <pair_oxrna2>` -
 * :doc:`pace <pair_pace>` - Atomic Cluster Expansion (ACE) machine-learning potential
 * :doc:`pace/extrapolation <pair_pace>` - Atomic Cluster Expansion (ACE) machine-learning potential with extrapolation grades
+* :doc:`pod <pair_pod>` - Proper orthogonal decomposition (POD) machine-learning potential
 * :doc:`peri/eps <pair_peri>` - peridynamic EPS potential
 * :doc:`peri/lps <pair_peri>` - peridynamic LPS potential
 * :doc:`peri/pmb <pair_peri>` - peridynamic PMB potential

doc/src/pair_sw.rst

@@ -176,9 +176,13 @@ are placeholders for atom types that will be used with other potentials.
 
 .. note::
 
-   When the *threebody off* keyword is used, multiple pair_coeff commands may
+   When the *threebody off* keyword is used, multiple pair_coeff
-   be used to specific the pairs of atoms which don't require three-body term.
+   commands may be used to specific the pairs of atoms which don't
-   In these cases, the first 2 arguments are not required to be \* \*.
+   require three-body term.  In these cases, the first 2 arguments are
+   not required to be \* \*, the potential parameter file is only read
+   by the first :doc:`pair_coeff command <pair_coeff>` and the element
+   to atom type mappings must be consistent across all *pair_coeff*
+   statements.  If not LAMMPS will abort with an error.
 
 Stillinger-Weber files in the *potentials* directory of the LAMMPS
 distribution have a ".sw" suffix.  Lines that are not blank or

doc/src/pair_table.rst

@@ -120,6 +120,13 @@ best effect:
 
 ----------
 
+Suitable tables in the correct format for use with these pair styles can
+be created by LAMMPS itself using the :doc:`pair_write <pair_write>`
+command.  In combination with the :doc:`pair style python <pair_python>`
+this can be a powerful mechanism to implement and test tables for use
+with LAMMPS.  Another option to generate tables is the Python code in
+the ``tools/tabulate`` folder of the LAMMPS source code distribution.
+
 The format of a tabulated file has an (optional) header followed by a
 series of one or more sections, defined as follows (without the
 parenthesized comments). The header must start with a `#` character

doc/src/pair_tracker.rst

@@ -174,8 +174,8 @@ the specified attribute.
 Restrictions
 """"""""""""
 
-This fix is part of the MISC package.  It is only enabled if LAMMPS
+This pair style is part of the MISC package.  It is only enabled if
-was built with that package.  See the :doc:`Build package
+LAMMPS was built with that package.  See the :doc:`Build package
 <Build_package>` page for more info.
 
 This pair style is currently incompatible with granular pair styles

doc/src/pair_ylz.rst

@@ -82,7 +82,7 @@ mixing as described below:
 
 * :math:`\epsilon` = well depth (energy units)
 * :math:`\sigma` = minimum effective particle radii (distance units)
-* :math:`\zeta` = tune parameter for the slope of the attractive branch
+* :math:`\zeta` = tuning parameter for the slope of the attractive branch
 * :math:`\mu` = parameter related to bending rigidity
 * :math:`\beta` = parameter related to the spontaneous curvature
 * cutoff (distance units)

doc/src/python.rst

@@ -8,14 +8,25 @@ Syntax
 
 .. parsed-literal::
 
-   python func keyword args ...
+   python mode keyword args ...
 
-* func = name of Python function
+* mode = *source* or name of Python function
-* one or more keyword/args pairs must be appended
+
+  if mode is *source*:
 
   .. parsed-literal::
 
-     keyword = *invoke* or *input* or *return* or *format* or *length* or *file* or *here* or *exists* or *source*
+     keyword = *here* or name of a *Python file*
+       *here* arg = inline
+         inline = one or more lines of Python code which defines func
+                  must be a single argument, typically enclosed between triple quotes
+       *Python file* = name of a file with Python code which will be executed immediately
+
+* if *mode* is the name of a Python function, one or more keywords with/without arguments must be appended
+
+  .. parsed-literal::
+
+     keyword = *invoke* or *input* or *return* or *format* or *length* or *file* or *here* or *exists*
        *invoke* arg = none = invoke the previously defined Python function
        *input* args = N i1 i2 ... iN
          N = # of inputs to function
@@ -24,7 +35,7 @@ Syntax
            SELF = reference to LAMMPS itself which can be accessed by Python function
            variable = v_name, where name = name of LAMMPS variable, e.g. v_abc
        *return* arg = varReturn
-         varReturn = v_name  = LAMMPS variable name which return value of function will be assigned to
+         varReturn = v_name  = LAMMPS variable name which the return value of the Python function will be assigned to
        *format* arg = fstring with M characters
          M = N if no return value, where N = # of inputs
          M = N+1 if there is a return value
@@ -38,10 +49,6 @@ Syntax
          inline = one or more lines of Python code which defines func
                   must be a single argument, typically enclosed between triple quotes
        *exists* arg = none = Python code has been loaded by previous python command
-       *source* arg = *filename* or *inline*
-         filename = file of Python code which will be executed immediately
-         inline = one or more lines of Python code which will be executed immediately
-                  must be a single argument, typically enclosed between triple quotes
 
 Examples
 """"""""
@@ -70,80 +77,105 @@ Examples
        lmp.command("pair_style lj/cut ${cut}")   # LAMMPS commands
        lmp.command("pair_coeff * * 1.0 1.0")
        lmp.command("run 100")
-    """
+   """
+
+   python source funcdef.py
+
+   python source here "from lammps import lammps"
+
 
 Description
 """""""""""
 
-Define a Python function or execute a previously defined function or
+The *python* command allows interfacing LAMMPS with an embedded Python
-execute some arbitrary python code.
+interpreter and enables either executing arbitrary python code in that
+interpreter, registering a Python function for future execution (as a
+python style variable, from a fix interfaced with python, or for direct
+invocation), or invoking such a previously registered function.
+
 Arguments, including LAMMPS variables, can be passed to the function
-from the LAMMPS input script and a value returned by the Python
+from the LAMMPS input script and a value returned by the Python function
-function to a LAMMPS variable.  The Python code for the function can
+assigned to a LAMMPS variable.  The Python code for the function can be included
-be included directly in the input script or in a separate Python file.
+directly in the input script or in a separate Python file.  The function
-The function can be standard Python code or it can make "callbacks" to
+can be standard Python code or it can make "callbacks" to LAMMPS through
-LAMMPS through its library interface to query or set internal values
+its library interface to query or set internal values within LAMMPS.
-within LAMMPS.  This is a powerful mechanism for performing complex
+This is a powerful mechanism for performing complex operations in a
-operations in a LAMMPS input script that are not possible with the
+LAMMPS input script that are not possible with the simple input script
-simple input script and variable syntax which LAMMPS defines.  Thus
+and variable syntax which LAMMPS defines.  Thus your input script can
-your input script can operate more like a true programming language.
+operate more like a true programming language.
 
 Use of this command requires building LAMMPS with the PYTHON package
 which links to the Python library so that the Python interpreter is
 embedded in LAMMPS.  More details about this process are given below.
 
 There are two ways to invoke a Python function once it has been
-defined.  One is using the *invoke* keyword.  The other is to assign
+registered.  One is using the *invoke* keyword.  The other is to assign
 the function to a :doc:`python-style variable <variable>` defined in
-your input script.  Whenever the variable is evaluated, it will
+your input script.  Whenever the variable is evaluated, it will execute
-execute the Python function to assign a value to the variable.  Note
+the Python function to assign a value to the variable.  Note that
-that variables can be evaluated in many different ways within LAMMPS.
+variables can be evaluated in many different ways within LAMMPS.  They
-They can be substituted for directly in an input script.  Or they can
+can be substituted with their result directly in an input script, or
-be passed to various commands as arguments, so that the variable is
+they can be passed to various commands as arguments, so that the
-evaluated during a simulation run.
+variable is evaluated during a simulation run.
 
-A broader overview of how Python can be used with LAMMPS is given on
+A broader overview of how Python can be used with LAMMPS is given in the
-the :doc:`Python <Python_head>` doc page.  There is an examples/python
+:doc:`Use Python with LAMMPS <Python_head>` section of the
-directory which illustrates use of the python command.
+documentation.  There also is an ``examples/python`` directory which
+illustrates use of the python command.
 
 ----------
 
-The *func* setting specifies the name of the Python function.  The
+The first argument of the *python* command is either the *source*
-code for the function is defined using the *file* or *here* keywords
+keyword or the name of a Python function.  This defines the mode
-as explained below. In case of the *source* keyword, the name of
+of the python command.
-the function is ignored.
+
+.. versionchanged:: TBD
+
+If the *source* keyword is used, it is followed by either a file name or
+the *here* keyword.  No other keywords can be used.  The *here* keyword
+is followed by a string with python commands, either on a single line
+enclosed in quotes, or as multiple lines enclosed in triple quotes.
+These Python commands will be passed to the python interpreter and
+executed immediately without registering a Python function for future
+execution.  The code will be loaded into and run in the "main" module of
+the Python interpreter.  This allows running arbitrary Python code at
+any time while processing the LAMMPS input file.  This can be used to
+pre-load Python modules, initialize global variables, define functions
+or classes, or perform operations using the python programming language.
+The Python code will be executed in parallel on all MPI processes.  No
+arguments can be passed.
+
+In all other cases, the first argument is the name of a Python function
+that will be registered with LAMMPS for future execution.  The function
+may already be defined (see *exists* keyword) or must be defined using
+the *file* or *here* keywords as explained below.
 
 If the *invoke* keyword is used, no other keywords can be used, and a
-previous python command must have defined the Python function
+previous *python* command must have registered the Python function
 referenced by this command.  This invokes the Python function with the
-previously defined arguments and return value processed as explained
+previously defined arguments and the return value is processed as
-below.  You can invoke the function as many times as you wish in your
+explained below.  You can invoke the function as many times as you wish
-input script.
+in your input script.
-
-If the *source* keyword is used, no other keywords can be used.
-The argument can be a filename or a string with python commands,
-either on a single line enclosed in quotes, or as multiple lines
-enclosed in triple quotes. These python commands will be passed
-to the python interpreter and executed immediately without registering
-a python function for future execution.
 
 The *input* keyword defines how many arguments *N* the Python function
-expects.  If it takes no arguments, then the *input* keyword should
+expects.  If it takes no arguments, then the *input* keyword should not
-not be used.  Each argument can be specified directly as a value,
+be used.  Each argument can be specified directly as a value, e.g. '6'
-e.g. 6 or 3.14159 or abc (a string of characters).  The type of each
+or '3.14159' or 'abc' (a string of characters).  The type of each
 argument is specified by the *format* keyword as explained below, so
 that Python will know how to interpret the value.  If the word SELF is
 used for an argument it has a special meaning.  A pointer is passed to
-the Python function which it converts into a reference to LAMMPS
+the Python function which it can convert into a reference to LAMMPS
-itself.  This enables the function to call back to LAMMPS through its
+itself using the :doc:`LAMMPS Python module <Python_module>`.  This
-library interface as explained below.  This allows the Python function
+enables the function to call back to LAMMPS through its library
-to query or set values internal to LAMMPS which can affect the
+interface as explained below.  This allows the Python function to query
-subsequent execution of the input script.  A LAMMPS variable can also
+or set values internal to LAMMPS which can affect the subsequent
-be used as an argument, specified as v_name, where "name" is the name
+execution of the input script.  A LAMMPS variable can also be used as an
-of the variable.  Any style of LAMMPS variable can be used, as defined
+argument, specified as v_name, where "name" is the name of the variable.
-by the :doc:`variable <variable>` command.  Each time the Python
+Any style of LAMMPS variable returning a scalar or a string can be used,
-function is invoked, the LAMMPS variable is evaluated and its value is
+as defined by the :doc:`variable <variable>` command.  The *format*
-passed to the Python function.
+keyword must be used to set the type of data that is passed to Python.
+Each time the Python function is invoked, the LAMMPS variable is
+evaluated and its value is passed to the Python function.
 
 The *return* keyword is only needed if the Python function returns a
 value.  The specified *varReturn* must be of the form v_name, where
@@ -153,8 +185,9 @@ numeric or string value, as specified by the *format* keyword.
 
 As explained on the :doc:`variable <variable>` doc page, the definition
 of a python-style variable associates a Python function name with the
-variable.  This must match the *func* setting for this command.  For
+variable.  This must match the *Python function name* first argument of
-example these two commands would be self-consistent:
+the *python* command.  For example these two commands would be
+consistent:
 
 .. code-block:: LAMMPS
 
@@ -163,21 +196,22 @@ example these two commands would be self-consistent:
 
 The two commands can appear in either order in the input script so
 long as both are specified before the Python function is invoked for
-the first time.
+the first time.  Afterwards, the variable 'foo' is associated with
+the Python function 'myMultiply'.
 
-The *format* keyword must be used if the *input* or *return* keyword
+The *format* keyword must be used if the *input* or *return* keywords
-is used.  It defines an *fstring* with M characters, where M = sum of
+are used.  It defines an *fstring* with M characters, where M = sum of
 number of inputs and outputs.  The order of characters corresponds to
 the N inputs, followed by the return value (if it exists).  Each
 character must be one of the following: "i" for integer, "f" for
-floating point, "s" for string, or "p" for SELF.  Each character
+floating point, "s" for string, or "p" for SELF.  Each character defines
-defines the type of the corresponding input or output value of the
+the type of the corresponding input or output value of the Python
-Python function and affects the type conversion that is performed
+function and affects the type conversion that is performed internally as
-internally as data is passed back and forth between LAMMPS and Python.
+data is passed back and forth between LAMMPS and Python.  Note that it
-Note that it is permissible to use a :doc:`python-style variable <variable>` in a LAMMPS command that allows for an
+is permissible to use a :doc:`python-style variable <variable>` in a
-equal-style variable as an argument, but only if the output of the
+LAMMPS command that allows for an equal-style variable as an argument,
-Python function is flagged as a numeric value ("i" or "f") via the
+but only if the output of the Python function is flagged as a numeric
-*format* keyword.
+value ("i" or "f") via the *format* keyword.
 
 If the *return* keyword is used and the *format* keyword specifies the
 output as a string, then the default maximum length of that string is
@@ -192,12 +226,13 @@ truncated.
 
 Either the *file*, *here*, or *exists* keyword must be used, but only
 one of them.  These keywords specify what Python code to load into the
-Python interpreter.  The *file* keyword gives the name of a file,
+Python interpreter.  The *file* keyword gives the name of a file
-which should end with a ".py" suffix, which contains Python code.  The
+containing Python code, which should end with a ".py" suffix.  The code
-code will be immediately loaded into and run in the "main" module of
+will be immediately loaded into and run in the "main" module of the
-the Python interpreter.  Note that Python code which contains a
+Python interpreter.  The Python code will be executed in parallel on all
-function definition does not "execute" the function when it is run; it
+MPI processes.  Note that Python code which contains a function
-simply defines the function so that it can be invoked later.
+definition does not "execute" the function when it is run; it simply
+defines the function so that it can be invoked later.
 
 The *here* keyword does the same thing, except that the Python code
 follows as a single argument to the *here* keyword.  This can be done
@@ -208,14 +243,15 @@ proper indentation, blank lines, and comments, as desired.  See the
 how triple quotes can be used as part of input script syntax.
 
 The *exists* keyword takes no argument.  It means that Python code
-containing the required Python function defined by the *func* setting,
+containing the required Python function with the given name has already
-is assumed to have been previously loaded by another python command.
+been executed, for example by a *python source* command or in the same
+file that was used previously with the *file* keyword.
 
-Note that the Python code that is loaded and run must contain a
+Note that the Python code that is loaded and run must contain a function
-function with the specified *func* name.  To operate properly when
+with the specified function name.  To operate properly when later
-later invoked, the function code must match the *input* and
+invoked, the function code must match the *input* and *return* and
-*return* and *format* keywords specified by the python command.
+*format* keywords specified by the python command.  Otherwise Python
-Otherwise Python will generate an error.
+will generate an error.
 
 ----------
 
@@ -225,19 +261,19 @@ LAMMPS.
 Whether you load Python code from a file or directly from your input
 script, via the *file* and *here* keywords, the code can be identical.
 It must be indented properly as Python requires.  It can contain
-comments or blank lines.  If the code is in your input script, it
+comments or blank lines.  If the code is in your input script, it cannot
-cannot however contain triple-quoted Python strings, since that will
+however contain triple-quoted Python strings, since that will conflict
-conflict with the triple-quote parsing that the LAMMPS input script
+with the triple-quote parsing that the LAMMPS input script performs.
-performs.
 
 All the Python code you specify via one or more python commands is
-loaded into the Python "main" module, i.e. __main__.  The code can
+loaded into the Python "main" module, i.e. ``__name__ == '__main__'``.
-define global variables or statements that are outside of function
+The code can define global variables, define global functions, define
-definitions.  It can contain multiple functions, only one of which
+classes or execute statements that are outside of function definitions.
-matches the *func* setting in the python command.  This means you can
+It can contain multiple functions, only one of which matches the *func*
-use the *file* keyword once to load several functions, and the
+setting in the python command.  This means you can use the *file*
-*exists* keyword thereafter in subsequent python commands to access
+keyword once to load several functions, and the *exists* keyword
-the other functions previously loaded.
+thereafter in subsequent python commands to register the other functions
+that were previously loaded with LAMMPS.
 
 A Python function you define (or more generally, the code you load)
 can import other Python modules or classes, it can make calls to other
@@ -264,12 +300,13 @@ outside the function:
      nvaluelast = nvalue
      return nvalue
 
-Nsteplast stores the previous timestep the function was invoked
+The variable 'nsteplast' stores the previous timestep the function was
-(passed as an argument to the function).  Nvaluelast stores the return
+invoked (passed as an argument to the function).  The variable
-value computed on the last function invocation.  If the function is
+'nvaluelast' stores the return value computed on the last function
-invoked again on the same timestep, the previous value is simply
+invocation.  If the function is invoked again on the same timestep, the
-returned, without re-computing it.  The "global" statement inside the
+previous value is simply returned, without re-computing it.  The
-Python function allows it to overwrite the global variables.
+"global" statement inside the Python function allows it to overwrite the
+global variables from within the local context of the function.
 
 Note that if you load Python code multiple times (via multiple python
 commands), you can overwrite previously loaded variables and functions
@@ -285,19 +322,39 @@ copy of the Python function(s) you define.  There is no connection
 between the Python interpreters running on different processors.
 This implies three important things.
 
-First, if you put a print statement in your Python function, you will
+First, if you put a print or other statement creating output to the
-see P copies of the output, when running on P processors.  If the
+screen in your Python function, you will see P copies of the output,
-prints occur at (nearly) the same time, the P copies of the output may
+when running on P processors.  If the prints occur at (nearly) the same
-be mixed together.  Welcome to the world of parallel programming and
+time, the P copies of the output may be mixed together.  When loading
-debugging.
+the LAMMPS Python module into the embedded Python interpreter, it is
+possible to pass the pointer to the current LAMMPS class instance and
+via the Python interface to the LAMMPS library interface, it is possible
+to determine the MPI rank of the current process and thus adapt the
+Python code so that output will only appear on MPI rank 0.  The
+following LAMMPS input demonstrates how this could be done. The text
+'Hello, LAMMPS!' should be printed only once, even when running LAMMPS
+in parallel.
 
-Second, if your Python code loads modules that are not pre-loaded by
+.. code-block:: LAMMPS
-the Python library, then it will load the module from disk.  This may
+
-be a bottleneck if 1000s of processors try to load a module at the
+   python python_hello input 1 SELF format p here """
-same time.  On some large supercomputers, loading of modules from disk
+   def python_hello(handle):
-by Python may be disabled.  In this case you would need to pre-build a
+       from lammps import lammps
-Python library that has the required modules pre-loaded and link
+       lmp = lammps(ptr=handle)
-LAMMPS with that library.
+       me = lmp.extract_setting('world_rank')
+       if me == 0:
+           print('Hello, LAMMPS!')
+   """
+
+   python python_hello invoke
+
+If your Python code loads Python modules that are not pre-loaded by the
+Python library, then it will load the module from disk.  This may be a
+bottleneck if 1000s of processors try to load a module at the same time.
+On some large supercomputers, loading of modules from disk by Python may
+be disabled.  In this case you would need to pre-build a Python library
+that has the required modules pre-loaded and link LAMMPS with that
+library.
 
 Third, if your Python code calls back to LAMMPS (discussed in the
 next section) and causes LAMMPS to perform an MPI operation requires
@@ -315,22 +372,21 @@ Python function is as follows:
 
 .. code-block:: python
 
-   def foo(lmpptr,...):
+   def foo(handle,...):
      from lammps import lammps
-     lmp = lammps(ptr=lmpptr)
+     lmp = lammps(ptr=handle)
      lmp.command('print "Hello from inside Python"')
      ...
 
-The function definition must include a variable (lmpptr in this case)
+The function definition must include a variable ('handle' in this case)
-which corresponds to SELF in the python command.  The first line of the
+which corresponds to SELF in the *python* command.  The first line of
-function imports the :doc:`"lammps" Python module <Python_module>`.
+the function imports the :doc:`"lammps" Python module <Python_module>`.
-The second line creates a Python object ``lmp`` which
+The second line creates a Python object ``lmp`` which wraps the instance
-wraps the instance of LAMMPS that called the function.  The "ptr=lmpptr"
+of LAMMPS that called the function.  The 'ptr=handle' argument is what
-argument is what makes that happen.  The third line invokes the
+makes that happen.  The third line invokes the command() function in the
-command() function in the LAMMPS library interface.  It takes a single
+LAMMPS library interface.  It takes a single string argument which is a
-string argument which is a LAMMPS input script command for LAMMPS to
+LAMMPS input script command for LAMMPS to execute, the same as if it
-execute, the same as if it appeared in your input script.  In this case,
+appeared in your input script.  In this case, LAMMPS should output
-LAMMPS should output
 
 .. parsed-literal::
 
@@ -344,8 +400,8 @@ The :doc:`Python_head` page describes the syntax
 for how Python wraps the various functions included in the LAMMPS
 library interface.
 
-A more interesting example is in the examples/python/in.python script
+A more interesting example is in the ``examples/python/in.python`` script
-which loads and runs the following function from examples/python/funcs.py:
+which loads and runs the following function from ``examples/python/funcs.py``:
 
 .. code-block:: python
 
@@ -495,24 +551,35 @@ Restrictions
 """"""""""""
 
 This command is part of the PYTHON package.  It is only enabled if
-LAMMPS was built with that package.  See the :doc:`Build package <Build_package>` page for more info.
+LAMMPS was built with that package.  See the :doc:`Build package
+<Build_package>` page for more info.
 
-Building LAMMPS with the PYTHON package will link LAMMPS with the
+Building LAMMPS with the PYTHON package will link LAMMPS with the Python
-Python library on your system.  Settings to enable this are in the
+library on your system.  Settings to enable this are in the
 lib/python/Makefile.lammps file.  See the lib/python/README file for
 information on those settings.
 
-If you use Python code which calls back to LAMMPS, via the SELF input argument
+If you use Python code which calls back to LAMMPS, via the SELF input
-explained above, there is an extra step required when building LAMMPS.  LAMMPS
+argument explained above, there is an extra step required when building
-must also be built as a shared library and your Python function must be able to
+LAMMPS.  LAMMPS must also be built as a shared library and your Python
-load the :doc:`"lammps" Python module <Python_module>` that wraps the LAMMPS
+function must be able to load the :doc:`"lammps" Python module
-library interface.  These are the same steps required to use Python by itself
+<Python_module>` that wraps the LAMMPS library interface.  These are the
-to wrap LAMMPS.  Details on these steps are explained on the :doc:`Python
+same steps required to use Python by itself to wrap LAMMPS.  Details on
-<Python_head>` doc page.  Note that it is important that the stand-alone LAMMPS
+these steps are explained on the :doc:`Python <Python_head>` doc page.
-executable and the LAMMPS shared library be consistent (built from the same
+Note that it is important that the stand-alone LAMMPS executable and the
-source code files) in order for this to work.  If the two have been built at
+LAMMPS shared library be consistent (built from the same source code
+files) in order for this to work.  If the two have been built at
 different times using different source files, problems may occur.
 
+Another limitation of calling back to Python from the LAMMPS module
+using the *python* command in a LAMMPS input is that both, the Python
+interpreter and LAMMPS, must be linked to the same Python runtime as a
+shared library.  If the Python interpreter is linked to Python
+statically (which seems to happen with Conda) then loading the shared
+LAMMPS library will create a second python "main" module that hides the
+one from the Python interpreter and all previous defined function and
+global variables will become invisible.
+
 Related commands
 """"""""""""""""
 

doc/src/read_data.rst

@@ -340,16 +340,20 @@ and are called "tilt factors" because they are the amount of
 displacement applied to faces of an originally orthogonal box to
 transform it into the parallelepiped.
 
-By default, the tilt factors (xy,xz,yz) can not skew the box more than
+The tilt factors (xy,xz,yz) should not skew the box more than half the
-half the distance of the corresponding parallel box length.  For
+distance of the corresponding parallel box length.  For example, if
-example, if xlo = 2 and xhi = 12, then the x box length is 10 and the
+:math:`x_\text{lo} = 2` and :math:`x_\text{hi} = 12`, then the :math:`x`
-xy tilt factor must be between -5 and 5.  Similarly, both xz and yz
+box length is 10 and the :math:`xy` tilt factor must be between
-must be between -(xhi-xlo)/2 and +(yhi-ylo)/2.  Note that this is not
+:math:`-5` and :math:`5`.  Similarly, both :math:`xz` and :math:`yz`
-a limitation, since if the maximum tilt factor is 5 (as in this
+must be between :math:`-(x_\text{hi}-x_\text{lo})/2` and
-example), then configurations with tilt = ..., -15, -5, 5, 15, 25,
+:math:`+(y_\text{hi}-y_\text{lo})/2`.  Note that this is not a
-... are all geometrically equivalent.  If you wish to define a box
+limitation, since if the maximum tilt factor is 5 (as in this example),
-with tilt factors that exceed these limits, you can use the :doc:`box tilt <box>` command, with a setting of *large*\ ; a setting of
+then configurations with tilt :math:`= \dots, -15`, :math:`-5`,
-*small* is the default.
+:math:`5`, :math:`15`, :math:`25, \dots` are all geometrically
+equivalent.  Simulations with large tilt factors will run inefficiently,
+since they require more ghost atoms and thus more communication.  With
+very large tilt factors, LAMMPS will eventually produce incorrect
+trajectories and stop with errors due to lost atoms or similar.
 
 See the :doc:`Howto triclinic <Howto_triclinic>` page for a
 geometric description of triclinic boxes, as defined by LAMMPS, and

doc/src/reset_atom_ids.rst

@@ -1,94 +0,0 @@
-.. index:: reset_atom_ids
-
-reset_atom_ids command
-======================
-
-Syntax
-""""""
-
-.. code-block:: LAMMPS
-
-   reset_atom_ids keyword values ...
-
-   * zero or more keyword/value pairs may be appended
-   * keyword = *sort*
-
-.. parsed-literal::
-
-   *sort* value = *yes* or *no*
-
-Examples
-""""""""
-
-.. code-block:: LAMMPS
-
-   reset_atom_ids
-   reset_atom_ids sort yes
-
-Description
-"""""""""""
-
-Reset atom IDs for the system, including all the global IDs stored
-for bond, angle, dihedral, improper topology data.  This will
-create a set of IDs that are numbered contiguously from 1 to N
-for a N atoms system.
-
-This can be useful to do after performing a "delete_atoms" command for
-a molecular system.  The delete_atoms compress yes option will not
-perform this operation due to the existence of bond topology.  It can
-also be useful to do after any simulation which has lost atoms,
-e.g. due to atoms moving outside a simulation box with fixed
-boundaries (see the "boundary command"), or due to evaporation (see
-the "fix evaporate" command).
-
-If the *sort* keyword is used with a setting of *yes*, then the
-assignment of new atom IDs will be the same no matter how many
-processors LAMMPS is running on.  This is done by first doing a
-spatial sort of all the atoms into bins and sorting them within each
-bin.  Because the set of bins is independent of the number of
-processors, this enables a consistent assignment of new IDs to each
-atom.
-
-This can be useful to do after using the "create_atoms" command and/or
-"replicate" command.  In general those commands do not guarantee
-assignment of the same atom ID to the same physical atom when LAMMPS
-is run on different numbers of processors.  Enforcing consistent IDs
-can be useful for debugging or comparing output from two different
-runs.
-
-Note that the spatial sort requires communication of atom IDs and
-coordinates between processors in an all-to-all manner.  This is done
-efficiently in LAMMPS, but it is more expensive than how atom IDs are
-reset without sorting.
-
-Note that whether sorting or not, the resetting of IDs is not a
-compression, where gaps in atom IDs are removed by decrementing atom
-IDs that are larger.  Instead the IDs for all atoms are erased, and
-new IDs are assigned so that the atoms owned by an individual
-processor have consecutive IDs, as the :doc:`create_atoms
-<create_atoms>` command explains.
-
-.. note::
-
-   If this command is used before a :doc:`pair style <pair_style>` is
-   defined, an error about bond topology atom IDs not being found may
-   result.  This is because the cutoff distance for ghost atom
-   communication was not sufficient to find atoms in bonds, angles, etc
-   that are owned by other processors.  The :doc:`comm_modify cutoff <comm_modify>` command can be used to correct this issue.
-   Or you can define a pair style before using this command.  If you do
-   the former, you should unset the comm_modify cutoff after using
-   reset_atom_ids so that subsequent communication is not inefficient.
-
-Restrictions
-""""""""""""
-none
-
-Related commands
-""""""""""""""""
-
-:doc:`delete_atoms <delete_atoms>`
-
-Default
-"""""""
-
-By default, *sort* is no.

doc/src/reset_atoms.rst

@@ -0,0 +1,283 @@
+.. index:: reset_atoms
+
+reset_atoms command
+===================
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   reset_atoms property arguments ...
+
+* property = *id* or *image* or *mol*
+* additional arguments depend on the property
+
+  .. code-block:: LAMMPS
+
+     reset_atoms id keyword value ...
+
+  * zero or more keyword/value pairs can be appended
+  * keyword = *sort*
+
+    .. parsed-literal::
+
+       *sort* value = *yes* or *no*
+
+  .. code-block:: LAMMPS
+
+     reset_atoms image group-ID
+
+  * group-ID = ID of group of atoms whose image flags will be reset
+
+  .. code-block:: LAMMPS
+
+     reset atoms mol group-ID keyword value ...
+
+  * group-ID = ID of group of atoms whose molecule IDs will be reset
+  * zero or more keyword/value pairs can be appended
+  * keyword = *compress* or *offset* or *single*
+
+    .. parsed-literal::
+
+       *compress* value = *yes* or *no*
+       *offset* value = *Noffset* >= -1
+       *single* value = *yes* or *no* to treat single atoms (no bonds) as molecules
+
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   reset_atoms id
+   reset_atoms id sort yes
+   reset_atoms image all
+   reset_atoms image mobile
+   reset_atoms mol all
+   reset_atoms mol all offset 10 single yes
+   reset_atoms mol solvent compress yes offset 100
+   reset_atoms mol solvent compress no
+
+
+Description
+"""""""""""
+
+.. versionadded:: TBD
+
+The *reset_atoms* command resets the values of a specified atom
+property.  In contrast to the set command, it does this in a
+collective manner which resets the values for many atoms in a
+self-consistent way.  This is often useful when the simulated system
+has undergone significant modifications like adding or removing atoms
+or molecules, joining data files, changing bonds, or large-scale
+diffusion.
+
+The new values can be thought of as a *reset*, similar to values atoms
+would have if a new data file were being read or a new simulation
+performed.  Note that the set command also resets atom properties to
+new values, but it treats each atom independently.
+
+The *property* setting can be *id* or *image* or *mol*.  For *id*, the
+IDs of all the atoms are reset to contiguous values.  For *image*, the
+image flags of atoms in the specified *group-ID* are reset so that at
+least one atom in each molecule is in the simulation box (image flag =
+0).  For *mol*, the molecule IDs of all atoms are reset to contiguous
+values.
+
+More details on these operations and their arguments or optional
+keyword/value settings are given below.
+
+----------
+
+*Property id*
+
+Reset atom IDs for the entire system, including all the global IDs
+stored for bond, angle, dihedral, improper topology data.  This will
+create a set of IDs that are numbered contiguously from 1 to N for a N
+atoms system.
+
+This can be useful to do after performing a "delete_atoms" command for
+a molecular system.  The delete_atoms compress yes option will not
+perform this operation due to the existence of bond topology.  It can
+also be useful to do after any simulation which has lost atoms,
+e.g. due to atoms moving outside a simulation box with fixed
+boundaries (see the "boundary command"), or due to evaporation (see
+the "fix evaporate" command).
+
+If the *sort* keyword is used with a setting of *yes*, then the
+assignment of new atom IDs will be the same no matter how many
+processors LAMMPS is running on.  This is done by first doing a
+spatial sort of all the atoms into bins and sorting them within each
+bin.  Because the set of bins is independent of the number of
+processors, this enables a consistent assignment of new IDs to each
+atom.
+
+This can be useful to do after using the "create_atoms" command and/or
+"replicate" command.  In general those commands do not guarantee
+assignment of the same atom ID to the same physical atom when LAMMPS
+is run on different numbers of processors.  Enforcing consistent IDs
+can be useful for debugging or comparing output from two different
+runs.
+
+Note that the spatial sort requires communication of atom IDs and
+coordinates between processors in an all-to-all manner.  This is done
+efficiently in LAMMPS, but it is more expensive than how atom IDs are
+reset without sorting.
+
+Note that whether sorting or not, the resetting of IDs is not a
+compression, where gaps in atom IDs are removed by decrementing atom
+IDs that are larger.  Instead the IDs for all atoms are erased, and
+new IDs are assigned so that the atoms owned by an individual
+processor have consecutive IDs, as the :doc:`create_atoms
+<create_atoms>` command explains.
+
+.. note::
+
+   If this command is used before a :doc:`pair style <pair_style>` is
+   defined, an error about bond topology atom IDs not being found may
+   result.  This is because the cutoff distance for ghost atom
+   communication was not sufficient to find atoms in bonds, angles, etc
+   that are owned by other processors.  The :doc:`comm_modify cutoff
+   <comm_modify>` command can be used to correct this issue.  Or you can
+   define a pair style before using this command.  If you do the former,
+   you should unset the *comm_modify cutoff* after using *reset
+   atoms id* so that subsequent communication is not inefficient.
+
+----------
+
+*Property image*
+
+Reset the image flags of atoms so that at least one atom in each
+molecule has an image flag of 0.  Molecular topology is respected so
+that if the molecule straddles a periodic simulation box boundary, the
+images flags of all atoms in the molecule will be consistent.  This
+avoids inconsistent image flags that could result from resetting all
+image flags to zero with the :doc:`set <set>` command.
+
+.. note::
+
+   If the system has no bonds, there is no reason to use this command,
+   since image flags for different atoms do not need to be
+   consistent.  Use the :doc:`set <set>` command with its *image*
+   keyword instead.
+
+Only image flags for atoms in the specified *group-ID* are reset; all
+others remain unchanged.  No check is made for whether the group
+covers complete molecule fragments and thus whether the command will
+result in inconsistent image flags.
+
+Molecular fragments are identified by the algorithm used by the
+:doc:`compute fragment/atom <compute_cluster_atom>` command.  For each
+fragment the average of the largest and the smallest image flag in
+each direction across all atoms in the fragment is computed and
+subtracted from the current image flag in the same direction.
+
+This can be a useful operation to perform after running longer
+equilibration runs of mobile systems where molecules would pass
+through the system multiple times and thus produce non-zero image
+flags.
+
+.. note::
+
+   Same as explained for the :doc:`compute fragment/atom
+   <compute_cluster_atom>` command, molecules are identified using the
+   current bond topology.  This will **not** account for bonds broken by
+   the :doc:`bond_style quartic <bond_quartic>` command, because this
+   bond style does not perform a full update of the bond topology data
+   structures within LAMMPS.  In that case, using the :doc:`delete_bonds
+   all bond 0 remove <delete_bonds>` will permanently delete such
+   broken bonds and should thus be used first.
+
+----------
+
+*Property mol*
+
+Reset molecule IDs for a specified group of atoms based on current
+bond connectivity.  This will typically create a new set of molecule
+IDs for atoms in the group.  Only molecule IDs for atoms in the
+specified *group-ID* are reset; molecule IDs for atoms not in the
+group are not changed.
+
+For purposes of this operation, molecules are identified by the current
+bond connectivity in the system, which may or may not be consistent with
+the current molecule IDs.  A molecule in this context is a set of atoms
+connected to each other with explicit bonds.  The specific algorithm
+used is the one of :doc:`compute fragment/atom <compute_cluster_atom>`
+Once the molecules are identified and a new molecule ID computed for
+each, this command will update the current molecule ID for all atoms in
+the group with the new molecule ID.  Note that if the group excludes
+atoms within molecules, one (physical) molecule may become two or more
+(logical) molecules.  For example if the group excludes atoms in the
+middle of a linear chain, then each end of the chain is considered an
+independent molecule and will be assigned a different molecule ID.
+
+This can be a useful operation to perform after running reactive
+molecular dynamics run with :doc:`fix bond/react <fix_bond_react>`,
+:doc:`fix bond/create <fix_bond_create>`, or :doc:`fix bond/break
+<fix_bond_break>`, all of which can change molecule topologies. It can
+also be useful after molecules have been deleted with the
+:doc:`delete_atoms <delete_atoms>` command or after a simulation which
+has lost molecules, e.g. via the :doc:`fix evaporate <fix_evaporate>`
+command.
+
+The *compress* keyword determines how new molecule IDs are computed.  If
+the setting is *yes* (the default) and there are N molecules in the
+group, the new molecule IDs will be a set of N contiguous values.  See
+the *offset* keyword for details on selecting the range of these values.
+If the setting is *no*, the molecule ID of every atom in the molecule
+will be set to the smallest atom ID of any atom in the molecule.
+
+The *single* keyword determines whether single atoms (not bonded to
+another atom) are treated as one-atom molecules or not, based on the
+*yes* or *no* setting.  If the setting is *no* (the default), their
+molecule IDs are set to 0.  This setting can be important if the new
+molecule IDs will be used as input to other commands such as
+:doc:`compute chunk/atom molecule <compute_chunk_atom>` or :doc:`fix
+rigid molecule <fix_rigid>`.
+
+The *offset* keyword is only used if the *compress* setting is *yes*.
+Its default value is *Noffset* = -1.  In that case, if the specified
+group is *all*, then the new compressed molecule IDs will range from 1
+to N.  If the specified group is not *all* and the largest molecule ID
+of atoms outside that group is M, then the new compressed molecule IDs will
+range from M+1 to M+N, to avoid collision with existing molecule
+IDs.  If an *Noffset* >= 0 is specified, then the new compressed
+molecule IDs will range from *Noffset*\ +1 to *Noffset*\ +N.  If the group
+is not *all* there may be collisions with the molecule IDs of other atoms.
+
+.. note::
+
+   Same as explained for the :doc:`compute fragment/atom
+   <compute_cluster_atom>` command, molecules are identified using the
+   current bond topology.  This will **not** account for bonds broken by
+   the :doc:`bond_style quartic <bond_quartic>` command, because this
+   bond style does not perform a full update of the bond topology data
+   structures within LAMMPS.  In that case, using the :doc:`delete_bonds
+   all bond 0 remove <delete_bonds>` will permanently delete such broken
+   bonds and should thus be used first.
+
+
+Restrictions
+""""""""""""
+
+The *image* property can only be used when the atom style supports bonds.
+
+Related commands
+""""""""""""""""
+
+:doc:`compute fragment/atom <compute_cluster_atom>`
+:doc:`fix bond/react <fix_bond_react>`,
+:doc:`fix bond/create <fix_bond_create>`,
+:doc:`fix bond/break <fix_bond_break>`,
+:doc:`fix evaporate <fix_evaporate>`,
+:doc:`delete_atoms <delete_atoms>`,
+:doc:`delete_bonds <delete_bonds>`
+
+Defaults
+""""""""
+
+For property *id*, the default keyword setting is sort = no.
+
+For property *mol*, the default keyword settings are compress = yes,
+single = no, and offset = -1.

doc/src/reset_mol_ids.rst

@@ -1,116 +0,0 @@
-.. index:: reset_mol_ids
-
-reset_mol_ids command
-=====================
-
-Syntax
-""""""
-
-.. parsed-literal::
-
-   reset_mol_ids group-ID keyword value ...
-
-* group-ID = ID of group of atoms whose molecule IDs will be reset
-* zero or more keyword/value pairs may be appended
-* keyword = *compress* or *offset* or *single*
-
-  .. parsed-literal::
-
-       *compress* value = *yes* or *no*
-       *offset* value = *Noffset* >= -1
-       *single* value = *yes* or *no* to treat single atoms (no bonds) as molecules
-
-Examples
-""""""""
-
-.. code-block:: LAMMPS
-
-   reset_mol_ids all
-   reset_mol_ids all offset 10 single yes
-   reset_mol_ids solvent compress yes offset 100
-   reset_mol_ids solvent compress no
-
-Description
-"""""""""""
-
-Reset molecule IDs for a group of atoms based on current bond
-connectivity.  This will typically create a new set of molecule IDs
-for atoms in the group.  Only molecule IDs for atoms in the specified
-group are reset; molecule IDs for atoms not in the group are not
-changed.
-
-For purposes of this operation, molecules are identified by the current
-bond connectivity in the system, which may or may not be consistent with
-the current molecule IDs.  A molecule in this context is a set of atoms
-connected to each other with explicit bonds.  The specific algorithm
-used is the one of :doc:`compute fragment/atom <compute_cluster_atom>`
-Once the molecules are identified and a new molecule ID computed for
-each, this command will update the current molecule ID for all atoms in
-the group with the new molecule ID.  Note that if the group excludes
-atoms within molecules, one (physical) molecule may become two or more
-(logical) molecules.  For example if the group excludes atoms in the
-middle of a linear chain, then each end of the chain is considered an
-independent molecule and will be assigned a different molecule ID.
-
-This can be a useful operation to perform after running reactive
-molecular dynamics run with :doc:`fix bond/react <fix_bond_react>`,
-:doc:`fix bond/create <fix_bond_create>`, or :doc:`fix bond/break
-<fix_bond_break>`, all of which can change molecule topologies. It can
-also be useful after molecules have been deleted with the
-:doc:`delete_atoms <delete_atoms>` command or after a simulation which
-has lost molecules, e.g. via the :doc:`fix evaporate <fix_evaporate>`
-command.
-
-The *compress* keyword determines how new molecule IDs are computed.  If
-the setting is *yes* (the default) and there are N molecules in the
-group, the new molecule IDs will be a set of N contiguous values.  See
-the *offset* keyword for details on selecting the range of these values.
-If the setting is *no*, the molecule ID of every atom in the molecule
-will be set to the smallest atom ID of any atom in the molecule.
-
-The *single* keyword determines whether single atoms (not bonded to
-another atom) are treated as one-atom molecules or not, based on the
-*yes* or *no* setting.  If the setting is *no* (the default), their
-molecule IDs are set to 0.  This setting can be important if the new
-molecule IDs will be used as input to other commands such as
-:doc:`compute chunk/atom molecule <compute_chunk_atom>` or :doc:`fix
-rigid molecule <fix_rigid>`.
-
-The *offset* keyword is only used if the *compress* setting is *yes*.
-Its default value is *Noffset* = -1.  In that case, if the specified
-group is *all*, then the new compressed molecule IDs will range from 1
-to N.  If the specified group is not *all* and the largest molecule ID
-of atoms outside that group is M, then the new compressed molecule IDs will
-range from M+1 to M+N, to avoid collision with existing molecule
-IDs.  If an *Noffset* >= 0 is specified, then the new compressed
-molecule IDs will range from *Noffset*\ +1 to *Noffset*\ +N.  If the group
-is not *all* there may be collisions with the molecule IDs of other atoms.
-
-.. note::
-
-   The same as explained for the :doc:`compute fragment/atom
-   <compute_cluster_atom>` command, molecules are identified using the
-   current bond topology.  This will not account for bonds broken by
-   the :doc:`bond_style quartic <bond_quartic>` command because it
-   does not perform a full update of the bond topology data structures
-   within LAMMPS.
-
-Restrictions
-""""""""""""
-none
-
-Related commands
-""""""""""""""""
-
-:doc:`reset_atom_ids <reset_atom_ids>`, :doc:`fix bond/react <fix_bond_react>`,
-:doc:`fix bond/create <fix_bond_create>`,
-:doc:`fix bond/break <fix_bond_break>`,
-:doc:`fix evaporate <fix_evaporate>`,
-:doc:`delete_atoms <delete_atoms>`,
-:doc:`compute fragment/atom <compute_cluster_atom>`
-
-Default
-"""""""
-
-The default keyword settings are compress = yes, single = no, and
-offset = -1.

doc/utils/sphinx-config/LAMMPSLexer.py

@@ -1,25 +1,33 @@
 from pygments.lexer import RegexLexer, words, include, default
 from pygments.token import *
 
-LAMMPS_COMMANDS = ("angle_coeff", "angle_style", "atom_modify", "atom_style",
+LAMMPS_COMMANDS = ("angle_coeff", "angle_style", "atom_modify",
-"balance", "bond_coeff", "bond_style", "bond_write", "boundary", "box",
+                   "atom_style", "balance", "bond_coeff", "bond_style",
-"clear", "comm_modify", "comm_style",
+                   "bond_write", "boundary", "clear", "comm_modify",
-"compute_modify", "create_atoms", "create_bonds", "create_box", "delete_atoms",
+                   "comm_style", "compute_modify", "create_atoms",
-"delete_bonds", "dielectric", "dihedral_coeff", "dihedral_style", "dimension",
+                   "create_bonds", "create_box", "delete_atoms",
-"displace_atoms", "dump_modify", "dynamical_matrix", "echo", "elif", "else",
+                   "delete_bonds", "dielectric", "dihedral_coeff",
-"fix_modify", "group2ndx", "hyper", "if", "improper_coeff",
+                   "dihedral_style", "dimension", "displace_atoms",
-"improper_style", "include", "info", "jump", "kim",
+                   "dump_modify", "dynamical_matrix", "echo", "elif",
-"kspace_modify", "kspace_style", "label", "labelmap", "lattice", "log",
+                   "else", "fix_modify", "group2ndx", "hyper", "if",
-"mass", "mdi", "message", "minimize", "min_modify", "min_style", "molecule",
+                   "improper_coeff", "improper_style", "include",
-"ndx2group", "neb", "neb/spin", "neighbor", "neigh_modify", "newton", "next",
+                   "info", "jump", "kim", "kspace_modify",
-"package", "pair_coeff", "pair_modify", "pair_style", "pair_write",
+                   "kspace_style", "label", "labelmap", "lattice",
-"partition", "prd", "print", "processors", "python", "quit", "read_data",
+                   "log", "mass", "mdi", "message", "minimize",
-"read_dump", "read_restart", "replicate", "rerun", "reset_ids",
+                   "min_modify", "min_style", "molecule", "ndx2group",
-"reset_timestep", "restart", "run", "run_style", "server", "set", "shell",
+                   "neb", "neb/spin", "neighbor", "neigh_modify",
-"special_bonds", "suffix", "tad", "temper", "temper/grem", "temper/npt", "then",
+                   "newton", "next", "package", "pair_coeff",
-"thermo", "thermo_modify", "thermo_style", "third_order", "timer", "timestep",
+                   "pair_modify", "pair_style", "pair_write",
-"units", "velocity", "write_coeff",
+                   "partition", "plugin", "prd", "print", "processors",
-"write_data", "write_restart")
+                   "python", "quit", "read_data", "read_dump",
+                   "read_restart", "replicate", "rerun", "reset_atoms",
+                   "reset_timestep", "restart", "run", "run_style",
+                   "server", "set", "shell", "special_bonds", "suffix",
+                   "tad", "temper", "temper/grem", "temper/npt", "then",
+                   "thermo", "thermo_modify", "thermo_style",
+                   "third_order", "timer", "timestep", "units",
+                   "velocity", "write_coeff", "write_data",
+                   "write_restart")
 
 #fix ID group-ID style args
 #compute ID group-ID style args

doc/utils/sphinx-config/false_positives.txt

@@ -228,6 +228,7 @@ Bartels
 Bartelt
 barycenter
 barye
+basename
 Bashford
 bashrc
 Baskes
@@ -256,6 +257,7 @@ berlin
 Berne
 Bertotti
 Bessarab
+bessel
 Beutler
 Bext
 Bfrac
@@ -447,6 +449,7 @@ checkbox
 checkmark
 checkqeq
 checksum
+chemflag
 chemistries
 Chemnitz
 Cheng
@@ -608,6 +611,7 @@ curv
 Cusentino
 customIDs
 cutbond
+cutghost
 cuthi
 cutinner
 cutlo
@@ -688,6 +692,7 @@ delocalized
 Delong
 delr
 deltaHf
+deltamu
 dem
 Dendrimer
 dendritic
@@ -738,6 +743,7 @@ diel
 Dietz
 differentiable
 diffusively
+diffusivities
 diffusivity
 dihedral
 dihedrals
@@ -749,6 +755,7 @@ dimensionality
 dimensioned
 dimgray
 dipolar
+dipoleflag
 dir
 Direc
 dirname
@@ -757,6 +764,7 @@ discretization
 discretized
 discretizing
 disp
+dispersionflag
 dissipative
 Dissipative
 distharm
@@ -819,6 +827,7 @@ du
 dU
 Ducastelle
 Dudarev
+Dufils
 Duin
 Dullweber
 dumpfile
@@ -906,6 +915,7 @@ elastance
 Electroneg
 electronegative
 electronegativity
+electroneutrality
 Eleftheriou
 ElementN
 elementset
@@ -971,7 +981,6 @@ equilization
 equipartitioning
 eradius
 erate
-erc
 Ercolessi
 Erdmann
 erf
@@ -1022,12 +1031,14 @@ evirials
 ew
 ewald
 Ewald
+ewaldflag
 excitations
 excv
 exe
 executables
 extep
 extrema
+extxyz
 exy
 ey
 ez
@@ -1094,8 +1105,10 @@ Fincham
 Fint
 fingerprintconstants
 fingerprintsperelement
+finitecutflag
 Finnis
 Fiorin
+fitpod
 fixID
 fj
 Fji
@@ -1138,6 +1151,7 @@ Forschungszentrum
 fortran
 Fortran
 Fosado
+fourbody
 fourier
 fp
 fphi
@@ -1271,6 +1285,7 @@ greenyellow
 Greffet
 grem
 gREM
+Grepl
 Grest
 Grigera
 Grimme
@@ -1510,6 +1525,9 @@ inumeric
 inv
 invariants
 inversed
+invertible
+invertibility
+ionicities
 ionizable
 ionocovalent
 iostreams
@@ -1629,6 +1647,7 @@ Kalia
 Kamberaj
 Kantorovich
 Kapfer
+Karhunen
 Karls
 Karlsruhe
 Karniadakis
@@ -1883,6 +1902,7 @@ ln
 localhost
 localTemp
 localvectors
+Loeve
 Loewen
 logfile
 logfreq
@@ -1934,6 +1954,7 @@ Mackrodt
 MacOS
 Macromolecules
 macroparticle
+Maday
 Madura
 Magda
 Magdeburg
@@ -2198,6 +2219,7 @@ msd
 msi
 MSI
 msm
+msmflag
 msse
 msst
 Mtchell
@@ -2248,6 +2270,7 @@ MxN
 myCompute
 myIndex
 mylammps
+myMultiply
 MyPool
 mysocket
 mySpin
@@ -2276,6 +2299,8 @@ nanometer
 nanometers
 nanoparticle
 nanoparticles
+nanopores
+nanostructures
 nanotube
 Nanotube
 nanotubes
@@ -2359,6 +2384,7 @@ Ng
 nghost
 Nghost
 Ngpu
+ngpus
 Ngyuen
 nh
 nharmonic
@@ -2384,6 +2410,7 @@ nktv
 nl
 nlayers
 nlen
+Nlimit
 nlines
 Nlines
 nlo
@@ -2473,7 +2500,9 @@ nsq
 Nstart
 nstats
 Nstep
-Nsteplast
+Nsteps
+nsteps
+nsteplast
 Nstop
 nsub
 Nsw
@@ -2503,7 +2532,7 @@ numpy
 Numpy
 Nurdin
 Nvalue
-Nvaluelast
+nvaluelast
 Nvalues
 nvc
 nvcc
@@ -2547,6 +2576,7 @@ Omelyan
 omp
 OMP
 oneAPI
+onebody
 onelevel
 oneway
 onlysalt
@@ -2628,6 +2658,7 @@ Pastewka
 pathangle
 pathname
 pathnames
+Patera
 Patomtrans
 Pattnaik
 Pavese
@@ -2769,6 +2800,7 @@ PowerShell
 ppme
 ppn
 pppm
+pppmflag
 Prakash
 Praprotnik
 prd
@@ -2923,6 +2955,7 @@ Rcmx
 Rcmy
 Rco
 Rcut
+rcut
 rcutfac
 rdc
 rdf
@@ -2948,6 +2981,7 @@ refactoring
 reflectionstyle
 Reinders
 reinit
+reinitflag
 relaxbox
 relink
 relres
@@ -3015,6 +3049,7 @@ Rij
 RIj
 Rik
 Rin
+rin
 Rinaldi
 Rino
 RiRj
@@ -3090,6 +3125,7 @@ Rutuparna
 rx
 rxd
 rxnave
+rxnbond
 rxnsum
 ry
 Ryckaert
@@ -3150,6 +3186,7 @@ sdpd
 SDPD
 se
 seagreen
+Searles
 Secor
 sectoring
 sed
@@ -3178,6 +3215,7 @@ setvel
 sfftw
 sfree
 Sg
+sgcmc
 Shan
 Shanno
 Shapeev
@@ -3280,6 +3318,7 @@ SPH
 spica
 SPICA
 Spickermann
+spinflag
 splined
 spparks
 Sprik
@@ -3465,6 +3504,7 @@ thermo
 thermochemical
 thermochemistry
 thermodynamically
+thermopotentiostat
 Thermophysical
 thermostatted
 thermostatting
@@ -3832,6 +3872,7 @@ workflows
 Workum
 Worley
 Wriggers
+writedata
 Wuppertal
 Wurtzite
 www

examples/PACKAGES/electrode/README

@@ -1,14 +1,14 @@
 These examples demonstrate the use of the ELECTRODE package for constant potential molecular dynamics.
 
 planar/
-  au-vac.data   -- gold electrodes with vacuum
+  data.au-vac -- gold electrodes with vacuum
   in.planar* -- comparison of gold electrodes with vacuum to theoretical capacitance of planar capacitor
     -- 5x, further labeled by long-range solver (ewald / pppm) and boundary correction (ew2d / ew3dc / ffield)
     -- the pppm-ew2d combination would not give correct results and will throw an error if selected
   test.sh -- run all in.planar files and check charge at 1.2V and %difference from theoretical (last column)
 
 graph-il/
-  graph-il.data -- graphene electrodes with electrolyte (coarse-grained BMIm-PF6)
+  data.graph-il -- graphene electrodes with electrolyte (coarse-grained BMIm-PF6)
   in.conp    -- reference run at constant potential
   in.etypes  -- type-based smart neighborlists
   in.ffield  -- finite field method with fully periodic cell
@@ -18,10 +18,22 @@ graph-il/
   in.thermo  -- thermalize electrolyte with thermopotentiostat instead of NVT
 
 au-aq/
-  au-aq.data -- gold electrodes with electrolyte (SPC water + NaCl)
+  data.au-aq -- gold electrodes with electrolyte (SPC water + NaCl)
   in.ffield  -- finite field method with fully periodic cell
   in.tf      -- Thomas-Fermi metallicity model with more delocalized charges
 
+madelung/
+  data.au-elyt -- tiny electrodes with two electrolyte atoms in between
+  settings.mod -- common settings
+  in.*         -- setup KSpace and fix electrode/conp
+  plate_cap.py -- compute reference energy and charges from Madelung style sum
+  eval.py      -- compare output of reference and Lammps job (used by test.sh)
+  test.sh      -- run all in.* files and check charge at 1 V and %difference from theoretical (last column)
+
+piston/
+  data.piston -- two electrodes with water
+  in.piston   -- equilibrate distance between rigid electrodes
+
 # future work:
   # in.cylinder -- comparison of carbon nanotube to theoretical induced charge for charge near circular conductor
 

examples/PACKAGES/electrode/au-aq/in.ffield

@@ -4,7 +4,7 @@
 boundary p p p # ffield uses periodic z-boundary and no slab
 include settings.mod # styles, groups, computes and fixes
 
-fix conp bot electrode/conp -1.0 1.805132 couple top 1.0 symm on ffield yes etypes 6*7
+fix conp bot electrode/conp -1.0 1.805132 couple top 1.0 symm on ffield yes etypes on
 
 thermo 50
 thermo_style custom step temp c_ctemp epair etotal c_qtop c_qbot c_qztop c_qzbot

examples/PACKAGES/electrode/au-aq/in.tf

@@ -6,7 +6,7 @@
 boundary p p p # ffield uses periodic z-boundary and no slab
 include settings.mod # styles, groups, computes and fixes
 
-fix conp bot electrode/conp -1.0 1.805132 couple top 1.0 symm on ffield yes etypes 6*7
+fix conp bot electrode/conp -1.0 1.805132 couple top 1.0 symm on ffield yes etypes on
 fix_modify conp tf 6 1.0 18.1715745
 fix_modify conp tf 7 1.0 18.1715745
 

examples/PACKAGES/electrode/au-aq/log.1Dec2022.au-aq-ffield.g++.1

@@ -1,4 +1,4 @@
-LAMMPS (24 Mar 2022)
+LAMMPS (3 Nov 2022)
 # electrodes with constant potential using finite field
 # for z-periodic gold-saline electrochemical cell
 
@@ -37,8 +37,8 @@ Finding 1-2 1-3 1-4 neighbors ...
      1 = max # of 1-3 neighbors
      1 = max # of 1-4 neighbors
      2 = max # of special neighbors
-  special bonds CPU = 0.002 seconds
+  special bonds CPU = 0.006 seconds
-  read_data CPU = 0.051 seconds
+  read_data CPU = 0.097 seconds
 
 group bot type 6
 1620 atoms in group bot
@@ -52,11 +52,12 @@ group electrolyte type 1 2 3 4 5
 
 fix nvt electrolyte nvt temp 298.0 298.0 241
 fix shake SPC shake 1e-4 20 0 b 1 2 a 1
+Finding SHAKE clusters ...
        0 = # of size 2 clusters
        0 = # of size 3 clusters
        0 = # of size 4 clusters
     2160 = # of frozen angles
-  find clusters CPU = 0.002 seconds
+  find clusters CPU = 0.006 seconds
 
 variable q atom q
 variable qz atom q*(z-lz/2)
@@ -67,12 +68,41 @@ compute qzbot bot reduce sum v_qz
 compute ctemp electrolyte temp
 
 
-fix conp bot electrode/conp -1.0 1.805132 couple top 1.0 symm on ffield yes etypes 6*7
+fix conp bot electrode/conp -1.0 1.805132 couple top 1.0 symm on ffield yes etypes on
 3240 atoms in group conp_group
 
 thermo 50
 thermo_style custom step temp c_ctemp epair etotal c_qtop c_qbot c_qztop c_qzbot
 run 500
+
+CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
+
+Your simulation uses code contributions which should be cited:
+
+- kspace_style pppm/electrode command:
+
+@article{Ahrens2021,
+author = {Ahrens-Iwers, Ludwig J.V. and Mei{\ss}ner, Robert H.},
+doi = {10.1063/5.0063381},
+title = {{Constant potential simulations on a mesh}},
+journal = {Journal of Chemical Physics},
+year = {2021}
+volume = {155},
+pages = {104104},
+}
+- fix electrode command:
+
+@article{Ahrens2022
+author = {Ahrens-Iwers, Ludwig J.V. and Janssen, Mahijs and Tee, Shern R. and Mei{\ss}ner, Robert H.},
+doi = {10.1063/5.0099239},
+title = {{ELECTRODE: An electrochemistry package for LAMMPS}},
+journal = {The Journal of Chemical Physics},
+year = {2022}
+volume = {157},
+pages = {084801},
+}
+CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
+
 PPPM/electrode initialization ...
   using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
   G vector (1/distance) = 0.24017705
@@ -82,9 +112,9 @@ PPPM/electrode initialization ...
   estimated relative force accuracy = 1.093542e-07
   using double precision MKL FFT
   3d grid and FFT values/proc = 472567 349920
-  generated 21 of 21 mixed pair_coeff terms from geometric mixing rule
+Generated 21 of 21 mixed pair_coeff terms from geometric mixing rule
 Neighbor list info ...
-  update every 1 steps, delay 10 steps, check yes
+  update: every = 1 steps, delay = 0 steps, check = yes
   max neighbors/atom: 2000, page size: 100000
   master list distance cutoff = 17
   ghost atom cutoff = 17
@@ -105,35 +135,35 @@ Neighbor list info ...
       pair build: skip
       stencil: none
       bin: none
-Per MPI rank memory allocation (min/avg/max) = 194.6 | 194.6 | 194.6 Mbytes
+Per MPI rank memory allocation (min/avg/max) = 194.8 | 194.8 | 194.8 Mbytes
    Step          Temp         c_ctemp         E_pair         TotEng         c_qtop         c_qbot        c_qztop        c_qzbot    
-         0   171.61215      298.06731     -39212.819     -35306.164      4.1391573     -4.1391573      78.718381      131.56372    
+         0   171.61215      298.06731     -39021.917     -35115.261      4.1391573     -4.1391573      78.718381      131.56372    
-        50   147.03139      255.37383     -39870.139     -36523.051      4.1312167     -4.1312167      78.563872      131.30255    
+        50   147.03139      255.37383     -39679.603     -36332.515      4.1312167     -4.1312167      78.563872      131.30255    
-       100   149.89027      260.33932     -39878.859     -36466.689      4.0217834     -4.0217834      76.482548      127.82573    
+       100   149.89027      260.33932     -39693.369     -36281.2        4.0217834     -4.0217834      76.482548      127.82573    
-       150   151.7382       263.54893     -39873.178     -36418.942      4.0469977     -4.0469977      76.967548      128.59855    
+       150   151.7382       263.54893     -39686.526     -36232.29       4.0469977     -4.0469977      76.967548      128.59855    
-       200   151.7508       263.57081     -39827.015     -36372.492      4.1830375     -4.1830375      79.554159      132.93925    
+       200   151.7508       263.57081     -39634.089     -36179.566      4.1830375     -4.1830375      79.554159      132.93925    
-       250   152.61146      265.06566     -39791.293     -36317.177      4.1835865     -4.1835865      79.56665       132.97185    
+       250   152.61146      265.06566     -39598.341     -36124.226      4.1835865     -4.1835865      79.56665       132.97185    
-       300   153.51486      266.63475     -39751.841     -36257.16       4.1571861     -4.1571861      79.061431      132.12905    
+       300   153.51486      266.63475     -39560.107     -36065.426      4.1571861     -4.1571861      79.06143       132.12905    
-       350   156.35115      271.561       -39754.955     -36195.708      4.3498059     -4.3498059      82.720202      138.28678    
+       350   156.35115      271.561       -39554.338     -35995.09       4.3498059     -4.3498059      82.720202      138.28678    
-       400   156.26118      271.40474     -39690.781     -36133.582      4.3444079     -4.3444079      82.619396      138.11873    
+       400   156.26118      271.40474     -39490.412     -35933.213      4.344408      -4.344408       82.619398      138.11874    
-       450   158.54164      275.36558     -39681.083     -36071.97       4.2020488     -4.2020488      79.912674      133.55185    
+       450   158.54163      275.36557     -39487.28      -35878.167      4.2020489     -4.2020489      79.912677      133.55186    
-       500   161.40138      280.33258     -39684.185     -36009.972      4.3021924     -4.3021924      81.807527      136.7464     
+       500   161.40137      280.33257     -39485.763     -35811.55       4.3021927     -4.3021927      81.807532      136.74641    
-Loop time of 246.197 on 1 procs for 500 steps with 9798 atoms
+Loop time of 146.959 on 1 procs for 500 steps with 9798 atoms
 
-Performance: 0.175 ns/day, 136.776 hours/ns, 2.031 timesteps/s
+Performance: 0.294 ns/day, 81.644 hours/ns, 3.402 timesteps/s, 33.336 katom-step/s
-356.3% CPU use with 1 MPI tasks x no OpenMP threads
+100.0% CPU use with 1 MPI tasks x no OpenMP threads
 
 MPI task timing breakdown:
 Section |  min time  |  avg time  |  max time  |%varavg| %total
 ---------------------------------------------------------------
-Pair    | 105.64     | 105.64     | 105.64     |   0.0 | 42.91
+Pair    | 69.832     | 69.832     | 69.832     |   0.0 | 47.52
-Bond    | 0.0010592  | 0.0010592  | 0.0010592  |   0.0 |  0.00
+Bond    | 0.00091634 | 0.00091634 | 0.00091634 |   0.0 |  0.00
-Kspace  | 37.643     | 37.643     | 37.643     |   0.0 | 15.29
+Kspace  | 33.817     | 33.817     | 33.817     |   0.0 | 23.01
-Neigh   | 5.8827     | 5.8827     | 5.8827     |   0.0 |  2.39
+Neigh   | 4.2067     | 4.2067     | 4.2067     |   0.0 |  2.86
-Comm    | 0.18181    | 0.18181    | 0.18181    |   0.0 |  0.07
+Comm    | 0.12212    | 0.12212    | 0.12212    |   0.0 |  0.08
-Output  | 0.0055762  | 0.0055762  | 0.0055762  |   0.0 |  0.00
+Output  | 0.0031896  | 0.0031896  | 0.0031896  |   0.0 |  0.00
-Modify  | 96.78      | 96.78      | 96.78      |   0.0 | 39.31
+Modify  | 38.92      | 38.92      | 38.92      |   0.0 | 26.48
-Other   |            | 0.06346    |            |       |  0.03
+Other   |            | 0.05687    |            |       |  0.04
 
 Nlocal:           9798 ave        9798 max        9798 min
 Histogram: 1 0 0 0 0 0 0 0 0 0
@@ -147,4 +177,4 @@ Ave neighs/atom = 842.63544
 Ave special neighs/atom = 1.3227189
 Neighbor list builds = 22
 Dangerous builds = 0
-Total wall time: 0:19:39
+Total wall time: 0:05:33

examples/PACKAGES/electrode/au-aq/log.1Dec2022.au-aq-ffield.g++.4

@@ -1,4 +1,4 @@
-LAMMPS (24 Mar 2022)
+LAMMPS (3 Nov 2022)
 # electrodes with constant potential using finite field
 # for z-periodic gold-saline electrochemical cell
 
@@ -39,7 +39,7 @@ Finding 1-2 1-3 1-4 neighbors ...
      1 = max # of 1-4 neighbors
      2 = max # of special neighbors
   special bonds CPU = 0.002 seconds
-  read_data CPU = 0.149 seconds
+  read_data CPU = 0.118 seconds
 
 group bot type 6
 1620 atoms in group bot
@@ -53,11 +53,12 @@ group electrolyte type 1 2 3 4 5
 
 fix nvt electrolyte nvt temp 298.0 298.0 241
 fix shake SPC shake 1e-4 20 0 b 1 2 a 1
+Finding SHAKE clusters ...
        0 = # of size 2 clusters
        0 = # of size 3 clusters
        0 = # of size 4 clusters
     2160 = # of frozen angles
-  find clusters CPU = 0.003 seconds
+  find clusters CPU = 0.002 seconds
 
 variable q atom q
 variable qz atom q*(z-lz/2)
@@ -68,12 +69,41 @@ compute qzbot bot reduce sum v_qz
 compute ctemp electrolyte temp
 
 
-fix conp bot electrode/conp -1.0 1.805132 couple top 1.0 symm on ffield yes etypes 6*7
+fix conp bot electrode/conp -1.0 1.805132 couple top 1.0 symm on ffield yes etypes on
 3240 atoms in group conp_group
 
 thermo 50
 thermo_style custom step temp c_ctemp epair etotal c_qtop c_qbot c_qztop c_qzbot
 run 500
+
+CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
+
+Your simulation uses code contributions which should be cited:
+
+- kspace_style pppm/electrode command:
+
+@article{Ahrens2021,
+author = {Ahrens-Iwers, Ludwig J.V. and Mei{\ss}ner, Robert H.},
+doi = {10.1063/5.0063381},
+title = {{Constant potential simulations on a mesh}},
+journal = {Journal of Chemical Physics},
+year = {2021}
+volume = {155},
+pages = {104104},
+}
+- fix electrode command:
+
+@article{Ahrens2022
+author = {Ahrens-Iwers, Ludwig J.V. and Janssen, Mahijs and Tee, Shern R. and Mei{\ss}ner, Robert H.},
+doi = {10.1063/5.0099239},
+title = {{ELECTRODE: An electrochemistry package for LAMMPS}},
+journal = {The Journal of Chemical Physics},
+year = {2022}
+volume = {157},
+pages = {084801},
+}
+CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
+
 PPPM/electrode initialization ...
   using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
   G vector (1/distance) = 0.24017705
@@ -83,9 +113,9 @@ PPPM/electrode initialization ...
   estimated relative force accuracy = 1.093542e-07
   using double precision MKL FFT
   3d grid and FFT values/proc = 138958 87480
-  generated 21 of 21 mixed pair_coeff terms from geometric mixing rule
+Generated 21 of 21 mixed pair_coeff terms from geometric mixing rule
 Neighbor list info ...
-  update every 1 steps, delay 10 steps, check yes
+  update: every = 1 steps, delay = 0 steps, check = yes
   max neighbors/atom: 2000, page size: 100000
   master list distance cutoff = 17
   ghost atom cutoff = 17
@@ -106,35 +136,35 @@ Neighbor list info ...
       pair build: skip
       stencil: none
       bin: none
-Per MPI rank memory allocation (min/avg/max) = 118.1 | 120.6 | 123.1 Mbytes
+Per MPI rank memory allocation (min/avg/max) = 118.2 | 120.7 | 123.2 Mbytes
    Step          Temp         c_ctemp         E_pair         TotEng         c_qtop         c_qbot        c_qztop        c_qzbot    
-         0   171.61215      298.06731     -39212.819     -35306.164      4.1391573     -4.1391573      78.718381      131.56372    
+         0   171.61215      298.06731     -39021.917     -35115.261      4.1391573     -4.1391573      78.718381      131.56372    
-        50   147.03139      255.37383     -39870.139     -36523.051      4.1312167     -4.1312167      78.563872      131.30255    
+        50   147.03139      255.37383     -39679.603     -36332.515      4.1312167     -4.1312167      78.563872      131.30255    
-       100   149.89027      260.33932     -39878.859     -36466.689      4.0217834     -4.0217834      76.482548      127.82573    
+       100   149.89027      260.33932     -39693.369     -36281.2        4.0217834     -4.0217834      76.482548      127.82573    
-       150   151.7382       263.54893     -39873.178     -36418.942      4.0469977     -4.0469977      76.967548      128.59855    
+       150   151.7382       263.54893     -39686.526     -36232.29       4.0469977     -4.0469977      76.967548      128.59855    
-       200   151.7508       263.57081     -39827.015     -36372.492      4.1830375     -4.1830375      79.554159      132.93925    
+       200   151.7508       263.57081     -39634.089     -36179.566      4.1830375     -4.1830375      79.554159      132.93925    
-       250   152.61146      265.06566     -39791.293     -36317.177      4.1835865     -4.1835865      79.56665       132.97185    
+       250   152.61146      265.06566     -39598.341     -36124.226      4.1835864     -4.1835864      79.56665       132.97185    
-       300   153.51486      266.63475     -39751.841     -36257.16       4.1571861     -4.1571861      79.061431      132.12905    
+       300   153.51486      266.63475     -39560.107     -36065.426      4.1571861     -4.1571861      79.06143       132.12905    
-       350   156.35115      271.561       -39754.955     -36195.708      4.3498059     -4.3498059      82.7202        138.28678    
+       350   156.35115      271.561       -39554.338     -35995.09       4.3498059     -4.3498059      82.720201      138.28678    
-       400   156.26118      271.40474     -39690.781     -36133.582      4.3444079     -4.3444079      82.619398      138.11873    
+       400   156.26118      271.40474     -39490.412     -35933.213      4.3444079     -4.3444079      82.619397      138.11873    
-       450   158.54163      275.36558     -39681.083     -36071.97       4.2020488     -4.2020488      79.912675      133.55185    
+       450   158.54163      275.36558     -39487.279     -35878.167      4.202049      -4.202049       79.912678      133.55186    
-       500   161.40138      280.33257     -39684.185     -36009.972      4.3021924     -4.3021924      81.807527      136.7464     
+       500   161.40137      280.33256     -39485.763     -35811.55       4.3021925     -4.3021925      81.807529      136.7464     
-Loop time of 111.902 on 4 procs for 500 steps with 9798 atoms
+Loop time of 97.2399 on 4 procs for 500 steps with 9798 atoms
 
-Performance: 0.386 ns/day, 62.168 hours/ns, 4.468 timesteps/s
+Performance: 0.444 ns/day, 54.022 hours/ns, 5.142 timesteps/s, 50.381 katom-step/s
-97.2% CPU use with 4 MPI tasks x no OpenMP threads
+87.0% CPU use with 4 MPI tasks x no OpenMP threads
 
 MPI task timing breakdown:
 Section |  min time  |  avg time  |  max time  |%varavg| %total
 ---------------------------------------------------------------
-Pair    | 21.816     | 31.136     | 40.866     | 166.5 | 27.82
+Pair    | 19.363     | 28.08      | 37.126     | 160.3 | 28.88
-Bond    | 0.00073413 | 0.00080346 | 0.00084203 |   0.0 |  0.00
+Bond    | 0.00094043 | 0.00096516 | 0.00098016 |   0.0 |  0.00
-Kspace  | 29.546     | 39.137     | 48.326     | 146.4 | 34.97
+Kspace  | 31.655     | 40.554     | 49.143     | 131.3 | 41.71
-Neigh   | 2.5867     | 2.5872     | 2.5877     |   0.0 |  2.31
+Neigh   | 2.2289     | 2.2294     | 2.2297     |   0.0 |  2.29
-Comm    | 0.33289    | 0.33603    | 0.33791    |   0.3 |  0.30
+Comm    | 0.5341     | 0.54079    | 0.5478     |   0.9 |  0.56
-Output  | 0.0022537  | 0.0030028  | 0.005192   |   2.3 |  0.00
+Output  | 0.0024141  | 0.0026943  | 0.0034176  |   0.8 |  0.00
-Modify  | 38.498     | 38.635     | 38.77      |   2.2 | 34.53
+Modify  | 25.6       | 25.755     | 25.908     |   2.8 | 26.49
-Other   |            | 0.06679    |            |       |  0.06
+Other   |            | 0.07733    |            |       |  0.08
 
 Nlocal:         2449.5 ave        2908 max        2012 min
 Histogram: 2 0 0 0 0 0 0 0 0 2
@@ -148,4 +178,4 @@ Ave neighs/atom = 842.63544
 Ave special neighs/atom = 1.3227189
 Neighbor list builds = 22
 Dangerous builds = 0
-Total wall time: 0:07:48
+Total wall time: 0:03:03

examples/PACKAGES/electrode/au-aq/log.1Dec2022.au-aq-tf.g++.1

@@ -1,4 +1,4 @@
-LAMMPS (24 Mar 2022)
+LAMMPS (3 Nov 2022)
 # electrodes with constant potential using finite field
 # for z-periodic gold-saline electrochemical cell
 # using Thomas-Fermi metallicity model: electrode q and qz will be
@@ -39,8 +39,8 @@ Finding 1-2 1-3 1-4 neighbors ...
      1 = max # of 1-3 neighbors
      1 = max # of 1-4 neighbors
      2 = max # of special neighbors
-  special bonds CPU = 0.010 seconds
+  special bonds CPU = 0.003 seconds
-  read_data CPU = 0.115 seconds
+  read_data CPU = 0.065 seconds
 
 group bot type 6
 1620 atoms in group bot
@@ -54,11 +54,12 @@ group electrolyte type 1 2 3 4 5
 
 fix nvt electrolyte nvt temp 298.0 298.0 241
 fix shake SPC shake 1e-4 20 0 b 1 2 a 1
+Finding SHAKE clusters ...
        0 = # of size 2 clusters
        0 = # of size 3 clusters
        0 = # of size 4 clusters
     2160 = # of frozen angles
-  find clusters CPU = 0.010 seconds
+  find clusters CPU = 0.002 seconds
 
 variable q atom q
 variable qz atom q*(z-lz/2)
@@ -69,7 +70,7 @@ compute qzbot bot reduce sum v_qz
 compute ctemp electrolyte temp
 
 
-fix conp bot electrode/conp -1.0 1.805132 couple top 1.0 symm on ffield yes etypes 6*7
+fix conp bot electrode/conp -1.0 1.805132 couple top 1.0 symm on ffield yes etypes on
 3240 atoms in group conp_group
 fix_modify conp tf 6 1.0 18.1715745
 fix_modify conp tf 7 1.0 18.1715745
@@ -77,6 +78,35 @@ fix_modify conp tf 7 1.0 18.1715745
 thermo 50
 thermo_style custom step temp c_ctemp epair etotal c_qtop c_qbot c_qztop c_qzbot
 run 500
+
+CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
+
+Your simulation uses code contributions which should be cited:
+
+- kspace_style pppm/electrode command:
+
+@article{Ahrens2021,
+author = {Ahrens-Iwers, Ludwig J.V. and Mei{\ss}ner, Robert H.},
+doi = {10.1063/5.0063381},
+title = {{Constant potential simulations on a mesh}},
+journal = {Journal of Chemical Physics},
+year = {2021}
+volume = {155},
+pages = {104104},
+}
+- fix electrode command:
+
+@article{Ahrens2022
+author = {Ahrens-Iwers, Ludwig J.V. and Janssen, Mahijs and Tee, Shern R. and Mei{\ss}ner, Robert H.},
+doi = {10.1063/5.0099239},
+title = {{ELECTRODE: An electrochemistry package for LAMMPS}},
+journal = {The Journal of Chemical Physics},
+year = {2022}
+volume = {157},
+pages = {084801},
+}
+CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
+
 PPPM/electrode initialization ...
   using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
   G vector (1/distance) = 0.24017705
@@ -86,9 +116,9 @@ PPPM/electrode initialization ...
   estimated relative force accuracy = 1.093542e-07
   using double precision MKL FFT
   3d grid and FFT values/proc = 472567 349920
-  generated 21 of 21 mixed pair_coeff terms from geometric mixing rule
+Generated 21 of 21 mixed pair_coeff terms from geometric mixing rule
 Neighbor list info ...
-  update every 1 steps, delay 10 steps, check yes
+  update: every = 1 steps, delay = 0 steps, check = yes
   max neighbors/atom: 2000, page size: 100000
   master list distance cutoff = 17
   ghost atom cutoff = 17
@@ -109,35 +139,35 @@ Neighbor list info ...
       pair build: skip
       stencil: none
       bin: none
-Per MPI rank memory allocation (min/avg/max) = 194.6 | 194.6 | 194.6 Mbytes
+Per MPI rank memory allocation (min/avg/max) = 194.8 | 194.8 | 194.8 Mbytes
    Step          Temp         c_ctemp         E_pair         TotEng         c_qtop         c_qbot        c_qztop        c_qzbot    
-         0   171.61215      298.06731     -39190.106     -35283.45       4.0804484     -4.0804484      79.075127      131.20697    
+         0   171.61215      298.06731     -39001.911     -35095.255      4.0804484     -4.0804484      79.075127      131.20697    
-        50   147.14308      255.56782     -39849.964     -36500.334      3.9990346     -3.9990346      77.497181      128.57759    
+        50   147.14308      255.56782     -39665.525     -36315.894      3.9990346     -3.9990346      77.497181      128.57759    
-       100   149.94935      260.44194     -39857.533     -36444.019      3.8613914     -3.8613914      74.82985       124.15315    
+       100   149.94935      260.44194     -39679.441     -36265.927      3.8613914     -3.8613914      74.82985       124.15315    
-       150   151.95924      263.93285     -39855.567     -36396.299      3.8677064     -3.8677064      74.957279      124.33201    
+       150   151.95924      263.93285     -39677.184     -36217.916      3.8677064     -3.8677064      74.957279      124.33201    
-       200   151.66737      263.42591     -39802.585     -36349.961      3.99842       -3.99842        77.491015      128.54496    
+       200   151.66737      263.42591     -39618.173     -36165.549      3.99842       -3.99842        77.491015      128.54496    
-       250   152.36874      264.64408     -39763.306     -36294.716      3.9925863     -3.9925863      77.379445      128.37226    
+       250   152.36874      264.64408     -39579.164     -36110.574      3.9925863     -3.9925863      77.379445      128.37226    
-       300   153.83916      267.19802     -39737.075     -36235.012      3.94995       -3.94995        76.553896      127.00395    
+       300   153.83916      267.19802     -39554.899     -36052.836      3.94995       -3.94995        76.553896      127.00395    
-       350   155.88897      270.75827     -39722.265     -36173.539      4.0598524     -4.0598524      78.679643      130.5394     
+       350   155.88897      270.75827     -39535.02      -35986.294      4.0598524     -4.0598524      78.679643      130.5394     
-       400   156.51993      271.85415     -39674.759     -36111.669      4.1312899     -4.1312899      80.060369      132.83599    
+       400   156.51993      271.85415     -39484.219     -35921.13       4.1312898     -4.1312898      80.060368      132.83598    
-       450   160.21339      278.26919     -39697.962     -36050.793      3.9068098     -3.9068098      75.713484      125.59216    
+       450   160.21339      278.26919     -39517.776     -35870.607      3.9068098     -3.9068098      75.713484      125.59216    
-       500   161.63639      280.74075     -39669.412     -35989.849      3.9261656     -3.9261656      76.0806        126.22255    
+       500   161.63639      280.74075     -39488.333     -35808.771      3.9261656     -3.9261656      76.080599      126.22255    
-Loop time of 280.183 on 1 procs for 500 steps with 9798 atoms
+Loop time of 185.804 on 1 procs for 500 steps with 9798 atoms
 
-Performance: 0.154 ns/day, 155.657 hours/ns, 1.785 timesteps/s
+Performance: 0.233 ns/day, 103.225 hours/ns, 2.691 timesteps/s, 26.366 katom-step/s
-341.9% CPU use with 1 MPI tasks x no OpenMP threads
+100.0% CPU use with 1 MPI tasks x no OpenMP threads
 
 MPI task timing breakdown:
 Section |  min time  |  avg time  |  max time  |%varavg| %total
 ---------------------------------------------------------------
-Pair    | 119.69     | 119.69     | 119.69     |   0.0 | 42.72
+Pair    | 91         | 91         | 91         |   0.0 | 48.98
-Bond    | 0.0010952  | 0.0010952  | 0.0010952  |   0.0 |  0.00
+Bond    | 0.0010315  | 0.0010315  | 0.0010315  |   0.0 |  0.00
-Kspace  | 42.137     | 42.137     | 42.137     |   0.0 | 15.04
+Kspace  | 40.32      | 40.32      | 40.32      |   0.0 | 21.70
-Neigh   | 6.5403     | 6.5403     | 6.5403     |   0.0 |  2.33
+Neigh   | 5.6505     | 5.6505     | 5.6505     |   0.0 |  3.04
-Comm    | 0.19411    | 0.19411    | 0.19411    |   0.0 |  0.07
+Comm    | 0.15304    | 0.15304    | 0.15304    |   0.0 |  0.08
-Output  | 0.0053644  | 0.0053644  | 0.0053644  |   0.0 |  0.00
+Output  | 0.0041829  | 0.0041829  | 0.0041829  |   0.0 |  0.00
-Modify  | 111.54     | 111.54     | 111.54     |   0.0 | 39.81
+Modify  | 48.607     | 48.607     | 48.607     |   0.0 | 26.16
-Other   |            | 0.07244    |            |       |  0.03
+Other   |            | 0.06807    |            |       |  0.04
 
 Nlocal:           9798 ave        9798 max        9798 min
 Histogram: 1 0 0 0 0 0 0 0 0 0
@@ -151,4 +181,4 @@ Ave neighs/atom = 842.7079
 Ave special neighs/atom = 1.3227189
 Neighbor list builds = 23
 Dangerous builds = 0
-Total wall time: 0:21:11
+Total wall time: 0:06:18

examples/PACKAGES/electrode/au-aq/log.1Dec2022.au-aq-tf.g++.4

@@ -1,4 +1,4 @@
-LAMMPS (24 Mar 2022)
+LAMMPS (3 Nov 2022)
 # electrodes with constant potential using finite field
 # for z-periodic gold-saline electrochemical cell
 # using Thomas-Fermi metallicity model: electrode q and qz will be
@@ -41,7 +41,7 @@ Finding 1-2 1-3 1-4 neighbors ...
      1 = max # of 1-4 neighbors
      2 = max # of special neighbors
   special bonds CPU = 0.002 seconds
-  read_data CPU = 0.091 seconds
+  read_data CPU = 0.114 seconds
 
 group bot type 6
 1620 atoms in group bot
@@ -55,11 +55,12 @@ group electrolyte type 1 2 3 4 5
 
 fix nvt electrolyte nvt temp 298.0 298.0 241
 fix shake SPC shake 1e-4 20 0 b 1 2 a 1
+Finding SHAKE clusters ...
        0 = # of size 2 clusters
        0 = # of size 3 clusters
        0 = # of size 4 clusters
     2160 = # of frozen angles
-  find clusters CPU = 0.001 seconds
+  find clusters CPU = 0.002 seconds
 
 variable q atom q
 variable qz atom q*(z-lz/2)
@@ -70,7 +71,7 @@ compute qzbot bot reduce sum v_qz
 compute ctemp electrolyte temp
 
 
-fix conp bot electrode/conp -1.0 1.805132 couple top 1.0 symm on ffield yes etypes 6*7
+fix conp bot electrode/conp -1.0 1.805132 couple top 1.0 symm on ffield yes etypes on
 3240 atoms in group conp_group
 fix_modify conp tf 6 1.0 18.1715745
 fix_modify conp tf 7 1.0 18.1715745
@@ -78,6 +79,35 @@ fix_modify conp tf 7 1.0 18.1715745
 thermo 50
 thermo_style custom step temp c_ctemp epair etotal c_qtop c_qbot c_qztop c_qzbot
 run 500
+
+CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
+
+Your simulation uses code contributions which should be cited:
+
+- kspace_style pppm/electrode command:
+
+@article{Ahrens2021,
+author = {Ahrens-Iwers, Ludwig J.V. and Mei{\ss}ner, Robert H.},
+doi = {10.1063/5.0063381},
+title = {{Constant potential simulations on a mesh}},
+journal = {Journal of Chemical Physics},
+year = {2021}
+volume = {155},
+pages = {104104},
+}
+- fix electrode command:
+
+@article{Ahrens2022
+author = {Ahrens-Iwers, Ludwig J.V. and Janssen, Mahijs and Tee, Shern R. and Mei{\ss}ner, Robert H.},
+doi = {10.1063/5.0099239},
+title = {{ELECTRODE: An electrochemistry package for LAMMPS}},
+journal = {The Journal of Chemical Physics},
+year = {2022}
+volume = {157},
+pages = {084801},
+}
+CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
+
 PPPM/electrode initialization ...
   using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
   G vector (1/distance) = 0.24017705
@@ -87,9 +117,9 @@ PPPM/electrode initialization ...
   estimated relative force accuracy = 1.093542e-07
   using double precision MKL FFT
   3d grid and FFT values/proc = 138958 87480
-  generated 21 of 21 mixed pair_coeff terms from geometric mixing rule
+Generated 21 of 21 mixed pair_coeff terms from geometric mixing rule
 Neighbor list info ...
-  update every 1 steps, delay 10 steps, check yes
+  update: every = 1 steps, delay = 0 steps, check = yes
   max neighbors/atom: 2000, page size: 100000
   master list distance cutoff = 17
   ghost atom cutoff = 17
@@ -110,35 +140,35 @@ Neighbor list info ...
       pair build: skip
       stencil: none
       bin: none
-Per MPI rank memory allocation (min/avg/max) = 118.1 | 120.6 | 123.1 Mbytes
+Per MPI rank memory allocation (min/avg/max) = 118.2 | 120.7 | 123.2 Mbytes
    Step          Temp         c_ctemp         E_pair         TotEng         c_qtop         c_qbot        c_qztop        c_qzbot    
-         0   171.61215      298.06731     -39190.106     -35283.45       4.0804484     -4.0804484      79.075127      131.20697    
+         0   171.61215      298.06731     -39001.911     -35095.255      4.0804484     -4.0804484      79.075127      131.20697    
-        50   147.14308      255.56782     -39849.964     -36500.334      3.9990346     -3.9990346      77.497181      128.57759    
+        50   147.14308      255.56782     -39665.525     -36315.894      3.9990346     -3.9990346      77.497181      128.57759    
-       100   149.94935      260.44194     -39857.533     -36444.019      3.8613914     -3.8613914      74.82985       124.15315    
+       100   149.94935      260.44194     -39679.441     -36265.927      3.8613914     -3.8613914      74.82985       124.15315    
-       150   151.95924      263.93285     -39855.567     -36396.299      3.8677064     -3.8677064      74.957279      124.33201    
+       150   151.95924      263.93285     -39677.184     -36217.916      3.8677064     -3.8677064      74.957279      124.33201    
-       200   151.66737      263.42591     -39802.585     -36349.961      3.99842       -3.99842        77.491015      128.54496    
+       200   151.66737      263.42591     -39618.173     -36165.549      3.99842       -3.99842        77.491015      128.54496    
-       250   152.36874      264.64408     -39763.306     -36294.716      3.9925863     -3.9925863      77.379445      128.37226    
+       250   152.36874      264.64408     -39579.163     -36110.574      3.9925863     -3.9925863      77.379445      128.37226    
-       300   153.83916      267.19802     -39737.075     -36235.012      3.94995       -3.94995        76.553896      127.00395    
+       300   153.83916      267.19802     -39554.899     -36052.835      3.94995       -3.94995        76.553896      127.00395    
-       350   155.88897      270.75827     -39722.265     -36173.539      4.0598524     -4.0598524      78.679643      130.5394     
+       350   155.88897      270.75826     -39535.02      -35986.294      4.0598523     -4.0598523      78.679642      130.5394     
-       400   156.51993      271.85415     -39674.759     -36111.669      4.1312899     -4.1312899      80.060369      132.83599    
+       400   156.51993      271.85415     -39484.219     -35921.13       4.1312897     -4.1312897      80.060366      132.83598    
-       450   160.21339      278.26919     -39697.962     -36050.793      3.9068098     -3.9068098      75.713485      125.59216    
+       450   160.21339      278.26919     -39517.776     -35870.607      3.9068099     -3.9068099      75.713486      125.59216    
-       500   161.63639      280.74075     -39669.412     -35989.849      3.9261654     -3.9261654      76.080597      126.22255    
+       500   161.63639      280.74075     -39488.333     -35808.771      3.9261657     -3.9261657      76.080602      126.22256    
-Loop time of 110.716 on 4 procs for 500 steps with 9798 atoms
+Loop time of 104.099 on 4 procs for 500 steps with 9798 atoms
 
-Performance: 0.390 ns/day, 61.509 hours/ns, 4.516 timesteps/s
+Performance: 0.415 ns/day, 57.833 hours/ns, 4.803 timesteps/s, 47.061 katom-step/s
-97.2% CPU use with 4 MPI tasks x no OpenMP threads
+87.7% CPU use with 4 MPI tasks x no OpenMP threads
 
 MPI task timing breakdown:
 Section |  min time  |  avg time  |  max time  |%varavg| %total
 ---------------------------------------------------------------
-Pair    | 21.17      | 30.449     | 39.65      | 164.9 | 27.50
+Pair    | 20.951     | 30.326     | 40.07      | 166.7 | 29.13
-Bond    | 0.0007313  | 0.00077537 | 0.00081477 |   0.0 |  0.00
+Bond    | 0.00098259 | 0.0010706  | 0.0011926  |   0.3 |  0.00
-Kspace  | 29.854     | 38.911     | 48.058     | 143.8 | 35.14
+Kspace  | 33.465     | 43.037     | 52.268     | 137.5 | 41.34
-Neigh   | 2.7206     | 2.7213     | 2.722      |   0.0 |  2.46
+Neigh   | 2.6007     | 2.6014     | 2.6021     |   0.0 |  2.50
-Comm    | 0.33023    | 0.33225    | 0.33384    |   0.2 |  0.30
+Comm    | 0.57766    | 0.58318    | 0.58875    |   0.7 |  0.56
-Output  | 0.0024528  | 0.0027565  | 0.0035754  |   0.9 |  0.00
+Output  | 0.0022277  | 0.0024765  | 0.0031841  |   0.8 |  0.00
-Modify  | 38.091     | 38.233     | 38.365     |   2.1 | 34.53
+Modify  | 27.292     | 27.47      | 27.647     |   3.1 | 26.39
-Other   |            | 0.06636    |            |       |  0.06
+Other   |            | 0.0787     |            |       |  0.08
 
 Nlocal:         2449.5 ave        2908 max        2017 min
 Histogram: 2 0 0 0 0 0 0 0 0 2
@@ -147,9 +177,9 @@ Histogram: 2 0 0 0 0 0 0 0 0 2
 Neighs:    2.06421e+06 ave  2.7551e+06 max 1.40237e+06 min
 Histogram: 2 0 0 0 0 0 0 0 0 2
 
-Total # of neighbors = 8256853
+Total # of neighbors = 8256852
-Ave neighs/atom = 842.708
+Ave neighs/atom = 842.7079
 Ave special neighs/atom = 1.3227189
 Neighbor list builds = 23
 Dangerous builds = 0
-Total wall time: 0:08:22
+Total wall time: 0:03:12

examples/PACKAGES/electrode/graph-il/.gitignore

@@ -1 +1,2 @@
 log.lammps*
+in.temp

examples/PACKAGES/electrode/graph-il/algo_test.sh

@@ -0,0 +1,28 @@
+#!/bin/bash
+
+LMP_BIN="$1"
+NP="${2:-1}"
+echo "MPI over $NP procs:"
+for feat in conp etypes tf
+do
+  echo "Using base input file in.$feat:"
+  echo "mat_inv, log excerpts:"
+  logfile="log.algo_test.$NP.$feat"
+  mpirun -np $NP $LMP_BIN -i in.$feat -l $logfile > /dev/null 2>&1
+  grep -A2 'Per MPI rank' $logfile
+  grep -B1 'Loop time' $logfile
+  rm $logfile
+  for cgtype in mat_cg cg
+  do
+    for tol in 1e-4 1e-5 1e-6
+    do
+      echo "$cgtype, tol = $tol, log excerpts:"
+      logfile="log.algo_test.$NP.$feat.$cgtype.$tol"
+      sed '/electrode/ s/$/ algo '"$cgtype"' '"$tol"'/' in.$feat > in.temp
+      mpirun -np $NP $LMP_BIN -i in.temp -l $logfile > /dev/null 2>&1
+      grep -A2 'Per MPI rank' $logfile
+      grep -B1 'Loop time' $logfile
+      rm $logfile
+    done
+  done
+done

examples/PACKAGES/electrode/graph-il/in.conp

@@ -3,9 +3,9 @@
 
 boundary p p f # slab calculation
 include settings.mod # styles, groups, computes and fixes
-kspace_modify slab 3.0
+kspace_modify slab 3.0 # amat twostep
 
-fix conp bot electrode/conp -1.0 1.979 couple top 1.0 symm on
+fix conp bot electrode/conp -1.0 1.979 couple top 1.0 symm on #algo mat_inv
 
 thermo 50
 thermo_style custom step temp c_ctemp epair etotal c_qbot c_qtop

examples/PACKAGES/electrode/graph-il/in.conq

@@ -5,14 +5,10 @@ boundary p p f # slab calculation
 include settings.mod # styles, groups, computes and fixes
 kspace_modify slab 3.0
 
-fix conq bot electrode/conq -1.0 1.979 couple top 1.0 etypes 5 # conq doesn't take symm option
+fix conq bot electrode/conq -1.0 1.979 couple top 1.0 etypes on # symm on
-
+variable dv equal f_conq[2]-f_conq[1] 
-# ask fix conq to output electrode potentials to internal variables
+# symm on and off give different electrode potentials, but identical potential difference
-variable vbot internal 0.0
-variable vtop internal 0.0
-fix_modify conq set v bot vbot
-fix_modify conq set v top vtop
 
 thermo 50
-thermo_style custom step temp c_ctemp epair etotal c_qbot c_qtop v_vbot v_vtop
+thermo_style custom step temp c_ctemp epair etotal c_qbot c_qtop f_conq[1] f_conq[2] v_dv
 run 500