diff --git a/.github/workflows/check-vla.yml b/.github/workflows/check-vla.yml
index 26f23cc33f0..ab89018a3de 100644
--- a/.github/workflows/check-vla.yml
+++ b/.github/workflows/check-vla.yml
@@ -27,9 +27,9 @@ jobs:
- name: Install extra packages
run: |
+ sudo apt-get update
sudo apt-get install -y ccache \
libeigen3-dev \
- libgsl-dev \
libcurl4-openssl-dev \
mold \
mpi-default-bin \
diff --git a/.github/workflows/compile-msvc.yml b/.github/workflows/compile-msvc.yml
index 03af27788b2..7560bc05493 100644
--- a/.github/workflows/compile-msvc.yml
+++ b/.github/workflows/compile-msvc.yml
@@ -11,6 +11,10 @@ on:
workflow_dispatch:
+concurrency:
+ group: ${{ github.event_name }}-${{ github.workflow }}-${{ github.ref }}
+ cancel-in-progress: ${{github.event_name == 'pull_request'}}
+
jobs:
build:
name: Windows Compilation Test
diff --git a/.github/workflows/full-regression.yml b/.github/workflows/full-regression.yml
index 73e1803bb64..a6b5353b9b0 100644
--- a/.github/workflows/full-regression.yml
+++ b/.github/workflows/full-regression.yml
@@ -30,8 +30,9 @@ jobs:
- name: Install extra packages
run: |
+ sudo apt-get update
sudo apt-get install -y ccache ninja-build libeigen3-dev \
- libgsl-dev libcurl4-openssl-dev python3-dev \
+ libcurl4-openssl-dev python3-dev \
mpi-default-bin mpi-default-dev
- name: Create Build Environment
diff --git a/.github/workflows/quick-regression.yml b/.github/workflows/quick-regression.yml
index 985177b2c11..88794bfa0a9 100644
--- a/.github/workflows/quick-regression.yml
+++ b/.github/workflows/quick-regression.yml
@@ -8,6 +8,10 @@ on:
workflow_dispatch:
+concurrency:
+ group: ${{ github.event_name }}-${{ github.workflow }}-${{ github.ref }}
+ cancel-in-progress: ${{github.event_name == 'pull_request'}}
+
jobs:
build:
name: Build LAMMPS
@@ -30,8 +34,9 @@ jobs:
- name: Install extra packages
run: |
+ sudo apt-get update
sudo apt-get install -y ccache ninja-build libeigen3-dev \
- libgsl-dev libcurl4-openssl-dev python3-dev \
+ libcurl4-openssl-dev python3-dev \
mpi-default-bin mpi-default-dev
- name: Create Build Environment
diff --git a/.github/workflows/style-check.yml b/.github/workflows/style-check.yml
new file mode 100644
index 00000000000..7be2c4fc46c
--- /dev/null
+++ b/.github/workflows/style-check.yml
@@ -0,0 +1,37 @@
+# GitHub action to run checks from tools/coding_standard
+name: "Check for Programming Style Conformance"
+
+on:
+ push:
+ branches:
+ - develop
+ pull_request:
+ branches:
+ - develop
+
+ workflow_dispatch:
+
+concurrency:
+ group: ${{ github.event_name }}-${{ github.workflow }}-${{ github.ref }}
+ cancel-in-progress: ${{github.event_name == 'pull_request'}}
+
+jobs:
+ build:
+ name: Programming Style Conformance
+ if: ${{ github.repository == 'lammps/lammps' }}
+ runs-on: ubuntu-latest
+
+ steps:
+ - name: Checkout repository
+ uses: actions/checkout@v4
+ with:
+ fetch-depth: 1
+
+ - name: Run Tests
+ working-directory: src
+ shell: bash
+ run: |
+ make check-whitespace
+ make check-permissions
+ make check-homepage
+ make check-errordocs
diff --git a/.github/workflows/unittest-linux.yml b/.github/workflows/unittest-linux.yml
index 366db25a99d..ce98fcea35e 100644
--- a/.github/workflows/unittest-linux.yml
+++ b/.github/workflows/unittest-linux.yml
@@ -11,6 +11,10 @@ on:
workflow_dispatch:
+concurrency:
+ group: ${{ github.event_name }}-${{ github.workflow }}-${{ github.ref }}
+ cancel-in-progress: ${{github.event_name == 'pull_request'}}
+
jobs:
build:
name: Linux Unit Test
@@ -27,9 +31,9 @@ jobs:
- name: Install extra packages
run: |
+ sudo apt-get update
sudo apt-get install -y ccache \
libeigen3-dev \
- libgsl-dev \
libcurl4-openssl-dev \
mold \
ninja-build \
diff --git a/.github/workflows/unittest-macos.yml b/.github/workflows/unittest-macos.yml
index b0bc4b27275..0d478a9d6bb 100644
--- a/.github/workflows/unittest-macos.yml
+++ b/.github/workflows/unittest-macos.yml
@@ -11,6 +11,10 @@ on:
workflow_dispatch:
+concurrency:
+ group: ${{ github.event_name }}-${{ github.workflow }}-${{ github.ref }}
+ cancel-in-progress: ${{github.event_name == 'pull_request'}}
+
jobs:
build:
name: MacOS Unit Test
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
index c68a9253241..1bd387b5b91 100644
--- a/cmake/CMakeLists.txt
+++ b/cmake/CMakeLists.txt
@@ -497,7 +497,7 @@ if((CMAKE_CXX_COMPILER_ID STREQUAL "Intel") AND (CMAKE_CXX_STANDARD GREATER_EQUA
PROPERTIES COMPILE_OPTIONS "-std=c++14")
endif()
-if(PKG_ATC OR PKG_AWPMD OR PKG_ML-QUIP OR PKG_ML-POD OR PKG_ELECTRODE OR BUILD_TOOLS)
+if(PKG_ATC OR PKG_AWPMD OR PKG_ML-QUIP OR PKG_ML-POD OR PKG_ELECTRODE OR PKG_RHEO OR BUILD_TOOLS)
enable_language(C)
if (NOT USE_INTERNAL_LINALG)
find_package(LAPACK)
@@ -572,7 +572,7 @@ else()
endif()
foreach(PKG_WITH_INCL KSPACE PYTHON ML-IAP VORONOI COLVARS ML-HDNNP MDI MOLFILE NETCDF
- PLUMED QMMM ML-QUIP SCAFACOS MACHDYN VTK KIM COMPRESS ML-PACE LEPTON RHEO EXTRA-COMMAND)
+ PLUMED QMMM ML-QUIP SCAFACOS MACHDYN VTK KIM COMPRESS ML-PACE LEPTON EXTRA-COMMAND)
if(PKG_${PKG_WITH_INCL})
include(Packages/${PKG_WITH_INCL})
endif()
diff --git a/cmake/Modules/Documentation.cmake b/cmake/Modules/Documentation.cmake
index 400109067ff..7b8f4a5ba06 100644
--- a/cmake/Modules/Documentation.cmake
+++ b/cmake/Modules/Documentation.cmake
@@ -4,6 +4,8 @@
option(BUILD_DOC "Build LAMMPS HTML documentation" OFF)
if(BUILD_DOC)
+ option(BUILD_DOC_VENV "Build LAMMPS documentation virtual environment" ON)
+ mark_as_advanced(BUILD_DOC_VENV)
# Current Sphinx versions require at least Python 3.8
# use default (or custom) Python executable, if version is sufficient
if(Python_VERSION VERSION_GREATER_EQUAL 3.8)
@@ -18,14 +20,6 @@ if(BUILD_DOC)
find_package(Doxygen 1.8.10 REQUIRED)
file(GLOB DOC_SOURCES CONFIGURE_DEPENDS ${LAMMPS_DOC_DIR}/src/[^.]*.rst)
- add_custom_command(
- OUTPUT docenv
- COMMAND ${VIRTUALENV} docenv
- )
-
- set(DOCENV_BINARY_DIR ${CMAKE_BINARY_DIR}/docenv/bin)
- set(DOCENV_REQUIREMENTS_FILE ${LAMMPS_DOC_DIR}/utils/requirements.txt)
-
set(SPHINX_CONFIG_DIR ${LAMMPS_DOC_DIR}/utils/sphinx-config)
set(SPHINX_CONFIG_FILE_TEMPLATE ${SPHINX_CONFIG_DIR}/conf.py.in)
set(SPHINX_STATIC_DIR ${SPHINX_CONFIG_DIR}/_static)
@@ -44,14 +38,32 @@ if(BUILD_DOC)
# configure paths in conf.py, since relative paths change when file is copied
configure_file(${SPHINX_CONFIG_FILE_TEMPLATE} ${DOC_BUILD_CONFIG_FILE})
- add_custom_command(
- OUTPUT ${DOC_BUILD_DIR}/requirements.txt
- DEPENDS docenv ${DOCENV_REQUIREMENTS_FILE}
- COMMAND ${CMAKE_COMMAND} -E copy ${DOCENV_REQUIREMENTS_FILE} ${DOC_BUILD_DIR}/requirements.txt
- COMMAND ${DOCENV_BINARY_DIR}/pip $ENV{PIP_OPTIONS} install --upgrade pip
- COMMAND ${DOCENV_BINARY_DIR}/pip $ENV{PIP_OPTIONS} install --upgrade ${LAMMPS_DOC_DIR}/utils/converters
- COMMAND ${DOCENV_BINARY_DIR}/pip $ENV{PIP_OPTIONS} install -r ${DOC_BUILD_DIR}/requirements.txt --upgrade
- )
+ if(BUILD_DOC_VENV)
+ add_custom_command(
+ OUTPUT docenv
+ COMMAND ${VIRTUALENV} docenv
+ )
+
+ set(DOCENV_BINARY_DIR ${CMAKE_BINARY_DIR}/docenv/bin)
+ set(DOCENV_REQUIREMENTS_FILE ${LAMMPS_DOC_DIR}/utils/requirements.txt)
+
+ add_custom_command(
+ OUTPUT ${DOC_BUILD_DIR}/requirements.txt
+ DEPENDS docenv ${DOCENV_REQUIREMENTS_FILE}
+ COMMAND ${CMAKE_COMMAND} -E copy ${DOCENV_REQUIREMENTS_FILE} ${DOC_BUILD_DIR}/requirements.txt
+ COMMAND ${DOCENV_BINARY_DIR}/pip $ENV{PIP_OPTIONS} install --upgrade pip
+ COMMAND ${DOCENV_BINARY_DIR}/pip $ENV{PIP_OPTIONS} install --upgrade ${LAMMPS_DOC_DIR}/utils/converters
+ COMMAND ${DOCENV_BINARY_DIR}/pip $ENV{PIP_OPTIONS} install -r ${DOC_BUILD_DIR}/requirements.txt --upgrade
+ )
+
+ set(DOCENV_DEPS docenv ${DOC_BUILD_DIR}/requirements.txt)
+ if(NOT TARGET Sphinx::sphinx-build)
+ add_executable(Sphinx::sphinx-build IMPORTED GLOBAL)
+ set_target_properties(Sphinx::sphinx-build PROPERTIES IMPORTED_LOCATION "${DOCENV_BINARY_DIR}/sphinx-build")
+ endif()
+ else()
+ find_package(Sphinx)
+ endif()
set(MATHJAX_URL "https://github.com/mathjax/MathJax/archive/3.1.3.tar.gz" CACHE STRING "URL for MathJax tarball")
set(MATHJAX_MD5 "b81661c6e6ba06278e6ae37b30b0c492" CACHE STRING "MD5 checksum of MathJax tarball")
@@ -97,8 +109,8 @@ if(BUILD_DOC)
endif()
add_custom_command(
OUTPUT html
- DEPENDS ${DOC_SOURCES} docenv ${DOC_BUILD_DIR}/requirements.txt ${DOXYGEN_XML_DIR}/index.xml ${BUILD_DOC_CONFIG_FILE}
- COMMAND ${DOCENV_BINARY_DIR}/sphinx-build ${SPHINX_EXTRA_OPTS} -b html -c ${DOC_BUILD_DIR} -d ${DOC_BUILD_DIR}/doctrees ${LAMMPS_DOC_DIR}/src ${DOC_BUILD_DIR}/html
+ DEPENDS ${DOC_SOURCES} ${DOCENV_DEPS} ${DOXYGEN_XML_DIR}/index.xml ${BUILD_DOC_CONFIG_FILE}
+ COMMAND Sphinx::sphinx-build ${SPHINX_EXTRA_OPTS} -b html -c ${DOC_BUILD_DIR} -d ${DOC_BUILD_DIR}/doctrees ${LAMMPS_DOC_DIR}/src ${DOC_BUILD_DIR}/html
COMMAND ${CMAKE_COMMAND} -E create_symlink Manual.html ${DOC_BUILD_DIR}/html/index.html
COMMAND ${CMAKE_COMMAND} -E copy_directory ${LAMMPS_DOC_DIR}/src/PDF ${DOC_BUILD_DIR}/html/PDF
COMMAND ${CMAKE_COMMAND} -E remove -f ${DOXYGEN_XML_DIR}/run.stamp
diff --git a/cmake/Modules/FindSphinx.cmake b/cmake/Modules/FindSphinx.cmake
new file mode 100644
index 00000000000..3718ecc5434
--- /dev/null
+++ b/cmake/Modules/FindSphinx.cmake
@@ -0,0 +1,29 @@
+# Find sphinx-build
+find_program(Sphinx_EXECUTABLE NAMES sphinx-build
+ PATH_SUFFIXES bin
+ DOC "Sphinx documenation build executable")
+mark_as_advanced(Sphinx_EXECUTABLE)
+
+if(Sphinx_EXECUTABLE)
+ execute_process(COMMAND ${Sphinx_EXECUTABLE} --version
+ OUTPUT_VARIABLE sphinx_version
+ OUTPUT_STRIP_TRAILING_WHITESPACE
+ RESULT_VARIABLE _sphinx_version_result)
+
+ if(_sphinx_version_result)
+ message(WARNING "Unable to determine sphinx-build verison: ${_sphinx_version_result}")
+ else()
+ string(REGEX REPLACE "sphinx-build ([0-9.]+).*"
+ "\\1"
+ Sphinx_VERSION
+ "${sphinx_version}")
+ endif()
+
+ if(NOT TARGET Sphinx::sphinx-build)
+ add_executable(Sphinx::sphinx-build IMPORTED GLOBAL)
+ set_target_properties(Sphinx::sphinx-build PROPERTIES IMPORTED_LOCATION "${Sphinx_EXECUTABLE}")
+ endif()
+endif()
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(Sphinx REQUIRED_VARS Sphinx_EXECUTABLE VERSION_VAR Sphinx_VERSION)
diff --git a/cmake/Modules/Packages/RHEO.cmake b/cmake/Modules/Packages/RHEO.cmake
deleted file mode 100644
index 7639acd8bca..00000000000
--- a/cmake/Modules/Packages/RHEO.cmake
+++ /dev/null
@@ -1,2 +0,0 @@
-find_package(GSL 2.6 REQUIRED)
-target_link_libraries(lammps PRIVATE GSL::gsl)
diff --git a/cmake/presets/mingw-cross.cmake b/cmake/presets/mingw-cross.cmake
index 100ce136320..413744b078f 100644
--- a/cmake/presets/mingw-cross.cmake
+++ b/cmake/presets/mingw-cross.cmake
@@ -67,6 +67,7 @@ set(WIN_PACKAGES
REACTION
REAXFF
REPLICA
+ RHEO
RIGID
SHOCK
SMTBQ
diff --git a/cmake/presets/most.cmake b/cmake/presets/most.cmake
index d01642f94db..05282eebddd 100644
--- a/cmake/presets/most.cmake
+++ b/cmake/presets/most.cmake
@@ -60,6 +60,7 @@ set(ALL_PACKAGES
REACTION
REAXFF
REPLICA
+ RHEO
RIGID
SHOCK
SPH
diff --git a/cmake/presets/windows.cmake b/cmake/presets/windows.cmake
index 403d40efa44..71241c559c0 100644
--- a/cmake/presets/windows.cmake
+++ b/cmake/presets/windows.cmake
@@ -60,6 +60,7 @@ set(WIN_PACKAGES
REACTION
REAXFF
REPLICA
+ RHEO
RIGID
SHOCK
SMTBQ
diff --git a/doc/src/Build_extras.rst b/doc/src/Build_extras.rst
index ac7edc7678e..8465bea8290 100644
--- a/doc/src/Build_extras.rst
+++ b/doc/src/Build_extras.rst
@@ -2251,28 +2251,38 @@ verified to work in February 2020 with Quantum Espresso versions 6.3 to
RHEO package
------------
-To build with this package you must have the `GNU Scientific Library
-(GSL) ` installed in locations that
-are accessible in your environment. The GSL library should be at least
-version 2.7.
+This package depends on the BPM package.
.. tabs::
.. tab:: CMake build
- If CMake cannot find the GSL library or include files, you can set:
-
.. code-block:: bash
- -D GSL_ROOT_DIR=path # path to root of GSL installation
+ -D PKG_RHEO=yes # enable the package itself
+ -D PKG_BPM=yes # the RHEO package requires BPM
+ -D USE_INTERNAL_LINALG=value # prefer internal LAPACK if true
+
+ Some features in the RHEO package are dependent on code in the BPM
+ package so the latter one *must* be enabled as well.
+
+ The RHEO package also requires LAPACK (and BLAS) and CMake
+ can identify their locations and pass that info to the RHEO
+ build script. But on some systems this may cause problems when
+ linking or the dependency is not desired. By using the setting
+ ``-D USE_INTERNAL_LINALG=yes`` when running the CMake
+ configuration, you will select compiling and linking the bundled
+ linear algebra library and work around the limitations.
.. tab:: Traditional make
- LAMMPS will try to auto-detect the GSL compiler and linker flags
- from the corresponding ``pkg-config`` file (``gsl.pc``), otherwise
- you can edit the file ``lib/rheo/Makefile.lammps``
- to specify the paths and library names where indicated by comments.
- This must be done **before** the package is installed.
+ The RHEO package requires LAPACK (and BLAS) which can be either
+ a system provided library or the bundled "linalg" library. This
+ is a subset of LAPACK translated to C++. For that, one of the
+ provided ``Makefile.lammps.`` files needs to be copied
+ to ``Makefile.lammps`` and edited as needed. The default file
+ uses the bundled "linalg" library, which can be built by
+ ``make lib-linalg args='-m serial'`` in the ``src`` folder.
----------
diff --git a/doc/src/Build_settings.rst b/doc/src/Build_settings.rst
index 2787560be5b..e4a53ddee75 100644
--- a/doc/src/Build_settings.rst
+++ b/doc/src/Build_settings.rst
@@ -229,8 +229,7 @@ can be used with the Intel or GNU compiler (see the ``FFT_LIB`` setting
above).
The NVIDIA Performance Libraries (NVPL) FFT library is optimized for NVIDIA
-Grace Armv9.0 architecture. You can download it from
-`https://docs.nvidia.com/nvpl/`_.
+Grace Armv9.0 architecture. You can download it from https://docs.nvidia.com/nvpl/
The cuFFT and hipFFT FFT libraries are packaged with NVIDIA's CUDA and
AMD's HIP installations, respectively. These FFT libraries require the
diff --git a/doc/src/Intro_authors.rst b/doc/src/Intro_authors.rst
index 78c85064214..84470ba3a04 100644
--- a/doc/src/Intro_authors.rst
+++ b/doc/src/Intro_authors.rst
@@ -56,7 +56,7 @@ lammps.org". General questions about LAMMPS should be posted in the
- SNL
- jmgoff at sandia.gov
- machine learned potentials, QEq solvers, Python
- * - Megan McCarthy
+ * - Meg McCarthy
- SNL
- megmcca at sandia.gov
- alloys, micro-structure, machine learned potentials
@@ -67,7 +67,7 @@ lammps.org". General questions about LAMMPS should be posted in the
* - `Trung Nguyen `_
- U Chicago
- ndactrung at gmail.com
- - soft matter, GPU package
+ - soft matter, GPU package, DIELECTRIC package, regression testing
.. _rb: https://rbberger.github.io/
.. _gc: https://enthalpiste.fr/
diff --git a/doc/src/fix_adapt.rst b/doc/src/fix_adapt.rst
index 03aef12a6c5..a44ce8e780c 100644
--- a/doc/src/fix_adapt.rst
+++ b/doc/src/fix_adapt.rst
@@ -319,25 +319,32 @@ all types from 1 to :math:`N`. A leading asterisk means all types from
:math:`N` (inclusive). A middle asterisk means all types from m to n
(inclusive).
-Currently *bond* does not support bond_style hybrid nor bond_style
-hybrid/overlay as bond styles. The bond styles that currently work
-with fix_adapt are
-
-+------------------------------------+------------+------------+
-| :doc:`class2 ` | r0 | type bonds |
-+------------------------------------+------------+------------+
-| :doc:`fene ` | k,r0 | type bonds |
-+------------------------------------+------------+------------+
-| :doc:`fene/nm ` | k,r0 | type bonds |
-+------------------------------------+------------+------------+
-| :doc:`gromos ` | k,r0 | type bonds |
-+------------------------------------+------------+------------+
-| :doc:`harmonic ` | k,r0 | type bonds |
-+------------------------------------+------------+------------+
-| :doc:`morse ` | r0 | type bonds |
-+------------------------------------+------------+------------+
-| :doc:`nonlinear ` | epsilon,r0 | type bonds |
-+------------------------------------+------------+------------+
+If :doc:`bond_style hybrid ` is used, *bstyle* should be a
+sub-style name. The bond styles that currently work with fix adapt are:
+
++---------------------------------------------------+---------------------------+------------+
+| :doc:`class2 ` | r0 | type bonds |
++---------------------------------------------------+---------------------------+------------+
+| :doc:`fene ` | k,r0 | type bonds |
++---------------------------------------------------+---------------------------+------------+
+| :doc:`fene/expand ` | k,r0,epsilon,sigma,shift | type bonds |
++---------------------------------------------------+---------------------------+------------+
+| :doc:`fene/nm ` | k,r0 | type bonds |
++---------------------------------------------------+---------------------------+------------+
+| :doc:`gromos ` | k,r0 | type bonds |
++---------------------------------------------------+---------------------------+------------+
+| :doc:`harmonic ` | k,r0 | type bonds |
++---------------------------------------------------+---------------------------+------------+
+| :doc:`harmonic/shift ` | k,r0,r1 | type bonds |
++---------------------------------------------------+---------------------------+------------+
+| :doc:`harmonic/restrain ` | k | type bonds |
++---------------------------------------------------+---------------------------+------------+
+| :doc:`mm3 ` | k,r0 | type bonds |
++---------------------------------------------------+---------------------------+------------+
+| :doc:`morse ` | r0 | type bonds |
++---------------------------------------------------+---------------------------+------------+
+| :doc:`nonlinear ` | epsilon,r0 | type bonds |
++---------------------------------------------------+---------------------------+------------+
----------
@@ -357,15 +364,34 @@ all types from 1 to :math:`N`. A leading asterisk means all types from
:math:`N` (inclusive). A middle asterisk means all types from m to n
(inclusive).
-Currently *angle* does not support angle_style hybrid nor angle_style
-hybrid/overlay as angle styles. The angle styles that currently work
-with fix_adapt are
-
-+------------------------------------+----------+-------------+
-| :doc:`harmonic ` | k,theta0 | type angles |
-+------------------------------------+----------+-------------+
-| :doc:`cosine ` | k | type angles |
-+------------------------------------+----------+-------------+
+If :doc:`angle_style hybrid ` is used, *astyle* should be a
+sub-style name. The angle styles that currently work with fix adapt are:
+
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`harmonic ` | k,theta0 | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`charmm ` | k,theta0 | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`class2 ` | k2,k3,k4,theta0 | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`cosine ` | k | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`cosine/periodic ` | k,b,n | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`cosine/squared/restricted ` | k,theta0 | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`dipole ` | k,gamma0 | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`fourier ` | k,c0,c1,c2 | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`fourier/simple ` | k,c,n | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`mm3 ` | k,theta0 | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`quartic ` | k2,k3,k4,theta0 | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
+| :doc:`spica ` | k,theta0 | type angles |
++--------------------------------------------------------------------+-----------------+-------------+
Note that internally, theta0 is stored in radians, so the variable
this fix uses to reset theta0 needs to generate values in radians.
diff --git a/doc/src/region.rst b/doc/src/region.rst
index 9d2af01de1d..3a27c4b5ff2 100644
--- a/doc/src/region.rst
+++ b/doc/src/region.rst
@@ -18,13 +18,13 @@ Syntax
*delete* = no args
*block* args = xlo xhi ylo yhi zlo zhi
xlo,xhi,ylo,yhi,zlo,zhi = bounds of block in all dimensions (distance units)
- xlo,xhi,ylo,yhi,zlo,zhi can be a variable
+ xlo,xhi,ylo,yhi,zlo,zhi can be a variable (see below)
*cone* args = dim c1 c2 radlo radhi lo hi
dim = *x* or *y* or *z* = axis of cone
c1,c2 = coords of cone axis in other 2 dimensions (distance units)
radlo,radhi = cone radii at lo and hi end (distance units)
lo,hi = bounds of cone in dim (distance units)
- c1,c2,radlo,radhi,lo,hi can be a variable (see below)
+ c1,c2,radlo,radhi,lo,hi can be a variable (see below)
*cylinder* args = dim c1 c2 radius lo hi
dim = *x* or *y* or *z* = axis of cylinder
c1,c2 = coords of cylinder axis in other 2 dimensions (distance units)
@@ -38,6 +38,7 @@ Syntax
*plane* args = px py pz nx ny nz
px,py,pz = point on the plane (distance units)
nx,ny,nz = direction normal to plane (distance units)
+ px,py,pz can be a variable (see below)
*prism* args = xlo xhi ylo yhi zlo zhi xy xz yz
xlo,xhi,ylo,yhi,zlo,zhi = bounds of untilted prism (distance units)
xy = distance to tilt y in x direction (distance units)
@@ -166,7 +167,7 @@ extending in the y-direction from -5.0 to the upper box boundary.
.. versionadded:: 4May2022
-For style *ellipsoid*, an axis-aligned ellipsoid is defined. The
+For style *ellipsoid*, an axis-aligned ellipsoid is defined. The
ellipsoid has its center at (x,y,z) and is defined by 3 axis-aligned
vectors given by A = (a,0,0); B = (0,b,0); C = (0,0,c). Note that
although the ellipsoid is specified as axis-aligned it can be rotated
@@ -206,9 +207,10 @@ parameters a,b,c for style *ellipsoid*, can each be specified as an
equal-style :doc:`variable `. Likewise, for style *sphere*
and *ellipsoid* the x-, y-, and z- coordinates of the center of the
sphere/ellipsoid can be specified as an equal-style variable. And for
-style *cylinder* the two center positions c1 and c2 for the location
-of the cylinder axes can be specified as a equal-style variable. For style *cone*
-all properties can be defined via equal-style variables.
+style *cylinder* the two center positions c1 and c2 for the location of
+the cylinder axes can be specified as a equal-style variable. For style
+*cone* all properties can be defined via equal-style variables. For
+style *plane* the point can be defined via equal-style variables.
If the value is a variable, it should be specified as v_name, where
name is the variable name. In this case, the variable will be
diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt
index 3a714ffdd6f..70d6b4e3232 100644
--- a/doc/utils/sphinx-config/false_positives.txt
+++ b/doc/utils/sphinx-config/false_positives.txt
@@ -141,6 +141,7 @@ arg
arge
args
argv
+Armv
arrhenius
Arun
arXiv
diff --git a/lib/linalg/dlauu2.cpp b/lib/linalg/dlauu2.cpp
new file mode 100644
index 00000000000..d90a84798d5
--- /dev/null
+++ b/lib/linalg/dlauu2.cpp
@@ -0,0 +1,77 @@
+#ifdef __cplusplus
+extern "C" {
+#endif
+#include "lmp_f2c.h"
+static doublereal c_b7 = 1.;
+static integer c__1 = 1;
+int dlauu2_(char *uplo, integer *n, doublereal *a, integer *lda, integer *info, ftnlen uplo_len)
+{
+ integer a_dim1, a_offset, i__1, i__2, i__3;
+ integer i__;
+ doublereal aii;
+ extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, integer *);
+ extern int dscal_(integer *, doublereal *, doublereal *, integer *);
+ extern logical lsame_(char *, char *, ftnlen, ftnlen);
+ extern int dgemv_(char *, integer *, integer *, doublereal *, doublereal *, integer *,
+ doublereal *, integer *, doublereal *, doublereal *, integer *, ftnlen);
+ logical upper;
+ extern int xerbla_(char *, integer *, ftnlen);
+ a_dim1 = *lda;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+ *info = 0;
+ upper = lsame_(uplo, (char *)"U", (ftnlen)1, (ftnlen)1);
+ if (!upper && !lsame_(uplo, (char *)"L", (ftnlen)1, (ftnlen)1)) {
+ *info = -1;
+ } else if (*n < 0) {
+ *info = -2;
+ } else if (*lda < max(1, *n)) {
+ *info = -4;
+ }
+ if (*info != 0) {
+ i__1 = -(*info);
+ xerbla_((char *)"DLAUU2", &i__1, (ftnlen)6);
+ return 0;
+ }
+ if (*n == 0) {
+ return 0;
+ }
+ if (upper) {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ aii = a[i__ + i__ * a_dim1];
+ if (i__ < *n) {
+ i__2 = *n - i__ + 1;
+ a[i__ + i__ * a_dim1] =
+ ddot_(&i__2, &a[i__ + i__ * a_dim1], lda, &a[i__ + i__ * a_dim1], lda);
+ i__2 = i__ - 1;
+ i__3 = *n - i__;
+ dgemv_((char *)"No transpose", &i__2, &i__3, &c_b7, &a[(i__ + 1) * a_dim1 + 1], lda,
+ &a[i__ + (i__ + 1) * a_dim1], lda, &aii, &a[i__ * a_dim1 + 1], &c__1,
+ (ftnlen)12);
+ } else {
+ dscal_(&i__, &aii, &a[i__ * a_dim1 + 1], &c__1);
+ }
+ }
+ } else {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ aii = a[i__ + i__ * a_dim1];
+ if (i__ < *n) {
+ i__2 = *n - i__ + 1;
+ a[i__ + i__ * a_dim1] =
+ ddot_(&i__2, &a[i__ + i__ * a_dim1], &c__1, &a[i__ + i__ * a_dim1], &c__1);
+ i__2 = *n - i__;
+ i__3 = i__ - 1;
+ dgemv_((char *)"Transpose", &i__2, &i__3, &c_b7, &a[i__ + 1 + a_dim1], lda,
+ &a[i__ + 1 + i__ * a_dim1], &c__1, &aii, &a[i__ + a_dim1], lda, (ftnlen)9);
+ } else {
+ dscal_(&i__, &aii, &a[i__ + a_dim1], lda);
+ }
+ }
+ }
+ return 0;
+}
+#ifdef __cplusplus
+}
+#endif
diff --git a/lib/linalg/dlauum.cpp b/lib/linalg/dlauum.cpp
new file mode 100644
index 00000000000..632bd4ba85d
--- /dev/null
+++ b/lib/linalg/dlauum.cpp
@@ -0,0 +1,101 @@
+#ifdef __cplusplus
+extern "C" {
+#endif
+#include "lmp_f2c.h"
+static integer c__1 = 1;
+static integer c_n1 = -1;
+static doublereal c_b15 = 1.;
+int dlauum_(char *uplo, integer *n, doublereal *a, integer *lda, integer *info, ftnlen uplo_len)
+{
+ integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
+ integer i__, ib, nb;
+ extern int dgemm_(char *, char *, integer *, integer *, integer *, doublereal *, doublereal *,
+ integer *, doublereal *, integer *, doublereal *, doublereal *, integer *,
+ ftnlen, ftnlen);
+ extern logical lsame_(char *, char *, ftnlen, ftnlen);
+ extern int dtrmm_(char *, char *, char *, char *, integer *, integer *, doublereal *,
+ doublereal *, integer *, doublereal *, integer *, ftnlen, ftnlen, ftnlen,
+ ftnlen);
+ logical upper;
+ extern int dsyrk_(char *, char *, integer *, integer *, doublereal *, doublereal *, integer *,
+ doublereal *, doublereal *, integer *, ftnlen, ftnlen),
+ dlauu2_(char *, integer *, doublereal *, integer *, integer *, ftnlen),
+ xerbla_(char *, integer *, ftnlen);
+ extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *,
+ ftnlen, ftnlen);
+ a_dim1 = *lda;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+ *info = 0;
+ upper = lsame_(uplo, (char *)"U", (ftnlen)1, (ftnlen)1);
+ if (!upper && !lsame_(uplo, (char *)"L", (ftnlen)1, (ftnlen)1)) {
+ *info = -1;
+ } else if (*n < 0) {
+ *info = -2;
+ } else if (*lda < max(1, *n)) {
+ *info = -4;
+ }
+ if (*info != 0) {
+ i__1 = -(*info);
+ xerbla_((char *)"DLAUUM", &i__1, (ftnlen)6);
+ return 0;
+ }
+ if (*n == 0) {
+ return 0;
+ }
+ nb = ilaenv_(&c__1, (char *)"DLAUUM", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
+ if (nb <= 1 || nb >= *n) {
+ dlauu2_(uplo, n, &a[a_offset], lda, info, (ftnlen)1);
+ } else {
+ if (upper) {
+ i__1 = *n;
+ i__2 = nb;
+ for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
+ i__3 = nb, i__4 = *n - i__ + 1;
+ ib = min(i__3, i__4);
+ i__3 = i__ - 1;
+ dtrmm_((char *)"Right", (char *)"Upper", (char *)"Transpose", (char *)"Non-unit", &i__3, &ib, &c_b15,
+ &a[i__ + i__ * a_dim1], lda, &a[i__ * a_dim1 + 1], lda, (ftnlen)5, (ftnlen)5,
+ (ftnlen)9, (ftnlen)8);
+ dlauu2_((char *)"Upper", &ib, &a[i__ + i__ * a_dim1], lda, info, (ftnlen)5);
+ if (i__ + ib <= *n) {
+ i__3 = i__ - 1;
+ i__4 = *n - i__ - ib + 1;
+ dgemm_((char *)"No transpose", (char *)"Transpose", &i__3, &ib, &i__4, &c_b15,
+ &a[(i__ + ib) * a_dim1 + 1], lda, &a[i__ + (i__ + ib) * a_dim1], lda,
+ &c_b15, &a[i__ * a_dim1 + 1], lda, (ftnlen)12, (ftnlen)9);
+ i__3 = *n - i__ - ib + 1;
+ dsyrk_((char *)"Upper", (char *)"No transpose", &ib, &i__3, &c_b15,
+ &a[i__ + (i__ + ib) * a_dim1], lda, &c_b15, &a[i__ + i__ * a_dim1], lda,
+ (ftnlen)5, (ftnlen)12);
+ }
+ }
+ } else {
+ i__2 = *n;
+ i__1 = nb;
+ for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__1) {
+ i__3 = nb, i__4 = *n - i__ + 1;
+ ib = min(i__3, i__4);
+ i__3 = i__ - 1;
+ dtrmm_((char *)"Left", (char *)"Lower", (char *)"Transpose", (char *)"Non-unit", &ib, &i__3, &c_b15,
+ &a[i__ + i__ * a_dim1], lda, &a[i__ + a_dim1], lda, (ftnlen)4, (ftnlen)5,
+ (ftnlen)9, (ftnlen)8);
+ dlauu2_((char *)"Lower", &ib, &a[i__ + i__ * a_dim1], lda, info, (ftnlen)5);
+ if (i__ + ib <= *n) {
+ i__3 = i__ - 1;
+ i__4 = *n - i__ - ib + 1;
+ dgemm_((char *)"Transpose", (char *)"No transpose", &ib, &i__3, &i__4, &c_b15,
+ &a[i__ + ib + i__ * a_dim1], lda, &a[i__ + ib + a_dim1], lda, &c_b15,
+ &a[i__ + a_dim1], lda, (ftnlen)9, (ftnlen)12);
+ i__3 = *n - i__ - ib + 1;
+ dsyrk_((char *)"Lower", (char *)"Transpose", &ib, &i__3, &c_b15, &a[i__ + ib + i__ * a_dim1],
+ lda, &c_b15, &a[i__ + i__ * a_dim1], lda, (ftnlen)5, (ftnlen)9);
+ }
+ }
+ }
+ }
+ return 0;
+}
+#ifdef __cplusplus
+}
+#endif
diff --git a/lib/linalg/dpotri.cpp b/lib/linalg/dpotri.cpp
new file mode 100644
index 00000000000..9c0a609e1b3
--- /dev/null
+++ b/lib/linalg/dpotri.cpp
@@ -0,0 +1,40 @@
+#ifdef __cplusplus
+extern "C" {
+#endif
+#include "lmp_f2c.h"
+int dpotri_(char *uplo, integer *n, doublereal *a, integer *lda, integer *info, ftnlen uplo_len)
+{
+ integer a_dim1, a_offset, i__1;
+ extern logical lsame_(char *, char *, ftnlen, ftnlen);
+ extern int xerbla_(char *, integer *, ftnlen),
+ dlauum_(char *, integer *, doublereal *, integer *, integer *, ftnlen),
+ dtrtri_(char *, char *, integer *, doublereal *, integer *, integer *, ftnlen, ftnlen);
+ a_dim1 = *lda;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+ *info = 0;
+ if (!lsame_(uplo, (char *)"U", (ftnlen)1, (ftnlen)1) && !lsame_(uplo, (char *)"L", (ftnlen)1, (ftnlen)1)) {
+ *info = -1;
+ } else if (*n < 0) {
+ *info = -2;
+ } else if (*lda < max(1, *n)) {
+ *info = -4;
+ }
+ if (*info != 0) {
+ i__1 = -(*info);
+ xerbla_((char *)"DPOTRI", &i__1, (ftnlen)6);
+ return 0;
+ }
+ if (*n == 0) {
+ return 0;
+ }
+ dtrtri_(uplo, (char *)"Non-unit", n, &a[a_offset], lda, info, (ftnlen)1, (ftnlen)8);
+ if (*info > 0) {
+ return 0;
+ }
+ dlauum_(uplo, n, &a[a_offset], lda, info, (ftnlen)1);
+ return 0;
+}
+#ifdef __cplusplus
+}
+#endif
diff --git a/lib/rheo/Makefile.lammps b/lib/rheo/Makefile.lammps
index ec587403700..5785f8978ba 100644
--- a/lib/rheo/Makefile.lammps
+++ b/lib/rheo/Makefile.lammps
@@ -1,14 +1,5 @@
-# Settings that the LAMMPS build will import when this package is installed
+# Settings that the LAMMPS build will import when this package library is used
-ifeq ($(strip $(shell pkg-config --version)),)
- # manual configuration w/o pkg-config/pkgconf
- # change this to -I/path/to/your/lib/gsl/include/
- rheo_SYSINC = -I../../lib/rheo/gsl/include/
-
- # change this to -L/path/to/your/lib/gsl/lib/
- rheo_SYSLIB = -L../../lib/rheo/gsl/lib/ -lgsl -lgslcblas
-else
- # autodetect GSL settings from pkg-config/pkgconf
- rheo_SYSINC = $(shell pkg-config --cflags gsl)
- rheo_SYSLIB = $(shell pkg-config --libs gsl)
-endif
+rheo_SYSINC =
+rheo_SYSLIB = -llinalg
+rheo_SYSPATH = -L../../lib/linalg$(LIBOBJDIR)
diff --git a/lib/rheo/Makefile.lammps.empty b/lib/rheo/Makefile.lammps.empty
new file mode 100644
index 00000000000..f71390299c5
--- /dev/null
+++ b/lib/rheo/Makefile.lammps.empty
@@ -0,0 +1,5 @@
+# Settings that the LAMMPS build will import when this package library is used
+
+rheo_SYSINC =
+rheo_SYSLIB =
+rheo_SYSPATH =
diff --git a/lib/rheo/Makefile.lammps.installed b/lib/rheo/Makefile.lammps.installed
new file mode 100644
index 00000000000..8900470077f
--- /dev/null
+++ b/lib/rheo/Makefile.lammps.installed
@@ -0,0 +1,5 @@
+# Settings that the LAMMPS build will import when this package library is used
+
+rheo_SYSINC =
+rheo_SYSLIB = -lblas -llapack
+rheo_SYSPATH =
diff --git a/lib/rheo/Makefile.lammps.linalg b/lib/rheo/Makefile.lammps.linalg
new file mode 100644
index 00000000000..5785f8978ba
--- /dev/null
+++ b/lib/rheo/Makefile.lammps.linalg
@@ -0,0 +1,5 @@
+# Settings that the LAMMPS build will import when this package library is used
+
+rheo_SYSINC =
+rheo_SYSLIB = -llinalg
+rheo_SYSPATH = -L../../lib/linalg$(LIBOBJDIR)
diff --git a/lib/rheo/README b/lib/rheo/README
index ae421b6e809..fe082797f18 100644
--- a/lib/rheo/README
+++ b/lib/rheo/README
@@ -1,7 +1,5 @@
-This directory has a Makefile.lammps file with settings that allows LAMMPS to
-dynamically link to the GSL library. This is required to use the RHEO package
-in a LAMMPS input script. If you have the pkg-config command available, it
-will automatically import the GSL settings. Otherwise they will have to be
-added manually.
-
-See the header of Makefile.lammps for more info.
+This directory has multiple Makefile.lammps variant files with settings that
+allows LAMMPS to link with a BLAS/LAPACK or compatible library or the bundled
+linalg library (which is subset of BLAS/LAPACK). Copy the suitable file
+to Makefile.lammps and edit, if needed.
+This is required to use the RHEO package in a LAMMPS input script.
diff --git a/src/CG-SPICA/angle_spica.cpp b/src/CG-SPICA/angle_spica.cpp
index e315e20f13b..913428cd9b0 100644
--- a/src/CG-SPICA/angle_spica.cpp
+++ b/src/CG-SPICA/angle_spica.cpp
@@ -522,3 +522,15 @@ double AngleSPICA::single(int type, int i1, int i2, int i3)
double tk = k[type] * dtheta;
return tk*dtheta + e13;
}
+
+/* ----------------------------------------------------------------------
+ return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleSPICA::extract(const char *str, int &dim)
+{
+ dim = 1;
+ if (strcmp(str, "k") == 0) return (void *) k;
+ if (strcmp(str, "theta0") == 0) return (void *) theta0;
+ return nullptr;
+}
diff --git a/src/CG-SPICA/angle_spica.h b/src/CG-SPICA/angle_spica.h
index 539512c0e98..5e590ba7a05 100644
--- a/src/CG-SPICA/angle_spica.h
+++ b/src/CG-SPICA/angle_spica.h
@@ -37,6 +37,7 @@ class AngleSPICA : public Angle {
void read_restart(FILE *) override;
void write_data(FILE *) override;
double single(int, int, int, int) override;
+ void *extract(const char *, int &) override;
protected:
double *k, *theta0;
diff --git a/src/CLASS2/angle_class2.cpp b/src/CLASS2/angle_class2.cpp
index 118179ad910..5000f9f629e 100644
--- a/src/CLASS2/angle_class2.cpp
+++ b/src/CLASS2/angle_class2.cpp
@@ -467,3 +467,17 @@ double AngleClass2::single(int type, int i1, int i2, int i3)
return energy;
}
+
+/* ----------------------------------------------------------------------
+ return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleClass2::extract(const char *str, int &dim)
+{
+ dim = 1;
+ if (strcmp(str, "k2") == 0) return (void *) k2;
+ if (strcmp(str, "k3") == 0) return (void *) k3;
+ if (strcmp(str, "k4") == 0) return (void *) k4;
+ if (strcmp(str, "theta0") == 0) return (void *) theta0;
+ return nullptr;
+}
diff --git a/src/CLASS2/angle_class2.h b/src/CLASS2/angle_class2.h
index f5fbd62b575..4ed6f344ae5 100644
--- a/src/CLASS2/angle_class2.h
+++ b/src/CLASS2/angle_class2.h
@@ -35,6 +35,7 @@ class AngleClass2 : public Angle {
void read_restart(FILE *) override;
void write_data(FILE *) override;
double single(int, int, int, int) override;
+ void *extract(const char *, int &) override;
protected:
double *theta0, *k2, *k3, *k4;
diff --git a/src/DIPOLE/angle_dipole.cpp b/src/DIPOLE/angle_dipole.cpp
index 6ad4a0fb4cd..a66f3e10426 100644
--- a/src/DIPOLE/angle_dipole.cpp
+++ b/src/DIPOLE/angle_dipole.cpp
@@ -263,3 +263,15 @@ double AngleDipole::single(int type, int iRef, int iDip, int /*iDummy*/)
return kdg * deltaGamma; // energy
}
+
+/* ----------------------------------------------------------------------
+ return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleDipole::extract(const char *str, int &dim)
+{
+ dim = 1;
+ if (strcmp(str, "k") == 0) return (void *) k;
+ if (strcmp(str, "gamma0") == 0) return (void *) gamma0;
+ return nullptr;
+}
diff --git a/src/DIPOLE/angle_dipole.h b/src/DIPOLE/angle_dipole.h
index 2e557226731..de0f958f989 100644
--- a/src/DIPOLE/angle_dipole.h
+++ b/src/DIPOLE/angle_dipole.h
@@ -36,6 +36,7 @@ class AngleDipole : public Angle {
void read_restart(FILE *) override;
void write_data(FILE *) override;
double single(int, int, int, int) override;
+ void *extract(const char *, int &) override;
protected:
double *k, *gamma0;
diff --git a/src/EXTRA-FIX/fix_ave_correlate_long.cpp b/src/EXTRA-FIX/fix_ave_correlate_long.cpp
index 738ae3ae4cd..abb1ad87dea 100644
--- a/src/EXTRA-FIX/fix_ave_correlate_long.cpp
+++ b/src/EXTRA-FIX/fix_ave_correlate_long.cpp
@@ -166,12 +166,12 @@ FixAveCorrelateLong::FixAveCorrelateLong(LAMMPS *lmp, int narg, char **arg) :
overwrite = 1;
iarg += 1;
} else if (strcmp(arg[iarg], "title1") == 0) {
- if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate/long title1", error);
+ if (iarg + 2 > nargnew) utils::missing_cmd_args(FLERR, "fix ave/correlate/long title1", error);
delete[] title1;
title1 = utils::strdup(arg[iarg + 1]);
iarg += 2;
} else if (strcmp(arg[iarg], "title2") == 0) {
- if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate/long title2", error);
+ if (iarg + 2 > nargnew) utils::missing_cmd_args(FLERR, "fix ave/correlate/long title2", error);
delete[] title2;
title2 = utils::strdup(arg[iarg + 1]);
iarg += 2;
diff --git a/src/EXTRA-MOLECULE/angle_cosine_periodic.cpp b/src/EXTRA-MOLECULE/angle_cosine_periodic.cpp
index 0b2a6d336d0..6aae8d3ff48 100644
--- a/src/EXTRA-MOLECULE/angle_cosine_periodic.cpp
+++ b/src/EXTRA-MOLECULE/angle_cosine_periodic.cpp
@@ -336,3 +336,16 @@ void AngleCosinePeriodic::born_matrix(int type, int i1, int i2, int i3, double &
du = prefactor * sin(m_angle) / s;
du2 = prefactor * (c * sin(m_angle) - s * cos(m_angle) * multiplicity[type]) / (s * s * s);
}
+
+/* ----------------------------------------------------------------------
+ return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleCosinePeriodic::extract(const char *str, int &dim)
+{
+ dim = 1;
+ if (strcmp(str, "k") == 0) return (void *) k;
+ if (strcmp(str, "b") == 0) return (void *) b;
+ if (strcmp(str, "multiplicity") == 0) return (void *) multiplicity;
+ return nullptr;
+}
diff --git a/src/EXTRA-MOLECULE/angle_cosine_periodic.h b/src/EXTRA-MOLECULE/angle_cosine_periodic.h
index f04ed047841..f63029919ee 100644
--- a/src/EXTRA-MOLECULE/angle_cosine_periodic.h
+++ b/src/EXTRA-MOLECULE/angle_cosine_periodic.h
@@ -36,6 +36,7 @@ class AngleCosinePeriodic : public Angle {
void write_data(FILE *) override;
double single(int, int, int, int) override;
void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override;
+ void *extract(const char *, int &) override;
protected:
double *k;
diff --git a/src/EXTRA-MOLECULE/angle_cosine_squared_restricted.cpp b/src/EXTRA-MOLECULE/angle_cosine_squared_restricted.cpp
index 2da31ef8932..c6d78ea1331 100644
--- a/src/EXTRA-MOLECULE/angle_cosine_squared_restricted.cpp
+++ b/src/EXTRA-MOLECULE/angle_cosine_squared_restricted.cpp
@@ -296,3 +296,15 @@ void AngleCosineSquaredRestricted::born_matrix(int type, int i1, int i2, int i3,
du2 = 2 * k[type] * numerator / denominator;
}
+
+/* ----------------------------------------------------------------------
+ return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleCosineSquaredRestricted::extract(const char *str, int &dim)
+{
+ dim = 1;
+ if (strcmp(str, "k") == 0) return (void *) k;
+ if (strcmp(str, "theta0") == 0) return (void *) theta0;
+ return nullptr;
+}
diff --git a/src/EXTRA-MOLECULE/angle_cosine_squared_restricted.h b/src/EXTRA-MOLECULE/angle_cosine_squared_restricted.h
index 674252b7d0f..b38b6bc4bde 100644
--- a/src/EXTRA-MOLECULE/angle_cosine_squared_restricted.h
+++ b/src/EXTRA-MOLECULE/angle_cosine_squared_restricted.h
@@ -36,6 +36,7 @@ class AngleCosineSquaredRestricted : public Angle {
void write_data(FILE *) override;
double single(int, int, int, int) override;
void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override;
+ void *extract(const char *, int &) override;
protected:
double *k, *theta0;
diff --git a/src/EXTRA-MOLECULE/angle_fourier.cpp b/src/EXTRA-MOLECULE/angle_fourier.cpp
index da1667c06f6..abcda6d0361 100644
--- a/src/EXTRA-MOLECULE/angle_fourier.cpp
+++ b/src/EXTRA-MOLECULE/angle_fourier.cpp
@@ -309,3 +309,16 @@ void AngleFourier::born_matrix(int type, int i1, int i2, int i3, double &du, dou
du = k[type] * (C1[type] + 4 * C2[type] * c);
}
+/* ----------------------------------------------------------------------
+ return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleFourier::extract(const char *str, int &dim)
+{
+ dim = 1;
+ if (strcmp(str, "k") == 0) return (void *) k;
+ if (strcmp(str, "C0") == 0) return (void *) C0;
+ if (strcmp(str, "C1") == 0) return (void *) C1;
+ if (strcmp(str, "C2") == 0) return (void *) C2;
+ return nullptr;
+}
diff --git a/src/EXTRA-MOLECULE/angle_fourier.h b/src/EXTRA-MOLECULE/angle_fourier.h
index 8fa5d14b266..c0e30c8e1a4 100644
--- a/src/EXTRA-MOLECULE/angle_fourier.h
+++ b/src/EXTRA-MOLECULE/angle_fourier.h
@@ -36,6 +36,7 @@ class AngleFourier : public Angle {
void write_data(FILE *) override;
double single(int, int, int, int) override;
void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override;
+ void *extract(const char *, int &) override;
protected:
double *k, *C0, *C1, *C2;
diff --git a/src/EXTRA-MOLECULE/angle_fourier_simple.cpp b/src/EXTRA-MOLECULE/angle_fourier_simple.cpp
index 6de7956ffa7..143a008039d 100644
--- a/src/EXTRA-MOLECULE/angle_fourier_simple.cpp
+++ b/src/EXTRA-MOLECULE/angle_fourier_simple.cpp
@@ -316,3 +316,16 @@ void AngleFourierSimple::born_matrix(int type, int i1, int i2, int i3, double &d
du2 = k[type] * C[type] * N[type] * (cos(theta) * sin(N[type] * theta)
- N[type] * sin(theta) * cos(N[type] * theta)) / pow(sin(theta),3);
}
+
+/* ----------------------------------------------------------------------
+ return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleFourierSimple::extract(const char *str, int &dim)
+{
+ dim = 1;
+ if (strcmp(str, "k") == 0) return (void *) k;
+ if (strcmp(str, "C") == 0) return (void *) C;
+ if (strcmp(str, "N") == 0) return (void *) N;
+ return nullptr;
+}
diff --git a/src/EXTRA-MOLECULE/angle_fourier_simple.h b/src/EXTRA-MOLECULE/angle_fourier_simple.h
index 3296ba60673..d37b3a83a87 100644
--- a/src/EXTRA-MOLECULE/angle_fourier_simple.h
+++ b/src/EXTRA-MOLECULE/angle_fourier_simple.h
@@ -36,6 +36,7 @@ class AngleFourierSimple : public Angle {
void write_data(FILE *) override;
double single(int, int, int, int) override;
void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override;
+ void *extract(const char *, int &) override;
protected:
double *k, *C, *N;
diff --git a/src/EXTRA-MOLECULE/angle_quartic.cpp b/src/EXTRA-MOLECULE/angle_quartic.cpp
index aade6b4534d..616c81c749c 100644
--- a/src/EXTRA-MOLECULE/angle_quartic.cpp
+++ b/src/EXTRA-MOLECULE/angle_quartic.cpp
@@ -325,3 +325,17 @@ void AngleQuartic::born_matrix(int type, int i1, int i2, int i3, double &du, dou
du2 = (2.0 * k2[type] + 6.0 * k3[type] * dtheta + 12.0 * k4[type] * dtheta2) / (s*s) -
(2.0 * k2[type] * dtheta + 3.0 * k3[type] * dtheta2 + 4.0 * k4[type] * dtheta3) * c / (s*s*s);
}
+
+/* ----------------------------------------------------------------------
+ return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleQuartic::extract(const char *str, int &dim)
+{
+ dim = 1;
+ if (strcmp(str, "k2") == 0) return (void *) k2;
+ if (strcmp(str, "k3") == 0) return (void *) k3;
+ if (strcmp(str, "k4") == 0) return (void *) k4;
+ if (strcmp(str, "theta0") == 0) return (void *) theta0;
+ return nullptr;
+}
diff --git a/src/EXTRA-MOLECULE/angle_quartic.h b/src/EXTRA-MOLECULE/angle_quartic.h
index 7de51b24d1f..3ff7f6f3e4a 100644
--- a/src/EXTRA-MOLECULE/angle_quartic.h
+++ b/src/EXTRA-MOLECULE/angle_quartic.h
@@ -36,6 +36,7 @@ class AngleQuartic : public Angle {
void write_data(FILE *) override;
double single(int, int, int, int) override;
void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override;
+ void *extract(const char *, int &) override;
protected:
double *k2, *k3, *k4, *theta0;
diff --git a/src/EXTRA-MOLECULE/bond_harmonic_shift.cpp b/src/EXTRA-MOLECULE/bond_harmonic_shift.cpp
index bd106c85671..6c87d47f5eb 100644
--- a/src/EXTRA-MOLECULE/bond_harmonic_shift.cpp
+++ b/src/EXTRA-MOLECULE/bond_harmonic_shift.cpp
@@ -228,3 +228,16 @@ void BondHarmonicShift::born_matrix(int type, double rsq, int /*i*/, int /*j*/,
du2 = 2 * k[type];
if (r > 0.0) du = du2 * dr;
}
+
+/* ----------------------------------------------------------------------
+ return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *BondHarmonicShift::extract(const char *str, int &dim)
+{
+ dim = 1;
+ if (strcmp(str, "k") == 0) return (void *) k;
+ if (strcmp(str, "r0") == 0) return (void *) r0;
+ if (strcmp(str, "r1") == 0) return (void *) r1;
+ return nullptr;
+}
diff --git a/src/EXTRA-MOLECULE/bond_harmonic_shift.h b/src/EXTRA-MOLECULE/bond_harmonic_shift.h
index 922f1ba00dc..68afc57bf74 100644
--- a/src/EXTRA-MOLECULE/bond_harmonic_shift.h
+++ b/src/EXTRA-MOLECULE/bond_harmonic_shift.h
@@ -36,6 +36,7 @@ class BondHarmonicShift : public Bond {
void write_data(FILE *) override;
double single(int, double, int, int, double &) override;
void born_matrix(int, double, int, int, double &, double &) override;
+ void *extract(const char *, int &);
protected:
double *k, *r0, *r1;
diff --git a/src/ML-IAP/mliap_data.cpp b/src/ML-IAP/mliap_data.cpp
index 5d847ee25e7..1ca57e27457 100644
--- a/src/ML-IAP/mliap_data.cpp
+++ b/src/ML-IAP/mliap_data.cpp
@@ -185,7 +185,6 @@ void MLIAPData::generate_neighdata(NeighList *list_in, int eflag_in, int vflag_i
int jtype = type[j];
const int jelem = map[jtype];
- lmp_firstneigh[ii][jj] = firstneigh[i][jj];
if (rsq < descriptor->cutsq[ielem][jelem]) {
pair_i[ij] = i;
jatoms[ij] = j;
@@ -193,6 +192,7 @@ void MLIAPData::generate_neighdata(NeighList *list_in, int eflag_in, int vflag_i
rij[ij][0] = delx;
rij[ij][1] = dely;
rij[ij][2] = delz;
+ lmp_firstneigh[ii][ninside] = firstneigh[i][jj];
ij++;
ninside++;
}
@@ -228,6 +228,7 @@ void MLIAPData::grow_neigharrays()
memory->grow(ielems, natomneigh, "MLIAPData:ielems");
memory->grow(itypes, natomneigh, "MLIAPData:itypes");
memory->grow(numneighs, natomneigh, "MLIAPData:numneighs");
+ memory->grow(lmp_firstneigh, natomneigh, nneigh_max, "MLIAPData:lmp_firstneigh");
natomneigh_max = natomneigh;
}
diff --git a/src/MOLECULE/angle_charmm.cpp b/src/MOLECULE/angle_charmm.cpp
index 1b66260c55a..11b5abd6990 100644
--- a/src/MOLECULE/angle_charmm.cpp
+++ b/src/MOLECULE/angle_charmm.cpp
@@ -309,3 +309,15 @@ double AngleCharmm::single(int type, int i1, int i2, int i3)
return (tk * dtheta + rk * dr);
}
+
+/* ----------------------------------------------------------------------
+ return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleCharmm::extract(const char *str, int &dim)
+{
+ dim = 1;
+ if (strcmp(str, "k") == 0) return (void *) k;
+ if (strcmp(str, "theta0") == 0) return (void *) theta0;
+ return nullptr;
+}
diff --git a/src/MOLECULE/angle_charmm.h b/src/MOLECULE/angle_charmm.h
index 2a77ac78641..cc4095e90f1 100644
--- a/src/MOLECULE/angle_charmm.h
+++ b/src/MOLECULE/angle_charmm.h
@@ -35,6 +35,7 @@ class AngleCharmm : public Angle {
void read_restart(FILE *) override;
void write_data(FILE *) override;
double single(int, int, int, int) override;
+ void *extract(const char *, int &) override;
protected:
double *k, *theta0, *k_ub, *r_ub;
diff --git a/src/MOLECULE/bond_fene_expand.cpp b/src/MOLECULE/bond_fene_expand.cpp
index c7821b18267..e115596eb1c 100644
--- a/src/MOLECULE/bond_fene_expand.cpp
+++ b/src/MOLECULE/bond_fene_expand.cpp
@@ -273,3 +273,18 @@ double BondFENEExpand::single(int type, double rsq, int /*i*/, int /*j*/, double
return eng;
}
+
+/* ----------------------------------------------------------------------
+ return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *BondFENEExpand::extract(const char *str, int &dim)
+{
+ dim = 1;
+ if (strcmp(str, "k") == 0) return (void *) k;
+ if (strcmp(str, "r0") == 0) return (void *) r0;
+ if (strcmp(str, "epsilon") == 0) return (void *) epsilon;
+ if (strcmp(str, "sigma") == 0) return (void *) sigma;
+ if (strcmp(str, "shift") == 0) return (void *) shift;
+ return nullptr;
+}
diff --git a/src/MOLECULE/bond_fene_expand.h b/src/MOLECULE/bond_fene_expand.h
index cdce710ea14..13524b0972c 100644
--- a/src/MOLECULE/bond_fene_expand.h
+++ b/src/MOLECULE/bond_fene_expand.h
@@ -36,6 +36,7 @@ class BondFENEExpand : public Bond {
void read_restart(FILE *) override;
void write_data(FILE *) override;
double single(int, double, int, int, double &) override;
+ void *extract(const char *, int &) override;
protected:
double *k, *r0, *epsilon, *sigma, *shift;
diff --git a/src/RHEO/Install.sh b/src/RHEO/Install.sh
index e34ca3a5554..07a439f44bc 100644
--- a/src/RHEO/Install.sh
+++ b/src/RHEO/Install.sh
@@ -47,6 +47,7 @@ if (test $1 = 1) then
sed -i -e 's/[^ \t]*rheo[^ \t]* //' ../Makefile.package
sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(rheo_SYSINC) |' ../Makefile.package
sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(rheo_SYSLIB) |' ../Makefile.package
+ sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(rheo_SYSPATH) |' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then
diff --git a/src/RHEO/README b/src/RHEO/README
index 4b6f2a162a4..15b642442a6 100644
--- a/src/RHEO/README
+++ b/src/RHEO/README
@@ -3,8 +3,9 @@ multiphase fluid systems. The authors include Joel Clemmer (Sandia), Thomas
O'Connor (Carnegie Mellon), and Eric Palermo (Carnegie Mellon).
Bond style rheo/shell, compute style rheo/property/atom, and fix style
-rheo/temperature all depend on the BPM package.
+rheo/temperature depend on the BPM package, so it is required to install
+the BPM package with RHEO.
-This package requires the GNU scientific library (GSL). We recommend version
-2.7 or later. To build this package, one must first separately install GSL in
-a location that can be found by your environment.
+This package requires the BLAS/LAPACK. This can be either a seperate installation
+or you can use the bundled "linalg" library. Please see the LAMMPS manual at
+https://docs.lammps.org/Build_extras.html#rheo for details.
diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp
index dd05901b32e..453f2b90284 100644
--- a/src/RHEO/compute_rheo_kernel.cpp
+++ b/src/RHEO/compute_rheo_kernel.cpp
@@ -37,10 +37,6 @@
#include "utils.h"
#include
-#include
-#include
-#include
-#include
using namespace LAMMPS_NS;
using namespace RHEO_NS;
@@ -50,6 +46,13 @@ using namespace MathExtra;
// max value of Mdim 1 + dim + dim * (dim + 1) / 2 with dim = 3
static constexpr int MAX_MDIM = 12;
+// declare LAPACK functions
+
+extern "C" {
+ void dpotrf_(const char *uplo, const int *n, double *a, const int *lda, int *info);
+ void dpotri_(const char *uplo, const int *n, double *a, const int *lda, int *info);
+}
+
/* ---------------------------------------------------------------------- */
ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) :
@@ -89,7 +92,7 @@ ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) :
comm_forward_save = comm_forward;
corrections_calculated = 0;
- gsl_error_flag = 0;
+ lapack_error_flag = 0;
}
/* ---------------------------------------------------------------------- */
@@ -156,9 +159,9 @@ void ComputeRHEOKernel::init_list(int /*id*/, NeighList *ptr)
int ComputeRHEOKernel::check_corrections(int i)
{
- // Skip if there were gsl errors for this atom
- if (gsl_error_flag)
- if (gsl_error_tags.find(atom->tag[i]) != gsl_error_tags.end()) return 0;
+ // Skip if there were lapack errors for this atom
+ if (lapack_error_flag)
+ if (lapack_error_tags.find(atom->tag[i]) != lapack_error_tags.end()) return 0;
// Skip if undercoordinated
if (coordination[i] < zmin) return 0;
@@ -558,19 +561,15 @@ void ComputeRHEOKernel::calc_dw_rk2(int i, double delx, double dely, double delz
void ComputeRHEOKernel::compute_peratom()
{
- gsl_error_flag = 0;
- gsl_error_tags.clear();
+ lapack_error_flag = 0;
+ lapack_error_tags.clear();
if (kernel_style == QUINTIC) return;
corrections_calculated = 1;
- int i, j, ii, jj, inum, jnum, a, b, gsl_error;
+ int i, j, ii, jj, inum, jnum, a, b, lapack_error;
double xtmp, ytmp, ztmp, r, rsq, w, vj, rhoj;
double dx[3];
- gsl_matrix_view gM;
-
- // Turn off GSL error handler, revert RK to Quintic when insufficient neighbors
- gsl_set_error_handler_off();
double **x = atom->x;
int *type = atom->type;
@@ -633,7 +632,7 @@ void ComputeRHEOKernel::compute_peratom()
}
} else if (correction_order > 0) {
- // Moment matrix M and polynomial basis vector cut (1d for gsl compatibility)
+ // Moment matrix M and polynomial basis vector cut (1d for LAPACK compatibility)
double H[MAX_MDIM], M[MAX_MDIM * MAX_MDIM];
for (ii = 0; ii < inum; ii++) {
@@ -647,7 +646,9 @@ void ComputeRHEOKernel::compute_peratom()
// Zero upper-triangle M and cut (will be symmetric):
for (a = 0; a < Mdim; a++) {
- for (b = a; b < Mdim; b++) { M[a * Mdim + b] = 0; }
+ for (b = a; b < Mdim; b++) {
+ M[a * Mdim + b] = 0;
+ }
}
for (jj = 0; jj < jnum; jj++) {
@@ -700,37 +701,50 @@ void ComputeRHEOKernel::compute_peratom()
// Populate the upper triangle
for (a = 0; a < Mdim; a++) {
- for (b = a; b < Mdim; b++) { M[a * Mdim + b] += H[a] * H[b] * w * vj; }
+ for (b = a; b < Mdim; b++) {
+ M[a * Mdim + b] += H[a] * H[b] * w * vj;
+ }
}
}
}
// Populate the lower triangle from the symmetric entries of M:
for (a = 0; a < Mdim; a++) {
- for (b = a; b < Mdim; b++) { M[b * Mdim + a] = M[a * Mdim + b]; }
+ for (b = a; b < Mdim; b++) {
+ M[b * Mdim + a] = M[a * Mdim + b];
+ }
}
// Skip if undercoordinated
if (coordination[i] < zmin) continue;
- // Use gsl to get Minv, use Cholesky decomposition since the
+ // Use LAPACK to get Minv, use Cholesky decomposition since the
// polynomials are independent, M is symmetrix & positive-definite
- gM = gsl_matrix_view_array(M, Mdim, Mdim);
- gsl_error = gsl_linalg_cholesky_decomp1(&gM.matrix);
+ const char uplo = 'U';
+ dpotrf_(&uplo, &Mdim, M, &Mdim, &lapack_error);
- if (gsl_error) {
- //Revert to uncorrected SPH for this particle
- gsl_error_flag = 1;
- gsl_error_tags.insert(tag[i]);
+ if (lapack_error) {
+ // Revert to uncorrected SPH for this particle
+ lapack_error_flag = 1;
+ lapack_error_tags.insert(tag[i]);
- //check if not positive-definite
- if (gsl_error != GSL_EDOM)
- error->warning(FLERR, "Failed decomposition in rheo/kernel, gsl_error = {}", gsl_error);
+ // check if not positive-definite
+ if (lapack_error > 0)
+ error->warning(FLERR, "Failed DPOTRF2 decomposition in rheo/kernel, info = {}",
+ lapack_error);
continue;
}
- gsl_linalg_cholesky_invert(&gM.matrix); //M is now M^-1
+ // M is now M^-1
+ dpotri_(&uplo, &Mdim, M, &Mdim, &lapack_error);
+
+ // make result matrix symmetric
+ for (int i = 0; i < Mdim; ++i) {
+ for (int j = i+1; j < Mdim; ++j) {
+ M[i * Mdim + j] = M[j * Mdim + i];
+ }
+ }
// Correction coefficients are columns of M^-1 multiplied by an appropriate coefficient
// Solve the linear system several times to get coefficientns
diff --git a/src/RHEO/compute_rheo_kernel.h b/src/RHEO/compute_rheo_kernel.h
index 20516255be1..8b70509e6a9 100644
--- a/src/RHEO/compute_rheo_kernel.h
+++ b/src/RHEO/compute_rheo_kernel.h
@@ -53,8 +53,8 @@ class ComputeRHEOKernel : public Compute {
private:
int comm_stage, comm_forward_save;
int interface_flag;
- int gsl_error_flag;
- std::unordered_set gsl_error_tags;
+ int lapack_error_flag;
+ std::unordered_set lapack_error_tags;
int corrections_calculated;
int kernel_style, zmin, dim, Mdim, ncor;
diff --git a/src/YAFF/angle_mm3.cpp b/src/YAFF/angle_mm3.cpp
index 3ff7df1653b..920041f7e94 100644
--- a/src/YAFF/angle_mm3.cpp
+++ b/src/YAFF/angle_mm3.cpp
@@ -327,3 +327,15 @@ void AngleMM3::born_matrix(int type, int i1, int i2, int i3, double &du, double
du = -k2[type] * df / s;
du2 = k2[type] * (d2f - df * c / s) / (s * s) ;
}
+
+/* ----------------------------------------------------------------------
+ return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *AngleMM3::extract(const char *str, int &dim)
+{
+ dim = 1;
+ if (strcmp(str, "k2") == 0) return (void *) k2;
+ if (strcmp(str, "theta0") == 0) return (void *) theta0;
+ return nullptr;
+}
diff --git a/src/YAFF/angle_mm3.h b/src/YAFF/angle_mm3.h
index 22f5bd746cd..126b275e725 100644
--- a/src/YAFF/angle_mm3.h
+++ b/src/YAFF/angle_mm3.h
@@ -36,6 +36,7 @@ class AngleMM3 : public Angle {
void write_data(FILE *) override;
double single(int, int, int, int) override;
void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override;
+ void *extract(const char *, int &) override;
protected:
double *theta0, *k2;
diff --git a/src/YAFF/bond_mm3.cpp b/src/YAFF/bond_mm3.cpp
index 31ce2dad3e9..b3e69881e12 100644
--- a/src/YAFF/bond_mm3.cpp
+++ b/src/YAFF/bond_mm3.cpp
@@ -238,3 +238,15 @@ void BondMM3::born_matrix(int type, double rsq, int /*i*/, int /*j*/, double &du
du = 2.0 * k2[type] * dr + 3.0 * K3 * dr2 + 4.0 * K4 * dr3;
du2 = 2.0 * k2[type] + 6.0 * K3 * dr + 12.0 * K4 * dr2;
}
+
+/* ----------------------------------------------------------------------
+ return ptr to internal members upon request
+------------------------------------------------------------------------ */
+
+void *BondMM3::extract(const char *str, int &dim)
+{
+ dim = 1;
+ if (strcmp(str, "k2") == 0) return (void *) k2;
+ if (strcmp(str, "r0") == 0) return (void *) r0;
+ return nullptr;
+}
diff --git a/src/YAFF/bond_mm3.h b/src/YAFF/bond_mm3.h
index ea89ac826d0..b9ebf464bbe 100644
--- a/src/YAFF/bond_mm3.h
+++ b/src/YAFF/bond_mm3.h
@@ -36,6 +36,7 @@ class BondMM3 : public Bond {
void write_data(FILE *) override;
double single(int, double, int, int, double &) override;
void born_matrix(int, double, int, int, double &, double &) override;
+ void *extract(const char *, int &) override;
protected:
double *r0, *k2;
diff --git a/src/angle_hybrid.cpp b/src/angle_hybrid.cpp
index a015882a15a..1261a78176a 100644
--- a/src/angle_hybrid.cpp
+++ b/src/angle_hybrid.cpp
@@ -320,6 +320,14 @@ void AngleHybrid::init_style()
if (styles[m]) styles[m]->init_style();
}
+/* ---------------------------------------------------------------------- */
+
+int AngleHybrid::check_itype(int itype, char *substyle)
+{
+ if (strcmp(keywords[map[itype]], substyle) == 0) return 1;
+ return 0;
+}
+
/* ----------------------------------------------------------------------
return an equilbrium angle length
------------------------------------------------------------------------- */
diff --git a/src/angle_hybrid.h b/src/angle_hybrid.h
index a6da29245eb..a84096b2976 100644
--- a/src/angle_hybrid.h
+++ b/src/angle_hybrid.h
@@ -42,8 +42,10 @@ class AngleHybrid : public Angle {
double single(int, int, int, int) override;
double memory_usage() override;
+ int check_itype(int, char *);
+
protected:
- int *map; // which style each angle type points to
+ int *map; // which style each angle type points to
int *nanglelist; // # of angles in sub-style anglelists
int *maxangle; // max # of angles sub-style lists can store
int ***anglelist; // anglelist for each sub-style
diff --git a/src/bond_hybrid.cpp b/src/bond_hybrid.cpp
index 307cbd72fd2..bd5badb54cc 100644
--- a/src/bond_hybrid.cpp
+++ b/src/bond_hybrid.cpp
@@ -385,6 +385,14 @@ void BondHybrid::init_style()
else map[0] = -1;
}
+/* ---------------------------------------------------------------------- */
+
+int BondHybrid::check_itype(int itype, char *substyle)
+{
+ if (strcmp(keywords[map[itype]], substyle) == 0) return 1;
+ return 0;
+}
+
/* ----------------------------------------------------------------------
return an equilbrium bond length
------------------------------------------------------------------------- */
diff --git a/src/bond_hybrid.h b/src/bond_hybrid.h
index ba520b81b41..d93b5c75580 100644
--- a/src/bond_hybrid.h
+++ b/src/bond_hybrid.h
@@ -44,6 +44,8 @@ class BondHybrid : public Bond {
double single(int, double, int, int, double &) override;
double memory_usage() override;
+ int check_itype(int, char *);
+
protected:
int *map; // which style each bond type points to
int has_quartic; // which style, if any is a quartic bond style
diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp
index cad157f2be3..c725707b299 100644
--- a/src/fix_adapt.cpp
+++ b/src/fix_adapt.cpp
@@ -15,8 +15,10 @@
#include "fix_adapt.h"
#include "angle.h"
+#include "angle_hybrid.h"
#include "atom.h"
#include "bond.h"
+#include "bond_hybrid.h"
#include "domain.h"
#include "error.h"
#include "fix_store_atom.h"
@@ -386,11 +388,15 @@ void FixAdapt::init()
if (utils::strmatch(force->pair_style,"^hybrid")) {
auto pair = dynamic_cast(force->pair);
- for (i = ad->ilo; i <= ad->ihi; i++)
- for (j = MAX(ad->jlo,i); j <= ad->jhi; j++)
- if (!pair->check_ijtype(i,j,pstyle))
- error->all(FLERR,"Fix adapt type pair range is not valid "
- "for pair hybrid sub-style {}", pstyle);
+ if (pair) {
+ for (i = ad->ilo; i <= ad->ihi; i++) {
+ for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) {
+ if (!pair->check_ijtype(i,j,pstyle))
+ error->all(FLERR,"Fix adapt type pair range is not valid "
+ "for pair hybrid sub-style {}", pstyle);
+ }
+ }
+ }
}
delete[] pstyle;
@@ -416,8 +422,16 @@ void FixAdapt::init()
if (ad->bdim == 1) ad->vector = (double *) ptr;
- if (utils::strmatch(force->bond_style,"^hybrid"))
- error->all(FLERR,"Fix adapt does not support bond_style hybrid");
+ if (utils::strmatch(force->bond_style,"^hybrid")) {
+ auto bond = dynamic_cast(force->bond);
+ if (bond) {
+ for (i = ad->ilo; i <= ad->ihi; i++) {
+ if (!bond->check_itype(i,bstyle))
+ error->all(FLERR,"Fix adapt type bond range is not valid "
+ "for pair hybrid sub-style {}", bstyle);
+ }
+ }
+ }
delete[] bstyle;
@@ -442,8 +456,16 @@ void FixAdapt::init()
if (ad->adim == 1) ad->vector = (double *) ptr;
- if (utils::strmatch(force->angle_style,"^hybrid"))
- error->all(FLERR,"Fix adapt does not support angle_style hybrid");
+ if (utils::strmatch(force->angle_style,"^hybrid")) {
+ auto angle = dynamic_cast(force->angle);
+ if (angle) {
+ for (i = ad->ilo; i <= ad->ihi; i++) {
+ if (!angle->check_itype(i,astyle))
+ error->all(FLERR,"Fix adapt type angle range is not valid "
+ "for pair hybrid sub-style {}", astyle);
+ }
+ }
+ }
delete[] astyle;
diff --git a/src/fix_ave_chunk.cpp b/src/fix_ave_chunk.cpp
index d9723cec9fd..6a3c2e20326 100644
--- a/src/fix_ave_chunk.cpp
+++ b/src/fix_ave_chunk.cpp
@@ -153,7 +153,7 @@ FixAveChunk::FixAveChunk(LAMMPS *lmp, int narg, char **arg) :
while (iarg < nargnew) {
if (strcmp(arg[iarg],"norm") == 0) {
- if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk norm", error);
+ if (iarg+2 > nargnew) utils::missing_cmd_args(FLERR, "fix ave/chunk norm", error);
if (strcmp(arg[iarg+1],"all") == 0) {
normflag = ALL;
scaleflag = ATOM;
@@ -166,13 +166,13 @@ FixAveChunk::FixAveChunk(LAMMPS *lmp, int narg, char **arg) :
} else error->all(FLERR,"Unknown fix ave/chunk norm mode: {}", arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"ave") == 0) {
- if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk ave", error);
+ if (iarg+2 > nargnew) utils::missing_cmd_args(FLERR, "fix ave/chunk ave", error);
if (strcmp(arg[iarg+1],"one") == 0) ave = ONE;
else if (strcmp(arg[iarg+1],"running") == 0) ave = RUNNING;
else if (strcmp(arg[iarg+1],"window") == 0) ave = WINDOW;
else error->all(FLERR,"Unknown fix ave/chunk ave mode: {}", arg[iarg+1]);
if (ave == WINDOW) {
- if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk ave window", error);
+ if (iarg+3 > nargnew) utils::missing_cmd_args(FLERR, "fix ave/chunk ave window", error);
nwindow = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
if (nwindow <= 0) error->all(FLERR,"Illegal fix ave/chunk number of windows: {}", nwindow);
}
@@ -180,21 +180,21 @@ FixAveChunk::FixAveChunk(LAMMPS *lmp, int narg, char **arg) :
if (ave == WINDOW) iarg++;
} else if (strcmp(arg[iarg],"bias") == 0) {
- if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk bias", error);
+ if (iarg+2 > nargnew) utils::missing_cmd_args(FLERR, "fix ave/chunk bias", error);
biasflag = 1;
id_bias = utils::strdup(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"adof") == 0) {
- if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk adof", error);
+ if (iarg+2 > nargnew) utils::missing_cmd_args(FLERR, "fix ave/chunk adof", error);
adof = utils::numeric(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"cdof") == 0) {
- if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk cdof", error);
+ if (iarg+2 > nargnew) utils::missing_cmd_args(FLERR, "fix ave/chunk cdof", error);
cdof = utils::numeric(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if ((strcmp(arg[iarg],"file") == 0) || (strcmp(arg[iarg],"append") == 0)) {
- if (iarg+2 > narg)
+ if (iarg+2 > nargnew)
utils::missing_cmd_args(FLERR, std::string("fix ave/chunk ")+arg[iarg], error);
if (comm->me == 0) {
if (strcmp(arg[iarg],"file") == 0) fp = fopen(arg[iarg+1],"w");
@@ -208,23 +208,23 @@ FixAveChunk::FixAveChunk(LAMMPS *lmp, int narg, char **arg) :
overwrite = 1;
iarg += 1;
} else if (strcmp(arg[iarg],"format") == 0) {
- if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk format", error);
+ if (iarg+2 > nargnew) utils::missing_cmd_args(FLERR, "fix ave/chunk format", error);
delete[] format_user;
format_user = utils::strdup(arg[iarg+1]);
format = format_user;
iarg += 2;
} else if (strcmp(arg[iarg],"title1") == 0) {
- if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk title1", error);
+ if (iarg+2 > nargnew) utils::missing_cmd_args(FLERR, "fix ave/chunk title1", error);
delete[] title1;
title1 = utils::strdup(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"title2") == 0) {
- if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk title2", error);
+ if (iarg+2 > nargnew) utils::missing_cmd_args(FLERR, "fix ave/chunk title2", error);
delete[] title2;
title2 = utils::strdup(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"title3") == 0) {
- if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk title3", error);
+ if (iarg+2 > nargnew) utils::missing_cmd_args(FLERR, "fix ave/chunk title3", error);
delete[] title3;
title3 = utils::strdup(arg[iarg+1]);
iarg += 2;
diff --git a/src/fix_ave_grid.cpp b/src/fix_ave_grid.cpp
index ca89c918ba6..1b69c5644c8 100644
--- a/src/fix_ave_grid.cpp
+++ b/src/fix_ave_grid.cpp
@@ -199,14 +199,14 @@ FixAveGrid::FixAveGrid(LAMMPS *lmp, int narg, char **arg) :
while (iarg < nargnew) {
if (strcmp(arg[iarg],"discard") == 0) {
- if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/grid command");
+ if (iarg+2 > nargnew) error->all(FLERR,"Illegal fix ave/grid command");
if (strcmp(arg[iarg+1],"yes") == 0) discardflag = DISCARD;
else if (strcmp(arg[iarg+1],"no") == 0) discardflag = KEEP;
else error->all(FLERR,"Illegal fix ave/grid command");
iarg += 2;
} else if (strcmp(arg[iarg],"norm") == 0) {
- if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/grid command");
+ if (iarg+2 > nargnew) error->all(FLERR,"Illegal fix ave/grid command");
if (strcmp(arg[iarg+1],"all") == 0) normflag = ALL;
else if (strcmp(arg[iarg+1],"sample") == 0) normflag = SAMPLE;
else if (strcmp(arg[iarg+1],"none") == 0) normflag = NONORM;
@@ -214,13 +214,13 @@ FixAveGrid::FixAveGrid(LAMMPS *lmp, int narg, char **arg) :
iarg += 2;
} else if (strcmp(arg[iarg],"ave") == 0) {
- if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/grid command");
+ if (iarg+2 > nargnew) error->all(FLERR,"Illegal fix ave/grid command");
if (strcmp(arg[iarg+1],"one") == 0) aveflag = ONE;
else if (strcmp(arg[iarg+1],"running") == 0) aveflag = RUNNING;
else if (strcmp(arg[iarg+1],"window") == 0) aveflag = WINDOW;
else error->all(FLERR,"Illegal fix ave/grid command");
if (aveflag == WINDOW) {
- if (iarg+3 > narg) error->all(FLERR,"Illegal fix ave/grid command");
+ if (iarg+3 > nargnew) error->all(FLERR,"Illegal fix ave/grid command");
nwindow = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
if (nwindow <= 0) error->all(FLERR,"Illegal fix ave/grid command");
iarg++;
@@ -228,19 +228,19 @@ FixAveGrid::FixAveGrid(LAMMPS *lmp, int narg, char **arg) :
iarg += 2;
} else if (strcmp(arg[iarg],"bias") == 0) {
- if (iarg+2 > narg)
+ if (iarg+2 > nargnew)
error->all(FLERR,"Illegal fix ave/grid command");
biasflag = 1;
id_bias = utils::strdup(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"adof") == 0) {
- if (iarg+2 > narg)
+ if (iarg+2 > nargnew)
error->all(FLERR,"Illegal fix ave/grid command");
adof = utils::numeric(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"cdof") == 0) {
- if (iarg+2 > narg)
+ if (iarg+2 > nargnew)
error->all(FLERR,"Illegal fix ave/grid command");
cdof = utils::numeric(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
diff --git a/src/region.h b/src/region.h
index f273485dce1..19fdec31c71 100644
--- a/src/region.h
+++ b/src/region.h
@@ -20,6 +20,8 @@ namespace LAMMPS_NS {
class Region : protected Pointers {
public:
+ enum { CONSTANT, VARIABLE };
+
char *id, *style;
Region **reglist;
int interior; // 1 for interior, 0 for exterior
diff --git a/src/region_block.cpp b/src/region_block.cpp
index efa3d8ca6ab..93760168438 100644
--- a/src/region_block.cpp
+++ b/src/region_block.cpp
@@ -23,8 +23,6 @@
using namespace LAMMPS_NS;
-enum { CONSTANT, VARIABLE };
-
static constexpr double BIG = 1.0e20;
/* ---------------------------------------------------------------------- */
diff --git a/src/region_cone.cpp b/src/region_cone.cpp
index dc37eeefe3b..401ed537359 100644
--- a/src/region_cone.cpp
+++ b/src/region_cone.cpp
@@ -27,8 +27,6 @@
using namespace LAMMPS_NS;
-enum { CONSTANT, VARIABLE };
-
static constexpr double BIG = 1.0e20;
/* ---------------------------------------------------------------------- */
diff --git a/src/region_cylinder.cpp b/src/region_cylinder.cpp
index 11783dc1258..2ad0ba82f5b 100644
--- a/src/region_cylinder.cpp
+++ b/src/region_cylinder.cpp
@@ -26,8 +26,6 @@ using namespace LAMMPS_NS;
static constexpr double BIG = 1.0e20;
-enum { CONSTANT, VARIABLE };
-
/* ---------------------------------------------------------------------- */
RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
diff --git a/src/region_ellipsoid.cpp b/src/region_ellipsoid.cpp
index daabd621c89..a0b4b9e544c 100644
--- a/src/region_ellipsoid.cpp
+++ b/src/region_ellipsoid.cpp
@@ -23,8 +23,6 @@
using namespace LAMMPS_NS;
-enum { CONSTANT, VARIABLE };
-
static double GetRoot2D(double r0, double z0, double z1, double g);
static double GetRoot3D(double r0, double r1, double z0, double z1, double z2, double g);
diff --git a/src/region_plane.cpp b/src/region_plane.cpp
index 154b072633c..6dc162eeadb 100644
--- a/src/region_plane.cpp
+++ b/src/region_plane.cpp
@@ -14,6 +14,9 @@
#include "region_plane.h"
#include "error.h"
+#include "input.h"
+#include "update.h"
+#include "variable.h"
#include
@@ -21,13 +24,48 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
-RegPlane::RegPlane(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
+RegPlane::RegPlane(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg),
+ xstr(nullptr), ystr(nullptr), zstr(nullptr)
{
+ xvar = yvar = zvar = 0.0;
+
options(narg - 8, &arg[8]);
- xp = xscale * utils::numeric(FLERR, arg[2], false, lmp);
- yp = yscale * utils::numeric(FLERR, arg[3], false, lmp);
- zp = zscale * utils::numeric(FLERR, arg[4], false, lmp);
+ if (utils::strmatch(arg[2], "^v_")) {
+ xstr = utils::strdup(arg[2] + 2);
+ xp = 0.0;
+ xstyle = VARIABLE;
+ varshape = 1;
+ } else {
+ xp = xscale * utils::numeric(FLERR, arg[2], false, lmp);
+ xstyle = CONSTANT;
+ }
+
+ if (utils::strmatch(arg[3], "^v_")) {
+ ystr = utils::strdup(arg[3] + 2);
+ yp = 0.0;
+ ystyle = VARIABLE;
+ varshape = 1;
+ } else {
+ yp = yscale * utils::numeric(FLERR, arg[3], false, lmp);
+ ystyle = CONSTANT;
+ }
+
+ if (utils::strmatch(arg[4], "^v_")) {
+ zstr = utils::strdup(arg[4] + 2);
+ zp = 0.0;
+ zstyle = VARIABLE;
+ varshape = 1;
+ } else {
+ zp = zscale * utils::numeric(FLERR, arg[4], false, lmp);
+ zstyle = CONSTANT;
+ }
+
+ if (varshape) {
+ variable_check();
+ RegPlane::shape_update();
+ }
+
normal[0] = xscale * utils::numeric(FLERR, arg[5], false, lmp);
normal[1] = yscale * utils::numeric(FLERR, arg[6], false, lmp);
normal[2] = zscale * utils::numeric(FLERR, arg[7], false, lmp);
@@ -54,9 +92,20 @@ RegPlane::RegPlane(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
RegPlane::~RegPlane()
{
+ delete[] xstr;
+ delete[] ystr;
+ delete[] zstr;
delete[] contact;
}
+/* ---------------------------------------------------------------------- */
+
+void RegPlane::init()
+{
+ Region::init();
+ if (varshape) variable_check();
+}
+
/* ----------------------------------------------------------------------
inside = 1 if x,y,z is on normal side of plane or on plane
inside = 0 if x,y,z is on non-normal side of plane and not on plane
@@ -113,3 +162,45 @@ int RegPlane::surface_exterior(double *x, double cutoff)
}
return 0;
}
+
+/* ----------------------------------------------------------------------
+ change region shape via variable evaluation
+------------------------------------------------------------------------- */
+
+void RegPlane::shape_update()
+{
+ if (xstyle == VARIABLE) xp = xscale * input->variable->compute_equal(xvar);
+
+ if (ystyle == VARIABLE) yp = yscale * input->variable->compute_equal(yvar);
+
+ if (zstyle == VARIABLE) zp = zscale * input->variable->compute_equal(zvar);
+}
+
+/* ----------------------------------------------------------------------
+ error check on existence of variable
+------------------------------------------------------------------------- */
+
+void RegPlane::variable_check()
+{
+ if (xstyle == VARIABLE) {
+ xvar = input->variable->find(xstr);
+ if (xvar < 0) error->all(FLERR, "Variable {} for region plane does not exist", xstr);
+ if (!input->variable->equalstyle(xvar))
+ error->all(FLERR, "Variable {} for region plane is invalid style", xstr);
+ }
+
+ if (ystyle == VARIABLE) {
+ yvar = input->variable->find(ystr);
+ if (yvar < 0) error->all(FLERR, "Variable {} for region plane does not exist", ystr);
+ if (!input->variable->equalstyle(yvar))
+ error->all(FLERR, "Variable {} for region plane is invalid style", ystr);
+ }
+
+ if (zstyle == VARIABLE) {
+ zvar = input->variable->find(zstr);
+ if (zvar < 0) error->all(FLERR, "Variable {} for region plane does not exist", zstr);
+ if (!input->variable->equalstyle(zvar))
+ error->all(FLERR, "Variable {} for region plane is invalid style", zstr);
+ }
+}
+
diff --git a/src/region_plane.h b/src/region_plane.h
index 2025586a7c5..0e4ecda6d46 100644
--- a/src/region_plane.h
+++ b/src/region_plane.h
@@ -28,13 +28,23 @@ class RegPlane : public Region {
public:
RegPlane(class LAMMPS *, int, char **);
~RegPlane() override;
+ void init() override;
int inside(double, double, double) override;
int surface_interior(double *, double) override;
int surface_exterior(double *, double) override;
+ void shape_update() override;
private:
double xp, yp, zp;
double normal[3];
+
+ int xstyle, xvar;
+ int ystyle, yvar;
+ int zstyle, zvar;
+ char *xstr, *ystr, *zstr;
+
+ void variable_check();
+
};
} // namespace LAMMPS_NS
diff --git a/src/region_sphere.cpp b/src/region_sphere.cpp
index cd20a697d42..f4499789381 100644
--- a/src/region_sphere.cpp
+++ b/src/region_sphere.cpp
@@ -22,8 +22,6 @@
using namespace LAMMPS_NS;
-enum { CONSTANT, VARIABLE };
-
/* ---------------------------------------------------------------------- */
RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) :
diff --git a/src/variable.cpp b/src/variable.cpp
index b2f6c2882c2..279c14d9997 100644
--- a/src/variable.cpp
+++ b/src/variable.cpp
@@ -613,7 +613,7 @@ void Variable::set(int narg, char **arg)
// unrecognized variable style
- } else error->all(FLERR,"Unknown variable keyword: {}", arg[1]);
+ } else error->all(FLERR,"Unknown variable style: {}", arg[1]);
// set name of variable, if not replacing one flagged with replaceflag
// name must be all alphanumeric chars or underscores
diff --git a/unittest/commands/test_variables.cpp b/unittest/commands/test_variables.cpp
index 2390b1b675b..c7686cbf12d 100644
--- a/unittest/commands/test_variables.cpp
+++ b/unittest/commands/test_variables.cpp
@@ -206,7 +206,7 @@ TEST_F(VariableTest, CreateDelete)
TEST_FAILURE(".*ERROR: Invalid variable loop argument: -1.*",
command("variable dummy loop -1"););
TEST_FAILURE(".*ERROR: Illegal variable loop command.*", command("variable dummy loop 10 1"););
- TEST_FAILURE(".*ERROR: Unknown variable keyword: xxx.*", command("variable dummy xxxx"););
+ TEST_FAILURE(".*ERROR: Unknown variable style: xxx.*", command("variable dummy xxxx"););
TEST_FAILURE(".*ERROR: Cannot redefine variable as a different style.*",
command("variable two string xxx"););
TEST_FAILURE(".*ERROR: Cannot redefine variable as a different style.*",
diff --git a/unittest/force-styles/tests/angle-charmm.yaml b/unittest/force-styles/tests/angle-charmm.yaml
index 52e78abaf6a..8ff51183901 100644
--- a/unittest/force-styles/tests/angle-charmm.yaml
+++ b/unittest/force-styles/tests/angle-charmm.yaml
@@ -15,7 +15,9 @@ angle_coeff: ! |
3 40.0 120.0 35.0 2.410
4 33.0 108.5 30.0 2.163
equilibrium: 4 1.9216075064457567 1.9425514574696887 2.0943951023931953 1.8936822384138476
-extract: ! ""
+extract: ! |
+ k 1
+ theta0 1
natoms: 29
init_energy: 85.42486388459771
init_stress: ! |2-
diff --git a/unittest/force-styles/tests/angle-class2.yaml b/unittest/force-styles/tests/angle-class2.yaml
index ae2e3ff5ee7..8901157d171 100644
--- a/unittest/force-styles/tests/angle-class2.yaml
+++ b/unittest/force-styles/tests/angle-class2.yaml
@@ -23,7 +23,11 @@ angle_coeff: ! |
3 ba 10.0 10.0 1.5 1.5
4 ba 0.0 20.0 1.5 1.5
equilibrium: 4 1.9216075064457565 1.9373154697137058 2.0943951023931953 1.8936822384138474
-extract: ! ""
+extract: ! |
+ k2 1
+ k3 1
+ k4 1
+ theta0 1
natoms: 29
init_energy: 46.44089683774903
init_stress: ! |2-
diff --git a/unittest/force-styles/tests/angle-cosine_periodic.yaml b/unittest/force-styles/tests/angle-cosine_periodic.yaml
index 5c8227fcbdd..ee3e5c14690 100644
--- a/unittest/force-styles/tests/angle-cosine_periodic.yaml
+++ b/unittest/force-styles/tests/angle-cosine_periodic.yaml
@@ -15,7 +15,10 @@ angle_coeff: ! |
3 50.0 -1 3
4 100.0 -1 4
equilibrium: 4 3.141592653589793 3.141592653589793 2.0943951023931957 2.356194490192345
-extract: ! ""
+extract: ! |
+ k 1
+ b 1
+ multiplicity 1
natoms: 29
init_energy: 1178.5476942873006
init_stress: ! |2-
diff --git a/unittest/force-styles/tests/angle-cosine_squared_restricted.yaml b/unittest/force-styles/tests/angle-cosine_squared_restricted.yaml
index 400babb3c03..341ccb39194 100644
--- a/unittest/force-styles/tests/angle-cosine_squared_restricted.yaml
+++ b/unittest/force-styles/tests/angle-cosine_squared_restricted.yaml
@@ -17,7 +17,9 @@ angle_coeff: ! |
3 50.0 120.0
4 100.0 108.5
equilibrium: 4 1.9216075064457567 1.9373154697137058 2.0943951023931953 1.8936822384138476
-extract: ! ""
+extract: ! |
+ k 1
+ theta0 1
natoms: 29
init_energy: 43.16721849625078
init_stress: ! |2-
diff --git a/unittest/force-styles/tests/angle-dipole.yaml b/unittest/force-styles/tests/angle-dipole.yaml
index 877ffa19c7f..3914973f0dd 100644
--- a/unittest/force-styles/tests/angle-dipole.yaml
+++ b/unittest/force-styles/tests/angle-dipole.yaml
@@ -20,7 +20,9 @@ angle_coeff: ! |
3 50.0 120.0
4 100.0 108.5
equilibrium: 4 1.9216075064457565 1.9373154697137058 2.0943951023931953 1.8936822384138474
-extract: ! ""
+extract: ! |
+ k 1
+ gamma0 1
natoms: 29
init_energy: 1003.6681304854917
init_stress: ! |2-
diff --git a/unittest/force-styles/tests/angle-fourier.yaml b/unittest/force-styles/tests/angle-fourier.yaml
index 61165c5a92b..275d62beca9 100644
--- a/unittest/force-styles/tests/angle-fourier.yaml
+++ b/unittest/force-styles/tests/angle-fourier.yaml
@@ -16,7 +16,11 @@ angle_coeff: ! |
3 50.0 0.0 0.0 1.0
4 100.0 0.3 0.3 0.3
equilibrium: 4 3.141592653589793 1.5707963267948966 1.5707963267948966 1.8234765819369754
-extract: ! ""
+extract: ! |
+ k 1
+ C0 1
+ C1 1
+ C2 1
natoms: 29
init_energy: 400.84036632010225
init_stress: ! |-
diff --git a/unittest/force-styles/tests/angle-fourier_simple.yaml b/unittest/force-styles/tests/angle-fourier_simple.yaml
index e1a394ee3a9..bd72e67912a 100644
--- a/unittest/force-styles/tests/angle-fourier_simple.yaml
+++ b/unittest/force-styles/tests/angle-fourier_simple.yaml
@@ -16,7 +16,10 @@ angle_coeff: ! |
3 50.0 1.0 3.0
4 100.0 -0.5 1.5
equilibrium: 4 3.141592653589793 1.5707963267948966 1.0471975511965976 2.0943951023931953
-extract: ! ""
+extract: ! |
+ k 1
+ C 1
+ N 1
natoms: 29
init_energy: 2474.0748013590646
init_stress: ! |-
diff --git a/unittest/force-styles/tests/angle-mm3.yaml b/unittest/force-styles/tests/angle-mm3.yaml
index 9fb94601837..381f43187ec 100644
--- a/unittest/force-styles/tests/angle-mm3.yaml
+++ b/unittest/force-styles/tests/angle-mm3.yaml
@@ -16,7 +16,9 @@ angle_coeff: ! |
3 50.0 120.0
4 100.0 108.5
equilibrium: 4 1.9216075064457565 1.9373154697137058 2.0943951023931953 1.8936822384138474
-extract: ! ""
+extract: ! |
+ k2 1
+ theta0 1
natoms: 29
init_energy: 44.72461548562619
init_stress: ! |2-
diff --git a/unittest/force-styles/tests/angle-quartic.yaml b/unittest/force-styles/tests/angle-quartic.yaml
index 6ded709e849..15ec06d82a2 100644
--- a/unittest/force-styles/tests/angle-quartic.yaml
+++ b/unittest/force-styles/tests/angle-quartic.yaml
@@ -16,7 +16,11 @@ angle_coeff: ! |
3 120.0 50.0 -9.5 -1.5
4 108.5 100.0 5.0 -2.0
equilibrium: 4 1.9216075064457565 1.9373154697137058 2.0943951023931953 1.8936822384138474
-extract: ! ""
+extract: ! |
+ k2 1
+ k3 1
+ k4 1
+ theta0 1
natoms: 29
init_energy: 41.0458477552901
init_stress: ! |2-
diff --git a/unittest/force-styles/tests/angle-spica.yaml b/unittest/force-styles/tests/angle-spica.yaml
index 7f88553c70c..46e82383496 100644
--- a/unittest/force-styles/tests/angle-spica.yaml
+++ b/unittest/force-styles/tests/angle-spica.yaml
@@ -20,7 +20,9 @@ angle_coeff: ! |
3 40.0 120.0
4 33.0 108.5
equilibrium: 4 1.9216075064457565 1.9425514574696887 2.0943951023931953 1.8936822384138474
-extract: ! ""
+extract: ! |
+ k 1
+ theta0 1
natoms: 29
init_energy: 38.36438529349082
init_stress: ! |2-
diff --git a/unittest/force-styles/tests/atomic-pair-meam_ms.yaml b/unittest/force-styles/tests/atomic-pair-meam_ms.yaml
index fff938d9407..d8dcd7b1eb4 100644
--- a/unittest/force-styles/tests/atomic-pair-meam_ms.yaml
+++ b/unittest/force-styles/tests/atomic-pair-meam_ms.yaml
@@ -2,7 +2,7 @@
lammps_version: 7 Feb 2024
tags: slow
date_generated: Wed Feb 28 17:07:42 2024
-epsilon: 2.5e-12
+epsilon: 2.5e-11
skip_tests:
prerequisites: ! |
pair meam/ms
diff --git a/unittest/force-styles/tests/atomic-pair-pedone.yaml b/unittest/force-styles/tests/atomic-pair-pedone.yaml
index ea97d9ee8c2..82c6405e65a 100644
--- a/unittest/force-styles/tests/atomic-pair-pedone.yaml
+++ b/unittest/force-styles/tests/atomic-pair-pedone.yaml
@@ -1,6 +1,6 @@
---
lammps_version: 7 Feb 2024
-tags:
+tags: unstable
date_generated: Tue Apr 9 07:44:34 2024
epsilon: 7.5e-13
skip_tests:
diff --git a/unittest/force-styles/tests/bond-fene_expand.yaml b/unittest/force-styles/tests/bond-fene_expand.yaml
index fc859d477cb..250f89af159 100644
--- a/unittest/force-styles/tests/bond-fene_expand.yaml
+++ b/unittest/force-styles/tests/bond-fene_expand.yaml
@@ -17,7 +17,12 @@ bond_coeff: ! |
4 650 2.4 0.015 1.2 0.15
5 450 2 0.018 1 0.09
equilibrium: 5 1.5550000000000002 1.117 1.321 1.3139999999999998 1.06
-extract: ! ""
+extract: ! |
+ k 1
+ r0 1
+ epsilon 1
+ sigma 1
+ shift 1
natoms: 29
init_energy: 5926.020859124294
init_stress: ! |-
diff --git a/unittest/force-styles/tests/bond-harmonic_shift.yaml b/unittest/force-styles/tests/bond-harmonic_shift.yaml
index 7a41c2c3cda..61212a468bd 100644
--- a/unittest/force-styles/tests/bond-harmonic_shift.yaml
+++ b/unittest/force-styles/tests/bond-harmonic_shift.yaml
@@ -17,7 +17,10 @@ bond_coeff: ! |
4 650.0 1.2 0.2
5 450.0 1.0 0.0
equilibrium: 5 1.5 1.1 1.3 1.2 1
-extract: ! ""
+extract: ! |
+ k 1
+ r0 1
+ r1 1
natoms: 29
init_energy: -9395.519982389222
init_stress: ! |-
diff --git a/unittest/force-styles/tests/bond-mm3.yaml b/unittest/force-styles/tests/bond-mm3.yaml
index eb7443f1c27..f5ba5c237c8 100644
--- a/unittest/force-styles/tests/bond-mm3.yaml
+++ b/unittest/force-styles/tests/bond-mm3.yaml
@@ -17,7 +17,9 @@ bond_coeff: ! |
4 650.0 1.2
5 450.0 1.0
equilibrium: 5 1.5 1.1 1.3 1.2 1
-extract: ! ""
+extract: ! |
+ k2 1
+ r0 1
natoms: 29
init_energy: 4.247265008273143
init_stress: ! |-
diff --git a/unittest/force-styles/tests/dihedral-cosine_squared_restricted.yaml b/unittest/force-styles/tests/dihedral-cosine_squared_restricted.yaml
index f67a093017e..3f4d217b9a3 100644
--- a/unittest/force-styles/tests/dihedral-cosine_squared_restricted.yaml
+++ b/unittest/force-styles/tests/dihedral-cosine_squared_restricted.yaml
@@ -1,8 +1,8 @@
---
lammps_version: 7 Feb 2024
-tags:
+tags:
date_generated: Sat Apr 13 11:41:16 2024
-epsilon: 2.5e-13
+epsilon: 1.0e-11
skip_tests:
prerequisites: ! |
atom full