ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
Stan Moore
2024-09-23 16:14:12 -06:00
parent 7e31a4f482 b95f7c2ab2
commit 80832ce763
93 changed files with 921 additions and 194 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
.github/workflows/check-vla.yml vendored
View File

@ -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 \

4
.github/workflows/compile-msvc.yml vendored
View File

@ -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

3
.github/workflows/full-regression.yml vendored
View File

@ -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

7
.github/workflows/quick-regression.yml vendored
View File

@ -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

37
.github/workflows/style-check.yml vendored Normal file
View File

@ -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

6
.github/workflows/unittest-linux.yml vendored
View File

@ -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 \

4
.github/workflows/unittest-macos.yml vendored
View File

@ -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

4
cmake/CMakeLists.txt
View File

@ -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()

48
cmake/Modules/Documentation.cmake
View File

@ -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

29
cmake/Modules/FindSphinx.cmake Normal file
View File

@ -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)

2
cmake/Modules/Packages/RHEO.cmake
View File

@ -1,2 +0,0 @@
find_package(GSL 2.6 REQUIRED)
target_link_libraries(lammps PRIVATE GSL::gsl)

1
cmake/presets/mingw-cross.cmake
View File

@ -67,6 +67,7 @@ set(WIN_PACKAGES
REACTION
REAXFF
REPLICA
RHEO
RIGID
SHOCK
SMTBQ

1
cmake/presets/most.cmake
View File

@ -60,6 +60,7 @@ set(ALL_PACKAGES
REACTION
REAXFF
REPLICA
RHEO
RIGID
SHOCK
SPH

1
cmake/presets/windows.cmake
View File

@ -60,6 +60,7 @@ set(WIN_PACKAGES
REACTION
REAXFF
REPLICA
RHEO
RIGID
SHOCK
SMTBQ

34
doc/src/Build_extras.rst
View File

@ -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) <https://www.gnu.org/software/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.<config>`` 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.
----------

3
doc/src/Build_settings.rst
View File

@ -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

4
doc/src/Intro_authors.rst
View File

@ -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 <tn_>`_
- 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/

78
doc/src/fix_adapt.rst
View File

@ -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
If :doc:`bond_style hybrid <bond_hybrid>` is used, *bstyle* should be a
sub-style name. The bond styles that currently work with fix adapt are:
+------------------------------------+------------+------------+
| :doc:`class2 <bond_class2>` | r0 | type bonds |
+------------------------------------+------------+------------+
| :doc:`fene <bond_fene>` | k,r0 | type bonds |
+------------------------------------+------------+------------+
| :doc:`fene/nm <bond_fene>` | k,r0 | type bonds |
+------------------------------------+------------+------------+
| :doc:`gromos <bond_gromos>` | k,r0 | type bonds |
+------------------------------------+------------+------------+
| :doc:`harmonic <bond_harmonic>` | k,r0 | type bonds |
+------------------------------------+------------+------------+
| :doc:`morse <bond_morse>` | r0 | type bonds |
+------------------------------------+------------+------------+
| :doc:`nonlinear <bond_nonlinear>` | epsilon,r0 | type bonds |
+------------------------------------+------------+------------+
+---------------------------------------------------+---------------------------+------------+
| :doc:`class2 <bond_class2>` | r0 | type bonds |
+---------------------------------------------------+---------------------------+------------+
| :doc:`fene <bond_fene>` | k,r0 | type bonds |
+---------------------------------------------------+---------------------------+------------+
| :doc:`fene/expand <bond_fene_expand>` | k,r0,epsilon,sigma,shift | type bonds |
+---------------------------------------------------+---------------------------+------------+
| :doc:`fene/nm <bond_fene>` | k,r0 | type bonds |
+---------------------------------------------------+---------------------------+------------+
| :doc:`gromos <bond_gromos>` | k,r0 | type bonds |
+---------------------------------------------------+---------------------------+------------+
| :doc:`harmonic <bond_harmonic>` | k,r0 | type bonds |
+---------------------------------------------------+---------------------------+------------+
| :doc:`harmonic/shift <bond_harmonic_shift>` | k,r0,r1 | type bonds |
+---------------------------------------------------+---------------------------+------------+
| :doc:`harmonic/restrain <bond_harmonic_restrain>` | k | type bonds |
+---------------------------------------------------+---------------------------+------------+
| :doc:`mm3 <bond_mm3>` | k,r0 | type bonds |
+---------------------------------------------------+---------------------------+------------+
| :doc:`morse <bond_morse>` | r0 | type bonds |
+---------------------------------------------------+---------------------------+------------+
| :doc:`nonlinear <bond_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
If :doc:`angle_style hybrid <angle_hybrid>` is used, *astyle* should be a
sub-style name. The angle styles that currently work with fix adapt are:
+------------------------------------+----------+-------------+
| :doc:`harmonic <angle_harmonic>` | k,theta0 | type angles |
+------------------------------------+----------+-------------+
| :doc:`cosine <angle_cosine>` | k | type angles |
+------------------------------------+----------+-------------+
+--------------------------------------------------------------------+-----------------+-------------+
| :doc:`harmonic <angle_harmonic>` | k,theta0 | type angles |
+--------------------------------------------------------------------+-----------------+-------------+
| :doc:`charmm <angle_charmm>` | k,theta0 | type angles |
+--------------------------------------------------------------------+-----------------+-------------+
| :doc:`class2 <angle_class2>` | k2,k3,k4,theta0 | type angles |
+--------------------------------------------------------------------+-----------------+-------------+
| :doc:`cosine <angle_cosine>` | k | type angles |
+--------------------------------------------------------------------+-----------------+-------------+
| :doc:`cosine/periodic <angle_cosine_periodic>` | k,b,n | type angles |
+--------------------------------------------------------------------+-----------------+-------------+
| :doc:`cosine/squared/restricted <angle_cosine_squared_restricted>` | k,theta0 | type angles |
+--------------------------------------------------------------------+-----------------+-------------+
| :doc:`dipole <angle_dipole>` | k,gamma0 | type angles |
+--------------------------------------------------------------------+-----------------+-------------+
| :doc:`fourier <angle_fourier>` | k,c0,c1,c2 | type angles |
+--------------------------------------------------------------------+-----------------+-------------+
| :doc:`fourier/simple <angle_fourier_simple>` | k,c,n | type angles |
+--------------------------------------------------------------------+-----------------+-------------+
| :doc:`mm3 <angle_mm3>` | k,theta0 | type angles |
+--------------------------------------------------------------------+-----------------+-------------+
| :doc:`quartic <angle_quartic>` | k2,k3,k4,theta0 | type angles |
+--------------------------------------------------------------------+-----------------+-------------+
| :doc:`spica <angle_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.

14
doc/src/region.rst
View File

@ -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 <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

1
doc/utils/sphinx-config/false_positives.txt
View File

@ -141,6 +141,7 @@ arg
arge
args
argv
Armv
arrhenius
Arun
arXiv

77
lib/linalg/dlauu2.cpp Normal file
View File

@ -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

101
lib/linalg/dlauum.cpp Normal file
View File

@ -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

40
lib/linalg/dpotri.cpp Normal file
View File

@ -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

17
lib/rheo/Makefile.lammps
View File

@ -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)

5
lib/rheo/Makefile.lammps.empty Normal file
View File

@ -0,0 +1,5 @@
# Settings that the LAMMPS build will import when this package library is used
rheo_SYSINC =
rheo_SYSLIB =
rheo_SYSPATH =

5
lib/rheo/Makefile.lammps.installed Normal file
View File

@ -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 =

5
lib/rheo/Makefile.lammps.linalg Normal file
View File

@ -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)

12
lib/rheo/README
View File

@ -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.

12
src/CG-SPICA/angle_spica.cpp
View File

@ -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;
}

1
src/CG-SPICA/angle_spica.h
View File

@ -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;

14
src/CLASS2/angle_class2.cpp
View File

@ -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;
}

1
src/CLASS2/angle_class2.h
View File

@ -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;

12
src/DIPOLE/angle_dipole.cpp
View File

@ -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;
}

1
src/DIPOLE/angle_dipole.h
View File

@ -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;

4
src/EXTRA-FIX/fix_ave_correlate_long.cpp
View File

@ -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;

13
src/EXTRA-MOLECULE/angle_cosine_periodic.cpp
View File

@ -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;
}

1
src/EXTRA-MOLECULE/angle_cosine_periodic.h
View File

@ -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;

12
src/EXTRA-MOLECULE/angle_cosine_squared_restricted.cpp
View File

@ -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;
}

1
src/EXTRA-MOLECULE/angle_cosine_squared_restricted.h
View File

@ -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;

13
src/EXTRA-MOLECULE/angle_fourier.cpp
View File

@ -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;
}

1
src/EXTRA-MOLECULE/angle_fourier.h
View File

@ -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;

13
src/EXTRA-MOLECULE/angle_fourier_simple.cpp
View File

@ -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;
}

1
src/EXTRA-MOLECULE/angle_fourier_simple.h
View File

@ -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;

14
src/EXTRA-MOLECULE/angle_quartic.cpp
View File

@ -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;
}

1
src/EXTRA-MOLECULE/angle_quartic.h
View File

@ -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;

13
src/EXTRA-MOLECULE/bond_harmonic_shift.cpp
View File

@ -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;
}

1
src/EXTRA-MOLECULE/bond_harmonic_shift.h
View File

@ -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;

3
src/ML-IAP/mliap_data.cpp
View File

@ -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;
}

12
src/MOLECULE/angle_charmm.cpp
View File

@ -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;
}

1
src/MOLECULE/angle_charmm.h
View File

@ -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;

15
src/MOLECULE/bond_fene_expand.cpp
View File

@ -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;
}

1
src/MOLECULE/bond_fene_expand.h
View File

@ -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;

1
src/RHEO/Install.sh
View File

@ -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

9
src/RHEO/README
View File

@ -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.

74
src/RHEO/compute_rheo_kernel.cpp
View File

@ -37,10 +37,6 @@
#include "utils.h"
#include <cmath>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
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

4
src/RHEO/compute_rheo_kernel.h
View File

@ -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<tagint> gsl_error_tags;
int lapack_error_flag;
std::unordered_set<tagint> lapack_error_tags;
int corrections_calculated;
int kernel_style, zmin, dim, Mdim, ncor;

12
src/YAFF/angle_mm3.cpp
View File

@ -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;
}

1
src/YAFF/angle_mm3.h
View File

@ -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;

12
src/YAFF/bond_mm3.cpp
View File

@ -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;
}

1
src/YAFF/bond_mm3.h
View File

@ -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;

8
src/angle_hybrid.cpp
View File

@ -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
------------------------------------------------------------------------- */

4
src/angle_hybrid.h
View File

@ -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

8
src/bond_hybrid.cpp
View File

@ -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
------------------------------------------------------------------------- */

2
src/bond_hybrid.h
View File

@ -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

40
src/fix_adapt.cpp
View File

@ -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<PairHybrid *>(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<BondHybrid *>(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<AngleHybrid *>(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;

22
src/fix_ave_chunk.cpp
View File

@ -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;

14
src/fix_ave_grid.cpp
View File

@ -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;

2
src/region.h
View File

@ -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

2
src/region_block.cpp
View File

@ -23,8 +23,6 @@
using namespace LAMMPS_NS;
enum { CONSTANT, VARIABLE };
static constexpr double BIG = 1.0e20;
/* ---------------------------------------------------------------------- */

2
src/region_cone.cpp
View File

@ -27,8 +27,6 @@
using namespace LAMMPS_NS;
enum { CONSTANT, VARIABLE };
static constexpr double BIG = 1.0e20;
/* ---------------------------------------------------------------------- */

2
src/region_cylinder.cpp
View File

@ -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) :

2
src/region_ellipsoid.cpp
View File

@ -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);

99
src/region_plane.cpp
View File

@ -14,6 +14,9 @@
#include "region_plane.h"
#include "error.h"
#include "input.h"
#include "update.h"
#include "variable.h"
#include <cmath>
@ -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);
}
}

10
src/region_plane.h
View File

@ -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

2
src/region_sphere.cpp
View File

@ -22,8 +22,6 @@
using namespace LAMMPS_NS;
enum { CONSTANT, VARIABLE };
/* ---------------------------------------------------------------------- */
RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) :

2
src/variable.cpp
View File

@ -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

2
unittest/commands/test_variables.cpp
View File

@ -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.*",

4
unittest/force-styles/tests/angle-charmm.yaml
View File

@ -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-

6
unittest/force-styles/tests/angle-class2.yaml
View File

@ -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-

5
unittest/force-styles/tests/angle-cosine_periodic.yaml
View File

@ -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-

4
unittest/force-styles/tests/angle-cosine_squared_restricted.yaml
View File

@ -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-

4
unittest/force-styles/tests/angle-dipole.yaml
View File

@ -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-

6
unittest/force-styles/tests/angle-fourier.yaml
View File

@ -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: ! |-

5
unittest/force-styles/tests/angle-fourier_simple.yaml
View File

@ -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: ! |-

4
unittest/force-styles/tests/angle-mm3.yaml
View File

@ -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-

6
unittest/force-styles/tests/angle-quartic.yaml
View File

@ -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-

4
unittest/force-styles/tests/angle-spica.yaml
View File

@ -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-

2
unittest/force-styles/tests/atomic-pair-meam_ms.yaml
View File

@ -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

2
unittest/force-styles/tests/atomic-pair-pedone.yaml
View File

@ -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:

7
unittest/force-styles/tests/bond-fene_expand.yaml
View File

@ -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: ! |-

5
unittest/force-styles/tests/bond-harmonic_shift.yaml
View File

@ -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: ! |-

4
unittest/force-styles/tests/bond-mm3.yaml
View File

@ -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: ! |-

4
unittest/force-styles/tests/dihedral-cosine_squared_restricted.yaml
View File

@ -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